Jiwei LI, Jingfeng WANG,*, Liming YANG, Chng SHU

a College of Aerospace Engineering, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China

b Department of Mechanical Engineering, National University of Singapore, Singapore 119260, Singapore

KEYWORDS Finite volume method;Fluid-structural-thermal;Hybrid lattice Boltzmann flux solver;Hypersonic;Improved switch function

Abstract In this paper, a hybrid Lattice Boltzmann Flux Solver (LBFS)with an improved switch function is proposed for simulation of integrated hypersonic fluid-thermal-structural problems. In the solver, the macroscopic Navier-Stokes equations and structural heat transfer equation are discretized by the finite volume method,and the numerical fluxes at the cell interface are reconstructed by the local solution of the Boltzmann equation.To compute the numerical fluxes,two equilibrium distribution functions are introduced. One is the D1Q4 discrete velocity model for calculating the inviscid flux across the cell interface of Navier-Stokes equations, and the other is the D2Q4 model for evaluating the flux of structural energy equation. In this work, a new dual thermal resistance model is proposed to calculate the thermal properties at the fluid-solid interface. The accuracy and stability of the present hybrid solver are validated by simulating several numerical examples,including the fluid-thermal-structural problem of cylindrical leading edge. Numerical results show that the present solver can accurately predict the thermal properties of hypersonic fluid-thermalstructural problems and has the great potential for solving fluid-thermal-structural problems of long-endurance high-speed vehicles.

1. Introduction

Hypersonic vehicles experience severe aero-thermal loads when they fly through the atmospheric environments.1-3The accurate prediction of aero-thermal loads is especially important to design the lightweight structures and the Thermal Protection Systems(TPS)for high-speed vehicles,because the choices of material and wall thickness mainly depend on the accurate estimation of temperature distribution in the vehicles structure.4Therefore,the Conjugate Heat Transfer(CHT)5between the fluid flow and the structural heat conduction is still an urgent issue to be solved for the prediction of aerodynamic thermal load and its effect on hypersonic vehicle. Because of the difficulty and high cost of wind tunnel tests,the numerical simulations have recently become the major tool to evaluate the hypersonic fluid-thermal-structural performance of the TPS.

At present, the traditional partition-coupled numerical approaches are widely used to simulate the hypersonic fluidthermal-structural problems. In these methods,6-9the external high-speed flow and internal thermal-structure are separated into two individual parts for independent simulations. Then,the interdisciplinary interactions and iterative process between the fluid and solid domains are generally performed through a two-way data transmission strategy of temperature and heat flux at the fluid-solid interfaces. Therefore, these partitioncoupled approaches can be easily implemented by utilizing the available existing numerical computation solvers of fluid dynamics and structural dynamics.

However, with the development of computational technology,now increasingly more researchers are coming to the conclusion that the partition-coupled method cannot accurately simulate the coupling characteristics of fluid-thermalstructural multi-physical fields, because the hypersonic CHT problem is a continuous and impartible physical process in both spatial and time domains. In Wieting’s researches,10,11it is especially emphasized that the interfacing partitioncoupled method is relatively imprecise and unreliable, which is attributed to the incompatible numerical calculation models between the different solvers. These mathematical models not only require series of coupled iteration strategies, such as the tightly coupling and loosely coupling strategies, but also preprocess the request of interfacing data transfer between the unmatched meshes of flow and solid domains. None of the above extensive processing can satisfy the physical change conditions of the whole continuous multi-physics,which will result in the reduction of computation efficiency and accuracy.Therefore,it is highly essential and desirable to integrate these multidisciplinary coupling methods into a unified fluidthermal-structural solver.In addition,various upwind schemes have been developed to evaluate inviscid flux in simulation of compressible flows in traditional Computational Fluid Dynamics (CFD) methods, such as Roe scheme and AUSM scheme. However, for the simulation of hypersonic flows, the Roe scheme may exhibit carbuncle phenomenon and induce numerical instability. And sometimes, in AUSM scheme, the numerical oscillations occur distinctly in shock regions.Consequently, how to implement this unified algorithm and overcome the drawback of the existing upwind schemes motivates the present work.

In recent years,the Lattice Boltzmann Method(LBM)12-14has drawn increasing attention as an alternative CFD approach,due to its unique advantages of kinetic nature, simplicity and scalability. Nonetheless, at present, the simulation of high-Mach number flows is still a challenge for LBM on account of the limitations of numerical schemes and compressible lattice Boltzmann models.So far,a number of studies have focused on developing compressible lattice Boltzmann models and some promising models have been proposed.For instance,Qu et al.15developed a series of compressible lattice Boltzmann models by using a simplified circular function to replace the Maxwellian distribution function. These models were applied successfully in simulating some strong shock waves flow. Dellar16also derived some compressible lattice Boltzmann models including one-dimensional(1D)unsplit and split seven-velocity models. However, the lattice velocities in these models are required to be artificially set in advance, which may not be suitable for hypersonic flows.

In view of this, Yang et al.17developed a moment conservation-based platform to derive non-free parameter compressible lattice Boltzmann model. The macroscopic Navier-Stokes equations can be locally recovered by using the compressible lattice Boltzmann model obtained from this platform. In this work, a non-free parameter D1Q4 model developed by Yang et al.18will be locally applied to compute the numerical flux at the cell interface in simulating hypersonic flows.

Traditionally, to simulate compressible flows by LBM, the Discrete Velocity Boltzmann Equation(DVBE)-based method is used to update the distribution functions in time.Qu et al.,19Li et al.,20Ubertini et al.21and Stiebler et al.22have made some efforts on this method.But it should be noted that the process of directly solving DVBE is tedious and the computational efficiency is relatively low. To develop a more efficient solver for compressible flows with the good properties of LBM,Shu and his coworkers17,18,23-26proposed the Lattice Boltzmann Flux Solver (LBFS). In LBFS, the conventional Finite Volume Method(FVM)is applied to discretize macroscopic governing equations, and the local solution of 1D compressible lattice Boltzmann model is used to reconstruct the inviscid flux at the cell interface. The good positive property and efficiency of the LBFS in simulating compressible flows with strong shock wave have been demonstrated. In the existing LBFS,the non-equilibrium component of the distribution function is used to introduce numerical dissipation for evaluating inviscid flux at the cell interface, which is important for capturing strong shock wave. However, the introduced numerical dissipation is too large to disturb the true solution in smooth regions such as in boundary layers, especially for accurately predicting the heat flux in hypersonic flow regime. To both capture strong shock and predict the heating rate accurately,a modified switch function will be proposed to control the numerical dissipation in the present work.

On the other hand,the LBM can also be applied in the thermal flow. A number of promising efforts have been made to solve this problem, such as, the Thermal Lattice Boltzmann Method (TLBM)27-30and the Thermal Lattice Boltzmann Flux Solver(TLBFS).31,32Inspired by these meaningful works,the application of LBM will be extended into structural heat conduction calculation in this paper. Compared with the TLBM, the TLBFS has been proven to be more efficient and can be applied to the non-uniform meshes, which is used popularly.Similar to the above LBFS,in TLBFS,the FVM is used in spatial discretization, and the numerical fluxes at the cell interface are reconstructed locally from the solution of the thermal lattice Boltzmann model. Therefore, the LBFS and TLBFS can be combined into one set of unified method for simulation of integrated hypersonic fluid-thermal-structural problems in this paper. In addition, a number of thermal LBE models are developed for incompressible thermal flow.27,28,33-35In these models, the D2Q4 model developed by Guo et al.35is relatively efficient and simple. And it is demonstrated that the D2Q4 has excellent numerical stability and accuracy by Huang et al.36and Wang et al.37Hence, the D2Q4 model will be adopted into the TLBFS to solve the structural heat conduction equation in the present work.

To sum up, the purpose of this paper is to develop a new unified solver for integrated hypersonic fluid-thermalstructural analysis based on a hybrid lattice Boltzmann flux solver. In this approach, the FVM is applied to discretize all the macroscopic equations. For the Navier-Stokes equations,the inviscid flux at the cell interface is calculated by using LBFS with D1Q4 model, and the viscous flux is calculated using conventional second-order central scheme.For the structural heat conduction equation,the TLBFS with D2Q4 model is used to evaluate the numerical flux. In addition, a new twothermal-resistance model is proposed to compute the thermal parameters for satisfying the continuity condition of the heat flux across the fluid-solid interface. To verify the feasibility and accuracy of the developed unified solver,several numerical examples are simulated and the simulation results are compared with the experimental data38and the computed results in the available references.

2. Fluid-thermal-structural equations

2.1. Aerodynamic flow

In this present work, the Reynolds average Navier-Stokes equations are solved by using the FVM based on cell-centred scheme for the compressible flow.In general,the integral form of Navier-Stokes equations without source term can be written as

To solve Eq. (1), the key process is to evaluate the convective fluxes Fcand diffusive fluxes Fv. In accordance with the elliptic property of the viscous flux, Fvis usually estimated by using smooth functions such as polynomial interpolation function.Therefore,the only remaining problem is to calculate the inviscid flux Fc.

2.2. Structural heat conduction

Generally, the macroscopic energy equation can be used to describe the thermal response of the structure as

where u and T are the velocity and temperature respectively.Within the solid domain, u can be fixed at zero. κ represents the thermal diffusivity and can be written as

where ρs,csand ksare the structural density,specific heat and thermal conductivity respectively.

3. Solution procedure

3.1. Fluid analysis

3.1.1. Relationship between Navier-Stokes equations and nonfree parameter D1Q4 lattice Boltzmann model

To effectively calculate the inviscid flux of Navier-Stokes equations, the LBFS is applied in this work. In LBFS, the FVM is applied to discretize the macroscopic Navier-Stokes equations, and the inviscid flux at the cell interface is reconstructed using the local solution of 1D compressible lattice Boltzmann model. The 1D compressible lattice Boltzmann model used in this work is the non-free parameter D1Q4 model, which is developed to simulate the hypersonic flows with strong shock waves.

As mentioned above,the viscous flux term Fvis usually calculated using second-order central difference method, which will not be described in detail in this article. In this work,the main step is evaluating the inviscid flux Fcby utilizing the LBFS with D1Q4 model.In LBFS,the non-free parameter D1Q4 model is applied into 1D velocity space at the cell interface. Therefore, the vectors of conservative variables Wand inviscid fluxes Fcin Eq.(2)need to be converted into an alternative expression with the normal and tangential velocities.At the cell interface, the vector of tangential velocity can be obtained from

From the physical conservation relations in recovering Euler equation, the vectors of conservative variables W and inviscid fluxes Fcalong the normal direction can be written as

where ep= [（3 - γ）/2]e is the particle potential energy. And fi（r ，t） is the distribution function at the cell-interface, which will be presented subsequently in detail.

Fig. 1 Discrete lattice velocities distribution of non-free parameter D1Q4 model.

From Eq. (13), we can see that the evaluations of the conservative variables W and inviscid fluxes Fcare also related with the tangential velocity, which cannot be determined directly using the LBFS with D1Q4 model. How to calculate approximately the contributed flux of the tangential velocity will be shown in the following section.

3.1.2. Evaluation of inviscid flux by local solution of D1Q4 lattice Boltzmann model

According to the Chapman-Enskog analysis,to evaluate inviscid flux, only the equilibrium part of the distribution function needs to be determined in principle and the non-equilibrium part should be taken as zero, which can provide the accurate solutions for the smooth flow domains with very little numerical dissipation. However, this scheme cannot be applied directly to capture the strong shock waves with stable solution.Consequently,both the equilibrium and non-equilibrium parts of the distribution function need to be considered at the cell interface in the present LBFS. The contribution of the nonequilibrium part will be actually retained straightforwardly as the numerical dissipation. As a result, the full distribution function at the cell interface located at r = 0 has the following expression:

where τ0=τ/δt is the dimensionless collision time and δt is the streaming time step. As previously stated, fineq（0 ，t） is kept as the numerical dissipation in the present LBFS. Thus,τ0can be viewed as the weight coefficient of the numerical dissipation.

By substituting Eq. (16) into Eq. (11), the inviscid flux at the cell interface in the normal direction can be expressed as

The last undetermined variable in Eq. (16) and Eq. (17) is the dimensionless collision time τ0, treated as the weight coefficient of the numerical dissipation. Thus, a switch function is constructed to determine the value of τ0for the control of the numerical dissipation. Inspired by the construction of limiter of Barth and Jespersen, Yang et al.25proposed a switch function of pressure based on hyperbolic tangent function to manage the numerical dissipation at the interface, which is defined as

where the subscripts‘‘L”and‘‘R”denote the left and the right grid cells respectively. PL, PRare the pressures at the left and right sides of cell interface. C is an constant amplification factor. NfLand NfRare the number of the faces of the left and right grid cells respectively. In this switch function, the hyperbolic tangent function tanh （x ） is used to calculate the dimensionless collision time τiat the interface, which can keep the value of τ0changing smoothly from 0 to 1 with the pressure difference of the left and right cells.

In Yang’s numerical study, the switch function has good performance in capturing strong shock waves and boundary layer for simulation of high-speed flows. However, this switch function is designed based on the uniform grids. In practical applications, computational meshes are rarely uniform, and especially,the boundary layer meshes have a large aspect ratio.

Fig. 2 Schematic diagram of uniform and non-uniform grids with different aspect ratios.

Fig. 2 illustrates the schematic diagram of uniform and non-uniform grids with different aspect ratios. If the same pressures are considered, Yang’s switch function may have the result of the same values of τ0for the uniform and nonuniform grids. Nevertheless, it is obvious that the values of τ0at the interface in the uniform and non-uniform grids should be different. So, this switch function is needed to be modified for the non-uniform grid.

Furthermore, the accuracy of calculating wall heat flux depends on the accurate prediction of both the shock wave and thermal boundary layer, which is very sensitive to the numerical dissipation. According to the kinetic theory, the dimensionless relaxation time τ0depends on both the pressure and the temperature.Therefore,based on the above analysis,a new improved switch function is proposed to precisely control the numerical dissipation, which is defined as

3.2. Structural heat transfer analysis

3.2.1. D2Q4 model for energy equation

To calculate the heat transfer equation, the local solution of LBE with the D2Q4 model proposed by Guo et al.35is utilized to calculate the numerical flux of energy Eq. (3) in this work.Fig. 3 shows the discrete particle velocities distributions at a cell interface in the non-free parameter D2Q4 model.

It can be seen that all the discrete velocity points are located on the coordinate axis.This makes the model be easily applied to evaluate the numerical flux of the energy equation.

From the Chapman-Enskog analysis and Yang’s work,41Eq. (3) can be recovered from the discrete velocity Boltzmann equation expressed as

where hαis the temperature distribution function and heqαis the equilibrium state approached by hαthrough particle collisions within a collision time scale τκ. eαis the vector of the particle velocity. In the non-free parameter D2Q4 model, the equilibrium distribution function heqα and the corresponding discrete particle velocity eαcan be written as

Fig. 3 Discrete lattice velocity distribution of D2Q4 model at cell interface.

where Δ l and Δ r are the shortest edge length of the left and right cells around the cell interface respectively.

3.3. Corresponding thermal relationship at fluid-solid interface

To satisfy the continuity requirement on the heat flux, the basic physical parameters (the temperature, the temperature gradient and the thermal conductivity) must be redefined at the fluid-solid interface due to the significant differences in thermal physical properties between the fluid field and solid field. This is an important step to accurately predict the wall heat flux in the CHT problem.

To clearly describe the coupling relationship between the fluid domain and solid domain, Fig. 4 illustrates the corresponding relationship between the left cell and right cell at the fluid-solid interface, where Ωfand Ωsrepresent the fluid and solid domains respectively.

As we know, the evaluation of temperature gradient is directly related to the calculation of heat-flux. In order to improve the calculation accuracy of temperature gradient, a modified method based on the Green-Gauss42approach to the evaluation of gradients is applied in this work, i.e.,

where ?Tland ?Trrepresent the temperature gradients approximated by Green-Gauss theorem at the centroids of the left and right cells respectively; Trand Tlare the temperatures at the centroids of the right and left cells respectively; rlrand Llrdenote the unit vector and distance between the cellcentroids of the left and right cells respectively.

Based on the continuity condition of heat flux, the temperature at the fluid-solid interface can be given by

where kland krdenote the thermal conductivities of the fluid domain (left cell) and solid domain (right cell) respectively.

In addition, the thermal conductivity is also one of the important parameters for evaluating the heat flux. Therefore,it is a critical step to calculate the thermal conductivity accurately at the fluid-solid interface for simulating CHT problem.Inspired by the concept of structural thermal resistance, a dual-thermal-resistance model is proposed to calculate the equivalent thermal conductivity at the interface for simulating the interphase heat transfer process accurately.

In general, the thermal resistance can be written as

where A is the heat transfer area,d is the heat transfer distance,and k is the thermal conductivity. The dual thermal resistance between the centroids of the left and right grid cells can be expressed as From Eq.(44)and Eq.(45),the following expression can be obtained:

where klris interface equivalent. Lland Lrare the distances from the centroids of the left and right cells to the interface respectively. Thus, the equivalent thermal conductivity at the fluid-solid interface can be given by

In order to explicitly describe the present hybrid solver,the detailed calculating procedure is illustrated in Fig. 5. Compared with the traditional partitioned approach, we can see that no additional data-exchange strategy is required in this solver, which will improve the calculation efficiency for simulating CHT problems.

4. Numerical examples

4.1. Case 1: Hypersonic flow around a circular cylinder

To investigate the capability of LBFS for simulation of hypersonic flows with strong shock waves and boundary layer, the hypersonic flow around a circular cylinder is simulated. This case is taken from the experiment done by Wieting,38where the freestream condition for the air is listed in Table 1.

Fig.4 Schematic diagram of coupling relationship at fluid-solid interface.

In this case, structured mesh with 320×320 cells is used,and the first-cell spacing sizes normal to the wall are determined by the cell Reynolds number43taken as Recell=ρ∞U∞-Δn/μ∞=1.835,where Δn is the mesh spacing of the first cell in the normal direction adjacent to the cylinder surface; ρ∞, U∞and μ∞are the density, velocity and dynamic viscosity of the freestream respectively.

Fig.6 shows the pressure,temperature and switch function contours of hypersonic flow around a circular cylinder obtained from LBFS. We can see that the strong shock can be captured stably without any oscillation in the post-shock region in Fig. 6(a) and (b). In addition, the switch function plays a desirable role well. The switch function is taken as a value close to zero in smooth regions and one around the strong shock wave in Fig. 6(c), as the expected theory.

Fig.7 presents the predicted pressure p/p0(upper)and heat flux q/q0(bottom) distributions normalized to the stagnation point values along the cylindrical surface. Also displayed in this figure are the experimental data,38and results calculated by AUSM+scheme44and gas-kinetic BGK scheme.45It can be observed that the pressure and heat flux distributions predicted by LBFS are compared well with these reference values.

In addition, Table 2 compares the computed pressure and heat flux at stagnation point with the relevant references. As shown in this table, the present prediction results agree well with the solutions computed by Xu et al.45and Yang et al.18This numerical example shows that the present LBFS can well capture strong shock waves and thin boundary layer of hypersonic viscous flows, so its accuracy is verified to be acceptable for the subsequent analysis.

4.2. Case 2: Thermal conduction for 2D-rectangle domain

A good benchmark case for investigating the capability and accuracy of present method for simulation of structural heat conduction is the thermal conduction for 2D-rectangle domain.There is an exact solution for this numerical example.The length and width of the two-dimensional rectangle are a and b respectively. The temperature boundary conditions are listed as follows:

Fig. 5 Calculation procedure of hybrid flux solver.

Table 1 Initial freestream conditions.

For the convenience of calculation, in this case, the geometry sizes of the 2D-rectangle domain are set as a=1 m and b=1 m.The thermal properties are set as T0=1 K,ρs=1-kg/m3, cs=1 J/(kg.K), ks=1 W/(m.K). In this simulation,both unstructured mesh and uniform structured grid with Nx-×Ny=20×20 cells are used.

Fig. 8 presents the comparative temperature contours between exact solution and numerical solutions using present method with unstructured and structured meshes respectively.It can be observed that the temperature contours obtained from the present scheme are in good agreement with the exact solution. In order to better observe the temperature variation,the predicted temperature distribution along the lines(Y=0.4 m and Y=0.8 m) within the rectangle domain is illustrated in Fig.9.It can be clearly seen that good agreements between numerical solutions and exact solution are achieved.

To study the effect of grid scale on calculation accuracy,the computational domain is divided by four different uniform grids Nx×Ny=20×20, 40×40, 60×60, and 80×80.Two points (Point A and Point B) are selected randomly to compare the detailed numerical error between exact value and the numerical solutions with the different meshes, as shown in Table 3. The coordinate of Point A is (0.9872,0.8123) and located in the boundary cell of the 2D-rectangle domain. The coordinate of Point B is (0.8875, 0.1875) and located within the internal cell.

From Table 3, it can be clearly seen that good agreements between the exact value and the numerical solutions predicted by the present method are achieved. And, grid refinement can contribute to reducing the numerical error. Given the above,the performance of the developed scheme is verified by its application to the thermal conduction within 2D-rectangle model and can be acceptable for use in the subsequent integrated fluid-thermal-structural analysis.

4.3. Case 3: Integrated fluid-thermal-structural analysis for cylindrical leading edge

Fig. 6 Pressure, temperature and switch function contours of hypersonic flow around a circular cylinder.

Fig.7 Comparison of pressure and heat flux distributions along cylindrical surface.

Table 2 Comparison of stagnation pressure and heat flux.

In this section, the solution to Mach number 6.47 flow over a cylinder is used to validate the accuracy of the hybrid lattice Boltzmann flux solver and study the fluid-thermal-structural behavior of aerodynamically heated leading edges. In 1987,the aerodynamic heating experiment of a stainless-steel cylinder was performed in NASA Langley High Temperature Tunnel by Wieting.38This experiment has been repeatedly utilized to verify the partition-coupled approaches for CHT analysis by a lot of researchers.11,46-50In the experiment, the cylinder is made of 321 stainless steel (1Crl8Ni9Ti) with outer radius Rout=0.0381 m and inner radius Rout=0.0254 m. The temperature-independent thermal properties of 321 stainless steel are listed in Table 4. In addition, the initial freestream conditions for the hypersonic flow are given in Table 5.

Fig. 10 shows the computational grid and thermal boundary conditions of the cylindrical leading edge. An integrative computational grid of the fluid domain (green zone) and solid domain(red zone)is used in the hybrid lattice Boltzmann flux solver for the integrated fluid-thermal-structural simulation.The mesh in the regions of shock wave and boundary layer is locally refined for sufficient resolution.

To accurately predict the wall heat flux without extra calculation burden, the normal first-cell spacing sizes at the fluidsolid interface within both the flow and cylinder domains are set as 1×10-6m with the corresponding cell Reynolds number43Recell=1.13. The extremely small grid-spacing size is able to capture the steep temperature gradient within the thin thermal layer.

4.3.1. Unsteady results

According to the physical characteristic time scale51of hypersonic flow, a fixed real-time step Δt=1×10-4s is used in the unsteady solution.

A major difficulty encountered in the present integrated fluid-thermal-structural analysis is to accurately predict the aerodynamic heating rates. This is because a good accurate resolution of the strong shock wave and high temperature gradient within the thin thermal layer at the fluid-solid interface is highly required. As mentioned above, an improved switch function is proposed to precisely control the value of relaxation time τ0for calculating wall heat flux accurately.

To exhibit the improved performance of the new switch function, Fig. 11 displays the switch function contours obtained from Yang’s method18and present work.We can find that there exists the significant difference between Fig. 11(a)and (b) in shock expansion region. In theory, the value of relaxation time τ0should be taken arbitrarily close to zero in the smooth region and one in the region of strong shock wave.Hence,the performance of Yang’s switch function is not good enough to control the value of τ0with the non-uniform grids.But this expectation can be well achieved using the improved switch function, as shown in Fig. 11(b). This is because, after the modification, the relaxation time τ0depends on both pressure and temperature. Furthermore, a control parameter δ is introduced into the improved switch function to reduce the influence of excessive aspect ratio of computational mesh.For high-speed flow, the normal first-cell spacing size at the wall within the boundary layer is usually extremely small to capture the temperature gradient and hence the aerodynamic heating rate,which leads to the negative effect of too large grid aspect ratio on the relaxation time τ0.

Fig. 8 Comparison of temperature contours within 2D-rectangle domain.

Fig. 9 Comparative temperature distributions between exact and computed solutions.

Table 3 Comparison of numerical error between exact value and numerical solutions with different meshes.

Table 4 Material thermal properties of cylinder.

Table 5 Initial freestream conditions.

In addition,Fig.12 shows the comparative surface heat flux distributions obtained from the present work, Yang’s method and AUSM+ scheme. As can be seen, the present heat flux distributions are in excellent agreement with the AUSM+scheme solution. But the result of Yang’s method is definitely over-high in the shock expansion zone, which is attributed to the large numerical dissipation introduced by the unmodified switch function.However,the difference will be reduced when the computational grid is enough uniform and refined.

Fig.10 Integrated computational grid and boundary conditions.

Fig. 11 Switch function contours for hypersonic flow around cylinder.

As illustrated in Fig.13,the predicted heating rate distribution normalized to the stagnation point heating rate value at time zero coincides well with the available experimental results.In this work, the predicted maximum heating rate value at stagnation point is 492.4 kW/m2, which is in excellent agreement with the results of the Fay-Riddell solution52(482.6 kW/m2) and a viscous shock-layer solution53(470.2 kW/m2). In order to adequately illustrate the accuracy of the calculated stagnation-point heating rate,Table 6 lists the computed values of stagnation-point heating rate at t=0 s in available literature.As we can see,the present predicted result is also compared well with the numerical solutions in the relevant references. However, all these numerical results of stagnation-point heating rate are lower than the experimental data38(670 kW/m2). This is because the aerodynamic heating rate is very sensitive to the high-speed flow turbulence,but the influence of the freestream turbulence on the aerodynamic heating rate in wind tunnel experiment is not considered in the numerical simulations.

Fig. 14 shows the time evolution of the predicted stagnation-point heating rate and its rate of change with time.As time advances, the stagnation-point heating rate decreases gradually, and its rate of change with time gradually approaches to zero after the rapid increase within a short initial period.These results can be explained as follows.First,the external heat generated by aerodynamic heating is continuously transmitted into the inside of the cylindrical body across the fluid-solid interface. Then, the wall temperature of the cylinder increases gradually, and the steep temperature gradient in the thin thermal layer of the stagnation point region decreases gradually, which leads to the gradual decline of the stagnation-point heating rate. Within two seconds, the stagnation-point heating rate is reduced by about 6.8% ,which is close to Dechaumphai’s11calculated result (approximately 8% ).

Fig. 12 Comparative heat flux distributions at t=0 s along cylinder surface.

Fig. 13 Comparison between present and experimental heating rate distributions along cylinder surface.

Fig.15 depicts the time evolution of computed temperature distribution along the symmetry-line (Y=0 m) within the fluid and cylinder domains. As can be seen, the flow temperature rises sharply from 241.5 K to about 2262.6 K after being heated by the strong bow shock wave. Then, the temperature of the heated flow rapidly drops to the initial wall temperature of 294.4 K at t=0 s within a very thin thermal boundary layer. This steep temperature gradient within the thin thermal layer, which is about 3.1% of the shock-layer thickness, produces a high stagnation-point heating rate.Over time,the flow temperature distribution along the centerline remains almost unchanged,only except within the thin thermal boundary layer near the fluid-solid interface. Since the aerodynamic heat is absorbed into the interior of the cylindrical body, the temperature distribution along the centerline gradually rises with time,and the stagnation-point temperature is always the highest. In addition, it is noted that the predicted bow shock wave position is located at X=-54.5 mm, which is compared well with the location (X=-54.4 mm) calculated by empirical formulae.54

In order to intuitively observe the continuous process of aerodynamic heat conduction, Fig. 16 illustrates the temperature field contours within the cylinder body at t=0.1 s,t=0.5 s,t=1.0 s and t=2.0 s.It can be seen obviously that,as time passes, the cylinder gets heated continuously by the external hypersonic flow.With the aerodynamic heat transferring gradually into the depth of cylinder body, the high temperature zone area expands continuously from the stagnation point region.

Fig. 17 shows the evolution of temperature distribution along the cylinder surface within the time period of 0 - 2.0 s.Accordingly, the time evolution of the predicted stagnationpoint temperature and its rate of change with time are shown in Fig. 18. As we can see, the predicted temperature distribution along the cylinder surface increases gradually with time and the increment of stagnation-point temperature is the largest, which is consistent with the change trend shown in Fig. 16. At t=2.0 s , the predicted maximum temperature ofabout 389.2 Koccurs at the stagnation point,which is in good agreement with Dechaumphai’s11calculated result of 388.8 K.Moreover, the comparison of stagnation-point temperature at t=2.0 s with the available references is also presented in Table 6. As shown in Fig. 18, the time evolution of stagnation-point temperature from 0 s to 2.0 s is compared well with the predicted results in the related references,49,50and the maximum temperature discrepancy at t=2.0 s is about approximately 3.2 K. It can also be seen that the stagnation-point temperature quickly increases at the beginning as a result of the initial severe aerodynamic heating,and then its rate of change with time also gradually tends to zero. This is because the gradually rising wall temperature causes the decrease of temperature gradient within the thermal layer, which weakens the external aerodynamic heating rate and further results in a gradual decline in stagnation-point temperature rise rate. The whole physics process accords with the real evolution of CHT problem.

Table 6 Comparison of stagnation-point temperature(t=2.0 s)and stagnation-point heating rate(t=0 s)with available references.

Fig. 14 Time evolution of predicted stagnation-point heat flux and change rate with time.

Fig.15 Temperature distribution with time along symmetry-line within fluid and solid domains.

To sum up,the present predicted results are in good agreement with those solutions of the traditional partitioned methods in the available references. Thus, the accuracy and reliability of the present integrated hybrid lattice Boltzmann flux solver are verified to be acceptable for simulation of hypersonic fluid-thermal-structural analysis.

4.3.2. Steady results

In the following, the steady solution of integrated fluidthermal-structural analysis for the cylindrical leading edge is performed and discussed. It should be noted that this present steady analysis is to simulate the thermal equilibrium state between the external hypersonic aerodynamic heating and the heat transfer within the cylinder body.

Fig.19(a)and(b)show the steady temperature contours of the hypersonic flow and the cylinder body respectively. And both the predicted steady temperature and heating rate distributions along the cylinder surface are illustrated in Fig.20.As can be observed from these figures, when the steady thermal equilibrium state is reached in the CHT phenomenon,the maximum temperature in the flow field is approximately 2264 K,which almost remains unchanged compared to the above unsteady results; the predicted maximum temperature within the cylinder body reaches about 648 K at the stagnation point.From the outer wall to the inner wall, the predicted temperature within the cylinder domain gradually decreases to 294.4 K owing to the thermal boundary of isothermal inner wall. In the steady simulation, the maximum temperature rise of the cylinder body is about 353.6 K from the initial value of 294.4 K. From Fig. 20, we can see that the wall heating rate and temperature along the fluid-solid interface consistently reduce from the stagnation point to the two ends.At the stagnation point, the maximum steady heating rate remains 327.7 kW/m2, which means that the aerodynamic heat constantly transfers into the cylinder body in the present steady state.

The comparison of the computed density contour (lower half)with the experimental schlieren photograph11(upper half)is shown in Fig.21(a).It can be clearly seen that the predicted shock shape and position are compared well with the experimental schlieren photograph, which is consistent with the results shown in Fig. 15. Therefore, the rationality and effectiveness of the assumption of a calorically perfect gas in the present work are validated by these comparisons. In addition,the excellent agreement between the predicted pressure distribution normalized stagnation-point pressure and the experimental data along the cylinder surface is also shown in Fig. 21(b). The relative error between the computed stagnation-point pressure and experimental result is within 3.9%. The differences between the calculated results obtained by Dechaumpai et al.11and Guo et al.47and the experimental data are approximately 5% and 7.08%respectively, so the predicted pressure presented in this study is correct and reliable.

Fig. 16 Temperature evolution within cylinder body.

Fig. 17 Temperature distribution with time along cylinder surface.

Fig. 19 Steady temperature contours within fluid and cylinder domains.

Fig. 20 Steady temperature and heat flux distributions along cylinder surface.

Fig.21 Comparison between predicted and experimental results.

In conclusion,the feasibility and accuracy of the hybrid lattice Boltzmann flux solver for simulating steady hypersonic CHT problem are also shown in this steady example. It is worth mentioning that the outstanding advantage of the present hybrid solver is that the steady temperature and heating rate distributions within both the fluid and solid domains can be obtained fast, which has the great potential for efficiently solving fluid-thermal-structural problems of longendurance high-speed vehicles. By contrast, the traditional interfacing partition-coupled methods are not applicable for solving these steady problems. A large number of alternating iterations in time domain between the independent interdisciplinary solvers with additional interfacing data transfer across the fluid-solid interface are required to achieve the thermal equilibrium state, which is uneconomical and inefficient for the steady CHT analysis. Therefore, the proposed integrated solver can effectively reduce the computing time for the steady case by eliminating the extra interdisciplinary coupling strategy.

5. Conclusions

In this paper,a new unified method based on the hybrid lattice Boltzmann flux solver is proposed for integrated hypersonic fluid-thermal-structural analysis. In this approach, the LBFS with D1Q4 model and TLBFS with D2Q4 model are used to evaluate the numerical fluxes of the Navier-Stokes equations and the structural heat conduction equation at the cell interface based on FVM respectively. In addition, to predict the heat flux of the hypersonic flows accurately,a modified switch function is introduced to precisely control the numerical dissipation in the present LBFS. Finally, by using a new dualthermal-resistance model for providing the continuity condition of the heat flux across the interface,the LBFS and TLBFS are combined simultaneously into one set of unified method for simulation of integrated hypersonic fluid-thermalstructural problems.

To demonstrate the accuracy, reliability and efficiency of the developed solver, several numerical examples are simulated, including hypersonic flow around a circular cylinder,2D structural-thermal conduction and the CHT on aerodynamically heated cylinder leading edges. Excellent agreements of the numerical results with the experiment and available data in the literature are achieved. These studies indicate that the present hybrid Boltzmann flux solver has great potential for providing insight into the features of the hypersonic fluidthermal-structural interactions and expanding the scope of application of LBFS in the CHT analysis.Furthermore,significant improvements in predicting steady CHT problems for long-endurance hypersonic vehicles can also be achieved in the present integrated method.

Acknowledgments

This study was co-supported by the Postgraduate Research &Practice Innovation Program of Jiangsu Province of China(No. KYCX17_0235), Nanjing University of Aeronautics and Astronautics PhD Short-term Visiting Scholar Project(No. 180602DF01), and National Numerical Wind Tunnel Project (Nos. NNW2018-ZT3B08 and NNW2019-ZT7B30).