- Email: [email protected]
Sandwich panels under interfacial debonding mechanisms
Accepted Manuscript Sandwich panels under interfacial debonding mechanisms Marco Francesco Funari, Fabrizio Greco, Paolo Lonetti PII: DOI: Reference:
S0263-8223(18)31357-6 https://doi.org/10.1016/j.compstruct.2018.06.113 COST 9907
To appear in:
Composite Structures
Received Date: Revised Date: Accepted Date:
12 April 2018 21 June 2018 27 June 2018
Please cite this article as: Funari, M.F., Greco, F., Lonetti, P., Sandwich panels under interfacial debonding mechanisms, Composite Structures (2018), doi: https://doi.org/10.1016/j.compstruct.2018.06.113
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Sandwich panels under interfacial debonding mechanisms Marco Francesco Funari, Fabrizio Greco, Paolo Lonetti Department of Civil Engineering, University of Calabria, Via P. Bucci, Cubo39B, 87030, Rende, Cosenza, Italy. ABSTRACT: The paper presents a nonlinear approach to investigate the behavior of composite sandwich structures with transversely compressible core, under static and dynamic loading conditions. The proposed model, formulated in the 2D framework, incorporates moving mesh cohesive modeling, crack initiation and nucleation at core/skin interfaces. Interface elements are used to predict debonding mechanisms, whereas shear deformable beams and two-dimensional plane stress elements identify skin and core behavior, respectively. In this framework, interfacial crack onset, layer kinematic and debonding propagation effects are correctly simulated. The moving mesh technique, combined with a multilayer formulation, ensures a reduction of the computational costs, required to predict crack onset and progressive evolution of debonding phenomena. Cohesive models for sandwich core/skin interfaces are calibrated by means of comparisons with numerical and experimental data with respect mode I and mode II configurations. Moreover, a parametric study to address the influence of the loading rate and sandwich characteristics on both static and dynamic frameworks is proposed.
Keywords: sandwich structures, debonding phenomena, moving mesh technique, dynamic fracture mechanics. Nomenclature
a
initial crack length
b B g kf
length of the internal discontinuities crack growth function
GI
energy release rate mode I
GII
energy release rate mode II
GIC
critical strain energy release rate mode I
GIIC
critical strain energy release rate mode II
g
width of the specimen
k f
crack growth function
GC L Ts
length of the specimen
Tc
thickness of the core
v cs
shear wave speed of the core
2∆
internal discontinuity
∆0
characteristic length parameter of the work of separation
Ω
L
per unit area work of separation
thickness of the skin
process zone on the left of the internal discontinuity
ΩR
process zone on the right of the internal discontinuity
Ωc
core domain
Ωs
skin domain
1
Γ S/C 1.
interface between skin and core
Introduction During the last decades, sandwich systems were object of studies and researches from both
manufacturing and design point of views. In particular, such structures are essentially based on two external high density materials, bonded to a thick core, made from low density material. Their configuration, coupled with an optimized material combination, offers a great variety of lightweight structural systems, which are utilized in aerospace, aircraft, marine industries and, recently, also in civil engineering applications [1]. However, the heterogeneity of the layered systems, in terms of geometry and material properties, may be the cause of macroscopic and microscopic damage phenomena. As a matter of fact, during the production process, the assembly of the structural constituents is an important procedure, in which moisture effects [2, 3], temperature variations and contact with ground support equipment [4, 5] are responsible of pre-manufacturing damage phenomena. Furthermore, during the life-structural service, sandwich structures are subjected to low-velocity impacts. This may cause failure mechanisms, with relevant discontinuities between skin and core [6]. In order to investigate the behavior of sandwich structures, several approaches have been proposed in literature. In particular, the analysis started by means of macroscale approaches, in which analytical and numerical formulations aimed to identify critical failure loads for design purpose were developed [7]. Such procedures were endorsed in equivalent plate/shell layer theories, which utilize, through the thickness, a single displacement approximation [8, 9]. Generalized models were also proposed by means multiscale approaches, working at macro-, meso- or microscales, mainly devoted to predict damage behavior of composite structures [10-12]. In particular, the local behavior, at the skin/core bond lines, requires proper micro-mechanical models to simulate the progressive evolution of the interfacial defects [13]. Moreover, at macroscale, nonlinear material models are utilized to predict shear failure mechanisms, transverse cracking, penetration 2
models and core deformations [14]. Typically, the crack growth is fixed along the interface region because of the material mismatch and low mechanical characteristics of the adhesive layer. However, the proposed model requires further extension to simulate kinking phenomena, by introducing crack growth angle criterion. The analysis of interfacial defects and their evolution on sandwich structures can be broadly classified into Cohesive Zone Models (CZMs) or Finite Fracture Mechanics (FFM). In particular, CZMs introduce interface elements at skin/core debonding line, with effective Traction Separation Law (TSL) constitutive relationships. The interlaminar stresses are described in terms of cohesive strength and relative displacements. Formulations available in literature differ in the definition of the constitutive TSLs, which could be based on non-potential or potential models [15-17]. FE models developed in this context range between continuum and discrete type of elements, in which proper TSLs reproduce the rod internal forces vs nodal separation displacements [18]. CZMs predict directly crack onset and propagation of interfacial delaminations, without introducing preexisting debonding length. In certain FE configurations, CZMs may produce numerical difficulties due to mesh dependence, computing inefficiency, sensitivity to the element aspect ratio [19]. Such computational complexities require advanced analyses, which prevent the use of CZMs for engineering practitioners. Alternatively, FFM allows to identify the crack growth by using Energy Release Rate (ERR) concept, which is determined numerically, e.g. based on the Virtual Crack Closure Technique (VCCT) [20] or the J-integral method [21]. Such analyses assume that an initial crack should be introduced prior of the analysis, in which a negligible size of the fracture process zone and the characteristic length of the structure is considered [22]. Although proper crack initiation criteria, based on energy and stress variables, should be utilized to overcome the requirement of an initial crack length, difficulties regarding the identification of the crack front during the mesh transition or the treatment of the singularity behavior of the stresses are well documented in the literature [23].
3
The above referred papers deal with the analysis of the sandwich structures under static loading conditions. However, an important issue to be achieved is the simulation of low or high speed impact conditions. In dynamics, most of the studies focus on the response of the whole structure, which is described in terms of vibration characteristics. Only few attempts discussed the behavior of sandwich structures, in presence of dynamic debonding phenomena. As a matter of fact, 3D analyses incorporating progressive damage models are carried out, in which low-velocity impact responses are investigated in terms of resistance curves and residual plastic deformations [24]. Free vibrations of sandwich panels are discussed in [25], in which a parametric study on the influence of size, location and number of debonding zones is proposed. Moreover, investigations on the dynamic response of debonded sandwich structures, in conjunction with kinematic constraints, governing contact-impact relationships, are proposed in [26]. Recently, dynamic debonding phenomena on sandwich structures were investigated in [27], in which a general model based on Extended HighOrder Sandwich Panel theory and cohesive theory is implemented to predict interfacial crack propagation under mixed model loading conditions. A common feature of the previous models is the treatment of the process zone ahead of the crack tip. In both cases a large number of computational points in proximity of the process zone is required, in order to capture singularity behavior of the stresses or node release effects during the crack growth. This circumstance is particularly relevant in the framework of sandwich structures because of the material heterogeneity. As a consequence, in order to reduce numerical difficulties in complex geometries during the crack advance, remeshing techniques should be utilized, with severe computational costs in the transition procedures [28]. In the present paper a generalized model to predict dynamic debonding phenomena based on the use of moving mesh method is proposed. In spite of other existing methods, the proposed strategy is able to introduce a low number of computational points in the whole geometry, reducing the complexity of the numerical model. This is achieved by using a moving mesh strategy based on Arbitrary-Lagrange Eulerian (ALE) formulation, which is able to reproduce mesh
4
movements according to the process zone motion, ensuring accuracy in the definition of the fracture variables. The outline of the paper is as follows. Section 2 presents the formulation of the governing equations for the ALE and interface approach, whereas in Section 3 the numerical implementation of the finite element model is reported. Finally, comparisons and parametric results to investigate static and dynamic behavior of debonding phenomena are proposed in Section 4.
2.
Formulation of the model The structure is modelled by means of an assembly of layers, which consists of an internal
core and two external skins. At skin/core bond lines, interlaminar interfaces are introduced to simulate initiation and growth of interfacial defects, which may be preexisting or produced during the load history. This is achieved by the use of interface elements based on moving mesh technique, which ensures an accurate description of the fracture variables and the application of cohesive interlaminar stresses in the process zone. Without loss of generality and for the sake of clarity, a sandwich structure is analyzed, in which fracture phenomena affect the upper skin/core interface only, producing two debonding lengths departing from Right (R) or Left (L) directions. A schematic representation of the model is reported in Fig.1.
2.1 Moving mesh cohesive model Moving interface elements between core and skins are introduced to simulate the onset and evolution of debonding phenomena. In particular, the crack growth is expressed as a function of two coordinate systems, i.e. referential and moving, which parametrize the motion of the process zone from the onset to the crack advance (Fig.2). From a mathematical point of view, a computational (Moving) mesh is implemented as the image of the grid points in the referential (Fixed) configuration. In particular, the prescribed motion is expressed in terms of the following Laplacebased equations developed for Static (S) or Dynamic (D) frameworks: 5
∆X
k d ,XX
∂ 2Φ k ( X ,t ) = =0 ∂X 2
∂ 3Φ k ( X ,t ) k =0 ∆ X d ,XX = ∂t ∂X 2
(S)
(D),
(1)
where k, with k=L, R, represents the index referred to the Left (L) or Right (R) process zones and
∆ X dk = X dk ( t ) − X k corresponds to the mesh displacement function [29]. Eqs. (1) require boundary and initial conditions (in dynamics), which depend from the motion of the process zone lengths. ALE interface elements are incorporated in the sandwich structure where debonding phenomena may occur (Fig.3). However, the position, where the crack growth is made possible, is identified on the basis of a proper crack growth function ( g kf ) , defined in terms of mixed mode ERR ratio [15] as follows:
( )
g kf X 1k
( GIC ,GIIC ) are
where
GI X k 1 = GIC
( )
r
2 GII X k 1 + GIIC
( )
r
2 − 1,
(2)
the total area under the traction separation laws, i.e. critical values, and ∆nc
( GI ,GII )
( )
are the individual energy release rates calculated as GI = ∫ Tn ∆n d ∆n and 0
∆tc
GII = ∫ Tt ( ∆t ) d ∆t with ∆n
and
∆s
representing the normal and tangential interlaminar
0
displacements respectively. For each mode components, the TSLs, namely Tn ( ∆n ) or Tt ( ∆t ) , could be described in terms of bilinear, trapezoidal, nonlinear relationships depending from the material interface characteristics. Without loss of generality, a linear-softening curve is considered in the present study. However, a generalization to other forms of damage laws would not require specific changes to the proposed modeling, which is independent from the assumed TSL. Similarly, a different crack onset criterion from Eq.(2) should be introduced to identify the position, in which debonding phenomena are activated.
6
The prescribed mesh motion introduces a rigid displacement of the process zone, which is identified on the basis of internal lengths for the left and right crack directions, namely Ω L or Ω R . Starting from the position ( X T ) in which onset condition is satisfied, two independent coordinates
(X
L T
, X TR ) are generated, which simulate the possibility for the crack tip to advance along the right
and left directions, respectively. The left and right debonding lengths are identified in terms of the positions
(X
L T
R
)
, X T , in which the fracture function g f assumes negative values, within a fixed
tolerance, as follows: L
Ω L = X TL − X T
R
Ω R = X TR − X T
with
g f ( X TL ) = 0
and
g f X T ≤ 0,
( )
(3)
with
g f ( X TR ) = 0 and
g f X T ≤ 0,
( )
(4)
L
R
Finally, the boundary conditions associated to Eq.(1), are required to simulate the displacements of the process zone for the left and right debonding lengths and at geometry edges for the mesh motion:
( )
(5)
( )
(6)
L
∆ X1L ( L ) = 0, ∆ X1L ( X TL ) = ∆ X1k X T = δ L , R
∆ X1R ( 0) = 0, ∆ X1R ( X TR ) = ∆ X1R X T = δ R ,
where (δ L ,δ R ) are the mesh displacements obtained by solving the following optimality conditions: +
+
+
+
δ L ≥ 0 ,g kf ( X TL ) ≤ 0 with δ L ⋅ g kf ( X TL ) = δ L ⋅ g kf ( X TL ) = 0, δ R ≥ 0 ,g kf ( X TR ) ≤ 0 with δ R ⋅ g kf ( X TR ) = δ R ⋅ g kf ( X TR ) = 0.
(7)
From the numerical point of view, the displacements of the debonding regions are determined on the basis of the values of the fracture function at the extremities of the debonding length. In particular, during the crack growth a null value of the fracture function is enforced at the internal extremity of the debonding length. Therefore, with reference to k-th defect, the crack tip displacement can be expressed by means of the following relationship: 7
∆X
kL ,R T
=0
g
kL ,R f
≤ 0, ∆ X
kL ,R T
=
g kf ( X 1kL ,R ) + toll g kf ( X 1kL ,R ) + toll ± g kf ( X 1kL ,R ± Ω kL ,R )
,R =0 Ω kL ,R g kL f
(8)
where toll is tolerance positive value able to trigger the crack growth. From the numerical point of view it has been chosen close to zero, i.e. 1E-2, to do not modify the interfacial cohesive laws [30]. 2.2 Kinematic formulation of the model
The skins follow a one dimensional model based on Timoshenko beam kinematics, whereas the core is modelled by using a 2D plane stress formulation. The materials of each layer are assumed to be linear elastic and hypotheses concerning small displacements and strains are considered in the analyses. Moving interface elements are introduced between the skin and core to reproduce potential debonding phenomena. The interfaces are initially assumed to be perfectly bonded to both upper and lower skin/core lines and thus moving coordinates match with the referential or material ones [30]. In this framework, dynamic equations of motion can be expressed by introducing equilibrium conditions along the horizontal and vertical directions and flexural rotation (counterclockwise) as follows: − µU2k + EAU 2k,X ,X + FX ( X ) + p x ( X ) = 0 , − µV2k + GA* (V2k,X − Φ 2k ) ,X + FY ( X ) + p y ( X ) = 0 , − I 0Φ2k + EIΦ 2k.XX − GA* (V2k,X − Φ 2k ) + m ( X ) = 0 ,
(
on Ω S , Ω S
div [ E∇U ] + f − µU = 0,
on
E33 (U ,YC + V,XC ) = FX ( X )
on
(E
C C 12U ,X + E22V,Y ) = FY ( X )
1
ΩC
(Γ
S1 / C
2
)
(9)
(10)
.Γ S2 / C
)
(11)
where, as shown in Fig.3, the superscript k (with k=S1,S2), here and in the following, refers to the lower (S1) or upper skins (S2), ( Ω X ,Γ X / Y ) are the domain of the X or the boundary between X and Y (with X,Y=S1, S2, C, ALE1, ALE2), ( I 0 , µ ) are the per unit length inertial rotation and mass of the skin, ( EI ,GA* ) are the flexural and shear stiffnesses, 8
( p , p ,m ) are the axial, vertical and flexural x
y
d(X) external loads applied to the skins, the notation X = indicates the time derivative of the dt
( )
function X, Eij (with i,j=1,2,3) is the 2D plane stress elastic matrix, U (with U T = U C V C ) is the vector containing the horizontal and vertical displacements of the core, f is the vector force per
unit volume. Moreover, the interlaminar stresses are defined on the basis of moving mesh coordinates and interlaminar relative displacements as follows:
FX ( X1k ) = Tn ∆nk ( X1k ) FY ( X1k ) = Tt ∆tk ( X1k )
,
on
Γ S1 / C and
Γ S2 / C
(12)
,
(13)
with
∆nS ( X ) = V C ( X , −0.5 ⋅ T C ) − V S ( X ) , 1
1
∆nS ( X ) = V S ( X ) − V C ( X ,0.5 ⋅ T C ) , 2
2
∆tS ( X 1R ) = U C ( X , −0.5 ⋅ T C ) − U S ( X 1R ) − Φ S ( X 1R )( 0.5 ⋅ T S ) 1
1
1
∆tS ( X 1R ) = U S ( X 1R ) + Φ S ( X 1R )( 0.5T S ) − U C ( X ,0.5 ⋅ T C ) 2
s
2
where T S and T C are the thickness of the skin and core, respectively. Previous equations, i.e. Eq.(9), should be considered in those cases, in which the interfaces are not affected by debonding phenomena. However, once crack initiation condition, defined by Eq.(2), is satisfied, at the position
X = X T debonding phenomena should be simulated from the right and left directions. In particular,
(
)
L R two different positions, i.e. X = X T , XT , which identify fictitious crack tips for the right and left
debonding lengths, are introduced. Such quantities differ from the internal length 2∆ , which refers to the initial defect predicted by the crack growth criterion. From the numerical point of view, they are assumed very small, in order to simulate the presence of an initial defect. For more details, the influence of such parameter was investigated in [30]. Therefore, the governing equations are expressed as follows:
9
− µU2k + EAU 2k,X ,X + FX ( X R ) + FX ( X L ) + px ( X ) = 0 , − µV2k + GA* (V2k,X − Φ 2k ) ,X + FY ( X R ) + FY ( X L ) + p y ( X ) = 0 , − I 0Φ2k + EIΦ 2k.XX − GA* (V2k,X − Φ 2k ) + m ( X ) = 0,
on Ω S 2
(14)
where the superscript k is k=S 2, ( X R ,X L ) are the effective moving coordinates for the left (L) and right (R) directions measured on the material domain. Moreover, the governing equations of the core, i.e. field and initial conditions, refer to a classical 2D plane stress modeling, which, for the sake of brevity, is not reported in detail in the paper. Boundary conditions concerning the interface regions between core and skins are defined by the following expressions:
E33 (U ,Y + V,X ) = FX ( X R ) + FX ( X L )
(E
U ,X + E22V,Y ) = FY ( X R ) + FY ( X L )
on Γ S2 / C
(15)
12
It is worth noting that interface traction forces are expressed in terms of the moving coordinate systems, i.e.
(X
R
, X L ) , for the right and left debonding phenomena, leaving unaltered the
governing equations of the structural elements, with respect to the standard material formulation [31].
3.
Numerical implementation
Governing equations presented in previous section are implemented by using a finite element approach. For the sake of conciseness, the implementation procedure is mainly referred to the ALE interface elements and moving traction forces acting at the skin/core interfaces only, since the governing equations of the structural model remain basically unaltered. In particular, Lagrange interpolation functions with cubic or quadratic orders are assumed for the Timoshenko or plane stress finite elements, whereas Lagrange linear interpolation functions are utilized for the moving interface variables:
10
NeS
q S = (U ,V ,Φ ) = ∑ NiS ( X ) (UiS ,Vi S ,Φ iS ) = N S1,2 U i =1 NeC
S1 ,2
(
on Ω S , ΩS
,
q = (U ,V ) = ∑ NiC ( X ) (U iC ,ViC ) = N C U , i =1 1,2
C
R
NeMM
X =
∑ F ( X ) ( X
R
i
i =1
) = F X
R
C
,X =
∑ F ( X ) ( X i
i =1
) = F X
L
(
1
(
1
, on Ω ALE ,Ω ALE
∑ F ( X ) ( Λ ) = F Λ , i
) (16)
L
NeMM
Λ=
2
on ΩC
NeMM
L
1
on Ω ALE ,Ω ALE
i
i =1
2
)
2
)
where • contains the nodal quantities of the [• ] function, N eS , N eC ,N eMM are the number nodes
(
(
)
)
in the skin, core and ALE elements, respectively, N S1 ,S2 ,N C ,F are the matrixes containing the
(
)
interpolation functions for the structural and ALE elements, q S1,2 ,q C collect nodal variables of the structural problem for the skin and core and
( X
R
, X L ) are the interpolation functions of the right
and left moving debonding length. Finally, Lagrange multiplier method is introduced as a result of weak constraint conditions by function Λ , which is assumed with a linear interpolation. Starting from Eq. (1) and Eq.(14), introducing weighing functions related to translational displacements (WUS,V ,WUC,V ) for the skins and core, rotation (WΦS ) fields for skins and positional variables (W k ) for the ALE variables, the weak forms for interface elements and the moving traction forces on the skin/core interfaces present the following relationships:
(∆X
k ,X
,W,Xk )
Ω ALE1
k C V Γ
+ (Tn ,W − W S V
(
= Λ1k δ ( X Tk ) ,W k
)
( ( ) ) k
Ω ALE1
k
+ Λ1 δ X T ,W k
+ Ω ALE1
k
)
S1 / C
T + Tt ,WUS + WΦS − WUC + 2 Γ S / C
(17)
1
(
+ ∆ X Tk − δ k δ ( X Tk ) ,WΛk 1
)
Ω ALE1
(
( ) )
k k + ∆ X T − δ k δ X T ,WΛ 2
=0 Ω ALE1
with Tnk = C nk ∆n ( X k ) ∆n ( X k ) , Tt k = Ctk ∆t ( X k ) ∆t ( X k ) 11
(18)
where the notation ( ⋅,⋅)Ω denotes ( ⋅,⋅)Ω = ∫ ( ⋅,⋅)dX , δ ( i ) is the delta dirac function, the superscript k Ω
(with k=R,L) refers to the Right and Left interfaces, and ( ∆s , ∆n ) are the interlaminar displacements defined in Eq.(13). Eqs.(18) are rearranged in vector notation as follows: T T = t Tn
Dt ( ∆t = 0
)
∆ t = D∆ Dn ( ∆n ) ∆n 0
(19)
where Dt ( ∆t ) ,Dn ( ∆n ) are the tangential and normal stiffnesses of the introduced to describe nonlinear TSLs. Moreover, the interlaminar displacement vector, ∆ , can be expressed by
(
introducing a mixed interpolation function matrix, namely N SC and two auxiliary matrixes Q1 ,Q2
)
:
T ∆t U S1 + Φ S1 − U C SC = = Q1 N S U S − Q2 N C U C = N SC U 2 ∆n V S1 − V C where U
SC
(with U
SC
C
(20)
S
= U ∪ U ) is the vector containing the nodal displacements of the core and
skin element. Since the interlaminar foces are evaluated on the moving domain, a mapping operator is required to identify the connection between referential and material coordinate systems[32]. To this end, a projector operator Φ : X → X d is introduced, which the relationship between the nodal
coordinates in the moving and referential configurations, i.e X d and X , with X d = Φ X . Substituting Eq,(16) into Eq.(17) and by using Eq.(20), the following discrete equations in which unknown quantities are represented by the position of the computational nodes, are derived: k
Ak X d + B1k λ1k + B2k λ2 + C CS ( k )U SC( k ) = 0 , k k k k k k B1 X d − C1 = 0 , B 2 X d − C 2 = 0 ,
with
12
(21)
X i +1 A = ∫ F,X ⊗ F,X dX ,B1 = δ ( X T ) F T ,B2 = δ ( X T ) F T X T X i
( )
C
CS
Xi +1 = ∫ ( N SC Φ X i
T
) D ( N
SC
Φ
) dX ,
(22)
( )(∆X
C1 = δ ( X T )( ∆ X T − δ f ) , C2 = δ X T
k T
− δ f ) ,
where ⊗ is the dyadic product between vectors and ( ) is the composition function. Eqs.(21), which are not modified from the classical definition [33], are implemented by using a customized FE subroutine in the framework of COMSOL Multiphysics, by means scripting capabilities of MATLAB® language [34]. In particular, the analysis is based on several steps, which are connected by using stop/restart procedures. In particular, the built-in Newton-Raphson or implicit transient analyses based backward differentiation formula (BDF) for static and dynamic analyses are employed. These are managed by means of external script files based on stop condition for crack initiation and restart/remeshing procedures for the activation of moving interface elements. A synoptic representation of the computational algorithm implemented in the FE environmental program is reported in Tab.1.
4.
Validation and numerical results
In this section, the proposed model is verified by means of comparisons with numerical and experimental data. The first step is developed with the purpose to analyze the consistency of the proposed formulation under mode I and mode II loading conditions. In particular, according to standard experimental methods, the static behavior of interfacial cracks for several sandwich debonding configurations is investigated. The main aim of the comparisons is to validate the proposed model and to examine its ability to describe debonding failure mechanisms in sandwich configurations. Subsequently, the structural behavior is investigated in both static and dynamic
13
frameworks to identify the influence of inertial effects on the crack growth phenomena, produced by different levels of the loading rate. 4.1 Mode I DCB loading scheme
At first, analyses are developed with reference to a classical DCB loading involving a pure mode I fracture analysis. The loading, boundary conditions and geometry of the specimen are illustrated in Fig.4, whereas, according to data recovered in [27, 35], mechanical properties assumed for the laminate, core and interfaces as well as the ones required by the cohesive elements, are summarized in Tab.2. The crack evolves in the interface region between core and face sheet (namely at Γ S2 / C ). In particular, a preexisting delaminated length is assumed in the upper interface, whereas the lower interface is considered as perfectly bonded. The numerical model is discretized by means of shear deformable beam elements with cubic interpolation functions for each face sheet, whereas the core is modeled by quadratic interpolation functions with linear elastic material behavior. Moreover, ALE elements with linear interpolation functions are introduced between each sheet-core interfaces to predict crack onset or debonding advance. The analysis is developed under a displacement control mode to ensure a stable crack propagation. Cohesive constitutive relationships are based on exponential laws, with smooth decay branches according to [36]. The numerical discretization utilized for the face sheet is assumed to be generally uniform with a length equal to ∆M / L = 1 / 30 , with ∆M the element length. The core has been discretized by means plane stress quadrilateral elements with maximum element length equal to ∆M / L = 1 / 30 . At the interfaces a coarse discretization is adopted, except at the process zone ( Ω R ,L ) where a mesh enrichment is considered, i.e. ∆M / L = 1 / 300 . In Fig. 5, comparisons in terms of loading curves as a function of opening end displacements, for two different core thickness configurations, i.e. Tc=15,20 mm, are reported. The results show that the proposed approach is able to reproduce correctly the behavior of the structure, since the loading curves are quite in agreement with 14
experimental [35] and analytical data available in [27]. The results show that an increment in the core thickness does not produce improvements in the loading curves, which are basically coincident in terms of stiffness, displacements except for a small variation of the peak load. Previous analyses developed essentially in statics are extended in a dynamic framework. The main aim is to investigate the influence of the loading rate on both resistance curve and crack tip speed evolution. In addition, the effect of an internal preexisting delamination region is also considered to verify the dynamic amplifications produced during the crack growth and their interaction with an external debonding length. The structure is subjected to an applied speed with a ramp curve until the time (t0), which is assumed proportional to the first period of vibration of the structure (t0=0.5T1); subsequently a fixed velocity is prescribed to the end displacement (v0). In Fig.6, comparisons in terms of resistance curve between static and dynamic results are developed, obtained by means several prescribed opening speeds with refers to low or high ranges. Moreover, the following configurations are considered in the analyses: •
external debonding length only with b=0 (C1);
•
internal and external discontinuities with different geometric lengths, i.e. b/a=0.2 (C2) or b/a=0.5 (C3).
The results show that under static loading the presence of an internal debonding length modifies load-displacement curves, increasing system instability during the crack growth and leading large oscillations with respect to the static case. The presence of an internal debonding, at low speeds, does not modify the evolution curve, which, apart little oscillations in the unstable path, is quite coincident with the C1 configuration. However, at large speeds, the dynamic effects modify the crack growth behavior, since C1 and C2 denote notable oscillations, which tend to be reduced, for increasing values of the opening displacements, leading to same predictions of the static solution. In Fig. 7, the influence of an internal discontinuity on the dynamic behavior of debonding phenomena is discussed. In particular, for a fixed loading rate, comparisons in terms of loading 15
curve are proposed for the C1, C2 and C3 configurations. The results show how the presence of an internal discontinuity modifies the equilibrium path, introducing a marked span-back phenomenon in the loading curve. However, once the crack tip overcomes the internal debonding length, the curves for the configurations with internal debonding lengths, i.e. C2 and C3, tend to the same prediction. Additional results are presented in Figs. 8-9, in which the relationship between speed and position of the crack tip is investigated. The analyses show that during the crack onset, the tip is affected by large accelerations, which tend to be reduced as far as the crack tip speed grows. However, at low loading rates, after a maximum value is reached, the crack tends to a steady state advance, since a constant speed is observed. In presence of an internal debonding region, the crack is affected by an oscillating behavior and, especially at low loading rates, crack arrest phenomena are observed. This phenomenon is mainly produced by the presence of the internal debonded region, which produces a local distribution of contact forces with respect to the case of perfectly debonded region. Large amplifications in the crack tip speeds are observed, especially, during the transition regions from the adhesive to the delaminated length.
4.2 Mode II CSB loading scheme
The analyses are extended to a loading scheme involving a mode II loading condition. Comparisons with numerical data are reported for a sandwich configuration based on a three-point bending loading scheme. As shown in Fig.10, the initial crack is located in the upper interface between the face-sheet and the core, whereas geometrical and mechanical characteristics are reported in Tab.3. The numerical discretization utilized for the face sheet is assumed to be generally uniform with a length equal to ∆M / L = 1 / 55 . The core has been discretized by means plane stress quadrilateral elements with maximum element length equal to ∆M / L = 1 / 55 . At the interfaces a coarse discretization is adopted, except at the process zone ( Ω R ,L ) where a mesh enrichment is considered, i.e. ∆M / L = 1 / 550 . The following configurations are considered in the analyses: 16
•
external debonding length, i.e. b=0, (C1);
•
internal and external discontinuities with b/a=0.5 (C2).
Comparisons in terms of applied loads as a function of mid-span vertical displacement for two different preexisting crack length configurations are reported. In both cases, results obtained by the proposed model are in agreement with the ones obtained by analytical approach based on pure cohesive model [27]. Moreover, in Fig. 12, comparisons are developed to verify the consistency of the solution with respect to the local distribution of interfacial stresses. Despite numerical data reported in the literature, the proposed model is able to reproduce the crack growth also for more extended tip displacements. The results point out in the load-displacement curves the presence of a knee once the position of the crack tip overcomes the one of the applied loads. The accuracy of the solutions is guaranteed also locally in terms of stress distribution. Good agreement with the values obtained by using a pure cohesive modeling, in which an accurate description in terms of mesh element length is required along all the interface regions is found. Contrarily in the proposed model, based on ALE approach, the mesh refinement is introduced along the process zone only, leading to a reduction of the computational costs involved in the analysis. Previous results developed essentially in statics are extended in a dynamic framework, in which applied speeds based on a ramp curve with different loading rates are considered. Resistance curves are reported in Figs. 13-14, which refer to configurations with or without an internal debonding length, i.e. C1 and C2, respectively. The results show that, at high loading rates, large oscillations and amplifications with respect to the static solution are observed. Similarly, to results obtained under pure mode I loading scheme, Mode II configuration is also affected by the presence of an internal discontinuity, since an instable behavior in the loading curve with small oscillations is observed with respect to the static equilibrium path. The stress distributions for two different structural scenarios related to C1 and C2 configurations are reported in Fig. 15. The results show that, in presence of an internal debonding length, the traction forces present a discontinuous behavior at the extremity points, leading to singular values of the interfacial stresses. In Fig.16, 17
comparisons in terms of crack tip speeds, denote how the presence of an internal discontinuity produces an amplifications of the tip speed also for low loading rates, leading to jumps of crack speeds once the tip reaches the debonded length. Moreover, as far as the crack tip overcomes the delaminate length; the crack tip speed oscillates leading to crack arrest phenomena during the crack evolution.
5.
Conclusions
The present study was aimed to investigate the behavior of sandwich structures in presence of debonding phenomena. The proposed model combines ALE formulation with a cohesive fracture approach to simulate debonding phenomena in presence of multiple delaminations. Based on the results obtained, the following conclusions can be drawn: •
the proposed methodology is able to reduce computational difficulties arising from the mesh dependence, convergence and accuracy of the solution, which typically affect cohesive approaches;
•
the results show the capabilities of the proposed model to reproduce correctly the behavior in both static and dynamic frameworks;
•
the parametric study denotes how debonding phenomena are quite influenced by the presence of internal material discontinuities, which produce strong oscillations and instabilities in the resistance curve and high speeds during the crack evolution;
•
this study could be very useful for developing modeling strategies to understand debonding phenomena in sandwich composites subject to low-velocity impact;
•
the proposed approach requires a further extension to predict kinking phenomena, since at this stage the crack position is fixed along the interfaces.
Acknowledgments
This work is supported by Italian Ministry of University and Research (P.R.I.N. National Grant 2015, B86J16002300001). 18
References
[1] Manalo A, Aravinthan T, Fam A, Benmokrane B. State-of-the-Art Review on FRP Sandwich Systems for Lightweight Civil Infrastructure. Journal of Composites for Construction. 2017;21(1). [2] Yu H, Zhou C. Sandwich diffusion model for moisture absorption of flax/glass fiber reinforced hybrid composite. Composite Structures. 2018;188:1-6. [3] Spadea S, Orr J, Ivanova K. Bend-strength of novel filament wound shear reinforcement. Composite Structures. 2017;176:244-53. [4] Ehsani M, editor ASCE Innovation Award Winner: Sandwich Construction Carbon FRP Pipe2017: American Society of Civil Engineers (ASCE). [5] Maletta C, Sgambitterra E, Niccoli F. Temperature dependent fracture properties of shape memory alloys: Novel findings and a comprehensive model. Scientific Reports. 2016;6. [6] Yu RC, Pandolfi A. 15 - Modeling of delamination fracture in composites: a review A2 Sridharan, Srinivasan. Delamination Behaviour of Composites: Woodhead Publishing; 2008. p. 429-57. [7] Triantafillou TC, Gibson LJ. Failure Mode Maps for Foam Core Sandwich Beams. Materials Science and Engineering. 1987;95:37-53. [8] Abrate S, Ferrero JF, Navarro P. Cohesive zone models and impact damage predictions for composite structures. Meccanica. 2015;50(10):2587-620. [9] Bruno D, Greco F, Lonetti P. Dynamic Mode I and Mode II Crack Propagation in Fiber Reinforced Composites. Mechanics of Advanced Materials and Structures. 2009;16(6):442-55. [10] Greco F, Leonetti L, Lonetti P, Nevone Blasi P. Crack propagation analysis in composite materials by using moving mesh and multiscale techniques. Computers and Structures. 2015;153:201-16. [11] Ascione L, Berardi VP, Giordano A, Spadea S. Pre-buckling imperfection sensitivity of pultruded FRP profiles. Composites Part B: Engineering. 2015;72:206-12. [12] Ciantia MO, Arroyo M, Butlanska J, Gens A. DEM modelling of cone penetration tests in a double-porosity crushable granular material. Computers and Geotechnics. 2016;73:109-27. [13] Rabinovitch O, Hamed E. Bending behavior of sandwich panels with a "soft" core and embedded rigid inserts. Sandwich Structures: Advancing with Sandwich Structures and Materials. 2005:261-70. [14] Mohammadi MS, Nairn JA. Balsa sandwich composite fracture study: Comparison of laminated to solid balsa core materials and debonding from thick balsa core materials. Composites Part B-Engineering. 2017;122:165-72. [15] Bruno D, Greco F, Lonetti P. Computation of energy release rate and mode separation in delaminated composite plates by using plate and interface variables. Mechanics of Advanced Materials and Structures. 2005;12(4):285-304. [16] Park K, Paulino GH. Cohesive Zone Models: A Critical Review of Traction-Separation Relationships Across Fracture Surfaces. Applied Mechanics Reviews. 2011;64(6). [17] Funari MF, Greco F, Lonetti P, Luciano R, Penna R. An interface approach based on moving mesh and cohesive modeling in Z-pinned composite laminates. Composites Part B-Engineering. 2018;135:207-17. [18] Xu G, Yan RJ. The Use of Sprint Interface Element Delamination Simulation of Sandwich Composite Beam. Applied Composite Materials. 2017;24(5):1049-60. [19] Rabinovitch O. Debonding analysis of fiber-reinforced-polymer strengthened beams: Cohesive zone modeling versus a linear elastic fracture mechanics approach. Engineering Fracture Mechanics. 2008;75(10):2842-59.
19
[20] Okada H, Higashi M, Kikuchi M, Fukui Y, Kumazawa N. Three dimensional virtual crack closure-integral method (VCCM) with skewed and non-symmetric mesh arrangement at the crack front. Engineering Fracture Mechanics. 2005;72(11):1717-37. [21] Lonetti P. Dynamic propagation phenomena of multiple delaminations in composite structures. Computational Materials Science. 2010;48(3):563-75. [22] Ramantani DA, Campilho RDSG, de Moura MFSF, Marques AT. Stress and Failure Analysis of Repaired Sandwich Composite Beams using a Cohesive Damage Model. Journal of Sandwich Structures & Materials. 2010;12(3):369-90. [23] Nishioka T. Computational dynamic fracture mechanics. International Journal of Fracture. 1997;86(1-2):127-59. [24] Foo CC, Chai GB, Seah LK. A model to predict low-velocity impact response and damage in sandwich composites. Composites Science and Technology. 2008;68(6):1348-56. [25] Burlayenko VN, Sadowski T. Dynamic behaviour of sandwich plates containing single/multiple debonding. Computational Materials Science. 2011;50(4):1263-8. [26] Kwon YW, Violette MA, McCrillis RD, Didoszak JM. Transient Dynamic Response and Failure of Sandwich Composite Structures under Impact Loading with Fluid Structure Interaction. Applied Composite Materials. 2012;19(6):921-40. [27] Odessa I, Frostig Y, Rabinovitch O. Modeling of interfacial debonding propagation in sandwich panels. International Journal of Solids and Structures. 2017 (in press). [28] Hou CY. Various remeshing arrangements for two-dimensional finite element crack closure analysis. Engineering Fracture Mechanics. 2017;170:59-76. [29] Funari MF, Greco F, Lonetti P. Dynamic debonding in layered structures: A coupled ALEcohesive approach. Frattura ed Integrita Strutturale. 2017;11(41):524-35. [30] Funari MF, Lonetti P. Initiation and evolution of debonding phenomena in layered structures. Theoretical and Applied Fracture Mechanics. 2017;92:133-45. [31] Funari MF, Greco F, Lonetti P. A moving interface finite element formulation for layered structures. Composites Part B: Engineering. 2016;96:325-37. [32] Greco F, Leonetti L, Lonetti P. A novel approach based on ALE and delamination fracture mechanics for multilayered composite beams. Composites Part B: Engineering. 2015;78:447-58. [33] Barbero EJ. Finite Element Analysis of Composite Materials using Abaqus: CRC Press, Taylor & Francis Group; 2013. [34] COMSOL. Multiphysics Reference Guide. 2014. [35] Prasad S, Carlsson LA. Debonding and Crack Kinking in Foam Core Sandwich Beams .2. Experimental Investigation. Engineering Fracture Mechanics. 1994;47(6):825-41. [36] Volokh KY, Needleman A. Buckling of sandwich beams with compliant interfaces. Computers & Structures. 2002;80(14-15):1329-35.
20
LIST OF TABS AND FIGURES
Table.1: Incremental-iterative procedure of the proposed algorithm Table.2: Mode I analysis: mechanical, geometrical and interface properties of the sandwich structure Table.3: Mode II analysis, mechanical, geometrical and interface properties of the sandwich structure. Fig.1: Synoptic representation of the sandwich structure. Fig.2: Moving and referential coordinate systems in the ALE description. Fig.3: Moving mesh interface elements, process zone and crack onset. Fig.4: Loading schemes and geometrical configuration for mode I analysis. Definition of the internal debonding geometry. Fig.5: Mode I analysis, normalized loading vs opening displacement: comparisons with numerical [27] and experimental [35] data in terms of normalized loading curve and opening displacement. Fig.6: Mode I analysis, normalized loading vs opening displacement: effect of the loading rate on the resistance curve: comparisons between static and dynamic loading curves. Fig.7: Mode I analysis, normalized loading vs opening displacement: Effect of the internal debonding length: comparisons between static and dynamic loading curves. Fig.9: Mode I analyses: crack tip speed vs crack displacement as a function of the loading rate and the internal debonding length for C1 and C3 configurations. Fig.10: Loading schemes and geometrical configuration for Mode II analysis. Definition of the internal debonding geometry. Fig.11: Mode II analysis, normalized loading vs opening displacement: comparisons with numerical [27] in terms of normalized loading curve and opening displacement. Fig.12: Mode II analysis: comparisons with numerical data arising from [27] in terms of traction forces for lower and upper interfaces. Fig.13: Mode II analysis, normalized loading vs opening displacement: comparisons between static and dynamic loading curves for C1 configuration. Fig.14: Mode II analysis, normalized loading vs opening displacement: comparisons between static and dynamic loading curves for C2 configuration. Fig.15: Mode II analysis: comparisons in terms of traction forces for different position of the debonding front. Fig.16: Mode II analysis: comparisons in terms of the crack tip speed between C1 and C2 configurations.
21
Table.1: Incremental-iterative procedure of the proposed algorithm START 0. Read the input data: geometry, material and interface characteristics 1. Loop for model coarse iteration 1.1. Initialize the nodal displacement vector 1.2. Loop for load increment 1.2.1.Determine the external force vector f , stiffness, mass and ALE (in case of onset condition) matrixes ( K , M ,W ) 1.2.2.Solve and compute nodal vector variables 1.2.3.Compute stress and strains for each element and evaluate the damage variable of the interfaces 1.2.3.1. Update the load increment counter; if crack onset condition is not satisfied go back to step 1.2 to solve current solution 1.2.3.2.
If crack onset condition is satisfied g kf ( X k ) = 0 with 0 ≤ X k ≤ L , perform STOP 1
1
Solver condition 1.2.3.3. In the X 1k coordinate insert the debonding length ( Ω L ,Ω R ) by modifying the interfacial geometry 1.2.3.4. Perform model refinement and transfer the interfaces damage variables to the new mesh by means a remeshing procedure 1.2.3.5. Activate interface ALE elements concerning the current debonding process 1.2.3.6. Prediction of the crack growth END
Table.2: Mode I analysis: mechanical, geometrical and interface properties of the sandwich structure. s [MPa] E 11
G12s [MPa]
νs
ρs [Kg/mc]
-
0.33
2700
-
Face - sheet aluminum 70 10
3
26 10
3
E 1c [MPa]
c [MPa] G12
νs
ρc [Kg/mc]
-
420
220
0.25
400
-
L [mm]
B [mm]
a [mm]
c [mm]
T s [mm]
Tc [mm]
152.4
25.4
50.8
27.1
2.2
15
Core - PMI AIRES R 90.400
Geometrical properties
22
-1 Gc [N mm ]
∆ 0 [mm]
-
-
-
0.550
0.12
-
-
-
Interface properties
Table.3: Mode II analysis, mechanical, geometrical and interface properties of the sandwich structure. s [MPa] E11
Es22 [MPa]
s [MPa] G12
νs
ρs [Kg/mc]
-
0.28
1800
-
Face - CFRP 135 10
3
9.750 10
3
6 10
3
E 1c [MPa]
c [MPa] G12
νs
ρc [Kg/mc]
-
-
105
42
0.25
75
-
-
L [mm]
B [mm]
a [mm]
c [mm]
Ts [mm]
Tc [mm]
550
250
100
25
2.25
25.7
-1 Gc [N mm ]
∆ 0 [mm]
-
-
-
-
0.385
0.20
-
-
-
-
Core - PMI RHOACELL 71RIST
Geometrical properties
Interface properties
Fig.1: Synoptic representation of the sandwich structure.
23
Fig.2: Moving and referential coordinate systems in the ALE description.
24
Fig.3: Moving mesh interface elements, process zone and crack onset.
Fig.4: Loading schemes and geometrical configuration for mode I analysis. Definition of the internal debonding geometry.
25
Fig.5: Mode I analysis, normalized loading vs opening displacement: comparisons with numerical [27] and experimental [35] data in terms of normalized loading curve and opening displacement.
26
Fig.6: Mode I analysis, normalized loading vs opening displacement: effect of the loading rate on the resistance curve: comparisons between static and dynamic loading curves.
27
Fig.7: Mode I analysis, normalized loading vs opening displacement: Effect of the internal debonding length: comparisons between static and dynamic loading curves.
Fig.8: Mode I analyses: crack tip speed vs crack displacement as a function of the loading rate and the internal debonding length for C1 and C2 configurations.
28
Fig.9: Mode I analyses: crack tip speed vs crack displacement as a function of the loading rate and the internal debonding length for C1 and C3 configurations.
Fig.10: Loading schemes and geometrical configuration for Mode II analysis. Definition of the internal debonding geometry.
29
Fig.11: Mode II analysis, normalized loading vs opening displacement: comparisons with numerical [27] in terms of normalized loading curve and opening displacement.
30
Fig.12: Mode II analysis: comparisons with numerical data arising from [27] in terms of traction forces for lower and upper interfaces.
31
Fig.13: Mode II analysis, normalized loading vs opening displacement: comparisons between static and dynamic loading curves for C1 configuration.
32
Fig.14: Mode II analysis, normalized loading vs opening displacement: comparisons between static and dynamic loading curves for C2 configuration.
Fig.15: Mode II analysis: comparisons in terms of traction forces for different position of the debonding front.
33
Fig.16: Mode II analysis: comparisons in terms of the crack tip speed between C1 and C2 configurations.
34
S0263-8223(18)31357-6 https://doi.org/10.1016/j.compstruct.2018.06.113 COST 9907
To appear in:
Composite Structures
Received Date: Revised Date: Accepted Date:
12 April 2018 21 June 2018 27 June 2018
Please cite this article as: Funari, M.F., Greco, F., Lonetti, P., Sandwich panels under interfacial debonding mechanisms, Composite Structures (2018), doi: https://doi.org/10.1016/j.compstruct.2018.06.113
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Sandwich panels under interfacial debonding mechanisms Marco Francesco Funari, Fabrizio Greco, Paolo Lonetti Department of Civil Engineering, University of Calabria, Via P. Bucci, Cubo39B, 87030, Rende, Cosenza, Italy. ABSTRACT: The paper presents a nonlinear approach to investigate the behavior of composite sandwich structures with transversely compressible core, under static and dynamic loading conditions. The proposed model, formulated in the 2D framework, incorporates moving mesh cohesive modeling, crack initiation and nucleation at core/skin interfaces. Interface elements are used to predict debonding mechanisms, whereas shear deformable beams and two-dimensional plane stress elements identify skin and core behavior, respectively. In this framework, interfacial crack onset, layer kinematic and debonding propagation effects are correctly simulated. The moving mesh technique, combined with a multilayer formulation, ensures a reduction of the computational costs, required to predict crack onset and progressive evolution of debonding phenomena. Cohesive models for sandwich core/skin interfaces are calibrated by means of comparisons with numerical and experimental data with respect mode I and mode II configurations. Moreover, a parametric study to address the influence of the loading rate and sandwich characteristics on both static and dynamic frameworks is proposed.
Keywords: sandwich structures, debonding phenomena, moving mesh technique, dynamic fracture mechanics. Nomenclature
a
initial crack length
b B g kf
length of the internal discontinuities crack growth function
GI
energy release rate mode I
GII
energy release rate mode II
GIC
critical strain energy release rate mode I
GIIC
critical strain energy release rate mode II
g
width of the specimen
k f
crack growth function
GC L Ts
length of the specimen
Tc
thickness of the core
v cs
shear wave speed of the core
2∆
internal discontinuity
∆0
characteristic length parameter of the work of separation
Ω
L
per unit area work of separation
thickness of the skin
process zone on the left of the internal discontinuity
ΩR
process zone on the right of the internal discontinuity
Ωc
core domain
Ωs
skin domain
1
Γ S/C 1.
interface between skin and core
Introduction During the last decades, sandwich systems were object of studies and researches from both
manufacturing and design point of views. In particular, such structures are essentially based on two external high density materials, bonded to a thick core, made from low density material. Their configuration, coupled with an optimized material combination, offers a great variety of lightweight structural systems, which are utilized in aerospace, aircraft, marine industries and, recently, also in civil engineering applications [1]. However, the heterogeneity of the layered systems, in terms of geometry and material properties, may be the cause of macroscopic and microscopic damage phenomena. As a matter of fact, during the production process, the assembly of the structural constituents is an important procedure, in which moisture effects [2, 3], temperature variations and contact with ground support equipment [4, 5] are responsible of pre-manufacturing damage phenomena. Furthermore, during the life-structural service, sandwich structures are subjected to low-velocity impacts. This may cause failure mechanisms, with relevant discontinuities between skin and core [6]. In order to investigate the behavior of sandwich structures, several approaches have been proposed in literature. In particular, the analysis started by means of macroscale approaches, in which analytical and numerical formulations aimed to identify critical failure loads for design purpose were developed [7]. Such procedures were endorsed in equivalent plate/shell layer theories, which utilize, through the thickness, a single displacement approximation [8, 9]. Generalized models were also proposed by means multiscale approaches, working at macro-, meso- or microscales, mainly devoted to predict damage behavior of composite structures [10-12]. In particular, the local behavior, at the skin/core bond lines, requires proper micro-mechanical models to simulate the progressive evolution of the interfacial defects [13]. Moreover, at macroscale, nonlinear material models are utilized to predict shear failure mechanisms, transverse cracking, penetration 2
models and core deformations [14]. Typically, the crack growth is fixed along the interface region because of the material mismatch and low mechanical characteristics of the adhesive layer. However, the proposed model requires further extension to simulate kinking phenomena, by introducing crack growth angle criterion. The analysis of interfacial defects and their evolution on sandwich structures can be broadly classified into Cohesive Zone Models (CZMs) or Finite Fracture Mechanics (FFM). In particular, CZMs introduce interface elements at skin/core debonding line, with effective Traction Separation Law (TSL) constitutive relationships. The interlaminar stresses are described in terms of cohesive strength and relative displacements. Formulations available in literature differ in the definition of the constitutive TSLs, which could be based on non-potential or potential models [15-17]. FE models developed in this context range between continuum and discrete type of elements, in which proper TSLs reproduce the rod internal forces vs nodal separation displacements [18]. CZMs predict directly crack onset and propagation of interfacial delaminations, without introducing preexisting debonding length. In certain FE configurations, CZMs may produce numerical difficulties due to mesh dependence, computing inefficiency, sensitivity to the element aspect ratio [19]. Such computational complexities require advanced analyses, which prevent the use of CZMs for engineering practitioners. Alternatively, FFM allows to identify the crack growth by using Energy Release Rate (ERR) concept, which is determined numerically, e.g. based on the Virtual Crack Closure Technique (VCCT) [20] or the J-integral method [21]. Such analyses assume that an initial crack should be introduced prior of the analysis, in which a negligible size of the fracture process zone and the characteristic length of the structure is considered [22]. Although proper crack initiation criteria, based on energy and stress variables, should be utilized to overcome the requirement of an initial crack length, difficulties regarding the identification of the crack front during the mesh transition or the treatment of the singularity behavior of the stresses are well documented in the literature [23].
3
The above referred papers deal with the analysis of the sandwich structures under static loading conditions. However, an important issue to be achieved is the simulation of low or high speed impact conditions. In dynamics, most of the studies focus on the response of the whole structure, which is described in terms of vibration characteristics. Only few attempts discussed the behavior of sandwich structures, in presence of dynamic debonding phenomena. As a matter of fact, 3D analyses incorporating progressive damage models are carried out, in which low-velocity impact responses are investigated in terms of resistance curves and residual plastic deformations [24]. Free vibrations of sandwich panels are discussed in [25], in which a parametric study on the influence of size, location and number of debonding zones is proposed. Moreover, investigations on the dynamic response of debonded sandwich structures, in conjunction with kinematic constraints, governing contact-impact relationships, are proposed in [26]. Recently, dynamic debonding phenomena on sandwich structures were investigated in [27], in which a general model based on Extended HighOrder Sandwich Panel theory and cohesive theory is implemented to predict interfacial crack propagation under mixed model loading conditions. A common feature of the previous models is the treatment of the process zone ahead of the crack tip. In both cases a large number of computational points in proximity of the process zone is required, in order to capture singularity behavior of the stresses or node release effects during the crack growth. This circumstance is particularly relevant in the framework of sandwich structures because of the material heterogeneity. As a consequence, in order to reduce numerical difficulties in complex geometries during the crack advance, remeshing techniques should be utilized, with severe computational costs in the transition procedures [28]. In the present paper a generalized model to predict dynamic debonding phenomena based on the use of moving mesh method is proposed. In spite of other existing methods, the proposed strategy is able to introduce a low number of computational points in the whole geometry, reducing the complexity of the numerical model. This is achieved by using a moving mesh strategy based on Arbitrary-Lagrange Eulerian (ALE) formulation, which is able to reproduce mesh
4
movements according to the process zone motion, ensuring accuracy in the definition of the fracture variables. The outline of the paper is as follows. Section 2 presents the formulation of the governing equations for the ALE and interface approach, whereas in Section 3 the numerical implementation of the finite element model is reported. Finally, comparisons and parametric results to investigate static and dynamic behavior of debonding phenomena are proposed in Section 4.
2.
Formulation of the model The structure is modelled by means of an assembly of layers, which consists of an internal
core and two external skins. At skin/core bond lines, interlaminar interfaces are introduced to simulate initiation and growth of interfacial defects, which may be preexisting or produced during the load history. This is achieved by the use of interface elements based on moving mesh technique, which ensures an accurate description of the fracture variables and the application of cohesive interlaminar stresses in the process zone. Without loss of generality and for the sake of clarity, a sandwich structure is analyzed, in which fracture phenomena affect the upper skin/core interface only, producing two debonding lengths departing from Right (R) or Left (L) directions. A schematic representation of the model is reported in Fig.1.
2.1 Moving mesh cohesive model Moving interface elements between core and skins are introduced to simulate the onset and evolution of debonding phenomena. In particular, the crack growth is expressed as a function of two coordinate systems, i.e. referential and moving, which parametrize the motion of the process zone from the onset to the crack advance (Fig.2). From a mathematical point of view, a computational (Moving) mesh is implemented as the image of the grid points in the referential (Fixed) configuration. In particular, the prescribed motion is expressed in terms of the following Laplacebased equations developed for Static (S) or Dynamic (D) frameworks: 5
∆X
k d ,XX
∂ 2Φ k ( X ,t ) = =0 ∂X 2
∂ 3Φ k ( X ,t ) k =0 ∆ X d ,XX = ∂t ∂X 2
(S)
(D),
(1)
where k, with k=L, R, represents the index referred to the Left (L) or Right (R) process zones and
∆ X dk = X dk ( t ) − X k corresponds to the mesh displacement function [29]. Eqs. (1) require boundary and initial conditions (in dynamics), which depend from the motion of the process zone lengths. ALE interface elements are incorporated in the sandwich structure where debonding phenomena may occur (Fig.3). However, the position, where the crack growth is made possible, is identified on the basis of a proper crack growth function ( g kf ) , defined in terms of mixed mode ERR ratio [15] as follows:
( )
g kf X 1k
( GIC ,GIIC ) are
where
GI X k 1 = GIC
( )
r
2 GII X k 1 + GIIC
( )
r
2 − 1,
(2)
the total area under the traction separation laws, i.e. critical values, and ∆nc
( GI ,GII )
( )
are the individual energy release rates calculated as GI = ∫ Tn ∆n d ∆n and 0
∆tc
GII = ∫ Tt ( ∆t ) d ∆t with ∆n
and
∆s
representing the normal and tangential interlaminar
0
displacements respectively. For each mode components, the TSLs, namely Tn ( ∆n ) or Tt ( ∆t ) , could be described in terms of bilinear, trapezoidal, nonlinear relationships depending from the material interface characteristics. Without loss of generality, a linear-softening curve is considered in the present study. However, a generalization to other forms of damage laws would not require specific changes to the proposed modeling, which is independent from the assumed TSL. Similarly, a different crack onset criterion from Eq.(2) should be introduced to identify the position, in which debonding phenomena are activated.
6
The prescribed mesh motion introduces a rigid displacement of the process zone, which is identified on the basis of internal lengths for the left and right crack directions, namely Ω L or Ω R . Starting from the position ( X T ) in which onset condition is satisfied, two independent coordinates
(X
L T
, X TR ) are generated, which simulate the possibility for the crack tip to advance along the right
and left directions, respectively. The left and right debonding lengths are identified in terms of the positions
(X
L T
R
)
, X T , in which the fracture function g f assumes negative values, within a fixed
tolerance, as follows: L
Ω L = X TL − X T
R
Ω R = X TR − X T
with
g f ( X TL ) = 0
and
g f X T ≤ 0,
( )
(3)
with
g f ( X TR ) = 0 and
g f X T ≤ 0,
( )
(4)
L
R
Finally, the boundary conditions associated to Eq.(1), are required to simulate the displacements of the process zone for the left and right debonding lengths and at geometry edges for the mesh motion:
( )
(5)
( )
(6)
L
∆ X1L ( L ) = 0, ∆ X1L ( X TL ) = ∆ X1k X T = δ L , R
∆ X1R ( 0) = 0, ∆ X1R ( X TR ) = ∆ X1R X T = δ R ,
where (δ L ,δ R ) are the mesh displacements obtained by solving the following optimality conditions: +
+
+
+
δ L ≥ 0 ,g kf ( X TL ) ≤ 0 with δ L ⋅ g kf ( X TL ) = δ L ⋅ g kf ( X TL ) = 0, δ R ≥ 0 ,g kf ( X TR ) ≤ 0 with δ R ⋅ g kf ( X TR ) = δ R ⋅ g kf ( X TR ) = 0.
(7)
From the numerical point of view, the displacements of the debonding regions are determined on the basis of the values of the fracture function at the extremities of the debonding length. In particular, during the crack growth a null value of the fracture function is enforced at the internal extremity of the debonding length. Therefore, with reference to k-th defect, the crack tip displacement can be expressed by means of the following relationship: 7
∆X
kL ,R T
=0
g
kL ,R f
≤ 0, ∆ X
kL ,R T
=
g kf ( X 1kL ,R ) + toll g kf ( X 1kL ,R ) + toll ± g kf ( X 1kL ,R ± Ω kL ,R )
,R =0 Ω kL ,R g kL f
(8)
where toll is tolerance positive value able to trigger the crack growth. From the numerical point of view it has been chosen close to zero, i.e. 1E-2, to do not modify the interfacial cohesive laws [30]. 2.2 Kinematic formulation of the model
The skins follow a one dimensional model based on Timoshenko beam kinematics, whereas the core is modelled by using a 2D plane stress formulation. The materials of each layer are assumed to be linear elastic and hypotheses concerning small displacements and strains are considered in the analyses. Moving interface elements are introduced between the skin and core to reproduce potential debonding phenomena. The interfaces are initially assumed to be perfectly bonded to both upper and lower skin/core lines and thus moving coordinates match with the referential or material ones [30]. In this framework, dynamic equations of motion can be expressed by introducing equilibrium conditions along the horizontal and vertical directions and flexural rotation (counterclockwise) as follows: − µU2k + EAU 2k,X ,X + FX ( X ) + p x ( X ) = 0 , − µV2k + GA* (V2k,X − Φ 2k ) ,X + FY ( X ) + p y ( X ) = 0 , − I 0Φ2k + EIΦ 2k.XX − GA* (V2k,X − Φ 2k ) + m ( X ) = 0 ,
(
on Ω S , Ω S
div [ E∇U ] + f − µU = 0,
on
E33 (U ,YC + V,XC ) = FX ( X )
on
(E
C C 12U ,X + E22V,Y ) = FY ( X )
1
ΩC
(Γ
S1 / C
2
)
(9)
(10)
.Γ S2 / C
)
(11)
where, as shown in Fig.3, the superscript k (with k=S1,S2), here and in the following, refers to the lower (S1) or upper skins (S2), ( Ω X ,Γ X / Y ) are the domain of the X or the boundary between X and Y (with X,Y=S1, S2, C, ALE1, ALE2), ( I 0 , µ ) are the per unit length inertial rotation and mass of the skin, ( EI ,GA* ) are the flexural and shear stiffnesses, 8
( p , p ,m ) are the axial, vertical and flexural x
y
d(X) external loads applied to the skins, the notation X = indicates the time derivative of the dt
( )
function X, Eij (with i,j=1,2,3) is the 2D plane stress elastic matrix, U (with U T = U C V C ) is the vector containing the horizontal and vertical displacements of the core, f is the vector force per
unit volume. Moreover, the interlaminar stresses are defined on the basis of moving mesh coordinates and interlaminar relative displacements as follows:
FX ( X1k ) = Tn ∆nk ( X1k ) FY ( X1k ) = Tt ∆tk ( X1k )
,
on
Γ S1 / C and
Γ S2 / C
(12)
,
(13)
with
∆nS ( X ) = V C ( X , −0.5 ⋅ T C ) − V S ( X ) , 1
1
∆nS ( X ) = V S ( X ) − V C ( X ,0.5 ⋅ T C ) , 2
2
∆tS ( X 1R ) = U C ( X , −0.5 ⋅ T C ) − U S ( X 1R ) − Φ S ( X 1R )( 0.5 ⋅ T S ) 1
1
1
∆tS ( X 1R ) = U S ( X 1R ) + Φ S ( X 1R )( 0.5T S ) − U C ( X ,0.5 ⋅ T C ) 2
s
2
where T S and T C are the thickness of the skin and core, respectively. Previous equations, i.e. Eq.(9), should be considered in those cases, in which the interfaces are not affected by debonding phenomena. However, once crack initiation condition, defined by Eq.(2), is satisfied, at the position
X = X T debonding phenomena should be simulated from the right and left directions. In particular,
(
)
L R two different positions, i.e. X = X T , XT , which identify fictitious crack tips for the right and left
debonding lengths, are introduced. Such quantities differ from the internal length 2∆ , which refers to the initial defect predicted by the crack growth criterion. From the numerical point of view, they are assumed very small, in order to simulate the presence of an initial defect. For more details, the influence of such parameter was investigated in [30]. Therefore, the governing equations are expressed as follows:
9
− µU2k + EAU 2k,X ,X + FX ( X R ) + FX ( X L ) + px ( X ) = 0 , − µV2k + GA* (V2k,X − Φ 2k ) ,X + FY ( X R ) + FY ( X L ) + p y ( X ) = 0 , − I 0Φ2k + EIΦ 2k.XX − GA* (V2k,X − Φ 2k ) + m ( X ) = 0,
on Ω S 2
(14)
where the superscript k is k=S 2, ( X R ,X L ) are the effective moving coordinates for the left (L) and right (R) directions measured on the material domain. Moreover, the governing equations of the core, i.e. field and initial conditions, refer to a classical 2D plane stress modeling, which, for the sake of brevity, is not reported in detail in the paper. Boundary conditions concerning the interface regions between core and skins are defined by the following expressions:
E33 (U ,Y + V,X ) = FX ( X R ) + FX ( X L )
(E
U ,X + E22V,Y ) = FY ( X R ) + FY ( X L )
on Γ S2 / C
(15)
12
It is worth noting that interface traction forces are expressed in terms of the moving coordinate systems, i.e.
(X
R
, X L ) , for the right and left debonding phenomena, leaving unaltered the
governing equations of the structural elements, with respect to the standard material formulation [31].
3.
Numerical implementation
Governing equations presented in previous section are implemented by using a finite element approach. For the sake of conciseness, the implementation procedure is mainly referred to the ALE interface elements and moving traction forces acting at the skin/core interfaces only, since the governing equations of the structural model remain basically unaltered. In particular, Lagrange interpolation functions with cubic or quadratic orders are assumed for the Timoshenko or plane stress finite elements, whereas Lagrange linear interpolation functions are utilized for the moving interface variables:
10
NeS
q S = (U ,V ,Φ ) = ∑ NiS ( X ) (UiS ,Vi S ,Φ iS ) = N S1,2 U i =1 NeC
S1 ,2
(
on Ω S , ΩS
,
q = (U ,V ) = ∑ NiC ( X ) (U iC ,ViC ) = N C U , i =1 1,2
C
R
NeMM
X =
∑ F ( X ) ( X
R
i
i =1
) = F X
R
C
,X =
∑ F ( X ) ( X i
i =1
) = F X
L
(
1
(
1
, on Ω ALE ,Ω ALE
∑ F ( X ) ( Λ ) = F Λ , i
) (16)
L
NeMM
Λ=
2
on ΩC
NeMM
L
1
on Ω ALE ,Ω ALE
i
i =1
2
)
2
)
where • contains the nodal quantities of the [• ] function, N eS , N eC ,N eMM are the number nodes
(
(
)
)
in the skin, core and ALE elements, respectively, N S1 ,S2 ,N C ,F are the matrixes containing the
(
)
interpolation functions for the structural and ALE elements, q S1,2 ,q C collect nodal variables of the structural problem for the skin and core and
( X
R
, X L ) are the interpolation functions of the right
and left moving debonding length. Finally, Lagrange multiplier method is introduced as a result of weak constraint conditions by function Λ , which is assumed with a linear interpolation. Starting from Eq. (1) and Eq.(14), introducing weighing functions related to translational displacements (WUS,V ,WUC,V ) for the skins and core, rotation (WΦS ) fields for skins and positional variables (W k ) for the ALE variables, the weak forms for interface elements and the moving traction forces on the skin/core interfaces present the following relationships:
(∆X
k ,X
,W,Xk )
Ω ALE1
k C V Γ
+ (Tn ,W − W S V
(
= Λ1k δ ( X Tk ) ,W k
)
( ( ) ) k
Ω ALE1
k
+ Λ1 δ X T ,W k
+ Ω ALE1
k
)
S1 / C
T + Tt ,WUS + WΦS − WUC + 2 Γ S / C
(17)
1
(
+ ∆ X Tk − δ k δ ( X Tk ) ,WΛk 1
)
Ω ALE1
(
( ) )
k k + ∆ X T − δ k δ X T ,WΛ 2
=0 Ω ALE1
with Tnk = C nk ∆n ( X k ) ∆n ( X k ) , Tt k = Ctk ∆t ( X k ) ∆t ( X k ) 11
(18)
where the notation ( ⋅,⋅)Ω denotes ( ⋅,⋅)Ω = ∫ ( ⋅,⋅)dX , δ ( i ) is the delta dirac function, the superscript k Ω
(with k=R,L) refers to the Right and Left interfaces, and ( ∆s , ∆n ) are the interlaminar displacements defined in Eq.(13). Eqs.(18) are rearranged in vector notation as follows: T T = t Tn
Dt ( ∆t = 0
)
∆ t = D∆ Dn ( ∆n ) ∆n 0
(19)
where Dt ( ∆t ) ,Dn ( ∆n ) are the tangential and normal stiffnesses of the introduced to describe nonlinear TSLs. Moreover, the interlaminar displacement vector, ∆ , can be expressed by
(
introducing a mixed interpolation function matrix, namely N SC and two auxiliary matrixes Q1 ,Q2
)
:
T ∆t U S1 + Φ S1 − U C SC = = Q1 N S U S − Q2 N C U C = N SC U 2 ∆n V S1 − V C where U
SC
(with U
SC
C
(20)
S
= U ∪ U ) is the vector containing the nodal displacements of the core and
skin element. Since the interlaminar foces are evaluated on the moving domain, a mapping operator is required to identify the connection between referential and material coordinate systems[32]. To this end, a projector operator Φ : X → X d is introduced, which the relationship between the nodal
coordinates in the moving and referential configurations, i.e X d and X , with X d = Φ X . Substituting Eq,(16) into Eq.(17) and by using Eq.(20), the following discrete equations in which unknown quantities are represented by the position of the computational nodes, are derived: k
Ak X d + B1k λ1k + B2k λ2 + C CS ( k )U SC( k ) = 0 , k k k k k k B1 X d − C1 = 0 , B 2 X d − C 2 = 0 ,
with
12
(21)
X i +1 A = ∫ F,X ⊗ F,X dX ,B1 = δ ( X T ) F T ,B2 = δ ( X T ) F T X T X i
( )
C
CS
Xi +1 = ∫ ( N SC Φ X i
T
) D ( N
SC
Φ
) dX ,
(22)
( )(∆X
C1 = δ ( X T )( ∆ X T − δ f ) , C2 = δ X T
k T
− δ f ) ,
where ⊗ is the dyadic product between vectors and ( ) is the composition function. Eqs.(21), which are not modified from the classical definition [33], are implemented by using a customized FE subroutine in the framework of COMSOL Multiphysics, by means scripting capabilities of MATLAB® language [34]. In particular, the analysis is based on several steps, which are connected by using stop/restart procedures. In particular, the built-in Newton-Raphson or implicit transient analyses based backward differentiation formula (BDF) for static and dynamic analyses are employed. These are managed by means of external script files based on stop condition for crack initiation and restart/remeshing procedures for the activation of moving interface elements. A synoptic representation of the computational algorithm implemented in the FE environmental program is reported in Tab.1.
4.
Validation and numerical results
In this section, the proposed model is verified by means of comparisons with numerical and experimental data. The first step is developed with the purpose to analyze the consistency of the proposed formulation under mode I and mode II loading conditions. In particular, according to standard experimental methods, the static behavior of interfacial cracks for several sandwich debonding configurations is investigated. The main aim of the comparisons is to validate the proposed model and to examine its ability to describe debonding failure mechanisms in sandwich configurations. Subsequently, the structural behavior is investigated in both static and dynamic
13
frameworks to identify the influence of inertial effects on the crack growth phenomena, produced by different levels of the loading rate. 4.1 Mode I DCB loading scheme
At first, analyses are developed with reference to a classical DCB loading involving a pure mode I fracture analysis. The loading, boundary conditions and geometry of the specimen are illustrated in Fig.4, whereas, according to data recovered in [27, 35], mechanical properties assumed for the laminate, core and interfaces as well as the ones required by the cohesive elements, are summarized in Tab.2. The crack evolves in the interface region between core and face sheet (namely at Γ S2 / C ). In particular, a preexisting delaminated length is assumed in the upper interface, whereas the lower interface is considered as perfectly bonded. The numerical model is discretized by means of shear deformable beam elements with cubic interpolation functions for each face sheet, whereas the core is modeled by quadratic interpolation functions with linear elastic material behavior. Moreover, ALE elements with linear interpolation functions are introduced between each sheet-core interfaces to predict crack onset or debonding advance. The analysis is developed under a displacement control mode to ensure a stable crack propagation. Cohesive constitutive relationships are based on exponential laws, with smooth decay branches according to [36]. The numerical discretization utilized for the face sheet is assumed to be generally uniform with a length equal to ∆M / L = 1 / 30 , with ∆M the element length. The core has been discretized by means plane stress quadrilateral elements with maximum element length equal to ∆M / L = 1 / 30 . At the interfaces a coarse discretization is adopted, except at the process zone ( Ω R ,L ) where a mesh enrichment is considered, i.e. ∆M / L = 1 / 300 . In Fig. 5, comparisons in terms of loading curves as a function of opening end displacements, for two different core thickness configurations, i.e. Tc=15,20 mm, are reported. The results show that the proposed approach is able to reproduce correctly the behavior of the structure, since the loading curves are quite in agreement with 14
experimental [35] and analytical data available in [27]. The results show that an increment in the core thickness does not produce improvements in the loading curves, which are basically coincident in terms of stiffness, displacements except for a small variation of the peak load. Previous analyses developed essentially in statics are extended in a dynamic framework. The main aim is to investigate the influence of the loading rate on both resistance curve and crack tip speed evolution. In addition, the effect of an internal preexisting delamination region is also considered to verify the dynamic amplifications produced during the crack growth and their interaction with an external debonding length. The structure is subjected to an applied speed with a ramp curve until the time (t0), which is assumed proportional to the first period of vibration of the structure (t0=0.5T1); subsequently a fixed velocity is prescribed to the end displacement (v0). In Fig.6, comparisons in terms of resistance curve between static and dynamic results are developed, obtained by means several prescribed opening speeds with refers to low or high ranges. Moreover, the following configurations are considered in the analyses: •
external debonding length only with b=0 (C1);
•
internal and external discontinuities with different geometric lengths, i.e. b/a=0.2 (C2) or b/a=0.5 (C3).
The results show that under static loading the presence of an internal debonding length modifies load-displacement curves, increasing system instability during the crack growth and leading large oscillations with respect to the static case. The presence of an internal debonding, at low speeds, does not modify the evolution curve, which, apart little oscillations in the unstable path, is quite coincident with the C1 configuration. However, at large speeds, the dynamic effects modify the crack growth behavior, since C1 and C2 denote notable oscillations, which tend to be reduced, for increasing values of the opening displacements, leading to same predictions of the static solution. In Fig. 7, the influence of an internal discontinuity on the dynamic behavior of debonding phenomena is discussed. In particular, for a fixed loading rate, comparisons in terms of loading 15
curve are proposed for the C1, C2 and C3 configurations. The results show how the presence of an internal discontinuity modifies the equilibrium path, introducing a marked span-back phenomenon in the loading curve. However, once the crack tip overcomes the internal debonding length, the curves for the configurations with internal debonding lengths, i.e. C2 and C3, tend to the same prediction. Additional results are presented in Figs. 8-9, in which the relationship between speed and position of the crack tip is investigated. The analyses show that during the crack onset, the tip is affected by large accelerations, which tend to be reduced as far as the crack tip speed grows. However, at low loading rates, after a maximum value is reached, the crack tends to a steady state advance, since a constant speed is observed. In presence of an internal debonding region, the crack is affected by an oscillating behavior and, especially at low loading rates, crack arrest phenomena are observed. This phenomenon is mainly produced by the presence of the internal debonded region, which produces a local distribution of contact forces with respect to the case of perfectly debonded region. Large amplifications in the crack tip speeds are observed, especially, during the transition regions from the adhesive to the delaminated length.
4.2 Mode II CSB loading scheme
The analyses are extended to a loading scheme involving a mode II loading condition. Comparisons with numerical data are reported for a sandwich configuration based on a three-point bending loading scheme. As shown in Fig.10, the initial crack is located in the upper interface between the face-sheet and the core, whereas geometrical and mechanical characteristics are reported in Tab.3. The numerical discretization utilized for the face sheet is assumed to be generally uniform with a length equal to ∆M / L = 1 / 55 . The core has been discretized by means plane stress quadrilateral elements with maximum element length equal to ∆M / L = 1 / 55 . At the interfaces a coarse discretization is adopted, except at the process zone ( Ω R ,L ) where a mesh enrichment is considered, i.e. ∆M / L = 1 / 550 . The following configurations are considered in the analyses: 16
•
external debonding length, i.e. b=0, (C1);
•
internal and external discontinuities with b/a=0.5 (C2).
Comparisons in terms of applied loads as a function of mid-span vertical displacement for two different preexisting crack length configurations are reported. In both cases, results obtained by the proposed model are in agreement with the ones obtained by analytical approach based on pure cohesive model [27]. Moreover, in Fig. 12, comparisons are developed to verify the consistency of the solution with respect to the local distribution of interfacial stresses. Despite numerical data reported in the literature, the proposed model is able to reproduce the crack growth also for more extended tip displacements. The results point out in the load-displacement curves the presence of a knee once the position of the crack tip overcomes the one of the applied loads. The accuracy of the solutions is guaranteed also locally in terms of stress distribution. Good agreement with the values obtained by using a pure cohesive modeling, in which an accurate description in terms of mesh element length is required along all the interface regions is found. Contrarily in the proposed model, based on ALE approach, the mesh refinement is introduced along the process zone only, leading to a reduction of the computational costs involved in the analysis. Previous results developed essentially in statics are extended in a dynamic framework, in which applied speeds based on a ramp curve with different loading rates are considered. Resistance curves are reported in Figs. 13-14, which refer to configurations with or without an internal debonding length, i.e. C1 and C2, respectively. The results show that, at high loading rates, large oscillations and amplifications with respect to the static solution are observed. Similarly, to results obtained under pure mode I loading scheme, Mode II configuration is also affected by the presence of an internal discontinuity, since an instable behavior in the loading curve with small oscillations is observed with respect to the static equilibrium path. The stress distributions for two different structural scenarios related to C1 and C2 configurations are reported in Fig. 15. The results show that, in presence of an internal debonding length, the traction forces present a discontinuous behavior at the extremity points, leading to singular values of the interfacial stresses. In Fig.16, 17
comparisons in terms of crack tip speeds, denote how the presence of an internal discontinuity produces an amplifications of the tip speed also for low loading rates, leading to jumps of crack speeds once the tip reaches the debonded length. Moreover, as far as the crack tip overcomes the delaminate length; the crack tip speed oscillates leading to crack arrest phenomena during the crack evolution.
5.
Conclusions
The present study was aimed to investigate the behavior of sandwich structures in presence of debonding phenomena. The proposed model combines ALE formulation with a cohesive fracture approach to simulate debonding phenomena in presence of multiple delaminations. Based on the results obtained, the following conclusions can be drawn: •
the proposed methodology is able to reduce computational difficulties arising from the mesh dependence, convergence and accuracy of the solution, which typically affect cohesive approaches;
•
the results show the capabilities of the proposed model to reproduce correctly the behavior in both static and dynamic frameworks;
•
the parametric study denotes how debonding phenomena are quite influenced by the presence of internal material discontinuities, which produce strong oscillations and instabilities in the resistance curve and high speeds during the crack evolution;
•
this study could be very useful for developing modeling strategies to understand debonding phenomena in sandwich composites subject to low-velocity impact;
•
the proposed approach requires a further extension to predict kinking phenomena, since at this stage the crack position is fixed along the interfaces.
Acknowledgments
This work is supported by Italian Ministry of University and Research (P.R.I.N. National Grant 2015, B86J16002300001). 18
References
[1] Manalo A, Aravinthan T, Fam A, Benmokrane B. State-of-the-Art Review on FRP Sandwich Systems for Lightweight Civil Infrastructure. Journal of Composites for Construction. 2017;21(1). [2] Yu H, Zhou C. Sandwich diffusion model for moisture absorption of flax/glass fiber reinforced hybrid composite. Composite Structures. 2018;188:1-6. [3] Spadea S, Orr J, Ivanova K. Bend-strength of novel filament wound shear reinforcement. Composite Structures. 2017;176:244-53. [4] Ehsani M, editor ASCE Innovation Award Winner: Sandwich Construction Carbon FRP Pipe2017: American Society of Civil Engineers (ASCE). [5] Maletta C, Sgambitterra E, Niccoli F. Temperature dependent fracture properties of shape memory alloys: Novel findings and a comprehensive model. Scientific Reports. 2016;6. [6] Yu RC, Pandolfi A. 15 - Modeling of delamination fracture in composites: a review A2 Sridharan, Srinivasan. Delamination Behaviour of Composites: Woodhead Publishing; 2008. p. 429-57. [7] Triantafillou TC, Gibson LJ. Failure Mode Maps for Foam Core Sandwich Beams. Materials Science and Engineering. 1987;95:37-53. [8] Abrate S, Ferrero JF, Navarro P. Cohesive zone models and impact damage predictions for composite structures. Meccanica. 2015;50(10):2587-620. [9] Bruno D, Greco F, Lonetti P. Dynamic Mode I and Mode II Crack Propagation in Fiber Reinforced Composites. Mechanics of Advanced Materials and Structures. 2009;16(6):442-55. [10] Greco F, Leonetti L, Lonetti P, Nevone Blasi P. Crack propagation analysis in composite materials by using moving mesh and multiscale techniques. Computers and Structures. 2015;153:201-16. [11] Ascione L, Berardi VP, Giordano A, Spadea S. Pre-buckling imperfection sensitivity of pultruded FRP profiles. Composites Part B: Engineering. 2015;72:206-12. [12] Ciantia MO, Arroyo M, Butlanska J, Gens A. DEM modelling of cone penetration tests in a double-porosity crushable granular material. Computers and Geotechnics. 2016;73:109-27. [13] Rabinovitch O, Hamed E. Bending behavior of sandwich panels with a "soft" core and embedded rigid inserts. Sandwich Structures: Advancing with Sandwich Structures and Materials. 2005:261-70. [14] Mohammadi MS, Nairn JA. Balsa sandwich composite fracture study: Comparison of laminated to solid balsa core materials and debonding from thick balsa core materials. Composites Part B-Engineering. 2017;122:165-72. [15] Bruno D, Greco F, Lonetti P. Computation of energy release rate and mode separation in delaminated composite plates by using plate and interface variables. Mechanics of Advanced Materials and Structures. 2005;12(4):285-304. [16] Park K, Paulino GH. Cohesive Zone Models: A Critical Review of Traction-Separation Relationships Across Fracture Surfaces. Applied Mechanics Reviews. 2011;64(6). [17] Funari MF, Greco F, Lonetti P, Luciano R, Penna R. An interface approach based on moving mesh and cohesive modeling in Z-pinned composite laminates. Composites Part B-Engineering. 2018;135:207-17. [18] Xu G, Yan RJ. The Use of Sprint Interface Element Delamination Simulation of Sandwich Composite Beam. Applied Composite Materials. 2017;24(5):1049-60. [19] Rabinovitch O. Debonding analysis of fiber-reinforced-polymer strengthened beams: Cohesive zone modeling versus a linear elastic fracture mechanics approach. Engineering Fracture Mechanics. 2008;75(10):2842-59.
19
[20] Okada H, Higashi M, Kikuchi M, Fukui Y, Kumazawa N. Three dimensional virtual crack closure-integral method (VCCM) with skewed and non-symmetric mesh arrangement at the crack front. Engineering Fracture Mechanics. 2005;72(11):1717-37. [21] Lonetti P. Dynamic propagation phenomena of multiple delaminations in composite structures. Computational Materials Science. 2010;48(3):563-75. [22] Ramantani DA, Campilho RDSG, de Moura MFSF, Marques AT. Stress and Failure Analysis of Repaired Sandwich Composite Beams using a Cohesive Damage Model. Journal of Sandwich Structures & Materials. 2010;12(3):369-90. [23] Nishioka T. Computational dynamic fracture mechanics. International Journal of Fracture. 1997;86(1-2):127-59. [24] Foo CC, Chai GB, Seah LK. A model to predict low-velocity impact response and damage in sandwich composites. Composites Science and Technology. 2008;68(6):1348-56. [25] Burlayenko VN, Sadowski T. Dynamic behaviour of sandwich plates containing single/multiple debonding. Computational Materials Science. 2011;50(4):1263-8. [26] Kwon YW, Violette MA, McCrillis RD, Didoszak JM. Transient Dynamic Response and Failure of Sandwich Composite Structures under Impact Loading with Fluid Structure Interaction. Applied Composite Materials. 2012;19(6):921-40. [27] Odessa I, Frostig Y, Rabinovitch O. Modeling of interfacial debonding propagation in sandwich panels. International Journal of Solids and Structures. 2017 (in press). [28] Hou CY. Various remeshing arrangements for two-dimensional finite element crack closure analysis. Engineering Fracture Mechanics. 2017;170:59-76. [29] Funari MF, Greco F, Lonetti P. Dynamic debonding in layered structures: A coupled ALEcohesive approach. Frattura ed Integrita Strutturale. 2017;11(41):524-35. [30] Funari MF, Lonetti P. Initiation and evolution of debonding phenomena in layered structures. Theoretical and Applied Fracture Mechanics. 2017;92:133-45. [31] Funari MF, Greco F, Lonetti P. A moving interface finite element formulation for layered structures. Composites Part B: Engineering. 2016;96:325-37. [32] Greco F, Leonetti L, Lonetti P. A novel approach based on ALE and delamination fracture mechanics for multilayered composite beams. Composites Part B: Engineering. 2015;78:447-58. [33] Barbero EJ. Finite Element Analysis of Composite Materials using Abaqus: CRC Press, Taylor & Francis Group; 2013. [34] COMSOL. Multiphysics Reference Guide. 2014. [35] Prasad S, Carlsson LA. Debonding and Crack Kinking in Foam Core Sandwich Beams .2. Experimental Investigation. Engineering Fracture Mechanics. 1994;47(6):825-41. [36] Volokh KY, Needleman A. Buckling of sandwich beams with compliant interfaces. Computers & Structures. 2002;80(14-15):1329-35.
20
LIST OF TABS AND FIGURES
Table.1: Incremental-iterative procedure of the proposed algorithm Table.2: Mode I analysis: mechanical, geometrical and interface properties of the sandwich structure Table.3: Mode II analysis, mechanical, geometrical and interface properties of the sandwich structure. Fig.1: Synoptic representation of the sandwich structure. Fig.2: Moving and referential coordinate systems in the ALE description. Fig.3: Moving mesh interface elements, process zone and crack onset. Fig.4: Loading schemes and geometrical configuration for mode I analysis. Definition of the internal debonding geometry. Fig.5: Mode I analysis, normalized loading vs opening displacement: comparisons with numerical [27] and experimental [35] data in terms of normalized loading curve and opening displacement. Fig.6: Mode I analysis, normalized loading vs opening displacement: effect of the loading rate on the resistance curve: comparisons between static and dynamic loading curves. Fig.7: Mode I analysis, normalized loading vs opening displacement: Effect of the internal debonding length: comparisons between static and dynamic loading curves. Fig.9: Mode I analyses: crack tip speed vs crack displacement as a function of the loading rate and the internal debonding length for C1 and C3 configurations. Fig.10: Loading schemes and geometrical configuration for Mode II analysis. Definition of the internal debonding geometry. Fig.11: Mode II analysis, normalized loading vs opening displacement: comparisons with numerical [27] in terms of normalized loading curve and opening displacement. Fig.12: Mode II analysis: comparisons with numerical data arising from [27] in terms of traction forces for lower and upper interfaces. Fig.13: Mode II analysis, normalized loading vs opening displacement: comparisons between static and dynamic loading curves for C1 configuration. Fig.14: Mode II analysis, normalized loading vs opening displacement: comparisons between static and dynamic loading curves for C2 configuration. Fig.15: Mode II analysis: comparisons in terms of traction forces for different position of the debonding front. Fig.16: Mode II analysis: comparisons in terms of the crack tip speed between C1 and C2 configurations.
21
Table.1: Incremental-iterative procedure of the proposed algorithm START 0. Read the input data: geometry, material and interface characteristics 1. Loop for model coarse iteration 1.1. Initialize the nodal displacement vector 1.2. Loop for load increment 1.2.1.Determine the external force vector f , stiffness, mass and ALE (in case of onset condition) matrixes ( K , M ,W ) 1.2.2.Solve and compute nodal vector variables 1.2.3.Compute stress and strains for each element and evaluate the damage variable of the interfaces 1.2.3.1. Update the load increment counter; if crack onset condition is not satisfied go back to step 1.2 to solve current solution 1.2.3.2.
If crack onset condition is satisfied g kf ( X k ) = 0 with 0 ≤ X k ≤ L , perform STOP 1
1
Solver condition 1.2.3.3. In the X 1k coordinate insert the debonding length ( Ω L ,Ω R ) by modifying the interfacial geometry 1.2.3.4. Perform model refinement and transfer the interfaces damage variables to the new mesh by means a remeshing procedure 1.2.3.5. Activate interface ALE elements concerning the current debonding process 1.2.3.6. Prediction of the crack growth END
Table.2: Mode I analysis: mechanical, geometrical and interface properties of the sandwich structure. s [MPa] E 11
G12s [MPa]
νs
ρs [Kg/mc]
-
0.33
2700
-
Face - sheet aluminum 70 10
3
26 10
3
E 1c [MPa]
c [MPa] G12
νs
ρc [Kg/mc]
-
420
220
0.25
400
-
L [mm]
B [mm]
a [mm]
c [mm]
T s [mm]
Tc [mm]
152.4
25.4
50.8
27.1
2.2
15
Core - PMI AIRES R 90.400
Geometrical properties
22
-1 Gc [N mm ]
∆ 0 [mm]
-
-
-
0.550
0.12
-
-
-
Interface properties
Table.3: Mode II analysis, mechanical, geometrical and interface properties of the sandwich structure. s [MPa] E11
Es22 [MPa]
s [MPa] G12
νs
ρs [Kg/mc]
-
0.28
1800
-
Face - CFRP 135 10
3
9.750 10
3
6 10
3
E 1c [MPa]
c [MPa] G12
νs
ρc [Kg/mc]
-
-
105
42
0.25
75
-
-
L [mm]
B [mm]
a [mm]
c [mm]
Ts [mm]
Tc [mm]
550
250
100
25
2.25
25.7
-1 Gc [N mm ]
∆ 0 [mm]
-
-
-
-
0.385
0.20
-
-
-
-
Core - PMI RHOACELL 71RIST
Geometrical properties
Interface properties
Fig.1: Synoptic representation of the sandwich structure.
23
Fig.2: Moving and referential coordinate systems in the ALE description.
24
Fig.3: Moving mesh interface elements, process zone and crack onset.
Fig.4: Loading schemes and geometrical configuration for mode I analysis. Definition of the internal debonding geometry.
25
Fig.5: Mode I analysis, normalized loading vs opening displacement: comparisons with numerical [27] and experimental [35] data in terms of normalized loading curve and opening displacement.
26
Fig.6: Mode I analysis, normalized loading vs opening displacement: effect of the loading rate on the resistance curve: comparisons between static and dynamic loading curves.
27
Fig.7: Mode I analysis, normalized loading vs opening displacement: Effect of the internal debonding length: comparisons between static and dynamic loading curves.
Fig.8: Mode I analyses: crack tip speed vs crack displacement as a function of the loading rate and the internal debonding length for C1 and C2 configurations.
28
Fig.9: Mode I analyses: crack tip speed vs crack displacement as a function of the loading rate and the internal debonding length for C1 and C3 configurations.
Fig.10: Loading schemes and geometrical configuration for Mode II analysis. Definition of the internal debonding geometry.
29
Fig.11: Mode II analysis, normalized loading vs opening displacement: comparisons with numerical [27] in terms of normalized loading curve and opening displacement.
30
Fig.12: Mode II analysis: comparisons with numerical data arising from [27] in terms of traction forces for lower and upper interfaces.
31
Fig.13: Mode II analysis, normalized loading vs opening displacement: comparisons between static and dynamic loading curves for C1 configuration.
32
Fig.14: Mode II analysis, normalized loading vs opening displacement: comparisons between static and dynamic loading curves for C2 configuration.
Fig.15: Mode II analysis: comparisons in terms of traction forces for different position of the debonding front.
33
Fig.16: Mode II analysis: comparisons in terms of the crack tip speed between C1 and C2 configurations.
34