Heat transport in 1d, two different materials at point contact with NDsolve -

I have two materials touching at a point x=0.

1) The first material lies within the interval [-L1, 0], has diffusion constant D1, and the heat conductivity kappa1. Its left end is kept on time dependent heat flux Cos[2 omega t]

2) The second material lies within the interval [0,L2], has diffusion constant D2, and the heat conductivity kappa2. Its right end is kept at constant temperature Temp.

Denoting temperatures in the first material as u1[x,t], and in the second as u2[x,t], I want the following continuity at x=0:

Continuity of Temperature: u1[0,t]==u2[0,t]

Continuity of Energy flow: kappa1 Derivative[1,0][u1][0,t]==kappa2 Derivative[1,0][u2][0,t]

The initial condition just requires that both u1[x,0] and u2[x,0] are equal to Temp.

I formalized my problem as follows:

heateqs = {
D[u1[x, t], t] == D1 D[u1[x, t], {x, 2}],
D[u2[x, t], t] == D2 D[u2[x, t], {x, 2}]};

bcs = {
Derivative[1, 0][u1][-L1, t] == Cos[2 omega t],
u1[0,t]==u2[0,t],
kappa1 Derivative[1, 0][u1][0, t] == kappa2 Derivative[1, 0][u2][0, t],
2[L2, t] == Temp
};

ics = {u1[x, 0] == Temp, u2[x, 0] == Temp};

sol = NDSolve[{heateqs, ics, bcs}, {u1, u2}, {x, L1, L2}, {t, 0,
1000}]


PROBLEM: mathematica wrote:

NDSolve::bcedge: "Boundary condition u1[0,t]==u2[0,t] is not specified on a single edge of the boundary of the computational domain."

Seems I can not connect two solutions u1[x,t] and u2[x,t] at x=0, in the naive (and logical) way as I did.

QUESTION: how to properly formalize my problem for NDSolve?

• Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour and check the faqs! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! Oct 27, 2019 at 19:06
• Why not just have one temperature field but define the equations so that they rely on a position dependent diffusion coefficient? Oct 27, 2019 at 19:06
• Mathematica can solve this if you define only one region and consider D as a function of x. The finite element method will automatically respects continuity of temperature and continuity of flux, provided that you express the equations in a certain form (the form eq1, and not the form eq2 here. Oct 27, 2019 at 19:43
• @Dunlop The problem is not at position dependent thermal Diffusion constant, but on different heat conductivities $\kappa_1$ and $\kappa_2$. Those do not enter directly the heat equation, but represents boundary condition. At the junction of two materials, apart of continuity of temeperature $u_1(x=0,t)=u_2(x,t)$ also the heat current should be conserved: $\kappa_1\partial_x u_1(x=0,t)=\kappa_2\partial_x u_2(x=0,t)$. This is a type of boundary-like condition at a mid point (in my case $x=0$). Oct 28, 2019 at 11:01
• Version 12.1 has a tutorial on Heat Transfer Modeling Mar 20, 2020 at 19:51

Update 2 Express Energy Balance in Coefficient Form

The most direct way to introduce $$\kappa$$ is to express the energy (not temperature) balance in coefficient form as described here and shown below.

$$m\frac{{{\partial ^2}}}{{\partial {t^2}}}u + d\frac{\partial }{{\partial t}}u + \nabla \cdot\left( { - c\nabla u - \alpha u + \gamma } \right) + \beta \cdot\nabla u + au - f = 0$$

In this case, only the $$d$$ and $$c$$ term are non-zero and given by:

$$d = \rho {{\hat C}_p} = \frac{\kappa }{D}$$ $$c=\kappa$$

This time, we will make the second phase a dense insulating material.

L1 = -1;
L2 = 1;
D1 = 1;
D2 = 1/160;
k1 = 1;
k2 = 0.1;
rhocp1 = k1/D1;
rhocp2 = k2/D2;
Temp = 0;
nv = NeumannValue[Cos[ Pi t], x == L1];
dc = DirichletCondition[u[t, x] == Temp, x == L2];
ic = u[0, x] == Temp;
pde = If[x < 0, rhocp1, rhocp2] D[u[t, x], t] +
D[-If[x < 0, k1, k2] D[u[t, x], x], x] == nv;
uif = NDSolveValue[{pde, dc, ic}, u, {t, 0, 10}, {x, L1, L2}];
imgs = Table[
Plot[uif[t, x], {x, L1, L2}, PlotRange -> {-0.55, 0.55},
ImageSize -> Medium], {t, 0, 10, 0.1}];
ListAnimate@imgs


Expressing PDEs in coefficient form make it easier to map to other solvers.

Again, Mathematica compares favorably to another solver giving me more confidence in the approach.

Update 1 To Perform Balance on Energy Instead of Temperature

If one wants to use a single state variable for thermal problems with distinct phases, then one should use enthalpy change instead of temperature as the field variable. The conservation equation in terms of enthalpy change becomes:

$$\frac{\partial }{{\partial t}}\rho {{\hat C}_p}T + \nabla \cdot \left( { - \alpha \nabla \rho {{\hat C}_p}T} \right) = 0$$

It is natural to set the field variable, $$u$$, to

$$u = \rho {{\hat C}_p}(T - {T_0})$$

Leading to a single pde for energy

$$\frac{\partial }{{\partial t}}u + \nabla \cdot \left( { - \alpha \nabla u} \right) = 0$$

This was the PDE that was solved previously. One has to think carefully about how to preserve temperature continuity. We should be able to convert back to temperature as a post-processing step by defining a piecewise function like:

$$T(t,x) = {T_0} + \frac{{u(t,x) - u(0,x)}}{{{\rho _1}{{\hat C}_{p1}}}};x < 0$$ $$T(t,x) = {T_0} + \frac{{u(t,0) - u(0,0)}}{{{\rho _1}{{\hat C}_{p1}}}} + \frac{{u(t,x) - u(t,0)}}{{{\rho _2}{{\hat C}_{p2}}}};Otherwise$$

At $$x=0$$, temperature continuity is preserved.

Comparison to Another Code

When possible, it is conducive to compare results among different codes. Here, I compare the Mathematica result to the Heat Transfer in Solids module in COMSOL where $$\rho {{\hat C}_p}$$ is held constant. The results look qualitatively similar.

Mathematica

imgs = Table[
Plot[uif[t, x], {x, L1, L2}, PlotRange -> {-0.5, 0.5},
ImageSize -> Medium], {t, 0, 10, 0.1}];
ListAnimate@imgs


COMSOL

As stated in the comments, you only need one field variable and a spatially diffusion coefficient. Since you did not provide parameter values, I created a case with test values.

L1 = -1;
L2 = 1;
D1 = 1;
D2 = 0.1;
Temp = 0;
nv = NeumannValue[Cos[ Pi t], x == L1];
dc = DirichletCondition[u[t, x] == Temp, x == L2];
ic = u[0, x] == Temp;
pde = D[u[t, x], t] + D[-If[x < 0, D1, D2] D[u[t, x], x], x] == nv;
uif = NDSolveValue[{pde, dc, ic}, u, {t, 0, 10}, {x, L1, L2}];
Plot3D[uif[t, x], {x, L1, L2}, {t, 0, 10}, PlotPoints -> 50,
PlotRange -> All]


• Thanks for your answer, but I think this is not an answer to my question. Your solution assumes different thermal diffusivity constants $D_1$ and $D_2$, but the same thermal conductivities $\kappa_1$ and $\kappa_2$. At least, I do not see them entering your code. In my case, at the point $x=0$ there is a breakdown and $\kappa_1 \neq \kappa_2$, hence according to the Fourier law one needs solution satisfying at $x=0$: $$\kappa_1 \partial_x u_1(x=0,t) = \kappa_2 \partial_x u_2(x=0,t)$$ plus continuity $$u_1(x=0,t)=u_2(x=0,t)$$ What I need is to know, how to implement those two conditions. Oct 28, 2019 at 7:56
• As far as I know, internal BCs are not generally supported in MMA. MMA does, however, support region dependent physical properties well. I updated my answer to show how to avoid the internal BC by performing the balance on energy versus temperature. You should be able to convert the energy back to temperature as a post-processing step using the piecewise function shown. Oct 28, 2019 at 15:53
• Thanks a lot Tim! What you have posted is, indeed, very helpful! This is a different strategy I was not thinking about. I will go for it. Oct 28, 2019 at 16:55
• @zelenysip You are welcome. Handling interfacial phenomena can be extremely tricky, especially in moving boundary problems, like solidification, with latent heat and hysteresis effects. Although MMA may not yet support a particular BC, one can often find a suitable work around with all the other features they have to offer. Good luck! Oct 28, 2019 at 20:09
• @andre314 I followed your recommended edit. Thanks! Oct 28, 2019 at 20:11