Walter M.B. Duval, David J. Chato, and Michael P. Doherty

Glenn Research Center, Cleveland, Ohio

Transient Convection Due to Imposed Heat Flux: 

Application to Liquid-Acquisition Devices

NASA/TM—2014-217442

June 2014

AIAA–2011–1321
NASA STI Program . . . in Profile

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA Scientific 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 Officer. 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 significant phase 

of research that present the results of NASA 

programs and include extensive data or theoretical 

analysis. Includes compilations of significant 

scientific 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. Scientific 

and technical findings 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. Scientific and 

technical findings by NASA-sponsored 

contractors and grantees.

• CONFERENCE PUBLICATION. Collected 

papers from scientific and technical 

conferences, symposia, seminars, or other 

meetings sponsored or cosponsored by NASA.

• SPECIAL PUBLICATION. Scientific, 

technical, or historical information from 

NASA programs, projects, and missions, often 

concerned with subjects having substantial 

public interest.

• TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientific 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 to help@sti.nasa.gov

• Fax your question to the NASA STI 

Information Desk at 443–757–5803

• Phone the NASA STI Information Desk at

443–757–5802

• Write to:

 STI Information Desk

 NASA Center for AeroSpace Information

 7115 Standard Drive

 Hanover, MD 21076–1320
Walter M.B. Duval, David J. Chato, and Michael P. Doherty

Glenn Research Center, Cleveland, Ohio

Transient Convection Due to Imposed Heat Flux: 

Application to Liquid-Acquisition Devices

NASA/TM—2014-217442

June 2014

AIAA–2011–1321

National Aeronautics and

Space Administration

Glenn Research Center

Cleveland, Ohio 44135

Prepared for the

49th Aerospace Sciences Meeting

sponsored by the American Institute of Aeronautics and Astronautics

Orlando, Florida, January 4–7, 2011
Acknowledgments

This work was performed in support of the Cryogenic Fluid Management program at the NASA Glenn Research Center. 

The computations were performed using the Columbia supercomputer (Silicon Graphics International Corp. Altix) 

at the NASA Ames Research Center’s NASA Advanced Supercomputing facility.

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. 
NASA/TM—2014-217442 1

Transient Convection Due to Imposed Heat Flux: 

Application to Liquid-Acquisition Devices

Walter M.B. Duval, David J. Chato, and Michael P. Doherty

National Aeronautics and Space Administration

Glenn Research Center

Cleveland, Ohio 44135

Summary

A model problem is considered that addresses the effect of heat load from an ambient laboratory 

environment on the temperature rise of liquid nitrogen inside an enclosure. This model has applications to 

liquid acquisition devices inside the cryogenic storage tanks used to transport vapor-free propellant to the 

main engine. We show that heat loads from Q = 0.001 to 10 W, with corresponding Rayleigh numbers 

from Ra = 109 to 1013, yield a range of unsteady convective states and temperature rise in the liquid. The 

results show that Q = 1 to 10 W (Ra = 1012 to 1013) yield temperature distributions along the enclosure 

height that are similar in trend to experimental measurements. Unsteady convection, which shows selfsimilarity in its planforms, is predicted for the range of heat-load conditions. The onset of convection 

occurs from a free-convection-dominated base flow that becomes unstable against convective instability 

generated at the bottom of the enclosure while the top of the enclosure is convectively stable. A number 

of modes are generated with small-scale thermals at the bottom of the enclosure in which the flow selforganizes into two symmetric modes prior to the onset of the propagation of the instability. These 

symmetric vertical modes transition to asymmetric modes that propagate as a traveling-wave-type motion 

of convective modes and are representative of the asymptotic convective state of the flow field. Intense 

vorticity production is created in the core of the flow field due to the fact that there is shear instability 

between the vertical and horizontal modes. For the higher Rayleigh numbers, 1012 to 1013, there is a 

transition from a stationary to a nonstationary response time signal of the flow and temperature fields with 

a mean value that increases with time over various time bands and regions of the enclosure. 

1.0 Introduction

The investigation of the effect of heat flux applied at the bottom boundary of an enclosure, which

contains cryogenic liquid nitrogen, to determine its temperature response requires knowledge of the 

dynamical state of the system when there is no applied heat flux; the enclosure is located inside an ambient 

laboratory environment. The no-applied-heat-flux case is necessary since it represents the effect of the heat 

flow from the background laboratory environment on the thermodynamic state of the liquid inside the 

enclosure. As part of a broader program to investigate heat entrapment effects in liquid acquisition devices,

Bolshinskiy et al. (Refs. 1 and 2) conducted an experiment to measure the effect of background heat flow 

from a laboratory environment on the temperature rise of liquid nitrogen inside an enclosed container 

(Dewar). The results showed that, when no heat flux was applied directly at the bottom boundary, temperature rise occurred at various positions along the vertical height of the enclosure. A model problem is 

addressed to obtain insight into the meaning of the temperature measurements (in no-screen, no-heater 

conditions, Ref. 1) for the sole condition of heat flow to the system from the background laboratory 

environment.

To determine a suitable range of parameters, we estimated a range of typical heat loads caused by the 

thermal background radiation, which yielded a conservative heat load of Q = 0.001 W to a typical maximum value of Q = 10 W. We approximated the problem as a rectangular cross section of a cylinder; this 

was convenient since, at the lowest heat load (Q = 0.001 W), the Rayleigh number Ra is on the order of 

109 and flow asymmetry occurs in this system for Ra on the order of 107

, as shown by Duval, Chato, and 

Doherty (Ref. 3). This implies a full transient calculation for a three-dimensional cylinder. However, a 

rectangular model provides an approximation to the system that allows the transient dynamics to be 
NASA/TM—2014-217442 2

investigated beyond asymmetry. The dynamics of the system for Ra = 1 to 108 has been considered for an 

aspect ratio of Ar = 1 (Ref. 3), and a short-time solution for Ra = 109 is presented for Ar = 3.24 (Ref. 3). 

Herein, we consider the asymptotic dynamics for Ra = 109 and Ar = 3.24, and we extend the calculations 

to Ra = 1013 for the highest heat-load condition. We show that at the highest heat load, the model predicts

temperature magnitude comparable to that in the experiments (Refs. 1 and 2). In addition, the model 

shows the dynamical state of the accompanying flow field, which was not accessible during the 

experiments since only temperatures inside the system were measured (Refs. 1 and 2).

In analogy to our model problem of a cryogenic liquid inside an enclosure that has a large height and 

that absorbs heat flux predominantly from its bottom and top boundaries, there have been a number of

investigations on the convective instability inside enclosures heated from below with prescribed temperature 

boundary conditions on the bottom and top boundaries, which lend insight into our model problem. 

Townsend (Ref. 4) showed that convection near a heated horizontal surface depends on the layer near the 

heated surface (independent of the geometrical height of the enclosure) and is determined by the heat flux, 

viscosity, and conductivity of the fluid. Townsend described the formation of thermals on the heated surface 

as rising columns of hot air from the edge of the conduction layer. A phenomenological theory was put forth 

by Howard (Ref. 5) who postulated that convection consists of a conductive phase, during which the 

conduction boundary layer is formed, followed by a short time interval, during which the thermals break off. 

Convective instability, or the formation of thermals, occurs after a time t

o that is shorter than the thermal 

diffusion time H2

/α, where H is the enclosure height and α is the thermal diffusivity. However, the thermals 

break away at a time that is short in comparison to the conductive period t

o

. In addition, Howard (Ref. 5)

showed that the mean temperature profile, based on the idea of marginal stability of the mean flow, depends

on the thermal boundary layer thickness δ when it becomes unstable (in which the local Rayleigh number is 

based on δ) and shows reasonable agreement with experimental data.

A number of two-dimensional numerical and experimental investigations of convective instability of 

flows inside enclosures that are heated from below consider the short-time evolution of the instability (Refs. 6 

to 9) instead of the statistically steady regime (Ref. 10) for long times. Elder (Ref. 6) showed that a HeleShaw cell can be used to simulate, experimentally and computationally, two-dimensional flow in a porous 

medium within the context of flow in a homogeneous horizontal slab that is heated from below. When the 

approximation of a Hele-Shaw cell was used, numerical prediction of unsteady convection (Ref. 7) inside an 

enclosure in which the width of the heating zone at the base can be adjusted showed that the first motion is a 

set of eddies growing on opposite ends of the local heating zones. As the instability evolves, additional eddies 

are formed in the central region of the heating zone. For a long heating zone that is approximately the horizontal length of the slab, the initial motion of a pair of end cells followed by successive growth of a string of 

cells above the heater is similar to the case of a short heating zone. These events are phenomenologically 

similar to the short-time prediction of our model for the initiation of convective instability.

For a gravitationally unstable layer of viscous fluid initially hot below and cold above a horizontal 

plane, Elder (Ref. 8) showed that the scale of the motion is set by the thickness of the interface δ and is 

independent of the vertical distance H of the horizontal plane in agreement with Townsend’s (Ref. 4)

findings. The amplification process of the disturbance in the interface is shown to be similar for the flow 

through a porous medium and for a viscous fluid inside a Hele-Shaw cell with impermeable boundaries. 

Experimental observation (Ref. 8) confirming numerical prediction shows that the instability develops as an 

array of blobs above a lower surface in a proto-sublayer. The disturbance remains embedded in the protosublayer until it reaches a finite amplitude. Similar events for short-time dynamics are shown to occur in our 

model for Rayleigh numbers on the order of 1012. The effect of various types of disturbance heat sources at 

the bottom boundary—namely, random, spatially varying, and constant—show similar developments of the 

instability (Ref. 8). This implies that the instability can be induced with a constant heat source as shown in 

our computational model. The development of the proto-sublayer is shown to be similar to the phenomenological description of the sublayer proposed by Howard (Ref. 5): that the mean temperature profile depends

on the boundary layer thickness of the sublayer in which conduction effects dominate.

Elder (Ref. 9) shows that using one-dimensional mean field equations to model the flow development of 

the proto-sublayer can simulate two-dimensional numerical models. The spatial means are taken over a 
NASA/TM—2014-217442 3

horizontal plane in the mean field equations, and the simulation is shown for Rayleigh numbers on the order 

of 105

. Even though Elder’s numerical simulation (Refs. 7 and 8) captures the growth of cells in the flow 

field as the convective instability evolves, the temperature field shows that thermals are roughly represented 

by blobs; this is also shown experimentally for two-dimensional flows simulated with a Hele-Shaw cell

configuration. The resolution of mushroomlike thermals exhibiting a blunted nearly hemispherical cap was 

captured in the experiments by Sparrow, Husar, and Goldstein (Ref. 11), which used water inside a threedimensional enclosure heated along a horizontal surface. They showed that thermals are generated 

periodically in time at fixed sites and that the breakup Rayleigh number based on the boundary layer 

thickness of the conduction layer is a constant, which validates a prediction by Howard (Ref. 5). Even 

though, for the short-time development of convective instability (in particular the proto-sublayer (Ref. 8)), 

the preceding works show trends that are similar to those of our model, for the long-time-scale asymptotic 

dynamics, we show herein that the dynamical state of the system is a traveling-wave-type motion of 

convective modes, which damps the evolution of thermals near the bottom boundary.

In the following sections, the formulation of the basic physics of the problem and discussion of the 

numerical solution are presented. Then, the global dynamics of the flow field—showing the growth of 

thermals, the transition to asymmetry, and the traveling-wave-type motion of the convective modes—are 

presented. Next, the flow field is probed at various fixed points to show time histories of the flow and 

temperature fields. Finally, we identify the range of heat load from the ambient laboratory environment,

showing predicted temperature magnitudes that are comparable to those of the experiments.

2.0 Formulation

The model, which addresses the effect of the external heat load Q from an ambient laboratory 

environment at temperature T∞ = 293 K on the temperature rise of saturated liquid nitrogen, initially at

TA = 77 K, inside an enclosure is shown in Figure 1. The heat load Q transmits constant heat flux q'' on 

the boundaries of the enclosure with equal heat flux absorption on the bottom and top (qb'' = qt''), whereas

the sidewall heat flux qs'' is assumed to be small, qs'' << qb'', since insulation was used at the side 

boundaries in these experiments (Refs. 1 and 2) to investigate heat entrapment effects inside liquid 

acquisition devices. The figure shows the locations of points inside the enclosure (xp, yp) that allow 

Figure 1.—Physical description of cryogenic fluid (ambient temperature, TA, 

and density, ρA) inside an enclosure absorbing surface heat flux, at the 

bottom, sides, and top (qb", qs", and qt"), from ambient environment heat 

flow, Q. Location of points in the flow field are denoted as (xp , yp).
NASA/TM—2014-217442 4

probing of the time history of the flow and temperature fields similar to the experimental measurement 

points. The maximum heat load absorbed from the ambient owing to background thermal radiation can be 

estimated to be on the order of 10 W, whereas a minimum conservative value is assumed to be 0.001 W. 

The absorption of an infinitesimal heat load Q through the boundaries of the large enclosure with height

H = 61.9 cm, approximately 3 times its horizontal length, L = 19.1 cm, drives the system off equilibrium.

A deviation from equilibrium conditions will cause an increase in temperature given by

DT 2T

Dt

=α∇ (1) 

which can give rise to fluid motion owing to its coupling with the body force (gy = ngo) through the 

density field. The body force is directed downward vertically, in which n represents the various 

gravitational levels of the Earth (n = 1), Moon (n = 1/6), or Mars (n = 1/3), and go is the gravitational 

acceleration on Earth. For an incompressible Boussinesq fluid (where V



is the velocity vector field),

∇⋅ = V 0  (2)

When the density ρ is independent of pressure but is a function of temperature, the density can be given 

as an equation of state:

ρ = ρ(TA)[1 – β∆T] (3)

where β = –(1/ρ)(∂ρ/∂T)p is the coefficient of thermal expansivity, the subscript p indicates pressure, and

∆T = T – TA is the temperature difference. The motion of the fluid satisfies

2 go

DV p V Tn

Dt

∇ = − + ν∇ − β∆

ρ

  (4) 

where t is time and p is pressure.

The dynamical equations of motion (Eqs. (1), (2), and (4)) satisfy the initial condition that liquid 

nitrogen is at saturation:

t = 0, T(x, y, 0) = TA (5)

as well as the no-slip boundary condition along the boundary Γ of the enclosure,

t > 0, V = Γ 0 on  (6)

Prescribed heat flux conditions are imposed along the bottom and top boundaries Γ of the enclosure:

0

0, , , b b

y

T q y x qq y =

′′ = ∀ = =− ′′ ′′ ∂

∂ κ (7) 

,, , t t

y H

T q yHx q q y =

′′ = ∀ = =− ′′ ′′ ∂

∂ κ (8)

and along the left and right walls,

0

0, , , s s x

T q x y qq x =

′′ = ∀ = =− ′′ ′′ ∂

∂ κ (9) 
NASA/TM—2014-217442 5

,, , s s x L

T q xLy q q x =

′′ = ∀ = =− ′′ ′′ ∂

∂ κ (10)

where κ is the thermal conductivity.

Since the top heat flux qt

′′ has been shown to stabilize the flow (Ref. 3), fluid motion is generated by 

either natural convection due to the wall heat flux qs

′′ or by convective instability due to the bottom heat 

flux . qb

′′ Since there is a threshold value on convective instability determined by qb

′′ prior to its onset, a 

natural convection base flow is always generated owing to the finite heat absorption s q′′ at the wall. This 

base flow is independent of the heat flux absorbed at the top and bottom boundaries.

3.0 Scaling Analysis and Computational Method

For a computational solution of the dynamical equations of motion (Eqs. (1), (2), and (4)), the vector 

momentum equation (Eq. (4)) can be transformed to a scalar equation for the vorticity field for twodimensional flow. A single vorticity component in the z direction is defined as

v u

x y

∂ ∂ ξ= − ∂ ∂ (11)

where v and u are the vertical and horizontal components of velocity, respectively. The resulting vorticity 

equation—when scaled using the characteristic length, time, and velocity (L, L2

/α, and α/L)—becomes

* * * 2* 2* * * *

* * * *2 *2 * v

T

u Pr Ra Pr

t x y xy x

∂ξ ∂ξ ∂ξ ∂ ξ ∂ ξ ∂   

+ + = ++    ∂ ∂ ∂ ∂∂ ∂    (12)

which implies that vorticity is generated by the horizontal temperature gradient that causes its advection 

and diffusion. Here, Pr is the Prandtl number and Ra is the Rayleigh number; a superscript asterisk (*) is 

used to denote a dimensionless quantity. The continuity equation (Eq. (2)) is satisfied identically when the 

velocity components are expressed in terms of the gradient of the scalar stream function ψ in which

u = ∂ψ/∂y and v = – ∂ψ/∂x; this transformation, when applied to the vorticity component equation

(Eq. (11)), yields Poisson’s equation:

2* 2* *

*2 *2 x y

∂ψ ∂ψ + = −ξ ∂ ∂ (13) 

which is used to solve for the flow field (u, v). The flow field is generated by the buoyancy source term in 

the vorticity field equation from the diffusion of heat along the boundary of the enclosure communicated 

through the density field. The buoyancy source term is coupled to the temperature field given by its 

dimensionless form:

* * * 2* 2* * *

* * * *2 *2 v

T T T TT

u

t x y xy

∂ ∂ ∂ ∂∂   

+ + =+    ∂ ∂ ∂ ∂∂    (14)

obtained using the following temperature scale:

* T TA T

q L

− = ′′ κ

(15)
NASA/TM—2014-217442 6

In addition to Ra and Pr in the vorticity equation (Eq. (12)), which are defined as

4 go and qn L Ra Pr β ν ′′ = = κνα α

(16a,b)

(where ν is the kinematic viscosity), the remaining parameters of the problem are obtained by nondimensionalizing the boundary conditions. This produces the heat flux ratios for the bottom and top boundaries,

fb and ft:

*

* * * 0

0, b

y

y f T

y =

= ∂ = − ∂ (17)

and

*

* * * , t

y Ar

y Ar T f y == ∂ = − ∂ (18) 

and the heat flux ratios of the sidewall boundaries, fs:

*

* * * 0

0, s x

x f T

x =

= ∂ = − ∂ (19)

and

*

* * * 1

1, s xx f T

x =

= ∂ = − ∂ (20) 

An additional four parameters are obtained. These are the aspect ratio Ar and dimensionless heat fluxes fb, 

ft, and fs defined as

H Ar

L = (21a)

b

b

q f

q

′′ = ′′ (21b)

t t

q f

q

′′ = ′′ (21c)

s s

q f

q

′′ = ′′ (21d)

The heat flux scale q" used in Equations (21b), (21c), and (21d) depends on the dominant prescribed 

heat flux at the enclosure boundary that causes convection. For this problem, it is the bottom heat flux qb

′′

that drives convective instability since qs

′′ << . qb

′′ In general, the parametric set Λ consists of six 

parameters:

Λ = Λ(Ra, Ar, Pr, fs, fb, ft) (22)

To simplify the parametric set, we consider a specific cryogenic fluid with fixed geometry such that 

Pr and Ar are constant. In addition, since we assume very low heat leakage at the sidewalls, fs is also 
NASA/TM—2014-217442 7

constant. The heat load from the ambient environment affects Ra, fb, and ft. Since fb = ft and the use of the 

heat flux scale b q′′ implies that fb = 1 and ft = 1, variation of the heat flux appears in fs. Thus, the reduced

parametric set Λ' becomes

Λ' = Λ'(Ra, fs) (23)

This implies that the effect of heat load from the ambient laboratory environment (ngo = 1go) on 

the temperature rise of liquid nitrogen inside an enclosure can be characterized by varying only two 

parameters: the Rayleigh number and the sidewall heat flux ratio. This simplification yields a reasonable 

set of parameters for computational parametric study.

To obtain a computational solution of the dynamical equations of motion for the reduced parametric 

set, we employ finite difference methods. First, we use a direct solver approach to solve Poisson’s 

equation (Eq. (13)) using a matrix inversion that yields the flow field. The input of the flow field is used 

to solve for the temperature field in Equation (14) via a flux-corrected transport method used to resolve 

small-scale thermals. Finally, since the horizontal temperature gradient drives vorticity production 

through buoyancy, we solve the vorticity equation (Eq. (12)). 

This direct solver approach requires the solution of Poisson’s equation for the first time step only. 

Subsequent time steps use the vorticity field (Eq. (12)) as input to the temperature field to obtain a timemarching solution of the flow and temperature fields. Since we are interested in time-accurate solutions, 

we use a third-order Adams-Bashforth time-discretization scheme with a typical time step of approximately one-thousandth of a second: ∆t = 0.00125 s. Such a small time step allows an asymptotic solution 

of the flow field to be computed. Since we are considering relatively high Rayleigh numbers, 109 to 1013, 

the smallest grid size employed is 1000 by 1000. Owing to the large vertical height of the enclosure, 

asymptotic dynamics of the flow field require long-time computations, which can take on the order of 500 

to 1000 hr on a Silicon Graphics International Corporation (SGI) Altix supercomputer system. Memory 

requirements for each case range from a low of 12 CPU to a midrange of 40 CPU and a high of 100 CPU 

on supercomputers with capacities of 508, 1024, and 2048 CPUs (1 CPU is approximately 2 GB). Thus,

substantial computational resources are required to obtain solutions for the reduced parametric set.

4.0 Numerical Results

4.1 Global Dynamics of the Flow Field for High Rayleigh Numbers

Owing to the large geometric length scales of the enclosure (L = 19.1 cm and H = 61.9 cm), computation at the lowest Rayleigh number, for the most conservative estimate of the heat load from the 

background laboratory environment, is on the order of 109

. For the range of Rayleigh numbers (109 to 

1013) considered to investigate the effect of background heat load Q from 0.001 to 10 W, we show that 

departure from thermodynamic equilibrium ensues from convective instability. The instability propagates 

from the bottom to the top of the enclosure and asymptotes to a traveling-wave-type motion of convective

modes toward dynamical equilibrium. The global dynamics of the flow field indicate self-similarity of the 

flow as an asymptotic state is approached. However, the approach to the asymptotic state depends on the 

magnitude of the Rayleigh number.

For the most conservative estimate of background heat load (Q = 0.001 W, which corresponds to 

Ra = 109

), Figure 2 shows a typical base flow from which convective instability ensues. In response to heat 

load from the ambient, a heat flux of 1 2 3.49 10 erg/cm -s q q b t ′′ ′′ = = × is transmitted through the top and bottom surfaces, and careful insulation of the sidewalls allows for a minimal heat flux of 2 0.001erg/cm -s. qs

′′ =

Throughout the parametric variation (the increase in background heat load), the absorption of heat flux 

through the sidewall is maintained constant at 2 0.001erg/cm -s. qs

′′ = Figure 2 shows that prior to convective 

instability, the sidewall heat flux fs causes natural-convection-dominated base flow with two large counterrotating cells. This typical base flow exists for all Rayleigh numbers. The temperature field shows the 

conduction boundary layer δ built up at the top and bottom of the enclosure owing to heat flux ft and fb. This 
NASA/TM—2014-217442 8

is similar to the conduction sublayer in the phenomenological theory of Howard (Ref. 5) and the 

proto-sublayer of Elder (Ref. 8).

Figure 3 shows the transition by which the base flow field becomes unstable against convective 

instability. The instability of the thermal boundary layer thickness δ is evident at t = 300 s and shows that 

the disturbance remains embedded in the proto-sublayer as pointed out by Elder (Ref. 8). Due to the fact that 

heat flux is absorbed at the top and bottom of the enclosure, convective instability unfolds at the bottom of 

the enclosure owing to an adverse temperature gradient stemming from the unstable density gradient driven 

by buoyancy, while the top of the enclosure is convectively stable since the density gradient is stably 

stratified. These results are in agreement with the findings of Townsend (Ref. 4) and Elder (Ref. 8) that 

convective motion at early times is dependent on the thermal boundary layer thickness and is independent of 

the height of the enclosure. As shown in Figure 3, the six cells that are generated initially are segregated (at t

= 300 s) toward the bottom wall corner, similar to the findings of Elder (Ref. 7). The growth of these cells 

yields small-scale thermals (at t = 432 s) (as shown by the temperature field) while the six cells self-organize 

into two larger cells with small secondary cells (as shown by the flow field). The amalgamation of the cells 

indicates a local increase in the wavelengths of the cells as also observed by Elder (Ref. 7). According to 

Elder (Ref. 8), the amalgamation of cells was also observed in laboratory experiments and represents a 

broadening of the spatial spectrum of the velocity field to lower wavenumbers or higher wavelengths.

The transition from short-time-scale events, in which thermals break off, to long-time-scale events, 

which mark the propagation of the instability, occurs for t > 432 s. The shear instability of the top 

convective cell, with the quiescent fluid on top, forms a local source that excites the growth of secondary 

vortices (t = 432 s). Self-repetition of the internal shear instability mechanism leads to the propagation of the 

local convective instability front toward the top of the enclosure (t = 1728 s). As the top of the enclosure is 

approached, this system of four symmetric counterrotating vortices becomes asymmetric, and the breakdown of symmetry in the system is shown for t = 1872 s. Regions of shear instability between the four cells 

bifurcate to a series of secondary cells, which leads to self-organization of the flow field into a set of four 

asymmetric cells. These cells are precursors to a traveling-wave-type motion of convective modes. The 

traveling-wave-type mode occurs after the computation time of 1872 s. However, as will be shown, 

increasing the background heat load Q to 0.01 W or Ra to 1010 increases the speed at which the instability 

propagates toward the top of the enclosure. Note that in Figure 3, and in the following figures, the local 

coordinates, xp and yp, are used for graphics which are normalized, for convenience, with respect to L and H, 

respectively. However, the normal dimensionless coordinates, x

* and y*

, which bound the internal region of 

the enclosure, are normalized with respect to the horizontal length scale L. 

Figure 2.—Base flow field, due to sidewall heat flux, which becomes unstable owing to local 

convective instability at the bottom boundary, where ψ is the stream function. Heat flux ratios at the 

bottom, top, and sidewall, fb = 1, ft = 1, and fs = 2.87×10–5

, respectively; aspect ratio, Ar = 3.24; 

Prandtl number, Pr = 2.27; Raleigh number, Ra = 109.
NASA/TM—2014-217442 9

Figure 3.—Global dynamics of the flow field, showing the growth of thermals over time, t, and the 

propagation of the flow field leading to incipient asymmetry; heat fluxes at the sidewall, bottom, 

and top, ′′ qs = 0.001 erg/cm2

-s and q q b t ′′ ′′ = = 3.49×101 erg/cm2-s, respectively; heat flux ratios at 

the bottom, top, and sidewalls, fb = 1, ft = 1, and fs = 2.87×10–5

, respectively; heat load, Q = 0.001 

W; aspect ratio, Ar = 3.24 (height and width of enclosure, H = 61.9 cm and L = 19.1 cm); Prandtl 

number, Pr = 2.27; and Raleigh number, Ra = 109

. Results are shown for increasing time.
NASA/TM—2014-217442 10

4.2 Traveling-Wave-Type Mode

The effect of increasing the absorption of heat on the boundary of the enclosure from the ambient 

environment is to increase the speed of transition to a traveling-wave-type motion of convective modes. 

For the same time scale as in Figure 3 (t = 1872 s), Figure 4 shows for Ra = 1010 the propagation of the 

basic four symmetric modes (t = 432 s) and their breakdown to asymmetric modes (t = 576 s), which 

leads to the traveling-wave-type motion of convective modes (t = 864 to 1872 s). Unlike for Ra = 109

, the 

transition to asymmetric modes (t = 576 s) occurs before the instability front reaches the top of the 

enclosure. The asymmetric front t = 576 s (t < 792 s) propagates toward the top of the enclosure through 

shear instability, which generates secondary modes until the top of the enclosure is reached (t = 792 s) 

and a fourth mode emerges. The fourth mode is a precursor to a traveling-wave-type mode (t = 864 s) that

propagates between the top and bottom of the enclosure as time increases. As shown by the vorticity field 

in Figure 4, the dynamics of this traveling-wave-type mode generates internal vorticity that mixes the 

liquid. Consequently, the temperature is uniform along the height of the enclosure even though, as will be 

shown, the flow and temperature fields are oscillating in time locally (t = 1872 s).

The characteristic observation time for the laboratory experiments for the case of no heater and no 

screen (Ref. 1) is t = 600 s. Figure 5 shows the characteristics of the flow and vorticity fields for a time 

scale of 720 s in order to gauge typical convective events. Dynamical events in Figure 5 are shown as a 

function of increasing heat absorption Q from 0.001 to 1.0 W, corresponding to increasing Ra from 109 to 

1012. The effect of increasing the background heat load Q is to increase the speed at which the instability 

propagates toward the top of the enclosure and its transition to a traveling-wave-type mode. Figure 5 

shows that the shear instability between the asymmetric modes generates the production of internal 

vorticity in which the intensity increases with increasing Rayleigh number. As Ra increases from 1010 to 

1012, the slight decrease in the maximum temperature near the top of the enclosure is evidence of the 

effectiveness of mixing due to the intense production of vorticity. On the other hand, there is a fast 

increase in the maximum velocity—from 0.83 to 7.48 cm/s. Even though the background heat load from 

the experiment (Ref. 1) is unknown, Figure 5 gives an indication of the effect of a range of heat-load

conditions on convective instability. For Ra = 109 and Q = 0.001 W, the instability penetrates only about 

20 percent of the height of the enclosure over the time scale of the experiment. Thus, it is unlikely that a 

heat load of Q = 0.001 W would be absorbed. The temperature measurements from the experiments

indicate that convective instability has penetrated the entire height of the enclosure, implying that

Q ≥ 0.1 W. Even though the depth of penetration of the instability helps to establish a lower bound on the 

heat absorption, it is the magnitude of the increase in temperature with time that dictates the neighborhood of the heat-load absorption. This will be shown by probing the flow field at various fixed points 

(xp, yp).

4.3 Small-Scale Flows

As the background heat load increases beyond Q = 0.01 W (Ra = 1010), it becomes challenging to 

obtain computational solutions owing to the increase in mode numbers, which characterizes the onset of 

convective instability. Figure 6 shows the bifurcation to higher mode numbers for increasing heat load or

Rayleigh number. The first bifurcation—from six to eight modes—occurs for Ra = 1011, which is characterized by intensive vorticity production near the bottom boundary. The second bifurcation—from 8 to 10

modes—occurs for a decade increase in Rayleigh number, Ra = 1012, or heat load, Q = 1 W. After another

decade increase in Rayleigh number, Ra = 1013, the mode numbers double to approximately 20 modes. 

This explosive increase in mode numbers required an increase in grid size to 1500 by 1500 for small-scale 

resolution. The explosive increase in mode numbers caused a considerable challenge in obtaining computational solutions of time-asymptotic dynamical states for Ra = 1013. Small-scale thermals erupt from the 

bottom wall on the same time scale as the growth of the instability; this is shown for Ra = 1013 and 1011 in

Figure 6.
NASA/TM—2014-217442 11

Figure 4.—Transition to asymmetry leading to the traveling-wave-type motion of convective modes; 

heat fluxes at the sidewall, bottom, and top, ′′ qs = 0.001 erg/cm2

-s and q q b t ′′ ′′ = = 3.49×102 erg/cm2-s, 

respectively; heat flux ratios at the bottom, top, and sidewalls, fb = 1, ft = 1, fs = 2.87×10–6

, 

respectively; heat load, Q = 0.01 W; aspect ratio, Ar = 3.24 (height and width of enclosure, 

H = 61.9 cm and L = 19.1 cm); Prandtl number, Pr = 2.27; and Rayleigh number, Ra = 1010. 

Results are shown for increasing times, t.
NASA/TM—2014-217442 12

Figure 5.—Effect of heat flux absorption at the sidewall, bottom, and top ( ′′ qs = 0.001 erg/cm2-s and 

′′ ′′ ′′ = = , qqq b t respectively) on vorticity production over a time scale of t = 720 s. Heat flux ratios 

fb = 1, ft = 1, fs = 2.87×10–5 to 2.87×10–8 (corresponding to heat loads of Q = 0.001 to 1.0 W); 

aspect ratio, Ar = 3.24 (height and width of enclosure, H = 61.9 cm and L = 19.1 cm); and Prandtl 

number, Pr = 2.27. Results are shown for increasing heat flux absorptions, q"; heat loads, Q; and 

Rayleigh numbers, Ra.
NASA/TM—2014-217442 13

Figure 6.—Bifurcation to higher mode numbers for increasing Rayleigh numbers (grid size increases from 

1000 by 1000 to 1500 by 1500 for Rayleigh number, Ra = 1013) and small-scale thermals. Aspect ratio, 

Ar = 3.24; and Prandtl number, Pr = 2.27. Results are shown for different times, t; heat fluxes, q"; heat 

loads, Q; heat flux ratios at sidewalls, fs; and Rayleigh numbers, Ra.

Since the resolution of small scales becomes important for the highest Rayleigh number, Ra = 1013, or 

heat load, Q = 10 W, considered, we show in Figure 7 the effect of grid size, which illustrates sensitivity to 

initial conditions. At Ra = 1013, the system is highly nonlinear and, as will be shown, becomes chaotic for 

Ra < 1013. A chaotic dynamical system implies sensitivity to initial conditions. For a time scale of t = 9 s,

the increase in grid numbers from 1000 by 1000 to 2500 by 2500 serves as a test for sensitivity to initial 

conditions. The increase in grid numbers represents a perturbation to the base flow solution and shows slight 

nuances in the dynamics of the modes; however, a grid size of 2000 by 2000 is adequate to resolve the 

short-time-scale dynamics of the growth of modes. The difference between 2000 by 2000 and 2500 by 2500 

is due to the fact that the dynamics of mode interactions is also sensitive to the grid number density. In 

addition, self-organization of the modes occurs rapidly, which can result in slight nuances in mode numbers.

This implies that a pattern similar to that for the 2000 by 2000 grid occurs for a 2500 by 2500 grid, but at

t < 9 s. Because of this sensitivity to the initial condition, the exact patterns are not be repeated at fixed time 

t = 9 s as the grid size increases. However, certain basic patterns remain the same: that is, the explosion to 

higher mode numbers and the increase in vorticity production near the bottom wall.

Figure 8 shows the self-similarity of the propagation of the convective instability for a 2000 by 2000

grid at Ra = 1013, which corresponds to the upper bound in heat load. This figure illustrates that, beyond 

the onset of convective instability at t = 8 s, in which there is an explosion in mode numbers, the growth 

of the instability, 9.6 s ≤ t ≤ 19.2 s, occurs through the self-organization of modes, which results in the 

basic four symmetric counterrotating vortices that also occur for the lowest Rayleigh number (Ra = 109

).

Beyond t > 19.2 s, this convective state follows trends similar to those shown for the lower Rayleigh 

numbers: that is, transition to asymmetry leading to a traveling-wave-type mode. As the flow propagates 

upwards, there is an increase in the small-scale internal production of vorticity in the bulk flow field.
NASA/TM—2014-217442 14

Figure 7.—Illustration of sensitivity to initial conditions showing the effect of grid size on the resolution of 

small scales at time, t = 9.0 s, for Rayleigh number, Ra = 1013; heat flux, q" = 3.49×105 erg/cm2

-s; heat 

load, Q = 10.0 W; heat flux ratios at the bottom, top, and sidewall, fb = 1, ft = 1, fs = 2.87×10–9

, 

respectively; aspect ratio, Ar = 3.24; and Prandtl number, Pr = 2.27.

Figure 8.—Evolution of instability at small scales, showing the self-organization of modes and vorticity 

production as the instability propagates from time, t, of 8 s to 19.2 s. Heat flux, q" = 3.49×105 erg/cm2

-s; 

heat load, Q = 10.0 W; aspect ratio, Ar = 3.24; Prandtl number, Pr = 2.27; Rayleigh number, Ra = 1013; 

and grid size, 2000 by 2000.
NASA/TM—2014-217442 15

4.4 Scalar Measures

For the time scale of the experiments (Ref. 1), Figures 9(a) and (b) show the characteristic maximum 

velocity and temperature due to heat absorption from the ambient background environment. There is a sharp 

increase in the maximum velocity as the heat absorption or the Rayleigh number increases (Fig. 9(a)). The 

fast increase in velocity is the cause for the decrease in time scale for the propagation of the convective 

instability along the height of the enclosure. An increase in the heat absorption from the environment causes 

a rise in the maximum temperature at the top boundary (Fig. 9(b)). For a short time (t = 360 s), the maximum temperature remains nearly constant, and for a long time (t = 720 s), the maximum temperature 

decreases as the Rayleigh number increases due to the fact that intensive convective motion mixes the fluid 

in the top boundary layer near the wall.

Figure 9.—Maximum velocity (Vmax) and temperature (Tmax) as a 

function of Rayleigh number, Ra. Aspect ratio, Ar = 3.24; Prandtl 

number, Pr = 2.27; heat flux ratios at the bottom and top, fb = 1 and 

ft = 1, respectively. (a) Maximum velocity in flow field as heat 

absorption from the ambient environment increases for a fixed time. 

(b) Effect of heat absorption from the ambient environment on 

maximum temperature inside the enclosure for a fixed time.
NASA/TM—2014-217442 16

Figure 10.—Temperature profile along the vertical axis at xp = 2.87 cm and 

fixed time t = 648 s for increasing Rayleigh numbers at ambient conditions

for two heat loads: heat flux ratios at the sidewall, fs = 2.87×10–5 (heat flux, 

q" = 3.48×101 erg/cm2

-s; heat load, Q = 0.001 W) and fs = 2.87×10–8

(q" = 3.49×104 erg/cm2

-s, Q = 1.0 W). Aspect ratio, Ar = 3.24; Prandtl 

number, Pr = 2.27; and heat flux ratio at the bottom and top, fb = 1 and 

ft =1, respectively.

The temperature profile along the vertical axis is contrasted for low Q = 0.001 W (Ra = 109

) to moderate Q = 1 W (Ra = 1012) heat absorption in Figure 10. Over the time scale of the experiment (Ref. 1) 

(t = 600 s), the propagation of convective instability penetrates about 20 percent of the height of the enclosure for Q = 0.001 W, whereas the instability propagates the entire height of the enclosure for Q = 1 W. For 

Ra = 1012, the temperature remains nearly constant throughout the bulk of the liquid, with a large temperature gradient in the stably stratified layer near the top wall and a well mixed layer near the bottom wall.

4.5 Local Dynamics at Fixed Points 

The response of the temperature and flow fields at a fixed point ( 0.15, 0.1) * * x p = y p = near the bottom 

of the enclosure for long (Ra = 109 to 1012) and short (Ra = 1013) time bands for various heat loads from the 

ambient laboratory environment is shown in Figure 11. These predictions indicate the response time for the 

onset of convective instability as the Rayleigh number or heat absorption increases. The response time of 

the instability at 10 percent of the enclosure height decreases from t = 400 s for Ra = 109 to t = 15 s for 

Ra = 1013 as more heat is absorbed from the environment. Both the temperature and velocity fields oscillate 

aperiodically from Ra = 109 to 1013, this is due in part to the large dimensions of the enclosure, Ar = 3.24 

(H = 61.9 cm and L = 19.1 cm). The velocity and temperature fluctuations increase as the heat load or 

Rayleigh number increase. The time series for the velocity field component is nearly stationary in the mean, 

varies about a fixed level for Ra = 109

, and varies slightly about the mean for a segment of the time series

for Ra = 1010 and 1011. However, for Ra = 1012 to 1013, the velocity field is nonstationary in the mean, does

not vary about a mean fixed level, and has an upward trend. This implies that the mean varies with time for 

high Rayleigh numbers. An estimation of an average value from such a time series should take into account

the variation of the mean with time. On the other hand, the temperature field is stationary for Ra = 109 to 

1010, and a nonstationary time series occurs for Ra = 1011 and 1012 for t > 200 s, which started out as 

stationary for the short time band t < 200 s. Since we show only a short time band of the time series for 

Ra = 1013, stationary behavior appears; however, a nonstationary trend would ensue for a time band on the 

order of 800 s, similar to Ra = 1012. These trends indicate that the averaging of temperatures for various 

points on a plane should take into account the nonstationary behavior of the time series so that the

temperature histories can be correlated properly with the convective states of the flow field.
NASA/TM—2014-217442 17

Figure 11.—Time history of nondimensional temperature, T*, and velocity fields (horizontal component of 

velocity, u) near the bottom of the enclosure ( . , . ) * * xp = 0 15 y p = 0 1 , showing an increase in intensity of the 

fluctuations as the Rayleigh number, Ra (heat absorption) increases for aspect ratio, Ar = 3.24 (height and 

length of enclosure, H = 69.1 cm and L = 19.1 cm); Prandtl number, Pr = 2.27; and heat flux ratios at 

bottom and top of enclosure, fb = 1 and ft = 1, respectively. Results are shown for different heat flux ratios 

at the side wall, fs, heat fluxes, q", heat loads, Q, and Rayleigh numbers, Ra.
NASA/TM—2014-217442 18

Since the temperature at fixed points (xp, yp) was measured experimentally along the height of the

enclosure (Ref. 1), Figure 12 provides the time history of the temperature field for various points along 

the enclosure height, contrasting low (Q = 0.01 W) and moderate (Q = 1.0 W) heat absorption. The 

experimental temperature measurements (Ref. 1) reported were averaged values over a plane consisting of 

five points at each * y p location—near 0.1, 0.35, 0.65, and 0.90—over a time span of t = 600 s. Figure 12

shows boundary layer locations at * y p = 0.001 and 0.999 near the bottom and top wall boundaries to 

illustrate the instantaneous response of the fluid so that it can be compared with experimental trends. For 

Ra = 1010, in response to heat absorption at the bottom boundary ( 0.001), * y p = there is a sharp rise in 

temperature over a time scale of t = 100 s, and the onset of oscillations follows for t > 100 s. As the 

convective instability propagates toward the top of the enclosure, there is a delay time τD for the onset of 

convective instability: for * y p = 0.1, τD = 150 s; whereas for * y p = 0.9, τD = 750 s. 

Similar to results for the bottom boundary, heat absorption at the top boundary increases the local 

temperature; however, owing to the stably stratified layer on top ( 0.999), * y p = the temperature continues 

to rise steadily past t = 100 s up to t = 800 s, beyond which the onset of convective instability begins. 

Based on the delay time of τD = 800 s for convective instability to propagate the entire length of the 

enclosure, a heat load of Q = 0.01 W is not sufficient. For Ra = 1010 (Q = 0.01 W), the delay time of 

τD = 750 s at * y p = 0.9 is too long for the propagation of convective instability in comparison to the 

experimental (Ref. 1) measurement time band of t = 600 s. This implies that the instability would not 

have reached the top of the enclosure. Thus, a background heat load of Q = 0.01 W is too small. However, 

a decade increase to Ra = 1011 (Q = 0.1 W) indicates a delay time of τD = 400 s at * y p = 0.9, which falls 

within the experimentally measured time band of t = 600 s. Thus, Ra = 1011 (Q = 0.1 W) can potentially 

serve as a lower bound for heat absorption from the ambient laboratory environment. However, when 

temperature magnitudes are compared with experimental values for Ra = 1011, they fall below the 

experimental measurements (Ref. 1). Another decade increase in heat load to Q = 1.0 W, Ra = 1012, 

shows a gradual approach to reported experimental values of average temperature and serves as a lower 

bound for heat absorption from the ambient environment. 

The dynamics of the temperature field at fixed points (xp, yp) spanning the height of the enclosure in 

Figure 12, shows a decrease in the penetrative convection response time as the Rayleigh number 

increases. In comparison to Ra = 1010, the delay response time of penetrative convective instability is 

reduced for Ra = 1012 as the height increases. The response time ranges from τD = 25 to 155 s, corresponding to * y p = 0.1 and 0.9, respectively. As the probing distance ( ) * y p increases from the bottom, the 

signal varies from a stationary time series about a constant mean value ( 0.001) * y p = to a nonstationary 

time series about a mean value that increases with time ( 0.1). * y p = This trend is maintained as the 

enclosure height is traversed to the top boundary ( 0.9), * y p = which indicates intermittency. In addition, 

the convective instability penetrates the top boundary layer at * y p = 0.999, showing intermittency through

periodic bursts of intensive convective motion. Owing to the lag in response time τD associated with 

locations from * y p = 0.1 to 0.9 (τD = 25 to 155 s), a comparison with experiments over the short time 

scale of t = 60 s is inadequate for 0.3 * y p > since the experimental measurements indicate nearly instantaneous rise in temperature as shown at * y p = 0.001 and 0.999. However, for a time band of t ≥ 200 s,

adequate comparisons can be made with experimental data. 

The computational model accounts for the response of the temperature field to a step input of heat 

absorption from the boundaries and represents a clearly posed initial value problem. However, there exists 

uncertainty associated with the initial condition of an experimental system. Since heat is absorbed 

instantaneously, the origin of time (t = 0) may be uncertain, considering the complexity of temperature 

measurements and the establishment of the initial saturation condition for liquid nitrogen. Considering the 
NASA/TM—2014-217442 19

Figure 12.—Comparison of low (heat load, Q = 0.01 W, and heat flux ratio at the sidewall, fs = 2.87×10–6) and 

moderate (Q = 1.0 W and fs = 2.87×10–8

) heat absorption from the ambient environment on penetrative 

convective instability along the height of the enclosure, as well as the typical dynamical state for Rayleigh 

number, Ra = 1012; aspect ratio, Ar = 3.24 (height and length of enclosure, H = 69.1 cm and L =19.1 cm); 

Prandtl number, Pr = 2.27; and heat flux ratios at bottom and top of enclosure, fb = 1 and ft = 1, 

respectively. T, temperature; t, time; q", heat flux; k, constant. Results are shown at different points inside 

the enclosure ( , ). * * xp y p
NASA/TM—2014-217442 20

complexity of an experiment, we conjecture that there could be uncertainty associated with the delay time

τD for penetrative convection that cannot be captured in an experimental system. Nevertheless, as is 

shown in Section 4.6, the temperature magnitude over the time scale of the experiment, t = 600 s,

indicates that the absorption of a heat load of Q = 1.0 W from the background environment can serve as a 

lower bound for a potential scenario of the actual experiment. These results imply that the temperature 

signals from an experimental system, prior to averaging, are likely to be nonstationary, with a mean value 

that varies with time; for the experimental system (Refs. 1 and 2) considered, only averaged values of 

temperature were reported. A nonstationary time series requires proper treatment of the mean value

(Ref. 12). Without proper averaging of a nonstationary time series, uncertainty can be introduced in the 

average value of temperatures. 

The time series shown for Ra = 1012 at 0.9 * y p = is quantified using the amplitude spectrum of the 

u component of the velocity field, defined as

2

0

1 ( ) ()

TR

if t

R

U f ute T

− π = ∫  (24)

In relation to Figure 12 showing the dynamical state for Ra = 1012, Equation (24) shows the square root of

the local kinetic energy per unit time, or the square root of power (Pu(f)) along the ordinate axis; the

abscissa represents the spectrum of the frequency of the time signal; and TR denotes the interval of the 

time series. A convolution relationship, U(f ) = WH(f ) * U (f ), using a Hanning spectral window WH(f ) is 

used to obtain a smooth estimate U(f ). The frequency decomposition of U(f ) in Figure 12 indicates a 

long-period event that has a low-frequency peak of 0.005 Hz and spans a broadband of 0.5 Hz. The 

broadband response falls below the critical Nyquist frequency of fc = 3.5 Hz, indicating that the effective 

bandwidth of the signal is free from aliasing. Note that the power spectrum, Pu(f ) = U(f ) × U*

(f )/2π, is

the product of U(f ) with its conjugate U*

(f ), which in terms of Equation (24) is the square of the 

magnitude of the Fourier-transformed velocity component u(t). The power spectrum gives results that are 

similar to those for the amplitude spectrum. However, the amplitude spectrum was found to have the 

advantage of showing the higher harmonics of the fundamental frequency on a linear scale. A

pseudophase space trajectory, obtained by lagging the time signal of u(t + k) by a constant value k, shows 

a dense trajectory. The combination of a broadband power spectrum and a dense pseudophase space

trajectory indicates a chaotic state. We conjecture that this chaotic dynamical state exists for the entire 

range of Rayleigh numbers considered—even for the lowest heat absorption (Q = 0.001 W) at Ra = 109

—

owing to the large size of the enclosure.

Related works on liquid helium by Behringer (Ref. 13) aimed at characterizing the onset of lowdimensional turbulence or chaos show that the trend of the critical Rayleigh number as a function of 

aspect ratio merges with the prediction of Charlson and Sani (Refs. 14 and 15) for long vertical walls, 

which shows an increase of the critical Rayleigh number for the onset of instability. In contrast,

enclosures with large horizontal extent show that turbulence occurs near the critical Rayleigh number. 

These observations parallel our conjecture of chaos at the lowest Rayleigh number since our enclosure 

has relatively long vertical and horizontal length scales.

4.6 Comparison of Model With Experiment

Figure 13 compares the mean fluctuating temperature T at fixed points (xp, yp) along the height of the 

enclosure with experimental data (Experiments (Exp.)—no screen, no heater) after Bolshinskiy et al.

(Ref. 1). The model uses the boundary condition of heat load from Q = 1.0 to 10 W (Ra = 1012 to 1013

) 

due to heat absorption from the ambient laboratory environment. The experimental data (o) correspond to 

average values of temperatures <T> for five points over a plane at each height (y) location. We compare 

the experimental data with the mean temperature T at a fixed location in order to show unfiltered 
NASA/TM—2014-217442 21

predicted temperatures without plane averaging; note the unfiltered temperatures in Figures 11 and 12 

reveal inherent fluctuations that may not be captured by averaging. The experimental data points are in 

close proximity to the computational model, and for comparison, the figures show data for a long time,

t = 600 s at Ra = 1012, and a short time, t = 50 s for Ra = 1013. Note that the delay time τD = 15 s for

Ra = 1013 is very short for high heat flux absorption and falls within the response time of the experimental 

data. For Ra = 1012, a heat load of Q = 1 W yields temperature magnitudes that are lower than those

measured in the bulk liquid; however, in the boundary layer near the top wall, the mean temperature,

T K = 79.5 , approaches the bulk experimental temperature, <T> = 80.8 K, of the liquid in the middle 

region. 

The model predicts the trend shown by the data: that is, a nearly uniform temperature in the central 

region of the enclosure, a nearly adverse temperature gradient near the bottom of the enclosure due to 

local convective instability, and a stable temperature gradient near the top region due to stable stratification. For an ambient heat load of Q = 10 W (Ra = 1013) in the neighborhood of the maximum value 

estimated from thermal radiation due to the background laboratory environment (T∞ = 293 K), the 

saturation temperature of liquid nitrogen, TA = 77 K, nearly matches the experimental data for short times. 

In this case, the temperature in the boundary layer near the bottom wall, T = 77.5K, nearly matches the 

highest temperature measured experimentally near the top region, <T> = 77.6 K. In contrast, the temperature in the boundary layer near the top wall, T = 83.5K , exceeds the highest measured temperature in the 

top region, <T> = 82 K, for the long duration time of t = 600 s. Given the potential uncertainty associated 

with establishing the proper initial condition that can occur for experimental systems, a minimum-tomaximum heat load range due to a background environment of Q = 1 to 10 W (Ra = 1012 to 1013) can be 

representative of actual experimental conditions. 

Figure 13.—Comparison of mean fluctuating temperature T at fixed points (xp, yp) in the flow field with experimental 

data (Exp.—no screen and no heater) (Ref. 1) for increasing Rayleigh numbers, Ra. (a) Time, t = 600 s. (b) Time, 

t = 50 s.
NASA/TM—2014-217442 22

5.0 Summary and Conclusion

A model problem is addressed on the effect of heat absorption from the background laboratory 

environment on the temperature rise of liquid nitrogen at saturation conditions (ρA, TA) inside a large 

enclosure (H = 61.9 cm and L = 19.1 cm). The model contributes toward an understanding of the flow 

field dynamics for an experimental system aimed at investigating heat entrapment effects in liquid 

acquisition devices. For estimated background heat load conditions, we consider a range from Q = 0.001 

to 10 W, corresponding to Rayleigh numbers from Ra = 109 to 1013. The estimated heat flux of the top 

and bottom boundaries (q q t b ′′ ′′ = ) ranges from 3.49×101 to 3.49×105 erg/cm2

-s, and a conservative 

estimate of qs

′′ = 0.001 erg/cm2

-s is used for the sidewalls owing to the use of insulation at the side 

boundary in the experiments. In response to heat load, the flow field dynamics for Ra = 109 to 1013 show 

self-similarity. The flow field is driven by local convective instability, against a free convection base 

flow, generated at the bottom of the enclosure. The flow field shows small-scale thermals for the short 

time scale and a traveling-wave-type motion of convective modes for the long time scale.

Finite heat leakage at the sidewalls generates the initial free-convection base flow, which consists of 

two counterrotating cells that become unstable against convective instability caused by an adverse 

temperature gradient at the bottom boundary, while the top boundary forms a stably stratified configuration. A number of modes arise, depending on the Rayleigh number, which self-organize into two modes 

owing to amalgamation prior to the propagation of the convective instability along the height of the 

enclosure. A number of symmetric vertical modes arise; these modes become asymmetric prior to 

reaching the top of the enclosure depending on the Rayleigh number. A traveling-wave-type motion of 

convective modes follows that generates aperiodic oscillations of the flow and temperature fields and is

representative of the asymptotic dynamics.

The time history of the temperature and velocity fields shows a characteristic delay time τD for the 

penetration of the convective instability along the height of the enclosure due to heat load absorption. The 

estimated heat loads, which serve as lower and upper bounds, owing to the background environment 

range from Q = 1 to 10 W, and show agreement with experimental trends (Ra = 1012 to 1013). For heat 

loads from Q = 1 to 10 W, the time signal is nonstationary over a time band and spatial location range,

with a mean that increases with time. Thus, cautionary measures should be taken for averaging of 

nonstationary temperature measurements. A chaotic dynamical state (stochastic time signal) seems to be 

the norm for the range of Rayleigh numbers considered, Ra = 109 to 1013, owing to the large size of the 

enclosure.

The computational model sheds light on temperature measurements of saturated liquid nitrogen inside 

an enclosure subjected to an unknown heat load from a background laboratory environment for application to liquid acquisition devices. To understand the effect of the applied heat flux at the bottom of the 

enclosure on the temperature rise of liquid nitrogen, it is necessary to establish the base condition of the 

experiment. This base condition, which was characterized experimentally by measurements of the 

temperature rise due to the background environment (Exp.—no screen, no heater) can be used to establish 

the boundary conditions. The computational model shows that the boundary conditions for the experiment 

lie between Q = 1 and 10 W (Ra = 1012 and 1013). These established boundary conditions can be used as 

input for the top heat flux qt

′′ along with the estimated side heat flux qs

′′ in order to analyze the next phase 

of the experiment, that is, the effect of independent variation of the bottom heat flux qb

′′ . This implies

that a bottom heat-load variation from a minimum of Q = 10 W to the maximum experimental input of 

Q = 104 W can be investigated. 
NASA/TM—2014-217442 23

Symbols

Ar aspect ratio

f heat flux ratio

fc critical Nyquist frequency

go gravitational acceleration

gy body force in y direction

H height of enclosure

k time lag constant

L horizontal length of enclosure

n coefficient representing gravitational level, where n = 1 represents the gravity on Earth

Pu(f ) power spectrum of u velocity component

Pr Prandtl number

p pressure

Q heat load from ambient

q'' imposed heat flux

Ra Rayleigh number

T temperature 

T mean fluctuating temperature

TA temperature of saturated liquid nitrogen

TR interval of the time series

T∞ ambient temperature of laboratory

t time

t

o characteristic time for formation of thermals

U( f ) smooth estimate of the amplitude spectrum of the u component of the velocity field (after 

Eq. (24))

U f ( )  amplitude spectrum of the u component of the velocity field (Eq. (24))

u horizontal component of velocity

V velocity magnitude

V

 velocity field 

v vertical component of velocity

WH(f ) Hanning spectral window

xp fixed horizontal point

yp fixed vertical point

α thermal diffusivity

β thermal expansion coefficient

Γ boundary of enclosure

δ thermal boundary layer thickness

κ thermal conductivity

ν kinematic viscosity

ξ single vorticity component in the z direction (Eq. (11))
NASA/TM—2014-217442 24

ρ density

ρA density of the saturated liquid nitrogen

τD delay time

ψ stream function

∀ for all

∇ gradient operator

Λ parametric set (Eq. (22))

Λ' reduced parametric set (Eq. (23))

 magnitude

Subscripts:

b bottom

p pressure (after Eq. (3))

s sidewall

t top

Superscripts:

* dimensionless quantity
NASA/TM—2014-217442 25

References

1. Bolshinskiy, L.G., et al.: Capillary Liquid Acquisition Device Heat Entrapment. NASA/TM—2007-

215074, 2007. http://ntrs.nasa.gov/

2. Bolshinskiy, L.G.; Hastings, L.J.; and Statham, G.: Cryogenic Capillary Screen Heat Entrapment. 

Cryogenics, vol. 48, no. 5, 2008, pp. 232–237.

3. Duval, W.M.B.; Chato, D.J.; and Doherty, M.P.: Heat Entrapment Effects Within Liquid Acquisition 

Devices. AIAA 2010–1298, 2010.

4. Townsend, A.A.: Temperature Fluctuations Over a Heated Horizontal Surface. J. Fluid Mech., 

vol. 5, no. 2, 1959, pp. 209–241.

5. Howard, L.N.: Convection at High Rayleigh Number. Applied Mechanics, Henry Gortler, ed.,

Springer-Verlag, Berlin, 1966, pp. 1109–1115.

6. Elder, J.W.: Steady Free Convection in a Porous Medium Heated From Below. J. Fluid Mech., 

vol. 27, no. 1, 1967, pp. 29–48. 

7. Elder, J.W.: Transient Convection in a Porous Medium. J. Fluid Mech., vol. 27, no. 3, 1967,

pp. 609–623.

8. Elder, J.W.: The Unstable Thermal Interface. J. Fluid Mech., vol. 32, no. 1, 1968, pp. 69–96.

9. Elder, J.W.: The Temporal Development of a Model of High Rayleigh Number Convection. J. Fluid 

Mech., vol. 35, no. 3, 1969, pp. 417–437.

10. Deardorff, J.W.; and Willis, G.E.: The Effect of Two-Dimensionality on the Suppression of Thermal 

Turbulence. J. Fluid Mech., vol. 23, no. 2, 1965, pp. 337–353.

11. Sparrow, E.M.; Husar, R.B.; and Goldstein, R.J.: Observation and Other Characteristics of Thermals. 

J. Fluid Mech., vol. 41, no. 4, 1970, pp. 793–800.

12. Wei, William W.S.: Time Series Analysis: Univariate and Multivariate Methods. Addison-Wesley

Publishing, 1990, pp. 67–84.

13. Behringer, R.P.: Rayleigh-Bénard Convection and Turbulence in Liquid Helium. Rev. Mod. Phys., 

vol. 57, no. 3, Part 1, 1985, pp. 657–687.

14. Charlson, G.S.; and Sani, R.L.: Thermoconvective Instability in a Bounded Cylindrical Fluid Layer. 

Int. J. Heat Mass Transfer, vol. 13, no. 9, 1970, pp. 1479–1496.

15. Charlson, G.S.; and Sani, R.L.: On Thermoconvective Instability in a Bounded Cylindrical Fluid 

Layer. Int. J. Heat Mass Transfer, vol. 14, no. 12, 1971, pp. 2157–2160.
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-06-2014

2. REPORT TYPE

Technical Memorandum

3. DATES COVERED (From - To)

4. TITLE AND SUBTITLE

Transient Convection Due to Imposed Heat Flux: Application to Liquid-Acquisition

Devices

5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S)

Duval, Walter, M.B.; Chato, David, J.; Doherty, Michael, 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-17597

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-2014-217442

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified-Unlimited

Subject Categories: 31 and 34

Available electronically at http://www.sti.nasa.gov

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

13. SUPPLEMENTARY NOTES

14. ABSTRACT

A model problem is considered that addresses the effect of heat load from an ambient laboratory environment on the temperature rise of

liquid nitrogen inside an enclosure. This model has applications to liquid acquisition devices inside the cryogenic storage tanks used to

transport vapor-free propellant to the main engine. We show that heat loads from Q = 0.001 to 10 W, with corresponding Rayleigh numbers

from Ra = 109

 to 1013, yield a range of unsteady convective states and temperature rise in the liquid. The results show that Q = 1 to 10 W

(Ra = 1012 to 1013) yield temperature distributions along the enclosure height that are similar in trend to experimental measurements.

Unsteady convection, which shows self-similarity in its planforms, is predicted for the range of heat-load conditions. The onset of

convection occurs from a free-convection-dominated base flow that becomes unstable against convective instability generated at the bottom

of the enclosure while the top of the enclosure is convectively stable. For the higher Rayleigh numbers, 1012 to 1013, there is a transition

from a stationary to a nonstationary response time signal of the flow and temperature fields with a mean value that increases with time over

various time bands and regions of the enclosure.

15. SUBJECT TERMS

Cryogenics; Thermal instability; Buoyancy-driven flow; Natural convection

16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF

 ABSTRACT

UU

18. NUMBER

 OF

 PAGES

32

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