[Python] Numerical calculation of two-dimensional thermal diffusion equation by ADI method (alternate direction implicit method)

Introduction

I am studying semiconductor physical properties at university. We are conducting experiments such as irradiating semiconductor materials with a laser and measuring thermal expansion due to light absorption with a piezoelectric element (element that converts expansion and contraction into voltage), and by reproducing the experimental results by calculation, thermal conductivity etc. I am trying to calculate the physical property value of.

I want to solve the thermal diffusion equation numerically for that calculation. One-dimensional was not enough to reproduce the experimental results, so we solved the two-dimensional thermal diffusion equation.

In this article, we will discuss the two-dimensional thermal diffusion equation. It is solved by the ADI method (Alternating direction implicit method). Please let me know if you make any mistakes or write a better program.

The equation you want to solve is as follows.

\frac{\partial T}{\partial t}=\lambda \left(\frac{\partial^2 T}{\partial x^2}+\frac{\partial^2 T}{\partial y^2}\right)

Where $ T $ is temperature, $ t $ is time, $ \ lambda $ is thermal diffusivity, and $ x, y $ are coordinates.

I studied the heat diffusion equation in Derivation of heat conduction equation.

In addition, the numerical calculation of the one-dimensional thermal diffusion equation is studied by Fluid simulation: Implementing the diffusion equation. Did.

You can calculate the temperature of a rod in 1D, a surface in 2D, and a solid in 3D.

Applying the Crank Nicholson method, which is often used in the numerical calculation of one-dimensional thermal diffusion equations, to two dimensions will increase the amount of calculation. On the other hand, the ADI method enables two-dimensional calculation while reducing the amount of calculation.

Image of ADI method

Assuming that all the temperatures at one time $ k $ are calculated by the ADI method, when the next $ k + 1 $ time is calculated, $ y $ coordinates $ j = 1,2, ... m $ as shown in the figure. The temperature of the entire surface is calculated by calculating $ i = 1,2, ... n $ in the $ x $ direction for each of. To find the next $ k + 2 $ time, calculate $ j = 1,2 ... m $ in the $ y $ direction at each of the $ x $ coordinates $ i = 1,2, ... n $. By doing so, the overall temperature is calculated. The $ k + 3 $ time is in the $ x $ direction, the $ k + 4 $ time is in the $ y $ direction, and so on. I think that it is called an alternating direction implicit method because it calculates the $ x $ direction and the $ y $ direction alternately.

ADI_image_r2.png

(Usually the origin is in the lower left, x-axis in the horizontal direction, y-axis in the upward direction, but when you open the solution ndarray in Spyder's variable explorer, there is 0 in the upper left, x in the downward direction, right Since there is a solution of y in the direction, it has become that image in me. I'm sorry it's hard to understand.)

ADI method calculation

Assuming that all $ T_ {i, j} ^ {k} $ at the $ k $ time has been calculated, $ T_ {i, j} ^ {k + 1} at the next $ k + 1 $ time To calculate , the thermal diffusion equation shown in "Introduction" is differentiated as follows. $ \frac{T_{i,j}^{k+1} - T_{i,j}^k}{\Delta t} = \lambda \left(\frac{T_{i-1,j}^{k+1} - 2 T_{i,j}^{k+1} + T_{i+1,j}^{k+1} }{\Delta x^2} + \frac{T_{i,j-1}^k - 2 T_{i,j}^k + T_{i,j+1}^k }{\Delta y^2}\right) $$

When this is transformed into $ k + 1 $ on the left side and $ k $ on the right side, it becomes as follows. $ T_{i-1,j}^{k+1}-\left(2+\frac{\Delta x^2}{\lambda \Delta t}\right)T_{i,j}^{k+1}+T_{i+1,j}^{k+1} = \left[2 \left( \frac{\Delta x^2}{\Delta y^2}\right)-\frac{\Delta x^2}{\lambda \Delta t}\right]T_{i,j}^k-\frac{\Delta x^2}{\Delta y^2} \left( T_{i,j-1}^k + T_{i,j+1}^k \right) $ Put this like this $ T_{i-1,j}^{k+1}-dT_{i,j}^{k+1}+T_{i+1,j}^{k+1} = BT_{i,j}^k-\frac{\Delta x^2}{\Delta y^2} \left( T_{i,j-1}^k + T_{i,j+1}^k \right) $

here, $ d = -\left(2+\frac{\Delta x^2}{\lambda \Delta t}\right), B=\left[2 \left( \frac{\Delta x^2}{\Delta y^2}\right)-\frac{\Delta x^2}{\lambda \Delta t}\right] $

If you set $ i = 1,2, ... m $ and express this as a matrix,


\left(
\begin{array}{cccc}
      d & 1 &  0 & \ldots  & \ldots & 0 \\
      1 & d & 1 & 0 & \ldots & 0 \\
      0  &1 & d & 1 & 0 & \ldots  \\
      \vdots & \ddots & \ddots & \ddots & \ddots & \ddots\\
      0  & 0 & \ddots & 1 & d & 1 \\
      0  & 0 & \ldots & 0 & 1 & d
    \end{array}
  \right)
  \left(
    \begin{array}{c}
      T_{1,j}^{k+1}  \\
      T_{2,j}^{k+1}  \\
      T_{3,j}^{k+1}    \\
      \vdots \\
      T_{m-1,j}^{k+1} \\
      T_{m,j}^{k+1} 
\end{array}
\right)

= \left( \begin{array}{c}
    B T_{1,j}^{k} -\frac{\Delta x^2}{\Delta y^2} \left( T_{1,j-1}^k + T_{1,j+1}^k \right)- T_{0,j}^k \\
      B T_{2,j}^{k} -\frac{\Delta x^2}{\Delta y^2} \left( T_{2,j-1}^k + T_{2,j+1}^k \right)  \\
      B T_{3,j}^{k} -\frac{\Delta x^2}{\Delta y^2} \left( T_{3,j-1}^k + T_{3,j+1}^k \right)  \\
      \vdots \\
      B T_{m-1,j}^{k} -\frac{\Delta x^2}{\Delta y^2} \left( T_{m-1,j-1}^k + T_{m-1,j+1}^k \right)  \\
      B T_{m,j}^{k} -\frac{\Delta x^2}{\Delta y^2} \left( T_{m,j-1}^k + T_{m,j+1}^k \right)-T_{m+1,j}^k
    \end{array} \right)

The $ m $ system of linear equations is obtained.

This transformation uses the Dirichlet boundary condition (the condition that the boundary temperature is constant). Since $ T_ {0, j} ^ k $ and $ T_ {m + 1, j} ^ k $ are the temperature of the boundary and are known numbers, the unknown number is $ m $.

Solving this system of equations for each case from $ j = 1 $ to $ j = n $ gives the temperature at the $ k + 1 $ time.

Next, assuming that all $ T_ {i, j} ^ {k} $ at the $ k + 1 $ time has been calculated, $ T_ {i, j} ^ at the next $ k + 2 $ time Differentiate as follows to calculate {k + 1} . $ \frac{T_{i,j}^{k+2} - T_{i,j}^{k+1}}{\Delta t} = \lambda \left(\frac{T_{i-1,j}^{k+1} - 2 T_{i,j}^{k+1} + T_{i+1,j}^{k+1} }{\Delta x^2} + \frac{T_{i,j-1}^{k+2} - 2 T_{i,j}^{k+2} + T_{i,j+1}^{k+2} }{\Delta y^2}\right) $$

In the same way $ T_{i,j-1}^{k+2}-\left(2+\frac{\Delta y^2}{\lambda \Delta t}\right)T_{i,j}^{k+2}+T_{i,j+1}^{k+2} = \left[2 \left( \frac{\Delta y^2}{\Delta x^2}\right)-\frac{\Delta y^2}{\lambda \Delta t}\right]T_{i,j}^{k+1}-\frac{\Delta y^2}{\Delta x^2} \left( T_{i-1,j}^{k+1} + T_{i+1,j}^{k+1} \right) $ Put this like this $ T_{i,j-1}^{k+2}-eT_{i,j}^{k+2}+T_{i,j+1}^{k+2} = CT_{i,j}^{k+1}-\frac{\Delta y^2}{\Delta x^2} \left( T_{i-1,j}^{k+1} + T_{i+1,j}^{k+1} \right) $

here, $ e=-\left(2+\frac{\Delta y^2}{\lambda \Delta t}\right), C=\left[2 \left( \frac{\Delta y^2}{\Delta x^2}\right)-\frac{\Delta y^2}{\lambda \Delta t}\right] $

When calculating the temperature of the next $ k + 2 $ time from the $ k + 1 $ time, the matrix can be drawn in the same way. (Omit)

Setting $ j = 1,2, ... n $ gives a system of linear equations of $ n $ for $ n $ unknowns. Solving the simultaneous equations in each case from $ i = 1 $ to $ i = m $ gives the temperature at the $ k + 2 $ time.

In the ADI method, the elements of simultaneous equations to be solved at one time can be $ m $ yuan or $ n $ yuan.

If we try to calculate in the same way by the Crank Nicholson method, we have to solve the $ m \ times n $ system of linear equations.

program

The calculation conditions are


# -*- coding: utf-8 -*-
"""
Calculation of two-dimensional thermal diffusion equation by alternating directional implicit method
"""

import numpy as np
import numpy.matlib
from decimal import Decimal, ROUND_HALF_UP
import scipy.sparse.linalg as spla
from scipy.sparse import csr_matrix
import matplotlib.pyplot as plt
import seaborn as sns
import time


dtnum = 5001 #How many times should the time be divided and calculated? It is better to set the ones digit to 1. It is related to the number at the end of the for statement at time t.
dxnum = 101 #How many divisions x and y should be calculated
dynum = 101

thick = 10 #Size in x direction
width = 10 #Size in y direction
time_calc = 500 #Time to calculate
beta = 0.1 #In the above formula, Lambda thermal diffusivity lambda=Thermal conductivity/(density*specific heat)

"""Calculation preparation"""
#Prepare an empty solution
solution_k1 = np.zeros([dxnum,dynum])
solution_k2 = np.zeros([dxnum,dynum]) #solution_The initial value of k2 becomes the initial condition


#boundary condition
irr_boundary = 100 #Surface boundary conditions(0,t)Temperature in
rear_boundary = 100 #Boundary conditions on the opposite side of the back(x,t)Temperature 0 was used as the reference temperature
side_0 = 100 #y=Temperature of 0
side_y = 100 #y=Boundary temperature opposite to 0

dt = time_calc / (dtnum - 1)
dx = thick / (dxnum - 1)
dy = width / (dynum - 1)

d = -(2+(dx**2)/(beta*dt))
B = (2*((dx/dy)**2)-((dx**2)/(beta*dt)))

e = -(2+(dy**2)/(beta*dt))
C = (2*((dy/dx)**2)-((dy**2)/(beta*dt)))


"""Ax=Prepare A of b"""
a_matrix_1 = np.identity(dxnum) * d \
            + np.eye(dxnum, k=1) \
            + np.eye(dxnum, k=-1)
            
a_matrix_2 = np.identity(dynum) * e \
            + np.eye(dynum, k=1) \
            + np.eye(dynum, k=-1)      
    
#Store sparse matrix csr method
a_matrix_csr1 = csr_matrix(a_matrix_1)
a_matrix_csr2 = csr_matrix(a_matrix_2)


#In the ADI method k+1 time and k+Since 2 times are calculated, the number of for statements is halved.
number = Decimal(str(dtnum/2)).quantize(Decimal('0'), rounding=ROUND_HALF_UP)
number = int(number)

solution = np.zeros([dxnum,dynum,number]) #Create a matrix to substitute the solution

#temp_temperature_array is for creating staggered columns and adding them
temp_temperature_array1 = np.zeros([dxnum,dynum+2])
temp_temperature_array1[:,0] = side_0
temp_temperature_array1[:,-1] = side_y

temp_temperature_array2 = np.zeros([dxnum+2,dynum])
temp_temperature_array2[0,:] = irr_boundary
temp_temperature_array2[-1,:] = rear_boundary

#ADI method calculation
for k in range(number): #About time t
    for j in range(dynum):#By calculating x and repeating the number of y, the time k+Calculate Tij of 1
        temp_temperature_array1[:,1:dynum+1] = solution_k2
        b_array_1 = B * temp_temperature_array1[:,j+1] \
            - ((dx/dy)**2) * (temp_temperature_array1[:,j] + temp_temperature_array1[:,j+2])
        
        #Boundary conditions at the beginning and end of b
        b_array_1[0] -= irr_boundary
        b_array_1[-1] -= rear_boundary
        
        #Find a solution
        temperature_solve1 = spla.dsolve.spsolve(a_matrix_csr1, b_array_1)#Solution for x
        solution_k1[:,j] = temperature_solve1
        
    for i in range(dxnum):#By calculating y and repeating the number of x, the time k+Calculate Tij of 2
        temp_temperature_array2[1:dxnum+1,:] = solution_k1
        b_array_2 = C * temp_temperature_array2[i+1,:] \
            - ((dy/dx)**2) * (temp_temperature_array2[i,:] + temp_temperature_array2[i+2,:])

        #Boundary conditions at the beginning and end of b
        b_array_2[0] -= side_0
        b_array_2[-1] -= side_y
        
        #Find a solution
        temperature_solve2 = spla.dsolve.spsolve(a_matrix_csr2, b_array_2)#Solution for y
        solution_k2[i,:] = temperature_solve2
    solution[:,:,k] = solution_k2


ax = sns.heatmap(solution[:,:,10], linewidth=0, vmin=0, vmax=100)
plt.show()

ax = sns.heatmap(solution[:,:,100], linewidth=0, vmin=0, vmax=100)
plt.show()

ax = sns.heatmap(solution[:,:,500], linewidth=0, vmin=0, vmax=100)
plt.show()

ax = sns.heatmap(solution[:,:,2000], linewidth=0, vmin=0, vmax=100)
plt.show()

If you execute it, it will take about 50 seconds on my computer to get the calculation result. The calculation result is as follows.

ADI_results.png

It is possible to calculate how heat is gradually transferred from the vicinity of the boundary and approaches 100 degrees with the passage of time from the initial condition of 0 degrees.

At the end

Is there an appropriate step size? What is the step size required when the size of an object is calculated on the order of micro or nanometer and the time is on the order of milliseconds? Should $ dx $ and $ dy $ be smaller than $ dt $ even in the implicit method? I would appreciate it if you could tell me.

References

Recommended Posts

[Python] Numerical calculation of two-dimensional thermal diffusion equation by ADI method (alternate direction implicit method)
[Scientific / technical calculation by Python] Numerical solution of two-dimensional Laplace-Poisson equation by Jacobi method for electrostatic potential, elliptic partial differential equation, boundary value problem
Numerical calculation of compressible fluid by finite volume method
[Scientific / technical calculation by Python] Numerical solution of one-dimensional and two-dimensional wave equations by FTCS method (explicit method), hyperbolic partial differential equations
[Numerical calculation method, python] Solving ordinary differential equations by Eular method
[Scientific / technical calculation by Python] Numerical solution of eigenvalue problem of matrix by power method, numerical linear algebra
[Scientific / technical calculation by Python] Fitting by nonlinear function, equation of state, scipy
[Scientific / technical calculation by Python] Sum calculation, numerical calculation
[Scientific / technical calculation by Python] Solving second-order ordinary differential equations by Numerov method, numerical calculation
[Scientific / technical calculation by Python] Numerical calculation to find the value of derivative (differential)
[Scientific / technical calculation by Python] Analytical solution to find the solution of equation sympy
[Scientific / technical calculation by Python] Solving one-dimensional Newton equation by the 4th-order Runge-Kutta method
[Scientific / technical calculation by Python] Lagrange interpolation, numerical calculation
Calculation of homography matrix by least squares method (DLT method)
Solving one-dimensional wave equation by finite difference method (python)
[Science / technical calculation by Python] Numerical solution of first-order ordinary differential equations, initial value problem, numerical calculation
[Scientific / technical calculation by Python] Numerical solution of second-order ordinary differential equations, initial value problem, numerical calculation
[Scientific / technical calculation by Python] List of matrices that appear in Hinpan in numerical linear algebra