Solve the heat equation for a semi infinite rod considering convection

Bug introduced in or after 10.3, persisting through 11.2.

I'm trying to solve following PDE (heat equation):

$$\begin{cases} u_t = a \, u_{xx} \\ u(x,0) = 0\\ \lim_{x\to \infty}u(x,t) =0\\ \alpha\, (\theta_0-u(0,t))+\dot{q}_0=-\lambda u_x(0,t) \end{cases}$$

Where basically I have an initial temperature of $0$ everywhere, a constant heat flux at the beginning of the rod, and convection between air and the rod at its beginning ($\theta_0$ is the air temperature which is assumed to be constant).

I found following analytical solution for the problem:

$$u(x,t) = \frac{\dot{q}_0+\alpha \, \theta_0}{\alpha}\left[ \mathrm{erfc}\, \left(\frac{x}{2\sqrt{a \, t}} \right) -\mathrm{exp}\,\left(\frac{\alpha}{\lambda}x+a \frac{\alpha^2}{\lambda^2}t \right)\mathrm{erfc}\,\left(\frac{x}{2\sqrt{a\, t}} + \frac{\alpha}{\lambda} \sqrt{a\, t} \right) \right]$$

which is physically meaningful. With mathematica, however, I get some meaningless results, probably due to my not that good mathematica skills.

This is what I tried to do:

λ = 0.8; c = 880; ρ = 1950; a = λ/(c ρ)
α = 15; θair = 0;

heqn = D[u[x, t], t] == a D[u[x, t], {x, 2}];
ic1 = u[x, 0] == 0;
bc1 = α (θair - u[0, t]) + 650 == -λ Derivative[1, 0][u][0, t];

sol = DSolveValue[{heqn, ic1, bc1}, u[x, t], {x, t}][[1, 1, 1]]

Which leads me to a complex solution (plotted below)

DensityPlot[sol, {t, 0, 2*3600}, {x, 0, 0.1},
PlotLegends -> Automatic, FrameLabel -> {t[s], x[m]}] [ Plot[Evaluate[Table[sol, {t, 3600, 7200, 3600}]], {x, 0, .1},
PlotRange -> All, Filling -> Axis, Axes -> True, AxesLabel -> {x[m], θ[°C]}] I'm aware of the fact that I didn't consider the condition at infinity. To do this, I tried to follow this answer without success. Also, mathematica finds the solution without this condition as soon as $\alpha = 0$.

This is the plot of the analytical solution I get:

Plot[{u[x, 600], u[x, 3600], u[x, 7200]}, {x, 0, .2}, Filling -> Axis,
Axes -> True, AxesLabel -> {x[m], θ[°C]}] • Do you want to obtain the same analytic solution with Mathematica, or just want to solve the problem numerically or analytically? – xzczd Jan 17 '18 at 11:08
• BTW, something seems to be wrong with the $\LaTeX$ formula, I failed to reproduce the graph with it: sola[t_, x_] = ( q + α θair)/α (Erfc[x/(2 Sqrt[a t])] - Exp[α/λ x + a α^2/λ^2 t] Erfc[ x/(2 Sqrt[α t]) + α/λ Sqrt[a t]]) – xzczd Jan 17 '18 at 11:26
• @xzczd at the moment I cannot check if I did some mistake by copying the analytical solution. I'll check it as soon as possible (probably in 1-2hours). Regarding your second question: I'm preparing a benchmark problem for a FEM simulation. A numerical solution is therefore in principle also ok. The only thing that seemed a bit strange to me is that Mathematica cannot solve a PDE which can be solved analytically (although I used Laplace transformations and this isn't probably the easiest thing for a software) – David Jan 17 '18 at 11:44
• @xzczd The latex formula is correct, you had a small mistake in your definition of the term Erfc[ x/(2 Sqrt[α t]) + α/λ Sqrt[a t]]. It should be x/(2 Sqrt[a t]) instead of x/(2 Sqrt[α t]). – David Jan 17 '18 at 13:35
• Oh, I see… I suggest adding the Mathematica code for the formula to the question so people can easily test it. BTW, thx for accepting, but you actually don't need to accept that fast, you can wait for 24 hours or so to see if someone will come up with a better answer :) . – xzczd Jan 17 '18 at 13:39

Semi-analytic Solution

The approach in the linked post can be used for solving your problem analytically. We just need an extra step i.e. making Laplace transform:

Clear@"`*"
heqn = D[u[x, t], t] == a D[u[x, t], {x, 2}];
ic1 = u[x, 0] == 0;
ic2 = α (θair - u[0, t]) + q == -λ Derivative[1, 0][u][0, t];

teqn = LaplaceTransform[{heqn, ic2}, t, s] /. Rule @@ ic1 /.
HoldPattern@LaplaceTransform[a_, __] :> a

tsol = Collect[u[x, t] /. First@DSolve[teqn, u[x, t], x], Exp[_], Simplify]
(*
E^((Sqrt[s] x)/Sqrt[a]) C + (E^(-((Sqrt[s] x)/Sqrt[
a])) (s^(3/2) λ C + Sqrt[a] (q + α (θair - s C))))/(
Sqrt[a] s α + s^(3/2) λ)
*)

tsol is the transformed solution, with one constant to be determined. Since the solution should vanish at $\infty$, C can only be 0:

tsolfunc[s_, x_] = tsol /. C -> 0

The last step is to transform back, but sadly InverseLaplaceTransform can't handle it. Nevertheless, we've indeed obtained an analytic solution involving integration: $$u=\frac{1}{2 \pi i } \int_{\gamma-i \infty }^{\gamma+i \infty } \frac{\sqrt{a} e^{-\frac{\sqrt{s} x}{\sqrt{a}}} (\alpha \theta_0+q)}{\sqrt{a} \alpha s+\lambda s^{3/2}} e^{st} \, ds$$

If we want to get the numeric result, we need to make use of numeric inverse Laplace transform package e.g. this:

λ = 8/10; c = 880; ρ = 1950; a = λ/(c ρ);
α = 15; θair = 0; q = 650;
(* Definition of FT isn't included in this post,
please find it in the link above. *)
Clear@sol; sol[t_, x_] = FT[tsolfunc[#, x] &, t];

Plot[sol[#, x] & /@ {600, 3600, 7200} // Evaluate, {x, 0, .2}, Filling -> Axis,
AxesLabel -> {x[m], θ[°C]}] Analytic Solution

Based on this answer, the following is one way to find the symbolic Laplace inversion:

coe = (q + α θair);
mid = Apart[LaplaceTransform[tsolfunc[s, x]/coe, x, S] /. s -> sqrts^2] /.
sqrts -> Sqrt[s]

mid2 = InverseLaplaceTransform[mid, s, t]; // AbsoluteTiming
(* {33.611455, Null} *)
final = FullSimplify[InverseLaplaceTransform[mid2[] // Expand, S, x], {t > 0, a > 0}]

coe final is the analytic solution:

coe final

$$u=\frac{(\alpha \theta_0+q)}{\alpha } \left(\text{erfc}\left(\frac{x}{2 \sqrt{a t}}\right)-e^{\frac{\alpha (a \alpha t+\lambda x)}{\lambda ^2}} \text{erfc}\left(\frac{2 a \alpha t+\lambda x}{2 \lambda \sqrt{a t}}\right)\right)$$

Check:

sola[t_, x_] = (q + α θair)/α (Erfc[x/(2 Sqrt[a t])] -
Exp[α/λ x + a α^2/λ^2 t] Erfc[x/(2 Sqrt[a t]) + α/λ Sqrt[a t]])

sola[t, x] == final coe // FullSimplify
(* True *)

Numeric Solution

NDSolve can also be used for solving the problem, of course. We just need to adjust the option a little:

inf = 0.2;

mol[tf : False | True, sf_: Automatic] := {"MethodOfLines",
"DifferentiateBoundaryConditions" -> {tf, "ScaleFactor" -> sf}}

nsol = NDSolveValue[{heqn, ic1, ic2, u[inf, t] == 0}, u, {x, 0, inf}, {t, 0, 7200},
Method -> mol[True, 100]];

np = Plot[nsol[x, #] & /@ {600, 3600, 7200} // Evaluate, {x, 0, inf}, Filling -> Axis,
AxesLabel -> {x[m], θ[°C]}}]

The result looks the same so I'd like to omit it here. To learn more about why the option needs to be adjusted, check this post.

• Thanks for your answer! I accepted it after getting the correct semi-analytical solution. However, now I'm trying the numerical one and I still get a NDSolveValue::ibcinc error (bc and ic are inconsistent), which, according to your other post should be avoided by using the method of lines. – David Jan 17 '18 at 16:34
• @david Well, the linked post is not about avoiding ibcinc, it's about getting an accurate enough result even if it pops up. Notice ibcinc is merely a warning, it doesn't necessarily mean the result will be unreliable. （You can compare np to the analytic solution. ） – xzczd Jan 18 '18 at 2:52
• Sure, you're right. thanks – David Jan 18 '18 at 7:39

Analytical solution seems to be wrong!

Just try

heqn = Derivative[0, 1][u][x, t] == a Derivative[2, 0][u][x, t];
ic1 = u[x, 0] == 0;
ic2 = \[Alpha] (\[Theta]air - u[0, t]) +qp0 == -\[Lambda] Derivative[1, 0][u][0, t];

sol = DSolve [{heqn, ic1, ic2}, u[x, t], {x, t}][]
(*{u[x, t] -> -((\[Alpha] (qp0 + \[Alpha] \[Theta]air) ((2 \[Alpha] +x \[Lambda]) Erfc[x/(2 Sqrt[a t])] -
2 E^((\[Lambda] (x \[Alpha] +a t \[Lambda]))/\[Alpha]^2) \[Alpha] Erfc[(
x \[Alpha] + 2 a t \[Lambda])/(2 Sqrt[a t] \[Alpha])]))/(2\[Lambda]^3))} *)

sol looks like the expected analytical solution!

u[infinity,t]->0 seems to be fullfilled:

Simplify[Normal[ Series[u[x, t] /. sol, {x, Infinity, 0}] ], {x > 0,t > 0, \[Alpha] > 0, \[Lambda] > 0, a > 0} ]
(*-((E^(-(x^2/(4 a t))) Sqrt[a t] \[Alpha] (qp0 + \[Alpha] \[Theta]air))/(Sqrt[\[Pi]] \[Lambda]^2))*)

because E^(-(x^2/(4 a t))) decreases with x!

The boundary condition u[x,0]==0 can be verified(a>0,x>0):

Simplify[Limit[u[x, t] /. sol[], t -> 0], {x > 0, a > 0}]
(* 0 *)

but the solution doesn't fullfill the pde

(E^(-(x^2/(4 a t))) Sqrt[a t] \[Alpha] (qp0 + \[Alpha] \[Theta]air))/\[Lambda] == 0

Also ic2 isn't true!

• Sadly it's not… you'll see the 2nd graph in the question if you plot it. – xzczd Jan 17 '18 at 12:15
• @xzczd hopefully it does. I edited my answer... – Ulrich Neumann Jan 17 '18 at 21:45
• Well, but it doesn't satisfy the b.c. at $x=0$, try test = {u[x, t] -> -((\[Alpha] (qp0 + \[Alpha] \[Theta]air) ((2 \[Alpha] + x \[Lambda]) Erfc[x/( 2 Sqrt[a t])] - 2 E^((\[Lambda] (x \[Alpha] + a t \[Lambda]))/\[Alpha]^2) \[Alpha] Erfc[( x \[Alpha] + 2 a t \[Lambda])/(2 Sqrt[a t] \[Alpha])]))/( 2 \[Lambda]^3))} /. (u[x, t] -> rest_) :> (u -> Function @@ {{x, t}, rest}); Simplify[ic2 /. test, t > 0]. You can also plot the solution and compare it to the graph in the question. – xzczd Jan 18 '18 at 3:06
• @xzczd u[x,0}==0 !I edited my answer... – Ulrich Neumann Jan 18 '18 at 8:17
• This is the initial condition… I mean the boundary condition i.e. ic2 in the code. – xzczd Jan 18 '18 at 8:22