(四)、传热学-热传导问题的数值解法

1、基本思想

把原来在时间、空间坐标系中连续的物理的场，用有限个离散点上的值的集合来代替，通过求解按一定方法建立起来的关于这些值的代数方程，来获得离散点上被求物理量的值。

2、基本步骤

(1)、建立控制方程及定解条件；

(2)、区域离散化；

(3)、建立节点物理量的代数方程；

(4)、设立迭代初场；

(5)、求解代数方程组；

(6)、解的分析。

3、内节点离散方程的建立方法

(1)、泰勒级数展开法

以节点(m,n)处的二阶偏导数为例。

(2)、热平衡法

对每个节点所代表的元体用傅里叶导热定律写出其能量守恒表达式。把节点看成是元体的代表。

通过元体的界面所传导的热流量可以对有关的两个节点应用傅里叶定律写出。

Φw=λAΔtΔxΦw=λΔy⋅1⋅tm−1,n−tm,nΔxΦw=λAΔtΔxΦw=λΔy⋅1⋅tm−1,n−tm,nΔx

类似的可以写出通过其它三个界面e、n、s，而传导给节点(m,n)的热量。对于所研究的问题，元体(m,n)的能量守恒方程为：

Φw+Φe+Φn+Φs=0Φw+Φe+Φn+Φs=0

带入得：

λΔytm−1,n−tm,nΔx+λΔytm+1,n−tm,nΔx+λΔxtm,n+1−tm,nΔy+λΔxtm,n−1−tm,nΔy=0λΔytm−1,n−tm,nΔx+λΔytm+1,n−tm,nΔx+λΔxtm,n+1−tm,nΔy+λΔxtm,n−1−tm,nΔy=0

4、边界节点离散方程的建立代数方程的求解

(1)、边界节点离散方程的建立

内部角点最具代表：

λΔytm−1,n−tm,nΔx+λΔy2tm+1,n−tm,nΔx+λΔxtm,n+1−tm,nΔy+λΔx2tm,n−1−tm,nΔy+3ΔxΔy4˙Φm,n+Δx+Δy2qw=0λΔytm−1,n−tm,nΔx+λΔy2tm+1,n−tm,nΔx+λΔxtm,n+1−tm,nΔy+λΔx2tm,n−1−tm,nΔy+3ΔxΔy4Φ˙m,n+Δx+Δy2qw=0

当Δx=ΔyΔx=Δy时有：

tm,n=16(2tm−1,n+2tm,n+1+tm,n−1+tm+1,n+3Δx22λ˙Φm,n+2Δxλqw)tm,n=16(2tm−1,n+2tm,n+1+tm,n−1+tm+1,n+3Δx22λΦ˙m,n+2Δxλqw)

其中qwqw可以按照绝热边界、不为0、对流边界，三种情况进行代入。

(2)、处理不规则区域得阶梯形逼近

(3)、求解代数方程的迭代法

5、非稳态导热问题的数值解法

(1)、时间-空间区域的离散化(泰勒展开法)

将函数t在节点(n,i+1)对点(n,i)作泰勒展开，可有：

t(i+1)n=t(i)n+Δτ∂t∂τ|n,i+Δτ22∂2t∂τ2|n,i+...t(i+1)n=t(i)n+Δτ∂t∂τ∣∣n,i+Δτ22∂2t∂τ2∣∣n,i+...

于是有：

∂t∂τ|n,i=t(i+1)n−t(i)nΔτ+O(Δτ)∂t∂τ∣∣n,i=t(i+1)n−t(i)nΔτ+O(Δτ)

式中，符号O(Δτ)O(Δτ)表示余项中ΔτΔτ的最低阶为一次。由上式可在点(n,i)处一阶导数的一种差分表示式。

∂t∂τ|n,i=t(i+1)n−t(i)nΔτ∂t∂τ∣∣n,i=t(i+1)n−t(i)nΔτ(向前差分)

另外可得向后差分和中心差分：

∂t∂τ|n,i=t(i)n−t(i−1)nΔτ∂t∂τ∣∣n,i=t(i)n−t(i−1)nΔτ(向后差分)

∂t∂τ|n,i=t(i+1)n−t(i−1)n2Δτ∂t∂τ∣∣n,i=t(i+1)n−t(i−1)n2Δτ(中心差分)

另外可得稳态和非稳态二阶中心差分：

∂2t∂x2=tn+1−2tn+tn−1Δx2∂2t∂x2=tn+1−2tn+tn−1Δx2(稳态二阶中心差分)

∂2t∂x2=tin+1−2tin+tin−1Δx2∂2t∂x2=tin+1−2tin+tin−1Δx2(非稳态二阶中心差分)

(2)、显式格式

对于一维非稳态导热方程，如扩散项取中心差分(并取i时刻之值)，非稳态项取向前差分，则有：

t(i+1)n−t(i)nΔτ=atin+1−2tin+tin−1Δx2t(i+1)n−t(i)nΔτ=atin+1−2tin+tin−1Δx2

显示差分格式：求解非稳态导热方程就是从已知的初始温度分布出发，根据边界条件依次求得以后各个时间层上的温度值。一旦i时层上各个节点的温度已知，可立即算出i+1时层上各内点的温度，而不必求解联立的方程。

优缺点：计算工作量小，对时间步长及空间步长有一定的限制，否则会出现不合理的震荡的解，称为稳定性问题。

(3)、隐式格式

对于一维非稳态导热方程，如扩散项取中心差分(并取i+1时刻之值)，非稳态项取向前差分，则有：

t(i+1)n−t(i)nΔτ=at(i+1)n+1−2t(i+1)n+t(i+1)n−1Δx2t(i+1)n−t(i)nΔτ=at(i+1)n+1−2t(i+1)n+t(i+1)n−1Δx2

隐式差分格式：已知的是i时层的值t(i)nt(i)n，而未知量有3个，因此不能直接算出t(i+1)nt(i+1)n之值，而必须求解(i+1)时层的一个联立方程组才能得出(i+1)时层各节点的温度。

优缺点：计算工作量大，但对步长没有限制，不会出现解的振荡现象。

(4)、边界节点的离散方程(能量守恒法)

无限大平板得右边界部分，其表面流体冷却，表面传热系数为h。

对元体应用能量守恒定律得：

Φright+Φleft=ΦchangΦright+Φleft=Φchang

λΔytiN−1−tiNΔx+hΔy(tf−tiN)=ρcV∂t∂τ=ρcΔx2Δyti+1N−tiNΔτλΔytiN−1−tiNΔx+hΔy(tf−tiN)=ρcV∂t∂τ=ρcΔx2Δyti+1N−tiNΔτ

(5)、数值稳定性

保证合理性，满足稳定性要求。