Kyushu University Global COE Program
Journal of Novel Carbon Resource Sciences, Vol. 3, pp. 11-16, Feb. 2011

Estimation of the Multi-phase Flow on the Vertical Pipe for the
Methane Hydrate Recovery
Shuji Hironaka *1, Yuki Fujisawa *2, Saki Manabe *2,
Gen Inoue*3, Yosuke Matsukuma*3, Masaki Minemoto*3*1

Graduate School of Engineering, Kyushu University*2
School of Engineering, Kyushu University*3
Faculty of Engineering, Kyushu University(Received December 13, 2010; accepted January 11, 2011)
In this study, the experiments with two vertical pipes of different diameters were conducted to obtain therelationship between the flow rate of the gas and liquid phase. The simulation model of three-phase flow was
proposed by the experimental result, which are the drag force between two-phase and the wall friction. The
simulation results are good agreement in the experimental results and it was clarified that the proposed simulation
model is useful to estimate the three-phase flow in the vertical pipe of the methane hydrate recovery system.

1. Introduction

Recent survey reveals there is methane gas estimatedmore than 100 times amount of domestic consumption in2007 on the seabed near Japan in the form of hydrate1).
It is methane hydrate (MH), which is stable crystal
like an ice under the condition of low temperature andhigh pressure and formed by water molecules’ trappinga methane molecule. The MH is expected as one of
new energy resource2) for Japan. The present platform
designed for oil or natural gas is, however, not applicablefor MH recovery because they are designed for fluid
resources, not for solid. Then, several kinds of systemhave been studied to recover the MH such as low pressure method3), heat injection method, inhibitor injection method, decomposition method, CO2 hydrate replacing method4) and gas lift method. Each of them has

Upward flow

both of advantages and disadvantages. For example, theheat injection method has disadvantages that heat-lossof hot water would be a serious problem though it has
been examined and most come into practical use at theMackenzie, Canada5). In the present paper, we focusedon the gas-lift method because it could have the advantage of lower running cost over the other methods6). It
is because of self-gas-lift effect of MH. It says that thedecomposition of MH might have additionally effects assame as injection of the methane gas, which is involved
by high temperature seawater and lower pressure as thehydrate is raised7, 8). Therefore less gas injection allows to
run in smaller cost. Enough knowledge is, however, notobtained for the multi-phase flow on a large scale vertical pipe, and it prevents us from constructing the actualcommercial plant. Figure 1 shows the conceptual imageof gas-lift system and rough principle. The up-stream
flow appears on the riser pipe by injection of methane
gas in the middle of riser pipe. Taking the force balance
between inside and outside of vertical pipe for this system, the following equation is obtained.
(Pressure at the bottom) - (Pressure at the Top) + (Pressure drop by wall friction)= (Hydrostatic pressure difference between top and bottom)(1)
Therefore,

Gasification of themethane hydrate
CH4 gasinjection(1－ε ) ρl gH + 4 f

1H
ρ U 2 = γρl gH
2D l l(2)

and we can estimate the superficial velocity of water by
the eq. (3), assuming other variables is constant.

Brokenhydrate lumps

Ul =(ε + γ －1) gD

Drill
Methane hydrate layer

Fig. 1 Conceptual image of gas-lift system.

2f

(3)

Here, γ is submerge ratio given by following equation.
γ = Ls/H(4)

Journal of Novel Carbon Resource Sciences, Vol.3, Feb. 2011

Figure 2 shows schematic of an experimental apparatus. The riser pipe (I) is made of acryl, downcomer (II) isof vinyl chloride (PVC), and both of them are connectedwith two acryl tanks at the top and bottom. Experiments
were conducted with two kind of diameter of riser pipe:
50 mm and 230 mm, and the height were 5.5 m and 5.73
m. These apparatus have the extra pipe (III) to measure
the static surface level for the submerge ratio. The air
and water were used as gas and liquid phase. The orificelike injection hole were set at 4.4 m deep for 50 mm riser
pipe, and 4.63 m for 230 mm from top of riser pipe. Air
is injected by the (IV) compressor through the (V) dryer
and it makes up-stream two-phase flow on the riser pipe.
Flow rate of air and water were recorded by mass flow
meter (VI) and electromagnetic flow meter (VII) respectively.

(3)

(2)

Extra pipe Downcommer

Ls

H(1) Riser pipe
(4) Compressor

Ul [m/sec ]2. Experiment

The volume flow rate of water and that of injected gas
were measured to obtain the superficial velocity.
According to eq. (3), water superficial velocity will
depend on the submerge ratio. Therefore, we paid
attention to keep it almost same in recording thesuperficial velocity by adjusting the amount of water on
the apparatus.
Figures 3 and 4 show the dependence of the watersuperficial velocity on the gas superficial velocity
when the submerge ratio was changed from 0.95 to
0.99 at 0.01 intervals for 50 mm and 230 mm diameter
pipe respectively. They also show the water superficial
velocity is influenced by the submerge ratio, because of
the larger hydrostatic pressure.

0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0

◆ γ= 0.95▲ γ= 0. 96■ γ= 0.97× γ= 0.98＊ γ= 0.99

0.00
Fig. 3Ul [m/sec ]It is the ratio of the length of recovery pipe under thesea surface to its total length. However, eq. (3) will not
be accurate as the amount of injected gas becomes larger
because the drag force between gas and liquid phases
is not negligible. This drag force is different from thesystem and various experiments and simulations9) havebeen conducted. On the other hand, from the view toapply the gas-lift system, large number of studies on
the multi-phase flow have been conducted10). They are,however, usually assumed that the phase transition will
not occur. In the paper, we will show results of experimentwith two vertical riser pipes of different diameters in thelaboratory scale to measure the relationship between
superficial velocity of the injection air and water. Next, thenumerical simulation was conducted to the experimental
system with one-dimensional multi-fluid model, andcompared with experimental result. And finally, theone-dimensional simulation of the three-phase flow was
carried out in the actual scale with the proper model
equations, including the heat transfer between hydrateand seawater, and hydrate’s decomposition.

0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0

0.02

(6)

0.08

0.10

Ul ’s dependence on Ug (50 mm Diameter).

◆ γ= 0.95▲ γ= 0. 96■ γ= 0.97× γ= 0.98＊ γ= 0.99

0.00

0.02

(5) Dryer

M Mass flowmeter

0.04
0.06
U g [m/sec ]Fig. 40.04
0.06
U g [m/sec ]

0.08

0.10

Ul ’s dependence on Ug (230 mm Diameter).

(7) Electromagnetic
flowmeter

3. Numerical simulation
3.1 Assumptions
Fig. 2 Schematic of experimental apparatus.

Basic equations for numerical simulation were derived
based on the multi-fluid model, which is to solve thepressure, volumetric fractions, and velocity of each

Estimation of the Multi-phase Flow on the Vertical Pipe for the Methane Hydrate Recovery

phase. The following assumptions were set to conduct
the calculation:I)

The multi-phase flow is assumed as the one-dimensional flow on the riser pipe.
II) Riser pipe does not have any flexures, thus its
cross sectional area is uniform.
III) The bubble form is spherical in the bubble flow
region.
IV) The contracted coefficient at the bottom of the riser pipe was 1.0.
V) All phases exist under the same pressure.
VI) Form of solid phase is spherical.
VII) Only gas density changes depending on the pressure and temperature, and other phase densities
are constant.
VIII) For the temperature of each phase, the gas andliquid's is equal, and the solid temperature isequilibrium temperature of the methane hydrate
under the surrounding pressure. Its distribution
through the solid sphere is uniform.
IX) The heat transfer is only considered for the liquid
phase.
X) The decomposed hydrate or ice is NOT reformed
again.

3.2 Equations

Concrete equations are given by following eqs. (5) ~
(8).

∂
∂
(α k ρk ) + (α k ρk vk ) = K
∂t
∂x

(5)

∂
∂
∂p
−α k
− ( fric, k ) vk
(α k ρk vk ) + (α k ρk vk2 ) =∂t∂x∂x− ∑ ( inter )kk ′ ( vk − vk ′ ) − ρ kα k g − Kvk(6)k ′≠ k∑α = 1

(7)

k

∂T 
 ∂Tl
+ vl l  =Qwall + Qls
∂t 
 ∂t

ρ lα l C p ()

Several terms appearing basic equations are semiexperimentally given by the following constitutive
equations. For wall friction term, it is given by following
eq. (12).

( fric, k ) =
Here,

2mf m
wk ρ kα k vkD

4.8 + 0.4 ln α g
m=
4.2

α g < 0.2
α g ≥ 0.2

16
 Rem
fm = 
0.0791 Re0.25mRem < 2300

(12)

(13)

(14)

Rem ≥ 2300

α g ≤ 0.9
0wg = 
2
(10α g − 9 ) ( 21α g − 20 ) α g ≥ 0.9(15)

wl = 1(16)

ws = 0(17)

Eq. (14) is the Blasius' equation, and it is only applicableto the single phase. Therefore, it is required to define theReynolds number to the mixture fluid to apply the gasliquid two-phase flow.

Rem =α g ρ g vg + α l ρl vl D

(α µ + α µ )
g

g

l

(18)

l

The drag force between gas and liquid phase is definedby eq. (19).

18( inter )lg = Cd ,l a ρl vg − vl

(19)

(8)

Eq. (5) is the equation of mass conservation, eq. (6)represents the momentum conservation, eq. (7) is thedefinition of volumetric fraction, and eq. (8) shows theheat conservation on the liquid phase. K denotes theamount of transition of the each phase, and they are given
by following equations respectively.

6α h (T − T )S=− s ⋅ ls l s
d s H decompθ p

3. 3 Constitutive Equations(9)L=− 1 − σ CH 4 θ p S(10)

G = −σ CH 4θ p S(11)

0.54α g−0.482Cd , l =  8  1 − α g   3  1 − α gs 

0 < α g ≤ 0.75

(20)0.75 < α g ≤ 1

 6α g ( 0.3 + 0.2α g )0 < α g ≤ 0.25

dg

a=
 4.5 α g − α gs 6α gs 1 − α g D ⋅ 1 − α + d ⋅ 1 − α 0.25 < α g ≤ 1
gs
g
gs


(21)

Moreover, the averaged gas volume fraction for theslug flow and the bubble diameter is given by following
equations.

Journal of Novel Carbon Resource Sciences, Vol.3, Feb. 2011( 0 ≤ α < 0.25)() (22)
(
)

g
α g

α gs 0.3929 − 0.5714α g 0.25 ≤ α g < 0.715
=

0.05
0.715 ≤ α g < 1(23)

d g = 0.95 DU g0.2

Similarly, the drag force between liquid and solid
phase is also defined by eq. (24).

( inter )ls = 6α 
1
Cd , s  s  ρl vl − vs8 ds 

(24)

calculation and experiment with γ changed for two
different diameters, 50 mm and 230 mm. Since theexperiment was conducted in the gas-liquid two-pahse
system, the solid phase is neglected. The experimental
results are dotted and simulations are lined. Simulation
results fairly agreed with the exepriment in the case ofdiameter D = 50 mm. Although the calculation results
underestimate liquid superficial velocity in the case of
D = 230 mm especially at the high gas velocity region, it
shows the similar qualitative tendency to the experiment.
In order to consider of the reason why the simulation
underestimated the liquid superficial velocity, the ratio
between the interface drag force and wall friction istaken with eq. (12) and (19).

( fric, l ) =( inter )And the drag force coefficient between liquid and solid,
and the solid diameter in eq. (24) are given by following
eq. (25) and (26).

24
 Re p
Cd , s = 10 0.5
 Re p0.44


2 < Rep < 500

(25)

0.9
0.8

500 < Rep

0.7

ρl vs − vl d s.
µl

The density of gas phase is given by the equation of
the state for the ideal gas because of the compressibility.

M CH 4 p(27)
RTl
Additionally, the temperature of the solid phase is sameto the equilibrium temperature of the methane hydrate asa function of the pressure.
Ts =8533.8
38.98 − ln  p5−3  (1.013 × 10 × 9.869 × 10 ) 

(28)

(31)

2

4U (Tl − Tsur )D

T L 
 6α 
Qls = −  s  hls (Tl − Ts ) + l ρlα l  d s 

4. Result and discussion0.6

γ=0.99,Exp.
γ=0.97,Exp.
◆ γ=0.95,Exp.
γ=0.99,Calc.
γ=0.97,Calc.
γ=0.95,Calc.
●
▲

0.5
0.4
0.3
0.2
0.1
0.0
0.00

0.02

0.04

0.06

0.08

0.10

U g [m/sec ]Fig. 5

1.2

Ul ’s dependence on Ug (50 mm Diameter).

γ=0.99,Exp.
γ=0.97,Exp.
◆ γ=0.95,Exp.
γ=0.99,Calc.
γ=0.97,Calc.
γ=0.95,Calc.
●
▲

1.0

The following equations represent the term of the heatexchange between the liquid and wall or solid phase in
the eq. (8) respectively.
Qwall = −Ul [m/sec ]

(26)

(29)

0.8

Ul [m/sec]α  3ds = ds0  s 
 α s0 

ρg =

 v 
DCd ,l a 1 − g vl Rep < 21

Here, Re p =16 f mα l

0.6
0.4

(30)

The simulation was compared with the experimentshowed in the section 2 to evaluate the simulation
model. Figures 5 and 6 show the comparison between

0.2
0.0
0.00

0.02

0.04

0.06

0.08

0.10

U g [m/sec]

Fig. 6

Ul ’s dependence on Ug (230 mm Diameter).

In eq. (31), the ratio between gas and liquid velocity,interfacial area density, and drag force coefficient
become more effective as the submerge ratio is smaller,
and the amount of the injection air is larger despite thatwall friction coefficient and liquid volume fraction arealmost constant. Therefore, eq. (31) shows that the interphase drag force is more dominant under the condition of
the smaller submerge ratio, larger air injection, and larger
pipe diameter. In conclusion, the present calculation
model is necessary to be improved if it is necessary to
calculate under the condition of the small submerge ratio.
Since it is obvious that the effect of the inter-phasedrag force is negligible to the present calculation model
under the condition of the larger submerge ratio around
0.99, the estimation of the actual system is conductedwith the submerge ratio set in 1.00. The effect of hydrate
decomposition was considered in this simulation. The
calculation condition is shown in Table 1.

in the previous figure, the gas expansion also involves
the increase of the velocity of the each phase to satisfy
the mass balance on the arbitrary control volume. Finally, Figure. 9 shows the decomposition ratio of the solid
phase on the left hand and the solid diameter on the right.
The decomposition rate is large in the case to recover thesolid of 0.01 m diameter compared with the result of the0.05 m diameter. It corresponds to the large mass flowrate of the solid phase in the case of small diameter as it
was previously denoted. Additionally, it shows the largesolid decomposes more slowly than the case of the small
diameter. Results in both cases suggest that the gas-liftsystem will be economical to run because the complete
decomposition of the methane hydrate allows to cut down
the amount of the injection CH4 gas.

Table 1 Calculation condition
Item
Value
Unit
Diameter of riser pipe
0.47[m]
Length500.0[m]
Injection depth
250.0
[m]
Heat conductivity of the
15.0[W/(m . K)]
pipe wall
Thickness of the pipe
0.019
[m]
Gas superficial velocity
0.5[Nm/sec]
Submerge ratio
1.00[-]
Temperature of liquid
277.0
[K]
phase at the surface
Density of liquid phase1024.0[kg/m3]
Density of solid phase2100.0[kg/m3]
Diameter of solid at the inlet 0.01, 0.05[m]Volume fraction of the
5.0[vol.%]
hydrate at the inlet
Decomposition enthalpy of 450,000[J/kg]
hydratePorosity
0.3[m3-pore/m3-sol.]
0.05[m3-CH4/m3-pore]
Saturation ratio of the CH4vg , vl , vvs [m/sec]

Estimation of the Multi-phase Flow on the Vertical Pipe for the Methane Hydrate Recovery

12.0
10.0
vg ( 0.01m )vg ( 0.05m )vl ( 0.01m )vl ( 0.05m )vs ( 0.01m )vs ( 0.05m )

8.0
6.0
4.0
2.0
0.0
0.0

Seabed

100.0

200.0

300.0

x [m]

400.0

500.0

Surface

Distribution of the velocity of the each phase.

1.0
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0

αg ( 0.01m )αg ( 0.05m )αl ( 0.01m )αl ( 0.05m )αg , αal [−]Fig. 7

0.0

100.0

200.0

300.0

x [m]

400.0

500.0

2.0
1.8
1.6
1.4
1.2
1.0
0.8
0.6
0.4
0.2
0.0

0.06
0.05
0.04
0.03
0.02

ds

Decomposition ratio
S [kg/(m3 . sec)]

Fig. 8 Distribution of the volume fraction of the gas
and liquid phase.

Solid diameter [m]Figures 7, 8 and 9 represent the distribution of the velocity, volume fraction, decomposition rate and diameterof the solid phase respectively. The x axis is directed tothe sea surface with set the origin at the bottom of theriser pipe (seabed) in all these 3 figures. In Figure 7, thevelocity of the gas phase appears at the deeper region because of the decomposition of the hydrate. Since the gas
hydrate completely decomposes before reaching at thesea surface, the solid velocity disappears at the middle
point. Moreover, this figure also shows the fundamental
qualitative tendency that the small solid is easier to recover than the larger solid. Figure 8 is the distribution of
the volume fractions of the each phase in two cases forthe solid diameter at the inlet. The gas phase suddenly
expands especially around the sea surface because thepressure decreases. As a result, its volume fraction becomes large, and the liquid's fraction is small. As we saw

14.0

S ( 0.01m )S ( 0.05m )ds ( 0.01m )ds ( 0.05m )

0.01

0.0

Fig. 9

100.0

200.0

300.0

x [m]

400.0

0.00
500.0

Distribution of the decomposition ratio and solid
diameter of the solid phase.

Journal of Novel Carbon Resource Sciences, Vol.3, Feb. 2011

5. Conclusion

The gas-lift experiment is conducted in the laboratory scale, and the simulation model is proposed to apply the one-dimensional two-fluid model. The numerical simulation fairly agreed with the experiment for the
relationship between the gas and liquid flow rate in thelarge submerge ratio. However, it is not able to estimate
the liquid superficial velocity under the small submerge
ratio, and it is because of the inter-phase drag force
model. Then, the multi-phase flow on the actual systemis estimated with the proposed model, and the decomposition of the hydrate considered. As a simulation result,the methane hydrate will completely decompose and it
showed the gas-lift system can be an economical runningsystem to recover the methane hydrate.
References

M. Sato, J. the Japan Institute of Energy, 80, 635 (2002).
S. Y. Lee, G.D. Holder, Methane hydrates potential as
future energy sources fuel processing technology, 71,
181 (2001).
3) Y. Masuda, Aquabiology, 23, 478 (2001).
4) I. Aya, K. Yamane, T. Kojima, S. Namie, The highpressure science and technology for the oceanic CO2
isolation and CO2 hydrate, 12- 1, 40 (2002).
5) R. Hamaguchi et al., Kagaku Kogaku Ronbunshuu,
31(1), 68 (2005).
6) R. Gregor, Geochimica et Cosmochimica Acta., 64, 285(2004).
7) S. K. Kelkar, Hydrate dissociation rates in pipelines
Fluid Phase Equilibria, 371 (1998).
8) M. Ishii, K. Mishima, I. Kataoka, G. Kocamustafaogullari,
Proc 9th U.S. National Congress of Applied Mechanics,
73 (1982).
9) S. Michio, N. Kawahara, T. Kimura, J. Nakao, Kikai
Gakkai Ronbunshuu(B) 74-742, 1263 (2008).
10) T. Yoshinaga, Y. Sato, Int. J. Multihphase Flow, 22(2),223 (1996).
1)2)Nomenclature
Lettersa: Interfacial area density
between liquid and gas phase: Heat capacity per unit massCp
at constant pressure
Cd,l : Drag force coefficient
between liquid and gas phase
Cd,s : Drag force coefficient
between liquid and solid phase[1/m][J/(K . kg)]
[-]
[-]

D
: Pipe diameter[m]: Bubble diameter[m]dg: Solid diameter[m]ds: Solid diameter at the inlet
[m]
ds0
[-]
f, fm : Wall friction coefficientg: Gravitational acceleration[m/sec2](fric,k) : The wall friction of the each phase [kg/m3 . sec]
H
: Pipe length[m]
Hdecomp : Decomposition enthalpy of themethane hydrate[J/kg]: Heat transfer coefficient between
hls
liquid and solid[W/(m . K)](inter) : Inter-phase drag force[kg/m3 . sec]
K
: Amount of the decomposed/generated mass
for the each phase (K= G, L, S)[kg/m3 . sec]
: Static surface level[m]
Ls[kg/mol]
MCH4 : Molar weight of the methanem: Modification factor of the wall friction [-]p: Pressure[Pa]: Heat exchanged between liquid
Qls
and solid[W/ms]
Qwall : Heat exchanged between liquid
and wall[W/ms]
R: Gas constant[J/(K . mol)]T: Temperature[K]
U: Overall heat transfer coefficientof the riser pipe
[W/(m . K)]: Superficial velocity[m/sec]Ukv: Velocity[m/sec]w: Distribution factor of the wall friction [kg/m3 . sec]α
: Volume fraction[-]: Averaged gas volume fraction
αgs
for the slug flow[-]
: Solid volume fraction at the inlet [-]α s0γ: Submerge ratio[-]ε: Gas volume fraction averaged
on the whole riser pipe[-]
: Porosity[-]
θpμ: Viscosity[Pa . sec]ρ
: Density[kg/m3]
σCH4 : Saturation ratio of the methane
to the porous part[-]
Subscriptionsk: Focused phase (k= g, l, s)k’: The other phaseg, l, s : Gas phase, liquid phase, and solid phase