翼伞和常规降落伞一样，在充气展开过程中，由于本身的柔性特性，在气流作用下柔性伞衣极易变形，工作特性发生改变，涉及到的流固耦合问题具有高度非线性、非定常特点，理论研究难度很大。随着计算机的发展和数值精度的提高，数值模拟逐渐成为研究降落伞流固耦合性能的重要手段。Purvis[4]、潘星和曹义华[5]建立了降落伞结构的弹簧阻尼多节点模型，与内部简化流场进行耦合计算。基于时空离散的MSD-CFD(Mass Spring Damper-Computational Fluid Dynamics)耦合方法[6-8]侧重于伞衣外形的计算，无法得到结构应力信息。Takizawa[9]、Stein[10]和Tezduyar[11]等将DSD/SST(Deforming-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 FFD六面体控制晶格 Fig. 1 FFD hexahedral control lattice |
| 图选项 |
本文采用六面体控制晶格，如图 1所示，建立局部坐标系WSH，W、S、H是局部坐标系轴矢量，沿坐标轴分别均匀布置l、m、n个控制节点，局部坐标系原点在全局坐标系下的坐标记为Q0。
FFD控制晶格各节点的全局坐标为
 | (1) |
式中：(i, j, k)为控制晶格各节点的局部标识，i=1, 2, …, l；j=1, 2, …, m；k=1, 2, …, n。
在FFD控制晶格中，翼伞任意一点的局部坐标是(w, s, h)，0≤w, s, h≤1，有
 | (2) |
该点的全局坐标为Q=Q0+wW+sS+hH。
本文采用Bernstein基函数建立控制晶格各节点位移ΔPi, j, k与ΔQ之间的映射函数关系：
 | (3) |
式中：Bli、Bmj、Bnk分别为(i, j, k)节点的l、m、n次Bernstein基函数。当晶格节点改变时，根据ΔPi, j, k求出ΔQ，从而得到变形后的翼伞外形。
1.2 流固耦合控制方程 流场域控制方程为
 | (4) |
式中：ρF为流场密度；t为时间；vi为流体速度；ei为Euler坐标；ωi=Δei/Δt为网格节点速度；fi为体积力；E为能量；σijF为流场的应力张量，σijF=-pδij+μ(vi, j+vj, i)，p为压强，μ为动力黏度，δij为Kronecker δ函数。
流场域采用无反射边界条件，以消除流场边界处的单元速度反射波干扰，防止边界对计算产生影响。
翼伞伞衣、伞绳和流场发生耦合，采用拉格朗日方法计算，控制方程为
 | (5) |
式中：ρS为结构密度；xiS为结构节点位移；σij, jS为结构的应力张量。
ALE网格运动方程为
 | (6) |
式中：Li为拉格朗日坐标。
对控制方程式(4)~式(6)进行全耦合计算，采用中心差分时间显式方法求解。对于每一个节点，流场和结构的位移x、速度u按式(7)进行更新：
 | (7) |
式中：Fint、Fext分别为内、外力矢量；M为质量对角矩阵。
1.3 流场域时步更新 翼伞充气过程是有限质量充气，为减小计算消耗，本文在每一时间步都对流场域进行更新，使其跟随载荷对象一起进行空间运动。首先在结构网格上任意选择不共线的3个节点A、B、C，坐标分别为KA、KB和KC，建立一个以A为原点的局部坐标系(见图 2)，其3个轴矢量分别为
 | (8) |
 |
| 图 2 流场域时间步更新 Fig. 2 Time-step update of flow field |
| 图选项 |
每个时间步后，局部坐标将跟随结构网格发生位移，根据局部坐标在每个时间步前后的位移获得变换矩阵T，则流场网格新的坐标系为
 | (9) |
式中：[e1* e2* e3* 1]为更新后的齐次坐标；[e1 e2 e3 1]为更新前的齐次坐标。流体域通过时间步更新技术实现运动与重构。
2 研究对象与模型验证 2.1 研究对象与计算工况 本文采用文献[22]的冲压式翼伞模型，如图 3所示，翼梢两侧有稳定幅，32根伞绳，平铺展长5.48 m，弦长2.74 m，绳长3.36 m，有7个气室，每个气室由一个非承载肋片分割成2个半气室。翼伞的折叠模型如图 4所示。
 |
| 图 3 冲压式翼伞结构示意图 Fig. 3 Schematic diagram of rammed parafoil structure |
| 图选项 |
 |
| 图 4 翼伞折叠模型 Fig. 4 Folding model of parafoil |
| 图选项 |
本文采用ALE流固耦合数值方法，将结构网格单元穿插于流场网格中，其流固耦合数值计算模型如图 5所示(规模为6WS×8WS×10WS，其中WS为伞衣弦长)，其中结构附近流场进行加密处理，参数设置如表 1所示。翼伞以0°迎角充气，初始充气速度为17 m/s。本文基于ANSYS平台计算，采用8核处理器，仿真0.6 s耗时约172 h。
 |
| 图 5 流固耦合数值计算模型 Fig. 5 Numerical calculation model of fluid-structure interaction |
| 图选项 |
表 1 流固耦合过程的数值模型信息 Table 1 Numerical model information of fluid-structure interaction process
| 参数 | 流场 | 伞衣 | 伞绳 |
| 单元数量 | 806 868 | 28 885 | 1 749 |
| 单元类型 | Solid | Shell | Beam |
| 材料类型 | Null | Fabric | Cable_Discrete_Beam |
| 密度/(kg·m-3) | 1.2 | 533.8 | 462.0 |
| 弹性模量/Pa | 4.3×108 | 9.7×1010 |
表选项
2.2 模型验证 为验证本文模型的合理性和准确性，采用文献[22]的空投试验数据进行验证计算，同时进行流场网格的无关性验证。本文采用3种不同网格密度的流场进行数值计算，对比结果如表 2所示。
表 2 网格无关性对比结果 Table 2 Comparison results of grid independence
| 网格数量 | 最大开伞动载/N | 误差/% | 计算时长/h |
| 587 000 | 4 661.30 | 6.08 | 123 |
| 849 000 | 4 760.10 | 4.09 | 172 |
| 1 189 000 | 4 764.13 | 4.01 | 238 |
| 注：文献[22]试验数据中，最大开伞动载为4 963.2 N。表中误差指文献[22]与本文最大开伞动载的误差。 | |||
表选项
随网格数量增加，计算误差减小，但消耗时长增加，因此综合考虑，本文选择网格数量为849 000的模型进行数值计算，得到的翼伞充满外形和开伞动载F(充气过程中作用在伞衣上的纵向载荷)与空投试验的对比情况如图 6、图 7所示。可以看出，翼伞充满外形与试验外形基本一致。本文开伞动载计算结果的变化趋势与空投试验变化规律相近，数值计算得到的开伞动载峰值出现时刻略早于文献[22]，数值为4 760.1 N，和文献[22]试验数据(4 963.2 N)比较，误差为4.09%，该误差满足柔性翼伞开伞计算的误差要求。综上认为，本文基于多气室双翼面冲压翼伞建立的数值计算模型可以有效模拟翼伞充气阶段工作过程。
 |
| 图 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为翼伞充气过程的外形变化。在充气初期，翼伞前缘切口进气，尾部最先膨胀展开。之后下翼面高压气流作用于稳定幅，拉动伞衣展向伸展，稳定幅产生应力集中现象。充气过程中，前缘切口约束较少，且在气室充满后由于切口角度的存在，始终承受气流吹袭，震颤频率大，产生持续的应力集中现象。因此，在翼伞设计中应加强稳定幅与前缘切口处的材料强度，避免伞衣破损。
 |
| 图 8 充气过程翼伞外形变化(范式等效应力云图) Fig. 8 Parafoil shape changes during inflation process (von Mises equivalent stress contours) |
| 图选项 |
为进一步探究翼伞充气过程中的结构及流场的非定常变化规律，选取a、b、c 3个截面(见图 9)进行研究，各个截面的结构及流场变化如图 10所示(各个时刻从上到下依次为a、b、c截面)，为刻画细节情况，各个时刻的图形缩比尺寸不同。
 |
| 图 9 选取截面示意图 Fig. 9 Schematic diagram of selected sections |
| 图选项 |
首先，通过观察图 10所示翼伞充气过程的弦向变化(a、b截面)：在翼伞初始折叠状态时，a、b截面的速度场分布相似，前缘产生涡。随翼伞逐渐充气，上翼面的涡逐渐后移至脱体，前缘处不断产生新的涡。在充气过程中，a截面随气室充气逐渐形成饱满翼型，内腔压力增加，伞衣刚度增加，类似刚性翼，外部流畅稳定。而b截面处于最外侧，伞衣由于约束较少变形明显，在翼伞充满后翼型饱满度不足，外部流场不稳定，伞衣变形状态具有一定的随机性。
 |
| 图 10 各截面的结构及流场变化 Fig. 10 Structure and flow field changes of each section |
| 图选项 |
图 10中翼伞充气过程的展向变化(c截面)表明：在翼伞初始折叠状态时，c截面周围流场产生对称绕流现象，此时伞衣折叠压缩率大，下翼面形成高压区(压力云图见图 11)。之后翼伞底部对称涡开始破裂，气流绕翼尖流出，翼伞展开充气。该过程上翼面流速快，下翼面流速慢，上下翼面形成的压力造成“翼尖上翘，中部凹陷”的翼伞尾流再附现象。因此在翼伞开伞过程中应加入引导伞改善充气状况，引导伞拉出伞衣中部减轻气室塌陷。此外，设计合理的稳定幅，阻挡下翼面到上翼面的绕流，消弱翼尖涡强度，减轻干扰影响。此后，伞衣随气室充满保持稳定充满翼型。
 |
| 图 11 充气过程压力云图 Fig. 11 Pressure contours of inflation process |
| 图选项 |
3.2 气室充气规律 为进一步研究翼伞的各气室充气规律，将翼伞沿弦向分为前(前缘)、中(中部)、后(尾缘)三部分，沿展向将7个气室逐次编号，两侧为气室1和气室7，中央气室为4号。图 12为各气室充满时间和充满宽度的变化规律。可以看出：
 |
| 图 12 各气室充满时间和充满宽度 Fig. 12 Inflation time and width of each air chamber |
| 图选项 |
1) 各气室的充气规律关于中央气室具有对称性，中央气室充气快且充满外形饱满，越靠近两侧充满速率越慢，充气效果越差。
2) 在高压气流推动作用下，进入气室的气体涌入尾部，翼伞尾缘最先展向打开，充气速率明显快于其他部分，前缘展开最慢。充满后中部充气效果最好，气室最为饱满。
3) 由于尾缘上下翼面缝合，翼尖绕流导致伞衣变形，内外气压不稳定，两侧气室的后部未完全充满。
充气过程的翼展变化见图 13。翼伞的前缘、中部、尾缘前、中、后三部分于0.39、0.37、0.25 s分别达到翼展峰值5.05、5.40、5.20 m。之后由于伞绳和伞衣的弹性作用，各部位充满后产生一定的回弹和气室“鼓包”现象，中部充满翼展稳定在5.0 m左右，略小于设计值。
 |
| 图 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。
 |
| 图 14 气动特性参数的变化 Fig. 14 Changes in aerodynamic characteristic parameter |
| 图选项 |
4 结论 本文针对非定常充气展开过程，数值模拟了翼伞的三维流固耦合动力学特性，得出：
1) 翼伞充气展开过程中稳定幅和前缘切口处产生应力集中现象，因此，在翼伞设计中应加强稳定幅与前缘切口处的材料强度。
2) 翼伞充气展开过程存在“翼尖上翘，中部凹陷”的翼伞尾流再附现象，造成伞衣内塌，导致系统稳定性受影响，因此在翼伞开伞过程中应加入引导伞以改善充气状况，设计合理的稳定幅可以减小翼尖涡的影响。
3) 各气室的充气规律关于中央气室具有对称性，中央气室充气快且充满外形饱满，两侧充满速率越慢，气室饱满度下降；翼伞尾部最先膨胀展开，前缘展开最慢，充满后存在一定的回弹和气室“鼓包”现象，充满后展长小于设计值。
4) 升阻力特性参数变化趋势一致，伞衣充满后翼伞滑翔比稳定在2.24，升力系数为1.56，阻力系数为0.70。
参考文献
| [1] | PATEL S, HACKETT N, JORGENSEN D, et al.Qualification of the guided parafoil air delivery system-light(GPADS-Light): AIAA-97-1493[R].Reston: AIAA, 1997. |
| [2] | BENNETT T, FOX R.Design, development & flight testing of the NASA X-387500 ft2 parafoil recovery system: AIAA-2003-2107[R].Reston: AIAA, 2003. |
| [3] | SHAW D O.A parafoil-based, hybrid airship design for extended Martian exploration[C]//AIAA SPACE Forum.Reston: AIAA, 2016: 1-7. |
| [4] | PURVIS J W. Theoretical analysis of parachute inflation including fluid kinetics[J]. Journal of Aircraft, 1982, 19(4): 290-296. DOI:10.2514/3.57392 |
| [5] | 潘星, 曹义华. 降落伞开伞过程的多结点模型仿真[J]. 北京航空航天大学学报, 2008, 34(2): 188-192. PAN X, CAO Y H. Simulation of parachute's opening process with multi-node model[J]. Journal of Beijing University of Aeronautics and Astronautics, 2008, 34(2): 188-192. (in Chinese) |
| [6] | STEEVES E C, BENNEY R J, STEIN K R.A computational model that couples aerodynamic and structural dynamic behavior of parachutes during the opening process: NASA-ADA264115[R].Washingtong, D.C.: NASA, 1993. |
| [7] | YU L, MING X. Study on transient aerodynamic characteristics of parachute opening process[J]. Acta Mechanica Sinica, 2007, 23(6): 627-633. DOI:10.1007/s10409-007-0112-3 |
| [8] | 彭勇, 张青斌, 秦子增. 降落伞主充气阶段数值模拟[J]. 国防科技大学学报, 2004, 26(2): 13-16. PENG Y, ZHANG Q B, QIN Z Z. Simulation of parachute final inflation phase[J]. Journal of National University of Defense Technology, 2004, 26(2): 13-16. DOI:10.3969/j.issn.1001-2486.2004.02.004 (in Chinese) |
| [9] | TAKIZAWA K, SPIELMAN T, MOORMAN C, et al. Fluid-structure interaction modeling of spacecraft parachutes for simulation-based design[J]. Journal of Applied Mechanics, 2012, 79(1): 010907. DOI:10.1115/1.4005070 |
| [10] | STEIN K, BENNEY R, KALRO V, et al. Parachute fluid-structure interactions:3-D computation[J]. Computer Methods in Applied Mechanics and Engineering, 2000, 190(3): 373-386. |
| [11] | TEZDUYAR T, BEHR M, MITTAL S, et al. A new strategy for finite element computations involving moving boundaries and interfaces-The deforming-spatial-domain/space-time procedure:Ⅱ.Computation of free-surface flows, two-liquid flows, and flows with drifting cylinders[J]. Computer Methods in Applied Mechanics and Engineering, 1992, 94(3): 353-371. DOI:10.1016/0045-7825(92)90060-W |
| [12] | KIM Y, PESKIN C S. 2-D parachute simulation by the immersed boundary method[J]. SIAM Journal on Scientific Computing, 2006, 28(6): 2294-2312. DOI:10.1137/S1064827501389060 |
| [13] | TUTT B, TAYLOR A.The use of LS-DYNA to simulate the inflation of a parachute canopy: AIAA-2005-1608[R].Reston: AIAA, 2005. |
| [14] | YU L, CHENG H, ZHAN Y N, et al. Study of parachute inflation process using fluid-structure interaction method[J]. Chinese Journal of Aeronautics, 2014, 27(2): 272-279. DOI:10.1016/j.cja.2014.02.021 |
| [15] | 陆伟伟, 张红英, 连亮. 大型翼伞的三维气动性能分析[J]. 航天返回与遥感, 2015, 36(3): 1-10. LU W W, ZHANG H Y, LIAN L. A three-dimensional analysis on aerodynamic performance of a large parafoil[J]. Spacecraft Recovery & Remote Sensing, 2015, 36(3): 1-10. DOI:10.3969/j.issn.1009-8518.2015.03.001 (in Chinese) |
| [16] | 朱旭, 曹义华. 翼伞平面形状对翼伞气动性能的影响[J]. 航空学报, 2011, 32(11): 1998-2007. ZHU X, CAO Y H. Numerical simulation of platform geometry effect on parafoil aerodynamic performance[J]. Acta Aeronautica et Astronautica Sinica, 2011, 32(11): 1998-2007. (in Chinese) |
| [17] | FOGELL N A, SHERWIN S, COTTER C J, et al.Fluid-structure interaction simulation of the inflated shape of ram-air parachutes: AIAA-2013-1326[R].Reston: AIAA, 2013. |
| [18] | KALRO V, ALIABADI S, GARRARD W, et al. Parallel finite element simulation of large ram-air parachutes[J]. International Journal for Numerical Methods in Fluids, 1997, 24(12): 1353-1369. DOI:10.1002/(SICI)1097-0363(199706)24:12<1353::AID-FLD564>3.0.CO;2-6 |
| [19] | 汪龙芳, 贺卫亮. 基于索膜有限元模型的翼伞气动变形仿真[J]. 北京航空航天大学学报, 2017, 43(1): 47-52. WANG L F, HE W L. Parafoil aerodynamic deformation simulation based on cable-membrane finite element model[J]. Journal of Beijing University of Aeronautics and Astronautics, 2017, 43(1): 47-52. (in Chinese) |
| [20] | 张春, 曹义华. 基于弱耦合的翼伞气动变形数值模拟[J]. 北京航空航天大学学报, 2013, 39(5): 605-609. ZHANG C, CAO Y H. Numerical simulation of parafoil aerodynamics and structural deformation based on loose coupled method[J]. Journal of Beijing University of Aeronautics and Astronautics, 2013, 39(5): 605-609. (in Chinese) |
| [21] | SEDERBERG T W, PARRY S R. Free-form deformation of solid geometric models[J]. ACM SIGGRAPH Computer Graphics, 1986, 20(4): 151-160. DOI:10.1145/15886.15903 |
| [22] | LEE C K, BUCKLEY J E. New technique for parafoil inflation control[J]. Journal of Aircraft, 2000, 37(3): 479-483. DOI:10.2514/2.2622 |
