Numerical study of delamination process of the CFRP composite

In recent years, progress in the analysis and development of composite materials has greatly impacted their application and use across various fields of engineering. As a result, scientific research on multi-domain contact problems holds significant promise for advancing applications in civil and mechanical engineering, as well as in the aeronautical and aerospace industries. This study focuses on analyzing delamination in interface contact issues within composite materials, which is also broadly relevant to engineering practice. A key challenge in designing and developing these structures is understanding the mechanisms of onset and progression of cracks $[1,2,3,4]$$[1,2,3,4]$[1,2,3,4][1,2,3,4]$[1,2,3,4]$.
Currently, the use of Cohesive Zone Models (CZM) to address fracture issues at interfaces is on the rise. These models are distinguished by a non-singular stress field near the crack tip, making CZM an effective method for simulating both the onset and progression of cracks [5,6,7].
In Damage Mechanics, interface damage is characterized by a non-dimensional damage parameter that reflects the irreversible process of decohesion. This concept has been utilized in interface contact models, where the degradation of adhesion, known as interface damage, is analyzed and discussed. The reduction in adhesion is governed by the irreversible nature of interface damage and is influenced by both normal and tangential deformations. Crack growth along the interface is driven by the energy release rate (ERR), which controls the advancement of crack progression [8,9,10].
To more accurately capture the physical behavior of the developed approach, the model incorporates viscoelastic properties [11]. The proposed method disregards inertial effects in its formulation, due to the slow nature of the external loading, meaning the model is treated as quasistatic. The theory of viscoelasticity accounts for the time-dependent relationship between stress and strain, allowing the model to be considered rate-dependent, while still being quasistatic. Quasistatic linear viscoelasticity theory offers a practical and reliable engineering approximation for numerous physical and mechanical applications involving composite materials.
Recently, the study and numerical analysis of interface degradation using computational contact mechanics have become critical for predicting failure loads in fiber-reinforced composites. One of the most effective methods for analyzing and simulating interface damage, particularly crack initiation and propagation, involves the use of CZM combined with frictional contact $[12,13]$$[12,13]$[12,13][12,13]$[12,13]$. However, numerically analyzing such nonlinear phenomena presents a complex and challenging problem that is not yet fully understood, with multiple approaches still being developed by researchers.
De Xie and Waas [14] have conducted similar research, utilizing a discrete cohesive zone model (DCZM) with the finite element method (FEM) to simulate fracture initiation and growth, particularly when material nonlinearities are significant. In their study, they addressed the issue of mesh size dependency by applying nodal forces to the rod elements, focused on a 2D analysis, and did not account for the viscosity parameter. A large number of numerical simulations of delamination using FEM are presented in [15,16,17]. Results from a numerical simulation of a double cantilever beam (DCB) under ModeI loading are provided in [18] to validate the modified adaptive cohesive model.
Research on fracture mechanics test methods aimed at determining delamination resistance or fracture toughness in fiberreinforced polymer-matrix composites remains a key area of investigation [19,20].
This study concentrates on the numerical analysis of delamination in interface contact problems within CFRP composite materials using the FEM and compares the numerical results for critical force with theoretical predictions.

In the plane of delamination, delaminated structure is divided into two sublaminates of thickness $h_M,h_N$$h_M,h_N$h_(M),h_(N)h_{M}, h_{N}$h_M,h_N$ (Fig.1), which are considered to be in the in-plane dimensions. The sublaminates are modeled as an assembly of in-plane elements connected by zero-thickness layer. The bond between both sublaminates is maintained by applying constraint equations, which are enforced using Lagrangian multipliers [10].

The ERR controls the progression of crack growth at the interface based on the given relationship

where $\zeta$$\zeta$zeta\zeta$\zeta$ expresses the damage parameter, which varies between one and zero at each interface point.

signifies that the adhesive is completely damaged at a specific point. i.e. complete delamination.

indicates an undamaged interface.
The total ERR in DCB sample can be formulated as

where $dU$$dU$dUd U$dU$ is stress increment, $da$$da$dad a$da$ is length increment (or crack increment) and $b$$b$bb$b$ is width of sample. With respect to the examination of the DCB sample using Mode I, the value of the ERR can be expressed using the Euler-Bernoulli theory as

where $EI$$EI$EIE I$EI$ is bending stiffness.
According to ASTM standards [21], there are several ways to determine the value of ERR, one of the methods is the socalled Modified Beam Theory (MBT). The $G_I$$G_I$G_(I)G_{I}$G_I$ value according to this theory is determined as

Due to the fact that the experiment is partially imperfect due to theoretical studies, a correction is introduced, which assumes that the delamination value is increased by the value $|\mathrm{\Delta }|$$|\mathrm{\Delta }|$|Delta||\Delta|$|\mathrm{\Delta }|$.
Thus, the total crack length is ( $\mathrm{a}+|\mathrm{\Delta }|$$\mathrm{a}+|\mathrm{\Delta }|$a+|Delta|\mathrm{a}+|\Delta|$\mathrm{a}+|\mathrm{\Delta }|$ ) and so the adjusted $G_I$$G_I$G_(I)G_{I}$G_I$ calculation looks as follows

where $\mathrm{\Delta }$$\mathrm{\Delta }$Delta\Delta$\mathrm{\Delta }$ is calculated in accordance with the ASTM standards.

If the laminate consists of $S$$S$SS$S$ sublaminates, each containing $n_q$$n_q$n_(q)n_{q}$n_q$ plies, where $q=1,\dots ,S$$q=1,\dots ,S$q=1,dots,Sq=1, \ldots, S$q=1,\dots ,S$, then components of in-plane stiffness matrix $\mathbf{A}$$\mathbf{A}$A\mathbf{A}$\mathbf{A}$ can be expressed as the sum of $S$$S$SS$S$ terms [22]

where $h_{k{n}_q}$$h_{k{n}_q}$h_(kn_(q))h_{k n_{q}}$h_{k{n}_q}$ is thickness and $E_{ij\left(k{n}_q\right)}$$E_{ij\left(k{n}_q\right)}$E_(ij(kn_(q)))E_{i j\left(k n_{q}\right)}$E_{ij\left(k{n}_q\right)}$ is component of transformed reduced stiffness matrix of the $k^{\text{th}}$$k^{\text{th }}$k^("th ")k^{\text {th }}$k^{\text{th }}$ layer of the $q^{\text{th}}$$q^{\text{th }}$q^("th ")q^{\text {th }}$q^{\text{th }}$ sublaminate, $\mathrm{and}{\left(A_{ij}\right)}_q$$\mathrm{and}{\left(A_{ij}\right)}_q$and(A_(ij))_(q)\operatorname{and}\left(A_{i j}\right)_{q}$\mathrm{and}{\left(A_{ij}\right)}_q$ is the component of the in-plane stiffness matrix of the $q^{\text{th}}$$q^{\text{th }}$q^("th ")q^{\text {th }}$q^{\text{th }}$ sublaminate.
For the two sublaminates $M$$M$MM$M$ and $N$$N$NN$N$ the combined in-plane stiffness matrix $\mathbf{A}$$\mathbf{A}$A\mathbf{A}$\mathbf{A}$ is

The components of bending-extension coupling stiffness matrix $\mathbf{B}$$\mathbf{B}$B\mathbf{B}$\mathbf{B}$ of the sublaminate can be written as

where ${\overline{z}}_k$${\overline{z}}_k$bar(z)_(k)\bar{z}_{k}${\overline{z}}_k$ is the $z$$z$zz$z$-coordinate of the mid-plane of the $k^{\text{th}}$$k^{\text{th }}$k^("th ")k^{\text {th }}$k^{\text{th }}$ layer assessed from the mid-plane of the combined laminate.

For the two sublaminates $M$$M$MM$M$ and $N$$N$NN$N$ the bending-extension coupling stiffness matrix $\mathbf{B}$$\mathbf{B}$B\mathbf{B}$\mathbf{B}$ can be written as

Similarly, the bending stiffness matrix $\mathbf{D}$$\mathbf{D}$D\mathbf{D}$\mathbf{D}$ for the combined laminate can be expressed as

Definition of visco-elasticity theory assess in the formulation dependence of stress and strain on time, such the model can be considered as the rate dependent (but quasistatic). The efficient methodology for FEM analysis of quasistatic visco-elastic delamination is considered in this research.
The behavior of polymers in the field of linear viscoelasticity can be studied using relaxation methods and creep tests [23]. In the first case, the result is a time dependence of the stress at a constant deformation, in the second case the time dependence of the relative deformation at constant stress. To describe this behavior, it is possible to use various rheological models, which are usually created by two basic elements: a Hooke's spring (H) with an elastic modulus $E$$E$EE$E$ and a viscous damper (N) with a liquid of viscosity $\eta$$\eta$eta\eta$\eta$. In the case of Maxwell model (Fig. 2), these elements are connected in series. Both the spring and damper transmit the same stress, while the overall deformation of the model in the loaded state is given by the sum of the deformations of both members. While the tension in the spring is controlled By Hooke's law, the stress in the damper is proportional to the rate of deformation. Maxwell model is suitable for stress relaxation description. The Kelvin-Voigt model (Fig. 2) can be used to describe creep, for which parallel connection of spring and damper is typical.

To establish the global constitutive relation of the Kelvin-Voight model we are first supposed to collect the corresponding geometric formulas following from the parallel connection and the physical relation $\sigma =E\epsilon$$\sigma =E\epsilon$sigma=E epsi\sigma=E \varepsilon$\sigma =E\epsilon$ and $\sigma =\eta \dot{\epsilon }$$\sigma =\eta \dot{\epsilon }$sigma=etaepsi^(˙)\sigma=\eta \dot{\varepsilon}$\sigma =\eta \dot{\epsilon }$, where $\eta$$\eta$eta\eta$\eta$ is dynamical viscosity.
Constitutive equation of the Kelvin-Voight model

The global physical relation of the Maxwell model is

The tangent of the phase angle is one of the basic measurable properties of the material and refers directly to its damping. Basically, it is a quantitative indicator of how effectively the material dissipates energy due to molecular displacements and internal friction. It is defined as a ratio loss and elastic modulus and thus is independent of the sample geometry

where $E^{\prime }$$E^{\prime }$E^(')E^{\prime}$E^{\prime }$ represents the real and $E^{\prime \prime }$$E^{\prime \prime }$E^('')E^{\prime \prime}$E^{\prime \prime }$ the loss (imaginary) component of complex modulus

The main idea of FEM is the discretization of the continuous, investigated area. Discretization is defined by creating a mesh that divides the continuous area into subareas called finite elements. According to the Ritz method applied to elastostatic problems, under the functional $\mathrm{\Pi }$$\mathrm{\Pi }$Pi\Pi$\mathrm{\Pi }$, we will consider potential energy.
In connection with the Ritz method, we approximate the field of displacements $\tilde{\mathbf{u}}$$\tilde{\mathbf{u}}$tilde(u)\tilde{\mathbf{u}}$\tilde{\mathbf{u}}$

where $\mathbf{N}$$\mathbf{N}$N\mathbf{N}$\mathbf{N}$ is the matrix of the shape functions, $\mathit{x}=(x,y,z)$$\mathit{x}=(x,y,z)$x=(x,y,z)\boldsymbol{x}=(x, y, z)$\mathit{x}=(x,y,z)$ is the position vector and $\mathbf{u}$$\mathbf{u}$u\mathbf{u}$\mathbf{u}$ is the finite element displacement vector.
We express the vector of stresses and strains as follows

where $\mathbf{E}$$\mathbf{E}$E\mathbf{E}$\mathbf{E}$ is the matrix of elastic constants, which we obtain by transforming from local to global coordinates using the transformation matrix $\mathbf{T}$$\mathbf{T}$T\mathbf{T}$\mathbf{T}$.

We express the total potential energy variation by the relation

of which applies

where $\mathbf{K}$$\mathbf{K}$K\mathbf{K}$\mathbf{K}$ is the stiffness matrix

${\mathbf{f}}_p$${\mathbf{f}}_p$f_(p)\mathbf{f}_{p}${\mathbf{f}}_p$ and ${\mathbf{f}}_q$${\mathbf{f}}_q$f_(q)\mathbf{f}_{q}${\mathbf{f}}_q$ are vectors of volume and surface external forces

If the components $\delta \mathbf{u}$$\delta \mathbf{u}$deltau\delta \mathbf{u}$\delta \mathbf{u}$ are independent, we get from equation (24) the linear equation system

The given equations valid for the finite element have an added index $E$$E$EE$E$. We express the internal energy of the finite element

where ${\mathbf{K}}_E$${\mathbf{K}}_E$K_(E)\mathbf{K}_{E}${\mathbf{K}}_E$ is the finite element stiffness matrix

Using the Boolean matrix ${\mathbf{L}}_E$${\mathbf{L}}_E$L_(E)\mathbf{L}_{E}${\mathbf{L}}_E$ is determined the correct position of each finite element in the stiffness matrix of the entire structure. For the entire structure, the vector of displacements of the finite element is placed into the vector of displacements using the relation

Then we obtain a system of equations by summation over all elements

After applying the boundary conditions, the matrix $\mathbf{K}$$\mathbf{K}$K\mathbf{K}$\mathbf{K}$ becomes a regular, i.e. positive definite matrix. We used the programs ADINA [24] for the numerical analysis of delamination process of the CFRP composite.

The mode I delamination of DCB within the FEM interface damage model is considered (Fig. 3). The DCB made of a laminate composite CFRP material is subjected to displacement $u$$u$uu$u$ at its unsupported ends. The laminate beam consists of two sublaminates. Each sublaminate consists of 8 UD laminate plies. Prescribed displacements are used as displacement-controlled loading for 100-time steps of size 0.01 s and full displacement of 7.15 mm at time 1.0. The automatic time-stepping method is used. Contact elements by means of springs are used in the model (Fig. 4). The mechanical viscoelastic characteristics of each layer, obtained from Impulse Excitation Technique (IET) are shown in Fig. 5 [11].
Cohesive properties:
$G_{\mathrm{Ic}}=0.262\mathrm{N}/\mathrm{mm}$$G_{\mathrm{Ic}}=0.262\mathrm{N}/\mathrm{mm}$G_(Ic)=0.262N//mmG_{\mathrm{Ic}}=0.262 \mathrm{~N} / \mathrm{mm}$G_{\mathrm{Ic}}=0.262\mathrm{N}/\mathrm{mm}$ [25], IL cohesive tensile strength $=75\mathrm{N}/{\mathrm{mm}}^2$$=75\mathrm{N}/{\mathrm{mm}}^2$=75N//mm^(2)=75 \mathrm{~N} / \mathrm{mm}^{2}$=75\mathrm{N}/{\mathrm{mm}}^2$, penalty stiffness $={10}^6\mathrm{N}/{\mathrm{mm}}^3$$={10}^6\mathrm{N}/{\mathrm{mm}}^3$=10^(6)N//mm^(3)=10^{6} \mathrm{~N} / \mathrm{mm}^{3}$={10}^6\mathrm{N}/{\mathrm{mm}}^3$.
Dimensions of DCB:
$L=100\mathrm{mm},b=25\mathrm{mm},b=2\mathrm{mm},a=50\mathrm{mm}$$L=100\mathrm{mm},b=25\mathrm{mm},b=2\mathrm{mm},a=50\mathrm{mm}$L=100mm,b=25mm,b=2mm,a=50mmL=100 \mathrm{~mm}, b=25 \mathrm{~mm}, b=2 \mathrm{~mm}, a=50 \mathrm{~mm}$L=100\mathrm{mm},b=25\mathrm{mm},b=2\mathrm{mm},a=50\mathrm{mm}$. The clamped edge is long 10 mm at the right side of DCB.
The 2-D plane strain elements were used for the problem by defining the orthotropic material. For numerical solution the cohesive interfaces were defined and the low-speed dynamics were specified. Delamination growth analysis is performed iteratively in time steps. Contact elements provide connection of elements in nodes by means of springs, which allow displacements in in-plane directions with prevention of layer penetration. From experiments, the actual mechanical properties can be obtained by the Resonalyser test within IET. UD carbon fibre-reinforced beam of 16 autoclaved prepreg carbon/epoxy layers is tested at ${22}^{\circ }\mathrm{C}$${22}^{\circ }\mathrm{C}$22^(@)C22^{\circ} \mathrm{C}${22}^{\circ }\mathrm{C}$, that is interpreted in [11].

For comparison, the results of critical force obtained from ADINA [24] numerical simulation and experiment [25] are shown in Table 1. Numerical results in chosen time steps are written in Table 2.

Amulti-layer laminate model has been introduced, which includes interface model grounded in fracture and contact mechanics.
The cohesive interface model within the FEM analysis was used to obtain the critical force away to the crack tip and, thus, to compare the numerical value of critical force in mode I delamination with experimental result. Contour plot of $z$$z$zz$z$-displacement in time 0.25 s and 1 s is shown in Figs. 6 and 7, respectively. Cohesive normal stress, maximum normal stress ${\sigma }_y$${\sigma }_y$sigma_(y)\sigma_{y}${\sigma }_y$ and cohesive damage in depending on time are presented in Figs. 8, 9 and 13, respectively. Load-displacement diagram is presented in Fig. 10. We can see in the Fig. 10 the nonlinear behaviour between load $F$$F$FF$F$ and the transverse displacement $\delta$$\delta$delta\delta$\delta$ after obtaining the crack initiation by opening displacement of 3.58 mm . In this nonlinear region, the interface strengthening is observed to the opening displacement of 4.44 mm . Up to reaching this point, dependence in loaddisplacement diagram is decreasing. This post-critical region has a large effect on the total fracture energy value compared to the critical value. Relative displacement and load in depending on crack length are presented in Figs. 11 and 12, respectively.

The critical load obtained from numerical modeling is 73.96 N and from experiment [25] is 77.14 N . Maximum of load, obtained from load-displacement diagram is 80.87 N . The delamination was stopped, when the crack length $a=90\mathrm{mm}$$a=90\mathrm{mm}$a=90mma=90 \mathrm{~mm}$a=90\mathrm{mm}$ and $F=50.24\mathrm{N}$$F=50.24\mathrm{N}$F=50.24NF=50.24 \mathrm{~N}$F=50.24\mathrm{N}$. The numerical example illustrates the effectiveness of the proposed model.

This work was supported by Slovak Research and Development Agency of the Ministry of Education, Research, Development and Youth of the Slovak Republic under Project APVV-23-0204 and projects VEGA 1/0307/23, VEGA 1/0363/21.

[1] Žmindák, M. and Dudinský, M. (2012). Finite Element Implementation of Failure and Damage Simulation in Composite Plates. In Eds, Ning HU: Composites and Their Properties. InTEch Publishing. DOI: 10.5772/48248.
[2] Vodička, R., Kormaníková, E. and Kšiňan, F. (2018). Interfacial debonds of layered anisotropic materials using a quasistatic interface damage model with Coulomb friction, International Journal of Fracture, 211(1-2), pp. 163-182. DOI: doi.org/10.1007/s10704-018-0281-z.
[3] Blázquez, A., Mantič, V., París, F. and McCartney, L.N. (2008). Stress state characterization cracks in [0/90] symmetric laminates by BEM. International Journal of Solids and Structures, 45, pp. 1632-1662. DOI: 10.1016/j.ijsolstr.2007.10.013.
[4] Kormaníková, E. and Kotrasová, K. (2017). Delamination modeling of laminate plate made of sublaminates. AIP Conf. Proc. 1863, 260006. DOI: 10.1063/1.4992420.
[5] Raimondo, A., Urcelay Oca, I. and Bisagni, C. (2021). Influence of interface ply orientation on delamination growth in composite laminates, Journal of Composite Materials, 55(27), pp. 1-18. DOI: 10.1177/00219983211031636.
[6] Knopp, A., Funck, E., Holtz, A. and Scharr, G. (2022). Delamination and compression-after-impact properties of zpinned composite laminates reinforced with circumferentially notched z-pins, Composite Structures, 285, 115188. DOI: 10.1016/j.compstruct.2022.115188.
[7] Zani, M., Fanteria, D., Catapano, A. and Montemurro, M. (2022). A consistent energy-based cohesive zone model to simulate delamination between differently oriented plies, Composite Structures, 282, 115042. DOI: 10.1016/j.compstruct.2021.115042.
[8] Saravanan, T.J., Gopalakrishnan, N. and Hari, B.K. (2019). Damage identification in structural elements through curvature mode shapes and nonlinear energy operator, Composite Materials and Engineering, 1(1), pp. 33-48. DOI: 10.12989/cme.2019.1.1.033.
[9] Tao, C.C., Zhang, C., Ji, H.L. and Qiu, J.H. (2021). Fatigue damage characterization for composite laminates using deep learning and laser ultrasonic, Composites Part B-Engineering, 216, 108816. DOI: 10.1016/j.compositesb.2021.108816.
[10] Barbero, E.J. (2014). Finite element analysis of composite materials, CRC Press, Taylor & Francis Group, DOI: 10.1201/b16295.
[11] Sol, H., Gu, J. and Kormanikova, E. (2022). Impulse Excitation Technique for the Determination of Elastic Engineering Constants of CFRP Composites, AIP Conference Proceedings, 2950(1), 195292, DOI: 10.1063/5.0181199.
[12] Saravanan T.J., Gopalakrishnan, N. and Karthick Hari, B. (2019). Damage identification in structural elements through curvature mode shapes and nonlinear energy operator, Composite Materials and Engineering, 1(1), pp. 33-48. DOI: 10.12989/cme.2019.1.1.033.
[13] Gall Trabal, G., Lau Verndal Bak, B., Chen, B., Mosbjerg Jensen, S. and Lindgaard, E. (2022). Delamination toughening of composite laminates using weakening or toughening interlaminar patches to initiate multiple delaminations: A numerical study, Engineering Fracture Mechanics, 273, 108730, DOI: 10.1016/j.engfracmech.2022.108730.
[14] Xie, D. and Waas, A.M. (2006). Discrete cohesive zone model for mixed-mode fracture using finite element analysis. Eng Frac Mech, 73(13), pp. 1783-1796.
[15] Huzni, S., Ilfan M., Sulaiman T., Fonna, M., Ridha, S. and Arifin, A.K. (2013). Finite Element Modeling of Delamination Process on Composite Laminate Using Cohesive Element, In: International Journal of Automotive and Mechanical Engineering (IJAME), 7, pp. 1023-1030, DOI: 10.15282/ijame.7.2012.18.0083.
[16] Rentala, M., Mohinoddin, M. and Venkatesh, S. (2023). Analysis of Delamination Studies in FRP Composites Using 3D Finite Element Simulation, B. Raj et al. (Eds.): ICETE 2023, AER 223, pp. 1088-1097, DOI: 10.2991/978-94-6463-252-1_109.
[17] Liu, Y.F., Qin, Z.Y., Chu, F.L. (2022). Investigation of magneto-electro-thermo-mechanical loads on nonlinear forced vibrations of composite cylindrical shells. Communications in Nonlinear Science and Numerical Simulation, 107, 106146. DOI: 10.1016/j.cnsns.2021.106146.
[18] Elmarakbi, A.M., Hu, N. and Fukunaga, H. (2009). Finite element simulation of delamination growth in composite materials using LS-DYNA, Composites Science and Technology, 69, pp. 2383-2391, DOI: 10.1016/j.compscitech.2009.01.036.
[19] Brunner, A.J., Blackman, B.R.K. and Davies, P. (2008). A status report on delamination resistance testing of polymermatrix composites, Engineering Fracture Mechanics, 75(9), pp. 2779-2794, DOI: 10.1016/j.engfracmech.2007.03.012.
[20] Higuchi, R., Warabi, S., Ishibashi, W. and Okabe, T. (2020). Experimental and numerical investigations on push-out delamination in drilling of composite laminates, Composites Science and Technology, 198, 108238, DOI: 10.1016/j.compscitech.2020.108238.
[21] ASTM D5528-01 (2002): Standard Test Method for Mode I Interlaminar Fracture Toughness of Unidirectional FiberReinforced Polymer Matrix Composites, Annual Book of ASTM Standard.
[22] Gürdal, Z., Haftka, R.T. and Hajela, P. (1999). Design and Optimization of Laminated Composite Materials, J. Wiley & Sons.
[23] Minárová, M. (2024). Rheology, Slovak University of Technology, Bratislava.
[24] Bathe, K.J. (2018). ADINA Primer, ADINA R & D, Inc. 1996 - 2018, USA.
[25] Kormaníková, E., Kšiňan, F. and Vodička, R. (2024). Computational and experimental approaches for Mode I delamination problems, International Journal of Solids and Structures, 300, 112926, DOI: 10.1016/j.ijsolstr.2024.112926.