基于CEL算法的水陆两栖飞机水上降落动力特性分析*

姚小虎 黄愉太 欧智成 刘晓明 胡伯仁

(1.华南理工大学 土木与交通学院，广东 广州510640;2.成都飞机工业(集团)有限责任公司 技术中心，四川成都610092)

水陆两栖飞机是一种能在水面上起飞和降落的特殊多用途飞机［1］，如水上轰炸机、水上侦查机、水上运输机、水上战斗机、海上救护机、灭火用水上飞机等，随着我国航空事业的快速发展，水陆两栖飞机也在飞机大家族中占据一席之地.

水陆两栖飞机的水上降落属于典型的流固耦合问题［2］，入水时水对机身船体的冲击使船体发生变形，而船体变形反过来又将影响水的流动.入水问题的研究早期以实验［3-4］和理论［5］为主，随着计算机的高速发展，数值模拟也越来越多地应用于入水问题.Brooks 等［6］最先使用LS-DYNA 模拟Apollo 指令舱撞击水面的过程. 卢鸣飞等［7］采用FEM-SPH 算法对飞机的入水进行了分析. 刘翔［8］建立了A320飞机模型，采用ALE/Penalty Coupling 方法成功模拟出飞机迫降的流固耦合作用，并将计算结果与实验进行了对比.张韬等［9］运用Dytran、采用一般耦合算法计算了飞机水上迫降的冲击荷载，并将结果与刚体飞机模型进行了对比.Chandra 等［10］对2009年在美国发生的一起飞机成功水上迫降的例子采用FEM-SPH 算法进行了分析，分析中将飞机视为刚体，结果证明采用该方法进行水上迫降计算是合理的.但是目前对于水上降落问题研究的大部分文献都是将机身处理成刚体［11］，且只给出定性结果而没有做定量的研究，而在实际入水过程中，考虑机身的冲击变形对提高结果的精确度是非常有必要的［12］.

文中采用耦合的欧拉－拉格朗日(CEL)算法对水上飞机降落的动力过程进行了研究，分析了某型水陆两栖飞机机头水上降落过程中不同姿态角和垂直速度对飞机结构动力响应的影响，并得到适合该机型水上降落的参考角度和速度，为该型水陆两栖飞机水上降落提供理论指导依据.

1 耦合欧拉－拉格朗日法

拉格朗日(Lagrangian)方法主要应用于固体力学中的受力和变形分析，在计算过程中，所分析物体的形状变化与有限元网格的变化完全保持一致，其优点在于它可以大大简化控制方程的求解过程，并可跟踪质点的运动轨迹，从而能准确地描述物体边界的运动;但是在涉及特大变形问题时，拉格朗日算法采用的单元网格将会出现严重的畸变现象，以致无法继续计算.欧拉方法(Eulerian)主要应用于流体力学的计算分析，其采用的有限元网格形状、大小和空间位置保持不变，在计算过程中材料可在固定的有限元网格中流动，因此欧拉算法能够处理物质的扭曲及一些大变形问题，但欧拉方法在捕捉物体边界信息上较为困难，不能够精确描述物质的边界［13］.

基于以上方法特点，Noh［14］首先提出了耦合的欧拉－拉格朗日算法(CEL)，它有机地结合拉格朗日方法和欧拉方法的优点，克服其中的缺陷，充分发挥各自的优势. 在CEL 方法中，固体变形小则采用拉格朗日单元，而流体变形大时采用欧拉单元. 因此，CEL 算法在解决物体的大位移时有较大优势，比如碰撞、流体动力学及流固耦合作用等.

2 水上飞机机头降落动力学模型

2.1 飞机机头入水有限元模型

机头包括雷达罩、驾驶舱、生活舱、船体4 部分，其中船体为入水首要冲击区域.建立机头入水模型，如图1 所示.机头蒙皮和腹板均采用壳单元，船体处加强件采用梁单元. 壳单元主要采用四边形和三角形单元，船体底部单元局部加密.水域模型为椭圆柱体，其长轴为9 m，短轴为6 m，高3 m，单元类型为六面体的欧拉单元，考虑到溅水，欧拉单元的区域要略大于初始水域，因此初始水域深2.7 m，欧拉单元表面在初始液面上0.3 m.开始时船底与水域相切，在入水位置附近加密网格，如图2 所示.

图1 机头模型Fig.1 Aircraft head model

2.2 材料本构

机头各部件均采用铝板，材料分布及物理参数如表1 所示，各型号铝材采用弹塑性双线性本构模型进行描述:

图2 机头入水模型Fig.2 Model of the head hitting the water

式中，σ 和ε 为应力和应变，σ0和ε0为屈服应力和屈服应变，E 为弹性模量，ET为强化模量.

表1 机头各部件材料物理参数1)Table 1 The material physical parameters of head

受冲击流体的材料行为极其复杂，可使用状态方程来描述材料的冲击力、体积和能量特性.各向同性冲击力可由Mie-Gruneisen 状态方程［15］得到，即

式中:Im为单位质量的内能;η =1 －ρ/ρ0，为名义体积压缩率，ρ0表示初始材料密度;Г0为Gruneisen 常数.激波速度Us与粒子速度Up之间有线性关系，Us=c+Up.在Us-UpHugoniot 模型中有PH=ρ0c2η·(1 －sη)－2，其中s 为极限的体积压缩率，c 为材料中的声速，当名义应变较小时，体积模量K=ρ0c2.

文中采用线性的Us-UpHugoniot 状态方程描述水的特性，取参数s=Г0=0，声速c=1483 m/s，密度ρ=1000 kg/m3，黏度μ=0.001 kg/(m·s).

3 机头入水结果分析

文中考虑机头在不同的入水角度和垂直入水速度下的落水情况，共5 种工况，如表2 所示.

表2 机头入水工况Table 2 Cases of head attacking the water

飞机降落时受到竖直向上的气动力2/3 g，与自重叠加，可近似认为飞机降落时有向下的1/3 g 加速度.水域施加一般重力场.水域边界为自由流入和无反射流出.所有工况的机头水平速度均为41 m/s，计算时间为0.1 s.

3.1 工况1 结果分析

图3 是工况1(入水角为7°，垂直入水速度为3m/s)下机头应力和水域在入水过程中的变化.从云图可看出，机头入水时受力区域主要为船体，机身部分基本无影响，因此重点分析船体的结构动力响应.机头从0s 开始与水面接触，在入水后的0.016 s 出现最大应力，为550MPa，位于船体骨架的纵向腹板处，超过材料屈服强度且接近抗拉强度，说明纵向腹板受力明显，危险性较高，机头其他部位受力相对较小，均未超其屈服强度.随着入水深度增加，机头船底的受冲击区域逐渐增大，入水过程中机头逐渐抬头.水域被划出一道水槽，影响区域逐渐增大，并伴随液体飞溅.

图3 工况1 下机头应力和水域的变化Fig.3 Head’s stress and water variation in case 1

3.2 不同垂直速度下入水的结果分析

当飞机入水角为7°，垂直速度分别为3、2、1m/s时，机头应力、船底冲击力和质心速度的变化情况如表3 和图4 所示.

表3 机头不同垂直速度入水结果Table 3 The results of head attacking the water at different vertical speeds

通过对3 种工况的比较可知，机头应力都是在入水初期(0.02 s 左右)出现最大值. 对于船壳最大应力，3 种工况相差不大，均保持在430 ～450MPa 之间，满足材料强度要求.机头船体骨架应力随垂直速度变化明显，在3 和2 m/s 时其最大应力甚至超过了船壳应力，分别为550 和529 MPa，接近材料抗拉强度，1m/s 时骨架应力迅速减小为337MPa，远小于抗拉强度，三者出现位置均位于船体骨架的纵向腹板处，可见腹板受垂直速度影响明显.随着垂直速度的减小，船体骨架的应力也减小，这是因为随着入水深度增加，速度越大，船体骨架尤其纵向腹板的受挤压程度也越大，容易发生局部屈曲和变形.当入水角度相同时，垂直速度越大，船底冲击力越小. 在入水初期，船底底面法向与水域垂直，此时底面受到的冲击力出现峰值，而后冲击力逐渐衰减，随着机头慢慢抬头，受冲击面法向发生变化，受水平速度影响，冲击力再次出现峰值.从质心速度减小量可看出，入水过程中机头整体保持稳定(水平晃动速度v2极小)，垂直速度(v3)越大，入水深度也越大，从而阻力增加，水平速度(v1)减小得也越多. 综合对比3 种工况，当入水速度为1 m/s时机头入水最安全.

3.3 不同角度下入水的结果分析

飞机在入水的垂直速度为3 m/s，入水角度分别为7°、6°、5°时的结果如表4 和图5 所示.

图4 不同垂直速度下入水的机头应力、船底冲击力和质心速度时程曲线Fig.4 Time-history curves of head stress，bottom pressure and mass centre of the head at different vertical speeds

表4 机头不同角度入水结果Table 4 The results of the head attacking the water at different angles

机头入水初期应力迅速增大，在0.02 s 左右出现最大值，并随入水角度减小而减小.3 种工况下，船壳的最大应力相当，变化不大，均在430 ～450MPa之间波动，满足材料强度要求，但对于船底骨架，最大应力变化明显，在7°和6°时超过船壳分别达到了550 和536 MPa，出现在纵向腹板处，接近抗拉强度，为最危险点.在入水角为5°时减小到365 MPa，满足强度要求.入水瞬间船底冲击力出现峰值，然后衰减，随着机头抬头，冲击面发生变化，冲击力再次出现峰值，当入水角度为6°时，船底的冲击力最大，这是因为船底为曲面，6°时船壳底面法向方向更接近竖直方向，因此入水的速度对于船底法向方向的分量越大，冲击力也越大.入水过程中，6°时船底底面与水面更接近平行，产生阻力相对较小，其水平速度(v1)较其他两工况减小量更少，对于垂直速度(v3)，由于入水深度大，7°时减小明显，3 种工况下机头整个过程都保持稳定(v2极小). 综上可得出，当飞机水上降落时入水角度为5°的安全性更高.

4 结论

文中建立了某型水陆两栖飞机机头水上降落流固耦合模型，采用CEL 算法分析了机头在不同入水角度和垂直入水速度下的结构响应. 船壳的应力在入水初期迅速增大，之后处于波动状态.各工况下船壳最大应力值相差不大，而船体骨架受力变化明显，其最大值大部分情况下甚至超过船壳，且都出现在纵向腹板处，因此在结构设计中应加强骨架尤其是纵向腹板处的结构强度. 机头入水瞬时船底冲击力即出现峰值.入水速度减小，机头应力和船底冲击力都随之减小;随着入水角度的减小，机头应力也减小，船底冲击力在6°时最大. 因此，根据计算结果，当入水角度为5°、垂直入水速度为1 m/s 时，飞机水上降落相对最安全.

图5 不同入水角度下机头应力、船底冲击力和质心速度时程曲线Fig.5 Time-history curves of head stress，bottom pressure and mass centre of the head at different attacking angles

［1］尹士达，林明.我国水上飞机的特点和主要性能［J］.江苏航空，1998(3):12-13.Yin Shi-da，Lin Ming. The main characteristics and performance of seaplane in China ［J］. Jiangsu Aviation，1998(3):12-13.

［2］张阿漫，戴绍仕.流固耦合动力学［M］.北京:国防工业出版社，2011:49-85.

［3］Thomas W L.Ditching investigation of a 1/20-scale model of the space shuttle orbiter ［M］. Bethpage:Grumman Aerospace Corporation，1975.

［4］Sawan A Shah. Water impact investigations for aircraft ditching analysis［D］.Melbourne:Royal Melbourne Institute of Technology，2010.

［5］Von Karman T.The impact of seaplane floats during landing［C］∥Proceedings of NACA TN 321，National Advisory Committee for Aeronautics. Washington:［s. n.］，1929:1-8.

［6］Brooks J R，Anderson L A. Dynamics of a space module impacting water ［J］. Spacecraft and Rockets，1994，31(3):509-515.

［7］卢鸣飞，华诚.基于FEM-SPH 耦合算法的飞机水上迫降的力学建模与分析［J］. 复旦学报:自然科学版，2013，52(5):583-591.Lu Ming-fei，Hua Cheng.Model and numerical simulations of water ditching by FEM-SPH coupling method ［J］.Journal of Fudan University:Natural Science，2013，52(5):583-591.

［8］刘翔.飞机水上迫降的运动特性分析和讨论［D］. 武汉:武汉理工大学理学院，2012.

［9］张韬，李书，江翔，等.民用飞机水上迫降分析模型和数值仿真［J］.南京航空航天大学学报，2010，42(3):392-394.Zhang Tao，Li Shu，Jiang Xiang，et al.Analysis model and numeral simulation for civil plane ditching［J］.Journal of Nanjing University of Aeronautics ＆ Astronautics，2010，42(3):392-394.

［10］Chandra C，Wong T Y，Bayandor J.Crash worthiness assessment in aircraft ditching incidents［C］∥Proceeding of 27th Congress of International Council of the Aeronautical Sciences.Nice:［s.n.］，2010:1-13.

［11］贺谦，陈效鹏，李磊，等.飞机水面降落的机身载荷研究［J］.强度与环境，2009，36(4):28-33.He Qian，Chen Xiao-peng，Li Lei，et al.Research on the load of aircraft landed on water［J］. Structure ＆ Environment Engineering，2009，36(4):28-33.

［12］宋长福.民机机身结构入水冲击问题数值仿真研究［D］.上海:上海交通大学航空航天学院，2011.

［13］徐文杰.基于CEL 算法的滑坡涌浪研究［J］.工程地质学报，2012，20(3):350-354.Xu Wen-jie. CEL algorithm study of reservoir surge induced by landslide［J］.Journal of Engineering Geology，2012，20(3):350-354.

［14］Noh W F.CEL:a time-dependent two-space-dimensional coupled Eulerian-Lagrangian code［M］.New York:Academic Press，1964:117-179.

［15］Shyue Keh-ming.A fluid-mixture type algorithm for compressible multicomponent flow with Mie-Gruneisen equation of state ［J］. Journal of Computational Physics，2001，171(2):678-707.