1 Introduction↩︎

Fiber reinforced polymer composites are increasingly adopted in the aerospace industry to meet lightweighting requirements. In particular, laminated stiffened panels made of thermoplastic composites can be fusion bonded, without the need of additional materials such as adhesives or bolts, resulting in more weight-savings, faster processing cycles and the possibility to manufacture composite parts with more complex geometries. The overmolding of a dissimilar material in a thermoplastic laminate is also an alternative that allows for high performance stiffened structural parts [1]. However, the failure behaviour of these fusion bonded/overmolded thermoplastic composites strongly depends on the processing conditions [1]–[4]. At present, the lack of efficient and reliable performance prediction tools forms an obstacle to the widespread use of thermoplastic composites.

Describing the mechanical behavior of composites through finite element simulations has been traditionally a challenge because of the complex anisotropic, inhomogeneous and multiscale nature of these materials. The most common approach in finite element modelling of composite materials is based on macroscale modelling, where layered shells, stacked solid or stacked shell elements are used to represent the laminate. However, such modelling strategies alone cannot explicitly capture interlaminar damage. To simulate delamination between laminate plies, cohesive elements are commonly introduced between adjacent layers [5], [6]. The cohesive elements are developed based on the cohesive zone model, proposed by Dugdale and Barenblatt [7], [8], in which a fracture process zone exists along the interface, ahead of the stress-free crack tip. A traction-separation relationship describes how the interfacial stresses and damage evolve with respect to the interfacial openings. Standard cohesive elements are usually developed for use between two solid elements to model their debonding [9]–[11].

Another common failure mode of laminated composite panels is skin-stiffener debonding due to relatively low interface strengths [12]. In thermoplastic composites, even if good welding occurs, the weaker zone moves to the laminate side at the interface between the matrix and the reinforcement, resulting in a weaker zone at the welding area [13]. Cohesive element models are also a versatile approach for the modeling of skin-stiffener separation in composite stiffened panels [12], [14], [15]. For instance, Balzani and Wagner [14] examined debonding between skin and stiffener with cohesive elements, implemented as surface interface elements between shells. Akterskaia et al. [12] developed a global-local methodology to evaluate skin–stiffener debonding with two accuracy levels: a fast global analysis to identify critical areas using a coarse mesh, followed by refined local submodels where cohesive elements simulate detailed damage propagation. The fracture strength of thermoplastic composite T-joints was investigated through numerical simulations in [15], [16]. While [15] evaluated mode I fracture strength using a mode-I cohesive model, [16] employed a framework combining a mixed-mode cohesive model with an anisotropic viscoplasticity model for the laminate.

Although widely used for delamination and skin–stiffener debonding modeling, standard cohesive elements require very fine meshes because the element size must be much smaller than the cohesive-zone length. High stress gradients develop within the cohesive zone during failure in composites [9], [11]. Therefore, a sufficiently refined mesh is required to accurately capture the solution. The fine mesh requirement of cohesive-element modeling has drawn the attention of many researchers in the past [17]–[23]. The so-called structural cohesive elements, which conform to \(C^1\) continuous structural elements, i.e., beams, plates and shells, are a recent approach that has demonstrated potential for delamination analysis with coarse meshes [22]–[24]. Shell models based on the Kirchhoff hypothesis have previously been adopted in the formulation of the structural cohesive elements [22], [23], [25]. The advantages of such models stem from the fact that the higher-order smooth representation of the shell mid-surface displacement fields allows to adopt relatively coarse discretizations without sacrificing solution accuracy. In multi-layer delamination modelling, compliant cohesive interfaces can still provide the transverse shear deformation expected in a laminate that isolated Kirchhoff models cannot capture. For instance, Balducci and Chen [22] developed a structural cohesive element for delamination problems which is compatible with TUBA3 plate elements [26]. Their results showed that the formulation allowed for the use of larger elements in delamination analysis. However, the curvature degrees of freedom (DoFs) of the element make it more complicated to set boundary conditions. Recently, Ai et al. [23] extended the Allman [27] triangular Kirchhoff plate element, which has only DoFs commonly used by engineers, for the modeling of composite plies. They also developed the corresponding structural cohesive element designed for delamination modeling. Even though the hybrid shell formulation employed in [23] does not actually give a smooth representation of the shell mid-surface transversal displacement in a strong sense, it does enforce higher-order continuity in a weak sense, which appears to be enough to give more accurate results than the standard models [23].

This article introduces a new structural cohesive element formulation specifically designed for the skin-stiffener debonding problem. First, the plate element from [23] is extended to 3D shell models by using the triangular membrane element with drilling DoFs developed by Bergan and Felippa [28]. This ensures higher-order approximations for the in-plane displacements and avoids membrane locking issues in case of in-plane bending of the stiffeners, while preserving a simple shell formulation based on a 3-node triangular element with DoFs commonly used to impose boundary conditions in structural analysis [29], [30]. Then, the novel structural cohesive element, namely the shell-to-shell cohesive line element, is proposed for the efficient modeling of skin-stiffener debonding in overmolded panels. The main novelty is the proposed kinematics that allows to obtain the jump vector over the skin-stiffener interface from the higher-order shell mid-surface displacement fields evaluated at the element edges. The model is suitable for analyzing debonding in skin-stiffener structures with non-constant damage through the stiffener thickness. It is designed to capture the overall skin–stiffener debonding response when the panel failure is governed solely by interface damage. Further developments are required to account for intralaminar damage or plasticity.

The paper is outlined as follows. The governing equations and finite element approximations for the shell element are briefly reviewed in Section 2. The proposed structural cohesive element is presented in Section 3. Numerical results are presented in Section 4. The model is first verified using mode I, mode II, and mixed-mode classical benchmarks. The performance of the proposed cohesive element is then demonstrated through comparisons with standard 3D cohesive elements in debonding problems. The debonding analysis of a complex stiffened panel is also presented for illustrative purposes. In Section 5, we draw the conclusions of this work.

2 Hybrid/mixed higher-order shell element↩︎

The three-node higher-order shell element is based on the membrane element with drilling DoFs developed in [28], [31], [32], and the Kirchhoff plate element developed in [23], [27]. Figure 1 shows the element and the nodal DoFs. A local coordinate system with an origin at the centroid is defined. The local axes are represented by \(x\), \(y\) and \(z\), to distinguish them from the global ones, \(X\), \(Y\) and \(Z\). The area of the element is \(A^{e}\). The anti-clockwise coordinate along the element boundary is denoted by \(s\), and \(n\) is the exterior normal. The angle between the normal \(n\) and the local axis \(x\) is \(\gamma\).

The in-plane displacements in the \(x\) and \(y\) directions defined over \(A^{e}\) are \(u(x, y)\) and \(v(x, y)\), respectively. Independent in-plane displacements \(\overline{u}_n(s)\) and \(\overline{u}_s(s)\) are assumed along the boundary \(\partial A^{e}\). The membrane DoFs are the nodal displacements \(\overline{u}\), \(\overline{v}\) and the nodal rotation \(\overline{\theta}_z\). The out-of-plane plate displacement defined over the domain is \(w(x, y)\). An independent out-of-plane boundary displacement, \(\overline{w}(s)\), and its compatible normal derivative, \(\partial \overline{w}/\partial n(s)\), are assumed along the boundary \(\partial A^{e}\). The plate DoFs at each node include the displacement \(\overline{w}\), and the two rotations \(\partial \overline{w}/\partial x =-\overline{\theta}_y\) and \(\partial \overline{w}/\partial y =\overline{\theta}_x\), defined according to Kirchhoff kinematics.

2.1 Membrane element with drilling DoFs↩︎

Consider a membrane surface domain \(A\) with boundary \(\partial A: \partial A_{\color{blue}N} \cup \partial A_{\mathrm{d}}\) in a plane stress condition. Membrane forces \(\mathbf{N}^{*} = \bigl\lfloor N_{nn}^{*},\;N_{ns}^{*} \bigr\rfloor^{\mathrm{T}}\) are prescribed on \(\partial A_{\color{blue}N}\), whereas in-plane displacements \(\mathbf{d}^{*}=\bigl\lfloor u_{n}^{*},\;u_{s}^{*} \bigr\rfloor^{\mathrm{T}}\) are prescribed on \(\partial A_{\mathrm{d}}\). The internal fields are in-plane displacements \(\mathbf{u}=\bigl\lfloor u,\;v \bigr\rfloor^{\mathrm{T}}\), membrane forces \(\mathbf{N}=\bigl\lfloor N_{\text{xx}},\;N_{\text{yy}},\;N_{\text{xy}} \bigr\rfloor^{\mathrm{T}}\), membrane strains \(\boldsymbol{\varepsilon}_{m}=\bigl\lfloor \varepsilon_{\text{xx}},\;\varepsilon_{\text{yy}},\;\gamma_{\text{xy}} \bigr\rfloor^{\mathrm{T}}\), and given forces per unit area \(\mathbf{q}^{*}= \bigl\lfloor q_{\text{x}}^{*},\;q_{\text{y}}^{*}\bigr\rfloor^{\mathrm{T}}\). The internal field equations are: \[\begin{align} {\color{blue}\boldsymbol{\varepsilon}_{\mathrm{m}}=\boldsymbol{\nabla}_{\mathrm{m}}\mathbf{u},\boldsymbol{N}=\mathbf{A} \boldsymbol{\varepsilon}_{\mathrm{m}},\boldsymbol{\nabla}_{\mathrm{m}}^{\mathrm{T}} \mathbf{N}+\mathbf{q}^{*}=\mathbf{0}\text{in A}} \end{align}\] where \[\begin{align} \nonumber \boldsymbol{\nabla}^{\mathrm{T}}_{\mathrm{m}}=\left[\begin{array}{ccc} \displaystyle \frac{\partial(\cdot)}{\partial x} & 0 & \displaystyle \frac{\partial(\cdot)}{\partial y} \\[1.0em] 0 & \displaystyle \frac{\partial(\cdot)}{\partial y} & \displaystyle \frac{\partial(\cdot)}{\partial x} \end{array}\right],\text{and}\mathbf{A}=\left[\begin{array}{ccc} A_{11} & A_{12} & A_{16} \\ A_{12} & A_{22} & A_{26} \\ A_{16} & A_{26} & A_{66} \end{array}\right] \end{align}\]

are the membrane strain operator and the membrane stiffness matrix, respectively. The hybrid/mixed parametrized energy functional of \(n\) membrane elements connected through their boundary displacements, as given by [32], can be written as \[\begin{align} \nonumber \Psi_{\mathrm{m}}(\tilde{\mathbf{u}}, \tilde{\mathbf{N}}, \bar{\mathbf{d}}) = \sum_{e=1}^{n} \left\{\frac{1}{2}(1-\gamma)\int_{A^{e}} \mathbf{N}^{u\mathrm{T}} \boldsymbol{\varepsilon}_{\mathrm{m}}^u \mathrm{~d} A -\frac{1}{2} \gamma \int_{A^{e}} \tilde{\mathbf{N}}^{T} \boldsymbol{\varepsilon}_{\mathrm{m}}^N \mathrm{~d} A \right.\\ \left.+\gamma \int_{A^{e}}\tilde{\mathbf{N}}^{T} \boldsymbol{\varepsilon}_{\mathrm{m}}^u \mathrm{~d} A-\mathbb{P}^{e}_{\mathrm{m}}\right\}, \end{align}\] where \(\gamma\) is a scalar, and \(\mathbb{P}^{e}_{\mathrm{m}}(\tilde{\mathbf{u}}, \tilde{\mathbf{N}}, \bar{\mathbf{d}})\) is the membrane forcing potential, analogous to the forcing potential of the d-generalized variational principle given by [32], which enforces in a weak sense \(\tilde{\mathbf{u}}\) and \(\bar{\mathbf{d}}\) to be equal at \(\partial A\).

For the finite element discretization, it is assumed that \[\begin{align} \tilde{\mathbf{u}}=\boldsymbol{\phi}_{u}\mathbf{q}_{u},\tilde{\mathbf{N}}=\boldsymbol{\phi}_{N}\mathbf{q}_{N}\text{in}A^{e}, \quad \text{and} \quad \bar{\mathbf{d}}=\boldsymbol{\phi}_{d} \mathbf{q}_{d}\text{in}\partial A^{e}. \label{Membrane95Approximations} \end{align}\tag{1}\] \[\begin{align} {\color{blue}\tilde{\mathbf{u}}=\left\{\begin{array}{l} \tilde{u} \\ \tilde{v} \end{array}\right\}=\sum_{i=1}^9 \boldsymbol{\phi}_{u_i} q_{u_i}=\boldsymbol{\phi}_{u}\mathbf{q}_{u}=\boldsymbol{\phi}_{r} \mathbf{q}_{\mathrm{r}}+\boldsymbol{\phi}_{c} \mathbf{q}_{\mathrm{c}}+\boldsymbol{\phi}_{h} \mathbf{q}_{\mathrm{h}},} \label{Displacement95decomposition} \end{align}\tag{2}\] where the rigid body, constant strain and higher-order modes \(\boldsymbol{\phi}_{r}\), \(\boldsymbol{\phi}_{c}\) and \(\boldsymbol{\phi}_{h}\) are given in [28]. The internal fields \(\tilde{\mathbf{u}}\) and \(\tilde{\mathbf{N}}\) may be discontinuous across elements. On the other hand, the boundary displacement field \(\bar{\mathbf{d}}\) have the same value on adjacent elements, as the approximation adopted for \(\bar{\mathbf{d}}\) on an edge separating two elements is uniquely interpolated by nodal values \(\mathbf{q}_{d}=\overline{\mathbf{u}}\): \[\begin{align} {\color{blue}\overline{\mathbf{u}}=\left\{\overline{u}_1, \overline{v}_1, \overline{\theta}_{z_1},\overline{u}_2, \overline{v}_2, \overline{\theta}_{z_2}, \overline{u}_3, \overline{v}_3, \overline{\theta}_{z_3} \right\}^{\mathrm{T}},} \label{Membrane95Dofs} \end{align}\tag{3}\] where \(\theta_{z}\) is given in geometrically linear fashion as \(\theta_{z}=\displaystyle \frac{1}{2}\left(\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}\right)\).

\[\begin{align} {\color{blue}\overline{\mathbf{u}}=\mathbf{G} \mathbf{q}_{u}= \left[\begin{array}{c} \mathbf{G}_{r} \quad \mathbf{G}_{c} \quad \mathbf{G}_{h} \end{array}\right] \mathbf{q}_{u}= \mathbf{G}_{r} \mathbf{q}_{r}+\mathbf{G}_{c} \mathbf{q}_{c}+\mathbf{G}_{h} \mathbf{q}_{h},} \label{au95to95qm95transformation} \end{align}\tag{4}\] where the (\(9\times3\)) matrices \(\mathbf{G}_{r}\), \(\mathbf{G}_{c}\) and \(\mathbf{G}_{h}\) are given in [28]. The resulting \(9\times9\) matrix \(\mathbf{G}\) is non-singular and may be inverted to give \[\begin{align} {\color{blue}\mathbf{q}_{u}= \left\{\begin{array}{ccc} \mathbf{q}_{r} & \mathbf{q}_{c} & \mathbf{q}_{h} \end{array}\right\}^{\mathrm{T}}= \mathbf{G}^{-1} \bar{\mathbf{u}}= \left[\begin{array}{ccc} \mathbf{H}_{r} & \mathbf{H}_{c} & \mathbf{H}_{h} \end{array}\right]^{\mathrm{T}} \bar{\mathbf{u}}.} \label{qm95to95au95transformation} \end{align}\tag{5}\]

Inserting the approximations into the functional \(\Psi_{\mathrm{m}}^{e}\) of an element, and making it stationary, yields the element equations. The internal displacement decomposition in Eq. (5 ) induces a partitioned version of the element equations, which can be reduced by using Eq. (2 ) and static condensation, as discussed in detail in [32]. The reduced element equations are given by \[\begin{align} {\color{blue}\left[\mathbf{K}_{\mathrm{b}}+(1-\gamma) \mathbf{K}_{\mathrm{h}}\right] \bar{\mathbf{u}}=\mathbf{f}_{\mathrm{m}} \quad \Rightarrow \quad \mathbf{K}_{\mathrm{m}} \bar{\mathbf{u}}=\mathbf{f}_{\mathrm{m}}} \label{Final95Membrane95System} \end{align}\tag{6}\]

\[\begin{align} {\color{blue}\mathbf{K}_{\mathrm{b}}=\frac{1}{A^{e}} \overline{\mathbf{L}} \mathbf{A} \overline{\mathbf{L}}^{\mathrm{T}}\text{and}\mathbf{K}_{\mathrm{h}}=\mathbf{H}_{\mathrm{h}}^{\mathrm{T}} \mathbf{K}_{qh} \mathbf{H}_{\mathrm{h}},\text{with}\mathbf{K}_{qh}=\int_{A^{e}} \mathbf{B}_{\mathrm{h}}^{\mathrm{T}} \mathbf{A} \mathbf{B}_{\mathrm{h}} \mathrm{~d} A.} \label{K95b95and95K95qh} \end{align}\tag{7}\] Expressions for the lumping matrix \(\overline{\mathbf{L}}\) and the generalized higher-order stiffness matrix \(\mathbf{K}_{qh}\) were explicitly derived in [28]. The external force vector of the membrane element is \[\begin{align} {\color{blue}\mathbf{f}_{\mathrm{m}} = \mathbf{f}_{N^{*}}+\mathbf{H}_{\mathrm{r}}^{\mathrm{T}} \mathbf{f}_{qr^{*}}+\frac{1}{A^{e}} \overline{\mathbf{L}} \mathbf{f}_{qc^{*}}+\mathbf{H}_{\mathrm{h}}^{\mathrm{T}} \mathbf{f}_{qh^{*}}.} \label{Membrane95Stiffnes95Matrices} \end{align}\tag{8}\]

2.2 Kirchhoff-plate element↩︎

Consider a plate surface domain \(A\) with boundary \(\partial A: \partial A_{V_n} \cup \partial A_{w}\) or \(\partial A:\partial A_{M_{n}} \cup \partial A_{w_{n}}\). Kirchhoff shear forces \(V_n^{*}\) and bending moments \(M_{nn}^{*}\) are prescribed on boundaries of the type \(\partial A_{V_n}\) and \(\partial A_{M_{n}}\), respectively, whereas out-of-plane displacements \(w^{*}\) and normal rotations \(\partial w^{*}/\partial n\) are prescribed on boundaries of the type \(\partial A_{w}\) and \(\partial A_{w_{n}}\), respectively. In addition, \(R_{N}^{*}\) are prescribed values of possible concentrated forces. The internal fields are the out-of-plane displacement \(w\), bending moments \(\mathbf{M}=\bigl\lfloor M_{\text{xx}},\;M_{\text{yy}},\;M_{\text{xy}} \bigr\rfloor^{\mathrm{T}}\), curvatures \(\boldsymbol{\kappa}=\bigl\lfloor \kappa_{\text{xx}},\;\kappa_{\text{yy}},\;\kappa_{\text{xy}} \bigr\rfloor^{\mathrm{T}}\), and given out-of-plane forces per unit area \(q_z^{*}\). The internal field equations are: \[\begin{align} {\color{blue}\boldsymbol{\kappa}=\boldsymbol{\nabla}_{\mathrm{p}}w,\boldsymbol{M}=\mathbf{D} \boldsymbol{\kappa},\boldsymbol{\nabla}_{\mathrm{p}}^{\mathrm{T}} \mathbf{M}+q_{z}^{*}=0 \text{in A}} \end{align}\] where \[\begin{align} \nonumber \boldsymbol{\nabla}^{\mathrm{T}}_{\mathrm{p}}=\left[ \begin{array}{ccc} \displaystyle \frac{\partial^{2}(\cdot)}{\partial x^{2}} & \displaystyle \frac{\partial^{2}(\cdot)}{\partial y^{2}} & 2\displaystyle\frac{\partial^{2}(\cdot)}{\partial x\,\partial y} \end{array} \right] \quad \text{and} \quad \mathbf{D}=\left[\begin{array}{ccc} D_{11} & D_{12} & D_{16} \\ D_{12} & D_{22} & D_{26} \\ D_{16} & D_{26} & D_{66} \end{array}\right] \end{align}\]

are the Kirchhoff-plate differential operator and the plate stiffness matrix. The hybrid energy functional of \(n\) plate elements connected through their boundary displacements can be written as [27] \[\begin{align} \nonumber \Psi_{\mathrm{p}}\left(\tilde{w}, \bar{w} \right) =\sum_{e=1}^{n} \left\{\int_{A^{e}} U^{w}_0 \mathrm{~d} A +\sum_{N=1}^{3}R^w_{N}\left(\bar{w}_N-\tilde{w}_{N}\right)\right. \\ \left.+\int_{\partial A^{e}} V^w_{n}\left(\bar{w}-\tilde{w}\right) \mathrm{~d} s - \int_{\partial A^{e}} M^w_{nn}\left(\frac{\partial \bar{w}}{\partial n}-\frac{\partial \tilde{w}}{\partial n}\right) \mathrm{~d} A -\mathbb{P}^{e}_{\mathrm{p}}\right\}, \label{Plate95Energy95Functional} \end{align}\tag{9}\] where \(\tilde{w}\) and \(\bar{w}\) are the domain and boundary out-of-plane plate displacements, respectively.

The strain energy density \(U^{w}_0\), computed from \(\tilde{w}\) for symmetric laminates, and the forcing potential \(\mathbb{P}^{e}_{\mathrm{p}}\) for the plate are presented in [23]. The bending moments \(\mathbf{M}^w=\bigl\lfloor M^w_{\text{xx}},\;M^w_{\text{yy}},\;M^w_{\text{xy}} \bigr\rfloor^{\mathrm{T}}\) give rise to a normal bending moment \(M^w_{nn}\), and resultant Kirchhoff shear force \(V^w_{n}\), on the element boundary, together with concentrated forces \(R^{w}_{N\,(N=1,2,3)}\) at the vertices [27].

For the finite element discretization, it is assumed that

\[\begin{align} \nonumber {\color{blue}\tilde{w}=A_1+A_2 x+A_3 y+\alpha_1 x^2+\alpha_2 x y+\alpha_3 y^2+\alpha_4 x^3}\\ {\color{blue}+\alpha_5 x^2 y+\alpha_6 x y^2+\alpha_7 y^3} \label{Plate95Approximations} \end{align}\tag{10}\]

\[\begin{align} \nonumber {\color{blue}\Psi_{\mathrm{p}}\left(\tilde{w}, \bar{w} \right) =\sum_{m} \left\{-\int_{A^{e}} U^{w}_0 \mathrm{~d} A +\sum_{N=1}^{3}R^w_{N}\bar{w}_N +\int_{\partial A^{e}} V^w_{n}\bar{w} \mathrm{~d} s\right.} \\ {\color{blue}\left.-\int_{\partial A^{e}} M^w_{nn}\frac{\partial \bar{w}}{\partial n} \mathrm{~d} s -\sum_{N=1}^{3}R_{N}^{*}\bar{w}_N -\int_{\partial A^{e}}V_n^{*} \bar{w} \mathrm{~d} s +\int_{\partial A^{e}} M_{nn}^{*} \frac{\partial \bar{w}}{\partial n} \mathrm{~d} s\right\},} \label{Modified95Plate95Energy95Functional} \end{align}\tag{11}\] \[\begin{align} {\color{blue}\int_{A^{e}} U^{w}_0 \mathrm{~d} x \mathrm{~d} y=\frac{1}{2} \boldsymbol{\alpha}^{\mathrm{T}} \mathbf{H} \boldsymbol{\alpha},} \label{Discrete95Strain95Energy} \end{align}\tag{12}\] \[\begin{align} {\color{blue}\mathbf{Q}^{w}=\mathbf{B}^{\mathrm{T}} \boldsymbol{\alpha},} \label{Generalized95Forces95alpha95relation} \end{align}\tag{13}\]

\[\begin{align} {\color{blue}\mathbf{Q}^{w}=\left\{R_1, R_2, R_3, V_n^{12}, V_n^{23}, V_n^{31}, M_{nn}^{12}, M_{nn}^{21}, M_{nn}^{23}, M_{nn}^{32}, M_{nn}^{31}, M_{nn}^{13}\right\}^{\mathrm{T}}} \label{Generalized95Forces} \end{align}\tag{14}\]

\[\begin{align} {\color{blue}\sum_{N=1}^3 R^{w}_N \bar{w}_N+\int_{\partial A^{e}} V^{w}_n \bar{w} \mathrm{~d} s-\int_{\partial A^{e}} M^{w}_{nn} \frac{\partial \bar{w}}{\partial n} \mathrm{~d} s=\mathbf{q}^{\mathrm{T}}\mathbf{Q}^{w}.} \label{Generalized95Forces95Work} \end{align}\tag{15}\]

\[\begin{align} {\color{blue}\overline{\mathbf{w}}=\left\{\bar{w}_1, \frac{\partial \bar{w}_1}{\partial x}, \frac{\partial \bar{w}_1}{\partial y}, \bar{w}_2, \frac{\partial \bar{w}_2}{\partial x}, \frac{\partial \bar{w}_3}{\partial y}, \bar{w}_3, \frac{\partial \bar{w}_3}{\partial x}, \frac{\partial \bar{w}_3}{\partial y}\right\}^{\mathrm{T}}} \label{plate95Dofs} \end{align}\tag{16}\] to the vector \(\mathbf{q}\) by means of: \(\mathbf{q}=\mathbf{T} \overline{\mathbf{w}}\). The derivation of \(\mathbf{T}\) assumes cubic approximation for \(\bar{w}\) and linear variation for \(\partial \bar{w} / \partial n\) along each side \(i\)-\(j\). Substituting \(\mathbf{q}=\mathbf{T} \overline{\mathbf{w}}\) and Eqs. (12 )-(15 ) into Eq. (11 ), the plate energy functional for a finite element under prescribed boundary loads is \[\begin{align} {\color{blue}\Psi^{e}_{\mathrm{p}}=-\frac{1}{2} \boldsymbol{\alpha}^{\mathrm{T}} \mathbf{H} \boldsymbol{\alpha}+\boldsymbol{\alpha}^{\mathrm{T}}(\mathbf{B T}) \overline{\mathbf{w}}-\mathbf{Q}^{* \mathrm{~T}} \mathbf{T} \overline{\mathbf{w}},} \label{Discrete95Plate95Energy95Functional} \end{align}\tag{17}\]

\[\begin{align} {\color{blue}\boldsymbol{\alpha}=\mathbf{H}^{-1}(\mathbf{B}\mathbf{T}) \overline{\mathbf{w}}} \label{alpha95Dofs95Relation} \end{align}\tag{18}\] by setting the coefficient of the arbitrary variation to zero [23]. Performing the variation of \(U^{e}\), with \(\boldsymbol{\alpha}\) substituted by Eq. (18 ), gives \[\begin{align} {\color{blue}\delta U^{e}=\delta \overline{\mathbf{w}}^{\mathrm{T}} \mathbf{K}_{\mathrm{p}} \overline{\mathbf{w}}, \quad \Rightarrow \quad \mathbf{K}_{\mathrm{p}}=(\mathbf{B T})^{\mathrm{T}} \mathbf{H}^{-1}(\mathbf{B}\mathbf{T}),} \label{Discrete95Internal95Work95Variation} \end{align}\tag{19}\]

\[\begin{align} {\color{blue}\delta W^{e}=\delta \overline{\mathbf{w}}^{\mathrm{T}} \mathbf{T}^{\mathrm{T}} \mathbf{Q}^*=\delta \overline{\mathbf{w}}^{\mathrm{T}} \mathbf{f}_{\mathrm{p}}, \quad \Rightarrow \quad \mathbf{f}_{\mathrm{p}}=\mathbf{T}^{\mathrm{T}} \mathbf{Q}^*.} \label{Discrete95External95Work95Variation} \end{align}\tag{20}\]

2.3 Shell element↩︎

For symmetric laminates, membrane–bending coupling is absent. Hence, the membrane and plate stiffness matrices can be assembled directly into the stiffness matrix of the flat shell element \[\begin{align} {\color{blue}\mathbf{K}_{\mathrm{s}}= \left[\begin{array}{cc} \mathbf{K}_{\mathrm{m}} & 0 \\ 0 & \mathbf{K}_{\mathrm{p}} \end{array}\right].} \label{Shell95Stiffness95Matrix} \end{align}\tag{21}\]

\[\begin{align} \nonumber {\color{blue}\mathbf{e}_{3}=\frac{\mathbf{x}_{12} \times \mathbf{x}_{13}}{ \left| \left| \mathbf{x}_{12} \times \mathbf{x}_{13} \right| \right|}, \quad \mathbf{e}_{1} = \frac{\mathbf{v}-(\mathbf{v} \cdot \mathbf{e}_{3})\mathbf{e}_{3}}{\left| \left|\mathbf{v}-(\mathbf{v} \cdot \mathbf{e}_{3})\mathbf{e}_{3}\right| \right|}, \quad \mathbf{e}_{2} = \mathbf{e}_{3} \times \mathbf{e}_{1},} \end{align}\] \[\begin{align} \nonumber {\color{blue}\mathbf{q}^{l}_{i}=\mathbf{L}_{q} \mathbf{q}^{g}_{i}, \quad \quad \quad \mathbf{q}^{g}_{i}=\mathbf{L}^{\mathrm{T}}_{q} \mathbf{q}^{l}_{i}, \quad \quad \quad \mathbf{L}_{q}=\operatorname{diag} \left[\mathbf{L}_{l}, \quad \mathbf{L}_{w\theta} \right]} \end{align}\] \[\begin{align} \nonumber {\color{blue}\mathbf{q}^{l}_{\mathrm{s}}=\mathbf{R} \mathbf{q}^{g}_{\mathrm{s}}, \quad \quad \quad \mathbf{q}^{g}_{\mathrm{s}}=\mathbf{R}^{\mathrm{T}} \mathbf{q}^{l}_{\mathrm{s}}, \quad \quad \quad \mathbf{R}=\operatorname{diag} \left[\mathbf{L}_{q}, \quad \mathbf{L}_{q}, \quad \mathbf{L}_{q} \right]} \end{align}\]

3 The shell-to-shell cohesive line interface model↩︎

3.1 Kinematics: Displacement jump and DoFs↩︎

The shell-to-shell cohesive line element must be kinematically compatible with two shell elements in the panel (\(p_1\) and \(p_2\)) and the stiffener element (\(s\)), as illustrated in Figure 2. The panel and the stiffener have constant thickness \(t^{p}\) and \(t^{s}\), respectively. The displacement jump over the cohesive surface \(\Omega_c=\left[-t^{s}/2,t^{s}/2\right] \times \Gamma_c\) is computed from the displacement jump over its central line \(\Gamma_c\), in addition to a contribution from the jump in rotations over \(\Gamma_c\)

\[\boldsymbol{\Delta}\left(\mathbf{x}\right)= \llbracket \mathbf{u}_{\Gamma_c} \rrbracket + z^{c} \llbracket \boldsymbol{\theta}_{\Gamma_c} \rrbracket,\mathbf{x} \in \Omega_c \label{JumpDecomposition}\tag{22}\]

where \(\boldsymbol{\Delta}\) is the displacement jump vector, written with respect to the cohesive local system \(\mathbf{e}_{i}^{c}\), and \(\llbracket \mathbf{u}_{\Gamma_c} \rrbracket\), \(\llbracket \boldsymbol{\theta}_{\Gamma_c} \rrbracket\) are, respectively, displacement and rotation jumps over the central line \(\Gamma_c\), which are defined as \[\llbracket \mathbf{u}_{\Gamma_c} \rrbracket=\mathbf{u}^{c+}_{\Gamma_c}-\mathbf{u}^{c-}_{\Gamma_c},\llbracket \boldsymbol{\theta}_{\Gamma_c} \rrbracket=\boldsymbol{\theta}^{c+}_{\Gamma_c}-\boldsymbol{\theta}^{c-}_{\Gamma_c}. \label{DispRotJumps}\tag{23}\]

The central line displacements \(\mathbf{u}^{c+}_{\Gamma_c}\) and \(\mathbf{u}^{c-}_{\Gamma_c}\) are given by \[\begin{align} \mathbf{u}^{c+}_{\Gamma_c} = \mathbf{L}_{cs}\mathbf{u}^{s}_{\Gamma_c}, \quad \mathbf{u}^{c-}_{\Gamma_c} = \mathbf{L}_{cp}\mathbf{u}^{p}_{\Gamma_c} \label{CentralDisplacementDefinition} \end{align}\tag{24}\] where the transformations \[\mathbf{L}_{c p} = \left[\begin{array}{lll} \mathbf{e}_{1}^{c} \cdot \mathbf{e}_{1}^{p} & \mathbf{e}_{1}^{c} \cdot \mathbf{e}_{2}^{p} & \mathbf{e}_{1}^{c} \cdot \mathbf{e}_{3}^{p}\\ \mathbf{e}_{2}^{c} \cdot \mathbf{e}_{1}^{p} & \mathbf{e}_{2}^{c} \cdot \mathbf{e}_{2}^{p} & \mathbf{e}_{2}^{c} \cdot \mathbf{e}_{3}^{p}\\ \mathbf{e}_{3}^{c} \cdot \mathbf{e}_{1}^{p} & \mathbf{e}_{3}^{c} \cdot \mathbf{e}_{2}^{p} & \mathbf{e}_{3}^{c} \cdot \mathbf{e}_{3}^{p} \end{array}\right],\mathbf{L}_{c s} = \left[\begin{array}{lll} \mathbf{e}_{1}^{c} \cdot \mathbf{e}_{1}^{s} & \mathbf{e}_{1}^{c} \cdot \mathbf{e}_{2}^{s} & \mathbf{e}_{1}^{c} \cdot \mathbf{e}_{3}^{s}\\ \mathbf{e}_{2}^{c} \cdot \mathbf{e}_{1}^{s} & \mathbf{e}_{2}^{c} \cdot \mathbf{e}_{2}^{s} & \mathbf{e}_{2}^{c} \cdot \mathbf{e}_{3}^{s}\\ \mathbf{e}_{3}^{c} \cdot \mathbf{e}_{1}^{s} & \mathbf{e}_{3}^{c} \cdot \mathbf{e}_{2}^{s} & \mathbf{e}_{3}^{c} \cdot \mathbf{e}_{3}^{s} \end{array}\right] \label{TransformationMatrices}\tag{25}\] are computed from the basis vectors: \(\mathbf{e}_{i}^{p}\), \(\mathbf{e}_{i}^{s}\) and \(\mathbf{e}_{i}^{c}\) in Fig. 2. The stiffener and panel displacement fields along \(\Gamma_c\) are defined based on the classical Kirchhoff plate kinematics \[\begin{align} & \mathbf{u}^{s}_{\Gamma_c}=\left\{\begin{array}{c} u^{s}|_{\Gamma_c} \\ v^{s}|_{\Gamma_c} \\ w^{s}|_{\Gamma_c} \end{array}\right\},\mathbf{u}^{p}_{\Gamma_c}=\left\{\begin{array}{c} u^{p}|_{\Gamma_{c}^{p}} \\ v^{p}|_{\Gamma_{c}^{p}} \\ w^{p}|_{\Gamma_{c}^{p}} \end{array}\right\}+\frac{t^{p}}{2}\left\{\begin{array}{c} -w_{\color{blue},x}^{p}|_{\Gamma_{c}^{p}} \\ -w_{\color{blue},y}^{p}|_{\Gamma_{c}^{p}} \\ 0 \end{array}\right\}, \end{align}\] in which \((\cdot)|_{\Gamma}\) stands for the internal fields evaluated at the edge \(\Gamma\), and the notation \(f_{\color{blue},x} = \partial f/\partial x\), \(f_{\color{blue},y} = \partial f/\partial y\) is adopted herein for derivatives.

The central line rotations \(\boldsymbol{\theta}^{c+}_{\Gamma_c}\) and \(\boldsymbol{\theta}^{c-}_{\Gamma_c}\) which cause motion out-of the stiffener plane are \[\begin{align} \tag{26} \boldsymbol{\theta}^{c+}_{\Gamma_c} &={\color{blue}\boldsymbol{\Theta}^{c+}_{\Gamma_c} \mathbf{e}_{3}} =\left\{\begin{array}{c} \phantom{-}\theta_x^{s}|_{\Gamma_c} \\ -\theta_y^{s}|_{\Gamma_c} \\ 0 \end{array}\right\}=\left\{\begin{array}{c} w_{\color{blue},y}^{s}|_{\Gamma_c} \\ w_{\color{blue},x}^{s}|_{\Gamma_c} \\ 0 \end{array}\right\}, \\ \tag{27} \boldsymbol{\theta}^{c-}_{\Gamma_c}&= {\color{blue}\boldsymbol{\Theta}^{c-}_{\Gamma_c} \mathbf{e}_{3}} =\left\{\begin{array}{c} \phantom{-}\theta_x^{p}|_{\Gamma_{c}} \\ -\theta_z^{p}|_{\Gamma_{c}} \\ 0 \end{array}\right\}=\left\{\begin{array}{c} w_{\color{blue},y}^{p}|_{\Gamma_{c}^{p}} \\ -\frac{1}{2}\left(v_{\color{blue},x}^{p}-u_{\color{blue},y}^{p}\right)|_{\Gamma_{c}^{p}} \\ 0 \end{array}\right\}. \end{align}\] where the skew-symmetric small rotation tensors (or spin tensors) \({\color{blue}\boldsymbol{\Theta}^{c+}_{\Gamma_c}}=\mathbf{L}_{cs}{\color{blue}\boldsymbol{\Theta}^{s}_{\Gamma_c}}\mathbf{L}^{T}_{cs}\) and \({\color{blue}\boldsymbol{\Theta}^{c-}_{\Gamma_c}}=\mathbf{L}_{cp}{\color{blue}\boldsymbol{\Theta}^{p}_{\Gamma_c}}\mathbf{L}^{T}_{cp}\) are computed with the transformation matrices (25 ), defined for the basis vectors illustrated in Fig. 2. Notice that \(\theta_x=w_{,y}\) and \(\theta_z=-1/2\left[v_{,x}-zw_{,yx}-\left(u_{,y}-zw_{,xy}\right)\right]\) in Eq.(27 ) do not depend on the local coordinate \(z\).

Considering Eqs. (23 )-(24 ) and (26 )-(27 ), the central line displacement and rotation jumps are given by \[\begin{align} \tag{28} &\llbracket \mathbf{u}_{\Gamma_c} \rrbracket =\left\{\begin{array}{c} v^{s}|_{\Gamma_c} - w^{p}|_{\Gamma_c^{p}} \\[4pt] u^{s}|_{\Gamma_c} - u^{p}|_{\Gamma_c^{p}}+\frac{t^{p}}{2}w_{\color{blue},x}^{p}|_{\Gamma_c^{p}} \\[4pt] -w^{s}|_{\Gamma_c}-v^{p}|_{\Gamma_c^{p}}+\frac{t^{p}}{2}w_{\color{blue},y}^{p}|_{\Gamma_c^{p}} \end{array}\right\},\\ &\llbracket \boldsymbol{\theta}_{\Gamma_c} \rrbracket =\left\{\begin{array}{c} w_{\color{blue},y}^{s}|_{\Gamma_c} - w_{\color{blue},y}^{p}|_{\Gamma_c} \\[4pt] w_{\color{blue},x}^{s}|_{\Gamma_c} + \frac{1}{2}\left(v_{\color{blue},x}^{p}-u_{\color{blue},y}^{p}\right)|_{\Gamma_{c}^{p}} \\[4pt] 0 \end{array}\right\}. \tag{29} \end{align}\]

Due to the incompatible internal fields of the adopted shell element, the panel fields \(u^{p}\), \(v^{p}\), \(v_{\color{blue},x}^{p}\), \(u_{\color{blue},y}^{p}\), \(w^{p}\), \(w_{\color{blue},x}^{p}\), \(w_{\color{blue},y}^{p}\) are herein defined as the average of the internal fields from panels 1 and 2, evaluated at their common edge. Notice that the local frame of reference of panel \(p_2\) corresponds to the local frame of \(p_1=p\), with its origin translated to the centroid of \(p_2\).

Evaluating the internal fields of the shell element along its edges requires the contribution of the DOFs associated with all its nodes. Therefore, the DoF vector of the cohesive element is defined as: \(\mathbf{q}_{\mathrm{c}}=\left[\overline{\mathbf{u}}_{p_1}^{\mathrm{T}}, \overline{\mathbf{w}}_{p_1}^{\mathrm{T}}, \overline{\mathbf{u}}_{p_2}^{\mathrm{T}}, \overline{\mathbf{w}}_{p_2}^{\mathrm{T}}, \overline{\mathbf{u}}_{s}^{\mathrm{T}}, \overline{\mathbf{w}}_{s}^{\mathrm{T}}\right]^{\mathrm{T}}\), where \(\overline{\mathbf{u}}_{\alpha}\) and \(\overline{\mathbf{w}}_{\alpha}\), \(\alpha = s, p_1, p_2\), are the membrane and plate DoF vectors for the first panel, second panel and stiffener elements, respectively.

The central task of cohesive element formulations is to find the kinematic operator \(\mathbf{B}_{\mathrm{c}}\) that relates the opening vector to the nodal DoFs \[\boldsymbol{\Delta}\left(\mathbf{x}\right)=\left[\Delta_{\mathrm{I}}, \Delta_{\mathrm{II}}, \Delta_{\mathrm{III}}\right]^{\mathrm{T}}=\mathbf{B}_{\mathrm{c}} \mathbf{q}_{\mathrm{c}}, \quad \mathbf{x} \in \Omega_c . \label{Jump95Dofs95relation95I}\tag{30}\]

Examining the expression of \(\boldsymbol{\Delta}\) in Eq. (22 ), it is possible to write

\[\boldsymbol{\Delta}=\left[\mathbf{B}_{\mathrm{c_d}} + z^{c} \mathbf{B}_{\mathrm{c_\theta}}\right]\mathbf{q}_\mathrm{c} = \mathbf{B}_{\mathrm{c}}\mathbf{q}_{\mathrm{c}} \label{Jump95Dofs95relation95II}\tag{31}\]

in which the operators \(\mathbf{B}_{\mathrm{c_d}}\) and \(\mathbf{B}_{\mathrm{c_\theta}}\) are defined from the shell internal fields, evaluated at \(\Gamma_c\) and \(\Gamma_c^{p}\). Therefore, they are functions solely of the parametric coordinates \(s^{s} \in \Gamma_c\) and \(s^{p} \in \Gamma_c^{p}\). From the expressions of \(\llbracket \mathbf{u}_{\Gamma_c} \rrbracket\) and \(\llbracket \boldsymbol{\theta}_{\Gamma_c} \rrbracket\) in Eqs. (28 )-(29 ), we can see that the \(\mathbf{B}_{\mathrm{c_d}}\) and \(\mathbf{B}_{\mathrm{c_\theta}}\) matrices shall be composed of sub-matrices that relate the approximation terms \[\begin{align} \nonumber &\tilde{u}^{s}, \tilde{v}^{s}, \tilde{w}^{s}, \tilde{u}^{p_1}, \tilde{v}^{p_1}, \tilde{w}^{p_1}, \tilde{w}_{\color{blue},x}^{p_1}, \tilde{w}_{\color{blue},y}^{p_1}, \tilde{u}^{p_2}, \tilde{v}^{p_2}, \tilde{w}^{p_2}, \tilde{w}_{\color{blue},x}^{p_2}, \tilde{w}_{\color{blue},y}^{p_2}, \\ \nonumber & \tilde{w}_{\color{blue},x}^{s}, \tilde{w}_{\color{blue},y}^{s}, \tilde{w}_{\color{blue},y}^{p_1}, \tilde{w}_{\color{blue},y}^{p_2}, \tilde{v}_{\color{blue},x}^{p_1}, \tilde{v}_{\color{blue},x}^{p_2}, \tilde{u}_{\color{blue},y}^{p_1}, \tilde{u}_{\color{blue},y}^{p_2}. \end{align}\] to the nodal DoFs in \(\mathbf{q}_{\mathrm{c}}\).

The in-plane displacements \(\tilde{u}^{\alpha}\), \(\tilde{v}^{\alpha}\) are those from Eq (2 ). These approximations, and their derivatives, can be expressed in terms of the membrane nodal Dofs \(\overline{\mathbf{u}}^{\alpha}\). The out-of-plane displacements \(\tilde{w}^{\alpha}\) and the rotations \(\tilde{w}_{\color{blue},x}^{\alpha}\) and \(\tilde{w}_{\color{blue},y}^{\alpha}\) are those derived from the cubic \(\tilde{w}^{\alpha}\) approximation in Eq. (10 ), which were explicitly expressed in terms of the plate nodal Dofs \(\overline{\mathbf{w}}^{\alpha}\) in [23]. All these approximations and the corresponding shape function matrices \(\mathbf{N}_{u}^{\alpha}\), \(\mathbf{N}_{v}^{\alpha}\), \(\mathbf{N}_{v_{\color{blue},x}}^{\alpha}\), \(\mathbf{N}_{u_{\color{blue},y}}^{\alpha}\), \(\mathbf{N}_{w}^{\alpha}\), \(\mathbf{N}_{w_{\color{blue},x}}^{\alpha}\), \(\mathbf{N}_{w_{\color{blue},y}}^{\alpha}\) are summarized in 6.

3.2 The \(\mathbf{B}_{\mathrm{c}}\) matrix↩︎

Once the shape functions have been explicitly defined in 6, the next step is to assemble the \(\mathbf{B}_{\mathrm{c}}\) matrix that relates the DoF vector \(\mathbf{q}_{\mathrm{c}}\) to the opening vector \(\boldsymbol{\Delta}\) \[\mathbf{B}_{\mathrm{c}}=\mathbf{B}_{\mathrm{c_d}}\left(s,s^{p}\right) + z^{c} \mathbf{B}_{\mathrm{c_\theta}}\left(s,s^{p}\right), \label{BCE95Matrix}\tag{32}\] with \(s \in \Gamma_c\), \(s^{p} \in \Gamma_c^{p}\) and \(z^{c} \in \left[-t^{s}/2,t^{s}/2\right]\).

Substituting the approximations in Eqs. (40 )-(45 ) to express the \(u\), \(v\), \(v_{\color{blue},x}\), \(u_{\color{blue},y}\), \(w\), \(w_{\color{blue},x}\) and \(w_{\color{blue},y}\) terms in Eqs. (28 )-(29 ), the expressions for \(\mathbf{B}_{\mathrm{c}_{d}}\) and \(\mathbf{B}_{\mathrm{c}_{\theta}}\) results \[\begin{align} \nonumber \mathbf{B}_{\mathrm{c}_{d}}= \left[\begin{array}{cccccc} \mathbf{0} & -\frac{1}{2}\mathbf{N}_{w}^{p_1}|_{\Gamma_c^{p}} & \mathbf{0} & -\frac{1}{2}\mathbf{N}_{w}^{p_2}|_{\Gamma_c^{p}} & \mathbf{N}_{v}^{s}|_{\Gamma_c} & \mathbf{0}\\[4pt] -\frac{1}{2}\mathbf{N}_{u}^{p_1}|_{\Gamma_c^{p}} & \frac{t^{p}}{4}\mathbf{N}_{w_{\color{blue},x}}^{p_1}|_{\Gamma_c^{p}} &-\frac{1}{2}\mathbf{N}_{u}^{p_2}|_{\Gamma_c^{p}} & \frac{t^{p}}{4}\mathbf{N}_{w_{\color{blue},x}}^{p_2}|_{\Gamma_c^{p}} & \mathbf{N}_{u}^{s}|_{\Gamma_c} & \mathbf{0}\\[4pt] -\frac{1}{2}\mathbf{N}_{v}^{p_1}|_{\Gamma_c^{p}} & \frac{t^{p}}{4}\mathbf{N}_{w_{\color{blue},y}}^{p_1}|_{\Gamma_c^{p}} &-\frac{1}{2}\mathbf{N}_{v}^{p_2}|_{\Gamma_c^{p}} & \frac{t^{p}}{4}\mathbf{N}_{w_{\color{blue},y}}^{p_2}|_{\Gamma_c^{p}} & \mathbf{0} & -\mathbf{N}_{w}^{s}|_{\Gamma_c} \end{array}\right] \end{align}\] and \[\begin{align} \nonumber \mathbf{B}_{\mathrm{c}_{\theta}}=\left[\begin{array}{cccccc} \mathbf{0} & -\frac{1}{2}\mathbf{N}_{w_{\color{blue},y}}^{p_1}|_{\Gamma_c^{p}} & \mathbf{0} & -\frac{1}{2}\mathbf{N}_{w_{\color{blue},y}}^{p_2}|_{\Gamma_c^{p}} & \mathbf{0} & \mathbf{N}_{w_{\color{blue},y}}^{s}|_{\Gamma_c}\\ \frac{1}{2}\mathbf{N}_{\theta_{z}}^{p_1}|_{\Gamma_c^{p}} & \mathbf{0} & \frac{1}{2}\mathbf{N}_{\theta_{z}}^{p_2}|_{\Gamma_c^{p}} & \mathbf{0} & \mathbf{0} & \mathbf{N}_{w_{\color{blue},x}}^{s}|_{\Gamma_c}\\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} \end{array}\right], \end{align}\] in which \(\mathbf{N}_{\theta_{z}}^{p_i}=1/2\left(\mathbf{N}_{v_{\color{blue},x}}^{p_i}-\mathbf{N}_{u_{\color{blue},y}}^{p_i}\right)\), \(i=1,2\), and \(\mathbf{N}|_{\Gamma}\) stands for the shape function matrix \(\mathbf{N}\) evaluated at the edge \(\Gamma\).

3.3 Turon mixed-mode damage cohesive model↩︎

The mixed-mode damage model used in the present work is the well-known Turon model [5]. The relationship between cohesive traction \(\mathbf{t}_{\mathrm{c}}=\left[\mathrm{t}_{\mathrm{I}}, \mathrm{t}_{\mathrm{II}}, \mathrm{t}_{\mathrm{III}}\right]^{\mathrm{T}}\) and the opening vector \(\boldsymbol{\Delta}=\left[\Delta_{\mathrm{I}}, \Delta_{\mathrm{II}}, \Delta_{\mathrm{III}}\right]^{\mathrm{T}}\) are given by

\[\begin{align} {\color{blue}\left\{\begin{array}{c} \mathrm{t}_{\mathrm{I}} \\ \mathrm{t}_{\mathrm{II}} \\ \mathrm{t}_{\mathrm{III}} \end{array}\right\}= \left[\begin{array}{ccc} \left(1-d_{\mathrm{I}}\right) K & 0 & 0 \\ 0 & (1-d) K & 0 \\ 0 & 0 & (1-d) K \end{array}\right]\left\{\begin{array}{c} \Delta_{\mathrm{I}} \\ \Delta_{\mathrm{II}} \\ \Delta_{\mathrm{III}} \end{array}\right\}} \label{Cohesive95Model} \end{align}\tag{33}\]

where \(K\) and \(d\) are the penalty stiffness and the damage variable of the cohesive model, respectively. The damage variable \(d\) in this work is updated by the bi-linear cohesive law proposed by Turon et al. [5], [6]. The damage variable \(d_{\mathrm{I}}\) under Mode I loading is distinguished from \(d\) to avoid interpenetration of the top and bottom surfaces under compression \[\begin{align} d_{\mathrm{I}}= \begin{cases}d, & \Delta_{\mathrm{I}} \geq 0 \\ 0, & \Delta_{\mathrm{I}}<0\end{cases}. \label{Interpenetration95Condition} \end{align}\tag{34}\] The parameters required by Turon’s model are the fracture energies \(G_{\mathrm{I}c}\), \(G_{\mathrm{II}c}\), the interface strengths \(\tau_{\mathrm{I}c}\), \(\tau_{\mathrm{II}c}\), the penalty stiffness \(K\) and the exponent \(\eta\) from the phenomenological relation between fracture toughness and mode ratio proposed by Benzeggagh and Kenane [33], from which the critical energy release rate \(G_{\mathrm{c}}\) for a mixed-mode ratio is given by \[\begin{align} {\color{blue}G_{\mathrm{c}}=G_{\mathrm{Ic}}+\left(G_{\mathrm{IIc}}-G_{\mathrm{Ic}}\right)\left(\frac{G_{\text{shear }}}{G_{\mathrm{T}}}\right)^\eta} \label{BenzeggaghandKenaneRelation} \end{align}\tag{35}\]

The proper choice of the penalty stiffness \(K\) is a key factor in cohesive models. The closed-form expression \[\begin{align} {\color{blue}K=\alpha_{p} \frac{E_3}{t}} \label{Penalty95Stiffness} \end{align}\tag{36}\]

3.4 Cohesive-interface element equation↩︎

The variational formulation of the stiffened panel model that includes the shell-to-shell cohesive line interface model may be obtained by adding the cohesive energy variation

\[\begin{align}\delta \Psi_{\mathrm{c}} &= \int_{\Omega_{c}} \delta {\boldsymbol{\Delta}}^{\mathrm{T}} \mathbf{t}_{\mathrm{c}} \mathrm{~d} \Omega= \int_{\Gamma_{c}} \int_{-\frac{t^{s}}{2}}^{\frac{t^{s}}{2}} \delta {\boldsymbol{\Delta}}^{\mathrm{T}} \mathbf{t}_{\mathrm{c}} \mathrm{~d} z^{c} \mathrm{~d} \Gamma \end{align}\] into the weak form of the equilibrium equations. Introducing the cohesive model (33 ) and the approximations (31 ) into the cohesive virtual work, and considering the arbitrary nature of the virtual DoFs \(\delta \mathbf{q}_{\mathrm{c}}\), the element equations result \[\begin{align}\delta \Psi_{\mathrm{c}} = \delta \mathbf{q}_{\mathrm{c}}^{\mathrm{T}}\mathbf{f}_{\mathrm{c}} =0 \quad \Leftrightarrow \quad \mathbf{f}_{\mathrm{c}}=\mathbf{0} \end{align}\] in which \[\begin{align}\mathbf{f}_{\mathrm{c}} = \int_{\Gamma_{c}} \int_{-\frac{t^{s}}{2}}^{\frac{t^{s}}{2}} \mathbf{B}_{\mathrm{c}}^{\mathrm{T}} \mathbf{t}_{\mathrm{c}} \mathrm{~d} z^{c}\mathrm{~d}\Gamma \label{f95CE} \end{align}\tag{37}\]

is the internal force vector of the cohesive element. The consistent tangent stiffness matrix of the proposed shell-to-shell cohesive line element results \[\begin{align}\mathbf{K}_{\mathrm{c}} &=\frac{\partial \mathbf{f}_{\mathrm{c}}}{\partial \mathbf{q}_{\mathrm{c}}}=\int_{\Gamma_{c}} \int_{-\frac{t^{s}}{2}}^{\frac{t^{s}}{2}} \mathbf{B}_{\mathrm{c}}^{\mathrm{T}} {\color{blue}\mathbf{D}_{\mathrm{c}}}\mathbf{B}_{\mathrm{c}} \mathrm{~d} z^{c} \mathrm{~d} \Gamma, \label{K95CE} \end{align}\tag{38}\]

where \({\color{blue}\mathbf{D}_{\mathrm{c}}}=\partial \mathbf{t}_{\mathrm{c}}/ \partial \boldsymbol{\Delta}\) is the consistent tangent matrix of Turon’s model, as derived in [34]. The consistency of the variational formulation is preserved by using the shell approximation fields in the definition of \(\Delta\).

The integrals in Eqs. (37 ) and (38 ) are difficult to evaluate analytically, mainly due to the complex nonlinear behavior of the cohesive traction vector \(\mathbf{t}_{\mathrm{c}}\) and tangent matrix \(\mathbf{D}_{\mathrm{c}}\) after damage initiation. Thus, numerical integration is applied with a standard Gaussian integration scheme. Earlier works have shown that using a higher number of quadrature points improves the solution accuracy and smoothness in delamination simulations [22]–[24], [35]. Accordingly, numerical integration is performed herein using a 13-point Gauss–Legendre quadrature along \(\Gamma_c\) and a 5-point Gauss–Legendre quadrature along \(z^{c} \in \left[-t^{s}/2,t^{s}/2\right]\).

The ordering of shell DoFs defined in Section 2.3 is also adopted for the DoFs of the panels and stiffener elements: \(\mathbf{q}^{l}_{\text{p}_1}\), \(\mathbf{q}^{l}_{\text{p}_2}\), \(\mathbf{q}^{l}_{\text{s}}\). Therefore \[\begin{align} \nonumber \mathbf{q}_{\mathrm{c}}^{l\mathrm{T}}=\left[\mathbf{q}_{\text{p}_1}^{l\mathrm{T}}, \mathbf{q}_{\text{p}_2}^{l\mathrm{T}}, \mathbf{q}_{\text{s}}^{l\mathrm{T}}\right]^{\mathrm{T}}. \end{align}\]

For this DoF vector, the following transformations hold \[\begin{align} \nonumber \mathbf{q}^{l}_{\mathrm{c}}=\mathbf{R}_{\mathrm{c}} \mathbf{q}^{g}_{\mathrm{c}}, \quad \quad \mathbf{q}^{g}_{\mathrm{c}}=\mathbf{R}_{\mathrm{c}}^{T} \mathbf{q}^{l}_{\mathrm{c}}, \quad \quad \mathbf{R}_{\mathrm{c}}=\left[\begin{array}{lll} \mathbf{R}_{\text{p}_1} & \mathbf{0} & \mathbf{0}\\ \mathbf{0} & \mathbf{R}_{\text{p}_2} & \mathbf{0}\\ \mathbf{0} & \mathbf{0} & \mathbf{R}_{\text{s}} \end{array}\right], \end{align}\]

in which \(\mathbf{q}^{g}_{\mathrm{c}}\) is the DoF vector of the cohesive element in the global coordinate system, and \(\mathbf{R}_{\text{p}_1}\), \(\mathbf{R}_{\text{p}_2}\), \(\mathbf{R}_{\text{s}}\) corresponds to the previously defined transformation \(\mathbf{R}\), defined for the first panel, second panel and stiffener frames of reference. The internal force vector and tangent stiffness matrix of the cohesive element, rewritten with respect to the new DoF ordering, are denoted as: \(\mathbf{f}^{l}_{\mathrm{c}}\) and \(\mathbf{K}^{l}_{\mathrm{c}}\). Applying the orthogonal transformation \(\mathbf{R}_{\mathrm{c}}\), it is possible to obtain the element force vector and tangent stiffness matrix written for the global coordinate system \[\begin{align} \mathbf{f}^{g}_{\mathrm{c}} = \mathbf{R}_{\mathrm{c}}^{T} \mathbf{f}^{l}_{\mathrm {c}} \text{,} \quad \quad \mathbf{K}^{g}_{\mathrm{c}} = \mathbf{R}_{\mathrm{c}}^{T}\mathbf{K}^{l}_{\mathrm{c}} \mathbf{R}_{\mathrm{c}}. \label{Global95Cohesive95System} \end{align}\tag{39}\]

4 Results↩︎

The proposed structural cohesive element (CE) is first verified using three classical unidirectional laminate benchmarks. A Mode I skin–stiffener debonding test is then used to assess the performance of the structural CEs relative to 3D solid modeling using standard CEs. Finally, a complex stiffened panel under three-point bending is analyzed to demonstrate the computational efficiency of the proposed approach. The analyses were performed in Jive, an open source C++ FEM library [36]. The CEs followed the Turon model [5], with the improved linearization proposed in [34]. All simulations were performed using a flexible path-following method that adopts both displacement control and energy-based arc-length control schemes [37], [38].

4.1 Model verification with benchmarks↩︎

The shell-to-shell cohesive line model was first verified on three classical unidirectional laminate benchmarks, namely the double cantilever beam (DCB), the end-notched flexure (ENF), and the mixed-mode bending (MMB) tests. Although these benchmarks are commonly used for standard delamination model validation, the setup of the shell-to-shell cohesive line models may also be applied to reproduce these tests, as discussed in Section 4.1.2. The structural CE analyses were conducted using the Free Formulation (FF) membrane element from Section 2, and standard Constant Strain Triangle (CST) membrane. A CST membrane is considered by replacing the \(\mathbf{N}_{u}, \mathbf{N}_{v}\) shape functions in the \(\mathbf{B}_{\mathrm{c}_{d}}\) and \(\mathbf{B}_{\mathrm{c}_{\theta}}\) definitions and suppressing the \(\theta_{z}\) DoF.

4.1.1 Description of the benchmark tests↩︎

The benchmarks studied herein are the same as those from [39], with the adaptations from [23]. The geometric parameters and boundary conditions are shown in Figure 3 and Table 1.

The detailed material properties are shown in Tables 2 and [tbl:IM747855295parameters].

10pt

The mixed mode ratio of the MMB specimen is \(50\%\) mode II (\(G_{II}/G_T = 0.5\)). The methodology proposed by [6] for accurately predicting the propagation of delamination under mixed-mode fracture with CEs was applied to the MMB models.

10pt

4.1.2 Description of the models↩︎

The setup of the shell-to-shell structural CE models used to reproduce the benchmarks is illustrated in Figure 4. The delamination benchmarks can be modeled with shells, adopting a horizontal mid-surface shell for the bottom

On the bottom side of the delamination, two horizontal element layers were used to create a central line of nodes for the insertion of shell-to-shell CEs. On the top side, orthogonally oriented element layers with unit thickness (\(t^{s}=1~\mathrm{mm}\)) were used. Shell-to-shell CEs - shown in red in Fig. 4 - were inserted between the two sides of the delamination, while the initial notch region \(a_0\) remained unmeshed. Any contact interaction between the upper ply (stiffener) and the bottom ply (skin) in the ENF test is avoided by simply constraining the vertical displacement of the top-right node of the stiffener. The analytical solutions for the benchmark problems were based on 2D beam models [41]. Due to the unit width of the bottom ply and the unit thickness of the top ply, the numerical responses of the shell-to-shell models are expected to match the mechanical behavior predicted by the analytical solutions, enabling a fair comparison.

For structural CE analyses using the FF membrane, the top ply of the benchmark tests was modeled with 1, 2, and 4 element layers. In contrast, for structural CE models employing the CST membrane, a larger number of layers was used to assess potential locking issues caused by in-plane bending. Specifically, the top ply of the DCB specimen was discretized using 2, 4, and 8 element layers, whereas the top ply of the ENF and MMB specimens was discretized with 4, 8, and 16 layers. For each modeling approach — structural CEs with FF membranes and structural CEs with CST membranes — different mesh refinements along the crack propagation direction were investigated, for fixed numbers of element layers.

4.1.3 Load-displacement results↩︎

The load–displacement curves obtained with the structural CE models are shown in Figures 5 and 6. These curves are plotted together with the analytical solutions from [41]–[43].

The results in Fig. 5 are used to compare the performance of the structural CE models with CST and FF membranes in terms of the number of element layers used for the top-ply discretization. The results show that structural CE - CST membrane models require significantly more element layers in the top-ply discretization than models based on the FF membrane. For the DCB, ENF, and MMB tests, the CST models required 4, 8, and 8 element layers, respectively, with mesh refinements of 0.5 mm, to avoid deviations from the analytical solution. In contrast, the corresponding FF models converged using only 2 layers and mesh refinements of 0.5 mm, 2.5 mm, and 0.5 mm, respectively.

The results in Fig. 6 focus on comparing the models with respect to mesh refinement in the propagation direction for fixed number of element layers for the top ply: 4 layers for the DCB CST-based models, 8 layers for the ENF and MMB CST-based models, and 2 layers for all the FF-based models. Mesh refinements along propagation direction ranging from 0.5 mm to 2.0 mm, 2.5 mm to 7.5 mm, and 0.5 mm to 2.5 mm were investigated for the DCB, ENF and MMB tests, respectively. For comparison purposes, it is interesting to notice that the cohesive zone length for delamination tests can be estimated from the material properties in Tables 2 and [tbl:IM747855295parameters], and the thickness of the plies, as discussed in [44] for idealized pure mode I and II fracture conditions. Therefore, the cohesize zone length expected for the DCB and ENF tests studied herein are approximately 1.7 mm and 8.0 mm, respectively.

The results confirm that, even with a reduced number of top-ply element layers, most of the FF-based structural CE models outperform the CST-based models, providing better approximations of the peak load, even when coarser meshes are employed, as well as less unstable responses during crack propagation. The MMB test is the only exception where the improved performance of the FF-based models is not directly evident. For this case, an additional FF-based model with eight element layers in the top-ply discretization was analyzed, as shown in Fig. 7 (b).

The comparison between CST- and FF-based models with 8 element layers shows a slight improvement when the FF membrane is employed. Very similar results to those shown in Fig. 6 (d) were obtained for FF-based structural CE models with a single top-ply element layer.

4.2 Mode I skin-stiffener debonding↩︎

The Mode I skin–stiffener debonding test was designed to ensure stable Mode I crack propagation along the interface. The geometry and boundary conditions of the 3D solid and shell models are given in Figure 8 and Table 4.

A pre-crack of length \(a_{0}\) was introduced at the loaded end to avoid instability, and a unit thickness \(t\) was assumed for both components. Two types of models with varying skin geometries are considered for the debonding evaluation: skin models with dimensions of 60 mm \(\times\) 50 mm (short debonding interface) and 120 mm \(\times\) 50 mm (long debonding interface).

The elastic properties of the skin are the same as those presented in Table 2, whereas the elastic properties of the stiffener are presented in Table 5. The interfacial fracture properties are given in Table 6. The interfacial properties follow those adopted in [23], except that the Mode I interfacial strength \(\tau_{\mathrm{I}c}\) was increased from \(30\) MPa to \(60\) MPa to reduce the fracture process zone (FPZ) size.

10pt

10pt

4.2.1 Description of the models↩︎

The 3D solid models employed Hex8 elements for both skin and stiffener, with a single layer of elements through the thickness, and linear 8-node cohesive elements at the interface. The shell models adopted an analogous approach using single-layer shell elements and the shell-to-shell CEs. Different mesh refinements were investigated for the 60 mm \(\times\) 50 mm models, with element sizes of 2.0 mm, 1.0 mm and 0.5 mm. The 3D solid models required a fine mesh of 0.5 mm to achieve accurate and stable cohesive crack propagation, whereas the structural shell models provided accurate peak loads and stable crack growth with a coarse mesh of 2.0 mm. Figure 9 (a) illustrates the most refined 3D solid model, whereas Figure 9 (b) illustrates the least refined structural shell model.

The 120 mm \(\times\) 50 mm skin–stiffener debonding tests were then used to assess the increase in computational cost with increasing cohesive interface length, considering only the coarsest mesh refinements that still ensure stable and accurate responses: 3D solid models with the 0.5 mesh and structural models with the 2.0 mm mesh. This enables a fair comparison of the CPU time required for failure analyses using 3D solid models with standard CEs and structural models with the proposed shell-to-shell CEs.

4.2.2 Load–displacement curves↩︎

The load-displacement curves obtained from the simulations of the 60 mm \(\times\) 50 mm models are shown in Figures 10 and 11. The designations “Solid CE” and “Structural CE” denote the results of the standard CE models and those of the proposed structural CE models, respectively.

The results for meshes of different element sizes are plotted together with a reference solution, obtained using a solid CE model with a fine mesh. Two values were considered for the penalty stiffness parameter \(\alpha_{p}\) in Eq. 36 : \(\alpha_{p}=50\), as suggested in [17], and \(\alpha_{p}=5\). Although not shown herein, higher values of \(\alpha_{p}\), such as \(\alpha_{p}=100\), led to increasingly unstable cohesive responses and larger errors in the predicted peak load.

Fig. 10 shows that the standard CE model requires element sizes below 0.5 mm to achieve accurate results. For meshes of 1.0 mm and 2.0 mm, the peak-load error exceeds 9% and 28%, respectively, and the post-peak response becomes highly unstable. In contrast, Fig. 11 shows that the proposed structural CE model remains in close agreement with the reference solution even with the 2.0 mm mesh, with peak-load errors below 8% and a stable post-peak response.

Minor post-peak oscillations observed for coarser meshes are attributed to the larger integration-point spacing. The post-peak oscillations are significantly reduced when \(\alpha_p=5\), indicating that lower penalty parameters than those recommended in [17] for standard CEs can be effectively used. The load-displacement curves obtained from the simulations of the 120 mm \(\times\) 50 mm models are shown in Figure 12. Despite minor peak-load differences, both models yield comparable responses, albeit at markedly different computational costs.

4.2.3 Cohesive stresses↩︎

The cohesive stresses along the interface central line \(\Gamma_c\) during crack propagation are presented in Figure 13 for the solid CE and structural CE models. The results show that the structural CE model with the coarse mesh, with an element size of 2.0 mm, delivered cohesive traction profiles similar to those provided by the most refined models. On the other hand, the solid CE model with the coarse mesh, with an element size of 2.0 mm, resulted in highly oscillatory traction profiles (see Figs. 13a, 13c and 13e).

It can be observed from Figure 13 that the numerically obtained cohesive zone length is smaller than 0.5 mm. The black vertical lines in Fig. 13 indicate the 2.0 mm spacing between the element boundaries adopted for the coarse mesh models. It is interesting to observe how the high-gradient traction profile in the cohesive region can be reproduced within a single high-order structural CE, which is the main reason for the superior performance of the proposed model on coarse meshes.

These results, together with the analytically obtained cohesive zone lengths for the benchmark tests and their corresponding structural responses in Figure 6, allow us to conclude that the proposed model is able to deliver very satisfactory responses, both in terms of load-displacement curves and cohesive traction profiles, for models with element sizes equal to or greater than the estimated cohesive zone length.

4.2.4 Computational performance↩︎

The proposed structural model is able to reduce the computational time of the debonding simulations considerably by allowing larger elements to be used. The results of the structural model are compared against those of the solid model. The CPU times for the analysis of the 60 mm \(\times\) 50 mm models are reported in Table 7.

A comparison between the 3D solid model with a 0.5 mm mesh and the structural model with a 2.0 mm mesh shows that the structural model reduces the computational time by more than 95%, while still providing accurate predictions of the peak load and the overall response. The CPU time for the analysis of the 120 mm \(\times\) 50 mm models are shown in Table 8.

Once again, the comparison between the 3D solid model and the structural model shows that the structural CE approach reduces the computational time of the debonding analysis by more than 95%.

4.2.5 Study on numerical integration↩︎

This section examines the influence of the number of integration points on the simulation results, using the error in the load–displacement curves relative to the reference solution as the evaluation metric. The analysis is restricted to the structural CE model with a 2.0 mm mesh and \(\alpha_{p}=5\). The corresponding load–displacement curves obtained with different numbers of integration points are shown in Figure 14. The numerical integration of the cohesive element equations, described in Section 3.4, may be performed using different numbers of integration points in both the through-thickness and crack-propagation directions. The influence of the number of integration points in each direction is shown in Figures 14 (a) and 14 (b), respectively.

Figure 14 (a) shows that, for this Mode I–dominated debonding test, increasing the number of integration points through the stiffener thickness has no effect on the predicted response. In contrast, the number of integration points in the crack-propagation direction has a minor influence, as illustrated in Figure 14 (b). Analyses with seven or fewer integration points in the propagation direction did not converge beyond the peak load. As reported in [23], increasing the number of integration points improves accuracy for coarse meshes but has negligible impact once the mesh is sufficiently refined, which appears to be the case for the 2.0 mm mesh considered here.

4.3 End-notched flexure test in a complex stiffened panel↩︎

The end-notched flexure (ENF) complex stiffened panel test is illustrated in Figure 15, with its geometrical parameters listed in Table 9. Figure 15 (a) shows a side view of the ENF-like test setup, while the complex stiffened configuration of the panel is shown in Figures 15 (b) in which the red dots indicate points of applied load.

A pre-crack of length \(a_{0}=50\) mm was introduced, cutting all skin–stiffener interfaces across the specimen width \(b\). The thickness \(t\) of the panel and the stiffeners was assumed to be 2.5 mm. The problem is governed by the out-of-plane bending of the panel and the in-plane bending of the stiffeners. The panel has a thickness-to-span ratio of \(t/2L=0.0167\), for which transverse shear effects are expected to be small. For thicker laminates, where transverse shear becomes significant, the multi-layer modelling approach with compliant cohesive interfaces from [23] can still provide the transverse shear deformation needed in the laminate. The elastic material properties adopted for the skin and the stiffener are those previously presented in Tables [tbl:IM747855295parameters] and 5, respectively. The interfacial fracture properties are also taken from Table [tbl:IM747855295parameters]. These fracture properties are assumed to be the same as those adopted in [23] for the original ENF delamination tests.

4.3.1 Description of the models↩︎

The structural finite element models employed shell elements for both the skin and stiffeners, with the proposed shell-to-shell CEs at the interfaces. Both components were discretized using a single layer of shell elements, and structural CEs were inserted in the uncracked region. Two mesh refinements were considered, with element sizes of 2.0 mm and 1.0 mm. Figure 16 illustrates the coarse and fine mesh models, which contain 8,485 and 36,110 nodes (corresponding to 50,910 and 216,660 DoFs), respectively.

4.3.2 Load–displacement curves↩︎

Figures 17 (a) and 17 (b) illustrate the load-displacement curves for the models with mesh elements of sizes 2.0 mm and 1.0 mm, respectively. The convergence of the structural response is observed by comparing the load-displacement responses of both models. Overall, the response resembles that of the original ENF test. However, several sharp snap-backs are observed, and the corresponding critical points are also illustrated in Fig. 17. The nature of these snapbacks is clarified by the damage evolution illustrated in Figures 18 and 19 for the 2.0 mm and 1.0 mm mesh models, respectively.

Even before the first peak load, significant debonding occurs, and the abrupt load drop immediately afterward is associated with the formation of two new damaged regions (Fig. 18a and Fig. 19a).

Similar mechanisms govern the second and third load peaks, where the corresponding load drops are driven by the emergence of symmetric damaged regions (Figs. 18b - 18c, and Figs. 19b - 19c). The subsequent peaks are followed by the appearance of new asymmetric damaged regions (Figs. 18d and 18e) for the 2.0 mm mesh model, and new symmetric damage regions (Figs. 19d and 19e) for the 1.0 mesh model. Beyond a deflection of 7 mm, the damage pattern returns to an approximately symmetric configuration (Fig. 18f and Fig. 19f).

The asymmetry observed in the numerical solution is initially triggered by asymmetry in the mesh. Minor asymmetry suffices for the extremely brittle failure events to occur sequentially on both sides of the panel instead of concurrently. For the model with the 1.0 mm mesh, the strong asymmetries previously observed in Figs. 18 (d) and 18 (e) did not occur, and the damage evolution was approximately symmetric throughout all load steps.

Additionally, Figure 20 shows a zoomed view of a region of interest, where it is possible to observe the cohesive zone - in which the damage field evolves from 0 to 1 - within the domain of a single element, as well as the non-constant damage values at opposite integration points through the stiffener thickness direction.

A more rigorous quantitative analysis of the stability and path-following sensitivity to improve the comprehension of the snap-back phenomena is carried out. To that end, the load-energy dissipation curves are presented in Figures 21 (a) and 21 (b) for the models with mesh element sizes of 2.0 mm and 1.0 mm, respectively.

The critical load steps previously presented in Fig. 17 are also indicated on those curves in order to clarify the physical nature of the snap-back phenomena. The energy dissipation is computed simply by summing up the work done by the cohesive tractions over the jump vector for all the integration points over all cohesive surfaces. The results presented in Figure 21 show that the snap-backs occurring at the critical points are associated with significant energy dissipation, showing that the snapbacks are not numerical artifacts but related to physical crack growth events.

Figure 22 illustrates the deformed configuration for the model with a 2.0 mm mesh at a load step near the end of the analysis.

4.3.3 Computational performance↩︎

The CPU times for the analysis of the models are shown in Table 10.

The CPU time required for the 1.0 mm model - 46.4 hours - is prohibitive for structural engineering practice. In contrast, the 2.0 mm model provides essentially the same structural response at a feasible computational cost of 2.7 hours.

4.4 Pull-out test in a cross-notched complex stiffened↩︎

In order to explore the robustness of the proposed element, the complex stiffened panel models from the previous section were studied under a different notch and load configuration. Instead of an initial edge rectangular notch and mid-span line load pattern, a central cross-notch within a circular region and a central point pull-out force were considered, with the panel fully constrained in the \(w\) direction. Figure 23 presents the notch and load configuration considered for the analysis, in which the cross-notch is indicated by the red dashed line, and the central pull-out force is indicated by the red dot.

The propagation pattern for this type of notch/load configuration is known to produced radial propagation from the initial circular notch in delamination applications [45]. A similar propagation pattern is also observed for the stiffener panel models studied herein, with cracks propagating in the radial direction from the initial circular region, as presented in the next subsection.

4.4.1 Load–displacement curves and damage propagation↩︎

Figures 24 (a) and 24 (b) illustrate the load-displacement curves for the models with mesh elements of sizes 2.0 mm and 1.0 mm, respectively. Very similar responses are obtained in terms of peak load for both models, with the advantage of the analysis with the coarse mesh model taking a fraction of the time to run compared to the analysis with the fine mesh model. The differences observed in the post-peak response may be attributed to mesh asymmetry effects, as discussed in the previous section. The structural response up to complete failure is highly complex, exhibiting several sharp snap-back events. The corresponding critical points associated with some of these snap-backs are also highlighted in Fig. 24.

Figure 25 presents the radial propagation pattern of interface damage for the mesh model with 1.0 mm elements.

The damage maps in Fig. 25 correspond to the damage states at the critical load steps highlighted in Fig. 24 (a). Additionally, Figure 26 shows a zoomed view of a region of interest, where it is possible to observe one more time the cohesive zone within the domain of a single element.

The damage evolution and load–displacement response of the models demonstrate the performance of the proposed formulation under a more challenging mixed-mode propagation scenario, characterized by nonlinear crack front evolution. Figure 27 illustrates the deformed configuration for the model with a 2.0 mm mesh at an intermediary load step.

5 Conclusions↩︎

A novel structural cohesive element is proposed for the efficient modeling of debonding in composite panels with overmolded stiffeners. Three-node higher-order hybrid/mixed shell elements based on the Kirchhoff hypothesis are used to model the panels and the stiffeners. The in-plane bending response of the element is enhanced by the Free Formulation (FF). The proposed cohesive element connects to T-jointed orthogonal shells. An appropriate definition of the displacement jump, consistent with shell kinematics, enables to compute the jump vector at any point over the cohesive surface from the shell displacement approximations evaluated at the element edges. The higher-order weakly continuous fields adopted for the evaluation of the jump vector enable reliable debonding analyses using relatively coarse meshes. The framework is suitable for analyzing debonding in skin–stiffener structures where non-constant damage through the stiffener thickness is expected.

The model is verified for mode I, mode II, and mixed-mode benchmark problems. The results show that the structural cohesive models based on the FF membrane delivered superior performance compared to those employing the Constant Strain Triangle (CST) membrane in most cases. A panel–stiffener debonding problem designed to promote stable crack propagation is also analyzed using both standard 3D cohesive elements and the proposed structural cohesive element with the FF membrane. The results demonstrate that the proposed models can employ coarser meshes than the standard cohesive models, achieving more than a 95% reduction in CPU time, while maintaining comparable accuracy. Consequently, the proposed element enables efficient analysis of stiffener debonding in laminated panels with complex overmolded stiffener grids. The model can be coupled in future studies with the cohesive element in [23] to simulate both delamination and panel–stiffener debonding in progressive failure analyses. It is also applicable to skin–stiffener debonding in overmolded panels with thickness-dependent interface properties, a feature of significant interest in thermoplastic applications.

Declaration of competing interest↩︎

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability↩︎

Data presented in this article will be available at the 4TU.ResearchData repository through https://doi.org/10.4121/9eaeb151-d8c1-42a7-b060-7fea8085e2d2.

Acknowledgements↩︎

This research was carried out as part of the project ENLIGHTEN (project number N21010 g) in the framework of the Partnership Program of the Materials innovation institute M2i (www.m2i.nl) and the Netherlands Organization for Scientific Research (www.nwo.nl).

6 Shape function approximations and derivatives↩︎

The approximations of the in-plane displacement approximations \(\tilde{u}\), \(\tilde{v}\) and their derivatives are expressed in terms of the membrane DoFs \(\overline{\mathbf{u}}\). The out-of-plane approximation \(\tilde{w}\) and its derivatives are expressed in terms of the plate DoFs \(\overline{\mathbf{w}}\).

6.1 Shape functions of \(\tilde{u}\), \(\tilde{v}\) and their derivatives↩︎

Referring back to Eq. (2 ) and the transformation (5 ), the membrane displacement approximations can be written in terms of the membrane DoFs \(\overline{\mathbf{u}}\) as \[\begin{align} \tilde{\mathbf{u}}= \left\{\begin{array}{l} \tilde{u} \\ \tilde{v} \end{array}\right\}= \left[\begin{array}{ccc} \boldsymbol{\phi}_{r}\mathbf{H}_{\mathrm{r}} + \boldsymbol{\phi}_{c}\mathbf{H}_{\mathrm{c}} + \boldsymbol{\phi}_{h}\mathbf{H}_{\mathrm{h}} \end{array}\right]\bar{\mathbf{u}}= \left[\begin{array}{c} \mathbf{N}_{\mathrm{u}} \\ \mathbf{N}_{\mathrm{v}} \end{array}\right]\bar{\mathbf{u}}, \label{Membrane95displacements95Dofs95II} \end{align}\tag{40}\]

where the (\(1\times9\)) matrices \(\mathbf{N}_{\mathrm{u}}\) and \(\mathbf{N}_{\mathrm{v}}\) are the shape functions of \(\tilde{u}\) and \(\tilde{v}\), respectively. The (\(2\times3\)) matrices \(\boldsymbol{\phi}_{r}\), \(\boldsymbol{\phi}_{c}\) and \(\boldsymbol{\phi}_{h}\) are the rigid body, constant strain and higher-order modes, as given in [28], and the (\(3\times9\)) matrices \(\mathbf{H}_{\mathrm{c}}\), \(\mathbf{H}_{\mathrm{r}}\) and \(\mathbf{H}_{\mathrm{h}}\) are defined in Eq. (5 ).

The derivatives of the approximation (40 ) are \[\begin{align} \tilde{\mathbf{u}}_{\color{blue},x}= \left\{\begin{array}{l} \tilde{u}_{\color{blue},x} \\ \tilde{v}_{\color{blue},x} \end{array}\right\}= \left[\begin{array}{ccc} \boldsymbol{\phi}_{r\color{blue},x}\mathbf{H}_{\mathrm{r}} + \boldsymbol{\phi}_{c\color{blue},x}\mathbf{H}_{\mathrm{c}} + \boldsymbol{\phi}_{h\color{blue},x}\mathbf{H}_{\mathrm{h}} \end{array}\right]\bar{\mathbf{u}}= \left[\begin{array}{c} \mathbf{N}_{\mathrm{u}_{\color{blue},x}} \\ \mathbf{N}_{\mathrm{v}_{\color{blue},x}} \end{array}\right]\bar{\mathbf{u}}, \label{Membrane95displacements95derivatives95x} \end{align}\tag{41}\] \[\begin{align} \tilde{\mathbf{u}}_{\color{blue},y}= \left\{\begin{array}{l} \tilde{u}_{\color{blue},y} \\ \tilde{v}_{\color{blue},y} \end{array}\right\}= \left[\begin{array}{ccc} \boldsymbol{\phi}_{r\color{blue},y}\mathbf{H}_{\mathrm{r}} + \boldsymbol{\phi}_{c\color{blue},y}\mathbf{H}_{\mathrm{c}} + \boldsymbol{\phi}_{h\color{blue},y}\mathbf{H}_{\mathrm{h}} \end{array}\right]\bar{\mathbf{u}}= \left[\begin{array}{c} \mathbf{N}_{\mathrm{u}_{\color{blue},y}} \\ \mathbf{N}_{\mathrm{v}_{\color{blue},y}} \end{array}\right]\bar{\mathbf{u}}, \label{Membrane95displacements95derivatives95y} \end{align}\tag{42}\]

where the derivatives \(\boldsymbol{\phi}_{k\color{blue},x}=\partial\boldsymbol{\phi}_{k}/\partial x\) and \(\boldsymbol{\phi}_{k\color{blue},y}=\partial\boldsymbol{\phi}_{k}/\partial y\), \(k=r,c,h\), are computed from the modes \(\boldsymbol{\phi}_{r}\), \(\boldsymbol{\phi}_{c}\) and \(\boldsymbol{\phi}_{h}\) defined in [28].

6.2 Shape functions of \(\tilde{w}\) and its derivatives↩︎

Starting from the cubic approximation (10 ), the displacement \(\tilde{w}\) and its derivatives \(\tilde{w}_{\color{blue},x} = \partial \tilde{w} / \partial x\) and \(\tilde{w}_{\color{blue},y} = \partial \tilde{w} / \partial y\) were expressed in terms of the plate DoFs \(\overline{\mathbf{w}}\) in [23]. The resulting approximations are given by \[\begin{align} \tag{43} \tilde{w} &= \left[\mathbf{S}^{\mathrm{T}} \mathbf{M}_A^{-1}\left(\mathbf{B}_A-\mathbf{M}_B \mathbf{C}\right) +\mathbf{R}^{\mathrm{T}} \mathbf{C}\right] \overline{\mathbf{w}}=\mathbf{N}_w \overline{\mathbf{w}} \\[4pt] \tag{44} \tilde{w}_{\color{blue},x} &=\left[\mathbf{S}_{\color{blue},x} \mathbf{M}_A^{-1}\left(\mathbf{B}_A-\mathbf{M}_B \mathbf{C}\right) +\mathbf{R}_{\color{blue},x}^{\mathrm{T}} \mathbf{C}\right] \overline{\mathbf{w}}=\mathbf{N}_{w_{\color{blue},x}} \overline{\mathbf{w}} \\[4pt] \tag{45} \tilde{w}_{\color{blue},y} &=\left[\mathbf{S}_{\color{blue},y} \mathbf{M}_A^{-1}\left(\mathbf{B}_A-\mathbf{M}_B \mathbf{C}\right) +\mathbf{R}_{\color{blue},y}^{\mathrm{T}} \mathbf{C}\right] \overline{\mathbf{w}}=\mathbf{N}_{w_{\color{blue},y}} \overline{\mathbf{w}} \end{align}\]

where the \((1 \times 9)\) matrices \(\mathbf{N}_w\), \(\mathbf{N}_{w_{\color{blue},x}}\) and \(\mathbf{N}_{w_{\color{blue},y}}\) are the shape functions of \(\tilde{w}\), \(\tilde{w}_{\color{blue},x}\) and \(\tilde{w}_{\color{blue},y}\), respectively. The matrices \(\mathbf{S}\), \(\mathbf{R}\), \(\mathbf{R}_{\color{blue},x}\) and \(\mathbf{R}_{\color{blue},y}\) are the ones that are not constants. They are defined as \[\begin{align} \nonumber \mathbf{S} &=\left[\begin{array}{lll} 1 & x & y \end{array}\right]^{\mathrm{T}},\mathbf{R}=\left[\begin{array}{lllllll} x^2 & xy & y^2 & x^3 & x^2y & xy^2 & y^3 \end{array}\right]^{\mathrm{T}}, \end{align}\]

\[\begin{align} \nonumber \mathbf{R}_{\color{blue},x} &=\left[\begin{array}{lllllll} 2x & y & 0 & 3x^2 & 2xy & y^2 & 0 \end{array}\right]^{\mathrm{T}}, \\[8pt] \nonumber \mathbf{R}_{\color{blue},y} &=\left[\begin{array}{lllllll} 0 & x & 2y & 0 & x^2 & 2xy & 3y^2 \end{array}\right]^{\mathrm{T}}. \end{align}\]

The constant matrix \(\mathbf{C}\), which depends on material and geometrical properties of plate element, is defined as \[\begin{align} \nonumber \mathbf{C}=(\mathbf{H})^{-1}(\mathbf{B}\mathbf{T}) \end{align}\] where the matrices \(\mathbf{H}\), \(\mathbf{B}\) and \(\mathbf{T}\) are the same as those from Eq. (18 ). The matrices \(\mathbf{M}_A\), \(\mathbf{M}_B\), \(\mathbf{B}_A\), \(\mathbf{S}_{\color{blue},x}\) and \(\mathbf{S}_{\color{blue},y}\) are also constant and given by \[\begin{align} \nonumber \mathbf{M}^{\alpha}_A=\left[\begin{array}{lll} 1 & x_1 & y_1 \\ 1 & x_2 & y_2 \\ 1 & x_3 & y_3 \end{array}\right],\mathbf{M}^{\alpha}_B=\left[\begin{array}{lllllll} 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 \\ 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 \\ x_3^2 & x_3 y_3 & y_3^2 & x_3^3 & x_3^2 y_3 & x_3 y_3^2 & y_3^3 \end{array}\right], \end{align}\] \[\begin{align} \nonumber \mathbf{B}_A=\left[\begin{array}{lllllllll} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \end{array}\right], \end{align}\] \[\begin{align} \nonumber \mathbf{S}_{\color{blue},x}=\left[\begin{array}{lll} 0 & 1 & 0 \end{array}\right]^{\mathrm{T}},\mathbf{S}_{\color{blue},y}=\left[\begin{array}{lll} 0 & 0 & 1 \end{array}\right]^{\mathrm{T}}, \label{Sx95Sy95matrices} \end{align}\tag{46}\] in which \(x_i,y_i\), \(i=1,2,3\), are the nodal coordinates of the element. Notice that \(\mathbf{R}_{\color{blue},x}\), \(\mathbf{R}_{\color{blue},y}\), \(\mathbf{S}_{\color{blue},x}\) and \(\mathbf{S}_{\color{blue},y}\) are just the derivatives of \(\mathbf{R}\) and \(\mathbf{S}\) with respect to \(x\) and \(y\).

References↩︎