基于Star-CCM+的滑翔器入水动力学特性研究
摘 要：为了研究滑翔器入水动力学特性的影响，文中基于 STAR-CCM+，并结合网格重构技术、重叠网格技术、 VOF 界面捕捉技术等建立了三维水下滑翔器入水数值计算模型，模拟了滑翔器入水全过程特性，并分析了其位移、 速度、流场等物理量的影响。模拟结果显示：所建立的数值模型能较精确地模拟滑翔器入水动力学特性问题，为滑 翔器入水问题提供技术支持。
关键词：STAR-CCM+；滑翔器；入水动力学
一、前言
结构物入水动力学问题广泛存在，入水过程时间历程短、 压力峰值大，易引起结构失效或破坏，其入水过程值得关注， 因此对结构物入水动力学特性进行研究具有十分重要的意 义。国内外学者围绕结构物入水问题开展了广泛研究。吴巧 瑞等[1]采用重叠网格法，对楔形体的入水砰击问题进行模拟。 随着入水深度的增加，气体逸出，出现二次峰值；当斜升角 大于等于 30°时，最大砰击压力出现在楔形体底端附近。速 度对二次砰击峰值的影响要比斜升角的影响大。吴云贵等[2] 对结构物垂直入水过程进行了仿真和试验研究。马林等[3]采 用两相 SPH 模型模拟三维结构物的入水过程。进行三维圆球 的垂直入水研究，给出入水过程中空泡演化过程以及空泡封 闭以后空泡的形态变化和泡内的压强变化，并实现了圆球斜 入水的模拟。刘双等[4]基于 STAR-CCM+和 Fluent，并结合 网格重构技术、重叠网格技术、VOF 界面捕捉技术等建立了 2 种可自相互验证的数值计算模型，对比分析了二维楔形体 的自由入水过程。模拟了不同参数（入水速度和底升角）下 楔形体的自由入水过程，并分析了其对压力、速度、流场等 物理量的影响。Von Karman 等人[5]对水上飞机浮箱降落过 程中撞击水面所受的砰击力进行了理论研究。
滑翔器是一种新型的水下机器人。由于工作效率高、续 航力大、可重复使用等特点被广泛的应用于军事及民用领域， 其投放投放过程容易受到入水速度、下落姿态角等参数的影 响，因此开展滑翔器入水动力学特性研究具有十分重要的意 义。文中基于 STAR-CCM+，并结合网格重构技术、重叠网 格技术、VOF 界面捕捉技术等建立了三维水下滑翔器入水数 值计算模型，模拟滑翔器入水过程，并分析了其位移、速度、 流场等物理量的影响。
二、数值计算模型
对于不可压缩黏性流体，采用连续性方程和雷诺平均方 程建立控制方程组。连续性方程和雷诺平均方程[6，7]分别为
式中：ρ 代表流体密度；t 为时间，u 、v 、w 分别表 示速度沿 XYZ 坐标轴 3 个方向上的分量：
式中：Fbx 、Fby 、Fbz 分别是质量力在 3 个方向上的分量； yx p 为流体内应力张量的分量。
滑翔器入水过程周围的流体会发生介质突变，使得流场 造成强烈的扰动，这种扰动在入水问题的数值求解过程中成 为传统网格法处理入水问题的难点。重叠网格技术可以将计 算与离散为具有任意方式相互重叠的多个不同网格的计算 域，能够高效耦合求解滑翔器六自由度运动方程，且允许各 个网格区域互相重叠、嵌套，在各个子网格间通过插值来传 递流畅信息。
三、几何模型和收敛性分析
1．滑翔器几何模型
如图 1 所示，文中的计算模型 spray 线型带双翼滑翔器， 翼型的选用 NACA 型号，以滑翔器重心为坐标轴原点。具体 参数如表 1 所示。
2．计算区域
如图 2 所示，计算区域的固定坐标系原点位于初始时刻 自由液面上，X 轴指向水平向右为正，Y 轴指向滑翔器左翼 为正，Z 轴垂直向上为正。滑翔器重心在距离水面 h =1m 的 高度处自由下落。定义滑翔器体长为 L，背景区域大小长度 方向的左边界距离计算区域原点 2L，右边界距离计算区域原 点 8L；宽度方向的边界距离计算区域原点 1.4L；顶部边界 距离计算区域原点 2L。
3．几何模型的网格划分与数值验证
当滑翔器落入水中时，在入水初期会激起自由液面（水 介质）的剧烈运动。因此，对自由液面处的网格进行加密， 运动轨迹方向以外的网格逐层递减，平稳地过渡到大的背景 网格。同时，对重叠区域做落体运动的轨迹进行加密，其尺 寸与重叠网格区域的尺寸相同，以确保重叠交界面插值的准 确性。如图 3 为网格划分示意图。
由于滑翔器入水相关研究较少，文中采用与本研究相似 的救生艇入水案例[9]开展收敛性与稳定性分析，具体参数详 见[8]，对案例的总网格数分别设置为 120、160、200 万， 如图 4（a）中所示，为以上网格总数下救生艇 X 轴速度-时 间的仿真结果与文献模拟结果图，从图中可以发现数值计算 中网格收敛性较好，权衡计算的精确性和需要的运算时间， 如图 4（b）为网格数量为 120 万时，采用时间步长分别为 0.0025、0.005、0.01s 下救生艇 Y 轴角速度模拟对比结果， 从图中可以发现数值计算中时间稳定性较好，因此本研究采 用时间步长为 0.01s 作为本案例时间步长。
四、计算结果与讨论
1．入水全过程特性研究
如图 5 为滑翔器以姿态角为 25°，2m/s 的初速度平抛 入水全过程运动变化曲线，可以从图 5（a）-（c）看出在 0.00~0.28s 时间段为自由抛落阶段，0.28~0.72s 时间段是 入水阶段，此时刻最大入水阶段为 0.309m，0.72~1.27s 时间段为第一段上浮阶段，此后滑翔器由于自身重力及惯性 矩，在 1.28s 重新进入下潜状态，达到平衡后，在 2.2s 进 行第二阶段上浮，直至 2.8s 达到稳定状态。从图 5（d）中 可以看出，在初始抛落时，姿态角略有减少；当入水阶段姿 态角受液面对滑翔器尾首部作用而使得滑翔器快速旋转，在 1.13s 内姿态角由 25°减少至-40°；当下潜到最深深度时， 上浮阶段姿态角从最小值-16.89°因水面的反作用力过大导 致滑翔器快速旋转至前端翘起。2s 时滑翔器翘起后头部因重 力向下快速坠落，第二次拍打水面。如图 6 为入水和上浮过 程自由液面变化情况，红色代表水域，蓝色代表空气域，从 图中可以清晰反映滑翔器入水和上浮的动态过程。
2．入水全过程速度变化
为了更清晰的了解滑翔器入水全过程动态特性，滑翔器 入水全过程速度变化曲线如图 7 所示，图中可以看出滑翔器 的自由抛落阶段是带有水平初速度的自由抛落运。自由抛落 阶段，滑翔器在重力加速度的作用下向下加速至-3.58m/s； 入水后，在滑翔器接触自由液面瞬间受到液面的反作用力， 在 0.71s 时垂向速度降为 0m/s；随后在上浮阶段受液面对 滑 翔器 的浮 力作 用， 使得 垂向 速度由 0m/s 加速 至 0.94m/s；随着上浮阶段结束后，滑翔器因自由液面的作用， 浮力减小与重力达到平衡，最终垂向速度变为 0m/s。
3．入水全过程受力变化
从上可知滑翔器在极短的时间内垂向速度发生了剧烈的 变化，因此探究滑翔器入水砰击过程中的受力情况，对于其 结构设计至关重要。如图 8 为入水砰击监测点布置图，在指 定位置选取 P1-P5 等五个不同位置的点作为监测点，其中 P1 点为艇体中性面上首先入水的点。图 9 所示为滑翔器入水 时的砰击压力时历曲线，随着滑翔器抛落入水，底板对液面 的挤压力不断增大，砰击压力也不断增大。待到液面抬升到 一定程度后出现射流现象，此时射流根部出现砰击压力峰值。 之后，射流与物面发生分离，砰击压力迅速减小。最大峰值 可达 0.11MPa。
五、结论
文中基于 STAR-CCM+，并结合网格重构技术、重叠网 格技术、VOF 界面捕捉技术等建立了三维水下滑翔器入水数 值计算模型，模拟滑翔器入水全过程，并分析了其位移、速 度、流场等物理量的影响。模拟结果显示：所建立的数值模 型能较精确地模拟滑翔器入水动力学特性问题，为滑翔器入 水问题提供技术支持。
参考文献
[1] 吴巧瑞，陈明辉，张珍等.基于重叠网格法的结构物入水 砰击研究[J].中国造船，2022，63（1）：102-112.
[2] 吴云贵，陈怀海，贺旭东.结构物垂直入水的冲击特性[J]. 振动.测试与诊断，2013，33（S2）：144-147+225.
[3] 马林.结构物入水的数值和实验研究[D].上海：上海交通 大学，2017.
[4] 刘双，武雨嫣，何广华等.物体自由入水的多参数影响分 析[J].哈尔滨工程大学学报，2020，41（1）：95-102.
[5] CARCATERRA A，CIAPPI E．Hydrodynamic Shock of Elastic Structures Impacting on the Water：Theory and Experiments[J].Journal of Sound and Vibration，2004，271 （1-2）：411-439.
[6] MENTEＲ F Ｒ．Zonal Two Equation Turbulence Models for Aerodynamic Flows[R]AIAA93－2906，1993.
[7] MENTEＲ F，SCHUTZE J，KUＲBATSKII K，et al ． Scale-resolving Simulation Techniques in Industrial CFD[C]//6th AIAA Theroretical Fluid Mechanics Conference，2011.
[8] Lyu，H. G.，Sun，P. N.，Miao，J. M.，& Zhang，A. M.3D multi-resolution SPH modeling of the water entry dynamics of free-fall lifeboats[J].Ocean Engineering，2022，257， 111648.
文章来源：中国水运
