- Email: [email protected]
Hydrogen particles in liquid helium
Cryogenics 44 (2004) 507–514 www.elsevier.com/locate/cryogenics
Hydrogen particles in liquid helium J. Xu a, K.M. Smith a, D. Celik b, M.Y. Hussaini a
a,*
, S.W. Van Sciver
b
School of Computational Science and Information Technology, Florida State University, 400 Dirac Science Library, Tallahassee, FL 32306, USA b National High Magnetic Field Laboratory, Florida State University, 1800 E. Paul Dirac Drive, Tallahassee, FL 32310, USA
Abstract Numerical investigations of phase changes during the generation of solid hydrogen particles in cryogenic helium, which include vaporization and re-condensation of helium around a hydrogen droplet and solidification of the droplet in helium, are conducted to unravel the basic mechanisms underlying the phase change phenomena. In a single-fluid, two-phase mixture framework, the simulation of vaporization/condensation involves a vapor–liquid exchange model that accounts for mass and energy transfers due to phase changes between liquid and vapor helium. The enthalpy-porosity, solidification–melting technique is employed to investigate the hydrogen solidification. The volume fraction of each phase in the solidification process is determined using a volume of fluid (VOF) two-phase model, along with a piecewise-linear approach to trace the liquid/solid interface. The results throw some light on the helium flow field resulting from the phase change and the motion of the droplet. They provide some insight into the evolution of liquid–vapor phase-exchange processes as well as the solidification evolution of hydrogen particles of various diameters. 2004 Elsevier Ltd. All rights reserved. Keywords: Solid hydrogen; Phase change; Solidification; Vaporization; Liquid helium; Single-fluid two-phase mixture
1. Introduction Hydrogen is the fuel used by the Space Shuttle Main Engine, as well as by many other hydrogen powered upper stages of rockets (Centaur, and those used by Delta, Atlas, and Ariane rocket upper stages, etc). However, due to its relative low density, larger storage tanks are required, which translates to higher structural weight and less payload. To reduce this unwanted effect, it has been proposed [1,2] that liquid hydrogen be refrigerated to near freezing temperatures and then pressurized to obtain dense, subcooled liquid fuel. Its density can be further increased by solidifying some portion of the liquid hydrogen in the form of small particles. Furthermore, atomic propellants, such as boron, carbon, or hydrogen can be embedded/stored in these small solid hydrogen particles, which can be stabilized and prevented from recombination [1–3], although such embedding techniques are yet to be developed. Thus, atomic hydrogen propellant feed systems may require the production and storage of hydrogen particles in liquid helium.
*
Corresponding author. Tel.: +1-850-6440601; fax: +1-8506451514. E-mail address: [email protected] (M.Y. Hussaini). 0011-2275/$ - see front matter 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.cryogenics.2004.02.017
Solid hydrogen particle formation in liquid helium has been experimentally investigated for years [1–4]. To create solid hydrogen particles in helium, droplets of liquid hydrogen at temperatures around 20 K are injected into a bath filled with liquid helium at 4 K. Since the temperature of the liquid hydrogen is higher than the boiling point of helium (4.2 K at 101,325 Pa [5,6]), and the temperature of ambient helium is lower than the solidification point of hydrogen (13.8 K at 101325 Pa [7]), vaporization, condensation and solidification occur simultaneously during the process. At the beginning, a thin layer of helium vaporizes around the hydrogen droplet, and the hydrogen droplet solidifies in the helium fluid. As the vaporized helium cools to the dew point of helium, it condenses in the helium bath, and reaches the temperature of the liquid helium bath in the end. Experimental efforts have been focused on the visual characterization of the particles, such as observing particle formation in liquid or vapor helium [4], video imaging of particle size [3], timing particle agglomeration [2,3], etc. Detailed experimental investigation of these phenomena is a formidable task due to very small time scales and other experimental design challenges. For example, Palaszewski [3] observed that in the hydrogen freezing process, large clouds of vapor ‘‘. . . would have also obscured the formation process, and
508
J. Xu et al. / Cryogenics 44 (2004) 507–514
thwarted efforts to see the final particles.’’ The authors encountered similar experimental difficulties. Numerical simulations, however, offer a relatively easy and less expensive alternative, and of course, they have to be validated by experimental observations. In addition, such simulations are becoming increasingly attractive due to advances in phase change theory [8–10], twophase models [11–15], and dynamic mesh techniques [16–20]. These numerical methodologies have been implemented in many phase change systems, but there are few numerical studies of the phenomena in helium– hydrogen cases. With the goal of better understanding the process of solid hydrogen particle generation, a series of computational studies have been conducted to investigate helium phase change around a hydrogen particle, and hydrogen droplet freezing in a helium bath. A vapor– liquid exchange model implemented in the framework of a single-fluid, two-phase mixture model is developed to study the vaporization/condensation process. An enthalpy-porosity, solidification–melting model is employed to investigate the solidification phenomena. Dynamic mesh techniques are introduced in the vaporization/condensation simulation to account for the movement of hydrogen droplets after injection. The effects of a phase changing, moving particle on the helium, as well as the evolution of the liquid–vapor helium phase exchange are studied. In addition, the time required for solidification of hydrogen particles with diameters ranging from 0.1 to 2.0 mm are investigated using the model developed.
2. Problem formulation and solution technique Formation of solid hydrogen particles in a helium bath involves a complex energy exchange processes between the two fluids resulting in vaporization of helium around the hydrogen droplet and its consequent condensation, and the simultaneous solidification of the liquid hydrogen droplet. The purpose of the present study is to study these processes in isolation and put them together in a later study, fully expecting that the sum will be greater than the whole as in any complex system. In other words, the ‘‘mushy layer’’ between the solidifying hydrogen droplet and the liquid helium is assumed as an infinitesimally thin interface and the interfacial conditions approximate/model the consequence of the mushy layer physics. Our study is the first cut at this extremely complex problem. 2.1. Vaporization and condensation of helium We assume that a hydrogen droplet of diameter r mm (r 2 ½0:1; 2:0) is injected vertically downwards with an initial velocity of ~ vy at the center of a rectangular bath of
Fig. 1. Simulation configuration. (The hydrogen droplet is shown as a red circular disk in the middle of the domain. The liquid helium is shown as blue, and is contained within the bath of width 8 mm and height 16 mm. Initially the helium temperature is 4 K, and the liquid hydrogen temperature is 20 K.)
dimensions 8 mm in width and 16 mm in depth, as shown in Fig. 1. The two-phase flow, comprising liquid and vapor, is modeled using a single-fluid two-phase approach [11,21]. In the two-phase approach, the phases are assumed to be interpenetrating, meaning that each phase has its own velocity vector field, and within any control volume there is a volume fraction of the primary (liquid) phase, and a volume fraction of the secondary (vapor) phase, which sum to one. 2.1.1. Governing equations The governing equations are the conservation equations for mass, momentum and energy of the mixture. It is expedient to consider these equations in a mesh system that moves with the hydrogen droplet. The continuity equation may be written as oqm þ r ðqm~ vÞ ¼ S_ q : ot
ð1Þ
ug , and ~ ug is the grid velocity of the Here, ~ v ¼~ vm ~ moving mesh, and ~ vm represents the mass-averaged velocity of the mixture, which is defined as ~ vm ¼
al ql~ vl þ av qv~ vv : qm
ð2Þ
The density of the mixture, qm , is defined to be qm ¼ al ql þ av qv ;
ð3Þ
where a represents the volume fraction, and subscripts l and v denote liquid and vapor phase, respectively. Following Engelberg-Forster and Greif [8,10], the mixture source/sink term, S_ q , in Eq. (1) is expressed as T Tb T Tb þ jc a v q v ; S_ q ¼ jv al ql Tb Tb
ð4Þ
J. Xu et al. / Cryogenics 44 (2004) 507–514
where T represents the temperature of ambient helium, Tb is the boiling point of helium, jv stands for the vaporization rate of liquid helium, while jc denotes the condensation rate of vapor helium. The momentum equation for the mixture is oqm~ v þ r ðqm~ vÞ þ qm~ v~ vÞ ¼ rp þ r ðlm r~ g ot ! X ak qk~ vdr;k~ vdr;k ; ð5Þ þ~ F þr k¼l;v
where lm ¼ al ll þ av lv , denotes the viscosity of the mixture; and ~ vdr;k ¼ ~ vkl ~ vm is the drift velocity for phase k, with ~ vkl ¼ ~ vk ~ vl the slip velocity of phase k, defined as the relative velocity of phase k over the primary phase l. Based on the assumption that a local equilibrium between the phases should be reached over spatial length scales, the slip velocity can be given by [21] ~ vkl ¼ skl~ a;
ð6Þ
where skl and ~ a are the relaxation time and acceleration of the secondary phase k, respectively. ~ F in Eq. (5) represents the body force, which is given by ~ F ¼ FD ð~ vm ~ ug Þ qm ;
ð7Þ
with the drag coefficient FD , given by a semi-empirical formula [4], FD ¼
24 ½1 þ 0:15Re0:687 ; p Rep
ð8Þ
where Rep is the particle Reynolds number. The energy equation for the mixture can be written in the form, X d X ðak qk hk Þ þ r ðak~ vk qk hk Þ ¼ r ðkeff rT Þ þ S_ h ; dt k¼l;v k¼l;v ð9Þ where hk is the sensible enthalpy of phase k, keff is the effective conductivity, and S_ h ¼ L Sq , with L the latent heat of helium. The volume fraction of the vapor (secondary) phase is governed by o ðav qv Þ þ r ðav qv~ vÞ ¼ r ðav qv~ vdr;v Þ þ S_ v ; ot where ( b ; if T > Tb ; vaporizes; jv al ql T T Tb S_ v ¼ T Tb jc av qv Tb ; if T < Tb ; condenses:
ð10Þ
ð11Þ
2.1.2. Initial and boundary conditions Initially, the helium bath and the liquid hydrogen droplet are at 4 and 20 K respectively, and the initial velocity of the droplet is 8 mm per second. During the
509
time the droplet freezes to 4 K, there is a complex interplay among the processes of vaporization, condensation and solidification. To study the vaporization/ condensation processes in isolation, we assume that (i) the shape and size of the droplet remains constant, (ii) the temperature evolution at the interface between helium and hydrogen droplet has a linear profile, T ¼ 20:0 c t, where c is the temperature change that the droplet suffers per unit time, and (iii) the droplet moves with y-velocity ~ vy ¼ 8:0 þ d t mm/s starting from the center of the bath. These assumptions are approximately supported by the results of droplet solidification simulation, which further determine c to be 0.4 K/s, and injection simulation, that gives d to be unity. Non-slip boundary conditions are applied at the walls of the bath. 2.1.3. Solution technique For computational purposes (easy management of the dynamic mesh), Eqs. (1), (5), (9) and (10) are transferred into a general integral form written as Z Z d qm /dV þ qm /~ v d~ A dt V Z ZoV ~ Cr/ dA þ S/ dV ; ð12Þ ¼ oV
V
where / represents a general transport quantity. In the continuity equation, / ¼ 1; in the momentum equation, / ¼~ v; and in the energy equation, / ¼ h (enthalpy). C is the diffusion, and S/ represents the source term for the quantity /, ~ A is the area vector, and V represents the control volume, which is identified with each mesh element. The time derivative term in the integral equations is discretized using a first-order backward difference formula, given by Z nþ1 n d ðq/V Þ ðq/V Þ ; ð13Þ q/ dV ¼ dt V Dt where the superscripts n and n þ 1 denote the respective quantity at the current and next time levels. The ðn þ 1Þth control volume V nþ1 is computed as V nþ1 ¼ V n þ
dV Dt; dt
ð14Þ
where dV =dt is the time derivative of the control volume. Following the grid conservation law, dV =dt is computed as Z nf X dV ~ ~ ¼ A¼ Aj ; ð15Þ ug d~ ug;j ~ dt oV j where nf is the number of faces on the control volume, ~ Aj is the j-face area vector and oV denotes the control volume surface. Dynamic mesh methodology [16–20] is employed to update the mesh where the shape of the
510
J. Xu et al. / Cryogenics 44 (2004) 507–514
domain is changing with respect to time due to the droplet motion. A linear spring-analogy smoothing method [19] and a local re-meshing method [18] are used in this simulation. A critical constraint on the moving mesh methodology is that the moving droplet should not advance more than one spatial element in a time step. In other words, if dy is the spatial scale in the direction of motion, dt is the time step, and vy is the moving droplet velocity, then the constraint dy > vy dt should be maintained at each time step. A second-order upwind scheme is used for spatial discretization of the momentum equation and the continuity equations. The multiphase flow equations are solved using the coupled semi-implicit method for pressure-linked equation algorithm [22]. 2.2. Solidification of hydrogen droplets 2.2.1. Governing equations For the problem of a static liquid hydrogen droplet immersed in a helium bath, convection and diffusion are the principal mechanisms of heat transfer. The R T energy equation in terms of the sensible enthalpy, h ¼ Tref cp dT , is oqh þ r ðq~ vhÞ ¼ r ðarhÞ þ Sh ; ot
ð16Þ
where q is the density of liquid or solid hydrogen, ~ v represents the velocity, where ~ v ¼ 0 within the solid and at the interface of the liquid and solid, a is the thermal diffusivity, and Sh is the source term that accounts for the latent heat of phase change. An enthalpy-porosity model proposed by Voller et al. [12–14] and modified by Brent et al. [15], is employed to represent the source term: Sh ¼
oðqDH Þ ; ot
ð17Þ
where DH , the latent heat content, can be written as a function of temperature T , 0; if T < Tm ; DH ¼ ð18Þ L; if T > Tm ; where L is the latent heat (of hydrogen), and Tm is the phase change temperature of hydrogen solidification/ melting. The latent heat is treated as a heat source or sink by assigning a nodal latent heat value to each computational element [13], according to the temperature of the cell. For the solid/liquid two-phase problem, a volume of fluid (VOF) model [11] is implemented to solve the volume fraction of each phase. The essential premise of the VOF model is the assumption that the two phases are not interpenetrating, which is true in the solidification case. Interface tracking between phases is accomplished by solving a continuity equation for the volume
fraction of each phase. The volume fraction is given by ai , where i ¼ s for the solid phase, and i ¼ l for the liquid phase. The volume fraction continuity equations are given by oas þ~ v ras ¼ 0: ot
ð19Þ
A single momentum equation is solved throughout the domain, given by o ðq~ vÞ þ r ðq~ v~ vÞ ¼ r ½lðr~ vÞ þ q~ g þ Sb ; ot
ð20Þ
where the term Sb ¼ q~ gbðT T0 Þ accounts for the buoyancy force of natural convection effects, and b is the thermal expansion coefficient. Volume fractions of phases for the element under consideration and its surrounding cells are used to determine the location and slope of the interface between two phases in each element, which is linearly approximated [23]. 2.2.2. Initial and boundary conditions Initially, the hydrogen droplet is at 20.0 K and at time t ¼ 0þ , the outermost boundary of the droplet in contact with the liquid helium suffers a temperature jump from 20.0 to 4.0 K, and is maintained at 4.0 K thereafter. 2.2.3. Solution technique The multiphase flow equations are solved by the phase-coupled, semi-implicit method for pressure-linked equation algorithm [22]. The initial time step is two orders of magnitude smaller than the characteristic time. Thereafter, larger time steps that satisfy the CFL (Courant–Friedrich–Levy) condition, may be used. At any time step, the solution is considered to have converged when (i) the normalized absolute residuals in the discrete equations summed over all elements are reduced by at least three orders of magnitude, and (ii) the absolute normalized mass imbalance summed over all elements is less than 10 3 .
3. Numerical results and discussion For the two-dimensional vaporization–condensation problem, the governing equations are solved by the finite volume method on an unstructured grid with 40,964 triangular elements for the case where the diameter of hydrogen droplet r is 2 mm. The grid was constructed with significant refinement so that the dynamic meshing scheme does not distort the grid much. For the twodimensional solidification problem, a mesh of 71,232 triangular elements is developed. A mesh of 30,091 tetrahedral elements is generated for the three-dimensional solidification case.
J. Xu et al. / Cryogenics 44 (2004) 507–514
511
Table 1 Physical parameters of helium and hydrogen [5–7] Parameter
Value Helium (liquid)
Helium (vapor)
Hydrogen (solid)
Hydrogen (liquid)
Density (kg/m3 ) Thermal conductivity (W/m K) Viscosity (·105 Pa s) Specific heat (J/kg K) Entropy (J/kg K) Boiling point (K) Latent heat (J/kg) Melting point (K) Surface tension (J/m2 ) Vapor pressure (Pa)
125.4 0.0177 0.1488 5170 3551 4.2 23,000 (vaporization)
0.7936 0.00387 0.0516 9033
92.0 0.18
71.0 1.82 0.84
14,304 20.3 350 (solidification) 13.8
0.00016 3129
The physical parameters relevant to the computation are listed in Table 1.
3.1. Helium vaporization and condensation Fig. 2 displays the helium vapor volume fraction contours at different instants of time. The first frame in Fig. 2 shows that a thin layer of helium vaporizes around the hydrogen droplet immediately after the introduction of the ‘‘hot’’ hydrogen droplet (20 K) into the bath. Frames 1–5 of Fig. 2 show that the helium vapor is quite dense initially around the droplet. As the droplet moves downward, helium vapor floats due to buoyancy forces and disperses in the helium bath. Fig. 3 shows the helium mixture temperature in the bath, corresponding to each frame in Fig. 2. It is apparent from two figures that the volume fraction of helium vapor decreases as it cools to the dew point and condenses into liquid helium. The mixture temperature finally reaches the temperature of the liquid helium bath.
In addition, the dispersion of helium vapor is strongly influenced by the helium flow field as well. Fig. 4 shows the flow field in the bath created by the wake of the moving droplet, vapor movement, and flow convection due to temperature gradients in the flow. The third frame in Fig. 4 demonstrates that there are two vortices following the moving droplet. As the droplet keeps moving, the helium flow is seriously disturbed, as shown in Frames 4 and 5 in Fig. 4. 3.2. Hydrogen droplet solidification Fig. 5 shows the evolution of the solidification process (of a spherical droplet of 2 mm diameter) from t ¼ 0:00 to 0.12 s, when the solidification is complete. The figure clearly demonstrates that (i) the phase change begins almost immediately after the hydrogen droplet contacts helium, and (ii) the phase front is a well-defined curve, virtually a circle. Fig. 6 displays the temperature contours within the hydrogen droplet. The contours are virtually circular too, meaning that the solidification
Fig. 2. The helium vapor volume fraction at different time.
512
J. Xu et al. / Cryogenics 44 (2004) 507–514
Fig. 3. Helium mixture temperature, corresponding to the times shown in Fig. 2.
Fig. 4. Velocity vectors within the helium fluid. Corresponding to the time steps shown in Figs. 2 and 3, this figure shows the velocity vectors of the helium fluid.
phase front is governed primarily by conduction effects, and convection has little influence. Fig. 7 shows the time required for solidification of hydrogen particles with different diameters. It is obvious that the smaller the particle, faster it freezes.
4. Summary and conclusions The studies show that a thin layer of helium vaporizes around the droplet immediately after it is injected into the helium medium. Due to the buoyancy force, helium vapor floats upwards. As the vaporized
helium cools, it condenses into liquid helium. The influence of the phase change and the droplet motion on the thermal and velocity fields is identified. Specifically, a graphic presentation of the temporal evolution of the helium volume fraction and mixture temperature is provided; also, snapshots of velocity vectors at various instants are provided to depict the evolution of the helium flow field. The studies on hydrogen solidification show that the outer-layer of the hydrogen droplet freezes almost immediately after it immerses in the liquid helium. Snapshots of the volume fraction of solidifying hydrogen and the associated temperature contours are shown at different time instants. Parametric study
J. Xu et al. / Cryogenics 44 (2004) 507–514
513
Fig. 5. Solidification of hydrogen droplet. Snapshots of the volume fraction of a solidifying hydrogen droplet are shown with times given in seconds. Red denotes liquid hydrogen, while blue represents solid hydrogen. It takes about 0.12 s for complete solidification of the droplet, whose diameter is 2.0 mm.
Fig. 6. Temperature contours within the solidifying hydrogen droplet. Elapsed time is given beneath each droplet, while the temperature scale is given at the bottom of the figure.
Solidification time (s)
0.14 0.12 0.1 0.08 0.06 0.04 0.02 0 0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
Diameter of droplets (mm)
Fig. 7. Solidification time for different size droplets.
2
2.2
shows that the smaller the hydrogen particle, the faster it freezes. A spherical hydrogen particle of diameter 2 mm freezes in 0.12 s, while it takes no more than 0.4 ms for complete solidification of a droplet of 0.1 mm. In conclusion, it has been demonstrated that the computational models discussed in this paper enable one to obtain a definitive qualitative picture of phase change phenomena associated with the freezing hydrogen particles in helium bath. These computational models can as well be used for quantitative analysis after validating with appropriate experimental results, and in fact, these models could help design proper experiments. Finally, there is room for improving these models. Our future goal is to couple the vaporization/condensation and
514
J. Xu et al. / Cryogenics 44 (2004) 507–514
solidification processes and check the validity of the current model. Acknowledgements This work was supported by the NASA Hydrogen Research program for Florida Universities (Contract No. NAG3-2751) and the Office of the Provost at The Florida State University. References [1] Palaszewski B. Launch vehicle performance with solid particle feed systems for atomic propellants. AIAA-98-3736 (NASA TM1998-208498). [2] Palaszewski B. Solid hydrogen experiments for atomic propellants. AIAA-2000-3855 (NASA TM-2001-211292). [3] Palaszewski B. Solid hydrogen experiments for atomic propellants: image analyses. AIAA-2001-3233 (NASA TM-2002211297). [4] Celik D, Van Sciver SW. Tracer particle generation in superfluid helium through cryogenic liquid injection for particle image velocimetry (PIV) applications. Exp Thermal Fluid Sci 2002;26(8): 971–5. [5] Van Sciver SW. Helium cryogenics. Plenum Press; 1986. [6] CRYODATA. Hepak, version 3.32. 1989. [7] Souers PC. Hydrogen properties for fusion energy. University of California Press; 1986. [8] Engelberg-Forster K. Heat transfer to a boiling liquid––mechanism and correlations. J Heat Transfer 1959;81:43–53. [9] Mikic BB, Rohsenow WM. A new correlation of pool-boiling data including the effect of heating surface characteristics. J Heat Transfer 1969;91:245–50. [10] Carey VP. Liquid–vapor phase-change phenomena: an introduction to the thermo-physics of vaporization and condensation processes in heat transfer equipment. Hemisphere Publishing Corporation; 1992.
[11] Stewart HB, Wendroff B. Two-phase flow: models and methods. J Comput Phys 1984;56:363–409. [12] Voller VR, Prakash C. A fixed-grid numerical modeling methodology for convection-diffusion mushy region phase-change problems. Int J Heat Mass Transfer 1987;30:1709–20. [13] Voller VR, Swaminathan CR. General source-based method for solidification phase change. Numer Heat Transfer B 1991; 19(2):175–89. [14] Voller VR, Markatos N, Cross M. Techniques for accounting for the moving interface in convection/diffusion phase change. In: Lewis RW, Morgan K, editors. Numerical methods in thermal problems (part 1), 1985. p. 595–609. [15] Brent AD, Voller VR, Reid KJ. Enthalpy-porosity technique for modeling convection–diffusion phase change: application to the melting of a pure metal. Numer Heat Transfer 1988;13:297– 318. [16] Zhao Y, Forhad A. A general method for simulation of fluid flows with moving and compliant boundaries on unstructured grids. Comput Methods Appl Mech Eng 2003;192:4439–66. [17] Lesoinne M, Farhat C. Geometric conservation laws for flow problems with moving boundaries and deformable meshes, and their impact on aeroelastic computations. Comput Methods Appl Mech Eng 1996;134:71–90. [18] Slone AK, Pericleous K, Bailey C, Cross M. Dynamic fluidstructure interaction using finite volume unstructured mesh procedures. Comput Struct 2002;80:371–90. [19] Degand C, Farhat C. A three-dimensional torsional spring analogy method for unstructured dynamic meshes. Comput Struct 2002;80:305–16. [20] Koobus B, Farhat C. Second-order time-accurate and geometrically conservative implicit schemes for flow computations on unstructured dynamic meshes. Comput Methods Appl Mech Eng 1999;170:103–29. [21] Manninen M, Taivassalo V, Kallio S. On the mixture model for multiphase flow. VTT Publications 288, Technical research center of Finland, Espoo, 1996. [22] Patankar SV. Numerical heat transfer and fluid flow. Hemisphere Publishing Corporation; 1980. [23] Youngs DL. Time-dependent multi-material flow with large fluid distortion. In: Morton KW, editor. Numerical methods for fluid dynamics. NY: Academic Press Inc.; 1982.
Hydrogen particles in liquid helium J. Xu a, K.M. Smith a, D. Celik b, M.Y. Hussaini a
a,*
, S.W. Van Sciver
b
School of Computational Science and Information Technology, Florida State University, 400 Dirac Science Library, Tallahassee, FL 32306, USA b National High Magnetic Field Laboratory, Florida State University, 1800 E. Paul Dirac Drive, Tallahassee, FL 32310, USA
Abstract Numerical investigations of phase changes during the generation of solid hydrogen particles in cryogenic helium, which include vaporization and re-condensation of helium around a hydrogen droplet and solidification of the droplet in helium, are conducted to unravel the basic mechanisms underlying the phase change phenomena. In a single-fluid, two-phase mixture framework, the simulation of vaporization/condensation involves a vapor–liquid exchange model that accounts for mass and energy transfers due to phase changes between liquid and vapor helium. The enthalpy-porosity, solidification–melting technique is employed to investigate the hydrogen solidification. The volume fraction of each phase in the solidification process is determined using a volume of fluid (VOF) two-phase model, along with a piecewise-linear approach to trace the liquid/solid interface. The results throw some light on the helium flow field resulting from the phase change and the motion of the droplet. They provide some insight into the evolution of liquid–vapor phase-exchange processes as well as the solidification evolution of hydrogen particles of various diameters. 2004 Elsevier Ltd. All rights reserved. Keywords: Solid hydrogen; Phase change; Solidification; Vaporization; Liquid helium; Single-fluid two-phase mixture
1. Introduction Hydrogen is the fuel used by the Space Shuttle Main Engine, as well as by many other hydrogen powered upper stages of rockets (Centaur, and those used by Delta, Atlas, and Ariane rocket upper stages, etc). However, due to its relative low density, larger storage tanks are required, which translates to higher structural weight and less payload. To reduce this unwanted effect, it has been proposed [1,2] that liquid hydrogen be refrigerated to near freezing temperatures and then pressurized to obtain dense, subcooled liquid fuel. Its density can be further increased by solidifying some portion of the liquid hydrogen in the form of small particles. Furthermore, atomic propellants, such as boron, carbon, or hydrogen can be embedded/stored in these small solid hydrogen particles, which can be stabilized and prevented from recombination [1–3], although such embedding techniques are yet to be developed. Thus, atomic hydrogen propellant feed systems may require the production and storage of hydrogen particles in liquid helium.
*
Corresponding author. Tel.: +1-850-6440601; fax: +1-8506451514. E-mail address: [email protected] (M.Y. Hussaini). 0011-2275/$ - see front matter 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.cryogenics.2004.02.017
Solid hydrogen particle formation in liquid helium has been experimentally investigated for years [1–4]. To create solid hydrogen particles in helium, droplets of liquid hydrogen at temperatures around 20 K are injected into a bath filled with liquid helium at 4 K. Since the temperature of the liquid hydrogen is higher than the boiling point of helium (4.2 K at 101,325 Pa [5,6]), and the temperature of ambient helium is lower than the solidification point of hydrogen (13.8 K at 101325 Pa [7]), vaporization, condensation and solidification occur simultaneously during the process. At the beginning, a thin layer of helium vaporizes around the hydrogen droplet, and the hydrogen droplet solidifies in the helium fluid. As the vaporized helium cools to the dew point of helium, it condenses in the helium bath, and reaches the temperature of the liquid helium bath in the end. Experimental efforts have been focused on the visual characterization of the particles, such as observing particle formation in liquid or vapor helium [4], video imaging of particle size [3], timing particle agglomeration [2,3], etc. Detailed experimental investigation of these phenomena is a formidable task due to very small time scales and other experimental design challenges. For example, Palaszewski [3] observed that in the hydrogen freezing process, large clouds of vapor ‘‘. . . would have also obscured the formation process, and
508
J. Xu et al. / Cryogenics 44 (2004) 507–514
thwarted efforts to see the final particles.’’ The authors encountered similar experimental difficulties. Numerical simulations, however, offer a relatively easy and less expensive alternative, and of course, they have to be validated by experimental observations. In addition, such simulations are becoming increasingly attractive due to advances in phase change theory [8–10], twophase models [11–15], and dynamic mesh techniques [16–20]. These numerical methodologies have been implemented in many phase change systems, but there are few numerical studies of the phenomena in helium– hydrogen cases. With the goal of better understanding the process of solid hydrogen particle generation, a series of computational studies have been conducted to investigate helium phase change around a hydrogen particle, and hydrogen droplet freezing in a helium bath. A vapor– liquid exchange model implemented in the framework of a single-fluid, two-phase mixture model is developed to study the vaporization/condensation process. An enthalpy-porosity, solidification–melting model is employed to investigate the solidification phenomena. Dynamic mesh techniques are introduced in the vaporization/condensation simulation to account for the movement of hydrogen droplets after injection. The effects of a phase changing, moving particle on the helium, as well as the evolution of the liquid–vapor helium phase exchange are studied. In addition, the time required for solidification of hydrogen particles with diameters ranging from 0.1 to 2.0 mm are investigated using the model developed.
2. Problem formulation and solution technique Formation of solid hydrogen particles in a helium bath involves a complex energy exchange processes between the two fluids resulting in vaporization of helium around the hydrogen droplet and its consequent condensation, and the simultaneous solidification of the liquid hydrogen droplet. The purpose of the present study is to study these processes in isolation and put them together in a later study, fully expecting that the sum will be greater than the whole as in any complex system. In other words, the ‘‘mushy layer’’ between the solidifying hydrogen droplet and the liquid helium is assumed as an infinitesimally thin interface and the interfacial conditions approximate/model the consequence of the mushy layer physics. Our study is the first cut at this extremely complex problem. 2.1. Vaporization and condensation of helium We assume that a hydrogen droplet of diameter r mm (r 2 ½0:1; 2:0) is injected vertically downwards with an initial velocity of ~ vy at the center of a rectangular bath of
Fig. 1. Simulation configuration. (The hydrogen droplet is shown as a red circular disk in the middle of the domain. The liquid helium is shown as blue, and is contained within the bath of width 8 mm and height 16 mm. Initially the helium temperature is 4 K, and the liquid hydrogen temperature is 20 K.)
dimensions 8 mm in width and 16 mm in depth, as shown in Fig. 1. The two-phase flow, comprising liquid and vapor, is modeled using a single-fluid two-phase approach [11,21]. In the two-phase approach, the phases are assumed to be interpenetrating, meaning that each phase has its own velocity vector field, and within any control volume there is a volume fraction of the primary (liquid) phase, and a volume fraction of the secondary (vapor) phase, which sum to one. 2.1.1. Governing equations The governing equations are the conservation equations for mass, momentum and energy of the mixture. It is expedient to consider these equations in a mesh system that moves with the hydrogen droplet. The continuity equation may be written as oqm þ r ðqm~ vÞ ¼ S_ q : ot
ð1Þ
ug , and ~ ug is the grid velocity of the Here, ~ v ¼~ vm ~ moving mesh, and ~ vm represents the mass-averaged velocity of the mixture, which is defined as ~ vm ¼
al ql~ vl þ av qv~ vv : qm
ð2Þ
The density of the mixture, qm , is defined to be qm ¼ al ql þ av qv ;
ð3Þ
where a represents the volume fraction, and subscripts l and v denote liquid and vapor phase, respectively. Following Engelberg-Forster and Greif [8,10], the mixture source/sink term, S_ q , in Eq. (1) is expressed as T Tb T Tb þ jc a v q v ; S_ q ¼ jv al ql Tb Tb
ð4Þ
J. Xu et al. / Cryogenics 44 (2004) 507–514
where T represents the temperature of ambient helium, Tb is the boiling point of helium, jv stands for the vaporization rate of liquid helium, while jc denotes the condensation rate of vapor helium. The momentum equation for the mixture is oqm~ v þ r ðqm~ vÞ þ qm~ v~ vÞ ¼ rp þ r ðlm r~ g ot ! X ak qk~ vdr;k~ vdr;k ; ð5Þ þ~ F þr k¼l;v
where lm ¼ al ll þ av lv , denotes the viscosity of the mixture; and ~ vdr;k ¼ ~ vkl ~ vm is the drift velocity for phase k, with ~ vkl ¼ ~ vk ~ vl the slip velocity of phase k, defined as the relative velocity of phase k over the primary phase l. Based on the assumption that a local equilibrium between the phases should be reached over spatial length scales, the slip velocity can be given by [21] ~ vkl ¼ skl~ a;
ð6Þ
where skl and ~ a are the relaxation time and acceleration of the secondary phase k, respectively. ~ F in Eq. (5) represents the body force, which is given by ~ F ¼ FD ð~ vm ~ ug Þ qm ;
ð7Þ
with the drag coefficient FD , given by a semi-empirical formula [4], FD ¼
24 ½1 þ 0:15Re0:687 ; p Rep
ð8Þ
where Rep is the particle Reynolds number. The energy equation for the mixture can be written in the form, X d X ðak qk hk Þ þ r ðak~ vk qk hk Þ ¼ r ðkeff rT Þ þ S_ h ; dt k¼l;v k¼l;v ð9Þ where hk is the sensible enthalpy of phase k, keff is the effective conductivity, and S_ h ¼ L Sq , with L the latent heat of helium. The volume fraction of the vapor (secondary) phase is governed by o ðav qv Þ þ r ðav qv~ vÞ ¼ r ðav qv~ vdr;v Þ þ S_ v ; ot where ( b ; if T > Tb ; vaporizes; jv al ql T T Tb S_ v ¼ T Tb jc av qv Tb ; if T < Tb ; condenses:
ð10Þ
ð11Þ
2.1.2. Initial and boundary conditions Initially, the helium bath and the liquid hydrogen droplet are at 4 and 20 K respectively, and the initial velocity of the droplet is 8 mm per second. During the
509
time the droplet freezes to 4 K, there is a complex interplay among the processes of vaporization, condensation and solidification. To study the vaporization/ condensation processes in isolation, we assume that (i) the shape and size of the droplet remains constant, (ii) the temperature evolution at the interface between helium and hydrogen droplet has a linear profile, T ¼ 20:0 c t, where c is the temperature change that the droplet suffers per unit time, and (iii) the droplet moves with y-velocity ~ vy ¼ 8:0 þ d t mm/s starting from the center of the bath. These assumptions are approximately supported by the results of droplet solidification simulation, which further determine c to be 0.4 K/s, and injection simulation, that gives d to be unity. Non-slip boundary conditions are applied at the walls of the bath. 2.1.3. Solution technique For computational purposes (easy management of the dynamic mesh), Eqs. (1), (5), (9) and (10) are transferred into a general integral form written as Z Z d qm /dV þ qm /~ v d~ A dt V Z ZoV ~ Cr/ dA þ S/ dV ; ð12Þ ¼ oV
V
where / represents a general transport quantity. In the continuity equation, / ¼ 1; in the momentum equation, / ¼~ v; and in the energy equation, / ¼ h (enthalpy). C is the diffusion, and S/ represents the source term for the quantity /, ~ A is the area vector, and V represents the control volume, which is identified with each mesh element. The time derivative term in the integral equations is discretized using a first-order backward difference formula, given by Z nþ1 n d ðq/V Þ ðq/V Þ ; ð13Þ q/ dV ¼ dt V Dt where the superscripts n and n þ 1 denote the respective quantity at the current and next time levels. The ðn þ 1Þth control volume V nþ1 is computed as V nþ1 ¼ V n þ
dV Dt; dt
ð14Þ
where dV =dt is the time derivative of the control volume. Following the grid conservation law, dV =dt is computed as Z nf X dV ~ ~ ¼ A¼ Aj ; ð15Þ ug d~ ug;j ~ dt oV j where nf is the number of faces on the control volume, ~ Aj is the j-face area vector and oV denotes the control volume surface. Dynamic mesh methodology [16–20] is employed to update the mesh where the shape of the
510
J. Xu et al. / Cryogenics 44 (2004) 507–514
domain is changing with respect to time due to the droplet motion. A linear spring-analogy smoothing method [19] and a local re-meshing method [18] are used in this simulation. A critical constraint on the moving mesh methodology is that the moving droplet should not advance more than one spatial element in a time step. In other words, if dy is the spatial scale in the direction of motion, dt is the time step, and vy is the moving droplet velocity, then the constraint dy > vy dt should be maintained at each time step. A second-order upwind scheme is used for spatial discretization of the momentum equation and the continuity equations. The multiphase flow equations are solved using the coupled semi-implicit method for pressure-linked equation algorithm [22]. 2.2. Solidification of hydrogen droplets 2.2.1. Governing equations For the problem of a static liquid hydrogen droplet immersed in a helium bath, convection and diffusion are the principal mechanisms of heat transfer. The R T energy equation in terms of the sensible enthalpy, h ¼ Tref cp dT , is oqh þ r ðq~ vhÞ ¼ r ðarhÞ þ Sh ; ot
ð16Þ
where q is the density of liquid or solid hydrogen, ~ v represents the velocity, where ~ v ¼ 0 within the solid and at the interface of the liquid and solid, a is the thermal diffusivity, and Sh is the source term that accounts for the latent heat of phase change. An enthalpy-porosity model proposed by Voller et al. [12–14] and modified by Brent et al. [15], is employed to represent the source term: Sh ¼
oðqDH Þ ; ot
ð17Þ
where DH , the latent heat content, can be written as a function of temperature T , 0; if T < Tm ; DH ¼ ð18Þ L; if T > Tm ; where L is the latent heat (of hydrogen), and Tm is the phase change temperature of hydrogen solidification/ melting. The latent heat is treated as a heat source or sink by assigning a nodal latent heat value to each computational element [13], according to the temperature of the cell. For the solid/liquid two-phase problem, a volume of fluid (VOF) model [11] is implemented to solve the volume fraction of each phase. The essential premise of the VOF model is the assumption that the two phases are not interpenetrating, which is true in the solidification case. Interface tracking between phases is accomplished by solving a continuity equation for the volume
fraction of each phase. The volume fraction is given by ai , where i ¼ s for the solid phase, and i ¼ l for the liquid phase. The volume fraction continuity equations are given by oas þ~ v ras ¼ 0: ot
ð19Þ
A single momentum equation is solved throughout the domain, given by o ðq~ vÞ þ r ðq~ v~ vÞ ¼ r ½lðr~ vÞ þ q~ g þ Sb ; ot
ð20Þ
where the term Sb ¼ q~ gbðT T0 Þ accounts for the buoyancy force of natural convection effects, and b is the thermal expansion coefficient. Volume fractions of phases for the element under consideration and its surrounding cells are used to determine the location and slope of the interface between two phases in each element, which is linearly approximated [23]. 2.2.2. Initial and boundary conditions Initially, the hydrogen droplet is at 20.0 K and at time t ¼ 0þ , the outermost boundary of the droplet in contact with the liquid helium suffers a temperature jump from 20.0 to 4.0 K, and is maintained at 4.0 K thereafter. 2.2.3. Solution technique The multiphase flow equations are solved by the phase-coupled, semi-implicit method for pressure-linked equation algorithm [22]. The initial time step is two orders of magnitude smaller than the characteristic time. Thereafter, larger time steps that satisfy the CFL (Courant–Friedrich–Levy) condition, may be used. At any time step, the solution is considered to have converged when (i) the normalized absolute residuals in the discrete equations summed over all elements are reduced by at least three orders of magnitude, and (ii) the absolute normalized mass imbalance summed over all elements is less than 10 3 .
3. Numerical results and discussion For the two-dimensional vaporization–condensation problem, the governing equations are solved by the finite volume method on an unstructured grid with 40,964 triangular elements for the case where the diameter of hydrogen droplet r is 2 mm. The grid was constructed with significant refinement so that the dynamic meshing scheme does not distort the grid much. For the twodimensional solidification problem, a mesh of 71,232 triangular elements is developed. A mesh of 30,091 tetrahedral elements is generated for the three-dimensional solidification case.
J. Xu et al. / Cryogenics 44 (2004) 507–514
511
Table 1 Physical parameters of helium and hydrogen [5–7] Parameter
Value Helium (liquid)
Helium (vapor)
Hydrogen (solid)
Hydrogen (liquid)
Density (kg/m3 ) Thermal conductivity (W/m K) Viscosity (·105 Pa s) Specific heat (J/kg K) Entropy (J/kg K) Boiling point (K) Latent heat (J/kg) Melting point (K) Surface tension (J/m2 ) Vapor pressure (Pa)
125.4 0.0177 0.1488 5170 3551 4.2 23,000 (vaporization)
0.7936 0.00387 0.0516 9033
92.0 0.18
71.0 1.82 0.84
14,304 20.3 350 (solidification) 13.8
0.00016 3129
The physical parameters relevant to the computation are listed in Table 1.
3.1. Helium vaporization and condensation Fig. 2 displays the helium vapor volume fraction contours at different instants of time. The first frame in Fig. 2 shows that a thin layer of helium vaporizes around the hydrogen droplet immediately after the introduction of the ‘‘hot’’ hydrogen droplet (20 K) into the bath. Frames 1–5 of Fig. 2 show that the helium vapor is quite dense initially around the droplet. As the droplet moves downward, helium vapor floats due to buoyancy forces and disperses in the helium bath. Fig. 3 shows the helium mixture temperature in the bath, corresponding to each frame in Fig. 2. It is apparent from two figures that the volume fraction of helium vapor decreases as it cools to the dew point and condenses into liquid helium. The mixture temperature finally reaches the temperature of the liquid helium bath.
In addition, the dispersion of helium vapor is strongly influenced by the helium flow field as well. Fig. 4 shows the flow field in the bath created by the wake of the moving droplet, vapor movement, and flow convection due to temperature gradients in the flow. The third frame in Fig. 4 demonstrates that there are two vortices following the moving droplet. As the droplet keeps moving, the helium flow is seriously disturbed, as shown in Frames 4 and 5 in Fig. 4. 3.2. Hydrogen droplet solidification Fig. 5 shows the evolution of the solidification process (of a spherical droplet of 2 mm diameter) from t ¼ 0:00 to 0.12 s, when the solidification is complete. The figure clearly demonstrates that (i) the phase change begins almost immediately after the hydrogen droplet contacts helium, and (ii) the phase front is a well-defined curve, virtually a circle. Fig. 6 displays the temperature contours within the hydrogen droplet. The contours are virtually circular too, meaning that the solidification
Fig. 2. The helium vapor volume fraction at different time.
512
J. Xu et al. / Cryogenics 44 (2004) 507–514
Fig. 3. Helium mixture temperature, corresponding to the times shown in Fig. 2.
Fig. 4. Velocity vectors within the helium fluid. Corresponding to the time steps shown in Figs. 2 and 3, this figure shows the velocity vectors of the helium fluid.
phase front is governed primarily by conduction effects, and convection has little influence. Fig. 7 shows the time required for solidification of hydrogen particles with different diameters. It is obvious that the smaller the particle, faster it freezes.
4. Summary and conclusions The studies show that a thin layer of helium vaporizes around the droplet immediately after it is injected into the helium medium. Due to the buoyancy force, helium vapor floats upwards. As the vaporized
helium cools, it condenses into liquid helium. The influence of the phase change and the droplet motion on the thermal and velocity fields is identified. Specifically, a graphic presentation of the temporal evolution of the helium volume fraction and mixture temperature is provided; also, snapshots of velocity vectors at various instants are provided to depict the evolution of the helium flow field. The studies on hydrogen solidification show that the outer-layer of the hydrogen droplet freezes almost immediately after it immerses in the liquid helium. Snapshots of the volume fraction of solidifying hydrogen and the associated temperature contours are shown at different time instants. Parametric study
J. Xu et al. / Cryogenics 44 (2004) 507–514
513
Fig. 5. Solidification of hydrogen droplet. Snapshots of the volume fraction of a solidifying hydrogen droplet are shown with times given in seconds. Red denotes liquid hydrogen, while blue represents solid hydrogen. It takes about 0.12 s for complete solidification of the droplet, whose diameter is 2.0 mm.
Fig. 6. Temperature contours within the solidifying hydrogen droplet. Elapsed time is given beneath each droplet, while the temperature scale is given at the bottom of the figure.
Solidification time (s)
0.14 0.12 0.1 0.08 0.06 0.04 0.02 0 0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
Diameter of droplets (mm)
Fig. 7. Solidification time for different size droplets.
2
2.2
shows that the smaller the hydrogen particle, the faster it freezes. A spherical hydrogen particle of diameter 2 mm freezes in 0.12 s, while it takes no more than 0.4 ms for complete solidification of a droplet of 0.1 mm. In conclusion, it has been demonstrated that the computational models discussed in this paper enable one to obtain a definitive qualitative picture of phase change phenomena associated with the freezing hydrogen particles in helium bath. These computational models can as well be used for quantitative analysis after validating with appropriate experimental results, and in fact, these models could help design proper experiments. Finally, there is room for improving these models. Our future goal is to couple the vaporization/condensation and
514
J. Xu et al. / Cryogenics 44 (2004) 507–514
solidification processes and check the validity of the current model. Acknowledgements This work was supported by the NASA Hydrogen Research program for Florida Universities (Contract No. NAG3-2751) and the Office of the Provost at The Florida State University. References [1] Palaszewski B. Launch vehicle performance with solid particle feed systems for atomic propellants. AIAA-98-3736 (NASA TM1998-208498). [2] Palaszewski B. Solid hydrogen experiments for atomic propellants. AIAA-2000-3855 (NASA TM-2001-211292). [3] Palaszewski B. Solid hydrogen experiments for atomic propellants: image analyses. AIAA-2001-3233 (NASA TM-2002211297). [4] Celik D, Van Sciver SW. Tracer particle generation in superfluid helium through cryogenic liquid injection for particle image velocimetry (PIV) applications. Exp Thermal Fluid Sci 2002;26(8): 971–5. [5] Van Sciver SW. Helium cryogenics. Plenum Press; 1986. [6] CRYODATA. Hepak, version 3.32. 1989. [7] Souers PC. Hydrogen properties for fusion energy. University of California Press; 1986. [8] Engelberg-Forster K. Heat transfer to a boiling liquid––mechanism and correlations. J Heat Transfer 1959;81:43–53. [9] Mikic BB, Rohsenow WM. A new correlation of pool-boiling data including the effect of heating surface characteristics. J Heat Transfer 1969;91:245–50. [10] Carey VP. Liquid–vapor phase-change phenomena: an introduction to the thermo-physics of vaporization and condensation processes in heat transfer equipment. Hemisphere Publishing Corporation; 1992.
[11] Stewart HB, Wendroff B. Two-phase flow: models and methods. J Comput Phys 1984;56:363–409. [12] Voller VR, Prakash C. A fixed-grid numerical modeling methodology for convection-diffusion mushy region phase-change problems. Int J Heat Mass Transfer 1987;30:1709–20. [13] Voller VR, Swaminathan CR. General source-based method for solidification phase change. Numer Heat Transfer B 1991; 19(2):175–89. [14] Voller VR, Markatos N, Cross M. Techniques for accounting for the moving interface in convection/diffusion phase change. In: Lewis RW, Morgan K, editors. Numerical methods in thermal problems (part 1), 1985. p. 595–609. [15] Brent AD, Voller VR, Reid KJ. Enthalpy-porosity technique for modeling convection–diffusion phase change: application to the melting of a pure metal. Numer Heat Transfer 1988;13:297– 318. [16] Zhao Y, Forhad A. A general method for simulation of fluid flows with moving and compliant boundaries on unstructured grids. Comput Methods Appl Mech Eng 2003;192:4439–66. [17] Lesoinne M, Farhat C. Geometric conservation laws for flow problems with moving boundaries and deformable meshes, and their impact on aeroelastic computations. Comput Methods Appl Mech Eng 1996;134:71–90. [18] Slone AK, Pericleous K, Bailey C, Cross M. Dynamic fluidstructure interaction using finite volume unstructured mesh procedures. Comput Struct 2002;80:371–90. [19] Degand C, Farhat C. A three-dimensional torsional spring analogy method for unstructured dynamic meshes. Comput Struct 2002;80:305–16. [20] Koobus B, Farhat C. Second-order time-accurate and geometrically conservative implicit schemes for flow computations on unstructured dynamic meshes. Comput Methods Appl Mech Eng 1999;170:103–29. [21] Manninen M, Taivassalo V, Kallio S. On the mixture model for multiphase flow. VTT Publications 288, Technical research center of Finland, Espoo, 1996. [22] Patankar SV. Numerical heat transfer and fluid flow. Hemisphere Publishing Corporation; 1980. [23] Youngs DL. Time-dependent multi-material flow with large fluid distortion. In: Morton KW, editor. Numerical methods for fluid dynamics. NY: Academic Press Inc.; 1982.