Study on Numerical Algorithm of the N-S Equation for Multi-Body Flows around Irregular Disintegrations in Near Space

Abstract

There has been a concern that the accurate numerical simulation of multi-body flow, which is caused by the multiple disintegrations of expired spacecraft re-entering into the near space, has a critical bottleneck impact on the falling area of the disintegrated debris. To solve this problem, an O-type grid topology method has been designed for the multi-body flow field of irregular debris formed by multiple disintegrations in near space, and a finite-volume implicit numerical scheme has been constructed for the Navier-Stokes equations to solve the aerodynamic interference characteristics of irregular multi-body flow, and further the N-S equation numerical algorithm has been established for the irregular multi-body flows in near space. The reliability of the method has been verified by the comparison of the present computation and the experiment of the low-density wind tunnel for the two-body flow of sphere, cylinder and square scripts. The objects of this study are from the multiple disintegrations of the Tiangong-1 spacecraft during uncontrolled re-entry into the atmosphere, including propelling cylinders and low-temperature lock cabinets. A series of simulations of multi-body flow mechanisms around different combinations have been carried out with varied shapes and spacing. As a result, it is found that when the distance of irregular debris (e.g., two propelling cylinders) in the near space is in the range of Δy < 3D or Δx < D, there is an obvious multi-body interference between debris, and the flow characteristics are obviously changed. When the distance between the debris in near space reaches a certain level, the influence of mutual interference can be ignored. For example, when the y-direction distance between multiple bodies is greater than 3D, the flow interference tends to be small and can be ignored, and we can regard them as two separate pieces to be carried out by the numerical prediction of flight track and falling area in engineering application. The results provide a practical design criterion for the integrated simulation platform which is used to simulate the multi-physics complex aerodynamics of space vehicles from the free-molecule flow of the outer space to the near-ground continuum flow.

1. Introduction

With the rapid development of aerospace science and technology, the research needs of related disciplines are becoming more and more urgent. The reentry process of spacecraft is an important issue in aerospace science and technology [1,2,3,4,5,6,7,8,9,10,11,12,13,14]. The full implementation of China’s manned space project and lunar exploration project has highlighted the importance of reentry technology in the space field. The flight process of spacecraft reentry into the dense atmosphere of the Earth can be divided into three categories in terms of function and destination: recovery, attack and fall [10,15,16,17,18]. The EDL (Entry, Descent and Landing) technology, which is currently attracting a lot of attention, is the technology related to the re-entry recycling process. The process includes that the spacecraft directly enters the celestial body to be landed along its orbit or changes its orbit away from its original orbit and enters the celestial body along the transformed orbit. If there is an atmosphere, it needs to pass through safely and to finally land on the surface of a celestial body smoothly by using atmospheric deceleration. The flight process of the terminal section (or including part of the middle section) of strategic and tactical ballistic missiles or some aerospace attack kinetic energy bombs is a typical case of reentry attacks. The spacecraft fallout reentry problem is neither an EDL technique nor an attack-oriented reentry/entry process, but a non-functional reentry situation such as the fallout of booster or launcher debris, failed satellites or other artificial objects in orbit after mission completion, which can be called as “non-conventional reentry problem”. The fall of a spacecraft is not a functional requirement and process relative to recovery or attack, it has thus received limited attention. The fallout area distribution and ground-based risk assessment of surviving debris from spacecraft fall are some of the areas of concern, which are the main driving forces for research on spacecraft fall analysis and forecasting. The reentry process of a large spacecraft with a truss structure at the end of its service life causes multiple disintegrations due to aerodynamic thermal effects. The fine simulation of the multi-body problems formed by the debris are the key foundations for predicting the dispersion range of the debris falling area, carrying out the numerical prediction of the airspace covered by the reentry disintegration, avoiding the risk of collision with other spacecraft, and assessing the hazard of ecosystem on the ground [15,16,17,18,19,20]. To this end, the National Key Basic Research Development Program (973 Program) project “Research on Key Basic Problems of Cross-Basin Aerodynamics and Flight Control of Space Vehicles” (Project No. 2014CB744100) combined the frontier basic research with the urgent national needs [10,17,21,22,23,24,25,26,27,28,29,30,31]. The project aims to establish a numerical prediction method for the unusual reentry trajectory of spacecraft by combining hypersonic aerodynamics, aerothermodynamics, disintegration and flight motion covering various flow regimes with an integrated platform of multi-physics field complex aerodynamics simulation of spacecraft from outer space free molecular flow to near ground continuous flow. By establishing a joint computational analysis mechanism for aerodynamic, aerothermal, structural disintegration, and flight ballistics, the research have formed a validation platform for the numerical prediction and simulation system of spacecraft de-orbiting and reentry disintegration falling zone at the end of service life [21,22,23,24,25,26,27,28,29,30,31].

The platform addresses the problem of numerical aerodynamic and thermal simulation of the flight and reentry process of large spacecraft debris by combining the Boltzmann equation collision integral physical analysis with computable modeling. The velocity distribution function equations of Boltzmann model that describes complex flow transport phenomena were established for various flow regimes during the reentry disintegration of expired spacecraft, the discrete velocity coordinate method and the improved Gaussian multiple integration method were developed, and the gas-kinetic unified algorithm (GKUA) for solving the Boltzmann model equations was established. The platform has been successfully used [17,22,25,28,32] for the study of hypersonic aerodynamic/thermal problems of the flow around multi debris generated during the reentry of large-scale spacecraft.

It was found that the disintegration of an expired spacecraft during reentry successively experiences from the highly rarefied free-molecule flow to the transition flow, slip flow and continuous flow regimes. It is a process of continuous decay of energy accumulation due to structural thermal response including deformation, melting and ablation by entering into the strong aerothermodynamics environment at the first cosmic velocity [21,25,32]. For example, the Tiangong-1 spacecraft failed after two and a half years of overdue service, and the altitude of its orbit decreased due to atmospheric drag. When it entered the atmospheric environment, it would face the reentry dis-integration and crash problems, which requires flight path prediction and falling area safety assessment and hazard analysis [16,17,18,19,20], involving many difficulties such as multiscale flow, various flow regimes, high temperature, aerodynamic heating in hypersonic flow [21,22,23,24,25]. The numerical forecasting indicates that the Tiangong-1 experienced the first disintegration at 110–105 km and multiple disintegrations at 95–56 km during reentry, and then the Mach number of the debris reduces to five, corresponding to a continuous flow regime. In this condition, the computational resources required to solve the Boltzmann model equation based on the six-dimensional phase space in physical and velocity using discrete velocity ordinate method with cost too much [32]. On the other hand, for numerical simulation of flow problems in continuous flow regime, the method based on Navier-Stokes equation is mainly used to solve the macroscopic flow problems. Among those methods, the RANS (Reynolds-Averaging Equations) based on the time-averaged form of the N-S equation, has been widely applied in engineering in continuous flow regime [33,34,35,36,37,38,39]. However, the studies on the mechanism of flow around multi-irregular debris generated by multiple disintegrations during the reentry of spacecraft in near space at the end of service life have been seldom reported.

During the process of reentry of spacecraft into near space and multiple disintegrations at the end of service, bow shock waves are generated near the leading edge of aircraft in the flow field of hypersonic dense gas with high temperature, resulting in the significant interference of strong shock waves on the surface of debris. In such a flow with high temperature, the fierce heat transfer process, including radiation heat transfer, has a dominant effect on the flow field. The extra-high temperature not only affects the aerodynamic force or momentum and aerothermodynamics characteristic of the flow, but also influences the static stability of the vehicle, thereby resulting in the multiple disintegrations of the vehicle structure through the deformation caused by the thermal stress [15,16,17,18,19,20,22,23,24,25,26,27,32,38]. To analyze the flow around debris in near space, a multi-block patched grid method is developed in this working to describe the computational region of flow around multi-irregular debris. Simultaneously, a finite volume implicit numerical scheme based on the Navier-Stokes equations for solving the problem is constructed. In addition, a numerical algorithm for the scheme is established to simulate the mechanism of flow around multibody composed of cylindrical and square cylinders under different incoming flow conditions with varying Reynolds numbers. The process of multiple disintegration of spacecraft during reentry into near space at the end of service life is shown in Figure 1, which is simplified to computational models with irregular shapes based on the features. After the extraction of flow characteristics, it is simplified to a computable model of multi-body flow around irregular shapes and different intervals, which is used as the research object of this working.

2. Physical Models and Numerical Methods

2.1. Governing Equations

During the process of off-orbit reentry and multiple disintegrations of large-scale spacecraft into near space at the end of service, due to the strong aerodynamic heat, the rapid thermal accumulation leads to the deformation and destruction of the dynamic response to the structure, and then the energy and flight speed decrease continuously with the role of aerodynamics. Specially, there is no dissociation nor ionization phenomena for the disintegrated flows around the irregular shapes. After several breakups into near space, the Mach number has become very less than 5, and the non-equilibrium effects of chemical reaction and dissociation ionization have become to be nonexistent. In terms of some practice of numerical prediction from the multiple disintegrated objects of the China’s Tiangong-1 target aircraft during uncontrolled fall reentry, it’s feasible to reveal the flow interference of multi-body flow for the propulsion gas cylinders and low-temperature lock cabinets. As the first step of the research, the Navier-Stokes (N-S) equations can be used to describe the multi-body flow around the irregular disintegrated bodies in near space. So, all of the computational simulations are carried out by an internal CFD code (Fortran language) developed by the authors. This code also develops the MPI (Message-Passing-Interface) parallel technology to improve calculation efficiency and uses the local time marching method to accelerate the convergence of the computation to a stable state.

Based on the three laws of mass conservation, momentum conservation and energy conservation followed by fluid motion, the Navier-Stokes (N-S) equations describing the multi-body flow can be deduced, and the following three assumptions are mainly introduced with the generalized Newton’s viscous stress formula, continuum hypothesis, and complete gas equation of state. If the external heat source and the volume force are not considered, the conservative form of the N-S equations in the Cartesian coordinate system can be written as [33,34,35,36,37,38]:

$$\frac{\partial \overrightarrow{Q}}{\partial t}+\frac{\partial \overrightarrow{F}}{\partial x}+\frac{\partial \overrightarrow{G}}{\partial y}+\frac{\partial \overrightarrow{H}}{\partial z}=\frac{\partial {\overrightarrow{F}}_v}{\partial x}+\frac{\partial {\overrightarrow{G}}_v}{\partial y}+\frac{\partial {\overrightarrow{H}}_v}{\partial z}$$

(1)

where, each variable is dimensionless according to the characteristic length, the incoming flow density and velocity of the research object. $\overrightarrow{Q}$ represents the conserved variable, $\overrightarrow{F}$, $\overrightarrow{G}$ and $\overrightarrow{H}$ represent the inviscid vector flux in the three coordinate directions of x, y and z, and ${\overrightarrow{F}}_v$, ${\overrightarrow{G}}_v$ and ${\overrightarrow{H}}_v$, respectively, represent the viscous vector flux in the three coordinate directions. The specific form is as follows:

$$\begin{gathered} \mathit{Q}=\left(\begin{array}{c}\rho \\ \rho u\\ \rho v\\ \rho w\\ \rho e\end{array}\right), \mathit{F}=\left(\begin{array}{c}\rho u\\ \rho u^2+p\\ \rho uv\\ \rho uw\\ (\rho e+p)u\end{array}\right), \mathit{G}=\left(\begin{array}{c}\rho v\\ \rho uv\\ \rho v^2+p\\ \rho vw\\ (\rho e+p)v\end{array}\right), \mathit{H}=\left(\begin{array}{c}\rho w\\ \rho uw\\ \rho vw\\ \rho w^2+p\\ (\rho e+p)w\end{array}\right) \\ {\mathit{F}}_v=\left(\begin{array}{c}0\\ {\tau }_{xx}\\ {\tau }_{xy}\\ {\tau }_{xz}\\ {\phi }_x\end{array}\right), {\mathit{G}}_v=\left(\begin{array}{c}0\\ {\tau }_{yx}\\ {\tau }_{yy}\\ {\tau }_{yz}\\ {\phi }_y\end{array}\right), {\mathit{H}}_v=\left(\begin{array}{c}0\\ {\tau }_{zx}\\ {\tau }_{zy}\\ {\tau }_{zz}\\ {\phi }_z\end{array}\right) \end{gathered}$$

(2)

where,

$$\begin{array}{l}{\phi }_x=u{\tau }_{xx}+v{\tau }_{xy}+w{\tau }_{xz}-q_x\\ {\phi }_y=u{\tau }_{yx}+v{\tau }_{yy}+w{\tau }_{yz}-q_y\\ {\phi }_z=u{\tau }_{zx}+v{\tau }_{zy}+w{\tau }_{zz}-q_z\end{array}$$

(3)

The components of the viscous stress term are:

$$\begin{array}{l}{\tau }_{xx}=\frac{2}{3}\mu \left(2\frac{\partial u}{\partial x}-\frac{\partial v}{\partial y}-\frac{\partial w}{\partial z}\right);\begin{array}{ccc}& & \end{array}{\tau }_{xy}={\tau }_{yx}=\mu \left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\\ {\tau }_{yy}=\frac{2}{3}\mu \left(2\frac{\partial v}{\partial y}-\frac{\partial u}{\partial x}-\frac{\partial w}{\partial z}\right);\begin{array}{ccc}& & \end{array}{\tau }_{yz}={\tau }_{zy}=\mu \left(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}\right)\\ {\tau }_{zz}=\frac{2}{3}\mu \left(2\frac{\partial w}{\partial z}-\frac{\partial u}{\partial x}-\frac{\partial v}{\partial y}\right);\begin{array}{ccc}& & \end{array}{\tau }_{zx}={\tau }_{xz}=\mu \left(\frac{\partial w}{\partial x}+\frac{\partial u}{\partial z}\right)\end{array}$$

(4)

where, $q_x$, $q_y$ and $q_z$ represent the heat flux in three directions, respectively, and the specific forms are as follows:

$$\begin{array}{l}{q}_x=-\frac{\mu }{\left(\gamma -1\right)\mathrm{Pr}}\frac{\partial T}{\partial x}\\ q_y=-\frac{\mu }{\left(\gamma -1\right)\mathrm{Pr}}\frac{\partial T}{\partial y}\\ q_z=-\frac{\mu }{\left(\gamma -1\right)\mathrm{Pr}}\frac{\partial T}{\partial z}\end{array}$$

(5)

where, ρ represents the fluid density, u, v and w represent the three-dimensional velocity components in Cartesian coordinate system, p represents the pressure, and e represents the total energy of gas per unit mass.

$$e=\frac{1}{\gamma -1}\frac{p}{\rho }+\frac{1}{2}\left(u^2+v^2+w^2\right)$$

(6)

where, μ represents the dynamic viscosity coefficient, which can be calculated by the Sutherland formula:

$$\frac{\mu }{{\mu }_0}={\left(\frac{T}{T_0}\right)}^{1.5}\left(\frac{T_0+T_s}{T+T_s}\right)$$

(7)

where, $T_0$ = 273.16 K. If the fluid is air, take ${\mu }_0$ = 1.7161 × 10−5 Pa s, $T_s$ = 110.56 K.

2.2. Finite Volume Method

The basic idea of finite volume method is to generate grids of computed flow field and integrate the differential equation to be solved into each control volume, so as to obtain a set of discrete equations. Compared with the finite difference method, the finite volume method has less dependence on the grid and is easier to deal with the flow with complex geometric shape. Therefore, the finite volume method has been widely used in computational fluid dynamics. Assuming that the grid does not change with time, the spatial derivative of differential N-S equation can be transformed into integral form by using Gauss divergence theory. In the curvilinear coordinate system, the conservative N-S equation has the following integral form:

$${\displaystyle \underset{\Omega }{∰}\frac{\partial Q}{\partial t}}dV+{\displaystyle \underset{S}{∯}(\mathit{F}-{\mathit{F}}_v)·d\mathbf{S}+}{\displaystyle \underset{S}{∯}(\mathit{G}-{\mathit{G}}_v)·d\mathbf{S}+}{\displaystyle \underset{S}{∯}(\mathit{H}-{\mathit{H}}_v)·d\mathbf{S}=0}$$

(8)

where, Ω represents the control body, S represents the control surface, and dS is the external normal area vector of each control surface of the control body.

In the calculation space under the general curvilinear coordinate system, the schematic diagram of finite volume grid element Ω (i, j, k) is in Figure 2 as follow:

Using the spatial discretization method, the control equations are discretized, and the semi discrete equations with Ω (i, j, k) lattice center physical quantities as unknowns can be obtained:

$$Vo{l}_{i,j,k}\frac{d{\mathit{Q}}_{i,j,k}}{dt}+\delta {(F+G+H-F_v-G_v-H_v)}_{i,j,k}=0$$

(9)

where, the latter item represents the net flux of fluid outflow Ω (i, j, k), in the specific form:

$$\begin{array}{ll}\delta {(\mathit{F}+\mathit{G}+\mathit{H}-{\mathit{F}}_v-{\mathit{G}}_v-{\mathit{H}}_v)}_{i,j,k}& ={\left((\mathit{F}-{\mathit{F}}_v)\cdot \mathrm{S}\right)}_{i+1/2,j,k}-{\left((\mathit{F}-{\mathit{F}}_v)\cdot \mathrm{S}\right)}_{i-1/2,j,k}\\ & +{\left((\mathit{G}-{\mathit{G}}_v)\cdot \mathrm{S}\right)}_{i,j+1/2,k}-{\left((\mathit{G}-{\mathit{G}}_v)\cdot \mathrm{S}\right)}_{i,j-1/2,k}\hfill \\ & +{\left((\mathit{H}-{\mathit{H}}_v)\cdot \mathrm{S}\right)}_{i,j,k+1/2}-{\left((\mathit{H}-{\mathit{H}}_v)\cdot \mathrm{S}\right)}_{i,j,k-1/2}\hfill \end{array}$$

(10)

In the calculation space, the interface area and volume of the grid element are 1, so the integral N-S equation based on the finite volume method of the element center can be written as follow:

$$\begin{array}{l}{\left(\frac{\partial \mathit{Q}}{\partial t}\right)}_{i,j,k}+{\left(\mathit{F}-{\mathit{F}}_v\right)}_{i+1/2,j,k}-{\left(\mathit{F}-{\mathit{F}}_v\right)}_{i-1/2,j,k}+{\left(\mathit{G}-{\mathit{G}}_v\right)}_{i,j+1/2,k}\\ \hfill -{\left(\mathit{G}-{\mathit{G}}_v\right)}_{i,j-1/2,k}+{\left(\mathit{H}-{\mathit{H}}_v\right)}_{i,j,k+1/2}-{\left(\mathit{H}-{\mathit{H}}_v\right)}_{i,j,k-1/2}=0\end{array}$$

(11)

2.3. Spatial Discrete Scheme

With high discontinuous resolutions, the original Roe’s FDS is one of the most famous CFD schemes. In this paper, the inviscid flux is divided by the Roe scheme. The Roe scheme belongs to the Godunov type method, which transforms the nonlinear Riemann problem into a linear problem. By solving the approximate solution of the linearized Riemann problem at each grid interface, the solution of the whole flow field is obtained. The key is to construct a linearized approximate matrix to approximately replace the Jacobi matrix of inviscid flux. Roe scheme has natural high resolution for contact discontinuity. Since the characteristics of viscous action are similar to contact discontinuity, roe scheme also has high viscous resolution, which is not available in general schemes.

The Jacobi matrix of inviscid flux can be diagonalized. Taking $\xi$ direction as an example, the Jacobi matrix of flux can be decomposed into its corresponding eigenvalues and eigenvectors:

$$A=R\Lambda R^{-1}$$

(12)

The numerical flux of Riemann problem can be written as:

$$\begin{array}{ll}{F}_{i+1/2,j,k}& =\frac{1}{2}\left[R(\Lambda +\left|\Lambda \right|)R^{-1}{Q}_{i,j,k}+R(\Lambda -\left|\Lambda \right|)R^{-1}{Q}_{i+1,j,k}\right]\\ & =\frac{1}{2}(R\Lambda R^{-1}{Q}_{i,j,k}+R\Lambda R^{-1}{Q}_{i+1,j,k})-\frac{1}{2}\left[R\left|\Lambda \right|R^{-1}(Q_{i+1,j,k}-Q_{i,j,k})\right]\\ & =\frac{1}{2}(F_{i,j,k}+F_{i+1,j,k})-\frac{1}{2}\left|A\right|(Q_{i+1,j,k}-Q_{i,j,k})\end{array}$$

(13)

Therefore, in the Roe format, the interface flux in the $\xi$ direction is [40]:

$$F_{i+1/2,j,k}=\frac{1}{2}[\tilde{F}(Q_R)+\tilde{F}(Q_L)-\left|\tilde{A}(Q_L,Q_R)\right|(Q_R-Q_L)]$$

(14)

where, the matrix $\tilde{A}$ is the Roe average Jacobian matrix, the $Q_L$ and $Q_R$ are the state variables on the left and right sides of the interface. The superscript “–” represents the Roe average, which is defined as follows:

$$\begin{array}{ll}\tilde{\rho }=\sqrt{{\rho }_L{\rho }_R}& \tilde{u}=(u_L+u_R\sqrt{{\rho }_R/{\rho }_L})/(1+\sqrt{{\rho }_R/{\rho }_L})\\ \tilde{v}=(v_L+v_R\sqrt{{\rho }_R/{\rho }_L})/(1+\sqrt{{\rho }_R/{\rho }_L})\hfill & \tilde{w}=(w_L+w_R\sqrt{{\rho }_R/{\rho }_L})/(1+\sqrt{{\rho }_R/{\rho }_L})\\ \tilde{H}=(H_L+H_R\sqrt{{\rho }_R/{\rho }_L})/(1+\sqrt{{\rho }_R/{\rho }_L})\hfill & {\tilde{c}}^2=(\gamma -1)[\tilde{H}-({\tilde{u}}^2+{\tilde{v}}^2+{\tilde{w}}^2)/2]\end{array}$$

(15)

After sorting, the dissipation term in the Roe format is given as follow:

$$\left|\tilde{A}\right|(Q_R-Q_L)=\left[\begin{array}{c}{\alpha }_4\\ \tilde{u}{\alpha }_4+{\widehat{\xi }}_x{\alpha }_5+{\alpha }_6\\ \begin{array}{l}\tilde{v}{\alpha }_4+{\widehat{\xi }}_y{\alpha }_5+{\alpha }_7\\ \tilde{w}{\alpha }_4+{\widehat{\xi }}_z{\alpha }_5+{\alpha }_8\end{array}\\ \tilde{H}{\alpha }_4+\tilde{\overline{u}}{\alpha }_5+\tilde{u}{\alpha }_6+\tilde{v}{\alpha }_7+\tilde{w}{\alpha }_8-\frac{{\tilde{c}}^2{\alpha }_1}{\gamma -1}\end{array}\right]$$

(16)

where,

$$\begin{array}{lll}{\alpha }_1=\left|S\right|·\left|\tilde{\overline{u}}\right|(∆\rho -\frac{∆p}{{\tilde{c}}^2}) \hfill & {\alpha }_4={\alpha }_1+{\alpha }_2+{\alpha }_3 & {\alpha }_7=\left|S\right|·\left|\tilde{\overline{u}}\right|(\tilde{\rho }∆v-{\widehat{\xi }}_y\tilde{\rho }\overline{u})\\ {\alpha }_2=\frac{1}{2{\tilde{c}}^2}\left|S\right|·|\tilde{\overline{u}}+\tilde{c}|(∆p+\tilde{\rho }\tilde{c}∆\tilde{u})\hfill & {\alpha }_5=\tilde{c}({\alpha }_2-{\alpha }_3)& {\alpha }_8=\left|S\right|·\left|\tilde{\overline{u}}\right|(\tilde{\rho }∆w-{\widehat{\xi }}_z\tilde{\rho }\overline{u})\\ {\alpha }_3=\frac{1}{2{\tilde{c}}^2}\left|S\right|·|\tilde{\overline{u}}+\tilde{c}|(∆p+\tilde{\rho }\tilde{c}∆\tilde{u})\hfill & {\alpha }_6=\left|S\right|·\left|\tilde{\overline{u}}\right|(\tilde{\rho }∆u-{\widehat{\xi }}_x\tilde{\rho }\overline{u})\end{array}$$

(17)

$$\begin{array}{lll}{\widehat{\xi }}_x=\frac{{\xi }_x}{\sqrt{{\xi }_x^2+{\xi }_y^2+{\xi }_z^2}}& {\widehat{\xi }}_y=\frac{{\xi }_y}{\sqrt{{\xi }_x^2+{\xi }_y^2+{\xi }_z^2}}& {\widehat{\xi }}_z=\frac{{\xi }_z}{\sqrt{{\xi }_x^2+{\xi }_y^2+{\xi }_z^2}}\\ \tilde{\overline{u}}=\frac{(\tilde{u}{\xi }_x+\tilde{v}{\xi }_y+\tilde{w}{\xi }_z)}{\sqrt{{\xi }_x^2+{\xi }_y^2+{\xi }_z^2}}& \left|S\right|=\frac{\sqrt{{\xi }_x^2+{\xi }_y^2+{\xi }_z^2}}{J}& \Delta (•)={(•)}_R-{(•)}_L\\ \tilde{\rho }=\sqrt{{\rho }_L{\rho }_R}& \tilde{D}=\sqrt{\frac{{\rho }_R}{{\rho }_L}}& (\widetilde{•})=\frac{{(•)}_L+{(•)}_R\tilde{D}}{1+\tilde{D}}\end{array}$$

(18)

The Roe scheme can finally be written as follow:

$${\mathit{F}}_{1/2}=\frac{1}{2}({\mathit{F}}_L+{\mathit{F}}_R)-\frac{1}{2}\mathit{R}\left|\widehat{\Lambda }\right|{\mathit{R}}^{-1}\Delta \mathit{Q}$$

(19)

where, ^ means the Roe average, R is the right eigenvector matrix and Λ is the diagonal matrix.

2.4. Time Marching Method

The time marching solution method used in the steady calculation in this working is the LU-SGS (Lower-Upper Symmetric Gauss Seidel) method proposed by Yoon and Jameson [41]. This method has become as the main implicit calculation method of CFD, and the stability of LU-SGS has been highly evaluated. Compared with the explicit method, the implicit time discretization method has obvious advantages in numerical stability and convergence efficiency.

The semi-discrete ordinary differential equation can be obtained by using the first-order difference for the time term in the governing equation, the implicit treatment for the inviscid flux and explicit treatment for the viscous flux, as follow:

$$\frac{\Delta {\mathit{Q}}_{i,j,k}}{\Delta t}+\delta {\left(\mathit{F}+\mathit{G}+\mathit{H}\right)}_{i,j,k}^{n+1}=\delta {\left({\mathit{F}}_v+{\mathit{G}}_v+{\mathit{H}}_v\right)}_{i,j,k}^n$$

(20)

where, ΔQ = Qn + 1 − Qn, n + 1 represents the next moment, n represents the current moment.

By linearizing the inviscid flux and omitting the second-order and higher-order terms, the following equation can be obtained:

$$\begin{array}{l}\frac{\Delta {\mathit{Q}}_{i,j,k}}{\Delta t}+{\left(A\Delta \mathit{Q}\right)}_{i+\frac{1}{2},j,k}-{\left(A\Delta \mathit{Q}\right)}_{i-\frac{1}{2},j,k}+{\left(B\Delta \mathit{Q}\right)}_{i,j+\frac{1}{2},k}-{\left(B\Delta \mathit{Q}\right)}_{i,j-\frac{1}{2},k}+{\left(C\Delta \mathit{Q}\right)}_{i,j+\frac{1}{2},k}\\ -{\left(C\Delta \mathit{Q}\right)}_{i,j-\frac{1}{2},k}=-\delta {\left(\mathit{F}+\mathit{G}+\mathit{H}\right)}_{i,j,k}^n+\delta {\left({\mathit{F}}_v+{\mathit{G}}_v+{\mathit{H}}_v\right)}_{i,j,k}^n={\mathit{R}}_{i,j,k}^n\end{array}$$

(21)

where, $R_{i,j,k}^n$ represents the residual at time n.

$$A=\frac{\partial \mathit{F}}{\partial \mathit{Q}}, B=\frac{\partial \mathit{G}}{\partial \mathit{Q}}, C=\frac{\partial \mathit{H}}{\partial \mathit{Q}}$$

(22)

By splitting the inviscid flux Jacobian matrix on the element interface according to the positive and negative eigenvalues, the following equations can be obtained,

$$\begin{array}{c}{\left(A\Delta \mathit{Q}\right)}_{i+\frac{1}{2},j,k}=A_{i,j,k}^{+}\Delta {\mathit{Q}}_{i,j,k}+A_{i+1,j,k}^{-}\Delta {\mathit{Q}}_{i+1,j,k}\\ {\left(A\Delta \mathit{Q}\right)}_{i-\frac{1}{2},j,k}=A_{i-1,j,k}^{+}\Delta {\mathit{Q}}_{i-1,j,k}+A_{i,j,k}^{-}\Delta {\mathit{Q}}_{i,j,k}\\ {\left(B\Delta \mathit{Q}\right)}_{i,j+\frac{1}{2},k}=B_{i,j,k}^{+}\Delta {\mathit{Q}}_{i,j,k}+B_{i,j+1,k}^{-}\Delta {\mathit{Q}}_{i,j+1,k}\\ {\left(B\Delta \mathit{Q}\right)}_{i,j-\frac{1}{2},k}=B_{i,j-1,k}^{+}\Delta {\mathit{Q}}_{i,j-1,k}+B_{i,j,k}^{-}\Delta {\mathit{Q}}_{i,j,k}\\ {\left(C\Delta \mathit{Q}\right)}_{i,j,k+\frac{1}{2}}=C_{i,j,k}^{+}\Delta {\mathit{Q}}_{i,j,k}+C_{i,j,k+1}^{-}\Delta {\mathit{Q}}_{i,j,k+1}\\ {\left(C\Delta \mathit{Q}\right)}_{i,j,k-\frac{1}{2}}=C_{i,j,k-1}^{+}\Delta {\mathit{Q}}_{i,j,k-1}+C_{i,j,k}^{-}\Delta {\mathit{Q}}_{i,j,k}\end{array}$$

(23)

Substituting into Equation (21), the following equation can be obtained:

$$\begin{array}{l}\left[I+\Delta t\left(A_{i,j,k}^{+}-A_{i,j,k}^{-}+B_{i,j,k}^{+}-B_{i,j,k}^{-}+C_{i,j,k}^{+}-C_{i,j,k}^{-}\right)\right]\Delta {\mathit{Q}}_{i,j,k}\\ +\Delta t\left(-A_{i-1,j,k}^{+}\Delta {\mathit{Q}}_{i-1,j,k}-B_{i,j-1,k}^{+}\Delta {\mathit{Q}}_{i,j-1,k}-C_{i,j,k-1}^{+}\Delta {\mathit{Q}}_{i,j,k-1}\right)\\ +\Delta t\left(A_{i+1,j,k}^{-}\Delta {\mathit{Q}}_{i+1,j,k}+B_{i,j+1,k}^{-}\Delta {\mathit{Q}}_{i,j+1,k}+C_{i,j,k+1}^{-}\Delta {\mathit{Q}}_{i,j,k+1}\right)=-\Delta t{\mathit{R}}_{i,j,k}^n\end{array}$$

(24)

Further, the LU decomposition of inviscid flux Jacobian matrix has the following form:

$$(L+D)D^{-1}(D+U)\Delta \mathit{Q}=-\Delta t{\mathit{R}}_{i,j,k}^n$$

(25)

where,

$$\begin{array}{l}L=-\Delta t(A_{i-1,j,k}^{+}+B_{i,j-1,k}^{+}+C_{i,j,k-1}^{+})\\ D=I+\Delta t\chi {({\sigma }_A+{\sigma }_B+{\sigma }_C)}_{i,j,k}I\\ U=\Delta t(A_{i+1,j,k}^{-}+B_{i,j+1,k}^{-}+C_{i,j,k+1}^{-})\end{array}$$

(26)

and,

$$A^{\pm }=\frac{1}{2}\left(A\pm \chi {\sigma }_{A}I\right);B^{\pm }=\frac{1}{2}\left(B\pm \chi {\sigma }_{B}I\right);C^{\pm }=\frac{1}{2}\left(C\pm \chi {\sigma }_{C}I\right)$$

(27)

The Jacobian matrix of viscous flux is complex and difficult to diagonalize. The influence of viscosity is usually not considered in implicit processing. However, in the region with significant viscous effect and dense grid, this simplification may cause instability of numerical calculation. To solve this problem, the following form of stability factor is introduced to improve the computational convergence:

$$v_{\xi }=\frac{2\mu \cdot M_{\infty }}{\rho \cdot Re}{\left|\nabla \xi \right|}^2, v_{\eta }=\frac{2\mu \cdot M_{\infty }}{\rho \cdot Re}{\left|\nabla \eta \right|}^2, v_{\zeta }=\frac{2\mu \cdot M_{\infty }}{\rho \cdot Re}{\left|\nabla \zeta \right|}^2$$

(28)

Thus, the modified form of diagonal operator D can be obtained:

$$D=I+\Delta t\chi {({\sigma }_A+{\sigma }_B+{\sigma }_C)}_{i,j,k}I+2\Delta t{\left(v_{\xi }+v_{\eta }+v_{\zeta }\right)}_{i,j,k}I$$

(29)

The Gauss Seidel iterative method is used to solve the time advance problem in two steps [41]. The L-Block operator and u-block operator are executed successively, and the whole field is scanned once before and after.

$$\Delta {\mathit{Q}}_{i,j,k}^{*}=\frac{\left\{-\Delta t{\mathit{R}}_{i,j,k}^n+\Delta t[{(A^{+}\Delta {\mathit{Q}}^{*})}_{i-1,j,k}+{(B^{+}\Delta {\mathit{Q}}^{*})}_{i,j-1,k}+{(C^{+}\Delta {\mathit{Q}}^{*})}_{i,j,k-1}]\right\}}{1+\Delta t{[\chi ({\sigma }_A+{\sigma }_B+{\sigma }_C)+2(v_{\xi }+v_{\eta }+v_{\zeta })]}_{i,j,k}}$$

(30)

$$\Delta {\mathit{Q}}_{i,j,k}=\Delta {\mathit{Q}}_{i,j,k}^{*}-\frac{\Delta t[{(A^{-}\Delta \mathit{Q})}_{i+1,j,k}+{(B^{-}\Delta \mathit{Q})}_{i,j+1,k}+{(C^{-}\Delta \mathit{Q})}_{i,j,k+1}]}{1+\Delta t{[\chi ({\sigma }_A+{\sigma }_B+{\sigma }_C)+2(v_{\xi }+v_{\eta }+v_{\zeta })]}_{i,j,k}}$$

(31)

Thus, the value of the conserved variable at n + 1 time is updated, as follow:

$${\mathit{Q}}_{i,j,k}^{n+1}={\mathit{Q}}_{i,j,k}^n+\Delta {\mathit{Q}}_{i,j,k}$$

(32)

According to the previous research experience, if the explicit scheme is adopted, the calculation time advance step is limited by the stability conditions of the scheme. When the shape of the aircraft is very complex, the non-uniform change of the measurement coefficient of the physical space grid will lead to a small calculation time step, especially for the calculation of the hypersonic flow around the irregular and complex aircraft in the continuous and near continuous flow region, Therefore, the research and development of implicit methods has become a more concerned problem. How to combine the physical processes of different scales for computable modeling and realize the application of implicit efficient algorithm in the near-space hypersonic flow of three-dimensional complex irregular aircraft has become an important development direction of implicit scheme. On the other hand, the use of a single grid system to represent complex shapes, whether the computational grid design is reasonable and whether the quality is high, will also directly affect the efficiency of computational convergence [26,27,42,43], and the physical space multi block docking grid technology can well solve this problem.

3. Model and Grid

3.1. Computable Model and Computational Conditions

The computable models in this work are selected as the side-by-side propulsion gas cylinders and low-temperature lock cabinets from the multiple disintegrated objects of the China’s Tiangong-1 target aircraft during uncontrolled fall reentry. As the real shapes of the disintegrated objects and debris are exceedingly complex and puzzling to be simulated by CFD. In order to reveal the flow interference of multi-body flow and find the engineering method to resolve the dispersion law of disintegration in near space, the original layouts of propulsion gas cylinders and low-temperature lock cabinets can be modeled as the computable two-body models shown in Figure 3b,d, respectively. The side-by-side propulsion gas cylinders are rotating bodies, which are composed of hemispheres and cylinders at both ends. The diameter of hemispheres at both ends is D = 450 mm and the length of cylinder section is L = 577 mm. The low-temperature lock cabinet is hexahedral, and the overall dimension is 450 mm × 450 mm × 495 mm. The computational conditions are taken as the uncontrolled reentry flight environment of the Tiangong-1 target aircraft as the research object at the altitude of 40247 m. Due to the continuous attenuation of its energy, the flight speed and Mach number are continuously reduced, which is 3.0, so that the thermochemical non-equilibrium effects weaken rapidly and there is no dissociation nor ionization phenomenon of air in this working.

3.2. Grid Generation Method for Multi-Body Flow Field around Irregular Disintegration

The reasonable design and high-quality generation of computational grid is the prerequisite of CFD calculation. The grid generation technology of complex shape has become one of the key factors affecting CFD calculation. The quality of grid directly affects the computational accuracy of numerical solutions, and this influence is decisive in many cases. There are many methods to generate the flow field grid around complex aircraft, such as structural grid, unstructured grid, Cartesian grid and so on. The structural grid is widely used because of its simple data structure, convenient program design, high computational efficiency and accuracy. The solution domain in general physical space is a region with complex geometry and known closed boundary, and the corresponding value of computational coordinates on the closed boundary is given. This is a typical boundary value problem of partial differential equation. The boundary value problem of elliptic partial differential equation is well posed. Therefore, by solving the elliptic Poisson equation and using the non-homogeneous term to adjust the grid density distribution in the solution domain, the corresponding flow field grid layout can be obtained, and the body fitted coordinate technology can be used to express the shape of the aircraft surface to generate the body fitted coordinate structured grid. However, it is difficult to generate a high-quality single block structure grid to discretize the physical space flow field due to the complex shape in engineering application, especially the irregular disintegration flying object formed during the falling and disintegration of the spacecraft at the end of service, which needs to be realized by multi block grid. Even for the multi-body flow around a single object with a simple shape, the multi-block grid technique is also needed to generate high-quality structural grids. Therefore, given the difficulty of generating body fitted structure grid in a single connected domain due to the multi-body flow around irregular disintegration, this paper adopts butt joint multi-block structure grid, automatically processes any butt joint relationship according to the transfer information between blocks, transfers flow field parameters between blocks through the butt joint surface, and treats the adjacent elements of the butt joint surface as internal points without interpolation, so as to automatically ensure the flux conservation of the whole flow field, It can effectively ensure the conservation of numerical flux in the flow field and is easy to adapt to the flow around aircraft with various complex shapes.

The grid structure is convenient to compute and organize data. However, constructing topology is the most important and difficult part of generating structural grid. After many seemingly simple shapes are combined, the topology becomes extremely complex. For example, it is relatively easy to generate the whole flow field grid of a single disintegrated object. However, in the whole flow field grid of multiple disintegrated objects, it is necessary to consider how to connect the grid of each disintegrated object. In the process of mesh generation, the O-type topology is adopted near the wall. On the one hand, it ensures the density and orthogonality of the mesh near the wall. On the other hand, it is conducive to control the mesh near the wall of each component in the O-type topology, which is no longer transmitted to the far field, so as to control the total amount of the mesh. Grid densification shall be carried out at key flow parts such as the boundary layer, while ensuring that the height of the normal grid of the first layer of the wall is y + < 1. The wall condition adopts the assumption of non-slip adiabatic wall. The overall grid meets the three elements of smooth transition, flow field adaptability and orthogonality, which is conducive to improving the calculation accuracy. The grid layout of symmetry plane and wall is shown in Figure 4. The division of the flow field around the computational model is shown in Figure 1 according to the butt-joint and multi-block structure grid, and the O-type topology near the wall is adopted to ensure the grid density and orthogonality near the wall, so as to obtain the computational grid for computing the multi-body flow field around the irregular disintegration in this working, as shown in Figure 5.

3.3. Grid Sensitivity Test

The grid is an important component in CFD (Computational Fluid Dynamics). It is necessary to perform the gird sensitivity test. Taking the paralleled and placed propulsion cylinder bottle and low-temperature lock cabinet with Δy = 4.5D as an example, a number of test cases with different gird resolutions and grid Reynolds numbers are tested as shown in

Table 1. Grid resolutions of different cases for the side-by-side propulsion cylinder (Ma = 3.0, H = 40 km).

	Case A (Coarse)	Case B (Medium)	Case C (Fine)
Grid Reynolds Number	200	100	50
Total Grid Number	560,651	989,901	1,845,371

, where the grid Reynolds number is defined as $R{e}_{cell}=\rho V∆l/\mu$, where $∆l$ is the normal length of the grid cell nearest the wall. For the same computational condition, the grid Reynold number can reflect the density of the grid near the wall. The schematic diagram of grid on wall and symmetry plane is shown in Figure 4 and Figure 5.

Figure 6 shows the dimensionless pressure contours of the paralleled and placed propulsion bottle and low-temperature lock cabinet with different grid resolutions. There is no much difference of macroscopic flow field information for different grid resolutions.

Table 2. Axial and normal force coefficients when the propulsion cylinders and low-temperature lock cabinets are placed and paralleled with different grids.

		Case A (Coarse)	Case B (Medium)	Case C (Fine)
Propulsion bottle	Axial force Coefficients	0.175	0.189	0.189
	Normal force Coefficients	0.186 × 10−5	0.197 × 10−5	0.197 × 10−5
Lock cabinet	Axial force Coefficients	0.079	0.088	0.088
	Normal force Coefficients	0.777 × 10−6	0.828 × 10−6	0.839 × 10−6

shows the axial and normal force coefficients with different grid resolutions. The force coefficients of Case A are relatively less than that of Case B and Case C. The force coefficients of Case B and Case C are almost the same. It is indicated that the grids of Case B and Case C satisfy the requirement of grid convergence, so the grid system corresponding to Case B is employed. To sum up the above, when the grid Reynolds number is less than 100, the computational grid can satisfy the requirement of grid convergence. Therefore, the grid Reynolds number of all the grid in this working is used to be less than 100, and the used total grid number is requested to be more than 989,901 to reliably solve the macroscopic flow field and aerodynamic coefficients around the propulsion cylinder bottles and low-temperature lock cabinets.

4. Numerical Simulation and Analysis on Flows around Irregular Multi-Debris in Near Space

4.1. Validation/Verification of Numerical Algorithm

In order to investigate the characteristics of flows around debris generated by multiple disintegration of the Tiangong-type vehicle during off-orbit reentry, the accuracy of the numerical method is firstly validated. The interference of aerodynamic properties of simplified debris models obtained by the low-density wind tunnel experiment [44] is analyzed and compared with the previous study. Based on the features, the debris generated by multiple disintegration of the Tiangong-type vehicle during reentry are classified and then simplified to three typical shapes: sphere, spherocylinder and cuboid. Subsequently, the experiment on the interference of aerodynamic properties of the debris is carried out. The figures of the experiment models are illustrated in Figure 7, with the diameter of sphere being 100 mm, the diameter of spherocylinder being 100 mm and the length of cuboid being 200 mm. the whole experiment is conducted at the Mach number of 9.96, the total temperature of 1040 K and the total pressure of 2.03 Mpa.

The pressure contours of interference between spheroid and sphere debris obtained through the numerical computation are shown in Figure 8, the changes of flow structures in the pressure contours at different relative distance d can be observed in the figure, where d indicates the upper and lower distances, and the interference model is located above the test model with the aligned head. When the spherocylinder and sphere are placed and paralleled, two smooth detached shock waves are generated in front of the spherical and spherocylindrical debris, forming a region with high pressure and temperature. The detached shock waves are intersected between the spherocylinder and sphere, generating a Mach disk. When the distance d of sphere and spheroid is less than 2.5 times the diameter D of the sphere, the detached shock waves are reflected by the Mach disk and then impinge on the surface of debris, seriously impacting on the flow field, as well as the stress and temperature distribution in the debris. With the increase of the distance between the sphere and spherocylinder, the location of the intersection of the shock waves moves downstream. In this regard, after reflection, the shock waves directly propagate downstream, and the interaction is weakened. When the distance between two debris is above 2.5 times the diameter of the sphere, both the Mach disk and the reflected shock waves have little influence on the debris, thus can be regarded as two independent objects and the numerical simulation on the flow around the debris can be carried out, respectively.

Figure 9 and Figure 10 display the changes of aerodynamic forces at different distance between two bodies when the interference model is placed above the measuring model. In the figures, H indicates that the interference model is aligned with the force measuring model, and V implies that they are placed up and down. The comparison of aerodynamic coefficients measured in the experiment and the computation through numerical simulation indicates the good agreement between numerical results and experiments, which further validates the accuracy and reliability of the numerical method for flows around multi-body with irregular shapes in near space established in this working. When the distance between spherical and sphero-cylindrical debris is small, the debris interact seriously with each other, the detached shock waves intersect and generate a Mach disk that reflect the shock waves and directly impinge on the surface of the upwards and downwards debris, which significantly affects the aerodynamic coefficients such as the axial force and normal force that changes with the distance between two debris. It can be indicated from Figure 9 and Figure 10 that the computational axial and normal force coefficients are in good agreement with that of the experiments from the low-density wind tunnel [44] for both spherocylinder—sphere and spherocylinder—spherecylinder flow interference, which verifies the validation and accuracy of the numerical algorithm for multi-body flows around irregular disintegrations in near space.

4.2. Numerical Simulation and Analysis on Flows around Irregular Multi-Debris in Near Space

4.2.1. Numerical Simulation and Analysis on Flow around Side-by-Side Propelling Cylinders in Near Space

In this section, the numerical simulation on the flows around two side-by-side propulsion cylinder bottles is carried out using the present numerical algorithm that focuses on the flow problems around the irregular multi-body based on N-S equation in near space. The Mach number contours are shown in Figure 11 with the main inflow Mach number of 3, and the flying altitude is 40 km with the distance Δy = 1.5D, 2D and 3.5D between the axes of two propulsion cylinders. As displayed in the figure, the detached shock waves are generated near the leading edge of each propulsion cylinder, which intersect with each other in the flow region between two propulsion bottles, resulting in complex interference between Mach disk and the flow such as the intersection, reflection, and fusion of shock waves.

The distance between two propulsion cylinders has significant impact on the characteristics of the flow field. When Δy = 1.5D, the detached shock waves generated near the leading edge of the propulsion bottles intersect with each other at X = −26 mm (the stagnation point of the bottle is located at X = 0 mm), generating a clear Mach disk. After intersection, transmitted shock waves impinge on the walls of the propulsion cylinder, inducing shock wave-boundary layer interaction, leading to a sudden increase in the pressure near the debris wall, thus changing the aerodynamic coefficients of the debris. When Δy = 2D, the detached shock waves are generated near the leading edge of each propulsion bottle and intersected with each other at X = 143 mm, the interference is still remarkable, generating obvious reflected shock waves that are similar as those in the hypersonic inlet. As Δy continue to increase, when it reaches 3.5D, the intersection of two shock waves moves further downstream. They intersect at X = 640 mm and the transmitted shock waves directly propagate to the downstream so that the interference of the shock waves on the debris is wakened. As is displayed in

Table 3. Minimum distance between shock wave and stagnation point of debris and location of intersection of shock waves.

Δy	1.5D	2D	3.5D
Minimum distance	109 mm	109 mm	109 mm
position of shocks intersection	X = −26 mm	X = 143 mm	X = 640 mm

, the minimum distance between the shock waves and the stagnation point of the propulsion bottle remains at 109 mm, indicating that the position of the axes of the detached shock waves in the front of the side-by-side propulsion bottles will not change with the variation of the distance between two propulsion bottles, in this regard, the interference of the flow field has no influence on the detached shock waves. When Δy is sufficiently small, almost equal to D, the two propulsion bottles are combined into a single object. In this way, the flow becomes similar to that of a single body. Conversely, when Δy is sufficiently large, the flow structures of two propulsion bottles will not affect each other, as a result, the two shock waves develop independently. Once the aerodynamic interference between the two propulsion cylinder bottles disappears, the two separated objects with mutual interference can be regarded as two separate ballistics computations for the numerical prediction [17,23,24,25,32] of the flight trajectory in the aerodynamics integration ballistics of the two propulsion cylinders.

In order to quantitatively analyze the flow around side-by-side propulsion bottles placed with different spacing, the axial and normal force coefficients on one of the bottles when Δy = 1.5D, 2D and 3.5D is listed in

Table 4. Axial and normal force coefficients of side-by-side propulsion bottles when Ma = 3 and the altitude is 40 km.

Δy	1.5D	2D	3.5D
Axial force coefficient	0.0880	0.0873	0.0755
Normal force coefficient	0.0367	0.0322	1.055 × 10−4

. As is shown in the table, with the increase of Δy, the magnitude of axial force coefficient remains almost the same, because it mainly influenced by the detached shock waves near the leading edge of the propulsion bottle. When y = 3.5D, the axial force coefficients decrease remarkably, this is because that when the position of the intersection of shock waves moves downstream, the transmitted shock wave directly enters the downstream and no longer affects the flow field near the wall of the propulsion bottles. When Δy = 1.5D or 2D, the shock wave interacts strongly between two propulsion bottles, resulting in a significant multi-body interference, thus improving the normal force coefficients. Since the normal force coefficients are mainly due to the interference of the side-by-side propulsion bottles, which is not obvious when Δy = 3.5D, the order of magnitude of the normal force coefficients in that condition is only 10⁻⁴. Moreover, the shock wave interaction results in a significant change in the characteristics of the flow field. When Δy is small, the mutual interference of two debris is strong, and the changes of the flow around the debris can’t be explained simply by analyzing the flow around a single debris. When Δy = 3.5D, the force coefficients of a single propulsion bottle, especially the normal force coefficients, are significantly reduced. The detached shock waves intersect with each other and generate two transmitted shock waves that directly propagate downstream. In this condition, the interference between two debris is substantially weakened so that its effect on the flow field can be neglected, thereby the forecast of the flight path of the debris can be conducted through the independent numerical simulation on each debris.

4.2.2. Numerical Simulation and Analysis on the Flow around Staggered Propulsion Cylinders in near Space

In this section, the analysis on the flow around two staggered propulsion cylinder bottles is conducted, the horizontal distance Δy keeps constant at 2D and the axial distance of the bottles is set as 0, D, and 3D, respectively. The dimensionless pressure contours and streamline structures are displayed in Figure 12a–c, the Mach number contours are illustrated in Figure 12d–f, and the dimensionless temperature fields are shown in Figure 12g–i. As is shown in the figures, a detached shock wave is generated in front of each propulsion bottle. For Δx = 0, two detached shock waves and the Mach disk of the wave intersection are clearly reproduced with severe flow interference, however, the symmetry distribution of up and down macroscopic flow variables and streamlines structure is a distinct feature shown in Figure 12a,d,g, respectively. For Δx = D, the above detached shock wave seriously impacts on the detached shock wave below, and the asymmetric flow structure presents in both streamlines and pressure, density and temperature distribution shown in Figure 12b,e,h, respectively. For Δx = 3D, the flow interference hardly exists and it is almost impossible for the above detached shock wave to impacts on the detached shock wave below, however, and the flow structure behind the upper propulsion cylinder is affected with small disturbance by the detached shock wave below shown in Figure 12c,f,i, respectively, only because the flow parameters of the detached shock wave below are still affected with some weak perturbation from the above detached shock wave. It has been indicated from Figure 12 in the comparison with Figure 11 that the disturbed region of flow interference between two staggered propulsion cylinders is much larger than that of two side-by-side propulsion bottles.

Specially, the cases when Δx = 0 have been analyzed in Section 4.2.1, after the intersection of two bow shock waves, the transmitted shock wave interacts with the boundary layer on the debris wall, the interference between two debris is significant. In the other hand, with the increase of Δx, the flow influenced by the propulsion bottle in the back moves downward, gradually entering the region influenced by the propulsion bottle in the front. As a result, the interference to the flow field near the propulsion bottle in the back becomes weaker. The trailing edge of the shock wave generated by the front propulsion bottle intersects with that generated by the rear propulsion bottle, affecting the flow field structure near the rear propulsion bottle. With the increase of Δx, the intensity of shock wave generated by the front propulsion bottle at the intersection becomes weaker, so its interference on the flow field of the rear propulsion bottle is suppressed. When Δx is sufficiently large, the whole rear propulsion bottle is submerged in the flow field influenced by the front propulsion bottle. The axial and normal force coefficients of the rear and front propulsion bottles are shown in

Table 5. Axial and normal force coefficient when propelling bottles are placed at different position.

Δx		0	D	3D
Propulsion bottle (front)	Axial force coefficient	0.0873	0.0872	0.0881
	Normal force coefficient	0.0322	1.147 × 10−4	5.3974 × 10−5
Propulsion bottle (back)	Axial force coefficient	0.0873	0.0920	0.1001
	Normal force coefficient	0.0322	0.0157	4.8940 × 10−3

. It can be seen from the table that when Δx increases, there is only a slight growth in the axial force coefficients of the front propulsion cylinder. However, as Δx increase from D to 3D, the normal force coefficients decrease significantly to approximately 5 × 10⁻⁵, indicating that the interference on the front propulsion bottle is extremely weak. The axial force coefficients of the rear propulsion bottle increase with the growth of Δx. this is because that the detached shock wave generated by the front propulsion bottle impinges on the bow shock wave near the leading edge of the rear propulsion bottle, affecting the structure of detached shock wave generated by the rear propulsion bottle, thereby increasing the axial force coefficients. When Δx = 3D, the rear propulsion bottle gradually enters the influenced flow field near the front propulsion bottle, thus the interference on the rear propulsion bottle is weakened. In this regard, the order of magnitude of the normal force coefficients of the rear gas cylinder is reduced to 10⁻³. It has been indicated from Table 5 in the comparison with Table 4 that when Δx is more than 3D, the force coefficients of a single propulsion bottle, especially the normal force coefficients, are significantly reduced to be minimal with 5 × 10⁻³ or 5 × 10⁻⁵, and the interference between two debris can be substantially neglected.

4.2.3. Numerical Simulation and Analysis on Flows around Paralleled Placed Propulsion Cylinder Bottles and Cryostats

In reality, the debris are often composed by irregular blocks with different shapes. In this section, taking the debris of propulsion cylinder bottles and cryostats as an example, the numerical simulation on the flows around debris in practical problems with the axial distance Δy = 1.5D, 3D, and 4.5D is carried out. As shown in the dimensionless pressure contours in Figure 13a–c, smooth detached shock waves are formed near the leading edge of both the propulsion cylinder bottle and the cryostat, generating the regions of high pressure. When Δy = 1.5D, the two detached shock waves intersect with each other and form the strong shock wave with a saddle shape between two pieces of debris. Similar to the flows around two paralleled propulsion cylinder bottles when Δy = 1.5D, the shock waves have a significant impact on the flow field. Nevertheless, with the increase of Δy, the location of intersection of the shock waves moves downstream and the interference between two debris is weakened.

In order to quantitatively analyze the influence of the relative position of the propulsion cylinder bottle and the cryostat on the flow, the axial and normal force coefficients when Δy equals to 1.5D, 3D and 4.5D are computed and listed in

Table 6. Axial and normal force coefficients when propulsion cylinder bottle and cryostat are placed and paralleled with different Δy.

Δy		1.5D	3D	4.5D
Propulsion bottle	Axial force coefficients	0.098	0.187	0.189
	Normal force coefficients	0.197 × 10−1	0.662 × 10−5	0.197 × 10−5
Cryostat	Axial force coefficients	0.188	0.087	0.088
	Normal force coefficients	0.217 × 10−1	0.688 × 10−4	0.839 × 10−6

. It can be indicated from Table 6 and Figure 13 that the effect of the detached shock wave is huge on the flow interference with the different spacing of Δy = 1.5D, 3D and 4.5D. For Δy = 1.5D, two detached shock waves form the saddle strong shock wave, which results in the great difference of flow variables from the front and rear of the propulsion cylinder bottle and the cryostat with the rear vacuum region in the rear of propulsion cylinder bottle and the cryostat, and then the different axial and normal force coefficients. From Δy = 3D to Δy = 4.5D, two detached shock waves form the Mach disk between the propulsion cylinder bottle and the cryostat and behind the cryostat, which results in the increased axial force coefficients from Ca = 0.098 to 0.187 and 0.189 for the propulsion cylinder bottle with longer length of 577 mm than that of Ca = 0.188 to 0.087 and 0.088 for the cryostat with the short length of 495 mm. As the axial force originates from the pressure difference between the front and back of the debris, the completely different pressure distribution of front and back flow determines the different axial force around the propulsion cylinder bottle and the cryostat, respectively. However, normal force originates from disturbing force perpendicular to surface and flow direction, and decreases greatly with increasing spacing from Δy = 1.5D to Δy = 3D and 4.5D, specially, the normal force coefficient has become as Cn = 0.197 × 10⁻⁵ and 0.839 × 10⁻⁶ for the propulsion cylinder bottle and the cryostat, respectively, when the spacing Δy = 4.5D, and this can be judged that the flow interference between two debris no longer exists.

5. Conclusion and Expectation

This study focuses on multi-body flow formed by multiple disintegrations of an expired spacecraft during off-orbit reentry. Firstly, an O-grid topology is designed using the multi-block patched grid technique to describe the flow field computation region. Then a finite-volume implicit scheme for solving the Navier-Stokes equation is constructed to study the aerodynamic interference characteristics of a multi-body flow formed by irregular debris. A numerical algorithm is then established to solve the Navier-Stokes equation with the butt-joint and multi-block structure grid technique for multi-body flows formed by irregular debris in near space. Finally, a simulation study on the flow mechanism of multi-body supersonic flow caused by different combinations of the Tiangong-1 target vehicle’s debris in near space is carried out. The results agree well with results of wind tunnel experiments, proving the reliability of the numerical method. The results show that there is strong interaction between irregular debris, e.g., propulsion cylinders, when their distance is Δy < 3D or Δx < D, and aerodynamic characteristics have significant changes. However, the interaction is negligible if the distance increases to a certain magnitude. When the distance between the debris in near space reaches a certain level, the influence of mutual interference can be ignored, and debris can be regarded as two separate pieces of disintegrated wreckage to be carried out the engineering application design of numerical prediction software for aerodynamic integrated with ballistic flight track and falling region.

This study fills the blank in the field of multi-body flow mechanism formed by irregular debris during reentry of expired spacecraft, and combined with the integrated simulation platform for solving complex aerodynamic problems with multiscale multi-physics fields from free molecular flow in outer space to near-ground continuous flow, makes up for deficiencies of macroscopic hydrodynamic methods in computing the aerodynamic interactions characteristics of a multi-body flow formed by irregular debris. The results provide support for applications of inertial-energy-inspired thermodynamic non-equilibrium effects in aerodynamic characteristics of multi-body flow mechanism formed by irregular debris during reentry of expired spacecraft. It lays a solid foundation and application support for accurate and reliable simulation of aerodynamic characteristics of multi-body flow around different flight heights, different disintegration flight speeds, and different flight angles of attack, which is furtherly considered of the thermodynamic non-equilibrium effect induced by internal energy.