5
$\begingroup$

The thermal conduction problem, described in polarcoordinates can be solved

t1 = 2; \[CapitalDelta]\[Phi] = 40 Degree; r1 = 0.161; r2 = 0.201;  
U = NDSolveValue[
{   Derivative[1, 0 , 0][u][t, r, \[CurlyPhi] ] +Derivative[0, 0 , 1][u][t,r, \[CurlyPhi] ] == \[Lambda] Laplacian[u[t, r, \[CurlyPhi]], {r, \[CurlyPhi]},"Polar"]+ 50 Boole[(0 < \[CurlyPhi] < \[CapitalDelta]\[Phi])]- 0.025 (u[t, r, \[CurlyPhi]] - 20)
, u[0, r, \[CurlyPhi]] == 100, u[t, r, 0] == u[t, r, 2 Pi]},
u,{r, r1, r2}, {\[CurlyPhi], 0, 2 Pi}, {t, 0, t1}
, Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid"}}] 

using Methodof Lines in an TensorProductGrid.

The solution

Quiet@Show[{Graphics[{Gray,Annulus[{0, 0}, r2 {1.05, 1.1}, {0, \[CapitalDelta]\[Phi]}], Line[{{0, 0}, {r2, 0}}],Line[{{0, 0}, r2 {Cos[\[CapitalDelta]\[Phi]], Sin[\[CapitalDelta]\[Phi]]}}]}],DensityPlot[U[#, Sqrt[x^2 + y^2], ArcTan[x, y]] , {x, -r2, r2}, {y, -r2,r2}, ColorFunction -> ColorData["ThermometerColors"] ,ColorFunctionScaling -> True , PlotLegends -> Automatic, RegionFunction -> Function[{x, y, z}, r1^2 < x^2 + y^2 < r2^2], PlotRange -> All, Exclusions -> None, Frame -> False]}] &[2]

enter image description here

looks fine

but NDSolve gives the message NDSolveValue::bcart: Warning: an insufficient number of boundary conditions have been specified for the direction of independent variable r. Artificial boundary effects may be present in the solution.

I tried to add two flux boundaries NeumannValue[0,r==r1] and NeumannValue[0,r==r2] to the pde but it doesn't work. Unfortunately Method->"FiniteElement" doesn't work.

My question: How can I avoid the error message and fullfill the two flux boundary conditions?

Thanks!

$\endgroup$
  • 2
    $\begingroup$ You did not define a $\lambda $. In the case of FEM, zero boundary conditions of Neumann type are automatically applied at boundaries $r=r_1,r_2$. So use Derivative[0, 1, 0][u][t, r1, \[CurlyPhi]] == 0, Derivative[0, 1, 0][u][t, r2, \[CurlyPhi]] == 0 $\endgroup$ – Alex Trounev Sep 6 '19 at 7:18
  • $\begingroup$ @AlexTrounev Thank you very much, it works! I only tried the NeumannValue, but couldn't use FEM because of the periodicty condition. $\endgroup$ – Ulrich Neumann Sep 6 '19 at 9:11
  • $\begingroup$ FEM also works well in this case. $\endgroup$ – Alex Trounev Sep 6 '19 at 9:32
  • $\begingroup$ @ AlexTrounev Unfortunately not. If I substitute TensorPorductGrid by FiniteElement NDSolve doesn't evaluate. Or would you suggest another way? $\endgroup$ – Ulrich Neumann Sep 6 '19 at 10:17
  • $\begingroup$ What is lambda equal to? $\endgroup$ – Alex Trounev Sep 6 '19 at 10:21
6
$\begingroup$

To solve the problem, we can apply FEM using explicit Euler in time. So that there were no messages, we organized a weak leak across the border.

Needs["NDSolve`FEM`"]; t1 = 2; \[CapitalDelta]\[Phi] = 
 40/180 Pi ; r1 = 0.161; r2 = 0.201; \[Lambda] = 1/2; 
f[x_] := Piecewise[{{50 , 0 <= x <= \[CapitalDelta]\[Phi]}, {0, 
    True}}]; mesh = ImplicitRegion[r1^2 <= x^2 + y^2 <= r2^2, {x, y}];
t0 = 1/20; n = 40;
U[0][x_, y_] := 100;
Do[U[t] = 
   NDSolveValue[(u[x, y] - U[t - t0][x, y])/
       t0 - \[Lambda] Laplacian[u[x, y], {x, y}] + (x D[u[x, y], y] - 
        y D[u[x, y], x]) - f[ArcTan[x, y]] + 0.025 (u[x, y] - 20) == 
     NeumannValue[10^-9 u[x, y], True], u, {x, y} \[Element] mesh, 
    Method -> {"FiniteElement", InterpolationOrder -> {u -> 2}, 
      "MeshOptions" -> {"MaxCellMeasure" -> 0.001, 
        "MeshOrder" -> 2}}];, {t, t0, n t0, t0}]
Table[DensityPlot[Evaluate[U[t][x, y]], {x, y} \[Element] mesh, 
  PlotLegends -> Automatic, PlotLabel -> Row[{"t = ", t*1.}], 
  ColorFunction -> ColorData["ThermometerColors"], PlotPoints -> 100, 
  Frame -> False, ColorFunctionScaling -> True], {t, 10 t0, n*t0,10 t0}]

Figure 1 For comparison, we used the 1D model, given the absence of gradients along r. We see that the data is consistent with FEM.

With[{r = (r1 + r2)/2}, 
 U1 = NDSolveValue[{Derivative[1, 0][u][t, \[CurlyPhi]] + 
      Derivative[0, 1][u][
       t, \[CurlyPhi]] == \[Lambda] D[
         u[t, \[CurlyPhi]], {\[CurlyPhi], 2}]/r^2 + f[\[CurlyPhi]] - 
      0.025 (u[t, \[CurlyPhi]] - 20), u[0, \[CurlyPhi]] == 100, 
    u[t, 0] == u[t, 2 Pi]}, u, {\[CurlyPhi], 0, 2 Pi}, {t, 0, t1}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", 
       "MinPoints" -> 200, "MaxPoints" -> 200, 
       "DifferenceOrder" -> 4}}]]

Table[DensityPlot[
  Evaluate[U1[t, ArcTan[x, y]]], {x, y} \[Element] mesh, 
  PlotLegends -> Automatic, PlotLabel -> Row[{"t = ", t*1.}], 
  ColorFunction -> ColorData["ThermometerColors"], PlotPoints -> 100, 
  Frame -> False, ColorFunctionScaling -> True], {t, 10 t0, n*t0, 
  10 t0}]

Figure 2 Also, in version 12, we can use FEM and time dependency. The result completely coincides with that obtained using explicit Euler:

U12 = NDSolveValue[{D[u[t, x, y], 
       t] - \[Lambda] Laplacian[
        u[t, x, y], {x, y}] + (x D[u[t, x, y], y] - 
        y D[u[t, x, y], x]) - f[ArcTan[x, y]] + 
      0.025 (u[t, x, y] - 20) == NeumannValue[10^-9 u[t, x, y], True],
     u[0, x, y] == 100}, u, {t, 0, t1}, {x, y} \[Element] mesh, 
   Method -> {"FiniteElement", 
     "MeshOptions" -> {"MaxCellMeasure" -> 0.001}}];
Table[DensityPlot[Evaluate[U12[t, x, y]], {x, y} \[Element] mesh, 
  PlotLegends -> Automatic, PlotLabel -> Row[{"t = ", t*1.}], 
  ColorFunction -> ColorData["ThermometerColors"], PlotPoints -> 100, 
  Frame -> False, ColorFunctionScaling -> True, PlotRange -> All], {t,
   10 t0, n*t0, 10 t0}]

Figure 3

The 2D "MethodOfLines" with periodic boundary conditions in polar coordinates gives a solution that does not have physical meaning. For example, the code below has no messages (unlike the one suggested by Ulrich Neumann)

U2 = NDSolveValue[{Derivative[1, 0, 0][u][t, r, \[CurlyPhi]] + 
     Derivative[0, 0, 1][u][t, 
      r, \[CurlyPhi]] == \[Lambda] Laplacian[
       u[t, r, \[CurlyPhi]], {r, \[CurlyPhi]}, "Polar"] + 
     f[\[CurlyPhi]] - 0.025 (u[t, r, \[CurlyPhi]] - 20), 
   u[0, r, \[CurlyPhi]] == 100, u[t, r, 0] == u[t, r, 2 Pi], 
   Derivative[0, 1, 0][u][t, r1, \[CurlyPhi]] == 0, 
   Derivative[0, 1, 0][u][t, r2, \[CurlyPhi]] == 0}, 
  u, {r, r1, r2}, {\[CurlyPhi], 0, 2 Pi}, {t, 0, t1}, 
  Method -> {"MethodOfLines", 
    "SpatialDiscretization" -> {"TensorProductGrid", 
      "MinPoints" -> 40, "MaxPoints" -> 100, 
      "DifferenceOrder" -> "Pseudospectral"}}]

Compare solution U2 and solution U1 (1D "MethodOfLines"). We see that in 2D the heat flux at the boundary $\phi =0=2\pi$ has a gap, which is devoid of physical meaning.

Plot[{U2[2, (r1 + r2)/2, phi], U1[2, phi]}, {phi, 0, 2 Pi}]

Figure 4

We add the condition of equality of heat fluxes at the boundary $2\pi=\phi =0$. In this case, the solution has the form

U2 = NDSolveValue[{Derivative[1, 0, 0][u][t, r, \[CurlyPhi]] + 
     Derivative[0, 0, 1][u][t, 
      r, \[CurlyPhi]] == \[Lambda] Laplacian[
       u[t, r, \[CurlyPhi]], {r, \[CurlyPhi]}, "Polar"] + 
     f[\[CurlyPhi]] - 0.025 (u[t, r, \[CurlyPhi]] - 20), 
   u[0, r, \[CurlyPhi]] == 100, u[t, r, 0] == u[t, r, 2 Pi], 
   Derivative[0, 0, 1][u][t, r, 0] == 
    Derivative[0, 0, 1][u][t, r, 2 Pi], 
   Derivative[0, 1, 0][u][t, r1, \[CurlyPhi]] == 0, 
   Derivative[0, 1, 0][u][t, r2, \[CurlyPhi]] == 0}, 
  u, {r, r1, r2}, {\[CurlyPhi], 0, 2 Pi}, {t, 0, t1}, 
  Method -> {"MethodOfLines", 
    "SpatialDiscretization" -> {"TensorProductGrid", 
      "MinPoints" -> 40, "MaxPoints" -> 100, 
      "DifferenceOrder" -> "Pseudospectral"}}]

Now we see an even greater difference between solutions of the equation obtained using 2D "MethodOfLines" and FEM (1D and FEM coincide on this scale).

Plot[{U2[2, (r1 + r2)/2, phi], U1[2, phi], 
  U12[2, Cos[phi] (r1 + r2)/2, Sin[phi] (r1 + r2)/2]}, {phi, 0, 2 Pi}, 
 PlotRange -> All, PlotLegends -> {"2D", "1D", "2D FEM"}, 
 AxesLabel -> {"\[CurlyPhi]", "u"}]

Figure 5

$\endgroup$
  • $\begingroup$ Thank you for your Help, I have to elaborate.... $\endgroup$ – Ulrich Neumann Sep 7 '19 at 15:36
  • $\begingroup$ What do you mean with: "Now we see an even greater difference with 1D and FEM (1D and FEM coincide on this scale)." Do you mean ....with 2D (TPG) and FEM... ? $\endgroup$ – user21 Sep 9 '19 at 7:02
  • 1
    $\begingroup$ @user21 Thank you, corrected. I do not have version 11.3, but in version 12 this code works. $\endgroup$ – Alex Trounev Sep 9 '19 at 10:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.