数值传热学——导热问题

数值传热学——导热问题

引用

CSDN

1.

https://blog.csdn.net/marhs/article/details/141965900

数值传热学是研究热能传递规律及其数值计算方法的一门学科，广泛应用于工程设计和科学研究中。本文将介绍如何使用有限差分法求解二维导热问题，包括显式格式和隐式格式的实现，并通过Python代码进行仿真计算。

一、控制方程

在笛卡尔坐标系下，二维导热问题的控制方程为：

$$
\frac{\partial}{\partial x}{(\lambda \frac{\partial T}{\partial x})} + \frac{\partial}{\partial y}{(\lambda \frac{\partial T}{\partial y})}= \rho c\frac{\partial T}{\partial \tau} \tag{1}
$$

其中，$\lambda$ 是导热系数，$\rho$ 是密度，$c$ 是比热容，$T$ 是温度，$\tau$ 是时间。

二、几何结构及边界条件

考虑一个边长为1的正方形几何体，如图1所示：

图1 几何结构

该几何体的初始条件为均匀的15℃。热物性参数如下：

边界条件均为第一类边界条件：

三、网格的划分

采用均匀正交网格，不再赘述。

四、方程的离散

空间步长采用中心差分。

1、显式格式

时间步长采用向前差分，为显式格式，时间步向前推进求解较为简单：

$$
T^+ = T + \left( \frac{T_{x^+} - 2T_x +T_{x^-}}{{\Delta x}^2} + \frac{T_{y^+} - 2T_y +T_{y^-}}{{\Delta y}^2} \right) \frac{\lambda}{\rho c} \Delta \tau \tag{2}
$$

采用Python实现整个算法：

import numpy as np
import matplotlib.pyplot as plt
import math

def cubic(t=3600*7.,dt=108,n=100,Lx=1.,Ly=1.):
    Ti=15
    Tb=50
    dx=Lx/n
    dy=Ly/n
    x=np.arange(0,1,dx)
    y=np.arange(0,1,dy)
    x,y=np.meshgrid(x,y)
    nt=math.floor(t/dt)
    dt=t/nt
    T=np.ones([n,n])*15
    T[:,0]=np.ones([n])*50
    Tn=T.copy()
    lamb=0.6
    rho=2600.
    c=1000.
    for it in range(0,nt):
        T[1:-1,1:-1] = Tn[1:-1,1:-1] + (
            (Tn[1:-1,:-2]-2*Tn[1:-1,1:-1]+Tn[1:-1,2:])/dx/dx+
            (Tn[:-2,1:-1]-2*Tn[1:-1,1:-1]+Tn[2:,1:-1])/dy/dy
        )*lamb/rho/c*dt
        Tn=T.copy()
    
    fig, ax = plt.subplots()
    ctf=ax.contourf(x,y,T,10,cmap='jet')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_aspect('equal')
    ax.set_xlim(0,1)
    ax.set_ylim(0,1)
    fig.colorbar(ctf)
    plt.show()

if __name__=='__main__':
    cubic()

计算结果如图2所示：


收敛条件的说明：

对于显示格式，当时间步长过大时，会存在不收敛的情况，该类问题的收敛条件为最大$Fo \leq 1$，即：

$$
max \left{ \frac{4 \alpha \Delta \tau}{(\Delta x)^2}, \frac{4 \alpha \Delta \tau}{(\Delta y)^2} \right} \leq 1 \tag{3}
$$

对于本问题：

$$
\Delta \tau \leq (\Delta x)^2 \rho c / \lambda /4 = 0.01^2 \times 2600 \times 1000/0.6 / 4 \approx 108
$$

因此，对于该问题来说，时间步长不可高于108。

2、隐式格式

时间步长采用向后差分，为隐式格式，需要迭代求解每个时间步：

$$
T^+ = T + \left( \frac{T_{x^+}^+ - 2T_x^+ +T_{x^-}^+}{{\Delta x}^2} + \frac{T_{y^+}^+ - 2T_y^+ +T_{y^-}^+}{{\Delta y}^2} \right) \frac{\lambda}{\rho c} \Delta \tau \tag{4}
$$

这里$T^+$为未来时间步的值，$T$为上一时间步。由此可知，中心点的值依赖于四面节点的值，常用追赶法求解，化简上式，得到迭代公式如下：

$$
T^+ = \frac{ T + \left( \frac{T_{x^+}^+ + T_{x^-}^+}{{\Delta x}^2} + \frac{T_{y^+}^+ + T_{y^-}^+}{{\Delta y}^2} \right) \frac{\lambda}{\rho c} \Delta \tau }{ 1+ \left( \frac{2}{{\Delta x}^2}+ \frac{2}{{\Delta y}^2} \right) \frac{\lambda}{\rho c} \Delta \tau } \tag{5}
$$

从而：

$$
T^+ = A T + BT_{x^+}^+ + CT_{x^-}^+ + DT_{y^+}^+ + ET_{y^-}^+ \tag{6}
$$

其中：

$$
A=\frac{1}{ 1+ \left( \frac{2}{{\Delta x}^2}+ \frac{2}{{\Delta y}^2} \right) \frac{\lambda}{\rho c} \Delta \tau } \
B=AFo_x \
C=B \
D=AFo_y \
E=D \
Fo_*=\frac{\lambda \Delta \tau}{\rho c {\Delta *}^2}
$$

求解流程：

• 设置几何参数
• 网格划分
• 初始化温度场
• 设置边界条件
• 开始一个时间步
• 基于式（6）对所有节点进行迭代
• 对比上一个迭代步的值与当前迭代步的值是否相同，计算误差，相对误差小于$1 \times 10^{-6}$，认为收敛
• 将当前迭代值作为该时间步的温度场开始下一轮迭代

采用Python实现整个算法，由于计算中采用了大量的迭代，向量化可以大幅度提高计算速度：

import numpy as np
import matplotlib.pyplot as plt
import math

def cubic(t=3600*7.,dt=1200,n=100,Lx=1.,Ly=1.):
    Ti=15
    Tb=50
    dx=Lx/n
    dy=Ly/n
    x=np.arange(0,1,dx)
    y=np.arange(0,1,dy)
    x,y=np.meshgrid(x,y)
    nt=math.floor(t/dt)
    dt=t/nt
    T=np.ones([n,n])*Ti
    T[:,0]=np.ones([n])*Tb
    Tn=T.copy()
    To=T.copy()
    lamb=0.6
    rho=2600.
    c=1000.
    A=1./(1.+(2./dx/dx+2./dy/dy)*lamb*dt/rho/c)
    Fox=lamb*dt/rho/c/(dx*dx)
    Foy=lamb*dt/rho/c/(dy*dy)
    B=A*Fox
    C=B
    D=A*Foy
    E=D
    for it in range(0,nt):
        for itr in range(0,100):
            To=Tn.copy()
            Tn[1:-1,1:-1]=T[1:-1,1:-1]*A+To[:-2,1:-1]*B+To[2:,1:-1]*C+To[1:-1,:-2]*D+To[1:-1,2:]*E
            #for i in range(1,n-1):                
                #for j in range(1,n-1):
                    #Tn[i,j]=T[i,j]*A+Tn[i-1,j]*B+Tn[i+1,j]*C+Tn[i,j-1]*D+Tn[i,j+1]*E
            if(np.max(np.abs(To-Tn)/Tn)<1e-4):
                break
        T=Tn.copy()
    
    fig, ax = plt.subplots()
    ctf=ax.contourf(x,y,T,10,cmap='jet')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_aspect('equal')
    ax.set_xlim(0,1)
    ax.set_ylim(0,1)
    fig.colorbar(ctf)
    plt.show()

if __name__=='__main__':
    cubic()

将时间步长设为1200，同样可以得到收敛的结果：


若使用传统的追赶法，则大量使用for循环，拖慢了计算速度，时间步长虽然变为原来的两倍，计算时间反而增大了。改为矩阵形式后，计算速度得到极大的提升。

五、扩展

隐式格式显然可以适应更大的时间步长，但是隐式格式在单个时间步内需要多次迭代。两种方法求解相同问题实际所需计算时间跟所采用的算法也有一定的关系，并且二者精度相当，隐式格式并不会因为求解复杂而提高算法的精度。