# How to modelize the temperature profile of a 2D plate heated on one edge?

I would like to modelize the temperature profile of a 2D plate, heated on one edge at temp T0=100. The external temp is Tinf=0.The Biot number is Bi=0.1. I've already discretised my plate like that :

I know the expression of the temperature on each point, as a function of the temperature of the previous point (for each family Color):

Parameters:

T0 = 100; Tinf = 0; Bi = 0.1;


Temperature function:

    T[i_, j_] =
If[i == 1, T0,
If[j == 1 && i != 1 && i != 21,
0.5  (0.5 (T[i, j + 1] + T[i, j - 1]) + T[2, j]),
If[j == 6 && i != 1 && i != 21,
1/(2 + Bi) (0.5 (T[6, j + 1] + T[6, j - 1]) + T[5, j] + Bi Tinf),
If[i == 21 && j != 1 && j != 6,
1/(2 (1 + Bi)) (T[i, 20] + 0.5 (T[i + 1, 21] + T[i - 1, 21]) +
Bi inf),
If[i == 21 && j == 1, 1/(2 + Bi) (T[1, 20] - T[2, 21] + Bi Tinf),
If[i == 21 && j == 6,
1/(2 (1 + Bi)) (T[6, 20] + T[5, 21] + Bi Tinf),
0.25 (T[i + 1, j] + T[i - 1, j] + T[i, j - 1] +
T[i, j + 1])]]]]]];


Matrix of temperatures :

M = Table[T[i, j], {i, 21}, {j, 6}]


But because my function is iterative, the matrix M won't show up : too long!

I know I should maybe use the Gauss-Seidel method, but I don't know how to... I would like to have something as follows, but obviously not that temperature profile :

Thank you for your help, I'm a beginner at Mathematica and I really don't know how to start.

• You need to put your conditions into a matrix and solve linear system. Feb 11, 2018 at 14:29
• I gave an answer to a related question that might give you some ideas! mathematica.stackexchange.com/questions/159906/… Feb 11, 2018 at 15:15
• Are you trying to solve heat equation? Feb 11, 2018 at 18:59
• @NavidRajil I think it's about the steady-state solution of the heat equation with Dirichlet boundary conditions. Feb 11, 2018 at 19:33
• @NavidRajil Now that I think of it: The OP probably aims at Dirichlet boundary conditions at the top boundary and Robin boundary conditions on the other boundaries... Feb 11, 2018 at 20:30

That's a good problem to investigate!

What you wrote down so far looks basically like the action of the Laplace operator on a temperature configuration T. You problem can be formulated as a linear system of the form A.x == b, where x and b are vectors of size n m (e.g., n=6, m=21 in your example) and the temperature matrix T is related to x by T = Partition[x,n]. One way to solve this could be to build the matrix A (preferably as a SparseArray) and the vector b (basically, b is zero for interior points and with the prescribed temperatures at the boundary points) and to use LinearSolve afterwards.

A second possibility is to write a function that calculates A.x and to apply iterative methods such as the conjugate gradient method, BiCStab, or GMRES to solve the equation A.x==b. These methods are often called Krylov methods. Lucky as we are, they are already implemented in Mathematica, but they are somewhat hidden. The method to use is SparseArrayKrylovLinearSolve. It takes as options basically all suboptions of "Krylov" in LinearSolve

This is how we can define the right hand side b and the action of the matrix A on a vector x; since the points of the computational domain lie on a regular grid, we merely have to use shift operations on rows and columns of T :

Biot = 0.1;
β = 1.;
α = Biot;
T0 = 100.;
T∞ = 0.;
m = 200/2;
hm = 20./m;
n = 60/2;
hn = 6./n;
b = ConstantArray[0., {m, n}];
b[[-1, All]] = α T∞;
b[[All, 1]] = α T∞;
b[[All, -1]] = α T∞;
b[[1, All]] = T0;
b = Flatten[b];
A = x0 \[Function] Module[{T = Partition[x0, n], y},
y = ConstantArray[0., {m, n}];
(*interior part*)
y[[2 ;; -2, 2 ;; -2]] = Subtract[
T[[2 ;; -2, 2 ;; -2]],
0.25 (T[[1 ;; -3, 2 ;; -2]] + T[[3 ;; -1, 2 ;; -2]] + T[[2 ;; -2, 1 ;; -3]] + T[[2 ;; -2, 3 ;; -1]])
];

(*Dirichlet conditions at top boundary*)
y[[1, 1 ;; -1]] = T[[1, 1 ;; -1]];
(*Robin conditions at left boundary*)
y[[2 ;; -1,1]] = (α + β/hn) T[[2 ;; -1, 1]] + (-β/hn) T[[2 ;; -1, 2]];

(*Robin conditions at right boundary*)
y[[2 ;; -1, -1]] = (α + β/hn) T[[2 ;; -1, -1]] + (-β/hn) T[[2 ;; -1, -2]];

(*Robin conditions at bottom boundary*)
y[[-1, 1 ;; -1]] -= (α + β/hm ) T[[-1,1 ;; -1]] + (-β/hm) T[[-2, 1 ;; -1]];

Flatten[y]
];


The solving step can be started with

T = Partition[
SparseArrayKrylovLinearSolve[A, b, Method -> "GMRES"],
n
];


Finally, we can plot the result:

ArrayPlot[T, ColorFunction -> "DarkRainbow"]


If you want it to look like a heated piece of iron, I would suggest the following color scheme:

ArrayPlot[T, ColorFunction -> "SunsetColors"]


Remark:

If one desires to increase m and n, one definitely has to think about a good preconditioner; otherwise SparseArrayKrylovLinearSolve` will take forever.