A翼型空化笔记结构化：OpenFOAM算例复现与σ扫描验证

简介面向水动力学、船舶推进与水下航行器方向的空化专题笔记围绕 NACA0015 二维水翼的云空化及其控制机理展开适合正在做空化数值模拟、需要梳理计算设置与研究思路的研究生和工程技术人员参考。内容从计算模型入手给出弦长 c0.1m、攻角 6°~10°、弧形进口边距前缘 5 倍弦长、出口距尾缘 10 倍弦长、上下边界相距 12 倍弦长的计算域布置并整理空化数取 0.9、进口速度 12m/s、汽化压强 0.044bar 等算例参数与时间步长设定。同时延伸到 Clark-Y 翼型的二维与三维建模尺寸、NACA66(MOD) 的弯度与厚度比参数以及汽化压强与空化数判据的说明。全包仅 1 个 docx 文档约 332KB以分条记录形式组织便于按小节查阅与摘录。已有 112 人学习可帮助读者快速搭建空化算例、理解参数选择依据并对照自己的模拟结果。1. 翻开一份 A翼型空化笔记先要解决的是口径问题翻到一份 2021–2022 年的《专题资料2021-2022年A翼型空化笔记.docx》里面大概率是几十条水洞试验或数值工况的手记空化数从 3.0 一路降到 0.8云空化在哪个点开始周期性脱落升力系数掉了几个百分点某几次还夹着出口压力调错、水温没记的坑。这类资料最麻烦的地方不是结论错而是它只躺在 Word 里——工况靠肉眼翻参数散在段落和表格之间想复算某个点得先把当时的边界条件重新猜一遍。真正省时间的方向有两条把 A翼型空化的物理量纲和判据固化成统一口径再用两相求解器把每个工况跑成可复现的算例最后把 docx 里的记录批量抽成结构化参数表。量纲、算例、结构化、验证这四步走下来一份陈年笔记才有二次利用价值。适合已经能跑基础 CFD、准备系统处理水翼空化数据的人。2. A翼型空化的量纲与判据从空化数到压力系数2.1 空化数与雷诺数必须同时记录A翼型这类二维水翼剖面的空化行为本质上由两个无量纲数主导。空化数衡量「当地压力距离汽化还有多远」σ (p∞ − p_v) / (½ρU∞²)雷诺数衡量边界层状态Re U∞·c / ν两者缺一个数据就没法横向比较。常见的水洞工况是弦长 c 0.1 m、来流 U∞ 10 m/s、水温 20 °C对应 ρ 998.2 kg/m³、ν 1.002×10⁻⁶ m²/s算下来 Re ≈ 1.0×10⁶。同样的 σ 在 Re 3×10⁵ 下是层流分离主导在 Re 1×10⁶ 下是湍流边界层空化起始位置能差出 20% 弦长。这里有个容易踩的坑σ 里的 p∞ 是绝对压力而水洞台架上的压力表往往读的是表压。笔记里如果只写「入口 0.2 MPa」而不注明表压绝压复算时 σ 会整体偏移 1 个大气压对应的量也就是约 1.0 的 σ 偏差足以把云空化工况误判成无空化。2.2 σ 从高往低降A翼型上依次出现三种流态对带攻角的 A翼型吸力面最小压力系数 Cp,min 在 α 6° 时通常在 −2.0 到 −2.8 之间。理论上的空化初生条件是 σ_i −Cp,min也就是 σ 降到 2 附近才会出现第一批气泡。继续往下扫流态大致这样演化σ/σ_i 区间流态典型特征脱落频率量级 1.0无空化只有单相湍流分离—0.9 ~ 1.0空化初生前缘零散游移气泡肉眼难辨无周期0.7 ~ 0.9片状空化附着型空泡长度随 σ 单调增长弱周期0.4 ~ 0.7云状空化尾部断裂成云团周期性脱落f·c/U∞ ≈ 0.1~0.2 0.4超空化空泡越过尾缘升力塌陷准稳态云状空化段的周期性脱落是笔记里最值得记的数据因为它同时对应严重的空蚀风险和升力系数折损。经验上脱落主频落在 f·c/U∞ ≈ 0.1~0.2 区间若用空泡长度 L_cav 做参考长度Strouhal 数会升到 0.25~0.35。写笔记时把参考长度一起记上否则两个数对不上后面做 FFT 验证时会怀疑是自己算错了。2.3 一条工况记录里必须落表的 8 个量笔记里最容易缺的恰恰是换算用的中间量。下面这张表可以直接当作记录模板字段名和 4.2 节的解析正则是一一对应的量符号单位示例备注来流速度U∞m/s10.0记实测值不记设定值弦长cm0.100二维剖面按单层处理攻角α°6.0与网格固定角保持一致水温T°C20决定 p_v不能省饱和蒸汽压p_vPa2339由 T 反算不抄手册参考压力p∞kPa(abs)50.0必须标注绝压空化数σ—0.955按 2.1 节定义最小压力系数Cp,min—−2.4无空化算例给出2.4 用一段 Python 把 σ、p_v 和初生点一次算准手算 σ 最容易错在 p_v 上20 °C 用 2339 Pa、25 °C 用 3169 Pa、30 °C 用 4246 Pa而很多人整套笔记都套 2339 Pa。下面这段把饱和蒸汽压按 Antoine 方程算出来再反推 σ 和初生判据import numpy as np def p_vapour(T_celsius): Antoine 方程估算水的饱和蒸汽压适用 1~100 degC返回 Pa A, B, C 8.07131, 1730.63, 233.426 # 系数对应 mmHg 和 degC p_mmHg 10 ** (A - B / (T_celsius C)) return p_mmHg * 133.322 # mmHg - Pa def sigma(p_inf_pa, U_inf, rho998.2, T20.0): 空化数p_inf_pa 必须传绝对压力 return (p_inf_pa - p_vapour(T)) / (0.5 * rho * U_inf ** 2) def p_inf_from_sigma(sigma_target, U_inf, rho998.2, T20.0): σ 扫描时反算需要的出口绝压Pa return sigma_target * 0.5 * rho * U_inf ** 2 p_vapour(T) U, c, alpha 10.0, 0.1, 6.0 print(Re , U * c / 1.002e-6) # 约 9.98e5 print(sigma , sigma(50e3, U)) # 约 0.955 print(p_inf sigma2.0 , p_inf_from_sigma(2.0, U) / 1e3, kPa)逻辑说明p_vapour用 Antoine 方程拟合20 °C 时返回约 2332 Pa与手册值 2339 Pa 相差 0.3%对 σ 的影响在千分之一量级完全可以接受。sigma强制要求传绝对压力参数名是为了在代码层面堵住表压/绝压混用。p_inf_from_sigma是后面第 5 章 σ 扫描的入口改 σ 只改一个参数U 和 T 保持不动这样每次变化的只有压力边界排错范围最小。参数说明ρ 默认取 20 °C 淡水密度 998.2 kg/m³做海水工况要改成 1025ν 取 1.002×10⁻⁶ m²/sT在 σ 计算里只影响 p_v但会通过影响 ρ 和 ν 间接改变 Re严谨做法是把三者都按实测水温查表。3. 用 OpenFOAM 跑通 A翼型空化的最小算例3.1 从翼型坐标到 C 型网格blockMesh 与 snappyHexMesh 的分工A翼型坐标到手通常是一份两列的 .dat 点集前缘在 (0 0 0)、弦向沿 x。别指望直接喂给 blockMeshblockMesh 只负责造背景六面体块翼型附近的贴体网格交给 snappyHexMesh。目录结构按 OpenFOAM 惯例铺开# 算例目录以 OpenFOAM v9 为例 A-foil-sigma1.0/ ├── 0/ U p alpha.water k omega nut ├── constant/ transportProperties cavitationProperties turbulenceProperties │ triSurface/A-foil.dat └── system/ blockMeshDict snappyHexMeshDict controlDict fvSchemes fvSolution背景域取半径 10c 的 C 型块弦向 400 个网格、法向 120 个前缘下游 1.5c 范围内再单独加密。执行顺序固定成下面四条日志单独落盘方便回头查是哪一步出的问题cd A-foil-sigma1.0 blockMesh log.blockMesh 21 surfaceFeatureExtract log.sFE 21 snappyHexMesh -overwrite log.sHM 21 checkMesh -allGeometry -allTopology | tee log.checkMesh-overwrite让 snappyHexMesh 直接覆盖 0/ 之外的网格目录省掉手工搬 polyMeshcheckMesh加-allGeometry会检查面法向和凹面空化算例里负体积网格经常表现为 α 场在某个角落突然冲到 1 又不消退根因基本都在这里。壁面第一层网格按 y 目标反算。u_τ ≈ U∞·√(C_f/2)C_f 用 0.058·Re^−0.2 估Re 1.0×10⁶ 时 u_τ ≈ 0.43 m/s。要 y ≈ 1 就得把首层厚度压到 2.3×10⁻⁶ m要 y ≈ 30 则是 7×10⁻ m。空化算例优先选前者配低雷诺数湍流模型因为空泡脱落位置对边界层分离点极其敏感壁面函数会把分离点往下游推。3.2 constant/ 下两个必改的文件相变模型参数集中在cavitationProperties这是整套算例里最需要按自己工况调的地方// constant/cavitationProperties字段名以本机 tutorials 为准 cavitationModel SchnerrSauer; SchnerrSauerCoeffs { pSat 2339; // 20 degC 水的饱和蒸汽压 [Pa]随水温改 n 1.6e13; // 单位液体体积内气核数 [1/m3] dNuc 2.0e-06; // 气核直径 [m] Cc 1.0; // 凝结系数 Cv 1.0; // 蒸发系数 }参数说明pSat必须和 2.4 节算出的 p_v 一致两处对不上会导致空化区大小系统性偏差n和dNuc决定初始气核含量水洞自来水取默认量级即可若用去气水做试验把 n 降到 1e12 更接近实际Cc和Cv是经验系数默认 1.0 时云空化脱落偏「黏」想加快空泡溃灭可把 Cc 调到 10~100但调过 100 容易让 α 场出现非物理振荡。transportProperties里定义两相物性注意 ν 和 ρ 的单位分别是 m²/s 和 kg/m³气相密度按 0.023 量级给别照抄空气的 1.2// constant/transportProperties phases (water vapour); water { transportModel Newtonian; nu 1.002e-06; rho 998.2; } vapour { transportModel Newtonian; nu 2.2e-05; rho 0.023; }3.3 system/ 里的时间步与界面捕捉时间步用adjustTimeStep自动控制但上限要卡住。弦向网格 400 等分时最小网格约 2.5×10⁻⁴ m取 Co 0.5 得到 dt ≈ 1.25×10⁻ s。云空化脱落主频按 f·c/U∞ 0.2 估U 10 m/s、c 0.1 m 对应 f ≈ 20 Hz一个脱落周期 0.05 s要抓频谱至少跑 10 个周期也就是endTime给到 0.5 s、约 4 万步。// system/controlDict 关键项 adjustTimeStep yes; maxCo 0.5; // 全局库朗数上限 maxAlphaCo 0.5; // 界面附近单独限制比 maxCo 更关键 deltaT 1e-05; // 初始步长 endTime 0.5; // 至少 10 个脱落周期 writeInterval 0.0025; // 每 5 个周期采一帧体积分数云图够用// system/fvSchemes 界面项别用一阶迎风 div(phi,alpha) Gauss vanLeer; div(phirb,alpha) Gauss interfaceCompression;// system/fvSolution 的 alpha.water 段 alpha.water.* { nAlphaCorr 2; // 每步界面代数修正次数 nAlphaSubCycles 1; cAlpha 1; // 界面压缩系数2 易造成界面破碎 MULESCorr yes; // 半隐式 MULES比显式稳定 nLimiterIter 10; solver smoothSolver; smoother symGaussSeidel; tolerance 1e-8; relTol 0; }参数说明maxAlphaCo单独限制界面附近的库朗数是云空化能不能稳定脱落的关键开关只调maxCo而放着它不管界面会在脱落瞬间被抹平。cAlpha控制界面锐度1 到 2 之间比较稳妥超过 2 会看到零散的气相斑点在尾迹里乱窜那是数值扩散造成的假空泡。nAlphaCorr给 2 一般够残差在界面处不降时可加到 3但计算量线性上升。3.4 首轮运行与结果可信度自查求解器用interPhaseChangeFoam前台跑一遍同时把日志留住interPhaseChangeFoam log.run 21 tail -f log.run | grep -E Time |alpha.water判断这一步结果的可靠性有三个即时可查的信号残差方面p 和 U 的初始残差应降到 1e-5 以下并保持平稳α 的残差在空化区本就不会太小看到 1e-3 量级不算异常网格方面checkMesh的Max aspect ratio控制在 200 以内C 型网格尾迹区拉长是常态但超过 500 会让脱落频率偏低场量方面吸力面前缘的 α_water 应在 0.2 到 0.9 之间连续过渡如果出现整片 0 或整片 1 的硬边界说明 cAlpha 或网格尺度至少有一个不对。提示不同 OpenFOAM 小版本里相变系数的文件名和字段名会有出入v9 之后部分字段从transportProperties挪到了cavitationProperties。最省事的确认方式是打开安装目录下tutorials/multiphase/interPhaseChangeFoam里任一算例对着自己的字典逐项核一遍比翻文档快。4. 把 2021–2022 的 docx 笔记转成可查的参数表4.1 先定字段再写正则docx 里的工况记录通常是自然语言句子比如「σ 1.20U 10 m/sT 20 ℃α 6°p∞ 45 kPa云空化周期性脱落明显」。要抽的是六类信息空化数、速度、温度、攻角、压力、流态描述。把正则按字段集中管理后面加字段只改字典PATTERNS { sigma: r[σs]\s*[:]\s*(-?\d(?:\.\d)?), U: rU\s*(?:∞|inf)?\s*[:]\s*(\d(?:\.\d)?)\s*m/s, T: rT\s*[:]\s*(\d(?:\.\d)?)\s*(?:°C|℃|C), alpha: r[αa]\s*[:]\s*(-?\d(?:\.\d)?)\s*(?:°|deg), p: rp\s*(?:∞|inf)?\s*[:]\s*(\d(?:\.\d)?)\s*kPa, state: r(无空化|空化初生|片状空化|云状空化|超空化), }用[σs]而不是直接写 σ是因为不少笔记会用字母 s 代替希腊字母state限定枚举值避免把「空化数」三个字误当成流态。这套写法对固定格式的手记命中率很高但对一句里塞两个工况的长句会失效那种情况只能拆句后再匹配。4.2 python-docx 按文档流抽段落和表格python-docx 的doc.paragraphs只返回段落表格得单独取两者顺序会丢。按 XML 子元素遍历才能还原真实阅读顺序下面这段是通用写法import re import pandas as pd from docx import Document from docx.table import Table from docx.text.paragraph import Paragraph def iter_blocks(doc): 按文档流顺序产出段落和表格 for child in doc.element.body.iterchildren(): tag child.tag.rsplit(}, 1)[-1] if tag p: yield Paragraph(child, doc) elif tag tbl: yield Table(child, doc) def parse_notes(path): doc Document(path) rows, table_idx [], 0 for block in iter_blocks(doc): if isinstance(block, Table): table_idx 1 for r in block.rows: cells [c.text.strip() for c in r.cells] if any(k in .join(cells) for k in (σ, 空化数)): rows.append({raw: | .join(cells), source: ftable{table_idx}}) continue text block.text.strip() if not text or (σ not in text and 空化数 not in text): continue row {raw: text, source: block.style.name} for key, pat in PATTERNS.items(): m re.search(pat, text) row[key] m.group(1) if m else None rows.append(row) return pd.DataFrame(rows) df parse_notes(专题资料2021-2022年A翼型空化笔记.docx) print(len(df), df[state].value_counts(dropnaFalse))逻辑说明iter_blocks直接遍历body的子节点并按命名空间后缀区分w:p和w:tbl这样段落和表格的先后关系被保留source字段能告诉你这条记录来自正文还是第几张表。表格里的行也走同一套正则是偷懒但有效的做法因为表格单元格里的文字格式和正文几乎一样。df[state].value_counts一眼就能看出 2021 和 2022 两年记的工况偏向哪种流态这往往是整理笔记时最先发现的规律。参数说明(σ not in text and 空化数 not in text)是预过滤把没有工况信息的段落直接跳过几十页的文档能在几秒内跑完如果笔记里用「sigma」英文写法需要在预过滤里补上。抽取完的 DataFrame 存成 parquet 比 csv 稳妥中文字段和空值都不会出问题。4.3 交叉校验把可疑记录挑出来抽出来的数字未必可信最有效的校验是拿 σ 反算压力——同一个工况下 σ、p∞、U、T 四个量是锁死的任意三个能推出第四个。偏离超过阈值就列为待核import numpy as np def qc(df, rho998.2, tol0.05): d df.dropna(subset[sigma, U, T, p]).copy() d[p_v] d[T].astype(float).map(p_vapour) d[sigma_check] ((d[p].astype(float) * 1e3 - d[p_v]) / (0.5 * rho * d[U].astype(float) ** 2)) d[dev] (d[sigma_check] - d[sigma].astype(float)).abs() d[qc] np.where(d[dev] tol, 核对压力单位/水温, OK) return d.sort_values(dev, ascendingFalse)判定规则和处置动作对应关系如下偏离量最可能原因处置恒定偏移约 1.0·(½ρU²)/... 量级p∞ 记成了表压加 101.3 kPa 后重算偏离随温度变化整份笔记套了同一个 p_v按各行 T 重算 p_v单条偏离极大小数点或单位抄错回原 docx 逐字核全部 ≤ 0.05记录自洽直接进算例生成流程注意docx 里嵌的截图、粘贴的 Excel 对象 python-docx 一律读不到doc.inline_shapes只能数出个数。真要把图也归档用zipfile打开 docx 列出word/media/下的文件按文档顺序重命名后单独存一份别指望文本解析能把图里的曲线还原成数据。5. σ 扫描与空化初生判定把笔记里的结论验证一遍5.1 六个工况点的批量生成复算的核心是只变一个量固定 U 10 m/s、T 20 °C、α 6°只改出口绝压。σ 取 3.0、2.5、2.0、1.5、1.0、0.8前半段用来定初生点后半段用来抓云空化脱落import shutil, os, re base, U, T A-foil-base, 10.0, 20.0 for s in [3.0, 2.5, 2.0, 1.5, 1.0, 0.8]: dst fA-foil-sigma{s:.1f} shutil.copytree(base, dst) p_out p_inf_from_sigma(s, U, TT) # 2.4 节的函数 f os.path.join(dst, 0, p) txt open(f).read() # 出口 fixedValue 段替换压力值其余边界保持 txt re.sub(r(outlet\s*\{[^}]*?value\suniform\s)[\d.eE-], rf\g1{p_out:.1f}, txt) open(f, w).write(txt) print(dst, p_outlet , round(p_out / 1e3, 2), kPa)逻辑说明copytree保证六个算例的网格、模型参数、时间步完全一致唯一差异是压力出口值任何结果差别都能归因到 σ 上。正则只替换outlet块里的value uniform不碰 inlet 的zeroGradient避免把速度入口误伤。六个算例可以并行提交用mpirun -np 4各跑各的互不干扰。5.2 三条判据一起用别只看云图只看体积分数云图判断「有没有空化」主观性太强。稳妥做法是三条判据同时满足才下结论一是时间平均升力系数相对无空化算例下降超过 5%这条对网格不敏感、可比性最好二是 α_vapour 最大值超过 0.5 且从第 3 个周期起连续出现说明空泡不是数值噪声三是压力监测点的频谱在预期频段出峰。第三条最能证伪用下面几行就能算出来import numpy as np c, U 0.1, 10.0 raw np.loadtxt(postProcessing/probe/0/p, skiprows1) # 列: t, p n0 int(0.3 * len(raw)) # 丢掉启动瞬态 t, p raw[n0:, 0], raw[n0:, 1] dt float(np.mean(np.diff(t))) freq np.fft.rfftfreq(len(p), dt) amp np.abs(np.fft.rfft(p - p.mean())) k int(np.argmax(amp[1:]) 1) print(f主频 {freq[k]:.2f} Hz, St_c {freq[k] * c / U:.3f})skiprows1对应 OpenFOAM 监测文件的首行表头丢掉前 30% 是为了剔除启动瞬态否则低频漂移会盖过真实的脱落峰。算出的 St_c 落在 0.1~0.2 之间说明抓到的确实是云空化脱落如果只有孤零零一个低频峰多半是出口压力反射造成的伪周期把出口位置再往下游推 5c 重跑即可。排错时按这个顺序查云空化完全不脱落先看maxAlphaCo是不是超了 0.5再看弦向网格是否不足 300空泡长度明显偏短先核对pSat和PSat是否一致再查dNuc是否被误设成 1e-3 量级六个 σ 算例结果几乎一样八成是0/p里的出口值没替换成功直接去算例目录里grep uniform 0/p确认。两年笔记里那些「不脱落」的记录相当一部分是时间步太大把界面抹平了值得按这套流程重跑一遍再下结论。本文还有配套的精品资源点击获取