# Projection (or fractional-step) methods Vs coupled method for incompressible Navier-Stokes

My question is in the context of the finite element method.

Incompressible Navier-Stokes equations can be solved using the coupled method or projection/fractional step methods. Each method has its advantages and disadvantages. Projection methods are said to be cost-effective because of smaller matrix systems; however, it is quite difficult to achieve higher-order (even second-order) temporal accuracies in the case of projection methods.

So, I am curious to know if there are comparisons in the literature between the two in terms of accuracy and computational cost for transient problems. I have not found any such literature so far. I would appreciate it if you could share your experience and/or share any relevant papers/reports showing the comparison in terms of cost and accuracy.

• I am pretty sure one could have discussions of religious fervor about this issue. Have you looked at Volker John's book, for example? Commented Apr 23, 2023 at 22:21
• @WolfgangBangerth Yeah. That book discusses the costs briefly. No data on computational cost is provided to understand it in detail. Commented Apr 23, 2023 at 22:43
• @ChennaK Some papers of Guermond et al. and Fehn et al. come into my mind + literature therein. You may also have a look at a third candidate, the artificial compressibilty method. Commented Apr 24, 2023 at 18:14
• @ConvexHull, Wouldn't dual time stepping needed for the artificial compressibility method make it more expensive? Commented Apr 25, 2023 at 9:29
• We tested the FEM implementation of several algorithms (including projection method) for solving transient 3D Navier-Stokes. As a reference solution we used exact solution from the paper onlinelibrary.wiley.com/doi/abs/10.1002/fld.1650190502 Commented Apr 27, 2023 at 15:35

Let us consider the system of equations describing viscous incompressible flow in a unit cube $$0\le x\le 1, 0\le y \le 1, 0 \le z \le 1$$, and in time interval $$0\le t \le 1$$, we have $$$$\label{NSE} \begin{array}{cc} \nabla.\bf{u}=0\\ \frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}.\nabla)\mathbf{u}+\frac{\nabla p}{\rho}=\nu\nabla ^2\mathbf{u} \end{array}$$$$ Here it is indicated:$$\rho=1$$ is density, $$\mathbf{u}=(u,v,w)$$ - flow velocity, $$\nu$$ - kinematic viscosity, $$p$$ - pressure. We put Dirichlet conditions on the boundary of cube in a form $$$$\label{NSE_bc} \mathbf{u}(x,y,z,t)=\mathbf{u}_e(x,y,z,t), p(x,y,z,t)=p_e(x,y,z,t)$$$$ and initial condition in the cube as follows $$$$\label{NSE_ic} \mathbf{u}(x,y,z,0)=\mathbf{u}_e(x,y,z,0), p(x,y,z,0)=p_e(x,y,z,0)$$$$ where $$\mathbf{u}_e=(u_e,v_e,w_e), p_e$$ is an exact (benchmark) solution described in the paper Exact fully 3D Navier–Stokes solutions for benchmarking, and given by $$$$\label{benchmarc_solution} \begin{array}{cc} u_e=-a \exp(-d^2 t) (\exp(a x) \sin(a y + d z) + \exp(a z) \cos(a x + d y)),\\ v_e =-a \exp(-d^2 t) (\exp(a y) \sin(a z + d x) + \exp(a x) \cos(a y + d z)),\\ w_e = -a \exp(-d^2 t) (\exp(a z) \sin(a x + d y) + \exp(a y) \cos(a z + d x)),\\ p_e = -\frac{a^2}{2}\exp(-2 d^2 t) (\exp(2 a x) + \exp(2 a y) + \exp(2 a z) +\\ 2 \sin(a x + d y) \exp(a (y + z)) \cos(a z + d x) +\\ 2 \sin(a y + d z) \exp(a (x + z)) \cos(a x + d y) +\\ 2 \sin(a z + d x) \exp(a (y + x)) \cos(a y + d z)) \end{array}$$$$ here $$a, d$$ are free parameters, and we put in numerical calculations $$a=d=1$$. To solve this problem we use nonlinear FEM numerical algorithm with an implicit time step $$$$\label{NSE_num} \begin{array}{cc} \nabla.\bf{u}(\mathbf{r},t)=0\\ \frac{\mathbf{u}(\mathbf{r},t)-\mathbf{u}(\mathbf{r},t-\tau)}{\tau}+(\mathbf{u}(\mathbf{r},t).\nabla)\mathbf{u}(\mathbf{r},t)+\nabla p(\mathbf{r},t)=\nu\nabla ^2\mathbf{u}(\mathbf{r},t) \end{array}$$$$ Here $$\tau$$ is a time step.
The projection method based on the system of equation $$$$\label{NSE_project} \begin{array}{cc} \frac{\mathbf{u}-\mathbf{u}_n}{\tau}+(\mathbf{u}_n.\nabla)\mathbf{u}=\nu\nabla ^2\mathbf{u}\\ \tau \nabla^2p-\nabla.\mathbf{u}=0,\\ \frac{u_{n+1}-u}{\tau}+\nabla p=0 \end{array}$$$$ here $$u_n,u,u_{n+1}$$ is velocity field on previous, intermediate and next step consequently, and p is a pressure. This linear system can be solved with linear FEM.
We can compare numerical solution with exact solution shown above at $$\nu =1$$. For this we use mesh of 1000 Hexahedron Element in the unit cub. In the first test we take $$\tau=1/40$$. Maximal absolute error of numerical solution $$\Delta u, \Delta v, \Delta w, \Delta p$$ computed with nonlinear FEM at 40 steps shown below for different $$z=z_i, i=1,2,...,8$$
Second test we run with $$\tau=1/40, 1/200, 1/400, 1/800$$ (first line, second line, third line, fourth line in the picture below) . Maximal absolute error of numerical solution $$\Delta u, \Delta v, \Delta w, \Delta p$$ computed with projection step at different $$z=z_i, i=1,2,...,8$$ Note, the velocity and pressure error decreases with $$\tau$$ decreasing, but for $$\tau =1/40$$ solution with projection step is not so precise as computed with nonlinear FEM. Computational time for the test with nonlinear FEM is about 267s (on my notebook), for projection step with $$\tau=1/40$$ the computational time is about 198s. All tests been implemented with Mathematica FEM. We can also implement projection method with using FFT - see my post here.