A novel TE-material based thermal protection structure and its performance evaluation for hypersonic flight vehicles

Aerospace Science and Technology 77 (2018) 458–470

Contents lists available at ScienceDirect

Aerospace Science and Technology www.elsevier.com/locate/aescte

A novel TE-material based thermal protection structure and its performance evaluation for hypersonic flight vehicles Chun-Lin Gong a,∗ , Jian-Jun Gou a , Jia-Xin Hu a , Feng Gao b a

Shaanxi Aerospace Flight Vehicle Key Laboratory, School of Astronautics, Northwestern Polytechnical University, Xi’an 710072, China State Key Laboratory of Solidification Processing, MIIT Key laboratory of Radiation Detection Materials and Devices, College of Material Science and Engineering, Northwestern Polytechnical University, Xi’an 710072, China b

a r t i c l e

i n f o

Article history: Received 21 November 2017 Received in revised form 15 January 2018 Accepted 17 March 2018 Available online 21 March 2018 Keywords: Thermal protection system Thermoelectric materials Hypersonic flight vehicle Numerical simulation

a b s t r a c t The traditional thermal protection system (TPS) of hypersonic flight vehicles is usually designed for thermal protection purpose only with low efficiency. In this paper, a thermoelectric material based multifunctional TPS structure concept is proposed and the evaluation approach of its mechanicalthermoelectric performance is developed based on a specific vehicle and a typical trajectory. The thermoelectric module in the structure can convert a certain amount of aerodynamic heat into electricity supply. The module consists of n-type Sr0.9 La0.1 TiO3 compound which is fabricated based on the solid state reaction method, and the widely used p-type Ca3 Co4 O9 . For a specific hypersonic flight vehicle with a typical trajectory curve, the aerodynamic heat is calculated by an engineering-based algorithm, the unsteady mechanical-thermoelectric characteristics of the structure is then analyzed based on a unit cell model and the thermoelectric conversion efficiency is finally evaluated. The results indicate that the multifunctional TPS structure has significant application potentials on the hypersonic flight vehicles. © 2018 Elsevier Masson SAS. All rights reserved.

1. Introduction For a flight vehicle, the increasing flying speed brings about a large amount of aerodynamic heat and leads to a very harsh thermal environment. For a supersonic vehicle, at some specific locations of leading edge, the temperature may reach to a value of higher than 1800 K [1–3]. Thus, a reliable and efficient thermal protection system (TPS) is essential to ensure the safety of the vehicle. The first input condition of design and evaluation of a TPS is the aerodynamic heat. There are two main methods to predict the aerodynamic heat, the engineering-based method and the numerical method. During the last decades, numerous relevant studies have been carried on to develop appropriate numerical method (mainly CFD) for the aerodynamic heat prediction of flow fields/vehicle with different configurations [4–8]. For instance, in Ref. [4], the capacity of DNS (direct numerical simulations), LES (large eddy simulations) and RANS (Reynolds-averaged Navier–Stokes) for the prediction of turbulent shock wave were evaluated based on cases of 2D compression and expansion-compression corners, 2D shock impingement, 3D single and double fins. At present, the numerical method still has some challenges, such as the high requirement

*

Corresponding author. E-mail address: [email protected] (C.-L. Gong).

https://doi.org/10.1016/j.ast.2018.03.028 1270-9638/© 2018 Elsevier Masson SAS. All rights reserved.

to mesh generation, especially the grid near the vehicle surface should be very carefully considered since the temperature gradient is very sensitive to the grids resolution. Also, the numerical simulation is a very resource-consuming method and has low efficiency. On the other hand, the engineering-based algorithms and programs, such as MINIVER [9] and LATCH [10] which are developed based on empirical correlations and engineering purposes are very popular. Such methods always contain a series of theories to adapt different flow conditions (e.g., subsonic or supersonic, laminar or turbulent). Thus, they have lower accuracy in comparison with numerical methods. However, they are much more efficient. In this work, a quick TPS evaluation approach is more concerned, and thus an engineering-based method is used to calculate the aerodynamic heat. A TPS always consists of a high-temperature layer and a below insulation layer. The high-temperature layer is the first layer to withstand the aerodynamic heating, and the candidate materials include C/SiC, C/C, ultra high temperature ceramic, etc. For the insulation layer, the basic requirement should be the low weight and low thermal conductivity [11,12]. Ma et al. [11] developed an integrated TPS with additional insulation layer, and different insulation material/structure including Mullite blanket, ceramic tile and aerogel are evaluated. In fact, the choosing of a TPS material is closely related to the material’s own characteristics and the TPS requirement. For instance, the TPS optional materials can be classified into

C.-L. Gong et al. / Aerospace Science and Technology 77 (2018) 458–470

459

Fig. 1. The vehicle and the trajectory.

two types in general, the ablative and non-ablative materials. For a structure with ablative materials [11,13–16], part of the aerodynamic heat will be taken away by its ablation process, and thus has higher efficiency. However, it is not applicable to a reusable vehicle due to its changing configuration during the flight mission. On the other hand, for a TPS based on non-ablative material [17–19], the surface radiation is the only dissipation of aerodynamic heat, and thus leaves some rooms for the improvement of its efficiency. In this work, a high-efficiency multi-functional TPS structure based on thermoelectric (TE) material is developed. The TE material and corresponding modules can convert aerodynamic heat into electric energy directly as the well-known Seebeck effect, and the conversion efficiency can reach a value of 5∼20% or even higher with some appropriate treatments like doping and nano-structuring [20]. According to its suitable temperature, the TE material can be classified as room-temperature (e.g., Bi2 Te3 ), mid-temperature (e.g., PbTe) and high-temperature ones (Oxide ones). A TE module is always a couple of n- and p-type materials and can be used as an electricity generator under a certain temperature difference. The performance of a TE generator is partly dependent on the temperature difference across the TE material which is related to the performance of additional heat exchangers, Lv et al. [21] studied three typical heat exchangers (air cooling, fluid cooling and heat pipe) and evaluated their influence on TE performance. Lan et al. [22] developed a dynamic model of TE generator for the waste heat recovery of automotive vehicles. For the application of flight vehicles, to the authors’ knowledge, the relevant researches are limited: Li and Wang [23] developed an integrated TE module with regenerative cooling system, and analyzed its performance based on the exergy analysis theory; Cheng et al. [24] developed a multi-stage TE module considering the suitable temperature of TE material, and analyzed its thermal protection performance. The very limited research activities motivate further studies. In this work, a TE material based multi-functional TPS structure is developed. For the p-type material, the widely used hightemperature material Ca3 Co4 O9 [25–27] is adopted, and for the n-type material, a Ti-doped ceramic Sr0.9 La0.1 TiO3 is fabricated by solid state reaction method. The evaluation approach of its mechanical-thermoelectric performance is developed based on a specific vehicle and a typical trajectory.

2. The vehicle, trajectory and aerodynamic heat It should be pointed out that under a limiting or an ideal condition, engineering structures can always have good performance, whereas relevant results have certain limitations. Therefore, the evaluation of an engineering system has to be accomplished according to a specific vehicle with specific mission trajectory. For the novel TPS structure proposed in this work, a complete evaluation process including the calculation of input thermal conditions, the mechanical-thermoelectric analysis and the TE conversion calculation is displayed in detail. Such approach can be referenced and used for other similar vehicles and trajectories. 2.1. The vehicle and trajectory The vehicle and the trajectory studied in this work should be introduced first. Fig. 1 shows the vehicle and the trajectory. The left part is the model of a hypersonic reusable launch vehicle designed by the authors. The black circle on the compression surface is a specific point M that will be considered in the calculation of this work. The size of the vehicle and the general distribution of structure materials are shown in Fig. 1. The vehicle has a length of 30 m, a wingspan of 15 m, and the width of the fuselage is 5 m. The compression surface has an area of about 30 m2 . The structure of the vehicle is divided into several zones of different materials mainly consists of Ti-alloy, Al-alloy and resin based composites, and more detailed descriptions can be found for a similar vehicle in the authors’ previous work [28]. The right part of Fig. 1 is its typical trajectory curve in a twodimensional pattern. The vertical axis is the altitude, the bottom horizontal axis is the flight time while the top one is the flight range. It should be noted that the flight range is not in a linear scale, and the largest value is about 6600 km and then the vehicle will return to the launch field. One can notice that the vehicle will reach 28 km height and its speed will increase from 0 to 8 Ma within 450 s, then it will fly to its highest altitude of about 50 km and accomplish its mission and return. The large amount of aerodynamic heat generated during the trans-atmospheric phase of its launch and re-entry gives rise of great challenges of TPS. In this work, the aerodynamic heat under the trajectory curve shown in Fig. 1 is calculated. The heat is then applied to analyze the thermal protection efficiency of the TPS structure proposed in this work.

460

C.-L. Gong et al. / Aerospace Science and Technology 77 (2018) 458–470

Table 1 Aerodynamic prediction methods for different regions. Configurations

Head Fuselage Wing

Flat Flat flat

2 < Ma < 6

Ma > 6

Windward

Leeward

Windward

Leeward

Tangent–Wedge Tangent–Cone Tangent–Cone

Dahlem–Buck Dahlem–Buck Dahlem–Buck

Modified Newtonian Modified Newtonian Modified Newtonian

Prandtl–Meyer Expansion Prandtl–Meyer Expansion Prandtl–Meyer Expansion

2.2. Aerodynamic heat In order to evaluate a TPS, the input condition which is the aerodynamic heat should be obtained first. For a non-ablative TPS system, the heat balance equation can be described by Eq. (1): 4 q = hconv ( T r − T w ) = σ ε T w + qcond

(1)

where q is the aerodynamic heat flux, hconv is the coefficient of convection, T r is the recovery wall temperature, T w is the wall temperature, σ = 5.67 × 10−8 W/(m2 K4 ) is the Stefan–Boltzmann constant, ε is the surface emissivity, and qcond is the heat flux conducted into the vehicle. The surface material of TPS is C/C-SiC and the emissivity is 0.8 [29] in this work. In this work, the aerodynamic heat is calculated by a fast engineering algorithm which is based on Prandtl’s boundary layer theory. The flow field is divided into inviscid flow which is external to the boundary layer and viscid field which is internal of the layer. The algorithm involves several theories of aerodynamic force and aerodynamic heat prediction. And three basic assumptions have to be stated first: a. In the calculation of aerodynamic force of the inviscid flow field, the head, the fuselage, and the wings can be considered as flat, and the corresponding methods are listed in Table 1. b. For the viscid compressible flow field internal of the boundary layer, it can be simulated by the incompressible theory with a compressibility correction, i.e., the reference temperature (enthalpy) method. c. For the situation of Ma < 2, the present method has low accuracy, thus the aerodynamic heat in that duration (0∼140 s) can be calculated by linear interpolation based on available data. For the inviscid flow field, the gas flow is considered as ideal and the aerodynamic parameters are calculated based on the isentropic conditions. The corresponding methods are listed in Table 1, and the detailed selection rationale (mainly based on speed and configuration) can be found in Ref. [30]. The calculation of inviscid field provides parameters of the outer edge of the boundary layer, which will be the input data in the viscid flow calculation. For the viscid flow internal of the boundary layer, the thermodynamic parameters are computed based on a reference temperature method, in which the reference temperature T ∗ is defined by:

T ∗ = T e + 0.58( T w − T e ) + 0.19( T r − T e )

( Mae < 5) T = 0.7T e + 0.58( T w − T e ) + 0.23( T r − T e ) (5 < Mae < 10) ∗

(2) where the subscript “e” indicates the outer edge of the boundary layer; T r is the recovery temperature which can be calculated based on the parameters (provided by the calculation of inviscid fields) of the outer edge of the boundary layer. The convective coefficient hconv in Eq. (1) can then be obtained based on the thermodynamic parameters. Based on iterative solving, the aerodynamic heat flux can be obtained. In the calculation, the radiation equilibrium condition

Fig. 2. The aerodynamic heat flux.

is assumed which means qcond in Eq. (1) is considered to be 0. Fig. 2 shows the heat flux at the location of specific point M. The horizontal axis is time and the left vertical axis is heat flux. It is well-known that the aerodynamic heat flux is closely related to the flight altitude and speed. In general, the lower altitude (larger air pressure and density) with higher speed lead to larger heat flux. From Fig. 2, one can find that for the vehicle and trajectory studied in this work, the aerodynamic heat flux increases to its highest value in a duration of about 450 s (trans-atmosphere phase), and then decreases to a valley value in a very short time (rarefied gas flow), and then fluctuates in a relatively lower value range. The heat flux in Fig. 2 will be used as input thermal boundary conditions in the later mechanical-thermoelectric coupling calculation. 3. The multi-functional TPS structure Fig. 3 is the schematic of multi-functional TPS structure proposed in this work. A three-layer structure is designed. The top layer is the high-temperature C/C-SiC satin woven composite. The middle layer consists two phases, the insulation material of Saffil Al-fiber and the TE materials. The bottom layer is Ti-6Al-4V alloy. The TE module is always composed of a p-type and a n-type composite which in this work are Ca3 Co4 O9 and Sr0.9 La0.1 TiO3 based materials, respectively. In a temperature gradient field, the n-type material will create an electric potential which is positive in the high temperature side while negative in the lower temperature side. This phenomenon is known specifically as the Seebeck effect. However, the p-type material has the exact opposite Seebeck effect, and thus a relative larger electric potential will be created by the couple of p and n materials. Both the C/C-SiC and the Ti-6Al-4V layer have a thickness of 5 mm, and each p and n material has a dimension of 3 × 3 × 15 mm. It should be noted that the transverse size of the structure can be an arbitrary value depends on its application. Fig. 4

C.-L. Gong et al. / Aerospace Science and Technology 77 (2018) 458–470

461

Fig. 3. The multi-functional TPS structure.

by Terasaki et al. in 1997 [31], metal oxide-based materials have been attracting continuous attentions. However, Na is volatile and the compound is instable in humid environment. On the other hand, another Co-based oxide Ca3 Co4 O9+δ has good TE performance ( Z T = 0.83 at 973 K for the single crystal) [32] with high thermal and chemical stabilities even up to 1200 K in air [25,33, 34]. In addition, the addition of Ag+ and Lu3+ can improve the Z T of the material (Ca2.85 Ag0.05 Lu0.15 Co4 O9+δ exhibits a Z T of up to 0.61 at 1118 K [35]). Therefore, the p material used in this work is a Ag+ and Lu3+ doped Ca3 Co4 O9 based compound. 3.2. Fabrication of n-type material and its TE properties

Fig. 4. The heat transfer path.

is the transfer path of aerodynamic heat from the vehicle top surface to the internal structure. If the aerodynamic heat generated in the structure surface is expressed as Q , before it enters the structure, a certain amount of heat Q rad will be dissipated by the radiation. When it transmits through the high-temperature layer and comes to the insulation layer, an amount of Q 1 will be converted into electric power by TE materials. And when the rest heat transferred through the Ti-6Al-4V, another certain amount Q 2 will be taken away by the convection heat transfer (possible bottom heat exchanger). Finally, the heat conducted into the vehicle internal structure is Q 3 . It is clear that Q = Q rad + Q 1 + Q 2 + Q 3 . The design of such a TPS structure has two key points, i.e., the choosing of p-type TE material, and the fabrication of n-type material and the measurement of its TE properties such as Seebeck coefficients, which will be discussed in the later subsections. 3.1. The choosing of p-type material By converting heat directly into electricity, TE generator offers a promising technology. The key to realize an efficient TE generator lies in finding materials with high TE performance: good durability and high robustness at high temperature in air. For a TE material, its conversion efficiency is often evaluated by the TE figure of merit ( Z T ):

ZT =

kS 2

λ

T,

S =−

E T

(3)

where S is the Seebeck coefficient, k is the electrical conductivity, λ is the thermal conductivity of the material,  E is the thermoelectric voltage and  T is the temperature difference. By far the most widely used TE materials are alloys of Bi2 Te3 , PbTe and SiGe, which are often harmful, expensive or have poor durability at high temperature. Metal oxides have been considered as an alternative to overcome these disadvantages. Since the discovery of large TE power in p-type NaCo2 O4 single crystals

A n-type TE material Sr0.9 La0.1 TiO3 is fabricated based on the solid-state reaction method in this work. Reagent-grade strontium carbonate (SrCO3 ), titanium oxide (TiO2 ), lanthanum oxide (La2 O3 ) and bismuth oxide (Bi2 O3 ) were used as starting materials. Firstly, the mixtures of SrCO3 , TiO2 , La2 O3 and Bi2 O3 were weighed and ball milled by zirconia balls (as grinding media) in ethanol for 12 h. The mixtures were then calcined at 1200 ◦ C and ball milled again in ethanol for 12 h. After drying, the powders were granulated by PVA (as a binder) and then pressed into pellets with the radius of 30 mm and thickness of 3 mm. The compacts were heated to 500 ◦ C for 2 h to burn out the binder and then sintered at 1450 ◦ C for 2 h in argon atmosphere. Then samples were cut into the size of 3 × 3 × 15 mm. In order to further reduce the electrical resistivity, an annealing process was taken at 1350 ◦ C for 8 h in argon and graphite atmosphere. The detailed fabrication process can be found in Fig. 5(a) and can also reference the authors’ previous work [36]. The Seebeck coefficient is the key property of TE materials. In this work, the Seebeck coefficient of the fabricated n-type material was measured by a 4-terminal method. Fig. 5(b) is the schematic of the measurement. A sample is set in a vertical position between the upper and lower blocks in the heating furnace, and there is an additional secondary heater in the lower block to provide a temperature gradient. The temperature difference  T and the thermal electromotive force  E can be measured by the thermocouples pressed against the side of the sample, and then the Seebeck coefficient can be obtained based on Eq. (3). Also, if a constant current I is applied to both ends of the sample, the voltage drop  V can be obtained, and the electric resistance can then be calculated. Fig. 6 shows the Seebeck coefficients and the electric resistivities of p- and n-type materials used in this work, in which that of p-type comes from Refs. [37,38]. Other physical properties are listed in Table 2.

462

C.-L. Gong et al. / Aerospace Science and Technology 77 (2018) 458–470

Fig. 5. The fabrication of n-type material and the measurement of Seebeck coefficient.

3.3. Physical properties of other materials The physical properties of C/C-SiC and Ti-6Al-4V alloy can be found in Ref. [39]. The elastic modulus and the Possion’s ratio of TE materials are 105 GPa and 0.25, respectively, and that for Saffil material are 20 MPa and 0.3, respectively. The thermal expansion coefficient of Saffil Al-fiber is assumed to be 100 × 10−6 K−1 in this work. 4. Numerical model 4.1. Governing and constitutive equations

Fig. 6. The Seebeck coefficient and electric resistivity.

Table 2 Physical properties of TE materials.

300 400 500 600 700 800 900 1000 1100

Thermal conductivity (p) (W/(m K))

Thermal conductivity (n) (W/(m K))

Thermal expansion coefficient (p–n) (1/K)

1.85 1.68 1.55 1.51 1.48 1.44 1.42 1.41 1.41

2.417 2.409 2.124 1.929 1.768 1.764 1.732 1.511 1.292

0.16 1.21 10.2 10.5 10.7 11.0 11.4 11.6 11.7

In this work, the mechanical-thermoelectric performance of the TPS structure is coupling analyzed by ANSYS Multiphysics which is based on Finite Element Method. In ANSYS, there are two types of coupled-field simulations, the load transfer and the direct coupling analyses. For the load transfer method, the fields are coupled by applying results from one field analysis as input loads in another field analysis. Such method is time saving while has lower accuracy. Therefore, in this work, the direct method is adopted. In order to save computational time, the thermoelastic and thermoelectric coupling are analyzed separately, and such treatment is reasonable, since to some extent the mechanical deformation and the electrical potential are independent to each other. The governing equations for the problem studied in this work can be listed as Eqs. (4), (5) and (6), respectively: Mechanical:

Equilibrium:

∂ τxz ∂ σxx ∂ τxy + + = 0; ∂x ∂y ∂z ∂ τ yx ∂ σ y y ∂ τ yz + + = 0; ∂x ∂y ∂z ∂ σzz ∂ τzx ∂ τzy + + =0 ∂x ∂y ∂z

C.-L. Gong et al. / Aerospace Science and Technology 77 (2018) 458–470

463

Fig. 7. The model, mesh and boundaries.

Strain-displacement: ∂u ∂v

εx =

; εy = ; ∂x ∂y ∂w ∂v ∂w εz = ; γ yz = + ; (4) ∂z ∂z ∂y ∂u ∂ w γxz = + ; ∂z ∂x ∂u ∂ v γxy = + ∂y ∂x where σxx , τxy . . . are stress components; εx , ε y , γ yz . . . are strain components; u, v and w are displacements in x, y and z directions, respectively. Electrical:

[k]∇ 2 φ = 0

(5)

The unit cell model (model formulation, model structure, mesh and boundaries) is shown in Fig. 7. As shown in the left picture of Fig. 7, the TPS structure is a typical translational symmetric structure along x- and y-axis, and thus a unit cell model which consists of only one p–n material couple can be formulated to represent the macro structure. It is clear that if the unit cell is translational moved along x- and y-axis, it will cover the whole structure. The unit cell model has a dimension of 12 × 6 × 25 mm. In comparison with the unit cell, the dimension of the TPS structure can be considered as macroscopic. The model is meshed by the couple-field 3D solid element SOLID226 and has about 18046 elements and 79825 nodes. SOLID226 has 20 nodes and each node can have different (up to 5) degrees of freedom: temperature and three displacement components for thermoelastic coupling, or temperature and electric potential for thermoelectric coupling.

where [k] is the electric conductivity matrix and is isotropic in this work; φ is the electric potential.

4.3. The boundary and initial conditions

Thermal:

For the top surface of the structure, the aerodynamic heat flux and the radiation dissipation are the boundary conditions. For the bottom surface, the assumed convection heat transfer which can be considered as active cooling is the boundary condition. The influence of convective coefficients on the TPS performance is studied in this paper. For a unit cell model to fully represent the macro structure, the accurate boundary conditions are essential. For the model in this work, the size-reducing occurs in transverse directions, thus, for the side planes of the structure, the boundary conditions should be considered more carefully. For a unit cell formulated by translational symmetries, the boundary conditions should be the so-called periodic conditions. The mechanical and thermal boundary conditions for the unit cell model are shown in Eqs. (9) and (10), respectively. In the equations, ∇ T x0 and ∇ T 0y are the macroscopic temperature gradients in x- and y-direction, respectively, and equal to 0 for the heat conduction problem in this work, K is the coefficient of convection heat transfer and qaerodynamic is the aerodynamic heat flux as shown in Fig. 2. In fact, the periodic boundary conditions in Eqs. (9) and (10) describe the relative displacement relations and the relative temperature relations between corresponding nodes, respectively. The corresponding nodes mean two nodes translational symmetric to each other, e.g., M = (x1 , y 1 , z1 ) and M  = (x1 + 12, y 1 , z1 ) which are translational symmetric along x-axis. In some commercial FEM software like ANSYS and ABAQUS, the boundary conditions of nodes on edges and vertices should be imposed separately to avoid the possible over constraint. In Eqs. (9) and (10), the boundary planes are expressed as P 1 , P 2 , P 3 , P 4 , P 5 and P 6 , and the edges and vertices are indicated by numbers. The boundaries of the model are shown in the rightmost picture of Fig. 7, and sim-

Thermoelasticity: ∂T ∂{ε } ρc = [λ]∇ 2 T − ( T ref + T off ){β}T ∂t ∂t Thermoelectrics:

(6a)

(6b) j 2x + j 2y + j 2z ∂T = [λ]∇ 2 T + ∂t k where ρ is the density; c is the heat capacity; [λ] is the thermal conductivity matrix; {β}T is the transpose of vector {β} = [C ]{α }, in which [C ] and {α } are the stiffness matrix and the thermal ex-

ρc

pansion coefficient vector, respectively; T ref (= 300 K in this work) is the reference temperature; T off (= 0 in this work) is the offset temperature from absolute zero; j x , j y and j z are the current density in x, y and z directions, respectively. It should be noted that the second terms of Eqs. (6a) and (6b) represent Piezocaloric effect for the thermoelastic coupling and Joule heating for the thermoelectric coupling, respectively. In addition to the above equations, the constitutive equations are:

{σ } = [C ]{ε } − {β}( T − T ref )   { j } = [k] {φi } − [ S ]{∇ T }

(7) (8)

where “[ ]” represents matrix, “{ }” represents vector, and “∇ ” means gradient. 4.2. The unit cell model and its mesh As discussed above, the TPS structure can have an arbitrary size in transverse direction (x- and y-), thus in order to save the computational cost, a unit cell model is formulated.

464

C.-L. Gong et al. / Aerospace Science and Technology 77 (2018) 458–470

ilar schematic can be found in the authors’ previous work [40]. In the picture, the white circles with numbers are vertices while the black circles with numbers are edges. The model has a dimension of a × b × h which is 12 × 6 × 25 mm in this work as stated above. It should be noted that the periodic boundary condition is the so-called Dirichlet boundary condition which is related to nodes, so Eqs. (9) and (10) describe the constraint conditions of nodes on planes, on edges and on vertices. However, some equations in Eq. (10) involve heat flux which is Neumann boundary condition, so it should be pointed out that such heat flux in fact is imposed on corresponding elements rather than on those nodes.

Planes:

Planes: P 1 –P 2 : u (0, y 1 ,z1 ) − u (a, y 1 ,z1 ) = aε

Fig. 8. Interfaces and relevant thermal contact resistances.

0 x,

P 1 –P 2 : T (0, y 1 ,z1 ) − T (a, y 1 ,z1 ) = 0

v (0, y 1 ,z1 ) − v (a, y 1 ,z1 ) = 0,

P 3 –P 4 : T (x1 ,0,z1 ) − T (x1 ,b,z1 ) = 0

w (0, y 1 ,z1 ) − w (a, y 1 ,z1 ) = 0

P 5 : q z = −h( T − T ∞ )

P 3 –P 4 : u (x1 ,0,z1 ) − u (x1 ,b,z1 ) = 0,



v (x1 ,0,z1 ) − v (x1 ,b,z1 ) = bε 0y ,

4 4 P 6 : q z = qaerodynamic − σ ε T w − T∞

w (x1 ,0,z1 ) − w (x1 ,b,z1 ) = 0

Edges:

P5 : u = v = w = 0

9–17 : T (0,0,z1 ) − T (a,0,z1 ) = a∇ T x0

Edges:

9–11 : T (0,0,z1 ) −

9–17 : u (0,0,z1 ) − u (a,0,z1 ) = aεx0 ,

Elements on 10, 13, 16, 19 : q z = − K ( T − T ∞ )

w (0,0,z1 ) − w (a,0,z1 ) = 0 v (0,0,z1 ) − v (0,b,z1 ) = bε 0y ,

Elements related to 1, 2, 3, 4 : q z = − K ( T − T ∞ )

9–14 : u (0,0,z1 ) − u (a,b,z1 ) = aεx0 ,

4 Elements related to 5, 6, 7, 8 : q z = qaerodynamic − σ ε T w

v (0,0,z1 ) − v (a,b,z1 ) = bε 0y , w (0,0,z1 ) − w (a,b,z1 ) = 0 10, 13, 16, 19 : u = v = w = 0

v (0,0,0) − v (a,0,0) = 0, w (0,0,0) − w (a,0,0) = 0 v (0,0,0) − v (a,b,0) = bε 0y , w (0,0,0) − w (a,b,0) = 0 1–4 : u (0,0,0) − u (0,b,0) = 0, v (0,0,0) − v (0,b,0) = bε 0y , w (0,0,0) − w (0,b,0) = 0 1–5 : u (0,0,0) − u (0,0,h) = 0, v (0,0,0) − v (0,0,h) = 0, w (0,0,0) − w (0,0,h) = hεz0 1–6 : u (0,0,0) − u (a,0,h) = aεx0 , v (0,0,0) − v (a,0,h) = 0, w (0,0,0) − w (a,0,h) = hεz0 1–7 : u (0,0,0) − u (a,b,h) = aεx0 , v (0,0,0) − v (a,b,h) = bε 0y , w (0,0,0) − w (a,b,h) = hεz0 1–8 : u (0,0,0) − u (0,b,h) = 0, v (0,0,0) − v (0,b,h) = bε 0y , w (0,0,0) − w (0,b,h) = hεz0



Vertices:

w (0,0,z1 ) − w (0,b,z1 ) = 0

1–3 : u (0,0,0) − u (a,b,0) = aεx0 ,



4 4 Elements on 12, 15, 18, 20 : q z = qaerodynamic − σ ε T w − T∞

9–11 : u (0,0,z1 ) − u (0,b,z1 ) = 0,

1–2 : u (0,0,0) − u (a,0,0) = aεx0 ,

(10)

T (0,b,z1 ) = b∇ T 0y

9–14 : T (0,0,z1 ) − T (a,b,z1 ) = a∇ T x0 + b∇ T 0y

v (0,0,z1 ) − v (a,0,z1 ) = 0,

Vertices:



(9)

For the TE materials involved in the multi-functional structure, the mechanical-thermoelectric performance is considered, while for other materials excluding TE ones, only thermal-mechanical problems are analyzed. For the TE material couple, the ground point is set on the lower surface of n material. The unsteady simulation is conducted for a duration of 0 to 1425 s. The initial temperature and electric potential are 300 K and 0 V, respectively. 4.4. Interfaces between different materials There are five different materials and eight typical interfaces involved in the present structure as shown in Fig. 8. The thermal contact resistance (TCR) at the interface will influence the heat transfer process and thus brings about some effects on the calculation results. TCR always has a value range of 10−6 to 10−4 (K m2 /W). In this work, the TCRs of interfaces a, d and g are integrated as TCR1, and that of c, h, and f are integrated as TCR2. TCR1 and TCR2 are simulated by two additional layers with a thickness of 0.1 mm for each. The thermal conductivity of the layer should equal to the thickness divided by the assumed TCR, which means that if the TCR is assumed to be 10−4 , the thermal conductivity of the layer would be 1 W/(m K). In this paper, the influence of TCR (both TCR1 and TCR2) is considered by values of 0 and 10−4 (K m2 /W). It should be noted that for the case of TCR equals 0, no additional layers are simulated in the structure. Thermal contact resistance is a very complex issue related to interface topography, contact deformation, etc., and no quantitative description can be obtained for the problem studied in this work, thus most discussions in later Section 5 is based on the situation of no TCR consideration.

C.-L. Gong et al. / Aerospace Science and Technology 77 (2018) 458–470

465

Fig. 9. Temperature fields (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)

5. Results and discussions The mechanical, thermal and electrical performances of the structure are analyzed along the typical trajectory (see Fig. 1) from 0 to 1425 s. The influence of TCR and the bottom active cooling is considered with the convective coefficient values of 10, 50, 100, 300 W/(m2 K) in this section. 5.1. Thermo-mechanical performance Figs. 9(a) and (b) are temperature distributions of the structure under convective coefficients of 10 and 100 W/(m2 K), respectively. From left to right are results under 180, 495, 900 and 1400 s, respectively, and the corresponding trajectory parameters for each time are summarized in the top of the figure. For the legend, only the lowest and highest values are shown. As shown in Fig. 9(a), the temperature on the top surface increases from 302 to 918 K and then decreases to 720 K at the time of 1400 s, while that on the bottom surface increases from 300 to 635 K in the duration of 1400 s. However, the temperature on the top surface in Fig. 9(b) increases from 302 to 917 K and then decreases to 647 K, while that on the bottom surface increases from 300 to 423 K and then decreases to 394 K in a duration of 1400 s. Also, if the location of the red isothermal lines is concerned, one can find that the line moves downward with the increasing flight time, and the situation of higher convective coefficient moves more quickly. Although there must be several temperature ups and downs between these four typical times, this

Fig. 10. Temperatures on the top and the bottom surface of the structure.

phenomenon reveals that the active cooling under the bottom of the TPS will cool the structure in a shorter time. A more detailed temperature distribution variation along with the flight time can be found in Fig. 10, in which four cases of different convective coefficients are calculated. The temperatures of the top and bottom surfaces which are represented by the center nodes of each surface are displayed by black and red lines, respectively. The temperature on the top surface has two obvious valleys

466

C.-L. Gong et al. / Aerospace Science and Technology 77 (2018) 458–470

Fig. 11. The distribution of y-displacement.

Fig. 12. Results comparison between different boundary conditions ( y-displacement).

after 495 s, while that on the bottom surface has much more gentle trends, which means that the temperature on the top surface is more sensitive to the variation of input heat flux. The black rectangles and the red circles in Fig. 10 are the temperature on the top and bottom surfaces, respectively, under the consideration of TCR of 10−4 (K m2 /W). The convective coefficient is 10 W/(m2 K) in this case. TCR will result in the increasing of temperature on the top surface while the decreasing of that on the bottom surface, however, the influence is very limited with a largest deviation of 5% and an average deviation of 0.9%. The bottom surface of such TPS structure faces the vehicle structure, hence the bottom temperature is the criterion to validate the reliability of the TPS. For the studied four convective coefficients of 10, 50, 100 and 300 W/(m2 K), the highest temperature of the bottom surface is 639, 497, 431 and 364 K, respectively. These values provide references for the TPS design of vehicles with different internal temperature control demands. It should be noted that the structures shown in Fig. 9 are the deformed volumes although they are very similar to the original ones. The thermal deformation calculated based on periodic mechanical boundary conditions should be discussed. Take

y-displacement as an instance as shown in Fig. 11 in which the convective coefficient is 100 W/(m2 K) and the four contour plots are results at 180, 495, 900 and 1400 s, respectively. It is clear that the distribution of y-displacement indicates the thermal deformation in y direction. According to the legend, one can find that the red region is the expansion in positive y-direction, while the blue one is in negative y-direction. Fig. 12 shows the difference of results obtained from nonperiodic boundary condition. In the right diagram, the planes in y-direction has no specific boundary conditions while that in x-direction are still in periodic constraint. The results of 495 s under convective coefficient of 100 W/(m2 K) are shown in the figure. The volume results and that on boundary planes of P 3 ( y = 0) and P 4 ( y = b) are displayed, respectively. It is clear that the non-periodic boundary conditions will result in non-periodic displacement fields on symmetric boundaries, and the contour lines are in up-and-down strip pattern. The thermal expansion in insulation region is much larger than that in C/C-SiC and Ti-6Al-4V region, this is obviously unreasonable if the neighboring unit cells are considered, i.e., the translational movement of such deformed unit cells could not compose a complete structure since the C/C-SiC

C.-L. Gong et al. / Aerospace Science and Technology 77 (2018) 458–470

467

Fig. 13. The distribution of electric potential.

(and Ti-6Al-4V) regions are not continuous. However, the deformed structures that based on periodic boundary conditions can form a complete deformed TPS structure. This reveals that the periodic boundary conditions are essential to obtaining reasonable results for such unit cell. 5.2. Thermoelectric performance Fig. 13 is the electric potential distribution of TE materials, and Figs. 13(a) and (b) are results under convective coefficients of 10 and 100 W/(m2 K), respectively. The electrical potential of the module increases from 0 at the bottom of n material to the highest value at the bottom of p material. Therefore, the largest value noted in the legend is the total electric potential of the module. For the convective coefficient of 10 W/(m2 K), the electric potential increases from 0.0003 to 0.144 and then decreases to 0.024 V at the time of 1400 s, while for the case of 100 W/(m2 K), the potential increases from 0.0003 to 0.149 and then decreases to 0.065 V. Detailed results are shown in Fig. 14, in which the black solid line, red dashed line, blue dash-dotted line and the pink dotted line are the total potential of cases 10, 50, 100 and 300 W/(m2 K), respectively. Two valleys can be found in each line after 495 s, at which the highest output is obtained. The maximum potentials of each line are 0.144, 0.147, 0.149 and 0.154 V, respectively, while the average values for each case in the 1400 s duration are 0.048, 0.063, 0.69 and 0.076 V, respectively. For a module of 0.01 m2 (100 × 100 mm) with 128 p–n couples, the maximum electric potential can be 18.4 to 19.7 V. It is clear that the higher convective coefficient obtains a larger electric potential, this should be due to the larger temperature difference across the TE materials. In

Fig. 14. The total electric potential.

Fig. 14, black rectangles represent the results considering the influence of TCR (10−4 (K m2 /W)) with the convective coefficient of 10 W/(m2 K). It reveals in the figure that TCR has little effects on the total electric potentials with a maximum deviation of 4% and an average deviation of 1%. Fig. 15 shows the temperature and electric potentials on the top and bottom surfaces of n and p materials. In the figure, the black solid lines and the red dashed lines are the temperature on the top surface and the bottom surface of TE materials, and the blue solid lines and the pink dashed lines are the total electrical po-

468

C.-L. Gong et al. / Aerospace Science and Technology 77 (2018) 458–470

Fig. 15. The thermo-electric results of p and n couple.

Fig. 16. The maximum output power and TE conversion efficiency.

tential values in Fig. 14 and that of n material, respectively. In the upper part of the figure, the black solid lines and the red dashed lines are the temperature on the top surface and bottom surface of TE materials, respectively. It should be noted that each line in fact contains two coincident lines of p and n materials, respectively, and it means that p and n have very close temperature distributions which can be due to that both of their thermal conductivities are much larger than that of insulation material. In Fig. 15, for the results obtained with consideration of TCR (10−4 (K m2 /W)), the black rectangles and the red circles represent the temperature on the top and the bottom surface of p-material, while the blue rectangles and the pink circles represent the electrical potentials of p- and n-material, respectively. The convective coefficient is 10 W/(m2 K). Very small deviation can be found between results with TCR and without TCR: the maximum and average deviation of temperature are 0.6% and 0.2%, and that of electrical potential are 4% and 1%, respectively. 5.3. Thermoelectric conversion efficiency The TE module can be used as a power supply, and its output power is related to the electrical load. In this work, the largest power for ideal conditions in which the load resistance equals to the resistance of TE materials are evaluated as shown in Fig. 16. In the figure, the solid lines indicate the output power, while the dashed lines are the TE conversion efficiency, and the black, red, blue and pink lines are results for different convective coefficients, respectively. The maximum output power is calculated by E 2 /(4r ),

where E is the total electric potential (see Fig. 14) while r is the electric resistance of the TE couple. In this work, r varies with flight time (changing temperature), and the values under average temperatures (( T top− p&n + T bottom− p&n )/2, see Fig. 15) are used in the calculation. The TE conversion efficiency is calculated by the aerodynamic heat divided by the output power, and the aerodynamic heat is account for the area of the unit cell model rather than only the TE materials. It has to be pointed out that in a real engineering condition, the electric load is always steady and larger than the resistance of TE materials, and thus the output power should be less than the ideal value shown in Fig. 16. As shown in Fig. 16, the largest power that a p–n couple can supply increases from 0 V to its peak value and then decreases to a smaller value during the mission. The TE conversion efficiency has very similar trend excluding the end of those lines, which are in a sharp increasing. This is because in an unsteady simulation as this work, when the input heat flux decreases in a very short time, the temperature difference that can be used to create electric energy will maintain for a while. The influence of TCR on the TE conversion can also be found in Fig. 16, in which the black rectangles and the black circles are the largest output power and the TE conversion efficiency, respectively. The maximum deviation is 9% and the average deviation is only 2%, which means that under the typical vehicle and trajectory studied in this work, the TCR at different material interfaces has limited influence on the conversion efficiency of the developed TE module. In addition, both the output power and the conversion efficiency increase with a larger conductive coefficient. As shown in Table 3, for cases of K = 10, 50, 100 and 300 W/(m2 K), the largest and the average values of output power and TE conversion efficiency are listed. The largest and the average power increase from 0.0182 to 0.0209 W and 0.00317 to 0.00713 W, respectively, while the largest and the average conversion efficiency increases from 2.63% to 4.09% and 0.237% to 0.607%, respectively. It means that the active cooling under the TPS bottom is very useful in improving the TE performance. It also worth noting that for a TPS structure (128 p–n couples) of just 0.01 m2 , the largest transient power and the average power will be 2.33 to 2.68 W and 0.406 to 0.913 W, respectively. The area of the compression surface of a typical reusable launch vehicle would be about 30 m2 as discussed in Section 1, which means that the largest transient output power and the average power during the mission for the studied vehicle and trajectory in this work can be higher than 6000 and 1000 W, respectively. Therefore, such kind of TPS structure has potentials to provide electrical power of thousands Watt and could be a part substitution of power supply especially for the vehicle that has a long time hypersonic flight. Conclusions In this work, a thermoelectric material based multifunctional TPS structure is proposed and the relevant evaluation approach is developed. The TE module in the structure can convert a certain amount of aerodynamic heat into electricity supply. The TE couple consists of Ca3 Co4 O9 based p-type material, and Sr0.9 La0.1 TiO3 based n-type compound which is fabricated by solid state reaction method. The mechanical-thermoelectric performance of the structure is studied based on a unit cell model. The aerodynamic heat is calculated by an engineering-based method based on a typical reusable launch vehicle with a typical trajectory curve, and is used as the input boundary conditions for the numerical simulation. The simulation approach developed in this work can be used in analysis of other trajectories of hypersonic vehicles. The results show some conclusions:

C.-L. Gong et al. / Aerospace Science and Technology 77 (2018) 458–470

469

Table 3 The maximum output power and the TE conversion efficiency. Convective coefficient (W/(m2 K))

10 50 100 300

Output power (W) Unit cell (single p–n couple)

TE conversion efficiency (%) TPS structure (128 p–n couples)

Maximum

Average

Maximum

Average

Maximum

Average

0.0182 0.0189 0.0195 0.0209

0.00317 0.00483 0.00591 0.00713

2.33 2.42 2.50 2.68

0.406 0.618 0.756 0.913

2.63 2.79 3.39 4.09

0.237 0.400 0.499 0.607

1. In the thermal-mechanical studies of a unit cell model that formulated based on translational symmetries, the periodic boundary condition is essential to obtaining reasonable results. 2. For the TPS structure proposed in this work, if a convection heat transfer with coefficients of 10, 50, 100 and 300 W/(m2 K) are provided under the bottom, in the studied trajectory duration the highest temperature of the bottom surface is 639, 497, 431 and 364 K, respectively, and this provide references for the TPS design of vehicles with different internal temperature control demands. 3. In the studied trajectory duration, a multi-functional TPS structure of 0.01 m2 can provide maximum average output powers of 0.406, 0.618, 0.756 and 0.913 W, respectively, and considering the large area of real vehicle, such TPS has potentials to provide a part substitution of power supply. 4. Active cooling is very useful in improving the thermoelectric performance including the output power and the TE conversion efficiency. Conflict of interest statement The authors declared that they have no conflict of interest to this work. Acknowledgements This study is supported by the China Civil Aerospace Program (No. D010403), the National Defense Fundamental Research Funds of China (No. JCKY2016204B102), the Fundamental Research Funds for the Central Universities (3102017OQD068) and the National Natural Science Foundation of China (Nos. 51672219, 51702259). References [1] B. Behrens, M. Müller, Technologies for thermal protection systems applied on re-usable launcher, Acta Astronaut. 55 (3–9) (2004) 529–536. [2] F. Gori, S. Corasaniti, W.M. Worek, W.J. Minkowycz, Theoretical prediction of thermal conductivity for thermal protection systems, Appl. Therm. Eng. 49 (2012) 124–130. [3] Y. Rong, Y. Wei, R. Zhan, Research on thermal protection by opposing jet and transpiration for high speed vehicle, Aerosp. Sci. Technol. 48 (Supplement C) (2016) 322–327. [4] D. Knight, H. Yan, A.G. Panaras, A. Zheltovodov, Advances in CFD prediction of shock wave turbulent boundary layer interactions, Prog. Aerosp. Sci. 39 (2) (2003) 121–184. [5] D. Knight, J. Longo, D. Drikakis, D. Gaitonde, A. Lani, I. Nompelis, B. Reimann, L. Walpot, Assessment of CFD capability for prediction of hypersonic shock interactions, Prog. Aerosp. Sci. 48 (Supplement C) (2012) 8–26. [6] Y. Kerboua, A.A. Lakis, Numerical model to analyze the aerodynamic behavior of a combined conical–cylindrical shell, Aerosp. Sci. Technol. 58 (Supplement C) (2016) 601–617. [7] P. Panagiotou, P. Kaparos, C. Salpingidou, K. Yakinthos, Aerodynamic design of a MALE UAV, Aerosp. Sci. Technol. 50 (Supplement C) (2016) 127–138. [8] P.D. Bravo-Mosquera, L. Botero-Bolivar, D. Acevedo-Giraldo, H.D. Cerón-Muñoz, Aerodynamic design analysis of a UAV for superficial research of volcanic environments, Aerosp. Sci. Technol. 70 (Supplement C) (2017) 600–614. [9] C.D. Engel, S.C. Praharaj, MINIVER Upgrade for the AVID System. Vol. 1: LANMIN User’s Manual, Nasa Cr, 1983. [10] E.V. Zoby, K. Sutton, Approximate convective-heating equations for hypersonic flows, J. Spacecr. Rockets 18 (1) (1981) 64.

[11] Y. Ma, B. Xu, M. Chen, R. He, W. Wen, T. Cheng, D. Fang, Optimization design of built-up thermal protection system based on validation of corrugated core homogenization, Appl. Therm. Eng. 115 (2017) 491–500. [12] S. Kumar, S.P. Mahulikar, Reconstruction of aero-thermal heating and thermal protection material response of a reusable launch vehicle using inverse method, Appl. Therm. Eng. 103 (Supplement C) (2016) 344–355. [13] R. Savino, M. De Stefano Fumo, D. Paterna, M. Serpico, Aerothermodynamic study of UHTC-based thermal protection systems, Aerosp. Sci. Technol. 9 (2) (2005) 151–160. [14] S. Liu, B. Zhang, Effects of active cooling on the metal thermal protection systems, Aerosp. Sci. Technol. 15 (7) (2011) 526–533. [15] R.J. Tobe, R.V. Grandhi, Hypersonic vehicle thermal protection system model optimization and validation with vibration tests, Aerosp. Sci. Technol. 28 (1) (2013) 208–213. [16] G. Xie, Q. Wang, B. Sunden, W. Zhang, Thermomechanical optimization of lightweight thermal protection system under aerodynamic heating, Appl. Therm. Eng. 59 (1) (2013) 425–434. [17] A. Riccio, F. Raimondo, A. Sellitto, V. Carandente, R. Scigliano, D. Tescione, Optimum design of ablative thermal protection systems for atmospheric entry vehicles, Appl. Therm. Eng. 119 (Supplement C) (2017) 541–552. [18] M. Cioeta, G. Di Vita, T. Signorelli Maria, G. Bianco, M. Cutroni, F. Damiani, V. Ferretti, A. Rotondo, Design, qualification, manufacturing and integration of IXV ablative thermal protection system, Acta Astronaut. 124 (Supplement C) (2016) 90–101. [19] J. Barcena, I. Garmendia, K. Triantou, K. Mergia, B. Perez, S. Florez, G. Pinaud, J-M. Bouilly, W.P.P. Fischer, Infra-red and vibration tests of hybrid ablative/ceramic matrix technological breadboards for Earth re-entry thermal protection systems, Acta Astronaut. 134 (Supplement C) (2017) 85–97. [20] C. Gayner, K.K. Kar, Recent advances in thermoelectric materials, Prog. Mater. Sci. 83 (Supplement C) (2016) 330–382. [21] S. Lv, W. He, Q. Jiang, Z. Hu, X. Liu, H. Chen, M. Liu, Study of different heat exchange technologies influence on the performance of thermoelectric generators, Energy Convers. Manag. 156 (Supplement C) (2018) 167–177. [22] S. Lan, Z. Yang, R. Chen, R. Stobart, A dynamic model for thermoelectric generator applied to vehicle waste heat recovery, Appl. Energy 210 (Supplement C) (2018) 327–338. [23] X. Li, Z. Wang, Exergy analysis of integrated TEG and regenerative cooling system for power generation from the scramjet cooling heat, Aerosp. Sci. Technol. 66 (Supplement C) (2017) 12–19. [24] K. Cheng, J. Qin, Y. Jiang, C. Lv, S. Zhang, W. Bao, Performance assessment of multi-stage thermoelectric generators on hypersonic vehicles at a large temperature difference, Appl. Therm. Eng. (2017), https://doi.org/10.1016/j. applthermaleng.2017.11.057. [25] E. Guilmeau, R. Funahashi, M. Mikami, K. Chong, D. Chateigner, Thermoelectric properties–texture relationship in highly oriented Ca3 Co4 O9 composites, Appl. Phys. Lett. 85 (9) (2004) 1490–1492. [26] R. Funahashi, M. Mikami, T. Mihara, S. Urata, N. Ando, A portable thermoelectric-power-generating module composed of oxide devices, J. Appl. Phys. 99 (6) (2006) 066117. [27] S. Urata, R. Funahashi, T. Mihara, A. Kosuga, S. Sodeoka, T. Tanaka, Power generation of p-type Ca3 Co4 O9 /n-type CaMnO3 module, in: 25th International Conference on Thermoelectrics, 2007. [28] C.L. Gong, B. Chen, L.X. Gu, Comparison study of RBCC powered suborbital reusable launch vehicle concepts, in: 20th AIAA International Space Planes and Hypersonic Systems and Technologies Conference, American Institute of Aeronautics and Astronautics, 2015. [29] V. Fernandez Villace, J. Steelant, The thermal paradox of hypersonic cruisers, in: 20th AIAA International Space Planes and Hypersonic Systems and Technologies Conference, American Institute of Aeronautics and Astronautics, 2015. [30] M. Moore, J. Williams, Aerodynamic prediction rationale for analyses of hypersonic configurations, in: 27th Aerospace Sciences Meeting, American Institute of Aeronautics and Astronautics, 1989. [31] I. Terasaki, Large thermoelectric power in NaCo2 O4 single crystals single crystals, Phys. Rev. B 56 (20) (1997) 12685–12687. [32] M. Shikano, R. Funahashi, Electrical and thermal properties of single-crystalline (Ca2 CoO3 )0.7 CoO2 with a Ca3 Co4 O9 structure, Appl. Phys. Lett. 82 (12) (2003) 1851–1853.

470

C.-L. Gong et al. / Aerospace Science and Technology 77 (2018) 458–470

[33] Y. Liu, Y. Lin, S. Zhan, C.W. Nan, Z. Shen, Preparation of Ca3 Co4 O9 and improvement of its thermoelectric properties by Spark Plasma Sintering, J. Am. Ceram. Soc. 88 (5) (2005) 1337–1340. [34] J.G. Noudem, M. Prevel, A. Veres, D. Chateigner, J. Galy, Thermoelectric Ca3 Co4 O9 ceramics consolidated by Spark Plasma sintering, J. Electroceram. 22 (1–3) (2009) 91–97. [35] N.N. Van, N. Pryds, S. Linderoth, M. Ohtaki, Enhancement of the thermoelectric performance of p-type layered oxide Ca3 Co4 O(9+δ ) through heavy doping and metallic nanoinclusions, Adv. Mater. 23 (21) (2011) 2484–2490. [36] C.L. Gong, G. Dong, J. Hu, Y. Chen, M. Qin, S. Yang, F. Gao, Effect of reducing annealing on the microstructure and thermoelectric properties of La–Bi co-doped SrTiO3 ceramics, J. Mater. Sci., Mater. Electron. 28 (19) (2017) 14893–14900.

[37] Y. Wang, Y. Sui, J. Cheng, X. Wang, W. Su, Comparison of the high temperature thermoelectric properties for Ag-doped and Ag-added Ca3 Co4 O9 , J. Alloys Compd. 477 (1) (2009) 817–821. [38] N. Van Nong, N. Pryds, S. Linderoth, M. Ohtaki, Enhancement of the thermoelectric performance of p-Type layered oxide Ca3 Co4 O9+δ through heavy doping and metallic nanoinclusions, Adv. Mater. 23 (21) (2011) 2484–2490. [39] J.J. Gou, X.J. Ren, Y.J. Dai, S. Li, W.Q. Tao, Study of thermal contact resistance of rough surfaces based on the practical topography, Comput. Fluids 164 (2018) 2–11, available online in 2016. [40] J.J. Gou, C.L. Gong, L.X. Gu, S. Li, W.Q. Tao, Unit cells of composites with symmetric structures for the study of effective thermal properties, Appl. Therm. Eng. 126 (Supplement C) (2017) 602–619.