Advances in Civil Engineering
Volume 2018 (2018), Article ID 8042965, 8 pages
https://doi.org/10.1155/2018/8042965

Research Article

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

Wangyang He and Zhennan Zhang

School of Naval Architecture, Ocean, and Civil Engineering, Shanghai Jiao Tong University, Shanghai 200240, China

Correspondence should be addressed to Zhennan Zhang; zhennanzhang@sjtu.edu.cn

Received 22 December 2017; Revised 31 January 2018; Accepted 13 February 2018; Published 5 April 2018

Academic Editor: Gaofeng Zhao

Copyright © 2018 Wangyang He and Zhennan Zhang. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Discretized virtual internal bond (DVIB) is a lattice model, which is composed of bond cells. Each bond cell has a finite number of bonds. The DVIB is used to model the creep fracture. It is done by introducing a viscous bond to the original hyperelastic DVIB. The hyperelastic bond is parallel coupled with a viscous bond together, forming a hybrid hyperelastic-Kelvin bond. The 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. The microbond parameters are calibrated based on the ideal cell approach. The simulation results suggest that this method can represent the typical features of creep and can simulate the creep fracture. The merit 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. The macrocreep fracture behavior is the natural response of the assembly of the micro hyperelastic-Kelvin bonds.

1. Introduction

The creep behaviors of rock play an important role in geotechnical engineering. Subjected to a constant load, the deformation of rock evolves with time. The 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–3] and the micromechanics-based models [4–12]. In these models, the constitutive relation is derived from some time-dependent mechanisms. The theoretical models are usually based on a 1D elasticity-viscosity coupled component. This 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. The 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. Thus, 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. The lattice model has many advantages in fracture simulation and has been extensively applied, for example, [13–21]. The 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 modified 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. The plastic and fracture energy-embedded 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.

2. 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 . So, the total strain energy of the bond cell is

where is the deformed bond length and is the bond number in a cell, defined as , with being the node number of a cell.

Figure 1: Modeling method of the DVIB.

By (1), the nodal force and tangent stiffness matrix of the bond cell are derived:

where is the displacement component.

The physical parameter of bond potential is calibrated as

where is Young’s modulus; is the cell volume; is the undeformed bond length; and is a coefficient, with , , and for 3D, plane stress, and plane strain cases, respectively.

3. Hyperelastic-Kelvin Bond-Based DVIB

3.1. Constitutive Model

The 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 is the hyperelastic potential of the bond; is the bond deformation ratio, ; and is the microviscosity coefficient, defined as , with being the internal force induced by the deformation rate of the bond.

Figure 2: Hyperelastic-Kelvin bond.

By (4), the nodal force vector is derived as

The tangent stiffness and viscous matrix are, respectively, derived as

3.2. Parameter Calibration

In the original DVIB [22], the ideal bond cell approach was proposed to calibrate the bond parameters. The initial bond stiffness was calibrated as for the 3D cases and for the in-plane stress case and 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. Therefore, the strain energy density of this cell is written as

where the operator denotes for the 3D cases and for the 2D cases. The 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 for 3D cases and for 2D cases.

By (7) and (8), the stress tensor is derived as

The tangent modulus tensor is

The 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 , Poisson’s ratio is fixed, namely, . 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 :

where 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

4. Numerical Implementation

4.1. Governing Equation

The general governing equation for the creep problem can be written as

where the restoring nodal force includes the effect of viscosity. To solve this equation, we use the central difference method. The relationship of quantities at the two immediate time steps complies with

where the subscript “0” means the current time step, while “1” means the next time step; means the time interval; and is a coefficient. When , this algorithm is unconditioned stable. In this paper, we take .

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

4.2. Computation of Bond Quantities

Take the single bond shown in Figure 3 as an example to show the computation of the stiffness matrix and viscous matrix. Let denote the bond vector in the reference and current configuration, respectively, , , with and being the coordinates of the node I and node II in the reference and current configuration, , . Let denote the displacement vector of the node I and node II. Then, are correlated by . The bond length in the reference and current configuration is, respectively, . Let denotes the displacement of two nodes. For a single bond, the nodal force is derived as

and the stiffness matrix is

Figure 3: A single bond in the reference and the current configuration.

The viscous matrix is

The derivatives of bond quantities are

5. 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 different meshing size schemes. The simulation specimen is shown in Figure 4(a), which is subjected to a constant axial tension, namely, . The specimen is discretized with irregular triangle cells, as shown in Figure 4(b). The material parameters are and . The analytical solution of this case is

(a)

(b)

(c)

Figure 4: Convergence check of the mechanical model: (a) simulation specimen; (b) meshing scheme; (c) comparison between the analytical and the simulated strain evolution with time (the size means the ratio of the average element area to the specimen area: , , , and ).

The elastic potential function used for this simulation is

where is the microbond parameter, calibrated as for the 3D case, for the 2D in-plane stress case, and for the 2D in-plane strain case.

The 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 the analytical results. Eventually, the simulation results almost overlap the analytical ones. This suggests that the present model is stable and convergent.

5.2. 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]. The 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. The material density of the bar is , Young’s modulus , Poisson’s ratio , and viscosity coefficient . The bar size is 50 mm  50 mm  1000 mm, whose aspect ratio can ensure the one-dimensional wave propagation. In the numerical simulation, different impact speeds, namely, , , and , are adopted. Here, we take the two-parameter cohesive force law as the hyperelastic bond potential, which is

where and are the microparameters and is calibrated through (3). is a scale parameter, defined as , with being the strain value at the peak stress of the uniaxial stress-strain curve.

(a)

(b)

(c)

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.

The 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 , which completely reduces to the hyperelastic situation. The 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. This 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). This example suggests that the present method can reflect the viscosity mechanism to great extent.

5.3. 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 and Young’s modulus , as shown in Figure 6(a). The dimension of the specimen is , and the crack inclination is . Take the two-parameter cohesive force law (30) as the hyperelastic bond potential.

(a)

(b)

(c)

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

To examine the effect of viscosity coefficient on the creep behaviors, we take , , and , respectively, and the compressive stress . The 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. This agrees with the general regularity of creep.

To examine the effect of load magnitude, we take , , , and , respectively, and . The simulated strain is shown in Figure 6(c). From Figure 6(c), it is seen that when the stress is small, say , the strain does not evolve further when it reaches a certain level. However, with the increase of stress, the strain evolves faster and faster. This also agrees with the typical features of creep. The simulated fracturing process with and is shown in Figure 7, which demonstrates that the specimen fails in the fracture mode.

(a)

(b)

(c)

Figure 7: Simulated creep fracturing process: (a) t = 4000 s; (b) t = 10500 s; (c) t = 11500 s.

6. 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. The 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. The 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.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (no. 11772190), which is gratefully acknowledged.

References

N. Challamel, C. Lanos, and C. Casandjian, “Creep damage modelling for quasi-brittle materials,” European Journal of Mechanics–A/Solids, vol. 24, no. 4, pp. 593–613, 2005.

I. G. Main, “A damage mechanics model for power-law creep and earthquake aftershock and foreshock sequences,” Geophysical Journal International, vol. 142, no. 1, pp. 151–161, 2000.

V. De Gennaro and J. Pereira, “A viscoplastic constitutive model for unsaturated geomaterials,” Computers and Geotechnics, vol. 54, pp. 143–151, 2013.

J. Shao, Q. Zhu, and K. Su, “Modeling of creep in rock materials in terms of material degradation,” Computers and Geotechnics, vol. 30, no. 7, pp. 549–555, 2003.

J. Shao, G. Duveau, N. Hoteit, M. Sibai, and M. Bart, “Time dependent continuous damage model for deformation and failure of brittle rock,” International Journal of Rock Mechanics and Mining Sciences, vol. 34, no. 3, pp. 285.e1–e285.e13, 1997.

H. Yoshida and H. Horii, “A micromechanics-based model for creep behavior of rock,” Applied Mechanics Reviews, vol. 45, no. 8, pp. 294–303, 1992.

H. Horii and S. Nemat-Nasser, “Compression-induced micro crack growth in brittle solids: axial splitting and shear failure,” Journal of Geophysical Research, vol. 90, no. B4, pp. 3105–3125, 1985.

A. Golshani, M. Oda, Y. Okui, T. Takemura, and E. Munkhtogoo, “Numerical simulation of the excavation damaged zone around an opening in brittle rock,” International Journal of Rock Mechanics and Mining Sciences, vol. 44, no. 6, pp. 835–845, 2007.

H. Konietzky, A. Heftenberger, and M. Feige, “Life-time prediction for rocks under static compressive and tensile loads: a new simulation approach,” Acta Geotechnica, vol. 4, no. 1, pp. 73–78, 2009.

N. Brantut, P. Baud, M. Heap, and P. Meredith, “Micromechanics of brittle creep in rocks,” Journal of Geophysical Research: Solid Earth, vol. 117, no. B8, 2012.

T. Xu, C. Tang, J. Zhao, L. Li, and M. Heap, “Modelling the time-dependent rheological behaviour of heterogeneous brittle rocks,” Geophysical Journal International, vol. 189, no. 3, pp. 1781–1796, 2012.

Y. Lu, D. Elsworth, and L. Wang, “A dual-scale approach to model time-dependent deformation, creep and fracturing of brittle rocks,” Computers and Geotechnics, vol. 60, pp. 61–76, 2014.

A. Hrennikoff, “Solution of problems of elasticity by the framework method,” Journal of Applied Mechanics, vol. 8, no. 4, pp. A169–A175, 1941.

G. Lilliu and J. G. M. Van Mier, “3D lattice type of fracture model for concrete,” Engineering Fracture Mechanics, vol. 70, pp. 927–941, 2003.

S. Kale and M. Ostoja-Starzewski, “Elastic-plastic-brittle transitions and avalanches in disordered media,” Physical Review Letters, vol. 112, no. 4, p. 045503, 2014.

M. Jirásek and Z. Bazant, “Macroscopic fracture characteristics of random particle systems,” International Journal of Fracture, vol. 69, no. 3, pp. 201–228, 1995.

G. Cusatis, Z. Bazant, and L. Cedolin, “Confinement-shear lattice model for concrete damage in tension and compression: I. Theory,” Journal of Engineering Mechanics, vol. 129, no. 12, pp. 1439–1448, 2003.

G. Wang, A. Al-Ostaz, A. H. D. Cheng, and P. R. Mantena, “Hybrid lattice particle modeling: theoretical considerations for a 2D elastic spring network for dynamic fracture simulations,” Computational Materials Science, vol. 44, no. 4, pp. 1126–1134, 2009.

S. F. Zhao and G. F. Zhao, “Implementation of a high order lattice spring model for elasticity,” International Journal of Solids and Structures, vol. 49, no. 18, pp. 2568–2581, 2012.

G. F. Zhao, “Developing a four-dimensional lattice spring model for mechanical responses of solids,” Computer Methods in Applied Mechanics and Engineering, vol. 315, pp. 881–895, 2017.

P. Grassl and M. Jirasek, “Meso-scale approach to modelling the fracture process zone of concrete subjected to uniaxial tension,” International Journal of Solids and Structures, vol. 47, no. 7-8, pp. 957–968, 2010.

Z. Zhang, “Discretized virtual internal bond model for nonlinear elasticity,” International Journal of Solids and Structures, vol. 50, no. 22-23, pp. 3618–3625, 2013.

Z. Zhang and Y. Chen, “Modeling nonlinear elastic solid with correlated lattice bond cell for dynamic fracture simulation,” Computer Methods in Applied Mechanics and Engineering, vol. 279, pp. 325–347, 2014.

Z. Zhang, Y. Yao, and X. B. Mao, “Modeling wave propagation induced fracture in rock with correlated lattice bond cell,” International Journal of Rock Mechanics and Mining Sciences, vol. 78, pp. 262–270, 2015.

Z. Zhang, J. F. Ding, A. Ghassemi, and X. R. Ge, “A hyperelastic-bilinear potential for lattice model with fracture energy conservation,” Engineering Fracture Mechanics, vol. 142, pp. 220–235, 2015.

L. Niu and W. Zhu, “Experimental investigation on loading rate effect of tensile strength influenced by sandstone viscosity,” Chinese Journal of Rock Mechanics and Rock Engineering, vol. 10, pp. 2466–2473, 2017.