Numerical prediction of flow characteristics of slush hydrogen in a horizontal pipe

Numerical prediction of flow characteristics of slush hydrogen in a horizontal pipe

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y x x x ( 2 0 1 6 ) 1 e1 2 Available online at www.sciencedirect.com ScienceDi...

2MB Sizes 0 Downloads 53 Views

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y x x x ( 2 0 1 6 ) 1 e1 2

Available online at www.sciencedirect.com

ScienceDirect journal homepage: www.elsevier.com/locate/he

Numerical prediction of flow characteristics of slush hydrogen in a horizontal pipe T. Jin a,c,*, Y.J. Li a, Z.B. Liang a, Y.Q. Lan b, G. Lei c, X. Gao c a

Institute of Refrigeration and Cryogenics of Zhejiang University/Key Laboratory of Refrigeration and Cryogenic Technology of Zhejiang Province, Hangzhou 310027, China b Beijing Institute of Aerospace Testing Technology, Beijing 100074, China c State Key Laboratory of Technologies in Space Cryogenic Propellants, Beijing 100028, China

article info

abstract

Article history:

Slush hydrogen has lower temperature, higher density and higher heat capacity than those

Received 10 April 2016

of liquid hydrogen, and is considered as a potential propellant for novel aerospace rockets

Received in revised form

with decreased spacecraft size and weight. In this study, an improved three-dimensional

8 September 2016

numerical model based on EulereEuler two-fluid model has been built to predict the

Accepted 9 September 2016

flow characteristics of slush hydrogen in a horizontal pipe. In this model, an effective

Available online xxx

viscosity of mixture, which takes the particle shape and size into consideration, is adopted to modify the drag law for interphase momentum exchange, and the wall boundary con-

Keywords:

ditions for the solid phase are based on JohnsoneJackson model which involves the friction

Slush hydrogen

and collision between the particle and the wall. The performance of the model has been

Two-phase flow

verified by the comparison between the calculated results and the experimental data from

Flow characteristics

the literatures and considered to be effective for slush hydrogen flow. The improved model is then used to analyze the effects of inlet velocity, solid fraction, particle size on the flow characteristics, including pressure gradient, solid volume fraction distribution and velocity distribution and abundant to predict the flow pattern shift. Moreover, the numerical results indicate that the pressure drops for subcooled liquid hydrogen, under some operating conditions, are greater than those of slush hydrogen, which are presented in some published experimental work. © 2016 Hydrogen Energy Publications LLC. Published by Elsevier Ltd. All rights reserved.

Introduction Slush hydrogen is the cryogenic suspension of solid hydrogen particles and subcooled liquid hydrogen, which has been considered as a potential propellant for novel aerospace rockets with decreased spacecraft size and weight due to its higher density and heat capacity than normal-boiling-point liquid hydrogen [1e4]. An important concern for the

application of slush hydrogen is the transfer characteristics, especially for high concentration slurries. Earlier experimental studies of NASA National Aerospace Plane (NASP) project showed that the slush hydrogen with 50% solid fraction could well flow through pipelines, valves, turbo-pumps and other flow components [5,6]. Another kind of cryogenic suspension, slush nitrogen, also have attracted attention for its potential as the coolant for high temperature

* Corresponding author. Institute of Refrigeration and Cryogenics of Zhejiang University/Key Laboratory of Refrigeration and Cryogenic Technology of Zhejiang Province, Hangzhou 310027, China. E-mail address: [email protected] (T. Jin). http://dx.doi.org/10.1016/j.ijhydene.2016.09.054 0360-3199/© 2016 Hydrogen Energy Publications LLC. Published by Elsevier Ltd. All rights reserved. Please cite this article in press as: Jin T, et al., Numerical prediction of flow characteristics of slush hydrogen in a horizontal pipe, International Journal of Hydrogen Energy (2016), http://dx.doi.org/10.1016/j.ijhydene.2016.09.054

2

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y x x x ( 2 0 1 6 ) 1 e1 2

superconductive cables [7e9]. Ohira K [8] found that the pressure drop reduction phenomenon (i.e., pressure drop for subcooled liquid nitrogen is greater than that of slush nitrogen) emerged clearly for high velocity flow from experimental results of slush nitrogen flow in horizontal pipe. With the development of computational technology, CFD (Computational Fluid Dynamics) method is becoming an effective way to predict the multiphase flow. In 1990s, a numerical code, called FLUSH, was developed to calculate the pressure drop and solid fraction loss for steady-state and one-dimensional flow of slush hydrogen through the pipe system in the frame of the NASP project and achieved good agreement with the preliminary experimental results [5]. Gamma et al. [10] used the finite-volume Navier-Stokes CFD solver to get a precise assessment of thermal fluid dynamics behavior of slush hydrogen for the Europe FESTIP program and took in consideration several situations, namely linear pipelines flow, Venturi duct flow and in-tank storage. However, so far the numerical model for slush hydrogen is still not effective for fully understanding the cryogenic flow mechanism. The CFD models employed for the pipe flows of ordinary solideliquid mixtures have been improved and applicable to determine the pressure gradient, volume fraction distribution and velocity profile. Nevertheless, these approaches, usually with dispersed phase of sand, spherical glass beads, ash, etc., are valid for fine and uniform particles, but not suitable for cryogenic slush hydrogen, since the latter carries nonspherical and coarser particles with an average size of 0.5e2 mm. Moreover, the modeling of cryogenic slurry flows should consider the interphase heat transfer and heat leak, which results in solid fraction loss. Ishimoto et al. [11] developed a two-dimensional model based on the unsteady thermal non-equilibrium EulerianeLagrangian approach to predict the flow characteristics of slush nitrogen flow in horizontal circular and converging-diverging pipes, but this approach is not applicable for the dense slurry flow. Several three-dimensional CFD simulations based on steady EulereEuler approach have also been conducted to determine the flow and heat transfer features of slush nitrogen [12,13] and slush hydrogen [14]. However, more attention needs to be paid to the effects of coarser cryogenic particles on the phase interaction for momentum and energy exchange, which are different from that of fine ambient particles. In the present study, an improved three-dimensional (3-D) computational model is built up and the basic equations are adopted for predicting the flow characteristics of slush hydrogen in a horizontal pipe. The novelty of the model resides in the consideration of the large-size solid hydrogen particles. By introducing the mixture effective viscosity, which is dependent on the characteristics of the solid particles, the effects of particles on the interfacial momentum exchange are taken into account. The wall boundary conditions for solid phase based on the JohnsoneJackson correlations [15] are employed to simulate the friction and collision between the particle and the wall, and the effective empirical parameters are adopted. The proposed model is validated by various sets of experimental data over a large scale operating conditions in the literatures, in terms of slurry type, mean flow velocity, particle size and solid concentration, and shows fairly good performance. Thus, the model is adopted to

calculate and analyze the effects of inlet velocity, solid fraction, particle size on the flow features, including pressure gradient, solid volume fraction distribution and velocity profiles.

Mathematical model The EulereEuler two-fluid approach allows for the modeling of multiple separate, yet interacting phases. In the 3-D numerical modeling of cryogenic slush hydrogen, the flows are considered incompressible and steady and each particle is inelastic with uniform and constant size and shape. The thermo-physical property values at the triple point adopted in the simulations for slush hydrogen are given in Table 1 [16]. The energy exchange between two phases is accounted, yet the variation in the properties and mass due to the temperature and pressure changes have not been considered, nor does the vaporization.

Governing equations Conservation of mass    v aq rq þ V$ aq rq ! v q ¼ m_ pq  m_ qp vt

(1)

where the subscripts p and q denote for either “s” referring to the solid phase or “l” referring to the liquid phase, respectively, and p is the contrary phase of q.

Conservation of momentum    v v q þ V, aq rq ! v q! v q ¼ aq Vp þ V,tq þ aq rq ! g  Vps aq rq ! vt ! ! ! þ F L;q þ F VM;q þ F D;q   þ m_ ! v  m_ ! v pq

pq

qp

qp

(2) where ps is the solids pressure composed of a kinetic term and a second term due to the particle collisions, as given by ps ¼ as rs Qs þ 2as 2 rs Qs ds g0;ss ð1 þ ess Þ. ess is the restitution coefficient for collisions between particles, taken as 0.7e0.99 for different types of particles. ds is the particle diameter put as 0.5e2 mm for slush hydrogen. Qs is the granular temperature and g0;ss is the radial distribution function. The lift force due to solideliquid velocity gradients is expressed as

Table 1 e Thermo-physical properties of hydrogen at triple point. Thermo-physical property Pressure Temperature Density Dynamic viscosity Specific heat Thermal conductivity Latent heat of fusion

Unit

Solid hydrogen

Liquid hydrogen

MPa K kg/m3 Pa$s kJ/kg$k W/(m$K) kJ/kg

e 13.957 86.59 e 2.6 0.129 62.8

0.00736 13.957 77.05 2.58E-05 7.0197 0.104 e

Please cite this article in press as: Jin T, et al., Numerical prediction of flow characteristics of slush hydrogen in a horizontal pipe, International Journal of Hydrogen Energy (2016), http://dx.doi.org/10.1016/j.ijhydene.2016.09.054

3

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y x x x ( 2 0 1 6 ) 1 e1 2

! F L;q ¼ CL as rl ð! vl ! v s Þ  ðV  ! v l Þ. CL is the lift coefficient, calculated by the SaffmaneMei model which is suitable for ffiC 0 , where the vorticity spherical solid particles: CL ¼ 2pp3ffiffiffiffiffiffi Rew L ! ! Reynolds number Rew ¼ rl j v l ml v s jds , and CL 0 ¼ 6:46[17]. The virtual mass force due to acceleration motion is !  ! ! v s . The virtual mass expressed as F VM ¼ CVM as rl d dtv l  d dt coefficients for the model are CVM ¼ 0:5[18]. The interfacial drag force between the solid and liquid phases, is expressed by the HuilineGidaspow model [19], which is a combination of the Wen and Yu model [20], and the Ergun equation [21]. The drag coefficient CD is given by Schiller and Naumann formula [22] 8   > < 24 1 þ 0:15Rep 0:687 CD ¼ Rep > : 0:44

Re  1000

(3)

Re > 1000

where Rep is the particle Reynolds number, defined as v q =mm . mm is the effective viscosity of the vp ! Rep ¼ rl ds ! mixture. There are some typical empirical formulas for the mixture effective viscosity [23e27], as shown in Table 2, using a fixed set of empirical coefficients as;max and ½h. The packing limit as;max specifies the maximum volume fraction for the granular phase, the intrinsic viscosity ½h characterizes the particles' contribution to the mixture viscosity. For rigid spheres, as;max ¼ 0:63 and ½h ¼ 2:5. Nevertheless, these two empirical parameters vary markedly from the particle shape size and type, which resulted in the complicated mathematical determinations and coordination. The method for the effective viscosity of mixture used in this study is an exponential formula developed by Cheng and Law [28], extended and derived from the Einstein's formula, as given by ( #) " 2:5 1  1 mm ¼ ml exp b ð1  as Þb

(9)

where b is the only parameter to account for the influence of particle shape, size and type, and its value varies from 0.95 to 3.9.

Granular temperature The transport equation derived for the solids phase based on the granular temperature was described by Ding and Gidaspow [29]:

  3 v ðas rs Qs Þ þ Vðas rs ! v s Qs Þ ¼  ps I þ ts 2 vt : V! v s þ V,ðkQs VQs Þ  gQs þ fls (12)

where the first term on the right side is the generation of energy by the solid stress tensor, and the second term describes the diffusive flux of granular energy. kQs is the diffusion coefficient, gQs is the collisional dissipation of energy within the particles, and fls is the energy exchange between liquid phase and solid phase. kQs and fls can be described as follows: kQs ¼

pffiffiffiffiffiffiffi

2 150rs ds pQ 6 1 þ as g0;ss ð1 þ ess Þ 5 384g0;ss ð1 þ ess Þ rffiffiffiffiffiffi Qs 2 þ 2as rs ds g0;ss ð1 þ ess Þ p

fls ¼ 3Kls Qs

(13) (14)

and gQs was described by Lun et al. [30]: gQs ¼

  12 1  ess 2 g0;ss pffiffiffi rs as 2 Qs 3=2 ds p

(15)

Turbulence modeling and wall boundary conditions In the present CFD analysis, the standard k  ε model is used for predicting the turbulent quantities. The turbulence kinetic energy k and its rate of dissipation ε are obtained from

  mt;l vðal rl kÞ þ V,ðal rl ! v l kÞ ¼ V, ml þ Vk þ Gk  al rl ε þ Pk vt sk (16)

  mt;l vðal rl εÞ ε ε2 þ V,ðal rl ! v l εÞ ¼ V, ml þ Vε þ C1ε Gk  C2ε al rl vt k sε k þ Pε (17) 2 rl Cm kε

is the turbulent viscosity, Gk represents the where mt;f ¼ generation of turbulence kinetic energy, Pk and Pε are the source terms which can be included to the model for the turbulent interaction between the dispersed phases and the continuous phase, the model constants are set to C1ε ¼ 1:44, C2ε ¼ 1:92, Cm ¼ 0:09, sk ¼ 1:0, and sε ¼ 1:3[27]. For each phase, the source terms are v pq ,! v dr Pk;q ¼ Cs aq Kpq !

(18)

εk Pε ¼ C3ε Pk k

(19)

The wall boundary condition for the particle phase is based on the JohnsoneJackson's correlation. The shear stress for the particles at the wall is in the following form

Table 2 e Empirical formulas for effective viscosity in the literatures. Literature

Empirical formulas

Einstein [23] Mooney [24]

mm ¼ ml ð1 þ ½has Þ (4)   s mm ¼ ml exp 1a½ha (5) =a s s;max

Thomas & Muthukumar [25]

mm ¼ ml ð1 þ 2:5as þ 10:05as 2 þ 0:00273 expð16:6as ÞÞ (6)

Mori et al. [26]

mm ¼ ml ð1 þ 3ðas 1  as;max 1 ÞÞ (7)

Barnes et al. [27]

mm ¼ ml expð1  as =as;max Þas;max ½h (8)

Please cite this article in press as: Jin T, et al., Numerical prediction of flow characteristics of slush hydrogen in a horizontal pipe, International Journal of Hydrogen Energy (2016), http://dx.doi.org/10.1016/j.ijhydene.2016.09.054

4

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y x x x ( 2 0 1 6 ) 1 e1 2

! p pffiffiffi as ts ¼  r g0;ss Qs 1=2 U sw 3f 6 as;max s

(20)

The general boundary condition of heat flux for the granular phase at the wall is calculated by p pffiffiffi as p pffiffiffi as ! ! ! qs ¼ r g0;ss Qs 1=2 U sw , U sw  3f 3 6 4 as;max s as;max   3=2 2 1  esw rs g0;ss Qs

(21)

! where U sw is the particle slip velocity parallel to the wall. f is the specularity coefficient between the particle and the wall, whose value ranges from 0 to unity depending on the surface conditions of both the wall and particles, set to 0.0001e0.1. esw is the restitution coefficient for the collisions between the particle and the wall, set to 0.7e0.99. The values of parameters for granular temperature is suggested to f ¼ 0:0001, ess ¼ 0:9 and esw ¼ 0:99 for slush nitrogen [13] and f ¼ 0:03, ess ¼ 0:9 and esw ¼ 0:7 for slush hydrogen [14].

Numerical domain and simulation set-up On the basis of the experimental setup for slush hydrogen of NASP project, the transfer duct modeled in this study is a horizontal circular pipe with the diameter D of 45 mm and the length L of 3 m which should meet the demand of L  50D for fully flow development. The computational grid of the analysis domain is an unstructured grid generated by ANSYS ICEM CFD 14.5 and the near-wall mesh is refined. The schematic of the computed system and the computational grid for the numerical analysis is given in Fig. 1. To determine the boundary conditions for the flow domain, the constant velocity and volume fraction for each phase are applied to the flow inlet. The average particle size of slush hydrogen is at the range of 0.5e2 mm, the inlet velocity is 1e4 m/s and the volume fraction for the solid phase is 5e35%. The outflow conditions are confined with an ambient outflow pressure. For the wall boundary conditions, no slip is assumed for the liquid phase, and JohnsoneJackson correlations are adopted for the granular phase. Referring to the suggested value [13,14], the specularity coefficient, particleeparticle and particle-wall restitution coefficients for slush hydrogen flow are f ¼ 0:01, ess ¼ 0:9 and esw ¼ 0:7, respectively. The effective viscosity parameter for slush hydrogen is b ¼ 1:0. Also, the adiabatic condition is applied for the wall.

A very fine grid scheme increases the accuracy and stability of multiphase flow but results in high computational expense. Thus, the grid independence investigation has been conducted, where the mesh is set with the cell quantities of 60 k, 378 k, 718 k, and 1809 k, respectively. As an example, the grid-dependence results for the case of inlet velocity of 1 m/s, solid fraction of 13% and particle size of 1 mm are given in Fig. 2. The velocity profiles of solid phase show the consistent trend and the maximum difference is 2.7% for the cell quantities over 378 k. Nevertheless, there is a front nose in the velocity distribution which indicates the particles are accelerated around the pipe center. It is considered to be due to interference among the particles and between particles and wall in the vicinity of pipe wall for flow at low velocity and light solid concentration. The maximum solid concentration is very close for cell quantities over 378 k with a relative deviation of 0.6%, and the greatest difference for solid fraction profiles is at the pipe bottom with a deviation of 4.6%. Hence, the mesh grid for the present calculation is set with a cell quantity of 378 k. The grid for the scheme is fined at the wall boundary layer with the minimum width near the pipe wall of 0.3 mm, and the grid in the flow direction is even. The standard k  ε model is used for predicting the turbulent quantities. The dimensionless wall distance for the cell adjacent to the wall yþ fulfills the requirement of standard wall function 30 < yþ < 70.

Results and discussion Model validation Due to the lack of experimental data on the flow characteristics of slush hydrogen, some investigations on the ordinary slurry flows from the literatures [8,31] are summarized in Table 3 for model validation. The reason to choose these experimental cases is that they include comprehensive consideration of the specific characteristics, such as pressure drop, solid fraction distribution and velocity profiles, regarding the cryogenic slurry. The cases are analyzed with the present mathematical model and the numerical results will then be compared with the corresponding experimental results. For Cases 1e3, the empirical parameters in the model adopted for sand/water, namely specularity coefficient, effective viscosity parameter, particleeparticle and particle-

Fig. 1 e Schematic of the computed system and the computational grid for the numerical analysis. Please cite this article in press as: Jin T, et al., Numerical prediction of flow characteristics of slush hydrogen in a horizontal pipe, International Journal of Hydrogen Energy (2016), http://dx.doi.org/10.1016/j.ijhydene.2016.09.054

5

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y x x x ( 2 0 1 6 ) 1 e1 2

0.5

0.5

0.4 0.3 0.2

0.3 0.2 0.1

0.0

y/D

y/D

0.1

60k 377k 718k 1809k

0.4

60k 377k 718k 1809k

-0.1

0.0 -0.1

-0.2

-0.2

-0.3

-0.3

-0.4

-0.4

-0.5

-0.5

0.0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

1.6

0.0

0.1

Velocity (m/s)

0.2

0.3

0.4

0.5

Solid fraction (%)

(a) velocity distribution of solid phase

(b) solid fraction profiles

Fig. 2 e Effect of the cell quantities on the numerical results for slush hydrogen flow.

Table 3 e Experimental slurry flow cases in the literatures. Ref.

Slurries

D (mm)

dp (mm)

rl (kg/m3)

rs (kg/m3)

C (%)

V (m/s)

Roco & Shook [31] Roco & Shook [31] Roco & Shook [31] Ohira [8] Ohira [8] Ohira [8]

Sand, water Sand, water Sand, water Slush nitrogen Slush nitrogen Slush nitrogen

51.5 51.5 51.5 15 15 15

0.165 0.165 0.165 1.36 1.36 1.36

1000 1000 1000 867.86 867.86 867.86

2650 2650 2650 1026.5 1026.5 1026.5

8.41 9.18 28.6 12.2 21.6 11e15 wt

1.66 3.78 4.33 1.68 4.1 1e4

Case ID C1 C2 C3 C4 C5 C6

1.0

1.0 C1 [Expt][28] C1 [CFD]

0.8

C1 [Expt][28] C1 [CFD] C3 [Expt][28] C3 [CFD]

0.8

C2 [Expt][28] C2 [CFD]

0.6

y/D

y/D

0.6 0.4

0.2

0.2 0.0

0.4

0

1

2

3

4

5

6

7

8

Velocity Distribution (m/s) (a) solid velocity distribution

0.0

0

5

10 15 20 25 30 35 40 45 50 55

Solid Fraction (%) (b) solid fraction profiles

Fig. 3 e Comparison between numerical data and experimental results of Cases 1e3.

wall restitution coefficients, are f ¼ 0:0001, b ¼ 3, ess ¼ 0:9 and ew ¼ 0:99, respectively. As presented in Fig. 3, the numerical results based on the present model are compared with the above experimental results to examine the applicability of model. Fig. 3(a) shows the solid velocity distribution results of Cases 1 and 2, with various inlet velocities. The results of CFD and experiments agree well. At the low velocity of 1.66 m/s in Case 1, the solid velocity at the upper part of pipe is significantly higher than that at the lower part, and the highest solid velocity locates at the upper part of pipe, which has good agreement with the experimental data. One of the possible reasons for the

asymmetry is that the forces to keep the solid suspending are slight at low velocity. At the high velocity of 3.78 m/s, the symmetry is obvious, and the highest solid velocity moves down to the nearby center of the pipe. Fig. 3(b) presents the solid fraction profiles of Case 1 and Case 3 to investigate the applicability of numerical model for the high concentration and velocity flow. The solid fraction and inlet velocity of Case 3 are much higher than the others’. For Case 1, the solid particles tend to gather at the bottom of the pipe, and the distribution of particles is non-uniform in the vertical direction. However, the flow of Case 3 is homogeneous with all particles suspending nearly uniformly. This is noted by both numerical results and

Please cite this article in press as: Jin T, et al., Numerical prediction of flow characteristics of slush hydrogen in a horizontal pipe, International Journal of Hydrogen Energy (2016), http://dx.doi.org/10.1016/j.ijhydene.2016.09.054

6

0.5 0.4

C4 [Expt][29] C4 [CFD] C5 [Expt][29] C5 [CFD]

0.3 0.2

y/D

0.1 0.0 -0.1 -0.2 -0.3 -0.4 -0.5 0.0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

V/Vmean

Pressure drop per unit length (kPa/m)

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y x x x ( 2 0 1 6 ) 1 e1 2

10 9

C6 [Expt][29]

8

C6 [CFD]

7 6 5 4 3 2 1 0 0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

4.5

Inlet velocity (m/s)

(a) normalized solid velocity

(b) pressure drop

Fig. 4 e Comparison between numerical data and experimental results of Cases 4e6.

Pressure drop characteristics Fig. 5 presents the numerical results of the pressure drop per unit length of the slush hydrogen and subcooled liquid hydrogen at the triple point. The calculation ranges are 1e4 m/ s for inlet flow velocity and 5e35% for solid fraction, and the

300

Pressure Drop (Pa/m)

experimental data. Nevertheless, there exists the data misfit at high velocity, i.e., the solid fraction gradient of numerical results is greater than that of experimental data, which indicates that the error of the model is considerable for predicting slurry flow at high flow rate or high solid concentration. The empirical parameters, namely, specularity coefficient, effective viscosity parameter, particleeparticle and particlewall restitution coefficients, adopted for slush nitrogen in Cases 4e5 are f ¼ 0:01, b ¼ 1:5, ess ¼ 0:9 and ew ¼ 0:95, respectively. Fig. 4 presents the comparison of numerical prediction and experimental results for Cases 4 and 5 of solid velocity distribution profiles at the inlet velocity of 1.68 m/s and 4.1 m/s, respectively. Here, V/Vmean is the normalized velocity at the cross-section of steady state flow. As shown in Fig. 4, the numerical results are in good agreement with the experimental results for both cases. For slush nitrogen, the effect of gravity on the flow is not obvious since the density ratio between solid and liquid phases is 1.18, resulting in the symmetrical distribution of solid velocity. However, for the low inlet velocity Case 4, the solid velocity values near the upper wall deviate about 8.9% from the experimental results. One of the possible reasons for the deviation is the deficit of lift modeling since the lift force is slight at low velocity, which needs further attention for model improvements. For Case 6, the numerical calculation data and experimental results of slush nitrogen flow at an initial velocity range of 1e4 m/s with solid fraction of 11e15 wt% in the horizontal circular pipe is presented in Fig. 4(b). The pressure drop per unit length is easily found to rise with the increased flow velocity for both the experimental and CFD results and the pressure drop values calculated by the 3D model agree fairly well with the experimental data for all the velocities considered in Case 6. Thus, the CFD model is considered effective for predicting the pressure drop of cryogenic slurries flow in horizontal circular pipe when the flow velocity is not too high.

Subcooled LH2 5 [%] 13 [%] 25 [%] 35 [%]

250 200 150 100 50 0

1

2

3

4

Inlet Velocity (m/s) Fig. 5 e Variation in pressure drop of slush hydrogen and subcooled liquid hydrogen for different inlet velocity and solid fraction.

particle diameter is 1 mm unless otherwise specified. The pressure drop data are obtained with the pipe length of 2e3 m, where the flow is regarded completely steady. Fig. 5 shows that the pressure drop for slush hydrogen rises with the increased flow velocity and solid fraction and that the increasing tendency is more obvious at higher velocity and concentration, however, the pressure drop of subcooled liquid hydrogen is greater than that of slush hydrogen at some cases of high velocity and low concentration. To better understand the pressure drop reduction results, the flow friction factors f of slush hydrogen and subcooled liquid hydrogen for different Reynolds numbers are given in Fig. 6. The dashed line denotes the result of the pffiffiffi pffiffiffi PrandtleKarman equation, i.e., 1= f ¼ 2:0  lgðRe f Þ  0:8. The Reynolds number is defined as Rel ¼ rl vD=ml , which is based on liquid flow. It can be seen that the friction factor of slush hydrogen is greater than that of liquid hydrogen at all Reynolds numbers with the solid fraction over 25%. On the other hand, for Rel exceeding 2.0  105 with solid fraction of

Please cite this article in press as: Jin T, et al., Numerical prediction of flow characteristics of slush hydrogen in a horizontal pipe, International Journal of Hydrogen Energy (2016), http://dx.doi.org/10.1016/j.ijhydene.2016.09.054

7

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y x x x ( 2 0 1 6 ) 1 e1 2

Friction Factor f

0.04 Prandtl-Karman Eq. Subcooled LH2 5 [%] 13 [%] 25 [%] 35 [%]

0.03

0.02

0.01

0.00

2.0E5

4.0E5

6.0E5

8.0E5

1.0E6

Renolds Number Rel Fig. 6 e Friction factor of slush hydrogen and subcooled liquid hydrogen for different Reynolds number.

rate for slush hydrogen is 45.5% under the condition with the inlet velocity of 4 m/s and the solid fraction of 5%. The pressure drop reduction phenomena may be a result of the gathering of solid particles in the center section of the pipe for slurry flow with high velocity and low concentration, accounting that solid particles in the pipe center region probably suppress the liquid turbulent flows near the pipe wall with decreasing turbulence kinetic energy. The influence of particle size on the slush hydrogen flow characteristics is presented in Fig. 8. The pressure drop of slush hydrogen for the particle diameters ds ¼ 0.5, 1, 2 mm with the inlet velocity of 1e4 m/s and the solid fraction of 13% are shown here and the pressure drops during tests are obtained at the pipe length of 2e3 m. As indicated in Fig. 8, the pressure drop for slush hydrogen rises with the increased particle size, and for the coarser particle, the rising rate is greater. Thus, the particle size is an important impact factor for the slush hydrogen flow. In order to manage the slush hydrogen flow pressure drop, the particle size should be controlled to be fine during the slurry production process.

0.04

Friction Factor f

0.03

0.02

0.01

0.00 0.0

2.0E5

4.0E5

6.0E5

8.0E5

1.0E6

Renolds Number Resl Fig. 7 e Friction factor of slush hydrogen and subcooled liquid hydrogen for different slush Reynolds number.

13% and for all Rel with the solid fraction of 5%, the friction factor of the liquid fluid is greater than that of the slush flow. Fig. 7 also presents the flow friction factors of slush hydrogen and subcooled liquid hydrogen with different Reynolds numbers. Nevertheless, the Reynolds number for the calculated case in Fig. 7 is the slush Reynolds number Resl ¼ rsl vD=mm , where rsl is the density of the slush hydrogen determined from rsl ¼ al rl þ as rs and mm is the effective viscosity of mixture. For the subcooled liquid hydrogen, Resl ¼ Rel . Apparently, the slush Reynolds number is more capable of characterizing the slurry flow. From Fig. 7, the friction factor for subcooled liquid hydrogen is greater than that of slush hydrogen for all concentrations with the Reynolds numbers of 1.35  105e5.0  105 and the difference between the friction factors for slush and liquid hydrogens rises with the decreased solid fraction. At the solid fraction over 13%, the friction factor for slush nitrogen drops as the Reynolds number rises and approaches to the PrandtleKarman equation line at the slush Reynolds number of 5.0  105. In the present calculation, the maximum friction factor reduction

Solid fraction distribution Figs. 9e11 show the simulated local solid volumetric concentration profiles and contours at the outlet cross-section of slush hydrogen flow. Fig. 9 presents the concentration profiles along the vertical centerline at the outlet for different inlet flow velocities (1e4 m/s) with the solid concentration of 13% (light) or 35% (high). Fig. 10 gives the solid concentration profiles for different inlet solid fractions (5e35%) with the velocity of 1 m/s (low) or 4 m/s (high). Generally, for all the inlet concentrations, the particles tend to gather at the lower part of the pipe because of the higher density of solid phase and the near-wall particle concentrations are minor due to the particle-wall interference. The highest value of solid fraction are about 45% for flow at light concentration (13%) and the gathering highest concentration point moves down along the centerline with the decreasing inlet velocity because of the significant suspension

200

Pressure Drop (Pa/m)

Prandtl-Karman Eq. Subcooled LH2 5 [%] 13 [%] 25 [%] 35 [%]

d=0.5 mm d=1 mm d=2 mm

150

100

50

0

1

2

3

4

Inlet Velocity (m/s) Fig. 8 e Pressure drop of slush hydrogen for different particle sizes and inlet velocity with the solid fraction of 13%.

Please cite this article in press as: Jin T, et al., Numerical prediction of flow characteristics of slush hydrogen in a horizontal pipe, International Journal of Hydrogen Energy (2016), http://dx.doi.org/10.1016/j.ijhydene.2016.09.054

8

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y x x x ( 2 0 1 6 ) 1 e1 2

0.5

0.5 U_in=1 m/s U_in=2 m/s U_in=3 m/s U_in=4 m/s

0.4 0.3 0.2 0.1

0.3 0.2 0.1

y/D

y/D

0.0

0.4

-0.1

U_in=1 m/s U_in=2 m/s U_in=3 m/s U_in=4 m/s

0.0 -0.1

-0.2

-0.2

-0.3

-0.3

-0.4

-0.4 -0.5

-0.5 0.0

0.1

0.2

0.3

0.4

0.0

0.5

0.1

0.2

0.3

0.4

0.5

0.6

Solid volumetric concentration

Solid volumetric concentration

(b) inlet solid fraction of 35%

(a) inlet solid fraction of 13%

Fig. 9 e Volumetric concentration along vertical centerline for different inlet velocity.

0.5 0.3 0.1

0.3 0.2 0.1

y/D

0.0 -0.1

0.0 -0.1

-0.2

-0.2

-0.3

-0.3

-0.4

-0.4

-0.5

-0.5

0.0

0.1

0.2

0.3

0.4

0.5

5[%] 13[%] 25[%] 35[%]

0.4

25[%] 35[%]

0.2

y/D

0.5

5[%] 13[%]

0.4

0.6

Solid volumetric concentration (a) Inlet flow velocity of 1m/s

0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

Solid volumetric concentration (b) Inlet flow velocity of 4 m/s

Fig. 10 e Volumetric concentration along vertical centerline at outlet cross-section for different solid fraction.

Fig. 11 e Volumetric concentration in vertical plane at outlet for different solid fractions with inlet flow velocity of 1 m/s.

force of solid phase at high velocity. However, the particle gathering exhibits a broader range and a greater maximum concentration value of 55% for flow at higher concentration (35%). On the other hand, at low velocity (1 m/s), the solid fraction distribution curves are more obviously affected by various inlet solid concentrations and the particles deposit and

gather more heavily at the lower part of the pipe. That is moving bed flow, where the particles tend to settle at the bottom of the pipe move by sliding along the bottom line. Also, the influence of velocity on the solid fraction distribution is more significant at high concentration, which is apparently shown in Fig. 11 by the means of solid fraction distribution contours.

Please cite this article in press as: Jin T, et al., Numerical prediction of flow characteristics of slush hydrogen in a horizontal pipe, International Journal of Hydrogen Energy (2016), http://dx.doi.org/10.1016/j.ijhydene.2016.09.054

9

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y x x x ( 2 0 1 6 ) 1 e1 2

Fig. 12 e Volumetric concentration in cross-sections along flow direction (solid fraction of 13%, inlet flow velocity of 4 m/s).

To better understand the solid fraction distribution in the flow development process, the local volumetric concentration in the typical cross-section along the flow direction is given in Fig. 12. The inlet solid fraction for the numerical case is 13% and the inlet flow velocity is 4 m/s. At the early flowdeveloping part, the particles are almost homogenously distributed in the cross-section. As the flow proceeds, it can be regarded fully developed at x ¼ 2.5 m, most particles gather in the lower central scope of the pipe and the flow is considered pseudo-homogeneous. The fact that the flow pressure drop reduction phenomenon exists for the cases with high velocity and low particle concentration can be attributed to the significant interreaction among the particles and between the particles and the liquid for pseudo-homogenous flow.

Velocity profiles Figs. 13 and 14 present the numerical solid phase velocity profiles and contours at the outlet cross-section of slush hydrogen flow. Fig. 13 shows the quantitative particle velocity distribution along the vertical centerline at the outlet for different inlet flow velocities (1e4 m/s) with a light concentration of 13% and a high concentration of 35%. For high concentration (35%), the particle velocity distribution exhibits some extent of symmetry and becomes more symmetrical as flow velocity is increased. Also, with the increased velocity, the particles in the lower half of the cross-section are accelerated with the corresponding de-acceleration in the upper half. These results indicate that the flow becomes pseudo-homogeneous

0.5

0.5

0.4

0.4

0.3 0.2

y/D

0.1 0.0

y/D

U_in=1 m/s U_in=2 m/s U_in=3 m/s U_in=4 m/s

-0.1

0.3

U_in=1 m/s

0.2

U_in=2 m/s

0.1

U_in=3 m/s

0.0

U_in=4 m/s

-0.1

-0.2

-0.2

-0.3

-0.3

-0.4

-0.4

-0.5 0.0

-0.5 0.2

0.4

0.6

0.8

1.0

V/V_mean (a) solid fraction of 13%

1.2

1.4

1.6

0.0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

1.6

V/V_mean (b) solid fraction of 35%

Fig. 13 e Solid velocity profiles along vertical centerline for different inlet flow velocity with solid fraction. Please cite this article in press as: Jin T, et al., Numerical prediction of flow characteristics of slush hydrogen in a horizontal pipe, International Journal of Hydrogen Energy (2016), http://dx.doi.org/10.1016/j.ijhydene.2016.09.054

10

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y x x x ( 2 0 1 6 ) 1 e1 2

0.5

0.5

0.4

0.4 5[%] 13[%] 25[%] 35[%]

0.3 0.2 0.0

5[%]

0.2

13[%] 25[%] 35[%]

0.1

y/D

y/D

0.1

0.3

-0.1

0.0 -0.1

-0.2

-0.2

-0.3

-0.3

-0.4

-0.4 -0.5

-0.5 0.0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

1.6

0.0

0.2

V/V_mean

0.4

0.6

0.8

1.0

1.2

1.4

V/V_mean

(a) flow velocity of 1 m/s

(b) flow velocity of 4 m/s

Fig. 14 e Solid velocity profiles along vertical centerline for different solid fraction with flow velocity.

Fig. 15 e Solid velocity profiles in vertical plane at outlet for solid fraction of 13% and 35% with inlet flow velocity of 1 m/s.

when the inlet velocity is high. Fig. 14 shows the quantitative particle velocity distribution for different inlet solid fractions (5e35%) with a low velocity (1 m/s) and a high velocity (4 m/s). At the low flow velocity (1 m/s), the solid particles are blocked and concentrated in the lower half of the pipe. The velocity of solid phase in the lower half of the cross-section is reduced at higher concentration. For the flow with high solid concentration and low velocity, it can be heterogeneous flow and the moving bed flow as shown in the velocity distribution contours of Fig. 15, which indicates that the particles deposit in sufficient quantity or move slowly at the bottom of the pipe. The solid velocity distribution contours at the typical cross-sections along the flow direction for an inlet solid fraction of 13% and a flow velocity of 4 m/s are given in Fig. 16. At the early flow-developing part, the particle velocity gradient throughout the cross-section is slight and as the flow proceeds, the velocity differences at the top and bottom parts become apparent, and the particles gradually concentrate and interreact with each other high-frequently. It can be noted from Fig. 16 that the particles at the lower part of the pipe move at high velocity which can avoid the deposition at the bottom, indicating that the flow is pseudo-homogenous flow.

Conclusion A 3-D numerical code based on EulereEuler two-fluid model was developed to study the flow characteristics of slush hydrogen in a horizontal circular pipe, where the effective viscosity of slush mixture, which takes the particle shape and size into consideration, is introduced to the model. The effective viscosity can help to modify the correlations for interphase momentum exchange. The other improvement in the wall boundary conditions for the solid phase are based on the JohnsoneJackson model [15] which involved the friction and collision between the particle and the wall. The specularity coefficient, effective viscosity parameter, particleeparticle and particle-wall restitution coefficients for the granular phase of slush hydrogen are f ¼ 0:01, b ¼ 1:0, ess ¼ 0:9 and ew ¼ 0:7, respectively. After the comprehensive validation with respect to the experimental data in the literatures, the numerical model was verified to be capable of predicting the flow characteristics of slush hydrogen when the flow velocity is not too high. Extensive cases at various working conditions (inlet solid fraction of 5e35%, inlet velocity of 1e4 m/s, particle size of 0.5e2 mm) were calculated to analyze the pressure

Please cite this article in press as: Jin T, et al., Numerical prediction of flow characteristics of slush hydrogen in a horizontal pipe, International Journal of Hydrogen Energy (2016), http://dx.doi.org/10.1016/j.ijhydene.2016.09.054

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y x x x ( 2 0 1 6 ) 1 e1 2

11

Fig. 16 e Solid velocity profiles on cross-sections at flow direction (solid fraction of 13%, inlet flow velocity of 4 m/s).

drop, solid concentration distribution and velocity profiles. The main conclusions that can be derived are: 1) The pressure drop and flow friction factor for slush hydrogen generally rise with the increased flow velocity and the rising tendency is more obvious at high velocity and concentration. The pressure drop reduction is confirmed in the numerical analysis since the friction factor for subcooled liquid hydrogen is greater than that of slush hydrogen in the cases for all the concentrations at slush Reynolds numbers of 1.35  105e5.0  105. At a solid fraction over 13%, the friction factor for slush nitrogen drops as the Reynolds number rises and approaches the PrandtleKarman Eq. line at slush Reynolds number of 5.0  105. The influence of particle size (ds ¼ 0:5; 1; 2 mm) on the slush hydrogen flow characteristics is elucidated. The pressure drop for slush hydrogen rises with the increased particle size and the rising rate is greater for the coarser particle. Thus, the particle size should be managed to be fine in the slurry production process. 2) Particles tend to gather at the lower parts of the pipe because of the higher density of solid phase, and the nearwall particle concentration is minor due to the particlewall interference. At a higher concentration, particle gathering exhibits a broader range, and the particle velocity distribution exhibits to be more symmetrical. Also, with the increased velocity, the particles in the lower half of the cross-section are accelerated while with the corresponding de-acceleration in the upper half, which results in the pseudo-homogeneous flow. 3) The pressure drop reduction phenomenon (i.e., the pressure drop for subcooled liquid flow is higher than that of slush nitrogen) is not yet fully understood. Detailed analysis about the basic fluid mechanism is

necessary for promoting the research on the flow characteristics of cryogenic slush hydrogen. Also, further experimental and numerical investigations should concern the influence of heat transfer, pipeline type and geometric structure parameters on the flow characteristics of slush hydrogen.

Acknowledgment This work is financially supported by the Zhejiang Provincial Natural Sciences Foundation (LZ14E060001) and the State Key Laboratory of Technologies in Space Cryogenic Propellants (SKLTSCP1408).

Nomenclature CD CL CVM D ds ess ew ! F D;q ! F L;q ! F VM;q Gk ! g g0;ss k kQs L

Drag coefficient Lift coefficient Virtual mass coefficient Pipe diameter, mm Particle diameter, mm Restitution coefficient for particles Restitution coefficient for particles and wall Drag force Lift force Virtual mass force Generation of turbulence kinetic energy Gravitational acceleration, m2/s Radial distribution function Turbulence kinetic energy Diffusion coefficient Pipe length, m

Please cite this article in press as: Jin T, et al., Numerical prediction of flow characteristics of slush hydrogen in a horizontal pipe, International Journal of Hydrogen Energy (2016), http://dx.doi.org/10.1016/j.ijhydene.2016.09.054

12 m_ p ps ! q Re Rep Rew ! U sw ! v

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y x x x ( 2 0 1 6 ) 1 e1 2

Mass flow, kg/s Pressure, Pa Solids pressure, Pa Heat flux, W/m2 Reynolds number Particle Reynolds number Vorticity Reynolds number Particle slip velocity parallel to the wall, m/s Velocity, m/s

Greek letters a Volume fraction, % The packing limit, % as;max b Effective viscosity parameter ε Rate of dissipation f Specularity coefficient Interphase energy exchange fls Collisional dissipation of energy gQs ½h Intrinsic viscosity, Pa$s m Fluid viscosity, Pa$s Turbulent interaction source term Pk ,Pε Granular temperature, m2/s2 Qs r Density, kg/m3 t Shear stress Subscripts l Liquid phase s Solid phase sl Slush p Either solid phase or liquid phase q The opposite phase of p

references

[1] Carney RR. “Slush hydrogen”production and handling as a fuel for space projects. In: Advances in cryogenic engineering. Springer; 1964. p. 529e36. [2] Smith EM. Slush hydrogen for aerospace applications. Int J Hydrogen Energy 1989;14(3):201e13.  lu TN, Sheffield JW. Review of [3] Gu¨rsu S, Sheriff SA, Vezirocg slush hydrogen production and utilization technologies. Int J Hydrogen Energy 1994;19(6):491e6. [4] Park YM. Literature research on the production, loading, flow, and heat transfer of slush hydrogen. Int J Hydrogen Energy 2010;35(23):12993e3003. [5] Hardy TL. FLUSH: a tool for the design of slush hydrogen flow systems. NASA-TM-102467, Lewis Research Center; 1990. [6] Hardy TL, Whalen MV. Slush hydrogen propellant production, transfer, and expulsion studies at the NASA KSite Facility. Cleveland, OH: Aiaa Journal; 1991. p. 1e7. [7] Takakoshi T, Murakami M, Ikeuchi M, Matsuo K, Tsukahara R. PIV measurement of slush nitrogen flow in a pipe. In: Proceedings of the Advances in Cryogenic Engineering: Transactions of the Cryogenic Engineering Conference-CEC. AIP Publishing; 2006. p. 1025e32. [8] Ohira K. Pressure drop reduction phenomenon of slush nitrogen flow in a horizontal pipe. Cryogenics 2011;51(7):389e96.

[9] Zhang P, Jiang YY. Forced convective heat transfer of slush nitrogen in a horizontal pipe. Int J Heat Mass Transf 2014;71(1):158e71. [10] Gamma F, del Monte L, Liberatore R. 1998. CFD-analysis of slush hydrogen flows in feed-systems of rocket-based combined cycles. AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit, Cleveland, OH. [11] Ishimoto J, Ono R. Numerical study of the two-phase flow characteristics of slush nitrogen. Cryogenics 2005;45(4):304e16. [12] Ohira K, Ota A, Mukai Y, Hosono T. Numerical study of flow and heat-transfer characteristics of cryogenic slush fluid in a horizontal circular pipe (SLUSH-3D). Cryogenics 2012;52(7):428e40. [13] Jiang YY, Zhang P. Numerical investigation of slush nitrogen flow in a horizontal pipe. Chem Eng Sci 2012;73:169e80. [14] Zhang P, Shi X. Numerical investigation of slush hydrogen flow in horizontal pipes. CIESC J 2014;65(z2):38e44. [15] Johnson PC, Jackson R. Frictional-collisional constitutive relations for granular materials, with application to plane shearing. J Fluid Mech 1987;176(1):67e93. [16] Leachman JW, Jacobsen RT, Penoncello SG, Lemmon EW. Fundamental equations of state for parahydrogen, normal hydrogen, and orthohydrogen. J Phys Chem Reference Data 2009;38(3):721e48. [17] Saffman PG. The lift on a sphere in a slow shear flow. J Fluid Mech 1965;22(2):385e400. [18] Roco MC. Particulate two-phase flow. Boston, U.S: Butterworth-Heinemann; 1993. [19] Lu H, Gidaspow D. Hydrodynamics of binary fluidization in a riser: CFD simulation using two granular temperatures. Chem Eng Sci 2003;58(16):3777e92. [20] Wen CY, Yu YH. Mechanics of fluidization. Chemengng Progsympser 1966;62:100e11. [21] Ergun S. Fluid flow through packed column. J Mater Sci Chem Eng 1952;48(2):89e94. [22] Schiller L, Naumann Z. A drag coefficient correlation. Vdi Ztg 1935;77(318):51. [23] Einstein A. Zur theorie der brownschen bewegung. Ann Der Phys 1906;324(2):371e81. [24] Mooney M. The viscosity of a concentrated suspension of spherical particles. J Colloid Sci 1951;6(2):162e70. [25] Thomas CU, Muthukumar M. Three-body hydrodynamic effects on viscosity of suspensions of spheres. J Chem Phys 1991;94(7):5180e9. [26] Mori Y, Oshio Y, Shintou H. Fundamental technology of the pipeline engineering. Tokyo: Keigaku Syuppan; 1980. [27] Barnes HA, Hutton JF, Walters K. An introduction to rheology. Amsterdam, Netherland: Elsevier; 1989. [28] Cheng NS, Law AWK. Exponential formula for computing effective viscosity. Powder Technol 2003;129(1e3):156e60. [29] Gidaspow D, Bezburuah R, Ding J, et al. 1992. Hydrodynamics of circulating fluidized beds: kinetic theory approach; Proceedings of the International Conference on Fluidization, Gold Coast, Australia. [30] Lun CKK, Savage SB, Jeffrey DJ, Chepurniy N. Kinetic theories for granular flow: inelastic particles in Couette flow and slightly inelastic particles in a general flowfield. J Fluid Mech 1984;140:223e56. [31] Roco MC, Shook CA. Modeling slurry flow: the effect of particle size. Can J Chem Eng 1983;61(4):494e503.

Please cite this article in press as: Jin T, et al., Numerical prediction of flow characteristics of slush hydrogen in a horizontal pipe, International Journal of Hydrogen Energy (2016), http://dx.doi.org/10.1016/j.ijhydene.2016.09.054