Space–time flow computation with boundary layer and contact representation: a 10-year history

Abstract

In computation of flow problems with moving solid surfaces, moving-mesh methods such as the space–time (ST) variational multiscale method enable mesh-resolution control near the solid surfaces and thus high-resolution boundary-layer representation. There was, however, a perception that in computations where the solid surfaces come into contact, high-resolution boundary-layer representation and actual-contact representation without leaving a mesh protection opening between the solid surfaces were mutually exclusive objectives in a practical sense. The introduction of the ST topology change (ST-TC) method in 2013 changed the perception. The two objectives were no longer mutually exclusive. The ST-TC makes moving-mesh computation possible even without leaving a mesh protection opening. The contact is represented as an actual contact and the boundary layer is represented with high resolution. Elements collapse or are reborn as needed, and that is attainable in the ST framework while retaining the computational efficiency at a practical level. The ST-TC now has a 10-year history of achieving the two objectives that were long seen as mutually exclusive. With the ST-TC and other ST computational methods introduced before and after, it has been possible to address many of the challenges encountered in conducting flow analysis with boundary layer and contact representation, in the presence of additional intricacies such as geometric complexity, isogeometric discretization, and rotation or deformation of the solid surfaces. The flow analyses conducted with these ST methods include car and tire aerodynamics with road contact and tire deformation and ventricle-valve-aorta flow. To help widen awareness of these methods and what they can do, we provide an overview of the methods, including those formulated in the context of isogeometric analysis, and the computations performed over the 10-year history of the ST-TC.

1 Introduction

In computation of flow problems with moving solid surfaces, high-resolution representation of boundary layers is very often a requirement for obtaining reliable solutions. Moving-mesh methods such as the Space–Time Variational Multiscale (ST-VMS) [1,2,3] method meet that requirement because they enable mesh-resolution control near the solid surfaces. There was, however, a perception that in computations where the solid surfaces come into contact, high-resolution boundary-layer representation and actual-contact representation without leaving a mesh protection opening between the solid surfaces were mutually exclusive objectives in a practical sense. The introduction of the ST Topology Change (ST-TC) method in 2013 [4] changed the perception. The two objectives were no longer mutually exclusive. The ST-TC makes moving-mesh computation possible even without leaving a mesh protection opening. The contact is represented as an actual contact and the boundary layer is represented with high resolution. Elements collapse or are reborn as needed, and that is attainable in the ST framework while retaining the computational efficiency at a practical level. The ST-TC now has a 10-year history of achieving the two objectives that were long seen as mutually exclusive.

With the ST-TC and other ST computational methods introduced before and after, it has been possible to address many of the challenges encountered in conducting flow analysis with boundary layer and contact representation, in the presence of additional intricacies such as geometric complexity, isogeometric discretization, and rotation or deformation of the solid surfaces. The methods used with the ST-TC in addressing the challenges include the ST-VMS and its precursor “ST-SUPS” [5,6,7], which serve as the core methods. They also include the special methods introduced in connection with the core methods. The special methods are the ST Slip Interface (ST-SI) method [8, 9], ST Isogeometric Analysis (ST-IGA) [1, 10, 11], syntheses of these ST methods, such as the “ST-SI-TC-IGA” [12], local length scales and stabilization parameters targeting IGA discretization [13, 14], element-based mesh relaxation (EBMR) [15], mesh relaxation and mesh moving based on fiber-reinforced hyperelasticity [16], Complex-Geometry IGA Mesh Generation (CGIMG) method [17], and NURBS Surface-to-Volume Guided Mesh Generation (NSVGMG) method [18]. The flow analyses conducted with these ST methods include car and tire aerodynamics with road contact and tire deformation [18] and ventricle-valve-aorta flow with leaflet–leaflet contact [19]. The computations had high-resolution boundary-layer representation with actual-contact representation of the tire–road contact and leaflet–leaflet contact. To help widen awareness of these methods and what they can do, we provide an overview of the methods and the computations performed over the 10-year history of the ST-TC.

1.1 Core ST methods: ST-SUPS and ST-VMS

The Deforming-Spatial-Domain/Stabilized ST (DSD/SST) method, introduced in 1990 [5, 20, 21] as a moving-mesh method, is now called “ST-SUPS.” The DSD/SST was intended to serve as a core method for flows with moving boundaries and interfaces (MBI). This is a wide class of problems that includes fluid–structure interactions (FSI), and in a more general sense, moving solid surfaces. The new acronym “ST-SUPS” reflects the stabilization components of the method. Those components are the Streamline-Upwind/Petrov–Galerkin (SUPG) method [22, 23] and the Pressure-Stabilizing/Petrov–Galerkin (PSPG) method [5, 24]. In a broader interpretation of the terminology, as pointed out in [25], the ST-SUPS includes the stabilization component coming from Least-Squares on the Incompressibility Constraint (LSIC), as how it was for the DSD/SST in its form in [6]. That makes three stabilization terms for the ST-SUPS: SUPG, PSPG, and LSIC. The ST-VMS is the VMS version of the DSD/SST, with the VMS components coming from the residual-based VMS (RBVMS) method [26,27,28]. The ST-VMS subsumes its precursor ST-SUPS and has two more stabilization terms beyond those the ST-SUPS has. With these two additional terms, the ST-VMS has better turbulence modeling features. Both core ST methods are residual-based. Even in flow computations without MBI, the higher-order accuracy of the ST context (see [1, 2, 29]) would make the ST-SUPS and ST-VMS favored methods.

The ST framework of the core ST methods, in its basic functionality, is comparable to the arbitrary Lagrangian–Eulerian (ALE) method, which is older and more commonly used. Its finite element version goes at least as far back as 1981 [30]. Its use in 3D finite element flow computations with residual-based stabilized methods like the SUPG, however, is newer compared to the ST-SUPS. The earliest 3D computation with the ALE-SUPS was the ram-air parachute FSI in 1999 [31]. Among the moving-mesh methods with the VMS stabilization, the ALE-VMS method [32] precedes the ST-VMS. The ST-SUPS, ALE-SUPS, ALE-VMS, and ST-VMS have much in common. Therefore the classes of problems computed with these methods since their inception also have much in common. We see that when we list the classes of problems computed.

The classes of problems computed with the ALE-SUPS, RBVMS, and ALE-VMS include wind turbines [33,34,35,36], turbomachinery [37,38,39], stratified flows [40, 41], bridges [42], marine applications [43, 44], free-surface flows [45, 46], two-phase flows [47,48,49], additive manufacturing [50], aircraft applications [51], hypersonic flows [52], parachutes [31], cardiovascular medicine [53], mixed ALE-VMS Immersogeometric Analysis (IMGA) computations [54] in the framework of the Fluid–Solid Interface-Tracking/Interface-Capturing Technique [55], and IMGA FSI and flow analysis [56].

The classes of problems computed with the ST-SUPS and ST-VMS include those summarized in [57] (all computed in 1993–2018), wind turbines [58], turbomachinery [59], ground vehicles and tires [60], fluid films [61], disk brakes [9], flapping-wing aerodynamics [62], spacecraft [63], parachutes [64], cardiovascular medicine [19], Taylor–Couette flow [65], U-ducts [66], higher-order temporal isogeometric discretization [29], boundary-layer mesh resolution studies [25], and long-wake flows [67]. The core ST methods, together with the special methods introduced in connection with them, enabled first-ever solutions in some of the most challenging classes of MBI problems (see [68] for an overview).

1.2 Mesh moving methods introduced in connection with the core ST methods

Moving-mesh methods require mesh update, which, most of the time, consists of mesh moving. The mesh moves to adapt to the boundary and interface motion and to retain the mesh resolution near the solid surfaces as close to how it was in the starting mesh as possible. If the element distortion goes beyond a level that we find acceptable, the mesh update involves remeshing. Overall, the first priority of a good mesh moving method should be to maintain good element quality near the solid surfaces where we care about the accuracy in representing the boundary layer. The second priority should be to reduce the need for remeshing, which is of course more time-consuming than mesh moving and is also disruptive. Section 4.7 of [69] provides an overview of the basic mesh moving concepts, including the concept of special- and general-purpose mesh moving. Since the introduction of the ST-SUPS in 1990, a significant number of special- and general-purpose mesh moving methods, with increasing levels of sophistication, have been introduced in connection with the core ST methods. Among those, the first general-purpose method was the linear-elasticity mesh moving with Jacobian-based stiffening [70,71,72], which is now called, for reasons explained in [19], “mesh-Jacobian-based stiffening (MJBS).” An overview of those methods can be found in a recent review article [73] on the subject, and we will discuss some of the most recent general-purpose mesh moving methods in Sect. 1.9.

1.3 ST-SI

The ST-SI accurately connects the parts of the solution obtained over two mesh zones with nonmatching meshes at the SI between the zones. It is a multifunction method. Its original function was to enable the use of the core ST methods as moving-mesh methods even in the presence of rotating solid surfaces, such as a car tire, so that we continue benefiting from using a moving-mesh method. The mesh covering the rotating solid surface rotates with it, retaining the high-resolution representation of the boundary layers, while the mesh on the other side of the SI remains independent of the solid-surface rotation. In that functionality, the ST-SI is essentially the ST version of the “sliding interface” method introduced in [74] in connection with the ALE-VMS. Accurately connecting the two sides of the solution is achieved by adding to the core ST formulation a set of integrals over the SI. These integrals account for the compatibility conditions for the velocity and stress at the SI. The ST-SI versions going beyond the original functionality of the ST-SI version with the “fluid–fluid SI” work in comparable ways and further increase the scope of the core ST methods. The ST-SI version with “fluid–solid SI,” also introduced in [8], is basically the ST counterpart of the weak-Dirichlet-condition method [75]. The ST-SI version with the “porosity SI” [8] is for a thin porous structure with fluid on its two sides. There are more ST-SI versions, with various functions, and they will be summarized in Sects. 1.5–1.7. More on the ST-SI and a partial list of the classes of problems computed with it can be found in [76].

1.4 Near contact, actual contact, and the ST-TC

In flow computations where the moving solid surfaces come into contact, the nature of the problem dictates if the contact can be represented as a “near contact” or needs to be represented as an actual contact. If the nearness is close enough for obtaining a good flow solution, then representing the contact as a near contact would still leave a narrow mesh protection opening between the solid surfaces. With no topology change in the flow domain, mesh update in moving-mesh methods becomes less challenging. Many classes of problems can be solved that way. For example, in particle-laden flows [77,78,79], the collision between the particles can be represented as a “near collision,” and in spacecraft parachute FSI [15, 80], the contact between the canopies of a parachute cluster can be represented as a near contact.

In some classes of problems, obtaining a meaningful flow solution requires the contact to be represented as an actual contact. That includes heart valve flow analysis, where the valve leaflets have to have an actual contact, wing-clapping aerodynamics, where the clapping wings have to have an actual contact, and tire aerodynamics with road contact, where the tire and the road have to have an actual contact. That is computationally challenging because we want to both represent the contact without leaving a mesh protection opening and represent the boundary layers with an accuracy we get from moving-mesh methods. The ST-TC has overcome this challenge. We will provide an overview of the ST-TC mechanism in Sect. 2.

1.5 Synthesis of the ST-SI and ST-TC: ST-SI-TC

With the ST-SI-TC [81], the elements collapse and are reborn independent of the nodes representing a solid surface. High-resolution flow representation is enabled even when parts of the SI are coinciding with a solid surface. The ST-SI-TC enables effective management of contact location change and contact sliding. More on the ST-SI-TC and a partial list of the classes of problems computed with it can be found in [76].

1.6 ST-IGA

The terminology “ST-IGA” represents, depending on the context we have, ST discretization with IGA basis functions in space or time or both space and time. Superior accuracy of discretization with IGA basis functions in space [32, 74, 82, 83] is now textbook knowledge [69, 84]; it enables more accurate representation of the geometry and increased accuracy in the flow solution. Because that superior accuracy is reached with fewer control points, the effective element sizes will be larger. Consequently, we can take larger time steps in the computation without exceeding the Courant number we target for good accuracy.

Discretization with IGA basis functions in time is uniquely offered by the ST framework. That leads to a number of advantages. Deriving full benefit from using IGA basis functions in space (see the stability and accuracy analysis in [1]) is among those advantages. With that specific advantage, the ST framework will increase the Courant number we should target for good accuracy. We will get that increase even if we use linear basis functions in time. We will get even more increase if we use higher-order B-splines in time. We mention two other advantages in using IGA basis functions in time. (i) As can be seen in [10] for example, it brings higher accuracy in representing the motion of a solid surface, a mesh motion consistent with that surface motion, and better efficiency in representing the mesh motion and in remeshing. Representation of a circular path will be exact with quadratic NURBS basis functions, achievable with as little as two patches. (ii) The “ST-C” [85], where the letter “C” means “continuous,” is a method for extracting time-continuous data from the computed data. It can work as a data compression method in dealing with large volumes of data (see, for example, [76]).

More on the ST-IGA and a partial list of the classes of problems computed with it can be found in [76].

1.7 Synthesis of the ST-SI, ST-SI-TC, and ST-IGA: ST-SI-IGA and ST-SI-TC-IGA

The ST-SI-IGA [11] is the synthesis of the ST-SI and ST-IGA, and the ST-SI-TC-IGA is the synthesis of the ST-SI-TC and ST-IGA. They extend the ST-SI and ST-SI-TC, with all the good properties discussed in Sects. 1.3–1.5, to computations with the ST-IGA. (i) The classes of problems requiring computation with the ST-SI-TC quite often involve narrow and curved gaps between solid surfaces. Its synthesis with the ST-IGA keeps element density required for good accuracy at a reasonable level. That in turn keeps the computational cost at a reasonable level. (ii) Sometimes we need to accurately connect the parts of the solution obtained over two mesh zones with nonmatching meshes not because of rotating solid surfaces requiring slip between the two mesh zones but because we want to be free from the matching requirement. That gives us more flexibility in IGA mesh generation. That is helpful considering that IGA mesh generation for complex geometries is a challenge. Multiple SIs in a computation, some for the rotating solid surfaces and some for the mesh generation flexibility, can be expected in ST flow computations, such as the ventricle-valve-aorta computation in [19]. (iii) Computations with the ST-SI-TC also quite often involve elements with high aspect ratios. Its synthesis with the ST-IGA brings more robustness to computations with such elements. (iv) With the special feature introduced in [61], the ST-SI-TC-IGA has a built-in Reynolds-equation limit. Computations involving fluid films between solid surfaces coming into contact can be carried out without the need to use separately a Reynolds-equation model. The ST-SI-TC-IGA can emulate the Reynolds-equation model with a comparable solution quality and computational cost, working also in regions of the problem domain beyond the scope of the Reynolds-equation model. More on the ST-SI-IGA and ST-SI-TC-IGA and a partial list of the classes of problems computed with them can be found in [76].

1.8 Local length scales and stabilization parameters targeting IGA discretization

Local length scales, also called element lengths, and stabilization parameters play an essential role in the SUPS, VMS, and other stabilized methods. In addition to having good stability and accuracy properties, they should have good invariance properties, for example, invariance with respect to element node reordering. There are typically two local length scales. They are the advection and diffusion length scales, used in the expressions for the advective and diffusive limits of the stabilization parameter. The advection length scales in use have always been in the flow direction. The diffusion length scales in use, however, have mostly been just element-geometry-dependent. The reasons for also making that direction-dependent so that the spatial variation of the solution is taken into account were articulated in [86]. Those were the reasons behind the introduction of the diffusion length scale in the solution-gradient direction, which was given in 2001 [87] without explicitly stating the reasons and identified as the diffusion length scale in [6]. We note that the integrals in the ST-SI also involve element length, and that is in the direction normal to the SI.

Until 2017, local length scales and stabilization parameters used in computations with IGA basis functions were those originally developed for computations with finite element basis functions. Local length scales and stabilization parameters targeting IGA discretization were introduced in [13] and further developed in [14]. Despite being relatively new, they have been used in computation of many classes of problems. The overview given in [76] on the local length scales and stabilization parameters includes how they started, how they evolved over the years, how the versions targeting IGA discretization were derived, and the classes of problems computed with these versions.

1.9 EBMR and mesh relaxation and mesh moving based on fiber-reinforced hyperelasticity

This subsection is an overview of some of the most recent general-purpose methods for mesh quality improvement and for low-distortion mesh moving. The EBMR improves the mesh quality. It can be used for improving the quality after the initial creation of a mesh, or occasionally during a mesh motion, or as part of the mesh motion. When used in connection with the mesh motion, it is less disruptive and less time-consuming than restoring the mesh quality by remeshing. In the EBMR, the new mesh has the same set of nodes and elements as before, but some of the nodes are moved slightly to improve the quality of some of the elements. The motion is determined from the nonlinear-elasticity equations of large-deformation mechanics and a locally-defined zero-stress state (ZSS) [88, 89]. The mesh relaxation and low-distortion mesh moving methods based on fiber-reinforced hyperelasticity and optimized ZSS were introduced in [16]. Although the methods targeted IGA discretization, they can of course also be used for, as a special case of IGA discretization, finite element discretization with linear basis function in each direction. Element distortion is reduced by stiffening the element in multiple directions, which is achieved by placing fibers in those directions. The ZSS is optimized by seeking orthogonality of the parametric directions, by mesh relaxation, and by making the ZSS time-dependent as needed More on the EBMR and mesh relaxation and low-distortion mesh moving based on fiber-reinforced hyperelasticity can be found in [76].

1.10 IGA mesh generation in complex-geometry flow computations

It is known that IGA mesh generation for complex-geometry flow problems has been a challenge. To overcome that challenge, two general-purpose, complex-geometry NURBS mesh generation methods were introduced. They are the Complex-Geometry IGA Mesh Generation (CGIMG)¹ method [17] and NURBS Surface-to-Volume Guided Mesh Generation (NSVGMG) method [18]. Special-purpose, complex-geometry NURBS mesh generation methods introduced for specific classes of problems include those for the aerodynamics of tires with near-actual tire geometry and road contact [91] and for heart valve flow analysis [92]. More on the CGIMG and a partial list of the classes of problems computed with it can be found in [76]. More on the NSVGMG can be found in [18].

2 ST-TC mechanism

Fig. 1

ST-TC mechanism. As two internal solid surfaces (red and blue) come into contact, the inner spatial mesh between them, consisting of two elements, contracts and the outer meshes expand. While the two elements eventually collapse, the ST elements they belong to remain active. As the servant points (white) of the ST-TC master–servant point system detach from the boundaries, they become free points, and the free points reaching a solid surface become servants. The numbers are in reference to the parent mesh. (Color figure online)

In this section, we provide an overview of the ST-TC-mechanism. Figure 1 illustrates the mechanism. When two solid surfaces come into contact, the spatial mesh between them contracts and eventually the elements collapse over the contact area. As an element collapses, the ST element it belong to remains active. For example, in Fig. 1, the spatial element with Points 1 and 2 has zero volume at the first time level. However, it has nonzero volume at the second time level, and therefore the corresponding ST element has nonzero volume and is active. As another example from Fig. 1, the spatial element with Points 3 and 4 has nonzero volume at the third time level and zero volume at the fourth time level, and thus the corresponding ST element has nonzero volume.

The ST-TC relies on a master–servant point system that maintains the connectivity of the “parent” mesh when there is a contact. Consequently, the distribution model in the parallel-computing environment does not change during the computations with moving meshes. That maintains the computational efficiency at a practical level. When a free point reaches a solid surface, it becomes a servant of a master point. When a servant point detaches from a solid surface, it becomes a free point.

3 Computations reported in 2013

Two classes of computations were reported in 2013 [4] with actual-contact and high-resolution boundary-layer representations, in the same article where the ST-TC was introduced. Both classes of computations were for 2D models, with the ST-SUPS, ST-TC, and finite element discretization. One was for a pair of symmetrically-flapping surfaces (“Flapping pair”). The other one was for two pairs of symmetrically-flapping surfaces with coordinated opening/closing action, intended to resemble the left ventricle. We summarize here only the flapping-pair computation.

Fig. 2

Flapping pair. Surface positions at $t= 0$ (also 18), 3, 6, 9, 12, and 15 s

Table 1 Flapping pair. Mesh data and the number of time steps per flapping cycle

Fig. 3

Flapping pair. Mesh 1, Mesh 2, and Mesh 3. Inner-zone part around the flapping surfaces

The pair of surfaces, with zero thickness, were undergoing a prescribed sinusoidal flapping, as shown in Fig. 2, with a period of $T = 18.0\,\hbox {s}$. The projected width of the flapping surfaces along the horizontal axis was 1.0 m. The inflow velocity was 0.1 m/s. The Reynolds number based on these length and velocity scales was 6000. The boundary conditions were no-slip on the flapping surfaces, uniform horizontal velocity at the inflow, zero-stress at the outflow, and slip at the upper and lower boundaries. Three different meshes were tested, each with a different time-step size, to see the influence of increased refinement in space and time. The meshes had a structured inner zone made of 4-node quadrilateral elements and an unstructured outer zone made of 3-node triangular elements. Table 1 shows the mesh data and the number of time steps per flapping cycle. Figure 3 shows the inner-zone part around the flapping surfaces. During the flapping motion, only the inner mesh was moving, with a special-purpose, algebraic mesh moving method. For clarity in illustrating how the elements collapse and are reborn, we show the results only for Mesh 1. Figure 4 shows the mesh around the flapping surfaces at different instants during the flapping cycle. Figure 5 shows the velocity magnitude at the same instants as in Fig. 4.

Fig. 4

Flapping pair. Mesh 1 at the same instants as in Fig. 2

Fig. 5

Flapping pair. Velocity magnitude (m/s) for Mesh 1 at the same instants as in Fig. 2

4 Computations reported in 2014

Three classes of computations were reported in 2014 [62, 93] with actual-contact and high-resolution boundary-layer representations, all in 3D. The first one was flow analysis of an aortic valve model with coronary arteries [93], the second one was flow analysis of a mechanical aortic valve [93], and the third one was wing-clapping aerodynamics [62].

4.1 Aortic valve with coronary arteries

Fig. 6

Aortic valve with coronary arteries. Model geometry. Aorta, leaflets, sinuses, and coronary arteries

Fig. 7

Aortic valve with coronary arteries. Leaflets at $t/T= 0.0, 0.2, 0.4, 0.5, 0.6,\, \hbox {and}\, 0.8$

The model, shown in Figs. 6 and 7, has three leaflets, two small outlets, corresponding to the coronary arteries, and one main outlet, corresponding to the beginning of the aorta. The inlet and main outlet diameters are 23 mm. The two coronary artery diameters are 2.9 mm. The motion of the leaflets was prescribed in the computation, with a period of $T = {0.6}\,\hbox {s}$. The velocity specified at the inlet was corresponding to a time-dependent flow rate that had an average of $5000\, \hbox {mL/min}$. The boundary conditions were no-slip on the arterial walls and valves, traction-free at the outflow boundaries, and uniform velocity at the inflow boundary. The computations were with the ST-SUPS for the first two nonlinear iterations of each time step, ST-VMS for the last two, ST-TC, and finite element discretization. The mesh had structured inner zones around the leaflets and an unstructured outer zone. The inner zones were made of tetrahedral, pyramid-shaped, and wedge-shaped elements, and the outer zone was made of tetrahedral elements. The mesh had 1,417,910 nodes and 4,184,614 elements. There were 79 layers of nodes extruded from both the upper and lower surfaces of the leaflets. During the prescribed motion, only the inner zones were moving with a special-purpose mesh moving method. The time-step size was $3.00{\times }10^{-3}\,\hbox {s}$. Figure 8 shows the velocity magnitude. Figures 9 and 10 show the wall shear stress (WSS) on the leaflet surfaces.

Fig. 8

Aortic valve with coronary arteries. Velocity magnitude (m/s) at the same instants as in Fig. 7

Fig. 9

Aortic valve with coronary arteries. WSS (Pa) on the lower surface of the leaflets at the same instants as in Fig. 7

Fig. 10

Aortic valve with coronary arteries. WSS (Pa) on the upper surface of the leaflets at the same instants as in Fig. 7

4.2 Mechanical aortic valve

Fig. 11

Mechanical aortic valve. Components of the model. The pipe and a valve plate

Fig. 12

Mechanical aortic valve. Model dimensions and plate positions

Fig. 13

Mechanical aortic valve. The plates at $\theta = 0.0^{\circ }, 9.43^{\circ }, 18.87^{\circ }, 28.30^{\circ }, 37.73^{\circ }$, and $47.17^{\circ }$

Fig. 14

Mechanical aortic valve. The rotation angle over a cycle

Figure 11 shows the components of the model, a pipe and two rotating plates. Figures 12, 13, and 14 show the model dimensions, plate positions, and how they rotate, with the rotation angle $\theta$ defined in Fig. 12. The period was $T = 1.597\,\hbox {s}$. The flow passages between the plates and between each plate and the pipe wall are all blocked when the valve is closed at $\theta = {47.17}^{\circ }$. The average flow rate was 5736 mL/min. The boundary conditions were no-slip on the pipe walls and plates and uniform velocity at the inflow boundary. At the outflow boundary, the method described in [94] was used to control the reverse flow. The computations were with the ST-VMS, ST-TC, and finite element discretization. The mesh was generated all manually and consisted of 1,241,163 nodes and 1,217,264 tetrahedral, pyramid-shaped, and wedge-shaped elements. The time-step size was $3.175{\times }10^{3}\,\hbox {s}$. The mesh moving was with a special-purpose mesh moving method. Figure 15 shows the velocity magnitude. Figure 16 shows the WSS on the upper and lower surfaces of the plates.

Fig. 15

Mechanical aortic valve. Velocity magnitude (m/s) at $t/T = 0.0, 0.2, 0.4, 0.5, 0.7,\, \hbox {and}\, 0.9$

Fig. 16

Mechanical aortic valve. WSS (Pa) vectors on the upper and lower surfaces of the plates at the same instants as in Fig. 15

4.3 Wing clapping

Fig. 17

Wing clapping. Model. The wing thickness is zero

Fig. 18

Wing clapping. Wing configurations at $t/T = 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8,\, \hbox {and}\, 0.9$

Figure 17 shows the model. The wings were flapping symmetrically, with prescribed motion, as shown in Fig. 18, with a period of $T = 0.0365\,\hbox {s}$. The free-stream velocity was 4.5 m/s. The Reynolds number based on that and the average chord length was 5423. Three cases were computed in [62], with the angle of attack $\alpha ={0}^{\circ }, 5^{\circ }$, and $10^{\circ }$. We present the results from the computation with $\alpha ={10}^{\circ }$. The boundary conditions were no-slip on the wings and body, uniform horizontal velocity at the inflow boundary, zero-stress at the outflow boundary, and slip at the upper, lower, and side boundaries. The computations were with the ST-SUPS for the first two nonlinear iterations of each time step, ST-VMS for the last two, ST-TC, and finite element discretization. The meshes had structured inner zones around the wings and an unstructured outer zone. Both the structured and unstructured zones were made of 4-node tetrahedral elements. The mesh for $\alpha ={10}^{\circ }$ had 1,152,367 nodes and 6,896,762 elements. The time-step size was $4.51 {\times } 10^{-4}\,\hbox {s}$. During the flapping, only the mesh in the inner zones was moving, with a special-purpose mesh moving method. Figure 19 shows the flow patterns during a flapping cycle. Figures 20 and 21 show the WSS on the body and wing surfaces.

Fig. 19

Wing clapping. Helicity isosurfaces ($\pm \,5$ and $\pm \,10\,\hbox {m}^{2}/\hbox {s}^{2}$) at the same instants as in Fig. 18. Blue is for negative values, and red for positive, with the darker shades representing the lower magnitudes. (Color figure online)

Fig. 20

Wing clapping. WSS (Pa) on the body and the wing surfaces at the same instants as in Fig. 18. The upper surface of the upper wing (left side) and the lower surface of the lower wing (right side)

Fig. 21

Wing clapping. WSS (Pa) on the body and the wing surfaces at the same instants as in Fig. 18. The lower surface of the upper wing (left side) and the upper surface of the lower wing (right side). The white regions are the closed parts of the wings. (Color figure online)

5 Computations reported in 2016

Two classes of computations were reported in 2016 [12, 81] with actual-contact and high-resolution boundary-layer representations. The computations reported in [81] were for tire aerodynamics with road contact, using 2D and 3D models that included transverse grooves. The 2D model was a tire section. The computations reported in [12] were for an aortic valve, with a more realistic 3D model geometry than the one in Sect. 4. We exclude the 2D computation from our summaries in this section.

5.1 Tire aerodynamics with road contact and transverse grooves

Fig. 22

Tire aerodynamics with road contact and transverse grooves. Model cross-section, which is invariant. The model dimension in the axial direction is 27.5 cm. The deformation region is the circular sector with central angle $36^{\circ }$

Fig. 23

Tire aerodynamics with road contact and transverse grooves. Boundary conditions. The left (red) and bottom (yellow) boundaries represent the inflow and road, where the velocity is set to the linear speed at the tire periphery. The innermost (blue) circle is the tire surface, which is rotating. The larger (pink) circle is the SI. The bottom of the SI is coinciding with the road and the tire–road interface. The conditions at the right (green) and top (cyan) boundaries are zero-stress and slip, respectively. The condition at the two side boundaries is also slip. (Color figure online)

The geometry of the model, which has 12 treads and an invariant cross-section, and the deformation pattern are shown in Fig. 22. The rotation speed was corresponding to a linear speed of 100 km/h at the tire periphery. Figure 23 shows the computational domain and boundary conditions. The computation was with the ST-SUPS, ST-SI-TC, and finite element discretization. The mesh had 2,370,840 nodes and 5,415,467 elements. The inner mesh around the tire was made of hexahedral elements. It was generated manually, rotated with the tire, and connected with an SI to the outer mesh. The outer mesh had a first layer of hexahedral elements, generated manually, next a layer of pyramid elements, also generated manually, and after that tetrahedral elements, generated automatically. There were 2880 time steps per rotation, equivalent to a time-step size of $2.474{\times }10^{-5}\;\hbox {s}$. Figure 24 shows the velocity magnitude.

Fig. 24

Tire aerodynamics with road contact and transverse grooves. Velocity magnitude (km/h)

5.2 Aortic valve with a realistic model geometry

Fig. 25

Aortic valve with a realistic model geometry. Aorta, leaflets, and sinuses. Leaflet identification: Leaflet 1 (red), 2 (green), and 3 (blue). (Color figure online)

Fig. 26

Aortic valve with a realistic model geometry. Leaflet motion. Time variation of $\theta$ for each of the three leaflets

Fig. 27

Aortic valve with a realistic model geometry. Inflow velocity (two cycles)

The model, shown in Fig. 25, has three leaflets and one main outlet, corresponding to the beginning of the aorta. The leaflet motion was prescribed in the computation. The computations reported in [12] were for two cases of leaflet motion: symmetric and asymmetric. In the symmetric case, the leaflet motions had 3-fold rotational periodicity, and in the asymmetric case, they did not. We summarize the asymmetric case here. Figure 25 also shows the leaflet identification. The leaflet positions were defined by means of a pseudo-time parameter $\theta$, with the values 0 and 1 corresponding to the fully open and fully closed positions. The prescribed motion was given through $\theta$ as shown in Fig. 26. The boundary conditions were no-slip on the arterial walls and the leaflets, zero-stress at the outflow boundary, and uniform velocity at the inflow boundary, with a time variation as shown in Fig. 27. The cycle period was 0.712 s. The no-slip condition on the arterial walls was enforced weakly. The computation was with the ST-SUPS, ST-SI-TC-IGA, and IGA discretization. The mesh, made of quadratic NURBS elements, was created with five SIs (see Fig. 28). Three of them connect the mesh sectors containing the leaflets in the valve region. The other two, the top and bottom circular planes in the figure, connect the independently created meshes in the inlet and outlet regions to the valve region. The mesh had 84,534 control points and 54,000 elements. The time-step size was $4.00{\times }10^{-3}\,\hbox {s}$. The mesh moving was with a special-purpose mesh moving method. Figure 29 shows the flow patterns. The viewing angle is as we see the leaflets in Fig. 25. Figure 30 shows the WSS on the leaflet surfaces. The viewing angle is as we see the leaflets in Fig. 31.

Fig. 28

Aortic valve with a realistic model geometry. The valve and the five SIs

Fig. 29

Aortic valve with a realistic model geometry. Isosurfaces corresponding to a positive value of the second invariant of the velocity gradient tensor, colored by the velocity magnitude (m/s). The frames are for $t =$ 0.804, 0.984, 1.028, 1.072, 1.080, and 1.252 s

Fig. 30

Aortic valve with a realistic model geometry. WSS (Pa) on the upper (left) and lower (right) leaflet surfaces. The frames are for the same instants as in Fig. 29

Fig. 31

Aortic valve with a realistic model geometry. Viewing angle for reporting the WSS. The leaflet identification is the same as in Fig. 25

6 Computations reported in 2018

The computations reported in 2018 [91] with actual-contact and high-resolution boundary-layer representations were for or related to tire aerodynamics with near-actual tire geometry, road contact, and tire deformation. The computations involved a simple 2D model and a 3D model. The 2D model was just for verification purposes, had no grooves, but the inner mesh still rotated with the tire and some of its elements collapsed and were reborn. The 3D model had a near-actual tire geometry and included longitudinal and transverse grooves. The deformation pattern was realistic; it was provided by the tire company. We exclude the 2D computation from our summary in this section.

Fig. 32

Tire aerodynamics with near-actual tire geometry, road contact, and tire deformation. Tire model in its deformed configuration

Fig. 33

Tire aerodynamics with near-actual tire geometry, road contact, and tire deformation. Preliminary and refined meshes near the tire surface. The checkerboard pattern is for differentiating between the NURBS elements

The tire model, in its deformed configuration, is shown in Fig. 32. The diameter is 1.03 m and the width is 260 mm. There are three longitudinal grooves and 72 transverse grooves. The deformation of the tire was provided by the tire company, at five uniformly spaced instants of a $5^{\circ }$ rotation (see [91]). Cubic NURBS basis functions were used in the temporal representation of the deformation, which was extracted with the ST-C (see Sect. 1.6) from the data at those five instants. The tire rotation speed was corresponding to a speed of 103 km/h at the periphery of the undeformed tire, and the translation speed was $U_0 = 100\,\hbox {km/h}$. The boundary conditions were $U_0$ at the inflow and bottom boundaries, zero-stress at the outflow boundary, and slip at the top and side boundaries. The computations were with the ST-VMS, ST-SI-TC-IGA, and IGA discretization. Two meshes, made of quadratic NURBS elements, were used. The preliminary mesh had 690,144 control points and 376,560 elements, and the refined mesh had 4,149,720 control points and 2,921,552 elements. Figure 33 shows the two meshes. The main complexity of the mesh generation was near the tire surface and that part of the mesh was generated manually. Although the rest of the mesh could have been generated with the CGIMG (see Sect. 1.10), because it was not so difficult, that was also done manually. The number of time steps per rotation was 1440 when the preliminary mesh was used, and four times that when the refined mesh was used. Figures 34 and 35 show the flow velocity near the tire–road contact area. Figures 36 and 37 show the flow patterns near the tire–road contact area.

Fig. 34

Tire aerodynamics with near-actual tire geometry, road contact, and tire deformation. Velocity magnitude (km/h) at three uniformly spaced instants, from the computation with the preliminary mesh. Contact area. Data displayed on planes perpendicular to the tire axis

Fig. 35

Tire aerodynamics with near-actual tire geometry, road contact, and tire deformation. Velocity magnitude (km/h) at three uniformly spaced instants, from the computation with the refined mesh. Contact area. Data displayed on planes perpendicular to the tire axis

7 Computations reported in 2019

Three classes of computations were reported in 2019 [19, 92, 95] with actual-contact and high-resolution boundary-layer representations. The computations reported in [95] were for or related to tire aerodynamics with near-actual tire geometry, road contact, tire deformation, road roughness, and fluid film. They involved 2D and 3D models. The 2D model was a tire section with grooves and was for a study on testing the built-in Reynolds-equation limit of the ST-SI-TC-IGA and different options for dealing with the fluid trapped in the grooves. The 3D model, like the 3D model in Sect. 6, had a near-actual tire geometry and longitudinal and transverse grooves, with the tire deformation provided by the tire company. The computational model also included road roughness and fluid film, and the built-in Reynolds-equation limit of the ST-SI-TC-IGA was utilized. The computation reported in [92] was for flow analysis with a realistic 3D model of a bioprosthetic heart valve. The computation reported in [19] was for flow analysis with a realistic 3D model of the ventricle, heart valve, and aorta sequence. We exclude the 2D computation from our summaries here.

7.1 Tire aerodynamics with near-actual tire geometry, road contact, tire deformation, road roughness, and fluid film

Fig. 36

Tire aerodynamics with near-actual tire geometry, road contact, and tire deformation. Isosurfaces corresponding to a positive value of the second invariant of the velocity gradient tensor, colored by the velocity magnitude (km/h), viewed from the bottom, at three uniformly spaced instants. Computed with the preliminary mesh. The gray zones are the contact areas. (Color figure online)

Fig. 37

Tire aerodynamics with near-actual tire geometry, road contact, and tire deformation. Isosurfaces corresponding to a positive value of the second invariant of the velocity gradient tensor, colored by the velocity magnitude (km/h), viewed from the bottom, at three uniformly spaced instants. Computed with the refined mesh. The gray zones are the contact areas. (Color figure online)

Fig. 38

Tire aerodynamics with near-actual tire geometry, road contact, tire deformation, road roughness, and fluid film. Road-roughness representation. Dents containing the fluid films were represented not on the road surface (top) but on the tire surface (bottom)

Fig. 39

Tire aerodynamics with near-actual tire geometry, road contact, tire deformation, road roughness, and fluid film. Road surface with dents

In representing the road roughness, it was understood that the mesh resolution on the tire surface would be higher than the mesh resolution on the road surface. Therefore, the dents containing the fluid films were represented not on the road surface but on the tire surface. Figure 38 illustrates this concept. The paths needed for deisolating the dents were created on the tire surface. The tire model was the same as in Sect. 6. The tire rotation speed was corresponding to a speed of 100 km/h at the periphery of the undeformed tire, and the translation speed was $U_0 = 98.86\;\hbox {km/h}$. The boundary conditions were the same as in Sect. 6. The road roughness was added in a random fashion. Figure 39 shows the dents. The computation was with the ST-SUPS, ST-SI-TC-IGA, and IGA discretization. The mesh was made of quadratic NURBS elements, generated manually. It had 706,872 control points and 342,296 elements. Figure 40 shows the mesh. The number of time steps per rotation was 1440. Figure 41 shows the flow patterns. Figure 42 shows the fluid friction on the tire surface.

Fig. 40

Tire aerodynamics with near-actual tire geometry, road contact, tire deformation, road roughness, and fluid film. The mesh around the tire. The mesh with the yellow–gray element pattern was fixed. The mesh with the other-color element patterns was rotating with the tire. The elements in the red–gray pattern are the special elements left for the fluid film gaps in the contact area. The elements in the green–gray pattern are those in the grooves. The blue–gray pattern indicates the rest of the elements. (Color figure online)

Fig. 41

Tire aerodynamics with near-actual tire geometry, road contact, tire deformation, road roughness, and fluid film. Isosurfaces corresponding to a positive value of the second invariant of the velocity gradient tensor, colored by the velocity magnitude (km/h), viewed from the bottom, at three uniformly spaced instants. The gray indicates the full-contact zones. (Color figure online)

Fig. 42

Tire aerodynamics with near-actual tire geometry, road contact, tire deformation, road roughness, and fluid film. Fluid friction on the tire surface in terms of the positional-averaged WSS (Pa), at three uniformly spaced instants. The vectors indicate just the direction of the WSS

7.2 Bioprosthetic-heart-valve flow analysis

Fig. 43

BHV flow analysis. Model. Leaflets, metal frame, and sinuses

Fig. 44

BHV flow analysis. Inflow velocity

Fig. 45

BHV flow analysis. Valve quadratic NURBS surfaces at $t =$ 0.175, 0.275, 0.315, 0.495, 0.585, 0.590, 0.595, 0.605, 0.635, and 0.785 s

The bioprosthetic-heart-valve (BHV) and arterial-surface motion, with a cardiac cycle $T = 0.86\,\hbox {s}$ and with a time-interval size $5.00 {\times } 10^{-3}\,\hbox {s}$, were coming from the IMGA FSI computation reported in [96]. The BHV model is shown in Fig. 43. The boundary conditions were no-slip on the valve and artery surfaces, zero-stress at the outflow, and uniform velocity at the inflow. The flow rate at the inflow was a modified version of the one in the IMGA FSI computation. The modification was for making sure that during the closed-valve part of the cardiac cycle, the inflow-boundary flow rate and the closed-space volume-change rate matched. The inflow velocity corresponding to the modified flow rate is shown in Fig. 44. The computation was with the ST-SUPS, ST-SI-TC-IGA, and IGA discretization. Both the valve and arterial surfaces were represented by quadratic NURBS. Figure 45 shows the valve NURBS surfaces. The mesh, made of quadratic NURBS elements, had 777,001 control points and 572,556 elements. It had three SIs and three parts (“Part 1,” “Part2,” and “Part 3”). Figure 46 shows the SIs and the parts. Part 1 was facing the SIs. It contained the elements that collapsed and were reborn as the leaflets were moving. Its motion was created with a special-purpose mesh moving method that took into account the contact. Part 2 was not changing during the leaflet motion. Part 3 was the rest of the mesh, between Part 1, Part 2, and the artery surface. It was changing during the computation. Its motion was determined with the EBMR (see Sect. 1.9) and MJBS (see Sect. 1.2). The motion was driven by the motion of Part 1, the valve motion, and the motion of the artery surface. The time-step size was $5.00{\times } 10^{-3}\,\hbox {s}$. Figure 47 shows the flow patterns during the second cardiac cycle, and Fig. 48 shows the corresponding valve WSS.

Fig. 46

BHV flow analysis. The three SIs and three parts: Part 1 (red), Part 2 (blue), and half of Part 3 (green). (Color figure online)

Fig. 47

BHV flow analysis. Isosurfaces corresponding to a positive value of the second invariant of the velocity gradient tensor, colored by the velocity magnitude (m/s), $t=$ 1.135, 1.355, 1.450, 1.455, 1.465, and 1.645 s

Fig. 48

BHV flow analysis. WSS (Pa) at the same instants as in Fig. 47. One-third of the valve is transparent

7.3 Ventricle-valve-aorta flow analysis

Fig. 49

Ventricle-valve-aorta flow analysis. Model. LV (blue), valve with the leaflets (orange) and sinuses (transparent), and aorta (green). Front view and view along the valve axis. (Color figure online)

Figure 49 shows the model for the left ventricle (LV), valve, and aorta sequence. The model does not include the mitral valve; it ends where the LV meets the left atrium, and that is the inflow boundary. The boundary condition at the inflow was zero-velocity for the closed position of the mitral valve, and the Constrained-Flow-Profile (CFP) Traction (see [19]) for the open position. The other boundary conditions were zero-stress at the main aorta outlet, specified-velocity at the three smaller outlets, an no-slip on all solid surfaces. At the traction-specified inflow and outflow boundaries, the conditions were from $p_\textrm{in} = 10.5\,\hbox {kPa}$ and $p_\textrm{out} = 0$. At the velocity-specified outflow boundaries, the estimated flow rates were from the Murray’s law [97]. The flow rate distribution between the three outlets was proportional to $D^3$, with D being the average diameter of the outflow cross-section (see [98]). With a cardiac cycle of $T = 0.9\,\hbox {s}$, the LV shrinks and expands, and the valve opens and closes. The aorta was assumed rigid. The computation was with the ST-SUPS, ST-SI-TC-IGA, and IGA discretization. The time-step size selected was equivalent to 320 time steps per cardiac cycle. The mesh, made of quadratic NURBS elements, had 777,001 control points and 572,556 elements. It was generated separately for the LV, valve, and aorta. These mesh parts, not matching at the interfaces between them, were connected with SIs.

Fig. 50

Ventricle-valve-aorta flow analysis. LV when it is most expanded and most shrunk. The blue line indicates roughly where the LV structure was truncated to create the segment used in the flow computation. (Color figure online)

The LV geometry and motion were built, in a tabular form, from CT scans of the LV at one instant in the cardiac cycle and other anatomically realistic data. Figure 50 shows the LV when it is most expanded and most shrunk. Figure 50 also shows where the LV structure was truncated to create the segment used in the flow computation. The starting flow computation mesh was generated with the CGIMG (see Sect. 1.10). The cardiac-cycle representation of the mesh motion, conforming to the LV motion, was built by using cubic B-splines in time, the ST-C (see Sect. 1.6), and the procedure [99, 100] that assures cycle periodicity. The mesh moving method was the linear-elasticity mesh moving with MJBS (see Sect. 1.2), with the parameters given in Section 10.1.1 of [19]. The valve leaflet motion was produced from the shell structural mechanics computation in [101], which was performed with a time-periodic, spatially uniform pressure difference between the upper and lower surfaces of the leaflets. The time profile of the pressure was similar to the one used in [102]. The valve mesh motion was created, by a transformation, from the mesh motion of the valve ST-SI-TC-IGA computation in [101], where the leaflet contact was represented as an actual contact and the boundary layers were represented with high-resolution meshes. The aorta geometry, represented by cubic T-splines, was based on a set of CT scans different than the ones for the LV. Given that geometry, the mesh was generated with the CGIMG.

Fig. 51

Ventricle-valve-aorta flow analysis. Mesh composed of the LV, valve, and aorta segments, plus the special-purpose element (brown) for the CFP Traction. The lines mark the element boundaries. (Color figure online)

Fig. 52

Ventricle-valve-aorta flow analysis. Seven SIs (green) of the mesh and color-coded boundary conditions. Red: CFP Traction inflow, blue: traction-specified outflow, yellow: velocity-specified outflow. (Color figure online)

Figure 51 shows the mesh, composed of the LV, valve, and aorta segments, plus the special-purpose element for the CFP Traction. The special-purpose element has 27 basis functions and is connected to the rest of the mesh with an SI. Figure 52 shows the seven SIs of the mesh and color-coded boundary conditions. Two of the SIs are between the valve and the other two mesh segments, making the mesh generation easier. Three are between the valve mesh segments containing each leaflet, making the mesh moving easier. One is between the LV segment and the special-purpose element, facilitating the CFP Traction. One is between the two LV mesh segments, making the mesh generation easier. Figure 53 shows the mesh at six instants in the cardiac cycle. Figure 54 shows the flow patterns.

Fig. 53

Ventricle-valve-aorta flow analysis. Mesh at $t/T =$ 0.003, 0.250, 0.625, 0.684, 0.894, and 0.981

Fig. 54

Ventricle-valve-aorta flow analysis. Isosurfaces corresponding to a positive value of the second invariant of the velocity gradient tensor, colored by the velocity magnitude (m/s), at the same instants as in Fig. 53

8 Computations reported in 2022

Two classes of computations were reported in 2022 [18, 76] with actual-contact and high-resolution boundary-layer representations, both in 3D. The computation reported in [18] was for car and tire aerodynamics with near-actual car body and tire geometries, road contact, and tire deformation. The computations reported in [76] were for high-resolution multi-domain car and tire aerodynamics with near-actual car body and tire geometries, road contact, and tire deformation. The aerodynamic influence of the car body on the tires and influence of the front-tire wake on the rear-tire aerodynamics were accounted for in the framework of the Multidomain Method [103].

Fig. 55

Car and tire aerodynamics with near-actual car body and tire geometries, road contact, and tire deformation. Car body and tire models. Tire (blue), wheel (green), and disk rotor (red) models. (Color figure online)

8.1 Car and tire aerodynamics with near-actual car body and tire geometries, road contact, and tire deformation

The car body and tire models are shown in Fig. 55. The car is 4.65 m long, 1.81 m wide, and 1.16 m high. The tire diameter and width are 632 mm and 211 mm. A steady-state structural mechanics computation was carried out to obtain the tire deformation, resulting in the shape shown in Fig. 55. The rotation speed was corresponding to a linear speed of 100 km/h at the undeformed tire periphery, with a car speed of $U_0 = 99.39\;\hbox {km/h}$. The tire rotation period, which we will use in reporting the results, is represented by the symbol $T_{\textrm{TR}}$. The boundary conditions were $U_0$ at the inflow and bottom boundaries, zero-stress at the outflow boundary, slip at the top and side boundaries, and no-slip on the car body, tires, wheels, and disk rotors.

The computation was with the ST-VMS, ST-SI-TC-IGA, and IGA discretization. The mesh, made of quadratic NURBS elements, had 1,193,580 control points and 707,372 elements. It was generated over five domain parts, separately, and the mesh parts were connected with SIs. Both the NSVGMG (see Sect. 1.10) and special-purpose mesh generation methods were used. The mesh relaxation method based on fiber-reinforced hyperelasticity (see Sect. 1.9) was also used. Figure 56 shows the mesh. Figure 57 shows the mesh around the tire and wheel. Although the tire model includes only longitudinal grooves, the mesh around the tire was still rotating to test the ST-TC mechanism. That was the only part of the mesh that was moving and deforming. The time-step size was $T_{\textrm{TR}}/144$. Time is measured from a point reached after a sufficiently long computation at full Reynolds number. The computation was carried out until $t = 5 T_{\textrm{TR}}$. Figure 58 shows the flow patterns during the last $T_{\textrm{TR}}/36$. Figures 59 and 60 show the flow patterns around the tire–road contact areas during the same period.

Fig. 56

Car and tire aerodynamics with near-actual car body and tire geometries, road contact, and tire deformation. Mesh. Entire domain and cross-sections along the width and flow directions. Dark gray indicates the car body. The checkerboard pattern is for differentiating between the elements. The colors are for differentiating between the patches. (Color figure online)

Fig. 57

Car and tire aerodynamics with near-actual car body and tire geometries, road contact, and tire deformation. Mesh around the tire and wheel. Dark gray indicates the tire and light gray indicates the wheel. The checkerboard pattern is for differentiating between the elements. The colors are for differentiating between the patches. (Color figure online)

Fig. 58

Car and tire aerodynamics with near-actual car body and tire geometries, road contact, and tire deformation. Isosurfaces corresponding to a positive value of the second invariant of the velocity gradient tensor, colored by the velocity magnitude (km/h), at 3 uniformly spaced instants during the last $T_{\textrm{TR}}/36$

Fig. 59

Car and tire aerodynamics with near-actual car body and tire geometries, road contact, and tire deformation. Isosurfaces corresponding to a positive value of the second invariant of the velocity gradient tensor, colored by the velocity magnitude (km/h), viewed from the bottom, at the same instants as in Fig. 58. Front right and left tires. The free-stream flow is from left to right. The dark gray zones are the contact areas. The top frames show the tires viewed, with the arrows indicating the viewing direction. (Color figure online)

Fig. 60

Car and tire aerodynamics with near-actual car body and tire geometries, road contact, and tire deformation. Isosurfaces corresponding to a positive value of the second invariant of the velocity gradient tensor, colored by the velocity magnitude (km/h), viewed from the bottom, at the same instants as in Fig. 58. Rear right and left tires. The free-stream flow is from left to right. The dark gray zones are the contact areas. The top frames show the tires viewed, with the arrows indicating the viewing direction. (Color figure online)

8.2 High-resolution multi-domain car and tire aerodynamics with near-actual car body and tire geometries, road contact, and tire deformation

Fig. 61

High-resolution multi-domain car and tire aerodynamics with near-actual car body and tire geometries, road contact, and tire deformation. Sequence of three local domains. The dimensions are in m. Blue indicates the body and red indicates the tires, wheels, and disk rotors. (Color figure online)

Fig. 62

High-resolution multi-domain car and tire aerodynamics with near-actual car body and tire geometries, road contact, and tire deformation. Boundary conditions for the three local domains. Blue: car body, red: ground, green: stress boundary condition, dark gray: tire, yellow: wheel, and transparent: prescribed-velocity boundaries. The two planes with orange frames are the inflow planes for those subdomains. (Color figure online)

The car body and tire models are the same as in Sect. 8.1. Taking the data from the computation described in Sect. 8.1 as the global-computation data, high-resolution computations for the left set of tires were carried out over the sequence of three local domains shown in Fig. 61. Figure 62 shows the boundary conditions for the three local domains. The prescribed-velocity boundary conditions were extracted from the global solution, except for the inflow planes of the second and third local domains. The velocities at those planes came from the solutions over the first and second local domains. The stress boundary conditions were extracted from the global solution. We note that the global data that was used as the extraction source was stored with the ST-C (see Sect. 1.6), using cubic B-splines in time. The time-step size for the cubic B-spline representation was $T_{\textrm{TR}}/48$, which was three times the time-step size used in the global computation.

Fig. 63

High-resolution multi-domain car and tire aerodynamics with near-actual car body and tire geometries, road contact, and tire deformation. Meshes over the three local domains. The checkerboard pattern is for differentiating between the elements. The colors are for differentiating between the patches

Fig. 64

High-resolution multi-domain car and tire aerodynamics with near-actual car body and tire geometries, road contact, and tire deformation. Local mesh around the tire and wheel. Dark gray indicates the tire and light gray indicates the wheel. The checkerboard pattern is for differentiating between the elements. The colors are for differentiating between the patches. (Color figure online)

Fig. 65

High-resolution multi-domain car and tire aerodynamics with near-actual car body and tire geometries, road contact, and tire deformation. Isosurfaces corresponding to a positive value of the second invariant of the velocity gradient tensor, colored by the velocity magnitude (km/h), at 3 uniformly spaced instants during the last $T_{\textrm{TR}}/36$

Fig. 66

High-resolution multi-domain car and tire aerodynamics with near-actual car body and tire geometries, road contact, and tire deformation. Isosurfaces corresponding to a positive value of the second invariant of the velocity gradient tensor, colored by the velocity magnitude (km/h), viewed from below the tires, at the same instants as in Fig. 65. The free-stream flow is from left to right. The dark gray zones are the contact areas. The top frames show the tires viewed, with the arrows indicating the viewing direction. (Color figure online)

The computations were with the ST-VMS, ST-SI-TC-IGA, and IGA discretization. The mesh over the first local domain had 1,973,195 control points and 1,487,632 elements. The mesh over the second local domain had 941,108 control points and 845,000 elements. The mesh over the third local domain had 2,079,835 control points and 1,576,232 elements. The meshes over the first and third local domains were generated over four parts each, and the meshes were connected with SIs. Both the CGIMG (see Sect. 1.10) and special-purpose mesh generation methods were used. The mesh relaxation method based on fiber-reinforced hyperelasticity (see Sect. 1.9) was also used. Figure 63 shows the meshes over the three local domains. Figure 64 shows the local mesh around the tire and wheel. Again, the mesh around the tire was rotating to test the ST-TC mechanism. The time-step size was $T_{\textrm{TR}}/432$. The computations over the first, second, and third local domains were started at $t=0$, $t = 0.59 T_{\textrm{TR}}$, and $t = 1.38 T_{\textrm{TR}}$. The initial conditions were from the global computation at those times. The computations were continued until $t = 5 T_{\textrm{TR}}$. We note that $0.59 T_{\textrm{TR}}$ and $1.38 T_{\textrm{TR}}$ are the approximate durations for the flow information advected from the inflow plane of the first local domain at approximately 100 km/h to reach the inflow planes of the second and third local domains. Figure 65 shows the flow patterns during the last $T_{\textrm{TR}}/36$. Figure 66 shows the flow patterns around the tire–road contact areas during the same period.

9 Concluding remarks

There was a perception that in computation of flow problems where moving solid surfaces come into contact, high-resolution boundary-layer representation by using a moving-mesh method and actual-contact representation without leaving a mesh protection opening between the solid surfaces were mutually exclusive objectives. That changed in 2013 with the introduction of the ST-TC. This is a method that enables, in the ST framework, moving-mesh computation with actual-contact and high-resolution boundary-layer representations. In the ST-TC, elements collapse or are reborn as needed, and that is attainable in the ST framework while retaining the computational efficiency at a practical level. The ST-TC now has a 10-year history of achieving the two objectives that were long seen as mutually exclusive.

With the ST-TC and other ST computational methods introduced before and after, it has been possible to address many of the challenges encountered in conducting flow analysis with boundary layer and contact representation, in the presence of additional intricacies such as geometric complexity, IGA discretization, and rotation or deformation of the solid surfaces. The methods used with the ST-TC in addressing the challenges include the ST-SUPS and ST-VMS, which serve as the core methods, and the special methods introduced in connection with the core methods. The special methods include the ST-SI, ST-IGA, synthesized ST methods like the ST-SI-TC-IGA, complex-geometry IGA mesh generation methods like the CGIMG, and mesh relaxation and mesh moving based on fiber-reinforced hyperelasticity. To help widen awareness of these ST methods and what they can do, we provided an overview of the methods and the computations performed over the 10-year history of the ST-TC. We summarized, in chronological order, 12 computations with actual-contact and high-resolution boundary-layer representations, one in 2D from 2013 and 11 in 3D from the subsequent years. The computations summarized included car and tire aerodynamics with road contact and tire deformation and ventricle-valve-aorta flow with leaflet–leaflet contact. The summaries demonstrate the level this class of ST computational flow analysis has reached over its 10-year history.