1. Introduction

A problem of significant interest and importance in solid mechanics is the analysis of delamination phenomena in joined solids. In this sense, composite laminates, materials with inclusions, masonry structures with mortar beds, may be considered as typical joined bodies. The structural response of joined solids strongly depends on the mechanical behaviour of the internal adhesive layers, where several microstructural failure mechanisms lead to progressive degradation of the mechanical properties, up to complete delamination (macrocrack). Moreover, under compressive loading conditions, macrocraks may close and produce frictional effects.

Adhesive layers are generally modelled as zero thickness surfaces by means of interface mechanical device, whose internal constitutive law is defined as a relation between traction and jump-displacement components.

Interface constitutive modelling may be formulated at mesoscale level, reproducing the average mechanical behaviour of an adhesive material portion, known as representative surface element (RSE) whose size is sufficiently large with respect to the microstructure characteristic dimensions.

In the literature, interface constitutive laws have been initially used in the pioneering cohesive-zone models of Dugdale and Barenblatt ^[1,2], to describe the progressive separation process in the fracture process zone (FPZ). These cohesive-zone models are defined in the framework of fracture mechanics, possibly coupled with damage mechanics or with plasticity theory. Interface models in the analysis of composite delamination are developed in the damage mechanics framework ^{[3,4,5,6,7,8,9,10,11]}, where plastic and frictional phenomena are neglected. Alternative formulations ^[12,13,14] are developed in the framework of non-associative softening-plasticity theory, allow to catch the residual frictional behaviour of the fully cracked material.

Cohesive damaging with residual frictional effects has been initially proposed in ^[15] for composite laminates and improved in ^[16,17]. In ^[18] and, by the same author in ^[19,20,21], an original mesoscale interpretation of the damage parameter allows to model the transition from the pure cohesive damaging behaviour of the sound material to the pure frictional behaviour of the fully cracked one. In ^[20], the cohesive-frictional constitutive model has been successfully used for the analysis of frictional effects in Mode Ⅱ delamination tests of composite material, focusing the attention on experimental tests, which were been performed with an increasing cycling loading path. In the present paper the theoretical and computational issues, related to the use of a cohesive-frictional model combined with the kinematic formulation of XFEM, are analysed in detail.

Most of the proposed constitutive models assume the initial behaviour of interface as elastic, by a linear relationship between traction and jump displacement. The normal and tangential elastic parameters are defined in function of the elastic behaviour of the adhesive medium keeping joined the two bodies. For problems of joined solids without any adhesive layer, where bodies are directly bonded one each other, the initial elastic behaviour is assumed as an approximate penalty approach to numerically simulate the mechanical constrain between the bodies.

From the kinematical point of view, the joining surface is modelled as jump displacement between the edges of the joined bodies to describe opening and relative sliding phenomena in the delamination phase. Jump displacements can be active also in the interface sound condition and considered as elastic deformation, which can produce material penetration in compressive state, or elastic delamination in tensile state. In order to bound such effects, a high initial elastic stiffness has to be assumed, but it can produce numerical problems such as spurious traction oscillations, as shown in ^[22].

In the framework of FEM resolution strategy, adhesive layers are generally modelled by means of interface elements, which are defined in a zero thickness domain (surface element for 3D mechanical problems or linear element for 2D mechanical problems) connecting, node by node, the meshes of the joined bodies. Interface finite elements are of simple implementation but, as a drawback, they require two different meshes for the connected bodies with coincident nodes on the connecting surface. Generally, the domain discretization problem by use of interface elements can be quit burdensome.

An alternative formulation is the extended finite element method (XFEM), which allows to define discontinuous displacement fields inside a generic domain, without the need of the mesh to be conforming to the discontinuity surface. The standard displacement based formulation is enriched by a discontinuous displacement field, which is governed by additional degrees of freedom. The spatial position of displacements discontinuity is implicitly defined by the level set method and it can be arbitrarily positioned inside the mesh.

The XFEM has been developed by Belytschko el al. ^[23,24] to model crack formation and propagation, voids and inhomogeneity. Several XFEM schemes have been developed and proposed in literature: for two dimensional linear elastic fracture mechanics (LEFM) in ^[23], with the step function enrichment in ^[24], for three dimensional crack propagation in ^[25], for material discontinuity problems in ^[26,27], with second-order XFEM formulation in ^[27], for cohesive crack models in ^[28,29,30]. The XFEM formulation has also been proposed for the modelling of physically existing discontinuity surface. In ^[31], XFEM and level set method are analysed for the modelling of interface material failure with a damage separation law for the interface, neglecting frictional effects. In ^[32], the XFEM formulation has been employed to model the frictional effects on the contact surface between two adjacent bodies, but neither decohesion nor delamination have been considered. In ^[33], the material microstructure has been modelled by XFEM formulation, including voids and material interface, in a finite element multiscale resolution strategy. In ^[34], material interface has been modelled in a higher order XFEM framework, including curved discontinuity surface, for which attention is focused on the integration procedure and convergence of numerical solutions.

The main target of the present paper is the modelling, in a unified framework, of several mechanical problems such as: frictional contact, delamination, decohesion and smooth transition from the initial cohesive behaviour up to the frictional residual on. Such a complex mechanical behaviour can be modelled, for a joining or discontinuity surface, in the kinematic framework of XFEM, for which the interface can be geometrically placed arbitrarily in a unique mesh, without the requirement of defining two or more meshes with coincident nodes at the joining surface.

The implemented constitutive model has been developed by the same author in ^[19]. The solution proposed in the present paper is limited to two-dimensional problems for which the connection surface is geometrically represented by a line. The finite element formulation is based on the triangular six nodes eXtended element published in ^[30], which allows to model displacement fields with a strong discontinuity at a generic line inside the mesh.

---

2. Extended finite element problem

Consider two solid domains $\Omega_a$ and $\Omega_b$ perfectly bonded and in mutual contact on the common internal surface $\Gamma_c$, with applied external tractions $\textbf{t}$ on the free boundary $\Gamma_t$ and with imposed displacement $\bar{\textbf{u}}$ on $\Gamma_u$, as represented in Figure 1a in a two-dimensional case. The hypotheses of homogeneous material and small strain are assumed. In Figure 1b the overall domain $\Omega=\Omega_a\cup\Omega_b$ is discretized in triangular finite elements and the mesh does not conform the internal connecting surface $\Gamma_c$.

The kinematic formulation developed in ^[30] allows to model a discontinuous displacement field for both the partially cracked finite element ("crack tip element") and for the completely cracked finite element. Moreover both a first-order scheme, with a three node triangular element and a second-order scheme, with a six node triangular element, are developed.

The displacement field of the whole domain $\Omega$ is defined as the additive combination of a continuous part $\textbf{u}_c$ and a discontinuous one $\textbf{u}_d$ :

$\textbf{u}\left(\textbf{x}\right) = \textbf{u}_c\left(\textbf{x}\right)+\textbf{u}_d\left(\textbf{x}\right)$

(1)

where the continuous part is defined by standard shape functions

$\textbf{u}_c\left(\textbf{x}\right) = \sum\limits_{i\in I}N_i\left(\textbf{x}\right)\textbf{u}_i$

(2)

and $I$ is the set of all nodes of the meshed domain and $\textbf{u}_i$ is the displacement vector of the i-th node. With reference to domain $\Omega_{enr}$, defined by the set of elements entirely intersected by the connecting surface $\Gamma_c$, and the relevant set on nodes $I_{enr}\subset I$, the discontinuous part of the displacement field is written in the following form

$\textbf{u}_d\left(\textbf{x}\right) = \sum\limits_{i\in I_{enr}}N_i\left(\textbf{x}\right)\!\Psi_i\left(\textbf{x}\right)\!\textbf{a}_i$

(3)

where $\Psi_i\left(x\right)$ is the enrichment function and $\textbf{a}_i$ contains the relevant degree of freedoms.

The spatial position of the connecting surface $\Gamma_c$ is defined by the level set method, that is

$\Gamma_c = \left\{\textbf{x}\in\Omega: \phi(\textbf{x}) = 0\right\}$

(4)

where function $\phi(\textbf{x})$ assumes different signs in the two domain $\Omega_a$ and $\Omega_b$ and whose iso-zero curve coincides with surface $\Gamma_c$, as shown in Figure 1a. A local reference with tangential and normal axes $(t, n)$ is defined at a generic point of the connecting surface $\Gamma_c$ and the normal axis $n$ is oriented towards the positive values of level set function.

For elements completely intersected by the surface $\Gamma_c$, the enrichment function is defined as

$\Psi_i\left(\textbf{x}\right) = \frac{1}{2}sign\left(\phi(\textbf{x})\right)-\frac{1}{2}sign\left(\phi_i\right), \ \ i\in I_{enr}, \ \ \textbf{x}\in \Omega_{enr}$

(5)

where $\phi_i = \phi(\textbf{x}_i)$ is the level set value at the i-th node.

Due to the specific form of the enrichment functions $\Psi_i(\textbf{x})$, the enrichment displacement field vanishes outside the enriched elements, as shown in Figure 2 for a three node triangular element.

In the finite element framework the interface deformation is defined as a jump displacement between the adjacent edges of the two connected bodies, lying on the surface $\Gamma_c$. The jump displacement across the connecting surface $\Gamma_c$ is defined as the difference of displacement between the two adjacent points on the positive and on the negative part of $\Gamma_c$, that is

$\left[\!\left[\mathbf{u}(\textbf{x})\right]\!\right] = \textbf{u}(\textbf{x}^+)-\textbf{u}(\textbf{x}^-), \ \ \ \ \textbf{x} \in\Gamma_c$

(6)

where $\textbf{x}^+$ and $\textbf{x}^-$ are points respectively on the positive side of $\Gamma_c$ and on its negative side, with reference to the sign of the level-set function. By substitution of Eqs 1, 2, 3 and 5 in Eq 6, it can be observed that jump displacement is function of the enriching parameters $\mathbf{a}_i$, that is

$\left[\!\left[\mathbf{u}(\textbf{x})\right]\!\right] = \sum\limits_{i\in I_{enr}}N_i\left(\textbf{x}\right)\!\textbf{a}_i, \ \ \ \ \textbf{x} \in\Gamma_c$

(7)

In the present paper, the initial behaviour of the interface is assumed elastic and the jump displacements can be different than zero also in the undamaged state of the connecting surface. The enriching degrees of freedom $\textbf{a}_i$ with $i\in I_{enr}$, which belong to nodes contained in the domain $\Omega_{enr}$ of elements crossed by the interface, are always unconstrained and some numerical problems, such as spurious traction oscillations, can emerge (see ^[22]).

In ^[31], two different formulations are considered respectively for the sound interface and for the cracked one. In the former case, the finite element model is enhanced with a weak discontinuity field, in order to model the kinematical behaviour at the initial state, when the two bodies are perfectly bonded with zero jump displacements and with independent deformation states on the two element parts. Finally, when in a finite element the internal interface attains the delamination condition, the element is enriched with a strong discontinuity field by releasing the relevant degrees of freedom. Nevertheless, in such formulation numerical problem of stress singularity or high stress concentration are observed.

---

2.1. Position of surface discontinuity

The position inside the mesh of the connecting surface is defined in the input file by means of linear or quadratic parametric equations, or circle equations. Inside a single finite element the connecting surface is linearised and the following intersection schemes are considered:

a) side to side, when connecting surface intersects two element sides (Figure 2);

b) node to side, when connecting surface contains one node and intersects its opposite side;

c) node to node, when connecting surface contains two vertex nodes;

d) one node, when connecting surface contains just one vertex node, without cutting the element internal domain.

In (a) and (b) schemes, the finite element is completely cut, whereas in scheme (c) the interface coincides with a finite element side.

---

3. Cohesive-frictional constitutive model for interface

The interface constitutive model proposed in ^[19] is adopted for the numerical analysis; the model is developed in two-dimensional space, with reference to a zero thickness linear interface and it is described in the present section.

Interface constitutive models typically link the displacement discontinuity $\left[\!\left[\mathbf{u}(\textbf{x})\right]\!\right]$ and the adhesion forces $\textbf{t}(\textbf{x})$. As described in previous section, the kinematics is defined by the displacement discontinuity across the interface, and it can be decomposed in two components: the normal one (Mode Ⅰ) and the tangential one (Mode Ⅱ). The interface equilibrium condition imposes continuity of traction $\textbf{t}$ across interface, therefore it is assumed $\mathbf{t} = \mathbf{t}^{+} = \mathbf{t}^{-}$, where $\mathbf{t}^{+}$ and $\mathbf{t^{-}}$ are the traction vectors, respectively, on the positive edge and on the negative one. Statics and kinematics of an infinitesimal element $dS$ of the connecting surface $\Gamma_c$ are represented in Figures 3a and b.

The proposed constitutive model is developed in the damage mechanics framework ^[35] considering the Representative Surface Element (RSE), whose size is assumed sufficiently large compared to the adhesive material inhomogeneities (Figure 4), and average damage in the RSE is defined as

$\omega = \frac{d S_{c}}{d S}$

(8)

where $dS_{c}$ measures the area of the micro-cracked RSE portion and $dS$ measures the RSE area.

The damage variable $\omega$ measures the specific areal crack density of the infinitesimal surface element $dS$ and it allows us to define, at the material point, the extensions of the two complementary fractions: the micro-cracked one $\omega\, dS$ and of the sound one $\left(1-\omega\right)dS$ (see Figure 4). Two different constitutive laws can be formulated for the two fractions and the damage variable governs the progressive transition of the adhesive constitutive model from the initial elastic-cohesive behaviour to the residual frictional one. In particular, a cohesive linear elastic law is assumed for the sound fraction and a non-associative elastic-plastic constitutive relation, with Coulomb yield surface, is assumed for the cracked fraction. The compressive behaviour of both the sound fraction and the damaged one, is modelled as linearly elastic.

The constitutive model formulation is developed considering the index $\mathit{s}$ for the sound fraction variables and index $\mathit{c}$ for the micro-cracked fraction ones. For the sound fraction the elastic deformation is define as

$\boldsymbol{\delta}_{s}^{e} = \left[\!\left[\mathbf{u}\right]\!\right]- \boldsymbol{\delta}_{s}^{p}$

(9)

where $\boldsymbol{\delta}_{s}^{p}$ measures the plastic contribution. For the cracked fraction, which is represented in Figure 5a, the variable $\boldsymbol{\delta}_{c}^{e}$ measures the asperities elastic deformation, represented in Figure 5b, and it is defined by the following relation

where $\boldsymbol{\delta}_{c}^{p}$ measures the plastic component (sliding and dilatancy represented in Figure 5c) and $\boldsymbol{\delta}_{c}^{d}$ is the detachment jump displacement (Figure 5d).

The vector $\boldsymbol{\delta}_{c}^{p}$ represents the jump displacement caused by the irreversible sliding of asperities among each others and it develops when the limit of frictional strength is attained. Moreover, since frictional behaviour is active only for compressive normal traction, positive normal traction produces loss of contact between the crack edges, which is measured by the detachment displacement $\boldsymbol{\delta}_{c}^{d}$.

The constitutive model is developed in a thermodynamic framework following the classic Coleman and Noll procedures (for more detail see ^[35,36,37]). The Helmholtz free energy density (for unit interface length) is defined as

$\begin{eqnarray} \Psi & = & \Psi_{el}\left(\boldsymbol{\delta}_{s}^{e}, \, \boldsymbol{\delta}_{c}^{e}, \, \omega\right)+ \Psi_{in}\left(\xi\right) = \nonumber \\ & = & \frac{1}{2}\left(1-\omega\right) \boldsymbol{\delta}_{s}^{eT}\mathbf{K}_{\mathit{s}}\, \boldsymbol{\delta}_{s}^{e}+\frac{1}{2}\omega \boldsymbol{\delta}_{c}^{eT}\mathbf{K}_{\mathit{c}}\, \boldsymbol{\delta}_{c}^{e}+\Psi_{in}\left(\xi\right) \label{eq:Helm} \end{eqnarray}$

(11)

where $\Psi_{el}$ and $\Psi_{in}$ are respectively the elastic part and the internal one; $\mathbf{K}_{\mathit{s}}$ and $\mathbf{K}_{\mathit{c}}$ are diagonal elastic matrices, respectively, of the sound fraction and of the micro-cracked one; $\xi$ is a scalar state internal variable, which governs damage hardening. The second thermodynamic principle, in the form of Clausius-Duhem inequality, imposes a pointwise positiveness of rate dissipation energy density given as the difference between the work done by external forces and the variation of the free energy, for any change of the state variables. This condition is mathematically written as

$D = \mathbf{t}^{T}\left[\kern-0.15em\left[\dot{\mathbf{u}} \right]\kern-0.15em\right]- \dot{\Psi} \geq0,$

(12)

where rate of Helmholtz free energy is defined as

$\begin{eqnarray} \dot{\Psi} & = & \frac{\partial\Psi_{el}}{\partial\omega}\dot{\omega} + \frac{\partial\Psi_{el}}{\partial\boldsymbol{\delta}_{s}^{e}}\cdot \dot{\boldsymbol{\delta}}_{s}^{e} + \frac{\partial\Psi_{el}}{\partial\boldsymbol{\delta}_{c}^{e}}\cdot \dot{\boldsymbol{\delta}}_{c}^{e}+\frac{\partial\Psi_{in}}{\partial\xi}\dot{\xi} = \label{eq:Helm rate} \\ & = & -\left(Y\dot{\omega}- \mathbf{t}_{s}^T\dot{\boldsymbol{\delta}_{s}^{e}}- \mathbf{t}_{c}^{T} \dot{\boldsymbol{\delta}_{c}^{e}} - \chi\dot{\xi}\right) \nonumber \end{eqnarray}$

(13)

$\mathbf{t}_{s}$ and $\mathbf{t}_{c}$ are conjugated variable of the elastic deformations $\boldsymbol{\delta}_{s}^{e}$ and $\boldsymbol{\delta}_{c}^{e}$, that is

$\mathbf{t}_{s} = \frac{\partial\Psi_{el}}{\partial\boldsymbol{\delta}_{s}^{e}} = \left(1-\omega\right)\mathbf{K}_{\mathit{s}}$

(14)

$\mathbf{t}_{c} = \frac{\partial\Psi_{el}}{\partial\boldsymbol{\delta}_{c}^{e}} = \omega\mathbf{K}_{\mathit{c}}\boldsymbol{\delta}_{c}^{e},$

(15)

$Y$ is the energy release rate and $\chi$ is the static hardening-like variable

$Y = -\frac{\partial\Psi_{el}}{\partial\omega} = \frac{1}{2}\boldsymbol{\delta}_{s}^{eT} \mathbf{K}_{s}\, \boldsymbol{\delta}_{s}^{e}- \frac{1}{2}\boldsymbol{\delta}_{c}^{eT}\mathbf{K}_{c} \, \boldsymbol{\delta}_{c}^{e}$

(16)

$\chi = \frac{\partial\Psi_{in}}{\partial\xi}.$

(17)

For an elastic loading step, where $\dot{\omega} = 0$, $\dot{\xi} = 0$, $\dot{\boldsymbol{\delta}_{s}^{p}} = \mathbf{0}$ and $\dot{\boldsymbol{\delta}_{c}^{p}} = \mathbf{0}$, rate elastic deformations are $\dot{\boldsymbol{\delta}_{s}^{e}} = \dot{\boldsymbol{\delta}_{c}^{e}} = \left[\kern-0.15em\left[\dot{\mathbf{u}} \right]\kern-0.15em\right]$, dissipation has to be null ($D = 0$) and after substitution Eq 12 reads

$D = \left(\mathbf{t}- \mathbf{t}_{s}- \mathbf{t}_{c}\right)^{T} \left[\kern-0.15em\left[\dot{\mathbf{u}} \right]\kern-0.15em\right] = 0.$

(18)

Since Eq 18 has to be verified for any elastic loading step $\left[\kern-0.15em\left[\dot{\mathbf{u}} \right]\kern-0.15em\right]$, the following internal equilibrium condition is obtained

$\mathbf{t} = \mathbf{t}_{s}+ \mathbf{t}_{c}$

(19)

and, after substitution of Eqs 14 and 15, the overall interface elastic traction-displacement relation is obtained

$\mathbf{t} = (1-\omega)\mathbf{K}_{s} \boldsymbol{\delta}_{s}^{e}+\omega\mathbf{K}_{c} \boldsymbol{\delta}_{c}^{e}.$

(20)

Equation 20 is assumed to be valid also for inelastic loading steps, in which dissipation is

$D = \mathbf{t}_{s}^{T}\dot{\boldsymbol{\delta}}_{s}^{p}+ \mathbf{t}_{c}^{T}\dot{\boldsymbol{\delta}}_{c}^{p}+ Y\dot{\omega}-\chi\dot{\xi}\geq0.$

(21)

where detachment $\boldsymbol{\delta}_{c}^{d}$ does not produce any dissipation.

Dissipation positiveness is ensured by means of suitable activation criteria, which govern the evolution of plastic and damaging phenomena. In particular, damage evolution is governed by the following yield function

$\phi_{d}\left(Y, \chi\right) = Y-\chi\left(\xi\right)-Y_{0}\leq0,$

(22)

with $Y_{0}$ assumed as initial yielding threshold. The relevant flow rules with the loading-unloading conditions read

$\begin{eqnarray} \dot{\omega} & = & \frac{\partial\phi_{d}}{\partial Y}\, \dot{\lambda}_{d} = \dot{\lambda}_{d}, \nonumber \\ \dot{\xi} & = & -\frac{\partial\phi_{d}}{\partial\chi}\, \dot{\lambda}_{d} = \dot{\lambda}_{d}, \label{eq:equazioni di flusso danno} \\ \dot{\lambda}_{d} & \geq & 0, \;\; \phi_{d}\dot{\lambda}_{d} = 0, \;\; \dot{\phi}_{d}\dot{\lambda}_{d} = 0. \nonumber \end{eqnarray}$

(23)

The damage hardening-like law is

$\chi\left(\xi\right): = \frac{1}{2}K_{N}^{s}\bar{u}_{e}^{2} \left[\left(\frac{\bar{u}_{f}}{\bar{u}_{f} \left(1-\xi\right)+\bar{u}_{e}\xi}\right)^{2}-1\right],$

(24)

where $\bar{u}_{e}$ and $\bar{u}_{f}$ are jumps displacement limit values, respectively, at the elastic threshold and at the unitary damage condition in a pure tensile state and $K_{N}^{s}$ is the normal elastic stiffness of the interface sound fraction. Hardening law in Eq 24 produce linear softening in the traction-jump displacement space and it produces the following simple Mode Ⅰ fracture energy

$G_{I} = \frac{1}{2}K_{N}^{s}\bar{u}_{e}\bar{u}_{f}.$

(25)

Plastic phenomena in the sound fraction are neglected (${\boldsymbol{\delta}_{s}^{p}} = \mathbf{0}$), whereas, in micro-cracked fraction, plastic phenomena evolution is governed by means of the classical Coulomb yield function, generally used for frictional materials,

$\phi_{p} = \left|t_{cT}\right|+\alpha\, t_{cN}\leq0.$

(26)

under the hypothesis of non-associative plasticity, by means of the following plastic potential

$\Omega_{p} = \left|t_{cT}\right|+\beta\, t_{cN},$

(27)

where $\alpha$ and $\beta$, with $\alpha \geq \beta$, are respectively the frictional coefficient and the dilatancy one. The plastic flow rules and loading/unloading conditions are

$\begin{eqnarray} \dot{\delta}_{cT}^{p} & = & \frac{\partial\Omega_{p}}{\partial t_{cT}}\, \dot{\lambda}_{p} = \textrm{sgn}\left(t_{cT}\right)\dot{\lambda}_{p}, \nonumber \\ \dot{\delta}_{cN}^{p} & = & \frac{\partial\Omega_{p}}{\partial t_{cN}}\, \dot{\lambda}_{p} = \beta\dot{\lambda}_{p}, \label{eq:leggi di flusso plastiche} \\ \dot{\lambda}_{p} & \geq & 0, \;\;\phi_{p}\dot{\lambda}_{p} = 0, \;\; \dot{\phi}_{p}\dot{\lambda}_{p} = 0. \nonumber \end{eqnarray}$

(28)

By substitution of Eqs 28 and 23 in Eq 21, dissipation assumes the following form

$D = \left(\left|t_{cT}\right|+\beta\, t_{cN}\right) \dot{\lambda}_p+ Y_0\dot{\lambda}_d\geq 0$

(29)

and positiveness is ensured for any damaging or plastic loading step. The cohesive-frictional interface constitutive model has been implemented in the open source finite element code FEAP ^[38]. The interface kinematics is implicitly defined inside the triangular six nodes extended finite element. The enrichment degrees of freedom are active for the nodes of all the finite elements crossed the joining surface $\Gamma _C$; on the contrary, the other enrichment degrees of freedom are restrained. Three Gauss integration points are considered for each finite element, where the two components of the jump displacement are computed through Eq 7 and passed to the subroutine of the interface constitutive model, which give back the tractions components and the constitutive stiffness matrix. Details on the implementation of both user defined finite elements and user defined constitutive models can be found in ^[38].

---

4. Numerical simulations

The proposed interface constitutive model and the triangular extended finite element have been implemented in the open source finite element code FEAP ^[38]. The numerical results of delamination phenomenon in joined bodies are presented.

The first simulation is the two-dimensional analysis of delamination between a square matrix the the circular inclusion, represented in Figure 6 with dimensions and boundary conditions. Matrix and inclusion are considered as elastic isotropic with the following constitutive parameters: $E_M = 40000\ N/mm^2$, $\nu_M = 0.2$, $E_I = 400000 \ N/mm^2$, $\nu_I = 0.2$. With reference to the proposed interface constitutive model, the elastic parameters of the sound fractions are: $K_N^s = 400000\ N/mm, K_T^s = 347830\ N/mm$; the elastic parameters of the cracked fractions are: $K_N^c = 400000\ N/mm, K_T^c = 6956.6\ N/mm$; the jump displacement limit values are $\bar{u}_{e} = 10^{-6}\ mm, \bar{u}_{f} = 5\times10^{-5}\ mm$; frictional and dilatancy coefficients are $\alpha = 0.364, \beta = 0.176$.

The domains of the two joined solids are numerically modelled by the not conforming mesh of XFEM triangular elements represented in Figure 6c; they are also modelled by quadrilateral isoparametric elements and interface elements in the mesh represented in Figure 6b, which conforms to the joining surface between matrix and inclusion. The numerical responses carried out by the two approaches, XFEM and classical interface elements, are compared in terms of force vs displacement in Figure 7 and the results are numerically coincident. Figures 8a and b plot the maps of vertical normal stress carried out, respectively, at the reference loading steps $a$ and $b$ represented in Figure 7.

The second proposed simulation is the two-dimensional sliding test between two elastic solids in mutual contact along a curvilinear surface, without any cohesive strength. The configuration test is represented in Figure 9a together with the adopted mesh, which is defined for the whole domain and it does not conform the contact surface. The hypothesis of small displacement is assumed. In the considered sliding test, the lower elastic solid is constrained along the lower boundary side and the upper elastic solid is subjected to a constant vertical force $F_y = 1000$ N and to a monotonically increasing horizontal displacement at the upper boundary side, as shown in Figure 9a. Rotation of the upper boundary side is restrained. The elastic parameters of the two solids are $E = 40000\ N/mm^2$, $\nu = 0.2$; for the contacts surface, frictional and dilatancy coefficients respectively are $\alpha = tan20^\circ$, $\beta = 0$. The map of vertical stresses, obtained by the numerical simulation for the imposed horizontal displacement $u_x = 3.0$ mm, is plotted in Figure 9b. The results of the numerical simulation are plotted in Figure 10 in terms of horizontal displacement vs horizontal force. The second numerical simulation shows that the proposed formulation allows to model problems of frictional contact with a very simple mesh generation, where a quite complex contact surface is defined in the input file as an isoparametric quadratic function.

The last numerical simulation reproduces the experimental results of a three points bending test performed on a carbon fibres-epoxy composite beam, whose dimensions are represented in Figure 11. The composite panel has been manufactured by hand lay-up and vacuum-bag techniques by a sequence of 12 unidirectional plies of carbon fibres. A Mylar foil was inserted at the mid-plane of the composite panels during the lay-up process in order to produce a starter delamination crack. Four specimens have been tested under displacement control and simultaneously the crack extension has been monitored by an high resolution camera. The experimental test results are plotted in Figure 12 in terms of imposed displacement vs vertical reaction. The fracture energy values obtained for the four experimental tests are reported in Table 1 with the relevant mean value. The three points bending test has been numerically simulated by using a two-dimensional mesh of $9 \times 110$ six nodes extended triangular elements in plane strain condition. The composite material is assumed as elastic isotropic with the following constitutive parameters: $E = 80580\ N/mm^2$ and $\nu = 0.23$. The interface constitutive parameters are: $K_N^s = K_T^s = 90000\ N/mm$, $K_N^c = 90000\ N/mm, K_T^c = 3000\ N/mm$; the jump displacement limit values are $\bar{u}_{e} = 6.9\times10^{-4} \ mm, \bar{u}_{f} = 0.014593 \ mm$; frictional and dilatancy coefficients respectively are $\alpha = tan35^\circ$, $\beta = 0$. The adopted constitutive parameters produces the fracture energy value $G_{II} = \frac{1}{2}K_N^s\cdot \bar{u}_{e} \cdot \bar{u}_{f} = 0.4531\ N/mm$. The results of the numerical simulation are compared to the experimental results in Figure 12 in terms of imposed displacement vs load P, showing a good agreement between the two results. The numerical simulation produces: the same stiffness in the elastic branch; the same values of load and imposed displacement at the delamination condition; a response close to the experimental data also in the post-delamination branch. Moreover, the set of parameters chosen for the cohesive law produce a mode Ⅱ fracture energy very close to the values obtained by the experimental data. The discretization of the pre-cracked specimen, by XFEM formulation, has been obtained by the generation the mesh of the overall rectangular domain and by defining, in the input file, the line of the pre-existing crack, with unitary damage, and the line of possible delamination. The discretization of the same problem with standard interface element would be extremely much more onerous.

Finally, the map of horizontal normal stress at the maximum loading condition is plotted in Figure 13 for the whole specimen together with a zoom around the crack tip.

The first numerical simulation proposed in the present paper has been discretized either by the standard interface elements and by the extended finite elements and the results of the two numerical simulations are practically coincident, showing the effectiveness of the implemented XFEM formulation.

Whereas interface finite elements are of very simple implementation in a finite element code, the most difficult in their use is the discretization phase, which requires a separate mesh for each connected body, with coincident nodes of the discretized bodies on the joining surface. The domain discretization phase, by use of interface elements, can be quite burdensome. Some alternative formulations of non-matching interface elements have been proposed in literature ^[39,40], which allow the discretization of the two connected bodies with independent meshes. The two meshes can be defined with non-coincident nodes on the joining surface $\Gamma_C$, with an easier discretization phase.

On the contrary, the principal advantage of the XFEM formulation in the modelling of different bodies connected by joining surface is the discretization. The whole domain can be discretized with a unique mesh and the position of joining surface can be analytically defined, with much lower time needed for such operation, with respect to the interface formulation. Moreover, the XFEM formulation allows also to model extrinsic cohesive law, with initial rigid behaviour. Disadvantages of the XFEM formulation is the implementation, which can be very onerous, and the mesh density of the joining surface, which cannot be easily controlled for each element.

---

5. Conclusions

The present paper proposes a unified approach for the modelling of friction and delamination phenomena on discontinuity surfaces, discretized in an overall mesh based on an XFEM formulation. The proposed approach allows us to model both the delamination phenomenon with crack propagation of two joined solids and the frictional contact phenomenon between two adjacent solids. The transition, from the initially cohesive behaviour with delamination phenomena, to the final behaviour with residual frictional strength is also considered.

The discontinuity surfaces has been arbitrarily defined inside the mesh in parametric formulation, a unique mesh can be defined for whole domain, including connected bodies and discontinuity surface. On the contrary, in a standard approach, separated meshes have to be defined for the connected or adjacent bodies and specific interface finite elements for the discontinuity surface. Moreover, the meshes of the adjacent or connected bodies have to be defined with spatially coincident nodes on the discontinuity surface.

---

6. Acknowledgments

Funding: The financial support of the Italian Ministry for University and Research (MIUR), under the grant PRIN-2015, project No.2015LYYXA8, "Multiscale mechanical models for the design and optimization of microstructured smart materials and metamaterials" is gratefully acknowledged.

---

7. Conflict of interest

The authors declare that they have no conflict of interest.

---