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

Authors: Wangyang He; Zhennan Zhang
Journal: Advances in Civil Engineering
 (2018)

---

## 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.


---

## Body

## 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 Φ=Φl. So, the total strain energy of the bond cell is(1)W=∑ΩΦl,where l is the deformed bond length and Ω is the bond number in a cell, defined as Ω=NN−1/2, with N 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:(2)Fi=∂W∂ui=∑∂Φl∂l⋅∂l∂ui,Kij=∂2W∂ui∂uj=∑∂2Φl∂ui∂uj=∑Φ″l⋅∂l∂ui⋅∂l∂uj+Φ′l⋅∂2l∂ui∂uj,where ui is the displacement component.The physical parameter of bond potential is calibrated as(3)Φ″l0=λEVΩl02,where E is Young’s modulus; V is the cell volume; l0 is the undeformed bond length; and λ is a coefficient, with λ=6, λ=3, and λ=16/5 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 Figure2. 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(4)ψ=Φl+η˜l˙l−l0,where Φl is the hyperelastic potential of the bond; l˙ is the bond deformation ratio, l˙=∂l/∂t; and η˜ is the microviscosity coefficient, defined as η˜=Fc/l˙, with Fc 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(5)Fi=∑∂ψ∂ui=∑Φ′l⋅∂l∂ui+η˜l˙⋅∂l∂ui+η˜l−l0⋅∂l˙∂ui.The tangent stiffness and viscous matrix are, respectively, derived as(6)Kij=∂Fi∂uj=∑Φ″l⋅∂l∂ui⋅∂l∂uj+Φ′l⋅∂2l∂ui∂uj+η˜l˙⋅∂2l∂ui∂uj+η˜⋅∂l∂ui⋅∂l˙∂uj+η˜l−l0⋅∂2l˙∂ui∂uj+η˜⋅∂l∂uj⋅∂l˙∂ui,Cij=∂Fi∂u˙j=∑η˜⋅∂l˙∂u˙j⋅∂l∂ui+η˜l−l0⋅∂2l˙∂ui∂u˙j.


### 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 Φ″l0=6EV/Ωl02 for the 3D cases and Φ″l0=3EV/Ωl02 for the in-plane stress case and Φ″l0=16EV/5Ωl02 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(7)W=1V∑ψ=1V∑Φl+η˜l˙l−l0=1VΦl+η˜l˙l−l0,where the operator ⋯ denotes ⋯=∫02π∫0π⋯Ω/4π⁡sin⁡θ⁡dθ⁡dϕ for the 3D cases and ⋯=∫02π⋯Ω/2πdθ for the 2D cases. The bond length and its ratio in terms of strain tensor are written as(8)l=l0+l0εijξiξj,l˙=l0ε˙ijξiξj.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(9)σij=∂W∂εij=1VΦ′l∂l∂εij+η˜∂l˙∂εijl−l0+η˜l˙∂l∂εij.The tangent modulus tensor is(10)Kijkl=∂σij∂εkl=1VΦ″l∂l∂εkl⋅∂l∂εij+Φ′l∂2l∂εij∂εkl+η˜∂2l˙∂εij∂εkll−l0+η˜∂l˙∂εij⋅∂l∂εkl+η˜∂l˙∂εkl⋅∂l∂εij+η˜l˙∂2l∂εij∂εkl.The tangent viscosity tensor is(11)Cijkl=∂σij∂ε˙kl=1VΦ″l∂l∂εkl⋅∂l∂ε˙ij+Φ′l∂2l∂εij∂ε˙kl+η˜∂2l˙∂εij∂ε˙kll−l0+η˜∂l˙∂εij⋅∂l∂ε˙kl+η˜∂l˙∂ε˙kl⋅∂l∂εij+η˜l˙∂2l∂εij∂ε˙kl.Substituting (8) into (10) and (11), we can obtain the tangent stiffness tensor and viscosity tensor under the undeformed state as(12)Kijkl=1VΦ″l0l02ξiξjξkξl,Cijkl=1Vη˜l02ξiξjξkξl.In the calibration of the initial bond stiffnessΦ″l0, Poisson’s ratio is fixed, namely, ν=0.25. Analogous to the initial bond stiffness calibration, we define the macroviscosity coefficient η as(13)σ=ηε˙.To relate the microviscosity to the macro one, we expand (13) to the 3D form, analogous to the elastic matrix with ν=0.25:(14)ΩC=25η311000131000113000000100000010000001,where ΩC complies with the relationship(15)σ11,σ22,σ33,τ12,τ23,τ31T=ΩCε˙11,ε˙22,ε˙33,γ˙12,γ˙23,γ˙31T.According to the tangent viscosity tensor, the viscosity matrix obtained by integrating (12) is(16)Ωb=Ωl0215V3η˜η˜η˜000η˜3η˜η˜000η˜η˜3η˜000000η˜000000η˜000000η˜.By equaling (14) to (16), the microviscosity coefficient is calibrated as(17)Ωl0215Vη˜=25η⁡⇒⁡η˜=6ηVΩl02.By the same calibration process, the parameters in 2D case are calibrated as(18)η˜=3ηVΩl02for  plane  stress16ηV5Ωl02for  plane  stress.


## 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 Figure2. 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(4)ψ=Φl+η˜l˙l−l0,where Φl is the hyperelastic potential of the bond; l˙ is the bond deformation ratio, l˙=∂l/∂t; and η˜ is the microviscosity coefficient, defined as η˜=Fc/l˙, with Fc 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(5)Fi=∑∂ψ∂ui=∑Φ′l⋅∂l∂ui+η˜l˙⋅∂l∂ui+η˜l−l0⋅∂l˙∂ui.The tangent stiffness and viscous matrix are, respectively, derived as(6)Kij=∂Fi∂uj=∑Φ″l⋅∂l∂ui⋅∂l∂uj+Φ′l⋅∂2l∂ui∂uj+η˜l˙⋅∂2l∂ui∂uj+η˜⋅∂l∂ui⋅∂l˙∂uj+η˜l−l0⋅∂2l˙∂ui∂uj+η˜⋅∂l∂uj⋅∂l˙∂ui,Cij=∂Fi∂u˙j=∑η˜⋅∂l˙∂u˙j⋅∂l∂ui+η˜l−l0⋅∂2l˙∂ui∂u˙j.


## 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 Φ″l0=6EV/Ωl02 for the 3D cases and Φ″l0=3EV/Ωl02 for the in-plane stress case and Φ″l0=16EV/5Ωl02 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(7)W=1V∑ψ=1V∑Φl+η˜l˙l−l0=1VΦl+η˜l˙l−l0,where the operator ⋯ denotes ⋯=∫02π∫0π⋯Ω/4π⁡sin⁡θ⁡dθ⁡dϕ for the 3D cases and ⋯=∫02π⋯Ω/2πdθ for the 2D cases. The bond length and its ratio in terms of strain tensor are written as(8)l=l0+l0εijξiξj,l˙=l0ε˙ijξiξj.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(9)σij=∂W∂εij=1VΦ′l∂l∂εij+η˜∂l˙∂εijl−l0+η˜l˙∂l∂εij.The tangent modulus tensor is(10)Kijkl=∂σij∂εkl=1VΦ″l∂l∂εkl⋅∂l∂εij+Φ′l∂2l∂εij∂εkl+η˜∂2l˙∂εij∂εkll−l0+η˜∂l˙∂εij⋅∂l∂εkl+η˜∂l˙∂εkl⋅∂l∂εij+η˜l˙∂2l∂εij∂εkl.The tangent viscosity tensor is(11)Cijkl=∂σij∂ε˙kl=1VΦ″l∂l∂εkl⋅∂l∂ε˙ij+Φ′l∂2l∂εij∂ε˙kl+η˜∂2l˙∂εij∂ε˙kll−l0+η˜∂l˙∂εij⋅∂l∂ε˙kl+η˜∂l˙∂ε˙kl⋅∂l∂εij+η˜l˙∂2l∂εij∂ε˙kl.Substituting (8) into (10) and (11), we can obtain the tangent stiffness tensor and viscosity tensor under the undeformed state as(12)Kijkl=1VΦ″l0l02ξiξjξkξl,Cijkl=1Vη˜l02ξiξjξkξl.In the calibration of the initial bond stiffnessΦ″l0, Poisson’s ratio is fixed, namely, ν=0.25. Analogous to the initial bond stiffness calibration, we define the macroviscosity coefficient η as(13)σ=ηε˙.To relate the microviscosity to the macro one, we expand (13) to the 3D form, analogous to the elastic matrix with ν=0.25:(14)ΩC=25η311000131000113000000100000010000001,where ΩC complies with the relationship(15)σ11,σ22,σ33,τ12,τ23,τ31T=ΩCε˙11,ε˙22,ε˙33,γ˙12,γ˙23,γ˙31T.According to the tangent viscosity tensor, the viscosity matrix obtained by integrating (12) is(16)Ωb=Ωl0215V3η˜η˜η˜000η˜3η˜η˜000η˜η˜3η˜000000η˜000000η˜000000η˜.By equaling (14) to (16), the microviscosity coefficient is calibrated as(17)Ωl0215Vη˜=25η⁡⇒⁡η˜=6ηVΩl02.By the same calibration process, the parameters in 2D case are calibrated as(18)η˜=3ηVΩl02for  plane  stress16ηV5Ωl02for  plane  stress.


## 4. Numerical Implementation

### 4.1. Governing Equation

The general governing equation for the creep problem can be written as(19)Mu¨+Fu,u˙=R,where the restoring nodal force Fu,u˙ 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(20)u1=u0+Δt1−θu˙0+θu˙1,u˙1=u˙0+Δt1−θu¨0+θu¨1,where the subscript “0” means the current time step, while “1” means the next time step; Δt means the time interval; and θ is a coefficient. 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(21)1θ2Δt2Mu1+F1u1,u˙1=R1+1θ−1R0+1θ2Δt2Mu0+1θ2ΔtMu˙0⁡−⁡1θ−1F0u0,u˙0.To solve (21), the Newton–Raphson method is adopted, whose iteration algorithm is(22)1θ2Δt2M+∂F1u1i,u˙1i∂u˙1i⋅1θΔt+∂F1u1i,u˙1i∂u1iΔu1i+1=R1+1θ−1R0+1θ2Δt2Mu0+1θ2ΔtMu˙0⁡−⁡1θ−1F0u0,u˙0−1θ2Δt2Mu1i+F1u1i,u˙1i,u1i+1=u1i+Δu1i+1,u˙1i+1=1θΔtu1i+1−1θΔtu0−1θ−1u˙0,and the quantities are updated at the next time step as(23)u˙1=1θΔtu1⁡−⁡1θΔtu0⁡−⁡1θ−1u˙0,u¨1=1θ2Δt2u1−1θ2Δt2u0−1θ2Δtu˙0⁡−⁡1θ−1u¨0.


### 4.2. Computation of Bond Quantities

Take the single bond shown in Figure3 as an example to show the computation of the stiffness matrix and viscous matrix. Let X¯,Y¯ denote the bond vector in the reference and current configuration, respectively, X¯=XII−XI, Y¯=YII−YI, with XI,YI and XII,YII being the coordinates of the node I and node II in the reference and current configuration, X¯=x¯1,x¯2,x¯3, Y¯=y¯1,y¯2,y¯3. Let uI,uII denote the displacement vector of the node I and node II. Then, X¯,Y¯ are correlated by Y¯=X¯+uII−uI. The bond length in the reference and current configuration is, respectively, l0=X¯TX¯,l=Y¯TY¯. Let u=u1,u2,u3,u4,u5,u6T denotes the displacement of two nodes. For a single bond, the nodal force is derived as(24)Fib=Φ′l⋅∂l∂ui+η˜l˙⋅∂l∂ui+η˜l−l0⋅∂l˙∂ui,and the stiffness matrix is(25)Kijb=Φ″l⋅∂l∂ui⋅∂l∂uj+Φ′l+η˜l˙⋅l0l∂2l∂ui∂uj+η˜⋅l0l∂l∂ui⋅∂l˙∂uj+∂l∂uj⋅∂l˙∂ui.Figure 3
A single bond in the reference and the current configuration.The viscous matrix is(26)Cijb=η˜l−l0lpij+l0l⋅∂l∂ui⋅∂l˙∂u˙j.The derivatives of bond quantities are(27)l02=x¯12+x¯22+x¯32,l2=y¯12+y¯22+y¯32=x¯1+u4−u12+x¯2+u5−u22+x¯3+u6−u32,l˙=y¯1lu˙4−u˙1+y¯2lu˙5−u˙2+y¯3lu˙6−u˙3,∂l∂ui=y¯1lδ4i−δ1i+y¯2lδ5i−δ2i+y¯3lδ6i−δ3i,∂l˙∂ui=1lu˙4−u˙1δ4i−δ1i+1lu˙5−u˙2δ5i−δ2i+1lu˙6−u˙3δ6i−δ3i−l˙l⋅∂l∂ui,∂l˙∂u˙j=∂l∂ui,∂2l∂ui∂uj=1lpij−1l⋅∂l∂uj⋅∂l∂ui,∂2l˙∂ui∂uj=−1l⋅∂l∂ui⋅∂l˙∂uj−1l⋅∂l∂uj⋅∂l˙∂ui−l˙l⋅∂2l∂ui∂uj,pij=δ4j−δ1jδ4i−δ1i+δ5j−δ2jδ5i−δ2i+δ6j−δ3jδ6i−δ3i,∂2l˙∂ui∂u˙j=1lpij−1l⋅∂l∂ui⋅∂l˙∂u˙j,∂l∂ui⋅∂l˙∂uj+∂l∂uj⋅∂l˙∂ui+l⋅∂2l˙∂ui∂uj+l˙⋅∂2l∂ui∂uj=0.


## 4.1. Governing Equation

The general governing equation for the creep problem can be written as(19)Mu¨+Fu,u˙=R,where the restoring nodal force Fu,u˙ 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(20)u1=u0+Δt1−θu˙0+θu˙1,u˙1=u˙0+Δt1−θu¨0+θu¨1,where the subscript “0” means the current time step, while “1” means the next time step; Δt means the time interval; and θ is a coefficient. 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(21)1θ2Δt2Mu1+F1u1,u˙1=R1+1θ−1R0+1θ2Δt2Mu0+1θ2ΔtMu˙0⁡−⁡1θ−1F0u0,u˙0.To solve (21), the Newton–Raphson method is adopted, whose iteration algorithm is(22)1θ2Δt2M+∂F1u1i,u˙1i∂u˙1i⋅1θΔt+∂F1u1i,u˙1i∂u1iΔu1i+1=R1+1θ−1R0+1θ2Δt2Mu0+1θ2ΔtMu˙0⁡−⁡1θ−1F0u0,u˙0−1θ2Δt2Mu1i+F1u1i,u˙1i,u1i+1=u1i+Δu1i+1,u˙1i+1=1θΔtu1i+1−1θΔtu0−1θ−1u˙0,and the quantities are updated at the next time step as(23)u˙1=1θΔtu1⁡−⁡1θΔtu0⁡−⁡1θ−1u˙0,u¨1=1θ2Δt2u1−1θ2Δt2u0−1θ2Δtu˙0⁡−⁡1θ−1u¨0.


## 4.2. Computation of Bond Quantities

Take the single bond shown in Figure3 as an example to show the computation of the stiffness matrix and viscous matrix. Let X¯,Y¯ denote the bond vector in the reference and current configuration, respectively, X¯=XII−XI, Y¯=YII−YI, with XI,YI and XII,YII being the coordinates of the node I and node II in the reference and current configuration, X¯=x¯1,x¯2,x¯3, Y¯=y¯1,y¯2,y¯3. Let uI,uII denote the displacement vector of the node I and node II. Then, X¯,Y¯ are correlated by Y¯=X¯+uII−uI. The bond length in the reference and current configuration is, respectively, l0=X¯TX¯,l=Y¯TY¯. Let u=u1,u2,u3,u4,u5,u6T denotes the displacement of two nodes. For a single bond, the nodal force is derived as(24)Fib=Φ′l⋅∂l∂ui+η˜l˙⋅∂l∂ui+η˜l−l0⋅∂l˙∂ui,and the stiffness matrix is(25)Kijb=Φ″l⋅∂l∂ui⋅∂l∂uj+Φ′l+η˜l˙⋅l0l∂2l∂ui∂uj+η˜⋅l0l∂l∂ui⋅∂l˙∂uj+∂l∂uj⋅∂l˙∂ui.Figure 3
A single bond in the reference and the current configuration.The viscous matrix is(26)Cijb=η˜l−l0lpij+l0l⋅∂l∂ui⋅∂l˙∂u˙j.The derivatives of bond quantities are(27)l02=x¯12+x¯22+x¯32,l2=y¯12+y¯22+y¯32=x¯1+u4−u12+x¯2+u5−u22+x¯3+u6−u32,l˙=y¯1lu˙4−u˙1+y¯2lu˙5−u˙2+y¯3lu˙6−u˙3,∂l∂ui=y¯1lδ4i−δ1i+y¯2lδ5i−δ2i+y¯3lδ6i−δ3i,∂l˙∂ui=1lu˙4−u˙1δ4i−δ1i+1lu˙5−u˙2δ5i−δ2i+1lu˙6−u˙3δ6i−δ3i−l˙l⋅∂l∂ui,∂l˙∂u˙j=∂l∂ui,∂2l∂ui∂uj=1lpij−1l⋅∂l∂uj⋅∂l∂ui,∂2l˙∂ui∂uj=−1l⋅∂l∂ui⋅∂l˙∂uj−1l⋅∂l∂uj⋅∂l˙∂ui−l˙l⋅∂2l∂ui∂uj,pij=δ4j−δ1jδ4i−δ1i+δ5j−δ2jδ5i−δ2i+δ6j−δ3jδ6i−δ3i,∂2l˙∂ui∂u˙j=1lpij−1l⋅∂l∂ui⋅∂l˙∂u˙j,∂l∂ui⋅∂l˙∂uj+∂l∂uj⋅∂l˙∂ui+l⋅∂2l˙∂ui∂uj+l˙⋅∂2l∂ui∂uj=0.


## 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 Figure4(a), which is subjected to a constant axial tension, namely, σ0=5⁡MPa. The specimen is discretized with irregular triangle cells, as shown in Figure 4(b). The material parameters are E=40⁡GPa and η=1000⁡GPa⋅s. The analytical solution of this case is(28)ε=σ0E1−exp−Eηt.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:size1=2.4×10−3, size2=7.1×10−4, size3=1.7×10−4, and size4=8.9×10−5).
(a)
(b)
(c)The elastic potential function used for this simulation is(29)Φl=12Al−l02,where A is the microbond parameter, calibrated as A=6EV/Ωl02 for the 3D case, A=3EV/Ωl02 for the 2D in-plane stress case, and A=16EV/5Ωl02 for the 2D in-plane strain case.The comparison between the analytical and simulated strain is shown in Figure4(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 ρ=2168.5⁡kg/m3, Young’s modulus E=8.4⁡Gpa, Poisson’s ratio v=0.29, and viscosity coefficient η=100⁡KPa⋅s. 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, V1=1.3⁡m/s, V2=2.0⁡m/s, and V3=2.5⁡m/s, are adopted. Here, we take the two-parameter cohesive force law as the hyperelastic bond potential, which is(30)Φ′l=Al−l0⁡exp−l−l0B,where A and B are the microparameters and A is calibrated through (3). B is a scale parameter, defined as B=εtl0, with εt being the strain value at the peak stress of the uniaxial stress-strain curve.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.
(a)
(b)
(c)The simulated results are shown in Figure5(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. 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σ0 and Young’s modulus E=40⁡GPa, as shown in Figure 6(a). The dimension of the specimen is 1.0⁡m×1.0⁡m, and the crack inclination is 45°. Take the two-parameter cohesive force law (30) as the hyperelastic bond potential.Figure 6
Simulation of creep fracture: (a) simulation specimen; (b) effect of viscosity on creep; (c) effect of external stress on creep.
(a)
(b)
(c)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. 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σ0=5.0⁡MPa, σ0=6.5⁡MPa, σ0=7.0⁡MPa, and σ0=7.5⁡MPa, respectively, and η=1000⁡GPa⋅s. The 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. This also agrees with the typical features of creep. The 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.Figure 7
Simulated creep fracturing process: (a)t = 4000 s; (b) t = 10500 s; (c) t = 11500 s.
(a)
(b)
(c)


## 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 Figure4(a), which is subjected to a constant axial tension, namely, σ0=5⁡MPa. The specimen is discretized with irregular triangle cells, as shown in Figure 4(b). The material parameters are E=40⁡GPa and η=1000⁡GPa⋅s. The analytical solution of this case is(28)ε=σ0E1−exp−Eηt.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:size1=2.4×10−3, size2=7.1×10−4, size3=1.7×10−4, and size4=8.9×10−5).
(a)
(b)
(c)The elastic potential function used for this simulation is(29)Φl=12Al−l02,where A is the microbond parameter, calibrated as A=6EV/Ωl02 for the 3D case, A=3EV/Ωl02 for the 2D in-plane stress case, and A=16EV/5Ωl02 for the 2D in-plane strain case.The comparison between the analytical and simulated strain is shown in Figure4(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 ρ=2168.5⁡kg/m3, Young’s modulus E=8.4⁡Gpa, Poisson’s ratio v=0.29, and viscosity coefficient η=100⁡KPa⋅s. 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, V1=1.3⁡m/s, V2=2.0⁡m/s, and V3=2.5⁡m/s, are adopted. Here, we take the two-parameter cohesive force law as the hyperelastic bond potential, which is(30)Φ′l=Al−l0⁡exp−l−l0B,where A and B are the microparameters and A is calibrated through (3). B is a scale parameter, defined as B=εtl0, with εt being the strain value at the peak stress of the uniaxial stress-strain curve.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.
(a)
(b)
(c)The simulated results are shown in Figure5(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. 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σ0 and Young’s modulus E=40⁡GPa, as shown in Figure 6(a). The dimension of the specimen is 1.0⁡m×1.0⁡m, and the crack inclination is 45°. Take the two-parameter cohesive force law (30) as the hyperelastic bond potential.Figure 6
Simulation of creep fracture: (a) simulation specimen; (b) effect of viscosity on creep; (c) effect of external stress on creep.
(a)
(b)
(c)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. 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σ0=5.0⁡MPa, σ0=6.5⁡MPa, σ0=7.0⁡MPa, and σ0=7.5⁡MPa, respectively, and η=1000⁡GPa⋅s. The 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. This also agrees with the typical features of creep. The 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.Figure 7
Simulated creep fracturing process: (a)t = 4000 s; (b) t = 10500 s; (c) t = 11500 s.
(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.



---

*Source: 8042965-2018-04-05.xml*