Convergence rate of solution scheme for incompressible Navier-Stokes equations

[Marc Fehling]

Dear community,

as the question title suggests, I'm having trouble verifying the convergence order of a solution scheme for the incompressible Navier-Stokes equations and I'm addressing this question particularly to those who are familiar with the verification of such a scheme.

We're currently working on a numerical solution scheme for the incompressible Navier-Stokes equations using Chorin's projection. Currently, we're using continuous Taylor-Hood elements. As mentioned in tutorial step-12, advection-like problems are not stable with CG methods, thus we're applying external stabilization, i.e. Taylor-Galerkin and grad-div-stabilization. The time marching scheme is chosen to be a semi-implicit scheme, basing on the implicit Euler scheme. We call it semi-implicit because we linearize the advection term, replacing the (u^* * nabla) u^* by (u^n * nabla) u^*.

Now, I want to verify the solution scheme with a convergence analysis for the flow velocity. I take the L2 error using the "integrate_difference(...)" function and compare the different values depending on the size of the time-step and the global refinement level. To get an error indicator, I take the root mean square over the L2 error at every timestep. The used function for verification is the non-trivial, two-dimensional solution to the incompressible Navier-Stokes (INS) equations made by McDermott (source). Since this function solves the INS equations intrinsically, no modification of the right hand side of the equations is made. Periodic boundaries are used to avoid the external imposition of boundary conditions and thus another error source.

If I run the simulation with different timesteps at a specified global level, I can reproduce the expected order in time, which is one, resulting from the implicit Euler scheme. But if I run the simulation at different refinement levels with a fixed timestep, I get convergence rates which are not consistent. Using second order elements for the velocity space, I would expect a convergence order of 3 for the L2 error. But the convergence rates I get are jumping wildly and seem to depend on the size of the chosen timestep. As an exmaple using small timesteps (dt=1e-4, dx<0.2, cfl<1e-3) I get convergence rates of roughly 3+/-0.6 between two global levels. With smaller timesteps (dt=1e-3), I get rates around 1.6+/-0.2. If I take a look at the convergence rates at different timepoints (not the root mean square one) for the smallest timestep (dt=1e-4), I see that convergence rates in the beginning are indeed the expected ones of 3, but are changing over time.

Why is the convergence rate in space inconsistent? Am I missing some crucial point?

Best regards,
Marc Fehling