连红博 喻东明 周昊 赵颖涛

摘 要：为研究十字型阻力伞的开伞特性， 利用ALE（Arbitrary Lagrange-Euler）流-固耦合算法， 对十字型阻力伞的充气过程进行了数值模拟研究。 本文讨论了计算域尺寸与网格密度对计算结果的影响， 并在此基础上详细分析了伞衣的外形变化及开伞动载情况。 结果显示， 不同开伞速度下， 最大动载的模拟结果与实验相吻合， 误差均在5%以内， 验证了模拟方法的可行性与准确性; 伞衣充气时首先由伞臂展开， 之后气流在伞衣顶部聚集， 待伞顶被完全撑开后， 开始逐渐向外扩展， 直至伞衣完全张满; 阻力伞的开伞动载受伞衣外形变化影响显著， 且当伞衣完全张满时， 阻力伞达到最大开伞动载。

关键词：     阻力伞; 流-固耦合; 数值模拟; 充气过程; 流场模型; 开伞动载

中图分类号：     TJ760; V211.3  文献标识码：    A 文章编号：     1673-5048（2021）06-0083-05

0 引  言

十字型阻力伞属于降落伞的范畴， 因具有开伞快、 稳定性好、 结构简单且动载小等特点， 常作为飞机着陆滑跑时的减速装置[1]。 当前评估阻力伞开伞性能的试验方法主要包括风洞试验、 水洞试验及火箭橇试验等， 这些试验方法的测试成本较高， 在实际情况下还存在开伞速度有限且难以精确控制等不足[2-4]。 因此， 数值模拟方法以其成本低、 周期短、 参数更改灵活等优点， 成为研究阻力伞充气展开过程的重要手段。

阻力伞的充气展开变形大且历时短。 在开伞过程中， 受伞衣初始状态、 透气性及自接触等因素影响， 周围的流场变化复杂， 属于典型的流-固耦合问题[5-6]。 近年来， ALE（Arbitrary Lagrange-Euler）算法被越来越多地应用于阻力伞开伞过程的模拟， 因为该算法综合了Lagrange算法与Euler算法的優点， 在准确描述耦合边界的同时， 避免了网格畸变的影响。

利用ALE算法， Tutt和Aquelet[7-8]针对TP8降落伞进行了有无透气性的对比分析， 发现伞衣材料的透气性对提升降落伞充气过程的稳定性具有显著效果。 国内针对阻力伞充气过程的数值研究起步较晚， 且所研究的伞型相对单一。  房明[9]利用LS-DYNA软件对远场边界下单、 双十字型伞的充气展开过程进行了仿真， 并通过对比验证了该算法与MSD（Mass Spring Damper）算法结果一致。  杨品等[10]采用最新的S-ALE（Structured ALE）算法模拟了飞机阻力伞从伞舱抛出至完全打开的全过程， 详细分析了阻力伞的开伞动载以及飞机的减速情况。  喻东明等[11]基于十字型阻力伞充满时间的实测值推算出了阻力伞的充满时间系数， 并利用数值仿真验证了其圆整结果的准确性。

目前， 国内外对于十字型阻力伞充气过程的仿真研究仍比较缺乏， 本文基于ALE流-固耦合算法对十字型阻力伞充气过程展开数值模拟， 详细给出了充气时伞衣结构形态及不同速度下开伞动载的变化情况。

1 基于ALE的流-固耦合方法

1.1 ALE方法介绍

如何描述伞衣与空气的受力及运动是模拟阻力伞充气过程的关键所在， 所以算法的选择十分重要。 目前有限元软件常用的算法有Lagrange、 Euler以及ALE算法， 三者之间的主要区别如表1所示。 ALE算法于1964年由Noh等[12]首次提出， 其是对Lagrange算法与Euler算法的综合考量。 同Euler算法类似， 其网格与物质材料相互独立， 不同的是其网格在计算时可以发生任意运动， 这不仅解决了网格的畸变问题， 也大大减小了物体发生大位移时ALE网格的计算域。

阻力伞充气过程涉及伞衣大变形与流体运动， 所以本文选择Lagrange算法来划分伞衣网格， ALE算法来描述空气域， 并将二者耦合实现阻力伞开伞过程的仿真。

问题—

连红博， 等： 十字型阻力伞充气过程的数值模拟

1.2 控制方程

阻力伞开伞过程中运动速度相对较小， 可将周围流体视作不可压缩。 采用ALE算法对阻力伞开伞过程进行流-固耦合模拟时， 周围流体遵循的控制方程仍然是Navier-Stokes方程[9]， 即

ρt=-ρvixi-ωiρxiρvit=σij，j+ρbi-ρωivixjρEt=σijvi，j+ρbivi-ρωjExjσij=-pδij+μvi，j+vj，i （1）

式中： ρ，vi，σij，p，E，bi，μ分别表示流体密度、 物质速度、 应力张量、 压力、 能量、 单位体积力、 动态粘性系数; ωi是相对速度， ωi=vi-ui， ui是网格速度; σij，j=σij/xj。

网格控制方程为[10]

fXi，tt=fxi，tt+wifxi，txi（2）

式中： Xi为拉格朗日网格坐标。

伞的结构控制方程为[13]

M+C+Kw=F（3）

式中： M， C， K分别指单元质量、 阻尼模量和弹性模量; F为膜单元所受合力。

1.3 流-固耦合及伞衣透气性算法

阻力伞开伞流-固耦合的重点在于伞衣与空气耦合面上节点信息的传递， 本文选用LS-DYNA进行数值解算， 其主要通过罚函数方法实现节点信息的传递。 罚函数方法会对Lagrange节点与Euler流体界面的位置进行检查， 当二者发生贯穿时， 会产生一对大小相等、 方向相反的耦合力。 耦合力的大小F=kd， 其中k是刚度系数， d为二者之间发生的相对位移。

阻力伞本身为织物材料， 所以在对其进行仿真时， 需充分考虑伞衣透气性对计算结果的影响[14]。 本文采用Ergun型多孔介质模型来描述阻力伞的透气性能， 其计算公式为[15]

Δpe=a（ε， μ）v+b（ε， ρ）v2（4）

式中： Δp，  e， v分别代表压力差、 伞衣厚度、 穿透速度; a（ε， μ）为粘性系数， 取值与空气动态粘性系数μ与伞衣孔隙率ε相关; b（ε，ρ）代表惯性系数， 可通过空气密度ρ以及伞衣孔隙率ε求得。

2 阻力伞模型

2.1 基本假设

伞衣的充气展开是一个十分复杂的过程。 为简化计算， 在不影响物理属性的情况下对仿真进行以下假设：

a. 伞衣初始为简单折叠的轴对称结构;

b. 伞绳充气前为拉直状态， 且与伞衣的初始应力为0;

c. 忽略重力影响。

2.2 阻力伞折叠模型

本文所研究阻力伞为十字型伞， 其平铺模型如图1所示。 其主要由主伞衣、 伞绳和连接带组成， 具体结构参数见表2。 其中， 主傘衣上还布有横向、 径向及对角加强带， 为伞衣打开时提供有力支撑。 伞臂上特意留有的伞缝和伞顶孔具有良好的排气效果。 32根伞绳与加强带直接相连， 在汇聚至连接带后挂在飞机连接环上， 从而起到阻力减速作用。 在仿真过程中， 整个模型取为轴对称结构， 伞衣采用壳单元， 网格尺寸约为61 mm， 加强带、 伞绳与连接带均采用一维安全带单元。 对于两种不同单元的重合处， 主要是通过共节点的方式进行连接。

图1为阻力伞的平铺状态， 而实际阻力伞放伞前， 通常是折叠状态。 为了更好地模拟阻力伞的开伞过程， 本文对平铺模型进行了简单折叠。 折叠过程参考文献[16]的方式，  文中不再赘述。

如图2所示。 首先将伞衣模型的1/4拆成两半分别进行“Z”折; 然后再把它们两个折在一起， 并且合并公共边上的节点; 利用已折叠的1/4模型， 便可通过镜像操作获得整个折叠伞衣模型; 最后建立伞绳、 连接带模型。 在伞衣折叠的过程中， 伞衣尺寸会发生一定的变化， 这会给阻力伞开伞后的外形带来影响， 所以需要对折叠处的单元进行微调。 同时对部分安全带单元进行松弛设置， 以使其在恢复至原始长度前不产生拉力作用。

3 流场无关性验证

阻力伞在被飞机抛出拉直后， 会在空气阻力的作用下快速充气展开， 因此流场的建立对于分析阻力伞的开伞过程十分关键。 考虑到计算域大小、 单元数量以及网格密度等因素的影响， 本文对不同情况下阻力伞开伞的计算结果进行对比分析， 以选择合适的仿真条件。

3.1 不同尺寸的计算域

在实际应用中， 该十字型伞稳定张满后的最大宽度一般在4.4～4.6 m之间， 故以张满宽度D=4.5 m为基础建立不同尺寸的流场计算域， 如图2所示。 其中， 保持计算域横向尺寸为6D， 并将垂直于横向的截面尺寸分别定为9D×9D， 7D×7D， 5D×5D和3D×3D。 此外， 对计算域中心与阻力伞的耦合区域进行加密， 四种不同计算域中加密区的伞衣与流场网格尺寸比均为1∶1。

流场入口速度设为250 km/h， 边界条件取为无反射边界， 得到不同尺寸计算域下阻力伞的开伞载荷-时间曲线， 如图3所示。 由图可知， 计算域的大小对计算结果和计算精度存在显著影响。 整体来看， 随着计算域的增大， 开伞达到的载荷峰值大致呈减小趋势。 但当计算域截面宽度大于5D后， 最大载荷已基本稳定在62.5 kN左右。 对比该型阻力伞试验数据， 发现其开伞速度为250 km/h时的最大载荷在61～67 kN之间， 认为数值计算结果与试验结果吻合良好[10]。 综合考虑计算精度和计算成本， 本文选用Mesh_7D作为后续仿真的计算域。

3.2 不同网格密度

保持计算域尺寸为Mesh_7D， 将加密区伞衣与流场网格尺寸比分别设为1∶1， 1∶1.5， 1∶2.5及1∶4， 如图4所示。 仿真得到对应的开伞载荷-时间曲线， 如图5所示。

从图5可以看出， 相比于Mesh_1∶1， 在伞衣单元尺寸与加密区单元尺寸出现较大偏差时， 开伞载荷曲线随时间的波动幅度变大， 且开伞最大载荷出现在伞衣顶部被完全撑开的瞬间， 这与实际不符。 不同的流场单元尺寸改变了伞衣与流体的耦合节点及二者之间的相对位移， 从而影响到整个开伞过程流-固耦合的计算精度。 故本文选取Mesh_1∶1的网格密度作为后续计算条件。

基于计算域尺寸及网格密度的综合考量， 本文选择计算域大小为6D×7D×7D、 伞衣与加密区网格尺寸比为1∶1的流场网格， 展开阻力伞充气过程的分析与计算。

4 结果及讨论

4.1 模拟结果准确性验证

利用以上仿真模型， 本文得到不同来流速度下开伞载荷-时间变化曲线， 如图6所示。 可以看出， 三种情况下， 伞衣所受载荷随时间的变化趋势基本一致， 开伞最大动载随来流速度的增大而增加， 且速度的增大也在一定程度上缩短了开伞完成时间。

为进一步验证数值模拟结果的准确性， 将阻力载荷的计算结果和该型阻力伞放伞过程的试验结果进行对比。 阻力伞载荷的工程计算公式为[17]

F=12ρv2CsAs（5）

式中： F表示伞衣阻力载荷; ρ为空气密度; v为开伞速度; Cs是阻力系数; As为伞衣面积。

对该型阻力伞进行不同开伞速度下的实物试验， 得到对应的最大开伞动载如图7所示。

对试验数据进行二次多项式拟合， 可得

F=0.001 012v2（6）

当来流速度为250 km/h， 270 km/h及300 km/h时， 仿真得到的最大开伞载荷分别为62.1 kN， 73.5 kN和89.2 kN。 由图可知， 3次数值模拟结果都在拟合曲线附近， 且误差均不超过5%， 验证了模拟方法的准确性。

4.2 阻力伞开伞外形分析

阻力伞初始为折叠状态， 为研究其充气过程中伞衣随时间的变化情况， 本文以来流速度250 km/h为例， 选取7个关键时间点来描述阻力伞的整个开伞过程， 如图8所示。

由图看出， 阻力伞在初始充气阶段， 首先自伞臂处展开， 外形变化缓慢; 从0.128 s至0.176 s， 空气向阻力伞顶部聚集， 该阶段伞衣外形变化剧烈， 直到气流将伞顶完全撑开; 从0.176 s至0.344 s， 聚集在伞顶的气流开始向四周扩散， 伞臂逐渐向外张开， 伞衣外形不断充满， 并在0.344 s时初步达到完全张满状态; 从0.344 s至0.528 s， 伞衣外形经历了多次调整后趋于稳定。

图9（a）为阻力伞开伞仿真得到的稳定后的伞衣张满外形（t=0.8 s）， 图9（b）为基于火箭橇试验平台获得的阻力伞实际张满状态[4]， 二者外形十分吻合。

4.3 阻力伞开伞动载分析

图10为自由来流250 km/h时阻力伞的开伞载荷-时间曲线。 结合图8分析可知， 在阻力傘初始充气阶段， 由于伞衣首先从伞臂张开， 伞衣阻力面积变化幅度不大，

开伞载荷增长缓慢; 而在t=0.128 s时， 开伞载荷开始急剧增大， 这是因为伞衣顶部开始在气流作用下展开， 伞衣外形变化剧烈; 直到0.176 s， 伞衣顶部被完全顶开， 此时开伞载荷达到约40.7 kN; 接下来一段时间， 由于前期伞衣展开太快而使得载荷发生了回落; 至t=0.2 s时， 伞衣内气流从顶部向径向流动， 使得伞衣不断向外张开， 开伞载荷也重新开始逐渐增大， 并于0.344 s时达到最大开伞力为62.1 kN， 计算结果和试验结果61～67 kN相符; 之后开伞载荷经历若干次波动， 幅度不断减小， 并在0.528 s时达到基本稳定状态， 稳定载荷约为49.2 kN。

5 结  论

（1） 提出一种采用ALE流-固耦合技术对十字型阻力伞充气过程进行数值模拟的方法，分析了计算域大小、 网格密度对计算精度的影响， 确定了适用的数值仿真条件， 并验证不同速度下最大动载的数值模拟结果与试验值拟合结果的误差均在5%以内。

（2） 十字型阻力伞在充气时首先由伞臂展开， 之后气流在伞衣顶部聚集， 使其膨胀， 然后膨胀部分向伞臂扩展， 直至伞衣完全张满。

（3） 阻力伞的开伞动载受伞衣外形的影响显著。 当伞衣完全张满时， 阻力伞达到最大开伞动载; 空气流速越快， 开伞动载越高， 所需要的开伞时间也越短。

参考文献：

[1] 韩雅慧， 杨春信， 肖华军， 等. 十字形伞的实验研究进展及展望[J]. 兵工自动化， 2013， 32（3）： 3-8.

Han Yahui， Yang Chunxin， Xiao Huajun， et al. Development History and Expectation of Cross Parachute[J]. Ordnance Industry Automation， 2013， 32（3）： 3-8.（in Chinese）

[2] Brocato B， Esteve L， Garcia D， et al. Experimental Study of Fluid-Structure Interactions on a Cross Parachute-Comparison of Wind Tunnel Data and Drop Data with CFD Predictions[C]∥15th Aerodynamic Decelerator Systems Technology Conference， 1999.

[3] Siefers T， Greene K， McLaughlin T， et al. Wind and Water Tunnel Measurements of Parachute Suspension Line[C]∥51st AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aero-space Exposition， 2013： 1015-1028.

[4] 余元元， 沈文波， 李磊子. 基于火箭橇试验的阻力伞试验平台研究[C]∥中国空气动力学会测控技术专委会第六届四次学术交流会论文集， 2013： 509-516.

Yu Yuanyuan， Shen Wenbo， Li Leizi. Research on the Drag Parachute Test Platform Based on the Rocket Sled Test[C]∥ Proceedings of the 6th Conference of Measurement and Control Professional Committee， Chinese Aerodynamics Research Society， 2013： 509-516.（in Chinese）

[5] 陈晨， 郭琪磊. 一种基于任意拉格朗日-欧拉方法的降落伞充气展开数值模型[J]. 科学技术与工程， 2021， 21（2）： 801-807.

Chen Chen， Guo Qilei. A Numerical Model of Parachute Deployment Inflation Process Based on Arbitrary Lagrange-Euler Method[J]. Science Technology and Engineering， 2021， 21（2）： 801-807.（in Chinese）

[6] Gao X L， Zhang Q B， Tang Q G. Fluid-Structure Interaction Ana-lysis of Parachute Finite Mass Inflation[J]. International Journal of Aerospace Engineering， 2016（1）： 1-8.

[7] Tutt B. The Application of a New Material Porosity Algorithm for Parachute Analysis[J]. International LS， 2006.

[8] Aquelet N， Tutt B. A New Fluid Structure Coupling[J]. European Journal of Computational Mechanics， 2007： 521-536.

[9] 房明. 火箭橇-阻力伞系统力学特性与流场分析[D]. 南京： 南京航空航天大学， 2018.

Fang Ming. Dynamic Characteristics and Flow Field Analyses on a Rocket Sled-Drag Parachute System[D]. Nanjing： Nanjing University of Aeronautics and Astronautics， 2018. （in Chinese）

[10] 杨品， 金楼， 胡雪鹏， 等. 基于LS-DYNA的飞机阻力伞放伞过程仿真研究[C]∥2019年中国航空科学技术大会论文集， 2019： 849-857.

Yang Pin， Jin Lou， Hu Xuepeng， et al. Simulation Study of Aircraft Drag Parachute Deployment Process Based on LS-DYNA [C]∥Proceedings of 2019 China Aviation Science and Technology Conference， 2019： 849-857.（in Chinese）

[11] 喻东明， 马啸民， 杨品. 亚声速条件下十字形伞充满时间系数的解算方法及仿真验证[J]. 南京航空航天大学学报， 2021， 53（2）： 177-181.

Yu Dongming， Ma Xiaomin， Yang Pin. Solution of Inflation Coeffcient of Cruciform Parachute under Subsonic Conditions and Simulation Verification[J]. Journal of Nanjing University of Aeronautics & Astronautics， 2021， 53（2）： 177-181.（in Chinese）

[12] Noh W F. Cel： A Time-Dependent， Two-Space-Dimensional， Coupled Eulerian-Lagrange Code[R]. Office of Scientific and Technical Information （OSTI）， 1963.

[13] 徐琳， 梁建， 段麗华， 等. 基于CSD/CFD舵面气动力流固耦合仿真分析[J]. 航空兵器， 2020， 27（2）： 47-52.

Xu Lin， Liang Jian， Duan Lihua， et al. Fluid Solid Coupling Simulation Analysis of Rudder Surface Aerodynamic Based on CSD/CFD[J]. Aero Weaponry， 2020， 27（2）： 47-52.（in Chinese）

[14] Zhang S Y， Yu L， Jia H， et al. Similarity Criteria for Canopy Porosity and Environmental Impact Analysis of Air Permeability[J]. Journal of Industrial Textiles， 2020.

[15] Skjetne E， Auriault J L. New Insights on Steady， Non-Linear Flow in Porous Media[J]. European Journal of Mechanics-B/Fluids， 1999， 18（1）： 131-145.

[16] 杨品， 胡雪鹏， 金楼. 基于LS-DYNA的某型阻力伞半折叠数值仿真研究[C]∥第八届中国航空学会青年科技论坛论文集， 2018： 774-781.

Yang Pin， Hu Xuepeng， Jin Lou. Numerical Simulation Study of a Certain Drag Parachute Based on LS-DYNA[C]∥ Proceedings of the 8th Youth Forum of Chinese Society of Aeronautics， 2018： 774-781.（in Chinese）

[17] GJB 67.4A—2008 军用飞机结构强度规范.第4部分： 地面载荷[S]. 北京： 中国人民解放军空军， 2008.

GJB 67.4 A—2008 Military Airplane Structural Strength Specification. Part 4： Ground Loads[S]. Beijing： Air Force of the PLA， 2008.（in Chinese）

Numerical Simulation on Inflation Process of Cruciform Drag Parachute

Lian Hongbo1， Yu Dongming2， Zhou Hao2， Zhao Yingtao1*

（1. School of Aerospace Engineering， Beijing Institute of Technology， Beijing 100081， China;

2. AVIC Aerospace Life-Support Industries， Ltd.， Xiangyang 441003，China）

Abstract：  In order to study the opening characteristics of the cruciform drag parachute， the ALE （Arbitrary Lagrange-Euler） fluid-solid coupling method is used to numerically simulate the inflation process of the cross drag parachute. The influence of the calculation domain size and grid densities on the calculation results is discussed， and on this basis， the shape changes of the canopy and the dynamic loading of the parachute are analyzed in detail. The results show that the simulation results of the maximum dynamic load at different opening speeds are consistent with the experiment， the error is within 5%. This can verify the feasibility and accuracy of the simulation method. When the canopy is inflated， the arms are firstly deployed， and then the airflow gathers on the top of the canopy. After being fully extended， it begins to gradually expand outwards until the canopy is completely full. The dynamic load of the drag parachute is significantly affected by the change in the shape of the canopy， and when the canopy is fully opened， the drag parachute reaches the maximum dynamic load.

Key words： drag parachute; fluid-solid coupling; numerical simulation; inflation process; flow field model; dynamic load during parachute opening