Editing Navier–Stokes equations (section)

==Incompressible flow==
The incompressible momentum Navier–Stokes equation results from the following assumptions on the Cauchy stress tensor:<ref name=Batchelor_142_148>Batchelor (1967) pp. 142–148.</ref>

<ul>
<li>the stress is '''[[Galilean invariance|Galilean invariant]]''': it does not depend directly on the flow velocity, but only on spatial derivatives of the flow velocity. So the stress variable is the tensor gradient <math display="inline">\nabla \mathbf{u}</math>.</li>
<li>the fluid is assumed to be [[isotropic]], as with gases and simple liquids, and consequently <math display="inline">\boldsymbol{\tau}</math> is an isotropic tensor; furthermore, since the deviatoric stress tensor can be expressed in terms of the [[dynamic viscosity]] <math display="inline">\mu</math>:
{{Equation box 1
|indent=:
|title='''Stokes' stress [[constitutive equation]]''' ''(expression used for incompressible elastic solids)''
|equation=:<math>\boldsymbol \tau = 2 \mu \boldsymbol \varepsilon</math>
|cellpadding
|border
|border colour = #0073CF
|background colour=#DCDCDC
}}
where
<math display="block">\boldsymbol{\varepsilon} = \tfrac{1}{2} \left( \mathbf{\nabla u} + \mathbf{\nabla u}^\mathrm{T} \right)</math>
is the rate-of-[[strain tensor]]. So this decomposition can be made explicit as:<ref name=Batchelor_142_148/>
{{Equation box 1
|indent=:
|title='''Stokes's stress constitutive equation''' ''(expression used for incompressible viscous fluids)''
|equation=:<math>\boldsymbol \tau = \mu \left[\nabla\mathbf{u} +  (\nabla\mathbf{u}) ^\mathrm{T}\right]</math>
|cellpadding
|border
|border colour = #0073CF
|background colour=#DCDCDC
}}
</li>
</ul>

This is constitutive equation is also called the '''[[Newtonian fluid#Newtonian law of viscosity|Newtonian law of viscosity]]'''.
Dynamic viscosity {{mvar|μ}} need not be constant – in incompressible flows it can depend on density and on pressure. Any equation that makes explicit one of these [[transport coefficient]] in the [[conservative variable]]s is called an [[equation of state]].<ref name="Batchelor 1967 p. 165">Batchelor (1967) p. 165.</ref>

The divergence of the deviatoric stress in case of uniform viscosity is given by:
<math display="block">\nabla \cdot \boldsymbol \tau = 2 \mu \nabla \cdot \boldsymbol \varepsilon = \mu \nabla \cdot \left( \nabla\mathbf{u} + \nabla\mathbf{u} ^\mathrm{T} \right) = \mu \, \nabla^2 \mathbf{u}</math>
because <math display="inline"> \nabla \cdot \mathbf{u} = 0</math> for an incompressible fluid.

Incompressibility rules out density and pressure waves like sound or [[shock wave]]s, so this simplification is not useful if these phenomena are of interest. The incompressible flow assumption typically holds well with all fluids at low [[Mach number]]s (say up to about Mach 0.3), such as for modelling air winds at normal temperatures.<ref name="Ach">See Acheson (1990).</ref> the incompressible Navier–Stokes equations are best visualized by dividing for the density:<ref>{{Cite journal |last1=Abdulkadirov |first1=Ruslan |last2=Lyakhov |first2=Pavel |date=2022-02-22 |title=Estimates of Mild Solutions of Navier–Stokes Equations in Weak Herz-Type Besov–Morrey Spaces |journal=Mathematics |language=en |volume=10 |issue=5 |page=680 |doi=10.3390/math10050680 |issn=2227-7390|doi-access=free }}</ref>

{{Equation box 1
|indent=:
|title='''[[Incompressible flow|Incompressible]] Navier–Stokes equations''' with uniform viscosity (''convective form'')
|equation=<math>\frac{D \mathbf{u}}{D t} = \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla) \mathbf{u} = \nu \,\nabla^2 \mathbf{u} - \frac{1}{\rho}\nabla p + \frac{1}{\rho} \mathbf{f} </math>
|cellpadding
|border
|border colour = #0073CF
|background colour=#DCDCDC
}}

where <math display="inline">\nu = \frac{\mu}{\rho}</math> is called the [[kinematic viscosity]]. 
By isolating the fluid velocity, one can also state:
{{Equation box 1
|indent=:
|title='''[[Incompressible flow|Incompressible]] Navier–Stokes equations''' with constant viscosity (''alternative convective form'')
|equation=<math> \left(\frac{\partial}{\partial t} + \mathbf{u} \cdot \nabla - \nu \,\nabla^2 \right) \mathbf{u} = - \frac{1}{\rho}\nabla p + \frac{1}{\rho} \mathbf{f}.</math>
|cellpadding
|border
|border colour = #0073CF
|background colour=#DCDCDC
}}

If the density is constant throughout the fluid domain, or, in other words, if all fluid elements have the same density, <math display="inline">\rho</math>, then we have

{{Equation box 1
|indent=:
|title='''[[Incompressible flow|Incompressible]] Navier–Stokes equations with constant density and viscosity''' (''convective form'')
|equation=<math>\frac{D \mathbf{u}}{D t} = \nu \,\nabla^2 \mathbf{u} - \nabla \frac{p}{\rho} + \frac{1}{\rho} \mathbf{f},</math>
|cellpadding
|border
|border colour = #0073CF
|background colour=#DCDCDC
}}
where <math display="inline">p/\rho</math> is called the unit [[pressure head]].

In incompressible flows, the pressure field satisfies the [[Poisson equation]],{{r|landau}}

:<math>\nabla^2 p = - \rho \frac{\partial u_i}{\partial x_k}\frac{\partial u_k}{\partial x_i} = - \rho \frac{\partial^2 u_iu_k}{\partial x_kx_i},</math>

which is obtained by taking the divergence of the momentum equations.

{{hidden
| A laminar flow example |
Velocity profile (laminar flow):
<math display="block">u_x = u(y), \quad u_y = 0, \quad u_{z} = 0</math>
for the {{mvar|x}}-direction, simplify the Navier–Stokes equation:
<math display="block">0 = -\frac{\mathrm{d} P}{\mathrm{d} x} + \mu\left(\frac{\mathrm{d}^2 u}{\mathrm{d} y^2}\right)</math>

Integrate twice to find the velocity profile with boundary conditions {{math|1=''y'' = ''h''}}, {{math|1=''u'' = 0}}, {{math|1=''y'' = −''h''}}, {{math|1=''u'' = 0}}:
<math display="block">u = \frac{1}{2\mu}\frac{\mathrm{d} P}{\mathrm{d} x} y^2 + Ay + B</math>

From this equation, substitute in the two boundary conditions to get two equations:
<math display="block">\begin{align}
0 &= \frac{1}{2 \mu}\frac{\mathrm{d} P}{\mathrm{d} x} h^2 + Ah + B \\
0 &= \frac{1}{2 \mu}\frac{\mathrm{d} P}{\mathrm{d} x} h^2 - Ah + B
\end{align}</math>

Add and solve for {{mvar|B}}:
<math display="block">B = -\frac{1}{2 \mu}\frac{\mathrm{d} P}{\mathrm{d} x} h^2</math>

Substitute and solve for {{mvar|A}}:
<math display="block">A = 0</math>

Finally this gives the velocity profile:
<math display="block">u = \frac{1}{2 \mu}\frac{\mathrm{d} P}{\mathrm{d} x} \left(y^2 - h^2\right)</math>
|style = border: 1px solid lightgray;
|headerstyle = text-align:left
}}

It is well worth observing the meaning of each term (compare to the [[Cauchy momentum equation]]):

<math display="block">\overbrace{ 
    \vphantom{\frac{}{}} 
    \underbrace{
        \frac{\partial \mathbf{u}}{\partial t}
    }_{\text{Variation}} 
    + 
    \underbrace{ 
        \vphantom{\frac{}{}} 
        (\mathbf{u} \cdot \nabla) \mathbf{u}
    }_{\begin{smallmatrix}
        \text{Convective}\\
        \text{acceleration}
    \end{smallmatrix}} 
}^{\text{Inertia (per volume)}} 
= 
\overbrace{ 
    \vphantom{\frac{\partial}{\partial}} 
    \underbrace{ 
        \vphantom{\frac{}{}} 
        -\nabla w
    }_{\begin{smallmatrix}
            \text{Internal}\\
            \text{source}
        \end{smallmatrix}
    } 
    + 
    \underbrace{ 
        \vphantom{\frac{}{}} 
        \nu \nabla^2 \mathbf{u}
    }_{\text{Diffusion}} 
}^{\text{Divergence of stress}} 
+ 
\underbrace{ 
    \vphantom{\frac{}{}} 
    \mathbf{g}
}_{\begin{smallmatrix}
        \text{External}\\
        \text{source}
    \end{smallmatrix}}.</math>

The higher-order term, namely the [[shear stress]] divergence <math display="inline">\nabla \cdot \boldsymbol{\tau}</math>, has simply reduced to the [[vector Laplacian]] term <math display="inline">\mu \nabla^2 \mathbf{u}</math>.<ref>Batchelor (1967) pp. 21 & 147.</ref> This Laplacian term can be interpreted as the difference between the velocity at a point and the mean velocity in a small surrounding volume. This implies that – for a Newtonian fluid – viscosity operates as a ''diffusion of momentum'', in much the same way as the [[heat conduction]]. In fact neglecting the convection term, incompressible Navier–Stokes equations lead to a vector [[diffusion equation]] (namely [[Stokes flow|Stokes equations]]), but in general the convection term is present, so incompressible Navier–Stokes equations belong to the class of [[convection–diffusion equation]]s.

In the usual case of an external field being a [[conservative field]]:
<math display="block"> \mathbf g = - \nabla \varphi </math>
by defining the [[hydraulic head]]:
<math display="block">h \equiv w + \varphi </math>

one can finally condense the whole source in one term, arriving to the incompressible Navier–Stokes equation with conservative external field:
<math display="block">\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla) \mathbf{u} - \nu \, \nabla^2 \mathbf{u} = - \nabla h.</math>

The incompressible Navier–Stokes equations with uniform density and viscosity and conservative external field is the '''fundamental equation of [[hydraulics]]'''. The domain for these equations is commonly a 3 or fewer dimensional [[Euclidean space]], for which an [[orthogonal coordinate]] reference frame is usually set to explicit the system of scalar partial differential equations to be solved. In 3-dimensional orthogonal coordinate systems are 3: [[Cartesian coordinate system|Cartesian]], [[Cylindrical coordinate system|cylindrical]], and [[Spherical coordinate system|spherical]]. Expressing the Navier–Stokes vector equation in Cartesian coordinates is quite straightforward and not much influenced by the number of dimensions of the euclidean space employed, and this is the case also for the first-order terms (like the variation and convection ones) also in non-cartesian orthogonal coordinate systems. But for the higher order terms (the two coming from the divergence of the deviatoric stress that distinguish Navier–Stokes equations from Euler equations) some [[tensor calculus]] is required for deducing an expression in non-cartesian orthogonal coordinate systems.
A special case of the fundamental equation of hydraulics is the [[Bernoulli's equation]].

The incompressible Navier–Stokes equation is composite, the sum of two orthogonal equations,
<math display="block">\begin{align}
\frac{\partial\mathbf{u}}{\partial t} &= \Pi^S\left(-(\mathbf{u}\cdot\nabla)\mathbf{u} + \nu\,\nabla^2\mathbf{u}\right) + \mathbf{f}^S \\
\rho^{-1}\,\nabla p &= \Pi^I\left(-(\mathbf{u}\cdot\nabla)\mathbf{u} + \nu\,\nabla^2\mathbf{u}\right) + \mathbf{f}^I
\end{align}</math>
where <math display="inline">\Pi^S</math> and <math display="inline">\Pi^I</math> are solenoidal and [[Conservative vector field|irrotational]] projection operators satisfying <math display="inline">\Pi^S + \Pi^I = 1</math>, and <math display="inline">\mathbf{f}^S</math> and <math display="inline">\mathbf{f}^I</math> are the non-conservative and conservative parts of the body force. This result follows from the [[Helmholtz decomposition|Helmholtz theorem]] (also known as the fundamental theorem of vector calculus). The first equation is a pressureless governing equation for the velocity, while the second equation for the pressure is a functional of the velocity and is related to the pressure Poisson equation.

The explicit functional form of the projection operator in 3D is found from the Helmholtz Theorem:
<math display="block">\Pi^S\,\mathbf{F}(\mathbf{r}) = \frac{1}{4\pi}\nabla\times\int \frac{\nabla^\prime\times\mathbf{F}(\mathbf{r}')}{|\mathbf{r}-\mathbf{r}'|} \, \mathrm{d} V', \quad \Pi^I = 1-\Pi^S</math>
with a similar structure in 2D. Thus the governing equation is an [[integro-differential equation]] similar to [[Coulomb's law|Coulomb's]] and [[Biot–Savart law|Biot–Savart's law]], not convenient for numerical computation.

An equivalent weak or variational form of the equation, proved to produce the same velocity solution as the Navier–Stokes equation,<ref>
{{citation
 | last = Temam | first = Roger
 | author-link = Roger Temam
 | title = Navier–Stokes Equations, Theory and Numerical Analysis
 | pages = 107–112
 | year = 2001 | publisher = AMS Chelsea
}}</ref> is given by,
<math display="block">\left(\mathbf{w},\frac{\partial\mathbf{u}}{\partial t}\right) = -\bigl(\mathbf{w}, \left(\mathbf{u}\cdot\nabla\right)\mathbf{u}\bigr) - \nu \left(\nabla\mathbf{w}: \nabla\mathbf{u}\right) + \left(\mathbf{w}, \mathbf{f}^S\right)</math>

for divergence-free test functions <math display="inline">\mathbf{w}</math> satisfying appropriate boundary conditions. Here, the projections are accomplished by the orthogonality of the solenoidal and irrotational function spaces. The discrete form of this is eminently suited to finite element computation of divergence-free flow, as we shall see in the next section. There, one will be able to address the question, "How does one specify pressure-driven (Poiseuille) problems with a pressureless governing equation?".

The absence of pressure forces from the governing velocity equation demonstrates that the equation is not a dynamic one, but rather a kinematic equation where the divergence-free condition serves the role of a conservation equation. This would seem to refute the frequent statements that the incompressible pressure enforces the divergence-free condition.

===Weak form of the incompressible Navier–Stokes equations===
====Strong form====
Consider the incompressible Navier–Stokes equations for a [[Newtonian fluid]] of constant density <math display="inline">\rho</math> in a domain
<math display="block"> \Omega \subset \mathbb R^d \quad (d=2, 3)</math>
with boundary
<math display="block"> \partial \Omega = \Gamma_D \cup \Gamma_N ,</math>
being <math display="inline">\Gamma_D</math> and <math display="inline">\Gamma_N</math> portions of the boundary where respectively a [[Dirichlet boundary condition|Dirichlet]] and a [[Neumann boundary condition]]  is applied (<math display="inline">\Gamma_D \cap \Gamma_N = \emptyset</math>):<ref name="Quarteroni">{{cite book |author-link1= Alfio Quarteroni |last1=Quarteroni |first1=Alfio |title=Numerical models for differential problems |publisher=Springer |isbn=978-88-470-5522-3 |edition=Second|date=2014-04-25 }}</ref>
<math display="block">
\begin{cases}
\rho \dfrac{\partial \mathbf{u}}{\partial t} + \rho (\mathbf{u} \cdot \nabla) \mathbf{u} - \nabla \cdot \boldsymbol{\sigma} (\mathbf{u}, p) = \mathbf{f} & \text{ in } \Omega \times (0, T) \\
\nabla \cdot \mathbf{u} = 0  & \text{ in } \Omega \times (0, T) \\
\mathbf{u} = \mathbf{g} & \text{ on } \Gamma_D \times (0, T) \\
 \boldsymbol{\sigma} (\mathbf{u}, p) \hat{\mathbf{n}} = \mathbf{h} & \text{ on } \Gamma_N \times (0, T) \\
\mathbf{u}(0)= \mathbf{u}_0 & \text{ in } \Omega \times \{ 0\}
\end{cases}
</math>
<math display="inline">\mathbf{u}</math> is the fluid velocity, <math display="inline">p</math> the fluid pressure, <math display="inline">\mathbf{f}</math> a given forcing term, <math>\hat{\mathbf{n}}</math> the outward directed unit normal vector to <math display="inline">\Gamma_N</math>, and <math display="inline">\boldsymbol{\sigma}(\mathbf{u}, p)</math> the [[viscous stress tensor]] defined as:<ref name="Quarteroni" /> 
<math display="block">\boldsymbol{\sigma} (\mathbf{u}, p) = -p \mathbf{I} + 2 \mu \boldsymbol{\varepsilon}(\mathbf{u}).</math>
Let <math display="inline">\mu</math> be the dynamic viscosity of the fluid, <math display="inline">\mathbf{I}</math> the second-order [[Identity matrix|identity tensor]] and <math display="inline">\boldsymbol{\varepsilon}(\mathbf{u})</math> the [[strain-rate tensor]] defined as:<ref name="Quarteroni" /> 
<math display="block">
\boldsymbol{\varepsilon} (\mathbf{u}) = \frac{1}{2} \left(\left( \nabla \mathbf{u} \right) + \left( \nabla \mathbf{u} \right)^\mathrm{T}\right).
</math>
The functions <math display="inline">
\mathbf{g}
</math> and <math display="inline">
\mathbf{h}
</math> are given Dirichlet and Neumann boundary data, while <math display="inline">
\mathbf{u}_0
</math> is the [[initial condition]]. The first equation is the momentum balance equation, while the second represents the [[conservation of mass|mass conservation]], namely the [[continuity equation]]. 
Assuming constant dynamic viscosity, using the vectorial identity
<math display="block"> \nabla \cdot \left( \nabla \mathbf{f} \right)^\mathrm{T} = \nabla ( \nabla \cdot \mathbf{f} )</math>
and exploiting mass conservation, the divergence of the total stress tensor in the momentum equation can also be expressed as:<ref name="Quarteroni" /> 
<math display="block">
\begin{align}
\nabla \cdot \boldsymbol{\sigma} (\mathbf{u}, p)
& = \nabla \cdot \left(-p \mathbf{I} + 2 \mu \boldsymbol{\varepsilon}(\mathbf{u}) \right) \\
& = - \nabla p + 2 \mu \nabla \cdot \boldsymbol{\varepsilon}(\mathbf{u}) \\
& = - \nabla p + 2 \mu \nabla \cdot \left [ \tfrac{1}{2} \left(\left(\nabla \mathbf{u} \right) + \left(\nabla \mathbf{u} \right)^\mathrm{T}\right) \right] \\
& = -\nabla p + \mu \left(\Delta \mathbf{u} + \nabla \cdot \left(\nabla \mathbf{u} \right)^\mathrm{T} \right) \\
& =  -\nabla p + \mu \bigl( \Delta \mathbf{u} + \nabla  \underbrace{(\nabla \cdot \mathbf{u})}_{=0} \bigr)
= -\nabla p + \mu \, \Delta \mathbf{u}.
\end{align}
</math>
Moreover, note that the Neumann boundary conditions can be rearranged as:<ref name="Quarteroni" /> 
<math display="block">
 \boldsymbol{\sigma}(\mathbf{u}, p) \hat{\mathbf{n}} = \left(-p \mathbf{I} + 2 \mu \boldsymbol{\varepsilon}(\mathbf{u})\right)\hat{\mathbf{n}} = -p \hat{\mathbf{n}} + \mu \frac{\partial \boldsymbol u}{\partial \hat{\mathbf{n}}}.
</math>

====Weak form====
In order to find the weak form of the Navier–Stokes equations, firstly, consider the momentum equation<ref name="Quarteroni" /> 
<math display="block">
\rho \frac{\partial \mathbf{u}}{\partial t}  - \mu \Delta \mathbf{u} + \rho (\mathbf{u} \cdot \nabla) \mathbf{u} + \nabla p  = \mathbf{f} 
</math>
multiply it for a test function <math display="inline">
\mathbf{v}
</math>, defined in a suitable space <math display="inline">
V
</math>, and integrate both members with respect to the domain <math display="inline">
\Omega
</math>:<ref name="Quarteroni" />
<math display="block">
\int \limits_\Omega \rho \frac{\partial \mathbf{u}}{\partial t}\cdot \mathbf{v}  - \int \limits_\Omega \mu \Delta \mathbf{u} \cdot \mathbf{v} + \int \limits_\Omega \rho (\mathbf{u} \cdot \nabla) \mathbf{u} \cdot \mathbf{v} + \int \limits_\Omega \nabla p \cdot \mathbf{v} = \int \limits_\Omega \mathbf{f} \cdot \mathbf{v}
</math>
Counter-integrating by parts the diffusive and the pressure terms and by using the Gauss' theorem:<ref name="Quarteroni" /> 
<math display="block">\begin{align}
-\int \limits_\Omega \mu \Delta \mathbf{u} \cdot \mathbf{v} &= \int_\Omega \mu \nabla \mathbf{u} \cdot \nabla \mathbf{v} - \int \limits_{\partial \Omega} \mu \frac{\partial \mathbf{u}}{\partial \hat{\mathbf{n}}} \cdot \mathbf{v} \\
\int \limits_\Omega \nabla p \cdot \mathbf{v} &= -\int \limits_\Omega p \nabla \cdot \mathbf{v} + \int \limits_{\partial \Omega} p \mathbf{v} \cdot {\hat{\mathbf{n}}}
\end{align}</math>

Using these relations, one gets:<ref name="Quarteroni" /> 
<math display="block">
\int \limits_\Omega \rho \dfrac{\partial \mathbf{u}}{\partial t}\cdot \mathbf{v}
 + \int \limits_\Omega \mu \nabla \mathbf{u} \cdot \nabla \mathbf{v}
 + \int \limits_\Omega \rho (\mathbf{u} \cdot \nabla) \mathbf{u} \cdot \mathbf{v}
 - \int \limits_\Omega p \nabla \cdot \mathbf{v}
= \int \limits_\Omega \mathbf{f} \cdot \mathbf{v}
 + \int \limits_{\partial \Omega} \left ( \mu \frac{\partial \mathbf{u}}{\partial \hat{\mathbf{n}}} - p \hat{\mathbf{n}}\right) \cdot \mathbf{v}
 \quad \forall \mathbf{v} \in V.
</math>
In the same fashion, the continuity equation is multiplied for a test function {{mvar|q}} belonging to a space <math display="inline">
Q
</math> and integrated in the domain <math display="inline">
\Omega
</math>:<ref name="Quarteroni" />
<math display="block">\int \limits_\Omega q \nabla \cdot \mathbf{u} = 0. \quad \forall q \in Q. </math>
The space functions are chosen as follows: 
<math display="block">\begin{align}
V = \left[H_0^1(\Omega) \right]^d &= \left\{ \mathbf{v} \in \left[H^1(\Omega)\right]^d: \quad \mathbf{v} = \mathbf{0} \text{ on } \Gamma_D \right\}, \\
Q &= L^2(\Omega)
\end{align}</math>
Considering that the test function {{math|'''''v'''''}} vanishes on the Dirichlet boundary and considering the Neumann condition, the integral on the boundary can be rearranged as:<ref name="Quarteroni" /> 
<math display="block">
\int \limits_{\partial \Omega} \left ( \mu \frac{\partial \mathbf{u}}{\partial \hat{\mathbf{n}}} - p \hat{\mathbf{n}} \right) \cdot \mathbf{v} 
= 
\underbrace{
    \int \limits_{\Gamma_D} \left ( \mu \frac{\partial \mathbf{u}}{\partial \hat{\mathbf{n}}} - p \hat{\mathbf{n}} \right) \cdot \mathbf{v}
}_{ \mathbf{v} = \mathbf{0} \text{ on } \Gamma_D \ } 
+ 
\int \limits_{\Gamma_N} \underbrace{ \vphantom{\int \limits_{\Gamma_N} }
    \left ( \mu \frac{\partial \mathbf{u}}{\partial \hat{\mathbf{n}}} - p \hat{\mathbf{n}} \right)
}_{= \mathbf{h} \text{ on } \Gamma_N} \cdot \mathbf{v}
= 
\int \limits_{\Gamma_N} \mathbf{h} \cdot \mathbf{v}.
</math>
Having this in mind, the weak formulation of the Navier–Stokes equations is expressed as:<ref name="Quarteroni" /> 
<math display="block">\begin{align}
&\text{find } \mathbf{u} \in L^2 \left(\mathbb R^+\; \left[H^1(\Omega)\right]^d\right) \cap C^0\left(\mathbb R^+ \; \left[L^2(\Omega)\right]^d\right) \text{ such that: } \\[5pt]
&\quad\begin{cases}
\displaystyle \int \limits_{\Omega}\rho \dfrac{\partial \mathbf{u}}{\partial t}\cdot \mathbf{v} + \int  \limits_{\Omega} \mu \nabla \mathbf{u} \cdot \nabla \mathbf{v} + \int \limits_{\Omega} \rho (\mathbf{u} \cdot \nabla) \mathbf{u} \cdot \mathbf{v} - \int \limits_{\Omega} p \nabla \cdot \mathbf{v} = \int \limits_{\Omega}\mathbf{f} \cdot \mathbf{v} +  \int \limits_{\Gamma_N} \mathbf{h} \cdot \mathbf{v} \quad \forall \mathbf{v} \in V, \\
\displaystyle  \int \limits_{\Omega} q \nabla \cdot \mathbf{u} = 0 \quad \forall q \in Q. 
\end{cases}\end{align}
</math>

===Discrete velocity===
With partitioning of the problem domain and defining [[basis function]]s on the partitioned domain, the discrete form of the governing equation is
<math display="block">\left(\mathbf{w}_i, \frac{\partial\mathbf{u}_j}{\partial t}\right) = -\bigl(\mathbf{w}_i, \left(\mathbf{u}\cdot\nabla\right)\mathbf{u}_j\bigr) - \nu\left(\nabla\mathbf{w}_i: \nabla\mathbf{u}_j\right) + \left(\mathbf{w}_i, \mathbf{f}^S\right).</math>

It is desirable to choose basis functions that reflect the essential feature of incompressible flow – the elements must be divergence-free. While the velocity is the variable of interest, the existence of the stream function or vector potential is necessary by the Helmholtz theorem. Further, to determine fluid flow in the absence of a pressure gradient, one can specify the difference of stream function values across a 2D channel, or the line integral of the tangential component of the vector potential around the channel in 3D, the flow being given by [[Stokes' theorem]]. Discussion will be restricted to 2D in the following.

We further restrict discussion to continuous Hermite finite elements which have at least first-derivative degrees-of-freedom. With this, one can draw a large number of candidate triangular and rectangular elements from the [[Bending of plates|plate-bending]] literature. These elements have derivatives as components of the gradient. In 2D, the gradient and curl of a scalar are clearly orthogonal, given by the expressions,
<math display="block">\begin{align}
\nabla\varphi &= \left(\frac{\partial \varphi}{\partial x},\,\frac{\partial \varphi}{\partial y}\right)^\mathrm{T}, \\[5pt]
\nabla\times\varphi &= \left(\frac{\partial \varphi}{\partial y},\,-\frac{\partial \varphi}{\partial x}\right)^\mathrm{T}.
\end{align}</math>

Adopting continuous plate-bending elements, interchanging the derivative degrees-of-freedom and changing the sign of the appropriate one gives many families of stream function elements.

Taking the curl of the scalar stream function elements gives divergence-free velocity elements.<ref>
{{Citation
 | last = Holdeman | first = J. T.
 | title = A Hermite finite element method for incompressible fluid flow
 | journal = Int. J. Numer. Methods Fluids
 | volume = 64 | pages = 376–408 | year = 2010
 | doi = 10.1002/fld.2154
 | issue = 4 |bibcode = 2010IJNMF..64..376H
| s2cid = 119882803
 }}</ref><ref>
{{Citation
 | last1 = Holdeman | first1 = J. T.
 | last2 = Kim | first2 = J. W.
 | title = Computation of incompressible thermal flows using Hermite finite elements
 | journal = Comput. Meth. Appl. Mech. Eng.
 | volume = 199 | pages = 3297–3304 | year = 2010
 | doi = 10.1016/j.cma.2010.06.036
 | issue = 49–52 |bibcode = 2010CMAME.199.3297H
}}</ref> The requirement that the stream function elements be continuous assures that the normal component of the velocity is continuous across element interfaces, all that is necessary for vanishing divergence on these interfaces.

Boundary conditions are simple to apply. The stream function is constant on no-flow surfaces, with no-slip velocity conditions on surfaces.
Stream function differences across open channels determine the flow. No boundary conditions are necessary on open boundaries, though consistent values may be used with some problems. These are all Dirichlet conditions.

The algebraic equations to be solved are simple to set up, but of course are [[#Nonlinearity|non-linear]], requiring iteration of the linearized equations.

Similar considerations apply to three-dimensions, but extension from 2D is not immediate because of the vector nature of the potential, and there exists no simple relation between the gradient and the curl as was the case in 2D.

===Pressure recovery===
Recovering pressure from the velocity field is easy. The discrete weak equation for the pressure gradient is,
<math display="block">(\mathbf{g}_i, \nabla p) = -\left(\mathbf{g}_i, \left(\mathbf{u}\cdot\nabla\right)\mathbf{u}_j\right) - \nu\left(\nabla\mathbf{g}_i: \nabla\mathbf{u}_j\right) + \left(\mathbf{g}_i, \mathbf{f}^I\right)</math>

where the test/weight functions are irrotational. Any conforming scalar finite element may be used. However, the pressure gradient field may also be of interest. In this case, one can use scalar Hermite elements for the pressure. For the test/weight functions <math display="inline">\mathbf{g}_i
</math> one would choose the irrotational vector elements obtained from the gradient of the pressure element.