Study on the Accuracy of Structural and FSI Heart Valves Simulations

The performance of heart valves, either native or artificial, can be evaluated by means of finite element analyses, either from a structural or a fluid–structure interaction (FSI) point of view. The latter captures the coupling between the valve leaflets and the blood in a more realistic way. The selection of the appropriate finite elements approach for the model is the first and fundamental step to achieve accurate simulations. The aim of this work is to investigate the influence of the type, formulation, size, and shape of the elements in heart valves simulations. The effects related to the choice of the finite elements-shell or solid- in structural and FSI simulations were analyzed. In particular, the analysis of grid convergence on both the structure and fluid domains, the influence of the element typology, formulation and damping factor in an idealized three-leaflets valve model loaded with physiological pressure conditions were investigated. Stress values and valve kinematics results confirmed the importance of performing a proper verification process for selecting the most appropriate elements with the optimal accuracy to computational cost ratio. In this regard, our results indicate the quadrangular shell with reduced integration and viscous hourglass control to be the best choice to model heart valves. If a solid discretization is required, quadratic hexahedral elements with full integration are also acceptable. Finally, our results show that the damping coefficient needs to be carefully selected in order to smooth out the high frequency modes of the structure without introducing excessive numerical artificial viscosity.


INTRODUCTION
Computational analysis has become one of the most widely used tools to investigate the human physiological and pathological conditions, to help in clinical decisions and interventional planning, and for the development and optimization of new devices. 16 Computer modeling can also be used for in silico clinical trials, where ''virtual'' treatments and clinical procedures are applied to ''virtual'' patients. 49 This recent line of research can potentially reduce the need for clinical trials by replacing them with reliable numerical models. Numerical models, if used as clinical predictors, need to be reliable and this entails the concept of verification and validation (V&V). 48 Therefore, the first step in defining and building a mathematical model is the verification process summarized by the following questions: ''Does the computational model accurately solve the underlying mathematical model?'' 49 and ''Does the software do what it is designed to do?. 48 The confidence assessment of a simulation is nothing new in the field of computational models. 30 In this regard, the American Society of Mechanical Engineering provides guidelines for verification and validation in computational solid mechanics 2,4 and fluid-dynamics. 3 Verification and validation are two different concepts normally used in conjunction to demonstrate the reliability and accuracy of an in silico model; the former aims to prove that the implemented model is correct, the latter that the problem is correctly modeled. As Roache 53 suggested, ''Verification deals with mathematics; validation deals with physics''. A number of reviews on this topic with recommendations and suggestions for the implementation of V&V analyses can be found in the literature. 6,45,46,53 In the field of cardiovascular applications, the fascination for numerical models to investigate heart valves is well known. 16,29,55 The recent literature involves a number of examples of numerical models of heart valves 39,65 and involving transcatheter valve technology. 60 In this context, the importance of verifying the ability of the numerical models to represent the reality is evident. 29 Without investigating whether the software solves the equations correctly, as most of the commercial solvers are already validated, issues related to the minimization of errors associated with mesh resolution, and differences in the quantities of interest obtained with different modeling strategies are important for investigation. Mesh resolutions, type, order and distribution of the elements are key factors determining the integrity and accuracy of the solution. In the past, models of heart valves have been discretized with triangular 12,13,35 or quadrangular shell, 1,7,14,15,25,27,42,43,47,58,62-64 tetrahedral 23 or hexahedral 9,37,38,41,56 elements with both linear and quadratic order, along with, both reduced and full integration formulation. Most of these studies correspond to finite element (FE) structural analysis in which the pressure exerted by the fluid is uniformly applied on the leaflets. However, the best numerical approach capable of reproducing the loading on the valve leaflets due to the fluid coupling is the fluidstructure interaction (FSI) approach. 36 The main issue related to the solution of an FSI algorithm is how to manage the interface between the solid structure and the fluid domain, which involves the solution of a non-linear problem. Different interface tracking techniques exist, which can be grouped in two typologies: the non-boundary fitted methods and the boundary-fitted methods. 9 In the non-boundary fitted methods, the interface movement is calculated through an interpolation technique without deforming the fluid mesh: the solid mesh results superimposed to the fluid mesh. The most widely used methods belonging to this group are the immersed boundary 51,52 and the fictitious domain. 24 They are well suited for problems where the structural domain undergoes large displacements, as is the case with biological or natural heart valves. On the contrary, in the boundary-fitted methods, the interface movement is obtained through the deformation of the fluid mesh, which continuously fits the solid domain. The Arbitrary Lagrangian-Eulerian (ALE) method 19,26,28 is the most common method in this group. Although, ALE guarantees very accurate results near the interface between the solid structures and the fluid, there are limitations associated with this methodology when used for heart valve simulations. 61 Apart from being computationally expensive due to the necessary remeshing step after each iteration, the deformation of the mesh may lead to excessive distortion of the fluid elements and the appearance of negative element volumes, which may cause convergence problems. In addition, ALE methods do not allow the fluid domain to be partitioned in two independent regions as happens during the closure of the valve.
In this context, the aim of this work is to verify and compare different technical details of heart valve simulations that are often neglected, namely, the convergence of the mesh, the finite element typology and formulations, and the damping factor. To the best of the authors' knowledge, this is the first study on the accuracy of heart valve simulations that simultaneously investigates all of these parameters. In this study, the verification process was performed with structural and FSI simulations on an idealized three-leaflets valve model, loaded with physiological pressure conditions. The aims are 3-fold: (i) To conduct a convergence grid analysis separately on both the structure solid domain (the valve) and the coupled fluid domain (the blood); (ii) To investigate the influence of the element typology and formulation, and of the damping factor on the kinematics and the stress field of the heart valve; and, (iii) To compare the corresponding structural and FSI simulations, highlighting their differences in terms of the load acting on the leaflets of the valve.

Governing Equations
The multi-physics FSI approach considers the interaction between structure and fluid domains. For the solid domain, based on the solution of the continuum mechanics equations, the finite element method discretizes in space the conservation of linear momentum equation: where q s is the density of the body, v s i the velocity of the material points, r s ij the Cauchy stress tensor, and f s i the body forces per unit volume. For the fluid domain, the computational fluid dynamic software solves the conservation of mass and the Navier-Stokes equations: where q f is the density of the fluid, v f i the velocity of the fluid particles at position x i , r f ij the Cauchy stress tensor, and f f i the fluid forces per unit volume.
The FSI algorithm entails that the fluid and the structures continuously couple forces and displacements through the wet surface. The partitioned approach solves the solid Eq. (1) and the fluid Eqs. (2-3) with two separate solvers, different from the monolithic approach characterized by a unique solver for both domains. Hence, a coupling algorithm to track the interface (the wet surface) is required. In addition, two additional boundary conditions are added at the interface: (i) continuity of the velocity field at the interface (no-slip condition), and (ii) the balance of the traction forces exchanged between the fluid and the structure domains: where t f ¼ r f n f and t s ¼ r s n s are the traction forces at the interface in the fluid and solid domains, respectively, with n s ¼ Àn f the normal to the wet surface.

Idealized Model
An idealized standard tri-leaflets aortic valve was adopted as a model. It was composed of three identical deformable leaflets and a rigid supporting corona (Fig. 1a). The internal and external diameters of the valve were 23 and 28 mm respectively, with a homogenous leaflets thickness of 0.4 mm; the total height of the device was 17 mm. The valve was considered to be made of pericardium which was assumed as linear elastic since, for the working strain range, the material behavior is well represented with the initial tangent modulus of the hyperelastic characteristic curve. 34 The material properties of the valve were 3 MPa for the elastic modulus, 0.49 for the Poisson's ratio and 1100 kg/m 3 for the density. 40 The fluid domain consisted of a rigid tube with a diameter of 27 mm, and a total length of 135 mm, (Fig. 1b). The blood was modeled as a Newtonian fluid, with a density of 1060 kg/m 3 and a dynamic viscosity of 3.5 cP. 37 A physiologic pressure difference was applied to the valve to reproduce two complete cardiac cycles each lasting 0.8 s (Fig. 1c). The ventricular pressure tracing registred 113.58 and 7.39 mmHg as maximum and minimum values, whereas the aortic pressure tracing registered 111.88 and 66.72 mmHg, respectively. The difference between upstream and downstream valve pressures presented a maximum value of 1.95 mmHg, during systole, and À 81.07 mmHg during diastole. 18 In the structural analyses, the pressure difference (discontinuous green line in Fig. 1c) was applied directly on the leaflets (aortic side). On the contrary, in the FSI analyses the pressure difference was applied to the inlet (ventricular side) whereas a zero-pressure condition was set to the aortic outlet section. A no-slip condition was applied to the external fluid nodes since the rigid ''container'' was simulated via nodal constraints. A self-contact with penalty algorithm among the three leaflets was defined. All of the numerical models were solved with an explicit algorithm.
All of the discretized models were created with Hypermesh 2017 (Altair Engineering, Inc., USA) and the simulations were performed on 16 CPUs of an Intel-MPI Xeon64 using the commercial finite element solver LS-Dyna 971 Release 10.1 (LSTC, Livermore CA, USA).

Mesh Sensitivity-Structural Analysis
Structural analyses were carried out to perform convergence studies for both of the considered element technologies: shells and bricks (specified in the model acronyms with the letters S and B, respectively).
Three different shell meshes, from coarse (S 1 ) to fine (S 3 ), were generated to perform a sensitivity analyses on the discretization density (Fig. 2a). For this preliminary analysis, a quadrangular shell element with reduced integration and viscous hourglass control was used. A mass proportional damping was assigned to the structure, with different values during the opening and closing phases. For the continuum solid models, three brick meshes (B 1 , B 2 , and B 3 ) were generated, based on the previous surface meshes S 1 , S 2 , and S 3 , with 4 elements across the leaflets thickness. A sensitivity analysis on the number of elements across the thickness was also performed, by considering the mesh B 2 with two (B 2 -T 2 ) or eight (B 2 -T 8 ) elements across the thickness. For the bricks, a fully integrated quadratic element with different values of mass damping in the cardiac cycle was used. The details of the element size and the total number of elements for each case are summarized in Table 1.
Geometric opening area (GOA) during the systolic opening phase and the average value of the maximum principal stresses in the central part of the leaflets were used to monitor the convergence of the mesh. The region where the first principal stresses were averaged was chosen to be far from the applied boundary conditions (commissure edges near the corona) and the contact area (free edges) ( Fig. 2a-white circle). FIGURE 2. Different surface meshes considered in the convergence study, from coarse (S 1 ) to fine (S 3 ) refinement with the region where the first principal stresses were averaged (white circle) (a); GOA during systolic phase for the shell (b) and the brick (d) models; first principal stress for the shell (c) and the brick models (e) in the entire cardiac cycle.

Element Formulation-Structural Analysis
Five shell element formulations were tested as reported in Table 2. The reference formulation (model S 2 ) from the grid sensitivity analysis was a Belytschko-Lin-Tsay (BLT) element with one-point integration 10 and viscosity hourglass control. 31 The BLT element with stiffness hourglass control, S 2 -HgS, was also tested for completeness. The computational efficiency of the BLT formulation is due to the mathematical simplification of the co-rotational and velocity-strain. The third element tested was based on the Belytschko-Leviathan formulation, S 2 -BL, which adds hourglass viscosity stresses to the physical stresses at element level. 22 The fourth element type was a fully integrated shell-type formulation, S 2 -FI, with an assumed strain interpolants used to alleviate the locking problem and enhance in-plane bending. 20 All the previous formulations are based on the Reissner-Mindlin kinematic assumption, in which mid-surface displacement with rotations are used to describe the shell deformation (5parameter shell model). 8 The last element formulation tested was a thickness enhanced one point integration formulation, S 2 -T, also called thick-thin 6-or 7parameter shell model. 8 In all these cases, only the element and hourglass control was changed, while keeping fixed all the other parameters of the simulation. For the shell models, mass proportional damping coefficients of 1.0 and 5.0 were set for the systolic and diastolic phases respectively. When weighted system damping is used, a force vector due to system damping, F n damp , is added to the external and internal loads where m and v are the mass and the velocity field of the body, and D s a damping constant. As D s affects the results, especially for the kinematics of the valve, two models based on the reference model S 2 were considered, a low damped model (S 2 -D 0.1 ) and a high damped model (S 2 -D 5 ) with values of D s , equal to 0.1 and 5, respectively.
For the solid models meshed with brick elements, three different element formulations were investigated as reported in Table 2. The reference formulation (model B 2 ), adopted also in the mesh sensitivity analysis, was a quadratic fully integrated element with nodal rotation. 50 This quadratic formulation is free of any locking problem. The second formulation was a tri-linear fully integrated solid element, B 2 -FI, in which the pressure is constant within the element to avoid pressure locking. 44 However, shear locking remains a problem since this formulation is not able to capture the pure bending mode of deformation. As a third formulation, an advanced version of the B 2 -FI element, the B 2 -FI Adv , was considered. In this improved formulation, the transverse shear locking is reduced at expense of a higher computationally cost. For solid elements, we also investigated the use of reduced integration with a viscous, B 2 -RI-HgV, or stiffness hourglass control with exact volume integration, B 2 -RI-HgS,. 22 For the solid models, mass proportional damping coefficients of 0.1 and 1.0 were set for the systolic and diastolic phases, respectively. In addition, for the B 2 reference model, the influence of the mass proportional damping was tested by setting the damping coefficient to a constant value of 0.05 (B 2 -D 0.05 ).
The average GOA during the systolic phase and the average maximum principal stress value of the central part of the leaflets along the cardiac cycle were evaluated as primary and integrated variables. The total computation time was also reported.

FSI Analysis
The ''operator split'' Lagrangian-Eulerian approach, a non-boundary fitted method implemented in the solver LS-Dyna, was adopted to define the interaction between the solid structure (valve) and the fluid (blood). In this approach, the Eulerian fluid Eqs. (2-3) are ''split'' into Lagrangian and advection steps. 11 In the Lagrangian step, the solid mesh deforms due to forces transmitted from the fluid, whereas during the advection step the velocity is remapped and interpolated in the fixed reference grid. This coupling algorithm prevents flow through the structure by applying penalty forces to the fluid and the structures.
After the convergence analysis of the solid part, the mesh for the valve is fixed and a convergent study on the fluid domain is performed. With the non-boundary fitted method, the nodes of the fluid and solid mesh do not necessarily coincide. However, a necessary condition is that the structural domain remains immersed in the fluid domain during the whole simulation. In addition, in order to avoid leakage through the valve leaflets, an appropriate number of coupling nodes is required. This condition is achieved if the element size of the fluid domain is nearly the same as, or smaller than, the solid elements where the coupling (wet surface) takes place. Three different fluid meshes were generated to perform a sensitivity analysis, from coarse (S 2 -FSI 1 ) to fine (S 2 -FSI 3 ) refinement (Table 3). The mesh of the fluid domain was generated considering the symmetry of the three valve leaflets. In addition, the fluid grid of the S 2 -FSI 2 model, which has the same element size for the valve and the fluid grids, was generated with the classical quadratic o-grid outline (model S 2 -FSI 2Q ), to investigate its influence on the results, especially on the pressure load transmitted to the structure from the fluid. Details about the element size and the total number of fluid elements for each model are listed in Table 3. Concerning the valve, the reference shell-type model S 2 (first row in Table 2) was used for the convergence analysis of the FSI simulations. As for the fluid domain, the fluid element formulation was set as Eulerian one-point ALE multimaterial elements. The coupling between the fluid and the structure was defined within a penalty-based algorithm. It consisted of the imposition of a resistance force on the slave nodes, typically the Lagrangian structure, which was proportional to their penetration on the master surface, the fluid. The GOA during the systolic opening phase and the average first principal stress value of the central part of the leaflets at the end of diastole were used to monitor the convergence of the mesh. To reduce the computational time, the bulk modulus of the blood was reduced to 1% of its real value. This modification, however, does not significantly affect the pressure and velocity fields. 5, 32,33 The structural analysis performed with model S 2 and the corresponding FSI analysis with the converged fluid grid were compared. According to the results from a previous study, 36 the pressure difference during the systolic phase was scaled by a factor of two (S 2 -FSI 2 -Sf 2 ) or three (S 2 -FSI 2 -Sf 3 ) in order to emphasize the difference in terms of pressure distribution on the leaflets between the structural and the FSI simulation. The diastolic pressure difference was, on the contrary, kept identical. The FSI analysis was also carried out with brick solid elements for the valve with both the physiological and scaled pressure boundary conditions (models B 2 -FSI 2 , B 2 -FSI 2 -Sf 2 and B 2 -FSI 2 -Sf 3 ). Table 4 summarizes the FSI models investigated. The comparison was performed in terms of average GOA during the systolic phase, blood flow rate, stroke volume, and total computing time. Table 5 summarizes the results obtained for all simulations in terms of maximum GOA at peak systole, average max principal stress at the central part of the leaflets at the end of diastole, and total computational time for one cardiac cycle. The three shell grids (S 1 , S 2 and S 3 ) showed an acceptable performance in terms of the monitored variables (see Figs. 2a, 2b, 2c). The main variable for evaluating the convergence was the kinematics of the leaflets, in terms of the averaged GOA of the valve during the systolic phase. The difference in percentage of the GOA between the medium (S 2 ) refinement and the coarse (S 1 ) was of 0.39%, while between the medium (S 2 ) and the fine (S 3 ) was of 3.45% (Fig. 2b). In terms of the cycle-averaged difference of the average maximum principal stress in the central part of the leaflets (white circle in Fig. 2a), between the medium (S 2 ) refinement and the coarse (S 1 ) was of 0.49%, whereas between the medium (S 2 ) and the fine (S 3 ) was of 10.77% (Fig. 2c). The fine mesh model (S 3 ) showed leaflet co-penetration during the closing phase, due to the excessive number of points in the penalty contact defined for the self-contact of the leaflets. This behavior leads to a larger difference in the results in terms of the stresses developed in the leaflets during the closing phase of the valve. For this reason, the medium mesh (S 2 ) was chosen for the remaining simulations. For the brick element meshes, the percentage difference in the GOA between the medium (B 2 ) and the coarse (B 1 ) meshes was of 10.70%; between the medium (B 2 ) and the fine (B 3 ) meshes, it was 0.54%; between the medium (B 2 ) and the medium B 2 with 2elements through thickness (B 2 -T 2 ), it was 60.94% and between the medium (B 2 ) and the medium B 2 with 8elements through thickness (B 2 -T 8 ), it was of 4.40% (Fig. 2c). Regarding the cycle-averaged maximum principal stress in the center of the leaflet, between the medium (B 2 ) and the coarse (B 1 ) meshes the percentage difference was of 2.01%; between the medium (B 2 ) and the fine (B 3 ) meshes, it was 15.63%; between the medium (B 2 ) and the medium B 2 with 2-element through thickness (B 2 -T 2 ), it was 16.96%; the simulation with the medium B 2 with 8-elements through thickness (B 2 -T 8 ) did not finish because of excessive leaflet co-penetration during the diastolic phase (Fig. 2d). This problem was also observed in the finer model B 3 . On the other hand, in the model with 2elements, through the thickness (B 2 -T 2 ) the valve did not fully open during the systolic phase. This was because the number of elements in the thickness was not enough to properly catch the bending behavior of the leaflet resulting in an overly stiff response. The mesh B 2 with 4-elements in the thickness was found to be a good compromise between the proper modeling of the bending behavior of the leaflet in systole and the prevention of contact problems during diastole. Therefore, it was chosen to perform the remaining simulations.

Element Formulation-Structural Analysis
The element formulation investigation showed how the different technical choices affect the kinematics of the valve and the stress field on the leaflets. The reference model (S 2 ) for the shell-based models was a reduced integration quadrilateral element with viscosity hourglass control. As regards the hourglass control formulation, the comparison between models S 2 , S 2 -BL and S 2 -HgS showed similar kinematic results (Fig. 3a), with model S 2 -BL resulting slightly stiffer (difference in systole-averaged GOA of less than 2%). The main difference was found in the cycle-averaged maximum principal stress in the center of the leaflet, with the S 2 -HgS reaching a value 23% larger with respect to the S 2 and S 2 -BL. This difference corresponded to much greater hourglass energy dissipation with respect to the S 2 and S 2 -BL models (Fig. 5a). The computing time was not influenced by the hourglass formulation (Table 5). The fully integrated shell model (S 2 -FI) and the thickness enhanced one-point integration model (S 2 -T) gave the same results as the reference  Sf scale factor. model (S 2 ), but with a significant increase in computing time (Table 5). As regards the influence of the damping factor, Fig. 3a demonstrates the substantial influence of this factor on the kinematics of the valve. In fact, in model S 2 -D 0.1 the valve opening time (VOT) reduces from 0.085 s, for the reference model, to 0.075 s, with a 2.4% larger maximum GOA. On the contrary, the damping effect in model S 2 -D 5 was found excessive since it prevented the valve from fully opening, with a maximum GOA of less than 60%, with respect to the reference model. This result is consistent with the increment in damping energy dissipation with the scale factor (Fig. 5d). For the solid cases, the reference model (B 2 ) was a quadratic fully integrated hexahedral element with four elements through the thickness. The tri-linear fully integrated model (B 2 -FI) gave a very stiff response, due to the shear locking problem that characterizes this type of formulation in bending dominated problems (Fig. 4a). The maximum GOA was of 47.69 mm 2 , with respect to 326.91 mm 2 of the B 2 model (see Table 5). The advanced tri-linear fully integrated formulation (model B 2 -FI Adv ) overcome this problem: the VOT was identical to the reference model, with a difference in the maximum value of  GOA of less than 9.50% (Fig. 4a). However, the computational cost increased five-fold with respect to the B 2 -FI model. Regarding the reduced integration models, B 2 -RI-HgV and B 2 -RI-HgS behaved very differently. Model B 2 -RI-HgV with the viscous hourglass control opened faster with a VOT of 0.07 s against the 0.085 s of the reference model but the same value of maximum GOA (see Table 5). On the contrary, the model with the stiffness hourglass control (B 2 -RI-HgS) behaved almost identically to the B 2 -FI model. The hourglass energy dissipation for the solidtype models resulted to be one order of magnitude larger than for the shell-type models (Fig. 5c). As regards the damping factor, it did not strongly influence the kinematic of the valve during systole (model B 2 -D 0.05 in Fig. 4a, damping energy dissipation curves in Fig. 5e). However, it caused instabilities in the contacts during the diastolic phase preventing the completion of the simulation. The comparison of the stress fields (Fig. 4b) demonstrated that the closing behavior of the valve was not strongly influenced by the element formulation; on the contrary, the computational time was affected. The tri-linear fully integrated model (B 2 -FI) was four times less expensive than the reference quadratic model (B 2 ), while the advanced linear fully integrated element resulted 50% more computationally expensive than the B 2 model. The reduced integration models (B 2 -RI-HgV and B 2 -RI-HgS) reduced the computational time with respect to the reference model (B 2 ) by 82 and 88%, respectively.
Finally, it was verified that a quasi-static condition in all the simulations was achieved as the ratio between the kinetic and the internal energy was less than 5% during all the simulated cardiac cycles (Fig. 5a).

FSI Analysis
The convergence analysis on the fluid meshes confirmed the general recommendation that the element size of the fluid and solid domains has to be nearly the same where coupling is to take place. In fact, the difference in the systole-averaged GOA between the reference FSI 2 model and the coarse FSI model (FSI 1 ) was 3.26%, whereas, between the FSI 2 model and the fine FSI model (FSI 3 ), the difference was only 0.04% (Fig. 6a). In addition, between model FSI 2 and the model meshed with the classic o-grid structures (FSI 2Q ) the difference in systole-averaged GOA was 2.99%. As regards the cycle-averaged maximum principal stress, results showed the same trend as the GOA, with percentage differences of 11.66% between the FSI 2 and the FSI 1 models, 9.97% between the FSI 2 and the FSI 3 models, and 9.60% between FSI 2 and FSI 2Q models (Fig. 6b). In particular, the comparison between models FSI 2 and FSI 2Q showed a non-negligible influence of the mesh topology on the results. Considering that models FSI 2 and FSI 2Q have the same element size, these results indicate that the mesh should respect the symmetry of the structure domain when possible in order to maximize computational efficiency and accuracy. Based on these results, the FSI 2 model with the same element size as the underlying structure grid was selected as the most appropriate.
The shell-type and brick-type models with the applied physiological and scaled pressure curves (Fig. 1c) were compared. In terms of the kinematics of the valve, the maximum opening area of the valve was globally lower in the FSI models than their equivalent structural models. In fact, the maximum value of GOA in the model S 2 -FSI 2 was 71.10 mm 2 , with respect to the 334.60 mm 2 of the reference structural model S 2 (Table 5). On the contrary, the scaled S 2 -FSI 2 -Sf 2 and S 2 -FSI 2 -Sf 3 models reached a value of 162.33 and 289.74 mm 2 , respectively. Accordingly, the resulting stress field in the leaflet (Fig. 7c), the blood flow rates (Fig. 7b), the velocity field, and the stroke volumes changed. In particular, the stroke volume varied from 4.95 mL for the S 2 -FSI 2 model, to 24.95 mL for the S 2 -FSI 2 -Sf 2 model, to 54.29 mL for the S 2 -FSI 2 -Sf 3 model. A similar result was found for the FSI models with the valve meshed using hexahedral elements. The maximum value of GOA in the model B 2 -FSI 2 was 11.42 mm 2 , with respect to the 326.91 mm 2 of the reference structural model B 2 (Fig. 7d). The scaled B 2 -FSI 2 -Sf 2 and B 2 -FSI 2 -Sf 3 models reached a value of 140.02 mm 2 and 285.69 mm 2 , respectively; the resulting stress field in the leaflet (Fig. 7f), blood flow rates (Fig. 7d), velocity field (Fig. 7g), and stroke volumes changed consistently, with the stroke volume varying from a value of 1.26 mL for the B 2 -FSI 2 model, to 14.92 mL for B 2 -FSI 2 -Sf 2 model and 46.37 mL for B 2 -FSI 2 -Sf 3 model. By considering the B 2 and B 2 -FSI 2 -Sf 3 models, the VOT increased from 0.085 s to 0.15 s, an increment of twofold. In general, with regard to the kinematics and the stress distribution, FSI models with a shell-based and brick-base structural meshes performed similarly (see Fig. 7). However, they greatly differ in terms of the computational time. For brickbased models, the increment in computational times with respect to shell-based modes was fourfold as shown in Table 5. This is, in great part due to the larger number of degrees of freedom associated with the brick type models.

DISCUSSION
When performing numerical simulations, the importance to assess the reliability and truthfulness of numerical models is of capital importance. 54 In this regard, the verification and validation (V&V) process is crucial. This is particularly true for heart valve simulations where different formulation and numerical techniques are used. 1,7,9,[12][13][14][15][16]23,25,27,29,[35][36][37][38][39][41][42][43]47,51,56,58,60,[62][63][64][65] The choice of the discretization technique is a fundamental step in the set-up phase of a numerical investigation; this includes the definition of the element topology, order of interpolation and type of integration, and the adequate number of nodes. Two main families of elements are commonly used to model heart valve: quadrilateral shell and hexahedral solid ele-ments. Quadrilateral shell elements and quadratic hexahedral solid elements are to be preferred, geometry permitting, since they can properly represent the bending mode of deformation which characterizes the valve mechanics. 17 Indeed, triangular and tetrahedral elements, which are constant-strain elements, do not properly model the bending behavior of the structure, in particular, tetrahedral elements which lead to an over stiff mechanical response. The same conclusion can be extended to tri-linear hexahedral elements which are prone to suffering shear-locking in bending dominated problems. These deficiencies associated with the tetrahedral and tri-linear hexahedral elements necessitate using discretization with a large number of elements across the thickness, which leads to a significant loss of computational efficiency against quadratic FIGURE 6. Different fluid meshes considered in the FSI convergence study, from coarse (FSI 1 ) to fine (FSI 3 ) refinements; GOA during the systolic phase (a) and the mean first principal stress for the valve in the central region of the leaflets for the entire cardiac cycle (see white circle in Fig. 2a) elements. Numerical techniques are also required to integrate various quantities over the volume of the element, e.g., to calculate the element stiffness matrix. In this regard, volume integration is carried out with Gaussian quadrature and can be full or reduced. Numerical issues affect both formulations. Full-integration with quadratic elements is computationally costly but results are adequate for bending dominated problems. On the contrary, reduced integration can eliminate the problem of shear-locking in tri-linear hexahedral elements but introduces zero energy modes, called hourglass modes, that need to be eliminated in the numerical formulation.
In the present study, the verification process was carried out in the cardiovascular numerical field, more precisely on biological heart valves numerical models. In literature, heart valves have been generally studied from a structural view point i.e., neglecting the inter-FIGURE 7. GOA during the systolic phase for the shell-type (a) and brick-type (d) FSI element formulation models; flow rate for the shell-type (b) and brick-type (e) FSI element formulation models; First principal stress field on the valve discretized with shell (c) and brick (f) elements at the moment of maximum valve opening; velocity field during the systolic phase in the brick-type element models (g). action with the blood. However, since their dynamic behavior is directly linked to hemodynamics, FSI analysis is gaining attention. 21,25,38,39,57,62 The aim of this study was to investigate how different technical details in the simulation of heart valves affected the numerical solution of the problem. We considered different scenarios of the heart valve modeling and verified which technical choices were suitable to accurately describe the problem. Our aim was not to compare the numerical model of the valve against in vitro tests since this was not a validation study. Structural and FSI simulations were performed, considering: different mesh typology, discretization, element formulation, damping factor, and applied boundary conditions. An idealized tri-leaflets linear elastic valve loaded with physiological pressure was considered.
Even though heart valve models discretized with triangular shell 12,13,35 and tetrahedral solid 23 elements are found in the literature, these element typologies were not investigated here because they are highly inefficient in bending dominated problems due to their constant strain characteristics. Hence, only quadrilateral shell and hexahedral elements were considered. The preliminary sensitivity analysis on the discretization refinement demonstrated that increasingly refining the mesh is not always the most efficient way to reach convergence when dealing with valve mechanics. While this could generally be true for the systolic phase, for the diastolic phase of the simulated cardiac cycle (closure of the valve), the finer meshes may exhibit problems of leaflets co-penetration during contact. This implies that excessively fine meshes may be counterproductive. In fact, in the penalty contact algorithm, for every slave node the solver searches the closest master node/segments, computes the orthogonal distance, and if penetration exists, then a force proportional to penetration depth is applied to both slave and master nodes. If the number of slave nodes is too large, the algorithm is unable to prevent the penetration between the slave and the master structure. Therefore, when dealing with valve mechanics, the optimal mesh size requires not only monitoring the variables of interest i.e., the GOA and the maximum principal stress but also verifying that the discretization adequately describes the contact between the leaflets during the diastolic phase.
For models based on shell elements with reduced integration, the viscous hourglass control resulted in the most appropriate model, since the stiffness hourglass control led to higher values of stress during the closure phase. The full integration and the thicknessenhanced models showed similar results, but with a significant increment in computational time. In fact, the greatest benefit of the reduced integration formulation is the significant saving in computational time. The cost of the fully integrated element may be justified by its higher reliability, and if used sparingly may actually increase the overall speed. Regarding the mass proportional damping, our results indicate that the damping coefficient has a strong impact on the kinematics of the valve. A little change in the damping coefficient value resulted in a significant change in the opening and closure time of the valve. For this reason, in shell-type heart valve simulations, a sensitivity analysis of this coefficient should be always conducted.
As regards the brick-based models, the trilinear hexahedral with full integration formulation showed shear locking since the opening and closing of the valve is dominated by bending. The advanced formulation of the trilinear fully integrated element overcomes this problem, but at additional computational cost. In this regard, the quadratic hexahedral element with full integration resulted in a good compromise between accuracy and computational cost because it was free of shear locking and allowed working with meshes with a contained number of nodes. The quadratic hexahedral element with reduced integration results were inadequate for simulating valve mechanics. The model based on quadratic hexahedral element with reduced integration resulted unstable during the closure of the valve when viscous hourglass control was used. On the contrary, when stiffness hourglass control was used, the model resulted as stiff as the model with trilinear hexahedral elements with full integration. Finally, unlike the shell-based models, the mass proportional damping did not influence the kinematics of the valve, due to the innate greater stability of the solid elements.
The chosen shell and brick element formulations were also analyzed for the FSI simulations with a nonboundary fitting method. The main advantage of this kind of method is that only the structure grid deforms, and the calculation of the variables is less expensive. 61 Thanks to the fixed fluid grid, this approach's results are well suited for problems where the structural domain undergoes large displacement as in the case of the heart valve dynamics. One of the weak points of this approach is that the results at the interfaces are less accurate. 9 This can be partially solved by refining the fluid mesh, but with a significant increase in computational time. In this regard, the results confirmed that at least the two grids (structural and fluid) should be of the same element size in order to adequately describe the valve mechanics. In addition, our results indicate that the mesh of the fluid domain should respect the symmetry of the structure domain when possible in order to maximize computational efficiency and accuracy of the load transfer between fluid and structure.
The direct comparison between structural and FSI analyses showed a great difference in terms of the GOA of the valve. This difference was already noticed in previous works. 33,36,59 The reason behind this difference is due to the pressure boundary conditions used for the structural simulations. On the contrary, in an FSI simulation, this pressure difference was applied between the inlet and outlet of the fluid domain in correspondence with actual boundary conditions. Therefore, the forces transmitted to the leaflets resulted not only lower for the FSI simulation than for the structural simulation but were also non-uniformly distributed on the leaflet. In other words, in the structural analysis, the pressure boundary conditions applied on the leaflets were overestimated. This explains not only the larger values of GOA obtained for the structural analysis, but also the faster valve opening with respect to FSI simulations conducted with the same element formulations for the solid domain. During the closure phase, all these differences were less evident due in part to the larger pressure gradient acting on the leaflets, but also to the fact that the leaflets enter in contact during this phase. A previous validation study from our group demonstrated that the structural analysis cannot fully reproduce a realistic loading on the valve. 36 Due to the presence of the fluid, which transmits the pressure to the structure, the FSI simulations are more appropriate than structural analyses to describe the real dynamic behavior of heart valves. In addition, the inertia of the fluid, not present in a structural analysis, delays the opening and the closure of the valve, as is the case in physiological conditions.
The study is not exempt from limitations, including the idealized, perfectly symmetric, linear elastic valve model to the rigid container of the fluid. However, these simplifications are reasonable for the study aim which was to show the influence of each investigated numerical aspects of a heart valve model. The obtained results can be generalized and applied to realistic and patient-specific heart valves; the selection of the best element type and size can also be adopted or more complex problems. The simplicity of the model was adopted for computational reasons since we wanted to carry out numerous simulations.
Moreover, during the simulations, the assumption that the material behavior is well reproduced with the first elastic modulus of the hyperelastic characteristic curve 34 was confirmed by the resulting strain field.
In conclusion, the goal of this study was to investigate the consequence of different numerical parameter choices. In this regard, modeling the valves with shells produced the best compromise between model complexity and computational efficiency. In particular, the quadrilateral shell with reduced integration and viscous hourglass control was the best adoption. If the geometry necessarily required a solid discretization technique, a quadratic hexahedral element with full integration can be used. Finally, the damping coefficient sensitivity analysis should be always performed in order to smooth the high frequency vibration of the structure without introducing numerical artificial viscosity.