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 iswhere is the deformed bond length and is the bond number in a cell, defined as , with being the node number of a cell.
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 aswhere 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 aswhere 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.
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 aswhere 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 aswhere 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 aswhere 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 withwhere 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 isand 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 asand the stiffness matrix is
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)
The elastic potential function used for this simulation iswhere 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 iswhere 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)
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)
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)
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.