**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( <sup>1</sup> <sup>2</sup>*π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

U<sup>w</sup>

(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.

<sup>2</sup>*π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*<sup>w</sup>

y x

z

prove crucial in subsequent developments of this problem.

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

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

*<sup>r</sup>* <sup>=</sup> *<sup>r</sup>*<sup>∗</sup> *a*

*<sup>u</sup>* <sup>=</sup> *<sup>u</sup>*<sup>∗</sup> *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.

; *<sup>z</sup>* <sup>=</sup> *<sup>z</sup>*<sup>∗</sup> *a*

> ; *<sup>w</sup>* <sup>=</sup> *<sup>w</sup>*<sup>∗</sup> *U*w

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

> *u ∂u ∂r* + *w ∂u <sup>∂</sup><sup>z</sup>* <sup>=</sup> <sup>−</sup> *<sup>∂</sup><sup>p</sup>*

*u ∂w ∂r* + *w ∂w <sup>∂</sup><sup>z</sup>* <sup>=</sup> <sup>−</sup> *<sup>∂</sup><sup>p</sup>*

One may now invoke the dyadic vector identity **<sup>u</sup>** · ∇**<sup>u</sup>** ≡ ∇( <sup>1</sup>

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*<sup>c</sup> = const *W*<sup>c</sup> cos( <sup>1</sup>

*<sup>W</sup>*c(<sup>1</sup> <sup>−</sup> *<sup>r</sup>m*)

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

<sup>2</sup>*πr*2)

; <sup>∇</sup> <sup>=</sup> *<sup>a</sup>*∇∗; *<sup>p</sup>* <sup>=</sup> *<sup>p</sup>*<sup>∗</sup>

Internal Flows Driven by Wall-Normal Injection 101

; **<sup>Ω</sup>** <sup>=</sup> **<sup>Ω</sup>**∗*<sup>a</sup> 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*<sup>2</sup> w

; *<sup>W</sup>*<sup>c</sup> <sup>=</sup> *<sup>W</sup>*<sup>∗</sup>

; *<sup>ψ</sup>* <sup>=</sup> *<sup>ψ</sup>*<sup>∗</sup> *a*2*U*<sup>w</sup> ;

> ; *<sup>L</sup>* <sup>=</sup> *<sup>L</sup>*<sup>∗</sup> *a*

*<sup>∂</sup><sup>z</sup>* <sup>=</sup> <sup>0</sup> (3a)

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

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

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

*<sup>∂</sup><sup>r</sup>* (3b)

*<sup>∂</sup><sup>z</sup>* (3c)

<sup>2</sup>**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*∗ <sup>0</sup> (*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*∗ <sup>c</sup> = *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*<sup>w</sup> = −*u*∗(*a*, *z*∗) is commensurate with propellant or fuel regression rates. In hybrids, *U*w can be appreciably smaller than *W*∗ <sup>c</sup> 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 <sup>r</sup> <sup>a</sup>

\* \*

Wc

equations will be performed.

**2.2 Mathematical idealization**

<sup>0</sup> (*r*∗). This could be given by

<sup>0</sup> (*r*∗) = *w*∗(*r*∗, *z*<sup>∗</sup> = 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

<sup>c</sup> = const uniform

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

<sup>c</sup> = *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*<sup>w</sup> = −*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

<sup>2</sup>*πr*∗2/*a*2) cosine

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

<sup>c</sup> due to slow fuel pyrolysis; in SRM

(1)

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

⎧ ⎪⎪⎪⎨

*W*∗

*W*∗ <sup>c</sup> cos( <sup>1</sup>

*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:

*<sup>∂</sup><sup>z</sup>* <sup>=</sup> *<sup>α</sup>*¯ *<sup>A</sup>* cos( <sup>1</sup>

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

*<sup>∂</sup><sup>r</sup>* <sup>=</sup> *rC*(*α*¯ *<sup>z</sup>* <sup>+</sup> *<sup>β</sup>*¯) cos( <sup>1</sup>

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 *β<sup>n</sup>* = 0, we recover Culick's original solution

*ψ*(*r*, *z*) = *z* sin( <sup>1</sup>

∞ ∑ *n*=0

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

*α<sup>n</sup>* sin[(*n* + <sup>1</sup>

*∂ψ*(0, *z*)

*∂ ψ*(1, *z*)

For convenience, we introduce *<sup>χ</sup><sup>n</sup>* <sup>≡</sup> <sup>1</sup>

*∂ψ*(1, *z*) *<sup>∂</sup><sup>z</sup>* <sup>=</sup>

> 1 *r*

*βn* 1 0

*∂ψ*(*r*, 0) *<sup>∂</sup><sup>r</sup>* <sup>=</sup> *<sup>π</sup>*

<sup>2</sup>*C*) = 0. This is satisfied by

and so cos( <sup>1</sup>

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

<sup>2</sup>*Cr*2) + *<sup>α</sup>*¯ *<sup>B</sup>* sin( <sup>1</sup>

<sup>2</sup>*Cr*2) 

<sup>2</sup> )*π*] = 1 or

<sup>2</sup>*πr*

∞ ∑ *n*=0

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

*ψ*(*r*, *z*) = *z* sin( <sup>1</sup>

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

<sup>2</sup>*πr* <sup>2</sup>) +

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*<sup>n</sup>* <sup>+</sup> <sup>1</sup>) cos<sup>2</sup> *<sup>χ</sup><sup>n</sup> r dr* <sup>=</sup> <sup>1</sup>

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

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

∞ ∑ *n*=0 <sup>2</sup>*Cr*2) 

*<sup>r</sup>*=<sup>1</sup> <sup>=</sup> 0; <sup>∀</sup> *<sup>z</sup>* <sup>∈</sup> **<sup>R</sup>**<sup>+</sup>

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

<sup>2</sup> )*πr*

∞ ∑ *n*=0

*α*<sup>0</sup> = 1 and *α<sup>n</sup>* = 0; ∀*n* �= 0 (23)

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

(*αnz* + *βn*) sin *χ<sup>n</sup>* (21)

*<sup>r</sup>*=<sup>0</sup> <sup>=</sup> <sup>0</sup> (17)

<sup>2</sup>] (20)

(−1)*nα<sup>n</sup>* <sup>=</sup> <sup>1</sup> (22)

<sup>2</sup>) (24)

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

*w*0(*r*) cos *χ<sup>n</sup> r dr* (27)

(2*n* + 1)*β<sup>n</sup>* cos *χ<sup>n</sup>* = *w*0(*r*) (26)

<sup>0</sup> (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* 

<sup>=</sup> [*F*(*ψ*)]*<sup>z</sup>* [*F*(*ψ*)]*<sup>r</sup>*

*∂*2*ψ <sup>∂</sup>r*<sup>2</sup> <sup>−</sup> <sup>1</sup> *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*)

*∂*<sup>2</sup> *ψ*(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( <sup>1</sup>

; *<sup>w</sup>* <sup>=</sup> <sup>1</sup> *r ∂ψ*

*<sup>∂</sup><sup>r</sup>* (9)

(11)

(12)

**Ω** = Ω*θ***e***<sup>θ</sup>* ≡ Ω**e***<sup>θ</sup>* (10)

<sup>=</sup> *<sup>ψ</sup><sup>z</sup> ψr*

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

*<sup>∂</sup><sup>z</sup>* <sup>=</sup> <sup>0</sup> (14a)

*<sup>∂</sup><sup>r</sup>* <sup>=</sup> <sup>0</sup> (14b)

*<sup>∂</sup><sup>z</sup>* <sup>=</sup> <sup>1</sup> (14c)

*<sup>∂</sup><sup>r</sup>* <sup>=</sup> *<sup>w</sup>*0(*r*) (14d)

*<sup>∂</sup><sup>z</sup>* <sup>=</sup> <sup>0</sup> (15a)

*<sup>∂</sup>r∂<sup>z</sup>* <sup>=</sup> <sup>0</sup> (15b)

<sup>2</sup>*Cr*2)] (16)

<sup>2</sup>*Cr*2) + *<sup>B</sup>* sin( <sup>1</sup>

(Ω/*r*)*<sup>r</sup>*

<sup>=</sup> *<sup>ψ</sup><sup>z</sup> ψr*

<sup>=</sup> 0 or (Ω/*r*)*<sup>z</sup>*

<sup>=</sup> *<sup>F</sup>ψψ<sup>z</sup> Fψψ<sup>r</sup>*

*<sup>u</sup>* <sup>=</sup> <sup>−</sup><sup>1</sup> *r ∂ψ ∂z*

Its substitution into the vorticity transport equation (5) yields

(Ω/*r*)*<sup>z</sup>* (Ω/*r*)*<sup>r</sup>*

*∂*2*ψ <sup>∂</sup>z*<sup>2</sup> <sup>+</sup>

 Ω *r* <sup>−</sup> *∂ψ ∂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( <sup>1</sup> <sup>2</sup>*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 *<sup>χ</sup><sup>n</sup>* <sup>≡</sup> <sup>1</sup> <sup>2</sup> (2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)*πr*<sup>2</sup> 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 *β<sup>n</sup>* = 0, we recover Culick's original solution

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

Other forms of *α<sup>n</sup>* will be discussed in Section 3. At present, we let *α*<sup>0</sup> = 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

*<sup>β</sup><sup>n</sup>* <sup>=</sup> <sup>4</sup>(−1)*nW*<sup>c</sup>

<sup>2</sup>) + <sup>4</sup>*W*<sup>c</sup> *π*2

> <sup>2</sup>) + <sup>4</sup>*W*<sup>c</sup> *π*

The character of (34) is illustrated in Figure 4. Using *W*<sup>c</sup> = *U*<sup>w</sup> = 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*<sup>c</sup> = 1 corresponds to a simulated solid

> *<sup>π</sup>* <sup>≡</sup> *<sup>W</sup>*h; *<sup>n</sup>* <sup>=</sup> <sup>0</sup> 0; otherwise

> > <sup>2</sup>*πr*

<sup>2</sup>*πr* 2);

> <sup>2</sup>*πr* 2)

Ω(*r*, *z*) = *π*2*rz* sin( <sup>1</sup>

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

∞ ∑ *n*=0

> ∞ ∑ *n*=0

> > <sup>2</sup>*πr*

(−1)*<sup>n</sup>*

(−1)*<sup>n</sup>*

*<sup>π</sup>*2(2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)<sup>2</sup> (33)

(2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)<sup>2</sup> sin *<sup>χ</sup><sup>n</sup>* (34)

(2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>) cos *<sup>χ</sup><sup>n</sup>* (35)

<sup>2</sup>) (36)

<sup>2</sup>) (38)

(37)

(39)

2005), and others (uniform flow).

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

*ψ*(*r*, *z*) = *z* sin( <sup>1</sup>

The axial velocity and vorticity may be easily determined to be

*w*(*r*, *z*) = *πz* cos( <sup>1</sup>

**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

*β<sup>n</sup>* =

⎧ ⎨ ⎩ *W*c

*ψ*(*r*, *z*)=(*z* + *W*h) cos( <sup>1</sup>

*w*(*r*, *z*) = *π*(*z* + *W*h) cos( <sup>1</sup>

<sup>Ω</sup>(*r*, *<sup>z</sup>*) = *<sup>π</sup>*2*r*(*<sup>z</sup>* + *<sup>W</sup>*h) sin( <sup>1</sup>

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.

<sup>2</sup>*πr*

<sup>2</sup>*π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 *β<sup>n</sup>* through (28). By way of confirmation, the classical Taylor–Culick solution with inert headwall may be readily recovered by setting *β<sup>n</sup>* = 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)*π*

<sup>2</sup>*πr* 2);

<sup>2</sup>) + *π*

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

<sup>2</sup>) + *π*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 *β<sup>n</sup>* through (28). By way of confirmation, the classical Taylor–Culick solution with inert headwall may be readily recovered by setting *β<sup>n</sup>* = 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*<sup>c</sup> = const uniform

*<sup>W</sup>*c(<sup>1</sup> <sup>−</sup> *<sup>r</sup>*2) parabolic

<sup>2</sup>*πr*2) half cosine

∞ ∑ *n*=0

> ∞ ∑ *n*=0

<sup>2</sup>*πr*

<sup>2</sup>*πr*

0.3

� 1 0 0.4 0.5 0.6 0.7 0.8 0.9

*w*0(*r*) cos *χ<sup>n</sup> r dr* (28)

�

(2*n* + 1)*β<sup>n</sup>* cos *χ<sup>n</sup>* (30)

(2*n* + 1)2*β<sup>n</sup>* sin *χ<sup>n</sup>* (31)

sin *χn* (29)

(32)

*w*0(*r*) cos *χ<sup>n</sup> r dr*

0.1

0.1

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

(2*n* + 1)*π*

∞ ∑ *n*=0

*<sup>β</sup><sup>n</sup>* <sup>=</sup> <sup>4</sup>

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

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

*<sup>u</sup>*(*r*) = <sup>−</sup>*r*−<sup>1</sup> sin( <sup>1</sup>

*w*(*r*, *z*) = *πz* cos( <sup>1</sup>

Ω(*r*, *z*) = *π*2*rz* sin( <sup>1</sup>

*w*0(*r*) =

⎧ ⎪⎪⎪⎨

⎪⎪⎪⎩

*W*<sup>c</sup> cos( <sup>1</sup>

this historical benchmark are illustrated in Figure 3.

**2.5 Axisymmetric headwall injection profiles**

injection pattern at *z* = 0, namely,

<sup>2</sup>*πr* <sup>2</sup>) +

*ψ*(*r*, *z*) = *z* sin( <sup>1</sup>

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*<sup>c</sup> = *U*<sup>w</sup> = 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*<sup>c</sup> = 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*<sup>c</sup> = 1.

<sup>2</sup>*πr*

<sup>2</sup>*πr*

<sup>2</sup>*πr*

*<sup>Q</sup>*(*r*, *<sup>z</sup>*) = �∇ × **<sup>u</sup>** <sup>×</sup> **<sup>Ω</sup>**� <sup>=</sup> <sup>−</sup> *<sup>∂</sup>*

*r*2 *∂ψ ∂z* + 1 *r ∂ψ ∂z ∂*Ω *<sup>∂</sup><sup>r</sup>* <sup>−</sup> <sup>1</sup> *r ∂ψ ∂r ∂*Ω

*ψ*(*r*, *z*) = *z* sin( <sup>1</sup>

*w*(*r*, *z*) = *πz* cos( <sup>1</sup>

Ω(*r*, *z*) = *π*2*rz* sin( <sup>1</sup>

straightforward to see that *Q* may be calculated from

In terms of the streamfunction and the vorticity, we have

*<sup>n</sup>ψ<sup>n</sup> r*

*∂ψn ∂z* + 1 *r ∂ψn ∂z ∂ ∂r* (*C*<sup>2</sup>

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

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

*<sup>Q</sup>*(*r*, *<sup>z</sup>*) = <sup>−</sup> <sup>Ω</sup>

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

0.13

0

**2.6 Nonlinear residual error**

residual. Using Ω = Ω*<sup>n</sup>* = *C*<sup>2</sup>

*<sup>Q</sup>* <sup>=</sup> <sup>−</sup> <sup>1</sup> *r*2

*Qn* <sup>=</sup> <sup>−</sup> *<sup>C</sup>*<sup>2</sup>

<sup>=</sup> <sup>−</sup> *<sup>C</sup>*<sup>2</sup> *<sup>n</sup>ψ<sup>n</sup> 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:

<sup>2</sup>) + <sup>8</sup>*W*<sup>c</sup> *π*3

> <sup>2</sup>) + <sup>8</sup>*W*<sup>c</sup> *π*2

> > <sup>2</sup>) + <sup>8</sup>*W*<sup>c</sup> *π 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

*<sup>n</sup>ψ<sup>n</sup>* + *<sup>C</sup>*<sup>2</sup> *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 <sup>∂</sup><sup>r</sup>* <sup>−</sup> <sup>1</sup> *r*

∞ ∑ *n*=0

*∂ψn ∂r*

∞ ∑ *n*=0

*∂*Ω*n*

*<sup>∂</sup><sup>z</sup>* (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*) <sup>−</sup> <sup>1</sup> *r ∂ψn ∂r*

(*u*Ω) <sup>−</sup> *<sup>∂</sup>*

*∂z*

*∂ ∂z* (*C*<sup>2</sup> *nrψn*)

*∂ψn <sup>∂</sup><sup>r</sup>* <sup>−</sup> *<sup>C</sup>*<sup>2</sup> *n ∂ψn ∂r*

1

1

1

0.50 0.63 0.75 0.88 1.0

1.1

(2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)<sup>3</sup> sin *<sup>χ</sup><sup>n</sup>* (41)

(2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)<sup>2</sup> cos *<sup>χ</sup><sup>n</sup>* (42)

(2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>) sin *<sup>χ</sup><sup>n</sup>* (43)

(*w*Ω) (44)

*<sup>∂</sup><sup>z</sup>* (45)

*∂ψn <sup>∂</sup><sup>z</sup>* <sup>=</sup> <sup>0</sup> (46)

Fig. 4. Streamlines corresponding to uniform headwall injection with *W*<sup>c</sup> = 1.

Fig. 5. Streamlines corresponding to cosine headwall injection with *W*<sup>c</sup> = 1.

#### **2.5.3 Parabolic injection**

For the parabolic, laminar-like profile, one may substitute *<sup>w</sup>*0(*r*) = *<sup>W</sup>*c(<sup>1</sup> <sup>−</sup> *<sup>r</sup>*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 *<sup>w</sup>*0(*r*) = *<sup>W</sup>*c(<sup>1</sup> <sup>−</sup> *<sup>r</sup>*2) into (28) and

*<sup>β</sup><sup>n</sup>* <sup>=</sup> <sup>8</sup>*W*<sup>c</sup>

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*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)3*π*<sup>3</sup> (40)

Fig. 4. Streamlines corresponding to uniform headwall injection with *W*<sup>c</sup> = 1.

0.13

Fig. 5. Streamlines corresponding to cosine headwall injection with *W*<sup>c</sup> = 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*<sup>c</sup> = 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 Ω = Ω*<sup>n</sup>* = *C*<sup>2</sup> *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

<sup>2</sup>**u** · **u** −

where *p*<sup>0</sup> = *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 <sup>∂</sup>r∂<sup>z</sup>* <sup>−</sup> *<sup>u</sup> 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

<sup>2</sup>*π*2*z*<sup>2</sup> <sup>−</sup> <sup>1</sup>

<sup>2</sup>*π*2*z*<sup>2</sup> <sup>−</sup> *<sup>W</sup>*c*π<sup>z</sup>* <sup>−</sup> <sup>1</sup>

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,

<sup>2</sup> *<sup>r</sup>*−<sup>2</sup> sin2( <sup>1</sup>

<sup>2</sup>)**e***<sup>r</sup>* + *<sup>π</sup>*(*<sup>z</sup>* + <sup>W</sup>*h*) cos( <sup>1</sup>

<sup>2</sup>*πr*

<sup>2</sup> *<sup>r</sup>*−<sup>2</sup> sin2( <sup>1</sup>

<sup>2</sup>*πr*

*∂*<sup>2</sup> *p <sup>∂</sup>r∂<sup>z</sup>* <sup>=</sup> *<sup>∂</sup>*<sup>2</sup> *<sup>p</sup> ∂z∂r*

 *u ∂w ∂r*

*dz* (57)

*<sup>∂</sup><sup>r</sup>* <sup>=</sup> <sup>0</sup> (59)

<sup>2</sup>) (60)

<sup>2</sup>)**e***<sup>z</sup>* (61)

<sup>2</sup>*πr*2) (62)

(58)

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

*<sup>p</sup>* <sup>=</sup> *<sup>p</sup>*<sup>0</sup> <sup>−</sup> <sup>1</sup>

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

*u ∂*2*w <sup>∂</sup>r*<sup>2</sup> <sup>+</sup> *<sup>w</sup>*

*<sup>p</sup>*(*r*, *<sup>z</sup>*) = *<sup>p</sup>*<sup>0</sup> <sup>−</sup> <sup>1</sup>

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

<sup>2</sup>*πr*

**<sup>u</sup>** <sup>=</sup> <sup>−</sup>*r*−<sup>1</sup> sin( <sup>1</sup>

*<sup>p</sup>*(*r*, *<sup>z</sup>*) = *<sup>p</sup>*<sup>0</sup> <sup>−</sup> <sup>1</sup>

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}$$

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

Furthermore, for the Taylor–Culick class of solutions, *α*<sup>0</sup> = 1 and *α<sup>n</sup>* = 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 (*β<sup>n</sup>* = 0, ∀*n*) or sinusoidal headwall injection profiles (*β<sup>n</sup>* = 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*<sup>2</sup> *<sup>n</sup>* <sup>−</sup> *<sup>C</sup>*<sup>2</sup> <sup>0</sup> = 0. In the case of no headwall injection, *β<sup>n</sup>* = *Dn* = 0, thus leading to an exact representation. As *Dn* → 0, the solutions become more accurate. Generally, *<sup>β</sup><sup>n</sup>* �<sup>=</sup> 0 and so *Dn* will only vanish when *<sup>C</sup>*<sup>2</sup> *<sup>n</sup>* = *C*<sup>2</sup> <sup>0</sup>. 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, *α*<sup>0</sup> = 1 and *α<sup>n</sup>* = 0, ∀ *n* �= 0. This leaves

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

*∂*Ω*n <sup>∂</sup><sup>z</sup>* <sup>=</sup> *<sup>∂</sup>*Ω<sup>0</sup>

∞ ∑ *n*=0

*<sup>n</sup>ψ<sup>n</sup>* + *<sup>C</sup>*<sup>2</sup> *nr ∂ψn*

> <sup>2</sup>*πr* <sup>2</sup>)*r* ∞ ∑ *n*=1

*Dn* cos � 1

*<sup>∂</sup><sup>r</sup>* <sup>=</sup> sin( <sup>1</sup>

<sup>2</sup>*πr* 2) ∞ ∑ *n*=1

*<sup>n</sup>* <sup>−</sup> *<sup>C</sup>*<sup>2</sup>

*Dn*,parabolic <sup>∼</sup> *<sup>n</sup>*(*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)

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

*Dn*,uniform <sup>∼</sup> (−1)*<sup>n</sup> <sup>n</sup>*(*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)

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*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)<sup>2</sup> −−−→ *<sup>n</sup>*→<sup>∞</sup>

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

*∂*Ω*n <sup>∂</sup><sup>r</sup>* <sup>−</sup> *<sup>C</sup>*<sup>2</sup> 0 ∞ ∑ *n*=0

*<sup>∂</sup><sup>z</sup>* <sup>=</sup> *<sup>α</sup><sup>n</sup>* sin *<sup>χ</sup><sup>n</sup>*

*<sup>∂</sup><sup>z</sup>* <sup>=</sup> *<sup>C</sup>*<sup>2</sup> 0*r ∂ψ*0

*<sup>n</sup>α<sup>n</sup>* sin *χ<sup>n</sup>*

*∂ψn ∂r*

*Cnβn*(*C*<sup>2</sup>

<sup>2</sup> (2*n* + 1)*πr*

<sup>0</sup> = 0. In the case of no headwall injection, *β<sup>n</sup>* = *Dn* = 0,

1

*<sup>n</sup>* = *C*<sup>2</sup>

�

*<sup>∂</sup><sup>r</sup>* (51)

*<sup>n</sup>* <sup>−</sup> *<sup>C</sup>*<sup>2</sup>

2 �

<sup>0</sup>) ≡ *n*(*n* + 1)(2*n* + 1)*β<sup>n</sup>* (54)

(2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>) −−−→ *<sup>n</sup>*→<sup>∞</sup> <sup>±</sup><sup>∞</sup> (56)

<sup>0</sup>. To illustrate this character, we

<sup>4</sup> (55)

*<sup>∂</sup><sup>z</sup>* (49)

<sup>0</sup>) cos *χ<sup>n</sup>* (52)

(48)

(50)

(53)

*∂*Ω*n <sup>∂</sup><sup>z</sup>* <sup>=</sup> *rC*<sup>2</sup>

*<sup>ψ</sup><sup>n</sup>* = (*αnz* <sup>+</sup> *<sup>β</sup>n*) sin *<sup>χ</sup>n*; *∂ψ<sup>n</sup>*

*<sup>∂</sup><sup>r</sup>* <sup>=</sup> *rCn*(*αnz* <sup>+</sup> *<sup>β</sup>n*) cos *<sup>χ</sup>n*;

<sup>2</sup>*πr* 2);

Ω*<sup>n</sup>* + 1 *r*

∞ ∑ *n*=0

*∂*Ω*n <sup>∂</sup><sup>r</sup>* <sup>=</sup> *<sup>C</sup>*<sup>2</sup>

*<sup>∂</sup><sup>z</sup>* <sup>=</sup> sin( <sup>1</sup>

� − 1 *r*2

where

us with

Finally, noting that

we expand (52) into

cosine profile for which *C*<sup>2</sup>

*<sup>Q</sup>*(*r*) = *∂ψ*<sup>0</sup> *∂z*

we retrieve

where

⎧ ⎪⎨

⎪⎩

*∂ψn <sup>∂</sup><sup>z</sup>* <sup>=</sup> *∂ψ*<sup>0</sup>

*<sup>Q</sup>* <sup>=</sup> *∂ψ*<sup>0</sup> *∂z*

> ∞ ∑ *n*=0

� *C*2 *<sup>n</sup>* <sup>−</sup> *<sup>C</sup>*<sup>2</sup> 0 � *∂ψn*

*Q*(*r*) = 4*π*3*r* sin( <sup>1</sup>

<sup>4</sup>*π*<sup>3</sup> *<sup>β</sup>n*(*C*<sup>2</sup>

parabolic injection, we find a quickly converging sequence, specifically

*Dn* <sup>≡</sup> *Cn*

*<sup>n</sup>* <sup>−</sup> *<sup>C</sup>*<sup>2</sup>

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

*∂ψ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*<sup>0</sup> = *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*−<sup>1</sup> sin( <sup>1</sup>

<sup>2</sup>*πr*

*w*uniform(0, *z*) = *w*parabolic(0, *z*) = *πz* + *W*<sup>c</sup> (64)

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

*nz*<sup>2</sup> <sup>−</sup> (2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)2*π*2*α<sup>n</sup> <sup>β</sup>nz* <sup>−</sup> <sup>1</sup>

<sup>2</sup>*π*2*z*<sup>2</sup> <sup>−</sup> *<sup>β</sup>*0*π*2*<sup>z</sup>* <sup>−</sup> <sup>1</sup>

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

<sup>2</sup>) = 0 (63)

<sup>2</sup>*π*2*z*<sup>2</sup> <sup>−</sup> *<sup>π</sup>zW*<sup>c</sup> (65)

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

2 *α*2

<sup>2</sup>*πr*

<sup>2</sup> *<sup>r</sup>*−<sup>2</sup> sin2( <sup>1</sup>

*nr*<sup>−</sup><sup>2</sup> sin2 *<sup>χ</sup><sup>n</sup>* (68)

<sup>2</sup>) (69)

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

*r*→0

*<sup>p</sup>*(0, *<sup>z</sup>*) = *<sup>p</sup>*<sup>0</sup> <sup>−</sup> <sup>1</sup>

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

*un* <sup>=</sup> <sup>−</sup>*αnr*−<sup>1</sup> sin *<sup>χ</sup>n*;

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

<sup>2</sup> (2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)2*π*2*α*<sup>2</sup>

*pn* <sup>=</sup> *<sup>p</sup>*<sup>0</sup> <sup>−</sup> <sup>1</sup>

∞ ∑ *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*, *<sup>z</sup>*) = *<sup>p</sup>*<sup>0</sup> <sup>−</sup> <sup>1</sup>

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 *<sup>n</sup>* = <sup>1</sup> <sup>2</sup> (*u*<sup>2</sup> *<sup>n</sup>* + *<sup>v</sup>*<sup>2</sup>

r

*En*(*r*, *θ*, *z*) = <sup>1</sup>

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

r

0

0

0

20

40

60

80 w <sup>W</sup><sup>c</sup>

20

40

60 w <sup>W</sup><sup>c</sup>

20

40

60 w <sup>W</sup><sup>c</sup>

Internal Flows Driven by Wall-Normal Injection 113

= 10

= 10

= 10

0 1

0 1

0 1

*<sup>n</sup>*) (71)

<sup>2</sup> )*πr*<sup>2</sup> (72)

r

r

r

(b) cosine

(d) parabolic

(f) uniform

*<sup>n</sup>* + *<sup>w</sup>*<sup>2</sup>

*<sup>χ</sup><sup>n</sup>* <sup>≡</sup> (*<sup>n</sup>* <sup>+</sup> <sup>1</sup>

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 <sup>W</sup><sup>c</sup>

10

20

30

40

50 w <sup>W</sup><sup>c</sup>

10

20

30

40

50

w <sup>W</sup><sup>c</sup>

= 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*<sup>c</sup> = 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*<sup>c</sup> = 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α<sup>n</sup>* <sup>=</sup> <sup>1</sup> (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 = <sup>1</sup> 6*π*3*L*<sup>3</sup>

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

*<sup>α</sup><sup>n</sup>* <sup>=</sup> (−1)*<sup>n</sup>*

(−1)*nα<sup>n</sup>* <sup>=</sup> <sup>1</sup>

*∂*G *∂λ* <sup>=</sup>

variables to obtain

**3.2 Critical length**

be accomplished by setting

and

G(*α*0, *<sup>α</sup>*1, *<sup>α</sup>*2,..., *<sup>λ</sup>*) = *<sup>E</sup>*<sup>V</sup> + *<sup>λ</sup>*

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*<sup>2</sup> (*an* + *π*−2*L*−2*dn*)−<sup>1</sup>

1 (*an* <sup>+</sup> *<sup>π</sup>*−2*L*−2*dn*) <sup>=</sup> *<sup>N</sup>*

∞ ∑ *i*=0

*αnan* +

∞ ∑ *n*=0

*<sup>α</sup><sup>n</sup>* <sup>=</sup> <sup>−</sup> <sup>6</sup>(−1)*n<sup>λ</sup>*

*<sup>λ</sup>* <sup>=</sup> <sup>−</sup> *<sup>π</sup>*3*L*<sup>3</sup>

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

(*an* <sup>+</sup> *<sup>π</sup>*−2*L*−2*dn*)*<sup>N</sup>* ; *<sup>N</sup>* <sup>=</sup>

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

∞ ∑ *n*=0

Some values of *α<sup>n</sup>* 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*<sup>V</sup> is completely determined.

Equation (74) can be normalized by *L*<sup>3</sup> 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<sup>∞</sup> = 2*π*/3. Granted this behavior, a critical aspect ratio *L*cr may be defined beyond

*N*

6 ∞ ∑ *n*=0 (−1)*nα<sup>n</sup>* <sup>−</sup> <sup>1</sup>

∇G(*αn*, *λ*) = 0; *n* ∈ **N**<sup>0</sup> (76)

+ (−1)*n<sup>λ</sup>* <sup>=</sup> <sup>0</sup> (77)

(−1)*nα<sup>n</sup>* <sup>−</sup> <sup>1</sup> <sup>=</sup> <sup>0</sup> (78)

*<sup>π</sup>*3*L*3(*an* <sup>+</sup> *<sup>π</sup>*−2*L*−2*dn*) (79)

1 *ai* + *π*−2*L*−<sup>2</sup>*di*

<sup>E</sup> <sup>=</sup> *<sup>E</sup>*<sup>V</sup> /*L*<sup>3</sup> (83)

(75)

(80)

(81)

*<sup>N</sup>* <sup>=</sup> <sup>1</sup> (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* <sup>−</sup><sup>1</sup> 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 *χ<sup>n</sup>*

*<sup>n</sup>π*−2*L*−2*dn*);

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

*nr*<sup>−</sup><sup>2</sup> sin2 *<sup>χ</sup><sup>n</sup>* + *<sup>π</sup>*2*α*<sup>2</sup>

50

100

150

<sup>200</sup> <sup>W</sup><sup>c</sup>

= 10

50

100

150

<sup>200</sup> <sup>W</sup><sup>c</sup>

= 10

0 1

0 1

*nz*<sup>2</sup>(2*<sup>n</sup>* + <sup>1</sup>)<sup>2</sup> cos<sup>2</sup> *<sup>χ</sup><sup>n</sup>*

*dn* <sup>=</sup> 3Cin[(2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)*π*] (74)

*<sup>r</sup>*<sup>2</sup> <sup>+</sup> *<sup>π</sup>*2*z*2(2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)<sup>2</sup> cos2 *<sup>χ</sup><sup>n</sup>*

*an* = (2*n* + 1)<sup>2</sup>

<sup>−</sup><sup>1</sup> 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*) = <sup>1</sup>

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

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

 *L* 0

∞ ∑ *n*=0

*nan* + *<sup>α</sup>*<sup>2</sup>

(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*<sup>3</sup> <sup>∞</sup> ∑ *n*=0 (*α*<sup>2</sup>

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

0.3, 0.5, 0.7, and 0.9.

*<sup>E</sup>*<sup>V</sup> =

where Cin(*x*) ≡

V

*<sup>E</sup>*<sup>V</sup> <sup>=</sup> <sup>1</sup>

 *x* 0

50

100

<sup>150</sup> <sup>W</sup><sup>c</sup>

= 1

50

100

<sup>150</sup> <sup>W</sup><sup>c</sup>

= 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 *α<sup>n</sup>* 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*<sup>V</sup> is completely determined.

#### **3.2 Critical length**

Equation (74) can be normalized by *L*<sup>3</sup> 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<sup>∞</sup> = 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)*<sup>n</sup>*

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*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)<sup>2</sup> sin *<sup>χ</sup><sup>n</sup>* �→ *<sup>r</sup>*

*<sup>n</sup>*=0(−1)*nα<sup>n</sup>*

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 α<sup>m</sup>* ∑*<sup>m</sup>*

Table 2. Convergence of the sidewall injection sequence *α<sup>n</sup>* 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

*<sup>ψ</sup>*(*r*, *<sup>z</sup>*) = <sup>8</sup>

0.333

0.500 0.667

0.833

0.857

<sup>2</sup>*z* (87)


Table 1. Convergence of the sidewall injection sequence *α<sup>n</sup>* 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

−<sup>1</sup>

1

 *an* ∞ ∑ *i*=0

(−1)*nα<sup>n</sup>* <sup>=</sup> <sup>8</sup>

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

: − ≤ <sup>∞</sup> <sup>∞</sup> L •• • cr cr 0.01

2 <sup>3</sup> • <sup>∞</sup> <sup>=</sup>

Ecr − E<sup>∞</sup> ≤ 0.01 E<sup>∞</sup> (84)

<sup>=</sup> <sup>8</sup>(−1)*<sup>n</sup>*

*<sup>π</sup>*2(2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)<sup>2</sup> (85)

(2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)<sup>2</sup> <sup>=</sup> <sup>1</sup> (86)

Table 1. Convergence of the sidewall injection sequence *α<sup>n</sup>* 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**

*<sup>L</sup>*→<sup>∞</sup> *<sup>α</sup><sup>n</sup>* = (−1)*<sup>n</sup>*

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 *α<sup>n</sup>* 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)*<sup>n</sup>* multiplier that appears in (89).

(−1)*nBq*

<sup>2</sup>*πr*2). In practice, profiles with *<sup>q</sup>* <sup>≥</sup> 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

*<sup>k</sup>*=0(2*<sup>k</sup>* <sup>+</sup> <sup>1</sup>)−*<sup>q</sup>* <sup>=</sup> (−1)*n*(2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)−*<sup>q</sup>*

*<sup>n</sup>* (*q*) =

q

2.5

2.6

2.7

4.5 4.75

*<sup>ζ</sup>*(*q*)(<sup>1</sup> <sup>−</sup> <sup>2</sup>−*q*) ; *<sup>q</sup>* <sup>≥</sup> 2 (Type I) (92)

0; otherwise (93)

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

(2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)*<sup>q</sup>* ; *<sup>q</sup>* <sup>≥</sup> <sup>2</sup> (94)

(2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)*<sup>q</sup>* <sup>=</sup> <sup>1</sup> (95)

2.1 2.4

(—) and 20 (−·−).

form of *αn*, namely,

*z* sin( <sup>1</sup>

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%

*<sup>n</sup>* (*q*) = (−1)*n*(2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)−*<sup>q</sup>* ∑<sup>∞</sup>

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

lim *<sup>q</sup>*→<sup>∞</sup> *<sup>α</sup>*<sup>−</sup>

*α*+

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

*<sup>n</sup>* (*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*<sup>2</sup> = 8/*π*<sup>2</sup> 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

*<sup>π</sup>*2(2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)<sup>2</sup> <sup>∼</sup> (−1)*nA*<sup>2</sup>

where *A*<sup>2</sup> = 8/*π*<sup>2</sup> 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)*<sup>n</sup>* (−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*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)<sup>2</sup> sin *<sup>χ</sup><sup>n</sup>* �→ *<sup>r</sup>*2*z w*<sup>−</sup>

(2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)<sup>2</sup> sin *<sup>χ</sup><sup>n</sup> <sup>w</sup>*<sup>−</sup>

*<sup>ζ</sup>*(*q*)(<sup>1</sup> <sup>−</sup> <sup>2</sup>−*q*)

*<sup>n</sup>* (*q*) = (−1)*nAq*

(2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)<sup>2</sup> (88)

(2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)*<sup>q</sup>* ; *<sup>q</sup>* <sup>≥</sup> <sup>2</sup> (89)

(2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)*<sup>q</sup>* <sup>=</sup> <sup>1</sup> (90)

∞ ∑ *k*=1

(−1)*<sup>n</sup>*

∞ ∑ *n*=0

ref <sup>+</sup> *<sup>W</sup>*<sup>c</sup> cos( <sup>1</sup>

∞ ∑ *n*=0

8*W*c *π*2

*k*−*<sup>q</sup>* (91)

(2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>) cos *<sup>χ</sup><sup>n</sup>* �→ <sup>2</sup>*<sup>z</sup>*

<sup>2</sup>*πr*2)

cos *χn* (2*n* + 1)<sup>2</sup>

(−1)*<sup>n</sup>* (2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>) cos *<sup>χ</sup><sup>n</sup>*

; *ζ*(*q*) =

ref <sup>≡</sup> <sup>8</sup> *π z* ∞ ∑ *n*=0

ref +

<sup>2</sup>) *w*<sup>−</sup>

(2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)<sup>3</sup> *<sup>w</sup>*<sup>−</sup>

4*W*c *π*

ref +

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

the generic Type I form

0 *ψ*−

*W*<sup>c</sup> *ψ*−

*<sup>W</sup>*c(<sup>1</sup> <sup>−</sup> *<sup>r</sup>*2) *<sup>ψ</sup>*<sup>−</sup>

*W*<sup>c</sup> cos( <sup>1</sup>

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

*<sup>n</sup>*=0(2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)−*<sup>q</sup>* <sup>=</sup> <sup>1</sup>

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

(−1)*<sup>n</sup>*

sin( <sup>1</sup> <sup>2</sup>*πr*

sin *χn*

∞ ∑ *n*=0

(−1)*<sup>n</sup>*

*Aq* <sup>=</sup> <sup>1</sup> ∑<sup>∞</sup>

> ∞ ∑ *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 <sup>≡</sup> <sup>8</sup> *π*2 *z*

<sup>2</sup>*πr*2) *<sup>ψ</sup>*<sup>−</sup>

ref +

*<sup>α</sup><sup>n</sup>* <sup>=</sup> <sup>8</sup>(−1)*<sup>n</sup>*

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}$$

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 the negative sign in the superscript of *α*− *<sup>n</sup>* . They can be bracketed between (87) and *ψ*(*r*, *z*) = *z* sin( <sup>1</sup> <sup>2</sup>*πr*2). In practice, profiles with *<sup>q</sup>* <sup>≥</sup> 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, respectively.

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

To capture solutions with energies that exceed that of Taylor-Culick's, a modified form of *αn* is needed. We begin by introducing

$$\alpha\_n^+(q) = \frac{B\_q}{(2n+1)^q}; \quad q \ge 2\tag{94}$$

The key difference here stands in the exclusion of the (−1)*<sup>n</sup>* multiplier that appears in (89). The remaining steps are similar. Substitution into (70) unravels

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

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

sin *χn*

Internal Flows Driven by Wall-Normal Injection 121

(−1)*<sup>n</sup>*

sin( <sup>1</sup> <sup>2</sup>*πr*

∞ ∑ *n*=0

(2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)<sup>2</sup> *<sup>w</sup>*<sup>+</sup>

(2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)<sup>3</sup> *<sup>w</sup>*<sup>+</sup>

0

Fig. 14. Effect of the kinetic energy power index on the two types of energy-based solutions.

Results are shown for (a) turn angle, (b) radial velocity, and (c) axial velocity.

4

8

10

w/z

(2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)<sup>2</sup> sin *<sup>χ</sup><sup>n</sup> <sup>w</sup>*<sup>+</sup>

sin *χn*

Table 4. Summary of solutions with most kinetic energy for various headwall injection

ref <sup>≡</sup> *<sup>π</sup>* C *z*

> 4*W*c *π*

ref +

ref +

<sup>2</sup>) *w*<sup>+</sup>

∞ ∑ *n*=0

∞ ∑ *n*=0

ref <sup>+</sup> *<sup>W</sup>*<sup>c</sup> cos( <sup>1</sup>

∞ ∑ *n*=0

8*W*c *π*2

most kinetic energy state

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

<sup>2</sup>*πr*2)

cos *χn* (2*n* + 1)<sup>2</sup>

(−1)*<sup>n</sup>* (2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>) cos *<sup>χ</sup><sup>n</sup>*

0 1

r

4 4 3

Type II

2

Type I least kinetic energy state Hart-McClure

(c) Axial velocity

3

q = 2

Taylor-Culick

∞ ∑ *n*=0

∞ ∑ *n*=0

> *W*c *π*

8*W*c *π*3

ref <sup>≡</sup> *<sup>z</sup>* C

> 4*W*c *π*2

ref +

ref +

<sup>2</sup> (2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)*πr*2.

<sup>q</sup> = 2 <sup>3</sup>

Taylor-Culick

<sup>2</sup> <sup>3</sup>

Type I

ref +

<sup>2</sup>*πr*2) *<sup>ψ</sup>*<sup>+</sup>

0 1

(a) Turn angle

4 3

Type I

0 1

(b) Radial velocity

r

Type II

Taylor-Culick

0 *ψ*<sup>+</sup>

*W*<sup>c</sup> *ψ*<sup>+</sup>

*<sup>W</sup>*c(<sup>1</sup> <sup>−</sup> *<sup>r</sup>*2) *<sup>ψ</sup>*<sup>+</sup>

*W*<sup>c</sup> cos( <sup>1</sup>

patterns. Here, *<sup>χ</sup><sup>n</sup>* <sup>≡</sup> <sup>1</sup>

0


u

r

0 q = 2

2

Type II

<sup>3</sup> <sup>4</sup>

u

30

60

90

Fig. 13. Comparison of the Taylor–Culick streamlines (—) and the Type II energy-maximized solution (*q* = 2) with stretched streamline curvature (−·−). Results are shown for an inert headwall.

or

$$B\_q = \frac{1}{\sum\_{n=0}^{\infty} (-1)^n (2n+1)^{-q}} = \frac{4^q}{\zeta(q, \frac{1}{4}) - \zeta(q, \frac{3}{4})};\tag{96}$$

where *ζ*(*q*, *a*) is the generalized Riemann zeta function given here as

$$\mathcal{Z}(q, a) = \sum\_{k=0}^{\infty} (k + a)^{-q}; \quad \forall a \in \mathbb{R} \tag{97}$$

Equation (94) yields the general structure of the Type II complementary family of solutions

$$\alpha\_n^+(q) = \frac{(2n+1)^{-q}}{\sum\_{k=0}^{\infty} (-1)^k (2k+1)^{-q}}$$

$$= \frac{4^q (2n+1)^{-q}}{\zeta(q, \frac{1}{4}) - \zeta(q, \frac{3}{4})}; \quad q \ge 2 \quad \text{(Type II)}\tag{98}$$

Note that the Type II solutions emerging from (98) dispose of kinetic energies that are higher than Taylor-Culick's. The variation of the solution with respect to *q* is embodied in the upper branch of Figure 12. According to this form of *α*<sup>+</sup> *<sup>n</sup>* , Taylor-Culick's model is recoverable asymptotically by taking the limit as *q* → ∞. Here too, most of the solutions exhibit energies that fall within one percent of Taylor-Culick's. The most interesting solutions are those corresponding to *q* = 2, 3, and 4 with energies that are 47.0, 8.08, and 2.4 percent larger than Taylor-Culick's. When the energy level is fixed at *q* = 2, a simplification follows for the Type II representation. Catalan's constant emerges in (98), namely,

$$\theta^{\ell} = \sum\_{k=0}^{\infty} (-1)^{k} (2k+1)^{-2} \simeq 0.915966\tag{99}$$

The Type II solution that carries the most energy at *q* = 2 is plotted in Figure 13 and listed in Table 4. In Figure 13, the Type II approximation is seen to overshoot the Taylor-Culick streamline curvature. In view of the two types of solutions with energies that either lag or surpass that of Taylor-Culick's, one may perceive the *q* → ∞ case as a saddle point to which other possible forms will quickly converge when their energies are shifted. Later in Section 3.12, we will use the entropy maximization principle to establish the Taylor–Culick model as a local equilibrium solution to which all other profiles will be attracted to.

26 Will-be-set-by-IN-TECH

0.250 0.375 0.500 0.625

*ζ*(*q*, <sup>1</sup>

<sup>4</sup> ) <sup>−</sup> *<sup>ζ</sup>*(*q*, <sup>3</sup>

4 )

(*<sup>k</sup>* <sup>+</sup> *<sup>a</sup>*)−*q*; <sup>∀</sup> *<sup>a</sup>* <sup>∈</sup> **<sup>R</sup>** (97)

; *q* ≥ 2 (Type II) (98)

*<sup>n</sup>* , Taylor-Culick's model is recoverable

(−1)*k*(2*<sup>k</sup>* <sup>+</sup> <sup>1</sup>)−<sup>2</sup> � 0.915966 (99)

0.750

0.875

L

; (96)

0.125

0.010

*Bq* <sup>=</sup> <sup>1</sup> ∑<sup>∞</sup>

*α*+

branch of Figure 12. According to this form of *α*<sup>+</sup>

Type II representation. Catalan's constant emerges in (98), namely,

∞ ∑ *k*=0

a local equilibrium solution to which all other profiles will be attracted to.

C =

where *ζ*(*q*, *a*) is the generalized Riemann zeta function given here as

*ζ*(*q*, *a*) =

*<sup>n</sup>* (*q*) = (2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)−*<sup>q</sup>* ∑<sup>∞</sup>

> <sup>=</sup> <sup>4</sup>*q*(2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)−*<sup>q</sup> ζ*(*q*, <sup>1</sup>

∞ ∑ *k*=0

0

headwall.

or

1

r

0 0.2L 0.4L 0.6L 0.8L

z

Fig. 13. Comparison of the Taylor–Culick streamlines (—) and the Type II energy-maximized solution (*q* = 2) with stretched streamline curvature (−·−). Results are shown for an inert

*<sup>n</sup>*=0(−1)*n*(2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)−*<sup>q</sup>* <sup>=</sup> <sup>4</sup>*<sup>q</sup>*

Equation (94) yields the general structure of the Type II complementary family of solutions

*<sup>k</sup>*=0(−1)*k*(2*<sup>k</sup>* <sup>+</sup> <sup>1</sup>)−*<sup>q</sup>*

4 )

Note that the Type II solutions emerging from (98) dispose of kinetic energies that are higher than Taylor-Culick's. The variation of the solution with respect to *q* is embodied in the upper

asymptotically by taking the limit as *q* → ∞. Here too, most of the solutions exhibit energies that fall within one percent of Taylor-Culick's. The most interesting solutions are those corresponding to *q* = 2, 3, and 4 with energies that are 47.0, 8.08, and 2.4 percent larger than Taylor-Culick's. When the energy level is fixed at *q* = 2, a simplification follows for the

The Type II solution that carries the most energy at *q* = 2 is plotted in Figure 13 and listed in Table 4. In Figure 13, the Type II approximation is seen to overshoot the Taylor-Culick streamline curvature. In view of the two types of solutions with energies that either lag or surpass that of Taylor-Culick's, one may perceive the *q* → ∞ case as a saddle point to which other possible forms will quickly converge when their energies are shifted. Later in Section 3.12, we will use the entropy maximization principle to establish the Taylor–Culick model as

<sup>4</sup> ) <sup>−</sup> *<sup>ζ</sup>*(*q*, <sup>3</sup>


Table 4. Summary of solutions with most kinetic energy for various headwall injection patterns. Here, *<sup>χ</sup><sup>n</sup>* <sup>≡</sup> <sup>1</sup> <sup>2</sup> (2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)*πr*2.

Fig. 14. Effect of the kinetic energy power index on the two types of energy-based solutions. Results are shown for (a) turn angle, (b) radial velocity, and (c) axial velocity.

*w*(*r*, 0) Ω−(*r*, *z*) Ω+(*r*, *z*)

Internal Flows Driven by Wall-Normal Injection 123

*W*<sup>c</sup> 0 Ω<sup>+</sup>

*<sup>W</sup>*c(<sup>1</sup> <sup>−</sup> *<sup>r</sup>*2) <sup>2</sup>*W*c*<sup>r</sup>* <sup>Ω</sup><sup>+</sup>

<sup>2</sup>*πr*2) <sup>Ω</sup><sup>+</sup>

distribution whereas the highest speed will emerge in the narrowest and most elongated profile connected with the state of most kinetic energy. Interestingly, although this profile slowly diverges at the centerline, it observes mass conservation. This may be explained by the

Having fully determined the velocity field, its vorticity companion may be determined from

This expression is evaluated for the least and most kinetic energy forms (*q* = 2) and provided

For the least kinetic energy solution (Type I, *q* = 2), the linear variation that accompanies the radial velocity as well as the uniformity of the axial velocity are characteristics of an irrotational motion. The vorticity in this case vanishes and the corresponding velocity field collapses into **u** = −*r***e***<sup>r</sup>* + 2*z***e***z*. This potential analogue of the Taylor–Culick velocity has been historically used by McClure et al. (1963) and Hart & McClure (1959) in modeling the internal

> ∞ ∑ *n*=0

By substituting *un* and *wn* from (72) into (4b), the pressure eigenmodes may be integrated.

(2*n* + 1)

*<sup>n</sup>* <sup>−</sup> <sup>1</sup> 2 *α*2 *n*

<sup>2</sup>*α*<sup>2</sup> *<sup>n</sup>* <sup>−</sup> <sup>1</sup> 2 ∞ ∑ *n*=0

<sup>2</sup> (2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)2*π*2*z*2*α*<sup>2</sup>

flow in SRMs. It is recovered here as an extreme state with the lowest kinetic energy.

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

(2*n* + 1)

∞ ∑ *n*=0 ref <sup>≡</sup> *<sup>π</sup>*<sup>2</sup>

<sup>2</sup><sup>C</sup> *rz* csc( <sup>1</sup>

ref

ref + 2*W*c*r*

ref <sup>+</sup> *<sup>π</sup>W*c*<sup>r</sup>* sin( <sup>1</sup>

<sup>2</sup>*πr* 2)

<sup>2</sup>*πr*2)

<sup>2</sup>*αnz* sin *χ<sup>n</sup>* (101)

*pn*(*r*, *z*) (102)

*α*2 *n*

*<sup>r</sup>*<sup>2</sup> sin2 *<sup>χ</sup><sup>n</sup>* (103)

*<sup>r</sup>*<sup>2</sup> sin2 *<sup>χ</sup><sup>n</sup>* (104)

0 0 Ω<sup>+</sup>

<sup>2</sup>*πr*2) *<sup>π</sup>W*c*<sup>r</sup>* sin( <sup>1</sup>

<sup>Ω</sup> = <sup>Ω</sup>*<sup>θ</sup>* = *<sup>π</sup>*2*<sup>r</sup>*

*W*<sup>c</sup> cos( <sup>1</sup>

*rw*+(*r*, *z*) = 0.

fact that lim*r*→<sup>0</sup>

**3.6.4 Vorticity**

in Table 5.

One gets

**3.6.5 Irrotational motion**

**3.7 Pressure evaluation**

One may approximate the pressure by taking

*pn* <sup>=</sup> *<sup>p</sup>*<sup>0</sup> <sup>−</sup> <sup>1</sup>

*<sup>p</sup>*(*r*, *<sup>z</sup>*) = *<sup>p</sup>*<sup>0</sup> <sup>−</sup> <sup>1</sup>

The total pressure is then determined by summing over all eigensolutions

2*π*2*z*<sup>2</sup> <sup>∞</sup> ∑ *n*=0

Table 5. Vorticity for least or most kinetic energy solutions.

#### **3.6 Behavior of the velocity and vorticity fields**

The sensitivity of the solution to the energy power index *q* is illustrated in Figure 14 where both components of the velocity are displayed in addition to the streamline turn angle *θ*.

#### **3.6.1 Turn angle**

This angle represents the slope of the local velocity measured from the radial injection direction. Making use of axial similarity, *θ* may be expressed as

$$\theta(r) = \frac{180}{\pi} \tan^{-1} \left( -\frac{1}{z} \frac{w}{u} \right) \tag{100}$$

The turn angle is shown in Figure 14(a) where, irrespective of *q*, the flow enters radially at the sidewall with *θ*(1) = 0. This feature confirms that the flow enters the chamber perpendicularly to the surface in fulfilment of the no-slip requirement. Conversely, for all cases considered, the establishment of strictly parallel motion along the centerline is reflected in *θ*(0) = 90. Crossing the region between the wall and the centerline, the *q* = 2 Type I case is accompanied by the sharpest change in the turn angle from 0 to 90 degrees. However, as we shift toward the state of most kinetic energy, the smoothing process causes the turn angle to change more gradually. This may be explained by the relative magnitudes of the radial and axial velocities. Specifically, for the Type I, *q* = 2 case, the axial velocity remains practically constant at any chamber cross-section, whereas the radial velocity magnitude increases with *r*. As we cross into the Type II region, the flow starts turning in the vicinity of the sidewall and progresses smoothly as the centerline is approached.

#### **3.6.2 Radial velocity**

The radial velocity is illustrated in Figure 14(b) for representative energy power indices. Starting with the Type II region, the *q* = 2 solution is seen to exhibit a maximum radial velocity overshoot of 16.5 percent relative to the sidewall injection speed. This overshoot reaches its peak at *r* = 0.66 and is required to compensate for the decreasing circumferential area (2*πrL*) normal to the injected stream. Recalling that the Taylor-Culick radial velocity exhibits a 7 percent overshoot at *r* = 0.861, the maximum overshoot calculated here is more than twice as large; it also occurs at a greater distance from the sidewall. Overall, the Type II solutions exhibit smoother curvatures as *q* is increased. In contrast, by examining the case of least kinetic energy in Figure 14(b), no radial overshoot is observed. Instead, the radial velocity displays its lowest absolute value by diminishing linearly from 1 at the wall to 0 at the centerline. This linear variation is accompanied by an essentially uniform axial velocity depicted for the Type I *q* = 2 case in Figure 14(c). At the outset, the locus of the overshoot varies between 0.66 < *r* < 1 as one moves from the Type II, *q* = 2 to the Type I, *q* = 2 case.

#### **3.6.3 Axial velocity**

In Figure 14(c), it is clear that the Type I axial velocities are initially blunt, with the flattest curve being the one corresponding to the top-hat profile at *q* = 2. As *q* is increased, all curves evolve into a sinusoid that approaches the Taylor–Culick model for *q* = 5 and above. Furthermore, as we cross into the Type II region, the centerline velocity continues to increase with increasing energy levels. Due to mass conservation, *Q* = 2*π* 1 0 *wr* d*r* = 2*πz*, and so the centerline speed at each power index is compelled to vary with its corresponding shape to preserve *Q*. The lowest centerline speed will thus accompany the spatially uniform


Table 5. Vorticity for least or most kinetic energy solutions.

distribution whereas the highest speed will emerge in the narrowest and most elongated profile connected with the state of most kinetic energy. Interestingly, although this profile slowly diverges at the centerline, it observes mass conservation. This may be explained by the fact that lim*r*→<sup>0</sup> *rw*+(*r*, *z*) = 0.

#### **3.6.4 Vorticity**

28 Will-be-set-by-IN-TECH

The sensitivity of the solution to the energy power index *q* is illustrated in Figure 14 where both components of the velocity are displayed in addition to the streamline turn angle *θ*.

This angle represents the slope of the local velocity measured from the radial injection

tan−<sup>1</sup> −1 *z w u* 

The turn angle is shown in Figure 14(a) where, irrespective of *q*, the flow enters radially at the sidewall with *θ*(1) = 0. This feature confirms that the flow enters the chamber perpendicularly to the surface in fulfilment of the no-slip requirement. Conversely, for all cases considered, the establishment of strictly parallel motion along the centerline is reflected in *θ*(0) = 90. Crossing the region between the wall and the centerline, the *q* = 2 Type I case is accompanied by the sharpest change in the turn angle from 0 to 90 degrees. However, as we shift toward the state of most kinetic energy, the smoothing process causes the turn angle to change more gradually. This may be explained by the relative magnitudes of the radial and axial velocities. Specifically, for the Type I, *q* = 2 case, the axial velocity remains practically constant at any chamber cross-section, whereas the radial velocity magnitude increases with *r*. As we cross into the Type II region, the flow starts turning in the vicinity of the sidewall

The radial velocity is illustrated in Figure 14(b) for representative energy power indices. Starting with the Type II region, the *q* = 2 solution is seen to exhibit a maximum radial velocity overshoot of 16.5 percent relative to the sidewall injection speed. This overshoot reaches its peak at *r* = 0.66 and is required to compensate for the decreasing circumferential area (2*πrL*) normal to the injected stream. Recalling that the Taylor-Culick radial velocity exhibits a 7 percent overshoot at *r* = 0.861, the maximum overshoot calculated here is more than twice as large; it also occurs at a greater distance from the sidewall. Overall, the Type II solutions exhibit smoother curvatures as *q* is increased. In contrast, by examining the case of least kinetic energy in Figure 14(b), no radial overshoot is observed. Instead, the radial velocity displays its lowest absolute value by diminishing linearly from 1 at the wall to 0 at the centerline. This linear variation is accompanied by an essentially uniform axial velocity depicted for the Type I *q* = 2 case in Figure 14(c). At the outset, the locus of the overshoot varies between 0.66 < *r* < 1 as one moves from the Type II, *q* = 2 to the Type I, *q* = 2 case.

In Figure 14(c), it is clear that the Type I axial velocities are initially blunt, with the flattest curve being the one corresponding to the top-hat profile at *q* = 2. As *q* is increased, all curves evolve into a sinusoid that approaches the Taylor–Culick model for *q* = 5 and above. Furthermore, as we cross into the Type II region, the centerline velocity continues to

and so the centerline speed at each power index is compelled to vary with its corresponding shape to preserve *Q*. The lowest centerline speed will thus accompany the spatially uniform

increase with increasing energy levels. Due to mass conservation, *Q* = 2*π*

(100)

 1 0

*wr* d*r* = 2*πz*,

**3.6 Behavior of the velocity and vorticity fields**

direction. Making use of axial similarity, *θ* may be expressed as

and progresses smoothly as the centerline is approached.

*<sup>θ</sup>*(*r*) = <sup>180</sup> *π*

**3.6.1 Turn angle**

**3.6.2 Radial velocity**

**3.6.3 Axial velocity**

Having fully determined the velocity field, its vorticity companion may be determined from

$$
\Omega = \Omega\_{\theta} = \pi^2 r \sum\_{n=0}^{\infty} (2n+1)^2 \alpha\_n z \sin \chi\_n \tag{101}
$$

This expression is evaluated for the least and most kinetic energy forms (*q* = 2) and provided in Table 5.

#### **3.6.5 Irrotational motion**

For the least kinetic energy solution (Type I, *q* = 2), the linear variation that accompanies the radial velocity as well as the uniformity of the axial velocity are characteristics of an irrotational motion. The vorticity in this case vanishes and the corresponding velocity field collapses into **u** = −*r***e***<sup>r</sup>* + 2*z***e***z*. This potential analogue of the Taylor–Culick velocity has been historically used by McClure et al. (1963) and Hart & McClure (1959) in modeling the internal flow in SRMs. It is recovered here as an extreme state with the lowest kinetic energy.

#### **3.7 Pressure evaluation**

One may approximate the pressure by taking

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

By substituting *un* and *wn* from (72) into (4b), the pressure eigenmodes may be integrated. One gets

$$p\_n = p\_0 - \frac{1}{2}(2n+1)^2 \pi^2 z^2 \alpha\_n^2 - \frac{1}{2} \frac{\alpha\_n^2}{r^2} \sin^2 \chi\_n \tag{103}$$

The total pressure is then determined by summing over all eigensolutions

$$p(r,z) = p\_0 - \frac{1}{2}\pi^2 z^2 \sum\_{n=0}^{\infty} (2n+1)^2 \alpha\_n^2 - \frac{1}{2} \sum\_{n=0}^{\infty} \frac{a\_n^2}{r^2} \sin^2 \chi\_n \tag{104}$$

23456

Fig. 16. Asymptotic behavior of the kinetic energy density for both Type I (- - -) and Type II (+

solutions bridges the gap between an essentially potential flow at *q* = 2 and a fully rotational field at *q* → ∞, thus yielding intermediate formulations with energies that vary across the

Using the absolute convergence and ratio tests, the series representations can be individually shown to be unconditionally convergent for *q* ≥ 2. The most subtle solutions to examine correspond to the Type II inert headwall case with maximum kinetic energy. The attendant velocity and vorticity forms require special attention. For the sake of illustration, we consider

(2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)*<sup>q</sup>* sin *<sup>χ</sup>n*; *<sup>χ</sup><sup>n</sup>* <sup>≡</sup> <sup>1</sup>

 ≤ ∞ ∑ *n*=0

> *e i* 1 <sup>2</sup> *<sup>π</sup>r*<sup>2</sup>

where the right-hand-side converges for *q* > 1. In evaluating quantities that require one or more differentiations (such as the vorticity), we find it useful to substitute, whenever possible, the closed-form analytical representations of the series in question. The equivalent finite expressions enable us to overcome the pitfalls of term-by-term differentiation which, for some infinite series, can lead to spurious results. The Type II axial velocity for the inert headwall configuration presents such an example at *q* = 2. This series can be collapsed into a

While term-by-term differentiation of the infinite series representation of *w*<sup>+</sup> diverges, the

<sup>2</sup><sup>C</sup> *<sup>r</sup>* csc( <sup>1</sup>

1

+ tanh−<sup>1</sup>

 *e* −*i* 1 <sup>2</sup> *<sup>π</sup>r*<sup>2</sup>

<sup>2</sup>*πr*2) (112)

(111)

Taylor-Culick

Internal Flows Driven by Wall-Normal Injection 125

q

3 12 <sup>∞</sup> • <sup>∞</sup> =

 Type I Type II

<sup>2</sup> (2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)*πr*<sup>2</sup> (109)

(2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)*<sup>q</sup>* (110)

2

∞ .

the Type II streamfunction, specifically

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

∞ ∑ *n*=0

The absolute convergence test may be applied to show that ∞ ∑ *n*=0

 

combination of inverse hyperbolic tangent functions by writing

(2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>) cos *<sup>χ</sup><sup>n</sup>* <sup>=</sup> <sup>1</sup>

derivative of the closed-form equivalent yields the correct outcome of

2C tanh−<sup>1</sup>

<sup>Ω</sup><sup>+</sup> <sup>=</sup> <sup>−</sup> *<sup>π</sup>*

*Bq*

1 (2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)*<sup>q</sup>* sin *<sup>χ</sup><sup>n</sup>*

**3.9 Convergence properties**

*w*<sup>+</sup> =

∞ ∑ *n*=0

*B*2

+ +) solutions.

range [0.81 <sup>−</sup> <sup>1</sup>]<sup>E</sup> <sup>∞</sup>

18%

47%

3

4

<sup>∞</sup>• 3.79

2.09

Fig. 15. Centerline pressure drop for the two types of energy-driven solutions.

where *p*<sup>0</sup> = *p*(0, 0). Interestingly, the pressure drop along the centerline collapses into

$$
\Delta p = p(0, z) - p\_0 = -\frac{1}{2}\pi^2 z^2 \sum\_{n=0}^{\infty} (2n+1)^2 a\_n^2 \tag{105}
$$

Equation (105) is plotted in Figure 15 for *q* = 2, 3, and ∞. Unsurprisingly, the largest pressure excursion is seen to accompany the Type II state with most kinetic energy while the smallest pressure loss is accrued in the least kinetic energy expression, specifically, in the *q* = 2 potential case.

#### **3.8 Asymptotic limits of the kinetic energy density**

When the large *L* approximation is employed with *q* = 2, the Type II kinetic energy density E <sup>+</sup> approaches a constant value of E <sup>+</sup> <sup>∞</sup> (2) = *<sup>π</sup>*5/(96<sup>C</sup> <sup>2</sup>) <sup>≈</sup> 3.79944. Note that the asymptotic value for Taylor–Culick's (i.e. when both *L* and *q* approach infinity), E <sup>∞</sup> <sup>∞</sup> <sup>≡</sup> *<sup>π</sup>*3/12 <sup>≈</sup> 2.5838, is recovered as *q* → ∞. In general, when *L* → ∞, the limit of the kinetic energy density can be written as

$$\mathcal{E}\_{\infty} = \frac{1}{12}\pi^3 \sum\_{n=0}^{\infty} (2n+1)^2 a\_n^2 = \mathcal{E}\_{\infty}^{\infty} \sum\_{n=0}^{\infty} (2n+1)^2 a\_n^2 \tag{106}$$

For the Type I solutions, substitution of (106) yields a closed-form expression,

$$\mathcal{E}\_{\infty}^{-}(q) = \mathcal{E}\_{\infty}^{\infty} \left[ \sum\_{k=0}^{\infty} \frac{1}{(2k+1)^q} \right]^{-2} \sum\_{n=0}^{\infty} (2n+1)^{2-2q} = \mathcal{E}\_{\infty}^{\infty} \frac{4^q - 4}{(2^q - 1)^2} \frac{\zeta(2q-2)}{\zeta(q)^2} \tag{107}$$

In like manner, for the Type II solutions, (106) leads to

$$\mathcal{E}\_{\infty}^{+}(q) = \mathcal{E}\_{\infty}^{\infty} \left[ \sum\_{k=0}^{\infty} \frac{(-1)^{k}}{(2k+1)^{q}} \right]^{-2} \sum\_{n=0}^{\infty} (2n+1)^{2-2q} = \mathcal{E}\_{\infty}^{\infty} \frac{4^{q}(4^{q}-4)\zeta(2q-2)}{[\zeta(q, \frac{1}{4}) - \zeta(q, \frac{3}{4})]^{2}} \tag{108}$$

As shown in Figure 16 both types approach E <sup>∞</sup> <sup>∞</sup> either from below or above, depending on *q*. The Taylor–Culick limit of 2.5838 is practically reached by both Type I and Type II solutions with differences of less than 0.287 and 0.265 percent at *q* = 6. The maximum range occurs at *q* = 2 while the total allowable excursion in energy that the mean flow can undergo may be estimated at E <sup>+</sup> <sup>∞</sup> (2) − E <sup>−</sup> <sup>∞</sup> (2) /E <sup>∞</sup> <sup>∞</sup> = 0.66. From an academic standpoint, the Type I family of 30 Will-be-set-by-IN-TECH

0 0.25L 0.5L 0.75L L

Equation (105) is plotted in Figure 15 for *q* = 2, 3, and ∞. Unsurprisingly, the largest pressure excursion is seen to accompany the Type II state with most kinetic energy while the smallest pressure loss is accrued in the least kinetic energy expression, specifically, in the *q* = 2

When the large *L* approximation is employed with *q* = 2, the Type II kinetic energy density

is recovered as *q* → ∞. In general, when *L* → ∞, the limit of the kinetic energy density can be

(2*n* + 1)

The Taylor–Culick limit of 2.5838 is practically reached by both Type I and Type II solutions with differences of less than 0.287 and 0.265 percent at *q* = 6. The maximum range occurs at *q* = 2 while the total allowable excursion in energy that the mean flow can undergo may be

*<sup>n</sup>* = <sup>E</sup> <sup>∞</sup> ∞ ∞ ∑ *n*=0

(2*n* + 1)2−2*<sup>q</sup>* = E <sup>∞</sup>

<sup>2</sup>−2*<sup>q</sup>* = E <sup>∞</sup> ∞

(2*n* + 1)2*α*<sup>2</sup>

0.3L 0.4L

where *p*<sup>0</sup> = *p*(0, 0). Interestingly, the pressure drop along the centerline collapses into

Fig. 15. Centerline pressure drop for the two types of energy-driven solutions.

<sup>Δ</sup>*<sup>p</sup>* <sup>=</sup> *<sup>p</sup>*(0, *<sup>z</sup>*) <sup>−</sup> *<sup>p</sup>*<sup>0</sup> <sup>=</sup> <sup>−</sup><sup>1</sup>

value for Taylor–Culick's (i.e. when both *L* and *q* approach infinity), E <sup>∞</sup>

For the Type I solutions, substitution of (106) yields a closed-form expression,

−<sup>2</sup> <sup>∞</sup> ∑ *n*=0

<sup>−</sup><sup>2</sup> <sup>∞</sup> ∑ *n*=0

12*π*<sup>3</sup> <sup>∞</sup> ∑ *n*=0

1 (2*k* + 1)*<sup>q</sup>*

(−1)*<sup>k</sup>* (2*k* + 1)*<sup>q</sup>* = 2

(2*n* + 1)2*α*<sup>2</sup>

<sup>∞</sup> (2) = *<sup>π</sup>*5/(96<sup>C</sup> <sup>2</sup>) <sup>≈</sup> 3.79944. Note that the asymptotic

(2*n* + 1)2*α*<sup>2</sup>

<sup>4</sup>*<sup>q</sup>* <sup>−</sup> <sup>4</sup> (2*<sup>q</sup>* <sup>−</sup> <sup>1</sup>)<sup>2</sup>

[*ζ*(*q*, <sup>1</sup>

<sup>∞</sup> = 0.66. From an academic standpoint, the Type I family of

∞

q

Taylor-Culick

z

2*π*2*z*<sup>2</sup> <sup>∞</sup> ∑ *n*=0

least kinetic energy state Type I

potential

3 3

2

*<sup>n</sup>* (105)

<sup>∞</sup> <sup>≡</sup> *<sup>π</sup>*3/12 <sup>≈</sup> 2.5838,

*<sup>n</sup>* (106)

*<sup>ζ</sup>*(*q*)<sup>2</sup> (107)

<sup>4</sup> )]<sup>2</sup> (108)

*ζ*(2*q* − 2)

<sup>4</sup>*q*(4*<sup>q</sup>* <sup>−</sup> <sup>4</sup>)*ζ*(2*<sup>q</sup>* <sup>−</sup> <sup>2</sup>)

<sup>∞</sup> either from below or above, depending on *q*.

<sup>4</sup> ) <sup>−</sup> *<sup>ζ</sup>*(*q*, <sup>3</sup>


potential case.

written as

E −

E <sup>+</sup>

estimated at

<sup>∞</sup> (*q*) = <sup>E</sup> <sup>∞</sup> ∞

<sup>∞</sup> (*q*) = <sup>E</sup> <sup>∞</sup> ∞

E <sup>+</sup>

<sup>∞</sup> (2) − E <sup>−</sup>




Δp

0

most kinetic energy state Type II

> -1.0 -0.8 -0.6 -0.4

**3.8 Asymptotic limits of the kinetic energy density**

E<sup>∞</sup> = <sup>1</sup>

In like manner, for the Type II solutions, (106) leads to

 ∞ ∑ *k*=0

 ∞ ∑ *k*=0

As shown in Figure 16 both types approach E <sup>∞</sup>

<sup>∞</sup> (2) /E <sup>∞</sup>

E <sup>+</sup> approaches a constant value of E <sup>+</sup>

Fig. 16. Asymptotic behavior of the kinetic energy density for both Type I (- - -) and Type II (+ + +) solutions.

solutions bridges the gap between an essentially potential flow at *q* = 2 and a fully rotational field at *q* → ∞, thus yielding intermediate formulations with energies that vary across the range [0.81 <sup>−</sup> <sup>1</sup>]<sup>E</sup> <sup>∞</sup> ∞ .

#### **3.9 Convergence properties**

Using the absolute convergence and ratio tests, the series representations can be individually shown to be unconditionally convergent for *q* ≥ 2. The most subtle solutions to examine correspond to the Type II inert headwall case with maximum kinetic energy. The attendant velocity and vorticity forms require special attention. For the sake of illustration, we consider the Type II streamfunction, specifically

$$\psi(r,z) = z \sum\_{n=0}^{\infty} \frac{B\_q}{(2n+1)^q} \sin \chi\_n; \quad \chi\_n \equiv \frac{1}{2}(2n+1)\pi r^2 \tag{109}$$

The absolute convergence test may be applied to show that

$$\sum\_{n=0}^{\infty} \left| \frac{1}{(2n+1)^q} \sin \chi\_n \right| \le \sum\_{n=0}^{\infty} \frac{1}{(2n+1)^q} \tag{110}$$

where the right-hand-side converges for *q* > 1. In evaluating quantities that require one or more differentiations (such as the vorticity), we find it useful to substitute, whenever possible, the closed-form analytical representations of the series in question. The equivalent finite expressions enable us to overcome the pitfalls of term-by-term differentiation which, for some infinite series, can lead to spurious results. The Type II axial velocity for the inert headwall configuration presents such an example at *q* = 2. This series can be collapsed into a combination of inverse hyperbolic tangent functions by writing

$$\log^{+} = \sum\_{n=0}^{\infty} \frac{B\_2}{(2n+1)} \cos \chi\_{\mathbb{R}^3} = \frac{1}{2\ell'} \left[ \tanh^{-1} \left( e^{i\frac{1}{2}\pi r^2} \right) + \tanh^{-1} \left( e^{-i\frac{1}{2}\pi r^2} \right) \right] \tag{111}$$

While term-by-term differentiation of the infinite series representation of *w*<sup>+</sup> diverges, the derivative of the closed-form equivalent yields the correct outcome of

$$
\Omega^+ = -\frac{\pi}{2\ell} r \csc(\frac{1}{2}\pi r^2) \tag{112}
$$

0 0.5 1

Type II <sup>2</sup> <sup>3</sup> Taylor-Culick

0 0.5 1

(b) *F*� (*r*)

Fig. 17. Comparison between analytical (—) and numerical (◦) solutions for (a) *F*(*r*), and (b)

(*r*) for Type I (blue) and Type II (red) solutions. Plots are shown for *q* = 2, 3 and ∞. Here,

To explore the physicality of our variational solutions, a second law analysis is helpful. To better understand the mechanisms responsible for the system to opt for one energy state over another, or one type of solution over another, the principle of entropy maximization may be referred to. This principle states that a system will tend to maximize entropy at equilibrium and may hence be applied to our problem by considering the different energy solutions as different states of the same system. As shown by Saad & Majdalani (2010), the second law analysis reveals that the volumetric entropy of the Type I family grows with successive increases in *q* but depreciates in the Type II case. So given an initial profile, the system may

If the system is initialized on the Type I branch, it will evolve toward the Taylor–Culick solution to the extent of maximizing its total entropy. While entropy could be further increased

Hart-McClure potential

Type I

(a) *F*(*r*)

Taylor-Culick

Hart-McClure potential

r

r

Type I

3

3

q = 2

2

3

Type II

Internal Flows Driven by Wall-Normal Injection 127

0

2

0

**3.12 Unphysicality of the Type II family of solutions**

evolve according to one of two scenarios that are described below.

*F*�

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

**3.12.1 Type I branching**

0.5

1.0

1.5

F*'*

q = 2

0.25

0.50

0.75

F

1

As it may be expected, the corresponding solution is accompanied by finite kinetic energy and mass flowrate despite its singularity at the centerline.

#### **3.10 Arbitrary headwall injection**

For T-burners, solid rocket motors with reactive fore-ends, and hybrid rocket chambers with injector faceplates, a model that accounts for headwall injection is required. For these problems, our analysis may be repeated assuming an injecting headwall with an axisymmetrically varying profile defined by (8). The streamfunction becomes

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

In the resulting expressions, *βn* does not vanish. As shown by Majdalani & Saad (2007b) and detailed in Section 2.1, orthogonality may be applied to obtain *βn* for an axisymmetric headwall injection pattern. Application of Lagrangian optimization in conjunction with the large *L* approximation yield identical results for *αn* as those obtained in (92) and (98). The streamfunction, axial velocity, and vorticity for several injection profiles are available through Tables 3, 4, and 5 where the least and most kinetic energy solutions are identified.

#### **3.11 Numerical verification**

Our analytical expansions may be verified by solving (13) using Runge-Kutta integration. We begin by introducing the transformation *ψ* = *z f*(*r*) through which (13) may be reduced to a second order ODE

$$F''(r) - \frac{1}{r}F'(r) + \mathcal{C}^2 r^2 F(r) = 0\tag{114}$$

In order to numerically capture the different variational solutions, the boundary conditions of (14) have to be carefully selected. Because our solutions are in series form, we first decompose *F*(*r*) into its eigenmode components by taking

$$F(r) = \sum\_{n=0}^{\infty} F\_n(r) \tag{115}$$

and so (114) becomes

$$F\_{n}^{\prime\prime}(r) - \frac{1}{r}F\_{n}^{\prime}(r) + \mathcal{C}\_{n}^{2}r^{2}F\_{n}(r) = 0; \quad n = 0, 1, \cdots, \infty \tag{116}$$

where *n* corresponds to the eigenmode associated with *Cn* = (2*n* + 1)*π*. Finally, the boundary conditions may be written as

$$\begin{cases} F\_{\mathfrak{n}}(0) = 0; & F(0) = \sum\_{n=0}^{\infty} F\_{\mathfrak{n}}(0) = 0 \\ F\_{\mathfrak{n}}(1) = (-1)^{n} a\_{\mathfrak{n}}; & F(1) = \sum\_{n=0}^{\infty} F\_{\mathfrak{n}}(1) = \sum\_{n=0}^{\infty} (-1)^{n} a\_{\mathfrak{n}} = 1 \end{cases} \tag{117}$$

Using 120 terms to reconstruct the series expansions, both numerical and analytical solutions for *F*(*r*) and *F*� (*r*) are displayed in Figures 17(a) and 17(b), respectively. This comparison is held at representative values of the kinetic energy power index corresponding to *q* = 2, 3, and ∞. It is gratifying that, irrespective of *q* and *n*, the variational solutions are faithfully simulated by the numerical data to the extent that visual differences between full circles (numerical) and solid lines (analytical) are masked.

32 Will-be-set-by-IN-TECH

As it may be expected, the corresponding solution is accompanied by finite kinetic energy and

For T-burners, solid rocket motors with reactive fore-ends, and hybrid rocket chambers with injector faceplates, a model that accounts for headwall injection is required. For these problems, our analysis may be repeated assuming an injecting headwall with an

(*αnz* + *βn*) sin[ <sup>1</sup>

In the resulting expressions, *βn* does not vanish. As shown by Majdalani & Saad (2007b) and detailed in Section 2.1, orthogonality may be applied to obtain *βn* for an axisymmetric headwall injection pattern. Application of Lagrangian optimization in conjunction with the large *L* approximation yield identical results for *αn* as those obtained in (92) and (98). The streamfunction, axial velocity, and vorticity for several injection profiles are available through

Our analytical expansions may be verified by solving (13) using Runge-Kutta integration. We begin by introducing the transformation *ψ* = *z f*(*r*) through which (13) may be reduced to a

(*r*) + *C*2*r*

In order to numerically capture the different variational solutions, the boundary conditions of (14) have to be carefully selected. Because our solutions are in series form, we first decompose

> ∞ ∑ *n*=0

where *n* corresponds to the eigenmode associated with *Cn* = (2*n* + 1)*π*. Finally, the boundary

∑ *n*=0

∑ *n*=0

Using 120 terms to reconstruct the series expansions, both numerical and analytical solutions

held at representative values of the kinetic energy power index corresponding to *q* = 2, 3, and ∞. It is gratifying that, irrespective of *q* and *n*, the variational solutions are faithfully simulated by the numerical data to the extent that visual differences between full circles (numerical) and

*Fn*(0) = 0

*Fn*(1) = <sup>∞</sup> ∑ *n*=0

(*r*) are displayed in Figures 17(a) and 17(b), respectively. This comparison is

*F*(*r*) =

*<sup>n</sup>*(*r*) + *<sup>C</sup>*<sup>2</sup>

*Fn*(0) = 0; *<sup>F</sup>*(0) = <sup>∞</sup>

*Fn*(1)=(−1)*nαn*; *<sup>F</sup>*(1) = <sup>∞</sup>

<sup>2</sup> (2*n* + 1)*πr*

<sup>2</sup>] (113)

<sup>2</sup>*F*(*r*) = 0 (114)

*Fn*(*r*) (115)

*nr*<sup>2</sup>*Fn*(*r*) = 0; *<sup>n</sup>* <sup>=</sup> 0, 1, ··· , <sup>∞</sup> (116)

(−1)*nα<sup>n</sup>* <sup>=</sup> <sup>1</sup>

(117)

axisymmetrically varying profile defined by (8). The streamfunction becomes

Tables 3, 4, and 5 where the least and most kinetic energy solutions are identified.

*<sup>F</sup>*��(*r*) <sup>−</sup> <sup>1</sup> *r F*�

∞ ∑ *n*=0

mass flowrate despite its singularity at the centerline.

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

**3.10 Arbitrary headwall injection**

**3.11 Numerical verification**

*F*(*r*) into its eigenmode components by taking

*F*�� *<sup>n</sup>* (*r*) <sup>−</sup> <sup>1</sup> *r F*�

second order ODE

and so (114) becomes

for *F*(*r*) and *F*�

conditions may be written as

⎧ ⎪⎪⎨

⎪⎪⎩

solid lines (analytical) are masked.

Fig. 17. Comparison between analytical (—) and numerical (◦) solutions for (a) *F*(*r*), and (b) *F*� (*r*) for Type I (blue) and Type II (red) solutions. Plots are shown for *q* = 2, 3 and ∞. Here, *ψ*(*r*, *z*) = *zF*(*r*).

#### **3.12 Unphysicality of the Type II family of solutions**

To explore the physicality of our variational solutions, a second law analysis is helpful. To better understand the mechanisms responsible for the system to opt for one energy state over another, or one type of solution over another, the principle of entropy maximization may be referred to. This principle states that a system will tend to maximize entropy at equilibrium and may hence be applied to our problem by considering the different energy solutions as different states of the same system. As shown by Saad & Majdalani (2010), the second law analysis reveals that the volumetric entropy of the Type I family grows with successive increases in *q* but depreciates in the Type II case. So given an initial profile, the system may evolve according to one of two scenarios that are described below.

#### **3.12.1 Type I branching**

If the system is initialized on the Type I branch, it will evolve toward the Taylor–Culick solution to the extent of maximizing its total entropy. While entropy could be further increased

effort to mimic the internal flow character in either SRM or hybrid rocket motors. Overall, we find that the flow field evolves to the self-similar Taylor–Culick sinusoid far downstream irrespective of the headwall injection pattern. Nonetheless, the details of headwall injection remain important in hybrid motors, short chambers, and T-burners where the foregoing approximations may be applied. In hybrid rockets, our models seem to capture the streamtube

Internal Flows Driven by Wall-Normal Injection 129

The other chief contribution of this chapter is the discussion of variational solutions that may be connected with the Taylor-Culick problem. Based on the Lagrangian optimization of the total volumetric energy in the chamber, we are able to identify two families of solutions with dissimilar energy signatures. These are accompanied by lower or higher kinetic energies that vary, from one end of the spectrum to the other, by up to 66 percent of their mean value.

up to Taylor-Culick's. The latter is asymptotically recovered in the limit of *q* → ∞, a case that corresponds to an equilibrium state with maximum entropy. In practice, most solutions become indiscernible from Taylor-Culick's for *q* ≥ 5. Indeed, those obtained with *q* = 2, 3, and 4 exhibit energies that are 18.9, 8.28, and 2.73 percent lower than their remaining counterparts. The least kinetic energy solution with *q* = 2 returns the classic, irrotational Hart-McClure profile. It can thus be seen that the application of the Lagrangian optimization principle to this problem leads to the potential form that historically preceded the Taylor-Culick motion. It can also be inferred that the Type I solutions not only bridge the gap between a plain potential representation of this problem and a rotational formulation, but also recover a continuous spectrum of approximations that stand in between. When the same analysis is repeated using

*<sup>n</sup>* <sup>∼</sup> (2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)−*q*; *<sup>q</sup>* <sup>≥</sup> 2, a complementary family of Type II solutions is identified with descending energy levels. These are shown to be purely academic, although they represent a class of exact solutions to the modified Helmholtz equation. Their most notable profiles correspond to *q* = 2, 3, and 4 with energies that are 47.0, 8.08, and 2.40 percent higher than Taylor-Culick's. Their entropies are also higher than that associated with the equilibrium state. Despite their dissimilar forms, both Type I and II solutions converge to the Taylor-Culick representation when their energies are incremented or reduced. Yet before using the new variational solutions to approximate the mean flow profile in porous tubes or the bulk gaseous motion in simulated rocket motors, it should be borne in mind that no direct connection exists between the energy steepened states and turbulence. For this reason, it is hoped that additional numerical and experimental investigations are pursued to test their physicality and the particular configurations in which they are prone to appear. As for the uniqueness of the Taylor-Culick equilibrium state, it may be confirmed from the entropy maximization principle and the Lagrangian-based solutions where, for a given set of boundary conditions, the equilibrium state may be asymptotically restored as *q* → ∞ irrespective of the form of *<sup>α</sup><sup>n</sup>* <sup>∼</sup> (−1)*n*(*p n* <sup>+</sup> *<sup>m</sup>*)−*q*, provided that the Lagrangian constraint <sup>∑</sup> (−1)*nα<sup>n</sup>* <sup>=</sup> 1 is faithfully

Lastly, we note that the collection of variational solutions that admit variable headwall injection increase our repertoire of Euler-based approximations that may be used to model the incompressible motion in porous tubes. For the porous channel flow analogue, the planar solutions are presented by Saad & Majdalani (2008a; 2009b). As for tapered grain

configuration, the reader may consult with Saad et al. (2006) or Sams et al. (2007).

sequence of Type I solutions is unraveled in ascending order, *α*−

*<sup>n</sup>* <sup>∼</sup> (−1)*n*(2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)−<sup>2</sup> yields the profile with least kinetic energy, a

*<sup>n</sup>* <sup>∼</sup> (−1)*n*(2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)−*q*; *<sup>q</sup>* <sup>&</sup>gt; 2,

motion quite effectively.

After identifying that *α*−

*α*+

secured.

by branching out to the Type II region, it may be argued that such development is not possible for two reasons. Firstly, the character of the two types of solutions is sharply dissimilar, especially in the expressions for *α*− *<sup>n</sup>* and *α*<sup>+</sup> *<sup>n</sup>* . Secondly, given that the Taylor–Culick solution maximizes the entropy for the Type I branch, it can be viewed as a local equilibrium state. As such, there is no necessity for the system to switch branches once it reaches the Taylor–Culick state.

#### **3.12.2 Type II branching**

If the system is initialized on the Type II branch, it will approach the solution with most vorticity (i.e. Type II, *q* = 2). Although this may be a mathematically viable outcome, it may not be physically realizable because it would be practically impossible to initialize a system with such a high level of vorticity without the aid of external work. The most natural flow evolution corresponds to an irrotational system originally at rest in which vorticity generation is initiated at the sidewall during the injection process. The ensuing motion will subsequently progress until it reaches the stable Taylor–Culick equilibrium state wherein it can settle with no further tendency to branch out.
