Python 有限体积法|极客笔记

Python 有限体积法

有限体积法是一种数值方法，用于求解偏微分方程组，特别适用于流体动力学、热传导等问题。在数值计算中，有限体积法将偏微分方程离散化为代数方程，通过在物理空间中定义有限体积，然后在网格中集成它们的守恒方程来进行计算。Python 是一种非常适合实现有限体积法的编程语言，其简洁的语法和丰富的科学计算库使得实现数值方法变得十分容易。

在本文中，我们将介绍如何使用 Python 实现有限体积法，并分步解释其原理和实现过程。

1. 一维热传导问题

我们首先考虑一个简单的一维热传导问题，即热方程：

$\frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2}$

其中 $u(x, t)$ 是温度分布函数， $x$ 是空间坐标， $t$ 是时间。我们将在一个有限的区域 $[0, L]$ 上求解该方程，假设边界条件为 $u(0, t) = 0, u(L, t) = 0$ ，初始条件为 $u(x, 0) = f(x)$ 。

我们将使用有限体积法离散化该方程，并通过数值求解得到温度分布的近似解。

2. 空间离散化

首先，我们将空间坐标 $x$ 离散化为 $N$ 个有限体积单元，每个单元之间的长度为 $\Delta x = \frac{L}{N}$ ，网格节点的坐标为 $x_i = i \cdot \Delta x$ ，其中 $i = 0, 1, \ldots, N$ 。

接下来我们定义每个有限体积单元的体积为 $V_i = \Delta x$ ，并通过积分得到离散化的守恒方程为：

$\frac{d}{dt} \int_{V_i} u \, dx = \frac{d}{dt} \left( u_i \cdot V_i \right) = -\int_{\partial V_i} F \cdot dS$

其中 $u_i$ 是节点 $x_i$ 上的温度值， $F$ 是热传导通量密度， $\partial V_i$ 是体积 $V_i$ 的边界面， $dS$ 是边界面的法向单位矢量。

3. 时间离散化

将时间坐标 $t$ 离散化为时间步长 $\Delta t$ ，我们可以使用前向差分法对时间导数进行近似，得到递推式：

$u_i^{n+1} = u_i^n + \frac{\Delta t}{V_i} \sum_{\text{faces}} F_{i+\frac{1}{2}} – F_{i-\frac{1}{2}}$

其中 $u_i^n$ 是节点 $x_i$ 上时间步 $n$ 的温度值， $F_{i+\frac{1}{2}}$ 和 $F_{i-\frac{1}{2}}$ 分别是体积单元 $[x_i, x_{i+1}]$ 和 $[x_{i-1}, x_i]$ 上的热传导通量。

4. 数值实现

接下来我们使用 Python 实现上述算法。首先，我们导入必要的库：

import numpy as np
import matplotlib.pyplot as plt

然后定义模拟参数和网格：

L = 1.0  # 区域长度
N = 100  # 网格数
dx = L / N
dt = 0.001  # 时间步长
T = 0.1  # 总时间
x = np.linspace(0, L, N+1)  # 网格节点坐标
u = np.zeros(N+1)  # 温度分布
u_new = np.zeros(N+1)  # 新一步的温度分布

下面我们定义初始条件和边界条件：

def initial_condition(x):
    return np.sin(np.pi * x)

def boundary_condition(u):
    u[0] = 0
    u[N] = 0

然后对时间进行迭代，并在每个时间步计算新的温度分布：

time = 0.0
while time < T:
    # 边界条件
    boundary_condition(u)

    # 更新温度分布
    u_new[1:-1] = u[1:-1] + dt/dx**2 * (u[:-2] - 2*u[1:-1] + u[2:])

    # 更新时间步
    time += dt
    u[:] = u_new

最后我们可以将结果可视化出来：

plt.plot(x, u)
plt.xlabel('x')
plt.ylabel('Temperature')
plt.title('One-dimensional Heat Conduction')
plt.show()

Python 有限体积法