耦合分布螺旋桨的大柔性机翼静气弹研究

吴轩,周洲,王正平

(西北工业大学 航空学院,陕西 西安 710072)

面对未来航空器高经济性、低能耗性、低噪声、低污染的“绿色航空发展要求”,采用分布式螺旋桨电推进系统的高空长航时无人机正被应用于各种场景[1]。此类飞机普遍具有大展弦比、低翼载、大柔性的特点,在飞行载荷作用下,机翼会产生较大的几何非线性变形,传统的线性气动弹性分析手段已经不再适用。同时,虽然单个螺旋桨相对于大展弦比机翼尺度较小,但沿展向分布的多个螺旋桨,其载荷和滑流都会对机翼的流场、载荷分布、变形分布有较大的干扰[2]。

针对大柔性机翼的几何非线性气弹问题,Patil[3]采用本征梁模型和有限状态入流模型系统研究了机翼结构的几何非线性对结构刚度、气动力、气动弹性的影响。Palacios等[4]采用应变梁模型与计算流体力学(CFD)相耦合的方法,研究了机翼的非线性静气弹响应问题。张驰[5]使用共旋转梁模型对非线性结构进行建模,在有限元的总体和单元坐标系中重新定义了飞行力学的地轴系、体轴系和风轴系,研究了耦合飞行力学下的气弹响应问题。针对螺旋桨对机翼气动弹性的影响研究,Hodges等[6]将螺旋桨推力定义为具有常值的集中力,研究了推力对机翼弯扭颤振的影响,发现推力对弯扭颤振同时具有促进和抑制作用。Agostinelli等[7]采用叶素理论和CFD方法,并结合线性有限元模型,分析了螺旋桨载荷、滑流诱导速度对弹性机翼的气动干扰。Teixeira等[8]使用黏性涡粒子法将螺旋桨滑流影响和陀螺效应引入到柔性机翼的气弹分析中,指出螺旋桨滑流增加了机翼变形且陀螺效应的影响可以忽略。谢长川等[9]使用动量理论和兰金涡核模拟了螺旋桨滑流对柔性机翼的气动影响,结果表明旋转的滑流会改变机翼当地绕流攻角,从而影响机翼气动力分布和变形分布。

从现有研究现状来看,对于柔性机翼在分布式螺旋桨影响下的非线性气弹响应问题,传统研究通常仅考虑螺旋桨的拉力,并将其作为集中力处理,并未考虑滑流的影响。文献[10]指出对于柔性飞机,未考虑螺旋桨滑流影响是仿真结果存在误差的重要原因之一。而对于螺旋桨滑流对机翼的影响,现有研究更注重二者之间的气动干扰问题,忽略了柔性机翼的几何非线性变形特征。同时,螺旋桨分布位置所带来的气动弹性影响也亟待研究。

对此,本文建立了一种考虑分布式螺旋桨载荷和滑流影响的柔性机翼非线性静气弹分析方法。首先,利用单个螺旋桨相对于大展弦比机翼尺度较小、转速较大的特性,采用滑流管模型模拟滑流对机翼的诱导速度,避免了螺旋桨复杂的尾迹计算。螺旋桨环量分布通过涡流叶素理论获取,机翼气动力则使用曲面涡格法进行计算。其次,在共旋转法的节点坐标系中,推导出展向不同位置处螺旋桨的局部柱坐标系,完成分布式螺旋桨与机翼结构在非线性大变形下的耦合建模。随后,对螺旋桨与机翼耦合的构型进行气动升力计算,将结果与现有方法对比,验证了所提气弹模型的准确性。最后,通过分布式螺旋桨与柔性机翼耦合的算例,对比分析了有无螺旋桨载荷和滑流以及螺旋桨分布位置对大柔性机翼静气弹响应的影响,再次表明在几何非线性气弹模型中引入分布式螺旋桨影响的重要性。

1 螺旋桨/机翼气动模型

1.1 涡流叶素理论

采用结合涡流模型的片条理论对小负荷螺旋桨进行气动建模,假定螺旋桨叶片与机翼一样,沿桨叶展向分布着环量变化的附着涡系Γ(r),在展向每一位置处均有自由涡溢出,并在后缘附近形成一个扭曲的涡面。对于每个单元中的叶素,其半径为r,弦长为b,受力和速度分解图如图1所示。

图1 叶素受力和速度分解图

在实际飞行中,叶素运动轨迹为螺旋线,前飞速度为V1,旋转角速度为Ω,桨盘平面内的切向速度V2=Ωr,则气流相对于叶素的几何合成速度Vgeo为

(1)

Vgeo与旋转平面之间的夹角即为初始入流角φ0

(2)

叶素受轴向诱导速度Va、切向诱导速度Vt影响,其实际入流速度Vreal为

(3)

经诱导后,气流的实际入流角φ为

(4)

令叶素的安装角为θ(叶素弦线方向与旋转平面之间的夹角),尾涡诱导产生的干涉角为β,则气流相对于叶素的当地迎角α为

α=θ-β-φ0

(5)

在迎角α下,当地升阻力系数和马赫数及雷诺数有关,即Cl=Cl(α,Re,Ma),Cd=Cd(α,Re,Ma)。根据涡流理论,经桨尖环量损失修正[11]后的当地叶素环量为

(6)

式中

式中:R为螺旋桨半径;NB为螺旋桨桨叶数目。

根据库塔-茹科夫斯基升力定理,当地环量分布与升力系数关系为

(7)

联立求解(6)～(7)式,可得螺旋桨环量分布Γ(r),进一步可以得到该单元处的拉力和扭矩

(8)

将单元处拉力dT和扭矩dM沿径向和周向积分,即可得到螺旋桨总拉力和总扭矩

(9)

式中,Rc为螺旋桨轮毂半径。

1.2 螺旋桨滑流管尾迹模型

对于本文研究的带有分布式螺旋桨的大展弦比机翼来说,机翼飞行速度远远小于螺旋桨旋转速度,因此可以忽略螺旋桨的非定常效应,同时假设机翼附近的螺旋桨尾迹不发生变形,从而避免了复杂的尾迹更新过程,提高气弹分析计算效率。

螺旋桨滑流对机翼的气动影响主要分为2个部分:①滑流的轴向诱导速度分量使得机翼当地流动速度增加、动压增大;②滑流的切向诱导速度分量引发局部上洗或者下洗流动。因此,如图2所示,本文采用Veldhuis[12]提出的滑流管模型,将螺旋桨尾涡分成轴向涡系和切向涡系2个部分。

图2 螺旋桨尾迹涡系模型

使用毕奥萨伐尔定律计算微元涡段对空间任一点Q(xQ,yQ,zQ)的诱导速度

(10)

式中,dVQ为微元涡段对Q点的诱导速度;Γ表示涡段环量大小;S表示点Q到涡段的向量;dl表示微元涡段向量。

图2中的轴向涡和切向涡可以表示为

(11)

式中:Γ(r)为桨叶的环量分布,由1.1节中给出;NB为桨叶数目;Ω表示螺旋桨旋转速度;U∞表示自由来流速度。

在螺旋桨局部柱坐标系中,对(10)式沿桨盘周向φ和纵向z积分,并结合(11)式,最终可得到切向涡对Q点的诱导速度为

(12)

轴向涡对Q点的诱导速度为

(13)

式中:b=rsinφ-xQ,c=-rcosφ-yQ。

桨盘处的诱导速度是个小量,可以忽略,因此滑流管模型对Q点的最终诱导速度为

(14)

1.3 机翼曲面涡格法气动力计算

鉴于螺旋桨与机翼的相对尺寸较小且螺旋桨旋转速度远高于机翼的飞行速度,故忽略机翼对螺旋桨的诱导速度。同时,螺旋桨之间距离较远,故忽略螺旋桨之间的气动干扰,着重关注螺旋桨载荷对柔性机翼气动的影响。

气动载荷使柔性机翼的升力面发生弯曲扭转变形从而形成空间曲面,使用曲面涡格法对其进行气动力计算。如图3所示,将气动面沿展向和弦向划分网格,共计M个单元。在每个气动单元内布置涡环基本解,控制点位于涡格1/4弦向处,此布置形式自然满足Kelvin条件和流体边界条件。

图3 考虑滑流影响的机翼曲面涡格

在当前构型下,升力面附着涡满足物面不可穿透条件,即在每个涡格的控制点处物面法方向nj上的诱导速度分量为零,即满足

(15)

2 螺旋桨/机翼结构模型

2.1 机翼几何非线性建模

为了耦合螺旋桨模型,采用基于位移变量的共旋转法对柔性机翼结构进行几何非线性分析,其物理意义更加明确。如图4所示,共旋转法通过局部节点坐标系和单元坐标系来描述结构的刚体运动和变形位形。

图4 共旋转法单元与节点坐标系

(16)

式中,R(α)是广义旋转位移对应的旋转矩阵。根据Rodrigues对旋转矩阵的参数化描述方法,可以得到旋转矩阵和旋转伪矢量的关系式

(17)

式中:α=|α|为旋转伪矢量的模;S(α)为斜对称矩阵

R(γ)=UTT

(18)

(19)

(20)

所有推导的节点坐标系和单元坐标系都是结构位形的随动坐标系,随结构变形实时更新。

根据(20)式,三维梁在局部单元坐标系下的位移为pl,总体坐标系下的位移为p:

(21)

对(21)式进行微分处理,即可得到总体坐标系下的位移增量δp和单元坐标系下的位移增量δpl的转换关系为

δpl=Fδp

(22)

式中,F为转换矩阵。

令梁单元在总体坐标系下的内力为q,在单元坐标系下的内力为ql,根据虚功原理,有

(23)

在单元坐标系下,共旋转法认为力与位移的平衡关系是线弹性的,同时根据(22)式,可以得到

q=FTql=FTKlpl

(24)

式中,Kl是线性刚度矩阵。

对(24)式进行变分处理即可得到增量平衡方程和对应的切向刚度矩阵

(25)

式中:FTKlF为材料刚度矩阵;Ktσ(ql)为几何刚度矩阵且为单元局部内力ql的函数,KT为共旋转法的切线刚度矩阵,具体表达式见文献[5]。

2.2 螺旋桨机翼耦合建模

上文1.2节中推导的滑流诱导速度是在螺旋桨局部柱坐标系下进行的,然而由于机翼的大变形,螺旋桨的位置和方向随机翼结构变形实时变化,需要考虑在结构模型中推导螺旋桨局部柱坐标系的表达形式,完成惯性坐标系、结构单元节点坐标系和螺旋桨柱坐标系的统一以及载荷变形等数据传递。

如图5所示,对于一般安装在机翼前缘的螺旋桨推进系统,在惯性系下,螺旋桨局部柱坐标系原点位于中心点Op处,平面极坐标rp(rp,φp)位于桨叶旋转平面,zp方向垂直于旋转平面与推力线方向OTP反向。在基于共旋转法的机翼非线性结构模型中,螺旋桨安装位置处机翼节点坐标系为Ui(Uix,Uiy,Uiz),节点坐标为Oi。

图5 螺旋桨柱坐标系与结构节点坐标系

(26)

(27)

忽略推进系统旋转产生的回转力矩效应,螺旋桨的推力Tp和转矩Mp在1.1节中可以通过涡流叶素理论获得。则在惯性系中,螺旋桨安装位置处结构单元节点所受到的力Fi和力矩Mi为

(28)

3 螺旋桨/机翼非线性静气弹分析

3.1 空间梁样条插值

采用空间梁样条完成在大柔性变形下,机翼结构节点与气动节点之间的位移与载荷插值。在共旋转梁模型中,结构节点平动位移和转动位移分别为Xs和Ψs。令气动节点到对应截面结构节点的矢量为ξa,假设机翼在变形过程中,截面保持刚性,则气动节点位移Xa可以通过结构节点位移插值得到

Xa=Xs+Us(Ψs)ξa=GXs

(29)

式中,Us为截面转动位移对应的转换矩阵;G为位移插值矩阵,具体表达式见文献[14]。

设结构界面的虚位移为δXs,气动界面的虚位移为δXa,气动界面的载荷为Fa,结构界面的载荷为Fs,根据虚功原理有

(30)

根据(29)式,可以得到气动结构节点的位移增量关系式

δXa=GδXs

(31)

将(31)式代入(30)式,得

(32)

根据虚位移的任意性,可以得到气动载荷与结构载荷之间的插值关系式为

Fs=GTFa

(33)

3.2 非线性静气弹分析框架

大柔性飞机从本质上来说是一个分布式的参数系统,全机不同位置处结构位移、气动响应呈现非线性变化。考虑柔性飞机的特点以及分布式螺旋桨的影响,本文建立如图6所示的非线性静气弹分析框架并完成相应各个模块的程序编写。

图6 非线性静气弹分析框架

整个分析框架分为2个部分:螺旋桨机翼气动计算模块和非线性结构计算模块。在气动计算模块,首先通过2.2节中非线性结构与螺旋桨的耦合建模,可以确定不同位置处螺旋桨的实时随动位形,并对空间气动曲面进行涡格划分。其次,通过1.1节中介绍的涡流叶素理论可以获取各个螺旋桨的载荷和环量分布,并结合1.2节和1.3节中的滑流管模型和曲面涡格法,完成考虑滑流影响的大柔性机翼气动特性计算。在结构计算模块,首先通过空间梁样条将螺旋桨载荷和气动载荷转换加载到结构有限元节点上。同时,利用迭代计算的方法求解2.1节中的增量平衡公式(25),直至结构达到变形收敛状态,最终完成考虑分布式螺旋桨影响的柔性机翼非线性静气弹分析并输出结果。

4 方法验证及算例分析

4.1 气弹模型验证

采用密歇根大学的高空长航时验证无人机[15](X-HALE)为研究对象,该机翼分为内翼段和外翼段,展长分别为2 m,1 m,弦长为0.2 m。机翼在翼根处有5°的安装角,外翼段上反10°,其他特性参数见表1。如图7所示,5个螺旋桨分布在中心点以及机翼半展长的±1/3,±2/3位置处,螺旋桨距离机翼前缘前伸0.2 m、下移0.028 m,相对翼根安装角为0°,螺旋桨型号为APC11×5.5E,参数见文献[16]。

图7 X-HALE机翼示意图

设定自由来流速度v∞=14 m/s,来流攻角为2°,机翼翼根处固支。设定螺旋桨转速为6 000 r/min,其拉力为1.70 N。根据收敛性研究,机翼弦向划分12个涡格,展向划分96个涡格。根据3.2节中的静气弹分析框架,计算机翼在刚性不变性和柔性大变形下的收敛位移和升力分布,计算结果如图8所示。

图8 X-HALE机翼位移与升力分布

图8a)为机翼受来流气动载荷和螺旋桨载荷共同作用下的垂向位移分布图,翼尖最大位移约占半展长的35%。本文结果和文献[15]结果基本一致,静变形误差基本可以忽略,验证了本文非线性静气弹分析模型的准确性。

图8b)为机翼的展向升力系数分布图,可以看出,滑流轴向诱导速度可以引起翼面加速流动从而使得局部升力变大,切向诱导速度在机翼附近引发上洗和下洗流动也会改变机翼的展向升力分布。机翼的整体升力系数如表2所示,变形后机翼的升力系数相对未变形构型降低17%左右,本文和文献[15]的结果相对误差在2%以内。

表2 机翼升力系数对比

升力分布以及升力系数的计算结果表明,本文采用的气动模型以及编写的程序是合理准确的,可以有效获取螺旋桨对机翼的气动影响。同时,计算时间相对于CFD等高阶方法大大降低,仅需30 s左右,本文气动模型可应用于此类飞机初始设计阶段的快速分析。

4.2 柔性机翼耦合分布式螺旋桨算例分析

选用一广泛被用作基准模型的大展弦比柔性机翼[17]作为研究对象,使用3.2节中的求解框架探究螺旋桨的载荷、滑流、分布位置对大柔性机翼的静气弹影响。图9所示为机翼半模示意图,10个螺旋桨沿展向对称布置在机翼半展长的3/8,4/8,5/8,6/8,7/8位置处。所有螺旋桨弦向距离前缘0.5m,垂向安装距离为0,桨轴与机翼共面,转向为正(机翼外侧为上洗区域)。

图9 耦合分布式螺旋桨的柔性机翼(半模)示意图

螺旋桨具体几何参数见文献[18],机翼的具体参数见表3。自由来流速度为25 m/s,空气密度为0.088 9 kg/m3,所有螺旋桨转速皆为4 400 r/min,机翼翼根处固支,安装角为4°。

表3 大柔性机翼参数

基于收敛性研究,机翼结构有限元模型共分为32个梁单元,机翼升力面弦向划分12个单元,展向划分128个单元,结构和气动的计算收敛误差皆在1%以内。

4.2.1 螺旋桨对结构变形的影响

机翼非线性弯曲、扭转变形计算结果对比如图10～11所示,在无螺旋桨的情况下,本文计算得到的垂向弯曲变形和Murua等[17]使用SHARP程序计算的结果基本一致,翼尖最大垂向位移误差为1.2%,表明本文气弹模型和代码具有一定的准确性。

图10 柔性机翼垂向弯曲变形图

图11 柔性机翼扭转变形图

在只考虑螺旋桨拉力扭矩载荷不考虑滑流影响的情况下,翼尖最大弯曲变形相对于干净机翼减小9.6%,翼尖最大扭转变形减小34.3%。这主要是由于较大的弯曲变形使螺旋桨拉力产生机翼结构负扭转,从而造成机翼局部迎角减小、升力减小,进一步使弯曲变形减小。

当引入滑流的影响时,翼尖最大弯曲变形相较于只考虑螺旋桨载荷的情况增加3.2%,翼尖最大扭转变形减小7.9%。滑流带来机翼局部入流速度增加,从而引起升力增加、机翼弯曲变形增加。同时,滑流引起气动俯仰力矩减小,产生低头力矩,从而减小扭转变形。

4.2.2 螺旋桨对气动特性的影响

图12给出了在不同迎角下,螺旋桨载荷和滑流对机翼升力系数曲线的影响。由于拉力产生结构负扭转,机翼局部迎角变小,升力系数相较于干净机翼降低10%左右。滑流使机翼局部来流速度增加,导致机翼升力系数相对于只考虑螺旋桨载荷的情形增加约2.5%左右。

图12 升力系数随迎角变化曲线

图13为机翼纵向静稳定裕度Cmα/Clα随迎角变化曲线,可以看出螺旋桨的载荷使静稳定裕度绝对值减小20%～35%,滑流使静稳定裕度绝对值增加2%～8%,差值随迎角增大而增大。

图13 静稳定裕度随迎角变化曲线

4.2.3 螺旋桨不同分布位置的影响

分别只保留机翼半展长3/8,4/8,5/8,6/8,7/8位置处对称的2个螺旋桨工作,关闭其他8个螺旋桨,计算条件不变,图14～15给出了螺旋桨位于不同位置处机翼升力系数和俯仰力矩系数随转速变化曲线。

图14 螺旋桨位于不同位置处升力系数随转速曲线

如图14所示,当螺旋桨的位置靠近翼根时,滑流的影响更大,导致升力增加;靠近翼尖时,拉力的影响更大,导致升力减小;当螺旋桨位于半展长的3/8到4/8之间时,载荷和滑流对升力系数的影响相互抵消。从曲线斜率可知,越靠近翼尖,转速对升力线斜率Clα影响越大。如图15所示,螺旋桨在各个位置均增大了俯仰力矩系数,造成机翼纵向静稳定裕度减小。同样,越靠近翼尖,转速对力矩系数影响越大。

图15 螺旋桨位于不同位置处俯仰力矩系数随转速曲线

5 结 论

1) 本文在基于涡流理论的滑流管模型以及共旋转法非线性结构模型的基础上,建立了考虑分布式螺旋桨载荷和滑流影响的大柔性机翼静气弹快速分析求解方法。将机翼升力分布和静气弹变形等结果与文献[15]对比,验证了本文模型具有一定的准确性,能够有效体现大柔性机翼的非线性变形特征以及分布式螺旋桨与机翼的耦合气动特性。

2) 分布式螺旋桨的载荷和滑流会对柔性机翼的弯曲扭转变形、气动力分布、静稳定裕度等产生较大的影响,使机翼的非线性效应更加显著。螺旋桨拉力会产生机翼结构负扭转,使当地有效迎角减小,造成气动升力损失和静稳定裕度降低;滑流则通过增加当地流速和改变绕流攻角,增加了升力和静稳定裕度。在此类飞机的初始设计阶段,必须考虑这些影响,以免造成设计偏差。

3) 螺旋桨在展向不同位置对大柔性机翼的性能影响不同,靠近翼根时增升,靠近翼尖时减升,且越靠近翼尖螺旋桨转速变化带来的气弹影响越显著,文献[19]也指出了分布式螺旋桨可以控制机翼升力分布以及抑制机翼弹性变形的潜力。在后续研究中,可以根据本文模型方法和分析结果,通过调节螺旋桨的转速、安装位置、分布形式等参数完成大柔性机翼与分布式螺旋桨的耦合优化设计。