The Submerged Nozzle Damping Characteristics in Solid Rocket Motor

Abstract

In this paper, the effects of the geometry of a submerged nozzle on the nozzle damping characteristics are studied numerically. Firstly, the numerical method is verified by the previous experimental data. Then, the mesh sensitivity analysis and the monitor position independence analysis are carried out. Thirdly, the effects of nozzle geometry on nozzle damping are systematically studied, and focuses are placed on the cavity size, convergent angle and divergent angle. The pulse decay method is utilized to evaluate the nozzle decay coefficient. Several important results are obtained: the submerged cavity with large volume leads to low frequency acoustic oscillations in the combustion chamber and corresponds to a small nozzle decay coefficient; then, as the nozzle convergent angle is decreased, the nozzle decay coefficient is increased. In addition, the nozzle divergent angle has a trivial effect on the nozzle decay coefficient; and lastly, the effects of the temperature on the nozzle damping capability are conducted. The results show that an increase of the working temperature leads to an increase of the nozzle decay coefficient; therefore, the damping force is increased.

1. Introduction

Combustion instability is an unwanted phenomenon in solid rocket motors [1]. It is generally raised by the coupling among the unsteady combustion process, the acoustic waves and the fluid dynamics [2]. It may lead to large-amplitude thrust oscillation, making the rocket motors unable to work normally or even to explode [3,4,5]. To investigate the primary mechanisms of combustion instability, extensive research has been performed [6,7,8,9,10]. According to the linear stability analysis, the thermoacoustic system response is determined by the sum of acoustic gain and damping factors [5]. The gain factors are mainly caused by vortex-acoustic coupling [11], thermo-acoustic coupling [12,13], and the distributed combustion [14], etc., while the damping factors mainly include the nozzle damping [15], particle damping [16], structural damping, and so on. When the acoustic energy gain is larger than the acoustic energy damping, the small disturbance in the combustion chamber is amplified, resulting in pressure oscillation and combustion instability. On the contrary, when the acoustic energy damping is larger than the acoustic energy gain, the disturbance in the combustion chamber is dissipated, and the solid rocket motor works stably. Therefore, increasing the acoustic energy damping is an effective method to suppress combustion instability. Among the main damping factors, nozzle damping can contribute more than 50% of the acoustic energy loss in the combustion chamber, where the acoustic energy is transmitted and radiated through the nozzle. Thus, studying the nozzle damping characteristics is very necessary and important.

Since World War II, significant achievements have been done on nozzle damping. Tsien [17] calculated the response of a choked nozzle under the influence of axial velocity perturbation and analyzed the transfer function of rocket nozzles. Crocco et al. extended Tsien’s work and introduced the nozzle admittance theory to study the influence of the nozzle on acoustic oscillations [18,19]. Crocco assumed that the mean flow is one-dimensional and that the acoustic oscillation in the nozzle is three-dimensional; in addition, he [20] first proposed the short-nozzle theory, in which the short nozzle refers to the nozzle whose convergent section is small compared to the chamber length. Zinn and Janardan et al. theoretically predicted the nozzle admittance based on the short-nozzle theory and experimentally studied the damping of axial instabilities of small-scale nozzles under cold-flow conditions [15,21,22,23,24,25]. Good agreement between the experimental data and the theoretical results are obtained. Buffum et al. [26,27] carried out several subscale tests to evaluate the nozzle decay coefficient, and three methods, including the steady-state resonance method, the pulse decay method, and the steady-state wave decay method, were used. Anthoine et al. [28,29] studied the effect of the nozzle cavity on a solid rocket motor with a submerged nozzle, both experimentally and numerically, and found that the amplitude of pressure oscillation is approximately linear with the nozzle cavity volume. Javed and Chakraborty [30] validated their computational model based on the experimental case reported by Buffum et al., then predicted the nozzle decay coefficient of a solid rocket motor with composite propellant. Su and Sun et al. [31,32,33] numerically studied the influence of nozzle geometries on nozzle damping via a steady-state wave decay method and evaluated the nozzle damping characteristics. Chen et al. [34] investigated the decay characteristics of pressure oscillation excited by pulses between finocyl-grain and axil-grain combustion chambers and indicated that the axil-grain is more stable than the finocyl-grain.

The present work is based on a solid rocket motor with a submerged nozzle, and the pulse decay method is used to study the effects of the geometry of a submerged nozzle and working temperature on the nozzle damping characteristics. In this study, the ANSYS software is utilized [35]. Specifically, the acoustic characteristics are calculated using ANSYS Mechanical, and the flow field is obtained using ANSYS Fluent. Firstly, the numerical method is verified by comparing the calculated results and experimental data. Then, the mesh sensitivity analysis and the monitor position independence analysis are carried out. Finally, the effects of nozzle geometry and working temperature on nozzle damping are studied and discussed in this paper. This work aims at providing information to help improve the passive control of the solid rocket motor with the submerged nozzle.

2. Analytical Model and Computational Methods

2.1. Theoretical Foundation

According to the liner stability analysis of combustion instability, the coefficient α is used to predict the growth or decay of a small disturbance in the combustion chamber, and further evaluate the stability of the solid rocket motor. If the coefficient $\mathsf{\alpha }>0$, the small disturbance will lead to pressure oscillations in the combustion chamber and the motor is unstable; if $\mathsf{\alpha }<0$, the small disturbance will be dissipated, and the motor is stable [36]. The coefficient α is expressed by the sum of all gain and damping coefficients, as follows:

$$\mathsf{\alpha }={\alpha }_{PC}+{\alpha }_{VC}+{\alpha }_{DC}+{\alpha }_N+{\alpha }_P+{\alpha }_{MF}+{\alpha }_G+{\alpha }_W+{\alpha }_{ST}$$

(1)

where the coefficients on the right side of the equation are the pressure coupling response coefficient, velocity coupling response coefficient, distributed combustion gain coefficient, nozzle damping coefficient, particle damping coefficient, average flow coefficient, gas phase damping coefficient, wall damping coefficient, and structural damping coefficient, respectively.

Based on previous studies [37,38], the linear stabilities of solid rocket motors are characterized by the exponential decay or growth rate of pressure oscillations as follows:

$$p^{\prime }=p_0{e}^{i\omega t+\alpha t}$$

(2)

where $p_0$ denotes the initial pressure amplitude, and $\omega$ is the radian frequency.

In this paper, the pulse decay method is numerically applied to evaluate the nozzle damping characteristics. During the transient calculation process, a sudden pulse is exerted on the combustion chamber, and then the decay coefficient can be obtained by processing the pressure oscillation decay curve. Suppose that there is no other source of acoustic energy existing in the chamber and there is no acoustic energy consumption. Thus, the acoustic energy is only affected by the nozzle and the behavior of the pressure oscillation in the chamber can be expressed as follows:

$$p^{\prime }=p_0{e}^{{\alpha }_{N}t}$$

(3)

Therefore, the nozzle decay coefficient ${\alpha }_N$ can be obtained from

$${\alpha }_N=\frac{\ln p_2-\ln p_1}{t_2-t_1}$$

(4)

where $p_1$ and $p_2$ are the values of two pressure peaks, and $t_1$ and $t_2$ are the respective time instants.

2.2. Finite Element Formulation of the Acoustic Wave Equation

In this study, the acoustic finite element analysis (FEA) is used for the acoustic simulation. According to the previous studies [34,39], the fluid momentum equations and continuity equations are simplified to get the acoustic wave equation by using the following assumptions:

• The fluid is compressible (density changes due to pressure variations);
• There is no body force;
• The fluid is irrotational;
• The pressure distribution in the combustion chamber is small;
• There is no mean flow;
• The gas is ideal, adiabatic, and reversible.

The linearized continuity equation is:

$$\nabla ·{\overrightarrow{v}}_a=-\frac{1}{{\rho }_0{c}^2}\frac{\partial p_a}{\partial t}+\frac{Q}{{\rho }_0}$$

(5)

The linearized momentum equation is:

$$\frac{\partial {\overrightarrow{v}}_a}{\partial t}=-\frac{1}{{\rho }_0}\nabla p_a+\frac{4\mu }{3{\rho }_0}\nabla \left(-\frac{1}{{\rho }_0{c}^2}\frac{\partial p_a}{\partial t}+\frac{Q}{{\rho }_0}\right)$$

(6)

where ${\overrightarrow{v}}_a$ is the acoustic velocity, $p_a$ denotes the acoustic pressure, $c$ represents the acoustic speed in fluid medium, $Q$ is the mass source, ${\rho }_0$ represents the mean fluid density, and $\mu$ is the dynamic viscosity.

By using Equations (5) and (6) and eliminating ${\overrightarrow{v}}_a$, the acoustic wave equation is given by:

$$\nabla ·\left(\frac{1}{{\rho }_0}\nabla \mathrm{p}\right)-\frac{1}{{\rho }_0{c}^2}\frac{{\partial }^{2}p}{\partial t^2}+\nabla ·\left[\frac{4\mu }{3{\rho }_0}\nabla \left(\frac{1}{{\rho }_0{c}^2}\frac{\partial p}{\partial t}\right)\right]=-\frac{\partial }{\partial t}\left(\frac{Q}{{\rho }_0}\right)+\nabla ·\left[\frac{4\mu }{3{\rho }_0}\nabla \left(\frac{Q}{{\rho }_0}\right)\right]$$

(7)

By integrating Equation (7) over the whole volume [35], the finite element formulation of acoustic wave equation is expressed as:

$$\begin{array}{c}\underset{{\mathsf{\Omega }}_F}{\overset{ }{{\displaystyle \iiint }}}\frac{1}{{\rho }_0{c}^2}w\frac{{\partial }^{2}p}{\partial t^2}\mathrm{d}V+\underset{{\mathsf{\Omega }}_F}{\overset{ }{{\displaystyle \iiint }}}\nabla \mathrm{w}·\left(\frac{4\mu }{3{\rho }_0^2{c}^2}\nabla \frac{\partial p}{\partial t}\right)\mathrm{d}V+\underset{{\mathsf{\Omega }}_F}{\overset{ }{{\displaystyle \iiint }}}\nabla \mathrm{w}·\left(\frac{1}{{\rho }_0}\nabla p\right)\mathrm{d}V\\ +\underset{{\mathsf{\Gamma }}_F}{\overset{ }{{\displaystyle ∯}}}w\widehat{n}·\frac{{\partial }^2{\overrightarrow{u}}_F}{\partial t^2}\mathrm{d}S=\underset{{\mathsf{\Omega }}_F}{\overset{ }{{\displaystyle \iiint }}}w\frac{1}{{\rho }_0}\frac{\partial Q}{\partial t}\mathrm{d}V+\underset{{\mathsf{\Omega }}_F}{\overset{ }{{\displaystyle \iiint }}}\nabla \mathrm{w}·\left(\frac{4\mu }{3{\rho }_0^2}\nabla Q\right)\mathrm{d}V\end{array}$$

(8)

where $w$ is the testing function, $V$ represents the volume of acoustic domain ${\mathsf{\Omega }}_F$, $S$ represents the surface of acoustic domain boundary ${\mathsf{\Gamma }}_F$, $\widehat{n}$ denotes the outward normal unit vector to the boundary, and ${\overrightarrow{u}}_F$ is the displacement of the fluid particle.

The finite element–approximating shape functions for the spatial variation of the pressure and displacement components are given by:

$$\mathrm{P}={\left\{N\right\}}^T\left\{P_e\right\},\mathrm{u}={\left\{N^{\prime }\right\}}^T\left\{u_e\right\}$$

(9)

where $\left\{N\right\}$ represents the element shape function for pressure, $\left\{N^{\prime }\right\}$ represents the shape function for displacements, $\left\{P_e\right\}$ is the nodal pressure vector, and $\left\{u_e\right\}$ is the displacement component vector.

The second time derivative of the pressure and displacement components and the virtual change in the pressure can be expressed as follows:

$$\frac{{\partial }^{2}p}{\partial t^2}={\left\{N\right\}}^T\left\{{\ddot{P}}_e\right\}$$

(10)

$$\frac{{\partial }^2}{\partial t^2}\left\{u\right\}={\left\{N^{\prime }\right\}}^T\left\{{\ddot{u}}_e\right\}$$

(11)

$$\delta P={\left\{N\right\}}^T\left\{\delta P_e\right\}$$

(12)

Therefore, the wave Equation (8) can be rewritten as:

$$\begin{array}{c}\underset{{\mathsf{\Omega }}_F}{\overset{ }{{\displaystyle \iiint }}}\frac{1}{{\rho }_0{c}^2}\left\{N\right\}{\left\{N\right\}}^T\mathrm{d}V\left\{{\ddot{p}}_e\right\}+\underset{{\mathsf{\Omega }}_F}{\overset{ }{{\displaystyle \iiint }}}\frac{4\mu }{3{\rho }_0^2{c}^2}{\left[\nabla N\right]}^T\left[\nabla N\right]\mathrm{d}V\left\{{\dot{p}}_e\right\}\\ +\underset{{\mathsf{\Omega }}_F}{\overset{ }{{\displaystyle \iiint }}}\frac{1}{{\rho }_0}{\left[\nabla N\right]}^T\left[\nabla N\right]\mathrm{d}V\left\{p_e\right\}+\underset{{\mathsf{\Gamma }}_F}{\overset{ }{{\displaystyle ∯}}}\left\{N\right\}{\left\{n\right\}}^T{\left\{N^{\prime }\right\}}^T\mathrm{d}S\left\{{\ddot{u}}_{e,F}\right\}\\ =\underset{{\mathsf{\Omega }}_F}{\overset{ }{{\displaystyle \iiint }}}\frac{1}{{\rho }_0}\left\{N\right\}{\left\{N\right\}}^T\mathrm{d}V\left\{\dot{q}\right\}+\underset{{\mathsf{\Omega }}_F}{\overset{ }{{\displaystyle \iiint }}}\frac{4\mu }{3{\rho }_0^2}{\left[\nabla N\right]}^T\left[\nabla N\right]\mathrm{d}V\left\{q\right\}\end{array}$$

(13)

where $\left\{n\right\}$ denotes outward normal vector at the fluid boundary, $\left\{q\right\}$ is the node mass source vector, and $\left\{\dot{q}\right\}$ represents the first time derivative of a nodal mass source vector.

Equation (13) can be written in matrix notation and a discretized wave equation is obtained:

$$\left[M_F\right]\left\{{\ddot{p}}_e\right\}+\left[C_F\right]\left\{{\dot{p}}_e\right\}+\left[K_F\right]\left\{p_e\right\}+{\rho }_0{\left[R\right]}^T\left\{{\ddot{u}}_{e,F}\right\}=\left\{f_F\right\}$$

(14)

where $\left[M_F\right]$ represents the acoustic fluid mass matrix, $\left[C_F\right]$ represents the acoustic fluid damping matrix, $\left[K_F\right]$ represents the acoustic fluid stiffness matrix, ${\left[R\right]}^T$ represents the acoustic fluid boundary matrix, and $\left\{f_F\right\}$ represents the acoustic fluid load vector. To solve the governing Equation (14), the ANSYS Mechanical software is utilized.

2.3. Governing Equations for the Flow-Field

In the flow field simulation, the flow is assumed to be compressible and rotational. The Reynolds averaging method is used, and the instantaneous variable $f$ is decomposed into the mean component $\overline{f}$ and the fluctuating component $f^{\prime }$, i.e., $f=\overline{f}+f^{\prime }$. The continuity, momentum, and energy conservation equations are given as follows:

$$\frac{\partial \overline{\rho }}{\partial t}+\frac{\partial }{\partial x_i}\left(\overline{\rho }{\overline{u}}_i\right)=0$$

(15)

$$\frac{\partial }{\partial t}\left(\overline{\rho }{\overline{u}}_i\right)+\frac{\partial }{\partial x_j}\left(\overline{\rho }{\overline{u}}_i{\overline{u}}_j\right)=-\frac{\partial \overline{p}}{\partial x_i}+\frac{\partial }{\partial x_j}\left({\overline{\tau }}_{ij}-\overline{\rho u_i^{\prime }{u}_j^{\prime }}\right)$$

(16)

$$\begin{array}{l}\frac{\partial }{\partial t}\left(\overline{\rho }\overline{E}\right)+\frac{\partial }{\partial x_j}\left(\overline{\rho }\overline{E}{\overline{u}}_j\right)\\ =-\frac{\partial }{\partial x_j}\left(\overline{p}{\overline{u}}_j\right)+\frac{\partial }{\partial x_j}\left[\left({\overline{\tau }}_{ij}-\overline{\rho u_i^{\prime }{u}_j^{\prime }}\right){\overline{u}}_i\right]+\frac{\partial }{\partial x_j}\left({\overline{\tau }}_{ij}-\overline{\rho u_i^{\prime }{u}_j^{\prime }}\right)\\ +\frac{\partial }{\partial x_j}\left(\lambda \frac{\partial \overline{T}}{\partial x_j}\right)\end{array}$$

(17)

where

$${\overline{\tau }}_{ij}=\mu \left(\frac{\partial {\overline{u}}_i}{\partial x_j}+\frac{\partial {\overline{u}}_j}{\partial x_i}-\frac{2}{3}{\delta }_{ij}\frac{\partial {\overline{u}}_l}{\partial x_l}\right)$$

(18)

$$-\overline{\rho u_i^{\prime }{u}_j^{\prime }}={\mu }_t\left(\frac{\partial {\overline{u}}_i}{\partial x_j}+\frac{\partial {\overline{u}}_j}{\partial x_i}-\frac{2}{3}{\delta }_{ij}\frac{\partial {\overline{u}}_k}{\partial x_k}\right)-\frac{2}{3}{\delta }_{ij}\overline{\rho }k$$

(19)

where $\overline{u}$ represents the mean velocity, $\overline{\rho }$ denotes the mean fluid density, $\overline{p}$ denotes the mean fluid pressure, $\overline{E}=\overline{\mathrm{e}}+\frac{1}{2}\overline{u_i{u}_i}$ and $\overline{\mathrm{e}}$ is the internal energy, $\lambda$ represents the heat diffusion coefficient, $\overline{T}$ denotes the mean temperature, ${\overline{\tau }}_{ij}$ represents the viscous stress tensor, $\mu$ is the molecular viscosity, $-\overline{\rho u_i^{\prime }{u}_j^{\prime }}$ represents the Reynolds stress, ${\mu }_t$ is the turbulent viscosity, $k=\frac{1}{2}\overline{u_i^{\prime }{u}_i^{\prime }}$ denotes the turbulence kinetic energy, and ${\delta }_{ij}$ represents the Kronecker symbol: if $i=j$, then ${\delta }_{ij}=1$; otherwise, ${\delta }_{ij}=0$.

The ideal gas law is added to the governing equations as follows:

$$p=\rho RT$$

(20)

where $R$ is the specific gas constant of the gas.

The standard $k-\epsilon$ turbulence model is used in this study, and therefore, the turbulent viscosity (${\mu }_t$) is computed by combining turbulence kinetic energy ($k$) and its dissipation rate ($\epsilon$). The turbulence kinetic energy and its dissipation rate are obtained from the following transport equations:

$$\frac{\partial }{\partial t}\left(\overline{\rho }k\right)+\frac{\partial }{\partial x_i}\left(\overline{\rho }k{\overline{u}}_i\right)=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_j}\right]+G_k-\rho \epsilon$$

(21)

$$\frac{\partial }{\partial t}\left(\overline{\rho }\epsilon \right)+\frac{\partial }{\partial x_i}\left(\overline{\rho }\epsilon {\overline{u}}_i\right)=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_{\epsilon }}\right)\frac{\partial \epsilon }{\partial x_j}\right]+C_{1\epsilon }\frac{\epsilon }{k}{G}_k-C_{2\epsilon }\rho \frac{{\epsilon }^2}{k}$$

(22)

$${\mu }_t=\rho C_{\mu }\frac{k^2}{\epsilon }$$

(23)

where $G_k$ represents the generation of turbulence kinetic energy due to the mean velocity gradients, $C_{1\epsilon }=1.44$, $C_{2\epsilon }=1.92$, and $C_{\mu }=0.09$, ${\sigma }_k$ and ${\sigma }_{\epsilon }$ are the turbulent Prandtl numbers for $k$ and $\epsilon$, respectively; they are set as ${\sigma }_k=1.0$ and ${\sigma }_{\epsilon }=1.3$.

In the flow simulation, the CFD solver ANSYS FLUENT is used, and a bounded central-differencing scheme is selected for momentum and energy equations. An implicit second-order method is used in an unsteady solver. The steady calculation is firstly carried out to get good convergent solutions. Then, a transient calculation is performed and the timestep is 1 × 10⁻⁵ s to capture the pressure oscillations. When good convergent results are again obtained, the mass flow rate at the inlet is increased suddenly in a short time as a pulse via user-defined function (UDF). After that, the mass flow rate turns back to the original value and the data on pressure oscillation is recorded.

2.4. Computational Domain and Boundary Conditions

In the acoustic modal analysis of the combustion chamber, with the exception of the throat and the nozzle divergent section, all the surfaces are defined as rigid walls with no displacement boundaries. In the flow-field simulation, the computational geometry and the boundary conditions of the solid rocket motor model with a submerged nozzle are illustrated in Figure 1. The total length of the combustion chamber is ${\mathrm{L}}_{\mathrm{c}}=0.96$ m and its radius is ${\mathrm{R}}_{\mathrm{c}}=0.1$ m. The nozzle throat radius is ${\mathrm{R}}_{\mathrm{t}}=0.02$ m and the cavity radius is expressed by ${\mathrm{R}}_{\mathrm{cavity}}$. ${\mathrm{L}}_{\mathrm{cavity}}$ represents the length of the cavity which is equal to ${\mathrm{R}}_{\mathrm{cavity}}$, and $\mathrm{h}=0.02$ m is constant. Therefore, the larger ${\mathrm{R}}_{\mathrm{cavity}}$ leads to the larger cavity volume. The left boundary is configured to be mass flow inlet, and gases are injected axially. The nozzle exit is defined as a pressure outlet, and the wall is characterized by no-slip boundary condition. The gas is simplified as an ideal gas and its physical properties used in the flow field simulation are shown in

Table 1. Physical properties for the cold flow.

Physical Properties	Value
$\mu$	Sutherland law
$C_p$	1700 J/(kg · K)
$M$	26 g/mol
$\gamma$	1.4
$\mathsf{\lambda }$	0.0242 w/(m · K)

.

Six monitoring points (as shown in Figure 1) are picked to record the pressure oscillations. The first five points are placed along the axis of symmetry and the last one is located in the submerged cavity. The coordinates of the monitoring point location are listed in

Table 2. Coordinates of monitoring point location.

Coordinate	P1	P2	P3	P4	P5	P6
X/m	0.001	0.240	0.480	0.720	0.960	0.960
Y/m	0.001	0.001	0.001	0.001	0.001	0.120

.

3. Verification of the Numerical Method

3.1. Validation of FEA in the Acoustic-Field Simulation

In this paper, the acoustic FEA method is performed to study the acoustic characteristics of the combustion chamber. To verify the accuracy of the method, the acoustic modes of the VKI model are simulated, and the results are compared with the experimental data [29]. As shown in

Table 3. Comparison between FEA results and experimental data.

Mode Order	FEA Results/Hz	Experimental Data/Hz	Errors
1st	417	408	2.1%
2nd	835	874	4.6%
3rd	1284	1285	0.1%
4th	1736	1744	0.5%

, the errors between the predicted results and experimental data for the first four acoustic frequencies are 2.1%, 4.6%, 0.1%, and 0.5%, respectively. Therefore, the frequencies calculated by this method are highly consistent with the experimental data. Although the model used in the verification is different from the model studied in this paper, the conclusions prove the accuracy of the FEA method used in this study.

3.2. Validation of the Flow Models

To confirm the reliability of the numerical method, the pulse decay method is utilized, and the experimental and theoretical results from Buffum tests [26] are introduced in Figure 2.

It can be found that the decay coefficients from experimental measurements are much higher than those from the numerical and theoretical methods. For the experimental results, except the nozzle losses, the decay coefficient also covers other damping contributions, such as viscous wall losses, heat conduction losses, and so on. In another words, the experimental values are equivalent to the sum of multiple damping coefficients, not just the nozzle damping coefficient. Although there are some obvious deviations between the experimental data and numerical results, the trends of the numerical prediction are closer to the experimental values than to the theoretical values. In summary, the method provides an effective way to evaluate the nozzle damping characteristics.

3.3. Mesh Sensitivity Analysis

In the flow-field simulation, to obtain high mesh qualities, the quadrilateral grids are generated in the computational domain. The grids are clustered toward the side wall, the nozzle wall, and the nozzle throat to capture large flow gradients, as shown in Figure 3.

To confirm the independence of the numerical results on the mesh qualities, it is necessary to perform the mesh sensitivity analysis before the systematical studies [33]. Here, six different computational meshes have been employed for the grid validation, ranging from 10,000 to 250,000 cell numbers. To evaluate the nozzle damping capability, the decay coefficient is used and therefore treated as the mesh quality evaluation criterion. The nozzle decay coefficients under different mesh cell numbers are calculated and the results are illustrated in Figure 4. As the mesh cell number increases, it can be seen that the decay coefficient increases firstly, and when the mesh cell number exceeds 150,000, the decay coefficient almost levels off. To save the computational time and maintain satisfactory accuracy, the mesh cell number of 150,000 is chosen in the next studies.

3.4. Independence Verification of Monitoring Points

From the acoustic point of view, the combustion chamber of the solid rocket motor can be treated as a closed-closed structure, and the two ends of the combustion chamber are the pressure antinodes [40]. Therefore, the pressure oscillation of the first-order acoustic mode at the two ends can reach the maximum at the same time but there is a phase difference of 180°. The first-order acoustic pressure oscillations and the phase difference calculated by the numerical method are illustrated in Figure 5. It can be seen that the pressure peak of point 1 (or point 2) is ahead of one half of the oscillation period compared with point 5 (or point 4 and point 6) because the pressure node is located between point 2 and point 4. It indicates that the pressure oscillation of point 1 and point 5 have an opposite phase which is in good agreement with the theoretical analysis.

To determine the effect of monitoring points on the nozzle decay coefficient, the pressure oscillations recorded at different monitoring points (excluding point 3) were processed and the results were shown in

Table 4. The decay coefficient of different monitoring points.

Point	P1	P2	P3	P4	P5	P6
${\alpha }_N/s^{-1}$	−8.59	−8.55	--	−8.63	−8.65	−8.65

. There are slight differences between the five values, which are mainly caused by artificial errors. Considering that point 1 and point 5 (or point 6) can describe the pressure oscillation best and that the pressure amplitude of point 1 is maximum, so the decay coefficient at point 1 is chosen as the final nozzle decay coefficient in the next studies.

4. Results and Analysis

In this paper, to study the effects of nozzle geometry and operating conditions on nozzle damping, 13 different cases with different nozzle parameters are designed for comparison. Detailed information is listed in

Table 5. Different cases studied in this paper.

Case	Cavity Radius/mm	Convergent Angle/deg	Divergent Angle/deg	Working Temperature/K
C0	40	50	30	300
C1	45	50	30	300
C2	50	50	30	300
C3	55	50	30	300
C4	40	30	30	300
C5	40	40	30	300
C6	40	60	30	300
C7	40	50	10	300
C8	40	50	20	300
C9	40	50	40	300
C10	40	50	30	900
C11	40	50	30	1500
C12	40	50	30	2100

.

4.1. Acoustic Field in Chamber

The three-dimensional model of case C0 is studied by the FEA method in this section. The first four acoustic mode shapes are shown in Figure 6. It can be seen that the two ends of the combustion chamber are the antinodes of the acoustic pressure and that the submerged nozzle has little effect on the distribution of the acoustic pressure. The effect of the submerged cavity on the acoustic field is to move the position of the acoustic pressure antinode from the nozzle throat to the end of the cavity, which is equivalent to extending the axial length of the acoustic chamber.

Based on the FEA method, the acoustic frequencies are obtained. Then, the spectra of the pressure signals recorded in the flow-field simulation are computed with a fast Fourier transform (FFT) and the frequency domain information can be obtained. The first fourth-order acoustic frequencies of different cases by the two methods are shown in Figure 7. It is evident that the errors of the results are very small except for the fourth-order acoustic frequencies of case C3, and the frequencies computed with CFD are higher than the frequencies obtained by FEA, which may be caused by the modification of the flow inside the combustor. Besides, according to the acoustic theory, the axial inherent acoustic frequency is inversely proportional to the length of the acoustic chamber. Thus, when the size of the submerged cavity increases gradually, the equivalent length of the acoustic chamber increases, and the frequency then decreases.

4.2. The Effect of Cavity Size on Nozzle Damping

The submerged nozzle is widely used in solid rocket motors with limited length, which can reduce the length of the combustion chamber and the mass of the motor. However, the submerged nozzle can also increase the risk of combustion instability. Figure 8 depicts the streamlines at the end of the motor, and it can be seen that there are three distinct vortices in the cavity, which have an important effect on the stability of the motor [41]. On the one hand, in the view of the acoustic gain, the structure of the submerged cavity is more prone to generate some eddies, and the energy transfer between the vorticity field and the acoustical field may occur due to the vortex-acoustic coupling. On the other hand, in the view of the acoustic damping, the submerged cavity weakens the nozzle damping.

To study the effect of cavity size on the nozzle damping characteristics, based on case C0, three other cases (C1, C2, and C3) with different cavity sizes are chosen for comparison. Using the same numerical method described above, the pressure decay curves after the mass pulse are illustrated in Figure 9. It can be found that although the pulse intensity for the different cases is consistent, the amplitude of the pressure oscillation is different. The larger size of cavity leads to the more severe pressure oscillation which had been verified by Anthoine’s experiment [28]. In addition, the decay rates of different cases can be observed visually; the specific values of the decay coefficient are illustrated in Figure 10.

Figure 10a depicts the fitting line of $\ln \left(p^{\prime }\right)-t$, where $p^{\prime }$ is the pressure peak value of the pressure decay curve in Figure 9, and the slope of the fitting line $k$ represents the value of the nozzle decay coefficient [42], which is calculated and illustrated in Figure 10b. It can be seen that the absolute value of the nozzle decay coefficient decreases with the increasing of cavity size, namely, the submerged cavity has negative effects on the nozzle damping. That is to say, the increase in cavity size will reduce the loss of acoustic energy, which is not beneficial to the improvement of nozzle damping. On average, for every 12.5% increase in the cavity radius compared to case C0, the nozzle decay coefficient decreases by about 6%. Therefore, the submerged nozzle should be carefully considered in the design of a solid rocket motor to avoid combustion instability.

4.3. The Effect of Convergent Angle on Nozzle Damping

The convergent angle of the submerged nozzle can affect nozzle damping, which is related to the reflection of the acoustic waves in the convergent section of the nozzle. At present, it is difficult to predict the effects of nozzle convergent angle on nozzle damping using theoretical analysis and empirical formula, so the numerical method is used in this section. In order to find out the effect of nozzle convergent angle on nozzle damping, four different convergent angles (C0, C4, C5, and C6) within 30°–60° with an increase of 10° are selected for the study, and the other parameters of the nozzle are consistent. According to the data processing method described in Section 4.2, the nozzle decay coefficients under different convergent angles are illustrated in Figure 11.

From Figure 11, it can be seen that the absolute value of the nozzle decay coefficient decreases with the increase of the nozzle convergent angle. On the one hand, the increase of the nozzle convergent angle will lead to the increase of combustion chamber volume, and then the acoustic energy contained in the chamber will increase. On the other hand, with the increase of the nozzle convergent angle, the equivalent reflection effects of the convergent section on the acoustic energy will increase, which decreases the loss of acoustic energy.

In addition, the axial velocity contours under different nozzle convergent angles are illustrated in Figure 12. It can be seen that the gas accelerates in the convergent section and reaches the speed of sound at the nozzle throat, but the axial velocity distributions in the convergent section under four convergent angles are slightly different. A large convergent angle leads to small regions enclosed by the same velocity contour in the convergent section, which further leads to less acoustic energy loss due to wave energy radiation and the convection of wave energy by the mean flow. Therefore, the increase in the convergent angle will result in a decrease in nozzle damping.

4.4. The Effect of Divergent Angle on Nozzle Damping

To examine the influences of divergent angle on the damping capability of the submerged nozzle, four different divergent angles (C0, C7, C8, and C9) are selected for comparison. The divergent angles are within 10°–40° with an increase of 10°, and the decay coefficients are illustrated in Figure 13. As the nozzle divergent angle is increased, it can be seen that the nozzle decay coefficient is slightly decreased and then increased. However, such a variation is trivial. Therefore, no matter for the classical short nozzle or the submerged nozzle, once the flow at the nozzle throat is sonic, the acoustic wave downstream cannot propagate upstream through the nozzle throat, and there is no reflection to occur when the acoustic wave reaches the throat. The incident wave is carried out of the nozzle by the mean flow. The nozzle divergent angle is not the factor that affects nozzle damping.

4.5. The Effect of Gas Temperature on Nozzle Damping

In the practical working process, the gas temperatures inside the combustion chamber of solid rocket motors are very high. However, most current experimental tests are carried out under cold flow conditions. To find out the effect of gas temperature on nozzle damping, the numerical method is used. In this section, different gas temperatures (C0, C10, C11, and C12) ranging from 300 K to 2100 K are studied, representing the transition from cold flow test to hot fire test. In this simulation, for simplicity, it is assumed that the gas temperature in the combustion chamber is distributed uniformly in the hot fire test, which is not true in the real engine.

Figure 14a depicts the decay coefficient under different working temperatures, it can be seen that the predicted nozzle decay coefficient shows a decreasing trend with the increasing temperature. Therefore, there are great differences in nozzle damping between cold flow and hot fire tests. In addition, this part considers two different conditions in the calculation: the gas thermodynamic properties change with temperature and the gas thermodynamic properties remain unchanged. The results show that there are small differences between the decay coefficients obtained under the two conditions. Therefore, the thermodynamic properties are assumed as constant in the following study. Figure 14b depicts the axial velocity of the motor under different working temperatures. It can be seen that the increase in working temperature leads to an increase in the axial velocity in the combustion chamber. Besides, the velocity at the nozzle throat also increases with increasing temperature which is caused by the increase of sound speed. Due to the increase of axial velocity, there would be more acoustic energy flowing out of the nozzle with the mean flow, which results in the increase of nozzle damping.

4.6. The Effect of Frequency on Nozzle Damping

When the chamber is excited by the pulse of a disturbance at the head end, multimode damping due to the presence of several frequencies in the combustion chamber occurs [30,43]. By filtering the pressure oscillation wave, the decay coefficients under certain frequencies can be obtained.

In this section, case C0 is taken as an example to be studied. At first, removing the submerged cavity of the C0 model and other geometric parameters of the model remains unchanged; then it turns into a simple solid rocket motor model. According to the numerical method and the data processing method described above, the time-frequency spectrogram of pressure oscillations and the decay coefficient under different frequencies are illustrated in Figure 15. It can be seen in Figure 15a that the pressure oscillation occurs at 100 ms when the chamber is excited by a pulse and the pressure oscillation is dominated by several axial acoustic modes. After filtering the pressure oscillation wave, the oscillations dominated by the first four axial acoustic modes are obtained and processed separately. The results are illustrated in Figure 15b; it can be seen that the absolute value of the nozzle decay coefficient increases with the increase of frequency. That is to say, the nozzle without a submerged cavity has greater damping effects on the high-frequency acoustic wave, so the axial pressure oscillation with low-frequency usually occurs in the combustion chamber, which is consistent with the phenomenon observed.

Then, the C0 model with the submerged nozzle is considered. The time-frequency spectrogram of the pressure oscillations and the nozzle decay coefficient are illustrated in Figure 16. Apparently, higher-order pressure oscillations decay relatively slowly in Figure 16a. However, in Figure 15a, it can be seen that the power of a high frequency decays faster than that of a low frequency, which represents the greater damping on the high-frequency acoustic wave. In addition, the absolute value of the decay coefficient is decreased with the frequency, which is contrary to that in Figure 15b. Therefore, the nozzle with a submerged cavity has fewer damping effects on the high-frequency acoustic wave.

Based on the above content, it can be deduced that the vortex-acoustic coupling occurring in the submerged cavity plays as a gain factor for the acoustic energy. Therefore, the nozzle damping of the submerged nozzle is equivalent to the sum of a nozzle damping and vortex-acoustic coupling gain, which is weakened by the gain compared to the nozzle without the submerged cavity. With the increase of acoustic frequency, both the damping and gain effects of the nozzle increase, but the gain effect increases more. Thus, the total effect is that the nozzle damping of the submerged nozzle decreases with the increase of frequency, which means that the nozzle with a submerged cavity has fewer damping effects on the high-frequency acoustic wave.

5. Conclusions

In this paper, numerical simulations are carried out to study the nozzle damping characteristics of the submerged nozzle via the pulse decay method. The effects of nozzle geometry and operating conditions on nozzle damping are studied, and the main conclusions are as follows:

• With the increase in cavity size, the acoustic frequency of the motor and the absolute value of nozzle decay coefficient decreases. On average, for every 12.5% increase in the cavity radius compared to case C0, the nozzle decay coefficient decreases by about 6%;
• As the nozzle convergent angle is decreased, the nozzle decay coefficient is increased; in addition, the nozzle divergent angle has a trivial effect on the submerged nozzle decay coefficient;
• There are great differences on nozzle damping between cold flow and hot fire tests, and the nozzle decay coefficient increases with the increase of gas temperature;
• Compared with the nozzle without a submerged cavity, the nozzle with a submerged cavity has fewer damping effects on the high-frequency acoustic wave.