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

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]


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) \
\[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.

• You can solve the equation in the case of concentric circles, and map the solution to the case you want with a bilinear transformation: math.stackexchange.com/questions/1197630/…
– mjw
May 31 at 2:10
• There is an entire section in the documentation on that, have you seen that? May 31 at 4:45
• It might also simplify you life if you use a HeatTransferPDEComponent May 31 at 5:09

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["NDSolveFEM"]

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) \

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}] • Thank you kindly for the detailed answer. Jun 2 at 8:51
• @ConfuzzledStudent You are welcome! Jun 2 at 8:52

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.

• Not my down vote, but this does not really answer the question. You could add a note at the beginning of your "answer" saying that this is an extended comment and not an answer. Just a thought. May 31 at 7:09
• What is defined in this community as "really answering" to an obviously wrong formula labeled as the "Heat Equation in3D"? Let's wait. May 31 at 11:56