- Email: [email protected]
Finite Element Modelling of Complex 3D Static and Dynamic Crack Propagation by Embedding Cohesive Elements in Abaqus
Acta Mechanica Solida Sinica, Vol. 23, No. 3, June, 2010 Published by AMSS Press, Wuhan, China
ISSN 0894-9166
FINITE ELEMENT MODELLING OF COMPLEX 3D STATIC AND DYNAMIC CRACK PROPAGATION BY EMBEDDING COHESIVE ELEMENTS IN ABAQUS Xiangting Su1,2
Zhenjun Yang2
Guohua Liu1
1
( College of Civil Engineering and Architecture, Zhejiang University, Hangzhou 310027, China) (2 Department of Engineering, the University of Liverpool, L69 3GQ, UK)
Received 23 November 2009; revision received 25 May 2010
ABSTRACT This study proposes an algorithm of embedding cohesive elements in Abaqus and develops the computer code to model 3D complex crack propagation in quasi-brittle materials in a relatively easy and efficient manner. The cohesive elements with softening traction-separation relations and damage initiation and evolution laws are embedded between solid elements in regions of interest in the initial mesh to model potential cracks. The initial mesh can consist of tetrahedrons, wedges, bricks or a mixture of these elements. Neither remeshing nor objective crack propagation criteria are needed. Four examples of concrete specimens, including a wedgesplitting test, a notched beam under torsion, a pull-out test of an anchored cylinder and a notched beam under impact, were modelled and analysed. The simulated crack propagation processes and load-displacement curves agreed well with test results or other numerical simulations for all the examples using initial meshes with reasonable densities. Making use of Abaqus’s rich pre/postprocessing functionalities and powerful standard/explicit solvers, the developed method offers a practical tool for engineering analysts to model complex 3D fracture problems.
KEY WORDS finite element method, cohesive elements, three-dimensional crack propagation, discrete crack model, concrete structures, Abaqus
I. INTRODUCTION Crack propagation is an inherent feature in quasi-brittle materials such as concrete due to the low material tensile strength, and it is a predominant source of nonlinearity and a main culprit for ultimate failure of structures made from these materials. Accurate understanding of the crack propagation behaviour so as to assess the structural load-carrying capacity, by either experimental studies or numerical modelling, remains a worldwide challenge, especially for three-dimensional (3D) problems with complex geometries and loading conditions. Numerical modelling of crack propagation has been an active research field since 1960s[1, 2]. Nowadays there exist a large number of numerical models. In terms of how the cracks are modelled geometrically, there are discrete crack models explicitly separating crack surfaces and modelling discontinuity, smeared crack models based on continuum mechanics, and more indirectly, lattice models. In terms of numerical methods, both the traditional finite element methods (FEM) and boundary element methods (BEM)
Corresponding author. Email: [email protected] Project supported by EPSRC UK (No. EP/F00656X/1). Xiangting Su’s one-year visit to the University of Liverpool was supported by the China Scholarship Council and the National Natural Science Foundation of China (No. 50579081).
· 272 ·
ACTA MECHANICA SOLIDA SINICA
2010
have been widely applied to crack propagation modelling. Extensive literature review of these models and methods in 2D problems has been given elsewhere[3–6] . It appears now that the discrete crack models, mostly based on the cohesive crack model (CCM) developed in ductile materials by Barenblatt[7] and Dugdale[8] and in quasi-brittle materials by Hillerborg[9] , are becoming more and more popular, because of their ability to model macroscopic cracks with strong discontinuity, the CCM’s capability of realistically representing the energy dissipation during fracture processes, and the ease of its implementation as cohesive interface elements (CIE) in the FEM and the BEM. In general, two types of approaches are used to model discrete cohesive crack propagation. The first type of approaches is based on sophisticated remeshing procedures that constantly change the meshes as cracks propagate[3, 10–12] . For problems with crack paths unknown a priori, objective crack propagation criteria are needed to judge when and in which direction a crack propagates. This usually involves calculation of stress intensity factors (SIFs) and stresses at crack tips, whose accuracy can only be ensured by fine crack-tip meshes or using singular elements. This in turn exacerbates the complexity of the remeshing procedure. It is also very difficult for the remeshing-based approaches to deal with complex crack propagation situations, such as multiple cracking and 3D problems. Indeed, the success of the remeshing-based approaches has been largely limited to 2D static fracture problems with single or a few cracks, although they are computationally efficient because a relatively small number of nonlinear CIEs are inserted into the meshes. Another type of approaches pre-insert or pre-embed CIEs between the finite elements in the initial meshes[6, 13–15] . Crack propagation is modelled by automatic opening and merging of the CIEs under applied loadings. Neither remeshing procedures nor crack propagation criteria are needed so complex crack propagation can be modelled. However, these approaches restrict the cracks to the finite element edges or surfaces and the predicted crack patterns may thus be mesh-dependent. Another disadvantage is the high computational cost due to the use of a large number of nonlinear CIEs. The CIEs can also be dynamically inserted when certain crack initiation criterion is satisfied[16, 17] so that the number of CIEs changes with the loading and the total number of CIEs remains relatively low. These procedures of dynamically inserting CIEs need objective crack initiation criteria and change meshes although the remeshing operation is relatively straightforward. In the last two decades, the partition of unity finite element method (PUFEM)[18, 19] and the extended finite element method (XFEM)[20–22] are becoming popular in modelling crack propagation. These methods introduce discontinuity into the finite elements. The crack growth is described independently of the mesh, so that remeshing is not needed. Meshless methods are also receiving much attention[23, 24] . These innovative methods, showing high potential in crack propagation modelling, generally need fine crack tip meshes like the traditional FEM to calculate accurate stresses or SIFs used in evaluating the crack propagation criterion. In this sense, fine initial meshes are needed if the crack paths are unknown in priori, which means high computational cost. The newly-developed semi-analytical scaled boundary finite element[25] is able to compute accurate SIFs without using fine meshes or singular elements. It has recently been applied to model discrete cohesive crack propagation[4, 5, 26], but so far only 2D problems were modelled. In summary, all the above crack models and numerical methods have their own advantages and disadvantages in terms of effectiveness, efficiency and applicability, and the selection largely depends upon the problems analysed and the preference and experience of the analysts. Most of the existing studies have modelled 2D crack propagation. 3D modelling of complex cohesive crack propagation has been sporadically conducted, and there is a long way to go before it becomes a routine exercise in research and practical design due to the inherent difficulties. In addition, most of the studies have used specialised in-house computer programs (e.g., Cornell University’s FRANC) that may be difficult for practical engineers to use. Most of the general-purpose commercial FEA packages, such as Ansys and Abaqus, are unable to model complicated crack propagation without user intervention or second development. This study aims at developing a practical finite element method for modelling complex 3D crack propagation in quasi-brittle materials under static and dynamic loadings. The method makes use of a special type of CIEs called ‘cohesive elements’, designed for modelling bonded interfaces in Abaqus V6.5 or higher[27] , to model potential cracks. An efficient algorithm was devised to insert these CIEs into concerned regions in FE meshes consisting of various types of solid elements. The method is simple to implement. The rich pre/post-processing functionalities and powerful implicit and explicit solvers
Vol. 23, No. 3
Xiangting Su et al.: Complex 3D Static and Dynamic Crack Propagation
· 273 ·
of Abaqus can be fully exploited. Its effectiveness is validated by modelling various fracture problems of concrete structures.
II. THE MODELLING METHOD The proposed method involves the following steps: (1) Meshing the domain using Abaqus/CAE and generating an input file. (2) Inserting cohesive elements into the regions of interest in the initial FE mesh using an in-house computer program and generating another input file. The initial mesh may consist of solid elements such as 4-noded tetrahedrons (C3D4), 6-noded wedges (C3D6) and 8-noded bricks (C3D8), or a mixture of these elements. Both 6-noded and 8-noded cohesive elements (COH3D6 and COH3D8) can be inserted. Figure 1 shows a mesh with solid elements only, the mesh after CIEs are inserted and the 3D CIEs, respectively. (3) Solving the problem using Abaqus standard or explicit solvers. This procedure can be readily automated by running a batch file. The following briefly discusses the key elements in the method.
Fig. 1. FE meshes before and after cohesive elements are inserted and 3D cohesive elements.
2.1. Cohesive Elements with Damage in Abaqus The cohesive crack model[7–9] assumes the existence of a fracture process zone (FPZ) in front of the real crack tip, in which energy dissipation occurs during fracture. In the FPZ, there exist tractions in the normal direction (tn ) and the two tangential or shear directions (ts and tt ) across the crack surfaces, resulting from mechanisms such as material bonding, aggregate interlocking and surface friction. Figure 2 shows a curve relating the normal traction tn and the crack opening displacement δn as an example. Similar traction-separation curves can be defined for shear tractions (ts and tt ) - crack sliding displacements Fig. 2 An example of traction-separation curve with softening (δs and δt ) relations. Before the crack initiates, a for cohesive elements. linear elastic ascending phase is assumed to model the initially un-cracked material. After the crack initiates, the traction decreases monotonically as functions of the corresponding separation, which is often termed tension or strain softening. The initial tensile stiffness kn0 and the initial shear stiffness ks0 and kt0 should be high enough to represent the un-cracked material, but not too high to cause numerical ill-conditioning. These initial stiffness values are determined by a trial and error approach. If δn is negative during loading increments or iterations, a compressive stiffness of magnitude equal to kn0 is assigned in order to prevent penetration of crack surfaces.
· 274 ·
ACTA MECHANICA SOLIDA SINICA
2010
The cohesive elements in Abaqus[27] are based on the cohesive crack model. The constitutive response of cohesive elements, defined in terms of traction-separation laws, assumes initially linear elastic behaviour followed by the initiation and evolution of damage. Several damage initiation criteria are available in Abaqus. This study used the quadratic nominal stress criterion in all the examples. The damage evolution criteria can be classified into three categories according to the relation between traction and separation: linear damage evolution, exponential damage evolution, and tabular damage evolution. All these criteria were used in this study. 2.2. Inserting Cohesive Elements into 3D FE Meshes Although inserting cohesive elements into 2D FE meshes is relatively straightforward[6], it is not a trivial task for 3D problems, especially if the FE mesh is not regular and consists of different types of solid elements. One challenge is how to robustly deal with the changes in the complicated 3D nodal and elemental connectivity due to the insertion of CIEs. This is tackled by an effective embedding algorithm implemented in a MATLAB[28] program with carefully designed data structures listed in Table 1. It has the following steps: Table 1. Data structures used in the algorithm of inserting cohesive elements in 3D meshes
Structure Name
Members
NODES index elemconn faceconn xcord ycord zcord isInDefinedArea
ELEMS index faceconn nodeconn type
FACES index nodeconn elemconn isBound isInDefinedArea
BALLS POINT isInDefinedArea
POINT elemindex xcord ycord zcord newnodeindex
(1) Reading the nodal coordinates and the nodal connectivity of the solid elements from the Abaqus input file. These data are then used to generate three arrays of type Structure: NODES, ELEMS and FACES. Most of the members in these structures are self-explanatory. For example, the members elemconn and faceconn in NODES structure are arrays containing the elements and faces the node is connected with, respectively. It is worth noting that only those nodes, elements and surfaces in the concerned or defined regions in the mesh are considered, in order to reduce the computational cost. This is indicated by the member isIndefinedArea. A consistent nodal ordering approach (clockwise or anti-clockwise) in nodeconn of FACES should be used to facilitate determination of correct nodal order when cohesive elements are created in the 4th step. (2) Inserting new nodes. For each existing node with isIndefinedArea=1, all the N number of solid elements connected with this node are first identified. N number of new nodes, each connected with one solid element, are then generated with the same coordinates of the existing node. All the new nodes are stored in an array of POINT, which is a member of a Structure BALLS associated with the existing node. The member elemindex in POINT represents the solid element the new node is connected with. An array of BALLS is then generated for all the existing nodes. The design of BALLS greatly facilitates subsequent generation of cohesive elements and is at the heart of the embedding algorithm. (3) Updating the arrays NODES and ELEMS. Since a large number of new nodes are generated, the elemental connectivity (elemconn) of NODES and the nodal connectivity (nodeconn) of ELMS are updated. It should be noted that the array FACES is not updated at this step. (4) Generating cohesive elements. This is based on the array FACES. There are two solid elements (elemconn) connected with each face. For each vertex of a face (nodeconn), a pair of new nodes connected to the two solid elements is found in POINT of BALLS. The three or four pairs of nodes found for this face comprise a cohesive element (COH3D6 or COH3D8). As the nodal order in FACES is carefully designed in the 1st step, the requirement on the nodal order of cohesive elements is readily satisfied. (5) Creating the Abaqus input file. The file contains the index and coordinates of POINT, the index and nodal connectivities of the updated ELEMS and the cohesive elements. To minimise the mesh-dependence for problems with crack paths unknown a priori, a reasonably fine initial mesh should be used. This means generally a large number of solid elements, nodes, surfaces
Vol. 23, No. 3
Xiangting Su et al.: Complex 3D Static and Dynamic Crack Propagation
· 275 ·
and CIEs are used. Computer implementation of the above inserting algorithm must be carried out with care, or the time spent on inserting CIEs may be unduly long. As an example, for an initial mesh with 46800 C3D8 elements and 52140 nodes, it takes about 30 minutes for the developed MATLAB program to insert 135,200 CIEs using a PC with an Intel Xeon [email protected] GHz and 3 GB physical memory.
III. NUMERICAL EXAMPLES Four 3D concrete structures subjected to static and dynamic loadings were modelled using the developed method as examples. A PC with two Intel Xeon CPUs @3.16GHz and 3 GB physical memory was used. Both Abaqus/standard and explicit solvers were used to solve the nonlinear equation systems. When the quasi-static explicit solver was used to model static problems, the analysis time was carefully chosen to avoid dynamic effects. Because the failure mode of all the examples is mainly tensile fracture and there is no experimental data about shear fracture resistance, the strength and the initial stiffness in shear are assumed as 10 times those in tension for the cohesive elements so that shear fracture failure does not occur. 3.1. Wedge-Splitting Test The first example is the wedge-splitting test of a concrete specimen CP250 carried out by Trunk[29] . The test was modelled in 2D by Trunk[29] and Feist[30] , and in 3D by Areias and Belytschko[22]. The geometry and the loading condition are shown in Fig.3. The concrete has a Young’s modulus E = 28300 MPa, the Poisson ratio ν = 0.2 and mass density ρ = 2.5 × 10−6 kg/mm3 . A bilinear softening curve shown in Fig.4 with parameters obtained from the experiment was used to model the normal tractionseparation relation of the CIEs. The area under the curve is the fracture energy Gf = 0.49 N/mm. The initial stiffness kn0 was assumed as 1 × 105 MPa·mm. The Abaqus/explicit quasi-static solver used an analysis time of 0.1 s.
Fig. 4. Bilinear softening law used for the wedgesplitting test. Fig. 3. Wedge-splitting test specimen (unit: mm).
Figure 5 shows a deformed FE mesh with 774 C3D6 and 1392 C3D8R solid elements, 1854 CIEs and 6852 nodes in total. It can be seen that the crack path was vertical as expected because of the symmetry of the mesh and loading conditions. Finer meshes were also modelled but they led to little difference from the predicted crack path and the load-crack mouth opening displacement (CMOD) curves. Figure 6 compares the load-CMOD curves from the experiment and modelling using different solvers and solid elements. The standard solver always failed by divergence at the post-peak stage, e.g., it failed at CMOD=3.2 mm for C3D8 as shown in Fig.6. This also often occurred in 2D simulations[6] . In contrast, the explicit solver always managed to model the whole fracture process when small enough time increments were used. The explicit solver using C3D8R with reduced integration and hourglass control (HC) resulted in nearly the same pre-peak responses as the standard solver using C3D8, while the explicit solver using C3D8R but without hourglass control led to too soft pre-peak responses. C3D8 was also tried with the explicit solver but the computational time nearly doubled when C3D8R was
· 276 ·
ACTA MECHANICA SOLIDA SINICA
2010
used. Therefore, the explicit solver was used in all the following simulations and C3D8R elements with hourglass control were chosen when 8-noded brick elements were used.
Fig. 5. Deformed mesh (CMOD=4 mm, scale=100).
Fig. 6. Load-CMOD curves: effect of solid element types and solvers.
Figure 7 compares the load-CMOD curves from 3D and 2D simulations. In the 2D simulation, a similar mesh to the in-plane mesh of the 3D mesh (Fig.5) was used assuming a plane stress condition. It can be seen that the peak load predicted from the 3D simulation is higher than that from the 2D, which reflects the constraining effect of the specimen thickness. A very similar mesh to Fig.5 but with a non-symmetric distribution of elements along the vertical central line was also modelled. The deformed mesh is shown in Fig.8. The crack path deviated from the central line, Fig. 7 Load-CMOD curves from 2D and 3D. indicating the effect of the initial mesh. The loadCMOD curves predicted from the two meshes are shown in Fig.9.
Fig. 8. Deformed mesh from a non-symmetric mesh (CMOD=4 mm, scale=100).
Fig. 9. Load-CMOD curves: effect of mesh arrangement.
3.2. Torsion Fracture Test The second example is a concrete bar under torsion tested by Brokenshire[31]. The geometry and the loading condition are illustrated in Fig.10. 3D modelling was carried out by Jefferson et al.[32] and Gasser and Holzapfel[33] . The three support points were assumed as fixed and the loading incremental procedure was controlled by the vertical displacement at the point C. The concrete has E = 35000 MPa, ν = 0.2 and
Vol. 23, No. 3
Xiangting Su et al.: Complex 3D Static and Dynamic Crack Propagation
· 277 ·
Fig. 10. A concrete bar under torsion (unit: mm).
ρ = 2.5 × 10−6 kg/mm3 . The same linear softening curve (Fig.11) as used in Ref.[32] was adopted to model the normal traction-separation relation with Gf = 0.08 N/mm. The initial stiffness was set as kn0 = 7.5 × 106 MPa·mm. An analysis time of 0.3 s was used in the explicit quasi-static procedure. CMOD is defined as the relative displacement between the points A and B normal to the Fig. 11 Linear softening law used for the torsion test. notch surface (ref. Fig.10). Two FE meshes with different solid element types were modelled to investigate the effect on the prediction of crack surfaces. Mesh A in Fig.12(a) has 10651 C3D4 solid elements and 9606 CIEs, and Mesh B in Fig.12(b) has 7027 C3D8R and C3D6 elements and 8055 CIEs. The predicted crack surfaces are compared with the test result in Fig.13. It is clear that Mesh A with C3D4 led to very realistic, tortuous fracture surfaces whereas Mesh B with C3D8R predicted an unrealistically flat and smooth crack path. The predicted load-CMOD curves are presented in Fig.14. It can be seen that both meshes led to curves close to the experimental data. Mesh A predicted a peak load slightly higher than Mesh B. This may be explained by the predicted crack surfaces in Fig.13, where the tortuous and unsmooth crack surfaces from Mesh A offer higher fracture resistance than the flat and smooth ones from Mesh B. Jefferson et al.[32] also modelled the same test using a similar mesh to Mesh B and the ‘Craft’ concrete
Fig. 12. Two FE meshes used for the torsion test.
· 278 ·
ACTA MECHANICA SOLIDA SINICA
2010
Fig. 13. Crack surfaces after failure for the torsion test.
Fig. 14. Load-CMOD curves for the torsion test.
Fig. 15. Pull-out test of a steel bar embedded in concrete (unit: mm).
model. The predicted load-CMOD curve is also shown in Fig.14 and it is very close to that from Mesh B in this study. 3.3. Pull-out Test of a Steel-Embedded Concrete Anchor This example models the pull-out behaviour of a steel anchor embedded in concrete. The detailed geometry, applied boundary and loading conditions of one-quarter specimen are illustrated in Fig.15. This example was also modelled in Ref.[19,22,23] using other methods. The anchor bolt, excluded from the FE model, was simulated by a boundary condition of vertical displacement imposed on the top
Fig. 16. Meshes used for the pull-out test.
Vol. 23, No. 3
Xiangting Su et al.: Complex 3D Static and Dynamic Crack Propagation
· 279 ·
surface of the groove. The concrete properties are: E = 30000 MPa, ν = 0.2 and ρ = 2.5×10−6 kg/mm3 . A linear softening curve with Gf = 0.1 N/mm, tensile strength=3 MPa and kn0 = 6 × 104 MPa·mm was used. An analysis time of 0.1 s was used in the explicit solver. To investigate the mesh-dependence of results, two meshes were modelled. The coarse mesh in Fig.16(a) has 31889 C3D4 solid elements and 60164 CIEs, and the numbers for the fine mesh in Fig.16(b) are 69798 and 138169, respectively. Figure 17 shows the load (F )-displacement (d) curves from the two meshes and other studies. It can be seen that the two meshes predicted very close responses, especially the same failure point when the load suddenly dropped at d = 0.25 mm.
Fig. 17. Load-displacement curves for the pull-out test.
Similar crack propagation processes and final crack patterns were predicted from the two meshes. Figure 18 illustrates four stages of the predicted crack propagation using the fine mesh, where the cracks are represented by grey areas consisting of highly-damaged cohesive elements (with the damage
Fig. 18. Predicted crack propagation process of the pull-out test.
· 280 ·
ACTA MECHANICA SOLIDA SINICA
2010
index D ≥ 0.99). At F = 310 kN (60% of the peak load Fu = 503 kN), tensile cracks initiated at the top edge of the groove, forming a cone around the groove (Fig.18(a)). As the load approached the peak, a splitting-mode longitudinal crack initiated at the top surface and propagated in the radial direction (Fig.18(b)). Beyond the peak load (Fig.18(c)), the longitudinal crack propagated rapidly, leading to a drastic decline on the load-displacement curve (Fig.17). After the complete propagation of the longitudinal crack (Fig.18(d)), the conical crack continued to advance with the load decreasing slowly. The simultaneous propagation of the longitudinal crack and the conical crack, observed in the experiments by Rots[34] , was well reproduced. This may explain the difference between the predicted load-displacement curves (Fig.17) by the present study and other studies[22, 23] which only modelled the conical crack. 3.4. Impact Test of a Concrete Beam The fourth example is a concrete beam under impact tested and modelled in 2D by Du et al.[35] . 3D modelling was reported by Belytschko et al.[36]. The geometry and the boundary condition are shown in Fig.19. The impact load history obtained from the test (Fig.20) was input as the loading condition. The concrete properties are: E = 34480 Fig. 19 Impact test of a concrete beam (unit: mm). MPa, ν = 0.2 and ρ = 2.5 × 10−6 kg/mm3 . An exponential softening curve (Fig.21) with the exponential law parameter α = 1.5 and the equivalent dynamic fracture energy Gf = 0.152 N/mm, suggested in Ref.[35], was adopted to model the CIEs. The initial stiffness kn0 = 1 × 105 MPa·mm was used. Figure 22 shows a deformed mesh at the post-peak stage when t = 1 ms. The mesh consists of 3645 solid elements and 606 cohesive elements. Figure 23 compares favourably the predicted loaddeflection curve with the test data. Figure 24 compares the measured and computed strain histories at three locations (SG01, SG02 and SG03). The strain history computed for SG03 agreed well with the measured. For the other two locations, the strain histories were underestimated, probably due to
Fig. 21. Exponential softening law for the impact test. Fig. 20. Load history used for the impact test.
Fig. 22. Deformed mesh for the impact test.
Fig. 23. Load-displacement curves for the impact test.
Vol. 23, No. 3
Xiangting Su et al.: Complex 3D Static and Dynamic Crack Propagation
· 281 ·
the negligence of shrinkage effects in the numerical modelling[37] . The crack-tip extension histories are compared in Fig.25.
Fig. 24. Strain history for the impact test.
Fig. 25. Crack-tip extension history for the impact test.
IV. CONCLUSIONS A simple yet effective finite element methodology has been developed to simulate 3D complex cohesive crack propagation in quasi-brittle materials, using Abaqus’s cohesive elements. An efficient algorithm to insert the cohesive elements into initial finite element meshes was devised and implemented by an in-house computer program. Several concrete examples were modelled to demonstrate the effectiveness of the method. It is concluded that, making use of the rich pre/post-processing functionalities and powerful standard/explicit solvers in Abaqus, the developed methodology seems promising to provide a practical simulation tool for modelling realistic 3D fracture in various industries.
References [1] Scordelis,A.C. and Ngo,D., Finite element analysis of reinforced concrete beams. Journal of the American Concrete Institute , 1967, 64: 152-163. [2] Rashid,Y.R., Ultimate strength analysis of prestressed concrete pressure vessels. Nuclear Engineering and Design, 1968, 7(4): 334-344. [3] Yang,Z.J. and Chen,J.F., Fully automatic modelling of cohesive discrete crack propagation in concrete beams using local arc-length methods. International Journal of Solids and Structures, 2004, 41(3-4): 801826. [4] Yang,Z.J., Fully automatic modelling of mixed-mode crack propagation using scaled boundary finite element method. Engineering Fracture Mechanics, 2006, 73(12): 1711-1731. [5] Yang,Z.J. and Deeks,A.J., Fully-automatic modelling of cohesive crack growth using a finite element-scaled boundary finite element coupled method. Engineering Fracture Mechanics, 2007, 74(16): 2547-2573. [6] Yang,Z.J., Su,X.T., Chen,J.F. and Liu,G.H., Monte Carlo simulation of complex cohesive fracture in random heterogeneous quasi-brittle materials. International Journal of Solids and Structures, 2009, 46(17): 32223234. [7] Barenblatt,G.I., The formation of equilibrium cracks during brittle fracture: general ideas and hypothesis, axially symmetric cracks. Applied Mathematics and Mechanics, 1959, 23: 622-636. [8] Dugdale,D.S., Yielding of steel sheets containing slits. Journal of Mechanics of Physics and Solids, 1960, 8: 100-104. [9] Hillerborg,A., Modeer,M. and Petersson,P., Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cement and Concrete Research, 1976, 6: 773-782. [10] Wawrzynek,P.A. and Ingraffea,A.R., An interactive approach to local remeshing around a propagating crack. Finite Elements in Analysis and Design, 1989, 5: 87-96. [11] Xie,M. and Gerstle,W.H., Energy-Based Cohesive Crack Propagation Modeling. Journal of Engineering Mechanics-ASCE, 1995, 121(12): 1349-1358. [12] Yang,Z.J. and Xu,X.F., A heterogeneous cohesive model for quasi-brittle materials considering spatially varying random fracture properties. Computer Methods in Applied Mechanics and Engineering, 2008, 197: 4027-4039.
· 282 ·
ACTA MECHANICA SOLIDA SINICA
2010
[13] Xu,X.P. and Needleman,A., Numerical simulations of fast crack growth in brittle solids. Journal of the Mechanics and Physics of Solids, 1994, 42(9): 1397-1434. [14] L´ opez,C., Carol,I. and Aguado,A., Meso-structural study of concrete fracture using interface elements. I: numerical model and tensile behavior. Materials and Structures, 2008, 41(3): 583-599. [15] L´ opez,C., Carol,I. and Aguado,A., Meso-structural study of concrete fracture using interface elements. II: compression, biaxial and Brazilian test. Materials and Structures, 2008, 41(3): 601-620. [16] Camacho,G.T. and Ortiz,M., Computational modelling of impact damage in brittle materials. International Journal of Solids and Structures, 1996, 33: 2899-2938. [17] Pandolfi,A. and Ortiz,M., An efficient adaptive procedure for three-dimensional fragmentation simulations. Engineering with Computers, 2002, 18: 148-159. [18] Melenk,J.M. and Babuska,I., The partition of unity finite element method: Basic theory and applications. Computer Methods in Applied Mechanics and Engineering, 1996, 139: 289-314. [19] Gasser,T.C. and Holzapfel,G.A., Modeling 3D crack propagation in unreinforced concrete using PUFEM. Computer Methods in Applied Mechanics and Engineering, 2005, 194: 2859-2896. [20] Moes,N., Dolbow,J. and Belytschko,T., A finite element method for crack growth without remeshing. International Journal for Numerical Methods in Engineering, 1999, 46: 131-150. [21] Dolbow,J.E., An Extended Finite Element Method with Discontinuous Enrichment for Applied Mechanics. Ph.D. dissertation, Northwestern University, 1999. [22] Areias,P.M.A. and Belytschko,T., Analysis of three-dimensional crack initiation and propagation using the extended finite element method. International Journal for Numerical Methods in Engineering, 2005, 63: 760-788. [23] Bordas,S.P, Rabczuk,T. and Zi,G., Three-dimensional crack initiation, propagation, branching and junction in non-linear materials by an extended meshfree method without asymptotic enrichment. Engineering Fracture Mechanics, 2008, 75: 943-960. [24] Sun,H.T. and Wang,Y.H., The meshless virtual boundary method and its applications to 2D elasticity problems. Acta Mechanica Solida Sinica, 2007, 20(1): 30-40. [25] Wolf,J.P. and Song,C.M., Finite-element Modelling of Unbounded Media, Wiley, Chichester, 1996. [26] Ooi,E.T. and Yang,Z.J., Modelling multiple cohesive crack propagation using a finite element-scaled boundary finite element coupled method. Engineering Analysis with Boundary Elements, 2009, 33: 915-929. [27] Abaqus 6.7 User Documentation, Dessault Systems, 2007. [28] Matlab R2008a User’s Guide, MathWorks, 2008. [29] Trunk,B., Einfluss der Bauteilgroesse auf die Bruchenergie von Beton. Aedificatio Publishers, 2000 (in German). [30] Feist,C., Numerical Simulations of Localization Effects in Plain Concrete. Ph.D. dissertation, University Innsbruck, 2003. [31] Brokenshire,D.R., A Study of Torsion Fracture Tests. Ph.D. Dissertation, Cardiff University, 1996. [32] Jefferson,A.D., Barr,B., Bennett,T. and Hee,S., Three dimensional finite element simulations of fracture tests using the Craft concrete model. Computers and Concrete, 2004, 1: 261-284. [33] Gasser,T.C. and Holzapfel,G.A., 3D Crack propagation in unreinforced concrete: A two-step algorithm for tracking 3D crack paths. Computer Methods in Applied Mechanics and Engineering, 2006, 195: 5198-5219. [34] Rots,J.G., Computational Modelling of Concrete Fracture. Ph.D. Dissertation, Delft University of Technology, 1988. [35] Du,J., Yon,J.H., Hawkins,N.M., Arakawa,K. and Kobayashi,A.S., Fracture process zone for concrete for dynamic loading. ACI Materials Journal, 1992, 89: 252-258. [36] Belytschko,T., Organ,D. and Gerlach,C., Element-free galerkin methods for dynamic fracture in concrete. Computer Methods in Applied Mechanics and Engineering, 2000, 187: 385-399. [37] Yon,J.H., Hawkins,N.M. and Kobayashi,A.S., Strain-rate sensitivity of concrete mechanical properties. ACI Materials Journal, 1992, 89: 146-153.
ISSN 0894-9166
FINITE ELEMENT MODELLING OF COMPLEX 3D STATIC AND DYNAMIC CRACK PROPAGATION BY EMBEDDING COHESIVE ELEMENTS IN ABAQUS Xiangting Su1,2
Zhenjun Yang2
Guohua Liu1
1
( College of Civil Engineering and Architecture, Zhejiang University, Hangzhou 310027, China) (2 Department of Engineering, the University of Liverpool, L69 3GQ, UK)
Received 23 November 2009; revision received 25 May 2010
ABSTRACT This study proposes an algorithm of embedding cohesive elements in Abaqus and develops the computer code to model 3D complex crack propagation in quasi-brittle materials in a relatively easy and efficient manner. The cohesive elements with softening traction-separation relations and damage initiation and evolution laws are embedded between solid elements in regions of interest in the initial mesh to model potential cracks. The initial mesh can consist of tetrahedrons, wedges, bricks or a mixture of these elements. Neither remeshing nor objective crack propagation criteria are needed. Four examples of concrete specimens, including a wedgesplitting test, a notched beam under torsion, a pull-out test of an anchored cylinder and a notched beam under impact, were modelled and analysed. The simulated crack propagation processes and load-displacement curves agreed well with test results or other numerical simulations for all the examples using initial meshes with reasonable densities. Making use of Abaqus’s rich pre/postprocessing functionalities and powerful standard/explicit solvers, the developed method offers a practical tool for engineering analysts to model complex 3D fracture problems.
KEY WORDS finite element method, cohesive elements, three-dimensional crack propagation, discrete crack model, concrete structures, Abaqus
I. INTRODUCTION Crack propagation is an inherent feature in quasi-brittle materials such as concrete due to the low material tensile strength, and it is a predominant source of nonlinearity and a main culprit for ultimate failure of structures made from these materials. Accurate understanding of the crack propagation behaviour so as to assess the structural load-carrying capacity, by either experimental studies or numerical modelling, remains a worldwide challenge, especially for three-dimensional (3D) problems with complex geometries and loading conditions. Numerical modelling of crack propagation has been an active research field since 1960s[1, 2]. Nowadays there exist a large number of numerical models. In terms of how the cracks are modelled geometrically, there are discrete crack models explicitly separating crack surfaces and modelling discontinuity, smeared crack models based on continuum mechanics, and more indirectly, lattice models. In terms of numerical methods, both the traditional finite element methods (FEM) and boundary element methods (BEM)
Corresponding author. Email: [email protected] Project supported by EPSRC UK (No. EP/F00656X/1). Xiangting Su’s one-year visit to the University of Liverpool was supported by the China Scholarship Council and the National Natural Science Foundation of China (No. 50579081).
· 272 ·
ACTA MECHANICA SOLIDA SINICA
2010
have been widely applied to crack propagation modelling. Extensive literature review of these models and methods in 2D problems has been given elsewhere[3–6] . It appears now that the discrete crack models, mostly based on the cohesive crack model (CCM) developed in ductile materials by Barenblatt[7] and Dugdale[8] and in quasi-brittle materials by Hillerborg[9] , are becoming more and more popular, because of their ability to model macroscopic cracks with strong discontinuity, the CCM’s capability of realistically representing the energy dissipation during fracture processes, and the ease of its implementation as cohesive interface elements (CIE) in the FEM and the BEM. In general, two types of approaches are used to model discrete cohesive crack propagation. The first type of approaches is based on sophisticated remeshing procedures that constantly change the meshes as cracks propagate[3, 10–12] . For problems with crack paths unknown a priori, objective crack propagation criteria are needed to judge when and in which direction a crack propagates. This usually involves calculation of stress intensity factors (SIFs) and stresses at crack tips, whose accuracy can only be ensured by fine crack-tip meshes or using singular elements. This in turn exacerbates the complexity of the remeshing procedure. It is also very difficult for the remeshing-based approaches to deal with complex crack propagation situations, such as multiple cracking and 3D problems. Indeed, the success of the remeshing-based approaches has been largely limited to 2D static fracture problems with single or a few cracks, although they are computationally efficient because a relatively small number of nonlinear CIEs are inserted into the meshes. Another type of approaches pre-insert or pre-embed CIEs between the finite elements in the initial meshes[6, 13–15] . Crack propagation is modelled by automatic opening and merging of the CIEs under applied loadings. Neither remeshing procedures nor crack propagation criteria are needed so complex crack propagation can be modelled. However, these approaches restrict the cracks to the finite element edges or surfaces and the predicted crack patterns may thus be mesh-dependent. Another disadvantage is the high computational cost due to the use of a large number of nonlinear CIEs. The CIEs can also be dynamically inserted when certain crack initiation criterion is satisfied[16, 17] so that the number of CIEs changes with the loading and the total number of CIEs remains relatively low. These procedures of dynamically inserting CIEs need objective crack initiation criteria and change meshes although the remeshing operation is relatively straightforward. In the last two decades, the partition of unity finite element method (PUFEM)[18, 19] and the extended finite element method (XFEM)[20–22] are becoming popular in modelling crack propagation. These methods introduce discontinuity into the finite elements. The crack growth is described independently of the mesh, so that remeshing is not needed. Meshless methods are also receiving much attention[23, 24] . These innovative methods, showing high potential in crack propagation modelling, generally need fine crack tip meshes like the traditional FEM to calculate accurate stresses or SIFs used in evaluating the crack propagation criterion. In this sense, fine initial meshes are needed if the crack paths are unknown in priori, which means high computational cost. The newly-developed semi-analytical scaled boundary finite element[25] is able to compute accurate SIFs without using fine meshes or singular elements. It has recently been applied to model discrete cohesive crack propagation[4, 5, 26], but so far only 2D problems were modelled. In summary, all the above crack models and numerical methods have their own advantages and disadvantages in terms of effectiveness, efficiency and applicability, and the selection largely depends upon the problems analysed and the preference and experience of the analysts. Most of the existing studies have modelled 2D crack propagation. 3D modelling of complex cohesive crack propagation has been sporadically conducted, and there is a long way to go before it becomes a routine exercise in research and practical design due to the inherent difficulties. In addition, most of the studies have used specialised in-house computer programs (e.g., Cornell University’s FRANC) that may be difficult for practical engineers to use. Most of the general-purpose commercial FEA packages, such as Ansys and Abaqus, are unable to model complicated crack propagation without user intervention or second development. This study aims at developing a practical finite element method for modelling complex 3D crack propagation in quasi-brittle materials under static and dynamic loadings. The method makes use of a special type of CIEs called ‘cohesive elements’, designed for modelling bonded interfaces in Abaqus V6.5 or higher[27] , to model potential cracks. An efficient algorithm was devised to insert these CIEs into concerned regions in FE meshes consisting of various types of solid elements. The method is simple to implement. The rich pre/post-processing functionalities and powerful implicit and explicit solvers
Vol. 23, No. 3
Xiangting Su et al.: Complex 3D Static and Dynamic Crack Propagation
· 273 ·
of Abaqus can be fully exploited. Its effectiveness is validated by modelling various fracture problems of concrete structures.
II. THE MODELLING METHOD The proposed method involves the following steps: (1) Meshing the domain using Abaqus/CAE and generating an input file. (2) Inserting cohesive elements into the regions of interest in the initial FE mesh using an in-house computer program and generating another input file. The initial mesh may consist of solid elements such as 4-noded tetrahedrons (C3D4), 6-noded wedges (C3D6) and 8-noded bricks (C3D8), or a mixture of these elements. Both 6-noded and 8-noded cohesive elements (COH3D6 and COH3D8) can be inserted. Figure 1 shows a mesh with solid elements only, the mesh after CIEs are inserted and the 3D CIEs, respectively. (3) Solving the problem using Abaqus standard or explicit solvers. This procedure can be readily automated by running a batch file. The following briefly discusses the key elements in the method.
Fig. 1. FE meshes before and after cohesive elements are inserted and 3D cohesive elements.
2.1. Cohesive Elements with Damage in Abaqus The cohesive crack model[7–9] assumes the existence of a fracture process zone (FPZ) in front of the real crack tip, in which energy dissipation occurs during fracture. In the FPZ, there exist tractions in the normal direction (tn ) and the two tangential or shear directions (ts and tt ) across the crack surfaces, resulting from mechanisms such as material bonding, aggregate interlocking and surface friction. Figure 2 shows a curve relating the normal traction tn and the crack opening displacement δn as an example. Similar traction-separation curves can be defined for shear tractions (ts and tt ) - crack sliding displacements Fig. 2 An example of traction-separation curve with softening (δs and δt ) relations. Before the crack initiates, a for cohesive elements. linear elastic ascending phase is assumed to model the initially un-cracked material. After the crack initiates, the traction decreases monotonically as functions of the corresponding separation, which is often termed tension or strain softening. The initial tensile stiffness kn0 and the initial shear stiffness ks0 and kt0 should be high enough to represent the un-cracked material, but not too high to cause numerical ill-conditioning. These initial stiffness values are determined by a trial and error approach. If δn is negative during loading increments or iterations, a compressive stiffness of magnitude equal to kn0 is assigned in order to prevent penetration of crack surfaces.
· 274 ·
ACTA MECHANICA SOLIDA SINICA
2010
The cohesive elements in Abaqus[27] are based on the cohesive crack model. The constitutive response of cohesive elements, defined in terms of traction-separation laws, assumes initially linear elastic behaviour followed by the initiation and evolution of damage. Several damage initiation criteria are available in Abaqus. This study used the quadratic nominal stress criterion in all the examples. The damage evolution criteria can be classified into three categories according to the relation between traction and separation: linear damage evolution, exponential damage evolution, and tabular damage evolution. All these criteria were used in this study. 2.2. Inserting Cohesive Elements into 3D FE Meshes Although inserting cohesive elements into 2D FE meshes is relatively straightforward[6], it is not a trivial task for 3D problems, especially if the FE mesh is not regular and consists of different types of solid elements. One challenge is how to robustly deal with the changes in the complicated 3D nodal and elemental connectivity due to the insertion of CIEs. This is tackled by an effective embedding algorithm implemented in a MATLAB[28] program with carefully designed data structures listed in Table 1. It has the following steps: Table 1. Data structures used in the algorithm of inserting cohesive elements in 3D meshes
Structure Name
Members
NODES index elemconn faceconn xcord ycord zcord isInDefinedArea
ELEMS index faceconn nodeconn type
FACES index nodeconn elemconn isBound isInDefinedArea
BALLS POINT isInDefinedArea
POINT elemindex xcord ycord zcord newnodeindex
(1) Reading the nodal coordinates and the nodal connectivity of the solid elements from the Abaqus input file. These data are then used to generate three arrays of type Structure: NODES, ELEMS and FACES. Most of the members in these structures are self-explanatory. For example, the members elemconn and faceconn in NODES structure are arrays containing the elements and faces the node is connected with, respectively. It is worth noting that only those nodes, elements and surfaces in the concerned or defined regions in the mesh are considered, in order to reduce the computational cost. This is indicated by the member isIndefinedArea. A consistent nodal ordering approach (clockwise or anti-clockwise) in nodeconn of FACES should be used to facilitate determination of correct nodal order when cohesive elements are created in the 4th step. (2) Inserting new nodes. For each existing node with isIndefinedArea=1, all the N number of solid elements connected with this node are first identified. N number of new nodes, each connected with one solid element, are then generated with the same coordinates of the existing node. All the new nodes are stored in an array of POINT, which is a member of a Structure BALLS associated with the existing node. The member elemindex in POINT represents the solid element the new node is connected with. An array of BALLS is then generated for all the existing nodes. The design of BALLS greatly facilitates subsequent generation of cohesive elements and is at the heart of the embedding algorithm. (3) Updating the arrays NODES and ELEMS. Since a large number of new nodes are generated, the elemental connectivity (elemconn) of NODES and the nodal connectivity (nodeconn) of ELMS are updated. It should be noted that the array FACES is not updated at this step. (4) Generating cohesive elements. This is based on the array FACES. There are two solid elements (elemconn) connected with each face. For each vertex of a face (nodeconn), a pair of new nodes connected to the two solid elements is found in POINT of BALLS. The three or four pairs of nodes found for this face comprise a cohesive element (COH3D6 or COH3D8). As the nodal order in FACES is carefully designed in the 1st step, the requirement on the nodal order of cohesive elements is readily satisfied. (5) Creating the Abaqus input file. The file contains the index and coordinates of POINT, the index and nodal connectivities of the updated ELEMS and the cohesive elements. To minimise the mesh-dependence for problems with crack paths unknown a priori, a reasonably fine initial mesh should be used. This means generally a large number of solid elements, nodes, surfaces
Vol. 23, No. 3
Xiangting Su et al.: Complex 3D Static and Dynamic Crack Propagation
· 275 ·
and CIEs are used. Computer implementation of the above inserting algorithm must be carried out with care, or the time spent on inserting CIEs may be unduly long. As an example, for an initial mesh with 46800 C3D8 elements and 52140 nodes, it takes about 30 minutes for the developed MATLAB program to insert 135,200 CIEs using a PC with an Intel Xeon [email protected] GHz and 3 GB physical memory.
III. NUMERICAL EXAMPLES Four 3D concrete structures subjected to static and dynamic loadings were modelled using the developed method as examples. A PC with two Intel Xeon CPUs @3.16GHz and 3 GB physical memory was used. Both Abaqus/standard and explicit solvers were used to solve the nonlinear equation systems. When the quasi-static explicit solver was used to model static problems, the analysis time was carefully chosen to avoid dynamic effects. Because the failure mode of all the examples is mainly tensile fracture and there is no experimental data about shear fracture resistance, the strength and the initial stiffness in shear are assumed as 10 times those in tension for the cohesive elements so that shear fracture failure does not occur. 3.1. Wedge-Splitting Test The first example is the wedge-splitting test of a concrete specimen CP250 carried out by Trunk[29] . The test was modelled in 2D by Trunk[29] and Feist[30] , and in 3D by Areias and Belytschko[22]. The geometry and the loading condition are shown in Fig.3. The concrete has a Young’s modulus E = 28300 MPa, the Poisson ratio ν = 0.2 and mass density ρ = 2.5 × 10−6 kg/mm3 . A bilinear softening curve shown in Fig.4 with parameters obtained from the experiment was used to model the normal tractionseparation relation of the CIEs. The area under the curve is the fracture energy Gf = 0.49 N/mm. The initial stiffness kn0 was assumed as 1 × 105 MPa·mm. The Abaqus/explicit quasi-static solver used an analysis time of 0.1 s.
Fig. 4. Bilinear softening law used for the wedgesplitting test. Fig. 3. Wedge-splitting test specimen (unit: mm).
Figure 5 shows a deformed FE mesh with 774 C3D6 and 1392 C3D8R solid elements, 1854 CIEs and 6852 nodes in total. It can be seen that the crack path was vertical as expected because of the symmetry of the mesh and loading conditions. Finer meshes were also modelled but they led to little difference from the predicted crack path and the load-crack mouth opening displacement (CMOD) curves. Figure 6 compares the load-CMOD curves from the experiment and modelling using different solvers and solid elements. The standard solver always failed by divergence at the post-peak stage, e.g., it failed at CMOD=3.2 mm for C3D8 as shown in Fig.6. This also often occurred in 2D simulations[6] . In contrast, the explicit solver always managed to model the whole fracture process when small enough time increments were used. The explicit solver using C3D8R with reduced integration and hourglass control (HC) resulted in nearly the same pre-peak responses as the standard solver using C3D8, while the explicit solver using C3D8R but without hourglass control led to too soft pre-peak responses. C3D8 was also tried with the explicit solver but the computational time nearly doubled when C3D8R was
· 276 ·
ACTA MECHANICA SOLIDA SINICA
2010
used. Therefore, the explicit solver was used in all the following simulations and C3D8R elements with hourglass control were chosen when 8-noded brick elements were used.
Fig. 5. Deformed mesh (CMOD=4 mm, scale=100).
Fig. 6. Load-CMOD curves: effect of solid element types and solvers.
Figure 7 compares the load-CMOD curves from 3D and 2D simulations. In the 2D simulation, a similar mesh to the in-plane mesh of the 3D mesh (Fig.5) was used assuming a plane stress condition. It can be seen that the peak load predicted from the 3D simulation is higher than that from the 2D, which reflects the constraining effect of the specimen thickness. A very similar mesh to Fig.5 but with a non-symmetric distribution of elements along the vertical central line was also modelled. The deformed mesh is shown in Fig.8. The crack path deviated from the central line, Fig. 7 Load-CMOD curves from 2D and 3D. indicating the effect of the initial mesh. The loadCMOD curves predicted from the two meshes are shown in Fig.9.
Fig. 8. Deformed mesh from a non-symmetric mesh (CMOD=4 mm, scale=100).
Fig. 9. Load-CMOD curves: effect of mesh arrangement.
3.2. Torsion Fracture Test The second example is a concrete bar under torsion tested by Brokenshire[31]. The geometry and the loading condition are illustrated in Fig.10. 3D modelling was carried out by Jefferson et al.[32] and Gasser and Holzapfel[33] . The three support points were assumed as fixed and the loading incremental procedure was controlled by the vertical displacement at the point C. The concrete has E = 35000 MPa, ν = 0.2 and
Vol. 23, No. 3
Xiangting Su et al.: Complex 3D Static and Dynamic Crack Propagation
· 277 ·
Fig. 10. A concrete bar under torsion (unit: mm).
ρ = 2.5 × 10−6 kg/mm3 . The same linear softening curve (Fig.11) as used in Ref.[32] was adopted to model the normal traction-separation relation with Gf = 0.08 N/mm. The initial stiffness was set as kn0 = 7.5 × 106 MPa·mm. An analysis time of 0.3 s was used in the explicit quasi-static procedure. CMOD is defined as the relative displacement between the points A and B normal to the Fig. 11 Linear softening law used for the torsion test. notch surface (ref. Fig.10). Two FE meshes with different solid element types were modelled to investigate the effect on the prediction of crack surfaces. Mesh A in Fig.12(a) has 10651 C3D4 solid elements and 9606 CIEs, and Mesh B in Fig.12(b) has 7027 C3D8R and C3D6 elements and 8055 CIEs. The predicted crack surfaces are compared with the test result in Fig.13. It is clear that Mesh A with C3D4 led to very realistic, tortuous fracture surfaces whereas Mesh B with C3D8R predicted an unrealistically flat and smooth crack path. The predicted load-CMOD curves are presented in Fig.14. It can be seen that both meshes led to curves close to the experimental data. Mesh A predicted a peak load slightly higher than Mesh B. This may be explained by the predicted crack surfaces in Fig.13, where the tortuous and unsmooth crack surfaces from Mesh A offer higher fracture resistance than the flat and smooth ones from Mesh B. Jefferson et al.[32] also modelled the same test using a similar mesh to Mesh B and the ‘Craft’ concrete
Fig. 12. Two FE meshes used for the torsion test.
· 278 ·
ACTA MECHANICA SOLIDA SINICA
2010
Fig. 13. Crack surfaces after failure for the torsion test.
Fig. 14. Load-CMOD curves for the torsion test.
Fig. 15. Pull-out test of a steel bar embedded in concrete (unit: mm).
model. The predicted load-CMOD curve is also shown in Fig.14 and it is very close to that from Mesh B in this study. 3.3. Pull-out Test of a Steel-Embedded Concrete Anchor This example models the pull-out behaviour of a steel anchor embedded in concrete. The detailed geometry, applied boundary and loading conditions of one-quarter specimen are illustrated in Fig.15. This example was also modelled in Ref.[19,22,23] using other methods. The anchor bolt, excluded from the FE model, was simulated by a boundary condition of vertical displacement imposed on the top
Fig. 16. Meshes used for the pull-out test.
Vol. 23, No. 3
Xiangting Su et al.: Complex 3D Static and Dynamic Crack Propagation
· 279 ·
surface of the groove. The concrete properties are: E = 30000 MPa, ν = 0.2 and ρ = 2.5×10−6 kg/mm3 . A linear softening curve with Gf = 0.1 N/mm, tensile strength=3 MPa and kn0 = 6 × 104 MPa·mm was used. An analysis time of 0.1 s was used in the explicit solver. To investigate the mesh-dependence of results, two meshes were modelled. The coarse mesh in Fig.16(a) has 31889 C3D4 solid elements and 60164 CIEs, and the numbers for the fine mesh in Fig.16(b) are 69798 and 138169, respectively. Figure 17 shows the load (F )-displacement (d) curves from the two meshes and other studies. It can be seen that the two meshes predicted very close responses, especially the same failure point when the load suddenly dropped at d = 0.25 mm.
Fig. 17. Load-displacement curves for the pull-out test.
Similar crack propagation processes and final crack patterns were predicted from the two meshes. Figure 18 illustrates four stages of the predicted crack propagation using the fine mesh, where the cracks are represented by grey areas consisting of highly-damaged cohesive elements (with the damage
Fig. 18. Predicted crack propagation process of the pull-out test.
· 280 ·
ACTA MECHANICA SOLIDA SINICA
2010
index D ≥ 0.99). At F = 310 kN (60% of the peak load Fu = 503 kN), tensile cracks initiated at the top edge of the groove, forming a cone around the groove (Fig.18(a)). As the load approached the peak, a splitting-mode longitudinal crack initiated at the top surface and propagated in the radial direction (Fig.18(b)). Beyond the peak load (Fig.18(c)), the longitudinal crack propagated rapidly, leading to a drastic decline on the load-displacement curve (Fig.17). After the complete propagation of the longitudinal crack (Fig.18(d)), the conical crack continued to advance with the load decreasing slowly. The simultaneous propagation of the longitudinal crack and the conical crack, observed in the experiments by Rots[34] , was well reproduced. This may explain the difference between the predicted load-displacement curves (Fig.17) by the present study and other studies[22, 23] which only modelled the conical crack. 3.4. Impact Test of a Concrete Beam The fourth example is a concrete beam under impact tested and modelled in 2D by Du et al.[35] . 3D modelling was reported by Belytschko et al.[36]. The geometry and the boundary condition are shown in Fig.19. The impact load history obtained from the test (Fig.20) was input as the loading condition. The concrete properties are: E = 34480 Fig. 19 Impact test of a concrete beam (unit: mm). MPa, ν = 0.2 and ρ = 2.5 × 10−6 kg/mm3 . An exponential softening curve (Fig.21) with the exponential law parameter α = 1.5 and the equivalent dynamic fracture energy Gf = 0.152 N/mm, suggested in Ref.[35], was adopted to model the CIEs. The initial stiffness kn0 = 1 × 105 MPa·mm was used. Figure 22 shows a deformed mesh at the post-peak stage when t = 1 ms. The mesh consists of 3645 solid elements and 606 cohesive elements. Figure 23 compares favourably the predicted loaddeflection curve with the test data. Figure 24 compares the measured and computed strain histories at three locations (SG01, SG02 and SG03). The strain history computed for SG03 agreed well with the measured. For the other two locations, the strain histories were underestimated, probably due to
Fig. 21. Exponential softening law for the impact test. Fig. 20. Load history used for the impact test.
Fig. 22. Deformed mesh for the impact test.
Fig. 23. Load-displacement curves for the impact test.
Vol. 23, No. 3
Xiangting Su et al.: Complex 3D Static and Dynamic Crack Propagation
· 281 ·
the negligence of shrinkage effects in the numerical modelling[37] . The crack-tip extension histories are compared in Fig.25.
Fig. 24. Strain history for the impact test.
Fig. 25. Crack-tip extension history for the impact test.
IV. CONCLUSIONS A simple yet effective finite element methodology has been developed to simulate 3D complex cohesive crack propagation in quasi-brittle materials, using Abaqus’s cohesive elements. An efficient algorithm to insert the cohesive elements into initial finite element meshes was devised and implemented by an in-house computer program. Several concrete examples were modelled to demonstrate the effectiveness of the method. It is concluded that, making use of the rich pre/post-processing functionalities and powerful standard/explicit solvers in Abaqus, the developed methodology seems promising to provide a practical simulation tool for modelling realistic 3D fracture in various industries.
References [1] Scordelis,A.C. and Ngo,D., Finite element analysis of reinforced concrete beams. Journal of the American Concrete Institute , 1967, 64: 152-163. [2] Rashid,Y.R., Ultimate strength analysis of prestressed concrete pressure vessels. Nuclear Engineering and Design, 1968, 7(4): 334-344. [3] Yang,Z.J. and Chen,J.F., Fully automatic modelling of cohesive discrete crack propagation in concrete beams using local arc-length methods. International Journal of Solids and Structures, 2004, 41(3-4): 801826. [4] Yang,Z.J., Fully automatic modelling of mixed-mode crack propagation using scaled boundary finite element method. Engineering Fracture Mechanics, 2006, 73(12): 1711-1731. [5] Yang,Z.J. and Deeks,A.J., Fully-automatic modelling of cohesive crack growth using a finite element-scaled boundary finite element coupled method. Engineering Fracture Mechanics, 2007, 74(16): 2547-2573. [6] Yang,Z.J., Su,X.T., Chen,J.F. and Liu,G.H., Monte Carlo simulation of complex cohesive fracture in random heterogeneous quasi-brittle materials. International Journal of Solids and Structures, 2009, 46(17): 32223234. [7] Barenblatt,G.I., The formation of equilibrium cracks during brittle fracture: general ideas and hypothesis, axially symmetric cracks. Applied Mathematics and Mechanics, 1959, 23: 622-636. [8] Dugdale,D.S., Yielding of steel sheets containing slits. Journal of Mechanics of Physics and Solids, 1960, 8: 100-104. [9] Hillerborg,A., Modeer,M. and Petersson,P., Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cement and Concrete Research, 1976, 6: 773-782. [10] Wawrzynek,P.A. and Ingraffea,A.R., An interactive approach to local remeshing around a propagating crack. Finite Elements in Analysis and Design, 1989, 5: 87-96. [11] Xie,M. and Gerstle,W.H., Energy-Based Cohesive Crack Propagation Modeling. Journal of Engineering Mechanics-ASCE, 1995, 121(12): 1349-1358. [12] Yang,Z.J. and Xu,X.F., A heterogeneous cohesive model for quasi-brittle materials considering spatially varying random fracture properties. Computer Methods in Applied Mechanics and Engineering, 2008, 197: 4027-4039.
· 282 ·
ACTA MECHANICA SOLIDA SINICA
2010
[13] Xu,X.P. and Needleman,A., Numerical simulations of fast crack growth in brittle solids. Journal of the Mechanics and Physics of Solids, 1994, 42(9): 1397-1434. [14] L´ opez,C., Carol,I. and Aguado,A., Meso-structural study of concrete fracture using interface elements. I: numerical model and tensile behavior. Materials and Structures, 2008, 41(3): 583-599. [15] L´ opez,C., Carol,I. and Aguado,A., Meso-structural study of concrete fracture using interface elements. II: compression, biaxial and Brazilian test. Materials and Structures, 2008, 41(3): 601-620. [16] Camacho,G.T. and Ortiz,M., Computational modelling of impact damage in brittle materials. International Journal of Solids and Structures, 1996, 33: 2899-2938. [17] Pandolfi,A. and Ortiz,M., An efficient adaptive procedure for three-dimensional fragmentation simulations. Engineering with Computers, 2002, 18: 148-159. [18] Melenk,J.M. and Babuska,I., The partition of unity finite element method: Basic theory and applications. Computer Methods in Applied Mechanics and Engineering, 1996, 139: 289-314. [19] Gasser,T.C. and Holzapfel,G.A., Modeling 3D crack propagation in unreinforced concrete using PUFEM. Computer Methods in Applied Mechanics and Engineering, 2005, 194: 2859-2896. [20] Moes,N., Dolbow,J. and Belytschko,T., A finite element method for crack growth without remeshing. International Journal for Numerical Methods in Engineering, 1999, 46: 131-150. [21] Dolbow,J.E., An Extended Finite Element Method with Discontinuous Enrichment for Applied Mechanics. Ph.D. dissertation, Northwestern University, 1999. [22] Areias,P.M.A. and Belytschko,T., Analysis of three-dimensional crack initiation and propagation using the extended finite element method. International Journal for Numerical Methods in Engineering, 2005, 63: 760-788. [23] Bordas,S.P, Rabczuk,T. and Zi,G., Three-dimensional crack initiation, propagation, branching and junction in non-linear materials by an extended meshfree method without asymptotic enrichment. Engineering Fracture Mechanics, 2008, 75: 943-960. [24] Sun,H.T. and Wang,Y.H., The meshless virtual boundary method and its applications to 2D elasticity problems. Acta Mechanica Solida Sinica, 2007, 20(1): 30-40. [25] Wolf,J.P. and Song,C.M., Finite-element Modelling of Unbounded Media, Wiley, Chichester, 1996. [26] Ooi,E.T. and Yang,Z.J., Modelling multiple cohesive crack propagation using a finite element-scaled boundary finite element coupled method. Engineering Analysis with Boundary Elements, 2009, 33: 915-929. [27] Abaqus 6.7 User Documentation, Dessault Systems, 2007. [28] Matlab R2008a User’s Guide, MathWorks, 2008. [29] Trunk,B., Einfluss der Bauteilgroesse auf die Bruchenergie von Beton. Aedificatio Publishers, 2000 (in German). [30] Feist,C., Numerical Simulations of Localization Effects in Plain Concrete. Ph.D. dissertation, University Innsbruck, 2003. [31] Brokenshire,D.R., A Study of Torsion Fracture Tests. Ph.D. Dissertation, Cardiff University, 1996. [32] Jefferson,A.D., Barr,B., Bennett,T. and Hee,S., Three dimensional finite element simulations of fracture tests using the Craft concrete model. Computers and Concrete, 2004, 1: 261-284. [33] Gasser,T.C. and Holzapfel,G.A., 3D Crack propagation in unreinforced concrete: A two-step algorithm for tracking 3D crack paths. Computer Methods in Applied Mechanics and Engineering, 2006, 195: 5198-5219. [34] Rots,J.G., Computational Modelling of Concrete Fracture. Ph.D. Dissertation, Delft University of Technology, 1988. [35] Du,J., Yon,J.H., Hawkins,N.M., Arakawa,K. and Kobayashi,A.S., Fracture process zone for concrete for dynamic loading. ACI Materials Journal, 1992, 89: 252-258. [36] Belytschko,T., Organ,D. and Gerlach,C., Element-free galerkin methods for dynamic fracture in concrete. Computer Methods in Applied Mechanics and Engineering, 2000, 187: 385-399. [37] Yon,J.H., Hawkins,N.M. and Kobayashi,A.S., Strain-rate sensitivity of concrete mechanical properties. ACI Materials Journal, 1992, 89: 146-153.