流体力学核心贯通课 II
2025–2026 学年第二学期
北京理工大学 空天科学与技术学院
任杰
算例实践：Couette 流动
参考书
- John D. Anderson Computational Fluid Dynamics:
The Basics with Applications
McGraw-Hill, 1995 - Chapter 9
真正理解 CFD 的最好方式，是自己实现数值算法并运行计算 - John D. Anderson
流动描述
库埃特流动定义如下：考虑两个相距垂直距离 \(D\) 的平行平板之间的粘性流动，如图所示。上平板以速度 \(u_e\) 运动，下平板静止。
- 两平板之间的流场完全由运动上平板施加在流体上的剪切应力驱动
- 库埃特流动可能是所有粘性流动中最简单的
- 将采用Crank–Nicolson（CN隐式）格式求解
控制方程
连续性方程（定常流）
\[ \frac{\partial (\rho u)}{\partial x}+ \frac{\partial (\rho v)}{\partial y} = 0 \]
考虑流向充分发展的流动：
\[ \frac{\partial}{\partial x} = 0 \Longrightarrow \frac{\partial (\rho v)}{\partial y} = 0 \]
考虑密度为常数的情况，根据边界条件：
\[ v = 0 \quad (y=0) \Longrightarrow v = 0 \quad \text{处处成立} \]
- 无法向速度（流线为平行直线）
控制方程
y方向动量方程
\[ \rho \frac{Dv}{Dt}= -\frac{\partial p}{\partial y}+ \frac{\partial \tau_{xy}}{\partial x}+ \frac{\partial \tau_{yy}}{\partial y}+ \frac{\partial \tau_{zy}}{\partial z}+ \rho f_y \]
其中，
\[ \rho \frac{Dv}{Dt}=\rho\left(\frac{\partial v}{\partial t}+ u \frac{\partial v}{\partial x}+ v \frac{\partial v}{\partial y}+ w \frac{\partial v}{\partial z}\right) \]
对于 Couette 流，方程简化为：
\[ 0 = -\frac{\partial p}{\partial y} + \frac{\partial \tau_{yy}}{\partial y}, \qquad \tau_{yy} = 2\mu \frac{\partial v}{\partial y} + \lambda \left(\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y}\right) = 0 \]
得到：
\[ \frac{\partial p}{\partial y} = 0 \]
控制方程
x方向动量方程
\[ \rho \frac{Du}{Dt}= -\frac{\partial p}{\partial x}+ \frac{\partial \tau_{xx}}{\partial x}+ \frac{\partial \tau_{yx}}{\partial y}+ \frac{\partial \tau_{zx}}{\partial z}+ \rho f_x \]
简化后的控制方程
\[ \frac{\partial}{\partial y} \left( \mu \frac{\partial u}{\partial y} \right) = 0 \]
若 \(\mu\) 为常数：
\[ \frac{\partial^2 u}{\partial y^2} = 0 \]
两次积分：
\[ u = c_1 y + c_2 \]
控制方程
x方向动量方程
两次积分：
\[ u = c_1 y + c_2 \]
根据边界条件：
- \(u(0) = 0\)
- \(u(D) = u_e\)
得到：
\[ u = \frac{u_e}{D} y \]
👉 线性速度分布
非定常 Couette 流
控制方程：
\[ \rho \frac{\partial u}{\partial t} = \mu \frac{\partial^2 u}{\partial y^2} \]
👉 抛物型方程（类似热传导方程）
无量纲化
\[ u' = \frac{u}{u_e}, \quad y' = \frac{y}{D}, \quad t' = \frac{t}{D/u_e} \]
得到：
\[ \frac{\partial u'}{\partial t'} = \frac{1}{Re_D} \frac{\partial^2 u'}{\partial y'^2}, \qquad Re_D = \frac{\rho u_e D}{\mu} \]
数值求解
为方便起见，以下考虑无量纲方程，不再标注上标
\[ \frac{\partial u}{\partial t} = \frac{1}{Re_D} \frac{\partial^2 u}{\partial y^2} \]
计算域 \(y\in[0,1]\) 边界条件
- \(u(0) = 0\)
- \(u(1) = 1\)
精确解 \(u=y\)
Crank–Nicolson 格式离散：
\[ \frac{u_j^{n+1} - u_j^n}{\Delta t}= \frac{1}{2 Re_D (\Delta y)^2} \left[ (u_{j+1}^{n+1} - 2u_j^{n+1} + u_{j-1}^{n+1}) + (u_{j+1}^{n} - 2u_j^{n} + u_{j-1}^{n}) \right] \]
数值求解
\[ \frac{u_j^{n+1} - u_j^n}{\Delta t}= \frac{1}{2 Re_D (\Delta y)^2} \left[ (u_{j+1}^{n+1} - 2u_j^{n+1} + u_{j-1}^{n+1}) + (u_{j+1}^{n} - 2u_j^{n} + u_{j-1}^{n}) \right] \]
整理得到三对角形式：
\[ A u_{j-1}^{n+1} + B u_j^{n+1} + A u_{j+1}^{n+1} = K_j \]
其中：
\[ A = -\frac{\Delta t}{2(\Delta y)^2 Re_D}, \qquad B = 1 + \frac{\Delta t}{(\Delta y)^2 Re_D}, \qquad K_j = \left[1 - \frac{\Delta t}{(\Delta y)^2 Re_D}\right] u_j^n+ \frac{\Delta t}{2(\Delta y)^2 Re_D}\left( u_{j+1}^n + u_{j-1}^n \right) \]
数值求解
边界条件 \(u_1\) 和 \(u_{N+1}\) 是已知的：
\[ u_1 = 0 \]
\[ u_{N+1} = 1 \]
因此方程组包含 \(N-1\) 个方程，对应 \(N-1\) 个未知量，即：
\[ u_2, u_3, \dots, u_N \]
数值求解
我们可以将该方程组展开为更具体的形式。
第一条方程：
\[ A u_1^{n+1} + B u_2^{n+1} + A u_3^{n+1} = K_2 \]
由于 \(u_1 = 0\)，则有：
\[ B u_2^{n+1} + A u_3^{n+1} = K_2 \]
最后一条方程：
\[ A u_{N-1}^{n+1} + B u_N^{n+1} + A u_{N+1}^{n+1} = K_N \]
由于 \(u_{N+1} = 1\)，则有：
\[ A u_{N-1}^{n+1} + B u_N^{n+1} = K_N - A \]
数值求解
由此，方程组可以写成矩阵形式：
\[ \begin{bmatrix} B & A & 0 & 0 & 0 & \cdots & 0 \\ A & B & A & 0 & 0 & \cdots & 0 \\ 0 & A & B & A & 0 & \cdots & 0 \\ 0 & 0 & A & B & A & \cdots & 0 \\ \vdots & & & \ddots & & & \vdots \\ 0 & \cdots & 0 & A & B & A \\ 0 & \cdots & 0 & 0 & A & B \end{bmatrix} \begin{bmatrix} u_2^{n+1} \\ u_3^{n+1} \\ u_4^{n+1} \\ \vdots \\ u_N^{n+1} \end{bmatrix}= \begin{bmatrix} K_2 - 0 \\ K_3 \\ K_4 \\ \vdots \\ K_{N-1} \\ K_N - A u_e \end{bmatrix} \]
求解该线性方程组后，可得到 \(u_2^{n+1},\; u_3^{n+1},\; \dots,\; u_N^{n+1}\)，即时间步 \(n+1\) 的速度分布。
随后：👉 重复该过程（时间推进）直到速度分布收敛到稳态
数值求解
\[ \frac{\partial u}{\partial t} = \frac{1}{Re_D} \frac{\partial^2 u}{\partial y^2} \]
回顾上节课内容：显式格式的稳定性条件为：
\[ \frac{1}{Re_D} \frac{\Delta t}{(\Delta y)^2} \le \frac{1}{2} \]
我们定义参数：
\[ E = \frac{\Delta t}{Re_D (\Delta y)^2} \]
数值求解
\[ E = \frac{\Delta t}{Re_D (\Delta y)^2} \]
将其代入系数表达式，
\[ A = -\frac{\Delta t}{2(\Delta y)^2 Re_D}, \qquad B = 1 + \frac{\Delta t}{(\Delta y)^2 Re_D}, \qquad K_j = \left[1 - \frac{\Delta t}{(\Delta y)^2 Re_D}\right] u_j^n+ \frac{\Delta t}{2(\Delta y)^2 Re_D}\left( u_{j+1}^n + u_{j-1}^n \right) \]
可得：
\[ A = -\frac{E}{2}, \qquad B = 1 + E, \qquad K_j = (1 - E) u_j^n + \frac{E}{2}(u_{j+1}^n + u_{j-1}^n) \]
数值求解
我们使用Crank–Nicolson隐式方法
👉 由作业二可知方法：无条件稳定（unconditionally stable）
- 可以使用较大的时间步长 👉 \(E\) 可以取任意值
- 加快收敛到稳态?
- 若要准确模拟瞬态过程：👉 必须选择较小的 \(\Delta t\)
💻 代码演示
Couette.m
关于参数 \(E\) 的选取
当我们简单地增大 \(\Delta t\)（即增大 \(E\)）时：
- 达到稳态所需的时间步数减少
- 体现了隐式方法的优势
但当 \(\Delta t\) 足够大（例如 \(E > 400\)）时：👉 情况发生反转：
- 所需时间步数反而增加
- 收敛效率变差
- 隐式方法的实际优势消失
关于参数 \(E\) 的选取
- Crank–Nicolson 方法：
- 无条件稳定
- 适用于稳态问题
- 但必须权衡：
\[ \text{精度} \quad \text{vs} \quad \text{计算效率} \]
👉 实际应用中：选择一个“适中”的时间步长最为关键
重要结论
稳态解：
- 稳态解与 \(Re_D\) 无关
- 但瞬态过程依赖于 \(Re_D\)
非定常过程：
- 初始剖面逐渐演化
- 最终趋于线性分布
数值现象：
- 大时间步 → 非物理振荡
- 小时间步 → 精度高但计算慢
一个研究方向：时间精度更高的隐式方法
👉 目的在于：将隐式方法更有效地应用于非定常流动问题
作业：Poiseuille 流动的数值计算与格式对比 | 5月9日提交打印版报告
在课堂 Couette 流动算例的基础上，进一步研究平行平板间的 Poiseuille 流动。
Poiseuille 流动由压强梯度驱动，其速度剖面在稳态时为抛物线分布，是验证粘性流动数值方法的重要经典算例。
控制方程
考虑两无限长平行平板之间的不可压粘性流动，板间距离为 \(D\)。
令无量纲坐标为：
\[ 0 \le y \le 1 \]
无量纲速度满足：
\[ \frac{\partial u}{\partial t}=\frac{1}{Re_D}\frac{\partial^2 u}{\partial y^2}+ S \]
作业：Poiseuille 流动的数值计算与格式对比 | 5月9日提交打印版报告
无量纲速度满足：
\[ \frac{\partial u}{\partial t}=\frac{1}{Re_D}\frac{\partial^2 u}{\partial y^2}+ S \]
其中：
- \(u\) 为无量纲流向速度
- \(Re_D\) 为基于板间距的雷诺数
- \(S\) 为无量纲压力梯度驱动项
本作业中取：
\[ Re_D = 100, \qquad S = \frac{8}{Re_D} \]
作业：Poiseuille 流动的数值计算与格式对比 | 5月9日提交打印版报告
边界条件与初始条件
采用无滑移边界条件：
\[ u(0,t)=0 \]
\[ u(1,t)=0 \]
初始条件可取（可以灵活选择）：
\[ u(y,0)=0 \]
作业：Poiseuille 流动的数值计算与格式对比 | 5月9日提交打印版报告
稳态精确解
当流动达到稳态时：
\[ 0=\frac{1}{Re_D}\frac{d^2 u}{dy^2}+ S \]
结合边界条件：
\[ u(0)=0, \qquad u(1)=0 \]
可得稳态精确解：
\[ u_{\text{exact}}(y)=\frac{Re_D S}{2} y(1-y) \]
根据 \[ S = \frac{8}{Re_D} \Longrightarrow u_{\text{exact}}(y)=4y(1-y) \]
作业：Poiseuille 流动的数值计算与格式对比 | 5月9日提交打印版报告
计算格式要求
请采用以下三种时间推进格式求解非定常 Poiseuille 流动。
- FTCS 显式格式
\[ \frac{u_j^{n+1}-u_j^n}{\Delta t}=\frac{1}{Re_D} \frac{u_{j+1}^{n}-2u_j^n+u_{j-1}^{n}}{(\Delta y)^2}+ S \]
整理得：
\[ u_j^{n+1}=u_j^n+E\left(u_{j+1}^{n}-2u_j^n+u_{j-1}^{n}\right)+\Delta t S \]
其中：
\[ E=\frac{\Delta t}{Re_D(\Delta y)^2} \]
作业：Poiseuille 流动的数值计算与格式对比 | 5月9日提交打印版报告
计算格式要求
注意显式格式稳定性条件为：
\[ E \le \frac{1}{2} \]
- BTCS 隐式格式
\[ \frac{u_j^{n+1}-u_j^n}{\Delta t}=\frac{1}{Re_D}\frac{u_{j+1}^{n+1}-2u_j^{n+1}+u_{j-1}^{n+1}}{(\Delta y)^2}+ S \]
整理为：
\[ -E u_{j-1}^{n+1}+(1+2E)u_j^{n+1}-E u_{j+1}^{n+1}=u_j^n+\Delta t S \]
作业：Poiseuille 流动的数值计算与格式对比 | 5月9日提交打印版报告
计算格式要求
- Crank–Nicolson 格式
\[ \frac{u_j^{n+1}-u_j^n}{\Delta t}=\frac{1}{2Re_D}\left[\frac{u_{j+1}^{n+1}-2u_j^{n+1}+u_{j-1}^{n+1}}{(\Delta y)^2}+\frac{u_{j+1}^{n}-2u_j^n+u_{j-1}^{n}}{(\Delta y)^2}\right]+ S \]
整理为：
\[ -\frac{E}{2}u_{j-1}^{n+1}+(1+E)u_j^{n+1}-\frac{E}{2}u_{j+1}^{n+1}=\frac{E}{2}u_{j-1}^{n}+(1-E)u_j^n+\frac{E}{2}u_{j+1}^{n}+\Delta t S \]
作业：Poiseuille 流动的数值计算与格式对比 | 5月9日提交打印版报告
任务 1：三种格式的程序实现
编程语言不限，计算报告中应给出关键代码片段
任务 2：与精确解对比
对三种格式分别计算到接近稳态，并与精确解比较：
\[ u_{\text{exact}}(y)=4y(1-y) \]
至少绘制：
- 四个不同典型时刻速度剖面图
- 最终稳态数值解与精确解对比图
- 误差随时间步数变化图
误差可定义为：
\[ error_{\infty}=\max_j\left|u_j-u_{\text{exact},j}\right| \]
作业：Poiseuille 流动的数值计算与格式对比 | 5月9日提交打印版报告
任务 3：研究参数 \(E\) 的影响
研究不同 \(E\) 对计算结果的影响。
其中：
- FTCS 格式只需测试满足稳定性条件的情况
- BTCS 和 CN 格式可测试较大 \(E\)
- 对比不同 \(E\) 下的收敛速度、误差水平和是否出现非物理振荡
作业：Poiseuille 流动的数值计算与格式对比 | 5月9日提交打印版报告
计算报告要求
提交一份计算分析报告，内容包括：
- 程序实现说明
- 网格与参数设置
- 数值结果图像
- 参数 \(E\) 对结果的影响
- 三种格式优缺点总结