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

Question

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. — Vsevolod A., Feb 11 '18 at 14:29
I gave an answer to a related question that might give you some ideas! https://mathematica.stackexchange.com/questions/159906/simulation-of-diffusion-in-a-grid/159908#159908 — Thies Heidecke, Feb 11 '18 at 15:15
@NavidRajil I think it's about the steady-state solution of the heat equation with Dirichlet boundary conditions. — Henrik Schumacher, Feb 11 '18 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... — Henrik Schumacher, Feb 11 '18 at 20:30

score 3 · Answer 1 · edited Jan 27 '23 at 12:27

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 SparseArray`KrylovLinearSolve. 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[
  SparseArray`KrylovLinearSolve[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 SparseArray`KrylovLinearSolve will take forever.

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

1 Answers1