Discontinuity problem with 3D cylindrical heat equation (possibly due to a conversion between Cartesian to Cylindrical coordinates)

Question

I have been working a certain type of 3D cylindrical diffusion equation for a bit now. I am trying to simulate a longitudinal diffusion process in a cylinder with a dislocation defect that will make the heat flow spiral upward as you can see underneath with these two animations :

Here you can find the code I used to simulate an upward diffusion process along a cylinder.

(*parameters*)
r0 = 0.5;
r1 = 1;
h = 10;
\[Beta] = 1;


(*arg[z_]/;Im[z]<0:=Arg[z]+2 \
Pi;arg[z_]/;Im[z]\[GreaterEqual]0:=Arg[z]  possible alternative to \
Mod[Artan[x,y],2Pi]*)

(*region definition*)
reg = Cuboid[{r0, 0., 0.}, {r1, 2 Pi, h}];

reg3D = ImplicitRegion[
   r0^2 <= x^2 + y^2 <= 1 && 0 <= z <= h, {x, y, z}];

(*equation + conditions*)
eq3 = D[u[t, r, \[Theta], z], 
    t] - (D[u[t, r, \[Theta], z], r, r] + (1/r)*
      D[u[t, r, \[Theta], z], r] + 
     1/(r^2)*D[
       u[t, r, \[Theta], 
        z], \[Theta], \[Theta]] + (r^2 + \[Beta]^2)/(r^2)*
      D[u[t, r, \[Theta], z], z, z] - (2 \[Beta]^2/r^2) D[
       u[t, r, \[Theta], z], z, \[Theta]]);


ic = u[0, r, \[Theta], z] == 1;
bc = DirichletCondition[u[t, r, \[Theta], z] == 0, z == 0];
nV = NeumannValue[0, r == r1] + NeumannValue[0, r == r0] + 
   NeumannValue[1, z == h];

(*solution computation*)
sol = NDSolveValue[{eq3 == nV, ic, bc}, 
   u, {t, 0, 5}, {r, \[Theta], z} \[Element] reg];

frames = Table[
   DensityPlot3D[
    sol[t, Sqrt[x^2 + y^2], 
     Mod[Arctan[x, y], 2 Pi],z], {x, y, z} \[Element] reg3D, 
    ColorFunction -> "Rainbow", OpacityFunction -> None, 
    Boxed -> False, Axes -> False, PlotRange -> All, PlotPoints -> 50,
     PlotLabel -> Row[{"t = ", t}], 
    ColorFunctionScaling -> False], {t, .05, 5, .05}];

Export["better view all.mov", frames]

As you can see with these animations, there seem to be a discontinuity on the left side of the cylinder. This discontinuity left me quite puzzled as I could not find any physical explanation. The equation itself does not hint at any kind of irregularities with respect to theta. The solution should be continuous around the cylinder, especially if we imagine that the solution represents temperature for instance, then such behavior would be prohibited.

After discussing this with some other people, I remembered a very simple fact that I had totally missed. The discontinuity I was trying hard to explain could simply come from the Arctan[x,y] function I was using when I tried to animate the whole diffusion process.

I tried using other alternatives, for instance, I tried using Mod[Arctan[y/x],2Pi] instead but it gave me the following animation.

At least I now know that it's the conversion between cartesian and polar coordinates that gives me such a peculiar discontinuity.

I also tried using the Arg function and even tried building one that goes from 0 to 2Pi following the answers found here. However it pretty much gave me the same result as with Mod[Arctan[x,y],2Pi].

I searched for some similar questions asked here regarding the Arctan function and found this question here showing us how to build a version of Arctan that could potentially help me here. However it was fruitless in the end as the simulation failed this time and gave me a blank video showing a white background.

Eventually I was very interested in how the discontinuity in the Arctan function was bypassed in this post by using the derivative of said function and then integrating when needed, however I am still working on that lead (I'm still a beginner when it comes to Mathematica alright...).

In the end, I'm fairly sure that it isn't a "physical discontinuity" but instead is induced by my own algorithm but I would like to know your opinion on that matter as I could very well be wrong from the start. And if so, is there a way for me to possibly get a nice simulation without such a discontinuity ?

I would really appreciate the help, even just tiny hints. Thanks in advance.

Edit : Thanks the advice of user xzczd, I revised my code so that It contains periodic boundary conditions, following the example given here using double pbcs and ghost vicinity. I now have the following code :

(*parameters*)
r0 = 0.5;
r1 = 1;
h = 5;
\[Beta] = .5;

(*Ghost vivinity parameter*)
epsilon = 0;


(*region definition*)
reg = Cuboid[{r0, -2 Pi - epsilon, 0.}, {r1, 2 Pi + epsilon, h}]; 

reg3D = ImplicitRegion[
   r0^2 <= x^2 + y^2 <= r1 && 0 <= z <= h, {x, y, z}];

(*Diffusion equation*)
eq = D[u[t, r, \[Theta], z], 
    t] - (D[u[t, r, \[Theta], z], r, r] + (1/r)*
      D[u[t, r, \[Theta], z], r] + 
     1/(r^2)*D[
       u[t, r, \[Theta], 
        z], \[Theta], \[Theta]] + (r^2 + \[Beta]^2)/(r^2)*
      D[u[t, r, \[Theta], z], z, z] - (2 \[Beta]^2/r^2) D[
       u[t, r, \[Theta], z], z, \[Theta]]);

(*Initial and boundary conditions*)
ic = u[0, r, \[Theta], z] == 1;
bc = DirichletCondition[u[t, r, \[Theta], z] == 0, 
   z == 0];
pbc={PeriodicBoundaryCondition[u[t,r,\[Theta],z],\[Theta]\
\[Equal]2Pi+epsilon&&0<z<h,TranslationTransform[{0,-2Pi,0}]\
],PeriodicBoundaryCondition[u[t,r,\[Theta],z],\[Theta]\[Equal]0-\
epsilon&&0<z<h,TranslationTransform[{0,2Pi,2Pi*\[Beta]0e[0, r == r0] + 
   NeumannValue[1, z == h];

(*solution computation*)
sol = NDSolveValue[{eq == nV, ic, bc(*,pbc*)}, 
   u, {t, 0, 15}, {r, \[Theta], z} \[Element] reg];


(*Animation of the diffusion process*)
frames = Table[
   DensityPlot3D[
    sol[t, Sqrt[x^2 + y^2], Mod[ArcTan[x, y], 2 Pi], 
     z], {x, y, z} \[Element] reg3D, ColorFunction -> "Rainbow", 
    OpacityFunction -> None, Boxed -> False, Axes -> False, 
    PlotRange -> All, PlotPoints -> 50, PlotLabel -> Row[{"t = ", t}],
     ColorFunctionScaling -> False], {t, .05, 15, .05}];

Export["beta=1.mov", frames]

This new code allows me to get a nice continuous animation, at the expense of the spiraling motion, however, which completely defeats the purpose...

I really don't know what to do now.

Answer 1

The end goal is to create a visualization of heat flow that spirals in a helical fashion due to a dislocation defect. The OP's initial attempt was to formulate the problem and cylindrical coordinates, which led to difficulty achieving the desired visual effects.

This answer proposes a Cartesian coordinate approach that may be irrelevant, but it will produce a spiraling heat diffusion effect.

Bilayered screw

I prototyped the model in SolidWorks by sweeping two rectangular slabs along a helical path as shown below.

If we step through a Z clip plane, we can see two semi disks spiraling. This is simple enough geometry that we should be able to construct a solvable mesh.

Helper functions

In addition to the FEM package, we will use MeshTools to extrude a quadrilateral disk mesh. We will apply a small twist to the mesh as we move up the Z direction.

(*Import required FEM package*)
Needs["NDSolve`FEM`"];
(*Install MeshTools*)
(*Uncomment if not installed*)
(*ResourceFunction["GitHubInstall"]["c3m-labs","MeshTools"]*)
Needs["MeshTools`"]
(*Wrap twisted mesh builder function in module*)
Clear[buildMeshes]
buildMeshes[rad_, nrad_, ht_, θ_, turns_] := 
 Module[{nelm, mesh2D, mesh3D, Ω3Dright, rmf, 
   regmarkerfn, inc, mean,
   regmarkers, newcrd, mesh, vis2D, vis3D, visTwisted},
  nelm = 360 ° turns/θ;
  mesh2D = DiskMesh[{0, 0}, rad, nrad];
  vis2D = Show[
    mesh2D["Wireframe"],
    Graphics[{Point[mesh2D["Coordinates"]]}]
    ];
  mesh3D = ExtrudeMesh[mesh2D, nelm, nelm];
  (*Right RegionMember Function*)
  Ω3Dright = Cuboid[{0, -rad, 0}, {rad, rad, nelm}];
  rmf = RegionMember[Ω3Dright];
  regmarkerfn = If[rmf[#], reg["right"], reg["left"]] &;
  (*Get mean coordinate of each hexa for region marker assignment*)
  inc = ElementIncidents[mesh3D["MeshElements"]];
  mean = Mean /@ GetElementCoordinates[mesh3D["Coordinates"], #] & /@ 
     inc // First;
  regmarkers = regmarkerfn /@ mean;
  inc = inc[[1]];
  (*Create element mesh*)
  mesh3D = 
   ToElementMesh["Coordinates" -> mesh3D["Coordinates"], 
    "MeshElements" -> {HexahedronElement[inc, regmarkers]}];
  vis3D = mesh3D[
     "Wireframe"["MeshElement" -> "MeshElements", 
       "MeshElementStyle" -> (Directive[Opacity[0.5], 
          FaceForm[#]] &  /@ {Red, White}), 
     ViewPoint -> {-1.5, 0.8, -3}, ViewVertical -> {0, 1, 0}, 
       PlotRange -> All]];
  newcrd = 
   mesh3D["Coordinates"] /. {{x_, y_, 
       z_} -> {x Cos[z θ] - y Sin[z θ], 
       y Cos[z θ] + x Sin[z θ], z/nelm ht}};
  (*Create twisted mesh*)
  mesh = ToElementMesh["Coordinates" -> newcrd, 
    "MeshElements" -> {HexahedronElement[inc, regmarkers]}];
  visTwisted = mesh[
     "Wireframe"["MeshElement" -> "MeshElements", 
       "MeshElementStyle" -> (Directive[Opacity[1], FaceForm[#](*, 
                 EdgeForm[]*)] &  /@ {Red, White}), 
     ViewPoint -> {2.5, -1.75, -1.7}, 
     ViewVertical -> {-0.2, 0.13, -0.97},
       PlotRange -> All]];
  <|"mesh2D" -> mesh2D, "mesh3D" -> mesh3D, "twistedMesh" -> mesh,
   "vis2D" -> vis2D, "vis3D" -> vis3D, "visTwisted" -> visTwisted|>
  ]

A small coarse test example

We will test the buildMeshes function on a single turn with a large twist angle to produce a candy cane type structure.

(*Association for Clearer Region Assignment*)
reg = <|"left" -> 1, "right" -> 2|>;
(*Set parameters*)
rad = 1;
ht = 1;
θ = 20 °;
turns = 1;
nrad = 6;
meshes = buildMeshes[rad, nrad, ht, θ, turns]
meshes["visTwisted"]

Larger model

Now, we will build a five-turn model with finer discretization both in the radial and theta direction.

(*Set parameters*)
rad = 1;
ht = 2;
θ = 5 °;
turns = 5;
nrad = 12;
meshes = buildMeshes[rad, nrad, ht, θ, turns];
mesh = meshes["twistedMesh"];

Heat transfer PDE system

The twistedMesh function produced a mesh that had element markers so that we can assign an element marker-based thermal conductivity. In this case, we will make one of the materials 10 times less conductive than the other. We will use the HeatTransferPDEComponent (new in 12.2) to set up the PDE system. We will apply a 1° temperature difference across the actual boundaries. For expediency, we will consider steady-state only.

(*Set material parameters*)
kleft = 1;
kright = 1/10;
k = 
  Evaluate[
   Piecewise[
    {
     {kleft, ElementMarker == reg["left"]},
      {kright, True}
     }]
   ];
(*Set up PDE system*)
vars = {Θ[x, y, z], {x, y, z}};
pars = <|"ThermalConductivity" -> {{k, 0, 0}, {0, k, 0}, {0, 0, k}}|>;
Γtop = 
  HeatTemperatureCondition[z == ht, vars, 
   pars, <|"SurfaceTemperature" -> 1|>];
Γbot = 
  HeatTemperatureCondition[z == 0, vars, 
   pars, <|"SurfaceTemperature" -> 0|>];
eqn = {HeatTransferPDEComponent[vars, pars] == 
   0, Γtop, Γbot}

Solving and postprocessing

Before we solve the PDE system we will first inspect the material conductivity along the z-axis.

(*Create thermal conductivity interpolation function*)
Short[values = 
  Function[{x, y, z}, Piecewise[{{(kleft + kright)/2, Abs[x] < 0.001},
      {kleft, x < 0},
      {kright, True}}]] @@@ meshes["mesh3D"]["Coordinates"]]
kif = ElementMeshInterpolation[{mesh}, values];
framesK = 
  SliceContourPlot3D[
     kif[x, y, z], {z == # ht}, {x, y, z} ∈ mesh, 
     ViewPoint -> {Infinity, Infinity, Infinity}, 
     ColorFunction -> "ThermometerColors", PlotLegends -> Automatic, 
     Contours -> 13] & /@ (Subdivide[0, 1, 50]~Join~
     Subdivide[1, 0, 50]);
ListAnimate@framesK

Now, we will look at the temperature solution. Since this is a steady-state solution, the temperature gradients are quite small so we will rescale the temperature at each clipping plane.

(*Solve and plot solution*)
Tfun = NDSolveValue[eqn, Θ, {x, y, z} ∈ mesh];
frames = SliceContourPlot3D[
     Tfun[x, y, z], {z == # ht}, {x, y, z} ∈ mesh, 
     ViewPoint -> {Infinity, Infinity, Infinity}, 
     ColorFunction -> "ThermometerColors", PlotLegends -> Automatic, 
     Contours -> 13] & /@ (Subdivide[0, 1, 50]~Join~
     Subdivide[1, 0, 50]);
ListAnimate@frames

It will be interesting to look at the heat flux along the axis. To do so, it will be necessary to create a derived variable for the temperature gradient and multiply it by our previously derived variable for thermal conductivity.

(*Set up derived temperature gradient and flux variable*)
Clear[gradT, flux]
gradT[x_, y_, z_] = {Derivative[1, 0, 0][Tfun][x, y, z], 
   Derivative[0, 1, 0][Tfun][x, y, z], 
   Derivative[0, 0, 1][Tfun][x, y, z]};
flux[x_, y_, z_] = -kif[x, y, z] gradT[x, y, z];
frames2 = 
  Show[SliceVectorPlot3D[
      flux[x, y, z], {z == # ht}, {x, y, z} ∈ mesh, 
      PlotLabel -> N[# ht], Axes -> False, Boxed -> False, 
      ViewPoint -> {Infinity, Infinity, -Infinity}, 
      VectorScaling -> Automatic, VectorSizes -> {0.1, 2}, 
      VectorMarkers -> Placed["Arrow3D", "Start"], 
      VectorAspectRatio -> 0.15, VectorPoints -> 25],
     SliceContourPlot3D[
      kif[x, y, z], {z == # ht}, {x, y, z} ∈ mesh, 
      ColorFunction -> "ThermometerColors", 
      Contours -> 13]] & /@ (Subdivide[0, 1, 50]~Join~
     Subdivide[1, 0, 50]);
ListAnimate@frames2

The heat flux is predominantly in the downward direction as expected. Maybe less expected, is the highest heat fluxes are confined to the interface and fall off rather quickly.