# Problem with boundary condition 2D heat transfer

I want solve this 2D Heat transfer:

So at Mathematica:

L = 1;
k = 237;
To = Tinf = 10;
α = 80*10^(-6);
q = 10;
h = 25;

s[vx_, vy_,
vt_] := (T[x, y, t] /.
Flatten[NDSolve[{1/\[Alpha]*D[T[x, y, t], t] ==
D[T[x, y, t], {x, 2}] + D[T[x, y, t], {y, 2}],
\!$$\*SuperscriptBox[\(T$$,
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[0, y, t] == 0,
\!$$\*SuperscriptBox[\(T$$,
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[L, y, t] == 0,
\!$$\*SuperscriptBox[\(T$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, 0, t] == -q/k, k*
\!$$\*SuperscriptBox[\(T$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, L, t] + h*T[x, L, t] == h*Tinf,
T[x, y, 0] == To},
T[x, y, t], {x, 0, L}, {y, 0, L}, {t, 0, 10}]]) /. {x -> vx,
y -> vy, t -> vt}


When I try to calculate T:

s[L/2, L/2, 0]

NDSolve::ibcinc: Warning: boundary and initial conditions are inconsistent.


I already tried to specify the Method of discretization the spacial derivatives but doesnt works.

• What's the question here? You don't want to see the warning, or the result is not expected, or something else? Have you read this post?: mathematica.stackexchange.com/a/127411/1871 – xzczd Jun 29 at 3:32
• I have created an "exemplary" package HeatTrans for 2D heat transfer model with convective boundary condition. With version 2.0.0 you can easily compare results of NDSolve and AceFEM. Maybe you find it interesting. – Pinti Jul 1 at 8:14

You need to use NeumannValue.

L = 1;
k = 237;
To = Tinf = 10;
q = 10;
h = 25;

R = DiscretizeRegion[Rectangle[{0, 0}, {L, L}], MaxCellMeasure -> 0.0001];

Tsol = NDSolveValue[{
D[T[x, y, t], t] == Laplacian[T[x, y, t], {x, y}]
+ NeumannValue[h Tinf / k - h T[x, y, t] / k, y == L]
+ NeumannValue[-q / k, y == 0]
+ NeumannValue[0, x == 0 || x == L],
T[x, y, 0] == To
}, T, {x, y} ~Element~ R, {t, 0, 10}]

Plot3D[Tsol[x, y, 10], {x, y} ~Element~ R]


The computer I'm using right now can't handle the alpha parameter, but you can just add that back in.

Also because of the way boundary conditions on 2D PDEs work, you might need to adjust the signs of the NeumannValue terms to get what you want.

• If you use ToElementMesh you get a second order mesh; DiscretizeRegion will create a first order mesh. – user21 Jun 29 at 10:56
• What do you mean change the signs of the NeumannValue? If I change the signs I will change the boundary condition. – Mateus Jun 29 at 17:51
• @Mateus Just that the NeumannValue is specified as the gradient of your function dotted with the outward unit normal of the surface. Depending on the side of the rectangle you're on this normal will be pointing in the negative direction. – ccosm Jun 29 at 20:25
• Oh right! I got it! Thanks. – Mateus Jun 29 at 21:15

Your flux boundary conditions are stated in terms of energy fluxes, but your dependent variable is temperature. To avoid confusion, you want to solve the energy balance equation versus the temperature equation. We will assume that temperature is proportional to energy via

$$u = \rho {\hat C_p}T$$

You should also note that having no flux conditions at the x extremes essentially makes this a 1-d problem. Your conductivity parameter suggest that material might be aluminum. We will assume that the density of Al is 2700 kg/cu m. We can now solve for the heat capacity to be 1097.22.

If you are trying to solve for the thermal penetration of 1 m Al block in 10 s, we should not expect the thermal front to penetrate too far. I turned ccosm's answer into a 1-d problem and converted it to an energy balance. You might also note that it is often best to convert your PDE into coefficient form, especially when considering Robin conditions (i.e., a heat transfer coefficient).

Needs["NDSolveFEM"]
\[Alpha] = 80*10^(-6);
rho = 2700;
cp = 1097.22;
k = 237;
L = 1;
To = Tinf = 10;
q = 10;
h = 25;
R = ToElementMesh[Interval[{0, L}], MaxCellMeasure -> 0.001];
op = D[u[t, x], t] + D[-\[Alpha] D[u[t, x], x], x];
rightflow = NeumannValue[h Tinf - h u[t, x]/rho/cp, x == L];
leftflow = NeumannValue[q, x == 0];
ic = u[0, x] == rho cp To;
ufunHeat =
NDSolveValue[{op == rightflow + leftflow, ic}, u,
x \[Element] R, {t, 0, 10}];
plots = Table[
Plot[ufunHeat[t, x]/rho/cp, {x, 0, 1},
PlotRange -> {9.9999, 10.0012}], {t, 0, 10, 0.1}];
ListAnimate[plots]


I verified the solution in COMSOL and the solutions seem similar as shown below.

• If you use ToElementMesh you get a second order mesh; DiscretizeRegion will create a first order mesh. – user21 Jun 29 at 10:56
• @user21 Excellent point. Typically, I would use ToElementMesh. I was fighting some numerical issues and they went away when I used ToElementMesh. I will update my post. – Tim Laska Jun 29 at 11:33
• Why you change the signal of the leftflow? Why q instead -q like the original boundary condition? – Mateus Jun 29 at 18:04
• @Mateus It is possible that I made a mistake. I will look into it in a bit. One thing to consider is that Neumann conditions by convention are positive when the flux is directed into the domain as I discussed in this answer. If you are heating leftflow should be positive and negative for cooling. – Tim Laska Jun 29 at 22:02

Analytic soln is tedious but straightforward.

pde = D[T[x, y, t], x, x] + D[T[x, y, t], y, y] - D[T[x, y, t], t]/α == 0;


Separate variables;

T[x_, y_, t_] = X[x] Y[y] Tt[t];

pde/T[x, y, t] // Expand
(*-(Tt'[t]/(α Tt[t]))+X''[x]/X[x]+Y''[y]/Y[y]==0*)


First use a zero constant

x0eq = X''[x]/X[x] == 0;

DSolve[x0eq, X[x], x] // Flatten;

x0 = X[x] /. % /. {C[1] -> c1, C[2] -> c2}
(*c1+c2 x*)


Use a negative constant for sinusoidal solutions

xeq = X''[x]/X[x] == -a^2;

DSolve[xeq, X[x], x] // Flatten;

x1 = X[x] /. % /. {C[1] -> c3, C[2] -> c4}
(*c3 Cos[a x]+c4 Sin[a x]*)


Now the y's

y0eq = Y''[y]/Y[y] == 0;

DSolve[y0eq, Y[y], y] // Flatten;

y0 = Y[y] /. % /. {C[1] -> c5, C[2] -> c6}
(*c5+c6 y*)

yeq = Y''[y]/Y[y] == -b^2;

DSolve[yeq, Y[y], y] // Flatten;
(*c7 Cos[b y]+c8 Sin[b y]*)

t0eq = -(Tt'[t]/(α Tt[t])) + X''[x]/X[x] + Y''[y]/Y[y] == 0 /. {X''[x]/X[x] -> 0, Y''[y]/Y[y] -> 0};

DSolve[t0eq, Tt[t], t] // Flatten;

Tt0 = 1;

teq = -(Tt'[t]/(α Tt[t])) + X''[x]/X[x] + Y''[y]/Y[y] == 0 /. {X''[x]/X[x] -> -a^2, Y''[y]/Y[y] -> -b^2};

DSolve[teq, Tt[t], t] // Flatten;

Tt1 = Tt[t] /. % /. C[1] -> 1 // Simplify
(*E^(α t (-(a^2+b^2)))*)


Put the solutions together

T[x_, y_, t_] = x0 y0 Tt0 + x1 y1 Tt1
(*E^(α t (-(a^2+b^2))) (c3 Cos[a x]+c4 Sin[a x]) (c7 Cos[b y]+c8 Sin[b y])+(c1+c2 x) (c5+c6 y)*)


Apply the bc's on x

(D[T[x, y, t], x] /. x -> 0) == 0
(*a c4 E^(α t (-(a^2+b^2))) (c7 Cos[b y]+c8 Sin[b y])+c2(c5+c6 y)==0*)


From which

c4 = 0;
c2 = 0;

T[x, y, t]
(*c3 Cos[a x] E^(α t (-(a^2+b^2))) (c7 Cos[b y]+c8 Sin[b y])+c1 (c5+c6 y)*)


Combine constants

c3 = 1;
c1 = 1;


The next bc on x

(D[T[x, y, t], x] /. x -> L) == 0
(*-a Sin[a L] E^(α t (-(a^2+b^2))) (c7 Cos[b y]+c8 Sin[b y])==0*)


Can satisfy this by

a = (n π)/L;

\$Assumptions = n ∈ Integers


The bc's on y

(D[T[x, y, t], y] /. y -> 0) == -(q/k)
(*b c8 Cos[(π n x)/L] E^(α t (-(b^2+(π^2 n^2)/L^2)))+c6==-(q/k)*)


From which

c8 = 0;
c6 = -(q/k);

((k D[T[x, y, t], y] + h T[x, y, t] /. y -> L) // FullSimplify) ==  h Tinf
(*c7 Cos[(π n x)/L] E^(α t (-(b^2+(π^2 n^2)/L^2))) (h Cos[b L]-b k Sin[b L])+c5 h-(q (h L+k))/k==h Tinf*)


We can set the first term to 0 by requiring:

beq = h Cos[b L] - b k Sin[b L] == 0;


so there will be an infinite number of b's, requiring an infinite sum in the final solution. save for later simplification

bRule = {Sin[b L] -> h/(b k) Cos[b L]};


Use the other terms to solve for c5

c5 = c5 /. Solve[c5 h - (q (h L + k))/k == h Tinf, c5][[1]];

T[x, y, t]
(*c7 Cos[b y] Cos[(π n x)/L] E^(α t (-(b^2+(π^2 n^2)/L^2)))+(h k Tinf+h L q+k q)/(h k)-(q y)/k*)


Now we only need to use the ic to solve for c7.

T[x, y, 0] == T0
(*c7 Cos[b y] Cos[(π n x)/L]+(h k Tinf+h L q+k q)/(h k)-(q y)/k==T0*)


Use orthogonality. Shift everything not containing c7 to the RHS, multiply by Cos[(π n x)/L]Cos[b y] and integrate.

c7*Integrate[Cos[(Pi*n*x)/L]^2*Cos[b*y]^2, {y, 0, L}, {x, 0, L}] ==
Integrate[(T0 - (h*k*Tinf + h*L*q + k*q)/(h*k) + (q*y)/k)*Cos[(Pi*n*x)/L]*Cos[b*y],
{y, 0, L}, {x, 0, L}]
(*(c7 L (2 b L+Sin[2 b L]))/(8 b)==0*)


We get 0 for general n. This means only n = 0 will contribute and there will be no x dependency in the final solution. This should not be surprising since the x derivatives or zero on the boundaries and there is not x dependency in the ic.

Try again

n = 0;

T[x, y, t]== T0
(*c7 Cos[b y]+(h k Tinf+h L q+k q)/(h k)-(q y)/k==T0*)


Shift sides. Multiply each side by cos(b y) and integrate.

c7*Integrate[Cos[b*y]^2, {y, 0, L}, {x, 0, L}] ==
Integrate[(T0 - (h*k*Tinf + h*L*q + k*q)/(h*k) + (q*y)/k)*Cos[b*y], {y, 0, L},
{x, 0, L}];

c7 = c7 /. Solve[%, c7][[1]] // Simplify;

T[x, y, t]
(*(E^(α (-b^2) t) Cos[b y] (4 h q Cos[b L]-4 (b k Sin[b L] (-(h T0)+h Tinf+q)+h q)))/(b h k (2 b L+Sin[2 b L]))+(h k Tinf+h L q+k q)/(h k)-(q y)/k*)


Check pde and bc's so far.

pde // Simplify
(*True*)

dTx = D[T[x, y, t], x]
(*0*)

D[T[x, y, t], y] /. y -> 0
(*-(q/k)*)

FullSimplify[k*D[T[x, y, t], y] + h*T[x, y, t] /. y -> L] == h*Tinf /. bRule
(*True*)


Now form a sum over the b's. We can drop the x dependency. Make it a finite sum so we can compute.

Temp[y_, t_, mm_] := (h*k*Tinf + h*L*q + k*q)/(h*k) - (q*y)/k +
(1/(h*k))*Sum[(E^(α*(-t)*b[[m]]^2)*Cos[y*b[[m]]]*(4*h*q*Cos[L*b[[m]]] -
4*(k*b[[m]]*(-(h*T0) + h*Tinf + q)*Sin[L*b[[m]]] + h*q)))/
(b[[m]]*(2*L*b[[m]] + Sin[2*L*b[[m]]])), {m, 1, mm}]


Plug in some numbers. The equation that determines the b's is a transcendental equation and must be solved numerically.

L = 1;
k = 237;
T0 = 10;
Tinf = 10;
\[Alpha] = 80 10^-6;
q = 10;
h = 25;

fb[b_] = beq[[1]]

Plot[fb[b], {b, 0, 4}]


Calculate the first 200 b's with FindRoot.

m = 200;
b = Table[0, {i, m}];
aa = FindRoot[fb[β] == 0, {β, 0.3},
WorkingPrecision -> 100, AccuracyGoal -> 80];
b[[1]] = β /. aa;
aa = FindRoot[fb[β] == 0, {β, 3.2},
WorkingPrecision -> 100, AccuracyGoal -> 80];
b[[2]] = β /. aa;
For[i = 3, i <= m,
aa = FindRoot[fb[β] == 0, {β, 2 b[[i - 1]] - b[[i - 2]]},
WorkingPrecision -> 100 , AccuracyGoal -> 80];
b[[i]] = β /. aa;
i++]
Clear[m]


Do some plots

Plot3D[Evaluate[Temp[y, t, 200]], {y, 0, L}, {t, 0, 2000},
AxesLabel -> {y, t, Temp}, PlotRange -> {9.9999, 10.02}]


gifs = Table[
Plot[Evaluate[Temp[y, t, 20]], {y, 0, L},
AxesLabel -> {y, Tempurature}, PlotRange -> {10, 10.02},
PlotLabel -> "t = " <> ToString[PaddedForm[N[t], {6, 1}]]], {t, 0,
2000, 20}];


We would get a little more action with a higher value for α.