Modeling stress wave propagation in rocks by distinct lattice spring model

Gaofeng Zhao

Centre for Infrastructure Engineering and Safety,School of Civil and Environmental Engineering,The University of New South Wales,Sydney,NSW 2052, Australia

Modeling stress wave propagation in rocks by distinct lattice spring model

Gaofeng Zhao*

Centre for Infrastructure Engineering and Safety,School of Civil and Environmental Engineering,The University of New South Wales,Sydney,NSW 2052, Australia

A R T I C L EI N F O

Article history:

Received 26 November 2013 Received in revised form

23 March 2014

Accepted 24 April 2014

Available online 24 June 2014

Distinct lattice spring model

In this paper,the ability of the distinct lattice spring model(DLSM)for modeling stress wave propagation in rocks was fully investigated.The inf l uence of particle size on simulation of different types of stress waves(e.g.one-dimensional(1D)P-wave,1D S-wave and two-dimensional(2D)cylindrical wave)was studied through comparing results predicted by the DLSM with different mesh ratios(lr)and those obtained from the corresponding analytical solutions.Suggested values oflrwere obtained for modeling these stress waves accurately.Moreover,the weak material layer method and virtual joint plane method were used to model P-wave and S-wave propagating through a single discontinuity.The results were compared with the classical analytical solutions,indicating that the virtual joint plane method can give better results and is recommended.Finally,some remarks of the DLSM on modeling of stress wave propagation in rocks were provided.

?2014 Institute of Rock and Soil Mechanics,Chinese Academy of Sciences.Production and hosting by Elsevier B.V.All rights reserved.

1.Introduction

Stress wave propagation in rocks is one of the most important issues in rock dynamics,e.g.the damage criteria of rock structures underdynamic loads are generallyregulated according tothreshold values of wave amplitudes:peak displacement,peak particle velocity or peak acceleration.Prediction of stress wave propagation in rocks is also the fundamental requirement in the study of mechanism of seismic events in earthquake.Up to now,a variety of theoretical,experimental and numerical studies have been conducted.For example,Schoenberg(1980)and Pyrak-Nolte et al. (1990)developed analytical solutions to predict the incident wave through a single dry or fully liquid-f i lled fracture using the displacement discontinuous models.Later,these equations were validated by laboratory experiments carried out by Myer et al. (1985)and Suárez-Rivera(1992),respectively.The analytical solutions to interface wave propagation alongside a single failure have been studied by Gu(1994)and Gu et al.(1995),which were also successfully validated by laboratory measurements.Stress wave propagation through a single discontinuity is simple and straightforward.However,stress wave propagation within a medium with multiple joints(a typical situation of rock mass in nature)is much more complex due to multiple ref l ections between separate failures.For this kind of situation,analytical solutions can only be derived under idealized conditions;examples can be found in Cai and Zhao(2000),Zhao et al.(2008)and Li et al.(2010).

To overcomethe limitation of analytical method,more and more numerical methods have been applied for the analysis of stress wave propagation in rocks,e.g.the f i nite element method(FEM) (e.g.Moran,1987),f i nite difference method(FDM)(e.g.Reeshidev and Mrinal,2008),boundary element method(BEM)(e.g.Demirel and Wang,1987),discrete element method(DEM)(e.g.Resende et al.,2010),discontinuous deformation analysis(DDA)(Jiao et al., 2007),discontinuous Galerkin method(DGM)(e.g.Park and Tassoulas,2002),and numerical manifold method(Fan et al., 2013;Zhao et al.,2014).Lattice spring model(LSM)can be viewed as a numerical model based on the concept of bottom-up andone-dimensional(1D)modelingconcept(Wang,2008; Rinaldi,2013).The LSM was originally developed by Hrennikoff (1941)to solve elasticity problems.However,due to computational limitations and the development of FEM,this method was underdeveloped.In recent years,many researchers have renewed their interests in this method due to its advantage in modeling solids fracturing.The LSMs are also adopted for stress wave propagation in rocks,e.g.O’Brien(2008)developed a visco-elastic LSM for seismic wave propagation,and Takekawa et al.(2013)proposeda similar model for stress wave propagation in solid.However,most of these models only considered the rocks as continuous media without addressing the joints/discontinuities.

In this work,the application of distinct lattice spring model (DLSM)(Zhao,2010;Zhao et al.,2011)to stress wave propagation through rocks is discussed.The main contributions of the DLSM are:(1)the restriction on the Poisson’s ratio in classical LSM was removed through a technique to evaluate spring deformation using the local strain technique rather than the particle displacements directly;(2)a close relationship among the spring parameters and the macro-elastic constants,Young’s modulus and Poisson’s ratio is established;and(3)the lattice structures can be both regular and irregular.A few examples of the DLSM on modeling of stress wave propagation through a continuous rock bar were described in Zhao et al.(2011).Verif i cation of DLSM on modeling 1D P-wave propagation through rock masses was studied by Zhu et al.(2011).In this context,a more comprehensive investigation on the ability of the DLSM to model stress wave propagation through rocks is presented,e.g.both 1D P-wave,1D S-wave and two-dimensional(2D) cylindrical wave will be covered.

2.Stress wave propagation by the DLSM

2.1.The model

In DLSM,the material is represented as particles bonded together by springs(see Fig.1).The equation of motion for the system is described as

where u is the vector of particle displacement,[K]is the stiffness matrix,[M]is the diagonal mass matrix,[C]is the damping matrix, and F(t)is the vector of external force.In DLSM,Eq.(1)was solved using the Newton’s Second Law.Details can be found in Zhao(2010) and Zhao et al.(2011).

The input elastic parameters in DLSM are the Young’s modulus and the Poisson’s ratio.During calculation,the spring parameters are calculated from the following equations:

whereknandksare the normal and shear spring stiffness,respectively;EiandEjare the Young’s moduli of the linked particles, respectively;νiandνjare the corresponding Poisson’s ratios;and α3Dis a microstructure geometry coeff i cient of the lattice model expressed as

whereliis the length of theith bond,andVis the volume of the model.

2.2.Viscous boundary condition

Stress wave propagation in a computational model with fi nite boundary causes the wave to be re fl ected and blended with the initial input.It is very dif fi cult to analyze the mixed results.To solve this problem,a non-re fl ection boundary was implemented into DLSM to simulate the computational model without fi nite boundaries.The viscous non-re fl ection boundary condition in DLSM is shown in Fig.2.Three dashpots were placed at particles on the arti fi cial boundary plane to minimize the re fl ected wave.Details on the implementation can be found in Zhao(2010).

Fig.1.Lattice structures in DLSM(Zhao et al.,2011).(a)Simple cubic lattice,(b)Simple cubic II lattice,(c)Simple cubic III lattice,(d)BCC lattice,(e)FCC lattice,and(f)Random lattice.

2.3.Joint/discontinuity in DLSM

The simplest way to represent a discontinuity is to treat it as a thin layer of material with weak mechanical properties(E＇，ν＇)as shown in Fig.3.There is no need to change the original DLSM code but only a few modif i cations on the pre-processor.The stiffness parameters of the discontinuity represented through this method can be obtained as

wheredis the thickness of the layer of weak material.

The idea of virtual joint plane method(Zhao,2010)is original from the smooth-joint contact model(Mas Ivars et al.,2008)as shown in Fig.4.A similar technique was adopted in DLSM to represent the discontinuity.The direction of the spring was changed into the normal direction of the joint plane when a spring was split by a joint plane,and in the meantime,the spring stiffness parameters were modif i ed according to

3.Examples

3.1.Particle size and 1D P-wave and S-wave

Inf l uence of particle size on the numerical accuracy of DLSM modeling of wave propagation problems was studied in this section.Similar work has been conducted for DEM and FEM,e.g.the mesh size inf l uence of UDEC onwave propagation has been studied by Zhao et al.(2008).In order to keep consistent with previous studies,the term called mesh ratio(lr)(ratio of the particle size to the wavelength of input wave)was also adopted as the control parameter.

Fig.4.The smooth-joint contact model(Mas Ivars et al.,2008).

Fig.5.DLSM models for 1D P-wave and S-wave propagation.

The numerical accuracy of 1D P-wave and S-wave representation in DLSM can be verif i ed by simulating a plane elastic wave propagation through an elastic medium.The used DLSM models areshowninFig.5.Thedimensionofthemodelwas 70 mm×140 mm×5 mm,and the diameter of spherical particles was 0.5 mm.The material properties were:density=2120 kg/m3, elastic modulus=27.878 GPa,Poisson’s ratio=0.298,shear wave velocityCs=2250m/s,andcompressionalwavevelocityCP=4200 m/s.A one-cycle sinusoidal wave with an amplitude of 100 mm/s was normally or tangentiallyapplied tothe top boundary and propagates along they-direction through the model.Seven measuring points were positioned in the specimen to record timehistories of the particle velocities(see Fig.5).For P-wave,the left and right side boundaries were f i xed inx-direction.

The wave frequencies of P-wave were taken different values as 0.1 MHz,0.2 MHz,0.5 MHz,1.0 MHz and 2.0 MHz to produce differentlrvalues as 1/82,1/41,1/17,1/8 and 1/4.The percentage error of DLSM on modeling the amplitude of P-wave(the ratio of the difference amplitude between transmitted-wave and inputwave)was used as the index to represent the accuracy of the numerical results.The results of 1D P-wave propagation are shown in Fig.6.It shows that the percentage error decreases with decreasing particle size and increases with the increasing distance from wave source.From these results,it was found that the percentage error will be less than 5%whenlris smaller than 1/41.For S-wave,the wave frequencies were selected as 0.2 MHz,0.1 MHz, 0.05 MHz and 0.025 MHz.The correspondinglrvalues are 1/22,1/ 45,1/90 and 1/180.The DLSM modeling results are shown in Fig.7. The suggestedlrvalue of DLSM modeling of S-wave propagation was obtained as 1/90.In DLSM,to obtain a precise wave form,a value oflr=1/180 is suggested.

Fig.6.Percentage error of wave amplitudes of DLSM modeling of P-wave propagation with differentlrvalues.

3.2.Inf luence of the mesh ratio on 2D wave propagation

The ability of the DLSM to model 2D wave propagation was studied by simulating stress wave propagation through a cylindrical cavity(see Fig.8).The analytical solution of the radial displacement,velocity and stress of the medium can be found in Graff(1979).The wave velocity attenuation ratio along the radial direction can be obtained as

Fig.7.Percentage error of wave amplitudes of DLSM modeling of S-wave propagation with differentlrvalues.

Fig.8.The problem of stress wave propagation from a cylindrical cavity.

Here,the wave attenuation ratio was used as the index to compare DLSM modeling results and the analytical ones.Fig.9 shows the DLSM used to simulate the stress wave propagation through the cylindrical cavity.A cavity with a radius of 10 mm exists in an inf i nite domain.A quarter symmetrical model with a dimension of 100 mm×100 mm×5 mm was used.The particle size was 0.5 mm and a total of 396,840 particles were used to build the computational model.The top and right boundaries were nonref l ection boundaries,while the left and the lower boundaries were symmetrical boundaries.A compressional harmonic wave with amplitude of 100 mm/s was applied at the boundary of the cavity. The wave frequencies were taken as 0.1 MHz and 0.2 MHz to representlrof 1/41 and 1/17,respectively.The mechanical parameters were taken as follows:elastic modulus is 27.878 GPa,Poisson’s ratio is 0.298 and the density is 2120 kg/m3.

In order to quantify the DLSM results,the error for detection point was given as

Fig.9.The used DLSM computational model of the stress wave propagation through cylindrical cavity.

Fig.10.The DLSM modeling results underlrof 1/17 and analytical solution of the wave propagation through cylindrical cavity.

whereAi(DLSM)is the attenuation value of the wave at theith monitoring point predicted by DLSM,andAi(analytical)is the correspondingvalueoftheanalyticalsolution.Theresultsof DLSM modeling and analytical solution are shown in Figs.10 and 11.

The average error was 10.86%for the DLSM withlrof 1/17 and 1.02%for thelrof 1/41.In this sense,the suggestedlrcan be also taken as 1/41 for 2D P-wave propagationproblems.The suggestedlrin the DLSM is smaller than that in the UDEC,e.g.thelrof 1/12 was suggested for the UDEC modeling of P-wave propagation by Zhao et al.(2008).One of the reasons is that the def i nitions of mesh size and particle size in UDEC and DLSM are different.One single element in UDEC includes four sub-triangle elements.In this case, the requirement in UDEC is actuallylr=1/24.For S-wave propagation problem,a strict requirement is set in DLSM(lr=1/90), while the UDEC can still uselr=1/24(the actual ratio).It is concluded that a stricter requirement on particle size is needed for DLSM to model wave propagation than that for UDEC.

3.3.P-wave/S-wave across a single fracture

The theoretical expression of the transmission coeff i cient for an incidentharmonicP-wave/S-waveacrossasinglelinearly deformable fracture in a continuous rock material can be calculated as(Zhao et al.,2008):

Fig.11.The DLSM modeling results underlrof 1/41 and analytical solution of the wave propagation through cylindrical cavity.

Fig.12.DLSMs for P-wave/S-wave incidence.

Fig.13.The modeling results of weak material layer method and analytical solution of P-wave/S-wave propagation through single discontinuity.

where|T1|is the transmission coeff i cient,kis the normal/shear fracturestiffness,ωisthe angularfrequencyof thewave,andzisthe P-wave/S-wave impedance(the product of P-wave/S-wave velocity andmaterialdensity).Itshouldbementionedthattheunitofnormal and shear stiffness of a three-dimensional joint is Pa/m,however,in this work it is taken as Pa for the 2D examples by assuming that the thickness in the analytical solution is 1 m.The fast Fourier transform(FFT)andinversefastFouriertransform(IFFT)canbeused to obtain the analytical solution of a half-cycle sinusoidal wave across a single fracture.Details can be found in Zhao et al.(2008).

The DLSMs for P-wave/S-wave across a single fracture are showninFig.12.Thedimensionofthesemodelswas 70 mm×140 mm×5 mm and the particle size used was 0.5 mm. The material parameters were:elastic modulus=27.878 GPa, Poisson’s ratio=0.298,and the density=2120 kg/m3.A half sinusoidal P-wave/S-wave with frequency of 20 kHz was applied at the top boundary of these models.Thelrwas 1/420 for P-wave propagation problem and 1/220 for S-wave case.

Firstly,the weak material layer method was used to represent the discontinuity.A small ratio of the elastic modulus of the base material is taken as the elastic modulus of the weak material layer, which was 0.005,0.01,0.02,0.05,0.08,0.1,0.3,and 0.5,to produce different normal and shear stiffnesses.The modeling results of theweak material layer method are shown in Fig.13.It points out that the difference between analytical solution and DLSM results was apparent.In order to provide a quantity comparison,the percentage errors between numerical and analytical solutions are listed in Table 1.It can be seen that the error decreases with increasing joint stiffness.The maximum error of weak material layer method is about 9%for P-wave and about 18%for S-wave.So this method is not a good choice for quantitative analysis of wave propagation through discontinuities.

Table 1The errors of weak material layer method on modeling P-wave/S-wave propagation through a single discontinuity.

Fig.14 shows the results of virtual joint plane method.It can be seen that better agreements are obtained.The percentage errors of the virtual joint plane method on modeling P-wave and S-wave propagation are given in Table 2.The maximum error for P-wave is 0.59%and 2.52%for S-wave.This means that the virtual joint plane method is better than the weak material layer method on modeling discontinuities for stress wave propagation problems.

Fig.14.The modeling results of virtual joint plane and analytical solution of P-wave/S-wave propagation through single discontinuity.

Table 2The errors of virtual joint plane method on modeling P-wave/S-wave propagation through a single discontinuity.

4.Conclusions

Abilities of the DLSM to model wave propagation have been studied.The inf l uence of particle size on the numerical error of DLSM modeling of 1D P-wave and 1D S-wave propagation was studied.The suggested mesh ratio(lr)for different conditions was provided.For DLSM modeling of wave problems,the suggestedlrfor P-wave is 1/41 and 1/90 for S-wave.The weak material layer method and virtual joint plane method were used to model P-wave and S-wave propagation through a single discontinuity and compared with the analytical solution.The virtual joint plane method is recommended for modeling discontinuities in the DLSM.

Conf l ict of interest

We wish to con fi rm that there are no known con fl icts of interest associated with this publication and there has been no signi fi cant fi nancial support for this work that could have in fl uenced its outcome.

Acknowledgments

This research is f i nancially supported by the Australian Research Council(Grant No.DE130100457).

Arfken G.Mathematical methods for physicists.3rd ed.Orlando:Academic Press; 1985.pp.604-10.

Cai JG,Zhao J.Effects of multiple parallel fractures on apparent attenuation of stress waves in rock masses.International Journal of Rock Mechanics and Mining Sciences 2000;37(4):661-82.

Demirel V,Wang S.An eff i cient boundary element method for two-dimensional transientwavepropagationproblems.AppliedMathematicalModelling 1987;11(6):411-6.

Fan LF,Yi XW,Ma GW.Numerical manifold method(NMM)simulation of stress wave propagation through fractured rock mass.International Journal of Applied Mechanics 2013.http://dx.doi.org/10.1142/S1758825113500221.

Graff KF.Wave motion in elastic solids.Ohio State University Press;1979.

Gu B,Nihei KT,Myer LR,Pyrak-Nolte LJ.Fracture interface waves.Journal of Geophysical Research 1995;101:827-35.

Gu B.Interface waves on a fracture in rock.PhD Thesis.Berkeley,USA:University of California;1994.

Hrennikoff A.Solution of problems of elasticity by the framework method.Journal of Applied Mechanics 1941;8:A619-715.

Jiao YY,Zhang XL,Zhao J,Liu QS.Viscous boundary of DDA for modeling stress wave propagation in jointed rock.International Journal of Rock Mechanics and Mining Sciences 2007;44(7):1070-6.

Li J,Ma G,Zhao J.An equivalent viscoelastic model for rock mass with parallel joints.Journal of Geophysical Research 2010;115(B3):B03305.

Mas Ivars D,Potyondy DO,Pierce M,Cundall PA.The smooth-joint contact model. In:Schref l er BA,Perego U,editors.The 8th World Congress on Computational Mechanics(WCCM8)/The 5th European Congress on Computational Methods in Applied Sciences and Engineering.Barcelona:International Center for Numerical Methods in Engineering(CIMME);2008.Paper No.a2735.

Moran B.A f i nite element formulation for transient analysis of viscoplastic solids with application to stress wave propagation problems.Computers and Structures 1987;27(2):241-7.

Myer LR,Hopkins D,Cook NGW.Effects of contact area of an interface on acoustic wave transmission characteristics.In:Rock Mechanics Symposium,vol.1. Boston;1985.pp.565-72.

O’Brien GS.Discrete visco-elastic lattice methods for seismic wave propagation. GeophysicalResearchLetters2008;35:L02302.http://dx.doi.org/10.1029/ 2007GL032214.

Park SH,Tassoulas JL.A discontinuous Galerkin method for transient analysis of wave propagation in unbounded domains.Computer Methods in Applied Mechanics and Engineering 2002;191(36):3983-4011.

Pyrak-Nolte LJ,Myer LR,Cook NGW.Transmission of seismic waves across single natural fractures.Journal of Geophysical Research 1990;95:17-38.

Reeshidev B,Mrinal KS.Finite-difference modelling of S-wave splitting in anisotropic media.Geophysical Prospecting 2008;56(3):293-312.

Resende R,Lamas LN,Lemos JV,Cal?ada R.Micromechanical modelling of stress waves in rock and rock fractures.Rock Mechanics and Rock Engineering 2010;43(6):741-61.

Rinaldi A.Bottom-up modeling of damage in heterogeneous quasi-brittle solids. Continuum Mechanics and Thermodynamics 2013;25:359-73.

Schoenberg M.Elastic wave behavior across linear slip interfaces.Journal of Acoustics Society of America 1980;68:1516-21.

Suárez-Rivera R.The inf l uence of thin clay layers containing liquids on the propagation of shear waves.PhD Thesis.Berkeley,USA:University of California;1992.

Takekawa J,Mikada H,Coto TN,Sanada Y,Ashida Y.Coupled simulation of seismic wave propagation and failure phenomena by use of an MPS method.Pure and Applied Geophysics 2013;170:561-70.

Wang JH.One-dimensional dynamical modeling of earthquakes:a review.Terrestrial Atmospheric and Oceanic Sciences 2008;19:183-203.

Zhao GF,Zhao XB,Zhu JB.Application of the numerical manifold method for stress wave propagation across rock masses.International Journal for Numerical and Analytical Methods in Geomechanics 2014;38(1):92-110.http:// dx.doi.org/10.1002/nag.2209.

Zhao XB,Zhao J,Cai JG,Hefny AM.UDEC modelling on wave propagation across fractured rock masses.Computers and Geotechnics 2008;35:97-104.

Zhao GF,Fang J,Zhao JA.3D distinct lattice spring model for elasticity and dynamic failure.International Journal for Numerical and Analytical Methods in Geomechanics 2011;35:859-85.

Zhao GF.Development of micro-macro continuum-discontinuum coupled numerical method.PhD Thesis.Switzerland:EPFL;2010.

Zhu JB,Zhao GF,Zhao XB,Zhao J.Validation study of the distinct lattice spring model(DLSM)on P-wave propagation across multiple parallel joints.Computers and Geotechnics 2011;38:298-304.

*Tel.:+61 2 9385 5022;fax:+61 2 9385 6139.

E-mail address:gaofeng.zhao@unsw.edu.au.

Peer review under responsibility of Institute of Rock and Soil Mechanics,Chinese Academy of Sciences.

Rock

Stress wave

Verif i cation