3、翼型几何建模：NACA 4位数字翼型的数学表达式、坐标生成方法、使用Python生成翼型坐标点

各位同学，咱们今天聊聊翼型建模。说实话，这是整个课程里最“硬核”的一章，但也是最有意思的。你想想看，一个简单的数字组合，比如“2412”，就能定义出一个完整的翼型形状——这本身就是空气动力学的魅力所在。

我个人习惯把NACA 4位数字翼型比作“空气动力学的乐高积木”。你掌握了它的生成逻辑，后面那些更复杂的5位、6位数字翼型，甚至超临界翼型，理解起来都会轻松很多。

3.1 NACA 4位数字翼型的数学表达式

NACA 4位数字翼型，说白了就是用四个数字来编码翼型的几何特征。比如NACA 2412：

• 第一位数字（2）：最大弯度占弦长的百分比。这里就是2%。
• 第二位数字（4）：最大弯度位置距前缘的距离，以弦长的十分之一表示。这里是0.4倍弦长位置。
• 最后两位数字（12）：最大厚度占弦长的百分比。这里是12%。

嗯，这里要注意：弯度和厚度是两个完全不同的概念。弯度决定了翼型的不对称程度，说白了就是产生升力的“底子”；厚度则影响阻力特性和结构强度。我在项目中遇到过有人把这两个搞混，结果设计出来的翼型飞起来像块砖头……

数学上，NACA 4位数字翼型由三部分组成：

1. 中弧线（Camber Line）：翼型上下表面的“中间线”。
2. 厚度分布（Thickness Distribution）：沿弦长方向对称分布的厚度函数。
3. 上下表面坐标：由中弧线加上或减去厚度得到。

中弧线的数学表达式是这样的：

对于前段（0 ≤ x ≤ p）：
  y_c = (m / p²) * (2p*x - x²)

对于后段（p ≤ x ≤ 1）：
  y_c = (m / (1-p)²) * (1 - 2p + 2p*x - x²)

其中：
  m = 最大弯度（第一位数字/100）
  p = 最大弯度位置（第二位数字/10）

厚度分布则统一使用这个公式：

y_t = (t / 0.2) * (0.2969*sqrt(x) - 0.1260*x - 0.3516*x² + 0.2843*x³ - 0.1015*x⁴)

其中：
  t = 最大厚度（最后两位数字/100）

这里有个坑，我必须要提醒你：厚度分布公式最后一项的系数是-0.1015，不是-0.1036。很多老教材上写的是-0.1036，但那是早期的手算近似值。用-0.1015才能保证后缘闭合。我曾经因为这个系数搞错过，画出来的翼型后缘像鱼尾巴一样张开……

3.2 坐标生成方法

有了数学表达式，生成坐标点就很简单了。基本流程如下：

1. 在弦长方向（x轴）上取一系列点，通常从0到1。
2. 计算每个x位置的中弧线高度y_c。
3. 计算每个x位置的厚度y_t。
4. 上下表面坐标分别为：
   • 上表面：x_u = x - y_t * sin(θ)，y_u = y_c + y_t * cos(θ)
   • 下表面：x_l = x + y_t * sin(θ)，y_l = y_c - y_t * cos(θ)
   其中θ是中弧线的切线角度，θ = arctan(dy_c/dx)。

你可能会问：“为什么要用sin和cos修正？直接加减不行吗？”

嗯，直接加减确实简单，但那样得到的翼型在中弧线曲率大的地方会失真。用垂直修正的方式，说白了就是沿着中弧线的法线方向加厚度，这样生成的翼型才是“干净”的。我在做高升力翼型设计时，就吃过这个亏——直接加减导致前缘形状畸变，CFD计算出来的升力系数差了10%。

3.3 使用Python生成翼型坐标点

好了，理论说完了，咱们直接上代码。我个人习惯用NumPy来处理数组运算，这样代码简洁，效率也高。

import numpy as np

def naca_4_digit(designation, n_points=200):
    """
    生成NACA 4位数字翼型坐标
    
    参数：
        designation: 翼型编号，如'2412'
        n_points: 生成的点数（默认200）
    
    返回：
        x_u, y_u: 上表面坐标
        x_l, y_l: 下表面坐标
    """
    # 解析翼型编号
    m = int(designation[0]) / 100.0      # 最大弯度
    p = int(designation[1]) / 10.0       # 最大弯度位置
    t = int(designation[2:]) / 100.0     # 最大厚度
    
    # 生成x坐标（前缘加密）
    beta = np.linspace(0, np.pi, n_points)
    x = (1 - np.cos(beta)) / 2.0
    
    # 厚度分布
    y_t = (t / 0.2) * (0.2969*np.sqrt(x) - 0.1260*x - 0.3516*x**2 + 
                       0.2843*x**3 - 0.1015*x**4)
    
    # 中弧线
    y_c = np.zeros_like(x)
    dy_c_dx = np.zeros_like(x)
    
    # 前段
    front = x <= p
    y_c[front] = (m / p**2) * (2*p*x[front] - x[front]**2)
    dy_c_dx[front] = (2*m / p**2) * (p - x[front])
    
    # 后段
    back = x > p
    y_c[back] = (m / (1-p)**2) * (1 - 2*p + 2*p*x[back] - x[back]**2)
    dy_c_dx[back] = (2*m / (1-p)**2) * (p - x[back])
    
    # 中弧线法线角度
    theta = np.arctan(dy_c_dx)
    
    # 上下表面坐标
    x_u = x - y_t * np.sin(theta)
    y_u = y_c + y_t * np.cos(theta)
    x_l = x + y_t * np.sin(theta)
    y_l = y_c - y_t * np.cos(theta)
    
    return x_u, y_u, x_l, y_l

# 生成NACA 2412翼型
x_u, y_u, x_l, y_l = naca_4_digit('2412', n_points=200)

# 保存到文件
with open('naca2412.dat', 'w') as f:
    f.write('NACA 2412 Airfoil\n')
    # 从上表面后缘开始，逆时针绕一圈
    for i in range(len(x_u)-1, -1, -1):
        f.write(f'{x_u[i]:.6f}  {y_u[i]:.6f}\n')
    for i in range(len(x_l)):
        f.write(f'{x_l[i]:.6f}  {y_l[i]:.6f}\n')

print('翼型坐标已生成，保存到 naca2412.dat')

3.4 知识体系总览

下面这张图，是我自己总结的NACA 4位数字翼型建模的核心逻辑。你把它理解了，整个流程就刻在脑子里了。

3.5 常见问题与避坑指南

好了，这一章的内容就到这里。你把这个NACA 4位数字翼型的生成逻辑吃透了，后面那些复杂的翼型建模，说白了就是在这个基础上加一些“花活”。代码我已经给了，建议你动手跑一遍，看看不同编号的翼型长什么样。嗯，动手才是最好的学习方式。