# Turbulent Flow Characteristics in a Model of a Solid Rocket Motor Chamber with Sidewall Mass Injection and End-Wall Disturbance

Authors: F. A. Albatati; A. M. Hegab; M. A. Rady; A. A. Abuhabaya; S. M. El-Behery
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

The present analytical, numerical, and experimental investigations are performed to study the flow field in acoustically simulated solid rocket motor (SRM) chamber geometry. The computational solution is carried out for a high Reynolds number and low Mach number internal flows driven by sidewall mass addition in a long chamber with end-wall disturbances. This kind of flow (transient, weakly viscous, and contains vorticity) have several features in common with a turbulent flow field. The numerical study is performed by solving the unsteady Reynolds-averaged Navier–Stokes equations along with the energy equation using the control volume approach based on a staggered grid system. The v2-f turbulence model has been implemented in the current study. A comparison of the SIMPLE and PISO algorithms showed that both algorithms provide identical results, and the computational time using the PISO algorithm is higher by about 6% than the corresponding value of the SIMPLE algorithm. A fair agreement has been obtained between the numerical, analytical, and experimental results. Moreover, the results showed that the complex turbulent internal flow patterns are induced inside the chamber due to the strong interaction of the sidewall injection with the traveling acoustic waves. Such a complex internal structure is shown to be dependent on the piston frequency and sidewall mass flux. The current study, for the first time, emphasizes the acoustic-fluid dynamics interaction mechanism and the accompanying unsteady rotational fields along with the effect of the generated turbulence on the unsteady vorticity and its impact on the real burning rate.


---

## Body

## 1. Introduction

The flow through porous channels with a sidewall injection and end-wall disturbance is used to mimic the unsteady mass addition due to irregular propellant burning and fluid dynamics in a solid rocket motor (SRM) chamber [1]. The prediction of the complex acoustic, turbulent, and unsteady fluid flow interaction provides rocket designers with the essential characteristics required for the optimum rocket design and operation. Therefore, such a complex fluid dynamics problem has received significant attention as in [2]. In such a configuration, the internal flow exhibits time-dependent, compressible, laminar-turbulent transition, and turbulent characteristics over a large portion of the SRM chamber [3].Many investigations to chamber flow turbulence modelling, based onk-ε, k-ω, and full Reynolds stress methods, have been studied. Liou and Lien [4] revealed that the turbulence models, such as k-ε and k-ω, that used to predict the turbulence level, are not very successful in predicting the transverse location of the turbulence intensity peak as well as show the overpredicting turbulence level than the measured values. For example, Beddini [5] employed a full Reynolds stress turbulence model to analyze the flow in a porous channel. At a higher Reynolds numbers, transition to a turbulent velocity profile is predicted. The comparison to the experimental data in cold-flow simulation by Yamada et al. [6] indicates agreement with the laminar description for the mean velocity profile. The results by Sabnis et al. [7] using the low Reynolds number k-ε turbulence model show greatly overpredicting turbulence levels but indicates agreement with the experimental data for the axial mean velocity profiles. In general, in spite of the overpredicting turbulence level and the unsuccessful in predicting the transverse location of the turbulence peak, the turbulence models proved satisfactory in the pretransition region of the mean flow.In 2011, El-Askary et al. [8] simulated the internal fluid dynamics inside an SRM chamber with and without nozzle configuration. Authors assessed the capability of various turbulence models to predict the flow characteristics. The validity and accuracy of these models were evaluated through comparison with experimental data collected from literature. Among the investigated models which proved capability predicting the flow characteristics reasonably were the Reynolds stress model (RSM) and the v2-f and shear stress transport (SST) k-ω models. Moreover, Balabel et al. [9] indicated that only the v2-f and SST k-ω models can predict the supersonic flow in the nozzle configuration. Although the case of supersonic flow has been investigated using large eddy simulation technique (LES) by Apte and Yang [10] and Rhea et al. [11], the complexity associated with LES makes it less favorable compared to the Reynolds-averaged Navier–Stokes (RANS) equations which are superior especially when the differences in the predicted average quantities are marginal.Relevant studies in [2, 12] revealed that transient injection from the sidewall in an SRM chamber is the source of acoustic instability and disturbances that propagate with a low axial Mach number (M), high Reynolds number (Re), and mean flow. Moreover, in these studies, the asymptotic analysis which is used to simulate the corrugation that arises from the burning of composite propellants has demonstrated that an interaction between the acoustic transients and the fluid injected from the sidewall causes transverse axial velocity gradients (vorticity) and transverse temperature gradients (heat transfer) to appear with large amplitudes at the sidewalls.Alternatively, acoustics in the SRM chamber can be generated by a vibrating piston at the end-wall boundary. In such cases, the opposite boundary conditions in the SRM chamber play a significant role in the behaviour of the internal fluid dynamics and acoustic characteristics [1, 13]. Moreover, the end-wall disturbance consideration has been investigated both analytically and experimentally by Nasr et al. [14]. Nevertheless, the effect of the turbulence transition inside the SRM chamber was not considered. Javed and Chakraborty [15] numerically evaluated nozzle damping in a solid rocket motor using ANSYS CFX. A sinusoidal pressure pulse was enforced at the closed end of a channel-fitted with nozzle, for a given period. The pulse was then removed and the pressure was monitored at the head end. The pressure decay increased as the throat-to-port area decreased. The study also reported noticeable differences between the laminar and turbulent flow predictions.In the present study, the effect of the generated end-wall disturbance on the internal flow in an SRM chamber with a sidewall mass injection is investigated numerically considering the turbulence characteristics of the induced flow. The numerical model was especially developed for the purpose of this study in FORTRAN environment. The developed model was validated using a wide range of computational fluid dynamics (CFD) problems in different aerodynamics applications with reacting and nonreacting flows [8, 16]. The study also presents RANS equations for time-dependent, compressible turbulent flows along with the implemented turbulence model. Also, the numerical procedure and associated boundary conditions are explored.


## 2. Physical Model

The considered physical configuration of the nozzleless solid rocket motor consists of a two-dimensional rectangular chamber with a length (L) and half height (H) along with the porous sidewall. The aspect ratio (δ) is specified to be L/H. The chamber is fitted at the head end by a piston with a specified frequency of oscillation (ω) and transient axial velocity amplitude (ε), and it is open to the atmosphere at the downstream end, as shown in Figure 1.Figure 1
The physical model of the SRM nozzleless chamber.


## 3. Mathematical Formulation

The governing conservation equations of the flow can be written for 2D, time-dependent, compressible, and turbulent flow as follows.Continuity equation:(1)∂ρ∂t+∂∂xiρui=0.Momentum equation:(2)∂ρui∂t+∂∂xjρuiuj=−∂p∂xi+∂∂xjμ∂ui∂xj+∂uj∂xi−23δij∂ul∂xl+∂∂xj−ρui′uj′¯.Energy equation:(3)∂ρe∂t+∂ujρe+p∂xj=∂∂xjλ+CpμtPrt∂T∂xj−uiρui′uj′¯,where ui denotes the mean velocities and ui′ and uj′ are the turbulent fluctuations. The total energy per unit volume is defined as ρe=ρCvT+0.5ρui2, where Cv is the specific heat at a constant volume and T is the temperature.The momentum equations contain additional terms, known as Reynolds stresses, which represent the effects of turbulence. These Reynolds stresses,−ρui′uj′¯, must be modelled to close the governing system of equations. This can be achieved through the appropriate selection of a turbulence model.For compressible flow, equations (1)–(3) can be interpreted as Favre-averaged conservation equations with the velocities representing mass-averaged values. The Reynolds-averaged conservation equations for time-dependent compressible turbulent flow along with a turbulence model are coupled with the equation of state p=ρRT to close the system of governing equations.


### 3.1. Turbulence Modelling

In industrial CFD applications, RANS modelling remains one of the main approaches when dealing with turbulent flows. The Reynolds-averaged approach to turbulence modelling requires appropriate modelling of the Reynolds stresses in equation (3). The conventional method to relate the Reynolds stresses to the mean velocity gradients is based on the Boussinesq hypothesis [17]:(4)−ρui′uj′=μt∂ui∂xj+∂uj∂xi−23ρk+μt∂uk∂xkδij,where δij is the Kronecker delta function (δij = 1 if i = j, and δij = 0 if i ≠ j) and k is the turbulent kinetic energy. One of the most challenging problems of turbulent flow is the estimation of the turbulent viscosity. In order to obtain the turbulent viscosity, other transport equations are needed. These equations differ from one model to another as mentioned in detail in our previous work [9]. The general transport equations for the adopted models are given below, while the different terms and coefficient of the turbulence models adopted are given in Table 1.Table 1
Coefficient for the adopted turbulence models.
STDk-ɛ modelETDk-ɛ modelv2˗f modelSSTk-ωβ11110β20000.09β30000.0828σk10.7511σɛ1.31.151.30σω0001.168Pk2vtSijSij2vtSijSij2vtSijSij2vtSijSijSij0.5 (ui,xj + uj,xi)0.5 (ui,xj + uj,xi)0.5 (ui,xj + uj,xi)0.5 (ui,xj + uj,xi)Pk∗PkPkPkmin (Pk,10, 10ɛ)C1ɛ1.41.151.4(1 + 0.05(k/v2)0.5)0C2ɛ1.921.91.90C3ɛ00.2500Tk/ɛk/ɛmax (k/ɛ, 6 (μ/ρɛ)0.5)1/ωμtρCμk2/ɛρCμk2/ɛρCμv2T0.31ρk/max(0.31ω, F(2SijSij)0.5)Cμ0.090.090.220L000.23 max (k1.5/ɛ, 70 (v0.75ɛ0.25))0F2000tanh (D1)2D1000((2k0.5/0.09ωy), 500v/ωy2)γ10000.4403FSST0002.336 (1 −F1) (1/ω) (k,xjω,xj)F1000tanh (D2)4D2000min(max(D1, 4.672ρk/D3y2))D3000max (2.336ρ(1/ω) (k,xjω,xj), 1·e−10)Thek-equation:(5)∂∂xjρujk=∂∂xjμ+μtσk∂k∂xj+ρPk∗−β1ϵ−β2kω.Theɛ-equation:(6)∂∂xjρujε=∂∂xjμ+μtσε∂ε∂xj+ρTC1εpk−C2εε+C3εPk2ε.Thev2-equation:(7)∂∂xjρujv2=∂∂xjμ+μtσk∂v2∂xj+ρkf−v2kε.Thef-equation:(8)L2∂2f∂xj2−f=1T−4.6v2k−0.83−0.3Pkk.Theɷ-equation:(9)∂∂xjρujω=∂∂xjμ+μtσω∂ω∂xj+ρ2ϑ1μt−ρβ3ω2+FSST.For the calculation of the turbulent viscosity, thev2-f turbulence model by Lien and Kalitzin [18] is applied. The distinguishing feature of the v2-f model is its use of the velocity scale, v2, instead of the turbulent kinetic energy, k, for evaluating the eddy viscosity. v2, which is the velocity fluctuation normal to the streamlines, provides the right scaling in representing the damping of turbulent transport close to the wall, a feature that k does not provide.This model requires the solution to three transport equations, where the transport equation of the turbulent kinetic energy and its dissipation rate can be written as follows [18]:(10)∂ρk∂t+∂∂xjρujk=∂∂xjμ+μtσk∂k∂xj+ρPr−ε,∂ρε∂t+∂∂xjρujε=∂∂xjμ+μtσε∂ε∂xj+ρC1εPrεk−C2εε2k+C3εPr2k,uj∂v¯′2∂xj=kf−6v¯′2εk+∂∂xjν+νt∂v¯′2∂xj,L2∂2f∂xj2−f=1TC1−6v¯′2k−23C1−1−C2Gk,T=maxkε,6νε,L=CLmaxk3/2ε,Cην3ε1/4,where Pr and ε represent the production rate and dissipation rate of the turbulent kinetic energy, k, respectively. The production rate is related to the mean strain of the velocity field through the Boussinesq assumption:(11)Pr=μtS2,where S is the strain rate magnitude and it can be defined as follows:(12)S=1ρ∂uj∂xi+∂ui∂xj−23δij∂uk∂xk∂ui∂xj.The turbulent viscosity is given by the following:(13)μt=ρCμv2T.The constants of the model are given in Table1 and can be extracted and written for the v2-f turbulence model as follows:(14)Cμ=0.22,σk=1,σε=1.3,C1ε=1.41+0.05kv2,C2ε=1.9,C1=1.4,CL=0.23,Cη=70.


## 3.1. Turbulence Modelling

In industrial CFD applications, RANS modelling remains one of the main approaches when dealing with turbulent flows. The Reynolds-averaged approach to turbulence modelling requires appropriate modelling of the Reynolds stresses in equation (3). The conventional method to relate the Reynolds stresses to the mean velocity gradients is based on the Boussinesq hypothesis [17]:(4)−ρui′uj′=μt∂ui∂xj+∂uj∂xi−23ρk+μt∂uk∂xkδij,where δij is the Kronecker delta function (δij = 1 if i = j, and δij = 0 if i ≠ j) and k is the turbulent kinetic energy. One of the most challenging problems of turbulent flow is the estimation of the turbulent viscosity. In order to obtain the turbulent viscosity, other transport equations are needed. These equations differ from one model to another as mentioned in detail in our previous work [9]. The general transport equations for the adopted models are given below, while the different terms and coefficient of the turbulence models adopted are given in Table 1.Table 1
Coefficient for the adopted turbulence models.
STDk-ɛ modelETDk-ɛ modelv2˗f modelSSTk-ωβ11110β20000.09β30000.0828σk10.7511σɛ1.31.151.30σω0001.168Pk2vtSijSij2vtSijSij2vtSijSij2vtSijSijSij0.5 (ui,xj + uj,xi)0.5 (ui,xj + uj,xi)0.5 (ui,xj + uj,xi)0.5 (ui,xj + uj,xi)Pk∗PkPkPkmin (Pk,10, 10ɛ)C1ɛ1.41.151.4(1 + 0.05(k/v2)0.5)0C2ɛ1.921.91.90C3ɛ00.2500Tk/ɛk/ɛmax (k/ɛ, 6 (μ/ρɛ)0.5)1/ωμtρCμk2/ɛρCμk2/ɛρCμv2T0.31ρk/max(0.31ω, F(2SijSij)0.5)Cμ0.090.090.220L000.23 max (k1.5/ɛ, 70 (v0.75ɛ0.25))0F2000tanh (D1)2D1000((2k0.5/0.09ωy), 500v/ωy2)γ10000.4403FSST0002.336 (1 −F1) (1/ω) (k,xjω,xj)F1000tanh (D2)4D2000min(max(D1, 4.672ρk/D3y2))D3000max (2.336ρ(1/ω) (k,xjω,xj), 1·e−10)Thek-equation:(5)∂∂xjρujk=∂∂xjμ+μtσk∂k∂xj+ρPk∗−β1ϵ−β2kω.Theɛ-equation:(6)∂∂xjρujε=∂∂xjμ+μtσε∂ε∂xj+ρTC1εpk−C2εε+C3εPk2ε.Thev2-equation:(7)∂∂xjρujv2=∂∂xjμ+μtσk∂v2∂xj+ρkf−v2kε.Thef-equation:(8)L2∂2f∂xj2−f=1T−4.6v2k−0.83−0.3Pkk.Theɷ-equation:(9)∂∂xjρujω=∂∂xjμ+μtσω∂ω∂xj+ρ2ϑ1μt−ρβ3ω2+FSST.For the calculation of the turbulent viscosity, thev2-f turbulence model by Lien and Kalitzin [18] is applied. The distinguishing feature of the v2-f model is its use of the velocity scale, v2, instead of the turbulent kinetic energy, k, for evaluating the eddy viscosity. v2, which is the velocity fluctuation normal to the streamlines, provides the right scaling in representing the damping of turbulent transport close to the wall, a feature that k does not provide.This model requires the solution to three transport equations, where the transport equation of the turbulent kinetic energy and its dissipation rate can be written as follows [18]:(10)∂ρk∂t+∂∂xjρujk=∂∂xjμ+μtσk∂k∂xj+ρPr−ε,∂ρε∂t+∂∂xjρujε=∂∂xjμ+μtσε∂ε∂xj+ρC1εPrεk−C2εε2k+C3εPr2k,uj∂v¯′2∂xj=kf−6v¯′2εk+∂∂xjν+νt∂v¯′2∂xj,L2∂2f∂xj2−f=1TC1−6v¯′2k−23C1−1−C2Gk,T=maxkε,6νε,L=CLmaxk3/2ε,Cην3ε1/4,where Pr and ε represent the production rate and dissipation rate of the turbulent kinetic energy, k, respectively. The production rate is related to the mean strain of the velocity field through the Boussinesq assumption:(11)Pr=μtS2,where S is the strain rate magnitude and it can be defined as follows:(12)S=1ρ∂uj∂xi+∂ui∂xj−23δij∂uk∂xk∂ui∂xj.The turbulent viscosity is given by the following:(13)μt=ρCμv2T.The constants of the model are given in Table1 and can be extracted and written for the v2-f turbulence model as follows:(14)Cμ=0.22,σk=1,σε=1.3,C1ε=1.41+0.05kv2,C2ε=1.9,C1=1.4,CL=0.23,Cη=70.


## 4. Numerical Procedure

The adopted model, including the conservation equations and turbulence model equations, is solved using the numerical method proposed by Versteeg and Malalasekera [19] based on the finite volume approach. The governing equations are discretized using a hybrid scheme for all variables except the density, which is interpolated using a first-order upwind scheme to obtain a linear system of algebraic equations. This system of equations is solved using the TDMA technique described in [19]. The second-order implicit discretization is applied to the unsteady term. An essential step of the numerical procedure is the linearization of the source terms that vary according to the considered equations.The solution is carried out using either SIMPLE or PISO algorithms, which are extended to compressible flow according to [19], and the program is considered to reach the final solution when the maximum normalized residual for all variables drops to 10−4 for a steady-state solution and to 10−8 for an unsteady solution during each time step. The integration time step Δt is assumed to be constant and small enough to achieve a stable and nonoscillatory solution, and it is selected based on the independent study. Different time steps are used in the simulation, and the time histories for pressure and velocity are compared. A time step of 10−6 s was found sufficient to produce independent results for all cases.


### 4.1. Computational Domain and Boundary Conditions

The computational domain is a two-dimensional rectangular channel with the same dimensions as those in [20] as in Figure 1. The channel height H is 10.3 mm, and its length L is 58.1 mm. A 250 × 250 computational grid is used based on a grid-independent study by [20]. Three kinds of boundary conditions are used in this study.The first boundary condition case is concerned with studying the internal flow field when only acoustic waves are generated at the closed end, while the other free end is open to the atmosphere. In this case, the no-slip condition is applied to the velocity components at the sidewalls (i.e.,u = v = 0). At the head end, where x = 0, the stream wise velocity is considered in the form u = ε sin (ωt), where ε and ω represent the wave amplitude and forced frequency, respectively. At the duct exit, a known pressure (atmospheric pressure) is applied, and the other parameters are extrapolated using a zero gradient.The second boundary condition case is the same as the experimental boundary condition of the experimental work conducted by Avalon et al. [21]. In this case, the air is injected uniformly across the bottom sidewall of the channel, where y = 0, and the head end is closed, while the other free end, x = L, is open to the atmosphere. At the head end, where x = 0, and the top wall, where y = H, the no-slip boundary conditions are applied, u = v = 0. At the duct exit, the gauge pressure is taken as 137 400 Pa, and the other parameters are extrapolated using the zero gradients at the exit ∂ϕ/∂x=0. At the bottom sidewall, the mass flux value of 2.619 kg/m2 s is used which is relevant to the real burning mass flux of the SRM propellant (Vinj ∗ ρg = ρs ∗ rb), where ρs, ρg, and rb represent the density of combustion product in the gas phase, the density of the solid propellant in the solid phase, and the real burning rate, respectively.In the third boundary condition case, the air is injected uniformly across the sidewalls of the chamber, and acoustic waves are generated at the head end, while the other free end is open to the atmosphere. The same boundary conditions for the first case are applied except that, at both sidewalls, the crosswise velocity is taken to be the injection velocityVinj.In the following sections, the samples of the obtained results are illustrated. The critical elements of the performed cases are the injection velocityVinj and the axial perturbation velocity of the vibrated piston up. Three different cases are considered, as summarized in Table 2Table 2
Summary of boundary conditions for the test cases.
VinjupCase I0.0u=3.4sinωtCase IIVinj=2.619/ρm/sec0.0Case IIIVinj=0.3m/secu=3.4sinωt


### 4.2. Analytical Solution

An analytical solution for the first case is developed for the reduced form of the governing system of equations. The analytical results provide important insight into the preliminary results of the problem under consideration, and they might be used as an indication of the logical level of both the numerical and experimental results. The analytical approach is based on the reduced form of the Navier–Stokes equations using asymptotic techniques [22–25], and the final solution for pressure and velocity can be written as follows.The nonresonance acoustic solution for axial velocity is as follows:(15)u¯x¯,t¯=ε¯sinω¯t¯+∑n=0∞2ε¯ω¯bn2−ω¯2ω¯bn2sinω¯t¯−sinbnt¯sinbnx¯,forω¯≠bn∗,where bn is the wavenumber, bn = (n + 1/2) π, and ω¯ is the dimensionless forced frequency. The perturbation pressure p¯x¯,t¯ is found using the following equation [22, 23]:(16)∂u¯x¯,t¯∂t¯=−1γM∂p¯x¯,t¯∂x¯.The boundary condition at the duct head end,x¯=1, reads p¯=1. Integrating equation (16) yields the pressure perturbation in the following form:(17)p¯x¯,t¯=1.0+γMε¯ω¯1.0−x¯cosω¯t¯+γM∑n=0∞2ε¯ω¯bn2−ω¯2ω¯2bn2sinω¯t¯−bnsinbnt¯cosbnx¯bn.The overbar in equations (15)–(17) indicates the dimensionless quantity. The following parameters are used for nondimensionalization:(18)u¯=uUzo,ε¯=εUzo,ω¯=ωLCzo,t¯=tCzoL,x¯=xL,p¯=ppo,where Czo = (γRTo)0.5 is the characteristic sound velocity, Uzo = MCzo is the characteristic axial velocity, and po is the pressure at the duct outlet.


## 4.1. Computational Domain and Boundary Conditions

The computational domain is a two-dimensional rectangular channel with the same dimensions as those in [20] as in Figure 1. The channel height H is 10.3 mm, and its length L is 58.1 mm. A 250 × 250 computational grid is used based on a grid-independent study by [20]. Three kinds of boundary conditions are used in this study.The first boundary condition case is concerned with studying the internal flow field when only acoustic waves are generated at the closed end, while the other free end is open to the atmosphere. In this case, the no-slip condition is applied to the velocity components at the sidewalls (i.e.,u = v = 0). At the head end, where x = 0, the stream wise velocity is considered in the form u = ε sin (ωt), where ε and ω represent the wave amplitude and forced frequency, respectively. At the duct exit, a known pressure (atmospheric pressure) is applied, and the other parameters are extrapolated using a zero gradient.The second boundary condition case is the same as the experimental boundary condition of the experimental work conducted by Avalon et al. [21]. In this case, the air is injected uniformly across the bottom sidewall of the channel, where y = 0, and the head end is closed, while the other free end, x = L, is open to the atmosphere. At the head end, where x = 0, and the top wall, where y = H, the no-slip boundary conditions are applied, u = v = 0. At the duct exit, the gauge pressure is taken as 137 400 Pa, and the other parameters are extrapolated using the zero gradients at the exit ∂ϕ/∂x=0. At the bottom sidewall, the mass flux value of 2.619 kg/m2 s is used which is relevant to the real burning mass flux of the SRM propellant (Vinj ∗ ρg = ρs ∗ rb), where ρs, ρg, and rb represent the density of combustion product in the gas phase, the density of the solid propellant in the solid phase, and the real burning rate, respectively.In the third boundary condition case, the air is injected uniformly across the sidewalls of the chamber, and acoustic waves are generated at the head end, while the other free end is open to the atmosphere. The same boundary conditions for the first case are applied except that, at both sidewalls, the crosswise velocity is taken to be the injection velocityVinj.In the following sections, the samples of the obtained results are illustrated. The critical elements of the performed cases are the injection velocityVinj and the axial perturbation velocity of the vibrated piston up. Three different cases are considered, as summarized in Table 2Table 2
Summary of boundary conditions for the test cases.
VinjupCase I0.0u=3.4sinωtCase IIVinj=2.619/ρm/sec0.0Case IIIVinj=0.3m/secu=3.4sinωt


## 4.2. Analytical Solution

An analytical solution for the first case is developed for the reduced form of the governing system of equations. The analytical results provide important insight into the preliminary results of the problem under consideration, and they might be used as an indication of the logical level of both the numerical and experimental results. The analytical approach is based on the reduced form of the Navier–Stokes equations using asymptotic techniques [22–25], and the final solution for pressure and velocity can be written as follows.The nonresonance acoustic solution for axial velocity is as follows:(15)u¯x¯,t¯=ε¯sinω¯t¯+∑n=0∞2ε¯ω¯bn2−ω¯2ω¯bn2sinω¯t¯−sinbnt¯sinbnx¯,forω¯≠bn∗,where bn is the wavenumber, bn = (n + 1/2) π, and ω¯ is the dimensionless forced frequency. The perturbation pressure p¯x¯,t¯ is found using the following equation [22, 23]:(16)∂u¯x¯,t¯∂t¯=−1γM∂p¯x¯,t¯∂x¯.The boundary condition at the duct head end,x¯=1, reads p¯=1. Integrating equation (16) yields the pressure perturbation in the following form:(17)p¯x¯,t¯=1.0+γMε¯ω¯1.0−x¯cosω¯t¯+γM∑n=0∞2ε¯ω¯bn2−ω¯2ω¯2bn2sinω¯t¯−bnsinbnt¯cosbnx¯bn.The overbar in equations (15)–(17) indicates the dimensionless quantity. The following parameters are used for nondimensionalization:(18)u¯=uUzo,ε¯=εUzo,ω¯=ωLCzo,t¯=tCzoL,x¯=xL,p¯=ppo,where Czo = (γRTo)0.5 is the characteristic sound velocity, Uzo = MCzo is the characteristic axial velocity, and po is the pressure at the duct outlet.


## 5. Code Validations

The analytical and experimental results by Hegab and Nasr [22] and Hegab et al. [23] are used to validate the computational codes that we constructed. The experimental study consisted of a piston which oscillates at the head end of similar cavity. The acoustic pressure fields are recorded at several locations along the tube by capacitive pressure transducers. These transducers are connected to a computer through a built-in data acquisition system. The motor frequencies are changed using voltage varies. The range of the generated frequencies from the available devices is small enough to catch the first natural frequencies (250 Hz). The pressure traces at x = 0.125, ε = 0.045, and ω = 0.52 (f ` = 63.44 Hz) was recorded.The validated results deal with the acoustic fields in the chamber without injection from the sidewall, using the experimental, computational, and analytical approaches. The analytical results are obtained to describe acoustic and viscous flow dynamics in an impermeable duct when an end-wall disturbance is used to derive axial waves. The results are used to show that resonance between forced and propagating acoustic wave modes may be a source of large acoustic instability as well as to assert how nonlinear effects alter the energy in acoustic eigenfunction. The experimental results are compared with computational results to validate the numerical technique.The first set of the results shows a comparison between experimental and computational pressure traces atx = 0.125, ε = 0.045, and ω = 0.52 (f′ = 63.44Hz), as shown in Figure 2.Figure 2
Comparison between the experimental and computational results for the pressure-time history at (x) = 0.125, ε = 0.045, and ω = 0.52 (f′ = 63.44 Hz).The result illustrates the wave oscillations (quasi-steady wave motion) after many wave cycles. A stable or bounded solution is observed at this frequency. A reasonable comparison between the experimental and computational approaches is noticed. The deviations are directed to the linearity and complexity of the experimental test rig design. As a result, undesirable harmonics are seen with the experimental work due to some leaks in its connection between the circular portion of the piston and the duct shape. In the experimental study, the acoustic wave was generated using the cylindrical piston and the measurements were carried out in a rectangular cross-section duct. Therefore, the transition from circular to rectangular cross-sections results in a nonsmooth pressure wave.Moreover, a comparison between the current numerical solution and the analytical results is presented in Figure3. The results are illustrated for the pressure-time history at the midlength of the chamber, ε = 0.1 and ω = 1.0 (f′ = 159Hz). The deviation between the analytical and the computational approaches appears reasonable for the period between 0 < t < 20, while the deviation increases as time increases. This deviation may be due to nonlinearity effects.Figure 3
Comparison between the computational and analytical results for the pressure-time history at the midlength of the chamber,ε = 0.1 and ω = 1.0 (f′ = 159 Hz).Coupling effects between the burning process and the flow characteristics inside the SRM chamber are associated with small-amplitude pressure oscillations, which can be described mathematically in terms of appropriate eigenfunctions. The co-existence of forced acoustic modes and propagating modes has a significant role in the energy exchange mechanism and acoustic instability inside the chamber of an SRM. One of the primary objectives of this research is to illustrate how acoustic wave patterns in a duct are generated by a simple time-dependent boundary condition at the head end in terms of an initial-boundary value formulation. The numerical solution to the full parabolized nonlinear, unsteady, and compressible Navier- Stokes equations is integrated with the analytical solutions using the perturbation theory to the weakly nonlinear analysis of the acoustic waves to examine the complex acoustic/fluid dynamics mechanism and their impact on the combustion of real solid fuel. Moreover, both approaches are used to show that the resonance between forced and propagating acoustics wave modes may be the source of large acoustic instability.Lower and higher-order calculations are formulated to find a weakly nonlinear solution within the core flow, which defines the eigenfunctions’ acoustic damping time. The lowest order acoustic equations are valid for the limitM ⟶ 0 and Re >> 1.Our recent study by Hegab et al. [23] revealed that a weakly nonlinear theory is required to describe the acoustic evolution for a longer time. Thus, the perturbation analysis for higher-order calculations is used to find a weakly nonlinear solution within the core flow to provide additional verification for the acoustic damping time of the eigenfunctions.Moreover, the higher-order calculations revealed that the quadratic nonlinear convective terms could not be a source of intermodal energy transfer at the O (M2) approximation level. As a result, energy in eigenfunction modes cannot be damped by quadratic convection effects when “open box” eigenfunctions are the basis set for the Fourier series. At the same time, the O (M3) cubic nonlinearities are responsible for nonlinearization in the cavity of the SRM chamber. It follows that wave-wave interactions (nonharmonic behaviour of eigenfunction amplitudes) will occur over a time scale long compared to the acoustic time t = O (1).Generally, we may conclude that Figure3 supports these remarks above by comparing a finite-difference solution for the 2D, nonlinear, unsteady, and parabolized Navier–Stokes equations and the one-dimensional linear acoustic theory equations.Finally, the comparison of our constructed numerical codes with both the experimental and numerical solutions shows a reasonable agreement and reflects the ability of the FORTRAN programmed codes to run for many wave cycles with stable and reliable solutions. Moreover, several cases, beside the above ones, are performed to validate the constructed numerical codes in Section6.


## 6. Results and Discussion

The investigation of the flow field inside the chamber of the solid rocket motor (SRM) has been carried out. The current study emphasizes the acoustic-fluid dynamics interaction mechanism and the accompanying unsteady rotational fields. Moreover, the effect of generated turbulence on these unsteady rotational fields is considered for the first time. Several cases are performed to validate the constructed numerical codes as follows.


### 6.1. Case 1

The first set of results is considered to compare the numerical solution with the analytical solution for the pure acoustic flow field inside the chamber with end-wall disturbances only. Figure4 illustrates the comparison of the numerical results obtained using both the SIMPLE and PISO turbulence algorithms and the weakly nonlinear analytical solution at ω¯ = 1.0 at several axial locations inside the SRM chamber.Figure 4
Pressure and velocity-time histories at different axial locations (x/L = 0.25, 0.5, and 0.75): comparisons of analytical solution (dashed-dot) and numerical predictions using SIMPLE (solid) and PISO (dashed) algorithms at ω¯ = 1.0.
(a)(b)(c)The comparisons showed a reasonable agreement between the numerical and analytical results for both the axial velocity and pressure-time histories. Additionally, no noticeable difference is obtained in the numerical results when either the SIMPLE or PISO algorithm is used.However, an important feature of the numerical results for the axial velocity and pressure profiles is the appearance of higher harmonics oscillations, which represent the eigenfunction mode contributions at different time intervals. This interesting phenomenon can be attributed to the nonlinear terms included in the numerical simulation and is neglected in the analytical analysis. Moreover, it is clearly seen from the results that the numerical solution using either SIMPLE or PISO algorithms are matched together, while the analytical solution is noticeably deviated. The reason behind these discrepancies is related to the nonlinearity considerations with the numerical treatments than that with the weakly nonlinear consideration for the analytical derivations. The importance of including the nonlinear terms is evident for the damping and attenuation in the velocity and pressure oscillations for the results obtained from the numerical techniques than that of the quasi-steady oscillations from analytical analysis, where it can be observed.The resonance occurs at a forced frequency ofω¯=n+0.5π. The code is also validated at a relatively high-forced frequency, ω¯ = 1.4, which is very close to the resonance frequency. Only the SIMPLE algorithm is used in this case.Figure5 shows a comparison between the time histories of the predicted velocity and pressure and the analytical solution at near resonance frequency, ω¯=1.4.Figure 5
Pressure and velocity-time histories at different axial locations: comparisons of analytical solution and numerical predictions at near resonance frequency,ω¯=1.4.
(a)(b)(c)Beats are observed for both the numerical and analytical solutions atω¯=1.4. The agreement between the numerical and analytical results is found at the beginning of the time, up to t = 0.05 s. As time increases, the deviations between the predicted results and the analytical solution are recognized. This deviation may be attributed to the weakly nonlinear assumptions that are considered with the analytical solution. The importance of these comparisons is to indicate the ability of the developed numerical code to predict the full features of the generated acoustic field that consists of the fundamental and higher harmonics inside the SRM chamber at a wide range of frequencies either away from or very close to the resonance inside the cavity.


### 6.2. Case 2

The second case represents the simulated SRM chamber as a rectangular channel with a sidewall mass injection. As we mentioned, the sidewall mass injection mimics the burning rate of the real solid rocket motor propellant. This case was thoroughly validated in our previous work [8] using different turbulence models. The RSM v2-f and SST k-ω models were in a close agreement with the measured experimental data from [25].In the present paper, we only explain some important features of this flow configuration, especially the turbulence characteristics. Due to the continuous sidewall mass addition, the axial flow velocity increases towards the duct exit, and the flow is subjected to a strong streamline curvature, as shown in Figure6.Figure 6
Streamlines and velocity vectors of the internal flow field in a simulated solid rocket motor chamber with sidewall injection.
(a)(b)This velocity gradient and streamlined curvature results in an increase in turbulent kinetic energy production. As a result, the turbulent kinetic energy increases towards the duct exit, as shown in Figure7. Similar results are obtained by Hegab et al. [23].Figure 7
Contours of turbulent kinetic energy (m2/s2) in a simulated solid rocket motor chamber with sidewall injection.Therefore, the possibility of a laminar-to-turbulent transition exists. One way to identify the transition from laminar to turbulent is to examine the wall friction coefficient. Figure8 shows the friction coefficient obtained using different turbulent models. The friction coefficient, Cf, is calculated as follows:(19)Cf=2τwρUm2,where Um=∫0Hudy/H is the local mean streamwise velocity.Figure 8
Distribution of the local friction coefficient along the wall of the solid rocket motor chamber with sidewall injection.The results shown in Figure8 reveal that the flow is laminar near the head end and transitions to turbulence at a particular axial distance. Dunlap et al. [26], Hegab et al. [23], and Griffond et al. [27] reported similar observations. Only the v2-f turbulence model can predict this phenomenon accurately. Therefore, the v2-f model is the appropriate choice for flow with a sidewall mass injection, and it will be used for the last case. The effect of turbulence on the steady-state vorticity is shown in Figure 9, which shows that the vorticity is very high near the duct exit for turbulent flow solution. These higher values of vorticity can cause erosive burning in a real SRM.Figure 9
Steady-state vorticity distribution: (a) laminar; (b) turbulent.
(a)(b)


### 6.3. Case 3

To indicate the coupling between the induced acoustics and turbulence, a channel with both an end-wall disturbance and sidewall injection is simulated. The numerical simulation is performed in the case of a laminar regime (without including a turbulence model) and, in the case of turbulence, modelling via thev2-f model. First, a converged solution is obtained without an end-wall disturbance. Once the internal flow field within the motor chamber is well established, a sinusoidal velocity is applied at the head end, and the pressure and velocity histories are recorded.Figure10 shows a comparison of the time histories for the axial velocity and pressure for the laminar and turbulent flow.Figure 10
Pressure and velocity-time histories at different axial locations for case 3: comparisons of laminar and turbulent predictions.
(a)(b)(c)(d)The left-hand side represents the velocity-time history, while the pressure-time history is illustrated on the right side. The deviation between the laminar and turbulent solution for the velocity-time history has an apparent significant change as the axial distance increases downstream. This difference can be attributed to the presence of turbulent flow near the duct exit. An examination of the velocity-time history atx/L = .75 shows nonsinusoidal behaviour for the laminar flow. This can be attributed to the excitation of subharmonic waves, which may be damped in the turbulent flow due to the higher effective viscosity. Nasr et al. [14] and Javed and Chakraborty [15] reported similar behaviour.Unlike the velocity-time history, the pressure-time history reveals a reasonably significant deviation between the two solutions over the whole channel length. Thus, the pressure is higher in the turbulent solution than in laminar solutions. To explain this behaviour, the means and root mean squares of the pressure are presented in Figure11, which shows that these values are higher in a turbulent than in a laminar flow solution due to the higher pressure drop associated with the turbulent flow. The mean pressure fields are identical to those obtained from the steady solution.Figure 11
Comparison of laminar and turbulent solutions via mean and root mean square pressures.
(a)(b)To indicate the effect of turbulence on the wave motion, the fluctuation pressure,p′, is calculated and presented in Figure 12, which shows that the fluctuating pressure obtained from the laminar and turbulent solutions is very close near the head end (i.e., at x/L = 0.125 and 0.25). As the flow approaches the duct exit, the two solutions are entirely different. Thus, the turbulent solution damps the excitation from the lower harmonic waves. Therefore, the v2-f turbulence model can predict the flow field in both the laminar (near the head end) and turbulent (near the duct exit) regions with fair accuracy.Figure 12
Time history of pressure fluctuation for laminar and turbulent simulations.
(a)The surface profiles for the unsteady vorticity at different times for the laminar and turbulent flow are presented in Figure13. The figure indicates that the turbulent flow decreases the oscillations in the vorticity due to the higher associated viscosity. This behaviour is identified near the duct exit.Figure 13
Vorticity surface plot at different times for laminar: (a) and turbulent (b) predictions.
(a)(b)


## 6.1. Case 1

The first set of results is considered to compare the numerical solution with the analytical solution for the pure acoustic flow field inside the chamber with end-wall disturbances only. Figure4 illustrates the comparison of the numerical results obtained using both the SIMPLE and PISO turbulence algorithms and the weakly nonlinear analytical solution at ω¯ = 1.0 at several axial locations inside the SRM chamber.Figure 4
Pressure and velocity-time histories at different axial locations (x/L = 0.25, 0.5, and 0.75): comparisons of analytical solution (dashed-dot) and numerical predictions using SIMPLE (solid) and PISO (dashed) algorithms at ω¯ = 1.0.
(a)(b)(c)The comparisons showed a reasonable agreement between the numerical and analytical results for both the axial velocity and pressure-time histories. Additionally, no noticeable difference is obtained in the numerical results when either the SIMPLE or PISO algorithm is used.However, an important feature of the numerical results for the axial velocity and pressure profiles is the appearance of higher harmonics oscillations, which represent the eigenfunction mode contributions at different time intervals. This interesting phenomenon can be attributed to the nonlinear terms included in the numerical simulation and is neglected in the analytical analysis. Moreover, it is clearly seen from the results that the numerical solution using either SIMPLE or PISO algorithms are matched together, while the analytical solution is noticeably deviated. The reason behind these discrepancies is related to the nonlinearity considerations with the numerical treatments than that with the weakly nonlinear consideration for the analytical derivations. The importance of including the nonlinear terms is evident for the damping and attenuation in the velocity and pressure oscillations for the results obtained from the numerical techniques than that of the quasi-steady oscillations from analytical analysis, where it can be observed.The resonance occurs at a forced frequency ofω¯=n+0.5π. The code is also validated at a relatively high-forced frequency, ω¯ = 1.4, which is very close to the resonance frequency. Only the SIMPLE algorithm is used in this case.Figure5 shows a comparison between the time histories of the predicted velocity and pressure and the analytical solution at near resonance frequency, ω¯=1.4.Figure 5
Pressure and velocity-time histories at different axial locations: comparisons of analytical solution and numerical predictions at near resonance frequency,ω¯=1.4.
(a)(b)(c)Beats are observed for both the numerical and analytical solutions atω¯=1.4. The agreement between the numerical and analytical results is found at the beginning of the time, up to t = 0.05 s. As time increases, the deviations between the predicted results and the analytical solution are recognized. This deviation may be attributed to the weakly nonlinear assumptions that are considered with the analytical solution. The importance of these comparisons is to indicate the ability of the developed numerical code to predict the full features of the generated acoustic field that consists of the fundamental and higher harmonics inside the SRM chamber at a wide range of frequencies either away from or very close to the resonance inside the cavity.


## 6.2. Case 2

The second case represents the simulated SRM chamber as a rectangular channel with a sidewall mass injection. As we mentioned, the sidewall mass injection mimics the burning rate of the real solid rocket motor propellant. This case was thoroughly validated in our previous work [8] using different turbulence models. The RSM v2-f and SST k-ω models were in a close agreement with the measured experimental data from [25].In the present paper, we only explain some important features of this flow configuration, especially the turbulence characteristics. Due to the continuous sidewall mass addition, the axial flow velocity increases towards the duct exit, and the flow is subjected to a strong streamline curvature, as shown in Figure6.Figure 6
Streamlines and velocity vectors of the internal flow field in a simulated solid rocket motor chamber with sidewall injection.
(a)(b)This velocity gradient and streamlined curvature results in an increase in turbulent kinetic energy production. As a result, the turbulent kinetic energy increases towards the duct exit, as shown in Figure7. Similar results are obtained by Hegab et al. [23].Figure 7
Contours of turbulent kinetic energy (m2/s2) in a simulated solid rocket motor chamber with sidewall injection.Therefore, the possibility of a laminar-to-turbulent transition exists. One way to identify the transition from laminar to turbulent is to examine the wall friction coefficient. Figure8 shows the friction coefficient obtained using different turbulent models. The friction coefficient, Cf, is calculated as follows:(19)Cf=2τwρUm2,where Um=∫0Hudy/H is the local mean streamwise velocity.Figure 8
Distribution of the local friction coefficient along the wall of the solid rocket motor chamber with sidewall injection.The results shown in Figure8 reveal that the flow is laminar near the head end and transitions to turbulence at a particular axial distance. Dunlap et al. [26], Hegab et al. [23], and Griffond et al. [27] reported similar observations. Only the v2-f turbulence model can predict this phenomenon accurately. Therefore, the v2-f model is the appropriate choice for flow with a sidewall mass injection, and it will be used for the last case. The effect of turbulence on the steady-state vorticity is shown in Figure 9, which shows that the vorticity is very high near the duct exit for turbulent flow solution. These higher values of vorticity can cause erosive burning in a real SRM.Figure 9
Steady-state vorticity distribution: (a) laminar; (b) turbulent.
(a)(b)


## 6.3. Case 3

To indicate the coupling between the induced acoustics and turbulence, a channel with both an end-wall disturbance and sidewall injection is simulated. The numerical simulation is performed in the case of a laminar regime (without including a turbulence model) and, in the case of turbulence, modelling via thev2-f model. First, a converged solution is obtained without an end-wall disturbance. Once the internal flow field within the motor chamber is well established, a sinusoidal velocity is applied at the head end, and the pressure and velocity histories are recorded.Figure10 shows a comparison of the time histories for the axial velocity and pressure for the laminar and turbulent flow.Figure 10
Pressure and velocity-time histories at different axial locations for case 3: comparisons of laminar and turbulent predictions.
(a)(b)(c)(d)The left-hand side represents the velocity-time history, while the pressure-time history is illustrated on the right side. The deviation between the laminar and turbulent solution for the velocity-time history has an apparent significant change as the axial distance increases downstream. This difference can be attributed to the presence of turbulent flow near the duct exit. An examination of the velocity-time history atx/L = .75 shows nonsinusoidal behaviour for the laminar flow. This can be attributed to the excitation of subharmonic waves, which may be damped in the turbulent flow due to the higher effective viscosity. Nasr et al. [14] and Javed and Chakraborty [15] reported similar behaviour.Unlike the velocity-time history, the pressure-time history reveals a reasonably significant deviation between the two solutions over the whole channel length. Thus, the pressure is higher in the turbulent solution than in laminar solutions. To explain this behaviour, the means and root mean squares of the pressure are presented in Figure11, which shows that these values are higher in a turbulent than in a laminar flow solution due to the higher pressure drop associated with the turbulent flow. The mean pressure fields are identical to those obtained from the steady solution.Figure 11
Comparison of laminar and turbulent solutions via mean and root mean square pressures.
(a)(b)To indicate the effect of turbulence on the wave motion, the fluctuation pressure,p′, is calculated and presented in Figure 12, which shows that the fluctuating pressure obtained from the laminar and turbulent solutions is very close near the head end (i.e., at x/L = 0.125 and 0.25). As the flow approaches the duct exit, the two solutions are entirely different. Thus, the turbulent solution damps the excitation from the lower harmonic waves. Therefore, the v2-f turbulence model can predict the flow field in both the laminar (near the head end) and turbulent (near the duct exit) regions with fair accuracy.Figure 12
Time history of pressure fluctuation for laminar and turbulent simulations.
(a)The surface profiles for the unsteady vorticity at different times for the laminar and turbulent flow are presented in Figure13. The figure indicates that the turbulent flow decreases the oscillations in the vorticity due to the higher associated viscosity. This behaviour is identified near the duct exit.Figure 13
Vorticity surface plot at different times for laminar: (a) and turbulent (b) predictions.
(a)(b)


## 7. Conclusions

The internal flow field in an SRMC with compressible flow transpiration and end-wall disturbance has been numerically investigated. Three different cases are considered. The first case of the study is directed to pure acoustic fields due to end-wall disturbance, while the second case considered a steady-state rotational flow due to a steady sidewall injection. The last set of the results was performed to account for the complete acoustic-fluid dynamics interaction generated from an end-wall disturbance and sidewall injection. The numerical results are validated through a comparison with the analytical analysis and the available experimental measurements for the first and the second case, respectively. The comparisons showed a general fair agreement for the axial velocity and pressure profiles at different locations in the SRMC. The effect of including the nonlinear terms could be visibly observed in predicting higher harmonics oscillation in the numerical results rather than the analytical results.In the third case, the effect of including turbulence models in predicting the internal flow characteristics is explained by a comparison between the laminar and turbulent simulation. The comparison showed coupling between the generated unsteady vorticity arising from acoustics-fluid dynamics interaction and the turbulence intensity that exists everywhere in the chamber. Hegab and Kassoy [25] extensively studied the former and showed how the acoustic-fluid interaction might affect the burning rate and, in turn, has a significant effect on erosive burning.In the current study, there is a novel essential coupling, considering the turbulence existence accompanying the acoustic-fluid dynamics interaction. This essential coupling may add another factor that controls erosive burning in a real SRM. In general, more intensive computational work is needed to extend the current research to account for the combustion of real SRM propellant.In the context of propellant combustion, the intensive unsteady vorticity generation and the associated shear stresses can be anticipated to have a direct impact on the stability of burning the real solid rocket propellant. Moreover, more intensive computations are required to account for the existence of the convergent-divergent nozzle along with the end-wall disturbances.



---

*Source: 9978102-2021-06-16.xml*