- Email: [email protected]
Geomechanical constitutive modelling of gas hydrate-bearing sediments by a state-dependent multishear bounding surface model
Journal Pre-proof Geomechanical constitutive modelling of gas hydrate-bearing sediments by a statedependent multishear bounding surface model Huolang Fang, Kenan Shi, Yang Yu PII:
S1875-5100(19)30371-3
DOI:
https://doi.org/10.1016/j.jngse.2019.103119
Reference:
JNGSE 103119
To appear in:
Journal of Natural Gas Science and Engineering
Received Date: 14 September 2019 Revised Date:
15 November 2019
Accepted Date: 18 December 2019
Please cite this article as: Fang, H., Shi, K., Yu, Y., Geomechanical constitutive modelling of gas hydrate-bearing sediments by a state-dependent multishear bounding surface model, Journal of Natural Gas Science & Engineering, https://doi.org/10.1016/j.jngse.2019.103119. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. 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. © 2019 Published by Elsevier B.V.
1
Title of paper
2
Geomechanical constitutive modelling of gas hydrate-bearing sediments by a state-dependent
3
multishear bounding surface model
4 5
Author 1
6
Huolang Fang
7
College of Civil Engineering and Architecture, Zhejiang University, 866 Yuhangtang Road,
8
Hangzhou 310058, China
9
Email: [email protected]
10 11
Author 2
12
Kenan Shi
13
College of Civil Engineering and Architecture, Zhejiang University, 866 Yuhangtang Road,
14
Hangzhou 310058, China
15
Email: [email protected]
16 17
Author 3
18
Yang Yu (Corresponding author)
19
Ocean College, Zhejiang University, 1 Zheda Road, Dinghai District, Zhoushan, Zhejiang 316021，
20
China
21
Email: [email protected]
1
22
Abstract
23
Natural gas hydrate in marine sediments and permafrost areas is considered as an important
24
potential energy source. Since hydrate dissociation will reduce the stability of gas hydrate-bearing
25
sediments (GHBS) and may cause wellbore failures and geological disasters during gas production,
26
it is necessary to reveal the mechanical behavior of GHBS for the safe exploitation of natural gas
27
hydrate. This paper proposes a geomechanical constitutive model of GHBS within the multishear
28
bounding surface framework. Following the slip theory of plasticity, a constitutive formulation is
29
obtained by splitting the macro constitutive response of sediments into a macro volume response
30
and a series of micro shear responses in spatial distributions related to virtual microshear structures.
31
Each microshear structure describes micro shear and dilatancy responses in three orthogonal
32
orientations. A micro stress–strain relationship and a micro stress–dilatancy relationship are
33
established for each orientation of the microshear structure. The model comprehensively describes
34
the consolidation, hardening, softening, dilatation, collapse, and non-coaxial characteristics of gas
35
hydrate-bearing sediments by introducing the multishear concept, state parameter, evolution law of
36
hydrate bonding and debonding, and collapse strain caused by hydrate dissociation. The
37
effectiveness of the model is confirmed by simulating the available published laboratory tests on
38
the samples of synthetic and natural GHBS under different pore pressures, temperatures, initial void
39
ratios, hydrate saturations, and initial effective confining stresses.
40
Keywords
41
gas hydrate, mechanical behavior, constitutive model, dilatancy, collapse deformation
2
42
1. Introduction
43
Natural gas hydrate is the crystalline solid where gas molecules are stored in ice-like crystal
44
lattices and is stable under high pressure and low temperature in marine sediments and permafrost
45
areas. The amount of natural gas trapped inside hydrate is roughly evaluated to be twice the existing
46
reserves of conventional fossil fuels in the earth. As a significant potential energy resource, the
47
exploitation of gas hydrate has attracted wide attention. At present, three possible ways to produce
48
gas from gas hydrate-bearing sediments (GHBS) are decompression, heat injection, and chemical
49
stimulation. Gas hydrate is usually located in unconsolidated sediments with low shear strength.
50
Hydrate dissociation will reduce the stability of GHBS due to hydrate exploitation. Particularly,
51
hydrate dissociation in the deep-sea marine sediments may result in wellbore failures and geological
52
disasters, such as casing damage, seafloor deformation, and submarine landslide. Therefore, it is
53
necessary to reveal the mechanical behavior of GHBS for the safe exploitation of natural gas
54
hydrate.
55
Over the past two decades, many experimental researches in the mechanical behavior of synthetic
56
GHBS have been performed (Ebinuma et al., 2005; Masui et al., 2005; Priest et al., 2005, 2009;
57
Clayton et al., 2005, 2010; Yun et al., 2007; Miyazaki et al., 2011; Ghiassian and Grozic, 2013;
58
Hyodo et al., 2013a, 2013b, 2014, 2017; Santamarina et al., 2015; Liu et al., 2016; Yoneda et al.,
59
2016; Kajiyama et al., 2017a, 2017b; Liu et al., 2017; Sultaniya et al., 2018; and Choi et al., 2018).
60
For example, Masui et al. (2005) performed several triaxial compression tests for samples of
61
synthetic GHBS involving different hydrate saturations, and concluded that the stiffness and
62
strength of samples were proportional to the saturation of hydrate. Priest et al. (2005, 2009) studied
63
the impacts of hydrate formation on the wave velocity in samples, revealing that the contact of
64
hydrate-cemented particles significantly enhanced the wave velocity. Miyazaki et al. (2011)
65
performed a set of creep tests in the triaxial apparatus for samples of synthetic GHBS under loading
66
and unloading at a given strain rate, and found that the samples first contracted and then dilated
67
owing to the existence of hydrate. Hyodo et al. (2013a, 2013b, 2014) conducted some triaxial
68
compression tests for the synthetic GHBS with Toyoura sands, and systematically investigated the
69
impacts of initial effective confining stress, pore pressure, hydrate saturation, temperature, and
70
density on the strength and deformation characteristics of GHBS. Kajiyama et al. (2017a, 2017b)
71
carried out several triaxial compression tests for the synthetic GHBS made with glass beads and
72
Toyoura sands, and studied the impact of grain properties on the strain and stress responses of
73
GHBS as well as the bonding and debonding mechanism of hydrate between particles. Hyodo et al. 3
74
(2017) conducted several triaxial compression tests for studying impacts of fines content as well as
75
density on the shear response of hydrate-free sediments and GHBS, and concluded that the increase
76
of fines content significantly enhanced the peak strength and dilative behavior. Choi et al. (2018)
77
experimentally investigated the mechanical response of synthetic GHBS through the triaxial tests of
78
single-stage and multi-stage, focusing on the strength, stiffness, and dilatancy of GHBS during
79
shearing.
80
In addition to the above experimental studies on synthetic GHBS, only a few experiments have
81
been conducted for samples of natural GHBS extracted from gas hydrate deposits due to the
82
difficulty and high cost of in-situ sampling. Winters et al. (2007) conducted several triaxial
83
compression tests for the GHBS extracted from the Mallik 2 L-38, Mackenzie Delta, Canada, and
84
indicated that the GHBS was higher than the natural hydrate-free sample in shear strength. Masui et
85
al. (2008) and Yoneda et al. (2015a, 2015b) investigated the mechanical behavior of GHBS
86
collected from the Eastern Nankai Trough, Japan. More recently, Yoneda et al. (2019) conducted
87
isotropic consolidation tests and triaxial compression tests on the samples of GHBS taken from the
88
Krishna-Godavari Basin, offshore India.
89
At the same time, the mechanical behaviors of GHBS have been investigated theoretically to
90
evaluate the wellbore stability or the submarine slope failure during hydrate dissociation (Sasaki,et
91
al., 2018; Yan et al., 2018; Malagar et al., 2019; Song et al., 2019; Yu et al., 2019). Up to now,
92
different types of constitutive models for GHBS have been presented, such as the nonlinear models
93
by Dai et al. (2012), Miyazaki et al. (2012), Pinkert and Grozic (2014), Pinkert et al. (2015), and
94
Yan et al. (2017), the elasto-plastic models by Klar et al. (2010), Sultan and Garziglia (2011),
95
Uchida et al. (2012, 2016), Sun et al. (2015), Lin et al. (2015, 2017), Shen et al. (2016), Gai and
96
Sánchez (2017), Sánchez et al. (2017), Yan and Wei (2017), and Zhou et al. (2018), the
97
elasto-visco-plastic model by Kimoto et al. (2010), and the hypoplastic model by Zhang (2017). For
98
example, Miyazaki et al. (2012) developed a nonlinear Duncan-Chang type of model for GHBS.
99
Although this model reasonably predicted the mechanical behavior of GHBS in the stiffness and
100
strength, it did not properly describe the volume variation of GHBS during loading and unloading.
101
Uchida et al. (2012, 2016) presented a Cam-clay type of model for GHBS, and verified its modeling
102
ability by comparing with the published test data. Sun et al. (2015) proposed a
103
thermodynamics-based model of GHBS with the concept of critical state, and confirmed that the
104
model well captured the stress softening and dilation behaviors of GHBS under different
105
morphologies and saturations of hydrate. Lin et al. (2015, 2017) proposed a subloading type of 4
106
model within the framework of spatial mobilized plane. The model capability was verified by
107
comparing with the triaxial compression test data of the synthetic and in-situ core specimens under
108
different hydrate saturations and initial effective confining stresses. Shen et al. (2016) developed a
109
state-dependent constitutive equation for GHBS. The simulation result showed that the model
110
reasonably predicted the constitutive response of GHBS at different hydrate saturations, initial
111
effective confining stresses, and void ratios. Gai and Sánchez (2017) and Sánchez et al. (2017)
112
presented the elasto-plastic model of GHBS in the hierarchical single-surface framework, and
113
concluded that the model had a good performance by comparing with the experimental data for the
114
synthetic and natural core samples with different hydrate morphologies, hydrate saturations, and
115
initial effective confining stresses. More recently, Sánchez et al. (2018) applied their constitutive
116
model and THMC to simulate laboratory experiments and field production tests.
117
The above works on GHBS are similar to the experimental, theoretical, and numerical
118
investigations on the fouling of nanoparticles (Nikkhah et al., 2015a, 2015b; Kamalgharibi et al.,
119
2016; Hemmat Esfe et al., 2018) and the CO2 sequestration (Hajirezaie et al., 2015; Rostamian et al.,
120
2016, 2019; Golkhou et al., 2019). They are closely related in thermodynamic properties. In
121
particular, when CO2 is injected into GHBS, the kinetics of CO2-CH4 hydrate in the reservoir
122
strongly depends on the physical mechanism of gas hydrate at the molecular level (Boswell et al.,
123
2017).
124
The existing constitutive models of GHBS are basically based on the coaxial assumption between
125
the directions of principal stress and plastic strain increment, and are mostly state-independent. To
126
overcome these shortcomings, this study presents a state-dependent multishear bounding surface
127
constitutive relation for GHBS. On the basis of plasticity slip theory, the formulation is obtained by
128
splitting a macro constitutive response into a macro volume response and a series of micro shear
129
responses in spatial distributions related to microshear structures. Each microshear structure
130
describes micro shear and dilatancy responses in three orthogonal orientations. A relation of the
131
micro stress–strain and a relation of the micro stress–dilatancy are established for each orientation
132
of the microshear structure. The hydrate-related state parameter and the evolution law of hydrate
133
bonding and debonding are introduced into the model. Moreover, in order to consider the impact of
134
hydrate dissociation on the deformation of GHBS, a simple collapse calculation method is proposed.
135
The validity of the constitutive relation is verified by comparing the calculated values with the test
136
results of GHBS under pore pressures of 8–15 MPa, temperatures of 1–10 °C, initial void ratios of
137
0.5–1.0, hydrate saturations of 0–0.8, and initial effective confining stresses of 1–5 MPa. 5
138
2. Model description
139
2.1. Formulation of multishear model
140
GHBS are assemblages of soil particles with large amounts of different sizes and shapes as well
141
as hydrates hosted in sediments. Soil particles interact with each other and endure external loads.
142
The role of hydrate differs from its accumulation habit. Fig. 1 gives the schematic diagram of the
143
multishear model through a unit sphere that represents the elementary volume of GHBS. There are a
144
number of randomly oriented microcontacts in space (Fig. 1(a)), which contribute to the macro
145
response of GHBS under loading. Following the concept of plasticity slip (Taylor, 1938), a macro
146
constitutive response of GHBS is split into a macro volume response and a series of micro shear
147
responses in spatial distributions related to microcontacts in different directions, in which the micro
148
shear contributes mainly to the macro deviatoric response (Fang et al., 2017; Fang et al., 2019). In
149
principle, this splitting approach somewhat resembles the volumetric-deviatoric decomposition used
150
in the microplane model (Caner and Bažant, 2013). Note that the micro shear here is not defined at
151
the real micro scale, but at the virtual equivalent scale smaller than the macro scale. As shown in
152
Fig. 1(b), for any microcontact colored in orange, a unit normal vector n on the sphere can be
153
found, which is consistent with the normal direction of the microcontact. For each microcontact, a
154
virtual microshear structure is assumed, which has the same direction as the microcontact (Fig. 1(c)).
155
The component in n is expressed with ni , in which the subscript i represents the global
156
coordinate xi . l and m refer to the orthogonal unit vectors, defining two tangential orientations
157
of the microshear structure. The components in l and m are denoted by li and mi , respectively.
158
To establish the relation of the micro shear and macro deviatoric responses, a macro strain is split
159
into the volumetric and deviatoric parts, i.e., ε ij = ε v δ ij / 3 + eij , where ε ij represents the macro
160
strain tensor, ε v = ε ijδ ij refers to the macro volumetric strain, eij represents the macro deviator
161
strain tensor, and δ ij denotes the delta of Kronecker. Generally, projecting a macro deviator strain
162
tensor onto the microcontact forms two tangential strain components corresponding to the l and
163
m directions and one deviatoric normal strain component corresponding to the n direction. As the
164
projected strain components of a macro deviator strain tensor, there is no difference among them
165
except for the direction. Hence, they can all be called micro shear strains, expressed as (see
166
Appendix A.1 for details or Fang et al. (2017))
167
γ ( k ) (n)=N ij( k ) (n )ε ij ( k =1, 2, 3) 6
(1)
168
where (k) denotes the kth component in the microshear structure; 1, 2, and 3 in k mean the
169
directions of l , m , and n in the microshear structure, respectively; γ (k ) (n ) refers to the kth
170
component of a micro shear strain in the microshear structure of the n direction; and N ij(k ) (n )
171
represents the projection tensor for the kth component of a micro shear strain in the microshear
172
structure of the n direction, which is given by equation (A4) in the appendix.
173
The kth component of a micro shear stress is expressed as τ (k ) (n ) in the microshear structure of
174
the n direction and its direction is assumed to be consistent with the corresponding micro shear
175
strain. In general, it is not possible for the micro strain and stress on the microcontact to be the
176
tensor projections for the macro strain and stress, respectively. Hence, a static equilibrium between
177
the macro and micro levels is approximately enforced in terms of the virtual work principle for all
178
strains and stresses. With the hypothesis that all microshear structures are mutually independent as
179
well as the micro strain and stress distribute uniformly in each microshear structure, the increment
180
of macro effective stress is derived as (see Appendix A.1 for details or Fang et al. (2017)) M
181
3
dσ ij = dpδij + ∑∑ 2w(m) Nij(mk )dτ (mk )
(2)
m=1 k =1
182
where w ( m ) denotes the weight coefficient in the mth microshear structure; σ ij represents the
183
macro effective stress tensor; p = σ ii /3 refers to the macro mean effective stress; m = 1, 2, ..., M
184
represent a collection of integral points, each of which corresponds to the microshear structure with
185
the direction of n ( m ) ; M denotes the total integral point number, which is set to 21 for an efficient
186
integration formula with an enough accuracy proposed by Bažant and Oh (1985); and (mk)
187
represents the kth component of a variable in the mth microshear structure with the direction of
188
n ( m ) , such as N ij( mk ) = N ij( k ) (n ( m ) ) and τ ( mk ) = τ ( k ) (n ( m ) ) . Equation (2) indicates that the increment
189
of macro stress consists of a direction-independent volumetric part and a deviatoric part of the
190
summation of micro shear stresses for all microshear structures.
191
2.2. Critical state and state parameter
192
The impact of hydrate on the sediment structure is dependent on the saturation and accumulation
193
habit of hydrate. There exist three main kinds of hydrate morphologies in sediments, i.e., the
194
pore-filling hydrate, the cementing hydrate, and the load-bearing hydrate (Masui et al., 2005), as
195
shown in Fig. 2. The pore-filling hydrate can significantly affect the permeability of sediments,
196
which also affects the mechanical behavior of sediments when the saturation of hydrate is high. 7
197
Particularly, the hydrate becomes the load-bearing hydrate if its saturation rises to the order of 0.3–
198
0.4. The cementing hydrate acts as the bonding material at inter-granular contacts and can
199
significantly increase the strength and stiffness of sediments by bonding near particles. The
200
load-bearing hydrate is a part of the sediment skeleton and can strengthen the skeleton structure of
201
sediments.
202
Hyodo et al. (2013a, 2013b) performed a set of triaxial compression tests for samples of synthetic
203
GHBS to study impacts of pore pressure, temperature, initial void ratio, hydrate saturation, and
204
initial effective confining stress on the deformation and strength behaviors of GHBS. The samples
205
with host sands were made through the partial water saturation approach, which belonged to the
206
cementing hydrate. Fig. 3 illustrates the relationship of the critical void ratio against mean effective
207
stress for synthetic GHBS, indicating that the location of critical state line is significantly affected
208
by the saturation of hydrate. Although the experimental data are somewhat scatted, the linear
209
relationship in the e − ln p plane can be approximately described for samples with any hydrate
210
saturation, where e refers to the void ratio. With the increase of the saturation of hydrate, the
211
critical state line moves upward. The relationship of the critical void ratio against mean effective
212
stress is described by
213
ec = eΓ − λc ln( p / pa )
214
where ec represents the void ratio at a critical state; λc refers to the material constant; pa
215
denotes the atmospheric pressure; and eΓ represents the critical void ratio parameter which
216
depends on the saturation of hydrate, defined as (Shen et al., 2016)
(3)
217
eΓ = eΓ 0 + ac Shnc
218
where eΓ 0 refers to the material parameter for host materials; ac and nc denote material
219
parameters; and S h represents the hydrate saturation. Fig. 4 illustrates the variation of eΓ with the
220
hydrate saturation obtained from Fig. 3, in which its fitting curve is also given.
(4)
221
Fig. 5 shows the relation of the critical stress ratio versus hydrate saturation for synthetic GHBS,
222
where pp and T represent the pore pressure and temperature, respectively. It can be seen from
223
this figure that although the test data are relatively scatted, the trend that the critical stress ratio is
224
directly proportional to the hydrate saturation is very obvious. Therefore, for simplicity, it is
225
assumed that the critical stress ratio increases linearly with the hydrate saturation, given by
226
M c = M c0 (1 + aM Sh )
8
(5)
227
where M c0 and M c represent the critical stress ratios of hydrate-free sediments and GHBS,
228
respectively; and aM denotes the material parameter, reflecting the influence of hydrate on the
229
critical stress ratio.
230
Some studies indicated that the stress and density affect the mechanical behavior of granular soils
231
(Been and Jefferies, 1985; Sun et al., 2008; Zhao and Guo, 2013). Similarly, to consider their
232
impacts in the model, the state parameter ψ for sands (Been and Jefferies, 1985) is applied. This
233
state parameter represents the interval of the critical state and current void ratios at the same mean
234
effective stress, expressed as
ψ = e − ec
235 236
(6)
2.3. Evolution law of hydrate bonding and debonding
237
Similar to cemented soils (Nguyen et al., 2014), the inter-particle bonding of hydrate can be
238
defined by a tensile strength of GHBS. This strength is called the bonding strength, which can be
239
calculated in terms of the experimental data. Fig. 6 illustrates the impact of saturation in hydrate on
240
the peak strength line of GHBS in the p − q plane, where q and qf are the deviator stress and
241
peak strength, respectively. From this figure, it can be observed that the mean effective stress at the
242
deviator stress of zero represents the initial value of the bonding strength, which is dependent on the
243
saturation of hydrate. Therefore, the relation of the initial value of the bonding strength and the
244
saturation of hydrate can be obtained from the test results as plotted in Fig. 6, which is
245
approximately expressed as
pt0 = at Sh
246
(7)
247
where pt0 represents the initial bonding strength; and a t denotes the material parameter. Fig. 7
248
shows the fitting result in the relation of the initial bonding strength versus hydrate saturation.
249
To describe the cementation and deterioration of hydrate during deformation, a simple evolution
250
law for the bonding strength is used, expressed as (Uchida et al., 2012, 2016; Sun et al., 2015; Lin
251
et al., 2015, 2017; Shen et al., 2016; Gai and Sánchez, 2017; Sánchez et al., 2017; Yan and Wei,
252
2017)
dpt = −kt pt dε qp
253
(8)
254
where p t refers to the current bonding strength; k t represents the model parameter reflecting the
255
p debonding rate; and ε q denotes the equivalent plastic shear strain. The integral form of equation
256
(8) is written as 9
pt = pt0 exp(−ktε qp )
257 258
(9)
2.4. Characteristic surfaces and lines
259
To consider the effect of bonding strength in the constitutive relation of GHBS, the macro
260
modified effective stress tensor and mean stress are expressed as σˆ ij = σ ij + ptδ ij and
261
pˆ = p + pt , respectively. Following the bounding surface theory, two types of characteristic
262
surfaces are described, one is conical surface and the other is flat cap. The conical surface depends
263
on the elasto-plastic loading owing to a change in the modified stress ratio defined as
264
rˆij = (σˆ ij − pˆ δij ) / pˆ , and the cap depends on the elasto-plastic loading owing to a change in the
265
modified mean effective stress at a unchanged modified stress ratio. Fig. 8 illustrates four
266
cone-shaped characteristic surfaces whose apexes are located at the origin. The bounding surface
267
represents the condition where the modified stress ratio rises to its peak. The surface of critical state
268
denotes the condition where any change in stress and void does not occur while the plastic shear
269
deformation unceasingly increases. The dilatancy surface refers to the condition where the
270
transformation from contraction to dilation takes place. The surface of maximum stress ratio
271
represents the condition where its maximum value occurs in loading histories. The shapes of four
272
conical surfaces are similar, given by the invariant of the modified stress ratio, the yielding shape
273
function, and their stress ratios at triaxial compression as (Fang et al. 2017)
Fi =Rˆ − Mi g (θˆ) = 0 (i=m, b, c, d)
274
(10)
275
where Fm , Fb , Fc , and Fd refer to the macro maximum stress ratio, bounding, critical state,
276
dilatancy surfaces, respectively; Rˆ represents the macro modified stress ratio, given by
277
Rˆ =
278
and dilatancy stress ratios under the condition of triaxial compression, respectively; and g (θˆ )
279
denotes the yielding shape function; θˆ refers to the angle of Lode. According to equation (10),
280
M m is equal to the maximum value of Rˆ / g (θˆ ) during loading. M d and M b are expressed as
281
(Li, 2002)
282 283
3 / 2rˆij rˆij ; M m , M b , M c , and M d mean the macro maximum, bounding, critical state,
M d = M c exp(ndψ ) , M b = M c exp(−nbψ ) where nd and nb are model constants.
10
(11)
284
To involve the influence of the intermediate principal stress in the model, the yielding criterion of
285
Matsuoka and Nakai (1974) or the yielding criterion of Lade and Duncan (1975) is applied. For the
286
yielding criterion of Matsuoka and Nakai, g (θˆ ) is written as (Matsuoka et al., 1999)
(
1 g (θˆ) = Rˆ 3 6
287
( Iˆ Iˆ − Iˆ ) ( Iˆ Iˆ − 9Iˆ ) − 1) 1 2
3
1 2
(12)
3
288
3 / 6 refer to where Iˆ1 = σˆ ij , Iˆ2 = (σˆ ii2 − σˆ rsσˆ sr ) / 2 , and Iˆ3 = σˆ ijσˆ jkσˆ ki / 3 − σˆ rsσˆ srσˆ mm / 2 + σˆ nn
289
invariants of the modified stress tensor. For the yielding criterion of Lade and Duncan, g (θˆ ) is
290
written as (Yao and Sun, 2000)
291 292
Rˆ 1 ˆ 3 g (θˆ) = I 3 /pˆ 1 − 3 pˆ 2
1 −1 ˆ 3 cos 3 cos ( − I 3 /pˆ )
−1
−1
(13)
The bounding surface of flat cap is expressed as (Fang et al. 2017)
Fp =pˆ − pˆ m = 0
293
(14)
294
where Fp refers to the bounding surface of flat cap; and pˆ m denotes the maximum modified
295
mean effective stress in loading histories.
296 297
Similarly, four characteristic lines for the kth component in the mth microshear structure, as illustrated in Fig. 9, are expressed as (Fang et al. 2017)
Fi ( mk ) =τ ( mk ) − ri( mk ) pˆ = 0 (i=m, b, c, d)
298
(15)
299
where Fm(mk ) , Fb(mk ) , Fc(mk ) , and Fd(mk ) denote the micro maximum stress ratio, bounding, critical
300
state, and dilatancy lines for the kth component in the mth microshear structure, respectively;
301
r (mk ) = τ (mk ) /pˆ represents the kth component of the micro stress ratio in the mth microshear
302
structure; and rm( mk ) , rb( mk ) , rc( mk ) , and rd( mk ) mean the micro maximum, bounding, critical state,
303
and dilatancy stress ratios for the kth component in the mth microshear structure, respectively.
304
These micro stress ratios correspond to different micro stress states. Their physical meanings are
305
similar to those corresponding to the macro modified stress ratios. According to equation (15),
306
rm( mk ) is set to the maximum value of r (mk ) in loading histories. ri ( mk ) (i=b, c, d) is given by
307
equation (A28) in the appendix.
11
308
2.5. Macro volumetric deformation
309
The macro volumetric strain without that due to dilatancy usually consists of elastic and plastic
310
parts, which depends on the modified mean effective stress. The incremental relationship of them is
311
given by (see Appendix A.2 for details or Fang et al. (2017))
dε v − dε vd =
312
dpˆ 1 1 dpˆ + h( pˆ − pˆ m ) dpˆ Ke Kp dpˆ
(16)
313
where K e denotes the elastic bulk modulus; K p represents the plastic bulk modulus; ε vd refers
314
to the macro volumetric strain due to dilatancy; h represents the Heaviside step function; and
315
refers to the Macaulay brackets.
316
The elastic bulk modulus is given by Ke =
317
1+ e
κ
ˆ a ) 0.5 (1 + aG S h )( pp
(17)
318
where aG refers to the material parameter; the term aG S h represents the influence of hydrate on
319
the modulus of GHBS; and κ denotes the material parameter related to the unloading in isotropic
320
compression, expressed as κ = κ 0 [(1 + e) / (2.97 − e)]2 where κ 0 represents the model constant
321
for host materials.
322 323
The plastic bulk modulus is defined as (Wang et al., 1990)
Kp =
Mc g(θˆ) 1+ e ˆ a )0.5 (1+ aGSh )( pp λ −κ Mc g(θˆ) − Rˆ
(18)
324
where λ represents the material parameter related to the isotropic compression, expressed as
325
λ = λ0 [(1 + e ) / (2.97 − e )]2 where λ0 denotes the model constant for host materials.
326 327
The macro volumetric strain due to dilatancy is the sum of the micro volumetric strains due to dilatancy for all microshear structures, written as (Fang et al., 2017) M
328
3
( mk ) dε vd = ∑ ∑ 2 w( m ) dε vd
(19)
m =1 k =1
329
where ε vd( mk ) denotes the kth component of the micro volumetric strain due to dilatancy in the mth
330
microshear structure. Following the cyclic torsional shear test result of sands, the micro stress–
331
dilatancy relation is given by (Fang et al., 2017)
332
dε vd( mk ) = d1 (±rd( mk ) − r ( mk ) )dγ p(mk )
12
(20)
333
where d1 represents the micro dilatancy parameter, which is assumed to be the same for all
334
microshear structures; ± is set as positive for dγ p(mk ) > 0 and negative for d γ p(mk ) < 0 ; and
335
γ p(mk ) denotes the kth component of the micro plastic shear strain in the mth microshear structure,
336
given by equation (A15) in the appendix.
337
2.6. Micro shear deformation
338
The micro shear strain usually consists of elastic and plastic parts. The former depends on the
339
micro shear stress, and the latter depends on the micro shear stress ratio and modified mean
340
effective stress. According to the plasticity theory of bounding surface (Wang et al., 1990), their
341
incremental correlation is expressed as (see Appendix A.3 for details or Fang et al. (2017))
342
dγ ( mk ) =
1
dτ ( mk ) + ( mk )
Ge
1 Gp( mk )
pˆ dr ( mk ) +
1
h( pˆ − pˆ m ) ( mk )
Hp
dpˆ ( mk ) r dpˆ dpˆ
(21)
343
where γ ( mk ) denotes the kth component of the micro shear strain in the mth microshear structure;
344
Gp(mk ) and H p(mk ) refer to the micro plastic shear moduli related to d r ( mk ) and dpˆ for the kth
345
component in the mth microshear structure, respectively; and Ge(mk ) represents the micro elastic
346
modulus for the kth component in the mth microshear structure, which depends on the macro elastic
347
shear modulus and is given by equation (A29) in the appendix. The macro elastic shear modulus
348
Ge is expressed as (Richart et al., 1970)
349
Ge = G0 (1 + aG S h )
(2.97 − e) 2 ˆ a ) 0.5 ( pp 1+ e
(22)
350
where G0 refers to the elastic shear modulus parameter for host materials; and G0 (1 + aG S h ) = Gh0
351
denotes the elastic shear modulus parameter for GHBS. Fig. 10 shows the Gh0 − S h relationship
352
from the experimental data (Hyodo et al., 2013a, 2013b).
353 354
The micro plastic shear modulus Gp(mk ) related to d r ( mk ) is defined as (Fang et al., 2017)
r ( mk ) ρ ( mk ) Gp( mk ) = h1Ge(mk ) b( mk ) 1( mk ) −1 rm ρ1
(23)
355
where h1 denotes the model parameter; ρ1( mk ) = r (mk ) − rr(mk ) represents the interval between the
356
current and back values of the kth micro stress ratio in the mth microshear structure; rr(mk ) denotes
357
the back value of the kth micro back stress ratio in the mth microshear structure, defined as the
358
origin during a virgin loading or the latest inflection location during a reverse loading; and
13
359
ρ1( mk ) = ± rm(mk ) − rr(mk ) represents the interval between the maximum and back values of the kth
360
micro stress ratio in the mth microshear structure, where ± is set as positive for d r ( mk ) > 0 and
361
negative for d r ( mk ) < 0 .
362
The micro plastic shear modulus H p(mk ) related to dpˆ is defined as (Fang et al., 2017)
H
363
( mk ) p
( mk ) ( mk ) c e ( mk )
= h2G
r r
ρ2 ρ2
(24)
364
where h2 denotes the model parameter; ρ 2 = pˆ − pˆ r refers to the interval between the current
365
and back values of the modified mean effective stress; pˆ r represents the back value of the
366
modified mean effective stress, defined as the origin during a virgin loading or the latest inflection
367
location during a reverse loading; and ρ 2 = pˆ m − pˆ r denotes the interval between the maximum
368
and back values of the modified mean effective stress where ρ 2 = pˆ r for d pˆ < 0 .
369
2.7. Collapse deformation due to hydrate dissociation
370
Klar et al. (2010) and Sánchez et al. (2017) proposed different methods to predict the collapse
371
deformation caused by hydrate dissociation. The Klar et al.’s method follows the concept of stress
372
relaxation. Due to the dissociation of hydrate, the effective stress in the hydrate skeleton gradually
373
releases, resulting in the decrease of the total effective stress in the skeleton of GHBS under
374
constant strain. In order to restore the released stress, an additional deformation is required. Because
375
the initial elastic modulus is applied in the Klar et al.’s method, the collapse deformation may be
376
underestimated. To overcome this drawback, an improved method is proposed. During hydrate
377
dissociation without external load, it can be assumed that the sediment is self-balanced through its
378
skeleton structure, and the effective stress remains unchanged from the time of t to the time of t+dt,
379
i.e.,
380
σ ijt = ( K eqε vδ ij + 2Geq eij )t = σ ijt +dt = ( K eqε vδ ij + 2Geq eij )t +dt
381
where the superscripts t and t+dt of a variable represent its values at the time of t and the time of
382
t+dt, respectively; and K eq and Geq denote the equivalent bulk and shear moduli, respectively.
383
Since σ ijt , ε vt , and eijt are known at the time of t, K eqt and Geqt can be approximately obtained
384
by using the linear elastic relationship between the stress and strain. By assuming that the
385
increments of the equivalent bulk and shear moduli are proportional to the variation of hydrate
386
saturation during the time increment of dt, K eqt + dt and Geqt + dt are approximately given by
14
(25)
387
K eqt +dt = K eqt +dK eq = K eqt (1+dSh )
(26)
388
Geqt +dt = Geqt +dGeq = Geqt (1+dSh )
(27)
389 390
From equations (25)–(27), the increment of the collapse strain produced by hydrate dissociation is given by
dε ijdis = ε ijt +dt − ε ijt = −[ pt / (3K eqt )δ ij + (σ ijt − p tδ ij ) / (2Geqt )]dSh / (1 + dSh )
391 392
where ε ijdis refers to the collapse strain caused by hydrate dissociation.
393
2.8. Constitutive formulation
394 395
(28)
The incremental relation of the stress and strain, including the effect of hydrate dissociation, is described as
dσ ij = Dijrs (dε rs − dε rsdis ) − dptδij
396
(29)
397
where Dijrs denotes the macro elasto-plastic stiffness tensor related to macro and micro quantities
398
(see Appendix A.4 for details), which enables the model to involve the influence of micro quantities
399
as well as their evolutions at different directions of microshear structures on the constitutive
400
response of GHBS.
401
3. Parameter calibration
402
The model has 17 parameters, 11 of which are associated with the mechanical behavior of host
403
materials, and the remaining 6 are associated with hydrate. The 11 model parameters for host
404
materials include two elastic modulus parameters ( G0 and κ 0 ), three critical state parameters
405
( M c0 , eΓ 0 , and λ c ), four plastic modulus parameters ( h1 , h2 , nb , and λ0 ), and two dilatancy
406
parameters ( d1 and nd ). These parameters for host materials are determined through the following
407
approaches: (1) G 0 is obtained by the relation between the deviator stress and axial strain for
408
triaxial tests. λ0 and κ 0 are determined by the e − p relations for consolidation tests. κ 0 is
409
also set by κ 0 = 3(1 − 2ν ) / [2(1 + ν ) G 0 ] , where ν refers to the ratio of Poisson; (2) M c0 , eΓ 0 ,
410
and
411
nd = ln( M d / M c ) / ψ d , where ψ d and M d take the values of ψ and stress ratio when the phase
412
transformation takes place, respectively. nb is determined by nb = ln( M c / M b ) / ψ b , where ψ b
413
and M b take the values of ψ and stress ratio when the peak of stress ratio occurs, respectively;
λc
are
directly
determined
by
triaxial
15
test
results;
(3)
nd
is
obtained
by
414
(4) d1 is obtained according to the volumetric and deviatoric strain data of triaxial tests; (5) h1
415
and h2 are set by the trial-and-error procedure to match optimally the model predictions with the
416
test results.
417
The additional 6 parameters for GHBS consist of one elastic modulus parameter ( aG ), three
418
critical state parameters ( ac , nc , and aM ), one bonding strength parameter ( a t ), and one strength
419
deterioration parameter ( k t ). aG is determined in terms of a variation in elastic shear modulus
420
with the saturation of hydrate, as shown in Fig. 10. ac and nc are obtained by the regression of
421
the critical void ratio with a variation in hydrate saturation, as illustrated in Fig. 4. aM is set by
422
defining that the critical stress ratio is dependent on the saturation of hydrate, as shown in Fig. 5.
423
a t for the initial bonding strength is given through the experimental peak strength data for GHBS,
424
as illustrated in Figs. 6 and 7. k t is calibrated by the trial-and-error procedure to simulate the
425
mechanical response of GHBS in tests.
426
4. Model performance
427
The model performance is evaluated by simulating the available experimental results published in
428
literatures. Five sets of triaxial compression tests for GHBS are used. Two sets are the tests for the
429
samples of synthetic GHBS (Masui et al., 2005; Hyodo et al., 2013a), and the other three sets are
430
the tests for the samples of natural GHBS taken from the Eastern Nankai Trough, Japan (Masui et
431
al., 2008; Yoneda et al., 2015b) and the samples of natural GHBS extracted from the Krishna–
432
Godavari Basin, offshore India (Yoneda et al., 2019). Moreover, the hypothetical torsional shear
433
tests during principal stress direction rotations are simulated to investigate the model ability for
434
predicting the non-coaxiality of GHBS. The simulation results are described as follows.
435
4.1. Triaxial compression tests for synthetic GHBS by Masui et al. (2005)
436
Masui et al. (2005) performed drained triaxial compression tests for synthetic cementing and
437
pore-filling types of GHBS, in which the pore-filling type of hydrate formed through the ice-seed
438
approach and the cementing type of hydrate formed through the partial water saturation approach.
439
The host materials were made of Toyoura sands. Tests were carried out under the pore pressure of 8
440
MPa, temperature of 5 °C, initial effective confining stress of 1 MPa. For the pore-filling type of
441
GHBS, the initial void ratio of the samples was 0.736 and the saturation of hydrate was in the extent
442
of 0.264–0.678. For the cementing type of GHBS, the initial void ratio of the samples was 0.605
16
443
and the saturation of hydrate was in the extent of 0.0–0.551. In the tests, the deviator stress response
444
data of all samples were successfully measured, but only the volume change data of the sample with
445
a hydrate saturation of 0.409 for the pore-filling hydrate and the samples with hydrate saturations of
446
0.0 and 0.407 for the cementing hydrate were obtained. The model parameters for host sands and
447
GHBS used in simulations are given in Table 1. Note that ac and nc are set in two regions
448
according to the properties of pore-filling hydrate.
449
Fig. 11 compares the predicted values with the experimental results for the pore-filling type of
450
GHBS, where e0 and p c0 refer to the initial values of the void ratio and effective confining stress,
451
respectively. As can be observed from this figure, most of the predicted deviator stresses and
452
volumetric strains agree well with those of the test data. Only the predicted deviator stress responses
453
for the hydrate saturation of 0.678 are larger than the experimental values when the axial strain
454
exceeds 8%. The experimental deviator stress for the hydrate saturation of 0.678 is quite different
455
from that in other cases. After the axial strain of 8%, its value decreases significantly. The reason
456
for this phenomenon was not explained in the paper of Masui et al. (2005), which may be that the
457
hydrate partially decomposed during the experiment of this case. Fig. 12 compares the calculated
458
values with the test data for the cementing type of GHBS. In the relation of the deviator stress and
459
axial strain, the calculated values are consistent with the test data before the axial strain of 9%, but
460
the calculated values deviate from the test data after the axial strain exceeds 9–12%. In the relation
461
of the volumetric and axial strains, there is a certain difference between the predicted and
462
experimental values, especially for the case with a hydrate saturation of 0.407. Moreover, it can be
463
observed from Figs. 11 and 12 that the peak strength, initial stiffness, and volumetric dilation for
464
the cementing samples are larger than those for the pore-filling samples although their saturations of
465
hydrate are almost the same. The calculated results show that the constitutive relation captures the
466
different behaviors of GHBS for two accumulated types of hydrates, particularly in the deviator
467
stress response, as observed in the tests.
468
4.2. Triaxial compression tests for synthetic GHBS by Hyodo et al. (2013a)
469
Hyodo et al. (2013a) conducted drained triaxial compression tests for synthetic GHBS. Toyoura
470
sands were employed as host materials. The initial void ratio of the samples was in the range of
471
0.629–0.842 and the saturation of hydrate was in the extent of 0.0–0.695. Tests were performed
472
under different pore pressures (5, 10 and 15 MPa), temperatures (1, 5, and 10 °C), and initial
473
effective confining stresses (1, 3, and 5 MPa). The model parameters for host sands and GHBS used
17
474
in simulations are given in Table 1. Fig. 13 compares the predicted values with the test results of
475
GHBS containing different hydrate saturations under the pore pressure of 10 MPa, temperature of
476
5 °C, and initial effective confining stress of 5 MPa. From this figure, it can be observed that the
477
model can generally well predict the variations of the deviator stress and volumetric strain with the
478
axial strain. The simulated and experimental results show that, in addition to the great dilative
479
behavior of the sample containing a high hydrate saturation of 0.531, the samples exhibit the
480
dominant compression in volume and the behavior of strain hardening. Furthermore, the strength
481
and initial stiffness are proportional to the saturation of hydrate.
482
Fig. 14 shows the comparison between the simulated results and test data for GHBS with almost
483
the same hydrate saturation and initial void ratio under different pore pressures (10 and 15 MPa),
484
temperatures (1, 5, and 10 °C), and initial effective confining stresses (1, 3, and 5 MPa). In the
485
relation of the deviator stress and axial strain, the calculated values well consist with the test data
486
for two cases under initial effective confining stresses of 1 and 5 MPa, while the calculated values
487
for the initial effective confining stress of 3 MPa are agreement with the test data under the
488
temperature of 1 °C. In the relation of the volumetric and axial strains, the predicted values are very
489
consistent with the test data under the initial effective confining stress of 1 MPa, while the
490
calculated values are quite different from the test data under initial effective confining stresses of 3
491
and 5 MPa. In addition, when the initial effective confining stress is 3 MPa, the experimental
492
deviator stress responses at 5 and 10 °C are very close, while the experimental deviator stress
493
response at 1 °C is significantly greater than that at other two temperatures; the experimental
494
volume responses at these three temperatures do not change much with temperature. These test
495
results show that the decrease of temperature leads to the nonlinear increase of strength, while the
496
volume deformation is less affected by temperature. In contrast, the simulation results with different
497
temperatures are very close under the same initial effective confining stress because the impact of
498
temperature is not directly considered in this model.
499
4.3. Triaxial compression tests for natural GHBS by Masui et al. (2008)
500
Masui et al. (2008) carried out drained triaxial compression tests for natural GHBS taken from
501
the Eastern Nankai Trough, Japan. The natural GHBS were mainly composed of silty sands. Four
502
samples of the natural GHBS were sheared using the triaxial compression apparatus. The initial
503
void ratios of four samples were 0.792, 0.961, 0.972, and 1.033, respectively. The corresponding
504
saturations of hydrate were 0.376, 0.225, 0.094, and 0.077, respectively. Tests were performed
18
505
under the pore pressure of 9 MPa, temperature of 5 °C, and initial effective confining stress of 1
506
MPa. The model parameters for host sands and GHBS used in simulations are given in Table 1. Fig.
507
15 compares the calculated values with the experimental data at different saturations of hydrate. As
508
can be observed from this figure, the simulated values generally consist with the experimental
509
results in the relation of the deviator stress and axial strain and in the relation of the volumetric and
510
axial strains except for the volumetric strain response of the sample containing a hydrate saturation
511
of 0.094. The experimental and predicted results indicate that the peak value in strength and the
512
initial value in stiffness for natural GHBS are directly proportional to the saturation of hydrate. The
513
higher the saturation of hydrate, the larger the volumetric change of the sample. The sample shows
514
a significant strength reduction and volume change when the saturation of hydrate is equal to 0.376.
515
In contrast, the sample shows a low level of softening and dilative response when the saturation of
516
hydrate is equal to 0.077. For the sample with a constant saturation of hydrate, softening takes place
517
after the peak in strength and its reduction is proportional to the saturation of hydrate.
518
4.4. Triaxial compression tests for natural GHBS by Yoneda et al. (2015b)
519
Yoneda et al. (2015b) performed drained triaxial compression tests for natural GHBS taken from
520
the Eastern Nankai Trough, Japan. The natural GHBS were mainly composed of silty sands. Three
521
samples of the natural GHBS and their reconstituted samples were sheared by using the triaxial
522
compression apparatus. The saturations of hydrate for three samples were 0.38, 0.7, and 0.79,
523
respectively. Two samples with hydrate saturations of 0.38 and 0.79 were successfully tested, but
524
the volume change data for the sample with a hydrate saturation of 0.7 were not obtained due to the
525
liquid leakage in the cell. Therefore, the test for the sample with a hydrate saturation of 0.7 is not
526
simulated because its deviator stress response was also influenced by the liquid leakage. The initial
527
void ratios of all the samples were very similar in the range of 0.645–0.669. Tests were conducted
528
under the pore pressure of about 13 MPa, temperature of 10 °C, and initial effective confining stress
529
of 1.5 or 1.6 MPa. The model parameters for natural GHBS and their reconstituted samples used in
530
simulations are given in Table 1. Fig. 16 compares the calculated values with the experimental data
531
on the responses in deviator stress and volumetric strain for natural GHBS and their reconstituted
532
samples. From this figure, it can be seen that the model performance is generally satisfactory with
533
good predictions in terms of the stiffness and strength behavior as well as in terms of the volume
534
response except for the volumetric strain of the sample containing a hydrate saturation of 0.79. The
535
initial stiffness, peak strength, and residual strength significantly depend on the saturation of
19
536
hydrate. The reconstituted sediments retain the same strength as host silty sands, which is basically
537
consistent with the in-situ strength of hydrate-free sediments. The volumetric strain responses also
538
strongly depend on the saturation of hydrate. The sample with a hydrate saturation of 0.79 exhibits a
539
significant dilative behavior.
540
4.5. Isotropic consolidation tests and triaxial compression tests for natural GHBS by Yoneda et al.
541
(2019a, 2019b)
542
Yoneda et al. (2019a, 2019b) conducted isotropic consolidation tests and drained triaxial
543
compression tests on the natural GHBS extracted from the Krishna–Godavari Basin, offshore India.
544
The natural GHBS were sediments rich in sands. Two types of isotropic consolidation tests for six
545
samples of the natural GHBS with hydrate saturations of 0.321–0.716 were performed by use of a
546
triaxial apparatus. The initial effective confining stress in situ was about 1 MPa. The results of the
547
first type of test are shown in Fig. 17(a), where samples were first loaded from 0.2 to 1 MPa, then
548
hydrates were dissociated, and finally samples were loaded to 12 MPa after hydrate dissociation.
549
The results of the second type of test are shown in Fig. 17(b), where samples were first loaded from
550
0.2 to 5 MPa, then hydrates were dissociated, and finally samples were loaded to 12 MPa after
551
hydrate dissociation. During the tests, unloading and reloading were carried out at different
552
effective confining stresses.
553
The drained triaxial compression tests for five samples of the natural GHBS were performed. The
554
hydrate saturations of five samples were 0.487, 0.552, 0.634, 0.734, and 0.822, respectively. The
555
corresponding initial void ratios were 0.712, 0.634, 0.831, 0.658, and 0.656, respectively. Tests
556
were carried out under the pore pressure of about 10 MPa, temperature of 4 °C, and initial effective
557
confining stress of 1 or 1.1 MPa. The experimental data of the deviator stress response and volume
558
change were successfully measured except for the volume change data of the sample with a hydrate
559
saturation of 0.634. Fig. 18 shows the experimental results of the deviator stress response and
560
volume variation with the axial strain. From this figure, it can be observed that the experimental
561
shear strength and dilation of GHBS were not always proportional to the saturation of hydrate when
562
the saturation of hydrate was relatively high. However, the previous experimental results (Miyazaki
563
et al., 2011; Hyodo et al., 2013a, 2013b) indicated that the strength and dilation of GHBS increased
564
with the enhancement of hydrate saturation. A possible reason for this phenomenon may be the
565
hydrate amount in the samples. The saturations of hydrate in these test samples are relatively high,
566
while the saturations of hydrate in the previous tests samples are usually low. The hydrate particles
20
567
enter into the voids at low saturation of hydrate, while the hydrate particles partially occupy the
568
voids and partially block the contacts of sediment particles at high saturation of hydrate. Therefore,
569
the distribution of hydrate particles in GHBS is related to the amount of hydrate, thereby changing
570
the equivalent skeleton void ratio, which in turn affect the shear strength and dilation of GHBS.
571
Although hydrates and fines are different kinds of materials, this characteristic may be somewhat
572
similar to that of the sand-fines mixtures (Thevanayagam et al., 2002).
573
To consider the above feature in simulations, the hydrate saturation is divided into two intervals,
574
and the equivalent skeleton void ratio parameters are set respectively. The model parameters of the
575
natural GHBS applied in simulations are illustrated in Table 1. Because the sediments of the sample
576
with a hydrate saturation of 0.321 is classified as silt (Yoneda et al., 2019a), κ 0 and λ 0 for this
577
sample are taken to be 0.01 and 0.03, respectively. Fig. 17 compares the model predictions with the
578
experimental data for isotropic consolidation tests. Note that the experimental void ratio at the mean
579
effective stress of 0.2 MPa for the two cases with hydrate saturations of 0.321 and 0.639 is quite
580
different from that in other cases. Therefore, for the simulations of these two cases, the void ratio
581
corresponding to the mean effective stress of 0.2 MPa is inversely determined according to the void
582
ratio at the mean effective stress of 1.0 or 1.1 MPa. It can be seen from Fig. 17 that the volume
583
response predicted by the model basically agree with that observed in the tests. Particularly, the
584
model can reasonably predict the collapse deformation caused by hydrate dissociation, indicating
585
the effectiveness of the proposed collapse calculation method. Fig. 18 shows the predicted and
586
experimental results for triaxial compression tests. From this figure, it can be seen that the
587
calculated values generally consist with the experimental data in the relation between the deviator
588
stress and axial strain and in the relation between the volumetric and axial strains.
589
4.6. Drained torsional tests during principal stress rotations
590
To investigate the ability of the model for predicting the non-coaxiality of GHBS, the drained
591
torsional shear tests under the pure rotation of principal stress directions are simulated. Since this
592
type of test for GHBS has not been reported, the test conditions for simulations are assumed
593
according to the existing triaxial compression tests of GHBS and torsional shear tests of sands. The
594
schematic diagram of the hollow cylindrical sample is showed in Fig. 19. The hydrate saturations of
595
three samples used in simulations are 0.0, 0.3, and 0.6, respectively. The initial void ratio and
596
effective confining stress are equal to 0.65 and 1 MPa, respectively. The pure stress rotation is
597
achieved by varying only the direction of the principal stress while maintaining the other conditions
21
598
unchanged. The mean effective stress, deviator stress, and intermediate principal stress ratio are
599
taken as 1 MPa, 1 MPa, and 0.5, respectively. The model parameters are set to be the same as those
600
applied in the simulation of the Yoneda et al.’s tests (Yoneda et al., 2015b), which are given in
601
Table 1. Fig. 19 shows the predicted variations of all strain components with the principal stress
602
direction during a cycle of continuous principal stress rotations. From this figure, it can be observed
603
that all strain components decrease with the increase of saturation in hydrate, indicating that the
604
strength and stiffness of sediments increase due to hydrate. Figs. 20 and 21 illustrate the predicated
605
stress paths and vectors of the total strain increment, as well as the predicated stress paths and
606
vectors of the plastic strain increment, respectively. The magnitudes of the increment vectors for
607
total and plastic strains in these two figures are define as
(d(ε a − ε θ )) 2 +(2dε aθ ) 2 / ds and
608
(d(ε ap − ε θp )) 2 +(2dε aθp ) 2 / ds , respectively, in which the stress increment ds is expressed as
609
ds = (d((σ a − σ θ ) / (σ a + σ θ ))) 2 +(d(2σ aθ / (σ a + σ θ ))) 2 . From these two figures, it can be observed
610
that with the enhancement of hydrate saturation, the increment in the plastic strain decreases and
611
inclines to the tangential direction along the circle, indicating that the non-coaxiality of GHBS
612
increases. Therefore, it can be said that the present constitutive relation is able to predict the
613
non-coaxial response of GHBS without new parameters.
614
5. Conclusions
615
This paper presented a state-dependent multishear constitutive formulation for predicting the
616
mechanical response of GHBS. The model follows the plasticity theory of bounding surface and
617
accords with the theory of critical state for granular soils. The effectiveness of the model was
618
confirmed by simulating different laboratory tests for synthetic and natural GHBS published in
619
literatures. The following conclusions are obtained:
620
(1) The model formulation was derived by splitting a macro constitutive response of GHBS into a
621
macro volume response and a series of micro shear responses in spatial distributions related to the
622
microshear structures. This splitting approach can provide a simple and effective means to predict
623
the complex three-dimensional constitutive response of GHBS by superposition of a macro
624
response in volume and individual responses of all microshear structures.
625
(2) A state parameter related to the hydrate saturation was introduced into the model, including
626
the joint impact of the effective confining stress and density on the constitutive responses of GHBS.
627
With this state parameter, the model not only accords with the theory of critical state for granular
22
628
soils, but also can capture the contractive and dilative behaviors for loose and dense GHBS in terms
629
of a single set of model parameters.
630
(3) The evolution law of hydrate bonding and debonding related to the plastic shear strain was
631
used to consider the effect of hydrate on the constitutive response of GHBS. This evolution law
632
characterizes the cementing mechanism of hydrate and the degradation of bonding during plastic
633
deformation, respectively.
634
(4) The simple method to calculate the collapse deformation produced by the dissociation of
635
hydrate was proposed. The validity of this method was verified by predicting the volume
636
compression of GHBS observed in the hydrate dissociation tests under constant stress.
637
(5) The comparisons of the experimental and simulated results for different drained laboratory
638
tests of synthetic and natural GHBS were performed. The result confirmed that the developed
639
model can not only predict the constitutive response of GHBS under isotropic consolidation and
640
triaxial compression, but also predict the non-coaxiality of GHBS during the pure rotation of
641
principal stress directions without new parameters. The model is suitable for the case of pore
642
pressures of 8–15 MPa, temperatures of 1–10 °C, initial void ratios of 0.5–1.0, hydrate saturations
643
of 0–0.8, and initial effective confining stresses of 1–5 MPa. However, since temperature is not
644
taken as an independent variable in the present model, and the influence of temperature on the
645
strength of GHBS cannot be ignored, temperature should be directly introduced into the constitutive
646
model in the future.
647
Acknowledgements
648
The authors would like to thank the reviewers for their creative suggestions and professional
649
comments. This work was financially supported by the National Natural Science Foundation of
650
China (No. 51878605 and No. 41807224) and the Zhejiang Provincial National Science Foundation
651
of China (LQ17D020001).
652
Notation
653
ac , aG , aM , at
654
d1
dilatancy parameter
655
Dijrs
stiffness tensor
656
ds
stress increment
657
e , ec , e0
current, critical, and initial void ratios
model parameters
23
658
eij
deviator strain tensor
659
eΓ , eΓ 0
critical state parameters
660
Fb , Fc , Fd , Fm , Fp
661
Fb( mk ) , Fc( mk ) , Fd( mk ) , Fm( mk ) micro characteristic lines
662
Ge , Geq
663
( mk ) G ( mk ) , Ge( mk ) , Gp
664
G0 , G h0
shear modulus parameters
665
g (θˆ )
yield shape function
666
h
Heaviside step function
667
h1 , h2
model parameters
668
H ( mk ) , H p( mk )
micro total and plastic shear moduli
669
Iˆ1 , Iˆ2 , Iˆ3
stress invariants
670
K , Ke , Keq , K p
671
K1
variable
672
kt
model parameter
673
l, m , n
unit vectors
674
M b , M c , M d , M m bounding, critical, dilatancy, and maximum stress ratios
675
M c0
critical stress ratio for host materials
676
n b , nc , nd
model parameters
677
N ij( mk )
projection tensor
678
p , pˆ
mean effective stress and modified mean effective stress
679
pa
atmospheric pressure
680
pc , pc0
current and initial effective confining stresses
681
pˆ m , pˆ r
maximum and back values of modified mean effective stress
682
pp
pore pressure
683
pt0 , pt
initial and current values of hydrate bonding strength
macro characteristic surfaces
elastic and equivalent shear moduli micro total, elastic, and plastic shear moduli
total, elastic, equivalent, and plastic bulk moduli
24
684
q , qf
deviator stress and peak strength
685
Qij , Q ( mk )
variables
686
Rˆ
macro modified stress ratio
687
rˆij
macro modified stress ratio tensor
688
r ( mk )
micro stress ratio
689
rb( mk ) , rc( mk ) , rd( mk ) , rm( mk )
690
rr( mk )
back value of micro stress ratio
691
Sh
hydrate saturation
692
t
time
693
T
temperature
694
w(m)
weight coefficient
695
α
principal angle
696
δ
variation
697
δ ij
Kronecker delta
698
ε ij , ε ijp
total and plastic strain tensors
699
ε ijdis
collapse strain tensor caused by hydrate dissociation
700
ε qp
equivalent plastic shear strain
701
ε v , ε ve , ε vp
total, elastic, and plastic volumetric strains
702
ε vd
volumetric strain caused by dilatancy
703
ε vd( mk )
micro volumetric strain caused by dilatancy
704
ε a , ε r , ε θ , ε aθ
705
ε ap , ε θp , ε apθ
706
γ (mk ) , γ e( mk ) , γ p( mk )
707
κ , κ0
bulk modulus parameters related to isotropic unloading
708
λ , λ0
bulk modulus parameters related to isotropic loading
709
λc
critical state parameter
710
ν
micro bounding, critical, dilatancy, and maximum stress ratios
total strain components
plastic strain components micro total, elastic, and plastic shear strains
Poisson ratio
25
711
θˆ
712
ρ1( mk ) , ρ1( mk ) mapping intervals for micro stress ratio
713
ρ2 , ρ2
714
σ a , σ r , σ θ , σ aθ stress components
715
σ ij , σˆij
effective stress tensor and modified effective stress tensor
716
σ1 , σ 2 , σ3
effective principal stresses
717
σˆ1 , σˆ 2 , σˆ3
modified effective principal stresses
718
τ ( mk )
micro shear stress
719
ψ
state parameter
720
Ω
hemisphere surface
721
Lode angle
mapping intervals for modified mean effective stress
Macaulay bracket
722 723
Appendix. Constitutive relation
724
A.1. Multishear model
725 726
By projecting the tensor of the macro deviator strain onto the microcontact, the strain components in the microshear structure of the n direction can be written as 1 2
1 2
727
γ (1) (n ) = (li n j + l j ni )eij = (li n j + l j ni )(ε ij − ε vδ ij / 3) = N ij(1) (n )ε ij
728
γ (2) (n ) = ( mi n j + m j ni )eij = ( mi n j + m j ni )(ε ij − ε vδ ij / 3) = N ij(2) (n )ε ij
729
1 2
1 2
γ (3) (n) = ni n j eij = ni n j (ε ij − ε vδ ij / 3) = N ij(3) (n )ε ij
(A1) (A2) (A3)
730
where γ (1) (n ) and γ (2) (n ) denote the tangential strain components corresponding to the l and
731
m directions, respectively; γ (3) (n ) represents the deviatoric normal strain corresponding to the
732
n direction; and N ij(1) (n ) , N ij(2) (n ) , and N ij(3) (n ) are expressed as
733 734 735
N ij(1) (n ) =
1 1 1 (li n j + l j ni ) , N ij(2) (n ) = (mi n j + m j ni ) , N ij(3) (n ) = ni n j − δ ij 2 2 3
(A4)
According to the virtual work principle, the stress vectors on all microshear structures and the macro stress tensor are variationally enforced, which is written as
26
736 737 738 739 740 741 742
3 4π 4π 1 dσ ij δε ij = dpδε v + 2∫ ∑ dτ ( k ) (n)δγ ( k ) (n)dΩ Ω 3 3 k =1 3
(A5)
where δ denotes a variation; and Ω represents the surface of hemisphere. By substituting δε v = δ ij δε ij and δγ ( k ) (n) = N ij( k ) (n)δε ij into equation (A5) and holding the variational equation for any variation δ ε ij , the macro effective stress increment is given as
dσ ij = dpδij +
3 M 3 1 (k ) (k ) N ( n )d τ ( n )d Ω ≈ d p δ + 2w(m) Nij(mk )dτ (mk ) ∑ ∑∑ ij ij ∫ 2π Ω k =1 m=1 k =1
(A6)
A.2. Stress–strain relation on macro volumetric deformation The increment of the volumetric strain is the sum of the elastic and plastic parts, i.e., dε v =dε ve +dε vp
743
(A7)
744
where ε ve refers to the elastic volumetric strain; and ε vp denotes the plastic volumetric strain that
745
involves the strain due to a variation in mean effective stress and the strain due to dilatancy.
746
Following the elasticity theory, dε ve is expressed as
748
751 752 753 754 755
756 757
(A8)
Following the plasticity theory of bounding surface, dε vp is given by
dεvp =
749 750
1 dpˆ Ke
dε ve =
747
dpˆ 1 h( pˆ − pˆm ) dpˆ + dεvd Kp dpˆ
(A9)
Combing equations (A7)–(A9) yields dε v =
dpˆ 1 1 dpˆ + h( pˆ − pˆ m ) dpˆ + dε vd Ke Kp dpˆ
(A10)
Equation (A10) is rewritten as
dε v =
1 dpˆ +dε vd K
(A11)
where the total bulk modulus K is expressed as 1 dpˆ 1 K = + h ( pˆ − pˆ m ) K dpˆ e Kp
−1
A.3. Stress–strain relation on micro shear deformation The increment of micro shear strain is the sum of the elastic and plastic parts, i.e.,
27
(A12)
dγ ( mk ) =dγ e( mk ) +dγ p( mk )
758
(A13)
759
where γ e( mk ) denotes the kth component of the micro elastic shear strain in the mth microshear
760
structure.
761
Following the elasticity theory, dγ e( mk ) is expressed as
dγe(mk) =
762 763
dγ p(mk ) =
1 Gp(mk )
pˆ dr (mk ) +
765
Combing equations (A13)–(A15) yields
766
dγ ( mk ) =
768
1 ( mk ) e
G
dτ ( mk ) +
1 ( mk ) p
G
1
h pˆ − pˆ m ) ( mk ) (
Hp
pˆ dr ( mk ) +
1 H
( mk ) p
dpˆ (mk ) r dpˆ dpˆ
h( pˆ − pˆ m )
dpˆ (mk ) r dpˆ dpˆ
(A15)
(A16)
A.4. Macro stress–strain relation Following the definition of r (mk ) , dτ ( mk ) is expressed as
dτ
769 770
(A14)
Following the plasticity theory of bounding surface, dγ p( mk ) is given by
764
767
1 dτ (mk) G (mk) e
(mk )
τ (mk ) = d pˆ = pˆ dr (mk ) + r(mk )dpˆ pˆ
(A17)
Combing equations (A16) and (A17) yields
dτ (mk ) = G(mk )dγ (mk ) + (1 − G(mk ) / H (mk ) )r (mk )dpˆ
771
(A18)
772
where H (mk ) and G (mk ) represent the total micro shear moduli related to d r ( mk ) and dpˆ for the
773
kth component in the mth microshear structure, respectively, define as
774
H
775 776
G
( mk )
1 1 = ( mk ) + ( mk ) Ge Gp
3
dε vd = ∑∑ 2w( m ) d1 (± rd( mk ) − r ( mk ) ){(1 − m =1 k =1
778
1 dpˆ 1 = ( mk ) + ( mk ) h ( pˆ − pˆ m ) Ge Hp dpˆ
−1
(A19)
−1
(A20)
From equations (19), (20), (A13), (A14), and (A18), dε vd is expressed as M
777
( mk )
G ( mk ) 1 G ( mk ) ( mk ) ( mk ) )d γ − (1 − )r dpˆ } Ge( mk ) Ge( mk ) H ( mk )
Substituting equations (A21) into equation (A11) yields
28
(A21)
M
3
dpˆ =K1{dε v − ∑∑ 2 w( m ) d1 (± rd( mk ) − r ( mk ) )(1 − G ( mk ) / Ge( mk ) ) d γ ( mk ) }
779
(A22)
m =1 k =1
780
where K1 is expressed as
K1 =
781
K M
3
1 − K ∑∑ 2 w d1 (± r ( m)
( mk ) d
m =1 k =1
782 783
(A23) )r
( mk )
(1 − G
( mk )
/H
( mk )
( mk ) e
)/G
From equations (1), (2), (A18), and (A22), the macro constitutive relation, including the effect of hydrate dissociation, is written as
dσ ij = Dijrs (dε rs − dε rsdis ) − dptδij
784 785
−r
( mk )
(A24)
where Dijrs is expressed as M
3
Dijrs = K1Qijδ rs + ∑∑ 2 w( m ) (−Q ( mk )Qij + G ( mk ) N ij( mk ) ) N rs( mk )
786
(A25)
m =1 k =1
787
in which M
788
3
Qij = δ ij + ∑∑ 2 w( m ) r ( mk ) (1 − G ( mk ) / H ( mk ) ) N ij( mk )
(A26)
Q ( mk ) = K1d1 ( ± rd( mk ) − r ( mk ) )(1 − G ( mk ) / Ge( mk ) )
(A27)
m =1 k =1
789 790 791
A.5. Relations of micro and macro parameters
792
By assuming that rb( mk ) , rc( mk ) , and rd( mk ) are the same for all microshear structures, they can be
793
determined through the macro model constants and yielding shape function as (Fang et al., 2017),
794
given by
ri ( mk ) =
795
2 3
M i g (θˆ) M
(i=b, c, d)
3
∑∑ 2w
(m)
m =1 k =1
N
(A28)
( mk ) 33
796 797 798 799
By assuming that all microshear structures have the same elastic shear modulus, Ge(mk ) is given by (Fang et al., 2017)
Ge( mk ) =
4 3
Ge M
3
∑∑ 2w
(m)
m =1 k =1
800
29
N
( mk ) 2 33
(A29)
Data Availability Datasets related to this article can be found at http.
References Bažant, Z.P., Oh, B.H., 1985. Microplane model for progressive fracture of concrete and rock. J. Eng. Mech. 111(4), 559–582. Boswell, R., Schoderbek, D., Collett, T.S., Ohtsuki, S., White, M., Anderson, B.J., 2017. The Iġnik Sikumi Field Experiment, Alaska North Slope: Design, operations, and implications for CO2-CH4 exchange in gas hydrate reservoirs. Energy & Fuels 31(1), 140–153. Been, K., Jefferies, M.G., 1985. A state parameter for sands. Géotechnique 35(2), 99–112. Caner, F.C., Bažant, Z.P., 2013. Microplane model M7 for plain concrete. I: Formulation. J. Eng. Mech. 139 (12), 1714–1723. Choi, J.H., Dai, S., Lin, J.H., Seol, Y., 2018. Multistage triaxial tests on laboratory-formed methane hydrate-bearing sediments. J. Geophys. Res. Solid Earth 123, 3347–3357. Clayton, C.R.I., Priest, J.A., Best, A.I., 2005. The effects of disseminated methane hydrate on the dynamic stiffness and damping of a sand. Géotechnique 55(6), 423–434. Clayton, C.R.I., Priest, J., Rees, E., 2010. The effects of hydrate cement on the stiffness of some sands. Géotechnique 60(6), 435–45. Dai, S., Santamarina, J.C., Waite, W.F., Kneafsey, T.J., 2012. Hydrate morphology: Physical properties of sands with patchy hydrate saturation. J. Geophys. Res. 117, B11205. Ebinuma, T., Kamata, Y., Minagawa, W., Ohuma, R., Nagao, J., Narita, H., 2005. Mechanical properties of sandy sediment containing methane hydrate. In Proceedings of the 5th International Conference on Gas Hydrate (pp. 958–961). Trondheim, Norway. Fang, H.L., Zheng, H., Zheng, J., 2017. Micromechanics-based multi-mechanism bounding surface model for sands. Int. J. Plast. 90, 242–266. Fang, H., Shen, Y., and Zhao, Y., 2019. Multishear bounding surface modelling of anisotropic sands accounting for fabric and its evolution. Comput. Geotech. 110, 57–70. Gai, X., Sánchez, M., 2017. A geomechanical model for gas hydrate-bearing sediments. Environ. Geotech. 4(2), 143-156. Ghiassian, H., Grozic, J.L.H., 2013. Strength behavior of methane hydrate bearing sand in undrained triaxial testing. Mar. Petrol. Geol. 43, 310–319.
30
Golkhou, F., Haghtalab, A., 2019. Measurement and thermodynamic modeling of carbon dioxide hydrate formation conditions using dry water through hydrophobic nano silica. J. Nat. Gas Sci. Eng. 68, 102906. Hajirezaie, S., Hemmati-Sarapardeh, A., Mohammadi, A.H., Pournik, M., Kamari, A., 2015. A smooth model for the estimation of gas/vapor viscosity of hydrocarbon fluids. J. Nat. Gas Sci. Eng. 26, 1452–1459. Hemmat Esfe, M., Tatar, A., Ahangar, M., Rostamian, H., 2018. A comparison of performance of several artificial intelligence methods for predicting the dynamic viscosity of TiOC2/SAE 50 nano-lubricant. Phys. E 96, 85–93. Hyodo, M., Yoneda, J., Yoshimoto, N., Nakata, Y., 2013a. Mechanical and dissociation properties of methane hydrate-bearing sand in deep seabed. Soils Found. 53, 299–314. Hyodo, M., Li, Y., Yoneda, J., Nakata, Y., Yoshimoto, N., Nishimura, A., Song, Y., 2013b. Mechanical behavior of gas-saturated methane hydrate-bearing sediments. J. Geophys. Res. Solid Earth 118, 5185–5194. Hyodo, M., Li, Y., Yoneda, J., Nakata, Y., Yoshimoto, N., Nishimura, A., 2014. Effects of dissociation on the shear strength and deformation behavior of methane hydrate-bearing sediments. Mar. Petrol. Geol. 51, 52–62. Hyodo, M., Wu, Y., Nakashima, K., Kajiyama, S., Nakata, Y., 2017. Influence of fines content on the mechanical behavior of methane hydrate-bearing sediments. J. Geophys. Res. Solid Earth 122, 7511–7524. Kajiyama, S., Hyodo, M., Nakata, Y., Yoshimoto, N., Wu, Y., Kato, A., 2017a. Shear behaviour of methane hydrate bearing sand with various particle characteristics and fines. Soils Found. 57, 176–193. Kajiyama, S., Wu, Y., Hyodo, M., Nakata, Y., Nakashima, K., Yoshimoto, N., 2017b. Experimental investigation on the mechanical properties of methane hydrate-bearing sand formed with rounded particles. J. Nat. Gas Sci. Eng. 45, 96–107. Kamalgharibi, M., Hormozi, F., Zamzamian, S.A.H., Sarafraz, M.M., 2016. Experimental studies on the stability of CuO nanoparticles dispersed in different base fluids: Influence of stirring, sonication and surface active agents. Heat Mass Transfer 52(1), 55–62. Kimoto, S., Oka, F., Tomohiko, F., 2010. A chemo-thermo-mechanically coupled analysis of ground deformation induced by gas hydrate dissociation. Int. J. Mech. Sci. 52(2), 365–376.
31
Klar, A., Soga, K., Ng, M.Y.A., 2010. Coupled deformation-flow analysis for methane hydrate extraction. Géotechnique 60(10), 765–776. Lade, P.V., Duncan, J.M., 1975. Elasto-plastic stress–strain theory for cohesionless soil. J. Geotech. Eng. 101(10), 1037–1053. Li, X.S., 2002. A sand model with state-dependent dilatancy. Géotechnique 52(3), 173–186. Lin, J.S., Seol, Y., Choi, J.H., 2015. An SMP critical state model for methane hydrate-bearing sands. Int. J. Numer. Anal. Meth. Geomech. 39(9), 969–987. Lin, J.S., Seol, Y., Choi, J.H., 2017. Geomechanical modeling of hydrate-bearing sediments during dissociation under shear. Int. J. Numer. Anal. Meth. Geomech. 41, 1523–1538. Liu, W., Zhao, T, Luo, Y., Li, Y., Song, Y., Zhu, Y., Liu, Y., Zhao, J., Wu, Z., Xu, X., 2016. Experimental study on the mechanical properties of sediments containing CH4 and CO2 hydrate mixtures. J. Nat. Gas Sci. Eng. 32, 20–27. Liu, Z., Dai, S., Ning, F., Peng, L., Wei, H., Wei, C., 2017. Strength estimation for hydrate-bearing sediments from direct shear tests of hydrate-bearing sand and silt. Geophys. Res. Lett. 45, 715– 723. Malagar, B.R.C., Lijith, K.P., Singh, D.N., 2019. Formation & dissociation of methane gas hydrates in sediments: A critical review. J. Nat. Gas Sci. Eng. 65, 168–184. Masui, A., Haneda, H., Ogata, Y., Aoki, K., 2005. The effects of methane hydrate formation on shear strength of synthetic methane hydrate sediments. In Proceedings of the 5th International Conference on Gas Hydrates (pp. 657–663). Trondheim, Norway. Masui, A., Miyazaki, K., Haneda, H., Ogata, Y., Aoki, K., 2008. Mechanical characteristics of natural and artificial gas hydrate bearing sediments. In Proceedings of the 6th International Conference on Gas Hydrates. Chevron, Vancouver, B. C., Canada. Matsuoka, H., Nakai, T., 1974. Stress-deformation and strength characteristics of soil under three difference principal stresses. Proc. JSCE 232, 59–70. Matsuoka, H., Yao, Y.P., Sun, D.A., 1999. The Cam-clay models revised by the SMP criterion. Soils Found. 39(1), 81–95. Miyazaki, K., Masui, A., Sakamoto, Y., Aoki, K., Tenma, N., Yamaguchi, T., 2011. Triaxial compressive properties of artificial methane-hydrate-bearing sediment. J. Geophys. Res. 116, B06102.
32
Miyazaki, K., Tenma, N., Aoki, K., Yamaguchi, T., 2012. A nonlinear elastic model for triaxial compressive properties of artificial methane-hydrate-bearing sediment samples. Energies 5(10), 4057–4075. Nikkhah, V., Sarafraz, M.M., Hormozi, F., 2015a. Application of spherical copper oxide (II) water nano-fluid as a potential coolant in a boiling annular heat exchanger. Chem. Biochem. Eng. Q. 29(3), 405–415. Nikkhah, V., Sarafraz, M.M., Hormozi, F., Peyghambarzadeh, S.M., 2015b. Particulate fouling of CuO-water nanofluid at isothermal diffusive condition inside the conventional heat exchanger-experimental and modeling. Exp. Therm. Fluid Sci. 60, 83–95. Nguyen, L.D., Fatahi, B., Khabbaz, H., 2014. A constitutive model for cemented clays capturing cementation degradation. Int. J. Plast. 56, 1–18. Pinkert, S., Grozic, J., 2014. Prediction of the mechanical response of hydrate-bearing sands. J. Geophys. Res. Solid Earth 119(6), 4695–4707. Pinkert, S., Grozic, J., Priest, J., 2015. Strain-softening model for hydrate-bearing sands. Int. J. Geomech. 15(6), 04015007. Priest, J.A., Best, A.I., Clayton, C.R.I., 2005. A laboratory investigation into the seismic velocities of methane gas hydrate-bearing sand. J. Geophys. Res. 110, B04102. Priest, J.A., Rees, E.V.L., Clayton, C.R.I., 2009. Influence of gas hydrate morphology on the seismic velocities of sands. J. Geophys. Res. Solid Earth 114, B11205. Richart, F., Hall, J, Woods, R., 1970. Vibrations of soils and foundations. International Series in Theoretical and Applied Mechanics. Englewood Cliffs, NJ, USA: Prentice-Hall. Rostamian, H., Lotfollahi, M.N., 2016. New functionality for energy parameter of redlich-kwong equation of state for density calculation of pure carbon dioxide and ethane in liquid, vapor and supercritical phases. Period. Polytech., Chem. Eng. 60(2), 93–97. Rostamian, H., Lotfollahi, M.N., 2019. A novel statistical approach for prediction of thermal conductivity of CO2 by response surface methodology. Phys. A 527, 121175. Sánchez, M., Gai, X., Santamarina, J.C., 2017. A constitutive mechanical model for gas hydrate bearing sediments incorporating inelastic mechanisms. Comput. Geotech. 84, 28–46. Sánchez, M., Santamarina, C., Teymouri, M., Gai, X., 2018. Coupled numerical modeling of gas hydrate-bearing sediments: From laboratory to field-scale analyses. J. Geophys. Res. Solid Earth 123(10), 326–348.
33
Santamarina, J.C., Dai, S., Terzariol, M., Jang, J., Waite, W.F., Winters, W.J., Suzuki, K., 2015. Hydro-bio-geomechanical properties of hydrate-bearing sediments from Nankai Trough. Mar. Petrol. Geol. 66, 434–450. Sasaki, T., Soga, K., Elshafie, M.E.B., 2018. Simulation of wellbore construction in offshore unconsolidated methane hydrate-bearing formation. J. Nat. Gas Sci. Eng. 60, 312–326. Shen, J., Chiu, C.F., Ng, C.W.W., Lei, G.H., Xu, J., 2016. A state-dependent critical state model for methane hydrate-bearing sand. Comput. Geotech. 75, 1–11. Song, B., Cheng, Y., Yan, C., Ly, Y., Wei, J., Ding, J., Li, Y., 2019. Seafloor subsidence response and submarine slope stability evaluation in response to hydrate dissociation. J. Nat. Gas Sci. Eng. 65, 197–211. Sultan N., Garziglia S., 2011. Geomechanical constitutive modelling of gas-hydrate bearing sediments. In Proceedings of the 7th International Conference on Gas Hydrates. Edinburgh, Scotland, United Kingdom. Sultaniya, A., Priest, J.A., Clayton, C.R.I., 2018. Impact of formation and dissociation conditions on stiffness of a hydrate-bearing sand. Can. Geotech. J. 55, 988–998. Sun, D., Huang, W., Yao, Y., 2008. An experimental study of failure and softening in sand under three-dimensional stress condition. Granular Matter 10, 187–195. Sun, X., Guo, X., Shao, L., Tang, H., 2015. A thermodynamics-based critical state constitutive model for methane hydrate bearing sediment. J. Nat. Gas Sci. Eng. 27, 1024–1034. Taylor, G.I., 1938. The plastic strain of metals. J. Inst. Metals 62, 307–324. Thevanayagam, S., Shenthan, T., Mohan, S., Liang, J., 2002. Undrained fragility of clean sands, silty sands, and sandy silts. J. Geotech. Geoenviron. Eng. 128(10), 849–859. Uchida, S., Soga, K., Yamamoto, K., 2012. Critical state soil constitutive model for methane hydrate soil. J. Geophys. Res. 117, B03209. Uchida, S., Xie, X.G., Leung, Y.F., 2016. Role of critical state framework in understanding geomechanical behavior of methane hydrate-bearing sediments. J. Geophys. Res. Solid Earth 121(8), 5580–5595. Wang, Z.L., Dafalias, Y.F., Shen, C.K., 1990. Bounding surface hypoplasticity model for sand. J. Eng. Mech. 116(5), 983–1001. Winters, W.J., Waite, W.F., Mason, D.H., Gilbert, L.Y., Pecher, I.A., 2007. Methane gas hydrate effect on sediment acoustic and strength properties. J. Petrol. Sci. Eng. 56(1–3), 127–135.
34
Yan, C., Li, Y., Cheng, Y., Wang, W., Song, B., Deng, F., Feng, Y., 2018. Sand production evaluation during gas production from natural gas hydrates. J. Nat. Gas Sci. Eng. 57, 77–88. Yan, G., Cheng, Y., Li, M., Han, Z., Zhang, H., Li, Q., Teng, F., Ding, J., 2017. Mechanical experiments and constitutive model of natural gas hydrate reservoirs. Int. J. Hydrogen Energy 42, 19810–19818. Yan, Y., Wei, C., 2017. Constitutive model for gas hydrate-bearing soils considering hydrate occurrence habits. Int. J. Geomech. 17(8), 04017032. Yao, Y., Sun, D., 2000. Application of Lade’s criterion to Cam-clay model. J. Eng. Mech. 126(1), 112–119. Yoneda, J., Masui, A., Konno, Y., Jin, Y., Egawa, K., Kidab, M., Ito, T., Nagao, J., Tenma, N., 2015a. Mechanical behavior of hydrate-bearing pressure-core sediments visualized under triaxial compression. Mar. Petrol. Geol. 66(2), 451–459. Yoneda, J., Masui, A., Konno, Y., Jin, Y., Egawa, K., Kida, M., Ito, T., Nagao, J., Tenma, N., 2015b. Mechanical properties of hydrate-bearing turbidite reservoir in the first gas production test site of the Eastern Nankai Trough. Mar. Petrol. Geol. 66(2), 471–486. Yoneda, J., Jin, Y., Katagiri, J., Tenma, N., 2016. Strengthening mechanism of cemented hydrate-bearing sand at microscales. Geophys. Res. Lett. 43(14), 7442–7450. Yoneda, J., Oshima, M., Kida, M., Kato, A., Konno, Y., Jin, Y., Jang, J., Waite, W.F., Kumar, P., Tenma, N., 2019a. Pressure core based onshore laboratory analysis on mechanical properties of hydrate-bearing sediments recovered during India's National Gas Hydrate Program Expedition (NGHP) 02. Mar. Petrol. Geol. (In press) Yoneda, J., Oshima, M., Kida, M., Kato, A., Konno, Y., Jin, Y., Tenma, N., 2019a. Consolidation and hardening behavior of hydrate-bearing pressure-core sediments recovered from the Krishna–Godavari Basin, offshore India. Mar. Petrol. Geol. (In press) Yoneda, J., Oshima, M., Kida, M., Kato, A., Konno, Y., Jin, Y., Jang, J., Waite, W.F., Kumar, P., Tenma, N., 2019b. Pressure core based onshore laboratory analysis on mechanical properties of hydrate-bearing sediments recovered during India's National Gas Hydrate Program Expedition (NGHP) 02. Mar. Petrol. Geol. (In press) Yu, T., Guan, G., Abudula, A., Yoshida, A., Wang, D., Song, Y., 2019. Application of horizontal wells to the oceanic methane hydrate production in the Nankai Trough, Japan. J. Nat. Gas Sci. Eng. 62, 113–131.
35
Yun, T.S., Santamarina, C.J., Ruppel, C., 2007. Mechanical properties of sand, silt, and clay containing tetrahydrofuran hydrate. J. Geophys. Res. 112, B04106. Zhang, X., Lin, J., Lu, X., Liu, L., Liu, C., Li, M., Su, Y., 2017. A hypoplastic model for gas hydrate-bearing sandy sediments. Int. J. Numer. Anal. Meth. Geomech. 42, 931–942. Zhao, J., Guo, N., 2013. Unique critical state characteristics in granular media considering fabric anisotropy. Géotechnique 63(8), 695–704. Zhou, M., Soga, K., Yamamoto, K., 2018. Upscaled anisotropic methane hydrate critical state model for turbidite hydrate-bearing sediments at East Nankai Trough. J. Geophys. Res. Solid Earth 123(8), 6277–6298.
36
Table 1 Model parameters Test name Triaxial compression tests for synthetic pore-filling type of GHBS (Masui et al., 2005)
Parameters for host materials G0 =150, κ0 =0.005
Triaxial compression tests for synthetic cementing type of GHBS (Masui et al., 2005)
G0 =150, κ0 =0.005
Triaxial compression tests for synthetic GHBS (Hyodo et al., 2013a)
G0 =200, κ0 =0.004
M c0 =1.47, eΓ 0 =1.126, λ c =0.132 h1 =0.37, h2 =0.35, nb =1.5, λ0 =0.01 d1 =0.7, nd =2.3 M c0 =1.47, eΓ 0 =1.126, λ c =0.132 h1 =0.66, h2 =0.35, nb =1.5, λ0 =0.01
Additional parameters for GHBS aG =3.5, a M =0.038 ac =0.04, nc =0.53 for S h < 0.3 ac =0.39, nc =0.53 for S h ≥ 0.3 a t =0.15MPa, k t =50
aG =2.66, a M =0.071 ac =0.6, nc =2.0 a t =0.7MPa, k t =50
d1 =0.7, nd =4.0 M c0 =1.20, eΓ 0 =1.126, λ c =0.132 h1 =0.55, h2 =0.35, nb =0.8, λ0 =0.008
aG =2.66, a M =0.226 ac =0.6, nc =2.0 a t =0.43MPa, k t =50
d1 =0.6, nd =3.5 Triaxial compression tests for natural GHBS (Masui et al., 2008)
G0 =200, κ0 =0.004 M c0 =1.48, eΓ 0 =1.126, λ c =0.132 h1 =0.55, h2 =0.35, nb =1.0, λ0 =0.008
aG =2.66, a M =0.071 ac =0.37, nc =0.0 a t =0.6MPa, k t =50
d1 =0.7, nd =1.5 Triaxial compression tests for natural GHBS (Yoneda et al., 2015b)
G0 =100, κ0 =0.0075 M c0 =1.44, eΓ 0 =1.126, λ c =0.132 h1 =0.55, h2 =0.35, nb =3.5, λ0 =0.015
aG =1.5, a M =0.157 ac =0.12, nc =0.3 a t =0.6MPa, k t =80
d1 =0.7, nd =3.0 Triaxial compression tests for natural GHBS (Yoneda et al., 2019a, 2019b )
G0 =100, κ0 =0.0075 M c0 =1.2, eΓ 0 =0.8, λ c =0.12 h1 =0.2, h2 =0.35, nb =1.0, λ0 =0.0225 d1 =2.0, nd =1.0
37
aG =1.4, a M =0.125 ac =452.0, nc =13.5 for S h < 0.6 ac =0.13, nc =-3.4 for S h ≥ 0.6 a t =1.7MPa, k t =100
Fig. 1 Illustration of conceptual framework for multishear model
Fig. 2 Schematic diagram of different morphologies for GHBS
38
1.2
Solid lines: regressions Test (Sh=0) Test (Sh=0.242-0.287) Test (Sh=0.311-0.387) Sh=0.439 (average) Test (Sh=0.419-0.450) ec=0.939-0.132lnp 2 Test (Sh=0.504-0.572) R =0.75
Sh=0.532 (average) ec=0.993-0.132lnp 2
R =0.23
Void ratio, e
1.0
0.8 Sh=0 ec=0.823-0.132lnp Sh=0.262 (average) ec=0.865-0.132lnp
2
0.6
R =0.76
2
R =0.55
Sh=0.342 (average) ec=0.894-0.132lnp 2
R =0.68 0.4 1
2
3
4
5
6
7
8
9 10
Mean effective stress, p (MPa)
Fig. 3 Critical state line in e − ln p plane (test data from Hyodo et al., 2013a)
2.0
Critical state void ratio parameter, eΓ
Regression Obtained by fitting lines in Fig.3
1.5
1.0
2
eΓ=1.126+0.6Sh 2
R =0.85 0.5
0.0 0.0
0.2
0.4
0.6
0.8
Hydrate saturation, Sh
Fig. 4 Relationship between critical state void ratio parameter and hydrate saturation (test data from Hyodo et al., 2013a)
39
1.6
Regression Test (pp=10MPa, T= 5 Test (pp= 5MPa, T= 5 Test (pp=15MPa, T=10 Test (pp=15MPa, T= 1
Critical stress ratio, Mc
1.4
) ) ) )
1.2 2
Mc=1.2(1+0.226Sh), R =0.32 1.0
0.8 0.0
0.2
0.4
0.6
0.8
Hydrate saturation, Sh
Fig. 5 Relationship between critical stress ratio and hydrate saturation (test data from Hyodo et al., 2013a)
15
Solid lines: regressions Test (Sh=0.242-0.296) Test (Sh=0.414-0.457) Test (Sh=0.504-0.572)
Deviator stress, q (MPa)
12
Sh=0.532 (average) qf=1.42p+0.35
9
Sh=0.432 (average) qf=1.38p+0.26
2
R =0.95
2
R =0.99
6 Sh=0.267 (average) qf=1.29p+0.11 3
2
R =0.92
0 -2
0
2
4
6
8
10
12
Mean effective stress, p (MPa)
Fig. 6 Peak strength lines of GHBS with different hydrate saturations (test data from Hyodo et al., 2013a)
40
0.4
Initial bonding strength, pt0 (MPa)
Regression Obtained by fitting lines in Fig. 6
0.3
2
pt0=0.43Sh, R =0.96 0.2
0.1
0.0 0.0
0.2
0.4
0.6
0.8
Hydrate saturation, Sh
Fig. 7 Relationship between initial bonding strength and hydrate saturation (test data from Hyodo et al., 2013a)
Fig. 8 Characteristic surfaces in macro stress space
41
Fig. 9 Characteristic lines in τ ( mk ) − pˆ plane
Elastic shear modulus parameter, Gh0
800 Regression Test
600
Gh0=200(1+2.66Sh) 2
R =0.61 400
200
0 0.0
0.2
0.4
0.6
0.8
Hydrate saturation, Sh
Fig. 10 Relationship between elastic shear modulus parameter and hydrate saturation (test data from Hyodo et al., 2013a)
42
10 Simulation Test Sh e0 None 0.678 0.736 None 0.501 0.736 0.409 0.736 None 0.264 0.736
Deviator stress, q (MPa)
8
pc0 (MPa) 1.0 1.0 1.0 1.0
6
4
2
Pore pressure: 8 MPa Temperature: 5
0 0
2
4
6
8
10
12
14
16
Axial strain, εa (%)
-10
Simulation Test Sh e0 None 0.678 0.736 None 0.501 0.736 0.409 0.736 None 0.264 0.736
Volumetric strain, εv (%)
-8
-6
pc0 (MPa) 1.0 1.0 1.0 1.0
Pore pressure: 8 MPa Temperature: 5
-4
-2
0
2 0
2
4
6
8
10
12
14
16
Axial strain, εa (%)
Fig. 11 Test and simulation results for triaxial compression tests of synthetic pore-filling type of GHBS (test data from Masui et al., 2005)
43
10
Simulation Test Sh None 0.551 0.407 None 0.257 0.000
Deviator stress, q (MPa)
8
e0 0.605 0.605 0.605 0.605
pc0 (MPa) 1.0 1.0 1.0 1.0
6
4
2
Pore pressure: 8 MPa Temperature: 5
0 0
2
4
6
8
10
12
14
16
Axial strain, εa (%)
-16
Simulation Test Sh None 0.551 0.407 None 0.257 0.000
Volumetric strain, εv (%)
-12
e0 0.605 0.605 0.605 0.605
pc0 (MPa) 1.0 1.0 1.0 1.0
-8
-4
0
Pore pressure: 8 MPa Temperature: 5
4 0
2
4
6
8
10
12
14
16
Axial strain, εa (%)
Fig. 12 Test and simulation results for triaxial compression tests of synthetic cementing type of GHBS (test data from Masui et al., 2005)
44
16
Deviator stress, q (MPa)
14 12 10 8 6 Simulation Test 4 Pore pressure: 10 MPa Temperature: 5
2
Sh 0.531 0.351 0.242 0.000
e0 0.669 0.645 0.565 0.650
pc0 (MPa) 5.0 5.0 5.0 5.0
0 0
2
4
6
8
10
12
14
16
12
14
16
Axial strain, εa (%)
-8 Simulation Test
Volumetric strain, εv (%)
-6 -4
Sh 0.531 0.351 0.242 0.000
e0 0.669 0.645 0.565 0.650
pc0 (MPa) 5.0 5.0 5.0 5.0
-2 0 2 4 Pore pressure: 10 MPa Temperature: 5
6 8 0
2
4
6
8
10
Axial strain, εa (%)
Fig. 13 Test and simulation results for triaxial compression tests of synthetic GHBS (test data from Hyodo et al., 2013a)
45
Simulation Test
20
e0
0.531 0.537 0.543 0.523 0.516
16
Deviator stress, q (MPa)
Sh
pc0 pp T (MPa) (MPa) ( ) 0.669 5.0 10.0 5.0 0.650 3.0 10.0 5.0 0.650 1.0 10.0 5.0 0.645 3.0 15.0 1.0 0.653 3.0 15.0 10.0
12
8
4
0 0
2
4
6
8
10
12
14
16
12
14
16
Axial strain, εa (%)
-10
Simulation Test
e0
0.531 0.537 0.543 0.523 0.516
-8
Volumetric strain, εv (%)
Sh
-6
pc0 pp T (MPa) (MPa) ( ) 0.669 5.0 10.0 5.0 0.650 3.0 10.0 5.0 0.650 1.0 10.0 5.0 0.645 3.0 15.0 1.0 0.653 3.0 15.0 10.0
-4
-2
0
2
4 0
2
4
6
8
10
Axial strain, εa (%)
Fig. 14 Test and simulation results for triaxial compression tests of synthetic GHBS (test data from Hyodo et al., 2013a)
46
7
Deviator stress, q (MPa)
6
5
4
3
2
Simulation Test Pore pressure: 9 MPa Temperature: 5
1
Sh 0.376 0.225 0.094 0.077
e0 0.792 0.961 0.972 1.033
pc0 (MPa) 1.0 1.0 1.0 1.0
0 0
2
4
6
8
10
12
14
16
12
14
16
Axial strain, εa (%)
-8
Simulation Test
Volumetric strain, εv (%)
-6
Sh 0.376 0.225 0.094 0.077
pc0 (MPa) 1.0 1.0 1.0 1.0
e0 0.792 0.961 0.972 1.033
Pore pressure: 9 MPa Temperature: 5
-4
-2
0
2 0
2
4
6
8
10
Axial strain, εa (%)
Fig. 15 Test and simulation results for triaxial compression tests of natural GHBS (test data from Masui et al., 2008)
47
14 Simulation Test Sample No. 9 No. 7 No. 9-2 No. 7-2
Deviator stress, q (MPa)
12
10
Sh 0.79 0.38 0.00 0.00
e0 pc0 (MPa) 0.650 1.6 0.790 1.5 0.748 1.6 0.887 1.5
8
6
4
2
Pore pressure: 13 MPa Temperature: 10
0 0
2
4
6
8
10
12
14
16
12
14
16
Axial strain, εa (%)
-9 Simulation Test Sample No. 9 No. 7 No. 9-2 No. 7-2
Volumetric strain, εv (%)
-6
Sh 0.79 0.38 0.00 0.00
e0 pc0 (MPa) 0.650 1.6 0.790 1.5 0.748 1.6 0.887 1.5
-3
0
3
6
Pore pressure: 13 MPa Temperature: 10
9 0
2
4
6
8
10
Axial strain, εa (%)
Fig. 16 Test and simulation results for triaxial compression tests of natural GHBS (test data from Yoneda et al., 2015b)
48
0.9
Simulation Test Sample Depth (m) Sh 16B-3P 273.81-273.89 0.321 16B-4P 277.22-277.30 0.716 16B-3P 274.13-274.23 0.601
0.8
Before dissociation
Void ratio, e
0.7
During dissociation 0.6
After dissociation
0.5
0.4
0.3 0.1
1
10
100
Effective confining stress, pc (MPa)
(a) Hydrate dissociation under effective confining stress of 1 MPa 0.9
0.8 Before dissociation
Void ratio, e
0.7
0.6
During dissociation After dissociation
0.5
0.4
0.3 0.1
Simulation Test Sample Depth (m) Sh 16B-7P 287.12-277.20 0.639 16B-3P 274.03-274.13 0.606 23C-11P 285.25-285.33 0.556 1
10
100
Effective confining stress, pc (MPa)
(b) Hydrate dissociation under effective confining stress of 5 MPa Fig. 17 Test and simulation results for isotropic consolidation tests of natural GHBS (test data from Yoneda et al., 2019a)
49
5
Pore pressure: 10 MPa Temperature: 4
Deviator stress, q (MPa)
4
3
2 Simulation Test 1
Sh 0.822 0.734 0.634 0.552 0.487
e0 0.656 0.658 0.825 0.634 0.712
pc0 (MPa) 1.1 1.0 1.0 1.0 1.0
0 0
2
4
6
8
10
12
14
16
12
14
16
Axial strain, εa (%)
-12
Simulation Test
Sh 0.822 0.734 None 0.634 0.552 0.487
Volumetric strain, εv (%)
-8
e0 pc0 (MPa) 0.656 1.1 0.658 1.0 0.825 1.0 0.634 1.0 0.712 1.0
-4
0
4
8
Pore pressure: 10 MPa Temperature: 4
12 0
2
4
6
8
10
Axial strain, εa (%)
Fig. 18 Test and simulation results for triaxial compression tests of natural GHBS (test data from Yoneda et al., 2019b)
50
εa
1.5
εθ
εaθ
εr Sh=0.0 Sh=0.3 Sh=0.6
Strains, εa, εθ, εaθ, εr (%)
1.0
0.5
0.0
-0.5 e0=0.65 pc0=1MPa
-1.0
-1.5 0
30
60
90
120
150
180
Principal stress direction, α (°)
Fig. 19 Predicted strain responses with principal stress direction for drained pure rotational shear tests
51
1.0
e0=0.65 pc0=1MPa
Strain increment Sh=0.0 Sh=0.3 Sh=0.6
σaθ (MPa), 2dεaθ/ds (%)
0.5
0.0
-0.5
1% -1.0 -1.0
-0.5
0.0
0.5
1.0
(σa-σθ)/2 (MPa), (dεa-dεθ)/ds (%)
Fig. 20 Predicted stress paths and strain increment vectors for drained pure rotational shear tests
1.0
e0=0.65 pc0=1MPa
Plastic strain increment Sh=0.0 Sh=0.3 Sh=0.6
p
σaθ (MPa), 2dεaθ/ds (%)
0.5
0.0
-0.5
0.2% -1.0 -1.0
-0.5
0.0
0.5
1.0
(σa-σθ)/2 (MPa), (dε -dεθ )/ds (%) p a
p
Fig. 21 Predicted stress paths and plastic strain increment vectors for drained pure rotational shear tests
52
Highlights •
A geomechanical constitutive model in the multishear bounding surface framework is proposed for gas hydrate-bearing sediments.
•
A macro constitutive response is decomposed into a macro volume response and a series of micro shear responses in spatial distributions.
•
The model captures the consolidation, hardening, softening, dilatation, collapse, and non-coaxial behavior of gas hydrate-bearing sediments.
1
Declaration of interests ☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. ☒The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:
None
S1875-5100(19)30371-3
DOI:
https://doi.org/10.1016/j.jngse.2019.103119
Reference:
JNGSE 103119
To appear in:
Journal of Natural Gas Science and Engineering
Received Date: 14 September 2019 Revised Date:
15 November 2019
Accepted Date: 18 December 2019
Please cite this article as: Fang, H., Shi, K., Yu, Y., Geomechanical constitutive modelling of gas hydrate-bearing sediments by a state-dependent multishear bounding surface model, Journal of Natural Gas Science & Engineering, https://doi.org/10.1016/j.jngse.2019.103119. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. 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. © 2019 Published by Elsevier B.V.
1
Title of paper
2
Geomechanical constitutive modelling of gas hydrate-bearing sediments by a state-dependent
3
multishear bounding surface model
4 5
Author 1
6
Huolang Fang
7
College of Civil Engineering and Architecture, Zhejiang University, 866 Yuhangtang Road,
8
Hangzhou 310058, China
9
Email: [email protected]
10 11
Author 2
12
Kenan Shi
13
College of Civil Engineering and Architecture, Zhejiang University, 866 Yuhangtang Road,
14
Hangzhou 310058, China
15
Email: [email protected]
16 17
Author 3
18
Yang Yu (Corresponding author)
19
Ocean College, Zhejiang University, 1 Zheda Road, Dinghai District, Zhoushan, Zhejiang 316021，
20
China
21
Email: [email protected]
1
22
Abstract
23
Natural gas hydrate in marine sediments and permafrost areas is considered as an important
24
potential energy source. Since hydrate dissociation will reduce the stability of gas hydrate-bearing
25
sediments (GHBS) and may cause wellbore failures and geological disasters during gas production,
26
it is necessary to reveal the mechanical behavior of GHBS for the safe exploitation of natural gas
27
hydrate. This paper proposes a geomechanical constitutive model of GHBS within the multishear
28
bounding surface framework. Following the slip theory of plasticity, a constitutive formulation is
29
obtained by splitting the macro constitutive response of sediments into a macro volume response
30
and a series of micro shear responses in spatial distributions related to virtual microshear structures.
31
Each microshear structure describes micro shear and dilatancy responses in three orthogonal
32
orientations. A micro stress–strain relationship and a micro stress–dilatancy relationship are
33
established for each orientation of the microshear structure. The model comprehensively describes
34
the consolidation, hardening, softening, dilatation, collapse, and non-coaxial characteristics of gas
35
hydrate-bearing sediments by introducing the multishear concept, state parameter, evolution law of
36
hydrate bonding and debonding, and collapse strain caused by hydrate dissociation. The
37
effectiveness of the model is confirmed by simulating the available published laboratory tests on
38
the samples of synthetic and natural GHBS under different pore pressures, temperatures, initial void
39
ratios, hydrate saturations, and initial effective confining stresses.
40
Keywords
41
gas hydrate, mechanical behavior, constitutive model, dilatancy, collapse deformation
2
42
1. Introduction
43
Natural gas hydrate is the crystalline solid where gas molecules are stored in ice-like crystal
44
lattices and is stable under high pressure and low temperature in marine sediments and permafrost
45
areas. The amount of natural gas trapped inside hydrate is roughly evaluated to be twice the existing
46
reserves of conventional fossil fuels in the earth. As a significant potential energy resource, the
47
exploitation of gas hydrate has attracted wide attention. At present, three possible ways to produce
48
gas from gas hydrate-bearing sediments (GHBS) are decompression, heat injection, and chemical
49
stimulation. Gas hydrate is usually located in unconsolidated sediments with low shear strength.
50
Hydrate dissociation will reduce the stability of GHBS due to hydrate exploitation. Particularly,
51
hydrate dissociation in the deep-sea marine sediments may result in wellbore failures and geological
52
disasters, such as casing damage, seafloor deformation, and submarine landslide. Therefore, it is
53
necessary to reveal the mechanical behavior of GHBS for the safe exploitation of natural gas
54
hydrate.
55
Over the past two decades, many experimental researches in the mechanical behavior of synthetic
56
GHBS have been performed (Ebinuma et al., 2005; Masui et al., 2005; Priest et al., 2005, 2009;
57
Clayton et al., 2005, 2010; Yun et al., 2007; Miyazaki et al., 2011; Ghiassian and Grozic, 2013;
58
Hyodo et al., 2013a, 2013b, 2014, 2017; Santamarina et al., 2015; Liu et al., 2016; Yoneda et al.,
59
2016; Kajiyama et al., 2017a, 2017b; Liu et al., 2017; Sultaniya et al., 2018; and Choi et al., 2018).
60
For example, Masui et al. (2005) performed several triaxial compression tests for samples of
61
synthetic GHBS involving different hydrate saturations, and concluded that the stiffness and
62
strength of samples were proportional to the saturation of hydrate. Priest et al. (2005, 2009) studied
63
the impacts of hydrate formation on the wave velocity in samples, revealing that the contact of
64
hydrate-cemented particles significantly enhanced the wave velocity. Miyazaki et al. (2011)
65
performed a set of creep tests in the triaxial apparatus for samples of synthetic GHBS under loading
66
and unloading at a given strain rate, and found that the samples first contracted and then dilated
67
owing to the existence of hydrate. Hyodo et al. (2013a, 2013b, 2014) conducted some triaxial
68
compression tests for the synthetic GHBS with Toyoura sands, and systematically investigated the
69
impacts of initial effective confining stress, pore pressure, hydrate saturation, temperature, and
70
density on the strength and deformation characteristics of GHBS. Kajiyama et al. (2017a, 2017b)
71
carried out several triaxial compression tests for the synthetic GHBS made with glass beads and
72
Toyoura sands, and studied the impact of grain properties on the strain and stress responses of
73
GHBS as well as the bonding and debonding mechanism of hydrate between particles. Hyodo et al. 3
74
(2017) conducted several triaxial compression tests for studying impacts of fines content as well as
75
density on the shear response of hydrate-free sediments and GHBS, and concluded that the increase
76
of fines content significantly enhanced the peak strength and dilative behavior. Choi et al. (2018)
77
experimentally investigated the mechanical response of synthetic GHBS through the triaxial tests of
78
single-stage and multi-stage, focusing on the strength, stiffness, and dilatancy of GHBS during
79
shearing.
80
In addition to the above experimental studies on synthetic GHBS, only a few experiments have
81
been conducted for samples of natural GHBS extracted from gas hydrate deposits due to the
82
difficulty and high cost of in-situ sampling. Winters et al. (2007) conducted several triaxial
83
compression tests for the GHBS extracted from the Mallik 2 L-38, Mackenzie Delta, Canada, and
84
indicated that the GHBS was higher than the natural hydrate-free sample in shear strength. Masui et
85
al. (2008) and Yoneda et al. (2015a, 2015b) investigated the mechanical behavior of GHBS
86
collected from the Eastern Nankai Trough, Japan. More recently, Yoneda et al. (2019) conducted
87
isotropic consolidation tests and triaxial compression tests on the samples of GHBS taken from the
88
Krishna-Godavari Basin, offshore India.
89
At the same time, the mechanical behaviors of GHBS have been investigated theoretically to
90
evaluate the wellbore stability or the submarine slope failure during hydrate dissociation (Sasaki,et
91
al., 2018; Yan et al., 2018; Malagar et al., 2019; Song et al., 2019; Yu et al., 2019). Up to now,
92
different types of constitutive models for GHBS have been presented, such as the nonlinear models
93
by Dai et al. (2012), Miyazaki et al. (2012), Pinkert and Grozic (2014), Pinkert et al. (2015), and
94
Yan et al. (2017), the elasto-plastic models by Klar et al. (2010), Sultan and Garziglia (2011),
95
Uchida et al. (2012, 2016), Sun et al. (2015), Lin et al. (2015, 2017), Shen et al. (2016), Gai and
96
Sánchez (2017), Sánchez et al. (2017), Yan and Wei (2017), and Zhou et al. (2018), the
97
elasto-visco-plastic model by Kimoto et al. (2010), and the hypoplastic model by Zhang (2017). For
98
example, Miyazaki et al. (2012) developed a nonlinear Duncan-Chang type of model for GHBS.
99
Although this model reasonably predicted the mechanical behavior of GHBS in the stiffness and
100
strength, it did not properly describe the volume variation of GHBS during loading and unloading.
101
Uchida et al. (2012, 2016) presented a Cam-clay type of model for GHBS, and verified its modeling
102
ability by comparing with the published test data. Sun et al. (2015) proposed a
103
thermodynamics-based model of GHBS with the concept of critical state, and confirmed that the
104
model well captured the stress softening and dilation behaviors of GHBS under different
105
morphologies and saturations of hydrate. Lin et al. (2015, 2017) proposed a subloading type of 4
106
model within the framework of spatial mobilized plane. The model capability was verified by
107
comparing with the triaxial compression test data of the synthetic and in-situ core specimens under
108
different hydrate saturations and initial effective confining stresses. Shen et al. (2016) developed a
109
state-dependent constitutive equation for GHBS. The simulation result showed that the model
110
reasonably predicted the constitutive response of GHBS at different hydrate saturations, initial
111
effective confining stresses, and void ratios. Gai and Sánchez (2017) and Sánchez et al. (2017)
112
presented the elasto-plastic model of GHBS in the hierarchical single-surface framework, and
113
concluded that the model had a good performance by comparing with the experimental data for the
114
synthetic and natural core samples with different hydrate morphologies, hydrate saturations, and
115
initial effective confining stresses. More recently, Sánchez et al. (2018) applied their constitutive
116
model and THMC to simulate laboratory experiments and field production tests.
117
The above works on GHBS are similar to the experimental, theoretical, and numerical
118
investigations on the fouling of nanoparticles (Nikkhah et al., 2015a, 2015b; Kamalgharibi et al.,
119
2016; Hemmat Esfe et al., 2018) and the CO2 sequestration (Hajirezaie et al., 2015; Rostamian et al.,
120
2016, 2019; Golkhou et al., 2019). They are closely related in thermodynamic properties. In
121
particular, when CO2 is injected into GHBS, the kinetics of CO2-CH4 hydrate in the reservoir
122
strongly depends on the physical mechanism of gas hydrate at the molecular level (Boswell et al.,
123
2017).
124
The existing constitutive models of GHBS are basically based on the coaxial assumption between
125
the directions of principal stress and plastic strain increment, and are mostly state-independent. To
126
overcome these shortcomings, this study presents a state-dependent multishear bounding surface
127
constitutive relation for GHBS. On the basis of plasticity slip theory, the formulation is obtained by
128
splitting a macro constitutive response into a macro volume response and a series of micro shear
129
responses in spatial distributions related to microshear structures. Each microshear structure
130
describes micro shear and dilatancy responses in three orthogonal orientations. A relation of the
131
micro stress–strain and a relation of the micro stress–dilatancy are established for each orientation
132
of the microshear structure. The hydrate-related state parameter and the evolution law of hydrate
133
bonding and debonding are introduced into the model. Moreover, in order to consider the impact of
134
hydrate dissociation on the deformation of GHBS, a simple collapse calculation method is proposed.
135
The validity of the constitutive relation is verified by comparing the calculated values with the test
136
results of GHBS under pore pressures of 8–15 MPa, temperatures of 1–10 °C, initial void ratios of
137
0.5–1.0, hydrate saturations of 0–0.8, and initial effective confining stresses of 1–5 MPa. 5
138
2. Model description
139
2.1. Formulation of multishear model
140
GHBS are assemblages of soil particles with large amounts of different sizes and shapes as well
141
as hydrates hosted in sediments. Soil particles interact with each other and endure external loads.
142
The role of hydrate differs from its accumulation habit. Fig. 1 gives the schematic diagram of the
143
multishear model through a unit sphere that represents the elementary volume of GHBS. There are a
144
number of randomly oriented microcontacts in space (Fig. 1(a)), which contribute to the macro
145
response of GHBS under loading. Following the concept of plasticity slip (Taylor, 1938), a macro
146
constitutive response of GHBS is split into a macro volume response and a series of micro shear
147
responses in spatial distributions related to microcontacts in different directions, in which the micro
148
shear contributes mainly to the macro deviatoric response (Fang et al., 2017; Fang et al., 2019). In
149
principle, this splitting approach somewhat resembles the volumetric-deviatoric decomposition used
150
in the microplane model (Caner and Bažant, 2013). Note that the micro shear here is not defined at
151
the real micro scale, but at the virtual equivalent scale smaller than the macro scale. As shown in
152
Fig. 1(b), for any microcontact colored in orange, a unit normal vector n on the sphere can be
153
found, which is consistent with the normal direction of the microcontact. For each microcontact, a
154
virtual microshear structure is assumed, which has the same direction as the microcontact (Fig. 1(c)).
155
The component in n is expressed with ni , in which the subscript i represents the global
156
coordinate xi . l and m refer to the orthogonal unit vectors, defining two tangential orientations
157
of the microshear structure. The components in l and m are denoted by li and mi , respectively.
158
To establish the relation of the micro shear and macro deviatoric responses, a macro strain is split
159
into the volumetric and deviatoric parts, i.e., ε ij = ε v δ ij / 3 + eij , where ε ij represents the macro
160
strain tensor, ε v = ε ijδ ij refers to the macro volumetric strain, eij represents the macro deviator
161
strain tensor, and δ ij denotes the delta of Kronecker. Generally, projecting a macro deviator strain
162
tensor onto the microcontact forms two tangential strain components corresponding to the l and
163
m directions and one deviatoric normal strain component corresponding to the n direction. As the
164
projected strain components of a macro deviator strain tensor, there is no difference among them
165
except for the direction. Hence, they can all be called micro shear strains, expressed as (see
166
Appendix A.1 for details or Fang et al. (2017))
167
γ ( k ) (n)=N ij( k ) (n )ε ij ( k =1, 2, 3) 6
(1)
168
where (k) denotes the kth component in the microshear structure; 1, 2, and 3 in k mean the
169
directions of l , m , and n in the microshear structure, respectively; γ (k ) (n ) refers to the kth
170
component of a micro shear strain in the microshear structure of the n direction; and N ij(k ) (n )
171
represents the projection tensor for the kth component of a micro shear strain in the microshear
172
structure of the n direction, which is given by equation (A4) in the appendix.
173
The kth component of a micro shear stress is expressed as τ (k ) (n ) in the microshear structure of
174
the n direction and its direction is assumed to be consistent with the corresponding micro shear
175
strain. In general, it is not possible for the micro strain and stress on the microcontact to be the
176
tensor projections for the macro strain and stress, respectively. Hence, a static equilibrium between
177
the macro and micro levels is approximately enforced in terms of the virtual work principle for all
178
strains and stresses. With the hypothesis that all microshear structures are mutually independent as
179
well as the micro strain and stress distribute uniformly in each microshear structure, the increment
180
of macro effective stress is derived as (see Appendix A.1 for details or Fang et al. (2017)) M
181
3
dσ ij = dpδij + ∑∑ 2w(m) Nij(mk )dτ (mk )
(2)
m=1 k =1
182
where w ( m ) denotes the weight coefficient in the mth microshear structure; σ ij represents the
183
macro effective stress tensor; p = σ ii /3 refers to the macro mean effective stress; m = 1, 2, ..., M
184
represent a collection of integral points, each of which corresponds to the microshear structure with
185
the direction of n ( m ) ; M denotes the total integral point number, which is set to 21 for an efficient
186
integration formula with an enough accuracy proposed by Bažant and Oh (1985); and (mk)
187
represents the kth component of a variable in the mth microshear structure with the direction of
188
n ( m ) , such as N ij( mk ) = N ij( k ) (n ( m ) ) and τ ( mk ) = τ ( k ) (n ( m ) ) . Equation (2) indicates that the increment
189
of macro stress consists of a direction-independent volumetric part and a deviatoric part of the
190
summation of micro shear stresses for all microshear structures.
191
2.2. Critical state and state parameter
192
The impact of hydrate on the sediment structure is dependent on the saturation and accumulation
193
habit of hydrate. There exist three main kinds of hydrate morphologies in sediments, i.e., the
194
pore-filling hydrate, the cementing hydrate, and the load-bearing hydrate (Masui et al., 2005), as
195
shown in Fig. 2. The pore-filling hydrate can significantly affect the permeability of sediments,
196
which also affects the mechanical behavior of sediments when the saturation of hydrate is high. 7
197
Particularly, the hydrate becomes the load-bearing hydrate if its saturation rises to the order of 0.3–
198
0.4. The cementing hydrate acts as the bonding material at inter-granular contacts and can
199
significantly increase the strength and stiffness of sediments by bonding near particles. The
200
load-bearing hydrate is a part of the sediment skeleton and can strengthen the skeleton structure of
201
sediments.
202
Hyodo et al. (2013a, 2013b) performed a set of triaxial compression tests for samples of synthetic
203
GHBS to study impacts of pore pressure, temperature, initial void ratio, hydrate saturation, and
204
initial effective confining stress on the deformation and strength behaviors of GHBS. The samples
205
with host sands were made through the partial water saturation approach, which belonged to the
206
cementing hydrate. Fig. 3 illustrates the relationship of the critical void ratio against mean effective
207
stress for synthetic GHBS, indicating that the location of critical state line is significantly affected
208
by the saturation of hydrate. Although the experimental data are somewhat scatted, the linear
209
relationship in the e − ln p plane can be approximately described for samples with any hydrate
210
saturation, where e refers to the void ratio. With the increase of the saturation of hydrate, the
211
critical state line moves upward. The relationship of the critical void ratio against mean effective
212
stress is described by
213
ec = eΓ − λc ln( p / pa )
214
where ec represents the void ratio at a critical state; λc refers to the material constant; pa
215
denotes the atmospheric pressure; and eΓ represents the critical void ratio parameter which
216
depends on the saturation of hydrate, defined as (Shen et al., 2016)
(3)
217
eΓ = eΓ 0 + ac Shnc
218
where eΓ 0 refers to the material parameter for host materials; ac and nc denote material
219
parameters; and S h represents the hydrate saturation. Fig. 4 illustrates the variation of eΓ with the
220
hydrate saturation obtained from Fig. 3, in which its fitting curve is also given.
(4)
221
Fig. 5 shows the relation of the critical stress ratio versus hydrate saturation for synthetic GHBS,
222
where pp and T represent the pore pressure and temperature, respectively. It can be seen from
223
this figure that although the test data are relatively scatted, the trend that the critical stress ratio is
224
directly proportional to the hydrate saturation is very obvious. Therefore, for simplicity, it is
225
assumed that the critical stress ratio increases linearly with the hydrate saturation, given by
226
M c = M c0 (1 + aM Sh )
8
(5)
227
where M c0 and M c represent the critical stress ratios of hydrate-free sediments and GHBS,
228
respectively; and aM denotes the material parameter, reflecting the influence of hydrate on the
229
critical stress ratio.
230
Some studies indicated that the stress and density affect the mechanical behavior of granular soils
231
(Been and Jefferies, 1985; Sun et al., 2008; Zhao and Guo, 2013). Similarly, to consider their
232
impacts in the model, the state parameter ψ for sands (Been and Jefferies, 1985) is applied. This
233
state parameter represents the interval of the critical state and current void ratios at the same mean
234
effective stress, expressed as
ψ = e − ec
235 236
(6)
2.3. Evolution law of hydrate bonding and debonding
237
Similar to cemented soils (Nguyen et al., 2014), the inter-particle bonding of hydrate can be
238
defined by a tensile strength of GHBS. This strength is called the bonding strength, which can be
239
calculated in terms of the experimental data. Fig. 6 illustrates the impact of saturation in hydrate on
240
the peak strength line of GHBS in the p − q plane, where q and qf are the deviator stress and
241
peak strength, respectively. From this figure, it can be observed that the mean effective stress at the
242
deviator stress of zero represents the initial value of the bonding strength, which is dependent on the
243
saturation of hydrate. Therefore, the relation of the initial value of the bonding strength and the
244
saturation of hydrate can be obtained from the test results as plotted in Fig. 6, which is
245
approximately expressed as
pt0 = at Sh
246
(7)
247
where pt0 represents the initial bonding strength; and a t denotes the material parameter. Fig. 7
248
shows the fitting result in the relation of the initial bonding strength versus hydrate saturation.
249
To describe the cementation and deterioration of hydrate during deformation, a simple evolution
250
law for the bonding strength is used, expressed as (Uchida et al., 2012, 2016; Sun et al., 2015; Lin
251
et al., 2015, 2017; Shen et al., 2016; Gai and Sánchez, 2017; Sánchez et al., 2017; Yan and Wei,
252
2017)
dpt = −kt pt dε qp
253
(8)
254
where p t refers to the current bonding strength; k t represents the model parameter reflecting the
255
p debonding rate; and ε q denotes the equivalent plastic shear strain. The integral form of equation
256
(8) is written as 9
pt = pt0 exp(−ktε qp )
257 258
(9)
2.4. Characteristic surfaces and lines
259
To consider the effect of bonding strength in the constitutive relation of GHBS, the macro
260
modified effective stress tensor and mean stress are expressed as σˆ ij = σ ij + ptδ ij and
261
pˆ = p + pt , respectively. Following the bounding surface theory, two types of characteristic
262
surfaces are described, one is conical surface and the other is flat cap. The conical surface depends
263
on the elasto-plastic loading owing to a change in the modified stress ratio defined as
264
rˆij = (σˆ ij − pˆ δij ) / pˆ , and the cap depends on the elasto-plastic loading owing to a change in the
265
modified mean effective stress at a unchanged modified stress ratio. Fig. 8 illustrates four
266
cone-shaped characteristic surfaces whose apexes are located at the origin. The bounding surface
267
represents the condition where the modified stress ratio rises to its peak. The surface of critical state
268
denotes the condition where any change in stress and void does not occur while the plastic shear
269
deformation unceasingly increases. The dilatancy surface refers to the condition where the
270
transformation from contraction to dilation takes place. The surface of maximum stress ratio
271
represents the condition where its maximum value occurs in loading histories. The shapes of four
272
conical surfaces are similar, given by the invariant of the modified stress ratio, the yielding shape
273
function, and their stress ratios at triaxial compression as (Fang et al. 2017)
Fi =Rˆ − Mi g (θˆ) = 0 (i=m, b, c, d)
274
(10)
275
where Fm , Fb , Fc , and Fd refer to the macro maximum stress ratio, bounding, critical state,
276
dilatancy surfaces, respectively; Rˆ represents the macro modified stress ratio, given by
277
Rˆ =
278
and dilatancy stress ratios under the condition of triaxial compression, respectively; and g (θˆ )
279
denotes the yielding shape function; θˆ refers to the angle of Lode. According to equation (10),
280
M m is equal to the maximum value of Rˆ / g (θˆ ) during loading. M d and M b are expressed as
281
(Li, 2002)
282 283
3 / 2rˆij rˆij ; M m , M b , M c , and M d mean the macro maximum, bounding, critical state,
M d = M c exp(ndψ ) , M b = M c exp(−nbψ ) where nd and nb are model constants.
10
(11)
284
To involve the influence of the intermediate principal stress in the model, the yielding criterion of
285
Matsuoka and Nakai (1974) or the yielding criterion of Lade and Duncan (1975) is applied. For the
286
yielding criterion of Matsuoka and Nakai, g (θˆ ) is written as (Matsuoka et al., 1999)
(
1 g (θˆ) = Rˆ 3 6
287
( Iˆ Iˆ − Iˆ ) ( Iˆ Iˆ − 9Iˆ ) − 1) 1 2
3
1 2
(12)
3
288
3 / 6 refer to where Iˆ1 = σˆ ij , Iˆ2 = (σˆ ii2 − σˆ rsσˆ sr ) / 2 , and Iˆ3 = σˆ ijσˆ jkσˆ ki / 3 − σˆ rsσˆ srσˆ mm / 2 + σˆ nn
289
invariants of the modified stress tensor. For the yielding criterion of Lade and Duncan, g (θˆ ) is
290
written as (Yao and Sun, 2000)
291 292
Rˆ 1 ˆ 3 g (θˆ) = I 3 /pˆ 1 − 3 pˆ 2
1 −1 ˆ 3 cos 3 cos ( − I 3 /pˆ )
−1
−1
(13)
The bounding surface of flat cap is expressed as (Fang et al. 2017)
Fp =pˆ − pˆ m = 0
293
(14)
294
where Fp refers to the bounding surface of flat cap; and pˆ m denotes the maximum modified
295
mean effective stress in loading histories.
296 297
Similarly, four characteristic lines for the kth component in the mth microshear structure, as illustrated in Fig. 9, are expressed as (Fang et al. 2017)
Fi ( mk ) =τ ( mk ) − ri( mk ) pˆ = 0 (i=m, b, c, d)
298
(15)
299
where Fm(mk ) , Fb(mk ) , Fc(mk ) , and Fd(mk ) denote the micro maximum stress ratio, bounding, critical
300
state, and dilatancy lines for the kth component in the mth microshear structure, respectively;
301
r (mk ) = τ (mk ) /pˆ represents the kth component of the micro stress ratio in the mth microshear
302
structure; and rm( mk ) , rb( mk ) , rc( mk ) , and rd( mk ) mean the micro maximum, bounding, critical state,
303
and dilatancy stress ratios for the kth component in the mth microshear structure, respectively.
304
These micro stress ratios correspond to different micro stress states. Their physical meanings are
305
similar to those corresponding to the macro modified stress ratios. According to equation (15),
306
rm( mk ) is set to the maximum value of r (mk ) in loading histories. ri ( mk ) (i=b, c, d) is given by
307
equation (A28) in the appendix.
11
308
2.5. Macro volumetric deformation
309
The macro volumetric strain without that due to dilatancy usually consists of elastic and plastic
310
parts, which depends on the modified mean effective stress. The incremental relationship of them is
311
given by (see Appendix A.2 for details or Fang et al. (2017))
dε v − dε vd =
312
dpˆ 1 1 dpˆ + h( pˆ − pˆ m ) dpˆ Ke Kp dpˆ
(16)
313
where K e denotes the elastic bulk modulus; K p represents the plastic bulk modulus; ε vd refers
314
to the macro volumetric strain due to dilatancy; h represents the Heaviside step function; and
315
refers to the Macaulay brackets.
316
The elastic bulk modulus is given by Ke =
317
1+ e
κ
ˆ a ) 0.5 (1 + aG S h )( pp
(17)
318
where aG refers to the material parameter; the term aG S h represents the influence of hydrate on
319
the modulus of GHBS; and κ denotes the material parameter related to the unloading in isotropic
320
compression, expressed as κ = κ 0 [(1 + e) / (2.97 − e)]2 where κ 0 represents the model constant
321
for host materials.
322 323
The plastic bulk modulus is defined as (Wang et al., 1990)
Kp =
Mc g(θˆ) 1+ e ˆ a )0.5 (1+ aGSh )( pp λ −κ Mc g(θˆ) − Rˆ
(18)
324
where λ represents the material parameter related to the isotropic compression, expressed as
325
λ = λ0 [(1 + e ) / (2.97 − e )]2 where λ0 denotes the model constant for host materials.
326 327
The macro volumetric strain due to dilatancy is the sum of the micro volumetric strains due to dilatancy for all microshear structures, written as (Fang et al., 2017) M
328
3
( mk ) dε vd = ∑ ∑ 2 w( m ) dε vd
(19)
m =1 k =1
329
where ε vd( mk ) denotes the kth component of the micro volumetric strain due to dilatancy in the mth
330
microshear structure. Following the cyclic torsional shear test result of sands, the micro stress–
331
dilatancy relation is given by (Fang et al., 2017)
332
dε vd( mk ) = d1 (±rd( mk ) − r ( mk ) )dγ p(mk )
12
(20)
333
where d1 represents the micro dilatancy parameter, which is assumed to be the same for all
334
microshear structures; ± is set as positive for dγ p(mk ) > 0 and negative for d γ p(mk ) < 0 ; and
335
γ p(mk ) denotes the kth component of the micro plastic shear strain in the mth microshear structure,
336
given by equation (A15) in the appendix.
337
2.6. Micro shear deformation
338
The micro shear strain usually consists of elastic and plastic parts. The former depends on the
339
micro shear stress, and the latter depends on the micro shear stress ratio and modified mean
340
effective stress. According to the plasticity theory of bounding surface (Wang et al., 1990), their
341
incremental correlation is expressed as (see Appendix A.3 for details or Fang et al. (2017))
342
dγ ( mk ) =
1
dτ ( mk ) + ( mk )
Ge
1 Gp( mk )
pˆ dr ( mk ) +
1
h( pˆ − pˆ m ) ( mk )
Hp
dpˆ ( mk ) r dpˆ dpˆ
(21)
343
where γ ( mk ) denotes the kth component of the micro shear strain in the mth microshear structure;
344
Gp(mk ) and H p(mk ) refer to the micro plastic shear moduli related to d r ( mk ) and dpˆ for the kth
345
component in the mth microshear structure, respectively; and Ge(mk ) represents the micro elastic
346
modulus for the kth component in the mth microshear structure, which depends on the macro elastic
347
shear modulus and is given by equation (A29) in the appendix. The macro elastic shear modulus
348
Ge is expressed as (Richart et al., 1970)
349
Ge = G0 (1 + aG S h )
(2.97 − e) 2 ˆ a ) 0.5 ( pp 1+ e
(22)
350
where G0 refers to the elastic shear modulus parameter for host materials; and G0 (1 + aG S h ) = Gh0
351
denotes the elastic shear modulus parameter for GHBS. Fig. 10 shows the Gh0 − S h relationship
352
from the experimental data (Hyodo et al., 2013a, 2013b).
353 354
The micro plastic shear modulus Gp(mk ) related to d r ( mk ) is defined as (Fang et al., 2017)
r ( mk ) ρ ( mk ) Gp( mk ) = h1Ge(mk ) b( mk ) 1( mk ) −1 rm ρ1
(23)
355
where h1 denotes the model parameter; ρ1( mk ) = r (mk ) − rr(mk ) represents the interval between the
356
current and back values of the kth micro stress ratio in the mth microshear structure; rr(mk ) denotes
357
the back value of the kth micro back stress ratio in the mth microshear structure, defined as the
358
origin during a virgin loading or the latest inflection location during a reverse loading; and
13
359
ρ1( mk ) = ± rm(mk ) − rr(mk ) represents the interval between the maximum and back values of the kth
360
micro stress ratio in the mth microshear structure, where ± is set as positive for d r ( mk ) > 0 and
361
negative for d r ( mk ) < 0 .
362
The micro plastic shear modulus H p(mk ) related to dpˆ is defined as (Fang et al., 2017)
H
363
( mk ) p
( mk ) ( mk ) c e ( mk )
= h2G
r r
ρ2 ρ2
(24)
364
where h2 denotes the model parameter; ρ 2 = pˆ − pˆ r refers to the interval between the current
365
and back values of the modified mean effective stress; pˆ r represents the back value of the
366
modified mean effective stress, defined as the origin during a virgin loading or the latest inflection
367
location during a reverse loading; and ρ 2 = pˆ m − pˆ r denotes the interval between the maximum
368
and back values of the modified mean effective stress where ρ 2 = pˆ r for d pˆ < 0 .
369
2.7. Collapse deformation due to hydrate dissociation
370
Klar et al. (2010) and Sánchez et al. (2017) proposed different methods to predict the collapse
371
deformation caused by hydrate dissociation. The Klar et al.’s method follows the concept of stress
372
relaxation. Due to the dissociation of hydrate, the effective stress in the hydrate skeleton gradually
373
releases, resulting in the decrease of the total effective stress in the skeleton of GHBS under
374
constant strain. In order to restore the released stress, an additional deformation is required. Because
375
the initial elastic modulus is applied in the Klar et al.’s method, the collapse deformation may be
376
underestimated. To overcome this drawback, an improved method is proposed. During hydrate
377
dissociation without external load, it can be assumed that the sediment is self-balanced through its
378
skeleton structure, and the effective stress remains unchanged from the time of t to the time of t+dt,
379
i.e.,
380
σ ijt = ( K eqε vδ ij + 2Geq eij )t = σ ijt +dt = ( K eqε vδ ij + 2Geq eij )t +dt
381
where the superscripts t and t+dt of a variable represent its values at the time of t and the time of
382
t+dt, respectively; and K eq and Geq denote the equivalent bulk and shear moduli, respectively.
383
Since σ ijt , ε vt , and eijt are known at the time of t, K eqt and Geqt can be approximately obtained
384
by using the linear elastic relationship between the stress and strain. By assuming that the
385
increments of the equivalent bulk and shear moduli are proportional to the variation of hydrate
386
saturation during the time increment of dt, K eqt + dt and Geqt + dt are approximately given by
14
(25)
387
K eqt +dt = K eqt +dK eq = K eqt (1+dSh )
(26)
388
Geqt +dt = Geqt +dGeq = Geqt (1+dSh )
(27)
389 390
From equations (25)–(27), the increment of the collapse strain produced by hydrate dissociation is given by
dε ijdis = ε ijt +dt − ε ijt = −[ pt / (3K eqt )δ ij + (σ ijt − p tδ ij ) / (2Geqt )]dSh / (1 + dSh )
391 392
where ε ijdis refers to the collapse strain caused by hydrate dissociation.
393
2.8. Constitutive formulation
394 395
(28)
The incremental relation of the stress and strain, including the effect of hydrate dissociation, is described as
dσ ij = Dijrs (dε rs − dε rsdis ) − dptδij
396
(29)
397
where Dijrs denotes the macro elasto-plastic stiffness tensor related to macro and micro quantities
398
(see Appendix A.4 for details), which enables the model to involve the influence of micro quantities
399
as well as their evolutions at different directions of microshear structures on the constitutive
400
response of GHBS.
401
3. Parameter calibration
402
The model has 17 parameters, 11 of which are associated with the mechanical behavior of host
403
materials, and the remaining 6 are associated with hydrate. The 11 model parameters for host
404
materials include two elastic modulus parameters ( G0 and κ 0 ), three critical state parameters
405
( M c0 , eΓ 0 , and λ c ), four plastic modulus parameters ( h1 , h2 , nb , and λ0 ), and two dilatancy
406
parameters ( d1 and nd ). These parameters for host materials are determined through the following
407
approaches: (1) G 0 is obtained by the relation between the deviator stress and axial strain for
408
triaxial tests. λ0 and κ 0 are determined by the e − p relations for consolidation tests. κ 0 is
409
also set by κ 0 = 3(1 − 2ν ) / [2(1 + ν ) G 0 ] , where ν refers to the ratio of Poisson; (2) M c0 , eΓ 0 ,
410
and
411
nd = ln( M d / M c ) / ψ d , where ψ d and M d take the values of ψ and stress ratio when the phase
412
transformation takes place, respectively. nb is determined by nb = ln( M c / M b ) / ψ b , where ψ b
413
and M b take the values of ψ and stress ratio when the peak of stress ratio occurs, respectively;
λc
are
directly
determined
by
triaxial
15
test
results;
(3)
nd
is
obtained
by
414
(4) d1 is obtained according to the volumetric and deviatoric strain data of triaxial tests; (5) h1
415
and h2 are set by the trial-and-error procedure to match optimally the model predictions with the
416
test results.
417
The additional 6 parameters for GHBS consist of one elastic modulus parameter ( aG ), three
418
critical state parameters ( ac , nc , and aM ), one bonding strength parameter ( a t ), and one strength
419
deterioration parameter ( k t ). aG is determined in terms of a variation in elastic shear modulus
420
with the saturation of hydrate, as shown in Fig. 10. ac and nc are obtained by the regression of
421
the critical void ratio with a variation in hydrate saturation, as illustrated in Fig. 4. aM is set by
422
defining that the critical stress ratio is dependent on the saturation of hydrate, as shown in Fig. 5.
423
a t for the initial bonding strength is given through the experimental peak strength data for GHBS,
424
as illustrated in Figs. 6 and 7. k t is calibrated by the trial-and-error procedure to simulate the
425
mechanical response of GHBS in tests.
426
4. Model performance
427
The model performance is evaluated by simulating the available experimental results published in
428
literatures. Five sets of triaxial compression tests for GHBS are used. Two sets are the tests for the
429
samples of synthetic GHBS (Masui et al., 2005; Hyodo et al., 2013a), and the other three sets are
430
the tests for the samples of natural GHBS taken from the Eastern Nankai Trough, Japan (Masui et
431
al., 2008; Yoneda et al., 2015b) and the samples of natural GHBS extracted from the Krishna–
432
Godavari Basin, offshore India (Yoneda et al., 2019). Moreover, the hypothetical torsional shear
433
tests during principal stress direction rotations are simulated to investigate the model ability for
434
predicting the non-coaxiality of GHBS. The simulation results are described as follows.
435
4.1. Triaxial compression tests for synthetic GHBS by Masui et al. (2005)
436
Masui et al. (2005) performed drained triaxial compression tests for synthetic cementing and
437
pore-filling types of GHBS, in which the pore-filling type of hydrate formed through the ice-seed
438
approach and the cementing type of hydrate formed through the partial water saturation approach.
439
The host materials were made of Toyoura sands. Tests were carried out under the pore pressure of 8
440
MPa, temperature of 5 °C, initial effective confining stress of 1 MPa. For the pore-filling type of
441
GHBS, the initial void ratio of the samples was 0.736 and the saturation of hydrate was in the extent
442
of 0.264–0.678. For the cementing type of GHBS, the initial void ratio of the samples was 0.605
16
443
and the saturation of hydrate was in the extent of 0.0–0.551. In the tests, the deviator stress response
444
data of all samples were successfully measured, but only the volume change data of the sample with
445
a hydrate saturation of 0.409 for the pore-filling hydrate and the samples with hydrate saturations of
446
0.0 and 0.407 for the cementing hydrate were obtained. The model parameters for host sands and
447
GHBS used in simulations are given in Table 1. Note that ac and nc are set in two regions
448
according to the properties of pore-filling hydrate.
449
Fig. 11 compares the predicted values with the experimental results for the pore-filling type of
450
GHBS, where e0 and p c0 refer to the initial values of the void ratio and effective confining stress,
451
respectively. As can be observed from this figure, most of the predicted deviator stresses and
452
volumetric strains agree well with those of the test data. Only the predicted deviator stress responses
453
for the hydrate saturation of 0.678 are larger than the experimental values when the axial strain
454
exceeds 8%. The experimental deviator stress for the hydrate saturation of 0.678 is quite different
455
from that in other cases. After the axial strain of 8%, its value decreases significantly. The reason
456
for this phenomenon was not explained in the paper of Masui et al. (2005), which may be that the
457
hydrate partially decomposed during the experiment of this case. Fig. 12 compares the calculated
458
values with the test data for the cementing type of GHBS. In the relation of the deviator stress and
459
axial strain, the calculated values are consistent with the test data before the axial strain of 9%, but
460
the calculated values deviate from the test data after the axial strain exceeds 9–12%. In the relation
461
of the volumetric and axial strains, there is a certain difference between the predicted and
462
experimental values, especially for the case with a hydrate saturation of 0.407. Moreover, it can be
463
observed from Figs. 11 and 12 that the peak strength, initial stiffness, and volumetric dilation for
464
the cementing samples are larger than those for the pore-filling samples although their saturations of
465
hydrate are almost the same. The calculated results show that the constitutive relation captures the
466
different behaviors of GHBS for two accumulated types of hydrates, particularly in the deviator
467
stress response, as observed in the tests.
468
4.2. Triaxial compression tests for synthetic GHBS by Hyodo et al. (2013a)
469
Hyodo et al. (2013a) conducted drained triaxial compression tests for synthetic GHBS. Toyoura
470
sands were employed as host materials. The initial void ratio of the samples was in the range of
471
0.629–0.842 and the saturation of hydrate was in the extent of 0.0–0.695. Tests were performed
472
under different pore pressures (5, 10 and 15 MPa), temperatures (1, 5, and 10 °C), and initial
473
effective confining stresses (1, 3, and 5 MPa). The model parameters for host sands and GHBS used
17
474
in simulations are given in Table 1. Fig. 13 compares the predicted values with the test results of
475
GHBS containing different hydrate saturations under the pore pressure of 10 MPa, temperature of
476
5 °C, and initial effective confining stress of 5 MPa. From this figure, it can be observed that the
477
model can generally well predict the variations of the deviator stress and volumetric strain with the
478
axial strain. The simulated and experimental results show that, in addition to the great dilative
479
behavior of the sample containing a high hydrate saturation of 0.531, the samples exhibit the
480
dominant compression in volume and the behavior of strain hardening. Furthermore, the strength
481
and initial stiffness are proportional to the saturation of hydrate.
482
Fig. 14 shows the comparison between the simulated results and test data for GHBS with almost
483
the same hydrate saturation and initial void ratio under different pore pressures (10 and 15 MPa),
484
temperatures (1, 5, and 10 °C), and initial effective confining stresses (1, 3, and 5 MPa). In the
485
relation of the deviator stress and axial strain, the calculated values well consist with the test data
486
for two cases under initial effective confining stresses of 1 and 5 MPa, while the calculated values
487
for the initial effective confining stress of 3 MPa are agreement with the test data under the
488
temperature of 1 °C. In the relation of the volumetric and axial strains, the predicted values are very
489
consistent with the test data under the initial effective confining stress of 1 MPa, while the
490
calculated values are quite different from the test data under initial effective confining stresses of 3
491
and 5 MPa. In addition, when the initial effective confining stress is 3 MPa, the experimental
492
deviator stress responses at 5 and 10 °C are very close, while the experimental deviator stress
493
response at 1 °C is significantly greater than that at other two temperatures; the experimental
494
volume responses at these three temperatures do not change much with temperature. These test
495
results show that the decrease of temperature leads to the nonlinear increase of strength, while the
496
volume deformation is less affected by temperature. In contrast, the simulation results with different
497
temperatures are very close under the same initial effective confining stress because the impact of
498
temperature is not directly considered in this model.
499
4.3. Triaxial compression tests for natural GHBS by Masui et al. (2008)
500
Masui et al. (2008) carried out drained triaxial compression tests for natural GHBS taken from
501
the Eastern Nankai Trough, Japan. The natural GHBS were mainly composed of silty sands. Four
502
samples of the natural GHBS were sheared using the triaxial compression apparatus. The initial
503
void ratios of four samples were 0.792, 0.961, 0.972, and 1.033, respectively. The corresponding
504
saturations of hydrate were 0.376, 0.225, 0.094, and 0.077, respectively. Tests were performed
18
505
under the pore pressure of 9 MPa, temperature of 5 °C, and initial effective confining stress of 1
506
MPa. The model parameters for host sands and GHBS used in simulations are given in Table 1. Fig.
507
15 compares the calculated values with the experimental data at different saturations of hydrate. As
508
can be observed from this figure, the simulated values generally consist with the experimental
509
results in the relation of the deviator stress and axial strain and in the relation of the volumetric and
510
axial strains except for the volumetric strain response of the sample containing a hydrate saturation
511
of 0.094. The experimental and predicted results indicate that the peak value in strength and the
512
initial value in stiffness for natural GHBS are directly proportional to the saturation of hydrate. The
513
higher the saturation of hydrate, the larger the volumetric change of the sample. The sample shows
514
a significant strength reduction and volume change when the saturation of hydrate is equal to 0.376.
515
In contrast, the sample shows a low level of softening and dilative response when the saturation of
516
hydrate is equal to 0.077. For the sample with a constant saturation of hydrate, softening takes place
517
after the peak in strength and its reduction is proportional to the saturation of hydrate.
518
4.4. Triaxial compression tests for natural GHBS by Yoneda et al. (2015b)
519
Yoneda et al. (2015b) performed drained triaxial compression tests for natural GHBS taken from
520
the Eastern Nankai Trough, Japan. The natural GHBS were mainly composed of silty sands. Three
521
samples of the natural GHBS and their reconstituted samples were sheared by using the triaxial
522
compression apparatus. The saturations of hydrate for three samples were 0.38, 0.7, and 0.79,
523
respectively. Two samples with hydrate saturations of 0.38 and 0.79 were successfully tested, but
524
the volume change data for the sample with a hydrate saturation of 0.7 were not obtained due to the
525
liquid leakage in the cell. Therefore, the test for the sample with a hydrate saturation of 0.7 is not
526
simulated because its deviator stress response was also influenced by the liquid leakage. The initial
527
void ratios of all the samples were very similar in the range of 0.645–0.669. Tests were conducted
528
under the pore pressure of about 13 MPa, temperature of 10 °C, and initial effective confining stress
529
of 1.5 or 1.6 MPa. The model parameters for natural GHBS and their reconstituted samples used in
530
simulations are given in Table 1. Fig. 16 compares the calculated values with the experimental data
531
on the responses in deviator stress and volumetric strain for natural GHBS and their reconstituted
532
samples. From this figure, it can be seen that the model performance is generally satisfactory with
533
good predictions in terms of the stiffness and strength behavior as well as in terms of the volume
534
response except for the volumetric strain of the sample containing a hydrate saturation of 0.79. The
535
initial stiffness, peak strength, and residual strength significantly depend on the saturation of
19
536
hydrate. The reconstituted sediments retain the same strength as host silty sands, which is basically
537
consistent with the in-situ strength of hydrate-free sediments. The volumetric strain responses also
538
strongly depend on the saturation of hydrate. The sample with a hydrate saturation of 0.79 exhibits a
539
significant dilative behavior.
540
4.5. Isotropic consolidation tests and triaxial compression tests for natural GHBS by Yoneda et al.
541
(2019a, 2019b)
542
Yoneda et al. (2019a, 2019b) conducted isotropic consolidation tests and drained triaxial
543
compression tests on the natural GHBS extracted from the Krishna–Godavari Basin, offshore India.
544
The natural GHBS were sediments rich in sands. Two types of isotropic consolidation tests for six
545
samples of the natural GHBS with hydrate saturations of 0.321–0.716 were performed by use of a
546
triaxial apparatus. The initial effective confining stress in situ was about 1 MPa. The results of the
547
first type of test are shown in Fig. 17(a), where samples were first loaded from 0.2 to 1 MPa, then
548
hydrates were dissociated, and finally samples were loaded to 12 MPa after hydrate dissociation.
549
The results of the second type of test are shown in Fig. 17(b), where samples were first loaded from
550
0.2 to 5 MPa, then hydrates were dissociated, and finally samples were loaded to 12 MPa after
551
hydrate dissociation. During the tests, unloading and reloading were carried out at different
552
effective confining stresses.
553
The drained triaxial compression tests for five samples of the natural GHBS were performed. The
554
hydrate saturations of five samples were 0.487, 0.552, 0.634, 0.734, and 0.822, respectively. The
555
corresponding initial void ratios were 0.712, 0.634, 0.831, 0.658, and 0.656, respectively. Tests
556
were carried out under the pore pressure of about 10 MPa, temperature of 4 °C, and initial effective
557
confining stress of 1 or 1.1 MPa. The experimental data of the deviator stress response and volume
558
change were successfully measured except for the volume change data of the sample with a hydrate
559
saturation of 0.634. Fig. 18 shows the experimental results of the deviator stress response and
560
volume variation with the axial strain. From this figure, it can be observed that the experimental
561
shear strength and dilation of GHBS were not always proportional to the saturation of hydrate when
562
the saturation of hydrate was relatively high. However, the previous experimental results (Miyazaki
563
et al., 2011; Hyodo et al., 2013a, 2013b) indicated that the strength and dilation of GHBS increased
564
with the enhancement of hydrate saturation. A possible reason for this phenomenon may be the
565
hydrate amount in the samples. The saturations of hydrate in these test samples are relatively high,
566
while the saturations of hydrate in the previous tests samples are usually low. The hydrate particles
20
567
enter into the voids at low saturation of hydrate, while the hydrate particles partially occupy the
568
voids and partially block the contacts of sediment particles at high saturation of hydrate. Therefore,
569
the distribution of hydrate particles in GHBS is related to the amount of hydrate, thereby changing
570
the equivalent skeleton void ratio, which in turn affect the shear strength and dilation of GHBS.
571
Although hydrates and fines are different kinds of materials, this characteristic may be somewhat
572
similar to that of the sand-fines mixtures (Thevanayagam et al., 2002).
573
To consider the above feature in simulations, the hydrate saturation is divided into two intervals,
574
and the equivalent skeleton void ratio parameters are set respectively. The model parameters of the
575
natural GHBS applied in simulations are illustrated in Table 1. Because the sediments of the sample
576
with a hydrate saturation of 0.321 is classified as silt (Yoneda et al., 2019a), κ 0 and λ 0 for this
577
sample are taken to be 0.01 and 0.03, respectively. Fig. 17 compares the model predictions with the
578
experimental data for isotropic consolidation tests. Note that the experimental void ratio at the mean
579
effective stress of 0.2 MPa for the two cases with hydrate saturations of 0.321 and 0.639 is quite
580
different from that in other cases. Therefore, for the simulations of these two cases, the void ratio
581
corresponding to the mean effective stress of 0.2 MPa is inversely determined according to the void
582
ratio at the mean effective stress of 1.0 or 1.1 MPa. It can be seen from Fig. 17 that the volume
583
response predicted by the model basically agree with that observed in the tests. Particularly, the
584
model can reasonably predict the collapse deformation caused by hydrate dissociation, indicating
585
the effectiveness of the proposed collapse calculation method. Fig. 18 shows the predicted and
586
experimental results for triaxial compression tests. From this figure, it can be seen that the
587
calculated values generally consist with the experimental data in the relation between the deviator
588
stress and axial strain and in the relation between the volumetric and axial strains.
589
4.6. Drained torsional tests during principal stress rotations
590
To investigate the ability of the model for predicting the non-coaxiality of GHBS, the drained
591
torsional shear tests under the pure rotation of principal stress directions are simulated. Since this
592
type of test for GHBS has not been reported, the test conditions for simulations are assumed
593
according to the existing triaxial compression tests of GHBS and torsional shear tests of sands. The
594
schematic diagram of the hollow cylindrical sample is showed in Fig. 19. The hydrate saturations of
595
three samples used in simulations are 0.0, 0.3, and 0.6, respectively. The initial void ratio and
596
effective confining stress are equal to 0.65 and 1 MPa, respectively. The pure stress rotation is
597
achieved by varying only the direction of the principal stress while maintaining the other conditions
21
598
unchanged. The mean effective stress, deviator stress, and intermediate principal stress ratio are
599
taken as 1 MPa, 1 MPa, and 0.5, respectively. The model parameters are set to be the same as those
600
applied in the simulation of the Yoneda et al.’s tests (Yoneda et al., 2015b), which are given in
601
Table 1. Fig. 19 shows the predicted variations of all strain components with the principal stress
602
direction during a cycle of continuous principal stress rotations. From this figure, it can be observed
603
that all strain components decrease with the increase of saturation in hydrate, indicating that the
604
strength and stiffness of sediments increase due to hydrate. Figs. 20 and 21 illustrate the predicated
605
stress paths and vectors of the total strain increment, as well as the predicated stress paths and
606
vectors of the plastic strain increment, respectively. The magnitudes of the increment vectors for
607
total and plastic strains in these two figures are define as
(d(ε a − ε θ )) 2 +(2dε aθ ) 2 / ds and
608
(d(ε ap − ε θp )) 2 +(2dε aθp ) 2 / ds , respectively, in which the stress increment ds is expressed as
609
ds = (d((σ a − σ θ ) / (σ a + σ θ ))) 2 +(d(2σ aθ / (σ a + σ θ ))) 2 . From these two figures, it can be observed
610
that with the enhancement of hydrate saturation, the increment in the plastic strain decreases and
611
inclines to the tangential direction along the circle, indicating that the non-coaxiality of GHBS
612
increases. Therefore, it can be said that the present constitutive relation is able to predict the
613
non-coaxial response of GHBS without new parameters.
614
5. Conclusions
615
This paper presented a state-dependent multishear constitutive formulation for predicting the
616
mechanical response of GHBS. The model follows the plasticity theory of bounding surface and
617
accords with the theory of critical state for granular soils. The effectiveness of the model was
618
confirmed by simulating different laboratory tests for synthetic and natural GHBS published in
619
literatures. The following conclusions are obtained:
620
(1) The model formulation was derived by splitting a macro constitutive response of GHBS into a
621
macro volume response and a series of micro shear responses in spatial distributions related to the
622
microshear structures. This splitting approach can provide a simple and effective means to predict
623
the complex three-dimensional constitutive response of GHBS by superposition of a macro
624
response in volume and individual responses of all microshear structures.
625
(2) A state parameter related to the hydrate saturation was introduced into the model, including
626
the joint impact of the effective confining stress and density on the constitutive responses of GHBS.
627
With this state parameter, the model not only accords with the theory of critical state for granular
22
628
soils, but also can capture the contractive and dilative behaviors for loose and dense GHBS in terms
629
of a single set of model parameters.
630
(3) The evolution law of hydrate bonding and debonding related to the plastic shear strain was
631
used to consider the effect of hydrate on the constitutive response of GHBS. This evolution law
632
characterizes the cementing mechanism of hydrate and the degradation of bonding during plastic
633
deformation, respectively.
634
(4) The simple method to calculate the collapse deformation produced by the dissociation of
635
hydrate was proposed. The validity of this method was verified by predicting the volume
636
compression of GHBS observed in the hydrate dissociation tests under constant stress.
637
(5) The comparisons of the experimental and simulated results for different drained laboratory
638
tests of synthetic and natural GHBS were performed. The result confirmed that the developed
639
model can not only predict the constitutive response of GHBS under isotropic consolidation and
640
triaxial compression, but also predict the non-coaxiality of GHBS during the pure rotation of
641
principal stress directions without new parameters. The model is suitable for the case of pore
642
pressures of 8–15 MPa, temperatures of 1–10 °C, initial void ratios of 0.5–1.0, hydrate saturations
643
of 0–0.8, and initial effective confining stresses of 1–5 MPa. However, since temperature is not
644
taken as an independent variable in the present model, and the influence of temperature on the
645
strength of GHBS cannot be ignored, temperature should be directly introduced into the constitutive
646
model in the future.
647
Acknowledgements
648
The authors would like to thank the reviewers for their creative suggestions and professional
649
comments. This work was financially supported by the National Natural Science Foundation of
650
China (No. 51878605 and No. 41807224) and the Zhejiang Provincial National Science Foundation
651
of China (LQ17D020001).
652
Notation
653
ac , aG , aM , at
654
d1
dilatancy parameter
655
Dijrs
stiffness tensor
656
ds
stress increment
657
e , ec , e0
current, critical, and initial void ratios
model parameters
23
658
eij
deviator strain tensor
659
eΓ , eΓ 0
critical state parameters
660
Fb , Fc , Fd , Fm , Fp
661
Fb( mk ) , Fc( mk ) , Fd( mk ) , Fm( mk ) micro characteristic lines
662
Ge , Geq
663
( mk ) G ( mk ) , Ge( mk ) , Gp
664
G0 , G h0
shear modulus parameters
665
g (θˆ )
yield shape function
666
h
Heaviside step function
667
h1 , h2
model parameters
668
H ( mk ) , H p( mk )
micro total and plastic shear moduli
669
Iˆ1 , Iˆ2 , Iˆ3
stress invariants
670
K , Ke , Keq , K p
671
K1
variable
672
kt
model parameter
673
l, m , n
unit vectors
674
M b , M c , M d , M m bounding, critical, dilatancy, and maximum stress ratios
675
M c0
critical stress ratio for host materials
676
n b , nc , nd
model parameters
677
N ij( mk )
projection tensor
678
p , pˆ
mean effective stress and modified mean effective stress
679
pa
atmospheric pressure
680
pc , pc0
current and initial effective confining stresses
681
pˆ m , pˆ r
maximum and back values of modified mean effective stress
682
pp
pore pressure
683
pt0 , pt
initial and current values of hydrate bonding strength
macro characteristic surfaces
elastic and equivalent shear moduli micro total, elastic, and plastic shear moduli
total, elastic, equivalent, and plastic bulk moduli
24
684
q , qf
deviator stress and peak strength
685
Qij , Q ( mk )
variables
686
Rˆ
macro modified stress ratio
687
rˆij
macro modified stress ratio tensor
688
r ( mk )
micro stress ratio
689
rb( mk ) , rc( mk ) , rd( mk ) , rm( mk )
690
rr( mk )
back value of micro stress ratio
691
Sh
hydrate saturation
692
t
time
693
T
temperature
694
w(m)
weight coefficient
695
α
principal angle
696
δ
variation
697
δ ij
Kronecker delta
698
ε ij , ε ijp
total and plastic strain tensors
699
ε ijdis
collapse strain tensor caused by hydrate dissociation
700
ε qp
equivalent plastic shear strain
701
ε v , ε ve , ε vp
total, elastic, and plastic volumetric strains
702
ε vd
volumetric strain caused by dilatancy
703
ε vd( mk )
micro volumetric strain caused by dilatancy
704
ε a , ε r , ε θ , ε aθ
705
ε ap , ε θp , ε apθ
706
γ (mk ) , γ e( mk ) , γ p( mk )
707
κ , κ0
bulk modulus parameters related to isotropic unloading
708
λ , λ0
bulk modulus parameters related to isotropic loading
709
λc
critical state parameter
710
ν
micro bounding, critical, dilatancy, and maximum stress ratios
total strain components
plastic strain components micro total, elastic, and plastic shear strains
Poisson ratio
25
711
θˆ
712
ρ1( mk ) , ρ1( mk ) mapping intervals for micro stress ratio
713
ρ2 , ρ2
714
σ a , σ r , σ θ , σ aθ stress components
715
σ ij , σˆij
effective stress tensor and modified effective stress tensor
716
σ1 , σ 2 , σ3
effective principal stresses
717
σˆ1 , σˆ 2 , σˆ3
modified effective principal stresses
718
τ ( mk )
micro shear stress
719
ψ
state parameter
720
Ω
hemisphere surface
721
Lode angle
mapping intervals for modified mean effective stress
Macaulay bracket
722 723
Appendix. Constitutive relation
724
A.1. Multishear model
725 726
By projecting the tensor of the macro deviator strain onto the microcontact, the strain components in the microshear structure of the n direction can be written as 1 2
1 2
727
γ (1) (n ) = (li n j + l j ni )eij = (li n j + l j ni )(ε ij − ε vδ ij / 3) = N ij(1) (n )ε ij
728
γ (2) (n ) = ( mi n j + m j ni )eij = ( mi n j + m j ni )(ε ij − ε vδ ij / 3) = N ij(2) (n )ε ij
729
1 2
1 2
γ (3) (n) = ni n j eij = ni n j (ε ij − ε vδ ij / 3) = N ij(3) (n )ε ij
(A1) (A2) (A3)
730
where γ (1) (n ) and γ (2) (n ) denote the tangential strain components corresponding to the l and
731
m directions, respectively; γ (3) (n ) represents the deviatoric normal strain corresponding to the
732
n direction; and N ij(1) (n ) , N ij(2) (n ) , and N ij(3) (n ) are expressed as
733 734 735
N ij(1) (n ) =
1 1 1 (li n j + l j ni ) , N ij(2) (n ) = (mi n j + m j ni ) , N ij(3) (n ) = ni n j − δ ij 2 2 3
(A4)
According to the virtual work principle, the stress vectors on all microshear structures and the macro stress tensor are variationally enforced, which is written as
26
736 737 738 739 740 741 742
3 4π 4π 1 dσ ij δε ij = dpδε v + 2∫ ∑ dτ ( k ) (n)δγ ( k ) (n)dΩ Ω 3 3 k =1 3
(A5)
where δ denotes a variation; and Ω represents the surface of hemisphere. By substituting δε v = δ ij δε ij and δγ ( k ) (n) = N ij( k ) (n)δε ij into equation (A5) and holding the variational equation for any variation δ ε ij , the macro effective stress increment is given as
dσ ij = dpδij +
3 M 3 1 (k ) (k ) N ( n )d τ ( n )d Ω ≈ d p δ + 2w(m) Nij(mk )dτ (mk ) ∑ ∑∑ ij ij ∫ 2π Ω k =1 m=1 k =1
(A6)
A.2. Stress–strain relation on macro volumetric deformation The increment of the volumetric strain is the sum of the elastic and plastic parts, i.e., dε v =dε ve +dε vp
743
(A7)
744
where ε ve refers to the elastic volumetric strain; and ε vp denotes the plastic volumetric strain that
745
involves the strain due to a variation in mean effective stress and the strain due to dilatancy.
746
Following the elasticity theory, dε ve is expressed as
748
751 752 753 754 755
756 757
(A8)
Following the plasticity theory of bounding surface, dε vp is given by
dεvp =
749 750
1 dpˆ Ke
dε ve =
747
dpˆ 1 h( pˆ − pˆm ) dpˆ + dεvd Kp dpˆ
(A9)
Combing equations (A7)–(A9) yields dε v =
dpˆ 1 1 dpˆ + h( pˆ − pˆ m ) dpˆ + dε vd Ke Kp dpˆ
(A10)
Equation (A10) is rewritten as
dε v =
1 dpˆ +dε vd K
(A11)
where the total bulk modulus K is expressed as 1 dpˆ 1 K = + h ( pˆ − pˆ m ) K dpˆ e Kp
−1
A.3. Stress–strain relation on micro shear deformation The increment of micro shear strain is the sum of the elastic and plastic parts, i.e.,
27
(A12)
dγ ( mk ) =dγ e( mk ) +dγ p( mk )
758
(A13)
759
where γ e( mk ) denotes the kth component of the micro elastic shear strain in the mth microshear
760
structure.
761
Following the elasticity theory, dγ e( mk ) is expressed as
dγe(mk) =
762 763
dγ p(mk ) =
1 Gp(mk )
pˆ dr (mk ) +
765
Combing equations (A13)–(A15) yields
766
dγ ( mk ) =
768
1 ( mk ) e
G
dτ ( mk ) +
1 ( mk ) p
G
1
h pˆ − pˆ m ) ( mk ) (
Hp
pˆ dr ( mk ) +
1 H
( mk ) p
dpˆ (mk ) r dpˆ dpˆ
h( pˆ − pˆ m )
dpˆ (mk ) r dpˆ dpˆ
(A15)
(A16)
A.4. Macro stress–strain relation Following the definition of r (mk ) , dτ ( mk ) is expressed as
dτ
769 770
(A14)
Following the plasticity theory of bounding surface, dγ p( mk ) is given by
764
767
1 dτ (mk) G (mk) e
(mk )
τ (mk ) = d pˆ = pˆ dr (mk ) + r(mk )dpˆ pˆ
(A17)
Combing equations (A16) and (A17) yields
dτ (mk ) = G(mk )dγ (mk ) + (1 − G(mk ) / H (mk ) )r (mk )dpˆ
771
(A18)
772
where H (mk ) and G (mk ) represent the total micro shear moduli related to d r ( mk ) and dpˆ for the
773
kth component in the mth microshear structure, respectively, define as
774
H
775 776
G
( mk )
1 1 = ( mk ) + ( mk ) Ge Gp
3
dε vd = ∑∑ 2w( m ) d1 (± rd( mk ) − r ( mk ) ){(1 − m =1 k =1
778
1 dpˆ 1 = ( mk ) + ( mk ) h ( pˆ − pˆ m ) Ge Hp dpˆ
−1
(A19)
−1
(A20)
From equations (19), (20), (A13), (A14), and (A18), dε vd is expressed as M
777
( mk )
G ( mk ) 1 G ( mk ) ( mk ) ( mk ) )d γ − (1 − )r dpˆ } Ge( mk ) Ge( mk ) H ( mk )
Substituting equations (A21) into equation (A11) yields
28
(A21)
M
3
dpˆ =K1{dε v − ∑∑ 2 w( m ) d1 (± rd( mk ) − r ( mk ) )(1 − G ( mk ) / Ge( mk ) ) d γ ( mk ) }
779
(A22)
m =1 k =1
780
where K1 is expressed as
K1 =
781
K M
3
1 − K ∑∑ 2 w d1 (± r ( m)
( mk ) d
m =1 k =1
782 783
(A23) )r
( mk )
(1 − G
( mk )
/H
( mk )
( mk ) e
)/G
From equations (1), (2), (A18), and (A22), the macro constitutive relation, including the effect of hydrate dissociation, is written as
dσ ij = Dijrs (dε rs − dε rsdis ) − dptδij
784 785
−r
( mk )
(A24)
where Dijrs is expressed as M
3
Dijrs = K1Qijδ rs + ∑∑ 2 w( m ) (−Q ( mk )Qij + G ( mk ) N ij( mk ) ) N rs( mk )
786
(A25)
m =1 k =1
787
in which M
788
3
Qij = δ ij + ∑∑ 2 w( m ) r ( mk ) (1 − G ( mk ) / H ( mk ) ) N ij( mk )
(A26)
Q ( mk ) = K1d1 ( ± rd( mk ) − r ( mk ) )(1 − G ( mk ) / Ge( mk ) )
(A27)
m =1 k =1
789 790 791
A.5. Relations of micro and macro parameters
792
By assuming that rb( mk ) , rc( mk ) , and rd( mk ) are the same for all microshear structures, they can be
793
determined through the macro model constants and yielding shape function as (Fang et al., 2017),
794
given by
ri ( mk ) =
795
2 3
M i g (θˆ) M
(i=b, c, d)
3
∑∑ 2w
(m)
m =1 k =1
N
(A28)
( mk ) 33
796 797 798 799
By assuming that all microshear structures have the same elastic shear modulus, Ge(mk ) is given by (Fang et al., 2017)
Ge( mk ) =
4 3
Ge M
3
∑∑ 2w
(m)
m =1 k =1
800
29
N
( mk ) 2 33
(A29)
Data Availability Datasets related to this article can be found at http.
References Bažant, Z.P., Oh, B.H., 1985. Microplane model for progressive fracture of concrete and rock. J. Eng. Mech. 111(4), 559–582. Boswell, R., Schoderbek, D., Collett, T.S., Ohtsuki, S., White, M., Anderson, B.J., 2017. The Iġnik Sikumi Field Experiment, Alaska North Slope: Design, operations, and implications for CO2-CH4 exchange in gas hydrate reservoirs. Energy & Fuels 31(1), 140–153. Been, K., Jefferies, M.G., 1985. A state parameter for sands. Géotechnique 35(2), 99–112. Caner, F.C., Bažant, Z.P., 2013. Microplane model M7 for plain concrete. I: Formulation. J. Eng. Mech. 139 (12), 1714–1723. Choi, J.H., Dai, S., Lin, J.H., Seol, Y., 2018. Multistage triaxial tests on laboratory-formed methane hydrate-bearing sediments. J. Geophys. Res. Solid Earth 123, 3347–3357. Clayton, C.R.I., Priest, J.A., Best, A.I., 2005. The effects of disseminated methane hydrate on the dynamic stiffness and damping of a sand. Géotechnique 55(6), 423–434. Clayton, C.R.I., Priest, J., Rees, E., 2010. The effects of hydrate cement on the stiffness of some sands. Géotechnique 60(6), 435–45. Dai, S., Santamarina, J.C., Waite, W.F., Kneafsey, T.J., 2012. Hydrate morphology: Physical properties of sands with patchy hydrate saturation. J. Geophys. Res. 117, B11205. Ebinuma, T., Kamata, Y., Minagawa, W., Ohuma, R., Nagao, J., Narita, H., 2005. Mechanical properties of sandy sediment containing methane hydrate. In Proceedings of the 5th International Conference on Gas Hydrate (pp. 958–961). Trondheim, Norway. Fang, H.L., Zheng, H., Zheng, J., 2017. Micromechanics-based multi-mechanism bounding surface model for sands. Int. J. Plast. 90, 242–266. Fang, H., Shen, Y., and Zhao, Y., 2019. Multishear bounding surface modelling of anisotropic sands accounting for fabric and its evolution. Comput. Geotech. 110, 57–70. Gai, X., Sánchez, M., 2017. A geomechanical model for gas hydrate-bearing sediments. Environ. Geotech. 4(2), 143-156. Ghiassian, H., Grozic, J.L.H., 2013. Strength behavior of methane hydrate bearing sand in undrained triaxial testing. Mar. Petrol. Geol. 43, 310–319.
30
Golkhou, F., Haghtalab, A., 2019. Measurement and thermodynamic modeling of carbon dioxide hydrate formation conditions using dry water through hydrophobic nano silica. J. Nat. Gas Sci. Eng. 68, 102906. Hajirezaie, S., Hemmati-Sarapardeh, A., Mohammadi, A.H., Pournik, M., Kamari, A., 2015. A smooth model for the estimation of gas/vapor viscosity of hydrocarbon fluids. J. Nat. Gas Sci. Eng. 26, 1452–1459. Hemmat Esfe, M., Tatar, A., Ahangar, M., Rostamian, H., 2018. A comparison of performance of several artificial intelligence methods for predicting the dynamic viscosity of TiOC2/SAE 50 nano-lubricant. Phys. E 96, 85–93. Hyodo, M., Yoneda, J., Yoshimoto, N., Nakata, Y., 2013a. Mechanical and dissociation properties of methane hydrate-bearing sand in deep seabed. Soils Found. 53, 299–314. Hyodo, M., Li, Y., Yoneda, J., Nakata, Y., Yoshimoto, N., Nishimura, A., Song, Y., 2013b. Mechanical behavior of gas-saturated methane hydrate-bearing sediments. J. Geophys. Res. Solid Earth 118, 5185–5194. Hyodo, M., Li, Y., Yoneda, J., Nakata, Y., Yoshimoto, N., Nishimura, A., 2014. Effects of dissociation on the shear strength and deformation behavior of methane hydrate-bearing sediments. Mar. Petrol. Geol. 51, 52–62. Hyodo, M., Wu, Y., Nakashima, K., Kajiyama, S., Nakata, Y., 2017. Influence of fines content on the mechanical behavior of methane hydrate-bearing sediments. J. Geophys. Res. Solid Earth 122, 7511–7524. Kajiyama, S., Hyodo, M., Nakata, Y., Yoshimoto, N., Wu, Y., Kato, A., 2017a. Shear behaviour of methane hydrate bearing sand with various particle characteristics and fines. Soils Found. 57, 176–193. Kajiyama, S., Wu, Y., Hyodo, M., Nakata, Y., Nakashima, K., Yoshimoto, N., 2017b. Experimental investigation on the mechanical properties of methane hydrate-bearing sand formed with rounded particles. J. Nat. Gas Sci. Eng. 45, 96–107. Kamalgharibi, M., Hormozi, F., Zamzamian, S.A.H., Sarafraz, M.M., 2016. Experimental studies on the stability of CuO nanoparticles dispersed in different base fluids: Influence of stirring, sonication and surface active agents. Heat Mass Transfer 52(1), 55–62. Kimoto, S., Oka, F., Tomohiko, F., 2010. A chemo-thermo-mechanically coupled analysis of ground deformation induced by gas hydrate dissociation. Int. J. Mech. Sci. 52(2), 365–376.
31
Klar, A., Soga, K., Ng, M.Y.A., 2010. Coupled deformation-flow analysis for methane hydrate extraction. Géotechnique 60(10), 765–776. Lade, P.V., Duncan, J.M., 1975. Elasto-plastic stress–strain theory for cohesionless soil. J. Geotech. Eng. 101(10), 1037–1053. Li, X.S., 2002. A sand model with state-dependent dilatancy. Géotechnique 52(3), 173–186. Lin, J.S., Seol, Y., Choi, J.H., 2015. An SMP critical state model for methane hydrate-bearing sands. Int. J. Numer. Anal. Meth. Geomech. 39(9), 969–987. Lin, J.S., Seol, Y., Choi, J.H., 2017. Geomechanical modeling of hydrate-bearing sediments during dissociation under shear. Int. J. Numer. Anal. Meth. Geomech. 41, 1523–1538. Liu, W., Zhao, T, Luo, Y., Li, Y., Song, Y., Zhu, Y., Liu, Y., Zhao, J., Wu, Z., Xu, X., 2016. Experimental study on the mechanical properties of sediments containing CH4 and CO2 hydrate mixtures. J. Nat. Gas Sci. Eng. 32, 20–27. Liu, Z., Dai, S., Ning, F., Peng, L., Wei, H., Wei, C., 2017. Strength estimation for hydrate-bearing sediments from direct shear tests of hydrate-bearing sand and silt. Geophys. Res. Lett. 45, 715– 723. Malagar, B.R.C., Lijith, K.P., Singh, D.N., 2019. Formation & dissociation of methane gas hydrates in sediments: A critical review. J. Nat. Gas Sci. Eng. 65, 168–184. Masui, A., Haneda, H., Ogata, Y., Aoki, K., 2005. The effects of methane hydrate formation on shear strength of synthetic methane hydrate sediments. In Proceedings of the 5th International Conference on Gas Hydrates (pp. 657–663). Trondheim, Norway. Masui, A., Miyazaki, K., Haneda, H., Ogata, Y., Aoki, K., 2008. Mechanical characteristics of natural and artificial gas hydrate bearing sediments. In Proceedings of the 6th International Conference on Gas Hydrates. Chevron, Vancouver, B. C., Canada. Matsuoka, H., Nakai, T., 1974. Stress-deformation and strength characteristics of soil under three difference principal stresses. Proc. JSCE 232, 59–70. Matsuoka, H., Yao, Y.P., Sun, D.A., 1999. The Cam-clay models revised by the SMP criterion. Soils Found. 39(1), 81–95. Miyazaki, K., Masui, A., Sakamoto, Y., Aoki, K., Tenma, N., Yamaguchi, T., 2011. Triaxial compressive properties of artificial methane-hydrate-bearing sediment. J. Geophys. Res. 116, B06102.
32
Miyazaki, K., Tenma, N., Aoki, K., Yamaguchi, T., 2012. A nonlinear elastic model for triaxial compressive properties of artificial methane-hydrate-bearing sediment samples. Energies 5(10), 4057–4075. Nikkhah, V., Sarafraz, M.M., Hormozi, F., 2015a. Application of spherical copper oxide (II) water nano-fluid as a potential coolant in a boiling annular heat exchanger. Chem. Biochem. Eng. Q. 29(3), 405–415. Nikkhah, V., Sarafraz, M.M., Hormozi, F., Peyghambarzadeh, S.M., 2015b. Particulate fouling of CuO-water nanofluid at isothermal diffusive condition inside the conventional heat exchanger-experimental and modeling. Exp. Therm. Fluid Sci. 60, 83–95. Nguyen, L.D., Fatahi, B., Khabbaz, H., 2014. A constitutive model for cemented clays capturing cementation degradation. Int. J. Plast. 56, 1–18. Pinkert, S., Grozic, J., 2014. Prediction of the mechanical response of hydrate-bearing sands. J. Geophys. Res. Solid Earth 119(6), 4695–4707. Pinkert, S., Grozic, J., Priest, J., 2015. Strain-softening model for hydrate-bearing sands. Int. J. Geomech. 15(6), 04015007. Priest, J.A., Best, A.I., Clayton, C.R.I., 2005. A laboratory investigation into the seismic velocities of methane gas hydrate-bearing sand. J. Geophys. Res. 110, B04102. Priest, J.A., Rees, E.V.L., Clayton, C.R.I., 2009. Influence of gas hydrate morphology on the seismic velocities of sands. J. Geophys. Res. Solid Earth 114, B11205. Richart, F., Hall, J, Woods, R., 1970. Vibrations of soils and foundations. International Series in Theoretical and Applied Mechanics. Englewood Cliffs, NJ, USA: Prentice-Hall. Rostamian, H., Lotfollahi, M.N., 2016. New functionality for energy parameter of redlich-kwong equation of state for density calculation of pure carbon dioxide and ethane in liquid, vapor and supercritical phases. Period. Polytech., Chem. Eng. 60(2), 93–97. Rostamian, H., Lotfollahi, M.N., 2019. A novel statistical approach for prediction of thermal conductivity of CO2 by response surface methodology. Phys. A 527, 121175. Sánchez, M., Gai, X., Santamarina, J.C., 2017. A constitutive mechanical model for gas hydrate bearing sediments incorporating inelastic mechanisms. Comput. Geotech. 84, 28–46. Sánchez, M., Santamarina, C., Teymouri, M., Gai, X., 2018. Coupled numerical modeling of gas hydrate-bearing sediments: From laboratory to field-scale analyses. J. Geophys. Res. Solid Earth 123(10), 326–348.
33
Santamarina, J.C., Dai, S., Terzariol, M., Jang, J., Waite, W.F., Winters, W.J., Suzuki, K., 2015. Hydro-bio-geomechanical properties of hydrate-bearing sediments from Nankai Trough. Mar. Petrol. Geol. 66, 434–450. Sasaki, T., Soga, K., Elshafie, M.E.B., 2018. Simulation of wellbore construction in offshore unconsolidated methane hydrate-bearing formation. J. Nat. Gas Sci. Eng. 60, 312–326. Shen, J., Chiu, C.F., Ng, C.W.W., Lei, G.H., Xu, J., 2016. A state-dependent critical state model for methane hydrate-bearing sand. Comput. Geotech. 75, 1–11. Song, B., Cheng, Y., Yan, C., Ly, Y., Wei, J., Ding, J., Li, Y., 2019. Seafloor subsidence response and submarine slope stability evaluation in response to hydrate dissociation. J. Nat. Gas Sci. Eng. 65, 197–211. Sultan N., Garziglia S., 2011. Geomechanical constitutive modelling of gas-hydrate bearing sediments. In Proceedings of the 7th International Conference on Gas Hydrates. Edinburgh, Scotland, United Kingdom. Sultaniya, A., Priest, J.A., Clayton, C.R.I., 2018. Impact of formation and dissociation conditions on stiffness of a hydrate-bearing sand. Can. Geotech. J. 55, 988–998. Sun, D., Huang, W., Yao, Y., 2008. An experimental study of failure and softening in sand under three-dimensional stress condition. Granular Matter 10, 187–195. Sun, X., Guo, X., Shao, L., Tang, H., 2015. A thermodynamics-based critical state constitutive model for methane hydrate bearing sediment. J. Nat. Gas Sci. Eng. 27, 1024–1034. Taylor, G.I., 1938. The plastic strain of metals. J. Inst. Metals 62, 307–324. Thevanayagam, S., Shenthan, T., Mohan, S., Liang, J., 2002. Undrained fragility of clean sands, silty sands, and sandy silts. J. Geotech. Geoenviron. Eng. 128(10), 849–859. Uchida, S., Soga, K., Yamamoto, K., 2012. Critical state soil constitutive model for methane hydrate soil. J. Geophys. Res. 117, B03209. Uchida, S., Xie, X.G., Leung, Y.F., 2016. Role of critical state framework in understanding geomechanical behavior of methane hydrate-bearing sediments. J. Geophys. Res. Solid Earth 121(8), 5580–5595. Wang, Z.L., Dafalias, Y.F., Shen, C.K., 1990. Bounding surface hypoplasticity model for sand. J. Eng. Mech. 116(5), 983–1001. Winters, W.J., Waite, W.F., Mason, D.H., Gilbert, L.Y., Pecher, I.A., 2007. Methane gas hydrate effect on sediment acoustic and strength properties. J. Petrol. Sci. Eng. 56(1–3), 127–135.
34
Yan, C., Li, Y., Cheng, Y., Wang, W., Song, B., Deng, F., Feng, Y., 2018. Sand production evaluation during gas production from natural gas hydrates. J. Nat. Gas Sci. Eng. 57, 77–88. Yan, G., Cheng, Y., Li, M., Han, Z., Zhang, H., Li, Q., Teng, F., Ding, J., 2017. Mechanical experiments and constitutive model of natural gas hydrate reservoirs. Int. J. Hydrogen Energy 42, 19810–19818. Yan, Y., Wei, C., 2017. Constitutive model for gas hydrate-bearing soils considering hydrate occurrence habits. Int. J. Geomech. 17(8), 04017032. Yao, Y., Sun, D., 2000. Application of Lade’s criterion to Cam-clay model. J. Eng. Mech. 126(1), 112–119. Yoneda, J., Masui, A., Konno, Y., Jin, Y., Egawa, K., Kidab, M., Ito, T., Nagao, J., Tenma, N., 2015a. Mechanical behavior of hydrate-bearing pressure-core sediments visualized under triaxial compression. Mar. Petrol. Geol. 66(2), 451–459. Yoneda, J., Masui, A., Konno, Y., Jin, Y., Egawa, K., Kida, M., Ito, T., Nagao, J., Tenma, N., 2015b. Mechanical properties of hydrate-bearing turbidite reservoir in the first gas production test site of the Eastern Nankai Trough. Mar. Petrol. Geol. 66(2), 471–486. Yoneda, J., Jin, Y., Katagiri, J., Tenma, N., 2016. Strengthening mechanism of cemented hydrate-bearing sand at microscales. Geophys. Res. Lett. 43(14), 7442–7450. Yoneda, J., Oshima, M., Kida, M., Kato, A., Konno, Y., Jin, Y., Jang, J., Waite, W.F., Kumar, P., Tenma, N., 2019a. Pressure core based onshore laboratory analysis on mechanical properties of hydrate-bearing sediments recovered during India's National Gas Hydrate Program Expedition (NGHP) 02. Mar. Petrol. Geol. (In press) Yoneda, J., Oshima, M., Kida, M., Kato, A., Konno, Y., Jin, Y., Tenma, N., 2019a. Consolidation and hardening behavior of hydrate-bearing pressure-core sediments recovered from the Krishna–Godavari Basin, offshore India. Mar. Petrol. Geol. (In press) Yoneda, J., Oshima, M., Kida, M., Kato, A., Konno, Y., Jin, Y., Jang, J., Waite, W.F., Kumar, P., Tenma, N., 2019b. Pressure core based onshore laboratory analysis on mechanical properties of hydrate-bearing sediments recovered during India's National Gas Hydrate Program Expedition (NGHP) 02. Mar. Petrol. Geol. (In press) Yu, T., Guan, G., Abudula, A., Yoshida, A., Wang, D., Song, Y., 2019. Application of horizontal wells to the oceanic methane hydrate production in the Nankai Trough, Japan. J. Nat. Gas Sci. Eng. 62, 113–131.
35
Yun, T.S., Santamarina, C.J., Ruppel, C., 2007. Mechanical properties of sand, silt, and clay containing tetrahydrofuran hydrate. J. Geophys. Res. 112, B04106. Zhang, X., Lin, J., Lu, X., Liu, L., Liu, C., Li, M., Su, Y., 2017. A hypoplastic model for gas hydrate-bearing sandy sediments. Int. J. Numer. Anal. Meth. Geomech. 42, 931–942. Zhao, J., Guo, N., 2013. Unique critical state characteristics in granular media considering fabric anisotropy. Géotechnique 63(8), 695–704. Zhou, M., Soga, K., Yamamoto, K., 2018. Upscaled anisotropic methane hydrate critical state model for turbidite hydrate-bearing sediments at East Nankai Trough. J. Geophys. Res. Solid Earth 123(8), 6277–6298.
36
Table 1 Model parameters Test name Triaxial compression tests for synthetic pore-filling type of GHBS (Masui et al., 2005)
Parameters for host materials G0 =150, κ0 =0.005
Triaxial compression tests for synthetic cementing type of GHBS (Masui et al., 2005)
G0 =150, κ0 =0.005
Triaxial compression tests for synthetic GHBS (Hyodo et al., 2013a)
G0 =200, κ0 =0.004
M c0 =1.47, eΓ 0 =1.126, λ c =0.132 h1 =0.37, h2 =0.35, nb =1.5, λ0 =0.01 d1 =0.7, nd =2.3 M c0 =1.47, eΓ 0 =1.126, λ c =0.132 h1 =0.66, h2 =0.35, nb =1.5, λ0 =0.01
Additional parameters for GHBS aG =3.5, a M =0.038 ac =0.04, nc =0.53 for S h < 0.3 ac =0.39, nc =0.53 for S h ≥ 0.3 a t =0.15MPa, k t =50
aG =2.66, a M =0.071 ac =0.6, nc =2.0 a t =0.7MPa, k t =50
d1 =0.7, nd =4.0 M c0 =1.20, eΓ 0 =1.126, λ c =0.132 h1 =0.55, h2 =0.35, nb =0.8, λ0 =0.008
aG =2.66, a M =0.226 ac =0.6, nc =2.0 a t =0.43MPa, k t =50
d1 =0.6, nd =3.5 Triaxial compression tests for natural GHBS (Masui et al., 2008)
G0 =200, κ0 =0.004 M c0 =1.48, eΓ 0 =1.126, λ c =0.132 h1 =0.55, h2 =0.35, nb =1.0, λ0 =0.008
aG =2.66, a M =0.071 ac =0.37, nc =0.0 a t =0.6MPa, k t =50
d1 =0.7, nd =1.5 Triaxial compression tests for natural GHBS (Yoneda et al., 2015b)
G0 =100, κ0 =0.0075 M c0 =1.44, eΓ 0 =1.126, λ c =0.132 h1 =0.55, h2 =0.35, nb =3.5, λ0 =0.015
aG =1.5, a M =0.157 ac =0.12, nc =0.3 a t =0.6MPa, k t =80
d1 =0.7, nd =3.0 Triaxial compression tests for natural GHBS (Yoneda et al., 2019a, 2019b )
G0 =100, κ0 =0.0075 M c0 =1.2, eΓ 0 =0.8, λ c =0.12 h1 =0.2, h2 =0.35, nb =1.0, λ0 =0.0225 d1 =2.0, nd =1.0
37
aG =1.4, a M =0.125 ac =452.0, nc =13.5 for S h < 0.6 ac =0.13, nc =-3.4 for S h ≥ 0.6 a t =1.7MPa, k t =100
Fig. 1 Illustration of conceptual framework for multishear model
Fig. 2 Schematic diagram of different morphologies for GHBS
38
1.2
Solid lines: regressions Test (Sh=0) Test (Sh=0.242-0.287) Test (Sh=0.311-0.387) Sh=0.439 (average) Test (Sh=0.419-0.450) ec=0.939-0.132lnp 2 Test (Sh=0.504-0.572) R =0.75
Sh=0.532 (average) ec=0.993-0.132lnp 2
R =0.23
Void ratio, e
1.0
0.8 Sh=0 ec=0.823-0.132lnp Sh=0.262 (average) ec=0.865-0.132lnp
2
0.6
R =0.76
2
R =0.55
Sh=0.342 (average) ec=0.894-0.132lnp 2
R =0.68 0.4 1
2
3
4
5
6
7
8
9 10
Mean effective stress, p (MPa)
Fig. 3 Critical state line in e − ln p plane (test data from Hyodo et al., 2013a)
2.0
Critical state void ratio parameter, eΓ
Regression Obtained by fitting lines in Fig.3
1.5
1.0
2
eΓ=1.126+0.6Sh 2
R =0.85 0.5
0.0 0.0
0.2
0.4
0.6
0.8
Hydrate saturation, Sh
Fig. 4 Relationship between critical state void ratio parameter and hydrate saturation (test data from Hyodo et al., 2013a)
39
1.6
Regression Test (pp=10MPa, T= 5 Test (pp= 5MPa, T= 5 Test (pp=15MPa, T=10 Test (pp=15MPa, T= 1
Critical stress ratio, Mc
1.4
) ) ) )
1.2 2
Mc=1.2(1+0.226Sh), R =0.32 1.0
0.8 0.0
0.2
0.4
0.6
0.8
Hydrate saturation, Sh
Fig. 5 Relationship between critical stress ratio and hydrate saturation (test data from Hyodo et al., 2013a)
15
Solid lines: regressions Test (Sh=0.242-0.296) Test (Sh=0.414-0.457) Test (Sh=0.504-0.572)
Deviator stress, q (MPa)
12
Sh=0.532 (average) qf=1.42p+0.35
9
Sh=0.432 (average) qf=1.38p+0.26
2
R =0.95
2
R =0.99
6 Sh=0.267 (average) qf=1.29p+0.11 3
2
R =0.92
0 -2
0
2
4
6
8
10
12
Mean effective stress, p (MPa)
Fig. 6 Peak strength lines of GHBS with different hydrate saturations (test data from Hyodo et al., 2013a)
40
0.4
Initial bonding strength, pt0 (MPa)
Regression Obtained by fitting lines in Fig. 6
0.3
2
pt0=0.43Sh, R =0.96 0.2
0.1
0.0 0.0
0.2
0.4
0.6
0.8
Hydrate saturation, Sh
Fig. 7 Relationship between initial bonding strength and hydrate saturation (test data from Hyodo et al., 2013a)
Fig. 8 Characteristic surfaces in macro stress space
41
Fig. 9 Characteristic lines in τ ( mk ) − pˆ plane
Elastic shear modulus parameter, Gh0
800 Regression Test
600
Gh0=200(1+2.66Sh) 2
R =0.61 400
200
0 0.0
0.2
0.4
0.6
0.8
Hydrate saturation, Sh
Fig. 10 Relationship between elastic shear modulus parameter and hydrate saturation (test data from Hyodo et al., 2013a)
42
10 Simulation Test Sh e0 None 0.678 0.736 None 0.501 0.736 0.409 0.736 None 0.264 0.736
Deviator stress, q (MPa)
8
pc0 (MPa) 1.0 1.0 1.0 1.0
6
4
2
Pore pressure: 8 MPa Temperature: 5
0 0
2
4
6
8
10
12
14
16
Axial strain, εa (%)
-10
Simulation Test Sh e0 None 0.678 0.736 None 0.501 0.736 0.409 0.736 None 0.264 0.736
Volumetric strain, εv (%)
-8
-6
pc0 (MPa) 1.0 1.0 1.0 1.0
Pore pressure: 8 MPa Temperature: 5
-4
-2
0
2 0
2
4
6
8
10
12
14
16
Axial strain, εa (%)
Fig. 11 Test and simulation results for triaxial compression tests of synthetic pore-filling type of GHBS (test data from Masui et al., 2005)
43
10
Simulation Test Sh None 0.551 0.407 None 0.257 0.000
Deviator stress, q (MPa)
8
e0 0.605 0.605 0.605 0.605
pc0 (MPa) 1.0 1.0 1.0 1.0
6
4
2
Pore pressure: 8 MPa Temperature: 5
0 0
2
4
6
8
10
12
14
16
Axial strain, εa (%)
-16
Simulation Test Sh None 0.551 0.407 None 0.257 0.000
Volumetric strain, εv (%)
-12
e0 0.605 0.605 0.605 0.605
pc0 (MPa) 1.0 1.0 1.0 1.0
-8
-4
0
Pore pressure: 8 MPa Temperature: 5
4 0
2
4
6
8
10
12
14
16
Axial strain, εa (%)
Fig. 12 Test and simulation results for triaxial compression tests of synthetic cementing type of GHBS (test data from Masui et al., 2005)
44
16
Deviator stress, q (MPa)
14 12 10 8 6 Simulation Test 4 Pore pressure: 10 MPa Temperature: 5
2
Sh 0.531 0.351 0.242 0.000
e0 0.669 0.645 0.565 0.650
pc0 (MPa) 5.0 5.0 5.0 5.0
0 0
2
4
6
8
10
12
14
16
12
14
16
Axial strain, εa (%)
-8 Simulation Test
Volumetric strain, εv (%)
-6 -4
Sh 0.531 0.351 0.242 0.000
e0 0.669 0.645 0.565 0.650
pc0 (MPa) 5.0 5.0 5.0 5.0
-2 0 2 4 Pore pressure: 10 MPa Temperature: 5
6 8 0
2
4
6
8
10
Axial strain, εa (%)
Fig. 13 Test and simulation results for triaxial compression tests of synthetic GHBS (test data from Hyodo et al., 2013a)
45
Simulation Test
20
e0
0.531 0.537 0.543 0.523 0.516
16
Deviator stress, q (MPa)
Sh
pc0 pp T (MPa) (MPa) ( ) 0.669 5.0 10.0 5.0 0.650 3.0 10.0 5.0 0.650 1.0 10.0 5.0 0.645 3.0 15.0 1.0 0.653 3.0 15.0 10.0
12
8
4
0 0
2
4
6
8
10
12
14
16
12
14
16
Axial strain, εa (%)
-10
Simulation Test
e0
0.531 0.537 0.543 0.523 0.516
-8
Volumetric strain, εv (%)
Sh
-6
pc0 pp T (MPa) (MPa) ( ) 0.669 5.0 10.0 5.0 0.650 3.0 10.0 5.0 0.650 1.0 10.0 5.0 0.645 3.0 15.0 1.0 0.653 3.0 15.0 10.0
-4
-2
0
2
4 0
2
4
6
8
10
Axial strain, εa (%)
Fig. 14 Test and simulation results for triaxial compression tests of synthetic GHBS (test data from Hyodo et al., 2013a)
46
7
Deviator stress, q (MPa)
6
5
4
3
2
Simulation Test Pore pressure: 9 MPa Temperature: 5
1
Sh 0.376 0.225 0.094 0.077
e0 0.792 0.961 0.972 1.033
pc0 (MPa) 1.0 1.0 1.0 1.0
0 0
2
4
6
8
10
12
14
16
12
14
16
Axial strain, εa (%)
-8
Simulation Test
Volumetric strain, εv (%)
-6
Sh 0.376 0.225 0.094 0.077
pc0 (MPa) 1.0 1.0 1.0 1.0
e0 0.792 0.961 0.972 1.033
Pore pressure: 9 MPa Temperature: 5
-4
-2
0
2 0
2
4
6
8
10
Axial strain, εa (%)
Fig. 15 Test and simulation results for triaxial compression tests of natural GHBS (test data from Masui et al., 2008)
47
14 Simulation Test Sample No. 9 No. 7 No. 9-2 No. 7-2
Deviator stress, q (MPa)
12
10
Sh 0.79 0.38 0.00 0.00
e0 pc0 (MPa) 0.650 1.6 0.790 1.5 0.748 1.6 0.887 1.5
8
6
4
2
Pore pressure: 13 MPa Temperature: 10
0 0
2
4
6
8
10
12
14
16
12
14
16
Axial strain, εa (%)
-9 Simulation Test Sample No. 9 No. 7 No. 9-2 No. 7-2
Volumetric strain, εv (%)
-6
Sh 0.79 0.38 0.00 0.00
e0 pc0 (MPa) 0.650 1.6 0.790 1.5 0.748 1.6 0.887 1.5
-3
0
3
6
Pore pressure: 13 MPa Temperature: 10
9 0
2
4
6
8
10
Axial strain, εa (%)
Fig. 16 Test and simulation results for triaxial compression tests of natural GHBS (test data from Yoneda et al., 2015b)
48
0.9
Simulation Test Sample Depth (m) Sh 16B-3P 273.81-273.89 0.321 16B-4P 277.22-277.30 0.716 16B-3P 274.13-274.23 0.601
0.8
Before dissociation
Void ratio, e
0.7
During dissociation 0.6
After dissociation
0.5
0.4
0.3 0.1
1
10
100
Effective confining stress, pc (MPa)
(a) Hydrate dissociation under effective confining stress of 1 MPa 0.9
0.8 Before dissociation
Void ratio, e
0.7
0.6
During dissociation After dissociation
0.5
0.4
0.3 0.1
Simulation Test Sample Depth (m) Sh 16B-7P 287.12-277.20 0.639 16B-3P 274.03-274.13 0.606 23C-11P 285.25-285.33 0.556 1
10
100
Effective confining stress, pc (MPa)
(b) Hydrate dissociation under effective confining stress of 5 MPa Fig. 17 Test and simulation results for isotropic consolidation tests of natural GHBS (test data from Yoneda et al., 2019a)
49
5
Pore pressure: 10 MPa Temperature: 4
Deviator stress, q (MPa)
4
3
2 Simulation Test 1
Sh 0.822 0.734 0.634 0.552 0.487
e0 0.656 0.658 0.825 0.634 0.712
pc0 (MPa) 1.1 1.0 1.0 1.0 1.0
0 0
2
4
6
8
10
12
14
16
12
14
16
Axial strain, εa (%)
-12
Simulation Test
Sh 0.822 0.734 None 0.634 0.552 0.487
Volumetric strain, εv (%)
-8
e0 pc0 (MPa) 0.656 1.1 0.658 1.0 0.825 1.0 0.634 1.0 0.712 1.0
-4
0
4
8
Pore pressure: 10 MPa Temperature: 4
12 0
2
4
6
8
10
Axial strain, εa (%)
Fig. 18 Test and simulation results for triaxial compression tests of natural GHBS (test data from Yoneda et al., 2019b)
50
εa
1.5
εθ
εaθ
εr Sh=0.0 Sh=0.3 Sh=0.6
Strains, εa, εθ, εaθ, εr (%)
1.0
0.5
0.0
-0.5 e0=0.65 pc0=1MPa
-1.0
-1.5 0
30
60
90
120
150
180
Principal stress direction, α (°)
Fig. 19 Predicted strain responses with principal stress direction for drained pure rotational shear tests
51
1.0
e0=0.65 pc0=1MPa
Strain increment Sh=0.0 Sh=0.3 Sh=0.6
σaθ (MPa), 2dεaθ/ds (%)
0.5
0.0
-0.5
1% -1.0 -1.0
-0.5
0.0
0.5
1.0
(σa-σθ)/2 (MPa), (dεa-dεθ)/ds (%)
Fig. 20 Predicted stress paths and strain increment vectors for drained pure rotational shear tests
1.0
e0=0.65 pc0=1MPa
Plastic strain increment Sh=0.0 Sh=0.3 Sh=0.6
p
σaθ (MPa), 2dεaθ/ds (%)
0.5
0.0
-0.5
0.2% -1.0 -1.0
-0.5
0.0
0.5
1.0
(σa-σθ)/2 (MPa), (dε -dεθ )/ds (%) p a
p
Fig. 21 Predicted stress paths and plastic strain increment vectors for drained pure rotational shear tests
52
Highlights •
A geomechanical constitutive model in the multishear bounding surface framework is proposed for gas hydrate-bearing sediments.
•
A macro constitutive response is decomposed into a macro volume response and a series of micro shear responses in spatial distributions.
•
The model captures the consolidation, hardening, softening, dilatation, collapse, and non-coaxial behavior of gas hydrate-bearing sediments.
1
Declaration of interests ☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. ☒The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:
None