， *， ， ， ，

(1.武汉理工大学 工程结构与力学系，武汉 430070；2.高速水动力航空科技重点实验室，荆门 448000)

1 引 言

水陆两栖飞机上半部按照飞机结构设计，下半部是双曲面的船型结构，能在陆地和水面起降，是进行海上巡逻、救援和大型灭火等便捷有效的工具。水陆两栖飞机水上起降较陆上起降更为频繁，尤其在海面复杂气象条件下降落时，复杂的水动载荷使得飞机降落时的滑行阻力增大，过载增高，纵摇运动及升沉运动变化剧烈，甚至会引起机身结构的变形破坏，对机体及机乘人员造成伤害。在规定的气象水文条件下起飞和降落的能力，即耐波性，是决定水陆两栖飞机性能的重要因素。

20世纪中期NASA进行了大量的水上飞机耐波性实验，研究了不同型号的水上飞机在水上起降时的耐波性能[1,2]。国内于20世纪70年代也开始了水上飞机的理论与模型物理水池试验的研究，对水上飞机耐波性分析方法等做了一定研究[3]。但是，由于飞机与波浪耦合作用的复杂性，理论研究并不能满足工程需要，物理水池试验研究成本高、周期长，随着数值计算方法发展日趋成熟，对于水上飞机耐波性的数值计算研究逐渐展开。刘沛清等[4]对边界元法、有限元法、光滑粒子水动力学法和有限体积法等在近期飞行器水上迫降中的具体应用进行了总结；Newman[5]基于线性波理论，研究了造波及消波的理论方法，建立了方形和圆形波浪水池数值模型；王硕等[6]基于VOF RANS求解器，对高速滑行艇在规则波中的升沉运动规律进行了数值研究。

本文基于LS-DYNA，通过数值试验的手段研究了水陆两栖飞机在波浪水面上降落时的耐波性能，首先采用基于微幅波理论的动边界数值造波法实现了对动态行进波浪的模拟，并通过数值造波算例对本文中所使用的数值造波方案进行了验证；随后，采用ALE及罚函数流固耦合方法，对水陆两栖飞机迎波条件下降落时的纵摇运动姿态角变化、升沉过载情况以及底部压力分布等运动学和动力学特性进行了研究；最后，分析了入水波浪的波高及波长等不同海况条件对水陆两栖飞机耐波性能的影响。

2 数值计算分析方法

参考真实物理水池中的推波板造波方法，图1为本文动边界造波法示意图，水域左侧动边界起到物理水池中推波板的作用，做冲程为S的简谐振动，扰动边界处水域，生成波浪并向前传播。右侧和底部为固壁边界，上部为自由液面。基于微幅波理论，对动边界造波法进行理论推导[7]，可以得到波浪波高和波长与动边界冲程之间的关系，表达式为

H=2S(cosh2kph-1)/(sin2kph+2kph)

σ2=gkptanhkph，L=(gT2/2π)tanhkph

(1)

由于计算成本限制，计算水域长度不能无限大，在水域右侧固壁边界会产生一定的反射波与正向波叠加，对规则波造波结果造成不良影响，进而影响波浪与飞机结构的耦合计算结果。波浪运动中水质点的振动属于低频运动，本文在水域右端 1～2倍波长的消波区域内的水域设置质量阻尼，水质点的运动速度衰减，从而达到吸收反射波浪能量的目的[8]。

(2)

式中k为罚函数刚度系数，c为罚函数阻尼系数，耦合力同时作用在主物质及从物质节点上，方向相反，保证了接触面上力的平衡，阻止相互穿透。在 LS -DYNA 中，接触刚度是以单位面积定义的，接触刚度和弹簧阻尼系统的运动方程表示为

(3，4)

合适的Pf及ξ值应保证接触面上的压力变化相对平滑，又不会显著改变接触系统的动力学特性，流体与结构的相互穿透量应足够小以至可以忽略不计。根据以上原则，本文确定Pf=0.06，阻尼因子ξ=0.0，即不考虑阻尼力作用。

图1 动边界造波法示意图

Fig.1 Illustration of wave -making method by plane wavemakers

3 数值计算模型

为了验证本文基于推板式动边界造波法的正确性，参考文献[11]的试验条件，建立了数值造波验证算例，得到推波板振动周期与波浪波高关系拟合曲线，如图2所示。数值计算结果与试验结果符合性良好，当推波板周期小于1 s时，数值造波所生成波浪波陡较大，以至波浪破碎；周期在1.0 s～2.5 s之间，误差小于10%；周期在2.5 s～3.0 s之间时，数值计算结果有下降趋势。水域高度及推波板冲程一定时，波浪波长随推波板周期增加而增加，波浪传播速度降低，传播相同波长个数，时间增加，波高会受到水域自身重力加速度影响，周期在2.5 s～3.0 s之间时波高有一定程度衰减。

参考文献[12]的试验条件建立算例模型，得出了推波板冲程与生成波浪波高的关系，数值计算值、试验值和理论值对比如图3所示，可以看出，本文数值结果与试验结果一致，误差在5%以内。

本文研究的水陆两栖飞机机身长为35 m，宽

图2 推波板振动周期与波高关系图(S=0.3 m)

Fig.2 Relationship between period of wavemaker and wave height (S=0.3 m)

图3 推波板冲程与波高关系图(T=2.0 s)

Fig.3 Relationship between stroke of wavemaker and wave height (T=2.0 s)

为3.22 m，翼展为38 m，起飞重量为50 t。相比支线客机，水陆两栖飞机有能够减少水阻力的尖劈形船底，保持飞机横向稳定性的翼下浮筒，方便水上飞机脱离水面的机腹断阶及方便运输的上单翼等结构，飞机结构特点及参数如图4所示。

飞机的计算有限元模型选择为实机尺寸的整体刚体模型，采用朗格朗日壳单元进行划分，在LS -DYNA中给定刚体飞机质量、质心位置及转动惯量，保证飞机刚体有限元模型与真机质量分布的一致性。

本文对气动载荷进行了等效简化处理，首先通过计算流体动力学软件(CFX)进行飞机气动力的计算，得到飞机在初始水平速度及姿态角下重心处的升力、阻力和俯仰力矩，如图5所示。接着，通过定义集中力载荷曲线，施加以上等效气动力在飞机水上降落数值计算模型中的飞机模型质心处，并考虑为在飞机触水前线性递减为0。这种简单加载方法的有效性在数值模拟工作中已经得到了验证，数值结果与试验结果保持了比较好的一致性[13]。

根据着水试验数据，飞机抗浪波高为1.2 m，起飞重量水陆两栖飞机抗浪波高参考值为1.6 m，飞机纵摇谐振波长在1～3倍机身长度范围内[3]，为研究波高及波长对飞机耐波性的影响，本文选取0.8 m，1.2 m和1.6 m三个有效波高等级工况，1、2和3倍机身长度L0波长工况，本文工况列入表1。

图4 水陆两栖飞机的几何外形特点

Fig.4 Geometric characteristics of amphibious aircraft

图5 机身的气动载荷作用示意图

Fig.5 Illustration of the aerodynamic load of aircraft

依据着水试验经验，飞机着水后滑行7～8倍波长距离即可判断飞机着水动力特性的优劣，且本文模拟波高在0.8 m～1.6 m之间，根据式(1)，确定水域高度为8 m，长度为8倍波长，其中消波区长度为1倍波长，作为空间预设的空物质域高度设为4 m。上单翼在着水过程中不触水，水域宽度设置为24 m。

合理的网格尺寸，首先应保证飞机气动外形不至失真，水域网格要与飞机网格尺寸大致相同才不至于在耦合计算过程中出现负体积等问题。其次，LS -DYNA计算采用显式算法，计算时间与单元数量成正比，与最小单元的尺寸成反比，本文中典型工况单元数约为160万，单元尺寸微小的改变会使得单元总数成千上万地增加，使得计算成本大大增加；而且，在满足以上原则的单元尺寸基础上减小单元尺寸后，计算结果没有太大差异，仅时间历程结果更加平滑，对结果分析没有太大意义[15]。

根据以上原则本文确定飞机壳单元最大尺寸为0.3 m，流体域网格约为0.5 m，流体域网格采用渐进式网格进行划分。按上述内容建立迎波条件下水陆两栖飞机耐波性数值计算模型，如图6和图7所示。波浪传播到消波区域后，飞机以水平速度50 m/s，垂直下降速度1.5 m/s降落，初始姿态角为7°。

表1 工况计算参数

Tab.1 Parameters of design conditions

工况编号动边界圆频率σ动边界冲程S/m有效波高H/m波长L/m1π/41.30.82L02π/42.01.22L03π/42.71.62L04π/41.41.2L05π/44.11.23L0

图6 耐波性数值计算模型

Fig.6 Numerical simulation model of seakeeping investigation

图7 水陆两栖飞机波浪面降落过程模拟

Fig.7 Simulation of amphibious aircraft landing on waves

4 结果分析

飞机姿态角的变化反映了机体纵摇运动的剧烈程度，图8为各工况飞机纵摇姿态角结果。飞机降落在1.2 m波高下，纵摇幅值大于0.8 m波高工况，1.6m波高工况飞机做掠过波峰的滑行运动；飞机在2L0工况下，纵摇幅度大于1L0和3L0工况，1L0和2L0工况下飞机跳波现象明显，跳波现象改变了飞机的遭遇频率，为所跨波峰个数与波浪频率之积。

飞机在波浪上滑行，运动响应在达到最大值时遭遇波浪的波长和频率称为谐振波长和谐振频率，与飞机固有频率、波幅和飞机速度相关。本文中飞机以50 m/s的速度降落滑行时，飞机纵摇运动谐振波长约为2倍机身长度，谐振频率为0.75 Hz左右，波长小于谐振波长时飞机做跳波运动。根据试验结果，飞机模型在该滑行速度下，飞机谐振波长在1.5～3.5倍机身长度范围内，谐振频率在0.8 Hz～1.5 Hz之间，本文计算结果与试验结果较为符合[3]。

飞机升沉运动过载也是判别飞机耐波性的重要指标，过高的升沉运动过载会导致机乘人员身体

图8 飞机纵摇姿态角变化结果

Fig.8 Results of attitude angle changes in pitching motions

不适及工作能力降低等。图9为飞机质心处垂向过载变化结果，其规律与飞机纵摇运动规律一致，升沉过载随波高增高而增大，升沉运动谐振波长约为2倍机身长度波长，谐振频率约为0.75 Hz。飞机由于受到水载荷产生的上升过载大于飞机在重力和伯努利效应作用下的下沉过载，过载最大值为1.11 g。

为研究水陆两栖飞机在波浪面上降落过程中机乘人员与物资所承受的过载情况，本文选取了机艏和机艉两个过载采集点，如图10所示。机艏和机艉处的加速度由飞机刚体平动及转动产生的过载叠加组成，按下式计算得出，

机艏：a艏τ=axcosθ-azsinθ+ω2L艏

a艏⊥=axsinθ+azcosθ+αL艏

机艉：a艉τ=axcosθ-azsinθ-ω2L艉

a艉⊥=axsinθ+azcosθ-αL艉

式中aτ和a⊥分别代表切向和垂向加速度，ax,az,ω和α分别为刚体飞机质心处的水平加速度、垂直加速度、转动角速度和转动角加速度。

机艏和机艉处的过载变化结果如图11所示，结果显示，机艏过载受波高影响不明显，受波长影

图9 飞机质心升沉过载变化结果

Fig.9 Results of overload changes in heaving motions

图10 机艏、机艉过载的测点位置

Fig.10 Station of fore -body & rear-body

图11 机艏、机艉处过载变化结果

Fig.11 Results of overload changes at fore -body & rear-body

图12 机身底部压力峰值分布

Fig.12 Pressure peak distribution at the bottom of the fuselage

响较大，波长越短，机艏过载越高。在波长小于谐振过载情况下，机艏过载随时间呈上升趋势，这是由于机艏破浪作用承受了大量的波浪能量，机艏当前过载与遭遇之前波峰时的残余过载叠加而产生上升趋势。

机艉过载在不同海况条件下的变化规律与机体纵摇和升沉运动结果规律一致，过载最大情况在谐振波长海况条件下出现，断阶的整流作用使得飞机未出现吸附后体的状况。

如图12所示，飞机底部冲击压力峰值分布受海况条件影响不明显，压力集中在船型机身前体扭曲处，尖劈形扭曲前体底部承受了大部分的波浪冲击能量作用，提高了飞机纵向安定性和喷溅特性。断阶消除了伯努利效应，使得后体负压转变成正压。同时，飞机刚体转动使得机艉底部触水，机艉底部也有一定压力分布。

5 结 论

水陆两栖飞机在波浪水面上降落的耐波性物理试验复杂，成本巨大，数值计算是目前最经济有效的分析方法。本文实现了推波板造波理论的数值计算，模拟了水陆两栖飞机在波浪水面上降落的过程，得到结论如下。

(1) 水陆两栖飞机在波浪面上滑行时，纵摇与升沉周期性运动在时间历程上除相位外具有同步性，当波浪波长小于飞机谐振波长时，飞机做跳波运动，运动轨迹波长逼近谐振波长。

(2) 机艏承受了大量的波浪能量，机艏过载与压力均较大，此处机乘人员与机体的承受能力是衡量飞机耐波性的重要指标。

(3) 本文水陆两栖飞机耐波性的数值计算方法可行，计算结果符合相关物理实验与数值计算结果，且高效可靠。

真实海洋环境中的波浪是非规则的、三维的，对飞机侧向平动及滚转运动有一定影响。飞机降落过程中空气动力影响着飞机在入水冲击及滑行过程中的运动特性，地效效应影响尤为重要。另外，弹性体飞机的变形吸能情况和飞机初始入水条件(水平速度、下降速度等)对飞机耐波性的影响也不容忽视。

:

[1] Hugli W C,Axt W C.Hydrodynamic investigation of a series of hull models suitable for small flying boats and amphibians[J].TechnicalReportArchive&ImageLibrary,1951,31(8):101-129.

[2] Benson J M ,Havens R F,Woodward D R.Landing characteristics in waves of three dynamic models of flying boats[J].TechnicalReportArchive&ImageLibrary,1952.

[3] 黄 淼,吴 彬,蒋 荣,等.水上飞机在波浪上运动响应特性试验研究[J].实验流体力学,2015，29(3):41-46.(HUANG Miao,WU Bin,JIANG Rong，et al.Experimental study on motion response of a seaplane on waves[J].JournalofExperimentsinFluidMechanics,2015,29(3):41-46.(in Chinese))

[4] 刘沛清,屈秋林,郭保东,等.数值计算技术在飞机水上迫降中的应用[J].力学与实践,2014,36(3):278-284.(LIU Pei-qing,QU Qiu-lin,GUO Bao -dong,et al.Application of computational fluid dynamics in the planned ditching of a transport airplane[J].MechanicsinEngineering,2014,36(3):278-284.(in Chinese))

[5] Newman J N.Analysis of wave generators and absorbers in basins[J].AppliedOceanResearch,2010,32(1):71-82.

[6] 王 硕,苏玉民,庞永杰，等.高速滑行艇在规则波中的纵向运动数值研究[J].哈尔滨工程大学学报,2014,35(1):45-52.(WANG Shuo,SU Yu-min，PANG Yong-jie,et al.Numerical study on longitudinal motions of a high-speed planing craft in regular waves[J].JournalofHarbinEngineeringUniversity,2014,35(1):45-52.(in Chinese))

[7] Dean R G,Dalrymple R A.WaterWaveMechanicsforEngineersandScientists[M].Prentice -Hall,1984.

[8] 韩 朋,任 冰,李雪临,等.基于VOF方法的不规则波数值波浪水槽的阻尼消波研究[J].水道港口,2009,30(1):9-13.(HAN Peng,REN Bing,LI Xue -lin,et al.Study on damping absorber for the irregular waves based on VOF method[J].JournalofWaterwayandHarbor,2009,30(1):9-13.(in Chinese))

[9] 张 苏,古 彪,曹东风,等.基于ALE算法的V形楔形体入水的水动力特性分析[J].固体力学学报,2014,35(s1):95-100.(ZHANG Su,GU Biao,CAO Dong-feng.Hydrodynamic behavior analysis of water entry of V-shaped wedge based on ALE method [J].ChineseJournalofSolidMechanics,2014,35(s1):95-100.(in Chinese))

[10] 浦锡锋,王仲琦,白春华,等.入水问题的Euler-Lagrange流固耦合数值模拟技术研究[J].计算力学学报,2013,30(s1):159-163.(PU Xi-feng,WANG Zhong-qi,BAI Chun-hua.An Euler-Lagrange coupling method for numerical simulation of water entry[J].ChineseJournalofComputationalMechanics,2013,30(s1):159-163.(in Chinese))

[11] 刘月琴,尹尚先.实验室造波设备及波浪模拟技术的探讨[J].华南理工大学学报(自然科学版),2000,28(6):105-109.(LIU Yue -qin,YIN Shang-xian.A tentative study of wave -making equipment and wave simulated technology in the laboratory[J].JournalofSouthChinaUniversityofTechnology(NaturalScienceEdition),2000,28(6):105-109.(in Chinese))

[12] 路 宽.实验室水槽造波控制系统的开发[D].天津大学,2007.(LU Kuan.The Development of Laboratory Wave Flume Control System [D].Tianjin University,2007.(in Chinese))

[14] 张松涛,梁利华,孙明晓,等.耦合船体横摇的减摇水舱流体运动和减摇效果预测[J].计算力学学报,2016,33(2):252-256,262.(ZHANG Song-tao,LIANG Li-hua.SUN Ming-xiao,et al.Fluid motion and stabilization effect prediction of anti-rolling tank coupled ship rolling [J].ChineseJournalofComputationalMechanics,2016,30(2):252-256,262.(in Chinese))

[15] 许 栋,徐 彬,David PAyet,等.基于GPU并行计算的浅水波运动数值模拟[J].计算力学学报,2016,33(1):113-120.(XU Dong,XU Bin,David PAyet.Numerical simulation of shallow water motion based on parallel computation using GPU[J].ChineseJournalofComputationalMechanics,2016,30(1):113-120.(in Chinese))

[16] Cadilha C,Inês M.Numerical Simulation of Aircraft Ditching of a Generic Transport Aircraft:Contribution to Accuracy and Effciency [D].Instituto Supe -rior Técnico,2014.