I'm having some problems with NDSolve and the problem of conduction of heat. Specifically I'm having problems with the continuity of heat flux. First, let me explain what is the physical problem:

enter image description here

Let´s say I want to study the heat conduction in a composite material with n layers (see figure above), each one presenting a different value for its thermal properties (thermal conductivies, $k_1, k_2,..., k_n$, specific heat, $c_{p_1},...,c_{p_n} $ and densities, $\rho_1, ..., \rho_n$). In this case we could write:

$ k[z] = k_1 + \sum_\limits{i=1}^{n-1} (k_{i+1}-k_i) H[z-z_i] $

$ \rho[z] = \rho_1 + \sum_\limits{i=1}^{n-1} (\rho_{i+1}-\rho_i) H[z-z_i] $

$ c_p[z] = c_{p_1} + \sum_\limits{i=1}^{n-1} (c_{p_{i+1}}-c_{p_{i}}) H[z-z_i] $

where $z$ indicates the depth and $H$ is a function that changes abrubtly (ideally is the UnitStep function). As you can see, there will be discontinuities in $k[z]$.

The PDE describing our problem is the heat equation:

$ \rho c_p \frac{\partial T}{\partial t} = \nabla.(k\nabla T) + Q_{ext} $

where $ \rho, c_p, k$ and $ Q_{ext} $ are the density, the specific heat, the thermal conductivity and heat generation rate of the external heat source.

In our problem, we have convection of heat at the top of layer 1 and a constant temperature at the bottom of the last layer. We also require the continuity of both the temperature and the heat flux. Therefore, the boundary conditions can be expressed as:

$ k_1 \frac{\partial T}{\partial z}_{z=0} = h (T_{env}-T_{z=0}) $

$ T_{z=z_n} = T_{wall} $

$ k_n \frac{\partial T}{\partial z}_{z\rightarrow z_n^+} = k_{n+1} \frac{\partial T}{\partial z}_{z\rightarrow z_n^-} \text{; for n = 1,...,n-1} $

where $h$, $T_{env}$ and $T_{wall}$ are just constants.

The last written set of boundary conditions indicates the continuity of the heat flux, $q'' = - k \frac{\partial T}{\partial z}$, between the interfaces of the different layers. They are the ones I don't know how to impose in NDSolveValue.

This is what I have been doing so far:

(*Required packages*)

(*Some paremeters*)
conds = {αp -> 0.8*10^-4, Tenv -> 22, Tcore -> 35, ha -> 12, 
   Ptotal -> 1.30, Vtotal -> 2.25*10^-6, rlaser -> 0.13, 
   Llaser -> 4.5*10^-3, lasert -> 3*60};
dvalues = {1*10^-3, 3*10^-3, 3*10^-3, 1*10^-3, 
  6*10^-3}; (*thicknesses*)
kcond = {0.235, 0.445, 0.445, 0.185, 
  0.51}; (*thermal conductivity of layers*)
ρdens = {1200, 1200, 1200, 1000, 1085}; (*densities*)
cspe = {3589, 3300, 3300, 2674, 3800}; (*specific heat*)
(*Depth of layers*)
zvalues = Table[Sum[dvalues[[j]], {j, 1, i}], {i, 1, Length[dvalues]}];
noflayers = Length[zvalues]; (*number of layers*)
(*Z dependence of thermal properties*)
kcondz[z_] := 
 kcond[[1]] + 
   Sum[(kcond[[j]] - kcond[[j - 1]])*(1/2 + 
       Tanh[(z - zvalues[[j - 1]])/αp]/2) , {j, 2, 
     noflayers}] /. conds;
ρxc[z_] := ρdens[[1]] cspe[[1]] + 
    Sum[(cspe[[j]] ρdens[[j]] - 
        cspe[[j - 1]] ρdens[[j - 1]])*(1/2 + 
        Tanh[(z - zvalues[[j - 1]])/αp]/2) , {j, 2, 
  noflayers}] /. conds;
Qlaser[t_, r_, 
   z_] := (Ptotal/Vtotal) UnitStep[
     rlaser - r] UnitStep[-t] Exp[-z/Llaser] /. conds;
(*Penne's bioheat equation in cylindrical coordinates*)
eqn = (D[kcondz[z], z]*D[u[t, r, z], z]) + 
    kcondz[z]*(D[u[t, r, z], z, 
        z] + (1/r) D[r D[u[t, r, z], r], r]) - (ρxc[z] D[
       u[t, r, z], t]) + Qlaser[t, r, z] /. conds;
(*Region studied*)
Ω = 
  ImplicitRegion[True, {{r, 0.001, 1}, {z, 0, zvalues[[noflayers]]}}];
sol = NDSolveValue[{eqn == 
    NeumannValue[(u[t, r, z] - Tenv) (ha)/(kcondz[z]) /. conds, 
     z == 0], 
   DirichletCondition[u[t, r, z] == Tcore /. conds, 
    z == zvalues[[noflayers]]], 
   DirichletCondition[u[t, r, z] == Tcore /. conds, r == 1], 
   u[-lasert, r, z] == Tcore /. conds}, {u, 
RowBox[{"0", ",", "0", ",", "1"}], ")"}],
MultilineFunction->None]\)}, {t, -lasert /. conds, 
   1500}, {r, z} ∈ Ω, 
  Method -> {"MethodOfLines", 
    "SpatialDiscretization" -> {"FiniteElement"}}]

NDSolveValue indeed finds a solution. The problem is, if you try to plot tdiff[z]*D[sol[t,r,z],z] (that is, the heat flux), you will find a lot of discontinuities (see figure below). Therefore this is not a physical solution to the problem.

Plot[-kcondz[z] Evaluate[sol[[2]][0, 0.5, z]], {z, 0, 
  zvalues[[noflayers]]}, PlotRange -> All]

Discontinuity of heat flux

Can somebody please show me a way to impose the boundary conditions concerning the heat flux continuity at the interfaces of the layers?


  • $\begingroup$ The continuity of heat flux can be imposed by making a change of variable. You just need to make the change $k ξ = z$. I'm a little busy at the moment so can't write an detailed answer this week. 8楼 of this post is an answer for a similar problem. (It's written in Chinese, but I think you can at least read the code. ) $\endgroup$
    – xzczd
    Jul 16, 2016 at 10:35
  • $\begingroup$ @xzczd I will give a look at this post. I will try to use Google Translator. Thanks! $\endgroup$
    – J.Edwards
    Jul 16, 2016 at 11:24
  • $\begingroup$ @xzczd I was thinking about the change of variable that you suggested but, since the relation between $\xi$ and $z$ is not unequivocal, I believe there will be problems if you have a composite material with more than two layers. For instance, in our problem $\xi$ assumes the value of 0.0221569 for both z = 0.0113 and z = 0.004102. Therefore, the condition $T_{z=0.0113}=T_{wall}$ cannot be expressed in terms of $\xi$. Am I right? $\endgroup$
    – J.Edwards
    Jul 17, 2016 at 11:28
  • $\begingroup$ I made a mistake in my last comment. The relationship between $z$ and $\xi$ should be $dz=k d\xi$ i.e. $z=k\xi + b$. (Notice $b$ is probably piecewise, too.) $\endgroup$
    – xzczd
    Jul 18, 2016 at 4:32
  • $\begingroup$ @xzczd Thanks! I will try to write something along these lines. $\endgroup$
    – J.Edwards
    Jul 18, 2016 at 10:10

1 Answer 1


This post contains several code blocks, you can copy them easily with the help of importCode.

Pre-v10 Soluion

Note: If you're in or after v10, it's recommended to turn to the solution in Update: "FiniteElement"-based Solution section of this answer. That one is much simpler and more efficient.

Finally I managed to find some time to write this answer. The key idea of this solution is, as said in the comment above, making the change $dz=kd\xi$ to ensure the continuity of heat flux.

…OK, I think I need to explain why the transformation $dz=kd\xi$ will ensure the continuity of heat flux first. Well, actually this has been mentioned in the Scope of document of NDSolve but a little vague. Anyway, have you ever noticed the behavior of NDSolve when it faces a piecewise coefficient? If not, check the following examples.

Example 1, 1st order ODE, with piecewise coeffcient

coef = Piecewise[{{1, x <= 2}, {-1, x > 2}}];
sol1 = NDSolveValue[{coef y'[x] == x, y[0] == 1}, y, {x, 0, 3}];
Plot[sol1[x], {x, 0, 3}]

enter image description here

As shown above, NDSolve outputs a 0 order continuous solution.

Example 2, 2nd order ODE, with piecewise coeffcient

sol2 = NDSolveValue[{coef y''[x] == x, y[0] == 1, y'[0] == -1}, y, {x, 0, 3}];
Plot[{sol2[x], sol2'[x]}, {x, 0, 3}]

enter image description here

As shown above, NDSolve outputs a 1st order continuous solution.

Example 3, 3rd order ODE, with piecewise coeffcient

sol3 = NDSolveValue[{coef y'''[x] == x, y[0] == 1, y'[0] == -1, y''[0] == 0}, 
   y, {x, 0, 3}];
Plot[{sol3[x], sol3'[x], sol3''[x]}, {x, 0, 3}]

enter image description here

As shown above, NDSolve outputs a 2nd order continuous solution.

Example 4, 1D heat conduction equation, with piecewise coeffcient for the $\frac{\partial ^2u}{\partial x^2}$ term

coef4 = Piecewise[{{1, x <= 3/2}, {3, x > 3/2}}];
sol4 = NDSolveValue[{D[u[t, x], t] == coef4 D[u[t, x], x, x], u[0, x] == x (3 - x), 
    u[t, 0] == 0, u[t, 3] == 0}, u, {t, 0, 2}, {x, 0, 3}, MaxStepSize -> 0.01];
Plot[{sol4[2, x], Derivative[0, 1][sol4][2, x]}, {x, 0, 3}]

enter image description here

As shown above, NDSolve outputs a solution that is 1st order continuous in x-direction.

To sum up, for differential equation involving n order derivative term with a piecewise coefficient, NDSolve will give a solution which is n-1 order continuous in the corresponding direction.

Back to your question, after the transform $dz=kd\xi$, the dependent variable of your equation will become

$$U(t, r, \xi)=u(t,r,z(\xi))$$

Notice I've used a new variable $U$ to strictly distinguish the function relationship.

According to the illustrations above, after the transformation, $\frac{\partial U(t,r,\xi)}{\partial \xi}$ will be continuous in the output of NDSolve, while

$$\frac{\partial U(t,r,\xi)}{\partial \xi}=\frac{\partial u(t,r,z(\xi))}{\partial z(\xi)}\frac{\partial z(\xi)}{\partial \xi}=k \frac{\partial u(t,r,z(\xi))}{\partial z(\xi)}$$

So, making the transformation is equivalent to impose the heat flux continuity condition.

OK, analysis finished, let's programme.

Notice I've made some tiny modifications to dvalues and zvalues just to make it conciser:

conds = {αp -> 0.8/10^4, Tenv -> 22, Tcore -> 35, ha -> 12, 
       Ptotal -> 1.3, Vtotal -> 2.25/10^6, rlaser -> 0.13, 
       Llaser -> 4.5/10^3, lasert -> 3 60}; 
dvalues = {0, 1/10^3, 3/10^3, 3/10^3, 1/10^3, 6/10^3}; 
kcond = {0.235, 0.445, 0.445, 0.185, 0.51}; 
ρdens = {1200, 1200, 1200, 1000, 1085}; 
cspe = {3589, 3300, 3300, 2674, 3800}; 
zvalues = Accumulate[dvalues]; 
Qlaser[t_, r_, z_] = (Ptotal UnitStep[rlaser - r] UnitStep[-t] 
            Exp[-(z/Llaser)])/Vtotal /. conds; 

Your original definition for kcondz and ρxc contains Tanh, which seems to be a smoothing for the function. As far as I can tell, this somewhat complicated definition isn't necessary for obtaining a solution and stops one from making an analytic transform. (One can turn to numeric transform instead, of course. ) So I chose the simpler Piecewise approach:

coefuncgenerator[cond_List] := 
 Piecewise[MapThread[{#1, #2[[1]] < z <= #2[[2]]} &, {cond, Partition[zvalues, 2, 1]}], 

kcondz[z_] = coefuncgenerator[kcond]; 
ρxc[z_] = coefuncgenerator[ρdens cspe]; 

Plot[ρxc[z], {z, 0, Last[zvalues]}, PlotRange -> All, Exclusions -> None]

enter image description here

Now let's find the relationship between z and ξ:

ξfunc = ξ /. First[DSolve[{ξ'[z] == 1/kcondz[z], ξ[0] == 0}, ξ, z]]; 

It's a pity that at least currently DChange can't do the change of variable for us, so I make the transformation in a more primitive way:

zfunc[ξ_] = 
 Piecewise@Transpose@{z /. Solve[ξ == #, z][[1]] & /@ ξfunc[z][[1, All, 1]], 
      ξfunc[z][[1, All, -1]] /. {z -> ξ, a_?NumericQ :> ξfunc@a}};

Plot[zfunc@ξ, {ξ, 0, ξfunc@Last@zvalues}]

enter image description here


You can also find relationship between z and ξ numerically in the following way:

ξfunc = 
  NDSolveValue[{ξ'[z] == 1/kcondz[z], ξ[0] == 0}, ξ, {z, 0, Last@zvalues}];

derivativeOnGrid = 
  Sequence @@ (Partition[#[[-2, -1]], #[[2, 3]] + 1]\[Transpose] // Rest) &;

zfunc = ListInterpolation[{#["Coordinates"][[1]], 1/derivativeOnGrid@#}\[Transpose], {#[
       "ValuesOnGrid"]}] &@ξfunc;

Notice this approach will make the last step i.e. differential equation solving slower.

ρxcnew[ξ_] = ρxc@zfunc@ξ // Simplify;
kcondznew[ξ_] = kcondz@zfunc@ξ // Simplify;

neweqn = kcondznew[ξ] (D[u[t, r, ξ], ξ, ξ]/kcondznew[ξ]^2 + 
            D[r D[u[t, r, ξ], r], r]/r) - ρxcnew[ξ] D[u[t, r, ξ], t] + 
       Qlaser[t, r, zfunc[ξ]]; 

Notice the D[kcondz[z], z]*D[u[t, r, z], z] in the original equation is taken away because it's zero in this case.

I'm not in v10 yet so can't use "FiniteElement" method. In v9 I added an extra condition Derivative[0, 1, 0][u][t, 10^-3, ξ] == 0 to suppress the warning NDSolveValue::bcart. Though NDSolveValue still spits out the warning NDSolveValue::ibcinc, it gives an answer:

   sol = NDSolveValue[{neweqn == 0, 
           ha (u[t, r, ξ] - Tenv) == D[u[t, r, ξ], ξ] /. ξ -> 0, 
           u[t, r, ξfunc[Last[zvalues]]] == 
               Tcore, Derivative[0, 1, 0][u][t, 10^(-3), ξ] == 0, 
           u[t, 1, ξ] == Tcore, 
           u[-lasert, r, ξ] == Tcore} /. conds, 
         {t, -lasert /. conds, 0}, {r, 1/10^3, 1}, {ξ, 0, ξfunc@Last@zvalues}, 
Method -> {"MethodOfLines", 
  "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> {25, 30}, 
    "MinPoints" -> {25, 30}, "DifferenceOrder" -> 2}, 
  "DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 1}}]; ]
(* {103.8961825, Null} *)

(* Temperature *)
Plot[sol[0, 0.5, ξfunc[z]], {z, 0, Last[zvalues]}, PlotRange -> All]

Mathematica graphics

(* Heat flux *)
Plot[Derivative[0, 0, 1][sol][0, 0.5, ξfunc[z]], {z, 0, Last[zvalues]}, 
   PlotRange -> All]

Mathematica graphics

Notice the option "DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 1} is necessary here, because NDSolve will automatically set "ScaleFactor" -> 0 for non-Dirichlet boundary condition, when i.c. and b.c. are inconsistent, this will make the non-Dirichlet inconsistent b.c. be severely ignored, for more information you may want to read this obscure tutorial and this answer.

Feel free to ask in the comment if you have difficulty in understanding the answer.

Update: "FiniteElement"-based Solution

The last written set of boundary conditions indicates the continuity of the heat flux, ……, They are the ones I don't know how to impose in NDSolveValue.

After several discussions (1, 2) in this site during past 2 years, it has been clear that "FiniteElement" method introduced since v10 is able to take care of heat flux continuity automatically, as long as the form of the equation is proper. (Check links above for more information. ) And your trial turns out to be almost successful, except that:

  1. When dealing with a region with internal boundaries, it's necessary to specify them when generating the mesh, or the quality of the solution will be bad. Here's a simpler example illustrating this issue.

  2. Your translation for the Robin type b.c. is incorrect, the correct one should be

    NeumannValue[(u - Tenv) ha, z == 0] 

    As to the reason, please check Details section of document of NeumannValue.

  3. The Tanh term aiming at smoothing the discontinuity seems to be troublesome.

OK, let's fix these issues.

Note: Definitions for zvalues, kcondz, etc. are the same as those in Pre-v10 Soluion section and omitted in this section.

We first generate a mesh with internal boundaries:

r0 = 0; rend = 1;
line = Flatten[{{{r0, #}, {rend, #}} & /@ 
     Flatten@{0, zvalues}, {{#, 0}, {#, zvalues[[-1]]}} & /@ {r0, rend}}, 1];
coord = Flatten[line, 1] // Union;
indexrule = MapIndexed[# -> #2[[1]] &, coord];
bmesh = ToBoundaryMesh["Coordinates" -> coord, 
   "BoundaryElements" -> (LineElement[{#}] & /@ (line /. indexrule))];
Show[bmesh[Wireframe], AspectRatio -> 1]
mesh = ToElementMesh[bmesh];
(* Show[mesh[Wireframe], AspectRatio -> 1] *)

Mathematica graphics

Then introduce the Robin b.c. in the correct way, I've also made some simplification in coding:

With[{u = u[t, r, z]}, 
 eqn = kcondz[z] (D[u, z, z] + (1/r) D[r D[u, r], r]) - ρxc[z] D[u, t] + 
   Qlaser[t, r, z];
 robinbc = NeumannValue[(u - Tenv) ha, z == 0];
 bc = DirichletCondition[u == Tcore, r == 1 || z == zvalues[[-1]]];
 ic = u == Tcore /. t -> -lasert;]

Finally, solve the equation and plot. We can see the result is the same as the previous, but NDSolveValue only takes about 45 seconds to finish:

tend = 1500;
sol = NDSolveValue[{eqn == robinbc, bc, ic} /. conds, {u, 
     Derivative[0, 0, 1][u]}, {t, -lasert, tend} /. conds, {r, z} ∈
     mesh]; // AbsoluteTiming
(* {45.6404, Null} *)
Plot[sol[[1]][0, 0.5, z], {z, 0, zvalues[[-1]]}, PlotRange -> All]

Mathematica graphics

Plot[kcondz[z] Evaluate[sol[[2]][0, 0.5, z]], {z, 0, zvalues[[-1]]}, PlotRange -> All]

Mathematica graphics

Method -> "FiniteElement" isn't set in NDSolveValue, because it'll be automatically chosen when NeumannValue is used to build the equation.

  • $\begingroup$ What does the new variable \[Xi] represent physically? Why would it cause the heat flux to be continous? $\endgroup$ Aug 1, 2016 at 15:14
  • $\begingroup$ In the neweqn with the change of variable I don't see the term corresponding to D[kcondz[z], z]*D[u[t, r, z], z] that was present in the original equation. $\endgroup$ Aug 1, 2016 at 15:16
  • $\begingroup$ @JackLaVigne The term D[kcondz[z], z]*D[u[t, r, z], z] is actually zero in this case so I simply omit it. (I believe OP has derived the equation in a somewhat blind way :) ). $\endgroup$
    – xzczd
    Aug 2, 2016 at 6:22
  • 1
    $\begingroup$ Actually I have included the term $\frac{d k}{dz} \frac{\partial u}{\partial z}$ because of the k dependence with z (remember it had a Tanh to smooth it). $\endgroup$
    – J.Edwards
    Aug 2, 2016 at 10:46
  • 1
    $\begingroup$ @J.Edwards Recently I revisit the tutorial for "MethodOfLines" and noticed an additional option adjustment is needed for obtaining a reasonable result. See my edit. $\endgroup$
    – xzczd
    Sep 28, 2016 at 6:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.