- Email: [email protected]
Pore-scale numerical simulation of mud erosion in the subsea sand–mud alternate layer using lattice Boltzmann method
Accepted Manuscript Pore-scale numerical simulation of mud erosion in the subsea sand–mud alternate layer using lattice Boltzmann method Takero Yoshida, Takuya Yamaguchi, Hiroyuki Oyama, Toru Sato, Georgios Fytianos, Yuki Kano, Jiro Nagao, Norio Tenma, Hideo Narita PII:
S1875-5100(16)30713-2
DOI:
10.1016/j.jngse.2016.10.003
Reference:
JNGSE 1849
To appear in:
Journal of Natural Gas Science and Engineering
Received Date: 28 April 2016 Revised Date:
2 October 2016
Accepted Date: 5 October 2016
Please cite this article as: Yoshida, T., Yamaguchi, T., Oyama, H., Sato, T., Fytianos, G., Kano, Y., Nagao, J., Tenma, N., Narita, H., Pore-scale numerical simulation of mud erosion in the subsea sand– mud alternate layer using lattice Boltzmann method, Journal of Natural Gas Science & Engineering (2016), doi: 10.1016/j.jngse.2016.10.003. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Pore-scale numerical simulation of mud erosion in the subsea sand–mud alternate layer using lattice Boltzmann method
RI PT
Takero Yoshidaa, Takuya Yamaguchib, Hiroyuki Oyamaa, Toru Satoa, Georgios Fytianosa, Yuki Kanoa,c, Jiro Nagaod, Norio Tenmae, Hideo Naritaf
Department of Ocean Technology, Policy, and Environment, University of Tokyo, 5-1-5
SC
a
Kashiwanoha, Kashiwa 277-8563, Japan
Department of System Innovation, University of Tokyo, 7-3-1 Hongo, Bunkyo-ku 113-8656,
M AN U
b
Japan c
Geological Survey of Japan, National Institute of Advanced Industrial Science and Technology
(AIST), 1-1-1 Higashi, Tsukuba 305-8567, Japan d
Methane Hydrate Project Unit, Research Institute of Energy Frontier, Department of Energy
and Environment, National Institute of Advanced Industrial Science and Technology (AIST), 2e
TE D
17 Tsukisamu-Higashi, Toyohira, Sapporo 062-8517, Japan
Methane Hydrate Project Unit, Research Institute of Energy Frontier, Department of Energy
and Environment, National Institute of Advanced Industrial Science and Technology (AIST), Tsukuba West, 16-1 Onogawa, Tsukuba 305-8569, Japan Department of Energy and Environment, National Institute of Advanced Industrial Science and
EP
f
AC C
Technology (AIST), 2-17 Tsukisamu-Higashi, Toyohira, Sapporo 062-8517, Japan
Corresponding author: T. Sato Tel: +81-4-7136-4726, Fax: +81-4-7136-4726, E-mail:[email protected] 1
ACCEPTED MANUSCRIPT
ABSTRACT
Commercial-scale extraction of methane gas from subsea methane hydrate is difficult. A major challenge is mud erosion at the interfaces of sand–mud alternate layers in the subsea strata.
RI PT
Without a clear understanding of the various parameters associated with mud erosion, there may be complications in developing methane hydrate deposits, such as reservoir clogging due to mud particles. This study aims to develop a numerical simulation method to predict mud erosion caused by water flow in pore-scale computational domain. Microscopic sand beds were
SC
generated numerically in the computational domain while the water flow through the pores of various-shaped sand grains was numerically simulated using the lattice Boltzmann method. Mud
M AN U
was treated as a continuum and eroded at its surface via the shear stress of water flow using the critical shear stress and erosion rate constant, which were estimated by fitting calculated erosion rates with those of a laboratory-scale experiment. Using this numerical model, simulations of long-term mud erosion were conducted. The simulation showed that when the average velocity of water in the pore space is larger than the critical value and smaller than ten times that value,
TE D
the erosion rate is almost zero after about 10 days.
KEYWORDS: Mud erosion; Sand–mud alternate layer; Lattice Boltzmann method; Subsea
AC C
EP
sediment; Critical shear stress; Erosion rate.
2
ACCEPTED MANUSCRIPT
1. INTRODUCTION Mud erosion by water and/or gas flows at the interface of sand–mud alternate layers may cause problems for subsea engineering processes such as oil and gas production. In particular, the production of methane gas by depressurization from methane hydrate, which is found in sub-
RI PT
seabed sand layers that alternate with mud layers in turbidite formations (Fujii et al., 2015), causes a two-phase flow through the sand layers that comprises dissociated gas and water. This can potentially result in mud erosion of and, possibly, blockage of the flow.
Several types of mud erosion for fine-grained sediments were examined by Amos et al.
SC
(1992). Winterwerp et al. (2004) classified the erosion models for cohesive sediments: floc erosion, entrainment, surface erosion, and mass erosion. In terms of erosion of sand–mud
M AN U
mixtures, Mitchener et al. (1996) stated that the critical shear stress, τc, increases when sand is added to mud, where τc is the smallest local shear stress to cause erosion. Moreover, they showed that τc for sand–mud mixtures with sand contents of 70%–100% depends on the size of the sand grains and on the cohesive properties of the mud materials. Detailed studies of erosion types for sand–mud mixtures were conducted by Jacobs et al. (2011). After investigating different soil samples with varying compositions, they categorized erosion of sand–mud mixtures into two
TE D
types: time-decreasing and time-independent erosion.
The above studies are mainly limited to mud erosion on sea floors or river bottoms. For mud erosion at the interface of sand–mud alternate layers, Brumby et al. (2015) conducted numerical simulations using a combination of the lattice Boltzmann method (LBM) (based on
EP
the approaches of Sugita et al. (2012) and Guo et al. (2013)) and the discrete element method (DEM) developed by Cundall et al. (1979). Based on erosion studies on cohesive sediments by
AC C
Partheniades (1965), Oyama et al. (2016) introduced an equation to express mud erosion rate:
(
α exp U − U c E= 0
)
U ≥ Uc
U < Uc
,
(1)
where E [kg/m2/s] is the erosion rate, αexp [kg/m3] is the erosion rate constant, U [m/s] is the average velocity of water, and Uc [m/s] is the critical velocity, which is the erosion threshold. The reason why the average water velocity was chosen over the Darcy velocity is that mud erosion is directly relevant to the shear stress on the mud surface caused by water flow. Thus, we 3
ACCEPTED MANUSCRIPT
decided that the average water velocity would be more appropriate to represent the shear stress than the Darcy velocity, which is affected by porosity. The erosion rate constant and critical velocity were estimated by Oyama et al. (2016), who performed experiments of mud erosion using artificial sand–mud alternate cores.
RI PT
The primary objective of this study was to develop a numerical simulation method to realize mud erosion in sand–mud alternate layers. We demonstrated pore-scale numerical
simulations of water-induced erosion using the LBM. First, we generated a numerical pore-scale domain that mimics a structure composed of several frame-sand grains. At the bottom of this
SC
microscopic computational domain, we filled the pore space with mud. Next, water flow through this porous domain was analyzed using the LBM. Then, based on the resultant water velocity, we
M AN U
calculated the shear stress acting on the surface of the mud layer and changed the shapes of the mud surface to represent the mud erosion in the LBM.
Although of great interest, the temporal change of the erosion rate over long periods of time was previously unknown; however, the results of this study’s long-time simulations show how the erosion rate changes with time. Therefore, the results of our numerical simulation may be able to offer important insights that will be valuable for the gas production process from
TE D
methane hydrate reservoirs.
2. NUMERICAL METHOD
2.1 Generation of the sand bed
EP
Sugita et al. (2012) extracted the shapes of the sand grains using computed tomography images of sand cores obtained from hydrate sedimentary layers off Japan’s coast. Thereafter, they
AC C
packed these grains in a pore-scale computational domain using a numerical growth method to attain the required porosity. Figure 1 illustrates an example of a generated sand bed with six grains, which has a porosity of 0.40 under a periodic boundary condition. Because of the limitations of our computational facilities, the length of each side of the domain was set to be 200 µm. Because the domain is so small that it cannot represent the characteristics of the whole sand layer, we generated 10 domains and took averages of the properties obtained from the calculations (details presented later in section 2.6).
4
RI PT
ACCEPTED MANUSCRIPT
Fig. 1 An example of a three-dimensional computational domain holding six sand grains with a
SC
periodic boundary condition.
2.2 Lattice Boltzmann method (LBM)
M AN U
Porous substances are structurally complicated, making the creation of a computational matrix difficult with conventional computational fluid dynamics (CFD) methods. The LBM is a powerful tool for the calculation of complex liquid flows with low Reynolds number. The principal advantages of using LBM for the simulation of flow through porous media have been confirmed by many researchers, e.g., Guo et al. (2002) confirmed the suitability of the LBM for the simulation of incompressible flow through porous media. The LBM has additionally been
TE D
used for a number of other applications, e.g., a study conducted by Jain et al. (2016) explored the conditions for transition to turbulence in intracranial aneurysms at spatial and temporal resolutions near the Kolmogorov microscales. LBMs for microfluidics were well reviewed by Zhang (2011). Sugita et al. (2012) also developed a simulation code to calculate water flows
(2012).
EP
within pore-scale porous media. In this study, we applied the LBM developed by Sugita et al.
AC C
For solid boundaries with curved surfaces, we adopted the boundary condition proposed
by Bouzidi et al. (2001), which can compensate for the staircase approximation in the LBM. There are various models for discretizing velocity, and in this study, we adopted the threedimensional (3D) fifteen-velocity model developed by Qian (1992) (D3Q15). Nash et al. (2014) confirmed the usage of this model in complex flow domains, and the method of Sugita et al. (2012) was very stable when using the D3Q15 stencils together with the Bhatnagar–Gross– Krook (BGK) model for the collision term (Bhatnagar et al., 1954). Using the BGK model for
5
ACCEPTED MANUSCRIPT
the collision term, we can express the particle distribution function fi with the discrete Boltzmann equation (MacNamara, 1988): ∂fi 1 + ci ⋅ ∇f i = − (f i − f ieq ), ∂t τ
RI PT
(2)
where ci indicates the virtual particle velocity in the direction towards i, which denotes each of the surrounding 14 lattice points, t is time, and τ is the relaxation time. The local equilibrium
is the general kinetic theory of fluids expressed as
c=
M AN U
3 9 3 2 f ieq = ρw i 1 + 2 c i ⋅ u + 4 (c i ⋅ u ) − 2 u ⋅ u , 2c 2c c
SC
distribution function fieq is an extension of the Maxwell distribution to a local equilibrium, which
∆x , ∆t
(3)
(4)
TE D
where ∆t [s] is the time step, ∆x [m] is the lattice unit scale, and wi is the weighting factor depending on the discrete velocity, ci. The fluid density ρ [kg/m3] and fluid velocity u [m/s] are
ρ = ∑f i , i
i
(5) (6)
AC C
u = ∑cifi .
EP
given by
To correspond to the Navier–Stokes equation, the relaxation time τ is set to satisfy the
following equation: 1 ν = τ − c s2 ∆t, 2
(7)
6
ACCEPTED MANUSCRIPT
where ν [m2/s] is the kinematic viscosity of water. Calculations are performed separately for the collision and streaming processes using the
1 f i ( x + c i ∆t, t + ∆t ) = fi ( x, t ) − [fi ( x, t ) − f ieq ( x, t )], τ f i ( x + c i ∆t , t + ∆t ) = fi (x, t ).
RI PT
following equations:
(8)
(9)
SC
2.3 Calculation of shear stress
Forces acting on the surface of the continuum mud can be calculated via the momentum exchange method (e.g., Mei et al., 2002; Ladd, 1994), which is a simple implementation for fluid
M AN U
flows, assuming that the surface of the continuum is a solid boundary. The force originating from momentum exchange per unit time between neighboring lattices is expressed as follows.
Fi =
1 [c i f i (x) − c −i f −i (x + c i ∆t) ], ∆t
(10)
TE D
where F [kg m/s2] is the force vector acting on the mud surface and −i denotes the reflected direction, which is the bounce-back direction, as opposed to the direction of i. To evaluate the shear stress on a corner mud lattice, which has more than one side exposed to water lattices, F on
EP
the corner mud lattice is calculated using (10) as the summation of Fi, which is toward the water lattices. The tangential force on the mud surface can be expressed using the following equation:
n=−
AC C
Ft = F − Fn = F − F ⋅ n ,
∑Q c
i i
i
2
,
(11)
(12)
∑Q c
i i
i
where Ft and Fn are the force components of the tangential and normal direction of the mud surface, respectively, n is the normal unit vector, and Qi is the amount of mud in a lattice cell 7
ACCEPTED MANUSCRIPT
with direction i. The shear stress can be obtained by dividing the tangential force with the area of the mud surface.
2.4 Validation of the shear stress calculations
RI PT
To determine the validity of the calculated shear stress values, numerical simulations were conducted for the flow between parallel flat plates, and the calculated shear stresses were
compared with a theoretical value of the Poiseuille flow. The exact value of the shear stress can
1 dP h, 2 dx
M AN U
τs =
SC
be expressed via the following equation:
(13)
where dP/dx [Pa/m] is the pressure gradient and h [m] is the distance between the two parallel plates. Meanwhile, the exact average velocity Ua [m/s] is given by
1 dP 2 h , 12 µ dx
(14)
TE D
Ua =
where µ [kg/m/s] is the viscosity of water.
The simulation conditions were as follows. The flat plate is shown in black in Fig. 2. The
EP
domain sizes were 50, 50, and 200 µm in the x, y, and z directions, respectively. A periodic boundary condition was imposed for the side boundaries. The surface of the flat plate had a noslip condition. The upper surface of the domain was a symmetric boundary so that the height of
AC C
the domain was equal to h/2. A pressure gradient of 1000 Pa/m was imposed in the x-direction in the whole domain.
Here, the scale of the computational lattices and the relaxation time were examined by
comparing them to the theoretical values calculated using (13) and (14). Figures 3 and 4 show the comparison with the average velocity and wall shear stress of the Poiseuille flow with various relaxation times and computational lattice sizes. According to Kruger et al. (2009), their LBM with a relaxation time of 0.9 provides small errors for the velocity and shear stress. Kruger et al. (2009) demonstrated the cases of Re = 1, 10, and 100, and Ma = 0.01, 0.04, and 0.20. The
8
ACCEPTED MANUSCRIPT
number of the lattice cells in their computations ranged from 16 × 18 × 18 to 128 × 130 × 130, and the total computational time was between 0.25 and 0.75 s. In our numerical simulations, Re and Ma were approximately 2 and 0.03. The number of the lattice cells in our computations
and 1.5 s. The water velocity was set to 0 as the initial condition.
RI PT
ranged from 10 × 10 × 40 to 25 × 25 × 100, and the total computational time was between 0.13
Figures 3 and 4 agree well with the suggestion by Kruger et al. (2009), although all the simulated results are in moderately good agreement with the theory. In this study, a relaxation time of 1.0 and a computational lattice with dimensions of 2.5 µm were adopted for the
EP
TE D
M AN U
shear stress can be computed with an error margin of <1%.
SC
following numerical simulation. Using these parameters, the average water velocity and the wall
Fig. 2 Computational domain for the Poiseuille flow. The lower part of a parallel flat plate is
AC C
shown in black. The upper surface is a symmetric boundary and a periodic condition is imposed on the side surfaces.
9
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
Fig. 3 Calculated average velocity compared with the theoretical values using different
AC C
EP
TE D
relaxation times and lattice sizes.
Fig. 4 Calculated wall shear stress compared with theoretical values using different relaxation times and lattice sizes. 10
ACCEPTED MANUSCRIPT
2.5 Numerical model of mud erosion Following the erosion-rate model, (1) used by Oyama et al. (2016), we also adopted a linear
α cal (τ s − τ c ) E= 0
τs ≥ τc
,
(15)
SC
τs < τ c
RI PT
relationship between the local erosion rate and the shear stress, given by
where E [kg/m2/s] is the local erosion rate, αcal [s/m] is the erosion rate constant, τs [Pa] is the
M AN U
shear stress, and τc [Pa] is the critical shear stress being applied on the mud. Mud erodes when the shear stress acting on its surface is larger than the critical value.
To match (15) with (1), the average erosion rate, E [kg/m2/s], in the numerical simulation is given by
1 E dA, A τ s ∫≥ τ c
(16)
TE D
E=
where A is the total area of the sand–mud interface, including the sectional areas of the sand grains that penetrate the interface. It is noted that (1) and (15) do not hold at the same time
EP
because the integration in (16) is only done when τ ≥ τ and such area varies depending on the average water velocity. Therefore, the linearity in (15) never assures that in (1). However, (1) is
AC C
also meaningful from an engineering viewpoint because it can be easily applied to actual gas or oil production operations.
When the volume of mud in a lattice cell is larger than half the volume of a lattice cell
before erosion and becomes smaller than half after the erosion, we consider that the lattice originally defined as a mud lattice has become a water lattice. In such a case, the distribution function of this lattice is set to 0, and consequently, the water velocity of this lattice becomes 0 because the mud was not originally moving before the erosion. Thus, local cell-wise mass conservation does not hold with this numerical treatment. However, the mud volume eroded is assumed to be very small, <0.1% of the initial volume after long-time simulation as shown later, 11
ACCEPTED MANUSCRIPT
and we consider that it is negligible.
2.6 Estimation of critical shear stress and erosion rate constant The critical shear stress and erosion rate constant in the numerical simulation were estimated
RI PT
based on experimental data. Oyama et al. (2016) conducted laboratory-scale experiments on mud erosion under varying water velocities. To create an artificial core that mimics sand–mud
alternate layers, they packed kaolinite as mud particles and Toyoura sand as large frame grains in each half of a cylinder of length 100 mm and an inner diameter of 49.8 mm. Water was injected
SC
into the core, and the discharged water containing eroded mud was sampled at set time intervals. After analyzing the results of the experiment, Oyama et al. (2016) obtained a critical velocity, Uc,
M AN U
of 1.81 × 10−5 m/s and an experimental erosion rate constant, αexp, of 0.00374 kg/m3 at a confining pressure of 4 MPa.
Although measuring shear stress directly on the mud surface is difficult, we need to input the critical shear stress in our mud erosion simulation. Here, numerical simulations were conducted under the condition that the average water velocity is set to the critical velocity, Uc, by adjusting the pressure gradient to identify the critical shear stress, τc. The computational
TE D
conditions are listed in Table 1. The mud layers were set at the bottom of the computational domains, as shown in Fig. 5. In these simulations, the mud was not eroded, and it was treated as a fixed solid with its surface having no no-slip condition. A periodic boundary condition was applied to the side surfaces. The upper surface is a boundary with a gradient of zero in the
EP
vertical direction.
AC C
Table 1. Computational condition
Scale of calculation area
Number of lattices Porosity
200 × 10−6 [m] 80 × 80 × 80 0.4
Lattice size
2.5 × 10−6 [m]
Time step
1.04 × 10−6 [s]
12
SC
RI PT
ACCEPTED MANUSCRIPT
M AN U
Fig. 5 Isometric view of a computational domain at the initial stage. Sand, mud, and water are shown as gray, black, and white, respectively. Mud is set in the pore space close to the bottom of the domain.
According to Sugita et al. (2012), the Darcy velocity was validated when 30 lattices were used for a glass-bead particle with a diameter of 100 µm. The sizes of the computational domain
TE D
and the sand grains used by Sugita et al. (2012) are identical to those used in the present simulation, and therefore, we consider that the lattice size is sufficiently small for an accurate analysis.
EP
Because our computational domain includes only six sand grains and pore-scale mud erosion depends on the shapes and locations of the grains, we conducted flow simulations using 10 different domains with different sand-packing patterns and took average of the resultant data
AC C
of the 10 calculations to determine meaningful parameters. Figure 6 shows the cumulative frequency of shear stress on the mud surface of the 10
computational domains. As the average water velocity was set to be Uc, there should be no mud erosion and the shear stress should be smaller than τc. Here, we determined τc to be the value at which the cumulative area of the mud surface is 99.73%, corresponding to the 3σ upper limit of the Gauss distribution. Eventually, τc is 0.0233 Pa, as shown in Fig. 6.
13
100
99.73 ％
80 60 40 20
0.0233 Pa
0 0
0.005
0.01 0.015 0.02 Shear stress [Pa]
0.025
RI PT
Cumulative frequency [%]
ACCEPTED MANUSCRIPT
0.03
SC
Fig. 6 Cumulative frequency of shear stress calculated using 10 sand patterns under the condition
M AN U
that the average water velocity is the same as the critical velocity.
Next, we examined the effect of the mud phase height in the computational domains, as shown in Fig. 7 (sand, mud, and water are displayed in gray, black, and white, respectively). Here, we imposed a pressure gradient of 1000 Pa. Figure 8 shows the relationship between the average of the calculated shear stress acting on the flat mud surface and the mud heights under the same pressure gradient; the mud height does not exceed 40 lattice cells and the shear stresses
TE D
are almost constant. It is therefore fair to say that the initial height of mud being two lattice cells
AC C
EP
in this mud erosion simulation is validated.
(a)
(b)
14
RI PT
ACCEPTED MANUSCRIPT
(c)
(d)
SC
Fig. 7 Examples of computational domains with different heights of mud phase. The sand, mud, and water are displayed in gray, black, and white, respectively. The number of mud
EP
TE D
M AN U
lattices are (a) 2, (b) 20, (c) 40, and (d) 60 in the vertical direction.
Fig. 8 Calculated shear stress with various mud lattice heights in the vertical direction, as shown
AC C
in Fig. 7. The plots are the average using 10 different sand patterns.
Finally, the erosion rate constant was estimated by matching the erosion rate with the
experimental data for kaolinite content of 100%. The average velocities of water in the numerical simulation were set to be the same as those of the experiment, i.e., water injection rates of 0.5, 0.7, 1, 2, and 3 g/min, by adjusting the pressure gradient. The simulations were conducted for each of these water velocities using 10 sand patterns. The errors in the erosion rate between the numerical simulation and the experiment were minimized by adjusting the erosion rate constant using the least squares method. The best-fit erosion rate constant, αcal, was 7.47 × 10−5 s/m. The
15
ACCEPTED MANUSCRIPT
erosion rates calculated using this value are shown in Fig. 9. Although the calculated and measured erosion rates are moderately identical, the general trend is that the calculated rates are slightly smaller than those in the experiments. This may be acceptable because the erosion rate
M AN U
SC
RI PT
totally depends on the sand grain arrangement in the domain.
Fig. 9 Comparison of the calculated erosion rates with the experimental values obtained by Oyama et al. (2016).
TE D
3. TEMPORAL CHANGE OF THE EROSION RATE OVER LONG TIME Now, using the values of αcal and τs obtained above, a decrease in the mud volume can be calculated for long time periods, e.g., several weeks. The computational conditions for long-term mud erosion are shown in Table 2.
EP
To conduct a long-term simulation after the mud surface reaches the bottom of the domain, we move the top lattice layer below the computational domain by filling mud into its
AC C
whole pore whenever a part of the eroded mud surface touches the bottom of the computational domain. This treatment poses no difficulty because the sand grains are packed in a periodic way in the vertical direction.
In addition, we enlarged the erosion rate constant, αcal, by 107 times unless the volume of
eroded mud per time step was still less than that of the computational cell. This is because too much time is needed to calculate erosion for a long time range and because the physics is theoretically the same if αcal × t is the same. As shown in Fig. 10, the results using one sand pattern indicate that using magnifications of 106 and 107 produce almost identical results. In contrast, using magnifications of 108 and 109 16
ACCEPTED MANUSCRIPT
produced different results compared with using magnifications of 106 and 107. As a result, using a larger magnification can result in an incorrect erosion rate numerically, although theoretically, the magnification should not have an effect on the result. This may be because if the volume of the mud eroded in one time step is as large as the lattice cell volume, the eroded volume in a cell
RI PT
many remain cubic-cell-shaped, which is an artifact. Based on the results in Fig. 10, we adopted
M AN U
SC
a magnification of 107 to save computational time.
Fig. 10 Temporal changes of the calculated erosion rate in 3 days using one sand pattern
TE D
with varying magnifications for the erosion rate constant, αcal: 106, 107, 108, and 109. Table 2. Numerical simulation conditions
EP
Scale of calculation area
200 × 10−6 [m]
Number of lattices
80 × 80 × 80
Porosity
0.4 2.5 × 10−6 [m]
Time step
1.04 × 10−6 [s]
Critical shear stress
0.0233 [Pa]
Erosion rate constant
7.47 × 10−5 [s/m]
Bulk density of mud (porosity = 0.4)
1110 [kg/m3]
Average velocity of water
2.5Uc and 7.5Uc
Total time into real scale
15 [days]
AC C
Lattice size
17
ACCEPTED MANUSCRIPT
Figure 11 shows the temporal change of the erosion rates, which are average values of the results using 10 sand patterns, for one day. In this and subsequent figures, the erosion rate and time are modified such that they are small and long by dividing and multiplying by the expansion ratio of αcal, 107, respectively. The erosion rates did not show significant changes with time, at
RI PT
least within several hours. Oyama et al. (2016) assumed that the erosion rate was constant during their experiment, which lasted up to several hours, to extract important erosion parameters.
(a) Average water velocity was 2.5Uc
M AN U
SC
Evidently, the present numerical simulation supports their assumption.
(b) Average water velocity was 7.5Uc
Fig. 11 Temporal changes of the calculated erosion rates in 1 day. The results are the average
TE D
values of 10 different sand patterns. The average water velocities were 2.5Uc and 7.5Uc, where Uc represents the critical velocity that induces mud erosion. Temporal changes in the calculated erosion rates, which are average values of the results
EP
using 10 sand patterns, over long periods of time are shown in Fig. 12. The erosion rate decreases monotonously with time with small fluctuations (almost zero). These fluctuations seem
AC C
to depend on the morphology of the sand and mud. To understand this temporal change in the erosion rate, the velocity vectors and mud locations, shown in black, are depicted in Fig. 13 on the 1st and 10th days on a vertical section of one of the 10 domains. The velocity vectors are drawn at every 3 computational lattices. On the 1st day, the mud was continuously eroded around point A in Fig. 13(a). On the 10th day, the water path around point A in Fig. 13(b) was broadened because of erosion, and the flow velocities right above the eroded mud surface decreased. When spaces between the sand grains and the mud surfaces became sufficiently large, the shear stress was smaller than τc along such channels, and there was no further erosion (Fig. 13(b)). In contrast, a small amount of the mud was eroded around point B in Fig. 13(a) and (b). 18
ACCEPTED MANUSCRIPT
This may be because point B is located in the dead wake of a sand grain and thus there is no
M AN U
SC
RI PT
potential for erosion. This is another mechanism through which erosion can cease.
EP
TE D
(a) Average water velocity was 2.5Uc
(b) Average water velocity was 7.5Uc
AC C
Fig. 12 Temporal changes of the calculated erosion rates during 15 days. The same graph with the time period ranging from the 5th to 15th day is enlarged (inset).
19
RI PT
ACCEPTED MANUSCRIPT
(b) 10th day
SC
(a) 1st day
Fig. 13 Vertical sections of the computational domain on the 1st and 10th day for one sand
M AN U
pattern. Sand, mud, and water are shown as gray, black, and white, respectively. The arrows indicate the velocity vectors of water flow [m/s].
Figure 14 shows the histograms of the calculated shear stress distribution when the average water velocity was 7.5Uc on the 0th and 15th days. The result shows that shear stresses larger than τc, 0.0233 Pa, disappeared by the 15th day. The decrease in the erosion rate may be
TE D
attributed to the water flow eroding the mud so that smooth flow channels are formed near the
AC C
EP
deformed mud surface, as can be seen in Fig. 13.
Fig. 14 Histogram of shear stress acting on mud surface in the case where average water velocity is 7.5Uc. The solid and dotted lines show the results of 0th and 15th day, respectively.
Overall, when the average velocity of water, U, is several times larger than the critical 20
ACCEPTED MANUSCRIPT
value, our results suggested that mud erosion occurs only at an early stage of the gas production process from methane–hydrate bearing layers (approximately within 10 days). Oyama et al. (2016) reported that microfractures are formed in sand layers when the average velocity is larger than 2.14 × 10−4 m/s, which is approximately 10Uc. According to Oyama et al. (2010), the water
RI PT
velocity should be less than approximately 10−4 m/s, even in the vicinity of the wellbore;
otherwise, the sand frame is broken and sand is produced. Here, let us assume that the radius of a wellbore is R and water velocity around it is 2 × 10−4 m/s, approximately 10Uc. We can then estimate water velocity at 10R from the wellbore as 2 × 10−5 m/s, approximately Uc. Thus, mud
SC
erosion is considered to occur within the circle with a radius of 10R over a time period of <10 days. Assuming that the radius of a production well is 0.15 m, the cumulative mass of the eroded
M AN U
mud in a radius of 10R within 10 days can be estimated approximately as less than 1 kg at an interface of a sand–mud alternate layer. As a result, this estimation suggests that the total mass of the eroded mud is trivial even if there are multiple sand–mud layers. In this study, the mud was treated as a continuum solid that changes its shape because of erosion using (15), and it is assumed that the eroded mud does not affect any flow properties because the concentrations of mud particles are supposedly very less.
TE D
The present numerical simulation has some drawbacks, namely, sand grains do not move after the mud bases have eroded. A possible reason for this is that as in Fig. 13, the eroded mud in the figure can be seen as approximately 50 µm in depth. Similarly, in most cases, the mud surfaces were eroded 50–100 µm in depth before the erosion stopped. This result suggests that
EP
the thickness of the eroded mud within 10 days is smaller than the size of the sand grains. In addition, the sands originally formed a frame and do not necessarily move. In future, we intend
AC C
to demonstrate the decomposition of frame-sand grains as well as the formation of microfractures in sand layers by incorporating the sand grains into the present mud erosion simulation.
4. CONCLUSIONS
A numerical simulation method using LBM was developed to estimate the mass of mud erosion at the interface of sand–mud alternate layers. The sand beds were numerically generated in microscopic computational domains and water flow was forced within the pore space by a constant pressure gradient. First, we confirmed the accuracy of this numerical method by 21
ACCEPTED MANUSCRIPT
accurately comparing the calculated wall shear stress with that of the Poiseuille flow. Then, mud erosion was simulated for the generated microscopic sand–mud layers using experimentally obtained erosion parameters. We studied the validity of the present method by comparing the calculated erosion rate with that of the experimentally obtained value.
RI PT
Next, we applied the present method to a long-term erosion simulation. The calculated results suggest that the erosion rate decreases to zero in less than 10 days under the condition that the average velocity is several times larger than the critical value for mud erosion. This may be because the mud surface with large shear stress is eroded first, resulting in moderately wide
SC
spaces appearing in the pore of the sand grains, which slows the water velocity to a value smaller
AC C
EP
TE D
M AN U
than the critical value.
22
ACCEPTED MANUSCRIPT
ACKNOWLEDGEMENTS
This study was financially supported by the Research Consortium for Methane Hydrate Resources in Japan (MH21 Research Consortium) of Methane Hydrate Research and
AC C
EP
TE D
M AN U
SC
RI PT
Development Program planned by the Ministry of Economy, Trade and Industry (METI).
23
ACCEPTED MANUSCRIPT
REFERENCES
Amos, C.L., Daborn, G.R., Christian, H.A., Atkinson, A., Robertson, A., 1992. In situ erosion
RI PT
measurements on fine-grained sediments from the Bay of Fundy. Mar. Geol. 108, 175-196.
Bhatnagar, P.L., Gross, E.P., Krook, M., 1954. A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems. Phys. Rev. 94, 511-525.
M AN U
fluid with boundaries. Phys. Fluids 13, 3452-3459.
SC
Bouzidi, M., Firdaouss, M., Lallemand, P., 2001. Momentum transfer of a Boltzmann-lattice
Brumby, P. E., Sato, T., Nagao, J., Tenma, N., Narita, H., 2015. Coupled LBM DEM micro-scale simulations of cohesive particle erosion due to shear flows. Transp. Porous Med. 109, 43-60.
Cundall, P.A., Strack, O.D.L., 1979. A discrete numerical model for granular assemblies.
TE D
Géotechnique 29, 47–65.
Fujii, T., Suzuki, K., Takayama, T., Tamaki, M., Komatsu, Y., Konno, Y., Yoneda, J., Yamamoto, K., Nagao, J., 2015. Geological setting and characterization of a methane hydrate reservoir distributed at the first offshore production test site on the daini-atsumi knoll in the
EP
eastern nankai trough. Japan. Mar. Petrol. Geol. 66, 310-322.
AC C
Guo, Z., Zhao, T. S., 2002. Lattice Boltzmann model for incompressible flows through porous media. Phys. Rev. E 66, 036304.
Guo, Z., Shu, C., 2013. Lattice Boltzmann method and its applications in engineering.World Scientific Publishing Co. Pte. Ltd., Singapore.
Jacobs, W., Hir, P. L., Kesteren, W. V., Cann, P., 2011. Erosion threshold of sand/mud mixtures. Continental Shelf Res. 31, 14-25.
24
ACCEPTED MANUSCRIPT
Jain, K., Roller, S., Mardal, K. A., 2016. Transitional flow in intracranial aneurysms - A space and time refinement study below the Kolmogorov scales using Lattice Boltzmann Method. Comput. Fluids 127, 36-46.
RI PT
Kruger, T., Varnik, F., Raabe, D., 2009. Shear stress in lattice Boltzmann simulations, Phys. Rev. E 79, 046704.
Ladd, A., 1994. Numerical simulations of particulate suspensions via a discretized Boltzmann
SC
equation. Part 1. Theoretical foundation, J. Fluid Mech. 271, 285-309.
M AN U
McNamara, G., Zanetti, G., 1988. Use of the Boltzmann equation to simulate lattice-gas automata. Phys. Rev. Lett. 61, 2332-2335.
Mei, R., Yu, D., Shyy, W., Luo, L-S., 2002. Force evaluation in the lattice Boltzmann method involving curved geometry. Phys. Rev. E 65, 041203.
TE D
Mitchener, H., Torfs, H., 1996. Erosion of mud/sand mixtures. Coast. Eng. 29, 1-25.
Nash, R. W., Carver, H. B., Bernabeu, M. O., Hetherington, J., Groen, D., Krger, T., Coveney, P. V., 2014. Choice of boundary condition for lattice-Boltzmann simulation of moderate-Reynolds-
EP
number flow in complex domains. Phys. Rev. E 89, 023303.
AC C
Oyama, H., Abe, S., Yoshida, T., Sato, T., Nagao, J., Tenma, N., Narita, H., 2016. Experimental study of mud erosion at the interface of artificial sand-mud alternate layer. J. Nat. Gas Sci. Eng. 34, 1106–1114.
Oyama, H., Nagao, J., Suzuki, K., Narita, H., 2010. Experimental analysis of sand production from methane hydrate bearing sediments applying depressurization method. J MMIJ 126, 497– 502 (in Japanese).
Partheniades, E., 1965, Erosion and deposition of cohesive soils. J. Hydraul. Div. 91, 105-139. 25
ACCEPTED MANUSCRIPT
Qian, Y. H., D’Humires, D., Lallemand, P., 1992. Lattice BGK models for the Navier-Stokes equation. Europhys. Lett. 17, 479-484.
RI PT
Sugita, T., Sato, T., Hirabayashi, S., Nagao, J., Jin, Y., Kiyono, F., Ebinuma, T., Narita, H., 2012. A pore-scale numerical simulation method for estimating the permeability of sand sediment. Transp. Porous Med. 94, 1-17.
SC
Winterwerp J. C., van Kesteren, W.G.M., 2004. Introduction to the physics of cohesive sediment
M AN U
in the marine environment. Chapter 9, Elsevier, Amsterdam.
Zhang, J., 2011. Lattice Boltzmann method for microfluidics: models and applications.
AC C
EP
TE D
Microfluid Nanofluid 10, 1-28.
26
ACCEPTED MANUSCRIPT
Highlights
A pore-scale numerical simulation method for mud erosion was developed.
•
Microscopic sand beds are generated numerically in the computational domain.
•
Water flow through the sand pore is simulated using a lattice Boltzmann method.
•
Numerically, mud is eroded by the water shear stress using erosion rate constant.
•
Long-term simulation suggests that erosion rate approaches zero in about 10 days.
AC C
EP
TE D
M AN U
SC
RI PT
•
S1875-5100(16)30713-2
DOI:
10.1016/j.jngse.2016.10.003
Reference:
JNGSE 1849
To appear in:
Journal of Natural Gas Science and Engineering
Received Date: 28 April 2016 Revised Date:
2 October 2016
Accepted Date: 5 October 2016
Please cite this article as: Yoshida, T., Yamaguchi, T., Oyama, H., Sato, T., Fytianos, G., Kano, Y., Nagao, J., Tenma, N., Narita, H., Pore-scale numerical simulation of mud erosion in the subsea sand– mud alternate layer using lattice Boltzmann method, Journal of Natural Gas Science & Engineering (2016), doi: 10.1016/j.jngse.2016.10.003. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Pore-scale numerical simulation of mud erosion in the subsea sand–mud alternate layer using lattice Boltzmann method
RI PT
Takero Yoshidaa, Takuya Yamaguchib, Hiroyuki Oyamaa, Toru Satoa, Georgios Fytianosa, Yuki Kanoa,c, Jiro Nagaod, Norio Tenmae, Hideo Naritaf
Department of Ocean Technology, Policy, and Environment, University of Tokyo, 5-1-5
SC
a
Kashiwanoha, Kashiwa 277-8563, Japan
Department of System Innovation, University of Tokyo, 7-3-1 Hongo, Bunkyo-ku 113-8656,
M AN U
b
Japan c
Geological Survey of Japan, National Institute of Advanced Industrial Science and Technology
(AIST), 1-1-1 Higashi, Tsukuba 305-8567, Japan d
Methane Hydrate Project Unit, Research Institute of Energy Frontier, Department of Energy
and Environment, National Institute of Advanced Industrial Science and Technology (AIST), 2e
TE D
17 Tsukisamu-Higashi, Toyohira, Sapporo 062-8517, Japan
Methane Hydrate Project Unit, Research Institute of Energy Frontier, Department of Energy
and Environment, National Institute of Advanced Industrial Science and Technology (AIST), Tsukuba West, 16-1 Onogawa, Tsukuba 305-8569, Japan Department of Energy and Environment, National Institute of Advanced Industrial Science and
EP
f
AC C
Technology (AIST), 2-17 Tsukisamu-Higashi, Toyohira, Sapporo 062-8517, Japan
Corresponding author: T. Sato Tel: +81-4-7136-4726, Fax: +81-4-7136-4726, E-mail:[email protected] 1
ACCEPTED MANUSCRIPT
ABSTRACT
Commercial-scale extraction of methane gas from subsea methane hydrate is difficult. A major challenge is mud erosion at the interfaces of sand–mud alternate layers in the subsea strata.
RI PT
Without a clear understanding of the various parameters associated with mud erosion, there may be complications in developing methane hydrate deposits, such as reservoir clogging due to mud particles. This study aims to develop a numerical simulation method to predict mud erosion caused by water flow in pore-scale computational domain. Microscopic sand beds were
SC
generated numerically in the computational domain while the water flow through the pores of various-shaped sand grains was numerically simulated using the lattice Boltzmann method. Mud
M AN U
was treated as a continuum and eroded at its surface via the shear stress of water flow using the critical shear stress and erosion rate constant, which were estimated by fitting calculated erosion rates with those of a laboratory-scale experiment. Using this numerical model, simulations of long-term mud erosion were conducted. The simulation showed that when the average velocity of water in the pore space is larger than the critical value and smaller than ten times that value,
TE D
the erosion rate is almost zero after about 10 days.
KEYWORDS: Mud erosion; Sand–mud alternate layer; Lattice Boltzmann method; Subsea
AC C
EP
sediment; Critical shear stress; Erosion rate.
2
ACCEPTED MANUSCRIPT
1. INTRODUCTION Mud erosion by water and/or gas flows at the interface of sand–mud alternate layers may cause problems for subsea engineering processes such as oil and gas production. In particular, the production of methane gas by depressurization from methane hydrate, which is found in sub-
RI PT
seabed sand layers that alternate with mud layers in turbidite formations (Fujii et al., 2015), causes a two-phase flow through the sand layers that comprises dissociated gas and water. This can potentially result in mud erosion of and, possibly, blockage of the flow.
Several types of mud erosion for fine-grained sediments were examined by Amos et al.
SC
(1992). Winterwerp et al. (2004) classified the erosion models for cohesive sediments: floc erosion, entrainment, surface erosion, and mass erosion. In terms of erosion of sand–mud
M AN U
mixtures, Mitchener et al. (1996) stated that the critical shear stress, τc, increases when sand is added to mud, where τc is the smallest local shear stress to cause erosion. Moreover, they showed that τc for sand–mud mixtures with sand contents of 70%–100% depends on the size of the sand grains and on the cohesive properties of the mud materials. Detailed studies of erosion types for sand–mud mixtures were conducted by Jacobs et al. (2011). After investigating different soil samples with varying compositions, they categorized erosion of sand–mud mixtures into two
TE D
types: time-decreasing and time-independent erosion.
The above studies are mainly limited to mud erosion on sea floors or river bottoms. For mud erosion at the interface of sand–mud alternate layers, Brumby et al. (2015) conducted numerical simulations using a combination of the lattice Boltzmann method (LBM) (based on
EP
the approaches of Sugita et al. (2012) and Guo et al. (2013)) and the discrete element method (DEM) developed by Cundall et al. (1979). Based on erosion studies on cohesive sediments by
AC C
Partheniades (1965), Oyama et al. (2016) introduced an equation to express mud erosion rate:
(
α exp U − U c E= 0
)
U ≥ Uc
U < Uc
,
(1)
where E [kg/m2/s] is the erosion rate, αexp [kg/m3] is the erosion rate constant, U [m/s] is the average velocity of water, and Uc [m/s] is the critical velocity, which is the erosion threshold. The reason why the average water velocity was chosen over the Darcy velocity is that mud erosion is directly relevant to the shear stress on the mud surface caused by water flow. Thus, we 3
ACCEPTED MANUSCRIPT
decided that the average water velocity would be more appropriate to represent the shear stress than the Darcy velocity, which is affected by porosity. The erosion rate constant and critical velocity were estimated by Oyama et al. (2016), who performed experiments of mud erosion using artificial sand–mud alternate cores.
RI PT
The primary objective of this study was to develop a numerical simulation method to realize mud erosion in sand–mud alternate layers. We demonstrated pore-scale numerical
simulations of water-induced erosion using the LBM. First, we generated a numerical pore-scale domain that mimics a structure composed of several frame-sand grains. At the bottom of this
SC
microscopic computational domain, we filled the pore space with mud. Next, water flow through this porous domain was analyzed using the LBM. Then, based on the resultant water velocity, we
M AN U
calculated the shear stress acting on the surface of the mud layer and changed the shapes of the mud surface to represent the mud erosion in the LBM.
Although of great interest, the temporal change of the erosion rate over long periods of time was previously unknown; however, the results of this study’s long-time simulations show how the erosion rate changes with time. Therefore, the results of our numerical simulation may be able to offer important insights that will be valuable for the gas production process from
TE D
methane hydrate reservoirs.
2. NUMERICAL METHOD
2.1 Generation of the sand bed
EP
Sugita et al. (2012) extracted the shapes of the sand grains using computed tomography images of sand cores obtained from hydrate sedimentary layers off Japan’s coast. Thereafter, they
AC C
packed these grains in a pore-scale computational domain using a numerical growth method to attain the required porosity. Figure 1 illustrates an example of a generated sand bed with six grains, which has a porosity of 0.40 under a periodic boundary condition. Because of the limitations of our computational facilities, the length of each side of the domain was set to be 200 µm. Because the domain is so small that it cannot represent the characteristics of the whole sand layer, we generated 10 domains and took averages of the properties obtained from the calculations (details presented later in section 2.6).
4
RI PT
ACCEPTED MANUSCRIPT
Fig. 1 An example of a three-dimensional computational domain holding six sand grains with a
SC
periodic boundary condition.
2.2 Lattice Boltzmann method (LBM)
M AN U
Porous substances are structurally complicated, making the creation of a computational matrix difficult with conventional computational fluid dynamics (CFD) methods. The LBM is a powerful tool for the calculation of complex liquid flows with low Reynolds number. The principal advantages of using LBM for the simulation of flow through porous media have been confirmed by many researchers, e.g., Guo et al. (2002) confirmed the suitability of the LBM for the simulation of incompressible flow through porous media. The LBM has additionally been
TE D
used for a number of other applications, e.g., a study conducted by Jain et al. (2016) explored the conditions for transition to turbulence in intracranial aneurysms at spatial and temporal resolutions near the Kolmogorov microscales. LBMs for microfluidics were well reviewed by Zhang (2011). Sugita et al. (2012) also developed a simulation code to calculate water flows
(2012).
EP
within pore-scale porous media. In this study, we applied the LBM developed by Sugita et al.
AC C
For solid boundaries with curved surfaces, we adopted the boundary condition proposed
by Bouzidi et al. (2001), which can compensate for the staircase approximation in the LBM. There are various models for discretizing velocity, and in this study, we adopted the threedimensional (3D) fifteen-velocity model developed by Qian (1992) (D3Q15). Nash et al. (2014) confirmed the usage of this model in complex flow domains, and the method of Sugita et al. (2012) was very stable when using the D3Q15 stencils together with the Bhatnagar–Gross– Krook (BGK) model for the collision term (Bhatnagar et al., 1954). Using the BGK model for
5
ACCEPTED MANUSCRIPT
the collision term, we can express the particle distribution function fi with the discrete Boltzmann equation (MacNamara, 1988): ∂fi 1 + ci ⋅ ∇f i = − (f i − f ieq ), ∂t τ
RI PT
(2)
where ci indicates the virtual particle velocity in the direction towards i, which denotes each of the surrounding 14 lattice points, t is time, and τ is the relaxation time. The local equilibrium
is the general kinetic theory of fluids expressed as
c=
M AN U
3 9 3 2 f ieq = ρw i 1 + 2 c i ⋅ u + 4 (c i ⋅ u ) − 2 u ⋅ u , 2c 2c c
SC
distribution function fieq is an extension of the Maxwell distribution to a local equilibrium, which
∆x , ∆t
(3)
(4)
TE D
where ∆t [s] is the time step, ∆x [m] is the lattice unit scale, and wi is the weighting factor depending on the discrete velocity, ci. The fluid density ρ [kg/m3] and fluid velocity u [m/s] are
ρ = ∑f i , i
i
(5) (6)
AC C
u = ∑cifi .
EP
given by
To correspond to the Navier–Stokes equation, the relaxation time τ is set to satisfy the
following equation: 1 ν = τ − c s2 ∆t, 2
(7)
6
ACCEPTED MANUSCRIPT
where ν [m2/s] is the kinematic viscosity of water. Calculations are performed separately for the collision and streaming processes using the
1 f i ( x + c i ∆t, t + ∆t ) = fi ( x, t ) − [fi ( x, t ) − f ieq ( x, t )], τ f i ( x + c i ∆t , t + ∆t ) = fi (x, t ).
RI PT
following equations:
(8)
(9)
SC
2.3 Calculation of shear stress
Forces acting on the surface of the continuum mud can be calculated via the momentum exchange method (e.g., Mei et al., 2002; Ladd, 1994), which is a simple implementation for fluid
M AN U
flows, assuming that the surface of the continuum is a solid boundary. The force originating from momentum exchange per unit time between neighboring lattices is expressed as follows.
Fi =
1 [c i f i (x) − c −i f −i (x + c i ∆t) ], ∆t
(10)
TE D
where F [kg m/s2] is the force vector acting on the mud surface and −i denotes the reflected direction, which is the bounce-back direction, as opposed to the direction of i. To evaluate the shear stress on a corner mud lattice, which has more than one side exposed to water lattices, F on
EP
the corner mud lattice is calculated using (10) as the summation of Fi, which is toward the water lattices. The tangential force on the mud surface can be expressed using the following equation:
n=−
AC C
Ft = F − Fn = F − F ⋅ n ,
∑Q c
i i
i
2
,
(11)
(12)
∑Q c
i i
i
where Ft and Fn are the force components of the tangential and normal direction of the mud surface, respectively, n is the normal unit vector, and Qi is the amount of mud in a lattice cell 7
ACCEPTED MANUSCRIPT
with direction i. The shear stress can be obtained by dividing the tangential force with the area of the mud surface.
2.4 Validation of the shear stress calculations
RI PT
To determine the validity of the calculated shear stress values, numerical simulations were conducted for the flow between parallel flat plates, and the calculated shear stresses were
compared with a theoretical value of the Poiseuille flow. The exact value of the shear stress can
1 dP h, 2 dx
M AN U
τs =
SC
be expressed via the following equation:
(13)
where dP/dx [Pa/m] is the pressure gradient and h [m] is the distance between the two parallel plates. Meanwhile, the exact average velocity Ua [m/s] is given by
1 dP 2 h , 12 µ dx
(14)
TE D
Ua =
where µ [kg/m/s] is the viscosity of water.
The simulation conditions were as follows. The flat plate is shown in black in Fig. 2. The
EP
domain sizes were 50, 50, and 200 µm in the x, y, and z directions, respectively. A periodic boundary condition was imposed for the side boundaries. The surface of the flat plate had a noslip condition. The upper surface of the domain was a symmetric boundary so that the height of
AC C
the domain was equal to h/2. A pressure gradient of 1000 Pa/m was imposed in the x-direction in the whole domain.
Here, the scale of the computational lattices and the relaxation time were examined by
comparing them to the theoretical values calculated using (13) and (14). Figures 3 and 4 show the comparison with the average velocity and wall shear stress of the Poiseuille flow with various relaxation times and computational lattice sizes. According to Kruger et al. (2009), their LBM with a relaxation time of 0.9 provides small errors for the velocity and shear stress. Kruger et al. (2009) demonstrated the cases of Re = 1, 10, and 100, and Ma = 0.01, 0.04, and 0.20. The
8
ACCEPTED MANUSCRIPT
number of the lattice cells in their computations ranged from 16 × 18 × 18 to 128 × 130 × 130, and the total computational time was between 0.25 and 0.75 s. In our numerical simulations, Re and Ma were approximately 2 and 0.03. The number of the lattice cells in our computations
and 1.5 s. The water velocity was set to 0 as the initial condition.
RI PT
ranged from 10 × 10 × 40 to 25 × 25 × 100, and the total computational time was between 0.13
Figures 3 and 4 agree well with the suggestion by Kruger et al. (2009), although all the simulated results are in moderately good agreement with the theory. In this study, a relaxation time of 1.0 and a computational lattice with dimensions of 2.5 µm were adopted for the
EP
TE D
M AN U
shear stress can be computed with an error margin of <1%.
SC
following numerical simulation. Using these parameters, the average water velocity and the wall
Fig. 2 Computational domain for the Poiseuille flow. The lower part of a parallel flat plate is
AC C
shown in black. The upper surface is a symmetric boundary and a periodic condition is imposed on the side surfaces.
9
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
Fig. 3 Calculated average velocity compared with the theoretical values using different
AC C
EP
TE D
relaxation times and lattice sizes.
Fig. 4 Calculated wall shear stress compared with theoretical values using different relaxation times and lattice sizes. 10
ACCEPTED MANUSCRIPT
2.5 Numerical model of mud erosion Following the erosion-rate model, (1) used by Oyama et al. (2016), we also adopted a linear
α cal (τ s − τ c ) E= 0
τs ≥ τc
,
(15)
SC
τs < τ c
RI PT
relationship between the local erosion rate and the shear stress, given by
where E [kg/m2/s] is the local erosion rate, αcal [s/m] is the erosion rate constant, τs [Pa] is the
M AN U
shear stress, and τc [Pa] is the critical shear stress being applied on the mud. Mud erodes when the shear stress acting on its surface is larger than the critical value.
To match (15) with (1), the average erosion rate, E [kg/m2/s], in the numerical simulation is given by
1 E dA, A τ s ∫≥ τ c
(16)
TE D
E=
where A is the total area of the sand–mud interface, including the sectional areas of the sand grains that penetrate the interface. It is noted that (1) and (15) do not hold at the same time
EP
because the integration in (16) is only done when τ ≥ τ and such area varies depending on the average water velocity. Therefore, the linearity in (15) never assures that in (1). However, (1) is
AC C
also meaningful from an engineering viewpoint because it can be easily applied to actual gas or oil production operations.
When the volume of mud in a lattice cell is larger than half the volume of a lattice cell
before erosion and becomes smaller than half after the erosion, we consider that the lattice originally defined as a mud lattice has become a water lattice. In such a case, the distribution function of this lattice is set to 0, and consequently, the water velocity of this lattice becomes 0 because the mud was not originally moving before the erosion. Thus, local cell-wise mass conservation does not hold with this numerical treatment. However, the mud volume eroded is assumed to be very small, <0.1% of the initial volume after long-time simulation as shown later, 11
ACCEPTED MANUSCRIPT
and we consider that it is negligible.
2.6 Estimation of critical shear stress and erosion rate constant The critical shear stress and erosion rate constant in the numerical simulation were estimated
RI PT
based on experimental data. Oyama et al. (2016) conducted laboratory-scale experiments on mud erosion under varying water velocities. To create an artificial core that mimics sand–mud
alternate layers, they packed kaolinite as mud particles and Toyoura sand as large frame grains in each half of a cylinder of length 100 mm and an inner diameter of 49.8 mm. Water was injected
SC
into the core, and the discharged water containing eroded mud was sampled at set time intervals. After analyzing the results of the experiment, Oyama et al. (2016) obtained a critical velocity, Uc,
M AN U
of 1.81 × 10−5 m/s and an experimental erosion rate constant, αexp, of 0.00374 kg/m3 at a confining pressure of 4 MPa.
Although measuring shear stress directly on the mud surface is difficult, we need to input the critical shear stress in our mud erosion simulation. Here, numerical simulations were conducted under the condition that the average water velocity is set to the critical velocity, Uc, by adjusting the pressure gradient to identify the critical shear stress, τc. The computational
TE D
conditions are listed in Table 1. The mud layers were set at the bottom of the computational domains, as shown in Fig. 5. In these simulations, the mud was not eroded, and it was treated as a fixed solid with its surface having no no-slip condition. A periodic boundary condition was applied to the side surfaces. The upper surface is a boundary with a gradient of zero in the
EP
vertical direction.
AC C
Table 1. Computational condition
Scale of calculation area
Number of lattices Porosity
200 × 10−6 [m] 80 × 80 × 80 0.4
Lattice size
2.5 × 10−6 [m]
Time step
1.04 × 10−6 [s]
12
SC
RI PT
ACCEPTED MANUSCRIPT
M AN U
Fig. 5 Isometric view of a computational domain at the initial stage. Sand, mud, and water are shown as gray, black, and white, respectively. Mud is set in the pore space close to the bottom of the domain.
According to Sugita et al. (2012), the Darcy velocity was validated when 30 lattices were used for a glass-bead particle with a diameter of 100 µm. The sizes of the computational domain
TE D
and the sand grains used by Sugita et al. (2012) are identical to those used in the present simulation, and therefore, we consider that the lattice size is sufficiently small for an accurate analysis.
EP
Because our computational domain includes only six sand grains and pore-scale mud erosion depends on the shapes and locations of the grains, we conducted flow simulations using 10 different domains with different sand-packing patterns and took average of the resultant data
AC C
of the 10 calculations to determine meaningful parameters. Figure 6 shows the cumulative frequency of shear stress on the mud surface of the 10
computational domains. As the average water velocity was set to be Uc, there should be no mud erosion and the shear stress should be smaller than τc. Here, we determined τc to be the value at which the cumulative area of the mud surface is 99.73%, corresponding to the 3σ upper limit of the Gauss distribution. Eventually, τc is 0.0233 Pa, as shown in Fig. 6.
13
100
99.73 ％
80 60 40 20
0.0233 Pa
0 0
0.005
0.01 0.015 0.02 Shear stress [Pa]
0.025
RI PT
Cumulative frequency [%]
ACCEPTED MANUSCRIPT
0.03
SC
Fig. 6 Cumulative frequency of shear stress calculated using 10 sand patterns under the condition
M AN U
that the average water velocity is the same as the critical velocity.
Next, we examined the effect of the mud phase height in the computational domains, as shown in Fig. 7 (sand, mud, and water are displayed in gray, black, and white, respectively). Here, we imposed a pressure gradient of 1000 Pa. Figure 8 shows the relationship between the average of the calculated shear stress acting on the flat mud surface and the mud heights under the same pressure gradient; the mud height does not exceed 40 lattice cells and the shear stresses
TE D
are almost constant. It is therefore fair to say that the initial height of mud being two lattice cells
AC C
EP
in this mud erosion simulation is validated.
(a)
(b)
14
RI PT
ACCEPTED MANUSCRIPT
(c)
(d)
SC
Fig. 7 Examples of computational domains with different heights of mud phase. The sand, mud, and water are displayed in gray, black, and white, respectively. The number of mud
EP
TE D
M AN U
lattices are (a) 2, (b) 20, (c) 40, and (d) 60 in the vertical direction.
Fig. 8 Calculated shear stress with various mud lattice heights in the vertical direction, as shown
AC C
in Fig. 7. The plots are the average using 10 different sand patterns.
Finally, the erosion rate constant was estimated by matching the erosion rate with the
experimental data for kaolinite content of 100%. The average velocities of water in the numerical simulation were set to be the same as those of the experiment, i.e., water injection rates of 0.5, 0.7, 1, 2, and 3 g/min, by adjusting the pressure gradient. The simulations were conducted for each of these water velocities using 10 sand patterns. The errors in the erosion rate between the numerical simulation and the experiment were minimized by adjusting the erosion rate constant using the least squares method. The best-fit erosion rate constant, αcal, was 7.47 × 10−5 s/m. The
15
ACCEPTED MANUSCRIPT
erosion rates calculated using this value are shown in Fig. 9. Although the calculated and measured erosion rates are moderately identical, the general trend is that the calculated rates are slightly smaller than those in the experiments. This may be acceptable because the erosion rate
M AN U
SC
RI PT
totally depends on the sand grain arrangement in the domain.
Fig. 9 Comparison of the calculated erosion rates with the experimental values obtained by Oyama et al. (2016).
TE D
3. TEMPORAL CHANGE OF THE EROSION RATE OVER LONG TIME Now, using the values of αcal and τs obtained above, a decrease in the mud volume can be calculated for long time periods, e.g., several weeks. The computational conditions for long-term mud erosion are shown in Table 2.
EP
To conduct a long-term simulation after the mud surface reaches the bottom of the domain, we move the top lattice layer below the computational domain by filling mud into its
AC C
whole pore whenever a part of the eroded mud surface touches the bottom of the computational domain. This treatment poses no difficulty because the sand grains are packed in a periodic way in the vertical direction.
In addition, we enlarged the erosion rate constant, αcal, by 107 times unless the volume of
eroded mud per time step was still less than that of the computational cell. This is because too much time is needed to calculate erosion for a long time range and because the physics is theoretically the same if αcal × t is the same. As shown in Fig. 10, the results using one sand pattern indicate that using magnifications of 106 and 107 produce almost identical results. In contrast, using magnifications of 108 and 109 16
ACCEPTED MANUSCRIPT
produced different results compared with using magnifications of 106 and 107. As a result, using a larger magnification can result in an incorrect erosion rate numerically, although theoretically, the magnification should not have an effect on the result. This may be because if the volume of the mud eroded in one time step is as large as the lattice cell volume, the eroded volume in a cell
RI PT
many remain cubic-cell-shaped, which is an artifact. Based on the results in Fig. 10, we adopted
M AN U
SC
a magnification of 107 to save computational time.
Fig. 10 Temporal changes of the calculated erosion rate in 3 days using one sand pattern
TE D
with varying magnifications for the erosion rate constant, αcal: 106, 107, 108, and 109. Table 2. Numerical simulation conditions
EP
Scale of calculation area
200 × 10−6 [m]
Number of lattices
80 × 80 × 80
Porosity
0.4 2.5 × 10−6 [m]
Time step
1.04 × 10−6 [s]
Critical shear stress
0.0233 [Pa]
Erosion rate constant
7.47 × 10−5 [s/m]
Bulk density of mud (porosity = 0.4)
1110 [kg/m3]
Average velocity of water
2.5Uc and 7.5Uc
Total time into real scale
15 [days]
AC C
Lattice size
17
ACCEPTED MANUSCRIPT
Figure 11 shows the temporal change of the erosion rates, which are average values of the results using 10 sand patterns, for one day. In this and subsequent figures, the erosion rate and time are modified such that they are small and long by dividing and multiplying by the expansion ratio of αcal, 107, respectively. The erosion rates did not show significant changes with time, at
RI PT
least within several hours. Oyama et al. (2016) assumed that the erosion rate was constant during their experiment, which lasted up to several hours, to extract important erosion parameters.
(a) Average water velocity was 2.5Uc
M AN U
SC
Evidently, the present numerical simulation supports their assumption.
(b) Average water velocity was 7.5Uc
Fig. 11 Temporal changes of the calculated erosion rates in 1 day. The results are the average
TE D
values of 10 different sand patterns. The average water velocities were 2.5Uc and 7.5Uc, where Uc represents the critical velocity that induces mud erosion. Temporal changes in the calculated erosion rates, which are average values of the results
EP
using 10 sand patterns, over long periods of time are shown in Fig. 12. The erosion rate decreases monotonously with time with small fluctuations (almost zero). These fluctuations seem
AC C
to depend on the morphology of the sand and mud. To understand this temporal change in the erosion rate, the velocity vectors and mud locations, shown in black, are depicted in Fig. 13 on the 1st and 10th days on a vertical section of one of the 10 domains. The velocity vectors are drawn at every 3 computational lattices. On the 1st day, the mud was continuously eroded around point A in Fig. 13(a). On the 10th day, the water path around point A in Fig. 13(b) was broadened because of erosion, and the flow velocities right above the eroded mud surface decreased. When spaces between the sand grains and the mud surfaces became sufficiently large, the shear stress was smaller than τc along such channels, and there was no further erosion (Fig. 13(b)). In contrast, a small amount of the mud was eroded around point B in Fig. 13(a) and (b). 18
ACCEPTED MANUSCRIPT
This may be because point B is located in the dead wake of a sand grain and thus there is no
M AN U
SC
RI PT
potential for erosion. This is another mechanism through which erosion can cease.
EP
TE D
(a) Average water velocity was 2.5Uc
(b) Average water velocity was 7.5Uc
AC C
Fig. 12 Temporal changes of the calculated erosion rates during 15 days. The same graph with the time period ranging from the 5th to 15th day is enlarged (inset).
19
RI PT
ACCEPTED MANUSCRIPT
(b) 10th day
SC
(a) 1st day
Fig. 13 Vertical sections of the computational domain on the 1st and 10th day for one sand
M AN U
pattern. Sand, mud, and water are shown as gray, black, and white, respectively. The arrows indicate the velocity vectors of water flow [m/s].
Figure 14 shows the histograms of the calculated shear stress distribution when the average water velocity was 7.5Uc on the 0th and 15th days. The result shows that shear stresses larger than τc, 0.0233 Pa, disappeared by the 15th day. The decrease in the erosion rate may be
TE D
attributed to the water flow eroding the mud so that smooth flow channels are formed near the
AC C
EP
deformed mud surface, as can be seen in Fig. 13.
Fig. 14 Histogram of shear stress acting on mud surface in the case where average water velocity is 7.5Uc. The solid and dotted lines show the results of 0th and 15th day, respectively.
Overall, when the average velocity of water, U, is several times larger than the critical 20
ACCEPTED MANUSCRIPT
value, our results suggested that mud erosion occurs only at an early stage of the gas production process from methane–hydrate bearing layers (approximately within 10 days). Oyama et al. (2016) reported that microfractures are formed in sand layers when the average velocity is larger than 2.14 × 10−4 m/s, which is approximately 10Uc. According to Oyama et al. (2010), the water
RI PT
velocity should be less than approximately 10−4 m/s, even in the vicinity of the wellbore;
otherwise, the sand frame is broken and sand is produced. Here, let us assume that the radius of a wellbore is R and water velocity around it is 2 × 10−4 m/s, approximately 10Uc. We can then estimate water velocity at 10R from the wellbore as 2 × 10−5 m/s, approximately Uc. Thus, mud
SC
erosion is considered to occur within the circle with a radius of 10R over a time period of <10 days. Assuming that the radius of a production well is 0.15 m, the cumulative mass of the eroded
M AN U
mud in a radius of 10R within 10 days can be estimated approximately as less than 1 kg at an interface of a sand–mud alternate layer. As a result, this estimation suggests that the total mass of the eroded mud is trivial even if there are multiple sand–mud layers. In this study, the mud was treated as a continuum solid that changes its shape because of erosion using (15), and it is assumed that the eroded mud does not affect any flow properties because the concentrations of mud particles are supposedly very less.
TE D
The present numerical simulation has some drawbacks, namely, sand grains do not move after the mud bases have eroded. A possible reason for this is that as in Fig. 13, the eroded mud in the figure can be seen as approximately 50 µm in depth. Similarly, in most cases, the mud surfaces were eroded 50–100 µm in depth before the erosion stopped. This result suggests that
EP
the thickness of the eroded mud within 10 days is smaller than the size of the sand grains. In addition, the sands originally formed a frame and do not necessarily move. In future, we intend
AC C
to demonstrate the decomposition of frame-sand grains as well as the formation of microfractures in sand layers by incorporating the sand grains into the present mud erosion simulation.
4. CONCLUSIONS
A numerical simulation method using LBM was developed to estimate the mass of mud erosion at the interface of sand–mud alternate layers. The sand beds were numerically generated in microscopic computational domains and water flow was forced within the pore space by a constant pressure gradient. First, we confirmed the accuracy of this numerical method by 21
ACCEPTED MANUSCRIPT
accurately comparing the calculated wall shear stress with that of the Poiseuille flow. Then, mud erosion was simulated for the generated microscopic sand–mud layers using experimentally obtained erosion parameters. We studied the validity of the present method by comparing the calculated erosion rate with that of the experimentally obtained value.
RI PT
Next, we applied the present method to a long-term erosion simulation. The calculated results suggest that the erosion rate decreases to zero in less than 10 days under the condition that the average velocity is several times larger than the critical value for mud erosion. This may be because the mud surface with large shear stress is eroded first, resulting in moderately wide
SC
spaces appearing in the pore of the sand grains, which slows the water velocity to a value smaller
AC C
EP
TE D
M AN U
than the critical value.
22
ACCEPTED MANUSCRIPT
ACKNOWLEDGEMENTS
This study was financially supported by the Research Consortium for Methane Hydrate Resources in Japan (MH21 Research Consortium) of Methane Hydrate Research and
AC C
EP
TE D
M AN U
SC
RI PT
Development Program planned by the Ministry of Economy, Trade and Industry (METI).
23
ACCEPTED MANUSCRIPT
REFERENCES
Amos, C.L., Daborn, G.R., Christian, H.A., Atkinson, A., Robertson, A., 1992. In situ erosion
RI PT
measurements on fine-grained sediments from the Bay of Fundy. Mar. Geol. 108, 175-196.
Bhatnagar, P.L., Gross, E.P., Krook, M., 1954. A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems. Phys. Rev. 94, 511-525.
M AN U
fluid with boundaries. Phys. Fluids 13, 3452-3459.
SC
Bouzidi, M., Firdaouss, M., Lallemand, P., 2001. Momentum transfer of a Boltzmann-lattice
Brumby, P. E., Sato, T., Nagao, J., Tenma, N., Narita, H., 2015. Coupled LBM DEM micro-scale simulations of cohesive particle erosion due to shear flows. Transp. Porous Med. 109, 43-60.
Cundall, P.A., Strack, O.D.L., 1979. A discrete numerical model for granular assemblies.
TE D
Géotechnique 29, 47–65.
Fujii, T., Suzuki, K., Takayama, T., Tamaki, M., Komatsu, Y., Konno, Y., Yoneda, J., Yamamoto, K., Nagao, J., 2015. Geological setting and characterization of a methane hydrate reservoir distributed at the first offshore production test site on the daini-atsumi knoll in the
EP
eastern nankai trough. Japan. Mar. Petrol. Geol. 66, 310-322.
AC C
Guo, Z., Zhao, T. S., 2002. Lattice Boltzmann model for incompressible flows through porous media. Phys. Rev. E 66, 036304.
Guo, Z., Shu, C., 2013. Lattice Boltzmann method and its applications in engineering.World Scientific Publishing Co. Pte. Ltd., Singapore.
Jacobs, W., Hir, P. L., Kesteren, W. V., Cann, P., 2011. Erosion threshold of sand/mud mixtures. Continental Shelf Res. 31, 14-25.
24
ACCEPTED MANUSCRIPT
Jain, K., Roller, S., Mardal, K. A., 2016. Transitional flow in intracranial aneurysms - A space and time refinement study below the Kolmogorov scales using Lattice Boltzmann Method. Comput. Fluids 127, 36-46.
RI PT
Kruger, T., Varnik, F., Raabe, D., 2009. Shear stress in lattice Boltzmann simulations, Phys. Rev. E 79, 046704.
Ladd, A., 1994. Numerical simulations of particulate suspensions via a discretized Boltzmann
SC
equation. Part 1. Theoretical foundation, J. Fluid Mech. 271, 285-309.
M AN U
McNamara, G., Zanetti, G., 1988. Use of the Boltzmann equation to simulate lattice-gas automata. Phys. Rev. Lett. 61, 2332-2335.
Mei, R., Yu, D., Shyy, W., Luo, L-S., 2002. Force evaluation in the lattice Boltzmann method involving curved geometry. Phys. Rev. E 65, 041203.
TE D
Mitchener, H., Torfs, H., 1996. Erosion of mud/sand mixtures. Coast. Eng. 29, 1-25.
Nash, R. W., Carver, H. B., Bernabeu, M. O., Hetherington, J., Groen, D., Krger, T., Coveney, P. V., 2014. Choice of boundary condition for lattice-Boltzmann simulation of moderate-Reynolds-
EP
number flow in complex domains. Phys. Rev. E 89, 023303.
AC C
Oyama, H., Abe, S., Yoshida, T., Sato, T., Nagao, J., Tenma, N., Narita, H., 2016. Experimental study of mud erosion at the interface of artificial sand-mud alternate layer. J. Nat. Gas Sci. Eng. 34, 1106–1114.
Oyama, H., Nagao, J., Suzuki, K., Narita, H., 2010. Experimental analysis of sand production from methane hydrate bearing sediments applying depressurization method. J MMIJ 126, 497– 502 (in Japanese).
Partheniades, E., 1965, Erosion and deposition of cohesive soils. J. Hydraul. Div. 91, 105-139. 25
ACCEPTED MANUSCRIPT
Qian, Y. H., D’Humires, D., Lallemand, P., 1992. Lattice BGK models for the Navier-Stokes equation. Europhys. Lett. 17, 479-484.
RI PT
Sugita, T., Sato, T., Hirabayashi, S., Nagao, J., Jin, Y., Kiyono, F., Ebinuma, T., Narita, H., 2012. A pore-scale numerical simulation method for estimating the permeability of sand sediment. Transp. Porous Med. 94, 1-17.
SC
Winterwerp J. C., van Kesteren, W.G.M., 2004. Introduction to the physics of cohesive sediment
M AN U
in the marine environment. Chapter 9, Elsevier, Amsterdam.
Zhang, J., 2011. Lattice Boltzmann method for microfluidics: models and applications.
AC C
EP
TE D
Microfluid Nanofluid 10, 1-28.
26
ACCEPTED MANUSCRIPT
Highlights
A pore-scale numerical simulation method for mud erosion was developed.
•
Microscopic sand beds are generated numerically in the computational domain.
•
Water flow through the sand pore is simulated using a lattice Boltzmann method.
•
Numerically, mud is eroded by the water shear stress using erosion rate constant.
•
Long-term simulation suggests that erosion rate approaches zero in about 10 days.
AC C
EP
TE D
M AN U
SC
RI PT
•