基于流固耦合方法的快速部署浮空器数值研究

龙飞， 郑亚明

(1.中国特种飞行器研究所， 荆门 448035； 2.南京航空航天大学航空学院， 南京 210016)

传统浮空器存在放飞条件限制多、升空过程风险大、部署时间长等问题，特别在一些外空间探索和军事应用场合，传统的浮空器并不适用。因此，快速部署浮空器应运而生，快速部署浮空器系统首先由载具(飞机、火箭等)运至投放区域，投放后由降落伞对浮空器系统进行减速，随后折叠态浮空器被拉出并进行快速充气，充满后与伞降系统分离。上述概念首先由Steinberg于20世纪60年代率先提出[1]，主要用于金星悬浮探测。随后，美国空军地球物理实验室(air force geophysics laboratory, AFGL)，启动模型气球系统(air-launched balloon system, ALBS)研制，由于当时数值模拟技术尚不成熟，其研究主要依赖大量地面试验和不同载具(飞机、直升机和火箭等)空投试验，但后续研究进展未见报道。而随着人类探索外空间的活动的日益频繁，快速部署浮空器再一次引起各国学者重视，并在近10年内逐渐重新启动相关研究[2-4]。

与传统浮空器相比，快速部署浮空器最复杂的工作阶段为空中快速充气展开阶段。高空或外空大气稀薄，伞降减速效果有限，必须快速充气产生浮力，而巨大冲击速度将使浮空器迅速变形，蒙皮位移又改变气源方向，进一步加剧内部流场变化不确定性，极易造成浮空器轻薄蒙皮发生破损。针对上述非常复杂的工程和科学问题，试验方法无疑是最直接的解决方案，但现有试验手段仅对蒙皮材料撕裂性能进行研究[5]，很难获得内部流场和蒙皮受力变化规律。而随着计算机硬件发展和数值方法的日益成熟，数值方案逐渐成为重要研究手段。比较有代表性的成果有：Borker等[6]提出基于势流理论的流固耦合方法求解飞艇静气动弹性问题；吴小翠等[7]采用松散耦合方式，结合量纲分析研究不同刚度构型对耦合特性影响；而张宇等[8]构建了显式动力学流固耦合分析框架，分析了某浮空器不同压差下结构响应。上述成果侧重于静气动弹性问题研究，而针对折叠态快速部署浮空器快速充气展开过程研究，相关数值研究成果几乎没有报道。

针对上述问题，以某型快速部署气球作为研究对象，采用边界附加压力模拟浮力环境，最终采用基于任意拉格朗日欧拉的多物质流固耦合方实现气球快速充气过程研究，并开展地面充气试验对设计方案和数值模型进行验证，分析充气展开过程结构场和流场变化规律，为后续快速部署浮空器的设计和优化提供一定借鉴。

1 数学模型

处于折叠状态的浮空器体存在大量不规则褶皱，因而很难建立结构内外贴体欧拉网格。因此，采用有限元描述织物及气体，拉氏网格最大优势在于可以准确跟踪材料边界，即自然满足质量守恒。此外，由于快速部署浮空器充气时间较短，且蒙皮具有不可透气性，暂不考虑不同气体间换热问题。因此，只需要考虑动量方程，其空间半离散方程为

(1)

而时间离散则采用显式计算常用的中心差分格式为

(2)

式(2)中：d为节点位移;t为时间。

由于采用逆向建模方式建立折叠态浮空器模型，四边形单元或者三角形高阶单元会造成充气前初始单元畸变，因此，采用最简单的3节点壳单元对结构域进行离散。此外，在充气展开过程中，结构域呈现小应变和大转动，因此，本文结构域采用Kirchhoff本构模型[式(3)Voigt矩阵形式]描述，通过式4转换为Cauchy应力后用于显式计算。

(3)

σ=J-1FSFT

(4)

式中:Sij为第二类Piolakirchhoff应力分量;Eij为Green应变分量;E为弹性模量;υ为泊松比;J为Jacobian行列式;F为变形梯度矩阵。而基于有限元描述的流体则不需要进行PK2应力和Cauchy应力间转换，其本构方程和状态方程为

σij=μ(vi,j+vj.i)-pδij

(5)

p=(1-γ)ρωint

(6)

式中：μ为动力黏性系数；γ为气体比热；ωint为单位质量内能。

本文中结构域和流场域都采用有限元描述，流固耦合便可以理解为两者拉氏描述材料间的接触。因此，两域中的耦合信息交互可以用接触算法实现。此外，基于拉氏描述的流体域在计算过程中，其弊端也十分明显，由于流体本构模型不同于固体，在剪切力作用下会造成较大变形，几乎在每一个时间步都会造成单元畸变，导致计算无法继续进行。因此，必须对畸变网格进行重构和流场信息的更新。而网格重构后，物质网格与参考网格拓扑关系不变，而流场域通量φ(包括质量、动量等流场信息)的求解更新[9-10]公式为

(7)

式(7)中：vi为流场域对流速度，其中对流项离散格式为

(8)

2 算例

2.1 结构与材料参数

以模型气球作为研究对象开展数值和试验研究，球体直接为1.425 m，外蒙皮由聚乙烯材料胶接而成(图1)。为防止充气前载荷重力造成球体破损，球体中间增加加强绳(1.425 m)，连接球体顶部和底部。同时为了防止高速气源在充气过程中造成蒙皮破损，在气源口周边增加内囊，内囊为高强聚乙烯复合材料，内囊上均匀布置10个直径为0.02 m的出气口。模型气球所用材料参数如表1所示。

图1 某型快速部署气球结构示意图Fig.1 Schematic figure of rapid deployment aerostat structure

表1 气球模型的材料参数Table 1 Parameter list of balloon model

为模拟空降过程的充气过程，球体顶部和底部均采用柔性约束。氦气通过充气管接入球体充气口，内囊迅速充盈后，氦气由内囊上10个出气口向球体充气(图2)。整个试验在一个大气压环境下进行，充气时间20 s，气源充气流量平稳，基本维持在0.1 kg/s，球内温度平均-30 ℃。

图2 地面充气试验Fig.2 Inflation test on ground

2.2 流固耦合模型建立

采用逆向折叠建模方法获得了折叠态模型气球模型，并建立了完整的流固耦合模型(图3)，其中结构域网格数为4.8万。此外，对球体顶部和底部施加固定约束代替加强绳作用。由于充气阀口及内囊上10个出气口面积非常小，若采用单元离散会造成局部单元尺寸过小，进而造成显式计算过程中临界时间步长Δtcrit过小和计算量急剧增加。因此，用质量点代替整个充气阀，采用点气源代替内囊上10个出气口。采用30万六面体网格对模型气球周边空间进行离散，整个流场域为高度5.5 m，直径4.6 m的圆柱体。为了模拟浮力背景，整个流场域又分为两个流场域，其中压力背景流场为最外层网格将内部流场域进行包裹，定义压力背景流场顶部单元上表面静水压力p0为一个标准大气压，自顶而下所有压力背景流场单元上表面静水压力为pi=p0+ρairghi，其中g为重力加速度，hi为第i个单元深度。

图3 模型气球充气计算网格模型Fig.3 Grid model of balloon inflation calculation

内囊材料强度远大于球体材料，前期充气试验也证明其不会发生破损，此外，研究重点在于球体充气。因此，在充气计算过程中，内囊不参与流固耦合计算，仅参与其他结构域的接触计算。内囊迅速充满后，通过内囊均布的10个出气口向球体充气，因此，在计算过程中在内囊内壁施加1 000 Pa压力，保证内囊迅速展开，点气源则根据设计在相应的单元节点上定义，而充气阀(附加质量点)和出气口(单元节点)定义矢量用于控制充气方向(图4)。

图4 内囊点气源Fig.4 Air supply position in internal capsule

3 结果与分析

图5为球体充气过程外形变化的计算结果与试验结果对比，可以发现两者较为吻合。从对比中可以发现整个充气过程大致可以分为三个阶段。

图5 外形结果变化对比Fig.5 Comparison of shape results

初期(0～2 s)：整个球体在该阶段的变化最为剧烈，氦气从内囊喷出后，造成底部迅速膨胀，随后氦气在浮力作用下，迅速流向顶部；由于气源口在气囊一侧，因此造成整个球体呈现非对称充气现象。

中期(2～6 s)：氦气在气球顶部大量聚集，顶部的褶皱开始逐渐展开，球体所受浮力克服其自重开始上浮，球体外形逐渐对称化。

后期(6～20 s)：球体顶部完全展开后，球体变化较为稳定，球体自顶向下逐渐充满，最终形成充满后外形。

图6分别为模型气球充气过程中的速度矢量图和压力云图。

图6 流场结果Fig.6 Flow field results

初期(0～2 s)：由于内囊非对称设计，球体在初始充气阶段仅造成球体局部褶皱展开。虽然出口氦气速度较大，但其密度相对于空气较小，未能形成明显的高压区，对远离气源的蒙皮影响较小。由于内囊附近蒙皮尚未形成有效约束，气源口随蒙皮变化而变化，进一步加剧流场的复杂性，流场非对称流动较为明显，球体呈现非对称展开。

中期(2～6 s)：聚集在顶部的氦气开始向未展开部分流动，顶部进一步充分展开，而球体底部受浮力影响，底部经线方向几乎完全展开，而纬线方向依然存在褶皱。不同于充气初期，该阶段的内部流场受球体约束和褶皱的影响较大，球体内部氦气会在顶部和底部形成高压区，但随着顶部褶皱的展开，顶部高压区消失，而底部高压区受底部球体约束并不消失。

后期(6～20 s)：随着球体外形逐渐对称变化，球内流场变化较为平缓，而内囊位置逐渐稳定，造成气源口附近区域形成较为稳定的高压区。

从图7所示的等效应力云图中，可以发现应力主要集中在褶皱区域，但应力集中区随着褶皱的展开而消失。同时，随着球体的充满和浮力的逐渐增加，球体顶部和底部受约束区域应力逐渐集中，特别是球体底部承受所有氦气产生的浮力，应力集中区域最为明显。

图7 蒙皮等效应力Fig.7 Von Mises stress contoursof skin

从图8所示气源所在蒙皮不同位置等效应力变化图中可以发现，气源附近蒙皮(单元S57067)受氦气冲击会造成应力急剧增加，随着球体展开，该区域应力迅速降低，但随着球体浮力增加，该处应力逐渐增加。而处于顶部的蒙皮(单元S19783和单元S17903)随着褶皱的展开，等效应力增加，但随着褶皱完全展开，其应力下降。中下部的蒙皮(单元S19007和单元S57656)虽然比顶部蒙皮率先膨胀，但该处的蒙皮并没有完全展开，而是随着球体的充满而缓慢展开，因此应力变化较为平缓。

图8 不同位置蒙皮等效应力变化Fig.8 Variation of equivalent stress of skin at different positions

4 结论

为了研究快速部署浮空器核心工作阶段的力学机理，设计并制作了用于快速部署的模型气球。并以此作为研究对象开展了地面充气试验和相对应的数值研究。其中采用任意拉格朗日欧拉方法实现了折叠态气球充气展开过程，数值结果与试验结果较为吻合，证明了本文数值方法的可行性和准确性。通过计算本文获得了试验研究很难获得的各个阶段结构和流场信息。从中可以发现整个充气展开过程大概可以分为三个阶段，初期存在非对称充气现象，导致不同位置蒙皮应力变化差别较大；随着展开过程进入中期，球体逐渐对称化；在充气后期，球体变化逐渐稳定。而整个充气过程中，位于底部的蒙皮受气源冲击和浮力影响，应力值要高于其他部位，在后续设计和制造过程中需要对底部进行结构增强。通过本文提出的数值方法可以为后续快速部署浮空器设计和优化提供参考。