Modeling Creep Fracture in Rock by Using Kelvin Discretized Virtual Internal Bond

Discretized virtual internal bond (DVIB) is a lattice model, which is composed of bond cells. Each bond cell has a finite number of bonds.-e DVIB is used to model the creep fracture. It is done by introducing a viscous bond to the original hyperelastic DVIB.-e hyperelastic bond is parallel coupled with a viscous bond together, forming a hybrid hyperelastic-Kelvin bond.-e hyperelastic bond reflects the microfracture mechanism, whereas the viscous bond reflects the creep mechanism. Based on this hyperelastic-Kelvin bond, the constitutive relation of a cell is derived. -e microbond parameters are calibrated based on the ideal cell approach. -e simulation results suggest that this method can represent the typical features of creep and can simulate the creep fracture.-emerit of this method lies in that the complicated 3D macrocreep problem is reduced to the 1D microbond creep problem. No creep law is previously derived.-emacrocreep fracture behavior is the natural response of the assembly of the micro hyperelastic-Kelvin bonds.


Introduction
e creep behaviors of rock play an important role in geotechnical engineering.Subjected to a constant load, the deformation of rock evolves with time.e rock fails with creep fracture propagating further.How to model the creep fracture has been an important and tough problem.So far, many theoretical creep models have been proposed, for example, the phenomenological models [1][2][3] and the micromechanicsbased models [4][5][6][7][8][9][10][11][12].In these models, the constitutive relation is derived from some time-dependent mechanisms.e theoretical models are usually based on a 1D elasticity-viscosity coupled component.is 1D component is usually assumed to subject to a fixed load, and then, the creep law is derived under the specific boundary conditions.To simulate the 3D creep behaviors, the 1D constitutive law of creep is extended to the 3D version with certain hypothesis.In the finite element method implementation, each element is enforced to follow this creep law.Although it can simulate some macrocreep behaviors, there are some limitations with this method.e macrocreep is defined as the strain evolution with time when the interesting body is subjected to special boundary conditions.However, for a given microcomponent, it is not subjected to the same creep boundary conditions.For example, for an element at the vicinity of a microcrack tip, it is not subjected to a fixed load although the whole interesting body is subjected to a fixed load on the macroscale level.With the crack opening, the load that the element bears varies.us, for a given microelement, it should not follow the creep law that is derived from special boundary conditions substantially.If we enforce each element to comply with the creep law, it restrains the degree of freedom of an element seriously.In addition, the extension from the 1D creep law to the 3D case is based on a certain simplification, which is not rigorous in theory.
To overcome these limitations of the conventional continuum method, we use the lattice approach to model the creep fracture in this paper.
e lattice model has many advantages in fracture simulation and has been extensively applied, for example, [13][14][15][16][17][18][19][20][21].e DVIB [22] is distinctively different from other lattice models in that its discrete structure is composed of bond cells.Each bond cell can take any geometry with any number of bonds.It is capable of capturing the mesostructural characteristics of rock since the rock is composed of mineral grains on mesoscale.Since the DVIB was proposed in 2013, the model has been well developed.For example, Zhang and Chen [23] used the modi ed Stillinger-Weber potential to describe the bond cell energy to make the DVIB capable of representing the various Poisson's ratios.Zhang et al. [24] developed the elastic-brittle DVIB to simulate the spalling fracture.e plastic and fracture energyembedded DVIB has been developed [25].In this paper, we extend the DVIB to the creep fracture cases.It is done by introducing the viscous bond into the DVIB.

Brief Introduction to the DVIB
DVIB [22] considers the material to consist of bond cells, as shown in Figure 1.Each cell can take any geometry with any number of bonds.Each bond is characterized by a bond potential Φ Φ(l).So, the total strain energy of the bond cell is where l is the deformed bond length and Ω is the bond number in a cell, de ned as Ω N(N − 1)/2, with N being the node number of a cell.By (1), the nodal force and tangent sti ness matrix of the bond cell are derived: where u i is the displacement component.
e physical parameter of bond potential is calibrated as where E is Young's modulus; V is the cell volume; l 0 is the undeformed bond length; and λ is a coe cient, with λ 6, λ 3, and λ 16/5 for 3D, plane stress, and plane strain cases, respectively.

Constitutive Model.
e original DVIB bond is hyperelastic.So, it cannot capture the creep properties of rock.To enable the DVIB to simulate the creep fracture, we introduce a viscous component, as shown in Figure 2. In the hybrid bond, the viscous component is linear, while the elastic one is hyperelastic.Here, the hybrid bond is called the hyperelastic-Kelvin bond, whose potential can be written as where Φ(l) is the hyperelastic potential of the bond; _ l is the bond deformation ratio, _ l zl/zt; and η is the microviscosity coe cient, de ned as η F c / _ l, with F c being the internal force induced by the deformation rate of the bond.
By (4), the nodal force vector is derived as ( e tangent sti ness and viscous matrix are, respectively, derived as Bond cell 2 Advances in Civil Engineering

Parameter Calibration.
In the original DVIB [22], the ideal bond cell approach was proposed to calibrate the bond parameters.e initial bond stiffness was calibrated as Φ ″ (l 0 ) � 6EV/(Ωl 2 0 ) for the 3D cases and Φ ″ (l 0 ) � 3EV/(Ωl 2 0 ) for the in-plane stress case and Φ ″ (l 0 ) � 16EV/(5Ωl 2 0 ) for the in-plane strain case.In this paper, the microbond viscosity is also calibrated based on the ideal bond cell approach.In an ideal cell, the bond number is large enough, and the bonds are uniformly distributed.
erefore, the strain energy density of this cell is written as where the operator for the 2D cases.e bond length and its ratio in terms of strain tensor are written as where ε is the strain tensor and ξ is the orientation vector of the bond, defined as ξ � [sin θ cos ϕ, sin θ sin ϕ, cos θ] T for 3D cases and ξ � [cos θ, sin θ] T for 2D cases.By ( 7) and ( 8), the stress tensor is derived as e tangent modulus tensor is e tangent viscosity tensor is Substituting ( 8) into ( 10) and ( 11), we can obtain the tangent stiffness tensor and viscosity tensor under the undeformed state as In the calibration of the initial bond stiffness Φ ″ (l 0 ), Poisson's ratio is fixed, namely, ] � 0.25.Analogous to the initial bond stiffness calibration, we define the macroviscosity coefficient η as To relate the microviscosity to the macro one, we expand (13) to the 3D form, analogous to the elastic matrix with ] � 0.25: where Ω C complies with the relationship According to the tangent viscosity tensor, the viscosity matrix obtained by integrating (12) is By equaling (14) to (16), the microviscosity coefficient is calibrated as By the same calibration process, the parameters in 2D case are calibrated as

Governing Equation.
e general governing equation for the creep problem can be written as where the restoring nodal force F(u, _ u) includes the effect of viscosity.To solve this equation, we use the central difference method.e relationship of quantities at the two immediate time steps complies with Advances in Civil Engineering where the subscript "0" means the current time step, while "1" means the next time step; Δt means the time interval; and θ is a coe cient.When θ ≥ 0.5, this algorithm is unconditioned stable.In this paper, we take θ 0.75.With (20), the governing equation of ( 19) is discretized into To solve (21), the Newton-Raphson method is adopted, whose iteration algorithm is and the quantities are updated at the next time step as

Computation of Bond
Quantities.Take the single bond shown in Figure 3 as an example to show the computation of the sti ness matrix and viscous matrix.Let X, Y denote the bond vector in the reference and current con guration, respectively, X X II − X I , Y Y II − Y I , with X I , Y I and X II , Y II being the coordinates of the node I and node II in the reference and current con guration, X (x 1 , x 2 , x 3 ), Y (y 1 , y 2 , y 3 ).Let u I , u II denote the displacement vector of the node I and node II.en, X, Y are correlated by Y X + u II − u I .e bond length in the reference and cur- T denotes the displacement of two nodes.For a single bond, the nodal force is derived as and the sti ness matrix is e viscous matrix is e derivatives of bond quantities are Figure 3: A single bond in the reference and the current con guration. 4 Advances in Civil Engineering (27)

Verification and Simulation Examples
5.1.Convergence Check.As a constitutive model, its simulation results should be convergent to the analytical solutions with the decrease of the element size.To check whether this model possesses this property, we simulate a macrocreep case with di erent meshing size schemes.e simulation specimen is shown in Figure 4(a), which is subjected to a constant axial tension, namely, σ 0 5 MPa.e specimen is discretized with irregular triangle cells, as shown in Figure 4(b).e material parameters are E 40 GPa and η 1000 GPa•s.e analytical solution of this case is  Advances in Civil Engineering e elastic potential function used for this simulation is where A is the microbond parameter, calibrated as A 6EV/(Ωl 2 0 ) for the 3D case, A 3EV/(Ωl 2 0 ) for the 2D in-plane stress case, and A 16EV/(5Ωl 2 0 ) for the 2D inplane strain case.e comparison between the analytical and simulated strain is shown in Figure 4(c).It is seen that, with the decrease of the element size, the simulated results approach to  is suggests that the present model is stable and convergent.

Simulation of the Impact Test.
To examine whether the present method can capture the viscosity mechanism, we simulate a long-bar impact test reported in [26].e test is shown in Figure 5(a), where a pendulum impacts a long bar and then the strain at a prescribed point is recorded by the sensor.e material density of the bar is ρ � 2168.5 kg/m 3 , Young's modulus E � 8.4 Gpa, Poisson's ratio v � 0.29, and viscosity coefficient η � 100 KPa•s.e bar size is 50 mm × 50 mm × 1000 mm, whose aspect ratio can ensure the onedimensional wave propagation.In the numerical simulation, different impact speeds, namely, V 1 � 1.3 m/s, V 2 � 2.0 m/s, and V 3 � 2.5 m/s, are adopted.Here, we take the twoparameter cohesive force law as the hyperelastic bond potential, which is where A and B are the microparameters and A is calibrated through (3).B is a scale parameter, defined as B � ε t l 0 , with ε t being the strain value at the peak stress of the uniaxial stress-strain curve.e simulated results are shown in Figure 5(b), which indicates that the simulated results are reasonably consistent with the tested one.To confirm that it is the viscosity, not the elasticity, that leads to these results, we simulate this example by setting the viscosity coefficient η � 0, which completely reduces to the hyperelastic situation.e simulated results are shown in Figure 5(c).It is seen that the simulated strain at the prescribed point oscillates and then drops down to zero quickly if no viscosity is considered. is is reasonable because there is no energy dissipation to resist the elastic deformation in the pure hyperelastic case.So, it is the viscosity that is responsible for the simulated behaviors in Figure 5(b).is example suggests that the present method can reflect the viscosity mechanism to great extent.

Creep Fracture Simulation.
To check the performance of this model in creep fracture simulation, we simulate a center-cracked plate subjected to a constant uniaxial compressive stress σ 0 and Young's modulus E � 40 GPa, as shown in Figure 6(a).
e dimension of the specimen is 1.0 m × 1.0 m, and the crack inclination is 45 °.Take the twoparameter cohesive force law (30) as the hyperelastic bond potential.
To examine the effect of viscosity coefficient on the creep behaviors, we take η � 500 GPa•s, η � 1000 GPa•s, and η � 1500 GPa•s, respectively, and the compressive stress σ 0 � 6.5 MPa. e simulated strain is shown in Figure 6(b).It is seen that the general profile of the strain-time curve presents the typical features of creep.With the increase of the viscosity coefficient, the acceleration deformation stage is delayed.is agrees with the general regularity of creep.
To examine the effect of load magnitude, we take σ 0 � 5.0 MPa, σ 0 � 6.5 MPa, σ 0 � 7.0 MPa, and σ 0 � 7.5 MPa, respectively, and η � 1000 GPa•s.e simulated strain is shown in Figure 6(c).From Figure 6(c), it is seen that when the stress is small, say σ 0 � 5.0 MPa, the strain does not evolve further when it reaches a certain level.However, with the increase of stress, the strain evolves faster and faster.is also agrees with the typical features of creep.e simulated fracturing process with η � 1000 GPa•s and σ 0 � 6.5 MPa is shown in Figure 7, which demonstrates that the specimen fails in the fracture mode.

Conclusions
By introducing a viscous component into a hyperelastic bond, a hyperelastic-Kelvin bond is constructed.Based on this bond, the constitutive relation of a bond cell is derived.
e hyperelastic part of the hybrid bond can reflect the microfracture mechanism, and the viscosity part can reflect the creep mechanism; the presented DVIB can simulate the creep fracture behaviors.e simulation results suggest that this method can represent the typical features of creep.Since the Kelvin model has some limitations to describe the creep, the present model is only used to simulate some simple creep fractures.To reflect a more complicated creep fracture, a more elaborate component should be embedded into the DVIB.Advances in Civil Engineering

Figure 5 :
Figure 5: Simulation of the long-bar impact test: (a) test system after Niu and Zhu [26]; (b) comparison between the simulated and tested strain in [26]; (c) comparison of the simulation results with and without viscosity.

Figure 6 :
Figure 6: Simulation of creep fracture: (a) simulation specimen; (b) e ect of viscosity on creep; (c) e ect of external stress on creep.