计算流体力学（CFD）公式转化与电脑求解原理

以下是关于流体力学公式如何转化为计算软件能用的形式，以及电脑如何求解的详细解释，用简单语言面向初学者。

1. 如何转化为计算软件使用的公式

流体力学方程是连续的微分公式，电脑只能处理数字，所以需要把它们“翻译”成离散的加减法形式。这个过程叫离散化。

1.1 把空间和时间切成小块

• 空间网格: 把模型（比如管道）分成小格子，像搭乐高积木。例如，一个管道分成10万个小立方体。
• 时间步长: 把时间分成小段，比如每0.001秒算一次，像拍电影一帧一帧。 比喻: 把水流拍成照片。

1.2 把微分变成差分

微分是连续变化，电脑用格子间的差值代替。

• 时间项（如 ∂e/∂t）: 用前后时间差近似：

  ∂e/∂t ≈ (e_iⁿ⁺¹ – e_iⁿ) / Δt

  • e_iⁿ: 当前内能，e_iⁿ⁺¹: 下一步内能，Δt: 时间步长。
• 空间项（如 ∇ · (ρu)）: 用相邻格子差值：

  ∇ · (ρu) ≈ ((ρu)_i+1ⁿ – (ρu)_i-1ⁿ) / (2 Δx)

  • Δx: 格子宽度。

1.3 处理复杂项

• 对流项（如 u · ∇e）: u_x (e_i+1 – e_i-1) / (2 Δx)
• 扩散项（如 ∇ · (k ∇T)）: (k (T_i+1 – T_i) – k (T_i – T_i-1)) / (Δx²)

1.4 生成方程组

每个格子一个公式，10万个格子就有10万个方程，形成一个大系统。 比喻: 把水流故事拆成一页页小纸条。

2. 电脑如何执行求解的原理

离散化后，方程变成代数形式，电脑通过以下步骤求解。

2.1 初始化：给个起点

• 假设初始值，比如 u = 0、p = 大气压。 比喻: 像迷宫游戏先选个起点。

2.2 迭代求解：一步步逼近

• 显式方法: 用当前值算下一步，简单但不稳定：

  e_iⁿ⁺¹ = e_iⁿ + Δt [-p_i (∇ · u) + ∇ · (k ∇T) + Φ_i – ρ (u · ∇e)] / ρ

• 隐式方法: 把 n+1 的值写进方程，解矩阵，更稳定。
• 比喻: 像猜谜语，反复调整直到猜对。

2.3 矩阵求解：处理大方程组

• 隐式方法生成矩阵 A x = b，x 是 n+1 的值。
• 迭代法: 从猜测值开始调整，如Gauss-Seidel。

2.4 检查收敛

检查值变化是否小于 10^-6，确保守恒。

2.5 输出结果

生成彩图或动画，如温度分布。 比喻: 把水流拍成视频。

3. 为什么可以转化？更基础的讲解

3.1 流体力学方程是什么

这些公式是“水流规则”，描述连续变化，但电脑只认数字。

3.2 为什么从“连续”变“离散”

• 连续: 微分是平滑曲线，如 ∂ρ/∂t。 比喻: 像一条蜿蜒的河流。
• 离散: 把时间和空间切成小块，用差值代替微分。 比喻: 把河流分成小段直线。
• 为什么行: 块越小，拼起来越接近原曲线。

3.3 例子：连续性方程

∂ρ/∂t + ∂(ρu)/∂x = 0

• 时间项: (ρ₂^现在 – ρ₂^0.1秒前) / 0.1
• 空间项: (ρ₃ u₃ – ρ₁ u₁) / Δx
• 新公式: ρ₂^现在 = ρ₂^0.1秒前 – 0.1 * (ρ₃ u₃ – ρ₁ u₁) / Δx

3.4 电脑怎么算

• 从猜测开始，用公式算下一步，反复调整直到稳定。

4. 具体离散公式与推导

4.1 连续性方程

∂ρ/∂t + ∂(ρu)/∂x = 0

离散公式:

ρ_iⁿ⁺¹ = ρ_iⁿ – (Δt / (2 Δx)) * [(ρu)_i+1ⁿ – (ρu)_i-1ⁿ]

推导:

1. 时间项: ∂ρ/∂t ≈ (ρ_iⁿ⁺¹ – ρ_iⁿ) / Δt
2. 空间项: ∂(ρu)/∂x ≈ ((ρu)_i+1ⁿ – (ρu)_i-1ⁿ) / (2 Δx)
3. 整合: 代入并整理。

4.2 纳维-斯托克斯方程

ρ (∂u/∂t + u ∂u/∂x) = -∂p/∂x + μ ∂²u/∂x²

离散公式:

u_iⁿ⁺¹ = u_iⁿ + (Δt / ρ) * [-(p_i+1ⁿ – p_i-1ⁿ) / (2 Δx) + μ (u_i+1ⁿ – 2 u_iⁿ + u_i-1ⁿ) / (Δx²) – ρ u_iⁿ (u_i+1ⁿ – u_i-1ⁿ) / (2 Δx)]

推导:

1. 时间项: ∂u/∂t ≈ (u_iⁿ⁺¹ – u_iⁿ) / Δt
2. 对流项: u ∂u/∂x ≈ u_iⁿ (u_i+1ⁿ – u_i-1ⁿ) / (2 Δx)
3. 压力项: -∂p/∂x ≈ -(p_i+1ⁿ – p_i-1ⁿ) / (2 Δx)
4. 粘性项: μ ∂²u/∂x² ≈ μ (u_i+1ⁿ – 2 u_iⁿ + u_i-1ⁿ) / (Δx²)

4.3 能量方程

ρ (∂e/∂t + u ∂e/∂x) = -p ∂u/∂x + k ∂²T/∂x²

离散公式:

e_iⁿ⁺¹ = e_iⁿ + (Δt / ρ) * [-p_iⁿ (u_i+1ⁿ – u_i-1ⁿ) / (2 Δx) + k (T_i+1ⁿ – 2 T_iⁿ + T_i-1ⁿ) / (Δx²) – ρ u_iⁿ (e_i+1ⁿ – e_i-1ⁿ) / (2 Δx)]

推导:

1. 时间项: ∂e/∂t ≈ (e_iⁿ⁺¹ – e_iⁿ) / Δt
2. 对流项: u ∂e/∂x ≈ u_iⁿ (e_i+1ⁿ – e_i-1ⁿ) / (2 Δx)
3. 压力项: -p ∂u/∂x ≈ -p_iⁿ (u_i+1ⁿ – u_i-1ⁿ) / (2 Δx)
4. 热传导项: k ∂²T/∂x² ≈ k (T_i+1ⁿ – 2 T_iⁿ + T_i-1ⁿ) / (Δx²)

5. 隐式方法的矩阵形式

5.1 连续性方程（隐式）

((ρ_iⁿ⁺¹ – ρ_iⁿ) / Δt) + (((ρu)_i+1ⁿ⁺¹ – (ρu)_i-1ⁿ⁺¹) / (2 Δx)) = 0

矩阵形式 (假设 u 已知):

[1 a 0] [ρ₁ⁿ⁺¹] [ρ₁ⁿ + a (ρu)₀ⁿ⁺¹]
[-a 1 a] [ρ₂ⁿ⁺¹] = [ρ₂ⁿ ]
[0 -a 1] [ρ₃ⁿ⁺¹] [ρ₃ⁿ – a (ρu)₄ⁿ⁺¹], a = (Δt / (2 Δx)) uⁿ⁺¹

5.2 纳维-斯托克斯方程（隐式）

ρ (((u_iⁿ⁺¹ – u_iⁿ) / Δt) + u_iⁿ ((u_i+1ⁿ⁺¹ – u_i-1ⁿ⁺¹) / (2 Δx))) = -((p_i+1ⁿ⁺¹ – p_i-1ⁿ⁺¹) / (2 Δx)) + μ ((u_i+1ⁿ⁺¹ – 2 u_iⁿ⁺¹ + u_i-1ⁿ⁺¹) / (Δx²))

矩阵形式 (假设 p 已知):

[b c 0] [u₁ⁿ⁺¹] [右边项],
[-a b c] [u₂ⁿ⁺¹] = [右边项],
[0 -a b] [u₃ⁿ⁺¹] [右边项],
a = (ρ u_iⁿ Δt) / (2 Δx) + (μ Δt) / (Δx²), b = ρ / Δt + (2 μ Δt) / (Δx²), c = -(ρ u_iⁿ Δt) / (2 Δx) + (μ Δt) / (Δx²)

5.3 能量方程（隐式）

ρ (((e_iⁿ⁺¹ – e_iⁿ) / Δt) + u_iⁿ ((e_i+1ⁿ⁺¹ – e_i-1ⁿ⁺¹) / (2 Δx))) = -p_iⁿ⁺¹ ((u_i+1ⁿ⁺¹ – u_i-1ⁿ⁺¹) / (2 Δx)) + k ((T_i+1ⁿ⁺¹ – 2 T_iⁿ⁺¹ + T_i-1ⁿ⁺¹) / (Δx²))

矩阵形式 (假设 p, u 已知):

[b c 0] [e₁ⁿ⁺¹] [右边项],
[-a b c] [e₂ⁿ⁺¹] = [右边项],
[0 -a b] [e₃ⁿ⁺¹] [右边项],
a = (ρ u_iⁿ Δt) / (2 Δx) + (k Δt) / (c_v Δx²), b = ρ / Δt + (2 k Δt) / (c_v Δx²), c = -(ρ u_iⁿ Δt) / (2 Δx) + (k Δt) / (c_v Δx²)

6. 为什么隐式方法用 n+1

6.1 显式 vs 隐式

• 显式: 只用 n 的值猜 n+1。 比喻: 走路只看脚下，容易摔。
• 隐式: 用 n+1 的值一起解。 比喻: 走路看全路，计划每步。

6.2 为什么用 n+1

隐式站在“未来角度”看问题。

• 显式问题: 只看现在，水流突变可能算错。 例子: 跑步踩坑摔倒。
• 隐式优势: 考虑未来，调整所有格子。 例子: 跑步提前避坑。

6.3 为什么更准确

不是 n+1 本身准，而是整体解更符合规律。

• 显式: 误差像滚雪球。 例子: 只看今天猜天气，错得离谱。
• 隐式: 所有格子一起算，误差不放大。 例子: 看趋势猜天气，更靠谱。

6.4 为什么更稳定

• 显式: Δt 大时结果“爆炸”。 例子: 跳格子跳太远摔出界。
• 隐式: 矩阵约束，Δt 大也稳。 例子: 拉着手跳，不摔。

6.5 例子

• 显式: ρ_iⁿ⁺¹ = 1000 – 变化量，可能变负数。
• 隐式: 矩阵解出 ρ_iⁿ⁺¹ = 998，合理且稳定。