Shot Peening Processes to obtain Nanocrystalline surfaces in metal alloys: S.C. Li et alii, Frattura ed Integrità Strutturale, 47 (2019) 1-16; DOI: 10.3221/IGF-ESIS.47.01 1 Numerical simulation method study of rock fracture based on strain energy density theory Shuchen Li, Tengfei Ma, Luchen Zhang, Qian Sun Geotechnical and Structural Engineering Research Center, Shandong University, Jinan, Shandong 250061, China 316159025@qq.com ABSTRACT. Many numerical methods are carried out to study the nonlinear failure behaviors of the rock; however, the numerical simulation methods for the failed rock are still in the research stage. This paper establishes the damage constitutive equation by combining the bilinear strain softening constitutive model with energy dissipation principles, as well as the energy failure criterion of mesoscopic elements based on the strain energy density theory. When the strain energy stored by an element exceeds a fixed value, the element enters the damage state and the damage degree increases with increasing energy dissipation. Simultaneously, the material properties of the damaged element change until it becomes an element with certain residual strength. As the load increases, the damage degree of an element increases. When the strain energy stored by an element exceeds the established value of the energy criterion, the element is defined to be failed. As the number of failed elements constantly increases, failed elements interconnect and form macrocracks. The rock fracture calculation program on the basis of the preceding algorithm is successfully applied to the fracture simulation process in Brazilian splitting, tensile tests with build-in crack and tunnel excavation. KEYWORDS. Strain energy density; Energy dissipation; Rock fracture; Flac; Numerical simulation; Bilinear softening Citation: Li S.C., Ma, T.F., Zhang, L.C., Sun, Q., Study on numerical simulation method based on strain energy density theory of rock fracture, Frattura ed Integrità Strutturale, 47 (2019) 1-16. Received: 22.07.2018 Accepted: 24.11.2018 Published: 01.01.2019 Copyright: © 2019 This is an open access article under the terms of the CC-BY 4.0, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. INTRODUCTION s a product of geological movements, the rock is a heterogeneous material with complex mechanical properties. The rock medium usually shows strong nonlinear characteristics [1] in the deformation process. Physical and mechanical properties will be irreversible during the failure process, and this irreversible process will cause various forms of energy dissipation. According to the law of thermodynamics, energy conversion is the essential feature during the material physical process, and material failure is a kind of instability driven by energy. Therefore, studying the energy change rules and establishing the relation between energy changes and strength as well as structural failure during the rock A http://www.gruppofrattura.it/VA/47/2136.mp4 S.C. Li et alii, Frattura ed Integrità Strutturale, 47 (2019) 1-16; DOI: 10.3221/IGF-ESIS.47.01 2 failure process will be more conducive to reflecting strength changes and essential characteristics [2] of structural failure of the rock under the action of external load. Since Huber firstly introduced the potential energy concept to define material damage, many scholars at home and abroad have described rock deformation behaviors through energy analysis and achieved tremendous progress [3-8]. Most work focused on tests and theoretical research and macroscopic failure behaviors of the rock were obtained through energy analysis. According to their studies, material failures are mainly caused by irreversible internal energy dissipation; meanwhile, the energy criterion is generally significant for determining the rock failure. The rock fracture is an entire process from damage, material progressive degradation, microcrack generation, expansion, and till penetration. Therefore, studying rock fracture from the micromechanics perspective, analyzing the damage rules of tiny rock elements, and studying element failure through element energy dissipation and energy criterion can systematically show the entire rock fracture process. Currently, the most representative strain energy failure criterion is distortional strain energy density theory (fourth strength theory) [9]. This is a better strength theory for plastic materials; however, it is applicable only to plastic materials with the same tension and compression properties rather than triaxial equivalent tension. Therefore, besides the distortional strain energy, the volume deformation energy also needs to be considered [10]. The application of strain energy density theory [11] can comprehensively consider the preceding problems and take the strain energy density that is the sum of the volume deformation energy density and shape change energy density as the criterion of material failure. The advantage of this theory is that it can be well applied to complex geometry, loading conditions, and development situations of mixed cracks [12]. According to the preceding analysis, nonlinear failure behaviors of the rock are considered to simulate the rock partial failure and the entire process from microcrack generation, expansion, and complete fracture. In addition, the bilinear strain softening constitutive model is adopted by referring to document [13], and the energy criteria for the damage and failure of mesoscopic rock elements according to the strain energy density theory and energy dissipation principle. When the strain energy stored by an element exceeds a fixed value, the element enters the damage status. The damage degree of the element is determined based on the strain energy density and the material properties of the damaged element change until it becomes an element with certain residual strength. For rock elements entering the damage state, the energy failure criterion of strain energy density is used to determine whether the element is damaged. As the load increases, the number of failed elements gradually grows, and the failed elements interconnect and form macrocracks, causing the structural failure of the rock specimen. During the numerical simulation process, the elastic modulus reduction of damaged elements after reaching the stress extreme value is discretized. The preceding method completes the nonlinear calculation process with linear calculation, avoids singularity of numerical calculation in element fracture, and simulates the rock post-peak fracture behaviors. In this calculation method, the rock fracture calculation program is developed with Fish language in the Flac. It is successfully applied to the fracture simulation process in Brazilian splitting, tensile tests with build-in crack and tunnel excavation, indicating the accuracy and feasibility of this method for simulating the rock fracture process. STRAIN ENERGY DENSITY THEORY he rock fracture mode is affected by many factors such as the loading type, geometry, and material properties. However, the strain energy density theory can comprehensively consider these factors. The rock is assumed to be a continuum composed of many tiny structural elements and each element contains per unit volume of the material. If the element deforms under the action of external force, strain energy will be stored inside the element. In this way, the energy stored by each element is called strain energy density ( / )dW dV . The strain energy density equation can be expressed as follows: 0 ( , ) ij ij ij dW d f T C dV   = +   (1) ij and ij are stress and strain components. T and C are the temperature and humidity variations, in general, when the temperature and humidity are basically constant, this part is negligible. So the strain energy density equation can be expressed as 0 / ij ij ijdW dV d   =  simply. For linear elastic loaded objects, the elastic strain energy is equal to the work done by the external force, and the strain energy stored in the material element depends only on the final value of the external force and the deformation, but has T S.C. Li et alii, Frattura ed Integrità Strutturale, 47 (2019) 1-16; DOI: 10.3221/IGF-ESIS.47.01 3 nothing to do with the loading sequence. Therefore, in the case of linear elasticity, the relationship between each principal stress and the corresponding principal strain remains linear, so the strain energy density of the triaxial stress-state is as follows: 1 1 2 2 3 3 1 ( ) 2 dW dV      = + + (2) According to the general Hooke's law: 1 1 2 3 2 2 3 1 3 3 1 2 1 [ ( )] 1 [ ( )] 1 [ ( )] v E v E v E              = − +    = − +   = − +  (3) Substituting formula (3) into formula (2), we can get it: 2 2 2 1 2 3 1 2 2 3 3 1 1 [ 2 ( )] 2 dW dV E          = + + − + + (4) According to Mohr stress circle theory of space state, the stress components on each surface are substituted into formula (4), in the case of linear elasticity, the general expression of strain energy density equation is as follows: 2 2 2 2 2 21 1 ( ) ( ) ( ) 2 x y z x y y z z x xy yz zx dW dV E E E               + = + + − + + + + + (5) Therefore, in the case of plane stress state, the elastic strain energy density equation of the element body is expressed as follows: 2 2 21 1 ( ) 2 x y x y xy dW dV E E E        + = + − + (6) Considering the stress-strain curves of materials under tensile conditions, as shown in the Fig. 1, Assuming that the stress continues to increase after reaching the yield stress Y , plastic deformation occurs. If unloading at point P, the unloading path will be along line PM, and the new loading path will be along line MPF. In the process of unloading and reloading, energy represented by area OAPM= ( / ) pdW dV is dissipated. Therefore, the effective energy for crack propagation * ( / ) c dW dV can be expressed in OAPM, or as the formula shows: * c c p dW dW dW dV dV dV       = −            (7) This formula represents the total energy required for unit volume when the material unit element fails. Strain energy density theory is a good failure criterion for predicting nonlinear damage phenomena. The failure of materials is generally the process of stable crack development until the global instability of the structure. Stable crack growth is a local or microstructural instability that can be predicted by the critical strain energy density ( / )cdW dV , the value can be obtained from the area under the complete stress-strain curve. No matter what stress the element bears, such as tension, compression, or shear stress, the strain energy density can comprehensively reflect the action of each stress component on the element. Each element can store limited strain energy S.C. Li et alii, Frattura ed Integrità Strutturale, 47 (2019) 1-16; DOI: 10.3221/IGF-ESIS.47.01 4 at a specified time and the strain energy stored by the element varies with different parts of the material. Therefore, the damage mode of the material can be evaluated based on the energy change process of the material from one element to another. Figure 1: Dissipative strain energy and recoverable strain energy under tensile condition Determination of the strain energy density function considering strain softening According to the strain energy density theory, the strain energy density function of each element in the constant humidity and temperature condition can be expressed as: 0 ij ij ij dW d dV   =  (8) According to the preceding formula, the density of strain energy stored in an element is determined by its stress ij and the deformation history of strain increment ijd . This theory determines the yield failure of the material element based on strain energy density ( / )dW dV . Limit value of strain energy density for failure of unit element, is determined by yielding test in uniaxial tension, and stipulate that, the total strain energy absorbed by the unit is equal to the energy absorbed ( / )cdW dV at fracture under uniaxial tension, the unit will yield failure. According to this theory, before the rock texture shows global instability under the action of load, partial failure and crack propagation already occur and seriously affect the macroscopic failure behaviors of the rock texture. The steady propagation of microcracks inside the rock will finally cause rock macroscopic fracture. The deformation process of each rock element is accompanied with energy dissipation that will cause material progressive damage, property deterioration, and strength loss [2]. The material mechanical damage can be described by strain softening. The bilinear strain softening constitutive model of the rock under uniaxial tension in document [13] is considered, that is, the uniaxial tension stress-strain relation of the rock is simplified, as shown in Fig. 2. Under the action of external force, the rock first shows elastic deformation and does not fail immediately after reaching the stress limit point U. Instead, it enters the strain softening stage and the material starts to damage. As shown in the Fig. 2, the density ( / )dW dV of strain energy absorbed by the material element at point A is area OUAC surrounded by the stress-strain curve. If the element is unloaded at point A, the unloading path will be along line AB and the new loading path will be along line BAF. In the process of unloading and reloading, energy represented by area OUAB is dissipated. Therefore, the density of strain energy absorbed by the material element is composed of the following two parts: ( ) ( )d e dW dW dW dV dV dV = + (9) In the preceding formula, ( / )ddW dV is the density of strain energy dissipated OUAB while ( / )edW dV is the density of strain energy recoverable BAC. Strain A P F O M E F’ p dV dW       * c dV dW       Y  Stress S.C. Li et alii, Frattura ed Integrità Strutturale, 47 (2019) 1-16; DOI: 10.3221/IGF-ESIS.47.01 5 Figure 2: Bilinear strain softening stress-strain curve If point A is on line OU in the elastic stage, the unloading path will be along the loading path and the material will not damage. For the material element without damage, its critical strain energy density ( / )cdW dV is equal to area OUF. However, energy dissipation occurs after the material element damages. The residual critical strain energy density after damage *( / )cdW dV is the density of strain energy regained after the element is unloaded (that is area BAF) and can be expressed as the following: * ( ) ( ) ( )c c d dW dW dW dV dV dV = − (10) The preceding formula indicates that a higher density of strain energy dissipated by the material element means more serious element damage and lower density of critical strain energy that can be borne. According to the preceding analysis, larger deformation of the material element under the action of external force means a higher density of strain energy absorbed ( / )dW dV and lower density of critical strain energy *( / )cdW dV . When * ( / ) ( / )cdW dV dW dV , cracks start to be generated. When ( / )dW dV is equal to the initial critical strain energy density ( / )cdW dV (that is, area OUF) of the material element, the element is completely fractured and cannot bear any load. The bilinear constitutive relation of the preceding rock element can be simply obtained through uniaxial tensile tests. The rock texture is usually under the action of complex external force and the internal rock elements are under the combined action of tension, compression, and shear stress. The strain energy density can be obtained based on the loading history of the rock element and it can comprehensively reflect the loading status of the element. The damage status of the rock element under complex stress conditions can be determined by comparing it with the strain energy density in different stages under uniaxial tension. Energy dissipation and damage constitutive model According to energy opinions, after the rock shows inelastic deformation, the inelastic deformation energy the rock can bear has been significantly reduced, that is, the rock constitutive energy has decreased. This is also an expression of rock performance deterioration caused by the changes of rock microstructure [15]. After the rock element reaches the peak point of the tensile stress, it enters the strain softening stage and shows inelastic deformation and decreased material strength. The strength decrease of the material element is defined as the elastic modulus reduction and is expressed by equivalent elastic modulus *E . As shown in Fig. 3, as the energy loses, the unloading peak strength decreases from point U to points G, H, I ... and the equivalent elastic modulus is * 1E , * 2E , * 3E , …and * nE respectively. For the sake of calculation, the equivalent elastic modulus is discretized into 20 different values: * ( 21 ) ( ) 20 n E n E − = (11) n=1, 2, ..., 20 S.C. Li et alii, Frattura ed Integrità Strutturale, 47 (2019) 1-16; DOI: 10.3221/IGF-ESIS.47.01 6 The following can be obtained according to the constitutive relation of continuous medium damage mechanics theory [16]: (1 )E D = − (12) E is the elastic modulus of the rock without damage. The damage parameter D of the rock element in this model is: * ( )E E n D E − = (13) Figure 3: Reduction of elastic modulus after energy dissipation n refers to the damage degree of the material element. For the discretization of the equivalent elastic modulus, a larger value of n means more accurate simulation results but remarkably lower calculation efficiency. n is set to 20 by taking both factors into consideration, which meets related requirements. According to the energy dissipation principle, energy dissipation is directly related to the damage and strength loss. The dissipation quantity reflects the reduction of the initial strength [2]. According to Fig. 2 and formula (9), both the strain energy density and the energy dissipation increase after the rock element reaches the stress limit point. Therefore, the strain energy density is used as the criterion to determine the rock element damage and failure. (1) When ( / ) ( / ) 1/ 2( )u u udW dV dW dV   = , the material element is in the elastic stage; no damage occurs; both the equivalent elastic modulus and critical strain energy density are the initial values of the element, that is, *E E= , * ( / ) ( / ) 1/ 2( )c c u fdW dV dW dV  = = . (2) When ( / ) ( / )udW dV dW dV= , the material element enters the damage stage. The discretized equivalent elastic modulus * ( )E n is regarded as the elastic modulus of the material element after damage. n refers to the damage degree of the material element ( 0 20)n  ) and its value is determined by the strain energy density. It also can be seen from Fig. 2 that the critical strain energy density *( / )cdW dV decreases with increasing ( / )dW dV after the material element enters the damage stage. (3) When *( / ) ( / )cdW dV dW dV= , the element fails, marking the beginning of crack generation inside the material. (4) When ( / ) ( / ) 1/ 2( )c u fdW dV dW dV  = = , according to formula (13), the damage degree of the material element reaches the largest (n = 20); both the equivalent elastic modulus and critical strain energy density of the element change to zero; the element is completely fractured and loses the bearing capacity. To maintain the integrity of the entire structure calculation model and element continuity, the element completely fractured will be given a very small residual modulus * 0.05cE E= instead of removing the element. Strain 1 * E 2 * E 3 * E U G H I nE * F Stress S.C. Li et alii, Frattura ed Integrità Strutturale, 47 (2019) 1-16; DOI: 10.3221/IGF-ESIS.47.01 7 Numerical implementation The FLAC2D finite difference numerical software is adopted to establish the rock element damage equation based on the preceding theory and develop a rock fracture calculation program with Fish language. Model loading adopts the load control stepwise loading method. In particular, equivalent tiny loads are accumulatively added to the model in turn. After calculation is balanced, the strain energy density of each element is calculated. Then, the preceding strain energy density criterion is used to determine the element damage degree and failure. After the first loading step is performed, all elements are in the elastic stage due to the small load and the strain energy density of each element is: 2 2 2 21 [ 2 2 ( )] 2 x y xy x y xy dW v dV E      = + + − − (14) After loading step i is performed, the strain energy density of each element is: 1 1 ( 1) ( 1) 1 1 1 1 1 ( ) ( ) ( ) ( ) ( )( ) 2 1 1 ( )( ) ( )( ) 2 2 i i i i i i i x x x x i i i i i i i i y y y y xy xy xy xy dW dW dW dW dV dV dV dV             − − − − − − − − = +  = + + − + + − + + − (15) In the formula, i≥2 i x , i y , and i xy are the element stress in loading step i while 1i x − , 1i y − , and 1i xy − are the element stress in the last loading step. Similarly, the strain has the same expression way. When the strain energy density of an element ( / ) ( / ) 1/ 2( )u u udW dV dW dV   = , the damage degree n is: ( / ) ( / ) 2( / ) 20 1 [ ( / ) ( / ) ] 20 u u u u f u u c u dW dV dW dV dW dV n dW dV dW dV       − − = = − − (16) n is set to the integer part of the calculation results on the right of the equation and ranges from 0 to 20. When 0n = , the element is still in the linear elastic stage and no damage occurs; when 20n = , the damage value is the maximum, indicating that the element is completely fractured and loses the bearing capacity. Meanwhile, the element elastic modulus * E is obtained by using formula (11). The density of strain energy dissipated by the element is: 2 2 2 2 * 1 ( / ) ( / ) [ 2 2 ( )] 2 d x y xy x y xydW dV dW dV v E      = − + + − − (17) Therefore, the critical strain energy density of the element is: * 2 2 2 2 * 1 ( / ) ( / ) ( / ) ( / ) 2 1 [ 2 2 ( )] 2 c c d u f x y xy x y xy dW dV dW dV dW dV dW dV v E         = − = − + + + − − (18) Fig. 4 shows the calculation process of the damage constitutive calculation model based on strain energy density. S.C. Li et alii, Frattura ed Integrità Strutturale, 47 (2019) 1-16; DOI: 10.3221/IGF-ESIS.47.01 8 Figure 4: Numerical simulation flow of rock failure based on strain energy density theory NUMERICAL SIMULATION OF ROCK FRACTURE Brazilian Splitting Numerical Simulation his case takes the Brazilian splitting with lateral concentrated load in document [17] as an example. Fig. 5 shows the calculation model. The disk diameter is 50mm and it is divided into 7860 elements. The model material adopts the marble in document [17] with the elastic modulus 66.8E GPa= , poisson's ratio 0.33v = , density 3 2620 /Kg m = T No , Yes/No Mark element failure. No No No Yes Yes Yes Yes , Yes/No Give residual modulus to the element. Does stepwise loading end? End calculation. Start calculation. Input model parameters. Impose load and balance FLAC program calculation. , Yes/No Calculate the element damage degree n and reduce modulus accordingly. Calculate the strain energy density of each element. S.C. Li et alii, Frattura ed Integrità Strutturale, 47 (2019) 1-16; DOI: 10.3221/IGF-ESIS.47.01 9 , element tension ultimate stress 20u MPa = , strain under ultimate stress 5 1.3 10u − =  , and strain when fracture occurs 4 1 10f − =  . The model top and bottom are imposed with load increased by 2KN in the negative y direction and positive y direction respectively. Figure 5: Numerical calculation model of Brazil splitting After the first 2KN load is imposed, no damage occurs inside the disk; each element is in the linear elastic stage; Fig. 5 shows the distribution of the stress inside the Brazilian disk. The results of comparison of distribution rules and theoretical results of stress distribution inside the disk in document [17] show that the simulated stress distribution rules basically agree well with the theoretical distribution. Simulated distribution Theoretical distribution yy xx xy Figure 6: Comparison between theoretical and numerical results of stress distribution in Brazil disc S.C. Li et alii, Frattura ed Integrità Strutturale, 47 (2019) 1-16; DOI: 10.3221/IGF-ESIS.47.01 10 As the load increases, the strain energy absorbed by each element gradually increases. The energy change process of elements in the middle of the disk top is recorded, as shown in Fig. 6. The strain energy density absorbed by the element ( / )dW dV gradually increases from zero while the critical strain energy density *( / )cdW dV gradually decreases from ( / ) 1000cdW dV = . When they two intersect, the element fails. Figure 7: Curve of element strain energy density versus load change. It can be seen from Fig. 7 that when calculation reaches the 7th loading step, that is, 14KN load is imposed, elements in the middle of the disk top have *( / ) ( / )cdW dV dW dV , indicating that elements at both ends of the disk start to fail. When calculation reaches the 15th loading step, that is, 30KN load is imposed, partial damage at both ends is suddenly developed into complete penetration in the middle of the specimen; meanwhile, there are partial failure areas on both ends of the disk. Fig. 8 shows the failure of the specimen when it is loaded to loading step 7, 12, 14, and 15. Step 7 Step 12 Step 14 Step 15 Figure 8: Element failure process of Brazil splitting under different loading steps Fig. 9 shows the damage status of each element when the specimen fails. The damage value n is expressed by the built-in variable ex_5 of the element and ranges from 1 to 20. It can be seen from the damage cloud chart of the element that the damage value of each element at the penetration point in the middle of the disk reaches the maximum value 20, indicating that elements in the middle of the disk are completely fractured and lose the bearing capacity. Elements on both ends of S.C. Li et alii, Frattura ed Integrità Strutturale, 47 (2019) 1-16; DOI: 10.3221/IGF-ESIS.47.01 11 the disk except a small number of elements are completely fractured, and surrounding elements also show different levels of damage. Figure 9: Element damage status in Brazilian disc in numerical calculation Documents [18-20] all conducted experimental studies on Brazilian disk imposed with lateral load. The results show that both ends of the disk firstly fail due to the stress concentration caused at the loading point during the splitting process of the brittle rock. As the load continues to increase, penetration occurs in the middle of the disk and the failure pie section appears on the loading points at both ends of the disk. Fig. 10 shows the sketch of splitting failure of the marble under different stepping bars in document [18]. Comparison of Fig. 8 and Fig. 10 shows that simulated results of numerical calculation in this paper agree well with test results in the document, indicating the feasibility of this method for simulating the rock fracture. Figure 10: Test results of Brazilian splitting for marble in document [18] Numerical Simulation of Tensile Test with Built-in Crack This case takes the tensile test with built-in crack in document [21] as an example. Fig. 11 shows the calculation model. The specimen size is 50 100mm mm ; the length of the built-in crack is 20mm ; the specimen is divided into 9772 elements. The specimen material adopts the rock mortar with the elastic modulus 5.17E GPa= , Poisson's ratio 0.192v = , density 32300 /Kg m = , element ultimate tensile stress 2.7u MPa = , strain under ultimate stress 4 5.8 10u − =  , and strain when fracture occurs 4 7 10f − =  . The model top and bottom are imposed with uniform tensile load increased by 100N respectively. Figure 11: Numerical calculation model of uniaxial tensile test with built-in crack. S.C. Li et alii, Frattura ed Integrità Strutturale, 47 (2019) 1-16; DOI: 10.3221/IGF-ESIS.47.01 12 The following calculation results are obtained: under the action of tensile load, stress concentration occurs on the crack tip; elements at both ends of the crack start to fail when calculation reaches the fourth loading step, that is, 400N load. Fig. 12 shows the strain energy density of the first element at the left end of the crack with changing load. Figure 12: Curve of element strain energy density versus load change When calculation reaches the fifth loading step, that is, 500N load, penetration occurs in the middle of the specimen in the crack direction. Fig. 13 shows the failure status of the model when calculation reaches loading steps 4 and 5. Step 4 Step 5 Figure 13: Element failure process of uniaxial tensile test under different loading steps Fig. 14 shows the damage status of each element when the specimen fails. It can be seen from the damage cloud chart of the element that each element surrounding the crack tip shows different levels of damage due to stress concentration and that the damage value of a series of elements from both ends of the crack to the model edge reaches 20, indicating that the specimen is completely fractured along the cross section of the crack. Figure 14: Element damage status in tensile test with built-in crack S.C. Li et alii, Frattura ed Integrità Strutturale, 47 (2019) 1-16; DOI: 10.3221/IGF-ESIS.47.01 13 Document [21] conducted uniaxial tensile tests on mortar materials with different prefabricated crack dips on the rigidity servo testing machine. The test results show that when the tensile stress reaches a specific value, the prefabricated crack starts to initiate and then rapidly expands, causing the structural failure of the specimen. When the crack is horizontal (that is, crack dip 90 = ), wing crack is generated at the front edge of the crack and expanded in the direction of the original crack until the specimen is fractured. Finally, the fracture trace is perpendicular to the axial tensile stress direction and the fracture surface is flat. Fig. 15 shows the tensile failure of horizontal crack in document [21]. Comparison of Fig. 13 and Fig. 15 shows that the simulated failure process and damage status in the built-in crack tensile tests of numerical calculation in this paper basically agrees with test results, indicating the accuracy and feasibility of this method in this paper for simulating the rock fracture. Figure 15: Test results of uniaxial tensile test in document [21] Numerical Simulation of Tunnel Excavation This case takes Liangshui Tunnel of Lanzhou-Chongqing Line as an example to simulate the overloading process of its similar model test. The horizontal direction of tunnel cross-section is x-axis, the vertical direction is y-axis, and the excavation direction is z-axis. The calculation range is as follows: 60 60x m−   ; 67.5 52.5y m−   ; the thickness of Z direction is 1m. The buried depth of the tunnel is 200 m and the lateral pressure coefficient is 0.54. The material parameters of surrounding rock are as follows: 2E GPa= , 0.35v = , 0.3c MPa= , 37 = , 1t MPa= , and 3 2620 /Kg m = . The damage judgement parameters are: ultimate tensile stress of element 3u MPa = , strain under the ultimate stress 31 10u − =  , and strain when fracture occurs 2 1 10f − =  . Fig. 16 shows the numerical calculation model. In order to simulate the tunnel buried depth of 200 meters, the vertical crustal stress at the buried depth of 147.5 meters and the corresponding horizontal stress are applied in the model. Figure 16: Numerical model of Liangshui tunnel FLAC3D 3.00 Itasca Consulting Group, Inc. Minneapolis, MN USA Settings: Model Perspective 17:36:11 Wed Apr 27 2011 Center: X: 0.000e+000 Y: -7.500e+000 Z: 5.000e-001 Rotation: X: 90.000 Y: 0.000 Z: 0.000 Dist: 3.918e+002 Mag.: 1.25 Ang.: 22.500 Block Group 1 2 3 4 S.C. Li et alii, Frattura ed Integrità Strutturale, 47 (2019) 1-16; DOI: 10.3221/IGF-ESIS.47.01 14 In the similar model test, overload test is carried out after tunnel excavation. Keeping the coefficient of lateral pressure unchanged, increase the buried depth of tunnel step by step with the gradient of 50 m buried depth. That is to say, vertical loads at a specific buried depth are applied on the top of the tunnel and horizontal loads at corresponding horizontal in- situ stress are applied on both sides of the tunnel to study the progressive failure characteristics of the tunnel. In the simulation calculation, the tunnel model is loaded in the same way. During the numerical simulation, the damage and failure of surrounding rock around the tunnel are as follows: (1) When loaded to 250 m buried depth, the damage and failure of surrounding rock is shown in Fig. 17. It can be seen that at this time, the surrounding rock has not been destroyed, and the surrounding rock in the tunnel top and bottom within about 1 times the diameter of the tunnel begins to damage, and a small part of the elements at the distance of about 1 times the diameter of the tunnel on both sides of the side wall also occurs damage. The maximum damage value of the element is 2, and it occurs at the vault and the inverted arch. Figure 17: Damage and failure contour of surrounding rock. (2) When loaded at 350 m burial depth, the elements closest to the free surface of the vault and inverted arch begin to fail, and the damage value of the element is the largest, reaching 16. The damage scope of surrounding rock at the top and bottom of the tunnel has expanded by two times, and the damage of the elements in the range of about 0.8 times of the tunnel diameter on both sides has also occurred. Fig. 18 shows the damage and failure of surrounding rock. Figure 18: Damage and failure contour of surrounding rock. (3) When loaded to the 600 m burial depth, the failure area of the vault and inverted arch continues to expand, extending to the tunnel spandrels, and some elements at the waist of the arch begin to fail. The damage value of the elements at the vault and inverted arch within about 1 times the diameter of the tunnel reaches the maximum value of 20, which indicates that they have completely broken and lost their bearing capacity, as shown in Fig. 19. Figure 19: Damage and failure contour of surrounding rock. FLAC3D 3.00 Itasca Consulting Group, Inc. Minneapolis, MN USA Step 7230 Model Perspective 10:50:26 Sat Apr 30 2011 Center: X: 0.000e+000 Y: -7.500e+000 Z: 5.000e-001 Rotation: X: 90.000 Y: 0.000 Z: 0.000 Dist: 3.918e+002 Mag.: 3.05 Ang.: 22.500 Block Contour of Zone Extra 2 0.0000e+000 to 1.0000e+000 1.0000e+000 to 2.0000e+000 2.0000e+000 to 2.0000e+000 Interval = 1.0e+000 FLAC3D 3.00 Itasca Consulting Group, Inc. Minneapolis, MN USA Step 7230 Model Perspective 10:51:44 Sat Apr 30 2011 Center: X: 0.000e+000 Y: -7.500e+000 Z: 5.000e-001 Rotation: X: 90.000 Y: 0.000 Z: 0.000 Dist: 3.918e+002 Mag.: 3.05 Ang.: 22.500 Block Group 4 FLAC3D 3.00 Itasca Consulting Group, Inc. Minneapolis, MN USA Step 19893 Model Perspective 21:05:41 Wed Mar 30 2011 Center: X: 0.000e+000 Y: -7.500e+000 Z: 5.000e-001 Rotation: X: 90.000 Y: 0.000 Z: 0.000 Dist: 3.918e+002 Mag.: 3.05 Ang.: 22.500 Block Contour of Zone Extra 2 0.0000e+000 to 1.0000e+000 2.0000e+000 to 3.0000e+000 4.0000e+000 to 5.0000e+000 6.0000e+000 to 7.0000e+000 8.0000e+000 to 9.0000e+000 1.0000e+001 to 1.1000e+001 1.2000e+001 to 1.3000e+001 1.4000e+001 to 1.5000e+001 1.6000e+001 to 1.6000e+001 Interval = 1.0e+000 FLAC3D 3.00 Itasca Consulting Group, Inc. Minneapolis, MN USA Step 19893 Model Perspective 21:03:44 Wed Mar 30 2011 Center: X: 0.000e+000 Y: -7.500e+000 Z: 5.000e-001 Rotation: X: 90.000 Y: 0.000 Z: 0.000 Dist: 3.918e+002 Mag.: 3.05 Ang.: 22.500 Block Group 4 fail FLAC3D 3.00 Itasca Consulting Group, Inc. Minneapolis, MN USA Step 30381 Model Perspective 21:27:26 Wed Mar 30 2011 Center: X: 0.000e+000 Y: -7.500e+000 Z: 5.000e-001 Rotation: X: 90.000 Y: 0.000 Z: 0.000 Dist: 3.918e+002 Mag.: 3.05 Ang.: 22.500 Block Contour of Zone Extra 2 0.0000e+000 to 1.0000e+000 2.0000e+000 to 3.0000e+000 4.0000e+000 to 5.0000e+000 6.0000e+000 to 7.0000e+000 8.0000e+000 to 9.0000e+000 1.0000e+001 to 1.1000e+001 1.2000e+001 to 1.3000e+001 1.4000e+001 to 1.5000e+001 1.6000e+001 to 1.7000e+001 1.8000e+001 to 1.9000e+001 2.0000e+001 to 2.0000e+001 Interval = 1.0e+000 FLAC3D 3.00 Itasca Consulting Group, Inc. Minneapolis, MN USA Step 30381 Model Perspective 21:26:27 Wed Mar 30 2011 Center: X: 0.000e+000 Y: -7.500e+000 Z: 5.000e-001 Rotation: X: 90.000 Y: 0.000 Z: 0.000 Dist: 3.918e+002 Mag.: 3.05 Ang.: 22.500 Block Group 4 fail S.C. Li et alii, Frattura ed Integrità Strutturale, 47 (2019) 1-16; DOI: 10.3221/IGF-ESIS.47.01 15 (4) When loaded to the 900 m burial depth, all the elements around the tunnel begin to fail, and the failed elements interconnect to form a whole region. Especially in the position of the vault and the inverted arch, a large area of elements has failed. The damage value of the elements around the tunnel reaches the maximum value of 20 except for a small part of the elements at the side wall, as shown in Fig. 20. Figure 20. Damage and failure contour of surrounding rock The loading test results of the scale model of Liangshui Tunnel are as follows: (1) When loaded to 350 m burial depth, the vault begins to crack and fall off. (2) When loaded to 600 m burial depth, the vault has larger deformation, and the side wall begins to crack gradually. (3) When loaded to 900 m burial depth, the vault material begins to peel off and collapse in large volume, the arch foot also appears obvious deformation, and the whole section of the tunnel has been seriously deformed. At this time, it is considered that the tunnel has been destabilized and destroyed. Fig. 21 shows the failure in the model overload test. Figure 21: Failure in the model overload test. The results of tunnel overload test and numerical simulation are basically the same. The only difference is that the failure occurs below the inverted arch during numerical simulation while no corresponding phenomenon is observed in the model test. This is due to the side friction effect in the model test on the one hand, and the failure of the inverted arch under the condition of excavation cannot be clearly observed on the other hand. CONCLUSIONS his numerical simulation method of rock fracture based on strain energy density theory completes the nonlinear calculation process with linear calculation, avoids singularity of numerical calculation in element fracture, and simulates the rock post-peak fracture behaviors. (1) The simple tensile bilinear strain softening constitutive model is adopted to analyze the energy and damage evolution process after the rock element reaches the stress limit point from the micromechanics perspective. In addition, the criterion is created to determine element failure based on strain energy density theory. This method can well simulate the entire process of rock damage. (2) The strain energy density of the rock element is obtained by accumulating the energy absorbed by the element after the load is imposed each time. Meanwhile, the modulus reduction of the damaged element is discretized to complete nonlinear calculation with linear calculation. A very small residual modulus is given to elements (regarded to be completely fractured) with the maximum damage values. Continuous calculation is used to implement discontinuous fracture behaviors of the rock to avoid singularity of numerical calculation in element fracture. FLAC3D 3.00 Itasca Consulting Group, Inc. Minneapolis, MN USA Step 38004 Model Perspective 21:23:49 Wed Mar 30 2011 Center: X: 0.000e+000 Y: -7.500e+000 Z: 5.000e-001 Rotation: X: 90.000 Y: 0.000 Z: 0.000 Dist: 3.918e+002 Mag.: 3.05 Ang.: 22.500 Block Contour of Zone Extra 2 0.0000e+000 to 1.0000e+000 2.0000e+000 to 3.0000e+000 4.0000e+000 to 5.0000e+000 6.0000e+000 to 7.0000e+000 8.0000e+000 to 9.0000e+000 1.0000e+001 to 1.1000e+001 1.2000e+001 to 1.3000e+001 1.4000e+001 to 1.5000e+001 1.6000e+001 to 1.7000e+001 1.8000e+001 to 1.9000e+001 2.0000e+001 to 2.0000e+001 Interval = 1.0e+000 FLAC3D 3.00 Itasca Consulting Group, Inc. Minneapolis, MN USA Step 38004 Model Perspective 21:22:57 Wed Mar 30 2011 Center: X: 0.000e+000 Y: -7.500e+000 Z: 5.000e-001 Rotation: X: 90.000 Y: 0.000 Z: 0.000 Dist: 3.918e+002 Mag.: 3.05 Ang.: 22.500 Block Group 4 fail T S.C. Li et alii, Frattura ed Integrità Strutturale, 47 (2019) 1-16; DOI: 10.3221/IGF-ESIS.47.01 16 (3) Numerical simulation of Brazilian splitting, tensile tests with built-in crack and tunnel excavation shows that the preceding method can well simulate the entire process from partial failure to structural fracture of the rock and obtain the damage evolution process of each element during the failure process. (4) Studying rock fracture from the energy perspective can comprehensively reflect the rock strength changes and essential characteristics of the structural failure. This paper develops a rock damage numerical calculation program with Fish language in the FLAC. Compared with the theory and analytic method, this program can directly resolve various complex problems and therefore has a wider application prospect. ACKNOWLEDGMENTS he authors would like to thank the support of the National Basic Research Program of China (Grant No. 2010CB732002), the National Natural Science Foundation of China (Grant No. 51179098、51379113), the Specialized Research Fund for the Doctoral Program of Higher Education of China (Grant No. 20120131110031). REFERENCES [1] Cai, M. (2002). Rock Mechanics and Engineering, Beijing: Science Press. [2] Xie, H., Ju, Y. and Li, L. (2005). Criteria for strength and structural failure of rocks based on energy dissipation and energy release principles, Chinese J. of Rock Mech. and Eng., 24(17), pp. 3003-3010. [3] Mikhalyuk, A.V. and Zakharov, V.V. (1997). Dissipation of dynamic-loading energy in quasi-elastic deformation processes in rocks, Appl. Mech. and Tech. Phys., 38(2), pp. 312-318. [4] Steffler, E.D., Epstein, J.S. and Conley, E.G. (2003). Energy partitioning for a crack under remote shear and compression, Int. J. Fracture, 120, pp. 563-580. [5] You, M., Hua, A. (2012). Energy analysis on failure process of rock specimens, Chinese J. of Rock Mech. and Eng., 21(6), pp.778-781. [6] Xie, H., Peng, R. and Ju, Y. (2014). Energy dissipation of rock deformation and fracture, Chinese J. of Rock Mech. and Eng., 23(21), pp. 3565–3570. [7] Yang, S., Xu, W. and Su C. (2006). Study on the deformation failure and energy properties of rock specimen in uniaxial compression, Acta Mechanica Solida Sinica, 27(2), pp. 213-216. [8] Carloni, C. and Nobile, L. (2002). Crack initiation behavior of orthotropic solids as predicted by the strain energy density theory, Theoretical and Applied Fracture Mechanics, (38), pp.109-119. [9] Liu, H. (1999). Mechanics of Materials (I), Beijing: High Education Press. [10] Yu, X. (2012). The revision of the fourth strength theorem, J.of Guangxi Uni.(Natural Science), 27(Supp.), pp.63-69. [11] Li, Q.M. (2011). Strain energy density failure criterion, Int. J. of Solids and Structures, (38), pp. 6997-7013. [12] Gdoutos, E.E. (1984). Problems of Mixed Mode Crack Propagation, Martinus Nijhoff Publ, [13] Wang, X., Hai, L. and Huang, M. (2014). Inhomogeneity analysis of damage under uniaxial tension (I) - basic theory, Chinese J. of Rock Mech. and Eng., 23(9), pp. 1446–1449. [14] Jeong, D.Y., Orringer, O. and Sih, G.C. (1994). Strain energy density approach to stable crack extension under net section yielding of aircraft fuselage, Theoretical and Appl. Fracture Mech., 22(2), pp.127-137. [15] Jin, F., Jiang, M. and Gao, X. (2014). Defining damage variable based on energy dissipation, Chinese J. of Rock Mech. and Eng., 23(12), pp.1976–1980. [16] Li, Z. (2012). Damage Mechanics and Its Application, Beijing: Science Press, pp.16–19. [17] Ye, J., Wu, F.Q. and Sun, J.Z. (2013). Estimation of the tensile elastic modulus using Brazilian disc by applying diametrically opposed concentrated loads, Int. J. of Rock Mech. & Mining Sci., 46, pp. 568-576. [18] Yang, T., Wang, B. and Sun, L. (2012). Effects of various spacer methods for rock split tests, Site Investigation Sci. and Tech., 1, pp.3-7. [19] Hudson, J.A., Brown, E.T. and Rummel, F. (1972). The controlled failure of rock discs and rings loaded in diametral compression, Int. J. Rock Mech. Min. Sci., 9, pp. 241–248. [20] Markides, Ch.F., Pazis, D.N. and Kourkoulis, S.K. (2010). Closed full-field solutions for stresses and displacements in the Brazilian disk under distributed radial load, Int. J. of Rock Mech. & Mining Sci., 47, pp.227-237. [21] Yang, L. (2012). Experimental and Numerical Simulation Study on the Tensile Fracture Mechanism of Three- dimensional Flaws in Fractured Rock Mass [Master Thesis], Ji Nan: Shandong University. T