Code-Verification Techniques for Hypersonic Reacting Flows in Thermochemical Nonequilibrium
CCode-Verification Techniques for Hypersonic Reacting Flows inThermochemical Nonequilibrium
Brian A. Freno a , Brian R. Carnes a , V. Gregory Weirs a a Sandia National Laboratories, PO Box 5800, MS 0828, Albuquerque, NM 87185
Abstract
The study of hypersonic flows and their underlying aerothermochemical reactions is particularly importantin the design and analysis of vehicles exiting and reentering Earth’s atmosphere. Computational physicscodes can be employed to simulate these phenomena; however, code verification of these codes is necessaryto certify their credibility. To date, few approaches have been presented for verifying codes that simulatehypersonic flows, especially flows reacting in thermochemical nonequilibrium. In this paper, we present ourcode-verification techniques for verifying the spatial accuracy and thermochemical source term in hypersonicreacting flows in thermochemical nonequilibrium. We demonstrate the effectiveness of these techniques onthe Sandia Parallel Aerodynamics and Reentry Code (
SPARC ). Keywords: code verification, hypersonic flow, thermochemical nonequilibrium, manufactured solutions
1. Introduction
Hypersonic flows are distinguished by flow velocities and stagnation enthalpies that are high enough toinduce chemical reactions and excitation of thermal energy modes. In particular, the time scales of thereactions and thermal excitation are comparable to the characteristic flow time, requiring a fully coupledmodeling approach. The study of hypersonic flows and their underlying aerothermochemical phenomena isparticularly important in the design and analysis of vehicles exiting and reentering Earth’s atmosphere [1, 2].For much of this flight regime, the gas can be modeled as a continuum. The relevant chemical species areeach tracked, such that the interactions within the flow field are accounted for in the evolution of thechemical state. Similarly, considering each internal energy mode allows the different modes to relax towardsequilibrium while accounting for the flow and chemical states. As the flow velocity or stagnation enthalpyincreases, more complex chemical kinetics and internal energy models must be included, and ionization,radiation, and other phenomena may become important.Verification and validation of the implementation and suitability of such models are necessary to developconfidence in the credibility of the simulations. Validation assesses how well the models instantiated inthe code represent the relevant physical phenomena. It is typically performed by comparing simulationresults with experimental results to assess the suitability of the models, the model error, and the practicalbounds of validity of the models. On the other hand, verification assesses the accuracy of the numericalsolutions the code produces, relative to the assumptions and expectations associated with the numericalmethods. Following Roache [3], Salari and Knupp [4], and Oberkampf and Roy [5], verification can be dividedinto code verification and solution verification. Solution verification focuses on numerical-error estimationfor a particular simulation, whereas code verification focuses on the correctness of the numerical-methodimplementation in the code. A review of code and solution verification is presented by Roy [6].This paper focuses on code verification. When numerically solving the underlying equations of theaforementioned models, the equations must be discretized using, for example, finite differences, finite volumes,or finite elements. Due to the finite nature of the discretization, the equations incur a truncation errorand, consequently, their solutions introduce a discretization error. As the discretization is refined, the
Email addresses: [email protected] (Brian A. Freno), [email protected] (Brian R. Carnes), [email protected] (V. Gregory Weirs) a r X i v : . [ phy s i c s . c o m p - ph ] J u l iscretization error should decrease. More rigorously, a code should achieve an expected order of accuracy:as the mesh is refined by a factor, the error should decrease at a rate that is an expected power of thatfactor. In practice, since the exact solution is generally unavailable, manufactured solutions are frequentlyemployed [7].Code verification has been performed on computational physics codes associated with several physicsdisciplines, including fluid dynamics [8, 9, 10, 11], solid mechanics [12], fluid–structure interaction [13], heattransfer in fluid–solid interaction [14], multiphase flows [15], radiation hydrodynamics [16], electrodynam-ics [17], and electromagnetism [18]. Though not as common, code-verification techniques for hypersonic flowshave been presented by Roy et al. [19] for a single-species perfect gas and by Gollan and Jacobs [20] for amulti-species gas in thermal equilibrium.In this paper, we discuss the code-verification techniques we have employed for flows without discontinu-ities. The most noteworthy contribution of this work is our approach to verifying hypersonic reacting flowsin thermochemical nonequilibrium. As the scope of this paper is limited to code verification, subsequentinstances of ‘verification’ are used to abbreviate ‘code verification’.To assess the spatial discretization, we employ manufactured and exact solutions on uniform and nonuni-form meshes and study the convergence of the spatial error norms. Because we consider smooth flows, weemploy the more rigorous L ∞ -norm, of which we demonstrate the greater effectiveness in detecting devi-ations in the expected accuracy. These tests identified lower-order boundary-condition implementations,which reduced the convergence rate and spatial accuracy before being corrected.However, while these techniques are effective for assessing the spatial discretization, they do not directlyreveal errors in the algebraic thermochemical source term. Therefore, to verify these terms, we comparewith an independently developed code for thousands of samples, spanning extreme conditions, with theexpectation that the values agree to machine precision. Additionally, we examine the convergence of thestatistics with respect to the number of samples to assess their sufficiency. This work exposed disagreementdue convergence criteria, as well as errors in model parameters, resulting in corrections to both.This paper is organized as follows. Section 2 describes the governing conservation, energy-exchange, andchemical-kinetics equations. Section 3 details our approach for verifying the spatial discretization throughthe use of manufactured and exact solutions. Section 4 demonstrates the effectiveness of the verificationtechniques for the spatial discretization. Section 5 presents our approach for verifying the thermochemicalsource term. Section 6 demonstrates the effectiveness of the verification techniques for the thermochemicalsource term. Section 7 provides conclusions and an outlook for future work.
2. Governing Equations
In this paper, we consider hypersonic reacting flows in thermochemical nonequilibrium. We make thefollowing approximations: (1) electronic energy is negligible, (2) vibrationally excited molecules can becharacterized by a single vibrational temperature T v , and (3) the translational and rotational temperaturesare in thermal equilibrium: T t = T r = T . With these approximations, we model the conservation of mass, momentum, and energy for a gas inthermochemical nonequilibrium [1, 2]: ∂ U ∂t + ∇ · F c ( U ) = −∇ · F p ( U ) + ∇ · F d ( U ) + S ( U ) , (1)where the conservative-variable state vector U , convective flux F c , pressure flux F p , diffusive flux F d , andthermochemical source term S are U = ρ ρ v ρEρe v , F c ( U ) = ρ v T ρ vv T ρE v T ρe v v T , F p ( U ) = p I p v T T , F d ( U ) = − J τ (cid:0) τ v − q − q v − J T h (cid:1) T (cid:0) − q v − J T e v (cid:1) T , ( U ) = ˙ w0 Q t − v + e Tv ˙ w . ρ = { ρ , . . . , ρ n s } T is the vector of the densities of the n s species, and ρ = P n s s =1 ρ s is the mixture density.For the particular case of air with n s = 5, which is considered later in this paper, the species are N ,O , NO, N, and O. v = { u, v, w } T is the velocity vector, p = P n s s =1 ρ s M s ¯ RT is the pressure, M s is themolecular weight of species s , ¯ R = 8314 .
47 J / kmol / K is the universal gas constant, and T is the translational–rotational temperature. J = [ ρ u , . . . , ρ n s u n s ] T is the diffusion flux of the species, and u s is the diffusionvelocity of species s . τ is the viscous stress tensor, q is the heat flux, q v is the vibrational heat flux, and h = { h , . . . , h n s } T is the vector of enthalpies per mass of the species. ˙ w = { ˙ w , . . . , ˙ w n s } T is the vector ofmass production rates of the species per volume, and Q t − v is the translational–vibrational energy exchange..At high temperatures, vibrational energy is internally excited within molecules, which can cause dissoci-ation or ionization within the gas [2]. The mixture vibrational energy per mass is e v = n s X s =1 ρ s ρ e v s , where e v = (cid:8) e v , . . . , e v ns (cid:9) T is the vector of vibrational energies per mass of the species, such that e v s = (P n vs m =1 e v s,m ( T v ) for molecules,0 for atoms.The vibrational energy per mass of mode m of species s is assumed to have a Boltzmann distribution overthe vibrational energy levels [1, 2]: e v s,m ( T ) = ¯ RM s θ v s,m exp (cid:0) θ v s,m /T (cid:1) − , where n v s is the number of vibrational modes of species s ( n v s = 0 for atoms), and θ v s,m is the characteristicvibrational temperature of mode m of species s . θ v s,m is typically obtained from experiments, and we retrievethe value from a lookup table, such as Table A2.The total energy per mixture mass is E = | v | n s X s =1 ρ s ρ ( c V s T + e v s + h os ) , where c V s is the specific heat at constant volume of species s , which is the sum of the translational androtational contributions. For diatomic molecules, c V s =
52 ¯ RM s and, for atoms, c V s =
32 ¯ RM s . h os is the formationenthalpy of species s . We use Fick’s first law to model the diffusion flux of the species : J = − ρ D ∇ ρ ρ , where D = diag { D , . . . , D n s } contains the effective diffusion coefficients of the species.Assuming the flow is Newtonian and satisfies Stokes’ hypothesis, we model the viscous stress tensor as τ = µ (cid:16) ∇ v + ( ∇ v ) T (cid:17) − µ ( ∇ · v ) I , where µ is the dynamic viscosity.To model the heat fluxes, we employ Fourier’s law: q = − κ ∇ T and q v = − κ v ∇ T v , where κ and κ v are,respectively, the translational–rotational and vibrational heat conductivities of the mixture. µ and κ for the mixture are computed using Wilke’s mixing rule [21].3 .3. Translational–Vibrational Energy Exchange Modeling The translational–vibrational energy exchange is computed using the Landau–Teller model [22]: Q t − v = n s X s =1 ρ s n vs X m =1 e v s,m ( T ) − e v s,m ( T v ) h τ s,m i , (2)where h τ s,m i is the translational–vibrational energy relaxation time for mode m of species s .To compute the relaxation time, we use the semi-empirical approach of Millikan and White [23], whichis extended to account for higher temperatures [24, 25]: h τ s,m i = n s X s =1 y s τ s,m,s ! − + N A n s X s =1 ρ s M s ! σ v s s π ¯ RTM s − , where the mole fraction of species s is y s = ρ s /M s P n s s =1 ρ s /M s , and the translational–vibrational energy relaxation time for mode m of species s when colliding with species s is modeled by τ s,m,s = exp (cid:2) a s,m,s (cid:0) T − / − b s,m,s (cid:1) − . (cid:3) p . (3) p is the pressure expressed in atmospheres. For many gases, a s,m,s and b s,m,s can be modeled by a s,m,s = 1 . × − µ / s,s θ / v s,m ,b s,m,s = 0 . µ / s,s , (4)where µ s,s = M s M s M s + M s is the reduced mass of species s and s . Additionally, the collision-limiting vibrationalcross section is modeled by σ v s = σ v s (cid:18) T (cid:19) , where N A = 6 . × kmol − is the Avogadro constant, and σ v s is the collision-limiting vibrationalcross section at 50,000 K. The mass production rate per volume for species s is modeled by˙ w s = M s n r X r =1 ( β s,r − α s,r ) ( R f r − R b r ) , (5)where α s,r and β s,r are, respectively, the reactant and product stoichiometric coefficients for species s inreaction r . R f r and R b r are the forward and backward reaction rates for reaction r .The reaction rates are defined by R f r = γk f r n s Y s =1 (cid:18) γ ρ s M s (cid:19) α s,r , (6) R b r = γk b r n s Y s =1 (cid:18) γ ρ s M s (cid:19) β s,r , (7)4uch that, in (6) and (7), k f r , k b r , and ρ s /M s are expressed using the centimeter–gram–second system of units(CGS), R f r and R b r are expressed using the meter–kilogram–second system of units, and γ = 1000 cm · kmolm · mol is the conversion factor. The forward and backward reaction-rate coefficients k f r and k b r are modeled usingthe approach of Park [26]: k f r ( T c ) = C f r T η r c exp ( − θ r /T c ) ,k b r ( T ) = k f r ( T ) K e r ( T ) , where C f r and η r are empirical parameters, and θ r is the activation energy of reaction r , divided by theBoltzmann constant. T c is the rate-controlling temperature; it is set to T c = √ T T v for dissociative reactionsand T c = T for exchange reactions. K e r is the equilibrium constant for reaction r , modeled by K e r ( T ) = exp " A r (cid:18) T (cid:19) + A r + A r ln (cid:18) T (cid:19) + A r T + A r (cid:18) T (cid:19) , (8)where A i r are empirical curve-fit coefficients.In this paper, K e r ( T ) is limited to [exp( − , exp(81)], and, when computing ˙ w , T and T v are increasedto 500 K if they are less than 500 K. We limit the scope of this paper to the five-species air model, which consists of N , O , NO, N, and O.These species can undergo the dissociation and exchange reactions listed in Table 1. Additional propertiesof the species and their reactions are provided in Appendix A. r Reaction Type of Reaction1–5 N + M (cid:11) N + N + M , M = { N , O , NO , N , O } Dissociation6–10 O + M (cid:11) O + O + M , M = { N , O , NO , N , O } Dissociation11–15 NO + M (cid:11) N + O + M , M = { N , O , NO , N , O } Dissociation16 N + O (cid:11) N + NO Exchange17 NO + O (cid:11)
N + O Exchange
Table 1: Five-species air model: Reactions.
We note that the five-species air model is just one of a few options. Gimelshein et al. [27, 28] provide acomparison of air models represented by different numbers of species. Nonetheless, the verification techniqueswe propose for the thermochemical source term are applicable to any number of species.
3. Verification Techniques for Spatial Discretization
We begin our approach to code verification by computing the spatial accuracy of the numerical discretiza-tion. To compute the spatial accuracy, we compare the solution to the discretized equations with the solutionto the continuous equations. For each of the flow variables, we compute error norms and compare the ratesat which the error norms decrease with respect to the rates at which the mesh size increases.
In a steady state, the governing system of partial differential equations (1) can be written generally as r ( U ) = , (9)where r is the residual vector, and U = U ( x ) is the state vector. To solve (9) numerically, it must bediscretized: r h ( U h ) = , (10)5here r h is the residual of the discretized system of equations, U h is the solution to the discretized equations,and h describes the size of the mesh employed by the discretization. To simplify the notation in this section,we assume appropriate mappings of the residuals and solutions onto continuous or discrete space.The truncation error is τ h ( V ) = r h ( V ) − r ( V ) . (11)Letting V = U h and adding (9), (11) becomes τ h ( U h ) = r h ( U h ) − r ( U h ) + r ( U ) = r ( U ) − r ( U h ) . When r is linear or linearized with respect to U , the discretization error e h = U h − U is related to thetruncation error by r ( e h ) = − τ h ( U h ) [29, 5].For a p th -order-accurate discretization, the truncation error is τ h = C r h p + O ( h p +1 ) , where C r is a function of derivatives of the state vector but is independent of h . Once the meshes arefine enough that the approximation is within the asymptotic region, h p +1 (cid:28) h p , then τ h ≈ C r h p and e h ≈ C U h p , where C U is also independent of h .The observed accuracy can be computed using two meshes. For example, let the coarser mesh be char-acterized by h . The second mesh, if q -times as fine in each dimension as the first, is characterized by h/q .For each scalar field α in U (e.g., α = { ρ , . . . , ρ n s , u, v, w, T, T v } ), each mesh has a discretization error: e αh ≈ C α h p and e αh/q ≈ C α ( h/q ) p . p can be approximated locally by p ≈ log (cid:12)(cid:12) e αh /e αh/q (cid:12)(cid:12) log q = log q (cid:12)(cid:12) e αh /e αh/q (cid:12)(cid:12) . (12) Computing the spatial accuracy p in (12) requires computing errors e h , which, in turn, require solutions U . Therefore, we employ exact and manufactured solutions. For limited cases, there exist exact solutions U Exact to (9), which can be directly compared with thecomputed solutions U h from (10) with negligible implementation effort. However, available exact solutionsonly span a small subset of the application space, so we require additional approaches to thoroughly test thecapabilities of the code. Manufactured solutions U MS enable us to develop solutions that exercise the features we intend to test.Unless the manufactured solutions are exact solutions, they will not satisfy (9): r ( U MS ) = . Therefore, aforcing term is added to (10) to account for the presence of the manufactured solutions: r h ( U h ) = r ( U MS ) . (13) r ( U MS ) in (13) is computed analytically since r and U MS are known.Because the equations are differential and the error is a function of derivatives of the state vector,the manufactured solutions should be smooth, continuously differentiable functions with generally nonzeroderivatives. Additionally, for the approximation to be in the asymptotic region without requiring especiallyfine meshes, variations over the domain should not be large.6 .3. Error Norms The spatial accuracy p can be computed at a single location in the domain, as is done in (12); however,this approach has two shortcomings: (1) for cell-centered schemes, the cell centroids of a coarser mesh onlycoincide with those arising after mesh refinement in very limited cases, and (2) in regions where the errorsvanish, the computed spatial accuracy is meaningless. Therefore, for each scalar field α , we use error norms ε α to quantify the spatial accuracy p : p = log q (cid:0) ε α h /ε α h/q (cid:1) . (14)Because the error norms are global quantities, they do not require coincident solution locations across meshes,and they will only be zero for trivial flows.We consider two norms in particular:1. L -norm: ε α = k α h ( x ) − α ( x ) k = Z Ω | α h ( x ) − α ( x ) | d Ω . (15)The L -norm enables us to compute the spatial accuracy based on the average error throughout thedomain without significant contamination from localized deviations. Such localized deviations can arisefrom discontinuities, such as shocks, as well as from boundary conditions discretized with lower-orderspatial accuracy than what is used for the domain interior.2. L ∞ -norm: ε ∞ α = k α h ( x ) − α ( x ) k ∞ = max x ∈ Ω | α h ( x ) − α ( x ) | . (16)Unlike the L -norm, the L ∞ -norm catches the aforementioned localized deviations by computing themaximum error throughout the domain. This is particularly useful when such deviations are unex-pected.For the smooth flows considered in this paper, the observed orders of accuracy computed by the L -normand L ∞ -norm are expected to be the same.
4. Spatial-Discretization Verification Results
We demonstrate the aforementioned and forthcoming code-verification techniques on the Sandia ParallelAerodynamics and Reentry Code (
SPARC ) [30, 31, 32, 33] presently being developed at Sandia NationalLaboratories.
SPARC is a compressible computational fluid dynamics code designed to model transonic andhypersonic reacting turbulent flows.
SPARC employs a cell-centered finite-volume discretization, and the simulations presented herein use theSteger–Warming flux-vector splitting scheme. The expectation is that
SPARC is second-order accurate ( p = 2)in space for flows without discontinuities. Therefore, because the midpoint rule has a discretization error of O ( h ), we integrate r ( U MS ) in (13) using the midpoint rule, without decreasing the order of accuracy. Weadditionally use approximations to compute the L -norm in (15) and the L ∞ -norm in (16): ε α ≈ n X i =1 (Ω i | α h ( x i ) − α ( x i ) | ) , (17) ε ∞ α ≈ max ≤ i ≤ n | α h ( x i ) − α ( x i ) | , (18)where x i and Ω i are, respectively, the centroid and volume of cell i , and n is the number of cells. For auniform mesh, (17) is the same as the discrete L -norm, when multiplied by the domain volume Ω. Equation(18) is the discrete L ∞ -norm. SPARC computes r ( U MS ) in (13) through the automatic differentiation tool Sacado [34, 35]. This avoids theneed for using an external package to compute r ( U MS ) through symbolic manipulation. A major drawbackto using symbolic manipulation is the large amount of output generated, which must be copied and formatted7nto a source-code file. For the work in this paper, a single term for one equation would have required severallines of code. On the other hand, the automatic-differentiation approach we employ only requires codingthe terms in the differential equation. By templating the exact solution, we can generate the equivalentsource term using only a few lines of code. The source term is never provided in symbolic form, but it can beefficiently evaluated. The Sacado implementation in SPARC has been unit tested by comparing with examplescomputed using symbolic manipulation. The results of those tests agree within machine precision.For these results, the goal is to ensure second-order accuracy is achievable. Therefore, limiters are disabledand boundary conditions are expected to be second-order accurate. For complex simulations, the accuracymay be deliberately reduced in favor of numerical stability. The implications of which can be assessedthrough solution verification, which is not the focus of this paper.Our verification of
SPARC begins with the simplest tests, and, as these tests are satisfied, we add com-plexity to the subsequent tests. We begin by verifying single-species supersonic inviscid flow, ultimatelyincreasing the complexity to multi-species hypersonic inviscid flow in thermochemical nonequilibrium.For many of our two-dimensional manufactured solutions, we use the following solution structures, orsubsets thereof: ρ N ( x, y ) = ¯ ρ N (cid:2) − (cid:15) sin (cid:0) πx (cid:1) (cid:0) sin (cid:0) πy (cid:1) + cos (cid:0) πy (cid:1)(cid:1)(cid:3) ,ρ O ( x, y ) = ¯ ρ O (cid:2) (cid:15) sin (cid:0) πx (cid:1) (cid:0) sin (cid:0) πy (cid:1) + cos (cid:0) πy (cid:1)(cid:1)(cid:3) ,ρ NO ( x, y ) = ¯ ρ NO (cid:2) (cid:15) sin (cid:0) πx (cid:1) (cid:0) sin (cid:0) πy (cid:1) (cid:1)(cid:3) ,ρ N ( x, y ) = ¯ ρ N (cid:2) (cid:15) sin (cid:0) πx (cid:1) (cid:0) cos (cid:0) πy (cid:1) (cid:1)(cid:3) ,ρ O ( x, y ) = ¯ ρ O (cid:2) (cid:15) sin (cid:0) πx (cid:1) (cid:0) sin (cid:0) πy (cid:1) + cos (cid:0) πy (cid:1)(cid:1)(cid:3) ,u ( x, y ) = ¯ u (cid:2) (cid:15) sin (cid:0) πx (cid:1) (cid:0) sin (cid:0) πy (cid:1) + cos (cid:0) πy (cid:1)(cid:1)(cid:3) ,v ( x, y ) = ¯ v (cid:2) − (cid:15) sin (cid:0) πx (cid:1) (cid:0) sin (cid:0) πy (cid:1) (cid:1)(cid:3) ,T ( x, y ) = ¯ T (cid:2) (cid:15) sin (cid:0) πx (cid:1) (cid:0) sin (cid:0) πy (cid:1) + cos (cid:0) πy (cid:1)(cid:1)(cid:3) ,T v ( x, y ) = ¯ T v (cid:2) (cid:15) sin (cid:0) πx (cid:1) (cid:0) sin (cid:0) πy (cid:1) + cos (cid:0) πy (cid:1)(cid:1)(cid:3) . (19)These solutions are shown in Figure 1 for ( x, y ) ∈ [0 ,
1] m × [0 ,
1] m. The inflow and outflow boundariesare located at x = 0 m and x = 1 m. The velocity field ensures the flow is tangential to the slip-wall(tangent-flow) boundaries located at y = 0 m and y = 1 m. The first set of tests consists of a single species ( n s = 1, ρ = ρ ) inviscid ( F d = 0) flow in thermochemicalequilibrium ( S = , T = T v ). For this type of flow, the reference velocity ¯ u is determined from a referenceMach number ¯ M : ¯ u = ¯ M q γR ¯ T , (20)where γ is the ratio of specific heats for air, and R is the specific gas constant of air. For this test, we simulate a simple, one-dimensional flow using manufactured solutions: ρ ( x ) = ¯ ρ [1 − (cid:15) sin( πx )] ,u ( x ) = ¯ u [1 − (cid:15) sin( πx )] ,T ( x ) = ¯ T [1 + (cid:15) sin( πx )] , with x ∈ [0 ,
1] m, ¯ ρ = 1 kg/m , ¯ M = 2 .
5, ¯ T = 300 K, and (cid:15) = 0 . p is computed from (14) for α = { ρ, u, T } . Five 1Dmeshes are used, consisting of 50, 100, 200, 400, and 800 elements.8 y . . . . . . . . . . . . . . ρ N / ¯ ρ N xy . . . . . . . . . . . . . . ρ O / ¯ ρ O xy . . . . . . . . . . . ρ NO / ¯ ρ NO xy . . . . . . . . . . . ρ N / ¯ ρ N xy . . . . . . . . . . . ρ O / ¯ ρ O xy . . . . . . . . . . u/ ¯ uxy − . − . − . − . − . . . . . . v/ ¯ v xy . . . . . . . . . . . . . . T/ ¯ T xy . . . . . . . . . . . . . . . . . . T v / ¯ T v Figure 1: 2D manufactured solutions from Equation (19) with (cid:15) = 0 . Table 2 shows the observed accuracy from the original state of the code using the L ∞ -norm and L -normof the error. The L -norm indicates second-order accuracy ( p = 2), whereas the L ∞ -norm indicates first-order accuracy ( p = 1). This example demonstrates the usefulness of the L ∞ -norm. The L -norm suggeststhe code, on average, is second-order accurate; however, the L ∞ -norm captures localized deviations in theaccuracy. These deviations are due to the supersonic-inflow and supersonic-outflow boundary-conditionimplementations being only first-order accurate. For this case, the accuracy reduction is limited to thevicinity of the boundaries. L ∞ -norm L -normMesh ρ u T ρ u T Table 2: 1D MMS, n s = 1, T v = T , ˙ w = : Observed accuracy p from original boundary conditions. ∞ -norm L -normMesh ρ u T ρ u T Table 3: 1D MMS, n s = 1, T v = T , ˙ w = : Observed accuracy p from corrected boundary conditions. . . . . . . . . n − − − − − − l og ( ε ∞ α / ¯ α ) , α = { ρ , u , T } ρuT O ( h ) O ( h ) (a) L ∞ -norm . . . . . . . . n − . − . − . − . − . − . − . − . l og ( ε α / ¯ α ) , α = { ρ , u , T } ρuT O ( h ) O ( h ) (b) L -norm Figure 2: 1D MMS, n s = 1, T v = T , ˙ w = : Norms of the error. Dashed lines denote original boundary conditions; solid linesdenote corrected boundary conditions. We corrected the two boundary-condition implementations to be second-order accurate, which is con-firmed in Table 3. Additionally, Figures 2a and 2b show the two error norms for each of the flow variables,before and after correcting the boundary conditions. As shown in Figure 2a, the maximum errors are reducedby orders of magnitude upon correcting the boundary conditions. Furthermore, Figure 2b shows that thecorrect boundary conditions reduce the average error by a factor of approximately three.
This test increases the complexity of the first test by including variations along a second dimension. Inaddition to the supersonic-inflow and supersonic-outflow boundary conditions of the first test, the slip-wallboundary condition is exercised.The manufactured solutions for this case for α = { ρ, u, v, T } are listed in (19) and shown in Figure 1.In (19) and in Figure 1, ρ = ρ N , ¯ ρ = ¯ ρ N , and ¯ v = ¯ u (20). The domain is a square with ( x, y ) ∈ [0 ,
1] m × [0 ,
1] m, and ¯ ρ = 1 kg/m , ¯ M = 2 .
5, ¯ T = 300 K, and (cid:15) = 0 . p is computed from (14). Five 2D meshes are used,consisting of 25 ×
25, 50 ×
50, 100 × × ×
400 elements. These meshes are chosen totest the spatial accuracy of the discretization for nonuniform meshes, and are created using the approachin Appendix B. The 50 ×
50 mesh is shown in Figure 3.Table 4 shows the observed accuracy using the L ∞ -norm and L -norm of the error. Both norms indicatefirst-order accuracy ( p = 1), despite the second-order-accuracy expectation. This inconsistency is due to thesupersonic-inflow, supersonic-outflow, and slip-wall boundary-condition implementations being only first-order accurate. Unlike the first case, the implications of the first-order-accurate boundary conditions areglobal for this case.The corrected boundary-condition implementations are confirmed to be second-order accurate ( p = 2) inTable 5. Figures 4a and 4b show the two error norms for each of the flow variables, before and after correctingthe boundary conditions. The correct boundary conditions reduce both the maximum and average error by10 y Figure 3: 2D MMS, n s = 1, T v = T , ˙ w = : Mesh with 50 ×
50 elements. L ∞ -norm L -normMesh ρ u v T ρ u v T Table 4: 2D MMS, n s = 1, T v = T , ˙ w = : Observed accuracy p from original boundary conditions. L ∞ -norm L -normMesh ρ u v T ρ u v T Table 5: 2D MMS, n s = 1, T v = T , ˙ w = : Observed accuracy p from corrected boundary conditions. orders of magnitude. For the subsequent results, we omit the L -norm and consider only the correctedboundary conditions. 11 . . . . . . . . . √ n − − − − − − − l og ( ε ∞ α / ¯ α ) , α = { ρ , u , v , T } ρuvT O ( h ) O ( h ) (a) L ∞ -norm . . . . . . . . . √ n − − − − − − − l og ( ε α / ¯ α ) , α = { ρ , u , v , T } ρuvT O ( h ) O ( h ) (b) L -norm Figure 4: 2D MMS, n s = 1, T v = T , ˙ w = : Norms of the error. Dashed lines denote original boundary conditions; solid linesdenote corrected boundary conditions. This test exercises the same boundary conditions exercised in Section 4.1.2, but for an exact solutionthat does not require an additional source term. The exact solution is a steady, isentropic vortex, which wesimulate in a quarter-annulus domain [36, 37, 38].The exact solutions for this case are ρ ( r ) = ρ i (cid:20) γ − M i (cid:18) − (cid:16) r i r (cid:17) (cid:19)(cid:21) γ − ,u r ( r ) = 0 ,u θ ( r ) = − a i M i r i r ,T ( r ) = T i (cid:20) γ − M i (cid:18) − (cid:16) r i r (cid:17) (cid:19)(cid:21) , with ρ i = 1, a i = 1, M i = 2 .
25, and T i = 1 / ( γR ). r is the distance from the center of the full annulus, whichis bounded between r i = 1 and r o = 1 . u i = v i = a i M i .Upon solving (10), the observed order of accuracy p is computed from (14) for α = { ρ, u, v, T } . Six 2Dmeshes are used, consisting of 32 ×
8, 64 ×
16, 128 ×
32, 256 ×
64, 512 × ×
256 elements. The64 ×
16 mesh is shown in Figure 6.Table 6 shows the observed accuracy, using the L ∞ -norm of the error, which indicates second-orderaccuracy ( p = 2). Figure 7 shows the L ∞ -norm for each of the flow variables.Mesh ρ u v T Table 6: 2D Exact, n s = 1, T v = T , ˙ w = : Observed accuracy p using L ∞ -norm of the error. y . . . . . . . . . . ρ/ρ i xy . . . . . . . . . . . u/u i xy . . . . . . . . . . . v/u i xy . . . . . . . . . . . T/T i Figure 5: 2D Exact, n s = 1, T v = T , ˙ w = : Exact solutions. xy Figure 6: 2D Exact, n s = 1, T v = T , ˙ w = : Mesh with 64 ×
16 elements. . . . . . . . . . √ n − . − . − . − . − . − . − . − . − . l og ( ε ∞ α / α i ) , α = { ρ , u , v , T } O ( h ) ρuvT Figure 7: 2D Exact, n s = 1, T v = T , ˙ w = : L ∞ -norm of the error. For this test, we consider a three-dimensional flow. The manufactured solutions for this case are ρ ( x, y, z ) = ¯ ρ (cid:2) − (cid:15) sin (cid:0) πx (cid:1) (sin( πy ) + cos( πy ))( sin( πz ) + cos( πz )) (cid:3) ,u ( x, y, z ) = ¯ u (cid:2) (cid:15) sin (cid:0) πx (cid:1) (sin( πy ) + cos( πy ))( sin( πz ) + cos( πz )) (cid:3) ,v ( x, y, z ) = ¯ v (cid:2) − (cid:15) sin (cid:0) πx (cid:1) (sin( πy ) )( sin( πz ) + cos( πz )) (cid:3) ,w ( x, y, z ) = ¯ w (cid:2) − (cid:15) sin (cid:0) πx (cid:1) (sin( πy ) + cos( πy ))( sin( πz ) ) (cid:3) ,T ( x, y, z ) = ¯ T (cid:2) (cid:15) sin (cid:0) πx (cid:1) (sin( πy ) + cos( πy ))( sin( πz ) + cos( πz )) (cid:3) , with ( x, y, z ) ∈ [0 ,
1] m × [0 ,
1] m × [0 ,
1] m, and ¯ ρ = 1 kg/m , ¯ M = 2 .
5, ¯ v = ¯ w = ¯ u (20), ¯ T = 300 K, and (cid:15) = 0 . p is computed from (14) for α = { ρ, u, v, w, T } .Five 3D meshes are used, consisting of 25 × ×
25, 50 × ×
50, 100 × × × × × ×
400 elements. These meshes are chosen to test the spatial accuracy of the discretization fornonuniform meshes, and are created using the approach in Appendix B. The 50 × ×
50 mesh is shown inFigure 8.Table 7 shows the observed accuracy, using the L ∞ -norm of the error, which indicates second-orderaccuracy ( p = 2). Figure 9 shows the L ∞ -norm for each of the flow variables.Mesh ρ u v w T Table 7: 3D MMS, n s = 1, T v = T , ˙ w = : Observed accuracy p using L ∞ -norm of the error. y z Figure 8: 3D MMS, n s = 1, T v = T , ˙ w = : Mesh with 50 × ×
50 elements. . . . . . . . . . √ n − . − . − . − . − . − . − . − . l og ( ε ∞ α / ¯ α ) , α = { ρ , u , v , w , T } O ( h ) ρuvwT Figure 9: 3D MMS, n s = 1, T v = T , ˙ w = : L ∞ -norm of the error. .2. Five-Species Inviscid Flow in Chemical Nonequilibrium This set of tests uses the five-species air model ( n s = 5, ρ = { ρ N , ρ O , ρ NO , ρ N , ρ O } ) mentioned inSection 2.5. The flow is inviscid ( F d = 0) and in chemical nonequilibrium ( ˙ w = ). For this type of flow,the reference velocity ¯ u is determined from a reference Mach number ¯ M :¯ u = ¯ M vuut γ n s X s =1 ρ s ρ ¯ RM s ! ¯ T . (21)While these tests do not directly test the thermochemical-source-term implementation (as described inSection 5), they do test the coupling of the source term with the differential terms, as well as the spatialdiscretizations of multiple species and temperatures.
For this test, the flow is in thermal equilibrium ( T = T v ). The manufactured solutions for this case for α = { ρ N , ρ O , ρ NO , ρ N , ρ O , u, v, T } are listed in (19) and shown in Figure 1.The domain is a square with ( x, y ) ∈ [0 ,
1] m × [0 ,
1] m, and ¯ ρ N = 0 .
77 kg/m , ¯ ρ O = 0 .
20 kg/m ,¯ ρ NO = 0 .
01 kg/m , ¯ ρ N = 0 .
01 kg/m , ¯ ρ O = 0 .
01 kg/m , ¯ M = 2 .
5, ¯ v = ¯ u (21), ¯ T = 3500 K and (cid:15) = 0 . p is computed from (14) for α = { ρ N , ρ O , ρ NO , ρ N ,ρ O , u, v, T } . Seven 2D meshes are used, consisting of 25 ×
25, 50 ×
50, 100 × × × × × L ∞ -norm of the error, which indicates second-orderaccuracy ( p = 2). Figure 10 shows the L ∞ -norm for each of the flow variables.Mesh ρ N ρ O ρ NO ρ N ρ O u v T Table 8: 2D MMS, n s = 5, T v = T , ˙ w = : Observed accuracy p using L ∞ -norm of the error. . . . . . . √ n − − − − − − − l og ( ε ∞ α / ¯ α ) , α = { ρ s , u , v , T } ρ N ρ O ρ NO ρ N ρ O uvT O ( h ) Figure 10: 2D MMS, n s = 5, T v = T , ˙ w = : L ∞ -norm of the error. .2.2. 2D Hypersonic Flow in Thermal Nonequilibrium using a Manufactured Solution For this test, the flow is in thermal nonequilibrium ( T = T v ). The manufactured solutions for this casefor α = { ρ N , ρ O , ρ NO , ρ N , ρ O , u, v, T, T v } are listed in (19) and shown in Figure 1.The domain is a square with ( x, y ) ∈ [0 ,
1] m × [0 ,
1] m, and ¯ ρ N = 0 . , ¯ ρ O = 0 . ,¯ ρ NO = 0 . , ¯ ρ N = 0 . , ¯ ρ O = 0 . , ¯ M = 8, ¯ v = ¯ u (21), ¯ T = 5000 K, ¯ T v = 1000K, and (cid:15) = 0 . p is computed from (14) for α = { ρ N , ρ O , ρ NO , ρ N ,ρ O , u, v, T, T v } . Seven 2D meshes are used, consisting of 25 ×
25, 50 ×
50, 100 × × × × × L ∞ -norm of the error, which indicates second-orderaccuracy ( p = 2). Figure 11 shows the L ∞ -norm for each of the flow variables.Mesh ρ N ρ O ρ NO ρ N ρ O u v T T v Table 9: 2D MMS, n s = 5, T v = T , ˙ w = : Observed accuracy p using L ∞ -norm of the error. . . . . . . √ n − − − − − − − − l og ( ε ∞ α / ¯ α ) , α = { ρ s , u , v , T , T v } ρ N ρ O ρ NO ρ N ρ O uvTT v O ( h ) Figure 11: 2D MMS, n s = 5, T v = T , ˙ w = : L ∞ -norm of the error.
5. Verification Techniques for Thermochemical Source Term
While measuring the spatial accuracy is an effective technique for assessing the discretization, it doesnot directly reveal errors in the coding of the source term S ( U ) in (1). With manufactured solutions, forexample, terms containing derivatives in (1) are evaluated numerically when computing r h ( U h ) in (13) andanalytically when computing r ( U MS ). The thermochemical source term is similarly evaluated on both sidesof (13); however, because the algebraic source-term evaluation does not depend on the spatial discretization,these evaluations use the same source code. Therefore, there are no differences in the evaluation on eachside, and, as a result, errors in the source-term implementation are not detected.To address this limitation, we independently developed a code to compute the terms within the sourceterm, specifically Q t − v ( ρ , T, T v ), e v ( ρ , T, T v ), and ˙ w ( ρ , T, T v ) in (1). We compute these terms for manysamples of { ρ , T, T v } and compare with those obtained from SPARC for a single-cell mesh when initialized to17hose values with no velocity. Additionally, we assess the sufficiency of the number of samples by performingconvergence studies on distribution properties of the values and on the differences between the two codes.Because this approach only computes the thermochemical source term, it can be considered an extensiveunit test. Examples of simpler, traditional unit tests include test computations of the forward reaction-rate coefficient k f r and the equilibrium constant K e r for each type of reaction, the translational–vibrationalenergy relaxation time τ s,m,s for each collision, and the dependent intermediate computations to obtain Q t − v and ˙ w .However, chemical-kinetics models provide individual reaction-rate parameters as a set, with the indi-vidual parameters determined in a manner such that a model best matches experimental data in the regimeof interest. The nonlinear interactions of the individual reactions can lead to much richer model behaviorthan the simple structure of the reaction-rate forms suggests. In addition to detecting errors, our samplingapproach attempts to confirm that the performance of the chemical-kinetics model is not sensitive to theparticular implementation. These remarks apply to the vibrational nonequilibrium model as well, but to alesser extent, as the models are less complicated and model parameters are more accessible to theoreticaldetermination.We have opted for sampling as an expedient approach that can be easily applied to other algebraicterms. However, it may be more efficient to analyze a chemical-kinetics model in detail to identify regions orsurfaces of high sensitivity, and use an adaptive sampling or optimization approach to choose the locationsto compare alternative implementations. Nonetheless, through our approach, we are able to quickly detecterrors.While it may be instinctive to dismiss these techniques as a typically low-rigor code-to-code comparison,we clarify the distinctive and rigorous features.1. This code is independently developed, using the same models and material properties expected to beemployed by SPARC but taken directly from the original references. Alternatively, external softwarecould be used, but, given the variety of published models and material properties, quantifying theagreement and, consequently, assessing the implementation becomes non-trivial.2. Because the models and material properties are the same, when computing the difference in the sourceterms, the required tolerance is tightened from what may typically be a few percent to near machineprecision.3. For numerical solutions to partial differential equations, code-to-code comparisons are typically em-ployed for a few canonical cases, through which it is difficult to identify and isolate errors in thenumerical-method implementation or attribute differences to specific sources. On the other hand, ourcode-to-code comparison targets the portion of the code that manufactured solutions do not assess,and heavily queries conditions covering the thermochemical model’s domain of validity.The effectiveness of these techniques is demonstrated in Section 6.
6. Thermochemical-Source-Term Verification Results
To assess the correctness of the thermochemical-source-term implementation, we generate n S Latin hy-percube samples using the ranges and spacings listed in Table 10, with n S = 2 i , for i = 0 , . . . ,
17. Atthese samples, we query
SPARC and the independent code described in Section 5 to compute Q t − v ( ρ , T, T v ), e v ( ρ , T, T v ), and ˙ w ( ρ , T, T v ) in (1).For each n S , Figures 12a, 13a, and 14a show the minimum, mean, and maximum of the translational–vibrational energy exchange Q t − v , the vibrational energies per mass e v , and the mass production rates pervolume ˙ w , as computed from the independent code. For the vector quantities e v and ˙ w , the elements arepooled. As n S is increased, these values are expected to converge; however, because the samples at each n S are independently determined, the convergence is not monotonic. Convergence of these values suggests thedistribution of the values is sufficiently resolved, and, therefore, the values are sufficiently represented.For n S = 2 = 131,072, Figures 15–17 show the distributions of | Q t − v | , e v , and | ˙ w | , as computed fromthe independent code. n q denotes the number of queries within the ranges on the abscissa. Since Q t − v and˙ w can be non-positive, their ranges, as well as those of e v are listed in Table 11. All of these values varydrastically in magnitude, and, with the exception of e v , in sign.18ariable Minimum Maximum Units Spacing ρ N − kg/m Logarithmic ρ O − kg/m Logarithmic ρ NO − kg/m Logarithmic ρ N − kg/m Logarithmic ρ O − kg/m Logarithmic T
100 15,000 K Linear T v
100 15,000 K Linear
Table 10: S ( U ), n s = 5, T v = T , ˙ w = : Ranges and spacings for Latin hypercube samples of ρ , T , and T v . For every sample, we compute a symmetric relative difference, defined by δ β = 2 | β SPARC − β || β SPARC | + | β | , (22)where β = (cid:8) Q t − v , e v N2 , e v O2 , e v NO , ˙ w N , ˙ w O , ˙ w NO , ˙ w N , ˙ w O (cid:9) , and the prime denotes computation by theindependent code.For each n S , Figures 12b, 13b, and 14b show the maximum relative difference across the samples. As withFigures 12a, 13a, and 14a, the maximum relative difference is expected to increase and converge, throughnot monotonically. Additionally, for n S = 2 , Figures 18a, 19a, and 20a show the relative differences in Q t − v , e v , and ˙ w .As mentioned in Section 5, the relative differences (22) are expected to be near machine precision; however,this is clearly not the case for the red curves in Figures 12b and 13b or the histograms in Figures 18a and19a. As shown in Figure 18a, for approximately 8.7% of the queries, δ Q t − v is greater than 10%, and, for 29%of the queries, δ Q t − v is greater than 1%. Additionally, as shown in Figure 19a, although δ e v is less than 10 − for 99% of the queries, δ e v is greater than 100% for a few of the queries. Even with n S = 1, the red curvein Figure 12b indicates δ Q t − v = 3 . δ e v is within machine precisionthrough n S = 32, but increases by orders of magnitude at n S = 64 and n S = 128. These observationsdemonstrate how thirty-two samples are not enough to represent a seven-dimensional space.These high relative differences were due to two causes.1. The lookup table used by SPARC contained incorrect values for the vibrational constants used in (3)for N and O when the colliding species is NO. These incorrect values introduced an error in Q t − v for all samples.2. The convergence criteria specified in the implementation of Newton’s method used to compute T v from ρe v was loose. Though sufficient for most values of T v , these criteria prove unsuitable for low values.These criteria introduced errors in Q t − v and e v for a few samples. T v also appears in T c for dissociativereactions, but the actual impact on ˙ w was quite small. For a converged, steady problem, however, theoriginal convergence criteria is not expected to affect the final solution.Upon correcting the lookup-table values and tightening the convergence criteria, we reran the SPARC simulations and recomputed δ Q t − v , δ e v , and δ ˙ w , which are shown in Figures 18b, 19b, and 20b for n S = 2 .These results are consistent with our expectations, as all δ Q t − v values are less than 10 − ; all δ e v values areless than 10 − ; and, with the exception of one query, which is slightly greater, all δ ˙ w values are less than10 − .Of the 131,072 δ Q t − v values, the forty-eight greater than 10 − occur when T and T v have a relativedifference of less than 0.2%. As a result, in the numerator of (2), e v s,m ( T ) and e v s,m ( T v ) share many of theleading digits; therefore, precision is lost when computing their difference.Of the 655,360 δ ˙ w values computed from the 131,072 samples, 109 are greater than 10 − . These slightlyelevated differences are a result of the precision loss that can occur from subtraction in (5).The blue curves in Figures 12b, 13b, and 14b show how the maximum relative differences vary withrespect to n S after correcting the lookup-table values and tightening the convergence criteria. As we expect,these curves are generally increasing and converging. The change in T c impacts δ ˙ w only about as much asthe precision loss. 19 n S l og | m ( Q t − v ) | m ( Q t − v ) = min Q t − v m ( Q t − v ) = ¯ Q t − v m ( Q t − v ) = max Q t − v (a) Minimum, mean, and maximum n S − − − − − − − − l og m a x δ Q t − v OriginalCorrected (b) Maximum relative difference
Figure 12: S ( U ), n s = 5, T v = T , ˙ w = : Convergence history of Q t − v sampling. n S − − − − − l og m ( e v ) m ( e v ) = min e v m ( e v ) = ¯ e v m ( e v ) = max e v (a) Minimum, mean, and maximum n S − − − − − − − − l og m a x δ e v OriginalCorrected (b) Maximum relative difference
Figure 13: S ( U ), n s = 5, T v = T , ˙ w = : Convergence history of e v sampling. n S − − − l og | m ( ˙ w ) | m ( ˙ w ) = min ˙ w m ( ˙ w ) = ¯˙ w m ( ˙ w ) = max ˙ w (a) Minimum, mean, and maximum n S − − − − − − − − l og m a x δ ˙ w OriginalCorrected (b) Maximum relative difference
Figure 14: S ( U ), n s = 5, T v = T , ˙ w = : Convergence history of ˙ w sampling. − − − δ Q t − v − l og n q Figure 15: S ( U ), n s = 5, T v = T , ˙ w = : Absolute values of the translational–vibrational energy exchange Q t − v for n S = 2 . . −
10 .. − − − − δ e v − l og n q Figure 16: S ( U ), n s = 5, T v = T , ˙ w = : Vibrational energies per mass e v for n S = 2 . . − − − δ ˙ w − l og n q Figure 17: S ( U ), n s = 5, T v = T , ˙ w = : Absolute values of the mass production rates per volume ˙ w for n S = 2 . −
10 .. − − − − − − − − − δ Q t − v − l og n q (a) Original . −
16 .. −
15 .. −
14 .. −
13 .. −
12 .. −
11 .. −
10 .log δ Q t − v − l og n q (b) Corrected Figure 18: S ( U ), n s = 5, T v = T , ˙ w = : Relative differences in the translational–vibrational energy exchange Q t − v for n S = 2 , using the original lookup table and convergence criteria (a), and the corrected lookup table with tighter convergencecriteria (b). Queries with δ Q t − v = 0 are placed in the lowest bin. . −
16 .. −
14 .. −
12 .. −
10 .. − − − − δ e v − l og n q (a) Original . − . − . − . − . − . δ e v − l og n q (b) Corrected Figure 19: S ( U ), n s = 5, T v = T , ˙ w = : Relative differences in the vibrational energies per mass e v for n S = 2 , using theoriginal convergence criteria (a), and the tighter convergence criteria (b). Queries with δ e v = 0 are placed in the lowest bin. . −
16 .. −
15 .. −
14 .. −
13 .. −
12 .. −
11 .. −
10 .log δ ˙ w − l og n q (a) Original . −
16 .. −
15 .. −
14 .. −
13 .. −
12 .. −
11 .. −
10 .. − δ ˙ w − l og n q (b) Corrected Figure 20: S ( U ), n s = 5, T v = T , ˙ w = : Relative differences in the mass production rates per volume ˙ w for n S = 2 , usingthe original convergence criteria (a), and the tighter convergence criteria (b). Queries with δ ˙ w = 0 are placed in the lowest bin. Q t − v − . × . × J/m /se v N2 . × − . × J/kg e v O2 . × − . × J/kg e v NO . × − . × J/kg˙ w N − . × . × kg/m /s ˙ w O − . × . × kg/m /s ˙ w NO − . × . × kg/m /s ˙ w N − . × . × kg/m /s ˙ w O − . × . × kg/m /s Table 11: S ( U ), n s = 5, T v = T , ˙ w = : Ranges for Q t − v , e v , and ˙ w for n S = 2 . Though the incorrect lookup-table values introduced high relative differences immediately, the conver-gence criteria required more samples to detect their looseness. The convergence of the distribution propertiesin Figures 12a, 13a, and 14a and the convergence and bounds of the errors in Figures 12b, 13b, and 14b pro-vide insight into how many samples are enough to assess the correctness of the thermochemical-source-termimplementation over the ranges listed in Table 10.To determine the impact of these modifications, we reran a high-enthalpy (20 MJ/kg), hypersonic, lam-inar double-cone flow case after correcting the lookup table values and tightening the convergence criteria.Relative to the original values and convergence criteria, we observed up to a 1.4% and 2.7% change in thepressure and heat flux, respectively, on the surface of the body.
7. Conclusions
In this paper, we presented our code-verification techniques for hypersonic reacting flows in thermochem-ical nonequilibrium.To assess the spatial accuracy, we employed manufactured and exact solutions with the L ∞ -norm and L -norm. These approaches revealed the impact of the lower-order boundary conditions.To assess the algebraic thermochemical source term, we queried an independent code, with the expectationthat agreement be near machine precision. We queried the independent code many times and studied theconvergence properties. This approach revealed the impact of erroneous lookup table entries and insufficientlytight convergence criteria.While the scope of this paper has been limited to flows in vibrational nonequilibrium with five speciesthat undergo dissociation and exchange reactions, these techniques could be analogously extended to addressmore complex flows in rotational and electronic nonequilibrium that contain additional species capable ofundergoing ionization reactions. However, it is important to perform a convergence study to determinewhether the number of samples sufficiently spans the ranges of the dependencies. Acknowledgments
The authors thank Derek Dinzl, Travis Fisher, Micah Howard, and Ross Wagnild, for their valuableassistance with
SPARC and the underlying models and properties. This paper describes objective technicalresults and analysis. Any subjective views or opinions that might be expressed in the paper do not nec-essarily represent the views of the U.S. Department of Energy or the United States Government. SandiaNational Laboratories is a multimission laboratory managed and operated by National Technology and En-gineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S.Department of Energy’s National Nuclear Security Administration under contract DE-NA-0003525.23 ppendix A. Five-Species Air Model Properties
The molecular weights M s and formation enthalpies h os for the five species are tabulated in Table A1. s Species M s [g/mol] h os , 0 K [J/kg]1 N . × . × . × Table A1: Five-species air model: Molecular weights M s and formation enthalpies h os [39]. The characteristic vibrational temperatures θ v s,m [39] and collision-limiting cross sections at 50,000 K σ v s [24, 25] are tabulated in Table A2. s Species n v s θ v s,m [K] σ v s [m ]1 N × − × − × − Table A2: Five-species air model: Characteristic vibrational temperatures θ v s,m [39] and collision-limiting cross sections at50,000 K σ v s [24, 25]. The vibrational constants a s,m,s and b s,m,s for Equation (3), associated with the reactions listed inTable 1, are listed in Table A3. N O NO s a s,m,s b s,m,s a s,m,s b s,m,s a s,m,s b s,m,s N Equation (4) Equation (4) 49.5 0.0420O Equation (4) Equation (4) 49.5 0.0420NO Equation (4) Equation (4) 49.5 0.0420N Equation (4) 72.4 0.0150 49.5 0.0420O 72.4 0.0150 47.7 0.0590 49.5 0.0420
Table A3: Five-species air model: Vibrational constants a s,m,s and b s,m,s [24]. α s,r and β s,r associated with the reactions listed in Table 1 are listed inTable A4. α s,r β s,r r N O NO N O N O NO N O1 2 1 22 1 1 1 23 1 1 1 24 1 1 35 1 1 2 16 1 1 1 27 2 1 28 1 1 1 29 1 1 1 210 1 1 311 1 1 1 1 112 1 1 1 1 113 2 1 1 114 1 1 2 115 1 1 1 216 1 1 1 117 1 1 1 1
Table A4: Five-species air model: Stoichiometric coefficients α s,r and β s,r . The reaction-rate-coefficient dependencies C f r , η r , and θ r associated with the reactions listed in Table 1are listed in Table A5. r C f r [CGS] η r θ r [K]1 7 . × − . . × − . . × − . . × − . . × − . . × − . . × − . . × − . . × − . . × − . . × . . × . . × . . × . . × . . × − . . × . Table A5: Five-species air model: Reaction-rate-coefficient dependencies C f r , η r , and θ r [26]. r A r A r A r A r A r . . . − . − . . . . − . . . . . − . − . . . . − . . − . − . − . − . − . Table A6: Five-species air model: Equilibrium-constant coefficients A r , A r , A r , A r , and A r [26]. Appendix B. Approach to Creating Nonuniform Meshes
To generate the nonuniform meshes in Figures 3 and 8, we employ the following nonlinear transformationfrom ( ξ, η, ζ ) ∈ [0 , × [0 , × [0 ,
1] to ( x, y, z ) ∈ [0 , L x ] × [0 , L y ] × [0 , L z ]: x = L ξ + A ˆ ξ , where x = { x, y, z } T , L = diag { L x , L y , L z } , ξ = { ξ, η, ζ } T , ˆ ξ = L (cid:0) − ξ (cid:1) , and A = α ( s ξ,η + s ξ,ζ ) − βs ξ,η βs ξ,ζ βs η,ξ α ( s η,ξ + s η,ζ ) − βs η,ζ − βs ζ,ξ βs ζ,η α ( s ζ,ξ + s ζ,η ) , with α = 1 − cos π , β = sin π , and s ξ ,η = sin( πξ ) sin( πη ). For the mesh in Figure 3, ζ = z = 0. References [1] P. A. Gnoffo, R. N. Gupta, J. L. Shinn, Conservation equations and physical models for hypersonic airflows in thermal and chemical nonequilibrium, Tech. Rep. NASA-TP-2867, NASA Langley ResearchCenter (1989).[2] J. D. Anderson, Jr., Hypersonic and High-Temperature Gas Dynamics, 2nd Edition, American Instituteof Aeronautics and Astronautics, 2006. doi:10.2514/4.861956 .[3] P. J. Roache, Verification and Validation in Computational Science and Engineering, Hermosa Publish-ers, 1998.[4] K. Salari, P. Knupp, Code verification by the method of manufactured solutions, Sandia ReportSAND2000-1444, Sandia National Laboratories (Jun. 2000). doi:10.2172/759450 .[5] W. L. Oberkampf, C. J. Roy, Verification and Validation in Scientific Computing, Cambridge UniversityPress, 2010. doi:10.1017/cbo9780511760396 .[6] C. J. Roy, Review of code and solution verification procedures for computational simulation, Journal ofComputational Physics 205 (1) (2005) 131–156. doi:10.1016/j.jcp.2004.10.036 .[7] P. J. Roache, Code verification by the method of manufactured solutions, Journal of Fluids Engineering124 (1) (2001) 4–10. doi:10.1115/1.1436090 .[8] C. J. Roy, C. C. Nelson, T. M. Smith, C. C. Ober, Verification of Euler/Navier–Stokes codes using themethod of manufactured solutions, International Journal for Numerical Methods in Fluids 44 (6) (2004)599–620. doi:10.1002/fld.660 . 269] R. B. Bond, C. C. Ober, P. M. Knupp, S. W. Bova, Manufactured solution for computational fluiddynamics boundary condition verification, AIAA Journal 45 (9) (2007) 2224–2236. doi:10.2514/1.28099 .[10] S. Veluri, C. Roy, E. Luke, Comprehensive code verification for an unstructured finite volume CFDcode, in: 48th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and AerospaceExposition, American Institute of Aeronautics and Astronautics, 2010. doi:10.2514/6.2010-127 .[11] T. Oliver, K. Estacio-Hiroms, N. Malaya, G. Carey, Manufactured solutions for the Favre-averagedNavier–Stokes equations with eddy-viscosity turbulence models, in: 50th AIAA Aerospace SciencesMeeting including the New Horizons Forum and Aerospace Exposition, American Institute of Aeronau-tics and Astronautics, 2012. doi:10.2514/6.2012-80 .[12] É. Chamberland, A. Fortin, M. Fortin, Comparison of the performance of some finite element dis-cretizations for large deformation elasticity problems, Computers & Structures 88 (11) (2010) 664 –673. doi:10.1016/j.compstruc.2010.02.007 .[13] S. Étienne, A. Garon, D. Pelletier, Some manufactured solutions for verification of fluid–structureinteraction codes, Computers & Structures 106-107 (2012) 56–67. doi:10.1016/j.compstruc.2012.04.006 .[14] A. Veeraragavan, J. Beri, R. J. Gollan, Use of the method of manufactured solutions for the verificationof conjugate heat transfer solvers, Journal of Computational Physics 307 (2016) 308–320. doi:10.1016/j.jcp.2015.12.004 .[15] P. T. Brady, M. Herrmann, J. M. Lopez, Code verification for finite volume multiphase scalar equationsusing the method of manufactured solutions, Journal of Computational Physics 231 (7) (2012) 2924–2944. doi:10.1016/j.jcp.2011.12.040 .[16] R. G. McClarren, R. B. Lowrie, Manufactured solutions for the p radiation-hydrodynamics equations,Journal of Quantitative Spectroscopy and Radiative Transfer 109 (15) (2008) 2590–2602. doi:10.1016/j.jqsrt.2008.06.003 .[17] J. R. Ellis, C. D. Hall, Model development and code verification for simulation of electrodynamic tethersystem, Journal of Guidance, Control, and Dynamics 32 (6) (2009) 1713–1722. doi:10.2514/1.44638 .[18] R. G. Marchand, The method of manufactured solutions for the verification of computational electro-magnetic codes, phdthesis, Stellenbosch (Mar. 2013).[19] C. J. Roy, M. A. McWherter-Payne, W. L. Oberkampf, Verification and validation for laminar hypersonicflowfields, part 1: Verification, AIAA Journal 41 (10) (2003) 1934–1943. doi:10.2514/2.1909 .[20] R. Gollan, P. Jacobs, About the formulation, verification and validation of the hypersonic flow solverEilmer, International Journal for Numerical Methods in Fluids 73 (1) (2013) 19–57. doi:10.1002/fld.3790 .[21] C. R. Wilke, A viscosity equation for gas mixtures, The Journal of Chemical Physics 18 (4) (1950)517–519. doi:10.1063/1.1747673 .[22] L. Landau, E. Teller, Zur theorie der schalldispersion, Physikalische Zeitschrift der Sowjetunion 10 (34)(1936).[23] R. C. Millikan, D. R. White, Systematics of vibrational relaxation, Journal of Chemical Physics 39 (12)(1963) 3209–3213. doi:10.1063/1.1734182 .[24] C. Park, Review of chemical-kinetic problems of future NASA missions, I: Earth entries, Journal ofThermophysics and Heat Transfer 7 (3) (1993) 385–398. doi:10.2514/3.431 .2725] C. Park, J. T. Howe, R. L. Jaffe, G. V. Candler, Review of chemical-kinetic problems of future NASAmissions, II: Mars entries, Journal of Thermophysics and Heat Transfer 8 (1) (1994) 9–23. doi:10.2514/3.496 .[26] C. Park, Nonequilibrium Hypersonic Aerodynamics, John Wiley & Sons, Inc., 1990.[27] S. F. Gimelshein, I. J. Wysong, Impact of the ionization reaction set in nonequilibrium hypersonic airflows, AIAA Journal 58 (3) (2019) 1255–1265. doi:10.2514/1.j058895 .[28] S. F. Gimelshein, I. J. Wysong, Applicability of 5, 7, and 11 species air models in nonequilibriumhypersonic reacting flows, AIAA SciTech Forum, American Institute of Aeronautics and Astronautics,2020. doi:10.2514/6.2020-2190 .[29] J. H. Ferziger, M. Perić, Computational Methods for Fluid Dynamics, Springer, 2002. doi:10.1007/978-3-642-56026-2 .[30] M. Howard, A. Bradley, S. W. Bova, J. Overfelt, R. Wagnild, D. Dinzl, M. Hoemmen, A. Klinvex,Towards performance portability in a compressible CFD code, in: 23rd AIAA Computational FluidDynamics Conference, American Institute of Aeronautics and Astronautics, 2017. doi:10.2514/6.2017-4407 .[31] B. Carnes, V. G. Weirs, T. Smith, Code verification and numerical error estimation for use in modelvalidation of laminar, hypersonic double-cone flows, in: AIAA SciTech 2019 Forum, American Instituteof Aeronautics and Astronautics, 2019. doi:10.2514/6.2019-2175 .[32] S. L. Kieweg, J. Ray, V. G. Weirs, B. Carnes, D. Dinzl, B. Freno, M. Howard, E. Phipps, W. Rider,T. Smith, Validation assessment of hypersonic double-cone flow simulations using uncertainty quantifi-cation, sensitivity analysis, and validation metrics, in: AIAA SciTech 2019 Forum, American Instituteof Aeronautics and Astronautics, 2019. doi:10.2514/6.2019-2278 .[33] J. Ray, S. L. Kieweg, D. Dinzl, B. Carnes, V. G. Weirs, B. Freno, M. Howard, T. Smith, I. Nompelis,G. V. Candler, Estimation of inflow uncertainties in laminar hypersonic double-cone experiments, in:AIAA SciTech 2019 Forum, American Institute of Aeronautics and Astronautics, 2019. doi:10.2514/6.2019-2279 .[34] E. Phipps, R. Pawlowski, Efficient expression templates for operator overloading-based automatic dif-ferentiation, in: Recent Advances in Algorithmic Differentiation, Springer Berlin Heidelberg, Berlin,Heidelberg, 2012, pp. 309–319. doi:10.1007/978-3-642-30023-3_28 .[35] R. A. Bartlett, D. M. Gay, E. T. Phipps, Automatic differentiation of C++ codes for large-scale scientificcomputing, in: Computational Science – ICCS 2006, Springer Berlin Heidelberg, Berlin, Heidelberg,2006, pp. 525–532. doi:10.1007/11758549_73 .[36] M. Aftosmis, D. Gaitonde, T. S. Tavares, Behavior of linear reconstruction techniques on unstructuredmeshes, AIAA Journal 33 (11) (1995) 2038–2049. doi:10.2514/3.12945 .[37] H. Luo, J. Baum, R. Lohner, An improved finite volume scheme for compressible flows on unstructuredgrids, in: 33rd AIAA Aerospace Sciences Meeting and Exhibit, American Institute of Aeronautics andAstronautics, 1995. doi:10.2514/6.1995-348 .[38] L. Krivodonova, M. Berger, High-order accurate implementation of solid wall boundary conditions incurved geometries, Journal of Computational Physics 211 (2) (2006) 492–512. doi:10.1016/j.jcp.2005.05.029doi:10.1016/j.jcp.2005.05.029