The equations of motion for a viscous fluid are called the*Navier-Stokes*equations.

We will derive these equations from Equations (2.25) and (2.27) in the following,
but before we do so we will make a further simplifying assumption. The fluids
of interest in computer graphics, such as smoke and liquids at speed far below
the speed of sound do not change their volumes easily, so it convenient to make
an assumption of*incompressible flow*. Incompressible flow is any flow that keeps
volumes constant,*i.e.* dV/dt =0. It follows from Reynold’s Transport theorem
that

which, following similar reasoning as used in the previous section, implies that

∇ ·**u**=0 (2.28)

which is called the *incompressibility constraint* and is the first equation in the
*incompressible*Navier-Stokes equations. Equation (2.28) is equivalent to assuming
that density is advected passively with the flow, *i.e.* D*ρ/*Dt =0, which is easily
seen from Equation (2.25).

To derive the other equation of incompressible Navier-Stokes, we need to specify a material model (also called a constitutive model). First, we will split the stress tensor into two parts

*σ*=−*pδ*+*τ* (2.29)

where*δ*is the identity tensor. We will use the pressure*p* to enforce the
incom-pressibility constraint and the viscous stress tensor*τ*to model viscosity. We will
use a linear stress-strain relation for*τ*

*τ*=2*µ"*+*λ*tr(")δ.

The dynamic viscosity coefficient *µ* controls the resistance to shearing and the
bulk viscosity coefficient*λ*controls viscosity due to compression and dilation. The
infinitesimal strain tensor*"*=^{1}_{2}(∇u+∇u^{T})is as we saw in Section 2.4. Intuitively,
the tracetr(")measures the local volume change, which is zero for incompressible
flow,*i.e.*tr(") =∇ ·**u**=0, so it is customary to set*λ*=0. If we further assume
that*µ*is constant in space and plug*σ*into Equation (2.27) we obtain

*ρ*Du

Dt +∇*p*−*µ∇*^{2}**u**=**f**

which is the*momentum equation*of Navier-Stokes. If we further expand the
ma-terial derivative D/Dt, divide through by*ρ*(recalling that**f**=*ρ***g**) and define the
kinematic viscosity*ν*=*µ/ρ* we get

*∂***u**

*∂t* + (**u**· ∇)**u**+ 1

*ρ*∇*p*−*ν∇*^{2}**u**=**g.** (2.30)
which is the form most commonly seen in computer graphics literature.

Equations (2.28) and (2.30) together form the incompressible Navier-Stokes equations. However, we still need to relate the pressure in the momentum equation to the incompressibility constraint, which we will do in the next section.

**2.7.1** **Pressure**

The pressure appearing in the momentum equation (Equation (2.30)) is somewhat
strange in the case of incompressible flow. The pressure itself is *not* a physical
quantity, but the gradient of pressure∇*p*is a force. In Chapter 4 we will be using
the pressure gradient to define an energy that measures deviations of a liquid
surface from the set of physically valid states, so it makes sense to spend a little
time here to try and understand pressure.

First, lets derive an expression for pressure. We will assume the
incompress-ibility constraint∇ ·**u**=0and see what we get. By taking the divergence of the
momentum equation we get

∇ · 1

*ρ*∇*p*=∇ · −*∂***u**/∂*t*−(**u**· ∇)**u**+*ν∇*^{2}**u**+**g.**

(2.31) If we now use the incompressibility assumption so (assuming derivatives commute)

∇ · ^{∂}_{∂}^{u}*t*

= _{∂}^{∂}*t*(∇ ·**u**) =0we obtain

∇ · 1

*ρ*∇*p*=∇ · −(**u**· ∇)**u**+*ν∇*^{2}**u**+**g**

(2.32) which is a (variable coefficient) Poisson equation.

We will now show that the pressure in Equation (2.32) *enforces* the
incom-pressibility constraint. Without assuming∇ ·**u**=0we substitute the expression
for pressure in Equation (2.32) into Equation (2.31) to obtain (after commuting
derivatives)

*∂*

*∂t*(∇ ·**u**) =0
If we integrate this equation in time we get

∇ ·**u**=*g*(**x**)

which is independent of time. If we now assume that the initial divergence is
zero this implies that *g*(**x**) =0for all time. This shows that the pressure in
Equa-tion (2.32) does indeed imply that incompressibility is enforced. In Appendix A
we show how Equation (2.32) can be solved numerically using methods from
Sec-tion 2.3 and be used to enforce the incompressibility constraint.

What we have shown is section is that*any* divergence-free velocity field
in-duces a pressure through the incompressible Navier-Stokes equations and that this
pressure ensures that the accelerations in Navier-Stokes are divergence-free
*guar-anteeing* that the velocity field remains divergence-free for*all* time. See Gresho
and Sani [1987] for a more detailed discussion.

**2.7.2** **Helmholtz-Hodge Decomposition Theorem**

Another way of looking at pressure is as a projection onto the set of divergence-free vector fields through the Helmholtz-Hodge decomposition.

If we let*Ω*be a subset ofR^{3}with smooth boundary*∂Ω* then the
Helmholtz-Hodge decomposition theorem states that*any*smooth vector field**u**:*Ω*→R^{3}can
be decomposed*uniquely*as

**u**=∇φ+∇×* Ψ*+

*(2.33)*

**γ**where*φ* :*Ω*→R is a scalar field, * Ψ* :

*Ω*→R

^{3}is a

*divergence-free*

^{1}vector field (meaning∇ ·

*=0) and*

**Ψ***:*

**γ***Ω*→R

^{3}is a

*harmonic*vector-valued field (meaning

∇ ·* γ*=0and∇×

*=0) [Cantarella et al. 2002].*

**γ**To extract the divergence-free part of**u**, let us take the divergence on both sides
of Equation (2.33)

∇ ·**u**=∇ · ∇φ=∇^{2}*φ*
which gives us a Poisson equation for*φ.*

We can now define a projection operatorP(**u**) =**u**− ∇φ, which extracts the
divergence-free part. Clearly,∇ ·P(**u**) =0andP(P(**u**)) =P(**u**)so Pis a proper
projection. Notice that these steps correspond exactly to the steps we took in
Section 2.7.1 with pressure taking the place of *φ* and the momentum equation
taking the place of**u**.

Alternatively, we can take the curl of Equation (2.33) and get

∇×**u**=∇×∇×**Ψ**

=∇:0
(∇ ·* Ψ)*− ∇

^{2}

**Ψ**=−∇^{2}**Ψ**

(2.34)

where the second equality follows from a vector calculus identity and the third
equality from the fact that ∇ ·* Ψ* = 0 in the decomposition. Notice that
Equa-tion (2.34) is a vector-Poisson equaEqua-tion for

**Ψ.**For*Ω*with simple topology, where we may assume that* γ*=0, we can again
define a projection operatorP(

**u**) =∇×

*to recover the divergence-free part of*

**Ψ****u**. Even when

*Ω*has non-trivial topology,

*is often small for the kinds of vector fields that arise from fluid flows [Tong et al. 2003]. This means that in many fluid simulation papers in computer graphics (*

**γ***e.g.*Ando et al. [2015]) the authors apply this method even when

*Ω*has non-trivial topology.

**2.7.3** **Boundary conditions**

So far we have only talked about what happens in the interior of the fluid and have carefully avoided any mention of boundaries. In this section we will answer which kinds of boundaries exists between a fluid and its surrounding environment, and what happens at these boundaries. This gives us a chance to talk about two important concept, namely surface tension and free surfaces.

In the following we will denote by*∂Ω*the boundary of the fluid. There are two
different kinds of boundaries to consider. First, there are *solid-fluid*boundaries

*∂Ω*Swhere the fluid is in contact with a rigid solid. Second, there are*fluid-fluid*

1Divergence-free velocity fields are also sometimes called*solenoidal*vector fields.

boundaries *∂Ω*F where two fluids are in contact. The two kinds of boundaries
together form the whole boundary of the fluid*∂Ω*=*∂Ω*S∪*∂Ω*F.

For the purposes of this thesis we will be assuming that when multiple fluids
interact, that the density of one of the fluids is significantly higher than the density
of the other fluids. For instance, this is the case of water and air that have a density
ratio of about one thousand. Under this assumption we may neglect the influence
of the lighter fluids on the motion of the heavier fluid. This approximation is
called the*free surface approximation*^{2}and is prevalent in computer graphics, since
it allows us to model only a single fluid with significant computational savings as
a result.

Even when we assume that fluid-fluid interfaces are free surfaces, there is
another effect that becomes important at such interfaces at small length scales and
that is*surface tension*. Surface tension is a cohesive force resulting from the fact
that fluid particles prefer to stick together with other fluid particles from the same
fluid. Surface tensions causes a fluid to tend towards reducing its surface area,
which causes extra tension at the surface.

We will not dwell too much on topic, but as was the case when we derived the equations of motion for the interior of a continuum in Section 2.6, what happens at the boundary is also given by a force balance. In this case the net force on the area element of the surface should be balanced which eventually leads to the following boundary condition (see Subramanian [2015] and Bridson [2008] for more details) for the stress

(σ1−*σ*2)**nˆ**=2H*γ***nˆ** on *∂Ω*F

where*σ*1and*σ*2are the stresses for the two fluids that form the interface,**nˆ** is
the surface normal to the interface and*H* is the mean curvature of the interface
(we will see how to compute this in Section 2.10). Surface tension is denoted by*γ*
and has units of force per unit length. For the interface between water and air at
room temperature it is approximately*γ*=0.073newton per meter. Recalling from
Equation (2.29) that the stress was defined as*σ*=−*pδ*+*τ, this means that for an*
*inviscid*fluid with no viscosity (*i.e.τ*=0) we have

(*p*_{1}−*p*_{2}) =2H*γ* on *∂Ω*F

since pressure only works in the normal direction (it was defined as a diagonal stress tensor). Under the free surface approximation we may simplify this further.

We could for instance set*p*_{2}to be the ambient air pressure. In fact, since pressure
is only ever evaluated as a gradient, constant offsets do not matter and we may
equivalently set*p*_{2}=0so (defining*p*:=*p*_{1})

*p*=2H*γ* on *∂Ω*_{F}

2Technically, a free surface is a surface subject to constant normal stress and zero shear stress. This prevents effects like the wind creating waves in the ocean, which is a shear effect.

Finally, in case we are mainly interested in large scales (*e.g.*if we are modelling
the ocean) we may neglect surface tension completely, which gives us

*p*=0 on *∂Ω*_{F}

At solid-fluid boundaries we certainly do not want fluid to flow into or out of the solid. This leads to the following boundary condition

**u**·**nˆ**=**u**_{solid}·**nˆ** on *∂Ω*_{S}.

where where**u**_{solid}is the velocity of the solid and**nˆ** is the normal vector to solid
obstacle. This boundary condition is called*free-slip*, since the fluid is unrestricted
in its movement tangentially to the solid. This is the correct boundary condition
for an inviscid fluid. For a viscous fluid it turns out that the correct boundary
condition is

**u**=**u**_{solid} on *∂Ω*S.

This boundary condition is called*no-slip*, since it forces the fluid to have the same
velocity as the solid at the boundary. Nevertheless, it is common for practioners
to use the free-slip boundary condition even for viscous fluids, since the no-slip
boundary condition, although accurate in the limit, is too inaccurate for the kinds
of coarse grids used in computer graphics.

**2.7.4** **Conservation form**

In Chapter 5 we will be using an alternative but equivalent version of
Navier-Stokes for convenience. This form of Navier-Navier-Stokes is called*conservation form*and
emphasizes the fact that quantities like momentum and mass are conserved.

*∂***q**
Here, **u** denotes the usual three-dimensional velocity and *σ* is Cauchy’s strain
tensor as defined in Equation (2.29). Notice that we include both the momentum
equation and the incompressibility constraint into*one*equation.

To see why Equation (2.35) is called conservation form, let**f**=**0**, integrate both
sides over a fixed volume*Ω*and apply the divergence theorem to obtain

*∂*

where**nˆ**is the outwards pointing normal vector to*∂Ω. This says that the change*
of**q**in time inside*Ω*is purely due to the flux**F**(**q**)through*∂Ω. SinceΩ*is arbitrary

this implies that**q**is (locally) conserved. Note, this is only true in the interior of the
fluid and separate zero-flux boundary conditions have to be applied to ensure that
the equations form a closed system so that the quantities are conserved globally.