A Unified Arbitrary Lagrangian–Eulerian Model for Fluid–Structure Interaction Problems Involving Flows in Flexible Channels

In this work a finite element-based model for analyzing incompressible flows in flexible channels is presented. The model treats the fluid–solid interaction problem in a monolithic way, where the governing equations for both sub-domains are solved on a single moving grid taking advantage of an arbitrary Lagrangian/Eulerian framework (ALE). The unified implementation of the governing equations for both sub-domains is developed, where these are distinguished only in terms of the mesh-moving strategy and the constitutive equation coefficients. The unified formulation is derived considering a Newtonian incompressible fluid and a hypoelastic solid. Hypoelastic constitutive law is based on the strain rate and thus naturally facilitates employing velocity as a kinematic variable in the solid. Unifying the form of the governing equations and defining a semi-Lagrangian interface mesh-motion algorithm, one obtains the coupled problem formulated in terms of a unique kinematic variable. Resulting monolithic system is characterized by reduced variable heterogeneity resembling that of a single-media problem. The model used in conjunction with algebraic multigrid linear solver exhibits attractive convergence rates. The model is tested using a 2D and a 3D example.

decades, it remains an area of active research to-date. Particularly important issue remains the balance between the robustness of the model and the computational efficiency, since for many practical applications lengthy simulations involving days or even weeks of computational time are unacceptable. In the present work we shall concentrate our attention on FSI problems involving flows in flexible channels, aiming at biomedical applications.
Two types of approaches to FSI modeling exist, namely the partitioned and the monolithic ones [1]. Partitioned approaches rely on solving the solid and the fluid problems using different solvers that are coupled via transferring data at a common interface. In monolithic approaches coupling conditions at interface are implicit to the solution procedure as they are included together with the governing equations of the sub-domains in a single discrete system. In spite of their robustness monolithic approaches generally result in poor conditioning of the linearized system describing the coupled problem due to different scaling of the variables entering the monolithic system. Solids are usually described in terms of displacements, incompressible fluid formulations are typically velocity-pressure-based and interface compatibility is enforced by Lagrange multipliers technique. This leads to a very large difference in the magnitude of coefficients entering different blocks of the monolithic generalized stiffness matrix (left-hand-side of the equation).
One way of tackling this challenge consists in developing or applying sophisticated linear solvers and/or pre-conditioners capable of dealing with given large and poorly conditioned systems. This approach was followed by various researchers in the field of monolithic FSI [2][3][4][5][6][7]. In this case computational efficiency of the model is determined predominantly by the selected linear solver.
A different approach to the efficiency problem relies on manipulating the monolithic system prior to passing it to the linear solver. Reducing variable heterogeneity and facilitating the FSI coupling can be achieved, among other techniques, by the so-called "unified approaches" proposed originally in [8]. Unified approaches [9][10][11][12] describe the entire FSI problem on a single mesh and define a unique kinematic degree of freedom for the entire domain, including the nodes where fluid and the solid are in contact . All the mentioned models rely on Lagrangian description of the entire problem and are advantageous for analyzing FSI problems involving free surface flows.
The use of unified Lagrangian approaches in 3D is rather limited due to their high computational cost associated to continuous re-meshing and time step size restrictions [13] Moreover, for modeling flow in flexible channels, Lagrangian approaches are not optimal as such problems involve flow continuously passing through a fixed control volume. Application of inlet and outlet boundary conditions in thus inconvenient as it requires injecting nodes at the inlet and erasing them at the outlet [14], which may cause additional computational difficulties. Clearly, for such problems, using non-Lagrangian frameworks is more convenient for describing not only the inlets and outlets (as proposed e.g. in [15]), but also the major part of the fluid domain.
In 2000s Turek and Hron [2,16] advocated the use of the Arbitrary Lagrangian/Eulerian (ALE) framework for FSI problems involving channel flows, particularly for biological systems (flow in blood vessels). They took advantage of using using different mesh-motion strategies in different sub-domains. Similar ideas have been later used for developing models for various flow problems in flexible channels in [17][18][19][20][21]. The above-mentioned setting allows treating the solid in a Lagrangian way, while modeling the major portion of the fluid in an ALE fashion, keeping the mesh fixed at the inlet/outlet boundaries.
In the present work we inherit the idea of using different mesh motion strategies for different parts of the domain and combine it with a unified approach (in the sense of both using single mesh and single kinematic variable) to the coupled problem. While commonly one uses an interface equation for ensuring compatibility between the displacement-based solid and the velocity-based fluid, we choose a single kinematic variable (velocity) for both sub-domains, thus naturally ensuring the kinematic and continuity of the velocity field across the interface. Unification also reduces the variable heterogeneity, which is beneficial for the computational efficiency of the model.
We strive to obtain a method that is characterized by good computational efficiency, due to reduced variable heterogeneity and reduced size of the coupled system, facilitating improved linear solver convergence. Thus, our approach may be viewed as an alternative to the monolithic FSI approaches where the efficiency improvement is obtained exclusively at the level of the linear solver. The developed approach has an additional benefit from the implementation point of view: a single code encompassing both fluid and solid elements (by varying materials coefficients) allows implementing the monolithic solver in a very similar way to that of a single-material problem (fluid-only or solid-only), using a standard finite element assembly procedure and other standard routines.
The paper is organized as follows. In Sect. 2 the rationale of the method is explained. Section 3 is devoted to the governing equations of the fluid-structure interaction system. Unified constitutive equation is presented and a unified form of the discrete monolithic FSI system is derived. Corresponding solution algorithm is presented with a certain emphasis given to the decoupled mesh smoothing step. In Sect. 4 numerical examples are solved. The model is validated by means of comparison with several references in 2D and 3D. In Sect. 5 concluding remarks are drawn.

Rationale of the Model
The fluid and the solid are discretized here using a unique mesh and are integrated in time as a single entity. In order to avoid the necessity of introducing the interface governing equations (ensuring continuity of the kinematic variables and normal stresses across the fluid-solid boundary), same kinematic variable for the fluid and the solid domains. During the monolithic step we propose not to solve any additional equation for defining the fluid mesh motion. These two features distinguish the present approach from [2] and other similar ALE models mentioned in the introduction, where difference in the kinematic variables between the solid and the fluid domain as well as the presence of ALE mesh velocity variable oblige introducing additional interface equations.
In order to illustrate the proposed mesh-moving strategy of the model, let us consider an example involving flow in a flexible pipe. The undeformed configuration is depicted in Fig. 1. The domain consists of a flexible solid (gray outer walls) and the fluid inside. A constant flux is prescribed on the left (inlet), the solid is fixed at its left and right extremes. The fluid flow deforms the solid, resulting in the configuration depicted in Fig. 1. One can see that the domain boundaries at the inlet and the outlet remain unchanged, while elsewhere the solid deforms in the direction predominantly orthogonal to the direction of the flow. Thus the fluid-solid interface can be tracked by a mobile mesh simply stretching the corresponding fluid elements. Once the monolithic FSI system (where mesh motion is completely defined by the deformation of the solid) is solved, a mesh smoothing step can be performed leading to an improved quality of the fluid mesh. This is indicated in Fig. 1.
The following combination of mesh-moving strategies is used when solving the monolithic FSI system: • Solid sub-domain: Lagrangian fashion (obtained in ALE model by setting the mesh velocity to be equal to the material velocity). Gray in Fig. 1. • Fluid sub-domain, except for those fluid elements that are encountered in contact with the solid: Eulerian fashion (fixed grid). White in Fig. 1. • Interface fluid elements (fluid elements in contact with the solid): the nodes shared with the solid move according to their convective velocity ("Lagrangian nodes"), while the rest of the nodes are maintained fixed ("Eulerian nodes"). Dashed in Fig. 1.
It is worth emphasizing that in the above-described basic mesh moving strategy employed at the monolithic solution step (Fig. 1b), there exist either "Eulerian" or "Lagrangian" nodes, meaning that the mesh velocity does not constitute an additional unknown.
During the separate mesh smoothing stage (Fig. 1c), position of the fluid nodes (nodes of "white elements") is updated by solving smoothing equation.
In the following the governing equations of the FSI problem are presented. As already mentioned, in order to facilitate monolithic coupling without additional interface equations the solid model is written adopting velocity as the kinematic variable so as to match the degrees of freedom of the fluid. Adopting hypoelastic model for the solid allows obtaining a unified FSI model where the fluid and the solid differ only in terms of the constitutive relation.

Conservation Equations
Let be a bounded domain containing a viscous incompressible fluid and a flexible solid. The evolution of the velocity v = v(x, t) and the pressure p = p(x, t) (the latter being defined only in the fluid sub-domain) is governed by momentum and mass conservation: where σ is the Cauchy stress, and, ρ and v are the density and the velocity vector, respectively and b is the body force. The | x symbol stands for the derivative with respect to the room-fixed coordinate system. The ALE differential form of the momentum equation is readily obtained from Eq. (1) replacing in the convective term the material velocity v with the convective velocity c = v − v m , where v m is the velocity of the moving mesh [22]. This results in where χ is the coordinate associated to the ALE reference system. Note that the right-hand side of the above equations is written in classical Eulerian form, while the arbitrary motion of the computational mesh is only reflected in the left-hand side [22,23].
Equations (2)  where v pr is the prescribed velocity, n is the outer unit normal to N , σ pr n is the prescribed traction vector. On the internal interfaces I , the coupling conditions are with n being the unit normal to int , and [[φ]] represents the jump of a quantity φ across the interface. Note that adopting a monolithic formulation in a single variable the interface conditions Eqs. (6), (7) are satisfied automatically.
The momentum Eq. (3) is independent of the material (i.e. it holds for both the fluid and the solid sub-domains), so we need to incorporate the information of the specific material to be modeled. This is done by means of specifying the material's constitutive equation, which links the stress tensor to the deformation. In the context of the present paper we treat the problems involving two materials: incompressible Newtonian fluid and isotropic hypoelastic solids. Accordingly, in the next section we discuss their respective constitutive equations.

Constitutive equation for incompressible Newtonian fluids
The constitutive equation of a Newtonian fluid is a linear relationship between the Cauchy stress σ and the deformation rate tensor D: where μ f and λ f are the Lamé coefficients and p is the thermodynamic (or hydrostatic) pressure of the fluid. By definition D = 1 2 (L + L T ) is the symmetric part of the velocity gradient tensor L = ∇v. Note that tr(D) = ∇ · v = ε v = 0 for incompressible fluids.

Constitutive equation for hypoelastic solids
An isotropic hypoelastic solid is characterized by a constitutive equation where a linear isotropic relationship between the stress rate and the deformation rate are defined as [24]: where τ is the Truesdell rate. This is evaluated as where τ = J σ andṠ are the Kirchhoff stress tensor and the material time-derivative of the second Piola-Kirchhoff stress. The latter is defined as S = F −1 · P = J F −1 · σ · F −T . F is the deformation gradient and J is its Jacobian.

Time-integrated constitutive equations
An implicit time discretization of the fluid's constitutive Eq. (8) can written as: where superindex n + 1 indicates the current time step. Time-discretization of the constitutive equation of an hypoelastic solid An implicit integration of Eq. (9) can be written as Replacing τ (n+1) by Eq. (10) into Eq. (9), the constitutive equation reads Equation (14) is obtained by applying the generalized midpoint rule forṠ (n+1) and replacing S by its definition.
The above equation introducing the so-called incremental deformation f (n+1) defined as f (n+1) · F (n) = F (n+1) can be written as follows which provides the final form of the time-integrated constitutive equation for the solid.
As mentioned in the Introduction, in the present work the solid will be described in the updated Lagrangian reference frame (as a limiting case of ALE).
In what follows the governing equations will be discretized. Depending on the reference configuration chosen for solving the discrete governing equations two options exist within the Updated Lagrangian framework. First one relies on choosing the known configuration corresponding to the time t n . In this case one obtains the unknowns solving the governing equations on a "frozen" domain of t n . Thus, the unknown and the reference configurations are linked via the deformation gradient (reflecting deformation undergone by the continuum during one time step, i.e. from t n to t n+1 ). The non-linear iterations in this case involves re-computation of the deformation gradient F and its functions. Another option, followed in the present work is to choose the unknown end-of-step (t n+1 ) configuration as the reference configuration, iterative approximating it.

Remark 1
Note that for the finite element implementation choosing current configuration as a reference implies the necessity to iteratively recompute the finite element matrices which change at each non-linear iteration due to the fact that the end-of-step domain configuration is not known apriori. Instead, the current-reference configuration is iterative updated until convergence is achieved, providing the true end-of-step configuration.

Remark 2
A similar idea of a unified constitutive equation has been proposed in [25] considering hyperelastic constitutive equation in the context of a fully Eulerian FSI method.
Taking the above into account and re-defining the Lamé constants asμ

A unified constitutive formulation for a Newtonian fluid and a hypoelastic solid
Now we make use of the similarity between the constitutive equation of the hypoelastic solid (16) and the constitutive equation of a Newtonian fluid (8) so as to write a unique constitutive equation for both materials following the idea presented in [26]. The difference between the two constitutive equations is restricted to the fact that the Lame coefficients of the hypoelastic solid are not constant. Consequently, we can write the following general expressions for both materials: where the terms are defined as shown in Table 1.

Governing Equations at Discrete Level
Next, the steps to be performed in order to obtain the discrete version of the problem in time and space are presented. For the sake of simplicity time integration is illustrated using the Backward Euler scheme. Nevertheless, in the implementation of the present method a version of Newmark-Bossak scheme was used [27]. Spatial discretization was done using mixed linear velocity-pressure finite elements. Note that in order to simplify the notation, the superindex n+1 will be neglected. Multiplying Eq. (3) and the natural boundary condition by a linear finite element test function N (note that N is the vector of shape functions) and integrating over the domain the weak form of the governing equations can be written as Applying the divergence theorem to the stress terms allows us to write the previous equation as Introducing the unified constitutive Eq. (17) into momentum conservation Eq. (20) gives Considering linear finite element linear approximations for the velocity and the pressure unknowns p (x) = N T (x)p and v i (x) = N T (x)v i wherep andv are the nodal values, leads to the following discrete form of the governing equations (note that for the sake of clarity of presentation an overbar indicating the nodal quantities will be omitted in the text below, i.e. instead ofp andv we will write p and v) The discrete governing system defined by Eq. (22) holds for both the fluid and the solid, with γ = 0 in s and γ = 1 in f .
The matrices and vectors presented above are defined as The discrete system Eq. (22) for the fluid is stabilized using the algebraic sub-grid scale method [28]. For the sake of simplicity, stabilization terms (S K , S G , S D , S p ) are denoted here symbolically. Their definition can be consulted in [29] where the stabilization used here is described in detail.

Remark 3
Note moreover that unifying the formulation prior to discretization has an additional benefit in comparison with the approaches that obtain a single kinematic variable posteriorly. By "posterior unification" we mean replacing displacements by velocity at the level of the discrete system. This latter option requires linearizing the stress tensor expressed in terms of displacements with respect to velocity (see e.g. [30], p. 931), which is not trivial. Even when expressing acceleration in terms of velocity in a displacement-based solid, one must take special care to ensure that the used relation is compatible with the employed time integration scheme. In case of using a solid formulation expressed in terms of velocity, these issues do not arise.

Solution Algorithm
Given the solution at t n , v n and p n , at the known configuration X , the procedure to find the solution v and p at time t n+1 using the proposed unified model is summarized in Algorithm 1. For improving the mesh quality after each monolithic step, a mesh smoothing algorithm has been implemented. Among the different options available [31][32][33][34], the Laplacian smoothing technique has been chosen due to its simplicity and efficiency [35][36][37]. The reader can find a comprehensive introduction to mesh-smoothing operators in [38] and detailed review in [33,39].
In the context of the present model, this approach consist in solving the Laplace equation for the increment of the nodal position after solving the conservation equations [Eq. (22)]. Introducing the mesh displacement variable d m,L = x L − x (where x L is the unknown improved position at t n+1 and x is the position obtained by solving Eq. (22)) the Laplacian smoothing equation can be written as: with the following boundary conditions: where sub-index L stands for the position obtained after applying the "Laplacian-smoothing".

Remark 4
Solving the mesh-update equation in a separate step has several benefits. First of all, including the mesh update equation inside the monolithic system requires introducing an additional degree of freedom (mesh displacement/velocity), resulting in an increased monolithic system size and variable heterogeneity. Moreover, in this case the left-hand-side of the monolithic system is more difficult to construct, as the linearization in this case must also include the derivatives of the governing equations with respect to these additional variables. Consequently, these linearization terms introduce new coupling blocks in the left-hand-side matrix. All the issues mentioned-above contribute to an increase of the computational cost of the solution of monolithic system. Some of these aspects are discussed in detail on p. 218 in [40].
Having said this, inclusion of the mesh-motion inside the monolithic system (sometimes called the "direct-coupling") may be beneficial from the point of view of capability to handle larger domain distortions. However, from the point of view of computational efficiency it increases the computational cost. Since our method aims at solving flows in flexible channels, the expected domain deformations are moderate. This justifies our selection of decoupling the mesh-update step.

Remark 5
In general, the "price to be paid" for the above-explained mesh motion decoupling is the necessity of estimating the "safe" time step to avoid excessive contraction of the fluid elements adjacent to the solid during the Lagrangian motion step. For such cases in the present work an estimator (that ensures that an interface node does not surpass a distance larger than the characteristic size of the corresponding adjacent fluid element) was implemented following the scheme described in [13]. It can be simplified to the following expression: where h is the characteristic element size and v norm int is the velocity of the interface nodes in the direction normal to the interface). We note, however, that in the numerical examples presented in the following section no time step adjustment was performed as the time steps of the benchmarks were sufficiently small so as not to lead to element degradation for the faced interface velocities.

Numerical Examples
The present model was implemented by the authors in the ULF Application branch of Kratos MultiPhysics code, an academic Open Source software [41,42]. For the solution of the linear system two solvers are used: the algebraic multigrid solver (AMGCL library [43]) and a biconjugate gradient solver (BiCG).
In the following, three numerical examples are solved. First one is a 2D FSI benchmark dealing with a flow in a cavity with a flexible bottom. This test is often used for the validation of ALE models. We also use this example for getting an insight of the computational efficiency of the model. Second example is another commonly used benchmark, a 3D test case dealing with a pressure propagation in a flexible pipe (its more complex version is solved in [44]).

Fluid-Structure Interaction in Driven Cavity with Flexible Bottom
To test the method and to get the insight of its computational efficiency, a two-dimensional lid driven cavity with flexible bottom is simulated. This example (with slight variations in geometry/boundary conditions) was previously used for testing partitioned and monolithic FSI by various authors ( [19,[45][46][47]    Vertical displacement at the middle of the flexible bottom was recorded (x = 0.5 m, y = 0.002 m in undeformed configuration). The comparison of the results obtained here using different meshes specified above are shown in Fig. 5a. One can see that the three results are nearly identical.
Next, we compare present results with the results found in literature: the partitioned strongly coupled FSI solution of Valdes [46]. In the reference the example was solved on a mesh consisting of 32 × 32 quad elements in the fluid domain. Figure 5b shows the comparison. The average displacement at the periodic state is estimated in [46] as approximately 0.235 m, while present work estimates it as 0.25 m, which corresponds to a discrepancy of ∼5 %. The amplitude is estimated by both models as approximately 0.032 m. Present results are also found to agree well with the very recent work [47] (partitioned scheme, solid modeled using membrane theory).
Computational efficiency Cavity example was used for estimating the computational efficiency of the monolithic ALE solver in [19], where the number of linear solver iterations was recorded for different mesh resolutions. In the afore-mentioned reference the ALE monolithic system included the governing equations for the fluid and the solid as well as the interface equations. The kinematic variables of the solid and the fluid are the displacement and the  [19]. All shown results were obtained using algebraic multigrid linear solver velocity, respectively. For solving the linearized system the authors of [19] employed algebraic multigrid (AM) solver. The Krylov space dimension was set to 50, GMRES solvers were used in conjunction with the ILU(0) pre-conditioner for each field. Convergence tolerance was set to ||r|| ||r 0 || < 10 −5 (where r is the residual). These settings were reproduced in the AM solver used in the present work. Additionally, we tested the model in conjunction with a simpler solver, namely the stabilized bi-conjugate gradient (BiCG) solver equipped with a diagonal pre-conditioner.
The non-linear solver stopped as soon as the velocity and the pressure fields in L 2 norms were below 10 −8 . Figure 6 shows the evolution of iterations number over time obtained using present method in conjuction with an AM linear solver. One can see that for both Mesh3 and Mesh3 our approach converges 2-3 times faster than that of [19]. It is particularly pronounced in case of Mesh3, where present method takes 200-250 iterations to converge, while the non-unified ALE model of [19] required around 800 iterations. The results of the tests performed using BiCG and AM solvers are summarized in Table 2. One can see that when using the unified model in conjunction with a BiCG solver, one obtains similar convergence rates to the ones obtained by a non-unified model with an AM solver. In fact, the number of linear iterations appear to grow as the Courant number of the problem grows (considering the mesh sizes used and the maximum velocity of the flow (2 m/s) it is simple to compute Courant number C for each case analyzed) . Independent of the mesh size, all the solvers show acceptable convergence rates for small Courant number (C < 1) considerably growing for C > 2. However, the present unified model exhibits a much smaller growth as C increases. It is well-known that for small Courant number the influence of the mass matrix term (which is a diagonal matrix) becomes stronger (in comparison with the convective and viscous matrices), which improves the conditioning of the linear system and facilitates convergence of iterative solvers.
Another interesting observation is that the convergence of the linear solver is clearly affected by the mesh distortion. It reaches maximum at the moment of maximum channel contraction, which results in elements with impoverished aspect ratio.

A Pressure Pulse Propagating in a Flexible Pipe
This benchmark deals with a propagation of a pressure wave in a viscous incompressible fluid enclosed into a flexible 3D pipe. The pipe geometry is as follows: length of 0.05 m, inner diameter of 0.01 m and wall thickness of 0.001 m. The geometry and material properties are taken from [20,21]. The material properties are summarized below: The pipe is fixed at both ends. At the inlet an external pressure p ext of 1333 Pa is prescribed during a time span from 0 till 0.003 s. At the outlet zero pressure is prescribed. The overall simulation time is set to 0.02 s. The model was discretized using ∼300.000 tetrahedral elements and the time step was set to 0.0001 s. Linear solver convergence tolerance were set identical to those of the previous test example (Sect. 4.1). Figure 7 shows the pressure wave propagation and the corresponding deformation of the solid in the longitudinal cross-section at t = 0.004 and t = 0.1 s. Figure 7 also displays the corresponding results presented in [20]. Figure 8 displays the pressure distribution and the domain configuration at t = 0.02 s. The domain displacement is magnified by a factor of 10 for clarity. The results are also nearly identical to those of [21]). For the sake of brevity the snapshots of the latter reference are not reproduced here. Figure 9 shows the comparison of the the present simulation results with the ones of [20,21]. Axial and radial displacements were measured at the inner tube wall at half the length of the pipe (x = 0.025 m). In order to perform a fair comparison the mesh was constructed ensuring existence of the nodes exactly at x = 0.025 m.

Results
The linear solver iterations history is shown in Fig. 10. One can see that for the time step originally employed (dt = 0.0001 s) BiCG solver converged in less than 100 iterations per step, while AM solver takes around 40 iterations. We note that the Courant number in this case approximately equals C = 0.1.Number of non-linear solver iterations varied between 2 and 3.
To test the solver at larger Courant number, the time step was increased to 0.002 (corresponding to C = 2). In this case BiCG solver iterations increase to 800 iterations per step in average. However, when using the algebraic multigrid solver the number of iterations grows much less, reaching 200 iterations per step only. This behavior is very similar to the convergence behavior observed in the previous example.

Summary and Conclusions
In the present paper a unified monolithic method was proposed for modeling flows in channels with flexible walls with small-to-moderate deformations. Different sub-domains of the coupled problem were treated using different mesh motion strategies according to the setting of the problem at hand within a unique framework. The unification we proposed had three main ingredients: (i) a single deforming boundary-conforming mesh for the entire FSI problem ii) a single kinematic variable in the model already at the continuum level iii) unified implementation (solid and fluid elements are implemented in a single generalized formulation, with the corresponding contributions distinguished by variable coefficients). This allowed reducing the variable heterogeneity in the monolithic system and permitted solving the coupled problems without defining any interface equation. Moreover, our approach allows for a simple implementation where a unique code can be used for the fluid and the solid elements with variable materials coefficients and the global system of equations can be build using standard finite element assembly procedure treating the fluid and the solid elements without any distinction. The basic mesh-moving strategy proposed here consisted in moving the solid and interface nodes following their material velocity, while keeping the rest of the nodes fixed, defining a "semi-Lagrangian" step. Mesh smoothing was applied in a separate step, posterior to solving the monolithic system.
The model was compared with the results of 2D and 3D test cases from the literature. It was confirmed that the use of the unified formulation based on a single kinematic variable facilitates convergence of the iterative linear solvers and results in a feasible easy-to-implement monolithic approach.
While the use of the model in conjunction with conventional Bi-conjugate gradient linear solver provided acceptable convergence rates only for time step sizes corresponding to Courant number C ≤ 1, the use of the model in conjunction with the algebraic multigrid solver led to attractive convergence rates even for larger time steps.

Availability of Data and Materials
The data will be made available upon request.

Conflict of interest
The authors declare no conflict of interest.

Code Availability
The code is an open source code available at https://github.com/KratosMultiphysics/.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.