- resource1: xfoil接口程序,能够求解截面翼型不同攻角气动力,并使用Viterna method外推得到180°性能曲线。
- resource2: BEMT程序,用于计算螺旋桨推力效率等。
- resource3: 遗传算法优化程序,用于优化螺旋桨的几何参数。
- resource4: 结合MS-GAN的推力优化程序。
- resource5: 结合深度预测网络的多目标优化程序,使用代理模型预测螺旋桨力学、噪声数据并进行优化。
王晨宇
经典BEM理论的完全推导由Glauert在1935年完成,详见 Glauert, H., “Airplane propellers,” Aerodynamic theory, Springer, 1935, pp. 169–360
来流系数包括轴向和径向。几何参数描述了旋翼和流管几何结构。性能参数表示测量旋翼性能的指标。下标0是远场,上标是气流尾迹参数。
Flow variables | Geometric variables | Performance parameters |
---|---|---|
|
R 转子半径 | T, |
|
r 径向距离 | Q,$C_Q$ 扭矩,扭矩系数 |
|
|
P,$C_P$ 动力(power),动力系数 |
|
S 转子扫掠面积 | J 推进比 |
|
|
|
|
|
|
W 入流速度和矢量 |
|
|
|
|
|
|
|
FM 品质因子 |
|
|
|
|
- 分析通过转子平面气流动量的改变量。规定控制面“流管”
- 不考虑旋翼几何,只考虑动量通量
对于轴向动量,上游初始速度
考虑牛顿第二定律:
忽略体力,考虑质量通量在面上的积分
r代表径向距离,$p_0 (S_0 - S') $是作用流管两端压力,
但是不可以转化为微分: $dT = 2\pi \rho u_x (u'_x -V_x) r dr $ ,这是由于微分形式需要保证微元间无相互作用
化简未知量下游尾迹速度
$$\begin{equation} P = \frac{1}{2} \dot{m_{air} } (u'x -V^2_x) = \frac{1}{2} \dot{m{air}}(u'x+V_x) = \frac{1}{2}(u'_x + V_x)T \end{equation}$$
忽略了旋转诱导的滑流移动
单位时刻做的有用功即4功率
假设转子面轴向速度为远场上游与尾流速度平均
或者通常用的归一化诱导因子
由上述两式,有转化
$$\begin{equation} dT = 4 \pi \rho u_x (u_x - V_x) r dr\end{equation}$$ 或
但是公式(9)不可以在0表示来流速度,此时
定义效率 $$ \begin{equation} \eta = \frac{P_{out}}{P_{in}} \end{equation}$$
输出功率是远场流动速度乘推力(无旋转动力损失情况)。考虑轴向损失下,有
- 考虑螺旋桨滑流旋转损失(之前忽略)。由于桨叶与空气分子相互作用,增加了转向动量。
- 点质量假设,有角动量
$L = u_\theta m r^2$ 也可以写为$L = \omega_\theta m r^2$
由角动量守恒,导出切向诱导速度与最终尾流速度
而扭矩是由角动量单位时间增长导致的
定义切向诱导系数
在转子面前后用不可压缩无黏伯努利方程,通过转子平面的压强变化为
进一步在周向面元上由无黏不可压缩假定,有
对单位时间扭矩做功(输入),单位时间推力做工(输出)积分,有
代入公式(17)得到 $$\begin{equation}\eta_1 \eta_2 = \frac{V_x T^*}{\Omega Q} = \frac{V_x \pi \rho (V_\theta - u_\theta)u_\theta R^2}{\Omega \pi \rho u_x u_{\theta} R^3} = \frac{V_x (V_\theta - \frac{1}{2}u_\theta)}{u_x V_\theta} \end{equation}$$
对比效率公式(11)可以得到转动诱导效率 $$\begin{equation}\eta_2 = \frac{V_\theta - \frac{1}{2}u_\theta}{V_\theta} \end{equation}$$
- 对于低保真(low-fidelity)转子模型可用,考虑转子几何形状。每个叶片是个二维薄翼
- 表现在叶素翼型上的升力阻力,依赖来流速度W,雷诺数,和攻角
$\alpha $
每个叶素翼型上的推力扭力表示为
升阻力系数表示为
W 为 $ W = \sqrt{u^2_x + (V_\theta - \frac{1}{2} u_\theta)^2 }$
轴向$u_x 看作 V_x$周向诱导速度为0。
B为机翼数量,当推力(22)与扭力(23)相等时
积分有
由公式(30)(31)得到速度分量
将无量纲数
引入入流角$\theta$
公式(32)(33)代入(36),已知推进比J定义为 $ J = \frac{V_\infty}{n D_p} = \frac{\pi V_x}{\pi n R} = \pi \frac{V_x}{V_\theta} $
得到
最终化简(37)得到残差迭代式
在公式(38)中,各参数均为弦长
Ning, S. A., “A simple solution method for the blade element momentum equations with guaranteed convergence,” Wind Energy, Vol. 17, No. 9, 2014, pp. 1327–1345. Ning, A., “Using blade element momentum methods with gradient-based design optimization,” Structural and Multidisciplinary Optimization, 2021, pp. 1–24
两文中探索了收敛条件,用二分法找根时
if $ \phi = 0 , C_x = C_l , C_\theta = C_d $ and
if $ \phi = \frac{\pi}{2} , C_x = -C_d , C_\theta = C_l $ and $\alpha = \theta - \frac{\pi}{2} $ ,化简公式(38)为
如果两点异号,满足
公式(41)在对称和正弯度翼型成立,而对于公式(42),在极少数
附表:
叶素动量理论总结
转子表现系数总结
程序计算流程
注: 本节流程来自Marius L.Ruh and John T.Hwang,具体推导详询 [email protected]
Marius L.Ruh , John T.Hwang,"Robust modeling and optimal design of rotors using blade element momentum theory" AIAA,2021,pp.1-8.
进一步BEM修正推荐阅读报告
Winarto, Hadi. "BEMT algorithm for the prediction of the performance of arbitrary propellers." Melbourne: The Sir Lawrence Wackett Centre for Aerospace Design Technology, Royal Melbourne Institute of Technology,2004.
由于使用BEM需要360°全角外推 resource1
文件夹,操作语句默认在此文件夹路径下进行。
使用方法:
如果输入naca翼型
python extrapolate --foil 'naca4412' --Re 1e5 --base_angle '-10 10' --save_file 'example_extrapolate'
如果输入翼型文件
python extrapolate --foil 'example_extrapolate/ara_d_20.dat' --Re 1e5 --base_angle '-10 10' --save_file 'example_extrapolate'
翼型文件格式参照example_extrapolate/ara_d_20.dat
参数说明:
输入参数 | 说明 |
---|---|
foil | 输入为naca翼型或翼型坐标文件 |
Re | 输入计算的雷诺数 |
base_angle | 输入xfoil的基础计算范围,对应 [-base_angle[0],base_angle[1]],注意此值在合理范围内才会有结果 |
save_file | 保存的文件夹 |
再次提醒,base_angle非常重要,尤其是大雷诺数时,需要谨慎选择。
输出结果
输出结果可以在文件夹save_file
查看(上例为example_extrapolate中),输出文件有两个,文件名格式为 xfoil_Re...
和extrapolate_Re...
后缀名为雷诺数与翼型名称,靠下划线“_”分割。
xfoil_Re
第一行为 Re,alpha,cl,cd,cm
下方按表格形式排列,可以用deal_extrapolate.py
中的extract_matrix
读取
extrapolate_Re 前十四行为说明文档,说明了外推后的曲线属性,详情见每行的文件说明。 外推后范围为-180°-180°,每隔1°输出一次。 xfoil计算原理详情见第一章,下面介绍外推原理。
使用Viterna method外推,这种方法是外推翼型数据的最著名方法,在
Mahmuddin, Faisal, et al. "Airfoil lift and drag extrapolation with viterna and montgomerie methods." Energy Procedia 105 (2017): 811-816.
中有详细介绍。 首先使用实验或数值计算生成较小攻角的升阻力系数,使用下列公式生成失速角并外推数据
其中系数 $$\begin{equation} C_{D_{max}} \backsimeq 1.11 + 0.018 AR \end{equation}$$ $$\begin{equation} A_1 = \frac{C_{D_{max}}}{2} \end{equation}$$ $$\begin{equation} B_1 = C_{D_{max}} \end{equation}$$ $$\begin{equation} A_2 = (C_{L_{Install}}- C_{D_{max}} \sin \alpha_{stall} \cos \alpha_{stall} )\frac{\sin \alpha_{stall} }{\cos^2 \alpha_{stall}}\end{equation}$$ $$\begin{equation}B_2 = \frac{C_{D_{stall}} - C_{D_{max}} \sin^2 \alpha_{stall} }{\cos \alpha_{stall}} \end{equation}$$
公式(45)中的AR系数代表纵横比(aspect ratio),可从有限叶片长度影响平板假设的BEM应用中获得,一般情况下定位 0.75位置处的展弦比 R/c 。
在多数假设中AR的推荐值为10,使用不同的AR值对最终结果影响不大。在外推$\alpha > 90 \degree
使用airfoilprep.py
完成,最小Cd定为0.013800
本节以pyBEMT程序为基础,编写完成快速计算已知旋翼的螺旋桨性能程序。配置环境与所需程序见resource2
文件夹,操作语句默认在此文件夹路径下进行。
- 旋翼各截面几何与环境说明文件,参考
prop.ini
文件 - 各截面翼型的外推力学参数库(可以从2.2.1的resource1文件中得到),参考
BEMT_program/airfoil
文件夹中的dat文件
如果翼型在BEMT_program/airfoil
中,循环v_inf 变量,从0.1到20,得到20个值:
python BEMT --para_file 'example_BEMT/vali1/prop.ini' --V '20 0.1 20' --save_file 'example_BEMT/result' --draw_pic False
如果翼型在BEMT_program/airfoil
中,循环rpm 变量,从1000到3200,得到20个值:
python BEMT --para_file 'example_BEMT/vali2/tmotor28.ini' --rpm '20 1000.0 3200.0' --save_file 'example_BEMT/result' --draw_pic
输入参数格式:
输入参数 | 说明 |
---|---|
para_file | prop说明文件所在位置 |
V/rpm | 需要计算的参数,三个数字代表计算个数与计算的下限和上限 |
save_file | 保存的文件路径 |
draw_pic | 是否画图:True/False |
输入文件格式:
[case] rpm = 5000 v_inf = 1.0 [rotor] nblades = 2 diameter = 0.3 radius_hub = 0.0125 section = naca0012 naca0012 naca0012 naca0012 naca0012 naca0012 naca0012 naca0012 naca0012 naca0012 naca0012 naca0012 naca0012 radius = 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.11 0.12 0.13 0.14 0.15 chord = 0.027622 0.03166421 0.03390364 0.03458931 0.03399109 0.03238396 0.03003231 0.0271742 0.02400563 0.02066482 0.01721652 0.01363621 0.00979446 pitch = 43.83725288 39.0653117 32.64704522 28.10653926 24.81339996 21.99531739 19.5775673 17.72758179 16.3399819 15.05472676 13.7582959 12.87508345 12.1134445 load_path = F:\graduate_student\T2_GANpropeller\test1\0_database\foil\result_extrapolate_1_re5e5
[fluid] rho=1.225 mu=1.81e-5
- 外推文件位置需要在load_path 中输入,对应文件夹需要包含所有section中需要的外推文件。
结果输出:
- 在
save_file
参数对应文件夹内,会有csv
文件保存,文件名为对应运行时间- 第一行为输入V或者rpm的列表
- 下三行分别为CP,CT,eta对应工况下的值
- 可以使用excel画图或者-draw_pic 参数作图结果如下
cd example_BEMT/vali1
python run_prop.py
实验数据来自Definition of a benchmark for low Reynolds number propeller aeroacoustics。 Casalino, Damiano, et al. "Definition of a benchmark for low Reynolds number propeller aeroacoustics." Aerospace Science and Technology 113 (2021): 106707.
基础螺旋桨 | 截面翼型 | 转子直径 | 实验范围 | 轴半径 |
---|---|---|---|---|
APC-96 | NACA4412 | 0.3m | J-0~0.8 | 1.25cm |
验证结果:
误差值:
CP | CT | eta | |
---|---|---|---|
MSE | 0.0859 | 0.1373 | 0.0439 |
max relative err | 0.28 | 0.23 | 0.015 |
cd example_BEMT/vali2
python run_tmotor.py
结果: 实验数据来自水下涡轮实验 Bahaj A.S., et al. “Power and thrust measurements of marine current turbines under various hydrodynamic flow conditions in a cavitation tunnel and a towing tank.” Renewable Energy 32.3 (2007): 407-426.
本节介绍了根据python程序实现遗传算法NSGA-II自动优化旋翼的计算程序,配置环境与所需程序见resource3
文件夹,操作语句默认在此文件夹路径下进行。
遗传算法(Genetic Algorithm, GA)起源于对生物系统所进行的计算机模拟研究。它是模仿自然界生物进化机制发展起来的随机全局搜索和优化方法,借鉴了达尔文的进化论和孟德尔的遗传学说。其本质是一种高效、并行、全局搜索的方法,能在搜索过程中自动获取和积累有关搜索空间的知识,并自适应地控制搜索过程以求得最佳解。
概念名称 | 概念含义 |
---|---|
基因型(genotype) | 性状染色体的内部表现 |
表现型(phenotype) | 染色体决定的性状的外部表现,或者说,根据基因型形成的个体的外部表现 |
进化(evolution) | 种群逐渐适应生存环境,品质不断得到改良。生物的进化是以种群的形式进行 |
适应度(fitness) | 度量某个物种对于生存环境的适应程度 |
选择(selection) | 以一定的概率从种群中选择若干个个体。一般,选择过程是一种基于适应度的优胜劣汰的过程 |
复制(reproduction) | 细胞分裂时,遗传物质DNA通过复制而转移到新产生的细胞中,新细胞就继承了旧细胞的基因 |
交叉(crossover) | 两个染色体的某一相同位置处DNA被切断,前后两串分别交叉组合形成两个新的染色体。也称基因重组或杂交 |
变异(mutation) | 复制时可能(很小的概率)产生某些复制差错,变异产生新的染色体,表现出新的性状 |
编码(coding) | DNA中遗传信息在一个长链上按一定的模式排列。遗传编码可看作从表现型到基因型的映射 |
解码(decoding) | 基因型到表现型的映射 |
个体(individual) | 指染色体带有特征的实体 |
种群(population) | 个体的集合,该集合内个体数称为种群的大小 |
以生物学步骤为蓝本,诠释为如下过程
- 评估每条染色体对应个体表现型的适应度
- 遵循适应度越高,选择概率越大的原则,从种群中选择两个个体作为父方与母方
- 抽取父母双方的染色体,进行交叉产生子代
- 对子代的染色体进行变异
- 重复2、3、4过程,直到新种群产生
以算法形式表现进化过程
生物进化 | 遗传算法 |
---|---|
个体 | 问题的一个解 |
个体的竞争力 | 适应函数 |
适者生存 | 适应值最大的值被保留概率最大 |
染色体 | 解的编码 |
基因 | 编码元素 |
群体 | 被选定的一组解 |
种群 | 根据适应函数选择的一组解 |
交叉 | 以一定方式由双亲产生后代的过程 |
变异 | 编码的某些分量产生变化的过程 |
具体例子可以自主查询著名的组合优化问题「背包问题」。
数学公式(待补充)
略(待补充)
根据读取BEM文件.ini
的翼型确定截面数目以及初始翼型形状,
以旋翼prop.ini为例
以"Definition of a benchmark for low Reynolds number propeller aeroacoustics." 中旋翼为例 弦长分布: 扭转角分布 模型图 (待补充)
程序将自动转化为0-1的几何参数,作为优化参数
基因型,表现型
几何约束满足对应截面位置几何条件的最大、最小值,同时几何约束应满足生成旋翼的光滑程度
推力满足整体大于最小约束,局部区域小于最大约束
由第二节的BEMT程序计算旋翼属性。 在满足推力约束的要求下,效率最小为最优,
若推力不满足约束,应当推力最大为最优。
在局部区域的推力最大值应当为最优先(涉及材料属性)
对应位置的几何约束条件为最优先(涉及装配与材料问题)
给予一个基础BEMT计算旋翼文件设置基本工况与基础旋翼几何,给定设计目标与约束目标(力与几何约束)文件,程序自动基于该翼型优化出符合条件的旋翼。
设置一个优化约束文件,给定约束文件,以适应不同需求
python OPTIMIZE --BEM_file 'example_OPTIMIZE/prop.ini' --restrict 'example_OPTIMIZE/restrict.ini' --save_file 'example_OPTIMIZE/test/' --get_stl False --pic 'example_OPTIMIZE/test/OPTIMIZE_fig/'
BEM计算文件参考 example_OPTIMIZE/prop.ini
上一节有所说明
控制文件参考 example_OPTIMIZE/restrict1.ini
翼型文件格式参照example_extrapolate/ara_d_20.dat
参数说明:
输入参数 | 说明 |
---|---|
BEM_file | 输入为naca翼型或翼型坐标文件 |
restrict | 输入计算的雷诺数 |
save_file | 输出的优化base和优化过程变化的图片(用处较小) |
get_stl | 是否生成stl模型 |
pic | 输出的结果路径(最终结果) |
控制文件参数说明: 见参考文件后注释
模块:
- NSGA2算法优化
- 目标优化 eta(max)or FM
- 曲线光滑
- 整体约束 min_T
- 局部约束 local_cl,local_pitch,local_cR
约束目标eta最小,满足推力约束
分段多项式拟合?最高点位置 (7)
翼尖固定!