- Email: [email protected]
Nonlinear viscoelastic medium equivalence for stress wave propagation in a jointed rock mass
International Journal of Rock Mechanics & Mining Sciences 50 (2012) 11–18
Contents lists available at SciVerse ScienceDirect
International Journal of Rock Mechanics & Mining Sciences journal homepage: www.elsevier.com/locate/ijrmms
Nonlinear viscoelastic medium equivalence for stress wave propagation in a jointed rock mass L.F. Fan a, G.W. Ma b,n, J.C. Li c a
School of Civil and Environmental Engineering, Nanyang Technological University, Singapore 639798, Singapore School of Civil and Resource Engineering, The University of Western Australia, Crawley, WA 6009, Australia c Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, Wuhan 430071, China b
a r t i c l e i n f o
abstract
Article history: Received 16 August 2010 Received in revised form 26 November 2011 Accepted 17 December 2011 Available online 30 January 2012
A nonlinear viscoelastic equivalent medium model is proposed for longitudinal stress wave propagation in a rock mass that contains equally spaced parallel joints. The joint deformation is assumed to satisfy the nonlinear Bandis–Barton (B–B) deformational model. By the simplification of the nonlinear B–B model into piecewise linear segments, the material parameters in the equivalent medium model are analytically derived. The nonlinear viscoelastic equivalent medium model is verified by comparison with the results based on the displacement discontinuity method. It is shown that for smaller amplitude stress waves, when the rock joints deform in the linear elastic region, the nonlinear viscoelastic equivalent medium model is reduced to its linear form, whereas for larger amplitude stress waves, the present nonlinear model describes the wave amplitude-dependence of wave velocity. Two important properties, i.e., the wave transmission coefficient and the effective wave propagation velocity, are discussed in detail. & 2011 Elsevier Ltd. All rights reserved.
Keywords: Equivalent medium method Nonlinear deformational joint Jointed rock mass Viscoelastic medium Wave propagation
1. Introduction The presence of rock joints often dominates the dynamic engineering properties and increases the complexity in modeling dynamic behavior of a jointed rock mass [1,2]. It is of interest to develop an efficient and accurate analytical model to describe the dynamic behavior of jointed rock masses and evaluate the effects of joints on stress wave propagation. It has been observed that the complete deformational behavior of a rock joint is generally nonlinear [2–6]. The linear elastic joint model can be regarded as a special case of a nonlinear model such as the B–B model [2] when the joint is dry and the magnitude of the stress wave is too small to mobilize the joint nonlinear deformation. Currently, there are two different theoretical approaches to evaluate the effect of rock joints on the stress wave propagation through a rock mass [1,7,8]. One is the displacement discontinuity method (DDM) [3,8–13] and the other is the equivalent medium method (EMM) [1,7,14–20]. In the DDM, the joints are assumed to be non-welded interfaces. The displacements across the joint interface are discontinuous, while the stresses across the joint interface are continuous. These two assumptions are used as two boundary conditions in the stress wave propagation equations. By using the DDM, a series of analyses of stress wave propagation
n
Corresponding author. E-mail address: [email protected] (G.W. Ma).
1365-1609/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijrmms.2011.12.008
through rock joints have been carried out [11,21,22], and it has been applied successfully to predict the seismic responses of jointed rock masses [3,11–13]. Since the DDM treats the rock mass as a discrete rock medium, it is difficult to give an explicit expression of stress wave propagation in a jointed rock mass, and it becomes very tedious and complex when the number of discontinuities increases [1]. On the other hand, the equivalent medium method treats the rock mass as a whole and predicts the effects of joint on the behavior of a rock mass system using a representative elementary volume (REV), so as to make a continuum medium analysis of practical problems. An elastic continuum medium is commonly used to represent the jointed rock mass [16–19]. Such conventional linear elastic equivalence is effective only if the stress wave frequency and amplitude effects, and the multiple wave reflections among the joints are ignored [1,8,23]. Li et al. [1] introduced a viscoelastic equivalent medium model for rock masses consisting of equally-spaced parallel joints. A virtual wave source (VWS) concept has been put forward to simulate the multiple wave reflections among the joints [1,7]. The wave dispersion, attenuation and frequency dependence of a jointed rock mass have been analytically investigated by solving an explicit stress wave propagation equation. However, the joints were assumed to deform linearly in the previous study. A nonlinear equivalent medium method is necessary in order to efficiently analyze stress wave propagation in a jointed rock mass where the rock joint deforms nonlinearly.
12
L.F. Fan et al. / International Journal of Rock Mechanics & Mining Sciences 50 (2012) 11–18
This paper proposes a nonlinear viscoelastic equivalent medium model for the study of stress wave propagation through a rock mass with equally spaced parallel rock joints. The deformation of the rock joints satisfies the B–B model [2]. The proposed equivalent medium model is an extended standard solid model with three frequency dependent piecewise-linear parameter sets. In order to derive the wave propagation equation in the viscoelastic medium efficiently, a method of characteristic lines is adopted. To theoretically verify the present model, a comparison with the results obtained from the displacement discontinuity method is performed. The amplitude dependence, which only exists in stress wave propagation through a nonlinear medium, and the frequency dependence of the present viscoelastic model are discussed.
2. Nonlinear viscoelastic medium equivalence
d¼
Fig. 1 illustrates an REV of a jointed rock mass and the definition of the corresponding equivalent medium. In the present study, the joint is assumed to be planar, large in extent and small in thickness. The length of the REV of the equivalent medium is defined to be the joint spacing S as shown in Fig. 1. The normal deformational behavior of a rock joint under quasi-static and cyclic loadings has been studied [2,5]. Fig. 2 illustrates a piecewise linear model to consider the nonlinear deformational behavior of a rock joint, which is compared with the linear joint model and the B–B model [2]. If the closure (opening) of the joint and the compression (tension) stress are assumed to be positive (negative), the normal constitutive relations between the stress and the closure can be described by the
s
ð1Þ
kn þ s=dmax
where d is the closure of the joint, s is the normal stress, dmax is the maximum allowable closure of the joint and kn is the normal stiffness of the joint at the initial state. The dashed curve in Fig. 2 gives a typical hyperbolic fitting of the B–B model [11]. A piecewise linear model is proposed to approximate such hyperbolic B–B model. As shown in Fig. 2 that the curve of the B–B model is a continuous hyperbolic dashed line, while the present piecewise linear model can be regarded as a polygonal curve consisting of a series of linear segments with their end points located on the hyperbolic B–B model. Therefore, the values of the closure and the stress at the two ends of each segment are assumed to satisfy the hyperbolic B–B model relation yields di ¼
2.1. Joint models
si kn þ si =dmax
,
i ¼ 1,2,3,
ð2Þ
where i is the segment number, di and si denote the discrete threshold closure and stress at the two ends of each segment, respectively. The difference of the stress and closure for every individual segment can be obtained as
Dsi ¼ si si1 , i ¼ 1,2,3,
ð3Þ
and
Ddi ¼ di di1 , i ¼ 1,2,3,
ð4Þ
, respectively. Therefore, the stiffness of the joint for segment i can be obtained from Eqs. (2)–(4) as Dsi kn di di1 , i ¼ 1,2,3, ki ¼ ¼ Ddi di di1 1ðdi =dmax Þ 1ðdi1 =dmax Þ ð5Þ Considering Eqs. (2) and (5), the stress–closure relation of the present piecewise linear model is described as
Joint Rock
d ¼ di1 þ
a
b
a
b
S
S
Fig. 1. Representative elementary volume. (a) Jointed rock mass; (b) Equivalent medium.
14
Piecewise linear model BB model Linear model
12 Normal stress σ (MPa)
B–B model as,
d
10
(di,σi)
8 ki
6
(di-1,σi-1)
Δσ
ssi1 ki
,
di1 rd o di
ð6Þ
where di is obtained from Eq. (2) and ki is obtained from Eq. (5). It can be obtained from Fig. 2 that for given stress s, the B–B model reveals that the closure will increase following the locus of the dashed hyperbolic curve to d, while the present piecewise model reveals that the closure will increase following the locus of the solid linear segments to d. When the closure difference Ddi is small enough, the present piecewise linear model approximates the B–B model with sufficient accuracy. The values of dmax and kn for the piecewise linear model and the B–B model can be determined from the laboratory measurements as obtained by Bandis et al. [2] and Barton et al. [24]. It can also be realized from Fig. 2 that the present piecewise linear model and the B–B model are able to describe the joint stiffness increase with the increase of the normal stress, however the joint stiffness of a linear model always keeps a constant. 2.2. Nonlinear viscoelastic equivalent medium method
Δd
4 2 0 0.0
0.1
0.2
0.3
0.4 0.5 0.6 Closure d (mm)
0.7
0.8
0.9
Fig. 2. Piecewise linear model, B–B model and linear model of a rock joint.
A nonlinear viscoelastic equivalent model is proposed to describe the overall dynamic mechanical behavior of a rock mass containing nonlinear deformational joints, which is shown in Fig. 3. A piecewise linear spring with a modulus of Evi is in parallel with a piecewise linear dashpot Zvi, then in series with another piecewise linear spring Eai. The piecewise linear quantities of Eai, Evi and Zvi are related to the stress–closure relation of the joint at segment i. For the ith segment, Eai, Evi and Zvi can be obtained from the coefficients
L.F. Fan et al. / International Journal of Rock Mechanics & Mining Sciences 50 (2012) 11–18
13
where the sign ‘‘ þ’’ means the rightward running waves and ‘‘ ’’ denotes the leftward running waves. The second solution of Eq. (11) is
Evi Eai
L¼M¼0 N a0
ð15Þ
Therefore, the characteristic line of the second solution is
ηvi
dx ¼ 0
ð16Þ
Fig. 3. Piecewise linear viscoelastic equivalent medium model.
of the phase shift and amplitude attenuation of stress wave propagation [1]. The constitutive equation of the nonlinear model shown in Fig. 3 is given as ðEai þ Evi Þs þ Zvi
@s @e Zvi Eai Eai Evi e ¼ 0 @t @t
ð7Þ
where e is the strain component, Eai and Evi are the piecewise elastic modulus of the springs, respectively, Zvi is the piecewise viscosity of the dashpot, t denotes the time variable, and Eai, Evi and Zvi vary with the change of the magnitude of joint stiffness ki. The motion equation of one-dimensional longitudinal wave propagation is governed by
r
@v @s ¼ @t @x
ð8Þ
where r is the density of the medium material, x is the Cartesian coordinate along the wave propagation path and v is the particle velocity at x. In addition to the constitutive Eq. (7) and the equation of motion (8), the continuity equation is known as @v @e ¼ @x @t
The compatibility condition along the characteristic line of Eq. (16) is de
ds ðE þEvi ÞsEai Evi e ai dt ¼ 0 rC 2vi rZvi C 2vi
ð17Þ
From Eqs. (13) and (16), the wave propagation in the viscoelastic medium can be described by the x–t plane with three sets of characteristic lines (Fig. 4). When the initial condition is given, it can be obtained from the so-called Hugoniot relations [25–27] that across the front of a strong discontinuity the stress wave propagates in a material satisfies the following equation:
s ¼ rC vi v ¼ Eai e
ð18Þ
Furthermore, since the locus of the strong discontinuity wave should coincide with the rightward running characteristic line P(0,0) P(i,i), it can been found from Eq. (14) that dv ¼
1
rC vi
ds þ
ðEai þ Evi ÞsEai Evi e
rZvi C 2vi
dx
ð19Þ
Using Eq. (18), the ordinary differential equation of the stress is obtained:
ð9Þ
Eqs. (7)–(9) comprise the governing equations. In order to solve the wave propagation problem in the viscoelastic medium, a characteristic line method is adopted by multiplying Eqs. (7)–(9) by undetermined coefficients N, M and L, respectively [25]. The three equations are then added together, it gives @ ðN Zv Eai LÞ @t@ þ0 @x e þ Mr @t@ þ L @x@ v @ @ þ NZvi M s þ NðEai þ Evi ÞsNEai Evi e ¼ 0 ð10Þ @t @x Along a characteristic line C(x,t), the undetermined coefficients L, M and N satisfy the following relations: dx 0 L M ¼ ¼ ¼ ð11Þ dt Cðx,tÞ L þ NZvi Eai Mr NZvi
ds Eai s þ ¼0 dx 2Zvi C vi
ð20Þ
Thus, along the characteristic line P(0,0) P(i,i), the stress is determined by E s ¼ sð0,tÞexp ai x ð21Þ 2Zvi C vi where s(0,t) denotes the stress at boundary x ¼0. By a similar derivation, the velocity and strain along the characteristic lines P(0,0) P(i,i) can be obtained as Eai v ¼ vð0,tÞexp x ð22Þ 2Zvi C vi t
Two solutions are obtained from Eq. (11). The first solution is L þN Zvi Eai ¼ 0
rM2 ¼ LNZvi
ð12Þ P
Based on Eq. (12), two sets of characteristic lines are determined to be sffiffiffiffiffiffi dx Eai ¼7 ð13Þ ¼ 7 C vi dt r
P
P P
P P
P
ð14Þ
P
P
P
where Cvi is defined as the wave velocity along the characteristic lines. And the corresponding two sets of compatibility conditions along the characteristic lines are 1 ðE þEvi ÞsEai Evi e dv ¼ 7 ds 7 ai dt rC vi rZvi C vi 1 ðE þ Evi ÞsEai Evi e ds þ ai dx ¼7 rC vi rZvi C vi 2
P P
Undisturbed Zone x
Fig. 4. Characteristic lines in x–t plane for wave propagation across a viscoelastic medium.
14
L.F. Fan et al. / International Journal of Rock Mechanics & Mining Sciences 50 (2012) 11–18
and e ¼ eð0,tÞexp
Eai x 2Zvi C vi
ð23Þ
where v(0,t) and e(0,t) are the velocity and strain at boundary x¼ 0, respectively. Based on the above derivation, it is found that there are two cases to determine the mechanical behavior of all the points on the plane x t. One is the triangle case, such as P(0,j)P(0,j 2)P(1,j 1) in Fig. 4, and the other is the diamond case, such as P(i,j)P(i 1,j 1) P(i,j 2)P(iþ1,j 1) in Fig. 4. For the triangle case, the velocity v(0,j) on the boundary line P(0,j) is derived from the incident wave. Furthermore, along the leftward running characteristic line P(1,j 1)P(0,j), if the time interval Dt is small enough, it can be determined from the differential part of Eq. (14) that
sð0,jÞ ¼ sð1,j1ÞrC vi ½vð0,jÞvð1,j1Þ ðE þEvi Þsð1,j1ÞEai Evi eð1,j1Þ Dt ai Zvi
sð0,jÞsð0,j2Þ eð0,jÞ ¼ eð0,j2Þ þ rC 2vi ðE þ Evi Þsð0,j2ÞEai Evi eð0,j2Þ þ ai 2Dt rZvi C 2vi
þ
1
1
rC vi
ð26Þ
½sði,jÞsði þ1,j1Þ
ðEai þ Evi Þsði þ 1,j1ÞEai Evi eði þ1,j1Þ Dt rZvi C vi ð27Þ
and
sði,jÞsði,j2Þ rC 2vi ðEai þ Evi Þsði,j2ÞEai Evi eði,j2Þ þ 2Dt rZvi C 2vi
eði,jÞeði,j2Þ ¼
ð28Þ
, respectively. Eqs. (26)–(28) can be rewritten in an explicit form, i.e., 1 2
sði,jÞ ¼ rC vi ½vði þ 1,j1Þvði1,j1Þ 1 ½sðiþ 1,j1Þ þ sði1,j1Þ 2 1 ðEai þ Evi Þ ½sði þ 1,j1Þ þ sði1,j1ÞDt 2 Zvi þ
þ
vði,jÞ ¼
1 Eai Evi ½eði þ 1,j1Þ þ eði1,j1ÞDt 2 Zvi
1 1 ½sðiþ 1,j1Þsði1,j1Þ 2 rC vi
Fig. 1 illustrates the REV of a jointed rock mass, in which S is the initial length of the representative element. When a single joint was considered, the transmitted wave of a sinusoid wave across a linearly deformable joint was derived as [1,23] vT ðx,tÞ ¼
The differential part of the compatibility conditions of the leftward running wave of Eq. (14) and the compatibility conditions of the upward running wave of Eq. (17) are expressed as vði,jÞvði þ 1,j1Þ ¼
ð31Þ
Eqs. (24), (25) and (29)–(31) show that the response at time j can be calculated from those at times j 1 and j 2. Therefore, considering the boundary condition Eqs. (21)–(23), and the initial condition v(0,j), any particle velocity, stress and strain at the x t plane can be determined by an iterative computation using the differential Eqs. (24), (25) and (29)–(31).
ð25Þ
½sði,jÞsði1,j1Þ
ðEai þ Evi Þsði1,j1ÞEai Evi eði1,j1Þ Dt rZvi C vi
eði,jÞ ¼ eði,j2Þ þ
2.3. Parameter determination
For the diamond case, such as at P(i,j)P(i 1,j 1)P(i,j 2) P(iþ1,j 1), when Dt is very small, it is derived from the differential part of the rightward running wave compatibility conditions of Eq. (14) that
rC vi
sði,jÞsði,j2Þ rC 2vi ðE þ Evi Þsði,j2ÞEai Evi eði,j2Þ þ ai 2Dt rZvi C 2vi
ð30Þ
ð24Þ
For the upward running characteristic line P(0,j 2)P(0,j), the differential part of Eq. (17) is expressed to be
vði,jÞvði1,j1Þ ¼
1 ½vðiþ 1,j1Þ þ vði1,j1Þ 2 1 Eai Evi ½eði þ 1,j1Þeði1,j1ÞDt þ 2 rZvi C vi 1 Eai þ Evi ½sðiþ 1,j1Þsði1,j1ÞDt 2 rZvi C vi þ
2ki =z Aexp iðotxo=C 0 Þ io þ2ki =z
ð32Þ
where ki is defined in Eq. (1) as the normal joint stiffness, z is the wave impendence and z¼ rC0, C0 is the longitudinal wave velocity in an intact rock and o denotes the angular frequency of the sinusoid incident wave. When x ¼S, the phase shift coefficient ar and the attenuation coefficient br in Eq. (32) are, respectively, given in Eq. (33) which has also been discussed in [1]
8 o 1 o > > < ar ¼ C 0 þ S arctan 2ki =z
ð33Þ 2ki =z 1 > p ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > b ¼ ln : r S 2 o2 þ ð2ki =zÞ
When the same incident wave is applied to the REV of a viscoelastic equivalent medium, the phase shift coefficient ae and the attenuation coefficient be can be obtained from [26]: 8 ( " #)1=2 1=2 > > E2ai þ E2ci o2 t2i Eai þ Eci o2 t2i ro2 > > a ¼ þ > e 2Eci Eai > 1 þ o2 t2i 1 þ o2 t2i < ð34Þ ( " #)1=2 1=2 > > E2ai þ E2ci o2 t2i Eai þ Eci o2 t2i > ro2 > 1 þ o2 t2 > be ¼ 2Eci Eai > 1 þ o2 t2i : i Define Eci ¼EaiEvi/(Eai þEvi) and ti ¼ Zvi/Evi as the time of retardation of the nonlinear Voiget element. When x¼ S, the transmitted wave is obtained from Eq. (33) for the jointed rock mass or from Eq. (34) for the corresponding equivalent medium. Equating Eqs. (33) and (34), the parameters of the nonlinear viscoelastic model are determined by 88 2 391=2 !1=2 > 2 2 < ro2 2 2 > Eai þEci o2 t2i 5= > o 1 o > 4 Eai þ Eci o ti þ ¼ þ arctan > > 2 2 2 2 C 2k =z S ; < :2Eci Eai 1 þ o ti 1 þ o ti 0 i 8 2 2 391=2 3 !1=2 > > < ro2 > E2ai þ E2ci o2 t2i Eai þ Eci o2 t2i 5= 1 6 2ki =z > 7 4 > ¼ ln4qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi5 > :2E E : S ; 1 þ o2 t2i 1 þ o2 t2i ci ai o2 þð2ki =zÞ2
ð29Þ
ð35Þ It is recommended in [1], that the Young’s modulus of the intact rock can be used to give Ea in the viscoelastic model. Therefore, the relations between Evi, Zvi and o need be
L.F. Fan et al. / International Journal of Rock Mechanics & Mining Sciences 50 (2012) 11–18
determined by the implicit expressions given in Eq. (35) which can be solved by the Newton’s method [28]. If ki is kept constant, Evi and Zvi degrade to those of the linear model as proposed in [1]. In the following, it is assumed that the rock density r is 2650 kg/m3, and wave velocity C0 in the intact rock material is 5830 m/s. The closure variation for each segment is equally set to be Dd ¼0.005 mm, the initial stiffness of the joint in Eq. (1) is kn ¼3.5 GPa/m and the maximum closure of the joint is dmax ¼1 mm. A non-dimensional joint spacing G is defined as
15
the ratio of the joint spacing S to the wavelength of the incident wave l. Evi and Zvi relies on the incident wave frequency, the segment number and the joint spacing. Fig. 5 shows such relations when G ¼1, where f denotes the frequency of the incident wave, and f ¼ o=2p. It is determined from Fig. 5 that with a fixed segment number, Evi and Zvi increases as the increase of the frequency of the incident wave. And for a fixed frequency of the incident wave, Evi and Zvi increases as the segment number increases.
3. Verification of results A displacement discontinuity method is used to verify the proposed nonlinear viscoelastic equivalent medium method. The transmitted velocity passing a single joint is derived to be vT ðx1 ,t j þ 1 Þ ¼
2vi ð0,tx1 =C 0 Þ2vt ðx1 ,t j Þ Dt z=ki þ vT ðx1 ,t j Þ, si r s o si þ 1 ,
i ¼ 1,2,3
ð36Þ
where vI and vT denote the incident and transmitted waves, respectively. x1 is the location of the joint and s is the normal stress applied onto the joint, which can be obtained as s ¼ rC0vI [11] and ki is the piecewise linear stiffness of the joint as defined in Eq. (1). It is different from the linear displacement method that ki increases (decreases) as the magnitude of stress wave increases (decreases). It is also different from the results obtained by Zhao and Cai [11] because ki varies due to the assumption of the piecewise linear model in the present study. Eq. (36) shows that the response at time step jþ 1 can be calculated from those at time step j. Therefore, considering the boundary and initial conditions, any particle velocity just passing through the piecewise linear deformational joint can be calculated by an iterative computation using Eq. (36). To verify the present nonlinear viscoelastic equivalent medium method, an incident wave propagates through the RVE of a jointed rock mass and the corresponding nonlinear viscoelastic equivalent medium is studied, and the transmitted waves are compared as shown in Fig. 6. Assume the incident wave applied on boundary a in Fig. 1 is in the form of ( vI ð0,tÞ ¼
Fig. 5. Relations of parameters Evi, Zvi and frequency f for segment i (G ¼ 1). (a) Evi f and (b) Zvi f.
Asinð2pf 0 tÞ,
when 0 rt r 1=ð2f 0 Þ
0,
when t 4 1=ð2f 0 Þ
ð37Þ
where A is the amplitude and f0 the frequency of the incident wave, which are as shown in Fig. 6a–d with f0 ¼50 Hz and amplitude A¼0.057, 0.151, 0.340 and 0.529 m/s, respectively. The transmitted waves obtained based on the displacement discontinuity method, the linear viscoelastic equivalent medium method and the present nonlinear viscoelastic equivalent medium method with respect to different incident wave amplitudes are shown in Fig. 6a–d, respectively. ‘‘LVEMM’’ and ‘‘NVEMM’’ denote the linear viscoelastic equivalent medium method and the present nonlinear viscoelastic equivalent medium method, respectively. Fig. 6 shows that the transmitted waves obtained using DDM and NVEMM agree very well, which proves that the present NVEMM can effectively describe the dynamic property of the jointed rock mass. It is also obtained from Fig. 6 that for the incident wave with a small amplitude, the transmitted wave determined by the NVEMM approximates the transmitted wave by the LVEMM. However, the difference of the transmitted waves obtained by these two methods increases as the incident wave amplitude increases.
16
L.F. Fan et al. / International Journal of Rock Mechanics & Mining Sciences 50 (2012) 11–18
0.06
0.16 Incident wave Transmitted wave using DDM Transmitted wave using NVEMM Transmitted wave using LVEMM
0.04 0.03 0.02
0.12
0.01
0.10 0.08 0.06 0.04 0.02
0.00 0.00
Incident wave Transmitted wave using DDM Transmitted wave using NVEMM Transmitted wave using LVEMM
0.14
Velocity (m/s)
Velocity (m/s)
0.05
0.00 0.01
0.02
0.03
0.00
0.04
0.01
Incident wave Transmitted wave using DDM Transmitted wave using NVEMM Transmitted wave using LVEMM
0.35
0.03
0.04
0.25 0.20 0.15
Incident wave Transmitted wave using DDM Transmitted wave using NVEMM Transmitted wave using LVEMM
0.5 Velocity (m/s)
Velocity (m/s)
0.30
0.02
Time (s)
Time (s)
0.4 0.3 0.2
0.10 0.1
0.05
0.0
0.00 0.00
0.01
0.02
0.03
0.04
0.00
0.01
Time (s)
0.02
0.03
0.04
Time (s)
Fig. 6. Comparison of transmitted waves based on DDM, NVEMM and LVEMM with different amplitude incident waves. (a) A ¼ 0.057 m/s, (b) A ¼ 0.151 m/s, (c) A ¼ 0.340 m/ s and (d) A ¼ 0.529 m/s. (a) Evi f and (b) Zvi f.
1.0
4.1. Wave attenuation
0.9
The transmission coefficient is used to assess the attenuation of the transmitted wave. When an incident wave given in the form of Eq. (37) propagates across a jointed rock mass with a joint spacing of S¼5.83 m, the relations of the transmission coefficient versus the incident wave amplitude and the incident wave frequency are investigated, respectively. The relation between the transmission coefficient and the amplitude of the incident wave based on the LVEMM and the present NVEMM for frequencies of f0 ¼100, 200, 300 and 400 Hz are shown in Fig. 7. It can be seen from Fig. 7 that the transmission coefficient obtained using the LVEMM is independent of the incident wave amplitude, while those obtained by the NVEMM is amplitude dependent. This is reasonable because a nonlinear model is adopted in the NVEMM. For the incident wave with a certain frequency, the transmission coefficient increases with the increase of the incident wave amplitude based on the present NVEMM. All the transmission coefficients converge to 1 if the amplitude of the incident wave is sufficiently large. This phenomenon can be explained that as the incident wave amplitude increases, the joint is compacted nonlinearly and thus the joint stiffness increases accordingly. The stiffness of the joint approaches the stiffness of the rock material when the joint is fully compacted. From Fig. 7, the transmission coefficient for the incident wave with a larger frequency is smaller than that with a smaller frequency, which coincides with
Transmission coefficient T
4. Discussion
LVEMM
NVEMM f =100
0.8
f =200 f =300
0.7
f =400
0.6 0.5 0.4 0.3 0.2 0.0
0.2 0.4 0.6 0.8 Amplitude of incident wave A (m/s)
1.0
Fig. 7. Effect of incident wave amplitude on transmission coefficient.
the geological investigation that a joint filters out stress waves with higher frequencies, but has less effect on the stress waves with lower frequencies [12]. Fig. 8 illustrates the relation between the transmission coefficient and the incident wave frequency with respect to different incident wave amplitudes, i.e., A¼0.05 m, 0.1 m, 0.3 m, 0.5 m and 0.8 m, respectively. It is seen from Fig. 8 that there is a flat low frequency region within which the transmission coefficient is nearly a constant or reduces slightly with frequency increasing.
L.F. Fan et al. / International Journal of Rock Mechanics & Mining Sciences 50 (2012) 11–18
17
6000
1.0
0.6 0.4
Effective Velocity (m/s)
Transmission Coefficient T
5500 NVEMM A=0.80 NVEMM A=0.50 NVEMM A=0.30 NVEMM A=0.10 NVEMM A=0.05 LVEMM
0.8
0.2
5000
NVEMM NVEMM NVEMM NVEMM NVEMM LVEMM
4500 4000 3500 3000 2500
0.0
2000 1500
-0.2 101
102
103
104
105
106
107
101
Frequency of incident wave f0 (Hz) Fig. 8. Effect of incident wave frequency on transmission coefficient.
5500 5000 Effective velocity (m/s)
A=0.80 A=0.50 A=0.30 A=0.10 A=0.05
LVEMM
NVEMM f =100
4500
f =200 f =300 f =400
4000 3500 3000 2500 0.0
0.2
0.4
0.6
0.8
1.0
Amplitude of incident wave (m) Fig. 9. Effect of incident wave amplitude on effective wave propagation velocity.
When the frequency is sufficiently large, the transmission coefficient becomes zero. Between these two extremes, there is a rapid decreasing region of the transmission coefficient. For a fixed incident wave amplitude, the incident waves with a larger frequency result in less transmission. For the same frequency, the incident waves of the larger amplitudes result in a larger transmission coefficient. It is also found from Fig. 8 that as the amplitude of the incident wave decreases, the relationship between the transmission coefficient and the incident wave frequency obtained by the NVEMM approaches to that given by the LVEMM. 4.2. Effective velocity The effective wave propagation velocity in an equivalent medium is defined as the wave traveling length over the corresponding wave traveling time. The wave traveling time is measured as the time difference between the peak velocity of the incident and transmitted waves. Considering the joint spacing S ¼5.83 m, the effective velocity versus the incident wave amplitude for different frequencies of f0 ¼100 Hz, 200 Hz, 300 Hz and 400 Hz are shown in Fig. 9. It is seen that, with a fixed wave frequency, the larger the incident wave amplitude, the larger the effective velocity. The effective velocity approaches the maximum
102 103 104 Frequency of incident wave f0 (Hz)
105
Fig. 10. Effect of incident wave frequency on effective wave propagation velocity.
value when the stiffness of the joint is close to the stiffness of the intact rock. In the LVEMM, however, the effective velocity is independent to the incident wave amplitude. The relations between the effective velocity and the incident wave frequency for different incident wave amplitudes, i.e., A¼0.05 m, 0.1 m, 0.3 m, 0.5 m and 0.8 m are shown in Fig. 10. There is a low frequency region for which the effective velocity is the minimum and it varies slightly in the region. And the effective velocity approaches the wave propagation velocity of the intact rock material at high frequencies, which is v¼5830 m/s in the present analysis. Between these two regions, there is a rapid increase in the effective velocity as the frequency increases. Similar results have been obtained by the displacement discontinuity method [23]. It is also observed from Fig. 10 that as the incident wave amplitude decreases, the relationship between the effective velocity and the incident wave frequency obtained by the NVEMM approaches to that by the LVEMM. The LVEMM can be used to replace the NVEMM approximately if the amplitude of the incident wave is insufficient to cause nonlinear deformation of the joint.
5. Conclusions A nonlinear viscoelastic equivalent medium model is proposed in the present study for the investigation of stress wave propagation in a rock mass consisting equally-spaced parallel joints with nonlinearly deformational behavior. The nonlinear behavior of the rock joint is approximated by a piecewise linear joint model. A method of characteristic line is adopted to solve the governing equations of stress wave propagation in the viscoelastic equivalent continuum. The parameters of the present viscoelastic equivalent model are theoretically derived and the frequencydependence and the amplitude-dependence of the transmission coefficient and effective wave propagation velocity are presented and discussed in detail. It is concluded that stress wave propagation through a nonlinear rock joint is amplitude-dependent. The conventional LVEMM is limited by the linear deformation assumption, and it is not able to describe the dependence of stress wave propagation with the incident wave amplitude. The present NVEMM not only determines accurately the transmission coefficient and the effective velocity when the amplitude of incident stress wave is sufficiently large to mobilize the nonlinear deformation of the
18
L.F. Fan et al. / International Journal of Rock Mechanics & Mining Sciences 50 (2012) 11–18
rock joint, but also degrades to exactly the LVEMM when the amplitude of the incident stress wave is small. Thus, the conventional LVEMM can be regarded as a special case of the present NVEMM. The proposed NVEMM has advantage over the DDM. The comparison of the results obtained by the DDM and the present NVEMM shows that the NVEMM can be used to study stress wave propagation in jointed rock masses as accurate as the DDM. However, because of its continuous nature in an equivalent medium model, the proposed NVEMM becomes more efficient in analyzing stress wave propagation in jointed rock masses when the joint number increases. It should be mentioned that the present study gives only some preliminary results when a jointed rock mass is approximated by a nonlinear viscoelastic equivalent medium. Further researches are still necessary to consider more complicated cases, for example, inclined impingement of a stress wave onto a rock joint, a rock mass with multiple joint sets, wave multiple-reflection effects, shear wave propagation, etc.
Acknowledgments The authors would like to acknowledge the support of the Collaborative Research Fund with Overseas, Hong Kong and Macau Scholars (Grant No 51028801) and the general project (Grant No 51178012) of the National Natural Science Foundation of China.
References [1] Li JC, Ma GW, Zhao J. An equivalent viscoelastic model for rock mass with parallel joints. J Geophys Res 2010;115:B03305. [2] Bandis SC, Lumsden AC, Barton NR. Fundamentals of rock joint deformation. Int J Rock Mech Min Sci Geomech Abstr 1983;20(6):249–68. [3] Cook NGW. Natural joint in rock: mechanical, hydraulic and seismic behaviour and properties under normal stress. Int J Rock Mech Min Sci Geomech Abstr 1992;29(3):198–223. [4] Ma GW, Li JC, Zhao J. Three-phase medium model for filled rock joint and interaction with stress waves. Int J Numer Anal Methods Geomech 2011;35(1): 97–110. [5] Li JC, Ma GW. Experimental study of stress wave propagation across a filled rock joint. Int J Rock Mech Min Sci 2009;46:471–8.
[6] Zhao J, Cai JG, Zhao XB, Li HB. Dynamic model of fracture normal behaviour and application to prediction of stress wave attenuation across fractures. Rock Mech Rock Eng 2008;41(5):671–93. [7] Perino A, Zhu JB, Li JC, Barla G, Zhao J. Theoretical methods for wave propagation across jointed rock masses. Rock Mech Rock Eng 2010;43:799–809. [8] Zhao J, Zhao XB, Cai JG. A further study of P-wave attenuation across parallel fractures with linear deformational behavior. Int J Rock Mech Min Sci 2006;43:776–88. [9] Miller RK. An approximate method of analysis of the transmission of elastic waves through a frictional boundary. J Appl Mech 1977;44:652–6. [10] Schoenberg M. Elastic wave behaviour across linear slip interfaces. J Acoust Soc Am 1980;68(5):1516–21. [11] Zhao J, Cai JG. Transmission of elastic P-waves across single fractures with a nonlinear normal deformational behavior. Rock Mech Rock Eng 2001;34(1): 3–22. [12] Pyrak-Nolte LJ Seismic visibility of fractures. PhD thesis. Berkeley, CA, USA: University of California; 1988. [13] Li JC, Ma GW, Zhao J. Analysis of stochastic seismic wave interaction with a slippery rock fault. Rock Mech Rock Eng 2011;44:85–92. [14] White JE. Underground sound. New York: Elsevier; 1983. [15] Schoenberg M, Sayers CM. Seismic anisotropy of fractured rock. Geophysics 1995;60(1):204–11. [16] Hudson JA. Wave speeds and attenuation of elastic waves in material containing cracks. Geophys J R Astron Soc 1981;64:133–50. [17] Amadei B, Goodman RE. A 3-D constitutive relation for fractured rock masses. In: Proceedings of the International Symposium on the Mechanical Behavior of Structural Media, Balkema; p. 249–68; 1981. [18] Schoenberg M. Reflection of elastic waves from periodically stratified media with interfacial slip. Geophys Prospect 1983;31:265–92. [19] Majer EL, McEvilly TV, Eastwood FS, Myer LR. Fracture detection using P- and S-wave vertical seismic profiling at The Geysers. Geophysics 1988;53(1): 76–84. [20] Sitharam TG. Equivalent continuum analyses of jointed rock mass: Some case studies. in: Ramamurthy T, editor. Engineering in Rocks for slopes, foundations and tunnels. New Delhi, India: Prentice-Hall of India; 2007. p. 518–44. [21] Zhao XB, Zhao J, Cai JG. P-wave transmission across fractures with nonlinear deformational behavior. Int J Numer Anal Methods Geomech 2006;30(11): 1097–112. [22] Li JC, Ma GW, Huang X. Analysis of wave propagation through a filled rock joint. Rock Mech Rock Eng 2010;43:789–98. [23] Pyrak-Nolte LJ, Myer LR, Cook NGW. Anisotropy in seismic velocities and amplitudes from multiple parallel fractures. J Geophys Res 1990;95(B7): 11345–58. [24] Barton NR, Bandis SC, Bakhtar K. Strength deformation and conductivity coupling of rock joints. Int J Rock Mech Min Sci Geomech Abstr 1985;22(3): 121–40. [25] Wang LL, Labibes K, Azari Z, Pluvinage G. Generalization of split Hopkinson bar technique to use viscoelastic bars. Int J Impact Eng 1994;15(5):669–86. [26] Kolsky H. Stress waves in solids. Oxford: Clarendon Press; 1953. [27] Wang LL. Foundations of stress waves. Beijing: National Defense Industry Press; 1985. [28] Kelley CT. Solving nonlinear equations with Newton’s method. Philadelphia: SIAM; 2003.
Contents lists available at SciVerse ScienceDirect
International Journal of Rock Mechanics & Mining Sciences journal homepage: www.elsevier.com/locate/ijrmms
Nonlinear viscoelastic medium equivalence for stress wave propagation in a jointed rock mass L.F. Fan a, G.W. Ma b,n, J.C. Li c a
School of Civil and Environmental Engineering, Nanyang Technological University, Singapore 639798, Singapore School of Civil and Resource Engineering, The University of Western Australia, Crawley, WA 6009, Australia c Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, Wuhan 430071, China b
a r t i c l e i n f o
abstract
Article history: Received 16 August 2010 Received in revised form 26 November 2011 Accepted 17 December 2011 Available online 30 January 2012
A nonlinear viscoelastic equivalent medium model is proposed for longitudinal stress wave propagation in a rock mass that contains equally spaced parallel joints. The joint deformation is assumed to satisfy the nonlinear Bandis–Barton (B–B) deformational model. By the simplification of the nonlinear B–B model into piecewise linear segments, the material parameters in the equivalent medium model are analytically derived. The nonlinear viscoelastic equivalent medium model is verified by comparison with the results based on the displacement discontinuity method. It is shown that for smaller amplitude stress waves, when the rock joints deform in the linear elastic region, the nonlinear viscoelastic equivalent medium model is reduced to its linear form, whereas for larger amplitude stress waves, the present nonlinear model describes the wave amplitude-dependence of wave velocity. Two important properties, i.e., the wave transmission coefficient and the effective wave propagation velocity, are discussed in detail. & 2011 Elsevier Ltd. All rights reserved.
Keywords: Equivalent medium method Nonlinear deformational joint Jointed rock mass Viscoelastic medium Wave propagation
1. Introduction The presence of rock joints often dominates the dynamic engineering properties and increases the complexity in modeling dynamic behavior of a jointed rock mass [1,2]. It is of interest to develop an efficient and accurate analytical model to describe the dynamic behavior of jointed rock masses and evaluate the effects of joints on stress wave propagation. It has been observed that the complete deformational behavior of a rock joint is generally nonlinear [2–6]. The linear elastic joint model can be regarded as a special case of a nonlinear model such as the B–B model [2] when the joint is dry and the magnitude of the stress wave is too small to mobilize the joint nonlinear deformation. Currently, there are two different theoretical approaches to evaluate the effect of rock joints on the stress wave propagation through a rock mass [1,7,8]. One is the displacement discontinuity method (DDM) [3,8–13] and the other is the equivalent medium method (EMM) [1,7,14–20]. In the DDM, the joints are assumed to be non-welded interfaces. The displacements across the joint interface are discontinuous, while the stresses across the joint interface are continuous. These two assumptions are used as two boundary conditions in the stress wave propagation equations. By using the DDM, a series of analyses of stress wave propagation
n
Corresponding author. E-mail address: [email protected] (G.W. Ma).
1365-1609/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijrmms.2011.12.008
through rock joints have been carried out [11,21,22], and it has been applied successfully to predict the seismic responses of jointed rock masses [3,11–13]. Since the DDM treats the rock mass as a discrete rock medium, it is difficult to give an explicit expression of stress wave propagation in a jointed rock mass, and it becomes very tedious and complex when the number of discontinuities increases [1]. On the other hand, the equivalent medium method treats the rock mass as a whole and predicts the effects of joint on the behavior of a rock mass system using a representative elementary volume (REV), so as to make a continuum medium analysis of practical problems. An elastic continuum medium is commonly used to represent the jointed rock mass [16–19]. Such conventional linear elastic equivalence is effective only if the stress wave frequency and amplitude effects, and the multiple wave reflections among the joints are ignored [1,8,23]. Li et al. [1] introduced a viscoelastic equivalent medium model for rock masses consisting of equally-spaced parallel joints. A virtual wave source (VWS) concept has been put forward to simulate the multiple wave reflections among the joints [1,7]. The wave dispersion, attenuation and frequency dependence of a jointed rock mass have been analytically investigated by solving an explicit stress wave propagation equation. However, the joints were assumed to deform linearly in the previous study. A nonlinear equivalent medium method is necessary in order to efficiently analyze stress wave propagation in a jointed rock mass where the rock joint deforms nonlinearly.
12
L.F. Fan et al. / International Journal of Rock Mechanics & Mining Sciences 50 (2012) 11–18
This paper proposes a nonlinear viscoelastic equivalent medium model for the study of stress wave propagation through a rock mass with equally spaced parallel rock joints. The deformation of the rock joints satisfies the B–B model [2]. The proposed equivalent medium model is an extended standard solid model with three frequency dependent piecewise-linear parameter sets. In order to derive the wave propagation equation in the viscoelastic medium efficiently, a method of characteristic lines is adopted. To theoretically verify the present model, a comparison with the results obtained from the displacement discontinuity method is performed. The amplitude dependence, which only exists in stress wave propagation through a nonlinear medium, and the frequency dependence of the present viscoelastic model are discussed.
2. Nonlinear viscoelastic medium equivalence
d¼
Fig. 1 illustrates an REV of a jointed rock mass and the definition of the corresponding equivalent medium. In the present study, the joint is assumed to be planar, large in extent and small in thickness. The length of the REV of the equivalent medium is defined to be the joint spacing S as shown in Fig. 1. The normal deformational behavior of a rock joint under quasi-static and cyclic loadings has been studied [2,5]. Fig. 2 illustrates a piecewise linear model to consider the nonlinear deformational behavior of a rock joint, which is compared with the linear joint model and the B–B model [2]. If the closure (opening) of the joint and the compression (tension) stress are assumed to be positive (negative), the normal constitutive relations between the stress and the closure can be described by the
s
ð1Þ
kn þ s=dmax
where d is the closure of the joint, s is the normal stress, dmax is the maximum allowable closure of the joint and kn is the normal stiffness of the joint at the initial state. The dashed curve in Fig. 2 gives a typical hyperbolic fitting of the B–B model [11]. A piecewise linear model is proposed to approximate such hyperbolic B–B model. As shown in Fig. 2 that the curve of the B–B model is a continuous hyperbolic dashed line, while the present piecewise linear model can be regarded as a polygonal curve consisting of a series of linear segments with their end points located on the hyperbolic B–B model. Therefore, the values of the closure and the stress at the two ends of each segment are assumed to satisfy the hyperbolic B–B model relation yields di ¼
2.1. Joint models
si kn þ si =dmax
,
i ¼ 1,2,3,
ð2Þ
where i is the segment number, di and si denote the discrete threshold closure and stress at the two ends of each segment, respectively. The difference of the stress and closure for every individual segment can be obtained as
Dsi ¼ si si1 , i ¼ 1,2,3,
ð3Þ
and
Ddi ¼ di di1 , i ¼ 1,2,3,
ð4Þ
, respectively. Therefore, the stiffness of the joint for segment i can be obtained from Eqs. (2)–(4) as Dsi kn di di1 , i ¼ 1,2,3, ki ¼ ¼ Ddi di di1 1ðdi =dmax Þ 1ðdi1 =dmax Þ ð5Þ Considering Eqs. (2) and (5), the stress–closure relation of the present piecewise linear model is described as
Joint Rock
d ¼ di1 þ
a
b
a
b
S
S
Fig. 1. Representative elementary volume. (a) Jointed rock mass; (b) Equivalent medium.
14
Piecewise linear model BB model Linear model
12 Normal stress σ (MPa)
B–B model as,
d
10
(di,σi)
8 ki
6
(di-1,σi-1)
Δσ
ssi1 ki
,
di1 rd o di
ð6Þ
where di is obtained from Eq. (2) and ki is obtained from Eq. (5). It can be obtained from Fig. 2 that for given stress s, the B–B model reveals that the closure will increase following the locus of the dashed hyperbolic curve to d, while the present piecewise model reveals that the closure will increase following the locus of the solid linear segments to d. When the closure difference Ddi is small enough, the present piecewise linear model approximates the B–B model with sufficient accuracy. The values of dmax and kn for the piecewise linear model and the B–B model can be determined from the laboratory measurements as obtained by Bandis et al. [2] and Barton et al. [24]. It can also be realized from Fig. 2 that the present piecewise linear model and the B–B model are able to describe the joint stiffness increase with the increase of the normal stress, however the joint stiffness of a linear model always keeps a constant. 2.2. Nonlinear viscoelastic equivalent medium method
Δd
4 2 0 0.0
0.1
0.2
0.3
0.4 0.5 0.6 Closure d (mm)
0.7
0.8
0.9
Fig. 2. Piecewise linear model, B–B model and linear model of a rock joint.
A nonlinear viscoelastic equivalent model is proposed to describe the overall dynamic mechanical behavior of a rock mass containing nonlinear deformational joints, which is shown in Fig. 3. A piecewise linear spring with a modulus of Evi is in parallel with a piecewise linear dashpot Zvi, then in series with another piecewise linear spring Eai. The piecewise linear quantities of Eai, Evi and Zvi are related to the stress–closure relation of the joint at segment i. For the ith segment, Eai, Evi and Zvi can be obtained from the coefficients
L.F. Fan et al. / International Journal of Rock Mechanics & Mining Sciences 50 (2012) 11–18
13
where the sign ‘‘ þ’’ means the rightward running waves and ‘‘ ’’ denotes the leftward running waves. The second solution of Eq. (11) is
Evi Eai
L¼M¼0 N a0
ð15Þ
Therefore, the characteristic line of the second solution is
ηvi
dx ¼ 0
ð16Þ
Fig. 3. Piecewise linear viscoelastic equivalent medium model.
of the phase shift and amplitude attenuation of stress wave propagation [1]. The constitutive equation of the nonlinear model shown in Fig. 3 is given as ðEai þ Evi Þs þ Zvi
@s @e Zvi Eai Eai Evi e ¼ 0 @t @t
ð7Þ
where e is the strain component, Eai and Evi are the piecewise elastic modulus of the springs, respectively, Zvi is the piecewise viscosity of the dashpot, t denotes the time variable, and Eai, Evi and Zvi vary with the change of the magnitude of joint stiffness ki. The motion equation of one-dimensional longitudinal wave propagation is governed by
r
@v @s ¼ @t @x
ð8Þ
where r is the density of the medium material, x is the Cartesian coordinate along the wave propagation path and v is the particle velocity at x. In addition to the constitutive Eq. (7) and the equation of motion (8), the continuity equation is known as @v @e ¼ @x @t
The compatibility condition along the characteristic line of Eq. (16) is de
ds ðE þEvi ÞsEai Evi e ai dt ¼ 0 rC 2vi rZvi C 2vi
ð17Þ
From Eqs. (13) and (16), the wave propagation in the viscoelastic medium can be described by the x–t plane with three sets of characteristic lines (Fig. 4). When the initial condition is given, it can be obtained from the so-called Hugoniot relations [25–27] that across the front of a strong discontinuity the stress wave propagates in a material satisfies the following equation:
s ¼ rC vi v ¼ Eai e
ð18Þ
Furthermore, since the locus of the strong discontinuity wave should coincide with the rightward running characteristic line P(0,0) P(i,i), it can been found from Eq. (14) that dv ¼
1
rC vi
ds þ
ðEai þ Evi ÞsEai Evi e
rZvi C 2vi
dx
ð19Þ
Using Eq. (18), the ordinary differential equation of the stress is obtained:
ð9Þ
Eqs. (7)–(9) comprise the governing equations. In order to solve the wave propagation problem in the viscoelastic medium, a characteristic line method is adopted by multiplying Eqs. (7)–(9) by undetermined coefficients N, M and L, respectively [25]. The three equations are then added together, it gives @ ðN Zv Eai LÞ @t@ þ0 @x e þ Mr @t@ þ L @x@ v @ @ þ NZvi M s þ NðEai þ Evi ÞsNEai Evi e ¼ 0 ð10Þ @t @x Along a characteristic line C(x,t), the undetermined coefficients L, M and N satisfy the following relations: dx 0 L M ¼ ¼ ¼ ð11Þ dt Cðx,tÞ L þ NZvi Eai Mr NZvi
ds Eai s þ ¼0 dx 2Zvi C vi
ð20Þ
Thus, along the characteristic line P(0,0) P(i,i), the stress is determined by E s ¼ sð0,tÞexp ai x ð21Þ 2Zvi C vi where s(0,t) denotes the stress at boundary x ¼0. By a similar derivation, the velocity and strain along the characteristic lines P(0,0) P(i,i) can be obtained as Eai v ¼ vð0,tÞexp x ð22Þ 2Zvi C vi t
Two solutions are obtained from Eq. (11). The first solution is L þN Zvi Eai ¼ 0
rM2 ¼ LNZvi
ð12Þ P
Based on Eq. (12), two sets of characteristic lines are determined to be sffiffiffiffiffiffi dx Eai ¼7 ð13Þ ¼ 7 C vi dt r
P
P P
P P
P
ð14Þ
P
P
P
where Cvi is defined as the wave velocity along the characteristic lines. And the corresponding two sets of compatibility conditions along the characteristic lines are 1 ðE þEvi ÞsEai Evi e dv ¼ 7 ds 7 ai dt rC vi rZvi C vi 1 ðE þ Evi ÞsEai Evi e ds þ ai dx ¼7 rC vi rZvi C vi 2
P P
Undisturbed Zone x
Fig. 4. Characteristic lines in x–t plane for wave propagation across a viscoelastic medium.
14
L.F. Fan et al. / International Journal of Rock Mechanics & Mining Sciences 50 (2012) 11–18
and e ¼ eð0,tÞexp
Eai x 2Zvi C vi
ð23Þ
where v(0,t) and e(0,t) are the velocity and strain at boundary x¼ 0, respectively. Based on the above derivation, it is found that there are two cases to determine the mechanical behavior of all the points on the plane x t. One is the triangle case, such as P(0,j)P(0,j 2)P(1,j 1) in Fig. 4, and the other is the diamond case, such as P(i,j)P(i 1,j 1) P(i,j 2)P(iþ1,j 1) in Fig. 4. For the triangle case, the velocity v(0,j) on the boundary line P(0,j) is derived from the incident wave. Furthermore, along the leftward running characteristic line P(1,j 1)P(0,j), if the time interval Dt is small enough, it can be determined from the differential part of Eq. (14) that
sð0,jÞ ¼ sð1,j1ÞrC vi ½vð0,jÞvð1,j1Þ ðE þEvi Þsð1,j1ÞEai Evi eð1,j1Þ Dt ai Zvi
sð0,jÞsð0,j2Þ eð0,jÞ ¼ eð0,j2Þ þ rC 2vi ðE þ Evi Þsð0,j2ÞEai Evi eð0,j2Þ þ ai 2Dt rZvi C 2vi
þ
1
1
rC vi
ð26Þ
½sði,jÞsði þ1,j1Þ
ðEai þ Evi Þsði þ 1,j1ÞEai Evi eði þ1,j1Þ Dt rZvi C vi ð27Þ
and
sði,jÞsði,j2Þ rC 2vi ðEai þ Evi Þsði,j2ÞEai Evi eði,j2Þ þ 2Dt rZvi C 2vi
eði,jÞeði,j2Þ ¼
ð28Þ
, respectively. Eqs. (26)–(28) can be rewritten in an explicit form, i.e., 1 2
sði,jÞ ¼ rC vi ½vði þ 1,j1Þvði1,j1Þ 1 ½sðiþ 1,j1Þ þ sði1,j1Þ 2 1 ðEai þ Evi Þ ½sði þ 1,j1Þ þ sði1,j1ÞDt 2 Zvi þ
þ
vði,jÞ ¼
1 Eai Evi ½eði þ 1,j1Þ þ eði1,j1ÞDt 2 Zvi
1 1 ½sðiþ 1,j1Þsði1,j1Þ 2 rC vi
Fig. 1 illustrates the REV of a jointed rock mass, in which S is the initial length of the representative element. When a single joint was considered, the transmitted wave of a sinusoid wave across a linearly deformable joint was derived as [1,23] vT ðx,tÞ ¼
The differential part of the compatibility conditions of the leftward running wave of Eq. (14) and the compatibility conditions of the upward running wave of Eq. (17) are expressed as vði,jÞvði þ 1,j1Þ ¼
ð31Þ
Eqs. (24), (25) and (29)–(31) show that the response at time j can be calculated from those at times j 1 and j 2. Therefore, considering the boundary condition Eqs. (21)–(23), and the initial condition v(0,j), any particle velocity, stress and strain at the x t plane can be determined by an iterative computation using the differential Eqs. (24), (25) and (29)–(31).
ð25Þ
½sði,jÞsði1,j1Þ
ðEai þ Evi Þsði1,j1ÞEai Evi eði1,j1Þ Dt rZvi C vi
eði,jÞ ¼ eði,j2Þ þ
2.3. Parameter determination
For the diamond case, such as at P(i,j)P(i 1,j 1)P(i,j 2) P(iþ1,j 1), when Dt is very small, it is derived from the differential part of the rightward running wave compatibility conditions of Eq. (14) that
rC vi
sði,jÞsði,j2Þ rC 2vi ðE þ Evi Þsði,j2ÞEai Evi eði,j2Þ þ ai 2Dt rZvi C 2vi
ð30Þ
ð24Þ
For the upward running characteristic line P(0,j 2)P(0,j), the differential part of Eq. (17) is expressed to be
vði,jÞvði1,j1Þ ¼
1 ½vðiþ 1,j1Þ þ vði1,j1Þ 2 1 Eai Evi ½eði þ 1,j1Þeði1,j1ÞDt þ 2 rZvi C vi 1 Eai þ Evi ½sðiþ 1,j1Þsði1,j1ÞDt 2 rZvi C vi þ
2ki =z Aexp iðotxo=C 0 Þ io þ2ki =z
ð32Þ
where ki is defined in Eq. (1) as the normal joint stiffness, z is the wave impendence and z¼ rC0, C0 is the longitudinal wave velocity in an intact rock and o denotes the angular frequency of the sinusoid incident wave. When x ¼S, the phase shift coefficient ar and the attenuation coefficient br in Eq. (32) are, respectively, given in Eq. (33) which has also been discussed in [1]
8 o 1 o > > < ar ¼ C 0 þ S arctan 2ki =z
ð33Þ 2ki =z 1 > p ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > b ¼ ln : r S 2 o2 þ ð2ki =zÞ
When the same incident wave is applied to the REV of a viscoelastic equivalent medium, the phase shift coefficient ae and the attenuation coefficient be can be obtained from [26]: 8 ( " #)1=2 1=2 > > E2ai þ E2ci o2 t2i Eai þ Eci o2 t2i ro2 > > a ¼ þ > e 2Eci Eai > 1 þ o2 t2i 1 þ o2 t2i < ð34Þ ( " #)1=2 1=2 > > E2ai þ E2ci o2 t2i Eai þ Eci o2 t2i > ro2 > 1 þ o2 t2 > be ¼ 2Eci Eai > 1 þ o2 t2i : i Define Eci ¼EaiEvi/(Eai þEvi) and ti ¼ Zvi/Evi as the time of retardation of the nonlinear Voiget element. When x¼ S, the transmitted wave is obtained from Eq. (33) for the jointed rock mass or from Eq. (34) for the corresponding equivalent medium. Equating Eqs. (33) and (34), the parameters of the nonlinear viscoelastic model are determined by 88 2 391=2 !1=2 > 2 2 < ro2 2 2 > Eai þEci o2 t2i 5= > o 1 o > 4 Eai þ Eci o ti þ ¼ þ arctan > > 2 2 2 2 C 2k =z S ; < :2Eci Eai 1 þ o ti 1 þ o ti 0 i 8 2 2 391=2 3 !1=2 > > < ro2 > E2ai þ E2ci o2 t2i Eai þ Eci o2 t2i 5= 1 6 2ki =z > 7 4 > ¼ ln4qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi5 > :2E E : S ; 1 þ o2 t2i 1 þ o2 t2i ci ai o2 þð2ki =zÞ2
ð29Þ
ð35Þ It is recommended in [1], that the Young’s modulus of the intact rock can be used to give Ea in the viscoelastic model. Therefore, the relations between Evi, Zvi and o need be
L.F. Fan et al. / International Journal of Rock Mechanics & Mining Sciences 50 (2012) 11–18
determined by the implicit expressions given in Eq. (35) which can be solved by the Newton’s method [28]. If ki is kept constant, Evi and Zvi degrade to those of the linear model as proposed in [1]. In the following, it is assumed that the rock density r is 2650 kg/m3, and wave velocity C0 in the intact rock material is 5830 m/s. The closure variation for each segment is equally set to be Dd ¼0.005 mm, the initial stiffness of the joint in Eq. (1) is kn ¼3.5 GPa/m and the maximum closure of the joint is dmax ¼1 mm. A non-dimensional joint spacing G is defined as
15
the ratio of the joint spacing S to the wavelength of the incident wave l. Evi and Zvi relies on the incident wave frequency, the segment number and the joint spacing. Fig. 5 shows such relations when G ¼1, where f denotes the frequency of the incident wave, and f ¼ o=2p. It is determined from Fig. 5 that with a fixed segment number, Evi and Zvi increases as the increase of the frequency of the incident wave. And for a fixed frequency of the incident wave, Evi and Zvi increases as the segment number increases.
3. Verification of results A displacement discontinuity method is used to verify the proposed nonlinear viscoelastic equivalent medium method. The transmitted velocity passing a single joint is derived to be vT ðx1 ,t j þ 1 Þ ¼
2vi ð0,tx1 =C 0 Þ2vt ðx1 ,t j Þ Dt z=ki þ vT ðx1 ,t j Þ, si r s o si þ 1 ,
i ¼ 1,2,3
ð36Þ
where vI and vT denote the incident and transmitted waves, respectively. x1 is the location of the joint and s is the normal stress applied onto the joint, which can be obtained as s ¼ rC0vI [11] and ki is the piecewise linear stiffness of the joint as defined in Eq. (1). It is different from the linear displacement method that ki increases (decreases) as the magnitude of stress wave increases (decreases). It is also different from the results obtained by Zhao and Cai [11] because ki varies due to the assumption of the piecewise linear model in the present study. Eq. (36) shows that the response at time step jþ 1 can be calculated from those at time step j. Therefore, considering the boundary and initial conditions, any particle velocity just passing through the piecewise linear deformational joint can be calculated by an iterative computation using Eq. (36). To verify the present nonlinear viscoelastic equivalent medium method, an incident wave propagates through the RVE of a jointed rock mass and the corresponding nonlinear viscoelastic equivalent medium is studied, and the transmitted waves are compared as shown in Fig. 6. Assume the incident wave applied on boundary a in Fig. 1 is in the form of ( vI ð0,tÞ ¼
Fig. 5. Relations of parameters Evi, Zvi and frequency f for segment i (G ¼ 1). (a) Evi f and (b) Zvi f.
Asinð2pf 0 tÞ,
when 0 rt r 1=ð2f 0 Þ
0,
when t 4 1=ð2f 0 Þ
ð37Þ
where A is the amplitude and f0 the frequency of the incident wave, which are as shown in Fig. 6a–d with f0 ¼50 Hz and amplitude A¼0.057, 0.151, 0.340 and 0.529 m/s, respectively. The transmitted waves obtained based on the displacement discontinuity method, the linear viscoelastic equivalent medium method and the present nonlinear viscoelastic equivalent medium method with respect to different incident wave amplitudes are shown in Fig. 6a–d, respectively. ‘‘LVEMM’’ and ‘‘NVEMM’’ denote the linear viscoelastic equivalent medium method and the present nonlinear viscoelastic equivalent medium method, respectively. Fig. 6 shows that the transmitted waves obtained using DDM and NVEMM agree very well, which proves that the present NVEMM can effectively describe the dynamic property of the jointed rock mass. It is also obtained from Fig. 6 that for the incident wave with a small amplitude, the transmitted wave determined by the NVEMM approximates the transmitted wave by the LVEMM. However, the difference of the transmitted waves obtained by these two methods increases as the incident wave amplitude increases.
16
L.F. Fan et al. / International Journal of Rock Mechanics & Mining Sciences 50 (2012) 11–18
0.06
0.16 Incident wave Transmitted wave using DDM Transmitted wave using NVEMM Transmitted wave using LVEMM
0.04 0.03 0.02
0.12
0.01
0.10 0.08 0.06 0.04 0.02
0.00 0.00
Incident wave Transmitted wave using DDM Transmitted wave using NVEMM Transmitted wave using LVEMM
0.14
Velocity (m/s)
Velocity (m/s)
0.05
0.00 0.01
0.02
0.03
0.00
0.04
0.01
Incident wave Transmitted wave using DDM Transmitted wave using NVEMM Transmitted wave using LVEMM
0.35
0.03
0.04
0.25 0.20 0.15
Incident wave Transmitted wave using DDM Transmitted wave using NVEMM Transmitted wave using LVEMM
0.5 Velocity (m/s)
Velocity (m/s)
0.30
0.02
Time (s)
Time (s)
0.4 0.3 0.2
0.10 0.1
0.05
0.0
0.00 0.00
0.01
0.02
0.03
0.04
0.00
0.01
Time (s)
0.02
0.03
0.04
Time (s)
Fig. 6. Comparison of transmitted waves based on DDM, NVEMM and LVEMM with different amplitude incident waves. (a) A ¼ 0.057 m/s, (b) A ¼ 0.151 m/s, (c) A ¼ 0.340 m/ s and (d) A ¼ 0.529 m/s. (a) Evi f and (b) Zvi f.
1.0
4.1. Wave attenuation
0.9
The transmission coefficient is used to assess the attenuation of the transmitted wave. When an incident wave given in the form of Eq. (37) propagates across a jointed rock mass with a joint spacing of S¼5.83 m, the relations of the transmission coefficient versus the incident wave amplitude and the incident wave frequency are investigated, respectively. The relation between the transmission coefficient and the amplitude of the incident wave based on the LVEMM and the present NVEMM for frequencies of f0 ¼100, 200, 300 and 400 Hz are shown in Fig. 7. It can be seen from Fig. 7 that the transmission coefficient obtained using the LVEMM is independent of the incident wave amplitude, while those obtained by the NVEMM is amplitude dependent. This is reasonable because a nonlinear model is adopted in the NVEMM. For the incident wave with a certain frequency, the transmission coefficient increases with the increase of the incident wave amplitude based on the present NVEMM. All the transmission coefficients converge to 1 if the amplitude of the incident wave is sufficiently large. This phenomenon can be explained that as the incident wave amplitude increases, the joint is compacted nonlinearly and thus the joint stiffness increases accordingly. The stiffness of the joint approaches the stiffness of the rock material when the joint is fully compacted. From Fig. 7, the transmission coefficient for the incident wave with a larger frequency is smaller than that with a smaller frequency, which coincides with
Transmission coefficient T
4. Discussion
LVEMM
NVEMM f =100
0.8
f =200 f =300
0.7
f =400
0.6 0.5 0.4 0.3 0.2 0.0
0.2 0.4 0.6 0.8 Amplitude of incident wave A (m/s)
1.0
Fig. 7. Effect of incident wave amplitude on transmission coefficient.
the geological investigation that a joint filters out stress waves with higher frequencies, but has less effect on the stress waves with lower frequencies [12]. Fig. 8 illustrates the relation between the transmission coefficient and the incident wave frequency with respect to different incident wave amplitudes, i.e., A¼0.05 m, 0.1 m, 0.3 m, 0.5 m and 0.8 m, respectively. It is seen from Fig. 8 that there is a flat low frequency region within which the transmission coefficient is nearly a constant or reduces slightly with frequency increasing.
L.F. Fan et al. / International Journal of Rock Mechanics & Mining Sciences 50 (2012) 11–18
17
6000
1.0
0.6 0.4
Effective Velocity (m/s)
Transmission Coefficient T
5500 NVEMM A=0.80 NVEMM A=0.50 NVEMM A=0.30 NVEMM A=0.10 NVEMM A=0.05 LVEMM
0.8
0.2
5000
NVEMM NVEMM NVEMM NVEMM NVEMM LVEMM
4500 4000 3500 3000 2500
0.0
2000 1500
-0.2 101
102
103
104
105
106
107
101
Frequency of incident wave f0 (Hz) Fig. 8. Effect of incident wave frequency on transmission coefficient.
5500 5000 Effective velocity (m/s)
A=0.80 A=0.50 A=0.30 A=0.10 A=0.05
LVEMM
NVEMM f =100
4500
f =200 f =300 f =400
4000 3500 3000 2500 0.0
0.2
0.4
0.6
0.8
1.0
Amplitude of incident wave (m) Fig. 9. Effect of incident wave amplitude on effective wave propagation velocity.
When the frequency is sufficiently large, the transmission coefficient becomes zero. Between these two extremes, there is a rapid decreasing region of the transmission coefficient. For a fixed incident wave amplitude, the incident waves with a larger frequency result in less transmission. For the same frequency, the incident waves of the larger amplitudes result in a larger transmission coefficient. It is also found from Fig. 8 that as the amplitude of the incident wave decreases, the relationship between the transmission coefficient and the incident wave frequency obtained by the NVEMM approaches to that given by the LVEMM. 4.2. Effective velocity The effective wave propagation velocity in an equivalent medium is defined as the wave traveling length over the corresponding wave traveling time. The wave traveling time is measured as the time difference between the peak velocity of the incident and transmitted waves. Considering the joint spacing S ¼5.83 m, the effective velocity versus the incident wave amplitude for different frequencies of f0 ¼100 Hz, 200 Hz, 300 Hz and 400 Hz are shown in Fig. 9. It is seen that, with a fixed wave frequency, the larger the incident wave amplitude, the larger the effective velocity. The effective velocity approaches the maximum
102 103 104 Frequency of incident wave f0 (Hz)
105
Fig. 10. Effect of incident wave frequency on effective wave propagation velocity.
value when the stiffness of the joint is close to the stiffness of the intact rock. In the LVEMM, however, the effective velocity is independent to the incident wave amplitude. The relations between the effective velocity and the incident wave frequency for different incident wave amplitudes, i.e., A¼0.05 m, 0.1 m, 0.3 m, 0.5 m and 0.8 m are shown in Fig. 10. There is a low frequency region for which the effective velocity is the minimum and it varies slightly in the region. And the effective velocity approaches the wave propagation velocity of the intact rock material at high frequencies, which is v¼5830 m/s in the present analysis. Between these two regions, there is a rapid increase in the effective velocity as the frequency increases. Similar results have been obtained by the displacement discontinuity method [23]. It is also observed from Fig. 10 that as the incident wave amplitude decreases, the relationship between the effective velocity and the incident wave frequency obtained by the NVEMM approaches to that by the LVEMM. The LVEMM can be used to replace the NVEMM approximately if the amplitude of the incident wave is insufficient to cause nonlinear deformation of the joint.
5. Conclusions A nonlinear viscoelastic equivalent medium model is proposed in the present study for the investigation of stress wave propagation in a rock mass consisting equally-spaced parallel joints with nonlinearly deformational behavior. The nonlinear behavior of the rock joint is approximated by a piecewise linear joint model. A method of characteristic line is adopted to solve the governing equations of stress wave propagation in the viscoelastic equivalent continuum. The parameters of the present viscoelastic equivalent model are theoretically derived and the frequencydependence and the amplitude-dependence of the transmission coefficient and effective wave propagation velocity are presented and discussed in detail. It is concluded that stress wave propagation through a nonlinear rock joint is amplitude-dependent. The conventional LVEMM is limited by the linear deformation assumption, and it is not able to describe the dependence of stress wave propagation with the incident wave amplitude. The present NVEMM not only determines accurately the transmission coefficient and the effective velocity when the amplitude of incident stress wave is sufficiently large to mobilize the nonlinear deformation of the
18
L.F. Fan et al. / International Journal of Rock Mechanics & Mining Sciences 50 (2012) 11–18
rock joint, but also degrades to exactly the LVEMM when the amplitude of the incident stress wave is small. Thus, the conventional LVEMM can be regarded as a special case of the present NVEMM. The proposed NVEMM has advantage over the DDM. The comparison of the results obtained by the DDM and the present NVEMM shows that the NVEMM can be used to study stress wave propagation in jointed rock masses as accurate as the DDM. However, because of its continuous nature in an equivalent medium model, the proposed NVEMM becomes more efficient in analyzing stress wave propagation in jointed rock masses when the joint number increases. It should be mentioned that the present study gives only some preliminary results when a jointed rock mass is approximated by a nonlinear viscoelastic equivalent medium. Further researches are still necessary to consider more complicated cases, for example, inclined impingement of a stress wave onto a rock joint, a rock mass with multiple joint sets, wave multiple-reflection effects, shear wave propagation, etc.
Acknowledgments The authors would like to acknowledge the support of the Collaborative Research Fund with Overseas, Hong Kong and Macau Scholars (Grant No 51028801) and the general project (Grant No 51178012) of the National Natural Science Foundation of China.
References [1] Li JC, Ma GW, Zhao J. An equivalent viscoelastic model for rock mass with parallel joints. J Geophys Res 2010;115:B03305. [2] Bandis SC, Lumsden AC, Barton NR. Fundamentals of rock joint deformation. Int J Rock Mech Min Sci Geomech Abstr 1983;20(6):249–68. [3] Cook NGW. Natural joint in rock: mechanical, hydraulic and seismic behaviour and properties under normal stress. Int J Rock Mech Min Sci Geomech Abstr 1992;29(3):198–223. [4] Ma GW, Li JC, Zhao J. Three-phase medium model for filled rock joint and interaction with stress waves. Int J Numer Anal Methods Geomech 2011;35(1): 97–110. [5] Li JC, Ma GW. Experimental study of stress wave propagation across a filled rock joint. Int J Rock Mech Min Sci 2009;46:471–8.
[6] Zhao J, Cai JG, Zhao XB, Li HB. Dynamic model of fracture normal behaviour and application to prediction of stress wave attenuation across fractures. Rock Mech Rock Eng 2008;41(5):671–93. [7] Perino A, Zhu JB, Li JC, Barla G, Zhao J. Theoretical methods for wave propagation across jointed rock masses. Rock Mech Rock Eng 2010;43:799–809. [8] Zhao J, Zhao XB, Cai JG. A further study of P-wave attenuation across parallel fractures with linear deformational behavior. Int J Rock Mech Min Sci 2006;43:776–88. [9] Miller RK. An approximate method of analysis of the transmission of elastic waves through a frictional boundary. J Appl Mech 1977;44:652–6. [10] Schoenberg M. Elastic wave behaviour across linear slip interfaces. J Acoust Soc Am 1980;68(5):1516–21. [11] Zhao J, Cai JG. Transmission of elastic P-waves across single fractures with a nonlinear normal deformational behavior. Rock Mech Rock Eng 2001;34(1): 3–22. [12] Pyrak-Nolte LJ Seismic visibility of fractures. PhD thesis. Berkeley, CA, USA: University of California; 1988. [13] Li JC, Ma GW, Zhao J. Analysis of stochastic seismic wave interaction with a slippery rock fault. Rock Mech Rock Eng 2011;44:85–92. [14] White JE. Underground sound. New York: Elsevier; 1983. [15] Schoenberg M, Sayers CM. Seismic anisotropy of fractured rock. Geophysics 1995;60(1):204–11. [16] Hudson JA. Wave speeds and attenuation of elastic waves in material containing cracks. Geophys J R Astron Soc 1981;64:133–50. [17] Amadei B, Goodman RE. A 3-D constitutive relation for fractured rock masses. In: Proceedings of the International Symposium on the Mechanical Behavior of Structural Media, Balkema; p. 249–68; 1981. [18] Schoenberg M. Reflection of elastic waves from periodically stratified media with interfacial slip. Geophys Prospect 1983;31:265–92. [19] Majer EL, McEvilly TV, Eastwood FS, Myer LR. Fracture detection using P- and S-wave vertical seismic profiling at The Geysers. Geophysics 1988;53(1): 76–84. [20] Sitharam TG. Equivalent continuum analyses of jointed rock mass: Some case studies. in: Ramamurthy T, editor. Engineering in Rocks for slopes, foundations and tunnels. New Delhi, India: Prentice-Hall of India; 2007. p. 518–44. [21] Zhao XB, Zhao J, Cai JG. P-wave transmission across fractures with nonlinear deformational behavior. Int J Numer Anal Methods Geomech 2006;30(11): 1097–112. [22] Li JC, Ma GW, Huang X. Analysis of wave propagation through a filled rock joint. Rock Mech Rock Eng 2010;43:789–98. [23] Pyrak-Nolte LJ, Myer LR, Cook NGW. Anisotropy in seismic velocities and amplitudes from multiple parallel fractures. J Geophys Res 1990;95(B7): 11345–58. [24] Barton NR, Bandis SC, Bakhtar K. Strength deformation and conductivity coupling of rock joints. Int J Rock Mech Min Sci Geomech Abstr 1985;22(3): 121–40. [25] Wang LL, Labibes K, Azari Z, Pluvinage G. Generalization of split Hopkinson bar technique to use viscoelastic bars. Int J Impact Eng 1994;15(5):669–86. [26] Kolsky H. Stress waves in solids. Oxford: Clarendon Press; 1953. [27] Wang LL. Foundations of stress waves. Beijing: National Defense Industry Press; 1985. [28] Kelley CT. Solving nonlinear equations with Newton’s method. Philadelphia: SIAM; 2003.