Mathematical and Numerical Modeling in Geotechnical EngineeringView this Special Issue
Fracture Analysis of Brittle Materials Based on Nonlinear FEM and Application in Arch Dam with Fractures
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.
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  and continuum damage mechanics . Some researchers combine both methods for numerical failure analysis [11, 12]. Crack propagation simulation could be concluded as smeared fracture model  and discrete fracture model . Besides, numerical tools have developed rapidly, including Finite Element Method (FEM) , Extended FEM , Element Free Method , Boundary Element Method , Discrete Element Method , Numerical Manifold Method , Discontinuous Deformation Analysis , and Fast Lagrangian Method .
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  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 , 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;
Equation (4) is known as the closest-point projection method (CPPM) , 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.
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 :
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.
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.  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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.View at: Google Scholar
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.View at: Google Scholar
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.View at: Google Scholar
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.View at: Google Scholar
G. R. Irwin, “Fracture dynamics,” in Fracturing of Metals, pp. 147–166, 1948.View at: Google Scholar
E. Orowan, “Fracture and strength of solids,” Reports on Progress in Physics, vol. 12, no. 1, pp. 185–232, 1948.View at: Google Scholar
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.View at: Google Scholar
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.View at: Google Scholar
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.View at: Google Scholar
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.View at: Google Scholar
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.View at: Google Scholar
G. H. Shi, “Manifold method,” in Proceedings of the 1st International Forum on DDA and Simulations of Discontinuous Media, Bekerley, Calif, USA, 1996.View at: Google Scholar
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.View at: Publisher Site | Google Scholar
R. Hill, The Mathematical Theory of Plasticity, Clarendon Press, Oxford, UK, 1950.View at: MathSciNet
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.View at: Publisher Site | Google Scholar | MathSciNet
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.View at: Google Scholar
L. M. Kachanov, “Time of the rupture process under creep conditions,” Izvestiya Akademii Nauk SSR Otdelenie Tekhnicheskikh Nauk, vol. 8, pp. 26–31, 1958.View at: Google Scholar
I. N. Rabotnov, Creep Problems in Structural Members, vol. 7, North-Holland, Amsterdam, the Netherlands, 1969.