Coupling of Level Set and Volume of Fluid Methods for Simulations of Transient Internal Flow Field in Solid Rocket Motors

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.


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][4][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][7][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][10][11][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][16][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.

Coupling of Level Set and Volume of Fluid Formulation.
The VOF method tracks the interface by defining the volume fraction F. The volume fraction is defined as the ratio of the vol-ume 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 F = 1, the mesh is occupied by solid; if F = 0, the mesh is occupied by fluid; and if 0 < F < 1, 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 n is solved by the method discussed below. The compression term works at the interface only when FðF − 1Þ is not equal to 0; the compressibility coefficient is usually c F = 1.
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: <0, in discrete phase, = 0, at the interface, >0, in continuous phase: The governing equation of ϕ is To ensure the sharpness of the volume fraction F value obtained by the VOF algorithm and the boundedness of the LS distance function, the obtained volume fraction F 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 Γ = 0:75Δx, and Δx represents the mesh cell size. When the LS distance function is expressed by the F value, in order to ensure the accuracy of the ϕ 0 value, if the volume fraction in a mesh is less than 0 or greater than 1 in a time step, the flux of F 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 ϕ 0 ðxÞ is the numerical solution of the LS function at a certain moment, and ϕ 0 ðxÞ = 0 is expressed as the interface. Construct a function ϕðx, tÞ as the signed distance function at this moment and make ϕðx, tÞ = 0 still represent the interface at this moment. Perform the iteration of 2 International Journal of Aerospace Engineering the following equation: Equation (6) is solved until the steady state, and ϕðx, tÞ is used as the new signed distance function. The interface normal vector n 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 F n of the current step (step n) and the velocity information of the flow field (2) Map the volume fraction F n to the initial distance function ϕ 0 n by solving Equation (6) (3) Obtain the final distance function ϕ n+1 obtained by reinitializing the algorithm to solve Equation (7) (4) Solve the interface position by the obtained distance function ϕ n+1 , then solve Equation (5) to get the volume fraction F n+1 for the next step 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 S m and energy source S e are added in the mesh to simulate the injection of gas. The flow governing equations are as follows [19]: Level set function Re-initialization eq.  International Journal of Aerospace Engineering where S represents the burning surface area inside the mesh, as shown in Figure 2, and ΔV denotes the volume of the mesh. T tot is the total temperature, T ref represents the reference temperature, and T ref = 298:15 K. D represents the viscous resistance coefficient, and C represents the inertial resistance coefficient.

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 k-ω 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 k-ω model are as follows [20]: where u i ≡ U = ½u 1 , u 2 , u 3 T; G k represents the production of turbulence kinetic energy; G ω represents the generation of ω; Γ k and Γ ω , respectively, represent the effective diffusivity of k and ω; Y k and Y ω , respectively, represent the dissipation of k and ω due to turbulence; D ω represents the cross-diffusion term; and S k and S ω are user-defined source terms. The eddy viscosity can be calculated by the following formula: where S * indicates the strain rate magnitude and y 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, α * = 1.

Numerical Discrete.
For the spatial discretization, a structured multiblock finite volume approach with cellcentered 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.

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 D and inertial resistance coefficient C 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   International Journal of Aerospace Engineering 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 vis-cous resistance coefficient D and inertial resistance coefficient C need to meet the terms below. And Figure 3 demonstrates the distribution of D and C in different computational domains.
(1) For the solid domain, the values of D and C are set to very large numbers (such as 10 10    The volume fraction F is used to construct a function to calculate D and C at the interface, and two parameters R s and R f are introduced to make the function meet the above three terms, where R s = 10 10 and R f = 10 2 . The relevant formula is as follows: C has the same value as D, and their distribution curve is shown in Figure 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).  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.

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 e and divide it into k segments, where the thickness of each section will be Δe = e/m (2) When the burning thickness of the charge is iΔe ði = 0, 1, 2,⋯Þ, the surface area of the burning surface at this time is expressed as S i (i = 0, 1, 2, ⋯) (3) The pressure under each combustion layer is calculated by the formula p c , i = ðρ p c * a 1 S i /A t Þ 1/ð1−bÞ , where c * is the characteristic velocity of charge and A t denotes nozzle throat area (a) t = 0.2s   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     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.
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.

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 D = 50 mm, charge length L = 100 mm, characteristic length l = 10 mm, star tip arc radius r 1 = 3 mm, star root arc radius r 2 = 3 mm, convergence half angle of nozzle β = 45°, expansion half angle of nozzle α = 15°, and nozzle area expansion ratio ξ = 6. The propellant parameters are shown in Table 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 15: Contour of temperature distribution. 8 International Journal of Aerospace Engineering Figure 13, from thus illustrating which the entire process of the regression of the burning surface of the star charge. 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 M is defined as The relationship between pressure p and Mach number M: The relationship between temperature T and Mach number M: The formula for the speed of sound c: where the specific heat ratio is γ = 1:4, the total pressure is p 0 = 9:0 MPa, the total temperature is T 0 = 3000 K, the gas constant is R = 286:69 J/(kg·K), the nozzle area expansion ratio is ξ = 6, and the Mach number of the throat is M 0 = 1. 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