Novel Analytical and Numerical Methods in Heat Transfer Enhancement and Thermal Management

1Dipartimento di Ingegneria Industriale, Universita degli Studi di Napoli Federico II, 80125 Napoli, Italy 2Laboratoire de Modelisation et Simulation Multi Echelle, Equipe Transferts de Chaleur et de Matiere, Universite PARIS-EST, 77454 Marne-la-Vallee Cedex 2, France 3School of Energy and Power Engineering, Xi’an Jiaotong University, Xi’an, Shaanxi 710049, China 4Laboratory of Steam Boilers andThermal Plants, School of Mechanical Engineering, National Technical University of Athens, Zografou, 15780 Athens, Greece 5Mechanical and Aerospace Engineering Department, Rutgers, the State University of New Jersey, Piscataway, NJ 08854-8058, USA

Heat transfer enhancement (HTE) and thermal management (THEMA) are very attractive issues in the research and industry fields. They play an important role in improving energy efficiency and developing high performance thermal systems. HTE and THEMA techniques are powerful tools to increase and improve heat transfer rate and thermal performance as well as to reduce the size of heat transfer systems in installation and operational costs.
The basic purpose of this special issue is to collect original research articles on the most recent analytical and numerical models applied in this field, with the purpose of providing guidelines for future research direction.
This special issue showed novel and interesting analytical and numerical procedures applied in heat transfer enhancement and thermal management.
The topics of the accepted papers cover a wide area of the heat transfer field, from the transient heat transfer to the heat transfer in porous media and to the system design and optimization in forced convection.
A brief overview of each manuscript selected for this special issue is presented in the following, in alphabetical order of the first author.
C. Devaraj et al. investigate numerically natural convection heat transfer in a two-dimensional square enclosure at various angles of inclination using a finite volume based computational procedure. The heat transfer is from a constant temperature heat source of finite length centered at one of the walls to the cold wall on the opposite side while the remaining walls are insulated. The authors analyze the effect of area ratio of the heat source, Rayleigh number, and angle of inclination of the enclosure on the flow field and the heat transfer characteristics. Moreover the paper shows an exhaustive verification and validation section and gives heat transfer correlation equations for each angle of inclination of enclosure investigated which are in good agreement with the numerical results.
I. Simões et al. present, in their paper, a set of fully analytical solutions, together with explicit expressions, in the time and frequency domain for the heat conduction response of homogeneous unbounded and of bounded rectangular spaces (three-, two-, and one-dimensional spaces) subjected to point, line, and plane heat diffusion sources. The authors pay particular attention to the case of spatially sinusoidal, harmonic line sources. This last solution, referred to in the literature as the 2.5D problem. Proposed Green's functions are combined using an image source technique to model a half space, a corner, a layer system, a laterally confined layer system, a solid rectangular column, a solid rectangular column with an end cross section, and a 3D parallelepiped inclusion. This is the first such derivation that is expected 2 Journal of Applied Mathematics to be efficient for formulating 3D thermodynamics problems using boundary elements and integral transforms. The solutions are validated by comparing computed results with those obtained directly in the time domain.
T.-W. Tu and S.-Y. Lee develop for the first time an analytical solution for the heat transfer in hollow cylinders with time-dependent boundary condition and timedependent heat transfer coefficient at different surfaces. The methodology employed by the authors is an extension of the shifting function method and permits transforming the system into a partial differential equation with homogenous boundary conditions. The transformed system is thus solved by series expansion theorem. Limiting cases of the solution are studied and numerical results are compared with those in the literature. The authors investigate the influence of physical parameters on the temperature distribution of a hollow cylinder along the radial direction, too.
X. Zhang and R. Li present two papers. In the first one they investigate the cooling and radiation of a vectoring nozzle and compute the gas spectral characteristic in infrared band, adopting a narrow band model. The radiative heat transfer between the hot gas and the wall is considered with an enclosure model whereas the calculation of film cooling is performed through a cooling effectiveness method. A coupled heat balance equation of heat flux and wall temperature is established on the multilayer structure of the nozzle, including the wall, heat shield, and outer shield, and a Newton-Raphson scheme is taken for solution. The results verification is carried out by means of the investigation of temperature on the expansion part of an experimental nozzle in NASA TN D-1988. The authors also investigate another vectoring nozzle with a multirow of film cooling.
The same authors investigate the similarity of plume radiation between reduced scaling solid rocket models and full scale ones in ground. In the second paper the authors compute flow and radiation of plume from solid rockets with scaling ratio from 0.1 to 1, using CFD code. The radiative transfer equation is solved by the finite volume method in the infrared band 2∼6 m. The spectral characteristics of plume gases have been calculated with the weighted-sum-of-graygas (WSGG) model, and those of the Al 2 O 3 particles have been solved by the Mie scattering model. The research shows a good similarity between spectral variations of plumes from different scaling solid rockets.

Introduction
Problems in thermodynamics can often be solved with the aid of formulas or expressions known as Green's functions. These functions, or fundamental solutions, relate the field variables (heat fluxes and temperatures) at some location in a solid body caused by thermodynamic sources placed elsewhere in the medium.
The fundamental solutions most often used are point sources in a three-dimensional (3D), infinite homogeneous space; line sources acting within two-dimensional (2D) spaces; and plane sources heating one-dimensional (1D) spaces. The reason for these choices is that these three fundamental solutions are known in closed-form in time domain and have a relatively simple structure [1].
They are frequently combined to simulate heat conduction in the time domain or in a transform space defined by the Laplace transform, in half-spaces, infinite plates, rectangular 2D spaces, wedges, and rectangular 3D spaces [1][2][3]. Solutions have also been proposed to deal with multilayer systems, and they include the matrix method [1], the thermal quadrupole method [3], the thin layer method [4], and methods based on the definition of potentials [5][6][7]. Chen et al. have described the use of image method to solve 2D and 3D problems in unbounded and half-space domains containing circular or spherical shaped boundaries [8][9][10].
This paper compiles alternative fundamental solutions in explicit form, specifically Green's functions for harmonic 2D and 3D and harmonic (steady state) line sources whose amplitude varies sinusoidally in the third dimension. This last solution, which is often referred to in the literature as the 2.5D problem, can be of significant value when formulating 3D thermodynamics problems via boundary elements together with integral transforms. In addition, the proposed Green's functions are combined using an image source technique to model a half-space, a corner, a layer system, a laterally confined layer system, a solid rectangular column, a solid rectangular column with an end cross section, and a 3D parallelepiped inclusion. To the best of our knowledge, this is the first such derivation that promises to be efficient for formulating 3D thermodynamics problems using boundary elements and integral transforms.

Fundamental Solution
The transient heat transfer by conduction in an infinite, homogeneous space can be described by the diffusion equation in Cartesian coordinates: in which is time, ( , , , ) is the temperature at a point ( , , ) in the domain, and is the thermal diffusivity defined by /( ), where is the thermal conductivity, is the density, and is the specific heat of medium.
The solution of (1) can be obtained in the frequency domain after the application of a Fourier transform in the time domain, which leads to the following equation: where = √ −1, 1 = √− / and is the frequency. Consider first an infinite, homogeneous space subjected at ( 0 , 0 , 0 ) to a harmonic point heat source of the form where = √( − 0 ) 2 + ( − 0 ) 2 . Consider next an infinite, homogeneous space subjected to a spatially varying line heat source of the form ( − 0 ) ( − 0 ) ( − ) , with being the wavenumber in . This source acts in one of the three coordinate directions, passes through ( 0 , 0 ), and varies sinusoidally in the (i.e., third) dimension. This type of source is often referred to in the literature as a 2.5D source. The response to this source can be obtained by applying a spatial Fourier transform in the direction to the equations for a point heat load.
Applying a Fourier transformation in the direction leads to the solutioñ( where 1 = √− / − 2 , 0 ( ) are Hankel functions of the second kind and order 0. The full 3D solution can then be achieved by applying an inverse Fourier transform in the domain. This inverse Fourier transformation can be expressed as a discrete summation if we assume the existence of virtual sources, equally spaced at along , which enables the solution to be obtained by solving a limited number of 2D problems, ( , , , ) with being the axial wavenumber given by = (2 / ) . The distance chosen must be big enough to prevent spatial contamination from the virtual sources.
Equation (4) can be further manipulated and written as a continuous superposition of heat plane phenomena, ( , , , ) where 1 = √− / − 2 − 2 and Im( 1 ) ≤ 0, and the integration is performed with respect to the horizontal wave number ( ) in the direction. Assuming the existence of an infinite number of virtual sources, we can discretize these continuous integrals. The integral in the above equation can be transformed into a summation if an infinite number of such sources are distributed along the direction, spaced at equal intervals . The above equation can then be written as where 0 = − /(2 ), = − 1 | | , = − ( ) , 1 = √− / − 2 − 2 and Im( 1 ) ≤ 0, and = (2 / ) , which can in turn be approximated by a finite sum of equations ( ). Note that = 0 is the 2D case,̃( , , ) = 0 ∑ =+∞ =−∞ ( / 1 ) with 1 = √− / − 2 . Next, the above Green's functions are combined so as to define Green's functions for a half-space, a corner, a single layer system, a U system, a solid rectangular pipe, a solid open box, and a 3D parallelepiped box. Expressions in frequency and time solutions are provided. The time solutions obtained after the application of inverse spatial and frequency Fourier transforms are compared with those given by Green's functions defined directly in the time domain.
Green's functions for the different spaces are determined using the image source method. By this method a distribution of virtual sources and sinks are combined so as to give null temperatures (Dirichlet boundary conditions) or heat fluxes on the required boundaries (Neumann boundary conditions). Other boundary conditions, such as Robin, are not studied in this paper. In the case of solid bodies bounded by two parallel surfaces the number of sources, placed perpendicular to the surfaces, is theoretically infinite. Table 1 = √ ( ) 2 + ( ) 2 + ( ) 2 0 = − 0 0 = − 0 0 = − 0 1 = + 0 − 2 1 1 = + 0 − 2 2 1 = + 0 − 2 3 2 = + 0 + 2 1 ( − 1) 2 = + 0 + 2 2 ( − 1) The superscripts , , and identify the position of the virtual sources along the , , and directions, respectively. The upper value of , , and is defined by the convergence criteria. Each value of , , and is associated with four possible source positions, which are identified by the subscripts , , and for the , , and directions, respectively. Thus, , , and may take the values of 1, 2, 3, and 4.
The use of complex frequencies allows the contribution of the sources placed at greater distances to vanish and so to limit the number of the virtual sources. The use of complex frequencies with a small imaginary part, taking the form = − (where = 0.7Δ and Δ is the frequency increment), has the additional effect of avoiding the aliasing phenomena. This shift in the frequency domain is subsequently taken into account in the time domain by means of an exponential window, , applied to the response.
Green's functions are validated assuming that the medium is subject to a Dirac delta source. This type of source would require the solution to be computed in the frequency domain [0.0, ∞] Hz. However, the response does not need to be computed for a very large number of frequencies since it decays very quickly as the frequency decays. Note that the static response for the frequency 0.0 Hz can be calculated thanks to the use of complex frequencies.
The number of virtual sources used depends directly on the predefined convergence criterion. As we move from one dimension to two dimensions and then to three dimensions, the number of sources grows significantly. Thus, although the method converges rapidly, the cost of computation grows significantly as we move from a one-dimensional to a threedimensional problem.

Green's Functions
Green's functions in the time and frequency domain will be grouped for the following three cases: (i) unbounded space, which includes Green's functions for 1D, 2D, and 3D sources; (ii) two-dimensional space, which contains Green's functions for a half-space, a space bounded by two perpendicular planes, a single layer system, a U system, and a solid rectangular pipe, when subjected to 2D and 3D sources; ( All results showed a good agreement between the different formulations for all cases.
In the examples provided the Dirac delta source is positioned at the coordinate ( 0 = 0. The notation in Table 1 is used.
Equations for an Unbounded Space, Assuming 1D, 2D, and 3D Heat Sources. 1D source is as follows:   2D source is as follows: , as the sum of plane sources.

Two-Dimensional Space
(a) Half-Space Defined by ≥ 0. Boundary conditions prescribed for the half-space (Cases 1 and 2) are shown in Figure 10.
Equations for the half-space (Cases 1 and 2), subjected to a 2D heat source, are as follows.
Case 2. Consider Equations for the half-space (Cases 1 and 2), subjected to a 3D heat source, are as follows.
as the sum of 2.5D sources.
as the sum of 2.5D sources.
(b) Bounded Space Defined by ≤ 1 and ≥ 0. Boundary conditions prescribed for the bounded space defined by ≤ 1 and ≥ 0 (Cases 1-4) are shown in Figure 11. Equations for the bounded space defined by ≤ 1 and ≥ 0 (Cases 1-4), subjected to a 2D heat source, are as follows.

Journal of Applied Mathematics
Equations defined for the bounded space defined by ≤ 1 and ≥ 0 (Cases 1-4) and subjected to a 3D heat source, are as follows.
as the sum of 2.5D sources.
as the sum of 2.5D sources.
as the sum of 2.5D sources.
as the sum of 2.5D sources.
(c) Horizontal Layer Defined by = 0 and = 2 . Boundary conditions prescribed for the horizontal layer (Cases 1-4) are shown in Figure 12.
Equations for the horizontal layer (Cases 1-4), subjected to a 2D heat source, are as follows.
as the sum of 2.5D sources. (27) as the sum of 2.5D sources. (28) as the sum of 2.5D sources. (29) as the sum of 2.5D sources. (30) (d) U System Bounded by ≥ 0, ≤ 2 , and ≤ 1 . Boundary conditions prescribed for the U system (Cases 1-4) are shown in Figure 13. Equations for the U system (Cases 1-4), subjected to a 2D heat source, are as follows. . (31) . (32) Journal of Applied Mathematics 13 Case 4. Consider Equations for the U system (Cases 1-4), subjected to a 3D heat source, are as follows.
as the sum of 2.5D sources.
as the sum of 2.5D sources. (36) as the sum of 2.5D sources.
as the sum of 2.5D sources.
Equations for the solid rectangular pipe (Cases 1-4), subjected to a 2D heat source, are as follows.
as the sum of 2.5D sources. (43)

16
Journal of Applied Mathematics as the sum of 2.5D sources. (44) as the sum of 2.5D sources. (45) as the sum of 2.5D sources. (46)

Three-Dimensional Space (a) Solid Open Box Defined by
as the sum of 2.5D sources. (47) as the sum of 2.5D sources. (48)
Equations for the 3D parallelepiped box (Cases 1-4), subjected to a 3D heat source, are as follows.
as the sum of 2.5D sources. (51)

20
Journal of Applied Mathematics as the sum of 2.5D sources. (52) as the sum of 2.5D sources. (53) Journal of Applied Mathematics as the sum of 2.5D sources. (54)

Conclusions
Fully analytical solutions for heat conduction for unbounded and rectangular spaces subjected to point, line, and plane sources have been presented. Two boundary conditions were assumed, namely, the Dirichlet and Neumann boundary conditions. Particular attention was given to the two-anda-half-dimensional fundamental solution or 2.5D Green's functions defined for spatially sinusoidal, harmonic line sources. The final expressions were validated by applying the equations to the problem of a Dirac delta source, for which the solutions in the time domain are known in analytical form. Excellent agreement was found between the numerical solutions given by Fourier synthesis and the exact solutions.

Introduction
As for base heating and thermopneumatic property of the solid rocket, plume radiation has been the focus of the investigation in thermal protection design of solid rockets in past decades [1]. Much related experimental and theoretical study has been carried out to study radiation characteristics of the solid rocket plume. As a limitation of ground test, most experiment study of rocket plume radiation used reduced scaling model, and also some of the theoretical studies used scaling model. Considering the application of research results with reduced scaling models in radiation knowledge and thermal protecting design of full scale rockets, similarity between plume radiation of reduced scaling model and full scale rocket must be obtained. However, to the authors' best knowledge, there has not been experimental study about that problem published till now. Only several numerical analysis studies can be found out. Rozanov and Lyapustin had studied a new integral form of similarity conditions to the error analysis of truncation techniques for the forward peak scattering and deduced the integral similarity method, that is, Delta-M [2]. Jun and Wenbing had studied the self-similarity of transmittance depth and radiation energy in a cool medium, with changing time, boundary temperature, and medium density [3]. Research of Duracz and Mccormick was focused on the ratio of two similarity parameters, radiation intensity and irradiance, and the influence of the single-scattering albedo and asymmetric coefficients within weakly absorbing clouds [4]. Mitrescu and Stephens had studied a new approach for determining the scaling parameter of the radiative transfer equation and proposed that the key assumption regarding the angular dependence of the radiative field is essential [5]. Bril et al. had studied the similarities between temperature field and concentration field in the near nozzle area; this study also proposed a method for determining a nondimensional radiance intensity with the outlet parameter, radiation wavelength, and temperature gradient of absorptivity [6].
Calculation of solid rocket plume radiation is highly complicated as the RTE is a high ordered nonlinear equation with differential term and integral terms, considering the strong spectral properties of plume gases and radiative absorption and scattering of groups of Al 2 O 3 particles with different temperatures and diameters. The plume gases consist of several polarity gases: H 2 O, CO 2 , CO, HCL, and OH, whose infrared radiative spectrum is made up with ten thousands of lines [7]. The common detailed calculating method for gas radiation uses a band model [8,9] or the weighted-sum-ofgray-gas (WSGG) model [10]. The band model has a higher precision for nonhomogeneous gases but is somewhat time lasting for large volume gases and is not very suitable for solid rocket plume gases. The WSGG model always shows a compromise of precision and calculation time, which has been widely used to calculate radiation of clouds or combustion gases. As for spectral property of Al 2 O 3 particles, effects of particle temperature, concentrations, and diameters are needed to be considered [11,12], and the Mie scattering theory for calculating particle radiation can also be used for Al 2 O 3 particles [13]. Many numerical solution methods have been developed to calculate radiation in nonhomogenous absorptive/emissive/scattering medium by solving RTE numerically, including Monte Carlo (MC) methods [14], the streaming method (SM) [15], discrete-ordinate method (DOM) [16], and finite volume methods (FVM) [17]. To improve the solution accuracy with less execution time, the FVM has often been used to calculate the spectral radiation of the high temperature two-phase plume with highly expanding volume.
To study the similarity of infrared radiation from plumes with different scaling ratio, the Trident D5 solid rocket and a series of reduced scaling models with similar geometry have been taken for investigation. CFD code has been used to compute the axial symmetric flowing field inside the rocket engine and in plume of the full scaled rocket and reduced scaling model, with same combustion chamber pressure and plume's flowing Mach number. The spectral absorption coefficient of plume gases: H 2 O, CO 2 , CO, HCL, and OH, will be computed with WSGG model. Eight groups of Al 2 O 3 particles with different diameters between 2 and 16 m have been taken into consideration. The Mie scattering model has been used to compute the spectral absorption, scattering coefficient, and phase function of Al 2 O 3 particles. To calculate the spectral radiance of the plume, a finite volume method has been used to solve the RTE. Integrating spectral radiance on exterior surface of plume volume in one direction will give the plume radiation intensity. Ratio of radiation intensity between reduced scaling model and full scaled one has been computed to study the similarity rules of plume radiation.

Calculation of Plume Flowing Field
The Trident D5 solid rocket used a composite propellant NEPE; there is 10% weight of Al in the composite propellant. Plume's flowing data of the full scaled rocket and reduced scaling models are calculated with the same input data for the rocket nozzle, pressure is 9 Mpa, and temperature is 3750 K. Input concentrations of major species are OH = 0.192, H 2 = 0.299, CO = 0.027, CO 2 = 0.146, and H 2 O = 0.0008, on the inlet of nozzle. The full scaled rocket has a nozzle with 1.55 m length, 0.35 m diameter, and 9.7 ratio of area expansion. The reduced scaling model has the same geometry but with a minimized diameter and length. To obtain the flowing data of plume, the CEA program has been adopted to compute the chemical-equilibrium compositions of the propellants and the inlet parameters.
Then, CFD code is used to compute the united flowing in rocket nozzle and plume. The time-marching method and the advection upstream splitting method (AUSM) have been used for numerical discretization of the 2D axis symmetric N-S equations. The -turbulence model is adopted to simulate blending of the plume and atmosphere. A finite rate chemical reaction model with 12 components and 17 reactions has been taken for consideration, the detailed reaction equations and related coefficients can be seen in [18]. For the Al 2 O 3 particles, the Lagrangian particle trajectory model has been used to simulate exchange of energy and momentum between particles and plume gases. The distribution of the particle diameter is determined by Braithwaite's proposed functions [19].
The plume's flowing data of the full scaled Trident D5 rocket and the 9 reduced scaling models have been calculated in our study. Only two groups of result are given in Figure 1. One is for the full scaled rocket and the other is for the 0.5 scaling model. Considering the axial symmetric characteristics of the plume's flowing field, only flowing field map on the upper -plane has been plotted. Figure 1 shows pressure and temperature of plume gases and volume fractions of the 3 major components: H 2 O, CO 2 , and CO. Volume fractions of HCL and OH molecule have not been shown here as their fractions are very small, less than 0.05. Temperature and number concentration of the 8 groups of Al 2 O 3 particles have been calculated in our study, only results of 3 groups with diameters being 6 m, 8 m, and 10 m are given here.
One can see from the figure that volume of the plume expands rapidly after the rocket nozzle, making the diameter of the high temperature core at = 20 m about 5 times that of the nozzle exit. All maps of the flowing map of the 0.5 scaling rocket model have close similarity with the full scaled rocket, with similar parameter distributions and a close data range. One can also find that the high temperature areas of gases and Al 2 O 3 particles are a continuous central stripe in the plume, but the dense concentration area of Al 2 O 3 particles number is quite discrete. Comparing volume fraction of CO 2 and CO, it can be found that the high volume fraction area of CO 2 is in the radial medium ring, at some axial distance after the nozzle exit, while that of CO is in the central thin stripe close behind the nozzle exit. For H 2 O, the high volume fraction area is continuous in the central stripe of the plume.

Calculation Method of Plume Radiation
A widely used equation of radiation transport in a nonhomogeneous absorptive/emissive/scattering medium can be written as variation of radiance ( ⃗ S) while passing a distance S along the path S in the medium [20,21] i (S)    where the subscript " " means spectral variable and the superscript " " means directional variable. The symbol Σ is the summed absorption coefficient of plume gases and Al 2 O 3 particles, and is the summed scattering coefficient of 8 groups of Al 2 O 3 particles, which are simply summed up mathematically. The symbol S means direction of the incident radiation, means the solid angle in S direction, and Φ(S, S ) means the scattering phase function between S and S . The above symbols: , , and Φ are all spectral variables, while the notations " " are omitted in (1) for abbreviation.
According to finite volume method, (1) can be solved to compute spectral radiance of plume by integrating the equation on a control volume. Considering the plume geometry, the cylindrical control volumes had been firstly used for solving (1), but it is found that numerical convergence of the developed coefficient matrix is very poor, being as the area difference between the two radial opposite faces of the control volume. For that reason, the cuboid control volume has finally been chosen in our study. The new computation domain for radiation computation has the same length with the flowing domain in Figure 1, but it has a square form end with side length being equal to the outer radius of the plume. The precision of plume radiance solution is mainly affected by flowing data: temperature, volume fraction, or number concentration of plume gases and Al 2 O 3 particles. Those flowing data do not change dramatically in the plume as seen in Figure 1, so the mesh size for radiance solution can be much larger than that used in CFD solution. For plume of the full scaled solid rocket, the meshing plot on lengthwise section of the computing domain is shown in Figure 2. The divided grids number is 80 × 20 in the upper -plane. Flowing data of control volume inside the plume takes the arithmetic average of data on CFD grids in the volume. That of control volume outside the plume is set to zero. In that case, all control volumes outside the plume would not actually influence the result of plume radiance, but the integral RTE equation is uniform for all the control volumes, and solution convergence can be improved.
Spectral radiance of each control volume is considered in a series of directions; every direction is corresponding to a set of angle, zenith angle, and azimuth angle. The zenith angle is defined as the angle between the plume central line ( -axis) and the direction vector, and azimuth angle is defined as angle between the projected direction vectors on -coordinate plane and the axis. For radiance on control volume in direction S , integrating (1) with Gaussian integral method, radiance of can be related to radiance of its neighbored control volumes in that direction. The detailed derivation procedure can be found in [19]; only the final integration equation is given here as follows: where subscripts , , , , , and mean the variables of the neighbored control volume and , , , , , and mean the variables on interface of neighboring volumes, as shown in Figure 3. The symbol Δ is the area of the interface, n is external normal vector of interface, and Δ is the volume of control volume. The upper equations can be also written as the following simplified form: The compound integral RTE equations of all control volumes have the form of vector equations Equation (4) is a nonlinear system of equations, as variable b in the right side is actually unknown and needs to be determined by radiance of control volume in all the other directions, S . So, the solution of (4) could not be completed with any iteration algorithm in one loop. A cyclic iteration algorithm with modification of b has been proposed in this study. In each cycle, b is firstly calculated with the radiance computed from the last cycle. Then, a Gauss iteration algorithm will be used to solve vector equations (4) to compute a new group of radiance for each control volume. The cyclic iteration will go on till the following convergence limitation is met: Considering the symmetric characteristics of geometry and flowing data of the plume, radiance of plume is also symmetric. So, we do not need to solve the vector equations of RTE on all control volumes in the plume, but only those control volumes with centers on -coordinate plane have been chosen for calculation. While radiance of neighbored control volumes is also required in (4), those control volumes are not centering on the -coordinate plane, and a radiance transformation algorithm on symmetric control volumes has been proposed in our study to calculate the radiance of other control volumes in the plume. The radiance transformation algorithm is derived based on rotation rules of vectors seen in Figure 4. If the radiance of control volume in direction S 1 is i (S 1 ), the radiance of control volume in direction S 2 is equal to i (S 1 ) when the following vector equation is satisfied:

Calculating Radiation Intensity from the Plume Surface
In the cases of studying base heating or heat protection of the solid rocket, spectral radiation intensity on the whole surface of the plume is concerned. Spectral radiation intensity is also a directional variable, which is the integration of spectral radiance i (S), on surface of bounding control volumes in direction S. If radiation intensity in some waveband is concerned, as for the 2∼6 m band, integration on the wavelength is also needed where is number of control volume surfaces on boundary of the plume, e S is unit vector of the radiance direction, and n Δ is unit vector of external normal on the surface.

Calculation of Spectral Property of Plume Gases and Al 2 O 3 Particles
The 3 major radiation gases of the solid rocket plume are H 2 O, CO 2 , and CO. To calculate the spectral absorption coefficient of each plume gas, the standard absorption coefficient ,STP at several temperatures and 1 atm pressure in [20,21] has been used. For plume gas with pressure being and temperature being , absorption coefficient can be modified from ,STP as follows: Absorption coefficient of the plume gases mixture is calculated with the WSGG model, and the weight is the volume fraction of each gas. If the volume fraction of one gas is , then the absorption coefficient of the gases mixture can be written with the WSGG model as The calculated absorption coefficient for plume gases H 2 O, CO 2 , and CO and the gases mixture in center of front end of the plume have been shown in Figure 5. One can find that absorption coefficient of CO 2 and CO is comparatively much larger than that of CO is 4.2∼5.5 m, the peak spectrum of CO 2 is 4.2∼5 m, and absorption coefficient of H 2 O is much smaller. Absorption coefficient of the gases mixture appears as an accumulated effect of the three gases.

Calculation of the Spectral Property of Al 2 O 3 Particles
In the calculation of spectral properties of Al 2 O 3 particles with Mie scattering model, particle size, number concentration, and complex refractive index must be determined firstly. As mentioned above, size and number concentration of For a single Al 2 O 3 particle, spectral properties including the scattering cross section , attenuation cross section , scattering factor , attenuation factor , and the scattering phase function Φ can be computed with Mie scattering model [13] as where and are Mie scattering coefficients and 1 and 2 are amplitude functions. Computation of the above  coefficients and functions has been illustrated in detail in [13], which is not repeated in this work. As for the spectral properties of the particle clouds, the single and independent effect of each particle is assumed and all the spectral properties of the particle clouds can be computed as the mathematic sum of all particles.
The calculated absorption coefficients of 6 groups of Al 2 O 3 particles with different diameters and the total value of particles cloud in the center of plume's front end have been shown in Figure 6. That of the other two groups with = 4 m, 12 m is not included in the figure as their concentration is zero at front end of the plume. One can see that absorption coefficient of Al 2 O 3 particles does not change significantly with wavelength like that of plume gases. The three groups of Al 2 O 3 particles with smaller diameters ( = 4 m, 6 m, 8 m) have larger absorption coefficient. Compared with Figure 5, it is found that the absorption coefficient of particles cloud is much larger than that of plume gases, so Al 2 O 3 particles are the major radiation composition in the rocket plume.
The calculated scattering coefficients of the 6 groups of Al 2 O 3 particles with different diameters and total value of particles cloud in the center of plume's front face have been shown in Figure 7. One can find that scattering coefficient for Al 2 O 3 particles with = 4 m is most strong, which is medium for particles with = 6 m and is very small for extra particles. Compared with Figure 6, it is found that scattering coefficient of particles cloud is much bigger than the absorption coefficient, and the scattering coefficient will grow up with the increasing wavelength.
The calculated scattering phase functions of Al 2 O 3 particles cloud for different wavelengths in the center of plume's front end have been shown in Figure 8. The scattering angle in the figure means the angle between the incident radiance and scattering direction. It can be seen that the scattering phase function in forward directions is much larger than

Results and Discussion
The calculation method of spectral property of plume gases and Al 2 O 3 particles and infrared radiation intensity has been verified with a reduced scaling solid rocket in a vacuum chamber. Infrared radiation intensity in wavebands 2.7∼ 2.95 m and 4.2∼4.45 m at three zenith angles = 60 ∘ , 90 ∘ , and 120 ∘ have been calculated and compared with test results. The maximum difference between calculated and test results is 20.9%, which arrives at = 90 ∘ , and the minimum difference is 10.4% which arrives at = 60 ∘ [22]. Studying the similarities of plume's radiation intensity from solid rocket with different scale and ratio of plume radiation between the reduced scaling and full size rockets in wavelengths = 2 m, 3 m, 4 m, and 5 m, and waveband 2∼6 m have been calculated and illustrated in Figure 9. The horizontal ordinate in the figure means the scaling ratio of rocket model, and the vertical ordinate / 0 means the ratio of plume's radiation intensity. The three dashed lines are corresponding to the 1.5, 2, and 2.5 power function of scaling ratio. Results of = 0 ∘ , 40 ∘ , 80 ∘ , 120 ∘ , and 160 ∘ , with = 30 ∘ , have been investigated. It can be seen in the figure that ratio of radiation intensity increases with the increasing scaling ratio, for all the four wavelengths and waveband 2∼6 m. The increasing rule of radiation intensity ratio for = 0 ∘ gets close to the curve of 1.5 power mostly, which gets closer to the 2power function of scaling ratio for = 40 ∘ and comes up to reach 2.5 power of scaling ratio for = 120 ∘ and 160 ∘ . This is caused by changing the length of high radiance section and  radiance attenuation of the cold gases in tail section of the plume. As plume's radiance mainly comes from the central stripe which has a higher temperature and concentration for both gases and Al 2 O 3 particles. Radiance attenuation of the tail cold gases increases more rapidly than the projected area on the surface of the plume with the growth of the scaling ratio, making growing rate of / 0 being smaller than that of surface area; the latter has a growing rate of 2 in small zenith angle. With the increase of zenith angle, the effect of radiance attenuation decreases greatly and causes the ratio of radiation intensity increasing to reach and surpass 2 .
In the viewing of spectral characteristics for solid rocket plume with different scaling ratios, the variation of radiation intensity with wavelength for the 9 reduced scaling rocket models and the full scaled one in two directions, = 30 ∘ and 90 ∘ with = 30 ∘ , has been plotted in Figure 10. One can see good similarity between spectral variations of different scaling rocket models. There are two peak spectra of radiation intensity, 2.7∼3.0 m and 4.3∼4.6 m. Compared with Figure 5, it is found that the former spectrum is the feature radiative spectrum of H 2 O; the latter is of CO 2 and CO. Besides, in the two peak spectra, spectral radiation intensity for plume of the 10 rockets decreases with the increasing wavelength in 2∼6 m waveband. Since the Al 2 O 3 particle is the major radiation composition in the plume, which has the maximum radiation for < 2.0 m as particle, the temperature is higher than 2500 K.
The spectral radiation intensity of plume in all the zenith angles with = 30 ∘ , for the 0.5 scaling rocket model and the full scaled one, has been shown in Figure 11; results of wavelength = 2 m, 3 m, and 4 m have been plotted. With the increase of the zenith angle, radiation intensity of the plume firstly increases, arrives at a maximum value in a middle angle, and decreases after that angle, since the projected area on plume surface changes with the zenith angle. Comparing radiation intensity of plume in the front and rear hemisphere, it's found that is stronger in the rear hemisphere. As read end surface of plume has a strong radiance in the rear hemisphere. Effect of radiance from the rear end also makes the direction with maximum radiation intensity deflects toward the rear hemisphere, which is remarkable for = 0.5, and the maximum radiation intensity arrives at = 50 ∘ in that case. The colder section in the tail of plume is shorter and the radiance from the rear end surface would be bigger for a smaller scaling rocket. Figure 12 shows map of spectral radiance on -coordinate plane in plume flow of the 0.5 scaling solid rocket and the full scaled ones, for wavelength = 3 m and 5 m, in = 0 ∘ , 45 ∘ , and 90 ∘ . The radiance map of the two solid rockets has similar shape and close magnitude. It is also found that the high light area with a strong radiance is located in the central stripe of the plume, where temperature and concentration of gases and Al 2 O 3 particles are very high. Compared the width of the high light stripe in radiance maps of the two wavelengths, it is narrower for = 3 m than = 5 m. As the radiation of Al 2 O 3 particles plays the main role in plume radiance of 3 m, the area with dense concentration of particles is smaller in radius than plume gases, causing high radiance area centralizing in the plume in that wavelength. Plume gases are the major radiation composition for = 5 m, which has a bigger expansion angle and radius size, and cause a wider stripe with high radiance. Viewing the axial length of the high light stripe, it is in the front part with half length of the plume for the 0.5 scaling model, but extends from the front to end part of plume for the full scaled rocket model.

Conclusion
Infrared radiation of plume from one full scaled solid rocket as well as the 0.1∼0.9 reducing scaled models in ground condition has been investigated to study the similarity of plume's infrared radiation for rocket models with different scaling ratios. Flowing field in the rocket and plume has been computed with CFD code, the omnidirectional and infrared radiance on symmetric plane in the plume has been calculated with the developed FVM code. The research shows that Al 2 O 3 particle is the major radiation composition in the rocket plume, whose scattering coefficient is much larger than its absorption coefficient. Ratio of radiation intensity from reduced scaling plume to that of the full scaled one increases with the scaling ratio. The power of the growth curve is 1.5 for azimuth angle = 0 ∘ , becomes 2 for = 40 ∘ , and arrives at 2.5 for = 120 ∘ and 160 ∘ . There is good similarity between spectral variations of plumes from different scaling solid rockets, and there are two peak spectra of radiation intensity in 2.7∼3.0 m and 4.3∼4.6 m. Radiation intensity of the plume in rear hemisphere is bigger than that in the front, and the effect of radiance from the rear end surface makes the direction with maximum radiation intensity deflects toward the rear hemisphere.

Introduction
Natural convection heat transfer in enclosures has been extensively studied by researchers, because of its practical significance in science and technology. Applications include heating and cooling of building spaces, solar energy collectors, heat exchangers, and effective cooling of electronic components and machinery. The fluid flow and heat transfer behaviour of such systems are analysed numerically and experimentally by a number of researchers, incorporating different geometrical and thermal boundary conditions, from low to high Rayleigh numbers. The studies conducted by different groups were extensively summarized by Wan et al. [1] and a set of benchmark quality data was also presented in this study for the natural convection range of 10 3 ≤ Ra ≤ 10 8 . On the other hand, several investigations [2][3][4] were carried out by the researchers on high Rayleigh numbers above 10 8 where flow becomes turbulent for which numerical solutions are obtained by applying suitable turbulence models. Mao and Zhang [5] have carried out numerical simulation of turbulent natural convection of compressible air in tall cavities to evaluate the accuracy of turbulent models for various temperature boundary conditions. This study reveals that -turbulence model is good in predicting velocity distribution while LES model is good in predicting temperature distribution at high Rayleigh numbers. Unsteady natural convection in a rectangular cavity, numerically investigated by Patterson and Imberger [6], shows that the transient flows strongly depend on the Prandtl number Pr and ultimately converge to identical steady state solutions. The numerical analysis conducted on convection in a cavity with variable viscosity fluid by Kandaswamy et al. [7] confirms that the heat transfer rate increases with the increase in viscosity of the fluid. Analytical studies [8,9] on natural convection in a cavity with volumetric heat generation exemplify the convergence behaviour of various computational methods. A finite method approach by Saha et al. [10] for natural convection in a square tilt open cavity reveals the strong influence of the inclination angle of the cavity on convection heat transfer in addition to the influence of thermal boundary 2 Journal of Applied Mathematics conditions. Corcione [11] has applied SIMPLER algorithm to determine the effect of thermal boundary conditions at the sidewalls on natural convection in a rectangular enclosure heated from below and cooled from above for the range of Rayleigh numbers between 10 3 and 10 6 . A study on laminar natural convection in shallow air filled rectangular cavities was presented by Samy and Elsherbiny [12] which showed strong dependency of natural convection heat transfer on Rayleigh number, enclosure aspect ratio, and angle of inclination of the enclosure. Aminossadati and Ghasemi [13] investigated the effects of orientation of an inclined enclosure with two adjacent sidewalls at different temperatures on natural convection heat transfer and observed that at high Rayleigh numbers the flow pattern and the rate of heat transfer depend on the angle of inclination of the enclosure. Kandaswamy and Eswaramurthi [14,15] have focused their study on the determination of the effect of density maximum of water around 4 ∘ C on natural convection in a square porous enclosure applying finite volume method. The influences of porosity, Grashof number, Darcy number, and internal heat generation on natural convection heat transfer were presented in detail in these studies. Numerical and experimental studies were conducted by different researchers by introducing geometrical features such as baffles, partition walls, fins, and conducting blocks proposing various heat transfer applications [16][17][18][19]. In each study, the configuration is described to affect the flow pattern within the enclosure which in turn influences the natural convection heat transfer. A number of studies on square enclosures with time dependent boundary conditions were found in the literature. For illustration, Aswatha et al. [20] applied finite volume method to study the effect of different thermal boundary conditions at bottom wall such as uniform/sinusoidal/linearly varying temperatures on natural convection in enclosures. The results presented from the abovementioned study show that uniform temperature gives higher Nusselt number compared to the sinusoidal and linearly varying bottom wall temperature cases.
In many practical situations, there is a need to view the natural convection in inclined enclosures where the heating component which can be treated as heat source does not spread over the entire length of the inner wall. The directions of gravitational and buoyancy forces remain the same, irrespective of the changes in the angle of inclination of the enclosure which alters the position of the heat source. The combined effect of size of the heat source and the angle of inclination causes additional changes in the flow pattern and the heat transfer characteristics. In the present study, natural convection phenomena in the enclosure for steady laminar flow are investigated for different values of area ratios of the heat source in the range of 0.2 to 1.0 and the angle of inclination varied from 0 ∘ to 360 ∘ . Each configuration is analysed for various laminar ranges of Rayleigh numbers, 10 3 ≤ Ra ≤ 10 7 .

Mathematical Formulation
The schematic diagram of the two dimensional square enclosure of size × , consisting of a discrete heat source of length at a constant temperature ℎ at one of its inner walls, is depicted in Figure 1. The wall opposite to the heat source acts as heat sink at constant temperature . Except the heat source and the sink, the remaining parts of the enclosure walls are assumed to be adiabatic. The origin of the coordinate system is fixed at one corner of the wall containing heat source. Numerical analyses are carried out for different area ratio of heat source / equal to 0.2, 0.4, 0.6, 0.8, and 1.0 assuming that the enclosure contains fluid which is incompressible and laminar. Except the density which is approximated by the Boussinesq model, the fluid properties are assumed to be constant as elaborated by Gray and Giorgini [21]. Heat is transferred from the heat source to the cold wall. Thermally induced density gradient develops a buoyant flow and the gravitational force always acts on the fluid particles in the vertically downward direction irrespective of the angle of inclination of the enclosure.

Governing Equations.
Conservation of mass, momentum, and energy equations are the governing equations to be solved for natural convection flow: Energy: The properties of the fluid and the acceleration due to gravity are assigned with the appropriate values to obtain the desired Rayleigh number which is given by The average Nusselt number over the cold wall, which is a measure of convection heat transfer in the enclosure, is estimated by using the relation

Computational Method and Grid Independence Study
In the present study, the finite volume based computational procedure is applied to solve the governing equations. The SIMPLE algorithm suggested by Patankar [22] is used for pressure velocity coupling and the second order upwind scheme is used for solving momentum and energy equations. Velocity gradients of the fluid would be high near the inner wall due to the boundary layer formation in comparison with the far field region. Rectangular form grid structure is created as illustrated in Figure 2 with fine grids close to the inner wall to capture the boundary layer accurately and with coarse grids near the centre of the enclosure to reduce the computational time.  Figure 3. Grid comparison indicates that grid independent results are obtained for grid sizes 40 × 40 and above. In order to ensure better accuracy of results, grid structure 60×60 is selected as the ideal one for all the analyses in the present study. The convergence criteria were set to 10 −6 for all the relative (scaled) residuals.

Validation Study
The methodology applied in the present study is validated rigorously by conducting similar studies reported by Wan et al. [1] as benchmark data for a wide range of Rayleigh numbers. Contour plots and -plots of various parameters are compared and found to produce similar results between the two studies. The comparison of nondimensional horizontal velocities ( ) at the midwidth ( = 0.5) is presented in Figure 4 to demonstrate the attainment of validation study. Pointwise measurements for the Rayleigh numbers 10 3 to 10 6 are perceived to be in good agreement between the two results. The lack of comparison in the pointwise measurement begins to appear for Rayleigh number equal to 10 7 at few corresponding points of the cavity. For example there is a maximum difference of 10% for Rayleigh number 10 7 at 0.1 unit distances from the top and bottom walls of the cavity due to high velocity gradients at these levels. However, the present study is focused on the laminar range of Rayleigh numbers not exceeding 10 7 .

Results and Discussion
The results of numerical investigation focused on the influence of area ratio of the heat source and the angle of inclination of the enclosure for Rayleigh numbers 10 3 to 10 7 are presented in this section. Contour plots are used for qualitative description considering typical cavities with heat source area ratios 0.4 and 0.8 and angles of inclination 0 ∘ and 60 ∘ . The influences on performance parameters are described quantitatively using -plots.

Influence of Area Ratio and Angle of Inclination on Flow
Pattern. Figure 5 shows the results associated with flow field in terms of streamlines for enclosure angles 0 ∘ and 60 ∘ with area ratios equal to 0. of these two parameters on temperature variation which in turn affects the rate of heat transfer. Thermal gradients are high near the hot and cold surfaces at all Rayleigh numbers. Isothermal lines are parallel to the hot and cold walls in the middle portion of the enclosure at low Rayleigh numbers and as the Rayleigh number increases these lines turn out to be almost perpendicular. The effect of angle of inclination of enclosure is noticeable in terms of isothermal lines plotted for 60 ∘ angle of inclination. The isothermal lines are closely packed near the adiabatic wall which is on the downstream side of the heat source.
The effect of area ratio for the same Rayleigh number and the angle of inclination can also be observed in Figures 5 and  6. Though the trends of the curves remain the same as the area ratio increases, additional concentration of the streamlines and isothermal lines is noticed, which indicates the increase in heat transfer rate with area ratio.

Variation of Average Nusselt Number.
The angle of inclination of the enclosure strongly influences the natural convection heat transfer particularly at high Rayleigh numbers, due to change in streamline and isothermal line patterns. For illustration, the quantitative variations of average Nusselt number with for different Rayleigh numbers with hot surface area ratio equal to 0.8 are shown in Figure 7. At low Rayleigh numbers 10 3 and 10 4 the flow circulation within the enclosure is weak and hence average Nusselt number is independent of its angle of inclination. At high Rayleigh numbers the average Nusselt number decreases as increases from 0 ∘ to 90 ∘ and attains its initial value corresponding to 0 ∘ when rises from 90 ∘ to 180 ∘ . When is increased beyond 180 ∘ and approaches 270 ∘ , the hot surface ascends above the centre of the enclosure and the cold surface moves towards the bottom level. The buoyancy force deviating away from the cold wall contributes additional drop in rate of heat transfer.
The lowest value of heat transfer rate is obtained at equal to 270 ∘ in which the hot and cold surfaces are on the extreme upper and lower sides, respectively, and the direction of the buoyancy force is exactly towards the hot wall. Increasing from 270 ∘ to 360 ∘ changes the average Nusselt number to its initial value which corresponds to 0 ∘ angle of inclination. The same trend of variation of Nu avg with is observed for all the Rayleigh numbers investigated with different area ratios of the heat source.

Heat Transfer Dimensionless Correlations.
In steady state condition all the heat flux ( ) from the heat source is absorbed by the cold wall and hence the Nu avg over the cold wall is a measure of heat transfer capacity of the enclosure. Based on the numerical results obtained in this study a correlation equation can establish the dependence of average Nusselt number on Rayleigh number and area ratio of the heat source, for a particular angle of inclination of the enclosure, as follows: The values of correlation coefficients , , and which are obtained by multiple regression analysis for all the investigated angles of inclination of the enclosure together with the maximum relative error and the average standard deviation avg (%) are detailed in Table 1 The values of and avg (%) vary between 0.0143 and 0.1523 and between 0.2216 and 3.567, respectively, which prove good agreement between the numerical results and the values obtained using correlation equation.

Conclusions
Natural convection heat transfer in a square enclosure heated from the discrete heat source and cooled from the opposite wall with the remaining walls insulated has been investigated. The studies are conducted for several area ratios of the heat source from 0.2 to 1.0, for the angles of inclination of the enclosure 0 ∘ to 360 ∘ , and for values of Rayleigh number in the range between 10 3 and 10 7 . Streamline and isothermal line patterns are found to be similar for all the cases at low Rayleigh numbers whereas considerable differences are observed at high Rayleigh numbers. Studies conducted on the influence of angle of inclination reveal that the average Nusselt number decreases drastically as the position of the heat source is moved above the horizontal centre line of the enclosure. Heat transfer correlation equations have 8 Journal of Applied Mathematics been developed for each angle of inclination of enclosure investigated which are in good agreement with the numerical results.

Nomenclature
: Area ratio of heat source ( × 1/ × 1) : Acceleration due to gravity : Length of the enclosure wall Nu: Nusselt number Nu avg : Average Nusselt number : Fluid static pressure Pr: Prandtl number : Heatflux Ra: Rayleigh number : Length of the heat source : T emperature : Temperature of cold wall ℎ : Temperature of the heat source , V: Velocitycomponentsin and directions , : Cartesian coordinates , , : Correlation coefficients.

Greek Symbols
: Thermaldiffusivity : Volume expansion coefficient : Angle of inclination of the enclosure ]: Kinematic viscosity : Fluid density : Maximum relative error avg : Average standard deviation.

Introduction
The problems of transient heat flow in hollow cylinders are important in many engineering applications. Heat exchanger tubes, solidification of metal tube casting, cannon barrels, time variation heating on walls of circular structure, and heat treatment on hollow cylinders are typical examples. It is well known that if the temperature and/or the heat flux are prescribed at the boundary surface, then the heat transfer system includes heat conduction coefficient only; on the other hand, if the boundary surface dissipates heat by convection on the basis of Newton's law of cooling, the heat transfer coefficient will be included in the boundary term.
For the problem of heat conduction in hollow cylinders with time-dependent boundary conditions of any kinds at inner and outer surfaces, the associated governing differential equation is a second-order Bessel differential equation with constant coefficients. After conducting a Hankel transformation, the analytical solutions can be obtained and found in the textbook byÖzisik [1].
For the heat transfer in hollow cylinders with mixed type boundary condition and time-dependent heat transfer coefficient simultaneously, the problem cannot be solved by any analytical methods, such as the method of separation of variable, Laplace transform, and Hankel transform. Few studies in Cartesian coordinate system can be found and various approximated and numerical methods were proposed. By introducing a new variable, Ivanov and Salomatov [2,3] together with Postol'nik [4] transformed the linear governing equation into a nonlinear form. After ignoring the nonlinear term, they developed an approximated solution, which was claimed to be valid for the system with Biot number being less than 0.25. Moreover, Kozlov [5] used Laplace transformation to study the problems with Biot function in a rational combination of sines, cosines, polynomials, and exponentials. Even though it is possible to obtain the exact series solution of a specified transformed system, the problem is the computation of the inverse Laplace transformation, which generally requires integration in the complex plane. Becker et al. [6] took finite difference method and Laplace transformation method to study the heating of the rock adjacent to water flowing through a crevice. Recently, Chen and his colleagues [7] proposed an analytical solution by using the shifting function method for the heat conduction in a slab with time-dependent heat transfer coefficient at one end. Yatskiv et al. [8] studied the thermostressed state of cylinder with thin near-surface layer having time-dependent thermophysical properties. They reduced the problem to an integrodifferential equation with variable coefficients and solved it by the spline approximation.
According to the literature, because of the complexity and difficulty of the methodology, none of any analytical solutions for the heat conduction in a hollow cylinder with time-dependent boundary condition and time-dependent heat transfer coefficient existed. This work extends the methodology of shifting function method [7,21,22] to develop an analytical solution with closed form for the heat transfer in hollow cylinders with time-dependent boundary condition and time-dependent heat transfer coefficient simultaneously. By setting the Biot function in a particular form and introducing two specially chosen shifting functions, the system is transformed into a partial differential equation with homogenous boundary conditions and can be solved by series expansion theorem. Examples are given to demonstrate the methodology and numerical results are compared with those in the literature. And last but not least, the influence of physical parameters on the temperature profile is studied.

Mathematical Modeling
Consider the transient heat conduction in heat exchanger tubes as shown in Figure 1. A fluid with time-varying temperature is running inside the hollow cylinder and the heat is dissipated by the time-dependent convection at the outer surface into an environment of zero temperature. The governing differential equation of the system is where is the temperature, is the space variable, is the thermal diffusivity, is the time variable, and and denote inner and outer radii, respectively. The boundary and initial conditions of the boundary value problem are Here, ( ) is a time-dependent temperature function at the inner surface, is the thermal conductivity, ℎ( ) is a timedependent heat transfer coefficient function, and 0 ( ) is an initial temperature function. For consistence in initial temperature field, (0) must be equal to 0 ( ). The above problem can be normalized by defining = , = , where is a constant reference temperature, and the dimensionless boundary value problem will then become = ( ) , at = , To keep the boundary condition of the third kind at outer surface in the following analysis, one sets the Biot function Bi( ) in the form of where and ( ) are defined as = Bi (0) , It is obvious that (0) = 0, and the boundary condition at = 1 can be rewritten as

Change of Variable.
To find the solution for the secondorder differential equation with time-dependent boundary condition and time-dependent heat transfer coefficient at different surfaces, the shifting function method [7,21,22] was extended by taking Substituting (11) into (4), (5), (10), and (7), one has the following equation: and the associated boundary and initial conditions now are Something worthy to mention is that (13) contains three functions, that is, V( , ) and ( ) ( = 1, 2), and hence it cannot be solved directly.

The Shifting Functions.
For convenience in the analysis, the two shifting functions are specifically chosen in order to satisfy the following conditions: Consequently, the shifting functions can be easily determined as Substituting these shifting functions and auxiliary time functions into (11) yields When setting = 1 in the equation above, one has the relation Therefore, two functions in governing differential equation (13) are integrated to one. With (16) and (18), (13) can be rewritten in terms of the function V( , ) as where ( ), ( = 1, 2) are defined as Meanwhile, the associated boundary conditions of the transformed function turn to homogeneous ones as follows: Since 2 (0) = − (0) (1, 0) and (0) = 0, hence, the associated initial condition can be simplified as

Series Expansion.
To find the solution for the boundary value problem of heat conduction, that is, (19)- (22), one uses the method of series expansion with try functions: satisfying the boundary conditions (21). Here the characteristic values are the roots of the transcendental equation The try functions have the following orthogonal property: where the norms are Now, one can assume that the solution of the physical problem takes the form of where ( ) ( = 1, 2, 3, . . .) are time-dependent generalized coordinates. Substituting solution from (27) into differential equation (19) leads to Expanding 1 ( ) and 1 ( ) on the right hand side of (28) in series forms we obtain where and are in which ( = 0, 1, 2, 3; = , ) are given as  After taking the inner product with try function ( ) and integrating from to 1, the resulting differential equation now iṡ( where and are and ( ) is The associated initial condition is As a result, the complete solution of the ordinary differential equation (33) subject to the initial condition (36) is where ( ) is After substituting (16), (18), (23), and (27) back to (11), one obtains the analytical solution of the physical problem where the summation is taken over all eigenvalues of the problem.

Constant
Heat Transfer Coefficient at = 1. When the heat transfer coefficient ℎ at = 1 is time-independent, the Biot function is a constant and ( ) = 0. The infinite series solution, (39), is reduced to where the generalized coordinates ( ) are The (0)'s for the problem under consideration are Introducing (42) in (41) and performing integration by parts, we can get Substituting (43) into (40) yields the temperature distribution: This solution is the same as that obtained via the integral transform method byÖzisik [1].

Verification and Example
To illustrate the previous analysis and the accuracy of the three-term approximation solution, one examines the following case. The time-dependent boundary condition ( ) considered at = is taken as and differentiating it with respect to leads tȯ where 1 and 1 are two arbitrary constants and 1 and 1 are two parameters. The Biot function considered at boundary = 1 is where 2 and 2 are two arbitrary constants and 2 and 2 are two parameters. According to (8)-(9), we obtain Journal of Applied Mathematics Consequently, the temperature distribution in the hollow cylinder is where the ( )'s are defined in (37). The associated ( ) now is To avoid numerical instability that occurred in computing ( ), (37) is rewritten as Since the initial conditions cannot have effect on the steady-state response, we consider only the heat conduction in a hollow cylinder with constant initial temperature 0 ( ) = 0 as prescribed in the previous sections. The (0)'s are now computed as For consistence in the temperature field, the constant 0 is taken as zero in the following examples.
In comparison with the literature, the example of constant Biot function is studied first. Bi( ) = 1 and time-dependent temperature function, ( ) = 1 − − , are chosen in the case. In Table 1, we find that the convergence of the present solution is faster than that ofÖzisik [1]. The error of threeterm approximation in present study is less than 0.4%; on the contrary, at least twenty-term approximation is required to get the same accuracy inÖzisik's [1] cases.
In the case of time-dependent boundary condition and time-dependent heat transfer coefficient at both surfaces, we consider the time-dependent temperature function, ( ) = 1 − − cos , and the Biot function, Bi( ) = 2 − − . From Table 2, one can find that the error of three-term approximation is less than 0.4%. Because of large values of Bi( ), the internal conductance of the hollow cylinder is small, whereas the heat transfer coefficient at the surface is large. In turn, the fact implies that the temperature distribution within the hollow cylinder is nonuniform. Therefore, we find that the larger the Biot function, that is, when approaches to 10 in Table 2, the more the iteration numbers. Figure 2 depicts the temperature profiles along the radial of the hollow cylinder at different times, = 0.1 and = 1. We find that the temperature at = 0.6 is higher than the temperature at = 1 and the temperature profile decreases at the negative slope for every case. It is clear since the heat    function ( ) the temperature in Bi( ) = 2 − −2 is less than that in Bi( ) = 2 − − as proceeds. That is to say, more heat will be dissipated into the surrounding environment for Bi( ) = 2 − −2 as goes. Figure 4 depicts the effect of the parameter 1 of temperature function ( ) upon the temperature variation of the hollow cylinder. It is found that, in the same temperature function ( ), the temperature for = 0.4 is less than that for = 1.0. Besides, as increases from 0.4 to 1.0, the difference between temperatures at 1 = 1 and at 1 = 2 becomes significant.
Periodical heat source versus time-varying Biot function is drawn to show the temperature variation of the hollow cylinder at = 0.8 and = 1.0 with respect to in Figures  5(a) and 5(b), respectively. Two cases of heat source ( ) = 1 − cos (solid lines) and ( ) = 1 − cos 2 (dash lines) are considered. At the same , the temperature of Bi( ) = 2− − is less than that of Bi( ) = 2− −0.1 for constant ( ). The reason is that more heat has been dissipated into the surrounding environment at the case of Bi( ) = 2 − − . It can be observed that as proceeds, in the beginning, the temperatures are nonsensitive with 1 parameters, as shown in Figure 5.

Conclusion
An analytical solution for the heat conduction in a hollow cylinder with time-dependent boundary conditions of different kinds at both surfaces was developed for the first time. The surface is subject to a time-dependent temperature field at inner surface, whereas the heat is dissipated by timedependent convection from outer surface into a surrounding environment at zero temperature. The methodology is an extension of the shifting function method and the present results are identical to those in the literature when constant Biot function is considered. Since the methodology does not use integral transform, it has a proven result. The proposed method can also be easily extended to various heat conduction problems of hollow cylinders with timedependent boundary conditions of different kinds at both surfaces.

Introduction
Research on high-performance aeroengines shows that the nozzle is always the most severely heated structure of the engine. The gas temperature of a vectoring nozzle under afterburning condition is higher than 2000 K [1]. To prevent heating damage on the nozzle wall, film cooling has been widely applied on the vectoring nozzle in the past decades. The infrared radiation characteristics of a deflection nozzle are important for the engine to achieve infrared stealth and radar tracking [2]. In the vectoring nozzle of the aeroengine, computation of the wall temperature is a complex nonlinear problem that involves multiple coupled variables, including gas radiation and wall temperature. The major radiative gases, H 2 O and CO 2 , are all strongly spectral dependent, and their infrared radiation is composed of millions of spectral lines. The waveband model of gas spectral radiation is generally computation intensive. The deflection of the nozzle breaks the axisymmetric distribution of gas flowing parameters and radiation, which also increases the difficulty and calculation of heat transfer analysis. Thus, the accurate coupled calculation of the film-cooling wall temperature and gas infrared radiation is a problem that must be solved in the heat protection optimization of a high-performance aeroengine.
The technology of applying film cooling on the convergent part of a nozzle was first studied on a rocket engine. Reference [3] demonstrated that when the pressure of a combustion chamber is higher than 5000 lbs/in, only 40% of the wall heat flux can be removed by regenerative cooling; the use of film cooling to enhance cooling also proposed. Reference [4] proposed the Hatch-Papell equation to calculate film cooling in a rocket engine and used empirical methods to calculate the radiative heat. Reference [5] experimentally researched the influence of the geometry of the film hole on the film-cooling heat transfer coefficient. Reference [6] studied a combustion chamber with a wiggle strip film-cooling construction using empirical correlations of gas emissivity to compute the wall temperature distribution. Reference [7] focused on the wall temperature distribution of multiple slots of film cooling on heat shields to account for the heat conduction at the conjunction surface by using gas emissivity to compute radiative heat flux. Reference [8] numerically studied the mixing of the film and main flows on a plane with a single-film slot. Reference [9] researched the coupled computation methods of flow and the heat transfer of gas, wall, and coolant flow in a thrust chamber of a high-thrust liquid rocket engine. Reference [10] adopted the commercial software FLUENT to simulate the flow field and wall temperature of a round to square nozzle. Reference [11] researched the heat transfer model of an axisymmetric divergent nozzle and considered the conductive, convective, and radiative heat. Lefebvre empirical methods were adopted to compute the gas radiation to the wall.
Based on the spectral characteristics of aeroengine gas radiation, [12] adopted a weighted-sum-of-gray-gas model to simulate gas radiation; it used finite volume methods to compute the radiation transport equation. Reference [13] developed a 2D simulation model of radiation that uses a spherical harmonic method of solving the directed integral form of radiance equation. Given the complexity of gas spectral radiation calculation, few published studies have considered gas spectral radiance in analyzing the engine wall temperature. The cooling of new high-performance aeroengines and infrared stealth technology require the accurate prediction of the radiative heat transfer between the nozzle wall and gas. The present work established a theoretical analysis model that focuses on the entire structure of an aeroengine axisymmetric vectoring exhaust nozzle (axisymmetric vectored nozzle) and the convective/radiative heat transfer of hot gas. The temperature distribution of the entire nozzle structure, including heat shield, nozzle wall, and outer shield, is computed, as well as the infrared radiation on the wall and outlet surface of the nozzle.

Geometric Model and Grid
The vectoring nozzle with a multirow of film cooling investigated in this work is a 0.624 m long axisymmetric nozzle with a 0.8 m diameter afterburning part. The convergent part of the nozzle is 0.184 m long. The heat shield is set on the nozzle, with a gap of 0.006 m from the wall. Three slots of film cooling are applied on the heat shield. The coolant flows through the slots and covers the heat shield surface, thereby forming the adherent film. A part of the coolant directly goes through the end of the gap, rushes to the divergent part of the nozzle, and quickly mixes with the main flow. The outer shield outside the nozzle has a space of 0.1 m from the nozzle wall. The secondary flow goes through the space and conveys some of the heat on outside surface of the nozzle wall. The entire structure of the nozzle is shown in Figure 1. The diagram of the three slot of film is shown in Figure 2, where the three slots are equally spaced along length of the heat shield. Height of each slot is 1.5 mm.
The calculation of the wall radiation heat and film-cooling heat transfer in this work is based on the flow field results of the computational fluid dynamics (CFD) software FLUENT. The CFD solution adopts an unstructured grid. We have done a series of calculations based on different grid resolution to get the grid-independent solution. The final grid resolution we adopted is about 107550 grids on -coordinate plane, spaced 0.008 m apart in nozzle axis. The near zone + value is about 30∼90. The grids system for CFD solution is not suitable for heat transfer and wall temperature computation, because of the irregular and too dense arrangement of nodes, particularly at the boundary layer of the wall. The density of grid nodes slows down the calculation of gas radiation. Therefore, we reset the grids for radiation calculation. The rearranged grids are divided with the equal space idea on 51 sections. The arrangement of the nodes in each section is shown in Figure 3. The gas properties at each node are obtained through a spatial interpolation method. The new set grids are quadrangular and regular, which facilitates the calculation of heat transfer. The final reset grids for computing heat transfer and wall temperature of the nozzle are shown in Figure 4. [14]. An enclosure radiation model is used to calculate the hot gas radiation to the wall. An enclosure is established; it consists of the included hot gases, nozzle wall, inlet surface, and outlet surface. The radiation of each wall unit of the enclosure is calculated. The net radiative heat transfer rate of each unit is given as follows:

Calculation Model of Gas Radiation
where is the emissivity of each unit. is the temperature of the wall unit. is the radiosity of each unit. The radiosity of unit consists of two parts: the radiation from the unit itself and the reflection of incident radiation . Thus, where is the wall reflectivity. Incident radiation comes from radiosity of all the other wall cells in the enclosure. Therefore, the radiosity of all wall cells should be correlated with one another and must be solved simultaneously.
On any of the wall cells, the incident radiation consists of the gas radiation of the incident hemisphere and the radiation from the other wall cells. To calculate , the integral of the incident radiation in each direction of the hemisphere is needed. The incident hemisphere of a wall unit Δ is divided into several solid angles . Each solid angle corresponds to a unit surface Δ on the wall. The incident radiation of each solid angle is calculated.
Their sum is used to obtain the total incident radiation of a wall unit. The incident radiation of each solid angle consists of gas radiation along the path and radiosity from the emitting wall unit Δ .
The incident radiation at each solid angle is the integral of the gas layer on the incident path, expressed as follows: where , represents the incident radiation along the connection of wall cells Δ and Δ . The symbol is the optical thickness between the two wall cells, * is the optical thickness from one specific calculation node to the incident wall, and * is the optical thickness within each calculation gas layer.
( * ) represents the black body radiation at the gas temperature. ( ) represents the gas transmittance of optical thickness .
is the gas layer number on the connection of the two wall cells (i.e., the number of gas layers in radiation integration). , , indicates the part of the radiosity from the wall unit Δ that goes through the intermediate gas to the wall unit Δ .
Considering the spectrum of gas radiation, the infrared radiation at any direction at the spectrum from 1 m to 5 m is divided into 41 discrete wavebands. Gas radiation in all the other spectrums is regarded as a transparent body. Thus, 42 groups of spectral radiosity , , can be integrated to give radiosity in all wavelength.
For the other waveband outside the infrared waveband, a simplified equation ignoring the gas radiation is given as follows: In each waveband, the integral of (3) in all the solid angles of the hemisphere of Δ is substituted into (2). Assume that the radiation angle factor is = cos cos Δ / 2 , so the radiosity of one wall cell Δ can be given as follows: Journal of Applied Mathematics where is the number of the infinitesimal solid angles in the incident hemisphere of the wall cell Δ . The radiosity equations of all the wall cells in the enclosure are a closed set of equations: The coefficient matrix of this equation is diagonal dominant and can be solved using the Gauss-Seidel method. The radiosity of all the wall cells in every waveband is obtained, and their sum is the total radiosity. Thus, the relation between the net radiation heat and wall temperature is obtained according to (1).

Calculation Model of Film Cooling.
A cooling effectiveness model is adopted to calculate the film cooling on the convergent part of the nozzle. The convective heat transfer rate of gas is computed according to the imperial correlation given by [14,15].
The convective heat transfer rate of film cooling is given as follows: where is the adiabatic wall temperature of film cooling, is the temperature of the nozzle wall, and ℎ is the convective heat transfer coefficient. The common solution to ℎ is given by the forced convection heat transfer Dittus-Boelter formula. The physical parameter used in the formula is the main flow parameter near the wall.
is given by the empirical correlation of the cooling effectiveness. Consider where is the temperature of mail flow in gas.
In this work, three cooling slots are applied on the heat shield. To calculate the compound cooling effectiveness of the three slots, the cooling effectiveness of multiple slots is calculated based on the result of a single slot. The cooling effectiveness equation of a single slot is given in [14] as follows: where = film film / is the film blowing ratio, is the distance to the slot, is the slot height, and 0 is function of and / . When the calculation point is close to the injection slot, the cooling effectiveness given by (10) is greater than 1, which is unreasonable. Under this condition, the Hatch-Papell equation is used instead, as recommended by [4]. In [4], the Hatch-Papell equation was applied on a rocket engine's remaining high precision. The Hatch-Papell equation is rewritten in metric form as follows: The compound cooling effectiveness of the three slots of film is given by [14] as follows: where 1 , 2 , and 3 are the cooling effectiveness of each slot.

Calculation Model of Wall
Temperature. Lumped parameter methods are adopted to calculate the distribution of the wall temperature of all the nozzle structures. A heat balance equation of convective and radiative heat transfer is established on every wall unit as follows: where the subscript represents the convective heat, means the radiative heat, means the inner side of the wall, and means the outer side of the wall. The calculation of the convective and radiative heat of the wall cells whose inner side is directly heated by gas is discussed in the section above. The convective heat of other wall cells is calculated using the Dittus-Boelter equation. The radiative heat is calculated using the following equation of two envelope surfaces: The convective terms of the heat balance equation of the wall cells not directly heated by gas can also be solved using the Dittus-Boelter equation.
Considering any wall cell, the radiative term computed with (13) in its heat balance equation involves the temperature of the other wall cells. Thus, these heat balance equations are coupled with one another. Assuming 1 , 2 , and 3 are the temperatures of heat shield, nozzle wall, and outer shield, Journal of Applied Mathematics 5 respectively, 1 , 2 , and 3 are their heat balance equations; then the coupled heat balance equations group of the three structures above can be given as follows: Equation (14) is a set of multivariate nonlinear equations. These equations can be established on every wall cell of a nozzle. These equations are solved using the Newton-Raphson methods in our program, which performs iteration in each coupled equation.
In the iteration of (14), the radiative term of the gas-heated wall correlates with its radiosity and temperature. Also the radiosity of a wall cell correlates with the temperature of the other wall cells. Hence, the calculation of the radiosity and temperature of the wall cells should be simultaneous. When the temperature of a wall unit changes, iteration of radiosity must be repeated to update the radiosity of all wall cells. However, the waveband model for gas radiation consumes considerable time to complete the radiosity computation. Thus, the real-time coupling of temperature and radiosity is difficult to achieve.
The solution we adopt is stepwise and round coupling. First, a set of initial temperatures of the wall cells is assumed. According to these initial temperatures, the heat balance equation is iterated, ignoring the gas radiation. A set of wall temperatures is obtained. Second, the iteration of radiosity is performed according to the wall temperature above. The gas radiative heat transfer rate of the wall cells is generated. Thirdly, according to the newly generated radiative heat, the iteration of the wall temperature is restarted. A set of wall temperatures that consider the gas radiation is obtained and regarded as the result of first round. A new round of iteration will be restarted if temperature difference between two rounds is greater than a minimum limit. Iteration errors of both heat balance and gas radiation get smaller in a new round of iteration.

Results and Analysis of the Verification Case.
To verify theoretical and numerical model of this work, a film-cooling aeroengine [6] and a rocket engine [4] are investigated. Reference [6] measured the temperature distribution from the fourth to the end of the sixth wiggle strip of a flame tube of Spey-25 engine. The flame tube of this engine used film cooling on the wiggle strip. The result of this work is compared with that of [4] in Figure 5, where nondimensional temperature is plotted against nondimensional distance. Nondimensional temperature is the ratio of the temperature of the calculation point and the maximum temperature. Nondimensional distance is the ratio of the axial distance to the inlet and the nozzle length. Figure 5 shows that the results of the reference and this work reflect the cooling effect of film cooling. The wall temperature near the injection point is low and increases downstream, where three waves are found near   three slots. Cooling effect of slot film is the strongest just at the injecting position, which will be gradually attenuated with distance increasing, as coolant will be heated in mixing with the gas. A new rise of cooling effect will come where a new slot appears.
The measured curve is more gradual in the said reference than in this work. The conductive heat transfer inside the wall is ignored in this work to simplify the calculation. However, high-temperature wall cells transfer heat to the neighboring low-temperature cells, thereby resulting in a gradual distribution of temperature. The difference between the results of the reference and this work is no more than 20%. Therefore, the result of this work can properly reflect the variation trend of the wall temperature under film cooling.
Reference [4] gives the measured temperature distribution of a rocket engine under ground test condition. One slot of film cooling was applied on a slight upstream of the nozzle throat. A comparison of this work and the wall temperature 30 s after ignition in [6] is shown in Figure 6. The two groups of results basically agree with each other. However, the result of this work has a maximum temperature that is closer to the throat. The axisymmetric vectored nozzle in this work is deflected after the throat, thus causing temperature to ascend after the throat because of the scour of hot gas.

Temperature Distribution of Axisymmetric Vectored Nozzle.
For the vectoring nozzle described in Section 2, this work calculates the temperature of its convergent part, as shown in Figure 4. The definition of nondimensional distance is the same as above. The temperature distribution of nonfilm cooling is given in Figure 7 for an intuitive comparison. Under the protection of film cooling, the wall temperature is significantly low. The fluctuation of wall temperature becomes great after applying the cooling film. The wall temperature near the coolant injection point can be 550 K lower when no film is applied. Under the isolation of the heat shield, the temperature of the convergent part of the nozzle is maintained at nearly 500 K, which is approximately 500 K lower than that of the heat shield. While it is cooled by the film coolant in the gap of the heat shield, the nozzle wall does not have direct contact with the gas. The heat shield, accordingly, plays a good part in cooling the nozzle wall. The temperature distribution on each structure of the nozzle is shown in Figure 8. The divergent part is more severely heated than the other parts of the nozzle. Without the protection of the heat shield, this part of the nozzle is directly exposed to the hot gas. Although the gas temperature after the nozzle obviously decreases, the wall temperature ascends to 940 K. Meanwhile, the average temperature of the outer shield increases from 440 K to 690 K.
The main effects of vectoring condition on temperature distribution are the uneven distribution of gas parameters after the throat and the occlusion of the throat to the upper wall cells after the throat. The wall temperature distribution of vectoring condition in three circumferential positions is shown in Figure 9, where is the circumferential angle to the deflection direction. The temperature of the deflection  direction ( = 0) is the lowest, whereas that of the opposite direction ( = ) is the highest. The opposite direction is directly heated by gas. This part of the wall cells is not blocked by the throat and is heated directly by the inlet surface radiation. The average wall temperatures of these two directions differ by 4.8%.

Comparison of Deflection and Nondeflection Conditions.
The wall temperature distribution of the nondeflection condition is also obtained in this work. Compared with the deflection of the 20 ∘ condition, as shown in Figure 10, the temperature distribution pattern at the divergent part (circumferential = 0) may be different. The temperature curve of the deflection condition is arched because of the scour of   Figure 11: Radiosity in infrared and other wavebands on the wall. hot gas, whereas the curve of the nondeflection condition exhibits a linear decline. The decline of the nondeflection condition is fast, and the maximum value appears at the throat. By contrast, the wall temperature of the deflection condition is maintained at a high value at the downstream of the throat, and the maximum temperature appears at a distance from the throat. Figure 11 shows the radiosity at an infrared waveband of 1 m to 5 m and other wavebands in the circumferential angle of = 0, /2, and under vectoring condition. The radiosity at a waveband of 1 m to 5 m is 10 times greater than radiosity at other wavebands. Therefore, the radiation of the nozzle is mainly concentrated in the infrared waveband. The radiosity in the convergent part is circumferentially uniform, whereas that in the divergent part has a minimal value when = 0.

Infrared Characteristics of Axisymmetric Vectored Nozzle.
The maximum value appears at = . The radiosity at the throat is lower than that on the wall on both sides because of the contraction of the flow area and the deflection of the nozzle, which block the radiation of nearby gas of the other side. The nearby gas has the most significant effect on the wall, thereby decreasing gas radiation near the throat. Figure 12 shows the spectral radiation at an infrared waveband of 1 m to 5 m on the outlet surface of the divergent part in the = 0, /2, and direction. The spectral radiation on the outlet surface shows a characteristic of gas spectral radiation. When the spectral peak of H 2 O and CO 2 radiation at 1.8, 2.7, and 4.3 m is reached, the spectral radiosity is enhanced. The circumferential difference of the radiosity has the same pattern as that of the wall temperature (i.e., the minimum value appears at = 0, and the maximum value appears at = ).

Radiation
Intensity at the Outlet. The radiation observed from the rear hemisphere of the nozzle outlet is also computed in this work. Figure 12 shows that the radiation intensity observed from the rear hemisphere is plotted against the emission angle along the upright direction. Figure 13 shows that directed radiation has a maximal value at an emission angle of 70 ∘ , which is the deflection angle of 20 ∘ . By contrast, the directed radiation distribution under nondeflection condition has the maximum value at the right direction of 90 ∘ . Thus, AVEN remarkably changes the emission angle of the directed radiation.

Conclusion
With the aim of designing the cooling structure design of the AVEN of a high-performance engine, this study established a heat transfer-spectral radiosity coupled heat balance equation model. A 3D all-structure nozzle simulation program of high accuracy was also developed by considering gas radiation heat transfer and film cooling. Our research obtained the following conclusions.
(1) The validation result of this work agrees well with the measured result in a reference, which demonstrates that this work may comprehensively reflect the heat transfer condition in the nozzle.
(2) The convergent part of the axial symmetric vectored nozzle studied in this work is shielded by the heat shield and cooling film. Thus, the wall temperature of this part is relatively low. By contrast, the wall temperature of the divergent part is high and should be the emphasis of the cooling design. The heat shield is the most severely heated structure of the nozzle, a phenomenon that must be considered by the design.
(3) The deflection of the axial symmetric vectored nozzle increases the wall temperature of the divergent part in the opposite direction of deflection and makes it last for a long distance at the downstream, which results in a severely heated area. The deflection also changes the distribution of the directed radiation at the outlet effectively.