审查日期：2026-10-23。本文记录连续方程、数值实现风险和需通过的检查，不代替具体算例的运行结果。旧 compute-wake.py、wake-verification.md 已阅读；旧例为固定圆柱，不是自由涡振。
设实验室坐标为 (X,Y)，结构位移为 y(t)，采用 ξ=X, η=Y−y(t)。网格节点在实验室中随结构整体平移，单元形状和面积不变。这是移动网格的一种特殊情形，不需要局部变形或重划分网格。
令相对流速 w=(u_lab,v_lab−ydot)。对恒密度不可压流体，坐标变换直接给出：
div(w) = 0
∂w/∂t + (w·∇)w = −∇p/ρ + ν∇²w − yddot e_y
m yddot + c ydot + k y = F_L
f_n = sqrt(k/m)/(2π), c = 2 ζ m (2π f_n)这里的 p 仍是物理压力；圆柱表面 w=0，来流 w=(U,−ydot)。将求得的壁面升力代回结构，再用新的运动更新流体边界和惯性源，即构成双向流固耦合。结构自由度的位移并未由正弦函数指定。
Leontini 等的低 Re 圆柱研究已经使用附着于圆柱的加速参考系，在动量方程中加入参考系加速度，并讨论流固时间推进的稳定性。因此该坐标路线有直接研究先例。原论文；作者公开 PDF。
整体平移下速度梯度不变，所以涡量可以直接使用随体速度的旋度。可视化若选择实验室坐标，须同时移动圆柱、网格和流场位置；粒子在实验室坐标积分时须使用 v_lab=w_y+ydot。仅把圆柱图标上下移动、尾迹仍完全固定于实验室坐标，会给出错误的运动关系。
本例体积力为 F=(0,−ρ yddot)。若格子时间步为 1，碰撞前分布对应的物理速度应为：
ρ w = Σ_i f_i c_i + F/2
S_i = (1−1/(2τ)) W_i [3(c_i−w)·F + 9(c_i·w)(c_i·F)]碰撞平衡态和源项都要使用同一个 w。忽略半步速度修正只是在分布上加力，不能按标准 Guo 格式解释。Guo、Zheng、Shi 原论文。
旧求解器采用 pull-stream 后 collide 的顺序，其返回数组是碰撞后的分布。同一个碰撞时刻的速度若从该数组恢复，应使用 ρw=Σf_post c−F/2；不能直接沿用碰撞前的正号。也可以改成统一的 collide-stream 顺序，但所有边界、导出和诊断都必须同步。
入口施加的是物理 w_y=−ydot，所以边界重建所需的原始动量速度为 Σf c_y/ρ=−ydot+yddot/2。出口若采用密度条件，应外推随体横速并保持相同半步约定。旧出口隐含的 w_y=0 不能照搬，否则相当于在实验室中强制出口横速为结构速度。基础入口/出口方法来源见 Zou 与 He 原论文。最终运行版已将入口和出口改为非平衡外推，具体按下文的源码核验记录解释。
圆柱固定在计算坐标中、内部节点不作为流体推进时，壁面反弹的动量交换提供壁面受力。保留物理压力和 −ρ yddot 源项的上述方程，并不要求人为再加或减 ρ×圆柱面积×yddot。后者可能来自另一种修改压力变量的推导，或某些包含虚拟内部流体的浸入边界算法，不能跨算法套用。
同理，若右侧 F_L 是完整 CFD 流体力，则结构左侧 m 是结构质量，不能再额外加入理论附加质量而重复计入流体惯性。界面中标注的 f_n 应明确为这一质量定义下的结构固有频率；流体中的实际响应频率可随耦合发生偏移。
惯性系中的移动固壁动量交换还涉及边界速度、节点相变等附加项；本方案固体掩码不变，不能直接套用所有移动掩码修正。Giovacchini 与 Ortiz 的力推导可用于区分这些算法情形。
源项、结构速度、结构加速度和反弹力应标明各自的时间层。简单一次更新属于松耦合，不应标为每步迭代收敛的强耦合。低质量比时尤其容易产生附加质量数值不稳定；若需用低质量比，优先迭代流体与结构直至残差满足容差，或至少给出时间步敏感性证据。
- 对同一流体、同一直径，保持运动黏度不变，升风速会同时提高 Re=UD/ν。若改变黏度以固定 Re，则必须在说明中明确，这是参数研究而非同一介质连续升速。
- 固定体脱涡参考频率 f_s0≈St U/D 随 U 增大；耦合后实际升力频率 f_L 和结构响应频率 f_y 应从记录识别，不能事先把它们强制设置成 f_n。
- f_s0 越过 f_n 提供寻找大响应区的线索。是否出现锁定要由多个持速工况中 f_L、f_y 和 f_n 的关系判断；单个幅值峰值不足以证明完整锁定区。
- 分别从静止/小扰动起算的多工况是离散风速扫描；沿用前一工况状态的是连续历史。网页切换预计算工况不能被描述为即时求解了连续变风速过程。
- 侧向周期边界仍表示同步运动的周期圆柱阵列。旧例间距为 10D，不可直接声称无限远场孤立圆柱基准得到验证。
- 低 Re 二维圆柱适合演示“尾涡—力—位移—流场反馈”。它不代表高 Re 桥梁真实断面的定量抗风结果。
1. 参考系符号检查：去掉圆柱，用空间均匀流场与规定 yddot(t) 检查 w_y+ydot≈0；可以检出惯性源符号、半步速度和时序错误。
2. 固定体回归：将结构锁住，检查是否恢复相同边界、同 Re、同网格的固定圆柱结果。
3. 结构独立检查：零外力自由衰减识别频率和对数衰减率；检查离散能量变化与 ∫F_L ydot dt−∫c ydot²dt 相容。
4. 每工况持速与测频：去除启动过程；比较后两个窗口的幅值，采用频谱和过零周期相互核验。保存实际分析窗口、记录长度和频率分辨率。
5. 数值敏感性：至少对大响应工况进行更细网格或更低格子速度的对照；检查密度范围、全部有限值、峰值速度和实际局部 Mach 数。质量密度稳定并不等于已完成流固收敛验证。
6. 结果接受规则：数值发散、持续异常密度漂移或未趋稳的工况不能包装成稳定幅值曲线。只报告确实运行并通过检查的结果。
Python 3.13 中 NumPy、Numba、SciPy 可导入；现有课程 LBM 因而可直接扩展。MATLAB R2022a 安装目录及 PDE、Simulink 工具箱目录存在，但未核验许可证或可用求解能力。PATH 中未发现 OpenFOAM/SU2 可执行文件；WSL 命令报告尚未安装发行环境。Python 中未发现 FEniCS/DOLFIN/DOLFINx、Taichi、JAX、PyVista、Gmsh。已连接工具中没有可直接求解流体或流固耦合的专用 CFD 服务。
这些结果只说明此次浅层检查没有找到可即时复用的成熟 CFD 求解安装，并不证明整台机器没有其他软件。当前最短可复现路线是扩展已有 NumPy/Numba LBM；若升级为研究级对比，应另用成熟 ALE/动网格求解器和独立基准交叉验证。
审计脚本为 verify-fsi.py，结果为 qa/fsi-independent-method-checks.json。下面记录的通过版本 SHA-256 是 d489f5cab2eb4df19b86609e2137ec933f532d56f0a92758592da615030c6e03；后续修改若涉及算法，需要重新运行。脚本导入生产求解器，不修改其源文件或计算数据。
初版显式结构更新在零阻尼自由振动 20 周期后能量增加约 28%，设置 ζ=0.01 时识别到约 0.0090184；这是数值反阻尼，不应解释为流体向结构传能。初版源项/出口时间层不一致还使无圆柱均匀加速流出现横速偏差。这一版已拒收，初始检查结果单独保存为 qa/fsi-independent-method-checks-rejected-initial.json，不得混入正式展示数据。
- 结构在一步中保持上一时刻 CFD 力，使用线性结构的隐式中点法。取 γ=c/m, ω²=k/m，保存的端点加速度为 a=(F_old−cv−ky)/m，则 Δv=(a−ω²v/2)/(1+γ/2+ω²/4)；y_new=y+v+Δv/2，v_new=v+Δv。
- 随体网格使用实际速度增量，当前步惯性源为 a_frame=−Δv。上一步与当前步的碰撞后半步力约定通过分布的一阶矩修正转换，避免人为引入半步速度跳变。
- 入口非平衡外推施加当前结构速度的相反数；出口从上一步邻点分布外推，并转换半步力约定。恒定体积力下不能多减一次源项，否则出口会滞后一个速度增量。平衡与非平衡分布分离、邻点外推的基础方法见 Guo 等 2002 年边界论文；变参考系半步转换是本实现的额外处理，不能把全部程序声称为原论文的逐项复现。
- 壁面反弹力仍由前一碰撞状态的动量交换得到，作为下一结构步的载荷。这是显式分区、旧力保持一步的时间耦合，未实施每步流固迭代，不称为强耦合。
无圆柱检查使用真实 step_flow、48×8 网格和 U=0.04，覆盖入口、出口和内部节点。规定网格运动后，正确实验室横速应保持零。
| 检查 | 最大实验室横速误差 | 解释 |
常加速度 +2×10⁻⁶，600 步 | 1.06×10⁻¹⁵ | 接近舍入误差 |
常加速度 −2×10⁻⁶，600 步 | 1.30×10⁻¹⁵ | 接近舍入误差 |
正弦加速度幅值 2×10⁻⁵，周期 200 步，600 步 | 4.15×10⁻⁸ | 约为来流速度的 1.04×10⁻⁶；不是机器精度不变性 |
正弦加速度幅值 5×10⁻⁶，周期 400 步，1200 步 | 6.26×10⁻⁹ | 约为来流速度的 1.57×10⁻⁷ |
结构检查直接调用生产 advance，使用全固体单格使流体力恒为零，从而隔离实际结构推进器；不是只核验另一套仿写公式。D=16、24，各运行 20 个结构周期，零阻尼最终/初始能量比均在 1±8×10⁻¹⁵ 内。ζ=0.01 分别识别为 0.0099999871、0.0099999910；ζ=0.05 分别识别为 0.0499998726、0.0499999436。自由衰减能量闭合误差相对于初始能量不超过 3.86×10⁻⁶（能量耗散项使用梯形积分）。
这些检查支持源项符号、时间层转换与结构质量/刚度/阻尼实现的一致性。它们没有检验圆柱壁面精度、附加质量空间误差或整场流固网格收敛。正式数据的独立频率、位移二阶差分残差和逐周期功平衡由 qa/fsi-independent-audit.py 另行处理。
七个 D=16 工况均已完成 45 个结构固有周期，独立审计使用 25≤f_n t≤45 的导出记录；全文数据与逐周期结果保存在 qa/fsi-main-data-audit.json。表中频率为独立 Hann 频谱峰值的局部抛物线估计；其原始频谱分辨率约 0.04994 f_n，不应把插值后的多位小数视作测频精度。
| 工况 | 最后10周期 A/D | 位移主频/fn（约） | 四阶位移方程残差/气动力RMS |
| Ur=3.5，Re=70，ζ=0.01 | 0.00799 | 0.573 | 0.390% |
| Ur=4.5，Re=90，ζ=0.01 | 0.09409 | 0.815 | 0.223% |
| Ur=5.5，Re=110，ζ=0.01 | 0.37405 | 1.001 | 1.193% |
| Ur=6.5，Re=130，ζ=0.01 | 0.06964 | 1.183 | 0.452% |
| Ur=8.0，Re=160，ζ=0.01 | 0.04487 | 1.557 | 0.705% |
| Ur=5.5，Re=110，ζ=0.08 | 0.14311 | 1.017 | 0.396% |
结果支持“频率接近时响应放大，随后远离时响应减小；增加阻尼降低响应”的课堂结论。当前五个风速点不足以确定完整锁定区上下边界。Ur=6.5 最后两个10周期幅值仍约差 7.6%，所有结果都按有限时间窗口呈现，不把它们一概称为严格稳态振幅。
峰值工况升力与位移主频分别约为 1.0009 fn 与 1.0007 fn，位移正向过零频率约为 1.0005 fn。原始升力含高次成分、在一个主周期中可能多次过零，因此其原始过零计数给出约 1.336 fn，不可将这一数值作为主频或作为与 FFT 吻合的证据。低速和高阻尼工况也存在升力过零计数偏差，界面应以频谱主峰为主。
数据配对检查覆盖每个工况的全部48帧，帧时间严格递增，流场和状态元组数量一致。每组有3个帧时刻恰与高频时程采样重合；对应 CL, CD, y, v, a 差值均精确为零。非固定工况每帧的端点运动方程残差亦在浮点误差量级。这与源码中“同一步保存场和结构元组”的实现一致，未发现将其他时刻位移、升力拼接进流场的问题。
结构峰值工况的二阶中心差分残差为气动力 RMS 的约8.44%；此时惯性项与刚度项接近抵消，二阶导数截断误差会被较小的气动力尺度放大。独立四阶差分后降为约1.19%；相对于惯性、阻尼、刚度和气动力各项 RMS 之和则约0.0158%。未用保存的加速度恒等式替代这一独立检查。
必须区分同刻物理功与实际算法旧力保持时序的功。直接用同刻 F_L(t)v(t) 计算窗口收支，Ur=8 工况相对于阻尼耗散的闭合偏差约13.9%；这项误差没有被删除。按隐式中点结构实际使用的上一步保持载荷，其等效连续载荷为 F_L(t−Δt/2)+O(Δt²)。对原始记录做三次样条重建、细化积分，六组非固定工况的对应算法收支偏差均不超过阻尼耗散的0.0211%。这是对实现时间层的一致性检查，不能据此说同刻物理传能已经精确或一阶耦合误差不存在。
每组检查保存有限值、密度和局部 Mach 数、等间隔记录、帧配对、端点方程、独立位移差分与算法能量闭合。七个主组均满足明确的教学数据一致性阈值；细网格及最终验收见下节。
Ur=5.5、Re=110、ζ=0.01 的 D=20 细网格（400×200）也已完成45周期。所有原始结果经 qa/fsi-independent-audit.py 统一复核，最终 qa/fsi-independent-audit.json 顶层为 passed: true，所需8组工况无缺项；qa/fsi-main-data-audit.json 是为总报告兼容保留的相同内容副本。
| 对照项目 | D=16，320×160 | D=20，400×200 |
| 最后10周期 A/D | 0.374045 | 0.410281 |
| 最后20周期 A/D | 0.374095 | 0.410752 |
| 位移频谱峰值插值/fn | 1.000654 | 0.998947 |
| 位移过零主频/fn | 1.000500 | 0.998991 |
| 四阶位移方程残差/气动力RMS | 1.193% | 1.456% |
最后10周期幅值相对细网格仍差8.83%（最后20周期为8.92%），位移频谱峰值插值相对差0.171%；两个网格的原始 FFT 主峰位于同一频率格。结论是主频较稳定，幅值仍有网格敏感性。单个细网格点、相同格子来流速度的对照同时改变空间分辨率和物理时间步，不能分别证明空间和时间已经收敛。
细网格的独立四阶残差相对于各力项 RMS 之和为0.0151%，按算法实际旧力时序重建的窗口能量收支偏差为阻尼耗散的0.000708%。48帧的数量、顺序和端点方程检查通过；其中2个帧时刻与时程采样重合，全部状态元组值完全相等。
最终 passed 的含义限于已明确标注局限的教学回放验收：所有数组有限、平均密度偏差小于1%、诊断局部 Mach 小于0.3、时程等间隔、帧状态配对精确、独立四阶差分残差低于气动力RMS的3%且低于各力项RMS之和的1%、实际算法时序收支残差低于阻尼耗散的1%。完整阈值逐项结果和输入文件 SHA-256 随 JSON 保存。这不是研究级孤立圆柱基准验证，也不是强耦合物理传能或幅值完全收敛的证明。
浏览器只读取 fsi-data.js，风速滑块选择已计算工况，未调用在线求解器。求解器数组及固体掩码固定于随体坐标；整体平移是映射回实验室坐标的解释，不应在教学界面混同外边界固定的局部 ALE 网格变形。
可视化曾用间距 0.5D、以圆柱中心为原点的线叠加在涡量图上。实际圆柱中心为 (4D,5.07D)，这些线不保证与原始均匀格点重合，所以原“抽稀网格”表述不准确。现已改称随体坐标参考线，默认关闭，并明确本页是预计算回放。