Consider the incompressible Navier–Stokes equations for a Newtonian fluid of constant density ρ ρ in a domain Ω Ω ⊂ ⊂ R d ( d = 2 , 3 ) with boundary ∂ ∂ Ω Ω = Γ Γ D ∪ ∪ Γ Γ N , being Γ Γ D and Γ Γ N portions of the boundary where respectively a Dirichlet and a Neumann boundary condition is applied ( Γ Γ D ∩ ∩ Γ Γ N = ∅ ∅ ): { ρ ρ ∂ ∂ u ∂ ∂ t + ρ ρ ( u ⋅ ⋅ ∇ ∇ ) u − − ∇ ∇ ⋅ ⋅ σ σ ( u , p ) = f in Ω Ω × × ( 0 , T ) ∇ ∇ ⋅ ⋅ u = 0 in Ω Ω × × ( 0 , T ) u = g on Γ Γ D × × ( 0 , T ) σ σ ( u , p ) n ^ ^ = h on Γ Γ N × × ( 0 , T ) u ( 0 ) = u 0 in Ω Ω × × { 0 } u is the fluid velocity, p the fluid pressure, f a given forcing term, n ^ ^ the outward directed unit normal vector to Γ Γ N , and σ σ ( u , p ) the viscous stress tensor defined as: σ σ ( u , p ) = − − p I + 2 μ μ ε ε ( u ) . Let μ μ be the dynamic viscosity of the fluid, I the second-order identity tensor and ε ε ( u ) the strain-rate tensor defined as: ε ε ( u ) = 1 2 ( ( ∇ ∇ u ) + ( ∇ ∇ u ) T ) . The functions g and h are given Dirichlet and Neumann boundary data, while u 0 is the initial condition. The first equation is the momentum balance equation, while the second represents the mass conservation, namely the continuity equation. Assuming constant dynamic viscosity, using the vectorial identity ∇ ∇ ⋅ ⋅ ( ∇ ∇ f ) T = ∇ ∇ ( ∇ ∇ ⋅ ⋅ f ) and exploiting mass conservation, the divergence of the total stress tensor in the momentum equation can also be expressed as: ∇ ∇ ⋅ ⋅ σ σ ( u , p ) = ∇ ∇ ⋅ ⋅ ( − − p I + 2 μ μ ε ε ( u ) ) = − − ∇ ∇ p + 2 μ μ ∇ ∇ ⋅ ⋅ ε ε ( u ) = − − ∇ ∇ p + 2 μ μ ∇ ∇ ⋅ ⋅ [ 1 2 ( ( ∇ ∇ u ) + ( ∇ ∇ u ) T ) ] = − − ∇ ∇ p + μ μ ( Δ Δ u + ∇ ∇ ⋅ ⋅ ( ∇ ∇ u ) T ) = − − ∇ ∇ p + μ μ ( Δ Δ u + ∇ ∇ ( ∇ ∇ ⋅ ⋅ u ) ⏟ ⏟ = 0 ) = − − ∇ ∇ p + μ μ Δ Δ u . Moreover, note that the Neumann boundary conditions can be rearranged as: σ σ ( u , p ) n ^ ^ = ( − − p I + 2 μ μ ε ε ( u ) ) n ^ ^ = − − p n ^ ^ + μ μ ∂ ∂ u ∂ ∂ n ^ ^ .