Heart valve isogeometric sequentially-coupled FSI analysis with the space–time topology change method

Heart valve fluid–structure interaction (FSI) analysis is one of the computationally challenging cases in cardiovascular fluid mechanics. The challenges include unsteady flow through a complex geometry, solid surfaces with large motion, and contact between the valve leaflets. We introduce here an isogeometric sequentially-coupled FSI (SCFSI) method that can address the challenges with an outcome of high-fidelity flow solutions. The SCFSI analysis enables dealing with the fluid and structure parts individually at different steps of the solutions sequence, and also enables using different methods or different mesh resolution levels at different steps. In the isogeometric SCFSI analysis here, the first step is a previously computed (fully) coupled Immersogeometric Analysis FSI of the heart valve with a reasonable flow solution. With the valve leaflet and arterial surface motion coming from that, we perform a new, higher-fidelity fluid mechanics computation with the space–time topology change method and isogeometric discretization. Both the immersogeometric and space–time methods are variational multiscale methods. The computation presented for a bioprosthetic heart valve demonstrates the power of the method introduced.


Introduction
In addressing the computational challenges of heart valve fluid-structure interaction (FSI) analysis with high-fidelity flow solutions, in this article we introduce an isogeometric sequentially-coupled FSI (SCFSI) method. The SCFSI method (see [1][2][3] and references therein) enables dealing with the fluid and structure parts individually at different steps of the solutions sequence, and also enables using different methods or different mesh resolution levels at different steps. In the isogeometric SCFSI analysis here, the first step is a previously computed (fully) coupled Immersogeometric Analysis (IMGA) FSI of the heart valve from [4], which has a reasonable flow solution. With the valve leaflet and arterial surface motion coming from that, we perform a new, higher-fidelity fluid mechanics computation with a space-time (ST) computational method composed of core and special ST methods. The core component is the ST Variational Multiscale (ST-VMS) method [5][6][7], which subsumes its precursor "ST-SUPS" (see Sect. 1.4) and shares the residual-based VMS (RBVMS) [8][9][10][11] feature with the IMGA [4]. Beyond the ST-VMS, the key components are the ST Topology Change (ST-TC) method [12,13], which is used in combination with the ST Slip Interface (ST-SI) method [14,15], and the ST Isogeometric Analysis (ST-IGA) [5,16,17]. Integration of these components, resulting in the ST-SI-TC [18] and ST-SI-TC-IGA [19,20] methods, gives us the increased scope and accuracy we want in FSI analysis when there is contact between the moving solid surfaces or other TC. With the SCFSI method based on the IMGA and ST-SI-TC-IGA, we conduct an FSI analysis of a bioprosthetic heart valve (BHV).

Moving-mesh and nonmoving-mesh methods
Using the terminologies and categorizations used in [21][22][23], a method for flows with moving boundaries and interfaces (MBI) can be an interface-tracking (moving-mesh) method or an interface-capturing (nonmoving-mesh) method, or possibly a combination of the two. In a moving-mesh method, as the interface moves and the fluid mechanics domain changes its shape, the mesh moves to adjust to the shape change and to follow (i.e. "track") the interface.
Moving-mesh methods require mesh update methods. Mesh update typically consists of moving the mesh for as long as possible and remeshing as needed. With the key objectives being to maintain the element quality near solid surfaces and to minimize the frequency of remeshing, a number of advanced mesh update methods [24][25][26][27][28] were developed to be used with the ST-SUPS, including those that minimize the deformation of the layers of small elements placed near solid surfaces. Some of these methods have also been used with other moving-mesh methods. The advanced mesh update methods developed more recently [12,13,16,[29][30][31][32] have been used mostly with the ST-VMS, and some of the methods are unique to the ST framework (see Sect. 1.8).
Moving the fluid mechanics mesh to follow a fluid-solid interface enables us to control the mesh resolution near the interface, have high-resolution representation of the boundary layers, and obtain accurate solutions in such critical flow regions. These desirable features do not come easily or do not come at all with the nonmoving-mesh methods. In these methods, the interface geometry is somehow represented over a nonmoving fluid mechanics mesh, with more accuracy in some methods than in some others, but the key point is that the fluid mechanics mesh does not move to follow the interface. Because the mesh is not following the interface, independent of how accurately the interface geometry is represented, the boundary layer resolution will be limited by the fluid mechanics mesh resolution where the interface is.
In many cases, finding the mesh update too challenging in an MBI problem is the reason for using a nonmoving-mesh method. Naturally, different researchers will have different thresholds for finding the mesh update too challenging. For some, even just a mesh moving without any remeshing could be too much. For some, the need for remeshing, no matter how small the associated cost is, could be too much. Some may find a near contact between solid surfaces or other near TC too much. Yet, some can bear it until the "end," giving up on the mesh update methods only when there is an actual contact or other TC, and even then, perhaps only partially. Some do not give up even then (see Sect. 1.6). What we need to keep in mind in all this is that, as pointed out in [33], "while it is understandable that fixed-mesh methods become more favored when the interface geometric complexity appears to be too high for a moving-mesh method, we need to remember that there is a difference between making the problem computable and obtaining good fluid mechanics accuracy near the interface." As mentioned in [22], one of course recognizes that certain classes of interfaces (such as free-surface and twofluid flows with splashing) might be too complex to handle with an interface-tracking technique and, therefore, for all practical purposes, require an interface-capturing technique. The Mixed Interface-Tracking/Interface-Capturing Technique (MITICT) [27] was introduced in 2000 for computation of MBI problems that involve both fluid-solid interfaces that can be accurately tracked with a moving-mesh method and fluid-fluid interfaces that are too complex or unsteady to be tracked. Such fluid-fluid interfaces are captured over the mesh tracking the fluid-solid interfaces.
Thinking similarly about MBI problems with an actual contact or other TC, the Fluid-Solid Interface-Tracking/Interface-Capturing Technique (FSITICT) was introduced in 2009 [3] as the FSI version of the MITICT. In the FSITICT, we track the interface we can with a moving mesh, and capture over that moving mesh the interfaces we cannot track, specifically the interfaces where we need to have an actual contact between the solid surfaces.

Computational challenges of heart valve FSI
Heart valve FSI analysis is one of the computationally challenging cases in cardiovascular fluid mechanics. The challenges include unsteady flow through a complex geometry, solid surfaces with large motion, and contact between the valve leaflets. Because the flow has to be completely blocked at contact, this is a case of actual contact. The heart valve IMGA FSI analysis reported in [4] was with actual contact and was conducted in the framework of the FSITICT.

SCFSI
The SCFSI method was introduced [34,35] as an approximate FSI method in the context of arterial FSI, with the name Sequentially-Coupled Arterial FSI (SCAFSI). In the SCAFSI, first we compute a "reference" (i.e. "base") arterial deformation as a function of time, driven only by the blood pressure, which is given as a function of time by specifying the pressure profile in a cardiac cycle. Then we compute a sequence of updates involving mesh motion, fluid dynamics calculation, and recomputing the arterial deformation. In [34,35] the method was in early stages of its development, the description was rather cursory, and the test computations were limited. A more extensive description of the method was provided in [1], together with a wider set of test computations.
Multiscale versions of the SCAFSI were introduced in [1], and the test computations were presented for the temporally multiscale version, using different time-step sizes for the structural and fluid mechanics parts. In the spatially multiscale versions proposed in [1], fluid mechanics meshes with different refinement levels are used at different stages of the FSI computation. We use a relatively coarse mesh at the early stages and reserve the highly-refined mesh for the stage where we plan to do the high-fidelity fluid mechanics computation. The spatially multiscale versions introduced in [1] included SCAFSI M1C. In that version, we first compute the time-dependent structure shape with the (fully) coupled FSI method and a relatively coarse fluid mechanics mesh, followed by mesh motion and fluid mechanics computation with a more refined mesh, obtaining a higher-fidelity fluid mechanics solution. Test computations with the spatially multiscale versions were first reported in a book chapter [36] and a conference paper [37], and then in journal articles [2,3]. The more general name "SCFSI" was introduced in [3,36,37], and SCFSI M2C was introduced in [3,37]. In SCFSI M2C, we first compute the time-dependent flow field with the (fully) coupled FSI method and a relatively coarse structural mechanics mesh, followed by a structural mechanics computation with a more refined mesh, obtaining a higher-fidelity structural mechanics solution.

ST-VMS and ST-SUPS
The Deforming-Spatial-Domain/Stabilized ST (DSD/SST) method [28,38,39] was introduced for computation of flows with MBI, including FSI. In flow computations with MBI, the DSD/SST functions as a moving-mesh method, possessing the associated desirable features. Because the stabilization components of the original DSD/SST are the Streamline-Upwind/Petrov-Galerkin (SUPG) [40] and Pressure-Stabilizing/Petrov-Galerkin (PSPG) [38] stabilizations, it is now called "ST-SUPS." The ST-VMS [5][6][7] is the VMS version of the DSD/SST. It 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 higher-order accuracy of the ST framework (see [5,6]), are desirable also in computations without MBI.
For completeness, we will include the ST-VMS in "Appendix A". The ST-SUPS, ALE-SUPS, RBVMS, ALE-VMS and ST-VMS all have some embedded stabilization parameters that play a significant role (see [22]). 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 [7,14,16,30,105,130,[136][137][138][139]. 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 [86,[140][141][142][143][144][145][146][147][148][149][150][151]. We will specify in "Appendix B" which stabilization parameters we use in the computation reported in this article.

ST-SI
The ST-SI was introduced in [14], 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 [52,53]. 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 [14]. 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 [15] 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 [14,63,64], thermo-fluid analysis of disk brakes [15], flow-driven string dynamics in turbomachinery [122][123][124], flow analysis of turbocharger turbines [17,[125][126][127][128], flow around tires with road contact and deformation [18,119,[129][130][131], fluid films [131,132], aerodynamic analysis of ram-air parachutes [133], and flow analysis of heart valves [19,20,116,118].
In the ST-SI version presented in [14] 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 [133]. The compressible-flow ST-SI methods were introduced in [134], 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 [134]. These, together with the compressible-flow ST SUPG method [152], 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 [134,135].
For completeness, we will include the ST-SI in "Appendix A". The interface terms in the ST-SI also involve element length, in the direction normal to the interface. We will specify in "Appendix B" which element length we use for that in the computation reported in this article.

ST-TC
The ST-TC [12,13] 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 [12]. 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 [12,13,19,20,107,116,118,119], wing clapping [12,32], flow around a rotating tire with road contact and prescribed deformation [18,119,[129][130][131], and fluid films [131,132].

ST-SI-TC
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 fluid-solid 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 [19,20,116,118], tire aerodynamics with road contact and deformation [18,119,[129][130][131]131], and fluid films [131,132].

ST-IGA
The success with IGA basis functions in space [42,52,70,153] motivated the integration of the ST methods with isogeometric discretization, which we broadly call "ST-IGA." The ST-IGA was introduced in [5]. Computations with the ST-VMS and ST-IGA were first reported in [5] 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 [5] 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 [5,6] and demonstrated in [16,29,109]. 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) [16,29,109], with the name coined in [30]. 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" [5,6,16,22] 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 [7,15,119,[122][123][124]154]. The STNMUM 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 [16,22,29,108], bioinspired MAVs [106,107,109,110] and wing-clapping [12,32], separation aerodynamics of spacecraft [101], aerodynamics of horizontal-axis [30,60,106,107] and vertical-axis [14,63,64] wind-turbines, thermo-fluid analysis of ground vehicles and their tires [7,119], thermo-fluid analysis of disk brakes [15], flow-driven string dynamics in turbomachinery [122][123][124], flow analysis of turbocharger turbines [17,[125][126][127][128], and flow analysis of coronary arteries in motion [120].

ST-SI-IGA and ST-SI-TC-IGA
The The turbocharger turbine flow [17,[125][126][127][128] and flowdriven string dynamics in turbomachinery [123,124] 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 [133] and spacecraft parachute compressible-flow analysis [135] 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 [19,20,116,118] 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 [130,131], the ST-SI-TC-IGA enables 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 keeps 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 [131,132], the ST-SI-TC-IGA enables solution with a computational cost comparable to that of the Reynolds-equation model for the comparable solution quality [132]. With that, narrow gaps associated with the road roughness [131] 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 [19,20,116,118], turbocharger turbines [17,[125][126][127][128], and spacecraft parachute compressibleflow analysis [135].

Outline of the remaining sections
In Sect. 2, we present the SCFSI analysis of the BHV. The concluding remarks are given in Sect. 3. The ST-VMS and ST-SI are given in "Appendix A", and the stabilization parameters in "Appendix B".

SCFSI analysis of a BHV
In the SCFSI analysis, the first step is a (fully) coupled IMGA FSI computation [4] of the BHV, with a cardiac cycle of T = 0.86 s. The computation was driven by a prescribed left-ventricular pressure at the inlet and a resistance boundary condition at the outlet, and the corresponding flow rate was obtained as part of the FSI solution. With the BHV and

Geometry
The model, shown in Fig. 1, has three leaflets and a metal frame. In the IMGA FSI computation [4], the BHV is represented with cubic T-splines, and the arterial surface with quadratic NURBS. Figure 2 shows the BHV and arterialsurface meshes.

Surface mesh in the ST-SI-TC-IGA computation
The structure model is of zero thickness in geometry. Even in the closed position of the valve, there are small gaps between the model surfaces, including the metal surfaces. However, the flow solver, because of the way the IMGA deals with zero-thickness geometries, in combination with the limited mesh refinement, could see the gaps as closed. In our case, we want the flow solver to see the model surfaces accurately so that the moving-mesh method can do what it is good at -enable high-resolution flow representation near those surfaces. Our method for closing the gaps is to add just enough thickness to the surface geometries coming from the IMGA computation. We select the thickness to be small so that our  [4] model geometry is as close to the IMGA geometry as possible. It is roughly 10 times smaller than the structure thickness used in the IMGA computation. We start with a thickness of 4.921×10 −3 cm, determined based on closing the gaps between the metal surfaces, which are not deforming. That being our maximum thickness, we reduce it locally to prevent overlap between the model surfaces, but only when we need to do so during the cardiac cycle, in a time-varying fashion.  Figure 3 shows the surfaces added on the two sides, giving the structure thickness. We project the surfaces created on the two sides to quadratic-NURBS representation, and then connect the surfaces with a rounding arch. The rounding arch and the surfaces are connected in a tangent form. With the rounding, as an added benefit, we can capture the small-scale flow patterns near the free edges of the leaflets. The structural displacements from the IMGA computation are projected to the original surface and that is applied to the corresponding points of the surfaces on the two sides. Figures 4 and 5 show the valve surfaces. Figure 6 shows the artery quadratic NURBS surfaces.

Mesh and inflow conditions
We create a template mesh with three SIs. The mesh has three parts, identified as Part 1, Part 2 and Part 3. Figure 7 shows the three SIs and the three parts of the mesh. The number of control points and number of elements are 429,780 and 289,452. Part 1 faces the SIs and will have the elements that will collapse and become reborn by the motion of the leaflets. Part 2 remains unchanged during the computation. Part 3 is the mesh between Part 1, Part 2 and the arterial wall. The motion of Part 1, following the leaflet motion, is generated with a method taking into account the contact. The motion of Part 3 is generated automatically, by solving the steady-state structural mechanics equations based on the neo-Hookean model with Jacobian-based stiffening [7,22,[24][25][26]162], following Part 1, the leaflet motion and the artery wall motion. Figures 8 and 9 show the mesh motion.
The boundary conditions are no-slip on the valve and the arterial wall, traction-free at the outflow boundary, and uniform velocity at the inflow boundary. We use the flow rate from the IMGA computation. When the valve is closed, the flow rate at the inflow boundary and the time derivative of the closed-space volume do not match. Figure 10 shows that discrepancy. We note that the volume we are time-differentiating is the volume of the mesh zone that covers the closed space, which we will call "time-differentiation volume." In dealing with the discrepancy, in the closed-valve time periods, we set the flow rate equal to the time derivative of the closed-space volume. In narrow time zones neighboring the closed-valve periods, we set the flow rate equal to a value blended between   Figure 11 shows the modified flow rate. Figure 12 shows the inflow velocity corresponding to the modified flow rate.

Computational conditions
We use the ST-SUPS (see "Appendix A", which contains also the ST-SI), with the stabilization parameters given in "Appendix B". The time-step size is 5.00×10 −3 s. The number of nonlinear iterations per time step is 3, and the number of GMRES iterations per nonlinear iteration is 300. Figures 13 and 14 show the flow patterns from the second cardiac cycle. The figures demonstrate that we are able to capture Fig. 12 Inflow velocity the sheet-like flow structure emanating from the leaflet edges and we have a reasonable flow field even when the leaflet surfaces come into contact. Figure 15 shows the corresponding wall shear stress (WSS) on the valve surfaces. The WSS is high around the leaflet edges and on the upstream-side of the leaflet surfaces.

Concluding remarks
We have introduced an isogeometric SCFSI method that can address some of the toughest computational challenges of heart valve FSI analysis with high-fidelity flow solutions. Beyond the challenges related to the unsteady flow through a complex geometry and the solid surfaces with large motion, we have addressed the challenges related to the contact between the valve leaflets with a SCFSI combination of immersogeometric and ST computational methods. Both the immersogeometric and ST methods are variational multiscale methods. Because the ST computational method is a moving-mesh method, it enables high-resolution representa-

A ST-VMS and ST-SI
For completeness, we include, mostly from [14,18], the ST-VMS and ST-SI methods.

A.1 ST-VMS
The ST-VMS is given as where are the residuals of the momentum equation and incompressibility constraint. Here, ρ, u, p, f, and h are the density, velocity, pressure, body force, and the stress specified at the boundary. The stress tensor is defined as σ σ σ = −pI+2με ε ε(u), where I is the identity tensor, μ = ρν is the viscosity, ν is the kinematic viscosity, and ε ε ε (u) = (∇ ∇ ∇u) + (∇ ∇ ∇u) T /2 is the strain-rate tensor. The test functions associated with the u and p are w and q. A superscript "h" indicates that the function is coming from a finite-dimensional space. The symbol Q n represents the ST slice between time levels n and n + 1, (P n ) h is the part of the slice lateral boundary associated with the boundary condition h, and n is the spatial domain at time level n. The superscript "e" is the ST element counter, and n el is the number of ST elements. The functions are discontinuous in time at each time level, and the superscripts "−" and "+" indicate the values of the functions just below and above the time level.

Remark 1
The ST-SUPS can be obtained from the ST-VMS by dropping the eighth and ninth integrations.

A.2.1 Two-side formulation (fluid-fluid SI)
In describing the ST-SI, labels "Side A" and "Side B" will represent the two sides of the SI. The ST-SI version of the formulation given by Eq. (1) includes added boundary terms corresponding to the SI. The boundary terms for the two sides are first added separately, using test functions w h A and q h A and w h B and q h B . Then, putting together the terms added to each side, the complete set of terms added becomes where Here, (P n ) SI is the SI in the ST domain, n is the unit normal vector, v is the mesh velocity, γ = 1, and C is a nondimensional constant. The element length h will be defined in "Appendix B2".

A.2.2 One-side formulation (fluid-solid SI)
On solid surfaces where we prefer weak enforcement of the Dirichlet conditions [49,51] for the fluid, we use the ST-SI version where the SI is between the fluid and solid domains. This version is obtained (see [14]) by starting with the terms added to Side B and replacing the Side A velocity with the velocity g h coming from the solid domain. Then the SI terms added to Eq. (1) to represent the weakly-enforced Dirichlet conditions become The element length h B will be given in "Appendix B2".

B.1 ST-VMS
There are various ways of defining the stabilization parameters τ SUPS and ν LSIC . Here, τ SUPS is mostly from [137]: The first and second components are given as and where r is the solution-gradient direction: Here G ST and G are the ST and space-only element metric tensors: wherê The ST and space-only Jacobian tensors are and where θ and ξ ξ ξ are the temporal and spatial parametric coordinates. The transformation tensor D ST is defined as The definitions used for D θ and D play an important role, especially for higher-order isogeometric discretization [137,139] and simplex elements [138]. However, in this article, we use D θ = 1 and D = I. The third component, originating from [7], is defined as where · F represents the Frobenius norm. The stabilization parameter ν LSIC is from [30]: where h LSIC is set equal to the minimum element length h MIN : For more ways of calculating the stabilization parameters in flow computations, see [86,136,[140][141][142][143][144][145][146][147][148][149][150][151].

B.2 ST-SI
The element length used in the ST-SI is given as