替友人 A 写一份气动弹性复现脚本：一次闲时 vibe coding 记录气动弹性复现脚本

这不是我自己的毕业课题。一次闲聊里，友人 A 提到他在看 2024 年西北工业大学的硕士论文《马赫反射二维壁板气动弹性理论分析》，想把图 3-11、3-12 和 5-7 重新画出来。我对这个方向并不熟，却正好想试试：一个非本领域的人，借助 vibe coding，能不能把论文中的公式整理成一份有用的复现脚本？

所以这个小项目从一开始就不是“我证明了论文结论”，而是“我替友人 A 做了一套可运行的计算草图”。A 在本文中匿名；代码由我在闲时组装、运行和检查，但大量模型翻译与绘图实现借助了生成式工具。这个前提很重要，因为代码能跑不等于科学上可信。

我先做的不是写界面，而是问 A 到底想看什么

最初的需求很像一句普通的 vibe coding 指令：把论文里的几张图做成可调参的 Python 程序，最好有 GUI，能导出图片。但真正开始时，问题很快从“怎么画”变成“凭什么算”：

• 论文中的无量纲参数是否全部给出？
• 马赫杆位置 ξT 如何进入两段气动力积分？
• 图线是静力平衡、线性稳定性，还是长时积分后的极限环？
• 原文没有写清的质量比、结构阻尼和面内力，脚本应该怎么处理？

我们最后把需求收缩为三件事：公式能被程序化，关键参数可显式调整，失败时不伪装成一条有效曲线。后来看，第三件做得还不够好，却也成了这次练习最有价值的部分。

和模型对话时，我必须自己守住方程的形状

项目从无量纲二维壁板方程出发，包含弯曲刚度、面内力、结构阻尼、惯性、气动力和静压差。马赫杆位置 ξT 把壁板分成两段：上游超声速区用线性活塞理论，下游亚声速区用压缩性修正势流。气动力在代码里被拆成 A1–A8，分别承担准定常项、气动阻尼和亚声速附加项。

位移用 Galerkin 展开：W(ξ, τ) = Σ qᵢ(τ) sin(iπξ)。连续偏微分方程因此变成有限个模态坐标 qᵢ 的非线性常微分方程。von Kármán 几何非线性又把所有模态的能量耦合回每个方程，所以“用两个正弦项近似”不代表后面就只剩线性叠加。

生成式工具很擅长把符号快速转成 NumPy 数组操作，也很容易在下标、分段积分或符号约定上“顺手补全”。我的工作不是逐字接受它的实现，而是不断把程序中的矩阵形状、对称性、量纲和极限情形拉回来核对。

三条计算路径，三种不同的失败

静力 Newton：没收敛就应该说没收敛

static_solve 用 Newton 迭代解模态残差，并提供显式 Jacobian 结构。它从较低气动力或接近零的初值出发，用于不同静压差下的静变形与静发散曲线。核心层会对收敛失败抛出异常；但绘图层曾将部分异常换成零数组。这个方便 GUI 不崩溃的处理，会让读者误以为“变形正好为零”。

时域积分：图画出来不代表积分已收敛

动力状态写成 [q, q̇]。核心模块实现固定步长四阶 Runge–Kutta，图 3-11 的极限环幅值则使用 SciPy solve_ivp 自适应积分。图 5-7 需要逐幅续算，脚本保留固定步长 RK4，并在数值发散时最多四次减半步长。这是一种实用保护，不是收敛性证明。

稳定性扫描：跨过零点只是数值判据

stability_matrix 在静平衡附近线性化系统，计算特征值实部。程序扫描无量纲动压 λ，找最大实部由负转正的区间，再线性插值得到颤振临界动压。这条路径在数值上很清晰，但仍依赖前面的平衡解、参数选择和模态截断。

我把第一批图发给 A 时，更像是发了一份检查清单

本文封面是直接调用项目 make_fig_311 生成的结果，不是概念插画。为控制计算时间，当时使用 N=2、C=0.01、λ=500，对 μ=0.1 / 0.5 / 1.0 和七个静压差点扫描。在这组输入下，临界动压对零静压差近似呈 U 形，质量比越大，临界值越高。

我发给 A 的同时也包括几个保留项：论文没有完整报告 μ、C、Rx 等输入；极限环幅值量级很小；当时还没做步长与模态阶数的收敛检查。因而它只能证明脚本在该参数组下产生了这些曲线，不能单独证明与论文严格一致。

图 3-12 教我警惕“很干净的零线”

六个子图比较整体静压差为 -100、-50、0、忽略、50、100 时的壁板挠度。黑线是较低动压下的静变形，蓝线是输入动压下的静发散解。正静压差工况中某些蓝线干净地落在零上，初看很像有明确物理意义；追到程序才发现，它可能是 Newton 失败后的零数组占位。

这个小错误很具有 vibe coding 特征：工具倾向于保证界面完整，而科学计算要求保证失败可见。更好的做法是把该工况标成 NOT CONVERGED，保留残差和迭代次数，而不是让零线伪装成结果。

图 5-7 没有被我们“调到好看”

振荡马赫杆让 ξT 随时间变化。脚本对振荡幅值 A 逐点推进，从上一稳态末端继续计算，丢弃瞬态后提取 ξ=0.75 处挠度的局部极值，再用 Poincaré 采样和极值聚类估计周期数。

论文工况 λ=500 属深度静发散，系统刚性强。高阶模态 N=4 很容数值发散，所以 GUI 默认改用 N=2，并允许缩短推进时间、减小步长或增大阻尼。我没有为了得到更像原图的分岔纹理而不断“调到好看”，因为那会把数值拟合冒充为物理复现。要讨论真实分岔，至少还需要步长收敛、模态截断收敛和不同积分器交叉验证。

作为一份闲时脚本，它已经有用；作为复现工程，它还没闭环

项目将 aeroelastic_panel.py、aeroelastic_figs.py 和 aeroelastic_gui.py 分开，核心函数基本不依赖 GUI，参数集中在字典中，这使得 A 可以快速调整条件、批量出图。从“替朋友做一份能用的工具”来说，它已经达到了最初目标。

但若要把它升级成可引用的复现工程，还有一份很具体的清单：

• 把表 3-2 和式 2-32 的对照结果写成可在 CI 重放的测试，而不只写在 GUI 文案中。
• 为每个参数记录单位、论文页码、公式号和“原文缺失”状态。
• 每张图同时保存参数 JSON、代码版本、收敛标志和积分统计。
• 将依赖版本、字体和导出环境锁定，避免换机器后图形或数值悄然变化。
• 删除或实现仍为空的 _cross_terms，不让函数名暗示已有的物理能力。

这次 vibe coding 给我留下的真正产物

最终产物不只是一个可打包的 Windows 窗口，也不是三张颜色很像论文的曲线。它更像一次边界训练：我可以借助模型快速进入一个陌生领域、拼出可运行工具，但必须对每个没有来源的参数、每次未收敛、每一条过于干净的曲线保持怀疑。

我替友人 A 写下的是一个能继续追问的起点，不是一张“复现完成”的证书。