Space–time VMS flow analysis of a turbocharger turbine with isogeometric discretization: computations with time-dependent and steady-inflow representations of the intake/exhaust cycle

Many of the computational challenges encountered in turbocharger-turbine flow analysis have been addressed by an integrated set of space–time (ST) computational methods. The core computational method is the ST variational multiscale (ST-VMS) method. The ST framework provides higher-order accuracy in general, and the VMS feature of the ST-VMS addresses the computational challenges associated with the multiscale nature of the unsteady flow. The moving-mesh feature of the ST framework enables high-resolution computation near the rotor surface. The ST slip interface (ST-SI) method enables moving-mesh computation of the spinning rotor. The mesh covering the rotor spins with it, and the SI between the spinning mesh and the rest of the mesh accurately connects the two sides of the solution. The ST Isogeometric Analysis enables more accurate representation of the turbine geometry and increased accuracy in the flow solution. The ST/NURBS Mesh Update Method enables exact representation of the mesh rotation. A general-purpose NURBS mesh generation method makes it easier to deal with the complex geometries involved. An SI also provides mesh generation flexibility in a general context by accurately connecting the two sides of the solution computed over nonmatching meshes, and that is enabling the use of nonmatching NURBS meshes in the computations. The computational analysis needs to cover a full intake/exhaust cycle, which is much longer than the turbine rotation cycle because of high rotation speeds, and the long duration required is an additional computational challenge. As one way of addressing that challenge, we propose here to calculate the turbine efficiency for the intake/exhaust cycle by interpolation from the efficiencies associated with a set of steady-inflow computations at different flow rates. The efficiencies obtained from the computations with time-dependent and steady-inflow representations of the intake/exhaust cycle compare well. This demonstrates that predicting the turbine performance from a set of steady-inflow computations is a good way of addressing the challenge associated with the multiple time scales.


Introduction
The challenges faced in computational flow analysis of turbocharger turbines include unsteady flow through a complex geometry, the need for high-resolution flow representation near the rotor surface, high Reynolds numbers, and multiscale flow behavior. The flow unsteadiness comes from the intake/exhaust cycle and the flow in the turbine. An additional challenge is associated with the multiple time scales involved. The computational analysis needs to cover a full intake/exhaust cycle, which is much longer than the turbine rotation cycle because of high rotation speeds. Therefore computations with time-dependent representation of the intake/exhaust cycle require long-durations in the turbine time scale.
The parts of these challenges faced in turbine computations in a more general context have been addressed also by other researchers, with approaches ranging from using a single blade with spatially-periodic boundary conditions (see, e.g., [1][2][3][4][5][6][7][8]) to "sliding interfaces" (see, e.g., [9][10][11][12]). The precursors of the work presented in this article were reported in [13][14][15][16]. The challenge associated with the multiple time scales was addressed in [16] by going through with the long-duration computation. Here we propose to address that challenge in a different way, by calculating the turbine efficiency from computations with steady-inflow representation of the intake/exhaust cycle.

ST-VMS and ST-SUPS
The deforming-spatial-domain/stabilized ST (DSD/SST) method [26][27][28] was introduced for computation of flows with moving boundaries and interfaces (MBI), including fluid-structure interactions (FSI). In MBI computations the DSD/SST functions as a moving-mesh method. Moving the fluid mechanics mesh to follow an interface enables mesh-resolution control near the interface and, consequently, high-resolution boundary-layer representation near fluidsolid interfaces. The stabilization components of the original DSD/SST are the Streamline-Upwind/Petrov-Galerkin (SUPG) [29] and Pressure-Stabilizing/Petrov-Galerkin (PSPG) [26] stabilizations, which are used widely. Because of the SUPG and PSPG components, the original DSD/SST is now called "ST-SUPS." The ST-VMS is the VMS version of the DSD/SST. The VMS components of the ST-VMS are from the residual-based VMS (RBVMS) method [30][31][32][33]. The ST-VMS has two more stabilization terms beyond those in the ST-SUPS, and the additional terms give the method better turbulence modeling features. The ST-SUPS and ST-VMS, because of the higher-order accuracy of the ST framework (see [17,18]), are desirable also in computations without MBI.
In the flow analysis presented here, the ST framework provides higher-order accuracy in a general context. The VMS feature of the ST-VMS addresses the computational challenges associated with the multiscale nature of the unsteady flow. The moving-mesh feature of the ST framework enables high-resolution computation near the rotor surface.

ST-SI
The ST-SI was introduced in [20], in the context of incompressible-flow equations, to retain the desirable moving-mesh features of the ST-VMS and ST-SUPS when we have spinning solid surfaces, such as a turbine rotor. The mesh covering the spinning surface spins with it, retaining the high-resolution representation of the boundary layers. The starting point in the development of the ST-SI was the version of the ALE-VMS for computations with sliding interfaces [45,46]. Interface terms similar to those in the ALE-VMS version are added to the ST-VMS to account for the compatibility conditions for the velocity and stress at the SI. That accurately connects the two sides of the solution. An ST-SI version where the SI is between fluid and solid domains was also presented in [20]. The SI in this case is a "fluid-solid SI" rather than a standard "fluidfluid SI" and enables weak enforcement of the Dirichlet boundary conditions for the fluid. The ST-SI introduced in [21] for the coupled incompressible-flow and thermaltransport equations retains the high-resolution representation of the thermo-fluid boundary layers near spinning solid surfaces. These ST-SI methods have been applied to aerodynamic analysis of vertical-axis wind turbines [11,12,20], thermo-fluid analysis of disk brakes [21], flow-driven string dynamics in turbomachinery [117][118][119], flow analysis of turbocharger turbines [13][14][15][16], flow around tires with road contact and deformation [113,[120][121][122], lubrication fluid dynamics [123], aerodynamic analysis of ram-air parachutes [124], and flow analysis of heart valves [110,114,115].
In the ST-SI version presented in [20] the SI is between a thin porous structure and the fluid on its two sides. This enables dealing with the porosity in a fashion consistent with how the standard fluid-fluid SIs are dealt with and how the Dirichlet conditions are enforced weakly with fluid-solid SIs. This version also enables handling thin structures that have T-junctions. This method has been applied to incompressibleflow aerodynamic analysis of ram-air parachutes with fabric porosity [124]. The compressible-flow ST-SI methods were introduced in [125], including the version where the SI is between a thin porous structure and the fluid on its two sides. Compressible-flow porosity models were also introduced in [125]. These, together with the compressible-flow ST SUPG method [127], extended the ST computational analysis range to compressible-flow aerodynamics of parachutes with fabric and geometric porosities. That enabled ST computational flow analysis of the Orion spacecraft drogue parachute in the compressible-flow regime [125,126].
In the computations here, the ST-SI enables spinning the mesh covering the solid surface together with the surface. This retains the high-resolution representation of the flow near the surface, and the SI between the spinning mesh and the rest of the mesh accurately connects the two sides of the solution.

ST-IGA and STNMUM
The ST-IGA is the integration of the ST framework with isogeometric discretization. It was introduced in [17]. Computations with the ST-VMS and ST-IGA were first reported in [17] in a 2D context, with IGA basis functions in space for flow past an airfoil, and in both space and time for the advection equation. Using higher-order basis functions in time enables getting full benefit out of using higher-order basis functions in space. This was demonstrated with the stability and accuracy analysis given in [17] for the advection equation.
The ST-IGA with IGA basis functions in time enables a more accurate representation of the motion of the solid surfaces and a mesh motion consistent with that. This was pointed out in [17,18] and demonstrated in [22][23][24]. It also enables more efficient temporal representation of the motion and deformation of the volume meshes, and more efficient remeshing. These motivated the development of the STN-MUM [22][23][24]; the name was given in [25]. The STNMUM has a wide scope that includes spinning solid surfaces. With the spinning motion represented by quadratic NURBS in time, and with sufficient number of temporal patches for a full rotation, the circular paths are represented exactly. A "secondary mapping" [17,18,22,37] enables also specifying a constant angular velocity for invariant speeds along the circular paths. The ST framework and NURBS in time also enable, with the "ST-C" method, extracting a continuous representation from the computed data and, in large-scale computations, efficient data compression [19,21,113,[117][118][119]128]. The STNMUM and the ST-IGA with IGA basis functions in time have been used in many 3D computations. The classes of problems solved are flapping-wing aerodynamics for an actual locust [22,23,37,101], bioinspired MAVs [24,99,100,102] and wing-clapping [103,104], separation aerodynamics of spacecraft [93], aerodynamics of horizontal-axis [10,25,99,100] and vertical-axis [11,12,20] wind-turbines, thermo-fluid analysis of ground vehicles and their tires [19,113], thermo-fluid analysis of disk brakes [21], flow-driven string dynamics in turbomachinery [117][118][119], and flow analysis of turbocharger turbines [13][14][15][16].
The ST-IGA with IGA basis functions in space enables more accurate representation of the geometry and increased accuracy in the flow solution. It accomplishes that with fewer control points, and consequently with larger effective element sizes. That in turn enables using larger time-step sizes while keeping the Courant number at a desirable level for good accuracy. It has been used in ST computational flow analysis of turbocharger turbines [13][14][15][16], flow-driven string dynamics in turbomachinery [118,119], ram-air parachutes [124], spacecraft parachutes [126], aortas [110,111], heart valves [110,114,115], tires with road contact and deformation [121,122], and lubrication fluid dynamics [123]. The image-based arterial geometries used in patient-specific arterial FSI computations do not come from the zero-stress state (ZSS) of the artery. A number of methods [35,37,[129][130][131][132][133][134][135][136][137][138] have been proposed for estimating the ZSS required in the computations. Using IGA basis functions in space is now a key part of some of the newest ZSS estimation methods [136][137][138][139] and related shell analysis [140].
In the flow analysis presented here, the ST-IGA enables more accurate representation of the turbine geometry, increased accuracy in the flow solution, and using larger timestep sizes. The STNMUM enables exact representation of the mesh rotation.

General-purpose NURBS mesh generation method
To make the ST-IGA use, and in a wider context the IGA use, even more practical in computational flow analysis with complex geometries, NURBS volume mesh generation needs to be easier and more automated. To that end, a generalpurpose NURBS mesh generation method was introduced in [14]. The method is based on multi-block-structured mesh generation with existing techniques, projection of that mesh to a NURBS mesh made of patches that correspond to the blocks, and recovery of the original model surfaces. The recovery of the original surfaces is to the extent they are suitable for accurate and robust fluid mechanics computations. The method is expected to retain the refinement distribution and element quality of the multi-block-structured mesh that we start with. Because there are ample good techniques and software for generating multi-block-structured meshes, the method makes general-purpose mesh generation relatively easy. Mesh-quality performance studies for 2D and 3D meshes, including those for complex models, were presented in [15]. A test computation for a turbocharger turbine and exhaust manifold was also presented in [15], with a more detailed computation in [16]. The mesh generation method was used also in the pump-flow analysis part of the flowdriven string dynamics presented in [119]. The performance studies and test computations demonstrated that the generalpurpose NURBS mesh generation method makes the IGA use in fluid mechanics computations even more practical. The general-purpose NURBS mesh generation method is used also in the turbocharger-turbine flow analysis presented here.

ST-SI-IGA
An SI provides mesh generation flexibility in a general context by accurately connecting the two sides of the solution computed over nonmatching meshes. This type of mesh generation flexibility is especially valuable in complex-geometry flow computations with isogeometric discretization, removing the matching requirement between the NURBS patches without loss of accuracy. This feature was used in the flow analysis of heart valves [110,114,115], turbocharger turbines [13][14][15][16], spacecraft parachute compressible-flow analysis [126], and flow around tires with road contact and deformation [122]. It is used also in the turbocharger-turbine flow analysis presented here.

Steady-inflow representation of the intake/exhaust cycle
In calculating the turbine efficiency from the computations with the steady-inflow representation of the intake/exhaust cycle, we first perform a set of steady-inflow computations at different flow rates. We then obtain the turbine efficiencies associated with those computations and tabulate them as a function of the tip-speed ratio. Next we calculate the tipspeed ratio for the moving-averaged inflow velocity coming from the time-dependent representation of the intake/exhaust cycle. The turbine efficiency corresponding to that timedependent tip-speed ratio is obtained by interpolation from the tabulated values of the steady-inflow turbine efficiencies.

Computations presented
For comparison purposes, we include from [16] some of the results reported for the computation with the time-dependent representation of the intake/exhaust cycle. The steady-inflow computations are carried out at 12 different flow rates.

Outline of the remaining sections
The governing equations are given in Sect. 2. For completeness, in Sect. 3 we describe the ST-VMS. The computations are presented in Sect. 4, and the concluding remarks in Sect. 5. In the Appendix, we describe the ST-SI and give the stabilization parameters used in the ST-VMS and ST-SI.

Governing equations
Let Ω t ⊂ R n sd be the spatial domain with boundary Γ t at time t ∈ (0, T ), where n sd is the number of space dimensions. The subscript t indicates the time-dependence of the domain. The Navier-Stokes equations of incompressible flows are written on Ω t and ∀t ∈ (0, T ) as where ρ, u and f are the density, velocity and body force. The stress tensor σ σ σ (u, p) = − pI + 2με ε ε(u), where p is the pressure, I is the identity tensor, μ = ρν is the viscosity, ν is the kinematic viscosity, and the strain rate ε ε ε(u) = ∇ ∇ ∇u + (∇ ∇ ∇u) T /2. The essential and natural boundary conditions for Eq. (1) are represented as u = g on (Γ t ) g and n · σ σ σ = h on (Γ t ) h , where n is the unit normal vector and g and h are given functions. A divergence-free velocity field u 0 (x) is specified as the initial condition.

ST-VMS
For completeness, we include, mostly from [17], the ST-VMS: where are the residuals of the momentum equation and incompressibility constraint. The test functions associated with the velocity and pressure are w and q. A superscript "h" indicates that the function is coming from a finite-dimensional space. The symbol Q n represents the ST slice between time levels n and n + 1, (P n ) h is the part of the lateral boundary of that slice associated with the stress boundary condition h, and Ω n is the spatial domain at time level n. The superscript "e" is the ST element counter, and n el is the number of ST elements. The functions are discontinuous in time at each time level, and the superscripts "−" and "+" indicate the values of the functions just below and just above the time level.

Remark 1
The ST-SUPS can be obtained from the ST-VMS by dropping the eighth and ninth integrations.
We also use in the computations the ST-SI, which we include in Appendix A. The stabilization parameters, embedded in the ST-VMS and ST-SI, are included in Appendix B.

Problem setup
The computation with the time-dependent representation of the full intake/exhaust cycle was reported in [16]. Here we carry out the computations with the steady-inflow representation of the intake/exhaust cycle, at 12 different values of the tip-speed ratio: r ω/U IN , where r ω and U IN are the tip speed and inflow velocity. We will compare the efficiencies extracted from the computations with these two representations. Figure 1 shows the turbocharger-turbine model used in the computations with the steady-inflow representation. The model consists of the volute, stator, turbine and the exhaust pipe. The computations with the time-dependent representation had an exhaust gas purifier, instead of an exhaust pipe, and it also had an exhaust manifold. When we compare the efficiencies from the computations with the two representations, we use the same flow region.
The rotor diameter is 30 mm and the rotor speed is 30,000 rpm, which translates to a turbine rotation period of T TR = 2.0×10 −3 s. The gas density and kinematic viscosity are 0.9 kg/m 3  We note that for the operational conditions we have here, incompressible-flow modeling is reasonable. Besides, comparing the laboratory test data from scaled-up models to data from actual-scale models, it was observed that the Mach number made little difference in the turbine efficiency

Efficiency definition
In reporting the turbine efficiency, we first define two quantities, one instantaneous and the other interval-based: where Γ INF and Γ OUTF are the inflow and outflow boundaries used in the efficiency calculation. Similarly, where P(t) is the instantaneous power extracted from the turbine. Then, a general form of the interval-based efficiency is written as We note that we selected the locations of Γ INF and Γ OUTF with the objective of having the smallest possible domain in the efficiency calculation so that we minimize the effect of the manifold and gas purifier.   Figure 2 shows the quadratic NURBS control mesh for the turbine model. The mesh was generated with the generalpurpose NURBS mesh generation method [14,15] described in Sect. 1.4. Table 1 shows the number of control points and elements in different parts of the mesh. Figure 3 shows the four SIs of the mesh. Two of the SIs have an actual slip. The other two are just for mesh generation purpose and connect nonmatching meshes. They are used for independent meshing in the volute region of the computational domain. The STNMUM enables exact representation of the mesh rotation.

Computational conditions
In temporal representation of the mesh rotation, we use quadratic NURBS basis functions. There are 90 time steps per rotation, which is equivalent to a time-step size of 2.22×10 −5 s. The number of nonlinear iterations per time step is 4, with 500, 500, 600 and 800 GMRES iterations

Computations with the steady-inflow representation
The boundaries Γ INF and Γ OUTF are shown in Fig. 4. Figure 5 shows, at the 12 tip-speed ratios, the turbine efficiency averaged over one rotation, based on Eq. (9), where t 3 = t 1 , t 4 = t 2 , and t 2 − t 1 = T TR .

Computation with the time-dependent representation
For completeness, we include from [16] some of the results reported for the computation with the time-dependent representation of the intake/exhaust cycle. The symbol T will denote the duration of the intake/exhaust cycle, and T = 6.0×10 −2 s. Figure 6 shows the velocity magnitude for the full cycle. Figure 7 shows the turbine efficiency η IB (0, T , t, t).
The turbine efficiency corresponding to the time-dependent tip-speed ratio is obtained by interpolation from the steadyinflow turbine efficiencies displayed in Fig. 5. We note from Figs. 8, 9 and 10 that if the averaging period is too short or too long, the efficiencies from the steady-inflow and time-dependent representations do not match well. However, if we choose the right averaging period, which is in this case the turbine rotation period, the efficiencies match very well. We also note that the effect of the manifold

Efficiency (%)
Efficiency from time-dependent representation Efficiency from steady-inflow representation and gas purifier is not significant, provided of course that we report the efficiency based on the same flow region, shown in Fig. 4. Overall, the turbine efficiency can be predicted from a set of computations with the steady-inflow representation of the intake/exhaust cycle, with little difference compared to the efficiency from the time-dependent representation.

Concluding remarks
One of the main computational challenges encountered in turbocharger-turbine flow analysis is the multiple time scales involved. The computational analysis needs to cover a full intake/exhaust cycle, which is much longer than the turbine rotation cycle because of high rotation speeds. Therefore computations with time-dependent representation of the intake/exhaust cycle require long-durations in the turbine time scale. In this article we focused on that challenge. Many of the other computational challenges were addressed earlier by an integrated set of ST computational methods, and we used those methods in the computations here. The core computational method is the ST-VMS. The ST framework provides higher-order accuracy in general, the VMS feature of the ST-VMS addresses the computational challenges associated with the multiscale nature of the unsteady flow, and the moving-mesh feature of the ST framework enables high-resolution computation near the rotor surface. The ST-SI enables moving-mesh computation of the spinning rotor. The mesh covering the rotor spins with it, and the SI between the spinning mesh and the rest of the mesh accurately connects the two sides of the solution. The ST-IGA enables more accurate representation of the turbine geometry and increased accuracy in the flow solution. The STNMUM enables exact representation of the mesh rotation. A general-purpose NURBS mesh generation method makes it easier to deal with the complex geometries involved. An SI also provides mesh generation flexibility in a general context by accurately connecting the two sides of the solution computed over nonmatching meshes, and that is enabling the use of nonmatching NURBS meshes in the computations.
The challenge associated with the multiple time scales was addressed earlier by going through with the longduration computation. As another way of addressing that challenge, we proposed here to calculate the turbine efficiency from computations with steady-inflow representation of the intake/exhaust cycle. For that, we first performed a set of steady-inflow computations at different flow rates. We then obtained the turbine efficiencies associated with those computations and tabulated them as a function of the tip-speed ratio. Next we calculated the tip-speed ratio for the moving-averaged inflow velocity coming from the timedependent representation of the intake/exhaust cycle. We obtained the turbine efficiency corresponding to that timedependent tip-speed ratio by interpolation from the tabulated values of the steady-inflow turbine efficiencies. The efficiencies obtained from the computations with time-dependent and steady-inflow representations of the intake/exhaust cycle compare well. This demonstrates that predicting the turbine performance from a set of steady-inflow computations is a good way of addressing the challenge associated with the multiple time scales.
Acknowledgements This work was supported in part by Grant-in-Aid for Challenging Exploratory Research 16K13779 from Japan Society for the Promotion of Science; Grant-in-Aid for Scientific Research (S) 26220002 from the Ministry of Education, Culture, Sports, Science and Technology of Japan (MEXT); Council for Science, Technology and Innovation (CSTI), Cross-Ministerial Strategic Innovation Promotion Program (SIP), "Innovative Combustion Technology" (Funding agency: JST); and Rice-Waseda research agreement. This work was also supported in part by Grant-in-Aid for Early-Career Scientists 19K20287 (first author) and ARO Grant W911NF-17-1-0046 and Top Global University Project of Waseda University (third author). 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. The HPC resources provided by Nagoya University High Performance Computing Research Project for Joint Computational Science also contributed to obtaining the results reported in the paper.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix A: ST-SI
In describing the ST-SI (see [20]), labels "Side A" and "Side B" will represent the two sides of the SI. The ST-SI version of the formulation given by Eq. (3) includes added boundary terms corresponding to the SI. The boundary terms for the two sides are first added separately, using the test functions w h A and q h A and w h B and q h B . Then, putting together the terms added to each side, the complete set of terms added becomes where and h is the element length at the SI (see Appendix B.4).
Here, (P n ) SI is the SI in the ST domain, v is the mesh velocity, γ = 1, and C is a nondimensional constant. For explanation of the added SI terms, see [20,120].
On solid surfaces where we prefer weak enforcement of the Dirichlet conditions [20,42] for the fluid, we use the ST-SI version where the SI is between the fluid and solid domains. Using the terminology of Sect. 1.2, that is a fluid-solid SI. This version is obtained (see [20]) by starting with the terms added to Side B and replacing the Side A velocity with the velocity g h coming from the solid domain. Then the SI terms added to Eq. (3) to represent the weakly-enforced Dirichlet conditions become

Appendix B: Stabilization parameters
We include from [141] the element metric tensor in space and in the ST framework. These are used in Appendices B.3 and B.4 in calculation of the stabilization parameters and element lengths.

Appendix B.1: Element metric tensor in space
Components of the Jacobian matrix Q are written as where ξ j is the parametric coordinate in jth direction. We first scale it with a matrix D to take into account the polynomial order or other factors such as the dimensions of the element domain in the parametric space: Remark 2 In the case of finite elements with Lagrange polynomials of order p i in ith direction, with a parametric space of − 1 ≤ ξ i ≤ 1, the scaling matrix can be chosen as WithQ and the direction r to which we want to associate the length, we define the element length (see [141]) as where Remark 3 From this derivation, what we get with D = I has been used in many methods for calculating the stabilization parameters (see, for example, [37]). In those methods, a scaling factor taking the polynomial order into account is applied to the element length, and here we do the scaling in the parametric space, for each of the parametric directions.
Sweeping over all the directions represented by r, we obtain the minimum and maximum element lengths: They are equivalent to and where λ max and λ min are the maximum and minimum eigenvalues of the argument matrix.

Remark 4
In the implementation, we take measures to keep the calculated element length between h MIN and h MAX .

Appendix B.2: Element metric tensor in the ST framework
The ST Jacobian matrix is where θ is the parametric coordinate in time, and the mesh The ST scaling matrix is given as and the scaling becomeŝ The ST metric tensor is defined as

Appendix B.4: ST-SI
The element length used in the ST-SI is given as These were introduced in [