Problem with boundary condition 2D heat transfer

Question

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.

Answer 1

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.

Answer 2

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 α.

Answer 3

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["NDSolve`FEM`"]
\[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.