import numpy as np
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')
# ============================================================================
# 1. 正确的翼型生成函数
# ============================================================================
def generate_naca_airfoil(m=0.02, p=0.4, t=0.12, n_points=200):
"""生成NACA 4位数翼型坐标"""
# 生成x坐标（从后缘到前缘再到后缘，逆时针）
beta = np.linspace(0, 2*np.pi, n_points)
x = 0.5 * (1 - np.cos(beta)) # 余弦分布，前缘和后缘密集
# 中弧线
yc = np.zeros_like(x)
dyc_dx = np.zeros_like(x)
# 前段 (0 <= x <= p)
idx1 = x <= p
yc[idx1] = (m / p**2) * (2 * p * x[idx1] - x[idx1]**2)
dyc_dx[idx1] = (2 * m / p**2) * (p - x[idx1])
# 后段 (p <= x <= 1)
idx2 = x > p
yc[idx2] = (m / (1 - p)**2) * ((1 - 2 * p) + 2 * p * x[idx2] - x[idx2]**2)
dyc_dx[idx2] = (2 * m / (1 - p)**2) * (p - x[idx2])
# 厚度分布
yt = (t / 0.2) * (0.2969 * np.sqrt(x) - 0.1260 * x -
0.3516 * x**2 + 0.2843 * x**3 - 0.1015 * x**4)
# 角度
theta = np.arctan(dyc_dx)
# 上表面和下表面
xu = x - yt * np.sin(theta)
yu = yc + yt * np.cos(theta)
xl = x + yt * np.sin(theta)
yl = yc - yt * np.cos(theta)
# 合并坐标（逆时针：从后缘开始，下表面→前缘→上表面→后缘）
# 但为了简单，我们只取一半点，然后对称
x_upper = xu[:n_points//2]
y_upper = yu[:n_points//2]
x_lower = xl[n_points//2:][::-1] # 反转以使从后缘开始
y_lower = yl[n_points//2:][::-1]
X = np.concatenate([x_lower, x_upper])
Y = np.concatenate([y_lower, y_upper])
# 确保闭合
X[0] = X[-1] = 1.0
Y[0] = Y[-1] = 0.0
return X, Y
# ============================================================================
# 2. 正确的涡面元法实现
# ============================================================================
def correct_vortex_panel_method(X, Y, alpha_deg, V_inf=1.0):
"""
正确的涡面元法实现
参考：Anderson "Fundamentals of Aerodynamics" 中的标准方法
"""
N = len(X) - 1 # 面板数量
# 面板几何参数
x = X.copy()
y = Y.copy()
# 面板中点（控制点）
xc = 0.5 * (x[:-1] + x[1:])
yc = 0.5 * (y[:-1] + y[1:])
# 面板长度和角度
S = np.sqrt((x[1:] - x[:-1])**2 + (y[1:] - y[:-1])**2)
theta = np.arctan2(y[1:] - y[:-1], x[1:] - x[:-1])
# 法向量（指向翼型外部）
nx = -np.sin(theta)
ny = np.cos(theta)
# 初始化影响系数矩阵
A = np.zeros((N, N))
# 计算影响系数
for i in range(N): # 控制点i
for j in range(N): # 面板j
if i == j:
# 自身影响：半平面贡献
A[i, j] = 0.5
else:
# 计算面板j对控制点i的影响
# 使用涡段的诱导速度公式
x1 = x[j]
y1 = y[j]
x2 = x[j+1]
y2 = y[j+1]
# 计算中间变量
dx1 = xc[i] - x1
dy1 = yc[i] - y1
dx2 = xc[i] - x2
dy2 = yc[i] - y2
r1_sq = dx1**2 + dy1**2
r2_sq = dx2**2 + dy2**2
# 避免除零
if r1_sq < 1e-12 or r2_sq < 1e-12:
A[i, j] = 0.0
continue
r1 = np.sqrt(r1_sq)
r2 = np.sqrt(r2_sq)
# 计算角度
theta1 = np.arctan2(dy1, dx1)
theta2 = np.arctan2(dy2, dx2)
dtheta = theta2 - theta1
# 调整角度范围
if dtheta > np.pi:
dtheta -= 2*np.pi
elif dtheta < -np.pi:
dtheta += 2*np.pi
# 计算对数项
log_term = np.log(r2/r1)
# 计算诱导速度在法向的分量
# 公式来自标准涡面元法教材
u_ind = (dtheta/(2*np.pi)) * np.cos(theta[i]) - (log_term/(2*np.pi)) * np.sin(theta[i])
v_ind = (dtheta/(2*np.pi)) * np.sin(theta[i]) + (log_term/(2*np.pi)) * np.cos(theta[i])
# 投影到法向
A[i, j] = u_ind * nx[i] + v_ind * ny[i]
# 右端向量：来流在法向的分量
alpha_rad = np.deg2rad(alpha_deg)
b = -V_inf * (np.cos(alpha_rad) * nx + np.sin(alpha_rad) * ny)
# 添加Kutta条件：γ_1 + γ_N = 0
A_aug = np.vstack([A, np.ones(N)]) # 添加一行全1
b_aug = np.append(b, 0.0) # 添加Kutta条件
# 使用最小二乘法求解（更稳定）
gamma, residuals, rank, s = np.linalg.lstsq(A_aug, b_aug, rcond=None)
# 计算表面速度
V_tau = V_inf * (np.cos(alpha_rad) * np.cos(theta) + np.sin(alpha_rad) * np.sin(theta))
# 涡诱导的切向速度
for i in range(N):
for j in range(N):
if i != j:
# 重新计算切向影响
x1 = x[j]
y1 = y[j]
x2 = x[j+1]
y2 = y[j+1]
dx1 = xc[i] - x1
dy1 = yc[i] - y1
dx2 = xc[i] - x2
dy2 = yc[i] - y2
r1 = np.sqrt(dx1**2 + dy1**2)
r2 = np.sqrt(dx2**2 + dy2**2)
theta1 = np.arctan2(dy1, dx1)
theta2 = np.arctan2(dy2, dx2)
dtheta = theta2 - theta1
if dtheta > np.pi:
dtheta -= 2*np.pi
elif dtheta < -np.pi:
dtheta += 2*np.pi
log_term = np.log(r2/r1)
u_ind = (dtheta/(2*np.pi)) * np.cos(theta[i]) - (log_term/(2*np.pi)) * np.sin(theta[i])
v_ind = (dtheta/(2*np.pi)) * np.sin(theta[i]) + (log_term/(2*np.pi)) * np.cos(theta[i])
V_tau[i] += gamma[j] * (u_ind * np.cos(theta[i]) + v_ind * np.sin(theta[i]))
# 后缘速度修正（强制Kutta条件）
V_te = 0.5 * (V_tau[0] + V_tau[-1])
V_tau[0] = V_tau[-1] = V_te
return gamma, V_tau, S, xc, yc, nx, ny, theta
# ============================================================================
# 3. 正确的气动系数计算
# ============================================================================
def correct_aerodynamic_coefficients(gamma, V_tau, S, nx, ny, xc, yc, alpha_deg, chord=1.0, V_inf=1.0):
"""正确的气动系数计算"""
# 1. 总环量（注意：gamma是单位长度的涡强）
Gamma = np.sum(gamma * S)
# 2. Kutta-Joukowski升力系数
C_L_KJ = 2 * Gamma / (V_inf * chord)
# 3. 压力系数
Cp = 1.0 - (V_tau / V_inf)**2
# 4. 压力积分法升力系数
alpha_rad = np.deg2rad(alpha_deg)
C_L_pressure = 0.0
C_m_LE = 0.0
for i in range(len(S)):
# 升力方向分量
lift_dir = -np.sin(alpha_rad) * nx[i] + np.cos(alpha_rad) * ny[i]
# 压力在升力方向的分量
C_L_pressure += -Cp[i] * lift_dir * S[i]
# 力矩（对前缘(0,0)）
dFx = -Cp[i] * nx[i] * S[i]
dFy = -Cp[i] * ny[i] * S[i]
C_m_LE += (xc[i] - 0.25) * dFy - yc[i] * dFx # 对1/4弦长点
# 无量纲化
C_L_pressure /= chord
C_m_LE /= (chord**2)
return C_L_KJ, C_L_pressure, C_m_LE, Cp, Gamma
# ============================================================================
# 4. 薄翼理论（用于对比）
# ============================================================================
def thin_airfoil_theory(alpha_deg, m=0.02, p=0.4):
"""薄翼理论计算"""
alpha_rad = np.deg2rad(alpha_deg)
# 零升攻角
alpha_L0 = -2 * m * (1 - p) / np.pi
alpha_L0_rad = alpha_L0
# 升力系数
C_L = 2 * np.pi * (alpha_rad - alpha_L0_rad)
# 零升力矩系数
C_m0 = -m * (1 - 2*p) / 4
# 对前缘力矩系数
C_m = C_m0 - C_L/4
return C_L, C_m
# ============================================================================
# 5. 任务5：单个攻角分析
# ============================================================================
def task5_single_analysis(alpha_deg=5.0):
"""任务5：单个攻角分析"""
print(f"\n{'='*60}")
print(f"任务5：单个攻角分析 (α={alpha_deg}°)")
print(f"{'='*60}")
# 生成翼型
X, Y = generate_naca_airfoil(m=0.02, p=0.4, t=0.12, n_points=150)
chord = 1.0
# 运行涡面元法
gamma, V_tau, S, xc, yc, nx, ny, theta = correct_vortex_panel_method(X, Y, alpha_deg)
# 计算气动系数
C_L_KJ, C_L_pressure, C_m, Cp, Gamma = correct_aerodynamic_coefficients(
gamma, V_tau, S, nx, ny, xc, yc, alpha_deg, chord)
# 薄翼理论
C_L_thin, C_m_thin = thin_airfoil_theory(alpha_deg)
print(f"计算结果:")
print(f" 总环量 Γ: {Gamma:.6f}")
print(f" Kutta-Joukowski法: C_L = {C_L_KJ:.6f}")
print(f" 压力积分法: C_L = {C_L_pressure:.6f}")
print(f" 薄翼理论: C_L = {C_L_thin:.6f}")
print(f" 力矩系数(1/4弦长): C_m = {C_m:.6f}")
print(f" 薄翼理论力矩: C_m = {C_m_thin:.6f}")
# XFOIL参考值
xfoil_data = {
0.0: 0.230,
2.0: 0.420,
5.0: 0.704,
8.0: 0.985,
10.0: 1.170
}
if alpha_deg in xfoil_data:
C_L_xfoil = xfoil_data[alpha_deg]
print(f"\nXFOIL参考值: C_L = {C_L_xfoil:.4f}")
print(f"误差分析:")
print(f" K-J法误差: {abs(C_L_KJ - C_L_xfoil)/C_L_xfoil*100:.2f}%")
print(f" 压力积分法误差: {abs(C_L_pressure - C_L_xfoil)/C_L_xfoil*100:.2f}%")
# 可视化
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
# 1. 翼型形状
axes[0, 0].plot(X, Y, 'k-', linewidth=2)
axes[0, 0].fill_between(X, Y, alpha=0.2, color='gray')
axes[0, 0].set_xlabel('x/c')
axes[0, 0].set_ylabel('y/c')
axes[0, 0].set_title('NACA 22012 翼型')
axes[0, 0].axis('equal')
axes[0, 0].grid(True, alpha=0.3)
# 2. 压力系数分布
sorted_idx = np.argsort(xc)
axes[0, 1].plot(xc[sorted_idx], -Cp[sorted_idx], 'b-', linewidth=2)
axes[0, 1].set_xlabel('x/c')
axes[0, 1].set_ylabel('-Cp')
axes[0, 1].set_title(f'压力系数分布 (α={alpha_deg}°)')
axes[0, 1].grid(True, alpha=0.3)
# 3. 升力系数对比
methods = ['K-J法', '压力积分', '薄翼理论']
values = [C_L_KJ, C_L_pressure, C_L_thin]
colors = ['blue', 'green', 'red']
bars = axes[1, 0].bar(methods, values, color=colors, alpha=0.7)
axes[1, 0].set_ylabel('升力系数 C_L')
axes[1, 0].set_title('升力系数对比')
axes[1, 0].grid(True, alpha=0.3, axis='y')
for bar, value in zip(bars, values):
height = bar.get_height()
axes[1, 0].text(bar.get_x() + bar.get_width()/2, height + 0.01,
f'{value:.4f}', ha='center', va='bottom')
# 4. 环量分布
axes[1, 1].plot(xc, gamma, 'r-', linewidth=2)
axes[1, 1].set_xlabel('x/c')
axes[1, 1].set_ylabel('涡强 γ')
axes[1, 1].set_title('涡强分布')
axes[1, 1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig(f'task5_correct_alpha_{alpha_deg}.png', dpi=300, bbox_inches='tight')
plt.show()
return C_L_KJ, C_L_pressure, C_L_thin, C_m, C_m_thin
# ============================================================================
# 6. 任务6：攻角扫描
# ============================================================================
def task6_alpha_sweep_analysis():
"""任务6：攻角扫描分析"""
print(f"\n{'='*60}")
print("任务6：攻角扫描分析")
print(f"{'='*60}")
# 生成翼型
X, Y = generate_naca_airfoil(m=0.02, p=0.4, t=0.12, n_points=100)
chord = 1.0
# 攻角范围
alphas = np.arange(-5, 11, 1.0)
n_alphas = len(alphas)
# 存储结果
C_L_KJ_list = np.zeros(n_alphas)
C_L_pressure_list = np.zeros(n_alphas)
C_m_list = np.zeros(n_alphas)
C_L_thin_list = np.zeros(n_alphas)
C_m_thin_list = np.zeros(n_alphas)
print(f"\n计算进度:")
print(f"{'攻角':^8} {'C_L(KJ)':^12} {'C_L(薄翼)':^12}")
print("-" * 40)
for i, alpha in enumerate(alphas):
try:
# 涡面元法
gamma, V_tau, S, xc, yc, nx, ny, theta = correct_vortex_panel_method(X, Y, alpha)
C_L_KJ, C_L_pressure, C_m, Cp, Gamma = correct_aerodynamic_coefficients(
gamma, V_tau, S, nx, ny, xc, yc, alpha, chord)
# 薄翼理论
C_L_thin, C_m_thin = thin_airfoil_theory(alpha)
C_L_KJ_list[i] = C_L_KJ
C_L_pressure_list[i] = C_L_pressure
C_m_list[i] = C_m
C_L_thin_list[i] = C_L_thin
C_m_thin_list[i] = C_m_thin
print(f"{alpha:^8.1f} {C_L_KJ:^12.4f} {C_L_thin:^12.4f}")
except Exception as e:
print(f"α={alpha}° 计算失败: {e}")
C_L_thin, C_m_thin = thin_airfoil_theory(alpha)
C_L_thin_list[i] = C_L_thin
C_m_thin_list[i] = C_m_thin
C_L_KJ_list[i] = np.nan
C_L_pressure_list[i] = np.nan
C_m_list[i] = np.nan
print("=" * 40)
# 计算升力线斜率
valid_mask = ~np.isnan(C_L_KJ_list)
if np.sum(valid_mask) > 2:
valid_alphas = alphas[valid_mask]
valid_C_L = C_L_KJ_list[valid_mask]
# 线性回归
coeff = np.polyfit(valid_alphas, valid_C_L, 1)
slope_deg = coeff[0] # 每度的斜率
slope_rad = slope_deg * 180/np.pi # 转换为每弧度
print(f"\n升力线斜率分析:")
print(f" 涡面元法: {slope_rad:.4f} /弧度")
print(f" 薄翼理论: {2*np.pi:.4f} /弧度")
print(f" 相对误差: {abs(slope_rad - 2*np.pi)/(2*np.pi)*100:.2f}%")
# 可视化
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# 1. 升力系数曲线
axes[0, 0].plot(alphas, C_L_KJ_list, 'bo-', linewidth=2, markersize=6, label='涡面元法')
axes[0, 0].plot(alphas, C_L_thin_list, 'r-', linewidth=2, alpha=0.7, label='薄翼理论')
# XFOIL数据点
xfoil_alphas = [0.0, 2.0, 5.0, 8.0, 10.0]
xfoil_cl = [0.230, 0.420, 0.704, 0.985, 1.170]
axes[0, 0].plot(xfoil_alphas, xfoil_cl, 'g^', markersize=8, label='XFOIL')
axes[0, 0].set_xlabel('攻角 α [度]')
axes[0, 0].set_ylabel('升力系数 C_L')
axes[0, 0].set_title('升力系数随攻角变化')
axes[0, 0].grid(True, alpha=0.3)
axes[0, 0].legend()
# 2. 力矩系数曲线
axes[0, 1].plot(alphas, C_m_list, 'bo-', linewidth=2, markersize=6, label='涡面元法')
axes[0, 1].plot(alphas, C_m_thin_list, 'r-', linewidth=2, alpha=0.7, label='薄翼理论')
axes[0, 1].set_xlabel('攻角 α [度]')
axes[0, 1].set_ylabel('力矩系数 C_m (1/4弦长)')
axes[0, 1].set_title('力矩系数随攻角变化')
axes[0, 1].grid(True, alpha=0.3)
axes[0, 1].legend()
# 3. 两种方法对比
axes[1, 0].plot(alphas, C_L_KJ_list, 'b-', linewidth=2, label='K-J法')
axes[1, 0].plot(alphas, C_L_pressure_list, 'g--', linewidth=2, label='压力积分法')
axes[1, 0].set_xlabel('攻角 α [度]')
axes[1, 0].set_ylabel('升力系数 C_L')
axes[1, 0].set_title('两种计算方法对比')
axes[1, 0].grid(True, alpha=0.3)
axes[1, 0].legend()
# 4. 升力系数差异
valid_idx = ~np.isnan(C_L_KJ_list)
if np.any(valid_idx):
diff = np.abs(C_L_KJ_list[valid_idx] - C_L_thin_list[valid_idx])
axes[1, 1].plot(alphas[valid_idx], diff, 'r-o', linewidth=2, markersize=6)
axes[1, 1].set_xlabel('攻角 α [度]')
axes[1, 1].set_ylabel('|C_L(涡面元) - C_L(薄翼)|')
axes[1, 1].set_title('与薄翼理论的差异')
axes[1, 1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('task6_correct_alpha_sweep.png', dpi=300, bbox_inches='tight')
plt.show()
# 保存数据
results = np.column_stack([alphas, C_L_KJ_list, C_L_pressure_list, C_m_list,
C_L_thin_list, C_m_thin_list])
np.savetxt('task6_correct_results.csv', results,
delimiter=',',
header='alpha_deg,C_L_KJ,C_L_pressure,C_m,C_L_thin,C_m_thin',
fmt='%.6f',
comments='')
return alphas, C_L_KJ_list, C_L_pressure_list, C_m_list, C_L_thin_list, C_m_thin_list
# ============================================================================
# 7. 基于解析解的分析（备选方案）
# ============================================================================
def analytical_analysis():
"""基于解析解的分析（当数值方法失败时使用）"""
print(f"\n{'='*60}")
print("基于解析解的分析（备选方案）")
print(f"{'='*60}")
# 攻角范围
alphas = np.arange(-5, 11, 1.0)
# 计算薄翼理论结果
C_L_thin = np.zeros_like(alphas)
C_m_thin = np.zeros_like(alphas)
for i, alpha in enumerate(alphas):
C_L_thin[i], C_m_thin[i] = thin_airfoil_theory(alpha)
# XFOIL参考数据
xfoil_alphas = np.array([0.0, 2.0, 5.0, 8.0, 10.0])
xfoil_cl = np.array([0.230, 0.420, 0.704, 0.985, 1.170])
xfoil_cm = np.array([-0.030, -0.130, -0.210, -0.290, -0.350])
print(f"\n薄翼理论计算结果:")
print(f"{'攻角':^8} {'C_L(薄翼)':^12} {'C_m(薄翼)':^12}")
print("-" * 40)
for i, alpha in enumerate(alphas):
print(f"{alpha:^8.1f} {C_L_thin[i]:^12.4f} {C_m_thin[i]:^12.4f}")
print("=" * 40)
# 计算升力线斜率
coeff = np.polyfit(alphas, C_L_thin, 1)
slope_deg = coeff[0]
slope_rad = slope_deg * 180/np.pi
print(f"\n升力线斜率分析:")
print(f" 薄翼理论斜率: {slope_rad:.6f} /弧度")
print(f" 理论值(2π): {2*np.pi:.6f} /弧度")
print(f" 相对误差: {abs(slope_rad - 2*np.pi)/(2*np.pi)*100:.2f}%")
# 可视化
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
# 升力系数曲线
axes[0].plot(alphas, C_L_thin, 'b-', linewidth=2, label='薄翼理论')
axes[0].plot(xfoil_alphas, xfoil_cl, 'ro', markersize=8, label='XFOIL数据')
axes[0].set_xlabel('攻角 α [度]')
axes[0].set_ylabel('升力系数 C_L')
axes[0].set_title('NACA 22012 升力系数')
axes[0].grid(True, alpha=0.3)
axes[0].legend()
# 力矩系数曲线
axes[1].plot(alphas, C_m_thin, 'b-', linewidth=2, label='薄翼理论')
axes[1].plot(xfoil_alphas, xfoil_cm, 'ro', markersize=8, label='XFOIL数据')
axes[1].set_xlabel('攻角 α [度]')
axes[1].set_ylabel('力矩系数 C_m')
axes[1].set_title('NACA 22012 力矩系数')
axes[1].grid(True, alpha=0.3)
axes[1].legend()
plt.tight_layout()
plt.savefig('analytical_results.png', dpi=300, bbox_inches='tight')
plt.show()
return alphas, C_L_thin, C_m_thin, xfoil_alphas, xfoil_cl, xfoil_cm
# ============================================================================
# 8. 主函数
# ============================================================================
def main():
"""主函数"""
print("="*60)
print("NACA 22012 翼型空气动力学分析")
print("="*60)
try:
# 任务5：单个攻角分析
print("\n>>> 执行任务5：单个攻角分析")
print("(使用5°攻角作为示例)")
task5_results = task5_single_analysis(alpha_deg=5.0)
if task5_results[0] is not None:
C_L_KJ, C_L_pressure, C_L_thin, C_m, C_m_thin = task5_results
print(f"\n任务5结果总结:")
print(f" Kutta-Joukowski法: C_L = {C_L_KJ:.4f}")
print(f" 压力积分法: C_L = {C_L_pressure:.4f}")
print(f" 薄翼理论: C_L = {C_L_thin:.4f}")
print(f" 力矩系数: C_m = {C_m:.4f}")
# 任务6：攻角扫描
print("\n>>> 执行任务6：攻角扫描分析")
print("攻角范围: -5° 到 10°, 间隔 1°")
task6_results = task6_alpha_sweep_analysis()
print("\n" + "="*60)
print("涡面元法分析完成!")
print("="*60)
except Exception as e:
print(f"\n涡面元法计算失败: {e}")
print("使用解析解作为备选方案...")
# 使用解析解
analytical_results = analytical_analysis()
print("\n" + "="*60)
print("基于解析解的分析完成!")
print("="*60)
# ============================================================================
# 9. 运行程序
# ============================================================================
if __name__ == "__main__":
main()
扫一扫
1201
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
