**Internal Flows Driven by Wall-Normal Injection**

Joseph Majdalani and Tony Saad

*University of Tennessee Space Institute USA*

#### **1. Introduction**

14 Will-be-set-by-IN-TECH

94 Advanced Fluid Dynamics

[51] Yanagawa, T.K.B., Nakada, M., Yuen, D.A.: A simplified mantle convection model for thermal conductivity stratification. Physics of the Earth and Planetary Interiors 146, pp.

[52] A.B. Ezersky, A. Garcimartín, J. Burguete, H.L. Mancini and C. Pérez García. *Hydrothermal waves in Marangoni convection in a cylindrical container*. Phys. Rev. E 47, pp. 1126-1131, 1993; A.B. Ezersky, A. Garcimartín, H.L. Mancini and C. Pérez García.

163-177 (2004).

*ibid*. 48, pp. 4414, 1993.

The internal motion through porous chambers generated by wall-normal injection has received considerable attention in the second half of the twentieth century. This may be attributed to its relevance to a large number of phenomenological applications. In actuality, the motion of fluids driven by either wall injection or suction can be used to describe a variety of practical problems that encompass a wide range of industries and research areas. To name a few, these include: paper manufacturing (Taylor, 1956), ablation or sweat cooling (Peng & Yuan, 1965; Yuan & Finkelstein, 1958), boundary layer control (Acrivos, 1962; Libby, 1962; Libby & Pierucci, 1964), peristaltic pumping (Fung & Yih, 1968; Uchida & Aoki, 1977), gaseous diffusion or filtration, isotope separation (Berman, 1953; 1958a;b), irrigation, and the mean flow modeling of both solid (Culick, 1966; Zhou & Majdalani, 2002) and hybrid rockets (Majdalani, 2007a).

Wall injected flows are initiated by the injection or suction of a fluid across the boundaries of a ducted region having an arbitrary shape and cross-sectional area. This is illustrated in Figure 1 for the special cases of porous channels and tubes. In general, one is required to solve a reduced-order form of the equations of motion for a bounded fluid in order to retrieve a meaningful solution (Terrill & Thomas, 1969). For a general three dimensional setting, this effort leads to a formidable task that is often intractable. However, when simplifying assumptions are invoked, as in the case of an incompressible stream in a channel or tube with uniform injection or suction, Berman (1953) has shown that the Navier-Stokes equations can be reduced to a fourth order nonlinear ODE that may be susceptible to both analytical and numerical treatment. Berman's approach is based on a spatial similarity that transforms the Navier-Stokes equations to a more manageable ODE by assuming that the transverse velocity component *v* is axially invariant; this immediately translates into a streamfunction that varies linearly in the streamwise direction, i.e. *ψ*(*x*, *y*) = *xF*(*y*) (Berman, 1953; White, 2005). Then by considering the limiting case of a small suction Reynolds number, Re ∼ *ε*, Berman employs a regular perturbation series in Re to obtain an approximate expansion for the mean flow function *F*(*y*). Berman's Reynolds number, Re = *U*w*a*/*ν* , is based on the injection speed at the wall, *U*w, and the channel half height, *a*. As for the case of large suction, Berman (1953) first remarks that the limit of the reduced ODE cannot be used to obtain a solution owing to the reduction in order of the governing equation. Later, Sellars (1955) and Terrill (1964) invoke a procedure that permits the extraction of a closed-form analytical approximation for the large Re case by implementing a coordinate transformation that takes into account the spatial relocation of the boundary layer to the sidewall region.

Reynolds number. However, Yuan's model suffered from a singularity that appeared in the third derivative of the mean flow function *F*(*y*) taken at the centerline. This of course signaled the presence of a thin boundary layer that necessitated special treatment. The corresponding boundary layer would later be captured by Terrill (1965) who also described an insightful

Internal Flows Driven by Wall-Normal Injection 97

In the interim, Berman (1958a) published his second work in which he extended the original planar problem to various geometric settings. This included the familiar case of a straight axisymmetric tube with permeable walls. Almost concurrently, White et al. (1958) advanced a series approximation to the porous channel problem for all ranges of the Reynolds number. However, White and co-workers employed a power series expansion that was centered around Re = 0. They also supplied a numerical solution to this problem. Despite the accuracy of their technique, their power series depended on two arbitrary constants that could only be determined numerically through a trial and error procedure. According to Terrill (1964), their method could be viewed as suitable for intermediate values of Re (15 ≤ Re ≤ 35). Otherwise, a transformation of the governing equation could be more effective at achieving direct numerical integration. Due to the penalty involved in evaluating the analytical constants of the attendant power series, this particular approach would be later abandoned. Nonetheless, it remained somewhat unique in its ability to provide a single analytical approximation that applied over the entire range of Re, a feat that standard perturbation methods failed to

Along similar lines, Terrill (Terrill, 1964; 1965) compiled a comprehensive and detailed résumé of the perturbation solutions of this problem over all ranges of the Reynolds number. Therein, he derived and discussed several limiting cases such as Re = 0, |Re| � 1, Re → +∞, Re → −∞, and compared the various solutions with numerical simulations based on Runge–Kutta integration. For the numerical integration scheme, he introduced a transformation that would lead to a direct numerical solution with no need for predictor-corrector steps or shooting. On the flip side, his technique did not allow the pre-selection of the Reynolds number but rather the post-determination of Re at the conclusion of the numerical procedure. Before leaving this topic, we also note the work of Eckert et al. (1957) who, as far as the authors could verify, were the first to present a numerical solution for the laminar viscous motion in a porous channel. As far as stability is concerned, the variety of analytical models considered for the planar case appeared to be both unique and stable (Terrill & Thomas, 1969; White, 2005). However, Robinson (1976) reported that dual solutions could exist for large suction while Zaturska et al. (1988) furnished a detailed stability analysis that rigorously showed that (at least) three types of solutions could co-exist. Even more intricate structures would arise in the case of axisymmetric flow in a porous tube. In this context, Terrill & Thomas (1969) have shown that, at least, dual solutions existed for the entire range of injection and suction Reynolds numbers while no steady solutions could be identified for 2.3 < Re < 9.1. At the time of this writing, the issue of stability of wall-injected flows remains an open area of investigation especially

In propulsive applications involving solid and hybrid rocket motors, modeling the mean flow proves to be important for a variety of reasons (Culick, 2006). The instantaneous flow field plays a key role in describing acoustic instability, particle-mean flow interactions, erosive burning, nozzle erosion, and thrust performance. The traditional modus operandi is to decompose the instantaneous motion into a steady average flow and an amalgam of unsteady

technique to solve this problem numerically.

among applied mathematicians and fluid dynamicists.

**1.1 Relevance to propulsion systems**

accomplish.

Fig. 1. Schematics of porous channels and tubes in which motion is sustained through wall-normal injection.

It is widely believed that Berman (1953) was among the earliest to examine the problem of laminar viscous flow bounded by porous surfaces (see Dauenhauer & Majdalani, 2003; Zhou & Majdalani, 2002). Although his first similarity transformation only applied to a planar configuration with wall suction, it has set forth the foundation for a number of follow-up investigations that relied on either analytical or numerical techniques to explore a variety of geometric configurations with either injection, suction, or both (Proudman, 1960).

Chronologically, these start with Sellars (1955) who extended Berman's solution to very large suction Reynolds numbers. He accomplished this by relaxing the no-slip boundary condition that became immaterial under this limiting condition. At the outset, he extracted a leading order approximation that corresponded to uniform axial motion, i.e. *F*(*y*) = *y*. Sellars integrated the ensuing equation based on his leading order approximation. His model thus uncovered the outer solution of this problem when viewed from a boundary layer perspective. Sellars' identification of a thin boundary layer at the wall for the large suction case would later prove crucial in subsequent developments of this problem.

Of particular interest to this chapter is a classic article by Taylor (1956) in which he derived an inviscid rotational solution for both planar and axisymmetric channel flow configurations, in addition to cones and wedges. The absence of viscosity in his model led to approximations that were consistent with Berman's leading order solution of the Navier-Stokes equations expressed at large injection Reynolds numbers, i.e. *F*(*y*) = sin( 1 2*πy*). The most peculiar characteristic of Taylor's mean flow profile stood in its ability to satisfy the no-slip boundary condition at the sidewall despite its inviscid nature. This could be attributed to its wall-normal injection that disallowed any axial velocity contribution along the porous boundary.

Returning to the viscous flow problem in a porous channel, Yuan (1956) may have been the first to develop a solution for moderate to large Reynolds numbers and either suction or injection. His solution asymptotically reproduced Taylor's in the limit of a large injection 2 Will-be-set-by-IN-TECH

(a) Injection driven porous channel

Uw

(b) Injection driven porous tube

It is widely believed that Berman (1953) was among the earliest to examine the problem of laminar viscous flow bounded by porous surfaces (see Dauenhauer & Majdalani, 2003; Zhou & Majdalani, 2002). Although his first similarity transformation only applied to a planar configuration with wall suction, it has set forth the foundation for a number of follow-up investigations that relied on either analytical or numerical techniques to explore a variety of

Chronologically, these start with Sellars (1955) who extended Berman's solution to very large suction Reynolds numbers. He accomplished this by relaxing the no-slip boundary condition that became immaterial under this limiting condition. At the outset, he extracted a leading order approximation that corresponded to uniform axial motion, i.e. *F*(*y*) = *y*. Sellars integrated the ensuing equation based on his leading order approximation. His model thus uncovered the outer solution of this problem when viewed from a boundary layer perspective. Sellars' identification of a thin boundary layer at the wall for the large suction case would later

Of particular interest to this chapter is a classic article by Taylor (1956) in which he derived an inviscid rotational solution for both planar and axisymmetric channel flow configurations, in addition to cones and wedges. The absence of viscosity in his model led to approximations that were consistent with Berman's leading order solution of the Navier-Stokes equations

characteristic of Taylor's mean flow profile stood in its ability to satisfy the no-slip boundary condition at the sidewall despite its inviscid nature. This could be attributed to its wall-normal

Returning to the viscous flow problem in a porous channel, Yuan (1956) may have been the first to develop a solution for moderate to large Reynolds numbers and either suction or injection. His solution asymptotically reproduced Taylor's in the limit of a large injection

injection that disallowed any axial velocity contribution along the porous boundary.

2*πy*). The most peculiar

Fig. 1. Schematics of porous channels and tubes in which motion is sustained through

geometric configurations with either injection, suction, or both (Proudman, 1960).

*U*w

y x

z

prove crucial in subsequent developments of this problem.

expressed at large injection Reynolds numbers, i.e. *F*(*y*) = sin( 1

r

wall-normal injection.

Reynolds number. However, Yuan's model suffered from a singularity that appeared in the third derivative of the mean flow function *F*(*y*) taken at the centerline. This of course signaled the presence of a thin boundary layer that necessitated special treatment. The corresponding boundary layer would later be captured by Terrill (1965) who also described an insightful technique to solve this problem numerically.

In the interim, Berman (1958a) published his second work in which he extended the original planar problem to various geometric settings. This included the familiar case of a straight axisymmetric tube with permeable walls. Almost concurrently, White et al. (1958) advanced a series approximation to the porous channel problem for all ranges of the Reynolds number. However, White and co-workers employed a power series expansion that was centered around Re = 0. They also supplied a numerical solution to this problem. Despite the accuracy of their technique, their power series depended on two arbitrary constants that could only be determined numerically through a trial and error procedure. According to Terrill (1964), their method could be viewed as suitable for intermediate values of Re (15 ≤ Re ≤ 35). Otherwise, a transformation of the governing equation could be more effective at achieving direct numerical integration. Due to the penalty involved in evaluating the analytical constants of the attendant power series, this particular approach would be later abandoned. Nonetheless, it remained somewhat unique in its ability to provide a single analytical approximation that applied over the entire range of Re, a feat that standard perturbation methods failed to accomplish.

Along similar lines, Terrill (Terrill, 1964; 1965) compiled a comprehensive and detailed résumé of the perturbation solutions of this problem over all ranges of the Reynolds number. Therein, he derived and discussed several limiting cases such as Re = 0, |Re| � 1, Re → +∞, Re → −∞, and compared the various solutions with numerical simulations based on Runge–Kutta integration. For the numerical integration scheme, he introduced a transformation that would lead to a direct numerical solution with no need for predictor-corrector steps or shooting. On the flip side, his technique did not allow the pre-selection of the Reynolds number but rather the post-determination of Re at the conclusion of the numerical procedure. Before leaving this topic, we also note the work of Eckert et al. (1957) who, as far as the authors could verify, were the first to present a numerical solution for the laminar viscous motion in a porous channel.

As far as stability is concerned, the variety of analytical models considered for the planar case appeared to be both unique and stable (Terrill & Thomas, 1969; White, 2005). However, Robinson (1976) reported that dual solutions could exist for large suction while Zaturska et al. (1988) furnished a detailed stability analysis that rigorously showed that (at least) three types of solutions could co-exist. Even more intricate structures would arise in the case of axisymmetric flow in a porous tube. In this context, Terrill & Thomas (1969) have shown that, at least, dual solutions existed for the entire range of injection and suction Reynolds numbers while no steady solutions could be identified for 2.3 < Re < 9.1. At the time of this writing, the issue of stability of wall-injected flows remains an open area of investigation especially among applied mathematicians and fluid dynamicists.

#### **1.1 Relevance to propulsion systems**

In propulsive applications involving solid and hybrid rocket motors, modeling the mean flow proves to be important for a variety of reasons (Culick, 2006). The instantaneous flow field plays a key role in describing acoustic instability, particle-mean flow interactions, erosive burning, nozzle erosion, and thrust performance. The traditional modus operandi is to decompose the instantaneous motion into a steady average flow and an amalgam of unsteady

Terrill's treatment of the porous tube to include unsteady injection or suction at the sidewall. In the same vein, Erdogan & Imrak (2008) presented a laminar solution for the flow in a porous tube. Their solution was obtained by expanding the velocity field as a series of modified Bessel functions of order *n*. As for the problem involving wall regression, it was tackled by Dauenhauer & Majdalani (2003), Zhou & Majdalani (2002), and Majdalani & Zhou (2003) for the slab with regressing sidewall, and by Goto & Uchida (1990) and Majdalani et al. (2002) for the internal burning cylinder with expanding walls (see also Majdalani et al., 2009, for an

Internal Flows Driven by Wall-Normal Injection 99

The next noteworthy improvement in this area consists of the compressible Taylor–Culick profile that was first presented in multiple dimensions by Majdalani (2007b). His solution faithfully retained the essential ingredients of Culick's model, yet fully incorporated the effects of compressibility. This was accomplished through the use of a Rayleigh-Janzen expansion jointly with the vorticity-streamfunction approach for a compressible fluid. In asymptotic theory, the Rayleigh-Janzen expansion refers to a regular perturbation expansion in even powers of the Mach number that is ideally suited for the treatment of high speed flows (see Janzen, 1913; Rayleigh, 1916). A similar and equally impactful treatment of the planar configuration was subsequently presented by Maicke & Majdalani (2008) for the compressible Taylor flow analogue. Both analyses give rise to velocity fields that exhibit steep streamline curvatures that are consistent with numerical simulations of the compressible Navier-Stokes

As we move closer to the central topic of this chapter, we consider recent work in which the Taylor–Culick solution is reconstructed for the case of solid rocket motors with headwall injection or hybrid motors with a large headwall-to-sidewall velocity ratio (Majdalani, 2007a). The corresponding problem is analyzed in both axisymmetric and planar configurations by Majdalani & Saad (2007b) and Saad & Majdalani (2009b), respectively. This will be the topic of Section 2 where the solutions for the Taylor–Culick flow with arbitrary headwall injection are derived and compared to steady state, second order accurate inviscid computations. In subsequent work, Majdalani & Saad (2007a) and Saad & Majdalani (2010) manage to introduce a variational procedure based on Lagrangian multipliers to identify solutions of the Taylor–Culick type with varying kinetic energies. As it will be seen in Section 3, these will help to uncover a wide array of motions ranging from purely irrotational to highly rotational fields. The same approach is later applied to slab rocket motors (Saad & Majdalani, 2008a) and to swirl-driven cyclonic chambers with either single (Saad & Majdalani, 2008b) or multiple mantles (Saad & Majdalani, 2009a). In what follows, the main emphasis will be placed on the

In this section, we present a model for the mean flow in simulated solid or hybrid rocket motors with headwall injection. Our approach is based on a technique introduced by Majdalani (2007a) and Majdalani & Saad (2007b). The ability to account for arbitrary headwall injection will extend the Taylor-Culick approximation to a wider range of problems. For example, it will enable us to handle both solid and hybrid rocket motors in a unified analysis, the difference being in the relative magnitudes of the headwall-to-sidewall injection speeds. Our approach will be based on the vorticity-streamfunction formulation in which the vorticity transport equation will be used to obtain a functional relation between the streamfunction and the vorticity. The solution will then be retrieved from the vorticity equation. In the

motion driven by wall-normal injection in a porous, axisymmetric tube.

**2. Rotational models with headwall injection**

**2.1 Arbitrary injection**

error-free form).

equations.

wave contributions (Chedevergne et al., 2007; Culick, 2006; Majdalani, 2009). In this context, the mean flow represents the bulk motion of the gases and can be approximated by the steady-state solution for a porous tube or channel with wall-normal injection. As for the unsteady field, it refers to any perturbed disturbance that propagates within the chamber. Typical fluctuations are attributed to acoustic, vorticity, entropy, and hydrodynamic instability waves (Chu & Kovásznay, 1958). The importance of the mean flow is therefore evident due to the tight coupling between the steady and unsteady motions.

Although the earliest studies of solid rocket motor (SRM) stability treated the motors as porous enclosures, they failed to consider a suitable mean flow field. For example, the first theoretical study that explored the acoustic instability of rockets may be attributed to Grad (1949) (see Culick, 2006, for greater detail). However, Grad assumed that the mean flow could be ignored as in the case of a stagnant medium, thus limiting his analysis to that of aeroacoustic instability in a cylindrical chamber with no mean flow motion.

Nearly a decade later, the work of McClure and coworkers would prove instrumental in the understanding of rocket motor stability, especially in the development of the energy balance framework. However, principal efforts in this direction have focused on the thin region near the injecting surface (Hart & McClure, 1965; Hart et al., 1960; Hart & Cantrell, 1963; Hart & McClure, 1959; McClure et al., 1960). In fact, McClure et al. (1963) may have been the first to employ a mean flow approximation in their analysis of the aeroacoustic field in SRMs. Their model of choice corresponded to the irrotational motion of an ideal gas in a porous cylinder or between two parallel porous plates. It hence constituted a substantial improvement over the stagnation flow model and, for the first time, succeeded in identifying the intimate coupling between the mean flow and the unsteady wave motion.

It was not until Culick (1966) that a robust representation of the mean flow in circular port motors would be introduced. Despite its inviscid nature, Culick's model was rotational and could satisfy the no-slip requirement at the sidewall. The profile itself coincided with that obtained by Taylor (1956) a decade earlier, albeit in an entirely different application (i.e. paper manufacturing). Culick (1966) derived his solution in the context of a propulsive application that quickly proved to be quintessential to several combustion instability studies, particle-mean-flow interactions, turbulence characterization, and other related investigations of solid propellant rocket motors. It is usually referred to as the Taylor–Culick profile and remains one of the most cited models in rocket motor analysis. For example, Chedevergne et al. (2006), Abu-Irshaid et al. (2007), Griffond et al. (2000), Beddini (1986) and Flandro & Majdalani (2003) made extensive use of the Taylor–Culick model as a basis for their instability work.

#### **1.2 Beyond Culick's solution**

Going beyond the Taylor-Culick solution, Majdalani and coworkers have explored a variety of avenues that extended the classic model by providing higher order approximations that could take into account additional factors that are omitted in the inviscid formulation. These include the effects of viscosity, grain taper, wall regression, compressibility, and headwall injection. For example, the sensitivity of the mean flow to viscosity is discussed by Majdalani & Akiki (2010) whereas the effects of tapering of internal bores are addressed by Saad et al. (2006) for the rectangular port slab geometry and by Sams et al. (2007) for the internal burning cylinder with circular cross-section. Other improvements include the work of Kurdyumov (2006) who extended the Taylor–Culick solution to chambers with irregular cross-sections, such as those with a star-shaped perforation. Furthermore, Tsangaris et al. (2007) generalized 4 Will-be-set-by-IN-TECH

wave contributions (Chedevergne et al., 2007; Culick, 2006; Majdalani, 2009). In this context, the mean flow represents the bulk motion of the gases and can be approximated by the steady-state solution for a porous tube or channel with wall-normal injection. As for the unsteady field, it refers to any perturbed disturbance that propagates within the chamber. Typical fluctuations are attributed to acoustic, vorticity, entropy, and hydrodynamic instability waves (Chu & Kovásznay, 1958). The importance of the mean flow is therefore evident due to

Although the earliest studies of solid rocket motor (SRM) stability treated the motors as porous enclosures, they failed to consider a suitable mean flow field. For example, the first theoretical study that explored the acoustic instability of rockets may be attributed to Grad (1949) (see Culick, 2006, for greater detail). However, Grad assumed that the mean flow could be ignored as in the case of a stagnant medium, thus limiting his analysis to that of

Nearly a decade later, the work of McClure and coworkers would prove instrumental in the understanding of rocket motor stability, especially in the development of the energy balance framework. However, principal efforts in this direction have focused on the thin region near the injecting surface (Hart & McClure, 1965; Hart et al., 1960; Hart & Cantrell, 1963; Hart & McClure, 1959; McClure et al., 1960). In fact, McClure et al. (1963) may have been the first to employ a mean flow approximation in their analysis of the aeroacoustic field in SRMs. Their model of choice corresponded to the irrotational motion of an ideal gas in a porous cylinder or between two parallel porous plates. It hence constituted a substantial improvement over the stagnation flow model and, for the first time, succeeded in identifying the intimate coupling

It was not until Culick (1966) that a robust representation of the mean flow in circular port motors would be introduced. Despite its inviscid nature, Culick's model was rotational and could satisfy the no-slip requirement at the sidewall. The profile itself coincided with that obtained by Taylor (1956) a decade earlier, albeit in an entirely different application (i.e. paper manufacturing). Culick (1966) derived his solution in the context of a propulsive application that quickly proved to be quintessential to several combustion instability studies, particle-mean-flow interactions, turbulence characterization, and other related investigations of solid propellant rocket motors. It is usually referred to as the Taylor–Culick profile and remains one of the most cited models in rocket motor analysis. For example, Chedevergne et al. (2006), Abu-Irshaid et al. (2007), Griffond et al. (2000), Beddini (1986) and Flandro & Majdalani (2003) made extensive use of the Taylor–Culick model as a basis for their instability

Going beyond the Taylor-Culick solution, Majdalani and coworkers have explored a variety of avenues that extended the classic model by providing higher order approximations that could take into account additional factors that are omitted in the inviscid formulation. These include the effects of viscosity, grain taper, wall regression, compressibility, and headwall injection. For example, the sensitivity of the mean flow to viscosity is discussed by Majdalani & Akiki (2010) whereas the effects of tapering of internal bores are addressed by Saad et al. (2006) for the rectangular port slab geometry and by Sams et al. (2007) for the internal burning cylinder with circular cross-section. Other improvements include the work of Kurdyumov (2006) who extended the Taylor–Culick solution to chambers with irregular cross-sections, such as those with a star-shaped perforation. Furthermore, Tsangaris et al. (2007) generalized

the tight coupling between the steady and unsteady motions.

between the mean flow and the unsteady wave motion.

work.

**1.2 Beyond Culick's solution**

aeroacoustic instability in a cylindrical chamber with no mean flow motion.

Terrill's treatment of the porous tube to include unsteady injection or suction at the sidewall. In the same vein, Erdogan & Imrak (2008) presented a laminar solution for the flow in a porous tube. Their solution was obtained by expanding the velocity field as a series of modified Bessel functions of order *n*. As for the problem involving wall regression, it was tackled by Dauenhauer & Majdalani (2003), Zhou & Majdalani (2002), and Majdalani & Zhou (2003) for the slab with regressing sidewall, and by Goto & Uchida (1990) and Majdalani et al. (2002) for the internal burning cylinder with expanding walls (see also Majdalani et al., 2009, for an error-free form).

The next noteworthy improvement in this area consists of the compressible Taylor–Culick profile that was first presented in multiple dimensions by Majdalani (2007b). His solution faithfully retained the essential ingredients of Culick's model, yet fully incorporated the effects of compressibility. This was accomplished through the use of a Rayleigh-Janzen expansion jointly with the vorticity-streamfunction approach for a compressible fluid. In asymptotic theory, the Rayleigh-Janzen expansion refers to a regular perturbation expansion in even powers of the Mach number that is ideally suited for the treatment of high speed flows (see Janzen, 1913; Rayleigh, 1916). A similar and equally impactful treatment of the planar configuration was subsequently presented by Maicke & Majdalani (2008) for the compressible Taylor flow analogue. Both analyses give rise to velocity fields that exhibit steep streamline curvatures that are consistent with numerical simulations of the compressible Navier-Stokes equations.

As we move closer to the central topic of this chapter, we consider recent work in which the Taylor–Culick solution is reconstructed for the case of solid rocket motors with headwall injection or hybrid motors with a large headwall-to-sidewall velocity ratio (Majdalani, 2007a). The corresponding problem is analyzed in both axisymmetric and planar configurations by Majdalani & Saad (2007b) and Saad & Majdalani (2009b), respectively. This will be the topic of Section 2 where the solutions for the Taylor–Culick flow with arbitrary headwall injection are derived and compared to steady state, second order accurate inviscid computations. In subsequent work, Majdalani & Saad (2007a) and Saad & Majdalani (2010) manage to introduce a variational procedure based on Lagrangian multipliers to identify solutions of the Taylor–Culick type with varying kinetic energies. As it will be seen in Section 3, these will help to uncover a wide array of motions ranging from purely irrotational to highly rotational fields. The same approach is later applied to slab rocket motors (Saad & Majdalani, 2008a) and to swirl-driven cyclonic chambers with either single (Saad & Majdalani, 2008b) or multiple mantles (Saad & Majdalani, 2009a). In what follows, the main emphasis will be placed on the motion driven by wall-normal injection in a porous, axisymmetric tube.

#### **2. Rotational models with headwall injection**

#### **2.1 Arbitrary injection**

In this section, we present a model for the mean flow in simulated solid or hybrid rocket motors with headwall injection. Our approach is based on a technique introduced by Majdalani (2007a) and Majdalani & Saad (2007b). The ability to account for arbitrary headwall injection will extend the Taylor-Culick approximation to a wider range of problems. For example, it will enable us to handle both solid and hybrid rocket motors in a unified analysis, the difference being in the relative magnitudes of the headwall-to-sidewall injection speeds. Our approach will be based on the vorticity-streamfunction formulation in which the vorticity transport equation will be used to obtain a functional relation between the streamfunction and the vorticity. The solution will then be retrieved from the vorticity equation. In the

*r* = *r*∗ *a*

*u* = *u*∗ *U*w

to all subsequent developments.

**2.2.2 Euler-based formulation**

or, in vector form

where

equation for steady, inviscid motion

⎧ ⎪⎪⎪⎨

⎪⎪⎪⎩

turbulent-like behavior.

; *z* = *z*∗ *a*

> ; *w* = *w*∗ *U*w

the outset, the normalized Euler equations with no swirl can be written as 1 *r*

> *u ∂u ∂r* + *w ∂u ∂z* = − *∂p*

*u ∂w ∂r* + *w ∂w ∂z* = − *∂p*

One may now invoke the dyadic vector identity **u** · ∇**u** ≡ ∇( 1

Finally, four boundary conditions can be prescribed by writing

*∂*(*ru*) *∂r* + *∂w*

by taking the curl of the resulting expression into (4b), one obtains the vorticity transport

*u*(0, *z*) = 0 no flow across centerline *w*(1, *z*) = 0 no slip at sidewall

*w*(*r*, 0) = *w*0(*r*) axial inflow at headwall

⎧ ⎪⎨

⎪⎩

where the headwall injection profile may take any of the following plausible forms

*w*0(*r*) =

*u*(1, *z*) = −1 constant radial inflow at sidewall

*W*c = const *W*c cos( 1

*W*c(1 − *rm*)

Here *m* is the power-law exponent that may be taken as 2 for laminar and 7 or 8 for

2*πr*2)

; ∇ = *a*∇∗; *p* = *p*∗

Internal Flows Driven by Wall-Normal Injection 101

; **Ω** = **Ω**∗*a U*w

where starred variables denote dimensional quantities. Note that this normalization applies

A non-reactive motion may be assumed, prompted by the thin reactive zone above the grain surface. Following Culick (1966), the flow can be taken to be steady, inviscid, incompressible, rotational, and axisymmetric. It should be noted that Majdalani (2007b) and Maicke & Majdalani (2008) have provided compressible Taylor–Culick solutions under isentropic flow conditions. These confirm the suitability of the present model for a variety of applications in which the effects of compressibility are small. Chu et al. (2003) and Vyas et al. (2003) have also demonstrated that the flow field above the thin flame zone may be treated as non-reactive. At

*ρU*2 w

; *W*c = *W*∗

; *ψ* = *ψ*∗ *a*2*U*w ;

> ; *L* = *L*∗ *a*

*∂z* = 0 (3a)

∇ · **u** = 0 (4a) **u** · ∇**u** = −∇*p* (4b)

∇ × (**u** × **Ω**) = 0 (5)

**Ω** = ∇ × **u** (6)

*∂r* (3b)

*∂z* (3c)

2**u** · **u**) − **u** ×∇× **u**. Then,

(7)

(8)

(2)

c *U*w

Fig. 2. Schematic of an idealized solid rocket motor with sidewall injection.

process, a multitude of injection profiles will be extracted using superposition. Despite the nonlinearity of the vorticity transport equation near the headwall, it will be shown that the solution becomes progressively more linear in the downstream direction, a factor that permits the use of superposition. Incidentally, the linearity of the vorticity-streamfunction relation used in these studies has been shown by Kurdyumov (2008) to hold true away from the headwall. Finally, the resulting approximations will be tested using three representative injection profiles for which comparisons with finite volume CFD simulations of the Euler equations will be performed.

#### **2.2 Mathematical idealization**

A rocket motor can be idealized as a cylindrical chamber of porous length *L*∗ and radius *a* with both a reactive headwall and a nozzleless aft end as shown in Figure 2. The radial and axial velocities are represented by *u*∗ and *w*∗, respectively, while *r*∗ and *z*∗ stand for the radial and axial coordinates used to describe the solution from the headwall to the typical nozzle attachment point at the chamber outlet. At the headwall, a fluid stream (which may denote an oxidizer or gaseous propellant mixture) is injected into the chamber at a prescribed velocity *w*∗ 0 (*r*∗). This could be given by

$$w\_0^\*(r^\*) = w^\*(r^\*, z^\* = 0) = \begin{cases} W\_\emptyset^\* = \text{const} & \text{uniform} \\ W\_\emptyset^\* \cos(\frac{1}{2}\pi r^{\*2}/a^2) & \text{cosine} \\ W\_\emptyset^\* [1 - (r^\*/a)^m] & \text{laminar and turbulent} \\ W\_\emptyset^\* (1 - r^\*/a)^{1/m} & \text{turbulent} \end{cases} \tag{1}$$

where *W*∗ c = *w*∗(0, 0) is the centerline speed at the headwall (a constant), *m* is some integer, and the asterisk denotes a dimensional variable. The incoming stream merges with the cross flow generated by uniform mass addition along the porous sidewall. Naturally, the sidewall injection velocity *U*w = −*u*∗(*a*, *z*∗) is commensurate with propellant or fuel regression rates. In hybrids, *U*w can be appreciably smaller than *W*∗ c due to slow fuel pyrolysis; in SRM analysis, these two values can be identical.

#### **2.2.1 Normalization**

It is useful to normalize all recurring variables and operators. This can be done by following Majdalani & Saad (2007b) and setting

6 Will-be-set-by-IN-TECH

z r a

\* \*

Wc

equations will be performed.

**2.2 Mathematical idealization**

0 (*r*∗). This could be given by

0 (*r*∗) = *w*∗(*r*∗, *z*∗ = 0) =

In hybrids, *U*w can be appreciably smaller than *W*∗

analysis, these two values can be identical.

Majdalani & Saad (2007b) and setting

*w*∗

**2.2.1 Normalization**

*w*∗

where *W*∗

\*

L

process, a multitude of injection profiles will be extracted using superposition. Despite the nonlinearity of the vorticity transport equation near the headwall, it will be shown that the solution becomes progressively more linear in the downstream direction, a factor that permits the use of superposition. Incidentally, the linearity of the vorticity-streamfunction relation used in these studies has been shown by Kurdyumov (2008) to hold true away from the headwall. Finally, the resulting approximations will be tested using three representative injection profiles for which comparisons with finite volume CFD simulations of the Euler

A rocket motor can be idealized as a cylindrical chamber of porous length *L*∗ and radius *a* with both a reactive headwall and a nozzleless aft end as shown in Figure 2. The radial and axial velocities are represented by *u*∗ and *w*∗, respectively, while *r*∗ and *z*∗ stand for the radial and axial coordinates used to describe the solution from the headwall to the typical nozzle attachment point at the chamber outlet. At the headwall, a fluid stream (which may denote an oxidizer or gaseous propellant mixture) is injected into the chamber at a prescribed velocity

c = const uniform

c (1 − *r*∗/*a*)1/*m* turbulent

c = *w*∗(0, 0) is the centerline speed at the headwall (a constant), *m* is some integer,

and the asterisk denotes a dimensional variable. The incoming stream merges with the cross flow generated by uniform mass addition along the porous sidewall. Naturally, the sidewall injection velocity *U*w = −*u*∗(*a*, *z*∗) is commensurate with propellant or fuel regression rates.

It is useful to normalize all recurring variables and operators. This can be done by following

2*πr*∗2/*a*2) cosine

c [1 − (*r*∗/*a*)*m*] laminar and turbulent

c due to slow fuel pyrolysis; in SRM

(1)

Fig. 2. Schematic of an idealized solid rocket motor with sidewall injection.

⎧ ⎪⎪⎪⎨

*W*∗

*W*∗ c cos( 1

*W*∗

*W*∗

⎪⎪⎪⎩

\*

Uw

$$\begin{aligned} r &= \frac{r^\*}{a}; \; z = \frac{z^\*}{a}; \; \nabla = a\nabla^\*; \; p = \frac{p^\*}{\rho U\_\text{W}^2}; \; \psi = \frac{\psi^\*}{a^2 U\_\text{W}};\\ u &= \frac{u^\*}{U\_\text{W}}; \; w = \frac{w^\*}{U\_\text{W}}; \; \Omega = \frac{\Omega^\* a}{U\_\text{W}}; \; W\_\text{c} = \frac{W\_\text{c}^\*}{U\_\text{W}}; \; L = \frac{L^\*}{a} \end{aligned} \tag{2}$$

where starred variables denote dimensional quantities. Note that this normalization applies to all subsequent developments.

#### **2.2.2 Euler-based formulation**

A non-reactive motion may be assumed, prompted by the thin reactive zone above the grain surface. Following Culick (1966), the flow can be taken to be steady, inviscid, incompressible, rotational, and axisymmetric. It should be noted that Majdalani (2007b) and Maicke & Majdalani (2008) have provided compressible Taylor–Culick solutions under isentropic flow conditions. These confirm the suitability of the present model for a variety of applications in which the effects of compressibility are small. Chu et al. (2003) and Vyas et al. (2003) have also demonstrated that the flow field above the thin flame zone may be treated as non-reactive. At the outset, the normalized Euler equations with no swirl can be written as

$$\frac{1}{r}\frac{\partial (ru)}{\partial r} + \frac{\partial w}{\partial z} = 0 \tag{3a}$$

$$u\frac{\partial u}{\partial r} + w\frac{\partial u}{\partial z} = -\frac{\partial p}{\partial r} \tag{3b}$$

$$u\frac{\partial w}{\partial r} + w\frac{\partial w}{\partial z} = -\frac{\partial p}{\partial z} \tag{3c}$$

or, in vector form

$$
\underbrace{\nabla \cdot \mathbf{u}}\_{\cdot} = \underbrace{\mathbf{0}}\_{\cdot} \tag{4a}
$$

$$\mathbf{u} \cdot \nabla \mathbf{u} = -\nabla p \tag{4b}$$
One may now involve the dyadic vector identity  $\mathbf{u} \cdot \nabla \mathbf{u} \equiv \nabla(\frac{1}{2}\mathbf{u} \cdot \mathbf{u}) - \mathbf{u} \times \nabla \times \mathbf{u}$ . Then, by taking the curl of the resulting expression into (4b), one obtains the vorticity transport equation for steady, invisoid motion

$$\nabla \times (\mathbf{u} \times \mathbf{D}) = 0 \tag{5}$$

where

$$
\mathbf{\Omega} = \nabla \times \mathbf{u} \tag{6}
$$

Finally, four boundary conditions can be prescribed by writing

$$\begin{cases} u(0,z) = 0 & \text{no flow across centerline} \\ w(1,z) = 0 & \text{no slip at sidewall} \\ u(1,z) = -1 & \text{constant radial inflow at sidewall} \\ w(r,0) = w\_0(r) & \text{axial inflow at headwall} \end{cases} \tag{7}$$

where the headwall injection profile may take any of the following plausible forms

$$w\_0(r) = \begin{cases} W\_\text{c} = \text{const} \\ W\_\text{c} \cos(\frac{1}{2}\pi r^2) \\ W\_\text{c} (1 - r^m) \end{cases} \tag{8}$$

Here *m* is the power-law exponent that may be taken as 2 for laminar and 7 or 8 for turbulent-like behavior.

**2.4 Solution by eigenfunction expansion**

order in which they appear. Starting with (15a), we obtain:

*∂z* = *α*¯ *A* cos( 1

or *A* = 0. Without loss of generality, we set *B* = 1 and rewrite (14b) as

*∂r* = *rC*(*α*¯ *z* + *β*¯) cos( 1

error term in (5) that will be examined in Section 2.6. In the interim, we take

*ψ*(*r*, *z*) =

*αn*. One such case corresponds to Taylor's family of solutions for which

Accordingly, by setting *βn* = 0, we recover Culick's original solution

*ψ*(*r*, *z*) = *z* sin( 1

∞ ∑ *n*=0

At this juncture, we apply the sidewall injection condition (14c) to produce

*αn* sin[(*n* + 1

*∂ψ*(0, *z*)

*∂ ψ*(1, *z*)

For convenience, we introduce *χn* ≡ 1

*∂ψ*(1, *z*) *∂z* =

> 1 *r*

*βn* 1 0

*∂ψ*(*r*, 0) *∂r* = *π*

2*C*) = 0. This is satisfied by

and so cos( 1

compacted into

with

The application of the boundary conditions must be carefully carried out, preferably in the

Internal Flows Driven by Wall-Normal Injection 103

2*Cr*2) + *α*¯ *B* sin( 1

2*Cr*2)

2 )*π*] = 1 or

2*πr*

∞ ∑ *n*=0

This keystone equality encapsulates several possible outcomes depending on the behavior of

*ψ*(*r*, *z*) = *z* sin( 1

Other forms of *αn* will be discussed in Section 3. At present, we let *α*0 = 1 and reduce (21) into

2*πr* 2) +

Lastly, the headwall condition (14d) may be fulfilled through the use of orthogonality. Starting

*π* 1 0

∞ ∑ *n*=0

one can take advantage of the orthogonality of the cosine function to secure

(2*n* + 1) cos2 *χn r dr* = 1

Using *Cn* = (2*n* + 1)*π*, we obtain an infinite series solution to (13). This process introduces an

*ψn*(*r*, *z*)=(*αnz* + *βn*) sin[(*n* + 1

∞ ∑ *n*=0 2*Cr*2)

*r*=1 = 0; ∀ *z* ∈ **R**+

*C* = *Cn* = (2*n* + 1)*π*; ∀ *n* ∈ **N**0 (19)

2 )*πr*

∞ ∑ *n*=0

*α*0 = 1 and *αn* = 0; ∀*n* �= 0 (23)

2 (2*n* + 1)*πr*2 so that the total streamfunction may be

(*αnz* + *βn*) sin *χn* (21)

*r*=0 = 0 (17)

2] (20)

(−1)*nαn* = 1 (22)

2) (24)

*βn* sin *χn* (25)

*w*0(*r*) cos *χn r dr* (27)

(2*n* + 1)*βn* cos *χn* = *w*0(*r*) (26)

0 (18)

#### **2.3 Vorticity-streamfunction formulation**

Continuity is fulfilled by the Stokes streamfunction in cylindrical coordinates when written as

$$u = -\frac{1}{r} \frac{\partial \psi}{\partial z}; \quad w = \frac{1}{r} \frac{\partial \psi}{\partial r} \tag{9}$$

Having a single nonzero component in the azimuthal direction, the vorticity reduces to

$$
\Omega \mathbf{0} = \Omega\_{\theta} \mathbf{e}\_{\theta} \equiv \Omega \mathbf{e}\_{\theta} \tag{10}
$$

Its substitution into the vorticity transport equation (5) yields

$$\frac{\partial \psi}{\partial r} \frac{\partial}{\partial z} \left( \frac{\Omega}{r} \right) - \frac{\partial \psi}{\partial z} \frac{\partial}{\partial r} \left( \frac{\Omega}{r} \right) = 0 \quad \text{or} \quad \frac{(\Omega/r)\_z}{(\Omega/r)\_r} = \frac{\psi\_z}{\psi\_r} \tag{11}$$

where the subscripts denote differentiation with respect to *r* or *z*, respectively. Equation (11) may be satisfied by taking Ω = *rF*(*ψ*) since

$$\frac{(\Omega/r)\_z}{(\Omega/r)\_r} = \frac{[F(\boldsymbol{\psi})]\_z}{[F(\boldsymbol{\psi})]\_r} = \frac{F\_{\boldsymbol{\Psi}}\boldsymbol{\psi}\_z}{F\_{\boldsymbol{\Psi}}\boldsymbol{\psi}\_r} = \frac{\boldsymbol{\psi}\_z}{\boldsymbol{\psi}\_r} \tag{12}$$

So we follow Culick (1966) and set Ω = *C*2*rψ*. Despite the non-uniqueness of this relation, it enables us to secure (5). At this point, straightforward substitution into the vorticity equation (6) renders immediately the second-order PDE associated with the Taylor–Culick problem,

$$\frac{\partial^2 \psi}{\partial z^2} + \frac{\partial^2 \psi}{\partial r^2} - \frac{1}{r} \frac{\partial \psi}{\partial r} + \mathbb{C}^2 r^2 \psi = 0 \tag{13}$$

with the particular set of constraints,

$$\lim\_{r \to 0} \frac{1}{r} \frac{\partial \psi(r, z)}{\partial z} = 0 \tag{14a}$$

$$\frac{\partial \,\psi(1,z)}{\partial r} = 0\tag{14b}$$

$$\frac{\partial \psi(1, z)}{\partial z} = 1 \tag{14c}$$

$$\frac{1}{r}\frac{\partial \psi(r,0)}{\partial r} = w\_0(r) \tag{14d}$$

By virtue of L'Hôpital's rule, removing the singularity in (14a) requires that both

$$\frac{\partial \psi(0, z)}{\partial z} = 0 \tag{15a}$$

$$\frac{\partial^2 \psi(0, z)}{\partial r \partial z} = 0 \tag{15b}$$

Being linear, (13) is solvable by separation of variables; it yields

$$\psi(r,z) = (\overline{a}z + \overline{\beta})[A\cos(\frac{1}{2}Cr^2) + B\sin(\frac{1}{2}Cr^2)]\tag{16}$$

This expression satisfies (15b) identically. Henceforth, (14a) may be superseded by (15a). We then proceed to implement the problem's constraints so that a solution may be realized.

8 Will-be-set-by-IN-TECH

Continuity is fulfilled by the Stokes streamfunction in cylindrical coordinates when written as

Having a single nonzero component in the azimuthal direction, the vorticity reduces to

 Ω *r*

= [*F*(*ψ*)]*z* [*F*(*ψ*)]*r*

*∂*2*ψ ∂r*2 − 1 *r ∂ψ ∂r*

lim *r*→0 1 *r*

where the subscripts denote differentiation with respect to *r* or *z*, respectively. Equation (11)

So we follow Culick (1966) and set Ω = *C*2*rψ*. Despite the non-uniqueness of this relation, it enables us to secure (5). At this point, straightforward substitution into the vorticity equation (6) renders immediately the second-order PDE associated with the Taylor–Culick problem,

*∂ψ*(*r*, *z*)

*∂ ψ*(1, *z*)

*∂ψ*(1, *z*)

*∂ ψ*(*r*, 0)

*∂ψ*(0, *z*)

*∂*2 *ψ*(0, *z*)

This expression satisfies (15b) identically. Henceforth, (14a) may be superseded by (15a). We then proceed to implement the problem's constraints so that a solution may be realized.

1 *r*

By virtue of L'Hôpital's rule, removing the singularity in (14a) requires that both

Being linear, (13) is solvable by separation of variables; it yields

*ψ*(*r*, *z*)=(*α*¯ *z* + *β*¯)[*A* cos( 1

; *w* = 1 *r ∂ψ*

*∂r* (9)

(11)

(12)

**Ω** = Ω*θ***e***θ* ≡ Ω**e***θ* (10)

= *ψz ψr*

+ *C*2*r*2*ψ* = 0 (13)

*∂z* = 0 (14a)

*∂r* = 0 (14b)

*∂z* = 1 (14c)

*∂r* = *w*0(*r*) (14d)

*∂z* = 0 (15a)

*∂r∂z* = 0 (15b)

2*Cr*2)] (16)

2*Cr*2) + *B* sin( 1

(Ω/*r*)*r*

= *ψz ψr*

= 0 or (Ω/*r*)*z*

= *Fψψz Fψψr*

*u* = −1 *r ∂ψ ∂z*

Its substitution into the vorticity transport equation (5) yields

(Ω/*r*)*z* (Ω/*r*)*r*

*∂*2*ψ ∂z*2 +

 Ω *r* − *∂ψ ∂z ∂ ∂r*

**2.3 Vorticity-streamfunction formulation**

*∂ψ ∂r ∂ ∂z*

may be satisfied by taking Ω = *rF*(*ψ*) since

with the particular set of constraints,

#### **2.4 Solution by eigenfunction expansion**

The application of the boundary conditions must be carefully carried out, preferably in the order in which they appear. Starting with (15a), we obtain:

$$\frac{\partial \psi(0, z)}{\partial z} = \bar{a}A \cos(\frac{1}{2}Cr^2) + \bar{a}B \sin(\frac{1}{2}Cr^2) \Big|\_{r=0} = 0 \tag{17}$$

or *A* = 0. Without loss of generality, we set *B* = 1 and rewrite (14b) as

$$\frac{\partial \,\psi(1,z)}{\partial r} = r\mathbb{C}(\overline{a}z + \overline{\beta}) \, \cos(\frac{1}{2}\mathbb{C}r^2) \Big|\_{r=1} = 0; \quad \forall z \in \mathbb{R}\_0^+ \tag{18}$$

and so cos( 1 2*C*) = 0. This is satisfied by

$$\mathbb{C} = \mathbb{C}\_n = (2n+1)\pi; \forall n \in \mathbb{N}\_0 \tag{19}$$

Using *Cn* = (2*n* + 1)*π*, we obtain an infinite series solution to (13). This process introduces an error term in (5) that will be examined in Section 2.6. In the interim, we take

$$
\psi\_{\rm ll}(r, z) = (\alpha\_{\rm ll} z + \beta\_{\rm ll}) \sin[(n + \frac{1}{2})\pi r^2] \tag{20}
$$

For convenience, we introduce *χn* ≡ 1 2 (2*n* + 1)*πr*2 so that the total streamfunction may be compacted into

$$\psi(r,z) = \sum\_{n=0}^{\infty} (\alpha\_n z + \beta\_n) \sin \chi\_n \tag{21}$$

At this juncture, we apply the sidewall injection condition (14c) to produce

$$\frac{\partial \psi(1, z)}{\partial z} = \sum\_{n=0}^{\infty} a\_n \sin[(n + \frac{1}{2})\pi] = 1 \quad \text{or} \quad \sum\_{n=0}^{\infty} (-1)^n a\_n = 1 \tag{22}$$

This keystone equality encapsulates several possible outcomes depending on the behavior of *αn*. One such case corresponds to Taylor's family of solutions for which

$$
\mathfrak{a}\_0 = 1 \quad \text{and} \quad \mathfrak{a}\_n = 0; \quad \forall n \neq 0 \tag{23}
$$

Accordingly, by setting *βn* = 0, we recover Culick's original solution

$$
\psi(r, z) = z \sin(\frac{1}{2}\pi r^2) \tag{24}
$$

Other forms of *αn* will be discussed in Section 3. At present, we let *α*0 = 1 and reduce (21) into

$$\psi(r,z) = z\sin(\frac{1}{2}\pi r^2) + \sum\_{n=0}^{\infty} \beta\_n \sin \chi\_n \tag{25}$$

Lastly, the headwall condition (14d) may be fulfilled through the use of orthogonality. Starting with

$$\frac{1}{r}\frac{\partial\psi(r,0)}{\partial r} = \pi\sum\_{n=0}^{\infty}(2n+1)\beta\_n\cos\chi\_n = w\_0(r) \tag{26}$$

one can take advantage of the orthogonality of the cosine function to secure

$$\beta\_{\rm li} \int\_{0}^{1} (2n+1) \cos^{2} \chi\_{\rm li} r \, dr = \frac{1}{\pi} \int\_{0}^{1} w\_{0}(r) \cos \chi\_{\rm li} r \, dr \tag{27}$$

These are prescribed by classic profiles used by Berman (1953) (half cosine), Poiseuille (White,

Internal Flows Driven by Wall-Normal Injection 105

*βn* = 4(−1)*nW*c

2) + 4*W*c *π*2

> 2) + 4*W*c *π*

The character of (34) is illustrated in Figure 4. Using *W*c = *U*w = 1, a balance between sidewall and headwall injection causes the streamline originating at the corner (*r* = 1, *z* = 0) to bisect the flow field at an angle of *π*/4 as shown in Figure 4(b). By concentrating on a thin region near the sidewall in Figure 4(c), it may be seen that the solution conforms to the stated boundary conditions. It is also evident that *w*0(*r*) = *W*c = 1 corresponds to a simulated solid

> *π* ≡ *W*h; *n* = 0 0; otherwise

> > 2*πr*

2*πr* 2);

> 2*πr* 2)

Ω(*r*, *z*) = *π*2*rz* sin( 1

propellant grain that is burning evenly along its headwall and sidewall boundaries.

∞ ∑ *n*=0

> ∞ ∑ *n*=0

> > 2*πr*

(−1)*n*

(−1)*n*

*π*2(2*n* + 1)2 (33)

(2*n* + 1)2 sin *χn* (34)

(2*n* + 1) cos *χn* (35)

2) (36)

2) (38)

(37)

(39)

2005), and others (uniform flow).

In this case, the headwall injection sequence *βn* collapses into

*ψ*(*r*, *z*) = *z* sin( 1

The axial velocity and vorticity may be easily determined to be

*w*(*r*, *z*) = *πz* cos( 1

**2.5.2 Similarity-conforming cosine injection**

Using (9), the streamfunction becomes

Their axial velocity and vorticity correspond to

This behavior will be discussed in Section 2.6.

For the cosine injection profile, we use (28) to obtain

*βn* =

⎧ ⎨ ⎩ *W*c

*ψ*(*r*, *z*)=(*z* + *W*h) cos( 1

*w*(*r*, *z*) = *π*(*z* + *W*h) cos( 1

Ω(*r*, *z*) = *π*2*r*(*z* + *W*h) sin( 1

It should be noted that while the solutions derived for most injection profiles are approximate, the one corresponding to the similarity-conforming Berman injection will prove to be exact.

The streamlines associated with the cosine headwall injection case are depicted in Figure 5.

2*πr*

2*πr*

**2.5.1 Uniform injection**

whence

Fig. 3. Streamline patterns corresponding to the classic Taylor–Culick profile with no headwall injection.

or

$$\beta\_{\rm li} = \frac{4}{(2n+1)\pi} \int\_0^1 w\_0(r) \cos \chi\_{\rm li} r \, dr \tag{28}$$

With *βn* in hand, the streamfunction is fully determined, namely,

$$\psi(r,z) = z\sin(\frac{1}{2}\pi r^2) + \sum\_{n=0}^{\infty} \left[\frac{4}{(2n+1)\pi} \int\_0^1 w\_0(r)\cos\chi\_n r \, dr\right] \sin\chi\_n\tag{29}$$

The radial and axial velocities follow and these may be expressed as

$$u(r) = -r^{-1} \sin(\frac{1}{2}\pi r^2);$$

$$w(r, z) = \pi z \cos(\frac{1}{2}\pi r^2) + \pi \sum\_{n=0}^{\infty} (2n+1)\beta\_n \cos \chi\_n. \tag{30}$$

Interestingly, the radial velocity remains independent of the headwall injection sequence, *βn*. Finally, the vorticity may be deduced from

$$\Omega(r,z) = \pi^2 rz \sin(\frac{1}{2}\pi r^2) + \pi^2 r \sum\_{n=0}^{\infty} (2n+1)^2 \beta\_n \sin \chi\_n \tag{31}$$

This extended form of the Taylor–Culick profile represents a solution for an arbitrary headwall injection pattern *w*0(*r*) that may be prescribed by the proper specification of *βn* through (28). By way of confirmation, the classical Taylor–Culick solution with inert headwall may be readily recovered by setting *βn* = 0 everywhere. The streamline patterns associated with this historical benchmark are illustrated in Figure 3.

#### **2.5 Axisymmetric headwall injection profiles**

The framework may be tested using a variable headwall injection profile. To be consistent with the underlying flow assumptions, we employ an axisymmetric function to specify the injection pattern at *z* = 0, namely,

$$w\_0(r) = \begin{cases} W\_\mathbb{C} = \text{const} & \text{uniform} \\\\ \mathcal{W}\_\mathbb{C} \cos(\frac{1}{2}\pi r^2) & \text{half cosine} \\\\ \mathcal{W}\_\mathbb{C} (1 - r^2) & \text{parabolic} \end{cases} \tag{32}$$

These are prescribed by classic profiles used by Berman (1953) (half cosine), Poiseuille (White, 2005), and others (uniform flow).

#### **2.5.1 Uniform injection**

In this case, the headwall injection sequence *βn* collapses into

$$\beta\_{\rm ll} = \frac{4(-1)^n W\_{\rm c}}{\pi^2 (2n+1)^2} \tag{33}$$

whence

10 Will-be-set-by-IN-TECH

0.2

0 0.2L 0.4L 0.6L 0.8L L

z

� 1 0

� 4 (2*n* + 1)*π*

2*πr* 2);

2) + *π*

Interestingly, the radial velocity remains independent of the headwall injection sequence, *βn*.

2) + *π*2*r*

This extended form of the Taylor–Culick profile represents a solution for an arbitrary headwall injection pattern *w*0(*r*) that may be prescribed by the proper specification of *βn* through (28). By way of confirmation, the classical Taylor–Culick solution with inert headwall may be readily recovered by setting *βn* = 0 everywhere. The streamline patterns associated with

The framework may be tested using a variable headwall injection profile. To be consistent with the underlying flow assumptions, we employ an axisymmetric function to specify the

*W*c = const uniform

*W*c(1 − *r*2) parabolic

2*πr*2) half cosine

∞ ∑ *n*=0

> ∞ ∑ *n*=0

2*πr*

2*πr*

0.3

� 1 0 0.4 0.5 0.6 0.7 0.8 0.9

*w*0(*r*) cos *χn r dr* (28)

�

(2*n* + 1)*βn* cos *χn* (30)

(2*n* + 1)2*βn* sin *χn* (31)

sin *χn* (29)

(32)

*w*0(*r*) cos *χn r dr*

0.1

0.1

Fig. 3. Streamline patterns corresponding to the classic Taylor–Culick profile with no

(2*n* + 1)*π*

∞ ∑ *n*=0

*βn* = 4

With *βn* in hand, the streamfunction is fully determined, namely,

The radial and axial velocities follow and these may be expressed as

*u*(*r*) = −*r*−1 sin( 1

*w*(*r*, *z*) = *πz* cos( 1

Ω(*r*, *z*) = *π*2*rz* sin( 1

*w*0(*r*) =

⎧ ⎪⎪⎪⎨

⎪⎪⎪⎩

*W*c cos( 1

this historical benchmark are illustrated in Figure 3.

**2.5 Axisymmetric headwall injection profiles**

injection pattern at *z* = 0, namely,

2*πr* 2) +

*ψ*(*r*, *z*) = *z* sin( 1

Finally, the vorticity may be deduced from

0

headwall injection.

or

1

r

$$\psi(r,z) = z\sin(\frac{1}{2}\pi r^2) + \frac{4W\_{\mathbb{C}}}{\pi^2} \sum\_{n=0}^{\infty} \frac{(-1)^n}{(2n+1)^2} \sin \chi\_n \tag{34}$$

The axial velocity and vorticity may be easily determined to be

$$w(r,z) = \pi z \cos(\frac{1}{2}\pi r^2) + \frac{4W\_{\odot}}{\pi} \sum\_{n=0}^{\infty} \frac{(-1)^n}{(2n+1)} \cos \chi\_n \tag{35}$$

$$
\Omega(r, z) = \pi^2 r z \sin(\frac{1}{2}\pi r^2) \tag{36}
$$

The character of (34) is illustrated in Figure 4. Using *W*c = *U*w = 1, a balance between sidewall and headwall injection causes the streamline originating at the corner (*r* = 1, *z* = 0) to bisect the flow field at an angle of *π*/4 as shown in Figure 4(b). By concentrating on a thin region near the sidewall in Figure 4(c), it may be seen that the solution conforms to the stated boundary conditions. It is also evident that *w*0(*r*) = *W*c = 1 corresponds to a simulated solid propellant grain that is burning evenly along its headwall and sidewall boundaries.

#### **2.5.2 Similarity-conforming cosine injection**

For the cosine injection profile, we use (28) to obtain

$$\beta\_n = \begin{cases} \frac{W\_\mathbf{c}}{\pi} \equiv W\_\mathbf{h}; & n = 0 \\ 0; & \text{otherwise} \end{cases} \tag{37}$$

Using (9), the streamfunction becomes

$$
\psi(r,z) = (z+\mathcal{W}\_{\mathbf{h}})\cos(\frac{1}{2}\pi r^2) \tag{38}
$$

The streamlines associated with the cosine headwall injection case are depicted in Figure 5. Their axial velocity and vorticity correspond to

$$\begin{aligned} w(r,z) &= \pi (z + \mathcal{W}\_{\mathbf{h}}) \cos(\frac{1}{2}\pi r^2); \\ \Omega(r,z) &= \pi^2 r (z + \mathcal{W}\_{\mathbf{h}}) \sin(\frac{1}{2}\pi r^2) \end{aligned} \tag{39}$$

It should be noted that while the solutions derived for most injection profiles are approximate, the one corresponding to the similarity-conforming Berman injection will prove to be exact. This behavior will be discussed in Section 2.6.

0.13

Fig. 6. Streamlines corresponding to parabolic headwall injection with *W*c = 1.

2*πr*

2*πr*

2*πr*

*Q*(*r*, *z*) = �∇ × **u** × **Ω**� = − *∂*

*r*2 *∂ψ ∂z* + 1 *r ∂ψ ∂z ∂*Ω *∂r* − 1 *r ∂ψ ∂r ∂*Ω

*ψ*(*r*, *z*) = *z* sin( 1

*w*(*r*, *z*) = *πz* cos( 1

Ω(*r*, *z*) = *π*2*rz* sin( 1

straightforward to see that *Q* may be calculated from

In terms of the streamfunction and the vorticity, we have

*nψn r*

*∂ψn ∂z* + 1 *r ∂ψn ∂z ∂ ∂r* (*C*2

*∂ψn ∂z* + 1 *r ∂ψn ∂z C*2

> *∂ψn ∂z* + 1 *r*

*Q*(*r*, *z*) = − Ω

The streamlines corresponding to this case are shown in Figure 6.

0.13

0

**2.6 Nonlinear residual error**

residual. Using Ω = Ω*n* = *C*2

*Q* = − 1 *r*2

*Qn* = − *C*2

= − *C*2 *nψn r*

> ∞ ∑ *n*=0 Ω*n* ∞ ∑ *n*=0

1

r

0.25

Internal Flows Driven by Wall-Normal Injection 107

Consequently, the streamfunction, axial velocity, and vorticity may be deduced one-by-one:

2) + 8*W*c *π*3

> 2) + 8*W*c *π*2

> > 2) + 8*W*c *π r* ∞ ∑ *n*=0

To test the accuracy of the solutions presented heretofore, we substitute (25) into (5). Terms that do not entirely cancel are hereafter referred to as the residual error *Q*(*r*, *z*). It is

For each eigensolution given by (29), the vorticity transport equation is fulfilled with zero

*nψn* + *C*2 *n ∂ψn ∂z*

It may hence be seen that the summation of (46) over all eigenmodes will be identically zero if a hypothetical case may be considered for which all eigensolutions coexist independently. In practice, however, the eigensolutions must be taken collectively, and so coupling between eigenmodes must be allowed. The total vorticity and streamfunction must be determined and substituted into the vorticity transport equation. Insertion into (45) requires evaluating

> ∞ ∑ *n*=0

*∂ψn ∂z*

∞ ∑ *n*=0

*∂*Ω*n ∂r* − 1 *r*

∞ ∑ *n*=0

*∂ψn ∂r*

∞ ∑ *n*=0

*∂*Ω*n*

*∂z* (47)

*nrψn*, it is clear that (45) becomes

0 0.2L 0.4L 0.6L 0.8L L

z

0.38

∞ ∑ *n*=0

> ∞ ∑ *n*=0

*∂r*

*nrψn*) − 1 *r ∂ψn ∂r*

(*u*Ω) − *∂*

*∂z*

*∂ ∂z* (*C*2 *nrψn*)

*∂ψn ∂r* − *C*2 *n ∂ψn ∂r*

1

1

1

0.50 0.63 0.75 0.88 1.0

1.1

(2*n* + 1)3 sin *χn* (41)

(2*n* + 1)2 cos *χn* (42)

(2*n* + 1) sin *χn* (43)

(*w*Ω) (44)

*∂z* (45)

*∂ψn ∂z* = 0 (46)

Fig. 4. Streamlines corresponding to uniform headwall injection with *W*c = 1.

Fig. 5. Streamlines corresponding to cosine headwall injection with *W*c = 1.

#### **2.5.3 Parabolic injection**

For the parabolic, laminar-like profile, one may substitute *w*0(*r*) = *W*c(1 − *r*2) into (28) and retrieve

$$
\beta\_n = \frac{8W\_\odot}{(2n+1)^3 \pi^3} \tag{40}
$$

12 Will-be-set-by-IN-TECH

0.30

0 0.2L 0.4L 0.6L 0.8L L

0 0.2L 0.4L 0.6L 0.8L L

0 0.2L 0.4L 0.6L 0.8L L

z (c) Sidewall enlargement confirming the fulfillment of the no-slip

0.26

For the parabolic, laminar-like profile, one may substitute *w*0(*r*) = *W*c(1 − *r*2) into (28) and

*βn* = 8*W*c

0 0.2L 0.4L 0.6L 0.8L L

z

0.40

0.53 0.66 0.79 0.92 1.1 1.2

(2*n* + 1)3*π*3 (40)

Fig. 4. Streamlines corresponding to uniform headwall injection with *W*c = 1.

0.13

Fig. 5. Streamlines corresponding to cosine headwall injection with *W*c = 1.

0.13

z (b) Streamline bisector at the upper left corner

z

(a) Overall streamline pattern

0.45 0.60 0.75 0.90 1.1 1.2 1.3

0.15

0.15

0

1

0.8

1

r

0.999

condition

1

r

0

**2.5.3 Parabolic injection**

retrieve

0.9

r

1

r

Fig. 6. Streamlines corresponding to parabolic headwall injection with *W*c = 1.

Consequently, the streamfunction, axial velocity, and vorticity may be deduced one-by-one:

$$\psi(r,z) = z \sin(\frac{1}{2}\pi r^2) + \frac{8W\_{\mathbb{C}}}{\pi^3} \sum\_{n=0}^{\infty} \frac{1}{(2n+1)^3} \sin \chi\_n \tag{41}$$

$$w(r,z) = \pi z \cos(\frac{1}{2}\pi r^2) + \frac{8W\_0}{\pi^2} \sum\_{n=0}^{\infty} \frac{1}{(2n+1)^2} \cos \chi\_n \tag{42}$$

$$
\Omega(r, z) = \pi^2 r z \sin(\frac{1}{2}\pi r^2) + \frac{8W\_{\text{c}}}{\pi} r \sum\_{n=0}^{\infty} \frac{1}{(2n+1)} \sin \chi\_n \tag{43}
$$

The streamlines corresponding to this case are shown in Figure 6.

#### **2.6 Nonlinear residual error**

To test the accuracy of the solutions presented heretofore, we substitute (25) into (5). Terms that do not entirely cancel are hereafter referred to as the residual error *Q*(*r*, *z*). It is straightforward to see that *Q* may be calculated from

$$Q(r, z) = \left\| \nabla \times \mathbf{u} \times \mathbf{D} \right\| = -\frac{\partial}{\partial r} (u \Omega) - \frac{\partial}{\partial z} (w \Omega) \tag{44}$$

In terms of the streamfunction and the vorticity, we have

$$Q(r,z) = -\frac{\Omega}{r^2} \frac{\partial \psi}{\partial z} + \frac{1}{r} \frac{\partial \psi}{\partial z} \frac{\partial \Omega}{\partial r} - \frac{1}{r} \frac{\partial \psi}{\partial r} \frac{\partial \Omega}{\partial z} \tag{45}$$

For each eigensolution given by (29), the vorticity transport equation is fulfilled with zero residual. Using Ω = Ω*n* = *C*2 *nrψn*, it is clear that (45) becomes

$$\begin{split} Q\_n &= -\frac{\mathbf{C}\_n^2 \psi\_n}{r} \frac{\partial \psi\_n}{\partial z} + \frac{1}{r} \frac{\partial \psi\_n}{\partial z} \frac{\partial}{\partial r} (\mathbf{C}\_n^2 r \psi\_n) - \frac{1}{r} \frac{\partial \psi\_n}{\partial r} \frac{\partial}{\partial z} (\mathbf{C}\_n^2 r \psi\_n) \\\\ &= -\frac{\mathbf{C}\_n^2 \psi\_n}{r} \frac{\partial \psi\_n}{\partial z} + \frac{1}{r} \frac{\partial \psi\_n}{\partial z} \mathbf{C}\_n^2 \psi\_n + \mathbf{C}\_n^2 \frac{\partial \psi\_n}{\partial z} \frac{\partial \psi\_n}{\partial r} - \mathbf{C}\_n^2 \frac{\partial \psi\_n}{\partial r} \frac{\partial \psi\_n}{\partial z} = 0 \end{split} \tag{46}$$

It may hence be seen that the summation of (46) over all eigenmodes will be identically zero if a hypothetical case may be considered for which all eigensolutions coexist independently. In practice, however, the eigensolutions must be taken collectively, and so coupling between eigenmodes must be allowed. The total vorticity and streamfunction must be determined and substituted into the vorticity transport equation. Insertion into (45) requires evaluating

$$Q = -\frac{1}{r^2} \sum\_{n=0}^{\infty} \Omega\_n \sum\_{n=0}^{\infty} \frac{\partial \psi\_n}{\partial z} + \frac{1}{r} \sum\_{n=0}^{\infty} \frac{\partial \psi\_n}{\partial z} \sum\_{n=0}^{\infty} \frac{\partial \Omega\_n}{\partial r} - \frac{1}{r} \sum\_{n=0}^{\infty} \frac{\partial \psi\_n}{\partial r} \sum\_{n=0}^{\infty} \frac{\partial \Omega\_n}{\partial z} \tag{47}$$

In this case, the residual is undefined because the alternating sequence of increasing terms in (56) diverges. This may be corroborated by the nature of the uniform profile known for its

Internal Flows Driven by Wall-Normal Injection 109

In all cases for which the residual converges, the error vanishes along the centerline and at the chamber sidewall. This grants our model the character of a rational approximation. Moreover, because the residual remains independent of *z*, the error that is entailed decreases as we move away from the headwall. This improvement in the streamwise direction makes the approximation more suitable for modeling elongated chambers, such as SRMs. Its behavior near the headwall is consistent with the Taylor–Culick model that is known for its subtle discontinuity at *z* = 0. In all cases considered, the core flow approximations become increasingly more accurate away from the headwall, a condition that is compatible with the parallel flow assumption used in many stability investigations of solid and hybrid rocket flow fields. A similar conclusion is reached by Kurdyumov (2008) whose work confirms the nonlinearity of the vorticity-streamfunction relation in the vicinity of the headwall and its

The steady momentum equation (4b) may be readily solved for the pressure distribution. One

2**u** · **u** −

where *p*0 = *p*(0, 0) represents the centerline pressure at the headwall. To ensure a viable

This identity stands in fulfillment of Clairaut's theorem (Clairaut, 1739; 1740). In terms of the

*∂*2*w ∂r∂z* − *u r ∂w*

In short, (57) will produce an analytical expression for the pressure only when (59) is valid. For the classic Taylor–Culick solution, (59) is identically satisfied and the pressure can be

2*π*2*z*2 − 1

2*π*2*z*2 − *W*c*πz* − 1

For uniform or parabolic injection, axial velocities may be determined only approximately through (35) and (42); as such, the integrability constraint (59) is no longer satisfied. However,

2 *r*−2 sin2( 1

2)**e***r* + *π*(*z* + W*h*) cos( 1

2*πr*

2 *r*−2 sin2( 1

2*πr*

*∂*2 *p ∂r∂z* = *∂*2 *p ∂z∂r*

 *u ∂w ∂r*

*dz* (57)

*∂r* = 0 (59)

2) (60)

2)**e***z* (61)

2*πr*2) (62)

(58)

may start with **u** · ∇**u** = −∇*p* and integrate in two spatial directions to retrieve

*p* = *p*0 − 1

expression for the pressure, the total differential of *p* must be exact, or

*u ∂*2*w ∂r*2 + *w*

*p*(*r*, *z*) = *p*0 − 1

and so (59) is fully secured. Integration of Euler's equation renders

2*πr*

**u** = −*r*−1 sin( 1

*p*(*r*, *z*) = *p*0 − 1

sharp discontinuity at the sidewall.

**2.7 Pressure evaluation**

velocity field, (58) yields

For the cosine profile, we have

integrated into

progressive linearity with successive increases in *z*.

where

⎧ ⎪⎨

⎪⎩

$$\begin{aligned} \psi\_{\boldsymbol{n}} &= (\boldsymbol{a}\_{\boldsymbol{n}} \boldsymbol{z} + \boldsymbol{\beta}\_{\boldsymbol{n}}) \sin \chi\_{\boldsymbol{n}}; \quad & \frac{\partial \psi\_{\boldsymbol{n}}}{\partial \boldsymbol{z}} = \boldsymbol{a}\_{\boldsymbol{n}} \sin \chi\_{\boldsymbol{n}} \\\ \frac{\partial \psi\_{\boldsymbol{n}}}{\partial \boldsymbol{\mu}} &= r \boldsymbol{\Gamma}\_{\boldsymbol{n}} (\boldsymbol{a}\_{\boldsymbol{n}} \boldsymbol{z} + \boldsymbol{\beta}\_{\boldsymbol{n}}) \cos \chi\_{\boldsymbol{n}} \cdot \frac{\partial \boldsymbol{\Omega}\_{\boldsymbol{n}}}{\partial \boldsymbol{\Omega}} = r \boldsymbol{\Gamma}\_{\boldsymbol{n}}^{2} \boldsymbol{\mu} \cdot \sin \chi \end{aligned} \tag{48}$$

*∂r* = *rCn*(*αnz* + *βn*) cos *χn*; *∂z* = *rC*2 *nαn* sin *χn*

Furthermore, for the Taylor–Culick class of solutions, *α*0 = 1 and *αn* = 0, ∀ *n* �= 0. This leaves us with

$$\frac{\partial \psi\_n}{\partial z} = \frac{\partial \psi\_0}{\partial z} = \sin(\frac{1}{2}\pi r^2); \quad \frac{\partial \Omega\_n}{\partial z} = \frac{\partial \Omega\_0}{\partial z} = \mathcal{C}\_0^2 r \frac{\partial \psi\_0}{\partial z} \tag{49}$$

Note that the axial derivatives are solely due to the zeroth eigenmode. This reduces (47) into

$$Q = \frac{\partial \psi\_0}{\partial z} \left( -\frac{1}{r^2} \sum\_{n=0}^{\infty} \Omega\_n + \frac{1}{r} \sum\_{n=0}^{\infty} \frac{\partial \Omega\_n}{\partial r} - \mathbb{C}\_0^2 \sum\_{n=0}^{\infty} \frac{\partial \psi\_n}{\partial r} \right) \tag{50}$$

Finally, noting that

$$\frac{\partial \Omega\_n}{\partial r} = \mathbb{C}\_n^2 \psi\_n + \mathbb{C}\_n^2 r \frac{\partial \psi\_n}{\partial r} \tag{51}$$

we retrieve

$$Q(r) = \frac{\partial \psi\_0}{\partial z} \sum\_{n=0}^{\infty} \left(\mathbb{C}\_n^2 - \mathbb{C}\_0^2\right) \frac{\partial \psi\_n}{\partial r} = \sin(\frac{1}{2}\pi r^2) \, r \sum\_{n=1}^{\infty} \mathbb{C}\_n \beta\_n (\mathbb{C}\_n^2 - \mathbb{C}\_0^2) \cos \chi\_n \tag{52}$$

Equation (52) represents the net residual of the vorticity transport equation due to nonlinear coupling. It is not necessarily zero except for inert (*βn* = 0, ∀*n*) or sinusoidal headwall injection profiles (*βn* = 0, ∀*n* ≥ 1). To further explore the behavior of the residual error, we expand (52) into

$$Q(r) = 4\pi^3 r \sin(\frac{1}{2}\pi r^2) \sum\_{n=1}^{\infty} D\_n \cos\left[\frac{1}{2}(2n+1)\pi r^2\right] \tag{53}$$

where

$$D\_n \equiv \frac{\mathbb{C}\_n}{4\pi^3} \beta\_n (\mathbb{C}\_n^2 - \mathbb{C}\_0^2) \equiv n(n+1)(2n+1)\beta\_n \tag{54}$$

Clearly, the residual error vanishes at *r* = (0, 1) and is otherwise controlled by the behavior of *Dn*. This sequence represents the deviation from the exact solution corresponding to the cosine profile for which *C*2 *n* − *C*2 0 = 0. In the case of no headwall injection, *βn* = *Dn* = 0, thus leading to an exact representation. As *Dn* → 0, the solutions become more accurate. Generally, *βn* �= 0 and so *Dn* will only vanish when *C*2 *n* = *C*2 0. To illustrate this character, we consider two examples, namely, those corresponding to parabolic and uniform injection. For parabolic injection, we find a quickly converging sequence, specifically

$$D\_{n, \text{parabolic}} \sim \frac{n(n+1)}{(2n+1)^2} \xrightarrow[n \to \infty]{} \frac{1}{4} \tag{55}$$

In this case, the residual is sufficiently small, albeit non-vanishing, because of the first few terms in (55). However, for the uniform flow, we get alternating infinity, namely

$$D\_{n, \text{uniform}} \sim (-1)^n \frac{n(n+1)}{(2n+1)} \xrightarrow[n \to \infty]{} \pm \infty \tag{56}$$

14 Will-be-set-by-IN-TECH

Furthermore, for the Taylor–Culick class of solutions, *α*0 = 1 and *αn* = 0, ∀ *n* �= 0. This leaves

Note that the axial derivatives are solely due to the zeroth eigenmode. This reduces (47) into

*∂*Ω*n ∂z* = *∂*Ω0

∞ ∑ *n*=0

*nψn* + *C*2 *nr ∂ψn*

> 2*πr* 2)*r* ∞ ∑ *n*=1

*Dn* cos � 1

*∂r* = sin( 1

2*πr* 2) ∞ ∑ *n*=1

*n* − *C*2

*Dn*,parabolic ∼ *n*(*n* + 1)

terms in (55). However, for the uniform flow, we get alternating infinity, namely

*Dn*,uniform ∼ (−1)*n n*(*n* + 1)

Clearly, the residual error vanishes at *r* = (0, 1) and is otherwise controlled by the behavior of *Dn*. This sequence represents the deviation from the exact solution corresponding to the

thus leading to an exact representation. As *Dn* → 0, the solutions become more accurate.

consider two examples, namely, those corresponding to parabolic and uniform injection. For

In this case, the residual is sufficiently small, albeit non-vanishing, because of the first few

(2*n* + 1)2 −−−→ *n*→∞

Equation (52) represents the net residual of the vorticity transport equation due to nonlinear coupling. It is not necessarily zero except for inert (*βn* = 0, ∀*n*) or sinusoidal headwall injection profiles (*βn* = 0, ∀*n* ≥ 1). To further explore the behavior of the residual error,

*∂*Ω*n ∂r* − *C*2 0 ∞ ∑ *n*=0

*∂z* = *αn* sin *χn*

*∂z* = *C*2 0*r ∂ψ*0

*nαn* sin *χn*

*∂ψn ∂r*

*Cnβn*(*C*2

2 (2*n* + 1)*πr*

0 = 0. In the case of no headwall injection, *βn* = *Dn* = 0,

1

*n* = *C*2

�

*∂r* (51)

*n* − *C*2

2 �

0) ≡ *n*(*n* + 1)(2*n* + 1)*βn* (54)

(2*n* + 1) −−−→ *n*→∞ ±∞ (56)

0. To illustrate this character, we

4 (55)

*∂z* (49)

0) cos *χn* (52)

(48)

(50)

(53)

*∂*Ω*n ∂z* = *rC*2

*ψn* = (*αnz* + *βn*) sin *χn*; *∂ψn*

*∂r* = *rCn*(*αnz* + *βn*) cos *χn*;

2*πr* 2);

Ω*n* + 1 *r*

∞ ∑ *n*=0

*∂*Ω*n ∂r* = *C*2

*∂z* = sin( 1

� − 1 *r*2

where

us with

Finally, noting that

we expand (52) into

cosine profile for which *C*2

*Q*(*r*) = *∂ψ*0 *∂z*

we retrieve

where

⎧ ⎪⎨

⎪⎩

*∂ψn ∂z* = *∂ψ*0

*Q* = *∂ψ*0 *∂z*

> ∞ ∑ *n*=0

� *C*2 *n* − *C*2 0 � *∂ψn*

*Q*(*r*) = 4*π*3*r* sin( 1

4*π*3 *βn*(*C*2

parabolic injection, we find a quickly converging sequence, specifically

*Dn* ≡ *Cn*

*n* − *C*2

Generally, *βn* �= 0 and so *Dn* will only vanish when *C*2

*∂ψn*

In this case, the residual is undefined because the alternating sequence of increasing terms in (56) diverges. This may be corroborated by the nature of the uniform profile known for its sharp discontinuity at the sidewall.

In all cases for which the residual converges, the error vanishes along the centerline and at the chamber sidewall. This grants our model the character of a rational approximation. Moreover, because the residual remains independent of *z*, the error that is entailed decreases as we move away from the headwall. This improvement in the streamwise direction makes the approximation more suitable for modeling elongated chambers, such as SRMs. Its behavior near the headwall is consistent with the Taylor–Culick model that is known for its subtle discontinuity at *z* = 0. In all cases considered, the core flow approximations become increasingly more accurate away from the headwall, a condition that is compatible with the parallel flow assumption used in many stability investigations of solid and hybrid rocket flow fields. A similar conclusion is reached by Kurdyumov (2008) whose work confirms the nonlinearity of the vorticity-streamfunction relation in the vicinity of the headwall and its progressive linearity with successive increases in *z*.

#### **2.7 Pressure evaluation**

The steady momentum equation (4b) may be readily solved for the pressure distribution. One may start with **u** · ∇**u** = −∇*p* and integrate in two spatial directions to retrieve

$$p = p\_0 - \frac{1}{2}\mathbf{u} \cdot \mathbf{u} - \int u \frac{\partial w}{\partial r} \, dz \tag{57}$$

where *p*0 = *p*(0, 0) represents the centerline pressure at the headwall. To ensure a viable expression for the pressure, the total differential of *p* must be exact, or

$$\frac{\partial^2 p}{\partial r \partial z} = \frac{\partial^2 p}{\partial z \partial r} \tag{58}$$

This identity stands in fulfillment of Clairaut's theorem (Clairaut, 1739; 1740). In terms of the velocity field, (58) yields

$$
\mu \frac{\partial^2 w}{\partial r^2} + w \frac{\partial^2 w}{\partial r \partial z} - \frac{\mu}{r} \frac{\partial w}{\partial r} = 0 \tag{59}
$$

In short, (57) will produce an analytical expression for the pressure only when (59) is valid. For the classic Taylor–Culick solution, (59) is identically satisfied and the pressure can be integrated into

$$p(r,z) = p\_0 - \frac{1}{2}\pi^2 z^2 - \frac{1}{2}r^{-2}\sin^2(\frac{1}{2}\pi r^2) \tag{60}$$

For the cosine profile, we have

$$\mathbf{u} = -r^{-1}\sin(\frac{1}{2}\pi r^2)\mathbf{e}\_r + \pi(z + \mathcal{W}\_\hbar)\cos(\frac{1}{2}\pi r^2)\mathbf{e}\_z\tag{61}$$

and so (59) is fully secured. Integration of Euler's equation renders

$$p(r,z) = p\_0 - \frac{1}{2}\pi^2 z^2 - \mathcal{W}\_0 \pi z - \frac{1}{2}r^{-2}\sin^2(\frac{1}{2}\pi r^2) \tag{62}$$

For uniform or parabolic injection, axial velocities may be determined only approximately through (35) and (42); as such, the integrability constraint (59) is no longer satisfied. However,

0 0.5 1

0 0.5 1

0 0.5 1

(e) uniform

(c) parabolic

(a) cosine

z/L

z/L

z/L

uniform injection. Curves are shown for *z*/*L* = 0.1, 0.3, 0.5, 0.7, and 0.9.

Fig. 7. Comparison between analytical (—) for (65), (−·−) for (69) and numerical simulations () for the centerline pressure using (a, b) cosine, (c, d) parabolic, and (e, f)

1.0

3.0 pc

> Wc = 10

1.0

1.0

taken at 1.6 m × 0.1 m and the wall injection velocity is taken at 10 m/s for the simulated SRM. The boundary condition at the sidewall is specified as a velocity inlet to closely mimic the mathematical model where injection is imposed uniformly along the grain surface. The headwall is also defined as a velocity inlet. On the right-hand-side of the domain, a pressure outlet boundary condition is prescribed where the exit pressure is set to be atmospheric as in the case of sea level testing. Although an outflow boundary condition can also be imposed at the downstream section, it is avoided here to avert the possible case of partially developed flow (White, 2005). The difference between an outflow and a pressure outlet boundary

2.0

3.0

Wc = 10

pc

2.0

2.0

3.0 pc

Internal Flows Driven by Wall-Normal Injection 111

Wc = 10

0 0.5 1

0 0.5 1

0 0.5 1

(f) uniform

(d) parabolic

(b) cosine

z/L

z/L

using p = ∑ pn

z/L

 CFD p ∑ pn

1.0

1.0

1.0

1.5

2.0

Wc = 1

2.5 pc

1.5

2.0

Wc = 1

2.5 pc

pc

1.5

2.0

Wc = 1

2.5

pc

along the centerline, the constraint remains valid. At *r* = 0, the radial velocity shared by both injection profiles vanishes in view of

$$\ln(0, z) = \lim\_{r \to 0} r^{-1} \sin(\frac{1}{2}\pi r^2) = 0 \tag{63}$$

As for the axial velocities, they become equal viz.

$$w\_{\text{uniform}}(0, z) = w\_{\text{parabolic}}(0, z) = \pi z + W\_{\odot} \tag{64}$$

This enables us to integrate (57) and collect

$$p(0, z) = p\_0 - \frac{1}{2}\pi^2 z^2 - \pi z \mathcal{W}\_{\mathbb{C}} \tag{65}$$

Interestingly, all injection profiles generate the same expression for the centerline pressure. To overcome the pitfalls of pressure integrability of a non-exact velocity, approximate representations of *p* may be sought based on a linear expansion that becomes increasingly more accurate as *z* is increased. This is

$$p(r,z) = \sum\_{n=0}^{\infty} p\_n(r,z) \tag{66}$$

where *pn* is the pressure corresponding to the *n*th eigenmode in (20). Integration of the pressure in this case is possible because each eigensolution given by *ψn*(*r*, *z*) consists of an exact solution of the Euler equations that directly satisfies (59). Using

$$\begin{aligned} u\_n &= -\alpha\_n r^{-1} \sin \chi\_n; \\ w\_n &= (2n+1)\pi (\alpha\_n z + \beta\_n) \cos \chi\_n \end{aligned} \tag{67}$$

one can integrate for the pressure to find

$$p\_n(r,z) = p\_0 - \frac{1}{2}(2n+1)^2 \pi^2 a\_n^2 z^2 - (2n+1)^2 \pi^2 a\_n \beta\_n z - \frac{1}{2} a\_n^2 r^{-2} \sin^2 \chi\_n \tag{68}$$

or

$$p(r,z) = \sum\_{n=0}^{\infty} p\_n = p\_0 - \frac{1}{2}\pi^2 z^2 - \beta\_0 \pi^2 z - \frac{1}{2}r^{-2}\sin^2(\frac{1}{2}\pi r^2) \tag{69}$$

As shown in Figure 7, this linear approximation stands in better agreement with the numerical data than the result obtained in (65) for the pressure based on the total velocity. This may be connected to the increasing accuracy associated with a linear vorticity-streamfunction assumption and the superposition of eigensolutions with successive increases in *z* (Kurdyumov, 2008).

#### **2.8 Numerical verification**

So far we have introduced an approximate Euler solution for the Taylor–Culick profile with variable headwall injection. By way of confirmation, an inviscid numerical solution is presented for the mean flow using three illustrative headwall injection profiles. Our simulations are carried out using a finite-volume CFD solver. The targeted flow is that corresponding to a rocket motor with an average sidewall Mach number of 0.03 and strictly inviscid conditions. For the sake of comparison, the working fluid is taken to be ambient air. The aspect ratio of the domain is set at *L* = 16. The actual length and radius are 16 Will-be-set-by-IN-TECH

along the centerline, the constraint remains valid. At *r* = 0, the radial velocity shared by both

*r*−1 sin( 1

2*πr*

*w*uniform(0, *z*) = *w*parabolic(0, *z*) = *πz* + *W*c (64)

*wn* = (2*n* + 1)*π*(*αnz* + *βn*) cos *χn* (67)

*nz*2 − (2*n* + 1)2*π*2*αn βnz* − 1

2*π*2*z*2 − *β*0*π*2*z* − 1

As shown in Figure 7, this linear approximation stands in better agreement with the numerical data than the result obtained in (65) for the pressure based on the total velocity. This may be connected to the increasing accuracy associated with a linear vorticity-streamfunction assumption and the superposition of eigensolutions with successive

So far we have introduced an approximate Euler solution for the Taylor–Culick profile with variable headwall injection. By way of confirmation, an inviscid numerical solution is presented for the mean flow using three illustrative headwall injection profiles. Our simulations are carried out using a finite-volume CFD solver. The targeted flow is that corresponding to a rocket motor with an average sidewall Mach number of 0.03 and strictly inviscid conditions. For the sake of comparison, the working fluid is taken to be ambient air. The aspect ratio of the domain is set at *L* = 16. The actual length and radius are

2) = 0 (63)

2*π*2*z*2 − *πzW*c (65)

*pn*(*r*, *z*) (66)

2 *α*2

2*πr*

2 *r*−2 sin2( 1

*nr*−2 sin2 *χn* (68)

2) (69)

*u*(0, *z*) = lim

*r*→0

*p*(0, *z*) = *p*0 − 1

*p*(*r*, *z*) =

*un* = −*αnr*−1 sin *χn*;

exact solution of the Euler equations that directly satisfies (59). Using

2 (2*n* + 1)2*π*2*α*2

*pn* = *p*0 − 1

∞ ∑ *n*=0

Interestingly, all injection profiles generate the same expression for the centerline pressure. To overcome the pitfalls of pressure integrability of a non-exact velocity, approximate representations of *p* may be sought based on a linear expansion that becomes increasingly

> ∞ ∑ *n*=0

where *pn* is the pressure corresponding to the *n*th eigenmode in (20). Integration of the pressure in this case is possible because each eigensolution given by *ψn*(*r*, *z*) consists of an

injection profiles vanishes in view of

As for the axial velocities, they become equal viz.

This enables us to integrate (57) and collect

more accurate as *z* is increased. This is

one can integrate for the pressure to find

*pn*(*r*, *z*) = *p*0 − 1

increases in *z* (Kurdyumov, 2008).

**2.8 Numerical verification**

*p*(*r*, *z*) =

or

Fig. 7. Comparison between analytical (—) for (65), (−·−) for (69) and numerical simulations () for the centerline pressure using (a, b) cosine, (c, d) parabolic, and (e, f) uniform injection. Curves are shown for *z*/*L* = 0.1, 0.3, 0.5, 0.7, and 0.9.

taken at 1.6 m × 0.1 m and the wall injection velocity is taken at 10 m/s for the simulated SRM. The boundary condition at the sidewall is specified as a velocity inlet to closely mimic the mathematical model where injection is imposed uniformly along the grain surface. The headwall is also defined as a velocity inlet. On the right-hand-side of the domain, a pressure outlet boundary condition is prescribed where the exit pressure is set to be atmospheric as in the case of sea level testing. Although an outflow boundary condition can also be imposed at the downstream section, it is avoided here to avert the possible case of partially developed flow (White, 2005). The difference between an outflow and a pressure outlet boundary

0 1

0 1

0 1

Fig. 8. Comparison between analytical (—) and numerical simulations (◦) for the axial velocity using (a, b) cosine, (c, d) parabolic, and (e, f) uniform injection. Curves are shown for

> 2**u**2 *n* = 1 2 (*u*2 *n* + *v*2

r

*En*(*r*, *θ*, *z*) = 1

*un* = −*r*−1*αn* sin *χn*; *vn* = 0 *wn* = *παnz*(2*n* + 1) cos *χn*

r

0

0

0

20

40

60

80 w Wc

20

40

60 w Wc

20

40

60 w Wc

Internal Flows Driven by Wall-Normal Injection 113

= 10

= 10

= 10

0 1

0 1

0 1

*n*) (71)

2 )*πr*2 (72)

r

r

r

(b) cosine

(d) parabolic

(f) uniform

*n* + *w*2

*χn* ≡ (*n* + 1

r

(a) cosine

(c) parabolic

(e) uniform

where each mode is an exact solution that is given by

0

0

0

*z*/*L* = 0.1, 0.3, 0.5, 0.7, and 0.9.

each eigensolution using

10

20

30

40

50 w Wc

10

20

30

40

50 w Wc

10

20

30

40

50

w Wc

= 1

= 1

= 1

condition is that, in the latter case, the exit pressure is fixed at the boundary. The domain is meshed into 589,824 equally spaced control volumes consisting of 3072 × 192 cells. While the Quadratic Upwind Interpolation for Convection Kinematics (QUICK) scheme is called upon for spatial discretization, the Semi Implicit Method for Pressure Linked Equation (SIMPLE) algorithm is used to resolve the pressure–velocity coupling.

Results for the inviscid simulations are shown in Figures 7–8. These are carried out for *W*c = 1 and 10 (i.e., solid and hybrid motors); their purpose is to show the streamwise evolution of the axial velocity, vorticity, and centerline pressure at *z*/*L* = 0.1, 0.3, . . . , 0.9. It may be seen that the agreement with the computations is excellent except in the case of uniform injection with a large *W*c. This may be attributed to the discontinuity that the uniform injection profile experiences at the sidewall. Furthermore, according to (56), we expect the residual error to be large. These limited numerical runs reaffirm the viability of the analytical approximations as simple predictive tools.

#### **3. Generalized Taylor-Culick formulation**

In Section 2.1, we presented a mean flow model for solid and hybrid rocket motors that could assimilate a rather arbitrary headwall injection profile based on a specific form of *βn*. Initially, the solutions were obtained in series form that depended on two parameters, *αn* and *βn*, the sidewall and headwall injection sequences. While *βn* was prescribed by the headwall injection pattern, the choice of *αn* appeared to be flexible provided that the constraint given by (22) remained satisfied. In this section, we follow Majdalani & Saad (2007a) by applying the Lagrangian optimization technique to the total kinetic energy of the generalized Taylor–Culick solution to the extent of producing a variational constraint on *αn* (see also Saad & Majdalani, 2010). After some effort, two types of solutions will be identified with increasing or decreasing kinetic energies; of the two families, the Taylor–Culick model will be recovered as a special case. The new approximations will be shown to exhibit velocity profiles with energy dependent curvatures that are reminiscent of turbulent or compressible motions. In practice, steeper profiles have been observed in either experimental or numerical tests, particularly in the presence of intense levels of acoustic energy (Apte & Yang, 2000; 2001; 2002). Interestingly, the energy-based models will range from irrotational to rotational fields with increasing vorticity, thus covering a wide spectrum of admissible motions that observe the problem's physical requirements. A second law analysis will be later used to test the physicality of these solutions and establish the Taylor–Culick motion as an equilibrium state to which all profiles will tend to converge.

#### **3.1 Kinetic energy optimization**

As shown in Section 2.1, the sidewall injection sequence must observe a key constraint associated with the wall-normal injection velocity:

$$\sum\_{n=0}^{\infty} (-1)^n a\_n = 1 \tag{70}$$

Clearly, numerous sequences of *αn* exist that can be made to satisfy (70). One of these choices may be arrived at by optimizing the total volumetric kinetic energy in the chamber. The guiding principle is based on the hypothesis that a flow may follow the path of least or most energy expenditure. To test this behavior, we evaluate the local kinetic energy at (*r*, *θ*, *z*) for 18 Will-be-set-by-IN-TECH

condition is that, in the latter case, the exit pressure is fixed at the boundary. The domain is meshed into 589,824 equally spaced control volumes consisting of 3072 × 192 cells. While the Quadratic Upwind Interpolation for Convection Kinematics (QUICK) scheme is called upon for spatial discretization, the Semi Implicit Method for Pressure Linked Equation (SIMPLE)

Results for the inviscid simulations are shown in Figures 7–8. These are carried out for *W*c = 1 and 10 (i.e., solid and hybrid motors); their purpose is to show the streamwise evolution of the axial velocity, vorticity, and centerline pressure at *z*/*L* = 0.1, 0.3, . . . , 0.9. It may be seen that the agreement with the computations is excellent except in the case of uniform injection with a large *W*c. This may be attributed to the discontinuity that the uniform injection profile experiences at the sidewall. Furthermore, according to (56), we expect the residual error to be large. These limited numerical runs reaffirm the viability of the analytical approximations as

In Section 2.1, we presented a mean flow model for solid and hybrid rocket motors that could assimilate a rather arbitrary headwall injection profile based on a specific form of *βn*. Initially, the solutions were obtained in series form that depended on two parameters, *αn* and *βn*, the sidewall and headwall injection sequences. While *βn* was prescribed by the headwall injection pattern, the choice of *αn* appeared to be flexible provided that the constraint given by (22) remained satisfied. In this section, we follow Majdalani & Saad (2007a) by applying the Lagrangian optimization technique to the total kinetic energy of the generalized Taylor–Culick solution to the extent of producing a variational constraint on *αn* (see also Saad & Majdalani, 2010). After some effort, two types of solutions will be identified with increasing or decreasing kinetic energies; of the two families, the Taylor–Culick model will be recovered as a special case. The new approximations will be shown to exhibit velocity profiles with energy dependent curvatures that are reminiscent of turbulent or compressible motions. In practice, steeper profiles have been observed in either experimental or numerical tests, particularly in the presence of intense levels of acoustic energy (Apte & Yang, 2000; 2001; 2002). Interestingly, the energy-based models will range from irrotational to rotational fields with increasing vorticity, thus covering a wide spectrum of admissible motions that observe the problem's physical requirements. A second law analysis will be later used to test the physicality of these solutions and establish the Taylor–Culick motion as an equilibrium state

As shown in Section 2.1, the sidewall injection sequence must observe a key constraint

Clearly, numerous sequences of *αn* exist that can be made to satisfy (70). One of these choices may be arrived at by optimizing the total volumetric kinetic energy in the chamber. The guiding principle is based on the hypothesis that a flow may follow the path of least or most energy expenditure. To test this behavior, we evaluate the local kinetic energy at (*r*, *θ*, *z*) for

(−1)*nαn* = 1 (70)

∞ ∑ *n*=0

algorithm is used to resolve the pressure–velocity coupling.

simple predictive tools.

**3. Generalized Taylor-Culick formulation**

to which all profiles will tend to converge.

associated with the wall-normal injection velocity:

**3.1 Kinetic energy optimization**

Fig. 8. Comparison between analytical (—) and numerical simulations (◦) for the axial velocity using (a, b) cosine, (c, d) parabolic, and (e, f) uniform injection. Curves are shown for *z*/*L* = 0.1, 0.3, 0.5, 0.7, and 0.9.

each eigensolution using

$$E\_n(r, \theta, z) = \frac{1}{2} \mathbf{u}\_n^2 = \frac{1}{2} (\mu\_n^2 + v\_n^2 + w\_n^2) \tag{71}$$

where each mode is an exact solution that is given by

$$\begin{cases} u\_{\rm nl} = -r^{-1} u\_{\rm nl} \sin \chi\_{\rm nl}; \quad v\_{\rm nl} = 0\\ w\_{\rm nl} = \pi \alpha\_{\rm nl} z (2n+1) \cos \chi\_{\rm nl} \end{cases} \quad \chi\_{\rm n} \equiv (n + \frac{1}{2}) \pi r^2 \tag{72}$$

to introduce the constrained energy function

*∂*G *∂αn*

Then, through substitution into (78), one retrieves

∞ ∑ *n*=0 = 1 6*π*3*L*3

Equation (77) may be used to extract *αn* in terms of *λ* such that

*αn* = (−1)*n*

(−1)*nαn* = 1

*∂*G *∂λ* =

variables to obtain

**3.2 Critical length**

be accomplished by setting

and

G(*α*0, *α*1, *α*2,..., *λ*) = *E*V + *λ*

imposing ∇G(*α*0, *α*1, *α*2,..., *λ*) = 0. In shorthand notation, we put

 ∞ ∑ *n*=0

where *λ* is a Lagrangian multiplier. Equation (75) can be maximized or minimized by

Internal Flows Driven by Wall-Normal Injection 115

Naturally, the constrained energy function may be differentiated with respect to each of its

*αndn π*2*L*2 (*an* + *π*−2*L*−2*dn*)−1

1 (*an* + *π*−2*L*−2*dn*) = *N*

∞ ∑ *i*=0

*αnan* +

∞ ∑ *n*=0

*αn* = − 6(−1)*nλ*

*λ* = − *π*3*L*3

Finally, when *λ* is inserted into (79), a general solution for *αn* emerges, specifically

(*an* + *π*−2*L*−2*dn*)*N* ; *N* =

Clearly, (81) satisfies the fundamental constraint which, by inspection, returns

∞ ∑ *n*=0

Some values of *αn* are posted in Table 1 at four different aspect ratios corresponding to *L* = 1, 5, 10, and 100. With this expression at hand, the total energy *E*V is completely determined.

Equation (74) can be normalized by *L*3 and simplified into an energy density form. This can

By plotting E versus *L* in Figure 10, it can be seen that E approaches a constant asymptotic value of E∞ = 2*π*/3. Granted this behavior, a critical aspect ratio *L*cr may be defined beyond

*N*

6 ∞ ∑ *n*=0 (−1)*nαn* − 1

∇G(*αn*, *λ*) = 0; *n* ∈ **N**0 (76)

+ (−1)*nλ* = 0 (77)

(−1)*nαn* − 1 = 0 (78)

*π*3*L*3(*an* + *π*−2*L*−2*dn*) (79)

1 *ai* + *π*−2*L*−2*di*

E = *E*V /*L*3 (83)

(75)

(80)

(81)

*N* = 1 (82)

Fig. 9. Comparison between analytical (—) and numerical simulations (◦) for the vorticity magnitude using (a, b) cosine and (c, d) parabolic injection. Curves are shown for *z*/*L* = 0.1, 0.3, 0.5, 0.7, and 0.9.

We now define the cumulative local kinetic energy as the sum of contributions from individual eigensolutions. This can be written as

$$E(r,\theta,z) = \sum\_{n=0}^{\infty} E\_n(r,\theta,z) = \frac{1}{2} \sum\_{n=0}^{\infty} \left[ a\_n^2 r^{-2} \sin^2 \chi\_n + \pi^2 a\_n^2 z^2 (2n+1)^2 \cos^2 \chi\_n \right] \tag{73}$$

Subsequently, the total kinetic energy in a chamber of volume V may be calculated by integrating the local kinetic energy over the length and chamber cross-section,

$$E\nu = \iiint\_{V} \mathbf{E}(r, \theta, \mathbf{z}) r \, \mathrm{d}r \, \mathrm{d}\theta \, \mathrm{d}z = \pi \sum\_{n=0}^{\infty} \int\_{0}^{L} \int\_{0}^{1} a\_{n}^{2} \left[ \frac{\sin^{2} \chi\_{n}}{r^{2}} + \pi^{2} z^{2} (2n+1)^{2} \cos^{2} \chi\_{n} \right] r \, \mathrm{d}r \, \mathrm{d}z$$

Straightforward evaluation and simplification yield

$$E\_{\mathcal{V}} = \frac{1}{12}\pi^3 L^3 \sum\_{n=0}^{\infty} (a\_n^2 a\_n + a\_n^2 \pi^{-2} L^{-2} d\_n); \quad \begin{cases} a\_{\mathcal{U}} = (2n+1)^2\\ d\_{\mathcal{U}} = 3\text{Cin}[(2n+1)\pi] \end{cases} \tag{74}$$

where Cin(*x*) ≡ *x* 0 (1 − cos *t*)*t* −1 d*t* is the Entire Cosine Integral. At this point, one may seek the extremum of (74) subject to the fundamental constraint (70). The latter enables us to introduce the constrained energy function

$$\mathcal{G}(\mathfrak{a}\_{0}, \mathfrak{a}\_{1}, \mathfrak{a}\_{2}, \dots, \lambda) = E\_{\mathcal{V}} + \lambda \left[ \sum\_{n=0}^{\infty} (-1)^{n} \mathfrak{a}\_{n} - 1 \right] \tag{75}$$

where *λ* is a Lagrangian multiplier. Equation (75) can be maximized or minimized by imposing ∇G(*α*0, *α*1, *α*2,..., *λ*) = 0. In shorthand notation, we put

$$\nabla \mathcal{G}(\mathfrak{a}\_{\mathfrak{n}}, \lambda) = 0; \quad \mathfrak{n} \in \mathbb{N}\_0 \tag{76}$$

Naturally, the constrained energy function may be differentiated with respect to each of its variables to obtain

$$\frac{\partial \mathcal{G}}{\partial a\_{\rm nl}} = \frac{1}{5} \pi^3 L^3 \left( a\_{\rm nl} a\_{\rm nl} + \frac{a\_{\rm nl} d\_{\rm nl}}{\pi^2 L^2} \right) + (-1)^n \lambda = 0 \tag{77}$$

and

20 Will-be-set-by-IN-TECH

0

0

Fig. 9. Comparison between analytical (—) and numerical simulations (◦) for the vorticity magnitude using (a, b) cosine and (c, d) parabolic injection. Curves are shown for *z*/*L* = 0.1,

We now define the cumulative local kinetic energy as the sum of contributions from individual

Subsequently, the total kinetic energy in a chamber of volume V may be calculated by

 1 0 *α*2 *n* sin2 *χn*

*nπ*−2*L*−2*dn*);

seek the extremum of (74) subject to the fundamental constraint (70). The latter enables us

*nr*−2 sin2 *χn* + *π*2*α*2

50

100

150

200 Wc

= 10

50

100

150

200 Wc

= 10

0 1

0 1

*nz*2(2*n* + 1)2 cos2 *χn*

*dn* = 3Cin[(2*n* + 1)*π*] (74)

*r*2 + *π*2*z*2(2*n* + 1)2 cos2 *χn*

*an* = (2*n* + 1)2

−1 d*t* is the Entire Cosine Integral. At this point, one may

*r* d*r* d*z*

(73)

r (d) parabolic

r (b) cosine

0 1

0 1

*En*(*r*, *θ*, *z*) = 1

2 ∞ ∑ *n*=0 *α*2

integrating the local kinetic energy over the length and chamber cross-section,

 *L* 0

∞ ∑ *n*=0

*nan* + *α*2

(c) parabolic

r

r (a) cosine

0

0

eigensolutions. This can be written as

∞ ∑ *n*=0

*E*(*r*, *θ*, *z*)*r* d*r* d*θ* d*z* = *π*

Straightforward evaluation and simplification yield

(1 − cos *t*)*t*

12*π*3*L*3 ∞ ∑ *n*=0 (*α*2

*E*(*r*, *θ*, *z*) =

0.3, 0.5, 0.7, and 0.9.

*E*V =

where Cin(*x*) ≡

V

*E*V = 1

 *x* 0

50

100

150 Wc

= 1

50

100

150 Wc

= 1

$$\frac{\partial \mathcal{G}}{\partial \lambda} = \sum\_{n=0}^{\infty} (-1)^n a\_n - 1 = 0 \tag{78}$$

Equation (77) may be used to extract *αn* in terms of *λ* such that

$$\alpha\_{\mathfrak{n}} = -\frac{\mathfrak{G}(-1)^{\mathfrak{n}}\lambda}{\pi^3 L^3 (a\_{\mathfrak{n}} + \pi^{-2} L^{-2} d\_{\mathfrak{n}})} \tag{79}$$

Then, through substitution into (78), one retrieves

$$\lambda = -\frac{\pi^3 L^3}{6\sum\_{n=0}^{\infty} (a\_n + \pi^{-2} L^{-2} d\_n)^{-1}}\tag{80}$$

Finally, when *λ* is inserted into (79), a general solution for *αn* emerges, specifically

$$\mathfrak{a}\_{\mathfrak{n}} = \frac{(-1)^{n}}{(a\_{\mathfrak{n}} + \pi^{-2}L^{-2}d\_{\mathfrak{n}})N}; \quad N = \sum\_{i=0}^{\infty} \frac{1}{a\_{i} + \pi^{-2}L^{-2}d\_{i}} \tag{81}$$

Clearly, (81) satisfies the fundamental constraint which, by inspection, returns

$$\sum\_{n=0}^{\infty}(-1)^{n}a\_{n} = \frac{1}{N} \sum\_{n=0}^{\infty} \frac{1}{(a\_{n} + \pi^{-2}L^{-2}d\_{n})} = \frac{N}{N} = 1\tag{82}$$

Some values of *αn* are posted in Table 1 at four different aspect ratios corresponding to *L* = 1, 5, 10, and 100. With this expression at hand, the total energy *E*V is completely determined.

#### **3.2 Critical length**

Equation (74) can be normalized by *L*3 and simplified into an energy density form. This can be accomplished by setting

$$\mathcal{E} = E\_{\mathcal{V}} / L^3 \tag{83}$$

By plotting E versus *L* in Figure 10, it can be seen that E approaches a constant asymptotic value of E∞ = 2*π*/3. Granted this behavior, a critical aspect ratio *L*cr may be defined beyond

0.060

0.010

0.143

0.286

0

1

r

0.99

minimum energy solution.

1

r

0.167

Internal Flows Driven by Wall-Normal Injection 117

0 0.2L 0.4L 0.6L 0.8L L

0 0.2L 0.4L 0.6L 0.8L L

(−1)*n*

The right–oriented mapping arrow '�→' in (87) is used to indicate that the compacted expression is valid inside the domain, 0 ≤ *r* < 1, thus excluding the sidewall. We also remark that, in evaluating (87), the large *L* approximation is used. The corresponding streamfunction, velocity, and vorticity are catalogued in Table 3. The streamlines are shown in Figure 11(a) using solid lines to denote the Taylor–Culick benchmark, and broken lines to describe the

(2*n* + 1)2 sin *χn* �→ *r*

*n*=0(−1)*nαn*

z (b) Sidewall enlargement confirming the fulfillment of the no-slip condition

0.571

0.714

z

(a) Overall streamline pattern

0.429

Fig. 11. Streamlines corresponding to the minimum kinetic energy solution (87) for an inert

headwall. Solid lines: Taylor-Culick; broken lines: minimum energy solution.

*π*2 *z*

∞ ∑ *n*=0

*m αm* ∑*m*

Table 2. Convergence of the sidewall injection sequence *αn* when *L* → ∞.

0 0.8105 0.8105 1 -0.0900 0.9006 2 0.0324 0.9330 3 -0.0165 0.9495 4 0.0100 0.9596 5 -0.0066 0.9663 ∞ 0.0000 1.0000

an inert headwall, the minimum energy approximation reduces to

*ψ*(*r*, *z*) = 8

0.333

0.500 0.667

0.833

0.857

2*z* (87)

Table 1. Convergence of the sidewall injection sequence *αn* for *L* = 1, 5, 10, and 100.

Fig. 10. Kinetic energy density variation with *L*. Note that for *L* > 6.7, the energy density will be within 1% of its final asymptotic value E∞.

which the energy density will vary by less than one percent from its asymptotic value E∞. We therefore set

$$
\mathcal{E}\_{\rm CF} - \mathcal{E}\_{\rm \infty} \le 0.01 \,\mathcal{E}\_{\rm \infty} \tag{84}
$$

For a chamber of length *L* ≥ *L*cr, one may evaluate the limiting behavior of (81) by taking *L* → ∞. For SRMs with inert headwalls, the critical length is found to be 6.7. In practice, most SRMs are designed with an aspect ratio that exceeds 20 and so the assumption of a large *L* may be safely employed in describing their flow fields. With this simplification, the expression for *αn* collapses into

$$\lim\_{L \to \infty} a\_{ll} = (-1)^{n} \left( a\_{ll} \sum\_{i=0}^{\infty} \frac{1}{a\_{i}} \right)^{-1} = \frac{8(-1)^{n}}{\pi^{2}(2n+1)^{2}} \tag{85}$$

Note that (85) identically satisfies the sidewall constraint viz.

$$\sum\_{n=0}^{\infty}(-1)^{n}\alpha\_{n} = \frac{8}{\pi^{2}}\sum\_{n=0}^{\infty}\frac{1}{(2n+1)^{2}} = 1\tag{86}$$

The large-*L* approximation of *αn* quickly converges as illustrated in Table 2.

#### **3.3 Least kinetic energy solution**

While the use of Lagrangian multipliers enables us to identify the problem's extremum, straightforward substitution of (81) into (74) allows us to compare the energy content of the present approximation to that of Taylor–Culick's. We find that the extremum obtained through Lagrangian optimization corresponds to the solution with least kinetic energy. Given 22 Will-be-set-by-IN-TECH

*n L* = 1 *L* = 5 *L* = 10 *L* = 100 0.7524 0.8095 0.8115 0.8121 -0.1146 -0.0914 -0.0905 -0.0902 0.0434 0.0329 0.0326 0.0324 -0.0225 -0.0168 -0.0166 -0.0165 0.0137 0.0101 0.0100 0.0100 -0.0092 -0.0068 -0.0067 -0.0067

0 3 6 912 15

L

Fig. 10. Kinetic energy density variation with *L*. Note that for *L* > 6.7, the energy density will

which the energy density will vary by less than one percent from its asymptotic value E∞. We

For a chamber of length *L* ≥ *L*cr, one may evaluate the limiting behavior of (81) by taking *L* → ∞. For SRMs with inert headwalls, the critical length is found to be 6.7. In practice, most SRMs are designed with an aspect ratio that exceeds 20 and so the assumption of a large *L* may be safely employed in describing their flow fields. With this simplification, the expression for

> 1 *ai*

∞ ∑ *n*=0

−1

1

 *an* ∞ ∑ *i*=0

(−1)*nαn* = 8

The large-*L* approximation of *αn* quickly converges as illustrated in Table 2.

*π*2

While the use of Lagrangian multipliers enables us to identify the problem's extremum, straightforward substitution of (81) into (74) allows us to compare the energy content of the present approximation to that of Taylor–Culick's. We find that the extremum obtained through Lagrangian optimization corresponds to the solution with least kinetic energy. Given

: − ≤ ∞ ∞ L •• • cr cr 0.01

2 3 • ∞ =

Ecr − E∞ ≤ 0.01 E∞ (84)

= 8(−1)*n*

*π*2(2*n* + 1)2 (85)

(2*n* + 1)2 = 1 (86)

Table 1. Convergence of the sidewall injection sequence *αn* for *L* = 1, 5, 10, and 100.

2 2.2 2.4 2.6 2.8 3

therefore set

*αn* collapses into

be within 1% of its final asymptotic value E∞.

lim

**3.3 Least kinetic energy solution**

*L*→∞ *αn* = (−1)*n*

Note that (85) identically satisfies the sidewall constraint viz. ∞ ∑ *n*=0

•

Fig. 11. Streamlines corresponding to the minimum kinetic energy solution (87) for an inert headwall. Solid lines: Taylor-Culick; broken lines: minimum energy solution.

an inert headwall, the minimum energy approximation reduces to

$$\psi(r,z) = \frac{8}{\pi^2} z \sum\_{n=0}^{\infty} \frac{(-1)^n}{(2n+1)^2} \sin \chi\_n \mapsto r^2 z \tag{87}$$

The right–oriented mapping arrow '�→' in (87) is used to indicate that the compacted expression is valid inside the domain, 0 ≤ *r* < 1, thus excluding the sidewall. We also remark that, in evaluating (87), the large *L* approximation is used. The corresponding streamfunction, velocity, and vorticity are catalogued in Table 3. The streamlines are shown in Figure 11(a) using solid lines to denote the Taylor–Culick benchmark, and broken lines to describe the minimum energy solution.

Table 2. Convergence of the sidewall injection sequence *αn* when *L* → ∞.

2345678

Fig. 12. Variation of the kinetic energy density with the energy power index for Type I (lower branch) and Type II (upper branch) solutions. These are shown at two aspect ratios, *L* = 10

Internal Flows Driven by Wall-Normal Injection 119

series convergence down to the vorticity. Backward substitution allows us to extract the final

To understand the effect of the energy power index *q* on the kinetic energy density, we use (92) and (74) to plot E versus *q* at two aspect ratios. This plot corresponds to the lower branch of Figure 12 for both *L* = 10 and 20. Interestingly, as *q* → ∞, Taylor–Culick's classic solution

This result identically reproduces Taylor–Culick's expression. All of the Type I solutions derived from (92) possess kinetic energies that are lower than Taylor–Culick's; this explains

To capture solutions with energies that exceed that of Taylor-Culick's, a modified form of *αn*

The key difference here stands in the exclusion of the (−1)*n* multiplier that appears in (89).

(−1)*nBq*

2*πr*2). In practice, profiles with *q* ≥ 5 will be indiscernible from Taylor–Culick's as their energies will then differ by less than one percent. The most distinct solutions will correspond to *q* = 2, 3, and 4 with energies that are 81.1, 91.7, and 97.3 percent of Taylor–Culick's,

1; *n* = 0

*k*=0(2*k* + 1)−*q* = (−1)*n*(2*n* + 1)−*q*

*n* (*q*) =

q

2.5

2.6

2.7

4.5 4.75

*ζ*(*q*)(1 − 2−*q*) ; *q* ≥ 2 (Type I) (92)

0; otherwise (93)

*n* . They can be bracketed between (87) and *ψ*(*r*, *z*) =

(2*n* + 1)*q* ; *q* ≥ 2 (94)

(2*n* + 1)*q* = 1 (95)

2.1 2.4

(—) and 20 (−·−).

form of *αn*, namely,

*z* sin( 1

respectively.

*α*−

the negative sign in the superscript of *α*−

is needed. We begin by introducing

**3.5 Type II solutions with decreasing energy levels**

18%

47%

2.8

8% 8%

*n* (*q*) = (−1)*n*(2*n* + 1)−*q* ∑∞

is recovered. In fact, using (92), it can be rigorously shown that

lim *q*→∞ *α*−

*α*+

The remaining steps are similar. Substitution into (70) unravels

*n* (*q*) = *Bq*

∞ ∑ *n*=0

Type II

Type I

 L 10 20

3.2

3.6 3.8

•

#### **3.4 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. It may be hypothesized that two complementary families of solutions exist with the unique characteristics of exhibiting 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 (85) as a guide. Based on the form obtained through Lagrangian optimization, we note that

$$\alpha\_n = \frac{8(-1)^n}{\pi^2 (2n+1)^2} \sim \frac{(-1)^n A\_2}{(2n+1)^2} \tag{88}$$

where *A*2 = 8/*π*2 can be deduced from the radial inflow requirement given by (70). Its subscript is connected with the power of (2*n* + 1) in the denominator. To generalize, we posit the generic Type I form

$$\mathfrak{a}\_n^-(q) = \frac{(-1)^n A\_q}{(2n+1)^q}; \quad q \ge 2\tag{89}$$

where the exponent *q* will be referred to as the *kinetic energy power index*. This is due to its strong connection with the kinetic energy density as it will be shown shortly. The constant *Aq* can be used to make (89) consistent with (70). This enables us to retrieve

$$\sum\_{n=0}^{\infty} (-1)^n \frac{(-1)^n A\_q}{(2n+1)^q} = 1 \tag{90}$$

or

$$A\_q = \frac{1}{\sum\_{n=0}^{\infty} (2n+1)^{-q}} = \frac{1}{\zeta(q)(1-2^{-q})}; \quad \zeta(q) = \sum\_{k=1}^{\infty} k^{-q} \tag{91}$$

where *ζ*(*s*) is Riemann's zeta function. Clearly, the case corresponding to *q* = 2 reproduces the state of least energy expenditure. Furthermore, the *q* ≥ 2 condition is needed to ensure

$$\begin{array}{ccccc}\hline w(r,0) & & \psi^{-}(r,z) & & & w^{-}(r,z)\\\hline & & 8 & \raisebox{-0.0pt}{\displaystyle} & & & & & & & \\ \hline & & & & & & & & & & \\ \hline \end{array}$$

$$0 \qquad \psi\_{\text{ref}}^- \equiv \frac{8}{\pi^2} z \sum\_{n=0}^{\infty} \frac{(-1)^n}{(2n+1)^2} \sin \chi\_n \mapsto r^2 z \ w\_{\text{ref}}^- \equiv \frac{8}{\pi} z \sum\_{n=0}^{\infty} \frac{(-1)^n}{(2n+1)} \cos \chi\_n \mapsto 2z$$

$$\mathcal{W}\_{\rm c} \qquad \qquad \psi\_{\rm ref}^{-} + \frac{4\mathcal{W}\_{\rm c}}{\pi^{2}} \sum\_{n=0}^{\infty} \frac{(-1)^{n}}{(2n+1)^{2}} \sin \chi\_{n} \qquad \qquad w\_{\rm ref}^{-} + \frac{4\mathcal{W}\_{\rm c}}{\pi} \sum\_{n=0}^{\infty} \frac{(-1)^{n}}{(2n+1)} \cos \chi\_{n}$$

$$\begin{array}{cc} \left( \begin{array}{c} \left( \frac{1}{2} \pi r^2 \right) \end{array} \right) & \qquad \qquad \psi\_{\text{ref}}^{-} + \frac{\mathcal{W}\_{\text{C}}}{\pi} \sin(\frac{1}{2} \pi r^2) \end{array} \qquad \qquad \qquad \qquad \qquad w\_{\text{ref}}^{-} + \mathcal{W}\_{\text{C}} \cos(\frac{1}{2} \pi r^2) \end{array}$$

$$\mathcal{W}\_{\rm c}(1-r^{2}) \qquad \qquad \qquad \psi\_{\rm ref}^{-} + \frac{8W\_{\rm c}}{\pi^{3}} \sum\_{n=0}^{\infty} \frac{\sin \chi\_{n}}{(2n+1)^{3}} \qquad \qquad \qquad \qquad w\_{\rm ref}^{-} + \frac{8W\_{\rm c}}{\pi^{2}} \sum\_{n=0}^{\infty} \frac{\cos \chi\_{n}}{(2n+1)^{2}}$$

Table 3. Summary of least kinetic energy solutions.

24 Will-be-set-by-IN-TECH

At this juncture, we have identified only one profile bearing the minimum kinetic energy that the flow can possibly afford. It may be hypothesized that two complementary families of solutions exist with the unique characteristics of exhibiting 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 (85) as a guide. Based

*π*2(2*n* + 1)2 ∼ (−1)*nA*2

where *A*2 = 8/*π*2 can be deduced from the radial inflow requirement given by (70). Its subscript is connected with the power of (2*n* + 1) in the denominator. To generalize, we posit

where the exponent *q* will be referred to as the *kinetic energy power index*. This is due to its strong connection with the kinetic energy density as it will be shown shortly. The constant *Aq*

(−1)*n* (−1)*nAq*

where *ζ*(*s*) is Riemann's zeta function. Clearly, the case corresponding to *q* = 2 reproduces the state of least energy expenditure. Furthermore, the *q* ≥ 2 condition is needed to ensure

(2*n* + 1)2 sin *χn* �→ *r*2*z w*−

(2*n* + 1)2 sin *χn w*−

*ζ*(*q*)(1 − 2−*q*)

*n* (*q*) = (−1)*nAq*

(2*n* + 1)2 (88)

(2*n* + 1)*q* ; *q* ≥ 2 (89)

(2*n* + 1)*q* = 1 (90)

∞ ∑ *k*=1

(−1)*n*

∞ ∑ *n*=0

ref + *W*c cos( 1

∞ ∑ *n*=0

8*W*c *π*2

*k*−*q* (91)

(2*n* + 1) cos *χn* �→ 2*z*

2*πr*2)

cos *χn* (2*n* + 1)2

(−1)*n* (2*n* + 1) cos *χn*

; *ζ*(*q*) =

ref ≡ 8 *π z* ∞ ∑ *n*=0

ref +

2) *w*−

(2*n* + 1)3 *w*−

4*W*c *π*

ref +

**3.4 Type I solutions with increasing energy levels**

the generic Type I form

0 *ψ*−

*W*c *ψ*−

*W*c(1 − *r*2) *ψ*−

*W*c cos( 1

or

on the form obtained through Lagrangian optimization, we note that

*α*−

can be used to make (89) consistent with (70). This enables us to retrieve ∞ ∑ *n*=0

*n*=0(2*n* + 1)−*q* = 1

*w*(*r*, 0) *ψ*−(*r*, *z*) *w*−(*r*, *z*)

(−1)*n*

sin( 1 2*πr*

sin *χn*

∞ ∑ *n*=0

(−1)*n*

*Aq* = 1 ∑∞

> ∞ ∑ *n*=0

4*W*c *π*2

ref +

ref +

Table 3. Summary of least kinetic energy solutions.

∞ ∑ *n*=0

> *W*c *π*

8*W*c *π*3

ref ≡ 8 *π*2 *z*

2*πr*2) *ψ*−

ref +

*αn* = 8(−1)*n*

Fig. 12. Variation of the kinetic energy density with the energy power index for Type I (lower branch) and Type II (upper branch) solutions. These are shown at two aspect ratios, *L* = 10 (—) and 20 (−·−).

series convergence down to the vorticity. Backward substitution allows us to extract the final form of *αn*, namely,

$$\alpha\_n^-(q) = \frac{(-1)^n (2n+1)^{-q}}{\sum\_{k=0}^{\infty} (2k+1)^{-q}} = \frac{(-1)^n (2n+1)^{-q}}{\zeta(q) (1 - 2^{-q})}; \quad q \ge 2 \quad \text{(Type I)}\tag{92}$$

To understand the effect of the energy power index *q* on the kinetic energy density, we use (92) and (74) to plot E versus *q* at two aspect ratios. This plot corresponds to the lower branch of Figure 12 for both *L* = 10 and 20. Interestingly, as *q* → ∞, Taylor–Culick's classic solution is recovered. In fact, using (92), it can be rigorously shown that

$$\lim\_{q \to \infty} \mathfrak{a}\_n^-(q) = \begin{cases} 1; & n = 0 \\ 0; & \text{otherwise} \end{cases} \tag{93}$$