一维稳态无源导热——从方程到代码一维稳态无源导热——从方程到代码知乎用户n6NGzC知乎用户n6NGzC​关注她45 人赞同了该文章1. 问题描述如图，有绝热等截面棒，长 �=0

一维稳态无源导热——从方程到代码 - 知乎 (zhihu.com)

一维稳态无源导热——从方程到代码

知乎用户n6NGzC

​关注她

45 人赞同了该文章

1. 问题描述

如图，有绝热等截面棒，长 �=0.5� ，截面 �=10×10−3�2 ，棒左端温度 �� 保持100K，右端温度 �� 保持500K。棒的材料导热系数 �=1000�/(�.�) 。求绝热棒在稳态时温度分布。

问题描述

2. 问题分析

一维稳态无源导热控制方程如下：

���(�����)=0

对于上述方程，利用有限体积法来进行离散求解。

离散求解PED方程的步骤为：

1. 离散控制域（网格划分）
2. 在每一个控制体上离散控制方程
3. 插值得到界面值，完成单元离散方程
4. 组装单元控制方程，形成整体控制方程组（Ax=b）
5. 求解代数方程组（直接，迭代）
6. 得到离散场变量

下面分别按照上面给出的六个步骤进行分析。

3. 问题求解

3.1 域离散

采用均匀网格将棒沿着长度方向离散，如图所示：

域离散

为了编写程序，我们需要容器来储存网格数据。对于一维问题，我们需要存储离散节点的位置坐标，单元的体积（在本例中为单元的长度乘以截面积）。

代码实现

import numpy as np

# 域离散
L = 0.5       # 棒总长
N = 5         # 控制单元数目
dx = L / N    # 控制单元长度
centriod = [i * dx + dx / 2 for i in range(N)]   # 控制单元中心坐标
A = [0.001] * N                                  # 控制单元中心对应的界面面积
k = [1000] * N                                   # 控制体单元中心处的导热系数
T = [0] * N                                      # 控制体单元中心温度初始值。可以任意给定
A = np.zeros((N, N))                             # 整体系数矩阵
b = np.zeros(N)                                  # 整体方程右端项
TA = 100
TB = 500

3.2 方程离散

在如图的控制体上，对控制方程进行积分得到：

∫�����(�����)=(�����)���−(�����)���=0

上式中间部分继续展开，得到：

����(����)�−����(����)�=0

在FVM中，材料属性和场变量都是存储在节点上的，但(3)式中需要到界面上的导热系数和温度梯度值。我们需要利用插值函数将节点值插值到界面上。采用不同的插值函数，就可以得到不同的离散方程，不同的离散方程其精度、收敛性和稳定性都不一样。这里不进行深入的讨论（如果讨论可能需要一大本书，并且有些问题也讨论不清楚。）

3.3 界面值插值

我们采用如下的差值方式：

对于导热系数（材料参数），常采用线性加权方法插值。如下：
��=����+(1−��)��
其中，gc​为一个几何系数，定义如下：

��=||��||||��||

对于均匀网格而言，gc​=1/2。所以，对于本例中，界面上的导热系数利用如下插值得到：

��=��+��2��=��+��2

对于界面上的梯度值，也有多种插值方法，我们采用简单的线性插值，如下：

(����)�=��−�����(����)�=��−�����

其中， ��� 表示C点到E点之间的距离。在本例中，采用均匀网格， ���=���=Δ�=�� ​。

对于界面面积： ��=��=�

3.4 单元方程整理

将导热系数和温度梯度在界面的插值表达式带入到(3)中，得到最后的单元离散方程如下：

��+��2����−��Δ�−��+��2����−��Δ�=0

整理得到如下简洁形式：

����+����+����=����=−(��+��)��=−��+��2Δ���=−��+��2Δ���=0

显然，对于每一个控制体，都用到了它自己的前后单元值。那么，我们就需要对边界处做特殊的处理。因为边界处的单元只有一个相邻单元。

代码实现

a_w = -(k[ii-1] + k[ii]) / 2 / dx
a_e = -(k[ii] + k[ii+1]) / 2 / dx
bc = 0

3.5 整体矩阵组装

通过上面的分析，我们成功得到了单元的离散形式的控制方程。因为单元控制方程在离散的时候要用到相邻单元的值，但是相邻单元的值也是未知的。因此这个时候我们可以如下的方法来处理：

1. 将相邻单元的温度值当做已知（一般给出一个猜测的场变量值），这样就可以显式的解得当前单元处的温度值，但是这样做要对全场进行多次迭代。
2. 将每个单元的控制方程组装，形成一个对应于全场未知量的系数矩阵，求解这个整体的代数方程组，统一解得场变量。

这两种方式本质其实是相同的。尤其是当整体系数矩阵较大的时候，CFD中控制方程的非线性导致直接求解方程组很困难，一般常采用迭代求解。那么对于整体方程组采用迭代求解就和第一种思路一样。

整体矩阵组装的时候，将控制体进行编号，未知温度值当做列向量，形成AT=b的形式。具体的组装方式在后面程序中给出。

代码实现：

# 方程离散和整体组装
for ii in range(1, N-1):  # 用来组装整体系数矩阵的整体循环，注意边界单元没有循环
    # C对应下标ii, W对应下标ii-1, E对应下标ii+1
    a_W = -(k[ii-1] + k[ii]) / 2 / dx
    a_E = -(k[ii] + k[ii+1]) / 2 / dx
    b_C = 0
    A[ii, ii-1] = a_W
    A[ii, ii+1] = a_E
    A[ii, ii] = -(a_W + a_E)
    b[ii] = b_C

3.6 边界处理

对于靠近A的边界单元：

从方程 ����(����)�−����(����)�=0 开始。此时，控制体示意图如下：

A端边界离散

��+��2���−��Δ�−�����−��Δ�/2=0(��+��2Δ�+��Δ�/2)��+(−��Δ�/2)��+(−��+��2Δ�)��=0����+����=����=−��+������=−��+��2Δ���=(��Δ�/2)��

代码实现：

# 第一个边界单元修正， 此时C对应下标0， k_A = k[0]
A[0, 1] = -(k[0] + k[1]) / 2 / dx
b[0] = 2 * k[0] / dx * TA
A[0, 0] = -A[0, 1] + b[0] / TA

对于靠近B的边界：

从方程 ����(����)�−����(����)�=0 开始。此时，控制体示意图如下：

B端边界离散

�����−��Δ�/2−��+��2���−��Δ�=0(��+��2Δ�+��Δ�/2)��+(−��Δ�/2)��+(−��+��2Δ�)��=0����+����=����=−��+������=−��+��2Δ���=(��Δ�/2)��

代码实现：

# 第二个边界单元修正,此时，C对应下标N-1, W对应下标N-2, k_B = k[N-1]
A[N-1, N-2] = -(k[N-2] + k[N-1]) / 2 / dx
b[N-1] = k[N-1] * 2 / dx * TB
A[N-1, N-1] = -A[N-1, N-2] + b[N-1] / TB

至此，我们已经得到了与原来由PDE加边界条件定义的问题近似一致的离散型代数方程组AT=b。对于这类代数方程组的求解，有许许多多的算法可用。我们的问题最终形成的是一个三对角矩阵。可以利用TDM方法求解。这里用python-numpy库自带的函数来求解。

如果本例中，我们将控制单元设置为5个，则可以得到最终的系数矩阵如下：

[[ 30000., -10000.,      0.,      0.,      0.],
 [-10000.,  20000., -10000.,      0.,      0.],
 [     0., -10000.,  20000., -10000.,      0.],
 [     0.,      0., -10000.,  20000., -10000.],
 [     0.,      0.,      0., -10000.,  30000.]]

右端列向量为：

[ 2000000.,        0.,        0.,        0., 10000000.]

3.7 矩阵求解

矩阵求解的方法很多。一般采用迭代求解。

# 代数方程组求解
T = np.linalg.solve(A, b)

3.8 结果可视化

import matplotlib.pyplot as plt
fig = plt.figure()
plt.plot(np.arange(N)/N*L, T, 'o')
plt.show()

离散数值解

3.9 结果分析

本问题有理论解，表达式如下：
T=800x+100
将数值解和理论解对比，如下：

import matplotlib.pyplot as plt
x = np.array([0, *centriod, L])
T_all = np.array([TA, *T, TB])
exact = 800 * x + 100
fig = plt.figure()
plt.plot(x, T_all, 'o', label='FVM')
plt.plot(x, exact, '-', label='Exact')
plt.legend()
plt.show()

数值解与精确解对比

4. 代码

完整代码如下：

import numpy as np

# 域离散
L = 0.5       # 棒总长
N = 5       # 控制单元数目
dx = L / N    # 控制单元长度
centriod = [i * dx + dx / 2 for i in range(N)]   # 控制单元中心坐标
A = [0.001] * N                                  # 控制单元中心对应的界面面积
k = [1000] * N                                   # 控制体单元中心处的导热系数
T = [0] * N                                      # 控制体单元中心温度初始值。可以任意给定
A = np.zeros((N, N))                             # 整体系数矩阵
b = np.zeros(N)                                  # 整体方程右端项
TA = 100
TB = 500


# 方程离散和整体组装
for ii in range(1, N-1):  # 用来组装整体系数矩阵的整体循环，注意边界单元没有循环
    A[ii, ii-1] = -(k[ii-1] + k[ii]) / 2 / dx
    A[ii, ii+1] = -(k[ii]+k[ii+1]) / 2 / dx
    A[ii, ii] = -(A[ii, ii-1] + A[ii, ii+1])
    b_C = 0
    b[ii] = b_C

# 第一个边界单元修正
A[0, 1] = -(k[0] + k[1]) / 2 / dx
b[0] = 2 * k[0] / dx * TA
A[0, 0] = -(A[0, 1]) + b[0] / TA

# 第二个边界单元修正
A[N-1, N-2] = -(k[N-2]+k[N-1]) / 2 / dx
b[N-1] = k[N-1] * 2 / dx * TB
A[N-1, N-1] = -(A[N-1, N-2]) + b[N-1] / TB

# 方程求解
T = np.linalg.solve(A, b)

# 结果可视化
import matplotlib.pyplot as plt
fig = plt.figure()
plt.plot(np.arange(N)/N*L, T, 'o')
plt.show()


# 结果分析
import matplotlib.pyplot as plt
x = np.array([0, *centriod, L])
T_all = np.array([TA, *T, TB])
exact = 800 * x + 100
fig = plt.figure()
plt.plot(x, T_all, 'o', label='FVM')
plt.plot(x, exact, '-', label='Exact')
plt.legend()
plt.show()