# Linear and Nonlinear Dynamic Analyses of Sandwich Panels with Face Sheet-to-Core Debonding

Authors: Vyacheslav N. Burlayenko; Tomasz Sadowski
Journal: Shock and Vibration
 (2018)

---

## Abstract

A survey of recent developments in the dynamic analysis of sandwich panels with face sheet-to-core debonding is presented. The finite element method within the ABAQUS™ code is utilized. The emphasis is directed to the procedures used to elaborate linear and nonlinear models and to predict dynamic response of the sandwich panels. Recently developed models are presented, which can be applied for structural health monitoring algorithms of real-scale sandwich panels. First, various popular theories of intact sandwich panels are briefly mentioned and a model is proposed to effectively analyse the modal dynamics of debonded and damaged (due to impact) sandwich panels. The influences of debonding size, form, and location and number of such damage incidents on the modal characteristics of sandwich panels are shown. For nonlinear analysis, models based on implicit and explicit time integration schemes are presented and dynamic responses gained with those models are discussed. Finally, questions related to debonding progression at the face sheet-core interface when dynamic loading continues with time are briefly highlighted.


---

## Body

## 1. Introduction

Sandwich panels have long been recognized as one of the efficient structural elements. Due to their intrinsic properties like a high bending stiffness at minimum mass, the capability to be tailored for specific uses, high damping properties, and a great potential for energy absorption, they have found applications in diverse kinds of modern industries and structures. Along with advantages that suggest the sandwich concept, the sandwich panels suffer the consequences of their constructive features. The large differences between thicknesses and elastic moduli of the constitutive layers make the sandwich panels susceptible to debonding of the face sheet from the core at their interface. The sources of debonding usually are a result of imperfections in the manufacturing process, the degassing of the foam core under direct sunlight, the ability for water absorption of cellular types of the core followed by repeated cycles of freezing and thawing at the face sheet-to-core interface, low-velocity impacts, stress concentrations due to localized loading, and so forth. The presence of debonding is often invisible, but it affects the dynamics and the strength of sandwich structures as well as it may become a reason of premature failure below the level of design loads [1]. Therefore, a correct use of sandwich panels in different engineering applications requires a better knowledge of their mechanical behaviour, in particularly, their dynamic response. This knowledge can also provide a basis for nondestructive health monitoring techniques.In modern sandwich panels, the core layer is often made of a soft and flexible material. This leads to the core compressibility and, as a result, to the change in the height of the core during deformation [2]. Thus, to provide a high-fidelity analysis of sandwich panels, the nonlinear deformation patterns in the core should be accounted for in prediction models. This makes the dynamic analysis of the sandwich panels more challenging than that of the laminates. Furthermore, high physical and geometrical mismatching is also an obstacle to handle sandwich panels using the equivalent single layer (ESL) two-dimensional (2D) finite element [3, 4] or dynamic stiffness element [5] models. In this regard, models adopting mixed theories or the layer-wise (LW) approach are preferable for sandwich structures [6]. The latter has given rise to a variety of high-order sandwich panels theories referred to as HSAPTs [7–9], in which the high-order effects in the sandwich panel do not result from some prior assumptions for the displacement field. These theories invite any plate/shell theory formulations for the face sheets and a 3D elasticity theory or equivalent one for the core [10–12].One of the other aspects that relates to the analysis issues of sandwich panels is a microscopically discrete structure of the core, for example, truss cores, honeycomb cores, and corrugated or folded cores of various shapes. The traditional analysis scheme of such sandwich panels is that, first, the core is simplified as an equivalent generally anisotropic homogeneous material and then modelled by using any theory and solution method. Thus, the analysis accuracy strongly depends on correct estimations of equivalent rigidities. Some examples of finite element analyses of sandwich panels with discrete cores can be found, for example, in [13–18].Keeping in mind the challenges in modelling sandwich structures, many efforts have been made to study their vibration responses using both linear and nonlinear models. In doing so, some studies have used the assumptions on geometrical nonlinearity in the dynamics of sandwich panels, for example, [19, 20], while the other ones have examined the influence of material nonlinearity. The effect of a viscoelastic core on damped forced oscillations has been considered, for example, in [21, 22]. The problem of the interlaminar slip between constitutive layers in the nonlinear free vibration has also been discussed, for example, in [23, 24]. Apart from these nonlinear problems, the dynamics of sandwich panels with debonding are another computationally challenging task, even when small displacements and a material linearity are assumed. The difficulties are concerned with modelling physical phenomena arising from debonding. First, the detached surfaces are free of shear and normal stresses and as a result the overall stiffness of the structure is reduced. Second, during loading or oscillations these surfaces may slip longitudinally one with respect to another and/or undergo normal compressive stress if contact between them exists. Such local changes within the debonded region give rise to alterations in the global dynamics of sandwich panels and, also, the stress fields accompanying this dynamic behaviour may result in fracture in the face sheet-to-core interface.To simplify the nature of the problem, the earliest solutions on free vibration were based on the split beam approach under assumptions that the decoupled layers either freely overlap each other or are constrained to move together [25–30]. Later, improved nonlinear models excluding penetration between the layers coming into contact were developed. In [31], to prevent overlapping, the dynamics of a delaminated beam was studied using a piecewise linear virtual spring model, whereas the kinematic contact conditions were established using the node-to-node contact model in [32]. To date, a large volume of the literature on this subject is available [33]. Some recent advanced studies point at the need to use models accounting for coupling between normal and flexural actions in vibrations [34]. As found, it is a reason of parametrically induced vibration of delaminated beams and plates [35–37]. Other studies confirmed the efficiency of the LW-based finite element models compared to those using the ESL approach [38, 39] and showed that such models provide a stress recovery for localized effects [40].Yet, the dynamic finite element analysis (FEA) contributes to developments of vibration control methods of sandwich structures [41, 42]. Alternatively, the results of the dynamic FEA of sandwich panels can be used for increasing the efficiency of the structural health monitoring (SHM). In the latter, natural frequencies, mode shapes, frequency response functions, and time or frequency domain data can be extracted from either linear or nonlinear dynamic FEA. For instance, in [43] natural frequencies extracted from the eigenvalue analysis have been used to detect and locate the saw cut within a composite plate, whereas mode shapes collected from a similar linear analysis have successfully been applied to the delamination location in a composite beam in [44]. Other authors have developed a strain-based damage index based on a linear model for the prediction of the delamination location in composite plates [45]. The curvatures of mode shapes have been utilized for quantifying the damage magnitude in a honeycomb beam in [46]. Also, a linear dynamic analysis is used in FRF-based damage detection techniques [47]. In [48] such approach has successfully been implemented for detecting debonding in a honeycomb beam. For the sake of debonding detection, the linear dynamic FEA has been carried out in [49] to highlight the relative changes of dynamic response between a healthy sandwich plate and a debonded one. More advanced techniques using modal dataset to detect debonding in sandwich structures can be found in some recent works, for example, [50–52].Motivated by the idea to examine the influence of debonding on the modal dynamics of sandwich panels in detail, a number of numerical and experimental studies have been performed during the last two decades. For instance, in [53–55] natural frequencies and corresponding vibration modes have been calculated for flexible debonded sandwich beams. The changes of the modal characteristics in sandwich plates caused by debonding depending on the debonding size and form and of boundary conditions and sandwich core properties have been presented in [56]. Similar researches using linear models to study the dynamic behaviour of sandwich plates that have undergone a postimpact damage or containing a multidebonding have been carried out in [57–59] and [60–62], respectively. A linear model has been also adopted to examine a dynamic stability of a delaminated beam under harmonic longitudinal loading in [63]. An experimental study on the modal dynamics of a delaminated plate has been done in [64]. Changes in the peaks and valleys of the FRFs due to delamination in a honeycomb sandwich beam have been proven by experimental tests in [65]. Effects of the debonding length on fatigue and vibration of sandwich composites have been tested in [66].Because, in general, the dynamics of composite structures are rather nonlinear, the numerical simulations could throw light on the nature of various nonlinear mechanisms, in particular, due to debonding and in turn this knowledge can make the SHM methodologies more reliable [67, 68]. Insight into the real-life dynamics of debonded sandwich panels can be gained by accounting for the “real contact” conditions between the debonded parts. Hence, nonlinear models have been developed to handle the nonlinearities caused by either of or both of the factors, namely, opening and closing of the debonding (breathing) and contact/impact between the debonded layers during vibration. In [69] the contact problem in delaminated surfaces has been modelled in terms of fictitious linear springs to study the dynamic transient behaviour of a delaminated plate. An analytical model of a 1D sandwich beam accounting for the real contact has been created in [70] to examine the influence of the contact phenomenon on the transient behaviour of the beam. In [71–73] an explicit 3D finite element contact model has been developed with ABAQUS for analysing the transient dynamics of an impacted sandwich panel. This model has also been used in [74, 75] to study the effect of the debonding size on the transient dynamic response and the stress state in sandwich plates. Other contact models addressed to discovering the role of interlaminar contact on the nonlinear dynamics of sandwich and laminated panels accounting for geometrical nonlinearity have been proposed in [76–78]. A nonlinear contact model within an implicit approach has been elaborated in [79, 80] to explore the dynamics of sandwich panels with debonding, which are subjected to harmonic forces. The research reported a complex dynamic behaviour of such panels highly dependent on the driving frequency that is inherent in structures with conditions of friction and contact [81]. A potential of such nonlinear models to identify damages in composite structures has been demonstrated in [82].The discussion above clearly demonstrates that high-fidelity dynamic finite element analysis methods of sandwich panels based on linear and nonlinear models are still needed for many engineering problems. It seems that research on the modelling aspects of the dynamics of sandwich panels with debonding is scattered across different models and diverse research areas. Often such studies exploit only one either linear or nonlinear model, and then it is not so easy to juxtapose them within the same problem being considered. Besides, many of the studies are focused on sandwich beams or unidirectional panels. This spatial reduction simplifies the nature of the physical phenomena, because it dictates a specific structural behaviour. To the best of our knowledge, there is a deficiency of 3D models for studying the dynamic behaviour of debonded sandwich plates. Some of them have been developed in our previous papers. The present work can be viewed as a comprehensive survey summarizing the recent developments in 3D finite element modelling of the dynamics of sandwich plates with debonding. Our previous findings in this research area are considered and reexamined under the light of our goal to characterize the dynamic behaviour of sandwich panels with debonding ranging from free vibration to nonlinear oscillations by both linear and nonlinear 3D finite element models. So, the present work takes a step towards the classification of computational models for quantitatively estimating the dynamics of sandwich panels with debonding and understanding the responses of such structures based on the numerical examples.


## 2. Outlines of Dynamic FEAs

In this section, basic aspects of the finite element procedures used in the current work to perform dynamic finite element analyses (FEAs) are briefly presented.


### 2.1. Equations of Motion

The finite element analysis of the system of elements provides the solution of the equations of motion corresponding to the nodal degrees of freedom, which, in general, can be written as follows [83]:(1)FtI+FtD+Ftint=Ftext,where t is a certain moment of time and FI, FD, Fint, and Fext stand for inertia, damping, and internally and externally applied forces, respectively. In the case of a linear problem, (1) reads as(2)MU¨+CU˙+KU=Ftext,where M, C, and K are, respectively, mass, damping, and stiffness matrices of the element assemblage and U are nodal points displacements, and U˙ and U¨ are their time derivatives referring to nodal velocities and acceleration, respectively. The matrices of the physical properties of the finite elements are related to the force vectors as FtI=MU¨, FtD=CU˙, and Ftint=KU. Under assumptions that initial conditions and boundary conditions applied to the system remain unchanged during loading, the solution of (2) completely describes a linear response of the system.In presence of contact between elements of the system, the system response is linear, if other nonlinearities except contact do not exist, only prior to the contact conditions met. Otherwise changes in boundary conditions at a certain load level occur. Then, (2) is to be supplemented by contact constraints imposed at a set of all prospective contact pairs over the loading time as follows:(3)gN≥0,tN≤0,tNgN=0,tNg˙N=0and the constraints appropriate for the prescription of a Coulomb friction law, for example, as [84](4)ΨT=tT-μtN≤0,LυgT=χtTtT,χ≥0,χΨT=0.In (3) and (4), gN and gT stand for normal and tangential “gap” functions defined for all active contact pairs, tN and tT are normal and tangential contact traction acting at those pairs, Lυ is a Lie derivative (a convective derivative that defines a frame invariant measure of rate of gT in the case of slipping), μ is the coefficient of friction, and χ is a nonnegative scalar. For conciseness, we skip details and refer to works in this area, for example, [84]. The contact traction contributes to equilibrium of the system as contact force vector FC acting on some nodal points at a certain moment of timet; that is, (2) takes the form(5)MU¨+CU˙+KU+FCU=Ftext.


### 2.2. Eigenvalue Problem

Dynamic characteristics of a linear structural system (2) are governed by the natural frequencies and the corresponding mode shapes. The determination of them requires the solution of an eigenvalue problem:(6)Kϕ=ω2Mϕ.Hereinω is an undamped circular frequency and ϕ is the corresponding vibration mode. Then,n uncoupled solutions of (6) can be composed in the form convenient for frequency extraction:(7)KΦ=MΦΩ2,where in the matrix Φ the columns are orthogonalized mass-normalized eigenvectors ϕ-1,…,ϕ-n and Ω2=diagω12⋯ωn2 is a diagonal matrix containing squared eigenfrequencies.


### 2.3. Mode-Based Dynamic Analysis

To provide predictions of the steady-state dynamic response of a linear system (2) due to harmonic excitation Ft=FeiΩ-t at a given driving frequency Ω-, a modal superposition method can be used. Then, in the modal subspace (2) takes the form(8)M~Q¨+C~Q˙+K~Q=ΦTFt.Here the nodal displacements U=ΦQ are expressed in terms of the vector of generalized displacements Q and the normal modes Φ being found from the eigenvalue analysis; the matrices M~, C~, and K~ are the projections of M, C, and K onto the modal subspace; that is, M~=ΦTMΦ and C~=ΦTCΦ and K~=ΦTKΦ and ΦTFt=F~eiΩ-t. In the case of linear analysis M~ is the identity matrix, K~=diagωi2 and in case of no damping coupling C~=diag2ωiξi, where ξi are modal damping rations. Hence, the equations in system (8) with respect to the generalized displacements are decoupled and for therth mode have the form(9)Qr¨+diag2ωrξrQr˙+diagωr2Qr=F~reiΩ-t.Once all Q are evaluated, the physical response of the original system in terms of nodal displacements, velocities, acceleration, and stresses can be recovered. The analysis is done as a frequency sweep by applying the loading at a series of different frequencies and recording the system response [85].


### 2.4. Direct Integration Analyses

The equations in system (5) cannot be solved independently; a direct integration over time is required. One would normally resort to numerical approximation procedures to do it. The explicit central difference time integrator used in ABAQUS/Explicit [85] presents (5) at a typical time increment Δt in the form(10)U¨t=M~-1Ftext-FDU˙t-FintUt+FCUt,where M~ is a lumped mass matrix obtained by transform of the consistent mass matrix and the velocities and the displacements are calculated at each time increment according to the expressions(11)U˙t+1/2Δt=U˙t-1/2Δt+ΔtU¨t,Ut+Δt=Ut+ΔtU˙t+1/2Δt.This explicit time integration scheme is conditionally stable only. An approximation to its stability limit can be achieved by evaluating the highest modal natural frequency ωmax in the FE mesh as Δt≤2/ωmax [83].The implicit Hilber-Hughes-Taylor (HHT) temporal integrator is used by ABAQUS/Standard [85]. This integration scheme applied to (5) yields(12)MU¨t+Δt+CU˙t+Δt+KUt+αΔt+FCUt+αΔt=Ft+αΔtext,where(13)Ut+αΔt=αUt+Δt+1-αUt,Ut+Δt=Ut+ΔtU˙t+Δt221-2βU¨t+2βU¨t+Δt,U˙t+Δt=U˙t+Δt1-γU¨t+γU¨t+Δt.The unconditional stability occurs for -1/3≤α≤0,  β=1/41-α2, and γ=1/2-α [83]. The linearized problem in (12) within thejth Newton-Raphson iteration has the form(14)∂∂UMU¨+CU˙+KU+FCUUt+αΔtjΔU=Ft+αΔtext-MU¨t+Δtj-CU˙t+Δtj-KUt+αΔtj-FCUt+αΔtjwith updating unknowns as Ut+αΔtj+1=Ut+αΔtj+ΔU.


### 2.5. Modelling Workflow

In Figure1 the modelling workflow followed in this work is summarized. The finite element analyses carried out with ABAQUS are presented in the order of their computational cost and modelling complexity, that is, the computational procedures mentioned above and used in those analyses and some approaches and features applied to the development of linear and nonlinear models. Thereby, for sandwich panels, first, the frequency eigenvalue analysis (6) and the modal-based harmonic analysis (8) are performed using linear models and, then, the general nonlinear dynamics under an impulse load and harmonic loading are worked out using explicit (9) and implicit (12) time integration schemes. Both a debonded sandwich panel and the same healthy one are simulated in the calculations. By tracking the differences between their dynamic responses, the effects of debonding on the linear and nonlinear dynamic behaviour of the sandwich panel are evaluated. Moreover, using both linear and nonlinear models for the same sandwich panel, the impact for each of them for highlighting dynamic effects can be clearly seen.Figure 1
Workflow followed in this work.


## 2.1. Equations of Motion

The finite element analysis of the system of elements provides the solution of the equations of motion corresponding to the nodal degrees of freedom, which, in general, can be written as follows [83]:(1)FtI+FtD+Ftint=Ftext,where t is a certain moment of time and FI, FD, Fint, and Fext stand for inertia, damping, and internally and externally applied forces, respectively. In the case of a linear problem, (1) reads as(2)MU¨+CU˙+KU=Ftext,where M, C, and K are, respectively, mass, damping, and stiffness matrices of the element assemblage and U are nodal points displacements, and U˙ and U¨ are their time derivatives referring to nodal velocities and acceleration, respectively. The matrices of the physical properties of the finite elements are related to the force vectors as FtI=MU¨, FtD=CU˙, and Ftint=KU. Under assumptions that initial conditions and boundary conditions applied to the system remain unchanged during loading, the solution of (2) completely describes a linear response of the system.In presence of contact between elements of the system, the system response is linear, if other nonlinearities except contact do not exist, only prior to the contact conditions met. Otherwise changes in boundary conditions at a certain load level occur. Then, (2) is to be supplemented by contact constraints imposed at a set of all prospective contact pairs over the loading time as follows:(3)gN≥0,tN≤0,tNgN=0,tNg˙N=0and the constraints appropriate for the prescription of a Coulomb friction law, for example, as [84](4)ΨT=tT-μtN≤0,LυgT=χtTtT,χ≥0,χΨT=0.In (3) and (4), gN and gT stand for normal and tangential “gap” functions defined for all active contact pairs, tN and tT are normal and tangential contact traction acting at those pairs, Lυ is a Lie derivative (a convective derivative that defines a frame invariant measure of rate of gT in the case of slipping), μ is the coefficient of friction, and χ is a nonnegative scalar. For conciseness, we skip details and refer to works in this area, for example, [84]. The contact traction contributes to equilibrium of the system as contact force vector FC acting on some nodal points at a certain moment of timet; that is, (2) takes the form(5)MU¨+CU˙+KU+FCU=Ftext.


## 2.2. Eigenvalue Problem

Dynamic characteristics of a linear structural system (2) are governed by the natural frequencies and the corresponding mode shapes. The determination of them requires the solution of an eigenvalue problem:(6)Kϕ=ω2Mϕ.Hereinω is an undamped circular frequency and ϕ is the corresponding vibration mode. Then,n uncoupled solutions of (6) can be composed in the form convenient for frequency extraction:(7)KΦ=MΦΩ2,where in the matrix Φ the columns are orthogonalized mass-normalized eigenvectors ϕ-1,…,ϕ-n and Ω2=diagω12⋯ωn2 is a diagonal matrix containing squared eigenfrequencies.


## 2.3. Mode-Based Dynamic Analysis

To provide predictions of the steady-state dynamic response of a linear system (2) due to harmonic excitation Ft=FeiΩ-t at a given driving frequency Ω-, a modal superposition method can be used. Then, in the modal subspace (2) takes the form(8)M~Q¨+C~Q˙+K~Q=ΦTFt.Here the nodal displacements U=ΦQ are expressed in terms of the vector of generalized displacements Q and the normal modes Φ being found from the eigenvalue analysis; the matrices M~, C~, and K~ are the projections of M, C, and K onto the modal subspace; that is, M~=ΦTMΦ and C~=ΦTCΦ and K~=ΦTKΦ and ΦTFt=F~eiΩ-t. In the case of linear analysis M~ is the identity matrix, K~=diagωi2 and in case of no damping coupling C~=diag2ωiξi, where ξi are modal damping rations. Hence, the equations in system (8) with respect to the generalized displacements are decoupled and for therth mode have the form(9)Qr¨+diag2ωrξrQr˙+diagωr2Qr=F~reiΩ-t.Once all Q are evaluated, the physical response of the original system in terms of nodal displacements, velocities, acceleration, and stresses can be recovered. The analysis is done as a frequency sweep by applying the loading at a series of different frequencies and recording the system response [85].


## 2.4. Direct Integration Analyses

The equations in system (5) cannot be solved independently; a direct integration over time is required. One would normally resort to numerical approximation procedures to do it. The explicit central difference time integrator used in ABAQUS/Explicit [85] presents (5) at a typical time increment Δt in the form(10)U¨t=M~-1Ftext-FDU˙t-FintUt+FCUt,where M~ is a lumped mass matrix obtained by transform of the consistent mass matrix and the velocities and the displacements are calculated at each time increment according to the expressions(11)U˙t+1/2Δt=U˙t-1/2Δt+ΔtU¨t,Ut+Δt=Ut+ΔtU˙t+1/2Δt.This explicit time integration scheme is conditionally stable only. An approximation to its stability limit can be achieved by evaluating the highest modal natural frequency ωmax in the FE mesh as Δt≤2/ωmax [83].The implicit Hilber-Hughes-Taylor (HHT) temporal integrator is used by ABAQUS/Standard [85]. This integration scheme applied to (5) yields(12)MU¨t+Δt+CU˙t+Δt+KUt+αΔt+FCUt+αΔt=Ft+αΔtext,where(13)Ut+αΔt=αUt+Δt+1-αUt,Ut+Δt=Ut+ΔtU˙t+Δt221-2βU¨t+2βU¨t+Δt,U˙t+Δt=U˙t+Δt1-γU¨t+γU¨t+Δt.The unconditional stability occurs for -1/3≤α≤0,  β=1/41-α2, and γ=1/2-α [83]. The linearized problem in (12) within thejth Newton-Raphson iteration has the form(14)∂∂UMU¨+CU˙+KU+FCUUt+αΔtjΔU=Ft+αΔtext-MU¨t+Δtj-CU˙t+Δtj-KUt+αΔtj-FCUt+αΔtjwith updating unknowns as Ut+αΔtj+1=Ut+αΔtj+ΔU.


## 2.5. Modelling Workflow

In Figure1 the modelling workflow followed in this work is summarized. The finite element analyses carried out with ABAQUS are presented in the order of their computational cost and modelling complexity, that is, the computational procedures mentioned above and used in those analyses and some approaches and features applied to the development of linear and nonlinear models. Thereby, for sandwich panels, first, the frequency eigenvalue analysis (6) and the modal-based harmonic analysis (8) are performed using linear models and, then, the general nonlinear dynamics under an impulse load and harmonic loading are worked out using explicit (9) and implicit (12) time integration schemes. Both a debonded sandwich panel and the same healthy one are simulated in the calculations. By tracking the differences between their dynamic responses, the effects of debonding on the linear and nonlinear dynamic behaviour of the sandwich panel are evaluated. Moreover, using both linear and nonlinear models for the same sandwich panel, the impact for each of them for highlighting dynamic effects can be clearly seen.Figure 1
Workflow followed in this work.


## 3. Aspects of Model Discretization

The first step in a finite element analysis is the discretization process, in which the actual continuous structure is idealized as an assemblage of discrete finite elements interconnected at node. The number and type of finite elements selected for this, in essence, depend on the physical complexity of the structure and the required accuracy in the displacements and stresses [83]. In accordance with the discussion above, various theories and, as a result, modelling techniques for the FEA of sandwich panels are possible. A summary of sandwich finite element models that are available for modelling, in particular with the ABAQUS code, is collected in Table 1.Table 1
Summary of techniques on finite element modelling of sandwich panels [11, 12].
No.
Modelling schemes
Comments
(1)
Homogenous shell elements based on either CPT1 or FSDT2
Simple efficient for general analysis. Transverse compressibility is neglected.
(2)
Layered shell elements with shear deformation
Computationally efficient. Transverse compressibility is ignored.
(3)
Plane stress/plane strain elements modelling sandwich cross-section
Applicable to 2D and axisymmetric problems only.
(4)
Plane stress/plane strain elements for the core and beam elements for skins
Skin is idealized by Euler-Bernoulli beam element.
(5)
Solid-brick elements for the core and conventional shell elements for skins
Compressibility of the core is accounted for. Conventional shell elements may suffer from shear-locking.
(6)
3D model with solid-brick elements for both the core and the skins
Computationally very expensive model.
(7)
3D model with layered solid cross-section
No restrictions on the core and skins kinematics.
(8)
Layer-wise models combining 3D solid and continuum shell elements
Core compressibility can be accounted for.
1Classical plate theory. 2First-order shear deformation theory.


### 3.1. Layer-Wise Based Model

To tackle the problem at hand, the model handling the core flexibility and the transverse shear variation and the geometry and location of an embedded debonded zone has been tailored in ABAQUS using the layer-wise shell-solid approach (number 5 in Table1). This modelling strategy has been used and tested for modelling various static and dynamic problems related to sandwich panels in our previous works, for example, [10, 11, 13, 56]. In the present paper this model is revised and unified for running both linear and nonlinear analyses by switching appropriate options and/or assigning new features and properties available in the preprocessing stage and solvers of ABAQUS [85]. The main emphasis of this finite element model is that it realizes the layer-wise HSAPT-like approach for the sandwich structure establishing no restrictions on the core kinematics and allowing independent moving of top and bottom skins with respect to each other. In this regard, the upper and lower skins are discretized with the eight-node 24-DOF quadrilateral continuum shell elements obeying the assumptions of the FSDT theory within the ESL approach, while the core is discretized with the eight-node 20-DOF hexahedral (brick) elements [85] as shown in Figure 2. The degrees of freedom (DOFs) of the brick and continuum shell elements are conveniently coupled with each other, so no additional compatibility conditions at the interface between skins and core is needed. Hence, the final governing equations of the sandwich structure assemble the displacements variables of the upper and lower skins, and the core automatically ensures the continuity of displacements at the interfaces between them. In the case of a heterogeneous core material (honeycomb), the core was replaced by an equivalent model of anisotropic material.Figure 2
Aspects of finite element discretization: (a) a 3D model of sandwich plate; (b) continuous shell finite element; and (c) solid-brick finite element.
(a)
(b)
(c)A sandwich panel that has one or more arbitrary debonded regions at one or two of its skin-to-core interfaces is geometrically idealized using the multipartition discretization. The 3D mesh of the debonded sandwich plate consists of four different parts connected with each other through the shared nodes. These parts are discretized by applying different mesh methods and mesh densities, Figure3. That is, a fine mesh is used for the debonded zone, the next zone surrounding the debonded one has a gradually decreasing mesh density, and a coarse mesh is introduced for the core to minimize a CPU time in calculations. In the model, the mesh size is suitable to be refined to get a required convergence. Debonding is modelled by creating a small gap (about 1% of the face sheet thickness) between the face sheet and the core. No artificial adjustment of either material properties or geometrical entities is made within the debonded region to ensure a physically real case as close as possible.Figure 3
A 3D finite element model of sandwich plate with a penny-shaped debonded region [10].As seen in Figure3, the finite element mode consists of two types of interfacial regions, namely, a “fully bonded zone” and a “debonded zone.” In general, the detached surfaces within the debonded region may contact each other; that is, they can slip one with respect to the other in any longitudinal direction or can come into contact/impact vertically either entirely or partially involving only some of their zones. Contact models developed using the ABAQUS code to handle the behaviour between the detached skin and core in linear and nonlinear dynamic FEAs are discussed below.


### 3.2. Models with Spring Element

In the modal dynamics of debonded sandwich panels, a spring element “SPRING2” is introduced between the double nodes within the debonded area (Figure4(a)). The stiffness of the spring element is assumed to be zero in the case of tension and to be a large value in compression, that is, when the relative transverse displacement between the face sheet and the core Δugoes to zero. Since any inelastic effects are deactivated for modal analyses [85], we approximate the behaviour of the spring element by its two discrete states as shown in Figure 4(b). In such way we are able easily to activate one of the two “contact” options in the modal dynamics. That is, setting a zero stiffness to the spring elements, the contact behaviour of the detached parts is described by the free delamination model; that is, the interfaces move freely, while a nonzero stiffness value of those elements realizes the constrained delamination model between them that restrains the face sheet and the core to move together. So, a contact problem between the face sheet and core in linear dynamic FEAs is reduced to using either “mode constrained” or “mode-free” delamination model [26–29]. It should be noticed that although the spring elements with given nonlinear constitutive laws enable modelling contact in a general dynamic analysis, they are less convenient for handling the debonding problem than the surface contact definition used in this work and considered below.Figure 4
A detailed view of the linear model sandwich plate with debonding [56]: (a) zoomed debonded zone and (b) constitutive behaviour of the spring element.
(a)
(b)


### 3.3. Models with “Real Contact”

A general dynamic analysis enables representing “real contact” and as a result a nonlinear finite element model is addressed to considerations. In this case, between the surfaces of parts (I) and (II) within the debonded zone (Figure3) contact definitions are imposed. The surface-to-surface contact formulation in terms of slave and master surfaces is used. The relative motion of the contacting surfaces in the contact simulation is described with small sliding kinematics that allow their separation and sliding, as shown in Figure 5. This is a reasonable assumption because of small oscillations in the debonded zone, but it is more computationally effective than a finite formulation that is especially important in the case of 3D simulations.Figure 5
Spatial behaviour of contacting surfaces: (a) full contact or closing; (b) partial contact; and (c) no contact or opening.
(a)
(b)
(c)The constitutive behaviour of the surfaces coming into contact in the normal direction is assumed to be governed by the “hard contact” model. The model implies that the interacting surfaces transmit no contact pressure unless the nodes of the slave surface contact the master surface and no penetration is allowed at each constraint location. In the case of friction, an isotropic Coulomb friction model is adopted to model frictional stick and slip phenomena. To resolve the contact constraints imposed on the mentioned surfaces during the dynamic analysis, appropriate contact algorithms are run in calculations. The kinematic predictor/corrector algorithm is used in the explicit dynamic FEA in the case of impulse loading, while the penalty algorithm processes contact interactions for the implicit dynamic FEA when harmonic excitations are assumed to occur. The former algorithm strictly enforces contact constraints and allows explicit updating contact forces in the calculations, Figure6(a). The latter scheme has a weaker enforcement of contact constraints, Figure 6(b), but it is less computationally expensive than other methods available in ABAQUS, for example, the method of Lagrange multipliers. A detailed description of these contact schemes can be found, for example, in [74, 75, 79, 80]. Furthermore, the finite element models accounting for “real contact” can recover stress distributions caused by both contact/impact actions and deformations. Hence, these models can provide a strength assessment of the face sheet-to-core interface and, in essence, are able to estimate a possibility of separation between the core and the face sheet along their interface, where a debonding failure often initiates. However, to track the debonding progression, additional modelling efforts are needed, as shown later.Figure 6
Schemes of contact algorithms [85]: (a) kinematic predictor/corrector; and (b) penalty enforcement.
(a)
(b)


### 3.4. Models of Impacted Sandwich Panels

Low-velocity impacts induce barely visible impact damage (BVID), which may be one of the reasons for debonding occurrence in sandwich panels [86]. Therefore, the modelling of the dynamics of sandwich panels with impact-induced damage directly falls into the objectives of this work. According to the experimental findings, a sandwich plate damaged by low-velocity impact has a combination of failure modes among which damage of the face sheet, core crushing, and face sheet-to-core debonding are primary. Figure 7(a) shows the key geometrical parameters of the representative cross-section of the sandwich specimen assumed to be impacted by a spherical object: the peak depth of the residual face sheet indentation, δdent, the peak depth associated with core crushing, δcr, the cavity (debonding), δcav, developed between the face sheet and the core as a consequence of differences in their indentation depths and the radii of the planar dimension of the impacted face sheet, Rdent, and the crushed core, Rcr [86].Figure 7
Impact-damaged region: (a) key geometrical parameters [86]; and (b) finite element discretization.
(a)
(b)The modelling strategy of sandwich panels which had been damaged by impact is the same as described above, except the postimpact zone itself. The damage imparted into the face sheet and the core because of an impact event is simulated by reducing elastic properties of the finite elements in the damaged regions. Appropriate reduction coefficients are used for this purpose. For the crushed core, the stiffness of elements is reduced gradually around the damaged site, as shown by a colour gradient in Figure7(b). The residual indentations of the impacted face sheet and the crushed core are modelled as curvilinear surfaces, and the cavity between them is presented as well. The spring elements and the contact models mentioned above are introduced between appropriate surfaces of the impacted zone in the case of linear and nonlinear dynamic analyses, respectively.


## 3.1. Layer-Wise Based Model

To tackle the problem at hand, the model handling the core flexibility and the transverse shear variation and the geometry and location of an embedded debonded zone has been tailored in ABAQUS using the layer-wise shell-solid approach (number 5 in Table1). This modelling strategy has been used and tested for modelling various static and dynamic problems related to sandwich panels in our previous works, for example, [10, 11, 13, 56]. In the present paper this model is revised and unified for running both linear and nonlinear analyses by switching appropriate options and/or assigning new features and properties available in the preprocessing stage and solvers of ABAQUS [85]. The main emphasis of this finite element model is that it realizes the layer-wise HSAPT-like approach for the sandwich structure establishing no restrictions on the core kinematics and allowing independent moving of top and bottom skins with respect to each other. In this regard, the upper and lower skins are discretized with the eight-node 24-DOF quadrilateral continuum shell elements obeying the assumptions of the FSDT theory within the ESL approach, while the core is discretized with the eight-node 20-DOF hexahedral (brick) elements [85] as shown in Figure 2. The degrees of freedom (DOFs) of the brick and continuum shell elements are conveniently coupled with each other, so no additional compatibility conditions at the interface between skins and core is needed. Hence, the final governing equations of the sandwich structure assemble the displacements variables of the upper and lower skins, and the core automatically ensures the continuity of displacements at the interfaces between them. In the case of a heterogeneous core material (honeycomb), the core was replaced by an equivalent model of anisotropic material.Figure 2
Aspects of finite element discretization: (a) a 3D model of sandwich plate; (b) continuous shell finite element; and (c) solid-brick finite element.
(a)
(b)
(c)A sandwich panel that has one or more arbitrary debonded regions at one or two of its skin-to-core interfaces is geometrically idealized using the multipartition discretization. The 3D mesh of the debonded sandwich plate consists of four different parts connected with each other through the shared nodes. These parts are discretized by applying different mesh methods and mesh densities, Figure3. That is, a fine mesh is used for the debonded zone, the next zone surrounding the debonded one has a gradually decreasing mesh density, and a coarse mesh is introduced for the core to minimize a CPU time in calculations. In the model, the mesh size is suitable to be refined to get a required convergence. Debonding is modelled by creating a small gap (about 1% of the face sheet thickness) between the face sheet and the core. No artificial adjustment of either material properties or geometrical entities is made within the debonded region to ensure a physically real case as close as possible.Figure 3
A 3D finite element model of sandwich plate with a penny-shaped debonded region [10].As seen in Figure3, the finite element mode consists of two types of interfacial regions, namely, a “fully bonded zone” and a “debonded zone.” In general, the detached surfaces within the debonded region may contact each other; that is, they can slip one with respect to the other in any longitudinal direction or can come into contact/impact vertically either entirely or partially involving only some of their zones. Contact models developed using the ABAQUS code to handle the behaviour between the detached skin and core in linear and nonlinear dynamic FEAs are discussed below.


## 3.2. Models with Spring Element

In the modal dynamics of debonded sandwich panels, a spring element “SPRING2” is introduced between the double nodes within the debonded area (Figure4(a)). The stiffness of the spring element is assumed to be zero in the case of tension and to be a large value in compression, that is, when the relative transverse displacement between the face sheet and the core Δugoes to zero. Since any inelastic effects are deactivated for modal analyses [85], we approximate the behaviour of the spring element by its two discrete states as shown in Figure 4(b). In such way we are able easily to activate one of the two “contact” options in the modal dynamics. That is, setting a zero stiffness to the spring elements, the contact behaviour of the detached parts is described by the free delamination model; that is, the interfaces move freely, while a nonzero stiffness value of those elements realizes the constrained delamination model between them that restrains the face sheet and the core to move together. So, a contact problem between the face sheet and core in linear dynamic FEAs is reduced to using either “mode constrained” or “mode-free” delamination model [26–29]. It should be noticed that although the spring elements with given nonlinear constitutive laws enable modelling contact in a general dynamic analysis, they are less convenient for handling the debonding problem than the surface contact definition used in this work and considered below.Figure 4
A detailed view of the linear model sandwich plate with debonding [56]: (a) zoomed debonded zone and (b) constitutive behaviour of the spring element.
(a)
(b)


## 3.3. Models with “Real Contact”

A general dynamic analysis enables representing “real contact” and as a result a nonlinear finite element model is addressed to considerations. In this case, between the surfaces of parts (I) and (II) within the debonded zone (Figure3) contact definitions are imposed. The surface-to-surface contact formulation in terms of slave and master surfaces is used. The relative motion of the contacting surfaces in the contact simulation is described with small sliding kinematics that allow their separation and sliding, as shown in Figure 5. This is a reasonable assumption because of small oscillations in the debonded zone, but it is more computationally effective than a finite formulation that is especially important in the case of 3D simulations.Figure 5
Spatial behaviour of contacting surfaces: (a) full contact or closing; (b) partial contact; and (c) no contact or opening.
(a)
(b)
(c)The constitutive behaviour of the surfaces coming into contact in the normal direction is assumed to be governed by the “hard contact” model. The model implies that the interacting surfaces transmit no contact pressure unless the nodes of the slave surface contact the master surface and no penetration is allowed at each constraint location. In the case of friction, an isotropic Coulomb friction model is adopted to model frictional stick and slip phenomena. To resolve the contact constraints imposed on the mentioned surfaces during the dynamic analysis, appropriate contact algorithms are run in calculations. The kinematic predictor/corrector algorithm is used in the explicit dynamic FEA in the case of impulse loading, while the penalty algorithm processes contact interactions for the implicit dynamic FEA when harmonic excitations are assumed to occur. The former algorithm strictly enforces contact constraints and allows explicit updating contact forces in the calculations, Figure6(a). The latter scheme has a weaker enforcement of contact constraints, Figure 6(b), but it is less computationally expensive than other methods available in ABAQUS, for example, the method of Lagrange multipliers. A detailed description of these contact schemes can be found, for example, in [74, 75, 79, 80]. Furthermore, the finite element models accounting for “real contact” can recover stress distributions caused by both contact/impact actions and deformations. Hence, these models can provide a strength assessment of the face sheet-to-core interface and, in essence, are able to estimate a possibility of separation between the core and the face sheet along their interface, where a debonding failure often initiates. However, to track the debonding progression, additional modelling efforts are needed, as shown later.Figure 6
Schemes of contact algorithms [85]: (a) kinematic predictor/corrector; and (b) penalty enforcement.
(a)
(b)


## 3.4. Models of Impacted Sandwich Panels

Low-velocity impacts induce barely visible impact damage (BVID), which may be one of the reasons for debonding occurrence in sandwich panels [86]. Therefore, the modelling of the dynamics of sandwich panels with impact-induced damage directly falls into the objectives of this work. According to the experimental findings, a sandwich plate damaged by low-velocity impact has a combination of failure modes among which damage of the face sheet, core crushing, and face sheet-to-core debonding are primary. Figure 7(a) shows the key geometrical parameters of the representative cross-section of the sandwich specimen assumed to be impacted by a spherical object: the peak depth of the residual face sheet indentation, δdent, the peak depth associated with core crushing, δcr, the cavity (debonding), δcav, developed between the face sheet and the core as a consequence of differences in their indentation depths and the radii of the planar dimension of the impacted face sheet, Rdent, and the crushed core, Rcr [86].Figure 7
Impact-damaged region: (a) key geometrical parameters [86]; and (b) finite element discretization.
(a)
(b)The modelling strategy of sandwich panels which had been damaged by impact is the same as described above, except the postimpact zone itself. The damage imparted into the face sheet and the core because of an impact event is simulated by reducing elastic properties of the finite elements in the damaged regions. Appropriate reduction coefficients are used for this purpose. For the crushed core, the stiffness of elements is reduced gradually around the damaged site, as shown by a colour gradient in Figure7(b). The residual indentations of the impacted face sheet and the crushed core are modelled as curvilinear surfaces, and the cavity between them is presented as well. The spring elements and the contact models mentioned above are introduced between appropriate surfaces of the impacted zone in the case of linear and nonlinear dynamic analyses, respectively.


## 4. Numerical Studies

An overview of the use of the linear and nonlinear finite elements models considered above concerning the dynamics of sandwich panels is presented in this section. In this regard, some computational examples taken from our previous works are revised and the main results are briefly discussed below.


### 4.1. Linear Dynamic Analyses

Example 1.
To make evident the effect of “soft” core on the free vibration kinematics of the sandwich panel, the free vibration analysis of a fully bonded sandwich plate is considered first. A squared simply supported plate with a sidea=200 mm and the total thickness h=20 mm has been analysed. The sandwich plate layout corresponds to a five-layer structure 0/90/core/0/90, where face sheets are made of a composite anisotropic material, while the core is isotropic. The ratio between the face sheet and core thicknesses is 10, whereas the mismatching in material properties is the order of 103 as shown in Table 2.Table 2
Material properties of sandwich panels used in the calculations.
Examples
Ref.
Constituents
Example1
[10]
Face sheets,ρf = 1627 kgm−3
E
1, GPa
E
2 = E3, GPa
G
12 = G13, GPa
G
23, GPa
ν
12 = ν13
ν
23
131.0
10.34
6.895
6.205
0.22
0.49
Core,ρc = 97 kgm−3
E
1 = E2, GPa
E
3, GPa
G
12 = G13, GPa
G
23, GPa
ν
12
ν
13 = ν23
6.890  ×  10−3
6.890  ×  10−3
3.450  ×  10−3
3.450  ×  10−3
0.0
0.0
Example2
[53]
Face sheets
ρ
f = 4400 kgm−3, Ef = 36.0 GPa, ν = 0.3
Core,ρc = 52 kgm−3
E
c = 50.0 MPa, Gc = 21.0 MPa
Examples3–5
[73]
Face sheets,ρf=1650 kgm−3
E
1 = E3, GPa
E
2, GPa
G
12 = G23, GPa
G
13, GPa
ν
12
ν
13 = ν23
19.3
3.48
1.65
7.70
0.05
0.25
Core,ρc = 52 kgm−3, Ec = 85.0 MPa, Gc = 30.0 MPa, ν = 0.42
Examples6–9
[56]
Face sheets,ρf = 1650 kgm−3
E
1, GPa
E
2 = E3, GPa
G
12 = G13, GPa
G
23, GPa
ν
12
ν
13 = ν23
140
10
4.6
3.8
0.25
0.25
Honeycomb core,ρc = 57.17 kgm−3
E
1 = E2, GPa
E
3, GPa
G
12, GPa
G
13, GPa
G
23, GPa
ν
12=ν13=ν23
0.461 × 10−3
1.464
0.194 × 10−3
341.7 × 10−3
192.1 × 10−3
0.31
ρ
c
H
100 = 100 kgm−3, Ec = 105.0 MPa, Gc = 45.0 MPa, ν= 0.32
ρ
c
H
200= 200 kgm−3, Ec = 230.0 MPa, Gc = 85.0 MPa, ν = 0.33In the finite element mesh, one continuum shell element with the two-layer cross-section scheme across the thickness of the face sheets and five solid elements across the thickness of the core are used. In plane, the plate is discretized by 50 elements along each of its edges. Apart from the bending modes, the sandwich plate with flexible core exhibits in-plane and shear (twisting) and symmetric out-of-plane (pumping) modes, Figure8. For the latter modes, the core is removed from the corresponding contour plots for the better presentation of this phenomenon. The pumping modes indicate that the two face sheets move vertically in the opposite direction that relates to the changes of the plate height. Therefore, the layer-wise shell-solid finite element model proposed for modelling sandwich plates with a flexible core is able to capture such nonclassical effects.Figure 8
Contour plots of the mode shapes of sandwich plate with flexible core: (a) an in-plane mode; (b) a shear mode; (c) and (d) pumping modes [10].
(a)
(b)
(c)
(d)Example 2.
Firstly, the linear model with spring elements has been verified by comparing the natural frequencies and corresponding mode shapes of a debonded sandwich beam with those that are known in the literature [53]. A sandwich beam cored by foam with rectangular cross-section contains a deboned zone at its middle span, Figure 9. Material properties of the beam are listed in Table 2.Figure 9
A debonded sandwich beam with foam core [53].A comparison of first six natural frequencies of the debonded sandwich beam is presented in Table3. The mode shapes associated with the appropriate natural frequencies are shown in Figure 10.Table 3
Comparison of natural frequencies (Hz) of intact and debonded foam cored sandwich beams.
Mode
Intact
Debonded
[53]
[55]
Present
[53]
[55]
Present
(1)
289.3
293.52
293.46
288.98
293.52
293.07
(2)
683.3
722.96
707.09
388.32
360.81
433.67
(3)
1096.9
1139.51
1106.7
1093.2
1139.31
1093.2
(4)
1151.6
1545.91
1495.8
1146.9
1146.32
1132.0
(5)
1778.2
-
1818.7
1771.3
1769.9
1769.9
(6)
1895.3
1863.62
1918.1
1842.2
1948.35
2080.2Figure 10
Mode shapes of the debonded sandwich beam corresponding to natural frequencies: (a) 293.07 Hz, (b) 433.67 Hz, (c) 1093.2 Hz, (d) 1132.0 Hz, (e) 1769.9 Hz, and (f) 2080.2 Hz.
(a)
(b)
(c)
(d)
(e)
(f)It follows from the analysis that debonding reduces the natural frequencies compared to those of the same intact beam and this reduction is different for each mode. Apart from the global modal forms, the mode shapes of the debonded beam demonstrate both local and coupled shapes of vibration as well.Example 3.
The free vibration analysis of a sandwich plate containing a penny-shaped debonded zone is further demonstrated. The plate with a 50 mm thick foam core and each of the 2.4 mm thick face sheets is simply supported. The planar dimensions of the plate are of 270 mm by 180 mm. The radius of the debonded zone,rd, is of 20 mm; that is, the damage parameter defined as D%=πrd2/Atotal is about 5% of the plate area, Atotal. The material properties of each plate layer are presented in Table 2.The influence of debonding on the modal parameters of the sandwich plate reveals, first, the change of the order of frequencies and the existence of new ones in the spectrum. This happens due to local or coupled modes relating to debonding. Second, the reduction of the low natural frequencies occurs, especially, for bending modes. Herewith, the decreasing strongly depends on the mode shape. In the case of the natural frequencies associated with in-plane modes, there are practically no changes, and the changes are small for shear modes. Finally, the debonding affects the high frequencies in another manner than the lower ones, but there is no defined trend as a mode number increases. In the high part of the spectrum, there exist frequencies which are even higher than those of the intact plate; that is, the thickening phenomenon occurs [64]. The shift in natural frequencies as a result of comparing the healthy plate with the debonded one is presented in Figure 11, and several local mode shapes of the debonded sandwich plate are illustrated in Figure 12.Figure 11
The shift in natural frequencies of the debonded sandwich plate.Figure 12
Local mode shapes in the debonded sandwich plate: (a)(1,1); (b) (1,2); and (c) (1,3).
(a)
(b)
(c)Example 4.
In the next step, the mode-base steady-state dynamic analysis of the intact and debonded sandwich plates from the previous example is carried out. A harmonic concentrated load is applied at the centre of the bottom (undamaged in the case of the debonded plate) face sheet and the driving frequency is swept within a given frequency interval from 800 Hz to 4000 Hz. This interval was chosen to encompass the biggest part of the frequency spectrum of the plates. Frequency response functions (FRFs) calculated at the central point of the top face sheet, which have been extracted from this analysis for both intact and debonded sandwich plates, are shown in Figure13. The difference between the FRF curves of those sandwich plates is very well visible. The shift in the resonance peaks and the existence of new ones on the FRF curve of the debonded plate are clearly demonstrated. Thus, this analysis visibly illustrates the influence of the debonding on the modal dynamics of sandwich panels; that is, both the frequency shifting and the presence of new frequencies in the frequency spectrum of the debonded plate are highlighted, as found in the previous frequency analysis.Figure 13
The FRFs obtained at the centre of the top face sheet of debonded and the same intact sandwich plates.Example 5.
A sandwich plate with material properties and geometry the same as in the previous example is assumed to be impacted by a hemispherical body. Thus, the sandwich plate with the impact-induced damage such thatRdent = Rcr = 30 mm as shown in Figure 7 is analysed. The influence of each key parameter of the impacted zone on the natural frequencies of the sandwich plate is examined. The calculations showed that the frequencies are practically not sensitive to the cavity depth, the face sheet indentation, and a level of the core crushing, and only the face sheet damage level causes small visible changes in the frequency spectrum, Figure 14. Herewith, the frequencies were changing in the manner mentioned above with increasing the plane size of the impacted zone, rd=Rdent only. That is, the latter parameter is a main factor that may manifest the impact damage within a sandwich plate if the linear model is used. It should be noted that mode shapes are more sensitive to changes of the key parameters of the impacted zone, but the differences between them are hardly recognized for each varying key parameter.Figure 14
The influence of the key parameters of the impacted zone on the natural frequencies [57].
(a)
(b)
(c)
(d)In this regard, further, the effect of the planar size of the debonded zone, its form, location within the plate area, and number of such zones (one or two) on natural frequencies and mode shapes is studied. The different boundary conditions and core properties are considered as well.Example 6.
A simply supported rectangular sandwich plate of 180 × 135 mm2 with a 7 mm thick honeycomb core and each 1 mm thick CFRP face sheet is considered. A penny-shaped debonded zone is located at the plate centre between the core and the top face sheet. The damage impaired by debonding is evaluated using the parameter D%. The material properties of the plate constituents are listed in Table 2. To study the influence of the core type, light H100 and heavy H200 PVC foams are used as a core material in calculations, Table 2. The main conclusion of those researches is that the bigger size of debonding, the larger its influence on the modal characteristics of sandwich plates regardless of the debonding form, the core type used, and the boundary conditions applied, while each of those mentioned parameters affects the modal dynamics in a different way. Moreover, the results obtained in the eigenvalue analysis illustrate the fact that with increasing the debonding size the frequency crossover and veering phenomena can be observed for certain vibrational modes, as shown in Figure 15.Figure 15
The influence of the debonding size on the natural frequencies [56].Example 7.
Next, numerical investigations are performed to study the effect debonding location. A simply supported rectangular sandwich plate of 270 mm length and 185 mm width with the face sheets the same as in Example2 and the core made of H100 PVC foam is examined. The thicknesses of the face sheets and the core are 2.4 mm and 50 mm, respectively. The plate contains either single or two equally sized penny-shaped debonded zones inflicting a total damage of D%=10%. Different positions of single debonding and two equally sized debondings within the plate area are considered, as shown in Figure 16.Figure 16
Location of penny-shaped debonded zones within the area of the sandwich plate:Di stands for a single debonding and DDi stands for two equally sized debondings, where i=1,4- is a position number.In Figure17(a) one can see that the first five natural frequencies of the debonded plates are shifted more from the intact one when debonding is situated near the corner (D4) and less for centrally located debonding (D1). The locations of debonding near the edges are intermediate cases between the mentioned two positions. Although debonding at the transverse edge more affects the frequencies, the higher modes follow the similar patterns of shifting with some exceptions for the several high modes.Figure 17
Natural frequency shifting for the debonded sandwich plate depending on location of the following: (a) single debonding and (b) two equally sized debondings.
(a)
(b)The comparative results for the debonded plate with two equally sized debondings are presented in Figure17(b). In general, the natural frequencies and their appropriate mode shapes highly depend on the location of the two debondings within the plate area. As seen from Figure 17(b) the two debondings, located along the longitudinal midline (DD2), most affect the frequencies of the plate. In the other cases of the debonding location, it influences on the frequency changes in different ways, but in a lesser extent than in DD2. Such tendency is retained for all modes from the lower to higher ones. More findings for this research can be found in detail in [56, 60], where this problem is better emphasized and extensively discussed with a large number of examples.


### 4.2. Explicit Dynamic Analysis

As mentioned earlier, the linear model is not able to account for contact and friction between the detached face sheet and core in the debonded zone, whereas nonlinear models can do it. The role of the contact behaviour on the dynamics of sandwich panels with debonding is evaluated in the next sections.Example 8.
An explicit transient dynamic analysis is carried out. In this regard, the sandwich plate with the impacted region such as in Example3 is excited with an impulse load applied to the centre of the health bottom face sheet. In the explicit analysis, the friction coefficient, μ, is accepted as 0.1, and the modal damping ratio is adopted as 1% of the critical value. The comparison of results obtained for the intact and the same debonded sandwich plates clearly demonstrated a high need of accounting for “real contact” modelling to properly capture the complex dynamic effects of debonded sandwich panels. The differences in transient time histories of the longitudinal strains calculated at the central point of the impacted top face sheet and the same point of the intact plate are shown in Figure 18(a). As one can see if contact is neglected, the calculations lead to sufficiently incorrect results, which mainly overestimate the amplitude of the dynamic response. Moreover, the importance of accounting for contact increases with enlarging the size of the impacted zone, as seen in Figure 18(b).Figure 18
The influence of contact on longitudinal strain [73]: (a) models with and without contact and (b) models with different sizes of the debonded zone.
(a)
(b)Besides, the computations revealed that varying the other key geometric parameters of the impacted zone (Figure7), visible changes can also be observed in the transient response of the damaged sandwich plate, Figure 19. That is, in contrast to the modal dynamics, the transient dynamic behaviour is sensitive to all the damage factors induced by impact. However, it is difficult to evaluate the severity of each parameter. In general, the time signals obtained for two different values of the parameter are not well-comparable between themselves. However, the debonding of bigger in-plane size can be easily differentiated from the transient signals for different sizes of the debonding.Figure 19
The transient behaviour of the impacted sandwich plate depending on [72] (a) indentation of face sheet; (b) cavity between face sheet and core; (c) level of core crushing; and (d) level of face sheet damage.
(a)
(b)
(c)
(d)In Figure20 one can see that the amplitude, time lag of free decay oscillations, and attenuation (viewed as trend-lines of oscillations) of displacement time curves calculated at the central plate point are clearly magnified with increasing the size of the debonded zone. This is an evidence of increasing the internal damping capacity in the debonded plate due to the presence of debonding. Thereby, the debonding introduces into the oscillations additional dissipative mechanisms caused by contact/impact and friction phenomena between the detached face sheet and core.Figure 20
Transient response for different sizes of debonding at the central plate point [75]: (a) displacement curves and (b) displacement curves trend-lines.
(a)
(b)Also, due to those interactions between the detached surfaces the contact-induced stresses arise in the debonded zone. These stresses may play a critical role in the localized effects and, as a result, in the global response of the sandwich plate, because they are able to contribute to development of the face sheet-to-core fracture mechanisms. Figure21 shows the transient time signals of the von Mises stresses calculated at the centre and the point of boundary of the debonded zone on the top face sheet. More detailed results, regarding the transient dynamic response and the dynamic stress state occurring in the debonded sandwich plate, can be found in our papers [10, 75, 79].Figure 21
Transient time history of the von Mises stress calculated at the points [79]: (a) centre of the top face sheet and (b) boundary of the debonded zone on the top face sheet.
(a)
(b)


### 4.3. Implicit Dynamic Analysis

Example 9.
To gain a deeper insight into nonlinear dynamics of sandwich plates with debonding, a long-term dynamic behaviour of a simply supported sandwich plate with a central penny-shaped debonded zone is studied. The geometrical and mechanical properties of the plate are adopted as such in Example6. The plate is subjected to harmonic loading with a driving frequency Ω at the central point of the health bottom face sheet. A wide range of the driving frequencies defined as a fraction of the fundamental frequency of the same intact sandwich plate, η=Ω/f0, is examined. One of the main trends, which could be identified from the simulations, is the fact that the dynamic response of the debonded sandwich plate is strongly frequency dependent, as shown in Figures 22–28. The time signals are presented at either the central point of the damaged top face sheet (N1), or the central points of the debonded region between the remaining plate (N2) and the detached skin (N3) or the points near boundary of the debonded region between the remaining plate (N4) and the detached skin (N5) in a time range, when the steady-state motion of the vibrating sandwich plate has been achieved.Figure 22
Steady-state dynamic response atη=1/2(Ω=500Hz) [80]: (a) displacement time histories at points N2 and N3; (b) phase portrait at point N1; and (c) frequency spectra of the displacement signals at point N1.
(a)
(b)
(c)Figure 23
Steady-state dynamic response atη=3/4(Ω=800Hz) [80]: (a) displacement time histories at points N2 and N3; (b) phase portrait at point N1; and (c) frequency spectra of the displacement signals at point N1.
(a)
(b)
(c)Figure 24
Steady-state dynamic response atη=2.0(Ω=2000Hz): (a) displacement time histories at points N2 and N3; (b) phase portrait at point N1; and (c) frequency spectra of the displacement signals at point N1.
(a)
(b)
(c)Figure 25
Steady-state dynamic response atη=3(Ω=3000Hz): (a) displacement time histories at points N2 and N3 [80]; (b) phase portrait at point N1; and (c) frequency spectra of the displacement signals at point N1.
(a)
(b)
(c)Figure 26
Steady-state dynamic response atη=4(Ω=4000Hz): ((a)-(b)) displacement and velocity time histories at points N2 and N3; (c) displacement time histories at points N4 and N5 [80]; (d) displacement time histories at point N1 [80]; (e) phase portrait at point N1 [80]; and (f) frequency spectra of the displacement signals at point N1.
(a)
(b)
(c)
(d)
(e)
(f)Figure 27
Steady-state dynamic response atη=5 (Ω=5000Hz): ((a)-(b)) displacement and velocity time histories at points N2 and N3; (c) frequency spectra of the displacement signals at point N1; (d) displacement time histories at point N1; (e) phase portrait at point N1 [80]; and (f) Poincaré section at point N1.
(a)
(b)
(c)
(d)
(e)
(f)Figure 28
Steady-state dynamic response atη=7 (Ω=7000Hz): ((a)-(b)) displacement and velocity time histories at points N2 and N3; (c) frequency spectra of the displacement signals at point N1; (d) displacement time histories at point N1; (e) phase portrait at point N1 [80]; and (f) Poincaré section at point N1 [80].
(a)
(b)
(c)
(d)
(e)
(f)Figure22 illustrates a periodic motion of the debonded sandwich plate with the frequency the same as the excitation one. This motion goes without contact between the detached face sheet and core. That is, there are no local effects in the debonded region and global oscillations of the plate occur only.Retaining the excitation amplitude and increasing the driving frequency lead to arising the contact behaviour between the detached segments in the sandwich plate. In Figure23 a periodic motion of the debonded sandwich plate with one contact per each excitation is demonstrated. It is clearly seen that the periodically colliding detached surfaces produce harmonics which are integer multiples of the forcing frequency.Although one contact per each excitation cycle also occurs in the case study shown in Figure24, the qualitative character of response of the debonded plate differs significantly from the previous case. The dynamic behaviour of the debonded plate is a period-one motion. And, in this case the oscillations of the debonded and intact plates are not coincident in the phase. This happens because impact-like contacts prevail in the interactions between the face sheet and the remaining part of the plate. Moreover, the orbit of phase trajectory of the debonded plate is remarkably deformed due to this contact.A period-two motion of the debonded sandwich plate has been detected at the driving frequencyη=3, Figure 25. The plate vibrates in such way that the detached skin and the remaining part come into contact twice per cycle of excitation. Such contacts generate both harmonics and halves of the harmonics of the excitation frequency in the spectral signals of the debonded plate.Figure26 presents a behaviour of the debonded plate that has no contact at the centre of the debonded zone, but there exist two contacts per each cycle of excitation at the points near the debonding boundary. It is worth noticing that the existence of such behaviour as partial contact (Figure 5(a)) has been established only due to 3D modelling of the dynamics of debonded plate. Because of such type of intermittent contact, the global dynamics of the debonded sandwich plate differs from the intact one. A phase shift between the displacement time signals of both plates is clearly seen in Figure 26(d).The general steady-state dynamic response of the debonded plate atη=5 is demonstrated in Figure 27. One can see a quasi-periodic motion of the plate at this driving frequency. The frequency spectrum along with the peaks at integer numbers of the excitation frequency and halves of the harmonics exhibits also sidebands consisting of at least two components clearly visible on the frequency-axis. These frequencies produced by almost regular contacts within the debonded zone are incommensurate with the excitation frequency of the external load.The harmonic excitation of the debonded sandwich plate at the driving frequency of about 7000 Hz, that is, (η≈7), is shown in Figure 28. One can see that contacts between the remaining plate part and the detached skin are not regular. The trajectory of the debonded plate in the phase plane does not repeat itself ever again for the excitation periods. Moreover, a set of points in the Poincaré section corresponds to irregular variations of the motion from one excitation cycle to another one. From the frequency spectrum graph, one can see that the spectral spikes are surrounded by a distribution of frequencies having the character of a broadband. Thereby, the oscillations of the debonded plate at the given driving frequency look like a chaotic motion.Although we went through each case and explained each of the dependencies found, one can highlight our papers [79, 80], which can be viewed for further insights and details in this area.


## 4.1. Linear Dynamic Analyses

Example 1.
To make evident the effect of “soft” core on the free vibration kinematics of the sandwich panel, the free vibration analysis of a fully bonded sandwich plate is considered first. A squared simply supported plate with a sidea=200 mm and the total thickness h=20 mm has been analysed. The sandwich plate layout corresponds to a five-layer structure 0/90/core/0/90, where face sheets are made of a composite anisotropic material, while the core is isotropic. The ratio between the face sheet and core thicknesses is 10, whereas the mismatching in material properties is the order of 103 as shown in Table 2.Table 2
Material properties of sandwich panels used in the calculations.
Examples
Ref.
Constituents
Example1
[10]
Face sheets,ρf = 1627 kgm−3
E
1, GPa
E
2 = E3, GPa
G
12 = G13, GPa
G
23, GPa
ν
12 = ν13
ν
23
131.0
10.34
6.895
6.205
0.22
0.49
Core,ρc = 97 kgm−3
E
1 = E2, GPa
E
3, GPa
G
12 = G13, GPa
G
23, GPa
ν
12
ν
13 = ν23
6.890  ×  10−3
6.890  ×  10−3
3.450  ×  10−3
3.450  ×  10−3
0.0
0.0
Example2
[53]
Face sheets
ρ
f = 4400 kgm−3, Ef = 36.0 GPa, ν = 0.3
Core,ρc = 52 kgm−3
E
c = 50.0 MPa, Gc = 21.0 MPa
Examples3–5
[73]
Face sheets,ρf=1650 kgm−3
E
1 = E3, GPa
E
2, GPa
G
12 = G23, GPa
G
13, GPa
ν
12
ν
13 = ν23
19.3
3.48
1.65
7.70
0.05
0.25
Core,ρc = 52 kgm−3, Ec = 85.0 MPa, Gc = 30.0 MPa, ν = 0.42
Examples6–9
[56]
Face sheets,ρf = 1650 kgm−3
E
1, GPa
E
2 = E3, GPa
G
12 = G13, GPa
G
23, GPa
ν
12
ν
13 = ν23
140
10
4.6
3.8
0.25
0.25
Honeycomb core,ρc = 57.17 kgm−3
E
1 = E2, GPa
E
3, GPa
G
12, GPa
G
13, GPa
G
23, GPa
ν
12=ν13=ν23
0.461 × 10−3
1.464
0.194 × 10−3
341.7 × 10−3
192.1 × 10−3
0.31
ρ
c
H
100 = 100 kgm−3, Ec = 105.0 MPa, Gc = 45.0 MPa, ν= 0.32
ρ
c
H
200= 200 kgm−3, Ec = 230.0 MPa, Gc = 85.0 MPa, ν = 0.33In the finite element mesh, one continuum shell element with the two-layer cross-section scheme across the thickness of the face sheets and five solid elements across the thickness of the core are used. In plane, the plate is discretized by 50 elements along each of its edges. Apart from the bending modes, the sandwich plate with flexible core exhibits in-plane and shear (twisting) and symmetric out-of-plane (pumping) modes, Figure8. For the latter modes, the core is removed from the corresponding contour plots for the better presentation of this phenomenon. The pumping modes indicate that the two face sheets move vertically in the opposite direction that relates to the changes of the plate height. Therefore, the layer-wise shell-solid finite element model proposed for modelling sandwich plates with a flexible core is able to capture such nonclassical effects.Figure 8
Contour plots of the mode shapes of sandwich plate with flexible core: (a) an in-plane mode; (b) a shear mode; (c) and (d) pumping modes [10].
(a)
(b)
(c)
(d)Example 2.
Firstly, the linear model with spring elements has been verified by comparing the natural frequencies and corresponding mode shapes of a debonded sandwich beam with those that are known in the literature [53]. A sandwich beam cored by foam with rectangular cross-section contains a deboned zone at its middle span, Figure 9. Material properties of the beam are listed in Table 2.Figure 9
A debonded sandwich beam with foam core [53].A comparison of first six natural frequencies of the debonded sandwich beam is presented in Table3. The mode shapes associated with the appropriate natural frequencies are shown in Figure 10.Table 3
Comparison of natural frequencies (Hz) of intact and debonded foam cored sandwich beams.
Mode
Intact
Debonded
[53]
[55]
Present
[53]
[55]
Present
(1)
289.3
293.52
293.46
288.98
293.52
293.07
(2)
683.3
722.96
707.09
388.32
360.81
433.67
(3)
1096.9
1139.51
1106.7
1093.2
1139.31
1093.2
(4)
1151.6
1545.91
1495.8
1146.9
1146.32
1132.0
(5)
1778.2
-
1818.7
1771.3
1769.9
1769.9
(6)
1895.3
1863.62
1918.1
1842.2
1948.35
2080.2Figure 10
Mode shapes of the debonded sandwich beam corresponding to natural frequencies: (a) 293.07 Hz, (b) 433.67 Hz, (c) 1093.2 Hz, (d) 1132.0 Hz, (e) 1769.9 Hz, and (f) 2080.2 Hz.
(a)
(b)
(c)
(d)
(e)
(f)It follows from the analysis that debonding reduces the natural frequencies compared to those of the same intact beam and this reduction is different for each mode. Apart from the global modal forms, the mode shapes of the debonded beam demonstrate both local and coupled shapes of vibration as well.Example 3.
The free vibration analysis of a sandwich plate containing a penny-shaped debonded zone is further demonstrated. The plate with a 50 mm thick foam core and each of the 2.4 mm thick face sheets is simply supported. The planar dimensions of the plate are of 270 mm by 180 mm. The radius of the debonded zone,rd, is of 20 mm; that is, the damage parameter defined as D%=πrd2/Atotal is about 5% of the plate area, Atotal. The material properties of each plate layer are presented in Table 2.The influence of debonding on the modal parameters of the sandwich plate reveals, first, the change of the order of frequencies and the existence of new ones in the spectrum. This happens due to local or coupled modes relating to debonding. Second, the reduction of the low natural frequencies occurs, especially, for bending modes. Herewith, the decreasing strongly depends on the mode shape. In the case of the natural frequencies associated with in-plane modes, there are practically no changes, and the changes are small for shear modes. Finally, the debonding affects the high frequencies in another manner than the lower ones, but there is no defined trend as a mode number increases. In the high part of the spectrum, there exist frequencies which are even higher than those of the intact plate; that is, the thickening phenomenon occurs [64]. The shift in natural frequencies as a result of comparing the healthy plate with the debonded one is presented in Figure 11, and several local mode shapes of the debonded sandwich plate are illustrated in Figure 12.Figure 11
The shift in natural frequencies of the debonded sandwich plate.Figure 12
Local mode shapes in the debonded sandwich plate: (a)(1,1); (b) (1,2); and (c) (1,3).
(a)
(b)
(c)Example 4.
In the next step, the mode-base steady-state dynamic analysis of the intact and debonded sandwich plates from the previous example is carried out. A harmonic concentrated load is applied at the centre of the bottom (undamaged in the case of the debonded plate) face sheet and the driving frequency is swept within a given frequency interval from 800 Hz to 4000 Hz. This interval was chosen to encompass the biggest part of the frequency spectrum of the plates. Frequency response functions (FRFs) calculated at the central point of the top face sheet, which have been extracted from this analysis for both intact and debonded sandwich plates, are shown in Figure13. The difference between the FRF curves of those sandwich plates is very well visible. The shift in the resonance peaks and the existence of new ones on the FRF curve of the debonded plate are clearly demonstrated. Thus, this analysis visibly illustrates the influence of the debonding on the modal dynamics of sandwich panels; that is, both the frequency shifting and the presence of new frequencies in the frequency spectrum of the debonded plate are highlighted, as found in the previous frequency analysis.Figure 13
The FRFs obtained at the centre of the top face sheet of debonded and the same intact sandwich plates.Example 5.
A sandwich plate with material properties and geometry the same as in the previous example is assumed to be impacted by a hemispherical body. Thus, the sandwich plate with the impact-induced damage such thatRdent = Rcr = 30 mm as shown in Figure 7 is analysed. The influence of each key parameter of the impacted zone on the natural frequencies of the sandwich plate is examined. The calculations showed that the frequencies are practically not sensitive to the cavity depth, the face sheet indentation, and a level of the core crushing, and only the face sheet damage level causes small visible changes in the frequency spectrum, Figure 14. Herewith, the frequencies were changing in the manner mentioned above with increasing the plane size of the impacted zone, rd=Rdent only. That is, the latter parameter is a main factor that may manifest the impact damage within a sandwich plate if the linear model is used. It should be noted that mode shapes are more sensitive to changes of the key parameters of the impacted zone, but the differences between them are hardly recognized for each varying key parameter.Figure 14
The influence of the key parameters of the impacted zone on the natural frequencies [57].
(a)
(b)
(c)
(d)In this regard, further, the effect of the planar size of the debonded zone, its form, location within the plate area, and number of such zones (one or two) on natural frequencies and mode shapes is studied. The different boundary conditions and core properties are considered as well.Example 6.
A simply supported rectangular sandwich plate of 180 × 135 mm2 with a 7 mm thick honeycomb core and each 1 mm thick CFRP face sheet is considered. A penny-shaped debonded zone is located at the plate centre between the core and the top face sheet. The damage impaired by debonding is evaluated using the parameter D%. The material properties of the plate constituents are listed in Table 2. To study the influence of the core type, light H100 and heavy H200 PVC foams are used as a core material in calculations, Table 2. The main conclusion of those researches is that the bigger size of debonding, the larger its influence on the modal characteristics of sandwich plates regardless of the debonding form, the core type used, and the boundary conditions applied, while each of those mentioned parameters affects the modal dynamics in a different way. Moreover, the results obtained in the eigenvalue analysis illustrate the fact that with increasing the debonding size the frequency crossover and veering phenomena can be observed for certain vibrational modes, as shown in Figure 15.Figure 15
The influence of the debonding size on the natural frequencies [56].Example 7.
Next, numerical investigations are performed to study the effect debonding location. A simply supported rectangular sandwich plate of 270 mm length and 185 mm width with the face sheets the same as in Example2 and the core made of H100 PVC foam is examined. The thicknesses of the face sheets and the core are 2.4 mm and 50 mm, respectively. The plate contains either single or two equally sized penny-shaped debonded zones inflicting a total damage of D%=10%. Different positions of single debonding and two equally sized debondings within the plate area are considered, as shown in Figure 16.Figure 16
Location of penny-shaped debonded zones within the area of the sandwich plate:Di stands for a single debonding and DDi stands for two equally sized debondings, where i=1,4- is a position number.In Figure17(a) one can see that the first five natural frequencies of the debonded plates are shifted more from the intact one when debonding is situated near the corner (D4) and less for centrally located debonding (D1). The locations of debonding near the edges are intermediate cases between the mentioned two positions. Although debonding at the transverse edge more affects the frequencies, the higher modes follow the similar patterns of shifting with some exceptions for the several high modes.Figure 17
Natural frequency shifting for the debonded sandwich plate depending on location of the following: (a) single debonding and (b) two equally sized debondings.
(a)
(b)The comparative results for the debonded plate with two equally sized debondings are presented in Figure17(b). In general, the natural frequencies and their appropriate mode shapes highly depend on the location of the two debondings within the plate area. As seen from Figure 17(b) the two debondings, located along the longitudinal midline (DD2), most affect the frequencies of the plate. In the other cases of the debonding location, it influences on the frequency changes in different ways, but in a lesser extent than in DD2. Such tendency is retained for all modes from the lower to higher ones. More findings for this research can be found in detail in [56, 60], where this problem is better emphasized and extensively discussed with a large number of examples.


## 4.2. Explicit Dynamic Analysis

As mentioned earlier, the linear model is not able to account for contact and friction between the detached face sheet and core in the debonded zone, whereas nonlinear models can do it. The role of the contact behaviour on the dynamics of sandwich panels with debonding is evaluated in the next sections.Example 8.
An explicit transient dynamic analysis is carried out. In this regard, the sandwich plate with the impacted region such as in Example3 is excited with an impulse load applied to the centre of the health bottom face sheet. In the explicit analysis, the friction coefficient, μ, is accepted as 0.1, and the modal damping ratio is adopted as 1% of the critical value. The comparison of results obtained for the intact and the same debonded sandwich plates clearly demonstrated a high need of accounting for “real contact” modelling to properly capture the complex dynamic effects of debonded sandwich panels. The differences in transient time histories of the longitudinal strains calculated at the central point of the impacted top face sheet and the same point of the intact plate are shown in Figure 18(a). As one can see if contact is neglected, the calculations lead to sufficiently incorrect results, which mainly overestimate the amplitude of the dynamic response. Moreover, the importance of accounting for contact increases with enlarging the size of the impacted zone, as seen in Figure 18(b).Figure 18
The influence of contact on longitudinal strain [73]: (a) models with and without contact and (b) models with different sizes of the debonded zone.
(a)
(b)Besides, the computations revealed that varying the other key geometric parameters of the impacted zone (Figure7), visible changes can also be observed in the transient response of the damaged sandwich plate, Figure 19. That is, in contrast to the modal dynamics, the transient dynamic behaviour is sensitive to all the damage factors induced by impact. However, it is difficult to evaluate the severity of each parameter. In general, the time signals obtained for two different values of the parameter are not well-comparable between themselves. However, the debonding of bigger in-plane size can be easily differentiated from the transient signals for different sizes of the debonding.Figure 19
The transient behaviour of the impacted sandwich plate depending on [72] (a) indentation of face sheet; (b) cavity between face sheet and core; (c) level of core crushing; and (d) level of face sheet damage.
(a)
(b)
(c)
(d)In Figure20 one can see that the amplitude, time lag of free decay oscillations, and attenuation (viewed as trend-lines of oscillations) of displacement time curves calculated at the central plate point are clearly magnified with increasing the size of the debonded zone. This is an evidence of increasing the internal damping capacity in the debonded plate due to the presence of debonding. Thereby, the debonding introduces into the oscillations additional dissipative mechanisms caused by contact/impact and friction phenomena between the detached face sheet and core.Figure 20
Transient response for different sizes of debonding at the central plate point [75]: (a) displacement curves and (b) displacement curves trend-lines.
(a)
(b)Also, due to those interactions between the detached surfaces the contact-induced stresses arise in the debonded zone. These stresses may play a critical role in the localized effects and, as a result, in the global response of the sandwich plate, because they are able to contribute to development of the face sheet-to-core fracture mechanisms. Figure21 shows the transient time signals of the von Mises stresses calculated at the centre and the point of boundary of the debonded zone on the top face sheet. More detailed results, regarding the transient dynamic response and the dynamic stress state occurring in the debonded sandwich plate, can be found in our papers [10, 75, 79].Figure 21
Transient time history of the von Mises stress calculated at the points [79]: (a) centre of the top face sheet and (b) boundary of the debonded zone on the top face sheet.
(a)
(b)


## 4.3. Implicit Dynamic Analysis

Example 9.
To gain a deeper insight into nonlinear dynamics of sandwich plates with debonding, a long-term dynamic behaviour of a simply supported sandwich plate with a central penny-shaped debonded zone is studied. The geometrical and mechanical properties of the plate are adopted as such in Example6. The plate is subjected to harmonic loading with a driving frequency Ω at the central point of the health bottom face sheet. A wide range of the driving frequencies defined as a fraction of the fundamental frequency of the same intact sandwich plate, η=Ω/f0, is examined. One of the main trends, which could be identified from the simulations, is the fact that the dynamic response of the debonded sandwich plate is strongly frequency dependent, as shown in Figures 22–28. The time signals are presented at either the central point of the damaged top face sheet (N1), or the central points of the debonded region between the remaining plate (N2) and the detached skin (N3) or the points near boundary of the debonded region between the remaining plate (N4) and the detached skin (N5) in a time range, when the steady-state motion of the vibrating sandwich plate has been achieved.Figure 22
Steady-state dynamic response atη=1/2(Ω=500Hz) [80]: (a) displacement time histories at points N2 and N3; (b) phase portrait at point N1; and (c) frequency spectra of the displacement signals at point N1.
(a)
(b)
(c)Figure 23
Steady-state dynamic response atη=3/4(Ω=800Hz) [80]: (a) displacement time histories at points N2 and N3; (b) phase portrait at point N1; and (c) frequency spectra of the displacement signals at point N1.
(a)
(b)
(c)Figure 24
Steady-state dynamic response atη=2.0(Ω=2000Hz): (a) displacement time histories at points N2 and N3; (b) phase portrait at point N1; and (c) frequency spectra of the displacement signals at point N1.
(a)
(b)
(c)Figure 25
Steady-state dynamic response atη=3(Ω=3000Hz): (a) displacement time histories at points N2 and N3 [80]; (b) phase portrait at point N1; and (c) frequency spectra of the displacement signals at point N1.
(a)
(b)
(c)Figure 26
Steady-state dynamic response atη=4(Ω=4000Hz): ((a)-(b)) displacement and velocity time histories at points N2 and N3; (c) displacement time histories at points N4 and N5 [80]; (d) displacement time histories at point N1 [80]; (e) phase portrait at point N1 [80]; and (f) frequency spectra of the displacement signals at point N1.
(a)
(b)
(c)
(d)
(e)
(f)Figure 27
Steady-state dynamic response atη=5 (Ω=5000Hz): ((a)-(b)) displacement and velocity time histories at points N2 and N3; (c) frequency spectra of the displacement signals at point N1; (d) displacement time histories at point N1; (e) phase portrait at point N1 [80]; and (f) Poincaré section at point N1.
(a)
(b)
(c)
(d)
(e)
(f)Figure 28
Steady-state dynamic response atη=7 (Ω=7000Hz): ((a)-(b)) displacement and velocity time histories at points N2 and N3; (c) frequency spectra of the displacement signals at point N1; (d) displacement time histories at point N1; (e) phase portrait at point N1 [80]; and (f) Poincaré section at point N1 [80].
(a)
(b)
(c)
(d)
(e)
(f)Figure22 illustrates a periodic motion of the debonded sandwich plate with the frequency the same as the excitation one. This motion goes without contact between the detached face sheet and core. That is, there are no local effects in the debonded region and global oscillations of the plate occur only.Retaining the excitation amplitude and increasing the driving frequency lead to arising the contact behaviour between the detached segments in the sandwich plate. In Figure23 a periodic motion of the debonded sandwich plate with one contact per each excitation is demonstrated. It is clearly seen that the periodically colliding detached surfaces produce harmonics which are integer multiples of the forcing frequency.Although one contact per each excitation cycle also occurs in the case study shown in Figure24, the qualitative character of response of the debonded plate differs significantly from the previous case. The dynamic behaviour of the debonded plate is a period-one motion. And, in this case the oscillations of the debonded and intact plates are not coincident in the phase. This happens because impact-like contacts prevail in the interactions between the face sheet and the remaining part of the plate. Moreover, the orbit of phase trajectory of the debonded plate is remarkably deformed due to this contact.A period-two motion of the debonded sandwich plate has been detected at the driving frequencyη=3, Figure 25. The plate vibrates in such way that the detached skin and the remaining part come into contact twice per cycle of excitation. Such contacts generate both harmonics and halves of the harmonics of the excitation frequency in the spectral signals of the debonded plate.Figure26 presents a behaviour of the debonded plate that has no contact at the centre of the debonded zone, but there exist two contacts per each cycle of excitation at the points near the debonding boundary. It is worth noticing that the existence of such behaviour as partial contact (Figure 5(a)) has been established only due to 3D modelling of the dynamics of debonded plate. Because of such type of intermittent contact, the global dynamics of the debonded sandwich plate differs from the intact one. A phase shift between the displacement time signals of both plates is clearly seen in Figure 26(d).The general steady-state dynamic response of the debonded plate atη=5 is demonstrated in Figure 27. One can see a quasi-periodic motion of the plate at this driving frequency. The frequency spectrum along with the peaks at integer numbers of the excitation frequency and halves of the harmonics exhibits also sidebands consisting of at least two components clearly visible on the frequency-axis. These frequencies produced by almost regular contacts within the debonded zone are incommensurate with the excitation frequency of the external load.The harmonic excitation of the debonded sandwich plate at the driving frequency of about 7000 Hz, that is, (η≈7), is shown in Figure 28. One can see that contacts between the remaining plate part and the detached skin are not regular. The trajectory of the debonded plate in the phase plane does not repeat itself ever again for the excitation periods. Moreover, a set of points in the Poincaré section corresponds to irregular variations of the motion from one excitation cycle to another one. From the frequency spectrum graph, one can see that the spectral spikes are surrounded by a distribution of frequencies having the character of a broadband. Thereby, the oscillations of the debonded plate at the given driving frequency look like a chaotic motion.Although we went through each case and explained each of the dependencies found, one can highlight our papers [79, 80], which can be viewed for further insights and details in this area.


## 5. Further Research

Although a number of studies on the dynamic behaviour of delaminated or debonded composite panels are widely reported in literature, research on the effects of delamination or debonding on strength and damage tolerance of composite panels subjected to dynamic loading is limited. The dynamic fracture analysis of debonded sandwich panels is challenging and demands expertise in damage and fracture mechanics and material science to guide this research. For better understanding and prediction of multiple failure mechanisms occurring in sandwich structures in a dynamic environment more reliable and efficient numerical methods should be devised. The most popular methods that can be potentially used and partially have been already been applied to for this purpose are briefly discussed with only few selected references below, because, as the authors may suppose, it is a matter of future more extensive research.


### 5.1. Virtual Crack Closure Technique (VCCT)

The delamination/debonding growth can be modelled by the VCCT. This method is based on the Griffith crack growth criterion within the framework of the linear elastic fracture mechanics (LEFM). That is, a debonding grows if the energy released at its propagation is equal to or larger than the energy required to create a new crack surface. This critical energy release rate or fracture toughness is derived from appropriate fracture tests. The current energy release is calculated under an assumption that the energy required to extend and to close the cracka over a length Δa are identical. A progressive delamination in laminated panels under static loading has been studied using the VCCT, for example, in [87–90], whereas studies on the delamination growth with VCCT under dynamic conditions can be found in few publications only, for example, in [91, 92].


### 5.2. Cohesive Zone Model (CZM)

A cohesive zone model allows for simulation of both the onset and growth of delamination/debonding in one approach. It uses cohesive elements placed at the interface between constitutive layers. The elements are endowed by a certain traction-separation law that combines a stress-based analysis to predict the onset of delamination/debonding and fracture mechanics based approach to govern its growth. In contrast to the VCCT, which uses one parameter, namely, the critical energy release rate, the CZM requires two parameters to be defined: fracture toughness and interfacial strength stresses. Modelling dynamic delamination growth in laminated panels with the CZM approach can be found, for example, in [93, 94], while the dynamic progression of debonding between the face sheet and the core in sandwich plates has been simulated with CZM, for example, in [95–97].


### 5.3. Extended Finite Element Method (XFEM)

The conventional finite element method is a versatile and powerful tool for the analysis of structural response of composite panels. Both the VCCT and CZM approaches can be implemented into the finite element formulation. However, the FEM uses a fixed mesh, within which the crack propagation is assumed as either releasing mesh nodes or removing finite elements. The XFEM no longer requires the mesh to be linked to discontinuities and their propagation. It uses enriched finite elements based on the partition of unity to make a complete or partial intersection of the finite element by the crack and introduces level sets to model the crack surface and the crack front, while maintaining the same mesh. Jump-like discontinuous functions are added to displacement interpolation functions for the completely intersected elements, while the elements where the crack front exists are supplemented by the asymptotic displacements fields resulting from the LEFM. The use of XFEM models for simulating quasi-static and dynamic delamination growth in laminated composites has been found, for example, in recent papers [98–100].


### 5.4. Peridynamics

Purpose of peridynamics is to unify the mechanics of continuous and discontinuous media within a single, consistent set of equations. This is reached with using integral equations, in contrast with the classical theory of continuum mechanics based on PDEs. In numerical implementation, the peridynamic approach is reduced to meshless and Lagrangian method, where the fracture analysis is carried out within the framework of the bond-based peridynamic theory. The latter incorporates damage in the pairwise force function allowing bonds to break when their elongation exceeds some prescribed value. After a bond breaks, it no longer sustains any force, and the endpoints are effectively disconnected from each other, and the force it was carrying is redistributed to other bonds that have not yet broken. This increased load makes it more likely that these other bonds will break too. The process of bond breakage and load redistribution, leading to further breakage, simulates cracks growth in the peridynamic model. Thus, the method enables predicting complex fracture phenomena such as spontaneous crack nucleation and crack branching, curving, and arrest. The model requires prescribed bond properties such as the stiffness constant and critical stretch, which are continuous functions of bond orientation in the principal material axes. An application of the peridynamics to predict delamination initiation and growth in laminated composites has been presented, for example, in [101, 102].


### 5.5. Phase Field Methods

A conceptually new computational approach to model fracture relying on phase field formulations has gained interest in recent years. Fracture phase field models make the use of a continuous scalar variable—the crack field order parameter to distinguish the fully broken and intact material phases. In such way, the sharp discontinuity is approximated through a smooth transition (diffuse transition zone) between the different values of the order parameter associated with the adjoining material phases. Thus, a crack phase field is set, while the evolution of this field due to loading models the fracture process. The description of fracture follows from the solution of a coupled nonlinear system of (quasi-static or dynamic) stress equilibrium equations and a gradient-type evolution equation for the crack phase field. Furthermore, the finite element implementation is straightforward in both two and three dimensions, using the classical FEM without any modification of initial meshes or shape functions. Examples of phase field models adopted for simulating fracture and damage in composite materials can be found, for example, in [103–105].The commercial finite element codes which may implement these techniques are ABAQUS, ANSYS, MSC Marc, MSC Nastran, and so forth. With this implementation, the performance and accuracy of those different approaches should be systematically investigated for both the numerical competence and the adequacy of the modelled behaviour to real structural response. Consequently, their application to the progressive delamination/debonding problems in practical structures is, as yet, limited and rather is a subject of future research activity.


## 5.1. Virtual Crack Closure Technique (VCCT)

The delamination/debonding growth can be modelled by the VCCT. This method is based on the Griffith crack growth criterion within the framework of the linear elastic fracture mechanics (LEFM). That is, a debonding grows if the energy released at its propagation is equal to or larger than the energy required to create a new crack surface. This critical energy release rate or fracture toughness is derived from appropriate fracture tests. The current energy release is calculated under an assumption that the energy required to extend and to close the cracka over a length Δa are identical. A progressive delamination in laminated panels under static loading has been studied using the VCCT, for example, in [87–90], whereas studies on the delamination growth with VCCT under dynamic conditions can be found in few publications only, for example, in [91, 92].


## 5.2. Cohesive Zone Model (CZM)

A cohesive zone model allows for simulation of both the onset and growth of delamination/debonding in one approach. It uses cohesive elements placed at the interface between constitutive layers. The elements are endowed by a certain traction-separation law that combines a stress-based analysis to predict the onset of delamination/debonding and fracture mechanics based approach to govern its growth. In contrast to the VCCT, which uses one parameter, namely, the critical energy release rate, the CZM requires two parameters to be defined: fracture toughness and interfacial strength stresses. Modelling dynamic delamination growth in laminated panels with the CZM approach can be found, for example, in [93, 94], while the dynamic progression of debonding between the face sheet and the core in sandwich plates has been simulated with CZM, for example, in [95–97].


## 5.3. Extended Finite Element Method (XFEM)

The conventional finite element method is a versatile and powerful tool for the analysis of structural response of composite panels. Both the VCCT and CZM approaches can be implemented into the finite element formulation. However, the FEM uses a fixed mesh, within which the crack propagation is assumed as either releasing mesh nodes or removing finite elements. The XFEM no longer requires the mesh to be linked to discontinuities and their propagation. It uses enriched finite elements based on the partition of unity to make a complete or partial intersection of the finite element by the crack and introduces level sets to model the crack surface and the crack front, while maintaining the same mesh. Jump-like discontinuous functions are added to displacement interpolation functions for the completely intersected elements, while the elements where the crack front exists are supplemented by the asymptotic displacements fields resulting from the LEFM. The use of XFEM models for simulating quasi-static and dynamic delamination growth in laminated composites has been found, for example, in recent papers [98–100].


## 5.4. Peridynamics

Purpose of peridynamics is to unify the mechanics of continuous and discontinuous media within a single, consistent set of equations. This is reached with using integral equations, in contrast with the classical theory of continuum mechanics based on PDEs. In numerical implementation, the peridynamic approach is reduced to meshless and Lagrangian method, where the fracture analysis is carried out within the framework of the bond-based peridynamic theory. The latter incorporates damage in the pairwise force function allowing bonds to break when their elongation exceeds some prescribed value. After a bond breaks, it no longer sustains any force, and the endpoints are effectively disconnected from each other, and the force it was carrying is redistributed to other bonds that have not yet broken. This increased load makes it more likely that these other bonds will break too. The process of bond breakage and load redistribution, leading to further breakage, simulates cracks growth in the peridynamic model. Thus, the method enables predicting complex fracture phenomena such as spontaneous crack nucleation and crack branching, curving, and arrest. The model requires prescribed bond properties such as the stiffness constant and critical stretch, which are continuous functions of bond orientation in the principal material axes. An application of the peridynamics to predict delamination initiation and growth in laminated composites has been presented, for example, in [101, 102].


## 5.5. Phase Field Methods

A conceptually new computational approach to model fracture relying on phase field formulations has gained interest in recent years. Fracture phase field models make the use of a continuous scalar variable—the crack field order parameter to distinguish the fully broken and intact material phases. In such way, the sharp discontinuity is approximated through a smooth transition (diffuse transition zone) between the different values of the order parameter associated with the adjoining material phases. Thus, a crack phase field is set, while the evolution of this field due to loading models the fracture process. The description of fracture follows from the solution of a coupled nonlinear system of (quasi-static or dynamic) stress equilibrium equations and a gradient-type evolution equation for the crack phase field. Furthermore, the finite element implementation is straightforward in both two and three dimensions, using the classical FEM without any modification of initial meshes or shape functions. Examples of phase field models adopted for simulating fracture and damage in composite materials can be found, for example, in [103–105].The commercial finite element codes which may implement these techniques are ABAQUS, ANSYS, MSC Marc, MSC Nastran, and so forth. With this implementation, the performance and accuracy of those different approaches should be systematically investigated for both the numerical competence and the adequacy of the modelled behaviour to real structural response. Consequently, their application to the progressive delamination/debonding problems in practical structures is, as yet, limited and rather is a subject of future research activity.


## 6. Conclusions

Debonding is an inherent potential reason for the structural failure in sandwich panels possessing the adhesively bonded structure. Hence, detailed knowledge about its influence on the dynamic behaviour of the sandwich composite is highly required. In this regard, reliable numerical prediction methods are of importance. In this paper the linear and nonlinear finite element models developed by the authors to examine the dynamics of sandwich panels with preexisting debonding are comprehensively reviewed and prospects of further investigations on this topic are briefly discussed. The structural problem at hand has been numerically solved using the ABAQUS package. Numerical results ranging from linear free and forced vibration of debonded sandwich panels to their general nonlinear dynamics accounting for contact and friction are presented. A typical sandwich plate configuration has been discretized with shell and solid finite elements in accordance with the HSAPT-like assumptions. In the modal dynamic analyses a spring finite element separated the face sheet from the core in the debonded region, while contact and friction laws representing the unique sources of nonlinearity are imposed between them in the case of general dynamic analysis.Observations from the present work suggest the following conclusions. First, one can point out that the pumping effect or symmetric vibration modes occurring in sandwich panels with soft flexible cores can be efficiently simulated using the proposed layer-wise shell-solid finite element models. This result is of importance for developing 3D models that are able to accurately predict the dynamic behaviour of debonded sandwich panels. Second, the 3D finite element models created using the mentioned modelling strategy in conjunction with spring elements are relevant for solving a range of problems of the modal dynamics of sandwich panels with both an artificial flaw embedded into the face sheet/core bond and a postimpact circular damage involving indentation, face sheet/core interface degradation and defect of core. The calculated natural frequencies and associated mode shapes have been in good compliance with results obtained by some other authors using other methods. The various parametric studies carried out in terms of debonding parameters such as size, location, and number may deserve an attention for applying these results to developing appropriate structural health monitoring techniques. The major findings from the present studies are general decreasing of natural frequencies of debonded plates compared to the same intact ones, possible thickening phenomenon, changing in the order of frequencies in the spectrum and/or appearance of new frequencies there, and potential existence of the crossover phenomenon. The steady-state mode-based analyses based on the same models highlighted the occurrence of new frequencies on the frequency response curve. Meanwhile, it was observed that the mode shapes may give less useful information in the sense of debonding detection since they are highly dependent on location of the debonded zone within the plate area.In the present paper it is also noticed that although many advances have been made in simulation and analysis of the modal dynamics of sandwich panels with damaged face sheet-to-core interface, the general dynamic analysis of such structures remains challenging and requires expertise in physical nature of the behaviour and accurate numerical methods applied to guide further theoretical research. Local debonding may lead to significant changes in the global dynamics such that the dynamic response may become nonlinear even under initially linear conditions. This is the next conclusion that follows from this research. Nonlinearities are mainly related to intermittent contact and friction between the detached face sheet and core during global oscillations of the plate. The present numerical studies clearly showed that, under long-term harmonic loading, the changes of the driving frequency leads to significant alterations in the type of interactions between the debonded face sheet and the remaining part and as a result the dynamic response of the sandwich plate may be either periodic, or only quasi-periodic, or even chaotic. In the case of an impulse load the major finding is that transient time history responses are very sensitive to the size and severity of the debonding. Thus, the nonlinear models’ methodology presented herein can be deemed to be a promising prediction approach for discovering a variety of physical phenomena occurring in debonded sandwich panels and for examining the effect of each of them on the dynamics at the design stage. The obtained results may give a possibility of devising new structural health monitoring techniques which are more sensitive than the ones using modal dynamics to identify the delamination/debonding damage on the structural level. Besides, the present work mentions that there is a need to extend the study of debonding on strength and damage tolerance of sandwich plates subjected to dynamic loads, because the debonding onset and the debonding evolution are critical issues to ensure a structural reliability of sandwich panels.



---

*Source: 5715863-2018-02-14.xml*