Thermo-Structural Response Caused by Structure Gap and Gap Design for Solid Rocket Motor Nozzles

The thermo-structural response of solid rocket motor nozzles is widely investigated in the design of modern rockets, and many factors related to the material properties have been considered. However, little work has been done to evaluate the effects of structure gaps on the generation of flame leaks. In this paper, a numerical simulation was performed by the finite element method to study the thermo-structural response of a typical nozzle with consideration of the structure gap. Initial boundary conditions for thermo-structural simulation were defined by a quasi-1D model, and then coupled simulations of different gap size matching modes were conducted. It was found that frictional interface treatment could efficiently reduce the stress level. Based on the defined flame leak criteria, gap size optimization was carried out, and the best gap matching mode was determined for designing the nozzle. Testing experiment indicated that the simulation results from the proposed method agreed well with the experimental results. It is believed that the simulation method is effective for investigating thermo-structural responses, as well as designing proper gaps for solid rocket motor nozzles.


Introduction
The solid rocket motor (SRM) nozzle always acts as an energy transformation device where the chemical energy of propellant turns into the kinetic energy of gas with extremely high temperature and pressure [1][2][3]. Normally, the temperatures can be 3000 K or higher and the pressure may reach about 10 MPa or higher for some special motor devices, especially those equipped with modern propellants that have better performance [4,5]. What is more, the temperatures and pressures of the gas are becoming higher and higher due to the increasing performance demands of SRMs. Overall, the nozzle is always working under extremely harsh circumstances [6]. Hence, the structure integrity is of great importance for the normal operation of the nozzle. In view of this, many researchers around the world have focused on several aspects of this subject.
Kumar et al. [7] investigated the thermo-structural response of a nozzle made of composite under thermal and mechanical loads, and the results showed that thermal stress played a more important role these results, some designed gap sizes were given for correction, which provided a reference for future nozzle design. Sun et al. [32] simulated the fully-coupled thermo-structural response of a nozzle with adhesive failure, and the results showed that treating the interface as contact status could reduce the stress intensity. Goyal et al. [33] simplified the traditional 3D method, and found that a 2D model could achieve acceptable results. This method can be applied to the simulation of several parts of the nozzle and the nozzle insert, as well as the thermo-structural simulation with consideration of elasticity modulus reduction caused by ablation. Wang [34] took an SRM nozzle as example to simulate the transient heat and stress fields, and showed that the connecting and bonding status should be carefully considered in structure design due to the significant effects of thermal expansion coefficients. Zheng [35] established a simulation method for the nozzle insert of C/C composite material, with an overall consideration of contact thermal resistance, ablation moving boundary, interface debonding and so on. Based on the method, interface debonding was considered as the key factor that led to nozzle failure. Zheng et al. [36,37] proposed an integrative computation method in Fluent for an SRM nozzle based on fluid-structure interaction, and realized the automatic identification of contact boundaries. Teng et al. [38] analyzed the transient temperature field of an axisymmetric nozzle with consideration of thermo-structural, geometry characteristics, material characteristics and boundary symmetry, and employed a special GAP element in MSC.NASTRAN to simulate the contact status of different materials, which provided useful basis for nozzle thermal-structural design. Li [39] regarded the interface between the nozzle insert and the back surface as a contact one because the interface would fail due to the high temperature and shear stress. Liu [40] confirmed the friction coefficient between C/C material and asbestos phenolic material, and provided an accurate friction coefficient for the contact simulation of thermo-structural response.
To sum up, a large amount of work has been done for investigating the thermo-structural response concerning the contact status. However, little attention has been paid to the gap design, especially the specific size of each gap. This paper mainly focuses on the contact simulation, and meanwhile a gap size optimization is included. Furthermore, some experimental results are given to evaluate the design of the nozzle. Finally, some conclusions are drawn in the end.

Geometry Structure
This paper takes a typical nozzle as example to conduct a thermo-structural analysis based on the contact status and the structure gap. Figure 1 shows the nozzle sketch. The nozzle usually has three layers, including ablation layer, heat insulating layer and support layer. The ablation layer, i.e., the convergent and divergent regions, is made of carbon-phenolic-resin, the heat insulating layer is of silica-phenolic-resin, and the support layer is of structural steel. The throat insert is made of carbon/carbon composite material. Table 1 gives the nozzle dimensions. nozzle with adhesive failure, and the results showed that treating the interface as contact status could reduce the stress intensity. Goyal et al. [33] simplified the traditional 3D method, and found that a 2D model could achieve acceptable results. This method can be applied to the simulation of several parts of the nozzle and the nozzle insert, as well as the thermo-structural simulation with consideration of elasticity modulus reduction caused by ablation. Wang [34] took an SRM nozzle as example to simulate the transient heat and stress fields, and showed that the connecting and bonding status should be carefully considered in structure design due to the significant effects of thermal expansion coefficients. Zheng [35] established a simulation method for the nozzle insert of C/C composite material, with an overall consideration of contact thermal resistance, ablation moving boundary, interface debonding and so on. Based on the method, interface debonding was considered as the key factor that led to nozzle failure. Zheng et al. [36,37] proposed an integrative computation method in Fluent for an SRM nozzle based on fluid-structure interaction, and realized the automatic identification of contact boundaries. Teng et al. [38] analyzed the transient temperature field of an axisymmetric nozzle with consideration of thermo-structural, geometry characteristics, material characteristics and boundary symmetry, and employed a special GAP element in MSC.NASTRAN to simulate the contact status of different materials, which provided useful basis for nozzle thermalstructural design. Li [39] regarded the interface between the nozzle insert and the back surface as a contact one because the interface would fail due to the high temperature and shear stress. Liu [40] confirmed the friction coefficient between C/C material and asbestos phenolic material, and provided an accurate friction coefficient for the contact simulation of thermo-structural response.
To sum up, a large amount of work has been done for investigating the thermo-structural response concerning the contact status. However, little attention has been paid to the gap design, especially the specific size of each gap. This paper mainly focuses on the contact simulation, and meanwhile a gap size optimization is included. Furthermore, some experimental results are given to evaluate the design of the nozzle. Finally, some conclusions are drawn in the end.

Geometry Structure
This paper takes a typical nozzle as example to conduct a thermo-structural analysis based on the contact status and the structure gap. Figure 1 shows the nozzle sketch. The nozzle usually has three layers, including ablation layer, heat insulating layer and support layer. The ablation layer, i.e., the convergent and divergent regions, is made of carbon-phenolic-resin, the heat insulating layer is of silica-phenolic-resin, and the support layer is of structural steel. The throat insert is made of carbon/carbon composite material. Table 1 gives the nozzle dimensions.

Material Model
Great efforts have been made to choose proper materials for the nozzle. Here, all properties of the material models are temperature-dependent. The specific properties are shown in Tables 2-4, where X represents the radial direction, Y represents the axial direction and Z represents the circumferential direction. By the way, several variables and their nomenclatures are used for saving space. The specific variables are as follows: Coefficient of thermal expansion (K´1) K Thermal conductivity (W/(m¨K)) C Specific heat (J/(kg¨K))  Table 3. Material properties of carbon-phenolic-resin [41].

Finite Element Model
It is noticed that the attached nozzle model is axially symmetrical. Normally, the axisymmetric finite element model can not only reduce the number of meshes, but also eliminate to a large extent the errors caused by mesh dissymmetry. Therefore, a 2D axisymmetric model is used in this paper for both the thermo-structural simulation and the flow simulation.

Meshing
A properly sized mesh can generate more accurate results for both the thermo-structural simulation and the flow simulation. Figure 2 shows the mesh generation result of the thermo-structural simulation.

Finite Element Model
It is noticed that the attached nozzle model is axially symmetrical. Normally, the axisymmetric finite element model can not only reduce the number of meshes, but also eliminate to a large extent the errors caused by mesh dissymmetry. Therefore, a 2D axisymmetric model is used in this paper for both the thermo-structural simulation and the flow simulation.

Meshing
A properly sized mesh can generate more accurate results for both the thermo-structural simulation and the flow simulation. Figure 2 shows the mesh generation result of the thermostructural simulation.

ANSYS Parameter Design Language in ANSYS Workbench
According to previous work [42,43], reliable results can be obtained with ANSYS. In this work, the simulation is conducted on ANSYS Workbench that does not have a function of directly coupled thermo-structural simulation by now, so that necessary ANSYS Parameter Design Language (APDL) codes have been developed by the authors ourselves. The thermo-structural response is simulated in the Mechanics of Transient Structure, which can only simulate the structure response. Firstly, an element-choosing procedure is conducted by the ET command for all the parts of the nozzle. Secondly, KEYOPT and RMODIF commands are executed for the contact definition in order to simulate the heat transfer. Considering the radiation, MP command is used to define the emission of the materials. Finally, the thermal boundary conditions are defined and applied by *DIM, *TREAD and SF commands.

Boundary Conditions
For the thermo-structural simulation, both thermal and structural boundary conditions are required. The specific settings for this simulation are as below.

ANSYS Parameter Design Language in ANSYS Workbench
According to previous work [42,43], reliable results can be obtained with ANSYS. In this work, the simulation is conducted on ANSYS Workbench that does not have a function of directly coupled thermo-structural simulation by now, so that necessary ANSYS Parameter Design Language (APDL) codes have been developed by the authors ourselves. The thermo-structural response is simulated in the Mechanics of Transient Structure, which can only simulate the structure response. Firstly, an element-choosing procedure is conducted by the ET command for all the parts of the nozzle. Secondly, KEYOPT and RMODIF commands are executed for the contact definition in order to simulate the heat transfer. Considering the radiation, MP command is used to define the emission of the materials. Finally, the thermal boundary conditions are defined and applied by *DIM, *TREAD and SF commands.

Boundary Conditions
For the thermo-structural simulation, both thermal and structural boundary conditions are required. The specific settings for this simulation are as below.

‚ Thermal Boundary Conditions
There exist two kinds of thermal boundary conditions for this simulation, i.e., the free convection heat transfer between the outermost structural steel and the air, and the forced convection from the inner flow. Many studies [44,45] simplified that it was adiabatic between the outmost and the atmosphere. Deng [17] accepted the probable coefficient scope of 1-10 W/(m 2¨K ) and discussed the effects of the coefficient variation on the simulation. As a result, an average value of 5 W/(m 2¨K ) from Deng's research is adopted in this paper. Meanwhile, this value is also an empirical free convection coefficient for one kind of default setting called Simplified Case in ANSYS Workbench.
For the forced convection, the coefficient can be calculated by the Bartz formulation [46]. When the hot gas enters the nozzle, the velocity is accelerated so fast that the turbulence effect is increased sharply, which enhances the convection heat transfer between the flow and the wall. For this special flow form, Bartz proposed an experimental formulation to predict the coefficient, and the predicted results show reasonable agreement with the experimental data except for a small section before the throat: where h is the forced convection coefficient; C = 0.026 for subsonic flow and C = 0.023 for supersonic flow; d t is the diameter of the throat; c p is the constant-volume specific heat of the gas; µ is the coefficient of viscosity of the gas; Pr is the Prandtl number; . m is the mass flow rate; A t is the area of the throat; r c is the curvature radius; A is the area of the local cross profile; σ is a modified parameter considering the vibration of the boundary layer; T w is the temperature of the wall; T 0 is the temperature of the total temperature of the gas; k is the ratio of specific heat; and Ma is the local Mach number.
To calculate heat transfer from the inner flow to the wall, the recovery temperature T r of the gas should be provided as the temperature of the gas that directly contacts the nozzle and the value may be a little bigger than the total temperature. The recovery temperature can be calculated as given equation: where T and Ma are the gas static temperature and Mach number of the local cross profile. For the Prandtl number Pr and coefficient of viscosity of the gas µ, they can be calculated using the equations below [46]: where M is the average molecular weight of the gas and k = 1.2 for hot gas in SRM.

Structure Boundary Conditions
In general, structure boundary conditions include fixed constraints, displacement, pressure and so on. For a nozzle at work, the fixed and displacement constrains are the main boundary conditions, and the inner flow pressure cannot be ignored either. The nozzle is often stuck to the chamber by bolts. Meanwhile, some devices are often installed on the outside of the nozzle in order to save space. To simplify the simulation procedure, firstly, it is assumed that the outmost structural steel can move along the axial direction but not the radial direction; secondly, the radial direction outlines at the head of the nozzle are constrained to move along the radial direction but not the axial direction; thirdly, the outlines at the axis movement are constrained to move along the axial direction but not the radial direction; lastly, the radial direction outlines at the end of the nozzle are also constrained to move along the radial direction but not the axial direction.

Initial Gap Definition
In this section, the initial gap size is defined by experience [29]. For the non-metal and metal fitting surface, the optimum gap size interval may be 0.1 to 0.2 mm, and the maximum shall be not larger than 0.25 mm. For the nonmetal and non-metal fitting surface, the designed gap size can be smaller. Larger gap size may bring about bad splicing results and cause damage to the structure such as the flame leaks mentioned in the Introduction section. For the chosen nozzle, two possible passageways are illustrated in Figure 3. In order to be in accordance with the experiment introduced in the Experimental Results section, the simulation will only focus on the left passageway, and the initial gap sizes are defined in Table 5.

Assumptions
Some reasonable assumptions are made as follows to simplify the mathematical model and accelerate the convergent course: 1. The complicated phenomena of erosion and pyrolysis behavior of the erosion and heat insulation materials are neglected; 2. The connecting part of the nozzle and the chamber has no heat exchange, and other parts experience free convection; 3. The contact thermal resistance is assumed to be a constant of 10 −3 m 2 •K/W [35], which ignores the variation of thermal resistance caused by environmental changes such as contact status, temperature and so on; 4. The contact thermal resistance is totally ignored for the bonding simulation; 5. During the frictional contact movement, energy dissipation caused by friction is totally transformed into heat, and the contact and the target respectively absorbs 50% of the heat; 6. There exists radiation heat transfer between the contact and the target. The Stefan-Boltzmann constant is 5.67 × 10 −8 W/(m 2 •K 4 ) and the heat emissivity for each material is 0.9, which is the same as that for C/C composite material [47]. Because the radiation from the hot gas to the nozzle wall is very low relative to the forced heat convection [31], the internal and external radiation heat transfer are neglected.

Assumptions
Some reasonable assumptions are made as follows to simplify the mathematical model and accelerate the convergent course: 1.
The complicated phenomena of erosion and pyrolysis behavior of the erosion and heat insulation materials are neglected; 2.
The connecting part of the nozzle and the chamber has no heat exchange, and other parts experience free convection; 3.
The contact thermal resistance is assumed to be a constant of 10´3 m 2¨K /W [35], which ignores the variation of thermal resistance caused by environmental changes such as contact status, temperature and so on; 4.
The contact thermal resistance is totally ignored for the bonding simulation; 5.
During the frictional contact movement, energy dissipation caused by friction is totally transformed into heat, and the contact and the target respectively absorbs 50% of the heat; Energies 2016, 9, 430 8 of 21 6. There exists radiation heat transfer between the contact and the target. The Stefan-Boltzmann constant is 5.67ˆ10´8 W/(m 2¨K4 ) and the heat emissivity for each material is 0.9, which is the same as that for C/C composite material [47]. Because the radiation from the hot gas to the nozzle wall is very low relative to the forced heat convection [31], the internal and external radiation heat transfer are neglected.

Flow Results Based on Quasi-1D Model
Theoretically, the flow in a nozzle can be generated in milliseconds. In comparison with the whole working time, this period of time can be ignored. In particular, this paper emphasizes the thermo-structural response regarding contact status, so that the generation of the flow is not very important. As a result, a Quasi-1D model is established to acquire the fluid properties along the nozzle axis and provide boundary conditions for the nonlinear structure simulation.
For one dimension isentropic flow, the relationship of temperature, pressure and Mach number at any section along the nozzle axis can be given by: where T 0 is the chamber temperature, P 0 is the chamber pressure, A t is the nozzle throat area, A is the local area along the nozzle axis, and k is the specific heat ratio of the gas. For the chosen nozzle, the chamber temperature T 0 is 3500 K, and the chamber pressure P 0 is 8 MPa. A t can be obtained from Table 1. A can be figured out during modeling of the nozzle sketch and k is determined by the gas property. Figure 4 gives the temperature and Mach number curves obtained from the quasi-1D model. And x is the current location along the axis of the nozzle.
where T0 is the chamber temperature, P0 is the chamber pressure, At is the nozzle throat area, A is the local area along the nozzle axis, and k is the specific heat ratio of the gas. For the chosen nozzle, the chamber temperature T0 is 3500 K, and the chamber pressure P0 is 8 MPa. At can be obtained from Table 1. A can be figured out during modeling of the nozzle sketch and k is determined by the gas property. Figure 4 gives the temperature and Mach number curves obtained from the quasi-1D model. And x is the current location along the axis of the nozzle.

Forced Convection Coefficient Processing
For thermo-structural simulation, the convection coefficient along the wall and the recovery temperature should be provided [29]. Based on Equations (1)-(3), the coefficient and recovery temperature can be calculated as shown in Figure 5. As can be seen, the maximum value of forced convection coefficient is about 18,000 W/(m 2 •K) and the minimum is about 1800 W/(m 2 •K). At the entrance of the nozzle, due to the turbulent layer effect, the recovery temperature Tr is a little higher than the total temperature T0. This phenomenon will disappear soon. It is noticed that the variation curve of forced convection coefficient presents a peak at the throat, which represents a tough ordeal for thermal protection design.

Forced Convection Coefficient Processing
For thermo-structural simulation, the convection coefficient along the wall and the recovery temperature should be provided [29]. Based on Equations (1)-(3), the coefficient and recovery temperature can be calculated as shown in Figure 5. As can be seen, the maximum value of forced convection coefficient is about 18,000 W/(m 2¨K ) and the minimum is about 1800 W/(m 2¨K ). At the entrance of the nozzle, due to the turbulent layer effect, the recovery temperature T r is a little higher Energies 2016, 9, 430 9 of 21 than the total temperature T 0 . This phenomenon will disappear soon. It is noticed that the variation curve of forced convection coefficient presents a peak at the throat, which represents a tough ordeal for thermal protection design.

Forced Convection Coefficient Processing
For thermo-structural simulation, the convection coefficient along the wall and the recovery temperature should be provided [29]. Based on Equations (1)-(3), the coefficient and recovery temperature can be calculated as shown in Figure 5. As can be seen, the maximum value of forced convection coefficient is about 18,000 W/(m 2 •K) and the minimum is about 1800 W/(m 2 •K). At the entrance of the nozzle, due to the turbulent layer effect, the recovery temperature Tr is a little higher than the total temperature T0. This phenomenon will disappear soon. It is noticed that the variation curve of forced convection coefficient presents a peak at the throat, which represents a tough ordeal for thermal protection design.

Coupled Nonlinear Contact Simulation
Thermo-structural simulation can not only verify the nozzle structure integrity, but also ascertain the weakest parts or points with stress concentration. In this paper, the finite element method is employed to simulate the temperature field in the thermo-structural modules, and the

Coupled Nonlinear Contact Simulation
Thermo-structural simulation can not only verify the nozzle structure integrity, but also ascertain the weakest parts or points with stress concentration. In this paper, the finite element method is employed to simulate the temperature field in the thermo-structural modules, and the 2-D 8-Node Coupled-Field Solid element (Plane223) is chosen. For other additional information readers can refer to the ANSYS Mechanical APDL Theory [48]. The calculation contains only one simulation with all necessary freedom elements, and the coupled analysis is carried out on all element matrixes or load vectors of the physical variables. The governing equations are as below, and the meaning of each symbol can be found in [48]: In this section, the left passageway of flame leak in Figure 3 (1Ñ2Ñ3Ñ4Ñ5) is analyzed and the contact status takes the initial gap sizes into account. All the other interfaces are treated as bonded ones. Further, two types of simulations are delivered.

Discussion on Erosion and Pyrolysis Behavior
As phenolic resin based materials, when heated, carbon-phenolic and silica-phenolic start to decompose at about 300˝C. Along with decomposition, vaporized SiO 2 , CO 2 and other kinds of gases will also be generated. When the temperature reaches about 800˝C and higher, the porosity factor will be constant, indicating the end of thermal decomposition [29].
Besides, carbon residue will form a carbon layer with porous media characteristics. The decomposition gas would also flow among the small holes and enter into the gaps (1Ñ2Ñ3Ñ4Ñ5 or 6Ñ7Ñ8Ñ9Ñ10). Except for thermal expansion stress and extrusion stress between different parts, the decomposition gas can lead to pressure load and cause stress responses on all the sides of gaps labeled 1Ñ10.
For silica/phenolic composites, Shi [49] investigated the ablation mechanism and thermo-mechanical behavior of silica/phenolic composites. By comparing the temperature distribution and pore pressure distribution of the material, a relationship between temperature and pore pressure can be obtained as shown in Figure 6. As stated in [49], the maximum gas pressure is located at the starting position of decomposition. This matches well with the phenomenon shown in Figure 6 that the maximum pore pressure is about 6.3 MPa when temperature reaches about 800 K.

9
For silica/phenolic composites, Shi [49] investigated the ablation mechanism and thermomechanical behavior of silica/phenolic composites. By comparing the temperature distribution and pore pressure distribution of the material, a relationship between temperature and pore pressure can be obtained as shown in Figure 6. As stated in [49], the maximum gas pressure is located at the starting position of decomposition. This matches well with the phenomenon shown in Figure 6 that the maximum pore pressure is about 6.3 MPa when temperature reaches about 800 K. Zhao [50] did research on ablation behavior of carbon/phenolic composites, and a similar relationship between temperature and pore pressure of carbon/phenolic can be obtained as shown in Figure 7. Like to silica/phenolic composites, the maximum pore pressure also arises at the starting position of decomposition of about 1.2 MPa at about 400 K. Zhao [50] did research on ablation behavior of carbon/phenolic composites, and a similar relationship between temperature and pore pressure of carbon/phenolic can be obtained as shown in Figure 7. Like to silica/phenolic composites, the maximum pore pressure also arises at the starting position of decomposition of about 1.2 MPa at about 400 K. For simulation concerning the contact status, such pressure could not be ignored since it will affect the contact behavior. As a result, to take the erosion and pyrolysis behavior into account as far as possible, such pressures will be imposed on the surfaces around the nozzle insert and back surface labeled 1 → 10. To simplify the pressure imposing process, maximum values of 6.3 MPa for silica/phenolic interfaces and 1.2 MPa for carbon/phenolic interfaces are directly applied to conduct the coupled nonlinear contact simulation. Figure 8 shows the temperature distribution graphs of the two types of interface treating methods. From these temperature distribution graphs, it can hardly figure out how contact status affects the heat transfer process.

Temperature Comparison
In order to analyze the results more clearly, another post-process method is used. Firstly, a new coordinate system is built. The point on the path 1→2→3→4→5 with the smallest y value is defined as the origin of the coordinate, and the coordinate axis value increases along the path 1→2→3→4→5. Meanwhile, d is defined as the current location along the coordinate axis (path 1→2→3→4→5) and For simulation concerning the contact status, such pressure could not be ignored since it will affect the contact behavior. As a result, to take the erosion and pyrolysis behavior into account as far as possible, such pressures will be imposed on the surfaces around the nozzle insert and back surface labeled 1Ñ10. To simplify the pressure imposing process, maximum values of 6.3 MPa for silica/phenolic interfaces and 1.2 MPa for carbon/phenolic interfaces are directly applied to conduct the coupled nonlinear contact simulation.  Figure 8 shows the temperature distribution graphs of the two types of interface treating methods. From these temperature distribution graphs, it can hardly figure out how contact status affects the heat transfer process. In order to analyze the results more clearly, another post-process method is used. Firstly, a new coordinate system is built. The point on the path 1Ñ2Ñ3Ñ4Ñ5 with the smallest y value is defined as the origin of the coordinate, and the coordinate axis value increases along the path 1Ñ2Ñ3Ñ4Ñ5. Meanwhile, d is defined as the current location along the coordinate axis (path 1Ñ2Ñ3Ñ4Ñ5) and D is the total length of the path. Figure 9 shows the temperature distribution curves of the passageway divided by time for the two types of interfaces. From the graphs, as time goes on, the difference between them becomes more and more obvious, and the temperatures of the bonded interface are higher than those of the frictional interface. This means that the thermal contact resistance has a significant effect on preventing heat transfer from the inner to the outer, which can effectively protect the heat insulation layer and support layer. Moreover, with regard to the small part closest to the inner hot gas, the temperature of the frictional interface is higher than that of the bonded interface. At the beginning of heat transfer, the thermal contact resistance plays an important role, and this small part would firstly absorb energy from the inner hot gas. Since the heat flow rate is the largest, the insert has the highest temperature rise, but the thermal contact resistance cuts off heat transfer from the insert to the convergent part. When the contact time reaches a certain point, this effect can be weakened.  Figure 10 shows the Von Mises stress distribution condition. The distribution character shows similar rules. There is also an apparent difference between the bonded and frictional results. The Von Mises stress level of the bonded interface is higher than that of the frictional one, especially at the interfaces around the nozzle insert and back surface. It can be inferred that the contact treatment is conductive to stress release. Figure 11 gives the Von Mises stress variation curves for the passageway. It can be easily found that frictional treatment can reduce the stress level at almost each point. At each time point, Von Mises stress of path 1 is basically at the same level. For path 2, the stress release effect of frictional treatment is the most efficient. For paths 3 to 5, the release effect is moderate. Figure 12 illustrates the axial stress variation curves of the passageway. As can be seen, there also exists considerable stress release effect, and the change tendencies are nearly the same for both interface treating methods. For most time and most section of the interface, the axial stress presents to be compressive stress. Figure 13 shows the radial stress variation curves of the passageway. The stress release effect is not as obvious as that of the axial stress, but the stress level indeed decreases, especially for path 2. The radial stress is also compressive stress for most time and most section of the interface. 12 Figure 12 illustrates the axial stress variation curves of the passageway. As can be seen, there also exists considerable stress release effect, and the change tendencies are nearly the same for both interface treating methods. For most time and most section of the interface, the axial stress presents to be compressive stress. Figure 13 shows the radial stress variation curves of the passageway. The stress release effect is not as obvious as that of the axial stress, but the stress level indeed decreases, especially for path 2. The radial stress is also compressive stress for most time and most section of the interface.     Table 6 shows the average release percentages of the three kinds of stress at different time points and different paths. From the table, it can be easily figured out that for axial stress and radial stress, the average release effect is not obvious as for the Von Mises stress. Examining the original data for Figures 11-13 directly, for most nodes along the paths, the release effect is positive, but for a small number of nodes, the value is negative and the absolute value is much larger. This is the main reason why the release effect shown in Table 6 is not as obvious as shown in Figures 11-13. All results show  Table 6 shows the average release percentages of the three kinds of stress at different time points and different paths. From the table, it can be easily figured out that for axial stress and radial stress, the average release effect is not obvious as for the Von Mises stress. Examining the original data for Figures 11-13 directly, for most nodes along the paths, the release effect is positive, but for a small number of nodes, the value is negative and the absolute value is much larger. This is the main reason why the release effect shown in Table 6 is not as obvious as shown in Figures 11-13. All results show acceptable agreement with [39], namely the simulation considering contact can produce more accurate results. To demonstrate the feasibility of the SRM nozzle design, it is necessary to carry out a structure integrity check for the two different types of interface treatment. It is found that treating the interfaces as bonded conditions would give too much unnecessary margin when choosing materials and designing shape. As a result, the contact nonlinear algorithm should be used when simulating the thermo-structural response of a nozzle.

Gap Size Optimization
Too large a gap size may lead to unexpected damage such as flame leaks [29,30]. In addition, the nozzle insert may deform under the combined loads of inner pressure, force from the neighborhood and its own thermal expansion. If the deformation is too large for other parts to hold the insert, it may fall off the nozzle. In the simulation work, this phenomenon can be demonstrated by the non-convergence of numerical calculations. To avoid flame leaks and prevent inserts from bursting out, a proper gap size should be determined for the given nozzle. In this study, the insert and the back surface are treated as bonded interfaces so that it is impossible for the insert to burst out. Hence, flame leaks will be the only factor to be discussed in this section.

Gap Size Definitions
The initial gap size has been given above. To investigate the effects of other sizes on the nozzle, more gap sizes should be defined. Based on Chen [29], it is assumed that all gap sizes are smaller than 0.25 mm. What is more, the carbon-phenolic-resin will be carbonized under too high temperatures [51]. Consequently, the hot gas from the flow into the structure gap must be appropriately dealt with. In view of this, the gap that directly contacts with the gas should be smaller than the inner ones so as to protect the heat insulation layer as far as possible.
This section mainly focuses on the left passageway 1Ñ2Ñ3Ñ4Ñ5. Nine gap size matching modes can be set up as shown in Table 7, and the same simulation will be carried out on them.

Fire-Puncture Confirmation and Gap Optimization
Firstly, flame leak criteria should be set up. If there is a period of time during which the contact number along the passageway is zero, namely no nodes from both sides are in contact with each other, then the inner gas can flow into the gaps. If the gaps are not closed before the gas reaches the support layer of structural steel, flame leaks will occur. If the time period for which the gaps are not closed is short, the risk of flame leak can be reduced. Figure 14 gives nine graphs corresponding to the gap size modes defined in Table 7. Based on the criteria, there is no flame leak risk for Modes A to C but a little risk for Modes D and E. For Modes G to I, the risk is much larger. Furthermore, it can be concluded that the gap size of interface 1 plays a vital in avoiding flame leaks. This suggests the important tip that workers should control the craft very strictly during manufacture of the nozzle so as to avoid unwanted damage to the nozzle. Finally, considering the stress level and flame leak confirmation, the best gap size matching for the chosen possible passageway can be determined as shown in Table 8.  Finally, considering the stress level and flame leak confirmation, the best gap size matching for the chosen possible passageway can be determined as shown in Table 8.

Experimental Results
In this section, the nozzle was designed based on the simulation considering gap size, and a test experiment was performed to verify the effectiveness of the design. The experiment is successful and the structure integrity of the nozzle is good. Here, the designers chose Mode E for path 1Ñ2Ñ3Ñ4Ñ5. Figure 15 shows the convergent region after the experiment. Under such high temperatures and pressures, the convergent region exhibits a tough appearance due to severe ablation and carbonization as shown in Figure 15a. The highlighted region of (a) is enlarged in Figure 16b. From this figure, a flame leak phenomenon occurs in paths 1 and 2, but no flame leak can be observed in path 3. For path 1, the flame leak is very obvious and the carbonization is as serious as that of the surface that directly contacts with the hot gas. It can be inferred that the hot gas transfers from the inner flow along path 1 and is prevented by path 2 before it reaches path 3. This means that path 2 is closed at one time point. This experimental result agrees well with the simulation result, namely there is a short period when the path is not closed (shown in Figure 14-Mode E).  Figure 16 shows the divergent region after the experiment. In comparison with the convergent region, less ablation and carbonization can be observed. Especially, a hole as highlighted in the red circle goes deep into the second layer at the gap with the insert, which illustrates a flame leak. This implies that the gap size of 0.2 mm chosen by designers may not be suitable. As a result, not only path 1→2→3→4→5, but also path 6→7→8→9→10 needs to be optimized.  Figure 16 shows the divergent region after the experiment. In comparison with the convergent region, less ablation and carbonization can be observed. Especially, a hole as highlighted in the red circle goes deep into the second layer at the gap with the insert, which illustrates a flame leak. This implies that the gap size of 0.2 mm chosen by designers may not be suitable. As a result, not only path 1Ñ2Ñ3Ñ4Ñ5, but also path 6Ñ7Ñ8Ñ9Ñ10 needs to be optimized. Figure 16 shows the divergent region after the experiment. In comparison with the convergent region, less ablation and carbonization can be observed. Especially, a hole as highlighted in the red circle goes deep into the second layer at the gap with the insert, which illustrates a flame leak. This implies that the gap size of 0.2 mm chosen by designers may not be suitable. As a result, not only path 1→2→3→4→5, but also path 6→7→8→9→10 needs to be optimized.

Conclusions
In summary, an axisymmetric model of a nozzle was established to analyze its quasi-1D flow and its thermo-structural response, and its effects on temperature field and stress field were investigated by using a non-linear algorithm when the interfaces between different materials were treated as bonded status and frictional status. Further, a test experiment was performed to verify the numerical results. Some conclusions can be drawn as follows: 1. When interfaces are treated as being in frictional status, the temperature rise can be decreased, which protects the heat insulation layer and the support layer from excessive temperatures. From the perspective of stress analysis, the frictional treatment can release the stress, especially the axial stress; 2. During shape design and thermal protection structure optimization of a nozzle, interfaces should be treated as contact ones to reduce the design margin; Figure 16. Divergent region after the experiment.

Conclusions
In summary, an axisymmetric model of a nozzle was established to analyze its quasi-1D flow and its thermo-structural response, and its effects on temperature field and stress field were investigated by using a non-linear algorithm when the interfaces between different materials were treated as bonded status and frictional status. Further, a test experiment was performed to verify the numerical results. Some conclusions can be drawn as follows: 1. When interfaces are treated as being in frictional status, the temperature rise can be decreased, which protects the heat insulation layer and the support layer from excessive temperatures. From the perspective of stress analysis, the frictional treatment can release the stress, especially the axial stress; 2. During shape design and thermal protection structure optimization of a nozzle, interfaces should be treated as contact ones to reduce the design margin; 3. Based on the defined flame leak criteria and the simulation results, the best gap matching mode can be determined for the nozzle; 4. Experimental results provide great support for verifying the effectiveness of the proposed method.
This article establishes a simplified and fast simulation procedure for verifying nozzle designs. This procedure considers a pore pressure imposed on the contact interfaces. For more detailed and accurate analysis, the fundamental phenomena of erosion and pyrolysis of the insulating materials should be included. Besides, experiments should be carried out to investigate in detail the material properties related to the simulation.
For nozzle design with gaps between different materials, another important factor that should be included is the matching of their thermal expansion coefficients. Inappropriate thermal expansion coefficient matching may result in excessive stress and destroy the material. Moreover, not every best gap matching mode obtained by simulation is available considering the manufacture options. As a result, some splicing materials can be used at the very part that directly contacts with the hot gas. In view of this, the simulation should take the property of the glue into account to simulate the debonding process.
Author Contributions: Lin Sun, Futing Bao and Weihua Hui discussed and aroused the idea; Lin Sun and Ning Zhang performed the simulation and drafted the manuscript; Lin Sun finalized the manuscript. Shaozeng Wang, Nan Zhang and Heng Deng did the experiment and supplied the experiment results.

Conflicts of Interest:
The authors declare no conflict of interest.

L
Length of the nozzle (mm) d t Diameter of throat (mm) x Current location along the axis of the nozzle (mm) ρ Density (kg/m 3 ) T Temperature (K) E Young's modulus (GPa) v Poisson's ratio G Shear modulus (GPa) α Coefficient of thermal expansion (K´1) K Thermal conductivity (W/(m¨K)) C Specific heat (J/(kg¨K)) h Forced convection coefficient (W/(m 2¨K )) c p Constant-pressure specific heat of the gas (J/(kg¨K)) µ The coefficient of viscosity of the gas (m 2 /s) Pr Prandtl number m Mass flow rate (kg/s) A t Area of the throat (m 2 ) r c Curvature radius (mm) A Area of the local cross profile (m 2 )

Ma
The local Mach number R Contact thermal resistance (m¨K/(W)) σ Stefan-Boltzmann constant (W/m 2¨K4 ) T 0 Temperature of the chamber P 0 Pressure of the chamber k Specific heat ratio of the gas d Current location along path 1Ñ2Ñ3Ñ4Ñ5 (mm) D The total length of the path 1Ñ2Ñ3Ñ4Ñ5 (mm)