Gas turbine computational flow and structure analysis with isogeometric discretization and a complex-geometry mesh generation method

A recently introduced NURBS mesh generation method for complex-geometry Isogeometric Analysis (IGA) is applied to building a high-quality mesh for a gas turbine. The compressible flow in the turbine is computed using the IGA and a stabilized method with improved discontinuity-capturing, weakly-enforced no-slip boundary-condition, and sliding-interface operators. The IGA results are compared with the results from the stabilized finite element simulation to reveal superior performance of the NURBS-based approach. Free-vibration analysis of the turbine rotor using the structural mechanics NURBS mesh is also carried out and shows that the NURBS mesh generation method can be used also in structural mechanics analysis. With the flow field from the NURBS-based turbine flow simulation, the Courant number is computed based on the NURBS mesh local length scale in the flow direction to show some of the other positive features of the mesh generation framework. The work presented further advances the IGA as a fully-integrated and robust design-to-analysis framework, and the IGA-based complex-geometry flow computation with moving boundaries and interfaces represents the first of its kind for compressible flows.


Introduction
The designs of commercial and military vehicles are continuously improved to deliver increased performance. The vehi-B Kenji Takizawa kenji.takizawa@tafsm.org Yuri Bazilevs yuri_bazilevs@brown.edu Tayfun E. Tezduyar tezduyar@tafsm.org 1 cles are often costly to build and maintain, and, as a result, optimal performance and operational reliability become the key aspects of the vehicle design in the commercial and military mobility applications. As the vehicle and engine designs get more sophisticated, the corresponding computational analysis methods must also mature to support advances in the mobility technology. Recent advances in geometry modeling, mesh generation, computational mechanics of fluids and structures, multiphysics modeling, and high-performance computing create a unique opportunity for the development of the next-generation predictive simulation methods and tools for general mobility applications involving complex geometry and fluid-structure interaction (FSI). These newgeneration methods and tools, described in some detail and demonstrated on a single-stage gas turbine design in this article, enable high-fidelity computational analysis with greatly improved representation of the complex geometries and multiphysics phenomena.

Moving-mesh methods
Using the terminologies and categorizations used in [1][2][3], a method for flows with moving boundaries and interfaces (MBI), including FSI, can be an interface-tracking (movingmesh) method or an interface-capturing (nonmoving-mesh) method, or possibly a combination of the two. In a movingmesh method, as the interface moves and the fluid mechanics domain changes its shape, the mesh moves to adjust to the shape change and to follow (i.e. "track") the interface. Moving the fluid mechanics mesh to follow a fluid-solid interface enables us to control the mesh resolution near the interface, have high-resolution representation of the boundary layers, and obtain accurate solutions in such critical flow regions. The arbitrary Lagrangian-Eulerian (ALE) method, with the ALE finite element method introduced in 1981 [4], is one of the earliest and most commonly used moving-mesh methods.

Compressible-flow SUPG method
The compressible-flow SUPG method was first introduced, in the context of conservation variables, in a 1982 NASA technical report [130], with a concise version published as a 1983 AIAA paper [131] and more thorough version with additional examples as a 1984 journal paper [132]. That 1982 method is now called "(SUPG) 82 ." Several SUPGlike compressible-flow methods were developed after that. Taylor-Galerkin method [133], for example, is very similar, and under certain conditions is identical, to one of the versions of (SUPG) 82 . Another example of the subsequent SUPG-like compressible-flow methods in conservation variables is the streamline-diffusion method [134].
When the flow field is expected to have a discontinuity, SUPG and VMS methods are often supplemented with a discontinuity-capturing (DC) term, which is also called "shock-capturing term" in the context of compressible flows. Supplementing the SUPG with a DC term started as early as 1986 [135,136]. The DC term played a significant role in the evolution of the compressible-flow SUPG. At its inception, (SUPG) 82 was not used with a DC term, and the need for additional stabilization at the shocks was obvious from the test computations. Then (SUPG) 82 was recast in entropy variables and supplemented with a DC term [137]. This, as expected, resulted in better shock profiles. In a 1991 ASME paper [138], (SUPG) 82 was supplemented with a very similar DC term. All the test computations presented in [138,139], which were among the most widely used test computations at that time, showed that, with the DC term, (SUPG) 82 was very comparable in accuracy to (SUPG) 82 recast in entropy variables. The stabilized and DC methods introduced in [136] for the advection-diffusion-reaction equation included safeguards to prevent "compounding", i.e., augmentation of the SUPG effect by the DC effect when the advection and discontinuity directions are not orthogonal.
The compressible-flow ST-SUPG method [140] is essentially the same as the compressible-flow DSD/SST method, but without necessarily implying a mesh motion. The compressible-flow DSD/SST is a straightforward mixture of the DSD/SST concept and the compressible-flow SUPG. The first 3D MBI computation with the compressible-flow ST-SUPG was reported in 1996 [141] for two high-speed trains passing each other in a tunnel.
In a recent article [142], the ALE-SUPG formulation of compressible flows using pressure-primitive variables with DC was developed and applied to the simulation of a single-stage gas turbine. Weakly-enforced essential boundary condition and sliding-interface formulations for the compressible-flow equations were added in [142] to enable gas turbine simulations with boundary-layer flow and statorrotor interaction. In [143], the formulation was applied to study off-design performance of the gas turbine stage, and in [144] it was used to simulate the full-rotorcraft aerodynamics. Here we build on the developments in [142] and introduce improvements to the DC, weakly-enforced essential boundary condition, and sliding-interface operators to increase the accuracy and robustness of the ALE-based compressible flow formulation, especially when it is used in combination with the Isogeometric Analysis (IGA) [21,44,45,145,146].

ALE Sliding-Interface method
The ALE sliding-interface (ALE-SI) method [21,22] was introduced, in the context of incompressible-flow equations, to retain the desirable moving-mesh features of the ALE-VMS and ALE-SUPS in computations involving spinning solid surfaces, such as a turbine rotor. The mesh covering the spinning surface spins with it, retaining the high-resolution representation of the boundary layers, while the mesh on the other side of the SI remains unaffected. This is accomplished by adding to the ALE-VMS formulation interface terms that account for the compatibility conditions for the velocity and stress at the SI, accurately connecting the two sides of the solution. The ALE-SI has been applied to wind turbine aerodynamics and FSI [28][29][30][31] (more specifically, VAWTs [33,34,37,38], floating wind turbines [39], wind turbines in atmospheric boundary layers [32][33][34], and fatigue damage in wind turbine blades [43]), hydrodynamics and FSI of a hydraulic arresting gear [61,62], hydrodynamics of tidal-stream turbines with free-surface flow [63], bioinspired FSI for marine propulsion [65,66], and compressible flows with emphasis on gas-turbine analysis [142,143] and fullrotorcraft simulation [144].
The ST Slip Interface (ST-SI) method (the acronym "SI" will also refer to "slip interface") was introduced in [88], also in the context of incompressible-flow equations. The interface terms added to the ST-VMS formulation are similar to those in the ALE-SI. An ST-SI version where the SI is between fluid and solid domains was also presented in [88]. The SI in that case is a "fluid-solid SI" rather than a standard "fluid-fluid SI" and enables weak enforcement of the Dirichlet boundary conditions for the fluid. The ST-SI introduced in [113] for the coupled incompressible-flow and thermal-transport equations retains the high-resolution representation of the thermo-fluid boundary layers near spinning solid surfaces. These ST-SI methods have been applied to aerodynamic analysis of VAWTs [33][34][35][36]88,89], thermo-fluid analysis of disk brakes [113], flow-driven string dynamics in turbomachinery [35,36,[114][115][116], flow analysis of turbocharger turbines [117][118][119][120][121], flow around tires with road contact and deformation [106,[122][123][124][125], fluid films [125,126], aerodynamic analysis of ram-air parachutes [41,42,127], and flow analysis of heart valves [57,58,102,[107][108][109] and ventricle-valve-aorta sequences [104]. In the ST-SI version presented in [88] 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 incompressible-flow aerodynamic analysis of ram-air parachutes with fabric porosity [41,42,127]. The compressible-flow ST-SI methods were introduced in [128], 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 [128]. These, together with the compressible-flow ST-SUPG, 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 [128,129].

Moving-mesh Isogeometric Analysis
The IGA [21,44,45,145,146] was proposed in an effort to "bridge the gap" between computer-aided design (CAD) and the finite element method. As such, at its foundation, the IGA, unlike other computational methods, has a strong link to CAD, which is routinely used for the geometric design of many objects, including vehicles and their mechanical components. For the discretization of the solution fields, the IGA uses B-splines, non-uniform rational B-splines (NURBS), and related basis-function technology, which are essentially piece-wise polynomial or rational functions joined with higher-order continuity. For many problems arising in both fluid and solid/structural mechanics, the higher-order approximation accuracy and continuity of the splines result in high per-degree-of-freedom accuracy necessary to resolve the spatially multiscale phenomena present in many physical systems. This yields superior robustness, necessary for engineering applications. Computational flow analysis 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.
The success with IGA basis functions in space motivated the integration of the ST methods with isogeometric discretization, which is broadly called "ST-IGA." The ST-IGA was introduced in [10]. Computations with the ST-VMS and ST-IGA were first reported in [10] 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 higherorder basis functions in time enables deriving full benefit from using higher-order basis functions in space. This was demonstrated with the stability and accuracy analysis given in [10] 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 [10,11] and demonstrated in [90,91,93]. It also enables more efficient temporal representation of the motion and deformation of the volume meshes, and more efficient remeshing. These motivated the development of the ST/NURBS mesh update method (STNMUM) [90,91,93], with the name coined in [85]. 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" [2,10,11,90] 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 [106,[112][113][114][115][116]147]. 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 [2,[90][91][92], bioinspired MAVs [86,87,93,94] and wing-clapping [95,96], separation aerodynamics of spacecraft [80], aerodynamics of horizontal-axis [30,[85][86][87] and vertical-axis [33][34][35][36]88,89] wind turbines, thermo-fluid analysis of ground vehicles and their tires [41,106,112], thermo-fluid analysis of disk brakes [113], flow-driven string dynamics in turbomachinery [35,36,[114][115][116], flow analysis of turbocharger turbines [117][118][119][120][121], and flow analysis of coronary arteries in motion [110].

Complex-geometry NURBS mesh generation method
Mesh generation in the IGA, such as NURBS meshes generation, is not as established and straightforward as mesh generation in the standard discretization methods such as the finite differences and finite elements. To make the IGAbased flow and FSI computations even more powerful and more practical, the mesh generation has to be less arduous, more straightforward, more automated, and basically as established as it is now for the finite difference and finite element methods. A complex-geometry NURBS mesh generation method was introduced in [118] to make the ST-IGA use, and in a wider context the IGA use, more practical in computational flow analysis with complex geometries. The method is based on 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 block-structured mesh that we start with. Because there are ample good techniques and software for generating block-structured meshes, the method makes complex-geometry mesh generation relatively easy. Mesh-generation tests and mesh-quality studies reported in [118,119] clearly demonstrated the good performance of the method. Practicality of the method in flow computations in real-world applications was also demonstrated in [118,119], with computational analysis of a turbocharger turbine with exhaust manifold, a disk-gap-band parachute, and a human aorta with branches. The method was used in more detailed computation of the aorta in [57,58,102,103] and turbocharger turbine in [120], with additional computational analysis in [121]. The method was shown to be very effective also in IGA-based computational flow analysis of other classes of problems, such as the pump-flow analysis part of the flow-driven string dynamics presented in [35,36,116], flow analysis of a tsunami-shelter VAWT [36,89], and ventricle-valve-aorta flow analysis [104].

ALE-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 context of the ST-IGA and under the name "ST-SI-IGA", in flow analysis of heart valves [57,58,102,[107][108][109], ventriclevalve-aorta flow analysis [104], turbocharger-turbine flow analysis [117][118][119][120][121], spacecraft parachute compressible-flow analysis [129], flow analysis around tires with road contact and deformation [124,125], and flow analysis of a tsunamishelter VAWT [36,89]. It is also used, in the context of the ALE-IGA and under the name "ALE-SI-IGA", in the gas turbine flow analysis presented here.

Stabilization and DC parameters and element lengths
In all the semi-discrete and ST stabilized and VMS methods discussed in the earlier subsections, an embedded stabilization parameter, almost universally known as "τ ", plays a significant role (see [2]). This parameter involves a measure of the local length scale (also known as "element length") and other parameters such as the element Reynolds and Courant numbers. The interface terms in the ALE-SI and ST-SI also involve element length, in the direction normal to the interface. Various element lengths and τ s were proposed, starting with those in [5,161] and [130][131][132], followed by the ones introduced in [135,136]. In many cases, the element length was seen as an advection length scale, in the flow-velocity direction. The set of τ s introduced in [130][131][132] in conjunction with (SUPG) 82 is now called "τ 82 ". The τ definition introduced in [136], which is for the advective limit and is now called "τ SUGN1 " (and the corresponding element length is now called "h UGN "), automatically yields lower values for higher-order finite element basis functions (see [162,163]).
The τ used in [138] with (SUPG) 82 was a slightly modified version of τ 82 . Subsequent minor modifications of τ 82 took into account the interaction between the DC and (SUPG) 82 terms in a fashion similar to how it was done in [136] for the advection-diffusion-reaction equation. Until 2004, all these slightly modified versions of τ 82 were always used with the same DC parameter, which was introduced in the 1991 ASME paper [138] and is now called "δ 91 ." This DC parameter was derived from the one given in [137] for the entropy variables.
Calculating the τ s based on the element-level matrices and vectors was introduced in [164] in the context of the advection-diffusion equation and the Navier-Stokes equations of incompressible flows. These definitions are expressed in terms of the ratios of the norms of the matrices or vectors. They automatically take into account the local length scales, advection field and the element Reynolds number. The definitions based on the element-level vectors were shown [164,165] to address the difficulties reported at small time-step sizes. A second element length scale, in the solution-gradient direction and called "h RGN ," was introduced in 2001 [8,166]. Recognizing this as a diffusion length scale, a new stabilization parameter for the diffusive limit, "τ SUGN3 ," was introduced in [8,167], to be used together with τ SUGN1 and "τ SUGN2 ," the parameters for the advective and transient limits. For the stabilized ST methods, "τ SUGN12 ," representing both the advective and transient limits, was also introduced in [8].
New ways of calculating the τ and DC parameter to be used with (SUPG) 82 were introduced in 2004 [ [167][168][169]. The new τ s, now categorized under the label "τ 04 ," have a matrix structure for viscous flows and reduce to a scalar for inviscid flows. The new DC parameters were of two types: one defined in a style the Discontinuity-Capturing Directional Dissipation [8,169,170] parameter was defined, and one that is now called "YZβ" DC parameter. The YZβ DC parameter is residual-based, and it is simpler than δ 91 . It has options for smoother or sharper computed shocks. A number of 2D and 3D test computations with YZβ DC were reported in [171][172][173]. These computations showed that in addition to being simpler than δ 91 , the YZβ DC parameter was superior in accuracy. The computations reported in [171][172][173] were based on the compressible-flow ST SUPG.
Some new options for the stabilization parameters used with the SUPS and VMS were proposed in [9,84,85,90,112]. These include a fourth τ component, "τ SUGN4 " [112], which was introduced for the VMS, considering one of the two extra stabilization terms the VMS has compared to the SUPS. They also include stabilization parameters [112] for the thermaltransport part of the VMS for the coupled incompressibleflow and thermal-transport equations.

Direction-dependent element lengths for isogeometric discretization
The stabilization and DC parameters and element lengths discussed in the previous subsection were all originally intended for finite element discretization but quite often used also for isogeometric discretization. The element lengths and stabilization and DC parameters introduced in [187] target isogeometric discretization but are also applicable to finite element discretization. They were introduced in the context of the advection-diffusion equation and the Navier-Stokes equations of incompressible flows. The direction-dependent element length expression was outcome of a conceptually simple derivation. The key components of the derivation are mapping the direction vector from the physical ST element to the parent ST element, accounting for the discretization spacing along each of the parametric coordinates, and mapping what has been obtained in the parent element back to the physical element. The test computations presented in [187] for pure-advection cases, including those with discontinuous solution, showed that the element lengths and stabilization parameters proposed result in good solution profiles. The test computations also showed that the "UGN" parameters give reasonably good solutions even with NURBS basis functions. The stabilization parameters given in [124], which were mostly from [187], were the latest ones designed in conjunction with the ST-VMS. The direction-dependent local-length-scale expressions introduced in [188] for B-spline meshes are outcome of a clear and convincing derivation and more suitable for element-level evaluation. The new expressions are based on a preferred parametric space, instead of the standard integration parametric space, and a transformation tensor that represents the relationship between the integration and preferred parametric spaces. We do not want the element splitting to influence the actual discretization, which is represented by the control or nodal points. Therefore, the local length scale should be invariant with respect to element split-ting. That invariance is a crucial requirement in element definition, because unlike the element definition choices based on implementation convenience or computational efficiency, it influences the solution. It was proven in [189] that the local-length-scale expressions introduced in [188] meet that requirement.

Outline of the remaining sections
The NURBS mesh generation method is described in Sect. 2. The governing equations of compressible flows, including their quasi-linear and ALE forms, are summarized in Sect. 3. The compressible-flow computational method, including the SUPG, DC, weak-boundary-condition and SI operators, is presented in Sect. 4. The gas turbine geometry is described in Sect. 5. The meshes for the flow and structure analyses are presented in Sect. 6, together with the mesh quality evaluations. The results from the flow simulations are presented in Sect. 7. Structural mechanics simulation of the rotor free vibration is presented in Sect. 8. In Sect. 9, as an example of calculating the direction-dependent local length scales in isogeometric discretization, we calculate the Courant number in the turbine computation. The concluding remarks are given in Sect. 10.

NURBS mesh generation method
For completeness, we include, mostly from [118,119], a shortened description of the NURBS mesh generation method. First we provide a brief introduction to the basic concepts and terminologies of NURBS meshes. Figure 1 shows, as an example, a NURBS control mesh for a stator blade of a turbocharger turbine, depicting the NURBS patches and control points. Figure 2 shows the corresponding physical mesh, depicting the physical patches and elements. The relationship between the control and physical meshes is given by NURBS basis functions, defined on each patch by the knot vectors of the patch, with each space direction having its own knot vector.

Basic method
We start with a block-structured mesh (grid). Such meshes are very common in finite difference and finite volume computations. The mesh quality can be measured in terms of the grid point distribution, grid line orthogonality, element aspect ratio, and the adjacent-element-length ratio. Generating high-  Fig. 1. The lines represent the element boundaries quality meshes requires some skills and experience, but there are ample good techniques and software for generating blockstructured meshes. We assume that the block-structured mesh we start with is the outcome of such techniques or software and is of high quality. The mesh consists of trilinear elements. We see each block as a precursor to a NURBS patch. The time it would take to generate the block-structured mesh would be the time the user would spend on generating the finite difference, finite volume or finite element mesh. Beyond that, the NURBS mesh generation process is fully automated and takes essentially no time.
The second step involves a sequence of projections. For those projections, we define a common parametric space between each block and the corresponding patch. We choose the parametric space to be ζ ζ ζ ∈ [0, 1] n pd , where n pd is the number of parametric dimensions. Then the knot vector of the patch for a given direction would be of the form where p is the polynomial order for the NURBS basis functions used in that direction. Two methods were proposed in [118,119] to determine the knot values. We refer the reader interested in knowing those methods to [118,119].

Projections
For each patch, the projections are done hierarchically. First the control points at the corners of the patch are set to the same locations as the corner grid points of the block. Then the edges between the corners are projected: Here χ χ χ h ∈ Ω ⊂ R n sd , where n sd is the number of space dimensions, is the position representation in the block, x h ∈ Ω is the position representation in the patch, and w h is the test function of the projection. The parametric coordinate ζ used in the integration is a general representation for all edges.
Next the surfaces between the edges are projected: Here, the combination of the parametric coordinates ζ 1 and ζ 2 is a general representation for all surfaces. The projection sequence is completed with the projection of the volumes between the surfaces:

Merging the patches
For two adjacent blocks, the projections described in Sect. 2.1.1 result in the same control point positions over the surface they share, provided that the two blocks have the same knot vectors over the shared surface. This might happen automatically or require some additional steps, depending on the method used in determining the knot values of the parametric space. The control point variables are declared to be the same between the adjacent patches over the shared surface, which results in C 0 continuity for the basis functions across the surface.

Remark 1
Alternatively, as proposed and used in [108,117], we can declare the control point variables over the shared surface to be separate, resulting in C −1 continuity, and connect them with an SI.

Element reduction
We know from our experience that, to generate good quality NURBS meshes with the method proposed, in many cases we may need more elements in the finite element mesh than the target number of NURBS elements. Therefore we apply a reduction factor r to the number of elements in the block, which we can possibly be different for each direction in To make the number of elements an integer, after applying the factor r , we raise the value to the next integer. After the element reduction, we use the modified versions of the methods for determining the knot values of the parametric space, and we refer the reader interested in knowing how the modifications are done to [118,119].

Methods for recovering the exact surfaces
In the NURBS mesh generation method we have, we can recover the exact surfaces, instead of just relying on the surfaces represented by the block-structured mesh.

Special method for arc surfaces
Suppose we have a finite element mesh with a boundary that is an arc in the exact representation (see Fig. 3). We first represent the arc by three control points and the corresponding quadratic NURBS basis functions. The arc template is shown in Fig. 4, and the corresponding knot vector and weights are and For details, see [2,11]. Figure 5 shows the NURBS representation of the arc, superimposed on the finite element mesh.
Then we insert to the knot vector given by Eq. (5) the knots of the initial NURBS mesh. From that, we obtain the new control points and weights representing the boundary. The new mesh can be seen in Fig. 6.

General method for CAD surfaces
If the CAD surface basis function space is a subset of our NURBS mesh basis function space, we can recover the exact surface by knot insertion to the CAD space to obtain our NURBS space. If not, the best we can do is to project the CAD surface to our NURBS space.

Governing equations of compressible flows
The Navier-Stokes equations of compressible flows are often expressed using the vector of conservation variablesŨ: where ρ is the density, u i is the i th velocity component, i = 1, . . . , n sd , e tot = e + u 2 /2, and e is the internal-energy density.
We also introduce the vector of pressure-primitive variables Y: where p and T are the pressure and temperature. Pressure, density, and temperature are related through the equation of state. Here we use the ideal-gas equation of state: where R is the ideal-gas constant. We assume a calorically perfect gas and define the internal-energy density as where c v = R/(γ −1) is the specific heat at constant volume, and γ is the heat capacity ratio. Throughout the article, we use (·) ,t to denote a partial time derivative holding the spatial position x fixed, and we use (·) ,i to denote the spatial gradient.
With the above definitions, the Navier-Stokes equations of compressible flows, which represent the conservation of mass, momentum, and energy, can be written in terms ofŨ as whereF adv i andF diff i are the vectors of advective and diffusive fluxes,S is the source term, and Res is the strong-form residual.
The fluxes are defined as where δ i j is the Kronecker delta. The viscous stress τ i j and heat flux q i are given as where μ is the dynamic viscosity, λ = −2μ/3 based on Stokes' hypothesis, and κ is the thermal conductivity.
Introducing the mass and momentum conservation into the energy equation of Eq. (11), we get where is the contribution of the stress power in the energy equation, and Res is the modified strong-form residual.

Quasi-linear and ALE forms of the governing equations
The quasi-linear form of Eq. (16) is Explicit expressions for the matrices appearing in the quasilinear forms above can be found in [142].
With the space-time Piola transformation and following the steps in [45], the convective ALE form of Eq. (16) is written as whereû i is the i th component of the domain velocityû and (·) ,t x denotes a partial time derivative holding the referential positionx fixed. The quasi-linear form of Eq. (23) in terms of U is and in terms of Y, The corresponding strong-form residual is For notational convenience, we define the matricesÂ ALE

Compressible-flow computational method
Let Ω ⊂ R n sd denote the domain in the current configuration, and let Γ be its boundary. In addition, let S and V denote the vector-valued pressure-primitive trial-and test-function spaces, respectively. The weak form of Eq. (25) is stated as follows: find Y ∈ S, such that ∀W ∈ V, where Here Γ h is the subset of Γ where the viscous stress and heat flux boundary conditions H are enforced weakly.

Remark 2
In compressible-flow computations, at the outflow boundary one often sets the pressure, viscous stress, and heat flux. One way to accomplish this is to set the nodal or control-point pressure to the desired values, and enforce the remaining conditions weakly by defining H as ( 3 0 ) An alternative approach, which we use here, is to relax the strong imposition of the outlet pressure and define H as Here n i is the ith component of outward unit surface normal vector n, and the subscript "out" denotes the prescribed outflow quantities.
Equations (27)- (29) are used as a point of departure to develop a formulation suitable for the discretization using the IGA or finite elements. In what follows, we present the additional method components, which are the SUPG stabilization, DC, weak enforcement of the essential boundary conditions, and the SI method.

SUPG operator
We assume the time-dependent domain Ω is divided into N el spatial finite elements each denoted by Ω e , and define the SUPG operator as This SUPG operator is an extension of (SUPG) 82 (see Sect. 1.2) from its original version with the scalar stabilization parameter τ 82 (see Sect. 1.7) to a version with a matrix stabilization parameter τ τ τ SUPG . When working with Y, the stabilization parameter is whereτ τ τ SUPG is the stabilization parameter when working with U, which was given in [190] aŝ Here Δt is the time-step size, C I is a positive constant derived from an appropriate element-wise inverse estimate [191], and G i j are the components of the element metric tensor. The expression in Eq. (34) requires calculation of the squareroot-inverse of a 5×5 matrix in 3D, which is done iteratively using a modified version of the Denman-Beavers algorithm [192,193]. More details can be found in [142].

DC operator
Following the approach for the SUPG operator, we first define the DC operator for working with U: whereK DC is the matrix DC parameter. Changing the variables we are working with from U to Y gives which, in turn, defines the matrix DC parameter for working with Y: We assume a diagonal form forK DC : wherê Here h is the element length, C C , C M , and C E are O(1) positive constants, Res w and ∇ ∇ ∇U w are the weighted norms Res w = c 2 |Res 1 | + u Res 2:n sd +1 + Res n sd +2 , κ cap is the maximum value of the DC parameter: where G i j = G i j −1 , and c is the acoustic speed.

Remark 3
Note that the weighted norm appropriately scales the components that have different dimensions.

Remark 4
The DC parameter upper bound or "cap",κ cap is a multi-dimensional generalization of the upwinding viscosity h(u + c)/2. While it is expected that, on average, the residual-based definition of the DC parameter will stay well below the upwinding limit, the division by the gradient norm can lead to local spikes, which are mitigated by the cap. The cap idea was introduced and successfully employed in [194], a recent reference focusing on developing residualbased shock-capturing methods for solids. We also note that the introduction of the cap reduces the degree of nonlinearity associated with the DC terms and, as a result, improves convergence of the Newton-Raphson iterations.

Remark 5
The DC parameter originates from the "CAU DC" [195], which, in turn, is just an extension of δ 91 (see Sect. 1.7), originally designed in the context of computing steady-state solutions, from its form based on the steady-state residual to its form based on the time-dependent residual. The DC parameter may also be viewed as YZβ DC with β = 1.

Weak-boundary-condition operator
Weakly-enforced essential boundary conditions act as nearwall models for under-resolved boundary-layer flows while converging to their strongly-enforced counterparts at optimal rate with mesh refinement. Here we develop an improved, relative to [142], weak-boundary-condition formulation for compressible flows. The essential boundary conditions for the velocity and temperature are enforced on Γ D ∈ Γ . Let W = [q w T w θ ] T be the vector of test functions with w = [w 1 w 2 w 3 ] T . The weak-boundary-condition operator is given as The first line on the right-hand side of Eq. (45) represents a convective contribution to the weak-boundary-condition operator. We note that the weak-boundary-condition operator is defined in terms of U and G = [ρ ρg T ρc v T b ] T , where g and T b are the prescribed velocity and temperature. Furthermore, {Â n } − is the "negative" part ofÂ n =Â i n i . It Γ n is computed, based on the eigenvalue decomposition Γ n T −1 associated with hyperbolic equation systems (see, e.g., [196]), using the expression where, defining u n = u −û · n, is the diagonal matrix of the eigenvalues ofÂ n , is the matrix constructed using the corresponding eigenvectors as columns, is the matrix of transformation from the vector [ρ u T p] T to U, and {Γ n } − is the negative part of Γ n , which is easily computed. The remaining terms in Eq. (45) correspond to the weak enforcement of essential boundary conditions coming form the viscous stress and heat-flux contributions as given in [142,144].

SI operator
The SI method used in this work follows the approach for the weak-boundary-condition operator. We consider two subdomains that are in relative motion and share an SI, denoted by Γ I . We use the subscripts 1 and 2 to represent the two sides of the SI. Compatibility of U, stresses, and heat fluxes at the SI is weakly-enforced with the following SI operator: where t(u, n) = μ ∇ ∇ ∇u + (∇ ∇ ∇u) T + λ(∇ ∇ ∇ · u)I n is the viscous stress operator. The convective contributions to the SI operator, given by the first two terms on the right-hand side of Eq. (50), make use of {Â n } − evaluated on each side of the SI. The remainder of the terms account for the viscous stress and heat flux contributions.

Semi-discrete formulation and time integration
The final semi-discrete form of the Navier-Stokes equations of compressible flows may be stated as: where S h and V h are the discrete counterparts of S and V. The generalized-α method [197,198] is used for time integration. At each time step, the nonlinear equation system is solved with the Newton-Raphson method. At each Newton-Raphson iteration the resulting linear equation system is solved iteratively using the GMRES search technique [199] with nodal-block-diagonal preconditioning [200].

Turbine geometry
The turbine geometry comes from the interactive geometry modeling platform described in [201] and was used in earlier finite element simulations reported in [142]. The geometry and dimensions are shown in Fig. 7. The axial length is 211.36 mm, the casing radius is 95.524 mm, and the shaft inner and outer radii are 38.862 and 77.724 mm. The stator has 24 blades, and the rotor 34. This is a smaller gas turbine design, similar to that used as part of a turboshaft for Black Hawk and Apache helicopters. We note that the gap between Fig. 7 Geometry, dimensions and problem setup the stator and rotor blades is larger than that in the design used earlier in [142].

Fluid mechanics mesh
We divide the domain into four parts and generate a blockstructured mesh for each part. Figure 8 shows the four parts.
We have five SIs in the mesh. One of them is with an actual slip, and the other four are just for mesh generation purposes, giving us the flexibility to have different mesh resolutions on the two sides of the SI. Figure 9 shows the five SIs. Due to the rotational periodicity, the meshes around the stator and rotor blades are created by generating only one stator blade mesh and one rotor blade mesh and repeating those in a rotationally periodic fashion. In the cylindrical parts, the NURBS weights are selected to represent the exact geometry. Figure 10 shows the control mesh. The number of control points is 716,706, and the number of elements is 393,274. Figures 11 and 12 show the mesh around the stator and rotor blades.

Rotor structural mechanics mesh
We divide the domain into three parts and generate a blockstructured mesh for each part. Figure 13 shows the three parts.
There is no SI. Due to the rotational periodicity, the mesh around the blades is created by generating only one blade  Figure 14 shows the control mesh. The number of control points is 246,738, and the number of elements is 126,922. Figure 15 shows the blade mesh.

Mesh quality evaluation
We evaluate the mesh quality based on the mesh line orthogonality and the aspect ratio. The mesh lines are defined along the parametric coordinates ξ α in each element: where α = 1, . . . , n pd . For a given pair of two directions α and β, we define the angle φ αβ between the two mesh lines as Since the choice of the parametric directions is arbitrary, in the mesh quality visualization, we report the maximum value of the mesh line angle. We also report the L p norm, which is defined in a standard fashion: To define the aspect ratio, we first define the length in each parametric direction: With that, for a given pair of two directions α and β, the aspect ratio r αβ is defined as In the mesh quality visualization, we again report the maximum value. The L p norm we report in this case is defined in a logarithmic fashion:

Remark 6
We note that here we use the "transformation tensor" definition D = I in calculating the local length scale (see [121] for the terminology and symbol definitions). There is no technical obstacle to using other definitions, such as "RQD-MAX" in [121]. Figures 16-19 show the mesh line angle and aspect ratio in the parts of the mesh around the stator and rotor blades. Tables 1 and 2 show the L 2 and max norms of the mesh line angle and aspect ratio.  Tables 3 and 4 show the L 2 and max norms of the mesh line angle and aspect ratio.

Problem setup
The problem setup is given in Fig. 7 We conduct two simulations, one using the quadratic NURBS mesh described in Sect. 6, and the other using a linear tetrahedral mesh similar to the one used in [142]. The meshes near the stator and rotor blades are shown in Fig. 22, and the mesh data is given in Table 5. The time-step size is 3×10 −7 s in both simulations.

Figures 23 and 24
show the global distribution of the flow speed, pressure, temperature, and Mach number for the computations with the NURBS and tetrahedral meshes. The results are in good agreement intrinsically, as well as with earlier simulations of this case. We note that the solution appears continuous across all five SIs of the NURBS mesh. Figure 25 shows the temperature distribution in the stator-rotor region. The quadratic NURBS mesh with a carefully constructed boundary layer mesh is able to much better represent the thin thermal boundary layer near the stator blade than the tetrahedral mesh. The more diffuse thermal boundary layer obtained with the tetrahedral mesh leads to much less flow acceleration behind the stator blade and slower flow in the rotor channels. This can be seen in Fig. 26, which shows, for both simulations, the relative flow speed and streamlines in the stator and rotor channels. The slower flow in the rotor channels in the tetrahedral-mesh simulation translates to lower suction pressure and underestimation of the shaft torque relative to the NURBS-mesh simulation. The shaft power obtained is 354.0 kW with the NURBS mesh and 242.3 kW with the tetrahedral mesh.
Gas turbine performance can also be assessed from the adiabatic efficiency of the turbine stage η ad , calculated using the expression (see [202]) Here, T RBO 0 and p RBO 0 are the total temperature and pressure at the rotor blade outlet, and T SBI 0 and p SBI 0 at the stator blade inlet. The total temperature and pressure are calculated as and where M is the local Mach number. The computed values of η ad are 98.6 and 98.9% in the simulations with the NURBS and tetrahedral meshes. This is higher than the adiabatic efficiency of 85% reported for a similar case in [143], however, as we pointed out earlier, the simulation here has slip boundary conditions on the shaft and casing and a larger gap between

Free-vibration analysis
We carry out a free-vibration analysis of the gas turbine rotor with the structural mechanics NURBS mesh described in Sect. 6. We use the material properties of Nimonic C-263,  Figure 27 shows the displacement magnitude corresponding to one of the lower and one of the higher natural frequencies. The lower frequency corresponds to a global deformation mode, while the higher frequency has deformation localized to the rotor blades. Color range from blue to red indicates the normalized-displacement range from low to high

An example of local-length-scale calculation in isogeometric discretization: Courant number
In the context of the fluid mechanics mesh and computation presented in the earlier sections, we give an example of calculating the direction-dependent local length scales in isogeometric discretization. We calculate the Courant number based on the local length scale in the direction of the flow velocity computed. The local length scale is based on the RQD-MAX version of the method in [121]. The time-step size, from Sect. 7, is 3×10 −7 s. Figures 28 and 29 show the Courant number.

Concluding remarks
We have made further advances in establishing the IGA as a fully-integrated and robust design-to-analysis framework. We have integrated a recently introduced NURBS mesh generation method for complex-geometry IGA and a state-of-the-art moving-mesh compressible-flow computation method. With that, we have successfully carried out IGA simulation of a gas turbine stage similar to the one used as part of a turboshaft for Black Hawk and Apache helicopters. In the course of this development, we have further improved the DC, weak-boundary-condition, and SI methods for compressible flows. A comparison between the IGA and finite element simulations of the gas turbine stage shows that the NURBS mesh generation method, which is as easy to use as a block-structured mesh generation method in finite volume or finite element analysis, enables higher-order and higher-continuity representation and superior accuracy. The NURBS discretization is able to represent better the thin thermal boundary layers and higher flow velocity on the suction side of the stator blade, leading to higher flow rate in the rotor channels and resulting in higher shaft power compared to the tetrahedral finite element discretization. A free-vibration analysis of the gas turbine rotor carried out using the structural mechanics NURBS crated with the mesh generation method shows that the method can be used very effectively also in structural mechanics analysis, producing informative global and local deformation patterns. The 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/.