02年世界杯韩国黑哨_曲棍球世界杯 - guanchafang.com

python——薄翼型理论及求解二维机翼压强系数

薄翼型理论及求解二维翼型压强系数

python——薄翼型理论及求解二维机翼压强系数薄翼型理论——翼型的攻角、弯度效应NACA四位数翼型翼型的攻角、弯度效应

翼型的厚度效应求解二维机翼压强系数二维机翼压强系数绘图部分

薄翼型理论——翼型的攻角、弯度效应

在无粘流中，对于相对厚度比较小的机翼和翼型，假设其厚度为零，而以其中弧线来代替，其升力效应则用分布涡来体现，这种理论称为薄翼理论。

NACA四位数翼型

NACA翼型是美国国家航空咨询委员会（NACA，NASA前身）开发的一系列翼型。 每个翼型的代号由“NACA”这四个字母与一串数字组成，将这串数字所描述的几何参数代入特定方程中即可得到翼型的精确形状。 NACA四位数翼型的表达方式为： 中弧线方程为： 其中f为中弧线最高点的纵坐标（单位为0.01），p为中弧线最高点的弦向位置（单位为0.1）。

翼型的厚度分布为： 其中t为厚度（单位为0.01）。

上下翼面的坐标为： 其中的

θ

\theta

θ 为中弧线在弦向位置x处切线的倾斜角。

代码如下：

import numpy as np

import sympy as sp

import matplotlib.pyplot as plt

import scipy

class ThinAirfoil:

def __init__(self, f, p, t, b, alpha, V_oo, acc): # f中弧线最高点纵坐标 p最高点的弦向位置 t厚度

self.f = f

self.p = p

self.t = t

self.b = b

self.alpha = alpha * 3.1415296 / 180

self.V_oo = V_oo

self.acc = acc

def naca_four(self): # f中弧线最高点纵坐标 p最高点的弦向位置 t厚度

x = sp.symbols('x')

self.f, self.p, self.t = 0.01 * float(self.f), 0.1 * float(self.p), 0.01 * float(self.t)

y_t = 5 * self.t * (

0.2969 * (x ** 0.5) - 0.126 * x - 0.3516 * x * x + 0.28433 * (x ** 3) - 0.10363 * (x ** 4)) # 厚度分布

if self.p != 0:

y_f = sp.Piecewise(((self.f / (self.p ** 2)) * (2 * self.p * x - x ** 2), x <= self.p),

((self.f / (1 - self.p) ** 2) * ((1 - 2 * self.p) + 2 * self.p * x - x ** 2), x > self.p)) # 中弧线

else:

y_f = -1e-16 * x * x

y_f_np = sp.lambdify(x, y_f, 'numpy')

theta = sp.atan(sp.diff(y_f, x))

x_upper = x - y_t * sp.sin(theta)

y_upper = y_f + y_t * sp.cos(theta)

x_lower = x + y_t * sp.sin(theta)

y_lower = y_f - y_t * sp.cos(theta)

x_upper_np = sp.lambdify(x, x_upper, 'numpy')

y_upper_np = sp.lambdify(x, y_upper, 'numpy')

x_lower_np = sp.lambdify(x, x_lower, 'numpy')

y_lower_np = sp.lambdify(x, y_lower, 'numpy')

x_vals = np.arange(0, 1, self.acc)

xu_val = x_upper_np(x_vals)

yu_val = y_upper_np(x_vals)

xl_val = x_lower_np(x_vals)

yl_val = y_lower_np(x_vals)

return (xu_val, yu_val), (xl_val, yl_val), y_f, y_f_np, y_t

绘制NACA4412翼型如下： 其他翼型在此不再赘述。

翼型的攻角、弯度效应

薄翼型理论的关键在于求解涡密度分布。通过三角变化得到的涡密度分布为： 变量置换的几何关系为： 代码如下：

def thinAirfoil(self, y_f): # alpha攻角 b弦长 # 攻角 弯度

x, theta = sp.symbols('x theta')

dy_dx = sp.diff(y_f, x)

dy_dx_theta = dy_dx.subs(x, 0.5 * self.b * (1 - sp.cos(theta)))

alpha_0 = sp.integrate((1 - sp.cos(theta)) * dy_dx_theta, (theta, 0, 3.1415926)) / 3.1415926 # 零升攻角

Cy = 2 * sp.pi * (self.alpha - alpha_0) # 升力系数

miu_0 = - sp.integrate((1 - sp.cos(2 * theta)) * dy_dx_theta, (theta, 0, sp.pi)) / 4

mz_0 = -2 * (miu_0 + sp.pi * alpha_0 / 4) # 零升俯仰力矩系数

mz = Cy / 4 + mz_0 # 俯仰力矩系数

xp = 0.25 + mz_0 / Cy # 压心相对位置

n = sp.symbols('n')

A_0 = self.alpha - (sp.integrate(dy_dx_theta, (theta, 0, sp.pi))) / sp.pi

A_n = 2 * sp.integrate(dy_dx_theta * sp.cos(n * theta), (theta, 0, sp.pi)) / sp.pi

v_y = self.V_oo * (dy_dx - self.alpha) # 下洗速度

gamma = 0.5 * self.b * self.V_oo * Cy # 总涡强

gamma_theta = 2 * self.V_oo * (A_0 * sp.cot(theta / 2) + sp.Sum(A_n * sp.sin(n * theta), (n, 1, 100))) # 涡密度

gamma_x = gamma_theta.subs(theta, sp.acos(1 - 2 * x / self.b))

x_vals = np.arange(0, self.b, self.acc * self.b)

gamma_func = sp.lambdify(x, gamma_x, 'numpy')

v_y_func = sp.lambdify(x, v_y, 'numpy')

gamma_vals = gamma_func(x_vals) # gamma/V_oo为压力系数

v_x_vals = gamma_vals / 2

C_p = gamma_vals / self.V_oo

v_y_vals = v_y_func(x_vals)

return gamma_vals, v_x_vals, v_y_vals, alpha_0

此代码不仅求解了薄翼型理论的涡密度分布，而且求解了零升攻角、零升俯仰力矩系数、压心相对位置等重要气动参数，具体公式不再赘述，可参考相关空气动力学书籍。

翼型的厚度效应

由于要求解翼型的压强系数，因此除了由攻角、弯度效应引起的压力分布，还要考虑到厚度效应引起的压力分布。

def thickPr(self, y_t): # 翼型的厚度问题

x, xi = sp.symbols('x xi')

dy_dx = sp.diff(y_t, x)

dy_dx = dy_dx.subs(x, xi)

q_x = 2 * self.V_oo * dy_dx # 密度分布

v_y = self.V_oo * dy_dx

v_y_func = sp.lambdify(xi, v_y, 'numpy')

x_vals_s = np.arange(0, self.b, self.acc * self.b)

v_y_vals = v_y_func(x_vals_s)

v_x_list = []

x_vals = [i for i in np.arange(self.acc, self.b, self.acc * self.b)]

for x_val in x_vals:

temp = q_x / ((x_val - xi) * 2 * sp.pi)

integrand = sp.lambdify(xi, temp, 'numpy')

eps = self.acc * 0.05

interval_1 = [0, x_val - eps]

interval_2 = [x_val + eps, self.b]

result_1 = scipy.integrate.quad(integrand, interval_1[0], interval_1[1])

result_2 = scipy.integrate.quad(integrand, interval_2[0], interval_2[1])

v_x = result_2[0] + result_1[0]

err = result_2[1] + result_1[1] # 积分误差

v_x_list.append(v_x)

return v_x_list, v_y_vals

同样，具体公式不再赘述，可参考相关空气动力学书籍。

求解二维机翼压强系数

二维机翼压强系数

由压强系数的公式可知，压强系数不满足叠加原理；但对于二维机翼采用小扰动线化理论，对压强系数只取一阶近似，那么压强系数Cp与旋涡产生的诱导速度vx成正比，因此可以直接叠加攻角、弯度产生的Cp与厚度产生的Cp。

def solveCp(self, v_x, v_x_list, v_y, v_y_list): # 求压强系数

v_y_list = list(v_y_list)

v_x = list(v_x)

v_y = list(v_y)

v_x.pop(0) # vx弯度&攻角 vx_list厚度

v_y.pop(0)

Cp_upper, Cp_lower = [], []

for i in range(0, len(v_x)):

# Cp_upper.append(-2*(v_x[i]+v_x_list[i])/V_oo-((v_y[i]+v_y_list[i])**2)/(V_oo**2)) # 二阶近似

# Cp_lower.append(-2*(-v_x[i]+v_x_list[i])/V_oo-((v_y[i]+v_y_list[i])**2)/(V_oo**2))

# Cp_upper.append(1-((v_x[i]+v_x_list[i])**2+(v_y[i]+v_y_list[i])**2)/(V_oo**2))

# Cp_lower.append(1-((v_x[i]-v_x_list[i])**2+(v_y[i]+v_y_list[i])**2)/(V_oo**2))

Cp_u = -2 * (v_x[i] + v_x_list[i]) / self.V_oo

if Cp_u <= 1:

Cp_upper.append(Cp_u) # 一次近似效果反而更好

else:

Cp_upper.append(1-((v_x[i]+v_x_list[i])**2+(v_y[i]+v_y_list[i])**2)/(self.V_oo**2))

Cp_l = -2 * (-v_x[i] + v_x_list[i]) / self.V_oo

if Cp_l <= 1:

Cp_lower.append(Cp_l) # 一次近似效果反而更好

else:

Cp_lower.append(1 - ((v_x[i] - v_x_list[i]) ** 2 + (v_y[i] + v_y_list[i]) ** 2) / (self.V_oo ** 2))

return Cp_upper, Cp_lower

与xflr5软件相比，不同攻角、来流速度及翼型情况的Cp分布如下：

一阶近似的效果在大部分情况下与xflr5软件的仿真匹配度较好。推测的原因在于线化理论对来流的压缩总体不足，但是舍弃的高次项又很好地弥补了这一点。

此外，本程序在前缘或后缘的压力分布估计误差较大，原因也许在于被积函数值趋近于无穷大，积分区间小，导致积分产生的误差较大，并会出现Cp>1的情况；而且在前缘处正压强梯度特别之大，边界层往往为紊流。因此，在机翼前缘部分采用了Cp的定义式进行计算。

绘图部分

def drawPlot(self):

a_vals = np.arange(0, 1, self.acc)

(xu_val, yu_val), (xl_val, yl_val), y_f, y_f_np, y_t = self.naca_four()

fig, ax = plt.subplots()

ax.plot(xu_val, yu_val)

ax.plot(xl_val, yl_val)

yf_val = y_f_np(a_vals)

ax.plot(a_vals, yf_val)

ax.set_aspect(1) # 绘制翼型

x_vals = np.arange(0, self.b, self.acc * self.b)

gamma_vals, v_x_vals, v_y_vals, _ = self.thinAirfoil(y_f)

fig, bx = plt.subplots()

bx.plot(x_vals, gamma_vals)

bx.plot(x_vals, v_y_vals) # 薄翼型的涡密度分布与下洗速度

v_x_list, v_y_list = self.thickPr(y_t)

Cp_upper, Cp_lower = self.solveCp(v_x_vals, v_x_list, v_y_vals, v_y_list)

x_vals = np.arange(self.acc, self.b, self.acc * self.b)

fig, cx = plt.subplots()

cx.plot(x_vals, Cp_upper)

cx.plot(x_vals, Cp_lower)

cx.invert_yaxis() # 绘制翼表面的压强系数

plt.show()

def main():

TA = ThinAirfoil(f=4, p=4, t=12, b=1, alpha=-4.3, V_oo=10, acc=0.005) # acc插值精度

TA.drawPlot()

if __name__ == '__main__':

main()

本程序十分粗糙，仅为本人为加深课程印象所写，仍有诸多不足，望能提出批评改进意见。

参考：《空气动力学》 钱翼稷编著