Ablative degradation of cryogenic thermal protection and fuel boil-off: Improvement of using graded density insulators

International Journal of Heat and Mass Transfer 54 (2011) 4864–4874

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

Ablative degradation of cryogenic thermal protection and fuel boil-off: Improvement of using graded density insulators Jaona Randrianalisoa ⇑, Remy Dendievel, Yves Bréchet SIMAP, CNRS, Grenoble INP, UJF, F-38402 Saint Martin d’Hères, France

a r t i c l e

i n f o

Article history: Received 30 March 2011 Received in revised form 24 June 2011 Accepted 27 June 2011 Available online 26 July 2011 Keywords: Heat transfer Cryogenic insulation Polymer Foams Ablation Boil-off Nucleate boiling

a b s t r a c t A simple modeling of the heat transfer through insulating bilayers of cryogenic reservoir subjected to surrounding environment or aerothermal flux is presented. The model permits to determine the instantaneous evolution of maximum fuel temperatures and the boil-off rate as well as the possible insulation ablation. To get a realistic estimation of the maximum fuel temperature and boil-off loss, the good knowledge of the internal convection coefficient, taking into account the possible fuel ebullition, is indispensable. The aerothermal flux provokes the insulation degradation and its ablation. Because of the ablation, the reservoir becomes less and less thermally protected, which may cause an important fuel boil-off. Compared to classical insulations with uniform density, those with spatially graded densities appear as potential insulations permitting to meet both light-weight and low boil-off requirements. Ó 2011 Elsevier Ltd. All rights reserved.

1. Introduction In many cryogenic applications, the design requirements impose to maintain a low temperature inside the reservoir, while the external atmosphere is at much higher temperature. In extreme situations such as rockets reservoir, the inner temperature must remain below the boiling of liquid hydrogen (LH2), i.e. 20 K at atmospheric pressure, while the external temperature can go up to several hundred Kelvins [1]. In transportation applications, the need to develop weight-light structures has led to foams as the insulating material while the structural function is carried out by a metallic shell [2]. The large dimensions of the reservoir (about several meters of diameter and several 10’s of meters of height) preclude any monolithic solution for the insulation, and therefore the insulating component consists of plates glued on the metallic shell. In highly demanding situations such as rocket launchers, the thermal insulation is submitted to degradation on the external side, due to severe aerothermal fluxes. As a consequence, its insulation capability decreases while time goes on. The design of these insulating plates must therefore fulfil conditions of low ebullition of the fuel in the reservoir even with a progressively disappearing protective layer. In this aim, the knowledge of the maximum fuel ⇑ Corresponding author. Current address: GRESPI Lab., Reims University, F-51687, France. E-mail address: [email protected] (J. Randrianalisoa). 0017-9310/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijheatmasstransfer.2011.06.042

temperature and the ebullition rate for a given insulation system and given external conditions is of crucial importance. A solution of this problem can be obtained using full numerical analysis (e.g. coupling computational fluid dynamics (CFD), for the fuel [3,4], and FE calculations, for the heat transfer through the thermal protection [5,6]). Such approaches are however computationally expensive and they are not necessary the best strategy for optimization purpose. In the current work, the main objective is to predict the maximum fuel temperature, which will determine the boil-off (or liquid fuel loss by boiling). A detailed analysis of the entire temperature field of the fuel would give such information. However, the knowledge of the fuel temperatures close to the warmer metallic surface can be sufficient due to the temperature continuity at the fuel/wall interface. The main concern is therefore the temperature distribution at any point of the internal metallic surface, which requires only solving the temperature field within the protecting layers. This simplification is especially welcome in a design spirit where the insulating material is a foam type which can span a large range of density, thickness, and thermal properties. In addition, the external conditions on the protective layer lead to degradation and possibly ablation processes of the polymer [7]. Since the continuous variable, namely the density, influences both thermal properties and the ablation rate of the foam, the optimal choice for this variable is by no way obvious and requires a detailed analysis. The final thickness of the insulating layer, driven by a minimum weight requirement, and the need to avoid or minimize ebullition

J. Randrianalisoa et al. / International Journal of Heat and Mass Transfer 54 (2011) 4864–4874

4865

Nomenclature C d D e fs Fsub g Gr  h, h H k kB K⁄, Ks L l⁄ _ m M Nu P Pr Qext R Ra RBO S t T x z

specific heat capacity (J/kg/K) diameter of the gas particles (m) cell size in a foam (m) thickness (m) solid volume fraction in a foam subcooling correction factor acceleration due to gravity (m/s2) Grashoff number exchange coefficient and average exchange coefficient, respectively (W/m2/K) height of cylindrical reservoir (m) thermal conductivity (W/m/K) Boltzmann constant (m2 kg/s2/K) extinction coefficients of foam and polymer, respectively (1/m) latent heat (J/kg) characteristic length in Eq. (18) (m) mass loss rate or boil-off rate (kg/s) initial mass of the fuel (kg) average Nusselt number pressure (N/m2) Prandtl number external heat flux (W/m2) radius of cylindrical reservoir (m) Rayleigh number mass loss ratio or boil-off ratio defined in Eq. (5) wetted reservoir surface (m2) time (s) temperature (K) radial axis of cylindrical coordinate system or space variable through the reservoir thickness (m) longitudinal axis of cylindrical coordinate system (m)

Greek symbols a thermal diffusivity (m2/s)

in the reservoir at any time, will also depend on the density of the selected foam. This paper intends to present a simple modeling of heat transfer through wall of cryogenic reservoirs taking into account the possible insulation ablative degradation. The paper is structured as follows. In Section 2, we will describe the problem configuration and the hypotheses used here. In Section 3, the detailed equations derived from our simplifying hypotheses will be presented. Special attention will be given to the formulation of an ablation model. In Section 4, the thermophysical properties are described. In Section 5, typical results of temperature profile through the reservoir wall, the ablative degradation of insulation, the maximum fuel temperature, and the fuel boiling–off ratio are shown. The critical foam density ensuring safe operating conditions taking into account the processes is analyzed, and the possibility of graded materials as far as their density is concerned is considered.

volumetric expansion coefficient (1/K) duration of flight (s) temperature difference (K) fraction of solid in a cell mean-free-path (m) m cinematic viscosity q, q⁄, hq⁄i density, kg/m3, relative density, and average relative density, respectively h tangential axis of cylindrical coordinate system (m) r Stefan–Boltzmann constant (kg/s3/K4), or surface tension (N/m2) b Dt DT / K

Subscripts cond refers to solid and gas thermal conductivity cr refers to heat flux or wall excess temperature at the end of nucleate boiling deg refers to degradation ext refers to quantities outside the reservoir f refers to fluid phase foam refers to foam properties i refers to insulating layer int refers to quantities inside the reservoir m refers to mean temperature NC refers to natural convection heat transfer NB refers to nucleate boiling heat transfer pad refers to time of ground phase rad refers to radiation contribution s refers to metallic shell or solid phase sat refers to saturation temperature or wall excess temperature v refers to vapor phase vap refers to vaporization w refers to foam cell wall thickness 0 refers to fuel initial temperature 1 refers to unperturbed fluid temperature

are LH2 and liquid oxygen (LO2). Since the functioning temperature of LO2 is much greater than that of LH2, the thermal insulation of the LH2 reservoir is most critical and is analyzed in this study. The LH2 compartment is filled up to a height of about

2. Geometry description and simplifying hypothesis A typical configuration of cylindrical cryogenic tank is described in Fig. 1. The main body, made of aluminum alloy, is primarily devoted to sustain pressure and all the other flight loads for space transport application. Domes are used to close the cylindrical main body at both ends. The external metallic surface is recovered by a layer of insulating panels to avoid heat penetration and hence fuel boiling. Classical cryogenic fuels or cryogens

Fig. 1. Lower-half part of a thermally insulated cylindrical tank containing cryogenic fuel at pressure P and temperature T and subjected to an external heat flux Qext.

4866

J. Randrianalisoa et al. / International Journal of Heat and Mass Transfer 54 (2011) 4864–4874

15 m and there is unfilled (vapor) space necessary for the fuel pressurization. The inner tank radius, namely R (about 5 m), the metallic shell thickness, namely es (about 2 mm) and the thickness of the insulating plates, namely ei, are supposed to be constant along the reservoir height. Because the wall thickness is much smaller than other dimensions, the major part of heat transfer will occur along the thickness direction. The problem is then solved by analyzing only the onedimensional (1D) heat transfer through a unit surface, for example a portion of the system comprising an insulating panel glued on a metallic plate, as depicted in Fig. 2. In this figure and in herein calculations, the presence of the adhesive layer is neglected. The heat exchange between the fuel and the tank wall is considered as a boundary condition of the problem with a convection exchange coefficient and the unperturbed fuel temperature as parameters (both may depend on the location and time) [8,9]. The system is subjected to external heat flux, namely Qext, which depends on the surrounding environment conditions. During the ground phase, convection at the external surface with free air and radiation transfer between the external surface and the surrounding medium are the main heat transfers. The radiation flux is linearized enabling us to treat the radiation transfer in a similar manner as the convection. During flight phase the aerothermal flux constitutes the major heat source. When the insulation reaches the polymer degradation temperature, ablation processes are assumed to take place. In fact, when a thin layer of the insulation undergoes degradation (by fusion, thermal decomposition or combustion), its mechanical strength can be strongly affected [10]. In presence of strong flowing air during flight, this leads to ablation processes [7,11]. Such mechanism is certainly very complex: for example, while being heated, a polymeric foam will successively soften, melt, then burn. In interaction with the air flow in contact with the external wall, each of these processes can lead to ablation. The exact rate at which this will occur depends both on the air flow, on the surface irregularities produced by the foam structural evolution, and on the remaining strength of the foam in any of the three situations listed above. Without entering the details of the process which would require an extensive and difficult experimental characterization, it is reasonable to model the ablation processes via a critical energy input necessary for the foam to reach a given state (‘‘softened’’, ‘‘molten’’ or ‘‘burned’’). This is more reasonable than assuming simply a critical temperature, since the amount of energy required to ‘‘melt’’ (or ‘‘burn’’) the foam depends on the quantity of matter involved. It is therefore more accurate to consider latent heats (of fusion, of combustion) and the foam density as the relevant physical parameters.

3. Governing equations The main unknown of the current problem is the temperature field across the insulating and metallic components. The energy balance equation in these two components along the x direction (according to Fig. 2) can be written as:

qC

  @T @ @T ¼ k @t @x @x

ð1Þ

with 0 6 x < es for the metallic wall and es 6 x 6 es + ei for the insulating layer. In (1), T is the temperature at the abscise x and at time t. q, k, and C are respectively the density, thermal conductivity and specific heat of the wall component at x and vary with the temperature. Within the insulating layer es 6 x 6 es + ei, the thermal conductivity k in (1) may include both pure conduction and radiative contributions as it will be discussed later. 3.1. Initial temperature profile A few minutes before launch, the tank is filled, then pressurized and brought to the flight pressure values. Typically, the LH2 functioning pressure P and temperature T are respectively, about 2.9 bars and 20.4 K [12]. This operation aims at stabilizing the tank structures, which will have to withstand significant mechanical loads during the ascent phase, and at providing a sufficient net positive suction pressure (difference between the static pressure in the tank and the saturation pressure in the pump) to prevent any cavitation in the turbo pumps during the draining (engine boost) phase. Generally, the fueling processes include a pre-cooling of the reservoir wall so that fuel temperature is not much altered with respect to the injection temperature (about 20.4 K for LH2). The initial condition we consider here, is when the fueling is achieved. At t = 0, the fuel can be thus assumed at injection temperature, T0, while a temperature gradient is established within the reservoir wall. The temperature field inside walls is not known a priori since it results from the cooling of the reservoir during the fueling. Theoretically, this temperature gradient depends on the external condition, the thermophysical properties of materials, and the fueling duration. It is however clear that to limit excessive boiling just after the fueling step, the temperature of the internal metallic surface should be as close as possible to the fuel temperature T0. In this study, the initial temperature distribution within walls is taken as the steady-state solution of (1) resulting from the following boundary conditions: imposed temperature T0 at internal wall due to the fueling processes and natural convection and room radiation [detailed in Eq. (6) below] at the external insulation surface. 3.2. Boundary conditions at the internal tank surface

t

Thickness reduction over time, t Cylinder axis z

Outward normal to the reservoir surface, x

x =0

External flux

hint

Q ext ei

Fuel at T int

es

Fig. 2. A portion of the cylindrical walls subjected to internal convection and external flux Qext provoking the insulating layer ablation at instant t. Dashed zone: reacted phase; doted area: functioning insulating layer; and hatched area: metallic wall. The adhesive layer is neglected.

Heat exchanges between the fuel and the internal metallic surface are assumed to be driven by convection. Two cases have to be distinguished depending on the temperature at the metallic surface T(x = 0).  When this temperature does not exceed the fuel vaporization temperature Tvap, free convection is assumed occurring in the reservoir. The liquid fuel is heated by the incoming flux through a natural fluid circulation. This can be represented by the following expression:

k

@T ¼ hint ðT  T int Þ at x ¼ 0 @x

ð2Þ

where hint is the exchange coefficient at the fuel/metallic wall interface and Tint is the fuel temperature far from the wall called also hereafter as fluid bulk temperature. Strictly speaking, Tint

J. Randrianalisoa et al. / International Journal of Heat and Mass Transfer 54 (2011) 4864–4874

varies with time due to the evolution of the internal wall temperature T. But, due to the large volume of the fuel that is contained inside the tank and the short duration storage, the bulk temperature Tint will not deviate far from the initial temperature T0. It is thus reasonable to assume in this study that Tint in (2) remains identical to T0.  When the wall temperature is greater than the liquid vaporization temperature Tvap, the incoming flux is released without heating the bulk fluid. Indeed, the fuel is heated through the contact in such a way that the supplied heat can only vaporize the fuel through pool boiling without increasing the temperature of the bulk liquid. The convection condition writes then:

k

@T ¼ hint ðT  T vap Þ at x ¼ 0 @x

ð3Þ

The expression of the exchange coefficient in (3) may include free and nucleate boiling mechanisms, as discussed later in Section 4.4. In this case, the mass loss due to ebullition is determined from the balance between the incoming heat flux and the energy released to vaporize the fuel:

_ vap ¼ Shint ðT  T vap Þ mL

ð4Þ

_ (in kg/s) is the mass loss of the fuel per unit time known as ‘‘mass m loss rate’’ or ‘‘boil-off rate’’. Lvap is the latent heat of vaporization of the fuel (Lvap = 446.6 kJ/kg for saturated LH2). S is the wetted surface of the reservoir. This surface decreases with time when the fuel undergoes ebullition and more particularly when the fuel flows out the tank (e.g. consumed by the cryogenic engine and converted to hydrogen gas for maintaining the vapor space pressure). Using (4), an estimate of how the initial fuel volume is affected by the heat source for a given insulation system and time can be determined. The boil-off rate is a key parameter to analyze the thermal insulation performance. How_ in (4) may fluctuate with time due to the variaever, the quantity m tion of the amount of heat flowing the walls and the wetted surface S. It is more appropriate here to use as boil-off parameter the mass loss ratio (or boil-off ratio), denoted by RBO :

Rt RBO ¼

0

_ 0 mdt M

ð5Þ

The numerator in (5) corresponds to the mass of liquid vaporized at time t while M refers to the fuel initial mass. Remember that cryogenic tanks for space launchers are namely designed for short lifetimes. As the fuel is consumed quickly (in the order of minutes), they furthermore require the propellant to be stored and maintained only for a relatively short period. In such application, the acceptable boil-off rate is about 1–2% of the initial mass per hour [13]. 3.3. External boundary conditions As evoked in Section 2, the incoming flux, Qext, depends on the external conditions that are ground and flight conditions. Let us denote by tpad the waiting time between the end of fueling and the launch.  For 0 < t 6 tpad during which the rocket is on the launch-pad for final inspections, natural convection and radiation are the main heat exchange with surrounding medium. Under the linearized radiation formulation, the overall heat flux can be expressed as:

Q ext ¼ hext ðT ext  TÞ at x ¼ es þ ei

ð6Þ

where Text is the surrounding environment temperature; T is the external insulation surface temperature; and hext is the overall exchange coefficient between the surrounding medium and the external surface taking into account the convection and the radiation contributions as detailed in Section 4.4.

4867

 For t P tpad where the rocket undergoes ascent followed by flight phases, different heat transfer mechanisms can take place at the external insulating surface. Among them, an aerothermodynamical heating generated by the friction of the insulation free surface with flowing air, the absorption of radiation from surrounding light sources (atmosphere, stars, sun...), and the emission of radiation by the warmed surface toward the surrounding environment. A detailed analysis of these heat sources (as function of the flight velocity, body orientation, and the protective layer properties for example) could enable us an altitude dependent formulation of Qext [14]; however such analysis is complicated and is beyond the scope of this study. The incoming flux Qext will be considered here constant during the flight as specified in Section 5.1. When the insulating layer does not reach the polymer degradation temperature, namely Tdeg, the energy balance at the external boundary writes:

Q ext ¼ k

@T @x

at x ¼ es þ ei

and for T < T deg

ð7Þ

Depending on magnitude of Qext, the flight duration, and the thermophysical properties of the thermal insulation, the temperature of external surface may increase and reach the degradation temperature of the polymer Tdeg provoking polymer melting or burning and finally ablation processes due to the shear force of flowing air. Tdeg can be the melting temperature if the polymer is crystalline or the decomposition (or combustion) temperature if the polymer is amorphous. A part of Qext is hence released to melt or burn the polymer through an endothermic transformation. The energy balance at the external boundary in presence of ablation according to Fig. 2 becomes [15,16]:

Q ext ¼ k

@T @x þ fs qs Ldeg @x @t

at x
ð8Þ

The product fsqs corresponds to the true polymer foam density (i.e. without accounting for the enclosed fluid density) in which fs and qs are the polymer volume fraction in the foam and the bulk polymer density, respectively; Ldeg denotes the latent heat of degradation of the polymer; and ox/ot refers to the velocity of the degradation front. Strictly speaking, ablation processes lead to Stefan type moving boundary problems, which are numerically very challenging [5,6,16]. In addition, the degradation front velocity ox/ot is not known a priori avoiding the use of the boundary condition (8) as in classical ablation problems [15]. In consequence, (1) should be solved at each time step with the boundary conditions (2) and (7). In this case, each time where a thin part of the insulation reaches the degradation temperature, it is removed by moving the external boundary at abscise at which the insulation temperature equals the degradation temperature [5,6,16]. It is clear that such method requires very fine time and space discretizations and therefore excessive computation time. In this work, we have chosen a simpler approach by reasoning continuously on the initial domain [0; es + ei] and by monitoring a ‘‘reaction front’’ according to the released energy rate associated with the degradation processes. Such an approach is inspired from the enthalpy method initially developed to solve moving boundary problems due to phase change either melting or solidification [17,18]. In this framework, the external flux condition (7) is always applied to the initial surface position, but as far as thermophysical properties are concerned, the thermal insulation undergoing ablation is treated as a bilayer material (see fig. 2). In doing so, we assign to the ablated region (‘‘reacted phase’’) a zero heat capacity and an infinite thermal conductivity according to:

4868

CðTÞ ¼

J. Randrianalisoa et al. / International Journal of Heat and Mass Transfer 54 (2011) 4864–4874

8 > C foam > > < > > > :

kðT; qÞ ¼

Ldeg fs qs DT deg

0 

for T < T deg  for T deg 

DT deg 2

DT deg 2

for T P T deg þ kfoam

for T < T deg

1

for T P T deg

6 T < T deg þ

DT deg 2

ð9Þ

DT deg 2

ð10Þ

The subscript ‘‘foam’’ refers to the properties of the solid foam material; DTdeg is the temperature interval (of mid temperature Tdeg) within which the degradation (melting or combustion) occurs. In such way, the ‘‘reacted matter’’ is considered as totally removed from the insulating layer, which is not unreasonable seeing the drastic shear stresses imposed by the flowing air. In the computation viewpoint, the problem consists in solving (1) subjected to boundary conditions (2) and (7) and using the specific heat and thermal conductivity given by (9) and (10). This can be performed without difficulty through classical numerical methods [19]. In this work, the Comsol FE solver [20] is chosen thanks to the simplicity of implementation.

compared to that of solid matrix), its effects on the foam thermal conductivity can not be neglected. However, if the confined gas is cooled at very low temperatures, it can liquefy or solidify and hence affects both foam thermal conductivity and specific heat. After storage during several tens of days, the blowing agent enclosed in cells diffuse out and is replaced by air at atmospheric pressure [25,26]. Due to our difficulty to retrieve air properties over the current temperature range and different pressures, we consider that nitrogen is the fluid enclosed in cells, which is a reasonable approximation. In addition, we consider that the confined nitrogen undergoes an isochoric transformation (due to cooling or heating from room temperature) neglecting thus the contraction or expansion of cells. According to the NIST database [21], the boiling point of nitrogen under isochoric cooling from room temperature is between 63 K and 65 K (which is lower than the boiling temperature of nitrogen about 77 K at atmospheric pressure). Below these temperatures, the nitrogen liquefies and may solidify if the temperature decreases continuously. The volumetric specific heat of foams, namely qCfoam, can be given by the following mixture rule:

qC foam ¼ ð1  fs Þqf C f þ fs qs C s 4. Thermophysical properties

ð11Þ

We focus only our attention to LH2 fuel. Hydrogen thermophysical properties from 20 K to 100 K at a constant (functioning) pressure of about 2.9 bars are taken from the National Institute of Standard and Technology (NIST) database [21]. Note that at the vaporization temperature (Tvap = 24.5 K at 2.9 bars), the hydrogen properties present an abrupt change due to the difference between liquid and gaseous hydrogen properties.

fs is equal to the foam relative density q as long as the enclosed fluid is in vapor state. The subscripts foam, f and s refer to foam, enclosed fluid (nitrogen here) and solid (PVC polymer here), respectively. The specific heat of amorphous PVC polymer, Cs, in the temperature range 20–300 K is available from the Advanced Thermal Analysis (ATHAS) database [27] while those above 300 K are measurements from literature [28]. The polymer density, qs, is taken equal to 1400 kg/m3. The density, qf, and specific heat of nitrogen, Cf, involved in (11) below the vaporization temperature 64 K are from literature [29] while those for higher temperatures are again from NIST database [21]. It is interesting to note that there are three magnitude orders between densities of liquid and vapor nitrogen. The overall thermal conductivity of closed-cell foams may be given by the following relations, including the contribution of (i) conduction within solid matrix (by lattice vibration) and gaseous substance (by gas molecule interactions); and (ii) radiation due to the thermal radiation exchanged between cell walls and through them:

4.2. The metallic shell

kfoam ¼ kcond þ krad

The thermophysical properties of the fuel, metallic shell, and insulating layer are all necessary. Their temperature dependence within the current temperature range is non-negligible and should be taken into account for realistic analyzes. Data used in the present calculations are essentially taken from literature as specified in the following. 4.1. The fuel

The metallic shell is made of aluminum alloy 2219T87 with a density of about 2800 kg/m3. The minimum thickness of the wall es is about 2 mm required to sustain the reservoir pressure and all the flight loads. This alloy is chosen thanks to its interesting characteristics for aerospace applications such as good weldability, corrosion resistance, and mechanical properties over a wide temperature range [22]. The specific heat and thermal conductivity of aluminum 2219T87 as function of temperature are obtained from reference [23]. 4.3. The insulating layer The choice of the cryogenic insulating material is constrained by several conditions such as light-weight, low thermal diffusivity, good resistance to permeability, and high degradation temperature. Among existing insulating material, low density PVC closedcells foams fulfil reasonably these requirements [24]. Such foams are characterized by relative density inferior to 0.3 [10] and can be delivered in form of rectangular panels of size about 1 m  1 m. The thermophysical properties of foam materials depend on the thermal properties and the relative density of constituting phases. Whereas the effects of enclosed gas on the foam specific heat is usually neglected (due to the low density of gas

ð12Þ

The solid and gas conductivities of foams kcond have been widely studied. These investigations have shown that the derivation of a unique conductivity model for porous materials including foams is very difficult because the heat transport is strongly affected by various factors (impurities, defects, cell anisotropy, foam tortuousity, foam density. . .). That explain why there exists several thermal conductivity models of porous materials derived from different approaches (field, resistor, phase averaging approaches. . .) [30]. Comparison between results of most relevant models and data of PVC foams [31] at room temperature shows that the best compromise over the low density range 0 < q⁄ < 0.35, is obtained from the Glicksman formula [32] below.

kcond ¼ ð1  fs Þkf þ

ð2  /Þfs ks 3

ð13Þ

where / is the fraction of solid contained in cell edges; and ks and kf are respectively the solid and fluid phase conductivities. We have checked that these last can be taken as the bulk polymer and continuum nitrogen conductivities showing that no size effect can be expected over the current temperature and gas pressure ranges. The temperature dependent thermal conductivity of PVC polymer are taken from [33] for T < 300 K and from [34] for T > 300 K. The isochoric thermal conductivity of nitrogen are retrieved from [35] for T < 64 K and from NIST database [21] for T > 64 K.

4869

J. Randrianalisoa et al. / International Journal of Heat and Mass Transfer 54 (2011) 4864–4874

The formulation of radiation within the foam analogous to Fourier’s law, which leads (12), requires that the radiative transfer is in the diffusion regime. For foam layers having thicknesses much greater than few millimeters and operate in a temperature range inferior to 600 K, it can be shown that the radiation diffusion regime prevails. In this case, the radiative conductivity krad can be expressed through the famous Rosseland diffusion relationship [36]:

K rad ¼

16 rT 3 3K 

ð14Þ

With r the Stefan–Boltzmann constant, T the temperature at which the conductivity is evaluated, and K⁄ the extinction coefficient of the foam. For closed-cell foams, refined calculation of K⁄ based on details of the foam microstructure can be performed as in [37] but, here, the simple expression proposed by Campo-Arnaiz suitable for thin cell walls is adopted according to K  ¼ 4:10ð/fs =DÞ1=2 þ ð1  /Þfs K s [38]. In this relation, D is the average cell size and Ks is the temperature dependent extinction coefficient of the polymer. For classical PVC foams, we have /  0.3 and D  200 lm. Ks can be approximated by the following quadratic increasing function: Ks = 0.0394T2  4.840T + 648.74 for 70 K < T < 500 K. Note that (i) with a typical wall thickness, namely ew, about 10 lm, the thin wall condition is fulfilled since Ksew  1 over the current temperature range; (ii) the radiation contribution is only expected to be significant at high temperature and for very low density foams; and (iii) the convection within cells is neglected because the cell size is generally smaller than few millimeters [1,10]. The degradation properties required in the current model are Tdeg, DTdeg, and Ldeg. They are intrinsic properties of the polymer constituting the foam and can be determined by calorimetric measurements on bulk polymer. The degradation kinetic analyzes of PVC have shown that the main degradation mechanism is the thermal decomposition and combustion [39,40]. The evidence for this degradation comes first from thermogravimetric (TGA) experiments where large decreases in mass occur beginning around 550 K [41]. The results of Dynamic Scanning Calorimetric (DSC) experiments present an endothermic peak around this temperature presumably because of the sample decomposition since it is severally charred. The main information about the degradation can be retrieved from this peak. For example, the degradation temperature corresponds to the temperature of the endothermic peak; the temperature interval corresponds to the peak width; and the energy of degradation can be given by the area under the peak. Based on ATHAS database [27], Tdeg and Ldeg of amorphous PVC are 546 K and 11 kJ/mol, respectively. DTdeg is about 46 K deduced from the DSC curves from [41]. 4.4. Heat transfer coefficients at internal and external wall surfaces Remember that the heat transfer between the surrounding medium and the external reservoir surface during the ground phase is modeled through the overall exchange coefficient hext,  NC , and ambient which includes here natural convection, namely h radiation, namely hrad . Also, the heat transfer between the internal metallic wall and the fuel being by natural convection characterized by the exchange coefficient hint. The radiation contribution to hext according to the linearized radiation flux formulation between the external wall and the surrounding medium writes [42]:

hrad ¼ erðT þ T ext ÞðT 2 þ T 2ext Þ

ð15Þ

Eq. (15) seems explicitly temperature dependent; however, it can be taken as constant because the wall temperature varies slightly with time before the launch. In this case, a typical value of hrad is estimated equal to 5 W/m2/K.

Convection transfer problems both inside and outside a reservoirs have attracted many attentions. It is now acknowledged that the heat exchange between the fluid and the walls occurs within a thermal boundary layer. The parameters characterizing such convection problems are the Grashof (Gr), Prandtl (Pr), and Rayleigh (Ra) numbers. For a vertical wall wetted up to height H, these numbers are defined by [42]:

Gr ¼

gbðT  T 1 Þ

m2

H3 ;

Pr ¼

m ; and Ra ¼ Gr  Pr a

ð16Þ

With g the acceleration due to gravity, b the fluid volumetric expansion coefficient, T the wall temperature, T1 the bulk fluid temperature (outside, T1 is the temperature Text while inside, it is the unperturbed fuel temperature Tint), m the fluid cinematic viscosity, a the fluid diffusivity. All the fluid properties involved in (16) are evaluated at the mean temperature Tm [42]. The value of Ra number determines generally the type of the flow regime. In the present study, estimation of Ra values both inside and outside shows that the flows are strongly turbulent (about 10+13 for external flow and 10+16 for internal flow). For external flow, thanks to the large curvature of the reservoir, the convection contribution to hext when the rocket is on the launch-pad can follow the well-known Chu and Churchill empiric law of convection on flat vertical walls [24,42,43]:

8 >  hNC H < Nu ¼ ¼ 0:86 þ h > kf :

0:367Ra

1=6

1 þ ð0:492=PrÞ9=16

92 > =

i8=27 > ;

ð17Þ

where Nu is the well-known Nusselt number averaged over the reservoir height H; kf is the fluid thermal conductivity. Based on Ra value of 7.86  10+13 and the Pr value of 0.72 at a mean air tem NC  6:04 W=m2 =K. It is perature Tm of 250 K, relation (17) predicts h interesting to note that the free convection and the ambient radiation are of the same magnitude order. Therefore, both contributions need to be considered in the boundary conditions during the ground phase. Finally, the overall exchange coefficient hext is about 11 W/m2/K. For the internal flow, single and two-phase natural flow can take place depending on the wall excess temperature defined by DTvap = T  Tvap. When DTvap < 0, single liquid phase flow prevails. In this case, the Ra number can reach a value of 7.6  10+16, which is quite classical in cryogenic tanks. It indicates that the internal flow is strongly turbulent. To the knowledge of the authors, none correlation of exchange coefficient in vertical cylinder with this high Ra number is available in literature. An alternative approach is adopted here by considering the fact that the reservoir size is extremely large (5 m diameter here) and is much larger than the boundary layer thicknesses. Thus, the boundary layer development can be assumed as if it occurs at an isolated vertical plate in an infinite medium. Based on this, we use again the above Chu and Churchill formula (17), which is generally suitable for any Ra value. It  NC  272 W=m2 =K for Ra = 7.6  10+16 and Pr = 1.32 at a predicts h mean fluid temperature of 22 K. Note however that this value of  NC corresponds to ground phase conditions during which the fuel h is not consumed. In real situation, the fuel quantity decreases during the flight as the fuel is consumed (not only by the cryogenic engine but also by the pressurization system as evoked previously) causing a decrease of wetted height H in (16). A simple estimation considering a typical fuel outflow of 0.615 m3/s shows that the exchange coefficient from (17) during the ground phase varies weakly with H. More precisely, it increases just of an amount of 2 W/m2/K for a 10 min flight time. When the wall temperature excesses the LH2 vaporization temperature (DTvap > 0), the heat is first removed by natural convection flow to the free surface and then by evaporation to the

4870

J. Randrianalisoa et al. / International Journal of Heat and Mass Transfer 54 (2011) 4864–4874

vapor space when the excess temperature is small (DTvap < 1 K for LH2 [44] depending on the vaporization pressure). The usual formulas of free convection such as (17) here can be applied but with fluid properties evaluated at Tvap instead of at Tm [42]. Now when DTvap excesses a certain value, nucleate then transition and film boiling regimes may take place [45]. For design purpose of cryogenic insulation, the amount of vaporized fuel should be as low as possible. It is thus sufficient to limit our analysis to the case where the wall excess temperature remains reasonable (i.e. allowing nucleate boiling regime). In this work, we adopt the Kutateladze nucleate pool boiling model thanks to its reasonable accuracy for cryogenic fluids whatever the geometry [45]. In the presence of subcooled liquid1 in the nucleate boiling regime, the Kutateladze correlation (indexed by NB) can be given by:

( !)1=3  2 2 P ½ql  qv 9 kl Lv  NB ¼ F h 3:37  10 sub  q2v rg l Cl

ð18Þ

The terms in bracket in (18) are the Kutateladze correlation for saturated liquid from reference [46]. g is the acceleration gravity; q is the heat flux crossing the walls; P is the pressure; q is the density; C is the specific heat; and r is the surface tension between the liquid and its own vapor. The subscripts l and v refer respectively to the liquid properties (evaluated at Tvap) and vapor properties (evaluated at the film mean temperature (T + Tvap)/2). These properties can be available from NIST data bank [21]. Fsub is the Sakurai correction factor [47] given by (19) and aiming to account for the fact that the fuel can be subcooled (i.e. the bulk temperature remains smaller than its saturation temperature with a temperature difference DTsub = Tvap  Tint > 0 although boiling occurs).

 0:65  1:5 q C f DT sub F sub ¼ 1 þ 0:87 l qv Lv

ð19Þ

In a log–log plot, the overall heat transfer coefficient hint varies with the wall excess temperature DTvap according to a hyperbolic function characterized by two asymptotes given by (17) and (18). At isobaric pressure 2.9 bars, it can be shown that free convection regime dominates the heat exchange when DTvap < 0.2 K while nucleate boiling takes over when 0.2 K < DTvap 6 2.8 K. The value DTvap = 2.8 K corresponds to the critical (maximum) heat flux for  NB may reach 7  10+4 W/m2/K. In fact, above this temperwhich h ature value, transition boiling starts to take place with a heat flux decreasing with the increases of excess wall temperature.

5. Results and discussion 5.1. Preliminary analysis concerning the incoming flux According to our simplifying hypothesis concerning the incoming external flux Qext, we adopt the following simple form which depends on time:

 Q ext ¼

1 Because the bulk fluid temperature remains smaller than its vaporization temperature.

3500

ð20Þ

for t pad < t 6 t pad þ Dt

with t = 0 the initial time following the fueling step; tpad = 5 min the ground phase duration; and Dt = 10 min the usual lifetime of the main cryogenic reservoir of launchers. The value of 3.5 kW/m2 in (20) is our estimation of the average external flux over the flight duration. Indeed, the external flux increases from the ground flux (<1 kW/m2) to a peak flux about 7 kW/m2 [48] during the launcher passage through the atmosphere. First of all, it seems necessary to study the sensitivity of the model to the external flux Qext, which is here the most uncertain input parameter. To do this, let us consider an insulation having a thickness of 2 cm and a density of 70 kg/m3 (or a relative density of q⁄ = 0.05), which are usual for cryogenic launcher insulation. For the cryogenic engine, a typical fuel outflow of 0.615 m3/s is considered throughout this study. It is required in order to evaluate the reduction of the reservoir wetted surface S with time [see Eq. (4)]. In Fig. 3 is plotted the boil-off ratio RBO at time t = tpad + Dt as function of Qext (bottom axis) or the relative deviation on Qext (top axis). For a given insulation system, the boil-off ratio RBO increases with Qext due to the increase of the incoming heat flux (more energy is available to vaporize the fuel). Based on the current calculation parameters, the boil-off ratio with Qext = 3.5 kW/ m2 and 15 min service time is 0.3% of the initial fuel mass or 1.2% per hour in term of boil-off rate. It falls into the admissible boil-off rate range 1–2% per hour [13] mentioned previously. This result tends to confirm that the Qext value of 3.5 kW/m2 taken as reference here seems a good estimate of the external flux over the current flight duration. 5.2. Analysis of heat transfer through the walls, fuel temperature, and boil-off For the sake of clarity, in the following, we will denote

q⁄ = (qfoam/qs) the foam relative density. In Fig. 4a is depicted the

4.5. Remarks

0.75

Mass loss ratio RBO, %

(i) For LH2 at 2.9 bars and an excess temperature less than 2.8 K, the soobcooling scaling factor Fsub is about 1.2. (ii) At the highest altitude of usual launchers as far as the cryogenic tank is concerned, the gravity acceleration (about 9.23 m/s2 at 200 km above sea level) is not much changed compared to the ground reference value 9.81 m/s2. Therefore, the reduction of gravity acceleration on fluid flows can be also neglected. (iii) The thermophysical properties described in this section are introduced in the Comsol program through polynomial functions, which best fit them. Moreover, piecewise functions are used whenever the properties present abrupt variations (e.g. for the hydrogen and nitrogen thermophysical properties, and the foam properties given by (9) and (10)).

hext ðT  T ext Þ for 0 6 t 6 t pad

-30

Uncertainty on Qext, % -10 0 10

-20

20

30

0.50

0.25

0.00 2.5

3.0

3.5 4.0 2 External flux Qext, kW/m

4.5

Fig. 3. Evolution of the total mass loss ratio as function of the external flux (bottom axis) or its uncertainty (top axis).

4871

J. Randrianalisoa et al. / International Journal of Heat and Mass Transfer 54 (2011) 4864–4874

Temperature field T(x), K

500 400 300 200 100 0 0.0

ρ* = 0.01 0 min 5 min 6 min 7 min 8 min 9 min 15 min ρ* = 0.05 0 min 5 min 6 min 7 min 8 min 9 min 15 min

0.4 0.8 1.2 1.6 distance from internal boundary x,cm

2.0

2.2

(a) ρ* = 0.01 0 min 5 min 6 min 7 min 8 min 9 min 15 min

35 30

ρ* = 0.05 0 min 5 min 6 min 7 min 8 min 9 min 15 min

25 20

0.0

ρ* = 0.01 0 min 5 min 6 min 7 min 8 min 9 min 15 min

540

0.2 0.4 0.6 distance from internal boundary x, cm

Temperature field T(x), K

Temperature field T(x), K

40

535 530

ρ* = 0.05 0 min 5 min 6 min 7 min 8 min 9 min 15 min

525 520 515

0.8

0.8

1.0 1.2 1.4 1.6 1.8 2.0 distance from internal boundary x, cm

2.2

(c)

(b)

Fig. 4. (a) Temperature fields at different instants after the fueling (t = 0 min) to the end of flight (t = 15 min). Solid curves with gray filled symbols: uniform insulation density q⁄ of 0.01; and solid curves with white filled symbols: uniform insulation density q⁄ of 0.05; (b) Zoom of the temperature fields near the internal boundary; (c) Zoom of the temperature fields near the external boundary during flight (t > 5 min).

temperature profiles across the tank walls at different instants starting after the fueling step to the reservoir releasing. The insulation system has uniform density: the solid curve and gray filled symbols are for an insulation of density q⁄ = 0.01; and the dash curve and white filled symbols are for an insulation of density q⁄ = 0.05. The thick solid and dash curves refer to initial temperature fields (t = 0 min) while the symbols correspond to temperature fields at different instants from the launch (t = 5 min) to end of the flight (t = 15 min) as far as the reservoir is concerned. The region 0.0 cm < x < 0.2 cm corresponds to the highly heat conductor aluminum component, which explains why the temperature in this region is quasi constant. The remaining region corresponds to the insulating layer and presents three distinct temperature slopes. (i) The small temperature slopes below 64 K correspond to the region of the insulation where the nitrogen enclosed within foam cells is in liquid state. The liquid nitrogen is denser and more conductor than the polymer. (ii) The large temperature gradients between 64 K and 523 K correspond to the region where the polymer is yet undamaged and the nitrogen is in gaseous state. Since the polymer conductivity is greater than the gas conductivity, the effects of the foam density on the temperature fields are more noticeable. And (iii) the region with quasi-constant temperatures above 523 K at instant t corresponds to the reacted phase of the insulation according to our simple ablation model although in reality, this melted (or burned) layer is instantaneously removed

by the flowing air. Fig. 4b gives details of temperature fields close to the internal boundary (x < 0.8 cm). The temperatures close to internal boundary increase with time up to 25 K, which is approximately the ebullition temperature of the LH2. The effect of the foam density is negligible during ground waiting phase while it is more noticeable during the flight phase because of the difference of the propagation of degradation fronts between low and high insulation densities (the position of the degradation front is given here by the abscissa of the inflection point at 523 K). Fig. 4c shows a zoom of the temperature fields close to the external boundary (x > 0.8 cm). It can be seen that the temperatures in the insulation remain inferior to Tdeg = 546 K because of degradation when T > Tdeg  DTdeg = 523 K (horizontal dash line) followed by the ablation processes. The advancement of the degradation fronts with time using data in Fig. 4c are shown in Fig. 5 by solid and dash curves. We can note that: (i) the beginning time of ablation depends on the insulation density. The lower the density, the quicker the initiation of ablation is. (ii) In a first time, the degradation front progresses quasi-linearly with time indicating that the front velocities are approximately constant. The lower the insulation density, the faster the degradation is. The front velocities2 are respectively around

2

Defined here by the slopes dx/dt with x the front position at the time t.

4872

J. Randrianalisoa et al. / International Journal of Heat and Mass Transfer 54 (2011) 4864–4874

Flight time, min -5. 0

-2. 5

0. 0

2 .5

5. 0

7 .5

-5.0

10. 0

25 Maximum fuel temperatures, K

Position of the degradation front, cm

2.2 2.0 1.8

dx/dt = 42 μm/s

1.6

dx/dt = 25 μm/s

ρ* = 0.01 ρ* = 0.05 linear law expo. law piecewise constant law

1.4 1.2

-2.5

Flight time, min 2.5 5.0

7.5

10.0

Vaporization temperature

24 23 22

ρ* = 0.01 ρ* = 0.05

21 20

1.0 0.0

2.5

5.0

7.5 10.0 time t, min

12.5

0.0

15.0

Fig. 5. Insulation thickness drop with time due to ablation for different insulating layer densities. Solid curve: uniform density q⁄ of 0.01; dash curve: uniform density q⁄ of 0.05; dot lines with black triangles: piecewise constant density; dot lines with white filled circles: linear density; and dot lines with gray squares: exponential density.

2.5

This is actually the metallic wall temperature at the abscissa x = 0 cm.

5.0

7.5 time t, min

10.0

12.5

15.0

Fig. 6. Instantaneous evolution of the maximum fuel temperatures for different insulating layer densities. Straight curve: uniform density q⁄ of 0.01; and dash curve: uniform density q⁄ of 0.05.

1.0

1.0

ρ* const = 0.01 ρ* const = 0.03 ρ* const = 0.05 linear law exponential law piecewise constant law

0.8

Boil-off ratio RBO, %

42 lm/s and 25 lm/s for the insulating layers of densities q⁄ = 0.01 and 0.05. (iii) After a certain time, the insulation layer damage is attenuated and then stops. In fact, when the functioning insulation becomes thin (about 1 cm here), its thermal resistance decreases drastically. As consequence, most of the flux Qext penetrates the walls and the rest of this heat flux is too small to provoke ablation. The evolutions of the maximum fuel temperature3 (hereafter referred to as MFT) and the boil-off ratio RBO with time for the two cases of uniform insulation densities are respectively shown in Figs. 6 and 7 by the solid (q⁄ = 0.01) and dash curves (q⁄ = 0.05). After the fueling, the fuel temperature is free of control whereas the reservoir is exposed to the ground thermal aggression (convection and radiation). During the first instants, the heat flux transferred into the fuel is little because the temperature difference between the fuel and the internal wall is small, and exchange coefficient is not very large. There is consequently energy accumulation within the walls. This is represented by the rapid increase of MFTs at short times. When the time goes on, the heat dissipates more easily into the fluid (because the temperature difference between the wall and bulk liquid increases). This results in a slow rise of MFTs. The fuel heating begins to be important few minutes after the launch and at the same time, the boil-off appears. In fact, the heat flux crossing the walls becomes intense due to the aerothermal heating but this heat is first dissipated slowly to the fuel because of the natural fluid circulation. This causes again an abrupt temperature elevation of the tank wall near the fuel/wall interface, thus the MFTs. The fuel boil-off occurs when the MFT excesses its vaporization temperature (indicated by gray line in Fig. 6). In a first time, i.e. during the short instant where the MFTs continue to increase rapidly, the boil-off is essentially due to fuel evaporation at free surface resulting from natural circulation of saturated fuel. After that, the nucleate boiling removes rapidly the heat from the internal wall provoking, on one hand, a rapid increase of boil-off ratio and, on the other hand, a slow down of the fuel heating. At long times, the heat transfer reaches the steady-state regime, which explains why the MFTs become quasi-constant. However, at the same time, the boil-off ratios RBO continues to increase _ rebut in few extents because the instantaneous mass loss rate m duces with time proportionally to the wetted surface S (in Eq. (4)).

3

0.0

0.6

0.8 0.6

0.4

0.4

0.2

0.2

0.0

0.0 5

6

7

8

9

10 11 time t, min

12

13

14

15

Fig. 7. Instantaneous evolution of the boil-off ratios for different insulating layer densities. Straight curve: uniform density q⁄ of 0.01; dash-dot curve: uniform density q⁄ of 0.03; dash curve: uniform density q⁄ of 0.05; dot lines with black triangles: piecewise constant density; dot lines with white filled circles: linear density; and dot lines with gray filled squares: exponential density.

The efficiency of each insulation is observed through the comparison of boil-off ratios RBO corresponding to different insulation densities. As far as uniform insulation densities are concerned, one can see that the lower the density, the higher the boil-off ratio is. With an insulation of density q⁄ = 0.05, the boil-off remains acceptable. However, when the density is decreased to 0.01, the boil-off ratio reaches quickly the value limit of 0.5% of initial mass (or 2% of initial mass per hour in term of boil-off rate) after 6 min flight. This last case should be prevented because it will decrease a lot the quantity of fuel and may provoke the launcher disfunctioning. It is very interesting to note that, if the boiling mechanism is neglected in the current modeling, i.e. only natural convection is accounted for, the maximum fuel temperature continues to increase during a long time and the boil-off ratio becomes underestimated. Moreover, the insulating layer degradation is falsified. This shows the crucial importance of capturing both single (free convection) and two-phase (free convection and nucleate boiling here) natures of internal flow.

J. Randrianalisoa et al. / International Journal of Heat and Mass Transfer 54 (2011) 4864–4874

Considering now both ground and flight phases, it is clear that the lowest insulation density (q⁄ = 0.01 here) is not really the best thermal insulation of cryogenic rockets due to its rapid degradation provoking rapid and important fuel boil-off. For satisfying the overall flight time with this type of insulation, either higher density (as q⁄ = 0.05 here) insulation or light and thick insulation should be chosen but in any way at the expanse of heavy structure. So far, the results with insulation systems having uniform density are discussed. To overcome the disadvantages for using very low density foams, the insulating foam can be recovered with a layer of thermal protection system (TPS). Note that the use of TPS is classical for reusable space vehicles due to the high aerothermal flux (greater than 100 kW/m2) during the re-entry phase [48], it is not however the best solution to meet the low weight and cost requirements for single-use cryogenic tanks as considered here. An alternative and much economical solution may be the use of insulating foam having graded density, i.e. low density sub-layers (bonded on the metallic shell surface) to ensure efficiently the thermal insulation and high density upper-layers to retard the degradation. This can be made, either by bonding various foam layers of different densities or directly, during the foam manufacturing by tuning for example the types of blowing agent and/or the processing temperatures. To investigate the improvements of using insulation with graded density, let us consider three cases of density profiles following respectively piecewise constant (21), linear (22) and exponential laws (23). In each case: the insulation thickness is fixed to ei = 2 cm; the highest relative densities are fixed to 0.05 (they are the relative densities of layers close to the external surface); and the average relative density (over the insulation thickness), denoted herein by hq⁄i, is fixed to 0.03. Therefore, the parameters involved in (21)–(23) with x in cm are respectively a1 = 0.01, b1 = 0.05; a2 = 0.02, b2 = 0.01; and a3 = 0.017, b3 = 0.539.

q ¼



a1

if 0 6 x < ei =2

b1

if ei =2 6 x < ei

ð21Þ

q ¼ a2 x þ b2 with 0 6 x 6 ei

ð22Þ

q ¼ a3 eb3 x with 0 6 x 6 ei

ð23Þ

The positions of the degradation front according to time corresponding to three graded density insulations are also shown in Fig. 5 by the dot curves with symbols. When t < 10 min, the insulation with piecewise constant density resists better to ablation because the dense upper-layer (local relative density of 5%) is not entirely ablated until this time. However when t > 11 min, the three insulations tend to ablate rapidly and similarly because the unablated layers have all low relative densities about 0.01. The evolution of corresponding boil-off ratios at different instants are also depicted in Fig. 7. For comparison, the result for an insulation with uniform density q⁄ of 0.03 is shown by dashdot curve. The lower and upper bounds of boil-off ratios correspond respectively to the cases of densest (q⁄ = 0.05) and lightest uniform density insulating layers (q⁄ = 0.01). It can be seen that for an insulation having identical average density (e.g. hq⁄i = 0.03 here), the graded density layers provide lower and satisfactory mass loss ratio (<0.5% of the initial mass). More particularly, the piecewise constant density variation seems better among the three graded density profiles considered here. It is clear from this analysis that the denser layer is the best solution as far as the boil-off constraint is concerned. However, in applications where the light-weight structures are required such as in space transports, the insulations with graded densities are most promising. Up the stage of this work, it is not obvious to identify the optimal density profile for the insulating layer. In fact, one

4873

should consider additional constraints such as the feasibility, mechanical reliability constrains, different flight durations. 6. Conclusion In this paper, a simple model for predicting both fuel boil-off ratio and thermal protection ablation of a cryogenic reservoir subjected to ground and flight external flux is presented. The following conclusions can be drawn: To get a realistic estimation of the maximum fuel temperature and fuel boil-off, a detailed modeling of the internal heat flow, which should include both single and two-phase flows, is necessary. It can be performed through exchange coefficients as shown here. With a typical insulating layer characterized by a 2 cm thick low-density PVC foam, it is shown that during the usual ground phase, the maximum fuel temperature remains well below the vaporization temperature, as far as heating by the surrounding environment through natural convection and radiation is considered. However, few times after the launch, the maximum fuel temperature increases rapidly above the vaporization temperature. Consequently, the flight conditions such as the flight duration and the magnitude of aerothermal flux rather than the ground conditions influence mainly the choice of insulation characteristics. The aerothermal flux during ascent phase provokes the insulation degradation then its ablation. While the insulation continues to act, the ablation rate is almost constant; otherwise, the ablation stops but the boil-off increases progressively with time. It is shown that using polymer foam as insulation, the lower the insulation density, higher the ablation velocity is. To satisfy the flight duration and the boil-off requirements, insulation capable to resist to ablation is required. Based on usual uniform density insulation, it leads to dense (and hence heavy) insulating layer. It is however demonstrated here that using insulating layers with spatially graded densities; it is possible to design light and anti-ablative insulation. Among the three graded-type density profiles having 3% average relative density investigated in this paper, the insulation with piecewise constant density appears the best solution. Indeed, it is shown to compete well against the usual cryogenic launcher insulation (i.e. with uniform density of 5%). The current study shows that there is probably optimum lightweight insulation either of uniform or graded density but its choice is not evident and requires additional investigations. Such further analyzes should take into account not only the thermal constraints, as considered here, but also other requirements such as the feasibility, the mechanical reliability of the metallic shell/insulating panels. For example, the modeling tool recently developed in reference [49], can be used for such reliability analysis. Thanks to its simplicity, the heat transfer model proposed here is without doubt straightforward to be implemented in an optimization code of cryogenic tank protection. Finally, so far, constant external heat flux is only considered. It is however interesting to note that if this parameter is known locally or over the area of an insulating panel, the current model can predict both local evolution of maximum fuel temperatures and the insulation ablation. In this case, the insulation dimensioning may lead to locally varied insulation thickness, which probably offers better weight optimization of the reservoir insulation. Acknowledgements The authors acknowledge gratefully the support of the French National Agency Research through the Project No. ANR-08MAPR-0009.

4874

J. Randrianalisoa et al. / International Journal of Heat and Mass Transfer 54 (2011) 4864–4874

References [1] C.L. Tien, C.R. Cunnington, Cryogenic insulation heat transfer, in: T.F. Irvine Jr., J.P. Hartnett (Eds.), Advances in Heat Transfer, vol. 9, Academic Press, New York, 1973, pp. 349–417. [2] F.M. Anthony, J.Z. Colt, R.G. Helenbrook, Development and Validation of Cryogenic Foam Insulation for LH2 Subsonic Transports, NASA Contractor Report, 1981 (CR-3404). [3] S.K. Mukka, M.M. Rahman, Analysis of fluid flow and heat transfer in a liquid hydrogen storage vessel for space applications, in: AIP Conference Proceedings in Space Technology and Applications International Forum, Albuquerque, NM, vol. 699, 2004, pp. 37–44. [4] S.H. Ho, M.M. Rahman, Three-dimensional analysis for liquid hydrogen in a cryogenic storage tank with heat pipe–pump system, Cryogenics 48 (2008) 31–41. [5] Y.K. Chen, F.S. Milos, Fully implicit ablation and thermal analysis program (FIAT), J. Spacecraft Rockets 36 (1999) 475–483. [6] J.A. Dec, R.D. Braun, Ablative thermal response analysis using the finite element method, in: Proceedings of AIAA 47th Aerospace Sciences Meeting and Exhibit, Orlando, FL, AIAA-2009-0259, 2009. [7] B. Laub , E. Venkatapathy, Thermal protection system technology and facility needs for demanding future planetary missions, in: A. Wilson (Ed.), Proceedings of the International Workshop Planetary Probe Atmospheric Entry and Descent Trajectory Analysis and Science, ESA SP 544, ESA Publications Division, Noordwijk, Netherlands, 2004, pp. 239–247. [8] W. Lin, S. Armfield, Long-term behaviour of cooling fluid in a vertical cylinder, Int. J. Heat Mass Transfer 48 (2005) 53–66. [9] I. Rodriguez, J. Castro, C.D. Pérez-Segarra, A. Oliva, Unsteady numerical simulation of the cooling process of vertical storage tanks under laminar natural convection, Int. J. Therm. Sci. 48 (2009) 708–721. [10] L.J. Gibson, M.F. Ashby, Cellular Solids: Structure and Properties, second ed., Cambridge University Press, Cambridge, UK, 1997. [11] J.R. Sharp, Ablation modeling of Ares-I upper stage thermal protection system using thermal desktop, in: Proceedings of 18th Thermal and Fluids Analysis Workshop, Cleveland, OH, TFAWS-07-1009, 2007. [12] J. Lacapere, M. Gardette, Space engineering activities at cryospace and air liquide, Fluent News Fall (2005) S4–S7. [13] S. Mital, J. Gyekenyesi, S. Arnold, R. Sullivan, J. Manderscheid, P. Murthy, Review of current state of the art and key design issues with potential solutions for liquid hydrogen cryogenic storage tank structures for aircraft applications, in: NASA Technical Memorandum, TM-214346, 2006. [14] H. Klinkrad, B. Fritsche, A. Kashkovsky, Prediction of spacecraft destruction during uncontrolled re-entries, in: C. Stavrinidis, A. Rolfo, E. Breitbach (Eds.), Proceedings of the European Conference on Spacecraft Structures, Materials and Mechanical Testing, ESA SP 468, European Space Agency, Noordwijk, Netherlands, 2001, pp. 485–491. [15] L. Dombrovsky, L. Schunk, W. Lipinski, A. Steinfeld, An ablation model for the thermal decomposition of porous zinc oxide layer heated by concentrated solar radiation, Int. J. Heat Mass Transfer 52 (2009) 2444–2452. [16] R.E. Hogan, B.F. Blackwell, R.J. Cochran, Application of moving grid control volume finite element method to ablation problems, J. Thermophys. Heat Transfer 10 (1996) 312–319. [17] A. Esen, S. Kutluay, A numerical solution of the Stefan problem with a Neumann-type boundary condition by enthalpy method, Appl. Math. Comput. 148 (2004) 321–329. [18] S.K. Wong, A. Walton, Numerical solution of single-phase Stefan problem using a fictitious material, Numer. Heat Transfer B 35 (1999) 211–223. [19] A.R. Gourlay, Hopscotch: a fast second-order partial differential equation solver, J. Inst. Math. Appl. 6 (1970) 375–390. [20] Comsol Multiphysics 3.5a, Comsol Inc., Burlington, MA, 2008. [21] E.W. Lemmon, M.O. McLinden, D.G. Friend, Thermophysical properties of fluid systems, in: P.J. Linstrom, W.G. Mallard (Eds.), NIST Chemistry WebBook, NIST Standard Reference Database 69, National Institute of Standards and Technology, Gaithersburg, MD, 2003. [22] J.R. Davis (Ed.), Aluminum and Aluminum Alloys, ASM Specialty Handbook, ASM International, Menlo Park, CA, 1996. [23] N.J. Simon, E.S. Drexler, R.P. Reed, Review of cryogenic mechanical and for the period thermal properties of AI–Li alloys and alloy 2219, Astronautics Laboratory (AFSC) Report, 1990 (AL-TR-90-064). [24] A.J. Colozza, Hydrogen Storage for Aircraft Applications Overview, NASA Contractor Report, 2002 (CR-211867).

[25] D. Eaves, Polymer Foams Trends in Use and Technology, Rapra Technology Ltd., Shawbury, UK, 2001. [26] Ch-j. Tseng, M. Yamaguchit, T. Ohmorit, Thermal conductivity of polyurethane foams from room temperature to 20 K, Cryogenics 31 (1997) 305–312. [27] M. Pyda, (Ed.), ATHAS Data Bank, 2005. . [28] M. Kok, K. Demirelli, Y. Aydogdu, Thermophysical properties of blend of poly(vinyl chloride) with poly(isobornyl acrylate), Int. J. Sci. Technol. 3 (2008) 37–42. [29] I.N. Kudryavtsev, K.E. Nemchenko, Lattice dynamics and heat capacity of solid nitrogen, in: Proceedings of the 10th International Conference on Phonon Scattering in Condensed Matter: Phonons 2001, Hanover, NH, Po-II 54, 2001. [30] R. Singh, Thermal conduction through porous systems, in: A. Öchsner, G.E. Murch, M.J.S. de Lemos (Eds.), Cellular and Porous Materials: Thermal Properties Simulation and Prediction, Wiley-VCH, Weinheim, 2008, pp. 199– 238. [31] CES 4, The Cambridge Engineering Selector, Granta Design Ltd, Cambridge, UK, 2002. [32] L.R. Glicksman, Heat transfer in foams, in: N.C. Hilyard, A. Cunningham (Eds.), Low Density Cellular Plastics: Physical Basis of Behaviour, Chapman & Hall, London, 1994, pp. 105–152. [33] L. Risegari, M. Barucci, E. Olivieri, G. Ventura, Low Temperature Thermal Conductivity of PVC, J. Low Temp. Phys. 144 (2006) 49–59. [34] Y.S. Touloukian, R.W. Powell, C.Y. Ho, P.G. Klemens, Thermal conductivity: nonmetallic solids, in: Y.S. Touloukian (Ed.), Thermophysical Properties of Matter, vol. 2, IFI/Plenum, New York, 1970. _ [35] Y.A. Freiman, V.V. Sumarakov, A. Jezowski, P. Stachowiak, J. Mucha, Thermal conductivity of solid oxygen, nitrogen, and their solid solutions, Low Temp. Phys. 22 (1996) 194–204. [36] R. Siegel, J.R. Howell, Thermal Radiation Heat Transfer, fourth ed., Taylor and Francis, Washington, DC, 2002. [37] A. Kaemmerlen, C. Voc, F. Asllanaj, G. Jeandel, D. Baillis, Radiative properties of extruded polystyrene foams: predictive model and experimental results, J. Quant. Spectrosc. Radiat. Transfer 111 (2010) 865–877. [38] R.A. Campo-Arnaiz, M.A. Rodrıguez-Perez, B. Calvo, J.A. de Saja, Extinction coefficient of polyolefin foams, J. Polym. Sci. B: Polym. Phys. 43 (2005) 1608– 1617. [39] N.M. O’Nara, Combustion of PVC, Pure Appl. Chem. 49 (1977) 649–660. [40] A. Marcilla, M. Beltran, Thermogravimetric kinetic study of poly(vinylchloride) pyrolysis, Polym Degrad Stability (1995) 117–124. [41] A. Manivannan, M.S. Seehra, Identification and quantification of polymers in waste plastics using differential scanning calorimetry, Polym. Preprint Am. Chem. Soc. Div. Fuel Chem. 42 (1997) 1028–1032. [42] F.P. Incropera, D.P. DeWitt, Fundamentals of Heat and Mass Transfer, John Wiley & Sons, Inc., New York, 1996. [43] K.-H. Kima, H.-J. Koa, K. Kima, Y.-W. Kimb, K.-J. Cho, Analysis of heat transfer and frost layer formation on a cryogenic tank wall exposed to the humid atmospheric air, Appl. Therm. Eng. 29 (2009) 2072–2079. [44] G.V. Brown, R.H. Jansen, J.J. Trudell, High Specific Power Motors in LN2 and LH2, NASA Technical Memorandum, 2007 (TM-215002). [45] L.I. Boehman, J.J. Schauer, J. Harold, T. Stevenson, Cryogenic liquid heat transfer analysis, Aerospace Structures Information and Analysis Center (1987). ASIAC887.1D. [46] I.L. Pioro, W. Rohsenow, S.S. Doerffer, Nucleate pool-boiling heat transfer. II: assessment of prediction methods, Int. J. Heat Mass Transfer 47 (2004) 5045– 5057. [47] A. Sakurai, K. Fukuda, Mechanisms of subcooled pool boiling CHFS depending on subcooling, pressure, and test heater configurations and surface conditions in liquids, in: Proceedings of ASME International Mechanical Engineering Congress & Exposition: IMECE’02, New Orleans, LA, IMECE200234066, 2002. [48] E. Cosson, F. Deneu, O. Le Couls, V. Peypoudat, P. Baiocco, Thermal architecture of hydrogen tank on RLV second stage, in: A. Wilson (Ed.,) Proceedings of the 4th European Workshop Hot Structures and Thermal Protection Systems for Space Vehicles, ESA SP 521, European Space Agency, Paris, France, 2003, pp. 239–245. [49] J. Randrianalisoa, R. Dendievel, Y. Bréchet, On the thermomechanical behavior of two-dimensional foam/metal joints with shear-deformable adherends – Parametric study. Composite Part B, doi:10.1016/ j.compositesb.2011.04.011.