Einführung

Ich studiere physikalische Eigenschaften von Halbleitern an der Universität. Wir führen Experimente wie die Bestrahlung von Halbleitermaterialien mit einem Laser und die Messung der Wärmeausdehnung aufgrund der Lichtabsorption mit einem piezoelektrischen Element (Element, das Expansion und Kontraktion in Spannung umwandelt) durch und reproduzieren die experimentellen Ergebnisse durch Berechnung, Wärmeleitfähigkeit usw. Ich versuche, den Wert der physikalischen Eigenschaft von zu berechnen.

Ich möchte die Wärmediffusionsgleichung für diese Berechnung numerisch lösen. Da eine Dimension nicht ausreichte, um die experimentellen Ergebnisse zu reproduzieren, haben wir die zweidimensionale Wärmediffusionsgleichung gelöst.

In diesem Artikel werden wir die zweidimensionale Wärmediffusionsgleichung diskutieren. Es wird durch die ADI-Methode (Alternating Direction Implicit Method) gelöst. Bitte lassen Sie mich wissen, wenn Sie Fehler machen oder ein besseres Programm schreiben.

Die Gleichung, die Sie lösen möchten, lautet wie folgt.

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

Wo $ T $ die Temperatur ist, ist $ t $ die Zeit, $ \ lambda $ die Wärmediffusion und $ x, y $ die Koordinaten.

Ich habe die Wärmediffusionsgleichung in Ableitung der Wärmeleitungsgleichung untersucht.

Darüber hinaus wird die numerische Berechnung der eindimensionalen Wärmediffusionsgleichung durch Fluidsimulation: Implementierung der Diffusionsgleichung untersucht. Hat.

Sie können die Temperatur eines Stabes in 1D, einer Oberfläche in 2D und eines Volumenkörpers in 3D berechnen.

Wenn die Crank Nicholson-Methode, die häufig bei der numerischen Berechnung eindimensionaler Wärmediffusionsgleichungen verwendet wird, auf zwei Dimensionen angewendet wird, wird der Rechenaufwand enorm. Andererseits ermöglicht die Verwendung der ADI-Methode zweidimensionale Berechnungen, während der Rechenaufwand reduziert wird.

Bild der ADI-Methode

Bei der ADI-Methode unter der Annahme, dass alle Temperaturen zu einem Zeitpunkt $ k $ erhalten werden, koordiniert $ y $ bei der Berechnung des nächsten Zeitraums von $ k + 1 $ $ j = 1,2, ... m $, wie in der Abbildung gezeigt. Die Temperatur der gesamten Oberfläche wird berechnet, indem $ i = 1,2, ... n $ in der $ x $ -Richtung für jede von berechnet wird. Um die nächste $ k + 2 $ Zeit zu finden, berechnen Sie $ j = 1,2 ... m $ in $ y $ Richtung für jede der $ x $ Koordinaten $ i = 1,2, ... n $. Auf diese Weise wird die Gesamttemperatur berechnet. Die $ k + 3 $ -Zeit ist in der $ x $ -Richtung, die $ k + 4 $ -Zeit ist in der $ y $ -Richtung und so weiter. Ich denke, dass es eine implizite Methode mit wechselnder Richtung genannt wird, weil es die $ x $ -Richtung und die $ y $ -Richtung abwechselnd berechnet.

(Normalerweise befindet sich der Ursprung unten links und die x-Achse ist horizontal und die y-Achse ist nach oben geschrieben. Da es eine Lösung von y in der Richtung gibt, ist es dieses Bild in mir geworden. Es tut mir leid, dass es schwer zu verstehen ist.)

Berechnung der ADI-Methode

Angenommen, alle $ T_ {i, j} ^ {k} $ zum Zeitpunkt $ k $ wurden berechnet, $ T_ {i, j} ^ {k + 1} zum nächsten Zeitpunkt $ k + 1 $ Zur Berechnung von $ wird die in "Einführung" gezeigte Wärmediffusionsgleichung wie folgt unterschieden. $ \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) $

Dies wird wie folgt in $ k + 1 $ auf der linken Seite und $ k $ auf der rechten Seite umgewandelt. $ 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) $ Setzen Sie das so $ 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) $

Hier, $ 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] $

Wenn Sie $ i = 1,2, ... m $ setzen und dies als Matrix ausdrücken,


\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 $ Man erhält ein lineares Gleichungssystem.

Diese Transformation verwendet die Diricre-Randbedingung (die Bedingung, dass die Grenztemperatur konstant ist). Da $ T_ {0, j} ^ k $ und $ T_ {m + 1, j} ^ k $ die Temperatur der Grenze sind und bekannte Zahlen sind, ist die unbekannte Zahl $ m $.

Das Lösen dieser simultanen Gleichung für jeden Fall von $ j = 1 $ bis $ j = n $ ergibt die Temperatur zum Zeitpunkt $ k + 1 $.

Als nächstes wird angenommen, dass alle $ T_ {i, j} ^ {k} $ zum Zeitpunkt $ k + 1 $ berechnet wurden, $ T_ {i, j} ^ zum nächsten Zeitpunkt $ k + 2 $ Differenzieren Sie wie folgt, um {k + 1} $ zu berechnen. $ \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) $

Auf die gleiche Weise $ 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) $ Setzen Sie das so $ 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) $

Hier, $ 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] $

Eine Matrix kann auf die gleiche Weise gezeichnet werden, wenn die Temperatur der nächsten $ k + 2 $ Zeit aus der $ k + 1 $ Zeit berechnet wird. (Auslassen)

Das Setzen von $ j = 1,2, ... n $ ergibt eine $ n $ elementare simultane lineare Gleichung für $ n $ Unbekannte. Das Lösen der simultanen Gleichungen in jedem Fall von $ i = 1 $ bis $ i = m $ ergibt die Temperatur zum Zeitpunkt $ k + 2 $.

Bei der ADI-Methode können die gleichzeitig zu lösenden Elemente simultaner Gleichungen $ m $ yuan oder $ n $ yuan sein.

Wenn Sie versuchen, auf die gleiche Weise mit der Crank Nicholson-Methode zu berechnen, müssen Sie die simultanen linearen Gleichungen $ m \ times n $ elementar lösen.

Programm

Die Berechnungsbedingungen sind

Ein Objekt mit einer Größe von 10 in der Länge und 10 in der Breite und einer Wärmediffusionsrate von 0,1
Berechnet als $ dx = dy = 0,1 $, $ dt = 0,1 $
Zum Zeitpunkt 0 beträgt die Temperatur des gesamten Objekts 0 Grad
Die Temperatur der vier Seiten des Objekts ist auf 100 Grad festgelegt


# -*- coding: utf-8 -*-
"""
Berechnung der zweidimensionalen Wärmediffusionsgleichung durch implizite Wechselrichtungsmethode
"""

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 #Wie oft muss die zu berechnende Zeit geteilt werden? Es ist besser, die Einerstelle auf 1 zu setzen. Sie bezieht sich auf die Zahl am Ende der for-Anweisung zum Zeitpunkt t
dxnum = 101 #Wie viele Teilungen x und y sollen berechnet werden
dynum = 101

thick = 10 #Größe in x-Richtung
width = 10 #Größe in y-Richtung
time_calc = 500 #Zeit zum Berechnen
beta = 0.1 #In der obigen Formel ist die Lambda-Temperaturdiffusionsrate Lambda=Wärmeleitfähigkeit/(Dichte*spezifische Wärme)

"""Berechnungsvorbereitung"""
#Bereiten Sie eine leere Lösung vor
solution_k1 = np.zeros([dxnum,dynum])
solution_k2 = np.zeros([dxnum,dynum]) #solution_Der Anfangswert von k2 wird zur Anfangsbedingung


#Randbedingung
irr_boundary = 100 #Oberflächenrandbedingungen(0,t)Temperatur in
rear_boundary = 100 #Randzustand auf der gegenüberliegenden Seite des Rückens(x,t)Als Referenztemperatur wurde die Temperatur 0 verwendet
side_0 = 100 #y=Temperatur von 0
side_y = 100 #y=Temperatur an der Grenze gegenüber 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=Bereiten Sie A von b vor"""
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)      
    
#Speichert die CSR-Methode mit spärlicher Matrix
a_matrix_csr1 = csr_matrix(a_matrix_1)
a_matrix_csr2 = csr_matrix(a_matrix_2)


#Bei der ADI-Methode k+1 Mal und k+Da 2 Mal berechnet werden, halbiert sich die Anzahl der for-Anweisungen.
number = Decimal(str(dtnum/2)).quantize(Decimal('0'), rounding=ROUND_HALF_UP)
number = int(number)

solution = np.zeros([dxnum,dynum,number]) #Erstellen Sie eine Matrix, um die Lösung zu ersetzen

#temp_temperature_Array dient zum Erstellen und Hinzufügen von versetzten Spalten
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

#Berechnung der ADI-Methode
for k in range(number): #Über die Zeit t
    for j in range(dynum):#Durch Berechnen von x und Wiederholen der Anzahl von y wird die Zeit k+Berechnen Sie Tij von 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])
        
        #Randbedingung am Anfang und Ende von b
        b_array_1[0] -= irr_boundary
        b_array_1[-1] -= rear_boundary
        
        #Finde eine Lösung
        temperature_solve1 = spla.dsolve.spsolve(a_matrix_csr1, b_array_1)#Lösung für x
        solution_k1[:,j] = temperature_solve1
        
    for i in range(dxnum):#Durch Berechnen von y und Wiederholen der Anzahl von x wird die Zeit k+Berechnen Sie Tij von 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,:])

        #Randbedingung am Anfang und Ende von b
        b_array_2[0] -= side_0
        b_array_2[-1] -= side_y
        
        #Finde eine Lösung
        temperature_solve2 = spla.dsolve.spsolve(a_matrix_csr2, b_array_2)#Lösung für 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()

Wenn Sie es ausführen, dauert es auf meinem Computer ungefähr 50 Sekunden, um das Berechnungsergebnis zu erhalten. Das Berechnungsergebnis ist wie folgt.

Es kann berechnet werden, dass Wärme im Laufe der Zeit ab dem Anfangszustand von 0 Grad allmählich aus der Nähe der Grenze übertragen wird und sich 100 Grad nähert.

Am Ende

Gibt es eine angemessene Schrittweite? Welche Schrittgröße ist erforderlich, wenn die Größe des Objekts in der Größenordnung von Mikro oder Nanometer berechnet wird und die Zeit in der Größenordnung von Millisekunden liegt? Sollten $ dx $ und $ dy $ auch bei der impliziten Methode kleiner als $ dt $ sein? Ich würde es begrüßen, wenn Sie es mir sagen könnten.

Verweise

Guo Shigeru Yamazaki, Einführung in die numerische Lösung partiell differenzierter Gleichungen, Morikita Publishing Co., Ltd. (1993).
http://nalab.mind.meiji.ac.jp/~mk/labo/text/heat-fdm-1.pdf
https://qiita.com/KQTS/items/97daa509991bb9777a4a