Continuity Equation (Mass Conservation)
Continuity Equation (Mass Conservation): Theoretical Foundations
Overview
Professor, the continuity equation often appears alongside the Navier-Stokes equations. What does it represent?
It's the differential form of the law of conservation of mass. It's the mathematical expression of the physics that "fluid does not appear out of nowhere or disappear into nowhere." It's the foundation of all CFD.
General Form (Compressible Flow)
The continuity equation for any fluid is as follows.
Expanding it gives,
How should each term be interpreted?
$\partial\rho/\partial t$ is the rate of density change within a control volume. $\nabla\cdot(\rho\mathbf{u})$ is the net outflow rate of mass flux through the control volume surface. Their sum being zero means mass is conserved.
Material Derivative Form
Using the material derivative $D/Dt$, it can be written in another form.
This means "the rate of density change observed while riding along with a fluid particle is proportional to the divergence of the velocity field."
Case of Incompressible Flow
For incompressible flow, $\rho$ is constant, right?
Exactly. Since $D\rho/Dt = 0$, the continuity equation becomes very simple.
In component form,
This means "the velocity field is divergence-free." It is also called the solenoidal condition.
Are there cases for incompressible flow where density differs by location (like multiphase flow)?
Sharp observation. For example, in buoyancy-driven flow due to temperature differences (Boussinesq approximation), density variation is considered only in the buoyancy term, and the continuity equation remains $\nabla\cdot\mathbf{u}=0$. However, when handling water and air with the VOF (Volume of Fluid) method, density changes discontinuously, so the compressible form of the continuity equation must be used.
Integral Form
The integral form, which is the basis of the finite volume method, is also important. For a control volume $V$,
In FVM, this integral form is discretized for each cell. The core of the SIMPLE algorithm is adjusting the pressure field so that the mass flux balance through cell faces becomes zero.
When an Artery Clogs, Flow Velocity Becomes 16 Times Faster
The scariest everyday example directly linked to the continuity equation (mass conservation) is arterial stenosis. When the cross-sectional area is halved, flow velocity doubles, but if the diameter is halved, the area becomes one-fourth, so velocity quadruples. So what happens if atherosclerosis reduces the diameter by half? Velocity becomes 4 times, and since wall shear stress is proportional to the square of velocity, it becomes over 16 times. If such stress is applied to the wall of an artery with plaque buildup—well, it's going to detach, right? The continuity equation may seem like a modest formula in textbooks, but cardiac surgeons use it every day.
Computational Methods for Continuity Equation (Mass Conservation)
Role of Pressure in Incompressible Flow
For incompressible flow, pressure isn't in the continuity equation, right? How is pressure determined?
A question that gets to the core. In the incompressible NS equations, since density is known, the continuity equation $\nabla\cdot\mathbf{u}=0$ acts as a constraint to determine pressure, replacing the equation of state. Pressure can be interpreted as a "Lagrange multiplier to keep the velocity field divergence-free."
Pressure-Velocity Coupling Methods
Let's organize the representative algorithms.
| Method | Overview | Application |
|---|---|---|
| SIMPLE | Pressure Correction Method. Provisional velocity → pressure correction → velocity correction | Steady / Unsteady |
| SIMPLEC | Convergence-accelerated version of SIMPLE. No need for correction under-relaxation. | Steady problems |
| PISO | Two pressure corrections within one step. Suitable for unsteady problems. | Unsteady problems |
| Fractional Step Method | Velocity prediction → pressure Poisson → velocity projection | Unsteady / Academic codes |
| Coupled | Solves velocity and pressure simultaneously | When fast convergence is needed |
I feel like SIMPLE is the one I see most often.
The pressure correction equation in SIMPLE is a Poisson-type equation derived from the continuity equation.
Here $p'$ is the pressure correction value, $\mathbf{u}^*$ is the provisional velocity (velocity obtained tentatively from the momentum equation), $a_P$ is the diagonal coefficient of the momentum equation. It corrects the residual $\nabla\cdot\mathbf{u}^* \neq 0$ using $p'$.
Verification of Mass Conservation
How can I check if mass conservation is satisfied in the calculation results?
Check the balance of mass flow rates across each face.
In Fluent, you can check the flow rate at each boundary face via Report > Fluxes > Mass Flow Rate. An imbalance below 0.1% of the total flow rate is considered sufficiently converged.
Rhie-Chow Interpolation
On collocated grids (where velocity and pressure are placed at the same location), checkerboard instability occurs. Rhie-Chow interpolation is a technique that adds a pressure gradient correction to the cell face velocity to resolve this issue, and is standardly adopted in modern finite volume method solvers.
The continuity equation seems simple at first glance, but implementing it numerically requires considerable ingenuity, doesn't it?
Yes. Most of the difficulty in incompressible CFD boils down to the problem of "how to satisfy the continuity equation," i.e., the pressure-velocity coupling problem.
All Flow Meters Operate on the Continuity Equation
The mechanism of differential pressure flow meters (orifice flow meters) used in factory piping is the continuity equation itself. By reducing the cross-sectional area at the constriction to increase flow velocity, pressure drops according to Bernoulli's theorem. Measuring the pressure difference before and after allows calculating the flow rate inversely. Aircraft airspeed indicators (pitot tubes) work on the same principle. In numerical analysis too, it's fundamental to check after calculation whether "inlet flow rate = outlet flow rate" holds. If the continuity equation is not satisfied, there is either a leak somewhere in the mesh or a contradiction in the boundary conditions.