Abstract
A unified cohesive zone model is proposed. By introducing a parameter which denotes the ductility of the material, the proposed model can be used to simulate both brittle and ductile fracture. No fitting experimental data is needed in calibrating the parameters of the proposed model, and this leads to a high predictive ability of the model. Several numerical examples for brittle and ductile materials are given using finite element method. The good agreement between numerical results and experimental data proves the applicability and accuracy of the proposed cohesive zone model.
1. Introduction
Fracture is one of main failure modes in engineering. Fracture mechanics techniques, including linear elastic fracture mechanics and nonlinear fracture mechanics, can be used to evaluate the size of a crack to cause catastrophic structural failure and the time a crack of known length reaches to this value in the current service environment. However, linear elastic fracture mechanics break down when the fracture process zone is too large relative to the crack length or other structural dimensions. Also, nonlinear fracture mechanics concepts like that of crack tip opening displacement (CTOD) or that of J-integral are not easy to apply in engineering analysis. One method to account for these problems is the cohesive zone model. In this method, a zone of nonzero traction is placed at the crack tip. This concept provides a convenient way of avoiding problems with crack tip singularities and can be used in many situations where conventional fracture mechanics do not work.
The idea of cohesive zone model can be traced to the strip yield models of Dugdale [1] and Barenblatt [2]. In this way, the unrealistic continuum mechanics stress singularity at the crack tip can be avoided. Several merits make the cohesive zone model popular in fracture mechanics. It seems to mimic many process zone effects in engineering materials. The fracture process, isolated from surrounding continuum constitutive material, is very simple: traction versus separation. It is easy to implement the cohesive zone model in the finite element method using interface elements or specified interface. The cohesive zone model has been successfully used in solving the complicated fracture problems such as mixed mode fracture [3], fatigue fracture [4, 5], and dynamic fracture [6, 7].
The choice of a constitutive law for the cohesive zone is the most delicate aspect of this approach to fracture mechanics. Due to the small scale of the cohesive zone in most materials, it is experimentally quite difficult to determine the precise nature of the constitutive behavior in the cohesive zone. Thus, it has become commonplace to postulate a phenomenological form for the cohesive zone constitutive model. Postulating the form of the cohesive zone model is mainly based on the fracture characteristic of the material, for example, ductile or brittle modes. In the literature, several cohesive zone constitutive models can be found. Typical cohesive zone laws include linear decreasing form suggested by Hillerborg et al. [8], cubic and exponential forms by Needleman [9, 10], constant form by Yuan et al. [11], bilinear form by Geubelle and Baylor [12], trilinear form by Tvergaard and Hutchinson [13], and polynomial form by Chaboche et al. [14].
The object of the present work is to provide a unified cohesive zone model for both brittle and ductile fractures. The model is limited to mode I crack. The paper begins with a description of the unified cohesive zone model and the parameter calibration of the model. After that, the numerical implementation of the cohesive zone model is introduced. This is then followed by solving some examples of brittle and ductile fractures to demonstrate the applicability and accuracy of the model. Finally, the concluding remarks are made based on the study.
2. Cohesive Zone Model
During crack propagation, a fracture process zone exists in front of the crack tip, where microvoids initiate, grow, and finally coalesce with the main crack, as shown in Figure 1 (a). The material response in this fracture process zone differs greatly from that of the surrounding bulk material. In fracture process zone, the material degradation is obvious. The conventional elastic-plastic constitutive relation used in bulk material is not suitable for this local material degraded region. When formulating analytic solutions to problems involving the fracture process zone, a new constitutive relation is necessary. Compared with the other characteristic dimensions of the region, the thickness of the fracture process zone is very small. Therefore, it is reasonable to assume that the fracture process zone can be approximated by a zone of cohesive tractions, which has zero thickness in the undeformed configuration. In this case, the fracture process zone is termed a cohesive zone which consists of two fictitious cohesive surfaces. The cohesive surfaces are connected by cohesive traction and can be separated during deformation as shown in Figure 1 (b). The relation between cohesive traction and displacement jump of cohesive surfaces is usually called traction separation law. In contrast to the standard concept in continuum mechanics, cohesive zone model allows for a displacement jump inside the material and for separation of cohesive surfaces which finally leads to failure of the material. By introducing cohesive zone model, the material's behavior is split into two parts: deformation and separation. For bulk material, deformation is accounted for by conventional continuum material model, whereas the fracture process zone is modeled by the traction separation law.

Fracture process zone (a) and schematic of the cohesive zone (b) in front of the crack tip under the loading condition.
The cohesive zone model is a phenomenological model and there is no direct evidence for the formulation of traction separation law in the model. Various forms have been proposed since its initiation, as shown in Figure 2 [15]. A challenge of using cohesive zone model is to choose a proper traction separation law for different materials. As mentioned above, there are various forms for different materials They can be classified into two categories based on the model features, for example, zero traction and nonzero traction at the initial zero separation. The typical cohesive zone models of these two categories are Needleman's and Hillerborg's models as shown in Figure 3. For ductile materials, the exponential form proposed by Needleman [10] is one of the most popular traction separation laws. In this model, the traction first increases drastically to a peak value T0, called the cohesive stress, and then decreases to zero at a critical separation Δ c . For brittle materials, a linear decreasing form suggested by Hillerborg et al. [8] is used widely. The model starts with an infinite stiffness until the cohesive stress T0 is reached. The area under the traction-separation curve is called the cohesive energy Γ. Recent atomistic simulations show that even for ductile materials [16], the cohesive zone model should have a nonzero critical traction. The separation cannot occur when the traction is below this critical traction. Jin and Sun [17] also conclude that the cohesive zone model must have a nonzero traction at the initial zero separation, and the cohesive traction should first increase quickly to the cohesive stress and finally decreases with further increasing separation. Based on this knowledge, we propose a new cohesive zone model, in which a parameter representing the ductility of the material is included to determine the critical traction. The proposed model is focused on the mode I crack problem.

Various traction separation laws for mode I crack [15]: (a) rigid plastic law, (b) linear elastic brittle law, (c) triangular law with initial linear elastic regime followed by a linear softening regime, (d) an arbitrary traction law, often approximated with a cubic function, (e) linear softening law, (f) nonlinear softening law.

For mode I fracture, the constitutive behavior of the cohesive zone is derived through a potential function Φ in the form
where T n is the normal traction and g n is the normal separation. The equation for the potential function Φ is proposed to be the form
with
where Δ0 is the cohesive length, which is the normal separation between the crack surfaces corresponding to the cohesive stress T0 for mode I crack problem. λ is a model parameter. The relation between the traction and separation is obtained by the substitution of (2) into (1):
The parameter λ can be defined to be the ratio of yield stress to tensile strength of the material. When the value of λ ranges between 0 and 1, the model is applicable to the materials with different ductilities. Figure 4 displays the cohesive zone model with different values of the parameter λ but the same cohesive stress and cohesive energy. By setting λ = 0, the cohesive zone model returns to with Needleman's cubic traction separation form model.

The cohesive zone models with different parameter λ but the same cohesive stress and cohesive energy Γ.
Until now, there is no unified method to determine the cohesive stress according to the author's knowledge. However, it is commonly considered to be dependent on the material property, that is, brittle or ductile. In this paper, the determination of the cohesive stress is based on the global measurement results for brittle material [18–20] and the micromechanical study of void growth for ductile material [21–26]. For brittle material, it fails without obvious plastic deformations, and the failure is dominated by the maximum principal stress. The material behaves linearly elastic before the maximum principal stress reaches the tensile stress. Once the tensile stress is exceeded, the material begins to fail. Based on this knowledge, the cohesive stress for the brittle materials can be taken as tensile strength of the material. For ductile material, both numerical and experimental results show that the cohesive stress strongly depends on the condition of mechanical constraint present at the crack tip [21]. The constraint state can be represented by the stress triaxiality, that is, the ratio of hydrostatic stress to equivalent stress χ = σ h /σ e [23]. The relationship between cohesive stress and the stress triaxiality can be determined by studying the material failure uncoupled from the crack tip problems using the so-called unit cell method [22]. In the unit cell method, the material is regarded as a periodic array of identical cells, and a representative unit cell of material under consideration is picked out to be investigated. The method is initially used to determine the constitutive relation of composite materials at the microscale [27] and recently used by researchers to determine the triaxiality dependent cohesive stress for cohesive zone model [28, 29]. To determine the relationship between cohesive stress and stress triaxiality, a representative unit cell in front of the crack tip can be studied separately from the macroscopic specimen. The constitutive equation of the material for volume element is based on the Gurson-Tvergaard-Needleman (GTN) model [24–26], which describes the void growth and nucleation process in ductile materials. Using the unit cell method, the maximum load carrying capacity can be captured during the void growth under the given stress triaxiality, and the value of the maximum load carrying capacity is equal to the cohesive stress.
The cohesive energy Γ for both ductile and brittle materials is identical to the fracture energy, that is, the energy needed to create unit area of new crack surface. It can be experimentally measured. Once the parameters, that is, ductility λ, cohesive stress T0, and cohesive energy Γ, are settled, the cohesive length Δ0 is subsequently determined.
3. Numerical Implementation
The cohesive zone model is usually implemented into finite element analysis through a zero-thickness interface element, that is, cohesive element. The formulation of cohesive element is described in the following. The cohesive element is made up of two surfaces: the upper and the lower. In undeformed condition, the two surfaces of the cohesive element lie together. Under external load, the two surfaces separate but are connected by cohesive traction. The traction is calculated at the integration points according to the traction separation law. When the separation reaches a critical separation, the element fails at this point, and the contribution of the integration point to the stiffness vanishes. In the study, the bulk material is modeled by 4-node quadrilateral elements. Based on the compatibility requirements in the finite element formulation, the corresponding cohesive element is a 4-node interface element with two nodes at each surface. The cohesive element condition and its integration point position are illustrated in Figure 5. When the separation across the crack reaches a critical separation Δ c , the solid elements connected by this cohesive element are disconnected, which implies that a real crack propagates to this element.

Cohesive element in two-dimension condition. The circles represent the element nodes, and the solid circles represent the Gauss points.
The stiffness matrix of the cohesive element and finite element equations can be derived from the principle of virtual work. The principle of virtual work reads
that is, the internal virtual work δΠ i is equal to the external virtual work δΠ e . The symbol δ denotes the variation of a quantity. The internal virtual work for the cohesive element is defined by
where
Inserting (7) into (6) leads to
with
where superscript ()T denotes a transpose operation.
The external virtual work δΠ e is defined by
where
with
where
Combining (5), (8), (9), and (11) gives
Since the cohesive traction is nonlinear functions of the separation, the incremental procedure is needed. The traction rate can be expressed in terms of nodal displacement velocities as
where
where g
t
and T
t
are the tangential separation and tangential traction, respectively. There are two degrees of freedom for two-dimensional problem, that is, normal and shear degrees. However, for mode I crack, all components related to shear degree are zero. Therefore, the cohesive stiffness
Substitution of (14) into (13) leads to the tangent stiffness matrix
with
4. Numerical Examples
In this section, four numerical examples of crack growth are carried out for two brittle materials and two ductile materials to verify the applicability and accuracy of the proposed cohesive zone model. Two brittle materials are two kinds of concrete and two ductile materials are steel 20 MnMoNi 55 and aluminium alloy AL5083. Cohesive zone model defined in (4) is implemented as the cohesive element through the user-defined subroutine UEL in the commercial finite element code ABAQUS.
4.1. Crack Growth in Brittle Materials
Two kinds of concrete, named as concrete 1 and concrete 2, are taken as examples of the brittle materials. Crack growth in three-point bending test is simulated. Experimental results taken from Fathy's work [30] are used for comparison with simulating results. The geometry, loading, and boundary conditions of three-point bending specimen are shown in Figure 6. The geometric parameters are S = 300 mm, L = 500 mm, H = 100 mm, a0 = 33 mm, and thickness B = 100 mm. The finite element meshes for the specimen are given in Figure 7, including 1368 elements and 1451 nodes. To improve the computational efficiency, only the potential damage zone along the expected crack path is covered by the cohesive elements, and the rest of the region is occupied by 4-node linear elements.

Geometry, loading, and boundary conditions of three-point bending.

Finite element meshes of the three-point bending specimen.
4.1.1. Crack Growth in Concrete 1
The material properties of concrete 1 are Young's modulus E = 47.7 GPa, Poisson's ratio ν = 0.2, tensile strength σ b = 6.18 MPa, and the cohesive energy Γ = 0.126 N/mm [30]. The parameter λ = 1 is employed. Cohesive stress T0 = 6.18 MPa which is equal to the tensile strength. The plane strain condition is assumed.
The experimental and simulating results of force-crack mouth opening displacement (F-CMOD) curves are summarized in Figure 8 (a). It can be seen that the numerical results based on the proposed model have a good agreement with the experimental data. Moreover, the prediction of Fathy et al. [30] is also displayed in Figure 8 (a). It can be found that the simulating results are very close to Fathy's results. Figure 8 (b) shows the proposed model and Fathy's model. It is obvious that the cohesive stress in Fathy's model is higher than the proposed model which causes the slight difference of the predictions in Figure 8 (a).

The simulating F-CMOD curve in comparison with the experimental as well as Fathy's results for concrete 1 [30] (a), and the comparison between proposed model and Fathy's model (b).
4.1.2. Crack Growth in Concrete 2
To further study the feasibility of the proposed model in brittle materials, another simulation has been done for concrete 2. The material properties of concrete 2 are Young's modulus E = 45.1 GPa, Poisson's ratio ν = 0.2, tensile strength σ b = 5.82 MPa, and the cohesive energy Γ = 0.121 N/mm [30]. Parameter λ = 1 and cohesive stress T0 = 5.82 MPa which is equal to the tensile strength. The same plane strain condition is adopted as concrete 1.
Both simulating and experimental F-CMOD curves are shown in Figure 9 (a), together with the prediction of Fathy et al. [30]. It can be found that the simulating results fit the experimental results very well in the initial part as well as the peak load value. For the postpeak response, Fathy's model gives a better result than the proposed model. Figure 9 (b) displays the comparison of the two models. It is clear to see that softening curve in Fathy's model decreases much steeper than the proposed model which could lead to the deviation of the simulating results in the postpeak response in Figure 9 (a).

The simulating F-CMOD curve in comparison with the experimental as well as Fathy's results for concrete 2 [30] (a), and the comparison between proposed model and Fathy's model (b).
4.2. Crack Growth in Ductile Materials
In order to study the feasibility of the proposed model in ductile materials, two different ductile materials, that is, steel 20 MnMoNi 55 and aluminium alloy AL5083, have been employed. Crack growth tests for both materials are carried out in standard compact tension (CT) specimen which is displaced in Figure 10. Figure 11 shows the finite element meshes used in the calculation, with 3877 elements and 3996 nodes. Same as the mesh arrangement of the three-point bending specimen in Section 4.1, the potential damage zone along the expected crack path is covered by the cohesive elements, and the rest of the region is occupied by 4-node linear elements. The experimental results are documented in [31] for 20 MnMoNi 55 and [32] for AL5083.

Geometry of standard CT specimen. V is the location of displacement measurement in the test, B is the thickness, and Bnet is the thickness of the net section.

Finite element meshes of the CT specimen.
4.2.1. Crack Growth in Steel 20 MnMoNi 55
The geometric parameters of the CT specimen for steel 20 MnMoNi 55 are W = 100 mm, a0 = 60 mm, B = 10 mm, and Bnet = 8 mm. The material properties are Young's modulus E = 210 GPa, Poisson's ratio ν = 0.3, yield stress σ y = 465 MPa, tensile strength σ b = 600 MPa, and the cohesive energy Γ = 120 N/mm [31]. The parameter λ calculated from the ratio of yield stress to tensile strength is 0.78. The plane strain condition is assumed in the simulation.
As mentioned above, the cohesive stress can be determined by studying the damage process at the microscale using the so-called unit cell method. At a certain point of the structure, the stress triaxiality is varied with the increase of external loading. For cohesive zone model, it results in a change of cohesive stress which depends on the stress triaxiality. The change of cohesive stress will not only make the simulation complicated but also lead to a convergency problem [29]. In this study, the maximum value of the stress triaxiality at the initial crack tip is used for evaluating the cohesive stress in unit cell method. The same treatment has been successfully used by other researchers [28, 33]. To determine the triaxiality dependent cohesive stress, the stress triaxiality at the initial point is studied first using the finite element method, in which no crack propagation is considered. The response of the stress triaxiality as the loading increases is shown in Figure 12 (a) and the parameter ρ represents the load fraction. The maximum value of stress triaxiality which represents the maximum constraint state that can endure during the whole plastic deformations is then used in the unit cell method. Once the stress triaxiality is settled, the unit cell method is used to determine the cohesive stress. Following the configuration in [29], a two-dimensional axisymmetric element is used in the unit cell method. For the material considered here, the GTN model parameters are taken from Pirondi's work [34]. Loading under various constant applied stress ratio, β = σ1/σ2, is considered. σ1 is applied in the horizontal direction while σ2 is applied in the vertical direction. The triaxiality of the element can be calculated as
For a given β, the applied stress and the vertical displacement at the top of the element can be evaluated. Since the maximum value of σ2 is taken as the cohesive stress, the relationship between the cohesive stress and the stress triaxiality can be obtained. Then the maximum value of the stress triaxiality at the initial crack tip can be used to determine the cohesive stress. The relationship of the applied stress and the displacement in the vertical direction is shown in Figure 12 (b), where u is the displacement in the vertical direction and D is the length scale associated with void growth in the unit cell method. As explained in Section 2, the value of the maximum load carrying capacity is the cohesive stress, namely; T0 = 3.44, σ y = 1600 MPa.

Determination of the cohesive stress for 20 MnMoNi 55. (a) the response of stress triaxiality as loading increases. (b) The relationship of the applied stress and the displacement in the vertical direction determined by unit cell method at χ = 2.2.
Once the cohesive parameters are determined, the crack propagation simulation is carried out. The simulating as well as the experimental force-displacement (F-V) curves are displayed in Figure 13 (a). The numerical prediction obtained by Cornec et al. [31] is also displayed. For whole force-displacement curve, the proposed model and Cornec's model [31] almost give the same values and match the experimental results accurately. Figure 13 (b) shows the proposed model and Cornec's model. The difference between the two models is obvious. It seems that for the ductile materials, the simulating results are not sensitive to the form of the cohesive zone model.

Comparison of experimental [31] and predicted F-V curves for steel 20 MnMoNi 55 (a), and comparison of proposed model and Cornec's model (b).
4.2.2. Crack Propagation in Aluminium Alloy AL5083
The CT specimen used for aluminium alloy AL5083 has different geometric parameters from that for steel 20 MnMoNi 55, which are W = 50 mm, a0 = 25 mm, and B = Bnet = 3 mm. The material properties of AL5083 are Young's modulus E = 68 GPa, Poisson's ratio ν = 0.33, yield stress σ y = 243 MPa, tensile strength σ b = 350 MPa, and the cohesive energy Γ = 10 N/mm [32]. In the simulation, plane stress condition is adopted. Parameter λ = 0.69 which is calculated by the ratio of yield stress to tensile strength. Using the same method mentioned in Section 4.2.1, the stress triaxiality is calculated and depicted in Figure 14 (a). The maximum value of stress triaxiality is 0.63, which will be used in the unit cell method. The GTN model parameters in the unit cell method for this material are determined by Gao et al. [35]. The relationship of the applied stress and the displacement in the vertical direction is determined by unit cell approach and shown in Figure 14 (b). From Figure 14 (b), the cohesive stress is determined as T0 = 2.3, σ y = 560 MPa.

Determination of the cohesive stress for AL5083. (a) The response of stress triaxiality with the increase of load fraction. (b) The relationship of the applied stress and the displacement in the vertical direction determined by unit cell method at χ = 0.63.
The comparison between the simulating and experimental results, as well as the prediction by Scheider et al. [32], is shown in Figure 15 (a). The comparison between proposed model and Scheider's model is given in Figure 15 (b). For whole force-displacement (F-V) curve, the proposed model presents a better agreement with the experimental results than Scheider's model. Same as in the case of steel 20 MnMoNi 55, the proposed model shows a big difference with Scheider's model, as shown in Figure 15 (b).

Comparison of experimental [32] and predicted F-V curves for AL5083 (a) and comparison of proposed model and Scheider's model (b).
5. Conclusions
In the present work, a unified cohesive zone model has been proposed, which is applicable for both brittle and ductile materials. To verify the feasibility of the proposed model, several numerical examples have been carried out, and the following conclusions are obtained.
Introducing a new parameter into conventional cohesive zone model to characterize the ductility of the material can lead to a unified cohesive zone model which is suitable for both brittle and ductile materials.
In the present investigation, all the parameters in cohesive zone model are calibrated independent on the cohesive zone model, that is, no fitting experimental data is needed, and this leads to a high predictive ability of the model.
The simulating results have a good agreement with the experimental data, which proves the applicability and accuracy of the proposed model. Moreover, the form of the cohesive zone model is not crucial for ductile materials because the different cohesive zone models can predict nearly same crack growth behavior for ductile materials.
Footnotes
Nomenclature
Acknowledgment
This study is based upon the work supported by the National Natural Science Foundation of China (Grant no. 51005020).
