Development of a Synchronization Method for Fluid-Thermal Study of Hypersonic Flow

, , ,

College of Aerospace Engineering, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, P.R. China

Abstract: A synchronization method is developed for the fluid-thermal study of hypersonic flow. Different from conventional loosely/tightly coupled methods which separately deal with the flow field and the structure temperature field, the presented method expresses the governing equations in a unified framework so that the two fields can be calculated simultaneously. For efficiently solving the unified equations, the finite volume method together with the dual-time stepping approach is employed. Like in the flow field, the local time step is also used in the temperature field, which is determined from thermal conductivity spectral radii. In order to treat the fluid-structure interface more conveniently, an expanded virtual boundary is introduced. For validation, several fluid-thermal hypersonic flow problems are simulated. The computed results are compared with those obtained from the coupled methods and the experiment. In the continuous heating problems, the stagnation temperatures predicted by both the coupled and synchronization methods are in good agreements with the experimental data. In the unsteady flow-thermal hypersonic flows, the stagnation heat fluxes predicted by the presented method and tightly coupled method are basically the same, which agree better with the experimental data than those predicted by the loosely coupled method. In terms of prediction of the stagnation temperature, the synchronization method shows better accuracy than the tightly coupled method.

Key words: hypersonic flow;aerothermal; synchronization method; expanded virtual boundary; conductivity spectral radii

Nomenclature

0 Introduction

In the research and development of hypersonic vehicles, aerodynamic heating is a major difficulty that needs to be overcome. The transient structural temperature field of a hypersonic vehicle requires accurate computation. In recent years, some coupled methods have been developed to simulate these structure-aerodynamic heating (i.e., flow-thermal) problems. These methods can be classified into two types, namely tightly coupled methods and loosely coupled methods. If a single interface iteration (coupling iteration) is performed at each time step (or stage of an implicit Runge-Kutta scheme), the algorithm is referred to as a loosely coupled method, and otherwise it is referred to as a tightly coupled method[1]. Another way to distinguish them is as follows. Loosely coupled methods presume that the time taken for the flow field to reach equilibrium can be neglected when comparing characteristic time step for heat conduction in the structure; tightly coupled methods presume that the flow field and the heat conduction in the structure should be uniformly solved by a smaller characteristic time step determined by the flow field[2]. Regardless of how these methods are distinguished, most researchers used the exchanging flow field and the structural temperature field boundary conditions to solve the coupling problem when simulating the structural temperature field and the flow field separately[3-9]. Geng et al.[10]proposed an integrated numerical method without iterative coupling of these two fields, which was similar to loosely coupled methods. The flow and heat conduction were described in unified equations. To eliminate huge differences in characteristic time steps between fluid and solid, the flow field was kept “frozen” until the structural surface temperature changed by a certain magnitude.

Most studies on the use of both loosely and tightly coupled methods showed that the former remedies the low computational efficiency of the latter, premised on ensuring the accuracy of the computation[2,7-9]. By using these methods, the information of the two fields needs to be exchanged at the interface, and when one field is being computed, the other is usually kept “frozen”. However, this situation is non-physical because both the flow field and the structural temperature field are changing simultaneously, and thus the solution accuracy may be reduced. In addition, the determination of the coupling time step is critical. A very small time step may increase the computational efforts while a very large value may lead to numerical errors. Duchaine et al.[11]proposed a direct coupling method of the fluid/solid domains in a large system of simultaneous equations that is handled by a “monolithic” solver, called parallel coupling strategy (PCS), and has been tested to deal with a fully non-stationary coupled problem with synchronization in physical time. There are still existing exchanges between fluid solver and solid solver.

To overcome these drawbacks and to satisfy engineering requirements in precision as well, a synchronization method to simulate both flow field and structural temperature field for hypersonic vehicles is developed. An expanded virtual boundary is applied at the interface between fluid and solid components. The local time step required for computations involving the solid component is determined from thermal conductivity spectral radii, which is similar to the viscous spectral radii of the fluid. The whole field then remains synchronized in the time domain. The numerical results of the developed synchronization method are compared with those of both loosely and tightly coupled methods and from experiment.

In what follows, we first describe the governing equations for the loosely/tightly coupled systems, and the numerical method in the gas domain, the solid domain and at the interface(Section2). In Section3, we present the synchronization method based on expanded virtual boundary and the thermal conductivity spectral radii. In the final part of the paper (Section4), the numerical results of the synchronization method are compared with those obtained from the loosely/tightly coupled methods, and experiment data.

1 Loosely/Tightly Coupled Methods

Traditional coupled methods include both tightly and loosely coupled methods, both of which contain a flow field solver and a thermal conduction solver that work independently but exchange information at the boundary. Since response of flow field is much faster than that of heat conduction, the steady-flow Navier-Stokes(NS) equations are used in loosely coupled methods and the time step is determined from the characteristics of the heat conduction, which differs from that for tightly coupled methods.

1.1 Flow field solver

Consider the three-dimensional NS equation described in the Cartesian coordinate system

(1)

whereWrepresents the conservative variable,f,g,qrepresent the convective flux vectors, andR,S,Tthe viscous flux vectors.

In the computation of hypersonic flows, the “advection upstream splitting method” (AUSM+) format has many clear and well-known advantages, such as positive pressure retention, good stability, and high resolution. Viscous fluxes are discretized by the cell central scheme. To improve the accuracy of the calculation, the third-order “monotonic upstream-centered scheme for conservation laws” (MUSCL) interpolation method is applied. As the flow field contains strong gradients, the limiter is required to suppress any numerical oscillation. The MUSCL interpolation method involves

(2)

(3)

whereεis set to 10-6.

For loosely coupled methods, the steady-flow NS equation is implemented in the flow field solver, and the implicit lower-upper symmetric Gauss-Seidel (LU-SGS) scheme is applied to temporal discretization; For tightly coupled methods, the unsteady-flow NS equation is applied to the flow field solver, and the dual time stepping method is applied to temporal discretization. The far field has non-reflecting boundary conditions imposed, and no-slip wall conditions are imposed at the interface between the solid and viscous flow. The conditions imposed for the heat flux at the isothermal boundary and the pressure gradient at the solid surface are

(4)

(5)

1.2 Thermal conduction solver

Assuming no heat sources, the heat equation is given by

(6)

The integral form of the thermal conduction control equation is then

(7)

whereΩrepresents the volume of the control cell. By Stokes’ theorem

(8)

whereΣrepresents the total area of the control cell.

Assuming that thermal properties of the solid (ks,Cs,ρs) are isotropic and homogeneous, Eq.(8) becomes

(9)

wherenis the unit normal vector on control surfacei,Tfacethe gradient of the temperature on each surface and obtained from

(10)

withΩI+1/2,J,Kthe volume of the virtual cell. From Eqs.(8)—(10), we obtain

(11)

which subsequently reduces to the form

(12)

The surface heat flux distribution is acquired by the flow field solver. The additional source term method is applied to define the heat flux as an equivalent source term for the control cell.

The method of determining the local time step in the heat conduction computation is the same as the synchronization method and will be introduced below.

1.3 Coupling procedure

With reference to the diagram (Fig.1), the computation procedure for both loosely and tightly coupled methods is summarized as follows.

(1) The computation is divided into several periods (coupling time steps).To initialize the coupled computations, a steady-state fluid computation is performed att1.

(2)Assuming the flow field is “frozen”, the structural temperature field steps fromt1(initial temperature field) tot2with the “frozen” rate of heating, as well as the interface temperature att2can be acquired using the thermal conduction solver.

(3)Assuming the structural temperature field is then “frozen”, the flow field steps fromt1tot2with the “frozen” interface temperature, and the rate of heating att2is determined using the steady flow field solver.

(4)The final temperature field can be acquired through several iterations between the flow field solver and thermal conduction solver until reaching the final time point (tn).

Fig.1 Procedure for traditional coupling methods

2 Synchronization Method

The traditional coupled methods need to exchange the information of flow field and structural temperature field at the interface, and when one field is being computed, the other is usually kept “frozen”. Moreover, a suitable coupling time step is required. Thus, a synchronization method is proposed to overcome these drawbacks and also to satisfy engineering requirements in precision.

2.1 Governing equations and discretization

For the fluid component, the discretization method is the same with the flow field solver in the traditional coupling method. For an ideal gas

(13)

For the solid component,u=v=w=0. Hence, the first three scalar equations in Eq.(1) are naturally satisfied. The fourth equation in Eq.(1) determines the energy

ρsEs=ρsCsTs

(14)

The discretization method for the heat equations is the same with Eqs.(6)—(12). In the following section, the determination of local time steps and treatments of the boundary will be presented.

2.2 Determination of the local time steps

Jameson’s dual time stepping method is used for temporal discretization, and the implicit LU-SGS scheme is applied to time stepping in the pseudo-iteration. Local time steps for the fluid can be determined using

(15)

(16)

(17)

wherePLandPTrepresent laminar and turbulent Prandtl numbers, respectively;crepresents the local speed of sound. The fluid heat conductivity is given by

(18)

or equivalently

(19)

(20)

According to the definition ofthe spectral radii of viscous flux for the fluid, and from Eq.(19), the spectral radii of thermal conduction is acquired from

(21)

The method of introducing the thermal conduction spectral radii to determine the local time steps is equivalent to the definition of “viscous” in the heat conduction. This ensures that the fluid thermal conduction is consistent with the solid thermal conduction during each pseudo-iteration process. Additionally, there is no need to determine the coupling time step or to freeze the field for either flow or structural temperature.

2.3 Initial and boundary conditions

The static initial flow field and temperature field att=0 s are required before the synchronous computation. The boundary in the far field is treated using the same condition in loosely/tightly coupled methods. In the following section, the implementation of the boundary conditions between fluid and solid will be presented. The initial coefficient of viscosity of the interface can be determined by the surface temperature and the Sutherland formula, and then the thermal conductivity coefficientkwof the initial surface can be obtained.

Because of differences in physical properties (such as thermal conductivity and coefficient of viscosity) between fluid and solid, an expanded virtual boundary is introduced to compute the viscous flux of the interface. “Expanded” means that the virtual cells and real cells have same physical properties. Specifically, as sketched in Fig.2, the cell markedK=0 represents the expanded virtual cell, which has the same properties as the cellK=1. These cells are not real interface cells, but presented for the computation of the boundary derivatives and called expanded virtual cells. These are used to compute the viscous flux, viscous coefficient, and thermal conductivity. Nevertheless, the temperature of interface should be determined by the real cells of the fluid and solid. The temperature of the expanded virtual cell for the fluid is obtained by interpolation

T0=2×Tw-T1

(22)

whereT0,T1represent the temperature of cellsK=0,1, respectively, andTwdenotes the temperature of the interface. It can be obtained by pseudo-time stepping using the continuity of heat flux at the interface[12]

(23)

wheredfanddsdenote the distances between the center of the first layer and the interface in the fluid and solid, respectively;TfandTsthe temperature in the first layer of the fluid and solid;ksandkfthe heat conductivities of solid and fluid parts, respectively, determined by the interface temperature. The pressure gradient normal to the wall and the velocity at the wall are zero. The interface density can be obtained using the ideal gas equation.

Fig.2 Expanded virtual boundary

At the interface, the coefficients for heat conductivity and viscous are considered as functions of temperature. Consequently, the implementation of the boundary treatments combines the computation of the flow field and the temperature field more tightly.

3 Applications

Loosely/tightly coupled methods and the synchronization method are assessed by simulating Mach-6.47 flow over a two-dimensional cylinder and Mach-9.86 flow over a three-dimensional biconical surface. The resultant flow behaviors are calibrated using experimental data. All the computations are performed on multiblock structured grids. All lengths are in meters unless otherwise stated.

3.1 Mach-6.47 flow over a two-dimensional cylinder

Similar to Ref.[13], the test model consists of a 12.7-mm-thick, 76.2-mm-diameter, and 321 stainless steel cylinder. Details of the experimental configurations, the tunnel flow conditions, and the experimental results were also given in Ref.[13]. The setting for the thermal properties of the solid and the aerodynamic parameters applied in this section are

Ma=6.47,α=0°,T=241.5 K,Tw=294.4 K

ρ=1.012×10-2kg/m3,p=7.018×102Pa

ρs=9030 kg/m3,ks=16 W/(m·K)

Cs=500 J/(kg·k)

The grids representing the flow domain and the cylinder (Fig.3) consist of 12 000 cells in the flow domain and 3 000 cells in the cylinder with the same discretization along the interface. Following Ref.[14], the smallest cell in the flow domain is only 0.010 16 mm in the direction normal to the cylinder surface. A detailed study of this model was presented in Ref.[14], and it has been used to validate algorithms by many researchers. In this paper, the numerical results are also compared with those in Ref.[14].

Fig.3 Fluid-solid finite volume model for flow over a 76.2-mm-diameter cylinder

The major difficulty in the fluid analysis is the prediction of the rates of aerodynamic heating, because a very accurate resolution of the flow temperature gradient normal to the cylinder surface is required. Both the synchronization method and loosely/tightly coupled methods need a static initial flow field att=0 s before computations. The computed flow temperature distribution att=0 s (Fig.4) along the flow symmetry line (y=0) clearly illustrates the sharp gradients that must be resolved. The free-stream temperature increases sharply from 241.5 K to about 2 200 K across the bow shock. Within a very thin layer at the flow stagnation point, the temperature drops sharply from 2 200 K to a surface temperature of 294.4 K. Thus, these computation results are in good agreement with Ref.[14]. Comparison of the density contours with those of Ref.[14] (Fig.5) shows that the shock wave is well captured.

Fig.4 Comparison of fluid temperature distributions between the reference and present model taken along the centerline of the initial flow

Fig.5 Comparison of density contours between reference and present model for Mach-6.47 flow over a 3-in-diameter cylinder

And then, we can continue the computation of coupled methods and synchronization method. For loosely coupled methods, the flow field and the structural temperature field are updated every 0.2 s. Physical time step of tightly coupled methods should be consistent with the synchronization method, determined by a trial computation of stagnation temperature att=2.5e-3 s with different physical time steps. The results are obtained by synchronization method with 100 pseudo-iterations (Table 1).

Table 1 Stagnation temperature at t=2.5e-3 s with

Fig.6 Comparisons of predicted temperature contours at t=2 s of the loosely coupled method, the tightly coupled method and the synchronization method with the result in Ref.[14]

Comparisons of the temperature contours att=2 s among the reference, the loosely/tightly coupled methods, and the synchronization method are shown in Fig.6, in which temperature values are cell centered. Since the dimension of temperature in Ref.[14] is given using the Rankine scale, the temperature has been converted to enable a clear comparison of contours. For the loosely coupled method, the stagnation temperature is 386.339 K (695.410 R) att=2 s; For the tightly coupled method, the stagnation temperature is 386.049 K (694.888 R) att=2 s; For the synchronization method, the stagnation temperature is 386.043 K (694.877 R) att=2 s; In Ref.[14], the stagnation temperature is 388.889 K (about 700 R) att=2 s. Obviously, the results of three methods are in good agreement with those from the experiment.

Comparisons of time histories of the stagnation heat flux between the experiment, the loosely/tightly coupled methods, and the synchronization method [Fig.7(a)] show that the stagnation heat flux drops sharply during the early period of the computation, but the trend gradually slows down in the time remaining. Due to the influence of heat flux, the time histories of stagnation temperature of these methods increase in the similar trend [Fig.7(b)].

In this case, the model is continuously heated by hypersonic flow, the stagnation temperatures with these three methods are all found to be in good agreement with experiment. The temperature contours, stagnation temperatures and heat fluxes of the three methods are almost same. Therefore, it is acceptable to apply the either the loosely/tightly coupled methods or the synchronization method to aerothermal analysis of continuously heated problem in engineering.

Fig.7 Comparisons of heat-transfer rates and temperatures at stagnation between loosely coupled method, tightly coupled method, and synchronization method

3.2 Mach-9.86 flow over a three-dimensional bicone

The model (Fig.8) used in the experiment was fabricated from stainless steel except for the nose tip, which was made from MACOR machinable glass-ceramic. For the convenience, the model presented here is uniquely machined from the MACOR. In accounting for the variation with temperature, the thermal properties of MACOR, taken from Ref.[15], are as follows.

For 300 K

ks=1.614 W/m·K

For 295 K

Cs=193.53+3.348 13×T-1.001 61×10-2×

T2+2.811 25×10-5×T3-4.248 05×

10-8×T4+2.334 33×10-11×T5J/(kg·K)

For 298 K

ρs=2 543.84-8.0×10-2Tkg/m3

The grids representing the flow domain and the bicone (Fig.9) consist of 1 436 974 cells in the flow domain and 1 048 174 cells in the bicone with the same discretization along the interface. The height of the first layer is determined from the grid Reynolds number[16]

(24)

where the feature sizenis usually taken as the normal height of the first layer of mesh near the wall. In accordance with Ref.[16], the smallest cell in the flow domain is 0.107 mm in the direction normal to the surface of stagnation. Details of the experimental configurations, the tunnel flow conditions, and the experimental results are given in Ref.[17]. The aerodynamic parameters applied in this section are

Ma=9.86,α=0°,T=48.88 K,Tw=300 K

ρ=4.271×10-3kg/m3,p=5.992×10 Pa

Fig.8 Details of the profile of the model (mm)

Fig.9 Symmetry plane for flow domain and bicone

As mentioned in Ref.[17], in the experiment,t=0 s is defined as the initiating time when the surface temperature of the model begins to increase as it is inserted into the flow from a sheltered position. According to the time history of heat flux, the time interval can be divided into two stages. During the first 0.35 s, the flow field is not stable, and the model can be considered as be accelerating in this sub-interval; the heat flux reaches a maximum att=0.35 s and begins decreasing slowly after the peak point. Assuming that the density and pressure of the free-stream remain constant, the velocity is considered to increase with time to simulate the accelerating period. Kang et al.[18]applied the Kemp-Riddell corrected formula

(25)

From the time history of the heat flux for the reference, and assuming that the wall temperature is a constant (300 K) during the accelerating period, the Mach number during the intervalt=0 s tot=0.35 s can be obtained from Eq.(25), which is displayed Table 2.

Table 2 Variation of Mach number from t=0 s to t=0.35 s

For the loosely coupled method, the flow field and the structural temperature field are updated every 0.1 s. For both the tightly coupled method and the synchronization method, the flow field and the structural temperature field are updated every 0.001 s with 400 steps pseudo-iteration. Actually, the experimental environment is more complex than that assumed in the numerical simulation because keeping the density and pressure of the free-stream constant is difficult. The difference between the numerical methods and experiment is attributed to free-stream turbulence in the test stream, which is not taken into account by the numerical method. Moreover, the whole computation concludes subsonic flow, transonic flow, supersonic flow, and hypersonic flow. The height of the first layer and the far field of the grid are designed for hypersonic flow, resulting in computation errors during the subsonic and transonic periods. Furthermore, differences between the assumption (a constant wall temperature of 300 K during the accelerating period) and the experiment may influence the evaluation of the Mach number for the free stream using Eq.(26). Additionally, apart from the MACOR nose tip, the model is fabricated from stainless steel. Obviously, as noted in Ref.[17], the entire model ideally should be made of MACOR glass-ceramic, but for this experiment it is not structurally or economically feasible.

Comparisons of the longitudinal heat transfer rate and temperature distribution att=1.0 s among the experiment, the loosely/tightly coupled methods, and the synchronization method (Fig.10) show that the heat fluxes and temperatures obtained by all these numerical methods are basically the same in the non-stagnation region. The synchronization method performs better than the other methods in calculation of stagnation heat flux and temperature.

Fig.10 Comparisons of the longitudinal heat transfer rates and temperature distributions at t=1.0 s between the experiment, the loosely coupled method, the tightly coupled method, and the synchronization method

Comparisons of temperature contours in solid domain att=0.35 s between the loosely coupled method, the tightly coupled method, and the synchronization method (Fig.11) show that the temperatures obtained by the synchronization method and the tightly coupled method, which are basically the same, are slight higher than those obtained by the loosely coupled method. That is why the surface heat fluxes obtained by the loosely coupled method are slight higher than those obtained by the other two methods att=0.35 s (Fig.12). Because of the increment of temperature, the surface heat flux att=1.5 s reduces evidently compared with that att=0.35 s as shown in Fig.12.

Fig.11 Comparisons of temperature contours in solid domain at t=0.35 s and t=1.5 s between the loosely coupled method, the tightly coupled method, and the synchronization method

Fig.12 Comparisons of surface heat flux contours at t=0.35 s and t=1.5 s between the loosely coupled method, the tightly coupled method, and the synchronization method

Time-history comparisons of the heat flux at stagnation among the experiment, the loosely/tightly coupled methods, and the synchronization method (Fig.13(a)) show that the heat fluxes obtained by the synchronization method and the tightly coupled method, which are basically the same, are closer to the experimental data than those obtained by the loosely coupled method. As reported in Ref.[17], the measured results for the temperature at stagnation are 402.3 K, 461.1 K, and 493.8 K att=0.5 s,t=1.0 s, andt=1.5 s, respectively. Similar comparisons for the temperature at stagnation (Fig.13(b)) show that the final stagnation temperatures att=1.5 s, computed by the loosely/tightly coupled methods, and the synchronization method are 482.08 K, 486.30 K, and 490.58 K, respectively. The stagnation temperature increases sharply from 300 K to about 500 K in 1.5 s. Similarly, the whole aerodynamic heating is composed of two stages: (1) The free-stream varies duringt=0 s tot=0.35 s, and the heat flux grows sharply to a maximum att=0.35 s (unsteady flow); (2) fromt=0.35 s tot=1.5 s, the rate of aerodynamic heating decreases slowly (continuous heating flow).

During the first stage,the Mach number increases from 0.37 to 9.86, and the stagnation temperature grows markedly by about 50 K in 0.35 s. Indeed, the stagnation temperature increases slightly during the first 0.2 s. In other words, the rise in stagnation temperature mainly occurs within 0.15 s (fromt=0.2 s tot=0.35 s). The unsteady flow effect of this period is very significant. Hence, the differences between the synchronization method and the loosely coupled method are mainly manifested during this period. The results (Fig.13(a)) during the accelerating period computed by the tightly coupled method and the synchronization method agree better than those by the loosely coupled method and experiment data in terms of accuracy. This is mainly caused by the unsteady flow field during this time period.

During the second stage, the model is continuously heated by hypersonic flow, similar to the two dimensional case mentioned above. At the beginning of the second stage, there still exists unsteady effect in flow field, resulting in the difference of stagnation temperature between the tightly coupled method and the synchronization method. The stagnation temperature of the synchronization method matches the experiment data better than that of the tightly coupled method during this period in terms of the accuracy. In the following computation, the trends of heat flux and temperature of these three methods are almost same.

Fig.13 Time history comparisons of heat transfer rates and temperatures at stagnation between the experiment, the loosely coupled method, the tightly coupled method, and the synchronization method

4 Conclusions

In this paper, an accurate and efficient synchronization method is developed for fluid-thermal study of hypersonic flows, through which three major conclusions are drawn as follows.

(1) In the developed method, the flow field and the structure temperature field are calculated simultaneously by solving the unified governing equations. As compared with conventional coupled methods, the synchronization method is physically more reasonable since fewer approximations are made. Numerically, without the need to combine different methods for solving the flow field and the structure temperature field, it is easier to be implemented.

(2) In the continuous heating problems, the stagnation temperatures and heat fluxes predicted by all the three methods are in good agreements with the results of the reference. It is acceptable to apply either loosely/tightly coupled methods or the synchronization method to aerothermal analysis of the continuous heating problems in engineering.

(3) In the unsteady flow-thermal hypersonic flows, the stagnation heat fluxes obtained from the synchronization method and the tightly coupled method are closer to the experimental data than those obtained from the loosely coupled method. Further, the stagnation temperature obtained from the synchronization method is closer to the experimental data than that obtained from the tightly coupled method. Thus, the presented synchronization method is a more accurate method for the unsteady hypersonic flow-thermal computations.

Acknowledgements

This work was supported by the National Natural Science Foundation of China(No.11872212). All numerical simulations have been conducted on the computers of Aerodynamic Department of College of Aerospace Engineering at Nanjing University of Aeronautics and Astronautics. The authors also gratefully acknowledge the technical insights of associate professor Wu Jie of Aerodynamic Department of College of Aerospace Engineering at Nanjing University of Aeronautics and Astronautics.