EFFICIENT AND HIGHER-ORDER ACCURATE SPLIT-STEP METHODS FOR GENERALISED NEWTONIAN FLUID FLOW

: In numerous engineering applications, such as polymer or blood flow, the dependence of fluid viscosity on the local shear rate plays an important role. Standard techniques using inf-sup stable finite elements lead to saddle-point systems posing a challenge even for state-of-the-art solvers and preconditioners. Alternatively, projection schemes or time-splitting methods decouple equations for velocity and pressure, resulting in easier to solve linear systems. Although pressure and velocity correction schemes of high-order accuracy are available for Newtonian fluids, the extension to generalised Newtonian fluids is not a trivial task. Herein, we present a split-step scheme based on an explicit-implicit treatment of pressure, viscosity and convection terms, combined with a pressure Poisson equation with fully consistent boundary conditions. Then, using standard equal-order finite elements becomes possible. Stability, flexibility and efficiency of the splitting scheme is showcased in two challenging applications involving aortic aneurysm flow and human phonation.


INTRODUCTION
Various engineering and industrial applications such as automotive design, wind or hydraulic power production, medical devices or synthetics manufacturing share incompressible viscous flow as a central element.The modeling and simulation of fluids has thus been of great interest even before the beginning of computer aided design.More often than not, such fluids are modeled assuming a linear relationship between shear rate and viscous stress via constant viscosity.As it turns out, this modeling assumption may be invalid in various scenarios, blood and polymer flows being practically relevant examples.A vast majority of numerical schemes, however, focuses on Newtonian fluids, neglecting these effects.Depending on the problem and specific flow regime, non-Newtonian characteristics can heavily impact the results obtained and conclusions drawn from them [1,2].The most popular approach to incorporate phenomena such as plug flow or shear thinning/thickening is to consider the viscosity depending on the shear rate, leading to so-called generalised Newtonian or quasi-Newtonian assumptions.
Driven by the ever increasing demand, numerical treatment of the Navier-Stokes equations for incompressible flows have become a staple in modern day computational engineering.But despite the enormous efforts invested, large-scale flow problems still challenge state-of-the-art Book of Extended Abstracts of the 6 th ECCOMAS Young Investigators Conference 7 th -9 th July 2021, Valencia, Spain high-performance computer architectures.The development of new algorithms and methods designed for such problems therefore remains an intense field of research.When employing finite elements, basis functions for velocity and pressure have to be chosen with caution, obeying the Ladyzhenskaya-Babuška-Brezzi (LBB) condition.Some extensions of classical workarounds for Newtonian fluids are readily available, ranging from penalty methods [3,4] to pressure Poisson stabilisation [5] and local pressure projection [2].Some residual-based stabilisations have been proposed [6,7], and we recently presented a novel one [8,9], eliminating spurious pressure boundary layers and poor conservation properties even in lowest-order discretisations.However, preconditioning the arising linear (block-) systems is a critical and often limiting task when developing numerical schemes, despite well-performing algorithms being available [10][11][12].
In view of these challenges, one might prefer projection or split-step schemes decoupling velocity and pressure [13,14], thereby decomposing the system into convection-diffusion, Poisson and simple mass matrix problems.Nonetheless, projection methods suffer from artificial pressure boundary conditions (refer to Guermond et al. [15] for an excellent overview), which often call for corrective measures [16,17].As an alternative, Liu [18] combines explicit treatment of the convective velocity with a pressure Poisson equation (PPE) equipped with consistent boundary conditions.While schemes of similar kind have been applied to challenging incompressible flow problems [19][20][21][22], the extension to the non-Newtonian case is in many aspects challenging.Deteix and Yakoubi [14] proposed the so-called shear rate projection scheme which, despite being accurate and simple, requires LBB-stable velocity-pressure pairs and the solution of an advection-diffusion equation, two Poisson problems and more than ten scalar mass matrix problems per time step.
By contrast, we focus herein on the recent extension of the PPE scheme [18] to generalised Newtonian fluids [23,24].This new framework allows for continuous equal-order finite elements, is higher-order accurate, iteration-free, and consists of an advection-diffusion equation, a single PPE, and two mass matrix solves to recover pressure Dirichlet data and viscosity.We focus on the full-traction variant, additionally including Galerkin least-squares (GLS) stabilisation [25] to counteract dominant convective terms and the popular three-element Windkessel model together with backflow stabilisation.

PROBLEM STATEMENT
As a starting point, let us consider mass and momentum balance equations for an incompressible fluid in Ω ⊂ R d , d = 2, 3 and a time interval from t = 0 to T : with a constant density ρ, velocity u, pressure p, volumetric body force f and viscous stress S.
For generalised Newtonian fluids the viscous stress S computes by where µ(x, t) ∈ R + denotes the variable dynamic viscosity and ∇ s u := 1 /2[∇u + (∇u) ⊤ ] is the symmetric part of the velocity gradient.System (1)-( 2) is supplemented by depending on the shear rate γ is usually formulated in terms of a map η : R + → R + \ {0}: A popular choice in the context of shear-thinning haemodynamics or polymeric flows is the well-established Carreau model [26] η with upper and lower limits η 0 and η ∞ , respectively, and further fitting parameters λ and n ≤ 1. Homogeneous Newtonian fluids are naturally included in this setting, e.g. for n = 1.

TIME-SPLITTING SCHEME
The time-splitting scheme is based on a consistently derived PPE equipped with suitable boundary conditions.So, let us start by taking minus the divergence of Eq. ( 1), to obtain We can further use ∇ • u = 0 and which simplifies to The Dirichlet condition for this auxiliary problem is obtained by dotting the traction boundary condition on Γ N (5) with the unit outward normal vector n: and similarly, dotting the momentum balance equation (1) with gives the Neumann condition for the PPE when restricted to Γ D .For a detailed derivation, the interested reader is referred to our recent work [23], while we herein focus directly on a weak formulation of the split-step scheme.Let us denote the L 2 (Ω) and L 2 (Γ D ) scalar products by ⟨•, •⟩ and ⟨•, •⟩ Γ D , respectively, and start off by multiplying the PPE (9) with a test function q ∈ H 1 (Ω), q| Γ N = 0, integrating by parts and inserting the Neumann boundary condition (11), thereby yielding Book of Extended Abstracts of the 6 th ECCOMAS Young Investigators Conference 7 th -9 th July 2021, Valencia, Spain Table 1: Coefficients α m j and β m j of order m = 2 for BDF and extrapolation [27].
which is rewritten using integration by parts once again as Note however, that second-order derivatives are still present, which we mend via omitting some of the details from [23] for brevity.Other vital ingredients of the split-step scheme are (i) the full decoupling of momentum balance and PPE through explicit treatment of the pressure gradient term in ( 1), (ii) projection of the PPE Dirichlet condition on Γ N (10), (iii) recovering the viscosity µ via an L 2 projection and (iv) improving conservation of mass using divergence damping [20,22].For the time integration, we consider variable time steps ∆t n = t n+1 − t n in higher-order accurate backward differentiation (BDF) and extrapolation formulae (indicated by ⋆ ) with coefficients α m j and β m j given in Tab.1: Then, given solutions from previous time steps, the split-step scheme reads 1. Momentum balance: 2. Project viscosity:

PPE Dirichlet condition:
Recover the continuous Book of Extended Abstracts of the 6 th ECCOMAS Young Investigators Conference 7 th -9 th July 2021, Valencia, Spain 5. Divergence damping: which will be used in the following time step.
Note here, that the Poisson problems ( 15)-( 16) can be combined, solving only one Poisson problem per time step (cf.[18,23]) and for Newtonian fluids, the viscosity projection step ( 14) is skipped.Also, the rheological law is easily replaced by swapping the right-hand side of ( 14) and was actually lumped in our experiments.Per time step, the scheme consists of solving a vector-valued advection-diffusion equation, an L 2 projection on Γ N , a lumped mass matrix and one Poisson problem in the auxiliary variable p := p + φ.

NUMERICAL EXAMPLES
The split-step algorithm ( 13)-( 16) is implemented in the open-source finite element library deal.II [28], using parallel algebraic multigrid (AMG) methods provided by Trilinos' ML package [29] for preconditioning the FGMRES and BiCGStab methods used to solve linear systems corresponding to fluid momentum and PPE, respectively.The versatility and computational performance of the scheme is showcased in two fundamentally different applications in biomechanics, the first being flow through an abdominal aortic aneurysm and the second example considering human phonation.

Abdominal aortic aneurysm
Aneurysms are pathological vessel malformations giving rise to deformed, bulging lumina, altering flow fields and triggering various critical health conditions.A physiological setup is created similar to [30] based on flow data and geometry provided in [31,32].This prototypical segment of the abdominal aorta with length L = 20 cm and inlet/outlet radius R = 1 cm is subject to periodic inflow and outlet pressure p depicted in Fig. 1.Starting from the quiescent state, i.e., u 0 = 0, we prescribe u = (u 1 , 0, 0) ⊤ smoothly ramped by with τ = 0.2 s and a quadratic velocity profile, matching the volumetric flow rate computed by the given mean velocity ū.Concerning the fluid parameters, we set ρ = 1060 kg/m 3 and η 0 = 56 mPas, η ∞ = 3.45 mPas, λ = 3.313 s and n = 0.3568 in (8) according to [33].Further modeling aspects such as three-element Windkessel models, backflow stabilisation and GLS stabilisation are included into the split-step scheme.These extensions, typical for haemodynamic applications, merely modify Neumann data h n+1 or add terms to the momentum equation, and a rigorous introduction is omitted for brevity.Moreover, we define the maximum element CFL and Reynolds numbers as with the number of elements N e and directional element size h i taken as the maximum vertex distance in direction i.Based on (18), we aim for CFL e ≤ 0.  shown in Fig. 2 at t ≈ 4.97 s.Consequently, viscosity spans the whole admissible spectrum η ∞ ≤ µ ≤ η 0 due to large variations in the local shear rate.All linear systems are solved reducing the residual by a factor of 10 −6 , taking the last timestep solution as the initial guess.Doing so, iteration counts for momentum balance (N u ) and PPE (N p ) stay below 20, while the projection of pressure Dirichlet data on Γ N needs a constant of 6 steps only for reaching convergence.Note here, how the former two mildly depend on the flow field as shown in Fig. 3, where we include the inlet velocity ū for reference.Fig. 3 also depicts the adapting time step size together with element CFL and Reynolds numbers, showing time steps decreasing from ≈ 0.015 to ≈ 0.002 shortly after peak inflow, yielding a maximum CFL e of ≈ 0.85 without repeating time steps.CFL e > 1 is admissible in the split-step scheme (cf.Pacheco et al. [23]), but only at the cost of increasing iteration counts in the momentum balance solve.

Human phonation
In a second numerical test, we aim to simulate human phonation, which is the process of vocal folds interacting with air from the lungs, creating the human voice.However, in this preliminary two-dimensional study, the setup inspired by Kniesburges et al. [34]   Starting again from a quiescent state (u 0 = 0) and ramping via ( 17) with τ = 0.01 s, we enforce a quadratic inflow profile.The maximum inlet velocity is prescribed as 80 cm/s, yielding an intraglottal maximum velocity of ūG ≈ 56 m/s and Re = ρū G H G /µ = O(10 3 ) being in the physiological range [34].On the outlet, a zero reference pressure is (approximately) set using h = 0. Regarding the solver settings, we choose an initial ∆t 0 = 10 −4 s, adapt the time step size such that CFL e ≤ 0.8 and reduce the residual by a factor of 10 −8 with the last time step's solution as initial guess.The resulting velocity field is characterised by a strong jet, triggering vortices which in return influence the jet direction.Moreover, the pressure field features fluctuations in the vicinity of the jet as shown in Fig. 5. Low iteration counts result over Comparing to the previous aneurysm example, a slight increase is seen, which is due to a combination of worsened element aspect ratios and higher Reynolds number, but also depends on the more strict convergence criterion.

CONCLUSION
Within this work, a time-splitting scheme suitable for incompressible (generalised) Newtonian fluids has been presented.Momentum and mass balance equations are decoupled using an implicit-explicit treatment of the pressure, viscosity and convection terms.Thus, only an advection-diffusion equation for momentum balance and a PPE with fully consistent boundary conditions are computationally relevant steps.Lower equal-order interpolation of velocity and pressure is also found admissible, while temporal accuracy is determined by suitable BDF and extrapolation formulae.Two challenging examples in biomedical context were tackled, namely, flow through an abdominal aortic aneurysm and human phonation, demonstrating the effectiveness and versatility of the presented approach.

given
Dirichlet data g on Γ D and Neumann data in terms of the full normal traction h on Γ N , where Γ D ∪ Γ N = ∂Ω and Γ D ∩ Γ N = ∅.The rheological law describing the viscosity µ Book of Extended Abstracts of the 6 th ECCOMAS Young Investigators Conference 7 th -9 th July 2021, Valencia, Spain

Figure 2 :
Figure 2: Strong recirculation and viscosity gradients in aneurysm at t ≈ 4.97 s (diastole), selected streamlines (left) and viscosity in cut domain together with selected velocity vectors (right).

Figure 3 :
Figure 3: Iteration counts of momentum balance N u , pressure boundary projection N ζ and PPE N p (left); maximum element CFL and Reynolds numbers and ∆t n (right) with inlet velocity ū for reference.

Figure 4 :
Figure 4: Computational domain for the phonation example with vocal folds in dark grey.

Figure 5 :
Figure 5: Snapshot of the solution at step 12000 (t = 13.576ms) in the fold region: selected velocity streamlines (left) and vectors (right) colored by |u| over pressure p in the background.