8
$\begingroup$

I try to model the Convection–diffusion equation for a rotating 3D annulus(small thickness):

enter image description here

The temperature u[t,x,y,z] is described in cartesian coordinates.

parameters

t1 = 2;
r1 = .161;r2 = .201;d = 0.02;\[CapitalDelta]\[CurlyPhi] = 20. Degree;
\[Rho]c = 2.419*^6;\[Lambda] = 172.;
p = 1.855 10^7;\[Alpha] = 90.;\[Omega] = 110.;T0 = 50;

The heat flux p (constant for reasons of simplification) is applied in a local range of the disc. Convection \[Alpha] is considered at the side z==d/2.

meshing

Because thickness is much smaller than the radii the automatic meshing leads to
huge element meshes(NDSolve can't start:::) I predefine a simple mesh (using MeshTools package package, thanks to @pinti))

Needs["NDSolve`FEM`"];
Needs["MeshTools`"]

mesh2D = AnnulusMesh[{0, 0}, {    r1 ,  r2 }, {0, 2 Pi}, {36, 10}];
mesh = ExtrudeMesh[mesh2D, d/2, 5];
scheibe = HexToTetrahedronMesh[mesh]

enter image description here

simulation

The system is simulated with MethodOfLines and FiniteElement because these methods allow the description of flux-boundaries using NeumannValue.

U = Monitor[
NDSolveValue[{ \[Rho]c ( Derivative[1, 0, 0, 0][u][t, x, y,z] + \[Omega] {-y, x, 0}.Grad[ u[t, x, y, z], {x, y, z}]) ==   \[Lambda]  Laplacian[u[t, x, y, z] , {x, y, z}] 
+ p NeumannValue[1, -x Tan[\[CapitalDelta]\[CurlyPhi]/2] <= y <= x Tan[\[CapitalDelta]\[CurlyPhi]/2] &&  z ==  d/2]
- \[Alpha]    NeumannValue[(u[t, x, y, z] - 0),z ==  d/2 && ! (-x Tan[\[CapitalDelta]\[CurlyPhi]/2] <= y <= x Tan[\[CapitalDelta]\[CurlyPhi]/2])   ]
, u[0, x, y, z] == T0},
u,{t, 0, t1}, Element[{x, y, z}, scheibe] 
, Method -> {"MethodOfLines", "TemporalVariable" -> t,"SpatialDiscretization" -> {"FiniteElement" }}, EvaluationMonitor :> (monitor= Row[{"t = ", CForm[Round[t, .01]]}])] , monitor]

Results of simulation

Show[{SliceDensityPlot3D[U[t, x, y, d/2], {"XStackedPlanes", Subdivide[0,t1, 5]},  {t, 0,t1}, {x, -r2, r2}, {y, -r2, r2},PlotRange -> {0, All},RegionFunction -> Function[{t, x, y}, r1^2 <= x^2 + y^2 <= r2^2],AxesLabel -> {Zeit,None,None}, Ticks -> {Automatic, None, None},ColorFunction -> (ColorData["TemperatureMap"][#] &), PlotLegends ->BarLegend[{(ColorData["TemperatureMap"][#] &), {0, 200}}]]}]

enter image description here

seem to be ok, but NDSolve gives a warning

NDSolveValue::femcsp: The computed Peclet number is 515.3080637051237 and is larger than the mesh order (1), and the result may not be stable. Adding artificial diffusion may help.`

Checking temperature over time (smooth curve expected)

Plot[U[t, 0, (r1 + r2)/2, d/2], {t, 0, t1},AxesLabel -> {time, temperature}]

enter image description here

indicates numerical problems.

My questions

How to avoid message peclet-number? Using automated meshgeneration inside NDSolve leads to huge mesh( NDSolve doesn't run)

How to force smooth time response?

Thanks!

$\endgroup$
  • $\begingroup$ @user21 Thanks for editing my question. I'm curious about your thoughts about this problem. $\endgroup$ – Ulrich Neumann Nov 25 '19 at 9:25
  • $\begingroup$ I had a brief look but could not come up with something. I'll hope to find some time to put it on a big machine to see if that would actually solve the issue. But might take some time until I get to that. $\endgroup$ – user21 Nov 27 '19 at 16:43
  • 1
    $\begingroup$ @user21 Thanks, I'll await your contribution! $\endgroup$ – Ulrich Neumann Nov 27 '19 at 17:36
  • 1
    $\begingroup$ Should these not overlap? Show[Region[ ImplicitRegion[-x Tan[\[CapitalDelta]\[CurlyPhi]/2] <= y <= x Tan[\[CapitalDelta]\[CurlyPhi]/2], {x, y}]], Region[RegionDifference[Disk[{0, 0}, r2], Disk[{0, 0}, r1]]]] $\endgroup$ – user21 Nov 28 '19 at 14:29
  • 1
    $\begingroup$ @user21 Yes the regions must be the same. Change the first to Region[ImplicitRegion[-x Tan[\[CapitalDelta]\[CurlyPhi]/2] <= y <=x Tan[\[CapitalDelta]\[CurlyPhi]/2] && r1^2 <= x^2 + y^2 <= r2^2, {x, y}]] $\endgroup$ – Ulrich Neumann Nov 28 '19 at 15:38
5
$\begingroup$

We can add artificial viscosity to correct the imperfections of the mesh. But for this problem, the artificial viscosity is 632 times greater than $\lambda$. The second option is to change the model - go to the reference system where the disk is stationary and the heat source is moving. In this case, there are also messages regarding the calculation of a compiled function. But this does not spoil the numerical solution.

t1 = 2;
r1 = .161; r2 = .201; d = 0.02; \[CapitalDelta]\[CurlyPhi] = 
 20. Pi/180;
\[Rho]c = 2.419*^6; \[Lambda] = 172.;
p = 1.855 10^7; \[Alpha] = 90.; \[Omega] = 110.; T0 = 50.;
Needs["NDSolve`FEM`"];
Needs["MeshTools`"]

mesh2D = AnnulusMesh[{0, 0}, {r1, r2}, {0, 2 Pi}, {36, 10}];
mesh = ExtrudeMesh[mesh2D, d/2, 5];
scheibe = HexToTetrahedronMesh[mesh]

U = NDSolveValue[{\[Rho]c Derivative[1, 0, 0, 0][u][t, x, y, 
       z] == \[Lambda] Laplacian[u[t, x, y, z], {x, y, z}] + 
      p NeumannValue[
        1, -x Tan[\[CapitalDelta]\[CurlyPhi]/2 - \[Omega] t] <= y <= 
          x Tan[\[CapitalDelta]\[CurlyPhi]/2 - \[Omega] t] && 
         z == d/2] -  \[Alpha] NeumannValue[(u[t, x, y, z] - 0.), 
        z == d/2 && ! (-x Tan[\[CapitalDelta]\[CurlyPhi]/
                2 - \[Omega] t] <= y <= 
            x Tan[\[CapitalDelta]\[CurlyPhi]/2 - \[Omega] t])], 
    u[0, x, y, z] == T0}, u, {t, 0, t1}, Element[{x, y, z}, scheibe], 
   Method -> {"MethodOfLines", "TemporalVariable" -> t, 
     "SpatialDiscretization" -> {"FiniteElement"}}] // Quiet

The result is different from what it was in the original version, but the solution here is smooth.

Show[{SliceDensityPlot3D[
   U[t, x, y, d/2], {"XStackedPlanes", Subdivide[0, t1, 5]}, {t, 0, 
    t1}, {x, -r2, r2}, {y, -r2, r2}, PlotRange -> {0, All}, 
   RegionFunction -> Function[{t, x, y}, r1^2 <= x^2 + y^2 <= r2^2], 
   AxesLabel -> {Zeit, None, None}, Ticks -> {Automatic, None, None}, 
   ColorFunction -> (ColorData["TemperatureMap"][#] &), 
   PlotLegends -> 
    BarLegend[{(ColorData["TemperatureMap"][#] &), {0, 2000}}]]}]

Plot[U[t, 0, (r1 + r2)/2, d/2], {t, 0, t1}, 
 AxesLabel -> {time, temperature}, PlotRange -> All]

Figure 1

It is necessary to normalize p and $\alpha$ to the real area. Then the scales coincide with Ulrich Neumann result.

p = 1.855 10^7 20/360; \[Alpha] = 90.*340/360;

Figure 2

| improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ Comments are not for extended discussion; this conversation has been moved to chat. $\endgroup$ – halirutan Nov 30 '19 at 17:18

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.