Consider the Navier stokes equation after the discretization with conforming finite elements with source term $f=0$. We have the algebraic structure of a saddle point problem:

$$M \dot{u} = f- Au -B^tp - C(u(t))u(t)$$ $$Bu=0$$

where $B$ is the $- \operatorname{div}$ matrix, $A$ is the one corresponding to the laplacian and $C$ corresponds to the trilinear form of the equation.

I've read that it's possible to write this system as an ODE in $(\mathbf{u},p)$, instead of a plain DAE. Basically the incompressibility constraint $Bu=0$ will be absorbed in the resulting system of ODEs. The big problem is that I can't see how this may be done! It should be a simple computation, but I've spent the whole day trying it and I am not able to find that out.

  • $\begingroup$ There is no time derivative of pressure in the second equation. Unless you modify the second equation by employing the artificial compressibility method, I don't see how it is possible to write the system as ODEs. $\endgroup$
    – Chenna K
    Jul 28, 2021 at 16:38
  • $\begingroup$ In my notes, the second equation of the system of ODEs (i.e. the one for $p$) reads $M \dot{p} =-B u$... so I think he adopted that approach. So the resulting system is ODEs is $$M \dot{u} = f- Au -B^tp - C(u(t))u(t)$$ $$M \dot{p} = -Bu$$, right? It seems kinda strange to me honestly. How can one derive such a system? @ChennaK $\endgroup$
    – FEGirl
    Jul 28, 2021 at 17:23
  • 1
    $\begingroup$ Writing it as $\kappa M\dot p = - Bu$ amounts to making the material slightly compressible, with a compressibility of $\kappa$. This is a commonly used trick, but it is not equivalent to the original model. $\endgroup$ Jul 28, 2021 at 17:53
  • $\begingroup$ Using the artificial compressibility method, the incompressibility constraint equation is modified as $\frac{\partial p}{\partial t} + \beta^2 \, \nabla \cdot \mathbf{u} = 0$, where $\beta$ is the weight factor, aka artificial compressibiity parameter. As $\beta \rightarrow \infty$, $\nabla \cdot \mathbf{u} \rightarrow 0$. For practical use large values of $\beta (>100)$ should be sufficient. Keep in mind that this becomes problematic for explicit time integration schemes. $\endgroup$
    – Chenna K
    Jul 28, 2021 at 20:04
  • $\begingroup$ Thanks to both of you for the comments! Could you please point me to some reference on the topic ? A beginner level would be better, as I am a novice on cfd. $\endgroup$
    – FEGirl
    Jul 28, 2021 at 22:31

1 Answer 1


One approach to convert this into an ODE is with index reduction methods. These allow you to convert high-index DAEs into low-index DAEs or ODEs. See section VII.2 of "Solving Ordinary Differential Equations II" by Hairer and Wanner.

Consider a generic, Hessenberg index-2 DAE $$ \begin{align} y' &= f(y, z) \\ 0 &= g(y) \end{align} $$ Indeed, your NS equation fits in this formulation. Differentiating the algebraic constraint with respect to time gives $$ 0 = g_y(y) f(y, z) $$ If we differentiate this again, we finally recover an ODE in z (pressure) $$ z' = -(g_y(y) f_z(y, z))^{-1} \left( g_{yy}(f(y, z), f(y, z)) + g_y(y) f_y(y, z) f(y, z) \right) $$ The obvious downside to this is you need to compute partial derivatives. Since $g(y)$ is linear in NS, $g_{yy}(y) = 0$, which helps to simplify this a bit.

If we apply an ODE solver to our new ODE in terms of $y'$ and $z'$, it is not guaranteed that $g(y) = 0$ or the "hidden" constraint $0 = g_y(y) f(y, z)$ is satisfied. This is known as the "drift-off" phenomenon. These can be addressed with stabilization techniques or by projecting onto a feasible solution manifold.

In summary, it is possible to convert your problem into an ODE, but index reduction has several limitations. I am not an expert in NS, so I can't say if practitioners commonly use this approach, but I do know there are more direct ways to integrate NS.

Edit: For completeness, here would be the reduced equations $$ \begin{align} M \dot{u} &= f- Au -B^tp - \mathcal{N}(u), \\ B M^{-1} B^T \dot{p} &= B M^{-1} (A + \mathcal{N}_u(u)) M^{-1} (f- Au -B^tp - \mathcal{N}(u)) \end{align} $$ where $\mathcal{N}(u) = C(u) u$.

  • $\begingroup$ With the need to compute the derivatives of the constraint, this approach would be impractical for the numerical solutions of incompressible NS equations. $\endgroup$
    – Chenna K
    Jul 28, 2021 at 22:06
  • 2
    $\begingroup$ I agree this is not particularly practical. It's the price to pay for having an ODE. The Jacobian of the nonlinear term might be problematic to compute, and there are significantly more linear solves involving the mass matrix. Getting the derivative of the constrain isn't so bad, though: $g_y(y) = B$. $\endgroup$ Jul 28, 2021 at 22:40
  • 1
    $\begingroup$ Thanks Steven for your answer :-) @StevenRoberts $\endgroup$
    – FEGirl
    Aug 5, 2021 at 18:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.