[Scientific / technical calculation by Python] Numerical solution of one-dimensional unsteady heat conduction equation by Crank-Nicholson method (implicit method) and FTCS method (positive solution method), parabolic partial differential equation

Introduction

Atomic vibration (lattice vibration) energy is transported by the temperature gradient in a substance. This is heat conduction [1]. As a result, the temperature of the substance changes over time and eventually does not change over time. It has reached a steady state. Heat conduction is a common phenomenon in everyday life, and its engineering application (heat transfer engineering) is one of the foundations that support today's life.

The unsteady heat conduction equations that describe temperature fluctuations are classified as parabolic partial differential equations. ** Elementary numerical solutions include the central finite difference method in space and the explicit difference method (FTCS) that uses forward finite difference for time evolution [2]. The FTCS method is easy to understand and implement, but the [numerical stability] of differential equations (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) is not high. ** **

** [Crank-Nicholson Method](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] is one of the implicit methods with excellent numerical stability (unconditional stability). It is necessary to solve simultaneous linear equations to calculate time evolution, and it is more difficult to implement than the FTCS method, but it is useful for solving parabolic partial differential equations with less error for time evolution in addition to stability. It is one of the methods. ** **

** Here, the one-dimensional unsteady heat conduction equation is solved by the Crank-Nicholson method. ** **

Contents

Internal heat generation $ q (x, t) $, thermal equation according to the temperature $ T (x, t) $ of an object with a constant thermal diffusivity $ 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 $ (initial condition) $ T (0, t) = 0 $ (boundary condition) $ T (100, t) = 50 $ (boundary condition)

To

(1) Solve by the crank-Nicholson method. Here, $ \ kappa $ is the thermal conductivity, $ \ rho $ is the density, and $ C_v $ is the equal volume specific heat.

(2) Solve by FTCS method.

Calculation code

(1) Crank-Nicholson method

Faithful implementation. ** I'm using numpy's linalg.solve method to solve a one-dimensional system of equations. ** **


"""
One-dimensional unsteady heat conduction
crank-Nicholson method
"""

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


Nx =100 #Grid score in x direction
Nt =5000#Number of grid points in the t direction
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])  #Function to be sought


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

#boundary condition
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

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

    Amat=np.zeros([Nx-2,Nx-2])  #Generation of coefficient matrix of simultaneous linear equations
    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=Generation of b vector of 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) #Solve simultaneous linear equations
    for i in range(Nx-2):
        uu[i+1,j+1]=uvec[i]
        
#for visualization
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

Plot of the solution.


Next, let's see the animation of the state leading up to the thermal equilibrium state.


%matplotlib nbagg 
from matplotlib.animation import ArtistAnimation #Import methods for creating animations

fig = plt.figure()

anim = [] #A list for storing the data of the para-para diagram drawn for 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) #Animation creation
plt.xlabel('x')
plt.ylabel('t')

fig.show() 

anim.save("t.gif", writer='imagemagick')   #Animation.Save as gif and create a gif animation file.



t.gif

As time goes by, the temperature approaches constant (steady state).

(2) FTCS method

"""
One-dimensional unsteady heat conduction
FTCS method
"""

#%matplotlib nbagg
#matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import ArtistAnimation #Import methods for creating animations

#constant
L = 0.01
D= 4e-6 #Thermal diffusivity
N=100 #Number of steps
del_L= L/N #Space step size
del_t=  0.0001#Time step width
dum = del_t/1000

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

T_low = 0.0
T_mid = 20.0
T_high=50.0

#Illustrated
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

#Main

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

while t < t_end :
    #Temperature calculation
   # 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
    
    #Plot with the selected t
    
    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

A plot of the temperature profile with some time selected.


References

[1] Regarding the derivation of the heat conduction equation, Qiita article: "Derivation of heat conduction equation" Is polite and easy to understand.

[2] Guo Shigeru Yamazaki, ["Introduction to Numerical Solving of Partially Differentiated Equations"](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.


Addendum

  1. Two-dimensional heat conduction equations can also be solved by the Crank-Nicholson method, but the number of times to solve simultaneous linear equations increases. ** Assuming that the number of spatial grids is $ Ns $ and the time grid is $ Nt $, it is necessary to solve approximately $ Nt \ times Ns ^ N $ simultaneous linear equations in order to solve the $ Ns $ dimensional heat conduction equation. As the number of dimensions increases, the calculation cost of the implicit method by the Crank-Nicholson method increases significantly. ** One of the methods to reduce the calculation cost is ** Alternating Direction Implicit method (ADI method) **.

  2. Examples of parabolic PDEs in physics include heat conduction equations and diffusion equations. % 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) and time-dependent Schledinger equation 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) and so on.

Recommended Posts

[Scientific / technical calculation by Python] Numerical solution of one-dimensional unsteady heat conduction equation by Crank-Nicholson method (implicit method) and FTCS method (positive solution method), parabolic partial differential equation
[Scientific / technical calculation by Python] Numerical solution of one-dimensional and two-dimensional wave equations by FTCS method (explicit method), hyperbolic partial differential equations
[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
[Scientific / technical calculation by Python] Numerical solution of one-dimensional harmonic oscillator problem by velocity Verlet method
[Scientific / technical calculation by Python] Numerical solution of second-order ordinary differential equations, initial value problem, numerical calculation
[Scientific / technical calculation by Python] Numerical solution of eigenvalue problem of matrix by power method, numerical linear algebra
[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
[Science / technical calculation by Python] Comparison of convergence speeds of SOR method, Gauss-Seidel method, and Jacobi method for Laplace equation, partial differential equation, boundary value problem
[Python] Numerical calculation of two-dimensional thermal diffusion equation by ADI method (alternate direction implicit method)
[Science / technical calculation by Python] Numerical solution of first-order ordinary differential equations, initial value problem, numerical calculation
[Scientific / technical calculation by Python] Sum calculation, numerical calculation
[Scientific / technical calculation by Python] Fitting by nonlinear function, equation of state, scipy
[Scientific / technical calculation by Python] Drawing animation of parabolic motion with locus, matplotlib
[Scientific / technical calculation by Python] Lagrange interpolation, numerical calculation
[Scientific / technical calculation by Python] Solving one-dimensional Schrodinger equation in stationary state by shooting method (1), well-type potential, quantum mechanics
[Scientific / technical calculation by Python] Solving one-dimensional Schrodinger equation in stationary state by shooting method (2), harmonic oscillator potential, quantum mechanics
[Scientific / technical calculation by Python] List of matrices that appear in Hinpan in numerical linear algebra
[Scientific and technical calculation by Python] Drawing of fractal figures [Sierpinski triangle, Bernsley fern, fractal tree]
Python --Differential equation numerical solution Euler method & central difference method & Runge-Kutta method
[Numerical calculation method, python] Solving ordinary differential equations by Eular method
[Scientific / technical calculation by Python] Basic operation of arrays, numpy
[Scientific / technical calculation by Python] Monte Carlo integration, numerical calculation, numpy
[Scientific / technical calculation by Python] Monte Carlo simulation of thermodynamics of 2D Ising spin system by Metropolis method
Scientific / technical calculation by Python] Drawing and visualization of 3D isosurface and its cross section using mayavi
[Scientific / technical calculation by Python] Numerical integration, trapezoidal / Simpson law, numerical calculation, scipy
[Scientific / technical calculation by Python] Solving simultaneous linear equations, numerical calculation, numpy
[Scientific / technical calculation by Python] 2D random walk (drunk walk problem), numerical calculation
[Science and technology calculation by Python] (Partial) differential, mathematical formula, sympy
Numerical calculation of partial differential equations with singularity (for example, asymptotic behavior analysis of Hardy-Hénon type heat equation)
[Scientific / technical calculation by Python] Calculation of matrix product by @ operator, python3.5 or later, numpy
[Scientific / technical calculation by Python] Solving ordinary differential equations, mathematical formulas, sympy
[Scientific / technical calculation by Python] Drawing of 3D curved surface, surface, wireframe, visualization, matplotlib