# Transient Heatflow in a Wheel Disc Brake

I try to model the Convection–diffusion equation for a rotating 3D annulus(small thickness): 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["NDSolveFEM"];
Needs["MeshTools"]

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


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

• 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. Nov 27 '19 at 16:43
• @user21 Thanks, I'll await your contribution! Nov 27 '19 at 17:36
• 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]]]] Nov 28 '19 at 14:29
• @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}]] Nov 28 '19 at 15:38

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["NDSolveFEM"];
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] 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;
` • Comments are not for extended discussion; this conversation has been moved to chat. Nov 30 '19 at 17:18