OPENFOAM / MODEL AUDIT / SOURCE-FIRST

wallBoiling：从 C++ 源码
重新理解壁面沸腾与 CHT

基于 OpenFOAM Foundation v14：逐项推导 RPI 热流、气泡闭式、质量源项、温度壁面迭代与跨区域耦合，并解释此前示意中可能误写的公式。

2026-10-10源码核验单文件 · 离线可读内置 SVG · 交互公式
00 · 一页读懂

先说结论：源码并不是简单的“三项热流相加”

OpenFOAM Foundation OpenFOAM-14 中的 multiRegion/CHT/wallBoiling，在同一案例中组合了固体导热、Euler–Euler 两流体、壁面核态沸腾、流体内部气液相变。模型沿用 RPI 思想，但实现层面存在面积系数的截断、壁面温度迭代、体积化的质量源项，以及通过湍流热扩散率边界条件回写热流等步骤。

研究对象

R12 竖直管流沸腾

轴对称楔形网格，管壁外侧施加热流，液体在管内向上流动；壁面生成蒸汽，流体内部还可继续蒸发或冷凝。

重要区别

两类相变不能混写

wallBoiling 计算壁面气泡生成；heatTransferLimitedPhaseChange 计算气液界面热量失衡产生的体积相变。它们有各自独立的质量源项。

q_{\mathrm{RPI}}^{\prime\prime}=q_{\mathrm{conv}}^{\prime\prime}+q_{\mathrm{quench}}^{\prime\prime}+q_{\mathrm{evap}}^{\prime\prime}
RPI 的概念性热流分解；下文给出 OpenFOAM-14 实际使用的 A₁、A₂ 与 A₂E。

核验基准：壁面模型核心代码 ↗ · 教程 fvModels ↗ · 区域求解器配置 ↗

01 · 物理边界

DEBORA 算例：固体外壁输入热流，而不是液体壁面直接给热流

模型针对 R12 的竖直管内亚冷流动沸腾。实际字典首先构造半径为 9.6 mm 的流体区域及 1 mm 厚的固体壁，再经 extrudeMesh 创建楔形网格、splitMeshRegions 分成 fluid 与 solid。

轴向剖面（原案例为绕 x 轴的 wedge 几何；图示旋转仅为阅读方便） 外壁均匀输入热流q″outer = 66,919.2 W/m²固体：k = 16.3 W/(m·K) 壁面核态沸腾成核 Nₐ / 脱离直径 d / 频率 f 体积气液相变Hₗ(Tₗ−Tₛₐₜ)+Hᵥ(Tᵥ−Tₛₐₜ) 流体侧双温度场T.liquid / T.gas（平均场） 入口液体 → 轴向向上流动 流体：Rin = 9.6 mm壁厚 1 mm出口 ↑几何加热段长度 3.5 m
物理示意非按比例绘制。气泡是模型中的统计平均量，不是被网格分辨的真实气液界面。
参数当前案例字典或物理含义
工质R12；摩尔质量 120.914气、液使用物性表
加热段长度3.5 mblockMeshDict
流体半径 / 固体外半径9.6 mm / 10.6 mm管壁厚 1 mm
入口液体温度341.67 KT.liquid
初始压力2.62 MPap
初始液相轴向速度1.75175 m/sU.liquid
固体外壁热流66,919.2 W/m²externalTemperature
固体导热率16.3 W/(m·K)solid/physicalProperties
q^{\prime\prime}_{\mathrm{out}}=q^{\prime\prime}_{\mathrm{in}}\,\frac{R_{\mathrm{in}}}{R_{\mathrm{out}}}
等热功率、等轴向长度下的圆柱壁面积换算，忽略端面传热。

代入 Rin = 0.0096 m、Rout = 0.0106 m、q″inner = 73,890 W/m²，恰得到外壁 66,919.2 W/m²。这个换算解释了为何固体温度边界与验证脚本中的热流数值不同。几何网格字典 ↗ · 外壁施热字典 ↗

02 · 软件架构

双区域 + 双相 + 两套相变：实际的数据依赖关系

foamMultiRun fluid → multiphaseEulersolid → solid 两流体方程αᵢ、Uᵢ、Tᵢ、p / 湍流与相间力 wallBoiling（wall patch）q″RPI、ṁw、气泡统计量；αₜ 边界 heatTransferLimitedPhaseChange气液界面热传递限制的体积相变 ṁlv 固体导热方程ρs、Cvs、κs、T 外壁：externalTemperature规定 q″ = 66,919.2 W/m² CHT界面 mappedWall + 两侧温度耦合：multiphaseCoupledTemperature ↔ coupledTemperature
这是求解模块之间的物理与数据依赖图，不代表内部唯一的函数调用顺序。
system/controlDict · 原案例选取的求解模块
regionSolvers
{
    fluid           multiphaseEuler;
    solid           solid;
}

两个区域的耦合不能被简写为“固体的 T 等于液体的 T”。solid/T 的 wall_inner 使用 coupledTemperature 且 Tnbr T.liquid；流体侧 T.liquid 与 T.gas 的 wall 均使用 multiphaseCoupledTemperature。实际边界函数以各相体积分数和有效导热率形成热导权重，并通过 mapped patch 交换邻域信息。

多相流体侧边界源码 ↗ · 固体侧温度边界源码 ↗ · 液体 T 字典 ↗

q_{\mathrm{solid}\to\mathrm{fluid}}^{\prime\prime}\approx-\kappa_s\,\nabla T_s\cdot\boldsymbol{n}\quad\mathrm{(interface\ flux\ balance)}
这是用于解释热量方向的连续体近似关系；实际界面采用 mixed 耦合边界离散，不是简单强制相温度相等。
03 · 源码级数学

RPI 热流分解：公式逐项核验

这部分以 wallBoiling.C::calcBoiling() 为准。代码的返回值是当前液体侧沸腾面所计算的总热流密度。尤其需要注意 A1、A2 与 A2E 不是同一个面积系数。

3.1 壁面润湿比例 Fw：Lavieville

partitioningModel_->wetFraction(lagProps.alphaLiquid) 用液相体积分数 α_l 计算液体润湿比例（并非蒸发面积率）。对于本教程 alphaCrit = 0.2：

F_{w}=1-0.5\exp[-20(\alpha_l-\alpha_{\mathrm{crit}})]\quad(\alpha_l\geq\alpha_{\mathrm{crit}})
αₗ ≥ αcrit 的分支。
F_{w}=0.5(\alpha_l/\alpha_{\mathrm{crit}})^{20\alpha_{\mathrm{crit}}}\quad(\alpha_l<\alpha_{\mathrm{crit}})
αₗ < αcrit 的分支，当前指数 20αcrit = 4。

在 α_l = 0.2 时 F_w = 0.5；在 α_l → 0 时润湿比例趋零；在 α_l → 1 时接近一。Lavieville::calculate() ↗

3.2 成核点密度 Na：Lemmert–Chawla

N_{a}=C_{n}N_{\mathrm{ref}}\left[\max\left(\frac{T_{w}-T_{\mathrm{sat}}}{\Delta T_{\mathrm{ref}}},0\right)\right]^{1.805}
Nₐ 单位为 m⁻²；仅在壁面温度高于当地饱和温度时得到正的成核点密度。

当前教程采用 Cn=1、NRef=3×10⁷ m⁻²、ΔTRef=10 K。源码中正指数严格为 1.805，不能用线性超热假设替代。LemmertChawla::calculate() ↗

3.3 脱离直径 ddep：Tolubinski–Kostanchuk

d_{\mathrm{dep}}=\max\left[d_{\min},\min\left(d_{\mathrm{ref}}\exp\left(\frac{T_l-T_{\mathrm{sat}}}{45\,\mathrm{K}}\right),d_{\max}\right)\right]
与局部液体过冷度相关，且受上、下界截断。

教程配置为 dRef=0.00024 m、dMax=0.0014 m、dMin=1e−6 m。特别注意公式使用的是 T_l 与 T_sat，而不是直接把固体温度当作液体温度。TolubinskiKostanchuk::calculate() ↗

3.4 脱离频率 fdep：Kocamustafaogullari–Ishii

f_{\mathrm{dep}}=\frac{C_f}{d_{\mathrm{dep}}}\left[\frac{\sigma|g|(\rho_l-\rho_v)}{\rho_l^2}\right]^{1/4}
源码使用四次方根，重力大小 |g|、表面张力 σ、气液密度均参与计算。

当前 Cf=1.18。这是一个量纲可闭合的频率关联式，不能写成简单的 sqrt(g/d)。频率模型 calculate() ↗

3.5 过冷 Jakob 数、影响面积与截断

Ja=\frac{\rho_l C_{p,l}\,\max(T_{\mathrm{sat}}-T_l,0)}{\rho_v L}
仅计入过冷部分；若 Tₗ ≥ Tₛₐₜ，则 max(...) = 0。
A_l=F_w\,4.8\exp\left[\min\left(-\frac{Ja}{80},\ln(\mathrm{vGreat})\right)\right]
源码使用了 4.8、80 以及 exp 上界保护；通常可理解为 Fw×4.8×exp(−Ja/80)。
X=\frac{\pi d_{\mathrm{dep}}^2 N_a A_l}{4},\quad A_2=\min(X,1),\quad A_1=\max(1-A_2,10^{-4})
A₂ 被截断到 1，A₁ 下界设为 10⁻⁴。
A_{2E}=\min(X,5)\quad\mathrm{(not\ a\ bounded\ area\ fraction)}
A₂E 截断上界是 5，不是 1；因此它不能被称为有界的几何表面积率。

必须区分三个系数。 A1 用在对流热流项；A2 用在淬冷热流项；A2E 用在气泡蒸发质量源项。将三者都替换为 A2，即便整体公式看上去符合传统 RPI 结构，计算结果也会和当前 OpenFOAM-14 源码不一致。

A1 / A2 / A2E 的实现 ↗

3.6 壁面产生的体积化质量源项

\dot m_w=\frac{1}{6}\,A_{2E}\,d_{\mathrm{dep}}\,\rho_v\,f_{\mathrm{dep}}\,\frac{S_f}{V_c}
S₍f₎ 为 wall patch 面积；V₍c₎ 为对应近壁单元体积。该表达式给出 kg/(m³·s)。
[\dot m_w]=[\rho_v][d][f][S_f/V_c]=\mathrm{kg\,m^{-3}\,s^{-1}}
量纲核验：不是 kg/(m²·s) 的面积通量。

源码直接写作 mDot = (1.0/6.0)*A2E*dDeparture*rhoVapour*fDeparture*AbyV。其中 AbyV=magSf/V。这也意味着对计算网格局部面积 / 体积的使用应单独检查。wallBoiling.C 第 482 行 ↗

3.7 淬冷换热系数与淬冷热流

a_l=\frac{\kappa_l}{\rho_l C_{p,l}}
这里的 aₗ 是热扩散率；不要与液相体积分数 αₗ 混淆。
h_Q=2\,\kappa_l\,f_{\mathrm{dep}}\sqrt{\frac{\eta_w}{\pi a_l f_{\mathrm{dep}}}}
ηw = bubbleWaitingTimeRatio，默认 0.8；本式为源码表达式在 fdep > 0 时的整理形式，源码还包含 max(fdep,small) 数值保护。
q_{\mathrm{quench}}^{\prime\prime}=A_2 h_Q\max(T_w-T_l,0)
淬冷发生在 A₂ 所表示的相应面积上；温差为当前壁面温度与估计液体温度之差。

因此，qQuenching 不能简单写成 hQ(Tw−Tc)。源码采用的是T_l，其默认可由热壁函数估算到 y+≈250 的相应液体温度。calcBoiling 淬冷热流段 ↗

3.8 蒸发热流与对流热流

q_{\mathrm{evap}}^{\prime\prime}=\dot m_w L\frac{V_c}{S_f}
蒸发项利用同一个 mDot 换算回 W/m²；体积/面积因子与质量源项中相反。
q_{\mathrm{conv}}^{\prime\prime}=A_1\,\alpha_{t,\mathrm{conv},l}\,C_{p,l}\,\delta\,\max(T_w-T_{c,l},\varepsilon T_{c,l})
δ 为壁面法向的 deltaCoeffs；ε 表示代码使用的 small 正则项。αₜ,conv,l 为 Jayatilleke 热壁函数计算的动态湍流热扩散率。

上述三部分相加，才是 calcBoiling() 返回的沸腾热流密度。这里 qconv 的量纲严格为 W/m²，因为 αt 的量纲是 kg/(m·s)，与 Cp、∇T 相乘后是热流密度。返回总热流的源码 ↗

04 · 数值实现

壁面温度不是随手代入：OpenFOAM 用二分法闭合壁面热流

当前案例液体侧为 multiphaseCoupledTemperature，在 C++ 中它属于 mixed 类型温度边界。沸腾源模型读取 mixed 边界所代表的热量约束，构造关于未知壁面液体温度 T_w 的残差函数，再用区间二分法求解。

R(T_w)=q_{\mathrm{RPI}}^{\prime\prime}(T_w)-\mathcal{Q}+hT_w=0
仅用抽象符号 ℚ 与 h 表示混合温度边界贡献：源码形式为 calcBoiling(T) − TLiquidHTaPlusQa + TLiquidH*T。
wallBoiling.C · 二分法逻辑（保留关键语句）
auto R = [&](const scalarField& T)
{
    return calcBoiling(mDot, lagProps, T)
         - TLiquidHTaPlusQa
         + TLiquidH*T;
};

mDot.boiling_ = neg(R(lagProps.Tsat));
scalarField T0(lagProps.Tsat);
scalarField T1
(
    max(TLiquid + (TLiquid - lagProps.Tsat),
        lagProps.Tsat*(1 + sqrt(tolerance_)))
);
// interval bisection on T0/T1 follows...

随后源码将求出的整体壁面热流回写为液体侧等效湍流热扩散率：这不是增加一个新的固体导热边界条件，而是让能量方程离散通量体现当前 RPI 沸腾热流。

\alpha_{t,\mathrm{boil},l}=\frac{q_{\mathrm{RPI}}^{\prime\prime}}{C_{p,l}\,\delta\,\max(T_w-T_{c,l},\varepsilon T_{c,l})\,\max(\alpha_l,\varepsilon_\alpha)}
实际代码还用 mDot.boiling_ 在沸腾区与普通对流区之间切换；此式展示沸腾区的等效量。

气相也有对应的热扩散率回写：气侧使用液体润湿比例的余部分，这是另一个不能将“总沌腾热流”直接全部给液体的原因。

\alpha_{t,\mathrm{boil},v}=\frac{1-F_w}{\max(1-\alpha_l,\varepsilon_\alpha)}\,\alpha_{t,\mathrm{conv},v}
该关系对应 wallBoiling.C 中的 alphatVapour_ 更新。

还存在一处容易忽略的参数差异：本模型的 wallBoiling::Prt_ 默认是 0.85；而本案例 thermophysicalTransport.liquid 中 RAS/Prt 是 1。两者进入不同的子模型或通量计算路径，不应未经核对就视为同一个全局常数。读取内部缺省参数 ↗

05 · 第二种相变

壁面成汽 ≠ 气液界面蒸发冷凝：相间模型要单独列出

heatTransferLimitedPhaseChange 假定流体内部气液相界面处满足局部饱和状态，以 Tsat(p) 确定界面温度，再使用气液两侧的传热系数计算界面热量失衡。核心源码如下：

\dot m_{lv}=\frac{H_l(T_l-T_{\mathrm{sat}})+H_v(T_v-T_{\mathrm{sat}})}{L}
简化表达忽略模型中显式使用的松弛：源码对 mDot 做旧值与新值的线性加权。
heatTransferLimitedPhaseChange.C · correctMDot()
const Pair<tmp<volScalarField>> Hs =
    solver_.heatTransfer.Hs(phase1_, phase2_, scalar(0));
const volScalarField::Internal& H1 = Hs.first();
const volScalarField::Internal& H2 = Hs.second();

mDot_ = (1 - f)*mDot_
      + f*(H1*(T1 - Tsat) + H2*(T2 - Tsat))/L;

当前 heatTransfer 字典在气泡分散于液体时，对液体侧使用 RanzMarshall，气侧使用 spherical；两个换热系数的构造不是由 RPI 的成核点密度直接给出的。

壁面源项

wallBoiling:mDot

位于沸腾 wall patch 的相邻单元；依赖 A2E、气泡直径、频率以及 Sf/Vc。

相间源项

heatTransferLimitedPhaseChange:mDot

在气液共存区域，由两侧温度、饱和温度、传热系数与潜热决定，可对应蒸发或冷凝。

两个源项的物理位置与闭合关系不同。不能把壁面 q_evap 与相间 mDot×L 不加区分地相加，作为一个统一的“壁面总热流”；否则容易造成错误的能量统计或重复计账。正确的全局守恒检查应基于各区域能量方程、面通量、相变源项的实际定义。

相间模型计算式 ↗ · 教程相变模型配置 ↗

06 · 交互验证

把两个闭式关系画出来：不用实际 CFD 计算也能检查趋势

下方曲线由当前教程的实际系数与 C++ 对应数学表达式在浏览器实时计算。这里只显示两个单独的经验关联式，不等于预测完整壁面热通量。

润湿比例 Fw(αl)

液体体积分数 αl = 0.80

计算值 Fw = 1.0
临界体积分数 0.2；公式见 Lavieville.C

脱离直径 ddep(ΔTsub)

过冷度 Tsat−Tl = 20.0 K

计算值 ddep = 0.15 mm
仅使用教程中的 dRef、dMax、dMin。

曲线为当前 OpenFOAM-14 子模型的静态公式，而非运行算例结果
07 · 更正记录

哪些常见简化会导致公式看着对、源码却对不上？

易误写的描述OpenFOAM-14 的实际行为
“RPI 三项热流都使用相同的气泡影响面积比例。”对流 A1、淬冷 A2、蒸发质量源项 A2E，上限分别不一样。
“wetFraction 是蒸发面积比例。”它是根据液体体积分数计算的润湿因子，先参与 Al，再进入面积系数。
“壁面蒸发率直接是 kg/(m²·s)。”当前 mDot_ 是 near-wall cell 中的体积质量源项，单位 kg/(m³·s)，带 Sf/Vc。
“直接用 Tw−Tcell 算所有换热项。”淬冷用 Tw−Tl；对流梯度用 Tw−TcLiquid，且 Tl 默认受液体温度壁函数影响。
“wallBoiling 自己把壁面热量传给固体。”固液区域热通量由 coupledTemperature / multiphaseCoupledTemperature 实现映射耦合；沸腾模型回写液体等效 alphat。
“当前 RPI 沸腾属于解析 VOF 气泡。”不是。它在 Euler–Euler 平均场框架下借助亚网格经验关联式表示气泡生成。
“所有蒸发冷凝都在 wallBoiling 里。”案例另有独立的 heatTransferLimitedPhaseChange，允许气液内部发生相间相变。
“各版本 wall boiling 边界名称一致。”OpenFOAM-14 使用 alphatPhaseChangeWallFunction；不能直接照抄其他分支的 alphatWallBoilingWallFunction。
08 · 前世今生

本案例的历史，以及为什么 v14 的公式应优先看 v14 源码

1991 · RPI / Kurul–Podowski

提出多维沸腾流动中的壁面热流分配框架。

2012 / 2019 · VTT 研究与模型验证

Peltola 等人的壁面沸腾模型研究构成 OpenFOAM 实现的重要文献基础。

2022-11-16 · 现有 CHT 案例的直接起点

VTT 的 Juho Peltola 贡献 CHT 版本，取代单区域 Euler–Euler 的旧 wallBoiling 案例。历史提交 ↗

2023 · 模块化多区域求解器

教程转入 multiRegion/CHT 目录，采用模块化求解与 foamMultiRun。

2025-05-21 · 相变模型重构为 fvModels

壁面沸腾与相间相变模型重新组织。历史提交 ↗

2026-02-05 · 壁面相变热扩散边界统一

沸腾与壁面冷凝共用 alphatPhaseChangeWallFunction；仅允许不同模型在非重叠的活跃面上生效。历史提交 ↗

文献关联：Kurul & Podowski (1991)；Peltola & Pättikangas (2012)；Peltola et al. (2019)。本博客的公式优先采用 OpenFOAM Foundation v14 中的 C++ 实现而非不同版本的二次转述。

09 · 实操路径

若要复现实验或迁移到新的相变工质，推荐怎样检查？

建议先原样运行官方案例并保存质量、能量平衡结果，再改工质。下面是当前 Allrun 的主要步骤，按原教程执行并保留验证过程：

当前教程核心运行流程
blockMesh
extrudeMesh
splitMeshRegions -cellZones all
decomposePar -allRegions
foamMultiRun -parallel       # 由 Allrun 中的 runParallel 包装执行
reconstructPar -latestTime -allRegions
# 之后执行 foamPostProcess 和 validation 脚本

迁移到氨热管时，务必重新建立氨的饱和曲线、液/汽物性、表面张力、流动条件与壁面沸腾经验关联式适用性。热管的蒸发器、绝热段、冷凝器、毛细芯与回流机制，也不能仅凭这一个 DEBORA 强制对流沸腾案例直接推导。本案例能够验证 CHT 接口与壁面相变计算路径，但不能单凭该算例保证热管模型可用。

最值得追踪的输出

wallBoiling:dDeparture、wallBoiling:fDeparture、wallBoiling:nucleationSiteDensity、wallBoiling:wetFraction、wallBoiling:mDot；同时记录流体进出口质量/焓通量，以及固体输入与流体接受热量之间的偏差。system/functions 已包含大量相关 functionObject 设置。

10 · 可复查证据

源码索引

所有主要公式均按 OpenFOAM Foundation OpenFOAM-14 当前仓库的 C++ 与 tutorial 文件核对。源文件链接指向 Foundation 官方仓库；历史沿革附有 Git commit 链接。没有将 openfoam.com（ESI/OpenCFD 分支）的实现直接套用到 Foundation-14。

[01] wallBoiling.C（RPI 热流、质量源项与迭代）

[02] wallBoiling.H（算法定义、模型引用）

[03] Lavieville.C（wetFraction）

[04] LemmertChawla.C（成核点密度）

[05] TolubinskiKostanchuk.C（脱离直径）

[06] KocamustafaogullariIshiiDepartureFrequency.C

[07] heatTransferLimitedPhaseChange.C（相间蒸发冷凝）

[08] multiphaseCoupledTemperatureFvPatchScalarField.C（流体侧多相耦合）

[09] coupledTemperatureFvPatchScalarField.C（固体侧耦合）

[10] alphatPhaseChangeWallFunctionFvPatchScalarField.C

[11] tutorial/constant/fluid/fvModels（本案例配置）

[12] tutorial/system/controlDict（区域求解器）

[13] tutorial/system/blockMeshDict（几何尺寸）

[14] tutorial/0/solid/T（外壁热流）

[15] tutorial/0/fluid/T.liquid（液体 CHT 温度）

[16] tutorial/Allrun（运行流程）

[17] tutorial/system/functions（输出与守恒检查）

[18] phaseChange.C（源项及潜热定义）

发布日期：2026-10-10。免责声明：本文做了静态源码核验与公式量纲核查，未在此环境中完成该算例的重新编译、运行及与实验曲线的数值复验。文中“当前结果”指源码公式的浏览器演算，而非 CFD 求解输出。

↑ 返回页首