Modeling and Validation of Thrust Prediction of Underwater Solid Rocket Motor

Shilin Hu¹, Chao Yin², Wei Kang^1*, Muyao Xue²

¹School of Astronautics, Northwestern Polytechnical University, Xi’an, China.
²Shanghai Space Propulsion Technology Research Institute, Shanghai, China.

DOI: 10.4236/wjet.2024.124069                

Abstract

The solid rocket motor driven system is one of the common ways for submarines to launch underwater missiles. It has significant advantages in improving the missile’s water exit speed, anti-interference capability, and enemy striking power. The prediction of the underwater loading is a preliminary factor for the power system design of the underwater vehicle. This paper presents a rapid prediction method and validated by the experimental study for the underwater thrust of the solid rocket motor. Based on the potential flow assumption of the water field, a model of the bubble and a one-dimensional quasi-steady model of the nozzle are established to directly solve the flow status of the nozzle. The aerodynamic thrust and hydrodynamic thrust have been calculated and analyzed. The calculation results are within 5% error of the experimental results. Moreover, a design platform to predict the underwater thrust of the solid rocket motor has been developed based on Python and the PyQt library, which shows excellent system adaptability and computational efficiency.

Keywords

Power System Design, Underwater Vehicle, Solid Rocket Motor, Thrust Prediction

1. Introduction

In modern naval warfare, the use of submarines to launch underwater missiles has the characteristics of strong concealment, large maneuvering range and strong offensive capability. Underwater launch is divided into unpowered launch and powered launch according to whether the carrier is powered. For unpowered launch, the missile maintains inertial speed after being ejected from the launch tube, and the trajectory must be stable. Since the cavitation formed during the launch process imposes a great influence on the ballistic stability, the general unpowered launch exhibits a low exit speed and poor anti-interference ability, and cannot adapt well to the variant combat environment in high-intensity naval battles. In order to improve the ability to raid and strike the enemy, powered launch with solid rocket motor (SRM) have gradually been developed and applied. Compared with unpowered launch, powered launch has significant advantages in exit speed, attitude stability, launch depth, anti-interference ability and enemy strike capability.

Throughout the world, nations have begun to equip with submarine-launched missiles with underwater ignition capabilities, such as Russia’s SS-N-21 submarine-launched missile, France’s SM-39 Exocet short-range submarine-to-ship missile, and the United States’ Tomahawk missile launched using torpedo tubes. In order to accurately obtain the trajectory of the missile during a powered launch and maintain its stability, the thrust and its change process during underwater ignition are particularly important. During underwater ignition, the density of water is much greater than that of the working gas, and the high-temperature combustion gas will interfere with the surrounding water field after being ejected from the tail nozzle. These problems will cause a drastic change in the surrounding pressure field. How to accurately and quickly predict the thrust change for the underwater weapon has become the key to the powered launch of the missile.

In the underwater powered system, the coupling between the water field and the gas bubble field is the source of the thrust oscillation of underwater SRM. How to model and predict this coupling relationship is the key to predicting the thrust. The prediction of oscillating thrust is based on the comprehensive study of the growth and shedding process of the bubbles. However, because this method is of complexity and has low computational efficiency, it cannot meet the rapid response requirements of the modern combat environment. Therefore, it is urgent to develop a rapid prediction method for underwater thrust.

Theoretical and experimental studies have been carried out on the underwater powered system. Lu et al. [1] used the time-stepping method to achieve the numerical solution of water-gas-bubble flow for the axisymmetric underwater loads problem. The gas in the nozzle was solved according to the one-dimensional quasi-steady model using the finite difference method, and the water field was solved using the boundary integral equation. The bubble was simplified to an isentropic bubble model, and the evolution of the bubble was simulated using the Euler-Lagrange method. The calculation results showed that the pressure in the bubble will quickly form a peak, and then the pressure will gradually become flat. Huang et al. [2] studied the underwater thrust prediction at different water depths from the perspective of fluid dynamics. They divided the solution system into three parts: water field solution, nozzle flow solution, and bubble solution. The water field solution used an incompressible potential flow model, the nozzle flow field solution used an unsteady solution model, and the bubble solution used a zero-dimensional calculation model based on mass and energy conservation. The research results showed that the water depth has a great influence on the pressure of the tail bubble. When the water depth is large, the pressure inside the bubble will oscillate, which will cause the hydrodynamic force to oscillate, and then have a great impact on the missile thrust. Wang et al. [3] used the same solution method for the underwater thrust prediction. They considered the influence of heat transfer and vaporization between high-temperature gas and water medium. The calculation results simulated the growth and shedding process of the bubble and revealed the basic law of flow inside the bubble. Shan et al. [4] unified the nozzle solution and the bubble solution based on the axisymmetric Euler equation. The potential flow theory was used for the water flow field. By solving the three-dimensional Hess-Smith equation, the development process of the bubble was obtained. Yi et al. [5] used the finite difference method to solve the N-S equation to obtain the water field parameter distribution near the bubble. The calculation results showed that a pressure peak is formed in the bubble in a short period of time.

During the underwater ignition of the engine, the coupling between the water field and the gas bubble field is the source of the underwater engine thrust oscilla-tion. How to model and predict this coupling relationship is the key to predicting the engine thrust [6] [7]. The paper adopts the assumption of potential flow in the water field, then the nozzle outlet status and the thrust of the underwater SRM are directly solved. The thrust oscillation is captured and compared with the experimental results. Finally, a underwater SRM thrust prediction platform for various working conditions is developed. The paper is organized as follows: Section 1 gives a review of the prediction of underwater SRM thrust. The mathematical models are presented in Section 2. The calculation results and the comparison with the experimental results are shown in Section 3. The main conclusions are given in Section 4.

2. Mathematical Models

2.1. Governing Equations for Bubble

The water is assumed incompressible and inviscid. The velocity potential of the water can be solved by the Laplace’s equation.

${\nabla }^2\phi =0$ (1)

The boundary condition at nozzle outlet and infinity:

$\{\begin{array}{l}\frac{\partial \phi }{\partial n}=0,\text{nozzle}\text{outlet}\\ \frac{\partial \phi }{\partial n}=0,\phi =0,\text{infinity}\end{array}$ (2)

By substitute water velocity by water velocity potential in the Bernoulli equation, the water velocity potential can be solved by:

$-gH+\frac{p_w}{{\rho }_w}+\frac{1}{2}{\left(\nabla \phi \right)}^2+\frac{\partial \phi }{\partial t}=\frac{p_l}{{\rho }_w}$ (3)

where $p_w,{\rho }_w$ is the static pressure and the density of water respectively. $p_l$ is the atmosphere pressure at sea level. The Equation (3) can be discretized:

$\begin{array}{l}\frac{\phi \left(t+\Delta t,r\left(t+\Delta t\right),z\left(t+\Delta t\right)\right)-\phi \left(t,r\left(t+\Delta t\right),z\left(t+\Delta t\right)\right)}{\Delta t}\\ =\frac{p_l+{\rho }_{w}gH-p_w}{{\rho }_w}-\frac{1}{2}{[}^{{(}^{\frac{\partial \phi }{\partial n}}}]\end{array}$ (4)

$\begin{array}{l}\phi \left(t,r\left(t+\Delta t\right),z\left(t+\Delta t\right)\right)\\ =\phi \left(t,r\left(t\right),z\left(t\right)\right)+\left\{\frac{p_l+{\rho }_{w}gH-p_w}{{\rho }_w}+\frac{1}{2}\left[{\left(\frac{\partial \phi }{\partial n}\right)}^2-{\left(\frac{\partial \phi }{\partial s}\right)}^2\right]\right\}\Delta t\end{array}$ (5)

Then the water potential $\phi =\phi \left(t_0,r\left(t_0\right),z\left(t_0\right)\right)$ at the interface between the bubble and the water at time $t_0$ can be obtained. The $\frac{\partial \phi }{\partial s}$ at time $t_0$ can also be solved readily( $\frac{\partial \phi }{\partial s}=\frac{{\phi }_{i+1}-\phi }{\Delta L}$ , $\Delta L$ is the distance between the point $i$ and $i+1$ ). Moreover, the value of $\frac{\partial \phi }{\partial n}$ is calculated by the boundary element methods.

The location of the bubble is controlled by:

$\{\begin{array}{l}\frac{dr}{dt}=\frac{\partial \phi }{\partial n}\sin \theta +\frac{\partial \phi }{\partial s}\cos \theta \\ \frac{dz}{dt}=\frac{\partial \phi }{\partial n}\cos \theta -\frac{\partial \phi }{\partial s}\sin \theta \end{array}$ (6)

where $\theta$ is the angle between the external normal direction and the positive z-axis direction. The $\frac{\partial \phi }{\partial n}$ at time $t_0$ can be solved by the boundary integral equation:

$C\phi ={\displaystyle {\int }_{\Gamma }\left[{\phi }^{*}\cdot \frac{\partial \phi }{\partial n}-\phi \cdot \frac{\partial {\phi }^{*}}{\partial n}\right]d\Gamma }$ (7)

where $\text{Γ}$ represents the interface between bubble and water. The fundamental solution of Laplace operator is represented by ${\phi }^{*}$ , which can be expressed by:

${\phi }^{*}=\frac{1}{2\pi }ln\frac{1}{r\left(P,Q\right)}$ (8)

where $r\left(P,Q\right)$ is the distance between point $P$ and the source point $Q$ . The value of $C$ is decided by the location of the point $P$ .

$\{\begin{array}{l}C=1,P\in \Omega \\ C=\frac{\theta \left(P\right)}{2\pi },P\in \Gamma \end{array}$ (9)

where $\theta$ is the solid angle of boundary point $P$ in three dimensions. The velocity potential $\phi$ and the related normal gradient $\frac{\partial \phi }{\partial n}$ of water can be solved by the boundary value.

2.2. Governing Equations for Nozzle

2.2.1. Nozzle Flow Status Model

After obtaining the governing equations of the bubble, it is also necessary to solve the mass flow at nozzle outlet $\stackrel{\dot{}}{M}$ in advance. In this paper, the flow is steady and quasi-one-dimensional throughout the nozzle.

Firstly, the critical back pressure of the nozzle is determined according to the nozzle profile size (nozzle outlet area $A_{out}$ , nozzle throat area $A_t$ ) and the total pressure parameters of the combustion chamber. The flow aerodynamic function $q\left({\lambda }_{out}\right)$ of the nozzle outlet of this configuration is obtained by the ratio of the throat cross-sectional area to the outlet cross-sectional area:

$q\left({\lambda }_{out}\right)=\frac{A_t}{A_{out}}={\left(\frac{\gamma +1}{2}\right)}^{\frac{1}{\gamma -1}}{\lambda }_{out}{\left(1-\frac{\gamma -1}{\gamma +1}{\lambda }_{out}{}^2\right)}^{\frac{1}{\gamma -1}}$ (10)

where ${\lambda }_{out}$ is the nozzle exit velocity coefficient, which has two solutions. Then the values of the velocity coefficient at the exit position in the supersonic part and the subsonic part ${\lambda }_{out,\text{super}}$ and ${\lambda }_{out,\text{sub}}$ can be obtained.

According to the one-dimensional isentropic flow theory, the outlet aerodynamic function $\pi \left({\lambda }_{out,\text{super}}\right)$ is obtained. The first critical back pressure $p_{a1}$ can be obtained:

$\pi \left({\lambda }_{out,\text{super}}\right)={\left(1-\frac{\gamma -1}{\gamma +1}{\lambda }_{out,\text{super}}{}^2\right)}^{\frac{\gamma }{\gamma -1}}$

$p_{a1}=\pi \left({\lambda }_{out,\text{super}}\right)\ast p_0$ (11)

where $p_0$ is the total pressure of the combustion chamber. Based on the subsonic part in the one-dimensional isentropic flow aerodynamic function table, the outlet aerodynamic function $\pi \left({\lambda }_{out,\text{sub}}\right)$ is obtained. The third critical back pressure $p_{a3}$ is obtained:

$\pi \left({\lambda }_{out,\text{sub}}\right)={\left(1-\frac{\gamma -1}{\gamma +1}{\lambda }_{out,\text{sub}}{}^2\right)}^{\frac{\gamma }{\gamma -1}}$

$p_{a3}=\pi \left({\lambda }_{out,\text{sub}}\right)\ast p_0$ (12)

The second critical back pressure $p_{a2}$ is obtained by calculating the flow status of the normal shock wave at the outlet position:

$p_{a2}=p_{a1}\times \frac{p_{a2}}{p_{a1}}$ (13)

The flow status inside the nozzle is determined by the ambient pressure $p_b$ :

・ When $0<p_b<p_{a1}$ , the nozzle is under-expanded.

・ When $p_b=p_{a1}$ , the nozzle is full-expanded.

・ When $p_{a1}<p_b<p_{a2}$ , the nozzle is over-expanded, the shock is outside the nozzle.

・ When $p_b=p_{a2}$ , the nozzle is over-expanded, the shock is at the nozzle exit.

・ When $p_{a2}<p_b<p_{a3}$ , the nozzle is over-expanded, the shock is inside the nozzle.

・ When $p_b=p_{a3}$ , the nozzle throat reaches the speed of sound, and the rest of the flow is subsonic.

・ When $p_b>p_{a3}$ , the flow inside the nozzle is completely subsonic.

Due to the relationship between the exit velocity coefficient and the Mach number, the exit Mach number $M{a}_{out}$ is obtained. The exit temperature $T_{out}$ and the exit pressure $p_{out}$ can also be obtained according to the isentropic relationship:

$M{a}_{out}=\sqrt{\frac{\frac{2}{\gamma +1}{\lambda }_{out}{}^2}{1-\frac{\gamma -1}{\gamma +1}{\lambda }_{out}{}^2}}$

$T_{out}=T_0\cdot \left(1-\frac{\gamma -1}{\gamma +1}{\lambda }_{out}{}^2\right)$ (14)

$p_{out}=p_0\cdot {\left(1-\frac{\gamma -1}{\gamma +1}{\lambda }_{out}{}^2\right)}^{\frac{\gamma }{\gamma -1}}$

$V_{out}=\sqrt{\gamma R_g{T}_{out}}\cdot M{a}_{out}$

When $p_b<p_{a3}$ , the nozzle throat reaches the speed of sound, and the mass flow at outlet is directly determined by the combustion chamber state parameters:

$\stackrel{\dot{}}{M}={\left(\frac{\gamma }{R_g}\right)}^{1/2}{\left(\frac{2}{\gamma +1}\right)}^{\frac{\gamma +1}{2\left(\gamma -1\right)}}\frac{p_0}{\sqrt{T_0}}{A}_t$ (15)

where $T_0$ is total temperature of the combustion chamber.

When $p_b>p_{a3}$ , the flow inside the nozzle is completely subsonic, the mass flow at outlet:

$\stackrel{\dot{}}{M}={\rho }_{out}{V}_{out}{A}_{out},{\rho }_{out}=\frac{p_{out}}{R_g{T}_{out}}$ (16)

2.2.2. Thrust Prediction Model of Underwater SRM

The force situation at the nozzle exit is shown in Figure 1. The thrust of underwater SRM is divided into two parts. The first part is the aerodynamic thrust acting on the bottom of the nozzle, which can be obtained by solving the nozzle outlet flow state

Figure 1. The schematic of underwater SRM thrust.

$F_{aero}=\stackrel{\dot{}}{M}\cdot V_{out}+A_{out}\cdot \left(p_{out}-p_{b0}\right)$ (17)

where $p_{b0}$ is the static pressure of the water. The second part is the hydrodynamic thrust acting on the bottom of the underwater SRM, which can be solved by the internal pressure of the bubble:

$F_{hydro}=A_{hydro}\cdot \left(p_b-p_{b0}\right)$ (18)

where $A_{hydro}$ is the bottom area of the underwater SRM.

The thrust of the underwater SRM:

$F=\stackrel{\dot{}}{M}\cdot V_{out}+A_{out}\cdot \left(p_{out}-p_b\right)+A_{hydro}\cdot \left(p_b-p_{b0}\right)$ (19)

3. Results and Discussion

The calculation of underwater thrust depends on the operating status of the combustion chamber. The variation in the combustion chamber pressure can be divided into several operating phases: 1) The whole operating process that includes the initial phase with a sharp rise, the stable working phase with a relatively flat pressure output, and the backfire phase with a sharp drop in pressure; 2) The stable state that directly links to the design condition, which accounts for the major operating time. Therefore, it can only consider the stable state when predicting the underwater thrust of the SRM for the design purpose. For situations that require consideration of the entire working process of the combustion chamber, the ballistic curve should be considered as an input for the underwater thrust prediction. This section will discuss the underwater thrust under the two mentioned states.

3.1. Thrust Prediction of Underwater SRM in the Stable State

The case of underwater SRM thrust calculation from [8] is selected. The related parameters of the SRM are shown in Table 1. The initial radius of the bubble is set to 0.06m, and the initial bubble pressure is set to 0.5Mpa. The calculation process lasts for 0.1s.

Table 1. The parameters of the underwater SRM.

Figure 2 and Figure 3 present the variations of the bubble pressure and the thrust over time. It can be observed that the calculated bubble pressure and thrust are in good agreement with the literature results. In the initial stage, the bubble pressure exhibited a peak that is six times higher than the local hydrostatic pressure at the water depth of 30 m, which has been confirmed in many theoretical and experimental studies. At the same time, a significant thrust peak is also observed. In the initial stage, the bubble pressure $p_b$ is greater than the local hydrostatic pressure $p_{b0}$ due to presence of the bubble. Therefore the hydrodynamic thrust term $A_{hydro}\cdot \left(p_b-p_{b0}\right)$ in Equation (19) is relatively large, causing the thrust peak. As the characteristic radius of the bubble increases, the bubble

Figure 2. The variation of the bubble pressure.

Figure 3. The variation of the thrust.

pressure drops rapidly and gradually fluctuates around the value of local hydrostatic pressure. The change in bubble pressure leads to a gradual decrease in the hydrodynamic thrust term, resulting in the oscillation of the thrust.

3.2. Thrust Prediction of Underwater SRM with Ballistic Curve

There are initial, steady, and backfire phase in actual combustion chamber operating conditions. The total pressure in the chamber is no longer a constant value, and the working state of the nozzle will change over time. For the governing equation of bubble in Section 2.1, the exit mass flow rate $\stackrel{\dot{}}{M}$ will depend on the flow state of the nozzle.

The nozzle profiles used in the experiment are shown in Table 2. The Adiabatic index of the gas is 1.16. The ballistic curves of the nozzles are shown in Figure 4.

Table 2. The parameters of nozzles.

Figure 4. The ballistic curves of the nozzles.

Four different models of nozzles are tested at water depths of 25 m, 50 m, 75 m, and 100 m, respectively. Underwater thrust prediction calculations and error analysis are conducted.

Figure 5 presents the thrust curve of No. 01 nozzle and its comparison with the experimental results. At the initial stage, the nozzle is in a subsonic working state since a large amount of gas in the nozzle is obstructed and cannot be expelled smoothly due to the inertia of water. As the total pressure rises rapidly, the obstructed gas is fully expanded. Therefore, the nozzle quickly transits from a subsonic state to a supersonic state. Figures 6-8 present the thrust curves of No. 02,

Figure 5. The thrust comparison with experimental results of No. 01 nozzle.

Figure 6. The thrust comparison with experimental results of No. 02 nozzle.

Figure 7. The thrust comparison with experimental results of No. 03 nozzle

Figure 8. The thrust comparison with experimental results of No. 04 nozzle.

No. 03, and No. 04 nozzles and their comparisons with experimental results, respectively. The trends are similar to those of No. 01 nozzle.

Figures 9-12 show the statistical error distribution of sample points for four types of nozzles. More than 70% of the sample points have an error within 10%. Only a few sample points have significant errors, such as the initial and tailing stages of the ballistic curves. The minimum average error of underwater thrust for the four types of nozzles is 0.50%, and the maximum average error is 3.53%, which demonstrates the reliability and accuracy of the calculation model.

Figure 9. The statistical error distribution of No. 01 nozzle.

Figure 10. The statistical error distribution of No. 02 nozzle.

3.3. Prediction Platform of the Underwater SRM Thrust

To meet the needs of multi-input-condition calculations for thrust of the underwater SRM, a platform is developed based on Python and the PyQt library. The interface is shown in Figure 13. The results can be illustrated in the platform as shown in Figure 14. This platform has the following functions and features:

・ Two working state of the combustion chamber can be choosed: calculation in stable state & calculation with ballistic curve.

Figure 11. The statistical error distribution of No. 03 nozzle.

Figure 12. The statistical error distribution of No. 04 nozzle.

・ The water depth can be set freely with a range of 0-300 meters. The simulation platform supports two calculation modes: fixed depth calculation & variable depth calculation.

・ The nozzle parameters and physical parameters can be freely set and changed.

・ The simulation platform supports automatically calculation and results output in multi-condition calculations.

・ All operations are interactively executed in a visual interface.

Figure 13. The prediction platform of the underwater SRM thrust.

Figure 14. Results illustration in the platform.

4. Conclusions

In this study, Bernoulli equation with velocity potential is established for the water filed with the assumption of incompressible and inviscid water. The nozzle flow status of the underwater SRM is solved based on one-dimensional flow theory. The aerodynamic thrust and hydrodynamic thrust acting on the underwater SRM are calculated respectively. A predictive model for the thrust of the underwater SRM is established. The following conclusions have been drawn:

・ The governing equation for the bubble and the coupling relationship between water field and bubble have been established and calculated.

・ The governing equation for the nozzle and the predictive model for the thrust of the underwater SRM is established and calculated.

・ Thrust predictions with ballistic curves and stable states with various water depths have been conducted. The computational results indicate that the model accurately captured the thrust peak and oscillation phenomena. The error between the computational and experimental results is within 5%.

・ A simulation platform for underwater SRM thrust has been developed based on Python and PyQt library. The platform has excellent adaptability to various working environments and can perform the automation of multiple working cases calculation efficiently.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Lu, C.-J., Chen, F., Fan, H., et al. (1992) The Fluid Dynamic Research on the Under-Water Ignition of Missile. Acta Aeronautica et Astronautica Sinica, No. 4, 124-130.
[2]	Hang, J.C., Ye, Q.Y. and Zhu, S.Q. (1994) Gas-Water Dynamic Calculation for the Underwater Ignition of a Missile at Different Depths. Chinese Journal of Applied Mechanics, No. 3, 19-24+138.
[3]	Wang, C., Ye, Q.Y. and He, Y.S. (1997) Calculation of an Exhausted Gas Cavity behind an Under-Water Launched Missile. Chinese Journal of Applied Mechanics, No. 3, 3-9+147.
[4]	Shan, X.-X., Yang, R.-G. and Ye, Q.-Y. (2001) Fluid Forces on a Missile with Control System of Vectorial Thrust. Journal of Shanghai Jiaotong University, No. 4, 625-629.
[5]	Yi, S.-Q., Yan, K., Xu, S., et al. (2002) Numerical Simulation of the Initial Stage of Noncondensing High-Speed Gas Jet in Liquid. Journal of Hydrodynamics, No. 4, 448-453.
[6]	Fu, Y., Wei, Y. and Zhang, J. (2009) Parametric Study on the Thrust of Bubbly Water Ramjet with a Converging-Diverging Nozzle. Journal of Hydrodynamics, 21, 591-599. https://doi.org/10.1016/s1001-6058(08)60189-4
[7]	Song, A., Wang, N., Li, J., Ma, B. and Chen, X. (2020) Transient Flow Characteristics and Performance of a Solid Rocket Motor with a Pintle Valve. Chinese Journal of Aeronautics, 33, 3189-3205. https://doi.org/10.1016/j.cja.2020.04.023
[8]	Zhang, Y.-W. (2007) Research of Working Characteristic of Solid Motor in Water. University of Science and Technology of China.