[Calcul scientifique / technique par Python] Solution numérique d'une équation de conduction thermique non stationnaire unidimensionnelle par méthode Crank-Nicholson (méthode implicite) et méthode FTCS (méthode de solution positive)

introduction

L'énergie de vibration atomique (vibration de réseau) est transportée par le gradient de température dans le matériau. C'est la conduction thermique [1]. En conséquence, la température de la substance change avec le temps et finalement ne change pas avec le temps. Il a atteint un état stable. La conduction thermique est un phénomène courant dans la vie de tous les jours, et son application d'ingénierie (génie du transfert de chaleur) est l'un des fondements qui soutiennent la vie d'aujourd'hui.

Les équations de conduction thermique non stationnaires qui décrivent les fluctuations de température sont classées comme des équations différentielles partielles paraboliques. ** Les solutions numériques élémentaires incluent la méthode des différences centrales pour l'espace et la méthode explicite (FTCS) qui utilise les différences directes pour l'évolution temporelle [2]. La méthode FTCS est facile à comprendre et à mettre en œuvre, mais la [stabilité numérique] de l'équation différentielle (https://ja.wikipedia.org/wiki/%E6%95%B0%E5%80%A4%E7%9A% 84% E5% AE% 89% E5% AE% 9A% E6% 80% A7 # .E6.95.B0.E5.80.A4.E5.BE.AE.E5.88.86.E6.96.B9.E7 .A8.8B.E5.BC.8F.E3.81.A7.E3.81.AE.E5.AE.89.E5.AE.9A.E6.80.A7) n'est pas élevé. ** **

** [Méthode Crank-Nicholson](https://ja.wikipedia.org/wiki/%E5%B7%AE%E5%88%86%E6%B3%95#.E3.82.AF.E3.83 .A9.E3.83.B3.E3.82.AF.E3.83.BB.E3.83.8B.E3.82.B3.E3.83.AB.E3.82.BD.E3.83.B3. E6.B3.95) [2] est l'une des méthodes implicites avec une excellente stabilité numérique (stabilité inconditionnelle). Il est nécessaire de résoudre des équations linéaires simultanées pour calculer l'évolution temporelle, et c'est plus difficile à mettre en œuvre que la méthode FTCS, mais en plus de la stabilité, l'erreur par rapport à l'évolution temporelle est faible et elle est utile pour résoudre l'équation différentielle partielle parabolique. C'est l'une des méthodes. ** **

** Ici, l'équation de conduction thermique non stationnaire unidimensionnelle est résolue par la méthode Crank-Nicholson. ** **

Contenu

Génération de chaleur interne $ q (x, t) $, équation thermique en fonction de la température $ T (x, t) $ d'un objet à taux de diffusion thermique constant $ D $,

$ \frac{\partial T(x,t)}{\partial t} = D \frac{\partial^2 T(x,t)}{\partial x^2} +q(x,t)$

D = \frac{\kappa} {\rho\  C_V}

$ T (x, 0) = 20 $ (condition initiale) $ T (0, t) = 0 $ (condition aux limites) $ T (100, t) = 50 $ (condition aux limites)

À

(1) Résoudre par la méthode crank-Nicholson. Ici, $ \ kappa $ est la conductivité thermique, $ \ rho $ est la densité et $ C_v $ est la chaleur spécifique à volume égal.

(2) Résoudre par la méthode FTCS.

Code de calcul

(1) Méthode Crank-Nicholson

Mise en œuvre fidèle. ** J'utilise la méthode linalg.solve de numpy pour résoudre une équation simultanée unidimensionnelle. ** **


"""
Conduction thermique non stationnaire unidimensionnelle
manivelle-Méthode Nicholson
"""

%matplotlib notebook
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np


Nx =100 #Score de la grille dans la direction x
Nt =5000#Nombre de points de grille dans la direction t
Lx =0.01
Lt =1.5
delta_x=Lx/Nx
delta_t=Lt/Nt
r=delta_t/(delta_x**2)
print("r=",r)

uu = np.zeros([Nx,Nt])  #Fonction à rechercher


#Conditions initiales
#for i in range(1,Nx-1):
uu[:,0] = 20   #Conditions initiales

#condition limite
for i in range(Nt):
    uu[0,i] = 0  
    uu[-1,i] = 50

p=np.ones([Nx,Nt])
for i in range(Nx):
    p[i,:] =4e-6

#print("stability=",p[0,0]*r)
q=np.zeros([Nx,Nt])

alpha=np.ones([Nx,Nt])
for i in range(Nx):
    alpha[i,:]= r*p[i,:]/2

#Principale
for j in range(Nt-1):

    Amat=np.zeros([Nx-2,Nx-2])  #Génération d'une matrice de coefficients d'équations linéaires simultanées
    for i in range(Nx-2):
        Amat[i,i] = 1/alpha[i,j] +2
        if i >=1 :
            Amat[i-1,i] = -1
        if i <= Nx-4 :
                Amat[i+1,i] = -1

    
    bvec=np.zeros([Nx-2]) # Ax=Génération de b vecteur de b
    for i in range(Nx-2):
        bvec[i] =  uu[i,j]+ (1/alpha[i+1,j] - 2)*uu[i+1,j]+uu[i+2,j]+q[i+1,j]
    bvec[0] += uu[0,j+1]
    bvec[Nx-3] += uu[-1,j+1]
    
    uvec = np.linalg.solve(Amat ,bvec) #Résoudre des équations linéaires simultanées
    for i in range(Nx-2):
        uu[i+1,j+1]=uvec[i]
        
#pour la visualisation
x=list(range(Nx))
y=list(range(Nt))

X, Y = np.meshgrid(x,y)

def functz(u):
    z=u[X,Y]
    return z

Z = functz(uu)
fig = plt.figure()
ax = Axes3D(fig)
ax.plot_wireframe(X,Y,Z, color='r')
ax.set_xlabel('t')
ax.set_ylabel('x')
ax.set_zlabel('T')

plt.show()

スクリーンショット 2017-08-24 19.41.20.png

Tracé de la solution.


Ensuite, regardons l'état d'atteindre l'équilibre thermique avec l'animation.


%matplotlib nbagg 
from matplotlib.animation import ArtistAnimation #Importer des méthodes pour créer des animations

fig = plt.figure()

anim = [] #Une liste pour stocker les données du diagramme para-para dessiné pour l'animation
for i in range(Nt):
    T=list(uu[:,i])
    x=list(range(Nx))
    if i % int(Nt*0.02) ==0: 
        im=plt.plot(x,T, '-', color='red',markersize=10, linewidth = 2, aa=True)
        anim.append(im)

anim = ArtistAnimation(fig, anim) #Création d'animation
plt.xlabel('x')
plt.ylabel('t')

fig.show() 

anim.save("t.gif", writer='imagemagick')   #Animation.Enregistrez-le en tant que gif et créez un fichier d'animation gif.



t.gif

Avec le temps, la température se rapproche de la constante (régime permanent).

(2) Méthode FTCS

"""
Conduction thermique non stationnaire unidimensionnelle
Méthode FTCS
"""

#%matplotlib nbagg
#matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import ArtistAnimation #Importer des méthodes pour créer des animations

#constant
L = 0.01
D= 4e-6 #Taux de diffusion de chaleur
N=100 #Nombre d'étapes
del_L= L/N #Taille de pas d'espace
del_t=  0.0001#Largeur du pas de temps
dum = del_t/1000

print("stability=",D*del_t/(del_L**2))

T_low = 0.0
T_mid = 20.0
T_high=50.0

#Illustré
t1 = 0.01
t2 = 0.1
t3 = 0.4
t4 = 1.0
t5 = 10.0 
t_end = t5 +dum 

#
T = np.empty(N+1)
T[0] = T_high
T[N] = T_low
T[1:N] = T_mid

Tp = np.empty(N+1)
Tp[0] = T_high
Tp[N] = T_low

#Principale

t = 0.0
c =  del_t*D/(del_L**2)

while t < t_end :
    #Calcul de la température
   # for i in range(1,N):
    #   Tp[i] = T[i] + c*(T[i+1]+T[i-1]-2*T[i])
    Tp[1:N] = T[1:N] + c*(T[0:N-1]+T[2:N+1]-2*T[1:N])

    T, Tp = Tp, T
    t += del_t
    
    #Tracer avec le t sélectionné
    
    if np.abs(t-t1) < dum :
        plt.plot(T, label='t = t1')
        
    if np.abs(t-t2) < dum :
        plt.plot(T, label='t = t2')
    if np.abs(t-t3) < dum :
        plt.plot(T, label='t = t3')    
    if np.abs(t-t4) < dum: 
        plt.plot(T, label='t = t4')
    if np.abs(t-t5) < dum :
        plt.plot(T, label='t = t5')

plt.xlabel('x', fontsize=24)
plt.ylabel('T', fontsize=24)
plt.legend(loc='best')


plt.show()

t.png

Un tracé du profil de température sélectionné à plusieurs reprises.


Les références

[1] Concernant la dérivation de l'équation de conduction thermique Article Qiita: "Dérivation de l'équation de conduction thermique" Est poli et facile à comprendre.

[2] Guo Shigeru Yamazaki, ["Introduction à la résolution numérique d'équations partiellement différenciées"](https://www.google.co.id/search?q=%E5%B1%B1%E5%B4%8E+%E5% 81% 8F% E5% BE% AE% E5% 88% 86 & oq =% E5% B1% B1% E5% B4% 8E +% E5% 81% 8F% E5% BE% AE% E5% 88% 86 & aqs = chrome. 69i57j0l5.2548j0j7 & sourceid = chrome & ie = UTF-8), Morikita Publishing Co., Ltd., 1993.


Addenda

  1. Les équations de conduction thermique bidimensionnelles peuvent également être résolues par la méthode Crank-Nicholson, mais le nombre de fois pour résoudre des équations linéaires simultanées augmente. ** En supposant que le nombre de grilles spatiales est $ Ns $ et que la grille temporelle est $ Nt $, il est nécessaire de résoudre approximativement $ Nt \ fois Ns ^ N $ équations linéaires simultanées afin de résoudre l'équation dimensionnelle de conduction thermique $ Ns $. À mesure que le nombre de dimensions augmente, le coût de calcul de la méthode implicite par la méthode Crank-Nicholson augmente considérablement. ** Une des méthodes pour réduire le coût de calcul est ** Méthode implicite de direction alternative (méthode ADI) **.

  2. En plus de l'équation de conduction thermique, les exemples qui apparaissent dans la physique des équations aux dérivées partielles paraboliques sont équations de diffusion. % A1% E6% 95% A3% E6% 96% B9% E7% A8% 8B% E5% BC% 8F & oq =% E6% 8B% A1% E6% 95% A3% E6% 96% B9% E7% A8% 8B% E5% BC% 8F & aqs = chrome..69i57j69i61l2.193j0j7 & sourceid = chrome & ie = UTF-8) et dépendant du temps [équation de Schledinger](https://ja.wikipedia.org/wiki/%E3%82%B7% E3% 83% A5% E3% 83% AC% E3% 83% BC% E3% 83% 87% E3% 82% A3% E3% 83% B3% E3% 82% AC% E3% 83% BC% E6% 96% B9% E7% A8% 8B% E5% BC% 8F) et ainsi de suite.

Recommended Posts

[Calcul scientifique / technique par Python] Solution numérique d'une équation de conduction thermique non stationnaire unidimensionnelle par méthode Crank-Nicholson (méthode implicite) et méthode FTCS (méthode de solution positive)
[Calcul scientifique / technique par Python] Solution numérique d'équations d'ondes unidimensionnelles et bidimensionnelles par méthode FTCS (méthode explicite), équations aux dérivées partielles bi-courbes
[Calcul scientifique / technique par Python] Solution numérique de l'équation de Laplace-Poisson bidimensionnelle pour la position électrostatique par la méthode Jacobi, équation aux dérivées partielles elliptiques, problème des valeurs aux limites
[Calcul scientifique / technique par Python] Solution numérique d'un problème d'oscillateur harmonique unidimensionnel par vitesse Méthode de Berle
[Calcul scientifique / technique par Python] Solution numérique d'une équation différentielle ordinaire du second ordre, problème de valeur initiale, calcul numérique
[Calcul scientifique / technique par Python] Solution numérique du problème des valeurs propres de la matrice par multiplication de puissance, algèbre linéaire numérique
[Calcul scientifique / technique par Python] Résolution de l'équation différentielle ordinaire du second ordre par la méthode Numerov, calcul numérique
[Calcul scientifique / technique par Python] Calcul numérique pour trouver la valeur de la dérivée (différentielle)
[Calcul scientifique / technique par Python] Solution analytique sympa pour résoudre des équations
[Calcul scientifique / technique par Python] Résolution de l'équation de Newton unidimensionnelle par la méthode Runge-Kutta du 4ème ordre
[Calcul scientifique / technique par Python] Comparaison des vitesses de convergence de la méthode SOR, de la méthode Gauss-Seidel et de la méthode Jacobi pour l'équation de Laplace, équation différentielle partielle, problème des valeurs aux limites
[Python] Calcul numérique d'une équation de diffusion thermique bidimensionnelle par méthode ADI (méthode implicite de direction alternative)
[Calcul scientifique / technique par Python] Solution numérique d'équations différentielles ordinaires du premier ordre, problème de valeur initiale, calcul numérique
[Calcul scientifique / technique par Python] Calcul de somme, calcul numérique
[Calcul scientifique / technique par Python] Ajustement par fonction non linéaire, équation d'état, scipy
[Calcul scientifique / technique par Python] Dessin d'animation de mouvement parabolique avec locus, matplotlib
[Calcul scientifique / technique par Python] Interpolation de Lagrange, calcul numérique
[Calcul scientifique / technique par Python] Résolution de l'équation de Schrödinger unidimensionnelle à l'état stationnaire par méthode de tir (1), potentiel de type puits, mécanique quantique
[Calcul scientifique / technique par Python] Résolution d'une équation de Schrödinger unidimensionnelle en régime permanent par méthode de tir (2), potentiel d'oscillateur harmonique, mécanique quantique
[Calcul scientifique / technique par Python] Liste des matrices qui apparaissent dans Hinpan en algèbre linéaire numérique
[Calcul scientifique et technique par Python] Dessin de figures fractales [Triangle de Shelpinsky, fougère de Bernsley, arbre fractal]
Python - Solution numérique d'équation différentielle Méthode d'Euler et méthode de différence centrale et méthode Rungekutta
[Méthode de calcul numérique, python] Résolution d'équations différentielles ordinaires par la méthode Eular
[Calcul scientifique / technique par Python] Fonctionnement de base du tableau, numpy
[Calcul scientifique / technique par Python] Intégration Monte Carlo, calcul numérique, numpy
[Calcul scientifique / technique par Python] Simulation de Monte Carlo par la méthode metropolis de la thermodynamique du système de spin ascendant 2D
Calcul scientifique / technique avec Python] Dessin et visualisation d'isoplans 3D et de leurs vues en coupe à l'aide de mayavi
[Calcul scientifique / technique par Python] Intégration numérique, loi trapézoïdale / Simpson, calcul numérique, scipy
[Calcul scientifique / technique par Python] Résolution d'équations linéaires simultanées, calcul numérique, numpy
[Calcul scientifique / technique par Python] Marche aléatoire 2D (problème de marche ivre), calcul numérique
[Calcul scientifique / technique par Python] différentiel (biaisé), formule mathématique, sympy
Calcul numérique d'équations aux dérivées partielles avec singularité (en utilisant les équations thermiques de type Hardy-Hénon à titre d'exemple)
[Calcul scientifique / technique par Python] Calcul du produit de la matrice par l'opérateur @, python3.5 ou supérieur, numpy
[Calcul scientifique / technique par Python] Résolution d'équations différentielles ordinaires, formules mathématiques, sympy
[Calcul scientifique / technique par Python] Dessin de surface courbe 3D, surface, fil de fer, visualisation, matplotlib