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, 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, -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 Neumann boundary conditions (`Inactive[Div][{{k[x]}}.Inactive[Grad][u[x, t]`

) that are applied over the interfaces to connect the solutions in each domain, are there other boundary conditions applied implicitly?