CN111027152B - 一种围绕水翼空化流固耦合数值模拟方法 - Google Patents

Abstract

本申请涉及一种围绕水翼空化流固耦合数值模拟方法，包括：确定围绕三维水翼的流体域，并对所述流体域进行网格划分；建立空化流动的流体运动方程组；根据所述空化流动的流体运动方程组和预设边界条件得到所述流体域的空化数值模拟结果；根据所述空化数值模拟结果，以及通过径向基函数的插值方法进行非定常流固耦合计算，得到所述三维水翼的结构振动特性与所述流体域的空化流动规律。该技术方案本申请实施例提供的方法，通过基于悬臂梁模态的三维水翼结构，建立了能反映复杂空化流动和三维水翼结构的时空演变规律的数值模拟方法，进而实现提高空化流固耦合模拟精度。

Description

一种围绕水翼空化流固耦合数值模拟方法

技术领域

本申请涉及流固耦合数值模拟领域，尤其涉及一种围水翼空化流固耦合数值模拟方法。

背景技术

空化是水力机械、海洋工程领域的一种极为复杂的非定常现象，伴随着空化初生、发展、脱落和溃灭，会造成水动力的时空演化，特别是空泡溃灭所产生的瞬态高压脉冲会对结构造成巨大的冲击载荷。反过来，结构形状的改变也会影响空化的流动特性。这是一种典型的流体与结构相互作用的流固耦合现象。与传统的金属材料相比，复合材料具有比重小、比强度大等优点，使其开始在水力机械、海洋工程领域中广泛应用。在这种情况下，空化流激振动成为影响结构安全和稳定运行的不可忽视的问题。

针对绕弹性水翼的空化流固耦合问题，很多学者将三维水翼的弯曲、扭转简化为二维水翼在平面上的平动、转动，无法真实的模拟空化所造成的水动力载荷时空分布，也无法反映水翼的三维扭转效应。因此，发展针对三维弹性水翼的高效率的空化流固耦合数值模拟方法具有重要性，也具有必要性。

发明内容

为了解决上述技术问题或者至少部分地解决上述技术问题，本申请提供了一种围绕水翼空化流固耦合数值模拟方法。

本申请实施例提供了一种围绕水翼空化流固耦合数值模拟方法，包括：

确定围绕三维水翼的流体域，并对所述流体域进行网格划分；

建立空化流动的流体运动方程组；

根据所述空化流动的流体运动方程组和预设边界条件得到所述流体域的空化数值模拟结果；

根据所述空化数值模拟结果，以及通过径向基函数的插值方法进行非定常流固耦合计算，得到所述三维水翼的结构振动特性与所述流体域的空化流动规律。

在一个可能的实施方式中，所述确定围绕三维水翼的流体域，包括：

确定所述三维水翼的结构；

根据所述三维水翼的结构建立所述三维水翼的流体域。

在一个可能的实施方式中，所述建立空化流动的流体运动方程组，包括：

根据质量守恒以及动量守恒建立第一控制方程组，所述第一控制方程组包括：连续性方程，动量方程以及含汽率输运方程；

基于空泡动力学方法建立空化模型；

对所述连续性方程以及所述动量方程分别进行平均得到第二控制方程组；

根据所述空化模型以及第二控制方程组建立所述流体运动方程组。

在一个可能的实施方式中，

式中，

为求偏导数，ρ为混合相密度，u_j为在j方向上的速度量，x_j为位置分量，t为时间；

所述动量方程为：

式中，

为求偏导数，u_i为在i方向上的速度量，u_j为在j方向上的速度量，p为压力，μ为流体粘度，x_j为坐标分量，x_i为坐标分量；

所述含汽率输运方程为：

式中，m⁺为空化的蒸发速率，m^-为空化的凝结速率，α_v为汽相的体积分数，ρ_v为汽相密度。

在一个可能的实施方式中，所述空化模型用于计算空化的蒸发速率以及凝结速率；

所述空化的蒸发速率m⁺的计算过程如下：

所述空化的凝结速率m^-的计算过程如下：

式中，R_b为空泡半径，p_v为蒸汽压力，N_b为空泡密度，ρ_l为液相密度。

在一个可能的实施方式中，所述第二控制方程组，包括：

式中，τ_ij为亚格子应力项，μ_t为涡粘系数。

在一个可能的实施方式中，所述根据所述空化流动的流体运动方程组和预设边界条件得到所述流体域的空化数值模拟结果，包括：

确定所述预设边界条件，所述预设边界条件包括：所述流体域的进口面设为来流速度，所述流体域的出口面设为压力边界，所述流体域的上下面、两侧面、三维水翼的表面均设置为无滑移壁面；

按照库朗数小于等于5确定时间步长，采用所述流体运动方程组以及所述预设边界条件进行非定常空化流场的数值计算，得到所述流体域的空化数值模拟结果。

在一个可能的实施方式中，所述方法还包括：

建立悬臂梁模态的所述三维水翼的结构运动方程；

其中，所述结构运动方程包括归一化振型以及对应的固有频率，所述结构运动方程如下：

其中，M为质量矩阵，K为刚度矩阵，η为广义位移，

为加速度、u为位移，u＝φη，φ为质量归一化振型矩阵，如下式：

φ_i(x)＝cosβ_ix-coshβ_ix-θ_i(sinβ_ix-sinhβ_ix)

式中，

杨氏模量E，结构密度ρ_s；

所述固有频率ω_i如下：

式中，I为惯性矩，S是三维水翼沿着流向的投影面积，l是水翼的展长。

在一个可能的实施方式中，所述根据所述空化数值模拟结果，以及通过径向基函数的插值方法进行非定常流固耦合计算，得到所述三维水翼的结构振动特性与所述流体域的空化流动规律，包括：

基于径向基函数建立所述三维水翼中结构节点与所述流体域中流体网格点的插值关系；

根据所述空化数值模拟结果获取分布于所述流体网格点的第一作用力；

根据虚功等效原理确定所述第一作用力与分布于所述结构节点上第二作用力的关系，并计算所述第二作用力；

对所述结构运动方程进行求解得到所述结构节点的变形位移；

基于所述结构节点与所述流体网格点的插值关系更新所述流体网格点的坐标；

基于空化流动的流体运动方程组更新当前时间步的空化数值模拟结果；

重复上述步骤依次迭代直至满足预设数据交换次数时，确认完成一个物理时间步的计算；

当所述物理时间步达到预设数量时，得到所述三维水翼的结构振动特性与所述流体域的空化流动规律。

本申请实施例提供的上述技术方案与现有技术相比具有如下优点：

本申请实施例提供的方法，通过基于悬臂梁模态的三维水翼结构，建立了能反映复杂流体域和三维水翼结构的时空演变规律的数值模拟方法，进而实现提高空化流固耦合模拟精度。

本申请实施例提供的方法，实现围绕水翼的空化流固耦合数值模拟，有助于对空化流固耦合问题的流动规律和结构响应特性进行深入的研究，能够应用于围绕水翼空化流固耦合领域，并能够解决相关工程问题。

附图说明

此处的附图被并入说明书中并构成本说明书的一部分，通过附图示出了符合本发明的实施例，并与说明书一起用于解释本发明的原理。

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，对于本领域普通技术人员而言，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本申请实施例提供的一种围绕水翼空化流固耦合数值模拟方法的流程图。

具体实施方式

为使本申请实施例的目的、技术方案和优点更加清楚，下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本申请的一部分实施例，而不是全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本申请保护的范围。

图1为本申请实施例提供的一种围绕水翼空化流固耦合数值模拟方法的流程图。如图1所示，该方法包括以下步骤：

步骤S11，确定围绕三维水翼的流体域，并对流体域进行网格划分；

步骤S12，建立空化流动的流体运动方程组；

步骤S13，根据空化流动的流体运动方程组和预设边界条件得到流体域的空化数值模拟结果；

步骤S14，根据空化数值模拟结果，以及通过径向基函数的插值方法进行非定常流固耦合计算，得到三维水翼的结构振动特性与流体域的空化流动规律。

本实施例基于悬臂梁模态的三维水翼结构，提出了能反映复杂流体域和三维水翼结构的时空演变规律的数值模拟方法，进而实现提高空化流固耦合模拟精度。

本实施例中，获取围绕三维水翼的流体域的具体实现方式如下：获取三维水翼的结构，根据三维水翼的结构获取三维水翼的包围参数，基于包围参数确定围绕所述三维水翼的流体域。

具体的，包围参数如下：流体域的进口面距离水翼前缘为n倍的水翼弦长，流体域的出口面距离水翼尾缘为m倍的水翼弦长，流体域的上下面分别距离水翼弦线中点为p,q倍的水翼弦长，流体域的一侧面与水翼根部同一个平面，流体域的另一侧面距离水翼尖端d。可采用PRO/E、SolidWorks、UG等几何建模软件。流体域大小(n,m,p,q,d)可根据试验段尺寸、计算精度和效率来综合确定。

对围绕三维水翼的流体域进行网格划分，并对靠近水翼的网格进行加密，使其满足空化流动的流体运动方程组的要求，以便于捕捉空化流动细节。空化流动细节包括准确预测空化的初生、发展、脱落和溃灭等关键过程，水翼升力和阻力瞬时特性，水翼壁面压力脉动等。可采用ICEM、Pointwise、Hypermesh、snapHexMesh等软件来实现。

本实施例中，基于所述流体域建立流体运动方程组，包括：

(一)根据质量守恒以及动量守恒建立第一控制方程组，第一控制方程组包括：连续性方程，动量方程以及含汽率输运方程；

连续性方程为：

式中，

为求偏导数，u_j为在j方向上的速度量，ρ为混合相密度，μ为流体粘度，x_j为位置分量，t为时间。

动量方程为：

式中，p为压力，u_i为在i方向上的速度量，u_j为在j方向上的速度量，μ为流体粘度，x_j为位置分量，x_i为位置分量。

含汽率输运方程为：

式中，m⁺为空化的蒸发速率，m^-为空化的凝结速率，α_v为汽相的体积分数，ρ_v为汽相密度。

(二)基于空泡动力学方法建立空化模型；

本实施例中的基于空泡动力学方法来推导空化模型，即根据Rayleigh-Plesset方程来计算单一空泡体积的变化速率。其中，空化源项包括蒸发速率和凝结速率，分别反映空化的生成和凝结。本发明采用Schnerr and Sauer空化模型，模型参数包括饱和蒸汽压力p_v，空泡半径R_b。

空化的蒸发速率m⁺的计算过程如下：

空化的凝结速率m^-的计算过程如下：

式中，R_b为空泡半径，p_v为蒸汽压力，N_b为空泡密度，ρ_l为液相密度。

(三)对连续性方程以及动量方程进行平均计算得到第二控制方程组；

其中对连续性方程以及动量方程进行Favre平均，得到大涡模拟的控制方程组(即第二控制方程组)，第二控制方程组包括：

式中，τ_ij为亚格子应力项，其计算过程如下：

μ_t为涡粘系数，其计算过程如下：

其中，

Δ_s＝min(kd,CsV^1/3)。

(四)根据空化模型以及第二控制方程组建立流体运动方程组。

本实施例中，在建立流体运动方程组后，根据流体运动方程组和预设边界条件得到流体域的空化数值模拟结果。具体实现方式如下：

获取预设边界条件，预设边界条件包括：流体域的进口面设为来流速度，流体域的出口面设为压力边界，流体域的上下面、两侧面、三维水翼的表面均设置为无滑移壁面；

按照库朗数小于等于5确定时间步长，采用空化流动的流体运动方程组以及预设边界条件进行非定常空化流场的数值计算，得到流体域的空化数值模拟结果。

本实施例在根据空化数值模拟结果进行非定常流固耦合计算之前，还需建立基于悬臂梁模态的三维水翼的结构运动方程，结构运动方程包括归一化振型以及对应的固有频率。其中，结构运动方程如下：

其中，M为质量矩阵，K为刚度矩阵，η为广义位移，

为加速度、u为位移，u＝φη，φ为质量归一化振型矩阵，如下式：

φ_i(x)＝cosβ_ix-coshβ_ix-θ_i(sinβ_ix-sinhβ_ix)

式中，

β₁l＝1.875,β₂l＝4.694,β₃l＝7.855，杨氏模量E，结构密度ρ_s；

所述固有频率如下：

式中，I为惯性矩，S是三维水翼沿着流向的投影面积，l是水翼的展长。

本实施例中，根据所述空化数值模拟结果，以及通过径向基函数的插值方法进行非定常流固耦合计算，得到所述三维水翼的结构振动特性与所述流体域的空化流动规律。具体包括以下步骤：

(一)基于径向基函数建立三维水翼中结构节点与流体域中流体网格点的插值关系；

由于流体域与三维水翼在交界面上的网格节点不是一一对应，所以两者需要通过插值的方式交换力与位移的信息，数据插值需要满足虚功等效，即水动广义力等于结构广义力。本发明的数据交换方法基于径向基函数(RBF)，如下式：

其中，X＝(x,y,z)，||x-x_i||²＝(x-x_i)²+(y-y_i)²+(z-z_i)²，紧支撑的径向基函数为：

其中，d_i＝||x-x_i||，R为支撑半径，表示RBF插值的径向影响范围。φ表示基函数，可选用Gaussian、Thin plate spline(TPS)、Hardy’s multiquadric(HMQ)、Wendland’s C0(C0)、Euclid’s Hat(EH)等，如下表所示：

表1径向基函数类型

本实施例中，以x方向的分量为例叙述。把N个结构节点坐标代入公式(8)，联合公式(9)，得到：

φ_s1s2＝φ(||x_s1-x_s2||)

其中，下标"s"表示结构，Δx_s表示结构节点的位移。求解方程(11)，得到系数λ，再把流体网格点坐标代入式(8)，就得到流体网格点的位移：

ΔX_a＝A_asλ

其中，下标"a"表示流体，M表示物面流体节点的数目。

(二)根据空化数值模拟结果获取分布于流体网格点的第一作用力，本实施例中，通过对三维水翼的结构表面的粘性力和压力积分得到流体作用力向量f_a，即第一作用力。

(三)根据虚功等效原理确定第一作用力与分布于结构节点上第二作用力的关系，并计算第二作用力；

把f_a插值到结构节点上得到f_s。力的插值需要满足虚功等效，即

但该方程不能直接求解。通过引入四个辅助计算的系数

则虚功等效方程可以重新表示为：

把式(11)和式(12)代入式(13)，得到

求解线性方程(14)，就可以得到结构节点上的力，即第二作用力f_s。线性方程(14)的前四行是：

这表明载荷的传递满足力守恒和力矩守恒。

(四)对结构运动方程进行求解得到结构节点的变形位移。

(五)基于结构节点与流体网格点的插值关系更新流体网格点的坐标；由于CSD物面的运动会引起CFD物面的变化，这就需要对CFD网格进行更新。常见的网格变形实现方法包括基于弹簧网络的网格光滑算法(Spring-based Method)、局部网格重构算法(LocalRemeshing Method)、PDE方法，Delaunay方法以及RBF方法等。

(六)基于空化流动的流体运动方程组更新当前时间步的空化数据值模拟结果。

(七)重复上述步骤(二)至步骤(六)依次迭代直至满足预设数据交换次数时，确认完成一个物理时间步的计算。

然后重新执行步骤(二)，进行下一个物理时间步的计算，当物理时间步达到预设数量时，得到三维水翼结构振动特性与流体域的空化流动规律。

进一步需要说明的是，在本文中，诸如“第一”和“第二”等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的过程、方法、物品或者设备中还存在另外的相同要素。

以上所述仅是本发明的具体实施方式，使本领域技术人员能够理解或实现本发明。对这些实施例的多种修改对本领域的技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所申请的原理和新颖特点相一致的最宽的范围。

Claims (4)

1.一种围绕水翼空化流固耦合数值模拟方法，其特征在于，包括：

获取基于悬臂梁模态的三维水翼的结构；

根据三维水翼的结构获取三维水翼的包围参数，基于包围参数确定围绕所述三维水翼的流体域，并对所述流体域进行网格划分；

根据质量守恒以及动量守恒建立第一控制方程组，所述第一控制方程组包括：连续性方程，动量方程以及含汽率输运方程；

基于空泡动力学方法建立空化模型；

对所述连续性方程以及所述动量方程分别进行平均得到第二控制方程组；

根据所述空化模型以及第二控制方程组建立所述流体运动方程组；

建立基于悬臂梁模态的所述三维水翼的结构运动方程；

其中，所述结构运动方程包括归一化振型以及对应的固有频率，所述结构运动方程如下：

其中，M为质量矩阵，K为刚度矩阵，η为广义位移，

为加速度、u为位移，u＝φη，φ为质量归一化振型矩阵，如下式：

φ_i(x)＝cosβ_ix-coshβ_ix-θ_i(sinβ_ix-sinhβ_ix)

式中，

其中，β₁l＝1.875,β₂l＝4.694,β₃l＝7.855，E为杨氏模量，ρ_s为结构密度；

所述固有频率ω_i如下：

式中，I为惯性矩，S是三维水翼沿着流向的投影面积，l是水翼的展长；

根据所述空化流动的流体运动方程组和预设边界条件得到所述流体域的空化数值模拟结果；

根据所述空化数值模拟结果，以及通过径向基函数的插值方法进行非定常流固耦合计算，得到所述三维水翼的结构振动特性与所述流体域的空化流动规律；

所述连续性方程为：

式中，ρ为混合相密度，

为求偏导数，u_j为在j方向上的速度量，x_j为位置分量，t为时间；

所述动量方程为：

式中，

为求偏导数，u_i为在i方向上的速度量，u_j为在j方向上的速度量，p为压力，μ为流体粘度，x_j为坐标分量，x_i为坐标分量；

所述含汽率输运方程为：

式中，m⁺为空化的蒸发速率，m^-为空化的凝结速率，α_v为汽相的体积分数，ρ_v为汽相密度；

所述空化模型用于计算空化的蒸发速率以及凝结速率；

所述空化的蒸发速率m⁺的计算过程如下：

所述空化的凝结速率m^-的计算过程如下：

式中，R_b为空泡半径，p_v为蒸汽压力，N_b为空泡密度，ρ_l为液相密度。

2.根据权利要求1所述的方法，其特征在于，所述第二控制方程组，包括：

式中，τ_ij为亚格子应力项，μ_t为涡粘系数。

3.根据权利要求2所述的方法，其特征在于，所述根据所述空化流动的流体运动方程组和预设边界条件得到所述流体域的空化数值模拟结果，包括：

确定所述预设边界条件，所述预设边界条件包括：所述流体域的进口面设为来流速度，所述流体域的出口面设为压力边界，所述流体域的上下面、两侧面、三维水翼的表面均设置为无滑移壁面；

按照库朗数小于等于5确定时间步长，采用所述流体运动方程组以及所述预设边界条件进行非定常空化流场的数值计算，得到所述流体域的空化数值模拟结果。

4.根据权利要求3所述的方法，其特征在于，所述根据所述空化数值模拟结果，以及通过径向基函数的插值方法进行非定常流固耦合计算，得到所述三维水翼的结构振动特性与所述流体域的空化流动规律，包括：

基于径向基函数建立所述三维水翼中结构节点与所述流体域中流体网格点的插值关系；

根据所述空化数值模拟结果获取分布于所述流体网格点的第一作用力；

根据虚功等效原理确定所述第一作用力与分布于所述结构节点上第二作用力的关系，并计算所述第二作用力；

对所述结构运动方程进行求解得到所述结构节点的变形位移；

基于所述结构节点与所述流体网格点的插值关系更新所述流体网格点的坐标；

基于空化流动的流体运动方程组更新当前时间步的空化数值模拟结果；

重复上述步骤依次迭代直至满足预设数据交换次数时，确认完成一个物理时间步的计算；

当所述物理时间步达到预设数量时，得到所述三维水翼的结构振动特性与所述流体域的空化流动规律。

Images