Isogeometric cohesive zone model for thin shell delamination analysis based on Kirchhoff-Love shell model

Abstract

We present a cohesive zone model for delamination in thin shells and composite structures. The isogeometric (IGA) thin shell model is based on Kirchhoff-Love theory. Non-Uniform Rational B-Splines (NURBS) are used to discretize the exact mid-surface of the shell geometry exploiting their C¹-continuity property which avoids rotational degrees of freedom. The fracture process zone is modeled by interface elements with a cohesive law. Two numerical examples are presented to test and validate the proposed formulation in predicting the delamination behavior of composite structures.

Keywords

cohesive zone model / IGA / Kirchhoff-Love model / thin shell analysis / delamination

Introduction

Fracture in material is usually related to the reduction of material strength which eventually leads to the growth and coalescence of cracks. To adequately represent this complex phenomenon, considerable efforts have been devoted leading to a great number of different computational strategies [¹]. The extended finite element method (XFEM) [²–⁴], generalized finite element method (GFEM) [⁴,⁵], some certain meshfree methods [⁶–⁹], the embedded finite element method [¹⁰], the strong discontinuity embedded approach [¹¹–¹⁴], and the peridynamics formulation for fracture [¹⁵–¹⁷] are among the most popular techniques to model crack propagation with minimal remeshing.

If the process zone is not significantly small compared to the structure dimensions, conventional linear elastic fracture mechanics cannot be applied. To account for the energy dissipation at post-localization, cohesive zone models (CZMs) are commonly used which relate the cohesive tractions at the crack surface to the jump in the displacement field. A cohesive layer is added as a middle surface between two neighboring layers in order to describe the nonlinear interfacial softening response with an additional constitutive relation [¹⁸], these correspond to interface models [¹⁹–²¹] or material interface debonding [²²,²³]. The methodology is applied basing on finite analysis is known as CZMs. The constitutive equation of the interface surface is the cohesive law which describes the traction-jumping displacement relation of the interface. CZMs have also been applied in composites to model delamination between neighbor layers and these layers can be modeled by inserting interface elements across the discontinuity [²⁴]. The law can be based on polynomial, bi-linear and exponential formulations [²³–²⁷] and continuum damage approach [²⁸–³²] to describe the softening behavior when the separation starts to propagate. An overview of different traction-separation laws can be found in Ref. [²¹].

On the other hand, the composite layers are commonly described by shell theories such as Reissner-Mindlin theory or Kirchhoff-Love (KL) theory for thin shells. The shell formulation in exactly capturing the structural geometric was originally proposed in Refs. [³³,³⁴]. Since then there have been a number of shell formulations which account for fracture analysis. However, most approaches were developed for through-the-thickness cracks if classical shell theory was exploited. For instance, the contributions in Refs. [³⁵–³⁷] propose a (local) partition of unity (PU) enrichment in the context of the XFEM [³⁸,³⁹] or the phantom node method [⁴⁰,⁴¹]. Alternative approaches to PU enrichment include efficient remeshing techniques [⁴²,⁴³], enhanced gradient and phase field approaches for fracture [⁴⁴–⁴⁹]. The challenge of applying KL theory in finite elements comes from the fact that C¹ continuity is difficult to achieve by using Lagrange polynomials. A thin shell formulation based on KL theory which does not require rotational degrees of freedom was proposed in the context of meshfree methods exploiting their higher order continuity [⁵⁰,⁵¹]. The formulations in Ref. [⁵⁰] were developed in statics as well as dynamics and also accounted for fracture through the PU enrichment. These formulations were extended to account for fluid-structure interaction, i.e., fluid flow through the opening of dynamically propagating cracks in Ref. [⁵²].

The higher-order continuous IGA formulations have two major advantages over the meshfree approaches: they require less computational effort and the parametrization of complex geometries seems easier though very efficient and interesting parametrizations of complex geometries have been proposed in the context of meshfree LME approaches [⁵³]. Note also that CAD models are commonly based on surface representations and therefore, the IGA concept of integration geometric design into computational analysis [⁵⁴] seems quite natural. We also would like to mention the IGA shell formulations based on Reissner-Mindlin theory [⁵⁵], the blended shells [⁵⁶]. The isogeometric shell solid shells model was proposed by Refs. [⁵⁷,⁵⁸], and the idea of avoiding rotational degrees of freedom with the higher order continuity of IGA basis functions such as Non-Uniform Rational B-Splines (NURBS) and PHT-splines was exploited also in Refs. [⁵⁹–⁶⁴]. The solid-like shell was applied to model propagating delaminations [²⁸,⁵⁷,⁶⁵] where the finite thickness is interpolated by double-knot insertion to reduce the interlayer continuity. However, when the failure phenomena occur in a very narrow area, the thin shell KL model becomes more efficient to describe the decohesion caused by the strong discontinuity between the thin surfaces. Moreover, applying smooth NURBS for KL formulation allows formulating exactly the thin shell surface which would give a benefit in computational efficiency when requiring fewer DOFs is required compared with solid-like shell model.

In this work, differing from the previous finite thickness shell formulations based on solid-like theory, an isogeometric KL thin shell formulation is proposed to study the delamination of the thin shell structures. The fracture of multilayer structures is analyzed in by using a damage based CZM with isotropic damage crossing the interface. In this work, the kinematics of the shell is derived and the implementation of the cohesive zone formulation in the context of IGA is presented. The performance and the capability of proposed formulation in predicting the delamination in thin shell structures are tested via two numerical examples.

Shell kinematics

Figure 1 shows the discretized mid-surface of the shell in the deformed con-figuration \Omega and reference configurations {{\displaystyle \Omega }}_{\mathrm{0}}, respectively. The shell thickness is assumed to be uniform and constant. The position vectors Q and q of material points are presented in the curvilinear coordinates as

(1)\begin{array}{l}\mathit{Q}\left({\mathit{\xi }}^{\mathrm{1}},{\mathit{\xi }}^{\mathrm{2}},{\mathit{\xi }}^{\mathrm{3}}\right)=\mathit{X}\left({\mathit{\xi }}^{\mathrm{1}},{\mathit{\xi }}^{\mathrm{2}}\right)+{\mathit{\xi }}^{\mathrm{3}}{{\displaystyle \mathit{G}}}_{\mathrm{3}}\left({\mathit{\xi }}^{\mathrm{1}},{\mathit{\xi }}^{\mathrm{2}}\right),\hfill \\ \mathit{q}\left({\mathit{\xi }}^{\mathrm{1}},{\mathit{\xi }}^{\mathrm{2}},{\mathit{\xi }}^{\mathrm{3}}\right)=\mathit{x}\left({\mathit{\xi }}^{\mathrm{1}},{\mathit{\xi }}^{\mathrm{2}}\right)+{\mathit{\xi }}^{\mathrm{3}}{{\displaystyle \mathit{g}}}_{\mathrm{3}}\left({\mathit{\xi }}^{\mathrm{1}},{\mathit{\xi }}^{\mathrm{2}}\right),\hfill \end{array}

where -{\displaystyle \frac{h}{\mathrm{2}}}\le {\mathit{\xi }}^{\mathrm{3}}\le {\displaystyle \frac{h}{\mathrm{2}}} with h is the shell thickness. The mid-surface basic vectors are the covariant vectors computed as partial differentiation of the position vectors with respect to the curvilinear coordinates

(2)\begin{array}{l}{{\displaystyle \mathit{G}}}_{\alpha }={\displaystyle \frac{\partial X}{{\partial \mathit{\xi }}^{\alpha }}}={{\displaystyle X}}_{,\alpha },\hfill \\ {{\displaystyle \mathit{g}}}_{\alpha }={\displaystyle \frac{\partial x}{{\partial \mathit{\xi }}^{\alpha }}}={{\displaystyle x}}_{,\alpha }.\hfill \end{array}

The Greek subscripts take the value 1 and 2. The covariant metric coefficients are written as

(3)\begin{array}{l}{{\displaystyle \mathit{G}}}_{\alpha \beta }={{\displaystyle \mathit{G}}}_{\alpha }·{{\displaystyle \mathit{G}}}_{\beta },\hfill \\ {{\displaystyle \mathit{g}}}_{\alpha \beta }={{\displaystyle \mathit{g}}}_{\alpha }·{{\displaystyle \mathit{g}}}_{\beta }.\hfill \end{array}

The contravariant metric coefficients are computed from the inverse of its covariant counterparts

(4)\begin{array}{l}{G}^{\alpha \beta }={\left({{\displaystyle G}}_{\alpha \beta }\right)}^{-\mathrm{1}},\hfill \\ g^{\alpha \beta }={\left({{\displaystyle g}}_{\alpha \beta }\right)}^{-\mathrm{1}}.\hfill \end{array}

The normal vectors of the middle surface are

(5)\begin{array}{l}{\mathit{G}}^{\mathrm{3}}={{\displaystyle \mathit{G}}}_{\mathrm{3}}={\displaystyle \frac{{{\displaystyle \mathit{G}}}_{\mathrm{1}}\times {{\displaystyle \mathit{G}}}_{\mathrm{2}}}{\Vert {{\displaystyle \mathit{G}}}_{\mathrm{1}}\times {{\displaystyle \mathit{G}}}_{\mathrm{2}}\Vert }},\hfill \\ \\ {\mathit{g}}^{\mathrm{3}}={{\displaystyle \mathit{g}}}_{\mathrm{3}}={\displaystyle \frac{{{\displaystyle \mathit{g}}}_{\mathrm{1}}\times {{\displaystyle \mathit{g}}}_{\mathrm{2}}}{\Vert {{\displaystyle \mathit{g}}}_{\mathrm{1}}\times {{\displaystyle \mathit{g}}}_{\mathrm{2}}\Vert }}.\hfill \end{array}

The Green-Lagrange strain tensor is defined as

(6)\mathit{E}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\mathit{F}}^{\mathrm{T}}·\mathit{F}-I\right)={{\displaystyle E}}_{\alpha \beta }{\mathit{G}}^{\alpha }\otimes {\mathit{G}}^{\beta }.

Applying Eqs. (2), (3), (4) to Eq. (6) and neglecting the higher order term, the strain tensor of the middle surface can be expressed as

(7)E_{\alpha \beta }={\epsilon }_{\alpha \beta }+{\mathit{\xi }}^3{\mathit{\kappa }}_{\alpha \beta }.

The first component {{\displaystyle \epsilon }}_{\alpha \beta } contributes to the membrane strain meanwhile the second component {\mathit{\kappa }}_{\alpha \beta } contributes to the bending behavior and they have the forms

(8)\begin{array}{l}{{\displaystyle \epsilon }}_{\alpha \beta }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({{\displaystyle \mathit{g}}}_{\alpha }·{{\displaystyle \mathit{g}}}_{\beta }-{{\displaystyle \mathit{G}}}_{\alpha }·{{\displaystyle \mathit{G}}}_{\beta }\right),\hfill \\ \\ {{\displaystyle \mathit{\kappa }}}_{\alpha \beta }={{\displaystyle \mathit{g}}}_{\alpha }·{{\displaystyle \mathit{g}}}_{\mathrm{3},\beta }-{{\displaystyle \mathit{G}}}_{\alpha }·{{\displaystyle \mathit{G}}}_{\mathrm{3},\beta }.\hfill \end{array}

The displacement field of a material point is

(9)\mathit{u}\left({\mathit{\xi }}^{\mathrm{1}},{\mathit{\xi }}^{\mathrm{2}}\right)=\mathit{X}\left({\mathit{\xi }}^{\mathrm{1}},{\mathit{\xi }}^{\mathrm{2}}\right)-\mathit{x}\left({\mathit{\xi }}^{\mathrm{1}},{\mathit{\xi }}^{\mathrm{2}}\right).

Substituting Eq. (9) into Eq. (8), after neglecting the higher-order term, the membrane strain tensor and the bending strain tensor are written as

(10){{\displaystyle \epsilon }}_{\alpha \beta }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({{\displaystyle \mathit{G}}}_{\alpha }·{{\displaystyle \mathit{u}}}_{,\beta }+{{\displaystyle \mathit{G}}}_{\beta }·{{\displaystyle \mathit{u}}}_{,\alpha }\right),

{{\displaystyle \mathit{\kappa }}}_{\alpha \beta }=-{{\displaystyle \mathit{u}}}_{\alpha ,\beta }·{{\displaystyle \mathit{G}}}_{\mathrm{3}}+{\displaystyle \frac{\mathrm{1}}{\sqrt{{{\displaystyle J}}_{\mathrm{1}}}}}\left[{{\displaystyle \mathit{u}}}_{,\mathrm{1}}·\left({{\displaystyle \mathit{G}}}_{\alpha ,\beta }\times {{\displaystyle \mathit{G}}}_{\mathrm{2}}\right)+{{\displaystyle \mathit{u}}}_{,\mathrm{2}}·\left({{\displaystyle \mathit{G}}}_{\mathrm{1}}\times {{\displaystyle \mathit{G}}}_{\alpha ,\beta }\right)\right]

(11)+{\displaystyle \frac{{{\displaystyle \mathit{G}}}_{\mathrm{3}}\times {{\displaystyle \mathit{G}}}_{\alpha ,\beta }}{\sqrt{{{\displaystyle J}}_{\mathrm{1}}}}}\left[{{\displaystyle \mathit{u}}}_{,\mathrm{1}}·\left({{\displaystyle \mathit{G}}}_{\mathrm{2}}\times {{\displaystyle \mathit{G}}}_{\mathrm{3}}\right)+{{\displaystyle \mathit{u}}}_{,\mathrm{2}}·\left({{\displaystyle \mathit{G}}}_{\mathrm{3}}\times {{\displaystyle \mathit{G}}}_{\mathrm{1}}\right)\right],

where \sqrt{{{\displaystyle J}}_{\mathrm{1}}}=\Vert {{\displaystyle \mathit{G}}}_{\mathrm{1}}\times {{\displaystyle \mathit{G}}}_{\mathrm{2}}\Vert is the length of director field in the current configuration. The resultant effective membrane stress {{\displaystyle \overline{n}}}_{\alpha \beta } and bending stress {{\displaystyle \overline{m}}}_{\alpha \beta } are calculated from the membrane strain {{\displaystyle \epsilon }}_{\alpha \beta } and bending strain {{\displaystyle \mathit{\kappa }}}_{\alpha \beta } as

(12)\begin{array}{l}{\overline{n}}_{\alpha \beta }=t{H}^{\alpha \beta \gamma \delta }{{\displaystyle \epsilon }}_{\gamma \delta },\hfill \\ {\overline{m}}_{\alpha \beta }={\displaystyle \frac{t^{\mathrm{3}}}{\mathrm{12}}}{H}^{\alpha \beta \gamma \delta }{{\displaystyle \mathit{\kappa }}}_{\gamma \delta },\hfill \end{array}

where t is the thickness of the shell and the elastic tensor is expressed in the global coordinate as

H^{\alpha \beta \gamma \delta }={\displaystyle \frac{E}{\mathrm{1}-v^{\mathrm{2}}}}[v\left(G^{\alpha \beta }{G}^{\gamma \delta }\right)+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\mathrm{1}-v\right)\left(G^{\alpha \gamma }{G}^{\delta \beta }\right)

(13)+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\mathrm{1}-v\right)\left(G^{\alpha \delta }{G}^{\gamma \beta }\right)].

The stress tensors are written in Voigt notation as

(14)\left[\begin{array}{l}{{\displaystyle \overline{n}}}^{\mathrm{11}}\hfill \\ {{\displaystyle \overline{n}}}^{\mathrm{22}}\hfill \\ {{\displaystyle \overline{n}}}^{\mathrm{12}}\hfill \end{array}\right]=t\mathit{H}\left[\begin{array}{l}{\epsilon }^{\mathrm{11}}\hfill \\ {\epsilon }^{\mathrm{22}}\hfill \\ {\epsilon }^{\mathrm{12}}\hfill \end{array}\right]\mathrm{,}\left[\begin{array}{l}{{\displaystyle \overline{m}}}^{\mathrm{11}}\hfill \\ {{\displaystyle \overline{m}}}^{\mathrm{22}}\hfill \\ {{\displaystyle \overline{m}}}^{\mathrm{12}}\hfill \end{array}\right]={\displaystyle \frac{t^{\mathrm{3}}}{\mathrm{12}}}\mathit{H}\left[\begin{array}{l}{\mathit{\kappa }}^{\mathrm{11}}\hfill \\ {\mathit{\kappa }}^{\mathrm{22}}\hfill \\ {\mathit{\kappa }}^{\mathrm{12}}\hfill \end{array}\right],

where H is the constitutive tensor which is computed from Eq. (13) for isotropic material as

(15)\mathit{H}={\displaystyle \frac{E}{\mathrm{1}-v^{\mathrm{2}}}}\left[\begin{array}{ccc}\hfill {\left(G^{\mathrm{11}}\right)}^{\mathrm{2}}\hfill & \hfill v{G}^{\mathrm{11}}{G}^{\mathrm{22}}+\left(\mathrm{1}-v\right){\left(G^{\mathrm{12}}\right)}^{\mathrm{2}}\hfill & \hfill G^{\mathrm{11}}{G}^{\mathrm{12}}\hfill \\ \hfill \mathrm{...}\hfill & \hfill {\left(G^{\mathrm{22}}\right)}^{\mathrm{2}}\hfill & \hfill G^{\mathrm{22}}{G}^{\mathrm{12}}\hfill \\ \hfill \mathrm{sym}\hfill & \hfill \mathrm{...}\hfill & \hfill {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\mathrm{1}-v\right)G^{\mathrm{11}}{G}^{\mathrm{12}}+\left(\mathrm{1}-v\right){\left(G^{\mathrm{12}}\right)}^{\mathrm{2}}\hfill \end{array}\right].

CZM for thin shell

In this part, the cohesive formulation for two shell layers with an adhesive interface connecting the separate regions (see Fig. 2) is presented. The formulation is suitable to be extended to the general case of multilayer shells. The boundary conditions are presented in Fig. 2 are Dirichlet boundary condition {{\displaystyle \mathrm{\Gamma }}}_{\mathrm{u}}, Neumann boundary condition boundary condition {{\displaystyle \mathrm{\Gamma }}}_{\mathrm{t}} and the adhesive interface between two shell layers in the reference configuration is denoted by {{\displaystyle \mathrm{\Gamma }}}_{\mathrm{d}}. With an assumption that the adhesive layer thickness is considerably smaller than the in-plan shell dimensions, this idealization allows using a pseudocohesive traction as a constitutive relation to describe the failure of the confined area [³¹,³²]. The kinematic of the interface is possible to be defined by a jumping displacement which enables stress transition from one to another surface. The principle idea of the multilayer shell model is illustrated in Fig. 3, simplified shown in 2D. The interface element inherits the control points and weights of its corresponding layers. The shape function of the interface element is expanded from the NURBS basis function which has the opposite sign for the layer below. The difference between the displacements on two sides of the interface defines the crack opening displacement W. The equilibrium is derived from a weak form of the principle of virtual work which is expressed as

(16)\delta W={{\displaystyle \delta W}}_{{{\displaystyle \mathrm{int}}}_{\mathrm{shell}}}+{{\displaystyle \delta W}}_{{{\displaystyle \mathrm{int}}}_{\mathrm{coh}}}+{{\displaystyle \delta W}}_{\mathrm{ext}}=\mathrm{0},

where {{\displaystyle \delta W}}_{{{\displaystyle \mathrm{int}}}_{\mathrm{shell}}} is the elastic potential of the shell, {{\displaystyle W}}_{{{\displaystyle \mathrm{int}}}_{\mathrm{coh}}} is the elastic potential energy of the cohesive layer and {{\displaystyle W}}_{\mathrm{ext}} is the potential energy of the external force. In the absence of body forces, after applying the Gauss’ theorem, the weak form of the equilibrium equation is rewritten in matrix-vector notation as

(17){{\displaystyle \int }}_{{\mathrm{\Omega }}_{\mathrm{0}}}\left({\mathit{n}}^{\mathrm{T}}\delta \epsilon +{\mathit{m}}^{\mathrm{T}}\delta \mathit{\kappa }\right)\mathrm{d}{\mathrm{\Omega }}_{\mathrm{0}}+{{\displaystyle \int }}_{{\mathrm{\Omega }}_{\mathrm{d}}}\left({\mathit{t}}^{\mathrm{T}}\delta w\right)\mathrm{d}{\mathrm{\Gamma }}_{\mathrm{d}}={{\displaystyle \int }}_{{{\displaystyle \mathrm{\Gamma }}}_{\mathrm{t}}}\left({\mathit{q}}^{\mathrm{T}}\delta u\right)\mathrm{d}{\mathbf{\Gamma }}_t,

where q is the distributed force, t is the traction force on the cohesive surface, w is the jump in the displacement field and n, m, \epsilon, \mathit{\kappa } were defined from Eqs. (10), (11), and (14). The traction force t in Eq. (17) is computed via a cohesive law given by

(18)\mathit{t}=\mathit{D}·\mathit{w},

where D is the elastic tensor of the cohesive zone which takes into account the inuence of damage. The detail of the cohesive relation employed in this study is presented in Appendix A.

Isogeometric element formulation

The total displacement field and the test functions are discretized as

(19)\begin{array}{l}\mathit{u}\left({\mathit{\xi }}^{\mathrm{1}},{\mathit{\xi }}^{\mathrm{2}}\right)={\displaystyle \sum }_{I=\mathrm{1}}^{NP}{{\displaystyle N}}_I\left({\mathit{\xi }}^{\mathrm{1}},{\mathit{\xi }}^{\mathrm{2}}\right){{\displaystyle \mathit{a}}}_I,\hfill \\ \delta \mathit{u}\left({\mathit{\xi }}^{\mathrm{1}},{\mathit{\xi }}^{\mathrm{2}}\right)={\displaystyle \sum }_{I=\mathrm{1}}^{NP}{{\displaystyle N}}_I\left({\mathit{\xi }}^{\mathrm{1}},{\mathit{\xi }}^{\mathrm{2}}\right){{\displaystyle \delta \mathit{a}}}_I,\hfill \end{array}

where {{\displaystyle N}}_I\left({\mathit{\xi }}^{\mathrm{1}},{\mathit{\xi }}^{\mathrm{2}}\right) are NURBS functions defined in the parametric coordinates, NP is the number of control points and a_I is the displacement degree of freedom. Applying Eq. (19) to Eqs. (10) and (11) gives the approximation of membrane and bending strains

(20a)\mathbf{}\mathit{ϵ}={\displaystyle \sum }_I^{NP}{{\displaystyle \mathit{B}}}_m^I\left({\mathit{\xi }}^{\mathrm{1}},{\mathit{\xi }}^{\mathrm{2}}\right){{\displaystyle \mathit{a}}}_I,\delta \mathit{ϵ}={\displaystyle \sum }_I^{NP}{{\displaystyle \mathit{B}}}_m^I\left({\mathit{\xi }}^{\mathrm{1}},{\mathit{\xi }}^{\mathrm{2}}\right){{\displaystyle \delta \mathit{a}}}_I,

(20b)\mathit{\kappa }={\displaystyle \sum }_I^{NP}{{\displaystyle \mathit{B}}}_b^I\left({\mathit{\xi }}^{\mathrm{1}},{\mathit{\xi }}^{\mathrm{2}}\right){{\displaystyle \mathit{a}}}_I,\delta \mathit{\kappa }={\displaystyle \sum }_I^{NP}{{\displaystyle \mathit{B}}}_b^I\left({\mathit{\xi }}^{\mathrm{1}},{\mathit{\xi }}^{\mathrm{2}}\right){{\displaystyle \delta \mathit{a}}}_I.

The matrices B^m and B^b can be computed from Eqs. (B-1) and (B-2) in Appendix B. The jump of the displacement field and the associated test functions are computed as

(21a)\mathit{w}={\displaystyle \sum }_{I=\mathrm{1}}^{NI}{{\displaystyle \mathit{B}}}_c^I\left({\mathit{\xi }}^{\mathrm{1}},{\mathit{\xi }}^{\mathrm{2}}\right){{\displaystyle \mathit{a}}}_I^{\mathrm{int}},

(21b)\delta \mathit{w}={\displaystyle \sum }_{I=\mathrm{1}}^{NI}{{\displaystyle \mathit{B}}}_c^I\left({\mathit{\xi }}^{\mathrm{1}},{\mathit{\xi }}^{\mathrm{2}}\right){{\displaystyle \delta \mathit{a}}}_I^{\mathrm{int}},

where {{\displaystyle a}}_I^{\mathrm{i}\mathrm{n}\mathrm{t}} is the displacement degree of freedom at the interface surface and {{\displaystyle \delta a}}_I^{\mathrm{i}\mathrm{n}\mathrm{t}} is its virtual counterpart. The shape function of traction force is evaluated as

(22){{\displaystyle \mathit{B}}}_c^I=\left[\begin{array}{cccc}\hfill -{{\displaystyle N}}_I\left({\mathit{\xi }}^{\mathrm{1}},{\mathit{\xi }}^{\mathrm{2}}\right)\hfill & \hfill \mathrm{0}\hfill & \hfill -{{\displaystyle N}}_I\left({\mathit{\xi }}^{\mathrm{1}},{\mathit{\xi }}^{\mathrm{2}}\right)\hfill & \hfill \mathrm{0}\hfill \\ \hfill \mathrm{0}\hfill & \hfill -{{\displaystyle N}}_I\left({\mathit{\xi }}^{\mathrm{1}},{\mathit{\xi }}^{\mathrm{2}}\right)\hfill & \hfill \mathrm{0}\hfill & \hfill -{{\displaystyle N}}_I\left({\mathit{\xi }}^{\mathrm{1}},{\mathit{\xi }}^{\mathrm{2}}\right)\hfill \end{array}\right],

where N_I is the number of interface control points. Substituting Eqs. (14), (18), (20a), (20b), (21a), and (21b) into Eq. (17), with the absence of body force b, the obtained equation must be hold for any admissible virtual {{\displaystyle \delta a}}_{\mathrm{1}}, this leads to

(23)\begin{array}{l}{\displaystyle \sum }_{I=\mathrm{1}}^{NP}{\displaystyle \sum }_{J=\mathrm{1}}^{NP}{{\displaystyle \int }}_{{{\displaystyle {\Omega }_{\mathrm{0}}}}_e}\left[t{({{\displaystyle \mathit{B}}}_m^I)}^{\mathrm{T}}\mathit{H}({{\displaystyle \mathit{B}}}_m^J)+{\displaystyle \frac{t^{\mathrm{3}}}{\mathrm{12}}}{({{\displaystyle \mathit{B}}}_b^I)}^T\mathit{H}({{\displaystyle \mathit{B}}}_b^J)\right]\mathrm{d}{{\displaystyle {\mathrm{\Omega }}_{\mathrm{0}}}}_{e}a\\ +{\displaystyle \sum }_{I=\mathrm{1}}^{NI}{\displaystyle \sum }_{J=\mathrm{1}}^{NI}{{\displaystyle \int }}_{{{\displaystyle {\mathrm{\Gamma }}_d}}_e}{({{\displaystyle \mathit{B}}}_c^I)}^{\mathrm{T}}({\mathit{R}}^{\mathrm{T}})\mathit{D}(R){{\displaystyle \mathit{B}}}_c^J\mathrm{d}{{\displaystyle {\mathrm{\Gamma }}_d}}_e{a}^{\mathrm{int}}={\displaystyle \sum }_{I=\mathrm{1}}^{NP}{{\displaystyle \int }}_{{{\displaystyle {\mathrm{\Gamma }}_t}}_e}{{\displaystyle \overline{\mathit{t}}}}^{\mathrm{T}}{{\displaystyle \mathit{N}}}_I\mathrm{d}{{\displaystyle {\mathrm{\Gamma }}_t}}_e,\end{array}

where R is the rotation matrix from the local coordinates to global coordinate. The calculation of this matrix for the mid-surface requires the triad of the interface surface, this is presented detailed in Appendix C. Equation (23) demonstrates a nonlinear behavior, which needs a linearization to solve, hence the Newton-Raphson method is applied

(24)\begin{array}{l}{\displaystyle \mathrm{\sum }}_{I=\mathrm{1}}^{NP}{\displaystyle \mathrm{\sum }}_{J=\mathrm{1}}^{NP}{{\displaystyle \int }}_{{{\displaystyle {\mathrm{\Omega }}_{\mathrm{0}}}}_e}\left[t{({{\displaystyle \mathit{B}}}_I^m)}^{\mathrm{T}}\mathit{H}({{\displaystyle \mathit{B}}}_J^m)+{\displaystyle \frac{t^{\mathrm{3}}}{\mathrm{12}}}{({{\displaystyle \mathit{B}}}_I^b)}^{\mathrm{T}}\mathit{H}({{\displaystyle \mathit{B}}}_J^b)\right]\mathrm{d}{{\displaystyle {\mathrm{\Omega }}_{\mathrm{0}}}}_{e}a\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle \mathrm{\sum }}_{I=\mathrm{1}}^{NI}{\displaystyle \mathrm{\sum }}_{J=\mathrm{1}}^{NI}\left[{{\displaystyle \int }}_{{{\displaystyle {\mathrm{\Gamma }}_d}}_e}t{({{\displaystyle \mathit{B}}}_c^I)}^{\mathrm{T}}({\mathit{R}}^{\mathrm{T}})+(\mathit{D}+{{\displaystyle \mathit{D}}}_q)(\mathit{R}){{\displaystyle \mathit{B}}}_c^J\right]\mathrm{d}{{\displaystyle {\mathrm{\Omega }}_d}}_e\Delta a^{\mathrm{i}\mathrm{n}\mathrm{t}}\\ ={\displaystyle \mathrm{\sum }}_{J=\mathrm{1}}^{NP}{{\displaystyle \int }}_{{{\displaystyle {\mathrm{\Gamma }}_t}}_e}{{\displaystyle \overline{t}}}^{\mathrm{T}}{{\displaystyle N}}_I\mathrm{d}{{\displaystyle {\mathrm{\Gamma }}_t}}_e-{\displaystyle \mathrm{\sum }}_{I=\mathrm{1}}^{NP}{\displaystyle \mathrm{\sum }}_{J=\mathrm{1}}^{NP}{{\displaystyle \int }}_{{{\displaystyle {\mathrm{\Omega }}_{\mathrm{0}}}}_e}\left[t{({{\displaystyle \mathit{B}}}_I^m)}^{\mathrm{T}}\mathit{H}({{\displaystyle \mathit{B}}}_J^m)+{\displaystyle \frac{t^{\mathrm{3}}}{\mathrm{12}}}{({{\displaystyle \mathit{B}}}_I^b)}^{\mathrm{T}}\mathit{H}({{\displaystyle \mathit{B}}}_J^b)\right]\mathrm{d}{{\displaystyle {\mathrm{\Omega }}_{\mathrm{0}}}}_{e}a\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-{\displaystyle \mathrm{\sum }}_{I=\mathrm{1}}^{NI}{\displaystyle \mathrm{\sum }}_{J=\mathrm{1}}^{NI}\left[{{\displaystyle \int }}_{{{\displaystyle {\mathrm{\Gamma }}_t}}_e}{({{\displaystyle \mathit{B}}}_c^I)}^{\mathrm{T}}({\mathit{R}}^{\mathrm{T}})(\mathit{D})(\mathit{R})({{\displaystyle \mathit{B}}}_J^b)\right]\mathrm{d}{{\displaystyle {\mathrm{\Gamma }}_t}}_e{a}^{\mathrm{i}\mathrm{n}\mathrm{t}}\end{array}

Matrix D and D_q comes from the linearization of the traction force (see Appendix A for a detailed calculation). Tracing a nonlinear equilibrium path requires an arc-length technique. In this work, the dissipation-based arc-length control [⁶⁶] is employed.

Numerical integration and stress oscillation phenomena

The full Gaussian integration scheme is applied for all the numerical examples in this work, it means a \left({{\displaystyle p}}_{\mathrm{1}}+\mathrm{1}\right)\times \left({{\displaystyle p}}_{\mathrm{2}}+\mathrm{1}\right) Gauss quadrature rule is employed for calculating the numerical integral of the mid-surface (where p₁ and p₂ are the order of shape functions). From Refs. [¹⁹,³⁰], an undesired spurious oscillation in the stress pattern was observed. A similar phenomena was also report in Ref. [⁶⁷]. This problem comes from the incompatibility between the stiffness of the structure and the interface [¹⁹]. A remedy to overcome this drawback is using a numerical integration basing on Newton-Cotes or Gauss-Lobatto scheme instead of Gaussian quadrature. It is reported in Ref. [⁶⁸] that using higher order continuous isotropic interface elements seems to aggravate the oscillations, and high values of the dummy stiffness is not the only one reason of the oscillatory response. However, in this work, the integrals for the interface and mid-surface KL shell are both evaluated in 2D parametric coordinates, this means there is no problem of dimension inconsistency. We employ standard Gaussian quadrature for both the interface and shell model. In our numerical examples, we show that the spurious oscillations although cannot be totally dismissed but can be reduced significantly when using appropriate fine mesh sizes.

Numerical examples

Double cantilever beam (DCB) test

A panel subjected with initial delamination under tension is considered first. The geometry and boundary conditions of the specimen are shown in Fig. 4. The thickness h of each shell is 0.5 mm with a length L = 10 mm and a width of b = 0.5 mm. An initial delamination length a₀ = 1 mm between two layers is set. The material parameters are: Young’s modulus E = 0.1 GPa, Poisson’s ratio v = 0.3, dummy stiffness K_n = 5 × 10¹⁰ N/m³. Due to the mode I domination, the input parameters for cohesive law are chosen as: the ultimate strength F_sh = F_t = 10⁶ N/m², fracture energy G_I = G_II = 100 N/m, the Benzzegahn-Kenane is 2.284. To simulate the delamination of the DCB, two sets of control points which share a same NUBRBS functions and weight factors are generated. The distance between the upper and lower layers are equal to the shell thickness h. Three discretizations are used, i.e., 13 × 23, 13 × 33, and 13 × 43 control points with quadratic NURBS functions. The loading displacement curves measured at point A for each shell layer is presented in Fig. 5 where the mesh convergence is confirmed. The simulation results have a good agreement with the results from the American Society for Testing and Materials (ASTM) standard [⁶⁹] as well as the reference result for 2D cohesive surface formulation with stress triaxiality [⁷⁰]. Moreover, spurious oscillations can be observe in the region u_z = [⁴,⁵] mm of the loading deection curve from mesh 1. However, this artifact decreases with finer mesh size and becomes insignificant for mesh 3. Figure 6 shows the deformation of the DCB at three loading states S₁, S₂, and S₃ (in Fig. 5) and the component ε_xx of the strain field. As expected, the maximum value of ε_xx evolves along the shell where the delamination propagates.

Mixed mode test

In this second numerical example, the mixed delamination in a curved panel structure is investigated. Similar to the first numerical example, two sets of control points are generated for two mid-surface shells. The inner layer has two radii R₁ = 50 mm and r = 10 mm as depicted in Fig. 7. The outer layer has the same geometry with radius R₂ = R₁ + h where h = 1.5 mm is the shell thickness (see Appendix D for the data of control points and its weigh factors). The potential of surface separation is triggered by an initial delamination with \alpha ={\displaystyle \frac{\mathrm{\pi }}{\mathrm{6}}}. A distributed force in x-direction is applied on the inner layer. The material parameters are: Young’s modulus E = 0.1 GPa, Poisson’s ratio v=\mathrm{0.3}, and dummy stiffness {{\displaystyle K}}_{\mathrm{n}}=\mathrm{3.333}\times {\mathrm{10}}^{\mathrm{11}}\mathrm{N}/{\mathrm{m}}^{\mathrm{3}}. The input parameters for the cohesive law are: ultimate strength {{\displaystyle F}}_{\mathrm{sh}}=\mathrm{5}{{\displaystyle F}}_{\mathrm{t}}={\mathrm{10}}^{\mathrm{7}}\mathrm{N}/{\mathrm{m}}^{\mathrm{2}}, fracture energy {{\displaystyle G}}_{\mathrm{I}}={{\displaystyle G}}_{\mathrm{II}}=\mathrm{400}\mathrm{N}/{\mathrm{m}}^{\mathrm{3}}; the Benzzegahn-Kenane is chosen 2.284. Three meshes are used with 18 × 18, 28 × 18, and 48 × 18 control points with quadratic NURBS functions (see Fig. 8 for the original and refined control mesh for the finest discretization). The load displacement curves for the different meshes show convergent results as illustrated in Fig. 9 which confirm the mesh insensitivity response. Figure 10 illustrates the strain field component ε_xx of the deformed structure (deformations are amplified by 20 times) at several loading steps.

Conclusions

An IGA-based CZM has been developed. Due to the fact that imperfections in the geometry have a very high sensitivity to the performance of the shell especially in thin structures, the proposed IGA formulation and its better geometry representation are beneficial in analyzing the surface separation in composites.

The performance and accuracy of the proposed formulation have been demonstrated through two numerical examples. In the first example of a benchmark (DCB) test, the obtained results show an excellent agreement to the semianalytical results provided in the ASTM standard and a conventional 2D cohesive beam model. The second more complex example shows the ability of the formulation to reliably and mesh independently predict mixed mode fracture of a curved panel with complex geometries.

References

[1] Rabczuk T. Computational methods for fracture in brittle and quasibrittle solids: State-of-the-art review and future perspectives. ISRN Applied Mathematics, 2013, 2013: 849231

[2] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. International Journal for Numerical Methods in Engineering, 1999, 45(5): 601–620

[3] Moës N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. International Journal for Numerical Methods in Engineering, 1999, 46(1): 131–150

[4] Belytschko T, Gracie R, Ventura G. A review of extended/generalized finite element methods for material modeling. Modelling and Simulation in Materials Science and Engineering, 2009, 17(4): 043001

[5] Strouboulis T, Copps K, Babuška I. The generalized finite element method. Computer Methods in Applied Mechanics and Engineering, 2001, 190(32–33): 4081–4193

[6] Belytschko T, Lu Y Y, Gu L. Element-free galerkin methods. International Journal for Numerical Methods in Engineering, 1994, 37(2): 229–256

[7] Zhuang X, Augarde C, Mathisen K. Fracture modeling using meshless methods and level sets in 3d: Framework and modeling. International Journal for Numerical Methods in Engineering, 2012, 92(11): 969–998

[8] Zhuang X, Zhu H, Augarde C. An improved meshless shepard and least squares method possessing the delta property and requiring no singular weight function. Computational Mechanics, 2014, 53(2): 343–357

[9] Zhuang X, Augarde C, Bordas S. Accurate fracture modelling using meshless methods, the visibility criterion and level sets: Formulation and 2d modelling. International Journal for Numerical Methods in Engineering, 2011, 86(2): 249–268

[10] Simo J C, Oliver J, Armero F. An analysis of strong discontinuities induced by strain-softening in rate-independent inelastic solids. Computational Mechanics, 1993, 12(5): 277–296

[11] Zhang Y, Zhuang X. Cracking elements: A self-propagating strong discontinuity embedded approach for quasi-brittle fracture. Finite Elements in Analysis and Design, 2018, 144: 84–100

[12] Zhang Y, Lackner R, Zeiml M, Mang H A. Strong discontinuity embedded approach with standard SOS formulation: Element formulation, energy-based crack-tracking strategy, and validations. Computer Methods in Applied Mechanics and Engineering, 2015, 287: 335–366

[13] Nikolić M, Do X N, Ibrahimbegovic A, Nikolić Ž. Crack propagation in dynamics by embedded strong discontinuity approach: Enhanced solid versus discrete lattice model. Computer Methods in Applied Mechanics and Engineering, 2018, 340: 480–499

[14] Nikolic M, Ibrahimbegovic A, Miscevic P. Brittle and ductile failure of rocks: Embedded discontinuity approach for representing mode I and mode II failure mechanisms. International Journal for Numerical Methods in Engineering, 2015, 102(8): 1507–1526

[15] Han F, Lubineau G, Azdoud Y, Askari A. A morphing approach to couple state-based peridynamics with classical continuum mechanics. Computer Methods in Applied Mechanics and Engineering, 2016, 301: 336–358

[16] Ren H, Zhuang X, Rabczuk T. Dual-horizon peridynamics: A stable solution to varying horizons. Computer Methods in Applied Mechanics and Engineering, 2017, 318: 762–782

[17] Rabczuk T, Ren H. A peridynamics formulation for quasi-static fracture and contact in rock. Engineering Geology, 2017, 225: 42–48

[18] Hillerborg A, Modéer M, Petersson P E. Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cement and Concrete Research, 1976, 6(6): 773–781

[19] Schellekens J, De Borst R. A non-linear finite element approach for the analysis of mode-I free edge delamination in composites. International Journal of Solids and Structures, 1993, 30(9): 1239–1253

[20] Allix O, Ladeveze P, Corigliano A. Damage analysis of interlaminar fracture specimens. Composite Structures, 1995, 31(1): 61–74

[21] Alfano G, Crisfield M. Finite element interface models for the delamination analysis of laminated composites: Mechanical and computational issues. International Journal for Numerical Methods in Engineering, 2001, 50(7): 1701–1736

[22] Needleman A. A continuum model for void nucleation by inclusion debonding. Journal of Applied Mechanics, 1987, 54(3): 525–531

[23] Tvergaard V. Effect of fibre debonding in a whisker-reinforced metal. Materials Science and Engineering A, 1990, 125(2): 203–213

[24] Harper P W, Hallett S R. Cohesive zone length in numerical simulations of composite delamination. Engineering Fracture Mechanics, 2008, 75(16): 4774–4792

[25] Camanho P P, Davila C G, De Moura M. Numerical simulation of mixed-mode progressive delamination in composite materials. Journal of Composite Materials, 2003, 37(16): 1415–1438

[26] Ortiz M, Pandolfi A. Finite-deformation irreversible cohesive elements for three-dimensional crack-propagation analysis. International Journal for Numerical Methods in Engineering, 1999, 44(9): 1267–1282

[27] Hohberg J M. A note on spurious oscillations in FEM joint elements. Earthquake Engineering & Structural Dynamics, 1990, 19(5): 773–779

[28] Nguyen V P, Kerfriden P, Bordas S P. Two- and three-dimensional isogeometric cohesive elements for composite delamination analysis. Composites. Part B, Engineering, 2014, 60: 193–212

[29] Sprenger W, Gruttmann F, Wagner W. Delamination growth analysis in laminated structures with continuum-based 3d-shell elements and a viscoplastic softening model. Computer Methods in Applied Mechanics and Engineering, 2000, 185(2–4): 123–139

[30] Borst R, Gutiérrez M A, Wells G N, Remmers J J, Askes H. Cohesive-zone models, higher-order continuum theories and reliability methods for computational failure analysis. International Journal for Numerical Methods in Engineering, 2004, 60(1): 289–315

[31] Sacco E, Lebon F. A damage-friction interface model derived from micromechanical approach. International Journal of Solids and Structures, 2012, 49(26): 3666–3680

[32] Freddi F, Sacco E. An interface damage model accounting for in plane effects. International Journal of Solids and Structures, 2014, 51(25–26): 4230–4244

[33] Simo J. On a stress resultant geometrically exact shell model. Part VII: Shell intersections with 56-DOF finite element formulations. Computer Methods in Applied Mechanics and Engineering, 1993, 108(3–4): 319–339

[34] Simo J C, Fox D D, Rifai M S. On a stress resultant geometrically exact shell model. Part III: Computational aspects of the nonlinear theory. Computer Methods in Applied Mechanics and Engineering, 1990, 79(1): 21–70

[35] Dolbow J, Moës N, Belytschko T. Discontinuous enrichment in finite elements with a partition of unity method. Finite Elements in Analysis and Design, 2000, 36(3–4): 235–260

[36] Rabczuk T, Bordas S, Zi G. On three-dimensional modelling of crack growth using partition of unity methods. Computers & Structures, 2010, 88(23–24): 1391–1411

[37] Asadpoure A, Mohammadi S. Developing new enrichment functions for crack simulation in orthotropic media by the extended finite element method. International Journal for Numerical Methods in Engineering, 2007, 69(10): 2150–2172

[38] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. International Journal for Numerical Methods in Engineering, 1999, 45(5): 601–620

[39] Moës N, Belytschko T. Extended finite element method for cohesive crack growth. Engineering Fracture Mechanics, 2002, 69(7): 813–833

[40] Song J H, Areias P M, Belytschko T. A method for dynamic crack and shear band propagation with phantom nodes. International Journal for Numerical Methods in Engineering, 2006, 67(6): 868–893

[41] Meer F P, Sluys L J. A phantom node formulation with mixed mode cohesive law for splitting in laminates. International Journal of Fracture, 2009, 158(2): 107–124

[42] Areias P, Rabczuk T. Finite strain fracture of plates and shells with configurational forces and edge rotations. International Journal for Numerical Methods in Engineering, 2013, 94(12): 1099–1122

[43] Areias P, Rabczuk T, Camanho P. Initially rigid cohesive laws and fracture based on edge rotations. Computational Mechanics, 2013, 52(4): 931–947

[44] Areias P, Rabczuk T, Msekh M. Phase-field analysis of finite-strain plates and shells including element subdivision. Computer Methods in Applied Mechanics and Engineering, 2016, 312: 322–350

[45] Thai T Q, Rabczuk T, Bazilevs Y, Meschke G. A higher-order stress-based gradient-enhanced damage model based on isogeometric analysis. Computer Methods in Applied Mechanics and Engineering, 2016, 304: 584–604

[46] Ambati M, De Lorenzis L. Phase-field modeling of brittle and ductile fracture in shells with isogeometric nurbs-based solid-shell elements. Computer Methods in Applied Mechanics and Engineering, 2016, 312: 351–373

[47] Zhou S, Zhuang X, Zhu H, Rabczuk T. Phase field modelling of crack propagation, branching and coalescence in rocks. Theoretical and Applied Fracture Mechanics, 2018, 96: 174–192

[48] Wu J Y. Robust numerical implementation of non-standard phase-field damage models for failure in solids. Computer Methods in Applied Mechanics and Engineering, 2018, 340: 767–797

[49] Wu J Y, Nguyen V P. A length scale insensitive phase-field damage model for brittle fracture. Journal of the Mechanics and Physics of Solids, 2018, 119: 20–42

[50] Rabczuk T, Areias P, Belytschko T. A meshfree thin shell method for non-linear dynamic fracture. International Journal for Numerical Methods in Engineering, 2007, 72(5): 524–548

[51] Amiri F, Millán D, Shen Y, Rabczuk T, Arroyo M. Phase-field modeling of fracture in linear thin shells. Theoretical and Applied Fracture Mechanics, 2014, 69: 102–109

[52] Rabczuk T, Gracie R, Song J H, Belytschko T. Immersed particle method for fluid-structure interaction. International Journal for Numerical Methods in Engineering, 2010, 81(1): 48–71

[53] Amiri F, Millán D, Arroyo M, Silani M, Rabczuk T. Fourth order phase-field model for local max-ent approximants applied to crack propagation. Computer Methods in Applied Mechanics and Engineering, 2016, 312: 254–275

[54] Hughes T J, Cottrell J A, Bazilevs Y. Isogeometric analysis: Cad, finite elements, nurbs, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering, 2005, 194(39–41): 4135–4195

[55] Benson D, Bazilevs Y, Hsu M C, Hughes T. Isogeometric shell analysis: The Reissner-Mindlin shell. Computer Methods in Applied Mechanics and Engineering, 2010, 199(5–8): 276–289

[56] Benson D, Hartmann S, Bazilevs Y, Hsu M C, Hughes T. Blended isogeometric shells. Computer Methods in Applied Mechanics and Engineering, 2013, 255: 133–146

[57] Hosseini S, Remmers J J, Verhoosel C V, De Borst R. Propagation of delamination in composite materials with isogeometric continuum shell elements. International Journal for Numerical Methods in Engineering, 2015, 102(3–4): 159–179

[58] Bouclier R, Elguedj T, Combescure A. Efficient isogeometric NURBS-based solid-shell elements: Mixed formulation and B-method. Computer Methods in Applied Mechanics and Engineering, 2013, 267: 86–110

[59] Benson D, Bazilevs Y, Hsu M C, Hughes T. A large deformation, rotation-free, isogeometric shell. Computer Methods in Applied Mechanics and Engineering, 2011, 200(13–16): 1367–1378

[60] Kiendl J, Bletzinger K U, Linhard J, Wüchner R. Isogeometric shell analysis with Kirchhoff-Love elements. Computer Methods in Applied Mechanics and Engineering, 2009, 198(49–52): 3902–3914

[61] Nguyen-Thanh N, Valizadeh N, Nguyen M, Nguyen-Xuan H, Zhuang X, Areias P, Zi G, Bazilevs Y, De Lorenzis L, Rabczuk T. An extended isogeometric thin shell analysis based on Kirchhoff-Love theory. Computer Methods in Applied Mechanics and Engineering, 2015, 284: 265–291

[62] Nguyen-Thanh N, Zhou K, Zhuang X, Areias P, Nguyen-Xuan H, Bazilevs Y, Rabczuk T. Isogeometric analysis of large-deformation thin shells using RHT-splines for multiple-patch coupling. Computer Methods in Applied Mechanics and Engineering, 2017, 316: 1157–1178

[63] Nguyen-Thanh N, Kiendl J, Nguyen-Xuan H, Wüchner R, Bletzinger K U, Bazilevs Y, Rabczuk T. Rotation free isogeometric thin shell analysis using PHT-splines. Computer Methods in Applied Mechanics and Engineering, 2011, 200(47–48): 3410–3424

[64] Wang P, Xu J, Deng J, Chen F. Adaptive isogeometric analysis using rational PHT-splines. Computer Aided Design, 2011, 43(11): 1438–1448

[65] Reinoso J, Paggi M, Blázquez A. A nonlinear finite thickness cohesive interface element for modeling delamination in fibre-reinforced composite laminates. Composites. Part B, Engineering, 2017, 109: 116–128

[66] Verhoosel C V, Remmers J J, Gutiérrez M A. A dissipation-based arc-length method for robust simulation of brittle and ductile failure. International Journal for Numerical Methods in Engineering, 2009, 77(9): 1290–1321

[67] Kaliakin V, Li J. Insight into deficiencies associated with commonly used zero-thickness interface elements. Computers and Geotechnics, 1995, 17(2): 225–252

[68] Vignollet J, May S, de Borst R. On the numerical integration of isogeometric interface elements. International Journal for Numerical Methods in Engineering, 2015, 102(11): 1733–1749

[69] ASTM D5528-13. Standard test method for mode I interlaminar fracture toughness of unidirectional fiber-reinforced polymer matrix composites. ASTM International, 2013

[70] Remmers J J C, Borst R, Verhoosel C V, Needleman A. The cohesive band model: A cohesive surface formulation with stress triaxiality. International Journal of Fracture, 2013, 181(2): 177–188

[71] Benzeggagh M L, Kenane M. Measurement of mixed-mode de lamination fracture toughness of unidirectional glass/epoxy composites with mixed-mode bending apparatus. Composites Science and Technology, 1996, 56(4): 439–449

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag GmbH Germany, part of Springer Nature