A modified virtual multi-dimensional internal bonds model for geologic materials

In fracture simulation field, developing new models with more realistic microstructure is a flourishing tendency. Xu and Needleman (1994) [1] developed the cohesive surface models by introducing a series of discrete cohesive surfaces during the finite element discretization. Those cohesive surfaces have limited cohesive strength and limited work of fracture. And cracks are formed and spread along element boundaries. Gao and Klein (1998) [2] proposed the virtual internal bond (VIB) model, which is a different approach by introducing virtual internal bonds with stretch stiffness attaching erratically dispersed material particles inside the solid and statistically averaging the stochastic spatial network of bonds. The VIB uses the Cauchy-Born rule of crystal elasticity to obtain macroscopic collective behaviour of the stochastic bond network inside solid material, and it only considers the stretch energy of virtual internal bond and has only one type of bond possessing only one kind of stiffness coefficient (stretch stiffness), which lead to free rotation of bonds. Volokh and Gao (2005) [3] proposed a modified virtual internal bond (MVIB) model with bonds owning both stretch and bend stiffness. In the MVIB model, the Green strain tensor is decomposed into spherical and deviatoric parts and the stretch deformation of bonds are decomposed into spherical dilatational deformation and deviatoric deformation to calculate bonds’ strain energy. Zhang and Ge (2005) [4, 5] proposed another modified version of VIB, which introduced an R-bond to restrict the rotation freedoms of material particles. The new model is called as virtual multi-dimensional internal bonds (VMIB) model, and it only can be used to material with Poisson ratio lower than 0.25. Zhang and Gao (2012) [6] proposed an augmented virtual internal bond (AVIB) model using the Xu-Needleman potential to describe bonds’ energy. The AVIB decomposed bond’s strain into stretch strain and shear strain, and then calculates the bond strain energy by the Xu-Needleman potential. Zhang et al. (2014) [7] used the Stillinger-Weber potential, which is a combination of twoand three-body interactions, to calculate bonds energy in VIB model. All the above-mentioned models could only be applied to material with Poisson ratio smaller than 0.25, which is a severely restriction to the further development of virtual internal bond theory.


Introduction
In fracture simulation field, developing new models with more realistic microstructure is a flourishing tendency.Xu and Needleman (1994) [1] developed the cohesive surface models by introducing a series of discrete cohesive surfaces during the finite element discretization.Those cohesive surfaces have limited cohesive strength and limited work of fracture.And cracks are formed and spread along element boundaries.Gao and Klein (1998) [2] proposed the virtual internal bond (VIB) model, which is a different approach by introducing virtual internal bonds with stretch stiffness attaching erratically dispersed material particles inside the solid and statistically averaging the stochastic spatial network of bonds.The VIB uses the Cauchy-Born rule of crystal elasticity to obtain macroscopic collective behaviour of the stochastic bond network inside solid material, and it only considers the stretch energy of virtual internal bond and has only one type of bond possessing only one kind of stiffness coefficient (stretch stiffness), which lead to free rotation of bonds.Volokh and Gao (2005) [3] proposed a modified virtual internal bond (MVIB) model with bonds owning both stretch and bend stiffness.In the MVIB model, the Green strain tensor is decomposed into spherical and deviatoric parts and the stretch deformation of bonds are decomposed into spherical dilatational deformation and deviatoric deformation to calculate bonds' strain energy.
Zhang and Ge (2005) [4,5] proposed another modified version of VIB, which introduced an R-bond to restrict the rotation freedoms of material particles.The new model is called as virtual multi-dimensional internal bonds (VMIB) model, and it only can be used to material with Poisson ratio lower than 0.25.Zhang and Gao (2012) [6] proposed an augmented virtual internal bond (AVIB) model using the Xu-Needleman potential to describe bonds' energy.The AVIB decomposed bond's strain into stretch strain and shear strain, and then calculates the bond strain energy by the Xu-Needleman potential.Zhang et al. (2014) [7] used the Stillinger-Weber potential, which is a combination of two-and three-body interactions, to calculate bonds energy in VIB model.
All the above-mentioned models could only be applied to material with Poisson ratio smaller than 0.25, which is a severely restriction to the further development of virtual internal bond theory.

Constitutive model for non-linear elastic material
In VIB theory, the solid materials are consisted of discrete and random mass particles in micro level, while they don't have to be atoms, and virtual internal bonds are added to describe the interactions between particles (Gao and Klein, 1998) [2].Then the deformation of the continuum displacement field is placed in the deformation of the crystallite correspondently through the Cauchy-Born rule, that is to say, the atomic positions are connected to the continuum fields via the local deformation gradient (Tadmor et al., 1996) [8].The local deformation gradient F is applied to the undeformed crystal lattice basis and rebuilding the crystal through the altered base vectors to obtain the deformed crystal structure (see Fig. 1).In this way, each continuum representative volume element is represented by infinite crystals undertaking uniform deformation.
Fig. 1 The Cauchy-Born rule In this modified model, the virtual internal bond possesses both stretch stiffness and shear stiffness.Fig. 2 shows the microstructure used in this modified model.Virtual internal bond refers to the radial link between two mass particles.Therefore each bond has a unique orientation.If all the virtual internal bonds within one unit of volume of a material are rearranged by their orientations in the sphere coordinate system (Fig. 3).And then a virtual internal bond's orientation ξ would be ex- pressed as: Fig. 3 Virtual internal bond's orientation and the unit sphere In this model, the deformation of each virtual internal bond is decomposed into three parts: the spherical stretch strain n  , the deviatoric stretch strain n  and the shear strain t  .The spherical stretch strain is the bond's normal strain owing to spherical dilatational deformation, the deviatoric stretch strain is the bond's normal strain owing to the deviatoric deformation and the shear strain is the bond's tangential strain due to the whole deformation.And in the case of small deformation, the expressions of the spherical stretch strain, the deviatoric stretch strain and the shear strain are:

T T T ξ ε εξ ξ εξ
(2) In the original virtual multi-dimensional internal bonds model, the stretch strain is the summation of the deviatoric stretch strain and the spherical stretch strain.And its shear parameter would be negative when material Poisson ratio is greater than 0.25.This is mainly caused by the excessive stretch energy.The spherical stretch strain is the major stretch strain of virtual internal bonds, so, in this modified model, it would be considered as an effective stretch strain for bond's stretch energy.
The strain energy stored in virtual internal bonds consists of stretch energy n U and shear energy t U .As- suming that the stretch stiffness is function of spherical stretch strain n  and shear stiffness is function of shear here k 0 is initial stretch stiffness and r 0 is initial shear stiffness.
Assuming that the material is isotropic, then the bond density and stiffness should be uniform in spatial.So the average bond energy for per surface area of the unit sphere (Fig. 3) within one unit of volume of a material is here k and r are bonds' average stretch and shear material parameters for per surface area of unit sphere, respectively.In quasi-continuum theory, the strain energy density must be equivalence in order to make a representative volume element of discrete microstructure mechanically equivalent to a continuum volume element.Therefore the macro strain energy density Φ should be equal to the summation of the entire bonds' stretch and shear energy (Fig. 3): Based on the work conjugate principle, the Caustress tensor ij  is (Ogden, 1984)  [9]: So that with Eq. ( 4), the derivatives of t W and n W would be: And the partial derivatives of n  and 2 t  in Eqs. ( 6) and ( 7) would be: , , sin cos sin sin cos

Incremental constitutive model under triaxial compression
Assuming that materials are isotropic elastic under triaxial compression incremental deformation, and tangent Young's modulus t E and tangent Poisson ratio t  are the material elastic parameters in each stage of deformation.These two parameters change value related with stages of deformation.In geologic materials, basing on the Duncan-Chang non-linear elastic model [10], , tt Ev are given as: , here i E is the initial tangent modulus; S is the stress level; f R is the failure ratio; i v is the initial tangent Poisson ratio; 13 ,  are the major and minor principal stresses; D is the material parameter.
In incremental deformation two tangent parame-ters , tt kr are used to replace , kr and: Substituting Eq. ( 11) into Eq.( 8): Basing on FEM analysis, the tensor ijmn C could be expressed in a matrix form: And in isotropic linear elastic material the stressstrain relation could be expressed as: Substituting Eqs. ( 9) and ( 12) into Eq.( 7), and then integrating Eq. ( 7) to obtain ijmn C : By combining Eq. ( 14) with Eq. ( 15), the relation between , tt kr and the linear elastic material parameters (tangent Young's modulus t E and tangent Poisson ratio t  ) could be obtained:

5
. 41 From Eq. ( 16), both the two material parameters , tt kr remain positive when tangent Poisson ratio t  is greater than 0.25.And the limit Poisson ratio of this model is consistent with the material limit Poisson ratio 0.5.

Applying to the uniaxial tension process in brittle rock materials
In consideration of simple loading case, the nonlinear elastic model can be applied to simulate uniaxial tension process in isotropic brittle rock materials.And in order to simulate this process, those functions     here constants A, B are used to fit material properties ( ,0 AB ).
The non-linear elastic constitutive Eqs. ( 6) and ( 7) are combined with the Eqs.( 8), ( 9) and ( 17) and then embed into the ANSYS.Those material constitutive equations are written into the ANSYS UserMat Subroutine.The Us-erMat subroutine is used to define a material's stress-strain relationship and applies to any analysis procedure involving mechanical behavior.The subroutine is called at every material integration point of the elements during the solution phase.The program passes in stresses, strains, and state variable values at the beginning of the time increment and strain increment at the current increment, then updates the stresses and state variables to the appropriate values at the end of the time increment [12].
With the UserMat subroutine, some numerical simulations of uniaxial tension in geologic materials are performed in ANSYS.In these simulations, a cylindrical specimen is fixed in one side and then applied with a tensile displacement in the other side.Axisymmetric solid plane182 element is used to mesh the specimen (Fig. 4 shows the loading conditions and the mesh of elements).Fig. 4 The loading conditions and the mesh of elements The real tensile tests of Tako Sandstone, Sanjome andesite and Kimachi Sandstone were done in article [13,14], and the test data were cited here to compare with the numerical simulation results.Kimachi Sandstone k = 3.531 GPa, r = 1.810GPa, A = 2.165e-4, B = 7.828e-7.As shown in Fig. 5, the stressstrain curves given by numerical simulation linearly increase at first and then concave upward to the peak strength point, after that the stress drops almost vertically, then a distinct transition appears and the stress starts to decrease slowly until it gradually levels off near a very low constant value.The curve of Sanjome andesite is very much close to the real test data, and it fits well with the experimental results.The curves of Tako Sandstone and Kimachi Sandstone are close to the test data except at postpeak part.

Conclusions
The modified virtual multi-dimensional internal bonds model proposed in this paper uses the spherical stretch strain and the shear strain to calculate the bonds' stretch energy and shear energy respectively, so as to overcome the excessive stretch energy in the VMIB model.By comparing with the Duncan-Chang non-linear elastic model, the micro material parameters , tt kr in the modified model are obtained and expressed by the macro material parameters (the tangent Young's modulus t E and the tan- gent Poisson ratio t  ).These two micro material parame- ters remain positive when the material Poisson ratio is greater than 0.25.The modified model successfully expands the limit Poisson ratio of the VMIB model to the material limit Poisson ratio 0.5.
A non-linear constitutive model is given in chapter 4.Under small deformation and simple loading, this non-linear constitutive model is applied to simulate the uniaxial tension process in brittle rock materials and it fits well with the experimental data except at the post-peak part.