当前位置：首页 > article >正文

有限体积法：基于一维稳态扩散问题及其程序实现

article 2024/9/17 3:55:17

文章目录

- - 有限体积法：基于一维稳态扩散问题及其程序实现
  - - 1. 问题描述
    - 2. 有限体积法的应用步骤
    - - 2.1 步骤一：生成离散网格
      - 2.2 步骤二：构造离散方程
      - 2.3 步骤三：解方程组
    - 3. 编程实现
    - 4. 总结

有限体积法：基于一维稳态扩散问题及其程序实现

在工程与科学计算中，有限体积法（Finite Volume Method, FVM）是一种常用的数值方法，广泛用于求解偏微分方程，特别是在计算流体动力学（CFD）中。本文将以一个经典的一维稳态导热问题为例，详细介绍有限体积法的应用步骤，并结合 Matlab 编程实现，展示其求解过程。

1. 问题描述

我们考虑一个无内热源的一维稳态导热问题，如下图所示：

棒的长度为 $0.5\ \text{m}$ ，截面积为 $\times 10^{-3}\ \text{m}^2$ ；
左端温度保持在 $T_A = 100^\circ\text{C}$ ，右端温度保持在 $T_B = 500^\circ\text{C}$ ；
棒的导热系数为 $1000\ \text{W}/(\text{m}\cdot\text{K})$ 。

求解该棒在稳态下的温度分布。
稳态导热棒

2. 有限体积法的应用步骤

有限体积法求解一维稳态导热问题的步骤可以分为以下三步：

2.1 步骤一：生成离散网格

首先，将求解域划分为若干个控制体积。这里将棒分为 5 个控制体积，每个体积对应一个节点。节点间的距离 $\delta x = 0.1\ \text{m}$ 。如下图所示，这些节点标号为 1 到 5。
绝热棒离散网络

2.2 步骤二：构造离散方程

根据导热问题的控制方程：

$\frac{d}{dx}\left( kA\frac{dT}{dx} \right) = 0$

在没有内热源的情况下，导热方程可以简化为：

$\frac{d^2 T}{dx^2} = 0$

在离散化过程中，节点 2、3 和 4 的离散方程如下：

$a_P T_P = a_W T_W + a_E T_E$

其中 $a_W = \frac{kA}{\delta x}$ ， $a_E = \frac{kA}{\delta x}$ ， $a_P = a_W + a_E$ 。

对于边界节点 1 和 5，离散方程的处理稍有不同。以节点 1 为例，其离散方程为：

$\left( \frac{k}{\delta x}A + \frac{2k}{\delta x}A \right) T_1 = \frac{k}{\delta x}A T_2 + \frac{2k}{\delta x}A T_A$

通过类似方法，节点 5 的离散方程可以表示为：

$\left( \frac{k}{\delta x}A + \frac{2k}{\delta x}A \right) T_5 = \frac{k}{\delta x}A T_4 + \frac{2k}{\delta x}A T_B$

将离散方程组合成一个线性方程组：

$\begin{aligned} 300 T_1 &= 100 T_2 + 200 T_A \\ 200 T_2 &= 100 T_1 + 100 T_3 \\ 200 T_3 &= 100 T_2 + 100 T_4 \\ 200 T_4 &= 100 T_3 + 100 T_5 \\ 300 T_5 &= 100 T_4 + 200 T_B \end{aligned}$
写成矩阵形式有

$\left[\begin{array}{ccccc} 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 \end{array}\right]\left\{\begin{array}{c} T_1 \\ T_2 \\ T_3 \\ T_4 \\ T_5 \end{array}\right]=\left\{\begin{array}{c} 200 T_A \\ 0 \\ 0 \\ 0 \\ 200 T_B \end{array}\right\}$

2.3 步骤三：解方程组

将 $T_A=100, T_B=500$ 代入, 解此方程组, 可得

$\left\{\begin{array}{l} T_1 \\ T_2 \\ T_3 \\ T_4 \\ T_5 \end{array}\right\}=\left\{\begin{array}{l} 140 \\ 220 \\ 300 \\ 380 \\ 460 \end{array}\right\}$

此两端固定温度值的绝热棒稳态分析解为线性温度分布：

$T = 800 x + 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