Boiling Model
Theory and Physics
Overview
Professor, how does CFD for boiling work? Can you simulate water boiling in a kettle?
Of course you can. However, boiling phenomena are extremely complex, involving multiphase flow problems where bubble generation/detachment/condensation at the wall and liquid film behavior are intertwined. In CFD, we use wall heat flux partitioning models to describe the heat transfer associated with boiling.
What is a wall heat flux partitioning model?
The most representative one is the RPI model (Rensselaer Polytechnic Institute model), proposed by Kurul & Podowski (1990). It's a concept that decomposes the total heat flux from the wall into three components.
Governing Equations
Please show me the equation for the RPI model.
The total wall heat flux $q_w$ is partitioned into the following three components.
The meaning of each component is as follows.
- $q_{fc}$: Single-phase convective heat transfer (portion where liquid covers the wall)
- $q_{quench}$: Quenching heat flux (transient heat transfer when cold liquid contacts the wall after bubble detachment)
- $q_{evap}$: Evaporative heat flux (heat directly consumed for bubble generation)
What are the specific equations for each?
The single-phase convection component is based on the area fraction of the wall in contact with the liquid $(1 - A_b)$.
The quenching component is related to the bubble detachment frequency $f$ and waiting time.
The evaporation component is determined from the active nucleation site density $N_a$ on the wall, bubble departure diameter $d_w$, and detachment frequency $f$.
$A_b$ is the fraction of the wall area covered by bubbles, right?
Yes. $A_b = \min\left(1,\; K \frac{\pi d_w^2}{4} N_a\right)$, where $K$ is an empirical constant. The bubble departure diameter $d_w$ often uses the Tolubinsky & Kostanchuk (1970) model.
How do you determine the active nucleation site density $N_a$?
The correlation by Lemmert & Chawla (1977) is commonly used.
Here, $\Delta T_{sup} = T_w - T_{sat}$ is the wall superheat. $n$ is typically 1.805, and $C$ is a parameter determined from experiments.
Boiling Regimes
Are there different types of boiling?
Boiling regimes transition according to wall superheat. Taking pool boiling as an example, it is organized by the Nukiyama curve (boiling curve).
| Regime | Wall Superheat | Characteristics |
|---|---|---|
| Natural Convection | $\Delta T_{sup} < 5$ K | No bubbles, single-phase convection |
| Nucleate Boiling | 5–30 K | Bubbles detach from wall, high heat transfer coefficient |
| Transition Boiling | 30–100 K | Unstable, alternating liquid and vapor film formation |
| Film Boiling | $> 100$ K | Vapor film covers wall, low heat transfer coefficient |
The transition point from nucleate to film boiling is CHF (Critical Heat Flux: Critical Heat Flux), which is the most important parameter in nuclear reactor safety assessment.
So it's dangerous to exceed CHF.
Exceeding CHF causes a rapid rise in wall temperature leading to burnout. In nuclear engineering, DNBR (Departure from Nucleate Boiling Ratio) is managed as a safety margin during design.
Nukiyama Curve—The Japanese Who Discovered Boiling's "Cliff"
In 1934, Shizuo Nukiyama of Tohoku University precisely measured the relationship between heat flux and wall temperature by varying the electrical power to a wire immersed in water. The result was an "S-shaped" curve (Nukiyama curve), where heat flux first increases to a maximum value (Critical Heat Flux, CHF), then wall temperature rises sharply before stabilizing again. This is because the moment CHF is exceeded, the wall becomes covered by a vapor film, drastically reducing the heat transfer coefficient. It represents a design limit that must never be exceeded in reactors or evaporators. This discovery became the starting point for boiling engineering and is still used worldwide 90 years later as a validation benchmark for CFD boiling models.
Physical Meaning of Each Term
- Temporal Term $\partial(\rho\phi)/\partial t$: Imagine the moment you turn on a faucet. At first, water comes out spluttering and unstable, but after a while, the flow becomes steady, right? This "during the change" is described by the temporal term. The pulsation of blood flow from a heartbeat, or the flow fluctuation each time an engine valve opens/closes—all are unsteady phenomena. So what is steady-state analysis? It looks only at "after sufficient time has passed and the flow has settled"—meaning setting this term to zero. Since computational cost drops significantly, solving first with a steady-state approach is a basic CFD strategy.
- Convection Term $\nabla \cdot (\rho \mathbf{u} \phi)$: If you drop a leaf into a river, what happens? It gets carried downstream by the flow, right? This is "convection"—the effect where fluid motion transports things. Warm air from a heater reaching the far side of a room is also because the "carrier," air, transports heat by convection. Here's the interesting part—this term contains "velocity × velocity," making it nonlinear. That is, as flow speed increases, this term rapidly strengthens, making control difficult. This is the root cause of turbulence. A common misconception: "Convection and conduction are similar" → They're completely different! Convection is carried by flow, conduction is transmitted by molecules. There's an order of magnitude difference in efficiency.
- Diffusion Term $\nabla \cdot (\Gamma \nabla \phi)$: Have you ever put milk in coffee and left it? Even without stirring, after a while it naturally mixes. That's molecular diffusion. Now a question—honey and water, which flows more easily? Obviously water, right? Honey has high viscosity ($\mu$), so it flows poorly. Higher viscosity strengthens the diffusion term, making the fluid move "sluggishly." In low Reynolds number flow (slow, viscous), diffusion dominates. Conversely, in high Re number flow, convection overwhelms and diffusion plays a supporting role.
- Pressure Term $-\nabla p$: When you push a syringe plunger, liquid shoots out forcefully from the needle, right? Why? The plunger side is high pressure, the needle tip is low pressure—this pressure difference provides the force pushing the fluid. Dam discharge works on the same principle. On a weather map, where isobars are tightly packed? That's right, strong winds blow. "Flow arises where there is a pressure difference"—this is the physical meaning of the pressure term in the Navier-Stokes equations. A point of confusion here: "Pressure" in CFD is often gauge pressure, not absolute pressure. If results go wrong immediately after switching to compressible analysis, mixing up absolute/gauge pressure might be the cause.
- Source Term $S_\phi$: Heated air rises—why? Because it becomes lighter (lower density) than its surroundings, so buoyancy pushes it up. This buoyancy is added to the equation as a source term. Other examples: chemical reaction heat from a gas stove flame, Lorentz force acting on molten metal in a factory's electromagnetic pump... These are all actions that "inject energy or force into the fluid from outside," expressed by the source term. What happens if you forget the source term? In a natural convection analysis, forgetting to include buoyancy means the fluid doesn't move at all—a physically impossible result like turning on a heater in a winter room but the warm air doesn't rise.
Assumptions and Applicability Limits
- Continuum Assumption: Valid for Knudsen number Kn < 0.01 (molecular mean free path ≪ characteristic length)
- Newtonian Fluid Assumption: Shear stress and strain rate have a linear relationship (non-Newtonian fluids require viscosity models)
- Incompressibility Assumption (for Ma < 0.3): Treat density as constant. For Mach number ≥ 0.3, consider compressibility effects
- Boussinesq Approximation (Natural Convection): Consider density change only in the buoyancy term, using constant density in other terms
- Non-applicable Cases: Rarefied gas (Kn > 0.1), supersonic/hypersonic flow (shock capturing required), free surface flow (VOF/Level Set, etc. required)
Dimensional Analysis and Unit Systems
| Variable | SI Unit | Notes / Conversion Memo |
|---|---|---|
| Velocity $u$ | m/s | When converting from volumetric flow rate for inlet conditions, pay attention to cross-sectional area units |
| Pressure $p$ | Pa | Distinguish between gauge and absolute pressure. Use absolute pressure for compressible analysis |
| Density $\rho$ | kg/m³ | Air: ~1.225 kg/m³ @20°C, Water: ~998 kg/m³ @20°C |
| Viscosity Coefficient $\mu$ | Pa·s | Be careful not to confuse with kinematic viscosity $\nu = \mu/\rho$ [m²/s] |
| Reynolds Number $Re$ | Dimensionless | $Re = \rho u L / \mu$. Indicator for laminar/turbulent transition |
| CFL Number | Dimensionless | $CFL = u \Delta t / \Delta x$. Directly related to time step stability |
Numerical Methods and Implementation
Details of Numerical Methods
When implementing a boiling model in CFD, what framework is used to solve it?
Boiling analysis is mainly based on the Euler-Euler two-fluid model, with the RPI model incorporated as a wall boundary condition. The continuity and momentum equations are solved for the liquid and vapor phases separately.
The transport equation for the vapor volume fraction $\alpha_v$ includes source terms $\dot{m}$ for evaporation and condensation.
The evaporation rate is determined from the RPI model's $q_{evap}$, right?
Exactly. The evaporative mass flux at the wall is as follows.
Condensation in the bulk is calculated from the interfacial heat transfer coefficient obtained via the Ranz-Marshall correlation, based on heat exchange with the subcooled liquid around the bubbles.
Bubble Force Models
How do bubbles generated by boiling move?
Modeling the forces acting on bubbles is important. The following interphase forces are considered.
| Force | Model | Role |
|---|---|---|
| Drag Force | Schiller-Naumann, Ishii-Zuber | Governs bubble velocity difference |
| Lift Force | Tomiyama | Lateral force due to velocity gradient |
| Wall Lubrication Force | Antal | Pulls bubbles away from the wall |
| Turbulent Dispersion Force | Lopez de Bertodano | Turbulent diffusion of bubbles |
| Virtual Mass Force | Auton | Acceleration effect |
Can Tomiyama's lift force change sign?
Good question. When the bubble diameter exceeds a critical value based on the Eötvös number $Eo$, the direction of the lift force reverses. Small bubbles move towards the wall, large bubbles move towards the pipe center. This is important physics determining wall-peaking and core-peaking void distributions.
Treatment of Wall Functions
Can normal wall functions be used on boiling surfaces?
They cannot. This is because bubble agitation significantly alters the near-wall flow structure compared to single-phase flow. Software like Fluent implements boiling-specific wall functions that perform wall temperature calculations consistent with the RPI model.
For wall meshing, it's more important that the first cell height is appropriate relative to the bubble departure diameter $d_w$ than constraints on $y^+$. Generally, it's recommended that the first cell height be at least $d_w$.
Time Step and Stability
Why does boiling analysis tend to diverge easily?
Because the volume expansion accompanying evaporation is rapid, generating large local volume source terms. The following countermeasures are effective.
- Gradually increase wall superheat (ramp-up)
- First obtain a single-phase steady solution, then activate boiling from there
- Set time step to $10^{-4}$ s or less
- Apply under-relaxation (0.3–0.5) to the volume fraction equation
Dissecting the RPI Wall Boiling Model—The De Facto Standard for Nucleate Boiling CFD
The de facto standard for analyzing nucleate boiling in CFD is the RPI model. It decomposes the heat flux from the wall to the liquid into q_evap (evaporation) + q_quench (quenching) + q_conv (single-phase convection), calculating each from nucleation site density, bubble departure diameter, and detachment frequency. This model is implemented in ANSYS Fluent/CFX, but since the coefficient in the nucleation site density equation strongly depends on experimental conditions, caution is needed when extrapolating to new conditions. One researcher warns that "the RPI model is a stack of three correlations, so the uncertainty of each correlation is amplified multiplicatively," and errors in CHF prediction exceeding 30% are not uncommon.
Upwind Difference (U
Related Topics
なった
詳しく
報告