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/REV1

July 2011

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/REV1

July 2011

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://www.sti.nasa.gov

Level of Review: This material has been technically reviewed by technical management. 

Document Change History

NASA/TM—2010-216789/REV1, July 2011

Heat Entrapment Effects Within Liquid Acquisition Devices

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

This printing replaces NASA/TM—2010-216789, September 2010, because extensive changes were made to the document.
NASA/TM—2010-216789/REV1 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-type motion of 

convective modes occur 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/REV1 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 screen 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 for the 

experimental conditions with partitions, the capillary screen meshes were impervious to natural 

convection. This is verified experimentally (Ref. 1) from comparison between data for meshes versus a 

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; 

this implies that the impermeable zero mass flux condition of the solid boundary models the boundary 

condition for the screen meshes. 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, of 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/REV1 3 

(Refs. 8 and 9) who considers 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 from the standpoint 

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 wavelike motion along the wall. The visual observation by Seki et al. (Ref. 12) has 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 the limiting 

case of 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 

Rayleigh-Benard 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 (Refs. 10 to 12) of a 

narrow vertical and horizontal cavity (Ref. 13) 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 

laboratory condition, and also to assess 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 sidewall 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/REV1 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-type mode prior to 

transition to asymmetry. Due to the large geometric length scale of the prototype, attainment of 

asymptotic states can be prohibitive. We employ dynamical similarity using a model with smaller 

geometric dimensions to show that similar dynamical events occur which facilitate 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 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 sidewall, 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 1A(b). 
NASA/TM—2010-216789/REV1 5 

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

”

 = qt

”

 = qb

”

 = 0; this 

condition is not likely met for a low temperature cryogenic liquid inside an enclosure. Thus offequilibrium condition need to satisfy the dynamical motion of the fluid govern 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

æ ö ¶ ¶ ¶ ¶ ¶¶ ç ÷

+ + =- +n + ç ÷

ç ÷ ¶¶¶ ¶ ç¶ ¶ ÷÷ è ø

(2)

2 2

2 2 0

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

æ ö ¶ ¶ ¶ ¶ ¶¶ ç ÷ + + = - + n + - bD ç ÷

ç ÷ ¶¶¶ ¶ ç¶ ¶ ÷÷ è ø

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

2 2

2 2

T T T TT

u v

t xy x y

æ ö ¶ ¶ ¶ ¶¶ ç ÷

+ + =a + ç ÷

ç ÷ ¶¶¶ ç¶ ¶ ÷÷ è ø

    (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, 

( )[1 ] r = r - bD T

A 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 b = - r ¶r ¶T

Superscript ~ denotes dimensional temperature and D=- 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 = =G (7) 

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

"" " 0 0, ,0 , ) / s xs x q q yH q k T

x = = = £ £ =- ¶ ¶

 (8)

"" " , ,0 , ) / s xL s x Lq q y H q k T

x = = = £ £ =- ¶ ¶

 (9)

"" " 0, ,0 , ) / b yb 0 y q q xL q k T

y = = = £ £ =- ¶ ¶

 (10)

"" " , ,0 , ) / t yH t y Hq q x L q k T

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

y y x ¶ ¶ ¶¶ = =- = - ¶ ¶ ¶¶ (12) 
NASA/TM—2010-216789/REV1 7 

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

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

" /

A T T

T

qL k

- =

 

(13)

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

2* 2* *

*2 *2

x y

¶y ¶y + = -z ¶ ¶ (14)

* * * 2* 2* * *

* * * *2 *2 * Pr Pr T u v Ra

t x y xy x

æ öæ ö ¶x ¶x ¶x ¶ x ¶ 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 f T

x = = £ £ =- ¶¶ (17) 

* * *

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

x = = £ £ =- ¶¶ (18) 

* * *

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

y = = £ £ =- ¶¶ (19) 

* * *

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

y = = £ £ =- ¶¶ (20) 

The set of parameters resulting from scaling is defined as, 

"4 " " "

" "" , Pr , , , , o sbt

sbt

q ng L q q q H Ra Ar f f f k L qqqb n

na a = = == = = (21) 

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

fluxes. Thus the parametric set, 

( , , Pr, , , ) L=L 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, 

respectively. The side conditions 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 
NASA/TM—2010-216789/REV1 8 

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 lead to convective stability in 

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

Figure 2(a). These trends show that stability of the flow field can be obtained, below the critical Rayleigh 

number, provided that there is a perfectly insulated sidewall fs = 0 for cases in which: (1) there is bottom 

heating only (fb = 1, ft = 0), (2) top heating only (fb = 0, ft = 1), and (3) 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 of an 

experimental system (Ref. 1) since there is finite heat leakage along the sidewall, i.e., fs  0, regardless of 

precautionary measures. Note that absolute stability occurs for case 2 regardless of the Rayleigh number 

magnitude, however, at a critical Rayleigh number cases 1 and 3 become unstable marked by the onset of 

convection. 

Perturbation of Stability 

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

in Figure 2(b) for Ra = 1.26106

 in order to isolate simultaneous effects of bottom and top heating, and to 

show the interaction of natural convection and convective instability. For a conservative heat leakage on 

the order of 1 μW projected on the sidewall surface of a cylinder fs = 210–6 (qs”

 = 0.001 erg/cm2-s), the 

stable base state (fb = 1, ft = 0) 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 (fb = 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), with a side heat leakage of fs = 210–6, we show that the top 

stratified layer (fb = 0, ft = 1, fs = 0) remains stable for the same heat leakage; however horizontal density 

gradient gives rise to a natural convection base flow with a symmetric double cell. These two cases 

illustrate that even though the top stratified layer is stable for all Rayleigh numbers when the sidewall is 

insulated fs = 0, the effect of finite side heat leakage is to drive a natural convection base flow. 
NASA/TM—2010-216789/REV1 9 
NASA/TM—2010-216789/REV1 10 

Effect of Side Heating 

Since side heat leakage is an effective source to excite instability and drive a natural convective base 

flow, 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 in Figure 2(c). These two cases 

contrast heat flow on the order of 1 μW (qs” = 0.001 erg/cm2

-s) and 10 mW (qs” = 500 erg/cm2-s) 

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 events that occur in applications to address the effect of 

thermal stratification for partially filled cryogenic storage tanks (Refs. 8 and 9). 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 heat leakage toward the lower limit 

of Ra = 2.5. 

Perturbation of Double Layer 

A typical heating condition with finite sidewall heating is illustrated in Figure 2(d) 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

 based on case 3. The destabilization of the conduction dominated 

regime (fb = 1, ft = 1, fs = 0), near the critical Rayleigh number, is contrasted 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 the growth of the stratified layer at the top region for early time t = 360 s remains stable. 

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 characteristic inside the 

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

of the 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; there is 

saturation in mode numbers between Ra = 1.26107

 and Ra = 1.26109. The growth of modes segregates 

near the walls as Ra increases and the wavelength between the cells decreases. This trend indicates that 

high Rayleigh number flows require resolution of small scales for computational accuracy. 

 
NASA/TM—2010-216789/REV1 11 
NASA/TM—2010-216789/REV1 12 

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

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

Rayleigh number considered (1.26109

), l = 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 (Fig. 2(c)), in this convective instability dominated situation (Fig. 4), the 

maximum temperature occurs near the mid-center of the cavity due to bottom heating since the sidewalls 

are nearly insulated (fs = 210–5). 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) selforganizes 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/REV1 13 

 
NASA/TM—2010-216789/REV1 14 
NASA/TM—2010-216789/REV1 15 

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

absorption occurs 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 partitioned co-rotating cells with the middle cell 

(tertiary mode) counter-rotating. 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/REV1 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. 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/REV1 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 for Ra > 107

. 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

 and ft < 1; 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, thus a reverse trend is predicted with an asymptotic 

approach to a lower temperature in comparison to ft = 1. If heat entrapment is based on the condition ft  1, 

indicating finite heat leakage at the divider boundary then heat entrapment is less effective for Ra < 107

when ft = 1, whereas for Ra > 107

 it becomes more effective. The trend for Ra > 107 is in agreement with the 

trend in experimental measurements of temperature increase near the bottom divider (Ref. 1), since Ra > 107 

for the experimental system. 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 

experiments (Ref. 1) 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 3B. For example the horizontal length scale of the prototype can be four times greater 
NASA/TM—2010-216789/REV1 18 

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 h using 4 Central Processing Units (CPU’s) now requires between 300 to 800 h using 

20 to 32 CPU’s for the prototype. Even though, direct numerical simulation of the prototype is the most 

exact, dynamical 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) a travelling wave-type mode which eventually 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) allows the resolution of small length scales 

such as the fine structures during 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 subsequently 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-type mode event that occurs as shown in Figure 7C. The travelling wave-type 

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 travelling wave-type mode propagates 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 further up the cavity vertical walls. 

Effect of Aspect Ratio 

The effect of increasing the height of the cavity to actual prototype dimension 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 a 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 
NASA/TM—2010-216789/REV1 19 

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/REV1 20 
NASA/TM—2010-216789/REV1 21 

 
NASA/TM—2010-216789/REV1 22 

Application of Dynamical Similarity 

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 L =L 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 Figures 7A, B, 

C 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 its 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/REV1 23 
NASA/TM—2010-216789/REV1 24 
NASA/TM—2010-216789/REV1 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–5), 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 aperiodically in the core region. 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.004 Hz with slight harmonics at (0.15, 0.90); 

this corresponds to a long period of 250 s. Note that near the top (0.002, 0.998) conduction effects 

dominate 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 

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.001 Hz and harmonics of the frequency, 2f, 3f, 4f

(4f = 0.004 Hz) which spans a range of period from 250 to 1000 s. Comparison of the model, showing a 

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 dynamical state 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 response, however the frequency (0.04 Hz) is ten times 

larger with a corresponding shorter period of 25 s, and higher harmonics indicating nonlinear effects. 

Such dynamical similar results indicate that a model can be used to infer the dynamical characteristics of 

a prototype approximately. 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 h to run may only require 40 h for a model. The major 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/REV1 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 h 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), for 

comparison with experimental measurements (Ref. 1) for the baseline condition of no screen, 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/REV1 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, provides 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 and temperature fields 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 a baseline condition for a nonequilibrium 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/REV1 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. These stages consist of a transition from short 

wavelength instability to nonlinear growth of thermals, propagation of the flow field to the top of the 

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

travelling wave-type motion of convective modes which occur 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 challenging 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 shows 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, indicates 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. Bolshinskiy, 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/REV1 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-07-2011 

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-1 

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-REV1 

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

This printing replaces NASA/TM--2010-216789, September 2010. 

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-type motion of convective modes occur 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