# Eliminating instabilities in a transient finite element solution at a discontinuity near t = 0

I have found a transient solution to the 1D heat equation where the initial condition is discontinuous. The results are accurate except for when time is small. The initial condition looks like this:

$\qquad$ Initial Condition

After solving with NDSolveValue using "PDEDiscretization" the solution at time 0 looks like this:

$\qquad$ Solution at Time t = 0

There is clearly some instability near the discontinuity that is resulting in inaccurate evaluation of the temperature near the simulated boundary. The solution stabalizes with very little time such that the solution at time 0.01 looks like this:

$\qquad$ Solution at Time t= 0.01

I am happy with the long term results but for my purposes i need an accurate solution within an extremely small time-frame (micro seconds, t ~ .000001). Can anyone please suggest ways to eliminate these issues that are being caused by the discontinuity?

<< "NDSolveFEM"

ts = .000250; tsl = .000250; (* m *)
ρs = 3980;  ρsl = 958; (* kg/m3 *)
ks = .035;  ksl = .00067; (* kW/m/K *)
cs = .75; csl = 4.22; (* kJ/kg/K *)

ρ = If[0 <= z < ts, ρs, ρsl];
k = If[0 <= z < ts, ks, ksl];
c =  If[0 <= z < ts, cs, csl];

eqn1 = k*D[T1[z], z] + 28;
Tbl = 100;
Subscript[Γ1, D] = DirichletCondition[T1[z] == Tbl, z == ts + tsl];

BCr = NDSolveValue[{eqn1 == 0, Subscript[Γ1, D]}, T1, {z, 0, ts + tsl}];

Plot[BCr[z], {z, 0, ts + tsl}, GridLines -> {{ts}, {0}}]

Ti[z_] :=  Piecewise[{{BCr[z], 0 <= z < ts}, {100, True}}];
Plot[Ti[z], {z, 0, ts + tsl}, GridLines -> {{ts}, {0}}]

eqn2 = ρ*c*D[T[t, z], t] - k*D[T[t, z], z, z];
Subscript[Γ, D] = DirichletCondition[T[t, z] == Tbl, z == ts + tsl];
Subscript[Γ, N] = NeumannValue[-28, z == 0];

soln =
NDSolveValue[
{eqn2 == Subscript[Γ, N], Subscript[Γ, D], T[0, z] == Ti[z]},
T, {t, 0, 1}, {z, 0, ts + tsl},
Method ->
{"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement"}}}]

Animate[
Plot[soln[t, z], {z, 0, ts + tsl},
GridLines -> {{ts, ts + .00005}, {0}},
PlotRange -> {99, 111}],
{t, 0, 1}]
Plot[soln[t, ts + .00005], {t, 0, 1},
GridLines -> {{0}, {107}},
PlotRange -> Full]
Plot[soln[0.01, z], {z, 0, ts + tsl}, GridLines -> {{ts}, {0}}]


• You can experiment with "SpatialDiscretization" -> {"FiniteElement", "MeshOptions" -> {"MeshOrder" -> 1}} or things like "SpatialDiscretization" -> {"FiniteElement", "MeshOptions" -> {"MaxCellMeasure" -> (ts + tsl)/1011}} Commented Nov 26, 2015 at 1:26
• @ThiesHeidecke, OP uses the FEM and equivalent of MinPoints is changing the MaxCellMeasure and the difference order reduction can be done with specifying a "MeshOrder". Commented Nov 26, 2015 at 1:28