CFD代理模型实战——翼型气动代理模型从采样到ONNX部署的全流程

从零搭建翼型气动代理模型——零CFD依赖即可完成90%开发。本文用解析函数做真值、Prandtl框架做训练、ONNX做部署，让你看到代理模型从采样到部署的完整链路。

一、目标：五步完成翼型代理模型

本文以薄翼型升力系数预测为例，展示CFD代理模型的完整开发流程：

1. 问题定义：输入参数化、输出定义
2. 采样设计：在输入空间选择训练点
3. 真值计算：用解析公式（替代CFD）生成标签
4. 模型训练：GP和MLP对比
5. 模型部署：导出ONNX用于生产

关键前提：开发代理模型框架时不需要任何CFD数据。解析函数做真值是验证框架逻辑最快、最可靠的方式——等框架验证通过后再对接真实CFD也不迟。

二、环境准备

pip install prandtl-cfd numpy matplotlib scikit-learn onnx onnxruntime

确认安装：

import prandtl as pr
print(f"Prandtl version: {pr.__version__}")

三、步骤1：问题定义

薄翼型升力公式

根据薄翼型理论，升力系数 $$C_L$$ 由迎角和弯度决定：

$$C_L = 2\pi(\alpha + 2h)$$

其中：

• $$\alpha$$：迎角（弧度）
• $$h$$：相对弯度（无量纲，如 NACA 2412 的 $$h = 0.02$$）

这个公式是势流理论的精确结果，适用于小迎角附着流。我们把它当作”真值函数”——在真实项目中这会是CFD求解器。

参数	符号	范围	物理意义
迎角	$$\alpha$$	[-5°, 15°]	来流与弦线夹角
相对弯度	$$h$$	[0.01, 0.10]	NACA四位数第2位（除以100）

四、步骤2：采样设计

用拉丁超立方采样（LHS）在二维输入空间中选择训练点：

from prandtl.sampling import lhs_sampling

# 100个训练点（纯解析函数，0秒生成）
bounds = [(-5.0, 15.0), (0.01, 0.10)]
X_train = lhs_sampling(bounds, n=100)

# 50个测试点（独立采样，评估泛化能力）
X_test = lhs_sampling(bounds, n=50)

print(f"训练集: {X_train.shape}")  # (100, 2)
print(f"测试集: {X_test.shape}")   # (50, 2)

LHS确保每个维度的一维投影都是均匀的，避免随机采样的聚集问题。

五、步骤3：生成真值

import numpy as np

def cl_true(alpha_deg, camber_pct):
    alpha = np.radians(alpha_deg)
    return 2.0 * np.pi * (alpha + 2.0 * camber_pct)

Y_train = cl_true(X_train[:, 0], X_train[:, 1]).reshape(-1, 1)
Y_test = cl_true(X_test[:, 0], X_test[:, 1]).reshape(-1, 1)

print(f"Y_train 范围: [{Y_train.min():.3f}, {Y_train.max():.3f}]")

在真实项目中，这里会调用CFD求解器。但从框架开发的角度，解析函数更快、更可控、支持无限次迭代。

六、步骤4：模型训练与对比

6.1 高斯过程回归（GP）

gp_model = pr.Surrogate(
    params=["alpha", "camber"],
    outputs=["CL"],
    method="gp"
)
gp_model.fit(X_train, Y_train)

gp_report = gp_model.validate(X_test, Y_test)
print(f"GP R-squared: {gp_report['CL']['r2']:.6f}")
print(f"GP RMSE: {gp_report['CL']['rmse']:.6f}")

预期结果：对线性函数，GP的R-squared应该 > 0.999。不是因为GP强大——是任务本身简单。

6.2 多层感知器（MLP）+ 物理约束

from prandtl.physics import Monotonicity, BoundaryValue

mlp_model = pr.Surrogate(
    params=["alpha", "camber"],
    outputs=["CL"],
    method="mlp"
)
mlp_model.fit(X_train, Y_train, physics=[
    Monotonicity(param_idx=0, sign=1, weight=0.1),
    # CL随迎角单调递增

    BoundaryValue({"alpha": 0.0}, {"CL": 0.0}, weight=10.0)
    # 对称翼型零迎角时升力为零
])

mlp_report = mlp_model.validate(X_test, Y_test)
print(f"MLP R-squared: {mlp_report['CL']['r2']:.6f}")
print(f"MLP RMSE: {mlp_report['CL']['rmse']:.6f}")

物理约束在这个简单线性函数上作用不明显——但在真实CFD中、数据稀少的区域，物理约束是防止模型做出物理上不可能预测的关键。

6.3 结果对比

print("方法      | R-squared | RMSE")
print("----------|-----------|--------")
print(f"GP        | {gp_report['CL']['r2']:.6f} | {gp_report['CL']['rmse']:.6f}")
print(f"MLP+物理  | {mlp_report['CL']['r2']:.6f} | {mlp_report['CL']['rmse']:.6f}")

两者都精确拟合了这个线性函数——证明框架核心逻辑正确。

七、步骤5：模型部署（ONNX导出）

MLP可以导出为ONNX格式——一个跨框架、跨平台的推理标准。GP不能导出（非参数模型，推理依赖全量训练数据）。

# 导出MLP为ONNX
mlp_model.export("cl_surrogate.onnx")
print(f"ONNX model saved: cl_surrogate.onnx")

加载并推理

import onnxruntime as ort
import numpy as np

# 加载ONNX模型
session = ort.InferenceSession("cl_surrogate.onnx")

# 推理：预测迎角5度、弯度2%时的升力系数
input_data = np.array([[5.0, 0.02]], dtype=np.float32)
cl_pred = session.run(None, {"input": input_data})[0]
cl_true_val = cl_true(5.0, 0.02)

print(f"ONNX prediction: {cl_pred[0, 0]:.4f}")
print(f"Ground truth:     {cl_true_val:.4f}")
print(f"Error:            {abs(cl_pred[0, 0] - cl_true_val):.6f}")

ONNX的优势

• 推理时不需要Python：可以用C++、C#、Java加载
• 零依赖：只需要onnxruntime，3MB的库
• 实时能力：单次推理 < 0.1ms，满足飞控400Hz要求
• 工业标准：所有主流ML框架都支持导出ONNX

八、从解析函数到真实CFD

框架验证通过后，对接真实CFD只需更换真值函数：

# 解析函数版（开发验证用）
Y = analytical_function(X)

# 替换为 OpenFOAM版（生产用）
import subprocess

def run_openfoam(alpha, camber):
    """生成翼型网格 -> 运行simpleFoam -> 提取CL"""

    # 1. 用参数化工具生成几何
    # 2. 自动网格划分（blockMesh/snappyHexMesh）
    # 3. 运行 simpleFoam
    # 4. 解析 forceCoeffs 输出
    pass

# 200次CFD计算，32核并行约31小时
from concurrent.futures import ProcessPoolExecutor
with ProcessPoolExecutor(max_workers=32) as executor:
    Y_train = list(executor.map(run_openfoam,
                   X_train[:, 0], X_train[:, 1]))

九、代理模型部署到实时仿真

训练好的ONNX模型最终部署到仿真器中。Pegasus Simulator是目前唯一可行的开源方案——它提供PX4 SITL桥接和PhysX集成，气动模块通过策略模式替换为你的代理模型：

全程推理 < 0.1ms，而等效的CFD求解需要数秒——差距超过四个数量级。

十、从翼型到旋翼飞行器

翼型代理模型是最简单的起点。对于完整的旋翼飞行器（如无人机），输入空间扩展到7-9维：

参数	符号	物理含义
马赫数	$$M$$	飞行速度
迎角	$$\alpha$$	机身迎角
侧滑角	$$\beta$$	机身侧滑
基准转速	$$\Omega_{\text{ref}}$$	旋翼转速
俯仰差动	$$\Delta\Omega_{\text{pitch}}$$	前后旋翼差速
滚转差动	$$\Delta\Omega_{\text{roll}}$$	左右旋翼差速
偏航差动	$$\Delta\Omega_{\text{yaw}}$$	对角旋翼差速

物理分解策略（最有效的减数据技巧）：

$$F_{\text{total}} = F_{\text{isolated rotors}} + F_{\text{fuselage}} + \Delta F_{\text{interference}}$$

孤立旋翼力用叶素动量理论（BEMT）解析计算，机身力用细长体理论——只让代理模型学习干扰修正项 $$\Delta F_{\text{interference}}$$。这通常占总力的10-30%，意味着预测任务更容易、数据需求更低。

十一、总结

阶段	内容	CFD次数	耗时
开发验证	解析函数 -> 训练 -> 验证	0	秒级
框架集成	单次CFD冒烟测试	1-2	小时级
批量训练	200次CFD并行计算	200	1-2天
部署	导出ONNX -> 嵌入仿真器	0	分钟级

核心理念：不要在CFD层调试框架——在解析函数层完成所有验证，确认代理模型的拟合能力没问题后，再花CFD的真金白银。

参考文献

1. Prandtl: CFD代理模型框架. https://github.com/goodisok/prandtl
2. ONNX Runtime. https://onnxruntime.ai/
3. McKay, M. D., Beckman, R. J., & Conover, W. J. (1979). A Comparison of Three Methods for Selecting Values of Input Variables in the Analysis of Output from a Computer Code. Technometrics, 21(2), 239-245. DOI: 10.1080/00401706.1979.10489755
4. Sobol, I. M. (1967). On the Distribution of Points in a Cube and the Approximate Evaluation of Integrals. USSR Computational Mathematics and Mathematical Physics, 7(4), 86-112.
5. Leishman, J. G. (2006). Principles of Helicopter Aerodynamics. Cambridge University Press.
6. Pegasus Simulator: 开源PX4-IsaacSim桥接. https://github.com/PegasusSimulator/PegasusSimulator