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

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)

(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()

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.

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()

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

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

[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