Heat transfer with functions defined on different domains

Question

I'm attempting to model a situation in which a polymer initially at a higher temperature is sandwiched between two cooler metallic mold pieces with conductive boundary conditions in between the polymer and mold. On either outer side of the mold, there is convection with another fluid. Since the boundary conditions at either side of the sandwiched polymer where contact is made with the mold are the same, I'm focusing on only one of the mold pieces, with the mold spanning from $0 \leq z \leq d$ and the polymer spanning from $-w \leq z \leq 0$. The following equations model the setup:

$$\frac{\partial T_{mold}(z,t)}{\partial t}=a_{mold}\frac{\partial^2T_{mold}(z,t)}{\partial z^2}, 0 \leq z \leq d$$ $$\frac{\partial T_{poly}(z,t)}{\partial t}=a_{poly}\frac{\partial^2T_{poly}(z,t)}{\partial z^2}, -w \leq z \leq 0$$

Boundary conditions: $$k_{mold}\frac{\partial T_{mold}(z,t)}{\partial z}|_{z=d}=h_c(T_{mold}(z,t)-T_c)$$ $$-k_{mold}\frac{\partial T_{mold}(z,t)}{\partial z}|_{z=0}=-k_{poly}\frac{\partial T_{poly}(z,t)}{\partial z}|_{z=0}$$ $$-k_{mold}\frac{\partial T_{mold}(z,t)}{\partial z}|_{z=-w}=-k_{poly}\frac{\partial T_{poly}(z,t)}{\partial z}|_{z=-w}$$

Initial conditions: $$T_{mold}(z,t=0)=T_m$$ $$T_{poly}(z,t=0)=T_p$$

I tried to approach solving this based on the wonderful iterative approach of @Alex Trounev from this post using the following code (the signs on NeumannValue might be incorrect due to my inexperience):

amold = 0.004;
Tm = 323;
apoly = 0.0001;
Tp = 433;
Tc = 303;
h = 5000;
d = 0.005;
w = 0.005;

poly[0][z_, t_] := Tp;
Do[mold = 
  NDSolveValue[{D[Tmold[z, t], t] - 
      amold*D[Tmold[z, t], {z, 2}] == 
     NeumannValue[-h (Tmold[z, t] - Tc)/
         kmold, z == d] + 
      NeumannValue[kpoly/kmold)*(D[poly[i - 1][z, t], z] /. z -> 0), z == 0], 
    Tmold[z, 0] == Tm}, 
   Tmold, {z, 0, d}, {t, 0, 600}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"FiniteElement", 
       "MeshOptions" -> {"MaxCellMeasure" -> {"Length" -> 0.001}}}}]; 
 poly[i] = 
  NDSolveValue[{D[Tpoly[z, t], t] - 
      apoly*D[Tpoly[z, t], {z, 2}] == 
     NeumannValue[kmold/kpoly)*(D[mold[z, t], z] /. z -> 0), z == 0] + 
      NeumannValue[kmold/kpoly)*(D[mold[z, t], z] /. z -> 0), z == -w], 
    Tpoly[z, 0] == Tp}, 
   Tpoly, {z, -w, 0}, {t, 0, 600}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"FiniteElement", 
       "MeshOptions" -> {"MaxCellMeasure" -> {"Length" -> 
            0.001}}}}];, {i, 1, 20}]

However, I get a very unstable solution, and the kernel crashes after ~30 iterations. I also tried solving the system of PDEs without an iterative approach below, but I get the error that the boundary condition at $z=-w$ is not specified on a single edge of the boundary of the computational domain:

sys1 = {D[Tmold[z, t], t] - 
     amold*D[Tmold[z, t], {z, 2}] == 
    NeumannValue[-h (Tmold[z, t] - Tc)/
       kmold, z == d], 
   Tmold[z, 0] == Tm, 
   D[Tpoly[z, t], t] == 
    apoly*
     D[Tpoly[z, t], {z, 2}], -kmold*(D[Tmold[z, t], z] /. z -> 0) == -kpoly*(D[Tpoly[z, t], z] /. z -> 0), -kmold*(D[Tmold[z, t], z] /. z -> 0) == -kpoly*(D[Tpoly[z, t], z] /. z -> -w), 
   Tpoly[z, 0] == Tp};

heat = NDSolveValue[sys1, {Tmold,Tpoly}, {z, -w, d}, {t, 0, 600}, 
  Method -> {"MethodOfLines", 
    "SpatialDiscretization" -> {"FiniteElement", 
      "MeshOptions" -> {"MaxCellMeasure" -> {"Length" -> 0.001}}}}]

I was wondering if anyone had insight into how I can improve either of these two approaches. Thank you again in advance for any help!

Answer 1

As I understood the question, the polymer is cooled on two sides by the mold. This implies a symmetry condition at the center of the polymer as shown in the image below.

Create mesh

In multi-material domains, it is sometimes useful to refer to the different domains by ElementMarker. The following workflow will create a two region mesh. The polymer is indicated by $\color{Red}{Red}$ and the mold by $\color{Blue}{Blue}$.

Needs["NDSolve`FEM`"]
(* User Supplied Parameters *)
amold = 0.004;
Tm = 323;
apoly = 0.0001;
Tp = 433;
Tc = 303;
h = 5000;
d = 0.005;
w = 0.005;
(*Put user supplied parameters into list*)
g = {w/2, d};(*thicknesses*)
gw = {0}~Join~Accumulate[g]; (* Nodes *)
λ = {0.4, 44.5};
ρ = {975, 7850};
cp = {1550, 475};
(*functions to build multi-domain mesh with markers*)
ClearAll[seg, appendCrdRight]
seg[thick_, nelm_, marker_] := Module[{crd, inc, marks},
  crd = Subdivide[0, thick, nelm];
  inc = Partition[Range[crd // Length], 2, 1];
  marks = ConstantArray[marker, inc // Length];
  <|"c" -> crd, "i" -> inc, "m" -> marks|>
  ]
appendCrdRight[a1_, a2_] := Module[{crd, inc, marks, len, lcrd},
  len = a1["c"] // Length;
  lcrd = a1["c"] // Last;
  inc = Join[a1["i"], a2["i"] + len - 1];
  crd = Join[a1["c"], Rest[a2["c"] + lcrd]];
  marks = Join[a1["m"], a2["m"]];
  <|"c" -> crd, "i" -> inc, "m" -> marks|>]
(*Build two region mesh with markers*)
a = Fold[appendCrdRight, MapIndexed[seg[#1, 200, First[#2]] &, g]];
mesh = ToElementMesh["Coordinates" -> Partition[a["c"], 1], 
   "MeshElements" -> {LineElement[a["i"], a["m"]]}, 
   "BoundaryElements" -> {PointElement[{{1}, {a["c"] // Length}}, {1, 
       2}]}];
Show[mesh["Wireframe"["MeshElementStyle" -> {Red, Blue}]], 
 PlotRange -> {-0.0001, 0.0001}]

Set up and solve PDE system

To use the HeatTransferPDEComponent as described in the Heat Transfer Tutorial, you will need the mass density, specific heat capacity, and thermal conductivity versus just using the diffusivities is alone. Above, I used typical material parameters for polyethylene and stainless steel for the polymer and mold, respectively. An example workflow is shown below:

(*Create region dependent physical properties*)
dens = Evaluate[Piecewise[{{ρ[[1]], ElementMarker == 1},
     {ρ[[2]], True}}]];
specheat = Evaluate[Piecewise[{{cp[[1]], ElementMarker == 1},
     {cp[[2]], True}}]];
thermcon = Evaluate[Piecewise[{{λ[[1]], ElementMarker == 1},
     {λ[[2]], True}}]];
(*Set up variables and parameters for PDE*)
vars = {T[t, x], t, {x}};
pars = <|"MassDensity" -> dens, "SpecificHeatCapacity" -> specheat, 
   "ThermalConductivity" -> {thermcon}, 
   "HeatTransferCoefficient" -> h, "AmbientTemperature" -> Tc|>;
(*Convective heat transfer boundary condition*)
Γc = HeatTransferValue[x == Last@gw, vars, pars];

ic1 = T[0, x] == If[x <= gw[[2]], Tp, Tm];
eqn = HeatTransferPDEComponent[vars, pars] ==
   Γc;
Tfun = NDSolveValue[{eqn, ic1}, T, {t, 0, 100}, {x} ∈ mesh];
frames = Plot[Tfun[#, x], {x} ∈ mesh, 
     PlotRange -> {Tc, 1.025 Tp}, 
     ColorFunction -> Function[{x, y}, If[x <= gw[[2]], Red, Blue]], 
     ColorFunctionScaling -> False] & /@ Subdivide[0, 50, 60];
ListAnimate@frames

Pre-Mathematica 12.2 versions

If you do not have at least Mathematica version 12.2, then you will not have access to the HeatTransferPDEComponent. In that case, you will need to replace the code above from pars on down with:

pars120 = {md -> dens, shc -> specheat, tc -> thermcon, htc -> h, 
   Ta -> Tc};
Γc = 
  NeumannValue[htc (Ta - T[t, x]), x == Last[gw]] /. pars120;
ic1 = T[0, x] == If[x <= gw[[2]], Tp, Tm];
eqn = Inactive[Div][{{-tc}} . Inactive[Grad][T[t, x], {x}], {x}] + 
    md*shc*D[T[t, x], {t}] /. pars120;
pde = eqn == Γc;
Tfun = NDSolveValue[{pde, ic1}, T, {t, 0, 100}, {x} ∈ mesh];
frames = Plot[Tfun[#, x], {x} ∈ mesh, 
     PlotRange -> {Tc, 1.025 Tp}, 
     ColorFunction -> Function[{x, y}, If[x <= gw[[2]], Red, Blue]], 
     ColorFunctionScaling -> False] & /@ Subdivide[0, 50, 60];
ListAnimate@frames

Verification of temperature and heat flux continuity across material interfaces

The following verifies that the above approach preserves both temperature and heat flux continuity across material interfaces without the need of prescribing NeumannValue.

(*Create region dependent thermal conductivity function*)
kfn = Function[{x}, 
   Piecewise[{{(λ[[1]] + λ[[2]])/2, 
      Abs[x - gw[[2]]] < 1.*^-8}, {λ[[1]], 
      x < gw[[2]]}, {λ[[2]], True}}]];
(*Set up derived temperature gradient and flux functions*)
Clear[gradT, flux]
gradT[t_, x_] := Derivative[0, 1][Tfun][t, x];
flux[t_, x_] := -kfn[x] gradT[t, x];
(*Create temperature and heat flux plotting functions*)
Clear[plotter]
plotter[pf_, time_, ymin_, ymax_, lbl_] := 
 With[{pft = pf /. {t -> time}, lblt = lbl /. {t -> time}, 
   plotOptions = Options[Plot]},
  (* WithCleanup requires version 12.2 or higher *)
  WithCleanup[
   SetOptions[
    Plot, {PlotLabel -> lblt, 
     PlotRange -> {{gw[[1]], gw[[-1]]}, {0.95 ymin, 1.05 ymax}}, 
     PlotStyle -> {Directive[Thickness[0.01]]}, 
     Epilog -> {Opacity[0.4], HatchFilling[Pi/4, .05, 5], 
       EdgeForm[{Black, Thick}], 
       Rectangle[{gw[[2]], ymin}, {gw[[3]], ymax}], 
       Inset[Style["Mold", Large, Bold, Background -> White], 
        Mean[{{gw[[2]], ymin}, {gw[[3]], ymax}}]]}, 
     ColorFunction -> "ThermometerColors"}],
   Plot[pft, {x, gw[[1]], gw[[-1]]}],
   SetOptions[Plot, plotOptions]
   ]
  ]
tpltfn = plotter[Tfun[t, x], #, Tc, Tp, 
    StringForm["Temperature @ t=``", NumberForm[N [t], {3, 2}]]] &;
fpltfn = plotter[flux[t, x], #, 0, 250000, "Heat Flux"] &;
imgs = ParallelTable[
   GraphicsColumn[{tpltfn[t], fpltfn[t]}], {t, 0, 3, 1/33}];
ListAnimate@imgs