CN109948301A - 基于网格控制的近水面滑行跳跃流固耦合数值预测方法 - Google Patents

Abstract

本发明涉及一种基于网格控制的近水面滑行跳跃流固耦合数值预测方法，属于船舶、近水面航行器技术领域。本发明通过三维流域模型及网格划分、建立计算流体力学模型、建立网格控制方法，进行近水面滑行跳跃流场计算，然后对计算结果进行后处理，获得近水面跳跃随时间变化的动态信息。充分考虑了两相流和结构运动的影响，提高了数值计算结果的可信度；且本发明能够实现对近水面弹跳现象进行高精度的数值预测，实现触水物体不同俯仰角和攻角的设计工况，节省实验成本和时间。

Description

基于网格控制的近水面滑行跳跃流固耦合数值预测方法

技术领域

本发明涉及一种基于网格控制的近水面滑行跳跃流固耦合数值预测方法，属于船舶、近水面航行器技术领域。

背景技术

近年来，随着科学技术应用的迅猛发展，近水面滑行跳跃技术相比于水上飞行和水下航行技术具有较强的灵活性和突防能力。相比于跨介质入水实验和水上航行实验，近水面滑行跳跃的实验难度大、成本高，致使其实验的成功率较低。20世界90年代以来，随着计算机设备的发展和计算技术的进步，逐渐由传统的实验研究，转向数值与实验相结合的研究模式。目前针对近水面跨介质运动的数值研究主要聚焦于物体的入水问题，比如小球入水、射弹入水等，有相对成熟和完善的数值计算模型。而对于近水面滑行跳跃的数值计算研究极其匮乏，以至于无法对触水物体周围复杂的流场环境做出进一步的分析和研究，同时也无法对触水物体的水动力性能和结构相应特性做出合理的解释。因此，对于近水面触水跳跃问题，有必要发展和完善考虑复杂流场和入水物体结构特性的数值预测方法，建立一种基于网格控制的近水面滑行跳跃流固耦合数值预测方法是全新的技术问题，更具实际工程价值和科学意义。

发明内容

本发明的目的是为了解决近水面滑行跳跃的实验难度大、成本高的问题，提出一种基于网格控制的近水面滑行跳跃流固耦合数值预测方法。该方法能够对实验工况的成功性进行评估预测，提高实验的成功率、降低实验危险系数和实验成本。

本发明的目的是通过下述技术方案实现的：

基于网格控制的近水面滑行跳跃流固耦合数值预测方法，具体步骤如下：

步骤一：建立三维流域模型。

根据确定的触水物体的几何形状，确认其中心点为整体流域的坐标原点。然后根据坐标原点依次划分三维流域：内域、中域和外域三部分。内域是将触水物体包围起来的圆形区域，用于触水物体周围网格的加密，同时圆形内域能够方便快捷的调整触水物体的俯仰角度。中域是内域与外域的交界区域，作用于内外域网格之间的过渡，保证外域网格的精确度。外域是将所有域包起来的长方体区域，且三个流域的宽度均大于触水物体宽度。

步骤二：划分三维流域网格。

对步骤一建立的三维流域进行网格划分，其中内域采用非结构化网格对触水物体周围的流域网格进行加密处理，以便捕捉非定常流动细节；中域采用结构化网格，将中域网格由外域往内域方向逐渐加密；外域模型中将中域的位置尺寸作为参数使用，并进行非结构化网格划分。

步骤三：建立计算流体力学模型。

为了能够对步骤二所建立的三维流域网格进行流场计算求解，需要先建立计算流体力学模型。计算流体力学模型包括气液界面捕捉模型、流场控制方程、湍流模型和刚体运动方程。

气液界面捕捉模型：

式中，u_t为速度对时间的一阶导，u＝(u,v,w)是流体速度；F代表体积力；ρ＝ρ(x,t)是流体密度；p为流体压力；μ＝μ(x,t)是介质粘度；D是粘性应力张量；式(1)中的最后一项代表集中在相界面上的表面张力，σ代表表面张力系数；κ代表相界面的曲率；δ是DiracDelta函数；d表示计算区域中的点与相界面的垂直距离；n表示相界面上法向朝外的单位向量

流场控制方程包括质量方程(3)和动量方程(4)：

式中，表示函数对相应的变量求偏导数，ρ为流体的密度，t为时间，u_i、u_j代表流体的速度分量，x_i、x_j代表流体的位置分量，p为流场入口处压强，μ_l和μ_t分别为流体的层流和紊流粘性系数。

湍流模型为SST(Shear-Stress Transport)k-ω湍流模型，包括湍流动力粘性系数μ_t(5)、湍动能k(6)和比耗散率ω(7)：

式中，为ρ为流体的密度，α₁为经验常数，S为剪切力张量的常数项；x_j为坐标向量；u_i为速度分量；Γ_k、Γ_ω分别代表k和ω的有效耗散系数；G_k为由于平均速度梯度而导致的湍流动能k产生项；G_ω为ω产生项；Y_k、Y_ω分别表示k和ω的湍流耗散项；D_ω为交叉耗散项；S_k、S_ω为源项(用户自定义)。

刚体运动方程为包括平动方程(8)和旋转方程(9)：

M＝Iα+ω×Iω＝M_fluid+M_ext (9)

式中，m为刚体质量，为刚体的质心加速度，F为所有外力的矢量和，F_fluid为流体作用力，g为重力加速度，α为角加速度，ω为加速度，M_fluid为流体力矩，M_ext为其他外力矩，(x_G,y_G,z_G)为刚体质心坐标。

步骤四：进行近水面滑行跳跃非定常流场的流固耦合数值计算。

采用步骤三建立的流体力学模型对三维流域近行非定常流场数值计算。在计算流体动力学(CFD)求解器中，对计算参数进行初始化；所述的对计算参数进行初始化的方法为：外域的流场入口、出口及前后边界给定开放条件，上下边界为无滑移、光滑壁面边界条件；外域、中域和内域的接触面均为交界面；内域中触水物体的壁面给定为流固耦合刚体模型，并建立刚体相对于外域x和y方向位移D_x和D_y的表达式和网格质量判断公式。

基于上述边界条件和初始条件，进行近水面滑行跳跃非定常流场的流固耦合数值计算的方法为：

步骤4.1，对流固耦合刚体进行计算，对刚体运动方程的平动方程(8)和旋转方程(9)进行离散求解，计算出流固耦合刚体的运动位移；

步骤4.2，把求得的刚体的运动位移作为流场网格的变形，对气液界面捕捉模型(方程(1)和方程(2))、流场控制方程(方程(3)和方程(4))和湍流模型(方程(5)、方程(6)和方程(7))进行离散求解，得出外域网格变形数据，即刚体在x和y方向位移；

步骤4.3，根据步骤4.2得到的外域网格变形数据和步骤四的网格质量判断公式，基于计算流体力学(CFD)求解器进行网格控制，即对外域网格进行更新；

步骤4.4，用步骤4.3更新后的外域网格替换旧网格，进行流场的非定常计算。

步骤4.5，重复步骤4.1至步骤4.4，直至达到预定求解时间或者数值计算结果收敛，得到近水面滑行跳跃的随时间演化的过程，即实现了对近水面滑行跳跃过程的预测。

步骤4.3所述的网格控制的方法，具体步骤如下：

步骤1)，将步骤4.2监测的当前时间步的刚体在x和y方向位移的几何参数输出；

步骤2)，将步骤1)输出的几何参数作为外域新的输入几何参数；

步骤3)，对步骤2)读取的当前时间步的几何参数进行网格重构，并生成新的外域几何网格；

获得近水面跳跃流场和触水物体随时间的演化信息。

对步骤四所得到的数值计算结果进行后处理，获得近水面跳跃随时间变化的动态信息。后处理方法为：提取流场区域内的流动参数，以及触水物体的位移、速度和加速度，则能够反应触水物体的运动轨迹和过载特性。所述流动参数包括速度、压力和水气界面，其中速度分布通过流线图表示，压力分布通过云图表示，水气界面通过等值面或云图表示。

步骤一所述触水物体的底面为弧形物体。

有益效果

1、本发明的一种基于网格控制的近水面滑行跳跃流固耦合数值预测方法，相比现有技术，不仅可以获得大量触水物体和流场动态信息，更重要的是可以实现触水物体不同俯仰角和攻角的设计工况，能够降低实验的危险系数，节省实验成本和时间；

2、本发明的一种基于网格控制的近水面滑行跳跃流固耦合数值预测方法，充分考虑了两相流和结构运动的影响，提高数值计算结果的精度和可信度，能够对实验的成功性进行评估预测；

3、本发明的一种基于网格控制的近水面滑行跳跃流固耦合数值预测方法，能够实现对近水面弹跳现象进行高精度的数值预测。

附图说明

图1是本发明一种基于网格控制的近水面滑行跳跃流固耦合数值预测方法流程图；

图2是本发明三维流域网格划分的简化示意图；

图3是本发明所采用的网格控制方法具体实现流程图；

图4是使用本发明获得的橄榄球形刚体的运动轨迹。

具体实施方式

下面结合附图，对本发明的具体实施方式进行详细说明。

实施例1

本实施例中选用的触水物体形状为椭圆形状(简称刚体)，刚体底面位于距离静水面80mm的位置处，俯仰角为0°，攻角为2°，入水初速度为100m/s。

一种基于网格控制的近水面滑行跳跃流固耦合数值预测方法的具体流程，如图1所示，主要由以下步骤实现：

步骤一：三维流域模型建立。

针对刚体的几何形状，确认其中心点为整体流域的坐标原点，并确定流域几何参数：内域是直径为2倍刚体长轴的圆形；中域是边长为3倍刚体长轴、内部为圆形的正方形；外域是将所有域包围起来的长方体区域，该长方体区域左端距刚体中心1.5倍长轴，长方体区域右端距刚体中心8.5倍长轴，长方体区域上端距刚体中心1.5倍长轴，长方体区域下端距刚体中心2.5倍长轴。三个流域的宽度均为刚体宽度的2倍，。

步骤二：三维流域网格划分。

对步骤一建立的三维流域进行网格划分，其中内域采用非结构化网格对刚体周围的流域网格进行加密处理，以便捕捉流动细节；中域采用结构化网格，利用节点伸缩设置中域网格由外域往内域方向逐渐加密；外域基于ANSYS Workbench平台调用geometry导入模型并将中域相对于外域的位置尺寸D_x和D_y作为参数使用，随后进行非结构化网格划分，如图2所示。

步骤三：建立计算流体力学模型。

为了能够对步骤二所建立的三维流域网格进行流场计算求解，需要先建立计算流体力学模型。计算流体力学模型包括气液界面捕捉模型、流场控制方程、湍流模型和刚体运动方程。

气液界面捕捉模型：

式中，u_t为速度对时间的一阶导，u_t＝(u,v,w)是流体速度；F代表体积力；ρ＝ρ(x,t)是流体密度；p为流体压力；μ＝μ(x,t)是介质粘度；D是粘性应力张量；式(11)中的最后一项代表集中在相界面上的表面张力，σ代表表面张力系数；κ代表相界面的曲率；δ是DiracDelta函数；d表示计算区域中的点与相界面的垂直距离；n表示相界面上法向朝外的单位向量

流场控制方程包括质量方程(13)和动量方程(14)：

式中，表示函数对相应的变量求偏导数，ρ为流体的密度，t为时间，u_i、u_j代表流体的速度分量，x_i、x_j代表流体的位置分量，p为流场入口处压强，μ_l和μ_t分别为流体的层流和紊流粘性系数。

湍流模型为SST(Shear-Stress Transport)k-ω湍流模型，包括湍流动力粘性系数μ_t(15)、湍动能k(16)和比耗散率ω(17)：

式中，为ρ为流体的密度，α₁为经验常数，S为剪切力张量的常数项；x_j为坐标向量；u_i为速度分量；Γ_k、Γ_ω分别代表k和ω的有效耗散系数；G_k为由于平均速度梯度而导致的湍流动能k产生项；G_ω为ω产生项；Y_k、Y_ω分别表示k和ω的湍流耗散项；D_ω为交叉耗散项；S_k、S_ω为源项(用户自定义)。

刚体运动方程为包括平动方程(18)和旋转方程(19)：

M＝Iα+ω×Iω＝M_fluid+M_ext (19)

式中，m为刚体质量，为刚体的质心加速度，F为所有外力的矢量和，

F_fluid为流体作用力，g为重力加速度，α为角加速度，ω为加速度，M_fluid为流体力矩，M_ext为其他外力矩，(x_G,y_G,z_G)为刚体质心坐标。

步骤四：基于ANSYS CFX进行近水面滑行跳跃非定常流场的流固耦合数值计算，获得刚体结构与流场结构的动态变化过程。

基于ANSYS CFX进行数值计算，采用步骤三建立的流体力学模型对三维流域进行非定常流场数值计算。对计算参数进行初始化：外域的流场入口、出口及前后边界给定开放条件，上下边界为无滑移、光滑壁面边界条件；外域、中域和内域的接触面均为交界面；内域中触水物体的壁面给定为流固耦合刚体模型，并建立刚体相对于外域x和y方向位移D_x和D_y的表达式和网格质量判断公式。

D_x＝x₀[m]+abs(ave(Total Centroid Displacement X)@fly) (21)

D_y＝y₀[m]+abs(ave(Total Centroid Displacement Y)@fly) (22)

Quality＝minVal(Orthogonality Angle)@out<50[degree] (23)

式中，D_x和D_y分别为刚体相对于外域的实时距离，x₀和y₀分别为刚体相对于外域的初始位置。命令“abs”表示绝对值，“ave”表示平均值，“Total Centroid Displacement”表示刚体实际运动的位移，“@”表示提取位置，“fly”为刚体的命名。Quality表示网格的质量，“minVal”表示最小值，“Orthogonality Angle”表示网格变形的角度，“out”为外域的命名。

采用High Resolution和二阶向后差分格式，设定计算时间步长为Δt＝1×10^-4s，总时间t＝1s。将网格质量的判断公式作为中断，建立configuration进行网格控制。通过建立的流固耦合刚体模型，对刚体运动方程的平动方程(8)和旋转方程(9)进行离散求解，计算出流固耦合刚体的运动位移Total Centroid Displacement X和Total CentroidDisplacement Y。利用气液界面捕捉模型(方程(1)和方程(2))、流场控制方程(方程(3)和方程(4))和湍流模型(方程(5)、方程(6)和方程(7))进行求解，根据刚体运动位移进行流场网格的变形，并通过网格质量判断公式(23)对外域网格的质量判断。将更新的外域网格替换旧网格，继续进行流场的非定常计算，直至到达预定计算时间1s，得到近水面滑行跳跃的随时间演化的过程，即实现了对近水面滑行跳跃过程的预测。。

步骤五：基于ANSYS Workbench平台的三维流域外域网格的网格控制。

以步骤四得到的刚体相对于外域x和y方向位移D_x和D_y作为初始条件进行三维流域外域网格的网格控制，流程如图3所示。当步骤四开始数值计算时，求解器开始工作，首先判断时间是否在设定的总时间1s内，若在设定时间内，则判断网格质量，若符合要求，则继续求解，若不符合要求，则中断启动，即进行网格重构。重构过程中先调用步骤二建立的三维流域外域的ANSYS Workbench平台，读取步骤四生成的结果文件，将监测的当前时间步刚体相对于外域x和y方向位移D_x(21)和D_y(22)几何参数输出，并作为更新几何的输入几何参数D_x’和D_y’。根据读取的几何参数进行模型重构，并自行生成新的几何网格。将生成的新的几何网格替换旧网格，新网格赋值。重复步骤四至步骤五，直至达到预定求解时间1s或者外域网格质量一直不会触发中断条件。

步骤六：基于ANSYS CFX-Post对计算结果进行后处理，获得近水面跳跃流场和刚体随时间的演化信息。改变刚体的俯仰角和攻角，重复步骤一到步骤五，获得不同角度下刚体近水面跳跃随时间变化的动态信息。提取流场区域内的流动参数(包括速度、压力和水气界面)，其中速度分布通过流线图表示，压力分布通过云图表示，水气界面通过等值面或云图表示；通过提取刚体的位移、速度和加速度等信息(图4为本实施例得到的刚体触水弹跳运动轨迹)，可分析刚体触水弹跳或是下沉的原因，得出实现弹跳的最佳角度。

该实施例应用一种基于网格控制的近水面滑行跳跃流固耦合数值预测方法，对近水面运动现象进行了预测。通过改变刚体的俯仰角、攻角等设计因素，可以实现不同工况下的近水面运动的数值预测，从中选择出能够实现刚体水面弹跳的最佳设计方案，并应用到工程实践当中。由此表明，一种基于网格控制的近水面滑行跳跃数值预测方法能够对实验工况的成功性进行评估预测，提高实验的成功率、降低实验危险系数和实验成本，更加具有工程实际的应用价值。

最后需要说明的是，以上仅用以说明本发明的技术方案，本领域的普通技术人员可以对本发明的技术方案进行修改或者等同替换。凡在本发明的精神和原则之内所作修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (4)

1.基于网格控制的近水面滑行跳跃流固耦合数值预测方法，其特征在于：具体步骤如下：

步骤一：建立三维流域模型；

根据确定的触水物体的几何形状，确认其中心点为整体流域的坐标原点；然后根据坐标原点依次划分三维流域：内域、中域和外域三部分；内域是将触水物体包围起来的圆形区域，用于触水物体周围网格的加密，同时圆形内域能够方便快捷的调整触水物体的俯仰角度；中域是内域与外域的交界区域，作用于内外域网格之间的过渡，保证外域网格的精确度；外域是将所有域包起来的长方体区域，且三个流域的宽度均大于触水物体宽度；

步骤二：划分三维流域网格；

对步骤一建立的三维流域进行网格划分，其中内域采用非结构化网格对触水物体周围的流域网格进行加密处理，以便捕捉非定常流动细节；中域采用结构化网格，将中域网格由外域往内域方向逐渐加密；外域模型中将中域的位置尺寸作为参数使用，并进行非结构化网格划分；

步骤三：建立计算流体力学模型；

为了能够对步骤二所建立的三维流域网格进行流场计算求解，需要先建立计算流体力学模型；计算流体力学模型包括气液界面捕捉模型、流场控制方程、湍流模型和刚体运动方程；

气液界面捕捉模型：

式中，u_t为速度对时间的一阶导，u＝(u,v,w)是流体速度；F代表体积力；ρ＝ρ(x,t)是流体密度；p为流体压力；μ＝μ(x,t)是介质粘度；D是粘性应力张量；式(1)中的最后一项代表集中在相界面上的表面张力，σ代表表面张力系数；κ代表相界面的曲率；δ是Dirac Delta函数；d表示计算区域中的点与相界面的垂直距离；n表示相界面上法向朝外的单位向量

流场控制方程包括质量方程(3)和动量方程(4)：

式中，表示函数对相应的变量求偏导数，ρ为流体的密度，t为时间，u_i、u_j代表流体的速度分量，x_i、x_j代表流体的位置分量，p为流场入口处压强，μ_l和μ_t分别为流体的层流和紊流粘性系数；

湍流模型为SSTk-ω湍流模型，包括湍流动力粘性系数μ_t(5)、湍动能k(6)和比耗散率ω(7)：

式中，为ρ为流体的密度，α₁为经验常数，S为剪切力张量的常数项；x_j为坐标向量；u_i为速度分量；Γ_k、Γ_ω分别代表k和ω的有效耗散系数；G_k为由于平均速度梯度而导致的湍流动能k产生项；G_ω为ω产生项；Y_k、Y_ω分别表示k和ω的湍流耗散项；D_ω为交叉耗散项；S_k、S_ω为源项；

刚体运动方程为包括平动方程(8)和旋转方程(9)：

M＝Iα+ω×Iω＝M_fluid+M_ext (9)

式中，m为刚体质量，为刚体的质心加速度，F为所有外力的矢量和，F_fluid为流体作用力，g为重力加速度，α为角加速度，ω为加速度，M_fluid为流体力矩，M_ext为其他外力矩，(x_G,y_G,z_G)为刚体质心坐标；

步骤四：进行近水面滑行跳跃非定常流场的流固耦合数值计算；

采用步骤三建立的流体力学模型对三维流域近行非定常流场数值计算；在计算流体动力学(CFD)求解器中，对计算参数进行初始化；所述的对计算参数进行初始化的方法为：外域的流场入口、出口及前后边界给定开放条件，上下边界为无滑移、光滑壁面边界条件；外域、中域和内域的接触面均为交界面；内域中触水物体的壁面给定为流固耦合刚体模型，并建立刚体相对于外域x和y方向位移D_x和D_y的表达式和网格质量判断公式；

基于上述边界条件和初始条件，进行近水面滑行跳跃非定常流场的流固耦合数值计算的方法为：

步骤4.1，对流固耦合刚体进行计算，对刚体运动方程的平动方程(8)和旋转方程(9)进行离散求解，计算出流固耦合刚体的运动位移；

步骤4.2，把求得的刚体的运动位移作为流场网格的变形，对气液界面捕捉模型，即式(1)和式(2)、流场控制方程式(3)和式(4)和湍流模型式(5)、式(6)和式(7)进行离散求解，得出外域网格变形数据，即刚体在x和y方向位移；

步骤4.3，根据步骤4.2得到的外域网格变形数据和步骤四的网格质量判断公式，基于计算流体力学求解器进行网格控制，即对外域网格进行更新；

步骤4.4，用步骤4.3更新后的外域网格替换旧网格，进行流场的非定常计算；

步骤4.5，重复步骤4.1至步骤4.4，直至达到预定求解时间或者数值计算结果收敛，得到近水面滑行跳跃的随时间演化的过程，即实现了对近水面滑行跳跃过程的预测。

2.如权利要求1所述的基于网格控制的近水面滑行跳跃流固耦合数值预测方法，其特征在于：步骤4.3所述的网格控制的方法，具体步骤如下：

步骤1)，将步骤4.2监测的当前时间步的刚体在x和y方向位移的几何参数输出；

步骤2)，将步骤1)输出的几何参数作为外域新的输入几何参数；

步骤3)，对步骤2)读取的当前时间步的几何参数进行网格重构，并生成新的外域几何网格。

3.如权利要求1所述的基于网格控制的近水面滑行跳跃流固耦合数值预测方法，其特征在于：获得近水面跳跃流场和触水物体随时间的演化信息的方法为：对步骤四所得到的数值计算结果进行后处理，获得近水面跳跃随时间变化的动态信息；后处理方法为：提取流场区域内的流动参数，以及触水物体的位移、速度和加速度，则能够反应触水物体的运动轨迹和过载特性；所述流动参数包括速度、压力和水气界面，其中速度分布通过流线图表示，压力分布通过云图表示，水气界面通过等值面或云图表示。

4.如权利要求1所述的基于网格控制的近水面滑行跳跃流固耦合数值预测方法，其特征在于：步骤一所述触水物体的底面为弧形物体。