# Dynamic Stability of a Structurally Damped Delaminated Beam Using Higher Order Theory

Authors: Tamás Pölöskei; András Szekrényes
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

The static and dynamic stability of the composite beam with a single delamination are investigated using the Timoshenko beam theory. The mechanical model is discretized using the finite element method and the equation of motion is obtained using Hamilton’s principle. The coefficients of the mass and stiffness matrix for the damping matrix are determined using experimental modal analysis. The effect of harmonic excitation on the dynamic stability of a single delaminated composite beam is investigated using Bolotin’s harmonic balance method. The stability boundaries of the damped and undamped system are compared for different static load values and delamination lengths on the excitation frequency-excitation force amplitude parameter field.


---

## Body

## 1. Introduction

The first application of composite material is approximately more than 6000 years old and during this period of time the emphasis from the utilisation as a building material moved greatly towards the engineering applications with the requirement for high performance material properties. In the 21st century the application is ranging between construction, transportation, electronics, sports, medical, and aerospace industry [1]. This broad application palette makes it important to create accurate mechanical models of these mechanical structures.Natural frequency and static buckling analyses of intact straight beams and plates [2] as well as circular arches [3–5] are a widely investigated topic. However laminated composite structures due to their heterogeneous material are able to become partially or totally damaged in the presence of fabrication defects, impact loading, or chemical corrosion [6]. In case of partial damage on the macroscopic scale the most prominent defect is the delamination.One of the first modelling approaches to investigate the delamination was introduced by Wang et al. [7] in which the membrane coupling between the intact part and delaminated part as well as the rotary inertia was neglected during the natural frequency calculations. Later Mujumdar and Suryanarayan [8] introduced the constrained mode model which resulted in better accordance with the measurements. Finally Shen and Grady [9] implemented the missing compatibility between the membrane displacements and carried out further investigations into the free vibration analysis. Shu [10] and Shu and Fan [11] investigated the effect of the delamination front displacements and the number of delaminations on the primary frequencies. Luo and Hanagud [12, 13] considered the mechanics of the delamination surfaces not as totally constrained or free but the mix of two and applied nonlinear springs to achieve this effect. Della and Shu in [14, 15] made further observations on beams with multiple delaminations and in [16] on the effect of material properties of bimaterial split beams. Up to this point the introduced models were using equivalent single layer theories; however another large group of modelling techniques is the layer-wise approach [17]. In the subsequent sections the equivalent single layer theory will be used entirely.Another approach for determining the natural frequencies in delaminated structures is the use of finite element method. Most of the articles were made using classical or first-order shear deformable theories. Timoshenko’s beam theory [18] which will be later used throughout this paper belongs to the latter group. It has several types of discretization formulation [19–23]; however we selected one which is based on the equilibrium equations. Krawczuk et al. [24] developed the finite element model of the cracked beam and plate using first-order shear deformable theories and the obtained results for the natural frequencies were in good agreement with the measured ones.A more complex analysis is required if one wishes to determine the dynamic stability boundaries of elastic structures under time periodic loads. The equation which derives the phenomenon is a Mathieu equation. It is known from literature that in case of one-degree-of-freedom systems in the absence of damping the Lindstedt-Poincare method [25] is applicable. However in bigger systems the use of Bolotin’s harmonic balance method [26] is easier and there is an algorithmic approach to increase the order of approximation. The method was successfully applied in many types of systems such as in case of beams [27–29], plates [30–32], and shells [33, 34] as well.


## 2. Mechanical Formulation

To formulate the finite element (FE) model of the delaminated beam first the displacement, strain, and stress fields of the 5 parts of a delaminated beam with constant thickness are derived (see Figure1). We begin our formulation with the simplest part regarding the mechanical behaviour, i.e., the delaminated part.Figure 1
Mechanical model and geometric parameters of the delaminated composite beam. ①: left uncracked portion; ②: left delamination tip; ③: delamination portion; ④: right delamination tip; ⑤: right uncracked portion.


### 2.1. Delaminated Part

There are two separate equivalent single layers (ESL) on this part of the delaminated beam, a top and a bottom one which are moving independently of each other. In this case the displacement field of the two can be described with the following:(1)uδTx,yδ=u0δx-yδθδxvδx0,where u0δ is the membrane displacement of the local cross section, θδ is the total cross section rotation, and vδ is the transverse displacement of the beam. δ can take the values t and b if the top or bottom part is taken into consideration, respectively. Using (1) we can formulate the linear and the nonlinear part of the strain tensor which will be later used to determine the elastic strain energy stored in the beam and in the derivation of the geometric stiffness matrix.(2)εδL=12∇⊗uδ+uδ⊗∇(3)εδNL=12∇⊗uδ·uδ⊗∇,where ∇ represent the differential operator whose elements in a 3D Cartesian coordinate system are ∂/∂x, ∂/∂yδ, and ∂/∂z and ⊗ represents the dyadic product operator. Applying the operations the nonzero elements of the linear strain tensor and the nonlinear strain tensor element later to be used are the following:(4)εxδL=u0δ′-yδθδ′γxyδL=vδ′-θδ(5)εxδNL=12vδ′2+u0δ′-yδθδ′2Applying the assumptions of a Timoshenko beam the stress field is in linear relationship with the strain field.(6)σxδk=Eδkεxδτxyδk=Gδkγxyδ,where Ek and Gk are the kth layer’s elastic and shear modulus, respectively. The linear elastic strain energy in a given length of the beam is calculated as follows assuming that the strain and stress values do not depend on the variable z:(7)UδL=12∫lb∫tδσxδkεxδL+τxyδkγxyδLdyδdx,where b is the thickness of the beam.Substituting (4) and (6) into (7) and performing the integration along the yδ axis the following is obtained:(8)UδL=12∫lA11δu0δ′2-2B11δu0δ′θδ′+D11δθδ′2dx+12∫lkA66δvδ′-θδ2dx,where A11δ, B11δ, and D11δ is the membrane, coupling, and bending stiffness, respectively [35]. They are calculated as follows:(9)A11δ,B11δ,D11δ=b∫tδ1,yδ,yδ2EδkdyδA66δ is calculated in a similar manner to A11δ; in the computation the Young modulus of thekth layer is exchanged to the shear modulus of thekth layer.The positions of the top and bottom reference plane can be selected arbitrarily; however it is worth selecting them so that the coupling stiffness vanishes. We will use this approach in the Results but not during the FE equation derivation.The shear correction coefficientk is introduced to compensate the overestimation of the potential energy connected to the shear stress of the first-order theory. Its value for a rectangle cross section is 5/6 [36].The kinetic energy is calculated by taking the time derivative of the displacement vector (1) and multiplying it with itself and with the density of the kth layer.(10)Tδ=12∫lb∫tδρδku˙0δ-yδθ˙δ2+v˙δ2dyδdx(11)Tδ=12∫lI0δu˙0δ2+v˙δ2-2I1δu˙0δθ˙δ+I2δθ˙δ2dx,where I0δ, I1δ, and I2δ is the zeroth, first, and second mass moment of inertia, respectively [37], and are calculated as follows:(12)I0δ,I1δ,I2δ=b∫tδ1,yδ,yδ2ρδkdyδ,where ρδk is the kth layer’s density.During the FE discretization the static equilibrium equation connected toθδ will be used which is determined by applying the Euler-equation [38] to (8). In this case the equilibrium equation takes the following form:(13)δθδ:D11δθδ′′+kA66δvδ′-θδ-B11δu0δ′′=0


### 2.2. Intact Part

On the undelaminated region the beam is cut into two parts by the delamination plane and the continuity of the displacement field on the interface plane is ensured by applying the system of exact kinematic conditions (SEKC) [39]. The general forms of the delaminated part’s and the intact part’s top and bottom ESL’s displacement fields are the same as the exception where in the latter case there is only one transverse displacement.(14)uδTx,yδ=u0δ-yδθδv0 If the delamination plane is above compared to the global reference plane the SEKC takes the following form; refer to Figure 2:(15)ubx,tb-ybR=utx,-ytRubx,yR-ybR=u0xSolving (15) for u0t and u0b then substituting the results into the original expression the axial displacements will have the following form:(16)ut=u0-ytRθt-tb-yRθb-ytθtub=u0-ybR-yRθb-ybθbFigure 2
Displacement field in the vicinity of the left crack tip.Using (16) the strain and stress fields can be derived in a similar manner as for the delaminated part. Moreover the potential and the kinetic energy for the top and bottom ESL is calculated using the same principles described earlier. The total potential and kinetic energy is determined by taking the sum of the top and bottom values. However the static equilibrium equations of the intact part differ significantly; therefore the ones which will have a role in the FE discretization are presented here in detail:(17)δθt:D11tθt′′+kA66tv′-θt′-ytRA11tu0′′-ytRθt′′-tb-yRθb′′-B11tu0′′-2ytRθt′′-tb-yRθb′′=0δθb:D11bθb′′+kA66bv′-θb′-tb-yRA11tu0′′-tb-yRθb′′-ytRθt-B11bu0′′-2ybR-yRθb′′=0


## 2.1. Delaminated Part

There are two separate equivalent single layers (ESL) on this part of the delaminated beam, a top and a bottom one which are moving independently of each other. In this case the displacement field of the two can be described with the following:(1)uδTx,yδ=u0δx-yδθδxvδx0,where u0δ is the membrane displacement of the local cross section, θδ is the total cross section rotation, and vδ is the transverse displacement of the beam. δ can take the values t and b if the top or bottom part is taken into consideration, respectively. Using (1) we can formulate the linear and the nonlinear part of the strain tensor which will be later used to determine the elastic strain energy stored in the beam and in the derivation of the geometric stiffness matrix.(2)εδL=12∇⊗uδ+uδ⊗∇(3)εδNL=12∇⊗uδ·uδ⊗∇,where ∇ represent the differential operator whose elements in a 3D Cartesian coordinate system are ∂/∂x, ∂/∂yδ, and ∂/∂z and ⊗ represents the dyadic product operator. Applying the operations the nonzero elements of the linear strain tensor and the nonlinear strain tensor element later to be used are the following:(4)εxδL=u0δ′-yδθδ′γxyδL=vδ′-θδ(5)εxδNL=12vδ′2+u0δ′-yδθδ′2Applying the assumptions of a Timoshenko beam the stress field is in linear relationship with the strain field.(6)σxδk=Eδkεxδτxyδk=Gδkγxyδ,where Ek and Gk are the kth layer’s elastic and shear modulus, respectively. The linear elastic strain energy in a given length of the beam is calculated as follows assuming that the strain and stress values do not depend on the variable z:(7)UδL=12∫lb∫tδσxδkεxδL+τxyδkγxyδLdyδdx,where b is the thickness of the beam.Substituting (4) and (6) into (7) and performing the integration along the yδ axis the following is obtained:(8)UδL=12∫lA11δu0δ′2-2B11δu0δ′θδ′+D11δθδ′2dx+12∫lkA66δvδ′-θδ2dx,where A11δ, B11δ, and D11δ is the membrane, coupling, and bending stiffness, respectively [35]. They are calculated as follows:(9)A11δ,B11δ,D11δ=b∫tδ1,yδ,yδ2EδkdyδA66δ is calculated in a similar manner to A11δ; in the computation the Young modulus of thekth layer is exchanged to the shear modulus of thekth layer.The positions of the top and bottom reference plane can be selected arbitrarily; however it is worth selecting them so that the coupling stiffness vanishes. We will use this approach in the Results but not during the FE equation derivation.The shear correction coefficientk is introduced to compensate the overestimation of the potential energy connected to the shear stress of the first-order theory. Its value for a rectangle cross section is 5/6 [36].The kinetic energy is calculated by taking the time derivative of the displacement vector (1) and multiplying it with itself and with the density of the kth layer.(10)Tδ=12∫lb∫tδρδku˙0δ-yδθ˙δ2+v˙δ2dyδdx(11)Tδ=12∫lI0δu˙0δ2+v˙δ2-2I1δu˙0δθ˙δ+I2δθ˙δ2dx,where I0δ, I1δ, and I2δ is the zeroth, first, and second mass moment of inertia, respectively [37], and are calculated as follows:(12)I0δ,I1δ,I2δ=b∫tδ1,yδ,yδ2ρδkdyδ,where ρδk is the kth layer’s density.During the FE discretization the static equilibrium equation connected toθδ will be used which is determined by applying the Euler-equation [38] to (8). In this case the equilibrium equation takes the following form:(13)δθδ:D11δθδ′′+kA66δvδ′-θδ-B11δu0δ′′=0


## 2.2. Intact Part

On the undelaminated region the beam is cut into two parts by the delamination plane and the continuity of the displacement field on the interface plane is ensured by applying the system of exact kinematic conditions (SEKC) [39]. The general forms of the delaminated part’s and the intact part’s top and bottom ESL’s displacement fields are the same as the exception where in the latter case there is only one transverse displacement.(14)uδTx,yδ=u0δ-yδθδv0 If the delamination plane is above compared to the global reference plane the SEKC takes the following form; refer to Figure 2:(15)ubx,tb-ybR=utx,-ytRubx,yR-ybR=u0xSolving (15) for u0t and u0b then substituting the results into the original expression the axial displacements will have the following form:(16)ut=u0-ytRθt-tb-yRθb-ytθtub=u0-ybR-yRθb-ybθbFigure 2
Displacement field in the vicinity of the left crack tip.Using (16) the strain and stress fields can be derived in a similar manner as for the delaminated part. Moreover the potential and the kinetic energy for the top and bottom ESL is calculated using the same principles described earlier. The total potential and kinetic energy is determined by taking the sum of the top and bottom values. However the static equilibrium equations of the intact part differ significantly; therefore the ones which will have a role in the FE discretization are presented here in detail:(17)δθt:D11tθt′′+kA66tv′-θt′-ytRA11tu0′′-ytRθt′′-tb-yRθb′′-B11tu0′′-2ytRθt′′-tb-yRθb′′=0δθb:D11bθb′′+kA66bv′-θb′-tb-yRA11tu0′′-tb-yRθb′′-ytRθt-B11bu0′′-2ybR-yRθb′′=0


## 3. Finite Element Formulation

The finite element model is built using 4 different elements. On the intact parts (1st and 5th part) an 8-DOF complex 2D beam element is used; in the left and right crack tips (2nd and 4th part) two similar 10-DOF elements are used which are capable of connecting the top and bottom part of the delamination with the intact part; and finally on the delaminated part a simple 6-DOF beam element is used; refer to Figure3.Figure 3
Finite element model of the delaminated beam using Timoshenko beam theory.


### 3.1. Delaminated Part

As before let us begin the FE element matrix derivation with the delaminated part. In this paragraph let the value ofi+1 on Figure 3 be 1. In this case the nodal DOFs of a general Timoshenko finite element are:(18)ueδT=u1δv1δθ1δu2δv2δθ2δAssuming constant shear strain along one element the displacement and the total cross section rotation are interpolated by the following functions:(19)ueδξ=b0δ+b1δξ(20)veδξ=a0δ+a1δξ+a2δξ2+a3δξ3(21)θeδξ=ddxveδξ-γ0δ=1lea1δ+2a2δξ+3a3δξ2-γ0δ,where le is the element’s length and ξ represents the nondimensional coordinate whose value is selected from the 0,1 domain. Equations (19)–(21) contain seven unknown parameters. The nodal conditions yield six equations and the seventh comes from the equilibrium equation. In (13) the exact displacements are exchanged with the interpolated ones and the global derivation is carried out in the element’s local coordinate system.(22)ueδ0=u1δveδ0=v1δθeδ0=θ1δueδ1=u2δveδ1=v2δθeδ1=θ2δ(23)kle3A66δγ0δ+6D11δa3δ=0Solving (22) and (23) for the interpolation function parameters and substituting them into (19)–(21), the nodal DOFs can be separated and the interpolating functions can be expressed as scalar products.(24)ueδ=Nuδ·ueδ,veδ=Nvδ·ueδ,θeδ=Nθδ·ueδ,γ0δ=Nγδ·ueδSubstituting the interpolated displacement fields found in (24) into the expression of the kinetic and potential energy found in (8)–(11) and changing the integration domain to the nondimensional interval, the aforementioned quantities can be expressed in a general form:(25)Teδ=12u˙eδT∑i=1NMmiδ∫01Miδledξu˙eδ=12u˙eδTMeδu˙eδUeLδ=12ueδT∑i=1NKkiδ∫01Kiδledξueδ=12ueδTKeδueδ,where the fact that the vectors of interpolating functions do not depend on time and Meδ and Keδ represent the element mass and stiffness matrix, respectively, is used.The elements of the summation can be found in Table1. It is worth noticing that in the first three cases the miδ and kiδ coefficients have a similarity which is expected observing (9) and (12). The integrands likewise show some resemblance in the first three cases. Namely, the terms of the dyadic products of Kiδ are obtained by taking the global derivative of the terms found in Miδ. Finally if the terms found on the left and right side of the dyadic product are different then the resulting matrix integrand is found by taking the arithmetic mean of the original and the reversed order dyadic product to obtain the symmetric element matrix.Table 1
Mass and stiffness matrix coefficients and matrix integrands for the delaminated part.
miδ M i δ k i δ Kiδ I0δ N u δ ⊗ N u δ A 11 δ Nuδ′⊗dNuδ′ -2I1δ N u δ ⊗ N θ δ - 2 B 11 δ Nuδ′⊗Nθδ′ I2δ N θ δ ⊗ N θ δ D 11 δ Nθδ′⊗Nθδ′ I0δ N v δ ⊗ N v δ A 66 δ Nγδ⊗NγδThe geometric stiffness matrices are calculated using the principal theory found in [40]. First the elastic energy connected to the nonlinear terms is determined:(26)UδNL=∫lb∫tδσδ0:εδNLdx,where σδ0 is the prestress calculated using only the linear part of the strain tensor and εδNL represents the nonlinear part of the strain tensor. In our case the prestress is the result of an axial force; therefore the initial strain tensor’s only nonzero term is the membrane strain in the x direction εxδ0. Substituting the interpolated fields into (26) the geometric stiffness matrix for one element is formed similarly to the mass and stiffness matrices:(27)UeδNL=12ueδT∑i=1NKGkiGδ∫01Nuδ′·ue0δKiGδledξueδ=12ueδTKeGδueδ,where ue0δ denotes the known nodal displacement vector and the coefficients and matrices of the sum can be found in Table 2.Table 2
Geometric stiffness matrix coefficient and matrix integrands for the delaminated part.
kiGδ KiGδ A11δ Nuδ′⊗dNuδ′ -2B11δ Nuδ′⊗Nθδ′ D11δ Nθδ′⊗Nθδ′ A11δ Nvδ′⊗dNvδ′


### 3.2. Intact Part

In this paragraph let the value ofi-1 in Figure 3 be 1. In this case the 8 DOFs of the FE using the ESL theory are(28)ue=u1v1θ1tθ1bu2v2θ2tθ2bThe interpolation of the axial and transverse displacements of the top and bottom ESL on the intact part is similar to the delaminated portion with the exception of the axial displacement in which case only one function is used to describe the fields found in (19)–(21).The interpolating functions contain 12 parameters which are determined in terms of the nodal parameters by solving the system consisting of the 2 equilibrium equations (17) and the 10 nodal conditions:(29)ue0=u1vet0=v1veb0=v1θet0=θ1tθeb0=θ1bue1=u2vet1=v2veb1=v2θet1=θ2tθeb1=θ2bAs can be seen it is only guaranteed in the elements’ nodes that the top and bottom ESL’s transverse displacement is the same; however this inaccuracy will be negligible in our computations.In order to be able to use the equilibrium equations for one element the exact displacement functions in (17) are swapped by the interpolated ones. In the equation corresponding to the top and bottom part’s equilibria vet and veb are substituted, respectively. The resulting equations for the two ESL are the following:(30)kle3A66δγ0δ+6αδta3t+αδba3b=0,where αδt and αδb are the following in case of the top and bottom ESL:(31)αtt=D11t+ytRytRA11t+2B11tαtb=tb-yRytRA11t+B11tαbt=tb-yRytRA11t+B11tαbb=D11b+tb-yR2A11t+ybR-yRybR-yRA11b+2B11bUsing these results the interpolating functions can be expressed as scalar products of the vector of interpolating functions and the vector of nodal DOFs:(32)ue=Nu·ueveδ=Nvδ·ueθeδ=Nθδ·ueγδ0=Nγδ·ueSubstituting the interpolated functions into the expression of kinetic and potential energy and taking into consideration the fact that the vectors of interpolating functions do not depend on time a previously introduced general summation form is obtained.Six mass and stiffness matrix coefficients and matrix integrand show resemblance similarly to the delaminated part; thus only the quantities connected to the mass matrix are presented in Table3. However there are two mass and stiffness matrices which are not alike; they are presented in Table 4.Table 3
Mass matrix coefficients and integrands for the intact part. Stiffness matrices are similar.
mi Mi I0t+I0b Nu⊗Nu -2ytRI0t+I1t Nu⊗Nθt -2tb-yRI0t+ybR-yRI0b-I1b Nu⊗Nθb ytR2I0t+2ytRI1t+I2t Nθt⊗Nθt 2tb-yRytRI0t+I1t Nθt⊗Nθb t b - y R 2 I 0 t + y b R - y R 2 I 0 b N θ b ⊗ N θ b +2ybR-yRI1b+I2bTable 4
Mass and stiffness matrix coefficients and integrands.
mi M i k i Ki I0t N v t ⊗ N v t k A 66 t Nγt⊗Nγt I0b N v b ⊗ N v b k A 66 b Nγb⊗Nγb


### 3.3. Left Transition Element

Let the value ofi in Figure 3 be 1 in this paragraph. In this case the vectors of nodal parameters for the left transition element are(33)ueδT=u1v1θ1tθ1bu2tv2tθ2tu2bv2bθ2bThe element has a top and bottom part whose displacement and interpolations are the same as those for the delaminated part (see (19)–(21)). According to that there are 14 unknown interpolation parameters. At the left side of the element the nodal conditions make, it possible to connect the intact part’s displacement field with the delaminated parts. At this location the conditions for the membrane displacements are determined by evaluating the first and second equation of (16) at the top and bottom reference planes, respectively (yt=0 and yb=0); the rest of the equations are self-evident:(34)uet0=u1-ytRθ1t-tb-yRθ1bvet0=v1θet0=θ1tueb0=u1-ybR-yRθ1bveb0=v1θeb0=θ1buet1=u2tvet1=v2tθet1=θ2tueb1=u2bveb1=v2bθeb1=θ2bThe unknown parameters are determined in terms of the nodal parameters using the 12 nodal conditions (see (34)) and the static equilibrium equations found in (23). Subsequently the interpolated displacements can be expressed as scalar products:(35)uet=Nut·uevet=Nvt·ueθet=Nθt·ueγt0=Nγt·ueueb=Nub·ueveb=Nvb·ueθeb=Nθb·ueγb0=Nγb·ueSubstituting the interpolated displacements for the top and bottom part into the expression of the potential and kinetic energy and the elastic energy connected to the nonlinear terms in (8), (11), and (26) the mass, stiffness, and geometric stiffness matrices can be formulated. The coefficients and matrices of the summation form can be found in Tables 1 and 2. The resulting element matrices for the transition element are found by taking the sum of the matrices determined for the top and bottom part.


### 3.4. Right Transition Element

The derivation of the right transition element matrices is analogous to the left one. Hence only the bare minimum is presented here. In Figure3 let the value of j-1 be 1. According to that the vector of nodal parameters are(36)ueδT=u1tv1tθ1tu1bv1bθ1bu2v2θ2tθ2bTo determine the vectors of interpolating functions 12 conditions are determined for the nodal values:(37)uet0=u1tvet0=v1tθet0=θ1tueb0=u1bveb0=v1bθeb0=θ1buet1=u2-ytRθ2t-tb-yRθ2bvet1=v2θet1=θ2tueb1=u2-ybR-yRθ2bveb1=v2θeb1=θ2b The equilibrium equations for the top and bottom part are used as well.


### 3.5. Equation of Motion

The equation of motion (EOM) and the structural FE matrices are derived using Hamilton’s principle, which in case of no external work is as follows [41]:(38)δ∫t1t2T-ULdt=0,where T and U represent the total kinetic and potential energy of the system, respectively. Carrying out the variation after the introduction of the FE matrices the EOM of the undamped system is obtained:(39)MU¨+KU=0,where M and K denote N by N structural mass and stiffness matrices, respectively, and U is the vector of nodal displacements.The natural circular frequencies and mode shapes of the undamped system are determined via the generalized eigenvalue-eigenvector problem:(40)-αi2M+Kϕi=0,where αi and ϕi are the ith natural circular frequency and mode shape.In order to take into account the effect of damping as a first approach the proportional damping model is chosen [42]. The model calculates the system’s damping as a linear combination of the mass and stiffness matrix:(41)C=cMM+cKKImplementing the effect of damping in the EOM and transforming it using the modal matrix containing mass matrix orthonormal mode shape vectors,N uncoupled second-order differential equation is obtained:(42)Y¨i+cM+cKαi2Y˙i+αi2Yi=0It is known from the analysis of 1-DOF spring-damper-mass system that the EOM is similar to (42); however the coefficient of the first time derivative is 2Diαi, where Di represents the Rayleigh damping. Using this relationship the Rayleigh damping can be expressed as(43)Di=12cMαi+αicKBased on the principles of experimental modal analysis the Rayleigh damping can be calculated as the function of the peak and notch frequencies of the measured frequency response function around the natural frequency:(44)Di=ω22-ω122ω22+ω12,where ω1 and ω2 correspond to the peak and notch circular frequencies on the real part of the FRF; refer to Figure 4.Figure 4
Real and imaginary part of the FRF close to a natural frequency.Based on the measured FRFs severalαi-Di pairs are determined from which the cM-cK pair is calculated using least squares method. Using these values the damped EOM in case of constant and harmonic conservative force excitation is as follows:(45)MU¨+CU˙+KU=F0+Fssin⁡ωt+Fccos⁡ωtThe solution of (45) is like the excitation; therefore it has the following general form:(46)U=U0+Asin⁡ωt+Bcos⁡ωtBy back-substitution of the solution to the EOM and applying the balance of similar terms the following matrix equation is obtained:(47)K000-ω2M+K-ωC0ωC-ω2M+KU0AB=F0FsFc


### 3.6. Stability Analysis

According to [40] the stability of the whole structure subjected to the previously defined excitation can be described by the following Mathieu-type equation:(48)MU¨+KG0+KGssin⁡ωt+KGccos⁡ωt+KU=0,where KG0, KGs, and KGc are the geometric stiffness matrices calculated using U0, A, and B, respectively.Based on the analysis of the one-DOF system using Floquet-theory the solution of (48) in case of nonzero excitation frequency is dynamically stable; i.e., the spectral range of the Floquet multipliers is exactly one (the amplitude of the vibration stays bounded) or dynamically unstable; i.e., the vibration amplitude is unbounded [43]; this concept was extended for multi-DOF systems as well [26].We seek the solution of (48) on the boundary of stability using two trigonometric function series. One corresponds to the so called T periodic and the other one to the 2T periodic solution [26]:(49)UT=b0+∑k=2,4…aksin⁡kωt2+bkcos⁡kωt2(50)U2T=∑k=1,3…aksin⁡kωt2+bkcos⁡kωt2Substituting (49) and (50) into (48) and applying trigonometric identities we can determine the matrix equation corresponding to the solution of the harmonic balance method. These equations due to their sizes are presented in detail in the Appendix; however their general form is shown here:(51)STcT=0S2Tc2T=0,where the system matrices ST and S2T contain combinations of the structural matrices and the system coefficient vectors cT and c2T contain the coefficients of the trigonometric functions found in (49) and (50): (52)cTT=b0a2b2a4b4⋯c2TT=a1b1a3b3⋯If the effect of damping is negligible and there is no sinusoidal excitation theak and bk coefficients will appear only on the side of sine and cosine functions, respectively; therefore the harmonic balance method results are not two but four system matrices [27]. The results of this approach will be used later as a basis of comparison for the damped case.Our aim is not to determine the solution of these homogeneous matrix equations but to define the parameter combinations in which case the solution is the nontrivial one; i.e., the system matrix’s determinant is zero. These parameter combinations will yield the curves or surfaces which separate the stable and unstable part of the parameter field.The system matrix determinant’s analytical form cannot be determined effectively in case of larger FE discretized mechanical systems such as ours; therefore the stability boundaries has to be determined by means of numerical method. The implemented algorithm was a 2D global root finding algorithm using the concepts of bisection method found in [44].Another problem which arises during the calculation of the system matrix determinant is the use of floating point arithmetic. The magnitude of the determinant is quickly increasing with the number of DOFs of the structure; nonetheless we are only interested in the sign of the determinant in the stability boundary determination algorithm. These are the reasons which lead to the use of the LU matrix decomposition and to use the properties of the factors [45]:(53)sgnS=sgnLU=sgnU=∏i=1NsgnUii,where L is a lower triangular matrix with identity in the diagonal and U is an general upper triangular matrix.If in (45) the excitation frequency is zero only the effect of the static conservative force acting at the free end has to be considered. It can be shown that in this case only divergence type of instability can occur and the critical load can be calculated by reducing the time dependent problem to a time independent one [46–48]. In this case the stability analysis is corresponding to the static stability analysis also known as static buckling analysis.In order to determine the critical load the vector of nodal displacements has to be calculated for a given type of “unit” load; then the geometric stiffness matrix has to be determined based on this displacement field. In this case the equation which describes the buckling phenomena is as follows:(54)K+λiKG0ψi=0,where λi determines the multiplier of the external load vector in which case the structure loses its stability and ψi represents vector of nodal displacements corresponding to the buckled shape.


## 3.1. Delaminated Part

As before let us begin the FE element matrix derivation with the delaminated part. In this paragraph let the value ofi+1 on Figure 3 be 1. In this case the nodal DOFs of a general Timoshenko finite element are:(18)ueδT=u1δv1δθ1δu2δv2δθ2δAssuming constant shear strain along one element the displacement and the total cross section rotation are interpolated by the following functions:(19)ueδξ=b0δ+b1δξ(20)veδξ=a0δ+a1δξ+a2δξ2+a3δξ3(21)θeδξ=ddxveδξ-γ0δ=1lea1δ+2a2δξ+3a3δξ2-γ0δ,where le is the element’s length and ξ represents the nondimensional coordinate whose value is selected from the 0,1 domain. Equations (19)–(21) contain seven unknown parameters. The nodal conditions yield six equations and the seventh comes from the equilibrium equation. In (13) the exact displacements are exchanged with the interpolated ones and the global derivation is carried out in the element’s local coordinate system.(22)ueδ0=u1δveδ0=v1δθeδ0=θ1δueδ1=u2δveδ1=v2δθeδ1=θ2δ(23)kle3A66δγ0δ+6D11δa3δ=0Solving (22) and (23) for the interpolation function parameters and substituting them into (19)–(21), the nodal DOFs can be separated and the interpolating functions can be expressed as scalar products.(24)ueδ=Nuδ·ueδ,veδ=Nvδ·ueδ,θeδ=Nθδ·ueδ,γ0δ=Nγδ·ueδSubstituting the interpolated displacement fields found in (24) into the expression of the kinetic and potential energy found in (8)–(11) and changing the integration domain to the nondimensional interval, the aforementioned quantities can be expressed in a general form:(25)Teδ=12u˙eδT∑i=1NMmiδ∫01Miδledξu˙eδ=12u˙eδTMeδu˙eδUeLδ=12ueδT∑i=1NKkiδ∫01Kiδledξueδ=12ueδTKeδueδ,where the fact that the vectors of interpolating functions do not depend on time and Meδ and Keδ represent the element mass and stiffness matrix, respectively, is used.The elements of the summation can be found in Table1. It is worth noticing that in the first three cases the miδ and kiδ coefficients have a similarity which is expected observing (9) and (12). The integrands likewise show some resemblance in the first three cases. Namely, the terms of the dyadic products of Kiδ are obtained by taking the global derivative of the terms found in Miδ. Finally if the terms found on the left and right side of the dyadic product are different then the resulting matrix integrand is found by taking the arithmetic mean of the original and the reversed order dyadic product to obtain the symmetric element matrix.Table 1
Mass and stiffness matrix coefficients and matrix integrands for the delaminated part.
miδ M i δ k i δ Kiδ I0δ N u δ ⊗ N u δ A 11 δ Nuδ′⊗dNuδ′ -2I1δ N u δ ⊗ N θ δ - 2 B 11 δ Nuδ′⊗Nθδ′ I2δ N θ δ ⊗ N θ δ D 11 δ Nθδ′⊗Nθδ′ I0δ N v δ ⊗ N v δ A 66 δ Nγδ⊗NγδThe geometric stiffness matrices are calculated using the principal theory found in [40]. First the elastic energy connected to the nonlinear terms is determined:(26)UδNL=∫lb∫tδσδ0:εδNLdx,where σδ0 is the prestress calculated using only the linear part of the strain tensor and εδNL represents the nonlinear part of the strain tensor. In our case the prestress is the result of an axial force; therefore the initial strain tensor’s only nonzero term is the membrane strain in the x direction εxδ0. Substituting the interpolated fields into (26) the geometric stiffness matrix for one element is formed similarly to the mass and stiffness matrices:(27)UeδNL=12ueδT∑i=1NKGkiGδ∫01Nuδ′·ue0δKiGδledξueδ=12ueδTKeGδueδ,where ue0δ denotes the known nodal displacement vector and the coefficients and matrices of the sum can be found in Table 2.Table 2
Geometric stiffness matrix coefficient and matrix integrands for the delaminated part.
kiGδ KiGδ A11δ Nuδ′⊗dNuδ′ -2B11δ Nuδ′⊗Nθδ′ D11δ Nθδ′⊗Nθδ′ A11δ Nvδ′⊗dNvδ′


## 3.2. Intact Part

In this paragraph let the value ofi-1 in Figure 3 be 1. In this case the 8 DOFs of the FE using the ESL theory are(28)ue=u1v1θ1tθ1bu2v2θ2tθ2bThe interpolation of the axial and transverse displacements of the top and bottom ESL on the intact part is similar to the delaminated portion with the exception of the axial displacement in which case only one function is used to describe the fields found in (19)–(21).The interpolating functions contain 12 parameters which are determined in terms of the nodal parameters by solving the system consisting of the 2 equilibrium equations (17) and the 10 nodal conditions:(29)ue0=u1vet0=v1veb0=v1θet0=θ1tθeb0=θ1bue1=u2vet1=v2veb1=v2θet1=θ2tθeb1=θ2bAs can be seen it is only guaranteed in the elements’ nodes that the top and bottom ESL’s transverse displacement is the same; however this inaccuracy will be negligible in our computations.In order to be able to use the equilibrium equations for one element the exact displacement functions in (17) are swapped by the interpolated ones. In the equation corresponding to the top and bottom part’s equilibria vet and veb are substituted, respectively. The resulting equations for the two ESL are the following:(30)kle3A66δγ0δ+6αδta3t+αδba3b=0,where αδt and αδb are the following in case of the top and bottom ESL:(31)αtt=D11t+ytRytRA11t+2B11tαtb=tb-yRytRA11t+B11tαbt=tb-yRytRA11t+B11tαbb=D11b+tb-yR2A11t+ybR-yRybR-yRA11b+2B11bUsing these results the interpolating functions can be expressed as scalar products of the vector of interpolating functions and the vector of nodal DOFs:(32)ue=Nu·ueveδ=Nvδ·ueθeδ=Nθδ·ueγδ0=Nγδ·ueSubstituting the interpolated functions into the expression of kinetic and potential energy and taking into consideration the fact that the vectors of interpolating functions do not depend on time a previously introduced general summation form is obtained.Six mass and stiffness matrix coefficients and matrix integrand show resemblance similarly to the delaminated part; thus only the quantities connected to the mass matrix are presented in Table3. However there are two mass and stiffness matrices which are not alike; they are presented in Table 4.Table 3
Mass matrix coefficients and integrands for the intact part. Stiffness matrices are similar.
mi Mi I0t+I0b Nu⊗Nu -2ytRI0t+I1t Nu⊗Nθt -2tb-yRI0t+ybR-yRI0b-I1b Nu⊗Nθb ytR2I0t+2ytRI1t+I2t Nθt⊗Nθt 2tb-yRytRI0t+I1t Nθt⊗Nθb t b - y R 2 I 0 t + y b R - y R 2 I 0 b N θ b ⊗ N θ b +2ybR-yRI1b+I2bTable 4
Mass and stiffness matrix coefficients and integrands.
mi M i k i Ki I0t N v t ⊗ N v t k A 66 t Nγt⊗Nγt I0b N v b ⊗ N v b k A 66 b Nγb⊗Nγb


## 3.3. Left Transition Element

Let the value ofi in Figure 3 be 1 in this paragraph. In this case the vectors of nodal parameters for the left transition element are(33)ueδT=u1v1θ1tθ1bu2tv2tθ2tu2bv2bθ2bThe element has a top and bottom part whose displacement and interpolations are the same as those for the delaminated part (see (19)–(21)). According to that there are 14 unknown interpolation parameters. At the left side of the element the nodal conditions make, it possible to connect the intact part’s displacement field with the delaminated parts. At this location the conditions for the membrane displacements are determined by evaluating the first and second equation of (16) at the top and bottom reference planes, respectively (yt=0 and yb=0); the rest of the equations are self-evident:(34)uet0=u1-ytRθ1t-tb-yRθ1bvet0=v1θet0=θ1tueb0=u1-ybR-yRθ1bveb0=v1θeb0=θ1buet1=u2tvet1=v2tθet1=θ2tueb1=u2bveb1=v2bθeb1=θ2bThe unknown parameters are determined in terms of the nodal parameters using the 12 nodal conditions (see (34)) and the static equilibrium equations found in (23). Subsequently the interpolated displacements can be expressed as scalar products:(35)uet=Nut·uevet=Nvt·ueθet=Nθt·ueγt0=Nγt·ueueb=Nub·ueveb=Nvb·ueθeb=Nθb·ueγb0=Nγb·ueSubstituting the interpolated displacements for the top and bottom part into the expression of the potential and kinetic energy and the elastic energy connected to the nonlinear terms in (8), (11), and (26) the mass, stiffness, and geometric stiffness matrices can be formulated. The coefficients and matrices of the summation form can be found in Tables 1 and 2. The resulting element matrices for the transition element are found by taking the sum of the matrices determined for the top and bottom part.


## 3.4. Right Transition Element

The derivation of the right transition element matrices is analogous to the left one. Hence only the bare minimum is presented here. In Figure3 let the value of j-1 be 1. According to that the vector of nodal parameters are(36)ueδT=u1tv1tθ1tu1bv1bθ1bu2v2θ2tθ2bTo determine the vectors of interpolating functions 12 conditions are determined for the nodal values:(37)uet0=u1tvet0=v1tθet0=θ1tueb0=u1bveb0=v1bθeb0=θ1buet1=u2-ytRθ2t-tb-yRθ2bvet1=v2θet1=θ2tueb1=u2-ybR-yRθ2bveb1=v2θeb1=θ2b The equilibrium equations for the top and bottom part are used as well.


## 3.5. Equation of Motion

The equation of motion (EOM) and the structural FE matrices are derived using Hamilton’s principle, which in case of no external work is as follows [41]:(38)δ∫t1t2T-ULdt=0,where T and U represent the total kinetic and potential energy of the system, respectively. Carrying out the variation after the introduction of the FE matrices the EOM of the undamped system is obtained:(39)MU¨+KU=0,where M and K denote N by N structural mass and stiffness matrices, respectively, and U is the vector of nodal displacements.The natural circular frequencies and mode shapes of the undamped system are determined via the generalized eigenvalue-eigenvector problem:(40)-αi2M+Kϕi=0,where αi and ϕi are the ith natural circular frequency and mode shape.In order to take into account the effect of damping as a first approach the proportional damping model is chosen [42]. The model calculates the system’s damping as a linear combination of the mass and stiffness matrix:(41)C=cMM+cKKImplementing the effect of damping in the EOM and transforming it using the modal matrix containing mass matrix orthonormal mode shape vectors,N uncoupled second-order differential equation is obtained:(42)Y¨i+cM+cKαi2Y˙i+αi2Yi=0It is known from the analysis of 1-DOF spring-damper-mass system that the EOM is similar to (42); however the coefficient of the first time derivative is 2Diαi, where Di represents the Rayleigh damping. Using this relationship the Rayleigh damping can be expressed as(43)Di=12cMαi+αicKBased on the principles of experimental modal analysis the Rayleigh damping can be calculated as the function of the peak and notch frequencies of the measured frequency response function around the natural frequency:(44)Di=ω22-ω122ω22+ω12,where ω1 and ω2 correspond to the peak and notch circular frequencies on the real part of the FRF; refer to Figure 4.Figure 4
Real and imaginary part of the FRF close to a natural frequency.Based on the measured FRFs severalαi-Di pairs are determined from which the cM-cK pair is calculated using least squares method. Using these values the damped EOM in case of constant and harmonic conservative force excitation is as follows:(45)MU¨+CU˙+KU=F0+Fssin⁡ωt+Fccos⁡ωtThe solution of (45) is like the excitation; therefore it has the following general form:(46)U=U0+Asin⁡ωt+Bcos⁡ωtBy back-substitution of the solution to the EOM and applying the balance of similar terms the following matrix equation is obtained:(47)K000-ω2M+K-ωC0ωC-ω2M+KU0AB=F0FsFc


## 3.6. Stability Analysis

According to [40] the stability of the whole structure subjected to the previously defined excitation can be described by the following Mathieu-type equation:(48)MU¨+KG0+KGssin⁡ωt+KGccos⁡ωt+KU=0,where KG0, KGs, and KGc are the geometric stiffness matrices calculated using U0, A, and B, respectively.Based on the analysis of the one-DOF system using Floquet-theory the solution of (48) in case of nonzero excitation frequency is dynamically stable; i.e., the spectral range of the Floquet multipliers is exactly one (the amplitude of the vibration stays bounded) or dynamically unstable; i.e., the vibration amplitude is unbounded [43]; this concept was extended for multi-DOF systems as well [26].We seek the solution of (48) on the boundary of stability using two trigonometric function series. One corresponds to the so called T periodic and the other one to the 2T periodic solution [26]:(49)UT=b0+∑k=2,4…aksin⁡kωt2+bkcos⁡kωt2(50)U2T=∑k=1,3…aksin⁡kωt2+bkcos⁡kωt2Substituting (49) and (50) into (48) and applying trigonometric identities we can determine the matrix equation corresponding to the solution of the harmonic balance method. These equations due to their sizes are presented in detail in the Appendix; however their general form is shown here:(51)STcT=0S2Tc2T=0,where the system matrices ST and S2T contain combinations of the structural matrices and the system coefficient vectors cT and c2T contain the coefficients of the trigonometric functions found in (49) and (50): (52)cTT=b0a2b2a4b4⋯c2TT=a1b1a3b3⋯If the effect of damping is negligible and there is no sinusoidal excitation theak and bk coefficients will appear only on the side of sine and cosine functions, respectively; therefore the harmonic balance method results are not two but four system matrices [27]. The results of this approach will be used later as a basis of comparison for the damped case.Our aim is not to determine the solution of these homogeneous matrix equations but to define the parameter combinations in which case the solution is the nontrivial one; i.e., the system matrix’s determinant is zero. These parameter combinations will yield the curves or surfaces which separate the stable and unstable part of the parameter field.The system matrix determinant’s analytical form cannot be determined effectively in case of larger FE discretized mechanical systems such as ours; therefore the stability boundaries has to be determined by means of numerical method. The implemented algorithm was a 2D global root finding algorithm using the concepts of bisection method found in [44].Another problem which arises during the calculation of the system matrix determinant is the use of floating point arithmetic. The magnitude of the determinant is quickly increasing with the number of DOFs of the structure; nonetheless we are only interested in the sign of the determinant in the stability boundary determination algorithm. These are the reasons which lead to the use of the LU matrix decomposition and to use the properties of the factors [45]:(53)sgnS=sgnLU=sgnU=∏i=1NsgnUii,where L is a lower triangular matrix with identity in the diagonal and U is an general upper triangular matrix.If in (45) the excitation frequency is zero only the effect of the static conservative force acting at the free end has to be considered. It can be shown that in this case only divergence type of instability can occur and the critical load can be calculated by reducing the time dependent problem to a time independent one [46–48]. In this case the stability analysis is corresponding to the static stability analysis also known as static buckling analysis.In order to determine the critical load the vector of nodal displacements has to be calculated for a given type of “unit” load; then the geometric stiffness matrix has to be determined based on this displacement field. In this case the equation which describes the buckling phenomena is as follows:(54)K+λiKG0ψi=0,where λi determines the multiplier of the external load vector in which case the structure loses its stability and ψi represents vector of nodal displacements corresponding to the buckled shape.


## 4. Results

In our numerical investigations a unidirectional fibre reinforced composite beam is considered with a single through width delaminated and dense ply structure; therefore an orthotropic material model is selected which can be described by 3 material properties. These values and the ones describing the geometry of the beam and the delamination can be found in Table5.Table 5
Common material and geometric parameter values.
E 11 [GPa] G [GPa] ρ  kg/m3 L [mm] b [mm] t [mm] tb/tt1 33 3 1330 180 20 6.2 13Let us select the local reference planes so that the coupling stiffness vanishes. Due to the homogeneous material properties some expressions get simplified:(55)yδR=tδ2yR=t2The FE model is built using elements with an average length of 5 mm. To eliminate the effect of the free mode model, i.e., the unacceptably large decrease in the natural frequencies and critical forces, the constrained mode model assumption is used which is achieved by setting the transverse displacement of the delaminated part to be the same on the top and bottom part.Moreover only small amplitude vibrations are considered, which is sufficient in case of linear stability analysis; therefore the pressure appearing on the delamination surfaces remains small compared to the external load; thus the effect of frictional force between the top and bottom part can be neglected.As a first step our FE model is validated with a model using higher order plate theory derived by Radu and Chattopadhyay in [31]. For comparison reasons the results determined using Euler-Bernoulli beam theory found in [49] are also given in Table 6. These natural frequencies and critical forces correspond to a thin beam with a single through width delamination found in the symmetry plane. It is seen that taking into account the effect of transverse shear deformation has an advantageous effect on the determined values. Compared to the EB beam theory every value is smaller by 5% approximately; however the difference between the HOT and the Timoshenko beam theory still remains above 10%.Table 6
Comparison of the Euler-Bernoulli (EB) and Timoshenko (TIM) beam theory with the higher order plate theory (HOT) for various delamination lengths.L=127 [mm], b=12.7 [mm], and h=1.016 [mm].
a / L [%] f N 1 [Hz] f N 2 [Hz] F c r 1 [N] EB TIM HOT EB TIM HOT EB TIM HOT 0 97.06 93.22 82.12 608.3 582.8 513.3 22.82 22.81 16.34 20 96.54 92.73 81.19 607.3 581.9 509.2 22.76 21.61 16.07 40 93.20 89.51 76.48 592.2 567.6 469.0 21.89 20.90 15.05 60 85.89 82.47 67.26 532.5 510.6 369.1 19.87 19.18 12.71 80 75.98 72.98 56.64 439.2 421.5 325.8 16.94 16.50 9.934It is known that there is friction between the top and bottom part of the delaminated part during static and dynamic loading which has a positive effect on the damping of the structure [50]. However in our model we do not want to implement this effect, only the structure’s inner material damping. Therefore in our modal analysis the test specimen was an intact beam with the same material properties and total beam length as the delaminated one. In the measurement setup the beam with fixed-free end was excited with a Brüel and Kjær type impact hammer in points A and B and the response of the beam end was measured with a single axis PCB type accelerometer in points C and D (see Figure 5). This way it was possible to separate the natural frequencies corresponding to the bending and torsional modes of which we are only interested in the first one. The averages of 5 measurements can be seen in Figure 6 where the first letter in the subscript represents the position where the beam was excited and the second one the point of measured acceleration signal. Using the imaginary parts of the four transfer functions the frequencies were selected at which every curve had distinct local maxima. Based on these measured natural frequencies the nearest local minima and maxima on the real part, i.e., the peak and notch frequency values, were also determined. Table 7 contains the measured natural frequencies and calculated Rayleigh damping using (44) corresponding to bending vibrations.Table 7
Measured natural frequencies and calculated Rayleigh damping values corresponding to bending vibrations.
TRF i 1 2 3 4 HAC f N i   H z 133.5 840 2281 4337 103Di1 13.18 6.839 8.003 16.18 HAD f N i   H z 133.5 842 2282 4340 103Di1 13.18 6.532 7.998 17.05 HBC f N i   H z 133 847 2301 - 103Di1 7.518 21.92 6.519 - HBD f N i   H z 133 844 2293 4382 103Di1 7.518 23.22 5.460 15.58Figure 5
Damping measurement setup. A and B: impact hammer excitation points; C and D: accelerometer positions.Figure 6
Transfer functions of the intact beam.To obtain the proportional damping coefficients the least squares method (LSM) was used on the average of each data set corresponding to the same natural frequency to achieve uniform weighting. It can be seen in case ofi=2 in Table 7 that the calculated Rayleigh damping from two groups based on the excitation point. Due to that two fitted curves were calculated; one takes into consideration the outlying points and the other does not. The resulting fitted parameters are found in Table 8 and the best fitting curves with the points are presented in Figure 7. The difference between the fitted model parameters is less than 10%; moreover it is seen in Figure 7 that the effect of the outlying points is not significant in the 0–500 Hz range, where the stability analysis will be carried out. In the subsequent paragraphs the parameter which results in a lower damping is selected because in this case the global stability boundaries will be more reliable.Table 8
Proportional damping coefficients determined from the measurement data with and without the outlying points using least squares method.
c M   1 / s c K   s with 18.78 1.208⋅10-6 without 17.02 1.100⋅10-6Figure 7
Measurement results for the natural frequency-Rayleigh damping pairs and fitted proportional damping model with and without the outlying points.◊: HAC; ×: HAD; □: HBC; +: HBD.For comparison reasons in Table9 the first four natural frequencies of the undamped system and first critical load are presented for different delamination lengths; the intact parts’ size can be found in [49]. Similar to before, the higher order beam theory results in lower natural frequencies and critical forces compared to the EB beam theory. The stability maps are determined for the case when the external force excitation acts at the free end along the beam axis and has the following general form:(56)Ft=F0+Fssin⁡ωt+Fccos⁡ωtThe resulting stability maps for the aforementioned excitation are not presented here for every delamination lengths because the difference is not significant in every case. Therefore only the case of the intact beam and the beam with the longest delamination length is presented (see Figures 8 and 9). In both figures the left and right boundaries of the unstable parameter regions (gray areas) were determined from the undamped EOM when four sets of system matrix are obtained through the harmonic balance method and the black curves are corresponding to the damped system; also the matrices were truncated after the second-order terms. Moreover if the sign of F0 is negative or positive compressive or tensile static load is acting on the beam, respectively. The complete stability chart is obtained by superimposing the boundary curves determined using the T and 2T periodic system matrices. The unstable and stable regions of the maps are intuitive, but numerical time integration using the Newmark-beta algorithm [51] was carried out in several points of the chart; the beam end’s transverse deflection time history is shown in Figure 10. It can be seen that in case of points “a” and “c” the solution remains bounded; however in case of points “b” and “d” the solution increases exponentially; therefore the beam is in unstable state; thus the dynamically unstable regions are the V shaped regions on the parameter field.Table 9
Natural frequencies and critical forces of the intact and delaminated beams using Timoshenko beam theory.
a [mm] f N 1 [Hz] f N 2 [Hz] f N 3 [Hz] f N 4 [Hz] F c r 1 [N] 0 153.5 944.1 2573 4878 994.8 60 153.1 941.6 2484 4764 992.0 80 152.7 938.2 2432 4670 989.6 100 151.9 929.2 2396 4570 986.2 120 150.9 918.2 2370 4504 980.6Figure 8
Stability boundaries of the proportionally damped (black curves) and unstable regions of the undamped (gray regions) intact beam for different static load values. (a)F0=0.5Fcr1, (b) F0=0Fcr1, (c) F0=-0.5Fcr1, a=0 [mm], and Fs=0 [N].
(a) (b) (c)Figure 9
Stability boundaries of the proportionally damped (black curves) and unstable regions of the undamped (gray regions) delaminated beam for different static load values. (a)F0=0.5Fcr1, (b) F0=0Fcr1, (c) F0=-0.5Fcr1, a=120 [mm], and Fs=0 [N].
(a) (b) (c)Figure 10
Time history of an intact beam’s end deflection in case of different excitation frequency and nondimensional force amplitudes: (a) (150, 0.25), (b) (150, 0.75), (c) (275, 0.25), and (d) (275, 0.75).
(a) (b) (c) (d)Looking individually at Figures8 and 9 by applying tensional static force on the beams the unstable regions become thinner which is expected because this kind of load works against the buckling phenomenon; as well the tip of the unstable regions are shifted towards higher frequencies; in both aspects the opposite can be said in case of compressive force. The effect of damping on the stability region is almost negligible in case of the first unstable region connected to the 2T type stability loss (rightmost unstable region). However for smaller frequencies this effect becomes more and more remarkable and the tip of the regions shifts towards higher dimensionless force amplitudes.One important aspect to notice is that the stability boundary finding algorithm based on the bisection method is highly dependent on the minimal element size of the stability map. In case of narrow crevice like boundary curve parts the method is not able to work properly. Therefore it is possible that the calculated lowest point of the damped system’s stability region is not the real one. This phenomenon is shown in Figure11 on which the same region is shown for different resolution level. Another problem which arises from the use of finite sized elements is that some stability regions will remain separated even though they correspond to the same unstable region (see Figure 8(c)  f ≈ 75 Hz). Despite these flaws the algorithm is still considered to be a fast and accurate method owing to the fact that far from the stability boundaries the resolution of the mesh can be several orders lower compared to the vicinity its and the required precision can be selected in advance.Figure 11
Effect of resolution level on the unstable region for the intact proportionally damped beam. The levels of resolution for the figures are (a) 1, (b) 3, (c) 5, and (d) 7.
(a) (b) (c) (d)Comparing the stability charts obtained from the two structures it can be observed that the existence of delamination does not modify significantly the shape of the unstable regions but rather shifts them along the frequency axis with a value proportional to the difference of two beams’ natural frequencies.


## 5. Conclusions

In this paper the global dynamic stability of laminated composite beam with a single through width delamination was investigated using the finite element method. During the mechanical formulation of the equivalent single layer theory and for the connection of the intact and delaminated part the system of exact kinematic conditions was used. The mechanical model and the element matrices were derived using the first-order shear deformable theory also known as Timoshenko’s beam theory with constant shear strain along one element and with the application of the static equilibrium equation. It was found that the higher order theory gives a better model of the cracked beam compared to the Euler-Bernoulli theory from the viewpoint of vibration and static and dynamic stability; however plate theories are still superior in this sense.The stability maps were determined assuming mass and stiffness proportional structural damping for which the coefficients were calculated using experimental data obtained from modal analysis. The stability of the equation of motion was determined using Bolotin’s harmonic balance method, the LU matrix decomposition, and the bisection method for theT and 2T periodic system response. The effects of damping and static force acting along the beam axis on the stability maps were investigated as well. It was found that the damping determined for this structure does not decrease significantly the unstable regions; on the other hand the compressional force decreases while tensional force increases the size of the stable regions. Comparing the maps corresponding to the intact and the delaminated beam it was observed that the delamination does not change the shape of the unstable domains; rather it shifts them along the frequency axis.



---

*Source: 2674813-2018-06-06.xml*