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

Authors: Ze Zhang; Feng Gao; Rui Lv; Yang Gao
Journal: International Journal of Aerospace Engineering
 (2021)

---

## 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.


---

## Body

## 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 Section2, 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 fractionF. 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]:
(1)∂F∂t+∇⋅UF=0,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:
(2)∂F∂t+∇⋅UF+∇⋅cFUnFF−1F=0,where n is solved by the method discussed below. The compression term works at the interface only when FF−1 is not equal to 0; the compressibility coefficient is usually cF=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:
(3)ϕx,t<0,indiscretephase,=0,attheinterface,>0,incontinuousphase.The governing equation ofϕ is
(4)∂ϕ∂t+V⋅∇ϕ=0.To ensure the sharpness of the volume fractionF 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:
(5)ϕ0=2F−1⋅Γ,where Γ=0.75Δx, and Δx represents the mesh cell size.When the LS distance function is expressed by theF 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ϕ0x is the numerical solution of the LS function at a certain moment, and ϕ0x=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 the following equation:
(6)∂ϕ∂t=Sϕ01−∇ϕ,(7)ϕx,0=ϕ0x,(8)Sϕ0=ϕ0ϕ02+ζ2.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:
(9)n=∇ϕ∇ϕ+ζ.Figure 1
Schematic diagram of the implementation process of the CLSVOF method.The main steps of CLSVOF method calculation are as follows:(1)
Calculate the volume fractionFn of the current step (step n) and the velocity information of the flow field(2)
Map the volume fractionFn to the initial distance function ϕ0n 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 Fn+1 for the next step(5)
After the distance functionϕn+1 and the volume fraction Fn+1 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 Lawr=a1pb, where p represents the pressure of the combustion chamber, a1 is the burning rate coefficient, and b 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 sourceSm and energy source Se are added in the mesh to simulate the injection of gas. The flow governing equations are as follows [19]:
(10)∂ρ∂t+∇ρU=Sm,∂ρE∂t+∇ρUE=−∇pU+∇τU+λ∇2T+Se,∂ρU∂t+∇ρUU=−∇p+μ∇2U+13μ∇∇U+Sv,p=ρRTg,Sm=ρp−ρSΔVr∇ϕ,Sv=−DμU−C12ρUU,Se=SmcpTtot−Tref,where S represents the burning surface area inside the mesh, as shown in Figure 2, and ΔV denotes the volume of the mesh. Ttot is the total temperature, Tref represents the reference temperature, and Tref=298.15 K. D represents the viscous resistance coefficient, and C represents the inertial resistance coefficient.Figure 2
Schematic diagram of burning surface in mesh.


#### 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 SSTk‐ω 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]:
(11)∂∂tρk+∂∂xiρkui=∂∂xiΓk∂k∂xi+Gk−Yk+Sk,∂∂tρω+∂∂xjρωuj=∂∂xjΓω∂ω∂xj+Gω−Yω+Dω+Sω,where ui≡U=u1,u2,u3T; Gk represents the production of turbulence kinetic energy; Gω represents the generation of ω; Γk and Γω, respectively, represent the effective diffusivity of k and ω; Yk and Yω, respectively, represent the dissipation of k and ω due to turbulence; Dω represents the cross-diffusion term; and Sk and Sω are user-defined source terms.The eddy viscosity can be calculated by the following formula:(12)μt=ρkω1max1/α∗,S∗F2/a1ω,F2=tanhΦ22,Φ2=max2k0.09ωy,500μρy2ω,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
(13)α∗=0.024+ρk/6μω1+ρk/6μω,where in the high-Reynolds number case, α∗=1.


#### 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)):
(14)Sv=−DμU−C12ρUU.The viscous resistance coefficientD 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 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 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 ofD and C 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 ofD and C 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 affectedFigure 3
Resistance coefficient in the computational domain.The volume fractionF is used to construct a function to calculate D and C at the interface, and two parameters Rs and Rf are introduced to make the function meet the above three terms, where Rs=1010 and Rf=102. The relevant formula is as follows:
(15)D=FF⋅1/Rs+1−F⋅1/Rf.C has the same value as D, and their distribution curve is shown in Figure 4.Figure 4
Distribution curve of resistance coefficient.


### 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 (Figure5).Figure 5
Flow chart showing the coupled computing scheme.


## 2.1. Coupling of Level Set and Volume of Fluid Formulation

The VOF method tracks the interface by defining the volume fractionF. 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]:
(1)∂F∂t+∇⋅UF=0,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:
(2)∂F∂t+∇⋅UF+∇⋅cFUnFF−1F=0,where n is solved by the method discussed below. The compression term works at the interface only when FF−1 is not equal to 0; the compressibility coefficient is usually cF=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:
(3)ϕx,t<0,indiscretephase,=0,attheinterface,>0,incontinuousphase.The governing equation ofϕ is
(4)∂ϕ∂t+V⋅∇ϕ=0.To ensure the sharpness of the volume fractionF 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:
(5)ϕ0=2F−1⋅Γ,where Γ=0.75Δx, and Δx represents the mesh cell size.When the LS distance function is expressed by theF 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ϕ0x is the numerical solution of the LS function at a certain moment, and ϕ0x=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 the following equation:
(6)∂ϕ∂t=Sϕ01−∇ϕ,(7)ϕx,0=ϕ0x,(8)Sϕ0=ϕ0ϕ02+ζ2.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:
(9)n=∇ϕ∇ϕ+ζ.Figure 1
Schematic diagram of the implementation process of the CLSVOF method.The main steps of CLSVOF method calculation are as follows:(1)
Calculate the volume fractionFn of the current step (step n) and the velocity information of the flow field(2)
Map the volume fractionFn to the initial distance function ϕ0n 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 Fn+1 for the next step(5)
After the distance functionϕn+1 and the volume fraction Fn+1 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 Lawr=a1pb, where p represents the pressure of the combustion chamber, a1 is the burning rate coefficient, and b 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 sourceSm and energy source Se are added in the mesh to simulate the injection of gas. The flow governing equations are as follows [19]:
(10)∂ρ∂t+∇ρU=Sm,∂ρE∂t+∇ρUE=−∇pU+∇τU+λ∇2T+Se,∂ρU∂t+∇ρUU=−∇p+μ∇2U+13μ∇∇U+Sv,p=ρRTg,Sm=ρp−ρSΔVr∇ϕ,Sv=−DμU−C12ρUU,Se=SmcpTtot−Tref,where S represents the burning surface area inside the mesh, as shown in Figure 2, and ΔV denotes the volume of the mesh. Ttot is the total temperature, Tref represents the reference temperature, and Tref=298.15 K. D represents the viscous resistance coefficient, and C represents the inertial resistance coefficient.Figure 2
Schematic diagram of burning surface in mesh.


### 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 SSTk‐ω 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]:
(11)∂∂tρk+∂∂xiρkui=∂∂xiΓk∂k∂xi+Gk−Yk+Sk,∂∂tρω+∂∂xjρωuj=∂∂xjΓω∂ω∂xj+Gω−Yω+Dω+Sω,where ui≡U=u1,u2,u3T; Gk represents the production of turbulence kinetic energy; Gω represents the generation of ω; Γk and Γω, respectively, represent the effective diffusivity of k and ω; Yk and Yω, respectively, represent the dissipation of k and ω due to turbulence; Dω represents the cross-diffusion term; and Sk and Sω are user-defined source terms.The eddy viscosity can be calculated by the following formula:(12)μt=ρkω1max1/α∗,S∗F2/a1ω,F2=tanhΦ22,Φ2=max2k0.09ωy,500μρy2ω,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
(13)α∗=0.024+ρk/6μω1+ρk/6μω,where in the high-Reynolds number case, α∗=1.


### 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.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 Lawr=a1pb, where p represents the pressure of the combustion chamber, a1 is the burning rate coefficient, and b 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 sourceSm and energy source Se are added in the mesh to simulate the injection of gas. The flow governing equations are as follows [19]:
(10)∂ρ∂t+∇ρU=Sm,∂ρE∂t+∇ρUE=−∇pU+∇τU+λ∇2T+Se,∂ρU∂t+∇ρUU=−∇p+μ∇2U+13μ∇∇U+Sv,p=ρRTg,Sm=ρp−ρSΔVr∇ϕ,Sv=−DμU−C12ρUU,Se=SmcpTtot−Tref,where S represents the burning surface area inside the mesh, as shown in Figure 2, and ΔV denotes the volume of the mesh. Ttot is the total temperature, Tref represents the reference temperature, and Tref=298.15 K. D represents the viscous resistance coefficient, and C represents the inertial resistance coefficient.Figure 2
Schematic diagram of burning surface in mesh.


## 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 SSTk‐ω 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]:
(11)∂∂tρk+∂∂xiρkui=∂∂xiΓk∂k∂xi+Gk−Yk+Sk,∂∂tρω+∂∂xjρωuj=∂∂xjΓω∂ω∂xj+Gω−Yω+Dω+Sω,where ui≡U=u1,u2,u3T; Gk represents the production of turbulence kinetic energy; Gω represents the generation of ω; Γk and Γω, respectively, represent the effective diffusivity of k and ω; Yk and Yω, respectively, represent the dissipation of k and ω due to turbulence; Dω represents the cross-diffusion term; and Sk and Sω are user-defined source terms.The eddy viscosity can be calculated by the following formula:(12)μt=ρkω1max1/α∗,S∗F2/a1ω,F2=tanhΦ22,Φ2=max2k0.09ωy,500μρy2ω,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
(13)α∗=0.024+ρk/6μω1+ρk/6μω,where in the high-Reynolds number case, α∗=1.


## 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)):
(14)Sv=−DμU−C12ρUU.The viscous resistance coefficientD 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 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 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 ofD and C 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 ofD and C 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 affectedFigure 3
Resistance coefficient in the computational domain.The volume fractionF is used to construct a function to calculate D and C at the interface, and two parameters Rs and Rf are introduced to make the function meet the above three terms, where Rs=1010 and Rf=102. The relevant formula is as follows:
(15)D=FF⋅1/Rs+1−F⋅1/Rf.C has the same value as D, and their distribution curve is shown in Figure 4.Figure 4
Distribution curve of resistance coefficient.


## 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 (Figure5).Figure 5
Flow chart showing the coupled computing scheme.


## 3. Numerical Tests and Discussion

### 3.1. Tubular Charge SRM Transient Internal Flow Field Simulations

#### 3.1.1. Physical Model

Figure6 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.Figure 6
Schematic diagram of the SRM model with tubular charge.Table 1
Propellant parameters.
ValuePropellant density (kg/m3)1800Burning rate coefficient0.015Burning rate pressure exponent0.3


#### 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 ise 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 isiΔei=0,1,2,⋯, the surface area of the burning surface at this time is expressed as Si (i=0,1,2,⋯)(3)
The pressure under each combustion layer is calculated by the formulapc,i=ρpc∗a1Si/At1/1−b, where c∗ is the characteristic velocity of charge and At denotes nozzle throat area(4)
The burning rate at stepi is ri=a1pc,i+pc,i+1/2bFigure7 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.Figure 7
Contours of pressure distribution at different time points.Figure 8
Comparison of combustion chamber pressure.As presented in Figure9, 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.Figure 9
Schematic diagram of burning surface regression.
(a)
t=0.2s(b)
t=0.8s(c)
t=1.6s(d)
t=2.0s(e)
t=2.6sFigure 10
Comparison of burning surface area.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. Figure11 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 r1=3 mm, star root arc radius r2=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.Figure 11
Cross-section of the star charge.Table 2
Propellant parameters.
ValuePropellant density (kg/m3)1800Burning rate coefficient0.015Burning rate pressure exponent0.3


#### 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 Figure12, 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.Figure 12
Comparison of burning surface area.Figure 13
Schematic diagram of burning surface regression.
(a)
t=0.01s(b)
t=0.4s(c)
t=0.8s(d)
t=1.2s(e)
t=1.6sFigure14 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.Figure 14
Contour of pressure distribution.Figure 15
Contour of temperature distribution.Figure 16
Contour of Mach number distribution.Figure 17
Contour of 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
(16)ξ=M0M1+γ−1/2M21+γ−1/2M02γ+1/2γ−1.The relationship between pressurep and Mach number M:
(17)pp0=1+γ−12M2−γ/γ−1.The relationship between temperatureT and Mach number M:
(18)TT0=1+γ−12M2−1.The formula for the speed of soundc:
(19)c=γRT,where the specific heat ratio is γ=1.4, the total pressure is p0=9.0 MPa, the total temperature is T0=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 M0=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 Table3. 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.Table 3
Propellant parameters.
ParameterTheoretical calculation valueNumerical calculation valueRelative errorPressure (MPa)1.431.451.39%Temperature (K)917.85937.462.13%Velocity (m/s)2044.141999.572.18%Mach number3.373.282.67%


## 3.1. Tubular Charge SRM Transient Internal Flow Field Simulations

### 3.1.1. Physical Model

Figure6 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.Figure 6
Schematic diagram of the SRM model with tubular charge.Table 1
Propellant parameters.
ValuePropellant density (kg/m3)1800Burning rate coefficient0.015Burning rate pressure exponent0.3


### 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 ise 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 isiΔei=0,1,2,⋯, the surface area of the burning surface at this time is expressed as Si (i=0,1,2,⋯)(3)
The pressure under each combustion layer is calculated by the formulapc,i=ρpc∗a1Si/At1/1−b, where c∗ is the characteristic velocity of charge and At denotes nozzle throat area(4)
The burning rate at stepi is ri=a1pc,i+pc,i+1/2bFigure7 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.Figure 7
Contours of pressure distribution at different time points.Figure 8
Comparison of combustion chamber pressure.As presented in Figure9, 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.Figure 9
Schematic diagram of burning surface regression.
(a)
t=0.2s(b)
t=0.8s(c)
t=1.6s(d)
t=2.0s(e)
t=2.6sFigure 10
Comparison of burning surface area.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.1.1. Physical Model

Figure6 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.Figure 6
Schematic diagram of the SRM model with tubular charge.Table 1
Propellant parameters.
ValuePropellant density (kg/m3)1800Burning rate coefficient0.015Burning rate pressure exponent0.3


## 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 ise 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 isiΔei=0,1,2,⋯, the surface area of the burning surface at this time is expressed as Si (i=0,1,2,⋯)(3)
The pressure under each combustion layer is calculated by the formulapc,i=ρpc∗a1Si/At1/1−b, where c∗ is the characteristic velocity of charge and At denotes nozzle throat area(4)
The burning rate at stepi is ri=a1pc,i+pc,i+1/2bFigure7 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.Figure 7
Contours of pressure distribution at different time points.Figure 8
Comparison of combustion chamber pressure.As presented in Figure9, 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.Figure 9
Schematic diagram of burning surface regression.
(a)
t=0.2s(b)
t=0.8s(c)
t=1.6s(d)
t=2.0s(e)
t=2.6sFigure 10
Comparison of burning surface area.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. Figure11 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 r1=3 mm, star root arc radius r2=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.Figure 11
Cross-section of the star charge.Table 2
Propellant parameters.
ValuePropellant density (kg/m3)1800Burning rate coefficient0.015Burning rate pressure exponent0.3


### 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 Figure12, 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.Figure 12
Comparison of burning surface area.Figure 13
Schematic diagram of burning surface regression.
(a)
t=0.01s(b)
t=0.4s(c)
t=0.8s(d)
t=1.2s(e)
t=1.6sFigure14 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.Figure 14
Contour of pressure distribution.Figure 15
Contour of temperature distribution.Figure 16
Contour of Mach number distribution.Figure 17
Contour of 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
(16)ξ=M0M1+γ−1/2M21+γ−1/2M02γ+1/2γ−1.The relationship between pressurep and Mach number M:
(17)pp0=1+γ−12M2−γ/γ−1.The relationship between temperatureT and Mach number M:
(18)TT0=1+γ−12M2−1.The formula for the speed of soundc:
(19)c=γRT,where the specific heat ratio is γ=1.4, the total pressure is p0=9.0 MPa, the total temperature is T0=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 M0=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 Table3. 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.Table 3
Propellant parameters.
ParameterTheoretical calculation valueNumerical calculation valueRelative errorPressure (MPa)1.431.451.39%Temperature (K)917.85937.462.13%Velocity (m/s)2044.141999.572.18%Mach number3.373.282.67%


## 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. Figure11 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 r1=3 mm, star root arc radius r2=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.Figure 11
Cross-section of the star charge.Table 2
Propellant parameters.
ValuePropellant density (kg/m3)1800Burning rate coefficient0.015Burning rate pressure exponent0.3


## 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 Figure12, 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.Figure 12
Comparison of burning surface area.Figure 13
Schematic diagram of burning surface regression.
(a)
t=0.01s(b)
t=0.4s(c)
t=0.8s(d)
t=1.2s(e)
t=1.6sFigure14 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.Figure 14
Contour of pressure distribution.Figure 15
Contour of temperature distribution.Figure 16
Contour of Mach number distribution.Figure 17
Contour of 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
(16)ξ=M0M1+γ−1/2M21+γ−1/2M02γ+1/2γ−1.The relationship between pressurep and Mach number M:
(17)pp0=1+γ−12M2−γ/γ−1.The relationship between temperatureT and Mach number M:
(18)TT0=1+γ−12M2−1.The formula for the speed of soundc:
(19)c=γRT,where the specific heat ratio is γ=1.4, the total pressure is p0=9.0 MPa, the total temperature is T0=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 M0=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 Table3. 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.Table 3
Propellant parameters.
ParameterTheoretical calculation valueNumerical calculation valueRelative errorPressure (MPa)1.431.451.39%Temperature (K)917.85937.462.13%Velocity (m/s)2044.141999.572.18%Mach number3.373.282.67%


## 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.



---

*Source: 2984992-2021-06-19.xml*