A unified approach for the solution of fluid-solid interaction problems with hyperelastic deformation in internal flows

In the single domain method for solving fluid–solid interaction (FSI) problems, a unified formulation is used for the entire computational domain. In such monolithic FSI solvers, all of the governing equations are solved simultaneously. In the present study, the single domain method is further extended to an interface-tracking FSI solver which accounts for mesh movement via an Arbitrary Lagrangian–Eulerian (ALE) description of the governing equations. The focus is on internal flow problems with large deformation. Pressure and velocity are selected as the dependent variables for both solid and fluid parts of the computational domain. A distinguishing feature of the proposed method is that the governing equations at the interface are discretized in a conservative manner. Interfacial boundary conditions are enforced via a pressure–velocity splitting method to convert the kinematic and dynamic conditions at the interface into pressure–velocity relations. A PISO-like procedure is used to solve the discretized equations. In order to evaluate the proposed solver, strongly-coupled FSI benchmark test cases are employed. The performance of the proposed method and computational results are also compared with those obtained by a conventional partitioned solver. The results show that the proposed solver provides more accurate results on a coarser mesh compared to the benchmark solutions. The proposed method is also capable of solving strongly coupled problems for which the partitioned solver fails to converge. K

Coupling techniques are divided into two major classes: partitioned and monolithic. Monolithic coupling often refers to the simultaneous solution of all governing equations in a fully-coupled fashion. In the partitioned approach, the equations of fluid and solid and mesh moving are solved sequentially. This enables the use of existing fluid and solid solvers, a significant motivation for adopting the partitioned approach. For several problems, the partitioned approach works well and is very efficient. However, it sometimes suffers from convergence difficulties, most commonly when the structure is light and the fluid is heavy, and when a weakly compressible fluid is enclosed by a very flexible solid. On the other hand, the main advantage of the monolithic approach is the robustness of the solver, especially for problems with high degrees of coupling, that is, strongly-coupled problems, due to implicit exchange of information between fluid and solid regions. This feature has motivated many studies, [3][4][5] as well as the current study, to adopt a monolithic approach for the coupling scheme.
Mathematical formulation of fluid and solid behaviors in FSI problems usually involves two separate, yet relevant, decisions regarding the selection of dependent variables and the coordinate systems. Having the fluid in mind, the desired variables are velocity and pressure, and an Eulerian approach is usually chosen to describe the fluid dynamics. On the other hand, deformation with respect to the initial shape, in the context of a Lagrangian framework, is preferred in solid mechanics. Therefore, different FSI solvers have been developed that are described as Lagrangian-Eulerian, [6][7][8][9] Lagrangian-ALE, 10,11 full Eulerian, [12][13][14][15][16] and full Lagrangian. [17][18][19] Using a unified coordinate framework for the entire solution domain makes it easier to describe fluid and solid domains as a single continuum with identical dependent variables. Unified approaches are known to be computationally more stable. For instance, in the full Eulerian method, one can define the velocity as the unknown variable in the unified governing equations, resulting in an implicit treatment of the coupling between fluid and solid parts of the solution domain. However, a special technique is required to capture the shared interface. 13,20,21 To deal with the large deformation of solid boundaries, a unified 'one-fluid' formulation in the full Lagrangian framework, 19 which benefits from interface-tracking strategies, is preferred. One of the issues with this latter approach is the high computational cost associated with the tracking of fluid particles.
Interface-tracking and interface-capturing techniques are two approaches to locating the fluid-solid interface throughout the solution process of FSI problems. Interface-capturing methods 13,22 employ additional equations to "capture" the interface location, and hence, eliminate the requirements for mesh motion; therefore they are the preferred choice for problems with a complex motion of the interface. On the other hand, in interface-tracking methods, 23 as the interface moves due to solid deformation, the spatial domain occupied by the fluid changes its shape in order to follow, that is, "track", the interface motion. For interfaces with reasonable geometric complexity, moving the fluid mesh to track a fluid-solid interface enables one to control the mesh resolution near that interface and obtain accurate solutions in such critical flow regions. According to this argument, the interface-tracking methods are more popular for addressing problems like flow passing a flexible pipe especially for high Reynolds conditions. Using identical unknowns for fluid and solid equations avoids large disparity between the variable values and improves the conditioning of the FSI solver. Having this in mind, some researchers have adopted the unified momentum approach, 12 where velocity appears in the momentum equation for the entire domain and other quantities such as deformation in solid and pressure in the fluid, are lagged in the iterative computational procedure. In another approach, called the "single domain method," fluid and solid domains are treated as a single continuum with different material properties but similar dependent variables, that is, pressure and velocity. [24][25][26][27] Among them, Karac 27 applied the finite volume version of this method to study interactions between fluid and its container under drop impact and reported a 2.5 times faster convergence rate compared to the traditional FSI solver with separate formulations for the fluid and solid parts. A similar approach was adopted by Greenshields et al. 26 to solve the case of pressure propagation in a flexible tube to prevent instability associated with strongly coupled problems. The results in their work reveal that this method can efficiently solve strongly coupled problems due to intrinsic implicit treatment of the interface condition.
The deformation of the solid is an inseparable part of an FSI problem. However, for solid materials such as elastomers, polymers, foams and biological tissues, special care is needed due to the large deformation and also the highly non-linear nature of the interactions. Among these materials, incompressible solids have been studied more. For example, Wheel 28 studied an incompressible linear elastic body using the SIMPLE method on a collocated grid with Rhie and Chow interpolation to calculate pressure and displacement in the solid domain. Such developments in computational solid dynamics were carried out for incompressible elastic materials by Bijenlonja et al. 29 Since for an incompressible material, = 0.5 and consequently λ tends to infinity, they concluded that employment of the commonly used constitutive relations could lead to numerical issues such as locking, where the calculated displacements were unrealistically small. To overcome these problems, another constitutive equation was suggested as = 2 − pI, where p is an 'equivalent pressure' in solid. Thus, incompressibility constraint can be enforced by ⋅ u = 0 where u denotes displacement vector. As mentioned before, the SIMPLE algorithm with Rhie and Chow interpolation, which is usually used for flow solutions, can be applied in solids to avoid non-physical oscillation of the pressure. In the case of large deformation of incompressible solids, Bijelonja et al. 30 developed displacement-pressure-based finite volume formulation for the simulation of hyperelastic materials undergoing large deformations. They used Mooney-Rivlin incompressible material model to simulate material deformation. To enforce the incompressibility constraint, an integral form of mass conservation was applied to the deformed geometry. Also, a segregated SIMPLE-like approach was used to solve the resulting set of nonlinear coupled algebraic equations with displacement and pressure as unknowns. Recently, a conservative fluid-like pressure-velocity-based formulation was proposed by Tandis and Ashrafizadeh in which large deformation of a weakly compressible hyperelastic solid was solved on an Arbitrary Lagrangian-Eulerian (ALE) framework. 31 The proposed algorithm in Reference 21 allows for the employment of various material models such as Mooney-Rivlin and Neo-Hookean constitutive laws and can be regarded as a proper basis for the development of a unified solver for the simulation of the interactions between a fluid and a hyperelastic solid.
Solution of FSI problems actually involves seeking unknowns so that they satisfy the governing equations for fluid and solid which are constrained by their relevant boundary conditions as well as kinematic and dynamic boundary conditions at the interface. In the previous development of the "single domain method," 25-27 the interface was treated as an internal face for which material properties were interpolated from fluid and solid sides. The type of interpolation, for example, linear, harmonic, and so forth, strongly affects the performance of the solver. 17 As a result of this approach, exchanging information between fluid and solid is treated implicitly which makes the algorithm more stable. In the authors' view, even though this type of treatment of the interface removes the need for the explicit exchange of the information between fluid and solid, proper enforcement of kinetic and dynamic boundary conditions is not necessarily ensured. In addition, the mass flux at the interface in this approach can be nonzero, which is obviously non-physical.
Providing implicit information exchange between flow and solid and hence more robust solution, various unified formulations have been suggested as discussed in the above-mentioned studies. The motivation for the current study, apart from avoiding issues arising from the available unified algorithms such as remeshing and mesh-to-mesh data transfer, 17,18,32 or solving additional phase equations, 33 is to develop a 3-D FVM-based unified formulation which allows for easy implementation within relevant open-source computational packages, for example, OpenFOAM, and alleviates stability issues that the available solvers suffer from. Additionally, this study is motivated to carry out the implicit data exchange between fluid and solid so that the mass and momentum are conserved while transferred across the shared interface.
The main objective of this paper is to develop a unified FSI solver for internal flow problems involving hyperelastic deformation of the solid boundary. Being computationally efficient for the interaction of weakly compressible fluid and linearly elastic solid, the "single domain method" [24][25][26][27] suffers from the limitation of small deformation. Therefore, the development of an algorithm for hyperelastic and potentially large deformations necessitates some special care for the deformation of the computational domains as well as enforcement of the boundary conditions at the fluid-solid interface. This work can be regarded as an improvement to the single domain method where the main purpose is to employ a unified formulation for the entire domain aiming to build a monolithic interface-tracking FSI solver which accounts for mesh motion in largely deformed domains. Here, the pressure-velocity formulation on the ALE framework across the entire domain is adopted and an iterative computational loop is used to deal with all non-linearities including those associated with geometrical deformations. In order to provide a unified pressure-velocity formulation, the solid is modeled as a hyperelastic material obeying the Mooney-Rivlin constitutive relation. 31 A PISO-like algorithm is then applied to solve the weakly compressible Navier-Stokes equations on 3D polygonal control volumes (CVs) across the entire continuum. Moreover, and as the main contribution of this paper, enforcement of the consistency conditions at the interface is handled via an innovative pressure-velocity formulation based on the conservation of momentum and mass. In order to validate the accuracy of the proposed solver, a challenging FSI benchmark case, that is, the movement of a pressure wave inside a tube, is employed. Finally, a two-dimensional internal flow test case is used to compare the performance of the proposed algorithm with a traditional partitioned solver.

MATHEMATICAL FORMULATION
In this section, a unified mathematical model is proposed for FSI problems. This model treats the fluid and solid as a single domain and provides the unified equations for the conservation of mass and momentum on an ALE framework, where pressure and velocity are the unknown variables. Afterwards, constitutive relation and equation of state for weakly compressible Newtonian fluid and Mooney-Rivlin solid are presented. The section is wrapped up with a discussion about fluid-solid coupling at the shared interface.

Governing equations
An arbitrary Lagrangian-Eulerian approach establishes relations between differential operators in the Lagrangian, Eulerian and 'referential' configurations whose associated coordinates are denoted by X, x, andX, respectively. Using the definition of 'referential' coordinateX, the time derivatives of a quantity in these three frameworks are related to each other as follows 34 : Whereû denotes the referential velocity. Using the referential framework, governing equations are formulated on a CV, which moves with respect to the continuum. Hence, conservation of mass and momentum in the absence of body forces and in a time-dependent continuous domain Ω (t), enclosed by the boundary (t), are described as follows: Here,û is the velocity of boundaries of the domain Ω (t), is the density of volume V, is Cauchy stress tensor, u is the velocity of the continuum, and n is the outward pointing unit normal to the control surface S. It should be noted that in the context of FVM, Ω (t) and (t) denote the CV and its control surfaces, respectively. Therefore,û is actually representing mesh velocity at control surfaces associated with the CVs. According to the space conservation law, the change of the volume Ω (t) and the velocity of the control surfaceû must satisfy the following constraint: This relation is actually used for the calculation ofû.

Constitutive laws
The distinction between fluid and solid, as two types of continuum, is manifested in their different constitutive laws. The constitutive equation for a weakly compressible Newtonian fluid, based on primitive variables, that is, velocity and pressure, is as follows: p and are pressure and dynamic viscosity of the fluid. For the solid part, traditionally, the deformation is used as the unknown variable. In order to develop a unified pressure-velocity formulation, a constitutive relation similar to the fluid's relation is developed. For a hyperelastic solid, the following relation has been proposed in Reference 31: where I represents the identity tensor; B is the left Cauchy Green strain tensor; J indicates the Jacobian of the deformation gradient tensor, that is, J = ; C 1 , C 2 , and D 1 are material constants.

Equation of state
In a compressible continuum, pressure is related to the density and internal energy through an equation of state. In the case of weakly compressible material, the equation of state is usually modified using the bulk modulus K as shown below 26 : This relation can be further simplified by the definition of the speed of sound, c, in the medium.
Combining Equations (8) and (9), a simple linear equation of state is derived for a weakly compressible fluid: Here 0 and p 0 are the density and pressure of the fluid at rest. This linear equation is valid for density variations less than 0.01 0 .To mimic the form of equations commonly used in traditional compressible flow solvers, Equation (9) is expressed using the isothermal compressibility, ψ = ( p ) T , as follows: where C is a constant as defined below: Substituting Equation (11) into the first term of Equation (2) results in the following equation: It should be noted that the second term on the right-hand side of Equation (13) is not zero for moving CVs.

Boundary conditions at the fluid-solid interface
In the previous studies in the unified pressure-velocity formulation of FSI, [25][26][27] quantities such as velocity and material properties at the fluid-solid interface are linearly interpolated. However, as Karac 27 reported, the type of interpolation, for example, linear or harmonic, for evaluating the mentioned quantities at the interface is very important and can influence the performance of the algorithm. From the authors' point of view, an alternative way to address the issue of interpolation is to define pressure and velocity at the interface in such a way that ensures the consistency conditions. In this section, enforcement of the kinetic and dynamic boundary conditions at the interface in the context of pressure-velocity formation is explained in more detail. This section can be regarded as the major contribution of this study since it aims at introducing a physically consistent formulation of pressure and velocity at the interface in an FSI problem. More specifically, a new pressure-velocity formulation is developed at the interface in such a way that not only the consistency conditions, that is, kinematic and dynamic conditions, are enforced, but also the conservation of mass, that is, zero mass-flux across the interface, is satisfied. Additionally, and from the numerical point of view, such a formulation allows for applying a unified system of the equations for FSI problems where pressure and velocity are the only unknowns. The proposed formulation is obtained by implementing the projection method 35 as discussed below. In order to extract the boundary conditions at the interface Γ, each face at the interface is modeled as a zero-thickness CV (Δ = 0), as shown in Figure 2. The pressure and velocity at the interfacial CVs are calculated such that the consistency conditions and the conservation of mass and momentum are satisfied simultaneously. This means that pressure and velocity are sought throughout the entire domain, including at the interface, so that not only the conservation of mass and momentum are enforced for all internal CVs, but also consistency conditions are satisfied for all interfacial CVs.
The dynamic boundary condition for each face at the interface is, in fact, the traction balance at that face. This balance can be written as the momentum balance for the imaginary CV, indicated by C Γ , as follows: Decomposing the stress tensor into hydrostatic and deviatoric parts gives, Applying Gauss's Theorem transforms the integration from the imaginary CV to its surrounded faces which yields, The subscripts Γf and Γs refer to the fluid and solid sides at the interface, respectively; S Γf and S Γs are the surface vectors corresponding to the fluid and solid sides of the interfacial face, respectively. In order to bring velocity into the formulation, the deviatoric part of the stress tensor is substituted by the associated constitutive relations for the fluid and solid sides. For the fluid side, substitution from Equation (5) brings the velocity into the formulation directly. However, for the solid part, there is no explicit relation between the velocity and stress; so, velocity is introduced into the equation by adding and subtracting a velocity term as follows: Here, pseudo is a parameter that mathematically has no effect on the formulation and the final results, since the associated terms cancel out each other. However, choosing the value of (2 s + s ) Δt for this parameter provides the maximum consistent implicit contribution to component-wise discretization. 25,36 Additionally, the first and third terms in the first integral in Equation (17), that is, − 2 3 ⋅ u I and ( u) T , are usually small for weakly compressible flow and can be safely dropped from the equation.
Having a more careful look at Equation (17) reveals that while the pressure terms only refer to the value of pressure at the interface, the velocity terms appear in the form of divergence and gradient expressions. Considering the definition of gradient and divergence from the numerical perspective, this means that the velocities of neighbor nodes also play a role in this equation. Therefore, this equation is actually coupled to the momentum equations in the fluid and solid sides of the interface.
In order to develop a constraint on the interfacial pressure and close the system of equations at the interface, an idea borrowed from the projection method is employed: splitting the velocity into the intermediate velocity and pressure-driven velocity components. Such splitting is achieved via rearranging the discrete momentum equation in a way that gives the velocity at each cell in terms of velocity at the neighbor cells, the gradient of the pressure at the current cell, and other source terms: where a P represents the coefficient for each individual component of the velocity vector at cell P; the H P (u) operator contains all the matrix coefficients of cell neighbors multiplied by their corresponding velocities, plus source terms at cell P other than those from the pressure gradient. This description is a de facto standard description of the discrete momentum equation in order to distinguish the role of pressure and has been frequently used in the literature. Following this approach, the velocities at the interfacial faces, corresponding to the fluid and solid sides of the imaginary interfacial cell in Figure 2, are expressed as follows: This equation, indeed, represents two equations associated with the fluid and solid sides of the interface which are denoted by subscript r. Here, H(u) is defined similar to H P (u) except it is evaluated based on the momentum equation at the imaginary cells, that is, it represents the contribution of the velocity at adjacent cells and all source terms other than the pressure gradient at the imaginary cells corresponding to the region r; a refers to the diagonal coefficient of interface velocity in the discrete momentum equation at the interface, that is, Equation (17). Assuming that the velocity of the mesh motion at the interface isû Γ , the mathematical statement corresponding to the zero mass-flux at the interface is as follows: Since the proposed method adopts an interface-tracking approach, this equation states that the velocity at interface u Γ and velocity of mesh motionû Γ must be equal to avoid mass transfer across the interface. Substituting velocity from Equation (18), the above equation becomes the following constraint on the interfacial pressure: Equations (17) and (20) are actually associated with boundary conditions at the interface and enforce both the conservation and consistency constraints. It should be noted that while the pressure in Equation (20) is free to hold different values at two sides of the interface, the velocity at the interface can only take one value. This means that it is allowed to use the following equality in the discretization procedure: Therefore, for each interfacial CV, there are three unknown variables, that is, u Γ , p Γf and p Γs , which are constrained to obey Equation (17) and two equations represented by Equation (20). Moreover,û Γ is the mesh velocity at the interface, which is lagged in the iterative computational procedure as explained later.
It is worth mentioning that the concept of pressure discontinuity at the interface of fluid and solid, that is, assuming two values of the pressure associated with two sides of the interface, has been employed by other researchers as well, for example, Reference 37.

NUMERICAL SOLUTION APPROACH
In the current study, the cell-centered finite volume method in conjunction with the PISO algorithm is used to develop a solver for a weakly compressible unified continuum in the context of an ALE framework. Two distinguishing features of the proposed solver merit particular attention. The first one is the flexibility of the ALE, which allows arbitrary mesh motion, from fixed Eulerian to moving Lagrangian, depending on the nature of the problem and the extent of the domain deformation. The second one is the use of a weakly compressible model which is more realistic compared to the assumption of strict incompressibility. 38 In order to apply the FV method, firstly, the computational domains, that is, time and spatial domains, are discretized into a set of time steps and CVs. As shown in Figure 3, each CV can have a general polygonal shape which is constructed by an arbitrary number of control surfaces f . The unknown quantities are defined at the center of CV, that is, node P. The vector d f connects node P to the node of neighbor CV, that is, node N, and n f is the normal vector at the control surface which is pointed outward of the CV.
The discretization procedure is briefly described below: • Implicit Euler scheme 39 for the time discretization • Second order upwind scheme 40 for the advection velocity at CV faces • Approximation of surface normal gradient by the following formula: where d f is the distance between two neighbor nodes; subscripts P and N represent current and neighbor cells, respectively; Surface normal vector n f is decomposed into a vector parallel to d f and its normal direction, that is, d and k, as follows: It should be noted that for the discretization of the gradient terms at either side of the interface face, d f is computed as the distance between the interface face and the adjacent cell in the corresponding region. That implies the pressure gradient and velocity gradient are not continuous quantities across the interface. Also, in this section, a PISO-like algorithm is introduced which is applied to solve the system of coupled equations throughout the entire solution domain. Even though similar dependent variables and coordinate systems are used for both fluid and solid parts of the solution domain, the final form of the governing equations for the fluid and solid are different due to the difference of constitutive relations. In order to employ a unified formulation, a source term is added to the momentum equation as follows: where, Here, pseudo is (2 s + λ s ) Δt where s and λ s are material constants for the solid part. The discretized form of the momentum equation can be written in the following compact form: Where operator H(u) was defined previously. The unified pressure equation is also formulated as follows: Solution of this equation is followed by correcting the velocity and calculation of the flux term: Figure 4 illustrates the iterative numerical solution procedure.
In the present study, OpenFOAM, which allows users to employ a rich group of existing solvers and to develop new computational tools, is used to implement the proposed FSI solver. This software benefits from the object-oriented capability of C++ to present various numerical tasks within modules, libraries, and classes and, consequently, is suitable for both users and developers. Users mostly employ applications, which consist of solvers and utility parts of this software, while various options in these applications are available for users via libraries. For instance, "simpleFoam" is a solver for the solution of steady-state incompressible flow problems and users of this solver can either choose a turbulent model from available libraries or develop a new library for a new turbulence model. Generally, not only turbulence models but also many other models and numerical procedures, for example, mesh motion and boundary conditions can be added as new libraries.
Considering the above explanations, a number of libraries were developed in C++ using existing libraries in Open-FOAM. The purpose of newly developed libraries is to carry out some tasks such as defining pressure and velocity boundary conditions for both traction and interface boundaries and also for adding constituent models such as the Mooney-Rivlin model (Algorithm 1). Algorithm 1. Part of the pseudocode, implemented in the framework of OpenFOAM is as follows A) Read data, for example, mesh, initial conditions, material properties, and so forth. indicates this error at the previous outer iteration. Parameter is set to a value so that the desired convergence level is achieved.
The movement of the fluid mesh is handled by solving the Laplace equation for the mesh velocity, that is,û, then multiplied by Δt to give mesh displacement. The boundary conditions for this equation are zero-velocity for all boundary faces except for those at the fluid-solid interface. For these faces, the velocity at the interface is obtained from the solution of the unified equations and then is applied as the boundary condition of the Laplace equation. For the solid mesh movement, two options are available: ALE and Lagrangian. For the ALE method, mesh velocities at all boundaries are obtained from the solution of the unified equations. The mesh velocity of the internal nodes is obtained by solving a Laplace equation in the solid domain, then, multiplied by Δt to prescribe the motion (displacement) of the internal points at the current time step. In the Lagrangian approach, the motion of the internal and boundary points in the solid region is also evaluated by mapping velocities from cells to points and then multiplying the point velocity by Δt.
It must be noted that since mesh motion for both regions is carried out within the outer loop, both fluid and solid domains are treated as unknowns. This is in contrast to the standard OpenFOAM solver for FSI, icoFsiElasticNonLinUL-Solid, where the domain at the previous time step is used for the discretization of the equations in the solid region (Updated Lagrangian). Therefore, the proposed method for the mesh movement may be called a unified ALE method. However, when we choose to move the solid mesh using the Lagrangian method, it is actually an ALE-DL method, where ALE refers to the fluid part and DL (Deformed Lagrangian) refers to the solid part.
It is worth mentioning that since for the solid part, the deviatoric part of the Cauchy stress is evaluated explicitly, other solid models, for example, the elasto-plastic model, can also be implemented and used in this solver. In the present study, only the hyperelastic model is employed.

VALIDATION TEST CASES
Three test cases are used to investigate the performance of the proposed solver. In order to validate the proposed solver, a challenging benchmark case for internal flows, namely wave propagation inside a flexible, fluid-carrying pipe, is used as the first test case. According to the study carried out by Gerbeau et al., 41 this case differs from other FSI cases found in the area of aeroelasticity since its numerical stability strongly depends on the accurate resolution of interfacial conditions. Therefore, this case seems to be an appropriate validation case for the proposed solver and the innovative formulation of the consistency conditions at the interface. The second case combines two singe-physic problems, namely, flow past a cylinder with a flag behind it and deformation of the cantilever beam due to loads exerted by a fluid. To validate and assess the stability of the solver in various coupling conditions, this case is solved in three sets of parameters for the fluid and solid parts which lead to both stable deformation and self-induced oscillation of the flag (beam) structure. The third case is a 2-D version of the first case; it uses a simpler geometry with structured mesh, yet possesses the main characteristic of the first test case, that is, interaction between internal flow and its surrounding wall. Such simplification for FSI cases in internal flows has been also used in the literature. 41,42 The aim of such simplification in this study, however, is to carry out a comparative study between the proposed method and a partitioned solver in terms of accuracy and robustness using various grid resolutions. The algorithm for the partitioned solver and its convergence criteria has been explained in detail within the literature. 43 The partitioned solver used in this study is a modified version of icoFsiElasticNonLinULSolid solver in OpenFOAM. It is called wcoFsiElasticNonLinULSolid, and has been developed by Tandis and Ashrafizadeh 38 to simulate weakly compressible fluids in flexible tubes. For all cases in this study, the Aitken scheme is applied for the enforcement of the coupling conditions for the partitioned approach. By trial and error, it is found that the acceptable tolerance of the normalized interface residual for this algorithm is 10 −5 for all cases since a residual lower than this value does not affect the final results. This solver can be used to solve a well-defined OpenFOAM case of interaction between weakly compressible fluid and hyperelastic solid with the option of applying various hyperelastic models.

Wave propagation in a fluid inside a flexible pipe
In this section, a transient, three-dimensional (3D) FSI problem, originally proposed by Formaggia, 42 is studied using the proposed method. Being the theoretical model for simulating blood flow in a vessel, this benchmark case has been frequently used in the literature to evaluate various FSI solvers. 38,42,[44][45][46] Therefore, this test case not only is a suitable problem for validating results but also can be used to assess the performance of the solvers in terms of the computational cost. It should be noted that only a limited number of experimental studies for this case are available in the literature. 47,48 In particular, as the flexibility of the elastic wall increases, one hardly finds experimental results regarding the impact of the moving wall on the flow field. The test case, shown in Figure 5, consists of an elastic tube with inner diameter and thickness of 10 and 1 mm, respectively, through which laminar fluid flow occurs due to the pressure pulse at the inlet. The fluid is characterized by kinematic viscosity of 3 mm 2 /s and density equal to 1000 kg/m 3 . To mimic the incompressibility condition, here, we set the bulk modulus to 10 GPa. Mechanical properties of the tube wall, modeled by St-Venant material, are E = 0.3 MPa, = 0.3, and = 1200 kg/m 3 . The tube is fixed at the left and right ends while the outer surface of the wall satisfies the zero-traction constraint. The inner surface of the tube wall is interacting with the flowing fluid. At t = 0, the fluid is at rest and pressure suddenly changes to 1333 Pa. The pressure is then kept at 1333 Pa for 3 ms and finally drops suddenly to zero Pascal at t = 3 ms. The simulation lasts for 0.02 s. Other boundary conditions include zero fluid pressure at the outlet and a moving wall for the shared interface.
Following a methodology similar to the previous studies, a 3D model of the problem is solved first and the results are compared to those reported in the literature as well as the results obtained by a partitioned solver.
For both partitioned and proposed approaches, the solution domain is discretized using 3D cells as shown in Figure 6. Table 1 provides information about the computational grids used in the grid refinement study. Both partitioned and the proposed unified solvers are employed to simulate this case on similar meshes and with the same time step (Δt = 10 −4 s). Figures 7-10 provide pressure contours and deformation of the tube wall for four time-steps, where the deformations are exaggerated by a factor of 15 for the sake of clarity. As seen in the contours of these figures, the motion of the longitude pressure wave inside the fluid is the main cause of the radial deformation of the solid, which behaves as a transverse wave. The interested reader might refer to Reference 38 for further information.
In Figure 11, axial and radial displacements obtained by the unified solver for two grid sizes are compared with the results from the partitioned solver as well as those reported in Reference 46. As shown in this figure, the accuracy of the results on grid level 2 is acceptable. In addition, a comparison of the results obtained by the unified FSI solver with those of both the partitioned solver and reference in Figure 11 reveals that the unified solver is capable of producing results with acceptable accuracy.
It must be noted that in this case, tight convergence for the unified pressure equation in the inner iterations loop leads to instability. In order to avoid that, a limited number of iterations, here 10, is specified by the user. This issue is discussed in more detail in Sections 4.3.

Flow past a 2-D cylinder attached to an elastic flag within a channel
This test case is a benchmark FSI problem originally suggested by Turek and Hron. 49 The schematic of the problem is illustrated in Figure 12. For this test case, we simulate the interaction between an elastic flag attached to a fixed cylinder mounted in a 2-D channel, where a weakly compressible laminar flow passes. Depending on the fluid and solid parameters, the interaction results in either stable deformation or self-induced oscillation of the flag structure with a particular amplitude and frequency. The boundary conditions for the fluid flow include: • Parabolic velocity profile at the left of the channel with a smooth increase in its value with respect to time as below: is velocity profile at the steady state condition, T is the oscillation period for the profile at the inlet, and V is the average inlet velocity.
• No-slip condition at top and bottom wall of the channel, and fluid-solid interface, that is, circle and flag.
• Non-reflective boundary condition 50 at the right side of the channel.
Here, the Reynolds number is defined by Re = 2rV v f , where f denotes the kinematic viscosity of the fluid and r is the radius of the cylinder. Numerical results are available for Re = 10, 100, and 200.
This test case is a challenging problem as it helps evaluate FSI solvers while dealing with some main issues of the fully coupled problems involving large deformation and high value of fluid-solid density ratio which usually result in critical instability. Therefore, performance, for example, robustness and numerical efficiency, of most of the fully coupled FSI solvers 33,51-53 were assessed via this benchmark problem by changing the fluid or solid characteristics. For instance, Hoffman and his colleagues employed a fully coupled FSI solver based on unified continuum modeling 33 to solve this case with the assumption of the Neo-Hookean material model for the incompressible structure, while for the original case, 52 St.Venant-Kirchhoff model was applied to describe constitutive equation for a compressible hyperelastic solid. Hoffman et al. used amplitude and frequency of y-displacement oscillation to compare their results with two references and found them in a reasonable match, that is, for amplitude within less than 1%-2% of the results reported by Hron and Turek 52 and 11% of those reported by Dunne et al., 53 and for frequency within 3% of results reported by both references. Hoffman justified these errors by noting that different constitutive models were applied in the references.
In this study, results given by the proposed algorithm are validated against those reported by Horn et al. 52 at various degrees of coupling. Control point A, shown in Figure 13, is used to measure quantities like x-and y-displacement of the end of the beam structure with respect to time, that is, x(t) and y(t). In addition, forces exerted on the whole structure, that is, lift and drag forces acting on the system of cylinder and beam, are the other quantities of interest which are measured with respect to time. The Geometric dimensions, as presented by Turek and Hron 49 and shown in Table 2, are the same for three different configurations.
As shown in Table 3, the average velocity at the inlet and material properties associated with both fluid and solid are changed to generate various degrees of coupling. It should be noted that for the original case, the bulk modulus K is, in fact, infinity as the fluid is assumed incompressible. Here, we use a large K to mimic the incompressibility condition; however, as mentioned, the flow in this study is not truly incompressible but weakly compressible.
To achieve the mesh-independent solution, three meshes with different resolutions are used. Figure 14 shows the mesh Level 2 and more details are present in Table 4. For the results presented for this case, mesh Level 3 is employed, as this size of mesh was found sufficiently fine for the validation purpose.
Despite being solved in many studies, the available results for this test case do not exactly match our needs for validation purposes. The reason is that in our study, the fluid is modeled as weakly compressible while, to the author's knowledge, fluid was assumed incompressible in all available studies. As a result, we have to modify this test case to be solvable by the proposed algorithm. One of these modifications is associated with the pressure condition at the outlet of the channel. More precisely, since we employ a weakly compressible solver, a consistent boundary condition for pressure must be adopted at the outlet. The purpose of such a boundary condition, which is different from zero gradient pressure at the outlet for the original problem, is to have a boundary condition that does not reflect pressure waves from the outlet. This type of boundary condition, called "pressureTransmissive" in OpenFOAM, is explained in Reference 50. Figures 15-17 present velocity contours and streamlines inside fluid regions as well as deformation of the flag for three configurations, that is, FSI1, FSI2, and FSI3 respectively. As seen in these figures, the FSI1 configuration leads to a stable solution with the smallest deformation of the flag. In addition, while the deformation of the FSI2 configuration is the largest among the others, wakes behind the cylinder in the FSI3 configuration are larger and move along the length of the flag.
For the FSI1 configuration, as reported by Horn et al., 52 the solution leads to a stable state for the solid. The final value for the displacement of control point A in the x-and y-directions and also lift and drag force are presented in Table 5. Figures 18 and 19 show the transient motion of the control point A in both x-and y-directions as well as variation of force components in terms of time for the FSI2 and FSI3 configurations. As clear, the configurations of FSI2 and FSI3 lead to oscillatory behavior. Comparing the shape of curves with those reported by Hron and Turek 52 shows the qualitative agreement between these two cases. Tables 6 and 7 give a quantitative comparison between results obtained by the unified algorithm and reference. In order to compare the results in more detail, these tables present the frequency, mean, and amplitude of the displacement and force oscillations. As seen, the proposed solver produces results with acceptable accuracy when compared to those presented by Hron and Turek. 52 The slight difference between the results given by the proposed solver and those reported

Flow between 2-D compliant plates
In the first two test cases, the proposed algorithm was validated and its capability of dealing with strongly coupled problems was proven. The current test case is actually aimed at comparing the partitioned and proposed unified solvers in more detail on a simplified geometry. This particular case, which can be regarded as a Cartesian version of the previous test case, is chosen due to its simplicity (generating rectangular mesh on the computational domain) and also its low computational cost. Figure 20 demonstrates a schematic of this case: Here, the boundary condition of the inlet and outlet and also dimensions of the geometry are the same as in the first test. In order to provide different degrees of coupling, three types of mechanical properties for solid walls are adopted. Table 8 gives details of this test case for three configurations. These three configurations have similar geometrical dimensions and fluid parameters and they only differ due to the difference in mechanical properties of the solid. The change in the solid properties provides different degrees of coupling. The degree of coupling is, in fact, proportional to the fluid-solid density ratio and inversely proportional to the modulus of elasticity. Therefore, the test cases here can be arranged in terms of the coupling intensity, that is, case3 > case2 > case1.
The total time of simulation is set to 0.01 s for all cases. In order to provide a fair comparison, various levels of discretization for time and space are applied. Table 9 presents details of six levels of mesh which are used with four time-steps, Δt = 100, 50, 25, and 10 μs, to carry out the computations.
For both partitioned and unified solvers, similar discretization schemes are used for both spatial and temporal terms, that is, Euler implicit for temporal terms, "linearUpwind" for the convective term, and "Green-Gauss" scheme for the explicit evaluation of gradient inside cells. Other required interpolations are basically geometrical in nature. A tight accuracy criterion is employed for the partitioned solver in terms of both residuals in individual regions as well as the satisfaction of coupling conditions. The simulations were carried out on a PC with Intel(R) Core(TM) i7-6800K CPU @ 3.40 GHz and 32 GB memory. Finally, horizontal and vertical displacement of point A (L/2, H + t), shown in Figure 20, in terms of time is used as desired output for the sake of comparison. Figure 21 shows the contour of the fluid pressure as TA B L E 8 Specification of different coupling intensities for the case of flow between 2-D compliant plates.

Model
Quantity Value well as deformation of the walls for several time instances in the coupling 1 case. It should be noted that the deformations shown in these contours are exaggerated by a factor of 5. In Figure 22, transient variations of vertical and horizontal displacements of point A, obtained by both solvers, are presented. It is clear that the final results obtained from partitioned and unified solvers are in good agreement. In addition, as shown in Figures 23 and 24, the results obtained by the unified solver do not significantly change as the mesh gets finer. In other words, the dependency of the results of the proposed algorithm on the mesh size vanishes at mesh level 2. On the other hand, mesh convergence for the partitioned solver occurs at mesh level 4. Therefore, the computational results from the unified solver seem to be more accurate on a given mesh. Figure 25 shows the deformation of the wall for coupling 2 at four time instances. In order to compare with the coupling 1 case, the deformation for this case is also exaggerated by a factor of 5. Comparing Figures 25 and  21, it is obvious that the deformation in the coupling 2 case is larger than in the case of coupling 1. This is due to the lower modulus of elasticity in the coupling 2 case. Another difference is associated with the lower speed of the vertical displacement wave in the case of coupling 2. This is also due to the lower shear modulus ( s ) as a result of the lower Young's modulus, since the shear modulus of the solid has a positive effect on the traveling speed of this wave. Figures 26-32 compare the results obtained from the partitioned and unified solvers for coupling level 2. According to Figure 26, the results of both solvers are close to each other when the finest sizes are used for the time step and the mesh. However, the partitioned solver becomes unstable for Δt = 10 −4 s when a mesh finer than level 3 is used. Also, compared to the coupling level 1, the partitioned solver needs smaller time steps with more solution time to converge. This behavior is actually due to the explicit exchange of the information in the partitioned approach. In comparison, the proposed unified solver converges for all time steps without a major change in the solution time. Moreover, the proposed solver converges on a coarser mesh for all time steps.
The degree of coupling in coupling level 3 is the strongest among all the cases in this section. The partition solver becomes unstable and diverges for all meshes, L1-L5, and all time steps, 1 × 10 −4 to 1 × 10 −5 . However, the unified solver is capable of solving the problem on all meshes and all time steps. Figure 33 shows the deformation of the wall at different time steps for coupling level 3. Note that the amount of the deformation is not much different from the previous case, coupling level 2, due to the fact that the solid in both cases possesses the same elasticity and the only difference is associated with the solid's density difference. However, since the density of the solid is an order of magnitude lower than in the previous case, the degree of coupling for this case is stronger. This explains the failure of the partitioned solver for this very strongly coupled test case.
As shown in Figures 34-37, the results of the proposed solver for all time steps are independent of the size of the mesh for this coupling level 3 case. In addition, Figure 38 compares the results for four time-steps on mesh level 4. Note that the time step of 2.5 × 10 −5 s is small enough to obtain the desired accuracy in the calculation of the vertical displacement.
As the results for the second case, flow between 2-D compliant plates, at a high degree of coupling suggests the proposed unified solver outperforms in terms of stability when it is compared to the performance of the conventional partitioned approach. In the authors' view, the underlying reason for such superiority is the simultaneous solution of the equations in both domains as well as innovative treatment of the consistency conditions which have become feasible via unification of the formulations and unknown variables throughout the entire computational domain, including interface. [Colour figure can be viewed at wileyonlinelibrary.com] As a result, the proposed solver is a monolithic type which is capable of solving strongly-coupled cases without applying relaxation schemes.
It has to be mentioned that a fair comparison between partitioned and unified solvers in terms of the CPU time is difficult to make because the computer codes written for these two solvers are not optimized similarly. However, in order to provide more insights into the computational time for both solvers, Table 10 presents the number of outer iterations as well as total computational time for different coupling conditions with various discretization sizes.
As clear from this table, compared to the unified algorithm, the partitioned solver demands less computational time for coupling 1; however, the performance of this solver degrades for coupling 2 and 3. On the other hand, for the unified solver, an opposite behavior is observed, that is, the unified solver is more robust for strongly coupled conditions. In other words, for the weakly coupled case, the partitioned approach seems to be more efficient; as the degree of coupling increases, the convergence rate of the partitioned solver dramatically decreases while that of the unified solver becomes relatively better. Additionally, the results for the unified algorithm suggest that applying smaller time steps reduces the overall solution time for coupling 1 and has an opposite effect on coupling 3. For coupling 2, changing the time step does not lead to a major change in the solution time. It should be noted that for coupling 2 and 3, the unified pressure equation is solved using 10 and 5 iterations, respectively; using a tighter tolerance for the solution of this equation could lead to divergence, as encountered in the previous test case. In the author's view, seeking an efficient convergence criterion for the inner PISO loop as well as the solution of the pressure equation might improve the convergence behavior. For the partitioned solver, on the other hand, using a very small relaxation factor, for example, 10 −5 s. could not avoid instability for the coupling degree 3 and the solver encounters divergence after a few outer iterations within the first time step.

Stability of the proposed solver
As mentioned, the proposed unified solver in the current form is prone to the instability issue associated with the exact solution of the unified pressure equation. The number of iterations for the pressure solver as well as that for the PISO loop not only could cause instability but also could affect the convergence rate and solution time. In this section, the underlying reason and proposed remedies are discussed.
In the authors' view, the observed instability is due to the explicit treatment of the coupling between the unified pressure equation and the interfacial boundary conditions. More precisely, as mentioned in Sections 2.4, the pressure and velocity at the interface are coupled together as well as to the pressure and velocity values at the adjacent cells. As a result, an appropriate solution strategy is needed to avoid instability. In this study, to avoid complexity in the implementation, TA B L E 10 Solution details of the partitioned and unified solvers for three coupling conditions for the case of flow between 2-D compliant plates and using various discretization sizes. the coupling between equations at internal cells and those deduced from the interface conditions at imaginary cells are not solved in a single equation system but simply as two separate equation systems whose coupled terms are lagged and updated within the outer loop. Therefore, using an internal PISO loop with a very tight tolerance might lead to instability. One solution to this issue is to simply solve internal and interface pressure equations in a single system; this approach of applying boundary conditions as additional equations inside the overall equation system is also adopted for fluid only 54 and solid only problems. 55 The explicit treatment of the coupling between pressure and velocity via the PISO algorithm also escalates the issue of instability. Therefore, as a second remedy, an alternative to the PISO algorithm can be applied in order to couple unified pressure and unified velocity equations. In fact, despite the popularity of the PISO due to its simple implementation, the coupled algorithms have shown computational benefits 54,56 and can be used in order to mitigate instability and convergence rate issues. That means the simultaneous solution of all equations together in a fully coupled manner, that is, unified pressure, unified momentum and interfacial equations solver. Therefore, it seems that as fluid-solid coupling becomes stronger, developing a fully implicit algorithm for the solution of the unified pressure-velocity formulation becomes more beneficial from the stability point of view.

CONCLUDING REMARKS
A unified formulation and solution algorithm was proposed for strongly-coupled fluid-solid interaction problems. Solution of FSI problems in detail. In the proposed solver, pressure and velocity equations were developed at the interface to conservatively enforce the consistency constraints. Two benchmark test cases with different coupling intensities were numerically simulated using an already developed partitioned solver in OpenFOAM and the proposed unified solver. A third test case was also introduced to examine and compare the stability of both solvers. The comparison revealed that: • The proposed algorithm accurately models strongly coupled fluid-structure interaction problems.
• The proposed solver provides more accurate results on a coarser mesh as compared to the partitioned solver.
• The partitioned solver has a better convergence speed for weakly coupled problems.
• The unified solver solves strongly coupled problems with different mesh resolutions and time steps for which the partitioned solver fails to converge.
This study showed that pressure-velocity formulation for fluid-solid interaction problems can be used as an alternative to the traditional partitioned approach particularly for relatively largely deformed computational domains. In particular, the proposed formulation allows for the simultaneous solution of the equations for fluid and solid regions as well as those at the interface, that is, a monolithic solver where information is implicitly exchanged between fluid and solid regions, and therefore potentially improves the solution stability for the strongly-coupled problems. However, in the current state of development, the convergence speed of the unified solver depends on the fine-tuning of the number of inner iterations in the PISO-type loop of the solver. Fully implicit treatment, that is, simultaneous solution of the pressure and velocity in fluid and solid regions as well as at the traction and interface boundaries, is suggested as the next step for this study and is expected to improve the stability and convergence rate.
Another suggestion for further development of the proposed method is to include elasto-plastic solid models in the solver. Since no major assumption is made for the constitutive relation in the unified momentum equation (stress is calculated explicitly), implementation of the majority of constitutive relations is possible.
Finally, it is worth mentioning that the time step is fixed throughout the solution process for cases studied in this paper. However, since fluid and solid regions have, in general, very different physical properties, appropriate time steps can be quite different for these regions. Therefore, methods such as adaptive time stepping 57 or locally varying time stepping 58 schemes can potentially improve the convergence rate and stability.

FUTURE WORK
The main objective of this article is to develop a unified pressure-velocity formulation for interaction between weakly compressible flow and a hyperelastic solid in the context of the finite volume method. The test cases in this paper demonstrate the new formulation combined with the PISO algorithm can accurately solve the FSI problems with strong coupling conditions, which proves the validity of the proposed formulation as well as its superiority over the classical partitioned solver. In the authors' view, adopting a coupled algorithm (solving both pressure and velocity equations inside all sub-domains and at the interface) and using proper preconditioning of the overall equation system can enhance the stability of the unified solver and, hence, is recommended for future studies.

ACKNOWLEDGMENT
Open access funding provided by IReL.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.