Abstract
This study entails the analysis of the working performance of solid rocket motors (SRM), featuring the essential element of internal ballistic analysis. Therefore, the internal flow field under the condition of burning surface regression needs to be calculated. The boundary of the internal flow field of the SRM moves with the combustion of the propellant; therefore, it is necessary to accurately track the mobile interface to provide boundary conditions for the flow field calculation. The coupling of the level set method and the volume fraction method is utilized to track the burning surface, and the porous media model is used to divide the fluid and solid calculation domains. The interface between the two calculation domains is used to characterize the burning surface, and then, the area of the burning surface is obtained by solving the area of the interface. The calculation and analysis are carried out for SRM with tubular charge and star charge. The results verify that the calculation model can accurately calculate the transient internal flow field of SRM under the condition of burning surface regression.
1. Introduction
Solid rocket motors (SRM) are widely used as power plants on missiles and spacecraft because of their superior performance. In the research and development of SRM, a series of experiments need to be carried out to predict performance, which often consume huge funds. With the development of computer technology, numerical calculation methods to simulate the working process of SRM can be used to supplement the experimental means, and they play a highly important role in reducing expenses and shortening the development cycle. Therefore, directions involving such methods have become one of the most active ones in the field of SRM. The main purpose of computer simulation calculation is to obtain the internal ballistic curve of the SRM, which is the law of working pressure in the combustion chamber changing with time and space. Moreover, pressure is a significant parameter in SRM operation, which determines a series of performance parameters such as thrust, working time, and structural integrity of the SRM; hence, it is necessary to accurately calculate the interior ballistic curve. The burning surface area of the grain is a decisive factor for this calculation; thus, it is necessary to use the corresponding algorithm to simulate the burning surface regression, so as to obtain the burning surface area change curve of the grain during the combustion of the propellant [1, 2].
For the numerical simulation of burning surface regression, accurate mathematical models and robust and effective numerical techniques are required to ensure that the burning surface regression is synchronized with the combustion of the propellant. Currently, there are three main methods for calculating the burning surface regression, including the following: (1) The solid modeling method mainly relies on computer-aided design (CAD) software for secondary development, which can obtain the area of the burning surface at any time; however, it is difficult to realize the coupling calculation with the flow field [3–5]. (2) The dynamic mesh method continuously updates the mesh through a certain algorithm to simulate the movement of the interface. This method easily realizes the coupling calculation with the flow field, but the dynamic interface needs to be explicitly defined in the calculation process. Therefore, when the geometry of the grain is fairly complex, this method will become very difficult [6–8]. (3) The interface tracking method regards the burning surface of the grain as the boundary of the flow field and tracks the burning surface through a certain algorithm [9–12]. Among the relevant methods, the widely used ones are the level set (LS) method [13] and the volume of fluid (VOF) method [14], each with its inherent advantages and limitations. The LS method can compute the calculation model of complex interface deformation or topological structure changes, and more importantly, regardless of the changes of flow field, the LS function always remains smooth. Nonetheless, this method cannot satisfy the conservation of mass during reinitialization; thus, it will cause the loss of physical quantity. The VOF method can ensure the conservation of physical quantities, but it is difficult to ensure accuracy when calculating curvature, normal vector, and geometric parameters, which will cause the interface to be easily damaged and difficult to reconstruct and cause calculation errors. Therefore, the VOF method will be very difficult when dealing with complex-shaped interfaces. Consequently, a new interface tracking coupling method was proposed in recent years, which is called the CLSVOF method [15–17]. This approach not only overcomes the shortcomings of the VOF method in that it is difficult to accurately calculate the normal and curvature of the interface but also solves the problem that LS is not a conservative method, i.e., addresses the loss of physical quantities in the calculation process.
The key problem for the coupled calculation of burning surface regression and transient flow field is to establish mathematical models and numerical techniques. The coupling process is realized through the data transmission between the models, and the algorithm needs to be optimized in the calculation, so that the accuracy, efficiency, and stability of the simulation calculation are realized at a lower calculation cost. Accordingly, the outline of the paper is as follows. In Section 2, firstly, we introduce the coupling process of the level set method and the volume of fluid method (Section 2.1). Secondly, we describe the numerical modeling of the flow field calculation, which mainly includes the governing equations, turbulence models, and the numerical discrete (Section 2.2). Thirdly, we detail the mathematical models making up the porous media model (Section 2.3). Finally, we elaborate on the coupling method of burning surface regression and flow field calculation (Section 2.4). In Section 3, we carry out numerical simulations of the burning surface regression of SRM with tubular charge (Section 3.1) and SRM with star charge (Section 3.2), then compare and analyze the numerical calculation results with the theoretical calculation results. The paper ends with conclusions in Section 4.
2. Numerical Modeling
2.1. Coupling of Level Set and Volume of Fluid Formulation
The VOF method tracks the interface by defining the volume fraction . The volume fraction is defined as the ratio of the volume of fluid in the cell to the volume of the cell and is solved by the VOF equation. The governing equation is as follows [18]: where if , the mesh is occupied by solid; if , the mesh is occupied by fluid; and if , there is a moving interface in the mesh.
When solving this volume fraction transport equation, in order to ensure the sharpness of its interface, a compression term is artificially added into Equation (1). The revised equation is as follows: where is solved by the method discussed below. The compression term works at the interface only when is not equal to 0; the compressibility coefficient is usually .
The LS method is a distance function method, where the two-phase interface is represented by the zero point of a high-order LS function , which uses algebraic values to distinguish the phases in the calculation area. The LS method can be expressed as a continuous function:
The governing equation of is
To ensure the sharpness of the volume fraction value obtained by the VOF algorithm and the boundedness of the LS distance function, the obtained volume fraction value is coupled to the level set distance function. Here, we use an explicit algorithm to limit the flux to solve Equation (2). The following is the implementation process of the CLSVOF method: where , and represents the mesh cell size.
When the LS distance function is expressed by the value, in order to ensure the accuracy of the value, if the volume fraction in a mesh is less than 0 or greater than 1 in a time step, the flux of needs to be modified to limit the volume fraction. Therefore, when solving Equation (1), several subcycles will be performed in a single step. The reinitialization of the LS function is then performed, so that the solution of LS function is always a signed distance function.
Suppose that is the numerical solution of the LS function at a certain moment, and is expressed as the interface. Construct a function as the signed distance function at this moment and make still represent the interface at this moment. Perform the iteration of the following equation:
Equation (6) is solved until the steady state, and is used as the new signed distance function. The interface normal vector can then be calculated through the function. The implementation process of the CLSVOF method is shown in Figure 1:
The main steps of CLSVOF method calculation are as follows: (1)Calculate the volume fraction of the current step (step ) and the velocity information of the flow field(2)Map the volume fraction to the initial distance function by solving Equation (6)(3)Obtain the final distance function obtained by reinitializing the algorithm to solve Equation (7)(4)Solve the interface position by the obtained distance function , then solve Equation (5) to get the volume fraction for the next step(5)After the distance function and the volume fraction have been solved, go to step 1 and start the next iteration
2.2. Flow Analysis
2.2.1. Evaluation and Simplification of Physical Problems
The combustion of the propellant causes burning surface regression, which is a complex process. Consequently, it is difficult to calculate the coupling of burning surface regression and transient internal flow field under propellant burning conditions, and thus problematic, it is difficult to realize the coupling calculation of the burning surface regression and the transient flow field of the combustion chamber under propellant burning conditions. In contrast, it is easier to simulate the transient flow field of the combustion chamber under the burning surface regression by using the burning surface regression calculation method instead of considering the burning mechanism. In order to simplify the calculation, the following assumptions are made: (1)There is uniformity for the physical and chemical composition of the propellant(2)The entire burning surface of the grains burns at the same time(3)The burning surface advances toward the inside of the grain along the normal direction at the burning rate of the local propellant(4)The gas produced by the combustion of propellant is an ideal gas with a single component(5)Erosive burning is not included in the calculation(6)The regression velocity of the burning surface is calculated through St. Robert’s Law , where represents the pressure of the combustion chamber, is the burning rate coefficient, and denotes the burning rate pressure exponent(7)The solid properties of the grain are characterized by limiting the flow velocity of the solid domain to zero through the porous media model
2.2.2. Governing Equations
Gravity and radiation heat transfer in the combustion chamber, as well as the chemical reaction process of propellant combustion, are ignored in the calculation. When the burning surface regression reaches a certain place, the propellant in that place is burnt and is directly converted into gas; that is, the mass source and energy source are added in the mesh to simulate the injection of gas. The flow governing equations are as follows [19]: where represents the burning surface area inside the mesh, as shown in Figure 2, and denotes the volume of the mesh. is the total temperature, represents the reference temperature, and K. represents the viscous resistance coefficient, and represents the inertial resistance coefficient.
2.2.3. Turbulence Model
The flow in the combustion chamber and nozzle of the SRM is turbulent, which comprises the following aspects: shock wave, boundary layer, shear layer, recirculation zone, and the interaction of all these. Taking into account the complexity of the flow field and the compressibility of the gas, in this paper, we use the SST model to calculate turbulence, which has shown good performance in calculating the distribution of recirculation zone, wall pressure distribution, temperature distribution, and velocity distribution. The equations of the SST model are as follows [20]: where ; represents the production of turbulence kinetic energy; represents the generation of ; and , respectively, represent the effective diffusivity of and ; and , respectively, represent the dissipation of and due to turbulence; represents the cross-diffusion term; and and are user-defined source terms.
The eddy viscosity can be calculated by the following formula: where indicates the strain rate magnitude and is the distance to the next surface. The coefficient damps the turbulent viscosity causing a low-Reynolds number correction and is given by where in the high-Reynolds number case, .
2.2.4. Numerical Discrete
For the spatial discretization, a structured multiblock finite volume approach with cell-centered data storage is adopted. The diffusion term is discretized by the central difference scheme, and the convection term is discretized by the second-order upwind scheme. For the time discretization, the fourth-order explicit Runge-Kutta method is adopted. The time step is chosen as the minimum value of maximum stable time step in all fluid cells, and the SIMPLE algorithm is applied for the calculation of the flow field.
2.3. Porous Media Model
Porous media are modeled by the addition of a momentum source term to the standard fluid flow equations [21]. The source term is composed of two parts: a viscous loss term (the first term on the right-hand side of Equation (14) and an inertial loss term (the second term on the right-hand side of Equation (14)):
The viscous resistance coefficient and inertial resistance coefficient are added through the porous media model to limit the flow velocity of the solid domain to zero. When the differential pressure between the solid domain and the fluid domain is large, the slow flow in the solid domain is allowed to balance the differential pressure, so that the calculation will not diverge because of the pressure jump. In order for the equation to meet the calculation requirements of the fluid domain as well as the solid domain, the viscous resistance coefficient and inertial resistance coefficient need to meet the terms below. And Figure 3 demonstrates the distribution of and in different computational domains. (1)For the solid domain, the values of and are set to very large numbers (such as 1010), so that the flow of the solid domain is blocked(2)For the fluid domain, the values of and are set to zero to ensure that there is no effect on the flow(3)The coefficients at the interface can be well transitioned, so that the calculation efficiency is not affected
The volume fraction is used to construct a function to calculate and at the interface, and two parameters and are introduced to make the function meet the above three terms, where and . The relevant formula is as follows:
has the same value as , and their distribution curve is shown in Figure 4.
2.4. Coupling of Burning Surface Regression and Flow Field Calculation
The combination of CLSVOF and porous media modeling technology is used to perform the coupled calculation of burning surface regression and SRM transient internal flow field. The CLSVOF is used to control the regression of the combustion surface and to provide the area of the burning surface and volume fraction for the porous media model for the calculation of the source term to simulate the generation of gas. Computational fluid dynamics are used to calculate the gas flow in the combustion chamber to obtain the pressure distribution in the flow field, so that the burning rate equation can be used to compute the burning surface regression velocity. The coupled calculation scheme is shown in the flow chart (Figure 5).
3. Numerical Tests and Discussion
3.1. Tubular Charge SRM Transient Internal Flow Field Simulations
3.1.1. Physical Model
Figure 6 presents the SRM model with tubular charge. The propellant parameters are shown in Table 1. The charge of this motor is tubular, the front-end surface is covered, and the inner surface and the rear-end surface are burning surfaces.
3.1.2. Calculation Results and Analysis
The length-diameter ratio of the tubular charge SRM is relatively small, the gas flow rate in the combustion chamber channel is low, and there is no obvious gap along the axis. Therefore, it is appropriate to use the zero-dimensional internal ballistic method to calculate the combustion chamber pressure. The results calculated by the zero-dimensional interior ballistic method and the CLSVOF method are then compared to verify the accuracy of the CLSVOF method. The zero-dimensional interior ballistic calculation process is as follows: (1)Suppose the thickness of the charge is and divide it into segments, where the thickness of each section will be (2)When the burning thickness of the charge is , the surface area of the burning surface at this time is expressed as ()(3)The pressure under each combustion layer is calculated by the formula , where is the characteristic velocity of charge and denotes nozzle throat area(4)The burning rate at step is
Figure 7 shows the contours of pressure distribution of the flow field at different time points, and the clear flow field boundary can be seen from the figure. The pressure of the combustion chamber at different times is extracted and compared with the calculation results of the zero-dimensional internal ballistic method, as shown in Figure 8. The pressure curve can be divided into three segments: ascending segment, stable segment, and descending segment. The calculation results of the two methods of ascending segment and descending segment are consistent, while the result of the CLSVOF method in the stable segment is greater than the result of the zero-dimensional interior ballistic method. Based on the analysis, this might be due to the following two reasons: (1) At the beginning of SRM operation, the stable flow in the combustion chamber is not established, the mass flow rate of gas injected into the combustion chamber from the burning surface is greater than the rate of mass flow from the nozzle, and the zero-dimensional interior ballistic method does not involve this effect in the calculation. Thus, the combustion chamber pressure calculated by the CLSVOF method is higher. (2) The CLSVOF method considers the actual three-dimensional flow in the calculation; therefore, the calculated size of the nozzle throat area is smaller than the geometric size, and this is also consistent with the actual working process of the SRM.
As presented in Figure 9, the position of the burning surface at different times is extracted, such that the entire process of the regression of the burning surface of the tubular charge can be seen. Next, the area of the burning surface is extracted during the regression of the burning surface and compared with the result calculated by the computer-aided design (CAD) method, as demonstrated in Figure 10, showing that the calculation results of the two methods are consistent.
(a)
(b)
(c)
(d)
(e)
In summary, the effectiveness and accuracy of the CLSVOF method were verified through the comparative analysis results of the combustion chamber pressure and the combustion surface area. Thus, this method can effectively calculate the transient internal flow field of the tubular charge SRM under the regression of the burning surface.
3.2. Star Charge SRM Transient Internal Flow Field Simulations
3.2.1. Physical Model
The calculation and analysis of star charge SRM are carried out to further verify the ability of the CLSVOF method to calculate the complex burning surface. Figure 11 shows the cross-section of the star charge, and the specific parameters are as follows: outside diameter mm, charge length mm, characteristic length mm, star tip arc radius mm, star root arc radius mm, convergence half angle of nozzle , expansion half angle of nozzle , and nozzle area expansion ratio . The propellant parameters are shown in Table 2.
3.2.2. Calculation Results and Analysis
The area of burning surface is extracted during the regression of the burning surface and compared with the result calculated by the CAD method, as seen in Figure 12, which shows that the calculation results of the two methods are consistent. The position of the burning surface at different times is extracted, as shown in Figure 13, from thus illustrating which the entire process of the regression of the burning surface of the star charge.
(a)
(b)
(c)
(d)
(e)
Figure 14 shows the contour of the pressure distribution of the flow field in the star charge SRM, Figure 15 shows the contour of the temperature distribution, Figure 16 shows the contours of the Mach number distribution, and Figure 17 shows the contour of the velocity distribution.
The numerical calculation results and the theoretical calculation results are compared and analyzed to verify the accuracy of the numerical calculation results earlier. The theoretical calculation formula of the parameters at the nozzle outlet is as follows.
The relationship between the nozzle area expansion ratio and Mach number is defined as
The relationship between pressure and Mach number :
The relationship between temperature and Mach number :
The formula for the speed of sound : where the specific heat ratio is , the total pressure is MPa, the total temperature is K, the gas constant is J/(kg·K), the nozzle area expansion ratio is , and the Mach number of the throat is .
The theoretical calculation results for parameters at the nozzle outlet is obtained by solving the above formula, and the relative error is obtained by comparison with the numerical calculation result, as shown in Table 3. The numerical calculation results for pressure and temperature are larger than those of the theoretical calculation, and the numerical calculation results for speed and Mach number are smaller than those of the theoretical calculation. The isentropic flow calculation formula is used in the theoretical calculation, and the gas viscosity is not considered. However, the viscosity is considered in the numerical simulation, which causes the velocity at the nozzle outlet to be lower than the theoretical value and the pressure and temperature to be higher than the theoretical value. The error between the numerical calculation and the theoretical calculation results is within a reasonable range, which shows that the calculation model in this paper can effectively calculate the transient internal flow field of the star charge SRM.
4. Summary
In this study, a numerical calculation model for SRM burning surface regression was developed and implemented with the purpose of accurately and efficiently calculating the performance of the SRM. It is necessary for the interface to be clearly tracked during the regression of the burning surface to provide boundary conditions for the flow field calculation. Then, the gas flow in the flow field is calculated numerically according to the obtained boundary conditions, and the pressure, temperature, and velocity of the gas in the flow field can be solved through computational fluid dynamics. The CLSVOF method is used to track the interface, which is a coupled LS and VOF method. This method not only overcomes the shortcomings of the VOF method in that it is difficult to accurately calculate the normal and curvature of the interface but also solves the problem of nonconservation in the calculation procedure of the LS method. The porous media model is used to divide the fluid and solid calculation domains, which avoids the problem of mesh reconstruction in the calculation. The interface between the fluid domain and the solid domain is used to characterize the burning surface, and then, the area of the combustion surface at different moments can be obtained by solving the area of the interface.
The calculation model of this paper is verified by the zero-dimensional interior ballistic method and the CAD method. The results confirm that the proposed calculation model can compute the internal ballistics of SRM with different shapes of charges. However, some remaining issues need to be further studied to achieve more accurate calculations, such as erosion combustion and fluid-structure coupling calculations, of the performance of SRM.
Nomenclature
| : | Flow velocity (m/s) |
| : | Burning surface regression velocity (m/s) |
| : | Time (s) |
| : | Level set function (−) |
| : | Volume fraction (−) |
| : | Interface normal vector (−) |
| : | Compression coefficient (Pa) |
| : | Initial distance function |
| : | Characteristic cell size (−) |
| : | Infinitesimal (−) |
| : | Density (kg/m3) |
| : | Gravity (m/s2) |
| : | Temperature (K) |
| : | Dynamic viscosity (Pa·s) |
| : | Pressure (Pa) |
| : | Thermal conductivity (W/(m·K)) |
| : | Tangential stress (Pa) |
| : | Internal energy (J) |
| : | Mass source (kg/(m3·s)) |
| : | Momentum source (N/m3) |
| : | Energy source (W/m3) |
| : | Gas constant (J/(mol·K)) |
| : | Propellant density (kg/m3) |
| : | Propellant burning rate (m/s) |
| : | Specific heat capacity at constant pressure (J/(kg·K)) |
| : | Turbulence kinetic energy (m2/s2) |
| : | Turbulence dissipation rate (m2/s3) |
| : | Specific dissipation rate (1/s) |
| : | Eddy viscosity (Pa·s) |
| : | Mach number (−). |
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.