ICMF 2010, Tampa, FL, May 30June 4, 2010
Direct numerical simulation of volcano eruption as an explosive bubbly channel
flow using onegridonebubble high resolution sharp interface method
Mingyu Sun* and Chihiro Mikami
Center for Interdisciplinary Research, Tohoku University, Sendai 9808578, 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 highpressure gas sealed
in highviscous liquid in a channel. A conservative VOFtype method has been developed for such a gasliquid
twophase 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 gasliquid 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 onedimensional 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 generalpurposed
solver for compressible twophase flows. The method
allows us to track subgridscale bubbles and droplets
with negligible overhead. The gasliquid interface can
be sharply captured without creating nonphysical 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
gasliquid flows on the eruption process by a 2D 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
gasliquid 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 30June 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 isodensity 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
PLICVOF 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 ~QAt 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 smallcell 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 nonequilibrium phases. These subgridscale
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 twofluid models. An analyt
ical solution has been found for the exchange of mo
mentum, volume and energy to the secondorder 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 nonphysical oscillations. The details
of the mathematical model and numerical method are re
ported in a separate paper (Sun 2010b).
Volumetracking 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 PLICVOF 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 secondorder
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 30June 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 cellcentered
and cellvertex 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 secondorder accuracy in time, the twostep
RungeKutta 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 solutionadaptive 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 solutionadaptive grid for pressure ratio
P,/Pot = 1000
Solutionadaptive unstructured grids In order to sim
ulate a realistic geography near a volcano, the numerical
method is developed for unstructured solutionadaptive
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 30June 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. 59. 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 30June 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. 101113, 2002
Engquist W. E. B., Li X., Ren W, and VandenEijnden,
E., Heterogenous multiscale methods: A review, Com
munications in Computational Physics, 2 (3), (2007)
367450.
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., Volumetracking 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 30June 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 PinPout5
Pin/Pout=lO
SPln/Pout=20
PinIPout=50
SPin/Pout100
Pil/Pout200
SPinPout=500
PidnPout1000
* 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.
offluid (VOF) and twofluid models: basic ideas and
1D 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 HLLRiemann 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 24, pp. 297313 (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. 192203
(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. 303315, 1999
