赵芸可，屈秋林，刘沛清

（北京航空航天大学 航空科学与工程学院，北京100083）

近年来，中国海洋事业在民用领域蓬勃发展，水上飞机因其有别于陆基飞机，具备在水面起降的特点，可以应对应急救援和军用领域的特殊需求，故而使中国对水上飞机开发研制的需求日渐迫切。

水上飞机水面起降的能力是其最为典型的特征，性能上要求其具备水面短距起降的能力、良好的稳定性和着水性能［1］，在设计过程中必须考虑气动力和水动力的共同作用。为得到良好的起降性能，水上飞机在设计阶段便投入大量研究，主要涉及多相流、自由面捕捉和刚体六自由度运动等复杂问题的相互耦合。

水上飞机的发展可追溯到20世纪，发展初期，理论解析与实验是研究水面起降等问题的主要方法。1929年，von Karman［2］将水上飞机降落时浮筒的冲击载荷简化为二维楔形体入水冲击问题，并提出附加质量力的概念。1932年，Wagner［3］基于von Karman的工作对小抬升角的二维楔形体的垂向入水冲击问题提出进一步的修正，考虑了边界条件、水面涌起与喷溅等条件的加入。此后的诸多理论研究都是在Wagner的基础上做进一步的修正，但是对附加质量公式的修正大都针对冲击阶段，过渡阶段及浸没阶段却不适用［4］。1945年，Mayo［5］首次引入了三维流动的切片理论，结合二维流动的理论解，提出水上飞机冲击与滑水过程的理论分析并与试验结果进行对比。1948年，Milwitzky［6］总结了水上飞机冲击与滑水的切片理论，并对飞机运动过程进行了分析。1952年，Miller［7］把纵倾角固定的水上飞机冲击与滑水的切片理论拓展到纵倾角自由情况下的切片理论。NACA早期开展了大量的水上飞机设计试验，研究其气动和水动特性，包括不同抬升角［8］、不同外形参数（如机身高度、船体断阶高度等）产生的影响［9］。中国特种飞行器研究所也做了大量实验，近年来王明振［10］对水上飞机的水载荷试验方法进行了探讨，并根据试验目的对比了六自由度和三自由度试验的优缺点。魏飞和许靖锋［11］介绍了飞机模型在拖曳水池开展水上迫降试验的原理和方法。

随着计算科学技术的快速发展，近年各种数值方法在水上飞机冲击与滑水过程的研究中得到应用。与经验公式的理论方法和费用昂贵、平台搭建时间较长的模型试验相比，数值模拟方法在成本和灵活性方面具备较大优势，不仅可以提取冲击载荷等重要参数，还可以详细地展示相应的流场结构、运动状态和压力分布，为飞机设计提供参考。

当前，主要应用在该领域的数值计算方法之一是以SPH 为代表的无网格方法。2004年，EADS-CASA 的 Pentecote和 Kohlgruber［12］采 用SPH方法研究了飞机（刚体）水上坠毁多场耦合问题，分析得到由飞机底部产生吸力对结果的影响是不能忽略的。2006年，上海交通大学吴卫等［13］使用SPH方法对二维水下块体下滑进行了数值模拟，得到了水面破碎和漩涡。SPH方法在自由面大变形问题的模拟中表现良好，但在模拟气垫效应和吸力行为时表现尚欠，这对物体运动状态和冲击力预报的准确性有直接影响。

结合两相流模型的有网格技术可以较好地模拟出气垫效应和吸力行为。以有限体积法（FVM）为代表，FVM 结合VOF自由面捕捉方法已经成为一种常用的技术，其主要优势在于能够计算飞机水上运动过程中所承受的气动力和水载荷，其中尾部吸力的捕捉对飞机运动姿态的影响尤为重要。2007年，Streckwall等［14］使用基于动量法的DITCH，基于RANS方程结合VOF捕捉自由水面的Comet，模拟不同尾部形状机身的着水冲击过程，并与实验结果进行了对比验证。2009年，北京航空航天大学陆士嘉实验室刘沛清、屈秋林等［15］使用FVM结合VOF方法，用Fluent软件模拟了ARJ21飞机在平静水面上的水上迫降过程，较好地捕捉到水面的变形并提出了最佳迫降姿态。2013年，上海交通大学邱良俊和宋文滨［16］将总阻力拆分为水动阻力和气动阻力两部分，分别用Fluent和CFX逐步分解计算，再将2种载荷的模拟结果合并，同时代入当前给定的发动机拉力，通过迭代模拟水上飞机的起飞过程。2015年，北京航空航天大学刘沛清、屈秋林［17］等用整体运动网格（Global Motion Mesh，GMM）方法模拟了二维圆柱垂向入水和飞机水上迫降运动，迫降过程同时伴随冲击入水和水面滑行运动，模拟得到的压力分布状况和运动状态变化历程与实验结果进行比对，结果对应良好，表明该方法解决此类问题具有较好的适用性。2018年，上海交通大学张浪等［18］采用VOF方法，用Fluent软件进行模拟，并在每一步代入当前的发动机拉力，模拟了姿态角固定情况下飞机的起飞过程。西北工业大学段旭鹏［19］等在开源平台OpenFOAM 两相流动态解算器中添加了激励盘来模拟螺旋桨的拉力作用，并求解了水上飞机以固定姿态角进行滑水运动的过程。

用数值方法模拟水上飞机水面降落的全过程存在以下几个难点：①从触水冲击、水面滑行到滑行停止，该过程运动距离长，对于传统的动网格方法而言计算域过大，而有网格方法本身又对网格精细程度有较高的要求，导致计算成本难以承受；②飞机在气动载荷和水载荷共同作用下做六自由度运动，模拟时对计算的收敛性有较高要求，且要求其结果可靠。目前，很多研究将该过程设定为固定姿态角或者受迫运动来进行模型简化，或者仅截取冲击或滑水某一阶段进行模拟。

本文采用了流体求解器、整体运动网格方法、VOF自由面捕捉和六自由度模型相结合的方法，模拟了水上飞机从水面上方冲击入水到滑行基本停止进入漂浮姿态的自由运动全过程。该方法的优点在于：①避免了传统动网格的变形重构过程，即避免了因运动姿态剧烈变化导致的网格畸变所造成的计算发散，同时大大降低了计算域即网格数量，降低计算成本的同时增强了收敛性；②传统的变形重构网格需要采用非结构网格，而整体运动网格理论上可以使用全结构网格，结构网格在使用VOF方法时模拟效果有明显优势，在模拟精度要求不变的基础上结构网格所需的网格数量可大幅减少。基于以上优点，本文选择整体运动网格方法模拟某型水上飞机水面降落的全过程。

1 整体运动网格方法及其验证

1.1 整体运动网格

整体运动网格的开发基于Fluent商业软件中的动网格方法，其特征是整个计算域跟随飞机做刚体运动。计算时采用地面坐标系，每一步的迭代网格节点坐标都会随着飞机的运动更新，故而不需要配合滑行距离而扩大计算域，即省去了传统变形网格的变形区域。其中水体并不跟随做刚体运动，而是始终维持地面坐标系下规定的水面高度，通过UDF使用VOF边界条件设置即可完成该设定。

VOF 模 型 由 Hirt、Nichols［20］和 Muzaferija等［21］提出，该方法可以捕捉多相流之间的自由交界面，其中各相流体之间不相容同一个网格中，各相的体积分数之和为1。假设第q相流体有一个体积分数αq，在某个网格中，当αq＝0时，则说明该网格中不存在第q相流体；当αq＝1时，则说明该网格中只有第q相流体；当0＜αq＜1时，则说明该网格中除了第q相流体外，还有其他相的流体存在。第q相流体的体积分数方程如下：

式中：ρq为流体密度；vq为流体速度；Sαq为质量源项；˙m为相间转化质量流率。计算时，初始设置水面高度为0，水面以上水相体积分数为0，水面以下为1，水气交界面为0.5。

六自由度模型通过求解平动和转动方程来得到质心的位置和物体的姿态角。平动方程在惯性坐标系（地面坐标系）下求解。

式中：m为飞机质量；vg为重心的平动速度；fg为飞机表面力与重力之和；下标g表示地面坐标系。

转动方程则在机体坐标系下求解：

式中：K为惯性张量；ω为飞机绕重心角速度；Mb为飞机所受绕重心力矩；下标b为机体坐标系。

1.2 计算格式

本文采用Fluent软件作为计算平台求解非定常RANS方程，采用增强壁面处理可实现k-ε湍流模型，壁面处采用增强壁面处理。压力速度耦合采用SIMPLEC算法，压力项采用体积力加权格式离散，动量方程中的对流项由三阶MUSCL格式离散。湍流方程中的对流项由二阶迎风格式离散，扩散项由二阶中心差分格式离散，非定常项由二阶隐式格式离散。本文中所有工况均采用此计算格式。

1.3 飞机水上迫降计算与实验结果的对比验证

为了验证整体运动网格在模拟飞机入水冲击过程中的精度和合理性，屈秋林等［17］针对NACATN-2929的模型F进行了水上迫降的数值预报，并与实验结果进行了详细的对比分析。实验NACA-TN-2929［22］研究了不同外形的机身后段对迫降的影响，其中模型F机身轴线尾部上翘为宽机身，如图1所示。

图1 实验NACA-TN-2929中的模型F机身模型［22］Fig.1 Fuselagemodel of Model F in NACA-TN-2929 experiment［22］

图2 模型F的表面网格［17］Fig.2 Surface grid of Model F［17］

图3 模型F的计算结果与实验结果的对比［17］Fig.3 Comparison of calculation results of Model F with experimental results［17］

模型F的飞机表面网格如图2所示［17］，采用半模计算，网格总量270万。飞机位于半立方体的计算域中央，模型对称面和计算域的对称面重合，前端机身轴线和计算域上下边平行，计算域上游边界距离机头2倍机身长度，计算域侧边界距离翼梢1倍机身长度。对称面采用对称边界条件，外场边界条件为速度入口。图3［17］显示了该工况的实验结果与COMET数值模拟方法、DITCH动量方法和整体运动网格方法模拟结果的对比，总体上数值计算的模拟结果与实验结果吻合度较高。其中，水平速度曲线显示实验结果的初始加速度很大，不同于模拟结果的加速度由小增大，而事实上模拟的结果更为接近真实情况。2007年，Streckwall等［14］进行模拟时也遇见了相同的问题，他们认为由于实验年代较早，受限于测量手段实验结果难免存在误差。造成误差的原因可能是水平速度的初始测量点选取在运动后的一段时间。

根据上述结果分析可见，整体运动网格技术在数值预报飞机水上迫降动态特性是合理可靠的。

2 某型水上飞机水面降落过程模拟

本节使用整体运动网格数值预报了某型水上飞机的水面降落过程。

2.1 模型信息和初始条件设置

如图4所示，某型水上飞机的机翼参考面积约156m2，襟翼为简单开缝式。起飞内外襟翼偏角20°，着水时内外襟翼偏角45°。机翼展弦比9，根稍比2.4，1／4弦线后掠角7.13°，船体前体长度14 m，断阶宽度2.7 m。本节分析了质量为49.8 t，转动惯量为1.94×106kg·m2，以飞机头部尖端为原点，重心位于前重心（15.91，0，1.15）m的某型水上飞机。设置初始运动状态水平速度44.72m／s，垂向入水速度1.5 m／s，初始姿态角5°，重心距水面高度4.25m，断阶最低点距离水面高度1.15m。

计算模型为半模，飞机位于半立方体的计算域中央，模型对称面和计算域的对称面重合，前端机身轴线和计算域上下边平行，计算域上游边界距离机头2倍机身长度，计算域侧边界距离翼梢1倍机身长度。对称面采用对称边界条件，外场边界条件为速度入口。飞机的中央翼盒、机翼和襟翼采用非结构网格，其余部分采用结构网格，可以确保水面基本处于结构网格中。网格数量约600万，计算迭代步长0.01 s，共30 s。

图4 某型水上飞机几何特征和计算网格Fig.4 Geometric features and computational grid of a certain seaplane

2.2 数值模拟结果

如图5所示，模拟结果展示了地面坐标系下的飞机姿态角、水平速度、重心距水面高度、垂向速度、过载和角加速度等变化历程。其中过载是无量纲参数，垂向过载不包含重力作用的部分，定义如下：

式中：Fx为飞机所受水平过载；fx为飞机所受水平合力；Fy为飞机所受垂向过载；fy为飞机所受垂向合力；g为重力加速度。本文主要依据飞机运动状态和过载的变化，将飞机整个降落过程分解为3个阶段，即冲击阶段—滑水阶段—漂浮阶段。

1）冲击阶段。飞机具备一定的初始动能入水，冲击造成的水载荷导致整体所受过载剧烈振荡，且随着动能的耗散振荡幅度逐渐减小；该阶段运动状态参数变化剧烈，垂向速度呈现大幅振荡，重心升沉幅度较大，期间可能发生弹跳（本文中即出现弹跳现象），并伴随姿态角剧烈俯仰。

2）滑水阶段。运动状态参数振荡幅度减小，飞机重心距水面高度变化先趋于稳定，后缓步下沉，此间不在发生弹跳或剧烈的俯仰。过载振荡幅度较小，滑水运动产生的水动力主导运动姿态变化。

3）漂浮阶段。所有参数随时间发展趋于稳定，飞机的水面姿态接近静水漂浮状态，静水浮力成为主要的外部载荷，垂向过载趋近于1（即外力作用主要为静水浮力，其值接近飞机自身重力）。

2.2.1 冲击阶段

冲击阶段是飞机在水面降落的初始阶段，该阶段包含了飞机从接近水面到触水冲击的过程，主要作用在该阶段的冲击载荷是垂向的入水冲击运动与水平方向的滑水运动共同产生的。由冲击引起的大过载会导致飞机运动姿态产生剧烈变化，其中飞机触水弹起现象是冲击阶段的典型特征之一。该阶段从开始时刻持续到约4.8 s，定义重心距水面高度位置变化“峰-谷-峰”为一次冲击过程，本工况冲击阶段共发生3次冲击，命名为首次冲击、第一次二次冲击、第二次二次冲击。图6展示了冲击阶段飞机的运动状态和过载的变化历程。

1）首次冲击

计算初期飞机自水面上方向水面迫近，触水前作用于飞机的外载荷为气动力。如图6所示，从初始时刻到0.4 s左右，垂向过载约为0.8不足以抵抗重力，故指向下的垂向速度不断增大，由初始的1.50m／s到最大2.53m／s；0.4～0.55 s尚未触水，而向上的垂向过载开始增加。本文认为发生该情况的原因是飞机逐渐迫近水面时，机身底部与水面之间形成了一层空气垫，气垫效应使飞机受到了向上的压力，故而垂向过载增加。0.55 s时姿态角为3.63°左右，而此前飞机已经展现出明显的低头趋势，可见初始的低头趋势主要是气动力造成的。

如图7所示，为显示清楚，压力云图的范围为正负0.3倍大气压（图中是相对标准大气压值，后同）。0.55 s时断阶最低点触水，由入水冲击产生的向上的垂向过载陡升，其峰值出现在飞机刚发生触水的片刻后0.61 s，为1.68，这与通常水上迫降的载荷系数量级是一致的。随着飞机的进一步下沉，下沉速度减小，垂向过载呈缓慢减弱趋势，0.9 s飞机下沉速度降为0m／s时，飞机重心距水面高度位置到达最低点3.71m，此时断阶距离水面高度0.39m。此刻之后在垂向载荷的作用下机身开始第一次反弹，随着重心被不断抬高，机身入水体积减少，静水浮力降低，垂向冲击过载陡降，但依旧可以维持飞机向上弹起。2.07 s时重心处于第一次弹起的最高值4.78m，而后开始回落，回落过程中冲击过载逐渐上升，开始第一次二次冲击。

2）第一次二次冲击

图7 冲击阶段0.55 s水面状况和机身底部压力云图Fig.7 Water surface condition and pressure contour of the bottom of fuselage at 0.55 s during impact stage

第一次二次冲击和首次冲击过程相似，但由于动能的耗散，过载和运动姿态的变化幅度都相较首次略有削弱。第一次二次冲击从2.07 s开始，期间垂向过载峰值发生在2.82s，为1.32。当3.1s时重心距水面高度回落到最低位置4.04m，此时断阶距离水面高度0.07m。此后重心又开始呈现出上升趋势，开始第二次弹起过程，第二次弹起过程与第一次类似，4.02 s到达第二次弹起的重心位置最高值4.29m，继而回落开始第二次二次冲击。

3）第二次二次冲击

第二次二次冲击的力度远小于首次冲击和第一次二次冲击，从4.02 s开始到4.8 s结束，期间最大载荷系数在1.11。考虑到重心距水面高度回升后再未弹起到断阶离开水面的高度，故该阶段并未包含“峰-谷-峰”的完整过程。此次冲击之后冲击阶段结束，进入滑水阶段。

在整个冲击阶段，角加速度和冲击过载变化趋势相似，故认为垂向载荷的作用对角加速度产生主要影响。图8展现了从机身触水到第一次弹起过程中，飞机机身底部压力云图。其中，前三张图在机身断阶范围处产生了明显的负压，即吸力。一般认为，断阶处发生的机身底部曲线大曲率变化是负压产生的主要原因，当流动绕过这里时因为加速而产生相对负压，这与水上迫降中得到的经验相同［17］。结合图6第三幅图和第四幅图分析可得，0.6 s冲击发生时正压主要作用位置在重心位置之前，断阶处则产生负压，水载荷作用在飞机机身底部产生了很大的抬头力矩，伴随飞机姿态角显著增大机尾逐渐下沉靠近水面；0.9 s时冲击产生的正压力减小，同时因为姿态角增大，机身下方正压作用位置后移，使抬头力矩减小；0.9 s后机身尾部开始入水，由于机尾距离重心较远而且机尾在重心后方，尽管机尾处的冲击过载较弱，仍然产生了较大的低头力矩，直到1.2 s左右低头力矩达到最大值，此后由于第一次反弹开始重心距水面高度上升，作用于机身尾部的正压力减小，低头力矩开始减小。从机身触水开始到1.61 s期间飞机姿态角始终维持增大的趋势，1.61 s时达到最大值10.72°；2.0 s时，飞机已经弹起离开水面，作用在机身底部的主要为气动力。

图8 冲击阶段0.6～2.0 s机身底部压力云图Fig.8 Pressure contour of the bottom of fuselage for 0.6-2.0 s during impact stage

2.2.2 滑水阶段

在冲击阶段的弹跳运动结束后，飞机垂向运动因为能量的耗散而不再对运动状态构成主要影响，此阶段飞机的水平速度依旧较大，故而滑水运动的影响开始占据主要地位。滑水阶段的姿态角在前半程基本维持不变，后半程姿态角开始减小的该时刻为分界点，滑水阶段被分为两个阶段，分别为机身后段未入水和机身后段入水，如图9所示，这同时也解释了姿态角变化的原因。

图9 滑水阶段运动状态和过载变化历程曲线Fig.9 Motion state and overload history curves during water skiing stage

1）机身后段未入水

在冲击阶段结束后，飞机进入滑水阶段。由滑水产生的水平阻力降低水平速度，从而气动力减小，飞机重心距水面高度下降，机身入水深度增加，水载荷开始承担主要角色。滑水阶段前半程主要滑行接触面为包括断阶部位的机身前段，此阶段姿态角基本维持在8.4°，一直持续到11.5 s。此时重心距水面高度位置缓慢下降，从4.09m下降到3.58m左右。

从图10所示的机身与水面位置关系可以看出，滑水阶段的前半段，机身尾部处于水面上方，飞机水平速度缓步降低但依旧保持向前滑行的运动，机身前方水面涌起并绕过机身，使上机身底部产生了一个冲击载荷区，该冲击载荷区的位置依旧在重心之前，伴随着重心位置的下沉逐渐向机头部位移动，但同时冲击载荷降低，使该阶段水载荷产生的抬头力矩与气动力产生的力矩恰巧维持了基本平衡，从而保持了姿态角的稳定。

图10 滑水阶段水面状况和机身底部压力云图（机身后段未入水）Fig.10 Water surface condition and pressure contour of the bottom of fuselage during water skiing stage（tail does not touch water）

2）机身后段入水

如图11所示，重心垂向位移呈现缓步降低的趋势，11.5 s时刻机身尾部开始入水，产生了较大的低头力矩，姿态角开始降低。同时机身后端入水使入水体积增加，导致12～13 s水平阻力增大，但逐步降低的水平速度又进一步限制了水平阻力的增加。

滑水阶段的最后，重心距水面高度位置下沉至2.66m，垂向过载趋近于1，主要用于抵抗重力。

2.2.3 漂浮阶段

当水平速度降低到一定程度，水平方向滑水运动产生的载荷逐渐失去主要位置，静水浮力开始成为水载荷的主要构成部分。漂浮阶段从17 s开始持续到最后。

图11 滑水阶段水面状况和机身底部压力云图（机身后段入水）Fig.11 Water surface condition and pressure contour of the bottom of fuselage during water skiing stage（tail touches water）

图12 漂浮阶段运动状态和过载变化历程曲线Fig.12 Motion state and overload history curves during floating stage

如图12所示，此阶段的垂向过载趋于1，说明此时垂向除重力外，受力基本为静水浮力；飞机保持惯性继续向前滑行，并在水平阻力的影响下稳步减缓滑行速度，当飞机水平速度进一步变小，作用在重心前端的冲击压力减弱，其提供的抬头力矩减弱，姿态角减小；随着飞机姿态角的减小，机身前部逐渐入水，冲击压力区逐渐前移，但总体上冲击压力已经很弱，抬头力矩不再随之增加。整个过程中姿态角由4.32°降低至2.71°，最终逐渐接近飞机静水漂浮时的姿态。图13展示了此阶段飞机的水面状况。

图13 漂浮阶段水面状况和机身底部压力云图Fig.13 Water surface condition and pressure contour of the bottom of fuselage during floating stage

3 结 论

本文利用整体运动网格技术结合VOF方法和六自由度模型，模拟了某型水上飞机水面降落的全过程，得到了较好的模拟结果，由此得到整体运动网格方法在解决水上飞机水面降落问题时具备良好的适应性。该方法计算成本低，稳定性好，模拟结果可靠，在处理长距离滑水问题时优点尤其突出。

本文的模拟结果展示了运动状态参数和受力状况的变化历程，并结合压力分布云图和水面状况，对水面降落过程进行定性与定量分析。本文将某型水上飞机在水面上降落过程划分为3个阶段，即冲击阶段、滑水阶段、漂浮阶段。各个阶段的主要特征如下：

1）冲击阶段。飞机垂向的冲击运动产生的过载主导该过程。该过程运动姿态变化剧烈，本文工况中出现的弹跳运动为该阶段的典型状况。

2）滑水阶段。垂向方向的动能随着冲击很快被水的阻尼耗散，而此时仍然具有较大的水平速度，因此水平方向滑水运动产生的过载主导该过程。该过程中飞机重心距水面高度稳步下沉，运动姿态不再剧烈变化，而是在各载荷的相互作用下趋向稳定。在本文工况中，机身后段入水前滑水阶段产生的稳定（相对于冲击阶段）冲击压力使得姿态角大小维持稳定，直到机身后段入水。机身后段入水后姿态角缓步减小。

3）漂浮阶段。随着水平速度的逐渐减小，滑水运动的作用减弱，静水浮力接替主导位置。该阶段冲击压力很弱，从而抬头力矩减弱，姿态角逐渐下降，重心距水面高度几乎不变，整体运动状态趋于最终的静水漂浮姿态。