**2. Boundary layer separation models and mean separation location**

Two physically-based models of shock-induced separation were developed in 2007 (Keanini & Brown, 2007). The models assume importance for two reasons. First, they provide verifiable insight into the physical processes underlying shock-separation of compressible turbulent boundary layers. Second, they allow prediction of the altitude-dependent mean separation line, crucial in determining both the location and magnitude of altitude-dependent random side loads (Keanini et al., 2011; Srivastava et al., 2010). This section describes the more refined of these two models, as well as Chapman's Free Interaction Model (Chapman et al., 1958). All three models rely on scaling analyses of the boundary layer and near-boundary flows, in and near the boundary layer separation zone.

#### **2.1 Time-average separation**

Two distinct separation processes have been identified in overexpanded rocket nozzles, *free shock separation,* in which the turbulent boundary layer separates without reattachment, and *restricted shock separation,* in which the separated boundary layer reattaches, forming a small, closed recirculation zone immediately downstream of the separation point; see, e.g., (Keanini & Brown, 2007; Ostlund, 2002) for recent reviews of this work. This Chapter focuses on free shock separation.

The time-average flow features associated with free shock separation in nozzles were first characterized by Summerfield et al. (Summerfield et al., 1954), and are depicted schematically in figures 1 and 2. As shown, the time average pressure along the nozzle wall increases from *Pi* at the incipient separation point, *xi*, to a peak value of *Pp* at *xp*. Depending on the nozzle and the shock location relative to the nozzle exit, *Pp* is typically on the order of 80 to 100 % of the ambient pressure, *Pa*. The time-average separation point, *xs*, lies immediately upstream of *xp*.

Of central importance in nozzle design is determining both the conditions under which separation will occur and the approximate separation location. A number of criteria have been proposed for predicting the nominal free shock separation point, *xs*; see, e.g., (Keanini & Brown, 2007; Ostlund, 2002) for reviews. Since the boundary layer pressure rise between *xi* and *xs* depends primarily on the inviscid flow Mach number, *Mi*, most criteria relate either a gross separation pressure ratio, *Pi*/*Pa*, or more recently, a refined ratio, *Pi*/*Pp*, to *Mi* (Ostlund, 2002). Given the separation pressure ratio, the separation location can then be determined using an appropriate model of flow upstream of separation.

Although the actual separation process is highly dynamic, the scaling model focuses on time average flow dynamics in the vicinity of the shock interaction zone. In order to provide physical context, we briefly review the dynamical features associated with free shock separation and note simplifying assumptions made. Shock motion over the shock interaction zone appears to be comprised of essentially two components: i) a low frequency, large scale motion produced by flow variations downstream of the separation point, and occurring over the length of the shock interaction zone, *lp* = *xp* − *xi*, at characteristic frequencies, *fs* [on the order of 300 to 2000 Hz in the case of compression ramp and backward facing step flows (Dolling & Brusniak, 1989)], and ii) a high frequency, low amplitude jitter produced by advection of vortical structures through the shock interaction zone (Dolling & Brusniak, 1989). The scaling models in (Keanini & Brown, 2007) limit attention to time scales that are long relative to *f* <sup>−</sup><sup>1</sup> *<sup>s</sup>* . In addition, the model assumes that the flow is statistically stationary and that the separation process is two-dimensional.

The time average pressure gradient over the shock interaction zone (*xi* ≤ *x* ≤ *xp*), given approximately by

$$\frac{\partial P}{\partial \mathbf{x}} \sim \frac{P\_p - P\_i}{l\_p} \tag{1}$$

in reality reflects the intermittent, random motion of the shock between *xi* and *xp*; see, e.g., (Dolling & Brusniak, 1989). As the shock-compression wave system oscillates randomly above (and partially within) the boundary layer, the associated pressure jump across the system is transmitted across the boundary layer on a time scale *τ<sup>s</sup>* ∼ *δi*/ <sup>√</sup>*kRTi*, where *<sup>δ</sup><sup>i</sup>* and *Ti* are the characteristic boundary layer thickness and temperature in the vicinity of *xi*. Under typical experimental conditions, *<sup>τ</sup><sup>s</sup>* is much shorter than the slow time scale, *<sup>f</sup>* <sup>−</sup><sup>1</sup> *<sup>s</sup>* (where *<sup>τ</sup><sup>s</sup>* <sup>≈</sup> 1 to 10*μs*); thus, the instantaneous separation point essentially tracks the random position of the shock-compression wave system, where the position of the separation point is described by a Gaussian distribution over the length of the interaction zone (Dolling & Brusniak, 1989).

#### **2.2 Scale analysis of shock-induced separation**

4 Will-be-set-by-IN-TECH

II) The stochastic boundary layer separation models developed in (Keanini et al., 2011; Srivastava et al., 2010) allow construction of stochastic side load models that, on one hand, are physically self-consistent, and on the other, are imbued with statistical properties that are fully consistent with available experimental observations. Our second objective focuses

III) Given physically consistent separation and side load models (Keanini et al., 2011; Srivastava et al., 2010), the effect of random nozzle side loads on rocket ascent can be computed. Our group recently developed (Keanini et al., 2011) a series of interconnected, analytical models describing: i) fast time-scale, altitude-dependent stochastic boundary layer separation, ii) associated short-time- and long(rocket-dynamics)-time-scale stochastic side load generation, and iii) stochastic, altitude-dependent rocket response. In addition, a high-fidelity numerical model which solved the full nonlinear, coupled equations of rocket rotational and translational motion, under the action of altitude-dependent random side loads, was also reported (Srivastava et al., 2010). Our third objective centers on outlining

IV) Our most recent work (Keanini et al., 2011) demonstrates that rocket pitch and yaw rates evolve as Ornstein-Uhlenbeck processes. Since stochastic pitch and yaw rate evolution determines not only the random evolution of pitch/yaw displacement, but also the stochastic evolution of the rocket's lateral velocity and displacement (Keanini et al., 2011), the rocket's rotational and translational dynamics, fundamentally, trace to a damped diffusion process. The last objective centers on highlighting the connection between

Two physically-based models of shock-induced separation were developed in 2007 (Keanini & Brown, 2007). The models assume importance for two reasons. First, they provide verifiable insight into the physical processes underlying shock-separation of compressible turbulent boundary layers. Second, they allow prediction of the altitude-dependent mean separation line, crucial in determining both the location and magnitude of altitude-dependent random side loads (Keanini et al., 2011; Srivastava et al., 2010). This section describes the more refined of these two models, as well as Chapman's Free Interaction Model (Chapman et al., 1958). All three models rely on scaling analyses of the boundary layer and near-boundary flows, in and

Two distinct separation processes have been identified in overexpanded rocket nozzles, *free shock separation,* in which the turbulent boundary layer separates without reattachment, and *restricted shock separation,* in which the separated boundary layer reattaches, forming a small, closed recirculation zone immediately downstream of the separation point; see, e.g., (Keanini & Brown, 2007; Ostlund, 2002) for recent reviews of this work. This Chapter focuses on free

The time-average flow features associated with free shock separation in nozzles were first characterized by Summerfield et al. (Summerfield et al., 1954), and are depicted schematically in figures 1 and 2. As shown, the time average pressure along the nozzle wall increases from *Pi* at the incipient separation point, *xi*, to a peak value of *Pp* at *xp*. Depending on the nozzle and the shock location relative to the nozzle exit, *Pp* is typically on the order of 80 to 100 % of

these analytical and numerical models, and on describing recent results.

stochastic side load-driven rocket response and Ornstein-Uhlenbeck diffusion.

**2. Boundary layer separation models and mean separation location**

near the boundary layer separation zone.

**2.1 Time-average separation**

shock separation.

on describing these new stochastic side load models.

In the vicinity of the separation point,*xs*, we recognize that a fluid particle's normal acceleration component within the separating boundary layer is determined by the normal component of the pressure gradient across the separating boundary layer. Thus, balancing these terms (in the Navier- Stokes equations) yields:

$$
\rho \frac{V\_s^2}{R} \sim \frac{\partial P}{\partial n} \tag{2}
$$

where *Vs* is the particle speed in the streamwise (*s*-)direction, and *R*−<sup>1</sup> is the local streamline curvature. The curvature can be evaluated by first defining the shape of the boundary layer's

Rocket Ascent 7

<sup>161</sup> Shock-Induced Turbulent Boundary Layer Separation in Over-Expanded Rocket Nozzles: Physics, Models, Random Side Loads, and the Diffusive Character of Stochastic Rocket Ascent

flow is determined by the streamwise rate of boundary layer growth. Using (8) in (5) and

*<sup>δ</sup><sup>s</sup>* <sup>∼</sup> *<sup>μ</sup>M*2*a*<sup>2</sup> cos *<sup>θ</sup>*

<sup>∼</sup> *<sup>τ</sup>*<sup>2</sup>

Recognizing that the resultant stress on a fluid particle near *xs* is approximately equal to the vector sum of the horizontally acting viscous shear stress, *τ*2, and the vertically-acting net pressure, Δ*P* = *Pp* − *P*1, then (10) shows that the resultant acts in the direction of the separating boundary layer, *θ*, as it must. Since *δs*/*ls* ∼ tan *θ*, then (10) is also consistent with Chapman's (Chapman et al., 1958) estimate for *τ*2/Δ*P*. Likewise, Chapman et al. (Chapman

*τ<sup>w</sup>* ∼ Δ*P*, then for nominally fixed *θ* (see Sec. II.D below), (10) is also consistent with

Returning to Eqn. (3), inserting (9) and the expression for *R*−<sup>1</sup> and rearranging leads to an

*Pp P*1

point,and where the approximation tan *θ* ≈ −*f* � has been used. Under typical conditions, e.g., those extant in experiments described below [*Mi* <sup>≈</sup> 5, *<sup>θ</sup>* <sup>≈</sup> <sup>16</sup>*o*, *To* <sup>=</sup> 310K, *Po* <sup>=</sup> 1.24 Mpa], *�* <sup>=</sup> *<sup>O</sup>*(10−3), i.e., *�* << 1 [where *<sup>f</sup>* ��|*max* <sup>≈</sup> (*dθ*/*dx*)|*max* <sup>=</sup> *<sup>O</sup>*(1)]. Thus, solving (11) for *Pp*/*P*<sup>1</sup>

∼ *G* − *�*

Importantly, this equation shows that *P*1/*Pp* ≈ *Pi*/*Pp* ∼ *P*1/*P*2, *demonstrating that the separation*

In order to close the approximate model embodied in (12), it is necessary to specify the shock angle, *β*, and the flow deflection angle, *θ*. Referring to earlier work, Summerfield et al. (Summerfield et al., 1954) used measured separation pressure ratios and Mach numbers in the oblique shock relations to infer *θ*; based on their data, they inferred a nominally fixed value, *<sup>θ</sup>* <sup>≈</sup> <sup>16</sup>*o*. By contrast, Ostlund (Ostlund, 2002), again using the same approach,argued that *θ* varies with *Mi*, albeit weakly; he fit his estimate with a linear relationship, *θ* = 1.678*Mi* + 9.347, valid for 2.5 ≤ *Mi* ≤ 4.5. For this range of *Mi*, however, the correlation indicates that *θ* only varies from 13.5*<sup>o</sup>* and 16.9*o*. Based on these indirect estimates, we assume that *θ* is constant; for simplicity, we will arbitrarily adopt the average value of *θ* indicated by Ostlund's correlation, *<sup>θ</sup>* <sup>≈</sup> 15.2*o*, nearly equal to Summerfield's (Summerfield et al., 1954)

2*G*

(*Pp* <sup>−</sup> *<sup>P</sup>*1)*tan<sup>θ</sup>* (9)

<sup>Δ</sup>*<sup>P</sup>* <sup>∼</sup> tan *<sup>θ</sup>* (10)

<sup>1</sup>/2), where *Cf* is the friction factor; since

+ *G*(1 + *�*) ∼ 0 (11)

), *μ*<sup>2</sup> is the gas viscosity near the separation

*<sup>G</sup>* <sup>−</sup> <sup>1</sup> (12)

solving for *δ<sup>s</sup>* finally yields an estimate for the boundary layer thickness near *xs* :

*μU*2/*δ<sup>s</sup> Pp* − *P*<sup>1</sup>

<sup>2</sup> <sup>−</sup> (<sup>1</sup> <sup>+</sup> *<sup>G</sup>*)

<sup>2</sup> *a*<sup>2</sup> cos *θμ*<sup>2</sup> *f* ��/(*P*<sup>1</sup> *f* �

*Pp P*1

where *U*<sup>2</sup> = *M*2*a*<sup>2</sup> cos *θ* and *a*<sup>2</sup> is the sound speed.

et al., 1958) argued that *<sup>δ</sup>*∗/*ls* <sup>∼</sup> *Cf* <sup>=</sup> *<sup>τ</sup>w*/(*ρ*1*u*<sup>2</sup>

 *Pp P*1

and neglecting terms smaller than *O*(*�*), we finally obtain

*pressure ratio essentially corresponds to the oblique shock pressure ratio.*

Chapman's estimate for *Cf* .

where *G* = *P*2/*P*1, *�* = *kM*<sup>3</sup>

**2.4 Shock and flow deflection angles**

expression of the following form:

Before proceeding, and as an aside, we rewrite the estimate in (9) as

outer-most streamline (i.e., a streamline in the vicinity of *δ* which is roughly parallel to the local displacement thickness) as *r*(*x*) = *f*(*x*), where *r*(*x*) is the radial distance from the nozzle centerline to the streamline, evaluated at axial position *x*. Thus, *R*−<sup>1</sup> = *f* ��( 1 + (*f* �)2)−3. Expressing *V*<sup>2</sup> *<sup>s</sup>* in terms of local cartesian velocity components,*V*<sup>2</sup> *<sup>s</sup>* = *u*<sup>2</sup> *<sup>s</sup>* + *v*<sup>2</sup> *<sup>s</sup>*, and estimating *us* and *vs* by their approximate free stream magnitudes downstream of the oblique shock (since again, at axial position *xs*, the boundary layer has passed through the compression system below the shock), we obtain *V*<sup>2</sup> *<sup>s</sup>* <sup>≈</sup> *<sup>U</sup>*<sup>2</sup> <sup>2</sup> <sup>+</sup> *<sup>V</sup>*<sup>2</sup> <sup>2</sup> . Replacing terms in (2) by their approximate magnitudes then leads to

$$
\rho\_2 \frac{\mathcal{U}\_2^2 + V\_2^2}{R} \sim \frac{P\_p - P\_2}{\delta\_\text{s}} \tag{3}
$$

where arguments given in (Keanini & Brown, 2007) lead to the necessary pressure gradient and density estimates.

#### **2.3 Boundary layer thickness**

In order to proceed, we must estimate the magnitude of the boundary layer thickness, *δs*, immediately upstream of the separation point, *xs*. First, note that at the wall between *xi* and *xs*, the x-momentum equation yields the approximate balance

$$
\mu \frac{\partial^2 \mu}{\partial y^2} \sim \frac{\partial P}{\partial x} \tag{4}
$$

where *μ* is the dynamic viscosity. Estimating the magnitude of each term in this equation leads to

$$\frac{\mathcal{U}\_2}{\delta\_s^2} \sim \frac{1}{\mu} \frac{(P\_p - P\_1)}{l\_s} \tag{5}$$

where, since we are focusing on the neighborhood of *xs*, *u* is approximated as *U*2, and where it is recognized that the streamwise pressure increases from approximately *P*<sup>1</sup> near *xi* to approximately *Pp* near *xs*. [Although *P* equals *Ps* at *xs*, due to the relatively small difference between *Ps* and *Pp*, for simplicity, we approximate *P*(*xs*) as *Pp*.] The x-length scale, *ls* = (*xs* − *xi*), is approximately equal to the length of the shock interaction zone, *lp* = *xp* − *xi*. Considering the continuity equation near *xs*, we recognize that since the boundary layer acquires a vertical velocity component as it travels toward and past *xs*, and since associated mass advection and volumetric dilatation terms, *<sup>ρ</sup>*−1**<sup>u</sup>** · ∇*ρ*, and ∇ · **<sup>u</sup>**, respectively, are of the same order, then

$$\frac{\partial u}{\partial x} \approx \frac{\partial v}{\partial y} \tag{6}$$

or in terms of orders of magnitude,

$$\frac{\mathcal{U}\_2}{l\_s} \sim \frac{V\_2}{\delta\_s} \tag{7}$$

where again the vertical velocity near *xs* is on the order of *V*2, the inviscid flow's vertical velocity component immediately downstream of the oblique shock. Thus, since *V*2/*U*<sup>2</sup> ∼ tan *θ*, we obtain the following estimate for *δs*/*ls* :

$$
\delta\_\text{s}/l\_\text{s} \sim \tan \text{ } \theta
\tag{8}
$$

Note that this relationship is analogous to one of the key assumptions underlying Chapman's (Chapman et al., 1958) free interaction model, viz, the displacement of the external inviscid flow is determined by the streamwise rate of boundary layer growth. Using (8) in (5) and solving for *δ<sup>s</sup>* finally yields an estimate for the boundary layer thickness near *xs* :

$$
\delta\_s \sim \frac{\mu M\_2 a\_2 \cos \theta}{(P\_p - P\_1) \tan \theta} \tag{9}
$$

where *U*<sup>2</sup> = *M*2*a*<sup>2</sup> cos *θ* and *a*<sup>2</sup> is the sound speed. Before proceeding, and as an aside, we rewrite the estimate in (9) as

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

outer-most streamline (i.e., a streamline in the vicinity of *δ* which is roughly parallel to the local displacement thickness) as *r*(*x*) = *f*(*x*), where *r*(*x*) is the radial distance from the nozzle

and *vs* by their approximate free stream magnitudes downstream of the oblique shock (since again, at axial position *xs*, the boundary layer has passed through the compression system

where arguments given in (Keanini & Brown, 2007) lead to the necessary pressure gradient

In order to proceed, we must estimate the magnitude of the boundary layer thickness, *δs*, immediately upstream of the separation point, *xs*. First, note that at the wall between *xi* and

where *μ* is the dynamic viscosity. Estimating the magnitude of each term in this equation

where, since we are focusing on the neighborhood of *xs*, *u* is approximated as *U*2, and where it is recognized that the streamwise pressure increases from approximately *P*<sup>1</sup> near *xi* to approximately *Pp* near *xs*. [Although *P* equals *Ps* at *xs*, due to the relatively small difference between *Ps* and *Pp*, for simplicity, we approximate *P*(*xs*) as *Pp*.] The x-length scale, *ls* = (*xs* − *xi*), is approximately equal to the length of the shock interaction zone, *lp* = *xp* − *xi*. Considering the continuity equation near *xs*, we recognize that since the boundary layer acquires a vertical velocity component as it travels toward and past *xs*, and since associated mass advection and volumetric dilatation terms, *<sup>ρ</sup>*−1**<sup>u</sup>** · ∇*ρ*, and ∇ · **<sup>u</sup>**, respectively, are of the

(*Pp* − *P*1) *ls*

∼ 1 *μ*

> *∂u <sup>∂</sup><sup>x</sup>* <sup>≈</sup> *<sup>∂</sup><sup>v</sup>*

*U*<sup>2</sup> *ls*

<sup>∼</sup> *<sup>V</sup>*<sup>2</sup> *δs*

where again the vertical velocity near *xs* is on the order of *V*2, the inviscid flow's vertical velocity component immediately downstream of the oblique shock. Thus, since *V*2/*U*<sup>2</sup> ∼

Note that this relationship is analogous to one of the key assumptions underlying Chapman's (Chapman et al., 1958) free interaction model, viz, the displacement of the external inviscid

*μ ∂*2*u <sup>∂</sup>y*<sup>2</sup> <sup>∼</sup> *<sup>∂</sup><sup>P</sup>*

*U*<sup>2</sup> *δ*2 *s*

<sup>2</sup> <sup>+</sup> *<sup>V</sup>*<sup>2</sup>

1 + (*f* �)2)−3.

(3)

(5)

(7)

*<sup>s</sup>*, and estimating *us*

*<sup>s</sup>* = *u*<sup>2</sup>

<sup>2</sup> . Replacing terms in (2) by their approximate

*<sup>s</sup>* + *v*<sup>2</sup>

*<sup>∂</sup><sup>x</sup>* (4)

*<sup>∂</sup><sup>y</sup>* (6)

*δs*/*ls* ∼ tan *θ* (8)

centerline to the streamline, evaluated at axial position *x*. Thus, *R*−<sup>1</sup> = *f* ��(

*<sup>s</sup>* in terms of local cartesian velocity components,*V*<sup>2</sup>

*<sup>s</sup>* <sup>≈</sup> *<sup>U</sup>*<sup>2</sup>

*ρ*2 *U*2 <sup>2</sup> <sup>+</sup> *<sup>V</sup>*<sup>2</sup> 2 *<sup>R</sup>* <sup>∼</sup> *Pp* <sup>−</sup> *<sup>P</sup>*<sup>2</sup> *δs*

*xs*, the x-momentum equation yields the approximate balance

Expressing *V*<sup>2</sup>

below the shock), we obtain *V*<sup>2</sup>

magnitudes then leads to

and density estimates.

leads to

same order, then

or in terms of orders of magnitude,

tan *θ*, we obtain the following estimate for *δs*/*ls* :

**2.3 Boundary layer thickness**

$$\frac{\mu \mathcal{U}\_2 / \delta\_s}{P\_p - P\_1} \sim \frac{\tau\_2}{\Delta P} \sim \tan \theta \tag{10}$$

Recognizing that the resultant stress on a fluid particle near *xs* is approximately equal to the vector sum of the horizontally acting viscous shear stress, *τ*2, and the vertically-acting net pressure, Δ*P* = *Pp* − *P*1, then (10) shows that the resultant acts in the direction of the separating boundary layer, *θ*, as it must. Since *δs*/*ls* ∼ tan *θ*, then (10) is also consistent with Chapman's (Chapman et al., 1958) estimate for *τ*2/Δ*P*. Likewise, Chapman et al. (Chapman et al., 1958) argued that *<sup>δ</sup>*∗/*ls* <sup>∼</sup> *Cf* <sup>=</sup> *<sup>τ</sup>w*/(*ρ*1*u*<sup>2</sup> <sup>1</sup>/2), where *Cf* is the friction factor; since *τ<sup>w</sup>* ∼ Δ*P*, then for nominally fixed *θ* (see Sec. II.D below), (10) is also consistent with Chapman's estimate for *Cf* .

Returning to Eqn. (3), inserting (9) and the expression for *R*−<sup>1</sup> and rearranging leads to an expression of the following form:

$$\left(\frac{P\_p}{P\_1}\right)^2 - (1+G)\frac{P\_p}{P\_1} + G(1+\epsilon) \sim 0\tag{11}$$

where *G* = *P*2/*P*1, *�* = *kM*<sup>3</sup> <sup>2</sup> *a*<sup>2</sup> cos *θμ*<sup>2</sup> *f* ��/(*P*<sup>1</sup> *f* � ), *μ*<sup>2</sup> is the gas viscosity near the separation point,and where the approximation tan *θ* ≈ −*f* � has been used. Under typical conditions, e.g., those extant in experiments described below [*Mi* <sup>≈</sup> 5, *<sup>θ</sup>* <sup>≈</sup> <sup>16</sup>*o*, *To* <sup>=</sup> 310K, *Po* <sup>=</sup> 1.24 Mpa], *�* <sup>=</sup> *<sup>O</sup>*(10−3), i.e., *�* << 1 [where *<sup>f</sup>* ��|*max* <sup>≈</sup> (*dθ*/*dx*)|*max* <sup>=</sup> *<sup>O</sup>*(1)]. Thus, solving (11) for *Pp*/*P*<sup>1</sup> and neglecting terms smaller than *O*(*�*), we finally obtain

$$\frac{P\_p}{P\_1} \sim G - \varepsilon \frac{2G}{G-1} \tag{12}$$

Importantly, this equation shows that *P*1/*Pp* ≈ *Pi*/*Pp* ∼ *P*1/*P*2, *demonstrating that the separation pressure ratio essentially corresponds to the oblique shock pressure ratio.*

#### **2.4 Shock and flow deflection angles**

In order to close the approximate model embodied in (12), it is necessary to specify the shock angle, *β*, and the flow deflection angle, *θ*. Referring to earlier work, Summerfield et al. (Summerfield et al., 1954) used measured separation pressure ratios and Mach numbers in the oblique shock relations to infer *θ*; based on their data, they inferred a nominally fixed value, *<sup>θ</sup>* <sup>≈</sup> <sup>16</sup>*o*. By contrast, Ostlund (Ostlund, 2002), again using the same approach,argued that *θ* varies with *Mi*, albeit weakly; he fit his estimate with a linear relationship, *θ* = 1.678*Mi* + 9.347, valid for 2.5 ≤ *Mi* ≤ 4.5. For this range of *Mi*, however, the correlation indicates that *θ* only varies from 13.5*<sup>o</sup>* and 16.9*o*. Based on these indirect estimates, we assume that *θ* is constant; for simplicity, we will arbitrarily adopt the average value of *θ* indicated by Ostlund's correlation, *<sup>θ</sup>* <sup>≈</sup> 15.2*o*, nearly equal to Summerfield's (Summerfield et al., 1954)

Rocket Ascent 9

<sup>163</sup> Shock-Induced Turbulent Boundary Layer Separation in Over-Expanded Rocket Nozzles: Physics, Models, Random Side Loads, and the Diffusive Character of Stochastic Rocket Ascent

absence of separation. For a range of pressure gradients observed in a number of different

from either (18) or (19); in the latter case, we follow Chapman (Chapman et al., 1958) and

)*kM*<sup>2</sup> *i* 

) (or *F*(*s*)), the predicted separation pressure ratio, *Pi*/*Pp*, can be determined

 *Cf i*

) at the effective separation point, *sp*[= (*xp* − *xi*)/(*xs* − *xi*)]. It is important to note

*<sup>i</sup>* − <sup>1</sup>)1/4

 + 1 −<sup>1</sup>

) [= 6.0; see (Ostlund, 2002)] is the value

(20)

<sup>√</sup>2(*M*<sup>2</sup>

that (Ostlund, 2002) has developed an alternative separation criterion which requires *a priori* specification of both the plateau pressure, *Pp*, and the friction coefficient, *Cf i* at *xi*. The criterion in (20) by contrast only requires information on *Cf i*. Fortuitously, and as originally shown by Chapman (Chapman et al., 1958), for shock-induced separation of *turbulent* boundary layers, the dependence of *Pi*/*Pp* on *Cf i* (or equivalently, on *Reδ*∗, the displacement thickness Reynolds number at *xi*) is weak, at least over the range of Mach numbers investigated (1.3 ≤ *Mi* ≤ 4.0), consistent with both the scale analysis above and previously developed correlations; see

A series of experiments were carried out in the Nozzle Test Facility at Marshall Space Flight Center (Keanini & Brown, 2007). The experiments were designed to investigate the role of boundary layer separation on nozzle side- loading and to examine fluid-solid interactions

A sub-scale, ideal-contour nozzle, having an area ratio of approximately 30:1 (exit to throat area) was operated under a range of cold-flow, overexpanded conditions. The nozzle was outfitted with a series of pressure taps, where tap spacing in the axial direction was 0.0254 m. Two sets of azimuthally spaced taps were also used, placed at two axial locations, at 45*<sup>o</sup>* intervals around the nozzle circumference. The axially-spaced taps allowed measurement of the instantaneous and time-average axial pressure distribution within the nozzle while the azimuthally distributed taps allowed examination of the instantaneous and time-average separation line (under conditions where the shock interaction zone coincides with either set of azimuthally distributed taps). The throat diameter was 0.0254 m and the design Mach number was 5.25. Pressures at all taps were sampled at 10 kHz, sufficiently high to allow study of the low-frequency, large-amplitude component of shock motion (Dolling & Brusniak, 1989; Keanini & Brown, 2007), but not sufficient to resolve small-scale, high frequency jitter.

The scaling relationship in (12) was fit to available data on separation in overexpanded nozzles. Data on free shock separation was obtained from a number of sources (Keanini & Brown, 2007), and represents flow in a variety of nozzle geometries, in both full-scale and sub-scale models, under both cold flow and hot fire conditions. Although the working fluid in most experiments was air, Bloomer's (Bloomer et al., 1961) hot fire measurements, which used a mixture of JP-4 rocket fuel and liquid oxygen (*k* = 1.2), are included in estimating best fit parameters. This approach is allowable due to the weak dependence between separation pressure ratio and *k*; see review in (Keanini & Brown, 2007). Since solutions for *P*1/*P*<sup>2</sup> at a

), are nearly identical (Ostlund, 2002).

nozzles, the two correlations, *F*(*s*) and *F*(*s*; *p*�

*Pi Pp* = *F*(*sp*; *p*�

where *ν*¯(*s*) is approximated as *ν<sup>i</sup>* and where *F*(*sp*; *p*�

(Keanini & Brown, 2007) for further details.

underlying oscillatory modes observed in side-loaded nozzles.

**2.6 Experimental measurements**

Given *F*(*s*; *p*�

of *F*(*s*; *p*�

**2.7 Results**

linearize (19) to obtain

estimate. Second, we follow Summerfield (Summerfield et al., 1954) and assume that the oblique shock relation

$$\tan \theta = \frac{2 \cot \beta (M\_1^2 \sin^2 \beta - 1)}{(k+1)M\_1^2 - 2(M\_1^2 \sin^2 \beta - 1)} \tag{13}$$

applies to the inviscid flow outside the separating boundary layer.

#### **2.5 Free interaction model**

Chapman's (Chapman et al., 1958) original analysis posited that thickening of the boundary layer displacement thickness displaced the inviscid flow above the boundary layer according to

$$P(x) - P\_i = \frac{\rho\_i u\_i^2}{\sqrt{M\_i^2 - 1}} d\theta \tag{14}$$

where *dδ*∗/*dx* = *dθ*, and where the subscripts refer to conditions at *xi*. He then estimated terms in the balance between the axial pressure gradient and cross-stream shear stress gradient as

$$\frac{P - P\_{\text{i}}}{I\_{\text{s}}} \sim \frac{\tau\_{\text{wi}}}{\delta^\*} \tag{15}$$

(where *τwi* is the wall shear stress near *xi*), then estimated *dδ*∗/*dx* as *δ*∗/*ls* (where *δ*<sup>∗</sup> ∼ *δs*), and finally combined and linearized (for small *P* − *Pi*) to obtain

$$\frac{P - P\_{\dot{i}}}{q\_{\dot{i}}} \sim \frac{\sqrt{C\_{fi}}}{(M\_{\dot{i}}^2 - 1)^{1/4}}\tag{16}$$

where *qi* = *ρiu*<sup>2</sup> *<sup>i</sup>* /2 and *Cf i* = *τwi*/*qi*. Dividing the left side of (16) by the right suggests that

$$f\left(\frac{P - P\_i}{q\_i}\right)\left(\frac{(M\_i^2 - 1)^{1/4}}{\sqrt{C\_{fi}}}\right) \sim f(\mathbf{x} - \mathbf{x}\_i) \tag{17}$$

i.e., that the term on the left depends only on position within the shock interaction zone. Erdos and Pallone (Erdos & Pallone, 1962) exploited this idea to develop a wall pressure correlation, *F*(*s*), which describes the self-similar pressure variation over the shock interaction zone, where

$$F(s) = \left(\frac{P - P\_i}{q\_i}\right) \left(\frac{(M\_i^2 - 1)^{1/4}}{\sqrt{2C\_{fi}}}\right) \tag{18}$$

and where *s* = (*x* − *xi*)/(*xs* − *xi*). Carriere et al. (Carriere et al., 1968) extended this work by developing a generalized version of (18), suitable for the non- uniform flows in nozzles. In this case, the self-similar pressure variation over the separation zone is described by

$$F(s; p') = \sqrt{\left(\frac{P - P\_i}{q\_i}\right) \left(\frac{\bar{\nu}(s) - \nu(s)}{\sqrt{C\_{fi}}}\right)}\tag{19}$$

where *p*� = (*δ*∗ *<sup>i</sup>* /*qi*)(*dP*/*dx*) is the normalized inviscid flow pressure gradient immediately upstream of *xi*, *ν*(*s*) is the Prandtl-Meyer function, and *ν*¯(*s*) is the value of the function in the absence of separation. For a range of pressure gradients observed in a number of different nozzles, the two correlations, *F*(*s*) and *F*(*s*; *p*� ), are nearly identical (Ostlund, 2002).

Given *F*(*s*; *p*� ) (or *F*(*s*)), the predicted separation pressure ratio, *Pi*/*Pp*, can be determined from either (18) or (19); in the latter case, we follow Chapman (Chapman et al., 1958) and linearize (19) to obtain

$$\frac{P\_i}{P\_p} = \left[ F(s\_p; p') k M\_i^2 \left( \frac{\sqrt{C\_{fi}}}{\sqrt{2} (M\_i^2 - 1)^{1/4}} \right) + 1 \right]^{-1} \tag{20}$$

where *ν*¯(*s*) is approximated as *ν<sup>i</sup>* and where *F*(*sp*; *p*� ) [= 6.0; see (Ostlund, 2002)] is the value of *F*(*s*; *p*� ) at the effective separation point, *sp*[= (*xp* − *xi*)/(*xs* − *xi*)]. It is important to note that (Ostlund, 2002) has developed an alternative separation criterion which requires *a priori* specification of both the plateau pressure, *Pp*, and the friction coefficient, *Cf i* at *xi*. The criterion in (20) by contrast only requires information on *Cf i*. Fortuitously, and as originally shown by Chapman (Chapman et al., 1958), for shock-induced separation of *turbulent* boundary layers, the dependence of *Pi*/*Pp* on *Cf i* (or equivalently, on *Reδ*∗, the displacement thickness Reynolds number at *xi*) is weak, at least over the range of Mach numbers investigated (1.3 ≤ *Mi* ≤ 4.0), consistent with both the scale analysis above and previously developed correlations; see (Keanini & Brown, 2007) for further details.

### **2.6 Experimental measurements**

A series of experiments were carried out in the Nozzle Test Facility at Marshall Space Flight Center (Keanini & Brown, 2007). The experiments were designed to investigate the role of boundary layer separation on nozzle side- loading and to examine fluid-solid interactions underlying oscillatory modes observed in side-loaded nozzles.

A sub-scale, ideal-contour nozzle, having an area ratio of approximately 30:1 (exit to throat area) was operated under a range of cold-flow, overexpanded conditions. The nozzle was outfitted with a series of pressure taps, where tap spacing in the axial direction was 0.0254 m. Two sets of azimuthally spaced taps were also used, placed at two axial locations, at 45*<sup>o</sup>* intervals around the nozzle circumference. The axially-spaced taps allowed measurement of the instantaneous and time-average axial pressure distribution within the nozzle while the azimuthally distributed taps allowed examination of the instantaneous and time-average separation line (under conditions where the shock interaction zone coincides with either set of azimuthally distributed taps). The throat diameter was 0.0254 m and the design Mach number was 5.25. Pressures at all taps were sampled at 10 kHz, sufficiently high to allow study of the low-frequency, large-amplitude component of shock motion (Dolling & Brusniak, 1989; Keanini & Brown, 2007), but not sufficient to resolve small-scale, high frequency jitter.

#### **2.7 Results**

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

estimate. Second, we follow Summerfield (Summerfield et al., 1954) and assume that the

Chapman's (Chapman et al., 1958) original analysis posited that thickening of the boundary layer displacement thickness displaced the inviscid flow above the boundary layer according

where *dδ*∗/*dx* = *dθ*, and where the subscripts refer to conditions at *xi*. He then estimated terms in the balance between the axial pressure gradient and cross-stream shear stress

(where *τwi* is the wall shear stress near *xi*), then estimated *dδ*∗/*dx* as *δ*∗/*ls* (where *δ*<sup>∗</sup> ∼ *δs*),

(*M*<sup>2</sup>

*<sup>i</sup>* <sup>−</sup> <sup>1</sup>)1/4 *Cf i*

i.e., that the term on the left depends only on position within the shock interaction zone. Erdos and Pallone (Erdos & Pallone, 1962) exploited this idea to develop a wall pressure correlation, *F*(*s*), which describes the self-similar pressure variation over the shock interaction

and where *s* = (*x* − *xi*)/(*xs* − *xi*). Carriere et al. (Carriere et al., 1968) extended this work by developing a generalized version of (18), suitable for the non- uniform flows in nozzles. In

> *P* − *Pi qi*

upstream of *xi*, *ν*(*s*) is the Prandtl-Meyer function, and *ν*¯(*s*) is the value of the function in the

(*M*<sup>2</sup>

*<sup>i</sup>* <sup>−</sup> <sup>1</sup>)1/4 2*Cf i*

 *ν*¯(*s*) − *ν*(*s*) *Cf i*

*<sup>i</sup>* /*qi*)(*dP*/*dx*) is the normalized inviscid flow pressure gradient immediately

<sup>∼</sup> *<sup>τ</sup>wi*

 *Cf i*

*<sup>P</sup>*(*x*) <sup>−</sup> *Pi* <sup>=</sup> *<sup>ρ</sup>iu*<sup>2</sup> 

> *P* − *Pi ls*

> > ∼

*P* − *Pi qi*

(*M*<sup>2</sup>

 *P* − *Pi qi*

this case, the self-similar pressure variation over the separation zone is described by

<sup>1</sup> <sup>−</sup> <sup>2</sup>(*M*<sup>2</sup>

<sup>1</sup> sin2 *<sup>β</sup>* <sup>−</sup> <sup>1</sup>)

*i M*<sup>2</sup> *<sup>i</sup>* − 1

<sup>1</sup> sin<sup>2</sup> *<sup>β</sup>* <sup>−</sup> <sup>1</sup>) (13)

*dθ* (14)

*<sup>δ</sup>*<sup>∗</sup> (15)

*<sup>i</sup>* <sup>−</sup> <sup>1</sup>)1/4 (16)

∼ *f*(*x* − *xi*) (17)

(18)

(19)

tan *<sup>θ</sup>* <sup>=</sup> 2 cot *<sup>β</sup>*(*M*<sup>2</sup>

applies to the inviscid flow outside the separating boundary layer.

and finally combined and linearized (for small *P* − *Pi*) to obtain

 *P* − *Pi qi*

*F*(*s*) =

*F*(*s*; *p*�

) = 

*<sup>i</sup>* /2 and *Cf i* = *τwi*/*qi*. Dividing the left side of (16) by the right suggests that

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

oblique shock relation

**2.5 Free interaction model**

to

gradient as

where *qi* = *ρiu*<sup>2</sup>

zone, where

where *p*� = (*δ*∗

The scaling relationship in (12) was fit to available data on separation in overexpanded nozzles. Data on free shock separation was obtained from a number of sources (Keanini & Brown, 2007), and represents flow in a variety of nozzle geometries, in both full-scale and sub-scale models, under both cold flow and hot fire conditions. Although the working fluid in most experiments was air, Bloomer's (Bloomer et al., 1961) hot fire measurements, which used a mixture of JP-4 rocket fuel and liquid oxygen (*k* = 1.2), are included in estimating best fit parameters. This approach is allowable due to the weak dependence between separation pressure ratio and *k*; see review in (Keanini & Brown, 2007). Since solutions for *P*1/*P*<sup>2</sup> at a

Rocket Ascent 11

<sup>165</sup> Shock-Induced Turbulent Boundary Layer Separation in Over-Expanded Rocket Nozzles: Physics, Models, Random Side Loads, and the Diffusive Character of Stochastic Rocket Ascent

> **Arens Farley J−2S Dumnov**

**Bloomer (k=1.2) Lawrence Model II**

**1.5 2 2.5 3 3.5 4 4.5 5 5.5**

**Mi**

Fig. 5. Comparison of Model II (Keanini & Brown, 2007) with separation measurements in

In order to use the free interaction separation criterion in (20), *Cf i* must be specified. As noted, and based on Chapman's (Chapman et al., 1958) observation that *Pi*/*Pp* is weakly dependent on *Cf i*, we assume that *Cf i* is constant. The assumed magnitude, *Cf i* = 0.00245, represents the characteristic value obtained from fitting the free interaction model to observed time-average shock interaction zone pressure variations, as described in (Keanini & Brown, 2007). In addition, this value is used in fitting the generalized quasi-one-dimensional flow model, necessary for computing the flow upstream of the separation zone, to our experimental

As shown in figure 6, over 1.75 ≤ *Mi* ≤ 5.5, predicted separation pressure ratios obtained via the free interaction model are quite similar to those obtained via the scaling analysis above. Given the reasonable agreement between model predictions and previous observations, this result simplifies Ostlund's (Ostlund, 2002) separation criterion by eliminating the need for *a priori* specification of *Pp*. Importantly, this result provides further evidence of the applicability of the free interaction model to separation in nozzles, and moreover, further indicates the

The models described in the previous section allow determination of the mean separation line position, *xs*(*t*) = *xs*(*H*(*t*)), within the nozzle, as a function of the instantaneous nozzle pressure ratio *NPR* = *NPR*(*H*(*t*)) = *Po*(*t*)/*Pa*(*H*(*t*)) (Keanini & Brown, 2007; Keanini et al., 2011; Srivastava et al., 2010), where *Po*(*t*) is the time-dependent combustion chamber pressure

**0.1**

rocket engine nozzles. The fitting constant equals 1.14.

**2.8 Separation pressure ratios via the free interaction model**

shock-free flow measurements; see (Keanini & Brown, 2007) for details.

**3. Physically consistent models of random nozzle side loads**

physical consistency of the scale analyses presented in (Keanini & Brown, 2007).

**0.2**

**0.3**

**0.4**

**Pi / P**

**p**

**0.5**

**0.6**

**0.7**

**0.8**

turning angle of *<sup>θ</sup>* <sup>=</sup> 15.2*<sup>o</sup>* do not exist for *Mi* <sup>&</sup>lt;<sup>≈</sup> 1.7, only data obtained at *Mi* <sup>≥</sup> 1.75 are used in the fitting procedure. For comparative purposes, however, the limited data available at *Mi* ≤ 1.75 are presented in the graphs below.

A comparison of separation pressure ratios predicted by the first scaling model (Keanini & Brown, 2007), model I in figures 4 and 5, with available data, shown in figure 4, indicates that the model provides reasonable predictions over the range 1.75 ≤ *Mi* ≤ 4.0. The least square fitting constant is found to be 1.52 . A similar comparison using model II, Eq. (12), shown in figure 5, likewise indicates reasonable agreement over 1.75 ≤ *Mi* ≤ 4.5, with significantly improved agreement for *Mi* > 4.0; the fitting constant in this case is 1.14. Comparing with Ostlund (Ostlund, 2002) and Frey's (Frey & Hagemann, 1998) correlations, which fit observed shock pressure ratios to the oblique shock pressure ratio (based on inferred shock and deflection angles), we note that their quoted ranges of validity were in both cases 2.5 ≤ *Mi* ≤ 4.5. At higher Mach numbers (*Mi* ≥ 4), the data suggests that *Pi*/*Pp* becomes

Fig. 4. Comparison of Model I in (Keanini & Brown, 2007) with separation measurements in rocket engine nozzles. The fitting constant equals 1.52.

largely independent of *Mi*. Although constancy of *Pi*/*Pp* is not inconsistent with separation remaining dominated by the oblique shock, since the asymptotic expression for *Pi*/*Pp* at large *Mi* is, from (12),

$$\frac{P\_i}{P\_p} \sim \frac{k+1}{2k\sin^2\beta M\_i^2}$$

then due to an 82 % variation in *M*<sup>2</sup> *<sup>i</sup>* over 4.0 ≤ *Mi* ≤ 5.4, the time-average deflection angle, *θ*, likely becomes moderately dependent on *Mi*.

Fig. 5. Comparison of Model II (Keanini & Brown, 2007) with separation measurements in rocket engine nozzles. The fitting constant equals 1.14.

#### **2.8 Separation pressure ratios via the free interaction model**

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

turning angle of *<sup>θ</sup>* <sup>=</sup> 15.2*<sup>o</sup>* do not exist for *Mi* <sup>&</sup>lt;<sup>≈</sup> 1.7, only data obtained at *Mi* <sup>≥</sup> 1.75 are used in the fitting procedure. For comparative purposes, however, the limited data available

A comparison of separation pressure ratios predicted by the first scaling model (Keanini & Brown, 2007), model I in figures 4 and 5, with available data, shown in figure 4, indicates that the model provides reasonable predictions over the range 1.75 ≤ *Mi* ≤ 4.0. The least square fitting constant is found to be 1.52 . A similar comparison using model II, Eq. (12), shown in figure 5, likewise indicates reasonable agreement over 1.75 ≤ *Mi* ≤ 4.5, with significantly improved agreement for *Mi* > 4.0; the fitting constant in this case is 1.14. Comparing with Ostlund (Ostlund, 2002) and Frey's (Frey & Hagemann, 1998) correlations, which fit observed shock pressure ratios to the oblique shock pressure ratio (based on inferred shock and deflection angles), we note that their quoted ranges of validity were in both cases 2.5 ≤ *Mi* ≤ 4.5. At higher Mach numbers (*Mi* ≥ 4), the data suggests that *Pi*/*Pp* becomes

> **Arens Farley J−2S Dumnov**

**Bloomer (k=1.2) Lawrence Model I**

**1.5 2 2.5 3 3.5 4 4.5 5 5.5**

**Mi**

Fig. 4. Comparison of Model I in (Keanini & Brown, 2007) with separation measurements in

largely independent of *Mi*. Although constancy of *Pi*/*Pp* is not inconsistent with separation remaining dominated by the oblique shock, since the asymptotic expression for *Pi*/*Pp* at large

> <sup>∼</sup> *<sup>k</sup>* <sup>+</sup> <sup>1</sup> 2*k* sin2 *βM*<sup>2</sup>

*i*

*<sup>i</sup>* over 4.0 ≤ *Mi* ≤ 5.4, the time-average deflection angle, *θ*,

*Pi Pp*

at *Mi* ≤ 1.75 are presented in the graphs below.

**0.1**

then due to an 82 % variation in *M*<sup>2</sup>

likely becomes moderately dependent on *Mi*.

*Mi* is, from (12),

rocket engine nozzles. The fitting constant equals 1.52.

**0.2**

**0.3**

**0.4**

**0.5**

**Pi / Pp**

**0.6**

**0.7**

**0.8**

In order to use the free interaction separation criterion in (20), *Cf i* must be specified. As noted, and based on Chapman's (Chapman et al., 1958) observation that *Pi*/*Pp* is weakly dependent on *Cf i*, we assume that *Cf i* is constant. The assumed magnitude, *Cf i* = 0.00245, represents the characteristic value obtained from fitting the free interaction model to observed time-average shock interaction zone pressure variations, as described in (Keanini & Brown, 2007). In addition, this value is used in fitting the generalized quasi-one-dimensional flow model, necessary for computing the flow upstream of the separation zone, to our experimental shock-free flow measurements; see (Keanini & Brown, 2007) for details.

As shown in figure 6, over 1.75 ≤ *Mi* ≤ 5.5, predicted separation pressure ratios obtained via the free interaction model are quite similar to those obtained via the scaling analysis above. Given the reasonable agreement between model predictions and previous observations, this result simplifies Ostlund's (Ostlund, 2002) separation criterion by eliminating the need for *a priori* specification of *Pp*. Importantly, this result provides further evidence of the applicability of the free interaction model to separation in nozzles, and moreover, further indicates the physical consistency of the scale analyses presented in (Keanini & Brown, 2007).

#### **3. Physically consistent models of random nozzle side loads**

The models described in the previous section allow determination of the mean separation line position, *xs*(*t*) = *xs*(*H*(*t*)), within the nozzle, as a function of the instantaneous nozzle pressure ratio *NPR* = *NPR*(*H*(*t*)) = *Po*(*t*)/*Pa*(*H*(*t*)) (Keanini & Brown, 2007; Keanini et al., 2011; Srivastava et al., 2010), where *Po*(*t*) is the time-dependent combustion chamber pressure

Rocket Ascent 13

<sup>167</sup> Shock-Induced Turbulent Boundary Layer Separation in Over-Expanded Rocket Nozzles: Physics, Models, Random Side Loads, and the Diffusive Character of Stochastic Rocket Ascent

Constituent displacements in the set of *N* displacements are assumed independent, and based on experimental observations (Dolling & Brusniak, 1989), gaussian. Thus, each *pI* is given by

exp

− *s*2 *I* 2*σ*<sup>2</sup> *s* 

*<sup>s</sup> δ*(*φ* − *φ*�

φ

**I**

**R(H(t))**

**2Ls**

Δφ

**Fs**(*t*) = **Fr**(*t*) + **Fx**(*t*) (24)

(22)

) (23)

<sup>2</sup>*πσ*<sup>2</sup> *s*

In moving to a continuous description of the separation line, (Srivastava et al., 2010) assumes

, *t*) >= *σ*<sup>2</sup>

Considering next the side load, we express the instantaneous force vector produced by

s **(H(t))x**

> **s(**φ**<sup>I</sup> ,t)**

Fig. 7. Model I (Srivastava et al., 2010) separation line model. The mean separation line position, *xs*(*H*(*t*)), moves down the nozzle axis, on the slow time scale associated with vertical rocket motion. By contrast, axial separation line motion about *xs*(*H*(*t*)), at any angular position, *φI*, is random, and takes place on a much shorter time scale. Rapid axial motion, in addition, is confined to the nominal shock-boundary layer interaction zone, again denoted by *Ls*. Pressures upstream and downstream of the instantaneous separation line, *Pi* = *Pi*(*H*(*t*)) and *P*<sup>2</sup> = *P*2(*H*(*t*)), respectively, are assumed to be spatially uniform within

asymmetric boundary layer separation, **Fs**(*t*), as a sum of radial and axial components

 = *σ*2,

In (Srivastava et al., 2010), the following *ad hoc* side load model was assumed:

*Fsz*<sup>−</sup> <sup>&</sup>lt; *Fsz* <sup>&</sup>gt; <sup>2</sup>

A) *Fsy* and *Fsz* are independent, gaussian random variables,

*pI*(*sI*) = <sup>1</sup>

< *s*(*φ*, *t*)*s*(*φ*�

φ

*Ls*. Adapted from (Srivastava et al., 2010).

 = 

B) < *Fsy* >= 0 and < *Fsz* >= 0,

*Fsy*<sup>−</sup> <sup>&</sup>lt; *Fsy* <sup>&</sup>gt; <sup>2</sup>

C)

**x**

that

Fig. 6. Comparison of free interaction model and model II (Keanini & Brown, 2007) with separation measurements in rocket engine nozzles.

and *Pa*(*H*(*t*)) is the altitude- dependent ambient pressure. Given *xs*(*t*), physically consistent models that describe the stochastic evolution of the separation line *shape* relative to *xs*(*t*), can be developed. This section first highlights the essential elements of the stochastic side load models developed in (Keanini et al., 2011; Srivastava et al., 2010) and then briefly outlines a new, physically consistent model of stochastic separation line evolution (Keanini et al., 2011). As detailed in (Keanini et al., 2011; Srivastava et al., 2010), and in response to the decaying altitude-dependent ambient pressure, the mean position of boundary layer separation line, *xs*(*t*), travels down the nozzle axis toward the nozzle exit, with motion taking place on a relatively slow time scale, *τ<sup>a</sup>* = Δ*xa*/*VR*, where Δ*xa* is the characteristic incremental altitude over which ambient pressure varies and *VR* is the characteristic rocket speed. Superposed on this slow motion is a fast, random, azimuthally homogeneous stochastic motion. Following (Keanini et al., 2011; Srivastava et al., 2010), the joint probability density, *ps*, associated with the instantaneous random separation line shape is given by

$$p\_s(s\_1, s\_2, \dots, s\_N) = \prod\_I p\_I = \frac{1}{(2\pi\sigma\_s)^{N/2}} \exp\left[-\frac{s\_1^2 + s\_2^2 + s\_3^2 + \dots + s\_N^2}{2\sigma\_s^2}\right] \tag{21}$$

where, as shown in figure 7, *sI* is the random axial displacement of the separation line at azimuthal angle *φI*, and *σ*<sup>2</sup> *<sup>s</sup>* is the (assumed) constant variance of local separation line displacements.

Constituent displacements in the set of *N* displacements are assumed independent, and based on experimental observations (Dolling & Brusniak, 1989), gaussian. Thus, each *pI* is given by

$$p\_I(s\_I) = \frac{1}{\sqrt{2\pi\sigma\_s^2}} \exp\left[-\frac{s\_I^2}{2\sigma\_s^2}\right] \tag{22}$$

In moving to a continuous description of the separation line, (Srivastava et al., 2010) assumes that

$$ = \sigma\_s^2 \delta(\phi - \phi')\tag{23}$$

Considering next the side load, we express the instantaneous force vector produced by

Fig. 7. Model I (Srivastava et al., 2010) separation line model. The mean separation line position, *xs*(*H*(*t*)), moves down the nozzle axis, on the slow time scale associated with vertical rocket motion. By contrast, axial separation line motion about *xs*(*H*(*t*)), at any angular position, *φI*, is random, and takes place on a much shorter time scale. Rapid axial motion, in addition, is confined to the nominal shock-boundary layer interaction zone, again denoted by *Ls*. Pressures upstream and downstream of the instantaneous separation line, *Pi* = *Pi*(*H*(*t*)) and *P*<sup>2</sup> = *P*2(*H*(*t*)), respectively, are assumed to be spatially uniform within *Ls*. Adapted from (Srivastava et al., 2010).

asymmetric boundary layer separation, **Fs**(*t*), as a sum of radial and axial components

$$\mathbf{F\_s}(t) = \mathbf{F\_r}(t) + \mathbf{F\_x}(t) \tag{24}$$

In (Srivastava et al., 2010), the following *ad hoc* side load model was assumed:

A) *Fsy* and *Fsz* are independent, gaussian random variables,

B) < *Fsy* >= 0 and < *Fsz* >= 0,

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

**Arens Farley J−2S Dumnov**

**Model II**

**Bloomer (k=1.2) Lawrence Free Int. Model**

**1 2 3 4 5 6**

**Mi**

Fig. 6. Comparison of free interaction model and model II (Keanini & Brown, 2007) with

and *Pa*(*H*(*t*)) is the altitude- dependent ambient pressure. Given *xs*(*t*), physically consistent models that describe the stochastic evolution of the separation line *shape* relative to *xs*(*t*), can be developed. This section first highlights the essential elements of the stochastic side load models developed in (Keanini et al., 2011; Srivastava et al., 2010) and then briefly outlines a new, physically consistent model of stochastic separation line evolution (Keanini et al., 2011). As detailed in (Keanini et al., 2011; Srivastava et al., 2010), and in response to the decaying altitude-dependent ambient pressure, the mean position of boundary layer separation line, *xs*(*t*), travels down the nozzle axis toward the nozzle exit, with motion taking place on a relatively slow time scale, *τ<sup>a</sup>* = Δ*xa*/*VR*, where Δ*xa* is the characteristic incremental altitude over which ambient pressure varies and *VR* is the characteristic rocket speed. Superposed on this slow motion is a fast, random, azimuthally homogeneous stochastic motion. Following (Keanini et al., 2011; Srivastava et al., 2010), the joint probability density, *ps*, associated with

**0.1**

separation measurements in rocket engine nozzles.

the instantaneous random separation line shape is given by

*I*

*pI* <sup>=</sup> <sup>1</sup>

(2*πσs*)*N*/2 exp

where, as shown in figure 7, *sI* is the random axial displacement of the separation line

− *s*2 <sup>1</sup> <sup>+</sup> *<sup>s</sup>*<sup>2</sup> <sup>2</sup> <sup>+</sup> *<sup>s</sup>*<sup>2</sup>

*<sup>s</sup>* is the (assumed) constant variance of local separation line

<sup>3</sup> <sup>+</sup> ... <sup>+</sup> *<sup>s</sup>*<sup>2</sup>

2*σ*<sup>2</sup> *s*

*N*

(21)

*ps*(*s*1,*s*2,...,*sN*) = ∏

at azimuthal angle *φI*, and *σ*<sup>2</sup>

displacements.

**0.2**

**0.3**

**0.4**

**0.5**

**Pi/Pp**

**0.6**

**0.7**

**0.8**

C) *Fsy*<sup>−</sup> <sup>&</sup>lt; *Fsy* <sup>&</sup>gt; <sup>2</sup> = *Fsz*<sup>−</sup> <sup>&</sup>lt; *Fsz* <sup>&</sup>gt; <sup>2</sup> = *σ*2,

Rocket Ascent 15

<sup>169</sup> Shock-Induced Turbulent Boundary Layer Separation in Over-Expanded Rocket Nozzles: Physics, Models, Random Side Loads, and the Diffusive Character of Stochastic Rocket Ascent

*<sup>σ</sup>*<sup>2</sup> exp

Theoretical determination of rocket response to side loads requires that the time correlation

this subsection, we propose(and physically justify) that local separation line dynamics can be modeled as an Ornstein-Uhlenbeck process. Once this assumption is made, then the second step rests on a rigorous argument showing that on the relatively long rocket dynamics time scale, the boundary layer separation line shape, and importantly, associated side load

�

<sup>−</sup> *<sup>A</sup>*<sup>2</sup> 2*σ*<sup>2</sup> 

)*Fsα*(*t*)�, be first determined. As detailed

)*Fsα*(*t*)� is developed in two steps. First, and as detailed in

*dsi*(*t*) = <sup>−</sup>*ksi*(*t*) + <sup>√</sup>*DsdW*(*t*) (31)

(30)

is the uniform probability density underlying the random direction *φs*, and

is the Rayleigh distribution for the amplitude *A*.

function for either side load component, �*Fsα*(*t*

in (Keanini et al., 2011), �*Fsα*(*t*

**3.2 Ornstein-Uhlenbeck model of separation line dynamics**

�

stochastic rotational rocket dynamics (Keanini et al., 2011).

e.g., (Dolling & Brusniak, 1989).

(Plotkin, 1975).

recirculation zone.

shock-separated flows near compression corners and blunt fins:

*pA*(*A*) = *<sup>A</sup>*

components, are all delta correlated in time. See (Keanini et al., 2011)for details.

The following simple, explicit stochastic model of separation line dynamics is proposed:

where *si* (*t*) = *s*(*φi*, *t*) is the instantaneous separation line position at *φi*, *k* and *Ds* are damping and effective diffusion coefficients, and *dW*(*t*) is a differential Weiner process. This equation, describing an Ornstein-Uhlenbeck process, allows straightforward, physically consistent calculation of statistical properties associated with separation line motion and, more importantly, serves as the first link in a chain that connects short-time-scale random separation line motion to short-time-scale random side loads, and in turn, to long-time-scale

The form of this equation is chosen based on the following experimental features, observed in

a) Under statistically stationary conditions, the feet of separation-inducing shocks oscillate randomly, up- and downstream, over limited distances, about a fixed mean position; see,

b) As observed in (Dolling & Brusniak, 1989)the distribution of shock foot positions within

c) The time correlation of shock foot positions, as indicated by wall pressure measurements within the shock-boundary layer interaction zone, decays rapidly for time intervals, �*t*, larger than a short correlation time, *τs*, a feature that can be inferred, for example, from

Physically, the damping term captures the fact that the shock sits within a pressure-potential energy well. Thus, downstream shock excursions incrementally decrease and increase, respectively, upstream and downstream shock face pressures; the resulting pressure imbalance forces the shock back upstream. A similar mechanism operates during upstream excursions. Introduction of a Weiner process models the combined random forcing produced by advection of turbulent boundary layer structures through the upstream side of the shock foot and pressure oscillations emanating from the downstream separated boundary layer and

We note that the proposed model is qualitatively consistent with Plotkin's model of boundary layer-driven shock motion near compression corners and blunt fins (Plotkin, 1975). Plotkin's

the shock-boundary layer interaction zone is approximately gaussian.

where, assuming ergodicity, < · > denotes either an ensemble or time average, and where the separation line model above is used to calculate the force variance *σ*2. The side load components *Fsy* and *Fsz* are expressed with respect to rocket-fixed coordinates; see (Keanini et al., 2011; Srivastava et al., 2010).

In order to demonstrate the physical consistency of the model introduced in (Srivastava et al., 2010), (Keanini et al., 2011) first shows that the *assumed* properties, A) -C), can be *derived* from the simple *separation line model* developed in (Srivastava et al., 2010). Second, and as shown in the next subsection, the side load model in A) - C) then leads to *experimentally observed* side load amplitude and direction densities (Keanini et al., 2011).

#### **3.1 Derivation of density functions for side load amplitude and direction**

Two important experimental and numerical observations concerning the side load, **Fs** (within rigid, axisymmetric nozzles) are first noted:


Both observations can be *derived*, starting from the simple statistical model of random side loads, A) - C), immediately above. Thus, given *A* and *φs*, the instantaneous side load components in body-fixed *y* and *z* directions are given by

$$F\_{sy} = A \cos \phi\_s \qquad \qquad \qquad F\_{sz} = A \sin \phi\_s$$

Following (Srivastava et al., 2010), write *Fsy* and *Fsz* as *Fsy* = *Y*¯ = *A* cos *φ<sup>s</sup>* and *Fsz* = *Z*¯ = *A* sin *φs*; thus, the joint probability density associated with *Fsy* and *Fsz* can be expressed as

$$p\_{\bar{Y}\bar{Z}}(\bar{Y},\bar{Z}) = p\_{\bar{Y}}(\bar{Y})p\_{\bar{Z}}(\bar{Z}) = \frac{1}{2\pi\sigma^2} \exp\left(-\frac{\bar{Y}^2 + \bar{Z}^2}{2\sigma^2}\right) \tag{25}$$

Following (Srivastava et al., 2010), we restate *pY*¯*Z*¯ in terms of *A* and *φ<sup>s</sup>* as,

$$p\_{A\phi\_{\bar{\*}}} = |I| p\_{\bar{\Upsilon}\mathcal{Z}}(\bar{Y}, \bar{Z}) \tag{26}$$

where *pA<sup>φ</sup><sup>s</sup>* (*A*, *φs*) is the joint pdf for the random amplitude and direction of **Fr**, and where the jacobian determinant is given by

$$||J|| = \begin{vmatrix} \frac{\partial \bar{Y}}{\partial A} & \frac{\partial \bar{Y}}{\partial \phi\_\*}\\ \frac{\partial \bar{Z}}{\partial A} & \frac{\partial \bar{Z}}{\partial \phi\_\*} \end{vmatrix} = A\tag{27}$$

Thus,

$$p\_{A\phi\_\*}(A,\phi\_\mathbf{s}) = \frac{A}{2\pi\sigma^2} \exp\left(-\frac{A^2}{2\sigma^2}\right) = \left(\frac{1}{2\pi}\right) \left[\frac{A}{\sigma^2} \exp\left(-\frac{A^2}{2\sigma^2}\right)\right] = p\_{\phi\_\*}(\phi\_\mathbf{s})p\_A(A) \tag{28}$$

where,

$$p\_{\phi\_s}(\phi\_s) = \frac{1}{2\pi} \qquad \qquad 0 < \phi\_s \le 2\pi \tag{29}$$

is the uniform probability density underlying the random direction *φs*, and

$$p\_A(A) = \frac{A}{\sigma^2} \exp\left(-\frac{A^2}{2\sigma^2}\right) \tag{30}$$

is the Rayleigh distribution for the amplitude *A*.

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

where, assuming ergodicity, < · > denotes either an ensemble or time average, and where the separation line model above is used to calculate the force variance *σ*2. The side load components *Fsy* and *Fsz* are expressed with respect to rocket-fixed coordinates; see (Keanini et

In order to demonstrate the physical consistency of the model introduced in (Srivastava et al., 2010), (Keanini et al., 2011) first shows that the *assumed* properties, A) -C), can be *derived* from the simple *separation line model* developed in (Srivastava et al., 2010). Second, and as shown in the next subsection, the side load model in A) - C) then leads to *experimentally observed* side

Two important experimental and numerical observations concerning the side load, **Fs** (within

a) the probability density of the random amplitude, *A* = |**Fs**|, is a Rayleigh distribution (Deck

b) the random instantaneous direction, *φs*, of **Fs** is uniformly distributed over the periphery of the nozzle, or *pφ<sup>s</sup>* (*φs*) = 1/2*π*, where *pφ<sup>s</sup>* is the pdf of the side load direction (Deck &

Both observations can be *derived*, starting from the simple statistical model of random side loads, A) - C), immediately above. Thus, given *A* and *φs*, the instantaneous side load

*Fsy* = *A* cos *φ<sup>s</sup> Fsz* = *A* sin *φ<sup>s</sup>* Following (Srivastava et al., 2010), write *Fsy* and *Fsz* as *Fsy* = *Y*¯ = *A* cos *φ<sup>s</sup>* and *Fsz* = *Z*¯ = *A* sin *φs*; thus, the joint probability density associated with *Fsy* and *Fsz* can be expressed as

where *pA<sup>φ</sup><sup>s</sup>* (*A*, *φs*) is the joint pdf for the random amplitude and direction of **Fr**, and where

 

*A*

*<sup>σ</sup>*<sup>2</sup> exp

 <sup>−</sup> *<sup>A</sup>*<sup>2</sup> 2*σ*<sup>2</sup>

*∂Y*¯ *∂A ∂Y*¯ *∂φ<sup>s</sup> ∂Z*¯ *∂A ∂Z*¯ *∂φ<sup>s</sup>*

<sup>2</sup>*πσ*<sup>2</sup> exp

−*Y*¯ <sup>2</sup> <sup>+</sup> *<sup>Z</sup>*¯ <sup>2</sup> 2*σ*<sup>2</sup>

*pA<sup>φ</sup><sup>s</sup>* <sup>=</sup> <sup>|</sup>*J*|*pY*¯*Z*¯(*Y*¯, *<sup>Z</sup>*¯) (26)

= *A* (27)

0 < *φ<sup>s</sup>* ≤ 2*π* (29)

= *pφ<sup>s</sup>* (*φs*)*pA*(*A*) (28)

(25)

al., 2011; Srivastava et al., 2010).

load amplitude and direction densities (Keanini et al., 2011).

components in body-fixed *y* and *z* directions are given by

rigid, axisymmetric nozzles) are first noted:

& Nguyen, 2004; Deck et al., 2002), and

Nguyen, 2004; Deck et al., 2002).

the jacobian determinant is given by

<sup>2</sup>*πσ*<sup>2</sup> exp

 <sup>−</sup> *<sup>A</sup>*<sup>2</sup> 2*σ*<sup>2</sup> = 1 2*π*

*<sup>p</sup>φ<sup>s</sup>* (*φs*) = <sup>1</sup>

*pA<sup>φ</sup><sup>s</sup>* (*A*, *<sup>φ</sup>s*) = *<sup>A</sup>*

Thus,

where,

**3.1 Derivation of density functions for side load amplitude and direction**

*pY*¯*Z*¯(*Y*¯, *<sup>Z</sup>*¯) = *pY*¯(*Y*¯)*pZ*¯(*Z*¯) = <sup>1</sup>

Following (Srivastava et al., 2010), we restate *pY*¯*Z*¯ in terms of *A* and *φ<sup>s</sup>* as,


2*π*

#### **3.2 Ornstein-Uhlenbeck model of separation line dynamics**

Theoretical determination of rocket response to side loads requires that the time correlation function for either side load component, �*Fsα*(*t* � )*Fsα*(*t*)�, be first determined. As detailed in (Keanini et al., 2011), �*Fsα*(*t* � )*Fsα*(*t*)� is developed in two steps. First, and as detailed in this subsection, we propose(and physically justify) that local separation line dynamics can be modeled as an Ornstein-Uhlenbeck process. Once this assumption is made, then the second step rests on a rigorous argument showing that on the relatively long rocket dynamics time scale, the boundary layer separation line shape, and importantly, associated side load components, are all delta correlated in time. See (Keanini et al., 2011)for details.

The following simple, explicit stochastic model of separation line dynamics is proposed:

$$ds\_i(t) = -k s\_i(t) + \sqrt{D\_s} dW(t) \tag{31}$$

where *si* (*t*) = *s*(*φi*, *t*) is the instantaneous separation line position at *φi*, *k* and *Ds* are damping and effective diffusion coefficients, and *dW*(*t*) is a differential Weiner process. This equation, describing an Ornstein-Uhlenbeck process, allows straightforward, physically consistent calculation of statistical properties associated with separation line motion and, more importantly, serves as the first link in a chain that connects short-time-scale random separation line motion to short-time-scale random side loads, and in turn, to long-time-scale stochastic rotational rocket dynamics (Keanini et al., 2011).

The form of this equation is chosen based on the following experimental features, observed in shock-separated flows near compression corners and blunt fins:


Physically, the damping term captures the fact that the shock sits within a pressure-potential energy well. Thus, downstream shock excursions incrementally decrease and increase, respectively, upstream and downstream shock face pressures; the resulting pressure imbalance forces the shock back upstream. A similar mechanism operates during upstream excursions. Introduction of a Weiner process models the combined random forcing produced by advection of turbulent boundary layer structures through the upstream side of the shock foot and pressure oscillations emanating from the downstream separated boundary layer and recirculation zone.

We note that the proposed model is qualitatively consistent with Plotkin's model of boundary layer-driven shock motion near compression corners and blunt fins (Plotkin, 1975). Plotkin's

Rocket Ascent 17

<sup>171</sup> Shock-Induced Turbulent Boundary Layer Separation in Over-Expanded Rocket Nozzles: Physics, Models, Random Side Loads, and the Diffusive Character of Stochastic Rocket Ascent

Fig. 8. In-nozzle stochastic side loads versus rocket altitude. Adapted from (Srivastava et al.,

Fig. 9. Single realizations of rocket center of mass trajectory under random side loads and with side loads suppressed. *Xo*, and (*Yo*, *Zo*) denote, respectively, the vertical, and (mutually orthogonal) lateral displacements of the rocket center of mass, relative to the launch location.

Adapted from (Srivastava et al., 2010).

2010).

model, which captures low frequency spectra of wall pressure fluctuations within these flows, corresponds to a generalized Ornstein-Uhlenbeck process in which a deterministic linear damping term is superposed with a non-Markovian random forcing term. We use an ordinary OU process model, incorporating a Weiner process, since again, it is consistent with the above observations and more particularly, since it allows much simpler calculation of statistical properties (Keanini et al., 2011).

### **4. Rocket response to random side loads**

This section highlights recent results on numerical simulation of the stochastic ascent of a sounding-rocket-scale rocket, subjected to altitude dependent random nozzle side loads (Srivastava et al., 2010). The numerical simulations solve the full, coupled, nonlinear equations of rotational and translational rocket motion. The simulations include the effects of altitude-dependent aerodynamic drag forces, random nozzle side loads, associated random torques, mass flux damping torques (Keanini et al., 2011), and time-varying changes in rocket mass and longitudinal moment of inertia. In order to clearly isolate the effects of nozzle side loads on rocket translational and rotational dynamics, random wind loads are suppressed. A simple scaling argument (Keanini et al., 2011) indicates that random winds: i) under most conditions, do not excite rotational motion, and ii) simply function as an *additive* source of variance in the rocket's *translational* motion. In other words, wind appears to have minimal influence on the stochastic, altitude-dependent evolution of *rotational* dynamics. Rather, (launch-site-specific) mean and random winds simply produce whole-rocket, random, *lateral translational* motion, superposed on a deterministic translational drift.

Given side load direction and amplitude densities in Eqs. (29) and (30), respectively, altitude- dependent side loads are simulated using a Monte Carlo approach; the nonlinear, coupled equations of translational and rotational motion are then solved using fourth order Runge-Kutta integration (Srivastava et al., 2010). In estimating altitude-dependent means and variances in translational and rotational velocities and displacements, an ensemble of 100 simulated flights are used; it is found that estimated statistics do not vary significantly when using ensembles of 40 and 100 flights (Srivastava et al., 2010). The parameters employed in the simulations are characteristic of medium sized sounding rockets (Srivastava et al., 2010). Characteristic results are presented in figures 8 through 13. The stochastic evolution of both lateral side load components, *Fsy* and *Fsz*, observed during a single simulated realization, is shown in figure 8. Initially, i.e., at launch, the pressure jump across the separation-inducing shock is relatively high (Srivastava et al., 2010), and is manifested by somewhat higher initial side load amplitudes. However, as the rocket gains altitude, ambient pressure decays, and the cross-shock pressure jump decreases - characteristic side load magnitudes become smaller. At the instant when the slowly traveling in-nozzle shock reaches the nozzle exit, the nozzle flow becomes shock free, boundary layer separation ceases, and side loading stops. For these simulations, this instant corresponds to a flight time of 10.85 seconds (Srivastava et al., 2010), or an altitude of approximately 3.75 km. Note that side load amplitudes are significant, on the order of 10 to 16 % of the rocket's initial weight.

Figure 9 shows a single realization of the rocket's trajectory, under the action of random side loads, and is compared against the trajectory taken when side loads are suppressed. [Here, and throughout, *Xo* and (*Yo*, *Zo*) denote, respectively, the vertical, and (mutually orthogonal) lateral displacements of the rocket center of mass, relative to the launch location.] It is clear, that absent active control, a rocket can exhibit significant lateral displacements relative to the predicted zero-side-load path. Scaling shows that the characteristic magnitudes of random 16 Will-be-set-by-IN-TECH

model, which captures low frequency spectra of wall pressure fluctuations within these flows, corresponds to a generalized Ornstein-Uhlenbeck process in which a deterministic linear damping term is superposed with a non-Markovian random forcing term. We use an ordinary OU process model, incorporating a Weiner process, since again, it is consistent with the above observations and more particularly, since it allows much simpler calculation of statistical

This section highlights recent results on numerical simulation of the stochastic ascent of a sounding-rocket-scale rocket, subjected to altitude dependent random nozzle side loads (Srivastava et al., 2010). The numerical simulations solve the full, coupled, nonlinear equations of rotational and translational rocket motion. The simulations include the effects of altitude-dependent aerodynamic drag forces, random nozzle side loads, associated random torques, mass flux damping torques (Keanini et al., 2011), and time-varying changes in rocket mass and longitudinal moment of inertia. In order to clearly isolate the effects of nozzle side loads on rocket translational and rotational dynamics, random wind loads are suppressed. A simple scaling argument (Keanini et al., 2011) indicates that random winds: i) under most conditions, do not excite rotational motion, and ii) simply function as an *additive* source of variance in the rocket's *translational* motion. In other words, wind appears to have minimal influence on the stochastic, altitude-dependent evolution of *rotational* dynamics. Rather, (launch-site-specific) mean and random winds simply produce whole-rocket, random, *lateral*

Given side load direction and amplitude densities in Eqs. (29) and (30), respectively, altitude- dependent side loads are simulated using a Monte Carlo approach; the nonlinear, coupled equations of translational and rotational motion are then solved using fourth order Runge-Kutta integration (Srivastava et al., 2010). In estimating altitude-dependent means and variances in translational and rotational velocities and displacements, an ensemble of 100 simulated flights are used; it is found that estimated statistics do not vary significantly when using ensembles of 40 and 100 flights (Srivastava et al., 2010). The parameters employed in the simulations are characteristic of medium sized sounding rockets (Srivastava et al., 2010). Characteristic results are presented in figures 8 through 13. The stochastic evolution of both lateral side load components, *Fsy* and *Fsz*, observed during a single simulated realization, is shown in figure 8. Initially, i.e., at launch, the pressure jump across the separation-inducing shock is relatively high (Srivastava et al., 2010), and is manifested by somewhat higher initial side load amplitudes. However, as the rocket gains altitude, ambient pressure decays, and the cross-shock pressure jump decreases - characteristic side load magnitudes become smaller. At the instant when the slowly traveling in-nozzle shock reaches the nozzle exit, the nozzle flow becomes shock free, boundary layer separation ceases, and side loading stops. For these simulations, this instant corresponds to a flight time of 10.85 seconds (Srivastava et al., 2010), or an altitude of approximately 3.75 km. Note that side load amplitudes are significant, on the

Figure 9 shows a single realization of the rocket's trajectory, under the action of random side loads, and is compared against the trajectory taken when side loads are suppressed. [Here, and throughout, *Xo* and (*Yo*, *Zo*) denote, respectively, the vertical, and (mutually orthogonal) lateral displacements of the rocket center of mass, relative to the launch location.] It is clear, that absent active control, a rocket can exhibit significant lateral displacements relative to the predicted zero-side-load path. Scaling shows that the characteristic magnitudes of random

*translational* motion, superposed on a deterministic translational drift.

properties (Keanini et al., 2011).

**4. Rocket response to random side loads**

order of 10 to 16 % of the rocket's initial weight.

Fig. 8. In-nozzle stochastic side loads versus rocket altitude. Adapted from (Srivastava et al., 2010).

Fig. 9. Single realizations of rocket center of mass trajectory under random side loads and with side loads suppressed. *Xo*, and (*Yo*, *Zo*) denote, respectively, the vertical, and (mutually orthogonal) lateral displacements of the rocket center of mass, relative to the launch location. Adapted from (Srivastava et al., 2010).

Rocket Ascent 19

<sup>173</sup> Shock-Induced Turbulent Boundary Layer Separation in Over-Expanded Rocket Nozzles: Physics, Models, Random Side Loads, and the Diffusive Character of Stochastic Rocket Ascent

Fig. 11. Time-(altitude-)dependent means and variances of the rocket's lateral center of mass position (*Yo*, *Zo*) and lateral center of mass velocity (*vyo*, *vzo*). Adapted from (Srivastava et

can be modeled as an Ornstein-Uhlenbeck process, wherein side loads function, at least on the long rocket dynamics time scale, as Weiner processes. Simultaneous to stochastic pitch/yaw amplification, deterministic pitch/yaw damping, produced by incremental changes in the nozzle mass flux vector, i.e., mass flux damping (Keanini et al., 2011), counteract amplification. Thus, during the side load period, pitch/yaw velocities grow, but at an ever-decreasing rate. In the second regime, which begins when the separation-inducing shock exits the nozzle and side loads cease, mass-flux damping continues, driving pitch and yaw velocities toward zero (Keanini et al., 2011). These qualitative features are likewise apparent in plots of time-(altitude-)dependent pitch/yaw means and variances; see figure 13. As detailed in (Keanini et al., 2011), all of these features can be rigorously explained using, e.g., a course grained description of short-time-scale side load statistics, along with an asymptotic model of

Starting with the nonlinear equations of rotational motion, (Keanini et al., 2011) use a rigorous argument to show that the equations governing evolution of pitch and yaw rates assume the

al., 2010).

rocket translational and rotational motion.

**5. Stochastic rocket ascent as a diffusion process**

Fig. 10. Single realizations of rocket center of mass velocity under random side loads and with side loads suppressed. *vx*,*<sup>o</sup>* and (*vy*,*o*, *vz*,*o*), denote, respectively, the vertical, and (mutually orthogonal) lateral velocities of the rocket center of mass, relative to the launch location. Adapted from (Srivastava et al., 2010).

lateral displacements, which can be on the order of 1-4 kilometers over the simulated flight time of 25 seconds, are fully consistent with displacements estimated using characteristic side load magnitudes (Keanini et al., 2011; Srivastava et al., 2010). The several order of magnitude difference between (vertical) thrust forces and the small vertical component of **Fs** explains the result shown in figure 8a.

Random side loads produce significant excitation of pitch and yaw dynamics (Keanini et al., 2011; Srivastava et al., 2010). Thus, for example, once a rocket begins a random pitching and yawing motion, that motion continues, even after side loading stops (at *t* = 10.85 s). Random pitch and yaw, in turn, produce random lateral velocities. As indicated in figures 10 and 11, and as discussed in detail in (Keanini et al., 2011), a certain amount of time - an induction period - must pass before side-load-induced pitch and yaws grow large enough for significant lateral thrust components to appear. Once the induction period has passed, however, large lateral thrusts begin to act and the rocket begins to experience large lateral velocities and displacements. Figures 10 and 11 capture these essential dynamical features. The work in (Keanini et al., 2011) provides a detailed, physically-based analysis of the complex dynamics a rocket experiences during side loading.

As shown in (Keanini et al., 2011), and as indicated in figure 12, stochastic pitch and yaw rotational dynamics are characterized by two qualitatively distinct regimes. During the side load period, *t* ≤ 10.85 s, the action of random side loads on pitch and yaw angular velocities 18 Will-be-set-by-IN-TECH

Fig. 10. Single realizations of rocket center of mass velocity under random side loads and with side loads suppressed. *vx*,*<sup>o</sup>* and (*vy*,*o*, *vz*,*o*), denote, respectively, the vertical, and (mutually orthogonal) lateral velocities of the rocket center of mass, relative to the launch

lateral displacements, which can be on the order of 1-4 kilometers over the simulated flight time of 25 seconds, are fully consistent with displacements estimated using characteristic side load magnitudes (Keanini et al., 2011; Srivastava et al., 2010). The several order of magnitude difference between (vertical) thrust forces and the small vertical component of **Fs** explains the

Random side loads produce significant excitation of pitch and yaw dynamics (Keanini et al., 2011; Srivastava et al., 2010). Thus, for example, once a rocket begins a random pitching and yawing motion, that motion continues, even after side loading stops (at *t* = 10.85 s). Random pitch and yaw, in turn, produce random lateral velocities. As indicated in figures 10 and 11, and as discussed in detail in (Keanini et al., 2011), a certain amount of time - an induction period - must pass before side-load-induced pitch and yaws grow large enough for significant lateral thrust components to appear. Once the induction period has passed, however, large lateral thrusts begin to act and the rocket begins to experience large lateral velocities and displacements. Figures 10 and 11 capture these essential dynamical features. The work in (Keanini et al., 2011) provides a detailed, physically-based analysis of the complex dynamics

As shown in (Keanini et al., 2011), and as indicated in figure 12, stochastic pitch and yaw rotational dynamics are characterized by two qualitatively distinct regimes. During the side load period, *t* ≤ 10.85 s, the action of random side loads on pitch and yaw angular velocities

location. Adapted from (Srivastava et al., 2010).

a rocket experiences during side loading.

result shown in figure 8a.

Fig. 11. Time-(altitude-)dependent means and variances of the rocket's lateral center of mass position (*Yo*, *Zo*) and lateral center of mass velocity (*vyo*, *vzo*). Adapted from (Srivastava et al., 2010).

can be modeled as an Ornstein-Uhlenbeck process, wherein side loads function, at least on the long rocket dynamics time scale, as Weiner processes. Simultaneous to stochastic pitch/yaw amplification, deterministic pitch/yaw damping, produced by incremental changes in the nozzle mass flux vector, i.e., mass flux damping (Keanini et al., 2011), counteract amplification. Thus, during the side load period, pitch/yaw velocities grow, but at an ever-decreasing rate. In the second regime, which begins when the separation-inducing shock exits the nozzle and side loads cease, mass-flux damping continues, driving pitch and yaw velocities toward zero (Keanini et al., 2011). These qualitative features are likewise apparent in plots of time-(altitude-)dependent pitch/yaw means and variances; see figure 13. As detailed in (Keanini et al., 2011), all of these features can be rigorously explained using, e.g., a course grained description of short-time-scale side load statistics, along with an asymptotic model of rocket translational and rotational motion.

#### **5. Stochastic rocket ascent as a diffusion process**

Starting with the nonlinear equations of rotational motion, (Keanini et al., 2011) use a rigorous argument to show that the equations governing evolution of pitch and yaw rates assume the

Rocket Ascent 21

<sup>175</sup> Shock-Induced Turbulent Boundary Layer Separation in Over-Expanded Rocket Nozzles: Physics, Models, Random Side Loads, and the Diffusive Character of Stochastic Rocket Ascent

Fig. 13. Time-(altitude-)dependent means and variances of the rocket's pitch and yaw

*Int. Congress of Applied Mech.,* Stanford Univ., Palo Alto, CA.

hysteresis regime, *AIAA J.,* Vol. 42, No. 9, pp. 1878-1888.

configurations, *J. Turbulence,* Vol. 3, No. 1, pp. 57-57(1).

*Transfer and Fluid Mechanics Institute,* Stanford Univ. Press.

Bloomer, H. E., Antl, R. J. & Renas, P. E. (1961). Experimental study of the effects of geometric

Carriere, P., Sirieix, M. & Solignac, J. L. (1968). Properties de similitude des phenomenes de

Chapman, D. R., Kuehn, H. K. & Larson, H. K. (1958). Investigation of separated flows in

Deck, S. & Nguyen, A. T. (2004). Unsteady side loads in thrust-optimized contour nozzle at

Deck, S., Garnier, E. & Guillen, P. (2002). Turbulence modeling applied to space launcher

Dolling, D. S. & Brusniak, L. (1989). Separation shock motion in fin, cylinder, and compression ramp-induced turbulent interactions, *AIAA J.,* Vol 27, No. 6, pp. 734-742. Erdos, J. & Pallone, A. (1962). Shock-boundary layer interaction and flow separation, *Proc Heat*

Flemming, E. L., Chandra, S., Schoeberl, M. R. & Barnett, J. J. (1988). Monthly Mean Global

Climatology of Temperature, Wind, Geopotential Height, and Pressure for 0-120 km,

variables on performance of conical rocket engine exhaust nozzles, NASA TN D-846.

decollement laminaires ou turbulents en ecoulement supersonic nonuniforme, *12th*

supersonic and subsonic streams with emphasis on the effect of transition, NACA

angular velocities. Adapted from (Srivastava et al., 2010).

**6. References**

Report 1356.

NASA TM-100697.

Fig. 12. Ensemble of 100 realizations of yaw angular velocity evolution. Side loads cease at *t* = 10.85 s. Adapted from (Srivastava et al., 2010).

form of an Ornstein-Uhlenbeck process:

$$d\omega\_{\pm} = -A(t)\omega\_{\pm} \pm \sqrt{D(t)}dW\_{\pm} \tag{32}$$

where *ω*<sup>+</sup> = *ω<sup>y</sup>* and *ω*<sup>−</sup> = *ω<sup>z</sup>* correspond, respectively, to yaw and pitch rate, *A*(*t*) and *D*(*t*) are, respectively, the time-(altitude-)dependent effective damping and diffusion coefficients, and *W*±(*t*) are Weiner processes.

Fundamentally, and as detailed in (Keanini et al., 2011), equation (32) provides the key to analyzing both the rotational and translational rocket response to side loading. Physically, *A*(*t*) is roughly proportional to both the squared moment arm from the rocket center of mass to the nozzle exit, *L*<sup>2</sup> *ce*, as well as the mass flux magnitude, and is inversely proportional to the lateral moment of inertia. Likewise, *D*(*t*) is proportional *L*<sup>2</sup> *ce*, as well as the squared (altitude-dependent) pressure difference between the interior and exterior of the nozzle, and the squared altitude-dependent position of the separation-inducing shock.

Practically, the detailed formulas for *A*(*t*) and *D*(*t*) in (Keanini et al., 2011) are related to both rocket-specific design parameters, as well as universal, non-specific parameters characterizing in-nozzle, shock-boundary layer separation. Thus, as discussed in (Keanini et al., 2011), the formulas allow straightforward identification of design criteria for, e.g., enhancing pitch/yaw damping and/or suppressing diffusive, i.e., stochastic growth of random pitch and yaw.

Fig. 13. Time-(altitude-)dependent means and variances of the rocket's pitch and yaw angular velocities. Adapted from (Srivastava et al., 2010).
