2
$\begingroup$

I am trying to solve heat diffusion equation on hollow cylinder with constant DirichletCondition on inner radius and zero NeumannCondition (no conduction) on outer radius. Initial condition is linearly with r decreasing temperature. My code:

h = 10;
cyl1 = Cylinder[{{0, 0, 0}, {0, 0, h}}, 2];
cyl2 = Cylinder[{{0, 0, 0}, {0, 0, h}}, 4];
cyl = DiscretizeRegion[RegionDifference[cyl2, cyl1]];
Laplacian[u[t, r, θ, z], {r, θ, z}, "Cylindrical"]
ClearAll[r]
ifun = NDSolveValue[{D[u[t, r, θ, z], t] == 
     Laplacian[u[t, r, θ, z], {r, θ, z}, 
       "Cylindrical"] + NeumannValue[0.1, (0 < z < h) && (r > 3.9)], 
    DirichletCondition[u[t, r, θ, z] == 50, r < 2.1], 
    u[0.001, r, θ, z] == 50 - ((50 - 10)/(4 - 2))*r}, 
   u, {t, 1, 5}, {r, θ, z} ∈ cyl, Method -> Automatic];
SliceContourPlot3D[
 ifun[4, r, θ, z], {r, θ, z} ∈ cyl, 
 ColorFunction -> "TemperatureMap", Boxed -> False, Axes -> None]

The question is why is NDSolveValue taking so long.

EDIT: How to convince Mathematica to interpret (r,theta,z) as cylindrical (not cartesian) coordinates of a region cyl?

$\endgroup$
3
  • 1
    $\begingroup$ Hi kular. On your machine, how long does it take to run? $\endgroup$ – QuantumDot Apr 24 '20 at 20:02
  • $\begingroup$ I'm running it for 2 hours, still no result. $\endgroup$ – fasdgr Apr 24 '20 at 20:07
  • 1
    $\begingroup$ @kular You made a typo with Cylinder[]. In the cylindrical coordinate it is not cylinder, but region $2\le r\le 2, 0\le \theta \le 2\pi, 0\le z\le h$ $\endgroup$ – Alex Trounev Apr 25 '20 at 19:21
6
$\begingroup$

Try

h = 10;
innerR = 2;
outerR = 4;
cyl1 = Cylinder[{{0, 0, 0}, {0, 0, h}}, innerR];
cyl2 = Cylinder[{{0, 0, 0}, {0, 0, h}}, outerR];
cyl = ToElementMesh[RegionDifference[cyl2, cyl1], MaxCellMeasure -> .5];
cyl["Wireframe"]

Mathematica graphics

pde = D[u[t, r, θ, z], t] == Laplacian[u[t, r, θ, z], {r, θ, z}, "Cylindrical"] + 
    NeumannValue[1/10, (0 < z < h) && (r > 39/10)];

ic = u[1/1000, r, θ, z] == 50 - ((50 - 10)/(4 - 2))*r;

ifun = NDSolveValue[{pde, DirichletCondition[u[t, r, θ, z] == 50, 
    r < 21/10], ic}, u, {t, 1/1000, 5}, {r, θ, z} ∈ cyl]

Mathematica graphics

Finishes right away.

$\endgroup$
1
  • $\begingroup$ There is a typo with Cylinder[]. In the cylindrical coordinates it is not cylinder, but Coboud[] region - see Tim Laska answer and my also. $\endgroup$ – Alex Trounev Apr 25 '20 at 23:02
4
$\begingroup$

In cylindrical space, the annular region would look like a cuboid. Since there is no $\theta$ dependence in your model, you may as well do an axisymmetric model. Fortunately, Mathematica has a heat transfer verification manual to give a guide on how to setup a verified case, from which, I will shamelessly copy and paste.

ClearAll[HeatTransferModel]
HeatTransferModel[T_, X_List, k_, \[Rho]_, Cp_, Velocity_, Source_] :=
  Module[{V, Q, a = k},
  V = If[Velocity === "NoFlow", 
    0, \[Rho]*Cp*Velocity.Inactive[Grad][T, X]];
  Q = If[Source === "NoSource", 0, Source];
  If[ FreeQ[a, _?VectorQ], a = a*IdentityMatrix[Length[X]]];
  If[ VectorQ[a], a = DiagonalMatrix[a]];
  (* Note the - sign in the operator *)
  a = PiecewiseExpand[Piecewise[{{-a, True}}]];
  Inactive[Div][a.Inactive[Grad][T, X], X] + V - Q]

HeatTransferModelAxisymmetric[T_, {r_, z_}, k_, \[Rho]_, Cp_, 
  Velocity_, Source_] :=
 Module[{V, Q},
  V = If[Velocity === "NoFlow", 
    0, \[Rho]*Cp*Velocity.Inactive[Grad][T, {r, z}]];
  Q = If[Source === "NoSource", 0, Source];
  1/r*D[-k*r*D[T, r], r] + D[-k*D[T, z], z] + V - Q]

TimeHeatTransferModel[T_, TimeVar_, X_List, k_, \[Rho]_, Cp_, 
  Velocity_, Source_] := \[Rho]*Cp*D[T, {TimeVar, 1}] + 
  HeatTransferModel[T, X, k, \[Rho], Cp, Velocity, Source]

TimeHeatTransferModelAxisymmetric[T_, TimeVar_, {r_, z_}, k_, \[Rho]_,
   Cp_, Velocity_, Source_] :=
 \[Rho]*Cp*D[T, {TimeVar, 1}] + 
  HeatTransferModelAxisymmetric[T, {r, z}, k, \[Rho], Cp, Velocity, 
   Source]

Now, we can setup a verified heat equation operator for an axisymmetric case and solve your problem.

h = 10;
rmin = 2;
rmax = 4;
Subscript[\[CapitalGamma], temp] = 
  DirichletCondition[u[t, r, z] == 50, r == rmin];
nv = NeumannValue[0.1, r == rmax];
ic = {u[0, r, z] == 50 - ((50 - 10)/(4 - 2))*r};
\[CapitalOmega] = Rectangle[{rmin, 0}, {rmax, h}];
tend = 5;
parmop = TimeHeatTransferModelAxisymmetric[u[t, r, z], t, {r, z}, 
   k, \[Rho], Cp, "NoFlow", "NoSource"];
op = parmop /. {k -> 1, \[Rho] -> 1, Cp -> 1};
pde = {op == nv, Subscript[\[CapitalGamma], temp], ic};
ifun = NDSolveValue[pde, 
  u, {t, 0, tend}, {r, z} \[Element] \[CapitalOmega]]

The solution returns quickly. We can visualize the solution using the example in the verification manual or from the example on the Wolfram website shown here.

uRange = MinMax[ifun["ValuesOnGrid"]];
legendBar = 
  BarLegend[{"TemperatureMap", uRange}, 50, 
   LegendLabel -> Style["[\[Degree]C]", Opacity[0.6`]]];
options = {PlotRange -> uRange, 
   ColorFunction -> ColorData[{"TemperatureMap", uRange}], 
   ContourStyle -> Opacity[0.1`], ColorFunctionScaling -> False, 
   Contours -> 30, AspectRatio -> 1,
PlotPoints -> 41, FrameLabel -> {"r", "z"}, 
   PlotLabel -> Style["Temperature Field: u(t,r,z)", 18], 
   AspectRatio -> Automatic, ImageSize -> 250};
nframes = 80;
frames = Table[
   Legended[
    ContourPlot[ifun[t, r, z], {r, z} \[Element] \[CapitalOmega], 
     Evaluate[options]], legendBar], {t, 0, tend, tend/nframes}];
frames = (Rasterize[#1, "Image", ImageResolution -> 80] &) /@ frames;
ListAnimate[frames, SaveDefinitions -> True]
cpfn = ContourPlot[ifun[#, r, z], {r, z} \[Element] \[CapitalOmega], 
    ColorFunction -> "Temperature"] &;
frames2 = 
  Rasterize[#1, "Image", 
     ImageResolution -> 
      100] & /@ (Show[
       RegionPlot3D[
        rmin^2 <= x^2 + y^2 <= rmax^2 && 0 <= z <= h, {x, -rmax, 
         rmax}, {y, -rmax, rmax}, {z, 0, h}, Boxed -> False, 
        Axes -> False, PlotPoints -> 40, PlotStyle -> {Opacity[0.2]}, 
        Mesh -> False], 
       Graphics3D[{EdgeForm[Red], FaceForm[Gray], 
         GraphicsComplex[{{rmin, 0, 0}, {rmax, 0, 0}, {rmax, 0, 
            h}, {rmin, 0, h}}, {Texture[
            Show[cpfn[#], Frame -> False, PlotRangePadding -> None]], 
           Lighting -> {{"Ambient", White}}, 
           Polygon[{{1, 2, 3, 4}}, 
            VertexTextureCoordinates -> {{{0, 0}, {1, 0}, {1, 1}, {0, 
                1}}}]}]}], ImageSize -> 200] & /@ 
     Subdivide[0, tend, 80]);
ListAnimate[frames2, SaveDefinitions -> True]

Verification Animation

New in 10 Animation

Here is an example of how to make it look more 3D with a SliceContourPlot3D

cyl1 = Cylinder[{{0, 0, 0}, {0, 0, h}}, 2];
cyl2 = Cylinder[{{0, 0, 0}, {0, 0, h}}, 4];
cyl = RegionDifference[cyl2, cyl1];
frames3 = 
  Rasterize@
     SliceContourPlot3D[
      ifun[#, Sqrt[x^2 + y^2], z], {x, y, z} \[Element] cyl, 
      PlotRange -> uRange, 
      ColorFunction -> ColorData[{"TemperatureMap", uRange}], 
      ContourStyle -> Opacity[0.5`], ColorFunctionScaling -> False, 
      Contours -> 30, Boxed -> False, Axes -> False, 
      PlotPoints -> 40] & /@ Subdivide[0, tend, 40];
ListAnimate@frames3

Slice Plot 3D

$\endgroup$
3
  • $\begingroup$ (+1) But this is actually 2D, not 3D. How it could be solve in a case of 3D? $\endgroup$ – Alex Trounev Apr 25 '20 at 19:07
  • $\begingroup$ @AlexTrounev Thanks! Yes, it is an axisymmetric 2D model. The OP question didn't have any $\theta$ dependence, so why not take advantage of symmetry? I was having difficulty coming up with an advantage to the approach in the OP unless it is to avoid some simple post processing steps. Your answer revealed something subtle, but important. Increasing the number of boundary surfaces going from cartesian (4) to cylindrical (6), increases the risk of boundary-boundary artefacts. Those can be difficult to fix. You may need to control PointElement markers to fix. $\endgroup$ – Tim Laska Apr 25 '20 at 19:47
  • $\begingroup$ I am wonder that it working in 3D with periodic boundary condition. Usually I solve this problem in cartesian coordinates, then I have no headache with boundary elements and PeriodicBoundaryCondition[]. $\endgroup$ – Alex Trounev Apr 25 '20 at 22:54
4
$\begingroup$

In cylindrical coordinates region Cylinder[] should be transformed into Cuboid[], so effectively we should solve this problem in the region {r, 2, 4}, {\[Theta], 0, 2 Pi}, {z, 0, h} with a periodic boundary condition on $\theta$. Code

Needs["NDSolve`FEM`"]; h = 10; reg = 
 ImplicitRegion[
  2 <= r <= 4 && 0 <= \[Theta] <= 2 Pi && 0 <= z <= h, {r, \[Theta], 
   z}]; 

mesh = ToElementMesh[reg]

mesh["Wireframe"]
pbc = PeriodicBoundaryCondition[u[t, r, \[Theta], z], \[Theta] == 0, 
   TranslationTransform[{0, 2 Pi, 0}]];

ifun = NDSolveValue[{D[u[t, r, \[Theta], z], t] - 
     Laplacian[u[t, r, \[Theta], z], {r, \[Theta], z}, 
      "Cylindrical"] == NeumannValue[0., True], 
   DirichletCondition[u[t, r, \[Theta], z] == 50, 
    r < 2.1 && 0 < \[Theta] < 2 Pi], 
   u[1, r, \[Theta], z] == 50 - ((50 - 10)/(4 - 2))*r}, 
  u, {t, 1, 5}, {r, \[Theta], z} \[Element] mesh]

SliceContourPlot3D[
 ifun[4, r, \[Theta], 
  z], "ZStackedPlanes", {r, \[Theta], z} \[Element] mesh, 
 ColorFunction -> "TemperatureMap", Boxed -> False, 
 AxesLabel -> Automatic]

Projection into {x,y,z} coordinates:

SliceContourPlot3D[
     ifun[4, Sqrt[x^2 + y^2], ArcTan[x, y] + Pi, 
      z], "ZStackedPlanes", {x, -4, 4}, {y, -4, 4}, {z, 0, h}, 
     ColorFunction -> "TemperatureMap", Boxed -> False, 
     AxesLabel -> Automatic]

Figure 1

$\endgroup$
5
  • $\begingroup$ Should there be any $\theta$ dependence? I would have thought there would be a series of straight lines. Any thoughts? $\endgroup$ – Tim Laska Apr 25 '20 at 17:56
  • $\begingroup$ @TimLaska Probably we should increase number of points on $\theta $ or may be it just issues of PeriodicBoundaryCondition. But I see that you post a right model. $\endgroup$ – Alex Trounev Apr 25 '20 at 18:41
  • $\begingroup$ I think there is an issue where you have two types nodal boundary conditions on adjacent edges competing for the same corner point. Mathematica makes a hidden decision about which one takes precedence. Resolution will help mitigate the issue, but it would be nice to change the precedence so that you could experiment. $\endgroup$ – Tim Laska Apr 25 '20 at 18:56
  • $\begingroup$ @TimLaska It was issues of Method->Automatic. With FEM it looks better. $\endgroup$ – Alex Trounev Apr 25 '20 at 20:10
  • $\begingroup$ Indeed, that looks much better! $\endgroup$ – Tim Laska Apr 25 '20 at 20:30

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.