Group Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Title: 2.4.3 - Direct numerical simulation of volcano eruption as an explosive bubbly channel flow using one-grid-one-bubble high resolution sharp interface method
ALL VOLUMES CITATION THUMBNAILS PAGE IMAGE ZOOMABLE
Full Citation
STANDARD VIEW MARC VIEW
Permanent Link: http://ufdc.ufl.edu/UF00102023/00054
 Material Information
Title: 2.4.3 - Direct numerical simulation of volcano eruption as an explosive bubbly channel flow using one-grid-one-bubble high resolution sharp interface method Computational Techniques for Multiphase Flows
Series Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Physical Description: Conference Papers
Creator: Sun, M.
Mikami, C.
Publisher: International Conference on Multiphase Flow (ICMF)
Publication Date: June 4, 2010
 Subjects
Subject: volcano eruption
bubbly channel flow
conservation
compressible flow
adaptive grid
 Notes
Abstract: Explosive volcanic eruption is driven by the rapid expansion of volatile matters such as water that are dissolved in magma at high temperature and pressure. This phenomenon is numerically modeled as high-pressure gas sealed in high-viscous liquid in a channel. A conservative VOF-type method has been developed for such a gas-liquid two-phase flow. Explosive eruption behavior and liquid fragmentation phenomena have been simulated successfully. The liquid mass flux induced by the eruption is quantitatively measured, and the effects of pressure ratio and void rate are investigated. The liquid mass flux is proportional to the pressure ratio and only weakly depends on gas volume fraction in the range of low void rate.
General Note: The International Conference on Multiphase Flow (ICMF) first was held in Tsukuba, Japan in 1991 and the second ICMF took place in Kyoto, Japan in 1995. During this conference, it was decided to establish an International Governing Board which oversees the major aspects of the conference and makes decisions about future conference locations. Due to the great importance of the field, it was furthermore decided to hold the conference every three years successively in Asia including Australia, Europe including Africa, Russia and the Near East and America. Hence, ICMF 1998 was held in Lyon, France, ICMF 2001 in New Orleans, USA, ICMF 2004 in Yokohama, Japan, and ICMF 2007 in Leipzig, Germany. ICMF-2010 is devoted to all aspects of Multiphase Flow. Researchers from all over the world gathered in order to introduce their recent advances in the field and thereby promote the exchange of new ideas, results and techniques. The conference is a key event in Multiphase Flow and supports the advancement of science in this very important field. The major research topics relevant for the conference are as follows: Bio-Fluid Dynamics; Boiling; Bubbly Flows; Cavitation; Colloidal and Suspension Dynamics; Collision, Agglomeration and Breakup; Computational Techniques for Multiphase Flows; Droplet Flows; Environmental and Geophysical Flows; Experimental Methods for Multiphase Flows; Fluidized and Circulating Fluidized Beds; Fluid Structure Interactions; Granular Media; Industrial Applications; Instabilities; Interfacial Flows; Micro and Nano-Scale Multiphase Flows; Microgravity in Two-Phase Flow; Multiphase Flows with Heat and Mass Transfer; Non-Newtonian Multiphase Flows; Particle-Laden Flows; Particle, Bubble and Drop Dynamics; Reactive Multiphase Flows
 Record Information
Bibliographic ID: UF00102023
Volume ID: VID00054
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: 243-Sun-ICMF2010.pdf

Full Text
ICMF 2010, Tampa, FL, May 30-June 4, 2010


Direct numerical simulation of volcano eruption as an explosive bubbly channel
flow using one-grid-one-bubble high resolution sharp interface method


Mingyu Sun* and Chihiro Mikami

Center for Interdisciplinary Research, Tohoku University, Sendai 980-8578, Japan
s11 1 oCii Lohoku jILc ip
Keywords: volcano eruption, bubbly channel flow, conservation, compressible flow, adaptive grid




Abstract

Explosive volcanic eruption is driven by the rapid expansion of volatile matters such as water that are dissolved in
magma at high temperature and pressure. This phenomenon is numerically modeled as high-pressure gas sealed
in high-viscous liquid in a channel. A conservative VOF-type method has been developed for such a gas-liquid
two-phase flow. Explosive eruption behavior and liquid fragmentation phenomena have been simulated successfully.
The liquid mass flux induced by the eruption is quantitatively measured, and the effects of pressure ratio and void rate
are investigated. The liquid mass flux is proportional to the pressure ratio and only weakly depends on gas volume
fraction in the range of low void rate.


Introduction

The explosive eruption of a volcano is believed as a sud-
den release of a pressurized liquid filled with gas bub-
bles (Mader et al. 1994). Magma is a silicate mixture
containing dissolved volatile gas, such as water, carbon
dioxide and sulfur oxides, and exists in underground at
high temperature and pressure. When it is open to the at-
mosphere by certain reasons, the magma is drivenby the
rapid expansion of volatile matters. The eruption pro-
cess can be modeled as an unsteady gas-liquid two phase
flow with phase changes under sudden reduced pres-
sures. Previous studies on this phenomenon have been
reported, both experimentally (e.g. Cagnoli et al. 2002,
Yamamoto et al. 2008) and numerically (e.g. Wilson
1980, Yoshida & Koyaguchi 1999). However, previous
numerical analysis is based on one-dimensional model-
ing of conduit flows, so that multidimensional character-
istics of the eruption, such as magma defragmentation,
have never been properly reproduced.
Recently, we have developed a general-purposed
solver for compressible two-phase flows. The method
allows us to track subgrid-scale bubbles and droplets
with negligible overhead. The gas-liquid interface can
be sharply captured without creating non-physical oscil-
lations. The advanced solution method based on Rie-
mann solvers has also been implemented in the solver to
resolve better wave dynamics in compressible materials.
The motivation of this study is to use the flow solver
for the simulation of volcano eruptions. The magma de-


fragmentation has been successfully captured. And we
further analyze the effect of various parameters of the
gas-liquid flows on the eruption process by a 2-D mod-
eling.


Numerical modeling and methods

Numerical model of explosive eruption of a volcano A
simplified numerical model is set for the explosive erup-
tion of a volcano, as shown in Fig. la. The liquid is filled
with a matrix of gas bubbles with a pressure higher than
atmosphere, the pressure ratio varying from 1.1 to 1000.
The liquid and the gas are simply taken as water and air
in the present computation. The air is assumed to be the
calorically perfect gas, and water follows the Tait equa-
tion of state. In this work, we are interested in the de-
fragmentation of the liquid phase, so that the viscosity is
neglected. A horizontal free surface separates the bubbly
water below and the atmosphere above. Two sequential
numerical results are shown in Fig. 1 for pressure ratio
of 100. The free surface starting from a straight line de-
forms due to the expansion of bubbles below, and even-
tually breaks up to yield small pieces.
Flow solver that can deal with subgrid interfacial
structures The bubbly liquid is highly pressurized with
gas-liquid interfaces. Both liquid and gas have to treated
as compressible, which imposes tremendous difficulty in
simulation.
In the single phase regions, the flow is governed by




ICMF 2010, Tampa, FL, May 30-June 4, 2010


000ooooo 00000
00000 000000
oo0, o o0000 o 0
00000 00000 00000
a. ^ ^ '^^-^ b ^^ ^ c. iaoi^iA

Figure 1: Modeling the explosive volcano eruption: the
sudden release of a matrix of pressurized underwater
bubbles with a free surface on the top. (a)-(c) are se-
quential iso-density contours.


the Euler equations,

aU OF aP
at s s (1)
where s represents the direction normal to the volume
faces. Conservative quantity U (p, pu,E)T repre-
sents specific mass, momentum and energy respectively.
The flux F is the product of normal velocity and conser-
vative quantity
F uU, (2)
where normal velocity u, = u s. P is surface terms
related to pressure,

P (0, -ps, -pU,)T.

The pressure and normal velocity at volume faces are
then approximated by using an approximate Riemann
solver, HLLC (Toro 1994), from two neighboring pure
phases. Denoting the solutions by and,*, we have

Fe ..* U*, (3)

Pe = (0, -,. s, -pM T. (4)
The mixed cells that contain both gas and liquid
phases are treated in two steps. The first Eulerian step
deals with the interactions between cells, or the length
scale of the grid spacing and larger. In this step, the
mixed cell is approximated by the single phase with the
larger volume (master phase), so that the whole flowfield
is covered only by the cells with pure phases with pos-
sible material interfaces exactly aligned along the grid
lines. If the upstream volume has no interface, the nu-
merical fluxes are then the same as the upwind scheme
for the single phase Euler equations. If the upstream
volume has an interface inside, we have to apportion the
velocity Un in (2) to two phases. This is done by the


Figure 2: Master and slave phases in interface cells. The
pressure and velocity at grid lines are determined only
by the master phases


PLIC-VOF method to be discussed later. Discretizing
(refeqn:qldeuler) with (3) by finite volume method leads
to

U2 + n+l Un U at Qn(a~ U )j,

Un+ + = U n ~Q-At Z(r U, ) jAt (P )j,

where original coefficient .. has been split into a..
and a*..* with


* + = 1.


It is seen that the phase with the smaller volume (slave
phase) in the cell does not affect the pressure and veloc-
ity, so that the small time step or the stiffness due to the
small volume, known as small-cell problem, is avoided.
Therefore, this method can handle fluid parcels smaller
than the grid spacing.
After conducting the Eulerian step, the pressure and
the velocity of two phases in the mixed cell are in gen-
eral not in mechanical equilibriums, which may become
the source of pressure oscillations that eventually pollute
numerical solutions far away. The second step, or the
subgrid modeling, deals with the interactions between
these two non-equilibrium phases. These subgrid-scale
interactions are modeled by two sets of Euler equations
for these phases with interfaces between. The continu-
ous pressure and velocity at the material interfaces are
assumed, and the cell is further imposed with the peri-
odic boundary condition, such that the resulting net flux
is zero. This is similar to the step of the constrained mi-
croscale simulation in heterogeneous multiscale meth-
ods (Engquist et al. 2007). It is similar to the relax-
ation procedure in many two-fluid models. An analyt-
ical solution has been found for the exchange of mo-






mentum, volume and energy to the second-order accu-
racy for both finite and instantaneous relaxations and for
general materials. For the resolved sharp interface, the
subgrid modeling actually forces two phases to develop
toward the equilibrium state with the same pressure and
velocity in the mixed cell. It turns out that an interface
even with a high density ratio can be resolved sharply
without producing non-physical oscillations. The details
of the mathematical model and numerical method are re-
ported in a separate paper (Sun 2010b).
Volume-tracking of subgrid particles The flow solver
requires a tracking method that can handle arbitrarily
small particles or fluid structures. The linear material
interface is reconstructed in the cell, according to the
phase volume and its surface normals (SNs), which are
given by integrating the evolution equations for surface
normals. This unique PLIC-VOF technique allows the
tracking of subgrid particles (Sun 2010a), which signif-
icantly enhances the spatial resolution for the simulation
of defragmentation. The volume of slave phase follows

(as)t + (amU*m)x = 0, (6)
which is the same as incompressible flows, and is up-
dated using the geometric VOF method. The master vol-
ume am is then givenby 1 as.
The equations for SNs are rewritten as

nt + (F1), + (F2)y S, (7)

where
F, (ul F2 v S -m +nvx + ,
F1 =M) F, S .
um vm mU '**

Consider a control volume 8 bounded by discrete faces
System (7) is discretized by the finite volume method.
Consider a control volume Q i bounded by discrete faces
with outward surface normal s = (s sy).

(ni)t- Si ,- 1E


where the numerical flux Fj is approximated by the up-
wind scheme, depending on the direction of normal ve-
locity Un u S,

S F, s + F sy, if ua > 0;
S .F s +s +F2sy, otherwise.

Velocity and surface normal at faces are required to de-
fine the numerical fluxes. For achieving second-order
accuracy in space, the surface normal and velocity at
grid interfaces are located at the center of face, r j They
are interpolated from the center of the volume, following
the MUSCL method,

M-'+ -= M '+ + (VM)'+ (rC r +),


ICMF 2010, Tampa, FL, May 30-June 4, 2010


where M represents both the velocity and surface nor-
mal required to define the numerical flux. Superscripts
+ indicates the values are defined from left and right
sides (or upstream and downstream) respectively. For
example, M '+ are those values located at r'+ re-
spectively. (r r~'+) is the distance between the
center of the grid interface and the location of values
defined. The method is general for both cell-centered
and cell-vertex data structures. In solving the hyper-
bolic conservation laws, the limiter is often used to sup-
press possible numerical oscillations around discontinu-
ities. The MINMOD slope limiter for limiting the gra-
dient (VM)'+ is also implemented. The source terms
are discretized using the central difference scheme. For
achieving second-order accuracy in time, the two-step
Runge-Kutta method is followed for the surface normal
equations. However, the interface reconstruction and the
volume flux evaluation are done only once based on the
surface normal at the last time step.


t=0s t=4.0ms t=5.0ms t=9.9ms
Figure 3: Deformation of material interfaces and cor-
responding solution-adaptive grid for pressure ratio
Pi./Pout = 5


II ,

Si II

t= Os t =2.0ms


I I

II

S4. 'ms t ill ms
t=4.lms t=6.lms


Figure 4: Deformation of material interfaces and cor-
responding solution-adaptive grid for pressure ratio
P,/Pot = 1000

Solution-adaptive unstructured grids In order to sim-
ulate a realistic geography near a volcano, the numerical
method is developed for unstructured solution-adaptive
grids based on quadrilaterals. A few example of these
grids, together with the material interface captured are
shown in Figs. 3 and 4. The finest grids are automati-
cally distributed around the material interfaces. The nu-


H4k


II II I I




ICMF 2010, Tampa, FL, May 30-June 4, 2010


medical algorithm can handle unstructured grids as well,
but not shown in this paper.

Numerical results

The defragmentation phenomenon due to the bubble ex-
pansion is investigated under various pressure ratios.
The sequential density and pressure distributions are
plotted in Figs. 5-9. The flowfields are truly multi-
dimensional, especially for density. The liquid front that
separates the atmosphere above and bubbly liquid below
deforms slowly for low pressure ratios, and the bubble
matrix remains regular. For higher pressures, the contin-
uous liquid soon breaks into small pieces. This defrag-
mentation is more violent at high pressures.


(a) Density distributions(a)t
(d)t=15ms, (e)t= N'.s


Density
[kg/m3]
984
902
820
738
656
574
492
410
328
246
164
82

(d) (e)
=0s, (b)t=5ms, (c)t=10ms,


(a) Density distributions(a)t
(d)t=15ms,(e)t= 2" .s


Density
[kg/m3]
984
902
820
738
656
574
492
410
328
246
164
82
0

(d) (e)
=0s, (b)t=5ms, (c)t=10ms,


Pressure
[MPa]
0.195
0.186
0.173
0.161
0.15
0.14
0.131
0.123
0.116
0 11
0.105


(a) (b) (c) (d) (e)
(3) Pressure distributions(a)t=0s, (b)t=5ms, (c)t=lOms,
(d)t= 15ms, (e)t= 2,, I-
Figure 5: Bubbly channel flows under the pressure ratio
of Pi/Pout=2

The location of leading liquid front (front tip) is mea-
sured for pressure ratios varying from 1.1 to 1000, and
plotted in Fig. 10. It is seen that the front almost moves
at a constant speed for all pressure ratios. By conduct-


20- Pr
[-



1.0 I "

o05 m

0 ) 9 i


essure
MPa]
0.89
0.77
0.66
0.56
0.47
0.39
0.32
0.26
0.21
0.17
0.14
0.12


[1 U 1
(a) (b) (c) (d) (e)
(3) Pressure distributions(a)t=0s, (b)t=5ms, (c)t=10ms,
(d)t=15ms, (e)t=2 '."
Figure 6: Bubbly channel flows under the pressure ratio
of Pi /Pot = 10


ing the linear curve fitting, the tip velocities are obtained,
and plotted in Fig. 11. The tip velocity always increase
with the pressure ratio, and the velocity increases much
more rapidly for low pressure ratios.
The total mass ejected in unit time is also measured,
and summarized in Fig. 12. Runs 1, 2 and 3 corre-
spond to volume fractions = 0.5, 0.13, 0.25 respec-
tively. The mass ejected increases with the pressure, but
slightly decreases at higher volume fractions. However,
even a very small void rate gives five times more mass
rate than a pure pressurized liquid for pressure ratio of
100, which is 6 x 103 (not shown in the figure).


Summary and conclusions

The liquid channel is filled with a matrix of gas bubbles
with a pressure varying from 1.1 to 1000bar. The liquid
and the gas are assumed as water and air in the com-
putation. A horizontal free surface separates the bubbly


II




ICMF 2010, Tampa, FL, May 30-June 4, 2010


2.0 Density
[kg/m3]
984
1.5- 902
820
738
10 656
574
492
0.5 410
328
246
o 164
82
0
-0.5
(a) (b) (c) (d) (e)
(a) Density distributions(a)t=0s, (b)t=5ms, (c)t=10ms,
(d)t=15ms,(e)t=2'". s

20 Pressure


15 4.2
3.6
3.05
1.0 2 55
2.1
[0 51.7
05 1.35
1.05
0.8
0 0.6
0.45
0.35
-5 0.3
(a) (b) (c) (d) (e)
(3) Pressure distributions(a)t=0s, (b)t=5ms, (c)t=10ms,
(d)t=15ms,(e)t=2' ", .
Figure 7: Bubbly channel flows under the pressure ratio
of P, /P,,ot=50


water below and the atmosphere above. The free surface
starting from a straight line deforms due to the expan-
sion of bubbles below, and eventually breaks up to yield
small pieces. The results are summarized as follows.

1. The defragmentation phenomenon is highly un-
steady, and the liquid fragments are launched in a
random manner for pressure ratio over 5.

2. The velocity of the leading front of the liquid in-
creases with the pressure ratio. The velocity is
nearly a constant for a given pressure ratio shortly
after the eruption.

3. The amount of liquid ejected increases with the
pressure ratio.

4. The amount of liquid ejected weakly depends on
the void ratio. However, a few percent of gas bub-
bles may launch about six times more than without
bubbles.


(a) Density distributions(a)t
(d)t=15ms,(e)t=2 N-'s


Density
[kg/m3]
984
902
820
738
656
574
492
410
328
246
164
82
0

(d) (e)
=0s, (b)t=5ms, (c)t=10ms,


20r Pressure
:[MPa]
19.9
1.5 1 17
11.75
9.5
1.0- 7.5
S5.75
S425
0.5 3.0
2.0
1.25
0 0.75
0.5

-05- 0.5
(a) (b) (c) (d) (e)
(3) Pressure distributions(a)t=0s, (b)t=5ms, (c)t=10ms,
(d)t=15ms,(e)t=2- '"s
Figure 8: Bubbly channel flows under the pressure ratio
of Pi,/P,,t=200


References

Cagnoli B, Barmin A, Melnik O, Sparks RSJ, Depres-
surization of fine powders in a shock tube and dynamics
of fragmented magma in volcanic conduits. Earth and
Planetary Science Letters, Vol. 204, pp. 101-113, 2002

Engquist W. E. B., Li X., Ren W, and Vanden-Eijnden,
E., Heterogenous multiscale methods: A review, Com-
munications in Computational Physics, 2 (3), (2007)
367-450.

Mader H., Zhang Y., Phillips JC, Sparks RSJ, Sturtevant
B., Stolper E., Experimental simulations of expolosve
degassing of magma, Nature 372 (1994) 85

Sun, M., Volume-tracking of subgrid particles, to appear
in Intl. J. for Numerical Methods in Fluids, (2010a).

Sun, M., A unified numerical approach for both volume-




ICMF 2010, Tampa, FL, May 30-June 4, 2010


(a) Density distributions(a)t=0s, (b)t
(d)t=15ms,(e)t=2 ,,.s


Density
[kg/m3]
984
902
820
738
656
574
492
410
328
246
164
82
0


:5ms, (c)t=10ms,


Pressure
[MPa]
82
70
59
49
40
32
25
19
14
12
10
7
6
5
4
3.5


S(a) (b) (c) (d) (e)
(3) Pressure distributions(a)t=0s, (b)t=5ms, (c)t=10ms,
(d)t=15ms, (e)t= 2", ,,
Figure 9: Bubbly channel flows under the pressure ratio
of Pi,/Po,,=1000


1.5


N
E 1 L


I* *
0.5

O !O
A*

0 5


0 PidPout 1.1
+PiniPout=2
A PinPout-5
Pin/Pout=lO
SPln/Pout=20
PinIPout=50
SPin/Pout-100
Pil/Pout200
SPinPout=500
PidnPout-1000

* A


10 15 20 25 30 35
t[ms]


Figure 10: Location histories of leading liquid surface
front for pressure ratios varying from 1.1 to 1000.



of-fluid (VOF) and two-fluid models: basic ideas and
1-D tests, submitted for publication, (2010b).


1.2 103


1 103

8 102

g 6 10

4 102

2 102

0
0 2 10 4 102 6 102 8 10 1 1031.2 103


Figure 11: Velocity of leading liquid surface (Tip ve-
locity) as a function of pressure ratio.


1.2 105

1 105


I*RunI
A Run 2
Rmi3


8 104


4104

2 104

'1


A*


0 200 400 600 800 1000 1200
P oP [-]
in out-]

Figure 12: Liquid mass ejected as a function of pressure
ratio.



Toro, E.E, Spruce, M., and Speares, W., "Restoration of
the contact surface in the HLL-Riemann solver", Shock
Waves, vol. 4, 24, 1994.

Wilson L. Relationships between pressure, volatile
content and ejecta velocity in three types of volcanic
explosion, Journal of Volcanology and Geothermal Re-
search, Vol. 8, Issues 2-4, pp. 297-313 (1980)

Yamamoto H., Kazuyoshi Takayama and Kenichi
Ishikawa, Model experiment on magma fragmentation
in explosive volcanic eruption, Journal of Mineralogical
and Petrological Sciences, Vol. 103, No. 3, pp. 192-203
(2008)

Yoshida S., Koyaguchi T, A new regime of volcanic
eruption due to the relative motion between liquid and
gas. Journal of Volcanology and Geothermal Research,
Vol 89, p. 303-315, 1999




University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - - mvs