On the four-dimensional lattice spring model for geomechanics

Go-Feng Zho,Xiodong Hu,Qin Li,Jijin Lin,Guowei M

aState Key Laboratory of Hydraulic Engineering Simulation and Safety,School of Civil Engineering,Tianjin University,Tianjin,300072,China

bSchool of Civil Engineering,Hebei University of Technology,Tianjin,300401,China

Keywords:Lattice spring model(LSM)Fourth-dimensional spatial interaction Fracturing Geomechanics

ABSTRACT Recently,a four-dimensional lattice spring model(4D-LSM)was developed to overcome the Poisson’s ratio limitation of the classical LSM by introducing the fourth-dimensional spatial interaction.In this work,some aspects of the 4D-LSM on solving problems in geomechanics are investigated,such as the ability to reproduce elastic properties of geomaterials,the capability of solving heterogeneous problems,the accuracy on modelling stress wave propagation,the ability to solve dynamic fracturing and the parallel computational efficiency.Our results indicate that the 4D-LSM is promising to deal with problems in geomechanics.

1.Introduction

Rock,concrete and soil are the main research subjects in geomechanics,which are commonly involved in various engineering activities,such as mining,tunnelling and underground engineering.Mechanical responses of these materials are the most concerns for researchers(Dai et al.,2010;Brown,2015;Feng et al.,2016).For decades,the predecessors,e.g.Dr.K.Terzaghi,Dr.E.Hoek and Dr.E.T.Brown,have made fruitful works.However,some fundamental problems are still unclear.For example,an accurate constitutive model for deformation and a rational model(Liu and Carter,2002)for three-dimensional(3D)crack propagation(Ingraffea,1987;Yang et al.,2017)are still under exploring.These fundamental questions are directly related to major engineering safety issues,e.g.tunnel collapse,slope instability,rockburst,and nuclear waste leakage.Geomaterials are usually heterogeneous and their mechanical responses are highly nonlinear,which made them too difficult to be described through pure analytical study.With the development of the computer science,especially the rapid advance of high-performance computing,numerical modelling provides a promising alternative solution for geomechanical researchers.Compared with the theoretical analysis,it is able to consider the complex geometries,the dynamic processes and the nonlinear constitutive responses involved in the related geotechnical phenomena.Compared with the experimental methods,numerical modelling is capable of performing a large number of parameter analyses with reasonable time and relatively low cost.Although many researchers still have doubts on the ability of the numerical modelling to solve practical engineering problems,the numerical simulation seems to be the only feasible solution to understand the complex geotechnical problems.For example,when a geotechnical hazard occurs,the numerical modelling is helpful to find the actual reason through reproducing the observed results with parametric analysis.Actually,Dr.E.T.Brown said that a breakthrough of computational methods for geomechanics is the most promising way to solve various hard problems faced in the actual geotechnical engineering.

Since the 1950s,many computational methods have been developed(Jing,2003;Zhu and Tang,2006;He et al.,2014;Lisjak and Grasselli,2014).The finite element method(FEM)is a typical technology that was firstly developed in the engineering application and then theoretically completed by mathematicians.Its development has experienced the initial popular,the middle trough,and the final mature(a typical growth curve of science and technology).Success of the commercial software of the FEM makes its usagequiteconvenient.Forexample,the ANSYSand ABAQUS are widely used in the geomechanics and geotechnical engineering.In addition to the FEM,the discrete element model(DEM)is another well-known computational method in the geomechanics,which was originally developed to solve the progressive failure and movement of rock masses(Cundall,1971).Now,it has been widely used in many disciplines such as material science and chemical engineering(Zhao,2015).In recent years,the scientific communities are interested in the DEM due to its ability to simulate high complex phenomena.Besides,due to its bottom-up philosophy,the basic principle of the DEM is very intuitive and easy to understand(Chen and Zhao,1998).Similar to the DEM,the lattice spring model(LSM)is a bottom-up computational method as well.However,its intrinsic limitation on representing the full range of Poisson’s ratios(Hrennikoff,1941)restricts the development of this method.In order to solve this problem,researchers have proposed many solutions,e.g.the multi-body shear spring(Zhao et al.,2011)and the nonlocal potential(Chen et al.,2014).In a recent work,this limitation was further solved by introducing the fourth-dimensional spatial interaction(Zhao,2017).

In this paper,we explore the ability of the four-dimensional LSM(4D-LSM)on solving specific geomechanical problems.We focus on applying some aspects related to problems in geotechnical mechanics,e.g.the Poisson’s ratio,material inhomogeneity,wave propagation,3D failure,as well as computational efficiency and parallelism.Its ability to naturally model heterogeneous problems is demonstrated through a numerical simulation on a two-phase bar compression problem.Capabilities of the 4D-LSM on stress wave propagation are verified against an analytical solution for the P-waveand S-wave propagations in one-dimensionalbar problems.Following this,the parallel computing efficiency and ability on handling 3D fracturing of the 4D-LSM are compared with the parallel distinct LSM(DLSM)(Zhao et al.,2013).Finally,we summarise and discuss the advantages and disadvantages of the 4D-LSM and further possible development.

2.Four-dimensional lattice spring model(4D-LSM)

The original LSM developed by Hrennikoff(1941)is intrinsically limited to solve elastic problems with a fixed Poisson’s ratio of 1/3.Its 3D version is onlyable to handle problems with a fixed Poisson’s ratio of 1/4.There are a number of solutions developed and a consensus is reached that noncentral/nonlocal interaction has to be introduced to solve the Poisson’s limitation of classical LSM.Nevertheless,recently,it has been shown that it is possible to overcome the Poisson’s limitationwith only central interaction,but we have to consider the fourth-dimensional interaction.Because our experience is based on the perception of the 3D space,it is hard to directly image the 4D space.The concept of parallel world,a common concept used in science fiction movies,is a straightforward way to explain the idea of the 4D-LSM.In the 4D-LSM,our world is assumed as a hyper-membrane made up from our visual 3D world and an invisible parallel world.Fig.1 illustrates the process of building a 4D computational model.First,build up a 3D lattice model and assign one additional dimension for each particle.Then,make a copy of this 3D model with an offset along the fourthdimensional direction(4D thickness).For regular cubic lattice,the lattice configuration can be viewed through a tesseract(a 4D cube).The interaction between two particles is given as

where k is the spring stiffness,unis the deformation of the spring and nijis the normal vector from particle i to particle j.Compared with classical 3D-LSM,the only difference is that the force and normal vectors have four components.In 4D-LSM,all springs representing the 3D interaction share the same stiffness(k3D),whereas differentspring stiffnesses were assigned for the fourthdimensional interaction,which is characterised by a ratioλ4D.To reproduce the isotropic elasticity,spring stiffnesses of the 4D-LSM have to be assigned through the following equation:

Fig.1.A regular packed 4D-LSM for the uniaxial compression test to extract the elastic properties.

where kα,kβ and kγare the specific fourth-dimensional stiffnesses.In the 4D-LSM,the central finite difference method is used to solve the system equation,which can be simply written as

where V3Dis the volume of the corresponding represented 3D model,l3D，iis the length of the 3D springs,E is the elastic modulus,and νis the Poisson’s ratio.More details of the 4D-LSM can be found in the work of Zhao(2017).

3.Numerical examples

3.1.In fluence of 4D thickness on the elastic prosperities

For the construction of 4D-LSM,we can consider that there is a parallel version of 3D-LSM in the fourth dimension using a parallel world concept.The distance between the 3D model and its parallel version is defined as 4D thickness ratio,which is a ratio of the distance between two models to the lattice length of 3D lattice(regular lattice).As shown in Fig.1,the model has 8000 particles,whose diameter is 1 mm and its elastic modulus is taken as 10 GPa.Auniaxial compression test is conducted for the model inwhichthe bottom surface is fixed in y direction and a loading velocity of 10 mm/s is applied on the upper surface.The calculation of elastic modulus is expressed as

where f is the force applied on the upper surface,Uyis the displacement of upper surface,and L is the length of the cube.

The corresponding Poisson’s ratio can be calculated by the following formula:

where Ux1and Ux2correspond to the displacements of two lateral surfaces,respectively.It is important to pointout that werecord the convergence values of the above parameters during the simulation as the corresponding numerically reproduced elastic modulus and Poisson’s ratio.The results of numerical simulation are shown in Figs.2 and 3.

The elastic modulus increase ratio is defined as the ratio of the elastic modulus obtained from current simulation to that of the corresponding lattice model with 3D interaction only(λ4D=0).Fig.2 shows the results of elastic modulus increase ratio calculated from the numerical simulation.As the 4D stiffness ratio increases,the elastic modulus increase ratio also increases.However,with the increase of the 4D thickness ratio,the increment of elastic modulus increase ratio becomes smaller and smaller.For any 4D stiffness ratio,as the 4D thickness ratio increases,the corresponding elastic modulus increase ratio decreases.When the 4D thickness ratio approaches in finity,the elastic modulus increase ratio approaches 1,which means that the 4D model degenerates into a 3D one.

Fig.2.In fluence of the 4D thickness ratio on the elastic modulus increase ratio.

Fig.3.In fluence of the 4D thickness ratio on the Poisson’s ratio.

Fig.3 shows the results of Poisson’s ratio from numerical simulation.As the 4D stiffness ratio increases,the Poisson’s ratio decreases.However,when the 4D thickness ratio is small,the decrement of Poisson’s ratio is not obvious.With the increase of 4D thickness ratio,the decrementof Poisson’s ratio becomes largerand larger.When the 4D thickness ratio is around 1,the decrement of Poisson’s ratio reaches its maximum.After this,with the increase of 4D thickness ratio,the decrement of Poisson’s ratio becomes smaller and smaller.

Fig.4.A random packed 4D-LSM for the uniaxial compression test to extract the Poisson’s ratio.

Fig.5.In fluence of the lattice number on the Poisson’s ratio.

Fig.6.In fluence of the 4D stiffness ratio on the Poisson’s ratio with different LNs.

Fig.7.Computational model of a two-phase bar under compression.

Fig.8.Comparison between displacement results along the two-phase bar predicted using the FEM and 4D-LSM.

Fig.9.Computational model for the P-and S-wave propagation in an intact bar.

Based on the above research,it can be concluded that the 4D thickness ratio has a great in fluence on the elastic properties.When the 4D thickness is small,the system stiffness increment is large(see Fig.2),and the in fluence of Poisson’s ratio is small.The system with large stiffness needs smaller time step,which is not recommended to use.When the 4D thickness ratio is large,the change of Poisson’s ratiois alsoverysmall.To have a flexiblerange of Poisson’s ratio,the recommended value of the 4D thickness ratio is around 1.

3.2.In fluence of lattice number on the Poisson’s ratio

Fig.10.Analytical solution and numerical prediction of the stress wave propagation problem.

Fig.11.Computational model of a pre-cracked Brazilian disk.

Fig.12.Computational time and speedup of the parallel 4D-LSM and DLSM codes.

In contrast with the regular packed 4D-LSM,a random packed lattice can also be used.As shown in Fig.4,the left one is a random packed particle model and the rightone is the corresponding lattice model.By setting different threshold values for spring bond formation between two particles,models with different lattice numbers(LNs)(the total lattice bonds divided by the total number of particles)can be generated.Through the numerical simulation,it can be seen that the Poisson’s ratio of the 4D model is affected by the LN(see Fig.5).With the increase of the LN,the Poisson’s ratio decreases conversely.The increase of the LN of the model makes the connection tighter.Consequently,the lateral deformation becomes smaller and the Poisson’s ratio becomes smaller correspondingly.When the LN is around 13,the Poisson’s ratio reaches its minimum.After this,as the LN continues to increase,the Poisson’s ratio no longer decreases and remains constant(around 0.25).Other important parameter that in fluences the Poisson’s ratio is the 4D stiffness ratio,as shown in Fig.6.As the LN increases,the Poisson’s ratio decreases.Furthermore,as the 4D stiffness ratio increases,the Poisson’s ratio also decreases,also due to the enhancement of the connections.

Based on the above simulation,it can be concluded that the LN and 4D stiffness ratio have great in fluences on the Poisson’s ratio.When the LN and 4D stiffness ratio are small,the Poisson’s ratio is relatively large.The increases of the parameters all make the decreasing Poisson’s ratio.From Fig.6,we can conclude that the 4D-LSM could cover a wide range of Poisson’s ratio of geomaterials.

Fig.13.Displacement contours along the x direction of the disk predicted by the 4D-LSM and the DLSM(unit:mm).

3.3.Deformation of a two-phase bar

For the strong meshless methods,heterogeneous simulation is a typical problem(Fang et al.,2009).Special processing methods are required to solve the case of two materials;however,there is no need of any special processing for the 4D-LSM.Here,the 4D-LSM is used to solve the deformation problem of a two-phase bar under compression(see Fig.7).The results of 4D-LSM will be compared with FEM to verify the ability of 4D-LSM to deal with heterogeneous model.The computational model is shown in Fig.7.The lengths of the model along y and z axes are 10 mm and the length along x axis is 40 mm.Following the x axis in the negative direction,there is a load on the right surface as 1.21 MPa.The left surface of the model is fixed in all directions(see Fig.7).The resultsare shown in Fig.8.The numerical results of the FEM and 4D-LSM are in good agreement except a little deviation for the condition of E1/E2=10.It should be caused by the different representations of mesh-and particle-based numerical methods,i.e.the material interface cannot be presented precisely as a zero thickness plane in a particle-based model.

Fig.14.Failure patterns of the pre-cracked Brazilian disk predicted by the 4D-LSM and DLSM.

3.4.Stress wave propagation through an intact bar

Fig.9 shows the computational model for the P-and S-wave propagation in an intact bar to analyse the relationship between the Poisson’s ratio and the wave velocity.There are two measuring points A(z=20 mm)and B(z=190 mm)along the z axis of the bar whose length is 200 mm.The lengthsof the barin x andy directions are both 50 mm.The right surface perpendicular to the z axis is non-reflection boundary,and the top and bottom surfaces are fixed in y direction.

We input the P-and S-waves at the left surface along the z axis,respectively,and the computational results are shown in Fig.10.The analytical solution shows that with the increase of Poisson’s ratio,the P-wave velocity increases slowly and the S-wave velocity decreases slowly.The numerical prediction of P-wave velocity by the 4D-LSM is very close to the analytical solution and the average deviation of S-wave velocity is about 3.3%.It can be concluded that the 4D-LSM could be used to deal with the stress wave propagation in geomechanics.

3.5.Crack propagation of a pre-cracked Brazilian disk

The computational model to study the crack propagation of a pre-cracked Brazilian disk is shown in Fig.11.The diameter of the disk is 50 mm and the thickness of the specimen along z axis is 10 mm.In the middle of the disk,there is a crack with length of 20 mm and width of 2 mm.The angle between the crack and the x axis is 30°.The bottom surface boundary of the disk is fixed in y direction and a loading velocity of 500 mm/s is applied at the top surface.The computational model established in 4D-LSM consists of particleswith finite size.The boundaryconditions(the upperand bottom surfaces)are applied to a single layer of particles(see Fig.11).In the model of a pre-cracked Brazilian disk,the numbers of particles assigned to the upper and bottom boundary conditions both are 14,which means that the corresponding width is 7 mm(see Fig.11).

The computer used to analyse this problem has 64 GB memory and its CPU is Intel Xeon E5-2630 v4@2.20 GHz with 20 cores and 40 threads.The DLSM and 4D-LSM are used to analyse the influences of threads number on the computational time and speedup of the computation,respectively.As shown in Fig.12,compared with the DLSM,the 4D-LSM has a shorter computational time and higher efficiency.However,the parallel computation efficiency of the 4D-LSM is lower than that of the DLSM,especially when the number of threads is relativelylarge.The probable reason should be that the efficiency of multi-core parallel computation for inverse matrix(the core computing task of the DLSM)is high.

Crack propagation and fracture of rock are represented as a series of breakage of spring bonds.The bond failure criterion is that deformation between two particles exceeds a given limit value.Fig.13 shows the displacement contour of the disk predicted by the 4D-LSM and the DLSM in steps 160 and 200,respectively.Comparing Fig.13a and b,we can find that the crack begins propagating at the tip of pre-crack at step 160 and the complete coalescence of crack occurs at step 200,both in the DLSM and 4D-LSM.The crack distribution of DLSM and 4D-LSM are very close to each other;however,the displacement value along the x direction is slightly different.

The failure patterns of the pre-cracked Brazilian disk predicted by the 4D-LSM and the DLSM are shown in Fig.14.The overall cracks of 4D-LSM are similar to those of DLSM.The crack paths are both from the tip of pre-crack to the loading end,which are similar to the experimental results of rock and gypsum materials.Nevertheless,both the methods are slightly different from the microscopic view.Generally,the crack surfaces of 4D-LSM are more smooth,which can be verified by 3D views of Fig.14a and b,respectively.In 3D view,the red particles represent the damaged particles that are still in a tension state(corresponding to the macroscopic cracks that naked eyes can see),while the green particles represent the damaged particles that are in a status of closure(corresponding to the macroscopic cracks that are not easy to observe).

4.Conclusions

This paper mainly focused on abilities of the 4D-LSM to recover the proper Poisson’s ratio,consider the material inhomogeneity,solve the wave propagation problem,and handle the dynamic fracturing.From numerical simulation,the suitable 4D thickness is suggested as the average spring length of the 3D lattice.Using irregular lattice model with various LNs,the 4D-LSM can reproduce a proper range of Poisson’s ratios for geomaterials.Capabilities of the 4D-LSM on solving stress wave propagation and progressive fracturing are also demonstrated.Our results also show that,compared with the DLSM,the 4D-LSM has a higher computational efficiency.However,to solve actual problems in geotechnical engineering,the 4D-LSM still needs to be incorporated with specific constitutive models that can consider the nonlinear deformation and failure characteristics of various geomaterials.

Conflicts of interest

We wish to confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have in fluenced its outcome.

Acknowledgments

This research was financially supported by the National Natural Science Foundation of China(Grant No.1177020290).