1. Introduction and motivation
In the present work, we have focused on a symmetry breaking bifurcation for the solutions of a fluidstructure interaction (FSI) problem in a planar contractionexpansion channel. We have chosen to tackle this problem because it leads to complex and rich dynamics for relatively small Reynolds numbers. In fact, the model is characterized by the socalled Coandă effect [coandaahmed], which in this setting causes the onset of an asymmetric flow that attaches to the nearby wall, despite the symmetry of both the geometry and the boundary conditions.
In particular, the application that motivated our analysis is a heart disease known as mitral valve regurgitation, which occurs when blood flows abnormally from the left ventricle to the left atrium because of the defective closure of the mitral valve (see Figure 1). In particular situations, this anomalous flow also exhibits the aforementioned wallhugging behavior, leading to incorrect readings of the blood volume through echocardiography, which is used to assess the severity of this pathology [echoassestment]. Since a channel with a sudden expansion preserves the essential properties of the blood flow through the two heart chambers, it is well suited to characterize the main features of this phenomenon.
In fact, the test case under consideration undergoes a bifurcation due to the nonlinearity of the NavierStokes system, and this results in a change of both the number of equilibrium solutions and their nature for a small variation of the Reynolds number . For low , the symmetric solution is unique and stable [Serrin, Shapira, Sobey]; however, as increases, the inertial effects become progressively more important, breaking the symmetry of the solution for a certain critical value [Battaglia]. Below this critical value, a physically unstable configuration of symmetric flow and two stable asymmetric ones coexist. The existence of multiple solutions for the same parameter values, which is typical of branching phenomena, gives rise to a bifurcation pattern known as supercritical pitchfork bifurcation.
The Coandă effect has been widely studied in the context of pure fluid dynamics (see Section 2); thus, we propose a generalization of the previous works by facing this phenomenon in a multiphysics context. To provide a simplified description of mitral regurgitation, we set up the FSI problem, which deals with the mutual inﬂuence between a deformable structure and the fluid flow around it. The effects of this interaction can be significant or even predominant due to the stress field exchanged at the interface. To this end, we used the Arbitrary LagrangianEulerian formulation [aleformulation], which is needed to solve the problem on a fixed reference configuration. Our investigation aims at understanding how the introduction of the structure can modify the bifurcating behavior. For this purpose, we carried out a comparison between the multiphysics case and the one where a nondeformable rigid body replaces the elastic structure. We further deepened our investigation by considering large scale deformations. Moreover, towards a better description of the physics involved, we introduced a nonlinear constitutive relation for the structure, i.e. the Saint VenantKirchoff model.
The study of the bifurcation implies being able to obtain solutions that share the same qualitative properties (which are said to belong to the same branch) while changing the value of the parameter, which in our case will be the kinematic viscosity of the fluid, strictly linked to the aforementioned Reynolds number. This task is not trivial since, after the bifurcation has occurred, we lose the uniqueness of the solution [seydel]. Consequently, we employed an algorithm combining the GalerkinFinite Element method [quarteroni2017numerical], the Newton method [ciarlet2013linear], and a continuation method [continuation]. This approach allows us to obtain a sequence of solutions associated with the same prescribed branch.
Despite this, since the highfidelity models usually require unaffordable computational costs due to fine mesh discretization and repeated evaluations, we also focused on the application of Reduced Order Models (ROM) [benner2017model, morepas2017]. In particular, we employed the Reduced Basis (RB) method [hesthaven2015certified, patera07:book], which has been already used to trace the bifurcating behavior of many problems in fluiddynamics [HESS2019379, pichistrazzullo, Pitton_2017], with further applications to other fields such as quantum mechanics [pichiquaini_old], computational mechanics [pichirozza] and biomathematics [coanda]. This method is based on projecting the nonlinear governing equations on a reduced space of much lower dimension. The latter is generated by a small set of RB functions, which can be obtained employing the Proper Orthogonal Decomposition (POD) [POD].
Our results suggest that although the bifurcation is caused by the nonlinearity of the constitutive equation for the fluid, the presence of the solid and the particular constitutive relation used to model it, scale the bifurcating behavior. Moreover, when dealing with large scale deformations, the bifurcation point translates into the parameter interval. This effect was observed for both the linear and nonlinear solid constitutive relation; however, it is amplified in the first case.
The work is organized as follows:

Section 2 provides the state of the art related to the Coandă effect in planar contractionexpansion channels;

Section 3 introduces the FSI problem, providing the variational formulation for the steady case;

Section 4 describes all the ingredients to tackle numerically a generic bifurcation problem held by a nonlinear PDE: the high order approximation, the Reduced Basis method, and the Newton method.

Section 5 shows the results for the FSI problem concerning the rigid body model, the linear elastic model, and the Saint VenantKirchoff nonlinear model.
2. Related Work
As anticipated in the introduction, the object of our investigation is a symmetrybreaking bifurcation for a planar incompressible flow through a contractionexpansion channel.
The appearance of an asymmetric behavior for laminar flows, undergoing a sudden section change, is a problem that has been widely studied both experimentally and numerically [Cherdron, Sobey, Tritton]. Such phenomenon is caused by the Coandă effect, which concerns velocity fluctuations that lead to a transversal pressure gradient capable of maintaining the flow asymmetry [Tritton, Wille].
Cherdron et al. [Cherdron] presented a first large experimental analysis that uses flow visualization and laserDoppler anemometry. This study explains the origin of the phenomenon through the disturbances generated at the corners of the expansion, and subsequently amplified by shear layers. The latter plays a fundamental role, as the Moffat eddies [Moffatt] forming inside them are interdependent due to confinement in the transversal direction.
Another important contribution was provided by Sobey & Drazin [Sobey], who dealt with the problem using an asymptotic, numerical and experimental approach. They understood the strong connection of the problem with bifurcation theory, being able to deduce the presence of a Hopf bifurcation characterized by the existence of a timeperiodic solution, whose
have been successively estimated by the work of Quaini et al.
[Quaini]. The very nature of the temporal dependence of the solution was then investigated by Fearn et al. [Fearn], who managed to trace the phenomenon to threedimensional effects through an experimental investigation.The results of all these studies highlighted the onset of the Coandă Effect for relatively low Reynolds numbers. This can be partially justified by the fact that although the Coandă effect occurs more strongly in turbulent flows, it is amplified in the twodimensional setting [Wille]. However, these early works were carried out working with fixed expansion ratio , or with few values. Battaglia et al. [Battaglia] and Drikakis [Drikakis] then generalized the previous discussion by studying the relation between and . Their results showed that the critical Reynold number increases as the expansion ratio decreases; however, subsequent investigations, that include a more complex case study, show a nontrivial behavior of the bifurcation with respect to the expansion ratio [pichi2021artificial]
The physical interpretation of the laminar flow’s dynamics inside the contractionexpansion channel was provided by Hawa & Rusak’s work [Hawa], who have shown that, for low Reynolds numbers, the flow is stabilized by viscous dissipation. As increases, this stabilization is not sufficient to counteract the destabilization by the convective effects upstream of the expansion, and precisely this imbalance leads to the instability of the symmetric solution.
Subsequent investigations mainly focused on different kinds of generalizations, such as the constitutive relationship [Mishra], geometry [Mizushima], or the link between numerical schemes and the critical Reynolds value [Quaini].
The analogy between the flow in the contractionexpansion channel and the anomalous jet due to mitral regurgitation (see Section 1), have prompted recent works [Pichi, coanda, Pitton_2017] to create a connection between the medical literature [Chandra, Little, Vermeulen] on this pathology and the numerical modeling of the cardiovascular system. In particular, these works have tackled the problem of obtaining efficient and realtime simulations, using reduced order models, to be exploited in a reallife scenario. For example, the work of Pitton et al. [coanda] allowed to reproduce the Coandă effect in a mock heart chamber.
The big limitation of these previous works lies in the fact that they have faced the problem by solving the NavierStokes system in the contractionexpansion channel with fixed straight walls. An important step towards more realistic test cases was made by Hess et al. [quainihess], by studying a contractionexpansion channel with curved walls.
Our work seeks to carry out this attempt at generalization. In fact we are not aware of studies on the effect that the coupling of the fluid with an elastic solid can have on the bifurcation. Therefore we have decided to tackle the problem in the context of fluidstructure interaction. In addition to providing a replicable benchmark for future investigation in the multiphysics context, we tested a reconstruction using the RB method, which has been proven to be effective for FSI problems [Bal2016, Ballarin2017, Nonino, fluids6060229].
3. Problem formulation
Let us consider the abstract strong form of a parametric PDE, which reads: given a parameter find such that
(1) 
being the functional space where we are seeking the solution , and its dual space. In order to retrieve a variational formulation of the problem, we introduce the following variational form:
where is a test function and denotes the duality pairing between and .
Then, the weak formulation of equation (1) reads: given find such that
(2) 
We are interested in studying the behavior of the solution as the value of the parameter varies. In particular, if we assume that is continuously differentiable, and we denote by its Fréchet derivative with respect to , then the following important result holds true [ma2005bifurcation]:
Theorem 3.1.
Let be a known solution of equation (1). Let be two balls of radius and around and , respectively. Let be a map. If is bijective then, there exist and a unique solution such that
Therefore, when the aforementioned assumptions fail, it can happen that even a slight a variation of the parameter leads to changes in the number of the existing solutions for a given value of the parameter [caloz1997numerical, Brezzi1980]. In particular, this means that multiple solutions could branch off from a reference solution . The value of the parameters for which this happens will be called bifurcation points of equation (1), and can be defined as done in [ambrosetti_malchiodi_2007]:
Definition 3.1.
A parameter value is a bifurcation point for equation (1) if there exists a sequence with such that:

;

.
The manifold represented by the solutions with the same qualitative behavior will be named branch and denoted by . The objective of our investigation is the reconstruction of the entire solution manifold, given by the union of all the branches:
(3) 
In particular, this is fundamental when trying to recover the socalled bifurcation diagram. The latter is an important tool, typical of the study of dynamical systems, which allows us to investigate the qualitative behavior of the solutions as a function of the parameters’ set taken into consideration. The idea is to trace the evolution of a significant scalar quantity , which can be considered as the output of the dynamical system. Typical examples in this regard are the norm of the solutions, or a pointwise evaluation of one of the solution’s components. The presence of multiple solutions corresponding to the same parameter value results in multiple branches in the bifurcation diagram.
In the following, we want to apply the former discussion to the fluidstructure interaction problem. This will be done by retrieving both the strong and weak form for the particular problem at hand. However, we need to first introduce the Arbitrary LagrangianEulerian formulation, which will be used to solve the problem encoding coherently the fluid’s flow and the displacement of the structure.
3.1. Preliminaries: configurations and ALE description
In this section, we are going to introduce the Arbitrary LagrangianEulerian (ALE) formulation [aleformulation] which is widely used for the description of the fluid motion in FSI problems [fsibook]. In particular, let be the physical domain, which in our case will be a subset of the twodimensional Euclidean space . We assume that the solid resolution domain and the fluid one are such that and they do not overlap^{1}^{1}1In the following, we will use the subscript (respectively ) to refer to solid (respectively fluid) ’s phase quantities..
We know that for a generic FSI problem we have to deal with a moving domain. As regards the description of the solid phase, the problem is solved by defining the equations in material form, i.e. in the reference domain . Typically the problem does not concern fluids as they adapt to the domain by assuming its form. However, in fluidstructure interaction problems, the fluid resolution mesh needs to be in motion in order to be consistent with the structural deformation. One possible approach to deal with this issue is to use the already mentioned ALE formulation, which consists in defining an intermediate reference system on which the fluid equations are pulled back. Therefore, it is necessary to introduce an additional field , which represents the displacement of the fluid mesh we are modeling; this displacement must necessarily satisfy a Dirichlet condition on the boundary with the solid resolution domain :
being the solid displacement. The crucial task is to extend the displacement inside the domain . To do so, as done in [aleformulation], we introduce a fictitious displacement which follows the Dirichlet datum prescribed on the boundary. This extension is said to be arbitrary because it only concerns the reference system used to describe the fluid motion, and not the motion itself. This led us to consider an extension problem, among which we have decided to use the harmonic extension. The latter consists in solving a Poisson problem with inhomogeneous Dirichlet boundary conditions on the fluidstructure interface:
(4) 
The harmonic extension has been formulated in the fixed reference configuration , which is the inverse image of the Eulerian fluid resolution domain through the ALE map :
(5) 
The time derivative in such a system, denoted with
, can be computed for any scalar or vector quantity
through the Reynolds transport formula:(6) 
3.2. Strong formulation for the FSI problem
In this section, we are going to introduce the fully coupled formulation for the FSI problem. Since in this setting the ALE description is used for the fluid and the Lagrangian one for the solid, we will need two distinct reference domains, namely one for the solid and one for the fluid (see Figure 2 for a visual representation). The unknowns in this setting are: the global displacement , the fluid velocity , and the pressure .
By combining the equations for the solid phase (momentum equation), the fluid phase (NavierStokes system) and the mesh displacement (harmonic extension) in a single system, we obtain the full formulation for the FSI problem, namely:
(7) 
where and are respectively the solid and fluid densities, whereas and represent the external volume forces for the two physics. Moreover, we denoted with the
Cauchy stress tensor
and with the second Piola tensor; their expressions are provided by means of the constitutive relations for Newtonian incompressible fluid and linear elastic solids, which read:(8) 
and are a couple of materialdependent quantities for the solid known as Lamè parameters, whereas is the linearized strain operator, which is defined as
(9) 
is instead a fluiddependent quantity, known as kinematic viscosity, which assumes a fundamental role in our analysis since it will be the bifurcation parameter under investigation.
In Equation (7) we have introduced a notation clarifying the coordinate system we use in each equation: the divergence in the momentum equation for the solid reads , because it needs to be taken w.r.t. the reference variable , whereas the time derivative in the NavierStokes equation is referred to as because it is the ALE time derivative which has been already introduced in Equation (6). The same notation is used also for fields, therefore and describe the same field, but the difference is in their dependence as a function of or ^{2}^{2}2It should be clear that we are considering the same field but expressed in different variables. Exactly because of this, one can obtain a ALE description of an Eulerian field or vice versa, using composition with or with .
There is still a problem with the former set of equations. In fact, even if the Eulerian time derivative has been substituted by the ALE one in the NavierStokes system, we have not pulled back those equations on the reference domain . In order to do so, we introduce:
(10) 
which are the gradient of the ALE map (see Equation (5)) and its determinant. We can perform a change of variables by exploiting the Piola transformation rule to express the fluid equations on the reference configuration, which leads to:
(11) 
with being the representation of in the reference configuration, namely:
(12) 
In order to deal with a wellposed problem, we need to specify some proper coupling and boundary conditions. As concerns the former, we supplemented system (7) by means of the following conditions on :
(13) 
While for the latter we imposed the following:
(14) 
In the previous equations is the normal vector outgoing from . See Figure 2 for a sketch of the different portions of the boundaries.
As for now, everything is expressed in the reference configuration , therefore, since there is no need to distinguish between different coordinate systems, we will simplify the notation removing the hat symbol.
3.3. Variational formulation for the steady case
We proceed now recovering the variational formulation of the problem (11), which is needed in order to solve the system by means of a Finite Elements discretization. First of all, since the model will be solved in a steady setting, let us proceed to recover its strong form. From the solid’s momentum equation, we observe that the acceleration in the reference domain is equal to 0. Whereas, for what concerns the NavierStokes system, one can proceed by annihilating the partial derivative in time for the Eulerian framework before doing the pullback on the reference configuration. This leads to the following (steady) system:
(15) 
Furthermore, also the coupling condition needs to be changed accordingly:
(16) 
Finally, if we define the functional setting as and , then we can multiply each equation of system (15) by an appropriate test function belonging to the corresponding functional space and apply the divergence theorem. This leads to the following weak formulation:

FLUIDSTRUCTURE MOMENTUM EQUATION:
(17) 
FLUID INCOMPRESSIBILITY EQUATION:
(18) 
EXTENSION EQUATION:
(19)
The previous system can be approached numerically using numerous strategies. Our numerical resolution was built upon a monolithic method, which consists in solving the fully coupled problem, without having to iterate between the two submodules (which is instead done with the segregated approach [degrote]). The advantages of this kind of treatment are the achieved stability properties [stabilitymonolithic], and the possibility to solve the whole problem employing a global nonlinear solver such as Newton linearization. The price to pay is in terms of computational cost, which is why we combined this approach with the use of a modeling reduction strategy (Section 4.2).
4. Numerical Approximation
In this section, we will present the ingredients needed to tackle numerically a generic bifurcation problem held by a nonlinear PDE. The starting point will be represented by the socalled full order (or highfidelity) discretization, which aims at providing a high accuracy approximation to the solution of the PDE. With this regard, we will present only the GalerkinFinite Element (FE) method which has been used for our simulation. However, one could employ our framework and change the underline full order discretization strategy.
We will present the algorithm we used for the reconstruction of the bifurcation diagram. The chosen approach is branchwise because it aims at reconstructing the different branches separately. It also requires coupling the numerical resolution with a Continuation method[continuation], which fosters convergence to solutions belonging to the branch to be approximated.
The mitigation of the computational cost required by the presented methodology will be tackled using a reduced order modelling technique.
4.1. Full order approximation
The GalerkinFE method is a numerical discretization technique that allows solving various PDEs by resorting to their variational formulation. Here, we used this method to solve the FSI problem (Section 3.3), but it will also assume a prominent role in the reduced order setting (Section 4.2).
We start by considering a finite dimensional subspace such that . The main idea is to solve the weak variational formulation (2) onto this generic subspace, that means: given , find in s.t.
(20) 
We will refer to the above problem as Discrete Weak Formulation ().
Since has finite dimension, i.e. then , i.e. every element of can be written as a linear combination of the basis functions . Thus, we can express the solution as .
Because of the former consideration, one can easily prove that solving () is equivalent to finding the solution of the following nonlinear system:
(21) 
The nonlinearity in the model forces us to consider a nonlinear solver for its approximation, thus we relied on the Newton method.
First, we introduce a parametrized linear variational form by making use of the Fréchet partial derivative of at w.r.t. as
(22) 
Then, one suitably chooses an initial guess , and proceed with a sequence of iterations composed by:

Solving the finite dimensional linear variational problem:
that can be cast in the following matrix formulation:
being the Jacobian matrix and the highorder residual vector:

Update the approximated solution:
(23)
As stopping criterion for the former algorithm one can use a control on the highfidelity residual and/or on the norm of the variation , combined with a check on the maximum number of iterations.
We now want to use the former numerical resolution to perform a branchwise reconstruction of the bifurcation phenomena [Pichi]. This means following a single branch while modifying the value of the parameter we are dealing with. However, after the bifurcation had occurred, we lose the uniqueness of the solution, and therefore this task is not trivial. This is why we will implement a continuation method, that is a procedure that allows to generate a sequence of solutions associated with the same branch [continuation]. We consider a discrete version of the parameter space , and let us imagine starting from a certain value of the parameter , for which the corresponding full order solution is known. The idea is to get the solution for the next value of the parameter , starting from a suitable initial guess given by (, ) . In the following, we will be using the simplest version of this methodology which consists in choosing directly (, ) as initial guess. However, there are more complex strategies which can be exploited to automatically discover new branches [classic_deflation, pintore2019efficient]. We remark that a lot of attention must be paid to the choice of the discretization step; if it is too large, we risk not to capture the bifurcating behaviour, while a toosmall step causes a possible waste of computational resources, especially in regions of the parameter space far from the bifurcation points. The former considerations can be summarized in the following pseudoalgorithm for the reconstruction of a branch (see [Pichi] for more details):
4.2. Reduced order model
When dealing with a class of parametrized problem in the form:
(24) 
we may be interested in looking for the solution for all the possible instances of the parameters’ value, that is finding for all . Basically, we would like to approximate the solution manifold (Equation (3)). So far, we have seen how to tackle the problem within a highorder approximation framework (Section 4.1
). However, given that the full order models present a considerable computational difficulty (which scales proportionally to the number of degrees of freedom
), various Reduced Order Models (ROMs) [benner2017model, morepas2017] have been proposed to tackle the mitigation of the computational cost.In the present work we have decided to use a particular ROM known as the Reduced Basis (RB) method [hesthaven2015certified, patera07:book, QuarteroniManzoniNegri2015]. In particular, the RB methodology is divided into two stages:

Offline phase: highfidelity numerical simulations for different parameter values are performed. These values must be representative, and several methods specialize in how that selection is made. This phase represents the most expensive one from the computational point of view, but it has to be done only once. The information obtained must be subsequently processed to be transmitted to the next phase.

Online phase The precomputed quantities of the offline phase are used to solve with reduced computational complexity new instances of the parameterized problem. Ideally, the new resolutions’ complexity should not depend on and should require resources that can be managed by lowpower devices such as smartphones.
The previous strategy involves building a suitable approximation of the solution manifold by exploiting the information collected during the offline phase. In particular, the RB method consists in constructing this approximation through a finite dimensional space generated by a suitable basis^{3}^{3}3Here, we will denote with the subscript all the reduced basis quantities.. In the next section, we will introduce the optimality condition satisfied by , together with a strategy for its computation.
4.2.1. PODGalerkin method
As discussed before, a crucial point when applying the RB method is the construction of a proper basis during the offline phase. In the present work, we use the Proper Orthogonal Decomposition (POD) [POD], which can be used to compress information in different settings by extracting a low dimensional representation of a given dataset.
One begins by introducing a discretization in the parameter domain, namely , whose cardinality will be denoted by . If we solve problem (24) for a sufficiently rich discretized parameter space , we obtain an approximation of the highorder solution manifold:
(25) 
given by the vector space
(26) 
identified by the span of the so called snapshots (i.e. truth solutions w.r.t. ). Then, the idea is to find the space that minimizes the following quantity
(27) 
among all the subspaces of with fixed dimension .
From a discrete point of view, the former procedure is equivalent to carrying out a Single Value Decomposition (SVD) on the snapshot matrix
(28) 
By using the first left singular vectors as expansion coefficients w.r.t. , one obtains the set of RB functions which span .
The procedure in question must be carried out independently for every field of the unknown variable. Once obtained the reduced basis space , we perform a Galerkin projection onto the subspace itself. This strategy is equivalent to the one presented in Section 4.1, the only difference lays in the basis that we are using in this setting, which will be the one obtained via POD. Thus, the reduced problem reads: given , find s.t.
(29) 
The latter reduces, similarly to the full order model, to a system of nonlinear equations, which can be solved by means of the Newton method, as presented in Section 4.1. In particular, one can prove that the th iteration of the Newton method can be cast in matrix form as
(30) 
where is the matrix whose columns denote the coefficients of the reduced basis functions in terms of the full order basis functions . With regard to the branchwise reconstruction, one can obtain a reduced order version of Algorithm 1, by substituting the full order quantities with their reduced order counterpart:
5. Numerical results
In this section, we will present the results of our simulations for the fluidstructure interaction problem. In particular, we will compare the multiphysics case and the rigid structure one, where a nondeformable rigid body replaces the elastic structure. We will also present the results for the reduction strategy implemented through the RB method, whose accuracy has been tested on never seen parametric instances of the branch we aim at approximating. To investigate the effect of large scale deformations, we will use the Saint VenantKirchoff nonlinear model for the solid constitutive relation.
Finally, all the proposed models will be compared, obtaining important insights on how the bifurcation point is modified by the introduction of the structure and the constitutive relation chosen to model it.
5.1. Preliminaries
Let us start by noting that we have used centimetergramsecond (CGS) unit system for all the equations. Consequently, the various quantities have the following units of measurement:
Physical quantity  Unit of measure 

lengths  cm 
barye (Ba)  
stokes (St) 
The computational domain is the same for the FSI problem and the rigid structure one, and it is shown in Figure 3.
In the case of the rigid structure, the red region in Figure 3 is removed from the computational domain and is used only to impose the boundary conditions. On the contrary, for the FSI problem, it represents the reference domain for the solid (see Figure 2), where we solve the elastic problem . As already anticipated in Section 4.1, as highfidelity approximation we adopted a FE discretization on unstructured meshes composed by triangular elements (see Table 2).
Problem  Grid type  Number of cells  Number of vertices 

FSI  Unstructured with triangles  36419  18584 
Rigid structure  Unstructured with triangles  36118  18460 
The numerical resolution was carried out using the FeniCS library [fenics], with the additional use of multiphenics [multiphenics], which facilitates the definition of multiphysics problems, and of RBniCS [rbnics], which instead implements several reduced order modeling techniques for parametrized problems.
5.2. Rigid structure
The rigid structure case involves solving the NavierStokes system in Eulerian coordinates. In fact, since the domain is fixed there is no need to use the ALE formulation, so that the system of equations we solve is:
(31) 
The boundary conditions on and are of Neumann type and concern the imposition of the inlet and outlet stress values. In particular, the condition on is stressfree being homogeneous. Therefore, the flow is exclusively animated by the value of , which encodes the gradient of the stresses between the inlet and the outlet. This choice represents an alternative to a Dirichlet boundary condition which directly imposes the velocity value at the inlet. This last approach has been excluded to avoid lifting technique and make the ROM approximation easier. For our simulations we set , and , while the parameter adopted for the bifurcation analysis is the kinematic viscosity^{4}^{4}4Since we have written the equations in the system, and being , the values of the kinematic viscosity and those for dynamic viscosity coincide (except for the dimensions). of the fluid . However, since the channel’s geometry is fixed, and the only physical parameter that varies is the viscosity, the dimensionless bifurcation parameter is the Reynolds number. For what concerns the FE discretization, we used the TaylorHood () elements to avoid stability issues [quarteroni2008numerical].
We varied the kinematic viscosity within the parametric interval , composed by equispaced points (which is equivalent to adopt a step ).
During the study of the system, we expect a behavior similar to the one investigated in [Pichi, coanda], that is a supercritical pitchfork bifurcation. This consists in having a stable symmetric solution above a critical value , while, below this value, two stable wallhugging solutions and an unstable symmetric one coexist. To provide the bifurcation diagram shown in Figure 4, we have chosen as a representative scalar quantity the vertical velocity at a point of the channel symmetry axis, that is , with . It should be noted that the procedure implemented to follow both the unstable symmetric branch and the two stable asymmetric branches is the same (see Algorithm 1); the only difference, within our approach, consists in the direction of exploration of the parameter space . The two branches (a), (b) of Figure 4 refer to the asymmetric wallhugging jets, respectively to the upper and lower walls. In order to reconstruct these two branches, we proceed by calculating the solution for which is obtained starting from the trivial null guess, and then we decrease the viscosity through the step. However, due to the unstable nature of branch (c), for its reconstruction it is necessary to start with and proceed by increasing instead of decreasing it.
As shown in Figure 4, we were able to completely reconstruct the bifurcation diagram and estimate that, for this particular setting, the bifurcation occurs at . Figures 5 and 6 show the velocity and pressure profiles belonging to the branches and for values of . It is important to note that although the inlet stress is constant, the velocity at the exit of the expansion orifice varies according to the viscosity. This is highlighted in Figure 7, which shows the velocity profile for the extreme values of the parametric range, namely
5.3. FSI problem with linear elasticity
In this section, we introduce the results related to the fluidstructure interaction problem. We will start by replacing the rigid structure with a linear elastic one composed by two leaflets. We recall that the FSI system we need to solve is
(32) 
which is completed by the constitutive relations, boundary conditions and coupling conditions introduced in Section 3.2. An important remark concerns the imposition of the coupling conditions, In fact we decided to weakly impose the constraints and , by means of a Lagrange multipliers [YU20051] approach, which implies introducing two new unknown . The space discretization is obtained using second order Lagrange FE for and , while first order Lagrange FE for and .
The solid’s constitutive relation is completely determined by the value of the Lamè constants. These two parameters have the same dimension as stress, and for the results of this section, they assume the following values:
Physical quantity  Value 

(Lamé’s first parameter)  Ba 
(Shear modulus)  Ba 
The setting for the bifurcation analysis is the same as in Section 5.2, that is we consider the same values for the inlet and outlet stresses and we adopt the discretization , composed by equispaced points in the interval . However, given the more complex nature of the problem, we have also recovered a RB approximation of the FSI bifurcating phenomena using the approach presented in Section 4.2. First of all the offline phase was perfomerd by computing a FE solution to system (32) for each in . Figures 8 and 9 show the highfidelity solutions for solid’s displacement, velocity, and pressure obtained during this phase.
It should be noted that each solid displacement field is associated with a fluid displacement field which represents its harmonic extension. Indeed, if we consider the solution for and plot and together we obtain the result shown in Figure 10. The structure is clearly that of a Poisson problem, and in fact , defined on , spreads from the interface in all .
We can now proceed as done in Section 5.2 and retrieve a bifurcation diagram for the fluidstructure interaction problem. The first approach relies on considering the fluid phase and taking again as characteristic scalar output the vertical velocity , with . If we limit ourselves to the lower stable branch, the result is shown in Figure (a)a.
At this point, we want to investigate how the structure undergoes the bifurcation. This question arises spontaneously given that the symmetry breaking bifurcation is driven by the NavierStokes nonlinearity. The first attempt consists in analyzing the behavior of the maximum displacement of the upper and lower leaflets. In fact, as long as the fluid flow is symmetric, we expect that the upper and lower displacements exhibit the same trend. However, when the flow becomes asymmetric, we expect the solid deformation field to become asymetric as well. Moreover the lower leaflet undergoes greater deformations, as the jet of fluid directed downwards creates a depression in the lower part of the channel. Therefore, if we track as a function of we obtain Figure(b)b.
The result obtained is justifiable in the light of the analysis of the results of Figure 8. In fact, an asymmetric flow downstream of the expansion causes the stress field to be asymmetric as well. In particular, an attempt can be made to trace back the deformations’ profile to the gradient of the stresses on the two interfaces of each leaflet. For example, let us consider the pressure drop downstream the expansion, that can be roughly quantified through the following quantities:
The trend of these indicators as a function of the kinematic viscosity is shown in Figure (b)b. Hence, the pressure gradient presents a profile capable of justifying an asymmetry in the field of displacements of the solid phase. A more precise analysis would require calculating the result of the complete stress tensor (which also includes the viscous component linked to the velocity gradient) and integrating on each leaflet’s surface.
In order to retrieve a bifurcation diagram for the solid phase, one can exploit the asymmetric trend of the upper and lower deformations below the critical value . If we indicate with the upper leaflet and with the lower one (Figure (a)a), we can then measure the difference between the deformations with the indicator If we plot the trend of this quantity as a function of we obtain the bifurcation diagram shown in Figure (a)a. Therefore, in this case, we arrive at a result that is similar to the bifurcation diagram for the fluid phase (Figure (a)a).
We want to remark that developing a bifurcation analysis as a function of the Reynolds number is possible. We can associate a Reynolds number to each viscosity value if we consider the expansion orifice’s width as the characteristic length. This can be done using the maximum characteristic velocity at the inlet of the channel’s width section. As anticipated at the end of Section 5.2, this velocity is not constant as varies, and its trend is shown in Figure (b)b. At this point, for each viscosity value, we compute the Reynolds number as ; and through it, we express the bifurcation diagrams as in Figure 14, where we have .
5.3.1. Reduced order model results
Now we will present the results of the reduction phase for the problem above. To this end, we used the snapshots coming from the continuation method presented in the previous section. Using the notation of Section 4.2.1, this means that is the discretized interval [0.5,2], and . We then processed the snapshots to extract a reduced basis for each of the fields involved in the FSI problem. In this regard, we remind that there are 6 unknown fields^{5}^{5}5In particular they are }, however only the first three have a physical meaning. The others are the Lagrange multipliers introduced to impose boundary conditions and the mesh displacement which is needed for the ALE formulation. and consequently, for each of them, it was necessary to perform a POD. However, the reduced velocity space was also expanded using a supremizer enrichment to avoid stability issues [bal].
We show in Figure 15 the exponential decay of the normalized singular values w.r.t. the basis cardinality. In particular, we notice that the two displacement fields decay faster than the pressure and the velocity of the fluid. In fact, we recall that the problems associated with and are both linear, and this justifies the lower complexity of their solution manifold.
A valuable tool for evaluating the effectiveness of the basis extracted from the POD is the projection errors’ trend as a function of the basis cardinality. Since we are interested in reconstructing the solution for all the values of , we consider the following average error:
(33) 
with being the projection operator of the highfidelity solution onto the reduced basis space. The trend of this indicator for the different fields is shown in Figure (a)a. It should be noted that the projection error represents the best approximation error, so one cannot obtain an error lower than that (Cèa’s Lemma).
We have analyzed the average error trend for some basis cardinalities in the case of solutions obtained through the Reduced Basis method. In this case, the difference between the high fidelity solutions and the reduced ones was measured as follows:
(34) 
In Figure (b)b we observe that, as expected, the errors are higher than their direct projection counterparts. In fact, the nonlinearity of the system, which makes its resolution dependent on the Newton method and the chosen initial guess, and its singularity at the bifurcation point, ensure that the error decay does not have a monotone behavior; however, this decay remains exponential. Furthermore, the pressure field is the one with the highest reconstruction errors, while the velocity field can be reconstructed with the same accuracy as the displacements.
At this point, we have decided to work with 16 reduced basis for each one of the fields. In this regard, to correctly test the effectiveness of reconstruction by the reduced basis, we have refined the discretization of the parameter interval by considering equally spaced points. Following the implementation of the reduced continuation method (see Algorithm 2), the reduced solutions for the lower stable branch are obtained. We then compared them with the full order counterpart through the relative errors for each field taken into consideration (see Figure 17).
We observe that these errors reflect the behavior already highlighted through the average errors (Figure (b)b), therefore the pressure is reconstructed less accurately than the other fields. Moreover, we can observe a slight peak for the error trend in correspondence with the bifurcation value and the low viscosity regime. Thus it seems that the presence of the structure somehow helps the regularity of the parametric manifold, allowing for a better RB approximation of the bifurcating phenomena w.r.t. the one for the rigid case [Nonino].
For some application it may be important also to understand the effectiveness of the reconstruction as a function of the spatial coordinate . For this reason, we have also computed the componentwise relative error field, which assumes the following expression:
(35) 
Figure 18 shows these relative error associated with the four fields for . We can observe how each field presents criticalities in different regions of the domain. In particular, the two displacement fields are reconstructed with less efficiency in correspondence with the interface. On the other hand, the pressure has the highest reconstruction error between the leaflets.
Finally, the error for the velocity field traces the magnitude of the field itself, with a morphology similar to the one of the deflected jet, also presenting an additional criticality at the outlet of the channel.
5.4. FSI problem with the nonlinear Saint VenantKirchoff model
As we have seen in the previous section, introducing a structure with a linear constitutive law does not seem to modify the phenomenon of bifurcation undergone by the fluid phase.
However, it is necessary to consider the fact that the previous results were obtained with a high structural stiffness, which results in a smallscale deformation field. Precisely because of this amplitude, it is reasonable to expect that the bifurcation point is only mildly affected by the parametric interval taken into consideration.
Therefore, now we want to reduce the structural stiffness and analyze the effect this has on the bifurcation. At the same time we will modify the solid phase’s constitutive relation, exploiting a nonlinear hyperelastic model. This is done because, for many materials, the linear elastic model fails to describe high strains’ behavior accurately. The adjective hyperelastic in this context refers to the fact that the relation between stress and strain can be obtained by deriving a scalar quantity known as the strain energy density function. The simplest hyperelastic model is the Saint VenantKirchoff (SVK) one, which represents a straightforward extension of the linear elastic model.
If we define the GreenLagrange strain tensor as
(36) 
then the strain energy turns out to be a function of the deformation gradient through an explicit dependence on , that is:
(37) 
where and are the Lamé’s constants. Finally, to find the constitutive relation that expresses the Piola tensor is given by
(38) 
One can replace the two Lamé’s constants by another pair of independent parameters, whose values characterize the material under consideration. A very popular choice consists in using the Young’s modulus () and the Poisson’s ratio (), such that we can write
(39) 
In the following we will keep the Poisson’s ratio fixed to , while we will consider different values for the Young’s modulus (see Table 4).
Case  E 

A)  1.87e5 
B)  9.34e4 
C)  5.60e4 
D)  4.67e4 
First of all, we can compare the maximum displacements of the upper and lower leaflets for the different cases, as shown in Figure 19.
As can be easily observed, we have reduced the structural stiffness until reaching deformations of the same order as the domain’s dimensions. However, the analysis carried out in Figure 19 does not yet allow us to understand if the bifurcation behavior is only rescaled due to the differences in the values of the Young’s modulus or if there are differences in response depending on the case considered. Because of the former consideration, it makes sense to use as an indicator the difference between the maximum deformations of the upper and lower leaflets , which was used in the previous section to compute a bifurcation diagram for the solid phase. The result is shown in Figure 20. In the same figure, we also present the trend of its dimensionless counterpart , defined as
(40) 
We can observe, it can be observed, through the trend of the normalized indicator , that the bifurcation behavior of the structure is only scaled by the stiffness (Young’s modulus). A final analysis with this respect consists in comparing the various bifurcation diagrams obtained for the fluid phase through the scalar output used previously, namely , with . Also in this case, the normalization with respect to the maximum value allows us to conclude that even for the fluid phase the behavior is simply scaled with respect to the value of (see Figure 21).
5.5. Model comparison for the FSI problem
Finally, we compare the results for the different models considered. We recall that the differences concern the modelling of the leaflets for which we used a rigid body model (Section 5.2), a linear elastic model (Section 5.3), and the Saint VenantKirchoff nonlinear model (Section 5.4). However, a rigorous comparison between the linear and the nonlinear model requires that the same values of the parameters in the constitutive relations are used. Therefore, we decided to add a test case (A) for the linear model, for which and (corresponding to the test case (D) of the former section). In particular, we recall that these values result in a structural deformation of the same order of magnitude of the domain’s dimensions. Table 5 shows the values of the parameters for the different models compared in this section.
Let us start the comparison between the linear model (A) and the SVK model. In Figure (a)a we examine the maximum displacements of the two leaflets. The first observation concerns the fact that the SVK model corresponds to more significant deformations.
Case  

linear (A)  4.67e4  0.21 
linear (B)  2.89e5  0.44 
SVK  4.67e4  0.21 
Furthermore, we can observe an asymmetry in the deformation of the two leaflets for the SVK model even before of the bifurcation point, while this does not occur for the linear counterpart. Another interesting fact is the behavior in the postbifurcation regime. In fact, the linear model predicts that the asymmetry of the two leaflets’ behavior is more accentuated than the counterpart represented by the SVK model. This observation can be quantified by comparing the diagrams related to the scalar output (see Figure (b)b). Therefore, the linear model, to which a field of minor deformation is associated, responds to the bifurcation in a more accentuated way, due to its difficulty to model large displacements. Finally, to investigate how this difference affects the behavior of the fluid phase, we compare the bifurcation diagrams related to the vertical velocity . The result is presented in Figure 22.
As expected, changing the structure’s constitutive relation has important effects on the fluid flow field, and this also affects the bifurcation point. In particular, the linear model (B) (which corresponds to a more stiff structure) seems to share the bifurcation point with the rigid structure. However, in the postbifurcation regime this model gives rise to a vertical flow of lower intensity compared to the other cases examined. The decrease in structural stiffness in the linear case (A) leads to a delay of the bifurcation with respect to the kinematic viscosity. However, after the bifurcation, this model results in a vertical velocity that closely approximates the rigid structure’s case.
Furthermore, the SVK model delays the bifurcation with respect to the case of the rigid structure but anticipates it with respect to the linear (A) counterpart.
Finally, we can note that the SVK model differs more significantly from the linear model precisely in correspondence with the bifurcation, while in the postbifurcation regime, there is a similar behavior.
6. Conclusions and perspectives
The goal of this work was to analyze a complex phenomenon of bifurcation in the multiphysics context of fluidstructure interaction problems. In particular, we have adopted a very general algorithm for the study of the phenomenon at hand. Although the Coandă effect and the consequent investigation of the bifurcation it causes have already been addressed in previous works, we are not aware of studies on the effect that the coupling of the fluid with an elastic solid can have on the bifurcation.
In line with our main objective, we compared the rigid structure case with that of elastic structure. In a first test case, conducted using a linear constitutive relationship for the solid phase, we observed that the structure passively undergoes a bifurcation, as shown by the leaflets’ displacements. Furthermore, the presence of the structure seems to rescale the bifurcation diagram. On this test case, we also conducted a reconstruction using the Reduced Basis method. The peculiarity of this procedure concerns that the setting of the problem faced is very complex: we are dealing with a nonlinear system with seven unknown fields, and which undergo a bifurcation in the parametric range of interest. Nevertheless, we were able to extract a basis able to accurately approximate the solution manifold, as demonstrated by the trend of construction errors.
We then decided to deepen our investigation by considering large scale deformations. For consistency with the physical reality of the phenomenon, this was done by introducing a nonlinear constitutive relationship for the structure, i.e. the Saint VenantKirchoff model.
In this case, we have observed that modifying the structural stiffness does not affect the bifurcation point’s position. However, the comparison between the case of the rigid structure and that of the linear model showed that the solid model actually affects the bifurcation point. We have, in fact, noticed that the introduction of a solid with a nonlinear constitutive law seems to delay the bifurcation, and this delay is amplified in the linear case with large scale deformations.
During the conclusion of this work, we have identified many interesting topics to be deepened. In the present work, we have not investigated the stability of the solutions. In particular, we intend to develop an eigenvalue analysis that makes it possible to have a quantitative indicator for the bifurcation points; this will allow us to deepen the comparison between the presented test cases.
Moreover, we have used a monolithic approach for both the high fidelity problem and the reduced order one. A possible extension concerns the use of partitioned schemes based on the Chorin–Temam decomposition. This approach could also be integrated with a semiimplicit scheme for the treatment of the coupling conditions between fluid and solid [Nonino]. Furthermore, the reduced method we used depends on the number of degrees of freedom of the high fidelity problem. Consequently, one possible development concerns the implementation of hyperreduction strategies that guarantee a more efficient decoupling between the offline and online phases, in order to obtain a more significant speedup.
An interesting perspective results in the integration of our framework with machine learning techniques. In particular,
[Pichi]presents an approach for the efficient detection of the bifurcation points, based on a neural network that exploits the Proper orthogonal decomposition (PODNN). This strategy can be used to obtain an efficient decoupling between offline and online stages, as an alternative to classical hyperreduction strategies. However, there are also other approaches of possible investigation, such as the physics informed neural networks (PINNs).
Finally, it is interesting to modify the physics of the problem we have focused on. This can be done for example by considering a setting where the bifurcation is governed by the nonlinearity of the constitutive relationship for the solid problem. For example, we could exploit the buckling results shown for certain nonlinear hyperelastic models. A further generalization results in the investigation of bifurcation phenomena for compressible fluid dynamics. This last problem is considerably complex and of strong interest, as it could also allow to investigate the relationship between bifurcations and turbulence.
7. Acknoledgments
We acknowledge the support by European Union Funding for Research and Innovation – Horizon 2020 Program – in the framework of European Research Council Executive Agency: Consolidator Grant H2020 ERC CoG 2015 AROMACFD project 681447 ”Advanced Reduced Order Methods with Applications in Computational Fluid Dynamics”. We also acknowledge the PRIN 2017 ”Numerical Analysis for Full and Reduced Order Methods for the efficient and accurate solution of complex systems governed by Partial Differential Equations” (NAFROMPDEs), and the ”GO for IT” program within a CRUI fund for the project ”Reduced order method for nonlinear PDEs enhanced by machine learning”.
We would also like to show our gratitude to Professor Simona Perotto (Politecnico di Milano) for her useful insights and suggestions which helped shape both the research and the manuscript.