I am trying to solve the heat balance in a 3 layer system in Mathematica 12.0.0. Each layer has different thermal properties. The layers are coupled via the boundary conditions (Neumann type) to ensure that heat flux is preserved over the interfaces with the central layer (layer 2). The system is schematically depicted here:
I came up with the following approach:
(*definition of parameters*)
km = 0.128; rhom = 975; Cpm = 1948;
ka = 0.024; rhoa = 1.292; Cpa = 1003;
kf = 0.33; rhof = 940; Cpf = 2100;
h = 5;
thickness = 0.18;
L1 = N[-thickness/2];
L2 = N[thickness/2];
Tair = 120;
timemax=2;
(*heat equation*)
heateqf1 =
rhof*Cpf*D[uf1[x, t], {t, 1}] - (10^6*kf)*D[uf1[x, t], {x, 2}] ==
NeumannValue[h*(u[x,t]- uf1[x, t]), x == L1]; (*miss the -infinity bc*)
heateqm =
rhom*Cpm*D[u[x, t], {t, 1}] - (10^6*km)*D[u[x, t], {x, 2}] ==
NeumannValue[h*(uf1[x, t] - u[x, t]), x == L1] +
NeumannValue[h*(u[x, t] - uf2[x, t]), x == L2];
heateqf2 =
rhof*Cpf*D[uf2[x, t], {t, 1}] - (10^6*kf)*D[uf2[x, t], {x, 2}] ==
NeumannValue[h*(u[x, t] - uf2[x, t]), x == L2] ; (*miss the +infinity bc*)
icf1 = uf1[x, 0] == T0;
icm = u[x, 0] == T0;
icf2 = uf2[x, 0] == T0;
bc1 = uf1[-10*thickness, t] == Tair;
bc2 = uf2[10*thickness, t] == Tair;
(*solution finding*)
sol1 = NDSolveValue[{{heateqf1, heateqm, heateqf2}, {icf1, icm, icf2,
bc1, bc2}}, {uf1, u, uf2}, {t, 0, timemax}, {x, L1 - Lf,
L2 + Lf},
Method -> {"MethodOfLines", "TemporalVariable" -> t,
"SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> {"Length" -> 0.001}}}}];
Besides, I am not sure that bc1 and bc2 are good definitions for the boundary conditions, I do not make it to the result.
NDSolve
cannot directly handle such kind of problem i.e. PDEs in different parts of the domain are different, but in many cases it's not too hard to circumvent. 2. Can you show us the b.c.s in traditional math notation? I suspect you haven't interpret them toNeumannValue
correctly. $\endgroup$