有限体积法：基于一维稳态扩散问题及其程序实现
有限体积法：基于一维稳态扩散问题及其程序实现
有限体积法（Finite Volume Method, FVM）是一种常用的数值方法，广泛应用于求解偏微分方程，特别是在计算流体动力学（CFD）中。本文将以一个经典的一维稳态导热问题为例，详细介绍有限体积法的应用步骤，并结合Matlab编程实现，展示其求解过程。
1. 问题描述
我们考虑一个无内热源的一维稳态导热问题，如下图所示：
- 棒的长度为L = 0.5 m
- 截面积为A = 10 × 10−3 m2
- 左端温度保持在TA = 100∘C
- 右端温度保持在TB = 500∘C
- 棒的导热系数为k = 1000 W/(m⋅K)
求解该棒在稳态下的温度分布。
2. 有限体积法的应用步骤
有限体积法求解一维稳态导热问题的步骤可以分为以下三步：
2.1 步骤一：生成离散网格
首先，将求解域划分为若干个控制体积。这里将棒分为5个控制体积，每个体积对应一个节点。节点间的距离δx = 0.1 m。如下图所示，这些节点标号为1到5。
2.2 步骤二：构造离散方程
根据导热问题的控制方程：
d/dx(kAdT/dx) = 0
在没有内热源的情况下，导热方程可以简化为：
kAd2T/dx2 = 0
在离散化过程中，节点2、3和4的离散方程如下：
aPTP = aWTW + aETE
其中aW = kA/δx，aE = kA/δx，aP = aW + aE。
对于边界节点1和5，离散方程的处理稍有不同。以节点1为例，其离散方程为：
(k/δxA + 2k/δxA)T1 = k/δxA T2 + 2k/δxA TA
通过类似方法，节点5的离散方程可以表示为：
(k/δxA + 2k/δxA)T5 = k/δxA T4 + 2k/δxA TB
将离散方程组合成一个线性方程组：
300T1 = 100T2 + 200TA
200T2 = 100T1 + 100T3
200T3 = 100T2 + 100T4
200T4 = 100T3 + 100T5
300T5 = 100T4 + 200TB
写成矩阵形式有
[300 -100 0 0 0
-100 200 -100 0 0
0 -100 200 -100 0
0 0 -100 200 -100
0 0 0 -100 300]
{T1
T2
T3
T4
T5} =
{200TA
0
0
0
200TB}
2.3 步骤三：解方程组
将TA = 100, TB = 500代入，解此方程组，可得
{T1
T2
T3
T4
T5} =
{140
220
300
380
460}
此两端固定温度值的绝热棒稳态分析解为线性温度分布：
T = 800x + 100
下图为分析解与数值解的图形比较，从中可以看出数值解有很高的计算精度（随着网格加密精度将进一步提高）。
3. 编程实现
以下是基于Matlab的编程实现：
function temperatureDistribution = computeHeatDistribution()
% 参数初始化
lengthRod = 0.5; % 棒的长度
thermalConductivity = 1000; % 导热系数
temperatureLeft = 100; % 左端温度
temperatureRight = 500; % 右端温度
numNodes = 5; % 节点数
deltaX = lengthRod / (numNodes - 1); % 控制容积宽度
% 系数矩阵和右端项初始化
coeffMatrix = zeros(numNodes, numNodes);
rhs = zeros(numNodes, 1);
% 内部节点的离散方程
for i = 2:numNodes-1
coeffMatrix(i, i-1) = thermalConductivity / deltaX;
coeffMatrix(i, i) = -2 * thermalConductivity / deltaX;
coeffMatrix(i, i+1) = thermalConductivity / deltaX;
end
% 边界条件
coeffMatrix(1, 1) = 3 * thermalConductivity / deltaX;
coeffMatrix(1, 2) = -thermalConductivity / deltaX;
rhs(1) = 2 * thermalConductivity / deltaX * temperatureLeft;
coeffMatrix(numNodes, numNodes-1) = -thermalConductivity / deltaX;
coeffMatrix(numNodes, numNodes) = 3 * thermalConductivity / deltaX;
rhs(numNodes) = 2 * thermalConductivity / deltaX * temperatureRight;
% 求解线性方程组
temperatureDistribution = coeffMatrix \ rhs;
% 计算分析解
x = linspace(0, lengthRod, 100);
analyticalSolution = 800 * x + 100;
% 绘图
figure;
plot(x, analyticalSolution, '--', 'LineWidth', 1.5, 'Color', [0, 0.6, 0.3]);
hold on;
plot(linspace(deltaX/2, lengthRod-deltaX/2, numNodes), temperatureDistribution, 'o', 'MarkerSize', 8 , 'LineWidth', 1.5, 'MarkerSize', 8, 'Color', [0.7, 0.2, 0.2]);
hold on;
% 添加边界点
plot([0, lengthRod], [temperatureLeft, temperatureRight], 'o', 'MarkerSize', 8, 'LineWidth', 1.5, 'Color', [0.7, 0.2, 0.2]);
xlabel('位置 (m)', 'FontSize', 12, 'FontWeight', 'bold');
ylabel('温度 (°C)', 'FontSize', 12, 'FontWeight', 'bold');
title('一维稳态导热问题的温度分布', 'FontSize', 14, 'FontWeight', 'bold');
legend('精确解', '数值解', 'Location', 'Best');
grid on;
set(gca, 'FontSize', 10);
hold off;
% 输出数值解
fprintf('数值解：\n');
disp(temperatureDistribution);
end
该代码通过构建系数矩阵和右端项，利用Matlab内置的线性方程求解函数，得到了每个节点的温度分布。最终的温度分布图如图2-5所示，数值解与理论解具有良好的吻合性，验证了有限体积法在该问题中的有效性。
4. 总结
本文通过详细解析一维稳态导热问题的有限体积法求解过程，结合一维稳态扩散问题展示了从问题描述到离散化再到数值求解的完整流程。通过该案例，读者可以更深入地理解有限体积法的基本原理和应用场景，为进一步的数值计算研究奠定基础。