## Abstract

The Taylor–Culick solution for a porous cylinder is often used to describe the bulk gas motion in idealized representations of solid rocket motors. However, other approximate solutions may be found that satisfy the same fundamental constraints. In this vein, steeper or smoother profiles may be observed in either experimental or numerical tests, particularly in the presence of intense levels of acoustic energy. In this study, we use the Lagrangian optimization principle to arrive at multiple solutions that can help to explain the practically observed motions. We then search for the extreme states that display either the least or the most kinetic energy. These are derived and found to be dependent on the chamber’s aspect ratio. By assuming a slender case, simple expressions are retrieved from which both the rotational Taylor–Culick and the irrotational Hart–McClure solutions are recovered as special cases. At the outset, a new family of flow approximations is derived extending from purely potential to highly rotational fields. These are constructed, verified and catalogued based on their kinetic energies. To help understand the tendency of a motion to shift energy states, a second law analysis is performed and used to rank these solutions based on their entropy content. Interestingly, the Taylor–Culick profile is found to represent an equilibrium state entailing the most energy and entropy among those starting from the rest. A formal extension of Kelvin’s minimum energy theorem to an open region is then pursued, followed by a direct implementation to the problem at hand. Three other flow motions with open boundaries are examined to further exemplify the application of Kelvin’s extended criteria. Our efforts illustrate the overarching harmony that exists between Lagrange’s minimization principle and Kelvin’s minimum energy theorem; both converge on the potential solution as the least kinetic energy bearer even when the velocity at the open barrier is finite.

## 1. Introduction

It may be argued that the Taylor–Culick model used to approximate the internal flow field in solid rocket motors (SRMs) stands at the foundation of a host of theoretical problems that are of fundamental interest to the propulsion community (Majdalani 2003). For example, in the study of aeroacoustic instability (Griffond *et al.* 2000; Majdalani *et al.* 2000; Chedevergne *et al.* 2006), the Taylor–Culick model has provided a mean-flow approximation about which fluctuations may be induced (Majdalani 2001*a*; Majdalani & Flandro 2002). In studying the effect of particle addition and particle mean-flow interactions, it has fallen at the epicentre of hydrodynamic instability theory (Féraille & Casalis 2003). In large-scale numerical simulations that involve particle burning and agglomeration, average speeds and accelerations within the chamber have routinely been estimated directly from the Taylor–Culick solution. One may refer in this regard to the work of Najjar *et al.* (2006), Balachandar *et al.* (2001), as well as others. In reactive flow simulations, the Taylor–Culick solution is so valuable in estimating the bulk gas motion that it has been either built into some codes or used as a benchmark to verify computations. These are illustrated by Chu *et al.* (2003) in their premixed propane–air simulation of solid propellant burning, and by Chedevergne *et al.* (2007) in their direct numerical simulations (DNS) of an idealized SRM.

The basic Taylor–Culick solution is steady, rotational, axisymmetric, inviscid and incompressible (Culick 1966). It has been modified by Majdalani & Flandro (2002) to account for viscous stresses and regressing walls, and by Majdalani & Saad (2007*b*) to comprise arbitrary headwall mass addition. By permitting variable headwall injection profiles, the extended Taylor–Culick approximation can be viewed as a viable model for idealized hybrid rocket chambers (Majdalani 2007*a*). Its accuracy has been verified in several investigations including computational (Baum *et al.* 1988), experimental (Dunlap *et al.* 1990; Barron *et al.* 2000) and theoretical studies carried out in both cylindrical and planar (slab) configurations (Majdalani & Van Moorhem 2001; Zhou & Majdalani 2002). Recently, it has been extended to non-circular cross sections by Kurdyumov (2006). It is clearly one of the most ubiquitous cold flow approximations for a full-length cylindrical motor (Apte & Yang 2000). This is especially true in applications that require an analytical mean-flow formulation. Among the relevant examples, one may cite those studies concerned with vortico-acoustic wave propagation (Majdalani & Roh 2000, 2001; Majdalani 2001*b*) and hydrodynamic instability treatment in porous chambers with (Féraille & Casalis 2003) and without particle interactions (Ugurtas *et al.* 2000; Fabignon *et al.* 2003; Chedevergne *et al.* 2006).

In this study, we revisit the procedure leading to the Taylor–Culick incompressible model. In the process, we unravel a family of approximate solutions that satisfy the problem’s fundamental constraints (see also Majdalani & Saad 2007*a*). Among those, we apply the Lagrangian optimization principle to identify the particular forms that require the most or the least kinetic energy to appear. The bracketing solutions are determined and catalogued according to their energy signature. In all cases, simple approximations are obtained assuming sufficiently long chambers. These solutions are verified and discussed in light of both an extended form of Kelvin’s minimum energy theorem and the entropy maximization principle.

## 2. Formulation

The basic rocket chamber can be modelled as a porous cylinder of length *L*_{0} and radius *a*. The headwall at the forward end is inert and the aft end is unrestricted. As shown in figure 1*a*, and stand for the radial and axial coordinates. Our solution domain extends from the headwall to the parallel, virtual nozzle attachment plane at the aft end. In what follows, we seek to approximate other solutions that may exist besides Taylor–Culick’s basic relation (Taylor 1956; Culick 1966). In particular, we hope to identify those specific solutions that require the least or most energy to excite.

### (a) Equations

A non-reactive flow may be assumed, prompted by the typically thin reactive zone above the grain surface. Following rote, the basic flow can be taken to be steady, inviscid, incompressible, rotational and axisymmetric. The inviscid equations of motion become
2.1
It should be noted that a compressible Taylor–Culick solution developed under isentropic flow conditions confirms the suitability of the present model for a variety of applications in which the effects of compressibility are small (Majdalani 2007*b*).

### (b) Boundary conditions

These are physically connected to

(i) axial symmetry and therefore no flow across the centreline, ;

(ii) vanishing parallel flow at the sidewall to secure the no-slip boundary condition, ;

(iii) an inert or open surface at the headwall, ; and

(iv) uniform injection at the cylindrical sidewall, .

### (c) Normalization

All recurring variables and operators may be normalized via 2.2 Here, and designate the fluid injection velocities at the headwall and sidewall, respectively. For steady inviscid motion, the vorticity transport equation reduces to 2.3 where 2.4 Similarly, the dimensionless boundary conditions take the form 2.5

### (d) Vorticity-stream function approach

The Stokes stream function may be introduced through
2.6
As usual, substitution into equation (2.3) requires the use of *ω*=*ω*_{θ}=*r**F*(*ψ*). One then follows tradition (Culick 1966) and selects the simplest relation between *ω* and *ψ*, namely
2.7
Despite the non-uniqueness of this expression, it permits securing equation (2.3). By inserting equation (2.7) into the vorticity equation (2.4), one eliminates *ω* and restores the PDE characteristic of this problem. This is
2.8
with the particular set of constraints being
2.9a
2.9b
2.9c
2.9d
By virtue of L’Hôpital’s rule, removing the singularity in equation (2.9a) requires securing both
2.10aand
2.10b
Equation (2.8) is then solved by the separation of variables; one finds
2.11
As equation (2.11) satisfies equation (2.10c) identically, it is sufficient to replace equation (2.9a) by (2.10a) in the forthcoming analysis.

## 3. Energy-driven solutions

### (a) Solution by eigenfunction expansions

The application of the boundary conditions must be carefully carried out, preferably in the order in which they appear. For example, equation (2.10a) gives
3.1
or *A*=0. Without the loss of generality, we set *B*=1 and rewrite equation (2.9b) as
3.2
and so ; this is satisfied by
3.3
Using *C*_{n}=(2*n*+1)*π* enables us to sum over eigenfunctions corresponding to wall suction and injection. We ignore negative integers to avoid self-cancellation. We write or
3.4
The headwall boundary condition (2.9c) implies that *β*_{n}=0, and so
3.5
Finally, the fourth condition (2.9d) becomes
3.6
Clearly, innumerable possibilities exist that will, in theory, satisfy equation (3.6), depending on the behaviour of *α*_{n}. Realizing that *α*_{n} is directly related to the sidewall injection boundary condition, we call it the ‘sidewall injection sequence’.

### (b) Kinetic energy optimization

One of the choices for *α*_{n} may be arrived at by optimizing the total kinetic energy in the chamber. The underlying principle projects that a flow may choose the path of least or most energy expenditure. To test this behaviour, we evaluate the local kinetic energy at (*r*,*θ*,*z*), for each eigensolution, using
3.7
where each mode is an exact solution bearing the form
3.8
where *η*≡(*n*+1/2)*π**r*^{2}. By assuming a system of eigensolutions with individual kinetic energies, their cumulative energy can be locally written as
3.9
Subsequently, the total kinetic energy in the chamber volume may be calculated by integrating the local kinetic energy over the length and chamber cross section
3.10
Straightforward evaluation and simplification yield
3.11
Here is the entire cosine integral. At this point, one may seek the extremum of the total kinetic energy subject to the fundamental constraint
3.12
Equation (3.12) enables us to introduce the constrained energy function
3.13
where *λ* is a Lagrangian multiplier. Equation (3.13) can be maximized or minimized by imposing . In shorthand notation, we put
3.14
Subsequently, the constrained energy function may be differentiated with respect to each of its variables to obtain
3.15
and
3.16
Equation (3.15) can be solved for *α*_{n} in terms of *λ* such that
3.17
Then, through substitution into equation (3.16), one retrieves
3.18
Finally, when *λ* is inserted into equation (3.17), a general solution for *α*_{n} emerges, specifically
3.19
With this expression at hand, the total energy given by equation (3.11) is fully determined.

### (c) Large L approximation

Given that the Taylor–Culick model is semi-infinite, it is useful to introduce a suitable form of energy density such as . Then by plotting versus *L* in figure 2, one is able to assess the energy requirements associated with our formulation as *L* is increased. We find that as the length of the chamber is incremented at fixed radius, approaches a constant asymptotic value of . A critical aspect ratio *L*_{cr} can therefore be conceived beyond which the kinetic energy will vary by less than 1 per cent from its final asymptotic value (i.e. ). For a chamber with *L*≥*L*_{cr}, one may assume an infinitely large aspect ratio in evaluating equation (3.19). For the case of SRMs with inert headwalls, *L*_{cr} is found to be approximately 6.7. Because the aspect ratio of most SRMs exceeds 20, the large *L* approximation may be safely used. A simple case may be illustrated for a simulated chamber with an aspect ratio that exceeds *L*_{cr}. In the limit as , equation (3.19) reduces to
3.20
Equation (3.20) identically satisfies the fundamental constraint expressed through equation (3.12).

### (d) Least kinetic energy solution

The optimization technique based on Lagrangian multipliers enables us to identify the problem’s extremum with no first-hand indication of whether the outcome corresponds to a minimum or a maximum. Nonetheless, the substitution of equation (3.19) into equation (3.11) enables us to compare the energy content of the present approximation to that of Taylor–Culick’s. We find that the strategy just pursued exposes the solution with least kinetic energy. For an inert headwall, the energy-minimized formulation that emerges from equation (3.20) collapses into
3.21
Note that the right-oriented arrow ‘↦’ in equation (3.21) is used to indicate that the reduced expression is valid inside the domain, 0≤*r*<1. The corresponding stream function and velocity associated with the least kinetic energy solutions are posted in table 1. Corresponding streamlines are illustrated in figure 3*a*. Using solid lines to denote the traditional Taylor–Culick’s, the steepened solution is shown using broken lines. The energy-minimized solution exhibits steep curvatures that are reminiscent of those associated with turbulent or compressible flow motions (Majdalani 2007*b*).

## 4. Generalization

### (a) Type I solutions with increasing energy levels

At this juncture, we have identified only one profile bearing the minimum kinetic energy that the flow can possibly afford. Pursuant to the original discussion by Majdalani & Saad (2007*a*), it may be hypothesized that two complementary families of solutions exist with the unique characteristics of exhibiting continually varying energy levels from which the Taylor–Culick model may be recovered. To this end, it may be useful to seek mean-flow solutions with either increasing or decreasing energies. It would also be instructive to rank the Taylor–Culick solution according to its energy content within the set of possible solutions. In the interest of simplicity, we consider long chambers and make use of equation (3.20) as a guide. Based on the form obtained through Lagrangian optimization, we note that
4.1
where *A*_{2}=8/*π*^{2} can be deduced from the radial inflow requirement given by equation (3.12). Its subscript is connected with the power of (2*n*+1) that stands in the denominator. To generalize, we posit the generic type I form
4.2
where *q*=2 reproduces the state of least energy expenditure. This relation can be made to satisfy equation (3.12) when
4.3
where is Riemann’s zeta function. Note that the *q*≥2 condition is needed to ensure series convergence down to the vorticity. Backward substitution enables us to collect the proper form of *α*_{n}, namely
4.4
The exponent *q* may be dubbed the ‘kinetic energy power index’. With the forms given by equations (4.4) and (3.11), we are able to plot the variation of the total kinetic energy versus the kinetic energy power index *q* at two chamber aspect ratios. The corresponding curves define the lower branch of figure 4. Interestingly, as , Taylor–Culick’s classic solution is recovered. In fact, using equation (4.4), it can be rigorously shown that
4.5
This result identically reproduces Taylor–Culick’s expression. All of the type I solutions derived from equation (4.4) possess kinetic energies that are lower than Taylor–Culick’s; this explains the negative sign in the superscript for the type I sidewall injection sequence. They can be bracketed between equation (3.21) and . In practice, all solutions with *q*≥5 will be Taylor–Culick-like as their energies will then differ by less than 1 per cent. The most distinct solutions will correspond to *q*=2, 3 and 4 with energies that are 81.1, 91.7 and 97.3 per cent of Taylor–Culick’s.

### (b) Type II solutions with decreasing energy levels

To capture solutions with energies that exceed that of Taylor–Culick’s, a modified formulation for *α*_{n} is in order. We begin by introducing
4.6
The key difference here stands in the exclusion of the (−1)^{n} multiplier that appears in equation (4.2). Unless this term is lumped into *B*_{q}, no solutions can be captured with energies exceeding that of Taylor–Culick’s. The remaining steps are similar. Substitution into equation (3.12) unravels
4.7
where *ζ*(*q*,*x*) is the generalized Riemann zeta function. Equation (4.7) yields
4.8
Note that the type II solutions emerging from equation (4.8) dispose of higher kinetic energies than Taylor–Culick’s. The variation of the solution with respect to *q* is embodied in the upper branch of figure 4. According to this form of , Taylor–Culick’s model is recoverable asymptotically by taking the limit as . Here too, most of the solutions exhibit energies that fall within 1 per cent of Taylor–Culick’s. The most interesting solutions are, in descending order, those corresponding to *q*=2, 3, and 4 with energies that are 47.0, 8.08 and 2.40 per cent larger than Taylor–Culick’s. When the energy level is fixed at *q*=2, a simplification follows for the type II representation. Catalan’s constant emerges in equation (4.8), namely in the form
4.9

The type II solution that carries the most energy at *q*=2 is plotted in figure 3*b* and listed in table 2. In figure 3*b*, the type II approximation is seen to overshoot the Taylor–Culick streamline curvature. In view of the two types of solutions with energies that either lag or surpass that of Taylor–Culick’s, one may perceive the case as a saddle point to which other possible forms will quickly converge when their energies are varied.

### (c) Behaviour of the velocity and vorticity fields

The sensitivity of the solution to the energy power index *q* is illustrated in figure 5 where both components of the velocity are displayed in addition to the streamline turn angle *θ*. This angle represents the slope of the local velocity measured from the radial injection direction. Making use of axial similarity, *θ* may be expressed as
4.10
Starting with the turn angle in figure 5*a*, it is clear that the flow is purely radial at the sidewall where *θ*(1)=0 irrespective of *q*. This feature confirms that the flow enters the chamber perpendicularly to the surface in fulfilment of the no-slip requirement. Conversely, for all cases considered, the establishment of strictly parallel motion along the centreline is reflected in *θ*(0)=90. Crossing the region between the wall and the centreline, the *q*=2 type I case is accompanied by the sharpest change in the turn angle from 0 to 90^{°}. However, as one moves towards the state of most kinetic energy, the smoothing process causes the turn angle to change more slowly, thus projecting smaller deviations from the radial inflow direction at intermediate positions. Physically, this behaviour is caused by a gradual increase in the radial velocity contribution with successive increases in energy content. Similar trends are depicted in figure 5*b* where the flow is the smoothest and the radial velocity is shown to be the highest at any radial position. Note that the *q*=2 type II solution also exhibits the maximum radial velocity overshoot of 16.5 per cent relative to the sidewall injection speed. This overshoot reaches its peak at *r*=0.66 and is required to compensate for the decreasing circumferential area (2*π**r**L*) normal to the injected flow. Recalling that the Taylor–Culick radial velocity exhibits a 7 per cent overshoot at *r*=0.861, the maximum overshoot calculated here is more than twice as large; it also occurs at a greater distance from the sidewall. By contrast, by examining the case of least kinetic energy in figure 5*b*, no radial overshoot is observed. Instead, the radial velocity displays its lowest absolute value by diminishing linearly from 1 at the wall to 0 at the centreline. This linear variation is accompanied by an essentially uniform axial velocity depicted for the type I *q*=2 case in figure 5*c*. Interestingly, the corresponding series approximation reduces to ** u**=−

*r*

*e*_{r}+2

*z*

*e*_{z}inside the domain. This classic potential flow solution is described by Culick (1966) and employed by Hart & McClure (1959) in modelling the internal flow in SRMs (see also McClure

*et al.*1963). It is recovered here as an extreme state with the lowest kinetic energy. It confirms Kelvin’s minimum-energy theorem which will be discussed in §7.

In figure 5*c*, it is clear that the type I axial velocities are initially blunt, with the flattest curve being the one corresponding to *q*=2. As *q* is increased, all curves evolve into a sinusoid that approaches the Taylor–Culick solution for *q*=5 and above. Further, as we cross into the type II region, the centreline velocity continues to increase with the energy content. Owing to mass conservation, , and so the centreline speed at each power index is compelled to vary with its corresponding shape to preserve *Q*. The lowest centreline speed will thus accompany the spatially uniform distribution whereas the highest speed will emerge in the narrowest and most elongated profile connected with the state of most kinetic energy. In departing from the extreme cases, *q* may be increased past *q*=5 to quickly restore the Taylor–Culick profile.

Having fully determined the velocity field, its vorticity companion may be deduced from 4.11 This expression is evaluated for the least and most kinetic energy forms and provided in table 3.

### (d) Asymptotic limits of the kinetic energy density

When the large *L* approximation is employed, with *q*=2, the type II kinetic energy density approaches a constant value of . Note that the asymptotic value for Taylor–Culick’s, , is recovered as . In general, the limit of the kinetic energy density can be written as
4.12
For the type I solutions, substitution of equation (4.4) yields a closed-form expression
4.13
In like manner, for the type II solutions, equation (4.8) leads to
4.14
Specific values of these limits are , and for type I, and , and for type II. Both types approach either from below or from above. Because these limits are quickly reached as the length of the chamber is increased, they may be suitable to simulate rocket motor flows. The energy associated with each kinetic energy power index may be inferred from figure 6. Note that the Taylor–Culick limit of 2.5838 is practically reached by both type I and type II solutions with differences of less than 0.287 and 0.265 per cent at *q*=6. Given the maximum range at *q*=2, the total allowable excursion in energy that the mean flow can undergo may be estimated at , an appreciable portion of . From an academic standpoint, the type I family of solutions bridges the gap between an essentially potential flow at *q*=2 and a highly rotational field at , thus yielding intermediate formulations with energies that vary across the range .

### (e) Pressure evaluation

One may approximate the pressure by taking
4.15
where each *p*_{n} corresponds to the pressure obtained from an eigensolution. Using equation (2.1) in normalized form, we retrieve
4.16
Then, by way of equation (3.8), the pressure may be integrated to obtain
4.17
where *p*_{0}=*p*(0,0). Interestingly, the pressure drop along the centreline collapses into
4.18
Equation (4.18) is plotted in figure 7 for *q*=2,3 and . Unsurprisingly, the largest pressure excursion is seen to accompany the type II state with most kinetic energy while the smallest pressure loss is accrued in the least kinetic energy expression, specifically, in the *q*=2 type I case (the irrotational approximation).

### (f) Arbitrary injection

For T-burners, SRMs with reactive fore-ends, and hybrid rocket chambers with injector faceplates, a model that accounts for headwall injection is required. For these problems, our analysis may be repeated assuming an injecting headwall with an axisymmetrically varying profile defined by
4.19
where *u*_{c}=*U*_{c}/*U*_{w} is the normalized centreline velocity. In the resulting expressions, *β*_{n} does not vanish. As shown by Majdalani & Saad (2007*b*), orthogonality may be applied to obtain, for an axisymmetric headwall injection profile *u*_{0}(*r*)
4.20
where *η*≡(2*n*+1)*π**r*^{2}/2. Equation (4.20) may be substituted into equation (3.4) to obtain the stream function and velocities:
4.21
Application of Lagrangian optimization in conjunction with the large *L* approximation yields identical results for *α*_{n} as those obtained in equations (4.4) and (4.8). The stream function, axial velocity and vorticity for several injection profiles are catalogued in tables 1–3 where the least and most kinetic energy solutions are identified.

### (g) Other geometric and flow configurations

The methodology presented above may be applied to other Taylor-like problems such as the inviscid rotational flow in a porous channel and the bidirectional vortex in a confined cylinder (Saad & Majdalani 2008*a*,*b*). These configurations have been thoroughly explored by Saad & Majdalani (2007*b*, 2009), Majdalani & Rienstra (2007) and Vyas & Majdalani (2006) using physical conditions that are furnished in figure 1*b,c*. Being a member of the same family of injection driven flows, these profiles exhibit physical characteristics that are similar to those ascribed to the axisymmetric problem at hand. For the porous channel flow analogue (Saad & Majdalani 2008*a*), the corresponding results are summarized in tables 4–6. In the same vein, those obtained for the bidirectional vortex in a confined cylinder are posted in table 7 (Saad & Majdalani 2008*b*). Note that all formulations in question are spawned by the Lagrangian optimization principle and exhibit two types of solutions with either increasing or decreasing energy.

## 5. Numerical verification

Our analytical expansions may be verified by solving equation (2.8) using Runge–Kutta integration. We begin by introducing the transformation *ψ*=*z**f*(*r*) through which equation (2.8) may be reduced to a second-order ordinary differential equation
5.1
To numerically reproduce the energy solutions, a judicious choice of boundary conditions for equation (5.1) is needed. It can be readily seen that
5.2
where *n* corresponds to the eigenmode associated with *C*_{n}=(2*n*+1)*π*. Subsequently, the standard boundary conditions in equation (2.9) become
5.3
Using 120 terms to reconstruct the series expansions, both numerical and analytical solutions for *f*(*r*) and *f*′(*r*) are displayed in figure 8*a,b*, respectively. This comparison is held at representative values of the kinetic energy power index corresponding to *q*=2,3 and . It is gratifying that, irrespective of the power index and mode numbers, the energy-based formulations are faithfully simulated by the numerical data to the extent that visual differences between circles (numerical) and solid lines (analytical) are masked.

## 6. Second law verification

In this section, we seek to reconcile our findings with the second law of thermodynamics and help to address questions that may arise concerning the physicality of our two families of approximations. Given the geometry at hand, we employ a control volume approach and note that, for a frictionless fluid, the entropy *s* remains invariant along individual streamlines. At the outset, one can put
6.1
As *ψ* is preset along streamlines, we express equation (6.1) in a form that is valid throughout the chamber by taking *s* to be proportional to the stream function. Choosing *σ* to denote this proportionality constant, we write
6.2
To solve for *σ*, we use the same average entropy *S*_{w} at the injecting wall as a common benchmark for all eigensolutions. With *σ*≠*σ*(*q*), we integrate equation (6.2) at the porous wall and set
6.3
The proportionality constant is thus determined as a function of the average entropy at the wall.

### (a) Entropy change for a given energy state

For a given energy state, we consider the control volume bounded by the inner surface of the porous cylinder depicted in figure 9. In the absence of heat transfer or mechanical work, the change of entropy as the fluid crosses may be determined from the net entropy flux
6.4
This result confirms that the inviscid flow corresponding to a fixed state of energy undergoes a reversible process. Note that in equation (6.4), the product of the sums is employed *viz.*
6.5
Using a different perspective, we pair the individual eigenmodes and use, assuming fully uncoupled total entropy and velocity,
6.6
Interestingly, evaluating the new difference between incoming and outgoing entropy fluxes leads to, ∀*q*≥2,
6.7
The entropy flux is thus conserved across the chamber irrespective of whether the eigensolutions remain strongly coupled or not.

### (b) Entropy change across energy states

So far, the analysis has given rise to a number of flow approximations exhibiting varying degrees of energy. A key principle that remains to be identified is the one responsible for the system to opt for one energy state over the other. In a viscous environment, entropy maximization may be invoked straightforwardly. In our case, a similar concept may be explored.

To start, we select a control volume bounded by the inner surface of the porous cylinder shown in figure 9. For the Taylor–Culick family of solutions, we evaluate the change of volumetric entropy as a function of the kinetic energy power index. Pursuant to the second law of thermodynamics, the equilibrium state is expected to carry the most entropy. This principle will be readily tested by integrating equation (6.2) and writing 6.8 For each of the family types, the variation of the total entropy with the energy power index is illustrated in figure 10. For the type I family, the least amount of entropy (78.5% of the equilibrium state) corresponds to the irrotational Hart–McClure profile. Conversely, the most entropy (114.8% of Taylor–Culick’s) marks the type II profile exhibiting the most kinetic energy. These observations confirm that the entropy in our problem is an increasing function of the kinetic energy.

#### (i) Type I

Now for the type I family of solutions, equation (6.8) reduces to
6.9
This enables us to evaluate the difference in total entropy between two energy states using
6.10
Accordingly, if one initializes the system with a given energy state *q*_{1}, Δ*S*^{−} will reach its maximum as The maximum change in entropy will thus occur between the Hart–McClure and the Taylor–Culick profiles for which
6.11
Clearly, from an entropy maximization perspective, the Taylor–Culick solution constitutes the equilibrium state among all possible type I approximations.

#### (ii) Type II

For the type II family, one may follow similar lines and write
6.12
Then using the form Δ*S*^{+}=*S*^{+}(*q*_{1})−*S*^{+}(*q*_{2})>0, *q*_{2}>*q*_{1}, the maximum change in volumetric entropy may be estimated to be
6.13
This confirms that the *q*=2 case entails the most entropy among the type II solutions.

### (c) Physicality of the type II family of solutions

The second law analysis shows that the volumetric entropy of the type I family grows with successive increases in *q* but depreciates in the type II case. Given an initial profile, the system may hence evolve according to one of the two scenarios that are described next.

#### (i) Type I branching

If the system is initialized with a type I solution, it will evolve, in the absence of other sources, to the Taylor–Culick profile with most entropy. While entropy could be further increased by switching to a type II motion, such a jump is unlikely to occur for two reasons. Firstly, the two solution branches bear widely dissimilar characteristics, including their steepness and discontinuity in the formulation of their sidewall injection sequences and . Secondly, given that the Taylor–Culick solution carries the most entropy for the type I branch, it may be viewed as a local, stable equilibrium state. Once a solution reaches this state, it will become trapped, with no tendency to continue switching branches.

#### (ii) Type II branching

If the system is initialized with a type II solution, it will tend, according to equation (6.13), to the solution with most vorticity (*q*=2). Although conditions may be conceived to trigger this mathematical outcome, it may be argued that the development may not be physically realizable as it would be difficult to initialize the flow with sufficient vorticity in the absence of external work being done on the system. The most natural flow evolution corresponds to an irrotational system originally at rest in which vorticity generation is initiated at the sidewall during the injection process. The ensuing motion will subsequently progress until it reaches its first equilibrium point that corresponds to the Taylor–Culick solution.

## 7. Kelvin’s minimum-energy theorem for an open region

Kelvin’s minimum-energy theorem predicts that the irrotational motion of an incompressible fluid in a simply connected region contains less kinetic energy than any other motion with the same normal velocity at its boundary (Thomson 1849; Hovgaard 1923). For a fluid extending to infinity, the theorem further requires a vanishing normal velocity at the far-field boundary (Batchelor 1967). From this perspective, a closed boundary will permit securing the normal velocity requirement at its surface, whereas an open boundary will not. In reference to our problem, the velocity constraint is observed everywhere except in the exit plane where an open boundary is located. On this boundary, the normal velocity varies, although its average remains constant by virtue of mass conservation. To the authors’ knowledge, the applicability of Kelvin’s theorem to problems with open boundaries has not been tested before. For this reason, we revisit Kelvin’s argument and show that it continues to hold in the presence of an open boundary provided that one of two complementary criteria is satisfied. We conclude by showing that, for the family of solutions presented here, these criteria are always met.

### (a) Proof in an incompressible medium

We let denote the potential velocity of an incompressible fluid in a volume of fluid . Then may be taken to represent the rotational contribution and difference between the velocity of any motion satisfying continuity and the potential solution . Both fields are incompressible and so
7.1
Pursuant to Kelvin’s argument, ** u** and must exhibit the same normal velocity along the boundary of or else vanish. In this context, a boundary along which the velocity requirement is breached may be called open. Using to define a surface that envelops the fluid, one may put
7.2
where and correspond to the closed and open sections of the boundary, respectively. As shown in figure 11,

**stands for the outward pointing normal unit vector. By definition, one must have, along each of these segments, 7.3 In the case of steady, homogeneous, incompressible flow,**

*n**T*and may be used to represent the specific kinetic energies associated with

**and . Then the increment of energy brought about by rotationality may be evaluated from 7.4 Recalling that and , substitution into equation (7.4) yields 7.5 The divergence theorem may be readily used to convert the last term in equation (7.5) into a surface integral. Two distinct integrals emerge and these are given by 7.6 In constructing Kelvin’s theorem, a closed region is assumed in which along the entire boundary. This permits setting and deducing that Δ**

*u**T*≥0 as for any rotational field. The identification of the potential flow as the minimum-energy bearing motion is straightforward. However, in the presence of an open boundary, the second member of equation (7.6) may be expanded and reduced to 7.7 To ensure that Δ

*T*≥0, it is necessary and sufficient to require that 7.8 Given that the first term in equation (7.8) cannot be negative, it is sufficient although not necessary to show that 7.9 Kelvin’s theorem may thus be extended to problems with open boundaries when either of the above two criteria is met. Their practical implementation may be carried out in the following order: 7.10aelse, if

*T*

_{o}<0, 7.10b

### (b) Discussion and implications

Special cases of interest may be considered for which the defining conditions in equation (7.10) are universally satisfied. For a non-vanishing rotational velocity outlet (), it is sufficient for *ϕ*>0 to be secured over the open boundary. The converse is true for an inlet where *ϕ* must be negative. As for the case of a velocity that comprises no components normal to *S*_{o}, the first criterion in equation (7.10) is met identically. Finally, if we consider a streamtube of purely rotational motion (with on the streamtube surface), then Kelvin’s theorem will be fulfiled so long as the velocity potential is either negative along all inlets or positive along all outlets (figure 12).

If, instead, we turn our attention to a multi-valued potential such as that occurring in multiply connected regions, the theorem cannot be applied except to the solution that carries the least kinetic energy among the possible potential solutions. Conversely, if we define a velocity potential as the difference between two dissimilar potential solutions bearing the same cyclic constant, *ϕ*=*ϕ*_{0}−*ϕ*_{1}, the theorem will then hold owing to the resulting potential becoming unique (Batchelor 1967).

In short, the present extension seems to indicate that Kelvin’s energy minimization is connected to the sign of an integral that depends on the rotational flux of the potential function over the open boundary. Based on the first criterion in equation (7.10), we may further infer that the theorem will be identically satisfied whenever the potential *ϕ* and the normal component of the rotational velocity exhibit the same sign at the boundaries. This extension takes us beyond the classic requirement of a vanishing or uniform velocity field at infinity. In fact, the present analysis enables us to accommodate arbitrary velocity variations at the fluid boundaries to the extent of granting Kelvin’s theorem applicability to more complex geometric settings and regions with irregular velocity distributions.

### (c) Verification

To test the validity of equation (7.10), we consider the family of solutions obtained in §3. As the type I–II velocities observe identical boundary conditions, they share the same velocity potential as well. In this case, the potential function may be returned from the Laplacian of *ϕ* over the domain bracketed by 0≤*r*≤1 and 0≤*z*≤*L*. Its essential boundary conditions translate into
7.11
The resulting *ϕ* common to all solutions is readily found to be
7.12
The corresponding rotational component may be expressed as
7.13
where *η*≡(2*n*+1)*π**r*^{2}/2. When evaluated over the boundary, equation (7.13) vanishes everywhere except in the exit plane at *z*=*L*. The criterion given by equation (7.10) may be tested by examining the sign of
7.14
It is straightforward to see that *T*_{o}≥0 as a direct consequence of
7.15
The test that *T*_{o}≥0 ensures that . The case of *T*_{o}=0 corresponds to the Hart–McClure type I solution with *q*=2 (for which ). The analysis above confirms that Kelvin’s theorem may be suitably extended to open regions and that the profile bearing the least kinetic energy is indeed the potential flow solution.

### (d) Other test cases

Several examples with open regions may be considered to which the extension of Kelvin’s theorem may be applied. From the multitude of available choices, we select the Poiseuille flow in a duct, the Taylor flow in a porous channel and the bidirectional vortex in a confined cylinder.

#### (i) Poiseuille flow in ducts of arbitrary cross sections

For this class of problems, the velocity potential corresponds to that of a uniform planar flow such that the velocity remains parallel to the channel walls. With the potential and rotational velocities being axially independent, and having a single inlet and outlet as open boundaries, it can be easily seen that 7.16 Specific velocity profiles may be obtained from Batchelor (1967) or White (1991). It is clear that the lower limit in the inequality (7.10a) is met.

#### (ii) Taylor flow in a porous channel

For the Taylor flow in a porous channel (Saad & Majdalani 2008*a*), we have
7.17
The purely rotational part may be extracted from
7.18
In this problem, the exit plane constitutes the only open boundary at *z*=*L*. Application of equation (7.10a) yields
7.19

#### (iii) Bidirectional vortex in a confined cylinder

For the bidirectional vortex (Saad & Majdalani 2008*b*), the velocity potential is a piecewise function that may be determined from
7.20
Vyas & Majdalani (2006) derived an inviscid rotational model given by
7.21
where, for simplicity, we have set *κ*=1 in the complex lamellar Vyas–Majdalani profile.

By evaluating equation (7.10a) along the open boundary at *z*=*L*, we recover
7.22
where denotes the sine integral function.

The sign of *T*_{o} confirms the applicability of Kelvin’s theorem to the problem at hand.

## 8. Convergence properties

Using the absolute convergence and ratio tests, our Fourier-like series representations can be individually shown to be unconditionally convergent for *q*≥2. The most subtle solutions to examine correspond to the type II inert headwall case with maximum kinetic energy. Its velocity and vorticity formulations require special attention. In this vein, we consider the type II stream function, specifically
8.1
The absolute convergence test may be applied to show that
8.2
where the right-hand side converges for *q*>1. In evaluating quantities that require one or more differentiations (such as the vorticity), we find it useful to substitute, whenever possible, the closed-form analytical representations of the series in question. The equivalent finite expressions enable us to overcome the pitfalls of term-by-term differentiation which, for some infinite series, can lead to spurious results. The type II axial velocity for the inert headwall configuration presents such an example at *q*=2. This series can be collapsed into a combination of inverse hyperbolic tangent functions by writing
8.3
While term-by-term differentiation of the infinite series representation of diverges, the derivative of the closed-form equivalent yields the correct outcome of
8.4
As it may be expected, the corresponding solution is accompanied by finite kinetic energy and mass flow rate despite its singularity at the centreline. The singularity in this particular case may be linked to the absence of viscosity in our model.

## 9. Conclusions

In the past four decades, the Taylor–Culick solution for a porous tube with an impermeable headwall has been extensively used in the propulsion and flow filtration communities. Recently, an extended form has been advanced in which variable headwall injection is incorporated (Majdalani & Saad 2007*b*). In this article, we show that other solutions may be obtained and these are accompanied by lower or higher kinetic energies that can vary, from one end of the spectrum to the other, by up to 66 per cent of their mean value. After identifying that yields the profile with least kinetic energy, similar type I solutions are unravelled in ascending order, , up to Taylor–Culick’s. The latter is asymptotically recovered in the limit of , a particular case that corresponds to an equilibrium state with maximum entropy. In practice, most solutions become indiscernible from Taylor–Culick’s for *q*≥5. Indeed, those obtained with *q*=2,3 and 4 exhibit energies that are 18.9, 8.28 and 2.73 per cent lower than their remaining counterparts. The least kinetic energy solution with *q*=2 returns the classic, irrotational Hart–McClure profile. It can thus be seen that the application of the Lagrangian optimization principle to this problem leads to the potential form predicted by Kelvin’s minimum energy theorem. It can also be seen that the type I solutions not only bridge the gap between a plain potential representation of this problem and a rotational formulation, but also recover a continuous spectrum of approximations that stand in between. When the same analysis is repeated using , a complementary family of type II solutions is identified with descending energy levels. These may be later proved to be purely academic, being representative of a class of exact solutions to the modified Helmholtz equation. Their most notable profiles correspond to *q*=2,3 and 4 with energies that are 47.0, 8.08 and 2.40 per cent higher than Taylor–Culick’s. These also exhibit higher entropy values than the equilibrium state. Here too, most type II solutions resemble the classic form that is engendered as . In fact, both type I and II solutions converge to the Taylor–Culick representation when their energies are augmented or reduced, respectively. Profiles belonging to these families of solutions with and *u*_{c}=1, *π*/2 and 2 for uniform, Berman and Poiseuille injection are numerically verified in a separate work by Majdalani & Saad (2007*b*). While the new solutions may be used to approximate the mean-flow profile in wall-injecting porous tubes and the bulk gaseous motion in simulated rocket motors, it should be borne in mind that no direct connection exists between their energy steepened states and turbulence. In this vein, it is hoped that additional numerical and experimental validations may be achieved to help elucidate their physicality and the particular configurations in which they are prone to appear. As for the uniqueness of the least kinetic energy state established here, we note that Kelvin’s theorem guarantees that the irrotational motion will carry less kinetic energy than any other motion in the same fluid region. This is especially true in a simply connected domain, such as ours, where the least kinetic energy profile reproduces the potential motion that is uniquely prescribed by Kelvin’s theorem and tacitly confirmed through Lagrange’s optimization principle. Along similar lines, the uniqueness of the equilibrium state may be inferred from the entropy maximization principle and the Lagrangian-based solutions in which, for a given set of equations and boundary conditions, the equilibrium state may be asymptotically restored as irrespective of the form of *α*_{n}∼(−1)^{n}(*p* *n*+*m*)^{−q}, provided that the Lagrangian constraint is faithfully withheld.

## Acknowledgements

This material is based on the work supported by the National Science Foundation through grant no. CMMI-0928762, Dr Eduardo A. Misawa, Program Director. The authors thank Professor Jie-Zhi Wu for his useful suggestion concerning the application of Lord Kelvin’s minimum-energy theorem to this problem. We are equally indebted to the anonymous referees for their constructive comments and for drawing our attention to the entropy maximization principle.

## Footnotes

- Received June 22, 2009.
- Accepted September 16, 2009.

- © 2009 The Royal Society