翼伞充气过程的流固耦合方法数值仿真

2020-07-31张思宇余莉刘鑫

北京航空航天大学学报 2020年6期

张思宇，余莉，，*，刘鑫

（1.南京航空航天大学飞行器环境控制与生命保障工业和信息化部重点实验室，南京210016；2.南京航空航天大学航空学院，南京210016）

冲压式翼伞相比于传统降落伞，升阻比高、滑翔稳定性能和操作性良好，广泛应用于航空航天等领域［1-3］，近年来成为各国研究重点。翼伞充气成功与否决定着翼伞系统工作的成败。

翼伞和常规降落伞一样，在充气展开过程中，由于本身的柔性特性，在气流作用下柔性伞衣极易变形，工作特性发生改变，涉及到的流固耦合问题具有高度非线性、非定常特点，理论研究难度很大。随着计算机的发展和数值精度的提高，数值模拟逐渐成为研究降落伞流固耦合性能的重要手段。Purvis［4］、潘星和曹义华［5］建立了降落伞结构的弹簧阻尼多节点模型，与内部简化流场进行耦合计算。基于时空离散的MSD-CFD（Mass Spring Damper-Computational Fluid Dynamics）耦合方法［6-8］侧重于伞衣外形的计算，无法得到结构应力信息。Takizawa［9］、Stein［10］和Tezduyar［11］等将DSD/SST（Deform ing-Spatial-Domain/Stabilized Space-Time）有限元计算方法用于三维降落伞流固耦合仿真计算研究中，该方法侧重于外部流场的数值计算，对初始充气段缺少研究。Kim和Peskin［12］采用浸入边界法进行降落伞充气过程的流固耦合模拟。Tutt［13］和Yu［14］等采用任意拉格朗日-欧拉（Arbitrary Lagrange-Euler，ALE）方法模拟了三维降落伞开伞过程。

但是，目前降落伞的流固耦合研究大多集中于平面圆形伞，针对翼伞充气过程的三维流固耦合数值研究一直没有突破性进展，常常做出一定简化，主要面向刚性翼伞［15-16］或者基于松散耦合方法进行稳态定常计算。例如，Fogell等［17］采用弱耦合方法分析了翼伞在稳定滑翔阶段的气动外形。Kalro等［18］采用并行有限元计算对翼伞初始充气阶段进行了数值模拟。汪龙芳和贺卫亮［19］对定常状况下翼伞的流固耦合变形问题进行了三维数值模拟。张春和曹义华［20］基于弱耦合方法分析了翼伞的气动变形。

翼伞充气过程大变形、大位移、多尺度流动的非定常流固耦合研究基本处于空白阶段。这主要是由于：翼伞为非展平曲面构成的多气室-双翼面复杂结构，之间由开孔翼肋隔成若干连通气室，前缘充气，后缘封闭，结构的复杂性导致其无法采用常规方法建立折叠模型。此外，多气室翼伞在充气过程中多方向同时展开，是一个典型的非稳态流场-柔性伞衣结构的剧烈耦合问题，规律复杂。上述问题作为翼伞研究领域的难点，一直备受学者关注。

为了解决上述问题，本文提出了基于自由曲面变形（Free Form Deformation，FFD）理论的多气室柔性冲压翼伞展向折叠建模方法，流体域通过时步更新实现随伞载系统运动，基于ALE方法开展折叠翼伞充气过程的非定常流固耦合数值计算，分析翼伞的非线性动力学行为，系统地探究非稳态流场-柔性伞衣结构之间的耦合机理，为翼伞安全设计提供理论依据。

1 计算模型

1.1 自由曲面变形理论

翼伞为非展平曲面所构成的多气室-双翼面结构，为解决其展向折叠问题，本文根据FFD理论［21］，提出了多气室柔性冲压翼伞的展向折叠建模方法。将待变形的翼伞模型嵌入FFD控制晶格中，通过移动晶格节点，将晶格变形传递给嵌入其中的翼伞模型，达到对其折叠建模的目的。

本文采用六面体控制晶格，如图1所示，建立局部坐标系WSH，W、S、H是局部坐标系轴矢量，沿坐标轴分别均匀布置l、m、n个控制节点，局部坐标系原点在全局坐标系下的坐标记为Q0。

FFD控制晶格各节点的全局坐标为

图1 FFD六面体控制晶格Fig.1 FFD hexahedral control lattice

式中：（i，j，k）为控制晶格各节点的局部标识，i＝1，2，…，l；j＝1，2，…，m；k＝1，2，…，n。

在FFD控制晶格中，翼伞任意一点的局部坐标是（w，s，h），0≤w，s，h≤1，有

该点的全局坐标为Q＝Q0+w W+s S+h H。本文采用Bernstein基函数建立控制晶格各节点位移ΔPi，j，k与ΔQ之间的映射函数关系：

1.2 流固耦合控制方程

流场域控制方程为

流场域采用无反射边界条件，以消除流场边界处的单元速度反射波干扰，防止边界对计算产生影响。

翼伞伞衣、伞绳和流场发生耦合，采用拉格朗日方法计算，控制方程为

ALE网格运动方程为

式中：Li为拉格朗日坐标。

对控制方程式（4）～式（6）进行全耦合计算，采用中心差分时间显式方法求解。对于每一个节点，流场和结构的位移x、速度u按式（7）进行更新：

式中：Fint、Fext分别为内、外力矢量；M 为质量对角矩阵。

1.3 流场域时步更新

翼伞充气过程是有限质量充气，为减小计算消耗，本文在每一时间步都对流场域进行更新，使其跟随载荷对象一起进行空间运动。首先在结构网格上任意选择不共线的3个节点A、B、C，坐标分别为KA、KB和KC，建立一个以A为原点的局部坐标系（见图2），其3个轴矢量分别为

图2 流场域时间步更新Fig.2 Time-step update of flow field

每个时间步后，局部坐标将跟随结构网格发生位移，根据局部坐标在每个时间步前后的位移获得变换矩阵T，则流场网格新的坐标系为

式中：［e*1e*2e*31］为更新后的齐次坐标；［e1e2e31］为更新前的齐次坐标。流体域通过时间步更新技术实现运动与重构。

2 研究对象与模型验证

2.1 研究对象与计算工况

本文采用文献［22］的冲压式翼伞模型，如图3所示，翼梢两侧有稳定幅，32根伞绳，平铺展长5.48m，弦长2.74m，绳长3.36m，有7个气室，每个气室由一个非承载肋片分割成2个半气室。翼伞的折叠模型如图4所示。

本文采用ALE流固耦合数值方法，将结构网格单元穿插于流场网格中，其流固耦合数值计算模型如图5所示（规模为6WS×8WS×10WS，其中WS为伞衣弦长），其中结构附近流场进行加密处理，参数设置如表1所示。翼伞以0°迎角充气，初始充气速度为17m/s。本文基于ANSYS平台计算，采用8核处理器，仿真0.6 s耗时约172 h。

图3 冲压式翼伞结构示意图Fig.3 Schematic diagram of rammed parafoil structure

图5 流固耦合数值计算模型Fig.5 Numerical calculation model of fluid-structure interaction

表1 流固耦合过程的数值模型信息Table 1 Numerical model information of fluid-structure interaction process

2.2 模型验证

为验证本文模型的合理性和准确性，采用文献［22］的空投试验数据进行验证计算，同时进行流场网格的无关性验证。本文采用3种不同网格密度的流场进行数值计算，对比结果如表2所示。

随网格数量增加，计算误差减小，但消耗时长增加，因此综合考虑，本文选择网格数量为849 000的模型进行数值计算，得到的翼伞充满外形和开伞动载F（充气过程中作用在伞衣上的纵向载荷）与空投试验的对比情况如图6、图7所示。可以看出，翼伞充满外形与试验外形基本一致。本文开伞动载计算结果的变化趋势与空投试验变化规律相近，数值计算得到的开伞动载峰值出现时刻略早于文献［22］，数值为4 760.1 N，和文献［22］试验数据（4 963.2 N）比较，误差为4.09%，该误差满足柔性翼伞开伞计算的误差要求。综上认为，本文基于多气室双翼面冲压翼伞建立的数值计算模型可以有效模拟翼伞充气阶段工作过程。

表2 网格无关性对比结果Table 2 Comparison results of grid independence

图6 数值计算结果与空投试验充满外形对比Fig.6 Comparison of canopy shape between numerical calculation result and airdrop test

图7 开伞动载随时间的变化曲线Fig.7 Curves of opening load over time

3 研究结果

3.1 伞衣外形及流场变化

翼伞充气展开过程为有限质量充气情况，许多随机参数都会对充气性能产生影响，本文未考虑风场等随机环境因素对充气过程的影响。图8为翼伞充气过程的外形变化。在充气初期，翼伞前缘切口进气，尾部最先膨胀展开。之后下翼面高压气流作用于稳定幅，拉动伞衣展向伸展，稳定幅产生应力集中现象。充气过程中，前缘切口约束较少，且在气室充满后由于切口角度的存在，始终承受气流吹袭，震颤频率大，产生持续的应力集中现象。因此，在翼伞设计中应加强稳定幅与前缘切口处的材料强度，避免伞衣破损。

为进一步探究翼伞充气过程中的结构及流场的非定常变化规律，选取a、b、c 3个截面（见图9）进行研究，各个截面的结构及流场变化如图10所示（各个时刻从上到下依次为a、b、c截面），为刻画细节情况，各个时刻的图形缩比尺寸不同。

图8 充气过程翼伞外形变化（范式等效应力云图）Fig.8 Parafoil shape changes during inflation process（von M ises equivalent stress contours）

图9 选取截面示意图Fig.9 Schematic diagram of selected sections

首先，通过观察图10所示翼伞充气过程的弦向变化（a、b截面）：在翼伞初始折叠状态时，a、b截面的速度场分布相似，前缘产生涡。随翼伞逐渐充气，上翼面的涡逐渐后移至脱体，前缘处不断产生新的涡。在充气过程中，a截面随气室充气逐渐形成饱满翼型，内腔压力增加，伞衣刚度增加，类似刚性翼，外部流畅稳定。而b截面处于最外侧，伞衣由于约束较少变形明显，在翼伞充满后翼型饱满度不足，外部流场不稳定，伞衣变形状态具有一定的随机性。

图10中翼伞充气过程的展向变化（c截面）表明：在翼伞初始折叠状态时，c截面周围流场产生对称绕流现象，此时伞衣折叠压缩率大，下翼面形成高压区（压力云图见图11）。之后翼伞底部对称涡开始破裂，气流绕翼尖流出，翼伞展开充气。该过程上翼面流速快，下翼面流速慢，上下翼面形成的压力造成“翼尖上翘，中部凹陷”的翼伞尾流再附现象。因此在翼伞开伞过程中应加入引导伞改善充气状况，引导伞拉出伞衣中部减轻气室塌陷。此外，设计合理的稳定幅，阻挡下翼面到上翼面的绕流，消弱翼尖涡强度，减轻干扰影响。此后，伞衣随气室充满保持稳定充满翼型。

图11 充气过程压力云图Fig.11 Pressure contours of inflation process

3.2 气室充气规律

为进一步研究翼伞的各气室充气规律，将翼伞沿弦向分为前（前缘）、中（中部）、后（尾缘）三部分，沿展向将7个气室逐次编号，两侧为气室1和气室7，中央气室为4号。图12为各气室充满时间和充满宽度的变化规律。可以看出：

1）各气室的充气规律关于中央气室具有对称性，中央气室充气快且充满外形饱满，越靠近两侧充满速率越慢，充气效果越差。

2）在高压气流推动作用下，进入气室的气体涌入尾部，翼伞尾缘最先展向打开，充气速率明显快于其他部分，前缘展开最慢。充满后中部充气效果最好，气室最为饱满。

图12 各气室充满时间和充满宽度Fig.12 Inflation time and width of each air chamber

3）由于尾缘上下翼面缝合，翼尖绕流导致伞衣变形，内外气压不稳定，两侧气室的后部未完全充满。

充气过程的翼展变化见图13。翼伞的前缘、中部、尾缘前、中、后三部分于0.39、0.37、0.25 s分别达到翼展峰值5.05、5.40、5.20m。之后由于伞绳和伞衣的弹性作用，各部位充满后产生一定的回弹和气室“鼓包”现象，中部充满翼展稳定在5.0m左右，略小于设计值。

图13 充气过程的翼展变化Fig.13 Changes in wingspan during inflation process

3.3 气动特性分析

翼伞的气动特性参数变化情况如图14所示。0.1 s之前，计算域尚未形成稳定绕流流场，气动参数不具有参考意义。随绕流流场形成，翼伞的升阻力特性参数变化趋势一致，在0.28 s时达到峰值，此时伞衣投影面积最大。之后由于呼吸回弹现象，伞衣面积减小，气动系数减小，翼伞充满稳定后气动系数略有回升。充满后伞载系统存在一定迎角，滑翔比CL／CD稳定在2.24，升力系数CL为1.56，阻力系数CD为0.70。