Higher-Order Discontinuous Galerkin Splitting Schemes for Fluids with Variable Viscosity

TL;DR

This work presents a matrix-free, higher-order discontinuous Galerkin discretization for the incompressible Navier–Stokes equations with variable viscosity, using SIPG for the viscous term. It compares monolithic (coupled) time stepping with projection-based splitting schemes, exploring multiple linearization strategies and preconditioners on hp-multigrid hierarchies. Numerical results demonstrate that projection/splitting approaches often deliver substantial throughput gains (roughly 1.3–3.5x) over monolithic solves across a range of viscosity contrasts and Reynolds numbers, while maintaining accuracy; however, fully implicit schemes can be more robust in some regimes. The findings underscore the value of splitting methods for variable-viscosity flows, especially when robust preconditioning is challenging, and point to the potential of improved preconditioners to narrow gaps with coupled solvers; BDF-3 stability issues advise sticking to lower-order implicit time integrators in practice.

Abstract

This article investigates matrix-free higher-order discontinuous Galerkin discretizations of the Navier--Stokes equations for incompressible flows with variable viscosity. The viscosity field may be prescribed analytically or governed by a rheological law, as often found in biomedical or industrial applications. The DG discretization of the adapted second-order viscous terms is carried out via the symmetric interior penalty Galerkin method, obviating auxiliary variables. Based on this spatial discretization, we compare several linearized variants of saddle point block systems and projection-based splitting time integration schemes in terms of their computational performance. Compared to the velocity-pressure block-system for the former, the splitting scheme allows solving a sequence of simple problems such as mass, convection-diffusion and Poisson equations. We investigate under which conditions the improved temporal stability of fully implicit schemes and resulting expensive nonlinear solves outperform the splitting schemes and linearized variants that are stable under hyperbolic time step restrictions. The key aspects of this work are i) a higher-order DG discretization for incompressible flows with variable viscosity, ii) accelerated nonlinear solver variants and suitable linearizations adopting a matrix-free $hp$-multigrid solver, and iii) a detailed comparison of the monolithic and projection-based solvers in terms of their (non-)linear solver performance. The presented schemes are evaluated in a series of numerical examples verifying their spatial and temporal accuracy, and the preconditioner performance under increasing viscosity contrasts, while their efficiency is showcased in the backward-facing step benchmark.

Figures (10)

• Figure 1: Spatial convergence study employing inf-sup stable finite element pairs with velocity degree $k_u$ and pressure degree $k_p = k_u-1$ in the projection and coupled solvers. Optimal orders of accuracy in the relative $L^2$ errors are achieved up to linear solver tolerance.
• Figure 2: Temporal convergence study employing uniform time steps and BDF-$m$ schemes $m=1,2,3$. Optimal orders of accuracy in the relative $L^2$ errors are achieved with the coupled solver. The projection solver yields identical results up to $m=3$. The known unstable behavior is also observed for the Newtonian splitting scheme ($\dag$) and is found to be more drastic for lower linear solver tolerance ($*$, $\epsilon_\mathrm{lin}^\mathrm{abs} = 10^{-16}$, $\epsilon_\mathrm{lin}^\mathrm{rel} = 10^{-8}$).
• Figure 3: Temporal convergence study employing uniform time steps and BDF-$2$ schemes under variation of the upper viscosity limit $\eta_0$, keeping $\eta_\infty=5\times10^{-2}~\text{m}^2/\text{s}$ fixed. Optimal orders of accuracy in the relative $L^2$ errors are achieved with the coupled solver. The projection solver yields optimal rates for the Newtonian case, but reduces earlier to linear convergence for high viscosity contrasts.
• Figure 4: Exemplary steady-state solution on the $xy$ symmetry plane of the cavity using $\eta_0 = 100 \, \eta_\infty = 10^{-1}~\mathrm{m}^2\mathrm{/s}$.
• Figure 5: Velocity magnitude $\text{\boldmath${u}$}$ at $t = 0.2, 0.4, 0.6, 0.8$ and $1.0$ s for $\mathrm{Re}=300$ in the backward-facing step benchmark with smoothly ramped inlet velocity. Legend applies to all solution snapshots.