# A Coupled Damage-Permeability Constitutive Model for Brittle Rocks Subjected to Explosive Loading

Authors: Qizhi Wang; Xuebin Su; Bangbiao Wu; Wei Wang; Wei Yuan
Journal: Advances in Civil Engineering
 (2018)

---

## Abstract

A dynamic constitutive model considering coupled damage-permeability effect is developed for rocks under explosion. The main purpose of this model is to provide reference to in situ leaching mining (include uranium ore, gold ore, etc.) to improve the permeability of rocks using blasting. The damage is induced and described by the summation of all principal tensile strains. This method does not need to determine (distinguish) the stress state of the material and thus does not need to consider the damage induced by differential stresses. A corresponding relation between the damage status and rock permeability is established through a logarithm function. Based on the theoretical foundation, this model is embedded in a three-dimensional code of ABAQUS, and its applicability is verified through simulating the experiments under the condition of explosive load, the results from numerical simulation are consistent with the experimental results. Moreover, the obvious advantage of the numerical constitutive model is that it can clearly distinguish permeability variation that cannot be done in laboratory test.


---

## Body

## 1. Introduction

In situ leaching mining of deeply buried minerals depends on the permeability of rocks. Methods to improve permeability usually generate a damage zone in which geotechnical and hydrogeological properties are deteriorated. This degradation process, which is due to microcracks growth, is generally accompanied by significant changes in flow and permeability properties. In the current solid metal industry, the increase of the permeability of reservoir rocks can be beneficial to the mining process and lead to economic efficiency. Different from the petroleum industry where large fractures are desirable for oil or gas production, the metal leaching method aims to increase the overall number of microcracks to expand chemical solution contact area by enhancing the general permeability for the purpose of recovery efficiency. Therefore, it is worth studying how to increase the permeability of reservoir rocks for mineral production. The main purpose of this article is to develop a damage-permeability model under explosion and to obtain the permeability-damage regularity from the numerical simulation and provide guidance for blasting practices.Different brittle rock dynamic damage models under explosion loading have been proposed and developed by Grady and Kipp [1], Kipp and Grady [2], Taylor et al. [3], Thorne [4, 5], Liu and Katsabanis [6], and Chen [7], respectively. These models constructed the damage criterion by means of the tensile damage. Based on mesoscopic damage mechanics and the framework of continuum damage mechanics, Zhou et al. [8] and Souley et al. [9] established the damage model coupled with permeability under static loadings. In addition, some phenomenological damage models were developed to explain the permeability evolution of brittle rocks mainly characterized by crack growth [10], density distribution [11], connectivity, and the relation to elastic wave velocities [12–15]. But few models have been reported considering the damage-permeability coupling effect under explosive load.Based on the work by Thorne [4, 5] and Souley et al. [9], the objective of the present study is to build a damage-permeability coupling model for rocks subjected to explosive loading. This model depicts the relationship between dynamic damage and permeability. In order to verify the reliability of the model, it is embedded in the three-dimensional code (ABAQUS) and used to simulate the borehole blasting. Comparison of the numerical results with experimental results validates the effectiveness of the model.


## 2. Damage Variable and Its Evolution Law

Different from existing damage criterion based on the volume strain, this model adopts the summation of the tensile principal strains as the criterion of damage accumulation, reflecting the material damage caused by differential stress. This model simplifies the description of the problem to make its physical meaning clearer and assumes that the permeability is isotropic; thus, the corresponding relation between permeability changes and the damage status of rocks is established.On the construction of the damage model, the summation of principal tensile strainθ and damage threshold θc is the basis that determines the nucleation and propagation of the damage. If the summation of the principal tensile strain of a point exceeds θc, the microcracks start to extend (i.e., the damage began to produce), with the decrease of the material stiffness. On the contrary, when θ<θc, the elastic-plastic model will be used in the simulation of the material dynamic response under explosion load.For simplicity, it is assumed that the modeled material is isotropic with randomly distributed cracks, which represents the damage accumulation and that the Poisson’s ratio is constant under explosive loading. Crack density can be expressed by the average ruler of microcracks as(1)Cd=γNa3,whereN is the number of activated cracks, γ is a proportionality constant, and a is the average size of microcracks. N is related to the tensile volumetric strain and obeys Weibull distribution [2].In the process of constructing the damage model in this paper,θ represents the sum of the principal tensile strain, and it is assumed that the propagation of microcracks starts when θ reaches its critical threshold θc, then number of activated cracks can be written as(2)N=kθ−θcm,(3)θ=∑εit,(4)θc=TcE,where θ is the sum of principal tensile strain, m and k are both material constants, εit=εi, εi is tensile strain, Tc is dynamic tensile strength, E is elasticity modulus, and Macaulay bracket  is the operator.(5)L=h=LL,where h=L is a step function, and h=L=0 when L≤0, h=L=1 when L>0.After getting the crack density, the influence of the crack density on the stiffness is obtained by statistical fracture mechanics. Denote the crack density asCd and the unit volume as V0 of the material, the probability for it not to contain a crack is(6)ps=e−CdV0.If there is one microcrack initiation in the unit volumeV0, the probability of the material to crack can be written as(7)pf=1−e−CdV0.In order to consider the effect of overlapping between microcracks, Englman and Jaeger [16] introduced the expression(8)F=1−e−αCd,where α=16/9. In this paper, a scalar damage variable D is introduced as(9)D=F=1−e−αCd.Differentiating (9) with time leads to(10)dDdt=αe−αCddCddt.Differentiating (1) with time gets(11)dCddt=γdNdta3+3γNa2dadt.Differentiating (2) with time gets(12)dNdt=mkθ−θcm−1dθdt.Substituting (11) into (12),(13)dCddt=mka2θ−θcm−1∗aγdθdt+3γθ−θcmdadt.Considering the change of strain rate, an expression for the average size of microcracks is proposed by Taylor el al. [3] and Kuszmaul [17]:(14)a=1220KICρcRmax2/3,where KIC is the fracture toughness, ρ is the density, and c is the P-wave velocity of the material, respectively. Rmax represents the experienced largest strain rate of the material. Assume a is time independent:(15)dadt=0.Combining (1), (2), and (14) gives(16)Cd=5γkθ−θcm2KICρcRmax2.Substituting (15) and (16) into (13) leads to(17)dCddt=5mrkθ−θcm−12KICρcRmax2dθdt.According to the study by Thorne [4], parameter m = 6. A damage evolution law is thus established.


## 3. Dynamic Constitutive Relation

Leopold [18] pointed out that the change of Poisson’s ratio of a material depends on the shape of the micro racks. When the material contains opening microcracks, the effective Poisson’s ratio is less than its inherent Poisson’s ratio. When the material contains a large number of closed microcracks, the effective Poisson’s ratio is greater than its inherent Poisson’s ratio.The developed model in this paper ignores the specific shapes and the opening or closing states of the microcracks and assumes that the probability of occurrence is the same for various types of microcracks. Therefore, it can be assumed that the Poisson’s ratio is constant during the damage process. According to the conclusion drawn by Yang et al. [19], the constant Poisson’s ratio can be used to maintain the constant relations among the elastic properties of the damaged material. Therefore, the incremental stress-strain relationship of the damaged material can be expressed as(18)dσij=2Gddεij+λddεijδij,(19)Gd=1−DE21+v,λd=1−DEν1+ν1−2ν,where E and v are elasticity modulus and Poisson’s ratio of the material, respectively.Microcracks begin to further extend only when the material is subjected to strain larger than previous maximum strain in the whole loading process, while the damage of material further accumulates. Therefore, in this model,θc is always substituted by the updated maximum value of the summation of the principal tensile strain θmax.By far, the stress-strain relationship is given in the case whenθ>θc; when θ≤θc, the elastoplastic servo hardening constitutive model is adopted to simulate the stress-strain relationship [20]:(20)φ=γk=12ξijξij−σy23=0,(21)ξij=sij−αij,σy=1+ε˙C11/pσ0,where ξij is a stress tensor, sij is the deviatoric stress tensor, αij is a tensor that represents the stress history, σ0 is the yield stress, C1 and p are material constants, and ε˙ is the strain rate.


## 4. Establishment of the Relationship between Permeability and Damage

Permeability of rock material is dominated by the rock mineral composition and the structure of the inner pore and fracture networks. When subjected to loadings, microcracks in rock would nucleate, grow, and coalesce, leading to the change of the permeability. Based on the comparison with the triaxial compression tests by Paterson [21], Souley et al. [22] made the following conclusions:The preliminary loadings will lead to the closure of existing initial pore and microcracks, resulting in the decrease of rock permeability. Then, the permeability nearly keeps constant in the early stable stage of crack extension because the hysteresis is obvious when the permeability is increasing. During the unstable crack propagation stage, the permeability increases substantially, which is reflected by the macroscopic fracture due to coalescence of the large number of microcracks. Although this conclusion is drawn from the static triaxial test, the nature of the damage associated with the increase of permeability is consistent, so the conclusion can be applied to the situation when subjected to dynamic load.Li [23]comprehensively summarized the equations showing the correlations between permeability and stress state. However, most of them were mainly focused on the study before material damage, which is inconvenient for engineering applications because the relationship is different for different materials and is also different of the same material under different stress states.This paper argues that the stress state of a material does not directly decide the permeability of the rock material, it is only one of the factors that induce the permeability change. Instead, the principal tensile strain can intuitively reflect the change of its permeability in some extent. Thus, based on the damage constitutive model by Souley et al. [9] revealing a logarithmic relationship between average microcrack size and permeability, it is assumed in this paper that the relationship between the change of the rock permeability and the corresponding damage status can be expressed as(22)logKK0=C2θθc−θsθc,θθc>θsθc,(23)K=K0,θθc≤θsθc,where parameter C2 is a material constant, K0 is the initial permeability of the material, and θs is defined as the critical value of strain which characterizes the hysteresis phenomenon of the permeability increase. The three parameters can be obtained from triaxial permeability tests. Therefore, the dynamic constitutive model of permeability-damage relations coupled with differential stress state is constructed by the closed (10), equations (17)–(20), (22), and (23).


## 5. Numerical Implementation

The damage model described above is implemented in the three-dimensional explicit finite program, ABAQUS. It is assumed that Poisson's ratio is constant under blasting and the material is isotropic. The main assumption associated with the damage tensor is that when the growth of microcracks occurs, the damage and stress tensors are assumed to be coaxial to the strain tensor, and both damage and stress tensors are adjusted according to the constitutive equations and the mechanical growth criteria.The implementation of the change in permeability induced by damage evolution is achieved as follows: (a) initial damage is computed by setting the initial crack density according to (9) and (10), and (b) with the continuous development of tensile strain, one or more subzones reach the threshold, the damage of material will be updated through (17) and (18), and the permeability coefficient required is then updated by (21) and (22) as shown in Figure 1.Figure 1
Flowchart of damage-permeability routine in ABAQUS mechanical calculations.


## 6. Comparison of the Modeling Results with Experimental Measurements

In order to verify the applicability of the numerical method, an explosion test is carried out to increase the permeability by blasting [24], using the typical five-borehole arrangement, as shown in Figure 2. Five-borehole layout can produce an obvious tensile stress region which is commonly used in engineering activities. After the ignition of the explosives, the propagation of shock wave results in instantaneous tangential tensile stress on the surrounding materials. When the dynamic stress is close to or greater than the rock strength, microcracks would nucleate and propagate. In the experiment, cement mortar is used to simulate the surrounding rock material [25] to increase the permeability by blasting, and the mortar pouring process is shown in Figure 3. Polyvinylidene fluoride (PVDF) pressure sensors are embedded around the boreholes to monitor the stress field and to establish the stress-time curve during the loading process. In situ transient pressure pulse testing method is used to measure the permeability, and the result is illustrated in Figure 4. The principle is to obtain the relationship of injected water pressure attenuation with time in the boreholes and then calculate the permeability from the attenuation characteristics, which is discussed by Wei et al. [26] in detail.Figure 2
Five blasting boreholes layout (cm).Figure 3
Implementation process of the experiment: mortar pouring and finished experimental setup.
(a)
(b)Figure 4
The explosive stress wave attenuation with distance curve (K = 1.79).The charging structure and detonating configuration are listed in Tables1 and 2. Figure 5 illustrates the experimental result of the permeability change, showing that the permeability is remarkably increased by two orders of magnitude after blasting (the initial permeability is 4.62 × 10−5 D) [16]. The permeability increasing distance is about 70 times of the charge radius. The experimental result also shows the relationship between the permeability and the decoupling coefficient: it is the same order of magnitude of the permeability when the decoupling coefficient is 1.79 and 2.57, respectively. But when the decoupling coefficient is greater than 3.29, the permeability decreases significantly, revealing the reduced effect of the blasting.Table 1
Physical-mechanical parameters and original permeabilities of the model.
Density (kg/m3)
Compressive strength (MPa)
Tensile strength (MPa)
Wave impedance (10−3/m2·s)
Poisson’s ratio
Permeability (D)
2140
21.2
1.36
6899
0.203
4.6 × 10−5Table 2
Charge structure and initiation sequence form.
Sample number
Charge (diameter × length) (cm)
Initiation sequence
K
1
1.4 × 24
1# and 2#, 3# and 4# blast holes millisecond blasting, respectively, 5# blast hole separate initiating
1.79
2
1.4 × 24
Single order initiation 1#-2#-3#-4#-5#
2.57
3
1.4 × 24
Single order initiation 1#-2#-3#-4#-5#
3.29
Sample numbers 1#, 2#, and 3# corresponding to 9#, 11#, and 12# in literature [24].Figure 5
Change of the permeability of rock with different decoupling coefficients.Three configurations of blasting design are used to compare the modeling results with the experimental ones. The only difference of the three blasting designs is the diameter of the blasting hole. The modeling parameters are as listed in Table3. Figure 6 illustrates a good agreement between the numerical predictions and the experiments on permeability change. The slight discrepancy might be due to the difference between the heterogeneity of the in situ blasting test and homogeneity of the numerical simulation model. However, the difference between the two situations is slight and the results obey the same regulations, which is acceptable.Table 3
Model parameter.
t
0
a
K
IC
0
MPa
m
v
T
c (MPa)
C
K
0 (D)
E (GPa)
θ
s
1.215
0.5
0.28
0.2
1.32
0.4
46.2
43
0.001
The parameters are derived from triaxial experiment and literature value [10, 22, 23, 27].Figure 6
Verification of numerical implementation with permeability measurements.Because of the hazards and limitations of blasting test, many experimental rules between parameters are hard to obtain, for example, the law of permeability changing with distance and the relationship between permeability variation range and crushing zone around boreholes. It is thus very necessary to develop suitable damage-permeability models under explosion on the basis of experiments.


## 7. Modeling Prediction of Damage and Permeability

Due to the limitation of testing methods, the relationship between crushing circle and permeability is difficult to obtain through experiments. However, numerical simulation can overcome the shortage of laboratory experiment using the damage-permeability model. In this section, the numerical model is carried out using the same size and configuration as the experimental ones in Figure2, and the material parameters are the same as in Tables 1 and 2.Figure7 illustrates the generation of the crushing area surrounding of boreholes with different decoupling coefficients, which are in a width range from 9 to 20 times of the radius. For example, in the blast hole 5#, the crushing circle range is about 13.33 cm (K = 1.79, cap radius = 0.7 cm), 8.65 cm (K = 2.57, cap radius = 0.7 cm), and 6.53 cm (K = 3.29, cap radius = 0.7 cm), respectively. The result shows that the smaller the value of K, the bigger the crushing zone. Although the crushing zone should be as small as possible to prevent the collapse, the filling material also absorbs energy which affects the transfer of energy to expand the permeability of surrounding rock. Therefore, each model has a reasonable borehole charge structure, to obtain a large permeability of surrounding rocks, but not to destroy the boreholes completely. The key parameters are summarized in Table 4, when the decoupling coefficient is 1.79, the boundary of the area where the permeability is 100 times of the initial value is about 64.29 cm from the center of borehole 5#, as marked in red in Figure 8.Figure 7
Crushing circle under the condition of decoupling coefficient: (a)K = 1.79; (b) K = 2.57; (c) K = 3.29. The crushing zone is defined as the area where the strain is greater than 1%.Table 4
Partial simulation results.
Model
Charge (cm × cm)
K
100 times region of influence (cm)
10 times region of influence (cm)
Radius of crushed zone (cm)
1
1.4 × 24
1.79
64.29
121.79
13.33
2
1.4 × 24
2.57
44.98
112.74
8.65
3
1.4 × 24
3.29
20.50
67.3
6.53Figure 8
Permeability variation: (a)K = 1.79; (b) K = 2.57; (c) K = 3.29; unit: μD. The initial permeability coefficient is 46.2 μD.Predicted damage of this model gives an affected area with the width 70 to 80 times of the charge radius at least, as shown in Figure8. The permeability changes dramatically along with the overall boundaries, where the boundary of the area with 100 times the initial permeability is about 64.29 cm from the borehole center, and the boundary of the area with 10 times the initial permeability is about 121.79 cm from the borehole center when the decoupling coefficient is 1.79. The value of the distance is 44.98 cm and 112.74 cm for 100 times and 10 times of the initial permeability, respectively, when the decoupling coefficient is 2.57.As an example, Figure9 presents the permeability coefficient with distance based on 4# hole, and Lines A, B, and C marked the boundaries of the permeability change of ten times for the three models, respectively, where A = 67.3 cm (K = 3.29), B = 112.74 cm (K = 2.57), and C = 121.79 cm (K = 1.79). The decoupling coefficient K obeys the arithmetic progression in the three models, as KA − KB = 0.72, KB − KC = 0.78. Comparing the region of influence and the size of the crushed zone, it is obvious that model 1 and model 2 have a similar range of permeability influence, which is larger than that of model 3, as has been discussed in the last paragraph. However, model 1 has a much larger crushed area than model 2 and model 3, which is not desirable for the blasting design. Therefore, the configuration in model 2 is better than that of model 1 and model 3.Figure 9
Permeability varies with distance.Figure10 compares the crush zone boundary and the 10 times initial permeability boundary of the three models. The data are normalized based on the result from model 3, and it is obvious that the slope of the crushing ratio curve of model 1 and model 2 is larger than the slope of the permeability change. The crushed zone area of model 1 is much larger than that of model 2, but the permeability influenced area is not obviously larger than model 2. In other words, model 2 (K = 2.57) achieves lower crush zone and higher performance-price ratio in increased permeability distribution range which is encouraged in engineering.Figure 10
Result normalized permeability boundaries and crushing range. The permeability coefficient is 462μD at all the permeability boundaries and initial permeability is 46.2 μD for all the three models.From the results in Table4, it can be noticed that the permeability coefficient of rock can be increased by two orders of magnitude in about 30–90 times charging radius and one order of magnitude in about 90–170 times charging radius (0.7 cm) based on the coupled damage-permeability dynamic constitutive model. And it is worth mentioning that crush boundary also needs to be considered to prevent the collapse of borehole in the pursuit of permeability coefficient in practice.


## 8. Conclusion

In this study, a rock dynamic constitutive model considering coupled damage-permeability effect is developed under explosive loading. The purpose is to provide a prediction for metallic ore mining using in situ leaching process by increasing the permeability in blasting. The constitutive model is controlled by the sum of all the tensile principal strains. The model was adopted to predict the crush zone and the permeability variation of a five-borehole blasting design, which is commonly used in mining industry. The comparison of the predicted results with the in situ experimental studies demonstrated the applicability and the effectiveness of the proposed dynamic constitutive model. The major contributions achieved in this study are as follows: First, the numerical prediction is proved to be valid since discrepancy of the permeability between the numerical prediction and the experimental test is slight and acceptable in practical applications. The main reason for the discrepancy is the heterogeneity of modeled material in the numerical prediction. Second, it is found that the permeability of rock can be increased two orders of magnitude in 80 times of explosive radius under the explosive load in proper decoupling coefficient. Third, the decoupling coefficient has a reasonable range to make optimized blasting effects—a smaller decoupling coefficient will lead to a larger crushed zone, while a larger decoupling coefficient affects the boundary of increasing permeability. Therefore, it is successful and useful to apply the proposed dynamic constitutive model for the prediction of the damage-permeability relationship, which is valuable for the optimization of blasting design.



---

*Source: 6816974-2018-07-22.xml*