基于显式动力学的软式飞艇流固耦合计算框架

张 宇, 王晓亮

(上海交通大学 航空航天学院，上海 200240)

飞艇能完成诸如通信、侦察、运输等任务，在军事和民事应用中具有巨大的前景.2009年3月，上海交通大学“致远一号”飞艇验证艇在苏州成功试飞，验证了导航控制、数传、测控等技术方案的可行性，随后于同年11月完成飞行试验[1].2014年，美国国家航空航天局(NASA)开展了“平流层飞艇设计”竞赛，目标是以平流层飞艇代替卫星[2].2015年，北京航空航天大学在内蒙古成功开展了平流层飞艇的长时留空试验[3].飞艇从结构上分为硬式飞艇、软式飞艇和半硬式飞艇[4].其中，软式飞艇不含内部骨架，艇身由柔性蒙皮材料组成，通过在气囊内注入氦气形成内外压差从而维持外形.与同体积的硬式飞艇相比，软式飞艇基于内外压差成型，易于制造和保养，且由于不含内部骨架、自重较轻、能携带更多的任务载荷，所以小型飞艇多为软式飞艇.但由于软式飞艇的结构特点，内外压差不易维持容易导致其流固耦合现象明显.

对软式飞艇而言，其流固耦合现象主要来源于蒙皮薄膜结构，针对薄膜结构的流固耦合分析方法已经获得了较好地发展.Lee等[5]通过求解三维Navier-Stokes(N-S)方程获得气动力，结构分析采用共旋有限元和共旋壳单元方法，分析了流固耦合效应对微型机器人薄膜扑翼推进性能的影响.De Nayer等[6]对一充气柔性膜结构在不同来流雷诺数下的瞬时响应进行了数值模拟，流场由大涡模拟获得.在飞艇的应用方面，Bessert等[7]通过流体求解器VSAERO和固体力学求解软件Abaqus，研究CL160飞艇在不同俯仰角下的升力系数变化.王晓亮等[8]通过Fluent和Abaqus软件，形成了非线性动态流体-结构交错积分耦合法，使用隐式方法对某平流层飞艇在突风环境下的结构响应进行了分析.Liu等[9]通过SIMPLE算法和非线性有限元方法对比了弹性与刚性飞艇模型之间的压力分布与流场参数差异.吴小翠等[10]研究了在不同刚度构型下，飞艇的定常流固耦合特性，并利用量纲分析证明了模型与实物间的相似律.由此可见，目前已经初步形成了关于飞艇的流固耦合分析方法，但上述方法主要针对飞艇的单向流固耦合特性进行研究，且由于隐式方法在求解固体结构响应时需要迭代计算，对于高度非线性问题可能难以收敛.

本文以Fluent与Abaqus作为求解器，使用双向流固耦合方法，通过求解非定常雷诺平均NS(URANS)方程获取艇身表面气动力，由显式动力学方法获得飞艇的结构响应，由径向基函数(RBF)方法进行气动力和位移数据交换，由Delaunay映射方法完成流场网格更新，由Fortran子程序完成气动载荷加载.最后，应用上述框架分析某软式飞艇在不同压差下的结构响应规律.

1 几何及环境描述

一般的硬式飞艇结构如图1所示.由于研究对象为软式飞艇，所以在硬式飞艇模型的基础上去除了吊舱和内部骨架，艇身外形由参数化截面(PARSEC)方法获得[11].为体现艇身的旋成体特点，在原有的11个参数上只需保留8个即可确定飞艇外形.PARSEC方法示意图如图2所示.其中：rh为头部半径；xd为最大横截面位置；rd为最大横截面半径；kd为最大半径处曲率；αt为尾部偏移角；βt为尾部张开角；dt为尾部厚度；ht为尾部高度.艇身外形所对应的各参数取值如表1所示.飞艇尾翼翼型为NACA0012.y坐标由如下的6次多项式方程确定[12]：

表1 艇身外形参数Tab.1 Shape parameters of hull

图1 硬式飞艇外形示意图Fig.1 Schematic diagram of rigid airship

图2 PARSEC方法Fig.2 PARSEC method

(1)

式中的未知量ai可由如下线性方程组获得：

(2)

式中：cte为无量纲翼型弦长，此处取为1.

软式飞艇所处的环境示意图如图3所示.其中：p为艇身内部压力；p∞为外界远场压力；L为艇身总长；D为最大直径；denv为蒙皮厚度；Eenv、ρenv和νenv分别为蒙皮弹性模量、密度和泊松比；g为重力加速度，方向垂直向下；“×”表示垂直纸面向里的横风.飞艇内外压差为Δp=p-p∞，飞行高度h=20 km，空气密度与动力黏度分别为ρair=0.089 kg/m3和μair=1.422 Pa·s.实际艇身通过PARSEC模型放大45倍而来，尾翼前缘距离原点0.75L，尾翼后掠角为30°，翼尖长度为0.08L，两十字尾翼翼尖距离为1.5D.艇身内部压力和外界远场压力的大小随高度的变化而变化，来流速度为v∞.上述有关变量的取值如表2所示.

图3 飞艇几何及其环境简图Fig.3 Diagram of airship geometry and its environment

表2 飞艇几何及蒙皮材料参数Tab.2 Parameters of airship geometry and envelope material

2 数值方法

2.1 计算流体力学(CFD)数值方法

使用连续性方程和N-S方程即可求解三维不可压缩流动问题，微分形式的连续性方程为

(3)

(4)

湍流模型采用Spalart-Allmaras模型，该模型适合求解航空外流场问题，使用壁面增强函数以便获得较好的近壁面流场信息.求解器采取不可压缩形式的压力基求解器，压力-速度耦合格式为“Coupled”，空间梯度离散方法为“Least Squares Cell Based”，压力项采用二阶格式离散，动量和修正湍流黏度采用二阶迎风格式离散.CFD时间增量参数为dtf，每一增量步内迭代步数设置为50以保证计算收敛.

流体计算域如图4(a)所示，速度入口距离飞艇中心20L=900 m，压力出口距离飞艇中心40L=1 800 m，远场边界距离飞艇中心20L=900 m.网格划分工具为Pointwise，在飞艇表面布置三角形网格，采取T-Rex方法生成边界层网格，远场流域采用非结构四面体网格填充.

图4 流体计算域及其网格Fig.4 Fluid computational domain and its meshes

2.2 可靠性验证

为降低流场网格因素带来的误差，有必要进行网格无关性检验.通过对艇身表面网格疏密程度的调整，达到既能保证结果精度又能节约计算资源的目的.为此，划分了4种疏密程度的网格，分别统计在v∞=12 m/s横风作用下飞艇z方向的阻力系数CDz，相关飞艇流场的网格信息如表3所示.由表3可知，中等疏密程度的网格已经能较好地满足计算结果的精度要求.为节约计算资源，选取中等网格划分策略进行后续分析，最终的流场网格如图4(b)所示.

表3 飞艇流场网格信息Tab.3 Mesh information of airship fluid

在上述数值方法及网格单元设置下，飞艇受到v∞=12 m/s横风作用时艇身表面与不同x剖面交线上的无量纲网格高度y+分布如图5所示.其中：

图5 不同交线上的y+分布Fig.5 Distribution of y+ on different intersection lines

径向坐标为艇身表面的y+值；圆周坐标θ为交线绕x轴的夹角.由图5可知，在不同艇身交线上的y+最大值均约为3，大部分处于0.5～2范围内.故所选取的数值方法及网格划分策略满足增强壁面函数及湍流模型对边界层质量的要求，能较好地捕捉艇身外流流场信息.

2.3 动网格更新方法

流场网格更新通过Delaunay映射方法实现，该方法针对中等变形问题具有很高的稳健性，避免了求解大型矩阵的高耗时缺点.关于该算法的具体原理参考文献[13].借助Fluent中的用户自定义函数(UDF)功能，将Delaunay方法集成在DEFINE_GRID_MOTION宏中，即可实现网格更新效果.

2.4 计算结构动力学(CSD)数值方法

对于结构动态响应，一般可通过隐式或显式动力学算法获得.隐式算法需在每个时间步内迭代求解平衡方程，当前后两次迭代值满足误差要求时才会进入下一时间步，而显式算法在任一时间步内无需迭代求解大型矩阵.相较于隐式算法，在结构设置合理的前提下，显式算法不存在收敛问题.因此，隐式算法一般只用于求解常规的动力学问题，而显式算法更适合高度非线性问题且能更好地捕捉动态响应细节.

显式动力学关于运动方程的中心差分形式可表示为

(5)

(6)

(7)

固体网格采用三角形壳单元划分，经过网格无关性分析后确定网格单元总数为 3 180.采取基于中心差分的显式动力学对运动方程进行积分求解.但中心差分算法是条件稳定的，其时间增量与应力波传播速度vd与网格特征尺寸Lc密切相关，在此设定CSD的时间增量参数为dts，其值一般要满足dts≤Lc/vd.计算时将十字尾翼固定，在飞艇内部施加均匀分布的法向压力作为内压，由CFD计算得到的气动力通过Fortran语言二次开发子程序施加到软式飞艇蒙皮外表面.

3 数据传递方法及验证

采取RBF实现气动力与位移的传递.RBF插值方法不依赖于网格拓扑结构，只须获取CFD和CSD的插值坐标即可进行插值.需要说明的是，进行数据传递时CFD与CSD的插值点分别为物面网格节点和壳单元中心点.

3.1 位移传递

假设已知S个结构网格单元中心的坐标ϑs和某一方向上的位移u，则可建立关于位移的插值格式为

u=Φω

(8)

式中：Φ为径向基函数矩阵；ω为系数向量.式(8)的展开形式为

(9)

式中：ξji为点ϑsj与点ϑsi的距离与径向基函数作用半径ref的比值，即

(10)

径向基函数选取 Wendland二阶光顺格式[14]，可表示为

(11)

由于矩阵Φ是对称正定的，其逆可通过Cholesky分解快速获得，所以可以获得系数向量ω的分量值，于是CFD网格节点的位移可由下式获得：

(12)

k=1,2,…,Sf

式中：m为结构壳单元编号；k为流体网格节点编号；Sf为流体物面网格节点个数.

3.2 气动力传递

气动力传递过程与上节类似，在此不再赘述.如果所面临的流体插值点数目远大于固体插值点数目的情况，可使用局部径向基函数插值方法来提高效率，具体过程参考文献[15].

3.3 流固耦合过程

关于软式飞艇的非定常流固耦合流程如图6所示.一个子循环步的耦合过程如下：① 给软式飞艇内部施加均匀分布的法向压力，并在外部施加幅值为0的气动力，此步骤用于形成因内部压力而变形的飞艇初始外形；② CSD模块在当前子循环步中以增量dts向前推进Δt，判断CSD模块向前推进的时间是否达到dtf，即是否满足|Δt-dtf|≤ε，若满足则进行下一步，否则继续进行CSD模块计算，其中ε为时间允许误差，在数值上等于dts；③ 将CSD模块的单元中心变形量通过RBF方法传递给CFD模块网格节点，并由Delaunay映射方法更新网格节点位置，开启CFD模块计算进程；④ 当CFD模块完成当前时间步的计算后，将气动力传递到CSD网格单元.重复上述过程即可完成双向流固耦合的交替计算，且从上述过程可知实际的耦合步长为dtf.

图6 流固耦合计算流程Fig.6 Calculation process of fluid-structure interaction

3.4 流固耦合方法验证

为验证所提方法的可靠性，以如图7所示的基于网格的并行化代码耦合接口(MpCCI)案例为研究对象，将密度为 1 000 kg/m3、弹性模量为0.1 GPa、泊松比为0.49的弹性板置于长方体导管中.入口速度为8 m/s，出口为压力出口，表压为0，管壁为无滑移壁面，弹性板的上端固定在管壁上，取弹性板下端顶点为控制点，其余参数及设置参考文献[16].分别利用MpCCI和本文所提方法(FSI-ED)耦合Fluent和Abaqus进行计算，最终获得控制点的振动幅值随时间的变化情况如图8所示.其中：ux为控制点在来流方向上的位移.由图8可知，MpCCI与FSI-ED计算所得的振动频率分别为44.8 Hz与45.0 Hz，其相对误差为0.45%，所得最大振幅的平均相对误差为6.1%.

图7 平板冲击案例计算域(m)Fig.7 Computational domain of plate impact case (m)

图8 FSI-ED与MpCCI计算结果对比Fig.8 Comparison of calculation results of FSI-ED and MpCCI

选取NACA0014机翼进行气动力传递的精度验证，其机翼表面网格如图9所示.来流马赫数为 0.839 5，攻角为3.06°，当地大气压为 100 311.75 Pa，当地温度为300.9 K.表4统计了通过积分获得的作用在流体和固体网格表面单元上3个方向的合力和力矩大小，其中：Fx、Fy和Fz为机翼在x、y和z方向受到的合力；Mx、My和Mz为机翼在x、y和z方向受到的合力矩.由表4可知，在气动力的传递过程中几乎没有能量损失，所有对比项的相对误差均在0.1%以下.

表4 传递前后的力与力矩对比Tab.4 Comparison of forces and moments before and after transmission

图9 NACA0014机翼表面网格Fig.9 Surface mesh of NACA0014 wing

综上，所提数值方法能较精确地对振动过程进行描述且能以较高的精度保证气动力的传递，可满足流固耦合问题的数值计算.

4 软式飞艇的动态结构响应

4.1 流固耦合时间步长

对于流固耦合分析而言，时间步长的选取至关重要，不当的时间步长甚至可能导致错误的结果.CSD模块的时间增量步dts与应力波传递速度vd和网格特征尺寸Lc相关，这两个参数直接反映了飞艇模型的固有频率.Abaqus在模型数据检查阶段会自动计算上述两个参数，并对显式时间步长进行估计.根据其估计量选取dts=5×10-5s.

对CFD模型而言，一般采取隐式算法，因此理论上其时间步长不受限制，但过小的时间步长会导致计算资源有所增加，过大的时间步长可能导致失去真实的结构响应细节.对软式飞艇而言，其振动频率大多低于5 Hz[17].出于对计算资源及计算精度的综合考量，选取dtf=1×10-2s.

4.2 结构振动响应

设定软式飞艇的内外压差Δp的范围为[208,852] Pa，增量为92 Pa.应用所提计算框架对不同压差下整体艇身的结构响应进行分析.

对能维持一定整体刚度的软式飞艇而言，在尾翼固定的情况下，艇身头部的振动幅值最大，能够较好地反映软式飞艇在横风作用下的振动特性[8]，故取飞艇头部顶点作为控制点进行监测.在不同压差情况下，飞艇控制点在时域上的响应曲线如图10所示，其中uz为控制点在z方向上的位移.从时域曲线可见，当Δp低于576 Pa时，艇身振动幅值包线随时间单调递减；当Δp高于668 Pa后，振动幅值包线不再单调变化，而是呈现出较为复杂的动态响应.通过快速Fourier变换对时域曲线进行转换得到如图11所示的频域响应曲线，其中：f为频率；A为幅值.Δp与fm的关系曲线如图12所示，通过对不同压差下的振动主频fm进行统计，可见振动主频与飞艇的内外压差之间呈近似线性关系：

fm=5.8×10-5Δp+3.667 9

(13)

图10 控制点的时域响应曲线Fig.10 Time domain response curves of control points

图11 控制点的频域响应曲线Fig.11 Frequency domain response curves of control points

图12 Δp与fm的关系曲线Fig.12 Δp versus fm

由式(13)可见，斜率仅为5.8×10-5，因此可认为在文中设定的压差范围内，振动主频fm几乎不受压差Δp的影响，即当软式飞艇受到的压差能基本维持艇身刚度之后，继续增加内压对振动频率的影响较小，这与相关试验结论吻合[17].

5 结语

本文基于URANS和显式动力学方程，通过Delaunay映射方法完成流场网格更新，采用径向基函数完成气动力和结构位移在流体单元和固体单元之间的传递，进而构建软式飞艇在非定常气动力激励下的显示动力学结构响应计算框架.该框架适合双向流固耦合计算，能较好地反映软式飞艇的结构响应规律，体现显式计算格式效率高的优点.飞艇艇身由PARSEC方法获得，应用上述框架对该软式飞艇的计算结果表明，当软式飞艇能基本维持艇身刚度后，继续增加内压对振动频率的影响较小，计算结果与相关试验结论吻合.本文形成的计算框架对于高空气球以及其他膜结构高空飞行器在瞬态气动力作用下的结构响应分析均适用.