Predicting Non-Linear Flow Phenomena through Different Characteristics-Based Schemes

The present work investigates the bifurcation properties of the Navier–Stokes equations using characteristics-based schemes and Riemann solvers to test their suitability to predict non-linear flow phenomena encountered in aerospace applications. We make use of a singleand multi-directional characteristics-based scheme and Rusanov’s Riemann solver to treat the convective term through a Godunov-type method. We use the Artificial Compressibility (AC) method and a unified Fractional-Step, Artificial Compressibility with Pressure-Projection (FSAC-PP) method for all considered schemes in a channel with a sudden expansion which provides highly non-linear flow features at low Reynolds numbers that produces a non-symmetrical flow field. Using the AC method, our results show that the multi-directional characteristics-based scheme is capable of predicting these phenomena while the single-directional counterpart does not predict the correct flow field. Both schemes and also Riemann solver approaches produce accurate results when the FSAC-PP method is used, showing that the incompressible method plays a dominant role in determining the behaviour of the flow. This also means that it is not just the numerical interpolation scheme which is responsible for the overall accuracy. Furthermore, we show that the FSAC-PP method provides faster convergence and higher level of accuracy, making it a prime candidate for aerospace applications.


Introduction
Despite the advances made in the field of computational fluid dynamics over the past decades, predicting flow patterns around aerodynamic shapes remains a challenge for aerospace applications.The flow around a wing can have a transonic behaviour which, at high angles of attack, may be supplemented by flow separation, strong crossflow gradients as well as a hysteresis in the lift slope [1,2].Traub [3] highlighted further that at low Reynolds number flows, laminar separation bubbles exist which have an inherently unsteady behaviour.These separation bubbles may reattach to the wing or transition into a fully turbulent flow, depending on the pressure gradient.Panaras [4] reviewed recent published studies on Reynolds-averaged Navier-Stokes (RANS) and large-eddy simulation (LES) results and pointed out that the turbulent kinetic energy in regions of reversed flow is low and can be considered to be almost laminar.Thus, it is important to capture the transition between laminar and turbulent flows as flow separation and reattachment ultimately have a strong influence on predicting the lift slope, stall angle and hysteresis.Furthermore, since RANS models are commonly based on the linear Boussinesq assumption, results may provide only a moderate accuracy while models based on non-linear theories may be more efficient in predicting non-isotropic flow behaviour.
From the considerations given above, it is important to realize that the non-linear term in the Navier-Stokes equations is the most sensitive part and an accurate numerical description is mandatory to capture those flow features.The literature on numerical schemes is vast and this article does not intend to give a comprehensive overview-interested readers are referred to reference [5,6]-however a brief overview is given.Initially the second-order central scheme was used and later augmented by artificial dissipation.This was necessary to account for the lost dissipation on the small scales due to an under-resolved computational grid, essentially mimicking the same approach RANS turbulence model take in order to provide turbulent, or physical, dissipation.The central scheme lacks transportiveness and so is not suitable to predict flows containing discontinuities.A range of high-resolution schemes were developed to address this issue, of which the monotonic upwind scheme for conservation laws (MUSCL), total variation diminishing (TVD), essential non-oscillatory (ENO) and weighted essential non-oscillatory (WENO) schemes are the most widely used schemes nowadays.These schemes all have the same goal; to capture the correct behaviour of the non-linear term.However, their closure is based on numerical considerations.In order to enhance these scheme and to provide highly-accurate physical values for the time integration procedure, two possible routes are available.The first one is commonly employed in the compressible community where an approximate Riemann solver (RS) is used.In-fact, by including the local Riemann problem, the values obtained through the numerical scheme are modified based on the local eigenstructure of the system of equations, thus presenting a hybrid approach to couple numerical and physical fluxes.The second approach is to take a characteristics-based (CB) scheme which shares many similarities of the local Riemann problem, although being rather different in the details.Here, the characteristics of the flow are found at each computational node where the primitive variables are updated along the compatibility equations, which are unique to each node.The CB scheme found widespread use in the 1960s and onwards for the calculation of compressible flows [7].
With the introduction of the Artificial Compressibility (AC) method of Chorin [8], which was the first hyperbolic method available in the framework on incompressible flows, it became possible to apply both the local Riemann problem as well as CB schemes to incompressible flows.Drikakis et al. [9] were among the first to introduce the CB scheme in conjunction with the AC method in the finite volume method, which was followed by a development in the finite element framework by Zienkiewicz and Codina [10] and Zienkiewicz et al. [11].The CB scheme of Drikakis [9] was developed for one-dimensional flows but can be extended to higher dimensions by applying the same equations to any coordinate direction in space.Thus, the approach is termed the single-directional characteristics-based, or SCB, scheme.The scheme was later revised by Neofytou [12] and tested by Su et al. [13], concluding that, although the revised scheme being more mathematical rigorous, there was little difference in the results using the original and revised SCB scheme.
A vastly different approach was taken by Razavi et al. [14] who proposed to use a multi-directional characteristics-based, or MCB, scheme.It was extended to handle curve-linear coordinate systems [15] as well as unstructured grids [16].Improved boundary conditions in combination with the MCB scheme was proposed by Hashemi and Zamzamian [17] for the far-field and Zamzamian and Hashemi [18] for solid boundary conditions.Fathollahi and Zamzamian [19] investigated the influence of the number of compatibility equations used at each node.They showed that increasing the number of compatibility equation did not affect the accuracy of the solution while increasing the convergence rate.The scheme has been tested for heat transfer and turbulent flows by Razavi and Adibi [20] and Razavi and Hanifi [21], respectively.In the MCB scheme, instead of tracking characteristic lines through each computational node as done in the SCB scheme, characteristic surfaces are multiplied into the governing equations and their solution is found by means of the system's eigenvalues and geometrical considerations.This creates a characteristic cone along which the compatibility equations are emanating.This allows the scheme to capture any anisotropic properties of the flow and therefore lends itself to capture non-linear flow features.
In the present work, our aim is to show the applicability of the MCB scheme to predict non-linear flow behaviours at low reynolds numbers.As a comparison, we also employ the SCB scheme and apply Rusanov's RS which provides a hybrid multi-directional, Godunov-type treatment for the non-linear term of the Navier-Stokes equations.We apply the scheme to a channel with a sudden expansion.This geometry was found to be particularly suitable to predict the onset of bifurcation and subsequently the reattachment point of the flow.Both phenomena are underlying features that are important to calculate turbulent flows and thus highly relevant to the discussion in the context of aerospace applications.
We use the AC method in order to apply the CB scheme and RS, which require a hyperbolic system of equations.Furthermore, we also employ a recent method developed by Könözsy [22], later introduced by Könözsy and Drikakis [23], which unifies Chorin's AC method with Chorin's and Temam's Fractional-Step, Pressure projection (FS-PP) method [24,25] which is termed the Fractional-Step, Artificial Compressibility with Pressure-Projection, or FSAC-PP method (see Section 2.2).It has been tested for classical benchmark cases for incompressible flows [22], multi-species and variable density flows in micro-channels [26], trapping and positioning of cryogenic propellants [27], forced separated flows over a backward facing step [28] and the vortex pairing problem [29].The FSAC-PP method was originally introduced using the SCB scheme and extended by Teschner et al. [30] to the MCB scheme.Smith et al. [31] removed the CB approach and used different Riemann solvers instead to obtain the inviscid fluxes.The FSAC-PP has shown to have superior convergence properties over the AC and FS-PP method, especially in low Reynolds number flows, while showing a similar level of accuracy.
The remaining structure of this article is as follows.In Section 2, we give numerical details on the AC and FSAC-PP method.Section 3 summarises the SCB and MCB scheme used in this study while the computational setup is presented in Section 4. We present the results for the sudden expansion geometry in Section 5 and highlight final remarks in Section 6.

Methodology
In this section, we briefly review the concepts of the AC and FSAC-PP method, which are employed throughout this study.Details can be found in [8] about the AC method and in [22,23] about the FSAC-PP approach.

The Artificial Compressibility Method
The absence of a thermodynamic functional relationship between the pressure and density has led to the notion of an Artificial Compressibility approach.Here, the density is replaced by the pressure in the time derivative of the continuity equation.Since it could be difficult to give a clear relationship between the two state variables, a numerical constant is introduced, β, which is a user-defined convergence parameter [8].Thus, the continuity equation loses any physical meaning, however, it provides a mechanism to predict the pressure in the momentum equation.Once a steady state is obtained, the time derivative vanishes to zero and so a divergence-free velocity field is obtained.The governing equations with the extensions of the AC method thus become where ρ and ν are the density and viscosity, respectively.The equations are iterated in pseudo-time due to the non-physical time derivative of pressure.In order to advance the system of Equations ( 1) and (2) in real time, a dual-time stepping procedure needs to be employed for which a real time derivative is added to the momentum equation.The time-step has to satisfy the following condition as

The Unified Fractional-Step, Artificial Compressibility and Pressure Projection Method
The FSAC-PP method unifies both the AC and FS-PP method into a single framework.In the first Fractional-Step, the continuity equation of the AC method is used along with the momentum equation of the FS-PP method, in which the pressure gradient is dropped.The governing equations of the numerical method can be written in a semi-discretised form by From Equations ( 4) and ( 5), we obtain an intermediate velocity field denoted by U * which is used in the second Fractional-Step of the numerical procedure to recover the pressure field as The Helmholtz-Hodge decomposition requires the velocity field at time level (n + 1) to become divergence-free.Thus, by taking the divergence of Equation ( 6), the following Poisson equation for the pressure is obtained by With an updated pressure value, we can recover the velocity at the next time level from Equation (6) as Equations ( 5)-( 8) are consistent with the FS-PP method of Chorin [24] and Temam [25] while Equation (4) provides an initial pressure field for the Poisson equation through the perturbed continuity equation of the AC method.Thus, the hyperbolic and elliptic features of the AC and FS-PP method are coupled through the Fractional-Step procedure which allows for a characteristics-based treatment of the convective fluxes while the pressure is stabilised through the elliptic Poisson equation.

Characteristic-Based Schemes for Incompressible Flows
While the FSAC-PP method was originally introduced using the SCB scheme [9], Teschner et al. [30] showed the extension to a multi-directional characteristics-based (MCB) scheme, which is based on the derivation given by Razavi et al. [14] for the AC method.

A Single-Directional Closure
The SCB scheme splits the governing equations into one-dimensional equations and assumes that the SCB scheme is equally applicable for any direction.The primitive variables are recovered in the context of the hyperbolic-type AC method as where The values of x and ỹ are set to x = 1, ỹ = 0 when considering the x-direction and x = 0, ỹ = 1 for the y-direction, respectively.The eigenvalues are found as The primitive variables with indices j = 0, 1, 2 are found through Godunov's RS where the intercell flux values U L i+1/2 and U R i+1/2 are found through a third-order interpolation procedure, see Section 4. In the FSAC-PP method, the characteristic pressure is set to zero as it is dropped in the momentum equation.From Equations ( 9) and ( 10) it can be seen, however, that the pressure provided by the continuity equation enters the characteristic velocity field which enforces a stronger coupling between velocity and pressure.More details on the SCB scheme can be found in [9,22,23].

A Multi-Directional Closure
In contrast to the SCB scheme, the multi-directional approach applies a characteristic surface to the governing equations.For a two-dimensional flow, the governing equations of the AC method are Multiplying Equations ( 19)-( 21) by d f and introducing the shorthand notation The determinant of the coefficient matrix corresponding to Equations ( 22)-( 24) yields to solution in the form of Ψ 1 = 0 and Ψ 2 = β/(ρn τ ), where n τ = f τ is the component of the unit normal vector and use has been made of the fact that the normal direction corresponds to the derivative of the characteristic surface in the form of n = ∇ f .In-fact, it can be shown that Ψ 1,2 corresponds to two independent characteristic surfaces [7,[32][33][34] which are called the stream-(Ψ 1 ) and wave-(Ψ 2 ) surface in accordance with the literature on the compressible MCB scheme.Thus, inserting Ψ 1,2 into Equations ( 23) and ( 24) yields two sets of compatibility equations along with Equation (22).Rusanov [7] showed that a mix of compatibility equations on Ψ 1 and Ψ 2 are necessarily needed to obtain the characteristic values on the cell interfaces.For the incompressible version of the MCB scheme, however, it is sufficient to only use the compatibility equations corresponding to Ψ 2 .This can be verified by computing the rank of the coefficient matrix of Equations ( 22)-( 24) for Ψ 2 which will result in a non-rank deficient matrix.Thus, three compatibility equations are found over four wave directions as The characteristic pressure is an average over all four wave directions and ( U − U j ) indicates that the velocity components, which are aligned with the wave directions, are to be used.The normal vector is where j = 1, 2, 3, 4 indices correspond to the intersection of the characteristic surfaces with the time level (n), which is illustrated in Figure 1, and the primitive variables at those locations (wave directions) are found through interpolation.Razavi et al. [14] introduced two versions of the MCB scheme, a first-and second-order scheme, where primitive variables at j = 1, 2, 3, 4 are set to the reconstructed intercell variables at i + 1 2 for the first-order scheme and are interpolated along the characteristic surfaces in the second-order scheme.Since the SCB scheme classifies as a first-order CB scheme under this definition, the first-order MCB scheme is used in the remaining investigation to compare both methods.The SCB scheme was introduced using Godunov's RS, Equation (18), which provided the transportiveness property of the SCB scheme.The MCB scheme, as introduced by Razavi et al. [14], did not feature any RS treatment.Thus, in order to ensure the transportiveness of the MCB scheme, a numerical interpolation scheme needs to be selected which provides this property.Alternatively, the MCB scheme can be supplemented with a RS to remove the necessity to provide transportiveness through the numerical interpolation scheme and make it a property of the MCB scheme itself.Another benefit of RS-based approaches is that different RS provide different level of numerical dissipation which may become important at high Reynolds number flows.This combination allows to use a relative cheap but high-order interpolation scheme for increased accuracy at reduced computational cost, where we make the properties traditionally linked to the numerical scheme now a property of the CB framework.Furthermore, the inclusion of the RS promotes a multi-directional Godunov-type characteristic framework for incompressible flows which has so far only been tested and introduced by Teschner et al. [30] for a simple channel geometry.Therefore, we make use of Rusanov's approximate RS [35] in the form of where F( U) i+1/2 is the inviscid flux at the cell interface from which the primitive variables are found directly through flux differentiation.We follow the approach of Davis [36] to obtain the signal velocity with Using this closure, we have extended the MCB scheme into the Godunov framework and will demonstrate that even a simple approximate RS has favourable influence on the accuracy of the overall procedure.

Computational Setup and Numerical Schemes
The intercell values of the inviscid fluxes are obtained through a third-order interpolation polynomial as The viscous fluxes, on the other hand, use a second-order reconstruction scheme which, for Cartesian coordinates, defaults to a finite difference discretization in the form of The time integration is carried out by a third-order Runge-Kutta TVD scheme in which each stage is updated by the preceding one as where Here, we set α = 1 for the AC method and α = 0 for the FSAC-PP method.The numerical convergence parameter is set to β = 1.0, the CFL number for the explicit Runge-Kutta integration scheme is kept at CFL = 0.8 while the SOR algorithm applied to the Poisson solver, Equation (7), uses an under-relaxation factor of ω = 0.7 with a total of ten Poisson iteration per time-step.Könözsy [22] found that the Poisson equation can be over-relaxed to gain further computational speed-up and suggested to use ω = 1.7 instead.Unlike the FS-PP method, in which the Poisson solver is the most computational expensive part to obtain a correct pressure field, the FSAC-PP method benefits from an initial prediction of the pressure field through the continuity equation of the AC method which substantially reduces the number of iteration needed for the Poisson equation.It was found that of the order of ten iterations is sufficient to speed up the convergence rate while keeping the same order of accuracy as obtained with the FS-PP method.The residuals are calculated based on where R U is normalised based on the residual at the first iteration and compared against the convergence criterion of ε = 10 −12 to ensure that the residuals have reduced by 12 orders of magnitude based on the divergence free constrain.

Results and Discussion
To investigate the performance of the MCB scheme using the Rusanov RS against the classical, SCB scheme, we simulate the flow inside a suddenly expanding channel with an expansion ratio of 3:1, see Figure 2. We have performed a grid convergence study and calculated the grid convergence index (GCI) based on the reattachment length according to Roache [37].Due to the Cartesian nature of the solver, the mesh cannot be locally refined, especially near walls where a good resolution is required.This, however, creates a challenging test case for both the AC and FSAC-PP method as strong velocity gradients are present close to the wall which also amplifies any differences between the two methods.Table 1 shows the results obtained for the grid dependency study where we give the reattachment length at the upper and lower wall, as well as its corresponding GCI value and the time and iteration it took to compute the results.The Cartesian mesh elements were divided by a factor of two for subsequent mesh levels, ensuring a refinement ratio of four in all cases.We also provide the expected value based on the Richardson extrapolation which can also be found in [37].At a Reynolds number of Re = 30, for which the GCI study has been carried out, the expected reattachment length as given by Oliveira [38] is X r1 , X r2 = 3.080.Looking at the results, we can see that the FSAC-PP method converges towards that results for the finest level, while the AC method fails to predict the results correctly.The FSAC-PP method also shows a lower level of GCI which is of the order of 10% while the AC method produces GCI values which are close to double of that result.In terms of relative change, the reattachment length changes less than 5% from the medium to the finest mesh level for both AC and FSAC-PP method.Although the AC method overpredicts the reattachment, in this case using no CB scheme and no RS, we will see later that the results are slightly better with other schemes.Hence, focusing on the FSAC-PP method, it would be desired to use the finest mesh level for all subsequent simulations.The computational time, however, becomes prohibitive in the case of the AC method if the finest mesh level was to be used, where a single simulation would last for more than 32 h which would cause excessive simulation times as several simulations at different Reynolds numbers are required.Furthermore, for very low Reynolds numbers, the AC method requires even more computational time [23,27].Thus, we chose the medium mesh with 22,400 elements as the next closest mesh for which results in a reasonable amount of time can be obtained.
To ensure that the convergence criterion of 10 −12 is sufficient, we furthermore carried out a convergence study where we have checked the influence of the convergence parameter from 10 −6 to 10 −12 for Equation (40).To judge convergence, we check by how much the reattachment length changes on the upper and lower wall.The results are summarised in Table 2.We can see that the chosen convergence threshold at 10 −12 is a sufficient indicator for convergence.The results for the FSAC-PP method do not change from 10 −11 while the results for the AC method do only differ for the fourth significant figure.Thus, we haven chosen 10 −12 as a convergence threshold for all subsequent simulations.Table 1.Grid dependency study for the AC and FSAC-PP method at a Reynolds number of Re = 30.The expected value of the reattachment length is X r1 , X r2 = 3.080 as given by Oliveira [38]. Figure 3 shows contours of the u-velocity component while Figure 4 presents the velocity profile at x/S = 5, where S indicates the step height, see Figure 2. We show results here obtained with the AC and FSAC-PP method for two Reynolds number.The first Reynolds number is Re = 34.6,based on the average inlet velocity and a characteristic lengthscale of unity, which is below the critical Reynolds number of Re crit ≈ 54 [38,39] and where a symmetrical flow pattern prevails.The second Reynolds number of Re = 80 is above the critical one and exhibits a breaking in the symmetry.All velocity profiles are compared against experimental data of Fearn et al. [39].

Method
Figure 4a shows the AC method at a sub-critical Reynolds number.Not using any characteristics or RS overpredicts the flow slightly which is similar to the MCB scheme by itself.Both the standard non-linear treatment in conjunction with the Rusanov RS, as well as the hybrid Rusanov MCB scheme show an under-prediction of the peak velocity.A similar observation can be made for the SCB scheme.Since the SCB scheme has been proposed together with Godunov's RS, we can see that all RS-based approaches are more dissipative than a non RS-based approach.Figure 4b shows the AC method with the same schemes above the critical Reynolds number.Here, only the non-CB and MCB approach were capable of predicting the bifurcation to a non-symmetric flow behaviour.All RS-based approaches do not predict accurately the onset of bifurcation which may be explained by the inherent numerical dissipation to the Riemann solvers themselves.The picture is fundamentally different for the FSAC-PP method.At Re = 34.6, Figure 4c, the results are independent of the method, except for the SCB scheme, which under-predicts the peak velocity slightly more than the other schemes.In Figure 4d we see the bifurcated state and all schemes, including the RS-based approaches, do predict the velocity profile correctly.We see here that not just the numerical scheme, but also the incompressible method used to perform the calculations makes a difference.This highlights that not just the numerical scheme that is used makes a difference, but also the incompressible closure assumption of a given method can have a significant influence.It should be kept in mind that the AC method is purely hyperbolic while the FSAC-PP method has a mixed hyperbolic-elliptic behaviour and it is not surprising that both method perform differently.However, it is surprising that the FSAC-PP method does predict the bifurcation correctly, regardless of the numerical scheme, while the AC method is highly depending on the numerical scheme.This is against the common believe that only the numerical scheme influences the accuracy of the solution while the numerical methods may only differ in the number of iterations.We make a case here and show that different mathematical behaviours (hyperbolic and elliptic) do have a significant effect on the final result.Since the current study is investigating the non-linear term of the Navier-Stokes equations, this is an important finding and it shows that for studies where the non-linear term is of importance, such as stall prediction around wings, it is worth considering not just a suitable scheme for the convective term but also a suitable incompressible method.The mixed hyperbolic-elliptic behaviour of the FSAC-PP method means that any pressure disturbances in the flow will instantaneously propagate through the domain using the elliptic Poisson solver.The purely hyperbolic nature of the AC method, on the other hand, means that pressure disturbances travel at a fixed signal velocity, determined by the local eigenstructure of the system.The onset of bifurcation does critically depend on the pressure boundary value at the convex corners which, in turn, are influenced by the recirculation area upstream.The FSAC-PP method allows that any change in the recirculation area will influence the upstream convex corner points and vice versa.The AC method, on the other hand, lags behind due to the finite information propagation speed and thus faces difficulties in predicting the onset of bifurcation.In-fact, the same argumentation was used by Teschner et al. [40] where a novel incompressible method was devised based on a parabolic pressure transport equation.It was argued that the pressure treatment plays a crucial role for incompressible flows and that a parabolic behaviour of the pressure matches the physical expected behaviour more closely.Thus, the elliptic part of the FSAC-PP method plays a crucial part in the current discussion as the pressure disturbances are handled in a different way as in the AC method.The outcome has a rather drastic impact, as can be seen from the results, where the bifurcation is predicted by the FSAC-PP method but not by the AC method.At the same time, however, we need to point out that the bifurcation phenomena can be predicted using the AC method, see for example Drikakis [41].We do not claim that it is not capable of predicting the non-linear behaviour correctly, however, the AC method may be more prone to different modelling approaches, for example, it may become necessary to use stretched grids near the wall to capture the correct flow physics.The FSAC-PP method does not show the same modelling issues which can be attributed, in parts, to the Poisson solver.Not only does it provide an elliptic behaviour which advantages have been discussed above, but also it provides a stabilisation and smoothing procedure for the pressure.A stabilised pressure field, in turn, will provide a more realistic velocity field.Alternative approaches, using a fully hyperbolic compressible (Euler/Navier-Stokes) or incompressible (AC) solver are, for example, to use artificial dissipation, to stabilise the velocity field.The FSAC-PP method contains the stabilising feature by default which is another indicator as to why the result show not such a drastic difference as the AC method and furthermore predict the correct results regardless of the numerical scheme employed.In essence this is what we aim to demonstrate in this study, that we can take classical properties of numerical scheme, such as transportiveness and stability, and make them properties of the CB scheme as well as the numerical method.In this way, we provide a robust framework which is capable of treating non-linear and anisotropic flow behaviours correctly while modelling errors due to different reconstruction schemes are removed.
The discussion above is also given quantitatively in Tables 3 and 4. Focusing on the L 1 norm, we can see the order of magnitude is similar for both methods and Reynolds numbers, where the correct flow was predicted and is between 2% and 4%.The inability of the AC method under the current set up to predict the bifurcation results in rather large errors while the FSAC-PP method shows similar errors for the bifurcated and non-bifurcated state.We can also see that the application of the Rusanov RS produces larger errors for the FSAC-PP methods and also fails to predict the bifurcation for the AC method, while the non Rusanov RS-based approached predicted the bifurcation.This can be explained by the numerical dissipation inherent to the Rusanov RS.It may seem as a disadvantage for the present results, however, its advantages become only apparent at high Reynolds numbers, where the added dissipation acts as physical dissipation where a non-dissipative interpolation scheme is used.For the present study, where laminar flows are concerned, it is not expected to see a large difference in the result.Rather, it is worth mentioning that the added dissipation does not contaminated the results significantly for low Reynolds numbers in the case of the FSAC-PP method, while the effect is felt rather strongly in the case of the AC method, where the bifurcation is no predicted.Thus, compared to the AC method, we show that our hybrid Rusanov scheme in conjunction with the MCB scheme does predict the bifurcation correctly and that the gains are significant.The numerical results could be improved by using a less dissipative RS like the Harten, Lax, van Leer (HLL) or HLL-Contact (HLLC) RS, however, Smith et al. [31] highlighted that the signal speed prediction becomes problematic.Thus, we accept the numerical dissipation by the Rusanov RS but gain conservativeness through the RS which allows to use a relative simple numerical interpolation scheme.At the same time, once real aerospace applications are of interest, the numerical dissipation provided by the Rusanov RS is sufficient to produce stable results.This, by itself, could be regarded as an alternative Implicit Large Eddy Simulation (ILES) approach where we do not take the numerical dissipation from the interpolation scheme but rather from the Riemann solver directly.
Figure 5 shows the bifurcation diagram for the AC (Figure 5a) and FSAC-PP (Figure 5b) method.Here, we have defined DX as the difference between the upper and lower reattachment length and show results obtained in the range of 10 ≤ Re ≤ 100.We further show numerical results by Oliveira et al. [38] who used the SIMPLEC scheme to solve the incompressible Navier-Stokes equations and provided tabulated data on the reattachment length.As has been already observed, the bifurcation was only predicted by the non-CB and MCB scheme for the AC method while all schemes bifurcated for the FSAC-PP method.Figure 5a shows that the AC method predicted the onset of bifurcation prematurely.The FSAC-PP method in Figure 5b, on the other hand, shows that the bifurcation was predicted slightly after the critical Reynolds number.The FSAC-PP method, however, does follow the experimental and numerical data more closely, especially for higher Reynolds numbers.The reattachment length in the AC method becomes weakly invariant to the Reynolds numbers at Re = 80 while the numerical data suggest that the difference between the upper and lower reattachment point should continue to grow.The behaviour is the same for the MCB and non-CB scheme.The hyperbolic nature and thus finite information propagation speed in the AC method could explain this phenomena where the distance between the reattachment point and the upstream convex corner points becomes too large so that disturbances are not propagated fast enough for them to have an effect.In the FSAC-PP method, on the other hand, the disturbances are propagated instantaneously and so even at larger distances, or higher values of DX, the disturbances are still felt upstream.It is interesting to note here that all RS-based approaches in the FSAC-PP method do initially predict the onset of bifurcation later than the non-CB or MCB scheme.It could be argued that the inherent numerical dissipation to the RS delays bifurcation.At higher Reynolds number and sufficiently far away from the critical one, all schemes eventually converge towards the same solution.The results in Figure 5 is summarised in Tables 5 and 6 for the AC and FSAC-PP method, respectively.Here we give both the upper and lower reattachment point and compare against the tabulated data of Oliveira [38].We further give the number of iterations it took to get converged results up to our convergence criterion of ε.We can confirm from the data given that the non-CB and MCB scheme may indeed match the reference data more closely, especially at high Reynolds numbers.The difference are, however, minute in the case of the FSAC-PP method.
The computational cost is shown in Figure 6, which shows an interesting trend.The AC method is know to have slow convergence properties at low Reynolds numbers, see for example [22,23].We see the same trend in Figure 6a, where the computational time required grows exponentially as the Reynolds number approaches zero.The non-CB and MCB scheme require substantially more CPU time than the SCB and RS-based approaches, however, the latter did not predict the bifurcation at all.The reason here is that the flow will always develop as a symmetrical flow, even for Reynolds number above Re crit .After the residuals have dropped to about ε = 10 −10 , the numerical fluctuations become small enough so that the physical fluctuations can promote a different and more stable equilibrium position.At this stage, the reattachment point at the upper and lower wall start to interact with the upstream pressure at the convex corner points and slight physical fluctuations-which obtain their energy through the non-linear term-promote the change to a non-symmetrical flow pattern, which is also discussed by Oliveira [38].Therefore, results obtained for the bifurcated state may in general require more computational time.At the critical Reynolds number, both the AC and FSAC-PP method peak in terms of the CPU time.It is here that the physical fluctuations become just important enough for the flow to register the change to a non-symmetrical state.At higher reynolds numbers, their effects may be felt stronger and earlier during the calculation which may force the bifurcation to occur earlier and thus with increasing Reynolds numbers, the computational time decreases again.When only considering the non-CB and MCB scheme of the AC method, which eventually bifurcate, it can be seen that the same simulation using the FSAC-PP method requires 2-3 times less CPU time sufficiently far away from the critical Reynolds number.This is consistent with findings of previous works in which it was highlighted that the FSAC-PP method generally performs faster than the AC method while retaining a high level of accuracy [22,23,28,30].At Reynolds number close to zero, the computational savings are even more significant.This was one of the reasons the FSAC-PP method was developed in the first place, to overcome the stiff nature of the AC and Pressure Projection method for low Reynolds number flows, for example in microfluidic applications, see [23] for a detailed discussion.At the onset of bifurcation, however, the AC method seems to be slightly more cost efficient than the FSAC-PP method.This can be explained by the elliptic behaviour of the Poisson solver, where any pressure change is propagated through the domain instantaneously.These pressure waves require longer to be damped while the hyperbolic behaviour of the AC method induces less pressure waves and hence a converged solution is found quicker.

Conclusions
In this work, we presented results to predict the highly non-linear behaviour produced by the Navier-Stokes equations at low Reynolds numbers inside a channel with a sudden expansion.In particular, we investigated the sub-and super-critical range of Reynolds numbers where the flow bifurcates from a symmetric to a non-symmetric state.We investigated the performance of a non-characteristic-based (CB) scheme, single-and multi-directional characteristics-based scheme (SCB/MCB), as well as the Rusanov Riemann solver (RS) and combinations of these schemes and tested all approaches with the Artificial Compressibility (AC) and Fractional-Step, Artificial Compressibility with Pressure-Projection (FSAC-PP) method.
We found that only the non-CB and the MCB scheme were capable of predicting the bifurcation using the AC method where the RS-based approaches showed too much numerical dissipation to correctly predict the flow patterns.A significant difference between the SCB and MCB scheme could be observed in the AC method, where we showed that the multi-directional nature of the MCB scheme is required to predict the bifurcation at all.The added Rusanov RS showed too much numerical dissipation for the current approach to predict the bifurcation and similar results were obtained using the hybrid Rusanov and MCB scheme.In the FSAC-PP method, however, all schemes correctly predicted the symmetry breaking and overall showed better agreement with reference data in terms of the reattachment length.The incompressible flow method itself could overcome the inherent numerical dissipation of the Rusanov RS, as well as the SCB scheme which uses Godunov's RS by default, showing that the mathematical classification of the method's partial differential equations play a dominant role.The fully hyperbolic behaviour of the AC method was not always capable of capturing the bifurcation phenomena correctly while the mixed hyperbolic-elliptic equations of the FSAC-PP method ensured always a physically correct solution.This comes at a slightly increased computational cost near the bifurcation, however, sufficiently far away from the critical Reynolds number, the FSAC-PP method required 2-3 times less CPU resources compared to the AC method.
Aerospace applications are presented with flow features that are highly non-linear, as discussed in Section 1.The bifurcation shares many similarities to the onset to turbulence and predicting both phenomena correctly is mandatory to gain highly accurate computations in which stall, strong crossflow gradients and lift slope hysteresis are of interest.We have demonstrated here that the success of predicting non-linear flow features is not just scheme dependent, but the incompressible method used for the computation also plays a dominant role.The MCB scheme by itself showed that its multi-directional nature was capable of predicting non-linear phenomena correctly, regardless of the incompressible flow method used, while the SCB scheme was only capable of predicting correctly the flow phenomena when the FSAC-PP method was used.Although only laminar flows were investigated in this study, at higher Reynolds numbers, as is characteristic of aerospace applications, the flow becomes turbulent.In these flow regimes, the numerical dissipation produced by different schemes and RS becomes important.The Rusanov RS provides a sufficient amount of inherent numerical dissipation to tackle high Reynolds number turbulent flows.In future work, we will investigate the proposed framework for high Reynolds number flows, however, the framework present in the current study is directly applicable to such flows for aerospace applications.

Figure 1 .
Figure 1.Intersection of the characteristic surfaces with time level (n).

Figure 2 .
Figure 2. Geometry of the channel with a sudden expansion.Boundary conditions are given for the inflow, outflow and solid (wall) boundary condition.

Figure 3 .
Figure 3. Contour plots of the axial velocity component using a non-CB treatment.

Figure 4 .
Figure 4. Velocity profiles at x/S = 5 for the AC (top) and FSAC-PP (bottom) method using different combinations of CB schemes and the Rusanov Riemann solver at Re = 34.6 (left) and Re = 80 (right).

Figure 5 .
Figure 5. Bifurcation diagram for the AC (a) and FSAC-PP (b) method using different combinations of CB schemes and the Rusanov Riemann solver.

Figure 6 .
Figure 6.Total number of iterations required for the AC (a) and FSAC-PP (b) method for different Reynolnds numbers.

Table 2 .
Convergence study for different level of convergence at a Reynolds number of Re = 30 based on the reattachment length X r1 and X r2 at the upper and lower wall.

Table 3 .
L 0 and L 1 error norm of the velocity profiles for the AC method using different combinations of CB schemes and the Rusanov Riemann solver at Re = 34.6 and Re = 80.

Table 4 .
L 0 and L 1 error norm of the velocity profiles for the FSAC-PP method using different combinations of CB schemes and the Rusanov Riemann solver at Re = 34.6 and Re = 80.

Table 5 .
Prediction of the reattachment length X r1 and X r2 and number of iterations required for convergence for different Reynolds numbers using the AC method for various combinations of CB schemes and the Rusanov Riemann solver.

Table 6 .
Prediction of the reattachment length X r1 and X r2 and number of iterations required for convergence for different Reynolds numbers using the FSAC-PP method for various combinations of CB schemes and the Rusanov Riemann solver.