High-resolution multi-domain space–time isogeometric analysis of car and tire aerodynamics with road contact and tire deformation and rotation

We are presenting high-resolution space–time (ST) isogeometric analysis of car and tire aerodynamics with near-actual tire geometry, road contact, and tire deformation and rotation. The focus in the high-resolution computation is on the tire aerodynamics. The high resolution is not only in space but also in time. The influence of the aerodynamics of the car body comes, in the framework of the Multidomain Method (MDM), from the global computation with near-actual car body and tire geometries, carried out earlier with a reasonable mesh resolution. The high-resolution local computation, carried out for the left set of tires, takes place in a nested MDM sequence over three subdomains. The first subdomain contains the front tire. The second subdomain, with the inflow velocity from the first subdomain, is for the front-tire wake flow. The third subdomain, with the inflow velocity from the second subdomain, contains the rear tire. All other boundary conditions for the three subdomains are extracted from the global computation. The full computational framework is made of the ST Variational Multiscale (ST-VMS) method, ST Slip Interface (ST-SI) and ST Topology Change (ST-TC) methods, ST Isogeometric Analysis (ST-IGA), integrated combinations of these ST methods, element-based mesh relaxation (EBMR), methods for calculating the stabilization parameters and related element lengths targeting IGA discretization, Complex-Geometry IGA Mesh Generation (CGIMG) method, MDM, and the “ST-C” data compression. Except for the last three, these methods were used also in the global computation, and they are playing the same role in the local computation. The ST-TC, for example, as in the global computation, is making the ST moving-mesh computation possible even with contact between the tire and the road, thus enabling high-resolution flow representation near the tire. The CGIMG is making the IGA mesh generation for the complex geometries less arduous. The MDM is reducing the computational cost by focusing the high-resolution locally to where it is needed and also by breaking the local computation into its consecutive portions. The ST-C data compression is making the storage of the data from the global computation less burdensome. The car and tire aerodynamics computation we present shows the effectiveness of the high-resolution computational analysis framework we have built for this class of problems.

These methods were addressing in [1] various computational challenges encountered in the car and tire aerodynamics. The ST context was bringing higher-order accuracy (see [2,3]). The VMS feature of the ST-VMS was bringing better representation of the multiscale, turbulent flow patterns. The ST context was also bringing what comes with moving-mesh methods, and that is the high-resolution flow computation near the moving solid surfaces. The ST-SI was making the moving-mesh computation possible with the tire rotating. The inner mesh around the tire rotates with the tire. The SI between the inner and outer meshes accurately connects the inner and outer parts of the solution. The SI was also, in some places, just connecting the solution parts represented over mesh zones with nonmatching meshes at the interface between the zones and the interface was being treated as an SI. The ST-TC was making the moving-mesh computation possible even with the TC created by the contact between the tire and the road. The contact was being represented without giving up on high-resolution flow representation near the tire. The ST-SI-TC [12,17], which is the integrated combination of the ST-SI and ST-TC, was making it possible to have high-resolution representation near the tire and road surfaces even when some parts of the SI were coinciding with those surfaces. It was also facilitating contact location change as well as contact sliding. The ST-IGA was bringing higher accuracy in representing the tire geometry and the flow solution. The ST-SI-IGA [10] and ST-SI-TC-IGA, which are the integrated combinations of the ST-IGA with the ST-SI and ST-SI-TC, were bringing that accuracy even when the ST-SI or ST-SI-TC was needed to be used. These integrated combinations were also making it possible to have a reasonable element density in the tire grooves and near the contact areas and therefore the computational cost was being kept at a reasonable level. The EBMR was improving the quality of the meshes generated, and the NSVGMG was making the NURBS mesh generation for the complex car and tire geometries less challenging.
The car and tire ST isogeometric analysis in [1] was carried with a reasonable mesh resolution. In this article, we are presenting a high-resolution ST isogeometric analysis. The focus is on the tire aerodynamics. The high resolution is in both space and time. The influence of the aerodynamics of the car body comes, in the framework of the Multidomain Method (MDM) [18], from the global computation in [1]. The high-resolution local computation, carried out for the left set of tires, takes place in a nested MDM sequence over indicates the tires, wheels, and disk rotors. Bottom: the NURBS mesh around the tire, wheel, and disk rotor. The checkerboard pattern is for differentiating between the elements and the colors are for differentiating between the patches three subdomains (see Fig. 1). The first subdomain contains the front tire. The second subdomain, with the inflow velocity from the first subdomain, is for the front-tire wake flow. The third subdomain, with the inflow velocity from the second subdomain, contains the rear tire. All other boundary conditions for the three subdomains are extracted from the global computation.
All the methods we listed in the first paragraph as the components of the global computational framework, except for the last one, are also in the local computational framework. They are playing in the local computation the roles they played in the global computation. Beyond that, the Complex-Geometry IGA Mesh Generation (CGIMG) method [19,20], MDM, and the "ST-C" [21] data compression are in the local computational framework. The CGIMG is making the IGA mesh generation for the complex geometries less arduous. The MDM is reducing the computational cost by focusing the high-resolution locally to where it is needed and also by breaking the local computation into its consecutive portions. The ST-C data compression is making the storage of the data from the global computation less burdensome. and moving boundaries and interfaces (MBI). It originated from and subsumes its precursor the Deforming-Spatial-Domain/Stabilized ST (DSD/SST) method [24][25][26]. The DSD/SST is mostly called "ST-SUPS," with the abbreviation "SUPS" denoting the stabilization components SUPG and PSPG, which stand for the Streamline-Upwind/Petrov-Galerkin [27] and Pressure-Stabilizing/Petrov-Galerkin [24]. The ST-SUPS, in broader interpretation of the terminology, includes the stabilization component coming from Least-Squares on the Incompressibility Constraint (LSIC), as the DSD/SST in its form in [25] did. The VMS components of the ST-VMS are from the residual-based VMS (RBVMS) method [28][29][30][31]. The increased accuracy associated with the ST framework (see [2,3,32]) makes the ST-SUPS and ST-VMS appealing also in flow computations without MBI. Furthermore, the framework, naturally, makes it possible to use IGA basis functions also in time [32].
The arbitrary Lagrangian-Eulerian (ALE) moving-mesh framework is older, though its use in 3D finite element flow computations with modern stabilized methods like the SUPG is somewhat newer compared to the ST-SUPS. The ram-air parachute FSI analysis in [33] was one of the earliest computations with the ALE-SUPS. In the category of moving-mesh methods with the stabilization components coming from the RBVMS, however, the ALE-VMS method [34-37] precedes the ST-VMS. The ST-SUPS, ALE-SUPS, ALE-VMS, and ST-VMS, as methods, have much in common, and so do the classes of problems computed with them since their inception.
The ST-SUPS, ALE-SUPS, ALE-VMS, and ST-VMS, like all moving-mesh methods, need to be complemented with mesh update methods in FSI and MBI computations. The mesh update most of the time consists of moving the mesh to accommodate the motion of the boundaries and interfaces and to control the mesh resolution near solid surfaces that are moving, and remeshing if the element distortion exceeds an acceptable level. We expect two things from a good mesh moving method: to reduce the need for remeshing and to give high priority to maintaining element quality near solid surfaces where accurate representation of the boundary layers matters. Since the inception of the ST-SUPS in 1990, a large number of special-and general-purpose mesh moving methods have been developed for computations with the ST-SUPS and ST-VMS. Some of them have also been used in computations with the ALE-SUPS and ALE-VMS. A recent article [141] on mesh moving methods provides an overview. The general-purpose methods include, as the first one, the linear-elasticity mesh moving method with mesh-Jacobian-based stiffening [158,161] introduced in 1992, and, as the most recent ones, the EBMR (see Sect. 1.7), where the mesh motion is determined by using the large-deformation mechanics equations and an element-based zero-stress state (ZSS), mesh relaxation and mesh moving methods [162] based on fiber-reinforced hyperelasticity and optimized ZSS, and the linear-elasticity mesh moving method with no cycleto-cycle accumulated distortion [155,163].

ST-SI
This subsection, included for completeness, is mostly from [22]. The "sliding interface" method was introduced in [164] in the context of the ALE-VMS. We will call that "ALE-SI." The ST-SI is the ST version of that. The acronym "SI" is implying both "sliding" and "slip," because, independent of the choice of the word, the SI serves the same purpose in both the ALE-SI and ST-SI, just in two different contexts. With the ALE-SI and ST-SI, the ST-SUPS, ALE-SUPS, ALE-VMS, and ST-VMS can be used even in the presence of a rotating solid surface, such as a car tire, benefiting from what comes with the moving-mesh methods. The mesh around the rotating solid surface and inside the SI, with higher refinement near the solid surface, rotates with it, sustaining the highresolution boundary layer representation. The mesh outside the SI does not rotate with the solid surface but could still be moving for some other reason. Accurately connecting the two sides of the solution is achieved by adding to the core formulation (ST-SUPS, ALE-SUPS, ALE-VMS, or ST-VMS) a set of integrals over the SI. The velocity and stress compatibility between the two sides of the SI is accounted for by the added integrals. Both the ALE-SI and ST-SI were originally formulated in the context of incompressible flows. Other ST-SI versions were introduced for purposes different than but comparable to the original purpose.
Together with the ST-SI version with the "fluid-fluid SI," a version with "fluid-solid SI" was introduced in [5]. One side of the SI becomes a solid surface. The integrals over the SI facilitate the weak enforcement of the Dirichlet boundary conditions for the fluid. This version is basically the ST construction of the weak-Dirichlet-condition method [165]. The ST-SI version with the "porosity SI" was also introduced in [5]. It has the SI between a thin porous structure and the fluid on its two sides and makes how the porosity is handled consistent with how the fluid-fluid and fluid-solid SIs are handled. The versions for the coupled incompressible-flow and thermal-transport equations, with both the fluid-fluid and fluid-solid SIs, were introduced in [6]. With those versions, thermo-fluid boundary layers near rotating solid surfaces can also have high-resolution representation and the weak enforcement of the Dirichlet boundary conditions is extended to thermo-fluid problems. The ST-SI versions introduced in [146] are the compressible-flow counterparts of the fluidfluid, fluid-solid, and porosity SIs. They were connected to the compressible-flow ST SUPG method [166] and the compressible-flow porosity models that were also introduced in [146].

ST-TC
This subsection, included for completeness, is mostly from [23]. In FSI and MBI problems, contact between moving solid surfaces can be an actual contact or "near contact." In the near contact, there is still some little separation between the solid surfaces and therefore there is no topology change in the fluid mechanics domain. In such cases, the nearness is close enough for obtaining physically reasonable solutions from the flow computation; in other words, computing with that little separation is basically good enough for solving the problem. With a good mesh moving method, no element needs to collapse and good boundary layer resolutions can be retained. Several classes of flow problems were computed with the ST-SUPS and ST-VMS under near-contact conditions with sufficient accuracy. Examples can be found in the references mentioned in [7].
In some classes of flow problems, however, it is essential to represent the contact as an actual contact. For example, in heart valve flow analysis, for obvious reasons, the contact between the valve leaflets needs to be represented as an actual contact, without any separation. As another example, in wing clapping aerodynamics of insects, the contact between the upper and lower wings needs to be an actual contact. The ST-TC was introduced to make ST moving-mesh computations possible even in flow problems that involve an actual contact. With that, we can both represent an actual contact and retain the good boundary layer resolution. Elements collapse as needed, which is viable in the ST context. The connectivity of the "parent" mesh, however, does not change during the process of element collapse or rebirth, and therefore the computational efficiency is not harmed.

ST-SI-TC
This subsection, included for completeness, is mostly from [128]. Some classes of problems will need both the ST-SI and ST-TC. In the ST-SI version with the fluid-fluid SI, we need elements on both sides of the SI. The SI might be between the solid surfaces coming into contact, or in a more general context, might be merging with a fluid-solid interface. The elements between the solid surface and the coinciding SI segment collapse in the ST-TC process, and the SI segment switches from the fluid-fluid SI to the fluid-solid SI, creating an SI that is a mixture of the two SI types. With the ST-SI-TC, i) the element collapse and rebirth process gains independence from the solid-surface nodes, ii) we can have high-resolution boundary layer representation near the fluidsolid interfaces even when an SI segment coincides with the solid surface, and iii) we can manage, in an effective fashion, contact location change and contact sliding.

ST-IGA
This subsection, included for completeness, is mostly from [22]. The IGA basis functions in space brought major accuracy increases in fluid and solid mechanics computations [34,93,164,168]. That made IGA basis functions appealing for the ST-SUPS and ST-VMS computations and led to the introduction of the ST-IGA, at the same time the ST-VMS was introduced. It is IGA discretization in the ST framework. The terminology "ST-IGA" implies, depending on the context, discretization with IGA basis functions in space or time or both. The test computations reported in [2], which were in 2D, were for flow past an airfoil and for pure advection of a scalar. The flow computation was with IGA basis functions in space, and the advection computations with IGA basis functions in both space and time, accompanied by a stability and accuracy analysis for the pure advection equation. The advection computations and stability and accuracy analysis showed what can naturally be expected, and that is, higher-order basis functions in space will deliver more if they are used together with higher-order basis functions in time. Keeping in mind that the increased accuracy the ST-IGA with IGA basis functions in space brings is attained with fewer control points, the effective element sizes will be larger. With that, larger time steps can be taken while still keeping the Courant number at or below the levels we target for good accuracy.
Using IGA basis functions in time is uniquely offered by the ST framework, and partly because of that the effort was focused on that track in the early years of the ST-IGA computations [2,3,9]. Taking advantage of that opportunity 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. The ST/NURBS Mesh Update Method (STNMUM) [9,123] was built around these positive attributes of the ST-IGA. Without loss of generality, the context of rotating solid surfaces can be used for highlighting the advantages of the STNMUM. Representation of the circular path is exact with quadratic NURBS basis functions, and that can achieved with as little as two patches. To have speeds that are invariant along the circular path, a constant angular speed can be prescribed with the aid of a secondary mapping (see [2,3,9,36]). The ST-C is another example of the good things that come out of the ST-IGA with IGA basis functions in time. The letter "C" in "ST-C" means "continuous." This is a method for extracting time-continuous data from the computed data, and it can work as a data compression method in dealing with large data volumes [4,6,21,[54][55][56][57][125][126][127]129,130,132]. The classes of problems computed by using the ST-IGA with IGA basis functions in time include wind turbines [5,52,54,55,123,124,128], turbomachinery [10,20,54,55,129-131], flapping-wing aerodynamics [7,9,36,[138][139][140][141], spacecraft cover separation aerodynamics [142], and higher-order temporal IGA discretization [32].

ST-SI-IGA and ST-SI-TC-IGA
This subsection, included for completeness, is mostly from Beyond that, the ST-SI-IGA addresses the mesh generation challenge in IGA discretization. This is accomplished with an SI that does not have a slip between the two sides. The SI just connects the parts of the solution obtained over two IGA mesh zones with nonmatching meshes at the SI between the zones. Because we are no longer constrained by a matching requirement, in computation of flow problems with complex geometries, the IGA discretization becomes more practical. In IGA-based computations with a thin porous structure embedded in the flow field, the ST-SI-IGA mechanism is essentially the same as what we described in Sect. 1.2 for the ST-SI. In some cases, the rotating solid surface has grooves or creates narrow spaces or the thin porous structure has gaps and slits. In computation of such flow problems, the ST-SI-IGA makes it possible to keep the element density, and consequently the computational cost, at an acceptable level. That makes computations even with such geometric complexities practical.
The ST-SI-TC-IGA, in flow computations with contact between moving solid surfaces, makes it possible to keep the element density in the narrow spaces close to the contact region at an acceptable level. While the solid surfaces come into contact, prior to the collapse of the elements between a solid surface and SI, we may have curved and complex boundaries and narrow spaces. These would need highaspect-ratio elements. The ST-SI-TC-IGA makes it possible to compute under such adverse conditions with an acceptable level of computational cost. With the enhancements introduced in [137], the ST-SI-TC-IGA acquired a built-in Reynolds-equation limit. With that, when the solid surfaces coming into contact have fluid films between them, we do not need to use separately a Reynolds-equation model in those regions. The ST-SI-TC-IGA can handle that with comparable solution quality and computational cost and also work in the other parts of the flow domain where the Reynolds-equation model would not work.

EBMR, ZSS, and fiber-reinforced hyperelasticity
This subsection, included for completeness, is mostly from [141,162]. It provides a short description of the EBMR, ZSS, and fiber-reinforced hyperelasticity concepts mentioned in Sect. 1.1.

EBMR
The EBMR restores the mesh integrity lost during the mesh motion, but does that without remeshing. The loss of mesh integrity in regions that we care more about does not happen so often because of the advanced mesh moving methods used with the ST-SUPS and ST-VMS, but could happen in computations with a high level of complexity. The FSI computations reported in [14, [186][187][188][189][190] for spacecraft parachute clusters, for example, had that type of complexity. As proposed in [14], when faced with a loss of mesh integrity, the EMBR can be used for relaxing the mesh without changing it at the fluid-structure interface, and the mesh integrity is restored to some extent. This is, as commented in [14], a less disruptive and less time-consuming alternative to remeshing. The EBMR does not change the number of elements or nodes, but just moves slightly some of the nodes to improve the quality of the elements in need. The motion is determined from the nonlinear-elasticity equations of large-deformation mechanics and an element-based ZSS (EBZSS). The EBZSS is essentially a shape generated for each element, and by design, the undeformed shape is made of "target elements" and is the shape we want to obtain by solving the nonlinear-elasticity equations. There are a number of options for building the target element shapes and they are given in [14]. The EMBR was successfully used in FSI computation of spacecraft parachute clusters (see [14]).

Locally-defined ZSS
The locally-defined ZSS started as an arterial ZSS estimation [105,[170][171][172][173]191,192]. It was formulated first as the EBZSS in the context of finite element discretization [191,192], next as the EBZSS in the context of isogeometric dis-cretization [170,171], and next as the integration-point-based ZSS (IPBZSS) in the context of isogeometric discretization [172,173]. In the EBZSS the ZSS is defined for each element by a set of positions. Nodes (or control points) from different elements mapping to the same node in the mesh do not have to have the same ZSS-defining positions. In the reference configuration, however, all elements are connected by nodes and the displacements are measured from that configuration. Formulating the structural mechanics problem in this fashion is referred to as "element-based total Lagrangian" (EBTL) in [191]. The EBTL is a key part of the EBMR [14]. In the IPBZSS, the way the EBZSS is defined is extended to its integration-point counterpart, with the ZSS represented in terms of the metric tensor. The IPBZSS has more parameters than the EBZSS. Therefore, while the conversion from the EBZSS representation to the IPBZSS representation is straightforward and will be exact, the reverse conversion, in general, will not be exact. Formulating the structural mechanics problem in this fashion is referred to as "integration-point-based total Lagrangian" (IPBTL) in [162].

Mesh relaxation and mesh moving based on fiber-reinforced hyperelasticity and optimized ZSS
The mesh relaxation and mesh moving methods based on fiber-reinforced hyperelasticity and optimized ZSS were introduced targeting IGA discretization. They of course also suitable for use with finite element discretization as a special case of IGA discretization. Element distortion during the mesh deformation is reduced by stiffening the element in multiple directions with the fibers placed in those directions. The ZSS is optimized by seeking, with mesh relaxation, orthogonality of the parametric directions and by making the ZSS time-dependent as needed. With the mesh relaxation, we improve the quality of the mesh after its initial creation and have an equilibrium state with the optimized ZSS, boundary conditions, and constitutive law. Preceding the use of the mesh relaxation in the global car and tire aerodynamics computation, the NURBS mesh used in the computational flow analysis reported in [124] for a tsunami-shelter vertical-axis wind turbine was obtained with the mesh relaxation method.

Stabilization parameters and local length scales targeting IGA discretization
This subsection, included for completeness, is mostly from  [194] for stabilization parameters were for separate τ SUPG and τ PSPG , and test computations with them were successful. Leaving that aside, a single parameter, "τ SUPS " [2,36], is used instead of two. The expressions for the stabilization parameters involve, among other constituents, local lengths scales, also called "element lengths." Precursors of the element lengths and stabilization parameters used with the SUPG, PSPG, ST-SUPS, ALE-SUPS, RBVMS, ALE-VMS, and ST-VMS today can be found in [27,[195][196][197]. The work on designing good element lengths and stabilization parameters continued as a significant part of the research on residual-based stabilized methods. That generated a good number favorite expressions for the element lengths and stabilization parameters (see, for example, [4,9,25,26,123,194,[198][199][200][201][202][203]), and until late 2017 they were all intended for finite element discretization. For about a decade, the element length was a local length scale in the flow direction, which is an advection length scale. A second local length scale, in the solution-gradient direction, was introduced in [200]) and was identified as the diffusion length scale in [25]. The element length also appears in some of the integrals over the SI in the ALE-SI and ST-SI, and the relevant direction in that case is the SI normal. Element lengths and stabilization parameters in the ST context were introduced in [25, 200], those specific to the VMS stabilization in [4], and those for coupled incompressible-flow and thermal-transport equations in [4]. Direction-dependent element lengths that have node-numbering invariance also for simplex elements were introduced in [203]. All these local length scales and stabilization parameters were, at their inception, intended for finite element discretization, yet they have also been in use in computations with IGA discretization.
Local length scales and stabilization parameters targeting IGA discretization were introduced in 2017 in [15]. It goes without saying that they would also be suitable for use in computations with finite element discretization as a special case of IGA discretization. The expression introduced in [15] for the direction-dependent local length scale is from a conceptually simple three-step derivation. In the first step, the direction vector is mapped from the physical element to the parent element; in the second step, the discretization spacing along each of the parametric coordinates is accounted for; in the third step, what has been obtained in the parent element is mapped back to the physical element. Although the derivation steps were in the ST framework, reducing them to space only is straightforward. The stabilization parameters given in [13], which are largely from [15], are the current ones that have been in use in ST-VMS and ST-SUPS computations. In deriving the expressions for the local length scales, we do not need to use the standard integration para-metric space. Instead, we can use a preferred parametric space that more effectively serves the purpose, which is obtaining good expressions. This idea was introduced and shown to work well in [16,203,204]. The expressions for the local local length scales include a transformation tensor that relates the two parametric spaces. Based on this idea, expressions for the direction-dependent local length scales targeting complex-geometry B-spline meshes were introduced in [16]. We require the local length scales to be invariant with respect to element splitting in B-spline meshes. Without that invariance, the flow solution would be influenced by something that it should not be influenced by. The expressions introduced in [16] meet that requirement, and the proof was presented in [204].

Complex-geometry IGA mesh generation
This subsection, included for completeness, is mostly from [23]. While the IGA offers superior accuracy, IGA mesh generation for complex geometries is significantly more challenging than finite element mesh generation. Widely available mesh generation software packages for finite element, finite volume, and finite difference methods encourage the usage of these methods. To make IGA-based flow computations more applicable to problems with complex geometries, and consequently more practical in computational analysis of real-world problems, the IGA mesh generation will have to be less challenging and more encouraging. The CGIMG and NSVGMG were introduced to that end. We will provide a brief overview of the CGIMG here, and for the NSVGMG, we refer the interested reader to [1].
The CGIMG consists of three steps. In the first step, a block-structured mesh is generated using existing techniques for such meshes. In the second step, that mesh is projected to a NURBS mesh that is built from patches corresponding to the blocks of the block-structured mesh. In the third step, the original model surfaces are recovered, to the extent the nature of the recovered surfaces does not impede the robustness or accuracy of the flow computations. The CGIMG is normally expected to preserve the element quality and refinement distribution of the block-structured mesh. Mesh generation and mesh quality tests were included in [19,20]. The tests showed that the CGIMG is a practical IGA mesh generation method with good performance. The CGIMG has been used in computing the following classes of problems: wind turbines [55,124,128]

MDM
This subsection, included for completeness, is mostly from [125,128]. The purpose of the MDM, when it was first introduced, was to compute the long-wake flow behind an object, with the final objective being the wake computation or the ultimate objective being the computation of the wake influence on a secondary object placed far downstream. The first one of a sequence of overlapping subdomains covers the primary object, with the subsequent subdomains covering the long wake and the last one covering the secondary object, if there is one. For the first subdomain, the inflow velocity is the free-stream velocity. For each subsequent subdomain, the inflow velocity is extracted from the subdomain before it. In some cases, the outflow boundary of a subsequent subdomain might also be in the subdomain before it. When that is the case, the outflow stress is also extracted from the preceding subdomain.
A number of 3D problems were computed when and soon after the MDM was introduced: wake influence on a wing placed downstream of a larger wing [18], a cylinder wake extending 300 diameters downstream [205], and aerodynamics [206] and FSI [207] of a parachute crossing an aircraft wake. The cylinder computation, carried out at Reynolds number 140, was able to capture the second phase of the Karman vortex street observed in laboratory experiments. In the parachute computations, the parachute subdomain was translating fully inside the subdomain before it.
The thermo-fluids computations reported in [4] for a truck and its rear tires were based on a spatially multiscale MDM version with global and local domains. The global domain, containing the entire truck, was the primary domain, and the local domain, containing the rear tires, the secondary. The secondary domain was fully inside the primary domain, with some of their boundaries coinciding. The thermo-fluid computation over the primary domain was carried out with a reasonable-resolution mesh. The inflow velocity and temperature were the free-stream values, the outflow stress and normal heat flux were zero, and the top and side computational boundaries had zero normal velocity, tangential stress, and normal heat flux. The time-history data from the computation, large in volume, was stored with the ST-C. Following that, the computation over the secondary domain was carried out with a higher-resolution mesh and consequently with increased accuracy in the thermo-fluids analysis and the tire heat transfer rates. In that computation, the velocity and temperature specified at the inflow, top, and side computational boundaries were extracted, at each time step, from the stored primary-computation data at the corresponding time. This was done by evaluating the temporal NURBS representation of the primary-computation velocity and temperature at that corresponding time. At the outflow boundary, the stress was extracted from the primary-computation data, and the normal heat flux was set to zero. Where the primary-and secondarydomain boundaries coincided, the conditions specified there for the secondary domain were from what was specified for the primary domain. We note that, because in general the nodal points of the secondary-domain boundaries were not nodal points in the primary domain, the data extraction was done with the least-squares projection.
In [51], the MDM was used in aerodynamic and FSI analysis of two wind turbines placed back to back in an atmospheric boundary layer flow. The two turbines acted as the primary and secondary objects, in the way the MDM first started. The velocity data from a plane located at 10 m downstream of the primary turbine was used in [125][126][127] as the inflow velocity in the IGA-based two-domain MDM computations of the wind turbine wake. The computational framework was presented in [125], studies on spatial and temporal resolution in [126], and studies on spatial-refinement directional preference in [127]. In [68,128], the MDM was used in computation of flow over a complex terrain.

Outline of the remaining sections
The computation settings are described in Sect. 2, the results are presented in Sect. 3, and the concluding remarks are given in Sect. 4.

Car and tire models
The car and tire models are the same as in [1]; we include them here for completeness. The models are shown in Figs. 2 and 3. The tire model is provided by YOKOHAMA RUBBER CO., LTD.; the car body, wheel, and disk rotor are not from an actual CAD data. The car is 4.65 m long, 1.81 m wide, and 1.16 m high. The tire, wheel, and disk rotor dimensions are given in Table 1. Figure 4 shows the tire model. It does not have transverse grooves. The groove widths are different between the two center and two side grooves. We carry out a steady-state structural mechanics computation to obtain the tire deformation, resulting in the shape shown in Fig. 3. All four tires have the same deformation.

Problem setup
The problem setup is the same as in [1]; we include it here for completeness. The rotation speed corresponds to a linear speed of 100 km/h at the undeformed tire periphery. The measured tire contact angle is 22 • . The corresponding car speed is U 0 = 99.39 km/h. The reference frame is attached to the car. The air density and kinematic viscosity are 1.205 kg/m 3 and 1.512 × 10 −5 m 2 /s. The tire rotation period, which we will use in reporting the results, is represented by the symbol T TR . Figure 5 shows the global computational domain from [1], which we will call "car global domain (CGD)," and the car body and the tires. The domain has a length 13 times the body length, width 20 times the body width, and height 18 times the body height. The car center is located at 5 times the body length from the inflow plane. The velocity is specified at the inflow and bottom boundaries, at U 0 . The outflow boundary condition is stress-free, and the lateral and top boundary conditions are slip. The car body is a no-slip boundary. The tires, wheels, and disk rotors are also no-slip boundaries, and the meshes around them are rotating with the tire. The velocity of the rotating surfaces is obtained from the motion.
Taking the data from the CGD computation in [1] as our global data, we perform high-resolution computations for the left set of tires, in the local domain shown in Fig. 6, which we will call "car local domain (CLD)." The shape of CLD is related to the mesh used over CGD.

CGD mesh
For relational context, we describe the CGD mesh used in [1]. It was generated over five domain parts, separately, and the mesh parts were connected with SIs. Figure 7 shows the domain parts. The outer part ("O") has the simplest shape. Inside that, is the body part ("B"), which excludes the tire parts. Each of the four tire parts consists of three parts: tire rotating part ("TR"), tire stationary part ("TS"), and the wheel part ("W"). The parts are identified with abbreviations in a combination form, such as "CGD-O" and "CGD-TR." The CGD-B mesh has controlled layers of elements near the car body. There are 13 layers of elements in the normal direction, with a first-layer thickness of 0.7 mm. Figure 8 shows the CGD-B mesh. In the CGD-TR mesh, there are 10 layers in the radial direction. The first-layer thickness is about 1.5 mm. There are 6 layers in the axial direction, with a first-layer thickness of about 5.0 mm. In the computation with the ST-SI-TC-IGA, only the CGD-TR mesh was deforming. Figure  9 shows the CGD-TR and CGD-W meshes. Although the CGD-TS meshes for the front and rear tires have the same number of control points and elements and mesh structure, the domain shapes are different. The number of control points and elements for all five parts of the CGD mesh are shown in Table 2.

CLD meshes
We define CLD as a subdomain that covers the CGD-TR, CGD-TS, and CGD-W meshes and, roughly, the inner part  of the CGD-B mesh shown in Fig. 10. We decompose CLD into three overlapping subdomains "CLD1," "CLD2," and "CLD3," and perform the high-resolution computations over those subdomains, as a nested MDM sequence. Figure 11 shows CLD and its three overlapping subdomains. The overlap between CLD2 and CLD1 is 5.0 % of the CLD1 length, and between CLD3 and CLD2 4.3 % of the CLD2 length.

CLD1
The CLD1 mesh is generated over the four parts of that subdomain: "CLD1-B," "CLD1-TS," "CLD1-TR," and "CLD1-W." The subdomain parts CLD1-TS, CLD1-TR, and CLD1-W are the same as the global-domain parts CGD-TS, CGD-TR, and CGD-W, and the meshes are obtained by knot insertion from the corresponding CGD meshes. We note that CLD1-TS is, of course, the front CGD-TS. Table 3 shows the number of control points and elements for CLD1. Figure 12 shows the CLD1-B mesh. It was generated by a starting with a block-structured finite element mesh from Pointwise, converting it to the NURBS mesh with the CGIMG, followed by mesh relaxation on that with the method described in Sect. 1.7.3. In the computation with the ST-SI-TC-IGA, as how it was for the CGD-TR mesh, only the CLD1-TR part of the CLD1 mesh is deforming. Figure 13 shows the CLD1-TR and CLD1-W meshes. Compared to the CGD-TR mesh, the CLD1-TR mesh has about 3 times as many elements in the circumferential direction. There are 11 elements in the radial direction, with a first-layer thickness of about 0.68 mm, and 12 elements in the axial direction, with a first-layer thickness of about 2.5 mm. In total, the CLD1-TR mesh has about 9 times as many elements as the CGD-TR mesh has. The number of elements in the CLD1-TS and CLD1-W meshes are twice as many in all three directions. Figure 14 shows the CLD2 mesh. The number of control points and elements are 941,108 and 845,000. It was generated by a starting with a block-structured finite element mesh from Pointwise, converting it to the NURBS mesh with the CGIMG.

CLD3
The CLD3 mesh is generated over the four parts of that subdomain: "CLD3-B," "CLD3-TS," "CLD3-TR," and "CLD3-W," in the same way how the CLD1 mesh was generated. We note that the CLD3-TR and CLD3-W meshes are the same as the CLD1-TR and CLD1-W meshes. We also note that, as mentioned in Sect. 2.3, the front and rear CGD-TS meshes, despite the different domain shapes, have the same mesh structures, and therefore so do the CLD1-TS (front) and CLD3-TS (rear) meshes. In the computation with the ST-SI-TC-IGA, as how it is for the CLD1 mesh, only the CLD3-TR part of the CLD3 mesh is deforming. Table 4 shows the number of control points and elements for CLD3. Figure 15 shows the CLD3-B mesh.

CLD1-TR and CLD3-TR mesh updates
The mesh update for the CLD1-TR and CLD3-TR meshes is done, as how it was done for the CGD-TR meshes, by a combination of mesh moving with a special-purpose method and element collapse on the SI plane with the ST-SI-TC.

Computational conditions
In all computations, the method is the ST-VMS, and the stabilization parameters are those given by Eqs. (4)- (9) in [13]. The element lengths appearing in the expressions for the stabilization parameters and in some of the integrals over the SI in the ST-SI are those targeting IGA discretization (see Sect. 1.8). In the CGD computation, the time-step size was T TR /144. In the CLD computations, it is T TR /432. The number of nonlinear iterations per time step is 3, and the number of GMRES iterations per nonlinear iteration is 300. In the CGD computation, time was measured from a point reached after a sufficiently long computation at full Reynolds number. The computation was carried out until t = 5T TR . We start the CLD1 computation at t = 0, CLD2 computation at t = 0.59T TR , and the CLD3 computation at t = 1.38T TR . We compute until t = 5T TR . The initial conditions are from the CGD computation at those times. We note that 0.59T TR and 1.38T TR are the approximate durations for the flow information advected from the CLD1 inflow plane with the speed U 0 to reach the CLD2 and CLD3 inflow planes. Figure 16 shows the boundary conditions for CLD1, CLD2, and CLD3. The prescribed-velocity boundary conditions are extracted from the CGD solution, except for the CLD2 and CLD3 inflow planes. The velocities at those planes come from the CLD1 and CLD2 solutions. The stress boundary conditions are extracted from the CGD solution. We note that the CGD data we are extracting from was stored with the ST-C, using cubic B-splines in time. The time-step size for the cubic B-spline representation is T TR /48, which is three times the time-step size used in the CGD computation. Figures 17 and 18 show, for CGD and CLD, the flow patterns during the last T TR /36. The CGD computation is from [1]. We see better-resolved vortex structures in the CLD solution. Figures 19 and 20 show, for CGD and CLD, the flow patterns around the tire-road contact areas during the last T TR /36. We see more resolved vortices in the CLD solution. Comparing the larger vortex structures from the CGD and CLD computations, our observations are different for the front and rear tires. The vortex structures are similar near the front tire, but quite different near the rear tire. In the CGD computation, the vortex structures near the rear and front tires are close. This Fig. 16 Boundary conditions for CLD1, CLD2, and CLD3. 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 implies that the advected small vortex structures influence the large-scale solution near the rear tire. Figures 21 and 22 show, for CGD and CLD, the positionalaveraged shear stress over the last T TR . We see that the shear stress magnitudes are, as can be expected, quite different between the CGD and CLD solutions. Furthermore, in the CLD solution, while the large vortex structures near the front and rear tires are quite different, the distributions of the shear stress magnitude are comparable. This suggests that if we need computations with even higher resolution to reach higher accuracy in flow properties in specific regions, that can be done over smaller, third-level MDM subdomains.

Concluding remarks
We have presented high-resolution ST isogeometric analysis of car and tire aerodynamics, with many of the complexities of the actual car and tire, such as the near-actual tire geometry, road contact, and tire deformation and rotation. We focused on the tire aerodynamics, with high-resolution in both space and time. The influence of the aerodynamics of the car body was included, in the framework of the MDM, from the global computation with near-actual car body and tire geometries, carried out earlier with a reasonable mesh resolution. The high-resolution local computation was for the left compression. Except for the last three, these methods were used also in the global computation, and they played the same role in the local computation. The ST-TC, for example, as in the global computation, made the ST moving-mesh computation possible even with contact between the tire and the road, thus enabled high-resolution flow representation near the tire. Of the three last methods, the CGIMG, which served as an alternative to the NSVGMG used in the global computation, made the IGA mesh generation for the complex geometries less arduous. The MDM reduced the computational cost by focusing the high-resolution locally to where it was needed and also by breaking the local computation into its consecutive portions. The ST-C data compression made the storage of the data from the global computation less burdensome. The car and tire aerodynamics computation presented show the effectiveness of the high-resolution computational analysis framework we have built for this class of problems. That includes effectiveness in comprehensive and detailed analysis of the vortex patterns near the tires, such as seeing that properly resolved small vortex structures in the front-tire wake influence the large-scale solution near the rear tire.
Acknowledgements This work was supported in part by Rice-Waseda research agreement and International Technology Center Indo-Pacific (ITC IPAC) Contract FA520921C0010. The work was also supported in part by ARO Grant W911NF-17-1-0046 and Contract W911NF-21-C-0030 (first and fourth authors), and Top Global University Project of Waseda University (fourth author). The tire geometry in Sect. 2.1 was provided by the YOKOHAMA RUBBER CO., LTD. The authors acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing HPC resources that have contributed to the research results reported within this paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.