4

I am solving the heat equation (diffusion only) over a 3-layer system (physical properties varying from layer to layer).

km = 0.128;
rhom = 925;
Cpm = 1550;

ka = 0.024;
rhoa = 1.292;
Cpa = 1003;


thickness = 2*0.136;

L1 = N[-thickness/2];
L2 = N[thickness/2];
T01 = 120;
T02 = 20;
T03 = 120;


Length1 = 0.5;


v= 40/60;


time1 = Length1/v;

For describing the properties change in space, I use the following equation:

slope = 1000;
SmoothedStepFunction[fL_, fmax_, fR_, tsL_, tsR_, m_] := 
 Function[t, (fL*Exp[tsL*m] + fmax*Exp[m*t])/(Exp[tsL*m] + 
     Exp[m*t]) - (fR*Exp[tsR*m] + fmax*Exp[m*t])/(Exp[tsR*m] + 
     Exp[m*t]) + fR];

rhoCp[x_] := 
  SmoothedStepFunction[rhoaCpa, rhomCpm, rhoa*Cpa, L1, L2, slope][x];
k[x_] := SmoothedStepFunction[10^6ka, 10^6km, 10^6*ka, L1, L2, 
    slope][x];

The actual heat balance is solved here:

heateq = rhoCp[x]*D[u[x, t], t] == 
   Inactive[Div][{{k[x]}}.Inactive[Grad][u[x, t], {x}], {x}];

ic[x_] := 
  Piecewise[{{T01, x < L1}, {T02, L1 <= x <= L2}, {T03, x > L2}}];


sol1 = First[
   NDSolve[{heateq, u[x, 0] == ic[x]}, 
    u, {x, -2thickness, 2thickness}, {t, 0, time1}, 
    Method -> {"MethodOfLines", 
      "SpatialDiscretization" -> {"FiniteElement", 
        "MeshOptions" -> {"MaxCellMeasure" -> 0.01}}}]];

Despite it works well, I still have a question: next to the heat flux continuity over the boundaries of materials (Inactive[Div][{{k[x]}}.Inactive[Grad][u[x, t]), are there other boundary conditions applied implicitly with this approach?

user21
  • 39,710
  • 8
  • 110
  • 167
Luigi
  • 1,301
  • 8
  • 14
  • 1
    Currently, Mathematica does not support Neumann values on interfaces. It will, however, maintain heat flux continuity across a material interface. – Tim Laska Mar 01 '21 at 13:29
  • in the formulation above, when does Inactive[] become active? – Luigi Mar 01 '21 at 16:32
  • Form the documentation Numerical Solution of Partial Differential Equations, "Inactive forms can be activated using Activate, but commands like NDSolve can also work with the inactive forms directly." It is usually best to express the FEM operator in Coefficient Form. Because Div/Grad evaluate immediately, the coefficient form is lost without Inactive, but NDSolve can process it. – Tim Laska Mar 01 '21 at 17:04
  • by using Inactive[] I can leave out two boundaries conditions. why? – Luigi Mar 01 '21 at 20:29
  • 3
    The default BC for a FEM problem is NeumannValue[0, (x == xmin || x == xmax)], which is an insulated wall. If you supply no BCs, the default will be assumed. Inactive does not change the default value. – Tim Laska Mar 01 '21 at 20:40

0 Answers0