Solving 3D heat equation with an off-center boundary condition

Question

So i have this code (albeit a simplified version but it'll do for this question)which solves a time dependant 3D heat equation on a cylinder.

r1 = 0.001;
r2 = 1;
h = 1;
m = 1.5; (*some parameter*)
reg3D = ImplicitRegion[
   r1^2 <= x^2 + y^2 <= r2^2 && 0 <= z <= h, {x, y, z}];


eq = D[T[t, r, \[Theta], z], t] - D[T[t, r, \[Theta], z], z, z] - 
   D[T[t, r, \[Theta], z], \[Theta]] - 
   1 /r^2*D[T[t, r, \[Theta], z], \[Theta], \[Theta]] - 
   1/r*D[T[t, r, \[Theta], z], r] - 
   1/(m^2* r) *D[T[t, r, \[Theta], z], r, \[Theta]] - (1 + m^2)/(
    2 *m^2)*D[T[t, r, \[Theta], z], r, r] ;

(*initial and boundary conditions*)
ic = T[0, r, \[Theta] , z] == 0;
bc = DirichletCondition[T[t, r, \[Theta], z] == 1, 
   0 < \[Theta] < 2 Pi && 0.1 <= r <= 0.5];
pbc = PeriodicBoundaryCondition[T[t, r, \[Theta], z], \[Theta] == 0 , 
   TranslationTransform[{0, 2 Pi, 0}]];

sol = NDSolveValue[{eq == 0, ic, bc, pbc}, 
  T, {r, 0.001, 1.}, {\[Theta], 0, 2 *Pi}, {z, 0, h}, {t, 0, 10}]

When I plot the result with the following code,

(*3d plot*)
SliceDensityPlot3D[
 sol[1, Sqrt[x^2 + y^2], Mod[ArcTan[x, y], 2 Pi], 
  z], {x^2 + y^2 == 0.98*r2^2, x^2 + y^2 == 1.01*r1^2, z == 0, 
  z == h}, {x, y, z} \[Element] reg3D, 
 ColorFunctionScaling  -> False, ColorFunction -> "Rainbow", 
 Boxed -> False, Axes -> False, PlotPoints -> {50, 50, 10}, 
 ViewPoint -> Above]

I get this :

Pretty much the result anyone could have expected given the look of our equation. However now I would like to "heat" the cylinder through an off-center circle, much like the following figure, where i want my circle in the middle to be either on the left or the right of the origin (on the right here):

However I can't seem to be able to do it at all. None of the solutions I tried worked. For instance I looked at the following post here and improvised a modified b.c. like that :

(*offcenter parameter*)
a = 0.5;
r0 = 0.5; (*circle centered on polar coordinates (r,\[Phi])=(r0,0) \
and of radius a=0.5*)
\[Alpha] = ArcSin[a/r0];
(*initial and boundary conditions*)
ic = T[0, r, \[Theta] , z] == 0;
bc = DirichletCondition[
   T[t, r, \[Theta], z] == 
    1, (0 < \[Theta] <= \[Alpha] || -\[Alpha] <= \[Theta] < 0) && 
    r <= r0*Cos[\[Theta]] + (a^2 - r0^2*Sin[\[Theta]])^0.5];

(I'm not even sure this might even be the correct equation for an off-center disk. I've seen lots of stuff for off center circles, but nothing on disks ... So i just tried putting a <= and see where that leads...)

Well, that didn't work at all : NDSolveValue::ndsz: At t == 0.07406818267524817`, step size is effectively zero; singularity or stiff system suspected.

So here i am. Anyone has an idea on how to achieve what i would like ? I would greatly appreciate it.

Have a good day.

Answer 1

Note, that periodic boundary conditions and Dirichlet type condition could not be apply in one border simultaneously as it shown in a picture above. Therefore we should replace of center b.c. from the line $\theta=0$. Also we need to restrict $a, r_0$ so that Dirichlet condition region is not crossing line $\theta =0$. With all this restriction we have

Needs["NDSolve`FEM`"]

r1 = 0.0;
r2 = 1;
h = 1;
m = 1.5; (*some parameter*)
reg = Cuboid[{0, 0, 0}, {r2, 2 Pi, h}]; mesh = 
 ToElementMesh[reg, MaxCellMeasure -> 1/2000]



(*offcenter parameter*)a = 0.4;
r0 = -0.5; (*circle centered on polar coordinates (r,\[Phi])=(r0,0) \
and of radius a=0.5*)

eq = D[T[t, r, \[Theta], z], t] - D[T[t, r, \[Theta], z], z, z] - 
   D[T[t, r, \[Theta], z], \[Theta]] - 
   1/r^2*D[T[t, r, \[Theta], z], \[Theta], \[Theta]] - 
   1/r*D[T[t, r, \[Theta], z], r] - 
   1/(m^2*r)*D[T[t, r, \[Theta], z], r, \[Theta]] - (1 + m^2)/(2*m^2)*
    D[T[t, r, \[Theta], z], r, r];

(*initial and boundary conditions*)
ic = T[0, r, \[Theta], z] == 0;
bc = DirichletCondition[T[t, r, \[Theta], z] == 1, 
   r^2 - 2 r r0 Cos[\[Theta]] + r0^2 <= a^2 && z == 1];
pbc = PeriodicBoundaryCondition[T[t, r, \[Theta], z], \[Theta] == 0, 
   TranslationTransform[{0, 2 Pi, 0}]];

sol = NDSolveValue[{eq == 0, ic, bc, pbc}, T, {t, 0, 2}, 
  Element[{r, \[Theta], z}, mesh]] 

We have message from NDSolve

NDSolveValue::femcsp: The computed Peclet number is 8.22826911771027` and is larger than the mesh order (2), and the result may not be stable. Refining the mesh or adding artificial diffusion may help.

Visualization

Table[DensityPlot[
  sol[2, Sqrt[x^2 + y^2], ArcTan[x, y] + Pi, z], {x, -1, 1}, {y, -1, 
   1}, ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
  PlotRange -> {0, 1}, Exclusions -> None, PlotPoints -> 100, 
  PlotLabel -> Row[{"z=", z}], Frame -> False], {z, 0, 1, .2}] 

Answer 2

It would be easier to decipher, if the equations are given conventionlly. What I see is

You can easily copy the cell by a right click, copy as LaTeX and paste it here between four $ \\\$\\\$ \ \ \\\$\\\$ $

$$\partial_t T =\frac{m^2+1}{2 m^2} \ \partial_{r,r} T + \frac{1}{r}\partial_r T + \frac{1}{r^2}\partial_{\theta \theta } T + \partial_{z,z} T + \partial_{r,r} T +\partial_\theta T + \frac{1}{m^2 r} \partial_{r,\theta} T$$

This is not the heat equation in cylinder coordinates at rest, at least not by the last two terms.