Y Guevel

T Allain

G Girault

J M Cadou

Numerical bifurcation analysis for 3-dimensional sudden expansion fluid dynamic problem

Keywords: asymptotic-numerical method, branch switching, continuation, high-performance computing, Navier-Stokes, nonlinear solvers, numerical bifurcation analysis, sudden expansion, 3-dimensional flow

This paper deals with bifurcation analysis methods based on the asymptotic-numerical method. It is used to investigate 3-dimensional (3D) instabilities in a sudden expansion. To do so, highperformance computing is implemented in ELMER, ie, an open-source multiphysical software. In this work, velocity-pressure mixed vectors are used with asymptotic-numerical method-based methods, remarks are made for the branch-switching method in the case of symmetry-breaking bifurcation, and new 3D instability results are presented for the sudden expansion ratio, ie, E = 3. Critical Reynolds numbers for primary bifurcations are studied with the evolution of a geometric parame-ter. New values are computed, which reveal new trends that complete a previous work. Several kinds of bifurcation are depicted and tracked with the evolution of the spanwise aspect ratio. One of these relies on a fully 3D effect as it breaks both spanwise and top-bottom symmetries. This bifurcation is found for smaller aspect ratios than expected. Furthermore, a critical Reynolds number is found for the aspect ratio, ie, A i = 1, which was not previously reported. Finally, pri-mary and secondary bifurcations are efficiently detected and all post-bifurcated branches are followed. This makes it possible to plot a complete bifurcation diagram for this 3D case.

INTRODUCTION

Numerical bifurcation analysis of flow problems governed by the discretized Navier-Stokes equations with a large number of degrees of freedom is still a challenging task. A review of the existing numerical bifurcation methods in fluid dynamics was proposed in the work of Dijkstra et al. [START_REF] Dijkstra | Numerical bifurcation methods and their application to fluid dynamics: analysis beyond simulation[END_REF] Toolbox exists for ordinary differential equations, such as MatCont [START_REF] Dhooge | MATCONT: a MATLAB package for numerical bifurcation analysis of ODES[END_REF] and AUTO-07p, [START_REF] Doedel | AUTO 94: Software for Continuation and Bifurcation Problems in Ordinary Differential Equations[END_REF] hence with a small number of degrees of freedom. [START_REF] Kłosiewicz | Using pseudo-arclength continuation to trace the resonances of the Schrödinger equation[END_REF] A large-scale toolbox for systems of discretized partial differential equations is the library of continuation algorithms LOCA [START_REF] Salinger | LOCA 1.0 Library of continuation algorithms: Theory and implementation manual[END_REF] as part of the Trilinos framework. [START_REF]The Trilinos Project[END_REF] An alternative to the well-known incremental-iterative methods [START_REF] Winters | Bifurcation and stability: a computational approach[END_REF][START_REF] Doedel | Numerical analysis and control of bifurcation problems (I): Bifurcation in finite dimensions[END_REF][START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF] is the asymptotic-numerical method (ANM). [START_REF] Cochelin | Méthode Asymptotique Numérique. Collection méthodes numériques[END_REF] This ANM is the association of a perturbation technique and a spatial discretization method, in our case the finite element method (FEM). It has been successfully applied in nonlinear solid mechanics [START_REF] Bekhoucha | Branch switching at Hopf bifurcation analysis via asymptotic numerical method: application to nonlinear free vibrations of rotating beams[END_REF][START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF] or in fluid mechanics. [START_REF] Guevel | Numerical comparisons of high-order nonlinear solvers for the transient Navier-Stokes equations based on homotopy and perturbation techniques[END_REF][START_REF] Medale | High performance computations of steady-state bifurcations in 3D incompressible fluid flows by asymptotic numerical method[END_REF][START_REF] Guevel | Parametric analysis of steady bifurcations in 2D incompressible viscous flow with high order algorithm[END_REF] The main advantage of the ANM is its ability to determine analytical nonlinear solutions [START_REF] Cochelin | A path-following technique via an asymptotic-numerical method[END_REF][START_REF] Cadou | ANM for stationary Navier-Stokes equations and with Petrov-Galerkin formulation[END_REF] with computational CPU time lower than with that of the classical incremental-iterative method. On these nonlinear solutions, some bifurcation indicators have been introduced. [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF][START_REF] Cadou | A numerical method for the computation of bifurcation points in fluid mechanics[END_REF][START_REF] Guevel | Automatic detection and branch switching methods for steady bifurcation in fluid mechanics[END_REF][START_REF] Jawadi | Asymptotic numerical method for steady flow of power-law fluids[END_REF] High-performance computation based on the ANM for flow problems was proposed in the work of Medale and Cochelin [START_REF] Medale | A parallel computer implementation of the asymptotic numerical method to study thermal convection instabilities[END_REF] as an in-house code.

Parametric studies were proposed for the 3-dimensional (3D) sudden expansion in other works. [START_REF] Durst | Low Reynolds number flow over a plane symmetric sudden expansion[END_REF][START_REF] Fearn | Nonlinear flow phenomena in a symmetric sudden expansion[END_REF][START_REF] Baloch | On two-and three-dimensional expansion flows[END_REF][START_REF] Battaglia | Bifurcation characteristics of flows in rectangular sudden expansion channels[END_REF] Numerical bifurcation analysis using the aspect ratio E = 3 is available in other works, [START_REF] Schreck | Numerical study of bifurcation in three-dimensional sudden channel expansions[END_REF][START_REF] Chiang | Side wall effects on the structure of laminar flow over a plane-symmetric sudden expansion[END_REF][START_REF] Chiang | Spanwise bifurcation in plane-symmetric sudden-expansion flows[END_REF][START_REF] Tsui | Side-wall effects on the bifurcation of the flow through a sudden expansion[END_REF] but only the first primary bifurcation was investigated. More recently, a new powerful feature of the ANM, [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF] which makes it possible to enhance the detection and continuation in the vicinity of a simple bifurcation, has been applied to study additional primary bifurcation for several spanwise aspect ratio. [START_REF] Medale | High performance computations of steady-state bifurcations in 3D incompressible fluid flows by asymptotic numerical method[END_REF] In this work, we propose to complete previous studies by tracking each kind of primary bifurcation when the aspect ratio evolves. Moreover, secondary bifurcations are investigated to plot a complete bifurcation diagram. Furthermore, new information are given for the spanwise symmetry-breaking bifurcation described in the works of Medale and Cochelin [START_REF] Medale | High performance computations of steady-state bifurcations in 3D incompressible fluid flows by asymptotic numerical method[END_REF] and Chiang et al. [START_REF] Chiang | Side wall effects on the structure of laminar flow over a plane-symmetric sudden expansion[END_REF] Regarding the ANM-based method, remarks are made for the branch-switching technique in the case of symmetry-breaking bifurcation. Lastly, previous ANM-based methods applied to fluid dynamics problems [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF][START_REF] Medale | High performance computations of steady-state bifurcations in 3D incompressible fluid flows by asymptotic numerical method[END_REF][START_REF] Guevel | Parametric analysis of steady bifurcations in 2D incompressible viscous flow with high order algorithm[END_REF][START_REF] Cadou | ANM for stationary Navier-Stokes equations and with Petrov-Galerkin formulation[END_REF][START_REF] Cadou | A numerical method for the computation of bifurcation points in fluid mechanics[END_REF][START_REF] Guevel | Automatic detection and branch switching methods for steady bifurcation in fluid mechanics[END_REF][START_REF] Jawadi | Asymptotic numerical method for steady flow of power-law fluids[END_REF][START_REF] Medale | A parallel computer implementation of the asymptotic numerical method to study thermal convection instabilities[END_REF][START_REF] Cadou | Numerical tools for the stability analysis of 2D flows: application to the two-and four-sided lid-driven cavity[END_REF] made no use of pressure information because static condensation was used. In this work, full velocity-pressure mixed vectors are used in the ANM-based methods, such that flows and specific vectors are evaluated with pressure information.

This paper is organized as follows. The numerical methods are described in the first part. Implementation of the ANM in ELMER is briefly mentioned in the second part. In the third part, bifurcation analysis is proposed for 3D flow in a sudden expansion. Within this part, information of interest are compared to the literature, the evolution of the first primary bifurcation is investigated, and bifurcation modes are depicted with velocity and pressure components. Lastly, a bifurcation diagram is proposed to study secondary bifurcations.

NUMERICAL METHODS

Governing equations

The steady Navier-Stokes equations for Newtonian and incompressible fluid in a domain (Ω) are as follows:

-2μ div D (u) + grad p + ρu • grad u = 0 in (Ω)

(1)

div u = 0 in (Ω) (2)
u = λu d on (Γ), (3)
with u and p being the velocity and pressure unknowns, respectively, and λ being a control parameter. A generic velocity profile u d is prescribed on the boundary (Γ). D (u) is the symmetric part of the velocity gradient tensor. Both density ρ and kinematic viscosity μ are constant in the fluid domain. The Reynolds number is defined as Re = λ ρu max d h μ , with h being a characteristic length of the model.

The weak form of the Navier-Stokes equations can be defined according to

L(U) + Q(U, U) = λF in (Ω), (4)
where the linear operator L (•) contains the diffusion, the pressure gradient, and the continuity terms. The nonlinear operator Q (•, •) contains the convection terms. Dirichlet boundary condition Equation (3) produces the right-hand-side (RHS) term, which is equivalent to a fictitious force. [START_REF] Cochelin | A path-following technique via an asymptotic-numerical method[END_REF][START_REF] Cadou | ANM for stationary Navier-Stokes equations and with Petrov-Galerkin formulation[END_REF]

Continuation with perturbation method

Pseudo-arc-length path-following technique [START_REF] Cochelin | A path-following technique via an asymptotic-numerical method[END_REF] based on the ANM is used to perform the continuation of the steady flow solutions. This has been successfully implemented for 2-dimensional (2D) flows in the case of steady incompressible Navier-Stokes equations in other works [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF][START_REF] Cadou | ANM for stationary Navier-Stokes equations and with Petrov-Galerkin formulation[END_REF][START_REF] Cadou | A numerical method for the computation of bifurcation points in fluid mechanics[END_REF][START_REF] Guevel | Automatic detection and branch switching methods for steady bifurcation in fluid mechanics[END_REF] and for non-Newtonian fluids in the work of Jawadi et al. [START_REF] Jawadi | Asymptotic numerical method for steady flow of power-law fluids[END_REF] Recently, 3D flow study has been proposed in the work of Medale and Cochelin. [START_REF] Medale | High performance computations of steady-state bifurcations in 3D incompressible fluid flows by asymptotic numerical method[END_REF] The ANM relies on a linearization technique where the unknowns of the nonlinear problem are sought as power series as follows:

X(a) = X 0 + N ∑ i=1 a i X i , ( 5
)
where X 0 = {U 0 λ 0 } ⊤ is a known regular solution with the mixed vector notation U 0 = {u 0 p 0 } ⊤ , a ∈ R is a perturbation parameter, and N is the truncation order of this polynomial representation. Although several definitions exist, the perturbation path parameter is chosen as follows:

a = ⟨u -u 0 ⟩u 1 + (λ -λ 0 )λ 1 , (6)
where ⟨•, •⟩ is the classic Euclidean dot product. Introducing Equation ( 5) in the nonlinear system Equation ( 1) and equating like powers of a, a set of linear algebraic systems is obtained. The same tangent operator is used for each of those linear systems and is defined as follows:

L 0 t (•) = L (•) + Q (U 0 , •) + Q (•, U 0 ). (7)
The first order reads

L 0 t (U 1 ) = λ 1 F (8a) ⟨u 1 , u 1 ⟩ + λ 2 1 = 1, (8b)
and the recurrence relation

for k ∈ [2, N] is L 0 t (U k ) = λ k F - k-1 ∑ i=1 Q(U i , U k-i ) (9a) ⟨u k , u 1 ⟩ + λ k λ 1 = 0. ( 9b
)
Once the series {U i , λ i } i=1,N are computed using Equation ( 8)-( 9), small part of the steady flow solution branch is approximated continuously with the polynomial representation Equation (5).

The range of validity of this approximation is defined as follows:

a maxpoly = ( η ||U 1 || ||U N || ) 1∕(N-1) , ( 10
)
with η being a user-defined tolerance. Other definitions may be used. [START_REF] Cochelin | A path-following technique via an asymptotic-numerical method[END_REF] Finally, the ANM requires 2 user parameters, ie, the first one being the truncation order N and the second one is the tolerance parameter η. The latter governs the accuracy of the computed nonlinear predictions. Lastly, new regular solution X 0 = X(a max ) is evaluated with Equation ( 5) and new continuation step might be performed. We recall that the approximation of the flow solution is continuous in a path-following continuation step. This makes it possible to evaluate steady flow solutions at an exact Reynolds number.

Usually, the ANM nonlinear prediction does not need correction at the end of a step. Nevertheless, for some values of the chosen parameter η, the accuracy of the solution obtained with the ANM cannot be satisfactory. It means that the computed residual is greater than a given accuracy. In such a case, very efficient and cheaper correctors can be used at the end of the continuation step to improve the quality of the ANM solution (see, for example, the work of Cadou and Potier-Ferry 31 ).

Steady bifurcation detection

Once steady nonlinear solution branch has been computed, it is worth determining the critical Reynolds numbers for which a steady bifurcation appears. These singularities often characterize the loss of uniqueness of the flow solution. In this study, only simple symmetry-breaking bifurcations of steady flow are considered. Nevertheless, these techniques might be applied for other kinds of bifurcation analysis. [START_REF] Bekhoucha | Branch switching at Hopf bifurcation analysis via asymptotic numerical method: application to nonlinear free vibrations of rotating beams[END_REF][START_REF] Medale | A parallel computer implementation of the asymptotic numerical method to study thermal convection instabilities[END_REF] Van Dyke [START_REF] Van Dyke | Analysis and improvement of pertubation series[END_REF] showed that power series analysis may reveal some features about singularity in its range of validity. Thus, the ANM series analysis has been proposed for steady bifurcation detection in the work of Cochelin and Medale. [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF] It has been applied to Newtonian incompressible fluid flows in the works of Cochelin and Medale [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF][START_REF] Medale | High performance computations of steady-state bifurcations in 3D incompressible fluid flows by asymptotic numerical method[END_REF] and for non-Newtonian fluids in the work of Jawadi et al. [START_REF] Jawadi | Asymptotic numerical method for steady flow of power-law fluids[END_REF] It is demonstrated in the work of Cochelin and Medale [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF] that the pseudo-arc-length continuation is perturbed in the neighborhood of a singularity. Geometric progression appears in the ANM series in the vicinity of a steady bifurcation point. The power series {X i } computed in the singularity vicinity using Equation ( 8)-( 9) is, in fact, a linear combination of a geometric progression series with a nearly flawless ones as follows:

{X i } = { Xi } flawless + {X i } geometric , for i = 1, N. ( 11
)
Determining a steady bifurcation consists in finding a geometric progression into the polynomial representation Equation (5). This is numerically detected via collinearity condition and relative error test applied to the last terms of the series as follows:

N-2 ∑ p=N-3 ( |α p | 1∕(N-p) |α n-1 | -1 ) 2 < ϵ gp1 , with α p = ⟨ X p , X N ⟩ ⟨X N , X N ⟩ (12)
N-1

∑ p=N-3 ( ‖ ‖ X ⟂ p ‖ ‖ ∕ ‖ ‖ X p ‖ ‖ ) < ϵ gp2 , with X ⟂ p = X p -α p X N . ( 13
)
The parameters ϵ gp1 = 10 -6 and ϵ gp2 = 10 -3 gave satisfactory results for steady bifurcation detection in this study. It should be noted that a large truncation order, around N = 30 in this study, is required to accurately detect the geometric progression. Once the common ratio 1∕α c = 1∕α N-1 and the associated scale factor α N c X N are determined, a nearly flawless enhanced series is computed as follows:

Xi = X i -α N-i c X N , for i = 1, N. ( 14
)
In the work of Cochelin and Medale, [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF] it is demonstrated that α c is exactly the arc-length distance to the critical point. Thus, the flow solution at this singular point is computed using the flawless enhanced series as follows:

X c = X(a = α c ) = X 0 + N-1 ∑ i=1 α i c Xi . ( 15
)
This constitutes an easy method to detect accurately critical solution because it is done while performing continuation steps with no additional computational effort.

Branch switching at simple bifurcation point

At simple steady bifurcation point, families of flow solution branches might be found using the classical ANM path-following technique by varying the ANM user parameters N and η. [START_REF] Cadou | ANM for stationary Navier-Stokes equations and with Petrov-Galerkin formulation[END_REF] It requires a lot of computations with no possible guess of which parameters lead to a specific branch. Guevel et al [START_REF] Guevel | Automatic detection and branch switching methods for steady bifurcation in fluid mechanics[END_REF] proposed a branch-switching method for the Navier-Stokes equations resolution. Then, approximations of nonlinear post-bifurcated branches emanating from the singular solution at a simple bifurcation point are computed using an adapted ANM continuation. This branch-switching technique relies on the traditional bifurcation theory. [START_REF] Doedel | Numerical analysis and control of bifurcation problems (I): Bifurcation in finite dimensions[END_REF][START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF][START_REF] Keller | Numerical solution of bifurcation and nonlinear eigenvalue problems[END_REF][START_REF] Brezzi | Finite dimensional approximation of nonlinear problems part III: simple bifurcation points[END_REF] It has been successfully coupled with the ANM in other works [START_REF] Bekhoucha | Branch switching at Hopf bifurcation analysis via asymptotic numerical method: application to nonlinear free vibrations of rotating beams[END_REF][START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF][START_REF] Vannucci | An asymptotic-numerical method to compute bifurcating branches[END_REF][START_REF] Boutyour | Bifurcation points and bifurcated branches by an asymptotic numerical method and Padé approximants[END_REF] in the solid mechanics framework. Because FEM operators in fluid mechanics problems are unsymmetric, an adaptation is proposed in the works of Medale and Cochelin [START_REF] Medale | High performance computations of steady-state bifurcations in 3D incompressible fluid flows by asymptotic numerical method[END_REF] and Guevel et al. [START_REF] Guevel | Automatic detection and branch switching methods for steady bifurcation in fluid mechanics[END_REF] First, tangents of the 2 intersecting branches at a simple bifurcation point are determined. By choosing one of those tangents, it makes it possible to switch from a branch to another. Those nonlinear post-bifurcated branches are computed with a modified ANM path-following technique as presented in the work of Guevel et al. [START_REF] Guevel | Automatic detection and branch switching methods for steady bifurcation in fluid mechanics[END_REF] Hence, the nonlinear post-bifurcated branches are sought as a power series representation as follows:

X bi (a) = X c + N ∑ j=1 a j X bi j , for i = 1, 2. ( 16
)
The starting point is the singular solution X c computed as in Equation (15). Injecting the polynomial representation Equation ( 16) in the nonlinear system Equation (1), a set of linear algebraic systems is obtained equating like powers of a. The tangent operator denoted as L c t is assembled as in Equation (7) with the solution at the bifurcation point. As it is singular, a special care is required to compute the post-bifurcated polynomial representation using this operator. Thus, tangents at the bifurcation point are the first term of this representation.

At a simple bifurcation point, the tangent operator has the following properties:

Ker ( L c t ) = {𝚽} , 𝚽 ∈ Ω, ⟨𝚽, 𝚽⟩ = 1 (17) Ker ( L c⊤ t ) = {𝚿} , 𝚿 ∈ Ω, ⟨𝚿, 𝚽⟩ = 1, (18)
with Φ being the so-called bifurcation mode and 𝚿 being the associated left bifurcation mode. Moreover, there exists a unique particular solution vector W such that

L c t (W) = F, W ∈ Ω, < W, 𝚽 >= 0. ( 19
)
Different strategies are available to determine those vectors. [START_REF] Keller | Numerical solution of bifurcation and nonlinear eigenvalue problems[END_REF] A way to compute those vectors is described in the following. Additionally, all those specific vectors are computed in this work as full velocity-pressure mixed vectors.

Tangents determination

At first order of the path parameter a, the tangents

X bi 1 = { U bi 1 λ bi 1 } i=1,2
are solutions of the following linear system:

L c t ( U bi 1 ) = λ 1 F (20a) ⟨ U bi 1 , U bi 1 ⟩ + ( λ bi 1 ) 2 = 1. ( 20b
)
The singularity at the bifurcation point is treated classically with the help of the Lyapunov-Schmidt reduction. [START_REF] Guevel | Automatic detection and branch switching methods for steady bifurcation in fluid mechanics[END_REF][START_REF] Keller | Numerical solution of bifurcation and nonlinear eigenvalue problems[END_REF][START_REF] Vanderbauwhede | Lyapunov-Schmidt method for dynamical systems[END_REF][START_REF] Golubitsky | Singularities and Groups in Bifurcation Theory[END_REF] The 2 tangents are written as follows:

U bi 1 = λ bi 1 W + η bi 1 𝚽, (21)
with λ bi 1 , η bi 1 ∈ R being 2 scalars to be determined. To do so, solve the following linear system defined at the second order:

L c t ( U bi 2 ) = λ 2 F -Q ( U bi 1 , U bi 1 ) (22a) ⟨ U bi 2 , U bi 1 ⟩ + λ bi 2 λ bi 1 = 0. (22b)
Projecting Equation (22a) on the left bifurcation mode 𝚿, Equation (18) leads to the well-known Algebraic Bifurcation Equation (ABE) as follows:

⟨ 𝚿, Q ( U bi 1 , U bi 1 )⟩ = 0. ( 23
)
Using the Lyapunov-Schmidt reduction Equation ( 21), the following quadratic ABE is obtained:

a b ( λ bi 1 ) 2 + b b λ bi 1 η bi 1 + c b ( η bi 1 ) 2 = 0. ( 24
)
This quadratic equation may be easily solved using Equation (20b), with the following coefficients:

a b = ⟨𝚿, Q (W, W)⟩ (25a) b b = ⟨ 𝚿, Q (Φ, W) ⟩ (25b) c b = ⟨𝚿, Q (Φ, Φ)⟩ . ( 25c
)
In Equation (25b), the operator

Q(a, b) is defined according to the expression Q(a, b) = Q(a, b) + Q(b, a).
The general case with restriction to η bi 1 ≠ 0 and a b ≠ 0 gives

t bi = λ bi 1 η bi 1 = -b b + (-1) i √ b 2 b -4a b c b 2a b . ( 26
)
Using the norm condition Equation (20b) and the property ⟨W, 𝚽⟩ = 0, it is obtained:

η bi 1 = 1 √ (t bi ) 2 ⟨W, W⟩ + ⟨𝚽, 𝚽⟩ + (t bi ) 2 , ( 27
)
where λ bi 1 is automatically deduced using Equation (26).

Remark on symmetry-breaking bifurcation case

Pitchfork bifurcation occurs at simple bifurcation point if flow symmetry is broken for one of the post-bifurcated family branch. Therefore, the coefficients defined in Equation (25a) are checked numerically. If a b = c b = 0 and b b ≠ 0, it implies that pitchfork bifurcation occurs at the detected critical point. [START_REF] Werner | The computation of symmetry-breaking bifurcation points[END_REF] One is easily convinced that the condition c b = 0 is equivalent to the fact that the bifurcation mode 𝚽 is the solution of the ABE proposed in Equation (23) and that it is one of the sought tangents. [START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF][START_REF] Werner | The computation of symmetry-breaking bifurcation points[END_REF] Moreover, a b = 0 implies that the other tangent is then collinear to the particular solution W.

Finally, in the case of pitchfork bifurcation, the 2 tangents are as follows:

X ba 1 = { U ba 1 λ ba 1 } = { 𝚽 0 } ( 28
)
X bs 1 = { U bs 1 λ bs 1 } = { λ bs 1 W λ bs 1 } (29)
with

λ bs 1 = 1 √ ⟨W, W⟩ + 1 . ( 30
)

Specific vectors computation

In order to evaluate the coefficients in Equation ( 25) or directly the pitchfork bifurcation tangents as in Equation ( 29)-( 28), the specific vectors 𝚽, W, and 𝚿 are required.

In other works, [START_REF] Cadou | A numerical method for the computation of bifurcation points in fluid mechanics[END_REF][START_REF] Guevel | Automatic detection and branch switching methods for steady bifurcation in fluid mechanics[END_REF][START_REF] Vannucci | An asymptotic-numerical method to compute bifurcating branches[END_REF] the bifurcation mode Φ is computed using a classic ANM linear stability analysis technique known as the "bifurcation indicator". It should be noted that this requires some computational efforts. Then, only W and Ψ remain to be computed using the definitions Equation ( 18)-( 19) by the mean of augmented systems presented in the following. This is done with no assumption of the symmetry-breaking simple bifurcation case.

The recent development proposed in the works of Cochelin and Medale [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF][START_REF] Medale | High performance computations of steady-state bifurcations in 3D incompressible fluid flows by asymptotic numerical method[END_REF] does not use the same Lyapunov-Schmidt reduction. Nevertheless, in the case of the symmetry-breaking bifurcations, vectors 𝚽 and W can be computed. The first tangent at the bifurcation point is obtained by differentiating the actual flawless enhanced power series Equation ( 16) at the critical point as follows:

X b1 1 = dX(a) da | | | |a=α c = N-1 ∑ i=1 iα i-1 c Xi . ( 31
)
A Gram-Schmidt technique is used to extract the orthogonal vector from the last term of the series as follows:

X b1 ⟂ 1 = X N - ⟨ X N , X b1 1 ⟩ ⟨ X b1 1 , X b1 1 ⟩ X b1 1 . ( 32
)
In the case of symmetric branch, the following equality holds: 𝚽 = X b1 ⟂ 1 . Contrariwise, if the continuation is performed on a post-bifurcated solution, from an antisymmetric state, then the bifurcation mode is the tangent computed as 𝚽 = X b1 1 . Moreover, the particular solution W might be computed using Equation (31) or Equation ( 32) because of the orthogonality relation. Instead, it is used as an augmented system, which makes it possible to compute W and is reused in the following to circumvent the singularity of the tangent operator. This does not require more computational effort. The vector W, which is the solution of the linear equation in Equation (19), is orthogonal to the bifurcation mode 𝚽 so that the augmented system to be solved is written as follows:

[ L c t 𝚽 𝚽 ⊤ 0 ] { W κ } = { F 0 } , ( 33
)
where κ acts as a Lagrangian multiplier, which would be as small as possible.

In the case of the Navier-Stokes equations, the left mode remains to be computed using either an augmented system, [START_REF] Guevel | Automatic detection and branch switching methods for steady bifurcation in fluid mechanics[END_REF] or with iterative techniques such as bordering technique as presented in the work of Cochelin and Medale, [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF] or with an inverse power method. [START_REF] Medale | High performance computations of steady-state bifurcations in 3D incompressible fluid flows by asymptotic numerical method[END_REF] Then, using the same consideration for the left bifurcation mode 𝚿 (see Equation ( 18)), the augmented system is written as follows:

[ L c ⊤ t 𝚽 𝚽 ⊤ 0 ] { 𝚿 κ } = { 0 1 } . ( 34
)
It is noticed that this augmented matrix is the transpose of the precedent one in Equation (33). Specific options might be used in linear solvers to reuse the already factorized matrix Equation (33). In this way, only one factorization is required to compute the specific vectors and for every linear system required for the post-bifurcated branches as described in the following.

Post-bifurcated branch continuation

Once tangents are determined as solutions of Equation ( 20), high-order terms of the power series Equation ( 16) are determined for each of the post-bifurcated branches as solution of

L c t ( U bi k ) = λ bi k F - k-1 ∑ j=1 Q ( U bi j , U bi k-j ) (35a) ⟨ U bi k , U bi 1 ⟩ + λ bi k λ bi 1 = 0. ( 35b
)
The following Lyapunov-Schmidt reduction [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF][START_REF] Guevel | Automatic detection and branch switching methods for steady bifurcation in fluid mechanics[END_REF] is used:

U bi k = λ bi k W + η bi k 𝚽 + V bi k . ( 36
)
Injecting Equation (36) in Equation ( 35), the new vector

V bi k is solution of L c t ( V bi k ) = - k-1 ∑ j=1 Q ( U bi j , U bi k-j ) (37a) ⟨ V bi k , 𝚽 ⟩ = 0. ( 37b
)
This system of equations is exactly the linear augmented system described in Equation ( 33), with different RHS vector. It makes it possible to compute V bi k in an efficient way because only a new RHS needs to be assembled. The coefficients λ bi k and η bi k are determined with the same procedure used for the tangents. Injecting the Lyapunov-Schmidt reduction in the order k + 1 system Equation (35a) projected on the left bifurcation mode 𝚿 and using the condition Equation (35b), the following system is obtained:

[ ⟨ 𝚿, Q ( W, U bi 1 ) ⟩ ⟨ 𝚿, Q ( 𝚽, U bi 1 ) ⟩ ⟨ W, U bi 1 ⟩ + λ bi 1 ⟨ 𝚽, U bi 1 ⟩ ] { λ bi k η bi k } = { -g k - ⟨ V bi k , U bi 1 ⟩ } ( 38
)
with the coefficient definition

g 2 = ⟨ 𝚿, Q ( V bi 2 , U bi 1 )⟩ (39)
and, for k > 2,

g k = ⟨ 𝚿, Q ( V bi k , U bi 1 )⟩ + k-1 ∑ j=2 ⟨ 𝚿, Q ( U bi j , U bi k-j+1 )⟩ . ( 40
)
When a complete power-series is obtained, the range of validity is evaluated. Then, for a chosen tangent, the associated post-bifurcated branches emanating from the bifurcation points are approximated using Equation ( 16), with a max or -a max in 2 directions. This procedure is repeated at wish for each of the computed tangents. Each new starting point on the post-bifurcated branches is either saved for a restart procedure or used directly for new continuation steps.

It should be noted that, in the case of simple steady bifurcation, 4 branches are approximated using 2 tangents and new regular solutions are saved for next restart.

Remark on symmetry-breaking bifurcation case for high-order terms

For the symmetric post-bifurcated branch, the tangent X bs 1 is orthogonal to the bifurcation mode, which leads to the following simplification:

λ bs k = - ⟨ V bs k , W ⟩ ⟨W, W⟩ + 1 (41a) η bs k = - g bs k b b . (41b)
Then, using the tangent of the antisymmetric branch X ba 1 , which is exactly the bifurcation mode, the sought coefficients are

λ ba k = - g ba k b b (42a) η ba k = 0, (42b)
leading to the following expression of the high-order terms:

U ba k = λ ba k W + V ba k . ( 43
)
In this particular case, each high-order term of the series is orthogonal to the bifurcation mode according to the pseudo-arc-length condition Equation (35b).

NUMERICAL IMPLEMENTATION

All the previous numerical methods were implemented in the open-source software ELMER, [START_REF]ELMER -Finite Element Solver for Multiphysical Problems[END_REF] which is a generic and multiphysics finite element code. Several tests were run with the UMFPACK linear direct solver, which is the default direct solver. For problems requiring memory space greater than 4 GB, the use of an alternative direct solver was mandatory.

In this study, the solver MUMPS [START_REF] Amestoy | A fully asynchronous multifrontal solver using distributed dynamic scheduling[END_REF] was chosen. Accordingly, it was compiled with the OpenBLAS linear algebra library [START_REF] Xianyi | Model-driven level 3 BLAS performance optimization on Loongson 3A processor[END_REF] to get multithreaded operations. Finally, ELMER was compiled to work with MUMPS. All the ANM-based methods and algorithms were implemented in ELMER as new modules and user-defined solvers. * Additional features are available with the solver MUMPS that makes it possible to compute both Equation (33) and Equation (34) using only one factorization of the augmented tangent operator.

NUMERICAL RESULTS

Bifurcation analysis is performed for the case of symmetry-breaking steady bifurcation. The critical Reynolds numbers are investigated for the 3D sudden expansion.

The ANM continuation of steady flow solution is performed. Then, the geometric progression criteria Equation ( 12) and Equation ( 13) are checked during the continuation. Possibly, steady bifurcations are detected. The critical Reynolds numbers and bifurcation modes are reported. Results obtained with the ANM methods in ELMER are compared with the bifurcation information of the literature.

Computations presented in this paper were performed on a DELL Precision T7500 using 2 Intel(R) Xeon(R) X5677 @ 3.47 GHz (16 threads) with 96 GB of memory or on a Dell PowerEdge R930 using 2 Intel(R) Xeon(R) CPU E7-8860 v4 @ 2.20 GHz (72 threads) with 1 TB of memory.

This section is organized as follows. Firstly, the 3D sudden expansion case is described. Secondly, the first primary bifurcation is described for several geometric ratios. Then, the other primary bifurcations are investigated. Finally, secondary bifurcations are presented for one geometric ratio.

Test case description

The 3D sudden expansion test case is well documented in other works. [START_REF] Durst | Low Reynolds number flow over a plane symmetric sudden expansion[END_REF][START_REF] Fearn | Nonlinear flow phenomena in a symmetric sudden expansion[END_REF][START_REF] Baloch | On two-and three-dimensional expansion flows[END_REF][START_REF] Battaglia | Bifurcation characteristics of flows in rectangular sudden expansion channels[END_REF] Bifurcation analysis using the aspect ratio E = 3 was proposed numerically in other works [START_REF] Schreck | Numerical study of bifurcation in three-dimensional sudden channel expansions[END_REF][START_REF] Chiang | Side wall effects on the structure of laminar flow over a plane-symmetric sudden expansion[END_REF][START_REF] Chiang | Spanwise bifurcation in plane-symmetric sudden-expansion flows[END_REF][START_REF] Tsui | Side-wall effects on the bifurcation of the flow through a sudden expansion[END_REF] and, more recently, in the work of Medale and Cochelin 14 using ANM-based methods. The fluid domain and the boundary conditions are reported in Figure 1. The length L is chosen long enough to capture the steady bifurcation phenomenon. Mesh features with several geometric ratios A i = W∕h and length L are presented in Table 1.

Linear interpolation H8/8C stabilized finite elements are used. [START_REF] Zienkiewicz | The Finite Element Method, Fluid Dynamics[END_REF] *https://github.com/YGuevel/ElmerMAN An established velocity profile is prescribed on the inlet. Many formulations can be found in the literature. The best-fit profile compared to a gravitational driven flow has been used in this study.

Firstly, Shah and London 44(eq. 332) and Titarev and Shakhov 45(eq. 3.3) proposed similar analytical formulation based on infinite sums of hyperbolic terms. Secondly, in the works of Baloch et al, 24(eq. 18) Theofilis et al, [START_REF] Theofilis | Viscous linear stability analysis of rectangular duct and cavity flows[END_REF](eq. 2.3) and Nicolas et al, 47(eq. 7) another formulation is obtained. Moreover, in the work of Spiga and Morino, 48(eq. 7) only trigonometric functions are used.

Finally, the analytic formulation used in this work is based on the work of Tanyeri et al. 49(S1) It uses infinite sum of hyperbolic terms in z for the spanwise direction multiplied with a classic parabolic profile for y, the height. It has been decided to modulate this latter parabolic profile as proposed in the work of Shah and London 44(eqs. 334-337) using a power exponent n. The inlet velocity profile that best fits a gravitational driven flow using ELMER is

u dx (𝑦, z) = C ( 1 - ( 2𝑦 h ) 2 ) n ∞ ∑ k=1,3,.. 1 k 3 ⎛ ⎜ ⎜ ⎜ ⎝ 1 - cosh ( kπz h ) cosh ( kπW 2h ) ⎞ ⎟ ⎟ ⎟ ⎠ , ( 44
)
where C is a normalization coefficient. Several tests were performed such that the values k ∈ [1, 15] and n = 0.95 have been used in this work.

Primary bifurcation

The ANM continuation is performed to follow the steady solution branches. The ANM parameters are N = 30 and η = 10 -16 . Steady bifurcation detection is performed with ϵ gp1 = 10 -6 and ϵ gp2 = 10 -3 . Using these parameters, the ANM method requires from 10 to 14 steps to detect the first bifurcation. Other parameters, for example, η = 10 -8 , allow to detect accurately the first primary critical values with less steps (around 5 to 7 step). However, it leads to degraded quality of solutions. Some computational cost information are given. The CPU factorization time presented in Table 1 is highly decreased because a multithreaded library has been used. Usually, in ANM, the tangent operator factorization represents the main part of the computational time. In this work, when using fine meshes with long bodies, it requires 3 "multithreaded LU" CPU time to perform one ANM step. The geometric progression detection is based on vector manipulations, which is negligible. Then, some specific vectors need to be computed. Due to MUMPS V5, it is possible to factorize only one augmented system to get all the specific vectors and to compute the post-bifurcated branches. Finally, the branch-switching technique allows to compute 4 parts of the post-bifurcated branches using only 1.66 ANM step. This represents half of a working day for standard meshes to a day or two for meshes with almost 10 7 degrees of freedom using 36 threads for the factorization step.

Readers may get further details of the ANM computational cost from other works. [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF][START_REF] Guevel | Numerical comparisons of high-order nonlinear solvers for the transient Navier-Stokes equations based on homotopy and perturbation techniques[END_REF][START_REF] Medale | High performance computations of steady-state bifurcations in 3D incompressible fluid flows by asymptotic numerical method[END_REF][START_REF] Guevel | Parametric analysis of steady bifurcations in 2D incompressible viscous flow with high order algorithm[END_REF][START_REF] Cadou | ANM for stationary Navier-Stokes equations and with Petrov-Galerkin formulation[END_REF][START_REF] Cadou | A numerical method for the computation of bifurcation points in fluid mechanics[END_REF][START_REF] Guevel | Automatic detection and branch switching methods for steady bifurcation in fluid mechanics[END_REF][START_REF] Jawadi | Asymptotic numerical method for steady flow of power-law fluids[END_REF] FIGURE 2 Bifurcation diagram for A i = 8. FB1 and FB2 denote fundamental branches, PB1 denotes the first primary bifurcation point, and BBr1n and BBr1p denote the post-critical unsymmetric nonlinear branches. Blue: fully symmetric solutions. Red: one ANM step using the singular initial solution. Orange: bifurcated branches with a broken symmetry but still not unsymmetrical

Bifurcation diagram

A velocity probe is used at coordinates (100,0,0) at the intersection of the vertical and horizontal symmetry planes. A 2D bifurcation diagram may be plotted using either u y or u z versus the Reynolds number. This choice depends on which symmetry breaking is investigated. A bifurcation diagram is proposed for the first primary bifurcation in Figure 2. From the fundamental branch FB1, a pitchfork bifurcation is detected. Then, the singular solution U c and the 3 specific vectors required by the branch-switching technique are computed. The red curve depicted one ANM step starting with a singular solution. A new regular point is chosen to follow an other branch.

Flow topology descriptions

Illustration of the 3D flow is depicted in Figure 3. Two Reynolds numbers are considered, ie, the critical value Re c = 103 and a higher value Re = 174 for which 3 solutions exist (see dashed line in Figure 2). The streamlines of the flow solutions allow us to identify recirculation zones and vortex lines.

The first primary bifurcation breaks the horizontal plane, and the vertical symmetry plane being still present.

In order to validate the flow topology, features are determined using skin friction lines. [START_REF] Délery | Skin Friction Lines Pattern and Critical Points[END_REF] It is depicted in skin friction lines on the floor, roof, and lateral wall in Figure 4 for the flow solution on the BBr1n branch. The kind of singularities and their positions are in good agreement with the pattern depicted in the work of Chiang et al. [START_REF] Chiang | Side wall effects on the structure of laminar flow over a plane-symmetric sudden expansion[END_REF]

Evolution of Re c vs A i

Critical Reynolds numbers for the first primary bifurcation are reported in Table 2 for several geometric aspect ratios A i . Values from the literature are also reported. Moreover, a transition from a symmetrical state for A i = 1.5 at Re = 314 to an unsymmetrical one for A i = 2 at Re = 297 was presented in the work of Tsui and Wang. [START_REF] Tsui | Side-wall effects on the bifurcation of the flow through a sudden expansion[END_REF] In the work of Chiang et al, 27 a transition was noted from a symmetrical state for A i = 3.5 at Re = 145 to unsymmetrical one for A i = 3.75 at Re = 143. Those transitions are reported respectively in orange and red in Figure 5. Additional information is given regarding critical Reynolds numbers found in the literature. In the work of Tsui and Wang, [START_REF] Tsui | Side-wall effects on the bifurcation of the flow through a sudden expansion[END_REF] inlet flow is uniform and mean velocity is used for the Reynolds number definition, so that coefficients are applied to match our definition. In Table 3, some ratios u max ∕ [START_REF] Doedel | AUTO 94: Software for Continuation and Bifurcation Problems in Ordinary Differential Equations[END_REF] 2 u mean are given. Coefficient for aspect ratios A i = 1.5 and A i = 2 are numerically evaluated using simulation of a flow in a rectangular duct. For higher aspect ratios, values are extracted from table 4 in the work of Chiang et al. [START_REF] Chiang | Side wall effects on the structure of laminar flow over a plane-symmetric sudden expansion[END_REF]

FIGURE 5

Critical Reynolds number associated to the first primary steady bifurcation for several A i using the asymptotic-numerical method (ANM)-based method. In black circles, the value from literature. In red and orange, the transition from symmetric flow (S) to antisymmetric flow (A) found in the literature Firstly, it can be observed that our results are in very good agreement with the literature. Secondly, between standard mesh and finer and/or longer mesh (see Table 1), critical values are very similar. This indicates that standard mesh with L = 30 is sufficient to accurately compute critical values for the first primary bifurcation.

To the best knowledge of the authors, no result was found in the literature for the ratio A i = 1. In the work of Tsui and Wang, [START_REF] Tsui | Side-wall effects on the bifurcation of the flow through a sudden expansion[END_REF] no primary bifurcation was found for Reynolds number lower than Re = 300. Using standard mesh, we were not able to compute accurately fundamental solutions for Reynolds number higher than Re = 450. However, with a finer and longer mesh, the first primary bifurcation was detected at Re cPB1 = 734.7.

These values are plotted in Figure 5. The trend for high values of aspect ratio is to converge slowly toward the critical Reynolds number Re c = 80.4 obtained for the 2D flow. [START_REF] Guevel | Parametric analysis of steady bifurcations in 2D incompressible viscous flow with high order algorithm[END_REF] This is in perfect agreement with critical values for higher ratios computed in the work of Medale and Cochelin. [START_REF] Medale | High performance computations of steady-state bifurcations in 3D incompressible fluid flows by asymptotic numerical method[END_REF] Finally, it is observed that the first primary bifurcation occurs for higher Reynolds number as the aspect ratio is decreasing.

Bifurcation vectors description

A spatial representation of the bifurcation modes 𝚽 and 𝚿 and of the particular solution W is depicted in Figure 6. In the same way, the pressure and the velocity parts of these 2 bifurcation modes are depicted in Figures 7 and8, respectively.

The velocity part of the bifurcation mode 𝚽 is in agreement with the work of Medale and Cochelin. [START_REF] Medale | High performance computations of steady-state bifurcations in 3D incompressible fluid flows by asymptotic numerical method[END_REF] Moreover, the vertical middle cut plane for the pressure part is similar to the 2D case. [START_REF] Battaglia | Bifurcation of low Reynolds number flows in symmetric channels[END_REF] It is observed that the left bifurcation mode 𝚿 drives the inlet flow. This behavior is in agreement with the work of Battaglia et al. 51 4 Critical Reynolds numbers and bifurcation type of second and third primary bifurcation for several geometric ratios. "Std" stands for standard mesh, "F/L" for finer or longer mesh, "NS" for not searched.

Re cPB2

Re

Evolution of additional primary bifurcations

Beyond the first primary bifurcation, fundamental branches are computed for various geometric ratios. Then, additional primary bifurcations are found. For each primary bifurcation, it was observed that bifurcation mode topologies were very closed to each others. Moreover, it is reported in Tables 4 and5 the values of the critical Reynolds number for the primary bifurcation with respect to the aspect ratio A i . Accordingly, it was proposed to track the evolution of the critical Reynolds numbers for each primary bifurcation with respect to the aspect ratio A i (see Figure 10). Whatever the rank of the primary bifurcation, it was observed that the greater the aspect ratio is, the lower the critical Reynolds numbers are. The different kinds of bifurcation modes for the primary bifurcation noted as ▴, □, and ⧫ are depicted in Figures 11 and12 for the case of A i = 8.

TABLE 5

Critical Reynolds numbers and bifurcation type for the fourth and fifth primary bifurcations versus the aspect ratio A i . "Std" stands for standard mesh, "F/L" for finer or longer mesh, "NS" for not searched, and "NF" for not found In Table 4, the critical Reynolds numbers are similar for standard and finer or longer meshes. Some bifurcations were not searched using standard meshes. The recirculation zones were not correctly captured with L = 30 because it was necessary to refine the mesh for higher Reynolds number. However, the recirculation zones were correctly computed using finer or longer meshes.

Re cPB4

For the ratio A i = 10, critical values issued from the work of Medale and Cochelin [START_REF] Medale | High performance computations of steady-state bifurcations in 3D incompressible fluid flows by asymptotic numerical method[END_REF] are Re cPB2 = 137.60, Re cPB3 = 210.86, and Re cPB4 = 280.05. Firstly, only two of these 3 additional bifurcations are detected. Even with a Newton corrector using a residual norm criterion prescribed at 10 -6 and a finer mesh, we were not able to detect the bifurcation at Re c = 210.86. However, the critical values and kinds of bifurcation are in agreement with the literature.

The primary bifurcation, noted as □, breaks both planar symmetries. The spanwise symmetry is broken. The flow is in the symmetry group Z 2 (y) ⊕ Z 2 (z) as described in the work of Chiang et al. [START_REF] Chiang | Spanwise bifurcation in plane-symmetric sudden-expansion flows[END_REF] More precisely, this is an X-axis symmetry. In Figure 9, a visual comparison between flows on post-bifurcated branches BBr1 (post •) and BBr2 (post □) is proposed. Firstly, in the case of the spanwise bifurcation, planes z = -25 and z = 25 are antisymmetric. Moreover, slice at middle plane shows a symmetric pattern. Whereas, the topology of the BBr1 flow for the 3 slices is in the same top-bottom direction. Secondly, contours make it possible to see differences between the BBr1 branch flow, which breaks the top-bottom symmetry and is in the planar symmetry group Z 2 (z), and the BBr2 branch flow, which breaks both plane of symmetry and is in the Z 2 ( y) ⊕ Z 2 (z) group. The latter flow is clearly symmetric according to the X-axis. This spanwise symmetry breaking was not found below A i = 30 in the work of Chiang et al. [START_REF] Chiang | Spanwise bifurcation in plane-symmetric sudden-expansion flows[END_REF] In the work of Medale and Cochelin, 14 this bifurcation was still present for ratios A i ≥ 10. In the present study, this spanwise bifurcation persists from A i = 7 to A i = 10. It is not found for smaller ratios. Another interesting behavior with A i = 7, is the permutation in the order of primary bifurcation. Using standard mesh, the spanwise bifurcation is found at Re c = 313.8 whereas the other bifurcation is found at Re c = 313.3. These 2 critical values are really closed to each other, but modes and flows on post-bifurcated branches are clearly different (see next subsection). Using a finer mesh, the spanwise bifurcation is found early at Re c = 308.3 against Re c = 314.8 for the other one.

In Table 5, the fourth and fifth primary bifurcations are reported. The critical values found differ about 2% depending of the meshes used. ⧫ bifurcation type is found at Re c = 568.8 for A i = 8 and Re c = 518.7 for A i = 10. However, for A i = 9, the detection of the ⧫ bifurcation type failed with the criterion proposed in the work of Cochelin and Medale [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF] and retained for this study. Nevertheless a pole appears in the rational representation of the series (as used in the work of Guevel et al [START_REF] Guevel | Automatic detection and branch switching methods for steady bifurcation in fluid mechanics[END_REF] ) indicating a bifurcation for Re c = 545. In this case, we were not able to determine the bifurcation type. Finer spatial discretization or finite element of higher degree might be needed to ensure a perfect detection.

Evolution of critical Reynolds number against the spanwise aspect ratio is plotted in Figure 10. The trend of the ▴ and ⧫ bifurcation type is similar to that of the first primary bifurcation described in Section 4.2.3. As the body is narrow, which corresponds to small geometric ratio A i , there is a stabilization effect on the bifurcation phenomenon. In the case of the spanwise breaking primary bifurcation □, this stabilization effect occurs for higher ratio. The trend is not the same for this □ primary bifurcation type. It might be explained by the fact that these kinds of bifurcation are comparable to a 2D bifurcation, whereas the spanwise symmetry-breaking primary bifurcation is a real 3D effect.

Primary bifurcation modes

The bifurcation types □, ▴, and ⧫ are now described for the ratio A i = 8 (see Figures 11 and12). Velocity components and pressure of bifurcation modes are proposed for the case A i = 8 using L = 55h.

The bifurcation modes 𝚽 PB3 (▴) and 𝚽 PB4 (⧫) are both similar in their topology with the first primary bifurcation 𝚽 PB1 seen in Figure 7. The bifurcation mode 𝚽 PB2 (□), which breaks the spanwise symmetry, has a different structure. The latter is in perfect agreement with the mode described in the work of Medale and Cochelin. [START_REF] Medale | High performance computations of steady-state bifurcations in 3D incompressible fluid flows by asymptotic numerical method[END_REF]

Case with secondary bifurcation for A i = 8

In previous subsections, only primary bifurcations were discussed. Now, results for the ratio A i = 8 with secondary bifurcation (SB) are presented. A secondary bifurcation appears on a nonfundamental post-bifurcated branch where a symmetry remains.

FIGURE 13

Three-dimensional bifurcation diagram obtained using the asymptotic-numerical method for the ratio A i = 8 using a fine and long mesh. Two distinct 2D bifurcation diagrams are also plotted as projections on the considered plane [Colour figure can be viewed at wileyonlinelibrary.com]

Critical Reynolds numbers

In Table 6, the critical Reynolds number of primary and secondary bifurcations depending on the initial fundamental branch part as presented in Figure 2 is reported.

Bifurcation diagram

A 3D bifurcation diagram is presented in Figure 13, where solution branches rely on u y and u z at probe as in Section 4.2.1. The first part of this diagram has been presented in Figure 2. The complexity of this diagram is clearly seen. Thus, a synthetic diagram is presented in Figure 14 with information on the state of symmetry of the branches. Nevertheless, the ANM-based method makes it possible to monitor the branch solutions and to perform branch switching very easily. A remark is made for the second primary bifurcation, which breaks the spanwise symmetry, [START_REF] Medale | High performance computations of steady-state bifurcations in 3D incompressible fluid flows by asymptotic numerical method[END_REF][START_REF] Chiang | Spanwise bifurcation in plane-symmetric sudden-expansion flows[END_REF] but does not appear on this diagram. This is due to the choice of the probe on the X-axis. The flow is symmetric only with respect to the X-axis so that both velocity components u y and u z at probe stay null.

This diagram is proposed in a synthetic version as presented in the work of Golubitsky and Schaeffer. [START_REF] Golubitsky | Singularities and Groups in Bifurcation Theory[END_REF] It gathered both antisymmetric solutions (such as BBr1n and BBr1p) into one family branch. Moreover, the symmetry group remaining to be broken is recalled. For example, Z 2 (z) is the planar symmetry regarding the vertical plane. In this case, the spanwise symmetry is still present. The double plane or X-axis symmetry is classically noted as Z 2 ( 𝑦) ⊕ Z 2 (z).

Secondary bifurcation modes

The secondary bifurcation modes are now described. Velocity components and pressure are proposed for the case A i = 8 using L = 55h. Those bifurcations lead to either a fully nonsymmetric flow, or let the previous symmetry unchanged after the bifurcation point.

In the case of the first primary bifurcation, the nontrivial post-bifurcated branch still conserves the spanwise symmetry. Both 𝚽 PB1SB1 and 𝚽 PB1SB2 in Figure 15 are similar in topology.

In the case of the second primary bifurcation, the nontrivial post-bifurcated branch still has a symmetry to be broken, ie, Z 2 (𝑦) ⊕ Z 2 (z). Here, 𝚽 PB2SB1 and 𝚽 PB2SB2 in Figure 16 are not similar in topology. 𝚽 PB2SB1 looks like 𝚽 PB3 as a top-bottom symmetry-breaking bifurcation mode.

The last secondary bifurcation mode 𝚽 PB3SB1 in Figure 17 has a complex structure without well-defined features.

CONCLUSIONS

A numerical bifurcation analysis for the 3D sudden expansion has been proposed. Continuation of steady flow solutions, bifurcation detection, and branch-switching techniques are implemented using the ANM. Power series analysis makes it possible to accurately detect bifurcation points. Specific continuation technique is used in the case of the path following of postcritical branch of steady flow solutions. These techniques are implemented as new modules and user-defined solver in ELMER. Implementation of MUMPS as serial direct linear solver using the multithread features of OpenBlas allows us to perform 3D flow simulation up to 9.7 million degrees of freedom in this present study. Large-scale problems are now possibly studied, either only for steady flow continuation or for a detailed bifurcation analysis. Three-dimensional sudden expansion critical values and flow topologies are correctly reproduced. The critical Reynolds numbers and modes for bifurcation are in good agreement with the literature. The proposed simplified branch-switching method at a pitchfork bifurcation is efficient. The stabilized FEM used in ELMER does not perturb the bifurcation analysis. Furthermore, while previous ANM studies [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF][START_REF] Medale | High performance computations of steady-state bifurcations in 3D incompressible fluid flows by asymptotic numerical method[END_REF][START_REF] Guevel | Parametric analysis of steady bifurcations in 2D incompressible viscous flow with high order algorithm[END_REF][START_REF] Cadou | ANM for stationary Navier-Stokes equations and with Petrov-Galerkin formulation[END_REF][START_REF] Cadou | A numerical method for the computation of bifurcation points in fluid mechanics[END_REF][START_REF] Guevel | Automatic detection and branch switching methods for steady bifurcation in fluid mechanics[END_REF][START_REF] Jawadi | Asymptotic numerical method for steady flow of power-law fluids[END_REF][START_REF] Medale | A parallel computer implementation of the asymptotic numerical method to study thermal convection instabilities[END_REF][START_REF] Cadou | Numerical tools for the stability analysis of 2D flows: application to the two-and four-sided lid-driven cavity[END_REF] did not use pressure information because static condensation was used, pressure component is successfully used in this work within the ANM bifurcation analysis framework. As full velocity-pressure mixed vectors are used, flows and specific vectors are presented in this study with pressure information.

New results are obtained. Firstly, evolution of the critical Reynolds number for several primary bifurcations has been determined for aspect ratio A i from 1 to 10. Those primary bifurcations are strongly stabilized with small spanwise aspect ratios. Different kinds of bifurcation were tracked and depicted. Secondly, we determined that the spanwise bifurcation described in the works of Medale and Cochelin [START_REF] Medale | High performance computations of steady-state bifurcations in 3D incompressible fluid flows by asymptotic numerical method[END_REF] and Chiang et al [START_REF] Chiang | Spanwise bifurcation in plane-symmetric sudden-expansion flows[END_REF] exists for ratios A i ≥ 7. Thirdly, a critical Reynolds number is found at Re cPB1 = 734 for the aspect ratio A i = 1. This was not found in the work of Tsui and Wang [START_REF] Tsui | Side-wall effects on the bifurcation of the flow through a sudden expansion[END_REF] for Re ≤ 300. To the best knowledge of the authors, this critical value was not previously reported in the literature. Finally, a bifurcation diagram and description of the primary and secondary bifurcation modes were proposed for the case A i = 8.

A generic tool is available for large-scale bifurcation analysis of the 3D Navier-Stokes equations. Eigensolver might be used in ELMER to determine the stability of flow solutions. Additional techniques based on the ANM may now be implemented, such as parametric analysis of steady bifurcation using homotopy, [START_REF] Guevel | Parametric analysis of steady bifurcations in 2D incompressible viscous flow with high order algorithm[END_REF] bifurcation analysis of non-Newtonian fluids, [START_REF] Jawadi | Asymptotic numerical method for steady flow of power-law fluids[END_REF] Hopf bifurcation detection, [START_REF] Cadou | A numerical method for the computation of bifurcation points in fluid mechanics[END_REF][START_REF] Marquet | Sensitivity analysis and passive control of cylinder flow[END_REF][START_REF] Brezillon | A numerical algorithm coupling a bifurcating indicator and a direct method for the computation of Hopf bifurcation points in fluid mechanics[END_REF][START_REF] Girault | An algorithm for the computation of multiple Hopf bifurcation points based on Padé approximants[END_REF][START_REF] Heyman | Computation of Hopf bifurcations coupling reduced order models and the asymptotic numerical method[END_REF] and transient nonlinear solvers. [START_REF] Guevel | Numerical comparisons of high-order nonlinear solvers for the transient Navier-Stokes equations based on homotopy and perturbation techniques[END_REF][START_REF] Razafindralandy | Time integration algorithm based on divergent series resummation, for ordinary and partial differential equations[END_REF] It is also challenging to study multiple bifurcations or to couple the reduced-order model techniques with the ANM-based methods for 3D flows.

Lastly, the series analysis detection method is based on the properties of singular tangent operators during a continuation procedure. Then, it is expected to be able to detect Hopf bifurcation using a specific continuation technique.

A.2 Comparisons for 2D and 3D sudden expansions

Combination of the iterative linear solver BiCGStab with precondtioner ILUt has been compared with the direct linear solver MUMPS. In the following, ANM truncation order has been set to N = 20. Only one ANM step is under study. Firstly, the same 2D geometry has been used with 3 mesh refinements. The CPU time and memory consumption are recalled in Table A3. The following 2 conclusions are made. MUMPS is faster. Moreover, the iterative process BiCGStab+ILUt uses approximately the same amount of memory than MUMPS.

Secondly, the flow in a 3D sudden expansion test case with the features E = 3, A i = 3, l = 3h, L = 30h has been studied with 3 mesh refinements. The H8/8C mixed finite element [START_REF] Zienkiewicz | The Finite Element Method, Fluid Dynamics[END_REF] has been used. An important remark is made regarding ILUt. It was not possible to perform all calculations with the combination BiCGStab+ILUt, whereas too much memory was required in ELMER. Thanks to the ELMER development team who advised us to use BiCGStab(l) with the parameter l = 4. The CPU time and memory usage are recalled in Table A4. In conclusion, MUMPS is faster than BiCGStab(l = 4) with no preconditioner, but more memory is required. Some remarks are made regarding the choice of a linear solver for the ANM. The choice of the correct iterative linear solver and the corresponding best preconditioner or other parameters might depend on the part of the branch that is evaluated. A more precise study has to be performed in the future to give more detailed explanations on those peculiar behaviors. Lastly, it appeared that the multithreaded MUMPS direct linear solver used in this study made it possible to largely decrease the computational time regarding the iterative solvers available in ELMER. On the other side, a large amount of memory was required.

FIGURE 1

 1 FIGURE 1 Sudden expansion. h = 10, H = 3h, l = 3h, L, W = A i h

FIGURE 3

 3 FIGURE 3 Streamlines for 3-dimensional steady flows, A i = 8, L = 30h. u c at Re c = 103, u FB2 and u BBr1n are both evaluated at Re = 174 [Colour figure can be viewed at wileyonlinelibrary.com]

FIGURE 4

 4 FIGURE 4 Skin friction for the postcritical BBr1n branch of the sudden expansion with A i = 6 at Re = 135. (Left) our study, made with ParaView with line integral convolution; (Right) extracted from the work of Chiang et al, 27 half-domain for roof and floor plane [Colour figure can be viewed at wileyonlinelibrary.com]

FIGURE 6 5 FIGURE 7

 657 FIGURE 6 Slice and streamlines for the singular tangent operator specific vectors at first steady primary bifurcation (PB1). Bifurcation mode 𝚽, left bifurcation mode 𝚿, and specific solution W. Case of A i = 5

FIGURE 8

 8 FIGURE 8 Particular solution W PB1 . Velocity components and pressure for E = 3A i = 5L = 30. Contours for W u are 5 • 10 -4 , 7 • 10 -5 , and 5 • 10 -5 . Red: positive; blue: negative

FIGURE 9

 9 FIGURE 9 Velocity comparison of flows on the post-bifurcated branches BBr1n (post • at Re = 444) and BBr2n (post □ at Re = 315) for the case E = 3, A i = 8, L = 55h. Up: crop of 3 planar slices, at z -= -25, z 0 = 0, and z + = 25, of velocity magnitude at different spanwise coordinates. Middle: front view of the u x velocity component signed contours: 0.2 (red) and -0.2 (blue). Down: overview of those contours with transparency for the positive value

FIGURE 10

 10 FIGURE 10Critical Reynolds number associated to primary bifurcation for several A i using ELMER coupled with the asymptotic-numerical method. In black circles, values from the literature [Colour figure can be viewed at wileyonlinelibrary.com]

FIGURE 11 FIGURE 12

 1112 FIGURE 11 Bifurcation modes 𝚽 PB2 □ (Left) and 𝚽 PB3 ▴ (right). Velocity components and pressure for E = 3, A i = 8, L = 55h. Contours for 𝚽 u are 10 -3 • 10 -4 , 10 -4 . Red: positive; blue: negative

FIGURE 14

 14 FIGURE 14 Synthetic bifurcation diagram. Critical Reynolds numbers are respected. Primary bifurcation type and symmetry state are given for each concerned branch

FIGURE 15

 15 FIGURE 15 Secondary bifurcation modes 𝚽 PB1SB1 (left) and 𝚽 PB1SB2 (right). Velocity components and pressure for E = 3A i = 8L = 55h. Contours for 𝚽 u are 10 -3 • 10 -4 , 10 -4 . Red: positive; blue: negative

FIGURE 16 FIGURE 17

 1617 FIGURE 16 Secondary bifurcation modes 𝚽 PB2SB1 (left) and 𝚽 PB2SB2 (right). Velocity components and pressure for E = 3A i = 8L = 55h. Contours for 𝚽 u are 10 -3 • 10 -4 , 10 -4 . Red: positive; blue: negative

TABLE 1

 1 Meshes for the case E = 3, l = 3h, and ratios A i = W∕h. Length L, number of H8 elements for entrance (Nh×NW×Nl) and body (NH×NW×NL), half-bandwidth (HBW), memory space, and CPU factorization time required by MUMPS are given. (v4) stands for mono-threaded computation using the MUMPS V4 library, otherwise the MUMPS V5 and OpenBlas libraries are used

			Number of H8 Elements	Number of	Number of
	A i	L	Entrance	Body	Nodes	DOF HBW	RAM	LU Thread
	1	30h	12×12×36	36×12×360	179 725	718 900	1131	14 GB	330 s 8T
	1	60h	22×22×66 66×22×1320	2 070 575	8 282 300	3611 368 GB 1550 s 36T
	1.25 30h	16×20×48	48×20×480	512 085	2 048 340	2415	75 GB	940 s 8T
	1.5	60h	12×18×36	36×18×720	515 755	2 063 020	1653	58 GB 1800 s 1T (v4)
	2	30h	12×24×36	36×24×360	345 625	1 382 500	2534	45 GB 1800 s 1T (v4)
	2	50h	16×32×48	48×32×800	1 322 145	5 288 580	4466 256 GB 1100 s 36T
	2.5	30h	12×30×36	36×30×360	428575	1 714 300	3416	68 GB 1000 s 8T
	3	30h	8×24×24	24×24×240	156 025	624 100	2162	17 GB	550 s 1T (v4)
	3	50h	14×42×42	42×42×700	1 323 239	5 292 956	6470 235 GB 1300 s 36T
	3.5	30h	8×28×24	24×28×240	180 989	723 956	2220	22 GB	860 s 1T (v4)
	3.75 30h	10×37×30	30×37×300	367 118	1 468 472	3414	59 GB	855 s 8T
	4	30h	10×40×30	30×40×300	396 101	1 584 404	3414	66 GB 1000 s 8T
	4	40h	14×56×42	42×56×560	1 757 463	5 643 684	6570 340 GB 1700 s 36T
	5	30h	10×50×30	30×50×300	492 711	1 970 844	3493	91 GB 1500 s 8T
	5	40h	14×70×42	42×70×560	1 757 463	7 029 852	6737 476 GB 2900 s 36T
	6	30h	8×48×24	24×48×240	305 809	1 223 236	2667	48 GB 3000 s 1T (v4)
	6	30h	15×90×45	45×90×450	1 953 406	7 813 624	9093 604 GB 3600 s 36T
	7	30h	9×63×27	27×63×270	502 912	2 011 648	3855	96 GB 1800 s 8T
	7	30h	15×105×45 45×105×450	2 275 396	9 101 584 10 473 749 GB 4900 s 48T
	8	30h	8×70×24	24×70×240	486 989	1 947 956	3785	89 GB 1500 s 8T
	8	55h	12×96×36	36×96×714	2 417 725	9 670 900	7647 563 GB 4300 s 48T
	9	30h	8×72×24	24×72×240	455 593	1 822 372	3867	83 GB 1500 s 8T
	9	40h	12×108×36 36×108×360	1 990 885	3 981 770	8535 592 GB 3300 s 36T
	10	30h	8×80×24	24×80×240	505 521	2 022 084	4267	96 GB 1800 s 8T
	10	40h	12×120×36 36×120×480	2 210 065	8 840 260	9423 662 GB 4300 s 36T

TABLE 2

 2 Critical Reynolds number associated to the first primary steady bifurcation for several A i . "FL" stands for finer mesh or longer body, "NS" stands for not searched

	This Study
	A i =

W∕h Re cPB1 Re cPB1 (FL) Literature Re c

	1	( > 450)	734.7	Tsui and Wang 29	> 300
	1.25	438.7	NS
	1.5	333.3	NS
	2	241.7	241.7
	2.5	196.1	NS
	3	170.6	170.2	Tsui and Wang 29	171
	3.5	149.9	NS	Chiang et al 27	145
	3.75	142.5	NS
	4	136.7	136.2
	5	121.2	120.7
	6	112.7	111.5	Schreck and Schäfer 26	113
	7	106.9	106.0
	8	103.1	102.4
	9	100.5	99.4
	10	98.3	97.3	Medale and Cochelin 14	98.3

TABLE 3

 3 Correspondence between u max and 3 2 u mean for a uniform velocity inlet in a rectangular duct. Values higher than A i = 3 are extracted from table 4 in the work of Chiang et al[START_REF] Chiang | Side wall effects on the structure of laminar flow over a plane-symmetric sudden expansion[END_REF]

	A i	Ratio U max ∕ 𝟑 𝟐 U mean
	1.5	1.576
	2	1.488
	3	1.242
	3.5	1.209
	3.75	1.195
	6	1.117

TABLE 6

 6 Critical Reynolds numbers for the steady bifurcations using finer or longer mesh with the ratio A i = 8. NS stands for not searched

	FB1	FB2	FB3	FB4
	PB 102.36 203.87 297.59 569.01
	SB1 442.06 315.08 653.05	NS
	SB2 602.93 602.20	NS	NS

TABLE A3

 A3 CPU time in seconds and memory usage to perform one ANM step using truncation order N = 20 for a 2-dimensional test case. Direct solver is MUMPS v4 with one thread. Iterative solver is BiCGStab with ILUt preconditioner with the notation [min, max] for memory usage

			MUMPS	BiCGStab+ILUt
	Mesh	DOF Time	RAM	Time	RAM
	2D-Q4-M1	73 155	1.4 s	1.04 GB	5.3 s	[0.21, 0.97 GB]
	2D-Q4-M2	289 155	6.9 s	3.69 GB	56.8 s	[1.77, 5.13 GB]
	2D-Q4-M3 1 149 699 36.1 s 14.62 GB 721.1 s [12.89, 26.42 GB]

TABLE A4

 A4 CPU time in seconds to perform one ANM step using truncation order N = 20 for a 3-dimensional test case. Direct solver is MUMPS v4 with one thread. Iterative solver is BiCGStab(l) using polynomial parameter l = 4 with no preconditioner

			MUMPS	BiCGStab(l = 4)
	Mesh	DOF	Time	RAM	Time	RAM
	3D-H8-M4	624 100	524.7 s 18.0 GB 1865.8 s	4.3 GB
	3D-H8-M5 1 197 964 1673.9 s 40.0 GB 4716.2 s	8.4 GB
	3D-H8-M6 1 584 400 2979.5 s 58.9 GB 7305.7 s 10.9 GB

ACKNOWLEDGMENTS

The authors would like to thank Obame Endamne Sylvain for the linear solvers comparisons proposed in Appendix A.1.

APPENDIX A : LINEAR SOLVERS COMPARISONS

In a context of large-scale problems analysis, question arises to use the most appropriate solver. It is known that iterative solvers are very efficient and computationally cheap compared with direct solvers. Nevertheless, it is briefly shown that ANM-based methods are tailored for direct solver because a sequence of algebraic systems is solved using the same matrix operator. Thus, a short comparison of direct and iterative linear solvers is proposed in this section. In the software ELMER FEM, the following linear solvers have been tested:

• Direct: UMFPACK [START_REF] Davis | Algorithm 832: UMFPACK v4.3-an unsymmetric-pattern multifrontal method[END_REF] and MUMPS 41

• Iterative: GMRES, [START_REF] Saad | GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF] CGS, [START_REF] Sonneveld | CGS: a fast Lanczos-type solver for nonsymmetric linear systems[END_REF] BiCGStab, [START_REF] Van Der Vorst | Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems[END_REF] TFQMR, [START_REF]A transpose-free quasi-minimal residual algorithm for non-Hermitian linear systems[END_REF] GCR, [START_REF] Eisenstat | Variational iterative methods for nonsymmetric systems of linear equations[END_REF] using the following iterative preconditioners: ILU(n) with n ∈ [0, 9] and ILUt. 63

A.1 Linear solvers' behavior

Firstly, comparisons have been performed using a 2D sudden expansion. Same definitions of the geometric parameters presented for the 3D test case have been used. Hence, the following geometric features are: h = 10, E = 3, l = 3h, L = 30h (see also the work of Guevel et al [START_REF] Guevel | Automatic detection and branch switching methods for steady bifurcation in fluid mechanics[END_REF] ). A Q4/4C mixed finite element [START_REF] Zienkiewicz | The Finite Element Method, Fluid Dynamics[END_REF] has been used in ELMER. The ANM parameters are N = 30 and η = 10 -14 . The same trivial initial solution at Re = 0 has been used. The ANM continuation procedure was stopped at Re ≈ 81.

A.1.1 Direct linear solvers

A simple continuation test has been performed to compare 2 direct linear solvers available in ELMER. The number of ANM steps and total CPU time are recalled in Table A1. In our study, MUMPS is faster than UMFPACK.

A.1.2 Iterative linear solvers

All combinations between the 5 iterative linear solvers, ie, GMRES, CGS, BiCGStab, TFQMR, GCR, and 11 preconditioners, ie, ILU n with n ∈ [0, 9] / ILU t , have been tested for one ANM step with a linear iteration tolerance parameter TOL = 10 -8 . The following conclusions have been reached: BiCGStab and CGS using either ILU (5) or ILUt were the most efficient in term of computational time.

Then, a simple ANM continuation procedure has been performed to compare those 2 iterative linear solvers using either ILU (5) or ILUt. Results are recalled in Table A2.

Using ILU( 5) the number of linear iteration per ANM step grows from 3 iterations, for the first step, to 14 iterations when reaching Re > 70. Twelve ANM steps are necessary to reach Re ≈ 81 with a time per step between 20 to 60 seconds. In this case, BiCGStab+ILU(5) was the most efficient combination.

When using ILUt, only one linear iteration per ANM step is performed for both BiCGStab and CGS linear solvers. Seven ANM steps are performed to reach the Reynolds number Re ≈ 81 with approximately 16.5 seconds per ANM step. In this case, results are comparable for both combinations.