Conservation of Mass: The Continuity Equation

Post 10 derived the material derivative \(D\mathbf{u}/Dt\), the \(a\) in \(F = ma\). Newton's law is not the only rule a flow obeys: mass also cannot appear or vanish. This post turns that rule into the continuity equation and shows that for water it becomes the condition \(\nabla \cdot \mathbf{u} = 0\) announced at the end of Post 10.

The path: count the mass crossing the faces of a box fixed in space, first in one direction for a gas and for water, then in two and three directions, where the count produces the divergence. The post ends with what kind of equation \(\nabla \cdot \mathbf{u} = 0\) is, and when a real fluid may be treated as obeying it.

Mass is never created or destroyed. To turn that into an equation, bolt an imaginary box in place, like the Eulerian speed sensor on the bridge pillar in Post 10. The box has no walls of its own; fluid passes freely through its faces. Since no mass can appear or vanish inside, the mass in the box can only change by crossing a face.

For example, 5 kg/s of fluid enter through the left face, 3 kg/s leave through the right face, and nothing crosses the other faces. The box gains 2 kg every second. In words, the rate of change inside is "in minus out":

Where \(M\) denotes the mass inside the box (kg), \(dM/dt\) denotes how fast it changes (kg/s), and "in" and "out" denote the total mass entering and leaving through the faces per second (kg/s). In the example, \(dM/dt = 5 - 3 = 2\) kg/s. The derivative is an ordinary \(d\) because \(M\) is a single number that depends only on time; the box does not move.

What enters either leaves again or piles up. Water is nearly impossible to squeeze, so whatever goes in comes straight back out: in = out. A gas can be squeezed, so the extra 2 kg/s piles up and the density inside rises, until the higher pressure pushes the gas out again (that push is a force and belongs to the next post). The rest of this post makes this rule precise, starting with its right-hand side: where does a number like "5 kg/s" come from?

Take a face of 0.5 m², with fluid moving straight through it at 2 m/s. In one second, every bit of fluid that was within 2 m upstream of the face makes it through; fluid further back does not reach the face in time. That fluid fills a slab 2 m long with the face as its cross-section: 2 × 0.5 = 1 m³. For water, at about 1000 kg/m³, that is 1000 kg every second.

In general, speed times area gives volume per second, m/s · m² = m³/s, and multiplying by the density in kg/m³ turns it into mass per second:

Where \(\rho\) denotes the density (kg/m³), \(A\) denotes the area of the face (m²), and \(u_\perp\) denotes the part of the velocity at right angles to the face (m/s). Units: kg/m³ · m/s · m² = kg/s.

The three factors multiply because each one scales the result proportionally: twice the speed gives a slab twice as long, twice the area a slab twice as wide, twice the density twice the kilograms in the same slab. And if any one of them is zero, nothing gets through.

Only \(u_\perp\) counts. Fluid sliding along a face never crosses it, like walking alongside a door instead of through it. If the fluid above moved 2 m/s through the face of 0.5 m² and also 1.5 m/s along it, the same 1 m³ would cross per second; the sideways part only shifts the slab sideways.

The product of the first two factors gets its own name. The mass flux \(\rho u_\perp\) is the mass crossing a face per second per square metre, with units kg/m³ · m/s = kg/(m²·s). Multiplying it by the face area gives back the mass per second. Applied to two opposite faces of a thin box, it gives the first real equation.

A gas flows along a straight pipe of constant cross-section \(A\). Its density \(\rho(x,t)\) and velocity \(u(x,t)\) change along the pipe. The box is the slice of pipe between \(x\) and \(x + \Delta x\).

A worked example at one frozen instant: faces at \(x = 1\) m and \(x = 1.1\) m, \(A = 0.5\) m². 3 kg/s enter through the left face and 2.8 kg/s leave through the right face, so the box gains 0.2 kg every second. That number alone says little about the gas, because a longer box would gain a different amount. Density is mass per volume, so dividing the gain by the box volume, \(A \cdot \Delta x = 0.5 \times 0.1 = 0.05\) m³, gives the rate at which the density rises: 0.2 / 0.05 = 4 kg/m³ every second.

The general version does the same with symbols. Mass per second is density × speed × area, so "in" is \(\rho u A\) at the left face and "out" is \(\rho u A\) at the right face. The gain divided by the box volume is

Where \((\cdot)\big|_{x}\) denotes a quantity evaluated at the face at position \(x\), \(\Delta x\) denotes the length of the box (m), and \(\partial \rho/\partial t\) denotes the rate at which the density in the box changes ((kg/m³)/s). The derivative is partial because the box stays at a fixed \(x\): it is the rate of change at a fixed spot, the \(\partial/\partial t\) of Post 10. The \(\approx\) is there because the density varies slightly inside the box; this error vanishes as the box gets thinner.

\(A\) is the same at both faces, so it cancels. Swapping the order of the difference and putting a minus sign in front gives the form "(right − left) / \(\Delta x\)":

This form is useful because, for a thin box, the difference of a quantity across the box is its slope times the width: \((\rho u)\big|_{x+\Delta x} - (\rho u)\big|_{x} \approx \frac{\partial (\rho u)}{\partial x}\,\Delta x\). This becomes exact as \(\Delta x \to 0\), and it is the same move Post 3 made with the tension at the two ends of a string segment. The \(\Delta x\) cancels, and what is left is the slope:

Where \(\partial (\rho u)/\partial x\) denotes how much the mass flux \(\rho u\) changes per metre along the pipe. The example fits: the mass flux is 3 / 0.5 = 6 kg/(m²·s) at the left face and 2.8 / 0.5 = 5.6 kg/(m²·s) at the right face, a slope of (5.6 − 6) / 0.1 = −4 kg/(m³·s), and the minus sign turns it into the +4 kg/m³ per second found above.

The minus sign is the "in minus out" of the bookkeeping principle. If the mass flux grows from left to right, more leaves through the right face than enters through the left, and the density falls. Units: the left side is (kg/m³)/s, the right side is kg/(m²·s) per metre, which is kg/(m³·s). They match.

Now let the fluid be water. Its density is treated as exactly constant, so the density in the box never changes: \(\partial \rho/\partial t = 0\), and the gas equation becomes \(0 = -\partial (\rho u)/\partial x\). Since \(\rho\) is the same at every \(x\), it is a constant factor and can be pulled out of the slope: \(\partial (\rho u)/\partial x = \rho\, \partial u/\partial x\). This step is allowed only because \(\rho\) is constant; a density that changed along \(x\) would contribute a slope of its own. Dividing by \(\rho\), which is positive and so never zero, leaves

In a pipe of constant width, water has the same speed everywhere along the pipe at any given instant. The whole column can still speed up or slow down together over time, for example when a tap opens, but at any point in time, the speed is the same at every \(x\). Take your time to look at the widget at the top of this blog post and change the fluid type to water.

We just found that water in a pipe has the same speed everywhere along it. But in Post 10, the river sped up where it got narrower. The reason: in the gas box, the area \(A\) was the same at both faces, so it cancelled. In a river, the width changes along the flow, so the area depends on position, \(A(x)\), and no longer cancels. The density \(\rho\) is still constant, so it can still be pulled out in front:

Where \(A(x)\) denotes the cross-sectional area at position \(x\) (m²), and \(V\) denotes the volume of the box (m³), which is no longer exactly \(A \cdot \Delta x\). For water, the left side is 0. Multiplying both sides by \(V\) and dividing by \(\rho\) leaves only the difference:

So speed × area, the volume per second in m³/s, is the same at both faces, and therefore at every cross-section of the river. Where the area shrinks, the speed has to grow. A river 5 m wide flowing at 1 m/s carries 5 m³/s per metre of depth; where it narrows to 2.5 m, the same 5 m³/s must pass, so it flows at 2 m/s.

Narrowing from 5 m to 2.5 m means each bank comes in by 1.25 m. Water flowing next to a bank cannot keep going straight, or it would run into the bank: it must also move sideways, toward the middle. Its velocity has a component \(u_y\) across the river.

The 1D model never tracked \(u_y\). \(A(x)\) was a trick in which the walls did the sideways bookkeeping for us: by counting everything that passes a whole cross-section, the model never had to ask how the water gets from the wide part into the narrow middle. To properly account for this the box needs faces in both directions.

The bookkeeping is the same as for the gas box, with two small changes. Water is incompressible, so we count volume instead of mass; the constant \(\rho\) would cancel anyway, as it did for the river. And the flow is now 2D, \((u_x, u_y)\), with nothing varying along \(z\) (3D follows below). The box is \(\Delta x \times \Delta y\) × 1 m deep. Its left and right faces are crossed only by \(u_x\), since \(u_y\) slides along them, and their area, which was \(A\) in the pipe, is now \(\Delta y\) × 1 m.

Take \(\Delta x = \Delta y = 0.1\) m, with \(u_x = 2.0\) m/s at the left face and 2.08 m/s at the right face. The slope is (2.08 − 2.0) / 0.1 = 0.8 (m/s)/m, so by slope × width the water leaves 0.8 × 0.1 = 0.08 m/s faster than it enters. Multiplying this speed difference by the face area gives the volume this box loses per second: 0.08 × 0.1 = 0.008 m³/s. Dividing by the box volume, as for the gas box, turns this property of the box into a property of the fluid: 0.008 / 0.01 = 0.8 /s.

In symbols:

Where \(\partial u_x/\partial x\) denotes how much \(u_x\) changes per metre along \(x\), and \(\Delta y\) denotes the height of the box (m). The box sizes cancel, and only the slope is left. A smaller box gives the same answer: with 0.05 m sides, it loses 0.8 × 0.05 × 0.05 = 0.002 m³/s and has a volume of 0.0025 m³, so 0.002 / 0.0025 = 0.8 /s again.

The unit 1/s is m³/s divided by m³: the fraction of the box's own volume lost per second. 0.8 /s means 80 % of the box's volume per second.

The top and bottom faces are the left and right faces turned by 90°: only \(u_y\) crosses them, and each has area \(\Delta x\) × 1 m. The same three steps (extra speed \(\frac{\partial u_y}{\partial y}\,\Delta y\), times face area \(\Delta x\) × 1 m, divided by the box volume) give \(\partial u_y/\partial y\). Adding both pairs gives the total fraction of the box's volume lost per second:

Where \(\partial u_y/\partial y\) denotes how much \(u_y\) changes per metre along \(y\) (1/s). For an incompressible fluid this must be 0: the box can neither lose water nor gain it. If the \(x\) pair drains 0.8 /s, the \(y\) pair must refill it: \(\partial u_y/\partial y = -0.8\) /s. In the narrowing river, this is the inward squeeze from the banks. With \(y\) measured across the river from its centre line, water above the centre line moves down toward it and water below moves up, so \(u_y\) decreases as \(y\) increases. A widening river does the opposite. Each row of the table sums to 0:

In 3D, a third pair of faces, at right angles to \(z\), adds \(\partial u_z/\partial z\) in exactly the same way. For example, 0.8 − 0.5 − 0.3 = 0: the flow stretches along \(x\), is squeezed along \(y\) and \(z\), and the box keeps its water. This sum is the divergence:

Where \(\partial u_z/\partial z\) denotes how much \(u_z\) changes per metre along \(z\), and \(\nabla \cdot \mathbf{u}\) denotes the divergence of \(\mathbf{u}\): the net fraction of a small box's volume flowing out per second (1/s). The notation is a dot product, as in Post 10. \(\nabla = (\partial/\partial x,\; \partial/\partial y,\; \partial/\partial z)\) is a vector of slope-takers, and "multiplying" a slope-taker by a component means applying it: \(\partial/\partial x\) times \(u_x\) is \(\partial u_x/\partial x\). Matching entries are combined and added, as in any dot product. For an incompressible fluid:

One plain number per place, and it must be zero everywhere.

We've had \(\nabla\) appear in multiple combinations now, and the position of the dot tells them apart:

Divergence does not measure how much motion there is, only whether a box sends out more than it receives. Drag and resize the box below in five 2D flows (coefficients in 1/s, positions in m). Each face arrow shows the water crossing per second: blue in, green out.

Spring and drain. In the spring, \(\mathbf{u} = (0.5x,\; 0.5y)\), the outer face of every box is faster than the inner one, so \(\nabla \cdot \mathbf{u} = 0.5 + 0.5 = +1\) /s. The drain, \(\mathbf{u} = (-0.5x,\; -0.5y)\), is the reverse: \(-1\) /s. Moving or shrinking the box changes the face amounts, never the divergence: the amounts belong to the box, the divergence to the flow.

Steady wind. \(\mathbf{u} = (2,\; 0)\) is the fastest flow, yet \(\nabla \cdot \mathbf{u} = 0\): the same 2 m/s crosses the left and right faces.

Whirlpool. \(\mathbf{u} = (-y,\; x)\) circles counterclockwise, with \(\nabla \cdot \mathbf{u} = \partial(-y)/\partial x + \partial x/\partial y = 0\). At (0, 1), \(u_x = -y\) depends only on height, and the left and right faces span the same heights, so 0.50 m³/s enters on the right and 0.50 m³/s leaves on the left. At (1, 0), the same holds for the top and bottom faces, since \(u_y = x\) depends only on \(x\). The water passes through each box like a lane of traffic.

Squeeze. \(\mathbf{u} = (0.8x,\; -0.8y)\) is the narrowing river from the table above. Here the faces of a pair do not balance: the \(x\) faces lose 0.20 m³/s, the \(y\) faces bring in 0.20 m³/s. \(\nabla \cdot \mathbf{u} = 0\) only asks for the sum.

A solver never sees a formula like \(\mathbf{u} = (0.5x,\; 0.5y)\). It only has values stored at grid points, so the slopes in the divergence have to be computed from neighbouring values. The slope at a grid point is approximated by the difference between its two neighbours divided by their distance, a central difference. In 2D:

Where \(i\) denotes the grid index along \(x\) and \(m\) the grid index along \(y\) (as in Post 8), \(u_x[i,m]\) and \(u_y[i,m]\) denote the stored velocity components at grid point \((x_i, y_m)\), \(\Delta x\) and \(\Delta y\) now denote the grid spacings, and \(2\Delta x\) and \(2\Delta y\) are the distances between the two neighbours used.

Check on the spring at (1, 1) with spacing 0.1 m: \(u_x = 0.5x\) is 0.55 at \(x = 1.1\) and 0.45 at \(x = 0.9\), and \(u_y = 0.5y\) is 0.55 at \(y = 1.1\) and 0.45 at \(y = 0.9\). So \((0.55 - 0.45)/0.2 + (0.55 - 0.45)/0.2 = 0.5 + 0.5 = 1\) /s, the exact value. This is the same kind of recipe as Post 2's three-point curvature stencil: a fixed combination of neighbouring values, here for first slopes instead of curvature.

When the density can vary, as in a gas, counting volume no longer works: a box can gain mass without gaining volume. So we count mass instead. Each pair of faces now carries the mass flux, \(\rho u_x\), \(\rho u_y\) and \(\rho u_z\), instead of the plain speed. The same three steps apply to each pair: extra mass flux at the far face as slope × width, times the face area to get kg/s, divided by the box volume to get kg/(m³·s). Summed over the pairs, this is the box's net mass outflow per m³ per second, and by the bookkeeping principle the density inside falls exactly that fast:

Where \(\rho(\mathbf{x}, t)\) now denotes the density field (kg/m³), and \(\nabla \cdot (\rho\mathbf{u}) = \frac{\partial (\rho u_x)}{\partial x} + \frac{\partial (\rho u_y)}{\partial y} + \frac{\partial (\rho u_z)}{\partial z}\) denotes the divergence of the mass flux. This is the continuity equation. With \(\nabla \cdot (\rho\mathbf{u})\) moved to the right side, it reads: the density at a place rises exactly as fast as mass flows in, net. With flow only along \(x\), it is the gas equation from the pipe.

Numbers: let \(\partial (\rho u_x)/\partial x = 0.3\) and \(\partial (\rho u_y)/\partial y = -0.1\) kg/(m³·s), with nothing varying along \(z\). Then \(\nabla \cdot (\rho\mathbf{u}) = 0.2\) and \(\partial \rho/\partial t = -0.2\) kg/(m³·s): more leaves than enters, and the gas thins out. Units: the mass flux is kg/(m²·s), its slope per metre is kg/(m³·s), and \(\partial \rho/\partial t\) is (kg/m³)/s. They match.

For water, \(\rho\) is constant. Because it is constant in time, \(\partial \rho/\partial t = 0\), and that term of the continuity equation drops out:

Because it is constant in space, \(\rho\) can be pulled out of the divergence:

Dividing by \(\rho\) (about 1000 kg/m³, never zero) gives \(\nabla \cdot \mathbf{u} = 0\), the same condition as from the box with four faces.

First, what \(\nabla \cdot \mathbf{u} = 0\) does to a piece of water. Take a dyed square of 0.1 × 0.1 m in a flow that squeezes it along \(x\) at \(\partial u_x/\partial x = -0.8\) /s. By slope × width, its right edge moves 0.8 × 0.1 = 0.08 m/s slower than its left edge, so after 0.1 s the width has shrunk by 0.008 m, to about 0.092 m. Since \(\nabla \cdot \mathbf{u} = 0\) forces \(\partial u_y/\partial y = +0.8\) /s, the same reasoning makes the height grow to about 0.108 m. The area stays 0.092 × 0.108 ≈ 0.01 m² (the small leftover difference shrinks for shorter time intervals). A gas could simply get denser; water cannot. Once the flow along one direction is given, the other is not free.

Every equation so far has had the form "time derivative = …": the heat equation \(\partial u/\partial t = \alpha\, \partial^2 u/\partial x^2\) (Post 2), the two halves of the wave equation \(\partial u/\partial t = v\) and \(\partial v/\partial t = c^2\, \partial^2 u/\partial x^2\) (Post 4), and the continuity equation of a gas, \(\partial \rho/\partial t = -\nabla \cdot (\rho\mathbf{u})\). That form is what lets us Euler-step an equation: the right side, computed from the current snapshot, produces the next one. For the gas:

Where \(\rho^{\,n}\) denotes the density at one point at time step \(n\) (the superscript is an index, not an exponent), \(\Delta t\) denotes the time step (s), and \(\big(-\nabla \cdot (\rho\mathbf{u})\big)^{n}\) denotes the rate computed from the snapshot at step \(n\).

\(\nabla \cdot \mathbf{u} = 0\) has no time derivative, so it cannot produce the next snapshot. It can only test a given one: the spring fails (\(\nabla \cdot \mathbf{u} = +1\) /s), the whirlpool passes. It is a constraint, not an evolution equation. It lost its time derivative because \(\rho\) was the only quantity the continuity equation evolved, and for water we froze \(\rho\). This puts \(\nabla \cdot \mathbf{u} = 0\) next to the Laplace equation of Post 2, the elliptic case in Post 3: no time derivative, and it holds everywhere at once.

The next post gives \(\mathbf{u}\) its evolution equation from Newton's law. Stepping it naively can produce snapshots that fail the test, so some force has to correct every snapshot. That force is pressure. The division of labour: Newton's law says how \(\mathbf{u}\) changes, the constraint says which \(\mathbf{u}\) are allowed, and pressure reconciles the two.

Every real fluid can be compressed. The question is whether this particular flow compresses it noticeably.

The yardstick is how a fluid gets out of the way. Air ahead of a car must move aside, and the message "move" travels as a pressure signal at the speed of sound \(c\), the wave speed of the wave-equation posts: about 340 m/s in air and 1480 m/s in water. A car much slower than \(c\) warns the air early, and it steps aside without being squashed. Near \(c\), the car almost catches its own warning, and the air piles up. The ratio of the two speeds is the Mach number:

Where \(U\) denotes the flow speed (m/s), for example the car's speed relative to the air; \(M\) is unitless. A rule of thumb gives the largest relative density change:

Where \(\Delta \rho\) denotes the largest density change the flow causes (kg/m³), so \(\Delta \rho / \rho\) is a unitless fraction.

In the widget, the rings are the pressure signal leaving the nose. Drag the speed up and watch them bunch up in front of the body as the warning's head start shrinks.

A car at 30 m/s has \(M \approx 0.088\), so the air at its bumper goes from 1.2 to about 1.2047 kg/m³, a change of 0.4 %. The squeeze is strongest at the bumper, where the air almost stops relative to the car, and fades smoothly over a few car lengths. Because \(\tfrac{1}{2}M^2\) grows quadratically, a boat in water changes the density by only 0.002 %, while an airliner at \(M = 0.74\) changes it by about 27 %.

This is a one-time modelling decision, made before choosing the equations. Below \(M \approx 0.3\), about 5 %, we use \(\nabla \cdot \mathbf{u} = 0\) with \(\rho\) exactly constant; above it, the full continuity equation with \(\rho\) as its own evolving field.

We now have the acceleration (Post 10) and the constraint (this post). The next post supplies the forces in \(F = ma\): pressure and viscous drag.