Abstract
The mixed freedom finite element method proposed for contact problems was extended to simulate the fracture mechanics of concrete using the fictitious crack model. Pairs of contact points were set along the potential developing path of the crack. The displacement of structure was chosen as the basic variable, and the nodal contact force in contact region under local coordinate system was selected as the iteration variable to confine the nonlinear iteration process in the potential contact surface which is more numerically efficient. The contact forces and the opening of the crack were obtained explicitly enabling the softening constitutive relation for the concrete to be introduced conveniently by the fictitious crack model. According to the states of the load and the crack, the constitutive relation of concrete under cyclic load is characterized by six contact states with each contact state denoting its own displacement-stress relation. In this paper, the basic idea of the mixed freedom finite element method as well as the constitutive relation of concrete under cyclic load is presented. A numerical method was proposed to simulate crack propagation process in concrete. The accuracy and capability of the proposed method were verified by a numerical example against experiment data.
1. Introduction
Concrete is a highly heterogeneous material with its properties varyong widely from point to point due to the presence of high strength aggregates, medium strength mortar, and weaker mortar-aggregate interfaces. The existence of geometric voids which act as stress raisers can also contribute to the changing property [1]. Due to the low tensile strength of concrete material, cracking is the dominant failure mechanism in concrete structures, and the mechanism of crack initiation and growth has received significant research effort.
In recent years, considerable amount of effort has been devoted to develop numerical models to simulate the fracture behavior of concrete used in civil engineering structures [2]. A number of modeling approaches have been implemented, such as plasticity models [3, 4], microplane models [5, 6], and fracture mechanics models. Generally, in the field of fracture mechanics, the numerical methods based on the finite element method (FEM) are classified into two groups [7]: “smeared crack approach” and “discrete crack approach.” In the smeared crack approach, the fracture is represented in a smeared manner in which an infinite number of parallel cracks with infinitely small opening are (theoretically) distributed (smeared) over the finite element [8]. Those cracks are usually modeled on a fixed finite element mesh and their propagation is simulated by the reduction of the stiffness and strength of the material. Most of the models in this approach treat concrete as a quasi-brittle material that exhibits strain-softening behavior under tensile loading [9]. The constitutive laws, defined by stress-strain relations, are nonlinear and involve a strain softening behavior which can be challenging in the analysis. The system of equations may become ill-posed [10] and can introduce localization instabilities and spurious mesh sensitivity to the finite element calculations [8]. In the field of fracture mechanics, there are a number of models that can deal with the nonlinear zone ahead of a crack tip in concrete, including the fictitious crack model (FCM) [11], the crack band model (CBM) [12], the two-parameter fracture model (TPFM) [13, 14], the size effect model (SZM) [15], the effective crack model (ECM) [16, 17], and the double-K fracture model (DKFM) [18–20].
The discrete approach is preferred when there is one crack, or a finite number of cracks, in the structure [2]. Among them, the cohesive zone model (CZM), which assumes that stresses act across a narrowly open crack, was developed by Hillerborg et al. [11] under the name FCM. In the literature, it is generally accepted that the FCM is the best simple fracture mechanics model as it provides reasonable approximation of the crack propagation in concrete using only simple calculations. The model is mainly applicable to mode-I fracture, but recently it has been extended to model shear failure mode-II and mode-III cracks [21]. The FCM has been used to investigate the fracture process in concrete under monotonic and cyclic loading. Various constitutive relationships for cyclic tension behaviour of concrete, such as the focal point model [22] and the continuous function model [23], have been developed. More recently, linear and bilinear constitutive relationships describing the unloading and reloading scenarios were developed [24, 25]. Although these FE-based models claim to accurately simulate the fracture behaviour of concrete structures, they are usually computationally expensive. This drawback makes them unsuitable for damage detection applications where adaptable and computationally efficient models are needed [9].
This paper presents an efficient approach of modeling the monotonic and cyclic flexural cracking behaviour of concrete using a FCM. Under the cyclic load, the constitutive relation of concrete fracture is very complex, the same displacement may correspond to three different kinds of stress states, and numerical conditions of the transformation between each stress state are not easy to define. In numerical analysis, it will lead to instability problems during the overlapping process of positive and negative stiffness. In the analysis of concrete crack propagation under cyclic load, the states of load (loading, unloading, and reloading) at joint interface should be defined first because different states correspond to different displacement-stress curves. And the softening load that fictitious crack plane transfer is a function of opening displacement on this plane based on fictitious crack model, but the opening displacement is an unknown quantity still to be determined, which means that the propagation calculation of concrete cracks under cyclic load is a nonlinear iterative process. This situation is very similar to contact problems in practical engineering in which the state of load process, contact forces, and displacements on interfaces cannot be known in advance. The solution is to make the equation fulfill the equilibrium condition in general and satisfy the constitutive relation between opening displacement and interface tension on interface from the determined states of load and crack.
The paper is structured as follows. In the next section, the basic idea of the mixed freedom finite element method for contact problem is given. Section 3 describes the main features of the proposed method, including the improved constitutive relationships for the crack propagation simulation under cyclic load. The numerical implementation and solution algorithms are given in Section 4. Section 5 reports the validation of the model and discusses the findings. Finally, concluding remarks are summarized in Section 6.
2. Basic Idea of the Mixed Freedom Finite Element Method
Pairs of contact points were set along the potential developing path of the crack. The mixed freedom FEM proposed for contact problem was extended to simulate the fracture mechanics of concrete using the FCM. The system of forces acting on the structure was divided into two parts: external forces and contact forces. The displacement of structure was chosen as the basic variable, and the nodal contact force in potential contact region under local coordinate system was chosen as the iteration variable, so that the nonlinear iteration process was only limited in the contact surface and much more efficient. The contact forces and the opening of the crack can be obtained explicitly enabling the softening constitutive relation for the concrete to be introduced conveniently by the fictitious crack model.
In Figure 1, the body is divided into two parts by the potential developing path of crack which was set in advance. Now we focus our attention on an arbitrary node pair constituted by points 1 and 2. is the displacement of the node pair, and is the contact force. Here, the superscript denotes the two body parts. is the local coordinate system on the interface with being the normal direction and , the tangent directions. Then, the crack opening displacement (COD) in normal direction between point 1 and 2 can be computed as .
Assuming that the analysis of th step has been finished, the incremental static equilibrium equation for th step is stated as follows by performing a finite element discretization of each body and : where is the global stiffness matrix; is the vector of displacement increment at step ; is the vector of total external load at step ; is the vector of total contact force at step ; is the strain-displacement matrix; is the Cauchy stress tensor; is the vector of incremental contact force. Clearly, the first term of the right-hand side (RHS) of the above equation, that is, the content in the bracket, stands for the external load increment of current load step since the contribution of to the system is already included into .
If a typical decomposition of global stiffness matrix was performed at the beginning of the analysis for one and only one time, (1) could be rewritten as
If is defined as and the matrix is introduced into the above equation, (3) is rewritten as where the matrix is the flexibility matrix which is defined on the possible contact boundary . An arbitrary component of , , represents the flexibility coefficient corresponding to the displacement at the freedom due to a unit force at the freedom . Here, or is only limited in the freedom of the region where contact is likely to take place.
It is important to be pointed out that (3) and (4) and the following equations in this section are only performed for the degree of freedoms (DOFs) on the potential contact surface, not for all the DOFs of the whole system.
Applying (4) to points 1 and 2 of a given node pair, respectively, gives
According to Newton’s third law, it is obvious that , and, moreover, incorporating flexibility matrix and into , so that if we subtract (5) from (6), we obtain
Equation (7) is the finite element compliance equation of the mixed finite element method for the static contact problems with friction and initial gaps proposed in this paper. In this equation, the second term in RHS, , stands for the difference of incremental displacement only due to the external load increment, as can be seen easily from (3); therefore, it has nothing to do with the current contact state and can be obtained by back substitution directly. However for the first term in RHS, , it denotes the difference of incremental displacement induced by both the external load increment and the contact force increment. From this point of view, we can see that both the RHS and the left-hand side (LHS) of the above equation are associated with the contact state. An iterative method taking into account different contact states is necessary to solve (7), and this will be given in detail in the next section.
3. FCM for the Crack Propagation under Cyclic Load
The strain softening behaviour can be applied to establish the so-called FCM according to Hillerborg [26] to simulate crack formation. In Figure 2, a zone with more microcracks called fracture process zone will appear before crack propagation because of the inhomogeneity of concrete. The fictitious cracks in this zone can transfer tensile stress, and the stress value can be obtained from the softening curve given its crack width. The nonlinear behavior of the propagating crack is described by a fictitious crack with cohesive forces acting between its interfaces [27]. At the crack tip, the maximum stress is equal to the tensile strength , and, for all other points along the crack path, the transmitted stress depends on the crack opening and is defined by the strain softening diagram. The area under the curve is equal to the specific fracture energy . The fracture energy is the ratio of the energy necessary for total fracture of a specimen to the projection of the fracture surface area. is related to material characteristic constants such as concrete mixture ratio, strength, aggregate type, particle diameters, and cement grade and will be determined by experiment. According to the widely adopted CEB-FIP Model Code 1990, can be determined by where MPa and is the maximum aggregate size. The value of is of certain relation with the aggregate size; for example, in CEB-FIP Model Code 1990, it is recommended that , is the parameter relative to the maximum aggregate size.
There are different kinds of loading/unloading model proposed by many researchers for crack propagation under cyclic load, such as the simple elastic loading and unloading model, Hordijk-Reinhardt model [28, 29], and the Toumi model [24]; see Figures 3(a)~3(c).
(a)
(b)
(c)
Based on the Toumi model [24], an improved constitutive model shown in Figure 4 is proposed for the crack propagation under cyclic load. The main improvement lies in that the loading and reloading curve equation of Toumi model is adjusted to ensure the continuity of displacement and stress during the crack propagation which is convenient for FEM implementation. The equations used in this improved curve are shown as follows.
(1) The softening curve of concrete is based on Cornelissen’s (softening) function, which can be expressed as [30] where coefficients , takes the recommended values 1.0 and 5.64, respectively. The ultimate crack developing width was determined from the assumption that the area under the softening curve is equal to the specific fracture energy and can be approximated by .
(2) The equation of the unloading curve at point shown in Figure 4 is where , represent the stress at point and displacement of interface opening during unloading and is a material characteristic parameter describing the relationship between the compressive stress required to reclose the crack with unloading stress and the tensile strength in fracture process zone. According to the experimental results of Reinhardt et al. [31], a recommended value of is used in this paper. In cracked areas, when unloading if , then . Otherwise, if , we have , and the stress can be calculated using (10).
(3) The equation of reloading curve at point in the proposed model is where is the shape parameter of reloading curve, which can be approximated by according to the experimental data of [31].
The stress value at point is governed by the reduction factor together with the stress value at point and can be determined using the following equation:
The value of takes as a constant of 0.05 according to Horri’s suggestion [29].
4. Numerical Implementation
In the propagation of the crack in concrete, there are three loading process states, that is, loading, unloading, and reloading. For model I crack, there are three interface gap states, that is, closure, opening, and virtual opening. In order to describe the whole loading/unloading process of model I crack in a physical model, six contact states are employed. In Figure 4, the corresponding positions of the six contact states in the constitutive relation of concrete are shown, respectively: ① the initial closed state in noncrack zone; ② fictitious crack extension state in fracture process zone while loading; ③ fictitious crack unloading state in fracture process zone and macro crack zone while unloading; ④ closed crack state in fracture process zone and macrocrack zone while being fully reclosed; ⑤ fictitious crack reloading state in fracture process zone while reloading; ⑥ the opening state in macro crack zone. For the crack of models II and III, we can add an extra state which is the sliding state of the fracture process zone.
When using mixed freedom FEM to simulate fracture process, the initial condition was taken from the converged result (contact state and contact force between crack surfaces) of former FEM step. Using the overall displacement field obtained from (1), the joint opening displacement can be determined and employed to solve the new contact force on joint interface iteratively. Assuming that the th contact iteration has been finished, the displacement and stress expression for th iteration is stated as follows: Here, the superscript denotes the contact iteration step while the subscript stands for the loading increment step. The displacement and stress from the above equation are prediction results and must be checked against the contact state after th iteration. Based on the calculated result from the former loading step and th contact iteration, as well as the joint opening displacement and contact stress after th iteration, the contact states were processed and updated which will be explained in detail as follows.
① If a contact state is at the initial closed state, where the initial condition is , , , . After th iteration, if the normal stress of the crack surfaces satisfied , the contact state is updated to fictitious crack extension state, and is set for the next iteration; otherwise, if was satisfied, the current state is changed to sliding state, and a new condition of is set accordingly. is the normal contact force component; , are two tangential components of the contact forces; is the friction coefficient, is the controlled area of pair nodes; is cohesive strength; is the tensile strength of the crack surface; is the residual tensile strength of the crack surface.
② If the current state is at the fictitious crack extension state, where the condition is , the normal stress must satisfy (9). If the opening width reaches the condition of , the current state is adjusted to fictitious crack unloading state, and contact force is set as , stress at the unloading point is set as , and the residual tensile strength of the contact pair is changed to . If is met, the current state is changed to opening state, and new values are set to make , . is the normal stress of crack at th loading step.
③ If the current state is at the fictitious crack unloading state, where the condition is , the normal stress must satisfy (10). If the opening width of crack after iteration meet , the current state is changed to fictitious crack reloading state, and contact force is set as , and stress at the reloading point is set as , . If is satisfied, the current state is changed to crack closed state, and the width is set to to avoid crack surface embedding.
④ If current state is at the crack closed state, where the condition is , , and . If the normal stress after iteration satisfies , the contact state is changed to fictitious crack reloading state and is set as , where is the stress when the crack is fully closed.
⑤ If the current state is at the fictitious crack reloading state, where the condition is , the normal stress is satisfied in (11) . If the normal stress after iteration satisfies , and the opening width of crack surface after iteration meets the condition , then the current state is changed to the fictitious crack extension state.
⑥ If the current state is at the opening state, where the condition is , . If the opening degree of joints after iteration meets , the current state is changed to the fictitious crack unloading state, and the stress at loading point is set .
⑦ If the current state is at the sliding state, where the condition is . If is met after iteration, the current state is changed to the initial closed state.
Having updated the new contact state of one loading increment step, the RHS of (7) is obtained and can be used to get the overall displacement field from (1). Then, the width of the crack surface can be updated, and the new contact force vector can be determined from the stress-displacement curve. The updated contact force and state are taken into the next iteration, and the above process is repeated until both the state and force have converged which will lead to the start of the next loading step.
5. Numerical Experiments
The simulation of crack propagation of Toumi’s [24] beam problem under cyclic load was carried out to verify the proposed scheme. The numerical setup is present in Figure 5 along with the mesh configurations. Similarly, an initial crack is set in BC section, and AB is the pair of contact points. The mesh size near the crack is 1 mm, and other parameters are length = 320 mm, depth = 80 mm, thickness = 50 mm, initial crack length = 40 mm, tensile strength = 5.2 MPa, fracture energy = 34.2 N/m, elastic modulus = 31.6 GPa, poisson’s ratio = 0.2.
The load-displacement curve under monotonous load was shown in Figure 6 and compared against the experimental data [24]. The simulated result lies between the two experimental curves and followed the trend closely. In Figure 7, the load-displacement curve under cyclic load condition was plotted and matches well the experimental data for each loading cycle.
Figure 8 shows the distributions of normal stress along the height of crack under three different loading conditions when the vertical displacement is 0.04 mm and 0.07 mm respectively. No abrupt stress change has been observed as the stress curve is in a continuous and smooth way, and the distribution of normal tensile stress of crack fits nicely with the concrete softening curve adopted. When unloading, the crack surface below the tip of fracture process zone transited gradually from tensile state to compressive state, and the compressive stress increases gradually towards the end. However, if the crack zone is long ( mm), and the opening width of the bottom crack surface is too large, the compressive stress tends to decrease. When reloading, the crack surface below the tip is mainly in the tensile state. Both findings match well general understanding and demonstrate great accuracy and capability of the current scheme.
(a)
(b)
(c)
(d)
(e)
(f)
For simulating the crack propagation, the size of mesh is a factor of great importance. To fully understand the influence of the mesh size on calculated results and investigate the mesh dependence of the suggested method, the load-displacement curves near the crack were plotted using different meshing sizes varying from 1 mm to 8 mm (about 1/40~1/5 of the crack height) in Figure 9. It can be seen that for a meshing size less than 1/10 of the crack height, the result of different meshes shows little change, while if the meshing size and crack height ratio is bigger than 1/5, the load-displacement curve is not smooth anymore. In this case, the trend of the curve in general agreed well with the result of fine mesh despite the big error in the far end.
(a)
(b)
6. Conclusion
(1)The mixed freedom finite element method proposed for contact problems is extended to investigate the crack propagation under cyclic load using the constitutive relation of model I crack. According to the states of the load and the crack, the constitutive relation is characterized by six contact states and solved iteratively. The current scheme is verified against a number of experimental cases and proved to be accurate and reliable. (2)Using the mixed freedom finite element method to simulate the crack propagation, the softening constitutive relation of concrete is implemented directly into the contact iteration. It will not lead to numerical instability during the overlapping process of positive and negative stiffness, which is convenient to adopt the softening constitutive relation of random shape. Moreover, the dependence of the result on the size of mesh is not obvious. (3)As there are limited resources available about the softening relation of mixed type crack, this paper only showed simulations on the propagation of model I type crack. However, the authors believe the proposed method is applicable for other type crack as well. Due to the fact that mixed freedom finite element method needs to preset pairs of contact points in the potential contact area, the application of the current method is confined to problems which the paths of the crack propagation are already known. Future work involves continuing the current work on different types of cracks, and a more detailed analysis of the calculation parameters of the constitutive relations is presented in this paper.
Conflict of Interests
The authors do not have any conflict of interests with the content of the paper.
Acknowledgments
The present work is supported by the National Natural Science Foundation of China (no. 51279050), the National High-tech R&D Program of China (863 Program) (no. 2012BAK10B04), the National Key Technology R&D Program in 12th Five-Year Plan (no. SS2012AA112507), and the Non-Profit Industry Financial Program of MWR (Ministry of Water Resources of China) (no. 201301058).