Compressible-flow geometric-porosity modeling and spacecraft parachute computation with isogeometric discretization

One of the challenges in computational fluid–structure interaction (FSI) analysis of spacecraft parachutes is the “geometric porosity,” a design feature created by the hundreds of gaps and slits that the flow goes through. Because FSI analysis with resolved geometric porosity would be exceedingly time-consuming, accurate geometric-porosity modeling becomes essential. The geometric-porosity model introduced earlier in conjunction with the space–time FSI method enabled successful computational analysis and design studies of the Orion spacecraft parachutes in the incompressible-flow regime. Recently, porosity models and ST computational methods were introduced, in the context of finite element discretization, for compressible-flow aerodynamics of parachutes with geometric porosity. The key new component of the ST computational framework was the compressible-flow ST slip interface method, introduced in conjunction with the compressible-flow ST SUPG method. Here, we integrate these porosity models and ST computational methods with isogeometric discretization. We use quadratic NURBS basis functions in the computations reported. This gives us a parachute shape that is smoother than what we get from a typical finite element discretization. In the flow analysis, the combination of the ST framework, NURBS basis functions, and the SUPG stabilization assures superior computational accuracy. The computations we present for a drogue parachute show the effectiveness of the porosity models, ST computational methods, and the integration with isogeometric discretization.

The DSD/SST method was introduced in [11,16,17], intended for computation of flows with moving boundaries and interfaces (MBI), including FSI. Because it functions as a moving-mesh method in MBI computations, the fluid mechanics mesh follows the fluid-solid interfaces, enabling mesh-resolution control and accurate flow representation near those interfaces. What made the original DSD/SST method "stabilized" were the Streamline-Upwind/Petrov-Galerkin (SUPG) [18] and Pressure-Stabilizing/Petrov-Galerkin (PSPG) [11] components, and for that the method is now also called "ST-SUPS." The ST Variational Multiscale (ST-VMS) method [14,15] is the VMS version of the DSD/SST method, with the VMS components coming from the residual-based VMS (RBVMS) method [19][20][21][22].
The core and special ST computational methods have been playing by far the most prevalent role in FSI analysis of spacecraft parachutes (as documented in [1][2][3][4][5][6][7][8][9][10] and references therein). Actually, parachute FSI analysis in 3D with the ST computational methods started about 10 years earlier than the FSI analysis of spacecraft parachutes. In that earlier 10-year period, the ST computational methods played the most prevalent role in FSI analysis of various types of personnel and cargo parachutes (as documented in the earlier references cited in [1][2][3][4][5][6][7][8][9][10]). However, computational FSI analysis of spacecraft parachutes involves challenges not only beyond those in a typical FSI analysis, but also beyond those in parachute FSI analysis, and even beyond those in large-parachute FSI analysis. One of those challenges is the "geometric porosity," a design feature created by the hundreds of gaps and slits that the flow goes through. Because FSI analysis with resolved geometric porosity would require resolving the flow that goes through the hundreds of gaps and slits as they change their shapes during the computation, it would be exceedingly time-consuming. That makes accurate geometric-porosity modeling essential. One of the special methods targeting spacecraft parachutes, the Homogenized Modeling of Geometric Porosity (HMGP) [3,98], was introduced to address this computational challenge. The HMGP was a key contributor to successful computational analysis and design studies of the Orion spacecraft parachutes since 2007 (see [1,2] and references therein, and [3][4][5][6][7][8][9]).
Until recently, ST computational analysis of spacecraft parachutes focused on the Orion spacecraft main parachutes, which are the parachutes used for landing, in the incompressible-flow regime, which is where the main parachutes operate. At the higher-altitudes, drogue parachutes will be used, and that will mostly be in the compressible-flow regime. These parachutes have a ribbon construction and 24 gores, with 52 ribbons in each gore, where a gore is the slice of the parachute canopy between two radial reinforcement cables running from the parachute vent to the skirt. This construction results in hundreds of gaps that the flow goes through, creating a geometric-porosity challenge similar to the one faced in FSI analysis of the main parachutes. Furthermore, there are three wider gaps along the gore, created by removing ribbons. Drogue FSI computations with the ST computational methods were first presented in [6,7,9], for the incompressible-flow part of the flight envelope.
Geometric-porosity models and ST computational methods for compressible-flow aerodynamics of spacecraft parachutes were introduced recently [10] in the context of finite element discretization. The key new component of the ST computational framework was the compressible-flow ST Slip Interface (ST-SI) method, introduced in conjunction with the compressible-flow ST SUPG method.
The compressible-flow ST SUPG method is essentially the same as the compressible-flow DSD/SST method, but without necessarily implying a mesh motion. The compressibleflow DSD/SST method is a straightforward mixture of the DSD/SST concept and the compressible-flow SUPG method. The first 3D computation with the compressibleflow DSD/SST method was reported in 1996 [99] for two high-speed trains passing each other in a tunnel. The interested reader can find in [10] a summary of when and in what context the compressible-flow SUPG method was introduced [100][101][102], how it evolved with the addition of a shockcapturing term [103,104] and with new stabilization and shock-capturing parameters [105][106][107], and what test computations [108][109][110] were reported.
The ST-SI method [72] was introduced to retain the desirable moving-mesh features of the ST-VMS method (and its reduced version, ST-SUPS) when we have spinning solid surfaces, such as a turbine rotor or a tire. With the ST-SI method, the mesh covering the spinning solid surface spins with it and we retain the high-resolution representation of the boundary layers. The starting point in the development of the ST-SI method was the version of the ALE-VMS method designed for computations with sliding interfaces [32,33]. In the ST-SI method, interface terms similar to those in the ALE-VMS version are added to the ST-VMS formulation to account for the compatibility conditions for the velocity and stress. The SI between the spinning mesh and the rest of the mesh accurately connects the two sides. An ST-SI version where the SI is between fluid and solid domains with weakly-imposed Dirichlet boundary conditions for the fluid was also presented in [72]. The ST-SI method introduced in [90] for the coupled incompressible-flow and thermal-transport equations addresses the challenge involved in high-resolution representation of the thermo-fluid boundary layers near spinning solid surfaces. These ST-SI methods have been successfully applied to aerodynamic analysis of vertical-axis wind turbines [72], thermo-fluid analysis of disk brakes [90], flowdriven string dynamics in turbomachinery [91], flow analysis of turbocharger turbines [92][93][94], flow around tires with road contact and deformation [95,96], aerodynamic analysis of ram-air parachutes [97], and heart valve flow analysis [84,86,87].
The ST-SI methods have additional desirable features. The SI provides mesh generation flexibility in a general context by accurately connecting nonmatching meshes. This feature was used in the flow analysis of a heart valve [84,87] and a turbocharger turbine [92][93][94]. 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 (see [93,94]). In another version of the ST-SI method presented in [72], the SI is between a thin porous structure and the fluid on its two sides. With this, the fabric porosity is dealt with in a fashion consistent with how the standard two-sided SIs are dealt with and how the Dirichlet conditions are enforced weakly. Furthermore, this version of the ST-SI method enables handling thin structures that have T-junctions. This method has been successfully used in incompressible-flow aerodynamic analysis of ram-air parachutes with fabric porosity [97].
The compressible-flow ST-SI methods were introduced in [10], 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 [10]. These, together with the compressible-flow ST SUPG method, extended the ST computational analysis range to compressible-flow aerodynamics of parachutes with fabric and geometric porosities. That enabled successful ST computational flow analysis of the Orion spacecraft drogue parachute in the compressible-flow regime [10]. The computations were in the context of finite element discretization.
The ST Isogeometric Analysis (ST-IGA), which is the integration of the ST methods with isogeometric discretization, was introduced in [14], inspired by the success of using IGA basis functions in space [23,32,48,111]. First computations with the ST-VMS method and ST-IGA were reported in [14] 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. The stability and accuracy analysis given [14] for the advection equation showed that using higher-order basis functions in time would be essential in getting full benefit out of using higher-order basis functions in space.
In the early stages of the ST-IGA, the emphasis was on IGA basis functions in time. As pointed out in [14,15] and demonstrated in [73,74,76], higher-order NURBS basis functions in time provide a more accurate representation of the motion of the solid surfaces and a mesh motion consistent with that. They also provide more efficiency in temporal representation of the motion and deformation of the volume meshes, and better efficiency in remeshing. That is how the ST/NURBS Mesh Update Method (STNMUM) was introduced and demonstrated in [73,74,76]. The name "STNMUM" was given in [69]. The STNMUM has a wide scope that includes spinning solid surfaces. With the spinning motion represented by quadratic NURBS basis functions in time, and with sufficient number of temporal patches for a full rotation, the circular paths are represented exactly, and a "secondary mapping" [2,14,15,73] enables also specifying a constant angular velocity for invariant speeds along the 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 [89][90][91]112]. The STNMUM and desirable features of the ST-IGA with IGA basis functions in time have been demonstrated in many 3D computations. The classes of problems solved are flapping-wing aerodynamics for an actual locust [2,[73][74][75], bioinspired MAVs [70,71,76,77] and wing-clapping [78,79], separation aerodynamics of spacecraft [4], aerodynamics of horizontal-axis [41,[69][70][71] and vertical-axis [72] wind-turbines, thermo-fluid analysis of ground vehicles and their tires [89], thermo-fluid analysis of disk brakes [90], flow-driven string dynamics in turbomachinery [91], and flow analysis of turbocharger turbines [92][93][94].
The ST-VMS method and ST-IGA with IGA basis functions in space have been successfully utilized in ST computational flow analysis of turbocharger turbines [92][93][94], ram-air parachutes [97], heart valves [84,86,87], and tires with road contact and deformation [96]. The turbocharger turbine analysis was based on the integration of the ST-SI method and ST-IGA. The IGA basis functions were used in the spatial discretization of the fluid mechanics equations and also in the temporal representation of the rotor and spinning-mesh motion. That enabled accurate representation of the turbine surfaces and rotor motion and increased accuracy in the flow solution. The meshes used in the turbine analysis presented in [93,94] was created by the generalpurpose NURBS mesh generation method introduced in [93] for complex-geometry flow computations. The ST-SI method is a key player in discretization with this generalpurpose mesh generation method; it removes, without loss of accuracy, the matching requirement between the NURBS patches. The ram-air parachute analysis was based on the integration of the ST-IGA, the ST-SI version that weakly enforces the Dirichlet conditions, and the ST-SI version that accounts for the porosity of a thin structure. The ST-IGA with IGA basis functions enabled, with relatively few number of unknowns, accurate representation of the parafoil geometry and increased accuracy in the flow solution. The volume mesh needed to be generated both inside and outside the parafoil, and the mesh generation inside was challenging near the trailing edge because of the narrowing space. Using IGA basis functions addressed that computational challenge and still kept the element density near the trailing edge at a reasonable level. The heart valve analysis was based on the integration of the ST-SI method, ST Topology Change (ST-TC) method [78,79,85,95], and the ST-IGA. The "ST-SI-TC-IGA," beyond enabling a more accurate representation of the surfaces and increased accuracy in the flow solution, kept the element density in the narrow spaces near the contact areas at a reasonable level. When solid surfaces come into contact, the elements between the surface and the SI collapse. Before the elements collapse, the boundaries could be curved and rather complex, and the narrow spaces might have high-aspect-ratio elements. With NURBS elements, it was possible to deal with such adverse conditions rather effectively. The computational flow analysis of tires with road contact and deformation presented in [96] was also based on the ST-SI-TC-IGA, with essentially the same desirable features as in the heart valve computational flow analysis.
In this article, we integrate the compressible-flow porosity models and ST computational methods with isogeometric discretization. We apply that to flow analysis of a drogue parachute. We use quadratic NURBS basis functions in the computations.
The governing equations are given in Sect. 2, and the porosity models in Appendix A. The compressible-flow ST SUPG and ST-SI methods are given in Appendices B and C. The computations are presented in Sect. 3 and the concluding remarks are given in Sect. 4.

Governing equations
Let Ω t ⊂ R n sd be the spatial domain with boundary Γ t at time t ∈ (0, T ). The subscript t indicates the timedependence of the domain. The symbols ρ, u and p will represent the density, velocity and pressure, respectively, and ε ε ε(u) = (∇ ∇ ∇u) + (∇ ∇ ∇u) T /2 is the strain-rate tensor. The stress tensor is defined as σ σ σ (u, p) = −pI + T, where I is the identity tensor, and T is the Newtonian viscous tensor: T = λ(∇ ∇ ∇ · u)I + 2με ε ε(u). Here λ and μ (= ρν) are the viscosity coefficients, ν is the kinematic viscosity, and it is assumed that λ = −2μ/3.
The Navier-Stokes equations of compressible flows can be written on Ω t and ∀t ∈ (0, T ) as where U = (ρ, ρu 1 , ρu 2 , ρu 3 , ρe) is the vector of conservation variables, e is the total energy per unit volume, and F i and E i are, respectively, the Euler and viscous flux vectors: Here δ i j are the components of I, q i are the components of the heat flux vector, and T i j are the components of T.
The equation of state typically corresponds to the ideal gas assumption. The term R represents all other components that might enter the equations, including the external forces.
Equation (1) can further be written in the form where The essential and natural boundary conditions for Eq. H are complementary subsets of the boundary Γ t , n i are the components of the unit normal vector n, and G and H are given functions. A function U 0 (x) is specified as the initial condition.
The porosity models we use in conjunction with the governing equations are given in Appendix A. The compressible-flow ST SUPG and ST-SI methods are given in Appendices B and C.

Parachute description
The Orion spacecraft drogue parachute is a variable porosity conical ribbon parachute [6,7,9] with a nominal diameter of 23 ft. It has 24 gores, each composed of 52 2-inch horizontal ribbons that are spaced and retained at close intervals by seven parallel, equidistant vertical tapes. The ribbon ends are stitched to the radial lines, which provide the primary longitudinal stiffness. The parachute also includes a vent band, which connects the 24 radial lines terminated at the vent. Figure 1 shows the parachute configuration. The spacing between the ribbons is varied in four levels. The 13 ribbons that are closest to the vent are spaced 0.3 inches apart, and the next three groups of 13 ribbons are spaced 0.4, 0.5 and 0.6 inches apart. Additionally, there are three locations along the gore where ribbons are removed, increasing the overall geometric porosity. These "missing" ribbons allow a localized increase in flow, which helps prevent the boundary layer from reattaching to the canopy, increasing parachute stability.
The ribbons are modeled with membrane elements. The upper, middle and lower ribbons have slightly different material properties. The various lines, tapes and bands are all modeled with cable elements, with material properties that vary by functionality. The material properties were obtained from NASA.

Flight conditions
The drogue parachute is designed to be used at a wide range of altitudes and Mach numbers. In [10] we had three altitudes, 10,000, 20,000 and 35,000 ft, and three free-stream Mach numbers, 0.3, 0.5 and 0.7. In this article, the free-stream Mach number is 0.3, and we have the same three altitudes as in [10]. We use the same notation as in [10], AM11, AM21, and AM31, where the first digit denotes the altitude, and the second digit the Mach number. Table 1 shows the flight conditions for the three cases.

Structural mechanics computations
Structural mechanics computations for the drogue parachute are conducted to obtain a deformed shape prior to the fluid mechanics computations.

Problem setup
The structural mechanics mesh, shown in Fig. 2, is a cubic NURBS mesh. It has 91,612 control points and is composed of 6576 membrane elements, 13,249 cable elements, and 1 payload element. The mesh resolution for the ribbons in the radial direction is 1 element.
As the fluid dynamics load on the canopy, we apply a uniform pressure difference equal to the dynamic pressure corresponding to the free-stream density and velocity, Δ p = 1 2 ρ ∞ u ∞

Results
Parachute configurations obtained from the structural mechanics computation in the three cases are shown in Figs. 3 and 4. The elongation of the three parachutes look quite similar. On the other hand, we see slight differences in the canopy shape. We see more roundness in AM31 than in AM11 because of the lower dynamic pressure.

One-gore fluid mechanics computations
For the purpose of calculating the porosity coefficients in modeling the geometric porosity in each of the three cases, we compute the flow field for a one-gore (15 • ) "slice" of the canopy configuration obtained from the structural mechanics computation.

Problem setup
The fluid mechanics domain for the one-gore slice is a slice of a cylinder with radius of 17.25 ft, height 57.5 ft, and an axial circular hole that has a very small diameter. The distance between the parachute skirt and the inflow boundary is approximately 15 ft. At the inflow, we specify ρ and u in the Euler fluxes, and we drop the terms that we identify as the 3rd, 4th and 5th terms in Eq. (50). At the outflow, we specify stress condition with atmospheric pressure and zero normal heat flux. On the side boundaries, we set the normal velocity weakly to zero and specify zero normal heat flux. On the parachute, when we have porosity, we use the conditions described in Appendix C.4, and when we do not have porosity, we use the conditions enforced by Eq. (53). We resolve the flow through the gaps as well as the wider gaps. In dealing with the fabric porosity, we use Eqs. (19) and (23). In Eq. (19), we set D μ equal to the fabric porosity, and drop the β term.  We use quadratic NURBS meshes. The canopy surface has 3432 control points and 1872 elements. The mesh resolution in the radial direction is 4 elements for the ribbons, 2 or 3 elements for the gaps, and 6 elements for the wider gaps. Figure 5 shows, as an example, the mesh for AM31. The mesh information is given in Table 2. In the case of AM21, we remesh at 0.33 s to improve the quality of the mesh. Before the remeshing, the mesh has the same connectivity and function space as the mesh for AM31.
The computations are carried out in two stages, first without porosity, and then with porosity. Both stages have 3 nonlinear iterations per time step, and in the GMRES iterations, we use a nodal-block-diagonal preconditioner. The time-step sizes and the number of GMRES iterations per nonlinear iteration are shown in Table 3. We set γ ACI = −1 in Eq. (49).  Figures 6 and 7 show the velocity magnitude and density for the three cases. We can see that we are resolving the flow through the gaps as well as the wider gaps.

Geometric-porosity modeling
We calculate the porosity coefficients based on the porosity model represented by Eqs. (19) and (23). The subscript G will refer to the geometric porosity, and the area is divided into four patches, as shown in Fig. 8. Figures 9 and 10 show the patch-averaged values of |ṁ| through the gaps and M 2 . In the periods with fabric porosity, these values decrease with increasing altitude for all patches, because the density, velocity and dynamic pressure all decrease with increasing altitude (see Table 1). Figure 11 shows the correlation between the time-averaged values of |ṁ| and M 2 displayed in Figs. 9 and 10, with the intervals used in the time-averaging shown in those two figures. The correlation is established by leastsquares curve fitting based on Eq. (19), which yields μ D and β that represent the geometric porosity of the gaps, and we call those μ D G and β G . Table 4 shows, for each patch, the μ D G and β G values, and the coefficient of determination in the curve fitting.

Modeled porosity
In the full-canopy fluid mechanics computations with the modeled porosity, for each patch, the relationship between M 2 and |ṁ| is given as where the subscript F refers to the fabric porosity, A F is the total fabric area for the patch, A 1 is the fluid surface area for the patch, taken as is set equal to the fabric porosity, taken as 40 CFM here, and β F = 0 here. Figure 12 shows that relationship for all the patches. Table 5 shows, for all the cases and patches, the values of μ D and β used in Eq. (5).

Computation
We do the full-canopy fluid mechanics computation for AM31. The domain is a cylinder with diameter 96.5 ft and height 96.5 ft. The distance between the parachute skirt and the inflow boundary is approximately 33 ft. The boundary conditions are the same as those described in 3.4.1 for the one-gore fluid mechanics computations.
We use a quadratic NURBS mesh. The number of control points and elements for the volume mesh are 270,296 and 190,464. The canopy surface has 6336 control points and 5040 elements. The mesh resolution in the radial direction for the wider gaps is 3 elements.  Figures 13 and 14 show the surface and volume meshes. The computations are carried out in three stages, first without porosity, and the second and the third with porosity. The timestep size for first and second stages is 5.0×10 −5 s, and for the third stage 5.0×10 −4 s. In all three stages, we have 3 nonlinear iterations per time step and 60 GMRES iterations per nonlinear iteration. In the GMRES iterations, we use a nodal-block-diagonal preconditioner. We set γ ACI = −1 in Eq. (49). Figures 15 and 16 show the velocity and density for AM31 after a settled solution is reached.

Concluding remarks
Geometric porosity, a design feature in spacecraft parachutes, is created by the hundreds of gaps and slits that the flow goes through. Accurate geometric-porosity modeling is essential. That is because FSI analysis with resolved geometric porosity would require resolving the flow that goes through the hundreds of gaps and slits as they change their shapes during the computation, which would be exceedingly time-consuming.
The geometric-porosity model introduced earlier in conjunction with the ST computational methods, the HMGP, enabled successful computational analysis and design studies of the Orion spacecraft main parachutes, which operate in the incompressible-flow regime. Recently, porosity models and ST computational methods were introduced, in the con- Fig. 12 Relationship between |ṁ| and M 2 as given by Eq. (5). The lines labeled "Fabric" represent the fabric-porosity part of that relationship (without the ratio AF A1 ), and the curves labeled "Geometric" come from Fig. 11   typical finite element discretization. In the flow analysis, the combination of the ST framework, NURBS basis functions, and the SUPG stabilization assures superior computational accuracy. The ST-SI plays a key role not only by enabling the porosity modeling with its version where the SI is between a thin porous structure and the fluid on its two sides, but also by removing, without loss of accuracy, the matching requirement between the NURBS patches. The porosity coefficients used in modeling the geometric porosity are determined in high-resolution computation of the flow field for a slice of the parachute canopy with the Mach number given and over a range of altitudes. In these high-resolution computations, the flow going though each gap is resolved. After performing those computations for the drogue parachute, using the

A Porosity models
Consider a porous media as shown in Fig. 17. The flow is only in the normal direction, and across the media the mass flow rate is invariant and the pressure is continuous. The temperature-related condition will be described later.
We assume that, compared to the fluxes, the term d dt is negligible. Then, the mass flow rate across the media, is constant. Here, u R is the velocity relative to the porous media, which is only in the normal direction. Fig. 17 Schematic representation of a porous media. The coordinates x B and x A represent the "B" ("below") and "A" ("above") sides of the media. The flow direction from the B side to the A side is taken as the positive flow direction. The flow is only in the normal direction, and across the media the mass flow rate is invariant and the pressure is continuous. The temperature-related condition will be described later We assume that the pressure gradient can be expressed as where S and L are model parameters. This is know as the Darcy-Forchheimer model. We assume a polytropic process in the media: where n is the exponent constant, and C is a constant. This is general enough to cover most processes.

A.1 Relationship between the fluid inside the media and the surrounding fluid
Multiplying Eq. (8) with the density, we obtain and using Eq. (9), we can integrate in the normal direction: With the transformed model parameters defined as the integration yields Substituting for C from Eq. (9), we obtain

A.2 Mass flux
We define M 2 as and because the definition translates to With that, we rewrite Eq. (15) as and this is the equation we solve for |ṁ|. We obtain for β = 0, and the form would be applicable also when β = 0. From that, we can geṫ Remark 1 Setting n = γ gives us M 2 for adiabatic process: Remark 2 Setting n = ∞ gives us M 2 for incompressibleflow process:

A.3 Momentum flux
The force acting on the fluid per unit area due to the media, h M , is expressed in terms of the momentum-conservation fluxes on the two sides:

A.4 Energy flux
Neglecting the energy exchange due to viscous forces, the heat leaving from the fluid to the media, q M , is expressed in terms of the energy-conservation fluxes on the two sides: where the unit normal vector n is pointing from the B side to the A side. Heat flux condition between a thin porous structure and the surrounding fluid specifies q M to a given value. As a special case of that, the adiabatic condition between a thin porous structure and the surrounding fluid specifies q M to zero.

B Compressible-flow ST SUPG method
The compressible-flow ST SUPG method is essentially the same as the compressible-flow DSD/SST method, but without necessarily implying a mesh motion. The compressibleflow DSD/SST method is a straightforward mixture of the DSD/SST concept and the compressible-flow SUPG method. The compressible-flow SUPG method [101][102][103][104] was introduced in 1982 and evolved over the years (see Sect. 1).
In the DSD/SST method, the finite element formulation is written over a sequence of N ST slabs Q n , where Q n is the slice of the ST domain between the time levels t n and t n+1 . The lateral boundary P n will have complementary subsets where essential and natural boundary conditions are enforced, just like how it is with Γ t . At each time step, the integrations are performed over Q n . The functions are continuous within an ST slab, but discontinuous from one ST slab to another, and the superscripts "−" and "+" will indicate the values of the functions just below and just above the time level. The trial solution and test function spaces are defined over Q n by using ST polynomials that are typically first-order, but sometimes higher-order. Each Q n is decomposed into elements Q e n , where e = 1, 2, . . . , (n el ) n . The subscript n used with n el is for the general case where the number of ST elements may change from one ST slab to another.
We assume that we have constructed some suitablydefined finite-dimensional trial solution and test function spaces (S h U ) n and (V h U ) n . The DSD/SST formulation [106][107][108][109][110][113][114][115] of Eq. (3) can be written as follows: given where τ τ τ SUPG is the SUPG stabilization matrix, and ν SHOC is the shock-capturing parameter. The stabilization is residualbased because the residual of the compressible-flow equations, R A (U h ), appears as a factor in the stabilization term. We start with (U h ) − 0 = U 0 (x) and apply the formulation sequentially to all ST slabs Q 0 , Q 1 , Q 2 , . . . , Q N −1 . The stabilization matrix is given from [106,[108][109][110] and [116] as where and τ e SUGN3 −1 = ν e r e r e : G, The symbols θ and ν e represent the temperature and the thermal diffusivity. See Appendix D for the definition of G ST and G. The parameter ν SHOC is defined as where Y is a diagonal scaling matrix constructed from the reference values of the components of U, and In the computations reported here, we use β = 1 and

C.1 ST-SI base version
First we define a new function and introduce a notation based on that: where v is the mesh velocity and v i is its ith component. In the ST-SI method associated with the formulation given by Eq. (27), there will be added boundary terms corresponding to the SI. We will use the labels "Side A" and "Side B" to represents the two sides of the SI. The boundary terms for the two sides will first be added separately, using test functions W h A and W h B . Then, putting together the terms added for each side, the complete set of terms added will be obtained. We give the boundary terms for only Side B: where Here, (P n ) SI is the SI in the ST domain, c is the acoustic speed, and C h is a tensor that will be defined later. Side A counterpart of Eq. (43) can be written by just interchanging subscripts A and B. For the definition of G, see Appendix D.

Remark 3
The first and second integrations set the Euler flux at the boundary to the Lax-Friedrichs flux.

Remark 4
The third integration contains the average viscous terms.

Remark 5
The fourth integration, with γ ACI = 1, is the adjoint consistency term introduced in the symmetricinterior-penalty discontinuous Galerkin method [117]. The other choice is γ ACI = −1, resulting in a method that is adjoint inconsistent, which is known as the nonsymmetricinterior-penalty discontinuous Galerkin method [118].

Remark 6
The fifth integration is a penalty-like term. Several forms of the tensor C h have been proposed and we use the one from [119]: where C is a nondimensional positive constant, which is 1.0 in the computations reported in this article.
Putting together the boundary terms added for each side, the complete set of terms added becomes

C.2 ST-SI version where the SI is a fluid-solid interface with weakly-imposed flow velocity and temperature conditions
The boundary terms added for Side B are given as where and g h and g h θ are given functions.

C.3 ST-SI version where the SI is a fluid-solid interface with weakly-imposed flow velocity and adiabatic conditions
The boundary terms added for Side B are given as where Note that U , j = ∂U ∂ x j .

C.4 ST-SI version where the SI is the interface between a thin porous structure and the surrounding fluid with weakly-imposed flow velocity and adiabatic conditions
In general, the adiabatic condition (q M = 0) between the thin structure and the surrounding fluid implies − κn · ∇ ∇ ∇θ | B = − κn · ∇ ∇ ∇θ | A .
As a special case of that, we might have

C.4.1 Special case
From Eq. (22) with Eq. (23), we obtain the mass flux as a function of the two conservation variables: The boundary terms added for Side B are given as Here the normal components of the Euler flux vectors are taken as where g e (U,ṁ B ) = e − 1 2 u 2 + 1 2 The normal components of the viscous flux vectors, not including the heat conduction flux, are taken as where h T (U) = (I − nn) (n · T) .
The vectors and tensors involved in the fourth and fifth integrations of Eq. (62) are given as

C.4.2 General case
The boundary terms added for Side B are given as where the normal component of the heat conduction part of the viscous flux vectors are taken as The vectors and tensors involved in the fourth and fifth integrations of Eq. (71) are given as

D Element metric tensor
Here provide the element metric tensor in space and in the ST framework from [120], which was based on [116].

D.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: With this vector, we define the element length (see [116]) as where G =Q −TQ−1 .
Remark 8 From this derivation, what we get with D = I has been used in many methods of calculating the stabilization parameters (see, for example, [2]). 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: h MAX ≡ 2 max r (rr : G) − 1 2 .
They are equivalent to and h MAX = 2 min r (rr : G) where λ max and λ min are the maximum and minimum eigenvalues of the argument matrix.

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

D.2 Element metric tensor in the ST framework
The ST Jacobian matrix is where θ is the parametric coordinate in time, and the mesh velocity v is v = ∂x ∂t ξ ξ ξ .
The ST scaling matrix is given as and the scaling becomeŝ The ST metric tensor is defined as