临近空间太阳能飞行器着陆状态阵风响应

2024-03-08吴金洋吴光辉魏小辉房兴波

海军航空大学学报 2024年1期

吴金洋，吴光辉，魏小辉，聂宏，房兴波，陈虎

（1.南京航空航天大学航空学院，江苏南京 210016；2.中国商用飞机有限责任公司，上海 200126；3.南京航空航天大学通用航空与飞行学院，江苏溧阳 213300）

0 引言

临近空间太阳能飞行器通常巡航飞行在20～100 km 的临近空间区域。与常规布局飞行器相比，该类飞行器往往具有尺寸大、质量轻、翼展长、机体柔性大等几何非线性结构特点，以及高升力、低阻力等气动性能，且着陆过程中飞行速度极低。这些因素极易在风切变、湍流等突发气流载荷条件下使机翼产生大的弯曲和扭转形变，引起机翼气动载荷分布的改变，进而对整个飞行器的气动特性产生较大影响，甚至导致飞行器解体坠毁事故。2016 年，美国Facebook 公司的Aquila 太阳能无人机在降落初期由于遭遇阵风载荷，最后导致坠毁解体。基于此，研究临近空间太阳能飞行器起降过程中的气动弹性阵风响应问题对其安全性具有重要意义。

针对临近空间太阳能飞行器的阵风响应问题，国内外研究人员开展了大量的研究工作，如使用非线性梁结构和大迎角气动力模型[1-2]，准模态法和非定常涡格法[3-4]，准模态法和修正Theodorsen 方法结合片条理论[5]等进行结构与气动建模，探讨了大柔性飞行器在不同阵风载荷下结构和气动特性的变化。近年来，基于欧拉方程或N-S 方程的CFD 与CSD 相互耦合解决非线性气动弹性问题的精确计算发展迅速[6]。在研究几何非线性大柔性飞行器阵风响应时，CFD/CSD耦合展现出良好的性能，能够准确反映飞行器在大位移及结构大变形条件下的动态特性[7]。

为了对飞行器阵风响应进行理论分析并设计相应的减缓控制系统，阵风模型也是1 个关键因素。早期研究者主要采用离散阵风模型[8]，如“锐边”阵风和“斜坡”阵风，它们可描述大气紊流中的峰值、飞行器尾流区的流动以及地形诱导的气流等[9]。随着阵风研究的深入，“锐边”和“斜坡”阵风逐渐被1-cos 型离散阵风取代，该类型阵风至今仍在使用，并被FAR-25、CCAR-25 等适航法规作为分析飞行器阵风载荷的理论模型。由于空间尺度上非均匀阵风的复杂性，上述阵风模型都基于一维均匀阵风场假设进行了简化，即认为飞行器阵风速度仅在飞行方向上变化，沿展向方向保持不变。然而，对于拥有较大展向尺寸的大展弦比柔性飞行器，采用一维均匀阵风场模型进行分析存在缺陷，需要进一步考虑展向阵风变化[10-11]。

为此，本文基于CFD/CSD 松耦合分析法，利用网格速度法引入阵风载荷，研究了二维阵风对临近空间太阳能飞行器着陆状态时阵风响应的影响，并将二维阵风响应与一维阵风响应进行对比。研究获得了临近空间太阳能飞行器着陆时，二维阵风沿飞行器翼展方向的变化对其结构和气动性能的影响规律。

1 临近空间太阳能飞行器模型

参考Facebook Aquila 无人机的设计，本文所计算的临近空间太阳能飞行器模型如图1所示。

图1 临近空间太阳能飞行器模型Fig.1 Model of near-space solar-powered aircraft

该太阳能飞行器采用内、外翼飞翼式布局；飞行器的展长为48.00 m，平均气动弦长为1.47 m，展弦比为32.65，总质量为280 kg。飞行器所用材料及参数见文献[12]。利用Ansys 软件进行模型构建。建模时，针对着陆状态，飞行器模型有一定迎角，且来流速度很低，空气密度接近海平面。鉴于计算对象具有对称特性，为提升计算效率，仅对一侧流场进行分析，为此选择半实体模型。为降低流体网格数量，构建流体域时使前后平面与飞行器前缘平行。借助Ansys Mechanical 进行流体和固体网格划分，对流固耦合面网格进行细化处理。图2展示了模型部分网格示意图。

图2 CFD局部网格示意图Fig.2 Schematic diagram of CFD local grid

2 计算方法

尽管所研究的临近空间太阳能飞行器巡航阶段低雷诺数效应显著，但其在起降过程中雷诺数超过百万，属于高雷诺数范围。据此，采用CFD/CSD 松耦合方法进行求解。在CFD 部分，利用基于SST 湍流模型的有限体积法对可压缩积分形式的雷诺平均Navier-Stokes（RANS）方程组进行计算。在空间离散过程中：对流通量项采用Roe 格式的二阶精度迎风格式；黏性通量项选用二阶精度的中心离散方式；时间推进则使用LU-SGS隐式时间推进技术。

2.1 气动模型

SST 湍流模型是将标准k-ω模型与标准k-ε模型通过混合函数相结合而形成[13]。其在继承k-ε模型的稳定性、计算时间少以及k-ω模型无需壁面函数即可进行壁面积分等优势的基础上，剔除了k-ε模型对逆压力梯度和边界层分离不敏感以及k-ω模型收敛困难等不足。因此，该模型在近壁区域和自由剪切层中展现出优良的数值特性，并且能够较好地模拟具有较大逆压梯度特征的流动。SSTk-ω湍流模型的通用控制方程为：

式（1）（2）中：t为计算时间（飞行时间）；ρ为密度；k为湍流动能；ω为湍流特殊耗散；uˉi、uˉj为湍流速度平均值；xi、xj为坐标分量；Γk、Γω为有效扩散系数；Fk、Fω为湍流生成项；Yk、Yω为k、ω的耗散项；Dω为扩散项。

2.2 网格速度法

阵风模型使用网格速度法引入，流体区域内的流场速度可表示为：

式（3）中：u、v和w为直角坐标系中的速度分量；xt、yt和zt分别为3个方向上的网格速度分量；i、j和k代表直角坐标系在各个方向上的单位向量。

当流体区域受到1 个向上、速度为wi的阵风影响时，z方向上的流体区域网格会出现1个速度为wi的突增，而实际上流体域网格未发生运动。此时，流体域相当于在网格不动的情况下整体都受到-wi的来流作用。据此，在阵风作用下，流体域中流场的速度可以表示为：

2.3 阵风计算模型

一维阵风模型采用1-cos工程阵风模型[14]，整个飞行器在竖直方向上受到如式（5）所示的阵风影响。

式（5）中：wg为阵风速度；Ude为阵风幅值；Lg为阵风长度；V为飞行器速度；t为计算时间（飞行时间）。

图3 为飞行器所遭遇的1-cos 型阵风速度型[15]，设飞行器平均气动弦长为c，无穷远处来流速度为V∞，定义无量纲时间S=2V∞t c。

图3 一维1-cos型阵风速度型[15]Fig.3 One-dimension 1-cos gust speed pattern

以只有竖直方向速度的1-cos阵风模型为基础，二维阵风模型考虑了阵风幅值沿展向的变化。因此，二维阵风速度型为[16-17]：

式（6）中：Lspan为飞行器翼展；a为阵风沿展向分布的因子；y为沿展向坐标。

图4 为a=1.0 时的阵风模型；图5 为不同a值下阵风沿展向的分布。由图可看出，随着a的增大，二维阵风向一维阵风逼近，当a接近无穷大时，二维阵风退化为一维阵风。

图4 二维阵风模型[17]Fig.4 Two-dimensional gust model

图5 二维阵风的展向分布[17]Fig.5 Two-dimensional gust spanwise distribution

2.4 计算流程

本研究中，即便降落时飞行器未遭遇阵风，机翼本身也处在一定变形下的平衡状态。为此，可建立如图6 所示的CFD/CSD 松耦合方式的计算流程。首先，通过CFD/CSD 松耦合迭代收敛求出飞行器降落时的静平衡状态；然后，在平衡后某一时刻添加阵风载荷；最后，通过CFD/CSD松耦合迭代，直到收敛为止。

图6 松耦合计算流程图Fig.6 Flowchart of loose coupling calculation

3 计算结果及分析

3.1 阵风维数的影响

取飞行器飞行高度为500 m，飞行速度为40 km/h（11.1 m/s），大气密度为1.167 kg/m3，降落时初始迎角为4°，先算出飞行器的静平衡状态，然后以此为基础，计算到10 s 时分别施加一维和二维阵风载荷，研究不同阵风维数对飞行器着陆时阵风响应的影响。其中，阵风幅值Ude1为6 m/s，阵风长度Lg为25 倍飞行器平均气动弦长[18]，二维阵风沿展向分布因子设为1.0。参照图4，在-0.25～0.25 飞行器翼展的跨度范围内，阵风速度方向竖直向上，而其余跨度范围内的阵风速度方向竖直向下。

图7 给出了在一维和二维阵风载荷作用下，飞行器达到静平衡状态后，受到阵风载荷时翼尖位移和扭转角的对比。可看出：在一维和二维阵风载荷作用下，翼尖变形均出现较大的响应峰值，但翼尖z向位移响应的峰值方向相反，且二维阵风载荷作用下得到的峰值比一维阵风载荷作用下的小，出现时间相对滞后；此外，达到峰值后，一维阵风载荷作用下得到的翼尖位移和扭转角恢复到平衡位置，而二维阵风载荷作用下得到的翼尖z向位移未恢复到静平衡状态，只是恢复到-0.09 m，翼尖扭转角则减小到-1.99°。

图7 翼尖位移和扭转角阵风响应计算结果对比Fig.7 Gust response calculation results comparison for wingtip displacement and torsional angle

图8 给出了在一维和二维阵风载荷作用下，飞行器达到静平衡状态后，受到阵风载荷时翼根所受的结构整体弯矩的对比。

图8 翼根所受结构弯矩阵风响应计算结果对比Fig.8 Gust response calculation results comparison for structural moment on the wing root

可看出：二维阵风载荷作用下得到的结构弯矩响应结果峰值比一维阵风载荷作用下的小，且方向相反，出现时间相对滞后；此外，二维阵风载荷作用下得到的弯矩在达到峰值后未恢复到静平衡状态，而是振荡恢复到-900 N·m左右，导致二维阵风载荷作用下得到的翼尖变形减小到负值。这是由于展向分布因子为1.0 时，翼尖部分所受的二维阵风速度方向竖直向下，与一维阵风速度方向相反所致。

图9 给出了在一维和二维阵风载荷作用下，飞行器达到静平衡状态后，受到阵风载荷时整体升力系数的对比。

图9 飞行器整体升力系数阵风响应计算结果对比Fig.9 Gust response calculation results comparison for lift coefficient of aircraft

可看出：二维阵风载荷作用下得到的结构整体升力系数响应峰值比一维阵风载荷作用下的小，且出现时间相对滞后；此外，二维阵风载荷作用下得到的整体升力系数在达到峰值后未恢复到静平衡状态，而是振荡恢复到0.43 左右。这是由于飞行器受到展向分布因子为1.0 的二维阵风载荷后，整体未恢复到静平衡状态，而是出现反方向扭转，使得飞行器整体迎角减小所致。

3.2 展向分布因子的影响

图10给出了飞行器达到静平衡状态后，展向分布因子分别取0.5、1.0、2.0 和5.0 时，二维阵风载荷作用下翼尖z向位移、扭转角、翼根所受结构弯矩和整体升力系数的对比。结合图7～9 可看出：随着展向分布因子的增加，飞行器在二维阵风载荷作用下的阵风响应逐渐接近于一维阵风作用下的阵风响应。当a≥2.0时，响应峰值随a的增加逐渐增大，接近受到一维阵风作用下的响应峰值；尤其当a=5.0 时，二维阵风的响应峰值和响应趋势与一维阵风的具有较好的一致性；而当a ＜2.0时，由于阵风存在速度方向竖直向下的部分，飞行器的阵风响应过程与受一维阵风载荷作用下的响应过程差异较大；当a=0.5时，所有的响应结果峰值较小，除扭转角响应结果出现2 个峰值外，翼尖位移、翼根结构弯矩以及升力系数响应结果峰值方向均与受一维阵风载荷作用下的方向相反；而当a=1.0时，所有的响应结果达到峰值后未恢复到平衡状态，且翼尖z向位移和翼根所受结构弯矩的响应结果的峰值方向均与受一维阵风载荷作用下的响应方向相反，并达到最大。