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[rhoa*Cpa, rhom*Cpm, rhoa*Cpa, L1, L2, slope][x];
k[x_] := SmoothedStepFunction[10^6*ka, 10^6*km, 10^6*ka, L1, L2, 

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, -2*thickness, 2*thickness}, {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?

  • 1
    $\begingroup$ Currently, Mathematica does not support Neumann values on interfaces. It will, however, maintain heat flux continuity across a material interface. $\endgroup$ – Tim Laska Mar 1 at 13:29
  • $\begingroup$ in the formulation above, when does Inactive[] become active? $\endgroup$ – Luigi Mar 1 at 16:32
  • $\begingroup$ 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. $\endgroup$ – Tim Laska Mar 1 at 17:04
  • $\begingroup$ by using Inactive[] I can leave out two boundaries conditions. why? $\endgroup$ – Luigi Mar 1 at 20:29
  • 3
    $\begingroup$ 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. $\endgroup$ – Tim Laska Mar 1 at 20:40

Your Answer

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

Browse other questions tagged or ask your own question.