Difficulty setting a time step or stabilizing the transient heat equation

Question

I have been trying to solve the transient heat equation using finite elements. When I set the thermal diffusivity to a realistic number I get instabilities in the solution. It seems that larger diffusivities work but smaller ones don't. This is backward from what one would typically expect with numerical solutions to the heat equation. I don't know if Mathematica is trying an implicit time integrator and I don't seem to be able to set one. Adding artificial diffusion at the level necessary to stabilize the problem is completely non-physical. I eventually want to create regions with different diffusivities in the domain, which is why there are extra lines in the code below.

boundarylayerthickness = 10^-5;
traythickness = 1*10^-2;
traywidth = 300*10^-3;
traydepth = 200*10^-3;
trayalpha = 14*10^-7;
oilalpha = 10^-7;
oilthickness = 10^-3;
objectthickness = 2*10^-2;
objectwidth = 100*10^-3;
objectdepth = 70*10^-3;
tray = ImplicitRegion[
x >= 0 && x <= traywidth && y >= 0 && y <= traydepth && z >= 0 && 
z <= traythickness, {x, y, z}];
oil = ImplicitRegion[
x >= 0 && x <= traywidth && y >= 0 && y <= traydepth && 
z > traythickness && z <= traythickness + oilthickness, {x, y, 
z}]; 
object = ImplicitRegion[
x >= (traywidth - objectwidth)/2 && 
x <= (traywidth + objectwidth)/2 && 
y >= (traydepth - objectdepth)/2 && 
y <= (traydepth + objectdepth)/2 && 
z > traythickness + oilthickness && 
z <= traythickness + oilthickness + objectthickness, {x, y, z}]; 
gammad[t_] = {DirichletCondition[u[t, x, y, z] == 200.*t, 
z == 0 && x >= 0 && x <= traywidth && y >= 0 && y <= traydepth]};
omega = RegionUnion[tray, oil, object];
uinit = NDSolveValue[{10^-7 \!\(
\*SubsuperscriptBox[\(\[Del]\), \({x, y, z}\), \(2\)]\(u[0, x, y, 
    z]\)\) == 0, gammad[0]}, 
  u[0, x, y, z], {x, y, z} \[Element] omega, 
  Method -> {"PDEDiscretization" -> {"FiniteElement", 
          "MeshOptions" -> {"BoundaryMeshGenerator" -> 
                "Continuation"}}}];
uif = NDSolveValue[{D[u[t, x, y, z], t] - 10^-7 \!\(
\*SubsuperscriptBox[\(\[Del]\), \({x, y, z}\), \(2\)]\(u[t, x, y, 
     z]\)\) == 0, gammad[t], u[0, x, y, z] == uinit}, 
u, {t, 0, 25}, {x, y, z} \[Element] omega];
SliceContourPlot3D[
uif[2, x, y, z], {x, 0, traywidth}, {y, 0, traydepth}, {z, 0, 
traythickness + oilthickness + objectthickness}, 
PlotTheme -> "Detailed"]

This comment doesn't concern the problem you have at the moment, but I prefer to warn you about another problem you will have when the diffusivities will be different : Be carefull of the continuity of flux at the interface between the materials. see this — andre314, Sep 23 '17 at 20:13
Thank you. I know that is going to be an issue. There is some Mathematica documentation on the boundary surfaces. I am hoping that I'll be able to follow the examples there. — tensor, Sep 23 '17 at 22:10
Sorry, I didn't notice the link. That will be really helpful. — tensor, Sep 24 '17 at 02:15
here is a link about some further considerations about the problem of flux — andre314, Sep 24 '17 at 12:44
The discontinuities are a tricky problem and I know almost nothing about FEM. When I have worked with discontinuities it has been in the context of 1D finite volume methods, which are of no use for my current problem with its strange geometries. — tensor, Sep 25 '17 at 21:16

Difficulty setting a time step or stabilizing the transient heat equation

0 Answers0