基于动网格的两栖车航行姿态数值模拟

赵彬，张敏弟，剧冬梅

（1.北京理工大学机械与车辆学院，北京100081；2.中国兵器科学研究院，北京100089）

基于动网格的两栖车航行姿态数值模拟

赵彬1，张敏弟1，剧冬梅2

（1.北京理工大学机械与车辆学院，北京100081；2.中国兵器科学研究院，北京100089）

针对两栖车水上航行姿态,构建描述两栖车水上运动的动力学模型,采用混合耦合算法和动网格技术研究静水直航状态下,两栖车航行姿态的变化规律,模拟结果与实验值吻合较好。研究结果表明：车体达到稳定航行姿态会经历大幅振荡调整阶段和平稳运动阶段,车体平衡后,升沉和纵倾仍会有小幅变化；两栖车由排水航行状态逐步增速到滑行状态的过程中,车体重心逐渐升高,动升力在支撑车重成分中所占的比例越来越大,静浮力则越来越小。

兵器科学与技术；两栖车；航行姿态；动网格技术；数值模拟

0　引言

两栖车在水上实际航行过程中伴随着升沉、摇摆等复杂运动，其水上航行特性对其安全性能的保证具有重大作用［1-4］。对于两栖车水上动态性能的数值计算，主要有两种计算方法：第一种是静态网格技术，第二种是动态网格技术。

关于静态网格技术，近年来学者已经做了大量的工作。韩占忠等［5］对两栖车的黏性阻力进行数值计算，并对车首防浪板的设计进行了计算分析；王涛等［6］分析两栖车各阻力成分随速度变化的规律，为分析两栖车水上快速性提供了依据；吴珂等［7］运用计算流体力学软件Fluent，对两栖车车轮收起前后的阻力进行了仿真分析。静态网格技术的航行姿态调整方法，主要是手动调整方法。对两栖车航行姿态的数值模拟主要是采用车体相对于流场静止不动，赋予来流一定的速度，从而得出车体的绕流场情况。根据车体受力调整航行姿态，直到航行姿态满足要求。这种模拟方法的不足是：工作量大，车体静止不动，没有考虑浮态的变化。

由于静态网格技术的局限性，动态网格技术应运而生。两栖车在水中的运动是一动态过程，采用动网格模型可以有效解决这一问题。相比于静态网格技术，两栖车动态网格技术的研究国内还较少。徐国英等［8］通过网格分域方法，对两栖车进行了阻力和浮态的数值模拟；孙伟等［9］通过应用动网格技术和VOF模型，计算了两栖车在一定速度的航态和水动力特性；万晓伟等［10］提出可以通过任意拉格朗日-欧拉方法来解决流-固动力耦合问题，从而获取车辆最终航态。动网格技术主要是对车体实施主动或被动运动，从而实现车体的动态运动过程，并在过程中监测车体姿态的变化，该方法可充分利用计算机的资源，根据每一个时间段的受力状态，基于软件对车体的姿态进行自动调整，但网格的精度不易于保证。

两栖车在静水中航行时，由于受到外力的作用，会产生6个自由度的运动，这几种运动之间存在非线性耦合作用。针对两栖车水上航行姿态，本文构建了两栖车水上航行数学模型，采用动网格技术研究了静水直航状态下，两栖车航行姿态的变化规律。

1　数学模型和数值模拟方法

1.1研究对象

本研究采用简化的两栖车车体为研究对象，车质量为6.5 t，在水中行驶时车轮处于提升状态，简化三维模型如图1所示。

图1　简化两栖车模型Fig.1 Modified model of amphibious vehicle

1.2湍流模型

两栖车水上航行属于典型的流-固耦合运动，不仅伴随水和空气两种黏性流体的相互作用，还伴随着水和车体、空气和车体之间的相互作用。车体在水中运动时，由于水和空气两种流体是有黏性的，所以车体受力满足黏性流体力学的基本方程。

湍流模型选取标准k-ε模型［11］，k和ε的输运方程为

式中：t为时间；uj和xj分别为速度分量和坐标分量；ρm和μ分别为密度和分子黏性系数；μt为湍流黏性系数，μt=Cμρmk2/ε，其中Cμ为经验常数，Cμ=0.09；σk和σε分别为k和ε的湍流普朗特数，σk=1.0，σε=1.3；Cε1和Cε2为ε方程常数，Cε1=1.44，Cε2=1.92；Pt为平均速度梯度引起的湍动能k产生项，

1.3两栖车动力学模型

满足两栖车升沉和纵倾两个自由度的三维模型见图2.图中：L表示车体y轴所受升力，y轴正向为正；M表示车体z轴所受转动力矩，正负方向遵守右手法则。

图2　满足h和θ两个自由度的两栖车三维模型Fig.2 A 3-DOF model representing h and θ degrees of freedom of an amphibious vehicle

升沉和纵摇的运动幅度通过重心升沉h和纵倾角θ来定义。整个系统的控制方程［12］为

式中：ms为车体的质量矩阵；Cs为车体的阻尼矩阵；Ks为车体的刚度矩阵；｛Mex｝为外部激励矩阵；为在结构求解器计算出网格变形后通过使用计算流体力学（CFD）求解器计算出来的非线性黏性流体表面力；｛x｝为车体的变形量矩阵。

1.4流-固耦合混合算法

流-固耦合问题研究的是流体和固体结构之间的相互作用，包括固体结构在流体载荷作用下产生的变形或运动及其对流场的影响，因此，求解流-固耦合问题，需要同时考虑流场和结构场的求解及其耦合［13-15］。一般用到4种数值解法［16-17］，分别为完全耦合（FC）算法、松式（LC）算法、紧式（TC）算法、混合耦合（HC）算法。FC算法是完全耦合算法，对流体和结构建立统一的耦合方程，在一个时间步内对流体域和固体域中所有的未知量进行求解，求解精度较高，但在工程应用上求解范围较窄；LC算法和TC算法将模型的结构和流体部分独立进行求解，通过交互程序实现结构求解器和流体求解器之间的数据交换，但求解误差较大；HC算法在（3）式左右两边增加附加质量力项，适当预测了流体对运动的惯性、阻尼以及刚度，减小了对物体运动位移的过预测，从而提高了数值算法稳定性。相比之下LC算法和TC算法需要更小的时间步长或者更长的计算时间，甚至计算发散。鉴于此，本次计算采用的是HC数值算法。HC算法离散后的方程为

式中：n为时间步增量；i为每一个时间步的子迭代步数；｛Mfluid｝是附加质量力项。｛Mfluid｝［16］假设为

式中：2b为两栖车的特征长度，2b=5.04 m；ρf为流体密度。

（5）式的求解采用2阶精度的半隐式Crank-Nicholson方法。

本次计算中，忽略车体本身的阻尼项，同时给出车体的质量矩阵和刚度矩阵，根据每一步计算得到的车体升力L和作用力矩M，得到车体升沉和纵倾的变化规律。

混合耦合算法的流程图如图3所示。其中：δ是一个介于0～1之间的调节系数，δ=0.5；ε是对于子迭代的容差值，ε=10-5；上标i表示相同时间步的子迭代步数；imax是子迭代的最大步数；下标n表示不同的时间步；x0和分别是初始变形量和变形速度。计算时，将imax值设为1在第一个子迭代步（i=1），附加质量力是从上一个时间步的最后一个子迭代步的变形算起的：

混合耦合算法具体步骤如下：

首先，从t=0，子迭代步为1开始计算，假设初始变形量x0，初始变形速度x·0，通过CFD软件计算出，并将其代入（5）式中进行迭代求解，解得下一子迭代步的变形量，变形速度以及变形加速度；然后，将相邻子迭代步求解的变形量进行容差比较，若满足要求，则输出变形结果，进行下一时间步的迭代计算，若不满足，则重新进行（5）式的子迭代计算得到，直到满足容差要求为止，得到，并输出结果；最后，当计算达到预设定总时间时停止计算，输出模型变形结果。

图3　混合耦合迭代算法具体流程Fig.3 Flow chart of hybrid algorithm

2　网格划分和边界条件设置

2.1网格划分及网格无关性验证

网格划分采用全结构化网格，网格总数约为100万。在车体周围进行网格加密，并设置边界层，网格的划分见图4.

图4　车体附近网格划分Fig.4 Mesh near amphibious vehicle

网格无关性从两栖车辆阻力系数Cd进行分析，Cd定义为

式中：R为车身所受阻力；ρ为水密度；v为车体航行速度；V为排水体积。通过计算，获得网格数量n与阻力特性关系如图5所示。可以看出，使用阻力系数作为网格无关性评判标准时，网格数量在100万以后，阻力特性变化不大，鉴于此，文中所使用的数值模拟网格约100万。

图5　网格无关性验证Fig.5 Grid-independent validation

2.2边界条件设置

如图6所示，计算区域以车体为中心向前延伸了1.5个车长，向后延伸了3个车长，向下延伸了10倍于车体静吃水深。边界条件：入口边界给定为速度入口；出口和周边外域给定为开放边界；车体给定为移动固壁边界。

图6　边界条件设置Fig.6 Boundary conditions

3　数值计算结果及分析

3.1相同弗劳德数Fr车体航行姿态分析

以车体Fr=0.93的计算过程为例进行分析。如图2，车体初始姿态为0.2 m的重心升沉，8°的纵倾角。不同时刻，升沉和纵倾的变化曲线见图7和图8.可知车体在运动响应计算过程中，经过一段时间的大幅振荡调整，找到了受力平衡点，达到动平衡状态。车体结束大幅振荡调整大概需要3.5 s的时间，初始大幅振荡调整阶段过后逐步微调进入稳定航行姿态，经过初期短暂振荡调整，车体后续航行状态基本稳定。

图7　车体升沉变化图Fig.7 Time-history of the predicted center-of-gravity

图8　车体纵倾变化图Fig.8 Time-history of the predicted pitch angle

为了说明车体航行姿态变化的原因，与此同时，对车体的单位排水量总阻力CD、单位排水量总升力CL和纵倾力矩系数CN也进行监测，从而确定车体达到平衡的时间和平衡时的具体状态。

3个无量纲系数的定义为

式中：D为静水时车体排水重力；N为车体所受实际绕z轴转矩。

3个系数的变化情况见图9.为了下文表述方便，将图中单位排水量总阻力CD、单位排水量总升力CL和纵倾力矩系数CN采用6个时间点进行划分，分别为t1～t6.

对比图7、图8和图9曲线可知，车体的运动分为两个阶段：第1个阶段（t1～t4时刻）为大幅振荡调整阶段，在此阶段中车身的运动幅度较大，车体的升沉变化量最大达到了0.25 m，纵倾变化量最大达到了0.13 rad，动态响应特性明显；第2个阶段（t4～t6时刻）为车体的平稳运动阶段，在此阶段中车身的运动幅度较小，升沉变化量最大只有0.08 m，纵倾变化量最大只有0.05 rad，动态响应特性减弱。以初始姿态运动时，在运动初期，航行中车体所受浮力大于静止时所受浮力，故车体重心会提升直到浮力和重力再度平衡；设计姿态航行时，纵倾力矩并非为0，故车体纵倾角会减小直到达成新的力矩平衡。经过这两个阶段之后，车体已进入平稳运动阶段。在平稳运动阶段中，车体的航行姿态也会略有改变，继而总升力和纵倾力矩也会改变，但变化幅度已不大。由此可以得出车体平衡实际是动态平衡的概念，从而弥补了静态网格技术车体姿态静止不动的不足。

为了进一步分析车体运动规律变化原因，表1中列出了典型时刻车体仰视图的压力变化云图和纵剖面速度矢量图，结合图9纵倾力矩系数的变化，选定这几个时刻分别为力矩最大值时刻t2，力矩最小值时刻t3以及平稳运动阶段中的时刻t5.

图9　不同时刻车体单位排水量总阻力、单位排水量总升力和力矩系数变化趋势图Fig.9 Resistance，lift and torque coefficients at different moments

表1　典型时刻车体压力云图及速度矢量图Tab.1 Contours of static pressures and velocity vectors at different moments

大幅振荡调整时期，纵倾力矩系数的两个峰值点，最大峰值t2及最小峰值t3时刻，对比其压力分布，可发现各处的压力分布与其水上兴波密切相关。在t2时刻的仰视图中，后底部的高压区域对车体产生了一个绕z轴正的转矩，仰视图中其余位置的压力分布几乎沿纵剖面对称；而t3时刻的仰视图中，首前部产生压力较高的区域，将产生一个绕z轴负的转矩。但是在平稳运动时期的时刻t5，车体表面的压力基本趋于稳定，不会导致大的波动出现。

观察纵剖面的速度矢量图，在t2和t3时刻，车体首部、尾部附近有大量旋涡存在，在车尾处水流有明显向上移动的趋势。旋涡的出现会增加涡流损失，增大两栖车在水中航行时的阻力，并且在旋涡区域内部会形成低压区，增大车体的黏压阻力，进而改变车体表面的压力分布情况，从而影响车体航行姿态的调整；而当车处于平稳运动t5时刻，车体首部、尾部附近只有少量兴波存在，不再有大尺度的旋涡出现，车体姿态趋于稳定。

3.2不同Fr数车体航行姿态对比

不同Fr下车体稳定后最终航行绕流形态图见表2。由表2车体绕流形态图得知，随着Fr数的增大车体重心位置逐渐升高，车体逐渐抬起。Fr= 0.93时，车体处于排水航行状态，静浮力是支持车重的主要成分，车体大部分在水面以下，车首兴波较小，尾部流场没有明显的“鸡尾流”形态［18］。Fr= 1.40时车体稳定后重心位置高于Fr=0.93工况，但仍有部分车体在水面以下，即车的重量由静浮力和动升力共同支持。这表明随着航速的提高，动升力的比重将越来越大，尾部流场开始出现“鸡尾流”形态。Fr=1.87时，大部分车体从水中抬起，车体基本处于在水面滑行的状态，车体首部兴波上卷较大，两侧兴波向后喷溅，尾部出现典型的“鸡尾流”形态，动升力为支持车重的主要成分。

3.3数值模拟结果与实验结果对比

实验水池长160 m，宽7 m，水深3.7 m；拖车为空腹梁结构，全数字式直流调速系统，速度范围0.01～8 m/s，速度精度为±1 mm/s；船模采用钢板制成，内外涂防腐漆。车模拖曳方式示意图如图10所示。

表2　不同Fr下车体稳定航行绕流形态图Tab.2 Water volume fractions round vehicle at different speeds

图10　车辆拖曳方式图Fig.10 Amphibious vehicle model in the towing tank

如图10所示，拖点位置为推力轴线和重心纵向位置线的交点；阻力仪用于记录车体航行阻力；首尾升沉记录装置用于记录车体航行时首尾升沉位移，并可由此计算出车体航行纵倾角；首尾导航板用于防止车体偏航。具体实验方法执行我国船舶行业标准CB/Z 244—1988滑行艇船模阻力测试方法

将自动调整方法模拟值与拖模实验值进行对比（车体最终稳定航行绕流形态和航行姿态），见表3和图11.由表3和图11可知，模拟结果与实验结果吻合较好，但是受网格尺度影响，无法描述车首周围向外抛出的喷溅水体与空气混合后的流动结构。

表3　数值模拟与实验的绕流形态图对比Tab.3 Comparison of simulation and experimental results for sailing

图11　升沉和纵倾角的模拟值与实验值对比Fig.11 Comparison of simulation and experimental results of the predicted center-of-gravity and the predicted pitch angle

4　结论

本文通过应用动网格技术，对两栖车静水直航的航行姿态进行了数值模拟，得出以下结论：

1）结合HC算法的动态网格技术可以有效求解两栖车的航行姿态，数值计算结果与实验结果的吻合程度较好。

2）不同Fr数下车体稳定后的航行姿态有显著差异，具体表现在随着Fr数的增大车体重心位置逐渐升高，动升力在支持车重成分中所占的比例越来越大，静浮力则越来越小。

（

）

［1］王涛，徐国英，郭齐胜.两栖车辆水上动态性能数值模拟方法及其应用［M］.北京：国防工业出版社，2009：1-7. WANG Tao，XU Guo-ying，GUO Qi-sheng.Amphibious vehicle numerical simulation of water［M］.Beijing：National Defense Industry Press，2009：1-7.（in Chinese）

［2］万晓伟，姚新民，王涛.基于CFD的两栖车辆水上阻力预测［J］.装甲兵工程学院学报，2013，27（3）：26-30. WAN Xiao-wei，YAO Xin-min，WANG Tao.Water drag prediction of amphibious vehicle based on CFD［J］.Journal of Academy of Armored Force Engineering，2013，27（3）：26-30.（in Chinese）

［3］王涛，郭齐胜，徐国英，等.基于CFD的两栖车辆水上耐波性能研究初探［J］.系统仿真学报，2009，21（10）：3142-3145. WANG Tao，GUO Qi-sheng，XU Guo-ying，et al.Preliminary study on seakeeping capacity of sailing amphibious vehicle based on CFD［J］.Journal of System Simulation，2009，21（10）：3142-3145.（in Chinese）

［4］袁益民，郭齐胜.两栖车辆水上操纵运动建模与仿真［J］.计算机仿真，2011，28（7）：336-339. YUAN Yi-min，GUO Qi-sheng.Modeling and simulation of maneuvering behavior of sailing amphibious vehicle［J］.Computer Simulation，2011，28（7）：336-339.（in Chinese）

［5］韩占忠，王国玉，闫为革.两栖车辆航行黏性阻力数值分析［J］.车辆与动力技术，2003（2）：6-10. HAN Zhan-zhong，WANG Guo-yu，YAN Wei-ge.Numerical simulation of viscosity resistance around a running amphibian vehicle［J］.Vehicle&Power Technology，2003（2）：6-10.（in Chinese）

［6］王涛，徐国英，姚新民，等.两栖车两相绕流场的模拟与水上快速性分析［J］.机械工程学报，2009，44（12）：168-172. WANG Tao，XU Guo-ying，YAO Xin-min，et al.Numerical simulation of two phases flow field around amphibious vehicle and analysis of power performance on water［J］.Journal of Mechanical Engineering，2009，44（12）：168-172.（in Chinese）

［7］吴珂，宋桂霞，赵又群.基于流体仿真的两栖车车轮收起前后阻力对比分析［J］.系统仿真技术，2007，3（4）：187-191.WU Ke，SONG Gui-xia，ZHAO You-qun.Comparative analysis on resistance of amphibious vehicle before and after retracting wheels based on fluid simulation［J］.System Simulation Technology，2007，3（4）：187-191.（in Chinese）

［8］徐国英，王俊，周景涛.基于CFD的两栖车辆阻力和浮态数值模拟［J］.舰船科学技术，2006，28（4）：22-25. XU Guo-ying，WANG Jun，ZHOU Jing-tao.Numerical simulation of the amphibious vehicle's drag force and attitude based on CFD［J］.Ship Science and Technology，2006，28（4）：22-25.（in Chinese）

［9］孙伟，韩建礼，刘西侠.基于动网格模型的两栖车辆数值模拟［J］.舰船科学技术，2009，31（1）：146-150. SUN Wei，HAN Jian-li，LIU Xi-xia.Numerical simulation of amphibious vehicle based on the dynamic-mesh model［J］.Ship Science and Technology，2009，31（1）：146-150.（in Chinese）

［10］万晓伟，王涛，姚新民.两栖车辆水上高速航态的数值仿真方法研究［J］.计算机仿真，2012，29（6）：323-327. WAN Xiao-wei，WANG Tao，YAO Xin-min.Study on numerical simulation method for high speed navigation state of amphibious vehicle［J］.Computer Simulation，2012，29（6）：323-327.（in Chinese）

［11］Launder B E，Spalding D B.The numerical computation of turbulent flows［J］.Computer Methods in Applied Mechanics and Engineering，1974，3（2）：269-289.

［12］叶正寅，张伟伟，史爱明，等.流固耦合力学基础及其应用［M］.哈尔滨：哈尔滨工业大学出版社，2010：37-51. YE Zheng-yin，ZHANG Wei-wei，SHI Ai-ming，et al.Fundamentals of fluid-structure coupling and its application［M］.Harbin：Harbin Institute of Technology Press，2010：37-51.（in Chinese）

［13］Dowell E H，Hall K C.Modeling of fluid-structure interaction［J］.Annual Review of Fluid Mechanics，2001，33：445-490.

［14］Munch C，Ausoni P，Braun O，et al.Fluid-structure coupling for an oscillating hydrofoil［J］.Journal of Fluids and Structures，2010，26（6）：1018-1033.

［15］Akcabay D T，Young Y L.Hydroelastic response and energy harvesting potential of flexible piezoelectric beams in viscous flow［J］.Physics of Fluids（1994-present），2012，24（5）：054106.

［16］Young Y L，Chae E J，Akcabay D T.Hybrid algorithm for modeling of fluid-structure interaction in incompressible，viscous flows［J］.Acta Mechanica Sinica，2012，28（4）：1030-1041.

［17］Young Y L，Chae E J，Akcabay D T.Transient hydroelasticreponse of a flexible hydrofoil in subcavitating and cavitaing flows［C］//29th Symposium on Naval Hydrodynamics.Gothenburg，Sweden：the US Office of Naval Research，2012：26-31.

［18］杨爽，王志东，吴贺贺，等.滑行艇纵向运动响应与静水航行数值模拟分析［J］.船舶工程，2012，34（4）：31-35. YANG Shuang，WANG Zhi-dong，WU He-he，et al.Numerical simulation and analysis on heave-pitch coupled motion response and motion in calm water of planing ship［J］.Ship Engineering，2012，34（4）：31-35.（in Chinese）

Numerical Simulation of Navigating Pose for Amphibious Vehicle based on Dynamic-mesh Model

ZHAO Bin1，ZHANG Min-di1，JU Dong-mei2
（1.School of Mechanical Engineering，Beijing Institute of Technology，Beijing 100081，China；2.Ordnance Science and Research Academy of China，Beijing 100089，China）

A dynamics model for the navigating pose of amphibious vehicle in water is established.The dynamic mesh technology and the hybrid coupled algorithm are used to study the change rule of navigating pose of amphibious vehicle in sailing，and the navigating pose of amphibious vehicles in still water are discussed.The numerical results are consistent with the experimental results.The result shows that the vehicle undergoes a sharp motion and a smooth motion during movement before achieving the ultimate navigating pose.After achieving the ultimate navigating pose，the center-of-gravity and the pitch angle slightly change.During the process of sailing from drainage state to coasting state，the body center of gravity gradually increases；the dynamic buoyancy to support the vehicle weight increases，and the static buoyancy decreases.

ordnance science and technology；amphibious vehicle；navigating pose；dynamic mesh technology；numerical simulation

U674.78

A

1000-1093（2015）03-0412-09

10.3969/j.issn.1000-1093.2015.03.005

2014-06-04

国家自然科学基金项目（51106009）

赵彬（1989—），男，硕士研究生。E-mail：qzb1989q@163.com；张敏弟（1971—），女，副教授，硕士生导师。E-mail：zmd1971@bit.edu.cn