# Heat balance in layer coupled via BC

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.

• 1. 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 to NeumannValue correctly. – xzczd Mar 10 at 2:45
• T0 and Lf are no longer defined. How you have sketched the problem, would it not lead to a trivial solution of all temperatures equaling T0 for all times? – Tim Laska Mar 10 at 17:50

As I commented in your previous question 240844, currently, Mathematica does not support internal NeumannValues. It does, however, allow you to create region-dependent material properties.

When you have multiple regions, it is often conducive to construct a mesh manually with each material region marked. Here is an example workflow.

# User supplied parameters

(*Import required FEM package*)
Needs["NDSolveFEM"];
(*OP defined 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];
Lf = 8/1000;
Tair = 120;
Length1 = 0.1;
time1 = 60*Length1/4;
(*T0 was not supplied*)
T0 = 0;
(*Create property arrays*)
λ = {kf, km, ka};
ρ = {rhof, rhom, rhoa};
cp = {Cpf, Cpm, Cpa};


# Construct ElementMesh

The mesh is constructed following my previous answer to a multiple material problem here 223034.

(*Create mesh with region markers*)
g = {Lf, thickness, Lf} // N;(*thicknesse*)
gw = {0}~Join~Accumulate[g] + (L1 - Lf);
bmesh = ToBoundaryMesh["Coordinates" -> Partition[gw, 1],
"BoundaryElements" -> {PointElement[{{1}, {2}, {3}, {4}}]}]; nrEle \
= 10; pt = Partition[gw, 2, 1];
regmarkers =
Transpose[{Partition[(Mean /@ pt), 1], {1, 2, 3},
Abs[Subtract @@@ pt]/nrEle}];
mesh = ToElementMesh[bmesh, "RegionMarker" -> regmarkers];


# Create region dependent physical properties

(*Create region dependent physical properties*)
rhocp = Evaluate[Piecewise[{{ρ[[1]] cp[[1]], ElementMarker == 1},
{ρ[[2]] cp[[2]], ElementMarker == 2},
{ρ[[3]] cp[[3]], ElementMarker == 3}}]];
k = Evaluate[Piecewise[{{λ[[1]], ElementMarker == 1},
{λ[[2]], ElementMarker == 2},
{λ[[3]], ElementMarker == 3}}]];


# Set up PDE system

We will use the HeatTransferPDEComponent and HeatTransferValue to set up the heat transfer PDE.

(*Set up PDE*)
vars = {Θ[t, x], t, {x}};
pars = <|"MassDensity" -> 1, "SpecificHeatCapacity" -> rhocp,
"ThermalConductivity" -> {k}|>;
pars["BC1"] = <|"HeatTransferCoefficient" -> h,
"AmbientTemperature" -> Tair|>;
Γleft =
HeatTransferValue[x == gw[[1]], vars, pars, "BC1"];
Γright =
HeatTransferValue[x == gw[[-1]], vars, pars, "BC1"];
ics = Θ[0, x] == T0;
eqn = HeatTransferPDEComponent[vars,
pars] == Γleft + Γright;


For those who do not have version 12.2, here is the input form of eqn:

Inactive[Div][{{-k}} . Inactive[Grad][Θ[t, x], {x}], {x}] +
rhocp*Derivative[1, 0][Θ][t, x] ==
NeumannValue[h*(Tair - Θ[t, x]), x == gw[[-1]]] +
NeumannValue[h*(Tair - Θ[t, x]), x == gw[[1]]]


# Solve and plot the PDE system

Tfun = NDSolveValue[{eqn, ics}, Θ, {t, 0,
500000}, {x} ∈ mesh];
Plot3D[Tfun[t, x], {t, 0, 100000}, {x, First@gw, Last@gw},
ColorFunction -> "ThermometerColors", PlotRange -> All,
MeshFunctions -> {#3 &}, Mesh -> 24]


• Well, I'm afraid you've posted an answer too fast. The NeumannValue in OP's code seems to suggest OP is trying to set b.c. like that in this question. (This should be clarified by OP, of course. ) – xzczd Mar 10 at 3:38
• @xzczd Perhaps I should have waited for additional clarification from the OP. You would agree, though, that with the standard out-of-the-box FEM, you cannot have an internal NeumannValue. Agree? – Tim Laska Mar 10 at 5:58
• Yeah, at least up to v12.2. – xzczd Mar 10 at 6:17
• indeed, I am trying to couple the heat balances in each layer via the BCs. Or in other words: I would like to impose either the value of the flux over the interface (layer1 and layer3) or that the flux is constant over the interface of the layers (layer 2). Also relevant to mention that I am using Mathematica 12.0.0, so I do not have access to the latest functions for heat transfer. – Luigi Mar 10 at 8:17
• @Luigi The dimensions of the Neumann value should be that of heat flux. Therefore, dividing by the thermal conductivity would not make it dimensionally correct. Given that you have a heat transfer coefficient, I presume you are trying to model fluid/solid heat transfer. This implies a flow perpendicular to the one dimension your modeling and the heat transfer coefficient should depend on both the Reynolds number and the Prandtl number. To get a better representation, you should consider a two-dimensional model. – Tim Laska Mar 16 at 2:43

This is an extended comment. You can have internal NeumannValues - but you will have to generate a mesh that actually has an internal boundary:

Needs["NDSolveFEM"]
bmesh = ToBoundaryMesh[
"Coordinates" -> {{0, 0}, {1, 0}, {2, 0}, {2, 1}, {1, 1}, {0, 1}},
"BoundaryElements" -> {LineElement[{{1, 2}, {2, 3}, {3, 4}, {4,
5}, {5, 6}, {6, 1}}], LineElement[{{2, 5}}]}];
mesh = ToElementMesh[bmesh];
mesh["Wireframe"]


if = NDSolveValue[{Laplacian[u[x, y], {x, y}] ==
1 + NeumannValue[u[x, y], {x == 1}]},
u, {x, y} \[Element] mesh];
Plot3D[if[x, y], {x, y} \[Element] mesh]


Update:

1D example:

bmesh = ToBoundaryMesh["Coordinates" -> {{0}, {1}, {2}},
"BoundaryElements" -> {PointElement[{{1}, {2}, {3}}]}];
mesh = ToElementMesh[bmesh];
eqn = {Laplacian[u[x], {x}] == 1 + NeumannValue[10, {x == 1}],
DirichletCondition[u[x] == 0, x == 0 || x == 2]};
if = NDSolveValue[eqn, u, {x} \[Element] mesh];
Plot[if[x], {x} \[Element] mesh]


• Interesting. How is the outward-facing unit normal $\overset{\rightharpoonup }{n}$ determined in this case? – xzczd Mar 10 at 12:45
• @xzczd, the normal vector is never determined for NeumannValue - regardless it is on an exterior boundary or an interior boundary. And that's the reason I do not quite understand when people make use of internal boundaries. The normal vector is only computed if you explicitly request it by using BoundaryUnitNormal[x,..] from the FEM context. – user21 Mar 10 at 13:05