Journal of Applied Mathematics

Volume 2013 (2013), Article ID 658160, 12 pages

http://dx.doi.org/10.1155/2013/658160

Research Article

Fracture Analysis of Brittle Materials Based on Nonlinear FEM and Application in Arch Dam with Fractures

Yuanwei Pan,¹ Yaoru Liu,¹ Zhixiong Cui,1,2 Xin Chen,³ and Qiang Yang¹

¹State Key Laboratory of Hydroscience and Hydraulic Engineering, Tsinghua University, Beijing 100084, China

²State Grid Xin Yuan Construction Co., Ltd., Beijing 100761, China

³China University of Mining & Technology, Beijing 100083, China

Received 7 June 2013; Revised 2 September 2013; Accepted 4 September 2013

Academic Editor: Pengcheng Fu

Copyright © 2013 Yuanwei Pan et al. 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.

Current fracture analysis models based on fracture mechanics or continuum damage mechanics are still limited in the application to three-dimensional structure. Based on deformation reinforcement theory coming from elastoperfect plastic theory, unbalanced force is proposed to predict initiation and propagation of cracks. Unbalanced force is the driving force of time-dependent deformation according to Perzyna’s viscoplasticity theory. It is also related to the damage driving force in viscoplastic damage model. The distribution of unbalanced force indicates cracks initiation area, while its direction predicts possible cracks propagation path. Uniaxial compression test of precrack specimen is performed as verification to this method. The trend and distribution of cracks are in good agreement with numerical results, proving that unbalanced force is feasible and effective for fracture analysis. The method is applied in fracture analysis of Xiaowan high arch dam, which is subjected to some cracks in dam due to the temperature control program. The results show that the deformation and stress of cracks and the stress characteristics of dam are insensitive to grouting of cracks. The existing cracks are stable and dam heel is still the most possible cracking position.

1. Introduction

Rock and concrete are heterogeneous anisotropic materials, containing numerous microcosmic voids and flaws. Fracture is a common and significant failure mode of geotechnical structure. Fracture evaluation for structure under certain loads, including crack initiation, propagation, and penetration, is still an unsolved problem. Thus fracture analysis method for rock and concrete structure is of significant importance in the sense of cracking prevention and global stability evaluation.

There is a key problem remaining in fracture analysis for brittle materials and structures, that is, a feasible fracture criterion. Common fracture criterions include stress criterions and energetic criterions [1–8]. The former state that failure occurs when the maximum principal stress in some local point exceeds the tensile strength. The latter are provided by linear elastic fracture mechanics, which covers some precise measurement, such as the stress intensity factor (SIF), energy release rate, and strain energy density. The stress criterions are acceptable for bodies without cracks, while the energetic criterions only work when a certain large crack already exists.

Current numerical methods on fracture analysis include fracture mechanics [9] and continuum damage mechanics [10]. Some researchers combine both methods for numerical failure analysis [11, 12]. Crack propagation simulation could be concluded as smeared fracture model [1] and discrete fracture model [13]. Besides, numerical tools have developed rapidly, including Finite Element Method (FEM) [14], Extended FEM [15], Element Free Method [16], Boundary Element Method [17], Discrete Element Method [18], Numerical Manifold Method [19], Discontinuous Deformation Analysis [20], and Fast Lagrangian Method [21].

The theories and methods mentioned above work well in planar analysis. There is still severe limitation on the applicability when extended to three-dimensional and complex structure. Besides, behavior of rock and concrete involves complex nonlinear overall deformation, which is beyond the capacity of common numerical methods.

This paper presents a new approach based on deformation reinforcement theory [22] to evaluate cracks growth in three-dimensional structure. Unbalanced force is proposed as a set of equivalent nodal forces of overstress beyond the yield surface. Unbalanced force drives time-dependent deformation, as well as damage evolution. The distribution of unbalanced force indicates crack initiation area, while its direction predicts potential crack propagation path. Uniaxial compression test of precrack specimen is performed, which is in a good agreement with the numerical results. The method is also applied in fracture analysis of Xiaowan high arch dam.

2. Fracture Analysis Method Based on Unbalanced Force

2.1. Deformation Reinforcement Theory

Most geotechnical structures are under complicated configurations and working conditions. Stability analysis of these engineering structure could be summarized as a boundary value problem with some basic equations, including kinematic admissibility, equilibrium condition, and constitutive equations. The classical elastoplastic theory aims at solving the displacement and stress fields that simultaneously satisfy all the aforementioned equations. However, the existence of such solution requires that the structure is stable. Structural instability occurs when action is greater than resistance. The difference between action and resistance is overstress, which is the key concept of the Deformation Reinforcement Theory (DRT) and defined as the unbalanced force.

For perfect elasto-plastic materials with associative flow rule, the constitutive equations can be stated as

where is defined as the plastic dissipation, . and are elastic and plastic strain rates, respectively. is the flexibility tensor. is the yield function and is an arbitrary stress state satisfying the yield criterion. Inequality (1) is known as principle of maximum plastic dissipation [23], which covers associative flow rule and the Kuhn-Tucker condition.

Assuming that is a pseudotime, the stress and plastic strain of material are and at time . With a strain increment in time interval , a new state and is achieved. Thus the linearized plastic strain rate is

where is termed the trial elastic stress while . is also an equilibrium stress field under given loads;

Applying (2) into Inequality (1), the process of solving real stress field is turned into the following minimization problem:

Equation (4) is known as the closest-point projection method (CPPM) [24], as shown in Figure 1. is an arbitrary stress field on the yielding surface, which represents the material resistance after previous minimization process. is an equilibrium stress field, which could be regarded as certain stress under external actions.

Figure 1: The closest-point projection method.

The minimization objective represents the difference between the plastic dissipations of the external action and the material resistance. It is defined as the volume density of the plastic complementary energy (PCE);

Thus, instability of a material point is equal to the statement that the external action is greater than the material resistance: , for all . As indicated by the minimization problem, if , the real stress state minimizes the resistance deficiency in the sense of PCE. In other words, material resistance capacity is fully developed.

Drucker-Prager yield criterion is adopted in the following nonlinear calculation:

where , , and and are material parameters. With the real stress field solved by (4), the following characteristic can be proved [22]:

Eventually the following analytical solution of can be derived:

where and are invariants of and and are the Lamé constants.

2.2. Demonstration in FEM

PCE and unbalanced force are implemented in FEM analysis. Since the process is restricted to displacement method, the kinematic admissibility is naturally satisfied.

In order to evaluate the global stability of structure under certain actions, two stress field sets, and , are defined for elastoperfect plastic problem. The two stress sets, respectively, satisfy the equilibrium condition and the yield criterion as

where is the displacement gradient matrix, is equivalent nodal force vector of external loads, and is structure volume. and , respectively, represent external actions and structural resistance. The real stress field must satisfy the yield criterion; . Equivalent nodal forces of the difference between the two stress fields, which is the driving force of structure deformation, can be defined as unbalanced force as

For a given load step, the unbalanced force is the driving force of deformation in the FE iterations. The principle of minimum PCE implies that elasto-plastic structures deform tending to the state where PCE is minimized under certain actions as

Unbalanced force is the measurement between external actions and structural resistance. The plastic complementary energy is a norm of unbalanced force. Once reaches its minimum value , the unbalanced force stays constant and the structure suffers steady plastic flow until failure occurs. Only if , the structure is stable.

In elastoplastic FEM iteration, PCE and unbalanced force can be calculated by the following steps. First, the nearest stable stress field can be achieved from an arbitrary equilibrium stress field . Then, the nearest equilibrium stress field is also obtained from . Finally, structure stress state tends to the two closest stress fields and , as shown in Figure 2. The iteration converges if the plastic complementary energy is reaching a steady value . For ideal model, is equal to its minimum value . If , the structure remains stable and the stress field is the real stress response. Otherwise, the structure is unstable and PCE is a magnitude estimation of its global instability. Unbalanced force reflects the structural failure behavior, including failure position and pattern.

Figure 2: Demonstration of the elastoplastic iteration in FEM analysis.

2.3. Explanation of Unbalanced Force in Viscoplastic Damage Model

Stress state beyond the yielding surface is unacceptable for perfect elastoplasticity. It exists only in the iteration process of FEM calculation. However, in viscoplasticity models, it is of explicit physical significance as the driving force of visco-plasticity deformation. The Perzyna associative visco-plasticity strain rate could be stated as [25, 26]

where is the visco-plasticity strain rate. is the viscosity parameter. is the yield function. is plastic potential function, and for associative flow rule, . is overstress function characterized by the yield function, commonly expressed as . is a reference constant of the same dimension with yield function .

According to Perzyna’s visco-plasticity theory, visco-plasticity strain rate is in proportion to overstress function. represents the direction of , while reflects the magnitude of . If the overstress function is restricted to yield function, plastic flow occurs only on the condition that the yield function value is positive. As the DRT states, unbalanced force is a set of the equivalent nodal forces of stress exceeding the yield function, which can be termed as the driving force of time-dependent deformation.

Unbalanced force is also involved in structural crack growth analysis based on damage mechanics methods. In the concept of continuum damage mechanics, the damage variable is usually related to equivalent plastic strain, that is, the accumulation of time-dependent deformation. A classical definition of damage variable is , where is total area and is intact area [27, 28]. Effective stress tensor in the undamaged configuration is determined by the damage variable and the nominal Cauchy stress tensor [29, 30] as

where is the effective stress tensor and is the nominal Cauchy stress tensor.

Darabi et al. [31] presented a thermoviscodamage model to describe the damage evolution law, which could be stated as

where is the effective total strain in the effective configuration. includes both viscoelastic and viscoplastic components. is the stress dependency parameter and is a material parameter. and are the reference damage viscosity parameter and the reference damage force, respectively. is a damage temperature function. is the damage driving force in the effective configuration, which can be assumed to have a modified Drucker-Prager-type form as

where and are the invariants of the effective stress tensor. gives the distinction of material behavior in compression and extension loading conditions, where implies that . Thus the damage driving force is expressed as

Damage evolution requires that ; that is, the equilibrium stress field exceeds the yield function. As stated above, the equivalent nodal force of stress exceeding the yield function is unbalanced force, which is the driving force of damage evolution. Unbalanced force is more explicit than equivalent plastic strain in the sense of physical significance.

2.4. Fracture Analysis Method

Unbalanced force is a set of equivalent nodal forces of plastic stress in the elements around a node, which reflects the difference between external actions and structural resistance. The process of FEM iteration is to find a set of additional forces to prevent failure from occurring, which is minimized in the sense of PCE. If the additional force, that is, unbalanced force, tends to be zero in the iteration, the cracks will not initiate or stay in a limit state. Otherwise, cracks will initiate and propagate, and the direction and distribution of unbalanced force indicate the potential failure area.

Thus, occurrence of unbalanced force can be the identification of local cracks initiation. Meanwhile, as unbalanced force is related to the damage driving force in the thermo-viscodamage model, the distribution of unbalanced force will expand gradually as material damage accumulates. This process corresponds to the cracks growth and propagation. In general, its direction predicts the possible path of crack’s propagation.

Current fracture analysis methods are based on planar problems, which provide good results dealing with a single crack or two cracks. When extended to 3D structures, the fracture criterion and nonlinear calculation efficiency remain to be solved. Unbalanced force is presented in this paper, which is of clear physical concept and could be effectively applied to elato-plastic FEM calculation.

3. Fracture Test and Numerical Analysis of Precrack Specimen

The strong correlation is found between the distribution of unbalanced force and the initial point of structural cracking, which is verified by the following model test and numerical results.

3.1. Specimen and Numerical Model

The precrack cubic specimen is made with gypsum, which is 15 cm wide, 15 cm high, and 5 cm thick, as shown in Figure 3. Steel slices are inserted during the specimen modeling process and pulled out once the initial set begins. By this means the precracks are manufactured. The angle between precrack and horizontal direction is 30°, and precrack length is 12 mm.

Figure 3: Sketch map of the precrack specimen.

Numerical model is exactly established according to the physical test. Precracks are simulated where the corresponding elements are set null. The precrack specimen and numerical model are illustrated in Figure 4.

(a)

(b)

Figure 4: The precrack specimen and numerical model.

3.2. Test and Calculation Results

Uniaxial compression test of the precrack specimen is performed with its bottom restrained. The pressure on the top is gradually increased from 0 MPa to 2.0 MPa by the increment of 0.1 MPa.

The precracks in the left part of specimen initiate first during the uniaxial compression test, when the pressure achieves 1.6 MPa. The cracks initiation direction is almost perpendicular to the precracks trend. With the pressure being gradually increased, cracks propagation path tends to be vertical, and precracks at the top right corner begin to initiate as well. Cracks propagate and eventually penetrate both sides of the specimen, where failure occurs with a compressive strength 2.0 MPa. The cracks growth process and final failure mode are shown in Figure 5.

(a)

(b)

(c)

(d)

Figure 5: Cracks growth and failure mode of specimen test.

Nonlinear numerical calculation is applied with the same constraint conditions and loading procedure. As clearly shown in (14), the damage evolution law is influenced by various factors, including stress, strain, strain rate, temperature, and damage history. The loading process duration of this test is relatively short, so the temperature effect is ignored. Since gypsum is quasi-brittle material, the time-dependent deformation of damage evolution process is simplified in the numerical calculation. For the region where the damage driving force is greater than the reference damage force , gypsum is regarded as damaged and elastic modulus , and shear strengths and are reduced by 50%. In this case, MPa. Table 1 shows the materials parameters.

Table 1: Materials parameters.

E/GPa
μ
/ kN·m⁻³
f
c/MPa

Gypsum
5
0.35
2
0.1
0.2

Damaged gypsum
2.5
0.35
2
0.05
0.1

Unbalanced force distribution and evolution process at crack tips are shown in Figure 6. When the specimen stays in a stress state of elasticity, there is no unbalanced force. With the pressure being increased, unbalanced force occurs at the precrack tips around the middle left part of the specimen. As the distribution area expands, the direction of unbalanced force growth tends to be perpendicular to the precracks. The unbalanced force eventually goes through the nearest crack tips, which indicates the most possible propagation path of specimen fracture.

(a)

(b)

(c)

(d)

Figure 6: Unbalanced force evolution process at crack tips.

Damaged area expanding process is shown in Figure 7. Damaged area appears around the tip of crack and gradually extends to neighbor elements in the perpendicular direction to the precracks. At last, the damaged area unites between the tips of the two parallel precracks.

(a)

(b)

(c)

(d)

Figure 7: Damaged area expansion process at crack tips.

The existence of unbalanced force indicates that the model is unable to balance the loads and fracture occurs, which agrees with the test results. As cracks propagate, a new structure is achieved and the stress field is redistributed. Since the crack propagation is a local and quasi-static failure process, the new structure retains its bearing capacity. Thus the model is able to sustain the loads and reaches a new equilibrium state. This process continues until the load reaches the compressive strength and structural failure occurs. The distribution of unbalanced force indicates cracks initiation area, while its direction predicts the possible cracks propagation path.

4. Fracture Analysis of Arch Dam with Cracks

4.1. Numerical Model

Xiaowan arch dam, with a height of 292 m, is subjected to some cracks in dam due to the temperature control program. Most of the cracks distribute on the middle and bottom elevations of the dam and tend tangentially to the arch. The average width of these cracks is 1 mm. Distribution of cracks on typical elevation and on typical dam section is shown in Figures 8 and 9, respectively.

Figure 8: Distribution of cracks on typical elevation.

Figure 9: Distribution of cracks in typical dam section.

The size of numerical model is as follows: upstream 300 m, downstream 900 m, and left and right banks 700 m each. Elevations of dam crest and model bottom are 1245 m and 595 m, and thus the model height is 650 m. The model includes 58989 nodes and 53850 elements, as shown in Figure 10. There are 11 major cracks that are numbered and simulated, as shown in Figure 11 and Table 2.

Table 2: Numbers of cracks in the dam.

Number
Crack

1
13-1

2
13-2

3
20-1

4
20-2

5
22-1

6
22-3

7
22-4

8
25-1

9
28-1

10
28-2

11
30-1

Crack 13-1 in the table means the 1st crack in dam section number 13, and so on.

Figure 10: Numerical model of Xiaowan arch dam.

Figure 11: Distribution of cracks in the dam.

The dam cracks are induced by the improper temperature control program during construction. Since concrete grouting in the cracks is operated, the dam cracks are simulated with thin layer elements. The material parameters are shown in Table 3.

Table 3: Material parameters.

E/GPa
μ
/kN·m⁻³
f
c/MPa

Foundation
18
0.28
28
1.52
1.33

Dam
23.1
0.2
24
1.4
1.6

Crack
5
0.3
27.5
1.12
0.9

4.2. Fracture Analysis

Structural fracture and parameters sensitivity analysis are presented. Three schemes are performed, including normal parameters (scheme 1), crack’s shear strengths and reduced by 50% (scheme 2), and crack’s shear strengths and reduced by 75% (scheme 3). Both the hydrostatic pressure and overload conditions are included in each scheme.

4.2.1. Dam Heel Cracking Analysis

The distribution of unbalanced force at dam heel under different loads in scheme 1 is shown in Figure 12. The other two schemes are of similar distribution. Under the hydrostatic pressure , unbalanced force rarely occurs and dam heel stays as a global steady state. As the hydrostatic pressure is overloaded, unbalanced force area expands and its magnitude increases rapidly, which is concentrated in the dam heel.

(a)

(b)

(c)

Figure 12: Distribution of dam heel unbalanced force of scheme 1.

Unbalanced forces of dam heel in three schemes are summarized in Table 4. In the overload conditions, unbalanced force increases rapidly in dam heel, which is the most possible location of cracks initiation. As the variation of unbalanced forces between the three schemes is relatively small, the dam heel cracking behavior is insensitive to the change of crack parameters.

Table 4: Unbalanced forces of dam heel (10⁴ N).

Scheme  

1.5
2.0

1
40249.56
257289.9
415601.5

2
40216.39
257269.2
415560.7

3
40126.59
256869.2
414973.0

4.2.2. Dam Body Cracking Analysis

The distribution of unbalanced force in dam under different loads in scheme 1 is shown in Figure 13. Under the hydrostatic pressure , there is a little unbalanced force existing in the right abutment of the arch dam. As the pressure is overloaded, unbalanced force occurs in both the dam abutments. Distribution of unbalanced force in dam is close to the foundation, where cracks initiate and failure may occur.

(a)

(b)

(c)

Figure 13: Distribution of unbalanced force in dam of scheme 1.

Distribution of dam body unbalanced force under overload in schemes 2 and 3 is shown in Figures 14 and 15, respectively. Dam unbalanced force is insensitive to the change of crack parameters. Crack number 10, that is, the second crack in dam section number 28, suffers from some unbalanced force in scheme 3, which indicates the most sensitive and dangerous crack in the dam. Meanwhile, the unbalanced force of cracks is much less than that of dam body, and the latter is still the emphasis of dam cracking prevention.

Figure 14: Distribution of dam body unbalanced force under overload in scheme 2.

Figure 15: Distribution of dam body unbalanced force under overload in scheme 3.

4.2.3. Analysis for Existing Cracks

The length of unbalanced force vectors in the cracks is enlarged 20 times since it is relatively less than that of dam body. Distribution of cracks unbalanced force under overload in the three schemes is shown in Figures 16, 17, and 18, respectively. Cracks initiation and propagation are very sensitive to the crack parameters. In other words, the stability of the existing cracks relies on the quality of concrete grouting in the cracks.

Figure 16: Distribution of cracks unbalanced force under overload in scheme 1.

Figure 17: Distribution of cracks unbalanced force under overload in scheme 2.

Figure 18: Distribution of cracks unbalanced force under overload in scheme 3.

Unbalanced force of cracks under overload in all three schemes is summarized in Table 5. Unbalanced force in dam heel is also listed in Table 5 as comparison. Unbalanced forces in dam heel increase earlier than the cracks in the dam. Dam heel contributes most of the unbalanced forces in all conditions. Namely, dam heel cracking occurs before any existing crack propagates. Among all existing cracks, 20-2 and 28-2 are the dominating cracks in the process of fracture propagation, while 30-1 and 13-2 are also possible initiation area, as shown in Figure 19. The cracks near the crown cantilever are stable.

Table 5: Unbalanced force of cracks under overload (10⁴ N).

Number
Crack
Scheme   1
Scheme   2
Scheme   3

1
13-1
0.00
0.00
12.54

2
13-2
0.02
0.87
221.69

3
20-1
0.00
0.00
2.83

4
20-2
0.03
0.70
7614.90 *

5
22-1
0.00
0.00
0.00

6
22-3
0.00
0.00
0.00

7
22-4
0.00
0.00
0.00

8
25-1
0.00
0.00
8.00

9
28-1
0.00
0.94
65.95

10
28-2
1.77
345.54
15810.91 *

11
30-1
0.02
6.04
687.85

Dam heel
415601.5
415560.7
414973

*Dominating cracks in the propagation process.

Figure 19: Local distribution of cracks unbalanced force (scheme 1, under overload ).

4.3. Comparison with Geomechanical Model Test

Results and conclusions of geomechanical model test are presented as comparison with numerical method. The final distribution of upstream cracks is shown in Figure 20. Dam heel cracking occurs as the overload increases to . There is no sign of existing cracks growth in the dam during the test. Instead, cracks that occur on the dam surface begin to extend after the overload reaches . Experimental results indicate that dam heel cracking, compared with existing cracks in dam, is the dominating problem of Xiaowan arch dam, which is corresponding to numerical results.

Figure 20: The final distribution of upstream cracks.

5. Conclusion

Unbalanced force is proposed based on deformation reinforcement theory to analyze fracture behavior, including initiation and propagation of cracks in 3D structures. Unbalanced force is a set of the equivalent nodal forces of stress exceeding the yield function, which can be termed the driving force of time-dependent deformation, as well as damage evolution.

The unbalanced force and damaged area are in good agreement with precrack specimen test results. The distribution of unbalanced force indicates cracks initiation area, while its direction predicts the possible cracks propagation path.

The method is applied in fracture analysis of Xiaowan high arch dam. Dam heel cracking occurs before any existing crack propagates, which is the most possible failure mode. Among all existing cracks, 20-2 and 28-2 are the dominating cracks in the process of fracture propagation.

Acknowledgments

The work reported here was supported by the State Key Laboratory of Hydroscience and Engineering of Hydroscience with Grant no. 2013-KY-2 and China National Funds for Distinguished Young Scientists with Grant no. 50925931.

References

A. Hillerborg, M. Modéer, and P. E. Petersson, “Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements,” Cement and Concrete Research, vol. 6, no. 6, pp. 773–781, 1976.

Z. P. Bazant and B. H. Oh, “Crack band theory for fracture of concrete,” Meterials and Structures, vol. 16, no. 3, pp. 155–177, 1983.

K. J. Willam and E. P. Warnke, “Constitutive models for the triaxial behavior of concrete,” in Proceedings of the International Association for Bridge and Structural Engineering, vol. 19, p. 174, ISMES, Bergamo, Italy, 1975.

J. Kemeny and N. G. W. Cook, “Effective moduli, non-linear deformation and strength of a cracked elastic solid,” International Journal of Rock Mechanics and Mining Sciences and, vol. 23, no. 2, pp. 107–118, 1986.

G. R. Irwin, “Fracture dynamics,” in Fracturing of Metals, pp. 147–166, 1948.

E. Orowan, “Fracture and strength of solids,” Reports on Progress in Physics, vol. 12, no. 1, pp. 185–232, 1948.

J. R. Rice, “A path independent integral and the approximate analysis of strain concentration by notches and cracks,” Journal of Applied Mechanics, vol. 35, pp. 379–386, 1968.

G. C. Sih and B. Macdonald, “Fracture mechanics applied to engineering problems-strain energy density fracture criterion,” Engineering Fracture Mechanics, vol. 6, no. 2, pp. 361–386, 1974.

W. Yang, Macro-Micro Fracture Mechanics, National Industry Press, Beijing, China, 1995 (Chinese).

M. G. D. Geers, R. De Borst, and R. H. J. Peerlings, “Damage and crack modeling in single-edge and double-edge notched concrete beams,” Engineering Fracture Mechanics, vol. 65, no. 2-3, pp. 247–261, 2000.

W. S. Blackburn, “Remeshing algorithm for three-dimensional crack growth and intersection with surfaces or cracks in non-coplanar planes,” Engineering Analysis with Boundary Elements, vol. 24, no. 4, pp. 343–350, 2000.

J. C. W. van Vroonhoven and R. de Borst, “Combination of fracture and damage mechanics for numerical failure analysis,” International Journal of Solids and Structures, vol. 36, no. 8, pp. 1169–1191, 1998.

Z. P. Bazant, “Crack band model for fracture of geomaterials,” in Proceedings of the 4th International Conference on Numerical Methods in Geomechanics, vol. 3, pp. 1137–1152, Alberta, Canada, 1982.

P. O. Bouchard, F. Bay, Y. Chastel, and I. Tovena, “Crack propagation modelling using an advanced remeshing technique,” Computer Methods in Applied Mechanics and Engineering, vol. 189, no. 3, pp. 723–742, 2000.

N. Moës and E. Béchet, “Modeling stationary and evolving discontinuities with finite elements,” in Proceedings of the 7th International Conference on Computational Plasticity (COMPLAS '03), CIMNE, Barcelona, Spain, 2003.

T. Belytschko, Y. Y. Lu, and L. Gu, “Element-free Galerkin methods,” International Journal for Numerical Methods in Engineering, vol. 37, no. 2, pp. 229–256, 1994.

M. F. Marji, H. Hosseini-Nasab, and A. H. Kohsary, “On the uses of special crack tip elements in numerical rock fracture mechanics,” International Journal of Solids and Structures, vol. 43, no. 6, pp. 1669–1692, 2006.

A. Baghbanan and L. Jing, “Stress effects on permeability in a fractured rock mass with correlated fracture length and aperture,” International Journal of Rock Mechanics and Mining Sciences, vol. 45, no. 8, pp. 1320–1334, 2008.

G. H. Shi, “Manifold method,” in Proceedings of the 1st International Forum on DDA and Simulations of Discontinuous Media, Bekerley, Calif, USA, 1996.

Y. H. Hatzor, A. A. Arzi, Y. Zaslavsky, and A. Shapira, “Dynamic stability analysis of jointed rock slopes using the DDA method: king Herod's Palace, Masada, Israel,” International Journal of Rock Mechanics and Mining Sciences, vol. 41, no. 5, pp. 813–832, 2004.

H. Hakami, “Rock characterisation facility (RCF) shaft sinking—numerical computations using FLAC,” International Journal of Rock Mechanics and Mining Sciences, vol. 38, no. 1, pp. 59–65, 2001.

Q. Yang, Y. Liu, Y. Chen, and W. Zhou, “Deformation reinforcement theory and its application to high arch dams,” Science in China E, vol. 51, no. 2, pp. 32–47, 2008.

R. Hill, The Mathematical Theory of Plasticity, Clarendon Press, Oxford, UK, 1950.

J. C. Simo, J. G. Kennedy, and S. Govindjee, “Nonsmooth multisurface plasticity and viscoplasticity. Loading/unloading conditions and numerical algorithms,” International Journal for Numerical Methods in Engineering, vol. 26, no. 10, pp. 2161–2185, 1988.

O. C. Zienkiewicz and I. C. Cormeau, “Viscoplasticity-plasticity and creep in elastic solids-a unified numerical solution approach,” International Journal for Numerical Methods in Engineering, vol. 8, no. 4, pp. 821–845, 1974.

P. Perzyna, “Fundamental problems in viscoplasticity,” Advances in Applied Mechanics, vol. 9, no. 2, pp. 243–377, 1966.

L. M. Kachanov, “Time of the rupture process under creep conditions,” Izvestiya Akademii Nauk SSR Otdelenie Tekhnicheskikh Nauk, vol. 8, pp. 26–31, 1958.

I. N. Rabotnov, Creep Problems in Structural Members, vol. 7, North-Holland, Amsterdam, the Netherlands, 1969.

U. Cicekli, G. Z. Voyiadjis, and R. K. Abu Al-Rub, “A plasticity and anisotropic damage model for plain concrete,” International Journal of Plasticity, vol. 23, no. 10-11, pp. 1874–1900, 2007.

R. K. Abu Al-Rub and G. Z. Voyiadjis, “Gradient-enhanced coupled plasticity-anisotropic damage model for concrete fracture: computational aspects and applications,” International Journal of Damage Mechanics, vol. 18, no. 2, pp. 115–154, 2009.

M. K. Darabi, R. K. Abu Al-Rub, E. A. Masad, C. W. Huang, and D. N. Little, “A thermo-viscoelastic-viscoplastic-viscodamage constitutive model for asphaltic materials,” International Journal of Solids and Structures, vol. 48, no. 1, pp. 191–207, 2011.