5
$\begingroup$

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!

Improve this question
$\endgroup$
2
  • 3
    $\begingroup$ Energy flux continuity should be preserved through material interfaces, so NeumannValues are not required at the interface. You will want to express the heat equation in coefficient form, which will require that you break the thermal diffusivity into density, heat capacity, and thermal conductivity parameters. A fairly extensive heat transfer tutorial can be found here. $\endgroup$ – Tim Laska May 17 at 3:08
  • $\begingroup$ Thank you @TimLaska! I'll try taking a look at that $\endgroup$ – Rpj May 17 at 14:05
7
$\begingroup$

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.

System description

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

Mesh domains

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

Simulation results

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

Temperature and flux continuity

Improve this answer
$\endgroup$
3
  • $\begingroup$ Thank you so much for your help and detailed explanation @Tim Laska!! When I try putting your code into Mathematica, I get the error NDSolveValue::derivs: No derivatives of dependent variables were found in the equations. NDSolveValue is designed to solve differential or differential algebraic equations. Use NSolve or FindRoot to numerically solve algebraic equations. Do you know off the top of your head what might be causing this? $\endgroup$ – Rpj May 17 at 23:49
  • 2
    $\begingroup$ @Rpj Evaluate $Version in Mathematica. I am on 12.2.0 for Microsoft Windows (64-bit) (December 12, 2020). This approach will only work on V12.2 or later. $\endgroup$ – Tim Laska May 18 at 0:20
  • 1
    $\begingroup$ That fixed the issue. Thanks again for your help, @Tim Laska!! $\endgroup$ – Rpj May 18 at 3:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.