Abstract

Multi-interface is a typical feature of solid rocket motor nozzles, and the interface evolution law is an important aspect to analyze the thermo-structural response of solid rocket motor nozzles. In this paper, based on the friction coefficient test at different temperatures, a coupled thermo-structural analysis model of solid rocket motor nozzle considering the variation of friction coefficient under operating conditions was established. Firstly, the friction coefficient was tested to model the variation at different temperatures. Then, by adopting structure gap, variable friction coefficient, thermal contact resistance, and friction heat production, a strongly coupled non-linear model was established. Simulations using the non-linear model and traditional model were performed, in which the stress of the throat insert was within the required stress range of the material. The ground firing test result demonstrated the validity of the analysis model, and the non-linear model was in a better agreement with the firing test than the traditional model. Therefore, it can be concluded that with a more specific friction coefficient to represent the friction behaviour of the different parts of the nozzle, the strongly coupled non-linear model established in this paper can reflect the essence of thermo-structural response of solid rocket motor nozzle.

1. Introduction

Nowadays, with the development of computer technology, numerical simulation of the solid rocket motor (SRM) working process has become more and more refined. In particular, SRM nozzles are subjected to extremely high temperature and pressure. Using high-energy propellant, the temperature can reach 3500 K or higher, and for some special SRM, the pressure may reach about 10 MPa or higher [1]. Thus, structural integrity is of great importance, and many researchers have focused on different relevant aspects, such as material fabrication [2], material intrinsic performance characterization [3], anti-ablation performance [4], and thermo-structural analysis methods [5].

For general numerical simulation methods, in the early days, thermo-structural analysis consists mainly of thermal and mechanical loads. Kumar et al. [6] took thermal and mechanical loads into account to simulate the thermo-structural response of a nozzle, which showed that mechanical loads were less important than thermal load. Mehta and Narayana [7] adopted temperature-dependent thermo-physical properties for the SRM to investigate the axisymmetric, transient, anisotropic heat conduction problems and simulate the thermal stress with two degrees of freedom. In terms of the liquid rocket motor, Hu [810] simulated the thermo-structural response under special conditions using a large expansion ratio nozzle that leads to the flow separation phenomenon. The simulation results were in great agreement with the test results. Gong [11] studied the difference in the fluences of thermal loading and mechanical loading, and the results showed that the thermal loading had an important influence on the stress of the throat insert, which was proved by a ground test. Later, the surface recession of the nozzle gradually became a research hotspot. Kuo and Keswani [12] established a full model containing both diffusion and chemical kinetics effects for evaluating the recession process of a carbon/carbon nozzle. Bianchi [13] took homogeneous phase reactions and material nonlinearity into account and analyzed the coupled case of ablation and heat flow. In view of radiation heat transfer, Zhang [14] studied the ablation of the C/C material, and predicted the safety performance of the rocket nozzle in the ignition stage. Qin [15] used mesoscopic/microscopic composite modelling technology to study the ablation process and ablation mechanism of carbon/carbon composites at a mesoscopic/microscopic scale. Liu [16] studied 4D braided carbon/carbon composite material ablation behaviour and found that the mechanical erosion effect was determined by a particle velocity threshold. Previous studies have shown that the thermo-structural analysis for SRM nozzle is always a significant project, mainly including the coupled analysis of heat transfer and structural response, contact thermal resistance, interface debonding, and constitutive model analysis. In addition, new concepts and novel analysis methods are highlighted. On that basis, it can be seen the interfaces between the different parts of the SRM nozzle gradually become a hot spot for predicting thermo-structural response.

Multi-interface is a typical feature of modern SRM nozzle and it is important to analyze the interface evolution during the SRM working process. As far as we know, few research have considered the interfacial contact status in the thermo-structural field, especially the exact friction coefficient that may affect thermo-structural response of the contact interface between several parts of the nozzle. The ground test results have proved that both thermal and mechanical load can lead to the failure of the bonded interface of the nozzle parts and cause the interface frictional. This phenomenon needs to be described by an accurate friction coefficient. It is undeniable that many researchers have attached great importance to the debonding process [1726]. However, previous studies mainly adopted the empirical value. For example, Wang et al. [27] carried out a multifield (flow-thermal-mechanical) coupled numerical simulation on the mesh-based parallel code coupled interface, and the results indicated that both the temperature field and stress field obtained by the coupled algorithm were slightly lower than the values obtained by the engineering algorithm. Although the algorithm was verified by the test results, the friction coefficient which used a constant value of 0.25 was not validated. Liu 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; however, no further study was performed. Subsequently, Nigar et al. [28] cited Liu’s test results to investigate the cracking behaviour by combining computational fluid dynamics simulation of the supersonic flow, finite element modeling of the thermal shock, and extended finite element method-based analysis. The results showed that most cracks initiate at the flow surface and propagate parallel to the surface in the form of Mode II cracks. Notably, Liu’s work mainly tested the carbon/carbon material, which is different from the materials used in Nigar’s paper. Thus, Nigar’s work is still not specific. Therefore, although researchers have paid great attention to interface frictional behaviour, there is still a lack of targeted studies.

To address this issue, the main objective of this paper is to establish a thermo-structural analysis model considering the variation of friction coefficient under operating conditions to reflect the essence of thermo-structural response of the SRM nozzle. First, a special test was used to investigate the friction coefficient between carbon/carbon material and silica/phenolic resin at different temperatures. The temperatures were set close to the values during firing test. Further, a thermo-structural analysis model considering the temperature-dependent friction coefficient was established. Nonlinear thermo-structural analysis for an SRM nozzle was performed using this model. Also, the thermo-structural analysis was performed when the friction coefficient was ignored. Finally, the ground tests were conducted, and the strain of the metal case showed that an accurate friction coefficient led to a more accurate simulation result. The results have demonstrated that the thermo-structural analysis model established in this paper is reasonable.

2. Friction Coefficient Test and Results

2.1. Principle and Process

During firing test, the temperature of the interface between several parts of the nozzle will increase, and the interface condition will change from a bonded state to a frictional state. Previous studies have shown that the change occurs very fast and can be ignored. Therefore, in this paper, from the beginning to the end of the firing test, the interface was considered to be frictional. To describe the frictional behaviour more accurately, a specific test of the friction coefficient at different temperatures should be performed.

Carbon/carbon material and silica/phenolic resin were chosen for testing on a high-temperature friction and wear testing machine, model HT-1000, which is shown in Figure 1. It is mainly composed of a computer control system, a high-temperature furnace, a mechanical friction system, and a loading mechanism. Its working principle is as follows: (1)Heat the samples in the high-temperature furnace to a set value. In this study, the temperature was set to 25°C (298.15 K), 50°C (323.15 K), 100°C (373.15 K), 200°C (473.15 K), 300°C (573.15 K), 400°C (673.15 K) and 500°C (773.15 K);(2)Load the demand load by the loading mechanism and drive the carbon/carbon samples on the sample plate rotating to rub with the silica/phenolic resin samples(3)Measure the actual temperature in the high-temperature furnace and the friction force by the computer(4)Calculate the friction coefficient by dividing the friction force by the loading force

The test process is as follows: (1)Cut the samples into blocks with a diameter of 25 mm and a height of 6 mm(2)Heat the samples at 10°C/min and when the temperature reaches the demand point, maintain the temperature for 10 min(3)Start the test and keep it running for 30 min. The load is 1000 g, the wear diameter is 8 mm, and the speed is 388 r/min(4)Collect the test data and repeat the test process at each temperature(5)Analyze the friction coefficient curve over time and get the friction coefficient for future use

2.2. Results and Analysis

Figure 2 shows the test results. As the temperature increases, the friction coefficient decreases with different slopes. In the beginning, the decline is slow. But when the temperature reaches about 500 K, the decreasing slope becomes larger. Then the slope slows down. The friction coefficient changes from about 0.156 to about 0.085, and the decline proportion is about 45.51%. After the test, the test samples were checked. Due to the high temperature, the surface of samples was carbonized and covered by lots of wear debris. This could be one of the reasons for the decrease in the friction coefficient because the wear debris can act as a lubricant and the friction characteristics of the contact surface are reduced.

It can be concluded that constant values from experience or values from reference may not reflect the actual contact behaviour, especially for those exposed to high temperatures. So, a variable friction coefficient between carbon/carbon material and silica/phenolic resin needs to be used, which is more targeted.

3. Establishment of the Thermo-Structural Analysis Model

Firstly, a non-linear strong coupled thermo-structural analysis model considering structure gap, variable friction coefficient, thermal contact resistance, and friction heat production of the interface was established, as shown in Figure 3. All the factors were implemented using ANSYS Parameter Design Language (APDL) in ANSYS Workbench platform. APDL codes were developed to realize the strongly coupled simulation function. The 2-D 8-Node Coupled-Field Solid element plane223 was chosen by ET command. For non-linear contact analysis, KEYOPT and RMODIF commands were executed to simulate contact functions, such as contact heat transfer. MP command was used to simulate the emission of the materials.

For element plane223, it usually involves just one analysis that uses a coupled-field element type containing all necessary degrees of freedom. Coupling was handled by calculating element matrices or element load vectors that contained all necessary terms. Its governing equations can be expressed as follows [29]: where is element mass matrix; is acceleration vector; is second time derivative of temperature; is element structural damping matrix; is element thermoelastic damping matrix; is element-specific heat matrix; is velocity vector; is time derivative of temperature; is element stiffness matrix; is element thermoelastic stiffness matrix; is element thermal conductivity matrix; is nodal displacement vector; is temperature vector; is the sum of the element nodal force vector; and is the sum of the element heat generation load and element convection surface heat flow vectors.

The detailed meaning of each symbol can also be found in ANSYS Mechanical APDL Theory [29].

The modeling method for the analysis can be described as follows: (1)Flow field parameters, such as temperature, pressure, and Mach, were calculated by a quasi-one-dimensional isentropic model(2)Detailed heat transfer coefficient and recovery temperature can be obtained using Bartz formula(3)The thermo-structural analysis was conducted by adopting pressure load and heat transfer load [21];(4)During the simulation, detailed nonlinear contact boundary, friction behaviour and thermal contact resistance of the interfaces between carbon/carbon nozzle insert and silica/phenolic resin part(s) were all taken into account, as shown in Figure 3. The variable friction coefficient obtained in the previous test is a highlight of the model

4. Computation Model

4.1. Working Condition

The nozzle used in the firing test is 2D axisymmetric, and a simple 2D sketch can fully describe the design details of the nozzle. Figure 4 shows the sketch of the nozzle. It consists of an insert made of carbon/carbon, a back surface made of silica-phenolic-resin, a convergent section made of carbon-phenolic-resin and silica-phenolic-resin, a divergent section made of carbon-phenolic-resin and silica-phenolic-resin, and a metal case made of titanium alloy. As a general SRM nozzle, this nozzle is also axially symmetric. Therefore, to simplify the analysis process such as modeling, mesh quantity, etc., a 2D axisymmetric model was adopted.

The nozzle was subjected to a total pressure of 12 MPa and a total temperature of 3534 K, and the operating time was 25 s. Besides, the environment pressure is 0.1 MPa and the environment temperature is 295.15 K.

4.2. Assumptions

There are still some factors that may affect the thermo-structural response of the nozzle but are not included in this paper. Therefore, the following assumptions were made [21, 30]: (1)Ablation of the materials is neglected for a fast solution(2)Structure gap appears as soon as the motor starts to work(3)The friction coefficient is only affected by temperature because as the temperature increases, the surface morphology will change due to material pyrolysis(4)Thermal contact resistance is a constant value, 10−3 m2·K/W [23];(5)Energy generated by the frictional movement of the surfaces is equally distributed to the contact and target surface

4.3. Material Model

The temperature of SRM gas could reach 3000 K or even higher. At this temperature level, material physical property parameters will vary, and to represent this, temperature-dependent material models were established for each material. The detailed properties are shown in Tables 13 [21].

ρ: Density (kg/m3)

T: Temperature (K/°C)

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))

4.4. Boundary Conditions
4.4.1. Principle, Formulas, and Results

(1) Quasi-One-Dimensional Isentropic Flow. Themo-structural analysis includes two types of boundary conditions, thermal and structural boundaries. These boundaries mainly arise from the internal flow. Previous research has proved that the quasi-one-dimensional isentropic flow analysis can be used to predict the boundary conditions in the thermo-structural simulation for SRM nozzle [20, 21]. Therefore, the both boundaries are calculated by quasi-one-dimensional isentropic flow field, as shown in Figure 3. The relationship of temperature, pressure and Mach number at any section along the nozzle axis can be expressed as [20, 21]: 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.

(2) Heat Transfer Boundary. By adopting the third type boundary condition in heat transfer theory, heat transfer coefficient and temperature near the nozzle wall can be provided. The heat transfer coefficient in the SRM nozzle can be calculated by the Bartz formula (1): where h is the forced convection coefficient; C =0.026 for subsonic flow and C =0.023 for supersonic flow; dt is the diameter of the throat; cp is the constant-volume specific heat of the gas; μ is the coefficient of viscosity of the gas; Pr is the Prandtl number; is the mass flow rate; At is the area of the throat; rc is the curvature radius; A is the area of the local cross profile; σ is a modified parameter considering the vibration of the boundary layer; Tw is the temperature of the wall; T0 is the chamber temperature; k is the ratio of specific heat; and Ma is the local Mach number.

The temperature near the nozzle wall is defined as recovery temperature and it can be calculated by the following equation: where T and Ma are the gas static temperature and Mach number of the local cross profile, respectively.

The Prandtl number Pr and the viscosity of the gas μ can be calculated by the following equations: where is the average molecular weight of the gas and k =1.2 for hot gas in SRM.

(3) Results. Temperature, heat transfer coefficient, and pressure along the nozzle axis were calculated by Eqs. (2)–(8), and the results are shown in Figure 5. From the figure, the heat transfer coefficient reaches the maximum value at the throat, indicating a tough heat transfer process. Due to the turbulence layer effect, the recovery temperature at the entrance of the nozzle throat is a little higher than the total temperature of the gas. From the nozzle inlet to outlet, the pressure first decreases quickly, and then the decreasing trend becomes smooth. The patterns of the temperature, heat transfer coefficient, and pressure curves are all consistent with the theory of SRM nozzle flow.

4.4.2. Thermal Boundaries

Thermal boundaries for this nozzle mainly include heat transfer between the outermost titanium alloy and the air, heat transfer between the inner wall of the nozzle, and the hot gas. The detailed thermal boundary conditions are shown in Figure 6. The heat transfer coefficient between the outermost titanium alloy and the air was set to the empirical value of 5 W/m2·K. The heat transfer coefficient between the inner wall of the nozzle and the hot gas was set to a constant value of 5893.5 W/m2·K for the inner wall which is labeled red triangles, because gas temperature and velocity around this location remain nearly constant. In contrast, for the inner wall labeled red squares, the convection coefficient and ambient temperature are shown in Figure 5(a).

4.4.3. Structural Boundaries

Structural boundaries for this nozzle mainly contain fixed support and pressure, as shown in Figure 7. Since the nozzle is always stuck to the combustion chamber of SRM, a fixed support type is used on several surfaces labeled blue lines. For pressure loads onto the inner wall of the nozzle and the wall labeled red triangles, a constant pressure value of 12 MPa is used since gas pressure around this location remain nearly constant; and for the inner wall labeled red squares, the pressure is shown in Figure 5(b).

4.5. Initial Structural Gap

Structural gaps are always exist due to high-temperature degeneration of the fitting surface and glue. This is the main cause of the non-linear response. To control the gap size, rules are defined in Reference [17]: (1)For the non-metal and metal fitting surface, the optimum gap size may be in the range of 0.1 to 0.2 mm, and the maximum is not greater than 0.25 mm(2)For the metal and non-metal fitting surface, the designed gap size can be smaller

Based on the above two principles, the sizes of the possible passageways that flame may leak along were designed, as shown in Figure 8 and Table 4.

4.6. Mesh Generation and Mesh Independence Study

Mesh quality can affect the precision of simulation results, and different mesh quantities may also lead to different simulation results. In this paper, three types of meshes were generated to check the mesh independence. The temperature and stress curves along the interface between the nozzle insert and back surface were selected for comparison.

For mesh independence investigation, three types of meshes, i.e., coarse (45,637 nodes), medium (59,524 nodes), and fine (65,594 nodes), were adopted in this study. Figure 9 shows the temperature distributions for three different types of meshes and Figure 10 shows the Von Mises stress distributions for three different types of meshes. Table 5 shows the differences in temperature and stress with three different types of meshes.

From Figures 9 and 10, and Table 5, we can conclude that with the medium mesh density, the results of both Von Mises stress and temperature are mesh independent. Figure 11 shows the mesh generation result of the thermo-structural simulation in this paper.

5. Results, Analysis and Verification

This section comes to an overview of the temperature and stress status of this nozzle. To be specific, two thermal-structural models were used for comparison, i.e., the non-linear strongly coupled thermo-structural analysis model established in this paper and a traditional model that ignores the non-linear behaviour.

5.1. Temperature Status Analysis

The temperature distribution produced by both models were almost identical. The results obtained by the non-linear strongly coupled thermo-structural analysis model were used as an example in the analysis. The temperature distribution at different instances during the firing test, shown in Figure 12, was used to analyze the heat transfer. Due to higher thermal conductivity, the temperature of the insert rose more rapidly than in other parts. During the process, the temperature of the other parts rose slowly. Until the end of the process, at 25 s, the temperature of the metal case remained the environment temperature, which indicated that the insulator of the nozzle played an effective role in suppressing heat transfer and protecting the structural integrity of the nozzle.

5.2. Stress Status Analysis

The Von Mises stress field graphs of the whole nozzle at different instances during firing test are shown in Figure 13. Overall, there are obvious differences between the Von Mises stress distributions produced by the non-linear model and the traditional model. First, the maximum stress value from the non-linear model is much larger than that from the traditional model. Second, according to the results of the non-linear model, it can be observed that several points in the contact parts around the insert have high stress. However, from the results of the traditional model, the stress values of these points are not prominent, which indicates that the traditional model cannot reflect the non-linear movement of the interfaces around the nozzle insert. In conclusion, the non-linear model is more suitable to predict the structure integrity.

The detailed analysis of stress distributions in different directions of the nozzle insert are shown in Figure 14. Obviously, using the non-linear model, the stress values of the nozzle insert are larger, but the maximum tensile or compression stress is mainly at the cusp point of the nozzle insert, which means that the stress is concentrated. From the results by the traditional model, the stress distribution graphs are uniform, which indicates that the stress concentration is not significant. The reason for the differences in stress concentration between the traditional model and the non-linear model is as follows: (1)In the non-linear model, all the parts can move freely during the working process of SRM. Due to the difference in the thermal expansion coefficient, thermal conductivity, and specific heat of each thermal protection material, each part deforms differently. For example, the surface of the nozzle insert is heated firstly and expands toward both sides, but the other parts will not deform with the insert and stress concentration will occur(2)In the traditional model, all the parts are bonded together and deform synchronously. As a result, the stress distribution graphs are uniform. The differences in the stress distribution by the two models are particularly significant in the locations with stress concentration in the non-linear model results

The stress extremum in each direction at the end of working time was obtained by the non-linear model. From the results, the radial stress range was [-31.346 MPa, 24.344 MPa], the axial stress range was [-155.48 MPa, 1.5234 MPa], and the circumferential stress range was [-26.288 MPa, 16.508 MPa]. In addition, the maximum stress was the cusp stress concentration. The actual stress and the permissible stress range are shown in Table 6. It is clear that the throat insert structure has high reliability.

5.3. Nonlinear Contact Results Analysis

Non-linear contact analysis can give a history of the number of contact for each contact pair. With the flame leak criteria (FLC) established in [21], this paper analyzed two possible flame leak passageways. Figure 15 gives the history of the contact number of fitting surfaces 1 to 8 defined in Table 4. It shows that after SRM firing, fitting surfaces 2, 5, 6, 7, and 8 close immediately. Based on FLC, for both possible flame leak passageways, although there are no contact numbers for a short period in the beginning, the gas flow will not leak into the passageways at a small duration. In other words, there is no flame leak risk for both possible flame leak passageways. Thus, the nozzle design is suitable for further manufacture and firing test.

5.4. Comparisons and Analysis of Numerical and Experimental Results

After the ground firing test, the nozzle was dissected, as shown in.

Figure 16 From the cross-section, the throat insert is structurally intact and no defects and cracks are observed. This indicated that the nozzle was in good condition and the above analysis based on simulation results were correct. Besides, during the firing test, the strain of the metal case was monitored at point A, as shown in Figure 17.

The strain results from the non-linear model, traditional model, and firing test are shown in Table 7. It can be found that the strain data obtained by the non-linear model was closer to the measured values in the ground firing test, while the errors in the traditional model results were much larger. The non-linear model developed in this paper is suitable to analyze the thermo-structural response of the nozzle.

Figure 18 shows an overview of the metal case. From the figure, there is no flame leakage at the two locations marked by the circles. It confirms that the nonlinear contact results analysis is valid.

6. Conclusion

In this paper, in order to predict the thermo-structural response of SRM nozzle, a strongly coupled non-linear thermo-structural analysis model considering the variation of friction coefficient under operating conditions was established based on friction coefficient test at different temperatures. In addition, similar to the established models in previous studies, structure gap, thermal contact resistance, and friction heat conduction were also considered as important factors. Then, a ground firing test was conducted to measure the strain of the point at the metal case. The contact non-linear analysis results were also validated by the test. Therefore, the model was proved to be effective, and some conclusions are drawn as follows: (1)Friction coefficient test showed that as the temperature increased, the friction coefficient slowly decreased in the range of [300 K, 500 K]. When the temperature reached the range of [500 K, 600 K], the friction coefficient decreased rapidly. The reason may be that the wear debris generated by the carbonized surface act as a lubricant(2)A strongly coupled non-linear analysis methodology was established by involving structure gap, variable friction coefficient, thermal contact resistance, and friction heat production, combined with quasi-one-dimensional isentropic flow results. Among these factors, the variation of friction coefficient under operating conditions plays an important role in the modeling of the non-linear contact simulation(3)The simulation results of the non-linear model and traditional model showed that the stresses of nozzle insert in different directions were all within the required stress range of the carbon/carbon material. The integrity of the structure can be guaranteed. Both models can be used to analyze the thermal structural integrity of nozzles(4)The analysis of the nonlinear contact results showed that both possible flame leak passageways close immediately after the SRM starts working. Based on the FLC in previous studies, there is no flame leak for this nozzle. The ground firing test showed no flame heating phenomenon at the metal case, which supports the simulation results(5)During the firing test, the strain response was monitored. Compared to the firing test, the non-linear result had the error of about 3.288% and the traditional model result had the error of about 11.99%. This indicated that the non-linear simulation results were more consistent with the firing test results than the traditional model. The strongly coupled non-linear thermo-structural analysis model is proved to be effective in predicting the thermo-structural response of SRM nozzle

Nomenclature

At:The nozzle throat area
A:The local area along the nozzle axis
c:Specific heat
C:A constant: C =0.026 for subsonic flow and C =0.023 for supersonic flow
cp:Constant-volume specific heat of the gas
dt:Diameter of the throat
E:Young’s modulus
G:Shear modulus
h:Convection coefficient
K:Thermal conductivity
k:Specific heat ratio of the gas
:Mass flow rate
Ma:Local Mach number
:The average molecular weight of the gas
P0:Chamber pressure
Pr:Prandtl number
rc:Curvature radius
SRM:solid rocket motor
T:Temperature
T0:Chamber temperature
Tw:Temperature of the wall
Tr:Recovery temperature
:Element structural damping matrix
:Element thermoelastic damping matrix
:Element specific heat matrix
:Sum of the element nodal force vector
:Element stiffness matrix
:Element thermoelastic stiffness matrix
:Element thermal conductivity matrix
:Element mass matrix
:Sum of the element heat generation load and element convection surface heat flow vectors
:Temperature vector
:Time derivative of temperature
:Second time derivative of temperature
:Nodal displacement vector
:Velocity vector
:Acceleration vector
ρ:Density
v:Poisson’s ratio
α:Coefficient of thermal expansion
μ:The coefficient of viscosity of the gas
σ:A modified parameter considering the vibration of the boundary layer.

Data Availability

Since most of the data in this manuscript are commercially confidential, I cannot make them fully available. In the future, I can share some data with reviewers and readers if necessary.

Conflicts of Interest

The authors declare that there is no conflict of interest in the publication of this paper.

Acknowledgments

This work was financially supported by the Shanghai Academy of Spaceflight Technology Innovation Fund (Program no. SAST 2019-107) and the Natural Science Basic Research Program of Shaanxi (Program no. 2020JQ-172).