📚 课程目录第2章气动系数插值与气动力矩计算

2. 气动系数插值与气动力矩计算

任务定义

本案例实现导弹气动系数随马赫数（Ma）和攻角（α）变化的二维插值函数，并基于插值结果计算导弹所受的气动力（升力、阻力、侧力）与气动力矩（俯仰、偏航、滚转）。输入为当前飞行状态（Ma、α、β、动压q、参考面积S、参考长度l），输出为体轴系下的气动力与力矩分量。

气动数据表结构

采用工程中常用的稀疏网格存储方式：每个气动系数（升力系数C_L、阻力系数C_D、侧力系数C_Y、俯仰力矩系数C_m、偏航力矩系数C_n、滚转力矩系数C_l）均以Ma为行索引、α为列索引的二维表存储。插值采用双线性插值，保证系数曲面连续且计算量小，适合实时仿真。

代码实现

import numpy as np
from scipy.interpolate import RegularGridInterpolator

class AeroCoeffTable:
    """气动系数二维插值表（Ma × α）"""
    def __init__(self, ma_grid, alpha_grid, coeff_data):
        """
        ma_grid: 1D数组，马赫数断点（升序）
        alpha_grid: 1D数组，攻角断点（度，升序）
        coeff_data: dict，键为系数名，值为2D数组 shape=(len(ma), len(alpha))
        """
        self.ma_grid = ma_grid
        self.alpha_grid = alpha_grid
        self._interp = {}
        for name, data in coeff_data.items():
            self._interp[name] = RegularGridInterpolator(
                (ma_grid, alpha_grid), data,
                bounds_error=False, fill_value=None
            )

    def get(self, name, ma, alpha_deg):
        """查询单个系数，ma和alpha_deg可为标量或数组"""
        pts = np.column_stack([np.atleast_1d(ma), np.atleast_1d(alpha_deg)])
        return self._interp[name](pts)

def compute_aero_forces_moments(aero, ma, alpha_deg, beta_deg, q, S, l):
    """
    计算体轴系气动力与力矩
    输入:
        aero: AeroCoeffTable实例
        ma: 马赫数
        alpha_deg: 攻角（度）
        beta_deg: 侧滑角（度）
        q: 动压 (Pa)
        S: 参考面积 (m^2)
        l: 参考长度 (m)
    返回:
        dict: 包含力(N)与力矩(N·m)的体轴系分量
    """
    # 插值获取基本气动系数（侧力/偏航/滚转系数随β近似线性，此处用β=0查表+线性修正）
    CL = aero.get('CL', ma, alpha_deg).item()
    CD = aero.get('CD', ma, alpha_deg).item()
    Cm = aero.get('Cm', ma, alpha_deg).item()

    # 侧向系数：用β=0查表值乘以β(rad)的线性比例（简化模型）
    beta_rad = np.radians(beta_deg)
    CY_beta = aero.get('CY_beta', ma, alpha_deg).item()  # 侧力导数 /rad
    Cn_beta = aero.get('Cn_beta', ma, alpha_deg).item()  # 偏航力矩导数 /rad
    Cl_beta = aero.get('Cl_beta', ma, alpha_deg).item()  # 滚转力矩导数 /rad

    CY = CY_beta * beta_rad
    Cn = Cn_beta * beta_rad
    Cl = Cl_beta * beta_rad

    # 气动力（体轴系：X前，Y右，Z下）
    Fx = -CD * q * S          # 阻力沿-X
    Fy = CY * q * S           # 侧力沿+Y
    Fz = -CL * q * S          # 升力沿-Z

    # 气动力矩（体轴系）
    Mx = Cl * q * S * l       # 滚转
    My = Cm * q * S * l       # 俯仰
    Mz = Cn * q * S * l       # 偏航

    return {
        'Fx_N': Fx, 'Fy_N': Fy, 'Fz_N': Fz,
        'Mx_Nm': Mx, 'My_Nm': My, 'Mz_Nm': Mz,
        'CL': CL, 'CD': CD, 'Cm': Cm,
        'CY': CY, 'Cn': Cn, 'Cl': Cl
    }

# ========== 示例数据与调用 ==========
if __name__ == '__main__':
    # 定义马赫数与攻角网格（度）
    ma_pts = np.array([0.3, 0.6, 0.9, 1.2, 1.6, 2.0])
    alpha_pts = np.array([-4, 0, 4, 8, 12, 16, 20])

    # 生成示例气动数据（实际工程中由CFD或风洞试验提供）
    # 此处用解析函数生成平滑数据用于演示插值逻辑
    A_mesh, AL_mesh = np.meshgrid(ma_pts, alpha_pts, indexing='ij')
    # 升力系数：随α近似线性，随Ma在跨声速区有波动
    CL_data = (0.08 * AL_mesh + 0.1) * (1 + 0.15 * np.sin(3 * A_mesh))
    # 阻力系数：随α二次增长，跨声速区阻力增大
    CD_data = 0.02 + 0.0004 * AL_mesh**2 + 0.05 * np.exp(-((A_mesh - 1.0)**2) / 0.1)
    # 俯仰力矩系数：静稳定导弹Cm随α负斜率
    Cm_data = -0.01 * AL_mesh - 0.02 * (A_mesh - 1.5)
    # 侧向导数（单位/rad）
    CY_beta_data = np.full_like(CL_data, -0.3)   # 简化：常数
    Cn_beta_data = np.full_like(CL_data, 0.08)   # 简化：常数
    Cl_beta_data = np.full_like(CL_data, -0.05)  # 简化：常数

    coeff_dict = {
        'CL': CL_data, 'CD': CD_data, 'Cm': Cm_data,
        'CY_beta': CY_beta_data, 'Cn_beta': Cn_beta_data, 'Cl_beta': Cl_beta_data
    }

    aero_table = AeroCoeffTable(ma_pts, alpha_pts, coeff_dict)

    # 测试插值：Ma=0.8, α=5°
    test_ma = 0.8
    test_alpha = 5.0
    print(f"插值测试: Ma={test_ma}, α={test_alpha}°")
    print(f"  CL = {aero_table.get('CL', test_ma, test_alpha).item():.4f}")
    print(f"  CD = {aero_table.get('CD', test_ma, test_alpha).item():.4f}")
    print(f"  Cm = {aero_table.get('Cm', test_ma, test_alpha).item():.4f}")

    # 计算气动力与力矩
    result = compute_aero_forces_moments(
        aero_table, ma=0.8, alpha_deg=5.0, beta_deg=2.0,
        q=45000.0, S=0.42, l=0.75
    )
    print("\n气动力与力矩计算结果:")
    for k, v in result.items():
        print(f"  {k} = {v:.4f}")

数据流与插值逻辑

运行与验证

依赖安装：需要NumPy和SciPy库，执行 pip install numpy scipy。

运行方式：直接执行脚本 python aero_interp.py，程序会输出插值测试点（Ma=0.8, α=5°）的系数值以及对应气动力/力矩计算结果。

验证方法：

插值正确性：将查询点改为网格节点（如Ma=0.9, α=8°），返回值应与表中该节点数值完全一致（因双线性插值在节点处精确通过）。
量纲验证：力单位为N（q·S = 45000×0.42 = 18900 N），力矩单位为N·m（再乘以l=0.75）。例如CL=0.5时，Fz = -0.5×18900 = -9450 N。
趋势验证：固定Ma=0.8，增大α从0°到10°，CL应近似线性增大，CD应呈二次增长，Cm应负向增大（静稳定导弹）。

预期输出（示例）：当Ma=0.8、α=5°、β=2°时，插值得到的CL约为0.512，CD约为0.041，Cm约为-0.072；计算得到Fz约-9690 N，My约-1020 N·m（具体数值取决于示例数据生成函数，此处仅示意量级）。