# Structural cohesive element for the modelling of delamination in composite
laminates without the cohesive zone limit

Xiaopeng Ai  Boyang Chen<EMAIL_ADDRESS>Christos Kassapoglou

###### Abstract

Delamination is a critical mode of failure that occurs between plies in a
composite laminate. The cohesive element, developed based on the cohesive zone
model, is widely used for modeling delamination. However, standard cohesive
elements suffer from a well-known limit on the mesh density—the element size
must be much smaller than the cohesive zone size. This work develops a new set
of elements for modelling composite plies and their interfaces in 3D. A
triangular Kirchhoff-Love shell element is developed for orthotropic materials
to model the plies. A _structural_ cohesive element, conforming to the shell
elements of the plies, is developed to model the interface delamination. The
proposed method is verified and validated on the classical benchmark problems
of Mode I, Mode II, and mixed-mode delamination of unidirectional laminates,
as well as on the single-leg bending problem of a multi-directional laminate.
All the results show that the element size in the proposed models can be ten
times larger than that in the standard cohesive element models, with more than
90% reduction in CPU time, while retaining prediction accuracy. This would
then allow more effective and efficient modeling of delamination in composites
without worrying about the cohesive zone limit on the mesh density.

###### keywords:

Composites , delamination , cohesive zone model , cohesive element

††journal: Computer Methods in Applied Mechanics and Engineering

[inst1]organization=Department of Aerospace Structures and Materials, Faculty
of Aerospace Engineering, Delft University of
Technology,addressline=Kluyverweg 1, city=Delft, postcode=2629 HS, country=The
Netherlands

##  1 Introduction

The accurate prediction of delamination is of critical importance for the
reliable design of fiber-reinforced composite structures. The cohesive element
(CE) is a widely used finite element technology to model delamination. CE is
developed based on the Cohesive Zone Model proposed by Dugdale and Barenblatt
[1, 2]. A fracture process zone, generally called the cohesive zone, exists
along the interface, ahead of the stress-free crack tip. A traction-separation
relationship, namely the cohesive law, describes how the interfacial stresses
and damage evolve with respect to the interfacial openings. Standard CEs are
usually developed for use between two solid elements to model their debonding
[3, 4, 5, 6, 7, 8, 9, 10].

While being a popular and versatile tool to model delamination, standard CEs
suffer from a well-known limit on the mesh density – the element size must be
much smaller than the Cohesive Zone Length (CZL) to accurately predict
delamination. According to the previous analysis [3, 4], high-stress gradients
could be produced within the cohesive zone during delamination in composites.
A very fine CE mesh must be used there to sufficiently capture the stress
gradients such that the internal virtual work of the CEs can be accurately
integrated. So far, there is no fixed rule on how fine the mesh should be in
the literature. In some research [5, 6, 7], the authors have demonstrated that
at least two or three CEs should be used inside the cohesive zone. With
coarser meshes, simulations would significantly over-predict the peak load
[8]. In a typical mode I delamination test of a unidirectional composites
coupon, i.e., the Double Cantilever Beam (DCB) test, the CZL is less than 1
mm.

The above-mentioned problem of cohesive zone limit on the mesh density of CEs
has drawn the attention of many researchers in the past. Turon et al. [8]
adopted an engineering method to solve the problem by reducing the material
strength to numerically extend the CZL. While shown to be promising in the
Mode I DCB case, in cases such as the pure mode II End-Notched-Flexure (ENF)
test, decreasing the cohesive strength can cause excessive under-prediction of
the overall strength [11, 12]. Yang [13] and Do [14]’s research work
demonstrated that larger CE sizes can be achieved by placing enough
integration points in the cohesive zone. Their method could predict the peak
load correctly for CE size up to 1.43 times the estimated CZL. However, its
first-order shape function limits the effectiveness of this method in further
expanding CE beyond the scale of the cohesive zone [15]. Guiamatsia et al.[16,
17, 18] used the beam on elastic foundation solution as an enrichment function
and tested it in the mixed mode delamination. However, the enrichment method
could lead to inaccuracies in interpolation for elements larger than 3 mm
[17]. Another enrichment approach with the piecewise linear shape functions
was proposed by Samimi [19, 20, 21]. However, the over-prediction problem of
peak load under large CE has not been solved. Van der Meer et al. [22] used
the level set method with a energy-based criterion to propagate delamination
without the cohesive zone limit. However, this method is limited to the case
of a single delamination. Lu et al. [23] proposed a adaptive version of the
Floating Node Method (FNM) [24] to adaptively refine the CEs in the cohesive
zone. However, this method is currently implemented for 2D problems. Alvarez
et al. [25] used quadratic CEs between quadratic solid elements and used
higher number of integration points, which increased the grid size to be
comparable to CZL but not beyond [15]. Mukhopadhyay and Bhatia [26] developed
a $hp$ refinement strategy to simulate the delamination between two solid
elements. However, it is only implemented in 2D in their work. Daniel [27]
developed an ERR-Cohesive method to simulate the delamination with large
elements by estimating the energy release rate by means of the virtual crack
closure technique (VCCT). However, this method is only presented in 2D, and
the VCCT approach relies on the existence of an initial crack and the
assumption of self-similar crack propagation.

Inspired by the earlier works in the literature, Russo and Chen[15] developed
a so-called structural CE, which conforms to Kirchhoff-Love structural
elements for the neighbouring plies. Their work was done in 2D. The Euler-
Bernoulli beam elements were used to model the plies. The structural CE,
sharing its nodes with the beam elements, was developed to model the
delamination. An adaptive integration scheme was used to place more
integration points in CEs containing the cohesive zone. Their results showed
that the structural CE could overcome the cohesive zone limit on mesh density,
allowing the element size to be ten times larger than that of the standard
linear CEs. Motivated by Russo and Chen’s work, Tosti Balducci and Chen [28]
extended the structural CE to 3D DCB problem by developing a structural CE
compatible with the TUBA3 plate elements [29]. Their results showed that the
TUBA3-based structural formulation of CE could overcome the cohesive zone
limit. However, the curvature degrees of freedom (DoFs) make it more
complicated to set boundary conditions on TUBA3 elements than on more common
elements such as solid, beam or Reissner-Mindlin shell elements, thereby
impeding the adoption of such elements by engineers in practice.

From the above reviews, we can see that a composites delamination model which:

  1. 1.

does not suffer from the cohesive zone limit on mesh density;

  2. 2.

works in 3D space;

  3. 3.

does not require artificial reduction of strengths;

  4. 4.

does not require a pre-crack and the self-similar propagation assumption;

  5. 5.

can model delamination along multiple interfaces; and

  6. 6.

does not require working with curvature DoFs,

is still missing. Based on the earlier works [15, 28], this article aims to
develop a new set of elements for the plies and interfaces of composites. The
idea is to develop a triangular Kirchhoff-Love shell element without curvature
DoFs to model the composite plies, then formulate the conforming structural CE
to be used between the plies to model their delamination. Allman [30] proposed
a simple triangle Kirchhoff-Love plate element named cubic displacement
element. There are only three DoFs at each node, one for out-of-plane
displacement and two for rotations. These DoFs are commonly used by engineers
and are much easier to handle than the curvature ones when setting boundary
conditions. Therefore, this article will extend this triangle plate element
for the modeling of composite plies and develop the corresponding structural
CE for the interfaces. If such a structural CE could overcome the cohesive
zone limit on its mesh density, then it would achieve all the requirements
listed above.

The rest of the paper will be structured in the following way. Section 2
presents the proposed element formulations in detail. Section 3 demonstrates
the performance of these elements on a series of benchmarks on delamination in
both unidirectional and multi-directional laminates. In the end, Section 4
draws the conclusions of this work and discusses some potential future work.

##  2 Method

###  2.1 Overall illustration

(a) Conventional CE between solids: interface thickness = distance between
node 4 and node 8

(b) Structural CE between shells: interface thickness $\neq$ distance between
node 4 and node 8

Figure 1: Geometrical comparison between conventional CE and structural CE

The overall idea of the proposed modelling approach is to represent the
composite plies by Kirchoff-Love shell elements and the interfaces by
structural CEs. Figure 1 shows the geometrical comparison between the
conventional modelling approach with solid elements and the proposed
structural approach with shells. Looking at the structural CE, the difference
from the conventional CE is that the nodes of the CE are placed at the mid-
plane of the two shell elements, which are not the actual surfaces of the
cohesive interface. Hence, the opening of the actual interface does not equate
to the distance between the upper nodes and lower ones of the structural CE.
This opening shall be calculated using the shell kinematics to be detailed in
later subsections. The rest of this section will firstly present the cubic
plate element formulation by Allman [30] and its adaptation for flat
composites shell, then move on to derive the structural CE formulation.

###  2.2 Cubic plate element

####  2.2.1 Geometric definitions of the plate element

Figure 2: Coordinate system and DoFs for the triangular cubic plate element
[30]

The plate part of the shell element in this work is based on the triangular
cubic plate element developed by Allman [30], as shown in Figure 2. A local
coordinate system is used in this element, with origin at the centroid of the
triangle element. The local axis are represented by lowercase letters, $x$ and
$y$, to distinguish them from the global ones, $X$ and $Y$. The area of this
triangular element is $A$. $s$ is the anti-clockwise coordinate along the
element boundary. $n$ is the exterior normal. The angle between the normal $n$
and the local axis $x$ is $\gamma$. The out-of-plane displacement defined over
the domain is ${w}(x,y)$. An independent out-of-plane boundary displacement,
$\overline{w}(s)$, and its compatible normal derivative,
$\nicefrac{{\partial\overline{w}}}{{\partial n}}(s)$, are assumed along the
boundary. The DoFs at each node include the displacement $\overline{w}$ and
two rotations $\nicefrac{{\partial\overline{w}}}{{\partial x}}$,
$\nicefrac{{\partial\overline{w}}}{{\partial y}}$.

The cubic plate element in Allman’s work only considered the case of isotropic
material [30]. This work extends the cubic plate element using the classical
laminate theory, such that symmetric composite laminates can also be modelled
by this plate element. The original cubic element formulation will be
presented in detail, with adaptations for composites specified along the way.

####  2.2.2 Modified potential energy for the plate element

The minimum potential energy principle is used to derive the finite element
formulation of the triangular plate element. The potential energy used is
referred to as the modified potential energy in Allman’s work [30]:

$\displaystyle\pi=\iint_{A}U_{0}\mathrm{d}x\mathrm{d}y+\sum_{N=1}^{3}R_{N}(\overline{w}_{N}-w_{N})+\int_{\partial
A}V_{n}(\overline{w}-w)\,\mathrm{d}s-\int_{\partial
A}M_{n}(\frac{\partial\overline{w}}{\partial n}-\frac{\partial w}{\partial
n})\mathrm{d}s$
$\displaystyle-\iint_{A}p^{*}w\mathrm{d}x\mathrm{d}y-\sum_{N=1}^{3}R^{*}_{N}\overline{w}_{N}-\int_{\partial
A}V^{*}_{n}\overline{w}\mathrm{d}s+\int_{\partial
A}M^{*}_{n}\frac{\partial\overline{w}}{\partial n}\,\mathrm{d}s$ (1)

where $U_{0}$ is the strain energy density. $M_{n}$ is the normal bending
moment resultant and $V_{n}$ is the Kirchhoff shear force distribution on the
element boundary, respectively. $R_{N\,(N=1,2,3)}$ are the concentrated forces
at the element vertices. In addition, $R^{*}_{N}$, $V^{*}_{n}$, and
$M^{*}_{n}$ are the values of prescribed concentrated force, Kirchhoff shear
force, and normal bending moment resultant, respectively. $p^{*}(x,y)$ denotes
the prescribed distributed pressure load on the element. This potential energy
essentially uses $R_{N}$, $V_{n}$, and $M_{n}$ as Lagrange multipliers to
enforce the compatibility between the two fields $w$ and $\overline{w}$ and
their normal derivatives along the boundary. The variational principle based
on this modified potential energy satisfies the equilibrium equations, the
boundary conditions, and the compatibility requirement between $w$ and
$\overline{w}$ [30].

The expression of $U_{0}$ for composites will be different from that for
isotropic materials in Allman’s work [30]. The classical laminate theory can
be used to describe the constitutive relationship of a composite laminate
under Kirchhoff-Love kinematic assumptions. For the case of a symmetric
laminate (the smallest of which would be a single composite ply) considered in
this work, there would be no membrane-bending coupling. Therefore, the
expressions of the moments are:

$\displaystyle M_{x}=-D_{11}\frac{\partial^{2}w}{\partial
x^{2}}-D_{12}\frac{\partial^{2}w}{\partial
y^{2}}-2D_{16}\frac{\partial^{2}w}{\partial x\partial y}$ $\displaystyle
M_{y}=-D_{12}\frac{\partial^{2}w}{\partial
x^{2}}-D_{22}\frac{\partial^{2}w}{\partial
y^{2}}-2D_{26}\frac{\partial^{2}w}{\partial x\partial y}$ $\displaystyle
M_{xy}=-D_{16}\frac{\partial^{2}w}{\partial
x^{2}}-D_{26}\frac{\partial^{2}w}{\partial
y^{2}}-2D_{66}\frac{\partial^{2}w}{\partial x\partial y}$ (2)

Using the above moment-curvature relations (equation 2.2.2), the strain energy
density $U_{0}$ can be written as:

$\displaystyle U_{0}$
$\displaystyle=\frac{1}{2}\left[D_{11}\left(\frac{\partial^{2}w}{\partial
x^{2}}\right)^{2}+D_{22}\left(\frac{\partial^{2}w}{\partial
y^{2}}\right)^{2}+2D_{12}\frac{\partial^{2}w}{\partial
x^{2}}\frac{\partial^{2}w}{\partial y^{2}}\right.$
$\displaystyle\left.+4D_{16}\frac{\partial^{2}w}{\partial
x^{2}}\frac{\partial^{2}w}{\partial x\partial
y}+4D_{26}\frac{\partial^{2}w}{\partial y^{2}}\frac{\partial^{2}w}{\partial
x\partial y}+4D_{66}\left(\frac{\partial^{2}w}{\partial x\partial
y}\right)^{2}\right]$ (3)

Figure 3: Sign conventions for shear force and bending moment resultants

The sign conventions for the shear force and bending moment resultants are
shown in Figure 3. For the triangular plate element, the normal bending moment
resultant $M_{n}$, Kirchhoff shear force distribution $V_{n}$ and concentrated
forces $R_{N\,(N=1,2,3)}$ can be expressed as [30]:

$\displaystyle M_{n}$
$\displaystyle=M_{x}\cos^{2}{\gamma}+M_{y}\sin^{2}{\gamma}+M_{xy}\sin{2\gamma}$
(4) $\displaystyle V_{n}$ $\displaystyle=Q_{n}+\frac{\partial M_{ns}}{\partial
s}$ (5) $\displaystyle R_{N}$ $\displaystyle=[M_{ns}]^{s^{+}}_{s^{-}}$ (6)

where

$\displaystyle M_{ns}=\frac{1}{2}(M_{y}-M_{x})\sin 2\gamma+M_{xy}\cos 2\gamma$
(7) $\displaystyle Q_{n}=\frac{\partial M_{n}}{\partial n}+\frac{\partial
M_{ns}}{\partial s}$ (8)

In this element, linear variation of the normal bending moment resultant and
constant Kirchhoff shear force distribution are assumed along the edges:

$\displaystyle M_{n}$ $\displaystyle=M^{12}_{n}(1-\xi)+M^{21}_{n}\xi$ (9)
$\displaystyle V_{n}$ $\displaystyle=V_{n}^{12}$ (10)

where $M_{n}^{12}$ and $M_{n}^{21}$ are the normal bending moment resultant
values at node $1$ and node $2$ of side 1-2 respectively. The coefficient
$\xi$ in equation 9 is defined as $\xi=s/l_{12}$. $V_{n}^{12}$ is the value of
the Kirchhoff shear force distribution along the side 1-2. The cyclic
permutation of superscripts 1, 2, and 3 in equations 9 and 10 produces the
shear force distribution and moment resultant on the other two sides (side 2-3
and side 3-1) of the triangular element.

The terms in brackets in equation 6 represent the difference in bending moment
resultant values at the element vertices. Therefore, the expressions of
$R_{N\,(N=1,2,3)}$ are:

$\displaystyle R_{1}=M_{ns}^{12}-M_{ns}^{13}$ $\displaystyle
R_{2}=M_{ns}^{23}-M_{ns}^{21}$ $\displaystyle R_{3}=M_{ns}^{31}-M_{ns}^{32}$
(11)

Finally, the directional derivatives in equations 8 are given by:

$\displaystyle\frac{\partial}{\partial s}=-\sin\gamma\frac{\partial}{\partial
x}+\cos\gamma\frac{\partial}{\partial y}$
$\displaystyle\frac{\partial}{\partial n}=\cos\gamma\frac{\partial}{\partial
x}+\sin\gamma\frac{\partial}{\partial y}$ (12)

####  2.2.3 Finite element approximation and the discretized equilibrium
equation

For finite element implementation, the potential energy $\pi$ should be
discretized. In this element, the internal out-of-plane displacement field
$w(x,y)$ is approximated by a cubic polynomial:

$w(x,y)=A_{1}+A_{2}\,x+A_{3}\,y+\alpha_{1}\,x^{2}+\alpha_{2}\,xy+\alpha_{3}\,y^{2}+\alpha_{4}\,x^{3}+\alpha_{5}\,x^{2}y+\alpha_{6}\,xy^{2}+\alpha_{7}\,y^{3}$
(13)

The three coefficients $A_{1}$, $A_{2}$, and $A_{3}$ represent the rigid body
motion and do not affect the value of the strain energy density.

With the above choice of $w$, applying the Green’s theorem [31] to transform
the strain energy integral to boundary integral gives:

$\displaystyle
2\iint_{A}U_{0}\,\mathrm{d}x\mathrm{d}y=\sum_{N=1}^{3}R_{N}w_{N}+\int_{\partial
A}V_{n}w\,\mathrm{d}s-\int_{\partial A}M_{n}\frac{\partial w}{\partial
n}\,\mathrm{d}s$ (14)

In this work, we assume that there is no distributed pressure load
$p^{*}(x,y)$. Then, the modified potential energy (equation 2.2.2) can be
rewritten more simply as:

$\displaystyle\pi=-\iint_{A}U_{0}\,\mathrm{d}x\mathrm{d}y+\sum_{N=1}^{3}R_{N}\overline{w}_{N}+\int_{\partial
A}V_{n}\overline{w}\,\mathrm{d}s-\int_{\partial
A}M_{n}\frac{\partial\overline{w}}{\partial n}\,\mathrm{d}s$
$\displaystyle-\sum_{N=1}^{3}R^{*}_{N}\overline{w}_{N}-\int_{\partial
A}V^{*}_{n}\overline{w}\,\mathrm{d}s+\int_{\partial
A}M^{*}_{n}\frac{\partial\overline{w}}{\partial n}\,\mathrm{d}s$ (15)

Substituting equation 13 into 2.2.2, the strain energy can be rewritten as:

$\iint_{A}U_{0}\,\mathrm{d}x\mathrm{d}y=\frac{1}{2}\boldsymbol{\mathrm{\alpha}}^{\rm{T}}\,\boldsymbol{\mathrm{H}}\,\boldsymbol{\mathrm{\alpha}}$
(16)

where the vector $\boldsymbol{\mathrm{\alpha}}$ is:

$\boldsymbol{\mathrm{\alpha}}=\\{\alpha_{1},\,\alpha_{2},\,\alpha_{3},\,\alpha_{4},\,\alpha_{5},\,\alpha_{6},\,\alpha_{7}\\}^{\rm{T}}$
(17)

The matrix $\boldsymbol{\mathrm{H}}$ is more complicated to derive for
symmetric composite laminates than for isotropic materials considered in
reference [30]. Its detailed derivation is presented in A.

After the strain energy, the $2^{\mathrm{nd}}$, $3^{\mathrm{rd}}$ and
$4^{\mathrm{th}}$ terms in equation 2.2.3 denote the work of the so-called
generalized forces:

$\displaystyle\sum_{N=1}^{3}R_{N}\overline{w}_{N}+\int_{\partial
A}V_{n}\overline{w}\,\mathrm{d}s-\int_{\partial
A}M_{n}\frac{\partial\overline{w}}{\partial
n}\,\mathrm{d}s=\boldsymbol{\mathrm{Q}}^{\rm{T}}\boldsymbol{\mathrm{q}}$ (18)

where $\boldsymbol{\mathrm{Q}}$ is the vector of the twelve generalized
forces:

$\boldsymbol{\mathrm{Q}}=\\{R_{1},R_{2},R_{3},V_{n}^{12},V_{n}^{23},V_{n}^{31},M_{n}^{12},M_{n}^{21},M_{n}^{23},M_{n}^{32},M_{n}^{31},M_{n}^{13}\\}^{\rm{T}}$
(19)

Substituting equations 9 and 10 into equation 18, the generalized
displacements corresponding to the generalized forces in equation 19 compose
the vector $\boldsymbol{\mathrm{q}}$:

$\displaystyle\boldsymbol{\mathrm{q}}=$
$\displaystyle\left\\{\overline{w}_{1},\,\overline{w}_{2},\,\overline{w}_{3},\;l_{12}\int_{0}^{1}\overline{w}\,\mathrm{d}\xi,\;l_{23}\int_{0}^{1}\overline{w}\,\mathrm{d}\xi,\;l_{31}\int_{0}^{1}\overline{w}\,\mathrm{d}\xi,\right.$
$\displaystyle\left.-l_{12}\int_{0}^{1}\frac{\partial\overline{w}}{\partial
n}(1-\xi)\,\mathrm{d}\xi,\;-l_{12}\int_{0}^{1}\frac{\partial\overline{w}}{\partial
n}\xi\,\mathrm{d}\xi,\;-l_{23}\int_{0}^{1}\frac{\partial\overline{w}}{\partial
n}(1-\xi)\,\mathrm{d}\xi,\right.$
$\displaystyle\left.-l_{23}\int_{0}^{1}\frac{\partial\overline{w}}{\partial
n}\xi\,\mathrm{d}\xi,\;-l_{31}\int_{0}^{1}\frac{\partial\overline{w}}{\partial
n}(1-\xi)\,\mathrm{d}\xi,\;-l_{31}\int_{0}^{1}\frac{\partial\overline{w}}{\partial
n}\xi\,\mathrm{d}\xi\right\\}^{\rm{T}}$ (20)

A matrix $\boldsymbol{\mathrm{B}}$ can be constructed to connect the vector
$\boldsymbol{\mathrm{Q}}$ and the vector $\boldsymbol{\mathrm{\alpha}}$:

$\boldsymbol{\mathrm{Q}}=\boldsymbol{\mathrm{B}}^{\rm{T}}\boldsymbol{\mathrm{\alpha}}$
(21)

With the moment-curvature relationship for composites in equation 2.2.2 and
the expression of $w$ in equaton 13, the matrix $\boldsymbol{\mathrm{B}}$ has
been derived for composites in this work. The details are shown in B.

A matrix $\boldsymbol{\mathrm{T}}$ can be defined to represent the
relationship between the nodal DoFs $\overline{\boldsymbol{\mathrm{W}}}$ and
the vector $\boldsymbol{\mathrm{q}}$:

$\displaystyle\boldsymbol{\mathrm{q}}=\boldsymbol{\mathrm{T}}\overline{\boldsymbol{\mathrm{W}}}$
(22)

where $\overline{\boldsymbol{\mathrm{W}}}$ is the vector of the nodal DoFs
defined by

$\displaystyle\overline{\boldsymbol{\mathrm{W}}}=\left\\{\overline{w}_{1},\frac{\partial\overline{w}_{1}}{\partial
x},\frac{\partial\overline{w}_{1}}{\partial
y},\overline{w}_{2},\frac{\partial\overline{w}_{2}}{\partial
x},\frac{\partial\overline{w}_{3}}{\partial
y},\overline{w}_{3},\frac{\partial\overline{w}_{3}}{\partial
x},\frac{\partial\overline{w}_{3}}{\partial y}\right\\}^{\rm{T}}$ (23)

With the assumptions of cubic line function for $\overline{w}$ and linear
variation for $\nicefrac{{\partial\overline{w}}}{{\partial n}}$ along each
edge of the element, the $\boldsymbol{\mathrm{T}}$ matrix can derived as [30]:

$\\+bmatrix[colspec={@{~{}}ccccccccc@{~{}}},colsep=2.2pt]1&00000000\\\
000100000\\\ 000000100\\\
\dfrac{l_{12}}{2}\scalebox{0.75}[1.0]{$-$}\dfrac{l^{2}_{12}}{12}\sin{\gamma_{12}}\dfrac{l^{2}_{12}}{12}\cos{\gamma_{12}}\dfrac{l_{12}}{2}\dfrac{l^{2}_{12}}{12}\sin{\gamma_{12}}\scalebox{0.75}[1.0]{$-$}\dfrac{l^{2}_{12}}{12}\cos{\gamma_{12}}000\\\
000\dfrac{l_{23}}{2}\scalebox{0.75}[1.0]{$-$}\dfrac{l^{2}_{23}}{12}\sin{\gamma_{23}}\dfrac{l^{2}_{23}}{12}\cos{\gamma_{23}}\dfrac{l_{23}}{2}\dfrac{l^{2}_{23}}{12}\sin{\gamma_{23}}\scalebox{0.75}[1.0]{$-$}\dfrac{l^{2}_{23}}{12}\cos{\gamma_{23}}\\\
\dfrac{l_{31}}{2}\dfrac{l^{2}_{31}}{12}\sin{\gamma_{31}}\scalebox{0.75}[1.0]{$-$}\dfrac{l^{2}_{31}}{12}\cos{\gamma_{31}}000\dfrac{l_{31}}{2}\scalebox{0.75}[1.0]{$-$}\dfrac{l^{2}_{31}}{12}\sin{\gamma_{31}}\dfrac{l^{2}_{31}}{12}\cos{\gamma_{31}}\\\
0\scalebox{0.75}[1.0]{$-$}\dfrac{l_{12}}{3}\cos{\gamma_{12}}\scalebox{0.75}[1.0]{$-$}\dfrac{l_{12}}{3}\sin{\gamma_{12}}0\scalebox{0.75}[1.0]{$-$}\dfrac{l_{12}}{6}\cos{\gamma_{12}}\scalebox{0.75}[1.0]{$-$}\dfrac{l_{12}}{6}\sin{\gamma_{12}}000\\\
0\scalebox{0.75}[1.0]{$-$}\dfrac{l_{12}}{6}\cos{\gamma_{12}}\scalebox{0.75}[1.0]{$-$}\dfrac{l_{12}}{6}\sin{\gamma_{12}}0\scalebox{0.75}[1.0]{$-$}\dfrac{l_{12}}{3}\cos{\gamma_{12}}\scalebox{0.75}[1.0]{$-$}\dfrac{l_{12}}{3}\sin{\gamma_{12}}000\\\
0000\scalebox{0.75}[1.0]{$-$}\dfrac{l_{23}}{3}\cos{\gamma_{23}}\scalebox{0.75}[1.0]{$-$}\dfrac{l_{23}}{3}\sin{\gamma_{23}}0\scalebox{0.75}[1.0]{$-$}\dfrac{l_{23}}{6}\cos{\gamma_{23}}\scalebox{0.75}[1.0]{$-$}\dfrac{l_{23}}{6}\sin{\gamma_{23}}\\\
0000\scalebox{0.75}[1.0]{$-$}\dfrac{l_{23}}{6}\cos{\gamma_{23}}\scalebox{0.75}[1.0]{$-$}\dfrac{l_{23}}{6}\sin{\gamma_{23}}0\scalebox{0.75}[1.0]{$-$}\dfrac{l_{23}}{3}\cos{\gamma_{23}}\scalebox{0.75}[1.0]{$-$}\dfrac{l_{23}}{3}\sin{\gamma_{23}}\\\
0\scalebox{0.75}[1.0]{$-$}\dfrac{l_{31}}{6}\cos{\gamma_{31}}\scalebox{0.75}[1.0]{$-$}\dfrac{l_{31}}{6}\sin{\gamma_{31}}0000\scalebox{0.75}[1.0]{$-$}\dfrac{l_{31}}{3}\cos{\gamma_{31}}\scalebox{0.75}[1.0]{$-$}\dfrac{l_{31}}{3}\sin{\gamma_{31}}\\\
0\scalebox{0.75}[1.0]{$-$}\dfrac{l_{31}}{3}\cos{\gamma_{31}}\scalebox{0.75}[1.0]{$-$}\dfrac{l_{31}}{3}\sin{\gamma_{31}}0000\scalebox{0.75}[1.0]{$-$}\dfrac{l_{31}}{6}\cos{\gamma_{31}}\scalebox{0.75}[1.0]{$-$}\dfrac{l_{31}}{6}\sin{\gamma_{31}}$
(24)

####  2.2.4 Stiffness matrix and force vectors of cubic plate element

Substituting equations 16, 21 and 22 into equation 2.2.3, the total modified
potential energy represented by the finite element method under prescribed
boundary loads is:

$\displaystyle\pi=-\frac{1}{2}\boldsymbol{\mathrm{\alpha}}^{\rm{T}}\boldsymbol{\mathrm{H}}\boldsymbol{\mathrm{\alpha}}+\boldsymbol{\mathrm{\alpha}}^{\rm{T}}(\boldsymbol{\mathrm{B}}\boldsymbol{\mathrm{T}})\overline{\boldsymbol{\mathrm{W}}}-\boldsymbol{\mathrm{Q}}^{*\rm{T}}\boldsymbol{\mathrm{T}}\overline{\boldsymbol{\mathrm{W}}}$
(25)

where the vector $\boldsymbol{\mathrm{Q}}^{*}$, whose components are obtained
by replacing the generalized forces in vector $\boldsymbol{\mathrm{Q}}$ with
the corresponding prescribed quantities, denotes all the external forces:

$\displaystyle\boldsymbol{\mathrm{Q}}^{*}=\\{R_{1}^{*},R_{2}^{*},R_{3}^{*},V_{n}^{*12},V_{n}^{*23},V_{n}^{*31},M_{n}^{*12},M_{n}^{*21},M_{n}^{*23},M_{n}^{*32},M_{n}^{*31},M_{n}^{*13}\\}^{\rm{T}}$
(26)

Thus, the internal work $U$ and the external work $W$ of the potential energy
can be expressed as:

$\displaystyle U$
$\displaystyle=-\frac{1}{2}\boldsymbol{\mathrm{\alpha}}^{\rm{T}}\boldsymbol{\mathrm{H}}\boldsymbol{\mathrm{\alpha}}+\boldsymbol{\mathrm{\alpha}}^{\rm{T}}(\boldsymbol{\mathrm{B}}\boldsymbol{\mathrm{T}})\overline{\boldsymbol{\mathrm{W}}}$
(27) $\displaystyle W$
$\displaystyle={\boldsymbol{\mathrm{Q}}^{*}}^{\rm{T}}\boldsymbol{\mathrm{T}}\overline{\boldsymbol{\mathrm{W}}}$
(28)

Based on equation 25, the minimum total potential energy principle gives:

$\displaystyle\delta\pi=\delta\boldsymbol{\mathrm{\alpha}}^{\rm{T}}[(\boldsymbol{\mathrm{B}}\boldsymbol{\mathrm{T}})\overline{\boldsymbol{\mathrm{W}}}-\boldsymbol{\mathrm{H}}\boldsymbol{\mathrm{\alpha}}]+\delta\overline{\boldsymbol{\mathrm{W}}}^{\rm{T}}[(\boldsymbol{\mathrm{B}}\boldsymbol{\mathrm{T}})^{\rm{T}}\boldsymbol{\mathrm{\alpha}}-\boldsymbol{\mathrm{T}}^{\rm{T}}\boldsymbol{\mathrm{Q}}^{*}]=0$
(29)

Setting the coefficient of the arbitrary variation
$\delta\boldsymbol{\mathrm{\alpha}}^{\rm{T}}$ to zero gives:

$\displaystyle\boldsymbol{\mathrm{\alpha}}=\boldsymbol{\mathrm{H}}^{-1}(\boldsymbol{\mathrm{B}}\boldsymbol{\mathrm{T}})\overline{\boldsymbol{\mathrm{W}}}$
(30)

Performing the variation of $U$ as expressed in 27 with
$\boldsymbol{\mathrm{\alpha}}$ substituted by equation 30, we arrive at
$\delta U$ for the element as:

$\displaystyle\delta U$
$\displaystyle=\delta\overline{\boldsymbol{\mathrm{W}}}^{\rm{T}}(\boldsymbol{\mathrm{B}}\boldsymbol{\mathrm{T}})^{\rm{T}}\boldsymbol{\mathrm{H}}^{-1}(\boldsymbol{\mathrm{B}}\boldsymbol{\mathrm{T}})\overline{\boldsymbol{\mathrm{W}}}=\delta\overline{\boldsymbol{\mathrm{W}}}^{\rm{T}}\boldsymbol{\mathrm{f}}_{\rm{int}}$
(31)

where $\boldsymbol{\mathrm{f}}_{\rm{int}}$ is the internal force vector of the
plate element:

$\displaystyle\boldsymbol{\mathrm{f}}_{\rm{int}}=(\boldsymbol{\mathrm{B}}\boldsymbol{\mathrm{T}})^{\rm{T}}\boldsymbol{\mathrm{H}}^{-1}(\boldsymbol{\mathrm{B}}\boldsymbol{\mathrm{T}})\,\overline{\boldsymbol{\mathrm{W}}}=\boldsymbol{\mathrm{K}}_{\mathrm{plate}}\,\overline{\boldsymbol{\mathrm{W}}}$
(32)

and $\boldsymbol{\mathrm{K}}_{\mathrm{plate}}$ is the stiffness matrix of the
plate element:

$\displaystyle\boldsymbol{\mathrm{K}}_{\mathrm{plate}}=(\boldsymbol{\mathrm{B}}\boldsymbol{\mathrm{T}})^{\rm{T}}\boldsymbol{\mathrm{H}}^{-1}(\boldsymbol{\mathrm{B}}\boldsymbol{\mathrm{T}})$
(33)

The variation of the external work $W$ gives the external force vector of the
element in the absence of distributed pressure load:

$\displaystyle\delta
W=\delta\overline{\boldsymbol{\mathrm{W}}}\,\boldsymbol{\mathrm{T}}^{\rm{T}}\boldsymbol{\mathrm{Q}}^{*}=\delta\overline{\boldsymbol{\mathrm{W}}}\,\boldsymbol{\mathrm{f}}_{\rm{ext}},\quad\Rightarrow\quad\boldsymbol{\mathrm{f}}_{\rm{ext}}=\boldsymbol{\mathrm{T}}^{\rm{T}}\boldsymbol{\mathrm{Q}}^{*}$
(34)

###  2.3 Formation of shell element

To form a flat shell element that considers membrane deformation and bending,
we superimpose a linear membrane element on top of the plate element developed
in the previous section. The displacement fields $u$, $v$ are defined in terms
of area coordinates $L_{1}$, $L_{2}$, $L_{3}$:

$\displaystyle u$ $\displaystyle=u_{1}L_{1}+u_{2}L_{2}+u_{3}L_{3}$ (35)
$\displaystyle v$ $\displaystyle=v_{1}L_{1}+v_{2}L_{2}+v_{3}L_{3}$ (36)

where $u_{i}$ and $v_{i}$ represent the displacements along $x$ and $y$ for
node $i$, respectively. The DoFs vector of the membrane element
$\overline{\boldsymbol{\mathrm{U}}}$ is:

$\displaystyle\overline{\boldsymbol{\mathrm{U}}}=\left\\{u_{1},\,v_{1},\,u_{2},\,v_{2},\,u_{3},\,v_{3}\right\\}^{\rm{T}}$
(37)

The stiffness matrix of the triangular membrane element is given by the
expression:

$\displaystyle\boldsymbol{\mathrm{K}}_{\mathrm{mem}}=\iint_{\Delta}\boldsymbol{\mathrm{B}}_{\mathrm{mem}}^{\rm{T}}\,\boldsymbol{\mathrm{D}}_{\mathrm{mem}}\,\boldsymbol{\mathrm{B}}_{\mathrm{mem}}\,t\,\mathrm{d}x\,\mathrm{d}y$
(38)

where $t$ is the thickness of the membrane element and
$\boldsymbol{\mathrm{B}}_{\mathrm{mem}}$ matrix in equation 38 is

$\displaystyle\boldsymbol{\mathrm{B}}_{\mathrm{mem}}=\frac{1}{2A}\begin{bmatrix}b_{1}&0&b_{2}&0&b_{3}&0\\\\[5.0pt]
0&c_{1}&0&c_{2}&0&c_{3}\\\\[5.0pt]
c_{1}&b_{1}&c_{2}&b_{2}&c_{3}&b_{3}\end{bmatrix}$ (39)

$A$ is the area of the membrane element and:

$\displaystyle b_{1}$ $\displaystyle=y_{2}-y_{3},\quad b_{2}=y_{3}-y_{1},\quad
b_{3}=y_{1}-y_{2}$ (40) $\displaystyle c_{1}$ $\displaystyle=x_{3}-x_{2},\quad
c_{2}=x_{1}-x_{3},\quad c_{3}=x_{2}-x_{1}$ (41)

The $\boldsymbol{\mathrm{D}}^{\mathrm{mem}}$ matrix in equation 38 is:

$\displaystyle\boldsymbol{\mathrm{D}}_{\mathrm{mem}}=\begin{bmatrix}\frac{E_{1}}{1-\nu_{12}\nu_{21}}\,&\frac{\nu_{12}E_{2}}{1-\nu_{12}\nu_{21}}\,&0\\\\[11.99998pt]
\,&\frac{E_{2}}{1-\nu_{12}\nu_{21}}\,&0\\\\[11.99998pt]
\text{symmetric}\,&\,&G_{12}\end{bmatrix}$ (42)

From equations 33 and 38, we have obtained the stiffness matrices of the plate
and membrane elements, respectively. They are assembled to form the stiffness
matrix of the shell element:

$\displaystyle\boldsymbol{\mathrm{K}}_{\mathrm{shell}}=\begin{bmatrix}\boldsymbol{\mathrm{K}}_{\mathrm{mem}}\,&0\\\\[15.00002pt]
0\,&\boldsymbol{\mathrm{K}}_{\mathrm{plate}}\end{bmatrix}$ (43)

The corresponding DoF vector of the shell element is simply:

$\displaystyle\boldsymbol{\mathrm{q}}_{\mathrm{shell}}$
$\displaystyle=\Bigl{\\{}\overline{\boldsymbol{\mathrm{U}}},\overline{\boldsymbol{\mathrm{W}}}\Bigl{\\}}$
(44)
$\displaystyle=\Bigl{\\{}u_{1},v_{1},u_{2},v_{2},u_{3},v_{3},{w}_{1},\frac{\partial\overline{w}_{1}}{\partial
x},\frac{\partial\overline{w}_{1}}{\partial
y},\overline{w}_{2},\frac{\partial\overline{w}_{2}}{\partial
x},\frac{\partial\overline{w}_{3}}{\partial
y},\overline{w}_{3},\frac{\partial\overline{w}_{3}}{\partial
x},\frac{\partial\overline{w}_{3}}{\partial y}\Bigl{\\}}^{\rm{T}}$

The contribution of the shell element to the overall force residual of the
finite element system equation is
$\boldsymbol{\mathrm{f}}_{\rm{ext}}-\boldsymbol{\mathrm{f}}_{\rm{int}}$. When
the applied loads are nodal forces, $\boldsymbol{\mathrm{f}}_{\rm{ext}}$ does
not need to be calculated explicitly in the element subroutine as the nodal
forces can be directly entered in the global external force vector of the
system. In this case, the residual contribution of this element can be
expressed as:

$\displaystyle\boldsymbol{\mathrm{f}}_{\mathrm{res}}=-\boldsymbol{\mathrm{f}}_{\rm{int}}=-\boldsymbol{\mathrm{K}}_{\text{shell}}\,\boldsymbol{\mathrm{q}}_{\mathrm{shell}}$
(45)

###  2.4 Structural CE between two plies

####  2.4.1 DoFs and opening vector

The structural CEs must be compatible with the top and bottom shell elements
kinematically. Thus, they need to share the same DoFs and displacement
interpolations along the interface. The DoF vector of the CE can be defined
as:

$\displaystyle\boldsymbol{\mathrm{q}}_{\rm{CE}}$
$\displaystyle=\left[\overline{\boldsymbol{\mathrm{U}}}_{\text{bot}}^{\rm{T}},\overline{\boldsymbol{\mathrm{W}}}_{\text{bot}}^{\rm{T}},\overline{\boldsymbol{\mathrm{U}}}_{\text{top}}^{\rm{T}},\overline{\boldsymbol{\mathrm{W}}}_{\text{top}}^{\rm{T}}\right]^{\rm{T}}$
(46)

where
$\overline{\boldsymbol{\mathrm{U}}}_{\text{bot/top}},\overline{\boldsymbol{\mathrm{W}}}_{\text{bot/top}}$
are the membrane and plate DoFs of the bottom/top ply, respectively.

Figure 4: Mode I opening

The Mode I opening displacement of the CE can be expressed as

$\displaystyle\Delta_{\text{I}}=w^{\text{topCE}}-w^{\text{botCE}}$ (47)

where $w^{\text{topCE}}$ and $w^{\text{botCE}}$ are the vertical displacements
of the top and bottom CE surfaces, respectively. From Figure 4, using
classical plate theory, we have:

$\displaystyle
w^{\text{topCE}}=w^{\text{top}}+\frac{h^{\text{top}}}{2}(1-\cos{\theta^{\text{top}}})$
(48) $\displaystyle
w^{\text{botCE}}=w^{\text{bot}}-\frac{h^{\text{bot}}}{2}(1-\cos{\theta^{\text{bot}}})$
(49)

$h^{\text{top}}$ and $h^{\text{bot}}$ are the thickness of top and bottom
plies, respectively. $w^{\text{top}}$ and $w^{\text{bot}}$ are the vertical
displacements of the neutral planes of top and bottom plies, respectively. In
this work, only geometrical linearity is considered. Thus, the rotations are
small such that $1-\cos{\theta}\approx 0$. The mode I opening is then simply
the relative displacement between the mid-planes of the plies:

$\displaystyle\Delta_{\text{I}}=w^{\text{top}}-w^{\text{bot}}$ (50)

Figure 5: Mode II opening

The Mode II opening of the CE shown in Figure 5 can be expressed as:

$\displaystyle\Delta_{\text{II}}=u^{\text{topCE}}-u^{\text{botCE}}$ (51)

where $u^{\text{topCE}}$ and $u^{\text{botCE}}$ are the displacements along
the x-axis of the top and bottom CE surfaces, respectively. Considering the
rotations of the shells’ neutral planes and the offsets of shell surfaces from
the neutral planes, they can be written as:

$\displaystyle
u^{\text{topCE}}=u^{\text{top}}+\frac{h^{\text{top}}}{2}\sin{\theta^{\text{top}}_{x}}$
(52) $\displaystyle
u^{\text{botCE}}=u^{\text{bot}}-\frac{h^{\text{bot}}}{2}\sin{\theta^{\text{bot}}_{x}}$
(53)

where $u^{\text{top}}$, $u^{\text{bot}}$ are the displacement $u$ of the
neutral planes of the top and bottom shells, respectively. For a small
rotation $\theta$ in geometrically linear problems,
$\sin{\theta}\approx\theta$. Then, the opening in mode II is:

$\displaystyle\Delta_{\text{II}}=u^{\text{top}}-u^{\text{bot}}+\frac{h^{\text{top}}}{2}\theta^{\text{top}}_{x}+\frac{h^{\text{bot}}}{2}\theta^{\text{bot}}_{x}$
(54)

The mode III opening is found from analogous kinematics in the yz-plane.
Hence:

$\displaystyle\Delta_{\text{III}}=v^{\text{top}}-v^{\text{bot}}+\frac{h^{\text{top}}}{2}\theta^{\text{top}}_{y}+\frac{h^{\text{bot}}}{2}\theta^{\text{bot}}_{y}$
(55)

The central task of CE formulation is to find the matrix,
$\boldsymbol{\mathrm{B}}_{\rm{CE}}$, that relates the opening vector to its
nodal DoFs:

$\displaystyle\boldsymbol{\mathrm{\Delta}}=\left[\Delta_{\text{I}},\Delta_{\text{II}},\Delta_{\text{III}}\right]^{\rm{T}}=\boldsymbol{\mathrm{B}}_{\rm{CE}}\,\boldsymbol{\mathrm{q}}_{\rm{CE}}$
(56)

Examining the expressions of
$\Delta_{\text{I}},\Delta_{\text{II}},\Delta_{\text{III}}$ in equations 50,
54, and 55, we can see that the $\boldsymbol{\mathrm{B}}_{\rm{CE}}$ matrix
shall be composed of sub-matrices that relate the following terms to the nodal
DoFs in $\boldsymbol{\mathrm{q}}_{\rm{CE}}$:

$\displaystyle
w^{\text{top}},\,w^{\text{bot}},\,u^{\text{top}},\,u^{\text{bot}},\,v^{\text{top}},\,v^{\text{bot}},\,\theta^{\text{top}}_{x},\,\theta^{\text{bot}}_{x},\,\theta^{\text{top}}_{y},\,\theta^{\text{bot}}_{y}$
(57)

In equation 35, the in-plane displacements $u$ and $v$ are already expressed
in the nodal membrane DoFs. However, the out-of-plane displacements $w$ and
the rotations $\theta$ in the above list remain to be explicitly expressed in
terms of nodal DoFs in $\boldsymbol{\mathrm{q}}_{\rm{CE}}$.

####  2.4.2 Shape functions of $w$

The displacement $w$ shall be expressed in terms of plate DoFs
$\overline{\boldsymbol{\mathrm{W}}}$ in this section. Referring back to
equation 13, the displacement $w$ is defined in the local Cartesian
coordinates ($x$ and $y$) by the coefficients $A_{i}$ and the vector
$\boldsymbol{\mathrm{\alpha}}$. The expression of $w$ is then divided into two
parts. The first part only contains $A_{i}$, and the second part only contains
$\alpha_{j}$:

$\displaystyle w=w_{A}+w_{\alpha}$ (58)

where

$\displaystyle w_{A}$ $\displaystyle=A_{1}+A_{2}\,x+A_{3}\,y$ (59)
$\displaystyle w_{\alpha}$
$\displaystyle=\alpha_{1}\,x^{2}+\alpha_{2}\,xy+\alpha_{3}\,y^{2}+\alpha_{4}\,x^{3}+\alpha_{5}\,x^{2}y+\alpha_{6}\,xy^{2}+\alpha_{7}\,y^{3}$
(60)

$w_{A}$ can be rewritten as:

$\displaystyle
w_{A}=\boldsymbol{\mathrm{S}}^{\rm{T}}\,\boldsymbol{\mathrm{A}},\quad\boldsymbol{\mathrm{S}}=[1,\,x,\,y]^{\rm{T}},\;\boldsymbol{\mathrm{A}}=[A_{1},\,A_{2},\,A_{3}]^{\rm{T}}$
(61)

Similarly, $w_{\alpha}$ can be rewritten as:

$\displaystyle
w_{\alpha}=\boldsymbol{\mathrm{R}}^{\rm{T}}\,\boldsymbol{\mathrm{\alpha}},\quad\boldsymbol{\mathrm{R}}=[x^{2},\,xy,\,y^{2},\,x^{3},\,x^{2}y,\,xy^{2},\,y^{3}]^{\rm{T}}$
(62)

Let us firstly look at $w_{\alpha}$. Based on equation 30,
$\boldsymbol{\mathrm{\alpha}}$ and $\overline{\boldsymbol{\mathrm{W}}}$ can be
related by multiplying the three matrices $\boldsymbol{\mathrm{H}}$,
$\boldsymbol{\mathrm{B}}$ and $\boldsymbol{\mathrm{T}}$ in section 2.2.4.
Here, for simplicity, we replace the product of these three matrices with
matrix $\boldsymbol{\mathrm{C}}$:

$\displaystyle\boldsymbol{\mathrm{C}}=\boldsymbol{\mathrm{H}}^{-1}(\boldsymbol{\mathrm{B}}\boldsymbol{\mathrm{T}})$
(63)

Thus, the vector $\boldsymbol{\mathrm{\alpha}}$ can be rewritten as:

$\displaystyle\boldsymbol{\mathrm{\alpha}}=\boldsymbol{\mathrm{C}}\overline{\boldsymbol{\mathrm{W}}}$
(64)

Substituting equation 64 in to 62, we can obtain the expression of
$w_{\alpha}$ in terms of $\overline{\boldsymbol{\mathrm{W}}}$:

$\displaystyle
w_{\alpha}=\boldsymbol{\mathrm{R}}^{\rm{T}}\boldsymbol{\mathrm{C}}\overline{\boldsymbol{\mathrm{W}}}$
(65)

Next, we move on to express $w_{A}$ in terms of
$\overline{\boldsymbol{\mathrm{W}}}$. From equations 58, 61 and 65, we have:

$\displaystyle
w_{A}=\boldsymbol{\mathrm{S}}^{\rm{T}}\,\boldsymbol{\mathrm{A}}=w-\boldsymbol{\mathrm{R}}^{\rm{T}}\boldsymbol{\mathrm{C}}\overline{\boldsymbol{\mathrm{W}}}$
(66)

Evaluating $w_{A}$ at the three nodes, we can obtain:

$\displaystyle
w_{A}(x_{i},y_{i})=\boldsymbol{\mathrm{S}}^{\rm{T}}(x_{i},y_{i})\,\boldsymbol{\mathrm{A}}=\overline{w}_{i}-\boldsymbol{\mathrm{R}}^{\rm{T}}(x_{i},y_{i})\,\boldsymbol{\mathrm{C}}\overline{\boldsymbol{\mathrm{W}}},\quad
i=1,2,3$ (67)

which give us the following matrix equation:

$\displaystyle\underbrace{\begin{bmatrix}1\,&x_{1}\,&y_{1}\\\\[5.0pt]
1\,&x_{2}\,&y_{2}\\\\[5.0pt]
1\,&x_{3}\,&y_{3}\end{bmatrix}}_{\boldsymbol{\mathrm{M}}_{A}}\boldsymbol{\mathrm{A}}=\begin{bmatrix}\overline{w}_{1}\\\\[5.0pt]
\overline{w}_{2}\\\\[5.0pt]
\overline{w}_{3}\end{bmatrix}-\underbrace{\begin{bmatrix}x_{1}^{2}&x_{1}y_{1}&y_{1}^{2}&x_{1}^{3}&x_{1}^{2}y_{1}&x_{1}y_{1}^{2}&y_{1}^{3}\\\\[5.0pt]
x_{2}^{2}&x_{2}y_{2}&y_{2}^{2}&x_{2}^{3}&x_{2}^{2}y_{2}&x_{2}y_{2}^{2}&y_{2}^{3}\\\\[5.0pt]
x_{3}^{2}&x_{3}y_{3}&y_{3}^{2}&x_{3}^{3}&x_{3}^{2}3_{2}&x_{3}y_{3}^{2}&y_{3}^{3}\end{bmatrix}}_{\boldsymbol{\mathrm{M}}_{\alpha}}\boldsymbol{\mathrm{C}}\overline{\boldsymbol{\mathrm{W}}}$
(68)

Since

$\displaystyle\begin{bmatrix}\overline{w}_{1}\\\\[5.0pt]
\overline{w}_{2}\\\\[5.0pt]
\overline{w}_{3}\end{bmatrix}=\underbrace{\begin{bmatrix}1&0&0&0&0&0&0&0&0\\\\[5.0pt]
0&0&0&1&0&0&0&0&0\\\\[5.0pt]
0&0&0&0&0&0&1&0&0\end{bmatrix}}_{\boldsymbol{\mathrm{B}}_{A}}\overline{\boldsymbol{\mathrm{W}}}$
(69)

we can obtain $\boldsymbol{\mathrm{A}}$ as:

$\displaystyle\boldsymbol{\mathrm{A}}=\boldsymbol{\mathrm{M}}_{A}^{-1}\left(\boldsymbol{\mathrm{B}}_{A}-\boldsymbol{\mathrm{M}}_{\alpha}\,\boldsymbol{\mathrm{C}}\right)\overline{\boldsymbol{\mathrm{W}}}$
(70)

Using equation 61, we obtain:

$\displaystyle
w_{A}=\boldsymbol{\mathrm{S}}^{\rm{T}}\boldsymbol{\mathrm{M}}_{A}^{-1}\left(\boldsymbol{\mathrm{B}}_{A}-\boldsymbol{\mathrm{M}}_{\alpha}\,\boldsymbol{\mathrm{C}}\right)\overline{\boldsymbol{\mathrm{W}}}$
(71)

Combining with equation 65, we obtain the expression of $w$:

$\displaystyle w$
$\displaystyle=\boldsymbol{\mathrm{S}}^{\rm{T}}\boldsymbol{\mathrm{M}}_{A}^{-1}\left(\boldsymbol{\mathrm{B}}_{A}-\boldsymbol{\mathrm{M}}_{\alpha}\,\boldsymbol{\mathrm{C}}\right)\overline{\boldsymbol{\mathrm{W}}}+\boldsymbol{\mathrm{R}}^{\rm{T}}\boldsymbol{\mathrm{C}}\overline{\boldsymbol{\mathrm{W}}}=\boldsymbol{\mathrm{N}}_{w}\,\overline{\boldsymbol{\mathrm{W}}}$
(72)

where the matrix $\boldsymbol{\mathrm{N}}_{w}$ ($1\times 9$) can be extracted
as:

$\displaystyle\boldsymbol{\mathrm{N}}_{w}=\boldsymbol{\mathrm{S}}^{\rm{T}}\boldsymbol{\mathrm{M}}_{A}^{-1}(\boldsymbol{\mathrm{B}}_{A}-\boldsymbol{\mathrm{M}}_{\alpha}\boldsymbol{\mathrm{C}})+\boldsymbol{\mathrm{R}}^{\rm{T}}\boldsymbol{\mathrm{C}}$
(73)

####  2.4.3 Shape functions of $\theta$

Using the cubic polynomial expression of $w$ from equation 13, $\theta_{x}$
can be expressed as:

$\displaystyle\theta_{x}=\frac{\partial w}{\partial
x}=A_{2}+2\alpha_{1}x+\alpha_{2}y+3\alpha_{4}x^{2}+2\alpha_{5}xy+\alpha_{6}y^{2}$
(74)

Similarly as in the case of $w$, the expression of the above rotations can
also be divided into two parts, one containing the $A$ coefficients only and
another the $\alpha$ coefficients only:

$\displaystyle\theta_{x}=A_{2}+\theta_{x\alpha}$ (75)

Since $A_{2}$ is part of the vector $\boldsymbol{\mathrm{A}}$ which has been
expressed in equation 70, we have:

$\displaystyle
A_{2}=\boldsymbol{\mathrm{B}}_{x}\boldsymbol{\mathrm{A}}=\boldsymbol{\mathrm{B}}_{x}\boldsymbol{\mathrm{M}}_{A}^{-1}(\boldsymbol{\mathrm{B}}_{A}-\boldsymbol{\mathrm{M}}_{\alpha}\boldsymbol{\mathrm{C}})\overline{\boldsymbol{\mathrm{W}}}$
(76)

where $\boldsymbol{\mathrm{B}}_{x}$ is a Boolean matrix:

$\displaystyle\boldsymbol{\mathrm{B}}_{x}=\begin{bmatrix}0&1&0\end{bmatrix}$
(77)

The expression $\theta_{x\alpha}$ can be written as:

$\displaystyle\theta_{x\alpha}=\boldsymbol{\mathrm{R}}_{x}^{\rm{T}}\boldsymbol{\mathrm{\alpha}}=\boldsymbol{\mathrm{R}}_{x}^{\rm{T}}\boldsymbol{\mathrm{C}}\overline{\boldsymbol{\mathrm{W}}}$
(78)

where equation 64 is used to express $\boldsymbol{\mathrm{\alpha}}$ and
$\boldsymbol{\mathrm{R}}_{x}$ is also a vector of local coordinates:

$\displaystyle\boldsymbol{\mathrm{R}}_{x}=\begin{bmatrix}2x&y&0&3x^{2}&2xy&y^{2}&0\end{bmatrix}^{\rm{T}}$
(79)

Substituting the expressions above into equation 74, the expression of
$\theta_{x}$ is

$\displaystyle\theta_{x}$
$\displaystyle=\boldsymbol{\mathrm{B}}_{x}\boldsymbol{\mathrm{M}}_{A}^{-1}(\boldsymbol{\mathrm{B}}_{A}-\boldsymbol{\mathrm{M}}_{\alpha}\boldsymbol{\mathrm{C}})\overline{\boldsymbol{\mathrm{W}}}+\boldsymbol{\mathrm{R}}_{x}^{\rm{T}}\boldsymbol{\mathrm{C}}\overline{\boldsymbol{\mathrm{W}}}$
(80)

The shape function $\boldsymbol{\mathrm{N}}_{x}$ of rotation $\theta_{x}$ is
therefore:

$\displaystyle\boldsymbol{\mathrm{N}}_{x}=\boldsymbol{\mathrm{B}}_{x}\boldsymbol{\mathrm{M}}_{A}^{-1}(\boldsymbol{\mathrm{B}}_{A}-\boldsymbol{\mathrm{M}}_{\alpha}\boldsymbol{\mathrm{C}})+\boldsymbol{\mathrm{R}}_{x}^{\rm{T}}\boldsymbol{\mathrm{C}}$
(81)

The derivation of shape function $\boldsymbol{\mathrm{N}}_{y}$ for rotation
$\theta_{y}$ is very similar to that for $\theta_{x}$ and is omitted here for
brevity.

####  2.4.4 $\boldsymbol{\mathrm{B}}_{\rm{CE}}$ matrix

Once we have determined the shape functions, the next step is to assemble the
$\boldsymbol{\mathrm{B}}_{\rm{CE}}$ matrix that relates the DoFs vector
$\boldsymbol{\mathrm{q}}_{\rm{CE}}$ to the opening vector
$\boldsymbol{\mathrm{\Delta}}$ (c.f. equation 56). Substituting the shape
functions in equations 73 and 81 to express the $w$ and $\theta$ terms in
equations 50, 54, and 55, the expression of
$\boldsymbol{\mathrm{B}}_{\rm{CE}}$ can be obtained as the following:

(82) (83) (84)

####  2.4.5 $\boldsymbol{\mathrm{D}}_{{\rm{CE}}}$ matrix

The relationship between traction $\boldsymbol{\mathrm{\tau}}$ and opening
vector $\boldsymbol{\mathrm{\Delta}}$ is expressed through the constitutive
matrix $\boldsymbol{\mathrm{D}}_{\rm{CE}}$:

$\displaystyle\boldsymbol{\mathrm{\tau}}=\boldsymbol{\mathrm{D}}_{\rm{CE}}\,\boldsymbol{\mathrm{\Delta}},\quad\boldsymbol{\mathrm{D}}_{\rm{CE}}=\begin{bmatrix}(1-d_{\text{I}})K&0&0\\\\[5.0pt]
0&(1-d)K&0\\\\[5.0pt] 0&0&(1-d)K\end{bmatrix}$ (85)

where $K$ and $d$ are the penalty stiffness and the damage variable of the CE,
respectively. The damage variable $d$ in this work is updated by the bi-linear
cohesive law proposed by Turon et al. [32, 33]. The damage variable
$d_{\text{I}}$ under Mode I loading is distinguished from $d$ to avoid
interpenetration of the top and bottom surfaces under compression:

$\displaystyle d_{\text{I}}=\begin{cases}d,\quad&\Delta_{\text{I}}\geq
0\\\\[5.0pt] 0,\quad&\Delta_{\text{I}}<0\end{cases}$ (86)

The penalty stiffness $K$ is set as:

$\displaystyle K=\alpha\frac{E_{3}}{t}$ (87)

where $E_{3}$ is the out-of-plane laminate Young’s modulus, $t$ is the
thickness of the laminate and $\alpha$ is a constant, here set to be 50 [8].

####  2.4.6 Stiffness matrix and residual vector

In this work, the secant stiffness matrix is used as the stiffness matrix of
the CE:

$\displaystyle\boldsymbol{\mathrm{K}}_{\rm{CE}}=\int_{\Gamma}\boldsymbol{\mathrm{B}}_{\rm{CE}}^{\rm{T}}\,\boldsymbol{\mathrm{D}}_{\rm{CE}}\,\boldsymbol{\mathrm{B}}_{\rm{CE}}\,\mathrm{d}\,\Gamma$
(88)

where $\Gamma$ represents the domain of the interface. The integral in
equation 88 is hard to calculate analytically. Thus, the Gaussian integration
scheme is applied to obtain the stiffness matrix numerically. Earlier works
have shown that using a higher number of quadrature points improves the
accuracy and smoothness of the load-displacement solutions of delamination
simulations [25, 15, 28]. Therefore, thirteen quadrature points are used in
this work for the integration of the structural CE, with their coordinates and
weights taken from the work of Cowper [34].

Assuming that no external distributed loads are applied to the cohesive
interfaces, only the internal force vector of the CE contributes to the
overall residual vector of the system. The residual vector contribution from
this CE can then be written as:

$\displaystyle\boldsymbol{\mathrm{f}}_{\mathrm{res}}=-\boldsymbol{\mathrm{f}}_{\rm{int}}=-\boldsymbol{\mathrm{K}}_{\rm{CE}}\,\boldsymbol{\mathrm{q}}_{\rm{CE}}$
(89)

##  3 Results

The structural elements section 2 have been implemented in the Abaqus user-
defined element subroutines. The structural CE model was firstly verified on
the three classical benchmarks, namely the double cantilever beam (DCB), the
end-notched flexure (ENF), and the mixed-mode bending (MMB) problems. As the
above benchmarks only pertain to unidirectional laminates, a multi-directional
laminate problem, namely the single-leg bending (SLB) problem [35, 36], was
simulated to demonstrate the 3D capacity of the model. All simulations of the
structural CE model were performed with the Quasi-Newton solver without
damping or viscosity. The reference solutions on these problems were obtained
with analytical equations and Abaqus standard FE simulations using 8-node
linear solid elements (C3D8I element in Abaqus) for the plies and 8-node
linear CEs for the interfaces. The linear CEs were implemented as user-defined
elements following the formulation in [5], with the cohesive law from Turon et
al. [32, 33]. The problem descriptions and modelling details will be presented
in this section, followed by comparisons of load-displacement curves and
computational time.

###  3.1 Unidirectional laminate benchmarks: DCB, ENF, and MMB

####  3.1.1 Description of the unidirectional laminate tests

(a) DCB

(b) ENF

(c) MMB

Figure 6: DCB, ENF, and MMB test specimens

The unidirectional laminate benchmarks are drawn from the work by Krueger
[36], except that the pre-crack length of the ENF model is increased to 35 mm
in this work to avoid the snapback in load-displacement response. The
geometric parameters and boundary conditions are shown in Figure 6 and Table
1. The detailed material properties are shown in Table 2 and 3.

Table 1: Geometric parameters for unidirectional benchmarks Model Param. (mm) | $2L$ | $a_{0}$ | $h$ | $b$ (width) | $c$
---|---|---|---|---|---
DCB | 150.0 | 30.5 | 1.50 | 25.0 | -
ENF | 101.6 | 35.0 | 2.25 | 25.4 | -
MMB | 100.8 | 25.4 | 2.25 | 25.4 | 41.3
Table 2: Material properties for DCB [36] T300/1076 Unidirectional
graphite/epoxy prepreg
---
$E_{11}$= 139.4 GPa |  $E_{22}$= 10.16 GPa |  $E_{33}$ = 10.16 GPa
$\nu_{12}$ = 0.30 |  $\nu_{13}$ = 0.30 |  $\nu_{23}$ = 0.436
$G_{12}$ = 4.6 GPa |  $G_{13}$ = 4.6 GPa |  $G_{23}$ = 3.54 GPa
Fracture toughness data |  |
$G_{\rm{Ic}}=0.170\,\mathrm{kJ/m^{2}}$ | $G_{\rm{IIc}}=0.494\,\mathrm{kJ/m^{2}}$ | $\eta=1.62$
Interfacial strength data [33, 12] |  |
$\tau_{\rm{Ic}}=30\,\mathrm{MPa}$ | $\tau_{\rm{IIc}}=60\,\mathrm{MPa}$ |
Table 3: Material properties for ENF and MMB [36] IM7/8552 Unidirectional
graphite/epoxy prepreg
---
$E_{11}$= 161 GPa |  $E_{22}$= 11.38 GPa |  $E_{33}$ = 11.38 GPa
$\nu_{12}$ = 0.32 |  $\nu_{13}$ = 0.32 |  $\nu_{23}$ = 0.45
$G_{12}$ = 5.2 GPa |  $G_{13}$ = 5.2 GPa |  $G_{23}$ = 3.9 GPa
Fracture toughness data |  |
$G_{\rm{Ic}}=0.212\,\mathrm{kJ/m^{2}}$ | $G_{\rm{IIc}}=0.774\,\mathrm{kJ/m^{2}}$ | $\eta=2.1$
Interfacial strength data [33, 12] |  |
$\tau_{\rm{Ic}}=30\,\mathrm{MPa}$ | $\tau_{\rm{IIc}}=60\,\mathrm{MPa}$ |

####  3.1.2 Description of the model

Since the plies of the benchmarks in this section are all $0\degree$, the
symmetric lay-up conditions in equation 2.2.2 are met. Therefore, the
structural model could simply use one layer of shell elements on each side of
the delamination, with one layer of structural CEs in between. The Abaqus
solid models are built according to the work of Krueger [36], except that CEs,
instead of VCCT, are used to model delamination.

####  3.1.3 Load-displacement curves

(a) DCB-Solid

(b) DCB-Structural

(c) ENF-Solid

(d) ENF-Structural

(e) MMB-Solid

(f) MMB-Structural

Figure 7: Results of the DCB, ENF, and MMB tests

The load-displacement curves obtained from the simulations are shown in Figure
7. The designations “-solid” and “-structural” denote the results of the
traditional solid element model and those of the proposed structural CE model,
respectively. The results on meshes of different element sizes are plotted
together with the analytical solutions [37, 38, 39].

The DCB results in Figure 7 indicated that the standard model with solid
elements requires that the element size should not exceed 0.5 mm. When the
element sizes are greater than or equal to 1 mm, the simulation results differ
greatly from the analytical solution. The error on peak load exceeds 30%, and
the post-peak curve also stays way above the analytical one. With the proposed
structural model, the results on 2-mm mesh remain in close agreement with the
analytical solution, with a less than 3% error on the peak load. The post-peak
curve also closely follows the analytical curve. The slight over-prediction of
the stiffness is expected as a result of neglecting transverse shear in the
Kirchhoff-Love shell elements. The analytical solution based on the corrected
beam theory, however, includes such transverse shear effect [37]. Even on a
5-mm mesh, the structural model can predict the peak load fairly well despite
the post-peak oscillations right above the analytical curve. Such oscillations
are due to the larger spacing between the integration points on coarser
meshes.

In the ENF case, the solid element model on 1-mm mesh could capture the peak
load correctly, thanks to the larger cohesive zone in Mode II delamination
than in Mode I. However, the result on 2.5-mm mesh already shows a clear drift
towards over-prediction. As the element size increases to 5 mm, the peak load
and post-peak response again become severely over-predicted. Correspondingly,
if the structural model is used on the 5-mm mesh, the predicted curve remains
close to the analytical solution throughout the loading history. Even the
7.5-mm structural model manages to capture the load-displacement response very
accurately. The slight under-prediction of the peak load is expected, as the
analytical curve is based on Linear Elastic Fracture Mechanics, which ignores
the presence of material softening (i.e. the cohesive zone) at the crack tip.

In the MMB case, very similar trends can be observed. The solid model on the
2.5-mm mesh cannot capture the correct load-displacement response, while the
structural model’s predictions on the 2.5-mm and 5-mm meshes oscillates
closely around the analytical curve. Even the 7.5-mm structural model predicts
the peak load correctly, albeit with bigger oscillations during the load drop
section due to the coarser distribution of integration points in larger
elements.

####  3.1.4 Computational performances

Table 4: Comparison of CPU time (unit: second) | DCB | ENF | MMB
---|---|---|---
Solid model (mesh size) | 5311.8 (0.5mm) | 4239.7 (1mm) | 5560.7 (0.5mm)
Structural model (mesh size) | 197.15 (5mm) | 285.33 (7.5mm) | 370.02 (7.5mm)
| 836.76 (2mm) | 736.41 (5mm) | 474.79 (5mm)
Reduction by structural | 96.3% | 93.3% | 93.3%
| 84.2% | 82.6% | 91.5%

By allowing larger elements to be used, the proposed structural model is able
to reduce the computational time of the delamination simulations considerably.
The results of the structural model are compared against those of the solid
element model. The comparison on CPU time is reported in Table 4. It can be
seen that the structural model can reduce the computational time in all three
problems by more than 90%, while retaining accurate predictions of the peak
loads and the overall load-drop curves.

###  3.2 Multi-directional laminate benchmark: SLB

####  3.2.1 Description of the SLB test

The SLB specimen is shown in Figure 8, with its geometrical parameters
specified in Table 5. Unlike in the previous benchmarks, the ply angles here
are no longer all $0\degree$. The material of the SLB model is C12K/R6376 and
its properties are shown in Table 6. The value of $\tau_{\rm{IIc}}$ cannot be
found in the literature. According to the review in Ref. [12], the ratio of
$\tau_{\rm{IIc}}$ to $\tau_{\rm{Ic}}$ ranges from 1.25 to 2. Hence, the value
of $\tau_{\rm{IIc}}$ here is estimated to be 38 MPa, 1.5 times that of
$\tau_{\rm{Ic}}$.

Figure 8: SLB specimen [35] Table 5: Geometric parameters for the SLB specimen Parameter | $2L$ | $a_{0}$ | $t_{1}$ | $t_{2}$ |  $b$ (width)
---|---|---|---|---|---
value (mm) | 177.8 | 60 | 2 | 2 | 25.4
Table 6: C12K/R6376 material properties for SLB specimen [36] $E_{11}$= 146.9 GPa |  $E_{22}$= 10.6 GPa |  $E_{33}$ = 10.6 GPa
---|---|---
$\nu_{12}$ = 0.33 |  $\nu_{13}$ = 0.33 |  $\nu_{23}$ = 0.33
$G_{12}$ = 5.45 GPa |  $G_{13}$ = 5.45 GPa |  $G_{23}$ = 3.99 GPa
Fracture toughness data |  |
$G_{\rm{Ic}}=0.34\,\mathrm{kJ/m^{2}}$ | $G_{\rm{IIc}}=1.286\,\mathrm{kJ/m^{2}}$ | $\eta=3.39$
Interfacial strength data |  |
$\tau_{\rm{Ic}}=25\,\mathrm{MPa}$ [40] |  $\tau_{\rm{IIc}}=38\,\mathrm{MPa}$* |
|  | * Estimated [12]

####  3.2.2 Description of the model

Figure 9: Layer-wise model for SLB

Since the SLB specimen is a multi-directional laminate and its lay-up is not
symmetric, we can no longer model the entire laminate with a single layer of
shell elements on each side of the delamination, as done in the unidirectional
models. Therefore, the SLB model uses one layer of shell elements for each ply
and one layer of structural CEs between every two plies (c.f., Figure 9).

####  3.2.3 Load-displacement curves

(a) SLB-Solid

(b) SLB-Structural

Figure 10: Load-displacement predictions of the SLB

In Figure 10, the load-displacement curves are reported for the solid and
structural models. Both are compared with the analytical solution derived in
C. All the load-displacement predictions by the solid model could not reach
numerical convergence within the allowable settings on solver iterations. As a
result, the load-decreasing section of the curve could not be obtained. From
the left figure, it can be seen that only the 0.5-mm solid model can capture
the correct peak load. For the 2.5-mm solid model, the peak load is much lower
than the analytical solution due to premature divergence of analysis. However,
the peak load of the model on the 5-mm mesh is much larger than the analytical
solution, in line with the solid model’s performance in the unidirectional
benchmarks. The right figure shows that all the structural models can obtain
converged solutions throughout the full loading history. Overall speaking, the
predicted curves of the structural models follow well the analytical solution.
On the 2.5-mm mesh, the error of peak load is 2%, the same level of accuracy
as that of the unidirectional structural model. On the 5-mm mesh, the
predicted peak load is larger than that of the analytical solution by 12.5%.
However, the post-peak part of the curve quickly converges to that of the
2.5-mm mesh result. It should be noted that the initial stiffness of the
structural models is 5% larger than that of the analytical solution. The
reason is that the compliance $C$ obtained by equation 204 is calculated by
simplifying the SLB model into a 2D beam model (c.f., Figure 1), where the
neutral axis is assumed to be aligned with the centerline of the uncracked
region for simplicity. However, the actual neutral axis of the SLB model
should be slightly above the assumed one, as only the top part of the cracked
region would be under bending. This discrepancy would cause the analytical
stiffness to be slightly lower than that in the actual situation. Another
discrepancy to note is that when the crack length $a$ is longer than $L$, the
load increases again with increasing displacement (c.f., C, curve DFE). For
this part, the predicted curves of the structural models are above the
analytical curve. This discrepancy is however expected. As the crack length
$a$ increases with delamination propagation, the mixed-mode ratio $B$ is found
to also increase, particularly when $a>L$ [41], which causes the critical
energy release rate $G_{c}$ to increase as well. However, the analytical
solution in C assumes a fixed mixed-ratio, which underestimates $G_{c}$ for
$a>L$, hence resulting in the analytical curve being below the numerical ones.
The analytical curve could be corrected, should an analytical expression of
$G_{c}$ in terms of $a$ becomes available.

In addition, it needs to be emphasized that the difficulty of convergence
leads to a significant increase in computational time. Due to this reason,
artificial viscosity would be needed to continue the simulation through the
convergence difficulties. However, this would introduce a fitting parameter
which generally requires trial and error to properly set its value. This
suggests that the proposed structural model not only allows coarser meshes
with faster computations but also ensures better numerical stability than the
solid element model.

##  4 Summary and conclusions

This work aims to establish a state-of-the-art numerical method to simulate
the delamination of composite laminates aimed at overcoming the cohesive zone
limit on mesh density. The Kirchhoff-Love triangular cubic plate element from
Allman [30] has been extended for the modelling of symmetric laminate shells.
A structural CE, conforming with the shell element, has been developed to
model delamination between the shells. The combination of the Kirchhoff-Love
shell element and the structural CE is shown to be a powerful new method that
overcomes the cohesive zone limit and models delamination with high accuracy
and efficiency under different loading conditions. This capability is expected
to make a strong impact in the composites modeling community, as the cohesive
zone limit has been a long-lasting and well-known problem for delamination
modeling.

This method has also been validated on the DCB, ENF, MMB, and SLB problems.
The results show excellent agreement with the analytical solutions. By
comparing the results with those from the solid element model, the structural
model has shown the following advantages: 1) it can accurately predict the
load-displacement curves with significantly coarser meshes than the solid
element model, allowing 5-mm elements to be used instead of the 0.5-mm
elements in the latter; 2) its computational speed is much faster than that of
the solid element model, achieving more than 90% reduction in CPU time; and 3)
it has better numerical stability than the solid element model, reaching
convergence without needing artificial viscosity for stabilization.

Future work includes the incorporation of intralaminar damage in the
structural model. Since Kirchhoff-Love thin shell elements are used to model
the plies, the transverse strains are ignored. This will impact the prediction
of intralaminar damage in the model because an accurate three-dimensional
damage criterion would require out-of-plane stresses and/or strains as inputs.
However, since the layer-by-layer modelling method used in this work employs
CEs between every two plies, the required out-of-plane stresses can be
obtained by extracting the tractions of the CEs and interpolating them between
the layers [15]. This is the subject of an ongoing work in the group.

## Acknowledgements

The authors would like to acknowledge the useful discussions with Giorgio
Tosti Balducci, PhD candidate in the same group. The first author would like
to acknowledge funding support from the China Scholarship Council
(No.201906290034) for this research.

## References

  * Barenblatt [1962] G. I. Barenblatt,  The mathematical theory of equilibrium cracks in brittle fracture,  Advances in applied mechanics 7 (1962) 55–129.
  * Dugdale [1960] D. S. Dugdale,  Yielding of steel sheets containing slits,  Journal of the Mechanics and Physics of Solids 8 (1960) 100–104.
  * Qiu et al. [2001] Y. Qiu, M. A. Crisfield, G. Alfano,  An interface element formulation for the simulation of delamination with buckling,  Engineering Fracture Mechanics 68 (2001) 1755–1776.
  * Yang and Cox [2005] Q. Yang, B. Cox,  Cohesive models for damage evolution in laminated composites,  International Journal of Fracture 133 (2005) 107–137.
  * Camanho and Dávila [2002] P. P. Camanho, C. G. Dávila,  Mixed-mode decohesion finite elements for the simulation of delamination in composite materials,  NASA Technical Reports 211737 (2002).
  * Moës and Belytschko [2002] N. Moës, T. Belytschko,  Extended finite element method for cohesive crack growth,  Engineering fracture mechanics 69 (2002) 813–833.
  * Xu and Needleman [1993] X.-P. Xu, A. Needleman,  Void nucleation by inclusion debonding in a crystal matrix,  Modelling and Simulation in Materials Science and engineering 1 (1993) 111.
  * Turon et al. [2007] A. Turon, C. G. Davila, P. P. Camanho, J. Costa,  An engineering solution for mesh size effects in the simulation of delamination using cohesive zone models,  Engineering fracture mechanics 74 (2007) 1665–1682.
  * Pirondi et al. [2014] A. Pirondi, G. Giuliese, F. Moroni, A. Bernasconi, A. Jamil,  Comparative study of cohesive zone and virtual crack closure techniques for three-dimensional fatigue debonding,  The Journal of Adhesion 90 (2014) 457–481.
  * Chen et al. [1999] J. Chen, M. Crisfield, A. Kinloch, E. P. Busso, F. L. Matthews, Y. Qiu,  Predicting progressive delamination of composite material specimens via interface elements,  Mechanics of composite materials and structures 6 (1999) 301–317.
  * Harper and Hallett [2008] P. W. Harper, S. R. Hallett,  Cohesive zone length in numerical simulations of composite delamination,  Engineering Fracture Mechanics 75 (2008) 4774–4792.
  * Lu et al. [2019] X. Lu, M. Ridha, B. Chen, V. Tan, T. Tay,  On cohesive element parameters and delamination modelling,  Engineering Fracture Mechanics 206 (2019) 278–296.
  * Yang et al. [2010] Q. Yang, X. Fang, J. Shi, J. Lua,  An improved cohesive element for shell delamination analyses,  International journal for numerical methods in engineering 83 (2010) 611–641.
  * Do et al. [2013] B. Do, W. Liu, Q. Yang, X. Su,  Improved cohesive stress integration schemes for cohesive zone elements,  Engineering Fracture Mechanics 107 (2013) 14–28.
  * Russo and Chen [2020] R. Russo, B. Chen,  Overcoming the cohesive zone limit in composites delamination: modeling with slender structural elements and higher-order adaptive integration,  International Journal for Numerical Methods in Engineering 121 (2020) 5511–5545.
  * Guiamatsia et al. [2009] I. Guiamatsia, J. Ankersen, G. Davies, L. Iannucci,  Decohesion finite element with enriched basis functions for delamination,  Composites Science and Technology 69 (2009) 2616–2624.
  * Guiamatsia et al. [2010] I. Guiamatsia, G. Davies, J. Ankersen, L. Iannucci,  A framework for cohesive element enrichment,  Composite structures 92 (2010) 454–459.
  * Guiamatsia et al. [2013] I. Guiamatsia, J. Ankersen, L. Iannucci,  A study of mixed-mode composite delamination using enriched interface elements,  The Aeronautical Journal 117 (2013) 959–967.
  * Samimi et al. [2009] M. Samimi, J. Van Dommelen, M. Geers,  An enriched cohesive zone model for delamination in brittle interfaces,  International journal for numerical methods in engineering 80 (2009) 609–630.
  * Samimi et al. [2011a] M. Samimi, J. Van Dommelen, M. Geers,  A self-adaptive finite element approach for simulation of mixed-mode delamination using cohesive zone models,  Engineering Fracture Mechanics 78 (2011a) 2202–2219.
  * Samimi et al. [2011b] M. Samimi, J. Van Dommelen, M. Geers,  A three-dimensional self-adaptive cohesive zone model for interfacial delamination,  Computer Methods in Applied Mechanics and Engineering 200 (2011b) 3540–3553.
  * Van der Meer et al. [2012] F. Van der Meer, N. Moës, L. J. Sluys,  A level set model for delamination–modeling crack growth without cohesive zone or stress singularity,  Engineering Fracture Mechanics 79 (2012) 191–212.
  * Lu et al. [2018] X. Lu, B.-Y. Chen, V. B. Tan, T.-E. Tay,  Adaptive floating node method for modelling cohesive fracture of composite materials,  Engineering Fracture Mechanics 194 (2018) 240–261.
  * Chen et al. [2014] B. Chen, S. Pinho, N. De Carvalho, P. Baiz, T. Tay,  A floating node method for the modelling of discontinuities in composites,  Engineering Fracture Mechanics 127 (2014) 104–134.
  * Álvarez et al. [2014] D. Álvarez, B. Blackman, F. Guild, A. Kinloch,  Mode i fracture in adhesively-bonded joints: A mesh-size independent modelling approach using cohesive elements,  Engineering Fracture Mechanics 115 (2014) 73–95.
  * Mukhopadhyay and Bhatia [2023] S. Mukhopadhyay, S. Bhatia,  Accurate coarse mesh simulation of delamination in composites using a novel hp-adaptive cohesive element,  Journal of Composite Materials 57 (2023) 2201–2218.
  * Daniel et al. [2023] P. M. Daniel, J. Främby, M. Fagerström, P. Maimí,  An efficient err-cohesive method for the modelling of delamination propagation with large elements,  Composites Part A: Applied Science and Manufacturing 167 (2023) 107423\.
  * Tosti Balducci and Chen [2024] G. Tosti Balducci, B. Chen,  Overcoming the cohesive zone limit in the modelling of composites delamination with tuba cohesive elements,  arXiv preprint arXiv:2403.09895 (2024).
  * Bell [1969] K. Bell,  A refined triangular plate bending finite element,  International journal for numerical methods in engineering 1 (1969) 101–122.
  * Allman [1976] D. Allman,  A simple cubic displacement element for plate bending,  International Journal for Numerical Methods in Engineering 10 (1976) 263–281.
  * Przemieniecki [1985] J. S. Przemieniecki, Theory of matrix structural analysis, Courier Corporation, 1985\.
  * Turon et al. [2006] A. Turon, P. P. Camanho, J. Costa, C. Dávila,  A damage model for the simulation of delamination in advanced composites under variable-mode loading,  Mechanics of materials 38 (2006) 1072–1089.
  * Turon et al. [2010] A. Turon, P. Camanho, J. Costa, J. Renart,  Accurate simulation of delamination growth under mixed-mode loading using cohesive elements: Definition of interlaminar strengths and elastic stiffness,  Composite structures 92 (2010) 1857–1864.
  * Cowper [1973] G. Cowper,  Gaussian quadrature formulas for triangles,  International Journal for Numerical Methods in Engineering 7 (1973) 405–408.
  * Davidson et al. [1995] B. Davidson, R. Krüger, M. König,  Three-dimensional analysis of center-delaminated unidirectional and multidirectional single-leg bending specimens,  Composites Science and Technology 54 (1995) 385–394.
  * Krueger [2015] R. Krueger,  A summary of benchmark examples to assess the performance of quasi-static delamination propagation prediction capabilities in finite element codes,  Journal of Composite Materials 49 (2015) 3297–3316.
  * Williams [1989] J. Williams,  The fracture mechanics of delamination tests,  The Journal of strain analysis for engineering design 24 (1989) 207–214.
  * AST [2007] ASTM D5528 /D5528M-21, Standard test method for mode I interlaminar fracture toughness of unidirectional fiber-reinforced polymer matrix composites, ASTM-International, West Conshohocken, PA, 2007.
  * AST [2022] ASTM D6671/D6671M-22, Standard test method for mixed mode I -mode II interlaminar fracture toughness of unidirectional fiber reinforced polymer matrix composites, ASTM-International, West Conshohocken, PA, 2022.
  * Bruyneel et al. [2009] M. Bruyneel, J.-P. Delsemme, P. Jetteur, F. Germain,  Modeling inter-laminar failure in composite structures: illustration on an industrial case study,  Applied Composite Materials 16 (2009) 149–162.
  * Krueger et al. [2020] R. Krueger, L. Deobald, H. Gu,  A benchmark example for delamination growth predictions based on the single leg bending specimen under fatigue loading,  Advanced Modeling and Simulation in Engineering Sciences 7 (2020) 1–25.

##  Appendix A $\boldsymbol{\mathrm{H}}$ matrix

Before we calculate $U_{0}$ in equation 2.2.2, we need to calculate six parts.
The first part is:

$\displaystyle\frac{D_{11}}{2}\left(\frac{\partial^{2}w}{\partial
x^{2}}\right)^{2}=D_{11}(2\alpha_{1}^{2}+12x\alpha_{1}\alpha_{4}+4y\alpha_{1}\alpha_{5}+18x^{2}\alpha_{4}^{2}+12xy\alpha_{4}\alpha_{5}+2y^{2}\alpha_{5}^{2})$
(126)

Based on equation 16, we can obtain:

$\displaystyle\iint_{A}\frac{D_{11}}{2}\left(\frac{\partial^{2}w}{\partial
x^{2}}\right)^{2}\,\mathrm{d}x\mathrm{d}y=\frac{1}{2}\boldsymbol{\mathrm{\alpha}}^{\rm{T}}\,\boldsymbol{\mathrm{H}}_{1}\,\boldsymbol{\mathrm{\alpha}}$
(127)

where $\boldsymbol{\mathrm{H}}_{1}$ is:

$\displaystyle\boldsymbol{\mathrm{H}}_{1}=\iint_{A}D_{11}\begin{bmatrix}4&0&0&12x&4y&0&0\\\\[5.0pt]
&0&0&0&0&0&0\\\\[5.0pt] &&0&0&0&0&0\\\\[5.0pt] &&&36x^{2}&12xy&0&0\\\\[5.0pt]
&\rm{Symmetric}&&&4y^{2}&0&0\\\\[5.0pt] &&&&&0&0\\\\[5.0pt]
&&&&&&0\end{bmatrix}\mathrm{d}x\mathrm{d}y$ (128)

The second part is:

$\displaystyle\frac{D_{22}}{2}\left(\frac{\partial^{2}w}{\partial
y^{2}}\right)^{2}=D_{22}(2\alpha_{3}^{2}+4x\alpha_{3}\alpha_{6}+12y\alpha_{3}\alpha_{7}+2x^{2}\alpha_{6}^{2}+12xy\alpha_{6}\alpha_{7}+18y^{2}\alpha_{7}^{2})$
(129)

Then, an equation which is similar as equation 127 in terms of
$\boldsymbol{\mathrm{H}}_{2}$ can be obtained:

$\displaystyle\iint_{A}\frac{D_{22}}{2}\left(\frac{\partial^{2}w}{\partial
x^{2}}\right)^{2}\,\mathrm{d}x\mathrm{d}y=\frac{1}{2}\boldsymbol{\mathrm{\alpha}}^{\rm{T}}\,\boldsymbol{\mathrm{H}}_{2}\,\boldsymbol{\mathrm{\alpha}}$
(130)

where $\boldsymbol{\mathrm{H}}_{2}$ is:

$\displaystyle\boldsymbol{\mathrm{H}}_{2}=\iint_{A}D_{22}\begin{bmatrix}0&0&0&0&0&0&0\\\\[5.0pt]
&0&0&0&0&0&0\\\\[5.0pt] &&4&0&0&4x&12y\\\\[5.0pt] &&&0&0&0&0\\\\[5.0pt]
&\rm{Symmetric}&&&0&0&0\\\\[5.0pt] &&&&&4x^{2}&12xy\\\\[5.0pt]
&&&&&0&36y^{2}\end{bmatrix}\mathrm{d}x\mathrm{d}y$ (131)

The third part is:

$\displaystyle D_{12}\frac{\partial^{2}w}{\partial
x^{2}}\frac{\partial^{2}w}{\partial y^{2}}=$ $\displaystyle
D_{12}(4\alpha_{1}\alpha_{3}+4x\alpha_{1}\alpha_{6}+12x\alpha_{3}\alpha_{4}+12y\alpha_{1}\alpha_{7}+4y\alpha_{3}\alpha_{5}$
$\displaystyle+12x^{2}\alpha_{4}\alpha_{6}+12y^{2}\alpha_{5}\alpha_{7}+36xy\alpha_{4}\alpha_{7}+4xy\alpha_{5}\alpha_{6})$
(132)

and $\boldsymbol{\mathrm{H}}_{3}$ can be obtained by using the same method:

$\displaystyle\boldsymbol{\mathrm{H}}_{3}=\iint_{A}D_{12}\begin{bmatrix}0&0&4&0&0&4x&12y\\\\[5.0pt]
&0&0&0&0&0&0\\\\[5.0pt] &&0&12x&4y&0&0\\\\[5.0pt]
&&&0&0&12x^{2}&36xy\\\\[5.0pt] &\rm{Symmetric}&&&0&4xy&12y^{2}\\\\[5.0pt]
&&&&&0&0\\\\[5.0pt] &&&&&&0\end{bmatrix}\mathrm{d}x\mathrm{d}y$ (133)

Therefore, the $\boldsymbol{\mathrm{H}}$ matrix of $4^{\mathrm{th}}$,
$5^{\mathrm{th}}$ and $6^{\mathrm{th}}$ part are expressed by the following
expressions:

$\displaystyle\boldsymbol{\mathrm{H}}_{4}=\iint_{A}D_{16}\begin{bmatrix}0&4&0&0&8x&8y&0\\\\[5.0pt]
&0&0&12x&4y&0&0\\\\[5.0pt] &&0&0&0&0&0\\\\[5.0pt]
&&&0&24x^{2}&24xy&0\\\\[5.0pt] &\rm{Symmetric}&&&16xy&8y^{2}&0\\\\[5.0pt]
&&&&&0&0\\\\[5.0pt] &&&&&&0\end{bmatrix}\mathrm{d}x\mathrm{d}y$ (134)

$\displaystyle\boldsymbol{\mathrm{H}}_{5}=\iint_{A}D_{26}\begin{bmatrix}0&0&0&0&0&0&0\\\\[5.0pt]
&0&4&0&0&4x&12y\\\\[5.0pt] &&0&0&8x&8y&0\\\\[5.0pt] &&&0&0&0&0\\\\[5.0pt]
&\rm{Symmetric}&&&0&8x^{2}&24xy\\\\[5.0pt] &&&&&16xy&24y^{2}\\\\[5.0pt]
&&&&&&0\end{bmatrix}\mathrm{d}x\mathrm{d}y$ (135)

and

$\displaystyle\boldsymbol{\mathrm{H}}_{6}=\iint_{A}D_{66}\begin{bmatrix}0&0&0&0&0&0&0\\\\[5.0pt]
&4&0&0&8x&8y&0\\\\[5.0pt] &&0&0&0&0&0\\\\[5.0pt] &&&0&0&0&0\\\\[5.0pt]
&\rm{Symmetric}&&&16x^{2}&16xy&0\\\\[5.0pt] &&&&&16y^{2}&0\\\\[5.0pt]
&&&&&&0\end{bmatrix}\mathrm{d}x\mathrm{d}y$ (136)

The $\boldsymbol{\mathrm{H}}$ matrix can be written as the sum of
$\boldsymbol{\mathrm{H}}_{i}$ $(i=1,2,...,6)$:

$\displaystyle\boldsymbol{\mathrm{H}}=\boldsymbol{\mathrm{H}}_{1}+\boldsymbol{\mathrm{H}}_{2}+\boldsymbol{\mathrm{H}}_{3}+\boldsymbol{\mathrm{H}}_{4}+\boldsymbol{\mathrm{H}}_{5}+\boldsymbol{\mathrm{H}}_{6}$
(137)

where the integral for each item of the $\boldsymbol{\mathrm{H}}$ matrix can
be calculated by using the formulas given in Reference [29].

##  Appendix B $\boldsymbol{\mathrm{B}}$ matrix

Plugging the expression of $w$ (equation 13) into equation 2.2.2, we can
express the moment resultants as:

$\displaystyle M_{x}=$
$\displaystyle-2D_{11}\alpha_{1}-2D_{16}\alpha_{2}-2D_{12}\alpha_{3}-6xD_{11}\alpha_{4}-(2yD_{11}+4xD_{16})\alpha_{5}$
$\displaystyle-(2xD_{12}+4yD_{16})\alpha_{6}-6yD_{12}\alpha_{7}$ (138)

which in matrix form is:

$\displaystyle
M_{x}=\boldsymbol{\mathrm{B}}_{M_{x}}^{\rm{T}}\,\boldsymbol{\mathrm{\alpha}}$
(139)

where $\boldsymbol{\mathrm{B}}_{M_{x}}^{\rm{T}}$ is:

$\displaystyle\boldsymbol{\mathrm{B}}_{M_{x}}^{\rm{T}}=[$
$\displaystyle-2D_{11},\,-2D_{16},\,-2D_{12},\,-6xD_{11},\,-(2yD_{11}+4xD_{16}),$
$\displaystyle-(2xD_{12}+4yD_{16}),\,-6yD_{12}]$ (140)

Similarly, for $M_{y}$, we have:

$\displaystyle
M_{y}=\boldsymbol{\mathrm{B}}_{M_{y}}^{\rm{T}}\,\boldsymbol{\mathrm{\alpha}}$
(141)

where $\boldsymbol{\mathrm{B}}_{M_{y}}^{\rm{T}}$ is:

$\displaystyle\boldsymbol{\mathrm{B}}_{M_{y}}^{\rm{T}}=[$
$\displaystyle-2D_{12},\,-2D_{26},\,-2D_{22},\,-6xD_{12},\,-(2yD_{12}+4xD_{26}),$
$\displaystyle-(2xD_{22}+4yD_{26}),\,-6yD_{22}]$ (142)

For $M_{xy}$, we can obtain:

$\displaystyle
M_{xy}=\boldsymbol{\mathrm{B}}_{M_{xy}}^{\rm{T}}\,\boldsymbol{\mathrm{\alpha}}$
(143)

where the matrix $\boldsymbol{\mathrm{B}}_{M_{xy}}^{\rm{T}}$ is:

$\displaystyle\boldsymbol{\mathrm{B}}_{M_{xy}}^{\rm{T}}=[$
$\displaystyle-2D_{16},\,-2D_{66},\,-2D_{26},\,-6xD_{16},\,-(4xD_{66}+2yD_{16}),$
$\displaystyle-(4yD_{66}+2xD_{26}),\,-6yD_{26}]$ (144)

From equation 2.2.2, we can calculate the derivative of the moment:

$\displaystyle\frac{\partial M_{x}}{\partial x}$
$\displaystyle=-D_{11}\frac{\partial^{3}w}{\partial
x^{3}}-D_{12}\frac{\partial^{3}w}{\partial x\partial
y^{2}}-2D_{16}\frac{\partial^{3}w}{\partial x^{2}\partial y}$
$\displaystyle\frac{\partial M_{x}}{\partial y}$
$\displaystyle=-D_{11}\frac{\partial^{3}w}{\partial x^{2}\partial
y}-D_{12}\frac{\partial^{3}w}{\partial
y^{3}}-2D_{16}\frac{\partial^{3}w}{\partial x\partial y^{2}}$
$\displaystyle\frac{\partial M_{y}}{\partial x}$
$\displaystyle=-D_{12}\frac{\partial^{3}w}{\partial
x^{3}}-D_{22}\frac{\partial^{3}w}{\partial x\partial
y^{2}}-2D_{26}\frac{\partial^{3}w}{\partial x^{2}\partial y}$
$\displaystyle\frac{\partial M_{y}}{\partial y}$
$\displaystyle=-D_{12}\frac{\partial^{3}w}{\partial x^{2}\partial
y}-D_{22}\frac{\partial^{3}w}{\partial
y^{3}}-2D_{26}\frac{\partial^{3}w}{\partial x\partial y^{2}}$
$\displaystyle\frac{\partial M_{xy}}{\partial x}$
$\displaystyle=-D_{16}\frac{\partial^{3}w}{\partial
x^{3}}-D_{26}\frac{\partial^{3}w}{\partial x\partial
y^{2}}-2D_{66}\frac{\partial^{3}w}{\partial x^{2}\partial y}$
$\displaystyle\frac{\partial M_{xy}}{\partial y}$
$\displaystyle=-D_{16}\frac{\partial^{3}w}{\partial x^{2}\partial
y}-D_{26}\frac{\partial^{3}w}{\partial
y^{3}}-2D_{66}\frac{\partial^{3}w}{\partial x\partial y^{2}}$ (145)

Using equation 4, we calculate $M_{n}$ for node $1$ and node $2$ on the side
1-2 as:

$\displaystyle
M_{n}^{12}=\cos^{2}{\gamma_{12}}\boldsymbol{\mathrm{B}}_{M_{x}^{1}}^{\rm{T}}\boldsymbol{\mathrm{\alpha}}+\sin^{2}{\gamma_{12}}\boldsymbol{\mathrm{B}}_{M_{y}^{1}}^{\rm{T}}\boldsymbol{\mathrm{\alpha}}+\sin{2\gamma_{12}}\boldsymbol{\mathrm{B}}_{M_{xy}^{1}}^{\rm{T}}\boldsymbol{\mathrm{\alpha}}$
(146) $\displaystyle
M_{n}^{21}=\cos^{2}{\gamma_{12}}\boldsymbol{\mathrm{B}}_{M_{x}^{2}}^{\rm{T}}\boldsymbol{\mathrm{\alpha}}+\sin^{2}{\gamma_{12}}\boldsymbol{\mathrm{B}}_{M_{y}^{2}}^{\rm{T}}\boldsymbol{\mathrm{\alpha}}+\sin{2\gamma_{12}}\boldsymbol{\mathrm{B}}_{M_{xy}^{2}}^{\rm{T}}\boldsymbol{\mathrm{\alpha}}$
(147)

where $\boldsymbol{\mathrm{B}}_{M_{x}^{j}}^{\rm{T}}$ means equation B
evaluated at the coordinates of node $j$ and
$\boldsymbol{\mathrm{B}}_{M_{y}^{j}}^{\rm{T}}$ means equation B evaluated at
the coordinates of node $j$, respectively. And the corresponding
$\boldsymbol{\mathrm{B}}$ matrix is:

$\displaystyle\boldsymbol{\mathrm{B}}_{M_{n}^{12}}^{\rm{T}}=\cos^{2}{\gamma_{12}}\boldsymbol{\mathrm{B}}_{M_{x}^{1}}^{\rm{T}}+\sin^{2}{\gamma_{12}}\boldsymbol{\mathrm{B}}_{M_{y}^{1}}^{\rm{T}}+\sin{2\gamma_{12}}\boldsymbol{\mathrm{B}}_{M_{xy}^{1}}^{\rm{T}}$
(148)
$\displaystyle\boldsymbol{\mathrm{B}}_{M_{n}^{21}}^{\rm{T}}=\cos^{2}{\gamma_{12}}\boldsymbol{\mathrm{B}}_{M_{x}^{2}}^{\rm{T}}+\sin^{2}{\gamma_{12}}\boldsymbol{\mathrm{B}}_{M_{y}^{2}}^{\rm{T}}+\sin{2\gamma_{12}}\boldsymbol{\mathrm{B}}_{M_{xy}^{2}}^{\rm{T}}$
(149)

Similarly, $M_{n}$ for the side 2-3 is:

$\displaystyle
M_{n}^{23}=\cos^{2}{\gamma_{23}}\boldsymbol{\mathrm{B}}_{M_{x}^{2}}^{\rm{T}}\boldsymbol{\mathrm{\alpha}}+\sin^{2}{\gamma_{23}}\boldsymbol{\mathrm{B}}_{M_{y}^{2}}^{\rm{T}}\boldsymbol{\mathrm{\alpha}}+\sin{2\gamma_{23}}\boldsymbol{\mathrm{B}}_{M_{xy}^{2}}^{\rm{T}}\boldsymbol{\mathrm{\alpha}}$
(150) $\displaystyle
M_{n}^{32}=\cos^{2}{\gamma_{23}}\boldsymbol{\mathrm{B}}_{M_{x}^{3}}^{\rm{T}}\boldsymbol{\mathrm{\alpha}}+\sin^{2}{\gamma_{23}}\boldsymbol{\mathrm{B}}_{M_{y}^{3}}^{\rm{T}}\boldsymbol{\mathrm{\alpha}}+\sin{2\gamma_{23}}\boldsymbol{\mathrm{B}}_{M_{xy}^{3}}^{\rm{T}}\boldsymbol{\mathrm{\alpha}}$
(151)

And the corresponding $\boldsymbol{\mathrm{B}}$ matrix is:

$\displaystyle\boldsymbol{\mathrm{B}}_{M_{n}^{23}}^{\rm{T}}=\cos^{2}{\gamma_{23}}\boldsymbol{\mathrm{B}}_{M_{x}^{2}}^{\rm{T}}+\sin^{2}{\gamma_{23}}\boldsymbol{\mathrm{B}}_{M_{y}^{2}}^{\rm{T}}+\sin{2\gamma_{23}}\boldsymbol{\mathrm{B}}_{M_{xy}^{2}}^{\rm{T}}$
(152)
$\displaystyle\boldsymbol{\mathrm{B}}_{M_{n}^{32}}^{\rm{T}}=\cos^{2}{\gamma_{23}}\boldsymbol{\mathrm{B}}_{M_{x}^{3}}^{\rm{T}}+\sin^{2}{\gamma_{23}}\boldsymbol{\mathrm{B}}_{M_{y}^{3}}^{\rm{T}}+\sin{2\gamma_{23}}\boldsymbol{\mathrm{B}}_{M_{xy}^{3}}^{\rm{T}}$
(153)

Finally, $M_{n}$ for the side 3-1 is:

$\displaystyle
M_{n}^{31}=\cos^{2}{\gamma_{31}}\boldsymbol{\mathrm{B}}_{M_{x}^{3}}^{\rm{T}}\boldsymbol{\mathrm{\alpha}}+\sin^{2}{\gamma_{31}}\boldsymbol{\mathrm{B}}_{M_{y}^{3}}^{\rm{T}}\boldsymbol{\mathrm{\alpha}}+\sin{2\gamma_{31}}\boldsymbol{\mathrm{B}}_{M_{xy}^{3}}^{\rm{T}}\boldsymbol{\mathrm{\alpha}}$
(154) $\displaystyle
M_{n}^{13}=\cos^{2}{\gamma_{31}}\boldsymbol{\mathrm{B}}_{M_{x}^{1}}^{\rm{T}}\boldsymbol{\mathrm{\alpha}}+\sin^{2}{\gamma_{31}}\boldsymbol{\mathrm{B}}_{M_{y}^{1}}^{\rm{T}}\boldsymbol{\mathrm{\alpha}}+\sin{2\gamma_{31}}\boldsymbol{\mathrm{B}}_{M_{xy}^{1}}^{\rm{T}}\boldsymbol{\mathrm{\alpha}}$
(155)

And the corresponding $\boldsymbol{\mathrm{B}}$ matrix is:

$\displaystyle\boldsymbol{\mathrm{B}}_{M_{n}^{31}}^{\rm{T}}=\cos^{2}{\gamma_{31}}\boldsymbol{\mathrm{B}}_{M_{x}^{3}}^{\rm{T}}+\sin^{2}{\gamma_{31}}\boldsymbol{\mathrm{B}}_{M_{y}^{3}}^{\rm{T}}+\sin{2\gamma_{31}}\boldsymbol{\mathrm{B}}_{M_{xy}^{3}}^{\rm{T}}$
(156)
$\displaystyle\boldsymbol{\mathrm{B}}_{M_{n}^{13}}^{\rm{T}}=\cos^{2}{\gamma_{31}}\boldsymbol{\mathrm{B}}_{M_{x}^{1}}^{\rm{T}}+\sin^{2}{\gamma_{31}}\boldsymbol{\mathrm{B}}_{M_{y}^{1}}^{\rm{T}}+\sin{2\gamma_{31}}\boldsymbol{\mathrm{B}}_{M_{xy}^{1}}^{\rm{T}}$
(157)

Before calculating $R_{N}$, we need to calculate $M_{ns}$ at first as shown in
equation 2.2.2. For node 1:

$\displaystyle
M_{ns}^{12}=\left[\frac{1}{2}\left(\boldsymbol{\mathrm{B}}_{M_{y}^{1}}^{\rm{T}}-\boldsymbol{\mathrm{B}}_{M_{x}^{1}}^{\rm{T}}\right)\sin
2\gamma_{12}+\boldsymbol{\mathrm{B}}_{M_{xy}^{1}}^{\rm{T}}\cos
2\gamma_{12}\right]\,\boldsymbol{\mathrm{\alpha}}$ (158) $\displaystyle
M_{ns}^{13}=\left[\frac{1}{2}\left(\boldsymbol{\mathrm{B}}_{M_{y}^{1}}^{\rm{T}}-\boldsymbol{\mathrm{B}}_{M_{x}^{1}}^{\rm{T}}\right)\sin
2\gamma_{31}+\boldsymbol{\mathrm{B}}_{M_{xy}^{1}}^{\rm{T}}\cos
2\gamma_{31}\right]\,\boldsymbol{\mathrm{\alpha}}$ (159)

Thus, plugging the above two terms into equation 2.2.2, $R_{1}$ is:

$\displaystyle R_{1}$ $\displaystyle=\left[\frac{1}{2}(\sin 2\gamma_{12}-\sin
2\gamma_{31})\left(\boldsymbol{\mathrm{B}}_{M_{y}^{1}}^{\rm{T}}-\boldsymbol{\mathrm{B}}_{M_{x}^{1}}^{\rm{T}}\right)+(\cos
2\gamma_{12}-\cos
2\gamma_{31})\boldsymbol{\mathrm{B}}_{M_{xy}^{1}}^{\rm{T}}\right]\,\boldsymbol{\mathrm{\alpha}}$
$\displaystyle=\boldsymbol{\mathrm{B}}_{R_{1}}\,\boldsymbol{\mathrm{\alpha}}$
(160)

where the matrix $\boldsymbol{\mathrm{B}}_{R_{1}}$ is:

$\displaystyle\boldsymbol{\mathrm{B}}_{R_{1}}^{\rm{T}}=\frac{1}{2}\left(\sin
2\gamma_{12}-\sin
2\gamma_{31}\right)\left(\boldsymbol{\mathrm{B}}_{M_{y}^{1}}^{\rm{T}}-\boldsymbol{\mathrm{B}}_{M_{x}^{1}}^{\rm{T}}\right)+(\cos
2\gamma_{12}-\cos 2\gamma_{31})\boldsymbol{\mathrm{B}}_{M_{xy}^{1}}^{\rm{T}}$
(161)

Similarly, the matrix $\boldsymbol{\mathrm{B}}_{R_{2}}$ related to the node 2:

$\displaystyle\boldsymbol{\mathrm{B}}_{R_{2}}^{\rm{T}}=\frac{1}{2}\left(\sin
2\gamma_{23}-\sin
2\gamma_{12}\right)\left(\boldsymbol{\mathrm{B}}_{M_{y}^{2}}^{\rm{T}}-\boldsymbol{\mathrm{B}}_{M_{x}^{2}}^{\rm{T}}\right)+(\cos
2\gamma_{23}-\cos 2\gamma_{12})\boldsymbol{\mathrm{B}}_{M_{xy}^{2}}^{\rm{T}}$
(162)

At last, the matrix $\boldsymbol{\mathrm{B}}_{R_{3}}$:

$\displaystyle\boldsymbol{\mathrm{B}}_{R_{3}}^{\rm{T}}=\frac{1}{2}\left(\sin
2\gamma_{31}-\sin
2\gamma_{23}\right)\left(\boldsymbol{\mathrm{B}}_{M_{y}^{3}}^{\rm{T}}-\boldsymbol{\mathrm{B}}_{M_{x}^{3}}^{\rm{T}}\right)+(\cos
2\gamma_{31}-\cos 2\gamma_{23})\boldsymbol{\mathrm{B}}_{M_{xy}^{3}}^{\rm{T}}$
(163)

In order to simplify the writing of the above formula, we define the following
notations:

$\displaystyle\overline{S}_{1}=\sin{2\gamma_{12}}-\sin{2\gamma_{31}},\quad\overline{C}_{1}=\cos{2\gamma_{12}}-\cos{2\gamma_{31}}$
$\displaystyle\overline{S}_{2}=\sin{2\gamma_{23}}-\sin{2\gamma_{12}},\quad\overline{C}_{2}=\cos{2\gamma_{23}}-\cos{2\gamma_{12}}$
$\displaystyle\overline{S}_{3}=\sin{2\gamma_{31}}-\sin{2\gamma_{23}},\quad\overline{C}_{3}=\cos{2\gamma_{31}}-\cos{2\gamma_{23}}$
(164)

The matrix $\boldsymbol{\mathrm{B}}_{R}^{\rm{T}}$ in equation 161, 162 and 163
can be rewritten as:

$\displaystyle\boldsymbol{\mathrm{B}}_{R_{1}}^{\rm{T}}=\frac{1}{2}\,\overline{S}_{1}\left(\boldsymbol{\mathrm{B}}_{M_{y}^{1}}^{\rm{T}}-\boldsymbol{\mathrm{B}}_{M_{x}^{1}}^{\rm{T}}\right)+\overline{C}_{1}\boldsymbol{\mathrm{B}}_{M_{xy}^{1}}^{\rm{T}}$
$\displaystyle\boldsymbol{\mathrm{B}}_{R_{2}}^{\rm{T}}=\frac{1}{2}\,\overline{S}_{2}\left(\boldsymbol{\mathrm{B}}_{M_{y}^{2}}^{\rm{T}}-\boldsymbol{\mathrm{B}}_{M_{x}^{2}}^{\rm{T}}\right)+\overline{C}_{2}\boldsymbol{\mathrm{B}}_{M_{xy}^{2}}^{\rm{T}}$
$\displaystyle\boldsymbol{\mathrm{B}}_{R_{3}}^{\rm{T}}=\frac{1}{2}\,\overline{S}_{3}\left(\boldsymbol{\mathrm{B}}_{M_{y}^{3}}^{\rm{T}}-\boldsymbol{\mathrm{B}}_{M_{x}^{3}}^{\rm{T}}\right)+\overline{C}_{3}\boldsymbol{\mathrm{B}}_{M_{xy}^{3}}^{\rm{T}}$
(165)

Next, we calculate $V_{n}$ according to equations 5 to 8. Plugging equations 7
and 8 into 5, we get:

$\displaystyle V_{n}=\frac{\partial M_{n}}{\partial n}+2\frac{\partial
M_{ns}}{\partial s}$ (166)

Plugging equation 2.2.2 into equation 166, we get:

$\displaystyle V_{n}=\cos\gamma\frac{\partial M_{n}}{\partial
x}+\sin\gamma\frac{\partial M_{n}}{\partial y}-2\sin\gamma\frac{\partial
M_{ns}}{\partial x}+2\cos\gamma\frac{\partial M_{ns}}{\partial y}$ (167)

The terms on the RHS above require eventually the evaluation of
$\frac{\partial M_{n}}{\partial x}$, $\frac{\partial M_{n}}{\partial y}$,
$\frac{\partial M_{ns}}{\partial x}$, and $\frac{\partial M_{ns}}{\partial
y}$.

For $V_{n}^{12}$ item, the derivation from equation 146 is:

$\displaystyle\frac{\partial M_{n}^{12}}{\partial
x}=\cos^{2}{\gamma_{12}}\frac{\partial\boldsymbol{\mathrm{B}}_{M_{x}^{1}}^{\rm{T}}}{\partial
x}\boldsymbol{\mathrm{\alpha}}+\sin^{2}{\gamma_{12}}\frac{\partial\boldsymbol{\mathrm{B}}_{M_{y}^{1}}^{\rm{T}}}{\partial
x}\boldsymbol{\mathrm{\alpha}}+\sin{2\gamma_{12}}\frac{\partial\boldsymbol{\mathrm{B}}_{M_{xy}^{1}}^{\rm{T}}}{\partial
x}\boldsymbol{\mathrm{\alpha}}$ (168) $\displaystyle\frac{\partial
M_{n}^{12}}{\partial
y}=\cos^{2}{\gamma_{12}}\frac{\partial\boldsymbol{\mathrm{B}}_{M_{x}^{1}}^{\rm{T}}}{\partial
y}\boldsymbol{\mathrm{\alpha}}+\sin^{2}{\gamma_{12}}\frac{\partial\boldsymbol{\mathrm{B}}_{M_{y}^{1}}^{\rm{T}}}{\partial
y}\boldsymbol{\mathrm{\alpha}}+\sin{2\gamma_{12}}\frac{\partial\boldsymbol{\mathrm{B}}_{M_{xy}^{1}}^{\rm{T}}}{\partial
y}\boldsymbol{\mathrm{\alpha}}$ (169)

where

$\displaystyle\frac{\partial\boldsymbol{\mathrm{B}}_{M_{x}^{1}}^{\rm{T}}}{\partial
x}$ $\displaystyle=\begin{bmatrix}0,\ 0,\ 0,\ -6D_{11},\ -4D_{16},\ -2D_{12},\
0\end{bmatrix}$ (170)
$\displaystyle\frac{\partial\boldsymbol{\mathrm{B}}_{M_{y}^{1}}^{\rm{T}}}{\partial
x}$ $\displaystyle=\begin{bmatrix}0,\ 0,\ 0,\ -6D_{12},\ -4D_{26},\ -2D_{22},\
0\end{bmatrix}$ (171)
$\displaystyle\frac{\partial\boldsymbol{\mathrm{B}}_{M_{xy}^{1}}^{\rm{T}}}{\partial
x}$ $\displaystyle=\begin{bmatrix}0,\ 0,\ 0,\ -6D_{16},\ -4D_{66},\ -2D_{26},\
0\end{bmatrix}$ (172)

And

$\displaystyle\frac{\partial\boldsymbol{\mathrm{B}}_{M_{x}^{1}}^{\rm{T}}}{\partial
y}$ $\displaystyle=\begin{bmatrix}0,\ 0,\ 0,\ 0,\ -2D_{11},\ -4D_{16},\
-6D_{12}\end{bmatrix}$ (173)
$\displaystyle\frac{\partial\boldsymbol{\mathrm{B}}_{M_{y}^{1}}^{\rm{T}}}{\partial
y}$ $\displaystyle=\begin{bmatrix}0,\ 0,\ 0,\ 0,\ -2D_{12},\ -4D_{26},\
-6D_{22}\end{bmatrix}$ (174)
$\displaystyle\frac{\partial\boldsymbol{\mathrm{B}}_{M_{xy}^{1}}^{\rm{T}}}{\partial
y}$ $\displaystyle=\begin{bmatrix}0,\ 0,\ 0,\ 0,\ -2D_{16},\ -4D_{66},\
-6D_{26}\end{bmatrix}$ (175)

Thus, the matrix $\boldsymbol{\mathrm{B}}_{M_{n,x}^{12}}^{\rm{T}}$ and
$\boldsymbol{\mathrm{B}}_{M_{n,y}^{12}}^{\rm{T}}$ are:

$\displaystyle\boldsymbol{\mathrm{B}}_{M_{n,x}^{12}}^{\rm{T}}=\cos^{2}{\gamma_{12}}\frac{\partial\boldsymbol{\mathrm{B}}_{M_{x}^{1}}^{\rm{T}}}{\partial
x}+\sin^{2}{\gamma_{12}}\frac{\partial\boldsymbol{\mathrm{B}}_{M_{y}^{1}}^{\rm{T}}}{\partial
x}+\sin{2\gamma_{12}}\frac{\partial\boldsymbol{\mathrm{B}}_{M_{xy}^{1}}^{\rm{T}}}{\partial
x}$ (176)
$\displaystyle\boldsymbol{\mathrm{B}}_{M_{n,y}^{12}}^{\rm{T}}=\cos^{2}{\gamma_{12}}\frac{\partial\boldsymbol{\mathrm{B}}_{M_{x}^{1}}^{\rm{T}}}{\partial
y}+\sin^{2}{\gamma_{12}}\frac{\partial\boldsymbol{\mathrm{B}}_{M_{y}^{1}}^{\rm{T}}}{\partial
y}+\sin{2\gamma_{12}}\frac{\partial\boldsymbol{\mathrm{B}}_{M_{xy}^{1}}^{\rm{T}}}{\partial
y}$ (177)

Using equation 7, the derivative of moment $M_{ns}$ on side 1-2 is:

$\displaystyle\frac{\partial M_{ns}^{12}}{\partial
x}=\left[\frac{1}{2}\left(\frac{\partial\boldsymbol{\mathrm{B}}_{M_{y}^{1}}^{\rm{T}}}{\partial
x}-\frac{\partial\boldsymbol{\mathrm{B}}_{M_{x}^{1}}^{\rm{T}}}{\partial
x}\right)\sin
2\gamma_{12}+\frac{\partial\boldsymbol{\mathrm{B}}_{M_{xy}^{1}}^{\rm{T}}}{\partial
x}\cos 2\gamma_{12}\right]\,\boldsymbol{\mathrm{\alpha}}$ (178)
$\displaystyle\frac{\partial M_{ns}^{12}}{\partial
y}=\left[\frac{1}{2}\left(\frac{\partial\boldsymbol{\mathrm{B}}_{M_{y}^{1}}^{\rm{T}}}{\partial
y}-\frac{\partial\boldsymbol{\mathrm{B}}_{M_{x}^{1}}^{\rm{T}}}{\partial
y}\right)\sin
2\gamma_{12}+\frac{\partial\boldsymbol{\mathrm{B}}_{M_{xy}^{1}}^{\rm{T}}}{\partial
y}\cos 2\gamma_{12}\right]\,\boldsymbol{\mathrm{\alpha}}$ (179)

Thus, the matrices $\boldsymbol{\mathrm{B}}_{M_{ns,x}^{12}}^{\rm{T}}$ and
$\boldsymbol{\mathrm{B}}_{M_{ns,y}^{12}}^{\rm{T}}$ are:

$\displaystyle\boldsymbol{\mathrm{B}}_{M_{ns,x}^{12}}^{\rm{T}}=\frac{1}{2}\sin
2\gamma_{12}\left(\frac{\partial\boldsymbol{\mathrm{B}}_{M_{y}^{1}}^{\rm{T}}}{\partial
x}-\frac{\partial\boldsymbol{\mathrm{B}}_{M_{x}^{1}}^{\rm{T}}}{\partial
x}\right)+\cos
2\gamma_{12}\frac{\partial\boldsymbol{\mathrm{B}}_{M_{xy}^{1}}^{\rm{T}}}{\partial
x}$ (180)
$\displaystyle\boldsymbol{\mathrm{B}}_{M_{ns,y}^{12}}^{\rm{T}}=\frac{1}{2}\sin
2\gamma_{12}\left(\frac{\partial\boldsymbol{\mathrm{B}}_{M_{y}^{1}}^{\rm{T}}}{\partial
y}-\frac{\partial\boldsymbol{\mathrm{B}}_{M_{x}^{1}}^{\rm{T}}}{\partial
y}\right)+\cos
2\gamma_{12}\frac{\partial\boldsymbol{\mathrm{B}}_{M_{xy}^{1}}^{\rm{T}}}{\partial
y}$ (181)

Finally, using equation 167, $\boldsymbol{\mathrm{B}}$ matrix for the
Kirchhoff shear force $V_{n}$ on side 1-2 is:

$\displaystyle\boldsymbol{\mathrm{B}}_{V_{n}^{12}}^{\rm{T}}=$
$\displaystyle\cos{\gamma_{12}}\boldsymbol{\mathrm{B}}_{M_{n,x}^{12}}^{\rm{T}}+\sin{\gamma_{12}}\boldsymbol{\mathrm{B}}_{M_{n,y}^{12}}^{\rm{T}}-2\sin{\gamma_{12}}\boldsymbol{\mathrm{B}}_{M_{ns,x}^{12}}^{\rm{T}}+2\cos{\gamma_{12}}\boldsymbol{\mathrm{B}}_{M_{ns,y}^{12}}^{\rm{T}}$
(182)

Similarly, $\boldsymbol{\mathrm{B}}_{V_{n}}$ matrix on side 2-3 is:

$\displaystyle\boldsymbol{\mathrm{B}}_{M_{n,x}^{23}}^{\rm{T}}=\cos^{2}{\gamma_{23}}\frac{\partial\boldsymbol{\mathrm{B}}_{M_{x}^{2}}^{\rm{T}}}{\partial
x}+\sin^{2}{\gamma_{23}}\frac{\partial\boldsymbol{\mathrm{B}}_{M_{y}^{2}}^{\rm{T}}}{\partial
x}+\sin{2\gamma_{23}}\frac{\partial\boldsymbol{\mathrm{B}}_{M_{xy}^{1}}^{\rm{T}}}{\partial
x}$ (183)
$\displaystyle\boldsymbol{\mathrm{B}}_{M_{n,y}^{23}}^{\rm{T}}=\cos^{2}{\gamma_{23}}\frac{\partial\boldsymbol{\mathrm{B}}_{M_{x}^{2}}^{\rm{T}}}{\partial
y}+\sin^{2}{\gamma_{23}}\frac{\partial\boldsymbol{\mathrm{B}}_{M_{y}^{2}}^{\rm{T}}}{\partial
y}+\sin{2\gamma_{23}}\frac{\partial\boldsymbol{\mathrm{B}}_{M_{xy}^{2}}^{\rm{T}}}{\partial
y}$ (184)

and

$\displaystyle\boldsymbol{\mathrm{B}}_{M_{ns,x}^{23}}^{\rm{T}}=\frac{1}{2}\sin
2\gamma_{23}\left(\frac{\partial\boldsymbol{\mathrm{B}}_{M_{y}^{2}}^{\rm{T}}}{\partial
x}-\frac{\partial\boldsymbol{\mathrm{B}}_{M_{x}^{2}}^{\rm{T}}}{\partial
x}\right)+\cos
2\gamma_{23}\frac{\partial\boldsymbol{\mathrm{B}}_{M_{xy}^{2}}^{\rm{T}}}{\partial
x}$ (185)
$\displaystyle\boldsymbol{\mathrm{B}}_{M_{ns,y}^{23}}^{\rm{T}}=\frac{1}{2}\sin
2\gamma_{23}\left(\frac{\partial\boldsymbol{\mathrm{B}}_{M_{y}^{2}}^{\rm{T}}}{\partial
y}-\frac{\partial\boldsymbol{\mathrm{B}}_{M_{x}^{2}}^{\rm{T}}}{\partial
y}\right)+\cos
2\gamma_{23}\frac{\partial\boldsymbol{\mathrm{B}}_{M_{xy}^{2}}^{\rm{T}}}{\partial
y}$ (186)

The $\boldsymbol{\mathrm{B}}$ matrix for $V_{n}^{23}$ is:

$\displaystyle\boldsymbol{\mathrm{B}}_{V_{n}^{23}}^{\rm{T}}=$
$\displaystyle\cos{\gamma_{23}}\boldsymbol{\mathrm{B}}_{M_{n,x}^{23}}^{\rm{T}}+\sin{\gamma_{23}}\boldsymbol{\mathrm{B}}_{M_{n,y}^{23}}^{\rm{T}}-2\sin{\gamma_{23}}\boldsymbol{\mathrm{B}}_{M_{ns,x}^{23}}^{\rm{T}}+2\cos{\gamma_{23}}\boldsymbol{\mathrm{B}}_{M_{ns,y}^{23}}^{\rm{T}}$
(187)

On side 3-1:

$\displaystyle\boldsymbol{\mathrm{B}}_{M_{n,x}^{31}}^{\rm{T}}=\cos^{2}{\gamma_{31}}\frac{\partial\boldsymbol{\mathrm{B}}_{M_{x}^{3}}^{\rm{T}}}{\partial
x}+\sin^{2}{\gamma_{31}}\frac{\partial\boldsymbol{\mathrm{B}}_{M_{y}^{3}}^{\rm{T}}}{\partial
x}+\sin{2\gamma_{31}}\frac{\partial\boldsymbol{\mathrm{B}}_{M_{xy}^{3}}^{\rm{T}}}{\partial
x}$ (188)
$\displaystyle\boldsymbol{\mathrm{B}}_{M_{n,y}^{31}}^{\rm{T}}=\cos^{2}{\gamma_{31}}\frac{\partial\boldsymbol{\mathrm{B}}_{M_{x}^{3}}^{\rm{T}}}{\partial
y}+\sin^{2}{\gamma_{31}}\frac{\partial\boldsymbol{\mathrm{B}}_{M_{y}^{3}}^{\rm{T}}}{\partial
y}+\sin{2\gamma_{31}}\frac{\partial\boldsymbol{\mathrm{B}}_{M_{xy}^{3}}^{\rm{T}}}{\partial
y}$ (189)

and

$\displaystyle\boldsymbol{\mathrm{B}}_{M_{ns,x}^{31}}^{\rm{T}}=\frac{1}{2}\sin
2\gamma_{31}\left(\frac{\partial\boldsymbol{\mathrm{B}}_{M_{y}^{3}}^{\rm{T}}}{\partial
x}-\frac{\partial\boldsymbol{\mathrm{B}}_{M_{x}^{3}}^{\rm{T}}}{\partial
x}\right)+\cos
2\gamma_{31}\frac{\partial\boldsymbol{\mathrm{B}}_{M_{xy}^{3}}^{\rm{T}}}{\partial
x}$ (190)
$\displaystyle\boldsymbol{\mathrm{B}}_{M_{ns,y}^{31}}^{\rm{T}}=\frac{1}{2}\sin
2\gamma_{31}\left(\frac{\partial\boldsymbol{\mathrm{B}}_{M_{y}^{3}}^{\rm{T}}}{\partial
y}-\frac{\partial\boldsymbol{\mathrm{B}}_{M_{x}^{3}}^{\rm{T}}}{\partial
y}\right)+\cos
2\gamma_{31}\frac{\partial\boldsymbol{\mathrm{B}}_{M_{xy}^{3}}^{\rm{T}}}{\partial
y}$ (191)

The $\boldsymbol{\mathrm{B}}$ matrix for $V_{n}^{31}$ is:

$\displaystyle\boldsymbol{\mathrm{B}}_{V_{n}^{31}}^{\rm{T}}=$
$\displaystyle\cos{\gamma_{31}}\boldsymbol{\mathrm{B}}_{M_{n,x}^{31}}^{\rm{T}}+\sin{\gamma_{23}}\boldsymbol{\mathrm{B}}_{M_{n,y}^{31}}^{\rm{T}}-2\sin{\gamma_{23}}\boldsymbol{\mathrm{B}}_{M_{ns,x}^{31}}^{\rm{T}}+2\cos{\gamma_{23}}\boldsymbol{\mathrm{B}}_{M_{ns,y}^{31}}^{\rm{T}}$

Therefore, the calculations of $\boldsymbol{\mathrm{B}}$ matrices for $V_{n}$
on sides 23 and 31 can be simplified through the use of the pre-calculated
matrices for side 12.

Now, all the components of $\boldsymbol{\mathrm{B}}$ matrix have been derived.
The assembled (7 x 12) matrix $\boldsymbol{\mathrm{B}}$ is:

$\displaystyle\boldsymbol{\mathrm{B}}=[$
$\displaystyle\boldsymbol{\mathrm{B}}_{R_{1}},\boldsymbol{\mathrm{B}}_{R_{2}},\boldsymbol{\mathrm{B}}_{R_{3}},\boldsymbol{\mathrm{B}}_{V_{n}^{12}},\boldsymbol{\mathrm{B}}_{V_{n}^{23}},\boldsymbol{\mathrm{B}}_{V_{n}^{31}},\boldsymbol{\mathrm{B}}_{M_{n}^{12}},\boldsymbol{\mathrm{B}}_{M_{n}^{21}},\boldsymbol{\mathrm{B}}_{M_{n}^{23}},\boldsymbol{\mathrm{B}}_{M_{n}^{32}},\boldsymbol{\mathrm{B}}_{M_{n}^{31}},\boldsymbol{\mathrm{B}}_{M_{n}^{13}}]$
(192)

##  Appendix C Analytical solution of SLB benchmark

Figure 1: Analytical solution of SLB model

Referring to Figure 1, the analytical load-displacement curve consists of
three parts. The first part, OB, is the linear elastic stage with the initial
crack length $a_{0}$. The second part comes from the curve ABC which
represents the load-displacement response during delamination propagation when
$a<L$ (the half length of the SLB model). When $a>L$, we can obtain the last
curve DFE. The full analytical curve would then be OBFE. The derivation
process of each curve is shown one by one in the following content.

Figure 2: Beam model of SLB specimen ($a<L$)

When the crack length $a$ is smaller than $L$, the analytical derivation of
the SLB benchmark follows from the work of Davidson et el. [35]. Referring to
the SLB model in Figure 8, it is evident that the bending stiffness in the
cracked region comes entirely from the part above the crack, which is midplane
symmetrical. The uncracked region is midplane anti-symmetric.

The SLB model ($a<L$) is simplified into a 2D beam model, as shown in Figure
2. The moment-curvature relationship for the SLB problem can be expressed as:

$\displaystyle
M=-b\,D\,\frac{\mathrm{d}^{2}w}{\mathrm{d}x^{2}}=b\,D\,\kappa_{x}$ (193)

where $M$ is the moment on a face along the x-axis, $b$ is the width of the
specimen, $w$ is the deflection along the $z$ direction. And $D$ is the
effective bending rigidity per unit width. Because the values of $D$ in the
uncracked region and cracked regions are different, $D_{0}$ represents the
effective bending rigidity of the uncracked region, and $D_{1}$ represents the
effective bending rigidity of the cracked region to distinguish them.

We assume that the classical laminate theory is applicable to the SLB problem.
In matrix form, the constitutive equation can be written as:

$\displaystyle\begin{bmatrix}N_{x}\\\\[5.0pt] N_{y}\\\\[5.0pt]
N_{xy}\\\\[5.0pt] \--\\\\[5.0pt] M_{x}\\\\[5.0pt] M_{y}\\\\[5.0pt]
M_{xy}\end{bmatrix}=\begin{bmatrix}A_{11}&A_{12}&A_{16}&|&B_{11}&B_{12}&B_{16}\\\\[5.0pt]
A_{12}&A_{22}&A_{26}&|&B_{12}&B_{22}&B_{26}\\\\[5.0pt]
A_{16}&A_{26}&A_{66}&|&B_{16}&B_{26}&B_{66}\\\\[5.0pt]
\--&\--&\--&|&\--&\--&\--\\\\[5.0pt]
B_{11}&B_{12}&B_{16}&|&D_{11}&D_{12}&D_{16}\\\\[5.0pt]
B_{12}&B_{22}&B_{26}&|&D_{12}&D_{22}&D_{26}\\\\[5.0pt]
B_{16}&B_{26}&B_{66}&|&D_{16}&D_{26}&D_{66}\end{bmatrix}\begin{bmatrix}\varepsilon_{x}^{0}\\\\[5.0pt]
\varepsilon_{y}^{0}\\\\[5.0pt] \gamma_{xy}^{0}\\\\[5.0pt] \--\\\\[5.0pt]
\kappa_{x}\\\\[5.0pt] \kappa_{y}\\\\[5.0pt] \kappa_{xy}\end{bmatrix}$ (194)

where $\varepsilon_{x}^{0},\,\varepsilon_{y}^{0},\,\gamma_{xy}^{0}$ are the
membrane strains of the neutral plane, and
$\kappa_{x},\,\kappa_{y},\,\kappa_{xy}$ are the curvatures, derived from the
out-of-plane displacement $w$:

$\displaystyle\kappa_{x}=-\frac{\partial^{2}w}{\partial
x^{2}},\quad\kappa_{y}=-\frac{\partial^{2}w}{\partial
y^{2}},\quad\kappa_{xy}=-\frac{\partial^{2}w}{\partial x\partial y}$ (195)

The inverted form of the constitutive relationship is:

$\displaystyle\begin{bmatrix}\varepsilon_{x}^{0}\\\\[5.0pt]
\varepsilon_{y}^{0}\\\\[5.0pt] \gamma_{xy}^{0}\\\\[5.0pt] \--\\\\[5.0pt]
\kappa_{x}\\\\[5.0pt] \kappa_{y}\\\\[5.0pt]
\kappa_{xy}\end{bmatrix}=\begin{bmatrix}\alpha_{11}&\alpha_{12}&\alpha_{16}&|&\beta_{11}&\beta_{12}&\beta_{16}\\\\[5.0pt]
\alpha_{12}&\alpha_{22}&\alpha_{26}&|&\beta_{12}&\beta_{22}&\beta_{26}\\\\[5.0pt]
\alpha_{16}&\alpha_{26}&\alpha_{66}&|&\beta_{16}&\beta_{26}&\beta_{66}\\\\[5.0pt]
\--&\--&\--&|&\--&\--&\--\\\\[5.0pt]
\beta_{11}&\beta_{12}&\beta_{16}&|&\delta_{11}&\delta_{12}&\delta_{16}\\\\[5.0pt]
\beta_{12}&\beta_{22}&\beta_{26}&|&\delta_{12}&\delta_{22}&\delta_{26}\\\\[5.0pt]
\beta_{16}&\beta_{26}&\beta_{66}&|&\delta_{16}&\delta_{26}&\delta_{66}\end{bmatrix}\begin{bmatrix}N_{x}\\\\[5.0pt]
N_{y}\\\\[5.0pt] N_{xy}\\\\[5.0pt] \--\\\\[5.0pt] M_{x}\\\\[5.0pt]
M_{y}\\\\[5.0pt] M_{xy}\end{bmatrix}$ (196)

If the plate constraint condition is “generalized plane stress”, $N_{y}$,
$N_{xy}$, $M_{y}$, and $M_{xy}$ are zero. When considering $N_{x}=0$, equation
196 gives:

$\displaystyle D=\frac{1}{\delta_{11}}$ (197)

For the plane strain condition
($\epsilon_{y}^{0}=\gamma_{xy}^{0}=\kappa_{y}=\kappa_{xy}=0$), equation 194
gives:

$\displaystyle D=D_{11}$ (198)

The expression of strain energy is:

$\displaystyle U=\int_{0}^{2L}\frac{M^{2}\mathrm{d}x}{2bD}$ (199)

From Figure 2, the boundary conditions of the SLB model can be considered
similar to three-point bending. So we can get the moment along the beam as:

$\displaystyle M=\begin{cases}-\dfrac{Px}{2},&0\leq x\leq L\\\\[10.00002pt]
\dfrac{Px}{2}-PL,&L\leq x\leq 2L\\\ \end{cases}$ (200)

Since the cross-section of the upper part of the cracked region $(0\leq x\leq
a)$ is different from that of the uncracked region, the moment $M$ is further
divided, and the expression for strain energy is the sum of three parts:

$\displaystyle\left.U\right|_{(a<L)}=\int_{0}^{a}\frac{M^{2}\mathrm{d}x}{2bD_{1}}+\int_{a}^{L-a}\frac{M^{2}\mathrm{d}x}{2bD_{0}}+\int_{L}^{2L}\frac{M^{2}\mathrm{d}x}{2bD_{0}}$
(201)

Substituting equation 200 into equation 201, we can get:

$\displaystyle\left.U\right|_{(a<L)}=\int_{0}^{a}\frac{\left({-\frac{Px}{2}}\right)^{2}\mathrm{d}x}{2bD_{1}}+\int_{a}^{L}\frac{\left({-\frac{Px}{2}}\right)^{2}\mathrm{d}x}{2bD_{0}}+\int_{L}^{2L}\frac{\left({\frac{Px}{2}-PL}\right)^{2}\mathrm{d}x}{2bD_{0}}$
(202)

Based on the expression of strain energy, the displacement $\delta$ can be
calculated by Castigliano’s second theorem:

$\displaystyle\left.\delta\right|_{(a<L)}$
$\displaystyle=\frac{\partial{\left.U\right|_{(a<L)}}}{\partial{P}}=\int_{0}^{a}\frac{Px^{2}\mathrm{d}x}{4bD_{1}}+\int_{a}^{L}\frac{Px^{2}\mathrm{d}x}{4bD_{0}}+\int_{L}^{2L}\frac{\left({Px^{2}-4PLx+4PL^{2}}\right)\mathrm{d}x}{4bD_{0}}$
$\displaystyle=\frac{Pa^{3}}{12bD_{1}}+\frac{2PL^{3}-Pa^{3}}{12bD_{0}}$ (203)

The above equation can be used to plot the $P-\delta$ curve for a certain
crack length $a$, which for the case of $a=a_{0}$ gives the first part of the
analytical curve, i.e., OB. The compliance $C$ during the elastic region is
defined as the center deflection $\delta$ divided by the center load $P$:

$\displaystyle\left.C\right|_{(a<L)}$
$\displaystyle=\frac{\left.\delta\right|_{(a<L)}}{P}=\frac{2L^{3}+a^{3}(R-1)}{12bD_{0}}$
(204)

where $R$ is the ratio of the bending rigidity of the uncracked region to that
of the cracked region:

$\displaystyle R=D_{0}/D_{1}$ (205)

The relationship between the critical energy release rate $G_{c}$ and the
derivative of compliance $C$ with respect to the crack length $a$ is:

$\displaystyle G_{c}=\frac{P^{2}}{2b}\,\frac{\partial C}{\partial a}$ (206)

Substituting equation 204 into equation 206, we can obtain:

$\displaystyle G_{c}=\frac{P^{2}a^{2}(R-1)}{8b^{2}D_{0}}$ (207)

Substituting equation 204 to equation 207, the critical energy release rate
$G_{c}$ can also be obtained by the critical load and deflection:

$\displaystyle
G_{c}=\frac{3Pa^{2}\left.\delta\right|_{(a<L)}}{2b}\,\frac{(R-1)}{[2L^{3}+a^{3}(R-1)]}$
(208)

The critical energy release rate $G_{c}$ in equation 208 is the value in mixed
mode. It can be calculated by:

$\displaystyle
G_{\mathrm{c}}=G_{\mathrm{Ic}}+(G_{\mathrm{IIc}}-G_{\mathrm{Ic}})B^{\eta}$
(209)

where $G_{\mathrm{Ic}}$ and $G_{\mathrm{IIc}}$ are the critical energy release
rates in pure mode I and mode II, respectively. The mixed mode ratio $B$ in
the SLB model is 0.4 [36]. Equation 207 can be used to calculate the critical
load $P$ at a certain crack length $a$. It can also be used to express $a$ in
terms of $P$ for a fixed $G_{c}$, which can be plugged into equation 208 to
obtain the relationship between $\delta$ and $P$ for a fixed $G_{c}$ for the
case of $a<L$, hence giving the curve ABC.

Figure 3: Beam model of SLB specimen ($a\geq L$)

When the crack length $a$ is greater than or equal to $L$, the analytical
derivation of the SLB benchmark can follow the same derivation process as
above. The SLB model ($a\geq L$) is simplified in Figure 3. Based on Figure 3,
the expression for strain energy is rewritten as:

$\displaystyle\left.U\right|_{(a\geq
L)}=\int_{0}^{L}\frac{M^{2}\mathrm{d}x}{2bD_{1}}+\int_{L}^{a}\frac{M^{2}\mathrm{d}x}{2bD_{1}}+\int_{a}^{2L}\frac{M^{2}\mathrm{d}x}{2bD_{0}}$
(210)

Since the boundary conditions have not changed, the expression of the moment
is consistent with equation 200. Substituting equation 200 into equation 210,
we can obtain:

$\displaystyle\left.U\right|_{(a\geq L)}=$
$\displaystyle\int_{0}^{L}\frac{P^{2}x^{2}\mathrm{d}x}{8bD_{1}}+\int_{L}^{a}\frac{\left(P^{2}x^{2}-4P^{2}Lx+4P^{2}L^{2}\right)\mathrm{d}x}{8bD_{1}}$
$\displaystyle+\int_{a}^{2L}\frac{\left(P^{2}x^{2}-4P^{2}Lx+4P^{2}L^{2}\right)\mathrm{d}x}{8bD_{0}}$
(211)

The displacement $\delta$ at the middle can also be found by Castigliano’s
second theorem:

$\displaystyle\left.\delta\right|_{(a\geq
L)}=-\frac{PL^{3}}{2bD_{1}}+\frac{Pa^{3}-6PLa^{2}+12PL^{2}a}{12bD_{1}}+\frac{8PL^{3}}{12bD_{0}}-\frac{Pa^{3}-6PLa^{2}+12PL^{2}a}{12bD_{0}}$
(212)

Introducing ratio $R$ in equation 205, the above equation can be rewritten as:

$\displaystyle\left.\delta\right|_{(a\geq
L)}=-\frac{(6R-8)PL^{3}}{12bD_{0}}+\frac{(R-1)(Pa^{3}-6PLa^{2}+12PL^{2}a)}{12bD_{0}}$
(213)

The expression of compliance $C$ under the condition $a\geq L$ is:

$\displaystyle\left.C\right|_{(a\geq L)}$
$\displaystyle=\frac{\left.\delta\right|_{(a\geq
L)}}{P}=-\frac{(6R-8)L^{3}}{12bD_{0}}+\frac{(R-1)(a^{3}-6La^{2}+12L^{2}a)}{12bD_{0}}$
(214)

Substituting equation 214 into equation 206, we have:

$\displaystyle G_{c}=\frac{P^{2}(R-1)(3a^{2}-12aL+12L^{2})}{24b^{2}D_{0}}$
(215)

where $G_{c}$ is the critical energy release rate given by equation 209. Using
equations 213 and 215, one can then vary $a$ to plot the third part of the
analytical curve DE.