INVESTIGATION OF THE MIXED-MODE FRACTURE IN DELAMINATION TESTS : NUMERICAL SIMULATIONS USING COHESIVE ZONE AND PARTITIONING METHODS

Original scientific paper Delamination (fracture) tests have been numerically investigated using various cohesive zone properties. The test utilises asymmetric and symmetric double cantilever beam specimens loaded with bending moment. Energy release rate contributions from mode I and mode II fracture are calculated using a global and local approach. Mode-mixities results are presented and analysed. The numerical partitioning results for different configurations are compared to two analytical partitioning theories, namely, after Williams and after Hutchinson and Suo. Opposite to these theories, partitioning is observed to be dependent on cohesive zone properties.


Introduction
Interlaminar fracture is one of the most important failure modes for many modern materials arranged in layers.Mixed-mode fracture is often observed in delamination, where fracture toughness can be vary depending on mode-mixity, i.e. amounts of different fracture modes are present (mode I and mode II are mostly considered).It is therefore of great importance for design consideration of these materials to define the interlaminar toughness for full range of mode-mixities, from pure mode I to pure mode II fracture.
The delamination tests for composite laminates and adhesive joints extensively utilise beam like geometries.The pure mode I and II fracture toughness calculations from the experimental results are pretty much straightforward, however, there is much confusion about analysis of test configurations with mixed mode fracture.To calculate contributions from mode I and mode II fracture one can employ analytical or numerical partitioning methods, each of which suffer from a number of uncertainties and can produce different results depending on choice of a theoretical approach, numerical model, etc.
The delamination test under consideration in this work utilises double cantilever beam specimen loaded with uneven bending moments (DCB-UBM) [1,2], also known as the fixed-ratio mode-mixity (FRMM) test.The test is simulated in the commercial finite element (FE) software package Abaqus [3] using cohesive zone model.The fracture energy (strain energy release rate) is calculated and partitioned numerically from the simulation results and mode-mixities are compared to the analytical partitioning theories from Williams [4] and Hutchinson and Suo [5].
The present work is part of an ongoing investigation of the mixed-mode fractures in beam-like geometries under the coordination of European Structural Integrity Society, Technical Committee 4 (ESIS TC4).Investigations in the project are conducted as Roundrobin exercises.The ultimate goal of this project is a new testing protocol with recommendations for the accurate determination of mode-mixity in all beam-like geometries.In the first phase, the numerical simulations are conducted without damage development and modemixity is calculated by employing virtual crack closure technique and interaction domain integral [6].Here, as a part of the second phase, damage development in simulations is involved.The present work is continuation and expansion of the previous works by authors [7,8].

Problem description
Total of eighteen FE configurations for simulation of delamination tests are set-up using different test and simulation parameters.These are summarised in Tab. 1 and more details about settings are given in the following subsections.

Test configuration
The DCB specimen geometry and the test configuration used in this work are illustrated in Fig. 1.The specimen consists of two arms (beams) with thickness h 1 for the upper arm and h 2 for the lower arm.The arms are joined together over a half of the length at the right-hand side and fixed at the right end.They are separated at the left-hand side, thus forming a pre-crack.The value for h 1 is kept constant (3 mm) and h 2 value is varied (0,3/3/30 mm) to provide different range of modemixities (symmetric and asymmetric DCB).The upper arm is loaded with a bending moment M, whereas the lower arm is left free.The beam material is linear elastic, isotropic with the modulus of elasticity 50 GPa and Poisson's ratio 0,38.

Cohesive zone model
Fracture process is simulated by using cohesive element formulation proposed in [9] with the capability of dealing with crack propagation under mixed-mode loading.Response of the elements is defined in terms of traction versus separation, providing coupled and uncoupled behaviour between the traction and separation components and different models and criteria for defining damage initiation and evolution.An overall evolution of damage is defined with the scalar damage variable, D(sdeg in Abaqus), representing a stiffness degradation.The variable monotonically evolves from 0 to 1 in a fracture process.
In this study the uncoupled initial linear elastic behaviour is used, which for a two-dimensional problem can be writen as: , 0 where t is the nominal traction vector (normal and shear tractions), K is the elasticity matrix and ε is the nominal strain vector.Since nominal strains are defined as separations in two directions ( , ) n t δ δ divided by the nominal thickness that is equal to 1 by default, nominal strains are equal to the separations.Uncoupled elasticity matrix is defined with arbitraly high set value of penalty stiffness with K nn = K tt = 10 15 MPa.
Element damage initiation is defined using the quadratic nominal stress criterion [10]: Dependence of the fracture energy on the mode mixity is given with the power law fracture criterion [11]: where G I and G II refer to the work done by the traction in the normal and the shear directions, respectively.Using power parameter α=1 Eq. (3) reduces to a linear form as employed in this research: The arbitrary choice of equal critical energy and equal inter-laminar strength for two fracture modes is made to avoid problems reported in [12,13] regarding the type of cohesive zone laws implemented in Abaqus under mixed-mode fracture.

FE model mesh and set-up
According to the test configuration given in Fig. 1, the DCB-FE model is made from two separate beams (parts) with coincident nodes connected along a half of the length (dashed line) with zero-thickness cohesive elements (having nominal thickness equal to 1).The other halves of the beams form the pre-crack and have unconnected coincident nodes.No surface interaction in pre-crack region is defined since the two pre-crack surfaces are separated immediately at the test initiation.Abaqus CPE4 (4-node bilinear plane strain quadrilateral) elements are used for modelling beams and COH2D4 (4node two-dimensional cohesive) elements for modelling cohesive zone elements.Uniform and non-uniform meshes are modelled for the different configurations where element sizes are chosen to ensure an adequate mesh density in a fracture process zone (FPZ), with 10 elements in the zone being the targeted minimum value [7].The FPZ refers to a zone comprising cohesive elements that have reached damage initiation criterion (2).A uniform mesh (throughout each beam) is used in configurations where a FPZ is larger and a non-uniform in configurations where it is smaller in order to rationalise number of elements in a model and CPU usage in a simulation.In the non-uniform mesh configurations, the finest mesh is used in the region around the initial crack tip with the length of the fine mesh region in the cohesive side being twice shorter than that in the pre-crack side.Bending moment is simulated via rotation applied incrementally at the end of the upper beam which is set to be rigid (nodes at the end line are connected into the rigid body).The rotation is chosen because it ensures constant moment loading during crack propagation and numerical stability of simulation.
Convergence difficulties are often encountered in numerical simulations of delamination in beam-like geometries due to softening behaviour and stiffness degradation of a model.To overcome these difficulties the viscous regularisation of the cohesive element constitutive equations is included in the simulations using the value of 10 −4 for the viscosity parameter [7].

Analytical partitioning
Extensive research has been carried out in the last 30 years in area of mixed mode fracture in layered materials producing several analytical partitioning theories that are used to predict mixed mode partitions as a function of geometry (h 1 /h 2 ) and loads, but there is still much confusion around their validity and application in practice (a more detailed review can be found in [6,14]).Two pioneering theories, namely, Williams analytical solution [4], obtained by global approach using beam theory, and Hutchinson and Suo semi-analytical solution [5], based on linear elastic fracture mechanics (LEFM) and local approach, are selected here for the analysis because they can be defined as two possible extremes [6,15].Here, the term 'global' refers to the entire region mechanically affected by the crack, and 'local' to the crack tip.
Fig. 2 shows the considered problem with the original notation from authors being modified to correspond to the notation used here (as in Fig. 1).It is a homogenous isotropic beam (modulus of elasticity E, moment of inertia I) containing a delamination (or crack) at specified distance from the top or bottom surface, separating it into two arms.The delamination is forced with the bending moments acting on the arms (in the two analyses the lower arm moment loading is applied in opposite sense).Both methods produce equal results for total energy release rate, although the calculation formulas are given in different forms.The fracture energy can be derived by using beam theory or alternatively J-integral [16] and is given here in the form [5]: where E' is the effective Young's modulus defined as: and ν is Poisson's ratio.

Williams solution
For the Williams solution [4], the strain energy release rate mode I and mode II partitions are given by: ( ) 1 16 1 ( ) where M I and M II represent the contributions from bending moments for the fracture modes given by 1 2 1 1 The contributions are determined using conditions that a pure mode II fracture is obtained when the curvature in the beams above and below the crack is the same (no crack opening between the crack faces), and that in pure mode I fracture moments acting on the beams are equal and opposite (no sliding between the crack faces).

Hutchinson and Suo solution
In this solution, the mode I and II fracture energy partitions are given by [5]: 1 cos sin , and where ( ) ( ) ( ) ( ) ( ) ) The above calculation is only valid for r<1, see Eq. ( 16), but it may be used for r>1 (as is considered here) simply by swapping the upper and bottom arm, i.e. values of h 1 and h 2 , as well as M 1 and M 2 [6], and by using 1/r in calculations that would correspond to r=h 2 /h 1 .

Numerical partitioning
Although energy release rate partitions are included in cohesive element damage evolution Eqs.(3), their individual values are not available as Abaqus simulation output variables.Hence, they must be calculated by integrating (numerically) outputs for stresses and strains (tractions and separations).The two implemented methods are here classified as 'local approach' and 'global approach' based on the size of a damage zone or model included in the calculations (one or more cohesive elements).Therefore, in the numerical partitioning, the term 'global' means 'integration of energy over fracture process zone', whereas 'the local approach' considers 'integration of energy going into a single cohesive element'.

Local approach
In the local approach, fracture energy calculation is based on tracking strain energy going into a single cohesive element by integration of traction-separation curves.The mode I and mode II energy component integrals are given by [9,17]: 0 0 , and where δ nm , δ tm are maximum (final) opening and shearing displacements, δ n , δ t are opening and shearing displacements and σ, τ are normal and shear stresses (tractions n t and s t ).These integrals are numerically calculated using history output for a cohesive element integration point via trapezoidal rule: ( ) ( ) , and 2 , 2 where n is the total number of simulation increments, σ, τ are stress components (Abaqus identifiers s22 and s12, respectively) and δ n , δ t nominal strain components (Abaqus identifiers ne22 and ne12, respectively) for a given increment i(i+1), Fig. 3.As already mentioned, the values of the strain components are equal to the values of the displacements since the nominal thickness of cohesive elements is set to 1.The cohesive element used in this study has two integration points (P 1 and P 2 ) and the final values of fracture energy components for a single element were obtained by averaging values from integration points which generally have different loading history (see Section 5): , and 2 2

Global approach
The global approach states that once damage region is fully developed, self-similar crack propagation will exist, i.e. every cohesive element will have the same loading history; hence, in the FPZ there will be cohesive elements from all stages of a fracture process (loading history).The following integrations can then be performed along the cohesive surfaces in order to obtain the global mode I and mode II energy release rates [7,15]: where l is the length of the FPZ, δ n , δ t are opening and shearing displacements of the cohesive elements, σ, τ are normal and shear stresses and axis x coincides with a crack propagation direction, Fig. 4. Mathematically, this is equal to the mode decomposed J-integral [12] in the absence of any crack tip singularity.
Above integrals are numerically calculated using field outputs for cohesive elements via trapezoidal rule: ( ) ( ) , and 2 , 2 where n is the number of cohesive element nodes on bottom or top surface belonging to cohesive elements included in integration, σ, τ are stress components and δ n , δ t are nominal strain components, Fig. 4.

Figure 4
Global approach numerical integration area and tractions (stresses) and separations (strains) in FPZ element nodes.
For practical reasons, integration is done over the almost whole crack path, excluding only few cohesive elements at the fixed end (due to the influence of the fixed boundary condition).Thus, contributions from the cohesive elements that are fractured or still not damaged, are zero or near zero and do not have significant influence on the integration results.Stresses and strains are averaged to the nodes between adjacent cohesive element integration points, and they are taken in the last increment before the cohesive element in the current crack tip fractures i.e. in the increment when it reaches maximum damage.

Results and discussions
The accuracy of all simulations is evaluated by analysing deviations of calculated total fracture energy release rates from prescribed critical fracture energy values and by monitoring reaction moment in the upper beam boundary, where the rotation is applied.The highest energy error calculated was 1,05 % with global and 1,00 % with local approach.In all simulations, constant trend of reaction moment after delimination onset was registered with very small oscillations, confirming steady-state crack propagation.Fig. 5 shows reaction moment change with rotation application (time) for three selected configurations in which constant moment trend after crack propagation onset can be observed.The other configurations have similar pattern of the moment vs. rotation curves.Values of bending moment during crack propagation are averaged and compared with the analytical prediction given by Eq. ( 7).Tab. 2 summarises this comparison, showing very good agreement beetwen results for all simulations.The numerical results also confirm that the moment depends on a geometry, beam modulus of elasticity and critical fracture energy, but not on other cohesive properties.Observed FPZ and crack propagation lenghts are given in Tab. 3. Element damage initiation and complete failure instances are registred using damage variable D. FPZs have nearly constant length throughout delamination process, as expected for the steady-state crack propagation.The growing trend of FPZs length is observed with increase in both critical fracture energy (G c ) and beam thickness ratio (h 1 /h 2 ), and decrease in inter-laminar strength (t o ), as shown in Fig. 6.The targeted number of 10 elements in the FPZ has not been achieved in only three configurations (6 being the lowest number, see Tab. 3) due to the very small damage zone sizes (less than 1 mm).Components of the fracture energy are calculated using local, Eq. ( 19) and (20), and global approach, Eq. ( 22), along the crack propagation direction, in order to monitor change in mode-mixity with crack propagation.Results for all simulations are given in Fig. 7, where mode-mixity ratio is given in relation to the crack propagation length; it refers to the results of local approach calculations for cohesive elements which collapse occurs at specific crack propagation length and the results of global approach calculations made at the instance (FE analysis time increment) before their collapse.
It can be noticed that mode-mixity is almost constant throughout crack propagation when global approach is used.The reason for this lies in the fact that the global approach calculations (integrations) are performed from the instance when the damage region is fully developed and the influence of the initial crack tip singularity can be neglected.On the other hand, the mode-mixity varies when the local approach is used; the highest values are at the initial crack tip with asymptotic convergence along the crack propagation path with different converging zone lengths to reach converged solution.This variation is a result of different loading histories of cohesive elements; elements near the initial crack tip are affected by the initial crack tip singularity, whereas the other elements become damaged inside the fracture process zone.Comparing the FPZ lengths given in Fig. 6, and the converging zone lengths (approximated from the modemixity curves) given in Fig. 7, one can observe that the two are related.The converging zone length is closer to the FPZ length for the configurations with mode-mixities having larger contribution of mode I.After the local approach results converge, the two approaches give very close results, agreeing more in the configurations with more cohesive elements, which provide better numerical accuracy.
Dependency of the mode partitioning on cohesive zone properties is observed, confirming the findings reported in [1,15].The dependency is more pronounced in asymmetric beam geometries (h 1 /h 2 ≠1), while for the symmetric case (h 1 /h 2 =1) mode-mixity is nearly independent on cohesive properties, as suggested by both analytical solutions.
Next to investigate are the mode-mixity histories for cohesive elements (or more precisely, their integration points) along the crack path against the damage evolution presented as a percentage of the total fracture energy, as shown in Fig 8 .The diagrams are given for six configurations with four representative cohesive elements: the initial crack tip element, an element from converging zone, and two elements from converged zone.In general, histories of the first and the second integration point of a cohesive element are different, which is expected since the points belong to the same element at different locations along the same line of the element where cohesive zone law is enforced.
As a consequence, the first integration point has greater mode I contribution than that of the second integration point.However, it can be shown that the average values at the full cohesive element damage are around steady-state value for all positions (dotted line in graphs).On the other hand, the histories of the corresponding integration points between different elements are very similar.In fact, they are identical for the elements belonging to the converged zone, meaning that all elements are going through the same damage evolution process.The only exceptions to aforementioned features are the initial crack tip elements and elements in the converging zone.Fig. 9 shows the mode I contribution for all simulation configurations together with two analytical solutions.As the global and the local numerical approach give very close results, only the global results are plotted.The mode-mixity curve for the considered case (Fig. 1) according to Williams is obtained by Eq. ( 12), and the one according to Hutchinson and Suo by Eq. ( 13) ÷ (17).It really seems that the two analytical solutions form the upper and lower bound for the mixed mode partitioning, confirming conclusions given in [1,15].As already said, exceptions are the symmetrical configurations for which the contribution of mode I is nearly independent on cohesive properties.It can also be observed that results for configurations with higher fracture energy and lower inter-laminar strength are closer to Williams solution, and those with lower fracture energy and higher cohesive strength to Hutchinson and Suo solution.Moreover, when comparing results from Fig. 6 and Fig. 9, one can observe that the position of the mode I contribution for a specific configuration between the two analytical curves is directly related to the FPZ length, with the configurations having larger values being closer to the Williams solution and those with smaller values to the Hutchinson and Suo solution.This captures a unique dependence of a mode-mixity for specific geometry (h 1 /h 2 ) on a FPZ length which is more clearly shown in Fig. 10.

Summary and conclusions
The global and local approach in numerical partitioning produce similar results.The local one provides better insight into change in local mode-mixity with crack propagation, but the global one provides results that are not dependent on crack propagation length, once a damage zone is fully developed.
Opposite to the analytical theories, numerical results show mode-mixity dependence on cohesive model properties (and FPZ length), except for the symmetrical geometries.This shows limitations of the exiting modemixity partitioning theories for application in practice and a need for a new, more accurate, property-dependent partitioning solution.
The implemented cohesive zone model has some limitations but it is believed that the observations made using it can be general.For more accurate numerical solutions, a physically more realistic cohesive model should be employed.

Figure 3
Figure 3Local approach numerical integration of traction-separation curve.

Figure 6
Figure 6 Effects of cohesive zone properties and DCB geometry on FPZ length

Figure 7 Figure 8
Figure 7 Mode-mixity with crack propagation for different configurations, according to the local and global approach

Figure 9 Figure 10
Figure 9 Mode-mixity for the investigated test obtained numerically (global approach) and analytically

Table 1
Test and simulation configurations a Top and b bottom beam, c For non-uniform mesh

Table 3
Crack and FPZ propagation in simulations