Heat balance in layer coupled via BC

Question

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 = 
  rhofCpfD[uf1[x, t], {t, 1}] - (10^6kf)D[uf1[x, t], {x, 2}] == 
   NeumannValue[h(u[x,t]- uf1[x, t]), x == L1]; (miss the -infinity bc*)
heateqm = 
  rhomCpmD[u[x, t], {t, 1}] - (10^6km)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 = 
  rhofCpfD[uf2[x, t], {t, 1}] - (10^6kf)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.

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 '21 at 17:50

Tim Laska · Accepted Answer · 2021-03-10T15:28:49.037

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["NDSolve`FEM`"];
(*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 '21 at 03: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 '21 at 05:58
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 '21 at 08:17
@Luigi As requested in the comment under the question, can you show us the required b.c. in traditional math notation (using $\LaTeX$, if you're having difficulty with $\LaTeX$ then screenshot or even photo is OK) to avoid possible ambiguity? — xzczd, Mar 10 '21 at 09:32
@xzczd, please see my comment on internal NeumannValues. You can have internal NeumannValues. — user21, Mar 10 '21 at 12:36
@Luigi, I added the eqn output as you indicated you do not have version 12.2. Also, in a transient problem, the fluxes will not be constant over the interface. However, you will see in this formulation, the fluxes across the interface tend towards a constant steady-state value. — Tim Laska, Mar 10 '21 at 15:35
@TimLaska: the boundary conditions NeumannValue[h*(Tair - Θ[t, x]), x == gw[[-1]]] should actually be NeumannValue[h*(Tair - Θ[t, x])/k, x == gw[[-1]]]? Otherwise, we are missing the conductivity in the heat flux. Agree? — Luigi, Mar 15 '21 at 15:23
@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 '21 at 02:43

user21 · Answer 2 · 2021-03-10T14:25:39.973

1

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["NDSolve`FEM`"]
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]

I still think Tim Laska's answer is what you want.

edited Mar 10 '21 at 14:25

answered Mar 10 '21 at 12:34

user21

39,710
8
110
167

Interesting. How is the outward-facing unit normal $\overset{\rightharpoonup }{n}$ determined in this case? – xzczd Mar 10 '21 at 12:45
my problem is 1D. I should first become familiar on how to make appropriate meshes. – Luigi Mar 10 '21 at 12:56
@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 '21 at 13:05
@Luigi, In 1D this would be adding points where you want this condition. – user21 Mar 10 '21 at 13:05
@user21 would it be possible to provide me with an example on how to connect it to my problem? – Luigi Mar 10 '21 at 14:17
@Luigi, see update. – user21 Mar 10 '21 at 14:26
@user21, for some reasons, it does not work on my machine (v12.0.0) – Luigi Mar 10 '21 at 15:02
I was not aware of the BoundaryUnitNormal option, although I see it in the documentation. However, I did not see any example usage beyond "Automatic." It would be interesting to see how it is used for curved lines and surfaces. It would be a useful feature to model parallel heat transfer mechanisms such as radiation or jump conditions in mass transfer. There is something strange with your example in that you have two oppositely directed fluxes converging at a point in a steady problem. There should be some accumulation. – Tim Laska Mar 10 '21 at 15:54
@TimLaska, here is an example of a non Automatic usage. The computation of the BoundaryUnitNormal is expensive, so when a BC operator returns it and the boundary normal is known, then things can be sped up a bit. – user21 Mar 11 '21 at 06:13
@TimLaska, if you have an examples with a BoundaryUnitNormal i'd be keen to take a look. – user21 Mar 11 '21 at 06:58
@TimLaska, how do you think the solution for an internal Neumann value should look like? – user21 Mar 11 '21 at 06:59
@user21 Thanks for the link! I will try some experiments with it and try to find you some examples. Regarding what the NeumannValue should look like for a steady problem without radiation, the temperature profile should be monotonic like this. It would be good to find some verification examples with internal NeumannValues to tune the setup of these types of problems. – Tim Laska Mar 11 '21 at 15:37