非均匀流场螺旋桨空泡数值模拟

图1 计算区域(m)

Fig.1 Computational field (m)

$式中:下标i 和 j 分别为坐标方向; 下标 l 、 v 分别指液相和气相; t 为时间; ρ_{m} 为混合流密度; u 和 p 分别为速度和压力; μ_{m} 为混合相的动力黏性系数; μ_{t, m} 为湍流黏性系数; α 为相体积分数 . ρ_{m} 和 μ_{m}$ 均通过体积加权平均后得到.采用可实现的k-ε湍流模型求解涡黏性系数,该模型适用于强旋转流动和剪切流,文献^[2-5,7]采用该湍流模型得到了较为满意的结果.

由于控制方程中增加体积分数这一未知量,因此需通过添加额外的方程进行封闭.相变过程可认为起始于溶解于水中的气核,单位控制体积内混合流的质量始终保持不变,根据气核的生长和溃灭速率可计算两相间的质量输运率和体积分数.基于Rayleigh-Plesset方程的Schnerr-Sauer空化模型能够对气液两相间的质量输运进行较好的描述^[11].体积分数的控制方程为

(5)

\begin{matrix} \frac{\partial (ρ_{v} α_{v})}{\partial t} + \frac{\partial (ρ_{v} α_{v} u_{i})}{\partial x_{i}} = {\overset{\cdot}{m}}^{+} + {\overset{\cdot}{m}}^{-} \end{matrix}

(6)

α_{v} = \frac{n_{0} \frac{4}{3} π R^{3}}{1 + n_{0} \frac{4}{3} π R^{3}}

${\overset{\cdot}{式中:m}}^{+} 和 {\overset{\cdot}{m}}^{-} 分别为汽化和凝结引起的质量变化率; n_{0} 为初始单位体积内汽泡个数 (本文研究蒸汽泡, 简称汽泡), 值取 10^{12}; R$ 为汽泡半径,初始半径为10^-6m.该空化模型认为单个汽泡的尺寸变化率与当地静压及饱和蒸气压相关,忽略表面张力和黏性的求解方程为

(7)

\frac{\partial R}{\partial t} = \sqrt[]{\frac{2 |p - p_{0}|}{3}}

式中:p₀为饱和蒸气压,取值 3850 Pa,试验时水温为28.4 ℃.随着空化逐渐发展,单位体积内汽泡个数N₀=α_ln₀,保证液相中汽泡密度始终保持不变.

当p<p₀时,液体汽化,质量输运速率为

(8)

\begin{array}{l} {\overset{\cdot}{m}}^{+} = N_{0} \frac{\partial (ρ_{v} 4 π R^{3} / 3)}{\partial t} = 4 π N_{0} ρ_{v} R^{2} \frac{\partial R}{\partial t} = \\ 4 π N_{0} ρ_{v} R^{2} \sqrt[]{\frac{2}{3} \max {(p_{0} - p), 0.005}} \end{array}

当p>p₀时,气体凝结,质量输运速率为

(9)

{\overset{\cdot}{m}}^{-} = - 4 π N_{0} ρ_{v} R^{2} \sqrt[]{\frac{2}{3} \max {(p - p_{0}), 0.01}}

2 研究对象、网格划分及可信性校验

2.1 研究对象

本文以某万吨级油船螺旋桨模型为研究对象,该桨模已在上海船舶运输科学研究所(SSSRI)空泡水筒做过大量试验^[12].螺旋桨模型直径(D)为248 mm,盘面比为0.4,螺距比_(0.35_D₎为0.826,桨叶数为4,旋向为右旋.

2.2 网格划分及可信性校验

计算区域与试验条件完全相同,SSSRI空泡水筒工作段长2.6 m,横截面呈正方形带0.1 m倒圆角,宽为0.6 m,高为0.6 m,螺旋桨位于距入口0.5 m处,见图1,原点坐标位于入口轴系中心处.

网格采用六面体核心网格Trim类型,该模型内置有一个用以控制表面模板单元格尺寸的上限,这个上限可避免严重的模型几何失真和生成质量差的网格.整个计算域分为静止区域和螺旋桨旋转区域,两个区域通过交界面连接.为能捕捉到梢涡空泡,需对桨叶梢附近(涡量较大的区域)进行网格加密.涡量大小的判断采用Q准则:

(10)

\begin{matrix} Q = \frac{1}{2} (‖Ω ‖^{2} - ‖ S ‖^{2}) \end{matrix}

式中: $Ω 为涡量; S$ 为漩涡变形率,均为二阶张量.图2给出了Q等值面图,从图中可看出,螺旋涡管从桨叶导边起始,由叶梢处脱出,与试验梢涡空泡起始的位置相同.再进行空化流场计算时发现,当叶背片空泡产生时,由于片空泡具有一定的厚度,此时涡管在桨叶上的脱出位置与无空化状态有所改变,产生向x轴负方向的微小位移.基于上述现象在对该涡管进行网格加密时,加密区域应包围此涡管,为降低网格数量,仅对旋转区域内螺旋涡管进行了加密.

图2

图2 Q等值面(Q=5×10⁴ s^-2)

Fig.2 Isosurface of Q(Q=5×10⁴ s^-2)

对桨的空化模拟过程进行不确定分析,共生成3套计算网格,分别为网格1、2、3.网格加细比为 $\sqrt[]{2}$ .对应网格数量分别为174 万、400 万和1000万.网格2具体划分如下:

(1) 静止区域最大网格单元尺寸为0.08D,最小尺寸为0.008D(交界面网格尺寸).为保证桨盘面处伴流场,从入口至下游0.7 m处,以轴系为中心1.62D直径范围内进行网格加密,网格尺寸为0.02D,网格单元总数约为123 万.

(2) 旋转区域最大网格单元尺寸为0.008D,最小网格尺寸为0.002D(桨叶表面),桨叶设置为无滑移壁面,边界层内第1层网格高度为2.76×10^-4D,共6层,总厚度为4.03×10^-4D,本文计算工况下浆叶表面无量纲法向距离y⁺均在30~100的范围内.为捕捉梢涡空化,对桨叶梢及螺旋管进行网格加密,尺寸为0.001D,网格单元总数约为277 万.网格划分见图3.

图3

图3 计算域网格(网格2)

Fig.3 Mesh of computational field (mesh 2)

计算入口边设置为速度入口,STAR-CCM+允许用户配置文件表设置入口流速分布.由于试验采用金属网格方法模拟船后伴流场,因此本文仅给定轴向 $(x 轴方向) 流速分布, 试验测量得到的桨盘面处无量纲轴向速度 W_{x} = 1 - u_{x} / {\bar{u}}_{x} (u_{x} 、 {\bar{u}}_{x}$ 分别为轴向瞬时速度和轴向平均速度)分布,见图4,图中角度为桨叶的相位角(视角为从船艉向船艏看),作为本文计算的输入条件,后文中θ亦指此角.出口设置为压力出口.水筒壁面、轴系及桨叶设置为无滑移壁面,交界面为滑移壁面.

图4

图4 试验测得的轴向伴流分布

Fig.4 Distribution of axial wake measured by the test

为对计算结果进行可信性校验,本文运用安全系数法^[13]对3套网格计算得到的螺旋桨平均推力系数 $K_{T} 进行不确定度估算, K_{T} = T / (ρ_{l} n^{2} D^{4}) (T 为螺旋桨推力, n 为螺旋桨转速) . 网格 1 、 2 、 3 计算得到的 K_{T}$ 分别为 0.1393、0.1388、0.1385,试验值为 0.1380,不确定度估算结果为4.464×10^-4.分析计算结果可知,随着网格数量的增加,推力系数逐渐降低且趋于试验值,但是网格1与2、网格2与3以及网格3与试验值之间的差异在变小.不确定度值相当小,可认为计算结果可信.图5给出了不同网格计算得到的空泡形态.分析图5可知网格1不能有效地捕捉梢涡空泡,网格2和3均能捕捉到梢涡空泡,叶背片空泡正投影面积和梢涡空泡形态一致.综合考虑计算消耗和模拟精度之后,选用网格2进行后续计算.

图5

图5 不同网格模拟的空泡形态与试验结果对比

Fig.5 Comparisons of cavity pattern between calculation and test results at different mesh sizes

3 结果分析

3.1 无空化流场数值模拟

螺旋桨数值计算条件、工况与试验条件完全相同^[12].试验中n为30 r/s,K_T为0.138,背景压力为p_ref,通过改变背景压力满足转速空泡数σ_n=(p_ref-p₀)/(0.5ρ_ln²D²)=2.96.在正式计算之前,无螺旋桨情况下,需确认在速度进口处给定的流速分布,流向桨盘面位置时速度分布是否发生改变.图6给出了入口处和桨盘面处的轴向伴流分布,分析图6可知在桨盘面处流场分布与入口几乎相同,因此可按该网格划分方法进行后续计算.

图6

图6 不同断面处轴向伴流分布

Fig.6 Distribution of axial wake at different cross-sections

为保证推力系数与试验相同,又能捕捉到梢涡空化,因此先进行无空化流场计算.采用滑移网格方法,不加入空化模型,计算时间步长设定为 2.7778×10^-4s,对应螺旋桨旋转角度为3°.背景压力设置为1个大气压,调整来流速度(轴向伴流分布保持不变)直至时间平均推力系数相等.无空化状态K_T为 0.1388,误差为0.6%.螺旋桨旋转一圈推力系数变化见图7,分析图7可知一个旋转周期内推力系数出现4个峰值,与桨叶数相对应.由图3可知从θ=45° 附近开始流场逐渐趋于均匀,推力系数曲线也在此出现波谷,无空化状态波谷θ分别为40°、130°、220°及310°.

图7

图7 螺旋桨旋转一圈K_T变化

Fig.7 Variations of K_T during one propeller rotation

3.2 空化流场数值模拟

当无空化流场计算稳定后,激活欧拉多相流模型和空化模型,调整环境压力使转速空泡数为2.96,计算时间步长为 9.25926×10^-5 s,对应螺旋桨旋转角度为1°.空化状态K_T为0.139,误差为0.72%,推力脉动高于无空化状态,最小推力系数出现在θ=45°,135°,225°,315°.为验证数值计算结果的准确性,与试验结果比较,对不同相位角桨叶空泡形态进行了对比分析,见图8.采用α_v=0.1等值面描述空泡形态.图中给出了桨叶进、出伴流区域空泡发展及溃灭的整个过程:当桨叶开始进入伴流区域时首先出现叶背片空泡,起始于桨叶导边附近,随着桨叶继续旋转,叶背片空泡面积逐渐增大,梢涡空泡开始脱出,最大叶背片空泡面积出现在θ=10° 的位置,随后片空泡面积逐渐减小,梢涡空化逐渐加强,随着桨叶逐渐出伴流区域,叶背片空泡消失,梢涡逐渐降弱直至消失.数值计算准确的预测了这一历程,片空泡初生、片空泡面积、梢涡初生及梢涡消失等现象与试验结果吻合.

图8

图8 不同θ时空泡形态计算结果与试验结果比较

Fig.8 Comparisons of cavity pattern between calculation results and test results observation at different θ values

但是对于梢涡空化而言,空化形态略有差异,例如θ=60° 时,试验照片显示桨叶脱出的梢涡并不平滑,刚从叶梢梢部脱落时,梢涡非常细,随后梢涡空化加强,涡管变粗且外围呈现雾状空化,涡管表面非常粗糙.这可能是由于梢涡脱出后进入高伴流区域,速度梯度变大,压力非常不稳定,致使空化呈现出强烈的非定常特性,在短时间内剧烈变化.计算中θ=60° 时梢涡空化沿螺旋涡管也是由弱变强,从叶梢附近脱出、消失再出现,但是梢涡空间结构更加平顺光滑.这主要是由于网格尺寸不能无限增加,再者RANS方法对该微小尺度涡不能够准确预测.即使采用LES方法,其梢涡涡管依然平顺^[7].

试验中空泡手描图能够对叶背片空泡的大小、形态及稳定性进行详尽描述,从θ=20° 开始梢涡在空间内呈现复杂的三维结构,手描图难以在二维平面内对其准确的描述,对于梢涡形态手描图仅进行了示意,因此仅对叶背片空泡进行了量化计算,分析了与试验结果的差异性.现将不同相位角空泡计算结果、试验手描图和试验照片进行图像处理,每张照片均进行图像重塑为正投影图片,分别测绘叶背空泡正投影平面面积和桨叶正投影平面面积.定义叶背片空泡面积百分比γ=A_S/A_B(A_S为叶背空泡面积,A_B为桨叶面积).表1给出了数值和试验叶背片空泡正投影面积百分比结果,分析表1可知数值计算结果与试验结果吻合,仅在θ=40° 时相差5%左右,其他情况下计算结果与试验结果的差值在2.5%以内.

表1 计算和试验叶背片空泡面积对比

Tab.1 Comparisons of back-sheet cavity area between calculation and test results

θ/(°)	γ/%
θ/(°)	试验值	计算值
340	6.18	8.60
350	15.98	17.71
0	21.39	21.90
10	25.87	25.16
20	25.50	23.02
30	18.68	20.58
40	5.73	10.46
50	1.54	1.98
60	0.00	0.00

新窗口打开| 下载CSV

4 结语

本文以较少的网格数量成功模拟了非均匀流场螺旋桨叶背片空泡和梢涡空泡.结果表明:Schnerr-Sauer空化模型和可实现的k-ε两层湍流模型能够对非均匀流场螺旋桨空泡形态进行准确预测,计算结果与试验结果吻合,该计算方法可为非均匀流场螺旋桨空泡模拟提供参考.

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

盛振邦, 刘应中. 船舶原理[M]. 上海: 上海交通大学出版社, 2005.

SHENG

Zhenbang

, LIU

Yingzhong

. Principles of shipping[M]. Shanghai: Shanghai Jiao Tong University Press, 2005.

DOI:10.5957/jsr.2003.47.1.48 URL [本文引用: 2]

[2]

YOUNG

Y L

, KINNAS

S A

Numerical modeling of supercavitating propeller flows

[J]. Journal of Ship Research, 2003, 47(1):48-62.

[3]

NIEDZWIEDZKA

, SCHNERR

G H

, SOBIESKI

Review of numerical models of cavitating flows with the use of the homogeneous approach

[J]. Archives of Thermodynamics, 2016, 37(2):71-88.

DOI:10.1515/aoter-2016-0013 URL [本文引用: 1]

[4]

HSIAO

C T

, MA

J S

, CHAHINE

G L

Multiscale tow-phase flow modeling of sheet and cloud cavitation

[J]. International Journal of Multiphase Flow, 2017, 90:102-117.

DOI:10.1016/j.ijmultiphaseflow.2016.12.007 URL [本文引用: 1]

[5]

朱志峰, 方世良, 王晓燕.

船舶螺旋桨黏性空化流场数值方法

[J]. 东南大学学报(自然科学版), 2010, 40(6):24-29.

ZHU

Zhifeng

, FANG

Shiliang

, WANG

Xiaoyan

Numerical method for viscous capitating flow around ship propeller

[J]. Journal of Southeast University (Natural Science Edition), 2010, 40(6):24-29.

[6]

LIU

Z H

, WANG

B L

, PENG

X X

, et al.

Calculation of tip vortex cavitation flows around three-dimensional hydrofoils and propellers using a nonlinear k-ε turbulence model

[J]. Journal of Hydrodynamics, Ser. B, 2016, 28(2):227-237.

[7]

胡健, 王雅楠, 王晴, 等.

基于螺旋加密网格的螺旋桨梢涡空化数值模拟

[J]. 华中科技大学学报(自然科学版), 2020, 48(3):30-34.

[本文引用: 3]

Jian

, WANG

Yanan

, WANG

Qing

, et al.

Numerical simulation of propeller tip vortex cavitation based on helical mesh encryption

[J]. Journal of Huazhong University of Science and Technology (Nature Science Edition), 2020, 48(3):30-34.

[本文引用: 3]

[8]

刘芳远, 傅慧萍, 李杰.

螺旋桨梢涡及梢涡空泡数值模拟

[J]. 船舶力学, 2019, 23(4):388-396.

LIU

Fangyuan

, FU

Huiping

, LI

Jie

Numerical si-mulation of propeller tip vortex and TVC

[J]. Journal of Ship Mechanics, 2019, 23(4):388-396.

DOI:10.1016/j.ijmultiphaseflow.2012.02.006 URL [本文引用: 1]

[9]

, LUO

X W

, PENG

X X

, et al.

Numerical analysis of cavitation evolution and excited pressure fluctuation around a propeller in non-uniform wake

[J]. International Journal of Multiphase Flow, 2012, 43:13-21.

[10]

傅慧萍, 李杰.

斜流中的螺旋桨空化及压力脉动计算

[J]. 中国造船, 2018, 59(3):1-12.

Huiping

, LI

Jie

Calculation of propeller cavita-tion and pressure pulse in oblique flow

[J]. Shipbuilding of China, 2018, 59(3):1-12.

[11]

SCHNERR

G H

, SAUER

Physical and numerical modeling of unsteady cavitation dynamics

[C]//4th International Conference on Multiphase Flow. New Orleans: ICMF-2001, 2001: 1-12.

[12]

崔健, 伍锐, 孙硕, 等.

螺旋桨空泡观测技术研究

[J]. 水动力学研究与进展A辑, 2020, 35(5):7-18.

CUI

Jian

, WU

Rui

, SUN

Shuo

, et al.

Research of experimental technique on propeller cavitation observation

[J]. Chinese Journal of Hydrodynamics. Ser. A, 2020, 35(5):7-18.

[13]

TAO

, FREDERICK

Closure to “Discussion of ‘factors of safety for Richardson Extrapolation’”

[J]. Journal of Fluids Engineering, 2011, 133(11):1-6.