What is the Boussinesq approximation?

It is an approximation that considers density variations only in the buoyancy term (the body force term), while treating density as constant in the mass conservation and inertial terms. It is applicable to natural convection with small temperature differences (β·ΔT ≪ 1) and allows buoyancy-driven flows to be handled with an incompressible solver, significantly reducing computational cost.

What happens when the Rayleigh number exceeds 10^9?

Around Ra > 10^9, buoyancy-driven flow transitions from laminar to turbulent. After this turbulent transition, turbulence analysis using RANS models (such as k-ε or k-ω SST) or LES becomes necessary, and mesh design that sufficiently resolves the thermal boundary layer near the wall becomes critical.

What is the most important point in mesh design for CFD of buoyancy-driven flow?

It is to sufficiently resolve the thermal and velocity boundary layers near the heated or cooled walls. Since the boundary layer thickness can be estimated as δ_T ∝ L·Ra^(-1/4), a minimum of 10 cell layers in the wall-normal direction is required, with the first cell height designed for y+ < 1.

How can you address convergence issues in natural convection CFD?

Key countermeasures include: (1) setting initial conditions to a zero-velocity field with a linear temperature distribution, (2) using a parameter continuation method to gradually increase the Rayleigh number, (3) reducing the relaxation factors for SIMPLE-type schemes (e.g., pressure ~0.2, momentum ~0.5), and (4) starting with a small time step in a transient simulation.

CFD Simulation of Buoyancy-Driven Flow

Q: How can you address convergence issues in natural convection CFD?

Key countermeasures include: (1) setting initial conditions to a zero-velocity field with a linear temperature distribution, (2) using a parameter continuation method to gradually increase the Rayleigh number, (3) reducing the relaxation factors for SIMPLE-type schemes (e.g., pressure ~0.2, momentum ~0.5), and (4) starting with a small time step in a transient simulation.

Category: 熱解析 > 自然対流 | Integrated 2026-04-12

Buoyancy-driven natural convection flow simulation showing temperature contours and velocity vectors in an enclosed cavity — 浮力駆動流れの温度コンター図と速度ベクトル場 — 加熱壁面近傍で上昇流、冷却壁面近傍で下降流が発生

Theory and Physics

What is Buoyancy-Driven Flow?

🧑‍🎓

Professor, is buoyancy-driven flow the phenomenon where hot bathwater rises to the top?

🎓

Yes, exactly that. Temperature differences cause changes in fluid density, and in a gravitational field, lighter (higher temperature) fluid rises while heavier (lower temperature) fluid descends. This is the essence of buoyancy-driven flow, or natural convection.

🧑‍🎓

Does it happen on much larger scales too, not just in bathwater?

🎓

The scale is truly vast. Important engineering applications include:

Natural ventilation in buildings — Temperature differences indoors create airflow, exhausting heat through windows and vents. The basis for ventilation tower design in passive houses.
Natural air cooling of electronic enclosures — Dissipating heat from IC chips without fans. Essential for IoT devices and LED lighting design.
Emergency cooling of nuclear reactors — Passive safety systems (ECCS) that maintain core cooling via natural circulation alone during power loss.
Data center thermal management — Interaction between thermal stratification between server racks and HVAC airflow.
Geophysics — Mantle convection, atmospheric convection cells (Hadley circulation), oceanic thermohaline circulation.

🧑‍🎓

Wow, even nuclear reactors! So even if pumps stop, cooling can be maintained just by buoyancy force?

🎓

Exactly. Since the Fukushima Daiichi accident, the design verification of passive safety systems has become more stringent. A natural circulation loop forms where coolant heated in the core rises and cools down in the steam generator before descending. Accurately predicting this flow with CFD is crucial for safety assessment.

Boussinesq Approximation

🧑‍🎓

When solving buoyancy-driven flow with CFD, do you calculate all the density changes due to temperature accurately? That seems tough...

🎓

Good question. That's where the Boussinesq Approximation comes in. It's the most standard approach for natural convection CFD, proposed by Joseph Boussinesq in 1903.

The concept is simple: density variation is considered only in the buoyancy term (body force term), while in all other terms (the continuity equation of mass conservation, the inertial term of the momentum equation), density is treated as constant $\rho_0$.

$$ \rho \approx \rho_0 \quad \text{(except for buoyancy term)} $$

$$ \rho = \rho_0 \left[1 - \beta (T - T_0)\right] \quad \text{(buoyancy term only)} $$

Here, $\beta$ is the volumetric expansion coefficient [1/K], $T_0$ is the reference temperature, and $\rho_0$ is the density at the reference temperature.

🧑‍🎓

I see, so it's only applicable when density changes are small?

🎓

Exactly. The applicability condition is $\beta \cdot \Delta T \ll 1$ (roughly $\beta \Delta T < 0.1$). For example, for air ($\beta \approx 1/300$ K$^{-1}$ at 300K), sufficient accuracy is achieved for temperature differences up to about 30K. However, in combustion fields ($\Delta T \sim 1000$K or more) or with cryogenic fluids, the Boussinesq approximation breaks down, requiring a fully compressible approach or a "non-Boussinesq" model that directly calculates density as a function of temperature.

Limits of Boussinesq Approximation and Alternative Methods

Volumetric expansion coefficient $\beta$: For ideal gases, $\beta = 1/T$ (inverse of absolute temperature). For liquids, it has temperature dependence; water shows unique behavior with $\beta=0$ at 4$^\circ$C (maximum density).
Non-Boussinesq model: Density is calculated directly from the equation of state $\rho = \rho(T, p)$. Essential for large temperature difference problems ($\beta \Delta T > 0.1$). In Ansys Fluent, the "Incompressible Ideal Gas" model is used.
Fully compressible solution: Also considers pressure fluctuations on the speed of sound scale. Used in fire simulations (FDS, etc.). Computational cost is highest but applicability is broad.

Dimensionless Numbers — Gr, Ra, Pr, Nu

🧑‍🎓

Are there many dimensionless numbers that characterize natural convection? Like the Reynolds number?

🎓

In forced convection, the Reynolds number $Re$ is the governing parameter, but in natural convection, since no external velocity is given, $Re$ cannot be defined. Instead, dimensionless numbers representing the strength of buoyancy take the lead.

Grashof Number — Ratio of buoyancy to viscous forces:

$$ \mathrm{Gr} = \frac{g \beta \Delta T \, L^3}{\nu^2} $$

Prandtl Number — Ratio of momentum diffusivity to thermal diffusivity:

$$ \mathrm{Pr} = \frac{\nu}{\alpha} = \frac{\mu c_p}{k} $$

Rayleigh Number — Governing parameter for natural convection:

$$ \mathrm{Ra} = \mathrm{Gr} \cdot \mathrm{Pr} = \frac{g \beta \Delta T \, L^3}{\nu \alpha} $$

Nusselt Number — Strength of convective heat transfer (ratio to conduction only):

$$ \mathrm{Nu} = \frac{h L}{k} $$

🧑‍🎓

Ra = Gr $\times$ Pr, so does a larger Rayleigh number mean more vigorous convection?

🎓

That's right. The Rayleigh number represents "whether buoyancy overcomes thermal diffusion (heat conduction)". For small Ra, heat moves only by conduction (Nu $\approx$ 1). As Ra increases, convection becomes dominant and Nu grows larger. Here's a practical numerical sense:

$\mathrm{Ra} < 10^3$: No convection (pure conduction), Nu $\approx$ 1
$10^3 < \mathrm{Ra} < 10^9$: Laminar natural convection
$\mathrm{Ra} > 10^9$: Turbulent transition

Typical Prandtl Numbers for Fluids

Fluid	Pr (room temp., atm. pressure)	Characteristics
Air	0.71	Temperature and velocity boundary layers are nearly the same thickness
Water	7.0 (20$^\circ$C)	Temperature boundary layer $<$ velocity boundary layer
Engine Oil	$\sim$1000	Temperature boundary layer is very thin
Liquid Metal (Na)	0.004〜0.01	Temperature boundary layer $\gg$ velocity boundary layer

Full Form of Governing Equations

🧑‍🎓

So, can you show me all the governing equations using the Boussinesq approximation?

🎓

It's the set of incompressible Navier-Stokes equations and the energy equation under the Boussinesq approximation. Solving these three equations simultaneously is the basis of CFD.

Continuity Equation (mass conservation):

$$ \nabla \cdot \mathbf{u} = 0 $$

Momentum Conservation Equation (Navier-Stokes equation + buoyancy term):

$$ \rho_0 \frac{\partial \mathbf{u}}{\partial t} + \rho_0 (\mathbf{u} \cdot \nabla)\mathbf{u} = -\nabla p' + \mu \nabla^2 \mathbf{u} + \rho_0 \beta (T - T_0) \mathbf{g} $$

Here, $p' = p - \rho_0 \mathbf{g} \cdot \mathbf{x}$ is the modified pressure (with hydrostatic component removed). The final term on the right side is the buoyancy term, the core part of the Boussinesq approximation.

Energy Equation:

$$ \rho_0 c_p \frac{\partial T}{\partial t} + \rho_0 c_p (\mathbf{u} \cdot \nabla) T = k \nabla^2 T + Q $$

🧑‍🎓

The last part of the momentum equation, $\rho_0 \beta (T - T_0) \mathbf{g}$, is the buoyancy term, right? Where temperature is high ($T > T_0$), an upward force is generated.

🎓

You understand perfectly. Without this buoyancy term, the fluid wouldn't move even with a temperature difference. In other words, the buoyancy term is the driving force of natural convection. And through this term, the momentum and energy equations are bidirectionally coupled. The temperature field drives the velocity field, and the velocity field transports the temperature field—this nonlinear coupling is what makes natural convection CFD both interesting and challenging.

Nu Number Correlations

🧑‍🎓

When validating CFD results, are there theoretical values for the Nu number or something similar?

🎓

There are several empirical correlations based on experimental data. Let me introduce the formulas to check first for CFD validation.

Laminar natural convection on a vertical flat plate (Churchill-Chu equation, full Ra range):

$$ \overline{\mathrm{Nu}}_L = \left\{ 0.825 + \frac{0.387 \, \mathrm{Ra}_L^{1/6}}{\left[1 + (0.492/\mathrm{Pr})^{9/16}\right]^{8/27}} \right\}^2 $$

Horizontal heated plate (top surface):

$$ \overline{\mathrm{Nu}} = 0.54 \, \mathrm{Ra}^{1/4} \quad (10^4 \leq \mathrm{Ra} \leq 10^7) $$

$$ \overline{\mathrm{Nu}} = 0.15 \, \mathrm{Ra}^{1/3} \quad (10^7 \leq \mathrm{Ra} \leq 10^{11}) $$

Enclosed cavity (vertical wall heating/cooling, de Vahl Davis benchmark solution):

Ra	Nu (Benchmark)	Application
$10^3$	1.118	Code verification (V&V)
$10^4$	2.243	Code verification
$10^5$	4.519	Code verification
$10^6$	8.800	Code verification

🧑‍🎓

The de Vahl Davis benchmark, I see it often in CFD textbooks! So you check if your code matches these values.

🎓

Exactly. Published in 1983, this benchmark solution is still used as the "litmus test" for natural convection CFD codes. Along with checking mesh convergence, it's standard practice to first confirm the solver's validity with this square cavity problem.

Turbulent Transition and Ra Number Threshold

🧑‍🎓

Earlier you said Ra $> 10^9$ is turbulent, is that like Re > 2300 for forced convection?

🎓

The concept is the same. However, the transition Rayleigh number depends on geometry:

Vertical flat plate: $\mathrm{Ra}_\mathrm{crit} \approx 10^9$ (Gr $\approx 10^9$)
Horizontal heated plate (top surface): $\mathrm{Ra}_\mathrm{crit} \approx 10^7$
Enclosed cavity: $\mathrm{Ra}_\mathrm{crit} \approx 10^8$ (depends on aspect ratio)

After turbulent transition, RANS, LES, or DNS becomes necessary. Especially in data center or large-space HVAC CFD, Ra often exceeds $10^{10}$, so the choice of turbulence model greatly influences result accuracy.

Coffee Break Trivia Corner

Rayleigh-Benard Convection and Beautiful Patterns

In a horizontal fluid layer heated from below and cooled from above, when the Rayleigh number exceeds a critical value ($\mathrm{Ra}_c = 1708$), regular convection cells (Benard cells) suddenly appear. Henri Benard discovered this experimentally in 1900, and Lord Rayleigh explained it theoretically in 1916. Deriving this critical Ra number is a classic problem in fluid dynamics stability theory and is still widely used as a CFD verification problem. Nature chooses the most beautiful hexagonal pattern at the minimum Ra number.

Numerical Methods and Implementation

Finite Volume Method Discretization

🧑‍🎓

For CFD of buoyancy-driven flow, which is more mainstream, the finite element method or the finite volume method?

🎓

In CFD, the Finite Volume Method (FVM) is overwhelmingly dominant. Ansys Fluent, STAR-CCM+, OpenFOAM—all are FVM-based solvers. The reason is simple: FVM strictly satisfies conservation laws (mass, momentum, energy) at the discrete level.

In FVM, the governing equations are integrated over a cell volume $V$ and transformed into surface integrals using Gauss's divergence theorem:

$$ \frac{\partial}{\partial t} \int_V \rho \phi \, dV + \oint_S \rho \phi \mathbf{u} \cdot d\mathbf{S} = \oint_S \Gamma \nabla \phi \cdot d\mathbf{S} + \int_V S_\phi \, dV $$

Here, $\phi$ is a generic scalar (velocity component, temperature, etc.), $\Gamma$ is the diffusion coefficient, and $S_\phi$ is the source term. In buoyancy-driven flow, the buoyancy term $\rho_0 \beta (T - T_0) g_i$ enters the $S_\phi$ of the momentum equation.

🧑‍🎓

Are there any points to note about discretizing the convection term? Any problems specific to natural convection?

🎓

Good point. Natural convection velocities are lower compared to forced convection, so the cell Peclet number is often on the order of O(1). In this regime:

First-Order Upwind (UDS) has excessive numerical diffusion, blurring the temperature boundary layer.
Central Differencing (CDS) is accurate but prone to oscillations.
Second-Order Upwind (SOU) or QUICK schemes offer a good balance.

In practice, if you don't use a scheme of second-order accuracy or higher, the Nu number can be underestimated by 20-30%. This is one of the most common mistakes in natural convection CFD.

Pressure-Velocity Coupling

🧑‍🎓

The incompressible Navier-Stokes equations don't have a pressure equation, so they require a special solution method, right?

🎓

That's right. In incompressible flow, density is constant so there's no equation of state to directly link pressure and density. Therefore, pressure is determined using the continuity equation ($\nabla \cdot \mathbf{u} = 0$) as a constraint. Representative methods:

SIMPLE (Semi-Implicit Method for Pressure-Linked Equations): Standard for segregated solvers. Suitable for steady-state calculations.
SIMPLEC: Modified version of SIMPLE. No pressure correction relaxation is needed, sometimes converging faster.
PISO: Suitable for transient calculations. Improves accuracy with multiple correction steps per time step.
Coupled Solver: Solves pressure and velocity simultaneously. Uses more memory but dramatically improves convergence for natural convection.

🧑‍🎓

Which one is recommended for natural convection?

🎓

In practice, I strongly recommend the Coupled Solver. The reason is that in natural convection, the coupling between the pressure and velocity fields is very strong, and segregated solvers (SIMPLE family) can take thousands of iterations to converge. The Coupled Solver often finishes in a few hundred iterations. However, memory usage increases by 1.5-2 times, so for large-scale problems, this trade-off needs to be considered.

Turbulence Model Selection

🧑‍🎓

For turbulent natural convection with Ra $> 10^9$, which turbulence model is appropriate?

🎓

This is a very important point. Turbulence in natural convection has a different character than in forced convection, so choosing the wrong model can lead to large errors.

Model