桥梁抗风 · 进阶实验室01 / 计算流体力学

CFD · 从时间采样理解离散化

导数，就是相邻记录的变化率

先用位移、速度和加速度理解原函数、一阶导数、二阶导数。拖动时间步 Δt，看连续曲线怎样变成计算机里的数值；空间导数与有限体积实验可在下方展开。

先看时间，再理解 CFD

每隔 Δt 记录一次：位移 → 速度 → 加速度

二维时间采样 · 可调步长

把桥梁某一点的位移记作 f(t)。计算机先存下一个个时刻的数值，再用相邻数值相减、除以时间步来近似导数。

实线：原函数○：采样值蓝 / 橙：前后两段平均速度灰虚线：精确导数

原函数 f(t) 位移 / m

一阶导数 f′(t) 速度 / (m/s)

二阶导数 f″(t) 加速度 / (m/s²)

三张图共用时间轴。平均速度画在各自的时间区间内；二阶差分画在中间时刻 ti。

01 · 原函数

先记录每一帧的位移

fi = f(ti)

下标 i 只是第几个采样点；相邻时刻相隔相同的 Δt。

02 · 一阶导数

相邻位移相减，除以 Δt

f′(ti) ≈ (fi+1 − fi) / Δt

右式是 [ti, ti+1] 内的平均速度 v̄+，用来近似 ti 的瞬时速度；这是前向一阶差分。

03 · 二阶导数

两段平均速度相减，再除以 Δt

f″(ti) ≈ (fi+1 − 2fi + fi−1) / Δt²

前段 v̄− = (fi−fi−1)/Δt；两段区间中点也相隔 Δt，所以 (v̄+−v̄−)/Δt 就得到上式。

公式条件：屏幕读数保留三位小数，计算使用未舍入的值。以上采用均匀时间步，二阶式使用 ti 前后各一点，是中心二阶差分。函数足够光滑时，前向一阶差分的截断误差为 O(Δt)，中心二阶差分为 O(Δt²)。匀加速例子的二阶差分恰好精确；切换“振动”再缩小 Δt，更容易看出近似误差。

和 CFD 有什么关系？

固定流场里的一个位置，每隔 Δt 记录当地风速 u，就能用 [u(t+Δt)−u(t)]/Δt 近似 ∂u/∂t。若 u 本身是风速，它的时间导数单位就是 m/s²；本实验先用位移，是为了借助熟悉的“速度、加速度”理解同一种差分运算。

以后只需把“相邻时刻、Δt”换成“同一时刻的相邻位置、Δx”，就得到空间差分的类似形式。理解时间差分后，再展开空间导数。

进阶：展开空间导数与空间—时间曲面（保留原实验）
固定 t，沿 x 读取同一条曲线
W
P
E
x / mφ
同一时刻 t = 0.600 s · 三点等距 Δx = 0.600 m
蓝：左割线橙：右割线灰：P 点切线（一阶导）
拖动旋转 · 滚轮缩放 · R 复位
W
P
E
STEP 1 / 4

先在同一个时刻，读出 W、P、E

把割线放到二维坐标中看清楚

先得到两段一阶斜率；它们所在的两个半点之间，相距 Δx。

只缩小间距，误差怎样变化？

真实流场有 x、y、z、t 四个自变量。 本页取其中一个空间方向，把 φ(x,t) 画成“x—时间—数值”曲面；竖直高度是 φ，不是第三个空间坐标。固定 t 看空间变化；固定 x 看时间变化。

导数公式依据与制作说明

NIST DLMF：三点二阶偏导与O(Δx²)精度；UCL：前向与中心差分的Taylor推导。本课使用与制作说明；原程序迁移与数学核查记录。

走向 CFD：为什么总会看到“右边减左边”？

在均匀网格、常扩散系数 Γ 下，右、左面梯度分别近似为上面的两段斜率。因此 Γ(me−mw)/Δx 恰好变成 Γ(φE−2φP+φW)/Δx²。有限差分直接近似导数；有限体积先积分守恒、再近似面通量，此时得到相同的内部模板。

展开后段实验：把导数装进守恒方程，真正求解一维对流扩散
一个单元，两个有效通量面
一维守恒模型 · 横向挤出示意CELL 10 / 20
x → 传输方向 · 拖动旋转 · 滚轮缩放 · 单击选择单元
左、右面共享同一组通量；横向无梯度，另外四个面的净通量为零。
对流 扩散 当前控制体
网格 PeΔ0.55PeL / N
中心点 RMS 误差—相对解析解的绝对误差
最大单元守恒残差—|Jₑ − Jw − SΔx|
数值解有界性—在边界值 [0, 1] 内？

数值解与解析解

同一组边界条件，所有网格均保留源项

橙色为中心格式，青色为迎风格式，深蓝为解析解；竖线标记选中单元。

网格收敛：N = 10 → 160

固定同一个物理问题，只改变网格

二阶中心格式趋近斜率 −2；迎风格式趋近 −1。未进入渐近区时，不能只凭两组网格判定阶次。

SELECTED CONTROL VOLUME

第 10 个控制体

J = F φf − Γ ∂φ/∂x
Je − Jw = SP Δx

项（正 x 为正）西面 w东面 e
对流 Fφf
扩散 −Γ∂φ/∂x
总通量 J

FROM FLUX TO MATRIX

这一行方程是怎样来的？

单元中心 数值 φP 解析 φ(xP)
展开最小求解过程
// 正 x 方向总通量；每个内部面只计算一次
J[f] = F * phiFace[f] - Gamma * gradient[f];
// 单元积分守恒：右出减左入 = 体源项
A * phi = sourceIntegral;
// 边界在面上，中心到边界的距离为 dx / 2
boundaryGradient = (phiP - phiBoundary) / (dx / 2);

网页的 physics.js 直接组装三对角线性方程，并独立重算通量与残差。小残差表示离散方程解得准确，不能代替网格误差检验。

读图之前

三件容易混淆的事

PeΔ 不是 Re。 本页 PeΔ = FΔx/Γ 衡量网格尺度上的对流与扩散；只有输运动量且 Γ = μ 时，才有 PeΔ = ReL Δx/L。

数值振荡不是涡脱落。 中心格式在 |PeΔ| > 2 时失去正系数保证；曲线越界可以来自离散误差，不能据此判断真实流动进入湍流。

有界不等于准确。 一阶迎风会增加数值扩散。应同时检查守恒、解析或基准比较、网格与时间步敏感性。

课堂任务 · 约 10 分钟

先预测，再操作，再解释

  1. 点击“触发数值振荡”。解为何越过边界的 [0,1]？检查 PeΔ 和邻点系数。
  2. 在相同网格上改用迎风。曲线变平滑后，RMS 误差一定更小吗？
  3. 保持 PeL = 11，逐次加密网格。两种格式误差的衰减速度有何区别？
  4. 切换源项问题。正的 φ″、零端值应给出什么符号的 φ？
源程序、纠错说明与阅读依据

由 wolai 抗风课程页中 cfd.mlx、cfd_fvm.mlx 的教学内容迁移。保留控制体与对流扩散主线，重新核对了源项符号、边界位置、解析解和网格收敛；详见 源程序审阅记录 与 制作说明。

本课不是桥梁绕流求解器。三维控制体是准确的一维模型挤出示意；桥梁计算域、网格和箭头只说明建模选择，未生成压力云图、流线、气动力系数或工程验算结果。