Royal Society Publishing

On steady rotational high speed flows: the compressible Taylor–Culick profile

Joseph Majdalani


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 Embedded Image, where Mw is the wall Mach number. We solve the ensuing equations to the order of Embedded Image 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, Ls. 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 Res 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 L0. We assume that the speed of the gas at the wall Uw is uniform and L0 can reach the sonic length. In SRMs, the sidewall velocity Embedded Image is commensurate with the propellant regression rate, Embedded Image. To justify a constant Uw, we recall that the streamwise depreciation in pressure and its wall-coupling effect can nearly offset the axial decrease in density. Here, Embedded Image 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≤xL and 0≤r≤1, where L=L0/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 Embedded Image, where a0 is the speed of sound at the origin. Although we concentrate on the interesting case of choked flow at x=Ls, our solution will still apply to an isobaric opening when x<Ls. In rocketry, the wall Mach number can be as high as 0.02.

Figure 1

Porous tube accommodating sidewall injection and inert headwall conditions.

(b) Formulation

We begin by normalizing all of the remaining variables and operators using standard axisymmetric descriptors. This can be done by settingEmbedded Image(2.1)where the subscript ‘0’ refers to conditions at the origin, whereas Embedded Image alludes to the injection constant along the sidewall. Owing to its frequent recurrence, we define the second-order linear operator D2 asEmbedded Image(2.2)As usual, the velocity, Stokes streamfunction, and vorticity are related byEmbedded Image(2.3)Moreover, for steady inviscid motion, the momentum equation can be manipulated and reduced toEmbedded Image(2.4)where γcp/cv. 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 areEmbedded Image(2.5)Embedded Image(2.6)Embedded Image(2.7)Given an isentropic flow of a calorically perfect gas, one also hasEmbedded Image(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 byEmbedded Image(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 Embedded Image, we expand each variable according toEmbedded Image(2.10)The leading and first-order velocity components can be readily obtained at even powers in Mw, specificallyEmbedded Image(2.11)Embedded Image(2.12)On substituting (2.10) into (2.5)–(2.8), we collect, at Embedded Image,Embedded Image(2.13)Embedded Image(2.14)Embedded Image(2.15)Embedded Image(2.16)Expansion of (2.9) yields, in turn:Embedded Image(2.17)and so on. At Embedded Image, we segregateEmbedded Image(2.18)Embedded Image(2.19)Embedded Image(2.20)Embedded Image(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=C20. This expression is then fed into (2.14) to obtain D2ψ0+C2r2ψ0=0. An incompressible solution that fulfils the problem's constraints ensues for C=π, namely, Embedded Image. With the replacement of ψ0ψ0 by Embedded Image in the last member of (2.15), the basic momentum equation can be integrated intoEmbedded Image(3.1)The companion density and temperature are tenable directly from equation (2.16). A simple division by γ yields, for example Embedded Image.

(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×∇p1 vanishes. We are left withEmbedded Image(3.2)The scalar projection of this vector yields a single component, namely,Embedded Image(3.3)We now insert equations (2.11)–(2.12) into (3.3) and use Ω0=π20. After expanding and simplifying, we collectEmbedded Image(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 letEmbedded Image(3.5)We find that the first-order vorticity transport equation will indeed be satisfied by eliminating the residualEmbedded Image(3.6)orEmbedded Image(3.7)Equation (3.7) represents a first-order PDE in Ω1c that can be suitably integrated. After some effort, we recoup the general formEmbedded Image(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 extractEmbedded Image(3.9)To make further headway, it is expedient to apply the transformation ψ1=F(x,η)sin η. After some algebra, we holdEmbedded Image(3.10)whereEmbedded Image(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 intoEmbedded Image(4.1)where Embedded Image and Embedded Image. Equation (4.1) will be referred to as the ‘exact’ solution. Upon closer scrutiny, we identify an approximate form Embedded Image that is practically equivalent in spatial behaviour due to the dominance of the third power in x. Thus, we proposeEmbedded ImageorEmbedded Image(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, Embedded Image will be carried alongside ψ1 during the upcoming development. In actuality, Embedded Image 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 ψ/Mw. 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 Mw. 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.

Figure 2

Effect of compressibility on representative streamlines. Unless otherwise stated, we use γ=1.4 and lines or dashes to denote incompressible and compressible solutions, respectively.

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 Mw is reduced by half in figure 2b, 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 Mw and the length to reach a level of development. The plot also suggests that compressibility effects may not be large in short chambers unless Mw 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 obtainEmbedded Image(4.3)Embedded Image(4.4)Then, using Ω0=π2rx sin η, the first-order vorticity may be expressed byEmbedded Image(4.5)Recalling the classical solution, u0=(u0, v0), the approximate velocity can be assembled intoEmbedded Image(4.6)and thusEmbedded Image(4.7)The total vorticity becomesEmbedded Image(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 fromEmbedded Image(4.9)The last two expressions differ by a mere Embedded Image. 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 2a 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 Mw=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).

Figure 3

Spatial evolution of axial and radial velocities for Mw=0.01. The compressible solution leads with respect to the Taylor–Culick marker lines.

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 Ls 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=Ls, 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 a0. 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 writeEmbedded Image(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 Embedded Image. Backward substitution into the perturbed pressure yieldsEmbedded Image(4.11)Along similar lines, we findEmbedded Image(4.12)where σ=πα−2β≃12.6229. It is clear that the pressure is dominated by Embedded Image. This series can be further confirmed by calculating the centreline pressure ratio pc. In fact, we find Embedded Image to be in perfect unison withEmbedded Image(4.13)The higher order density and temperature can be evaluated from (2.21). The latter may be arranged intoEmbedded Image(4.14)where T2 is given byEmbedded Image(4.15)orEmbedded Image(4.16)The temperature, here again, is dominated by Embedded Image. Centreline temperatures yield, in turn, Embedded ImageandEmbedded Image(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.Embedded Image(4.18)Clearly, Embedded Image controls the solution. Using Embedded Image or Embedded Image, 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 haveEmbedded Image(4.19)In the same manner, by retaining two terms in the velocity, we can putEmbedded Image(4.20)orEmbedded Image(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, Embedded Image, orEmbedded Image(4.22)Note that the spurious singularity as γ→1 can be overcome using l'Hôpital's rule, with the limit being Embedded Image. 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 Embedded Image, 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 Embedded Image, 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 findEmbedded Image(4.23)or equivalently,Embedded Image(4.24)These supplement the expanded form of Χ5, namely,Embedded Image(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 Mw. 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 Ls=Χ4 or Χ3 for small and large γ, respectively. The most precise estimate Χ5 for six representative values of Mw is given in table 1. In that respect, we note that for Mw<0.005, Ls 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.

View this table:
Table 1

Critical sonic length, Ls.

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 Embedded Image, where u is the Pythagorean sum of u and v. Several lines of constant Mach number can thus be produced and displayed in figure 4a for γ=1.4 and Mw=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, Ls, differences remain indiscernible. The same attempt is repeated using, this time, γ→1 and Embedded Image. Differences are found to be so minor that they do not warrant further attention. Hence, figure 4a 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=Ls 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=Ls 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.

Figure 4

Evolution of the local Mach number over a range of Mw and γ=1.4. (a) Contour plots extend slightly past the choking distance. (b) Centreline Mach numbers and area-averaged values are compared to the one-dimensional compressible flow formula by Gany & Aharon (1999).

The centreline Mach number Embedded Image can be readily obtained from equations (4.17) and (4.18), then plotted in figure 4b. A crude but compact approximation based on Embedded Image and the first-order temperature are also shown. We implyEmbedded Image(4.26)One may infer from figure 4b that Embedded Image 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, Mm, can be numerically integrated from Embedded Image. As shown in the graph, Mm is close to 0.7 and does not reach unity until x=Lm. Following a parametric analysis, we find Lm to be a weak function of γ and practically independent of Mw. It diminishes from Lm=1.257Ls to 1.204Ls as γ is escorted across its full range. Along this excursion, Mc 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 Mm would still serve a useful purpose as it enables us to judiciously compare with one-dimensional theory. According to the latter,Embedded Image(4.27)Equation (4.27) is superimposed as the chained line in figure 4b; clearly it follows Mm 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 Mw.

(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/Ls. One may substitute back into equation (4.26) and rearrange. As confirmed by equations (4.19)–(4.24), the ensuing product πMwLs=Γ(γ)∼1 is a function of γ, particularly one in which several approximations have been unravelled in increasing order of accuracy. It can be seen thatEmbedded Image(4.28)The influence of Mw 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 revisitingEmbedded Image(4.29)Note that dependence on Mw 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).

Figure 5

Streamwise evolution of centreline pressures and temperatures using available theories. From outside moving inward, γ=1, 1.2, 1.4, 1.67. The present solution refers to the compressible Taylor–Culick approximation.

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 pc is steeper than Taylor's or Gany & Aharon's except in the last 15% stretch, as shown in figure 5a. 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 5b. 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 fromEmbedded Image(4.30)Note that Gany & Aharon (1999) project a critical pressure ratio of Embedded Image (at x=Ls) although nozzle theory foresees a critical ratio of Embedded Image. Ours predict an intermediate quantity. For example, at the sonic point, equations (4.13) and (4.17) reduce toEmbedded Image(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, Embedded Image. The most accurate areEmbedded Image(4.32)Both renderEmbedded Image(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 Embedded Image 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.

View this table:
Table 2

Critical pressure and temperature ratios. Superscripts correspond to (−) Gany & Aharon (1999); (+) one-dimensional nozzle theory.

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 6a 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.

Figure 6

Comparison of theoretical pressures and temperatures to Navier–Stokes computations obtained from two turbulent flow models.

When temperatures are compared in figure 6b, 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 findEmbedded Image(4.34)andEmbedded Image(4.35)Although not shown on the graphs, pm and Tm 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 Ls 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 Ls.

Figure 7

Sonic distance versus wall Mach number.

Considering a chamber of length Ls 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.9Ls. 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.

Figure 8

Spatial evolution of axial and radial velocity profiles up to the sonic point. Although we use a wall Mach number of 0.01 and Ls=28, relative proportions are universal. The inset rescales the axial velocity by its centreline value and compares it to CFD data (hollow squares for the kω model). The gradual flattening of the compressible profile past Ls/2 is corroborated by both theory and experiment.

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 9a, the streamwise velocity, normalized by πMwx, 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.4Ls, 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.4Ls. 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.

Figure 9

(a) Gradual steepening of the streamwise velocity as choking is approached. (b) and (c) Computational curves at the critical location. The normalized Taylor profile is axially invariant.

As shown in figure 9b, 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 Embedded Image becomes a direct measure of relative amplification. For example, the theoretical amplification at x=Ls 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 9c 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/u0. At any cross-section, the maximum local ratio occurs along the centreline, where we define ϑc(x)≡u(x, 0)/u0(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 10a 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 Mw 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=us/u0(Ls, 0), based on us=uc(Ls). In figure 10b, 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 thatEmbedded Image(4.36)By insertion of Γ, we collect Embedded Image. The asymptotic amplification for X=1 may be readily estimated fromEmbedded ImageIn applications that feature Mw>0.02, an additional correction is required, specificallyEmbedded Image(4.37)These quantities enable us to directly measure the propensity of fluid compression and its bearing on the mean flow.

Figure 10

(a) The effect of varying Mw on the compressible centreline velocity ratio as choking is axially approached. (b) The maximum velocity ratio at the sonic point.

(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, atEmbedded Image(4.38)If xϵa<(1−q)L0, 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 Embedded Image that must be surpassed for this condition to hold. This isEmbedded Image(4.39)Thus, if this motor has Embedded Image, 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 Mw=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.

View this table:
Table 3

Components of the compressible Taylor–Culick solution.

View this table:
Table 4

Approximate representation of the compressible Taylor–Culick solution.

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 considerEmbedded Image(A1)whereEmbedded Image(A2)To be specific, we letEmbedded Image(A3)and soEmbedded Image(A4)The problem is ready to be handled. The next step is to introduce a solution of the formEmbedded Image(A5)Substitution into (A 1) reveals that solvability is possible, when A0=A2=0. We recover two simultaneous ordinary differential equations,Embedded Image(A6)Embedded Image(A7)where the prime denotes differentiation with respect to η. As it is linear, the resulting set is straightforward to manage. We obtainEmbedded Image(A8)Embedded Image(A9)Here Si and Ci are the fast converging sine and cosine integral functions. They are given byEmbedded Image(A10)

Embedded Image(A11)

Our solution Embedded Image contains a set of six unknown constants, (A1,A3,C1,C2,C3,C4). This is not an overly determined system, but rather consistent with equation (2.9). Out of the four existing boundary conditions,Embedded Image(A12)only three are useful. The fourth, u1(0, r)=0, is identically satisfied by ψ1 at x=0. The three remaining statements (ac) 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,Embedded Image(A13)Its evaluation begets two equations at Embedded Image and Embedded Image that must cancel simultaneously. The first one readsEmbedded Image(A14)while the second leads to C1=0.

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


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.


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


View Abstract