简介这份文档围绕不同工况下螺旋桨空泡的水动力性能展开面向船舶与海洋工程、水下航行器推进方向的研究生、科研人员及螺旋桨设计工程师帮助读者理解空泡产生机理、数值预报方法与参数化改善思路。包内仅含1个docx文件整体约1.57MB正文以文字与公式推导为主配合网格模型、计算域与敞水性能曲线等图示便于直接阅读与引用。文档重点梳理了国内外在网格生成、加密位置与湍流模型选择上的研究进展如结构网格与非结构网格的对比、叶梢尖端螺旋加密对梢涡空泡的捕捉并给出以MAU4-40标准桨为对象的计算流域设置、边界条件与1 046 147网格量的划分策略。数学模型部分涵盖均质混合流连续方程、动量方程、RNG k-ε湍流模型及Zwart-Gerber-Belamri空化模型还推导了临界转速与空泡初生条件讨论了纵斜角对空泡的改善作用以及AUV自航非均匀伴流场的影响。目前已有105人学习适合需要系统掌握空化数值模拟流程与工程改善思路的读者参考。1. 螺旋桨空泡在不同工况下的水动力性能先把工况这个词拆开同一只螺旋桨装在船上从压载航行切到满载低速进速系数从 0.7 掉到 0.3叶梢空泡先拉长再在前缘铺开成片空泡最后在叶背形成云状空泡。推力系数 K_T 在敞水曲线上还算平滑可一旦进到空泡工况K_T 会突然偏离扭矩系数 K_Q 偏得更狠敞水效率 η_0 直接掉几个点。拿敞水曲线去反算轴功率误差落在 5% 到 15% 区间是常事这就是为什么桨的功率预报必须自带空泡工况。不同工况下螺旋桨空泡的水动力性能研究要回答的正是这条偏离从哪来空泡数 σ 和进速系数 J 各自变化时片空泡长度、梢涡空泡轨迹、空泡体积随时间的响应是什么形态这些响应又怎么折算成 K_T、K_Q、η_0 和叶面脉动压力。做桨型设计、船舶快速性预报、空蚀与噪声评估的人需要它第一次用 CFD 碰空泡的工程师同样绕不开它。2. 螺旋桨空泡水动力计算的相似参数、空化模型与湍流模型选型2.1 四个必须锁死的相似参数J、σ_n、Re 与 Fr_n空泡算例最常见的错误是把工况设成转速固定、只改进速。真实场景里船体阻力、螺旋桨转速、桨轴浸深处的静压是联动的低速重载时转速未必下降但进速掉下来叶吸力面低于饱和蒸汽压的区域会急剧扩大。所以先用四个无量纲量把一个工况钉死再去调 CFD 的边界值。参数定义式物理含义敞水桨典型范围进速系数 JV_A / (nD)单位转速下桨前进的距离等价于载荷水平0.1 ~ 1.0空泡数 σ_n(p_0 − p_v) / (0.5 ρ n²D²)远场压力余量相对动压越小越易空化0.5 ~ 5.0雷诺数 Re_0.7Rc_0.7R · V_R / ν边界层状态与湍流尺度1×10⁵ ~ 3×10⁶傅汝德数 Fr_nn²D / g自由面兴波与吸气影响决定是否带自由面V_A 是桨盘面进速n 是转速D 是桨径p_0 取桨轴中心处的总压p_v 是当地温度下的饱和蒸汽压。注意 σ_n 用的是 nD 而不是 V_A螺旋桨空泡的驱动量是叶梢相对来流的合速度所以工程上习惯以转速为基准定义空泡数。Re 在模型尺度下往往比实船低一到两个量级若不做任何修正空泡初生点会提前这一点在把 CFD 结果和船模试验对比时要提前说清。2.2 空化模型Schnerr-Sauer、ZGB 与 Kunz 各自适合哪类空泡空化模型决定液汽两相之间的质量传递速率也是空泡算例里最容易调参调歪的地方。下面四个是绕不开的选项。模型相变速率来源主要经验系数适合的空泡类型收敛性Schnerr-Sauer气泡数密度 n 与半径 R 的解析关系n、dNuc、Cc、Cv片空泡、云空泡中等Zwart-Gerber-Belamri成核点体积分数 α_nucα_nuc、r_b、C_vap、C_cond梢涡空泡、空泡初生偏难Kunz直接给质量源不引入成核点C_dest、C_prod稳定片空泡好Merkle直接给质量源C_vap、C_cond与 Kunz 类似好// constant/phaseChangeProperties —— 两相空化相变源项配置 phaseChangeModel SchnerrSauer; SchnerrSauerCoeffs { pSat 2340; // 20 摄氏度水的饱和蒸汽压 Pa温度一变这个值必须跟着改 n 1.6e09; // 单位体积成核点数密度 1/m3调大更容易触发空化 dNuc 1.0e-06; // 成核点直径 m与 n 共同决定初始气泡半径 Cc 1.0; // 凝结系数控制空泡溃灭速率 Cv 1.0; // 蒸发系数大于 1 会让空泡长得更快 }pSat是整个文件里唯一不能靠手感凑的参数它由水温决定20 摄氏度取 2340 Pa若水温按 15 摄氏度算就要改成 1705 Pa。n和dNuc决定初始气泡半径 R₀ (3/(4πn))^(1/3) 的一半量级R₀ 越小空化越难发生。Cv与Cc才是通常用来做敏感性分析的旋钮Cv调大→空泡区体积分数偏高Cc调大→空泡尾部溃灭更急促、脉动压力偏大。做工况对比研究时一定要在同一套系数下扫 σ 和 J否则不同工况之间的差异会混进模型参数的影响。2.3 湍流模型与近壁处理SST k-ω 是基线DES 留给梢涡空泡// constant/turbulenceProperties simulationType RAS; RAS { RASModel kOmegaSST; // 空泡算例的默认基线兼顾逆压梯度与分离 turbulence on; printCoeffs on; }目标推荐模型壁面 y 要求敞水性能 K_T / K_Qk-ω SST30 ~ 100用壁面函数片空泡长度与厚度k-ω SST 局部加密前缘 1 ~ 5梢涡空泡轨迹、噪声IDDES / SST-SAS1 以下需要解析涡核大范围云空泡k-ω SST30 ~ 100SST k-ω 在逆压梯度下的分离点预测比标准 k-ε 稳得多做片空泡时前缘附近的分离线是否算准直接决定空泡长度。若关注梢涡空泡RANS 会因涡粘过大把涡核抹掉这时候才需要换成 IDDES 之类的尺度解析方法代价是网格量翻三到五倍。近壁第一层网格高度按 y 30 反推时桨叶弦长 0.05 m、来流合速度 20 m/s 的算例第一层大约取 3×10⁻⁵ m 量级。2.4 两相求解器怎么选interFoam 相变与 interPhaseChangeFoam新版 OpenFOAM 的主流做法是在 interFoam 上挂相变模型把两相设成 water 和 waterVapour相变通过phaseChangeProperties打开。老版本里对应的求解器是 interPhaseChangeFoam两者物理一致、字典组织不同。选型上记住三条带自由面就一定要用 VOF 类求解器纯空泡、无自由面可以用 interCondensatingEvaporatingFoam商业软件里对应的是 STAR-CCM 的 VOF Schnerr-Sauer 组合或 Fluent 的 mixture 多相模型 ZGB。做工况对比研究我更倾向 interFoam Schnerr-Sauer因为相变系数的物理含义相对直观工况之间的可比性好控制。3. 用 OpenFOAM 跑通螺旋桨空泡水动力算例的最小可复现流程3.1 计算域与旋转域划分MRF 起步、滑移网格收尾把计算域切成一个圆柱形旋转域和外面的静止域。旋转域直径取 1.2D ~ 1.5D包住整个桨叶并留出足够余量静止域上游留 3D、下游留 5D 以上侧向边界离桨轴 5D 以上出口压降才不会被边界反射污染。# 目录结构螺旋桨空泡算例的最小骨架 propCav/ ├── constant/ │ ├── triSurface/propeller.stl # 桨叶表面单位 m │ ├── phaseProperties # 相种类、表面张力、粘度 │ ├── phaseChangeProperties # 2.2 节的空化模型配置 │ ├── transportProperties │ ├── turbulenceProperties │ ── dynamicMeshDict # 旋转域的动网格/滑移界面设置 ├── system/ │ ├── blockMeshDict # 背景六面体网格 │ ├── snappyHexMeshDict # 桨叶表面贴合与加密 │ ├── surfaceFeatureExtractDict │ ├── controlDict │ ├── fvSchemes │ ├── fvSolution │ └── setFieldsDict # 初始化水相分布 └── 0.orig/ # 带占位符的初始场模板 ├── U p_rgh alpha.water k omega nut旋转域和静止域之间的交界面用 AMI 处理先用 MRF 跑几百步拿到一个稳定初场再切到滑移网格做非定常。空泡的一阶响应跟叶频强耦合直接用 MRF 算会平滑掉梢涡空泡的周期性脱落工况对比里这一点会被误判成这个 σ 下没有梢涡空泡。3.2 网格生成snappyHexMesh 与 yblockMesh # 生成背景网格 surfaceFeatureExtract # 提取桨叶前缘/随边特征线 snappyHexMesh -overwrite # 贴合 加密 边界层结果覆盖上层网格 checkMesh -allGeometry -allTopology # 必跑重点看 max non-orthogonality 与负体积snappyHexMeshDict里真正影响空泡结果的是refinementSurfaces与addLayersControls两组。桨叶面加密等级取 2 ~ 3 级、前缘再追加一级refinementRegions局部加密nLayer取 8 ~ 12expansionRatio取 1.2finalLayerThickness按 y 目标反算。checkMesh报出的最大非正交角超过 70 度时fvSchemes里必须补非正交修正否则压力场会出现锯齿。层覆盖失败通常发生在随边太薄的位置featureAngle从默认 60 度放到 75 度、或把minThickness调小比整体加密划算得多。3.3 边界条件与初始场远场压力按 σ_n 反算场入口出口侧向/远场桨叶面AMI 界面UfixedValue (V_A,0,0)inletOutletslipmovingWallVelocityAMIp_rghzeroGradientfixedValue p_outzeroGradientzeroGradientAMIalpha.waterfixedValue 1inletOutlet 1zeroGradientzeroGradientAMIkfixedValueinletOutletzeroGradientkqRWallFunctionAMIomegafixedValueinletOutletzeroGradientomegaWallFunctionAMInutcalculatedcalculatedcalculatednutkWallFunctioncalculatedp_out 由 σ_n 反算先把 σ_n 代进p0 p_v σ_n · 0.5ρ(nD)²得到桨轴中心处的总压再加上桨轴浸深的水柱压差换算到出口高度最后扣掉出口处的动压。这一步算错整条 σ 曲线会整体平移空泡初生点跟着偏。入口的 k 与 omega 按湍流强度 1% 和特征长度计算别直接给默认 0.24 和 0.1这类拍脑袋值。初始场用setFields把水相填到指定水位以下桨前区域可以局部抬高一点避免启动瞬间大范围空化。3.4 求解控制fvSchemes、fvSolution 与时间步// system/fvSchemes —— 空泡算例里对流项的选取 ddtSchemes { default CrankNicolson 0.9; } // 空泡响应快全隐式会让空泡滞后 gradSchemes { default Gauss linear; } divSchemes { default none; div(phi,U) Gauss limitedLinearV 1; // 速度受限线性抑制桨叶前后压力振荡 div(phi,alpha) Gauss vanLeer; // 界面锐化空泡边界才不糊 div(phirb,alpha) Gauss interfaceCompression; div(phi,k) Gauss limitedLinear 1; div(phi,omega) Gauss limitedLinear 1; } laplacianSchemes { default Gauss linear corrected; }// system/fvSolution —— PIMPLE 与 alpha 求解的关键项 alpha.water.* { nAlphaCorr 2; // 每步重算界面通量的次数云空泡可以提到 3 nAlphaSubCycles 1; cAlpha 1; // 界面压缩系数1 为默认调大界面更锐但更易振 MULESCorr yes; // 用 MULES 做有界修正保证 alpha 落在 [0,1] nLimiterIter 10; solver smoothSolver; smoother symGaussSeidel; tolerance 1e-8; relTol 0; } PIMPLE { momentumPredictor no; // 空泡算例关掉动量预测器稳定性显著提升 nOuterCorrectors 1; nCorrectors 3; nNonOrthogonalCorrectors 1; }时间步按每步转 1 度取n 25 r/s 时 Δt 1/(25×360) ≈ 1.111×10⁻⁴ s。maxCo卡在 3 ~ 5 之间maxAlphaCo卡在 1 ~ 2 之间别为了赶工把 maxAlphaCo 放到 10界面会直接糊掉。启动阶段先用固定时间步跑 500 步等 K_T 的滑动平均稳定后再开自适应步长。setFields # 初始化水相 decomposePar -force mpirun -np 32 interFoam -parallel log.interFoam 21 4. 不同工况的批量扫描从 J-σ 矩阵到 K_T、K_Q 与空泡体积响应4.1 工况矩阵怎么排先扫 σ 再扫 J工况矩阵不是越密越好。螺旋桨空泡对 σ 的敏感度远高于 J所以先用 4 个 J 值搭骨架σ 在每根骨架上取 3 ~ 5 个点把空泡初生点夹在两个相邻 σ 之间。初生点一旦定位再在它附近追加两三个 σ这样一条曲线上花掉算例最少。工况编号Jσ_n预期空泡形态用途C010.23.0无空泡 / 初生基准校 K_T 与试验敞水值C020.21.5片空泡重载工况主查点C030.20.8片空泡 云空泡效率跌落与脉动压力C040.63.0无空泡巡航工况基准C050.61.5梢涡空泡涡空泡轨迹对比C060.83.0无空泡轻载高速验证网格上限4.2 用 Python 批量生成并提交算例# batch_cav.py —— 按 J-σ 矩阵批量生成并提交 OpenFOAM 空泡算例 import itertools, subprocess, shutil from pathlib import Path D 0.25 # 桨模直径 m rho 998.2 # 20 摄氏度水密度 kg/m3 pSat 2340.0 # 饱和蒸汽压 Pa n 25.0 # 转速 r/s hShaft 0.30 # 桨轴中心浸深 m J_list [0.2, 0.4, 0.6, 0.8] sigma_list [0.8, 1.5, 3.0] base Path(templateCase) for J, sigma in itertools.product(J_list, sigma_list): case Path(fcases/J{J:.2f}_sig{sigma:.2f}) if case.exists(): shutil.rmtree(case) # 重跑时清掉旧结果避免读到污染场 shutil.copytree(base, case) Va J * n * D # 进速 m/s p0 pSat sigma * 0.5 * rho * (n * D) ** 2 # 桨轴中心处总压 Pa pOut p0 rho * 9.81 * hShaft # 折算到出口静压 Pa # 模板里 0.orig/U 用 __VA__、0.orig/p_rgh 用 __POUT__ 占位 for f, key, val in [(0.orig/U, VA, Va), (0.orig/p_rgh, POUT, pOut)]: p case / f p.write_text(p.read_text().replace(f__{key}__, f{val:.6f})) (case / caseInfo.txt).write_text( fJ{J}\nsigma{sigma}\nVa{Va:.4f}\nn{n}\nD{D}\np0{p0:.2f}\n) subprocess.run( blockMesh snappyHexMesh -overwrite setFields decomposePar -force mpirun -np 32 interFoam -parallel log.interFoam 21, shellTrue, cwdcase, checkTrue)p0与pOut的换算是整个脚本的核心sigma定义在桨轴中心而边界条件写在出口面上两者差一个 ρgh高差 0.3 m 对应约 2.9 kPa在 σ 0.8 的工况下已经占到动压的 1% 以上省掉这一项会让初生点预测偏后。itertools.product保证每个 (J, σ) 组合只出现一次不会因为列表顺序变化而重复提交。shutil.rmtree必须在复制前执行否则上一次的postProcessing目录会被保留后处理时会读到两段不同时间的力数据。4.3 从 forces.dat 提取 K_T、K_Q、η_0// system/controlDict 中追加用于记录桨叶受力 functions { forces { type forces; libs (libforces.so); patches (propeller); // 桨叶边界名与 snappyHexMeshDict 一致 rho rhoInf; rhoInf 998.2; CofR (0 0 0); // 力矩参考点取桨轴中心 writeControl timeStep; writeInterval 1; } }# kt_kq.py —— 从 forces.dat 的后两转提取 K_T、K_Q 与敞水效率 import numpy as np D, rho, n 0.25, 998.2, 25.0 t, Fx, My [], [], [] with open(postProcessing/forces/0/forces.dat) as fp: for line in fp: s line.strip() if not s or s.startswith(#): continue v s.replace((, ).replace(), ).split() t.append(float(v[0])); Fx.append(float(v[1])); My.append(float(v[4])) t, Fx, My np.array(t), np.array(Fx), np.array(My) tEnd t[-1] mask t tEnd - 2.0 / n # 只取最后两转丢掉启动瞬态 T_mean, Q_mean Fx[mask].mean(), My[mask].mean() KT T_mean / (rho * n**2 * D**4) KQ Q_mean / (rho * n**2 * D**5) info dict(l.split() for l in open(caseInfo.txt).read().strip().split(\n)) J float(info[J]) eta0 J * KT / (2 * np.pi * KQ) print(fJ{J} KT{KT:.4f} 10KQ{10*KQ:.4f} eta0{eta0:.4f}) # 加上 Fx[mask].std() 就是推力脉动的量级重载工况下它的增长比 K_T 跌落更早出现去掉启动瞬态这一步不能省。空泡算例的启动阶段通常要跑 3 ~ 5 转才进入准周期前面几转的 K_T 会虚高 10% 以上直接全时段平均会把空泡工况的效率跌落平均掉。若周期性强可以再做一次相位平均把叶频分量分离出来。4.4 空泡体积与空泡面积的定量提取// system/controlDict 中追加统计旋转域内的蒸汽相体积 vapourVolume { type volFieldValue; libs (libfieldFunctionObjects.so); fields (alpha.waterVapour); operation volIntegrate; // 对体积分数做体积分得空泡体积 m3 regionType cellZone; name rotorZone; // 限定在旋转域排除远场自由面干扰 writeFields false; writeControl timeStep; writeInterval 5; log true; }有了逐时刻的空泡体积曲线就能把水动力性能和空泡形态绑在一起看同一条 V_cav(t) 曲线上均值反映空泡总量峰值反映云空泡的瞬时爆发方差反映空泡的周期性有多强。σ 从 3.0 降到 1.5 时往往 V_cav 的均值只涨了一点、方差却翻了倍这时 K_Q 的脉动也跟着放大说明工况已经从稳态片空泡进入非稳态阶段。这个判据比单纯看 K_T 曲线敏感得多。5. 螺旋桨空泡水动力结果的网格无关性验证与发散处置5.1 网格与时间步无关性GCI 怎么算网格无关性不能靠看起来差不多了。在 σ 1.5、J 0.4 这个工况上做三套网格取基础尺寸加密比 r 1.3把 K_T 和空泡体积一起作为观测量用 Richardson 外推求表观收敛阶 p再算网格收敛指数。网格单元数前缘 yK_TV_cav (cm³)与上一级相对偏差粗320 万680.18723.42—中700 万420.18313.052.2% / 12.1%细1500 万260.18182.910.7% / 4.8%GCI Fs · |e| / (r^p − 1)安全因子 Fs 取 1.25。这里要盯的是空泡体积收敛得比 K_T 慢得多粗网格上 V_cav 偏高 17%说明空泡尾部的大尺度结构还没被解析出来。工程上的做法是让 K_T 的 GCI 控制在 2% 以内、V_cav 的 GCI 控制在 8% 以内再去做工况对比。时间步的无关性同理把 Δt 从 1 度/步减到 0.5 度/步若 K_Q 的变化小于 1%就可以认为时间离散够用。5.2 空泡算例发散的四类典型原因与对策现象常见原因处置残差在 10⁻³ 附近卡住alpha 场出现棋盘格界面压缩过强或网格非正交过大cAlpha 从 1 降到 0.5补 nNonOrthogonalCorrectors 到 2头几步就出现负 alpha初始场里桨叶附近已有大量蒸汽相用 setFields 把桨前水相填满首 200 步固定小步长压力场出现尖刺、p_rgh 出现负值跳变pSat 设置过高或 Cv 过大核对 pSat 与水温一致Cv 先回到 1.0算到十几转后突然炸掉梢涡区域网格太粗导致涡核耗散异常梢涡轨迹沿线追加 refinementRegions或换 IDDES排错顺序建议从物理量入手先看 alpha.water 的极值有没有越界再看 p_rgh 的最小值是不是远低于 pSat最后才去动离散格式。绝大多数发散都能在前两者里找到原因改格式只是掩盖问题。5.3 一个省时间的小技巧用非空泡解做空泡算例的初场空泡算例启动阶段的迭代量往往比稳态阶段多好几倍。把phaseChangeProperties里的 pSat 临时设成一个极小值比如 −1×10⁶ Pa相变就不会被触发此时求解器退化成纯两相敞水计算跑 1500 步左右把 U、k、omega 场收到稳定。然后改回真实 pSat重启继续算。这样能省掉大约三分之一的机时尤其是 σ 较高的工况几乎直接进入准周期状态。唯一要注意的是切换瞬间alpha.water场是敞水分布桨叶吸力面的低压区会在几十步内快速长出空泡这一步要盯着 V_cav 曲线确保它是单调爬升而不是振荡发散。本文还有配套的精品资源点击获取