CN105955928A - 基于cfd预报船舶阻力的计算方法 - Google Patents

Abstract

本发明提供了一种基于CFD预报船舶阻力的计算方法，包括：S1：采用两相不可压RANS作为控制方程，包括连续方程和动量方程；S2：使用单相level set方法处理自由面；S3：使用有限体积法将计算域离散成一系列小单元，计算流场信息储存在网格单元重心；S4：根据单元重心的插值得到单元表面的值；S5：对单元表面的值求和得到单元体的体积积分；S6：对于离散后得到的压力速度耦合方程，采用PISO隐式分离算法进行迭代求解，其中，每一时间步的计算包括动量预估、压力求解和速度修正，最终得到收敛的数值解。本发明方法使用CFD计算得到船舶在一定运营速度下的总阻力系数，进而通过简单的公式计算得到船舶总阻力。本方法高效准确，具备很强的工程实用性。

Description

基于CFD预报船舶阻力的计算方法

技术领域

本发明涉及船舶技术领域，具体地，涉及一种对船舶阻力进行预报的计算方法。

背景技术

船舶静水航行的阻力性能预报和兴波情况模拟一直是船型优化和设计的一个热门研究课题。提高阻力计算的精确度和阻力变化内在规律的研究具有十分重要的意义。现有的研究方法有模型试验、势流理论计算及计算流体力学(CFD)方法求解等方法。

随着计算机运算速度的提高，使得CFD方法计算具有高度非线性运动特征的船舶航行阻力成为可能。在CFD计算方法中，对绕船体自由面周围黏性流场进行数值模拟的方法主要有直接数值模拟(DNS)，大渦模拟(LES)。上述两种方法均存在计算量大、精度不高的问题。

发明内容

针对现有现有技术中存在的技术问题，本发明的目的是提供一种基于CFD的预报船舶阻力的计算方法，相比其他两种方法，它不太关注湍流的脉动量，而是通过某种模型在时均化方程中体现出来，从而避免了计算量大的问题。在实际工程应用中，本发明方法计算量相对较小，精度高。当船舶的运动具有高度非线性特点和复杂的尾流时，使用本发明方法能更快、更准确的预报船舶的阻力。

为达到上述目的，本发明所采用的技术方案如下：

1、一种基于CFD的预报船舶阻力的计算方法，其特征包括：

S1：采用两相不可压缩RANS方程作为控制方程，包括连续方程和动量方程；

S2：使用SST k-ω两方程湍流模型封闭RANS方程；

S3：使用单相level set方法处理自由面；

S4：使用有限差分法离散计算域；

S5：对于离散后得到的压力速度耦合方程，采用PISO隐式分离算法进行迭代求解，其中，每一时间步的计算包括动量预估、压力求解和速度修正，最终得到收敛的数值解。

S6：作用在船体表面上的总阻力由船体表面所受应力沿船体湿表面的积分得到。

所述S1，两相不可压缩RANS方程为：

$\▿\·U=0---\left(1\right)$

$\frac{\∂\mathrm{\ρ}U}{\∂t}+\▿\·\left(\mathrm{\ρ}\right(U-U_g\left)\right)U=-\▿p_d-g\·x\▿\mathrm{\ρ}+\▿\·({\mathrm{\μ}}_{eff}\▿U)+(\▿U)\·\▿{\mathrm{\μ}}_{eff}+f_{\mathrm{\σ}}+f_s---\left(2\right)$

其中：

U代表速度场，U_g代表网格移动速度，p_d为动压力，ρ为液体密度，g为重力加速度，μ_eff表示有效动力黏性，f_σ为表面张力，f_s为消波的源项。

为了适应复杂物体形状，计算在曲线坐标系下进行，即控制方程由物理(x,y,z,t)域转化为计算域(ξ,η,τ)，则坐标转换后的连续方程和动量方程表达式如下：

$\frac{1}{J}\frac{\∂}{\∂{\mathrm{\ξ}}^j}\left(b_i^j{U}_i\right)=0---\left(3\right)$

$\frac{\∂U_i}{\∂\mathrm{\τ}}+a_{U_i}^k\frac{\∂U_i}{\∂{\mathrm{\ξ}}^k}-\frac{1}{R_{eff}}{g}^{ii}\frac{{\∂}^2{U}_i}{\∂{\mathrm{\ξ}}^i\∂{\mathrm{\ξ}}^j}=-\frac{1}{J}{b}_i^k\frac{\∂P}{\∂{\mathrm{\ξ}}^k}+s_{U_i}---\left(4\right)$

其中，J和分别是Jacobian系数和转化矩阵。

所述S2，SST k-ω两方程为：

$\frac{\∂k}{\∂t}+(U_j-{\mathrm{\σ}}_k\frac{\∂{\mathrm{\ν}}_t}{\∂x_j})\frac{\∂k}{\∂x_j}-\frac{1}{R_k}{\▿}^{2}k+s_k=0---\left(5\right)$

$\frac{\∂\mathrm{\ω}}{\∂t}+(U_j-{\mathrm{\σ}}_{\mathrm{\ω}}\frac{\∂{\mathrm{\ν}}_t}{\∂x_j})\frac{\∂\mathrm{\ω}}{\∂x_j}-\frac{1}{R_{\mathrm{\ω}}}{\▿}^2\mathrm{\ω}+s_{\mathrm{\ω}}=0---\left(6\right)$

${\mathrm{\ν}}_t=\frac{k}{\mathrm{\ω}}---\left(7\right)$

$\▿k\·n=\▿\mathrm{\ω}\·n=0---\left(8\right)$

其中，有效雷诺数及湍流产生项定义如下：

k和ω的源项如下：

s_k＝R_k(-G+β^*ωk) (9)

$s_{\mathrm{\ω}}=R_{\mathrm{\ω}}\[-\mathrm{\γ}\frac{\mathrm{\ω}}{k}G+{\mathrm{\β\ω}}^2+2(1-F_1){\mathrm{\σ}}_{\mathrm{\ω}2}\frac{1}{\mathrm{\ω}}\frac{\∂k}{\∂x_j}\frac{\∂\mathrm{\ω}}{\∂x_j}\]---\left(10\right)$

$G={\mathrm{\τ}}_{ij}\frac{\∂U_i}{\∂x_j}---\left(11\right)$

(10)式中，F1是混合函数，用来表达k-ω模型和k-ε模型混合使用的情况，充分利用两种模型的优势，即k-ω模型在固体壁面处边界条件简单，k-ε模型对自由面边界条件没有依赖性，模型中的其他常量也是关于F₁的函数，即φ＝F₁φ₁+(1-F₁)φ₂,φ₁和φ₂分别是标准k-ω模型和改进k-ε模型的常量。

所述S3，Level set函数在整个流场中连续，并满足以下运动方程：

采用单相Level set方法求解自由面，即计算一相流体的流场，另一相流体流动通过速度扩展方法确定流场速度；

求解Level set函数所满足的速度边界条件：

$\▿\mathrm{\ν}\·n=0---\left(13\right)$

其中，

满足的压力边界条件：

$p=\frac{z}{{\mathrm{Fr}}^2}---\left(14\right)$

其中，Fr为佛汝德数，z为垂向位置。

所述S4，连续方程采用有限差分法离散，对流项采用二阶迎风格式离散，黏性项采用二阶中心差分离散格式：

所有变量时间项采用二阶向后差分格式：

将时间和空间的离散表达式带入坐标转换后的RANS方程，得动量方程离散形式如下：

$A_{ijk}{U}_i^n+\underset{nb}{\Σ}{A}_{nb}{U}_{i,nb}^n=S_{U_i}-\frac{1}{J}{b}_i^k\frac{\∂p^n}{\∂{\mathrm{\ξ}}^k}---\left(18\right)$

其中，A_ijk和A_nb分别是动量离散方程中网格单元中心点和其邻点变量前的系数，

将计算域离散成一系列小单元，计算流场信息储存在网格单元重心，根据单元重心的插值得到单元表面的值，对单元表面的值求和得到单元体的体积积分。

所述S5，不可压缩RANS方程采用PISO隐式分离算法求解，首先将离散形式的动量方程带入连续方程式中，形成压力-泊松方程：

$\frac{1}{J}\frac{\∂}{\∂{\mathrm{\ξ}}^j}\left(\frac{{\mathrm{Jb}}^{jk}}{A_{ijk}}\frac{\∂P^{*}}{\∂{\mathrm{\ξ}}^k}\right)=\frac{1}{J}\frac{\∂}{\∂{\mathrm{\ξ}}^j}{b}_i^j{\hat{U}}_i---\left(19\right)$

然后通过预测、修正、再修正三步求解，即速度由上一步给定的压力场P^n-1带入动量离散方程得到：

$U_i^{*}=-\frac{1}{A_{ijk}}(\underset{nb}{\Σ}{A}_{nb}{U}_{i,nb}^{*}+S_i-\frac{1}{J}{b}_i^k\frac{\∂p^{n-1}}{\∂{\mathrm{\ξ}}^k})---\left(20\right)$

由连续方程得到的速度初步修正：

${\hat{U}}_i=\frac{1}{A_{ijk}}(S_i-\underset{nb}{\Σ}{A}_{nb}{U}_{i,nb}^{*})---\left(21\right)$

最后将带入压力-泊松方程求得p^*进行速度再修正得到：

$U_i^{**}={\hat{U}}_i-\frac{1}{{\mathrm{JA}}_{ijk}}{b}_i^k\frac{\∂P^{*}}{\∂{\mathrm{\ξ}}^k}---\left(22\right)$

所述S6，作用在船体表面上的总阻力由船体表面所受应力沿船体湿表面的积分得到，具体表达式如下：

$F=\underset{s}{\∫}{\mathrm{\τ}}_{ij}{n}_{i}ds---\left(23\right)$

作用在船体表面流体应力张量由法向压应力和切向黏性应力组成：

${\mathrm{\τ}}_{ij}=-{\mathrm{p\δ}}_{ij}+\frac{1}{\mathrm{Re}}(\frac{\∂U_i}{\∂x_j}+\frac{\∂U_j}{\∂x_i})---\left(24\right)$

其中n_i是船体表面单元的法向量，dS_i是湿面积单元，模拟船舶在沿船体纵向给定均匀流中的流场情况，因此与船体纵向的总阻力和摩擦阻力相关，无因次化后总阻力系数如下：

$C_t=\frac{F}{\frac{1}{2}{\mathrm{\ρSU}}^2}---\left(25\right)$

与现有技术相比，本发明具有如下的有益效果：

1、本发明使用单相level set方法处理自由面，相对于采用模型试验方法和势流理论计算得到的波形下求解的流场，本方法是隐式的跟踪运动界面，克服了一般波前追踪方法难于处理复杂界面，及其发生拓扑变形情况的弱点。而且计算一相流体的流场，另一相流体流动通过速度扩展方法确定流场速度，求解RANS方程的计算只在水中进行，在空气中只布置少许网格点满足界面跳跃条件，从而节省了大量的计算资源。

2、本发明采用PISO隐式分离算法进行迭代求解，与传统的SIMPLE方法相比，PISO算法在SIMPLE法的每一个时间步添加了多步矫正的算法，其更适于非定常流动的计算，在每一时间步达到收敛的计算时间更少，且数值稳定性好，结果更加准确可靠。

附图说明

通过阅读参照以下附图对非限制性实施例所作的详细描述，本发明的其它特征、目的和优点将会变得更明显：

图1是本发明计算方法一个时间步内整个计算域数值计算流程图；

图2是本发明一实施例的船模整体网格图；

图3是本发明一实施例船模计算的船表面网格。

具体实施方式

下面结合具体实施例对本发明进行详细说明。以下实施例将有助于本领域的技术人员进一步理解本发明，但不以任何形式限制本发明。应当指出的是，对本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进。这些都属于本发明的保护范围。

图1所示为本发明计算方法一个时间步内整个计算域数值计算流程图。

图2、图3所示，本发明计算方法应用于计算不同佛汝德数(Fr)下DTMB5415船模的总阻力系数的一实施例。

DTMB船模主要参数如表1所示：

表1DTMB船模主要参数

船模计算的整体网格如图2。船模计算的船表面网格如图3。

使用本发明的方法计算DTMB5415船模阻力结果与试验结果的比较,计算多个佛汝德数(Fr)下的船舶阻力，以验证本方法的有效性，计算结果如表2所示。

表2船模阻力计算结果与试验结果的比较

从阻力计算结果与试验结果的比较可以看出，总阻力系数从低于试验值0.58％至高于试验值4.55％之间，误差处在合理范围内，满足工程需要。因此，本发明计算方法能够比较准确的预报船舶的总阻力。

以上对本发明的具体实施例进行了描述。需要理解的是，本发明并不局限于上述特定实施方式，本领域技术人员可以在权利要求的范围内做出各种变形或修改，这并不影响本发明的实质内容。

Claims (7)

1.一种基于CFD预报船舶阻力的计算方法，其特征包括：

S1：采用两相不可压缩RANS方程作为控制方程，包括连续方程和动量方程；

S2：使用SST k-ω两方程湍流模型封闭RANS方程；

S3：使用单相level set方法处理自由面；

S4：使用有限差分法离散计算域；

S5：对于离散后得到的压力速度耦合方程，采用PISO隐式分离算法进行迭代求解，其中，每一时间步的计算包括动量预估、压力求解和速度修正，最终得到收敛的数值解。

S6：作用在船体表面上的总阻力由船体表面所受应力沿船体湿表面的积分得到。

2.根据权利要求1所述的基于CFD预报船舶阻力的计算方法，其特征在于，所述S1，两相不可压缩RANS方程为：

▽·U＝0 (1)

其中：

U代表速度场，U_g代表网格移动速度，p_d为动压力，ρ为液体密度，g为重力加速度，μ_eff表示有效动力黏性，f_σ为表面张力，f_s为消波的源项。

为了适应复杂物体形状，计算在曲线坐标系下进行，即控制方程由物理(x,y,z,t)域转化为计算域则坐标转换后的连续方程和动量方程表达式如下：

其中，J和分别是Jacobian系数和转化矩阵。

3.根据权利要求1所述的基于CFD预报船舶阻力的计算方法，其特征在于，所述S2，SSTk-ω两方程为：

▽k·n＝▽ω·n＝0 (8)

其中，有效雷诺数及湍流产生项定义如下：

k和ω的源项如下：

s_k＝R_k(-G+β^*ωk) (9)

(10)式中，F1是混合函数，用来表达k-ω模型和k-ε模型混合使用的情况，充分利用两种模型的优势，即k-ω模型在固体壁面处边界条件简单，k-ε模型对自由面边界条件没有依赖性，模型中的其他常量也是关于F₁的函数，即φ＝F₁φ₁+(1-F₁)φ₂,φ₁和φ₂分别是标准k-ω模型和改进k-ε模型的常量。

4.根据权利要求1所述的基于CFD预报船舶阻力的计算方法，其特征在于，所述S3，Level set函数在整个流场中连续，并满足以下运动方程：

采用单相Level set方法求解自由面，即计算一相流体的流场，另一相流体流动通过速度扩展方法确定流场速度；

求解Level set函数所满足的速度边界条件：

▽ν·n＝0 (13)

其中，

满足的压力边界条件：

其中，Fr为佛汝德数，z为垂向位置。

5.根据权利要求2所述的基于CFD预报船舶阻力的计算方法，其特征在于，所述S4，连续方程采用有限差分法离散，对流项采用二阶迎风格式离散，黏性项采用二阶中心差分离散格式：

所有变量时间项采用二阶向后差分格式：

将时间和空间的离散表达式带入坐标转换后的RANS方程，得动量方程离散形式如下：

其中，A_ijk和A_nb分别是动量离散方程中网格单元中心点和其邻点变量前的系数，

将计算域离散成一系列小单元，计算流场信息储存在网格单元重心，根据单元重心的插值得到单元表面的值，对单元表面的值求和得到单元体的体积积分。

6.根据权利要求1所述的基于CFD预报船舶阻力的计算方法，其特征在于，所述S5，不可压缩RANS方程采用PISO隐式分离算法求解，首先将离散形式的动量方程带入连续方程式中，形成压力-泊松方程：

然后通过预测、修正、再修正三步求解，即速度由上一步给定的压力场P^n-1带入动量离散方程得到：

由连续方程得到的速度初步修正：

最后将带入压力-泊松方程求得p^*进行速度再修正得到：

。

7.根据权利要求1所述的基于CFD预报船舶阻力的计算方法，其特征在于，所述S6，作用在船体表面上的总阻力由船体表面所受应力沿船体湿表面的积分得到，具体表达式如下：

作用在船体表面流体应力张量由法向压应力和切向黏性应力组成：

其中n_i是船体表面单元的法向量，dS_i是湿面积单元，模拟船舶在沿船体纵向给定均匀流中的流场情况，考虑船体纵向的总阻力和摩擦阻力，无因次化后总阻力系数如下：

。