水下航行体充气上浮仿真方法研究

张晓光， 李斌， 党会学， 温金鹏， 孙潘

(1.西北工业大学 航空学院， 陕西 西安 710072； 2.长安大学 建筑工程学院， 陕西 西安 710061；3.中国工程物理研究院 总体工程研究所， 四川 绵阳 621900)

0 引言

水下无人航行体是指能在水下自主航行的智能化装置。航行体的相关装置成本很高，在科研试验或日常操演中，为了使航行体在航行结束后可靠上浮，以便安全打捞、回收，必须设计可靠的助浮装置。充气式上浮装置可以提供较大正浮力，且气囊折叠后所占空间小，质量轻，调试、使用、维修方便，可重复使用，是一种理想的助浮装置。

助浮装置对航行体运动轨迹有较大影响，为给予助浮装置的研制提供参考，需进行航行体6自由度(6DOF)运动的仿真模拟。对水下航行体充气上浮的研究，目前多采用试验方法。吕汝信[1]曾对气囊式充气上浮装置进行过系统研究，包括装置的系统构成、导流罩解脱技术等，并进行了水中拖曳试验、导流罩干扰试验以及实航试验。甄文强等[2]利用试验设计(DOE)方法研究了无人航行器的充气上浮过程，重点研究了气囊充气过程对潜航器运动的影响。李春雨等[3]研究了浮囊整流罩的安全解脱设计及试验验证方法。但试验研究受试验投入、场地和测试设备等条件的影响，成本相对较高，并且水下试验过程的可测数据往往非常有限，难以满足精细分析的需要。随着计算流体动力学以及动网格技术的发展，水下航行体充气上浮动力学仿真的相关文献陆续发表。王栋等[4]通过有限元方法，建立三维造波水池模型，研究了浮囊式操雷在不同强度波浪作用下的漂浮过程，但不能模拟展开上浮过程。程文鑫等[5]、叶慧娟等[6]建立了鱼雷用助浮环型气囊的充气展开动力学模型，研究了不同水深条件下气囊单纯充气展开行为，但未研究气囊体积膨胀对鱼雷航行姿态的影响。

要同时进行水下航行体气囊展开和上浮过程的紧耦合行为分析，需要将气囊展开动力学模型和航行体的水动力学模型同步进行紧耦合求解。实现这一紧耦合求解，一方面要模拟复杂的折叠气囊在充入气流作用下的非规则展开过程，还需应用动网格技术模拟水下航行体的6DOF运动，以实现多学科协同仿真分析。现有公开发表的航行体水动力模型的研究有入水问题[7-8]、潜射导弹出水问题[9-10]、水下航行体的水下航行问题[11-12]等。如潘光等[13]利用有限体积法与非结构动网格技术，运用协调一致的压力耦合方程组半隐式 (SIMPLEC) 方法和k-ε湍流模型(k为湍动能，ε为湍流耗散率)，对无人潜航器进行水下对接过程的数值模拟，得到了在不同攻角下无人水下航行体流体动力随着距离的变化规律。Nair等[14]则研究了轴对称体进出水的计算流体力学(CFD)方法。

通过以上总结发现，目前对水下航行体充气上浮的研究多集中于试验以及气囊的研究，相关仿真研究主要涉及刚体的水动力学研究，同时考虑气囊展开和航行体上浮姿态变化的仿真研究尚未有人开展。但是这一问题对于确保航行体的安全可靠回收不可避免，若气囊体积膨胀带来的增浮效应不够，导致航行体持续下沉，并进而增加气囊外部水压，阻碍气囊的后续展开，则可能导致回收失效。而本文研究的航行体上浮过程，不仅需要仿真航行体的刚体运动，还需要仿真气囊的局部大变形。制约这一问题进行紧耦合分析的难点在于：现有算法无法同时考虑气囊展开内流场和外部水流场的相互作用，以及气囊非规则大幅变形时动态外流场网格的同步生成和多重流体与固体(简称流固)耦合过程的迭代求解。根据对现有仿真软件的调研分析表明，现有软件无法直接实现上述过程的仿真。

本文综合考虑到仿真的可实现性和计算结果的有效性，提出一种两阶段的松耦合求解方法：第1阶段为折叠气囊展开过程的模拟。此时不考虑航行体的运动姿态，仅关注定深水压作用下气囊的折叠与展开，获得助浮气囊的折叠展开过程中的体积膨胀率时程曲线。第2阶段为航行体上浮过程的姿态运动模拟。航行体运动过程主要受气囊展开体积影响，第2阶段仿真以第1阶段所获得每个时间步的体积膨胀为输入，忽略气囊内凹折叠等非规则变形的影响，采用本文提出的动网格技术同步表征第1阶段计算的气囊展开过程，并实时更新水动力学模型的流场网格，采用单向流固耦合与航行力学多学科协同分析方法来实现航行体上浮姿态的分析。

1 充气展开过程的建模分析

1.1 有限元模型的建立

整个水下航行体几何尺寸如图1(a)所示。几何模型如图1(b)所示，图中红色区域表示气囊所在位置。航行体长度为l，直径为d，气囊所在区域宽度为0.072l，气囊所在区域中部距航行体头部0.146l.

图1 水下航行体几何尺寸与几何模型Fig.1 Geometric dimensions and model of underwater vehicle

整个航行体的刚度相对较大，在本文分析航行体充气上浮过程的动力学行为时忽略航行体自身的弹性。实际充气上浮的过程包括：发动机停车—接收上浮指令—气囊充气—抛整流罩—预定膨胀体积达到—排气稳压维持浮力—航行体上浮至水面。引言部分已经指出，完全真实地考虑航行体水下环境，以紧耦合方式来仿真模拟航行体的气囊充气展开及上浮过程是不现实的。本文建模分析的第1阶段是忽略航行体姿态变化对折叠气囊展开过程的影响，仅考虑定水深恒压环境，将气囊展开子系统单独分离出来，并应用控制体积方法来建模分析折叠气囊的充气展开过程。有关折叠气囊的充气展开建模方法在相关文献已经有较为系统的描述[5-6]，在此简要交代本文模型的建立，并给出计算结果，作为第2阶段建模分析的输入。

将展开子系统的各组分模型在三维建模软件CATIA中建立几何模型，划分网格后在有限元前后处理软件LS-PrePost中合并得到展开子系统的有限元模型。对于展开子系统，其有限元模型主要由折叠环形气囊、气瓶、充气管以及水下航行体部分刚壳段4部分组成，如图2所示。由于水下环形气囊为三维气囊，其一般采用较为复杂的环向折叠方式安装在水下航行体上，几何属于不可直接展平气囊，需采用初始矩阵法进行有限元网格模型的建立[15-17]。对于该环形气囊，为得到单元节点编号、连接方式一一对应的参考网格和映射网格，考虑到该环形气囊的参考网格和映射网格均为规则周期结构，本文采用数学分析软件MATLAB自编程序实现两套网格的建立，如图3所示。对于深水压力的影响，通过设定不同的环境压强进行模拟，其计算公式为p=pe+ρgh，pe为一个标准大气压，取101 300 Pa，ρ为水密度，h为气囊所在水下深度，g为重力加速度，g=9.8 m/s2.

图2 展开子系统有限元模型Fig.2 Expanded system finite element model

图3 环形气囊参考网格与映射网格及折叠方式示意图Fig.3 Reference grid, mapping grid and folding mode of annular airbag

1.2 气囊展开过程仿真结果

对气囊在水下50 m深度的充气展开过程进行仿真分析，得到其充气展开过程变形图，如图4所示。随着气瓶充气，囊内压强不断升高，整流罩打开，气囊不断向外展开，直至2.325 s时气囊充气完成，气囊的体积不再发生变化。充气过程中气囊体积膨胀率与时间的关系曲线如图5所示。

图4 气囊充气展开过程变形图Fig.4 Deformation diagram of airbag during deploying

图5 气囊体积膨胀率随时间变化曲线Fig.5 Changing curve of volume expansion rate of airbag

2 上浮过程的仿真分析方法

对于航行体充气上浮运动的仿真问题，主要是要解决气囊局部膨胀运动与航行体刚体运动姿态变化的耦合运动问题，本文提出采用流固耦合和航行力学的多学科协同分析方法来解决这一问题。为避免气囊初始展开阶段的非规则内凹几何带来的网格负体积问题，本文提出采用规则弧线的旋成体来等效模拟每个时间步的气囊体积膨胀行为。具体计算时，先根据图5曲线提供的每个时间步气囊体积，确定当前时间步等效旋成体的子午线，计算气囊处网格点的运动；然后求解Navier-Stokes方程，得到膨胀过程每一时刻的流场解；再沿着航行体表面对压力进行积分求得航行体所受的力(矩)，将已求得的力(矩)代入6DOF刚体运动方程，求解出航行体下一时刻的运动学物理量；利用动网格技术中的弹簧光顺法和局部网格重构法，根据计算出的边界值自动计算出航行体移动后的网格节点位置。重复上述过程，即可对气囊展开、航行体上浮的整个过程进行数值模拟。本文在CFD软件Fluent的基础上，利用用户自定义函数(UDF)进行编程，实现了这一耦合求解。求解时每一时间步计算分4个步骤：1)计算气囊膨胀的等效运动；2)求解非定常流场控制方程；3)求解航行体6DOF刚体运动；4)更新流场动网格。

2.1 气囊体积等效膨胀的动网格计算方法

考虑到动网格实现的可行性，本文将气囊处网格的等效膨胀简化处理为截面为圆弧的形状，气囊膨胀过程的体积时间曲线由1.2节气囊的膨胀过程仿真结果可得。计算气囊上网格点的运动时需要输入圆弧对应的半径，所以需要推导出气囊体积Va(mm3)与截面圆弧半径r(mm)的关系式：

当圆弧为劣弧时，

(1)

当圆弧为优弧时，

(2)

具体计算气囊处网格点的移动时，需要先定义一个结构体数组，里面分别存储编号、网格点的三坐标(未膨胀前气囊处网格点坐标已知)以及网格点在圆弧上的角度比例，再进行旋转角度的计算，确定它在真实坐标系下的坐标。

2.1.1 角度比例的计算

图6为气囊处点坐标运动示意图，其中：C′D′段为气囊所在位置，A、B分别为航行体的前后端点，计算时以C和D点的轴向位置为依据，根据每一个点在CD段上的距离，确定其在角度区域内的比例因子；N′点在CD段上的对应位置为N，其在气囊膨胀过程中，位置不断更新，分别到了N″、N‴、N‴′；E′点在CD段上的对应位置为E，其在气囊膨胀过程中，分别到了E″、E‴、E‴′.

图6 气囊处点坐标运动示意图Fig.6 Schematic diagram of motion of airbag’s point coordinates

2.1.2 气囊处网格点坐标的更新计算

在计算网格点坐标时，需要用到罗德里格旋转公式：设v是一个三维空间向量，k是旋转轴的单位向量，则v在右手螺旋定则意义下绕旋转轴k旋转角度θ得到的，向量可以由3个不共面的向量v、k和k×v构成的标架表示：

vrot=cosθv+(1-cosθ)(v·k)k+sinθk×v.

(3)

使用注意事项：所有的向量(转轴、被旋转向量)，均默认向量的起始点为原点(0 m, 0 m, 0 m)，并且要求旋转向量和被旋转向量交于原点。如果要求真正的坐标，还要加上原来的点坐标。如果被旋转向量和转轴向量不相交，则要把被旋转向量分解成两个向量(被旋转向量起点和转轴向量起点构成一个向量；被旋转向量终点和转轴向量起点构成一个向量)，然后分别计算向量两个点旋转后的终点坐标，再把这两个终点坐标重新构成向量。把二维转角映射到三维上去。在图6所示的二维图中，首先确定的向量为AN′、BN′，然后计算它们矢量积的向量AN′×BN′，再确定移动向量N′N″. 把BN′的B点平移到N′点，可构成图7所示的三角形。

图7 平移BN′得到的三角形Fig.7 Triangle by translating BN′

以此三角形为依据，向量N′N″减去向量BN′就是剩下一边的向量。利用余弦定理，可以确定相应的角度θ，这样就在二维平面条件下确定了夹角θ.

在三维条件下，对于每个网格结点位置，都有对应的AN′、BN′，则利用矢量积AN′×BN′，然后绕着这个单位矢量，把BN′旋转角度θ(在二维中已经计算得到)，则可以得到沿着N′N″的向量方向；确定了矢量长度|BN′|/|N′N″|的值，则可以直接确定N′N″量各个分量值；然后利用ON″=ON′+N′N″，则可以确定N″的点坐标。同样的方法可求出其余点坐标。

2.2 流场求解方法

展开系统上浮运动的模拟是考虑刚体运动与气囊体积变化的非定常流体动力学问题。在模拟计算过程中，需要对空间离散格式、压力速度耦合格式、时间推进格式、湍流模型、考虑气囊体积变化及刚体运动的动网格等进行相应的设置和调整，以实现大范围、大变形动态运动过程中，能准确模拟气囊展开、水下航行体的上浮过程。

对于存在边界运动的控制体，一般标量φ守恒方程的积分形式，在任意控制体积V上，其边界移动都可以写成：

(4)

式中：u为流速矢量；ug为动网格的网格速度；Γ为扩散系数；Sφ为φ的源项；∂V表示控制体积V的边界；A是面向量。

在方程中出现的时间导数用2阶时间精度向后差分表示为

(5)

式中：n+1、n、n-1分别表示第n+1个时刻、第n个时刻、第n-1个时刻。控制体积由(6)式进行迭代计算：

(6)

为满足守恒方程(4)式，体积对时间的导数由(7)式进行计算：

(7)

以上对积分形式的Navier-Stokes方程，基于有限体积方法用2阶迎风格式进行了空间离散；然后采用2阶欧拉向后隐式格式进行时间推进；最后采用Transition SST湍流模型来封闭Navier-Stokes方程。

2.3 水下航行体受力分析

水下航行体在完成任务后，调整航行的速度、姿态角，发动机停车，气囊展开产生浮力带动航行体上浮。上浮过程中其受力情况如图8所示。图8中R为航行体和气囊受到的阻力，Ff为航行体与气囊产生的浮力，mg为航行体自身重力。计算公式分别为

(8)

式中：D、Df分别为航行体与气囊的阻力系数；S、Sf分别为航行体与气囊的浸湿面积；v为水下航行体的速度；Cx为航行体的x轴方向力系数；Cy为航行体y轴方向的力系数；Z为航行体运动过程受到的侧向力；Cz为航行体z轴方向的力系数；Vtot为航行体与气囊的体积。

图8 水下航行体受力示意图Fig.8 Force diagram of underwater vehicle in the process of floating movement

气囊展开过程中会产生逐渐增大的正浮力，使得水下航行体所受到的合力不断变化，其运动姿态也相应地随时间改变。现有CFD软件中的6DOF模型主要用于模拟计算域中的刚体在流体力作用下的轨迹及姿态，相关计算均建立在物体边界为刚性边界的前提条件下。因此，为了模拟水下航行体由于气囊展开带动的上浮过程，需用自编UDF对航行体壁面积分，计算流体作用及重力作用共同产生的合力(对质心的3个分力)和合力矩(绕体轴的3个分力矩)，再代入到6DOF刚体动力学方程中，计算出航行体的运动轨迹及姿态。本文动网格设置与刚体6DOF运动程序经美国阿诺德工程发展中心(AEDC)[18]标准弹翼分离算例验证，准确有效。整个仿真计算的流程图如图9所示。

图9 航行体充气上浮仿真分析流程图Fig.9 Flow chart of simulation of underwater vehicle’s floating

3 水下航行体充气上浮仿真计算

3.1 几何模型和网格生成

水下航行体几何模型如图1所示，航行体模型为壳体，航行体质心距头部距离为0.466l，质量为245 kg，航行体浮心距头部0.464l. 气囊充气完成后，航行体与气囊的总体积为0.294 m3，浮心距航行体头部距离为0.407l.

在网格划分方面，航行体采用结构化网格划分，而在航行体与外流场之间又划分了一块椭球形区域，在航行体与椭球形之间的区域采用结构化六面体网格划分，来进行气囊展开的计算，且椭球域与航行体是绑定在一起的，所以航行体上浮时，椭球域与航行体做相同运动。对于航行体运动所需的空间区域，则采用四面体非结构化网格进行剖分，以充分发挥四面体网格在动网格方面的优势；在六面体和四面体网格之间，采用五面体(金字塔)网格来实现过渡，如图10所示。

图10 计算域网格示意图Fig.10 Schematic diagram of computational domain grid

具体进行网格划分时，计算域尺寸的选取相当重要，椭球域过大会导致外流场区域过大，椭球域过小会使气囊处网格膨胀时挤压六面体网格导致网格出现负体积。外流场区域过大会导致计算量成倍增加，外流场区域太小会导致计算出现“回流”现象，从而导致计算不收敛。经过参考文献和计算验证，椭球长度方向的直径定为10l，宽度和高度方向的直径定为5l. 外流场的尺寸定为100l×29l×57l，如图11所示。整个模型的网格单元数大约为600万。

图11 计算域外流场示意图Fig.11 External flow field of computational domain

本文算例航行体航行水深为50 m，发动机停车时，航行体前行速度为5 m/s. 进行瞬态计算时，在保证计算收敛和计算精度的前提下，为提高计算效率，需加大时间步长以减少计算时间，本文通过多次收敛性检验计算后，设定计算时间步长为1×10-3s. 在计算气囊膨胀的过程中，即使对气囊进行了规则化处理，气囊膨胀时对周围的结构化六面体网格挤压还是很严重，在气囊体积膨胀至52%时，六面体网格由于扭曲率过大而出现负体积。所以本文在出现负体积的时间步重新划分图10(c)椭球域内的六面体网格，进行重启动计算。因为椭球域内的六面体网格主要进行气囊膨胀的计算，所以重新划分此区域内的网格并不会影响此后的计算结果。

3.2 计算结果分析

3.2.1 航行体运动轨迹与姿态分析

图12为航行体上浮过程的网格截面示意图，其直观展示了航行体充气上浮过程中网格的变化过程，椭球域与航行体绑定在一起，它们具有相同的运动轨迹与姿态。椭球域内的结构化六面体网格用于气囊膨胀的计算，随着航行体的转动与上升，椭球域周围的五面体网格及外流场的四面体网格也逐步变化。当网格变形严重扭曲率过大时，采用网格光顺与局部网格重构的方式重新提高网格质量。

仿真获得的气囊膨胀和航行体上浮过程如图13所示，图14、图15分别为航行体质心运动的平面轨迹投影以及三维轨迹。计算结果表明，在当前工况条件下，随着气囊的不断膨胀，航行体所受浮力逐渐增大，气囊带动航行体由水平位置逐渐转动到竖直位置，气囊充气展开提供的增浮作用可有效实现航行体的上浮回收。

图13 航行体充气上浮过程Fig.13 Aeration and floating process of vehicle

图14 航行体质心位移随时间的变化曲线Fig.14 Displacement of center of mass of vehicle over time

图15 航行体质心在三维空间中的运动轨迹Fig.15 Trajectory of center of mass of vehicle

图16 航行体质心线速度随时间的变化曲线Fig.16 Linear velocity of center of mass of vehicle over time

由图14x轴方向质心位移曲线可知，在大约2.3 s以前，航行体一直沿x轴正向方向运动，而2.3 s后航行体又向x轴负向方向运动。航行体质心线速度如图16所示，由图可知，这是因为航行体x轴方向的速度逐渐减小，并且在2.5 s后，x轴方向的速度由0 m/s逐渐向负方向增加。航行体x轴方向的速度之所以会向负方向增加，是由于航行体上浮时受到水动力的作用，此时航行体又具有较大俯仰角，导致航行体出现“回旋”，航行体向后运动。由航行体质心在水深方向(y轴方向)的位移时间曲线可知，在前1.1 s由于航行体重力大于所受浮力，航行体一直处于下降趋势，而后随着气囊膨胀使航行体浮力增大，航行体逐渐上升。航行体开始上升时的气囊体积膨胀率为62%，由于气囊膨胀速度相对较快，所以航行体浮力增加很快，航行体在水中并未下降很深便开始上浮。航行体上浮过程中已向z轴负方向偏移了0.4 m，这是由于航行体受到沿z轴负方向力的作用。

图17为航行体运动过程中欧拉角随时间的变化曲线，从中可以看出偏航角在前1.5 s几乎没有什么变化。偏航角在1.5 s后逐渐向负向增加，而后又逐渐沿正向增加。说明航行体先向左偏转，向左偏转了约19°后又向右偏转了20°. 横滚角开始时变化不大，在3 s后上浮过程才出现较大滚转。虽然开始时浮力还很小，航行体还处于下降趋势，但由于气囊膨胀，浮心在质心前方不断向航行体头部方向移动，浮力产生了抬头力矩，航行体从开始时俯仰角就一直在增大，直至增大到90°. 而后由于水动力的作用，航行体还会接着向后俯仰，最终达到106°后又逐渐减小。结合偏航角的变化，可以看出航行体上浮时处于螺旋上升过程。

图17 航行体欧拉角随时间的变化曲线Fig.17 Euler angle of vehicle over time

3.2.2 航行体合力分量的分析与讨论

图18 航行体在3方向上的力系数Fig.18 Force coefficients of vehicle in three directions

图19 初始倾角-5°时航行体上浮过程流线图Fig.19 Streamline diagrams of floating process of vehicle at an initial inclination of -5°

如图18所示，分别为航行体上浮时合力在地轴系上的分力的无量纲值。图18(a)为航行体合力在x轴方向分量的变化曲线，开始时航行体所受x轴方向力急剧增大，大约0.5 s时，力增大到最大，随后逐渐减小。因为开始时航行体具有“停车”速度，并且开始时气囊膨胀较快，导致航行体所受阻力急剧增大。而后航行体逐渐减速，并且气囊膨胀速度减慢，这使得航行体所受x轴方向力逐渐减小。图18(b)为航行体的y轴方向力变化曲线，航行体在充气上浮过程中，所受的y轴方向力先增大后减小。在前2.3 s，y轴方向力一直增大，并且在1.1 s以前，y轴方向力一直为负。这是因为开始时，航行体所受重力大于浮力，航行体在竖直方向的合力为负。随着气囊不断膨胀，航行体俯仰角逐渐增大，1.1 s后航行体y轴方向力沿正向逐渐增大。结合图11，航行体在1.1 s后便逐渐上浮，1.1 s时总浮力小于重力，而此时攻角为正，说明航行体上浮时受到水动力提供的竖向力作用。大约2.3 s后，y轴方向力逐渐减小。结合气囊体积膨胀率曲线，在2.3 s时，气囊充气完成，航行体所受浮力不变。而此时航行体在y轴方向的阻力随着航行体速度的增大而增大，所以航行体在y轴方向的合力逐渐减小。图18(c)为航行体z轴方向力的变化曲线，航行体开始运动时z轴方向几乎不受力，而后在上浮过程中逐渐受到一定侧向力作用，这也是为什么航行体上浮过程并不是沿Oxy面完全对称的原因。但相较于x轴方向、y轴方向分力来说，z轴方向分力的大小相对较小，所以航行体上浮过程中在z轴方向的位移较小，如图15所示。

3.2.3 不同时刻流场分析

航行体充气上浮过程所受总力可分为体力与流体力，体力包括航行体自身时不变重力和由气囊体积确定的时变浮力，流体力是由航行体的运动与气囊的充气引起的。为了深入了解流体力变化的机理，绘制了航行体在不同初始俯仰角下充气上浮过程中的空间流线变化图，如图19、图20所示。

图20 初始倾角5°时航行体上浮过程流线图Fig.20 Streamline diagrams of floating process of vehicle at an initial inclination of 5°

从图19可知：t=0 s时，航行体以-5°入射角移动，由于总力向下，航行体将向下运动；在t=0.075 s时，俯仰角速度很小，因此流线在z=0 m平面上是对称的，在t=0.20 s时，充气气囊“源”产生较大的俯仰力矩，但入射角相对较小，虚拟气囊结构的外倾角不足以产生流体分离，因此流体形态在z=0 m平面上也是对称的，在t=0.40 s时，外倾角增大，航行体背风面产生对称分离涡对，在这3个时间点上，俯仰角速度没有产生尾翼涡脱落，因此流型由平动决定。在t=0.65 s时，由于足够大的俯仰角速度，非对称分离涡从尾翼脱落，这种加速的俯仰角速度抑制了气囊背风面的分离涡；在t=1.00 s时，流型相似，但流动的非对称性引入了横滚与偏航角速度，从而在航行体头部和气囊上观察到流动的非对称性。气囊背风面涡对的抑制作用可以维持到俯仰角速度明显减速为止；在t=1.40 s时，由于俯仰运动引起的水阻尼，俯仰角速度明显减小，因此附着的涡对与气囊表面分离，形成马蹄涡[19]。这种马蹄涡对与尾翼脱落的非对称涡相互作用，当它们顺流而下时，马蹄涡对消失，在气囊背风侧只能看到两个分离的涡。由于航行体尾翼涡的非对称性，背风涡在运动过程中(t=2.00 s、t=2.35 s和t=3.00 s)保持非对称，并持续使航行体产生横滚与偏航运动。

从图20初始倾角5°时流线的变化过程可以看出，流型的演化过程与-5°的情况类似，不同之处在于：与-5°情况相比，5°时航行体向上运动导致俯仰运动的角加速度增大，从而导致水流分离不对称提前。

3.2.4 试验验证

根据本文仿真计算的工况，开展了一系列试验以验证本文的仿真结果。气囊在陆地上的展开过程如图21所示，展开过程与仿真结果基本相同，气囊均匀向周边展开。进行水下试验时利用电磁计时器对气囊展开过程进行计时，研究了气囊在不同深度下的展开过程，仿真与试验展开时间对比如表1所示，从中可以看出误差均在6 %以内，说明仿真结果的正确性。

图21 气囊陆地试验展开过程Fig.21 Deployment process of airbag on the land

进行航行体上浮试验的测试时，由于测试设备的限制，只能利用电磁计时器测量航行体从气囊开始充气至上浮完成的时间，结果如表1所示，从中可以看出误差均在7 %以内，可验证仿真的正确性。

表1 气囊展开与航行体上浮时间仿真与试验对比

4 结论

本文主要进行了水下航行体利用气囊充气上浮的仿真方法研究，先由LS-DYNA软件利用控制体积算法，模拟气囊在水下的展开过程，得到气囊展开过程的体积时间曲线。考虑到CFD计算模型的需要，保持气囊的体积膨胀率同步等效，用规则的环型旋成体等效模拟气囊的展开过程。利用动网格技术和自编UDF程序实现气囊膨胀过程的表征，然后基于流体力学控制方程与刚体6DOF运动方程求解航行体的运动过程。得出以下主要结论：

1) 在本文的计算初始条件下，气囊2.3 s左右就充气完成，航行体在1.1 s后便开始上浮，在2 s时航行体俯仰角就达到90°，此后航行体便处于螺旋上升过程。

2) 气囊膨胀时，航行体所受x轴方向阻力会急剧增大，而后逐渐减小。在上浮过程中，航行体会因为漩涡结构的不对称性而受到一定侧向力的作用，上浮时并不是沿Oxy平面的对称上浮运动，所以若用二维方法模拟航行体充气上浮运动会不准确。

3) 航行体上浮时会受到水流提供的竖向力作用，所以上浮前应尽量调整航行体的攻角为正，以便充分利用航行体的残余速度，从而加快上浮。

该方法能较好地模拟航行体充气上浮的三维运动过程，为水下航行体上浮的工程设计提供参考，同时本文提出的刚体运动与局部大变形耦合的仿真分析方法也为相似问题的解决提供了参考。