2
$\begingroup$

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 : centered boundary condition result

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

enter image description 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.

$\endgroup$
3
  • 2
    $\begingroup$ 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/… $\endgroup$
    – mjw
    May 31, 2023 at 2:10
  • $\begingroup$ There is an entire section in the documentation on that, have you seen that? $\endgroup$
    – user21
    May 31, 2023 at 4:45
  • $\begingroup$ It might also simplify you life if you use a HeatTransferPDEComponent $\endgroup$
    – user21
    May 31, 2023 at 5:09

2 Answers 2

2
$\begingroup$

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}] 

Figure 1

$\endgroup$
2
  • $\begingroup$ Thank you kindly for the detailed answer. $\endgroup$ Jun 2, 2023 at 8:51
  • $\begingroup$ @ConfuzzledStudent You are welcome! $\endgroup$ Jun 2, 2023 at 8:52
0
$\begingroup$

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

tradtional Form of eq

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.

$\endgroup$
4
  • $\begingroup$ 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. $\endgroup$
    – user21
    May 31, 2023 at 7:09
  • $\begingroup$ What is defined in this community as "really answering" to an obviously wrong formula labeled as the "Heat Equation in3D"? Let's wait. $\endgroup$
    – Roland F
    May 31, 2023 at 11:56
  • $\begingroup$ This answer is useful! Adding a positive vote. $\endgroup$
    – mjw
    May 31, 2023 at 12:17
  • $\begingroup$ @Roland F The choice of tag was perhaps wrongly done, sorry about that. What I wanted to truly ask was how to implement off-center bc with NDsolve to solve an equation that just so happens to ressemble a heat equation. $\endgroup$ May 31, 2023 at 13:31

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.