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.

(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

An object with a size of 10 in length and 10 in width and a thermal diffusivity of 0.1
Calculated as $ dx = dy = 0.1 $, $ dt = 0.1 $
At time 0, the temperature of the entire object is 0 degrees
The temperature of the four sides of the object is fixed at 100 degrees


# -*- 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.

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

Koji Yamazaki, Introduction to Numerical Solving of Partially Differentiated Equations, Morikita Publishing Co., Ltd. (1993).
http://nalab.mind.meiji.ac.jp/~mk/labo/text/heat-fdm-1.pdf
https://qiita.com/KQTS/items/97daa509991bb9777a4a