Fluid–solid coupled simulation of the ignition transient of solid rocket motor

Acta Astronautica 110 (2015) 180–190

Contents lists available at ScienceDirect

Acta Astronautica journal homepage: www.elsevier.com/locate/actaastro

Fluid–solid coupled simulation of the ignition transient of solid rocket motor Qiang Li n, Peijin Liu, Guoqiang He Science and Technology on Combustion, Internal Flow and Thermal-Structure Laboratory, Northwestern Polytechnical University, Xi’an, Shaanxi 710072, PR China

a r t i c l e in f o

abstract

Article history: Received 8 December 2014 Received in revised form 24 January 2015 Accepted 28 January 2015 Available online 4 February 2015

The first period of the solid rocket motor operation is the ignition transient, which involves complex processes and, according to chronological sequence, can be divided into several stages, namely, igniter jet injection, propellant heating and ignition, flame spreading, chamber pressurization and solid propellant deformation. The ignition transient should be comprehensively analyzed because it significantly influences the overall performance of the solid rocket motor. A numerical approach is presented in this paper for simulating the fluid– solid interaction problems in the ignition transient of the solid rocket motor. In the proposed procedure, the time-dependent numerical solutions of the governing equations of internal compressible fluid flow are loosely coupled with those of the geometrical nonlinearity problems to determine the propellant mechanical response and deformation. The wellknown Zeldovich–Novozhilov model was employed to model propellant ignition and combustion. The fluid–solid coupling interface data interpolation scheme and coupling instance for different computational agents were also reported. Finally, numerical validation was performed, and the proposed approach was applied to the ignition transient of one laboratory-scale solid rocket motor. For the application, the internal ballistics were obtained from the ground hot firing test, and comparisons were made. Results show that the integrated framework allows us to perform coupled simulations of the propellant ignition, strong unsteady internal fluid flow, and propellant mechanical response in SRMs with satisfactory stability and efficiency and presents a reliable and accurate solution to complex multi-physics problems. & 2015 IAA. Published by Elsevier Ltd. All rights reserved.

Keywords: Computational fluid dynamics Finite elements method Fluid–solid interaction Coupling instance

1. Introduction The ignition transient of solid rocket motor (SRM) operation is complex, and it involves time-dependent processes that advance in time in a strongly coupled manner [1,2]. For example, in a real-world SRM with complex propellant configuration, the hot igniter jet injected from the igniter flows down the chamber bore and heats up the solid

n

Corresponding author. Fax: þ86 2988494163. E-mail address: [email protected] (Q. Li).

http://dx.doi.org/10.1016/j.actaastro.2015.01.017 0094-5765/& 2015 IAA. Published by Elsevier Ltd. All rights reserved.

propellant. Once the propellant surface temperature reaches the critical temperature, the propellant is ignited locally, and high-temperature gaseous combustion products are released from the propellant surface and injected into the combustion chamber. The time-dependent internal flow field, propellant heating and ignition, and flame spreading are fully coupled because the development of the internal flow field in the combustion chamber governs the convective and radiative heat transfer from the high-temperature combustion products to the propellant solid and the propellant ignition, which determines the flame spreading speed and the propellant burning surface ignition sequence, and in turn,

Q. Li et al. / Acta Astronautica 110 (2015) 180–190

the evolution of the internal flow field. Furthermore, in the ignition transient, the motor chamber pressurization loadings applied on the propellant surface cause deformations to the solid propellant that change the geometrical configurations of both the fluid subdomain and the solid subdomain. Thus, the development of internal flow field and the mechanical response of the propellant are tightly coupled. Previous studies that have investigated the ignition transient of SRMs can be classified into two main categories. In the first category, the internal flow fields in the motor ignition transient were numerically simulated with propellant ignition models of different complexities, and the mechanical responses of the solid propellant were not considered [1–3]. In the second category, the unsteady numerical calculations of the motor ignition transient were not conducted, and the mechanical responses of the solid propellant were simulated with different ignition pressurization loadings [4–6] that were applied on the propellant burning surface. These previous works serve as valuable references values by establishing novel scientific ideas that future studies could apply both in their methodology and theoretical analysis. However, few studies have been reported on the fluid–solid fully coupled simulation of the motor ignition transient. The key problems for fluid–solid fully coupled simulation of the SRM ignition transient consist of a precise and robust solver for the strong unsteady compressible flow; an ignition model with satisfactory accuracy for predicting the ignition sequence of the propellant burning surface and flame spreading in the chamber; a mesh smoothing technique for maintaining and improving the quality of unstructured volume meshes with severe deformations; an interpolation scheme with satisfactory accuracy and low cost for interpolating data conservatively across the fluid–solid coupling interfaces; and a coupling procedure for handling the coupling of different agents (fluid flow, solid mechanics, propellant ignition, and fluid–solid coupling interface data transfer) at different time and spatial scales as well as for improving the stability and efficiency of the simulations [7]. This study focuses on the mathematical models and numerical strategies for the fluid–solid fully coupled simulation of the ignition transient of SRM. The remainder of this paper is organized as follows. Section 2 revisits the governing equations and numerical methods employed for the compressible flow and solid mechanics in SRMs. Section 3 introduces the ignition model used. Section 4 documents the fluid–solid coupling procedure. Section 5 reports the validation of the numerical strategies and the application and numerical results of a laboratory-scale SRM.

181

The application of the arbitrary-Lagrangian–Eulerian (ALE) formulation allows computation of these equations. The set of governing equations for the compressible viscous flow in the ALE formulation can be expressed as [8–11] Z I I ∂ U dΩ þ Fc dS ¼ Fv dS ð1Þ ∂t Ω ∂Ω ∂Ω where Ω is the control volume, ∂Ω is the boundary of Ω, U is the vector of the conservative variables, and Fc is the convective flux vector and can be written as follows: 2 3 ρ 6 ρu 7 6 7 6 7 7 U¼6 ð2Þ 6 ρv 7 6 ρw 7 4 5 ρE 2

ρV

3

6 ρuV þ nx p 7 6 7 6 7 ρvV þ ny p 7 Fc ¼ 6 6 7 6 ρwV þ n p 7 z 5 4 ρHV þV g p

ð3Þ

where ρ and p are the fluid density and pressure; u; v; w are the fluid velocity components; and V is the velocity of fluid flow relative to the computational mesh and can be expressed by the following equation: V ¼ nx u þ ny v þnz w  V g

ð4Þ

where nx, ny, and nz are the unit normal vector components of the control volume face; V g is the grid velocity in the normal direction of the control volume face; E is the specific total energy of the fluid; and H is the specific stagnation enthalpy of the fluid. The viscous flux vector can be written as follows: 2 3 0 6 nx τxx þny τxy þ nz τxz 7 7 6 6 7 nx τyx þ ny τyy þnz τyz 7 Fv ¼ 6 ð5Þ 6 7 6 n τ þn τ þn τ 7 y zy z zz 5 4 x zx nx θx þ ny θy þ nz θz where τij represents the viscous stresses. It can be written as follows:     ∂vi ∂vj 2μ ∂vk τij ¼ μ þ δ ð6Þ  3 ∂xk ij ∂xj ∂xi Moreover, ∂T θx ¼ uτxx þvτxy þ wτxz þ k ∂x

∂T ∂y ∂T θz ¼ uτzx þ vτzy þ wτzz þ k ∂z

θy ¼ uτyx þ vτyy þ wτyz þk 2. Governing equations and numerical methods 2.1. Fluid subdomain In the fluid–solid coupled simulation of the motor ignition transient, the governing equations for the fluid flow must be formulated to address the compressible flow on moving meshes caused by the dynamic changes in the fluid subdomain geometry.

ð7Þ

are terms describing the work of viscous stresses and heat conduction in the fluid. In this paper, μ and k are the viscosity and thermal conductivity of the combustion gases in the SRM, T is the fluid temperature and δij is the Dirac function. For the governing equations, the Spalart–Allmaras one-equation turbulence model, which is of satisfactory

182

Q. Li et al. / Acta Astronautica 110 (2015) 180–190

accuracy for wall-bounded compressible flow, was used for model closure. In the present analysis, finite volume approach with cell-centered data storage was adopted because the method is known for easily satisfying conservation laws. For the convective fluxes in Eq. (1), the finite volume code uses the second-order Harten–Lax–van Leer contact wave scheme to address strong transient flows, whereas for the viscous fluxes, second-order central difference scheme is used for accurate solutions of the viscous terms. The discrete equations are integrated in time by using a fourstep, fourth-order explicit Runge–Kutta method. The time step is chosen as the minimum value of maximum stable time step in all fluid cells. 2.2. Solid subdomain In the ignition transient of SRMs, the mechanical responses of the solid propellant are controlled by the ignition pressurization loadings and solid propellant material properties. Transient response analysis is the most general method for evaluating mechanical response with timevarying loadings. For that all the modern solid propellant grains utilize an elastomeric binder, which is filled with a high level of solid particles. The application of a load causes different complicated mechanisms to take place in the binder, the filler or the interface between them. Under these influences solid propellant grains exhibit very complex behavior including features associated with time and rate effects, temperature and pressure dependence, large deformations and strains, etc. Therefore the attempt to represent all aspects of solid propellant grain behavior would result in a very complicated constitutive model [12–14]. In the present paper, a number of previous investigations have been concerned with certain features only and special attentions are paid to the nonlinearity of the problem confronted. In the ignition transient of SRM, boundary nonlinearities could be neglected because the solid propellant is tightly bounded to the motor case, and the bonding strength is strong enough and no contact problem is presented. Meanwhile, material nonlinearities can also be ignored because the time period of ignition transient is fairly short, the inherent viscoelastic feature of the solid propellant cannot be exhibited, and a linear elastic constitutive model is highly acceptable although modern solid propellant is complicated in both chemical compositions and constitutive models. Geometric nonlinearities must be formulated in the numerical model because the displacements of solid propellant subjected to ignition pressurization loadings may be large, and the small displacement assumptions inherent in the equations of linear analysis may be invalid. Assembling all the element, the finite element equations for the geometrical nonlinear problems read as follows [5,15,16]: Mq€ ðtÞ þ Cq_ ðtÞ þKqðtÞ ¼ PðtÞ

ð8Þ

where qðtÞ is the displacement vector, the single and double dots over the displacement vector qðtÞ denote derivatives with respect to time for velocity and acceleration, PðtÞ is the external loading acting on the structure, M is the mass matrix, C is the damping matrix and is approximated by

Rayleigh proportional damping method, and K is the stiffness matrix. The Poisson's ratio of the solid propellant involved in this paper is 0.496. Therefore, the displacement-based finite element formulation Eq. (8) can be employed. The initial conditions are as follows: qð0Þ ¼ q0

ð9Þ

_ qð0Þ ¼ dq0 =dt

ð10Þ

0

0

where q is the initial displacement vector, and dq =dt is the initial velocity vector. Different numerical algorithms have been proposed for the nonlinear finite element equations [17]. In this study, the well-known Newmark algorithm [17] is employed to solve the initial condition problem. Given the time interval ½t; t þ Δt , on knowing the value and derivative of displacement vector qðtÞ, the objective is to determine _ þΔtÞ, and q€ ðt þΔtÞ, which satisfy Eq. (8). The qðt þ ΔtÞ, qðt following equations can be formulated by using the truncated Taylor series expansion.   q_ t þ Δt ¼ q_ t þ ð1 βÞq€ t þ βq€ t þ Δt Δt ð11Þ qt þ Δt ¼ qt þ q_ t Δt þ q€ t þ Δt ¼

   1  α q€ t þ αq€ t þ Δt Δt 2 2

  1  1 1 _t  q  1 q€ t q  q  t þ Δt t αΔt 2α αΔt 2

ð12Þ

ð13Þ

where α and β are two constants, and the subscript t and t þ Δt represent the values at time t and t þ Δt, respectively. In this paper, the calculation involves obtaining the unknowns by using an incremental Newton–Raphson approach [17], and the geometrical nonlinearity is modeled by calculating the shape function matrix and the stiffness matrix for the deformed configuration at each load step. 3. Solid propellant ignition model The ignition of solid propellant is extremely complex. It involves complex physiochemical processes, such as propellant heating, decomposition, gasification, and surface gas phase chemical reactions [3]. Under rapid pressurization rate in the motor ignition transient, the propellant dynamic burning rate can be significantly different from the burning rate predicted by the quasi-steady combustion model because of the thermal relaxation effect in the solid [14,18]. Ignition models of different complexity have been presented and used previously in numerical simulations [14,18–20]. The simplest model switches on the burning boundary condition when the temperature of hot gases flowing over the propellant surface reaches a critical ignition temperature and an empirical formulation was utilized to calculate the mass flow rate of the combustion products. This model actually works well with the exception that it underestimates the ignition delay due to propellant heating and combustion initiation. The more detailed ignition models are extremely complex, such as the heterogeneous reaction models [12,14], which consider the condensed-phase energy conservation, gas-phase species diffusion and conservation, gas–solid interfacial boundary conditions, ambient oxidizer species

Q. Li et al. / Acta Astronautica 110 (2015) 180–190

183

diffusion, and the detailed heterogeneous reactions at the propellant surface reaction zone. The well-known and widely accepted Zeldovich–Novozhilov (ZN) model [21] was employed in this study to simulate the ignition process. ZN model involves the utilization of a one-dimensional unsteady heat conduction equation to determine the propellant surface temperature. This equation is formulated in a semi-infinite domain in the normal direction of the propellant surface. With a uniform temperature as an initial condition and an unsteady heat flux as a boundary condition, the one-dimensional unsteady heat conduction equation was numerically solved by the third-order central difference scheme. Propellant surface temperature ignition criterion was used. This criterion refers to a point on the propellant surface that is ignited when the temperature at that point reaches a critical ignition temperature. The one-dimensional unsteady heat conduction equation reads as follows [22,23]

4. Fluid–solid coupling strategies

∂T s ∂T s ∂2 T s þ r_ ¼ αp 2 ∂t ∂x ∂x

where Γ is the fluid–solid coupling interface, the superscript f and s represent the values associated with the fluid and solid, respectively, n is the unit outward normal of fluid and solid subdomain, and σ is solid stress, as shown in Fig. 1, ρp is the propellant density, ui is the velocity component in the i-th direction [26,27]. In the equilibrium conditions, the action of the shear force applied on the fluid–solid coupling interface by the gas flow is neglected because this shear force is small compared with the pressure force [28]. Thus, the duration of ignition transient of SRM is short, and the regression distance caused by propellant combustion is not included although the steady state regression rate of solid propellant is up to 48.0 mm/s.

T s ðx; 0Þ ¼ T is  λp

∂T ð0; tÞ ¼ hðT  T b Þ ∂x

T s ð  1; t Þ ¼ T 0s

ð14Þ

where T s is the temperature distribution in the solid propellant, r_ is the propellant burning rate, and αp is the propellant heat diffusion coefficient, λp is the propellant heat transfer coefficient, T is the gas temperature in the adjacent cell, T b is the propellant surface temperature. The superscript i and 0 represent initial and far-field condition. A comprehensive description of the numerical method on Eq. (14) can be referred to [24]. In numerical simulation, boundary conditions must be provided to solve Eq. (14). In this study, the semi-empirical convective heat transfer and radiative heat transfer coefficients reported by Johnston [25] were used to evaluate the heat flux. The local convective heat transfer model used in this study is as follows:  0:2 hc ¼ ð0:023ÞP r  2=3 C p μ=D ðρVÞ0:8 where P r is the Prandtl number, C p is the gas specific heat, μ is the gas viscosity, V is the gas velocity in the adjacent cell, D is the reference length. For surfaces in the bore of the motor, the local bore diameter is used for the reference length D. The radiative heat transfer model used in this study is as follows: hr ¼ C sf σðT 2 þ T 2b ÞðT þ T b Þ where σ is the Stefan–Boltzman constant and C sf is an empirical coefficient accounting for transfer factors and set to 0.25. ZN model is intrinsically one dimensional (into the propellant) and formulated in the normal direction of the propellant surface. In calculation, this model is employed independently at each fluid cell face on the burning propellant surface, thereby making this model effectively three dimensional.

4.1. Fluid–solid coupling interface boundary condition For the fluid–solid coupled system, the geometrical compatibility conditions and the equilibrium conditions on the coupling interface must be satisfied. In this study, the boundary conditions applied on the fluid–solid coupling interface are as follows: (i) compatibility conditions:





qf ¼ qs

Γ

ð15Þ

Γ

(ii) equilibrium conditions:







pnfj þ σ sji nsj þ ρu2i ρ=ρp  1 ¼ 0 Γ

Γ

i ¼ 1; 2; 3

ð16Þ

4.2. Fluid–solid interface data transfer In fluid–solid coupled simulations, determining ideal interpolation scheme to interpolate quantities across the coupling interface is of great importance because the fluid subdomain mesh and the solid subdomain mesh are generated independently by different software products to satisfy the special requirements of the solvers. Although the fluid interface and the solid interface share the same geometry where the two sub-domains interact with each other, the fluid interface mesh and the solid interface mesh are potentially non-matching, as shown in Fig. 1. An accurate interpolation scheme with low computational cost to interpolate quantities from one side of coupling interface to another with satisfactory accuracy is important [29,30]. Interpolating quantities across the coupling interface requires calculating quantities at each target mesh points based on the data associated with the source mesh points. In this paper, solution interpolation involves three steps for the target mesh.

Fig. 1. Coupling interface boundary conditions.

184

Q. Li et al. / Acta Astronautica 110 (2015) 180–190

First, the variables associated with the source mesh are interpolated to cell vertices by weight, with the weighting function for vertex i of cell k is defined as follows: wi ¼

Ai Ak

ð17Þ

where Ai is the area of the triangle formed by the center and two other vertices of cell k, and Ak is the cell k area. The interpolated data is assigned as the corresponding value of vertex i associated with the cell k. For each target mesh cell vertex j, the cell containing the target mesh point in the source mesh is located by the octree method, and the data associated with the source mesh cell vertices are interpolated to the target mesh cell vertex by weight. The weighting function is defined in the same way as in Eq. (17) with the coordinates of vertex j as the cell center. Finally, for the target mesh, the cell data are interpolated from the cell vertices to the cell centers by weight. The weighting function is defined in the same way as in Eq. (17), and other variables are then inferred. The interpolation scheme used in this study is the classic scheme widely employed in the finite element method [17]. This scheme is a linear interpolation and first order in accuracy. 4.3. Coupling procedure In fluid–solid coupled simulation of the ignition transient of SRM, the two-way coupling instance is of special importance. Given that the nonlinear governing equations of the fluid and solid subdomains are potentially different, these two sets of equations are extremely difficult to discretize and solve in the same computational frame. In addition, these two subdomains provide boundary conditions for each other. Thus, none of the set of equations can be solved independent of the other. Moreover, this type of simulation often involves several modules that interact with one another. The instances at which information is transferred from one module to the other and the arrangement of these modules advance in time have a deep effect on the stability and efficiency of the coupling method and on the reliability and accuracy of the results. The coupling method used to satisfy the geometrical compatibility and the equilibrium conditions on the coupling interface between the two subdomains and the coupling procedure are tough challenges confronted. The two-way coupling method used in fluid–solid coupled simulations can be classified into three main categories [31–33]. The first category is known as the monolithic approach or direct coupling method. In this method, the coupled equations are solved in a direct way, and the variables of the coupled system are updated simultaneously. The main drawbacks of this method consist of complexity in the model derivation, solver coding, less flexibility of numerical schemes available, high-memory requirement, and computation cost. The second category is termed as the partitioned methods or the iterative method. In this method, separate governing equations for the fluid and structure are discretized separately in each subdomains and are coupled by using a synchronization procedure both

in time and space by subiteration. In applications, the partitioned method may suffer from misconvergence and instability, and a fully implicit formulation must be used for stability reasons. This formulation usually increases the memory requirement and computation cost dramatically. The third category is the most practical approach, which is named as the loose coupling method. In this method, the governing equations for the fluid and solid are discretized independently in each subdomains by using the preferred numerical schemes individually, and the information is passed in two ways by using a synchronization procedure both in time and space without subiteration. Loose coupling method is the simplest and the best-known coupling algorithm and has also been extensively used in various fluid– structure-interaction problems in different engineering disciplines in the past 10 years. In this paper, a loose coupling instance is presented, as shown in Fig. 2 for its obvious advantages in both the accuracy and easy implementation of the algorithm. First, a system time step dt is prescribed empirically during calculation, which is usually approximately 10  7  10  6 s, to evolve the system from time level n to the new time level nþ 1. For the propellant ignition, fluid flow, and solid mechanics solvers, the time steps adopted internally are determined independently by the local Courant–Friedrichs–Lewy conditions within the corresponding applications. In one system time step, the calculation begins with the propellant ignition solver. This solver takes one or more internal steps to evolve the solutions to the new time level n þ1. Then, the combustion product mass flux mb and stagnation temperature T n , which is calculated from solutions at time level n þ1, are directly copied to the fluid solver for that the mesh entities of the propellant ignition module are gained directly from the fluid subdomain surface meshes and that they match well on the coupling interface. With the new applied boundary conditions, the fluid solver then takes one or more internal steps to advance the fluid subdomain solution to the new time level nþ1. Then the new surface pressure p is interpolated as coupling interface boundary conditions to the solid mechanics solver. The solid mechanics solver then takes several internal steps to evolve the solid subdomain solution to the new time level nþ1, and all subdomain solutions are now at the new time level nþ1. Finally, the fluid–solid coupling interface displace ment vector U x ; U y ; U z is interpolated from the solid mechanics solver to the fluid solver. Meanwhile, the combustion product pressure p and the computed convective and radiation heat fluxes qc , qr are copied from the fluid

Fig. 2. Schematic of the coupling method.

Q. Li et al. / Acta Astronautica 110 (2015) 180–190

185

solver to the propellant ignition solver. Then, the preparation is completed for the next system time step. During motor ignition transient, the solid subdomain deforms under the action of fluid flow pressure. Consequently, the fluid subdomain surface mesh also moves accordingly because the two subdomains must be geometrically compatible on the coupling interface. Given that the interior nodes of fluid subdomain are intact, which may give rise to mesh skewness and distortion, fluid volume mesh smooth is necessary. In this study, the local volume mesh smoothing and the corresponding volume data interpolation scheme presented in [11] are used. 5. Numerical results and discussions In this section, numerical validation of the numerical method and the coupling procedure was performed, and then the application was conducted to a laboratory-scale SRM to demonstrate the effectiveness and accuracy of the procedure.

Fig. 3. Comparison of the shockwave structures. (a) Experimental result and (b) numerical result.

5.1. Numerical validation In this subsection, numerical validation was performed against the experiment conducted by Giordano [34]. Details about the experimental setup and the parameters can be referred to [34]. In the validation in this paper, the length and the thickness of the panel are 50.0 mm and 1.0 mm, respectively. The shock wave moves from the inlet at Mach number of 1.21 in air at rest (1.105 Pa and 293 K). Both the fluid subdomain and the solid subdomain are discretized by tetrahedrons. The solid subdomain is composed of 113,500 cells, whereas the fluid subdomain is composed of 426,500 cells. Fig. 3 shows the shockwave structures of the experimental and numerical results, and Fig. 4 shows a comparison of the time–history curves of the panel tip displacement. The numerical results match well with the experimental results. This finding indicates that the numerical method and the corresponding coupling procedure presented in this paper are reliable and accurate.

Fig. 4. Time–history curves of the panel tip displacement.

5.2. Numerical simulations In this subsection, multiphysics simulation of the ignition transient of a laboratory-scale test SRM was conducted. Fig. 5 shows the schematic of the test motor, which consists of igniter, thick steel case including the head-end plate and cylinder, solid propellant, combustion chamber, and nozzle assembly. The igniter is enclosed in the head end plate and designed in such a way so that the initial velocity of the igniter jet is parallel to the motor axis. The nozzle throat is made up of carbon/carbon composite and ablation can be totally neglected in the ignition transient. The configuration of the solid propellant grain is simple and can be divided into three main sections. The first section, which is close to the motor head end, is cylindrical. The second section is conical, and the third section is cylindrical, which is close to the aft end. The propellant configuration presented is slightly simple compared with that of a real-world SRM. However, the simulation of this

Fig. 5. Schematic of the laboratory-scale motor.

laboratory-scale motor essentially poses all of the major challenges encountered in a real-world SRM given that the strong unsteady flow, convective, and radiative heat transfer in cylindrical and conical channels, propellant ignition, and flame spreading are all included. The constraints of the propellant by the motor case are also modeled. The computational subdomains are built from the geometry of the experimental motor setup. This setup features the combustion chamber, motor head and aft end, the igniter, and the nozzle. The computational subdomains are discretized by unstructured meshes, and the discrete

186

Q. Li et al. / Acta Astronautica 110 (2015) 180–190

geometrical model is composed of 1,106,200 tetrahedrons for the fluid subdomain and 123,400 tetrahedrons for the solid subdomain. To define the boundary conditions, the igniter exit is defined as an inlet with a prescribed mass flow rate and gas temperature. The propellant burning surface is also defined as an inlet, and both the mass flow rate and the temperature of the combustion products are calculated with the propellant ignition model. The initial velocity of combustion products is specified as normal to the burning surface. The nozzle exit is defined as an outlet, and the ambient pressure is applied. The motor case is thick and strong. Consequently, the deformation of the motor case is neglected, and the outer surface of the solid propellant grain is fixed in all directions. In this paper, the equilibrium compositions of the combustion products were calculated by the Chemical Equilibrium with Application Code, a free software for calculating chemical equilibrium products based on the free-energyminimization principle [35]. The estimated chamber pressure (22.0 Mpa) was used to calculate chemical equilibrium. Results show that the combustion products in the combustion chamber comprise various species, as shown in Table 1. Other parameters such as chamber pressure, total temperature and the propellant combustion law are shown in Table 2. In SRM, the internal ballistics, which refers to the curve of the chamber pressure with time, is a key test parameter that can be accurately measured during hot firing test. During calculation, the pressure time–history data at the SRM head end and aft end are both probed and recorded, also the head end pressure curve was predicted theoretically [36]. The head end theoretical predication, numerical data and experimental data with standard error bars are plotted in Fig. 6 along with the aft end numerical data. The head end theoretical predication, numerical and experimental data curves evidently match well, and the numerical data lie within the error bars of the experimental data. The calculated maximum relative error between the head end numerical data and experimental data is approximately 12.6% in the time interval from 4.2 ms to 8.0 ms, during which the recorded result is larger than the calculated result. The main factors can be attributed to the errors in the evaluation of heat transfer coefficient. In this paper, the calculation involves the evaluation of the convective heat transfer coefficients by a semi-empirical expression [25], which is essentially for the quasi-steady state turbulent convective heat transfer in a tube. In the

Table 1 Chemical compositions of SRM combustion products in the combustion chamber. Species

Mole fraction (mol%)

H2 (g) CO (g) HCl (g) H2O (g) CO2 (g) N2 (g) Al2O3 (g)

33.158 25.682 12.729 8.583 8.360 7.204 7.219

Table 2 Other parameters of the test SRM. Combustion law Stagnation temperature (K) Propellant density (Kg/m3) Propellant Young's modulus (MPa) Propellant Poisson's ratio Propellant ignition temperature (K) Propellant thermal conductivity (cal/cm s K)

r_ ¼ apn 3540.0 1780.0 2.1 0.498 760.0 0:9  10  3

Fig. 6. Time–history curves of the head end and aft end pressure.

motor ignition transient, the internal flow is strong unsteady and highly turbulent. The convective heat transfer coefficient of the internal flow is larger than that of a steady flow. Another factor is the formula for evaluating the radiative heat transfer [25], in which the heat transfer coefficient Csf is an empirical constant and is usually set to 0.25. This value is obtained from time-averaged static firing data and may not be particularly suitable for ignition transient simulation. For the test motor involved, the motor head end pressure increases from approximately 0.1 MPa to 23.6 MPa and eventually reaches an equilibrium state with a time interval of approximately 18.0 ms. The gradient of chamber pressure variation is steep, and no ignition peak or pressure wave is found in both the experimental curve and the numerical curve. The pressure difference between the head end and the aft end of the test motor is remarkable. The simulated data curves plotted in Fig. 6 show that the motor head end pressure is unequal with the pressure from the aft end. The axial pressure difference becomes evident at 4.9 ms, and then increases with time for up to 3.1 MPa at 18.0 ms. This phenomenon is generally known as pressure drop and has been reported and studied in research of the Space Shuttle Reusable Solid Rocket Motor [37]. The axial pressure drop is caused by a complicated combination of geometrical and physical features, including a large length to bore diameter ratio, low port to throat area ratio, mass addition along the bore, and flow restrictions. In this study, the length-tobore diameter ratio L/D of the test motor is 11.13, and the port-to-throat area ratio is 1.71. Thus, sufficient surface mass flow causes choking and meanwhile the aft end of the motor experiences the high Mach numbers and lower static pressure, which creates an axial static pressure drop.

Q. Li et al. / Acta Astronautica 110 (2015) 180–190

187

Fig. 7. Temperature distribution in the motor chamber at different times (K).

The evolution of the internal flow field temperature is of great importance in the ignition transient of SRMs because of its fully decisive role in determining the heat transfer from the hot combustion products to the propellant solid, the ignition of the propellant, and the flame propagation along the propellant burning surface. The time interval of motor ignition transient of the smalland moderate-scale SRMs is a small fraction of one second, during which the motor chamber pressure rises from an initial state to an equilibrium state. This time interval is conventionally divided into three successive stages for analytical purpose, namely, the ignition time lag, flamespreading interval, and chamber-filling interval [26]. Fig. 7 shows the temperature distributions in the motor chamber in the ignition transient. Confined by the propellant geometry, the igniter jet flows downward and fills most of the combustion chamber at 3.0 ms. The propellant burning surface is first ignited near the middle of the chamber bore at 3.42 ms, which marks the end of the ignition time lag. Then, the flame spreads to both the upward and downward directions of the chamber separately, and over 64% and 92% of the burning surface is ignited at 5.3 ms and 8.3 ms, respectively. The flamespreading interval finishes at 9.4 ms when the burning surface is fully ignited, and the chamber-filling interval ends at 18.0 ms when the chamber pressure achieves a quasi-steady state and the whole ignition transient completes. Compared with the ignition transient data reported in previous works [1,2,5,6], the ignition transient of the motor in this study is quite short. This finding is caused by the following reasons. On the one hand, the motor is small in scale. On the other hand, both the mass flow rate of the igniter jet, which is up to 490.0 g/s, and the propellant burning rate are comparatively high. The propellant deformation occurs primarily in the motor axial direction. In the ignition transient of an SRM, the deformation of the solid propellant under rapid pressurization is one of the most concerned key factors to determine the security of the motor operation. First, the configuration of the grain with deformation has a strong effect on the actual propellant burning surface area, which determines exactly the total mass flow rate of the combustion products and the corresponding internal ballistics of the motor. Second, SRM consists of many parts that are made of different materials. In the ignition transient, all parts of SRM undergo deformations of different degrees under the rapid pressurization rate. These deformations must meet the conformance requirements for a successful operation. For

Fig. 8. Displacement on the propellant burning surface at different times. (a) Axial displacement and (b) radial displacement.

that it is extremely different to measure the propellant deformation in the hot firing of SRMs, the numerical results play a crucial role in assisting the designer in making exact evaluations and estimations. Fig. 8 shows the propellant burning surface axial and radial displacements at different times. The axial displacement curves show that no deformation is noticeable before 5.0 ms, at which the maximum deformation is less than 0.41 mm and can be neglected. Then, the propellant geometry undergoes significant deformation, and at 16.0 ms, the axial displacement obtains its maximal negative value. This value is approximately  4.01 mm, which indicates that near this point, the propellant moves in the

188

Q. Li et al. / Acta Astronautica 110 (2015) 180–190

negative direction of motor axis. The grain obtains its maximal positive axial displacement at 18.0 ms. This value is approximately 9.62 mm at a distance of 42.60 mm from the motor head end. At this point, the axial displacement is positive, thereby indicating that the displacement occurs in motor axial direction. During the whole ignition transient, both the quantity and the variation of the radial displacement along the motor axis are small compared with those of the axial displacement. This finding indicates that the axial deformation of the motor involved in this study is predominant. The propellant deformation is essentially attributed to the axial pressure drop. Fig. 8 illustrates that the propellant axial displacements at different times increase at a distance of approximately 400.02 mm from the motor head end and then vary in diverse ways with increasing axial distance. However, the radial displacement curves vary in a similar manner by first increasing from 0.0 mm to approximately 1.71 mm at a distance of approximately 69.92 mm from the motor head end, then decreasing slowly with the increase of axial distance, and even becoming negative. The negative radial displacements indicate that the central chamber bore is locally narrowed. This type of distribution is caused by two reasons. The first reason is constraint. The head end of the propellant is spatially fixed in computation with zero displacement. Thus, the computed axial and radial displacements are absolute zero at the location where x¼0.0 mm. The second reason is related to the force applied on the propellant burning surface. Two kinds of forces are applied on the propellant burning surface, namely, pressure force and shear force. Shear force is fairly small compared with pressure force and is neglected as mentioned earlier. The simulation results show that in the chamber of the motor, the axial pressure drop is remarkably great. The pressure force near the motor head end is remarkably greater than that near the aft end and this pressure force imbalance causes the propellant to be deformed in the axial direction and even to narrow the chamber central bore. In turn, the narrowing of the chamber central bore increases the gas pressure force imbalance by increasing the local Mach numbers and by lowering the corresponding static pressure near the motor aft end. The narrowing of the chamber central bore is clearly reflected by the negative radial displacement near the motor aft end. The strain of the solid propellant grain under rapid pressurization is a basic element for systematic and comprehensive assessment of the reliability and safety of the SRM operation. The reason is that the strain of the propellant is closely related to the initiation and nucleation of the small-scale cracks in propellant geometry. The growth and propagation of these cracks have a potential effect on the overall motor performances. Moreover, if the growth and propagation of these cracks become unmanageable, largescale cracks in the grain that degrade the adaptability of the motor may form. These large-scale cracks may even lead to motor operation failure. Fig. 9 shows the axial normal strain curves on the propellant burning surface at different times. The curves show that the propellant grain is locally stretched or compressed alternately with the increased strength along the motor axial direction. This particular type of axial normal strain distribution is related directly to the pressure force

Fig. 9. Axial normal strain curves on the propellant burning surface at different times.

applied on the propellant burning surface. In this study, the pressurization rate of the motor is rather high, and the continuous and rapid rise of the pressure force applied on the pressure burning surface creates a series of mechanical wave propagating along the motor axial direction. The propellant is compressed at the crest and stretched at the trough of the mechanical wave. Moreover, the strength of the axial pressure imbalance is great and increases with time. Thus, the propellant axial deformation amplitude increases with time and so does the axial normal strain. A close-up view of the propellant axial displacement and normal strain curves shows that wave inversion occurs. For the test motor studies, the inner surface of the propellant conical section is an interface across which the media density variations are comparatively large. When the mechanical wave traveling in the propellant approaches and reaches the propellant conical section surface, a small fraction of the energy possessed by the mechanical wave is transmitted to the surrounding combustion products and the majority of the energy is reflected and the wave is inverted upon reflecting off. This process is generally referred as inversion. In inversion, a positive displacement reflects off the interface and returns with a negative displacement or vice versa. Thus, the positive displacement at 10.0 ms changes to negative displacement at 12.0 ms upon reaching and reflecting off the interface. This phenomenon occurs several times from 12.0 ms to 18.0 ms. In the calculation, one line, which is parallel to the motor axis, is drawn on the fluid–solid coupling interface. The pressure data from the fluid subdomain and the interpolated data from the solid subdomain are both extracted at 16.1 ms and plotted in Fig. 10 with the relative error e. The relative error e is calculated as follows:

org

p  pint

e ¼

porg

where porg and pint are the original data from the fluid side and interpolated data from the solid side, respectively. The pressure and the interpolated data match well. In the motor chamber, the pressure varies from 23.78 MPa at the head end to 19.89 Mpa at the aft end. In addition, the maximal relative error is smaller than 1.0%. These results indicate that the relative errors of the interpolation scheme presented in

Q. Li et al. / Acta Astronautica 110 (2015) 180–190

Fig. 10. Distributions of the flow parameter and interpolated data.

this paper are allowable for engineering applications. The calculated relative error is comparable with that reported in [11] where the same interpolation scheme was employed for the three-dimensional volume mesh data interpolation. This finding indicates that the presented interpolation scheme is robust and adaptable, although it is quite simple. 6. Conclusions In this paper, an integrated framework is presented for the coupled simulations of the multiphysics process involving propellant heating and ignition, combustion chamber strong unsteady fluid flow, and the response of the solid propellant in the ignition transient of SRM. The results show that the interpolation scheme employed to interpolate data through the fluid–solid coupling interface provides satisfactory accuracy for engineering applications. The integrated framework allows us to perform coupled simulations of the propellant ignition, strong unsteady internal fluid flow, and propellant mechanical response in SRMs with satisfactory stability and efficiency and presents a reliable and accurate solution to complex multi-physics problems. Moreover, the outstanding progresses in the solid propellant combustion made in the last ten years have strongly promoted our knowledge in this area, the condensed-phase energy conservation, gas-phase species diffusion and conservation, gas–solid interfacial boundary conditions, ambient oxidizer species diffusion, and the detailed heterogeneous reactions at the propellant surface reaction zone have all been comprehensively studied and formulated [38], all these are the improvements and contents that will be made in our future studies. References [1] G.D. Luke, M.A. Eagar, H.A. Dwyer, Ignition Transient Model for Large Aspect Ratio Solid Rocket Motor, AIAA paper 96-3273. [2] W.A. Johnston, Solid rocket motor internal flow during ignition, J. Propuls. Power 11 (1995) 489–496. [3] W.H. Andersen, Model of transient ignition of self-sustained burning, Combust. Sci. Technol. 5 (1972) 75–81. [4] E.W. Price, H.H. Bradley Jr, G.L. Dehority, M.M. Ibiricu, Theory of ignition of solid propellants, AIAA J. 4 (1966) 1153–1181.

189

[5] W.C. Shiang, Dynamic analysis of solid propellant grains subjected to ignition pressurization loading, J. Sound Vib. 268 (2003) 465–483. [6] W.C. Shiang, Studies of poisson's ratio variation for solid propellant grains under ignition pressure loading, Int. J. Press. Vessels Pip. 80 (2003) 871–877. [7] N.N. Smirnov, I.D. Dimitrienko, Convective combustion of porous compressible propellants, Combust. Flame 89 (1992) 260–270. [8] E. Lefrançois, G. Dhatt, D. Vandromme, Fluid–structural interaction with application to rocket engines, Int. J. Numer. Methods Fluids (1999) 865–895. [9] G.P. Guruswamy, Development and applications of a large scale fluids–structures simulation process on clusters, Comput. Fluids 36 (2007) 530–539. [10] L. Garelli, R.R. Paz, A. Mario, Fluid–structure interaction study of the start-up of a rocket engine nozzle, Comput. Fluids 39 (2010) 1208–1218. [11] Q. Li, G.Q. He, P.J. Liu, J. Li, Coupled simulation of fluid flow and propellant burning surface regression in a solid rocket motor, Comput. Fluids 93 (2014) 146–152. [12] W.B. Merrill, Karthik Puduppakkam, Piyush Thakre, Vigor Yang, Modeling of combustion and ignition of solid-propellant ingredients, Prog. Energy Combust. Sci. 33 (2007) 497–551. [13] E.W. Price, J.H. Bradley, G.L. Dehority, M.M. Ibiricu, Theory of ignition of solid propellants, AIAA J. 4 (1966) 1153–1181. [14] A.K. Kulkarni, M. Kumar, K.K. Kuo, Review of solid-propellant ignition studies, AIAA J. 20 (1982) 243–244. [15] S. Rugonyi, K.J. Bathe, On the analysis of fully-coupled fluid flows with structural interactions – a coupling and condensation procedure, Int. J. Comput. Civil Struct. Eng. 1 (2000) 29–41. [16] K.J. Bathe, H. Zhang, M.H. Wang, Finite element analysis of incompressible and compressible fluid flows with free surfaces and structural interactions, Comput. Struct. 56 (1995) 193–213. [17] O.C. Zienkiewicz, R.L. Taylor, The Finite Element Methods, Butterworth-Heinemann, Burlington, 2000. [18] A. Linan, A. Crespo, An asymptotic analysis of radiant and hypergolic heterogeneous ignition of solid propellants, Combust. Sci. Technol. 6 (1972) 223–232. [19] T. Niioka, Heterogeneous ignition of a solid fuel in a hot stagnationpoint flow, Combust. Sci. Technol. 18 (1978) 207–215. [20] V.N. Vilyunov, V.E. Zarko, Ignition of Solids, Elsevier Science Publishers, New York, 1989. [21] Zel' dovich, B. Ya, On the theory of propellant combustion, Zhurnal Eksp. I Theoret. Fiz. 12 (1942) 11–12. [22] Zel' dovich, B. Ya, On a burning rate under nonsteady pressure, Zhurnal Prikl. Mekhaniki I Tekh. Fiz. 3 (1964) 752–758. [23] B.V. Novozhilov, Theory of nonsteady burning and combustion stability of solid propellants by the Zeldovich-Novozhilov methodNonsteady Burning and Combustion Stability of Solid Propellants, . [24] W.E. William, Modeling the Unstable Combustion of Solid Propellants with Detailed Chemistry, . [25] W.A.J. Solid Rocket, Motor internal flow during ignition, J. Propuls. Power 11 (1995) 489–496. [26] N.N. Smirnov, Combustion of porous dispersing fuels, Combust. Explos. Shock Waves 27 (1991) 52–58. [27] N.N. Smirnov, I.D. Dimitrienko, Investigation of a convective combustion in compressible solid propellants with channels, Combust., Explos. Shock Waves 26 (1990) 14–22. [28] G.P. Sutton, B. Biblarz, Rocket Propulsion Elements, JOHN WILEY & SONS, New York, 2000. [29] X. Jiao, M.T. Heath, Common-refinement-based data transfer between nonmatching meshes in multiphysics simulations, Int. J. Numer. Methods Eng. 61 (2004) 2402–2427. [30] X. Jiao, M.T. Heath, Overlaying surface meshes, Part I: algorithms, Int. J. Comput. Geom. Appl. 14 (2004) 379–402. [31] S. Piperno, C. Farhat, B. Larrouturou, Partitioned procedures for the transient solution of coupled aroelastic problems part I: model problem, theory and two-dimensional application, Comput. Methods Appl. Mech. Eng. 124 (1995) 79–112. [32] F.J. Blom, A monolithical fluid–structure interaction algorithm applied to the piston problem, Comput. Methods Appl. Mech. Eng. 167 (1998) 369–391. [33] X. Jiao, G. Zheng, P.A. Alexander, M.T. Campbell, et al., A system integration framework for coupled multiphysics simulations, Eng. Comput. 22 (2006) 293–309. [34] J. Giordano, G. Jourdan, Y. Burtschell, et al., Shock wave impacts on deforming panel, an application of fluid–structure interaction, Shock Waves 14 (2005) 103–110. [35] S. Gordon, B.J. McBride, Computer Program for Calculation of Complex Chemical Equilibrium Compositions and Applications

190

Q. Li et al. / Acta Astronautica 110 (2015) 180–190

I: Analysis, NASA Lewis Research Center: NASA Reference Publication, 1994, 25–32 (report number NASA-RP-1311). [36] I.G. Assovskii, A.G. Merzhanov, Validity of experimental and theoretical modeling of combustion of energetic materials, Combust., Explos. Shock Waves 49 (2013) 264–272.

[37] A.L. Brian, Internal Flow Analysis of Large Scale Solid Rocket Motors, AIAA Paper 2000-3803. [38] I.G. Assovskiy, Physics of Combustion and Interior Ballistics, Nauka, Moscow, 2005.