Solution to PDE with "MethodOfLines" is not returning any values

Question

I have been working on this code for quite a while now trying different types of solutions and I cannot get any actual values out of the NDSolveValue function. I have included my code so you can see what I am trying to do and also see how there are no values in the solutions.

<< "NDSolve`FEM`"

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

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


eqn1 = k*\!\(
\*SubscriptBox[\(\[PartialD]\), \(z\)]\(T1[z]\)\) + 28;
Tbl = 100;
Subscript[\[CapitalGamma]1, D] = 
  DirichletCondition[T1[z] == Tbl, z == ts + tsl];

BCr = NDSolveValue[{eqn1 == 0, Subscript[\[CapitalGamma]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 = \[Rho]*c*\!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(T[t, z]\)\) - k*\!\(
\*SubscriptBox[\(\[PartialD]\), \(z, z\)]\(T[t, z]\)\);
Subscript[\[CapitalGamma], D] = 
  DirichletCondition[T[t, z] == Tbl, z == ts + tsl];

soln[t_, z_] = 
 NDSolveValue[{eqn2 == 0, Subscript[\[CapitalGamma], 
   D], (D[T[t, z], z] /. z -> 0) == -28/ks, T[0, z] == Ti[z]}, 
  T, {t, 0, .1}, {z, 0, ts + tsl}, 
  Method -> {"PDEDiscretization" -> {"MethodOfLines", 
      "SpatialDiscretization" -> {"TensorProductGrid", 
        MinStepSize -> .00000005, MaxStepSize -> .00000005}}}]

Animate[Plot[soln[t, z], {z, 0, ts + tsl}, 
  GridLines -> {{ts}, {0}}], {t, 0, .1}]

Can someone explain how I can get an actual solution out of this. I believe there is an issue that is occurring due to the discontinuity at position .00025 and have tried to use the "DiscontinuityProcessing" options to no avail. Any suggestions on how to get the solution are extremely appreciated.

Answer 1

How about this:

tEnd = 10.;
solnFEM = 
 NDSolveValue[{eqn2 == NeumannValue[-28, z == 0], 
   Subscript[\[CapitalGamma], D], T[0, z] == Ti[z]}, 
    T, {t, 0, tEnd}, {z, 0, ts + tsl}, 
    Method -> {"PDEDiscretization" -> {"MethodOfLines", 
            "SpatialDiscretization" -> {"FiniteElement"}}}]

Seems to get reasonably close:

Plot[solnFEM[tEnd, z], {z, 0, ts + tsl}]

solnFEM[tEnd, 0] - BCr[0]
0.0602549455763608`

Answer 2

I have managed to get a solution and it appears somewhat correct (not entirely physically accurate), I don't know why this works but by not defining the independent variables t and z in for the solution achieves results. I switched the following line for the new line of code given below.

Original Code:

soln[t_, z_]  = NDSolveValue[{eqn2 == 0, Subscript[\[CapitalGamma], 
   D], (D[T[t, z], z] /. z -> 0) == -28/ks, T[0, z] == Ti[z]}, 
  T, {t, 0, 100}, {z, 0, ts + tsl}, 
  Method -> {"PDEDiscretization" -> {"MethodOfLines", 
      "SpatialDiscretization" -> {"TensorProductGrid", 
        MinStepSize -> .00000005, MaxStepSize -> .00000005}}}]

Working Code:

soln = NDSolveValue[{eqn2 == 0, Subscript[\[CapitalGamma], 
   D], (D[T[t, z], z] /. z -> 0) == -28/ks, T[0, z] == Ti[z]}, 
  T, {t, 0, 100}, {z, 0, ts + tsl}, 
  Method -> {"PDEDiscretization" -> {"MethodOfLines", 
      "SpatialDiscretization" -> {"TensorProductGrid", 
        MinStepSize -> .00000005, MaxStepSize -> .00000005}}}]

The issue is that the solution as time goes to infinity does not approach the steady state solution that is found initially. Given that the steady state solution and the transient solution as time gets very large should be identical, I believe there is some issue with the variable k, \[Rho], andcand them not functioning properly for the second differential solution.