Full Text:   <3023>

Summary:  <1848>

CLC number: V21

On-line Access: 2014-03-04

Revision Accepted: 2014-01-17

Crosschecked: 2014-02-20

Cited: 5

Clicked: 6329

Citations:  Bibtex RefMan EndNote GB/T7714

 Journal of Zhejiang University SCIENCE A 2014 Vol.15 No.3 P.185-196 http://doi.org/10.1631/jzus.A1300341

Numerical simulation of aerodynamic heating and stresses of chemical vapor deposition ZnS for hypersonic vehicles*

 Author(s):  Yuan-chun Liu1, Yu-rong He1, Jia-qi Zhu2, Jie-cai Han2, Dong-liang Quan3 Affiliation(s):  1. School of Energy Science & Engineering, Harbin Institute of Technology, Harbin 150001, China; more Corresponding email(s):   rong@hit.edu.cn Key Words:  Chemical vapor deposition (CVD) ZnS, Infrared window material, Thermal and stress responses, Hypersonic vehicles Share this article to： More <<< Previous Article|Next Article >>>

Yuan-chun Liu, Yu-rong He, Jia-qi Zhu, Jie-cai Han, Dong-liang Quan. Numerical simulation of aerodynamic heating and stresses of chemical vapor deposition ZnS for hypersonic vehicles[J]. Journal of Zhejiang University Science A, 2014, 15(3): 185-196.

@article{title="Numerical simulation of aerodynamic heating and stresses of chemical vapor deposition ZnS for hypersonic vehicles",
author="Yuan-chun Liu, Yu-rong He, Jia-qi Zhu, Jie-cai Han, Dong-liang Quan",
journal="Journal of Zhejiang University Science A",
volume="15",
number="3",
pages="185-196",
year="2014",
publisher="Zhejiang University Press & Springer",
doi="10.1631/jzus.A1300341"
}

%0 Journal Article
%T Numerical simulation of aerodynamic heating and stresses of chemical vapor deposition ZnS for hypersonic vehicles
%A Yuan-chun Liu
%A Yu-rong He
%A Jia-qi Zhu
%A Jie-cai Han
%A Dong-liang Quan
%J Journal of Zhejiang University SCIENCE A
%V 15
%N 3
%P 185-196
%@ 1673-565X
%D 2014
%I Zhejiang University Press & Springer
%DOI 10.1631/jzus.A1300341

TY - JOUR
T1 - Numerical simulation of aerodynamic heating and stresses of chemical vapor deposition ZnS for hypersonic vehicles
A1 - Yuan-chun Liu
A1 - Yu-rong He
A1 - Jia-qi Zhu
A1 - Jie-cai Han
A1 - Dong-liang Quan
J0 - Journal of Zhejiang University Science A
VL - 15
IS - 3
SP - 185
EP - 196
%@ 1673-565X
Y1 - 2014
PB - Zhejiang University Press & Springer
ER -
DOI - 10.1631/jzus.A1300341

Abstract:
hypersonic vehicles subjected to strong aerodynamic forces and serious aerodynamic heating require more stringent design for an infrared window. In this paper, a finite element analysis is used to present the distributions of thermal and stress fields in the infrared window for hypersonic vehicles based on flowfield studies. A theoretical guidance is provided to evaluate the influence of aerodynamic heating and forces on infrared window materials. The aerodynamic heat flux from Mach 3 to Mach 6 flight at an altitude of 15 km in a standard atmosphere is obtained through flowfield analysis. The thermal and stress responses are then investigated under constant heat transfer coefficient boundary conditions for different Mach numbers. The numerical results show that the maximum stress is higher than the material strength at Mach 6, which means a failure of the material may occur. The maximum stress and temperatures are lower than the material strength and melting point under other conditions, so the material is safe.

## 1.  Introduction

Aircraft and spacecraft structures designed for supersonic and hypersonic flights are subjected to severe aerodynamic heating during the launch and reentry phases of their operations, which is caused by the air in the boundary layer being progressively slowed down (Ruan et al., ). Consequently, all external surfaces on the vehicle are heated. This in turn leads to non-uniform transient temperatures that produce dynamic thermal stresses and deformations. Therefore, high heating associated with shocks at leading edges are significant issues in the vehicle design. In addition to surface melting and ablation, the flight aerodynamics can be perturbed, leading to unacceptable flight trajectory deviations. Another problem is signal refraction passing through the shocked hot gas layer in front of the vehicle’s head (Saravanan et al., ). In recent years, there have been significant investments in the development of hypersonic vehicle technologies. Hypersonic flight began in February 1949 when a Women’s Army Corps (WAC) corporal rocket was ignited from a US captured V-2 rocket (Sun and Wu, ). Later, extensive numerical analysis (Jain and Hayes, ; Di Clemente et al., ; Gerdroodbary and Hosseinalipour, ) obtained the pressure, heat transfer, and surface temperatures on steady or unsteady heat transfer in boundary layers of hypersonic flow. Although some flight experiments are also conducted (Di Clemente et al., ; Marini et al., ) to acquire aerodynamic heating data in the flight condition to evaluate these prediction methods, flight data are not adequate for complete validations. The coupled approach between the external and internal fields in hypersonic conditions has been studied in recent years (Culler and McNamara, ). Di Clemente et al. () studied the surface temperature measurements in a plasma wind tunnel experiment and compared these results with the results of the test numerical rebuilding carried out through an integrated procedure coupling the external aerodynamic field to the internal thermal state of the structure.

As one of the key components of a hypersonic vehicle, an infrared window is used to transmit infrared signals, and to keep the aerodynamic shape and protect its imaging system. High levels of aerodynamic heating can cause a malfunction or even damage to the delicate onboard equipment (White, ; Gnemmi et al., ). Excessive heating can cause ablation to the vehicle surface material. Combined with the presence of high-pressure loads, severe heating can cause a complete material failure and affect the normal performance of other instruments underneath the window (Heubner et al., ). The window design requires (Russell et al., ): (1) accommodating the severe highly transient aeroheating and resulting thermal gradients; (2) providing an effective seal to prevent gas flow into the laser guidance electronics volume; (3) minimizing aerodynamic loads by surface flush; (4) minimizing dynamic effects on the window due to launch and flight shocks; (5) providing an adequate surface finish and transmissibility for the laser wavelength of interest.

Zinc sulphide (ZnS) is one of the most common materials for infrared windows (Harris et al., ). The most common technique to deposit cubic polycrystalline ZnS thick films for optical window applications is based on a chemical vapor deposition (CVD) process (Harris, ). This technique is faster and cheaper than the epitaxial methods (molecular beam epitaxy (MBE) and atomic layer epitaxy (ALE)), and has better purity than the wet chemical method (Tropf et al., ).

The main objective of this paper is to present a theoretical guidance for the design of the infrared window (CVD ZnS) for hypersonic vehicles. To accomplish this, a process for predicting the aeroheating and thermal response of typical vehicle configurations in high-speed flows is described first. The aerodynamic heat flux from Mach 3 to Mach 6 at an altitude of 15 km in a standard atmosphere is obtained through flowfield analysis. Then, a finite element analysis using a coupling method is described to study the distribution of thermal and stress fields of the infrared window for hypersonic vehicles. The remainder of this paper is arranged as follows. In the next section, the mathematical model used is illustrated including the flowfield model, the structural model, and the Von Mises theory. Next, the main results of the simulation study are discussed. The paper finalizes with the main conclusions and recommendations for future studies.

## 2.  Mathematical model

### 2.1.1.  Basic physical models for fluid flow

For all flows, the conservation equations for mass and momentum are solved. For flows involving heat transfer or compressibility, an additional equation for energy conservation is solved. Additional transport equations are also solved when the flow is turbulent.

The equation for conservation of mass, or continuity equation, can be written as follows (Shahmardan et al., ; Tang et al., ): $$\frac{{\partial \rho }}{{\partial t}} + \nabla \cdot (\rho v) = 0$$, where ρ is the density, t is the time, and v is the velocity.

For 2D geometries, the continuity equation is given by $$\frac{{\partial \rho }}{{\partial t}} + \frac{\partial }{{\partial x}}(\rho {v_x}) + \frac{\partial }{{\partial y}}(\rho {v_y}) = 0$$.

The conservation of momentum in an inertial (non-accelerating) reference frame is described by $$\frac{\partial }{{\partial t}}(\rho v) + \nabla \cdot (\rho vv) = - \nabla p + \nabla \cdot (\overline{\overline \tau } ) + \rho g + F$$, where p is the pressure, g is the gravitational acceleration (=9.80665 m/s2), and F is the force vector. The stress tensor $$\overline{\overline \tau }$$ is given by $$\overline{\overline \tau } = \mu \left[ {(\nabla v + \nabla {v^T}) - \frac{2}{3}\nabla \cdot \nu I} \right]$$, where μ is the molecular viscosity, and I is the unit tensor.

For 2D geometries, the axial and radial momentum conservation equations are given by $$\begin{array}{*{20}{c}} {\frac{\partial }{{\partial t}}\left( {\rho {v_x}} \right) + \frac{\partial }{{\partial x}}\left( {\rho {v_x}{v_x}} \right) + \frac{\partial }{{\partial y}}\left( {\rho {v_y}{v_x}} \right) = - \frac{{\partial p}}{{\partial x}} + \frac{\partial }{{\partial x}}\left[ {\mu \left( {2\frac{{\partial {v_x}}}{{\partial x}} - \frac{2}{3}\left( {\nabla \cdot v} \right)} \right)} \right] + \frac{\partial }{{\partial y}}\left[ {\mu \left( {\frac{{\partial {v_x}}}{{\partial y}} + \frac{{\partial {v_y}}}{{\partial x}}} \right)} \right] + {F_x},} \\ {\frac{\partial }{{\partial t}}\left( {\rho {v_y}} \right) + \frac{\partial }{{\partial x}}\left( {\rho {v_x}{v_y}} \right) + \frac{\partial }{{\partial y}}\left( {\rho {v_y}{v_y}} \right) = - \frac{{\partial p}}{{\partial y}} + \frac{\partial }{{\partial x}}\left[ {\mu \left( {\frac{{\partial {v_y}}}{{\partial x}} + \frac{{\partial {v_x}}}{{\partial y}}} \right)} \right] + \frac{\partial }{{\partial y}}\left[ {\mu \left( {2\frac{{\partial {v_y}}}{{\partial y}} - \frac{2}{3}\left( {\nabla \cdot v} \right)} \right)} \right] + {F_y},} \end{array}$$ where $$\nabla \cdot v = \frac{{\partial {v_x}}}{{\partial x}} + \frac{{\partial {v_y}}}{{\partial y}}$$.

The energy equation for 2D geometries described by temperature can be written as follows: $$\rho {C_p}\left( {\frac{{\partial T}}{{\partial t}} + {v_x}\frac{{\partial T}}{{\partial x}} + {v_y}\frac{{\partial T}}{{\partial y}}} \right) = \lambda \left( {\frac{{{\partial ^2}T}}{{\partial {x^2}}} + \frac{{{\partial ^2}T}}{{\partial {y^2}}}} \right)$$, where C p is the heat capacity at a constant pressure, T is the temperature, and λ is the thermal conductivity.

The continuity equation, momentum conservation equation, and energy equation can be described by the general differential equation: $$\frac{\partial }{{\partial t}}\left( {\rho \phi } \right) + \nabla \cdot \left( {\rho \phi v} \right) = \nabla \cdot \left( {{\Gamma _\phi }\nabla \phi } \right) + {S_\phi }$$, where ϕ is the general variable, and stands for the solving variables, such as vx , vy , and T. Γϕ is the diffusion coefficient for ϕ, and Sϕ is the source of ϕ.

For each control volume V, the discretization of the general differential equation can be written as follows: $$\int_V {\frac{{\partial (\rho \phi )}}{{\partial t}}dV} + \oint {\rho \phi v \cdot dA} = \oint {{\Gamma _\phi }\nabla \phi \cdot A} + \int_V {{S_\phi }dV}$$, where A is the surface area vector.

For 2D geometries, the discretization for the cell face f on a given cell yields is $$\frac{{\partial (\rho \phi )}}{{\partial t}}V + \sum\limits_f^{{N_{faces}}} {{\rho _f}{v_f}{\phi _f}} \cdot {A_f} = \sum\limits_f^{{N_{faces}}} {{\Gamma _\phi }\nabla {\phi _f}} \cdot {A_f} + {S_\phi }V$$, where N faces is the number of faces enclosing cell.

### 2.1.2.  Turbulence model

The shear-stress transport (SST) k-ω model was developed by Menter () to effectively blend the robust and accurate formulation of the k-ω model in the near-wall region with the free-stream independence of the k-ω model in the far field.

The SST k-ω model has a similar form to the standard k-ω model. The turbulence kinetic energy, k, and the specific dissipation rate, ω, are obtained from the following transport equations: $$\begin{array}{*{20}{c}} {\frac{\partial }{{\partial t}}\left( {\rho k} \right) + \frac{\partial }{{\partial x}}\left( {\rho k{v_x}} \right) + \frac{\partial }{{\partial y}}\left( {\rho k{v_y}} \right)\; = \frac{\partial }{{\partial x}}\left( {{\Gamma _k}\frac{{\partial k}}{{\partial x}}} \right) + \frac{\partial }{{\partial y}}\left( {{\Gamma _k}\frac{{\partial k}}{{\partial y}}} \right) + {{\tilde G}_k} - {Y_k} + {S_k},} \\ {\frac{\partial }{{\partial t}}\left( {\rho \omega } \right) + \frac{\partial }{{\partial x}}\left( {\rho \omega {v_x}} \right) + \frac{\partial }{{\partial y}}\left( {\rho \omega {v_y}} \right)\; = \frac{\partial }{{\partial x}}\left( {{\Gamma _\omega }\frac{{\partial \omega }}{{\partial x}}} \right) + \frac{\partial }{{\partial y}}\left( {{\Gamma _\omega }\frac{{\partial \omega }}{{\partial y}}} \right) + {{\tilde G}_\omega } - {Y_\omega } + {D_\omega } + {S_\omega }.} \end{array}$$

The effective diffusivities for the SST k-ω model are given by $$\begin{array}{*{20}{c}} {{\Gamma _k} = \mu + {\mu _t}/{\sigma _k},} \\ {{\Gamma _\omega } = \mu + {\mu _t}/{\sigma _\omega },} \end{array}$$ where σk is the turbulent Prandtl numbers for k, σω is the turbulent Prandtl numbers for ω, and μ t is the turbulent viscosity.

The term $${\tilde G_k}$$ represents the production of turbulence kinetic energy, and is defined as $${\tilde G_k} = \min ({G_k},\;10\rho {\beta ^ * }k\omega )$$, where $${\beta ^ * } = \beta _i^*\left[ {1 + {\zeta ^*}F({M_t})} \right]$$, where $$F({M_t})$$ is the compressibility function, and $$\beta _i^* = \beta _\infty ^*\left( {\frac{{4/15 + {{\left( {R{e_t}/{R_\beta }} \right)}^4}}}{{1 + {{\left( {R{e_t}/{R_\beta }} \right)}^4}}}} \right),R{e_t} = \rho k/(\mu \omega ),{\zeta ^*} = 1.5, {R_\beta } = 8, \beta _\infty ^* = 0.09$$.

The term $${G_\omega }$$ represents the production of ω, which is given by $${G_\omega } = (\iota /{v_t}){\tilde G_k}$$, where ι is a coefficient, in the high-Reynolds number form, ι=1, and ν t is the turbulent kinematic viscosity.

The term Yk represents the dissipation of turbulence kinetic energy, $${Y_k} = \rho {\beta ^ * }k\omega$$.

The term Yω represents the dissipation of ω, $${Y_\omega } = \rho \beta {\omega ^2}$$, where $$\beta = {\beta _i}\left[ {1 - (\beta _i^*/{\beta _i}){\zeta ^*}F({M_t})} \right]$$, $${\beta _i} = 0.072$$.

The introduction of a cross-diffusion term Dω is defined as $${D_\omega } = 2(1 - {F_1})\rho {\sigma _{\omega ,2}}\frac{1}{\omega }\left( {\frac{{\partial k}}{{\partial x}}\frac{{\partial \omega }}{{\partial y}} + \frac{{\partial k}}{{\partial y}}\frac{{\partial \omega }}{{\partial x}}} \right)$$, where F 1 is the blending function, and $${\sigma _{\omega ,2}} = 1.168$$.

### 2.2.  Structural model

2D geometries are used for the axisymmetric hypersonic flow in the fluid model. Therefore, the schematic quarter 3D model is used in the structural model to predict the temperature and stress distribution.

Based on the principle of energy balance, the solid wall temperature can be described by a differential relation. In a Cartesian coordinate system, the unsteady heat conduction differential equation without internal heat source is expressed as $$\rho {C_p}\frac{{\partial T}}{{\partial t}} = \lambda \left( {\frac{{{\partial ^2}T}}{{\partial {x^2}}} + \frac{{{\partial ^2}T}}{{\partial {y^2}}} + \frac{{{\partial ^2}T}}{{\partial {z^2}}}} \right)$$.

The boundary condition on the internal surface is adiabatic, while the boundary condition on the external surface is set as a constant thermal flux expressed as $$- \lambda {\left. {\frac{{\partial T}}{{\partial n}}} \right|_{{\tau _{\text{c}}}}} = h({T_s} - {T_a})$$, where $${\left. {\frac{{\partial T}}{{\partial n}}} \right|_{{\tau _c}}}$$ is the normal temperature gradient of a solid surface, h is the heat transfer coefficient, T a is the ambient temperature, and T s is the temperature of the window.

For the finite element analysis, the functional equation established by the weighted residual quantity is $$\int_V {\left( {\rho {C_p}\frac{{\partial T}}{{\partial t}} - \lambda {\nabla ^2}T} \right)\delta \varpi dV} - \int_S {\left( {\lambda \frac{{\partial T}}{{\partial n}} + {{\bar q}_n}} \right)\delta \varpi dS} = 0$$, where $${\bar q_n} = h({T_s} - {T_a})$$, ϖ is weight function, V is the volume of the element, and S is the boundary of the element.

Transforming the formula by Gauss-Green, integration is obtained as $$\int_V {\left[ {\left( {\rho {C_p}\frac{{\partial T}}{{\partial t}}\delta \varpi + \lambda \nabla T \cdot \nabla \delta \varpi } \right)} \right]dV} - \int_S {{{\bar q}_n}\delta \varpi dS} = 0$$.

The discretization of the integral domain into elements and the temperature elements are interpolated as follows: $$T = N{T_e},\varpi = N{\varpi _e},\delta \varpi = N\delta {\varpi _e}$$, where N is the shape function matrix.

Introducing the boundary condition to the functionality, we can obtain: $$\sum\limits_e {(\delta \varpi )_e^T} \left( {{M_e}{{T'}_e} + {K_e}{T_e} - {Q_e}} \right) = 0$$, where $${M_e} = \int_{\Delta V} {\rho {C_p}{N^T}NdV} ,{K_e} = \int_{\Delta V} {\lambda \nabla {N^T} \cdot \nabla NdV} + \int_{\Delta S} {\eta {N^T}N{\delta _{eb}}dS} ,{Q_e} = {\delta _{eb}}\int_{\Delta S} {\eta {T_b}{N^T}dS}$$, where η and T b are known quantities; for boundary elements, $${\delta _{eb}} = 1$$; for internal elements, $${\delta _{eb}} = 0$$.

The aerodynamic heating on the surface of the window is from the flowfield.

Meanwhile, the heat radiated into the low temperature space from the solid wall is ignored because it has little effect on the overall heat transfer for a relatively small Mach number.

There are two basic conditions for thermal stresses: temperature change and constraint. According to the Newton motion law, the aerospace mechanics balance equation for a 3D model is expressed as follows: ${\kern 0pt} \begin{matrix} \frac{{\partial {\sigma _{xx}}}}{{\partial x}} + \frac{{\partial {\tau _{yx}}}}{{\partial y}} + \frac{{\partial {\tau _{zx}}}}{{\partial z}} + {f_x} = \rho \frac{{\partial {u_x}}}{{\partial t}}, \hfill \\ \frac{{\partial {\tau _{xy}}}}{{\partial x}} + \frac{{\partial {\sigma _{yy}}}}{{\partial y}} + \frac{{\partial {\tau _{zy}}}}{{\partial z}} + {f_y} = \rho \frac{{\partial {u_y}}}{{\partial t}}, \hfill \\ \frac{{\partial {\tau _{xz}}}}{{\partial x}} + \frac{{\partial {\tau _{yz}}}}{{\partial y}} + \frac{{\partial {\sigma _{zz}}}}{{\partial z}} + {f_z} = \rho \frac{{\partial {u_z}}}{{\partial t}}, \hfill \\ \end{matrix}$ where σxx , σyy , and σzz are the normal stresses, τxy , τyx , τxz , τzx , τyz , and τzy are the shearing stresses, fx , fy , and fz are the components of the body force, ux , uy , and uz are the displacements, and $$\rho (\partial {u_x}/\partial t),\rho (\partial {u_y}/\partial t),\rho (\partial {u_z}/\partial t)$$ are dynamic terms.

The physical equation is given by (Ying et al., ): $$\left\{ {\left. {\begin{array}{*{20}{c}} {{\sigma _x}} \\ {{\sigma _y}} \\ {{\sigma _z}} \\ {{\tau _{xy}}} \\ {{\tau _{yz}}} \\ {{\tau _{zx}}} \end{array}} \right\}} \right. = \frac{{E(1 - \mu )}}{{(1 + \mu )(1 - 2\mu )}}\left[ {\begin{array}{*{20}{c}} 1\;{{A_1}}\;{{A_1}}\;0\;0\;0 \\ {{A_1}}\;1\;{{A_1}}\;0\;0\;0 \\ {{A_1}}\;{{A_1}}\;1\;0\;0\;0 \\ 0\;0\;0\;{{A_2}}\;0\;0 \\ 0\;0\;0\;0\;{{A_2}}\;0 \\ 0\;0\;0\;0\;0\;{{A_2}} \end{array}} \right] \cdot \left[ {\begin{array}{*{20}{c}} {{\varepsilon _x} - (1 + \psi ){\alpha _x}\Delta T} \\ {{\varepsilon _y} - (1 + \psi ){\alpha _y}\Delta T} \\ {{\varepsilon _z} - (1 + \psi ){\alpha _z}\Delta T} \\ {{\gamma _{xy}}} \\ {{\gamma _{yz}}} \\ {{\gamma _{zx}}} \end{array}} \right]$$, where ε is the linear strain, γ is the shear strain, A 1=ψ/(1–ψ), A 2=(1–2ψ)/[2(1–ψ)], ψ is the Poisson’s ratio, E is the Young’s modulus, α is the thermal expansion, and ΔT is the increment of temperature.

The strain compatibility equations are expressed as (Ma et al., ) ${\kern 0pt} \begin{matrix} \frac{{{\partial ^2}{\varepsilon _x}}}{{\partial {y^2}}} + \frac{{{\partial ^2}{\varepsilon _y}}}{{\partial {x^2}}} = \frac{{{\partial ^2}{\gamma _{xy}}}}{{\partial x\partial y}}, \hfill \\ \frac{\partial }{{\partial x}}\left( {\frac{{\partial {\gamma _{zx}}}}{{\partial y}} + \frac{{\partial {\gamma _{xy}}}}{{\partial z}} - \frac{{\partial {\gamma _{yz}}}}{{\partial x}}} \right) = 2\frac{{{\partial ^2}{\varepsilon _x}}}{{\partial y\partial z}}, \hfill \\ \frac{{{\partial ^2}{\varepsilon _y}}}{{\partial {z^2}}} + \frac{{{\partial ^2}{\varepsilon _z}}}{{\partial {y^2}}} = \frac{{{\partial ^2}{\gamma _{yz}}}}{{\partial y\partial z}}, \hfill \\ \frac{\partial }{{\partial y}}\left( {\frac{{\partial {\gamma _{xy}}}}{{\partial z}} + \frac{{\partial {\gamma _{yz}}}}{{\partial x}} - \frac{{\partial {\gamma _{zx}}}}{{\partial y}}} \right) = 2\frac{{{\partial ^2}{\varepsilon _y}}}{{\partial z\partial x}}, \hfill \\ \frac{{{\partial ^2}{\varepsilon _z}}}{{\partial {x^2}}} + \frac{{{\partial ^2}{\varepsilon _x}}}{{\partial {z^2}}} = \frac{{{\partial ^2}{\gamma _{zx}}}}{{\partial z\partial x}}, \hfill \\ \frac{\partial }{{\partial z}}\left( {\frac{{\partial {\gamma _{yz}}}}{{\partial x}} + \frac{{\partial {\gamma _{zx}}}}{{\partial y}} - \frac{{\partial {\gamma _{xy}}}}{{\partial z}}} \right) = 2\frac{{{\partial ^2}{\varepsilon _z}}}{{\partial x\partial y}},\hfill \\ \end{matrix}$ where the linear strain and the shear strain are defined as follows: ${\kern 0pt} \begin{matrix} {\varepsilon _x} = \frac{{\partial {u_x}}}{{\partial x}}, {\gamma _{xy}} = \frac{{\partial {u_y}}}{{\partial x}} + \frac{{\partial {u_x}}}{{\partial y}},\; \hfill \\ {\varepsilon _y} = \frac{{\partial {u_y}}}{{\partial y}}, {\gamma _{yz}} = \frac{{\partial {u_z}}}{{\partial y}} + \frac{{\partial {u_y}}}{{\partial z}}, \hfill \\ {\varepsilon _z} = \frac{{\partial {u_z}}}{{\partial z}}, {\gamma _{zx}} = \frac{{\partial {u_x}}}{{\partial z}} + \frac{{\partial {u_z}}}{{\partial x}}. \hfill \\ \end{matrix}$

According to the Von Mises’s criterion, failure occurs when the strain energy of the distortion reaches a critical value. In this case, the equivalent stress σ eq is defined as (Carrera and Giunta, ; Jin et al., ) ${\sigma _{eq}} = \sqrt {\frac{1}{2}\left[ {{{({\sigma _1} - {\sigma _2})}^2} + {{({\sigma _2} - {\sigma _3})}^2} + {{({\sigma _3} - {\sigma _1})}^2}} \right]} \leqslant [\sigma ]$, where [σ] is the permissible stress.

### 2.3.  Physical model

To generalize the study of window structures, a hypersonic vehicle configuration is used in a hypersonic flow, as shown in Fig. 1. The overall length of the configuration is 3 m and it consists of two main sections: the nose and guidance section and the motor section. The nose and guidance section is a wedge-shaped body, such that the window located on the surface is inclined with respect to the freestream flow. The half-angle of the nose is 16°, with the distance from the nose tip being 1.046 m. The simulations are carried out in the Fluent 14.0® environment by a finite volume method.

Fig.1
Hypersonic vehicle configuration

Fig. 2 shows the schematic quarter model used in the simulation with symmetry conditions imposed on the two middle faces. The outermost is a titanium constraint frame to support the window with Grafoil™ at the interface between the frame and the window to provide a seal and minimize stresses. The innermost is the infrared functional window CVD ZnS material. The structural dimensions are shown in Fig. 3. A finite element analysis by ANSYS 14.0® is presented to study the distribution of thermal and stress fields using the coupled cell Solid 5.

Fig.2
Physical model of the window structure

Fig.3
Structural dimensions of window at XY plane (a) and YZ plane (b) (unit: mm)

The thermal and mechanical properties of the window are listed in Table 1. The Grafoil™ and CVD ZnS properties were taken from vendor data.

#### Table 1

Thermal and mechanical properties
 Material Density (kg/m3) Thermal conductivity (W/(m·K)) Specific heat (J/(kg·K)) Thermal expansion (1/K) Young’s modulus (GPa) Poisson’s ratio Mean strength (MPa) Titanium 4510 15.24 520 10.0×10−6 102.01 0.30 650 Grafoil™ 1121 5.19 712 2.7×10−5 0.20 0.30 0.7 CVD ZnS 4090 27.20 515 6.5×10−6 74.50 0.28 93

### 2.4.  Mesh generation

The flowfield model uses quadrilateral elements which are refined in the boundary layer and the number of grids on the X and Y directions are 1406 and 199, as shown in Fig. 4a. An implicit formulation is constructed at all grid points. Gradient is discretized using a Green-Gauss node based scheme. Concerning the flow, the turbulent kinetic energy and a specific dissipation rate are discretized using a second-order upwind scheme. The convergence precisions for the flowfield model are set to 1×10−5.

Fig.4
Discretization of the computational domain for flowfield model (a) and structural model (b)

The structural model consists of thermal and stress fields using hexahedron elements with an initial temperature of 298 K and the number of grids for the structural model on X, Y, and Z directions are 110, 77, and 29, respectively (Fig. 4b). The preconditioned conjugate gradient (PCG) solver is applied for simulating the finite element model, and the convergence precisions for the structural model are set to 1×10−8.

The initial conditions used in fluid models are constant Mach numbers from 3 to 6 flight at an altitude of 15 km in a standard atmosphere. Consequently, the freestream flow conditions are P =12 112.0 Pa and T =216.8 K, and the wall temperature is 298 K.

Aerodynamic heating loads from the fluid models are applied to the skin of the structural models using a convection heating boundary condition, where the convection heat transfer coefficient on the outer surface and the adiabatic wall temperature on the inner surface are the inputs. The window is initially stress-free and has a uniform temperature of 298 K.

## 3.  Results and discussion

### 3.1.1.  Aerodynamic heating model verification

To validate the numerical solver, a case from the experimental investigation by Muylaert et al. () is selected. In the experiment, the surface heat flux referred to as ELECTRE was measured. The freestream flow conditions (293 s) were v =4230.0 m/s, p =53.0 Pa, and T =265.0 K. To test the grid sensitivity, three grids with different resolutions were examined prior to validating their results. The coarse, medium, and fine grids contained 5.160, 1.650, and 0.451 y plus, respectively. The profiles of heat flux variation along the body surface are plotted in Fig. 5. It is shown that the numerical solution based on the fine grid produces heat flux predictions that agree very well with experimental data.

Fig.5
Comparison of heat flux predictions between the propose method and Muylaert et al. ()

### 3.1.2.  Aerodynamic heating studies

Temperature contours of Mach 3, Mach 4, Mach 5, and Mach 6 are shown in Figs. 6a, 6b, 6c and 6d, respectively. Obviously, there is an attached oblique shock at the forebody leading edge at the hypersonic flow. The temperature indicates a large variation in the shock, especially near the solid wall. The maximum heating for all the Mach numbers occurs at the tip of nose (stagnation point) where the flow is usually decelerated and brings it down to zero velocity. In this region, boundary layers are very thin. In addition, the flowfield temperature increases with the increasing Mach number.

Fig.6
Temperature contours of the numerical flowfield at Mach 3 (a), Mach 4 (b), Mach 5 (c), and Mach 6 (d)

The surface heat flux distributions along the body surface at Mach 3, Mach 4, Mach 5, and Mach 6 are shown in Fig. 7. For a constant Mach number, the maximum heat flux occurs near the stagnation zone of the model, and beyond the maximum heating point (stagnation point), the heat flux gradually decreases. After the axial location of 0.2 m (i.e., from the stagnation point), the heat flux is nearly constant, which is clearly seen in Fig. 7. For different Mach numbers, the heat flux increases gradually from Mach 3 to Mach 6. The location of the window is 0.55 m from the nose tip, where it reaches the minimum heat flux value. The heat transfer coefficient and ambient temperature at the window location at different Mach numbers are listed in Table 2.

Fig.7
Heat flux along the body surface at different Mach numbers

#### Table 2

Value of heat transfer coefficient on the point that the window located
 Mach number Heat transfer coefficient (W/(m2·K)) Ambient temperature (K) 3 3084.57 361.78 4 7940.36 458.74 5 14848.22 605.65 6 23440.18 815.18

### 3.2.1.  Aerothermoelastic model verification

To determine the accuracy of the proposed method, the results are compared with those from Russell et al. (). From Fig. 8, we can see that the temperature is in good agreement with Russell et al. (). Therefore, the following work will be done using the proposed method.

Fig.8
Temperature results comparison between the proposed method and Russell et al. ()

### 3.2.2.  Aerothermoelastic studies

To determine the dangerous temperature and thermal stress that is likely to lead to a failure of infrared window, the heating loads from the fluid models (Table 2) are applied to the surface of the window to simulate the temperature and stress response at 10 s.

The thermal and stress responses at Mach 3, Mach 4, Mach 5, and Mach 6 are shown in Figs. 9 and 10, respectively.

Fig. 9 shows that the maximum temperature occurs at the center of the outer surface and the minimum occurs at the surrounding of the inner surface. The temperature decreases from inside to outside and from center to surrounding. It is caused by the severe aerothermal load on the outer surface of the window.

Fig.9
Thermal response at 10 s at Mach 3 (a), Mach 4 (b), Mach 5 (c), and Mach 6 (d)

Fig. 10 shows that the stress at the step corner is high and the stress at the center of the outer surface is low. There also is a stress concentration on the step corner where the dangerous point of the window material lies.

Fig.10
Stress response at 10 s at Mach 3 (a), Mach 4 (b), Mach 5 (c), and Mach 6 (d)

Distributions of maximum temperature, minimum temperature, and maximum stress at different conditions with time are shown in Figs. 11 and 12.

Fig.11
Temperature variation with time

Fig.12
Von Mises stress variation with time

Fig. 11 shows the distributions of maximum and minimum temperatures with time. As illustrated in Fig. 11, the temperature increases with increasing Mach number. At the same Mach number, the maximum temperature which occurs at the outer surface increases sharply within 1 s. After that, the temperature rising rate decreases significantly. Two seconds later, the temperature almost reaches a constant value. Meanwhile, the maximum temperature increases more quickly than the minimum which occurs at the inside surface at first and then the increasing rate stays nearly the same, so the temperature difference between the outside and inside surfaces increases first and then remains nearly the same. In addition, the difference between the outside and inside surfaces at Mach 3 is small, approximately 26 K, but for Mach 6, the difference increases to about 239 K.

It can be seen from Table 3 that the maximum temperatures over the range of Mach 3, 4, 5, 6 are 345.96 K, 446.71 K, 596.34 K, and 806.96 K, respectively, which are all very far from the melting point 2103 K of CVD ZnS, so the material meets the thermal response requirements of the window.

#### Table 3

Thermal and structural responses comparison*
 Mach No. Max temperature (K) Max stress (MPa) Melting point (K) Mean strength (MPa) Result 3 345.96 12.52 2103 93 + 4 446.71 40.72 + 5 596.34 83.94 + 6 806.96 149.21 −

• *+: safety; −: failure

•  Fig. 12 provides the distributions of Von Mises stress with time. Similarly to the temperature, the stresses increase with increasing Mach number. At the same Mach number, the stress increases sharply in an extremely short time and then decreases. The reason for this is that the temperature increasing rate is large at first. About 2.5 s later, the stress increases gradually and finally increases slowly with time because there is a stress concentration on the step corner. At last, the stress almost reaches a constant value.

It can be seen from Table 3 that the maximum stresses at Mach 3, Mach 4, and Mach 5 are 12.52 MPa, 40.72 MPa, 83.94 MPa, respectively, which are less than the mean strength of 93 MPa, so the material is safe in this case. The maximum stress at Mach 6 is 149.21 MPa which is beyond the mean strength of 93 MPa, so the material does not meet the working requirements in this case.

## 4.  Conclusions

After conducting flowfield studies, the temperature and stress distributions for an infrared window are obtained by a finite element analysis method. The judgment of the material’s safety is based on the Von Mises theory of strength.

1. The maximum heating for all the Mach numbers occurs at the tip of the nose (stagnation point). The location of the window is 0.55 m from the nose tip, where the minimum heat flux value is located.

2. The temperature decreases from inside to outside and from center to surrounding. The maximum temperature occurs at the center of the outer surface and the minimum occurs at the edges of the inner surface.

3. The stress decreases from outside to inside and from the edges to the center. The maximum stress occurs at the corner of the step and the minimum stress occurs at the center of the outer surface.

4. The maximum temperatures over the range of Mach 3 to 6 are all very far from the melting point 2103 K of CVD ZnS, so the material meets the requirements of the thermal response of the window.

5. The maximum stresses at Mach 3, Mach 4, and Mach 5 are less than the mean strength of the material, so the material is safe in these cases. However, the maximum stress at Mach 6 is far beyond the mean strength of the material, so it does not meet the working requirements in this case.

Therefore, a suitable material selection is necessary for the design of an infrared functional window for hypersonic vehicles. More work will be done to investigate the temperature and stress distributions of different materials and under different conditions.

* Project supported by the Foundation for Innovative Research Groups of the National Natural Science Foundation of China (No. 51121004), and the Fundamental Research Funds for the Central Universities (No. HIT.BRETIV.201315), China

## References

[1] Carrera, E., Giunta, G., 2008. Hierarchical models for failure analysis of plates bent by distributed and localized transverse loadings. Journal of Zhejiang University-SCIENCE A, 9(5):600-613.

[2] Culler, A.J., McNamara, J.J., 2010. Studies on fluid-thermal-structural coupling for aerothermoelasticity in hypersonic flow. AIAA Journal, 48(8):1721-1738.

[3] Di Clemente, M., Rufolo, G., Battista, F., 2007. An extrapolation from flight methodology for a re-entry vehicle wing leading edge test in a plasma wind tunnel facility. , Proceedings of 39th AIAA Thermophysics Conference, AIAA 2007-3895, Miami, FL, :

[4] Di Clemente, M., Marini, M., Di Benedetto, S., 2009. Numerical prediction of aerothermodynamic effects on a re-entry vehicle body flap configuration. Acta Astronautica, 65(1-2):221-239.

[5] Di Clemente, M., Rufolo, G., Ianiro, A., 2013. Hypersonic test analysis by means of aerothermal coupling methodology and infrared thermography. AIAA Journal, 51(7):1755-1769.

[6] Gerdroodbary, M.B., Hosseinalipour, S.M., 2010. Numerical simulation of hypersonic flow over highly blunted cones with spike. Acta Astronautica, 67(1-2):180-193.

[7] Gnemmi, P., Srulijes, J., Roussel, K., 2003. Flowfield around spike-tipped bodies for high attack angles at Mach 4.5. Journal of Spacecraft and Rockets, 40(5):622-631.

[8] Harris, D.C., 1995. Frontiers in infrared window and dome materials. , Proceedings of SPIE 2552, Infrared Technology XXI, San Diego, CA, 325-335. :325-335.

[9] Harris, D.C., Boronowski, M., Henneman, L., 2008. Thermal, structural, and optical properties of Cleartran multispectral zinc sulfide. Optical Engineering, 47(11):114001-114011.

[10] Heubner, L.D., Mitchell, A.M., Boudreaux, E.J., 1995. Experimental results on the feasibility of an aerospike for hypersonic missiles. , 33rd Aerospace Sciences Meeting and Exhibit, AIAA 95-0737, Reno, NV, :

[11] Jain, A.C., Hayes, J.R., 2004. Hypersonic pressure, skin-friction, and heat transfer distributions on space vehicles: planar bodies. AIAA Journal, 42(10):2060-2069.

[12] Jin, Z.J., Jiang, C.H., Wan, X.P., 2010. Plastic limit load analysis for pressure pipe with incomplete welding defects based on the extended Net Section Collapse Criteria. Journal of Zhejiang University-SCIENCE A (Applied Physics & Engineering), 11(6):440-448.

[13] Ma, T., Chen, Y., Zeng, M., 2012. Stress analysis of internally finned bayonet tube in a high temperature heat exchanger. Applied Thermal Engineering, 43(101-108):101-108.

[14] Marini, M., Di Benedetto, S., Rufolo, G., 2007. Test design methodologies for flight relevant plasma wind tunnel experiments. , Proceedings of West-east High Speed Flow Field Conference, Russian Academy of Science, Moscow, Russia, 1-37. :1-37.

[15] Menter, F.R., 1994. Two-equation eddy-viscosity turbulence models for engineering applications. AIAA Journal, 32(8):1598-1605.

[16] Muylaert, J., Walpot, L., Hauser, J., 1992. Standard model testing in the European high enthalpy facility F4 and extrapolation to flight. , AIAA 17th Aerospace Ground Testing Conference, AIAA 92-3905, Nashville, TN, 1-16. :1-16.

[17] Ruan, J.L., Feng, X., Zhang, G.B., 2010. Dynamic thermoelastic analysis of a slab using finite integral transformation method. AIAA Journal, 48(8):1833-1839.

[18] Russell, G., Stephen, C., Jones, M., 2003. Army hypersonic compact kinetic energy missile laser window design. , Window and Dome Technologies VIII, Proceedings of SPIE 5078, 28-39. :28-39.

[19] Saravanan, S., Jagadeesh, G., Reddy, K.P.J., 2009. Investigation of missile-shaped body with forward-facing cavity at Mach 8. Journal of Spacecraft and Rockets, 46(3):577-591.

[20] Shahmardan, M.M., Norouzi, M., Kayhani, M.H., 2012. An exact analytical solution for convective heat transfer in rectangular ducts. Journal of Zhejiang University-SCIENCE A (Applied Physics & Engineering), 13(10):768-781.

[21] Sun, M., Wu, J.H., 2003. Aerodynamic force generation and power requirements in forward flight in a fruit fly with modeled wing motion. Journal of Experimental Biology, 206(17):3065-3083.

[22] Tang, K., Yu, J., Jin, T., 2013. Influence of compression-expansion effect on oscillating-flow heat transfer in a finned heat exchanger. Journal of Zhejiang University-SCIENCE A (Applied Physics & Engineering), 14(6):427-434.

[23] Tropf, W.J., Thomas, M.E., Klocek, P., 1996. Infrared optical materials. , Proceedings of SPIE CR64, Inorganic Optical Materials, 137-169. :137-169.

[24] White, J.T., 1993. Application of Navier-Stokes flowfield analysis to the aerothermodynamic design of an aerospike-configured missile. , AIAA/AHS/ASEE Aerospace Design Conference, Irvine, CA, :

[25] Ying, J., Lv, C.F., Lim, C.W., 2009. 3D thermoelasticity solutions for functionally graded thick plates. Journal of Zhejiang University-SCIENCE A, 10(3):327-336.