[Python] Calcul numérique d'une équation de diffusion thermique bidimensionnelle par méthode ADI (méthode implicite de direction alternative)
introduction
J'étudie les propriétés physiques des semi-conducteurs à l'université. Nous menons des expériences telles que l'irradiation de matériaux semi-conducteurs avec un laser et la mesure de la dilatation thermique due à l'absorption de la lumière avec un élément piézoélectrique (élément qui convertit l'expansion et la contraction en tension), et en reproduisant les résultats expérimentaux par calcul, conductivité thermique etc. J'essaye de calculer la valeur de propriété physique de.
Je veux résoudre numériquement l'équation de diffusion thermique pour ce calcul. Unidimensionnel n'était pas suffisant pour reproduire les résultats expérimentaux, nous avons donc résolu l'équation de diffusion thermique bidimensionnelle.
Dans cet article, nous discuterons de l'équation de diffusion thermique bidimensionnelle. Il est résolu par la méthode ADI (méthode implicite de direction alternative). S'il vous plaît laissez-moi savoir si vous faites des erreurs ou écrivez un meilleur programme.
L'équation que vous souhaitez résoudre est la suivante.
\frac{\partial T}{\partial t}=\lambda \left(\frac{\partial^2 T}{\partial x^2}+\frac{\partial^2 T}{\partial y^2}\right)
Où $ T $ est la température, $ t $ le temps, $ \ lambda $ est la diffusion thermique et $ x, y $ sont les coordonnées.
J'ai étudié l'équation de diffusion de chaleur dans Dérivation de l'équation de conduction thermique.
De plus, le calcul numérique de l'équation de diffusion thermique unidimensionnelle est étudié par Simulation fluide: Implémentation de l'équation de diffusion. Fait.
Vous pouvez calculer la température d'une tige en 1D, d'une surface en 2D et d'un solide en 3D.
Si la méthode Crank Nicholson, qui est souvent utilisée dans le calcul numérique des équations de diffusion thermique unidimensionnelles, est appliquée à deux dimensions, la quantité de calcul devient énorme. D'autre part, l'utilisation de la méthode ADI permet des calculs bidimensionnels tout en réduisant la quantité de calcul.
Image de la méthode ADI
Dans la méthode ADI, en supposant que toutes les températures à la fois $ k $ sont obtenues, lorsque le temps $ k + 1 $ suivant est calculé, $ y $ coordonnées $ j = 1,2, ... m $ comme indiqué sur la figure. La température de toute la surface est calculée en calculant $ i = 1,2, ... n $ dans la direction $ x $ pour chacun de. Pour trouver le temps $ k + 2 $ suivant, calculez $ j = 1,2 ... m $ dans la direction $ y $ pour chacune des coordonnées $ x $ $ i = 1,2, ... n $. Ce faisant, la température globale est calculée. Le temps $ k + 3 $ est dans la direction $ x $, le temps $ k + 4 $ est dans la direction $ y $, et ainsi de suite. Je pense que cela s'appelle une méthode implicite de direction alternée car elle calcule alternativement la direction $ x $ et la direction $ y $.
(Habituellement, l'origine est en bas à gauche, et l'axe des x est écrit horizontalement et l'axe des y est écrit vers le haut. Puisqu'il y a une solution de y dans la direction, c'est devenu cette image en moi. Je suis désolé, c'est difficile à comprendre.)
Calcul de la méthode ADI
En supposant que tout $ T_ {i, j} ^ {k} $ au temps $ k $ a été calculé, $ T_ {i, j} ^ {k + 1} au temps $ k + 1 $ suivant Pour calculer
Ceci est transformé en $ k + 1 $ sur le côté gauche et $ k $ sur le côté droit comme suit.
ici,
Si vous définissez $ i = 1,2, ... m $ et l'exprimez sous forme de matrice,
\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)
$ m $ On obtient un système d'équations linéaires.
Cette transformation utilise la condition aux limites de Diricre (la condition selon laquelle la température aux limites est constante). Puisque $ T_ {0, j} ^ k $ et $ T_ {m + 1, j} ^ k $ sont la température de la frontière et sont des nombres connus, le nombre inconnu est $ m $.
La résolution de cette équation simultanée pour chaque cas de $ j = 1 $ à $ j = n $ donne la température au temps $ k + 1 $.
Ensuite, en supposant que tout $ T_ {i, j} ^ {k} $ au temps $ k + 1 $ a été calculé, $ T_ {i, j} ^ au temps $ k + 2 $ suivant Différenciez comme suit pour calculer {k + 1}
De la même manière
ici,
Une matrice peut être dessinée de la même manière lors du calcul de la température du temps $ k + 2 $ suivant à partir du temps $ k + 1 $. (Omettre)
Mettre $ j = 1,2, ... n $ donne une équation linéaire simultanée élémentaire $ n $ pour $ n $ inconnus. La résolution des équations simultanées dans chaque cas de $ i = 1 $ à $ i = m $ donne la température au temps $ k + 2 $.
Dans la méthode ADI, les éléments d'équations simultanées à résoudre en une seule fois peuvent être $ m $ yuan ou $ n $ yuan.
Si vous essayez de calculer de la même manière avec la méthode Crank Nicholson, vous devez résoudre les équations linéaires simultanées élémentaires $ m \ times n $.
programme
Les conditions de calcul sont
- Un objet d'une taille de 10 en longueur et 10 en largeur et un taux de diffusion thermique de 0,1
- Calculé comme $ dx = dy = 0,1 $, $ dt = 0,1 $
- Au temps 0, la température de l'objet entier est de 0 degré
- La température des quatre côtés de l'objet est fixée à 100 degrés
# -*- coding: utf-8 -*-
"""
Calcul de l'équation de diffusion thermique bidimensionnelle par méthode implicite de direction alternée
"""
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 #Combien de fois diviser le temps pour calculer Il est préférable de mettre le chiffre des unités à 1 Il est lié au nombre à la fin de l'instruction for au temps t
dxnum = 101 #Combien de divisions x et y doivent être calculées
dynum = 101
thick = 10 #Taille dans la direction x
width = 10 #Taille dans la direction y
time_calc = 500 #Temps de calculer
beta = 0.1 #Dans la formule ci-dessus, le taux de diffusion de la température lambda lambda=Conductivité thermique/(densité*chaleur spécifique)
"""Préparation des calculs"""
#Préparez une solution vide
solution_k1 = np.zeros([dxnum,dynum])
solution_k2 = np.zeros([dxnum,dynum]) #solution_La valeur initiale de k2 devient la condition initiale
#condition limite
irr_boundary = 100 #Conditions aux limites de surface(0,t)Température en
rear_boundary = 100 #Condition aux limites sur le côté opposé du dos(x,t)La température 0 a été utilisée comme température de référence
side_0 = 100 #y=Température de 0
side_y = 100 #y=Température à la limite opposée à 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=Préparer A de 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)
#Stocke la méthode CSR de matrice clairsemée
a_matrix_csr1 = csr_matrix(a_matrix_1)
a_matrix_csr2 = csr_matrix(a_matrix_2)
#Dans la méthode ADI k+1 fois et k+Étant donné que 2 fois sont calculés, le nombre d'instructions for est divisé par deux.
number = Decimal(str(dtnum/2)).quantize(Decimal('0'), rounding=ROUND_HALF_UP)
number = int(number)
solution = np.zeros([dxnum,dynum,number]) #Créer une matrice pour remplacer la solution
#temp_temperature_array sert à créer des colonnes décalées et à les ajouter
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
#Calcul de la méthode ADI
for k in range(number): #À propos du temps t
for j in range(dynum):#En calculant x et en répétant le nombre de y, le temps k+Calculer Tij de 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])
#Condition aux limites au début et à la fin de b
b_array_1[0] -= irr_boundary
b_array_1[-1] -= rear_boundary
#Trouver une solution
temperature_solve1 = spla.dsolve.spsolve(a_matrix_csr1, b_array_1)#Solution pour x
solution_k1[:,j] = temperature_solve1
for i in range(dxnum):#En calculant y et en répétant le nombre de x, le temps k+Calculer Tij de 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,:])
#Condition aux limites au début et à la fin de b
b_array_2[0] -= side_0
b_array_2[-1] -= side_y
#Trouver une solution
temperature_solve2 = spla.dsolve.spsolve(a_matrix_csr2, b_array_2)#Solution pour 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()
Si vous l'exécutez, il faudra environ 50 secondes sur mon ordinateur pour obtenir le résultat du calcul. Le résultat du calcul est le suivant.
Il est possible de calculer comment la chaleur est progressivement transférée à partir du voisinage de la frontière avec le passage du temps à partir de la condition initiale de 0 degré et s'approche de 100 degrés.
À la fin
Y a-t-il une taille de pas appropriée? Quelle est la taille de pas nécessaire lorsque la taille de l'objet est calculée de l'ordre du micro ou du nanomètre et que le temps est de l'ordre de la milliseconde? Est-ce que $ dx $ et $ dy $ devraient être plus petits que $ dt $ même dans la méthode implicite? Je vous serais reconnaissant si vous pouviez me le dire.
Les références
- Guo Shigeru Yamazaki, Introduction à la résolution numérique d'équations partiellement différentielles, Morikita Publishing Co., Ltd. (1993).
- http://nalab.mind.meiji.ac.jp/~mk/labo/text/heat-fdm-1.pdf
- https://qiita.com/KQTS/items/97daa509991bb9777a4a
Recommended Posts