## Abstract

High-fidelity numerical experiments and theoretical modelling are used to study the dynamics of a sounding-rocket-scale rocket, subject to altitude-dependent random wind and nozzle side loads and deterministic aerodynamic loading. This paper completes a series of studies that showed that Ornstein–Uhlenbeck (OU) rotational dynamics arise when random nozzle side loads dominate wind and aerodynamic loading. In contrast to the earlier work, this paper elucidates that under conditions where aerodynamic, wind and nozzle side loads are comparable, the rocket behaves as stochastic Brownian oscillator. The Brownian oscillator model allows straightforward interpretation of the complex rotational dynamics observed: three dynamical regimes—each characterized by differing balances between nozzle-side-load-induced torques, spring-like aerodynamic torques and mass flux damping torques—characterize rocket ascent. Further, the paper illuminates that in the limit where wind and aerodynamic loads are small, random mass flux variations *exponentially amplify* side-load-induced rotational stochasticity. In this practical limit, pitch/yaw dynamics are described by a *randomly damped* OU process; an exact solution of the associated Fokker–Planck equation can be obtained and used to compute, e.g. time-dependent pitch/yaw rate means and variances.

## 1. Introduction

Modern rockets, under the competing demands of cost and performance, are typically designed to achieve high thrust-to-weight ratios. This objective is often pursued using high area ratio, reduced length nozzles and/or optimized nozzle designs (Östlund & Muhammad-Klingmann 2005). While modern nozzles typically provide high in-vacuum performance, during low altitude flight, unsteady detachment of the over-expanded supersonic flow within the nozzle can lead to the generation of often dangerous nozzle side loads (Summerfield *et al.* 1954; Sunley & Ferriman 1964; Nave & Coffey 1973; Dumnov 1996; Alziary 2001; Brown *et al.* 2002; Östlund *et al*. 2004; Östlund & Muhammad-Klingmann 2005). The deleterious and sometimes catastrophic effects of these loads have been documented, for example, in J2, SSME Block I, Ariane-5, Chinese CZ-3B, Fastrac, Vulcain, and Japanese H-II and LE-7 engines (Nave & Coffey 1973; Chang 1996; Newman 2001; Wang 2004). In large engines, side loads can reach extreme magnitudes, for example, of the order of 250 000 pounds in Apollo Saturn V rockets (Brown *et al.* 2002).

Flow separation and associated side load generation results from complex, three-dimensional, unsteady, shock-turbulent boundary layer interactions within the nozzle. The fundamental problem of shock-boundary layer interaction has been well-studied (Zukoski 1967; Dolling & Murphy 1983; Zheltovodov 1996; Sinha *et al*. 2003). The more specialized problem of flow separation and side load generation in rocket nozzles has also attracted significant attention; Östlund & Muhammad-Klingmann (2005) and Östlund *et al*. (2004) provide comprehensive reviews. Original work on in-nozzle shock structure, side load generation and, for example, rocket structural response can be found, for example, in Onofri & Nasuti (1999), Shimizu *et al*. (2008), Frey & Hagemann (1999, 2000), Pekkari (1994), Schwane & Xia (2005) and Wang (2004).

While rocket dynamics comprises a vast, long-studied area of research (Meirovitch 1970; Brusch 1977; Ge & Cheng 1982; Calise *et al.* 1998; Banerjee 2000; Eke & Mao 2000; Dukeman 2002), prior to the recent work of Srivastava *et al.* (2010) and Keanini *et al.* (2011), the significant question concerning rocket response to random, in-nozzle side loads remained open.

The last two studies, which provide the foundation for the present investigation, lead to the following results:

— physical models of the stochastic, altitude-dependent evolution of the in-nozzle shock-boundary layer separation zone were developed;

— probability densities for the resulting side load amplitude and direction were theoretically derived and found in accordance with known densities;

— rigorous physical connections were established between stochastic separation line dynamics and associated side load generation—both taking place on a fast boundary layer time scale—and tied to slow, rocket-dynamics-time-scale rocket response;

— simplified linear rocket response models were derived and numerically validated; and

— analytical, numerically-validated expressions for altitude-dependent means and variances of sounding-rocket rotational and translational response were obtained.

Significantly, these studies focused on the practical, simpler limit in which random wind loads and (deterministic) aerodynamic lift, drag and side loads remained small relative to nozzle side loads.

The objective of the present paper thus centres on studying the effects of all three load features—random nozzle side loads, random winds and aerodynamic forces—on the ascent of sounding-rocket-scale rockets. As in Srivastava *et al.* (2010) and Keanini *et al.* (2011), our approach employs Monte Carlo ascent simulations, based on the full, coupled, nonlinear equations of motion, combined with development of simple models for interpreting observed numerical results.

The paper is organized as follows. Section 2 briefly outlines the procedures used in our numerical simulations. Section 3 highlights the two elements that are new to the simulation, the model used to represent the altitude-dependent random wind and the aerodynamic load model.

Section 4 derives a simplified linear stochastic model of the rocket's rotational response to the combined action of nozzle-side-load-induced torques, and wind and aerodynamic torques. Because the model has the same mathematical form, and captures physical features analogous to those underlying molecular-scale Brownian oscillators, we present a brief argument that ascending sounding rockets, under the action of aerodynamic, wind and nozzle side loads, can be viewed as macroscale analogues to Brownian oscillators. Practically speaking, the model provides a useful framework for interpreting observed numerical results.

Section 5 introduces a significant, though untreated question in rocket ascent dynamics: What are the dynamical effects of random variations in the mass flux vector—essentially the thrust—on rocket ascent? Such variations—owing, for example, to combustion instability (Sutton & Biblarz 2010)—are unrelated to random nozzle side loads, but can be treated using much of the same theoretical machinery. Thus, as in our treatment of nozzle side loads (Srivastava *et al.* 2010; Keanini *et al.* 2011), we suppress the complicating effects of aerodynamic and wind loading, and focus on the case where random mass flux variations act in concert with random nozzle side loads.

Finally, §6 presents and interprets the results of our numerical simulations, while §7 summarizes our principal findings.

## 2. Launch vehicle dynamics model

Numerical experiments are based on a model of rigid-body rocket motion, combined with a control volume formulation of the rocket's time-varying mass and internal flow (Srivastava *et al.* 2010). Essential kinematic and dynamic features are illustrated in figure 1. As shown, and following convention, three coordinate systems are defined: an Earth-fixed system, denoted *XYZ*, a rocket-fixed system, *xyz*, having origin at the rocket's instantaneous centre of mass, *O*, and a set of wind axes, *x*_{w}*y*_{w}*z*_{w}, having origin at the instantaneous centre of pressure, *Q*.

In addition to gravity, three external forces act on the rocket:

— A random, altitude-dependent nozzle side load acts normally to the nozzle interior (see inset on figure 1). As described in Srivastava

*et al.*(2010) and depicted in figure 2 in the electronic supplementary material, the axial location of the load is taken as the altitude-dependent, mean axial position of the in-nozzle boundary layer separation line,*x*_{s}(*t*). Using a simple Ornstein–Uhlenbeck model to describe evolution of the azimuthally varying random separation line component, Keanini*et al*. (2011) rigorously derived experimentally consistent probability densities for the side load amplitude and direction; these densities, which were assumed in Srivastava*et al.*(2010), then allow Monte Carlo simulation of the side load.— A random, altitude-dependent lateral wind load, directed in the Earth-fixed

*Y*- and*Z*-directions, is assumed to act at the rocket centre of pressure,*Q*. As described in §3, the wind load is produced by superposition of an altitude-varying, mean lateral wind and an altitude-independent, laterally homogeneous turbulent wind fluctuation. For simplicity, and for the sake of checking the consistency of numerical simulations, the*Y*- and*Z*-components of the mean profile are taken as identical.— Mach number- and rocket-orientation-dependent aerodynamic drag, lift and side forces,

*F*_{d},*F*_{1}and*F*_{c}, resolved respectively in the*x*_{w},*y*_{w}and*z*_{w}directions, also act at the rocket centre of pressure. These are also described in §3.

Additionally, the following assumptions are made:

— The rocket body is, at all times, axisymmetric.

— The internal flow of burnt products is axisymmetric and steady.

— At all times, the instantaneous mass centre lies on the rocket longitudinal axis,

*x*, i.e. the axis of symmetry. For the solid fuel sounding rocket model considered here, this is a reasonable assumption.— The exhaust gas flow is axisymmetric, uniform and steady.

Using the approach detailed in Srivastava *et al.* (2010), rocket centre-of-mass translational dynamics are governed by
2.1
2.2
and
2.3The terms (*p*_{e}−*p*_{a}) *A*_{e} and correspond, respectively, to the net pressure force at the nozzle exit and the inertial thrust force produced by expulsion of mass from the nozzle; here *p*_{e}, *p*_{a}, *A*_{e}, and *v*_{ex} are the nozzle exit and altitude-dependent atmospheric pressures, the nozzle exit area, the nozzle mass flow rate and the relative velocity of exiting nozzle flow. Nozzle side loads, resolved into the rocket-fixed *y*- and *z*-directions, are expressed as *F*_{sy} and *F*_{sz}. As discussed in Keanini *et al.* (2011), the terms involving correspond to reactive forces produced by pitch- and yaw-induced variations in the mass flux vector, these act perpendicularly to the instantaneous longitudinal rocket axis, *x*, *against* instantaneous pitch and yaw motion (figure 2 in the electronic supplementary material).

Rotation matrices between the three coordinate systems used in this problem give rise to the angles appearing in equations (2.1) through (2.3); the matrices taking the rocket-frame, *xyz*, to the Earth frame, *XYZ*, are given, for example, in Srivastava *et al.* (2010). Rotation matrices relating the rocket frame and wind frame, not considered in Srivastava *et al.* (2010), are given in the appendix A; as shown there, the matrices are stated in terms of the instantaneous angle of attack, *α*(*t*), and sideslip angle, *β*(*t*). Likewise, Euler angles for roll, pitch and yaw, denoted respectively by *φ*, *θ* and *ψ*, are related to the instantaneous rocket rotation vector, ** ω**, via another rotation matrix, which is also given in the appendix A.

### (a) Rotational dynamics

Following the method developed by Eke & Mao (2000) and Tran & Eke (2004), Srivastava *et al.* (2010) obtained the equations governing rocket rotational dynamics:
2.4
2.5
and
2.6

Moments of inertia with respect to rocket-fixed coordinates are denoted with *I*'s (*I*_{yy}=*I*_{zz}=*I*), and *L*, *b* and *R*_{e}, are the rocket length, half-length and nozzle exit radius, respectively.

Referring to figure 2 in the electronic supplementary material, the external moment, given by
2.7is determined by the torque produced by the instantaneous random nozzle side loads as well as the instantaneous aerodynamic load, *F*_{A}. Here, *r*_{OQ} is the moment arm from the centre of mass, *O*, to the centre of pressure, *Q*.

Although somewhat obscured at this point, equations (2.5) and (2.6), governing evolution of the yaw and pitch angles, *ψ*(*t*) and *θ*(*t*), respectively, each contain a torque damping term, a weak torque amplification term, ; an aerodynamic spring torque, *r*_{OQ}×*F*_{A}; and a stochastic, nozzle-side-load-induced forcing torque, corresponding to the first term on the right-hand side of (2.7). (Here, *α* corresponds to either *y* or *z*.) In addition, in order to isolate the dynamical effects of nozzle-side-load-induced torques, and wind and aerodynamic torques, we continue to focus on non-spinning rockets (Srivastava *et al.* 2010; Keanini *et al*. 2011). Because the initial roll rate, *ω*_{x}(0), is zero and because roll-inducing external torques remain small, the nonlinear spin-coupled terms containing *ω*_{x} in (2.5) and (2.6) thus likewise remain small.

As described in §4, these features, when combined with a simplified aerodynamic load model, lead to a physically transparent picture of stochastically forced pitch and yaw motion.

## 3. Aerodynamic and wind loads

This section briefly outlines the aerodynamic load and random wind models. Details of the random nozzle side load model can be found in Srivastava *et al*. (2010) and Keanini *et al*. (2011).

### (a) Aerodynamic forces

The aerodynamic lift, drag and side forces are given in generic form by
3.1where *S*=(1/2)*ρ*_{a}*A*_{R}|*v*_{rel}|^{2}, *A*_{R}=2*LR*_{r} is the rocket's projected area and *v*_{rel} is the rocket's instantaneous relative velocity
3.2and *C*_{η} represents the lift, drag or side load coefficient. Here, *v*_{wy} and *v*_{wz} are the absolute instantaneous wind velocity components, where the wind model is described in §3*b*.

Following Ferris (1967), Fuller & Foster (1963), Sutton & Biblarz (2010) and Shaver & Hull (1990), the drag and lift coefficients are modelled as:
3.3while the aerodynamic side force coefficient is modelled as (Fuller & Foster 1963; Stevens & Lewis 2003):
3.4The model in (3.3) captures the dominant dependence of *C*_{d,l} on angle of attack, *α*, and rocket Mach number, *M*_{R}, as well as the weak dependence on sideslip angle, *β*. Similarly, equation (3.4) reflects the weak dependence of *C*_{c} on *M*_{R} and *α* (Fuller & Foster 1963), and its strong dependence on *β*. (Note, *k*=−1/8 per degree of sideslip angle.) Detailed forms of the drag and lift coefficients are given, respectively, by
3.5aand
3.5bwhere the coefficient functions of *M*_{R} are determined by fits to representative experimental data given by Sutton & Biblarz (2010) and Shaver & Hull (1990).

We also assume that the centre of pressure, *Q* in figures 1 and, in the electronic supplementary material, figure 2, remains at a fixed position on the rocket's longitudinal axis. In reality, the centre of pressure varies, depending on the instantaneous angle of attack, rocket altitude, sideslip angle and Mach number. However, this assumption—which is supported by computational and experimental studies of high Mach number flow over finless, rocket-shaped projectiles (DeSpirito & Heavey 2006) and finned booster rockets (Sippel & Klevanski 2004)—can be argued using the following simple scaling analysis.

First, note that because pitch and yaw angles, *θ* and *ψ*, or equivalently, the angle of attack and sideslip angles, *α* and *β*, respectively, remain small during ascent (Busse & Kraft 1966; Hall *et al.* 1998), the COP, to *O*(*α*, *β*), remains on the rocket's longitudinal axis. Second, taking into account the altitude varying atmospheric density, and using a simple force balance, it is easily shown that the characteristic, dynamic pressure, , remains approximately fixed throughout the simulated ascent. Third, and as described later, because the drag coefficient likewise remains nominally fixed, then from the approximate relationship,
and the fact that the surface pressure distribution, *P*, is everywhere of the order of , we recognize that the centre of pressure thus remains nominally fixed. Hence, the moment arm from the centre of mass to the centre of pressure, |*r*_{OQ}|=*δ*_{c}, is fixed, and is here treated as a parameter.

Finally, an accurate model of altitude-dependent atmospheric pressure and density is required in order to capture mean separation line motion within the over-expanded nozzle, as well as altitude-dependent aerodynamic loads. Thus, a NASA atmospheric model is used in the present numerical simulations; see Srivastava *et al.* (2010) for details. The wind model is described in the electronic supplementary material.

## 4. Analogue self-propelled Brownian oscillator

### (a) Overview of §4 and 5

As will be shown, our numerical experiments expose complex rotational dynamics. In order to develop a simple physical framework for interpreting observed results, we use §4 to derive a simplified model of stochastic pitch and yaw evolution during ascent. It is found that the model has the same form as a micro- or molecular-scale, self-propelled Brownian oscillator, where self-propulsion is provided by random nozzle side loads.

Section 5, representing both a simplification *and* extension of the development in §4, focuses on a question that has not been treated: how do random variations in the nozzle mass flux vector, , affect rocket rotational dynamics? We consider the limit where nozzle-side-load-induced torques dominate wind- and aerodynamic torques (Srivastava *et al.* 2010; Keanini *et al.* 2011), focusing on the case where random mass flux variations and random nozzle side loads act in concert. Using the associated limiting form of the self-propelled Brownian oscillator model, we solve the corresponding Fokker–Planck equation to determine transition densities for altitude-dependent pitch and yaw rates, *ω*_{z} and *ω*_{y}. It is shown that under conditions where side load generation and mass flux variations take place on distinct time scales, a physically meaningful *average* transition density, averaged over the ensemble of randomly evolving mass flux vectors, , can be computed.

As an aside, we note that rocket body vibrations, particularly those involving the nozzle or its attachments, can play a significant role in nozzle side load generation (Brown *et al*. 2002). Because side loads, in turn, represent an important dynamical feature, the impact of rocket vibration on ascent constitutes an important open question.

### (b) Model development

In the absence of roll-inducing torques, i.e. torque components along the *x*-axis, and given an initial roll rate, *ω*_{x}(0)=0, (2.4) again implies that the roll rate remains zero, *ω*_{x}(*t*)=0. Assuming that pitch and yaw angles, *θ* and *ψ*, remain small, and likewise, that the angle of attack and sideslip angles, *α* and *β*, also remain small, then equations (2.5) and (2.6) simplify to the following forms:
4.1and
4.2where the damping coefficient, *A*(*t*), is given as
4.3aand where
4.3b(In detail, equations (4.1) and (4.2) follow from (2.5) and (2.6) by first noting that, for non-spin-stabilized rockets, *ω*_{x}(*t*)=0. Second, in equation (2.7), (where, for stability, *δ*_{c}<0, i.e. the centre of pressure is behind the centre of mass—figure 1) and . Finally, expressing in terms of , using the transformation matrices in equation (A1), and using the small angle approximation for *α* and *β*, yields equations (4.1) and (4.2)).

Considering first torque terms involving the side load components, *F*_{sz} and *F*_{sy}, because the latter act on short time scales relative to the long rocket dynamics time scale, and because each is a zero mean, Gaussian random process, delta correlated in time, these can be represented, on the long time scale, as Wiener processes (Keanini *et al.* 2011):
4.4where the diffusion coefficient is given by
4.5Referring to figure 2 in the electronic supplementary material, *R*(*t*) is the nozzle radius at the axial location of the mean boundary layer separation line, *x*_{s}(*t*), *H*(*t*) is the instantaneous rocket altitude, *P*_{i}(*H*(*t*)) is the in-nozzle pressure immediately upstream of the boundary layer separation line (at altitude *H*(*t*)), *P*_{a}(*H*(*t*)) is the corresponding atmospheric pressure, *σ*_{s} is the characteristic (axial) width of the in-nozzle shock-boundary layer interaction zone, *τ*_{c} is the correlation time for axial displacements of the stochastically evolving separation line and *ε* is the ratio of the short (boundary layer) time scale to the long (rocket dynamics) time scale. See Keanini *et al.* (2011) for a detailed discussion.

From (4.5), we observe that the degree of stochasticity introduced into the rocket's pitch and yaw motions, driven by nozzle-side-load-induced torques, increases, for example, as the square of the pressure difference, *P*_{i}(*H*(*t*))−*P*_{a}(*H*(*t*)), across the nozzle interior and exterior, and as the square of size of shock-boundary layer interaction zone, 2*πR*(*t*)*σ*_{s}. Likewise, rotational stochasticity is amplified as rotational inertia, *I*(*t*), decreases.

Turning to the torque produced by aerodynamic loads, the last terms on the right-hand sides of equations (4.1) and (4.2), assuming small angular deflections, and focusing on the drag coefficient given by equation (3.5a), we neglect, in (3.5a), linear and quadratic terms in the angle of attack, *α*. Likewise, the relative velocity, *v*_{rel}, in equation (3.2) is approximated as (where both approximations follow using straightforward order of magnitude estimates). Finally, using (A3*a*) and (A3*b*) (which relate angle of attack and sideslip angle, *α*(*t*) and *β*(*t*), to rocket and wind velocities, as well as pitch, yaw and roll angles), again introducing the small angle approximation, and introducing the assumption that , we obtain the following linear model for the stochastic evolution of pitch angle:
4.6A similar set of steps applied to equation (4.2) likewise leads to a linear evolution equation for yaw:
4.7

### (c) Analogue self-propelled Brownian oscillator

The simplified equations governing pitch and yaw evolution, equations (4.6) and (4.7), have the same form as those governing micro- or molecular-scale rotational or translational dynamics, driven by molecular-scale thermal forcing; see, e.g. Coffey *et al.* (1996). Here, the oscillator is obviously not immersed in a thermal bath; rather, stochastic forcing is produced by nozzle side loads. (Note, in deriving (4.6) and (4.7), owing to the dominance of the axial velocity, , over the random wind component, *u*′, random aerodynamic torques do not appear.)

Nevertheless, and with reference to D'Anna *et al*. (2003), we argue that the present system can be viewed as a macroscopic analogue of a *self-propelled Brownian oscillator*. D'Anna *et al*. (2003) created a macroscopic Brownian oscillator by immersing cone-shaped and cylindrically shaped macroscopic Brownian ‘particles’ in a bed of vibrationally fluidized glass beads; the vertically forced beads induced stochastic torsional motion on the Brownian particle. Here, taking the rocket nozzle as a structural element—within a larger Brownian particle—subjected to Brownian forcing, and recognizing that a random nozzle side load represents the summed, instantaneous aggregate of local, randomly distributed, net pressure loads acting around the nozzle periphery, it is clear that the present system shares many of the essential ingredients characterizing Brownian oscillators.

Three differences are apparent, however. First, stochastic forcing is produced by *virtual* particles, whose dynamical effect corresponds to the local, random pressure forces acting on the nozzle interior. Second, these virtual particles do not reside within an equilibrium thermal reservoir. Rather, and in analogy with D'Anna *et al*. (2003), where particle kinetic energy (driving Brownian motion) was supplied by mechanical (not thermal) means, the energy driving the present system is supplied chemically, via combustion upstream of the nozzle. Third, only a portion of the Brownian particle, i.e. the nozzle, experiences stochastic forcing.

In marked contrast to the analogue oscillator introduced by D'Anna *et al*. (2003), we note that in the present analogue, we have detailed understanding of (Srivastava *et al.* 2010; Keanini *et al.* 2011): (i) the physical origins and dynamical characteristics of the (virtual) Brownian particles producing Brownian oscillatory motion, and (ii) how the short time-scale dynamics of these virtual particles becomes rectified, on the long rocket dynamics time scale, into (slow time scale) rotational motion.

Practically speaking, equations (4.6) and (4.7) provide a useful basis for interpreting the rotational dynamics observed in our numerical experiments. They also provide the jumping off point for investigating stochastic rocket ascent driven by the combined action of random nozzle side loads and random mass flux variations, as detailed in §5.

As an aside, we note that a Green's function for the Brownian oscillator is available, applicable when the damping and spring coefficients—the coefficients in second and third terms on the left-hand sides of (4.6) and (4.7)—are fixed (Coffey *et al.* 1996). While these coefficients are time-dependent, over several-kilometre-scale segments of the rocket trajectory, they remain nominally fixed. Thus, Green's function can be used to construct segmented solutions that accommodate random initial conditions, and, for example, qualitative changes in a randomly varying mass flux,

## 5. Stochastic mass flux damping

Mass flux damping of pitch and yaw motion, again, represents an inertial effect in which rotation-induced changes in the nozzle mass flux vector, directed against the direction of rotation, always opposes rotation. The numerical model in §§2 and 3, and more transparently, the simplified model in §4 show that mass flux damping represents a true energy damping mechanism, analogous to a viscous dashpot attached to a driven spring-loaded mass. While this feature is well-known (Eke & Mao 2000; Tran & Eke 2004; Srivastava *et al.* 2010), the effects of stochastic variations in produced, for example, by random, asymmetric combustion variations upstream of the nozzle, have not been investigated. Mechanically, such variations are equivalent to attaching a randomly varying dash-pot to a stochastically forced spring-mass system.

Noting that stability requires small aerodynamic torques, we focus on the limit where wind and aerodynamic torques are negligible compared with those produced by nozzle side loads and mass flux damping. Using (4.6) and (4.7) as the starting point, with aerodynamic torque terms removed, we obtain a generic evolution equation for the pitch and yaw rate:
5.1where *ω* represents either *ω*_{z} or *ω*_{y}. The damping coefficient is expressed as a superposition of a mean, time-dependent component, *A*_{o}(*t*), and a random component, *A*′(*t*), where the former varies on the slow rocket rotational time scale, *τ*_{R}, and the latter varies on a much shorter time scale. From (4.3a), *A*_{o}(*t*) and *A*′(*t*) are in turn related to the mean and random mass fluxes:
5.2aand
5.2b

### (a) Pitch/yaw rate evolution under combined random nozzle side loading and mass flux fluctuations

Equation (5.1), which is of Ornstein–Uhlenbeck form, applies during the period when nozzle side loads act, 0≤*t*≤*t*_{s}. For times exceeding *t*_{s}, the last term on the right-hand side of (5.1) drops out and a separate treatment, described in §5*b*, is required.

The Fokker–Planck equation associated with (5.1) is given by
5.3where *p*(*ω*,*t*|*ω*_{o},*t*_{o}) is the transition density for the stochastic process, *ω*=*ω*(*t*), evolving according to (5.1), and *A*(*t*)=*A*_{o}(*t*)+*A*′(*t*).

Using the method described in Keanini (2011), equation (5.3) can be solved in closed form:
5.4where
5.5aand where
5.5bEquation (5.4) gives the transition density associated with a *single* realization of the random damping term, *A*(*t*); the result assumes that the time scale associated with random mass flux variations, , is significantly longer than that associated with side load generation. For mass flux vector fluctuations associated with, e.g. low-frequency combustion instabilities, i.e. those of the order of 10–50 Hz (Sutton & Biblarz 2010), this is a reasonable assumption because the characteristic frequency band for side load generation is of the order of 300–400 Hz (Keanini & Brown 2007).

In order to compute meaningful averages and variances for the pitch/yaw rate, we average, over *A*′, single-realization means and variances, 〈*ω*(*t*)〉 and *var* *ω*(*t*) (McComb 1990; Keanini 2011):
5.6and
5.7where expectations over the Weiner side loads are denoted by 〈⋅〉 and expectations over *A*′ are expressed as 〈⋅〉_{A}.

Under conditions where *A*′ is Gaussian and stationary, the following formulae, stated in terms of the correlation function, *R*(*τ*)=〈*A*′(*t*+*τ*)*A*′(*t*)〉, can be derived (Keanini 2011):
5.8and
5.9where and where *ω*_{o} is a deterministic, zero or non-zero initial pitch/yaw rate.

#### (i) Example: rotational response under short correlation time mass flux variations

For illustration, we focus on the case where the correlation time, for damping coefficient variations is short relative to the rocket rotational dynamics time scale, *τ*_{R}. In this limit, the time correlation for *A*′ can be approximated as
5.10The rotational dynamics time scale, *τ*_{R}, can be estimated either empirically using
where Δ*θ* and *ω*_{c} are the characteristic pitch/yaw angle and rate, as observed in our numerical simulations, or physically, by balancing the characteristic rotational inertia against side-load-induced torque:
where *I*_{o} is the initial moment of inertia. Both approaches lead to *τ*_{R}=*O* (0.2*s*).

Thus, the present limit holds when the characteristic frequency for mass flux variations, *f*_{M}, satisfies:
5.11where the characteristic side load frequency, *τ*^{−1}_{F}=*O*(300–400 *Hz*) (Keanini & Brown 2007), and *τ*^{−1}_{R}=*O*(1–10 *Hz*). In terms of mass flux variations driven by combustion instabilities, for example, this limit might arise during chugging and pogo instabilities (Sutton & Biblarz 2010).

In order to obtain easily interpreted results, we assume that *A*_{o}(*t*) remains constant over 0≤*t*≤*t*_{s}, where *t*_{s}≅10.8 s. This is a reasonable approximation because *A*_{o}(*t*), given by equation (5.2a), varies by less than 25 per cent. Using (5.10) in (5.8) and (5.9) and approximating *A*_{o}(*t*) as constant, one then obtains
5.12aand
5.12bDetails of these calculations (applied to randomly forced Burger's vortex sheets) are given in Keanini (2011).

Importantly, we first note that in the limit where random damping fluctuations are small, , or zero, and where the side load period is short, , equations (5.12a) and (5.12b) simplify to well-known mean and variance expressions for constant coefficient OU processes (Gardiner 2004); the latter are obtained by setting in (5.12a) and (5.12b). More generally, when is not small, (5.12*a*) and (5.12*b*) show that stochastic variations in the damping coefficient, *A*(*t*):

— inhibit the decay of the pitch/yaw rate to

*ω*_{o}, and—

*exponentially amplify*, via the term , pitch/yaw rate variance.

Physically, both effects are tied to the fact that tangentially acting differential mass flux vectors have large moment arms from the rocket's centre of mass (considering stable flight configuration/scenarios).

Equation (5.2*b*) also shows that rotational stochasticity is proportional to the side-load diffusion coefficient, *λ*_{f}. From (4.5), *λ*_{f}, in turn, depends on a number of parameters, including the nozzle radius, *R*(*t*), at the axial location of the mean boundary layer separation line, *x*_{s}(*t*), the instantaneous rocket altitude, *H*(*t*), the pressure difference between the nozzle interior and exterior (evaluated in the vicinity of *x*_{s}(*t*)), *P*_{i}(*H*(*t*))−*P*_{a}(*H*(*t*)), the characteristic (axial) width of the in-nozzle shock-boundary layer interaction zone, *σ*_{s}, the correlation time for axial displacements of the stochastically evolving separation line, *τ*_{c}, and the ratio of the short (boundary layer) time scale to the long (rocket dynamics) time scale, *ε*. Thus, as in the case of rotational stochasticity driven by nozzle side loads alone, the physically transparent form of equation (4.5) allows identification of practical strategies (Keanini *et al.* 2011) for minimizing rotational stochasticity when both random side loads and random mass flux variations act in tandem.

### (b) Pitch/yaw evolution following the nozzle side load period, *t*>*t*_{s}

Once side loading ceases, or more generally, when only vector mass flux variations take place, the Weiner forcing term in (5.1) turns off. Shifting the time origin to *t*=*t*_{s}, allowing for time-varying *A*_{o}(*t*), and integrating (5.1) then leads to:
5.13and
5.14where *τ*=*t*−*t*_{s}.

Equations (5.13) and (5.14) are general and only assume that *A*(*t*), or more specifically, the *y*- and *z*-components of , are independent, stationary, and Gaussian. We consider again the case where mass flux variations take place on time scales that are short relative to the rocket rotational time scale, . It is important to note that there is now no upper limit on the characteristic frequency associated with mass flux variations; the formulae in (5.15) and (5.16) below for the average pitch/yaw rate and associated variance hold for variations produced, for example, by high frequency and screeching instabilities (Sutton & Biblarz 2010).

Introducing (5.10) in (5.13) and (5.14) and defining then yields: 5.15and 5.16Thus, during the post-side-load period of ascent, or again, when only mass flux variations occur, we observe from (5.13) and (5.14) the following:

— small (or non-existent) mass flux variations, allow exponential damping of pitch/yaw rate means and variances towards zero;

— for , the mean pitch/yaw rate decays exponentially to zero while variance grows exponentially. This range of conditions, however, corresponding to extreme variations in , is not expected to arise except during catastrophic or near-catastrophic ascent; and

— for , both the mean pitch/yaw rate and variance grow exponentially; again, this corresponds to an extreme flight condition.

## 6. Results and discussion

Numerical experiments were performed by integrating equations (2.1) to (2.3) using the procedures outlined in Srivastava *et al.* (2010). Model parameters, given in table 1 in the electronic supplementary material, are representative of those associated with sounding rockets, e.g. the Peregrine (Dyer *et al.* 2008) and Black Brant (NASA Sounding Rocket Handbook 1999). The main features of the simulation are as follows:

— the altitude-dependent position of the mean separation line,

*x*_{s}(*H*(*t*)), is determined using boundary layer separation criteria developed in Keanini & Brown (2007);— given

*x*_{s}(*H*(*t*)), an instantaneous random separation line shape is generated relative to*x*_{s}(*H*(*t*)) (Srivastava*et al.*2010);— given the instantaneous separation line shape and the difference,

*P*_{i}(*H*(*t*))−*P*_{a}(*H*(*t*)) between the in-nozzle pressure immediately upstream of*x*_{s}(*H*(*t*)) and the ambient pressure, the instantaneous nozzle side load components,*F*_{sy}(*t*) and*F*_{sz}(*t*), are computed;— aerodynamic loads are computed as outlined in §3

*a*,*b*; and— kinematic relationships between instantaneous pitch, roll, and yaw angles

*θ*,*φ*and*ψ*, sideslip angle,*β*, angle of attack,*α*, roll, yaw and pitch rates,*ω*_{x},*ω*_{y}and*ω*_{z}, wind velocity components,*v*_{wy}and*v*_{wz}, and rocket centre of mass velocity, are given in the appendix A.

The first three steps comprise the instantaneous side load calculation and are performed only when the separation-inducing shock lies within the nozzle. Once the shock exits, here, at an approximate altitude of 3.17 km, side loads cease.

Aside: The following approximations, valid in the limit of small angular displacements, are useful in connecting Euler pitch and yaw angles, *θ* and *ψ*, pitch and yaw rates, *ω*_{z} and *ω*_{y}, and the angle of attack and sideslip angle, *α* and *β*:
6.1aand
6.1bIn addition, in our numerical experiments, estimated ensemble averages and variances are determined using a collection of 100 individual ascent realizations. It is found that time-dependent averages and variances computed over ensembles of 40 and 100 ascents vary by less than 6 per cent.

### (a) Three dynamic regimes

In interpreting numerical results, it is useful to refer to the simplified model in equations (4.6) and (4.7) and to break ascent into three regimes:

— Regime I is characterized by side-load-driven, under-damped, oscillatory flight. During this period, 0≤

*t*≤*t*_{s}, torques produced by nozzle side loads either dominate or are comparable to those produced by aerodynamic loads, (*L*−*b*−*x*_{s}(*t*))*F*_{sy,sz}(*t*)∼*F*_{A}(*t*)*δ*_{c}, while the mass flux damping torque is small. Here,*F*_{A}(*t*) represents either the aerodynamic drag, lift or side force. (Note,*t*_{s}=10.85 s in our simulations.)— Regime II corresponds to a period of aerodynamics-dominated, under-damped flight. It begins when nozzle side loads cease,

*t*=*t*_{s}, and lasts as long as atmospheric density remains high enough to allow aerodynamic spring forces to compete with mass flux damping, i.e. as long as*δ*_{c}*F*_{A}(*t*)>*I*(*t*)*A*(*t*)*ω*(*t*). As discussed later (and as indicated in figure 8, for example), Regime II ends roughly at*t*=20 s.— Once atmospheric density decays to small magnitudes, mass flux damping becomes dominant and oscillatory pitch and yaw rates are driven towards zero. This is Regime III, during which the right-hand side of the above inequality becomes much larger than the left.

The first two dynamical regimes are apparent in our numerical results, figures 2–7, while the third, again, can be ascertained in figure 8.

Figures 2 and 3 show single realizations of rocket pitch and yaw rate, and angle of attack and sideslip angle, respectively, under conditions where nozzle side loads act and do not act. Three features are apparent: (i) nozzle side loading amplifies both the rates and magnitudes of angular motion, (ii) owing to the relative weakness of the random wind component, i.e. , the mean wind profile (figure 1 in the electronic supplementary material) forces largely deterministic rotational motion once side loads cease (at *t*=*t*_{s}=10.85 s) and (iii) for *t*>*t*_{s}, aerodynamic torques, driven by the mean wind, quickly suppress rotational motion forced by random side loads. Consistent with equation (2.4) and as shown in figure 2*a*, owing to the absence of roll-inducing torques, the roll rate remains near its zero initial value.

Although the rocket is forced by both random nozzle side loads and random wind, the aerodynamic spring force plays a strongly organizing role in rotational dynamics. This feature is captured in figures 4–6, which show, respectively, the ensemble of 100 observed pitch and yaw rate histories, corresponding estimated means and standard deviations of the pitch/yaw rate, and estimated means and standard deviations of the angle of attack and sideslip angle. Oscillatory order, again driven by strongly deterministic lateral winds, dominant aerodynamic spring torque and under-damped response, almost certainly becomes less prominent for ascent into more turbulent winds.

Once nozzle side loads cease at *t*_{s}=10.85 s, the mean angle of attack and sideslip angle, figure 6, undergo small oscillations about zero; likewise, variances remain small and oscillatory. The initial, relatively large negative dip in *α* and *β* observed in figure 6 reflects, on the one hand, the rocket's zero initial velocity, and on the other, the existence of a small, non-zero lateral wind; refer to equations (A3*a*,*b*). Physically, the ground-level wind imposes a small initial tipping torque.

### (b) Qualitative insight using a simplified Brownian oscillator model

This section proposes a simplified version of the Brownian oscillator model, (4.6) and (4.7), as a basis for gaining qualitative insight into experimentally observed rotational dynamics (figures 2–6). Thus, consider a generic, simplified, averaged version of (4.6) and (4.7):
6.2where, 〈*η*(*t*)〉, via (6.1*a*,*b*), represents either the ensemble average angle of attack, 〈*α*(*t*)〉, or sideslip angle, 〈*β*(*t*)〉.

Here, averaging is over the Weiner nozzle side load, the damping coefficient, *A*(*t*), is given by (4.3*a*) and the *generic* spring coefficient, *K*(*t*) is given by
6.3and *C*_{o} is a characteristic magnitude for the lift, drag and aerodynamic side load coefficients. Thus, the simplified model in (6.2) and (6.3) circumvents speed- and orientation-dependent aerodynamic torques.

A representative fixed magnitude of *C*_{o}=0.15 is chosen for the aerodynamic coefficients, *C*_{l}, *C*_{d} and *C*_{c}. While these coefficients exhibit strong dependencies on *α* and *β* (Sutton & Biblarz 2010), in our experiments, owing to small sideslip angles and angles of attack, these dependencies are weak. Likewise, at small *α* and *β*, the dependence of *C*_{l}, *C*_{d} and *C*_{c} on the rocket Mach number, *M*_{R}, is, save for the transonic flow region, weak (Shaver & Hull 1990; Sutton & Biblarz 2010). For transonic *M*_{R}, i.e. *M*_{R} in the vicinity of 1, the drag coefficient, *C*_{d}, undergoes a delta-function-like variation; because transonic flow conditions exist only for approximately 1 s of any given simulated flight, however—refer to figure 11*a*—it is reasonable to neglect this effect.

Equation (6.2) is integrated using two further simplifications: (i) the rocket velocity component, , is approximated as ; figure 11*a*. (ii) The initial value of 〈*η*〉 is taken as −3^{°}, necessary in order to capture the initial wind-induced rocket tilt noted earlier. The other time-dependent parameters, *ρ*_{a}(*t*) and *I*(*t*), correspond to those used in the full numerical simulation.

Figure 7 compares estimated ensemble average angles of attack and sideslip angles, obtained in our numerical experiments, and shown in figure 6, versus those computed using the simplified model, equation (6.2). While there are clear differences, the relative simplicity of the Brownian oscillator model and the physical insight it provides argues its utility. Specifically:

— observed periodicity in averaged experimental results is qualitatively captured by the simplified model; in the latter case, the time-varying frequency is given by 6.4Note

*K*(*t*)≫*A*(*t*) from a trivial order analysis;— weak damping, which is more apparent in numerical experiments, is likewise captured; and

— the simplified model predicts sideslip and angle of attack amplitudes that are roughly of the observed magnitudes.

### (c) Further results: rotational dynamics

Pitch/yaw rate variance evolution observed here is compared in figure 8 against that predicted when aerodynamic torques are zero; as shown in Keanini *et al.* (2011), the latter well-predicts numerically observed zero-aerodynamic-torque results (Srivastava *et al*. 2010). The comparison shows the following:

— during the nozzle side load period, 0<

*t*≤*t*_{s}(Regime I), mean variance evolution under aerodynamic loading is qualitatively similar to that seen and predicted when aerodynamic loads are neglected (Keanini*et al.*2011). The Brownian oscillator model provides a simple interpretation: owing to the relative smallness of the random wind component, rotational stochasticity during the side load period is driven primarily by nozzle side loads. Strongly deterministic aerodynamic torque merely introduces an oscillatory, superposed component;— the variance increase observed between 15 and 20 s (Regime II), by contrast, is driven by the random wind component. Clearly, rotational dynamics remain under-damped, as confirmed by simple order of magnitude estimates of the mass flux damping and aerodynamic torques; and

— the variance decay seen for

*t*>∼20 s (Regime III) reflects ascent to altitudes where ambient density becomes so small that mass flux damping dominates the aerodynamic spring force.

While spring-like aerodynamic torques impose oscillatory order on the rocket's stochastic motion, the degree to which this occurs depends on the rocket's speed. This is apparent in both figure 9, which shows observed variance evolution for the angle of attack and sideslip angle, as well as from equations (6.2) and (6.3). As shown in figure 9, rotational variance, early in Regime I, grows explosively, driven by random nozzle side loads. However, once rocket speed, essentially given by in (6.3), becomes large enough, the strongly deterministic spring torque dominates, suppressing side-load-driven variance growth. (In the present model, aerodynamic suppression of side-load-induced stochasticity takes place over . Once side loads cease (at *t*_{s}), the rocket enters Regime II, where weak atmospheric turbulence sustains nominally fixed, near-zero levels of rotational variance.

### (d) Translational dynamics

Observed single realization and ensemble average translational dynamics are shown in figures 10–13. As expected, because the rocket is subjected to identical mean wind profiles in the lateral *y*_{o} and *z*_{o} directions, and because the random wind component is small relative to the local mean, lateral displacements, *y*_{o}(*t*) and *z*_{o}(*t*), and lateral velocities, *y*_{o} and *z*_{o} for single realizations (figures 10 and 11), and over the ensemble of 100 flights (figure 12), approximately mirror one another. Simple scaling arguments show that observed lateral displacements and velocities are consistent with the action of characteristic lateral wind loads extant during flight.

Inspection of figures 10–11, and comparison with results reported in Srivastava *et al.* (2010) shows that vertical velocities and displacements, and *x*_{o}(*t*), are largely unaffected by wind, aerodynamic and nozzle side loads. This reflects: (i) the large, approximate two order of magnitude difference between the vertical thrust force and the wind, aerodynamic and nozzle side loads, and (ii) the fact that dynamical angles remain small during ascent.

By contrast, nozzle side loads play a dominant role in lateral translational motion, *throughout* ascent. During the side load period, Regime I, the effect is *dynamic*, wherein nozzle side loads suppress rather than amplify lateral displacements and velocities; refer to figures 10 and 11. This counterintuitive result reflects three features: (i) Random, laterally acting nozzle side loads are uniformly distributed around the nozzle's interior periphery (Srivastava *et al.* 2010; Keanini *et al.* 2011). (ii) Nozzle side load magnitudes are approximately an order of magnitude larger than those of lateral wind and aerodynamic loads. (iii) Nozzle side loads act on short time scales relative to the rocket rotational time scale, *τ*_{R}. Together, these features produce a centring effect: weak, slow-time-scale, wind-induced lateral displacements are counteracted by large, fast-acting, uniformly distributed nozzle side loads.

During Regimes II and III, the effect of nozzle side loads on lateral translational dynamics is purely *kinematic*. Because lateral accelerations produced by wind and aerodynamic loads are small, of the order of and , respectively, the divergences observed over *t*>*t*_{s} between the ‘side-load-on’ and ‘side-load-off’ displacement and velocity histories, figures 10 and 11, reflect differing initial conditions extant at *t*=*t*_{s}, as determined by whether or not side loads act over 0<*t*≤*t*_{s}. Specifically, because post-side-load lateral accelerations are small, lateral velocities, for *t*>*t*_{s}, remain nominally fixed at the magnitudes observed at *t*=*t*_{s}:
6.5Importantly, equation (6.5) suggests that the time-varying ratio on the left-hand side remains approximately fixed throughout the post-side-load period, at the value observed at *t*=*t*_{s}, a prediction bore out in our single realization and ensemble average results.

Whereas random nozzle side loads play a dominant role in *translational motion*, throughout ascent, the role of side loads on *rotational motion* is largely limited to Regime I. Once side loading ceases and the rocket enters Regime II, aerodynamic torques, which dominate mass flux damping torques, quickly annihilate all memory of side-load-driven rotational dynamics.

The strongly organizing effect of aerodynamic torques, acting to pin the sideslip angle and angle of attack near zero—see figure 6—likewise suppresses stochasticity in translational motion. Figure 13 compares translational displacement and velocity variance growth, observed in the current set of numerical experiments, versus those observed in previous numerical experiments in which aerodynamic loads were neglected (Srivastava *et al.* 2010, Keanini *et al.* 2011). As shown, variances grow at much slower rates when aerodynamic torques are acting. The origin of post-side-load displacement and velocity variance growth in both experiments traces to a single dominant term (Keanini *et al*. 2011):
6.6representing variance produced by lateral components of rocket thrust, where, *F*_{T}(*t*) and *M*(*t*) are the time-dependent (deterministic) thrust force and rocket mass, and *ψ*_{α+}(*t*)=*ψ*(*t*) and *ψ*_{α−}(*t*)=*θ*(*t*) are the random yaw and pitch angles, respectively. Because, from (6.1*a*,*b*), *ψ*(*t*) and *θ*(*t*) are approximately equal to the angle of attack, *α*(*t*), and sideslip angle, *β*(*t*), then because the latter are, in the present experiments, aerodynamically forced towards zero, variance production by lateral thrust is likewise aerodynamically suppressed.

## 7. Summary and conclusions

The rotational and translational dynamics of sounding rocket-scale rockets have been studied under conditions where altitude-dependent random wind loads, speed, orientation and altitude-dependent random aerodynamic loads, and altitude-dependent random nozzle side loads act during rocket ascent. Numerical experiments and physical modelling reveal the following features:

— under (fairly typical) conditions where ascent takes place into a strongly deterministic wind, with nozzle side loads appearing in an over-expanded nozzle, and where the centre of pressure remains in the vicinity of the centre of mass, rocket rotational dynamics can be qualitatively described using a model analogous to that governing Brownian oscillators. Although the present macroscopic analogue does not meet a key criterion defining

*Brownian motors*,*viz.*, cyclic breaking of an underlying symmetry (Hänggi & Marchesoni 2009), the current system's rotational dynamics are clearly self-forced by the random nozzle side load;— numerical experiments, interpreted in light of the Brownian oscillator model, show that rotational dynamics can be characterized by three distinct regimes: (i) In Regime I, nozzle-side-load-driven torques force under-damped, aerodynamic oscillations. (ii) In Regime II, nozzle side loads cease while under-damped, wind-driven aerodynamic oscillations continue. (iii) In Regime III, atmospheric density decays to the point that aerodynamic spring torques are overcome by mass flux damping torques;

— under the typical circumstance where wind is strongly deterministic, stochastic nozzle-side-load-induced rotational dynamics are quickly annihilated (during Regime II) by wind-driven aerodynamic torques. Strongly deterministic, post-side-load, wind-driven aerodynamics thus serves as an organizing agent, suppressing variance growth in all rotational and translational dynamical variables;

— in contrast to the somewhat limited role played by nozzle side loads in determining rotational dynamics—limited largely to the side load period—side loads play a significant role, throughout ascent, in determining rocket translational motion. During the side load period, nozzle side loads

*suppress*the effect of lateral winds on lateral motion. Subsequent to the side load period—owing to typically weak lateral wind loads—translational motion is largely determined by the side-load-produced lateral velocity components extant at the end of Regime I; and— theoretical consideration of the effects of random mass flux variations on rotational dynamics, applicable in the practical limit where aerodynamic torques are small relative to side-load-driven torques, leads to the following results:

(i) An analytical solution to the associated Fokker–Planck equation can be obtained. The solution, giving the transition density for the pitch/yaw rate, applies to an OU process having a

*random*damping coefficient, driven by a Weiner forcing term; the random damping coefficient arises from random mass flux variations, while the Weiner term derives from the random nozzle side load torque.(ii) The analytical pitch/yaw rate transition density, in turn, allows analytical determination of, e.g. the time-varying ensemble average pitch/yaw rate and variance.

(iii) It is shown that during the side load period, random mass flux variations (which, physically, do not have the same origin as nozzle side loads)

*inhibit*mass flux damping of mean pitch/yaw rate and*exponentially amplify*pitch/yaw rate variance growth.(iv) Likewise, during the post side load period, stochastic mass flux variations again

*inhibit*damping of pitch/yaw rate means and variances. Variance amplification, however, can only occur under the extreme condition where the random mass flux variation exceeds the time-varying mean.

- Received May 6, 2012.
- Accepted August 16, 2012.

- This journal is © 2012 The Royal Society