Vibration of a laminated beam with a delamination including contact effects

Certain results are presented in this paper on damped vibration of a laminated cantilever beam with a single closing delamination. In order to investigate this task the finite element method has been applied in the current study. For modelling the beam higher order shear deformation beam finite elements have been used. The vibration of the beam is investigated in the time domain using a dynamic contact algorithm developed by the authors. The algorithm is based on the Newmark method and also incorporates a Newton-Raphson based procedure for resolving the equation of motion. The time series obtained from solving the equation of motion have been subsequently analysed in the frequency domain by using FFT (Fast Fourier Transform). The vibration responses of the beam due to various harmonic and impulse excitations, at different delamination locations, and for different delamination lengths, as well as changes in the dissipation of damping energy due to the delamination, have all been considered in the paper.


Introduction
In current engineering practice composite materials are widely and commonly used for numerous and various applications in areas such like civil engineering, mechanical engineering, aircraft industry, and space technologies, everywhere where high strength-to-weight ratios are required.However, structural elements made of composite materials, similar to those made of isotropic materials, are exposed to many different types of damage which may occur during their routine work and service.In the case of composite materials the most important types of damage are cracking of matrix or reinforcing fibres, or delamination of layers.
For structural elements made of laminated composite materials (laminates) delamination is one of the most important and most dangerous failure modes.Delamination introduced during manufacturing processes, being an effect of impact or another operational risk, may significantly influence and reduce the structural stiffness and the buckling load of a structure.As a consequence of that its static and dynamic characteristics can be substantially altered leading even to the total destruction of the structure.For these reasons extensive work has been carried out in past decades in this field in many scientific and research institutes and has been very well documented in the literature.
The influence of delamination on the buckling and post-buckling deformation, and the influence of different geometrical parameters like loading conditions, material properties, and boundary conditions, on the delamination growth were studied in the past by Yin et al. [22], Chen [2], Yin and Jane [23], and Jane and Chen [6].Shu [18] presented analytical and experimental results on the influence of a double delamination on vibration of sandwich beams.Similarly, the effect of multiple delaminations on the dynamics of axially compressed beams was studied by Lee et al. [12].Various aspects of the delamination influence on the dynamic characteristics of composite beams and plates were analysed numerically and experimentally by Mujumdar and Suryanarayan [13], Han et al. [4], Lai and Young [11], Krawczuk et al. [9], Gadelrab [3], Saravanos and Hopkins [17], Żak et al. [24], and Ostachowicz and Kaczmarczyk [15].
Another very important aspect of the work carried out in this field is damage detection based on known, simulated or experimentally measured, changes in dynamic characteristics of structural elements made of composite materials.Using known and analytically obtained dynamic characteristics Valoor and Chandrashekhara [20] used the method of neural networks for detection of the delamination in composite beams.Tadoroki et al. [19] presented results of an experimental investigation on the delamination detection in composite plates based on changes in their resistance.Palacz and Krawczuk [16] studied in their work the usefulness of different vibration-based damage detection criteria, however, the analytically obtained dynamic characteristics they used were for an isotropic cracked cantilever beam.A thorough review was presented by Zou et al. [26] on different delamination indicators useful in health monitoring of composite elements of structures based on various vibration characteristics.
However, very few papers related to the subject have considered other dynamic characteristics than those based on the natural frequencies or modes of vibration, and which could be useful in damage detection.In their review Zou et al. [26] proposed damping as potentially a very sensitive and attractive parameter as a delamination indicator.Zhang and Hartwig [25] studied the relation between damping and damage in composite materials, but their study was focused on matrix and fibre cracking rather than delamination of layers.Similarly, Kisa and Brandon [8] investigated the influence of damage closure on the dynamic behaviour of a cracked composite beam.
It appears clearly from the above review that the influence of the delamination on various dynamic characteristics, as well as its detection in structural elements made of composite materials has been well documented in the literature and a vast number of papers has been published on this subject.However, in the majority of them their authors neglected the effect of contact between delaminated layers primarily due to its complex nature, as well as the necessity of solving the non-linear equation of motion, which highly complicates the studied dynamic problem.In this work, however, the effect of contact in the area of delamination has been taken into account and this can significantly improve understanding of the dynamic behaviour of damaged composite structures.The results presented are intended to fill in that gap in the present literature related to this subject offering a better insight into the complex dynamical behaviour of such structures.
In order to investigate damped vibration of a laminated cantilever beam with a single closing delamination the finite element method has been applied.Higher order shear deformation theory finite elements have been used by the authors for modelling the beam and a new dynamic contact algorithm developed by the authors has been utilised for solving the equation of motion in the time domain.The contact algorithm is based on the Newmark method and also incorporates a Newton-Raphson based procedure for solving the equation of motion.The time series obtained from solving the equation of motion are subsequently analysed in the frequency domain by using FFT (Fast Fourier Transform).The vibration responses of the beam due to various harmonic and impulse excitations, at different delamination locations, and for different delamination sizes, as well as changes in the dissipation of damping energy due to the delamination, have all been considered in the paper.

Laminated cantilever beam with a delamination
Figure 1 shows a laminated cantilever beam with a single delamination.The beam has the following dimensions: the length L is 1000 mm, the width B is 20 mm, and the thickness H is 10 mm.All the material layers within the beam are located symmetrically with respect to the x-z plane.It has been assumed that the beam is made out of glass/epoxy laminate and consists of 10 viscoelastic material layers.The ply stacking sequence is [±45 • ] 5 and the relative volume fraction of the reinforcing glass fibres is 0.5.Selected mechanical properties of the beam material components are shown in Table 1.It should be emphasised here that in this work all the material properties have been assumed as frequency independent.The single delamination of a length a is arbitrarily located within the beam and extends to its full width.The centre of the delamination it located at the distance L 1 from the free end of the beam, and at the distance H 1 from the upper surface of the beam.The beam is modelled by the use of beam finite elements and the higher order shear deformation theory is employed, as presented by Ochoa and Reddy [14].Each element has 3 nodes and 3 degrees of freedom at each node (i.e.longitudinal and transverse displacements and an independent rotation).
According to the higher order shear deformation theory, and for the theory of small displacements, the displacement field within a single lamina of each beam finite element can be expressed in the following, well-known, form: where u(x, z) and w(x, z) are the longitudinal and transverse displacements of a point (x, z) within the lamina, u 0 (x) and w 0 (x) are the displacements of a point in the mid-plane of the beam, and θ(x) is the independent rotation of the transverse normal about the z axis.
The strain field within a single lamina is expressed according to the theory of small deformations by using the following relations: The inertia M and stiffness K matrices of the whole element can be easily evaluated by applying standard finite element formulae -see Appendix for details.These matrices are expressed as sums of the inertia M and stiffness K matrices calculated for each individual lamina, as presented by Żak et al. [24].
Along the delamination interface the beam is split into two distinct layers of finite elements.These layers, modelling the delamination, have been connected together at the delamination start and end at the nodes of appropriate finite elements.The present method of the delamination modelling corresponds to the rigid connector presented by Shu [18] (i.e. the cross-section of the beam perpendicular to the mid-plane remains perpendicular to the deformed mid-plane).In this way no other additional boundary conditions are required at the interface.
A standard convergence analysis has been performed to verify the accuracy of the chosen modelling technique.For this case, however, the simply-supported type of boundary condition of the beam has been selected, rather than 2 First-order shear deformation theory [21].
the cantilever type, due to the availability of appropriate analytical formulas.It has been also assumed that the delamination is not present.Results obtained for the natural frequencies of the beam, as well as the analytical results obtained from the literature [21], are presented in Table 2.It can be seen that a sufficient accuracy of the results can be achieved (the highest relative error around 10% for the mode VI) when the beam is divided only into 10 beam finite elements.However, in order to correctly model the dynamic phenomena in the delamination area, where the contact occurs, as well as to consider short delaminations below 10% of the total beam length, it has been decided that in the following results presented in this work the beam will be divided into 40 beam finite elements, or more, depending on the case.

Contact modelling and solution procedure
For the purpose of the present investigation a discrete contact model between contacting element nodes has been utilised.Additionally, it has been assumed that there is no friction in the delaminated area and also due to significant computational simplifications the lumped form of the element mass matrix M, rather than the consistent one, has been chosen.
In order to solve the equation of motion the implicit Newmark algorithm has been applied.During the computations the equation of motion is solved in the following form: where M, C and K are the mass, damping and stiffness matrices, respectively.The symbol Q denotes the vector of external loads depending on time t, and F is the vector of internal forces.The symbol ∆t is a time increment, while ∆r, ν, and a denote vectors of the generalised nodal displacement increment, nodal velocities and accelerations.Equation ( 3) must be considered with additional compatibility conditions for the transverse components of the nodal displacements y i , y j and velocities u i , u j , for an unknown pair of contact nodes i and j in contact, in an unknown contact time instance t c .These conditions express the absence of penetration between the nodes in contact as well as the conservation of momentum.
When on the time step ∆t penetrations p(t) for pairs of contacting nodes are detected the algorithm adjusts the current time step calculating a new time step ∆t c , for which these penetrations are zero (i.e.y i (t c ) − y j (t c ) = (0), and only one pair of contacting nodes remains in contact.These penetrations are defined as follows: where the symbols y i and y j are the transverse components of nodal displacements of contacting nodes i and j.The symbol δ is a small negative constant assumed as equal to −1 × 10 −6 .Because the coefficients of the Newmark method are functions of the time step ∆t, the new time step ∆t c is calculated iteratively using a Newton-Raphson based procedure.In this way the compatibility between the transverse components of the nodal displacements of the contacting nodes is achieved.Additional compatibility conditions must be imposed on the nodal velocities by utilising the conservation of momentum in the transverse direction: where m i and m j are the masses of contacting nodes i and j, respectively -similarly to Kwon and Aygunes [10].
The symbols u i , u j are the transverse components of nodal velocities calculated at the end of the new time step ∆t c .The symbols u + i , u + j and u + are the values coming up after applying the compatibility conditions, and k is a coefficient of restitution.Once the compatibility conditions (4) have been applied the algorithm returns to the original time step ∆t and the calculations continue.The flow diagram of the discussed algorithm is also presented in Fig. 2.
It should be noted here that the implicit Newmark method is, in general, characterised by unconditional stability, however in the case of the vibration problem studied in this work, certain numerical difficulties occurred when too a large time step was selected.This was particularly the case when a finite element model contained many contacting nodes.In this situation the selected time step had to be small enough to allow the algorithm to detect each particular instance of contact -numerical simulations suggested that the time step should be one or two orders higher than that resulting from, for example, the excitation frequency.
The effectiveness of the presented algorithm is shown in Fig. 3.In this case the beam has been excited at its free end by a transverse force impulse of 100 N, acting over a period of time 0.5 ms.A short delamination (a/L = 0.1) has been located in the middle of the beam length (L 1 /L = 0.5), and in the mid-plane of the beam (H 1 /H = 0.5).The beam is modelled by 44 finite elements, which includes 8 finite elements in the delaminated area (4 elements for the upper layer and 4 elements for the lower layer).The coefficient of restitution k has been assumed to be unity in this test, and the same value of the coefficient has been used for the following results presented in the paper.
The transverse displacements of the beam, shown in Fig. 3, have been calculated for the upper and lower nodes at the delamination centre.

Numerical results
In this section of the present work different dynamic responses of the considered beam with a delamination to impulse or harmonic excitations have been presented.These responses have been obtained for a variety of locations and lengths of the delamination, as well as for different excitation locations and excitation frequencies -in the case of the harmonic excitation.The time series obtained from solving the equation of motion have been analysed in the frequency domain by using FFT (Fast Fourier Transform).
As a first, however, two results from linear analysis are presented showing the influence of the delamination location and length on the first five bending modes of the beam -see also Zou et al. [26].
It can be seen from the results presented in Fig. 4 that the influence of the delamination location on the natural frequencies of the beam increases with the mode number, but the present increase has oscillatory characteristics strongly depending on the delamination location.A similar conclusion can be drawn based on the results shown in Fig. 5.The influence of the delamination length is greater on the higher modes of vibration and rapidly increases with the delamination length.However it should be noticed that for relatively short delaminations (a/L = 0.1) the natural frequencies are reduced by only 5% or even less.This can cause great difficulties in detection of short delaminations especially when based on changes in the natural frequencies of the beam.On the other hand using other dynamic characteristics than the natural frequencies or modes of vibration (i.e.amplitudes of forced vibration, damping) may help to solve that problem.The following results presented in this paper have been obtained from a non-linear analysis, when the contact in the delamination has been taken into account.It has been also decided that two types of beam excitations are considered here: a harmonic excitation by a transverse force of 100 N at 1 kHz or 2 kHz frequency, and a smooth impulse excitation by a transverse force of 100 N acting over 0.5 ms period of time.It has been assumed that two different delamination lengths are considered, i.e. 5% and 10% of the total beam length (a/L = 0.05 or 0.1).In each presented case the delamination is located at the centre of the beam (L 1 /L = 0.5), and positioned either in the mid-plane of the beam or one layer towards the upper surface (H 1 /H = 0.5 or 0.4).Vibration responses of the beam have been calculated in the form of amplitude-frequency or force-displacement diagrams, and are shown in graphs for three different locations along the beam length: 50%, 75% and 100% of the total beam length (L 1 /L = 0.5, 0.75, and 1.0).
The results of numerical calculations indicate a strong relation between the excitation and the delamination length and location.However it has been found for the beam excited at its free end, and for the delamination up to 10% of the total beam length, that practically no changes are observed in the obtained amplitude-frequency diagrams in comparison with the amplitude-frequency diagrams obtained for the beam without the delamination.No additional peaks resulting from delamination closing and opening can be seen as long as the delamination is located in the mid-plane of the beam.This is independent of the excitation frequency and the same results are obtained for the 1kHz and 2 kHz harmonic excitations.In contrast, when the delamination moves outwards from the mid-plane of the beam evident changes in the amplitude-frequency diagrams may be seen.The magnitude of these changes depends  on both the delamination length and the excitation frequency.The behaviour described above is well illustrated by Figs 6-9.It can be clearly seen that due to the non-linear effects visible additional peaks can be seen, which are multiples of the excitation frequencies, and which result from periodical closing and opening of the delamination.In the case of the longer 10% delamination these additional peaks also propagate from the delamination area and are observed along the whole considered length of the beam.
Obviously the beam response significantly changes when the beam is excited in the delamination area.It has been established that the presence of the additional peaks, resulting from the closing and opening of the delamination is in this case independent of the delamination length and location, as well as independent of the excitation frequency.However, it should be mentioned here that the intensity of the noticeable non-linear response is a function of these parameters, as shown in Figs 10-11.
Also very interesting results have been obtained in the case when the beam is excited to vibration by the impulse excitation.In all considered cases here of the delamination locations and lengths, as well as the excitation location, no visible changes resulting from the non-linear effects due to closing and opening of the delamination can be seen in the amplitude-frequency diagrams.Certain results obtained for the investigated beam are presented in Figs 12-13.As clearly seen from Figs 12-13 the vibration response of the beam weakens strongly while propagating from the excitation area towards the free end of the beam and contains only resonant peaks.No significant differences can be seen between the response of the beam without the delamination and when the delamination length is considerable and reaches even 10% of the total beam length.The authors believes that the absence of the non-linear response (additional peaks) in the amplitude-frequency diagrams results from strong viscous material damping in this case.
The presence of the delamination changes not only the amplitude-frequency characteristics of the investigated beam.Because of the fact that the beam material has viscoelastic properties the delamination also strongly affects, and significant alters, hysteresis loops of the beam due to the harmonic excitation.These changes can be seen in both the sizes of the hysteresis loops and in their tangential slopes.Moreover, detectible changes can be seen for all considered lengths and positions of the delamination, as well as the excitation frequency.This is clearly seen from the results presented in Fig 14 -18.
The observed changes in the sizes and slopes of the hysteresis loops are the effect of the changes in the stiffness of the beam resulting from the presence of the delamination.Furthermore, the non-linear effects arising from the delamination closing and opening additionally alter the stiffness of the beam.Due to the periodical nature of this alteration the shapes of the hysteresis loops may be additionally distorted when compared with the shapes the loops obtained in the absence of the delamination -see Figs 17-18 for details.

Conclusions
In this paper certain results have been presented on damped vibration of a laminated cantilever beam with a single closing delamination.The vibration responses of the beam studied in this work allow for the following, general, conclusions to be drawn: 1.The results obtained from a linear analysis show that the presence of the delamination reduces the natural frequencies of the beam, as well as alters its modes of vibration.These changes strongly depend on the delamination location and length, but in general the observed changes are small for short delaminations.For that reason the usage of these parameters as damage indicators is very limited.2. In the case of the non-linear analysis additional resonant peaks, resulting from periodical closing and opening of the delamination, may be seen in the vibration responses of the beam.The presence of these additional resonant peaks, as well as their magnitude and intensity, are strongly linked to the delamination location between material layers and the delamination length, as well as the type of excitation.It has been found that the presence of the additional resonant peaks only weakly depends on the delamination location along the beam length.3. When the beam is excited to vibration by a harmonic force applied at the delamination area the additional resonant peaks may be clearly seen regardless of the delamination location and length.When the excitation is applied at the free end of the beam the additional resonant peaks are present in the beam's responses only if the delamination is located outside the mid-plane of the beam.In certain cases the additional resonant peaks may be detected along the whole beam length -not only in the delamination area.This is especially well seen for higher frequency excitations and longer delaminations.4. The additional resonant peaks do not appear in the vibration responses of the beam when the beam is excited to vibration by a force impulse.The responses of the beam in each considered case contain only "linear" resonant peaks.This behaviour, in the author's opinion, is directly related to strong viscoelastic damping of the beam's material.5.When hysteresis loops of the transverse displacements for the beam are investigated noticeable changes in their sizes and slopes can be observed due to the delamination presence and material damping.These changes are very sensitive to both the delamination location and length, as well as the excitation frequency.Moreover, due to the periodical nature of the contact in the delamination area, for longer delaminations and higher excitation frequencies, the shape of these loops may be additionally distorted.

Appendix
The laminated beam investigated in the paper has been modelled by beam finite elements and the higher order shear deformation theory [14] have been applied.Each beam finite element has 3 nodes and 3 degrees of freedom at each node (i.e.longitudinal and transverse displacements and an independent rotation).
The displacement field within a single lamina of each beam finite element is expressed by Eq. ( 1).The mid-plane longitudinal and transverse displacements u 0 (x) and w 0 (x), and the independent rotation θ(x), can be approximated by second order polynomials:    u 0 (x) = a 1 + a 2 x + a 3 x 2 w 0 (x) = a 4 + a 5 x + a 6 x 2 θ(x) = a 7 + a 8 (x) + a 9 x 2 (A.1)By using the conditions at the nodes of the element the unknown coefficients of the approximating polynomials a 1 -a 9 can evaluated and expressed in terms of the element generalised nodal displacements q 1 -q 9 as: where the matrix A has the following form: Next from the relationships between the coefficients of the approximating polynomials a 1 -a 9 and the mid-plane displacements u 0 (x) and w 0 (x), and the rotation θ(x), the matrix of element shape functions N can be easily found: where the matrix X is defined as: In a similar manner the strain filed within a single lamina of each beam finite element can be expressed in terms of the element generalised nodal displacements q 1 -q 9 by taking into account Eq. (2) as well as the relations (A.1) and (A.2): where the matrix Y is defined as: Finally, the inertia matrix M and the stiffness K of the whole finite element can be calculated from the following relations by applying the standard finite element procedure:

Fig. 2 .
Fig. 2. The flow diagram of the algorithm.

Fig. 3 .Fig. 4 .
Fig.3.A time response of the transverse displacements of a laminated cantilever beam with a delamination to an impulse excitation at the delamination centre (a/L = 0.1, L 1 /L = 0.5, H 1 /H = 0.4).

Fig. 5 .
Fig.5.The influence of the relative delamination length a/L on relative natural frequencies of a laminated cantilever beam (L 1 /L = 0.5, H 1 /H = 0.5).

Fig. 6 .
Fig.6.An amplitude-frequency diagram of a laminated cantilever beam without a delamination, excited by 1 kHz harmonic force at its free end.

Fig. 17 .Fig. 18 .
Fig.17.Hysteresis loops of the transverse displacements for a laminated cantilever beam without a delamination, excited by 2 kHz harmonic force at half length of the beam.

Table 2
The convergence of the finite element results for a simply-supported laminated beam without a delamination.