On the Use of Drilling Degrees of Freedom to Stabilise the Augmented Finite Element Method

Abstract

The augmented finite element method (AFEM) embeds cracks within solid elements. These cracks are modelled without additional degrees of freedom thanks to a dedicated static condensation process. However, static condensation can induce a lack of constraint problem, resulting in singular stiffness matrices. To address this issue, we propose a new method called the stabilised augmented finite element method (SAFEM), which produces non-singular stiffness matrices. We conducted 2D experiments involving stationary traction-free cracks and propagating cohesive discontinuities to compare the performance of the SAFEM with the AFEM. The SAFEM outperforms the AFEM in modelling traction-free cracks.

1. Introduction

In-plane rotational degrees of freedom (DOFs) are often referred to as “drilling degrees of freedom” in membrane finite elements. Their development started in the 1960s and served two purposes. Firstly, drilling DOFs improve the performance of linear membrane elements by allowing the manufacture of elements of intermediate computational cost and accuracy between linear and quadratic elements [1,2,3,4]. Secondly, when classical membrane elements and standard plate-bending elements are combined to form a shell element, each node possesses three displacement DOFs and two rotational DOFs in the local coordinate system of the shell. The presence of a third rotational DOF is attractive because it prevents the singularity of the stiffness matrix that occurs if shells meeting at a node are co-planar as outlined in Chap. 13.5 of Ref. [5].

In this paper, we explore a novel application of drilling DOFs: stabilising the augmented finite element method (AFEM). The AFEM was developed by Yang and coworkers as a numerical technique to embed weak or strong discontinuities within standard finite elements [6,7]. By embedding discontinuities inside elements, the AFEM simplifies their modelling because they no longer need to conform to the elements boundaries as required in the finite element method (FEM). The AFEM has been widely used in structural analysis to model cracks, with successful applications, including unstable crack propagation [8], crack growth under thermomechanical loading [9], dynamic crack propagation [10], three-dimensional studies of composite materials [11] or the large deformation of cracked shells [12]. The method was implemented as a user element in Abaqus and proved to be 50 times faster than the phantom-node method (PNM) natively implemented in this software [13]. The AFEM offers a significant advantage over other partition of unity methods, such as the PNM or the extended finite element method [14,15] (XFEM): it does not require adding global unknowns to model discontinuities. This allows the AFEM to represent an arbitrary number of growing cracks while keeping the size of the global stiffness matrix constant. However, the AFEM does have a major drawback: it requires handling stiffness matrices that become singular in some situations [6,16]. The AFEM originators briefly discussed this issue in Ref. [6], but no satisfactory solution was provided. We thus previously proposed a dedicated procedure to handle these singular matrices [16], which is recalled in Section 3.4. As far as we know, our proposed strategy is the only solution for using augmented elements with traction-free cracks published to date.

The use of drilling DOFs proposed in this work improves the AFEM by preventing the occurrence of singular stiffness matrices, eliminating the need for a dedicated procedure to handle them. Additionally, drilling DOFs substantially increases the accuracy of the AFEM. The improved AFEM with drilling DOFs will be referred to as the stabilised augmented finite element method (SAFEM) in this work. Numerical experiments conducted with an in-house code [16,17] demonstrate that the SAFEM outperforms the AFEM when cracks are traction-free. The results obtained with cohesive cracks are of similar accuracy.

An outline of this paper is as follows. Section 2 provides a brief review of the formulation of finite elements with drilling DOFs and includes some accuracy checks. Section 3 presents a summary of the general AFEM framework for modelling cohesive cracks and demonstrates how drilling DOFs can be introduced within augmented elements. In Section 4, we perform multiple experiments to evaluate the accuracy of the SAFEM and compare it to the original AFEM.

2. Finite Elements with Drilling Degrees of Freedom

According to MacNeal and Harder, membrane finite elements with drilling degrees of freedom (DOFs) were proposed several times between 1960 and 1980, but these attempts were mostly unsuccessful [3]. The distinct approaches taken by Bergan and Felippa [18] and Allman [19], as well as the mathematical framework proposed by Hughes and Brezzi [20] paved the way for the successful formulations used today (see, for example, [21,22,23]). In this work, we will consider the approach proposed by Allman in 1984 to design membrane elements with in-plane rotational DOFs [19]. We chose Allman’s elements because of their simplicity, but other formulations could have been considered as well. To keep the exposition to a reasonable length, only 2D membrane elements will be investigated, but Allman’s approach is equally applicable to 3D elements, as exemplified in Refs. [24,25]. We begin by recalling the formulation of Allman’s elements and check their accuracy. A word of caution is in order: Allman derived other elements with drilling DOFs in 1988 [2], which are different from those considered here; indeed, they are incompatible. Hence, one must bear in mind that this work only employs Allman’s 1984 compatible elements.

The goal of the present section is to benchmark Allman’s elements against “classical" finite elements. It is a fundamental step. Indeed, the SAFEM uses Allman’s elements and will inherit their properties, such as the convergence rate and sensitivity to distortions. To the best of our knowledge, the results we are about to present regarding Allman’s 1984 elements were yet to be found in the literature.

2.1. Design and Formulation of Allman’s Elements

Allman derived a triangular element with drilling DOFs through intuitive geometrical concepts [19]. Cook and coworkers re-interpreted Allman’s work as the result of applying a coordinate transformation to elements with midside nodes and applied it to quadrilateral, tetrahedral and hexahedral elements [1,24,25]. The approach can be summarised as follows: one starts with a higher-order isoparametric element with straight edges and midside nodes (e.g., the six-node linear strain triangle or the eight-node Serendipity quadrangle). Then, one “transforms” the midside translational DOF into drilling DOFs located at the element’s corners through master–slave constraints. Felippa called the resulting operation a “retrofitting” procedure in Ref. [21].

This section focuses on the retrofitting fabrications of the following: (i) a three-node triangle with drilling DOFs, coined T3A (where “A” stands for Allman), from a six-node linear strain triangle (T6) and (ii) a four-node quadrilateral with drilling DOFs (Q4A) from an eight-node biquadratic quadrilateral (Q8), see Figure 1.

Allman’s elements are obtained by choosing specific polynomial orders for the components of the displacement vector along each edge of finite elements with midside nodes: the normal and tangential components of the displacement field along the edges, $u_n$ and $u_t$, are quadratic and linear, respectively. Following Allman’s approach, the expressions for $u_n$ and $u_t$ are [19]:

$$\begin{array}{c}{u}_n\left(s\right)=(1-\frac{s}{l_{ij}})u_{n_i}+\frac{s}{l_{ij}}{u}_{n_j}+\frac{4s}{l_{ij}}(1-\frac{s}{l_{ij}})u_{n_{ij}}\end{array}$$

(1)

$$\begin{array}{c}{u}_t\left(s\right)=(1-\frac{s}{l_{ij}})u_{t_i}+\frac{s}{l_{ij}}{u}_{t_j}\end{array}$$

(2)

where $l_{ij}$ is the length of the edge between nodes i and j, s is the curvilinear abscissa along the edges of the elements and i and j are indexes associated with the element nodes (cf. Figure 1). For instance, $u_{n_j}$ is the degree of freedom associated with the normal component of the displacement field at node j. The use of both indexes i and j means that one considers the degree of freedom associated with a midside node, e.g., $u_{t_{ij}}$, is the DOF associated with the tangential component of the displacement field at the midside node located between nodes i and j. The degrees of freedom supported by the midside nodes (i.e., $u_{n_{ij}}$ in Equation (1)) are then replaced by rotational DOFs supported by the apex nodes. Vertex rotations are introduced for that purpose and the following relation is made use of:

$$\begin{array}{c}\hfill u_{n_{ij}}=\frac{l_{ij}}{8}({\mathsf{\Theta }}_j-{\mathsf{\Theta }}_i)\end{array}$$

(3)

where ${\mathsf{\Theta }}_i$ are the so-called drilling DOFs representing the vertex rotation at node i. Despite the use of drilling DOFs, Allman’s formulation does not seek to interpolate a rotation field. Instead, it only employs vertex rotations to enhance the displacement field. This is in contrast to elements with drilling DOFs, which interpolate the rotation field and displacement field independently (as seen in [20,23]). Furthermore, it is important to note that the vertex rotations introduced in Allman’s approach are not true rotations as defined in continuum mechanics, but they are closely related. Indeed, the following relation holds for the T3A [19]:   

$$\begin{array}{cc}\hfill {\mathsf{\Omega }}_i-\frac{1}{3}({\mathsf{\Omega }}_1+{\mathsf{\Omega }}_2+{\mathsf{\Omega }}_3)& =\frac{3}{4}\left({\mathsf{\Theta }}_i-\frac{1}{3}({\mathsf{\Theta }}_1+{\mathsf{\Theta }}_2+{\mathsf{\Theta }}_3)\right)\mathrm{for}i\in [1,2,3]\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill \mathrm{with}\mathsf{\Omega }(x,y)& =\frac{1}{2}\left(\frac{\partial v(x,y)}{\partial x}-\frac{\partial u(x,y)}{\partial y}\right)\hfill \end{array}$$

(5)

where $\mathsf{\Omega }$ is the true rotation field; ${\mathsf{\Omega }}_i$ is the value of the true rotation field at node i; and u and v are the horizontal and vertical components of the displacement field, respectively. Equation (4) reveals that the differences between the true rotations at the vertices (${\mathsf{\Omega }}_i$) and the average true rotation are directly proportional to the discrepancy between the vertex rotations (${\mathsf{\Theta }}_i$) and their mean. Equations (1)–(3) highlight how the translational DOF supported by the midside nodes of higher-order elements can be turned into vertex rotation DOFs. Consequently, the midside nodes and their associated DOF are no longer required (cf. Figure 1) and the following relations hold:

$$\begin{array}{c}{\{u_1,v_1,\dots ,u_6,v_6\}}^t=\left[T_{\mathrm{T}6\text{-}\mathrm{T}3\mathrm{A}}\right]{\{u_1,v_1,{\mathsf{\Theta }}_1,\dots ,u_3,v_3,{\mathsf{\Theta }}_3\}}^t\end{array}$$

(6)

$$\begin{array}{c}{\{u_1,v_1,\dots ,u_8,v_8\}}^t=\left[T_{\mathrm{Q}8\text{-}\mathrm{Q}4\mathrm{A}}\right]{\{u_1,v_1,{\mathsf{\Theta }}_1,\dots ,u_4,v_4,{\mathsf{\Theta }}_4\}}^t\end{array}$$

(7)

where $\left[T_{\mathrm{T}6\text{-}\mathrm{T}3\mathrm{A}}\right]$ and $\left[T_{\mathrm{Q}8\text{-}\mathrm{Q}4\mathrm{A}}\right]$ are the so-called transformation matrices that solely depend on the coordinates of the element vertices [1,3]. Closed-form expressions for these matrices are provided in Appendix A. Once the transformation matrices are known, the stiffness matrices of Allman’s elements can be readily computed from the stiffness matrices of the higher-order elements they are based on. Specifically, if $\left[K_{\mathrm{T}6}\right]$, $\left[K_{\mathrm{Q}8}\right]$, $\left[K_{\mathrm{T}3\mathrm{A}}\right]$ and $\left[K_{\mathrm{Q}4\mathrm{A}}\right]$ are the stiffness matrices of the T6, Q8, T3A and Q4A, respectively, then:

$$\begin{array}{c}\left[K_{\mathrm{T}3\mathrm{A}}\right]={\left[T_{\mathrm{T}6\text{-}\mathrm{T}3\mathrm{A}}\right]}^t\left[K_{\mathrm{T}6}\right]\left[T_{\mathrm{T}6\text{-}\mathrm{T}3\mathrm{A}}\right]\end{array}$$

(8)

$$\begin{array}{c}\left[K_{\mathrm{Q}4\mathrm{A}}\right]={\left[T_{\mathrm{Q}8\text{-}\mathrm{Q}4\mathrm{A}}\right]}^t\left[K_{\mathrm{Q}8}\right]\left[T_{\mathrm{Q}8\text{-}\mathrm{Q}4\mathrm{A}}\right]\end{array}$$

(9)

The elements resulting from the use of Equations (8) and (9) are compatible, satisfy the patch tests, employ second-order shape functions and only involve DOFs located at corner nodes. Allman’s elements can be integrated using the same quadrature rules as their underlying higher-order elements. In this study, we integrate the stiffness matrices of the T3A and the T6 with three integration points, while the Q4A and the Q8 make use of a 3 × 3 Gaussian integration scheme. These quadratures fully integrate the corresponding higher-order elements, but the associated Allman’s elements are rank deficient. There exists a spurious deformation mode with identical vertex rotations (i.e., ${\mathsf{\Theta }}_1={\mathsf{\Theta }}_2=\dots ={\mathsf{\Theta }}_n$) whatever the chosen integration scheme [3]. Several strategies have been proposed to deal with it: (i) constraining a single drilling DOF of the entire mesh to prevent the appearance of the spurious mode or (ii) adding a penalty matrix to the stiffness matrices in order to restore their full rank [3,4,24,25]. These two options will be investigated. The penalty stiffness used to prevent rank deficiencies of Allman’s elements was designed by Macneal and coworkers and is given by (see Equations (46)–(49) of Ref. [3]):

$$\begin{array}{c}\left[K_{ij}^{penalty}\right]=\gamma VG\frac{{\partial }^2\left({\mathsf{\Theta }}_{\Delta }^2\right)}{\partial q_i\partial q_j}\end{array}$$

(10)

$$\begin{array}{c}\mathrm{with}{\mathsf{\Theta }}_{\Delta }=\frac{1}{n}\sum _{i=1}^n\left({\mathsf{\Theta }}_i-{\mathsf{\Omega }}_0\right)\mathrm{and}{\mathsf{\Omega }}_0=\frac{1}{2}{\left(\frac{\partial v(x,y)}{\partial x}-\frac{\partial u(x,y)}{\partial y}\right)}_0\end{array}$$

(11)

where $\gamma$ is the dimensionless penalty parameter, V is the element volume, G is the shear modulus, n is the number of nodes of Allman’s elements, ${\mathsf{\Theta }}_i$ is the drilling DOF at node i, $q_i$ is the ith component of the DOF vector (e.g., $q_1=u_1$, $q_2=v_1$ or $q_3={\mathsf{\Theta }}_1$, cf. Equation (6)) and ${\mathsf{\Omega }}_0$ is the true rotation computed at the centroid of the element thanks to the shape function derivatives. Unlike the stiffness matrices, the penalty matrix is computed without using quadrature rules and instead employs a closed-form expression. Equations (10) and (11) show that the penalty matrix elastically ties the weighted average of the corner rotations to the in-plane rotation at the centre of the element, with no effect on the translational DOF. Choosing an appropriate value for the dimensionless coefficient $\gamma$ is crucial, as it involves a compromise between competing effects. If $\gamma$ is too large, any accuracy improvement resulting from the inclusion of rotations will be lost. On the other hand, if it is too low, the resulting matrices will become rank deficient. In Ref. [3], Macneal and coworkers suggested a value of ${10}^{-6}$ for $\gamma$.

2.2. Accuracy of Allman’s Elements

Allman’s elements have been employed in various studies, as documented in the literature [1,2,3,4,19]. However, the influence of the dimensionless coefficient $\gamma$, which controls the level of stabilisation provided to these elements, has received little attention. Moreover, to the best of our knowledge, the convergence rate of Allman’s elements in the energy norm has never been examined. To gain more insight into the accuracy of these elements, we conduct two series of numerical experiments. Material and geometric parameters will be given without units, as is customary in computational mechanics [18,21,26].

2.2.1. Short Cantilever Beam Test

The short cantilever beam test is employed to study the coarse mesh accuracy of Allman’s elements, their sensitivity to distortions and the effect of the dimensionless coefficient $\gamma$. The short cantilever beam test consists of applying a parabolic distribution of shear stresses at the tip of a beam while the left-hand side is fully clamped (Figure 2). Following Felippa [21], the problem is considered in a non-dimensional form, and we use Young’s modulus $E=2660$ and a Poisson’s ratio of $\nu =0.2$ and apply a vertical force $P=100$, so that the exact vertical displacement at the centre of the end-loaded cross section is $v_{analytical}=\frac{{L_{beam}}^{3}P}{3E{I}_Z}+\frac{(4+5\nu )\left(L_{beam}P\right)}{2E{H}_{beam}}=10$, where $L_{beam}$ and $H_{beam}$ are the beam’s dimensions, provided in Figure 2, and $I_z=\frac{L_{beam}{H}_{beam}^3}{12}$. A thickness of 1 and a plane stress state are considered.

The numerical experiment is run with Allman’s elements and the underlying higher-order finite elements depicted in Figure 1. For comparison purposes, we also benchmark the $2\times 2$ Gauss-integrated bilinear quadrilateral (Q4) and the constant strain triangle (T3). Regular and distorted meshes with one element through the thickness are employed (Figure 3).

The resulting vertical displacements at the centre of the end-loaded cross section are compared with the analytical solution value in

Table 1. Comparison of the computed vertical displacements at the tip of the short cantilever beam.

	Regular Mesh	Distorted Mesh	Number of DOFs
exact solution	10.	10.	-
T3	2.55017	1.4409	20
T3A $\gamma =0.$ T3A	7.65647	5.60141	30
$\gamma =10{.}^{-6}$
T6	9.85754	9.35355	54
Q4	6.85744	4.02667	20
Q4A $\gamma =0.$	9.45564	9.15315	30
Q4A $\gamma =10{.}^{-6}$	9.45511	9.15279	30
Q8	9.88888	9.91559	46

.

The following conclusions can be drawn from the numerical results: (i) the accuracy of elements with drilling DOFs lies in between those of the underlying higher-order elements and the linear elements, (ii) the use of drilling DOFs substantially increases the performance of classical linear elements, (iii) Allman’s elements are less sensitive to distortions than linear elements and (iv) the penalty stiffness does not affect the results when the recommended value $\gamma ={10}^{-6}$ is used.

These results are highly significant and indicate a positive outlook for the incorporation of Allman’s elements in the AFEM framework. As stressed in Section 3, the stabilised augmented finite element method (SAFEM) will rely on Allman’s elements in arbitrarily distorted configurations. Therefore, had the elements been highly sensitive to distortions, they would have been unsuitable for the method.

2.2.2. Convergence Rate in the Energy Norm

To the best of the authors’ knowledge, no previous studies have investigated the rate of convergence of Allman’s elements. As a result, the speed at which error measures decrease with mesh refinement when using these elements remains unknown. In this section, we focus on the error in the energy norm, as this measure is crucial to several fundamental properties of the FEM, such as Galerkin orthogonality (see Chap. 4.3 of Ref. [27]). The objective of this section is to investigate the convergence behaviour of Allman’s elements in the energy norm through numerical experiments. We will also compare the performances of both triangular and quadrilateral elements with the finite element method. Before presenting the numerical experiment, we introduce the error in the energy norm and some related quantities. Let $\mathbf{u}$ be the exact solution of a mechanical problem and $\widehat{\mathbf{u}}$ an approximate solution, and the error $\mathbf{e}$ is:

$$\mathbf{e}=\mathbf{u}-\widehat{\mathbf{u}}$$

(12)

Let $\mathsf{\Omega }$ be the volume occupied by the structure of interest, and the relative error in the energy norm is:

$$\left|\left|\mathbf{e}\right|\right|={\left(\frac{{\int }_{\mathsf{\Omega }}(\mathit{ϵ}-\widehat{\mathit{ϵ}}):\mathit{C}:(\mathit{ϵ}-\widehat{\mathit{ϵ}})d\mathsf{\Omega }}{{\int }_{\mathsf{\Omega }}\mathit{ϵ}:\mathit{C}:\mathit{ϵ}d\mathsf{\Omega }}\right)}^{\frac{1}{2}}$$

(13)

where $\mathit{ϵ}$ is the exact strain field, $\widehat{\mathit{ϵ}}$ is the finite element strain field and $\mathit{C}$ is the elasticity stiffness tensor. Galerkin orthogonality allows for simplifying the computation of the relative error in the energy norm:

$$\left|\left|\mathbf{e}\right|\right|={\left(\frac{{\int }_{\mathsf{\Omega }}\mathit{ϵ}:\mathit{C}:\mathit{ϵ}d\mathsf{\Omega }-{\int }_{\mathsf{\Omega }}\widehat{\mathit{ϵ}}:\mathit{C}:\widehat{\mathit{ϵ}}d\mathsf{\Omega }}{{\int }_{\mathsf{\Omega }}\mathit{ϵ}:\mathit{C}:\mathit{ϵ}d\mathsf{\Omega }}\right)}^{\frac{1}{2}}$$

(14)

To investigate the convergence rate of Allman’s elements, we perform a test involving a circular beam subjected to end shear in a state of plane stress. The geometry and loading of the problem are shown in Figure 4. The displacements on the boundaries are prescribed:

$$u_{\theta }(r,\theta ={90}^{\circ })=u_r(r=a,\theta ={90}^{\circ })=0\mathrm{and}{u}_r(r,\theta =0^{\circ })=u_0$$

(15)

Following Zienkiewicz (example 2.4 of Ref. [5]), the following geometric and material properties are considered: $a=5$, $b=10$, $E=\mathrm{10,000}$ and $\nu =0.25$. A displacement of magnitude $u_0=-0.01$ is imposed on the right end of the circular beam. The above values lead to the following closed-form strain energy [5]:

$${\int }_{\mathsf{\Omega }}\mathit{ϵ}:\mathit{C}:\mathit{ϵ}d\mathsf{\Omega }=\frac{1}{\pi }\left(log\left(2\right)-0.6\right)$$

(16)

The finite element solution to the problem is obtained using uniform meshes made of Allman’s elements and the underlying higher-order finite elements: the Q8 and the T6. For comparison, we also benchmarked the $2\times 2$ Gauss-integrated bilinear quadrilateral (Q4) as well as the constant strain triangle (T3) (cf. Figure 5). Elements with drilling DOFs are sometimes referred to as “higher-order” ones [21,28] because they derive from second-order elements. However, as shown in Figure 5, Allman’s elements converge with an error of $\mathcal{O}\left(h\right)$ in the energy norm, just like linear finite elements. This convergence rate could be expected because the tangential displacement along the edges of Allman’s elements is only linear (cf. Equation (2)). Thus, while Allman’s elements are more accurate than standard finite elements that share the same topology, they cannot be considered as "higher-order" elements. One can also notice that the penalty stiffness does not affect the results when the recommended value $\gamma ={10}^{-6}$ is employed. This value will be used consistently throughout this work.

3. Augmented Finite Elements with Drilling Degrees of Freedom

This self-contained section presents the derivation of the augmented finite element method (AFEM), the need for stabilising the method and the incorporation of drilling DOFs for achieving stabilisation. In Section 4, the resulting stabilised augmented finite element method (SAFEM) will be benchmarked against the AFEM. In this study, we only consider strong discontinuities because stabilisation is not required for modelling weak discontinuities [6,16,29].

3.1. Strong Form

The reference situation to be considered is schematised in Figure 6. Let $\mathsf{\Omega }$ be the domain occupied by a solid. A material point inside the domain is labelled as $\mathbf{x}\in \mathsf{\Omega }$. A cohesive strong discontinuity surface ${\mathsf{\Gamma }}_c={\mathsf{\Gamma }}_c^{+}\cup {\mathsf{\Gamma }}_c^{-}$ splits $\mathsf{\Omega }$ into two subdomains, ${\mathsf{\Omega }}^{+}$ and ${\mathsf{\Omega }}^{-}$. The prescribed external tractions ${\mathbf{t}}_{\mathbf{ext}}^{+}$ and ${\mathbf{t}}_{\mathbf{ext}}^{-}$ are applied on boundary ${\mathsf{\Gamma }}_t={\mathsf{\Gamma }}_t^{+}\cup {\mathsf{\Gamma }}_t^{-}$, whereas the displacements ${\overline{\mathbf{u}}}^{+}$ and ${\overline{\mathbf{u}}}^{-}$ are imposed on boundary ${\mathsf{\Gamma }}_u={\mathsf{\Gamma }}_u^{+}\cup {\mathsf{\Gamma }}_u^{-}$. The domains on both sides of the discontinuity are assumed to be elastic and homogeneous, yet ${\mathsf{\Omega }}^{+}$ and ${\mathsf{\Omega }}^{-}$ can be made of different materials. We further assume small strain and displacement conditions.

In the absence of body forces, the field equations governing the boundary value problem obey the following relations:

$$\begin{array}{cccccccc}\hfill & \nabla {\mathit{\sigma }}^{+}\left(\mathbf{x}\right)=\mathbf{0}\hfill & \hfill & \mathbf{x}\in {\mathsf{\Omega }}^{+}\hfill & \hfill & \nabla {\mathit{\sigma }}^{-}\left(\mathbf{x}\right)=\mathbf{0}\hfill & \hfill & \mathbf{x}\in {\mathsf{\Omega }}^{-}\hfill \end{array}$$

(17)

$$\begin{array}{cccccccc}\hfill & {\mathit{\sigma }}^{+}\left(\mathbf{x}\right)·{\mathbf{n}}^{+}\left(\mathbf{x}\right)={\mathbf{t}}_{\mathbf{ext}}^{+}\left(\mathbf{x}\right)\hfill & \hfill & \mathbf{x}\in {\mathsf{\Gamma }}_t^{+}\hfill & \hfill & {\mathit{\sigma }}^{-}\left(\mathbf{x}\right)·{\mathbf{n}}^{-}\left(\mathbf{x}\right)={\mathbf{t}}_{\mathbf{ext}}^{-}\left(\mathbf{x}\right)\hfill & \hfill & \mathbf{x}\in {\mathsf{\Gamma }}_t^{-}\hfill \end{array}$$

(18)

$$\begin{array}{cccccccc}\hfill & {\mathbf{u}}^{+}\left(\mathbf{x}\right)={\overline{\mathbf{u}}}^{+}\left(\mathbf{x}\right)\hfill & \hfill & \mathbf{x}\in {\mathsf{\Gamma }}_u^{+}\hfill & \hfill & {\mathbf{u}}^{-}\left(\mathbf{x}\right)={\overline{\mathbf{u}}}^{-}\left(\mathbf{x}\right)\hfill & \hfill & \mathbf{x}\in {\mathsf{\Gamma }}_u^{-}\hfill \end{array}$$

(19)

$$\begin{array}{cccccccc}\hfill & {\mathbf{t}}_{\mathbf{coh}}^{+}\left(\mathbf{x}\right)={\mathit{\sigma }}^{+}\left(\mathbf{x}\right)·{\mathbf{n}}^{+}\left(\mathbf{x}\right)\hfill & \hfill & \mathbf{x}\in {\mathsf{\Gamma }}_c^{+}\hfill & \hfill & {\mathbf{t}}_{\mathbf{coh}}^{-}\left(\mathbf{x}\right)={\mathit{\sigma }}^{-}\left(\mathbf{x}\right)·{\mathbf{n}}^{-}\left(\mathbf{x}\right)\hfill & \hfill & \mathbf{x}\in {\mathsf{\Gamma }}_c^{-}\hfill \end{array}$$

(20)

${\mathit{\sigma }}^{+}$ (resp. ${\mathit{\sigma }}^{-}$) is the stress field in ${\mathsf{\Omega }}^{+}$ (resp. ${\mathsf{\Omega }}^{-}$), ${\mathbf{t}}_{\mathbf{coh}}^{+}$ and ${\mathbf{t}}_{\mathbf{coh}}^{-}$ are the tractions along the discontinuity surface ${\mathsf{\Gamma }}_c={\mathsf{\Gamma }}_c^{+}\cup {\mathsf{\Gamma }}_c^{-}$ and ${\mathbf{n}}^{+}$ (resp. ${\mathbf{n}}^{-}$) is the outward-pointing normal of ${\mathsf{\Omega }}^{+}$ (resp. ${\mathsf{\Omega }}^{-}$). Tractions ${\mathbf{t}}_{\mathbf{coh}}^{+}$ and ${\mathbf{t}}_{\mathbf{coh}}^{-}$ are related to the cohesive strong discontinuity opening, denoted as $\Delta \mathbf{u}$ and given by:

$$\Delta \mathbf{u}\left(\mathbf{x}\right)={\mathbf{u}}^{+}\left(\mathbf{x}\right)-{\mathbf{u}}^{-}\left(\mathbf{x}\right)\mathbf{x}\in {\mathsf{\Gamma }}_c$$

(21)

Equilibrium imposes the following equalities:

$${\mathbf{t}}_{\mathbf{coh}}^{-}\left(\mathbf{x}\right)=-{\mathbf{t}}_{\mathbf{coh}}^{+}\left(\mathbf{x}\right)={\mathbf{t}}_{\mathbf{coh}}\left(\mathbf{x}\right)\mathbf{x}\in {\mathsf{\Gamma }}_c$$

(22)

The constitutive law and the strain–displacement equations for the two subdomains are:   

$$\begin{array}{cccccccc}\hfill & {\mathit{\sigma }}^{+}\left(\mathbf{x}\right)={\mathbf{C}}^{+}:{\mathit{ϵ}}^{+}\left(\mathbf{x}\right)\hfill & \hfill & \mathbf{x}\in {\mathsf{\Omega }}^{+}\hfill & \hfill & {\mathit{ϵ}}^{+}\left(\mathbf{x}\right)=\frac{1}{2}({\nabla }^T{\mathbf{u}}^{+}\left(\mathbf{x}\right)+\nabla {\mathbf{u}}^{+}\left(\mathbf{x}\right))\hfill & \hfill & \mathbf{x}\in {\mathsf{\Omega }}^{+}\hfill \end{array}$$

(23)

$$\begin{array}{cccccccc}\hfill & {\mathit{\sigma }}^{-}\left(\mathbf{x}\right)={\mathbf{C}}^{-}:{\mathit{ϵ}}^{-}\left(\mathbf{x}\right)\hfill & \hfill & \mathbf{x}\in {\mathsf{\Omega }}^{-}\hfill & \hfill & {\mathit{ϵ}}^{-}\left(\mathbf{x}\right)=\frac{1}{2}({\nabla }^T{\mathbf{u}}^{-}\left(\mathbf{x}\right)+\nabla {\mathbf{u}}^{-}\left(\mathbf{x}\right))\hfill & \hfill & \mathbf{x}\in {\mathsf{\Omega }}^{-}\hfill \end{array}$$

(24)

where ${\mathbf{C}}^{+}$ and ${\mathbf{C}}^{-}$ are the stiffness tensors of the subdomains ${\mathsf{\Omega }}^{+}$ and ${\mathsf{\Omega }}^{-}$, respectively. In this work, a weakly coupled and bilinear elasto-damage law governs the relationship between cohesive tractions ${\mathbf{t}}_{\mathbf{coh}}$ and openings $\mathbf{\Delta }\mathbf{u}$. It is important to recognise that the SAFEM and the AFEM can be employed with any shape of cohesive law (e.g., exponential or multi-linear). The traction separation law employed in this article is given by

$$\begin{array}{c}\hfill {\mathbf{t}}_{\mathbf{coh}}=\left[{\mathbf{D}}_{\mathbf{coh}/\mathbf{sec}}(d_n,d_t)\right]\mathbf{\Delta }\mathbf{u}\end{array}$$

(25)

where we introduced the scalar damage variables $d_n$ and $d_t$ associated with normal and tangential directions to the discontinuity surface. The matrix $\left[{\mathbf{D}}_{\mathbf{coh}/\mathbf{sec}}\right]$ is the secant stiffness matrix associated with the employed cohesive law. Because any cohesive law can be embedded within an augmented element, the detailed description of the employed traction–opening relationship is deferred to Appendix B. Recent work has demonstrated that the AFEM also allows interphases to be embedded within solid elements rather than interfaces [30].

Equations (17)–(25) define the problem strong form. The AFEM differs significantly from other strong discontinuity approaches, such as the numerous embedded finite element methods (EFEM) or the extended finite element method (XFEM) (see, e.g., Refs. [15,30,31,32,33,34,35,36,37]). Indeed, the cohesive internal tractions are treated similarly to the external tractions, and no discontinuous displacement or strain field is introduced. The above strong form actually shares many similarities with that of cohesive zone models, see, e.g., [38].

3.2. Discretised Weak Form

The spaces of the admissible continuous displacement fields ${\mathcal{U}}^{+}$ and ${\mathcal{U}}^{-}$ are defined by

$$\begin{array}{c}\hfill {\mathcal{U}}^{+}=\{{\mathbf{u}}^{+}\in H1:{\mathbf{u}}^{+}={\overline{\mathbf{u}}}^{+}\mathrm{on}{\mathsf{\Gamma }}_u^{+}\}{\mathcal{U}}^{-}=\{{\mathbf{u}}^{-}\in H1:{\mathbf{u}}^{-}={\overline{\mathbf{u}}}^{-}\mathrm{on}{\mathsf{\Gamma }}_u^{-}\}\end{array}$$

(26)

where $H1$ is the space of functions with square-integrable derivatives, i.e., the Sobolev space of degree 1 defined by:

$$\begin{array}{c}\hfill H1=\{\mathbf{w}:\mathbf{w}\in L2,D\mathbf{w}\in L2\}\mathrm{with}L2=\{\mathbf{w}:\int {\mathbf{w}}^{2}dx<\infty \}\end{array}$$

(27)

where $\mathbf{w}$ is a function and $D\mathbf{w}$ denotes its partial derivative.

Let us introduce the spaces of continuous test functions ${\mathcal{V}}^{+}$ and ${\mathcal{V}}^{-}$ defined by

$$\begin{array}{c}\hfill {\mathcal{V}}^{+}=\{{\mathbf{v}}^{+}\in H1:{\mathbf{v}}^{+}=\mathbf{0}\mathrm{on}{\mathsf{\Gamma }}_u^{+}\}{\mathcal{V}}^{-}=\{{\mathbf{v}}^{-}\in H1:{\mathbf{v}}^{-}=\mathbf{0}\mathrm{on}{\mathsf{\Gamma }}_u^{-}\}\end{array}$$

(28)

The weak forms of the equilibrium equations in the two subdomains ${\mathsf{\Omega }}^{+}$ and ${\mathsf{\Omega }}^{-}$ are stated as follows: find ${\mathbf{u}}^{+}\in {\mathcal{U}}^{+}$ and ${\mathbf{u}}^{-}\in {\mathcal{U}}^{-}$ such that

$$\begin{array}{c}\hfill {\int }_{{\mathsf{\Omega }}^{+}}{\mathit{\sigma }}^{+}\left({\mathbf{u}}^{+}\left(\mathbf{x}\right)\right):{\mathit{ϵ}}^{+}\left({\mathbf{v}}^{+}\left(\mathbf{x}\right)\right)d\mathsf{\Omega }={\int }_{{\mathsf{\Gamma }}_t^{+}}{\mathbf{t}}_{\mathbf{ext}}^{+}\left(\mathbf{x}\right)·{\mathbf{v}}^{+}\left(\mathbf{x}\right)d\mathsf{\Gamma }-{\int }_{{\mathsf{\Gamma }}_c^{+}}{\mathbf{t}}_{\mathbf{coh}}\left(\mathbf{x}\right)·{\mathbf{v}}^{+}\left(\mathbf{x}\right)d\mathsf{\Gamma }\forall \mathbf{v}\in {\mathcal{V}}^{+}\end{array}$$

(29)

$$\begin{array}{c}\hfill {\int }_{{\mathsf{\Omega }}^{-}}{\mathit{\sigma }}^{-}\left({\mathbf{u}}^{-}\left(\mathbf{x}\right)\right):{\mathit{ϵ}}^{-}\left({\mathbf{v}}^{-}\left(\mathbf{x}\right)\right)d\mathsf{\Omega }={\int }_{{\mathsf{\Gamma }}_t^{-}}{\mathbf{t}}_{\mathbf{ext}}^{-}\left(\mathbf{x}\right)·{\mathbf{v}}^{-}\left(\mathbf{x}\right)d\mathsf{\Gamma }+{\int }_{{\mathsf{\Gamma }}_c^{-}}{\mathbf{t}}_{\mathbf{coh}}\left(\mathbf{x}\right)·{\mathbf{v}}^{-}\left(\mathbf{x}\right)d\mathsf{\Gamma }\forall \mathbf{v}\in {\mathcal{V}}^{-}\end{array}$$

(30)

The above equations are referred to as the principle of virtual work (see, e.g., p. 78 of Ref. [39]). The left-hand sides of Equations (29) and (30) are the internal virtual work, and the right-hand sides are the virtual work carried out by the external forces and the tractions along the strong discontinuity surface. The application of the (Bubnov–)Galerkin method to Equations (29) and (30) leads to the following matrix equations (see, e.g., Chap. 2.8 [39]):   

$$\begin{array}{ccc}\hfill \left[L^{+}\right]\left\{d^{+}\right\}=\left\{\begin{array}{c}fex{t}^{+}\\ -fco{h}^{+}\end{array}\right\}& \hfill & \hfill \left[L^{-}\right]\left\{d^{-}\right\}=\left\{\begin{array}{c}fex{t}^{-}\\ -fco{h}^{-}\end{array}\right\}\end{array}$$

(31)

where vector $\left\{d^{+}\right\}$ (resp. $\left\{d^{-}\right\}$) contains the DOF in ${\mathsf{\Omega }}^{+}$ (resp. ${\mathsf{\Omega }}^{-}$). $\left[L^{+}\right]$ and $\left[L^{-}\right]$ are the finite element stiffness matrices of subdomains ${\mathsf{\Omega }}^{+}$ and ${\mathsf{\Omega }}^{-}$, respectively. If, for example, ${\mathsf{\Omega }}^{+}$ is a triangular subdomain, then $\left[L^{+}\right]$ would correspond to the stiffness matrix of a T3 element in the AFEM formulation or the stiffness matrix of a T3A element (as defined by Equation (8)) in the SAFEM formulation. The external forces, defined later in this section, appear on the right-hand side of Equation (31). The displacement field in ${\mathsf{\Omega }}^{+}$ (resp. ${\mathsf{\Omega }}^{-}$) is expressed thanks to the shape function matrix $\left[N^{+}\right]$ (resp. $\left[N^{-}\right]$):

$$\begin{array}{c}\hfill {\mathbf{u}}^{+}\left(\mathbf{x}\right)=\left[N^{+}\left(\mathbf{x}\right)\right]\left\{d^{+}\right\}\mathbf{x}\in {\mathsf{\Omega }}^{+}{\mathbf{u}}^{-}\left(\mathbf{x}\right)=\left[N^{-}\left(\mathbf{x}\right)\right]\left\{d^{-}\right\}\mathbf{x}\in {\mathsf{\Omega }}^{-}\end{array}$$

(32)

The interpolations presented in Equation (32) are only valid for an element crossed by a strong discontinuity, i.e., an augmented element. Crack-free elements are standard finite elements that use a single shape function matrix $\left[N\right]$.

To better illustrate the formulation, we consider a triangular domain $\mathsf{\Omega }$ divided into a quadrilateral subdomain (${\mathsf{\Omega }}^{-}$) and a triangular subdomain (${\mathsf{\Omega }}^{+}$), as depicted in Figure 7. With this configuration, the shape function matrix $\left[N^{+}\right]$ (resp. $\left[N^{-}\right]$) corresponds to a T3A (resp. Q4A) element with the SAFEM and a T3 (resp. Q4) element with the AFEM. The DOF vectors in Equation (32) are partitioned between those associated with the strong discontinuity surface, denoted as $\left\{dco{h}^{+}\right\}$ and $\left\{dco{h}^{-}\right\}$, and those associated with the surface of external tractions, $\left\{dex{t}^{+}\right\}$ and $\left\{dex{t}^{-}\right\}$, as follows:

$$\begin{array}{ccc}\hfill \left\{d^{+}\right\}=\left\{\begin{array}{c}dex{t}^{+}\\ dco{h}^{+}\end{array}\right\}& \hfill & \hfill \left\{d^{-}\right\}=\left\{\begin{array}{c}dex{t}^{-}\\ dco{h}^{-}\end{array}\right\}\end{array}$$

(33)

The DOF vector $\left\{dco{h}^{+}\right\}$ (resp. $\left\{dco{h}^{-}\right\}$) allows for expressing the strong discontinuity surface opening thanks to the shape functions $N^{+}$ (resp. $N^{-}$). The opening–displacement matrix of the cohesive surface $\left[B^{coh}\right]$ writes [38]

$$\mathbf{\Delta }\mathbf{u}\left(\mathbf{x}\right)=\left[\begin{array}{cc}-N^{-}\left(\mathbf{x}\right)& N^{+}\left(\mathbf{x}\right)\end{array}\right]\left\{\begin{array}{c}dco{h}^{-}\\ dco{h}^{+}\end{array}\right\}=\left[B^{coh}\left(\mathbf{x}\right)\right]\left\{\begin{array}{c}dco{h}^{-}\\ dco{h}^{+}\end{array}\right\}\mathbf{x}\in {\mathsf{\Gamma }}_c$$

(34)

Matrices $\left[L^{+}\right]$ and $\left[L^{-}\right]$ in Equation (31) are the stiffness matrices of the two linear elastic subdomains and are given by

$$\begin{array}{c}\hfill \left[L^{+}\right]=\left[\begin{array}{cc}{L}_{11}^{+}& L_{12}^{+}\\ L_{21}^{+}& L_{22}^{+}\end{array}\right]={\int }_{{\mathsf{\Omega }}^{+}}{\left[B^{+}\left(\mathbf{x}\right)\right]}^t\left[{\mathbf{C}}^{+}\right]\left[B^{+}\left(\mathbf{x}\right)\right]d\mathsf{\Omega }\end{array}$$

(35)

$$\left[L^{-}\right]=\left[\begin{array}{cc}{L}_{11}^{-}& L_{12}^{-}\\ L_{21}^{-}& L_{22}^{-}\end{array}\right]={\int }_{{\mathsf{\Omega }}^{-}}{\left[B^{-}\left(\mathbf{x}\right)\right]}^t\left[{\mathbf{C}}^{-}\right]\left[B^{-}\left(\mathbf{x}\right)\right]d\mathsf{\Omega }$$

(36)

where we introduced the strain–displacement matrix $\left[B^{+}\right]$ (resp. $\left[B^{-}\right]$), which contains the derivatives of the shape functions in ${\mathsf{\Omega }}^{+}$ (resp. ${\mathsf{\Omega }}^{-}$). Vectors $\left\{fex{t}^{+}\right\}$ and $\left\{fex{t}^{-}\right\}$ in Equation (31) are the external force vectors induced by the external tractions ${\mathbf{t}}_{\mathbf{ext}}$:

$$\begin{array}{ccc}\hfill \left\{fex{t}^{+}\right\}={\int }_{{\mathsf{\Gamma }}_t^{+}}{\left[N^{+}\left(\mathbf{x}\right)\right]}^t{\mathbf{t}}_{\mathbf{ext}}^{+}\left(\mathbf{x}\right)d\mathsf{\Gamma }& \hfill & \hfill \left\{fex{t}^{-}\right\}={\int }_{{\mathsf{\Gamma }}_t^{-}}{\left[N^{-}\left(\mathbf{x}\right)\right]}^t{\mathbf{t}}_{\mathbf{ext}}^{-}\left(\mathbf{x}\right)d\mathsf{\Gamma }\end{array}$$

(37)

Vectors $\left\{fco{h}^{+}\right\}$ and $\left\{fco{h}^{-}\right\}$ in Equation (31) are the equivalent force vectors induced by the tractions exerted on the strong discontinuity and are given by   

$$\begin{array}{c}\hfill \left\{\begin{array}{c}fco{h}^{-}\\ fco{h}^{+}\end{array}\right\}=\left[K^{coh/sec}\right]\left\{\begin{array}{c}dco{h}^{-}\\ dco{h}^{+}\end{array}\right\}\end{array}$$

(38)

where $\left[K^{coh/sec}\right]$ is the secant finite element matrix associated with the cohesive law Equation (25):   

$$\left[K^{coh/sec}\right]=\left[\begin{array}{cc}{K}_{11}^{coh/sec}& K_{12}^{coh/sec}\\ K_{21}^{coh/sec}& K_{22}^{coh/sec}\end{array}\right]={\int }_{{\mathsf{\Gamma }}_c}{\left[B^{coh}\left(\mathbf{x}\right)\right]}^t\left[{\mathbf{D}}^{\mathbf{coh}/\mathbf{sec}}\left(\mathbf{x}\right)\right]\left[B^{coh}\left(\mathbf{x}\right)\right]d\mathsf{\Gamma }$$

(39)

In this study, we utilise classical cohesive elements equipped with displacement degrees of freedom only, as discussed in [38]. Although cohesive elements with drilling degrees of freedom do exist and have been reported in [28], their performance has not been thoroughly studied in the literature. Therefore, we have chosen to stick with the classical cohesive element formulation for this investigation. The use of cohesive elements with drilling degrees of freedom will be considered for future studies.

The substitution of Equations (33), (35), (36) and (38) into (31) leads to

$$\begin{array}{cc}\hfill \left[\begin{array}{cc}{L}_{11}^{-}& 0\\ 0& L_{11}^{+}\end{array}\right]\left\{\begin{array}{c}dex{t}^{-}\\ dex{t}^{+}\end{array}\right\}+\left[\begin{array}{cc}{L}_{12}^{-}& 0\\ 0& L_{12}^{+}\end{array}\right]\left\{\begin{array}{c}dco{h}^{-}\\ dco{h}^{+}\end{array}\right\}-\left\{\begin{array}{c}fex{t}^{-}\\ fex{t}^{+}\end{array}\right\}=& \left\{\begin{array}{c}0\\ 0\end{array}\right\}\hfill \end{array}$$

(40)

$$\begin{array}{cc}\hfill \left[\begin{array}{cc}{L}_{21}^{-}& 0\\ 0& L_{21}^{+}\end{array}\right]\left\{\begin{array}{c}dex{t}^{-}\\ dex{t}^{+}\end{array}\right\}+\left[\begin{array}{cc}{L}_{22}^{-}+K_{11}^{coh/sec}& K_{12}^{coh/sec}\\ K_{21}^{coh/sec}& L_{22}^{+}+K_{22}^{coh/sec}\end{array}\right]\left\{\begin{array}{c}dco{h}^{-}\\ dco{h}^{+}\end{array}\right\}=& \left\{\begin{array}{c}0\\ 0\end{array}\right\}\hfill \end{array}$$

(41)

Equations (40) and (41) represent the AFEM discretised equilibrium equations. The problem to solve can be summed up as follows: find the degrees-of-freedom vectors $\left\{dco{h}^{-}\right\}$, $\left\{dco{h}^{+}\right\}$, $\left\{dex{t}^{+}\right\}$ and $\left\{dex{t}^{-}\right\}$ knowing $\left\{fex{t}^{+}\right\}$ and $\left\{fex{t}^{-}\right\}$ such that Equations (40) and (41) are satisfied.

3.3. Condensed Discretised Equilibrium Equations and Nested Global/Local Solving Procedure

One critical aspect of augmented finite element methods lies in the nested global/local procedure adopted to solve the coupled equilibrium Equations (40) and (41). Indeed, this procedure allows for the condensation of DOFs related to the strong discontinuity surface, namely, $\left\{dco{h}^{+}\right\}$ and $\left\{dco{h}^{-}\right\}$. Consequently, cracks can be embedded in elements without requiring additional global DOFs. The condensation process is straightforward: assuming that Equations (40) and (41) have been solved (i.e., their right-hand sides are null), the substitution of Equation (41) into Equation (40) yields:   

$$\begin{array}{c}\underset{=\left[K_{AFEM}^{sec}\right]}{\underbrace{\left(\left[\begin{array}{cc}{L}_{11}^{-}& 0\\ 0& L_{11}^{+}\end{array}\right]-\left[\begin{array}{cc}{L}_{12}^{-}& 0\\ 0& L_{12}^{+}\end{array}\right]{\left[\begin{array}{cc}{L}_{22}^{-}+K_{11}^{coh/sec}& K_{12}^{coh/sec}\\ K_{21}^{coh/sec}& L_{22}^{+}+K_{22}^{coh/sec}\end{array}\right]}^{-1}\left[\begin{array}{cc}{L}_{21}^{-}& 0\\ 0& L_{21}^{+}\end{array}\right]\right)}}\left\{\begin{array}{c}dex{t}^{-}\\ dex{t}^{+}\end{array}\right\}\hfill \\ \hfill -\left\{\begin{array}{c}fex{t}^{-}\\ fex{t}^{+}\end{array}\right\}=\left\{\begin{array}{c}0\\ 0\end{array}\right\}\end{array}$$

(42)

The (secant) stiffness matrix $\left[K_{AFEM}^{sec}\right]$ of an element is dependent on the state of the cohesive crack, but it does not explicitly reference the associated degrees of freedom $\left\{dco{h}^{\pm }\right\}$. Thus, it is commonly referred to as a condensed stiffness matrix, and Equation (42) is known as the condensed discretised equilibrium equation. The matrix-vector products $\left[K^{coh/sec}\right].{\left\{\begin{array}{c}dco{h}^{-}dco{h}^{+}\end{array}\right\}}^t$ (Equation (41)) and $\left[K_{AFEM}^{sec}\right].{\left\{\begin{array}{c}dex{t}^{-}dex{t}^{+}\end{array}\right\}}^t$ (Equation (42)) are similar to the internal force vector $\int {\left[B\right]}^t\sigma dV$ computed with classical finite elements. In the AFEM framework, the global DOFs associated with cracked and uncracked elements are identical. Therefore, low-level aspects of finite element codes, such as DOF numbering and matrix assembly, are not affected by the use of augmented elements.

Let us now shift our focus to the incremental-iterative strategy adopted to solve Equation (42). In this regard, we introduce indices n, i and j related to increments, local iterations and global iterations, respectively. The problem aims to determine ${\left\{dex{t}^{\pm }\right\}}^{n+1}$, given ${\left\{fex{t}^{\pm }\right\}}^{n+1}$ and ${\left\{dco{h}^{\pm }\right\}}^n$ such that Equation (42) holds. The solution to the problem follows a two-stage approach, referred to as the global and local stages, and is elaborated in detail in Algorithm 1. Firstly, we compute a trial value of the external degrees of freedom ${\left\{dex{t}^{\pm }\right\}}^j$ associated with the known load vectors ${\left\{fex{t}^{\pm }\right\}}^{n+1}$, under the assumption that crack openings remain constant, i.e., $\{\Delta d_{coh}^{\pm }\}=\left\{0\right\}$. The condensed tangent stiffness matrix of augmented elements is introduced for that purpose and is given by

$$\left[K_{AFEM}^{tan}\right]=\left[\begin{array}{cc}{L}_{11}^{-}& 0\\ 0& L_{11}^{+}\end{array}\right]-\left[\begin{array}{cc}{L}_{12}^{-}& 0\\ 0& L_{12}^{+}\end{array}\right]{\left(\left[\begin{array}{cc}{L}_{22}^{-}& 0\\ 0& L_{22}^{+}\end{array}\right]+\left[K^{coh/tan}\right]\right)}^{-1}\left[\begin{array}{cc}{L}_{21}^{-}& 0\\ 0& L_{21}^{+}\end{array}\right]$$

(43)

where $\left[K^{coh/tan}\right]$ is the tangent stiffness matrix associated with the cohesive interface, as given in Appendix B. The computation of the trial values of ${\left\{dex{t}^{\pm }\right\}}^j$ does not require the assembly of quantities related to the cohesive interface: vectors ${\left\{dco{h}^{\pm }\right\}}^i$ are treated as internal variables. This one-step operation is the global stage.

The local stage consists of finding ${\left\{dco{h}^{\pm }\right\}}^i$ such that Equation (41) is satisfied. The external degrees of freedom remain constant during this iterative stage (i.e., ${\{\Delta d_{ext}^{\pm }\}}^j=\left\{0\right\}$). Two markedly different strategies can be adopted during the local stage. The one mostly employed by Yang and coworkers considers cracks as being local to the elements, i.e., cohesive surface-related quantities ($\left\{dco{h}^{\pm }\right\}$, $\left[K_{AFEM}^{tan}\right]$, $\left[K_{AFEM}^{sec}\right]$ and $\left\{fco{h}^{\pm }\right\}$) are not assembled but are instead treated at the element level. Algorithm 1 then becomes similar to a nonlinear constitutive law with the difference being that the iterations are performed at the element level and not at the integration point level. These key aspects enable the efficient implementation of augmented elements in any industrial code equipped with user-defined elements because Equation (41) is solved independently within each element. This strategy nevertheless introduces a significant drawback: the interelement compatibility of the displacement field is lost and, to the best of the authors’ knowledge, the equivalence between the problem’s strong and weak forms remains to be demonstrated. The second possible strategy maintains interelement compatibility of the displacement field. To do so, one follows Algorithm 1 and assembles the cohesive surface-related quantities before solving Equation (41). Because this approach requires DOF management and matrix assembly, it cannot be easily implemented in (closed-source) industrial software. Yang and coworkers coined this second strategy the conforming augmented finite element method (C-AFEM) [40,41].

Algorithm 1 Calculate ${\left\{dex{t}^{\pm }\right\}}^{n+1}$ given ${\left\{fex{t}^{\pm }\right\}}^{n+1}$, ${\left\{dex{t}^{\pm }\right\}}^n$ and ${\left\{dco{h}^{\pm }\right\}}^n$

${\left\{\begin{array}{c}fex{t}^{-}\hfill \\ {\mathit{fext}}^{+}\hfill \end{array}\right\}}^j\leftarrow {\left\{\begin{array}{c}fex{t}^{-}\\ fex{t}^{+}\end{array}\right\}}^{n+1}$       ${\left\{\begin{array}{c}dex{t}^{-}\\ dex{t}^{+}\end{array}\right\}}^j\leftarrow {\left\{\begin{array}{c}dex{t}^{-}\\ dex{t}^{+}\end{array}\right\}}^n$       ${\left\{\begin{array}{c}dco{h}^{-}\\ dco{h}^{+}\end{array}\right\}}^i\leftarrow {\left\{\begin{array}{c}dco{h}^{-}\\ dco{h}^{+}\end{array}\right\}}^n$ while${∥{\left\{\begin{array}{c}re{s}_{global}\end{array}\right\}}^j∥}_{L2}={∥{\left\{\begin{array}{c}fex{t}^{-}\\ fex{t}^{+}\end{array}\right\}}^j-\left[K_{AFEM}^{sec}\left({\left\{\begin{array}{c}dco{h}^{\pm }\end{array}\right\}}^i\right)\right]{\left\{\begin{array}{c}dex{t}^{-}\\ dex{t}^{+}\end{array}\right\}}^j∥}_{L2}>{ϵ}_{global}$  do     ${\left\{\begin{array}{c}\Delta dex{t}^{-}\\ \Delta dex{t}^{+}\end{array}\right\}}^j\leftarrow {\left[K_{AFEM}^{tan}\left({\left\{\begin{array}{c}dco{h}^{\pm }\end{array}\right\}}^i\right)\right]}^{-1}{\left\{\begin{array}{c}re{s}_{global}\end{array}\right\}}^j$     ${\left\{\begin{array}{c}dex{t}^{-}\\ dex{t}^{+}\end{array}\right\}}^j\leftarrow {\left\{\begin{array}{c}dex{t}^{-}\\ dex{t}^{+}\end{array}\right\}}^j+{\left\{\begin{array}{c}\Delta dex{t}^{-}\\ \Delta dex{t}^{+}\end{array}\right\}}^j$     while ${∥{\left\{\begin{array}{c}re{s}_{local}\end{array}\right\}}^i∥}_{L2}={∥\left[\begin{array}{cc}{L}_{21}^{-}& 0\\ 0& L_{21}^{+}\end{array}\right]{\left\{\begin{array}{c}dex{t}^{-}\\ dex{t}^{+}\end{array}\right\}}^j+\left(\left[\begin{array}{cc}{L}_{22}^{-}& 0\\ 0& L_{22}^{+}\end{array}\right]+\left[\begin{array}{c}{K}^{coh/sec}\left({\left\{\begin{array}{c}dco{h}^{\pm }\end{array}\right\}}^i\right)\end{array}\right]\right){\left\{\begin{array}{c}dco{h}^{-}\\ dco{h}^{+}\end{array}\right\}}^i∥}_{L2}>$ ${ϵ}_{local}$ do         ${\left\{\begin{array}{c}\Delta dco{h}^{-}\\ \Delta dco{h}^{+}\end{array}\right\}}^i\leftarrow {\left(\left[\begin{array}{cc}{L}_{22}^{-}& 0\\ 0& L_{22}^{+}\end{array}\right]+\left[K^{coh/tan}\left({\left\{\begin{array}{c}dco{h}^{\pm }\end{array}\right\}}^i\right)\right]\right)}^{-1}{\left\{\begin{array}{c}re{s}_{local}\end{array}\right\}}^i$         ${\left\{\begin{array}{c}dco{h}^{-}\\ dco{h}^{+}\end{array}\right\}}^i\leftarrow {\left\{\begin{array}{c}dco{h}^{-}\\ dco{h}^{+}\end{array}\right\}}^i+{\left\{\begin{array}{c}\Delta dco{h}^{-}\\ \Delta dco{h}^{+}\end{array}\right\}}^i$     end while end while ${\left\{\begin{array}{c}dex{t}^{-}\\ dex{t}^{+}\end{array}\right\}}^{n+1}\leftarrow {\left\{\begin{array}{c}dex{t}^{-}\\ dex{t}^{+}\end{array}\right\}}^j$       ${\left\{\begin{array}{c}dco{h}^{-}\\ dco{h}^{+}\end{array}\right\}}^{n+1}\leftarrow {\left\{\begin{array}{c}dco{h}^{-}\\ dco{h}^{+}\end{array}\right\}}^i$

The second (conforming) strategy is theoretically more appealing because it ensures weak- and strong-form equivalence. Still, we will focus on the first (nonconforming) approach to authorise the implementation of the method in industrial software. Furthermore, the nonconforming strategy has been observed to outperform well-established methods such as the FEM or the XFEM in some benchmarks [13,16].

3.4. Treatment of Stiffness Matrices Singularities

The inverse of some stiffness submatrices appears in the AFEM discretised equilibrium equation (Equation (42)). The invertibility of these submatrices depends on two factors: (i) the parent element used, such as the constant strain triangular element (T3) or the bilinear quadrilateral element (Q4), and (ii) the crack’s location inside the element. Essongue and colleagues have shown that using the T3 parent elements leads to the singularity of these submatrices when cohesive tractions are absent [16]. To help clarify, assume that the cohesive discontinuity is completely damaged and opened, meaning that the values for $d_n$ and $d_t$ in Equation (25) are both equal to 1, and ${\delta }_n$ in Equation (A3) is greater than 0. This causes the cohesive finite element matrix, defined by Equation (39), to become a null matrix, $\left[K^{coh/sec}\right]=\left[0\right]$. To use Equation (42), it is necessary for submatrices $\left[L_{22}^{+}\right]$ and $\left[L_{22}^{-}\right]$ to be invertible. However, when adopting the T3 parent elements, $\left[L_{22}^{+}\right]$ becomes singular. To understand the invertibility of submatrices $\left[L_{22}^{\pm }\right]$, it is helpful to consider their mechanical interpretation. $\left[L_{22}^{+}\right]$ (and $\left[L_{22}^{-}\right]$) is the stiffness matrix for subdomain ${\mathsf{\Omega }}^{+}$ (and ${\mathsf{\Omega }}^{-}$) with a fully constrained DOF $\left\{dex{t}^{+}\right\}$ (and $\left\{dex{t}^{-}\right\}$) [16]. This is illustrated in Figure 7 for both the triangular augmented element (AT3) and the stabilised augmented triangular element (AT3A).

As can be seen in Figure 7, the upper triangular region of the AT3 is free to rotate, resulting in a singular $\left[L_{22}^{+}\right]$. Drilling DOFs are used to prevent this rigid-body rotation and ensure full rank of $\left[L_{22}^{+}\right]$. The effectiveness of this approach is demonstrated in the numerical experiments presented in Section 4.1.1.

Although the originators of the AFEM noted that “the nonzero stiffness in cohesive law helps to regularize the matrix condition” [6], they did not provide guidance on how to handle a fully softened cohesive zone. To our knowledge, the only published strategy for dealing with this situation is to compute the pseudoinverse of $\left[L_{22}^{+}\right]$, as proposed in our previous work [16]. Our numerical experiments have demonstrated that this method converged with mesh refinement and competed with the FEM or XFEM. Therefore, we use this approach with AT3 elements. The AT3A does not require pseudoinverse computations because the introduction of drilling DOFs suppresses the occurrence of singular stiffness submatrices (as shown in Section 4.1.1).

The AFEM is not the only embedded crack method that produces singular matrices. Several embedded finite element methods (EFEMs) also share this drawback. However, stabilisation strategies used in EFEMs cannot be directly applied to the AFEM. This is because EFEMs use prescribed crack kinematics (such as a constant or linear crack opening) for stabilisation [42,43,44], unlike the AFEM, which does not compute a crack opening directly. In the AFEM, the subdomains’ displacements are first solved, and then the crack opening is computed (using Equations (34),(40) and (41)). Moreover, the two subdomains of an augmented element have different kinematics (Equation (32)), a feature not shared by most EFEMs [31,32,43]. Therefore, stabilisation procedures that have been developed for EFEMs cannot be directly applied to the AFEM due to their differences.

3.5. Working with Various Element Shapes

While the upcoming numerical experiments will primarily concentrate on 2D triangular elements, it is important to note that the SAFEM is applicable in 3D scenarios and with various element shapes. The only prerequisite is the ability to construct the stiffness matrix associated with subdomains ${\mathsf{\Omega }}^{+}$ and ${\mathsf{\Omega }}^{-}$ (Equation (42)), which might involve some additional effort. Suppose, for instance, that the initial crack-free elements are Allman quadrilaterals (Q4A). If the Q4A is augmented due to the nucleation of a crack that splits the element into triangular and pentagonal subdomains, the triangular domain is to be modelled with the T3A presented in Section 2. One then needs to derive a pentagonal element with drilling DOFs, which would be termed P5A, to model the pentagonal subdomain. To the best of the authors’ knowledge, this element is absent from the literature. Nevertheless, the methodology to derive an Allman element from an underlying high-order element is straightforward (cf. Section 2). Hence, in the present example, one could use the high-order pentagonals proposed in Ref. [45] or [46] to compute the stiffness matrix of the P5A. Importantly, the SAFEM algorithm remains unmodified.

4. Numerical Experiments

This section compares the performances of the AFEM and the SAFEM. To allow for detailed observations, the numerical experiments are divided into two sections: traction-free strong discontinuities are investigated in Section 4.1 while cohesive discontinuities are considered in Section 4.2. To prompt the need for stabilisation, the computations will be performed exclusively on 2D triangular meshes. The material and geometric parameters will be given without units as is customary in computational mechanics [18,21,26].

4.1. Traction-Free Strong Discontinuities

While the AFEM was primarily designed to simulate the propagation of cohesive cracks, studying its ability to represent traction-free discontinuities is also of paramount importance. Indeed, traction-free situations put to light some AFEM deficiencies that (i) may remain undetected when cohesive tractions are large enough but (ii) could result in the termination of a simulation when interfaces are fully damaged. The next two sections investigate the singularities of the matrices arising in the AFEM formulation, as well as its convergence rate in traction-free situations.

4.1.1. The Stabilising Effects of Drilling Degrees of Freedom

As highlighted in Section 3.4, the stiffness matrices of augmented elements may become singular due to an internal rigid-body rotation (Figure 7). To prevent this, drilling DOFs were introduced into the AFEM framework giving rise to the SAFEM. However, Allman’s drilling DOFs are not true rotations (as explained in Section 2.1). Therefore, it is reasonable to question whether they effectively suppress the internal rigid-body rotation that occurs in the AFEM framework. To investigate this, we conducted numerical experiments where we computed the lowest eigenvalue of the submatrix $\left[L_{{22}^{+}}\right]$ of the AT3 and AT3A elements, denoted as ${\lambda }_{AT3}$ and ${\lambda }_{AT3A}$, for a large number of geometries. These geometries were generated by varying the coordinates $\{x_1,y_1\}$ of the first node of the triangular domain illustrated in Figure 8.

The following isotropic material properties are (arbitrarily) chosen: Young’s modulus is $E={10}^7$ and Poisson’s ratio is $\nu =0.3$. The values of ${\lambda }_{AT3}$ and ${\lambda }_{AT3A}$ computed for ${10}^6$ different sets of coordinates $\{x_1,y_1\}$ are plotted in Figure 9.

${\lambda }_{AT3}$ oscillates around 0, illustrating the singularity of $\left[L_{{22}^{+}}\right]$ that occurs when AT3 elements are used. Conversely, ${\lambda }_{AT3A}$ is always positive and evolves smoothly with the coordinates of Node 1, even in the limited case of a triangular domain collapsed into a line (e.g., $\{x_1,y_1\}=\{1,0\}$). This is a significant outcome that demonstrates the SAFEM’s ability to produce non-singular element stiffness submatrices $\left[L_{{22}^{+}}\right]$ as required by the condensed discretised equilibrium equation (Equation (42)). Consequently, the computation of the pseudoinverse of $\left[L_{{22}^{+}}\right]$, which is necessary in the AFEM, is no longer needed. This simplifies and improves the efficiency of augmented elements because closed-form solutions for the inverse of symmetric matrices can be used advantageously. The next section investigates the impact of drilling DOF incorporation on the accuracy of augmented elements, which is essential to further validate the SAFEM’s performance.

4.1.2. Convergence of the AFEM and the SAFEM in the Energy and L2 Norms

The goal of this section is to compare the AFEM and the SAFEM convergence behaviour in the energy and L2 norms. To achieve this, we will perform numerical experiments using known analytical solutions. We conduct a test case encountered in the XFEM literature [26,47]: an infinite plate cut by a traction-free and horizontal crack is loaded by a remote stress field and one successively considers the Mode I and Mode II loadings of the crack. Plane strain is assumed and the material is homogeneous and isotropic with Young’s modulus E and Poisson’s ratio $\nu$. The analytical expression of the displacement field under Mode I loading, denoted as ${\mathit{u}}_{\mathit{I}}=\{u_{Ix},u_{Iy}\}$, is given by [48]:

$$\begin{array}{cc}\hfill u_{Ix}=& \frac{K_I}{E}\sqrt{\frac{r}{2\pi }}(1+\nu )cos\left(\frac{\theta }{2}\right)\left(3-4\nu -cos\left(\theta \right)\right)\hfill \end{array}$$

(44)

$$\begin{array}{cc}\hfill u_{Iy}=& \frac{K_I}{E}\sqrt{\frac{r}{2\pi }}(1+\nu )sin\left(\frac{\theta }{2}\right)\left(3-4\nu -cos\left(\theta \right)\right)\hfill \end{array}$$

(45)

where $K_I$ is the Mode I stress intensity factor and $(r,\theta )$ are the polar coordinates associated with a reference frame centred at the crack tip as depicted in Figure 10.

The closed-form solution of the displacement field under Mode II loading, denoted as ${\mathit{u}}_{\mathit{II}}$, is [48]:

$$\begin{array}{cc}\hfill u_{IIx}=& \frac{K_{II}}{E}\sqrt{\frac{r}{2\pi }}(1+\nu )sin\left(\frac{\theta }{2}\right)\left(5-4\nu +cos\left(\theta \right)\right)\hfill \end{array}$$

(46)

$$\begin{array}{cc}\hfill u_{IIy}=& \frac{K_{II}}{E}\sqrt{\frac{r}{2\pi }}(1+\nu )cos\left(\frac{\theta }{2}\right)\left(-1+4\nu -cos\left(\theta \right)\right)\hfill \end{array}$$

(47)

The numerical experiments are performed with a square domain $\mathsf{\Omega }=[0,5]\times [-2.5,2.5]$ cut by a crack ${\mathsf{\Gamma }}_c=[0,2.5]\times \left\{0\right\}$, see Figure 10. Young’s modulus and Poisson’s ratio are E = 200,000 and $\nu =0.3$, respectively, and the imposed stress intensity factors are $K_I=K_{II}=2802,5$. Following Laborde and coworkers [47], we impose the closed-form displacement field on the boundary of the plate, as depicted in Figure 10. Eight homogeneous meshes with element sizes ranging from $h=0.025$ to $h=1.25$ are employed to discretise the square domain and compare the relative accuracy of the AFEM and the SAFEM under the Mode I and Mode II loadings of the crack. The expression of the relative error in the energy norm was introduced in Section 2.2.2 (cf. Equation (13)) and is not repeated here. The relative error in the L2-norm is given by:

$${\left|\left|\mathbf{e}\right|\right|}_{L2}={\left(\frac{{\int }_{\mathsf{\Omega }}(\mathit{u}-\widehat{\mathit{u}})·(\mathit{u}-\widehat{\mathit{u}})d\mathsf{\Omega }}{{\int }_{\mathsf{\Omega }}\mathit{u}·\mathit{u}d\mathsf{\Omega }}\right)}^{\frac{1}{2}}$$

(48)

where $\mathit{u}$ is the exact displacement field ($\mathit{u}={\mathit{u}}_{\mathit{I}}$ in Mode I and $\mathit{u}={\mathit{u}}_{\mathit{II}}$ in Mode II) and $\widehat{\mathit{u}}$ is the displacement field provided by the chosen numerical method (i.e., the AFEM or the SAFEM).

The energy error of the SAFEM and the AFEM under the Mode I loading of the crack is plotted in Figure 11. One observes that the convergence rates of the AFEM and the SAFEM are identical. This convergence rate value ($\left|\right|e\left|\right|=\mathcal{O}\left(h^{0.5}\right)$) is well understood: it is caused by the crack-induced singularity of the stress field [49]. One also notices that the SAFEM outperforms the AFEM for a given element size and a given number of DOFs.

The error in the energy norm under the Mode II loading of the crack is plotted in Figure 12. The AFEM and the SAFEM share the same convergence rate, $\left|\right|e\left|\right|=\mathcal{O}\left(h^{0.43}\right)$, and are of similar accuracy. As stated previously, a convergence with an error of $\mathcal{O}\left(h^{0.5}\right)$ is expected in this situation. Numerical experiments performed in Ref. [29] revealed that the AFEM was overly soft under Mode II loadings, leading to suboptimal convergence rates in the energy norm. The similarities between the two curves plotted in Figure 12 suggest that the SAFEM shares this drawback.

Although the error levels shown in Figure 11 and Figure 12 may appear significant, it is worth noting that such errors are commonly observed when modelling singularities using the FEM or the XFEM. This observation is supported by the comparison between the FEM and the XFEM in Ref. [47] (refer to Figure 9). Likewise, the quantitative comparison between the AFEM and the FEM presented in [29] (refer to Figures 16 and 17) indicates similar levels of error and a virtually identical performance of the AFEM and the FEM under the Mode I loading of a crack. Therefore, the fact that the SAFEM is either as performant as (under Mode II loading) or outperforms the AFEM (under Mode I loading) makes it a highly attractive option for researchers and practitioners in the field.

Let us now consider the error in the L2-norm, also denoted as the displacement error. A convergence with an error of $\mathcal{O}\left(h\right)$ is expected in this situation [50]. Figure 13 compares the relative performances of the SAFEM and the AFEM under the Mode I loading of the crack. The differences between the two methods are considerable: the SAFEM converges at $\mathcal{O}\left(h^{1.32}\right)$, while the AFEM is suboptimal with an $\mathcal{O}\left(h^{0.52}\right)$ convergence rate. A similar situation is also observed under the Mode II loading of the crack where the SAFEM outperforms the AFEM, see Figure 14. In that case, the SAFEM convergence rate, $\mathcal{O}\left(h^{0.74}\right)$, is suboptimal but remains higher than the AFEM one, $\mathcal{O}\left(h^{0.5}\right)$.

Additional post-treatments shed light on the underperformance of the AFEM in the L2-norm. To allow for a detailed comparison between the AFEM and the SAFEM, the displacement error at the element level under the Mode I loading of the crack is plotted in Figure 15. This plot reveals that the displacement error is not governed by the crack-tip singularity but by the incompatibility of the interelement displacement field, inherent to the use of augmented elements. The SAFEM reduces this incompatibility, which explains its better accuracy.

This aspect is quantitatively illustrated in Figure 16 where one compares the exact and the approximate normal opening of the crack under Mode I loading. The normal opening of the crack, $\Delta u_n\left(r\right)$, is defined as:

$$\Delta u_n\left(r\right)=u_{Iy}(r,\theta =\pi )-u_{Iy}(r,\theta =-\pi )r\in [0,2.5]$$

(49)

where $u_{Iy}$ is the y-component of the closed-form displacement field, see Equation (44).

As shown in Figure 16, (i) the crack opening computed with the AFEM oscillates strongly, (ii) these oscillations do not significantly decrease with mesh refinement and (iii) the crack openings do not fluctuate around the exact values. In contrast, the crack openings computed with the SAFEM consistently converge toward the exact values. The same trends are observed under the Mode II loading of the crack (not shown for the sake of brevity).

The SAFEM and the AFEM exhibit similar performances in terms of the errors in the energy norm, with a slight advantage for the SAFEM. This error measure is dominated by the singularity of the stress field at the crack tip (cf. [49] or Figure 18 of Ref. [16]), where standard finite elements, specifically the T3 or the T3A, are used. The T3 and the T3A have identical convergence rates in the energy norm (Figure 5), which explains the similar performances of the AFEM and the SAFEM in this norm. However, when examining the displacement error, it becomes apparent that the interelement discontinuities of the augmented elements, namely, the AT3 or the AT3A, dominate the error (Figure 15). The SAFEM predicts the crack opening much more accurately than the AFEM (Figure 16). The stabilisation brought by the drilling DOFs is thus responsible for the superior performance of the SAFEM over the AFEM.

Despite the increased problem size caused by the use of drilling DOFs, the improved accuracy of the SAFEM more than compensates for it. Therefore, it can be concluded that the SAFEM is superior to the AFEM in modelling traction-free cracks for the reasons stated in

Table 2. Benefits and drawbacks of the SAFEM over the AFEM when modelling traction-free cracks.

Pros	Cons
slightly more accurate in the energy norm, under Mode I loading of a crack (Figure 11)	additional implementation effort required due to the use of Allman’s elements
faster convergence in the L2-norm (Figure 13 and Figure 14)
crack openings converge toward expected values (Figure 16)

.

4.2. Cohesive Strong Discontinuities

The AFEM was originally developed for cohesive crack propagation [6]. In this section, we shift our focus to cohesive crack modelling and present two sets of tests to compare the performance of the AFEM and SAFEM. The first set of tests examines the satisfaction of a discontinuous patch test, detailed in Section 4.2.1. The second set involves a delamination test case, discussed in Section 4.2.2.

4.2.1. A Discontinuous Patch Test

Patch tests are numerical problems that have exact solutions that lie within the span of a chosen numerical method. Passing patch tests requires that the numerical results match the exact solution. Such tests are often performed during the design of new finite elements as they may provide necessary and sufficient requirements for convergence (cf. Chap. 3 of Ref. [5]). To the best of our knowledge, it has not yet been demonstrated whether passing or failing a patch test ensures or prevents the convergence of the numerical methods that deal with embedded cohesive cracks, such as the AFEM. Therefore, the patch test presented in this section should not be considered as a means to rigorously study the convergence of the AFEM, but rather as a way to highlight its behaviour in simple situations. The discontinuous patch test was first proposed by Armero and coworkers [43] and subsequently adopted by several authors working on finite elements with embedded strong discontinuities [44,51]. The structure considered for the patch test is the $1\times 1$ isotropic block depicted in Figure 17. It is initially discretised with crack-free standard finite elements, the T3 and the T3A, and subjected to displacement boundary conditions inducing a uniaxial stress state in the entire domain. The displacement boundary conditions shown in Figure 17 are given by:

$$\begin{array}{c}-u(x=0)=u(x=1)=u_x\end{array}$$

(50)

$$\begin{array}{c}v(x=0,y=0)=v(x=1,y=0)=0\end{array}$$

(51)

where u and v are the horizontal and vertical displacement fields, respectively. Once the material critical tensile stress ${\sigma }^c$ is reached, a cohesive discontinuity is inserted through the block’s height, dividing it into two equally sized parts. The AFEM or SAFEM are utilised to embed the discontinuity within the finite elements in the numerical experiments. The structure undergoes progressive softening and eventually becomes stress-free as the crack opens. The loading direction is then reversed leading to the closure of the crack and a compressive stress state.

Plane stress conditions are assumed. The bulk material properties are as follows: Young’s modulus is $E=1000$ and Poisson’s ratio is $\nu =0.2$. A bilinear weakly coupled cohesive law (cf., Figure A1) further detailed in Appendix B is employed with the following material properties: the normal and tangential critical tractions are ${\sigma }^c=100$ and ${\tau }^c=10$, the normal and tangential critical openings are ${{\delta }_n}^c={{\delta }_s}^c=0.001$ and the normal and tangential final openings are ${{\delta }_n}^f={{\delta }_s}^f=0.3$.

The closed-form solution for this test, which demonstrates that the block experiences a uniform and uniaxial stress state ${\sigma }_{xx}$, is provided in Ref. [43]. We introduce an equivalent stress, $\overline{{\sigma }_{xx}}$, to assess the distance between the numerical and exact solutions. It is computed thanks to the horizontal reaction force, $f_x$ (cf. Figure 17), measured through the numerical experiment and given by $\overline{{\sigma }_{xx}}=f_x/S_x$, where $S_x$ is the surface of the block normal to the loading direction (equal to 1 in the present case). This equivalent stress makes it easier to compare the numerical and exact solutions as well as to analyse the influence of mesh refinement. If the numerical method matches the exact solution, $\overline{{\sigma }_{xx}}={\sigma }_{xx}$, and the h-convergence of the numerical method is simply represented by the curves $\overline{{\sigma }_{xx}}$ and ${\sigma }_{xx}$ getting closer. Furthermore, because the evaluation of $\overline{{\sigma }_{xx}}$ involves computing reaction forces, it provides access to the total energy dissipated by the structure, a key quantity in cohesive crack modelling.

Another crucial aspect of cohesive crack modelling is the integration scheme of the cohesive forces (Equations (29) and (30)). This is a widely debated topic, with varying suggestions found in the literature (see, e.g., [38,52,53,54,55]). To study the influence of the integration scheme of the cohesive forces, we employed seven integration schemes (ISs), described in Chap. 5.9 of Ref. [5], including one-point ISs; two-, three- and four-point Gauss ISs; and two-, three- and four-point Newton–Cotes ISs.

Three meshes made of randomly oriented triangular elements have been used to study the influence of spatial discretisations on the results (Figure 18). The intersections between the strong discontinuity surface and the triangular elements are arbitrarily close from the mesh nodes.

In Figure 19, we present the results obtained using the coarsest spatial discretisation and one-point as well as Gauss integration schemes. Some notable observations can be made regarding both the AFEM and the SAFEM. Firstly, they fail the patch test with all the tested integration schemes. Secondly, a sharp stress drop occurs once the strong discontinuity is introduced. Thirdly, the compressive stress state is not accurately simulated due to the interelement discontinuities inherent in the use of augmented elements. Instead of exerting pressure on the crack lips, the triangular subdomains rotate, as shown in Figure 19. This flaw was found to gradually disappear with mesh refinement. Still, we recommend using Newton–Cotes integration schemes, which induce nonzero compressive stresses, as demonstrated in Figure 20. Hence, we will focus on Newton–Cotes integration schemes in the subsequent analysis.

The comparison between the AFEM and the SAFEM reveals some differences that hold for all the integration schemes (ISs), as shown in Figure 19 and Figure 20. These differences include (i) the SAFEM exhibiting a softer response in compression than the AFEM, (ii) the SAFEM leading to larger stress drops than the AFEM and (iii) the AFEM producing a non-uniform crack opening under traction-free conditions, which is in stark contrast with the expected constant opening (see Step 3 of Figure 17).

It is worth noting that flaws (i) and (ii) disappear with mesh refinement, as shown in Figure 21 and Figure 22. However, erroneous crack openings computed with the AFEM are not corrected by finer discretisations as illustrated in Figure 21. Therefore, once again, the SAFEM exhibits superior performance over the AFEM.

Furthermore, our tests demonstrate that the three-point Newton–Cotes IS provides optimal results in this test case, as it is both similar to the four-point Newton–Cotes IS and more accurate than the two-point Newton–Cotes IS. We thus chose the three-point Newton–Cotes scheme to visualise the impact of mesh refinement on the results in Figure 22. Interestingly, we found that further mesh refinement does not improve the solutions significantly, as medium and fine meshes yield similar results. This conclusion holds for other integration schemes as well.

To evaluate the discrepancy between the exact and numerical solutions quantitatively, we chose to compare the dissipated energy. The exact dissipated energy is the product of the cohesive surface fracture energy under the Mode I loading $G_I^c$ (as given by Equation (A6)) and the crack surface area $S_x$, which is equal to 1 in this specific case. Thus, $E_{exact}=G_I^c\times S_x$. The dissipated energies associated with the numerical experiments are computed as $E_{numerical}=2\times {\int }_{u=0}^{max\left(u_x\right)}{f}_x\left(u\right)du$, where $f_x$ is the horizontal reaction force and $max\left(u_x\right)$ is the maximum imposed displacement. The signed relative error in the dissipated energy, denoted as $E_{\%}$, is plotted in Figure 23. This graphic illustrates the superior performances of the three- and four-point Newton–Cotes ISs. It also confirms that neither the AFEM nor the SAFEM converge toward the exact solution with mesh refinement.

4.2.2. Mode I Delamination Test

The previous section has highlighted the possible discrepancies between the computed and expected dissipated energies that may occur when augmented elements are made use of. The present section investigates the impact of this flaw on the force/displacement curves obtained in simulations involving multiaxial stress states. A test case with a known analytical solution [56] has been chosen for that purpose: the double-cantilever beam (DCB), see Figure 24. The DCB is used to simulate delamination propagation under pure Mode I loadings. The simulations exactly replicate those of Turon and coworkers [56]. Plane stress conditions are assumed and orthotropic material properties representative of carbon-fiber-reinforced epoxy composites are used: the longitudinal Young’s modulus is $E_{xx}=120$ GPa, Young’s modulus in the transverse direction is $E_{yy}=10.5$ GPa, the in-plane shear modulus is $G_{xy}=5.25$ GPa and the in-plane Poisson’s ratios are ${\nu }_{xy}={\nu }_{yx}=0.3$. As in Ref. [56], the simulations are performed under displacement control with no specific path-following strategy nor inertia-based stabilisation [8].

Delamination modelling throughout the DCB test is usually performed with a priori inserted cohesive elements located on the expected delamination plane. The use of augmented elements removes this meshing burden because strong discontinuity surfaces are embedded within the elements. The specimen’s initial crack (of coordinates $\{x,y\}\in [0,a_0]\times \{H/2\}$) as well as the expected delamination plane ($\{x,y\}\in [a_0,L]\times \{H/2\}$) are not explicitly meshed but modelled with augmented elements. A traction-free discontinuity is used to represent the initial crack, while the delamination plane is modelled with a weakly coupled bilinear cohesive law (cf. Figure A1), further detailed in Appendix B, and the material properties taken from Ref. [56]. The Mode I and Mode II fracture energies are $G_I^c=0.26$ N/mm and $G_{II}^c=1.002$ N/mm, the normal and tangential elastic stiffnesses of the cohesive surface are $K_n=K_s={10}^6$ N/mm${}^3$ and the normal and tangential critical stresses are ${\sigma }^c=30$ MPa and ${\tau }^c=120$ MPa. The cohesive forces are integrated with a three-point Newton–Cotes IS. Three levels of mesh refinement (coarse, medium and fine) are considered in the numerical experiments (Figure 25). The intersections between the strong discontinuity surface and the triangular elements are arbitrarily close from the mesh nodes.

Figure 26 compares the numerical and exact force/opening curves obtained when running the DCB test. The observed oscillations are typical in such simulations and tend to disappear with mesh refinement [52,54]. The SAFEM is less prone to oscillations than the AFEM in the present case. The overall results provided by the AFEM and the SAFEM are similar, and both methods converge toward the exact solution with mesh refinement. The convergence is nevertheless deemed “slow”. Indeed, the use of mesh-conforming cohesive elements with a $0.3$ mm Q4 provides force/opening curves that are closer to the analytical solution than those obtained in this study with the $0.1$ mm augmented elements [56].

The reason behind this slow convergence can be grasped by an inspection of the longitudinal stress field ${\sigma }_{xx}$. As shown in Figure 27, ${\sigma }_{xx}$ is too low within the augmented elements: they should support the highest longitudinal stress, due to the bending of the two arms, but interelement discontinuities limit the amount of transmitted stresses. This is in line with the observations made in Ref. [29] and highlights the overly soft behaviour of the augmented elements when the crack lips are loaded. Figure 27 puts to light the advantages of augmented elements with drilling DOFs: their use results in a consistent opening of the crack. In contrast, the crack opening computed with the AFEM is inaccurate.

The initial stiffness of the AFEM simulations (represented by the slope of the force/opening curves in Figure 26) remains close to the expected value, even with the coarsest mesh. In contrast, the SAFEM produces an overly soft result with the coarsest mesh. The reason behind this behaviour is straightforward: the AT3 is too soft due to insufficient transmitted stresses (Figure 27), while the T3 element is too stiff (cf. Table 1). These two negative aspects compensate well in the current situation. The T3A is not overly stiff (cf. Table 1), while the AT3A is too soft for the same reason as the AT3. As a result, SAFEM structures are softer than those modelled with the AFEM. The differences disappear with mesh refinement. The benefits and the drawbacks of the SAFEM over the AFEM when modelling cohesive cracks are gathered in

Table 3. Benefits and drawbacks of the SAFEM over the AFEM when modelling cohesive cracks.

Pros	Cons
smoother force/displacement results in DCB test (Figure 26)	additional implementation effort required due to the use of Allman’s elements
consistent crack openings (Figure 21 and Figure 27)

.

5. Conclusions

The AFEM is a validated method that has proven successful in modelling dynamic crack propagation or conducting three-dimensional studies of damage in composite materials. As well known, some combinations of augmented element geometry and crack orientation give rise to singular stiffness matrices when modelling traction-free cracks. This work successfully introduced a method to suppress these singularities, named the stabilised augmented finite element method (SAFEM). The SAFEM relies on the use of bulk elements possessing rotational DOFs. The so-called Allman elements with drilling DOFs have been used for that purpose.

Numerical experiments were conducted to compare the performance of the SAFEM with the traditional AFEM in modelling traction-free cracks and cohesive discontinuities in Mode I and Mode II. The results showed that the AFEM suffers from degraded convergence in the L2-norm and can produce erroneous (cf. Figure 16) and irregular (cf. Figure 27) crack openings. These limitations restrict the applicability of the AFEM in situations where accurate crack opening estimation is essential, such as chemo-mechanical couplings in ceramic matrix composites [57,58]. In contrast, the proposed SAFEM successfully addresses these issues, limiting interelement discontinuities and achieving an accurate estimation of crack openings when modelling traction-free cracks (Figure 16). The improved accuracy of the SAFEM compensates for the computational cost incurred by the introduction of drilling DOFs (cf. Figure 13 and Figure 14).

The SAFEM’s superiority over the AFEM is less pronounced when it comes to modelling cohesive discontinuities. The presence of interelement discontinuities can lead to a premature failure of the cohesive surface, which negatively impacts the accuracy of both the AFEM and the SAFEM. Although both methods have their merits, there is still room for improvement in achieving completely satisfactory results in this area.

This study has demonstrated, for the first time, that even the most straightforward drilling DOF formulation can noticeably enhance the AFEM performance. In future research, efforts will be made to develop tailored drilling DOFs to reduce interelement discontinuities when dealing with cohesive discontinuities.