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!

Improve this question
$\endgroup$
7
  • 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
    Commented Sep 6, 2019 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
    Commented Sep 6, 2019 at 9:11
  • $\begingroup$ FEM also works well in this case. $\endgroup$
    – Alex Trounev
    Commented Sep 6, 2019 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
    Commented Sep 6, 2019 at 10:17
  • $\begingroup$ What is lambda equal to? $\endgroup$
    – Alex Trounev
    Commented Sep 6, 2019 at 10:21

1 Answer 1

Reset to default
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

Improve this answer
$\endgroup$
3
  • $\begingroup$ Thank you for your Help, I have to elaborate.... $\endgroup$
    – Ulrich Neumann
    Commented Sep 7, 2019 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
    Commented Sep 9, 2019 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
    Commented Sep 9, 2019 at 10:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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