This work relies on a compressible biglobal stability approach to describe the wave dynamics in a planar rocket chamber modeled as a porous channel. At first, the effectiveness of the compressible formulation is demonstrated by reproducing, in the absence of a mean flow, the Helmholtz frequencies and mode shapes. Next, the unsteady vorticity fluctuations, which intensify near the walls, are shown to be consistent with those associated with parietal vortex shedding. In this context, the penetration depth of vorticoacoustic waves is found to be strongly dependent on the penetration number. The latter gauges the cubic power of the injection speed to the product of kinematic viscosity, chamber half-height, and frequency squared. As for the strictly hydrodynamic modes, they seem to develop at the porous walls and grow in the core region, where the mean flow velocity is most appreciable. The ensuing modal analysis enables us to predict both longitudinal and transverse modes for several test cases, thus illustrating the tendency of hydrodynamic modes to intensify at higher injection speeds and longer chambers. Furthermore, by repeating the analysis with an active mean flow, one finds that successive increases in the injection speed gradually reduce the predicted frequencies relative to the eigenmodes obtained in a quiescent medium. Finally, recognizing that the spectral analysis is capable of recovering both longitudinal and transverse modes induced by acoustic and hydrodynamic disturbances, their coupled interactions, which often lead to specifically amplified frequencies in static tests, are robustly captured, namely, without resorting to any particular wave decomposition.

## I. INTRODUCTION

Flow instability is commonly observed in various combustors, rockets, and gas turbine engines, where it is often referred to as “resonant combustion” or “combustion instability.” Characterized by internal pressure oscillations that appear at characteristic chamber frequencies, combustion instability is not exclusively associated with the instability of the combustion process itself. In fact, several equally influential sources exist that can support the transfer of energy into resonant modes, and these include fluctuations in mass, momentum and energy fluxes, vorticity generation, unsteady heat transfer, entropy production, and others.^{1} Naturally, one of the principal objectives of combustion stability quantification has been the understanding of the various mechanisms that trigger and drive these self-excited oscillations. Chief among the underlying mechanisms, the study of hydrodynamic instability and its bearing on chamber acoustics has often occupied the central stage.^{2}

On this note, one can refer to Varapaev and Yagodkin,^{3} who may have been the first to introduce a one-dimensional streamfunction formulation (SFF) to investigate the stability of porous channel flows that mimic the bulk gaseous motion in rocket motors. Their work is followed by several notable studies that may be worth enumerating. These include Sviridenkov and Yagodkin^{4} and Beddini,^{5} who explore the mean flow breakdown and turbularization in an idealized rocket chamber configuration; it is then followed by the works of Casalis *et al.*,^{6} who introduce the local-non-parallel (LNP) approach using a primitive variable formulation (PVF); Griffond *et al.*,^{7} who focus on the incompressible Taylor–Culick flow stability in a semi-infinitely long porous cylinder;^{8} and Griffond and Casalis,^{9,10} who consider the planar counterpart, namely, the stability of the incompressible Taylor motion in a semi-infinitely long porous channel.^{11} Complementary studies include those by Féraille and Casalis,^{12} who take into account the effect of particle entrainment on the porous channel flow stability, Abu-Irshaid *et al.*,^{13} who examine the effect of headwall injection on the incompressible Taylor–Culick profile, and those by Chedevergne *et al.*,^{14} who introduce the incompressible biglobal stability approach. The latter can be viewed as being a major development wherein the LNP's spatial emphasis on tracking local, incompressible, and slightly viscous disturbances is abandoned in favor of a temporal scheme that can be globally applied. In this vein, and using a vorticity-streamfunction technique, Chedevergne *et al.*^{15} are able to confirm that direct numerical simulations of the incompressible Taylor–Culick motion instability agree reasonably well with the oscillatory wave structure obtained by consolidating the biglobal stability predictions of the incompressible hydrodynamic waves with those determined analytically from a vorticoacoustic formulation.^{16} Follow-on studies by Casalis and co-workers^{17–19} continue this line of inquiry and connect the onset of parietal vortex shedding (PVS) to the presence of simple wall defects in realistic motor cases. Other investigations of PVS coupling include those by Chakravarthy and Chakraborty^{20} and Li *et al.*^{21} In the interim, Bouyges *et al.*^{22} manage to characterize the impact of small deviations from a circular-port, such as those associated with star-shaped propellant grains, on the resulting flow motion and stability. Along similar lines, Elliott and Gibson^{23} examine the effect of swirl on mean flow stability whereas Li *et al.*^{24} investigate the effect of grid refinement on the stability of wave structures in a long cylindrical chamber with porous walls. In the case of a chamber with non-porous walls, Yeddula and Morgans^{25} and Shelton and Majdalani^{26} develop models to describe the effects of arbitrary varying cross-sectional areas and temperature profiles on the acoustic mode shapes; the latter is further extended to capture the unsteady pressure field in a Rijke tube.

However, most of the aforementioned studies resort to a linear stability framework in which the incompressible Navier–Stokes equations are perturbed. As a result, the hydrodynamic modes obtained from the corresponding eigensolvers do not return the chamber's acoustic frequencies. To overcome this limitation, the present article will focus on the development of a fully compressible biglobal stability framework that can be applied, using Cartesian coordinates, to the compressible flow problem in a porous channel. The latter admits a closed-form analytical solution that can be used as the mean flow.^{27} One advantage of using a compressible stability formulation stands in its ability to predict the coupled hydrodynamic and acoustic response with no need for unsteady flow decomposition at the forefront of the analysis.^{28}

In fact, most stability analyses subdivide the unsteady flow disturbances into three oscillation modes that owe their origins to acoustic, vortical, and entropy-based sources.^{29} According to the Helmholtz–Hodge decomposition technique, these distinct fields can be associated with irrotational, solenoidal, and harmonic oscillations, with the latter being both curl-free and divergence-free.^{30} For more detail on this subject, Ewert and Schröder^{1} provide a judicious overview of the different decomposition techniques and their relations to sound prediction.

In applying the Helmholtz–Hodge decomposition technique to a fluid disturbance, Chu and Kovàsznay^{29} introduce a subdivision of the resulting wave structures into three types: pressure, vorticity, and entropy waves. While the pressure component describes the production, propagation, and absorption of acoustic waves, the vorticity and entropy waves remain connected to the production, convection, and dissipation of rotational and thermal fluctuations, respectively. The acoustic wave, typically represented by a traveling pressure wave, appears most prominently at resonant frequencies. As for the vorticity and entropy waves, they appear to be complementary in nature because of their inherent coupling to the acoustic waves either through the unsteady boundary layer or through the unsteady combustion process. These two waves act as pathways for energy to enter or exit the resonant acoustic modes. Self-excitation can subsequently occur when vorticity and entropy waves couple with the natural acoustics of the chamber to the extent of modulating the ongoing energy transfer. On this note, a Rijke tube can be referred to as a classic representation of a thermoacoustic device where resonance occurs as a by-product of entropy and acoustic waves.^{26} Fundamentally, a Rijke tube consists of a hollow vertical cylinder with a heat source positioned in its lower half. In this case, acoustic waves are induced by amplified thermal fluctuations over a wide range of frequencies, including those that excite the natural frequencies of the enclosure. The resulting instability waves are termed “thermoacoustic.” In a similar manner, vorticity waves can couple with the acoustics to generate “vorticoacoustic” waves. The latter can be triggered by the production and convection of vorticity. It can, therefore, be seen that PVS waves, which are generated hydrodynamically, are not contingent on the presence of a chemically reactive medium.^{18,19} Yet, of the variety of waves that can develop in a rocket chamber, our emphasis here will be placed on the hydrodynamic response of the compressible mean flow motion, with special attention being given to the much discussed vorticoacoustic coupling mechanism.

Among the earliest analytical studies associated with vorticoacoustic instabilities in injection-driven flowfields, one may cite Hart and McClure^{31,32} and Culick,^{8} who introduce expressions for the energy accumulation due to unsteady wave motion in chambers with sidewall injection. In a series of successive studies, Majdalani and Roh^{33} and Fabignon *et al.*^{34} manage to provide continual refinements to the vorticoacoustic formulations developed for geometrically simple configurations. Along similar lines, Brownlee^{35} conducts several illuminating experiments intended at measuring the effects of acoustic oscillations on the burning rate. He notes that longitudinal waves, which remain parallel to the injecting surface, can have a substantial bearing on the production of vorticity. Price^{36} explains that the key behind this mechanism can be traced back to the vortical and turbulent flow processes, which can be affected by velocity and pressure fluctuations near the wall. As a matter of fact, elucidating the complex coupling between unsteady vortical and acoustic waves remains a fundamental objective in the combustion stability community.

One way of accounting for the interactions between the shearing and the compressing processes is to enforce the no-slip boundary condition at the burning surface. Through experimental studies, Brown *et al.*^{37} and Dunlap *et al.*^{38} confirm the oscillatory behavior of the pressure and velocity waves in a porous chamber. It becomes evident, at an early stage, that the inviscid acoustic field, alone, cannot satisfy the velocity-adherence requirement. Vorticity must be generated at the wall to counteract the sweeping motion of the acoustic wave, thus giving rise to unsteady vorticity. Manifested at resonant states, this process controls the manner by which the acoustic field modulates the vorticity generation to draw energy from the main stream. Vuillot and Avalon^{39} use this technique to study the development of the oscillatory Stokes-like boundary layer above an injecting porous wall. Working on the same problem, Majdalani^{40–42} identifies a key parameter that controls the penetration depth of the vorticoacoustic coupling region. To maintain a constant penetration depth,^{43} this parameter provides the necessary balance between convection and diffusion, the two mechanisms responsible for the transport of the vorticoacoustic boundary layer, more properly termed “penetration depth.” Furthermore, Majdalani^{40–42} evaluates the sensitivity of the vortical wave to the mean flow injection pattern using asymptotic tools.

Although most of the original formulations have been motivated by the need to understand solid rocket motor stability, liquid rocket engines have also exhibited resonant combustion that can be modeled using similar approaches.^{44–49} Owing to their geometric configuration and direction of injection, transverse acoustic waves tend to produce steep-fronted shearing amplitudes that can be detrimental to engine hardware. The term “transverse” refers in this context to the direction normal to the axis of the engine; the undesirable emergence of such waves in the F-1/Saturn V engine remains one of the most cited cases of combustion instability. This engine suffered from intense transverse oscillations, requiring over 2700 full scale firings and redesigns, including the addition of multiple forward baffles, to finally reduce the severity of reported oscillations to manageable levels.^{50} In much the same way, Clayton and co-workers^{51–53} reported steep-fronted pressure disturbances of high amplitude tangential oscillations in their liquid rocket testing investigation. Their experimental studies recorded pressure spikes larger than the mean chamber pressure by one order of magnitude. In the same vein, Ando *et al.*^{54} used numerical techniques to study the transverse wave evolution in a pulse detonation engine. Their results showed peaks where the transverse waves collided.

It may be safe to say that only a few theoretical studies have directly targeted the description of transverse vorticoacoustic waves, with the exception of Maslen and Moore,^{55} Crocco,^{56} and several of their co-workers. Most research in this area has remained empirically based. Using a simplified engine configuration, Haddad and Majdalani^{57,58} have provided analytical approximations for the transverse oscillations in the presence of simple mean flow profiles. Their methods have relied on asymptotic techniques that leverage the flow injection Mach number as a perturbation parameter. Along similar lines, Kovacic *et al.*^{59} have managed to extract the acoustic transverse modes numerically, thus extending the validity of previous approaches to include high injection speeds and arbitrary chamber configurations. Their results have also compared favorably with the limited asymptotic solutions developed by Haddad and Majdalani.^{57,58}

In this work, a compressible biglobal stability formulation will be constructed in Cartesian coordinates and then applied to the porous channel flow problem assuming two-dimensional planar conditions. The resulting eigenvalue problem will be shown to provide accurate predictions of the oscillatory frequencies and mode shapes associated with the two-dimensional Taylor motion, thus helping to reconcile between hydrodynamic stability projections and other findings in the literature. These include the ability to reproduce “in one swoop” the rich structures that accompany vorticoacoustic wave motion, particularly, with no need for unsteady flow decomposition and later reconstruction. As for the nature of the oscillations themselves, we seek to capture those that evolve at steady-state operation, notwithstanding initial flow transients.

The paper is organized as follows. In Sec. II, the compressible biglobal framework is formulated, starting with a detailed description of the geometry, underlying assumptions, boundary conditions, and normalizing variables. The linearized forms of the Navier–Stokes, energy, and state equations, which prescribe the unsteady disturbances, are also provided, along with a description of the mean flow profile for the compressible Taylor motion in a porous channel. After reproducing the biglobal stability equations in planar geometry, the numerical approach used to solve the resulting eigenvalue problem is presented in Sec. III. This includes a discussion of the spectral method that we choose to implement in concert with an *N*-point Chebyshev discretization scheme. In Sec. IV, a close examination of the eigensolutions for the vorticoacoustic waves is undertaken both with and without a mean flowfield. In the absence of a mean flow, it is first hypothesized and then shown that the compressible solver can faithfully reproduce the strictly acoustic Helmholtz tones arising in an impermeable and quiescent channel, namely, with no sidewall injection. Conversely, in the presence of a mean flow, it is speculated and then demonstrated that the compressible biglobal framework can generate a broad spectrum of frequencies and growth rates that not only reproduce the chamber's hydrodynamic modes, but also contain the chamber's longitudinal, transverse, and mixed vorticoacoustic modes. In this process, the solution's sensitivity to the discretization size is carefully examined to ensure convergence. To further demonstrate that the framework can adequately resolve the effects of surface boundaries on the generation of unsteady vorticity, Sec. V is used to compare the biglobal eigensolutions to analytical predictions of vorticoacoustic waves in porous channels including the underlying rotational depths of penetration. The resulting comparison enables us to ascertain the strong dependence of the numerical solution on the penetration number, a non-dimensional parameter that controls the exponential decay rate of the analytical formulation of the vorticoacoustic wave. The sensitivity of the solution to various factors is also explored, including variations in the chamber aspect ratio, mean flow compressibility, viscosity, and wall injection speed. This is followed, in Sec. VI, by a discussion of the wave structure *vis-à-vis* comparisons to both numerical simulations and experimental measurements acquired by other researchers. Having established that the compressible biglobal stability framework is capable of producing reliable predictions of vorticoacoustic and hydrodynamic frequencies and mode shapes, several closing remarks are offered in Sec. VII.

## II. FORMULATION OF THE COMPRESSIBLE BIGLOBAL EQUATIONS

### A. Geometry

The internal flow corresponds to a rectangular channel with porous walls through which gas is injected. The permeable walls allow mass to be introduced in a manner that mimics the injection across the burning surface of a propellant grain in a simulated slab rocket motor.^{60} Since the scope of this analysis is to predict the instabilities that evolve after the gases have been produced,^{61} a constant *V _{w}* is taken to represent the normal speed across the non-regressing wall. A schematic diagram of the planar problem is given in Fig. 1 using both two- and three-dimensional representations. The sketch depicts a chamber that is bounded by two opposing porous walls, an inert wall to the left, and an open section to the right. The simulated gaseous products are injected from the porous walls and accelerated in the streamwise direction. Assuming symmetry about the midsection plane, half of the chamber is investigated using a domain that extends horizontally from $x*=0$ to

*L*and vertically from $y*=0$ to

_{w}*a*.

### B. Normalization

To parameterize the solution, we begin by converting the governing equations to non-dimensional forms. This enables us to identify several non-dimensional parameters that affect the characteristics of the solution. As we proceed, the flow variables are normalized according to

In the above, the asterisk marks dimensional variables and all spatial coordinates are referenced to *a*, the half-height of the chamber. The flow velocities, denoted here by the vector $u\u0303$, are normalized using the speed of sound at the wall, *c _{w}*, whereas $(\rho \u0303,T\u0303)$ are used to denote the density and temperature normalized by their reference values at the wall, $(\rho w,Tw)$. The conservation equations are provided for a non-reacting fluid, with no internal heat generation or radiation. Using

*μ*,

_{w}*k*,

_{w}*C*, and

_{p}*γ*to denote the dynamic viscosity, thermal conductivity, constant pressure specific heat, and ratio of specific heats, one obtains $Rea=\rho wcwa/\mu w$ for the acoustic Reynolds number, and $Pr=cp\mu w/kw$ for the Prandtl number. The normalized system of equations becomes

^{43}

where $\Phi $ stands for the mechanical dissipation function.

### C. Linearization

A linearized system can be obtained by recognizing the presence of two unique scales: the first comprising mean flow transients in *t*, and the second resulting from the oscillatory fluctuations that evolve over a faster timescale, *ωt*. This two-scale structure allows the decomposition of the velocity, pressure, density, and temperature, thus leading to a linearized set of equations. All instantaneous quantities can be separated in terms of a mean flow variable, $G\xaf$, and its unsteady component $g1$,

As alluded to earlier, fluctuations can be further decomposed into three unsteady fields. The first is the compressible, irrotational acoustic wave. Being strictly irrotational, this wave does not satisfy the no-slip condition at the wall, unless the solution is augmented using a vortical wave. Acting as a viscous correction to the compressible field, the vortical wave is generally coupled with the acoustic wave to form the so-called vorticoacoustic wave.^{33} As for the combined, composite wave, it manifests itself at the natural frequencies that are prescribed by the resonant states of the system in the longitudinal or transverse directions.^{57,58}

In this work, we hypothesize that an additional contributor to the instantaneous oscillations can be hydrodynamic in nature. Unlike the vorticoacoustic waves, hydrodynamic fluctuations occur over a spectrum of frequencies that are connected to mean flow instability. The resulting fluctuations may be perceived as being periodic, non-resonant components of the decomposed field.

#### 1. Unsteady variable decomposition

In the present framework, the equations of motion account for both compressibility and rotationality; as such, they have the potential to capture both hydrodynamic and vorticoacoustic waves simultaneously. From this perspective, the requirement to decompose the unsteady variables into their acoustic, vortical, and entropic constituents is circumvented. Using the standard Cartesian nomenclature for the velocity field, the separation of the mean flow and its unsteady components takes the form

At this point, the *O*(1) mean flow functions, $(U,V,W)\u2261(U\xaf,V\xaf,W\xaf)/Mw$, are specifically normalized using the wall injection speed, *V _{w}*; this causes the wall Mach number,

*M*, to appear as a secondary perturbation parameter, while still producing a total mean flow velocity that is normalized by the speed of sound,

_{w}*c*. After inserting the instantaneous variables back into Eqs. (2)–(5) and collecting terms that operate on similar scales, two sets of equations emerge. The first represents the nonlinear steady system describing the compressible mean flowfield. The second corresponds to the first-order collection of unsteady equations. As usual, higher-order terms can be ignored in a linear stability framework. After some algebra, the interaction equations for the fluctuating variables can be extracted and written as

_{w}and

#### 2. Stability criteria

According to the theory of dynamical systems, a flow is considered to be stable if all perturbations about the mean flow decay over time.^{62} From this perspective, stability becomes intertwined with the temporal evolution of perturbations. Mathematically, the flow will be stable if

This criterion will be used at the basis of the upcoming stability analysis.

### D. Mean flowfield

In selecting a suitable mean flow, our strategy is to adopt a known solution that satisfies the steady flow equations. For the rectangular channel configuration, a two-dimensional profile is available.^{27} Being closed-form, its flow variables can be readily incorporated into the present framework and later encoded into any numerical solver with essentially no computational overhead. Moreover, the solution consists of a compressible, fully two-dimensional profile, which satisfies the fundamental assumptions made in this study: the compressibility needed to capture the acoustic propagation and the two-dimensionality needed to identify both longitudinal and transverse fluctuations.

In their work, Maicke and Majdalani^{27} use perturbation tools to derive a closed-form approximation that describes the mean flowfield up to the point where fully choked conditions are reached. Their variables are expanded using a regular Rayleigh-Janzen expansion with respect to $Mw2$. Their solution, which can be readily incorporated into the present framework, is summarized below for the unitary streamfunction and velocity components,

As for the thermodynamic properties, we have

and

### E. Compressible biglobal equations

The biglobal approach is based on a two-dimensional spatial ansatz of the form,

where $g1$ represents any first-order unsteady disturbance in the velocity, vorticity, or thermodynamic variables; these include $u1=u(x,y)ei(nz\u2212\omega t),\u2009v1=v(x,y)ei(nz\u2212\omega t),\u2009\Omega 1=\Omega (x,y)ei(nz\u2212\omega t),\u2009p1=p(x,y)ei(nz\u2212\omega t),\u2009\rho 1=\rho (x,y)ei(nz\u2212\omega t),\u2009T1=T(x,y)ei(nz\u2212\omega t)$, etc. Moreover, the generic function $g(x,y)$ refers to the spatial eigensolution in the two primary dimensions, namely, the longitudinal and transverse directions; *n* denotes the mode number in the third spanwise dimension; and $\omega =\omega r+i\omega i$ stands for the complex frequency. The latter encapsulates both the circular frequency, $\omega r$, and the wave's temporal growth rate, $\omega i$. In principle, the wave amplitudes will decay with the passage of time so long as $\omega i<0$. As such, $\omega i$ can be used at the basis of the temporal stability criterion described briefly in Sec. II C 2. The corresponding modal form, given by Eq. (19), no longer restricts the spatial mode shape to be strictly sinusoidal in any one of the two principal directions *x* or *y*. However, Eq. (19) still assumes a periodic motion in the secondary *z*-direction with a mode number of *n*. With this ansatz in hand, the biglobal form for each unsteady variable can be inserted into the linearized form of the Navier–Stokes, energy, and state equations, thus leading to the following formidable set:

Concerning boundary conditions, one can suppress all three velocity components at the physical boundaries and set a non-reflective extrapolation condition at the outflow section.^{14} This is done to facilitate the smooth convection of vortices at the downstream section.^{63–65} As for the pressure disturbances, one may follow Poinsot and Lelef ^{66} by specifying a reflective condition on *p* that promotes acoustic wave reflections. This is done by imposing $n\xb7\u2207p=0$ on all boundaries. Temperatures and densities can be similarly specified by taking $n\xb7\u2207T=n\xb7\u2207\rho =0$ at all boundaries.

## III. NUMERICAL APPROACH

This section provides a brief summary of the numerical approach and eigensolver setup.

### A. Choice of a spectral method

Wave propagation problems require both high spatial and temporal resolutions that are particularly important for computing hydrodynamic modes that exhibit small scale vortex structures. A proper solution will therefore require a high order, numerically stable approach to properly capture the oscillatory nature of both vorticoacoustic and hydrodynamic waves. In this vein, a spectral collocation method is relied upon to provide a distinct advantage over other discretization techniques. Since this work is focused on analyzing a simple planar geometry, the use of a spectral method can be fairly effective even when taken over relatively coarse grids.

Our spectral approach is based on a global discretization of the solution that is capable of achieving a high order approximation. The resolution of this approximation is made realizable through either of two conventional approaches: Galerkin and collocation methods. On the one hand, the residual of the approximating function in the Galerkin approach must be orthogonal to the basis function. For the present problem with six unknowns, inner products can become laborious to the extent of diminishing the possibility of a successful outcome. On the other hand, the residual in a collocation method vanishes at a given set of points, which happen to be the zeroes of the interpolating polynomials; this renders the solution exact at these specific locations. Given the finite geometry of the porous channel flow domain, the choice of basis functions can be narrowed down to either Chebyshev or Legendre polynomials. Both types of polynomials yield identical rates of convergence and are expressible over the $[\u22121,1]$ interval. However, the maximum pointwise error of a Legendre series exceeds that of the Chebyshev series of the same order.^{67} As such, it has become customary in the combustion stability community, including the present work, to adopt Chebyshev polynomials as the basis functions for the spectral decomposition effort.^{15}

### B. Chebyshev polynomials

With its non-uniformly spaced collocation points, Chebyshev polynomials hold unique advantages when compared to equally spaced grid points, such as Fourier spectral methods, which can be susceptible to the Gibbs phenomenon, or Lagrange interpolation polynomials, which can experience the Runge phenomenon.^{68} The Chebyshev polynomials satisfy the differential equation

whose solution can be represented as

As in Fourier analysis, one can expand a given function in terms of Chebyshev polynomials in lieu of trigonometric functions. This expansion enables us to approximate the function to a desired degree *N*. In this work, the discrete polynomial approximation of a function *f* (*x*) at the collocation points $\xi i$ is expressed as

Note that the weight function $\lambda i(\xi )$ for Chebyshev interpolation takes the form^{69}

where

The foregoing discretization scheme results in a system of *N* equations with *N* unknowns, which needs to be solved simultaneously in order to retrieve the values of *f* at the collocation points.

### C. Pseudo-spectral derivatives

Compared to finite difference methods, collocation techniques are generally attractive because of the enhanced accuracy of their differentiation schemes. Since the residual of the approximating function vanishes at the collocation points, the derivative is rendered locally exact at all Chebyshev points. Within this construct, the differentiation of a polynomial approximation given by Eq. (28) can be written as

As seen from Eq. (31), the calculation of derivatives requires the intermediate evaluation of $\lambda \u2032i$, the derivative of the weight function. The latter can be calculated beforehand and stored in the pseudo-spectral differentiation matrix *D*. Then, with the differentiation matrix in hand, the derivative of a function can be readily expressed as

In this manner, given a finite number of collocation points *N*, the differentiation matrix can be fully resolved. Higher-order derivatives can be determined straightforwardly by raising the matrix to the corresponding power.

### D. Eigenvalue problem

The partial differential equations entailed in the present formulation give rise to a generalized eigenvalue problem of the form

where *ω* and $fi$ define the eigenvalue and solution eigenvector, while $Aij$ and $Bij$ denote the operator matrices. The term “generalized” stems from the matrix $Bij$, which returns the identity matrix in an eigenvalue problem containing a single matrix only.

It may be helpful to recall that, being formulated in a two-dimensional context, additional treatment is needed in resolving partial differential equations. Among the first steps taken in rearranging bounded problems, a variable transformation is applied to map the solution over the interval $[\u22121,1]$ in each direction. In a problem with two spatial variables, *x* and *y*, the general transformations can be taken as

where *A*, *B*, *C*, and *D* mark the boundaries of the domain. The corresponding derivatives become

At this juncture, the domain can be meshed using a two-dimensional grid on directionally independent Chebyshev points. To compute the operator matrix, the technique consists of building a single operator matrix from a two-dimensional collocation. This is made possible through the Kronecker product. More specifically, the Kronecker product of two matrices, $Aij$ and $Bij$, generates a block matrix comprising $aijBij$ terms that can be aligned according to

In fact, the Kronecker product can be used to define the pseudo-spatial differentiation matrices. The derivatives with respect to *x* take the form $(DN)\u2297(IN)$ while derivatives with respect to *y* take the form $(IN)\u2297(DN)$. To calculate higher-order derivatives, the differentiation matrix $(DN)$ is simply raised to the corresponding power.

In this work, each operator matrix is conveniently subdivided into sets that reproduce the coefficients of five dependent variables in five equations. Note that both the pressure variable and the state equation can be eliminated by direct substitution. Once the density and temperature are found, the pressure can be straightforwardly deduced from the state equation. Moreover, each subset encompasses the corresponding variable multipliers evaluated at each grid point. For example, given *N *×* N* collocation points, the operator matrices $Aij$ and $Bij$ will consist of $25N4$ elements. In the interest of clarity, a graphical representation of the operator matrices is provided below:

## IV. VORTICOACOUSTIC WAVE MODELING

In this section, the results from the present formulation are compared and validated against well-established solutions that happen to be available both in the longitudinal and transverse directions of propagation. The prevailing direction of propagation depends on the type of problem, geometrical dimensions, and other physical parameters. The ability to solve for both the longitudinal and transverse modes simultaneously enables us to compare their modal growth rates for the same enclosure. It also facilitates performing parametric trade studies of simulated slab motors where the influence of several key parameters on the oscillatory motion can be evaluated.

### A. Acoustic wave representation

The decomposition of the oscillatory component of the flowfield has been implemented in stability analyses based on the method used by Chu and Kovàsznay.^{29} Introduced briefly in Sec. I, this decomposition reduces the governing equations into three simplified sets. The first set describes the production, propagation, and absorption of pressure waves, whose compression and expansion lead to the onset of acoustic waves. The second set captures the production, convection, and dissipation of vorticity fluctuations. Finally, the entropy wave stems from the production, convection, and diffusion of heat. This decomposition offers valuable insight into the physical nature of the problem.^{1}

Understanding the coupling between acoustic and vortical wave motions has been shown to be possible and contingent upon a proper representation of the flow near the wall, specifically, by enforcing the no-slip requirement.^{15,19} At the outset, the surface-driven vortical waves become coupled to the strictly acoustic motion, independently of the additional rotational contributions associated with hydrodynamic flow breakdown. Being independent of resonant modes, these vortical waves, which remain unaccounted for in the asymptotically derived vorticoacoustic formulations,^{40–42} will be referred to as hydrodynamic waves. In addition to its ability to capture the hydrodynamic wave structures, the present framework obviates the need to decompose the oscillatory component of the flowfield at the forefront of the analysis. As a result, it can reliably capture the vorticoacoustic and hydrodynamic waves while taking into account their possible interactions.

Although it may be speculated that the compressible formulation captures acoustic waves, the framework can be tested by exploring a basic case. Let us consider the one-dimensional acoustic wave equation,

The analytical solution to Eq. (37) for a simple channel can be compared to the present eigensolution under the same quiescent conditions. The purpose of this simplified case is to verify that the framework can reliably predict the acoustic waves and their theoretical Helmholtz frequencies. In practice, these frequencies and mode shapes are influenced by the mean flow, which is fully accounted for in the present formulation. However, because well-known benchmark cases that describe the resonant acoustic frequencies with mean flow effects are not widely available, it is helpful to first establish the effectiveness of the present scheme at reproducing the expected modes. To concentrate on pure acoustic waves, the solver is initialized for a channel with a length to half-height ratio of $L=Lw/a=5$ in a flow with no wall injection.

Forthwith, our computations are seen to capture the full spectrum of frequencies appearing in the chamber over the entire range of spatial and temporal scales associated with turbulence generation. Numerically, the resulting frequencies are found to coincide with the theoretically predicted resonant frequencies for a closed rectangular chamber where the Hertzian frequency is given by

Here, *l* and *m* appear as positive integers representing the wave numbers in the longitudinal and transverse directions, respectively. The ensuing wave structures follow the expected sinusoidal forms, as shown in Fig. 2 for the pressure; the latter has been conveniently normalized by its maximum amplitude.

### B. Eigenvalues in porous channel and slab motor flow simulations

Using the compressible Taylor flow given by Eq. (15), the solver can be used to extract a spectrum of eigenvalues that contain both real and imaginary components for a specific set of parameters. As stated earlier, the real part, $\omega r$, represents the frequency of oscillations for each eigenmode, and the imaginary part, $\omega i$, determines the stability growth rate. As usual, unstable modes are identified by a positive *ω _{i}*. In the present simulations of the porous channel flow configuration, the tangential mode number is taken to be

*n*=

*0, thus imposing no changes in the*

*z*direction. In what follows, Fig. 3 is used to illustrate the spectrum of frequencies that are obtained for a simulation with $Mw=0.05$ and $Rea=2000$ in a chamber with an aspect ratio of 5. The graph reports a collection of acoustic frequencies with the first few longitudinal modes labeled as L1-L4. Also shown are the first few transverse modes, labeled as T1–T4, with trailing lines stemming from each. The latter represent the corresponding mixed modes. The modes appearing at the bottom left corner of the graph can be identified as being the hydrodynamic modes; these are known to propagate at much lower frequencies. In this basic case, the hydrodynamic modes remain temporally damped with deceasing growth rates. Section VI discusses these particular modes in more detail.

Since the vertical axis represents the temporal growth rate, the graph can be split into stable and unstable regions. In this case, only two eigenmodes undergo temporal amplification and are thus considered unstable. These are the first and second transverse modes, T1–T2, which correspond to wave propagation that is perpendicular to the porous wall. Compared to the theoretical acoustic frequencies calculated from Eq. (38), we find the frequencies to be slightly shifted to the left. More precisely, the frequencies predicted by the numerical solver are found to be consistently lower than their Helmholtz values obtained in a quiescent medium. The deviations from the Helmholtz values can be attributed to the effect of the mean flow, which only increases at higher speeds and in longer chambers. In practice, this slight shift in the frequencies is also observed in various experiments, including those reported by Casalis and Vuillot,^{71} and in other geometric configurations.^{28}

### C. Solution sensitivity to *N*

Given that the precision of the numerical solver is controlled to some extent by the number of collocation points, it is helpful to examine the sensitivity of the solution to the parametric value *N*. When using Chebyshev polynomials, at least *π* points per wavelength are needed on average to ensure convergence.^{68} For example, to obtain converged results for up to the tenth oscillation mode in each direction, *N* needs to be larger than 5 wavelengths and so we take $N>5\pi \u224815.71$. Table I displays the relative error in the numerical results compared to classic theory. For consistency, our computations are performed for a basic case with *L *=* *5 and no mean flow. Note that Table I catalogs the longitudinal acoustic frequencies only. This comparison confirms that the results for *N *=* *15 are accurate up to mode 7. Starting at mode 8, the error starts to grow beyond 3.7%. Although the error grows non-monotonically, it remains rather small. Now when the collocation points are increased to *N *=* *20, the results become accurate up to mode 11. To be conservative, the remaining simulations are performed with *N *=* *60, a value that leads to a well-resolved two-dimensional mesh consisting of 3600 grid points. We find that increasing *N* further does not affect the results in the range of frequencies of interest, namely, up to the twentieth mode. As for the unsteady, Stokes-like boundary layer at the wall, decreasing the wall Mach number at a fixed value of *Re _{a}* (or viscosity) leads to a reduction in the penetration depth. This in turn would require a larger value of

*N*to resolve the rapid changes that occur in the vicinity of the wall. In this vein, a very large value of

*N*can become computationally prohibitive, and the alternative option can be to decrease

*M*while increasing

_{w}*Re*or

_{a}*L*. In Sec. V C, we show that the depth of penetration of unsteady vorticity remains essentially unchanged when the penetration number, $Mj=ReaMw3cw2/(\omega r*2a2)=ReaMw3cw2/(4\pi 2a2f2)$, is held constant.

^{43}

. | . | N = 15
. | N = 20
. | ||
---|---|---|---|---|---|

Mode . | Theoretical . | Numerical . | Error . | Numerical . | Error . |

1 | 0.6283 | 0.6283 | 0.0000% | 0.6283 | 0.0000% |

2 | 1.2566 | 1.2566 | 0.0000% | 1.2566 | 0.0000% |

3 | 1.8850 | 1.8850 | 0.0000% | 1.8850 | 0.0000% |

4 | 2.5133 | 2.5132 | 0.0040% | 2.5133 | 0.0000% |

5 | 3.1416 | 3.1416 | 0.0000% | 3.1416 | 0.0000% |

6 | 3.7699 | 3.7670 | 0.0769% | 3.7699 | 0.0000% |

7 | 4.3982 | 4.3990 | 0.0182% | 4.3982 | 0.0000% |

8 | 5.0265 | 4.8390 | 3.7302% | 5.0275 | 0.0199% |

9 | 5.6549 | 5.3830 | 4.8082% | 5.6459 | 0.1592% |

10 | 6.2832 | 6.2832 | 0.0000% | 6.2832 | 0.0000% |

. | . | N = 15
. | N = 20
. | ||
---|---|---|---|---|---|

Mode . | Theoretical . | Numerical . | Error . | Numerical . | Error . |

1 | 0.6283 | 0.6283 | 0.0000% | 0.6283 | 0.0000% |

2 | 1.2566 | 1.2566 | 0.0000% | 1.2566 | 0.0000% |

3 | 1.8850 | 1.8850 | 0.0000% | 1.8850 | 0.0000% |

4 | 2.5133 | 2.5132 | 0.0040% | 2.5133 | 0.0000% |

5 | 3.1416 | 3.1416 | 0.0000% | 3.1416 | 0.0000% |

6 | 3.7699 | 3.7670 | 0.0769% | 3.7699 | 0.0000% |

7 | 4.3982 | 4.3990 | 0.0182% | 4.3982 | 0.0000% |

8 | 5.0265 | 4.8390 | 3.7302% | 5.0275 | 0.0199% |

9 | 5.6549 | 5.3830 | 4.8082% | 5.6459 | 0.1592% |

10 | 6.2832 | 6.2832 | 0.0000% | 6.2832 | 0.0000% |

Finally, in what concerns the choice of the characteristic Reynolds number, our *Re _{a}* differs from the wall injection Reynolds number, $Rew=\rho wVwa/\mu w$, which is widely used in the hydrodynamic stability community.

^{14,15}The two parameters are related through the wall Mach number, namely, $Rew=MwRea$. For example, using (

*M*,

_{w}*Re*) combinations of (0.05, 2000) or (0.025, 20 000) will lead to $Rew$ values of 100 and 500, respectively. Similar ranges are used by Boyer

_{a}*et al.*

^{18,19}Recalling that the present simulation is focused on a cold-flow setup, an even smaller value of $Rew$ is obtained when reducing the wall Mach number further to 0.003 at a fixed value of $Rea=2000$. Hence, for $Mw=0.003$, which is characteristic of solid rocket motors, simulating a small wall Reynolds number of 100 would require the use of $Rea=33\u2009333$. The latter would demand a large value of

*N*that could possibly exceed the user's hardware capabilities. Practically, however, this situation should not be concerning because wall Reynolds numbers that are characteristic of solid rocket motors often exceed $Rew=2000$. In what follows, the effects of different parameters on the retrieved eigenmodes and eigensolutions will be systematically examined and compared “whenever possible,” to predictions from other available models and experiments.

## V. PARAMETRIC STUDIES AND COMPARISONS

### A. Chamber aspect ratio and mean flow effects on frequency

We begin by exploring the effect of varying the aspect ratio on the spectral frequencies. We note, at the outset, that by increasing the chamber length-to-half-height *L*, the acoustic frequency is decreased while the streamwise velocity is correspondingly increased; as such, the effect of the mean flow on the frequency is enhanced with successive increases in chamber length for a fixed injection Mach number. To isolate the effect of the mean flow itself, the viscosity, and hence the acoustic Reynolds number, as well as the injection Mach number at the wall, are held constant at 2000 and 0.007, respectively.

Figure 4 displays the dimensional frequency spectra for three aspect ratios of *L *=* *10, 25, and 50. Here, solid circles are used to denote modes of hydrodynamic nature, while solid triangles are used to identify the vorticoacoustic modes that resonate with the acoustic modes. The purely acoustic or Helmholtz frequencies in an identical chamber with impermeable walls are superimposed using vertical dashed lines. A cursory comparison of the three cases in question confirms the strong dependence of the vorticoacoustic frequencies on the mean flow. Going from Figs. 4(a) to 4(c), the hydrodynamic frequencies may be seen to remain essentially unchanged. However, the vorticoacoustic frequencies shift to lower values in longer chambers where the mean flow velocity is permitted to grow proportionately larger. This result is not surprising: it is well established from a simple mass balance that the mean flow will increase linearly with the length of the chamber due to mass being added uniformly along the porous walls.^{72} In contrast, the effect of increasing *L* on the hydrodynamic frequencies appears to be weak. This may be expected by virtue of the primary mode of initiation of hydrodynamic waves being from shear layer instabilities, and not from coupling with the fast traveling acoustic waves; only the latter prove to be strongly dependent on the chamber's geometric characteristics. More specifically, as we alternate from Figs. 4(a) to 4(c), the hydrodynamic modes do not seem to shift, unlike their vorticoacoustic counterparts, which appear at the chamber's resonant frequencies and will, therefore, change with the geometry.

Interestingly, the present simulations and Fig. 4 help to confirm the suitability of the outflow boundary conditions imposed on the time-dependent velocity and, by extension, vorticity fields. In theory, unsteady vortices can convect through the open section without affecting the motion upstream. This condition is satisfied because the modal frequencies remain independent of the location assigned to the outflow boundary, be it *x *=* *10, 25, or 50. As for the stability companion, the growth rate *ω _{i}* continues to move (with

*L*) closer to the unstable region, thus leading to lower stability margins in longer chambers. This behavior is consistent with the mean flow character because longer chambers allow more mass to be introduced and accelerated in the axial direction. In summary, the compressible eigensolver predicts a monotonic decline in stability growth rates and a small, albeit persistent decrease in vorticoacoustic frequencies as

*L*is increased.

### B. Mean flow effect on the acoustic pressure waveform

Before undertaking comparisons to existing models that do not take into account the effects of the mean flow on the acoustic mode shapes, it may be helpful to examine the behavior of the latter in the presence of wall injection. Given that the unsteady vorticity waves are driven at the boundaries by the sweeping acoustic motion, any spatial distortion of the acoustic mode shapes will have a direct bearing on the intimately coupled vorticoacoustic waveform. By turning the wall injection speed on and off, it is possible to verify through the compressible eigensolver that the presence of a non-negligible mean flow will indeed shift the acoustic pressure oscillations away from their well-established sinusoidal shapes.

One benefit of solving the complete set of biglobal equations for the coupled disturbances is the ability to predict the actual pressure waveform associated with a given vorticoacoustic mode without imposing any restrictions on the pressure distribution. In fact, analytical formulations of vorticoacoustic waves have not yet managed to incorporate the skewness effect that the mean flow can impart on the standing pressure or velocity mode shapes.^{40–42} In existing asymptotic approximations, a purely sinusoidal form of the acoustic pressure is assumed. In practice, this proves to be a reasonable first approximation only, as we illustrate below using a typical test case.

For a case of *L *=* *50, $Mw=0.007$, and $Rea=2000$, we readily extract the unsteady pressure eigensolution for the first vorticoacoustic mode. We then display its temporal evolution in Fig. 5 at eight equally spaced timelines corresponding to $t=0,18\tau ,28\tau ,\u2026,78\tau $, where *τ* stands for the cyclical period. For the sake of comparison, the standing pressure waveform based on classic acoustic theory is superimposed in Fig. 5 using a dashed curve. When the two models are compared, the motion of the compressible eigensolution may be seen to slightly deviate from the strictly sinusoidal motion of a standing wave. The time-dependent sloshing behavior that ensues does not appear when wall injection is suppressed. As such, the distortion can be strictly attributed to the convective mean flow effect on the acoustic mode shape. So far, these flow-induced deviations in the mode shapes have eluded analytical approximations as they could only be predicted by solving the compressible equations of motion computationally.

### C. Comparisons to analytical vorticoacoustic waveforms

Another avenue for verification can be pursued by comparing the velocity eigensolutions predicted from the compressible biglobal stability formulation to a closed-form asymptotic approximation for the vorticoacoustic wave in a porous channel undergoing internal pressure oscillations.^{40–42} In the asymptotic work, the unsteady field is subdivided into two distinct fluctuations in a manner that is consistent with the Helmholtz–Hodge decomposition process:^{1} an inviscid, irrotational, compressible disturbance and a viscous, rotational, incompressible correction that is generated at the boundaries by the sweeping pressure oscillation. Then using perturbations in the normalized pressure amplitude and wall Mach number, the time-dependent rotational wave is determined based on the coupling at the boundaries with a simple standing acoustic waveform for the compressible disturbance evaluated at the chamber's Helmholtz frequencies. The resulting analytical solution for the oscillatory axial velocity, given by Eq. (3.38) in Majdalani,^{42} employs a Wentzel–Kramers–Brillouin (WKB) expansion to obtain

where

Here, *S* denotes the Strouhal number and *F*(*y*) represents the mean flow characteristic function; as shown by Berman,^{73} a self-preserving motion in porous channels with constant cross sections can be expressed as

For a direct comparison, a simple incompressible mean flow profile can be chosen, such as $F(y)=y$ or $F(y)=sin\u2009(12\pi y)$, which lead to the porous channel's irrotational and rotational mean flow profiles, respectively. These are given by^{43}

Interestingly, the irrotational motion corresponds to the planar flow analog of the incompressible Hart–McClure profile,^{31,32} which preceded the use of Taylor's rotational profile^{11} in rocket stability calculations,^{6,9,10} including those incorporating the effects of particle entrainment,^{12} grain taper,^{74} and headwall injection.^{75}

It may be instructive to also note that, based on Eq. (39), a characteristic parameter, $Mj=ReaMw3cw2/(\omega r*2a2)$, is found to control the penetration depth of the rotational wave oscillations measured from the sidewall. In essence, $Mj$ controls the exponential decay rate of the wave amplitude as $y\u21920$. Following the theory of time-dependent laminar flows,^{76} the penetration depth is defined to be the distance from the wall to the point where 99% of the rotational wave amplitude would have vanished.^{33} When $Mj$ is fixed, the penetration depth remains essentially unchanged.^{43} For example, by varying both the Reynolds number and the wall Mach number in such a way to preserve $Mj$, a similar penetration depth of the axial velocity fluctuations is obtained for the three cases illustrated in Fig. 6. In fact, the identification of the penetration depth control parameter $Mj$ may be viewed as a compelling example of the valuable insight that can be gained from the complementary asymptotic expansion of this problem.^{43} However, the latter remains restricted to simple geometries and incompressible mean flowfields.

To make further headway, the vorticoacoustic wave is computed and compared in Fig. 7 using two cases that explore the effects of varying the Reynolds and wall injection Mach numbers. We recall that these two parameters control the injection velocity and viscous effects, which directly impact the unsteady boundary-layer thickness, as discussed in Sec. V C. Here, the computed and theoretical approximations of the oscillatory axial velocity are shown in Fig. 7(a) for $Rea=2000$ and $Mw=0.05$, and then in Fig. 7(b) for $Rea=10\u2009000$ and $Mw=0.04$. All calculations are carried out midway through the chamber at *x *=* *2.5, namely, at the acoustic velocity node, assuming pressure oscillations at the first harmonic frequency. Forthwith, a qualitative assessment of the plots reflects a fair agreement, particularly in the similarity of the wave shape and spatial periodicity. However, distinctive differences can be seen in both the amplitude and the penetration depth of the wave. These differences can be directly attributed to the effect of the mean flow, which causes the frequency to shift below its Helmholtz value, and the mode shape to deviate from its standing sinusoidal form. Both of these variations occur because the frequency of oscillation happens to be a primary controlling factor in determining the penetration number and, as such, the penetration depth. Since $Mj$ stays inversely proportional to $\omega r*2$, higher frequency oscillations are anticipated to promote a shallower propagation of the vorticoacoustic boundary layer. Moreover, given that successive increases in the mean flow velocity will serve to lower the frequency, the vorticoacoustic boundary layer will be expected to slightly increase in size and extend deeper into the chamber, as predicted by the compressible biglobal solver, when a more appreciable mean flow is prescribed. Recalling that the asymptotic model in Eq. (39) does not capture the frequency shift caused by mean flow variations, it may be safely argued that these effects are more accurately represented by the compressible eigensolver and the solid lines in Fig. 7.

### D. Wave structure in compressible and incompressible mean flowfields

In former biglobal stability investigations of the flowfield in simulated rocket chambers,^{14,15,19} incompressible models have been routinely used to represent the mean flow profile. Since the present eigensolver is fully compressible, it would be helpful to examine its sensitivity to small deviations in the mean flow description. Specifically, it would be instructive to determine whether the use of an incompressible mean flow as the basic motion is sufficient, or whether adopting a compressible model is warranted. This task can be accomplished by comparing the spectra corresponding to an incompressible versus a compressible mean flow profile. In the case of the former, a porous channel solution may be traced back to Taylor.^{11} Despite its simplicity, Taylor's profile is implemented in the context of solid rocket motor instability first by Varapaev and Yagodkin^{3} and, later, by Beddini,^{5} Couton *et al.*,^{77} Avalon *et al.*,^{78} Kourta,^{79} Griffond and Casalis,^{9,10} Féraille and Casalis,^{12} Venugopal *et al.*,^{80} and many others. It can be reproduced from the normalized mean flow streamfunction,

Alternatively, a compressible form of the Taylor profile can be adopted. Based on a Rayleigh-Janzen expansion in the wall injection Mach number (squared), one actually gets^{27}

In this work, using these two profiles as comparators, simulations are performed assuming $Mw=0.05$ and $Rea=2000$, with a length to half-height ratio of 5. Although typical wall Mach numbers for solid rocket motors revolve around $O(10\u22123)$, a value of $Mw=0.05$ is deliberately chosen to enhance the effects of compressibility. Nonetheless, based on the modal spectra reported in Fig. 8, only small differences in the eigenfrequencies and stability growth rates may be observed; particularly, a minor reduction in both *ω _{r}* and

*ω*is realized when density variations are permitted. These differences suggest that the retention of dilatational properties in the mean flow will have a slight stabilizing effect. It will also lead to a slight decrease in frequency, thus further justifying the deviations observed in Fig. 7 between the asymptotic and compressible biglobal predictions. In practice, however, these deviations can be ignored when $Mw=O(10\u22123)$ or smaller.

_{i}Using the same basic properties and midway through the chamber, Figs. 9(a)–9(d) compare typical axial velocity and vorticity distributions stemming from both profiles for the first vorticoacoustic mode. The nearly indiscernible plots help to ascertain that compressibility in the mean flowfield has a negligible effect on the unsteady velocity and vorticity fields. We recall that Maicke and Majdalani^{27} show that compressibility effects become more noticeable in the downstream sections of long porous ducts. Practically, they become only discernible beyond 40% of the critical length *L _{s}* of a porous channel, where

*L*is sufficiently large to induce sonic conditions in the absence of a converging-diverging nozzle. For this reason, we also display in Figs. 9(b)–9(d) the axial and vorticity profiles at the aft-end of the domain. Despite the increased disparity between the two curves in the downstream direction, the differences between them remain relatively small. Then, given the short length

_{s}*L*in these particular cases in relation to the critical (sonic) length

*L*, the need to account for compressibility in the mean flow seems to be secondary, at least insofar as $Mw\u226a1$ and entropy effects are not considered.

_{s}In this vein, and since no compressible flow study is complete without due consideration of heating effects, the unsteady entropy in this problem can be estimated from a perturbation of the thermodynamic equations. In their dimensional form, these can be written as

Letting *s*^{*} = $s\u0303*$/*c _{p}*, Eq. (45) can be readily normalized and expanded using Eqs. (1), (7) and (19). The fluctuating component of the entropy can be subsequently retrieved, namely,

By way of illustration, *s*-profiles are evaluated for both the compressible and incompressible mean flowfields and shown in Fig. 10 at two axial locations. Unlike the behavior of the velocity and vorticity profiles, these two cases exhibit substantial differences. This behavior can be attributed to the transfer of energy between the base flow and the unsteady disturbances. We recall that the energy used in the compression and expansion processes in the steady field manifests itself through fluctuations in the density and temperature. This characterization is consistent with the Helmholtz–Hodge decomposition that enables us to split the unsteady disturbances into an irrotational (acoustic) fluctuation and a solenoidal (vortical) counterpart. As shown by Chu and Kovàsznay,^{29} entropy variations will accompany pressure fluctuations as a result of energy dissipation (or absorption) at the walls. Moreover, the diffusion of hot spots, which is caused by entropy fluctuations, will only generate a weak velocity disturbance, which is necessary to satisfy the conservation of mass when density is varied. For the case of an incompressible base flow, entropy fluctuations become more concentrated halfway between the midsection plane and the wall, and vanish categorically along the midsection plane. They may be ascribed to the flow turning effect,^{81} which is locally intensified and captured through the density variations that are still permitted in the unsteady component of the present formulation.

Because entropy fluctuations are significantly influenced by compressibility, contour plots for both cases are provided in Fig. 11. These graphs help to identify the areas where entropy production is most concentrated for the case of an incompressible fluid. Moreover, for a fully compressible motion, the gradual increase in dilatational effects as the flow advances downstream is visually demonstrated, thus confirming the results of Maicke and Majdalani.^{27} Because of the appreciable differences, particularly in the density and temperature fluctuations caused by dilatational effects in the mean flow profile, only compressible base flows will be considered hereafter.

### E. Effect of viscosity

The effect of viscosity on the unsteady variables can be explored by computing the eigensolutions at three characteristic Reynolds numbers while maintaining the same injection speed. Our chosen parameters consist of $Mw=0.03$ and three acoustic Reynolds numbers of *Re _{a}* = 200, 1000, and 2000. We recall that

*Re*is based on the speed of sound at the wall irrespective of the injection speed. As such, changing the Reynolds number and maintaining the same thermodynamic conditions at the wall implies a change in viscosity.

_{a}Forthwith, Fig. 12 displays the frequency spectra for the three cases at hand. The spectral data helps to identify two sets of modes rather distinctly, namely, those of acoustic and those of hydrodynamic origins. With respect to the acoustic modes, their oscillation frequencies remain virtually unaltered for all three cases. However, their growth rates can be seen to increase with successive increases in the acoustic Reynolds number, which is inversely proportional to the kinematic viscosity. From a physical perspective, this behavior suggests that the system will gradually become less stable as viscosity is decreased. Unsurprisingly, viscous dissipation may be viewed as a source of damping that has a rather stabilizing effect on wave propagation. As for the hydrodynamic modes, they remain strongly dependent on the rotationality of the flow, which is further promoted by increasing levels of viscosity. As a result, the hydrodynamic frequencies shift leftward while becoming more stable with successive increases in viscosity.

The damping role of viscosity can also be inferred from Fig. 13, where the axial velocity profiles are showcased in the normal direction for the first vorticoacoustic mode. In undertaking this comparison, all wave amplitudes obtained from our linear eigensolver are normalized by their peak values. Qualitatively, the damping rate can be gauged by the number of spatial undulations, which increase as the acoustic Reynolds number is raised, starting with a single undulation at $Rea=200$ and incrementing to three fluctuations at $Rea=2000$. Also noteworthy is the gradual increase in the penetration depth of vorticoacoustic waves with successive reductions in viscosity, namely, as $Rea$ is increased at a fixed value of $Mj$. The resulting trend is fully consistent with the asymptotic behavior of Eq. (39).

Given the direct connection between velocity gradients and vorticity computations, similar effects may be seen to extend to the iso-vorticity contours depicted in Fig. 14. These are calculated for the first vorticoacoustic mode at $Mw=0.03$ and three decreasing levels of fluid friction corresponding to *Re _{a}* = 200, $1000$, and $2000$. As usual, vorticity magnitudes are derived from the appropriate gradients of both axial and normal components of the fluctuating velocities.

The first characteristic that may be inferred from this comparison relates to the number of vortices that ripple across the domain. For the case of $Rea=200$ in Fig. 14(a), it may be observed that only one main vortex is generated at the wall, which then dissipates rather quickly without penetrating deeply into the chamber. When viscosity is reduced to $Rea=1000$ and $2000$ in Figs. 14(b) and 14(c), multiple vortices are formed and transported more prominently across the chamber. This behavior can be attributed to the convective effect of the mean flow on the unsteady disturbances.

### F. Effect of wall injection

To study the effect of varying the injection speed on the unsteady variables, three different cases are simulated. These correspond to $Mw=0.01$, 0.03, and 0.05 at a Reynolds number of $Rea=2000$. Here too, the eigenmodes for the three cases are extracted and presented in Fig. 15. In order to examine the effects of the wall injection Mach number on flow stability, the temporal growth rates are closely compared. Interestingly, we find that, at the higher injection speeds, longitudinal modes tend to exhibit lower growth rates, in contrast to the transverse modes, which experience increased growth rates. It can be, therefore, argued that a higher *M _{w}* can stabilize the longitudinal modes but destabilize the transverse and hydrodynamic modes. However, as far as the actual oscillation frequencies are concerned, they may be seen to shift slightly leftward to lower values. Nonetheless, the shift observed in Fig. 15 remains minimal and will be examined in more detail in Sec. VI.

To further examine the effect of *M _{w}* on the wave structure, iso-vorticity contour graphs are produced for the first vorticoacoustic mode and then presented in Fig. 16 at three increasing values of the wall Mach number. At the lowest injection speed considered in Fig. 16(a), only one rotating vortex structure can be seen to develop and remain confined to the vicinity of the sidewall. However, as the injection speed is gradually raised in Figs. 16(b) and 16(c), the vortex structures seem to intensify and extend deeper into the chamber. This behavior may be viewed as being consistent with the acoustic boundary-layer transport mechanisms affecting this problem. In this context, the underlying diffusion and convection mechanisms remain strongly controlled by the kinematic viscosity and the injection speed, respectively. Presently, the viscosity, and thus dissipation, remain constant for the three cases in question. The comparison also illustrates how the vorticity-laden boundary layer can be transported by way of mean flow convection. To control the level of transport and depth of penetration of wall-generated vorticity waves, a balance between convection and dissipation is clearly needed. This balance is prescribed analytically by the penetration number,

^{42}which is discussed in more detail in Sec. V C.

## VI. HYDRODYNAMIC WAVE STRUCTURE

In the context of biglobal stability analysis, the term “hydrodynamic” is often used to denote a non-resonant, periodically fluctuating field variable. On the one hand, the corresponding wave type is mainly caused by mean flow breakdown and the generation of vorticity during the mean flow transitioning into a turbulent field. On the other hand, the mechanism of parietal or surface vortex shedding^{19} stands among the chief contributors to the development and triggering of resonant-state oscillations or, more precisely, excited vorticoacoustic waves. In the present framework, the corresponding waves are captured using a compressible, time-dependent viscous model, namely, one that accounts for both hydrodynamic and vorticoacoustic waves with no prior requirement for flow decomposition.

To proceed with this effort, the fully compressible modes will be resolved and then compared to their theoretical and experimental counterparts obtained in other investigations, such as those conducted at the *Office National d'Études et de Recherches Aérospatiales* (ONERA). More specifically, after applying a filter that eliminates spurious modes, our theoretical predictions will be compared to both computational and experimental measurements that are independently obtained by researchers at ONERA. In fact, led by Casalis and co-workers, the French Aerospace Laboratory has been at the forefront of modeling idealized rocket chambers using both one-dimensional^{9,10,12,13} and two-dimensional hydrodynamic stability theories.^{14,15,19,22} By way of verification, our compressible stability findings will be compared to those reported in two pertinent studies by Casalis and co-workers^{17–19} as well as others.^{6,71,78,82} In the former, comparative data will be excerpted from an incompressible solver that can be robustly relied upon to capture the hydrodynamic waves.^{18,19} In the latter, experimental measurements acquired from a porous channel flow apparatus named, “*Veine d'Essai de la Couche Limite Acoustique*” (VECLA), will be used to validate our spectral predictions.^{6,71,78,82}

### A. Comparison to numerical stability simulations

We begin by comparing the filtered modes obtained from the compressible biglobal solver to those produced from an incompressible (finite element) code developed by Boyer *et al.*^{18,19} In their work, Boyer *et al.*^{18,19} investigate the implications of minor wall defects at the grain surface of a simulated solid propellant on the generation of hydrodynamic disturbances. The latter are characterized alongside the effects of varying the Reynolds number on the predicted frequencies and growth rates. We recall that, although the solver in question is based on the biglobal stability approach, it differs from the present framework in its reliance on the linearized form of the incompressible Navier–Stokes equations.

For the sake of illustration, eigenvalues are extracted for the case of a chamber aspect ratio of $L=8,$ an acoustic Reynolds number of $Rea=2000$, a wall Mach number of $Mw=0.003$, and an injection speed of *V _{w}* = 1 m/s. The comparison is provided in Fig. 17, where the stability growth rate

*ω*is displayed against the oscillatory frequency

_{i}*ω*. Note that the results from the present solver are designated by solid circles while those from Boyer

_{r}*et al.*

^{18,19}are represented by hollow triangles. Overall, the comparison suggests a favorable agreement between the two models in pinpointing the oscillation frequencies, notwithstanding the small discrepancies that affect their temporal growth rates. Another difference can be detected in the high frequency range where the compressible framework produces a fewer number of modes than the incompressible solver. Upon closer examination, these underlying disparities can be actually attributed to two main factors.

First, it may be speculated that the dilatational energy that is consumed by the compressible waves that are particular to the present framework can lead to fewer modes than projected otherwise by an incompressible framework. Since the compression process absorbs energy, it is likely that only a subset of the incompressible modes will manifest themselves in a compressible medium. Conversely, a highly resolved incompressible framework will have a tendency to overpredict the number of modes by generating spurious modes, especially in the treatment of compressible motion.

Second, it may be argued that the number of modes generated by an eigensolver will remain somewhat dependent on the number of collocation points *N* that is used in the discretization process. To confirm that the number of modes at higher frequencies can be affected by the grid resolution, results can be extracted at gradually increasing numbers of collocation points. In the numerical solver, one finds that increasing the number of collocation points does not affect the frequency spectrum in the practical (low frequency) range of interest. It does, however, continue to provide additional eigenmodes in the high frequency range. Clearly, the improved resolution that can be achieved using an even larger number of collocation points enhances the framework's effectiveness at capturing the smaller spatial and temporal scales that typically accompany higher frequency disturbances. When all of these factors are accounted for to justify the deviations seen in Fig. 17, the foregoing comparison becomes corroborative of the ability of a compressible biglobal framework to reproduce, given a sufficient resolution, the hydrodynamic modes ascribed to a conventional incompressible solver.

Before leaving this topic, it may be instructive to pause and discuss the amplified mode that appears leftward of the theoretical acoustic chamber frequency in Fig. 17. Here, the Helmholtz frequency is identified on the graph using a vertical dashed line whereas the amplified mode is labeled “1L,” being related to the first longitudinal acoustic mode. In practice, this mode emerges as a by-product of the coupling that takes place between shear layer instabilities in the flow and chamber acoustics, i.e., a mechanism that can lead to the triggering of the first vorticoacoustic mode. This mode differs from the theoretical acoustic frequency, although a less amplified hydrodynamic mode that is closer to the Helmholtz value also exists (at the bottom right).

In another study by the same group, Casalis *et al.*^{17} provide two-dimensional contour plots for the streamwise and normal velocities using a similar configuration to ours. In Figs. 18 and 19, snapshots of their data are presented side-by-side with the wave patterns extracted from the compressible eigensolver. Graphically, the results exhibit a striking resemblance in the wave structure and the approximate positions of their peaks and troughs. Bearing in mind that the waves are periodic in time, the foregoing comparison may be viewed as providing a qualitative verification of the waveforms associated with the hydrodynamic modes.

Along similar lines, unsteady velocities are extracted for a representative case with $Rea=2000,\u2009Mw=0.05$, and *L *=* *5. The corresponding velocity vector plots are shown in Fig. 20 for the third and fifth hydrodynamic modes. These plots display distinct patterns of vortex structures, particularly in the downstream section of the channel, that share similar qualitative features to those produced by Casalis and co-workers.

### B. Comparison to experimental measurements

In searching for usable experimental data, two cold-flow experimental apparatuses developed by Avalon *et al.*^{82–84} come to mind: the porous channel flow configuration introduced as VECLA in Sec. VI, and the so-called “*Veine Axisymétrique pour Limiter le Développement des Oscillations*” (VALDO). Although other experimental facilities exist,^{37,38} VECLA and VALDO provide the advantage of being specifically designed to investigate the effects of parietal vortex shedding on pressure oscillations in both planar and axisymmetric chamber configurations. Both facilities simulate the gaseous mass injection along a solid propellant grain by forcing air into a chamber through a uniformly porous wall made of sintered aluminum.

In this section, the experimental measurements reported by Casalis and Vuillot^{71} in the context of a VECLA investigation will be qualitatively compared to simulation data retrieved from the compressible eigensolver. The corresponding VECLA apparatus consists of a rectangular parallelepiped whose height can be modified. Air injection is maintained through the lower perforated surface while the upper wall is kept impermeable, thus mimicking the midsection plane of a channel with two opposing porous walls. The measurement of instantaneous velocities is subsequently carried out at any location within the chamber using hot wire anemometers. The experimental setup resembles the present model except for the upper impervious wall, which corresponds to the plane of symmetry in the present formulation. This condition can be reproduced by simply modifying the boundary conditions on the unsteady variables and adjusting the mean flowfield.

Although various experimental results from the VECLA facility exist, the two representative cases selected here correspond to distinct chamber heights of 10 and 20 mm, respectively. The length of the nozzleless chamber is maintained at 580 mm and the wall injection speed is set at 1.36 m/s. Instantaneous velocity measurements are then taken at a distance of 1 mm away from the solid wall using a hot wire. The signal is typically treated and decomposed into a mean value and a fluctuating component. A Fourier transform is finally applied to the latter in an effort to characterize the spectral dependence of the fluctuating velocity. These findings are reflected through the power spectral density (PSD) plots displayed using a solid line in Fig. 21.

In mirroring the experimental conditions, two cases of the nonreactive gaseous injection in a nozzleless chamber are simulated using a length-to-height ratio of *L *=* *29 and 58, respectively. A suitable wall injection Mach number is set at 0.004, as it corresponds to the injection speed used in the experiments. When these values are incorporated into the compressible biglobal solver, the associated spatial and temporal scales become retrievable from the eigenvalues and eigenvectors. Moreover, to obtain the signal corresponding to the instantaneous velocity measured by the hot wire probe, some processing of the results is required. Knowing that the governing equations are linear by nature, any linear combination of the eigenvectors also leads to a valid solution. Consequently, in order to fully reconstruct the spectrum found in the signal, the eigenvectors provided by the solver can be consolidated. This can be accomplished by evaluating the weighted sum of the modal contributions using

where *n* denotes the mode number and *N _{T}* represents the total number of eigenvalues. However, this technique requires a resolution of the weighting functions, ${Wk$; $k=0,1,2,\u2026}$, which are not known beforehand. Typically, the weights

*W*can be computed by evaluating the solution at the initial state when

_{k}*t*=

*0. Presently, the initial state of the unsteady flow consists of a small perturbation that manifests itself at all possible spatial and temporal scales. Since our interest remains focused on capturing the amplified frequencies and their growth rates (notwithstanding their actual amplitudes), this issue becomes immaterial, especially that the wave amplitudes remain arbitrary in a linear solver. Furthermore, because the velocity is measured at one point for a given duration of time, the growth rate itself becomes a major contributor to the extent that a simple approximation of the function proves sufficient to extract the amplified frequencies in these simulations. Without loss of generality, we therefore set*

*W*= 1 while carrying out our modal computation.

_{k}In this manner, based on the collective modal content, amplified frequencies are readily extracted and displayed in Fig. 21(a) using dashed lines; these are provided side-by-side with the experimental results that are retraced using a solid line. Additionally, theoretically calculated Helmholtz frequencies are inserted using vertically dotted lines. These acoustic frequencies are calculated according to Eq. (38), which depends mainly on the geometric characteristics of the chamber. As one may infer graphically, the present formulation captures multiple amplifications that correspond rather favorably to the experimental measurements. These also agree with the longitudinal harmonic frequencies of the cavity. Overall, this case helps to support the hypothesis that the amplified frequencies are acoustically driven and that the hydrodynamic modal contributions remain minimal.

The second case corresponds to a much larger length-to-height ratio. The attendant chamber elongation provides the means for a stronger mean flow and, by the same token, parietal vortex shedding. Results from this case are displayed in Fig. 21(b) along with the experimental and theoretical acoustic measurements. The compressible computations are found to predict quite accurately the amplified frequencies measured in the experiments. Although these do not precisely coincide with the pure acoustic tones, they do correspond to actual vorticoacoustic modes with slightly shifted frequencies. This comparison also stands as a clear example of the importance of incorporating mean flow effects when analyzing acoustic fields in high-speed flows. Evidently, a high-fidelity mean flow profile will be useful to produce accurate and reliable simulations. At present, the sample comparisons to experimental measurements may be viewed as being supportive of the notion that a compressible eigensolver can accurately predict the chamber's amplified frequencies while providing insight into the nature of the vorticoacoustic modes.

### C. Unsteady vorticity comparisons

Evidently, it would be helpful to differentiate the contributions of hydrodynamic modes from those of vorticoacoustic modes on flow instability and breakdown structures. This can indeed be accomplished in Fig. 22 by undertaking a comparison between results that employ only vorticoacoustic modes and those that also include hydrodynamic modes. To this end, Fig. 22(a) is used to depict the vorticity contours generated from the linear superposition of the first 30 vorticoacoustic modes. These encompass longitudinal, transverse, and mixed modes affiliated with the chamber's harmonics. Graphically, one can identify a clear concentration of vortices near the wall, which appear to be of the vortex shedding parietal type. In essence, these vortices are generated near acoustic frequencies and then convected downstream by virtue of the mean flow motion. In fact, this particular behavior corresponds to what Vuillot^{85} refers to as the coupling between vortex shedding and the acoustics of the chamber. The graph provides a particularly illuminating visual depiction of how acoustic motion, propagating at the speed of sound, can modulate the frequency of the vortex shedding process, to draw energy from the mean flowfield.

To isolate the effects of hydrodynamic modes, we now move to superimpose the eigensolutions for 50 hydrodynamic modes in addition to the 30 vorticoacoustic modes used earlier. The corresponding vorticity fluctuations are then computed and presented in the contour plots of Fig. 22(b). Apart from the vortices that develop near the wall, as in Fig. 22(a), we draw the reader's attention to the particular vortical structures that emerge in the core region and continue to grow as the aft end of the chamber is approached. Since hydrodynamic modes are caused by mean flow breakdown, they tend to intensify with the mean flow speed. In this vein, a closer look at the mean flow can help to justify the appearance of these structures.

In deriving the compressible Taylor profile in a porous channel, Maicke and Majdalani^{27} discuss in detail the mean flow evolution, which reaches its maximum in the midsection plane. Therein, the flow speed increases progressively as the flow advances downstream. Consistently with these theoretical findings and the behavior of Eq. (43), the hydrodynamically generated vortices in Fig. 22(b) appear to be most prominent in the core region of the chamber, where the mean flow component is most significant, and become even more appreciable in the downstream direction.

## VII. CONCLUSION

One of the long-standing challenges in modeling the oscillatory motion in the presence of an injection-driven flowfield stands in properly characterizing the interactions that occur between the acoustic and mean flow structures, especially at the walls, and the attendant coupling with the vortical waveforms that evolve inside the chamber. It is well known that the rotational part of the wave can serve as the mechanism by which energy may be transferred from the mean flowfield for the purpose of feeding the vorticoacoustic oscillations. Typically, this coupling is accomplished by insisting on the vortical wave observing the velocity-adherence condition at the wall to the extent of producing an acoustically generated vortical field.

In this work, by solving the fully compressible and viscous equations of motion, the intrinsic coupling of the aforementioned waves is systematically captured. The modal decomposition enables us to resolve a wide range of temporal and spatial scales, including those evolving from low to high speed. When visualizing the frequency spectrum, the coupled modes distinguish themselves quite compellingly from the rest of the group by a clear amplification of their stability growth rates. These carefully identified vorticoacoustic modes appear near the natural harmonics of the cavity in a manner that is quite similar to what is observed in real firings. Granted that the bulk fluid motion carries the energy responsible for flow breakdown and self-excitation, the underlying frequency shifts are attributable, in part, to the inclusion of mean flow effects in a fully compressible framework. In this context, the presence of a mean flow slightly alters the acoustic mode shapes, which are traditionally assumed to be purely sinusoidal standing waves of the Helmholtz type.

Overall, it is gratifying that the vorticoacoustic waveforms predicted by the compressible eigensolver compare favorably with existing analytical approximations;^{42} these have been effective at identifying a characteristic parameter that controls the penetration depth of unsteady vorticity waves.^{43} By balancing the convective and dissipative terms, the non-dimensional parameter, $Mj=Vw3/(a\nu \omega r*2)$, oversees the two mechanisms responsible for the transport of the rotational boundary-layer disturbances across the mean flow.^{42} While the theoretical formulation enables us to gain deeper insight into the physical aspects of the problem, the present approach helps to extend our framework to more elaborate flow configurations with irregular boundaries with no need for flow decomposition at the forefront of the analysis.

In what concerns hydrodynamic breakdown and shear layer instabilities, it appears that a sufficiently refined compressible framework is capable of promoting an accurate resolution of most temporal and spatial scales. Consequently, the present formulation enables us to predict the modes of oscillation that are not solely related to the acoustics of the chamber. In effect, some are directly influenced by the actual characteristics of the mean flow itself. Comparisons with incompressible stability studies help to corroborate that the present solver can capture the aforementioned modes, which may become important when combustion models are appended to the solver's capabilities. Moreover, our results seem to favorably compare to experimental, albeit limited, cold-flow data acquired using the VECLA facility at ONERA; therein, the biglobal stability findings appear to outperform the LNP predictions in locating the amplified frequencies in the chamber. In essence, using cold-gas injection, the facility in question records velocity fluctuations from which the characteristic oscillatory frequencies can be deduced. According to our comparisons, their amplified frequencies seem to be vorticoacoustic in nature. In hindsight, these observations and physical insights have only become possible because of the fully compressible coupling that has been achieved between vortical and acoustic waves.

Finally, let us recall that the formulation of the vorticoacoustic wave strongly depends on the proper satisfaction of the no-slip requirement at the wall. While the overarching concept leads to an accurate representation of velocity profiles with small viscosity, pressure waves remain essentially unaffected. In this work, interesting new characteristic properties of the pressure disturbances are unraveled. More specifically, the direct coupling between acoustic and vorticity fields has enabled us to predict pressure waves that exhibit spatial shifting with respect to the routinely assumed sinusoidal form. Similar interactions caused by hydrodynamic fluctuations also affect the distribution of unsteady vorticity. These give rise to vortical structures near the midsection plane that become more appreciable in the aft region of the chamber. Since the latter constitutes the segment with the highest mean flow velocity, these structures can be associated with the early development of turbulence.

## ACKNOWLEDGMENTS

This work was supported partly by the National Science Foundation through Grant No. CMMI-1761675 and partly by the Hugh and Loeda Francis Chair of Excellence, Department of Aerospace Engineering, Auburn University. The authors are deeply indebted to Dr. Martin J. Chiaverini, Director of Propulsion Systems at Sierra Space Corporation, for numerous technical exchanges and for his unwavering support of our stability investigations. We also thank Donald Brenner, Dr. Brian Pomeroy, and Arthur Sauer for valuable discussions and technical exchanges.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## References

_{2}-O

_{2}-diluent detonations