Finite deformation implementation of a mixed-mode single-integral type cohesive zone with reorienting surfaces of separation

Abstract

To model material ductile failure and crack propagation, cohesive zone elements can be embedded along potential fracture paths in a finite element simulation. When damage criteria are met, elements in the mesh decohere, simulating the formation and propagation of a crack. In this paper, we present a novel computational algorithm based on finite deformation theory, essential to modeling crack initiation and growth in solids undergoing large deformations. This new algorithm was formulated within a Lagrangian frame of reference to extend previous cohesive zone algorithms to include modeling crack growth in finite deformation contexts. The local coordinate system, necessary for defining an embedded cohesive zone, is constructed based upon the current configuration and is updated within the nonlinear iteration process, thereby resulting in the convergence of the solution for a growing crack in a large deformation quasi-static setting. The model’s accuracy was demonstrated by comparing finite element model simulation results with the analytic case of a constant surface separation, as shown in the verification examples. The power and efficacy of the algorithm to capture large deformations during crack growth were then demonstrated with a double cantilever beam example case. It indicates that the model can be applied to a variety of physical circumstances for predicting crack initiation and growth with delamination and fracture.

1. Introduction

Cohesive zone models have been developed and used to model computational fracture mechanics, with many variations in the approach proposed to capture the hardening and softening phases for both brittle fracture and ductile crack growth. However, the concept of cohesive forces to model cracks predates numerical algorithms for cohesive zones [1,2] and forms the foundation of modern cohesive zone models for brittle and ductile fracture. In particular, for ductile crack growth, cohesive zone models have demonstrated considerable success in predicting both loading history and rate dependence of crack growth. These complexities cannot normally be modeled using conventional elastic fracture models [3].

To address experimentally observed complexities, the cohesive zone model employed in this work has a convolution integral [4] relating the crack face tractions to the history of the crack opening displacement, which inherently includes both loading rate and history dependence within the model formulation. This nonlinear cohesive zone model can predict ductile crack growth and has been constructed in an incremental form suitable for time-stepping finite element algorithms [4,5].

Previously, the cohesive zone model [4] did not include the complexities necessary for handling solids undergoing finite deformations. Alternatively, several other approaches have been proposed to account for crack growth in the presence of large deformations. These include phase-field modeling based upon the energy release due to fracture [6,7], XFEM for ductile crack growth [8], J-integral and dynamic M-integral approaches applied in an arbitrary Lagrangian Eulerian (ALE) frame of reference for quasi-brittle fracture [9], cohesive zone models within a Lagrangian frame of reference with the use of Nanson’s formula [10], finite band methods [11], peridynamics [12], element failure [13], and GraFEA [14].

There are also ongoing efforts to use hybrid approaches to obviate limitations of cohesive zone models, such as the crack location being known a priori and the difficulty in characterizing model parameters. Hybrid approaches include combining methods such as phase-field with cohesive zones [6,15–17] and XFEM with cohesive zones [18–20]. Characterizing model parameters is seeing headway by incorporating machine learning algorithms [21–23]. In addition, questions exist for multi-field problems and how a cohesive zone can represent reality in these cases. These challenges include the coupling mechanics with thermodynamics [15], chemistry [15], and electromagnetics [15,24,25].

Other open questions in cohesive zone modeling include the proper treatment of dynamic and history-dependent effects. The rate dependence exhibited by materials can be cumbersome, and the cohesive zone model used in this work is only one of several approaches proposed to capture this behavior [26,27]. Similar challenges arise for history-dependent behavior [28,29], particularly in the context of fatigue failure, which is a critical concern for aging structural design and must be addressed to maintain safety over the life span of a structure’s service life.

This study extends the above-mentioned cohesive zone model by Allen and Searcy [4] for application to include the prediction of crack growth in solids undergoing finite deformations. Several cohesive zone models have been previously proposed, with at least one focused on the softening branch of the model [30]. The issue to be addressed herein is the reorientation of the temporally extending surfaces exposed to the traction–separation relationship. Unlike prior work that focused on finding the normal direction to the crack face via Nanson’s Formula [10], the approach proposed herein focused on the physical orientation of the crack faces within a total Lagrangian frame of reference [31]. Consequently, the temporally changing orientation of a generic cohesive zone is defined within the algorithm, and the global displacement jump vector is mapped into a rotating frame of reference that is consistent with the current configuration of the material. By continuously using this reorientating procedure (as opposed to using Nanson’s formula), both the normal and tangential directions within a generic cohesive zone are simultaneously defined for the purpose of modeling different directions of crack propagation based on the applied stress, termed mode 1 (opening), mode 2 (in-plane sliding), and mixed-mode (both opening and sliding simultaneously) fracture.

Accurately simulating these fracture modes is essential for predicting crack propagation in complex materials such as biologic tissues. The importance of modeling such phenomena cannot be overstated, as it directly relates to the health of living things. The proposed extension to the cohesive zone model herein shows its validity in large deformation contexts. This validity, in conjunction with the single-integral nature of this model, allows for capturing the behavior of materials separating while handling the complexities associated with materials that are governed by multiple time scales, something that the proposed model lends itself to with viscoelastic-based kernel functions. This capability is one that other modern cohesive zone models lack. The continuous consideration of multiple time scales of a material via a single-integral formulation makes this work uniquely poised to address the separation of highly viscoelastic materials such as those found within biologic entities.

2. Description of the initial boundary value problem

The first step in posing the initial boundary value problem (IBVP) is to define the spatial domain, including both the interior and external surfaces, as illustrated in Figure 1. Within this context, consider a two-dimensional domain with volume ( $\mathrm{\Omega }({\mathrm{x}}_{\mathrm{i}},\mathrm{t})$ ) with one-dimensional curvilinear surfaces for the external boundary (Γ^e(x_i,t)) and internal boundaries (Γⁱ(x_i,t)). We note, however, that the same procedure can be extended to three dimensions.

Figure 1.

Problem domain.

Next, capturing crack surfaces that temporally vary requires separating the external boundary, denoted by superscript e, into up to two subsections, denoted by subscripts 1 and 2. Subsection 1 denotes a surface with traction boundary conditions and cohesive zone tractions ( $t_i\left({\mathrm{\Gamma }}_1^e\right)$ ), whereas subsection 2 denotes a surface with displacement boundary conditions ( $u_i\left({\mathrm{\Gamma }}_2^e\right)$ ). Following this notation, the boundary description (equation (1)) is:

$${\mathrm{\Gamma }}_1^e\cap {\mathrm{\Gamma }}_2^e={\mathrm{\Gamma }}^e,\mathrm{where}{\mathrm{\Gamma }}_1^e\cap {\mathrm{\Gamma }}_2^e=0$$

(1)

Furthermore, note that there can be any arbitrary number of interior boundaries, denoted by superscript i in a problem that will be simulated using the new cohesive zone model. Including a subscript to denote each potential crack, the total of all interior boundaries is described by equation (2).

$$\sum _{k=1}^n{\mathrm{\Gamma }}_k^i={\mathrm{\Gamma }}^i$$

(2)

Next, the equations that govern the IBVP, as defined above, are stated below here. The Lagrangian form of the conservation of linear momentum (equation (3)) is solved under the quasi-static assumption, where the Green–Lagrange strain (equation (4)) is employed. For demonstrative purposes, a standard isotropic linear elastic constitutive model (equation (5)) is used to describe the bulk material behavior, although a more intricate constitutive model can be used, depending on the constitution of the bulk material behavior that may arise in application. A previously proposed single-integral cohesive zone model [4] (equations (6) and (7)) and an accompanying damage evolution law (equation (8)) are used herein for the internal boundaries that are in contact with one another. The Lagrangian frame of reference makes it such that the partial derivative of the opening parameter is equivalent to the total derivative, and thus, they are interchangeable. The externally applied loading may result in the separation of some or all of the predefined surfaces. The nonlinear nature of the chosen form of conservation of linear momentum and the chosen strain measure, both introduce complexity to the physical and mathematical nature of the problem. It should also be noted that the deployment of a cohesive zone model requires that the potential crack path be known a priori.

$$\nabla (S^T·F^T)+\rho f=0\mathrm{in}\mathrm{\Omega }$$

(3)

$$\epsilon =\frac{1}{2}(\nabla u+{(\nabla u)}^T+\nabla u·{(\nabla u)}^T)\mathrm{in}\mathrm{\Omega }$$

(4)

$$S=C\epsilon \mathrm{in}\mathrm{\Omega }$$

(5)

$$t_i^c=\frac{1}{\lambda }\frac{u_i^c}{u_i^{*}}(1-\alpha \left(t\right)){\int }_0^{t}E(t-\tau )\frac{\partial \lambda }{\partial \tau }d\tau \mathrm{in}{\mathrm{\Gamma }}_c^{\mathrm{e}}\mathrm{and}{\mathrm{\Gamma }}^i$$

(6)

$$\lambda ={[{\left(\frac{u_n}{u_n^{*}}\right)}^2+{\left(\frac{u_r}{u_r^{*}}\right)}^2]}^{\frac{1}{2}},\mathrm{where}{u}_n\ge 0\mathrm{in}{\mathrm{\Gamma }}_c^{\mathrm{e}}\mathrm{and}{\mathrm{\Gamma }}^i$$

(7)

$$\frac{d\alpha }{\mathrm{dt}}={\alpha }_1{\left(\lambda \right)}^n\mathrm{when}\frac{\partial \lambda }{\partial t}>0\mathrm{and}\alpha <1;\frac{d\alpha }{\mathrm{dt}}=0\mathrm{when}\frac{\partial \lambda }{\partial t}\le 0\mathrm{or}\alpha =1\mathrm{in}{\mathrm{\Gamma }}_c^{\mathrm{e}}\mathrm{and}{\mathrm{\Gamma }}^i$$

(8)

In the above equations, the variables are the second Piola–Kirchoff stress tensor (S), the deformation gradient (F), density (ρ), body force (f), Green–Lagrange strain (ε), global displacement (u), linear elastic stiffness matrix (C), cohesive zone traction (t ^c ), cohesive zone length parameter (u^*), relaxation modulus (E(t)), non-dimensional cohesive zone opening (λ), cohesive zone normal displacement (u _n ), cohesive zone tangent displacement (u _r ), damage function (α), and the damage evolution parameters (α ₁ and n).

The highly nonlinear nature of the separation of the surfaces governed by the single-integral cohesive zone model is the primary focus of this research. As such, the cohesive zone model naturally accounts for crack formation and growth within the time-stepping algorithm using an iteration procedure. The cohesive zone model (equations (6)–(8)) predicts the statically equivalent interfacial tractions to hold the faces together within the boundary subdomains modeled by the cohesive zone model until the damage function (α) reaches a value of 1. At this point, the tractions holding the faces together become zero and the faces are no longer subject to cohesive loading, thereby converting the internal boundary to traction-free surfaces. Successively increasing external loading will then induce crack growth wherever the damage reaches the limiting value of 1 within the domain of interest. During the damage process, energy is dissipated primarily through the damage parameter that has a micromechanical basis in fibrillation, but the algorithmic implementation is phenomenological, as represented by equation (8). The energy from the accumulation of damage is unrecoverable, as is the case for most real materials. In addition, energy is dissipated through viscosity effects that come from the relaxation modulus in equation (6).

The proposed algorithm is capable of accounting for crack growth in bodies experiencing large deformations. Toward this end, the required inputs are the domain geometry ( $\mathrm{\Omega }({\mathrm{x}}_{\mathrm{i}},t),{\mathrm{\Gamma }}^e({\mathrm{x}}_{\mathrm{i}},t),\mathrm{and}{\mathrm{\Gamma }}^i({\mathrm{x}}_{\mathrm{i}},t))$ , the material properties for the bulk medium (C), the cohesive zone material properties (E(t), $u_i^{*}$ , ${\alpha }_1$ , and n), and the boundary conditions (u _i and t _i on Γ^e). This concludes the description of the IBVP.

3. Cohesive zone reorientation methodology

For the purpose of tracking the reorientation of the internal crack surfaces, the weak form of conservation of linear momentum (equation (3)) is solved herein. Within each input load increment, displacement increments are predicted based upon the initial configuration and the total displacement accumulated from all previous load steps [31]. For improved accuracy, nine-node quadrilateral elements are used. With these elements, the position and displacement fields are represented with bi-linear quadratic Lagrange interpolation polynomials and the consistency of the cohesive zones is enforced via nodal area weighting. With the quadratic shape functions, the result is that one-sixth of the surface area is given to each corner node and two-thirds is given to the middle node. Direct iteration is used to solve the nonlinear system of equations. This concludes the numerical procedure used herein for solving the IBVP.

A necessary extension of a previous algorithm [4] is to account for the reorientation of the cohesive zones throughout the time-stepping procedure. A numerical scheme for accomplishing this is, therefore, proposed herein to account for this phenomenon in each cohesive zone within the spatial domain. A local coordinate system is defined in the reference configuration in terms of the normal and tangential directions within each of the cohesive zones, as shown in Figure 2.

Figure 2.

Reference configuration of cohesive zones between elements.

Therefore, a mapping is necessary as the edge reorients during each loading step. These normal and tangent directions are necessary to solve equations (6) through (8). An isoparametric-based mapping of this reorientation was used to calculate the local coordinate directions as a function of the global configuration of the elements. By considering only the edges connected by a cohesive zone (red circles shown in Figure 2), an arbitrarily reoriented global configuration is mapped into the local frame based on the element configuration (shown in Figure 3).

Figure 3.

Reoriented cohesive zones between elements.

The current configuration within a generic cohesive zone is the sum of the reference configuration and the total displacement vector therein, which results in the current configuration of a generic cohesive zone within the spatial domain. To account for this, the notation used herein is as follows: subscripts 1 and 2 correspond to the separate nodes within a single cohesive zone, resulting in a global opening displacement (equation (9)) within the cohesive zone denoted by a subscript c, thereby resulting in the following:

$$\overrightarrow{u_c}=(\overrightarrow{x_2}+\overrightarrow{u_2})-(\overrightarrow{x_1}+\overrightarrow{u_1})$$

(9)

This global displacement must then be mapped into the reoriented local basis directions of the cohesive zone. This is achieved through an isoparametric mapping, wherein the Jacobian defined in equations (10) and (11) is used to find the global orientation of the element edges.

$$J=\left[\begin{array}{cc}\frac{\partial x_1}{\partial \xi }\hfill & \frac{\partial x_2}{\partial \xi }\hfill \\ \multicolumn{1}{c}{\frac{\partial x_1}{\partial \eta }}& \frac{\partial x_2}{\partial \eta }\hfill \end{array}\right]$$

(10)

$$\left[\begin{array}{c}{\mathrm{dx}}_1\hfill \\ \multicolumn{1}{c}{{\mathrm{dx}}_2}\end{array}\right]=J^T\left[\begin{array}{c}d\xi \hfill \\ \multicolumn{1}{c}{d\eta }\end{array}\right]$$

(11)

The orientation of one edge is used to define the cohesive zone orientation. Note that in the event of a corner node, the orientation of neighboring elements is averaged. Furthermore, the latter statement converges with mesh refinement as the orientation of neighboring elements approaches equality, wherein the edge length of the elements approaches zero.

The orientation equation varies slightly, depending on which isoparametric variable (ξ,η) is deployed. The following mapping is for a ξ = constant edge. The process begins by defining the angle of the cohesive zone in the global frame (see Figure 4).

Figure 4.

Global angle.

This angle (equation (12)) can then be calculated using equations (10) and (11), noting that dξ is zero along the edge. The angle can then be expressed in terms of Jacobian indices.

$$\theta =\arctan \left(\frac{{\mathrm{dx}}_1}{{\mathrm{dx}}_2}\right)=\arctan \left(\frac{J_{21}}{J_{22}}\right)$$

(12)

A similar result is obtained when considering an edge where η = constant, but in terms of J ₁₁ and J ₁₂ . This orientation requires a correction depending on the signs of the Jacobian to have a correct global orientation for an arbitrary configuration. Note that a phase shift is employed to ensure correct global orientation for the displacement mapping when Jacobian indices are negative.

Once the orientation is determined, the global displacement vector defined in equation (9) can be mapped into the (n,r) local coordinate system (equation (13)). The rotation matrix (R) is defined by which of the four edges (ξ = −1, ξ = 1, η = −1, or η = 1) is currently considered, but note that all the edges follow the following equation:

$$\left[\begin{array}{c}{u}_n\hfill \\ \multicolumn{1}{c}{u_r}\end{array}\right]=R\left[\begin{array}{c}{u}_1\hfill \\ \multicolumn{1}{c}{u_2}\end{array}\right]$$

(13)

The incrementalized form of the cohesive zone tractions can now be calculated in terms of the reoriented cohesive zone surface. Finally, note that the above process must be carried out for both the case of the total displacement and the displacement increment during the iterative solving process to reach proper convergence. This is due to the dependence in the incrementalized form [4] of equation (6) on both the displacement at the beginning of the step and during the step.

4. Verification of the large deformation algorithm

The verification case is a prescribed constant rate of opening ( $\frac{\partial u_n}{\partial t}$ ), including a generalized standard linear solid material model deployed within the cohesive zone. This test case has been modeled in three versions to test the functionality of the algorithm: single face separation, multi-body separation, and large rotation separation, respectively. The first case is intended to verify that the cohesive zone itself is performing correctly, the second is to ensure that multi-body behavior is accounted for via multiple cohesive zones at a single point, and the third is to show that the reorientation algorithm is working properly. The algorithm’s output is compared to the analytic solution [4] written as follows:

$$t_n\left(t\right)=\frac{\frac{\partial u_n}{\partial t}t}{u_n^{*}\lambda }[1-\frac{{\alpha }_1{\left(\frac{\partial u_n}{\partial t}\right)}^n}{{\left(u_n^{*}\right)}^n(n+1)}{t}^{n+1}][E_{\infty }t+E_1{\rho }_1(1-e^{\frac{-t}{{\rho }_1}})]\left(\frac{\frac{\partial u_n}{\partial t}}{u_n^{*}}\right)$$

(14)

It is important to note how a lambda equal to zero is handled in the algorithm (i.e., when there is no opening). Due to the physical nature of the problem, it is known that if a crack is under no loading, no traction will be generated within the cohesive zone. Therefore, when there is no opening, traction is set to zero within the algorithm.

Separation of a single face is redundant when viewing the results of the multi-body separation and the large rotation separation cases and is thus omitted from the following figures. In addition, due to the nature of the analytic solution, the tangential opening case is mathematically identical to the normal opening cases shown and is therefore omitted from the figures for brevity. For both verification cases shown herein, the configuration at three distinct times will be elucidated: the initial configuration, a timestep near the middle of the simulation, and the final timestep, whereby full separation has occurred (α = 1).

The analytic solution is identical to that predicted by the algorithm proposed herein for both the multi-body separation and rotating-body separation cases. Note that the color of the line representing the cohesive zone in the length scale of separation images is defined by the value of the damage parameter at that point in time. The yellow in Figure 5 denotes that 1/3 < α < 2/3, and the red in Figure 6 denotes that 2/3 ≤ α ≤ 1. The above results are in complete agreement with the anticipated results, thereby demonstrating the accuracy of the algorithm proposed herein.

Figure 5.

Multi-body separation verification case.

Figure 6.

Full rotation (360 degrees) verification case.

Figure 6 illustrates the advantages of the finite deformation algorithm developed in this work. In contrast, a small-deformation algorithm cannot adequately account for the changes in orientations induced by large rotations. The finite deformation mapping scheme is a requirement to capture arbitrarily reorienting cohesive zone surfaces.

5. Modeling demonstration via example cases

5.1. Uniaxial separation with necking

The numerical example in this section consists of a uniaxial bar subjected to large deformations that induce necking in the simulation. Cohesive zone elements are inserted along the longitudinal midline of the bar (vertical red line in Figure 7). The left end of the bar is fixed, and the right end is subjected to a constant rate of displacement. Two cases were run, one with Poisson’s ratio equal to 0 to inhibit necking and one with Poisson’s ratio equal to 0.499 to promote necking. The purpose of this example is to demonstrate the importance of using the finite strain formulation for problems involving large deformations. In particular, it illustrates the effect of cross-sectional area reduction that comes from necking, which accelerates the separation compared to a case without area deduction. Because necking cannot be accurately captured by small-deformation formulations, the earlier failure of the bar underscores the importance of the large deformation formulation used in this work.

Figure 7.

Initial configuration of the uniaxial bar case.

The simulation was run with consideration of the symmetry along the vertical midline of the bar, meaning the results will only show the top half of the bar. Only the necking deformed configuration is shown in Figure 8, as the non-necking case is trivial and has no area reduction.

Figure 8.

Deformed configuration of necking case.

The importance of the necking case is twofold. First, the large deformations lead to a non-homogeneous strain field, inducing spatially varying damage rates in the cohesive zones (the red line in Figure 8), with the bar failing first at the top surface. Second, the reduction in area caused by necking requires the cohesive zone to produce larger tractions to maintain equilibrium, which in turn causes the surface to separate more. Together, these factors lead to a larger damage rate and substantially earlier failure. Figure 9 shows the damage progression of the system with and without necking; the necking case shown is for the top cohesive zone, where failure first occurs.

Figure 9.

Damage progression with and without necking.

With necking, the system fails at less than half the right-end displacement of the bar compared to the case without necking. This shows the importance of the large deformation formulation for predicting ductile failure of a material, as small-deformation formulations cannot model necking and therefore severely overpredict a structure’s ability to withstand loading.

5.2. Double cantilever beam delamination

The numerical example in this section consists of a double cantilevered beam with cohesive zone elements inserted along the horizontal midline (blue line in Figure 10) with an initial separation (horizontal red line in Figure 10) at the left end. The top and bottom of the free end of the beam are loaded vertically with forces of equal magnitude and opposite direction to induce separation along the cohesive zones, as shown in Figure 10. Cohesive zones are embedded along the unseparated part of the midline, and the load is increased until a significant reorientation has occurred, meaning the tangential component of the cohesive zone opening is no longer negligible. The initial goal of the simulation is to produce a deformed orientation of the free end of the beam that is greater than 15 degrees.

Figure 10.

Initial configuration of the double cantilever beam case.

The goal of this example is to demonstrate the capability of the cohesive zone reorientation algorithm, together with the total Lagrangian finite deformation solving algorithm. The following results use a fourth-order polynomial curve fit along the top surface of the symmetric separation. This polynomial has been chosen for comparison based upon the order of displacement found via the Euler–Bernoulli beam theory and has resulted in a perfect correlation (R² = 1) when applied to the separation surfaces. A purple line along the top surface of separation presents this, and both the equation and correlation coefficient are included.

Figure 11 shows snapshots of damage progression within the cohesive zone at load step intervals of 100. These snapshots represent the progression of the damage evolution law (α(t) shown in equation (8)) throughout loading, where it should be noted that damage equal to unity implies crack separation.

Figure 11.

Damage progression within the cohesive zones at load step intervals of 100.

This monotonic increase in the damage evolution law, with the opening end experiencing the largest damage before the onset of crack growth, is consistent with expected physical behavior. Next, a sequence of damaged configurations is shown in Figure 12. It shows the progression of deformation, the trendline and correlation value associated with the separation surface, and a comparison of the reorientation algorithm (blue line) with the polynomial curve fit (orange dotted line). These configurations demonstrate the capability and accuracy of the reorientation algorithm developed in this study.

Figure 12.

(a) Results after load step 1, (b) results after load step 100, (c) results after load step 200, (d) results after load step 300, and (e) results after load step 400.

As shown in Figure 12, the progression of the interpolation is a consistent representation of the deformed configuration, which is confirmed by a strong agreement between the slope of the curve fit and the representation of the physical space predicted by the finite element model. Furthermore, the reorientation algorithm agrees with the interpolation.

6. Conclusion

In this work, a new cohesive zone model was proposed that works with time-stepping finite element methods to predict crack initiation and propagation in solids experiencing finite deformations. This study developed a new approach focused on the physical orientation of crack faces within a total Lagrangian framework. This method defines the changing orientation of cohesive zones and maps displacement vectors into a rotating frame consistent with the material’s current configuration, allowing for the simultaneous definition of normal and tangential directions. Using this approach, the cohesive zone model can accurately capture various fracturing modes (i.e., the opening mode, sliding mode, and mixed modes) in large deformation applications.

This new cohesive zone model was verified using analytical cases involving constant surface separation. The verification cases showed that the model accurately predicts the event of multiple bodies separating and that the accuracy holds for arbitrary reorientation of crack surfaces. A further demonstration of the modeling with an example case (i.e., a double cantilevered beam largely deformed and delaminated) confirms the capability and the accuracy of the cohesive zone model with reorienting surfaces of separation. The numerical example also highlights the algorithm’s effectiveness in predicting crack initiation and growth, illustrating its strengths for use across physical scenarios.