Ventricle-valve-aorta flow analysis with the Space–Time Isogeometric Discretization and Topology Change

We address the computational challenges of and presents results from ventricle-valve-aorta flow analysis. Including the left ventricle (LV) in the model makes the flow into the valve, and consequently the flow into the aorta, anatomically more realistic. The challenges include accurate representation of the boundary layers near moving solid surfaces even when the valve leaflets come into contact, computation with high geometric complexity, anatomically realistic representation of the LV motion, and flow stability at the inflow boundary, which has a traction condition. The challenges are mainly addressed with a Space–Time (ST) method that integrates three special ST methods around the core, ST Variational Multiscale (ST-VMS) method. The three special methods are the ST Slip Interface (ST-SI) and ST Topology Change (ST-TC) methods and ST Isogeometric Analysis (ST-IGA). The ST-discretization feature of the integrated method, ST-SI-TC-IGA, provides higher-order accuracy compared to standard discretization methods. The VMS feature addresses the computational challenges associated with the multiscale nature of the unsteady flow in the LV, valve and aorta. The moving-mesh feature of the ST framework enables high-resolution computation near the leaflets. The ST-TC enables moving-mesh computation even with the TC created by the contact between the leaflets, dealing with the contact while maintaining high-resolution representation near the leaflets. The ST-IGA provides smoother representation of the LV, valve and aorta surfaces and increased accuracy in the flow solution. The ST-SI connects the separately generated LV, valve and aorta NURBS meshes, enabling easier mesh generation, connects the mesh zones containing the leaflets, enabling a more effective mesh moving, helps the ST-TC deal with leaflet–leaflet contact location change and contact sliding, and helps the ST-TC and ST-IGA keep the element density in the narrow spaces near the contact areas at a reasonable level. The ST-SI-TC-IGA is supplemented with two other special methods in this article. A structural mechanics computation method generates the LV motion from the CT scans of the LV and anatomically realistic values for the LV volume ratio. The Constrained-Flow-Profile (CFP) Traction provides flow stability at the inflow boundary. Test computation with the CFP Traction shows its effectiveness as an inflow stabilization method, and computation with the LV-valve-aorta model shows the effectiveness of the ST-SI-TC-IGA and the two supplemental methods.


Introduction
We address the computational challenges encountered in and present results from ventricle-valve-aorta flow analysis. We include the left ventricle (LV) in the model because that makes the flow into the valve anatomically more realistic, which, in turn, makes the flow into the aorta more realistic. The challenges include (i) accurate representation of the boundary layers near the valve leaflets as they come into contact, (ii) computation with high geometric complexity, (iii) anatomically realistic representation of the LV motion, and (iv) flow stability at the inflow boundary, where we have a traction condition.
The ST-discretization feature of the ST-SI-TC-IGA provides higher-order accuracy compared to standard discretization methods. The VMS feature addresses the computational challenges associated with the multiscale nature of the unsteady flow in the LV, valve and aorta. The movingmesh feature of the ST framework enables high-resolution computation near the leaflets. The ST-TC enables movingmesh computation even with the TC created by the contact between the leaflets, dealing with the contact while maintaining high-resolution representation near the leaflets. The ST-IGA provides smoother representation of the LV, valve and aorta surfaces and increased accuracy in the flow solution. The ST-SI connects the separately generated LV, valve and aorta NURBS meshes, enabling easier mesh generation, and connects the mesh zones containing the leaflets, enabling a more effective mesh moving. It also helps the ST-TC deal with leaflet-leaflet contact location change and contact sliding and helps the ST-TC and ST-IGA keep the element density in the narrow spaces near the contact areas at a reasonable level.
The ST-SI-TC-IGA is supplemented with two other special methods. a) A structural mechanics computation method generates the LV motion from the CT scans of the LV and anatomically realistic values for the LV volume ratio. In this method, the structural mechanics computations, performed in different ways for the diastole and systole, generate a "table" of LV volumes and shapes. From that and the volume ratio given, we obtain the cardiac-cycle representation of the LV motion by using cubic B-splines in time and the ST-C [12]. b) The Constrained-Flow-Profile (CFP) Traction provides flow stability at the inflow boundary. This is done by placing adjacent to the inflow boundary a special-purpose element consisting of 27 basis functions. An SI connects the flow solutions over that element and the rest of the mesh. The special-purpose element, with only one unspecified controlpoint velocity at the inflow, results in a constrained flow profile, which is quadratic. The solution obtained for the unspecified velocity, together with the quadratic profile, represents the flow rate generated by the traction conditions specified at the inflow and outflow boundaries.
We first present a 2D-channel-flow test computation with the CFP to show how it performs compared to the analytical solution. We then present a computation with an LV-valve-aorta model created based on the CT scans of the LV and aorta.
In Sect. 2, we provide an overview of the ST-VMS and ST-SUPS. The overviews of the ST-SI, ST-TC, ST-SI-TC,  ST-IGA and ST-SI-TC-IGA are provided in Sects. [3][4][5][6][7]. The CFP Traction is presented Sect. 8, and the 2D-channel-flow test computation with it in Sect. 9. The LV-aorta-valve flow analysis is presented in Sect. 10, with the special structural mechanics computation method that generates the LV motion described in the first subsection. The concluding remarks are given in Sect. 11. The mesh moving method used in the LVaorta-valve flow computation is described in "Appendix A".

ST-VMS and ST-SUPS
This section, included for completeness, is mostly from [13]. The Deforming-Spatial-Domain/Stabilized ST (DSD/SST) method [14][15][16] was introduced for computation of flows with moving boundaries and interfaces (MBI), including fluid-structure interaction (FSI). In MBI computations the DSD/SST functions as a moving-mesh method. Moving the fluid mechanics mesh to follow an interface enables meshresolution control near the interface and, consequently, highresolution boundary-layer representation near fluid-solid interfaces. Because the stabilization components of the original DSD/SST are the Streamline-Upwind/Petrov-Galerkin (SUPG) [17] and Pressure-Stabilizing/Petrov-Galerkin (PSPG) [14] stabilizations, it is now called "ST-SUPS". The ST-VMS is the VMS version of the DSD/SST. The VMS components of the ST-VMS are from the residual-based VMS (RBVMS) method [18][19][20][21]. The ST-VMS has two more stabilization terms beyond those in the ST-SUPS, and the additional terms give the method better turbulence modeling features. The ST-SUPS and ST-VMS, because of the higherorder accuracy of the ST framework (see [3,4]), are desirable also in computations without MBI.
The ST-SUPS, ALE-SUPS, RBVMS, ALE-VMS and ST-VMS all have some embedded stabilization parameters that play a significant role (see [26]). These parameters involve a measure of the local length scale (also known as "element length") and other parameters such as the element Reynolds and Courant numbers. There are many ways of defining the stabilization parameters. Some of the newer options for the stabilization parameters used with the SUPS and VMS can be found in [2,5,6,10,92,93,[127][128][129][130]. Some of the earlier stabilization parameters used with the SUPS and VMS were also used in computations with other SUPG-like methods, such as the computations reported in [70,[131][132][133][134][135][136][137][138][139][140][141][142]. We will specify which ones we use here when we describe the computations in Sects. 9 and 10.
For more on the ST-VMS and ST-SUPS, see [26]. In the flow analyses presented here, the ST discretization provides higher-order accuracy compared to standard discretization methods. The VMS and SUPS features of the ST-VMS and ST-SUPS address the computational challenges associated with the multiscale nature of the unsteady flow in the LV, valve and aorta. The moving-mesh feature of the ST framework enables high-resolution computation near the leaflets.

ST-SI
This section, included for completeness, is mostly from [13]. The ST-SI was introduced in [6], in the context of incompressible-flow equations, to retain the desirable moving-mesh features of the ST-VMS and ST-SUPS in computations involving spinning solid surfaces, such as a turbine rotor. The mesh covering the spinning surface spins with it, retaining the high-resolution representation of the boundary layers, while the mesh on the other side of the SI remains unaffected. This is accomplished by adding to the ST-VMS formulation interface terms similar to those in the version of the ALE-VMS for computations with sliding interfaces [34,35]. The interface terms account for the compatibility conditions for the velocity and stress at the SI, accurately connecting the two sides of the solution. An ST-SI version where the SI is between fluid and solid domains was also presented in [6]. The SI in that case is a "fluid-solid SI" rather than a standard "fluid-fluid SI" and enables weak enforcement of the Dirichlet boundary conditions for the fluid. The ST-SI introduced in [7] for the coupled incompressible-flow and thermal-transport equations retains the high-resolution representation of the thermo-fluid boundary layers near spinning solid surfaces. These ST-SI methods have been applied to aerodynamic analysis of vertical-axis wind turbines [6,46,47], thermo-fluid analysis of disk brakes [7], flow-driven string dynamics in turbomachinery [114][115][116], flow analysis of turbocharger turbines [11,[117][118][119][120], flow around tires with road contact and deformation [2,13,109,121,122], fluid films [13,123], aerodynamic analysis of ram-air parachutes [124], and flow analysis of heart valves [1,106,108,110,111].
In the ST-SI version presented in [6] the SI is between a thin porous structure and the fluid on its two sides. This enables dealing with the porosity in a fashion consistent with how the standard fluid-fluid SIs are dealt with and how the Dirichlet conditions are enforced weakly with fluid-solid SIs. This version also enables handling thin structures that have T-junctions. This method has been applied to incompressibleflow aerodynamic analysis of ram-air parachutes with fabric porosity [124]. The compressible-flow ST-SI methods were introduced in [125], including the version where the SI is between a thin porous structure and the fluid on its two sides.
Compressible-flow porosity models were also introduced in [125]. These, together with the compressible-flow ST SUPG method [143], extended the ST computational analysis range to compressible-flow aerodynamics of parachutes with fabric and geometric porosities. That enabled ST computational flow analysis of the Orion spacecraft drogue parachute in the compressible-flow regime [125,126].
For more on the ST-SI, see [6,7]. In the CFP Traction, the ST-SI connects the special-purpose element placed adjacent to the inflow boundary to the rest of the mesh, and we will describe the CFP Traction in Sect. 8. In the LV-aorta-valve flow analysis, the ST-SI connects the separately generated LV, valve and aorta NURBS meshes, enabling easier mesh generation, connects the three mesh zones of the valve, enabling a more effective mesh moving, and connects the two separately generated parts of the LV mesh, again enabling easier mesh generation. By integrating it with the ST-TC, we deal with leaflet-leaflet contact location change and contact sliding, and we will describe the ST-SI-TC in Sect. 5. By integrating it with the ST-TC and ST-IGA, we keep the element density in the narrow spaces near the contact areas at a reasonable level, and we will describe the ST-IGA and ST-SI-TC-IGA in Sects. 6 and 7.

ST-TC
This section, included for completeness, is mostly from [13]. The ST-TC [8,9] was introduced for moving-mesh computation of flow problems with TC, such as contact between solid surfaces. Even before the ST-TC, the ST-SUPS and ST-VMS, when used with robust mesh update methods, have proven effective in flow computations where the solid surfaces are in near contact or create other near TC, if the nearness is sufficiently near for the purpose of solving the problem. Many classes of problems can be solved that way with sufficient accuracy. For examples of such computations, see the references mentioned in [8]. The ST-TC made moving-mesh computations possible even when there is an actual contact between solid surfaces or other TC. By collapsing elements as needed, without changing the connectivity of the "parent" mesh, the ST-TC can handle an actual TC while maintaining high-resolution boundary layer representation near solid surfaces. This enabled successful moving-mesh computation of heart valve flows [1,8,9,95,106,[108][109][110][111], wing clapping [8,100], flow around a rotating tire with road contact and prescribed deformation [2,13,109,121,122], and fluid films [13,123].
For more on the ST-TC, see [8,9]. In the LV-aorta-valve flow analysis, the ST-TC enables moving-mesh computation even with the TC created by the contact between the leaflets, dealing with the contact while maintaining high-resolution representation near the leaflets.

ST-SI-TC
This section, included for completeness, is mostly from [13]. The ST-SI-TC is the integration of the ST-SI and ST-TC. A fluid-fluid SI requires elements on both sides of the SI. When part of an SI needs to coincide with a solid surface, which happens for example when the solid surfaces on two sides of an SI come into contact or when an SI reaches a solid surface, the elements between the coinciding SI part and the solid surface need to collapse with the ST-TC mechanism. The collapse switches the SI from fluid-fluid SI to fluidsolid SI. With that, an SI can be a mixture of fluid-fluid and fluid-solid SIs. With the ST-SI-TC, the elements collapse and are reborn independent of the nodes representing a solid surface. The ST-SI-TC enables high-resolution flow representation even when parts of the SI are coinciding with a solid surface. It also enables dealing with contact location change and contact sliding. This was applied to heart valve flow analysis [1,106,108,110,111], tire aerodynamics with road contact and deformation [2,13,109,121,122], and fluid films [13,123].
For more on the ST-SI-TC, see [1,121]. In the computational analyses presented here, the ST-SI-TC enables contact location change and contact sliding between the leaflets.

ST-IGA
This section, included for completeness, is mostly from [13]. The success with IGA basis functions in space [23, 34,53,144] motivated the integration of the ST methods with isogeometric discretization, which we broadly call "ST-IGA." The ST-IGA was introduced in [3]. Computations with the ST-VMS and ST-IGA were first reported in [3] in a 2D context, with IGA basis functions in space for flow past an airfoil, and in both space and time for the advection equation. Using higher-order basis functions in time enables deriving full benefit from using higher-order basis functions in space. This was demonstrated with the stability and accuracy analysis given in [3] for the advection equation.
The ST-IGA with IGA basis functions in time enables a more accurate representation of the motion of the solid surfaces and a mesh motion consistent with that. This was pointed out in [3,4] and demonstrated in [10,96,98]. It also enables more efficient temporal representation of the motion and deformation of the volume meshes, and more efficient remeshing. These motivated the development of the ST/NURBS Mesh Update Method (STNMUM) [10,96,98], with the name coined in [93]. The STNMUM has a wide scope that includes spinning solid surfaces. With the spinning motion represented by quadratic NURBS in time, and with sufficient number of temporal patches for a full rotation, the circular paths are represented exactly. A "secondary mapping" [3,4,10,26] enables also specifying a constant angular velocity for invariant speeds along the circular paths. The ST framework and NURBS in time also enable, with the "ST-C" method, extracting a continuous representation from the computed data and, in large-scale computations, efficient data compression [5,7,12,109,[114][115][116]. The STN-MUM and the ST-IGA with IGA basis functions in time have been used in many 3D computations. The classes of problems solved are flapping-wing aerodynamics for an actual locust [10,26,96,97], bioinspired MAVs [94,95,98,99] and wing-clapping [8,100], separation aerodynamics of spacecraft [87], aerodynamics of horizontal-axis [43,[93][94][95] and vertical-axis [6,46,47] wind-turbines, thermo-fluid analysis of ground vehicles and their tires [5,109], thermo-fluid analysis of disk brakes [7], flow-driven string dynamics in turbomachinery [114][115][116], flow analysis of turbocharger turbines [11,[117][118][119][120], and flow analysis of coronary arteries in motion [112].
The ST-IGA with IGA basis functions in space enables more accurate representation of the geometry and increased accuracy in the flow solution. It accomplishes that with fewer control points, and consequently with larger effective element sizes. That in turn enables using larger time-step sizes while keeping the Courant number at a desirable level for good accuracy. It has been used in ST computational flow analysis of turbocharger turbines [11,[117][118][119][120], flowdriven string dynamics in turbomachinery [115,116], ram-air parachutes [124], spacecraft parachutes [126], aortas [106][107][108], heart valves [1,106,108,110,111], coronary arteries in motion [112], tires with road contact and deformation [2,13,122], and fluid films [13,123]. Using IGA basis functions in space is now a key part of some of the newest arterial zero-stress-state (ZSS) estimation methods [108,[145][146][147][148][149][150] and related shell analysis [151].
For more on the ST-IGA, see [11,26,96,124]. In the flow analyses presented here, the ST-IGA provides smoother representation of the LV, valve and aorta surfaces and increased accuracy in the flow solution.

ST-SI-IGA and ST-SI-TC-IGA
This section, included for completeness, is mostly from [13]. The ST-SI-IGA is the integration of the ST-SI and ST-IGA, and the ST-SI-TC-IGA is the integration of the ST-SI, ST-TC and ST-IGA. The turbocharger turbine flow [11,[117][118][119][120] and flow-driven string dynamics in turbomachinery [115,116] were computed with the ST-SI-IGA. The IGA basis functions were used in the spatial discretization of the fluid mechanics equations and also in the temporal representation of the rotor and spinning-mesh motion. That enabled accurate representation of the turbine geometry and rotor motion and increased accuracy in the flow solution. The IGA basis functions were used also in the spatial discretization of the string structural dynamics equations. That enabled increased accuracy in the structural dynamics solution, as well as smoothness in the string shape and fluid dynamics forces computed on the string.
The ram-air parachute analysis [124] and spacecraft parachute compressible-flow analysis [126] were conducted with the ST-SI-IGA, based on the ST-SI version that weakly enforces the Dirichlet conditions and the ST-SI version that accounts for the porosity of a thin structure. The ST-IGA with IGA basis functions in space enabled, with relatively few number of unknowns, accurate representation of the parafoil and parachute geometries and increased accuracy in the flow solution. The volume mesh needed to be generated both inside and outside the parafoil. Mesh generation inside was challenging near the trailing edge because of the narrowing space. The spacecraft parachute has a very complex geometry, including gores and gaps. Using IGA basis functions addressed those challenges and still kept the element density near the trailing edge of the parafoil and around the spacecraft parachute at a reasonable level.
The heart valve flow analysis [1,106,108,110,111] was conducted with the ST-SI-TC-IGA. The method, beyond enabling a more accurate representation of the geometry and increased accuracy in the flow solution, kept the element density in the narrow spaces near the contact areas at a reasonable level. When solid surfaces come into contact, the elements between the surface and the SI collapse. Before the elements collapse, the boundaries could be curved and rather complex, and the narrow spaces might have high-aspect-ratio elements. With NURBS elements, it was possible to deal with such adverse conditions rather effectively.
In computational analysis of flow around tires with road contact and deformation [2,13], the ST-SI-TC-IGA enabled a more accurate representation of the geometry and motion of the tire surfaces, a mesh motion consistent with that, and increased accuracy in the flow solution. It also kept the element density in the tire grooves and in the narrow spaces near the contact areas at a reasonable level. In addition, we benefit from the mesh generation flexibility provided by using SIs.
In computational analysis of fluid films [13,123], the ST-SI-TC-IGA enabled solution with a computational cost comparable to that of the Reynolds-equation model for the comparable solution quality [123]. With that, narrow gaps associated with the road roughness [13] can be accounted for in the flow analysis around tires.
An SI provides mesh generation flexibility in a general context by accurately connecting the two sides of the solution computed over nonmatching meshes. This type of mesh generation flexibility is especially valuable in complex-geometry flow computations with isogeometric discretization, removing the matching requirement between the NURBS patches without loss of accuracy. This feature was used in the flow analysis of heart valves [1,106,108,110,111], turbocharger turbines [11,[117][118][119][120], and spacecraft parachute compressible-flow analysis [126].
For more on the ST-SI-TC-IGA, see [1,2]. In the computations presented here, the ST-SI-TC-IGA is used for the reasons given and as described in the earlier paragraphs of this section.

Constrained-Flow-Profile traction
In computational flow analysis, typically, we specify the velocity at the inflow boundary, and the traction at the outflow boundary. Sometimes, however, we may want to specify the traction at both the inflow and outflow. To discuss how that would work, we first write the expressions for the flow rate, influx energy at the inflow, and the outflux energy at the outflow: Here u is the velocity, n is the unit normal vector, Γ in and Γ out are the inflow and outflow boundaries, ρ is the density, and p is the pressure. We note that a spatially-constant pressure value of p 0 can be added to or subtracted from the pressure field without changing the energy balance. Even though the flow rate would correspond to the pressure difference between the inflow and outflow, without specifying the full velocity profile at the inflow, the energy coming into the domain would not be under control. With the energy influx not being under control, for the same flow rate we could get different flow fields. An uncontrolled energy influx can also occur at the outflow if the flow reverses, and an outflow stabilization method was introduced [36] to remedy that. The method basically removes the incoming portion of the energy flux.
Targeting the inflow boundary, we introduce a simple method: Constrained-Flow-Profile (CFP) traction. Figure 1 illustrates the concept. We place adjacent to the inflow boundary a special-purpose element consisting of 3 n sd basis functions, where n sd is the number of space dimensions. An SI connects the flow solutions over that element and the rest of the mesh. The special-purpose element, with only one unspecified control-point velocity at the inflow, results in a constrained flow profile, which is quadratic. The solution obtained for the unspecified velocity, together with the quadratic profile, represents the flow rate generated by the traction conditions specified at the inflow and outflow bound- aries. As shown in [123], having a single quadratic NURBS element across a narrow channel is enough to represent the flow, which will be a low-Reynolds-number flow, and what we do here was inspired by that.

Test computation with the CFP
Using a 2D channel flow as the test case, we compare the solution obtained with the CFP to the analytical solution and the solution obtained with the outflow stabilization method [36] applied to the inflow, which we will identify with the acronym "OS." For completeness, we provide from [36] the outflow stabilization term, written in the ST framework: Here P OS is the part of the ST-slice lateral boundary where we apply the outflow stabilization method, v is the velocity of the lateral boundary, w is the test function, and a superscript "h" indicates that the function is coming from a finite-dimensional space. The function {A} − = A−|A| 2 . The parameter β ≥ 1 2 for stability [27].

Problem setup and analytical solution
The symbols H and L represent the channel height and length, and p in and p out the inflow and outflow pressures (see Fig. 2). At low Reynolds number, the solution is known to be of the form where μ is the viscosity and − H 2 x 2 H 2 , and the pressure gradient can be expressed as From that, the flow rate is and the energy influx where α = 54 35 , which is the kinetic energy factor. Depending on the flow conditions, the factor will change as a function of the Reynolds number The total energy in the domain is and the average over the domain

Meshes and boundary conditions
We set H = 1.3×10 −2 m, L = 1.365×10 −1 m, ρ = 1050 kg/m 3 , μ = 4.2×10 −3 Pa s, p in = 1.31 Pa, and p out = 0. The boundary conditions are no-slip on the channel walls and traction conditions corresponding to p in and p out at the inflow and outflow. We use two kinds of meshes, a regular mesh for the OS, and the kind shown in Fig. 1 for the CFP.   Except for the part of the domain where the CFP mesh has only one element, the two meshes have identical refinements, with increasing mesh refinement in the normal direction as we get closer and closer to the walls. Figure 3 shows the meshes. The number of control points and elements for the two meshes are given in Table 1.

Computational conditions
We start the computations with zero velocity, and suddenly apply the traction boundary conditions corresponding to p in and p out . We use the ST-VMS, with the stabilization parameters given by Eqs. (4)-(9) in [2]. The time-step size is 8.6×10 −3 s. The number of nonlinear iterations per time step is 3, and the number of GMRES iterations per nonlinear iteration is 300.

Results
Without the CFP or OS method, the solution quickly diverges. Figure 4 shows the velocity magnitude for the OS and CFP. We note that in both the CFP and OS computations, at the outflow boundary we use the outflow stabilization method [36] with β = 0.5. The maximum velocity from the analytical solution is 4.84×10 −2 m/s. With the CFP, we achieve, as expected, a parabolic profile at the inflow boundary, which is the form in the analytical solution. Figure 5 shows p total . Figure 6 shows the cross-sectional average pressure. In the OS computation the outflow stabilization term cancels the kinetic energy influx, and that causes the substantial discrepancy between the computed pressure and p in . In the CFP computation the discrepancy is very small. The discrepancy could be because of the sudden transition from a mesh with a single element to a mesh with many elements, creating a large discrepancy between the integration accuracies on the two sides of the SI. This can be remedied by having additional SIs and thus having a milder transition.

Ventricle-valve-aorta flow analysis
The flow computation model consists of the LV, aortic valve with sinuses, and the aorta. We do not include the mitral valve in the model. The boundary between the LV and the left atrium becomes our inflow boundary. We use the CFP Traction there when the mitral valve is open, and zero-velocity when it is closed. The aorta main outlet is treated as a regular outflow boundary with traction condition, and the three smaller outlets are treated as prescribed-velocity outflow boundaries.

Geometry
The entire model is shown in Fig. 7. The quadratic NURBS meshes for the three parts are generated separately and the SIs connect the three solution parts. The LV shrinks and expands with a cardiac cycle of T = 0.9 s, the valve opens and closes, and the aorta remains stationary. We describe the three parts more in the next three subsections.

LV geometry and motion
We build the LV geometry and motion based on CT scans of the LV at an instant in the cardiac cycle and anatomically realistic values for the volume ratio defined as (15) and ejection fraction (EF) defined as The volume ratio is given in Fig. 8 (from [152]), EF = 70 % (from [153]), and we set V max = 107 m .
By starting from the CT geometry and using a shell model, we perform a sequence of structural mechanics computations to obtain the cardiac-cycle LV shapes used in building the motion. In the diastole, using the CT geometry as the ZSS configuration, we apply a gradually-increasing spatiallyuniform pressure until the LV volume reaches V max . In the structural mechanics computations performed in the systole, we start from the CT geometry with zero-pressure load. The ZSS configurations are defined by applying to the CT geometry a factor expressed as λ −2 , where λ is the uniform stretch, with gradually-decreasing values of λ. This is done by applying the factor to the metric tensor of the ZSS (see Section 2.2 in [149]). We decrease the stretch until the LV volume comes down to V min .
The structural mechanics computations generate a "table" of LV volumes and shapes. From that and the volume ratio given in Fig. 8, we obtain the cardiac-cycle representation of the LV motion by using cubic B-splines in time and the ST-C [12]. The cycle periodicity is assured with the procedure followed for the same purpose in [98,112].
The shell formulation is from [151], and the thickness is variable, as shown in Fig. 9. Fung material model is used with the exponent factor D 2 = 8.365, and the pressure applied is scaled as p * = p 2D 1 D 2 . In the computations, we fix some of the control points, shown in Fig. 10, to obtain anatomically realistic shapes. Figure 11 shows the LV at its most expanded stage. Figure 12 shows how the LV volume varies with the applied pressure. Figure 13 shows the LV at its most shrunk stage. Figure 14 shows how the LV volume varies as the stretch decreases. We truncate the LV structure, roughly where the arrows in Fig. 10 indicate, to create the actual LV part used in the fluid mechanics computation. The truncated segment still provides us some guidance in building the valve part and the segment of the aorta part connecting to the valve.   The starting fluid mechanics volume mesh is generated using the method described in [117]. Then the mesh is moved to conform to the LV shapes obtained in the structural mechanics computations. The mesh moving method is the one described in "Appendix A", with (J M ) 0 = 1 and χ = 1. The mesh moving generates a "table" of LV volumes and fluid mechanics volume meshes. From that and Fig. 15 Valve-mesh motion with the SSDM. Template mesh (left) and simple shape (right) prior to the transformation. The lines are the element boundaries in the template mesh and control lines in the simple shape Fig. 16 Valve-mesh motion with the SSDM. Simple shape (left) and template mesh (right) after the transformation. The lines are the element boundaries in the template mesh and control lines in the simple shape the volume ratio given in Fig. 8, we obtain the cardiac-cycle representation of the fluid mechanics mesh motion by using cubic B-splines in time and the ST-C.

Valve motion
The valve-mesh motion is obtained by transformation from the mesh motion used in the valve computation reported in [106]. We note that the mesh motion in [106] is based on the ST-SI-TC-IGA, which deals with the TC created by the contact between the leaflets while maintaining highresolution representation near them. The transformation is achieved with the simple-shape deformation model (SSDM) [96], and in that context the mesh motion from [106] serves as a template. The process, illustrated in Figs. 15 and 16, starts with the template mesh. A simple shape, represented by far fewer control points, is created to enclose the template mesh. Then the simple shape transforms to have its upper and lower surfaces match the aorta and LV geometries. From that, the transformation of the template mesh is obtained by projection. The transformation with the SSDM is applied also to the motion of the valve leaflets, which are modeled in the template mesh, generating the cardiac-cycle representation of the leaflet motion by linear functions in time. The  Fig. 17 The three mesh parts and the special-purpose element (brown) added for the CFP Traction. The lines are the element boundaries projection takes places as many times as the number of flowcomputation time steps in the cardiac cycle.

Aorta
The aorta geometry is based on a different set of CT scans and is represented by cubic T-splines. From that, the fluid mechanics volume mesh is generated using the method described in [117]. The aorta part remains rigid during the fluid mechanics computation.

Mesh, boundary conditions, and blood properties
To the mesh composed of the three parts, we add, for the CFP Traction, an SI and a special-purpose element consisting of 27 basis functions. Figure 17 shows the entire mesh. The number of control points and elements are given in Table 2. Figure 18 shows the SIs and boundary conditions. There are altogether 7 SIs. Two connect the three main parts of the mesh, one connects the special-purpose element to the rest of the mesh, three connect the three zones of the valve mesh, and one is for mesh generation convenience in the LV. The boundary condition on the valve surfaces and arterials walls is no-slip. The tractions specified at the inflow and outflow boundaries correspond to p in = 10.5 kPa and  p out = 0. The estimated flow rates at the prescribed-velocity outflow boundaries come from the Murray's law [154], a flow rate distribution proportional to D 3 , where D is the average diameter of the outflow cross-section (see [107]). We set ρ = 1050 kg/m 3 and μ = 4.2×10 −3 Pa s. Figure 19 shows the mesh at different instants in the cardiac cycle. Figures 20  and 21 show the mesh in the valve at the same instants.

Computational conditions
We use the ST-SUPS, with the stabilization parameters given by Eqs. (4)-(9) in [2], and the ST-SI-TC-IGA. The time-step size is 2.81×10 −3 s. The number of nonlinear iterations per time step is 3, and the number of GMRES iterations per nonlinear iteration is 300.   Fig. 19. The right pictures are the zoomed views. The checkerboard pattern is for differentiating between the NURBS elements, and the colors are for differentiating between the NURBS patches Figure 22 shows the flow patterns. The flow patterns demonstrate that we are able to capture the spiral flow through the valve and in the aorta and we have a reasonable flow field even when the leaflets come into contact. Figure 23 and 24 show the magnitude of the wall shear stress (WSS) on the leaflet lower and upper surfaces. The WSS is, as expected, high on the lower surface.

Concluding remarks
We have addressed the computational challenges encountered in LV-valve-aorta flow analysis and presented results from the computation performed. We included the LV in the tion, and connects the mesh zones containing the leaflets, enabling a more effective mesh moving. 7. The ST-SI also helps ST-TC deal with leaflet-leaflet contact location change and contact sliding and helps ST-TC and ST-IGA keep the element density in the narrow spaces near the contact areas at a reasonable level. 8. The special structural mechanics computation method generates the LV motion from the CT scans of the LV and anatomically realistic values for the LV volume ratio. 9. The CFP Traction provides flow stability at the inflow boundary. This is done by placing adjacent to the inflow boundary a special-purpose element, and an SI connects the flow solutions over that element and the rest of the mesh. The special-purpose element, with only one unspecified control-point velocity at the inflow that we solve for, results in a constrained flow profile that represents the flow rate generated by the traction conditions specified at the inflow and outflow boundaries.
We first presented a low-Reynolds-number 2D-channelflow test computation with the CFP Traction and showed that the results were very close to the analytical solution. This demonstrated the effectiveness of the CFP Traction as an inflow stabilization method. We then presented a computation with an LV-valve-aorta model created based on the CT scans of the LV and aorta. We were successful in capturing the spiral flow through the valve and in the aorta and we had a reasonable flow field even when the leaflets came into contact. This demonstrated the effectiveness of the ST-SI-TC-IGA and the two supplemental methods in a highly challenging computational cardiovascular flow analysis. 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://creativecomm ons.org/licenses/by/4.0/.

A Mesh moving method and Mesh-Jacobian-based stiffening
In the mesh moving method introduced in 1992 [ [155][156][157], the motion of the internal nodes is determined by solving the equations of linear elasticity. As the boundary condition, the normal velocity of the mesh at the boundaries and interfaces is specified to match the velocity of the fluid. With the notation of [26], the mesh moving formulation can be written as Here ε ε ε is the strain tensor, σ σ σ is the stress tensor,ŷ h is the displacement from the reference configuration X h to the current configuration x h = X h +ŷ h , andŷ h t is the displacement from X h to thet configuration x h t = X h +ŷ h t we compute the mesh motion from. The most common selection ist = t n when we are computing the mesh at t n+1 .
In [155][156][157] the mesh deformation is dealt with selectively based on the sizes of the elements. Selective treatment based on element sizes is attained by altering the way we account for the Jacobian of the transformation from the element domain to the physical domain: where ξ ξ ξ is the parametric coordinate, and the subscript "M" is to clarify that this is the mesh Jacobian. The objective is to stiffen the smaller elements, which are typically placed near solid surfaces, more than the larger ones. When this method was first introduced in [155][156][157], it consisted of simply dropping J M from Eq. (17). That is equivalent to modifying Eq. (17) as where (J M ) 0 was originally intended to be a global scaling value and is now seen as a free parameter. The modification results in the smaller elements being stiffened more than the larger ones, and the element stiffening factor is (J M ) 0 J M . The method was given the name "Jacobian-based stiffening" in [158]. It was also augmented in [158] to a more extensive kind by introducing a stiffening power χ that determines the degree by which the smaller elements are rendered stiffer than the larger ones: This approach, when χ = 1, would be identical to the method first introduced in [155], and most of the time χ = 1.

Remark 1
To also clarify that the "Jacobian" in the name of the method is the mesh Jacobian, we are renaming the method "Mesh-Jacobian-based stiffening" (MJBS).
From earlier tests with differentt settings, we know that settingt = t n gives reasonably good performance for most problems. We believe that is partly because when we compute the mesh motion from a configuration close to the current configuration, J M provides the method feedback on which elements to protect.

Remark 2
While computing from thet = t n configuration works well for most problems, there are two closely related drawbacks: (i) the method is path-dependent, (ii) once elements accumulate in some region, it is hard to undo that. The path-dependence leads to non-cyclic results in FSI computations that we expect to have cyclic or near-cyclic results. Using the large-deformation elasticity equations (see, for example, [85]) is one way of addressing this issue. For FSI problems with cyclic or near-cyclic solutions, as an alternative, we introduce here the "Back-cycle-based mesh moving" (BCBMM) method. With T being the cycle period, in the kth cycle, as we move from t n to t n+1 , we compute the mesh motion from the configuration att = t n − (k − k FORW − 1)T ort = t n+1 − (k − k FORW − 1)T , where k ≥ 2. We start with k FORW = 0. With that, in any cycle, we compute the mesh motion based on the first cycle. In later cycles, as needed, we can set k FORW to higher values to compute the mesh motion based on the recent past cycles. That would be a good way if there are significant differences between the FSI solutions at different cycles before the solution becomes cyclic or nearcyclic.
We note that in practical computations the element stiffening method proposed in [159] is approximately equivalent to the MJBS method with χ = 1. To show that, we first approximate (J M ) 0 J M as where V is the element volume and V 0 is a global scaling value. The approximation becomes exact for simplex elements. In [159], the stiffening factor, represented by "1+τ m ", is given as where V min and V max are the minimum and maximum element volumes. This can also be written as In practical computations, V max >> V min , and for most of the elements that matter, Therefore Because V max − V min is a global scaling value just like (J M ) 0 is, the stiffening factor given by Eq. (25) is equivalent to the one given by Eq. (21).

Remark 3
We note that a stiffening-factor expression in terms of the mesh Jacobian gives us the option of evaluating the factor at the integration points, without the need for any precalculation. We also note that such an expression is naturally applicable to isogeometric discretization.