## Abstract

We consider the compressible flow analogue of the well-known Taylor–Culick profile. We first present the compressible Euler equations for steady, axisymmetric, isentropic flow assuming uniform injection of a calorically perfect gas in a porous chamber. We then apply the Rayleigh–Janzen expansion in powers of , where *M*_{w} is the wall Mach number. We solve the ensuing equations to the order of and apply the results up to the sonic point in a nozzleless chamber. Area averaging is also performed to reconcile with one-dimensional theory. Our solution agrees with the existing theory to the extent that it faithfully captures the steepening of the Taylor–Culick profile with downstream movement. Based on the closed-form expressions that we obtain, the main flow attributes are quantified parametrically and compared to the existing incompressible and quasi-one-dimensional theories. Verification by computational fluid dynamics is also undertaken. Comparison with two turbulent flow models shows excellent agreement, particularly in retracing the streamwise evolution of the velocity. Regardless of the Mach number, we observe nearly identical trends in chambers that are rescaled by the (critical) sonic length, *L*_{s}. Using a suitable transformation, we prove the attendant similarity and provide universal criteria that can be used to assess the relative importance of gas compression in solid rocket motors. Owing to sharper velocity gradients at the wall, we find that an incompressible model underestimates the skin friction along the wall and underpredicts the centreline speed by as much as 13% at the sonic point. In practice, such deviations become appreciable at high-injection rates or chamber aspect ratios.

## 1. Introduction

The injection-driven flow of an incompressible fluid was first described by Taylor (1956), who not only resolved this problem in porous tubes, but also in channels, wedges and cones. In lieu of using Taylor's integral equation, other investigators, such as Yuan (1956), Yuan & Finkelstein (1956), Proudman (1960), Wageman & Guevara (1960), Terrill (1964, 1965) and Culick (1966), used different routes and recovered the same cosine-shaped profile found by Taylor in the inviscid limit. Some of these studies relied on the vorticity–streamfunction approach or asymptotic tools to perturb the similarity-reduced equation derived by Berman (1953) for steady laminar flow in porous channels.

In addition to its sheer simplicity, the Taylor–Culick solution would later prove surprisingly accurate in modelling the drainage of watery suspensions across porous sheets, in surface ablation and sweat cooling (Yuan 1959; Peng & Yuan 1965), and in simulating the mean flow in solid rocket motors (SRMs). This could be attributed, in part, to its streamlines being observant of the no slip requirement at the sidewall, a feature that has inspired a few to term it ‘quasi-viscous’. It quickly resembled the viscous solution with successive increases in the injection Reynolds number (*Re*).

Following the inventive experiments by Taylor (1956) under laminar conditions, Proudman (1960), Wageman & Guevara (1960) and Yamada *et al*. (1976) confirmed the same behaviour under *Re*s that were sufficiently high to induce turbulence. Turbulence appeared to have little bearing in altering the mean flow structure except in the downstream segment of the chamber, where flow steepening was observed. In this part of the domain, compressibility intensified to the extent that it became difficult to discern which of the two mechanisms was more responsible for the steepening effect. The flattened shape was further reported in cold-flow experiments by Dunlap *et al*. (1974) and Traineau *et al*. (1986); these employed nitrogen and air injection, respectively. The spatial evolution was also parametrically substantiated via numerical simulations by Beddini (1986), Baum *et al*. (1988), Sabnis *et al*. (1989), Liou & Lien (1995) and Apte & Yang (2000, 2001). In summary, these studies confirmed the existence of a laminar segment that could stretch over the entire length of the chamber depending on the injection *Re*. They also suggested that in sufficiently long chambers, transition from the classic solution typically occurred past the midsection plane, although transition of the mean velocity was delayed to the fully developed section, where the axial profile became fuller. Both compressibility and turbulence were posited as plausible candidates for causing these departures.

The ability of the Taylor–Culick equation to faithfully imitate the motion of reactive gases ejected from the regressing wall of a solid propellant rocket has vast implications and, as such, has been the subject of much scrutiny. In one study, Majdalani *et al*. (2002) have shown that unless the burning rate is exceedingly high, the slow expansion of the sidewall may be discounted in favour of a quasi-steady model. Moreover, the idea of a uniform burning rate along the grain and the attendant assumption of constant mass addition at the sidewall have been reiterated in both nozzle-adapted and nozzleless configurations. In the latter case, Gany & Aharon (1999) have pointed out that the axially accentuated erosive-burning effects on the regression rate are generally offset by the decreasing pressure in the downstream direction. In the same vein, the nearly non-reactive environment outside the thin propellant flame zone has been repeatedly simulated and shown to be appropriately treatable by cold-flow models (Chu *et al*. 2003). This can be ascribed to the weak coupling with the energy equation, where the temperature is almost uniform (Vuillot 1995). Such factors have placed Taylor's mean flow formula at the epicentre of a large body of theoretical investigations focused on modelling rocket performance, flame zone analysis, particle–mean flow interactions and the ever-daunting task of predicting aeroacoustic instability. Relevance to internal ballistics and performance prediction in rockets is briefly outlined by Majdalani (2003). By way of example, the accuracy of Taylor's formula has given Varapaev & Yagodkin (1969) the impetus to use it as the basis of their investigation of hydrodynamic instability. In the same context, it has been employed by Griffond *et al*. (2000) to investigate the evolution of turbulence in straight porous cylinders. It has been essential not only to Féraille & Casalis (2003) in the evaluation of particle–mean flow interactions, but also to Balachandar *et al*. (2001) in establishing the connection between experimentally observed swirl and small asymmetries in the injection-driven flow.

The need for refined analytical solutions has continued to grow hand-in-hand with computer muscle for two fundamental reasons. First, to gain deeper physical insight into an ‘eternally complex problem’ and second, to provide the much needed limiting process validations for numerical simulations of chemically reactive rocket chambers. As argued by Wasistho *et al*. (2004), a team of researchers at the Center for Simulation of Advanced Rockets (CSAR) has determined that analytical models often present the only valuable resource for checking numerical results. This is partly caused by the difficulty in acquiring specific experimental data and partly due to the harsh environment erupting in rocket chambers, particularly one that is hostile to proper instrumentation. In the same study, Wasistho *et al*. (2004) have suggested the need to employ a compressible mean flow model to promote better agreement with their full-scale computations. Their steady and unsteady flow results were compared, respectively, to the numerical solution of the integral equation obtained by Balakrishnan *et al*. (1992), and to the time-dependent expression for wave motion by Majdalani & Van Moorhem (1998). Favourable agreement was noted on both counts. Today, as regulatory requirements for validation and verification continue to mount in the propulsion industry, the demands for improved analytical tools also increase. In this study, we focus on the compressible Taylor flow analogue.

The quest for a compressible mean flow expression serves multiple objectives. First, it is needed to ameliorate our representation of the undisturbed inviscid field and its impact on the ensuing acoustic wave motion and flame zone analysis. For instance, it will help to obtain a more consistent solution for the total internal flow field by insisting that both mean and unsteady components retain a compressible ingredient. Existing expressions for the time-dependent flow field are based on an incompressible steady-flow representation in both planar and axisymmetric chambers (cf. Majdalani & Roh 2000; Majdalani 2001; Majdalani & Flandro 2002). By the same token, a compressible outer solution can be indispensable in analysing the thin flame zone that forms above the surface of a solid propellant or a hybrid fuel, where pyrolysis and vaporization must be carefully modelled. It has been often reported that the exothermic viscous layer treatment cannot be completed without a compressible outer approximation (Balakrishnan *et al*. 1991). Second, it is needed to upgrade our analysis of hydrodynamic instability, which so far has been limited to incompressible equations. In simulated rocket chambers, the modelling of travelling rotational waves by Griffond (2002) is based on the oscillatory solution by Majdalani (2001). The latter relies on an incompressible mean flow representation. The use of a compressible formula at the baseline of such development can markedly improve the actual methodology and its potential outcome. Venugopal (2003) shows that a minor increase in mean flow velocity caused by density variations along the wall can engender large excursions in the growth of disturbances in a simulated SRM. Third, it is needed to provide a more precise platform for calculating the acoustic growth rate factors inherited from volume integrals of internal flow components (Majdalani *et al*. 2005). Finally, it is needed to facilitate the characterization of nozzleless rocket motors, whose increased volumetric efficiency has made them prime candidates for ramjet-boosting applications. The absence of a nozzle has enabled investigators to focus on specific features of the flow, thus providing an attractive vehicle for academic studies that employ either experimentation or computation. The operation of such motors has been examined from one- and two-dimensional perspectives by Traineau *et al*. (1986), King (1987) and Gany & Aharon (1999). It has been recently considered by Balakrishnan *et al*. (1992), who followed Traineau *et al*. (1986) in presenting an inviscid, rotational and compressible integral equation that can be solved numerically in two dimensions. Using scaling arguments, their radial momentum equation was discarded and a perturbation expansion was implemented using the reciprocal of the aspect ratio. In the process, pressure and grain regression were related by the Saint-Robert power law. The resulting ‘pseudo-two-dimensional’ approximation has proved adequate in treating long motors with arbitrary cross-section and wall-injection rate.

In this study, we reconsider Taylor's original problem in an axisymmetric, constant-area chamber with uniform wall injection. Using a Rayleigh–Janzen expansion in the wall Mach number (see Currie (2002), or the original papers by Janzen (1913) and Rayleigh (1916)), we expand the system of equations up to fourth order and extract the compressible flow analogue of Taylor's solution.

## 2. Mathematical model

### (a) Geometry

We consider the steady, inviscid and non-heat conducting flow of an ideal gas in the domain bounded by the porous sidewall of a tube of radius *a* and finite length *L*_{0}. We assume that the speed of the gas at the wall *U*_{w} is uniform and *L*_{0} can reach the sonic length. In SRMs, the sidewall velocity is commensurate with the propellant regression rate, . To justify a constant *U*_{w}, we recall that the streamwise depreciation in pressure and its wall-coupling effect can nearly offset the axial decrease in density. Here, stand for the streamwise and radial coordinates (with the overbar denoting a dimensional quantity), as shown in figure 1. We normalize these variables by *a* and select a curvilinear coordinate system, whose origin is centred at the headwall. Axial symmetry reduces the field investigation to the region 0≤*x*≤*L* and 0≤*r*≤1, where *L*=*L*_{0}/*a*. The tube is closed at *x*=0, corresponding to a zero inlet profile in Berman's equation (figure 1). SRMs exhibit wall injection Mach numbers of order , where *a*_{0} is the speed of sound at the origin. Although we concentrate on the interesting case of choked flow at *x*=*L*_{s}, our solution will still apply to an isobaric opening when *x*<*L*_{s}. In rocketry, the wall Mach number can be as high as 0.02.

### (b) Formulation

We begin by normalizing all of the remaining variables and operators using standard axisymmetric descriptors. This can be done by setting(2.1)where the subscript ‘0’ refers to conditions at the origin, whereas alludes to the injection constant along the sidewall. Owing to its frequent recurrence, we define the second-order linear operator D^{2} as(2.2)As usual, the velocity, Stokes streamfunction, and vorticity are related by(2.3)Moreover, for steady inviscid motion, the momentum equation can be manipulated and reduced to(2.4)where *γ*≡*c*_{p}/*c*_{v}. At this stage, the curl of equation (2.4) may be taken to eliminate the pressure and obtain the vorticity transport equation. After grouping and rearranging, one identifies the set of equations that must be solved. These are(2.5)(2.6)(2.7)Given an isentropic flow of a calorically perfect gas, one also has(2.8)

Our auxiliary conditions are prescribed by the continuity of the flow across the centreline and the radial inflow at the sidewall, where the streamwise component of the velocity must vanish. This assortment can be expressed by(2.9)

Equations (2.5)–(2.8) must be solved sequentially according to the following strategy. In order to satisfy the vorticity transport equation, a relation between *Ω* and *ψ* must be first determined. This relation is needed in the vorticity equation to eliminate *Ω* and render a second-order partial differential equation (PDE) in *ψ*. A solution that satisfies (2.9) may then be obtained and returned to the momentum equation, so that the pressure distribution can be deduced. Isentropic relations could then be invoked to calculate the density and temperature. This procedure must be applied after perturbing the system.

### (c) Perturbation expansion

Using a Rayleigh–Janzen series in powers of , we expand each variable according to(2.10)The leading and first-order velocity components can be readily obtained at even powers in *M*_{w}, specifically(2.11)(2.12)On substituting (2.10) into (2.5)–(2.8), we collect, at ,(2.13)(2.14)(2.15)(2.16)Expansion of (2.9) yields, in turn:(2.17)and so on. At , we segregate(2.18)(2.19)(2.20)(2.21)For this system, we invoke homogenous boundary conditions to avoid disturbing the leading order approximation, which, alone, must bear the four chief requirements.

## 3. Solution

### (a) Basic analysis

Equations (2.13)–(2.17) can be solved consecutively to yield the classic Taylor–Culick solution. To sketch this, we first recall that the vorticity transport relation may be satisfied by *Ω*_{0}=*C*^{2}*rψ*_{0}. This expression is then fed into (2.14) to obtain D^{2}*ψ*_{0}+*C*^{2}*r*^{2}*ψ*_{0}=0. An incompressible solution that fulfils the problem's constraints ensues for *C*=*π*, namely, . With the replacement of *ψ*_{0}∇*ψ*_{0} by in the last member of (2.15), the basic momentum equation can be integrated into(3.1)The companion density and temperature are tenable directly from equation (2.16). A simple division by *γ* yields, for example .

### (b) First-order vorticity–streamfunction relation

We begin by determining a pivotal relation between *Ω*_{1} and *ψ*_{1} that satisfies the vorticity transport equation. Since the momentum equation is only used to solve for the pressure, it is vital that equation (2.18) be first secured. In view of the isentropic relation given by equation (2.16), the baroclinic term ∇*ρ*_{1}×∇*p*_{1} vanishes. We are left with(3.2)The scalar projection of this vector yields a single component, namely,(3.3)We now insert equations (2.11)–(2.12) into (3.3) and use *Ω*_{0}=*π*^{2}*rψ*_{0}. After expanding and simplifying, we collect(3.4)At first glance, equation (3.4) appears to be staggeringly complex, thereby intractable by known methods. However, taking a cue from the leading order relation for the vorticity, we realize that a solution may mirror the incompressible form within a linear correction. So we let(3.5)We find that the first-order vorticity transport equation will indeed be satisfied by eliminating the residual(3.6)or(3.7)Equation (3.7) represents a first-order PDE in *Ω*_{1c} that can be suitably integrated. After some effort, we recoup the general form(3.8)where *G*(*ψ*_{0}) is arbitrary.

### (c) First-order vorticity equation

With *Ω*_{1} at hand, we turn our attention to equation (2.19). With appropriate substitutions, we extract(3.9)To make further headway, it is expedient to apply the transformation *ψ*_{1}=*F*(*x*,*η*)sin *η*. After some algebra, we hold(3.10)where(3.11)What makes this PDE interesting to solve is the freedom in selecting *G*(*x* sin *η*); it must be chosen in a manner to permit the fulfilment of the problem's boundary conditions. For the reader's convenience, this effort is consigned to the appendix.

## 4. Results and discussion

The compressible correction derived in the appendix may be rearranged and simplified with the net product taking the form of a rather straightforward expression. It collapses into(4.1)where and . Equation (4.1) will be referred to as the ‘exact’ solution. Upon closer scrutiny, we identify an approximate form that is practically equivalent in spatial behaviour due to the dominance of the third power in *x*. Thus, we proposeor(4.2)Equation (4.2) represents the simplest axisymmetric form that satisfies the problem's four boundary conditions, while exhibiting a high degree of precision in observing first principles. It accrues a negligible error in comparison with *ψ*_{1} to the extent of being graphically indiscernible from equation (4.1). Furthermore, the error continues to diminish with the distance from the headwall. Hence, it is ideally suited to model long porous tubes such as those used to simulate rocket chambers. In making our choice, the exact solution will be utilized in graphical depictions although, for added confirmation, will be carried alongside *ψ*_{1} during the upcoming development. In actuality, may serve as a crucial replacement in applications where simplicity of the mean flow formula is essential in preventing a purely numerical outcome (e.g. those involving integration).

### (a) Representative streamlines

Figure 2 compares the flow turning behaviour with and without compressibility at two-wall Mach numbers and chamber lengths. Here, we plot the characteristic *ψ*/*M*_{w}. The first observation is that compressibility causes the fluid to turn more sharply at a given location and the streamline curvature is steepened with successive increases in *M*_{w}. This effect has important implications on flow stability (Griffond 2002). In fact, it appears that fluid expansion in the downstream direction has a similar effect to that of increasing the viscosity. As illustrated earlier (e.g. Majdalani *et al*. 2002), decreasing the wall *Re* has a similar injunction on the flow. Incorporating both density changes and viscosity may therefore lead to more pronounced steepening. Viscosity, which has often been neglected in SRM simulations, may become an important factor in hybrid rocket analysis, where the wall *Re* is considerably smaller. We mention in passing that for simulated SRMs, wall expansion has been shown to exhibit the reverse effect on curvature, albeit too small to be considered except for extremely large regression rates. A purely hypothetical case occurs when the sidewall regresses at the same speed as the fluid entering the tube. Under these auspices, the expansion process can offset the effect of injection to such a degree that streamlines become normal to the sidewall.

The second observation in figure 2 concerns the striking similarity between parts *a* and *b* and the connection with the corresponding ratio of Mach numbers and chamber lengths. As *M*_{w} is reduced by half in figure 2*b*, it takes the streamlines twice the distance to reach the same level of development and disparity with respect to the Taylor solution. This trend insinuates that a direct relation may exist between *M*_{w} and the length to reach a level of development. The plot also suggests that compressibility effects may not be large in short chambers unless *M*_{w} is sufficiently high, particularly, exceeding a threshold value that may be helpful to determine. Finding the choking distance is another critical requirement for the purpose of delimiting the subsonic domain.

### (b) Velocity and vorticity fields

Velocity corrections can be derived starting from equation (2.12). We obtain(4.3)(4.4)Then, using *Ω*_{0}=*π*^{2}*rx* sin *η*, the first-order vorticity may be expressed by(4.5)Recalling the classical solution, **u**_{0}=(*u*_{0}, *v*_{0}), the approximate velocity can be assembled into(4.6)and thus(4.7)The total vorticity becomes(4.8)As usual, the vorticity vanishes along the centreline and by observing surface adherence, it is maximum at the sidewall. It can be calculated from(4.9)The last two expressions differ by a mere . This reflects the inherent accuracy that accompanies our compact formula. As it builds up with successive increases in *x*, *Ω*_{w} reaches its highest value in the exit plane.

Typical profiles corresponding to figure 2*a* are illustrated in figure 3. Here, we plot both streamwise and radial components of velocity at six axial stations, *γ*=1.4, and a characteristic Mach number of 0.01. We choose the latter because of its proximity to *M*_{w}=0.0095, a value that has been fixed in two separate studies. These include Apte & Yang (2000, 2001) who added much insight to the elegant work of Traineau *et al*. (1986). Incidentally, both groups have focused on cold-flow compressible simulations and laboratory experiments of the planar, porous channel case with air injection. We hereby select the same values in the interest of uniformity and portability. Lower injection rates from 0.0018 to 0.0036 have also been investigated by Dunlap *et al*. (1990).

In figure 3, it is clear that the discrepancy with respect to the leading order solution increases with the distance from the headwall. In both plots, this difference is hardly noticeable before crossing nearly one-third of the chamber span. Subsequently, the error in the Taylor–Culick model becomes progressively more pronounced, especially after passing *x*=15. Both streamwise and radial profiles begin to flatten, precisely as predicted by numerical simulations and experimental evidences. When *x*=25, the error between both models, taken at the centreline, reaches about 10%.

In order to better understand the factors that influence our solution, we find it imperative to calculate the point where the Mach number reaches unity. To this end, we define the length *L*_{s} to be the distance along the centreline to the point where a sonic condition is met. We base our definition on the idea of a local Mach number instead of an area-averaged value. The latter is necessitated by one-dimensional or pseudo-two-dimensional models, where information is often lumped at a given cross-section. Bearing this in mind, we remark that when *x*=*L*_{s}, only the centreline velocity component would have reached the speed of sound. At this junction, the area-averaged Mach number will be smaller than unity. In order to precisely determine the sonic point, we realize that the temperature correction must be evaluated to accurately estimate the local *a*_{0}. In this manner, the momentum equation must first be solved for the higher order pressure correction. We do this for the sake of completeness, although higher order pressures and temperatures are not likely to be large contributors except in long chambers.

### (c) Higher order pressures and temperatures

In seeking the higher order pressure term, we return to equation (2.20) whose right-hand side is now fully determined. We write(4.10)Partial integration must be carefully carried out to identify common terms in both radial and streamwise directions. It is expedient to start with the approximate solution, as the process can help in verifying the final formulation. After some algebra, we obtain . Backward substitution into the perturbed pressure yields(4.11)Along similar lines, we find(4.12)where *σ*=*πα*−2*β*≃12.6229. It is clear that the pressure is dominated by . This series can be further confirmed by calculating the centreline pressure ratio *p*_{c}. In fact, we find to be in perfect unison with(4.13)The higher order density and temperature can be evaluated from (2.21). The latter may be arranged into(4.14)where *T*_{2} is given by(4.15)or(4.16)The temperature, here again, is dominated by . Centreline temperatures yield, in turn, and(4.17)These variables are crucial for estimating the maximum Mach number at a given cross-section. It may be interesting to remark that approximate results at *η*=0 may be deduced from the exact expressions by replacing *σ* with 8.

### (d) Critical chamber length

Knowing the temperature, we can calculate the centreline Mach number based on the local speed, viz.(4.18)Clearly, controls the solution. Using or , we can solve for the streamwise coordinate *x*=Χ, at which a sonic condition is first detected. Depending on the level of approximation in the speed and temperature, we obtain estimates of increasing accuracy. At the leading order, we have(4.19)In the same manner, by retaining two terms in the velocity, we can put(4.20)or(4.21)Interestingly, a simple expression can be obtained by ‘borrowing’ the well-known critical pressure ratio from the one-dimensional nozzle theory. By choosing only two terms in equation (4.13), we can equate at the sonic point, which stands for the hypothetical nozzle throat, , or(4.22)Note that the spurious singularity as *γ*→1 can be overcome using l'Hôpital's rule, with the limit being . However, it is not simple to render the truncation order of this quantity because of the uncertainty associated with utilizing the one-dimensional theory to evaluate the critical pressure ratio. No throat of a nozzle exists here. While equation (4.13) exhibits a truncation error of , the global error accrued in criss-crossing with the one-dimensional nozzle theory remains an open question. At this point, a comparison with the most precise value of *Χ* may be illuminating.

The most accurate representation of , say *Χ*_{5}, may be derived by mustering every term that arises in equations (4.17) and (4.18). The lengthy formula we obtain is omitted here, as it may be relegated to symbolic programming. Instead, we realize that within four-digit accuracy, an equivalent expression is attainable by employing the dominant members of the temperature and velocity. We find(4.23)or equivalently,(4.24)These supplement the expanded form of *Χ*_{5}, namely,(4.25)

As it turns out, equations (4.22), (4.23) and (4.25) are the closest to *Χ*_{5}. This trend is illustrated in table 1, where the critical lengths are catalogued at typical Mach numbers and a range of specific heats. We note that the errors we entail are commensurate with the order in which they are introduced. Interestingly, the pseudo-one-dimensional approximation *Χ*_{2} leads to a surprisingly valid estimate that keeps improving at higher *M*_{w}. This may explain the ubiquitous acceptance of one-dimensional analyses, which, in some work, have been reported to exceed the expectations in predicting compressible flow behaviour. This is especially true concerning pressure distributions. Depending on the desired tolerances, we recommend the use of *L*_{s}=*Χ*_{4} or *Χ*_{3} for small and large *γ*, respectively. The most precise estimate *Χ*_{5} for six representative values of *M*_{w} is given in table 1. In that respect, we note that for *M*_{w}<0.005, *L*_{s} becomes very large so as to render it impractical in propulsive applications that exclude headwall injection. This may explain why choking through a nozzle is the rule rather than the exception in rocket-based applications.

Past the sonic point, it is uncertain whether our solution will continue to hold. Unless area expansion is permitted, weak shocks are likely to form and these are normally accompanied by irreversibilities. Of course, the latter tend to invalidate our underlying assumptions. Barring these incidences, it may be argued that the solution may continue to hold up to the point, where the area-averaged Mach number has reached unity. In fact, according to Tollmien (1941) and Kaplan (1946), the divergence of the Rayleigh–Janzen development does not occur when the local speed of sound is first exceeded at the centreline. Although we do not wish to venture further downstream, it may be useful to examine the evolution of the area-averaged Mach number for two main reasons. First, to set an upper bound on the range of validity and second, to provide supplementary tools that may be needed in sketching a meaningful parallelism with one-dimensional theory.

### (e) Area averaging

The local Mach number may be calculated from the total velocity and temperature using , where *u* is the Pythagorean sum of *u* and *v*. Several lines of constant Mach number can thus be produced and displayed in figure 4*a* for *γ*=1.4 and *M*_{w}=0.01. A few representative streamlines are also shown in the plot. What is most intriguing is perhaps our attempt to recreate the same plots using widely dissimilar Mach values of 10^{−3} and 10^{−4}. As long as the streamwise coordinate is rescaled by the appropriate sonic length, *L*_{s}, differences remain indiscernible. The same attempt is repeated using, this time, *γ*→1 and . Differences are found to be so minor that they do not warrant further attention. Hence, figure 4*a* is characteristic of the expected isocontours. In fact, they bear a striking resemblance to those illustrated in figure 3 of Balakrishnan *et al*. (1992). We show our streamlines and Mach numbers past the choking distance to illustrate their spatial evolution in the event that a shock is delayed. As mentioned earlier, only the centreline velocity would have reached sonic speed at *x*=*L*_{s} and there is little certainty beyond that point. We recall that instead of deteriorating past *M*=0.3, the Taylor model continues to hold further downstream. The compressible flow analogue may exhibit comparable resilience. Furthermore, the area-averaged Mach number at *x*=*L*_{s} may be calculated and shown to vary quite gradually and specifically, between 0.696 and 0.710 as *γ* is reduced from 5/3 to 1. This will be expounded next.

The centreline Mach number can be readily obtained from equations (4.17) and (4.18), then plotted in figure 4*b*. A crude but compact approximation based on and the first-order temperature are also shown. We imply(4.26)One may infer from figure 4*b* that is an adequate representation except for the slight overshoot that it displays near the aft end. As *γ* is decreased, this discrepancy becomes less conspicuous and vanishes, eventually, near unity.

To set the stage for comparisons with one-dimensional theory, the area-averaged Mach number, *M*_{m}, can be numerically integrated from . As shown in the graph, *M*_{m} is close to 0.7 and does not reach unity until *x*=*L*_{m}. Following a parametric analysis, we find *L*_{m} to be a weak function of *γ* and practically independent of *M*_{w}. It diminishes from *L*_{m}=1.257*L*_{s} to 1.204*L*_{s} as *γ* is escorted across its full range. Along this excursion, *M*_{c} crawls from 1.341 to 1.390. We suspect that full choking in a straight porous tube may be delayed an extra 20–25% of the sonic length. If the opposite were proved true, then the evaluation of *M*_{m} would still serve a useful purpose as it enables us to judiciously compare with one-dimensional theory. According to the latter,(4.27)Equation (4.27) is superimposed as the chained line in figure 4*b*; clearly it follows *M*_{m} except for the sudden sprint in the last 10% stretch. This result confirms the ability of our solution to mimic one-dimensional analysis except near the sonic point, where radial flow variability comes into play. Our results also support the idea of a local Mach number being independent of *M*_{w}.

### (f) Key geometric similarity

Inspired by the form of equation (4.27), a more convincing proof of this behaviour may be pursued by introducing a rescaled variable, *X*=*x*/*L*_{s}. One may substitute back into equation (4.26) and rearrange. As confirmed by equations (4.19)–(4.24), the ensuing product *πM*_{w}*L*_{s}=*Γ*(γ)∼1 is a function of *γ*, particularly one in which several approximations have been unravelled in increasing order of accuracy. It can be seen that(4.28)The influence of *M*_{w} is clearly negligible. Furthermore, the dependence on *γ* is also feeble in view of the slowly varying filter *Γ*. The latter shifts by a total of 11%, from 0.8846 to 0.7860, over the full range of *γ*. In rocketry, *γ* seldom exceeds 1.2; the corresponding 3.78% deviation in *Γ* is marginal, being from 88.5 to 85.1% only. As *γ* barely changes in a given application, the behaviour of our solution with respect to *X* is nearly ‘frozen’. This finding is pivotal as it enables us to render the results in a more universal fashion, namely, independently of the wall Mach number. The same idea applies to the streamfunction and its derivatives. This can be demonstrated by revisiting(4.29)Note that dependence on *M*_{w} disappears in a suitably normalized variable. Being a parent function, all related family members inherit this property. This also explains and settles our earlier bewilderment regarding the nature of figure 3, where in hindsight, the geometric scaling required to achieve similarity was met (unintentionally). The critical lengths for the two cases were actually 26 and 52, respectively.

### (g) Critical pressure and temperature ratios

Unlike the local Mach number that can undergo wide excursions at a given cross-section, we find the pressures and temperatures to be quasi-uniform. This can be quickly inferred from the favourable comparisons that may be drawn between any two of the centreline, sidewall or area-averaged properties. Owing to their similarities, we focus hereafter on centreline behaviour (see figure 5). Our axisymmetric results can also be examined in the context of existing theories. Among those, we choose the often-cited one-dimensional model of Gany & Aharon (1999), the pseudo-two-dimensional model of Traineau *et al*. (1986) and the incompressible formula of Taylor (1956). In addition, owing to the variability that exists among available data, we opt to conduct our own numerical simulations of the porous tube. In the interest of uniformity, we replicate conditions that mimic those implemented by Traineau *et al*. (1986) and subsequently by Apte & Yang (2001).

Our porous tube is chosen to be 28 cm long with a 1 cm radius and a wall Mach number of 0.0095. We apply a mass flux of 13 kg m^{−2} s^{−1} at an injection temperature of 260 K and *γ*=1.4; the fluid is air with a molecular weight of 29 kg kmol^{−1} and a dynamic viscosity of *μ*=1.66×10^{−5} kg m^{−1} s^{−1}. Unlike Traineau *et al*. (1986), we do not attach a divergent nozzle at the aft end. Our chief concern is with the faithful replication of the geometry used in the mathematical derivation. We recognize that area expansion can expedite convergence, but opt to adhere with Taylor's model, albeit at a small computational expense. Our mesh consists of (100, 280) elements in the radial and axial directions, respectively. Grid independency is achieved by ensuring consistency after doubling the mesh resolution.

The streamwise variation of *p*_{c} is steeper than Taylor's or Gany & Aharon's except in the last 15% stretch, as shown in figure 5*a*. It becomes more pronounced with successive increases in *γ*. The incremental steepening reflects inevitable departures from the hypothetical case for which *γ*→1. The same can be said for the temperature variation in figure 5*b*. In both cases, the one-dimensional model accrues a smaller drop along a major portion of the chamber. It then briskly accelerates in the aft 15% of the length on its way to reclaiming its choke-point value. Its evolution may be traced from(4.30)Note that Gany & Aharon (1999) project a critical pressure ratio of (at *x*=*L*_{s}) although nozzle theory foresees a critical ratio of . Ours predict an intermediate quantity. For example, at the sonic point, equations (4.13) and (4.17) reduce to(4.31)where *Γ* may be directly extracted from equations (4.23), (4.24), or yet (4.22), if a blend with one-dimensional theory is to be entertained. The latter yields the simplest, . The most accurate are(4.32)Both render(4.33)Interestingly, the Taylor series expansion of *Γ* matches the *γ*-dependence of the highest order expression for *Χ*_{5} encapsulated by equation (4.25). Furthermore, the term in equation (4.31) may be safely ignored without affecting the three-digit accuracy in the algebraic outcome. For the reader's convenience, a summary of critical pressure and temperature ratios is given in table 2.

Over the range of *γ*, our pressure ratio depreciates from 0.583 to 0.459. This constitutes a 16.6–22.4% increase in pressure recovery with respect to one-dimensional theory. The increased pressure recovery may be attributed to the absence of viscosity in our model. If viscosity had been accounted for, a larger pressure drop would have been entailed. The one-dimensional model appears to be less vulnerable to friction and attendant irreversibilities: it integrates those away in the process of injecting the flow directly along the axis, where viscous damping is the least prominent. If it had been viscous, the present model would have been subject to irreversibilities that consume their share of thermal energy and increase both the temperature and internal energy of the flow. The pressure drop, which is partly a measure of these irreversibilities, would have been larger. The accentuated sensitivity of our axisymmetric representation to viscous losses can, therefore, explain the larger pressure recoveries observed relative to viscous simulations or one-dimensional approximations.

For example, using *γ*=1.4, our results are compared in figure 6*a* to computed pressure curves obtained under two widely accepted turbulent flow models. While the *k*–*ω* curve nearly coincides with our analytical solution over a significant portion of the chamber, it begins to slightly diverge near the aft end to the extent of closing at 0.400. The Spalart–Allmaras simulation is slightly higher near the headwall, but is terminated by the same value. The increased pressure drop compared to ours at 0.501 may be linked to viscous damping. Naturally, the computed value is nearest to the one-dimensional critical ratio (0.417) proposed by Gany & Aharon (1999). However, until the exit section, the present solution adequately performs in matching the *k*–*ω* curve.

When temperatures are compared in figure 6*b*, an equally interesting behaviour is captured. The turbulent flow computations seem to fall, as predicted, slightly above our theoretical curve; they specifically coincide with the one-dimensional solution over a major portion of the chamber. In contrast to the pressure behaviour, a better agreement with our model is seen near the aft end, where the numerical curves terminate at 0.845 and 0.843, respectively. This puts them nearest to our theoretical value of 0.835. A better accord with numerical simulations would thus require the participation of viscosity.

Before leaving this subject, it may be instructive to evaluate the area-averaged quantities as they can provide expeditious estimates. We find(4.34)and(4.35)Although not shown on the graphs, *p*_{m} and *T*_{m} vary in concert with their counterparts taken along the axis. This reaffirms the frail pressure and temperature gradients in the radial direction. It also justifies the appropriateness of earlier comparisons between centreline properties and area-averaged one-dimensional projections.

### (h) Characteristic effects of compressibility

Having established the importance of *L*_{s} as a geometric scaling parameter, we can use it to delineate the zone of validity for the subsonic region. We show this in figure 7 at four representative values of *γ*. The minute shifts in the sonic curves reflect the weak sensitivity of the solution to the ratio of specific heats. This reinforces the notion of recouping a nearly frozen state after rescaling with respect to *L*_{s}.

Considering a chamber of length *L*_{s} as our baseline case, the spatial evolution of the velocity field across the entire domain may be captured independently of the Mach number. This is illustrated in figure 8, where both streamwise and radial components are plotted in the upper and lower parts of the graph, respectively. These represent scaled-down versions that preserve the relative proportions of actual magnitudes. The former plot is repeated in the middle part, where it is magnified for better clarity (using the local centreline speed) and compared to computational fluid dynamic (CFD) predictions. These are shown at nine, equally spaced axial stations corresponding to 0.1, 0.2, 0.3,…, and 0.9*L*_{s}. Although we use operating parameters that co-tail our computational model, we emphasize the geometric similarity that renders this plot universal. Its response to changing *γ* is so minor that it may be ignored. This is verified by modifying the Mach number from 0.0001 to 0.01 with no distinct influence on the solution. Changing *γ* from 1 to 5/3 has some influence, but it is too small to deserve particular attention.

As expected, these profiles bear a striking resemblance to the numerical results and to both laboratory and computational experiments obtained by Traineau *et al*. (1986), Balakrishnan *et al*. (1992), Apte & Yang (2000, 2001) and others. Specifically, we note that the streamwise velocity becomes much fuller than Taylor's sinusoidal profile as choking is approached. This is accompanied by the classic linearization of the radial velocity. The evolution into a blunter, turbulent-like or pseudo-one-dimensional plug flow is conformant to both theory and experiment. It faithfully captures the increased gradients at the sidewall and these can have important implications in mean flow-related analyses. To the author's knowledge, this represents the first explicit analogue of Taylor's compressible flow that exhibits the correct steepening behaviour. Its agreement with CFD predictions is excellent despite the presence of small viscosity in the numerical simulation.

In order to quantify the steepening effect and magnification caused by compressibility, theoretical and CFD profiles are overlaid in figure 9. In figure 9*a*, the streamwise velocity, normalized by *πM*_{w}*x*, is shown at four evenly spaced distances from the sonic point. Since the leading order sinusoidal profile varies between 0 and 1, it can be conveniently adopted as a benchmark. The amplification of the streamwise velocity can thus be inferred relative to Taylor's model. We note that below 0.4*L*_{s}, the compressible solution and Taylor's model overlap (the discrepancy between them being less than 2%). This justifies the use of Taylor's model in chambers whose actual length is roughly *L*≤0.4*L*_{s}. The precise cut-off point is a matter of conjecture, but we shall attempt to define it based on the relative error with respect to Taylor's model. So first, this relative error must be characterized.

As shown in figure 9*b*, the maximum difference between the axial velocities, which appears at the centreline, can be evaluated. Being commensurate with the local impact of fluid compression, the ratio becomes a direct measure of relative amplification. For example, the theoretical amplification at *x*=*L*_{s} can be determined to be 1.112 versus a computed value of 1.133; the CFD limit is converged upon by each of the constituent models. Despite selecting the (most) critical location, a nearly perfect match with the turbulent profiles is realized in both fullness and extent. The turbulent simulations are performed using two standard schemes, namely, the *k*–*ω* and Spalart–Allmaras (B. A. Maicke, personal communication, 2005). These yield such similar distributions that their differences are graphically imperceptible. The concurrence of the analytical solution with turbulent simulations is gratifying, especially that the model is often applied to problems involving large cross-flow *Re*. In fact, laminar simulations with a very small viscosity are also obtained and compared in figure 9*c* to analytical and *k*–*ω* predictions. Although the same velocity overshoot is realized at the centreline by both laminar and turbulent models, the bluntness of the turbulent simulation more closely resembles the analytical profile. The affinity with turbulent predictions is not an artefact, but rather an attribute of inviscid flows.

Overall, the amplification at a radial position can be measured from the ratio of velocities, *ϑ*(*x*, *r*)≡*u*/*u*_{0}. At any cross-section, the maximum local ratio occurs along the centreline, where we define *ϑ*_{c}(*x*)≡*u*(*x*, 0)/*u*_{0}(*x*, 0). The smallest increase in *ϑ*_{c} can have an appreciable impact on core flow ingredients, including the growth rate of disturbances. As recently shown in a compressible flow simulation of a model rocket motor, a roughly 10% increase in centreline speed (caused by density variations) can lead to 70% overshoot in the growth of oscillatory wave amplitudes at a streamwise location of 40 (cf. Venugopal 2003, pp. 61–65). The evolution of *ϑ*_{c} with distance from the headwall is illustrated in figure 10*a* over a wide range of operating parameters and up to the sonic point. The latter occurs sooner at higher *γ*. Since the critical length depends on the gas compression ratio, dotted lines are used in the buffer region that is bracketed by the upper and lower bounds for *γ*. As *M*_{w} is lowered, it appears that for the chosen range of Mach numbers, the relative amplification ratio increases to a pure constant at fixed *γ*. This observation prompts us to define the maximum overall velocity amplification, *ϑ*_{s}=*u*_{s}/*u*_{0}(*L*_{s}, 0), based on *u*_{s}=*u*_{c}(*L*_{s}). In figure 10*b*, a plot of *ϑ*_{s} confirms our hypothesis. We find that the sonic amplification ratio does indeed asymptote to a constant value in the practical range of Mach numbers. The asymptotic limit varies between 1.130 for *γ*→1 and 1.103 for *γ*=5/3. Analytically, the corresponding magnification may be obtained from equation (4.18) such that(4.36)By insertion of *Γ*, we collect . The asymptotic amplification for *X*=1 may be readily estimated fromIn applications that feature *M*_{w}>0.02, an additional correction is required, specifically(4.37)These quantities enable us to directly measure the propensity of fluid compression and its bearing on the mean flow.

### (i) Extension to solid rockets

One may wonder under what circumstances density variations will become important in SRMs. The present study can help to answer this question given a desired level of precision set by the designer. Suppose an error that exceeds *ϵ*=4% over more than *q*=1/4 of the chamber length is deemed unacceptable. Our solution enables us to calculate the distance to reach such degree of disparity with respect to Taylor's equation. This is met when *ϑ*_{c}=1+*ϵ*, or to good approximation, at(4.38)If *x*_{ϵ}*a*<(1−*q*)*L*_{0}, then volumetric expansions must be accounted for. By applying this formula in reverse to a chamber with a relatively large aspect ratio, say *L*=50, we find the (minimum) cut-off Mach number that must be surpassed for this condition to hold. This is(4.39)Thus, if this motor has , a 4% error will affect more than its aft quarter length. Otherwise, it may be safely treated using Taylor's formula. Considering that many SRMs operate under *M*_{w}=0.0042 and have lengths that are shorter than *L*=50, the suitability of Taylor's model remains unchallenged in those particular instances. The present solution becomes necessary, however, in applications that exhibit accentuated sensitivity to the mean flow. Examples include the full-scale numerical simulations conducted by CSAR, the treatment of the Stokes acoustic layer with sidewall injection (cf. Wasistho *et al*. 2004), the investigation of instability (cf. Venugopal *et al*. 2001), and the analysis of the exothermic boundary layer (cf. Balakrishnan *et al*. 1992). Surely, when precision is required to mitigate unwanted deviations, the present solution may be resorted to. Overall, it increases our arsenal of physical approximations for the treatment of high speed compressible flow.

## 5. Conclusions

In this study, the compressible flow analogue to the Taylor–Culick problem is solved analytically and presented in two closed forms: one satisfying all first principles and one approximate, retaining the essence of the solution (see tables 3 and 4). These key findings lead to a formidable leap in our physical understanding of compressible flow behaviour in multi-dimensional settings. In propulsion-rooted applications, they permit the re-evaluation of many ballistic performance measures, which, so far, have been mostly confined to one-dimensional platforms.

The most significant result is perhaps the advent of a compact analytical expression for the streamwise and wall-normal velocity profiles, which unlike other models, can be calculated easily. It is gratifying that our solution can reproduce the fullness and extent of Navier–Stokes profiles computed using two reliable models of turbulence up to the sonic point. The second significant result is the discovery of choking length as the appropriate scale leading to nearly frozen behaviour irrespective of the wall Mach number. This has enabled us to recast the solution in a more portable, parametric-free form. In the wake of this generalization, simple design criteria that quantify the relative importance of gas compression have been deduced and shown to realistically project actual behaviour. Owing to the particular role it plays, the sonic length has been derived and presented using several approximations of increasing precision. The amplification of the centreline velocity, a key parameter in core flow analysis, has also been characterized and tacitly secured.

By circumventing the need to compute the mean flow at each point and set of operating variables, our compact solutions greatly facilitate parametric trade analyses. As it is compressible, these can be used in lieu of the classic Taylor–Culick profile to provide a better foundation for studying the growth of unsteady disturbances in injection-driven enclosures. Applications include the investigation of thermo-acoustic and hydrodynamic instability mechanisms in solid rockets. The fact that our streamlines are steepened at higher Mach numbers and fixed length lends support to Griffond's and Venugopal's recent hypotheses; these ascribe the increased stability sensitivity to steeper streamline curvatures. In light of this consensus regarding the destabilizing role of fluid dilatation, further exploratory studies may be warranted.

On the computational side, the present formulation enables us to devise a quick-starting numerical scheme through which the analytical solution may be used at start-up to initialize the problem's incumbent variables. The resulting simulations are likely to exhibit accelerated convergence. A similar methodology may be employed in problems with axial asymmetry and variable cross-section.

The small differences we find in pressure and temperature estimations are linked to the absence of viscosity in our model. This motivates the quest for a two-pronged viscous and compressible mean flow solution. Incorporating friction may be useful in modelling hybrid rocket core flows in which the wall *Re* is small in comparison to SRMs. The use of a refined compressible model is hoped to improve our strategy in analysing unsteady vortico-acoustic wave propagation and stability in porous enclosures. Repeating the analysis with constant mass-flux surface conditions may also provide a more appropriate model vis-à-vis some existing experimental and numerical simulations. Finally, having demonstrated the viability of the present methodology, the work may be extended to other geometric settings and non-isentropic flow environments in which the energy equation may be employed. In short, a plethora of research topics are now open to investigation.

## 6. Appendix

We consider(A1)where(A2)To be specific, we let(A3)and so(A4)The problem is ready to be handled. The next step is to introduce a solution of the form(A5)Substitution into (A 1) reveals that solvability is possible, when *A*_{0}=*A*_{2}=0. We recover two simultaneous ordinary differential equations,(A6)(A7)where the prime denotes differentiation with respect to *η*. As it is linear, the resulting set is straightforward to manage. We obtain(A8)(A9)Here Si and Ci are the fast converging sine and cosine integral functions. They are given by(A10)

(A11)

Our solution contains a set of six unknown constants, (*A*_{1},*A*_{3},*C*_{1},*C*_{2},*C*_{3},*C*_{4}). This is not an overly determined system, but rather consistent with equation (2.9). Out of the four existing boundary conditions,(A12)only three are useful. The fourth, *u*_{1}(0, *r*)=0, is identically satisfied by *ψ*_{1} at *x*=0. The three remaining statements (*a*–*c*) yield six constraints in integer powers of *x* that must vanish independently. Therefore, two additional constants (one could have expected four) are essential. As noted by Balakrishnan *et al*. (1992), this problem is quasi-parabolic. To illustrate this behaviour, we first consider (A 12)*a*, specifically,(A13)Its evaluation begets two equations at and that must cancel simultaneously. The first one reads(A14)while the second leads to *C*_{1}=0.

To ensure no slippage and uniform inflow along the sidewall, we must enforce *u*_{1}(*x*, 1)=*v*_{1}(*x*, 1)=0. In terms of the streamfunction, we now write(A15)This set engenders four equations in descending orders of *x*. By cancelling the , and terms, we retrieve *A*_{3}=4/3, *C*_{2}=−*π*^{2}/4 and . Finally, at we reap(A16)Equation (A 14) may be revisited to pick . This enables us to calculate the last constant from (A 16), namely, . Thus, we conclude our derivation of an axisymmetric solution that fully satisfies equations (2.18) and (2.19).

## Acknowledgements

This project was funded by the National Science Foundation through Grant No. CMS-0353518. The author is deeply grateful to Drs Mario A. Rotea and Masayoshi Tomizuka, NSF Program Directors, for supporting this work.

## Footnotes

- Received January 7, 2006.
- Accepted June 20, 2006.

- © 2006 The Royal Society