Nonlinear Modeling Study of Aerodynamic Characteristics of an X38-like Vehicle at Strong Viscous Interaction Regions

Strong viscous interaction and multiple flow regimes exist when vehicles fly at high altitude and high Mach number conditions. The Navier–Stokes(NS) solver is no longer applicable in the above situation. Instead, the direct simulation Monte Carlo (DSMC) method or Boltzmann model equation solvers are usually needed. However, they are computationally more expensive than the NS solver. Therefore, it is of great engineering value to establish the aerodynamic prediction model of vehicles at high altitude and high Mach number conditions. In this paper, the hypersonic aerodynamic characteristics of an X38-like vehicle in typical conditions from 70 km to 110 km are simulated using the unified gas kinetic scheme (UGKS), which is applicable for all flow regimes. The contributions of pressure and viscous stress on the force coefficients are analyzed. The viscous interaction parameters, Mach number, and angle of attack are used as independent variables, and the difference between the force coefficients calculated by UGKS and the Euler solver is used as a dependent variable to establish a nonlinear viscous interaction model between them in the range of 70–110 km. The evaluation of the model is completed using the correlation coefficient and the relative orthogonal distance. The conventional viscous interaction effect and rarefied effect are both taken into account in the model. The model can be used to quickly obtain the hypersonic aerodynamic characteristics of X38-like vehicle in a wide range, which is meaningful for engineering design.


Introduction
The viscous interaction effect, which describes the mutual interaction process between the boundary layer and the outer inviscid flow, is one of the three main effects [1] on hypersonic vehicles for ground-to-flight extrapolation. Depending on the degree of feedback from the inviscid flow on the boundary layer, strong viscous interaction and weak viscous interaction can be defined.
Traditionally, a similarity parameter, χ = M 3 ∞ √ C/ √ Re, is used to ascertain whether an interaction region is strong or weak. Re = ρ e U e x/µ e is the conventional Reynolds number based on properties, ρ e , U e and µ e at the outer edge of the boundary layer: C = µ w ρ w /(µ e ρ e ). Large values of χ correspond to the strong interaction and small values of χ indicate a weak region. For pressure and force coefficients on simple configurations such as a flat plate or a sharp cone, a different correlation parameter, ν ∞ = M ∞ √ C/ √ Re, is usually used. The study of viscous interaction correlation for force coefficients derived from the space shuttle program has identified a modified viscous interaction parameter ν ∞ [2], which has been widely used in the literature to correlate the aerodynamic characteristics obtained by different means such as wind tunnel, flight, or numerical calculation. ν ∞ is defined as ν ∞ = M ∞ √ C / √ Re L∞ . Re L∞ is the Reynolds number based on the characteristic length of the vehicle. C = µ T ∞ /(µ ∞ T ). T /T ∞ is the ratio of the reference temperature in the boundary layer to the incoming flow temperature.
with DSMC results, Titarev also evaluated the BGK and Shakhov model equations as applied to hypersonic flows for both aerodynamics and heat transfer in [14][15][16]. Apart from the unstructured mesh technique used in Nesvetay, another efficient approach based on an adaptive octree velocity mesh is proposed by Baranger [17]. The octree mesh contains many fewer points than a traditional Cartesian mesh. In 2010, Xu [22] proposed the unified gas kinetic scheme (UGKS) method which is based on the integral solution of the model equation. The NS solution can be recovered from the UGKS in the hydrodynamic limit [23]. Good agreements between UGKS and DSMC results have also been achieved in rarefied regime [22][23][24][25][26][27]. Many advanced techniques such as the adaptive velocity method [28,29], implicit method [30][31][32][33][34][35], multigrid method [36], and memory-saving method [37] have been implemented. UGKS has been widely used in the simulation of flow fields from low speed to high speed, from continuum flow to rarefied flow [22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40]. For hypersonic validations and applications, Jiang [31] has conducted a UGKS simulation and verified its accuracy by comparing the pressure, stress, and heat flux distributions on an M = 25 cylinder for different regimes with DS2V results. Li [40] has conducted a kinetic blind comparative study on the aerodynamic characteristics of a complex-scaled X38-like vehicle which is the same as the one under study in the current paper. The free-stream Mach number is 8 with four different Knudsen numbers, 0.00275, 0.0275, 0.275, and 2.75. Two in-house kinetic solvers are used based on the DSMC method and UGKS method, respectively. Despite having different methods (statistical vs. deterministic) and different meshes (unstructured vs. structured), both UGKS and DSMC solvers gave similar and reasonably consistent results. The average relative errors for the lift and drag coefficients are only 0.98% and 2.01%, respectively.
Based on the above understanding and our practical experiences with UGKS in the past decade, the goal of this paper was to establish a viscous interaction model applicable to all regimes. An in-house UGKS solver was used to predict the aerodynamic characteristics of a complex X38-like configuration at high altitude (70-110 km) and high Mach number (≥10). The viscous interaction correlation method derived from the space shuttle program [2] was used for reference.
The difference between the aerodynamic characteristics obtained by UGKS and the solution of inviscid Euler equations was used as the dependent variable. A prediction model relating the difference and the viscous interaction parameter was proposed. The model was evaluated using the concepts of correlation coefficient and relative orthogonal distance. Some new cases were selected and calculated by UGKS and the prediction model to verify the accuracy of the model prediction results.

The Inviscid Solver
The governing equations are the three-dimensional compressible Euler equations in general curvilinear coordinates. The equations are discretized based on the finite volume method and solved by the implicit LUSGS method. See [41] for more details.

The Viscous Solver
The governing equation is the Shakhov model equation [42] which can be written in non-dimensional form: Here f is the distribution function which is a function of the space x, the particle velocity u, and time t. τ is the collision time. g M is the local Maxwellian distribution function. The second term in f + is a correction term based on the original BGK model Entropy 2022, 24, 836 4 of 18 equation in order to obtain a reasonable Prandtl number, Pr.c and q are the random velocity vector and the heat vector, respectively. µ, ρ, and p are the non-dimensional viscosity, density, and pressure, respectively. Re ∞ and M ∞ are the free stream Reynolds number and Mach number, respectively. Dimensional free stream density ρ ∞ , velocity modulus |U ∞ |, temperature T ∞ , and viscosity µ ∞ are used to obtain the non-dimensional macroscopic quantities in the following way The superscript ' * ' denotes dimensional quantities. The power-law intermolecular interaction µ = T ω is assumed. The total length of the vehicle L re f is used as scale of Unless explicitly specified, all variables in the following are non-dimensional. The relations between the macroscopic conserved quantities Q, the stress P, the heat q and the distribution function are where ψ is the vector of moments and dΞ = dudvdw is the volume element in the phase space.
In UGKS, at the cell interface (i + 1/2,j,k) an integral solution of the Shakhov model in the following form is used to construct the solution: where f + = g + g + will be approximated separately. The subscripts i,j,k and l,m,n denote the indexes in three structured physical mesh directions and three Cartesian velocity mesh directions, respectively. x = x i+1/2 − u l (t − t ) is the particle trajectory and f 0 is the initial gas distribution function at the beginning of each time step around the cell interface x i+1/2 at particle velocity u = (u l , v m , w n ). As the distribution function inside each control volume is known at the beginning of each time step. f 0 can be obtained using TVD reconstruction.
The equilibrium state g around the cell interface x i+1/2 can be expanded with two slopes where H[x] is the Heaviside function. g 0 is a local Maxwellian distribution located at the cell interface. It can be determined by the corresponding macroscopic flow variables. a L , a R , and A are related to the derivatives of a Maxwellian distribution in space and time. For details to obtain g 0 , a L , a R and A, see [22,25,26].
With the determination of equilibrium state and the heat flux at the cell interface, the additional term g + in the Shakhov model can be determined. Substituting Equations (6) and (7) into Equation (5), the gas distribution function at the cell interface with particle velocity (u l , v m , w n ) can be expressed as f i+1/2,j,k,l,m,n x j+1/2 , y j , z k , t, u l , v m , w n = 1 − e −t/τ (g 0 + g + ) From the cell interface distribution function we can obtain the distribution function flux and macroscopic flux. We will update the macroscopic variables first with the macroscopic fluxes. Subsequently, we can immediately obtain the local Maxwellian g ζ+1 M and the additional term f +,ζ+1 at ζ + 1 time step inside each cell. Therefore, based on Equation (1) the update of distribution function in UGKS becomes where f f is the distribution function flux across the interface and S is the interface area. The trapezoidal rule has been used for time integration of the collision time. Equation (9) can be rearranged as This is the original explicit UGKS in [22,25].
To accelerate the convergence for steady flow, the authors of [34] introduced the implicit discrete ordinate method for an unstructured physical mesh [12,13] into UGKS. A brief introduction is given below.
Rewriting Equation (1) for f with a particle velocity Treating the loss term of collision integral semi-implicitly and the gain term explicitly we can find where R is the net cell flux averaged over the evolution time step, which can be expressed as The evolution time step ∆ t is determined by the following where ∆t min is the minimum marching time step determined by the stability condition. CFL is the CFL number. Equation (12) can be further written as where the subscript ii indicates the six faces of the physical cell (i,j,k). S i,j,k,ii is the area of the iith face. V i,j,k is the cell volume. The subscript (i1,j1,k1) indicates the cell which shares the iith face with cell (i,j,k). n ii is the outer normal vector of the iith face.
Substituting Equation (16) into Equation (15), we can obtain After a simple deformation, it can be written as Writing in matrix form where (I + ∆t · Z l,m,n ) is a seven-diagonal matrix. NI, NJ, and NK are the total points in the i, j, and k directions of a block in the structured physical mesh, respectively. Applying the LU decomposition yields L l,m,n , U l,m,n are both matrices.
The final form of the implicit UGKS is By performing direct and backward substitutions in a structured physical mesh, (∆ f ) i,j,k,l,m,n can be found. We can then obtain the distribution function f i,j,k,l,m,n at time step ς + 1. After that, macroscopic variables can be obtained with Equations (3) and (4).
The tests [34] on the flows over a cylinder with different free stream Mach numbers showed that the above implicit method can give the same result as the original explicit method with a properly chosen evolving time step. Meanwhile, the computational efficiency can be improved by 1~2 orders.
Due to the explicit treatment of f + i,j,k,l,m,n in the above method, slow convergence exists in small Knudsen number cases. To further accelerate the convergence, Zhu et al. [35] proposed a macroscopic variable prediction technique to deal with f + i,j,k,l,m,n in their implicit UGKS, which is proved to be efficient in all flow regimes.
Under the support of the National Numerical Wind Tunnel Program, an aerodynamic characteristics prediction software applicable for multiple flow regimes called NNW-UGKS [38] has been established, and the viscous flow in the current paper was simulated by this software. Decomposition both in the physical and velocity meshes is applied for MPI parallelism, which is similar to the one in [39]. The composite Newton-Cotes quadrature formula which can be used for any kinds of flow simulation including the current hypersonic or highly non-equilibrium flows, was chosen for integration.
The diffusive reflection wall boundary condition and perfect gas assumption was used.

Results and Modeling
As a demonstrator of the Crew Return Vehicle (CRV), the X38 vehicle has a number of advantages, such as relatively high lift-to-drag ratio and volumetric efficiency [43]. Although the X38 project has long been terminated, research on similar shapes still continues.
The sketch of the vehicle is shown in Figure 1. The reference length of the vehicle, L re f , is 4.67 m. The sketch of the vehicle is shown in Figure 1. Free-stream conditions are given in Table 1. A total number of 24 cases and 4 cases are simulated by viscous and inviscid solvers, respectively. The structured physical mesh is illustrated in Figure 2. The number of cells is 334,434 for altitudes lower than 110 km. For 110 km, the outer boundary is not large enough and additional 82,656 cells were added. The minimum distance near the wall is 1.67 mm which is nearly two times and 0.3% of the free stream mean free paths of 70 km and 110 km, respectively. The velocity mesh is 65 × 65 × 65 and 81 × 81 × 81 for M = 10 and M = 15, respectively, ranging from 2.5 U   to 2.5 U  .
To conduct a thorough mesh convergence for such a problem is almost impossible. As is shown in our previous paper [40] for a 1:16.7 scaled model, good agreements with the DSMC results can be obtained for four different free stream conditions. For the DSMC method, the cell size should be adjusted according to the free stream condition to be smaller than the local mean free path of particles. While for UGKS method, the same physical structured mesh can be used for different free stream conditions. This may be due to the coupling mechanism of the particle transport and collision in UGKS method. The cell size can be larger than the mean free path of particles.  Free-stream conditions are given in Table 1. A total number of 24 cases and 4 cases are simulated by viscous and inviscid solvers, respectively. The structured physical mesh is illustrated in Figure 2. The number of cells is 334,434 for altitudes lower than 110 km. For 110 km, the outer boundary is not large enough and additional 82,656 cells were added. The minimum distance near the wall is 1.67 mm which is nearly two times and 0.3% of the free stream mean free paths of 70 km and 110 km, respectively. The velocity mesh is 65 × 65 × 65 and 81 × 81 × 81 for M = 10 and M = 15, respectively, ranging from −2.5|U ∞ | to 2.5|U ∞ |.   Figure 3 shows the pressure contour of the flow field and the velocity vector on the symmetry plane at two altitudes. For sake of clarity, the grid in the vector diagram is one out of three. The viscous boundary layer can be clearly distinguished from the figure. With the increase in altitude, the shock stand-off distance and the thickness of boundary layer increase, and the wall slip velocity, increases obviously. To conduct a thorough mesh convergence for such a problem is almost impossible. As is shown in our previous paper [40] for a 1:16.7 scaled model, good agreements with the DSMC results can be obtained for four different free stream conditions. For the DSMC method, the cell size should be adjusted according to the free stream condition to be smaller than the local mean free path of particles. While for UGKS method, the same physical structured mesh can be used for different free stream conditions. This may be due to the coupling mechanism of the particle transport and collision in UGKS method. The cell size can be larger than the mean free path of particles.            The local Knudsen number is defined [44] as

Flow Field Characteristics
where l m f p is the local mean free path. The Knudsen number of this form has a great physical meaning. Traditionally, different flow regimes are defined according to the Knudsen numbers [10]. For continuum regime, Kn is smaller than 0.01. For a transitional regime, Kn ranges between 0.01 and 10. When Kn is larger than 10, the flow is considered as free-molecular. Thus, when Kn GLL is much less than unity the flow can be regarded as locally slightly perturbed from equilibrium which is a fundamental assumption of the NS equations. Therefore, it is an appropriate parameter to indicate the degree of non-equilibrium.
where mfp l is the local mean free path. The Knudsen number of this form has a great physical meaning. Traditionally, different flow regimes are defined according to the Knudsen numbers [10]. For continuum regime, Kn is smaller than 0.01. For a transitional regime, Kn ranges between 0.01 and 10. When Kn is larger than 10, the flow is considered as free-molecular. Thus, when GLL Kn is much less than unity the flow can be regarded as locally slightly perturbed from equilibrium which is a fundamental assumption of the NS equations. Therefore, it is an appropriate parameter to indicate the degree of non-equilibrium. Figure 6 shows the local Knudsen number comparison along the y = 500 mm line in front of the vehicle. The local Knudsen number is large inside the bow shock which usually locates in the first peak from left, and near the wall which has been marked on the right. Even at 70 km, the local Knudsen number near the wall and inside the shock is on the order of 0.01, where the continuum assumption may break down. Thus, it is necessary to use UGKS for simulation.     Figure 7 shows the comparison of the pressure distribution on the centerlines. The pressure distribution on the windward centerline shows an increasing trend with the increase in altitude. While on the leeward centerline it increases first and then decreases with the increase in altitude. The magnitude is about an order smaller than that on the windward centerline.  Figure 7 shows the comparison of the pressure distribution on the centerlines. The pressure distribution on the windward centerline shows an increasing trend with the increase in altitude. While on the leeward centerline it increases first and then decreases with the increase in altitude. The magnitude is about an order smaller than that on the windward centerline.      In early research on simple configurations, figures similar to Figure 8 have been frequently given and linear relationships have been obtained. For the current complex vehicle, in the range of 70 km~85 km and X/L = 0.1~0.5, there is a good linear relationship between the pressure change and the viscous interaction parameter on the windward side for 20 degrees angle of attack. In other areas and the whole leeward side, no good linear relationships can be seen. Figure 9 shows the aerodynamic force coefficients computed by the Euler and UGKS solvers. For the UGKS results, the contributions of the pressure and friction are separated. With the increase in the viscous interaction parameter, the axial force and the normal force coefficients increase, and the pressure part and viscous part also increase at the same time. In early research on simple configurations, figures similar to Figure 8 have been frequently given and linear relationships have been obtained. For the current complex vehicle, in the range of 70~85 km and X/L = 0.1~0.5, there is a good linear relationship between the pressure change and the viscous interaction parameter on the windward side for 20 degrees angle of attack. In other areas and the whole leeward side, no good linear relationships can be seen. Figure 9 shows the aerodynamic force coefficients computed by the Euler and UGKS solvers. For the UGKS results, the contributions of the pressure and friction are separated. With the increase in the viscous interaction parameter, the axial force and the normal force coefficients increase, and the pressure part and viscous part also increase at the same time. For the axial force, the viscous part increases rapidly as the altitude increases, from 34% at 70 km to 87% at 110 km. At 80 km and above, the viscous part exceeds the pressure part. For the normal force, the pressure part is dominant, decreasing from 95% at 70 km to 74% at 110 km, and the viscous part is relatively small. between the pressure change and the viscous interaction parameter on the windward side for 20 degrees angle of attack. In other areas and the whole leeward side, no good linear relationships can be seen. Figure 9 shows the aerodynamic force coefficients computed by the Euler and UGKS solvers. For the UGKS results, the contributions of the pressure and friction are separated. With the increase in the viscous interaction parameter, the axial force and the normal force coefficients increase, and the pressure part and viscous part also increase at the same time. For the axial force, the viscous part increases rapidly as the altitude increases, from 34% at 70 km to 87% at 110 km. At 80 km and above, the viscous part exceeds the pressure part. For the normal force, the pressure part is dominant, decreasing from 95% at 70 km to 74% at 110 km, and the viscous part is relatively small.

Aerodynamic Characteristics and Viscous Interaction Modelling
(a) Axial force coefficient (b) Normal force coefficient   Figure 10 shows the viscous force coefficients with the third viscous interaction parameter. Note that the viscous force coefficient is defined as the quantity due to viscous interaction [45] which is equal to the difference between UGKS and Euler solutions. Thus, it is different from the viscous part of UGKS.
rameter. Note that the viscous force coefficient is defined as the quantity due to viscous interaction [45] which is equal to the difference between UGKS and Euler solutions. Thus, it is different from the viscous part of UGKS.
At an altitude less than 100 km where ' From the expressions of A a , A b , and A c , we can conclude that the effect of Mach number can be ignored under the condition of high Mach number. In fact, the effect of Mach number may be mainly reflected in the viscous interaction parameter. Usually, Equation (27) is called a viscous interaction model for the axial force coefficient. Models for other viscous force coefficients can be obtained in a similar way.  At an altitude less than 100 km where ν ∞ is about 0.33 (M = 10), the viscous axial force coefficient has a weak linear relationship with the third viscous interaction parameter. The higher the altitude is, the more serious the deviation from the linear relationship is. In order to correlate the results for all ranges of calculation, it is assumed that the change in the aerodynamic coefficients due to the viscous interaction satisfies the following relationship: As a preliminary study, it is further assumed that (26) According to the calculated aerodynamic force coefficients and the parameters in Equation (26), the following expression of the viscous axial force coefficient can be obtained by fitting with the least square method, From the expressions of a A , b A , and c A , we can conclude that the effect of Mach number can be ignored under the condition of high Mach number. In fact, the effect of Mach number may be mainly reflected in the viscous interaction parameter.
Usually, Equation (27) is called a viscous interaction model for the axial force coefficient. Models for other viscous force coefficients can be obtained in a similar way. Figure 11 shows the correlation curve between the viscous interaction model prediction data which is represented by suffix '_ model' and the numerical simulation data which is represented by suffix '_ UGKS'. The data are basically distributed near the correlation line at different angles of attack and Mach numbers. It can be seen that the correlation between the data is good. tropy 2022, 24, x FOR PEER REVIEW Figure 11 shows the correlation curve between the viscous interaction m tion data which is represented by suffix '_ model' and the numerical sim which is represented by suffix '_ UGKS'. The data are basically distributed relation line at different angles of attack and Mach numbers. It can be seen th lation between the data is good. The Pearson Correlation Coefficient r, which is widely used in statistics characterize the degree of correlation between the aerodynamic prediction numerical simulation data, and its expression is    n x x y y   Figure 11. Correlation between viscous interaction model and numerical simulation results. The Pearson Correlation Coefficient r, which is widely used in statistics, is chosen to characterize the degree of correlation between the aerodynamic prediction data and the numerical simulation data, and its expression is The closer r is to 1, the better agreement between the predicted values of the model and the numerical results we obtain. The Pearson correlation coefficients of axial force, normal force and pitching moment are 0.999996, 0.999973, and 0.999863, respectively. They are all very close to 1, indicating that the correlation between the predicted data and the numerical simulation data is very good.
In order to further assess the viscous interaction model, the relative orthogonal distance, dr i , is defined to characterize the relative degree of deviation of the data from the correlation curve, as shown in the following, The dr i of the viscous axial force is shown in Figure 12. The maximum fitting deviation is only 1.8%. Finally, the accuracy of the prediction model is preliminarily evaluated. The UGKS and the viscous interaction models are used to calculate two new cases with altitude equal to 80 km and 90 km, respectively. The angle of attack is 30 degrees with a Mach number of 15. The results and relative errors are shown in Table 2. The relative error o viscous axial force is small partially due to its large magnitude. While the error of viscou pitching moment is large due to its small magnitude compared with the viscous axia force. However, the relative error of the pitching moment itself is small. Taking the 80 km case as an example, the relative error of predicted viscous pitching moment is 9.87%. How ever, the pitching moments obtained by UGKS simulation and predicted by the model are -0.2158 and -0.2180, respectively, resulting in a relative error of only 1.01%.  Finally, the accuracy of the prediction model is preliminarily evaluated. The UGKS and the viscous interaction models are used to calculate two new cases with altitudes equal to 80 km and 90 km, respectively. The angle of attack is 30 degrees with a Mach number of 15. The results and relative errors are shown in Table 2. The relative error of viscous axial force is small partially due to its large magnitude. While the error of viscous pitching moment is large due to its small magnitude compared with the viscous axial force. However, the relative error of the pitching moment itself is small. Taking the 80 km case as an example, the relative error of predicted viscous pitching moment is 9.87%. However, the pitching moments obtained by UGKS simulation and predicted by the model are −0.2158 and −0.2180, respectively, resulting in a relative error of only 1.01%.

Conclusions
Hypersonic viscous and inviscid flow fields around the X38-like vehicle are simulated by UGKS solver and Euler solver, respectively. Viscous force coefficients at different altitudes, Mach numbers, and attack angles are obtained by subtracting the two solutions and correlated by the third viscous interaction parameter. A nonlinear viscous interaction model of force coefficients is established, and some preliminary conclusions are as follows, (1) For the X38-like vehicle, the contribution of the viscous part to the axial force coefficients increases rapidly with altitude, and reaches 87% at 110 km for the typical conditions, with Ma = 10 and AOA = 20. The contribution of the viscous part to the normal force coefficients is small, and can only reach 26% at 110 km.
(2) For complex configurations such as the current X38-like vehicle, the changes of wall pressure and aerodynamic coefficients due to viscous interaction cannot be expressed linearly with the viscous interaction parameters in the whole flow field.
(3) A viscous interaction model can be established by taking the viscous interaction parameters as the independent variables combined with the inviscid solution and the viscous solution, which is helpful to quickly obtain the aerodynamic characteristics at moderate to high altitudes and has certain application value in engineering design.
In this paper, the idea of modeling the viscous interaction based on UGKS solver is applied to the X38-like vehicle, and a satisfactory result has been achieved. The prediction model can take into account both the viscous interaction effect and rarefied gas effect. However, the Cartesian velocity mesh in our UGKS solver causes huge waste both in computation and memory. The next step is to introduce an unstructured velocity mesh into our solver to reduce the cost and give a more accurate prediction model for more complex configurations.