W.M.B. Duval, D.J. Chato, and M.P. Doherty

Glenn Research Center, Cleveland, Ohio

Heat Entrapment Effects Within

Liquid Acquisition Devices

NASA/TM—2010-216789

September 2010

AIAA–2010–1298
NASA STI Program . . . in Profi le

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA Scientifi c and Technical Information (STI) 

program plays a key part in helping NASA maintain 

this important role.

The NASA STI Program operates under the auspices 

of the Agency Chief Information Offi cer. It collects, 

organizes, provides for archiving, and disseminates 

NASA’s STI. The NASA STI program provides access 

to the NASA Aeronautics and Space Database and 

its public interface, the NASA Technical Reports 

Server, thus providing one of the largest collections 

of aeronautical and space science STI in the world. 

Results are published in both non-NASA channels 

and by NASA in the NASA STI Report Series, which 

includes the following report types:

• TECHNICAL PUBLICATION. Reports of 

completed research or a major signifi cant phase 

of research that present the results of NASA 

programs and include extensive data or theoretical 

analysis. Includes compilations of signifi cant 

scientifi c and technical data and information 

deemed to be of continuing reference value. 

NASA counterpart of peer-reviewed formal 

professional papers but has less stringent 

limitations on manuscript length and extent of 

graphic presentations.

• TECHNICAL MEMORANDUM. Scientifi c 

and technical fi ndings that are preliminary or 

of specialized interest, e.g., quick release 

reports, working papers, and bibliographies that 

contain minimal annotation. Does not contain 

extensive analysis.

• CONTRACTOR REPORT. Scientifi c and 

technical fi ndings by NASA-sponsored 

contractors and grantees.

• CONFERENCE PUBLICATION. Collected 

papers from scientifi c and technical 

conferences, symposia, seminars, or other 

meetings sponsored or cosponsored by NASA.

• SPECIAL PUBLICATION. Scientifi c, 

technical, or historical information from 

NASA programs, projects, and missions, often 

concerned with subjects having substantial 

public interest.

• TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientifi c and 

technical material pertinent to NASA’s mission.

Specialized services also include creating custom 

thesauri, building customized databases, organizing 

and publishing research results.

For more information about the NASA STI 

program, see the following:

• Access the NASA STI program home page at 

http://www.sti.nasa.gov

• E-mail your question via the Internet to help@

sti.nasa.gov

• Fax your question to the NASA STI Help Desk 

at 443–757–5803

• Telephone the NASA STI Help Desk at

 443–757–5802

• Write to:

 NASA Center for AeroSpace Information (CASI)

 7115 Standard Drive

 Hanover, MD 21076–1320
W.M.B. Duval, D.J. Chato, and M.P. Doherty

Glenn Research Center, Cleveland, Ohio

Heat Entrapment Effects Within

Liquid Acquisition Devices

NASA/TM—2010-216789

September 2010

AIAA–2010–1298

National Aeronautics and

Space Administration

Glenn Research Center

Cleveland, Ohio 44135

Prepared for the

48th Aerospace Sciences Meeting

sponsored by the American Institute of Aeronautics and Astronautics

Orlando, Florida, January 4–7, 2010
Acknowledgments

This work is performed in support of the Cryogenic Fluid Management (CFM) program at NASA Glenn Research Center.

Available from

NASA Center for Aerospace Information

7115 Standard Drive

Hanover, MD 21076–1320

National Technical Information Service

5301 Shawnee Road

Alexandria, VA 22312

Available electronically at http://gltrs.grc.nasa.gov

Level of Review: This material has been technically reviewed by technical management. 
NASA/TM—2010-216789 1 

Heat Entrapment Effects Within Liquid Acquisition Devices

W.M.B. Duval, D.J. Chato, and M.P. Doherty

National Aeronautics and Space Administration

Glenn Research Center

Cleveland, Ohio 44135

Abstract

We introduce a model problem to address heat entrapment effects or the local accumulation of 

thermal energy within liquid acquisition devices. We show that the parametric space consists of six 

parameters, namely the Rayleigh and Prandtl numbers, the aspect ratio, and heat flux ratios for the 

bottom, side, and top boundaries of the enclosure. For the range of Ra considered 1 to 109

, beyond Ra on 

the order of 105

, convective instability is the dominant mode of convection in comparison to natural 

convection. The flow field transitions to asymmetric modes at Ra on the order of 107

. Direct numerical 

simulation of a large geometric length scale prototype for Ra on the order of 109 shows that the flow field 

evolves from small wavelength instability which gives rise to nonlinear growth of thermals, propagation

of the instability occurs via growth of secondary and tertiary modes, and a travelling wave mode occurs 

prior to asymmetry. The effect of a large aspect ratio is to increase the number of modes in the vertical 

direction. Due to the slow diffusion of heat in the prototype, asymptotic states are not readily attained, we 

show that dynamical similarity can be used for a model which allows the attainment of asymptotic states 

and that transition to a chaotic state occurs for Ra on the order of 109 via a broadband power spectrum. 

These dynamical events show that for the baseline condition in which heat is absorbed from background 

laboratory environment, higher heat flux is absorbed at the top and bottom boundaries of the enclosure 

than a nominal value of 34.9 erg/cm2

-s.

Nomenclature

Ra Rayleigh number

Ar aspect ratio

Pr Prandtl number

fs sidewall heat flux ratio

fb bottom heat flux ratio

ft top heat flux ratio

q” imposed heat flux

H cavity height

L cavity width

g gravitational acceleration

k thermal conductivity

α thermal diffusivity

β thermal expansion coefficient 

ρ density

Γ boundary of cavity
NASA/TM—2010-216789 2 

Introduction

Capillary liquid acquisition devices are used inside cryogenic storage tanks in order to supply vaporfree propellant to the engine. The liquid acquisition devise (LAD) interfaces directly with the feed system, 

for the supply of propellant, which can serve as a significant source of heat leak. In order to avoid 

excessive boil-off and ensure sufficient propellant subcooling to prevent cavitation and vapor formation 

during propellant outflow to the engine, the control of heat leak is necessary. Hence the localized 

accumulation of thermal energy, defined as heat entrapment, within the LAD flow channels becomes an 

issue of prime importance. We consider a model problem that addresses the effect of heat source at the 

boundary of an enclosure which contains typical cryogenic liquids (nitrogen, hydrogen) to lead insight 

into heat entrapment effects. This model problem has relevance towards the experimental works of 

Bolshinskiy et al. (Ref. 1), who considers heat entrapment effects modeled as controlled bottom heating 

of a cylindrical enclosure with and without partitions totally filled with a fluid. The partitions, which 

consist of typical meshes used in LADs as well as a solid divider, serve to emulate typical experimental 

conditions encountered during operations. Their results (Ref. 1) show that convection occurs over a time 

scale much greater than the time for transport through the screen meshes, thus the convection problem can 

be decoupled from that of transport. This is verified experimentally (Ref. 1) from comparison between 

data for meshes versus solid boundary which show that the temperature rise are not significantly different. 

Hence a solid boundary can be used as an approximation of the boundary condition to address heat 

entrapment effects. In our model we employ a solid boundary condition to address convection effects 

encountered in LAD applications. 

The basic physics of heat entrapment effects or the localized accumulation of thermal energy is 

addressed by broadly considering the effect of heat flux absorption on the boundary of an enclosure, due 

to controlled heat flux on the bottom and unknown heat leakage at the top and side boundaries, on an 

enclosed cryogenic liquid initially at saturation condition which totally fills the volume of the enclosure. 

Of basic interest is the rise in temperature due to various external heat flux conditions. The quantification 

of the temperature field can be used to assess optimum operating conditions for propellant flow to the 

engine. Related to the broad problem of the effect of heat absorption on the temperature rise of a 

cryogenic fluid are also issues on effect of thermal stratification on temperature rise due to side heating of 

partially filled cavities with a free surface as addressed by Tellep and Harper (Ref. 2), Vliet et al. (Ref. 3), 

and Vliet (Ref. 4) with applications to cryogenic storage tanks. Approximate boundary layer analysis 

using an integral approach coupled with experiments (Refs. 2 and 3) show that the effect of dominant side 

heating is to transport warmer fluid along the wall to the interface region which becomes accumulated at 

the free interface thus forming a stratified layer. Vliet (Ref. 4) addressed the effect of bottom heating 

using an approximate integral approach and show agreement with experimental data; the results show that 

as bottom heating increases there is a faster approach toward asymptotic conditions. Evans et al. (Ref. 5)

have addressed the effect of non uniform heating on thermal stratification of water inside a rectangular 

enclosure and shown that thermal stratification can be reduced by distributing the heat flux along the side 

wall such that the maximum heat flux is concentrated near the bottom and decreases toward the top. The 

primary mechanism for the reduction of thermal stratification was postulated to be boundary separation at 

the wall. The bottom of the enclosure (Ref. 5) was insulated while the free interface was coated with a 

surfactant film of stearic acid to prevent vaporization; though the condition at the free surface is not 

specified, the effect of the surfactant film may be approximated as an insulating boundary condition based 

on the asymptotic profile of the rise in temperature of the enclosed fluid.

The experimental works of Evans & Stefany (Ref. 6) addressed transient heat transfer effects inside 

totally filled cylindrical enclosures. They show that in response to a step-change in wall temperature, the 

mixed-mean temperature of the fluid exhibits an exponential rise in temperature which approaches the 

wall temperature asymptotically. In addition, the heat transfer coefficient oscillates in time as steady-state 

is approached. This finding was confirmed by the numerical and experimental works of Barakat & Clark

(Ref. 7) for a partially filled cylindrical enclosure. Numerical works have been put forth by Polezhaev
NASA/TM—2010-216789 3 

(Refs. 8 and 9) which consider natural convection inside partially filled cavities for rectangular (Ref. 8) as 

well as cylindrical (Ref. 9) cavities for various boundary conditions.

In contrast to the above works, our problem is somewhat unique in that there is no free surface and 

also that the fluid occupy the entire volume of the enclosure. The experimental works of Evans & Stefany

(Ref. 6) share the similarity of total volume filling, however the boundary conditions are different, a stepchange in temperature is prescribed. In our model we consider the Neumann prescribed heat flux 

conditions along the bottom, top and side of the enclosure. Thus depending on the dominant heat flux 

condition, there is coupling between natural convection driven by side heating and convective instability 

driven by bottom heating. The top heating serves as a stabilizing mechanism, due to imposed stable 

stratification, in contrast to bottom heating in which there is unstable stratification due to an increase of 

density with height. The dominant mechanism, natural convection or convective instability, depends on 

the external heat flux conditions during LAD operation. These two mechanism differ in that in natural 

convection problems there is no stability threshold, convection occurs for any temperature variation 

across the boundary of the cavity whereas convective instability has a stable solution. The characteristics 

of natural convection inside vertical cavity due to differentially imposed temperature, at the vertical 

boundary, have been addressed by Elder (Refs. 10 and 11) for laminar (Ref. 10) and turbulent (Ref. 11) 

free convection. Transition from unicellular to secondary and tertiary flows have been shown to occur 

during laminar free convection, whereas the turbulent regime is characterized via travelling wave-like 

motion along the wall. The visual observation by Seki et al. (Ref. 12) have confirmed the findings of 

Elder, and also demonstrated the effect of Prandtl number. The mechanism of convective instability has 

been demonstrated for differentially heated horizontal cavity, with a forced flow component, whose 

bottom temperature is greater than the top by Akiyama et al. (Ref. 13) They show that for forced 

convection between horizontal plates, the limit of zero Reynolds number (no forced flow) yields 

agreement with the critical Rayleigh number criterion of 1708 corresponding to the classical RayleighBenard problem (Ref. 14). Secondary flow patterns in the form of longitudinal vortex rolls were observed

to occur near the critical Rayleigh number. The two basic configurations of a narrow vertical and 

horizontal cavity serve as basis by which the interaction between natural convection and convective 

instability mechanism can be understood.

We introduce a model problem which considers the basic physics of heat entrapment inside an 

enclosure subjected to prescribed heat flux conditions to shed light on the convective state of the flow 

field in relation to the temperature field. This model problem is intended to lead insight into the 

temperature measurements (Ref. 1) made for liquid nitrogen inside an enclosure which absorbs heat from 

ambient condition as well as the effect of imposed heat flux condition at the bottom boundary. Our main 

focus is to investigate the effect of heat flow from the ambient laboratory environment to the enclosure in 

order to establish a baseline to address the effect of imposed heat flux on the bottom boundary. The 

experiments were conducted using double encapsulation along the sidewall, thus the heat flux absorbed 

along the side wall is minimal; we assume a value of heat flow on the order of 1 μW. However, other 

unknowns are the heat absorbed on the top and bottom boundary for the experiments carried-out under 

ambient laboratory conditions. We investigate a range of heat flux condition to lead insight into its effects 

on the temperature rise of the fluid. This approach is necessary, since the baseline condition is not known 

and it is needed in order to assess properly the effect of bottom heating. 

We show that the effect of heat absorption is to destabilize the equilibrium condition of a fluid inside 

an enclosure and give rise to flow. The dynamical motion of the flow field obeys the evolutionary 

equations of mass, momentum and energy. For an incompressible Boussinesq fluid we numerically solve 

the coupled set of governing equations consisting of the Navier-Stokes coupled with the scalar 

temperature field using finite-difference methods. We consider the range of Rayleigh numbers (Ra) of 

order 1 to 109 for heating conditions comprising of prescribed heat flux at the sidewall, top, and bottom of 

the cavity. For the sidewall heat source we show that there is no stability threshold, incipient convection 

occurs. The flow field is dominated by two counter-rotating cells for low Rayleigh numbers and 

secondary flows can occur as Ra increases. However, for the bottom heating condition there is a stability 

threshold for which the onset of convection occurs which is on the order of 105

. As the Rayleigh number 
NASA/TM—2010-216789 4 

increases, the horizontal wavelength which characterizes the onset of convection decreases. There is a 

transition to asymmetric flow for Ra on the order of 107

. Direct numerical simulation of the full scale 

prototype shows that for Ra on the order of 109 the onset of convective instability proceeds from short 

wavelength to nonlinear growth of thermals and culminates with a travelling wave mode prior to 

transition to asymmetry. Due to the large geometric length scale of the prototype, asymptotic states can be 

prohibitive. We use dynamical similarity using a model with smaller geometric dimensions to show that 

similar dynamical events occur which facilitates the attainment of asymptotic states. As the Rayleigh 

number increases from 106 to 109 the flow field undergoes a transition from periodic to chaotic 

oscillations. The parametric study helps to identify the basic mechanism that drives the flow and lays the 

foundation to identify the baseline condition for heat entrapment effects.

The paper is organized as follows: we postulate the basic physics of heat entrapment inside of an 

enclosure. We derive the parametric space of the problem through scaling analysis. We solve the coupled 

set of nonlinear partial differential equations using finite difference methods coupled with a FluxCorrected Transport method to resolve the small scale structure of the temperature field. We show that

direct numerical simulation coupled with dynamical similarity can be used advantageously to model the 

large prototype enclosure used in the experiments. We show the basic convective state of the flow field 

for a relevant parametric space and that transition to asymmetric flow is relevant to this problem. 

Formulation

The baseline condition for heat entrapment is shown fundamentally in Figure 1A(a) and consists of a 

cryogenic liquid at saturation condition represented by ρA which can absorb various amounts of heat at its 

boundaries comprising of surface heat fluxes qs

” , qt” , qb” denoting respectively the wall, top, and bottom 

boundaries. The enclosure in Figure 1A can be approximated geometrically as the mid-plane of a cylinder 

with equivalent height (H) and diameter (L). This basic configuration represents the baseline condition 

investigated experimentally (Ref. 1) for liquid nitrogen inside a cylindrical Dewar under ambient 

laboratory condition; for which the heat flux absorb at the sidewall is controlled by external means such 

as double encapsulation of the wall region to minimize heat flow. Regardless of precautionary measures a 

finite amount of heat leakage occurred as indicated by the rise of temperature (Ref. 1) of liquid ρA. This 

finding indicates that under laboratory environment, equilibrium condition is not attained for a given 

cryogenic liquid ρA. One of the main focus of our model is to investigate the effect of relative magnitude 

of the heat fluxes (qs

”

, qt

”

, qb

”

) on the dynamical state of the liquid ρA at saturation temperature TA

 when 

the system is driven off-equilibrium. This approach is necessary as it eliminates one of the unknown of 

the problem, if one of the heat flux is prescribed such as qb

” in the experiments (Ref. 1). Thus for a given 

qb

”

, optimal values for qs

” and qt” can be estimated based on ambient condition and experimental 

arrangement. Given the double encapsulation used on the sidewall of the experiments (Ref. 1), it is 

approximated that qs

” << qt”

. Hence for the baseline condition the dominant heat fluxes are qb

” and qt”

. 

More generally the effect of heat entrapment is investigated experimentally (Ref. 1) by placing a 

divider at mid-height, consisting of stainless-steel screen meshes or an aluminum solid divider, of the 

cavity as shown in Figure 1A(b) in which ρB = ρA. The effect of heat entrapment is quantified 

experimentally by measuring the temperature above and below the divider which is used to compare 

against the baseline experiment illustrated in Figure 1A(a). The findings (Ref. 1) indicate that 

temperatures below (bottom)/above (top) the divider position were higher/lower with the divider installed 

indicating heat entrapment/reduced stratification. The effects of the divider may be quantified by 

approximating the local heat flux qi

” into or out of the top region. Hence the problem may be simplified 

into two sub-problems as in Figure 1B in which heat flows into or out of the top/bottom region shown in 

(a) and (b). This simplification effectively reduces the aspect ratio by 1/2 and the problem may be 

similarly addressed as the baseline problem in Figure 1A(a) in which qb

” and qt” can take the respective 

values of qi

”

. This approach effectively reduces the parametric range of the problem which can be used to 

investigate the more general set of conditions as represented in Figure 1B.
NASA/TM—2010-216789 5 

In general, as shown in Figure 1A(a), a condition for equilibrium to exist is for qs

” = qt” = qb” = 0; it is 

not likely for a cryogenic fluid inside an enclosure to meet this condition. Thus off-equlibrium condition 

need to satisfy the dynamical motion of the fluid given by the continuity, momentum and energy 

equations as:

0 u v

x y

 

    (1)

2 2

2 2

u u u p uv u v

t xy x x y

        

      

        

(2)

2 2

2 2 0

v v v p vv u v T ng txy y x y

                   

        

 (3)
NASA/TM—2010-216789 6 

2 2

2 2

T T T TT

u v

t xy x y

       

     

       

    (4) 

In the above equations the density field is independent of pressure due to incompressibility of the liquid, 

however it is a function of temperature and can be linearized and given as an equation of state,

     ( ) T

A [1 ] T (5)

We consider a Boussinesq fluid in which variation of density is significant only in the body force term as 

shown in Equation (3) where the coefficient of thermal expansivity is given as 1 .  P      T

Superscript ~ denotes dimensional temperature and  T TTA   . The thermophysical properties of the 

liquid, its density ρ, kinematic viscosity ν, and thermal diffusivity α are taken to be at saturation condition 

for the specific cryogen, for liquid nitrogen TA

 = 77 K. The gravitational acceleration gy is denoted by 

ng0 in which n is a ratio by which the standard gravitational acceleration go can be reduced via low 

gravitational environments such as the International Space Station (ISS), Moon, and Mars surface.

The description of the initial and boundary conditions are made in reference to Figure 1A(a). The 

initial condition is taken as uniform temperature throughout the liquid,

t Txy T   0, ( , , 0) A   (6) 

The momentum Equations (2) and (3) satisfy the no-slip condition along the boundary of the cavity 

denoted by Γ and given as,

u x y t v x y t on ( , , ) 0, ( , , ) 0   (7) 

Whereas the energy Equation (4) satisfy the Neumann condition along the boundary Γ prescribed as,

" " " 0, ,0 , ) / s x s 0 x q q yH T q k x         

 (8)

" " " , ,0 , ) / s xL s x Lq q y H T q k x         

 (9)

" " " 0, ,0 , ) / b y b 0 y q q xL T q k y         

 (10)

" " " , ,0 , ) / t yH t y Hq q x L T q k y         

 (11)

Scaling Analysis and Computational Method

The governing set of Equations (1) to (4) may be simplified by transforming the momentum 

Equations (2) and (3) into a vorticity equation using the following definition,

, , v uu v

y x xy

             (12)
NASA/TM—2010-216789 7 

The resulting set of equations is scaled using the characteristics length, time, and velocity as 

2 (, / , / ) LL L   , respectively. The dimensionless temperature is scaled as

" /

A T T

T

qL k

 

 

(13)

The scaled dimensionless (denoted by *) set of equations and boundary conditions become, 

2* 2* *

*2 *2

x y

      

(14)

* * * 2* 2* * *

* * * *2 *2 * Pr Pr T

u v Ra

t x y xy x

              

   

       

(15)

2 2 * *

* * * *2 *2

T T T TT

u v

t x y xy

          

  

      

(16)

The boundary conditions become, 

* * *

* 0 0, 0 , ) s x x y Ar T f

x         (17)

* * *

* 1 1, 0 , ) s x x y Ar T f

x         (18)

* * *

* 0 0, 0 1, ) b y y x T f

y         (19)

* * *

* , 0 1, ) t y Ar y Ar x T f

y         (20)

The set of parameters resulting from scaling is defined as,

" 4 " " "

" "" , Pr , , , , o sbt

sbt

q ng L H q q q Ra Ar f f f k L qqq 

       (21)

which represents respectively the Rayleigh and Prandtl numbers, aspect ratio, and dimensionless heat 

fluxes. Thus the parametric set, 

( , , Pr, , , )  Ra Ar f f f sb t (22)

consists of six parameters that can be varied independently and represents a formidable parametric set. 

Simplification of the parametric range is obtained by addressing the relevant limiting conditions for a 

given aspect ratio, fixed cryogenic liquid, and typical heat flux scale q”

. Since the predominant control 

variable is the bottom heat flux qb

” which implies fb = 1, then fs, ft can take a range of values from zero to 

one. The conditions ft = 0 and ft = 1 signify insulated and heat absorption top boundary conditions. The 
NASA/TM—2010-216789 8 

side condition fs = 0 and fs = 1 indicate a perfectly insulated sidewall and dominant sidewall heating. 

Since the sidewall has double encapsulation in the experiments (Ref. 1), we consider values of fs on the 

order of 10–6 which is approximately 1 μW of heat leakage. Although this heat leakage is insignificant in 

comparison to fb and ft, it serves as a perturbation parameter to excite flow instability.

The governing set of Equations (14) to (16) with boundary conditions (Eqs. (7) and (17) to (20)) are

solved using finite-difference methods with a direct solver approach. The solution methodology consists 

of solving the stream function Poisson’s Equation (14) using matrix inversion which is followed by the 

solution of the energy Equation (16) and vorticity Equation (15). The vorticity field is used subsequently 

for each time step. Since we are interested in time evolution of the flow field dynamics, a time accurate 

methodology is employed which uses third order Adams Bashforth for time discretization. Since the 

bottom heat flux is dominant, the convective instability mechanism becomes important in comparison to 

natural convection, which requires high resolution of the temperature field. We employed the Flux 

Corrected Transport (FCT) method (Ref. 15) to resolve sharp interior gradients; in addition this method 

eliminates false numerical oscillation by ensuring positivity of the temperature scalar component 

throughout the domain of computation.

Numerical Results

Stabilization of Flow Field

The dynamical state of fluid motion inside the enclosure shown in Figure 1A(a) is first investigated 

for ideal cases in order to show the fundamental mechanism of convective motion. We first consider how 

stable states can be obtained and the cause of departures from stability. This is illustrated using liquid 

hydrogen as an example with Prandtl number Pr = 1.03. The destabilization of the flow field for bottom 

and side heat flux is summarized in Figure 2. The heating prescriptions that leads to convective stability 

in which there is no fluid motion even though conduction dominates is illustrated for Ra = 1.26×103

. 

These trends show that stability of the flow field can be obtained provided that there is a perfectly 

insulated wall fs = 0 for either bottom heating only (fb = 1, ft = 0), top heating only (fb = 0, ft = 1) or both 

top and bottom absorbing heat (fb = 1, ft = 1). The case when (fb = 1, ft = 1) represents actual experimental 

condition, since neither top ft = 1 nor bottom heating fb = 1 can be isolated. However these cases represent 

idealization since there is finite heat leakage along the wall, i.e., fs>0, regardless of precautionary 

measures. 

Perturbation of Stability

The perturbation of the instability for bottom (fb = 1, ft = 0) and top heating (fb = 0, ft = 1) is illustrated 

for Ra = 1.26×106 in order to isolate simultaneous effects. For a conservative heat leakage on the order of 

1 μW projected on the sidewall surface of a cylinder fs = 2×10–8, the stable base state becomes unstable 

due to unstable stratification caused by heavy fluid overlaying lighter fluid due to decrease of local 

density near the bottom of the enclosure caused by bottom heating (ft = 1). This case corresponds to qb

” = 

500 erg/cm2

-s which can be obtained from heat flow from ambient laboratory environment on the order of 

10 mW. In this case convective instability driven by bottom heating is dominant over natural convection 

driven by side heating. Since we are operating near the threshold of instability the growth of four modes 

occur. In contrast to instability of the bottom stratification (fb = 1, ft = 0, fs = 0) due to side heat leakage of 

fs = 2×10–8, we show that the destabilization of the top stratified layer (fb = 0, ft = 1, fs = 0) due to the same 

heat leakage gives rise to a symmetric double cell. This case illustrates that even though the top stratified 

layer is stable for all Rayleigh numbers, the effect of finite side heat leakage is to drive convection.
NASA/TM—2010-216789 9 
NASA/TM—2010-216789 10

Effect of Side Heating

Since side heat leakage is an effective source to drive instability, we isolate the effect of side heating 

by considering the case when there is no heat leakage from the bottom (fb = 0) and top (ft = 0) contrasting 

Ra = 2.5 and Ra = 1.26×106

. These two cases contrast heat flow on the order of 1 μW (fs = 2×10–8) and

10 mW (fs = 1) corresponding to Ra = 2.5 and Ra = 1.26×106

, respectively. For Ra = 2.5 natural 

convection becomes dominant, it generates a very weak flow field which does not disturb the conduction 

layer near the vertical walls. However, when there is significant side heat leakage, Ra = 1.26×106

, the 

cellular structure becomes unstable which results in growth of impending secondary modes due to shear 

instability. There is a transition from conduction dominated regime to a boundary layer regime indicated 

by the temperature field. The flow field transports warm fluid near the walls to the interface region which 

forms a stratified layer. The scenario of sidewall dominated heating is similar to the events that occur for 

applications to address the effect of thermal stratification in partially filled cryogenic storage tanks

(Refs. 7 and 8). The case for Ra = 1.26×106

 (fb = 0, ft = 0) is demonstrated for illustrative purposes; the 

precautionary measures taken to insulate the sidewall in the experiments (Ref. 1) would indicate a much 

smaller side leakage toward the lower limit of Ra = 2.5.

Perturbation of Double Layer

A typical heating condition is illustrated for the case when there is simultaneous heat absorption from 

the top and bottom boundaries which represents perturbation of the double layer for Ra = 1.26×106

. The 

destabilization of the conduction dominated regime (fb = 1, ft = 1, fs = 0) is illustrated for fs = 2×10–6 and

fs = 1. When fs = 2×10–6 convective instability dominates, there occurs growth of four modes at the bottom 

region while there is growth of a stratified layer at the top region for early time t = 360 s. However, if 

there is significant sidewall heat leakage (fb = 1, ft = 1, fs = 1) natural convection dominates and the 

characteristic of convection switches from an instability mechanism to a boundary layer dominated flow 

shown by the convective cells along the wall and secondary cells near the bottom. In this case the 

intensity of convection modifies the top stratified layer and there is a potential for local superheating near 

the top corners. In summary, the mechanism by which the double layer (fb = 1, ft = 1, fs = 0) becomes 

unstable is illustrated by perturbation of the stability with fs = 2×10–6, and the effect of dominant side 

heating on the stability is illustrated for fs = 1. This indicates that there can be a wide range of dynamical 

events that occurs for 2×10–6<fs<1, even though we fix fs = 2×10–6 for the parametric study, it is a key 

parameter to characterize heat entrapment effects. The results also show that when the double layer is 

perturbed convective mixing becomes important. These limiting conditions serve to illustrate the 

interaction of convective instability and natural convection that can occur for the baseline condition 

shown in Figure 1A(a).

Mechanism of Instability 

For the given experimental conditions, it is likely that the convective characteristics inside the 

enclosure are dominated by convective instability (fs = 2×10–6). In Figure 3A we illustrate the mechanism 

of instability as the Rayleigh number increases. The results show that from the threshold of instability

Ra = 1.26×105 there is a transition from long to short wavelength as Ra increases to 1.26×109

; the number 

of modes does not change beyond Ra = 1.26×107

. The growth of modes segregates near the walls as Ra 

increases and the wavelength between the cells decreases. This trend indicates that high Rayleigh 

numbers require resolution of small scales for computational purposes.
NASA/TM—2010-216789 11
NASA/TM—2010-216789 12

The characteristic wavelength ( */ )  L as a function of the Rayleigh number in Figure 3B

shows that  is inversely proportional to Ra, thus decreases sharply as Ra increases. For the highest 

Rayleigh number considered (1.26×109

),   0.075 which is on the order of 4 mm. This finding indicates 

that the length scale which needs to be resolved gets smaller and smaller for Ra > 109 and also illustrates 

the challenge posed for computations of high Rayleigh numbers.

Dynamical Events

We now turn to the dynamical events as the flow field approaches asymptotic states. This is 

illustrated for Ra = 1.26×106 in Figure 4, which shows growth of a single thermal, and Ra = 1.26×107 in 

Figure 5A which shows growth of double thermals and the breakdown of symmetry in Figure 5B. In 

contrast to natural convection dominated flows in which the maximum temperature occurs near the

sidewalls due to side heating, in this convective instability dominated situation, the maximum temperature 

occurs near the mid-center of the cavity due to bottom heating since the sidewalls are nearly insulated 

(fs = 2×10–6). The mechanism of the instability is due to the release of potential energy from unstable 

stratification caused by heavy fluid overlaying lighter fluid and its conversion to kinetic energy to drive 

fluid motion. In this case, the warmer fluid transported upwards and cooler fluid moving downwards 

result in the release of potential energy; the inner cells grow faster than the outer cells since the inner 

region has the maximum temperature. The inner cells grow as time increases, thus the instability self

propagates until the top wall is reached. The initial bottom four cells (t = 336 s) self-organizes into two 

rows of secondary cells (t = 408 s) due to the nonlinear growth of the thermal. Since the top of the cavity 

is insulated (ft = 0), the intensity of convection increases as the top of the cavity is approached. Buoyancy 

effects due to growth of the thermal propel the two inner cells toward the top of the cavity (t = 360 s) 

while the outer cells near the walls segregate toward the center of the cavity (t = 384 s). The maximum 

temperature decreases with time due to local mixing induced by the increase in convective intensity for 

Ra = 1.26×106

. 
NASA/TM—2010-216789 13
NASA/TM—2010-216789 14
NASA/TM—2010-216789 15

We consider the dynamical evolution of the flow field when there is six modes in which there is heat 

absorption at the top boundary as shown in Figures 5A and B for Ra = 1.26×107

. There is similarity of 

events in comparison to the top insulated case for Ra = 1.26×106

. Since the maximum temperature occurs 

toward the mid-section of the cavity, the inner cells grow faster than the outer cells, the growth of six 

modes leads to formation of two symmetric thermals with respect to the center-line. The instability 

produces buoyancy which lifts the two central cells thus forming tertiary modes (t = 130 s); the rows of 

cells from the center toward the wall form co-rotating cells with the middle cell (tertiary mode) counterrotating. The tertiary modes propagate toward the top of the cavity (t = 135 s) from a bifurcation of a local 

double homoclinic orbit which signals transient effects. There is subsequent merging of the smaller cells 

(tertiary modes) along the bottom, forming a single homoclinic orbit, into single cells (t = 140 s). The 

growth of the stratified layer from the top limits the penetration of the modes as shown in Figure 5B (t = 

960 s); thus the system transitions to four modes. These four modes become asymmetric (t = 1080 s) and 

the growth of the asymmetry leads to merging of the inner and outer cells (secondary modes). Thus two 

asymmetric modes are formed at t = 1200 s. This finding is significant as it implies that beyond the 

critical Rayleigh number the flow field becomes asymmetric. This implies that axi-symmetric models 

would be limited in their application to high Rayleigh numbers since this occurs for Ra = 1.26×107

. 
NASA/TM—2010-216789 16

Local Measures

The local measures which characterize the behavior of the fluid as a function of system parameters 

are summarized in Figures 6A, B, and C. The effect of increasing the Rayleigh number is considered in 

Figures 6A and B by showing the characteristic maximum velocity magnitude and temperature as a 

function of Ra contrasting the effect of insulated (ft = 0) and top absorbing (ft = 1) boundaries at t = 600 s. 

The results in Figure 6A show that Vmax remains constant for both ft = 1 and ft = 0 for Ra < 105 indicating 

stability or no convection. Instability occurs for Ra > 105 for which the solution remains the same up to 

107

. There is faster rise in velocity for the case ft = 0 for 107 < Ra < 109. This represents the range for
NASA/TM—2010-216789 17

which the stabilizing effect of the top boundary ft = 1 becomes effective and lowers the convective 

velocity. The trend of the maximum temperature as a function of Ra in Figure 6B shows that Tmax is 

higher for ft = 0 when Ra < 107 whereas for Ra > 107

 Tmax becomes greater for the stabilizing top heat 

absorption ft = 1 case. These trends indicate that local convective mixing is much more effective at 

lowering Tmax for the insulated top as opposed to the case when there is heat absorption on top. In terms 

of heat entrapment (localized increase in temperature), the trend from ft = 1 to ft = 0 is representative of 

expectations, that is an increase in local maximum temperature for Ra < 107

; since the effect of the 

divider is to limit heat flow, the region bounded by 0<ft<1 represents potential conditions where Tmax

increases beyond the limiting condition ft = 1. The effect of the divider approximates the trend for ft<1, 

thus ft = 0 represents a limiting condition for heat entrapment. For Ra > 107 the increase in convective 

effects for ft = 0 induces local mixing which decreases heat entrapment effects. In Figure 6C, we show the 

transient response of a cryogenic fluid along the central axis (0.5, y) of the cavity for ft = 1 and ft = 0. The 

trend of an adverse temperature gradient at the bottom of the enclosure is not greatly affected by the top 

conditions. The effect of top heating (ft = 1) is to create a stably stratified layer which increases with time. 

Whereas, when the top is insulated ft = 0, there is no stratified layer; the temperature becomes nearly 

uniform near the top as the asymptotic state is approached. 

Direct Numerical Simulation/Dynamical Similarity

We demonstrated the characteristic of the flow field and the response of the temperature field for a 

fixed Pr, Ar, and a range of Rayleigh number for a model. We now consider the prototype used in the 

experiments1 which has much greater geometric length scales, that is L = 19.1 cm, H = 61.9 cm. Due to 

the large length scale employed in the prototype it requires a much larger grid density to resolve the small 

length scales in the flow and temperature fields that occurs for large Rayleigh numbers Ra ≥ 108 as shown 

in Figure 1A(a). For example the horizontal length scale of the prototype can be four times greater than 

that of the model; for a typical grid size used for the model of 200×200, the prototype would require a 

grid density four times as large of 800×800 to resolve the same length scale as the model. This can require 

more than doubling the computational resources; a typical case for a model which requires between 30 to 

72 hr using 4 Central Processing Units (CPU’s) now requires between 300 to 800 hr using 20 to 32 CPU’s

for the prototype. Even though, direct numerical simulation of the prototype is the most exact, dynamical 
NASA/TM—2010-216789 18

similarity can be used to model the prototype. This leads to considerable savings in computational 

resources. We explore the application of dynamical similarity by showing that a model with smaller 

geometric length scale, requires matching of the parametric set ΛP (given by Eq. (22)) for the prototype to 

that of the model ΛM; we also show that dynamical events are similar with the exception of the time scale. 

The issue being that due to reduction of length scale of the model, the equivalent heat flux diffusing 

through the model takes much shorter time. The model allows the dynamical events to be accelerated, 

hence the savings in computational resources.

In order to address direct numerical simulation of the prototype (for a geometric length scale

L = 19.1 cm, H = 61.9 cm), we consider Ar = 1 as a base case, Ar = 1.32 (H/2) for the half scale model in

Figure 1B, and Ar = 3.26 (H) for the full length of the prototype in Figure 1A. The external condition for 

the prototype is absorption of heat on the order of 1 mW for the top and bottom boundaries, due to 

ambient background laboratory condition, which is equivalent to its projection on the surface of a circle 

with heat flux of q” = 34.9 erg/cm2

-s for liquid nitrogen. The condition at the sidewall is heat absorption 

on the order of 1μW or a heat flux q” = 0.001 erg/cm2

-s. The characteristics of the dynamical evolution of 

the flow field is illustrated for Ar = 1 in Figures 7A, B, and C which show three main events: (a) 

transition to short wavelength instability and growth of thermals, (b) vertical propagation of the 

instability, and (c) travelling wave mode which transitions to asymmetry.

Even though the cavity size for the prototype in Figure 7A is much larger than a typical model shown 

in Figure 3A (Ra = 1.26×109

), for nearly equivalent Rayleigh numbers and slightly different Prandtl 

number, similar transition instability occurs in Figure 3A at t = 12 s as compared to Figure 7A at t = 288 s

for Ra = 1.0×109 and Pr = 2.27. In this case there is a deceleration of the instability since the size of the 

cavity is so large for a given heat flux input. There is a similar growth of thermals to that shown in

Figure 5A, however the higher grid density used (500×500) allow the resolution of small length scales 

such as the fine structures of the interaction between thermal pairs shown in Figure 7A at t = 432 s. The 

propagation of the instability toward the top of the cavity occurs via a repetitive sequence of growth of 

secondary cells in Figure 7B. The mechanism of the instability occurs through the formation of 

homoclinic orbits (t = 648 s and t = 720 s) which serve to increase the size of the cells and form 

secondary cells inside the main cell which detaches and form secondary cells as shown from the transition 

between t = 720 s to t = 792 s. Once the flow field penetrates near the top of the cavity, there is a 

travelling wave mode event that occurs as shown in Figure 7C. The travelling wave mode occurs through 

the exchange of energy between tertiary and secondary modes as shown for t = 864 s, t = 1296 s, and

t = 1368 s. As the wave travels up and down the cavity there is a breakdown and reconnection of the 

vertical modes; there is also boundary layer entrainment along the vertical walls. The flow field 

subsequently transitions to asymmetry. Note in comparison to Figures 5A and B the boundary layer 

extends much further up the cavity vertical walls.

Effect of Aspect Ratio

The effect of increasing the height of the cavity to actual prototype model Ar = 1.62 (H/2) and Ar = 

3.24 (H) for long and short time is shown in Figure 8 contrasting insulated ft = 0 and absorbing top ft = 1 

boundary conditions. For the case of a taller cavity Ar = 1.62, the growth of tertiary modes is formed 

between two co-rotating secondary cells. The tertiary mode has a homoclinic orbit which signals 

propagation of travelling wave inside a fluid. This case is contrasted with Ar = 3.72, the prototype used in

the experiments (Ref. 1). Such a large cavity requires high grid density to resolve small scales; we show 

direct numerical simulation of early time evolution of the instability using 1000×1000 grids for short 

times. The results show the onset of instability which consists of six modes segregated toward the wall 

with small length scales. In this case, the characteristic time for the instability to occur increases

(t = 288 s) due to the large increase in vertical length scale as compared to a model in Figure 3A which 

shows the instability at t = 12 s; note the instability develops faster with a lower aspect ratio prototype as 

in Figure 7A at t = 288 s.
NASA/TM—2010-216789 19
NASA/TM—2010-216789 20
NASA/TM—2010-216789 21
NASA/TM—2010-216789 22

We illustrate the application of dynamical similarity using a model in Figure 9 with smaller length 

scale for the same heat flux boundary condition as the prototype in Figure 7A. In this case,

( , , Pr, , , ) ( , , Pr, , , ) M bP b s s t t  Ra Ar f f f Ra Ar f f f   (23)

note the matching of the parametric space yields similar dynamical events as compared to Figure 7A, 

except with the model the characteristic time for the instability to develop is much shorter t = 20 s as 

compared to t = 288 s for the prototype in Figure 7A. 

Local Flow Dynamics at Fixed Points

The local dynamics of the flow and temperature fields at fixed points (xp, yp) as denoted in

Figure 1A(b) is shown in Figure 10. The decomposition of the time history signal into it frequency

components for the flow field is obtained by estimating its power spectrum using,

2

2

0

1 ( ) ()

Tb

if t

V

b

P f Vte T

− π = ∫  (24)

in which V(t) is the time history of the corresponding component of the velocity field at the time interval 

of the data denoted by Tb. The smooth estimate PV(t) is obtained from a convolution relationship 

( ( ) ( ) ( )) Pf W f Pf V HV = ∗  using a Hanning spectral window (WH(f)). The dynamical characteristics of 
NASA/TM—2010-216789 23
NASA/TM—2010-216789 24
NASA/TM—2010-216789 25

flow and temperature fields at various fixed points (xp, yp) for the prototype are shown in Figure 10. In 

relation to Figure 1A(b), the series of points represent location near the left wall which proceeds from the 

bottom toward the top for a time interval Tb = 1500 s. As shown in Figure 10, the time history of the 

temperature and flow fields in the bottom region at (0.002, 0.002) shows that for a typical q” =

34.9 erg/cm2

-s in which (fb = 1, ft = 1, fs = 2.9×10–6), the temperature increases initially by conduction 

effects up to t = 300 s while the flow field remains quiescent; for t > 300 s the onset of instability occurs 

in which the flow and temperature fields oscillate periodically. As the distance increase from (0.15, 0.1) 

to (0.15, 0.35) at the bottom, it takes longer (t = 400 s) for the temperature and flow fields to respond and 

there is also an increase in their aperiodicity. As the distance is further increased toward the top of the 

cavity, from (0.15, 0.65) to (0.002, 0.998), the flow and temperature oscillations transition from aperiodic 

to a nearly periodic (mono-frequency) oscillation. This trend is due to absorption of heat (ft = 1) at the top 

boundary which stabilizes the flow. The power spectrum shows that there is a dominant mono-frequency 

component of 0.008 Hz with slight harmonics; this corresponds to a long period of 125 s. Note that near 

the top (0.002, 0.998) conduction effects dominates and it takes on the order of t = 800 s for the onset of 

convective instability. The behavior of the wall temperature at the mid-section of the bottom (0.5, 0.0) 

and top (0.5, 1.0) shows that at the bottom the onset of convective instability (t = 400 s) results in a 

decrease in temperature due to cooling by the flow field whereas at the top the temperature increases 

beyond the onset of oscillations (t = 800 s). These trends are representative of the behavior of the flow 

and temperature fields for the prototype for Ra = 1.0×109

. 

We now consider the behavior of the flow and temperature fields for a model near the threshold of 

instability at Ra = 1.26×106 as well as a case for Ra = 1.0×109 in order to compare with the prototype. The 

application of dynamical similarity (ΛM = ΛP) using a model with smaller geometric length scale, showing 

time asymptotic behavior near the top of the cavity at (0.15, 0.90), is shown in Figure 11. Since the cavity 

size of the model is smaller, the equivalent amount of heat input to the bottom diffuses much faster as 

compared to the prototype, therefore the instability occurs much sooner. Near the threshold of instability 

at Ra = 1.26×106 there is oscillation in both the temperature and flow fields; the power spectrum of the 

flow field shows oscillations with frequency of f = 0.002 Hz and harmonics of the frequency, 2f, 3f, 4f

(4f = 0.008 Hz) which spans a range of period from 125 to 500 s. Comparison of the model, 

corresponding to time asymptotic state, and prototype in Figure 10 at location (0.15, 0.9) shows that the 

model exhibits aperiodic behavior as time increases for Ra = 1.0×109

. The time asymptotic behavior of 

the velocity component shows intermittent periodic bursts which contain a broadband frequency content 

extending to 0.4 Hz as indicated by the power spectrum Pu; this indicates a transition to chaos. We 

compare a sub-interval of the time history to the prototype in Figure 10 at (0.15, 0.90) and the remarkable 

results show similar behavior for both the temperature and flow fields. The initial response time of the 

model occurs near t = 10 s whereas the prototype corresponds to t = 200 s, and the onset of convective 

instability occurs near t = 60 s whereas the prototype indicates t = 700 s. Similar to the prototype, the 

power spectrum indicates periodic frequency of 0.02 Hz, corresponding to a period of 50 s, with a higher 

harmonic. Such dynamical similar results indicate that a model can be used to infer the dynamical 

characteristics of a prototype. This finding is practical as it represents an alternate means to study the 

behavior of a prototype which requires significantly less computational resources; for example a typical

prototype that requires 800 hr to run may only require 40 hr for a model. The difference between the 

model and prototype is the characteristic time for the onset of instability, the model requires much shorter 

time, hence the savings in computational time.
NASA/TM—2010-216789 26

The effect of aspect ratio on the prototype for Ar = 1.6 and Ar = 3.24 is contrasted for long and short 

times in Figure 12. For Ar = 1.6, we show a case for Ra = 1.0×109 corresponding to subcooling at the top 

boundary (ft = –0.1) which is relatively low. The time history of the temperature and flow fields near the 

bottom of the enclosure at (0.0025, 0.0017) is similar to that in Figure 10 shown at (0.002, 0.002). Thus low

subcooling on the top boundary does not have a significant effect on the dynamical characteristics of the 

flow field. We show a full scale simulation of the prototype for Ar = 3.24 in which the grid density is 

increased from 400×600 (Ar = 1.6) to 1000×1000; this is demonstrated for the baseline condition in which 

the same amount of heat flux is absorbed at both the bottom and top boundary (fb = 1, ft = 1) for short times. 

The results show that it takes on the order of 300 s for the onset of instability; this corresponds to 

approximately 800 hr of computational time. To compute the time asymptotic behavior requires significant 

resources to reach the asymptotic state. The wall temperature at the bottom (0.5, 0.0) shows a monotonic 

increase without fluctuations; this occurs because for the short time considered conduction effects are 

dominant. In comparison to the parametric range of the experiment (Ref. 1), the upper limit of imposed heat 

flow at the bottom boundary is 104 W which corresponds to approximately qb

” = 3.49×106 erg/cm2

-s. In 

order to find the baseline condition to estimate the heat flux absorbed from the top and bottom boundary we 

assumed a heat flow of approximately 1 mW, from background environment, which corresponds to a heat 

flux of qb

” = 3.49×101 erg/cm2

-s. The scaling of the temperature magnitude using Equation (13), indicates 

that the heat flux absorbed at the boundary is greater than the assumed value of qb

” = 3.49×101 erg/cm2

-s and 

lies in the range of 3.49×101 erg/cm2

-s < qb

” < 3.49×106 erg/cm2

-s. Using a combination of dynamical 
NASA/TM—2010-216789 27

similarity and full scale simulation of the prototype, it remains to investigate the effect of higher heat flux 

absorption at the boundaries. However, the parametric range considered for low heat flux condition at the 

top and bottom boundaries provide a foundation to investigate heat entrapment effects at higher Rayleigh 

numbers.

Summary and Conclusions

We address heat entrapment effects with application to liquid acquisition device (LAD) inside 

cryogenic storage tanks. We use dynamical similarity for the model along with direct numerical 

simulation of the prototype with much larger geometric length scale to address the global dynamics of the 

flow field and local dynamics at fixed points. From scaling analysis, we show that the parametric space 

depends on six parameters namely (Ra, Pr, Ar, fs, fb, ft) the Rayleigh and Prandtl numbers, the aspect ratio, 

and the corresponding heat flux ratios for the side, bottom, and top boundaries. Due to the large 

parametric set involved, we show a selected parametric range based on typical experimental conditions 

applicable to LAD. The parametric set entails Ra between 1 and 1.26×109

, Ar between 1 and 3.24,

Pr = 1.03 and 2.27, (fs, fb) between 0 and 1, and ft between –0.1 to 1. We address particularly the effect of 

applied heat flux at the boundaries to establish baseline condition for a non-equilibrium cryogenic fluid 

inside an enclosure absorbing heat due to background laboratory environment. This serves to establish the 

conditions likely to occur at the top and side boundaries which can be used for simulation with known 

input heat flux to the bottom boundary.
NASA/TM—2010-216789 28

The results show that for typical conditions convection is driven via the interaction between natural 

convection driven by heat flux absorbed along the vertical walls and convective instability due to heat 

flux absorbed along the horizontal bottom boundary. Heat flux absorbed at the top boundary serves as a 

stabilizing mechanism. The predominant mechanism for convection is convective instability for which the 

longest wavelength occurs for Ra on the order of 105 which is near the threshold of instability. The 

maximum number of modes that occur is six for Ra of order 107

; these modes segregate toward the wall 

as Ra increases to 1.26×109

. The wavelength of the horizontal modes decreases as Ra increases which 

governs the length scale that needs to be resolved for computations of high Rayleigh numbers. For 

Rayleigh numbers near the threshold of instability the flow field remains symmetric, however for Ra on 

the order 107 there is a transition to asymmetric flow. 

Direct numerical simulation for large geometric length scale prototype for Ra = 1.0×109 shows that 

the global dynamics of the flow field evolve in three stages. There occurs a transition from short 

wavelength instability to nonlinear growth of thermals; the flow field propagates to the top of the cavity 

via a repetitive sequence of destabilization which gives rise to secondary and tertiary modes, and a 

travelling wave mode occurs prior to transition to asymmetry. The effect of increasing the aspect ratio is 

to increase the number of modes in the vertical direction. Due to the slow diffusion of heat for large 

aspect ratio enclosures, the onset of instability occurs much later in comparison to a model with smaller 

geometric dimensions and computation of asymptotic dynamics of the flow field becomes prohibitive for 

the prototype. We use dynamical similarity to show that similar events can be obtained with a model 

which facilitates the attainment of asymptotic states. The time history at specific points in the flow field 

show that for Ra on the order of 106 the temperature and flow fields oscillate periodically, however for Ra 

on the order of 109 there is aperiodic behavior with intermittency in the flow field which is indicative of a 

transition to a chaotic state. Scaling of the temperature magnitude for the parametric range of the model 

and prototype indicate that there is a higher heat flux absorption on the boundaries than the nominal value 

of 34.9 erg/cm2

-s or approximately 1 mW of background ambient heat flow. However, the parametic 

study shed light on the basic flow field dynamics associated with heat entrapment and can be used in 

future studies to address effect of increased heat flux at the boundaries. One of the highlights of the 

parametric study is the finding that transition to asymmetry occurs for Rayleigh numbers relevant to heat 

entrapment, this implies that axi-symmetric models are to be used with caution, since they are applicable 

to a very narrow parametric range near the critical Rayleigh number.

References

1. Boslshinskiy, L.G., Hastings, L.J., Statham, G., “Cryogenic Capillary Screen Heat Entrapment,” 

Cryogenics, Vol. 48, 2008, pp. 232-237.

2. Tellep, D.M. and Harper, E.Y., “Approximate Analysis of Propellant Stratification,” AIAA Journal, 

Vol. 1, No. 8, 1963, pp. 1954-1956.

3. Vliet, G.C., Brogan, T.J., Sheppard, T.S., Morse, F.H., and F.L. Hines, “Stratified Layer Flow Model 

– A Numerical Approach To Temperature Stratification In Liquids Contained In Heated Vessels,” 

Aerospace Sciences Meeting, N.Y., N.Y. Jan. 20-22, Preprint No. 64–37, 1964.

4. Vliet, G.C., “Stratification with Bottom Heating,” J. Spacecraft, Vol. 3, No.7, 1966, pp. 1142–1144.

5. Evans, L.B., Matulevicius, E.S., and Pan, P.C., “Reduction of Thermal Stratification by Nonuniform 

Wall Heating,” AIChE Journal, Vol. 17, No. 3, 1971, pp. 752–753.

6. Evans, L.B., and Stefany, N.E., “An Experimental Study of Transient Heat Transfer To Liquids In 

Cylindrical Enclosures,” Chem. Eng. Prog. Symposium Ser. Vol. 62, No. 64, 1966, pp. 209–215. 

7. Barakat, H.Z., Clark, J.A., “Analytical and Experimental Study of The Transient Laminar Natural 

Convection Flows in Partially Filled Liquid Containers,” Proceedings of the Third International Heat 

Transfer Conference, Vol. II, Chicago, Illinois, Aug. 7–12, 1966, Papers 41–80. 

8. Polyzhaev, V.I., “Nonstationary Laminar Heat Convection In A Closed Domain For a Given Heat 

Flux,” Izvestiya Akademii Nauk SSSR, Meckanika Zhidkosti I Gaza, No. 4, 1970, pp. 109–117.
NASA/TM—2010-216789 29

9. Polyzhaev, V.I., “Convective Interaction In a Cylindrical Vessel Partially Filled With Liquid And 

With Heat Supplied To The Side Surface, Free Surface, And Bottom,” Izvestiya Akademii Nauk 

SSSR, Meckanika Zhidkosti Gaza, No. 4, 1972, pp. 81–88.

10. Elder, J.W., “Laminar Free Convection In A Vertical Slot,” J. Fluid Mech., Vol. 23, Part 1, 1965,

pp. 77–98.

11. Elder, J.W., “Turbulent Free Convection In A Vertical Slot,” J. Fluid Mech., Vol. 23, Part 1, 1965, 

pp. 99–111.

12. Seki, N., Fukusako, S., Inaba, H., “Visual Observation of Natural Convection Flow In a Narrow 

Vertical Cavity,” J. Fluid Mech., Vol. 84, Part 4, 1978, pp. 695–704.

13. Akiyama, M., Hwang, G.J., Cheng, K.C., “Experiments On The Onset Of Longitudinal Vortices in 

Laminar Forced Convection Between Horizontal Plates,” J. Heat Transfer, Vol. 93, 1971, pp. 335–

341.

14. Chandrasekhar, S., Hydrodynamic and Hydromagnetic Stability, Oxford University Press, Oxford, 

UK, 1961, pp. 9–71.

15. Boris, J.P., and Book, D.L., “Flux-Corrected Transport I: SHASTA-A Fluid Transport Algorithm 

That Works,” J. Comp. Phys., Vol. 11, Jan. 1973, pp. 139–174.
REPORT DOCUMENTATION PAGE Form Approved 

OMB No. 0704-0188 

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining the 

data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing this 

burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. 

Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB 

control number. 

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS. 

1. REPORT DATE (DD-MM-YYYY)

01-09-2010 

2. REPORT TYPE

Technical Memorandum 

3. DATES COVERED (From - To)

4. TITLE AND SUBTITLE

Heat Entrapment Effects Within Liquid Acquisition Devices 

5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S)

Duval, W.M.B.; Chato, D., J.; Doherty, M., P. 

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

WBS 095240.04.03.03.01.03 

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration 

John H. Glenn Research Center at Lewis Field 

Cleveland, Ohio 44135-3191 

8. PERFORMING ORGANIZATION

 REPORT NUMBER

E-17188 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration 

Washington, DC 20546-0001 

10. SPONSORING/MONITOR'S

 ACRONYM(S)

NASA 

11. SPONSORING/MONITORING

 REPORT NUMBER

NASA/TM-2010-216789 

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified-Unlimited 

Subject Categories: 34, 29, and 31 

Available electronically at http://gltrs.grc.nasa.gov 

This publication is available from the NASA Center for AeroSpace Information, 443-757-5802

13. SUPPLEMENTARY NOTES

14. ABSTRACT

We introduce a model problem to address heat entrapment effects or the local accumulation of thermal energy within liquid acquisition 

devices. We show that the parametric space consists of six parameters, namely the Rayleigh and Prandtl numbers, the aspect ratio, and heat 

flux ratios for the bottom, side, and top boundaries of the enclosure. For the range of Ra considered 1 to 109, beyond Ra on the order of 105, 

convective instability is the dominant mode of convection in comparison to natural convection. The flow field transitions to asymmetric 

modes at Ra on the order of 107. Direct numerical simulation of a large geometric length scale prototype for Ra on the order of 109 shows 

that the flow field evolves from small wavelength instability which gives rise to nonlinear growth of thermals, propagation of the instability 

occurs via growth of secondary and tertiary modes, and a travelling wave mode occurs prior to asymmetry. The effect of a large aspect ratio 

is to increase the number of modes in the vertical direction. Due to the slow diffusion of heat in the prototype, asymptotic states are not 

readily attained, we show that dynamical similarity can be used for a model which allows the attainment of asymptotic states and that 

transition to a chaotic state occurs for Ra on the order of 109 via a broadband power spectrum. These dynamical events show that for the 

baseline condition in which heat is absorbed from background laboratory environment, higher heat flux is absorbed at the top and bottom 

boundaries of the enclosure than a nominal value of 34.9 erg/cm2-s.

15. SUBJECT TERMS

Convective instability; Cryogenics; Buoyancy-induced flow; Natural convection 

16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF

 ABSTRACT

UU 

18. NUMBER

 OF 

 PAGES

35 

19a. NAME OF RESPONSIBLE PERSON

STI Help Desk (email:help@sti.nasa.gov) 

a. REPORT

U 

b. ABSTRACT

U 

c. THIS 

PAGE

U 

19b. TELEPHONE NUMBER (include area code)

443-757-5802 

Standard Form 298 (Rev. 8-98)

Prescribed by ANSI Std. Z39-18