# A Finite Element Formulation for Bending-Torsion Coupled Vibration Analysis of Delaminated Beams under Combined Axial Load and End Moment

Authors: Mir Tahmaseb Kashani; Seyed M. Hashemi
Journal: Shock and Vibration
 (2018)

---

## Abstract

Free vibration analysis of beams with single delamination undergoing bending-torsion coupling is made, using traditional finite element technique. The Galerkin weighted residual method is applied to convert the coupled differential equations of motion into to a discrete problem, where, in addition to the conventional mass and stiffness matrices, a delamination stiffness matrix, representing the extra stiffening effects at the delamination tips, is introduced. The linear eigenvalue problem resulting from the discretization along the length of the beam is solved to determine the frequencies and modes of free vibration. Both “free mode” and “constrained mode” delamination models are considered in formulation, and it is shown that the continuity (both kinematic and force) conditions at the beam span-wise locations corresponding to the extremities of the delaminated region, in particular, play a great role in “free mode” model formulation. Current trends in the literature are examined, and insight into different types of modeling techniques and constraint types are introduced. In addition, the data previously available in the literature and those obtained from a finite element-based commercial software are utilized to validate the presented modeling scheme and to verify the correctness of natural frequencies of the systems analyzed here. The paper ends with general discussions and conclusions on the presented theories and modeling approaches.


---

## Body

## 1. Introduction

Materials or components built up in layers are used in various manufacturing sectors ranging from shipbuilding to aerospace. The ever-increasing application of such layered structural elements is primarily due to their many attractive features such as high specific stiffness, high specific strength, good buckling resistance, and formability into complex shapes, to name a few. Delamination is the most common failure mode in layered structural elements and may be introduced during processing or subsequent service conditions, leading to the reduction or loss of structural strength and stiffness [1], which, in turn, will modify the vibration characteristics of the structure, i.e., natural frequency and mode shape. Natural frequencies will decrease caused by the stiffness reduction, which may subsequently lead to resonance, should the reduced frequency match the working frequency. Therefore, it is imperative to understand the effect of the delamination and its influence on the vibration characteristics of structures. Depending on the number of delamination, size, and location, the defective sublaminate generally displays different vibration frequencies and modes.One of the earliest models for vibration analysis of composite beams with delaminations was proposed by Kulkarni and Fredericks [2] where they studied a one-dimensional (1D) problem with a single delamination. They investigated a defective, circumferentially symmetric (circular) cylindrical shell, with a small length delamination located at the middle surface. In their analysis, it was assumed that the flexural rigidity of the delamination region can be taken as the sum of the flexural rigidities of the all individual delaminated layers. This analysis resulted in system natural frequencies considerably lower than those found experimentally. Ramkumar et al. [3] investigated the vibration characteristics of defective composite beams, where a single through-the-width delamination was modeled as four Timoshenko beams connected at delamination edges. The predicted natural frequencies were consistently found to be lower than those measured experimentally. This discrepancy was attributed to the contact between the delaminated free surfaces during vibrations, and the authors suggested that the inclusion of the contact effect might improve the analytical prediction. However, it was later discovered that underpredicted natural frequencies for off midplane delaminations, resulting from the free mode model, are in fact due to unrestricted penetration of the beams into each other. Considering free vibration characteristics of isotropic defective beams with through-the-width split, Wang et al. [4] developed an improved theoretical model to investigate the effect of delamination on the system’s frequencies. Assuming each beam segment, modeled as an Euler beam, vibrates freely without touching each other, they included the coupling between flexural and axial vibrations of the delaminated sublaminates in their model. Although this improved analytical model resulted in frequency data closer to experimental results, it predicted a few physically inadmissible opening modes, addressed by Mujumdar and Suryanarayan [5]. The “constrained mode” model [5] assumes that the delaminated beam segments undergo identical transverse displacements and predicts vibration behaviour more accurately for off midplane delamination configurations. However, opening delamination modes, commonly seen in experimental analysis [6, 7], cannot be captured using the constrained model. Tracy and Pardoen [8] also used the constrained model to assess the effect of delamination on natural frequency of symmetric laminated Euler beams containing midplane delamination. Based on both the classical and high-order shear deformation beam theories, and using the FEM, Nagesh Babu and Hanagud [9] studied the effects of delamination on the system’s natural frequencies. Performing experiments, Shen and Grady [7] observed the presence of opening modes, undetected by the constrained model. By exploiting the Galerkin method, they also investigated the effects of delamination on the natural frequency and mode shape of composite laminated beams. Introducing coupling between longitudinal and bending motion, Stamos et al. [10] presented a delaminated composite beam model and furthermore suggested an inverse method to determine position and size of delamination based on the degradation of the first two natural frequencies. Lee et al. [11] analyzed the free vibration of multi-delaminated composite beam-columns subjected to axial compression load analytically and experimentally. In their work, the characteristic equation of the defective system is obtained by dividing the multi-delaminated beam-columns into segments and by imposing recurrence relation from the continuity conditions on each subbeam-column. To verify their results, experimental results are obtained for isotropic single delaminated beam-columns and confirmed that the sizes, location, and number of delaminations have significant effect on the system’s vibration and stability. Chen et al. [12] developed an analytical model for free vibration of a delaminated layered composite beam in the prebuckled state. They also presented a new constrained model, including both compressive axial force and bending-extension coupling effects. Saravanos and Hopkins [13] investigated the damped-free vibration of unstressed defective laminated composite beams, both analytically and experimentally, where they considered the kinematic perturbations induced by a delamination crack as additional degrees of freedom. In a later work, Chrysochoidis and Saravanos [14] assessed the effects of delamination on the dynamic response of composite beam specimens with surface attached piezoelectric actuators and sensors using experiments and same analytical method developed earlier in [13]. Later, they developed an analytical model using coupled linear layerwise laminate theory and a FE model to investigate both static and dynamic response of flexible defective composite beams, including delaminations and active piezoelectric sensors [15]. Della and Shu [16, 17] reported an analytical solution method, and more recently, Erdelyi and Hashemi used Dynamic Stiffness Matrix (DSM) [18] and presented a FEM-based novel assembly technique [19] to investigate the behaviour of a delaminated slender beam and compared their results with those reported by Wang et al. [4] and Lee [20]. In a recent work, Szekrényes [21] performed stability analysis on delaminated composite beams undergoing coupled flexural-longitudinal vibration. Starting with Timoshenko beam theory and then reducing the model to Euler–Bernoulli solution, they presented a formulation, considering both free and constrained models. In addition, the equality of axial forces in the top and bottom beams was derived and shown in an exact way. Also, the continuity of the effective bending moments was related to the equilibrium equations, and it was also concluded and reported that if the normal force is compressive in one of the half-periods of the vibration and reaches a critical value, then it can lead to delamination buckling [21–23].In two recent publications, the authors exploited the classical [24] and frequency-dependent (dynamic) FEM [25] formulations for the free vibration of intact prestressed beams. The free vibration/buckling of intact and delaminated axially-loaded beams has also been investigated before (e.g., [17–23]). However, to the authors’ best knowledge, the vibration of prestressed delaminated beams, subjected to combined axial load and end moment, has not been reported in the open literature.In the recent years, layered, sandwich and composite elements have found their place and have been integrated in some commercial software, such as ANSYS®, used to perform static and dynamic structural analysis. However, it is worth noting that modeling of delaminated configurations in such environments in not a straightforward procedure can involve cumbersome, complex, time-consuming, and error-prone processes. It requires manual model creation, involving the use of, for example, Multipoint Constraint Rigid Link (ANSYS element MPC184) to enforce the displacement and slope continuity at the edges of delamination region [26]. Therefore, in the present study, the free vibration behaviour of a delaminated prestressed beam subjected to combined axial force and end moment is investigated using the finite element method (FEM). By implementing the Galerkin weighted residual method, the governing differential equations of motion are discretized, leading to element matrices. As it will be explained later in this paper, in addition to the mass and stiffness matrices commonly resulting from the conventional FEM formulation, the presented procedure also leads to three extra stiffness matrices: a geometric stiffness matrix caused by the axial load, a delamination stiffness matrix, resulting from the relationships between the two delamination region’s extremities, and a coupling geometric stiffness matrix caused by the applied end moment. The assembly of the element matrices and application of relevant boundary conditions then leads to the system’s linear eigenvalue problem. The resulting eigenproblem is then solved to extract the system’s natural frequencies and modes. The correctness of the proposed FEM method is first validated, where the free vibration analysis of a simple delaminated beam (in the absence of prestress) is conducted, and the frequency results are compared with those presented in the literature. The application of the presented model and formulation is then further demonstrated through the illustrative examples of unstressed delaminated [09/0]2s composite beams, as well as bending-torsion coupled beams, subjected to various prestress load combinations (i.e., axial forces and end moments). The paper concludes with a brief discussion on the presented method and numerical results. While the FEM model presented in this paper assumes isotropic materials, further research is underway to extend it to sandwich and fibre-reinforced laminated composite beams, characterized by an extensional response coupled with flexural/torsional and coupled bending-torsion vibration.


## 2. Mathematical Model

Figure1 shows a prestressed two-layered beam of length L and thickness H1 having a single through-the-width delamination with a length L2, starting from x=x1. This delamination divides the beam into four segments, which are assumed to have a slenderness ratio of greater than 10 and, as a result will be analyzed as four interconnected Euler–Bernoulli beams. Beam segments 2 and 3 have H2 and H3 thicknesses, respectively. The structure is subjected to an axial force of P0 and constant end moment of Mzz, applied at its ends, x=0 and x=L.Figure 1
Schematic illustration of a beam with single delamination under axial load and end moment.Assuming that each uniform beam segment is made of linearly elastic, homogeneous, isotropic material with constant mechanical, material, and geometric properties, it can be shown that the partial differential equations of motion, governing free, linear, bending-torsion coupled vibrations of theith Euler–Bernoulli beam segment subjected to axial force and end moment are as follows [24]:(1)EIiw″″i+Pi0wi″+Mzz,iϕi″−ρAiw¨i=0,i=1−4,(2)GJiϕi″+Pi0IP,iϕi″+Mzz,iwi″−ρAiIP,iϕ¨i=0,i=1−4,where w is the lateral displacement along z direction, ϕ stands for the torsional twist about x-axis, E is Young’s modulus, G represents the shear modulus, J and Ip are, respectively, the torsion constant and the polar moment of inertia, and subscript i represents the beam segment’s number, where i=1,2,3,4. Pi0 and Mzz,i represent axial force and end moment for ith beam, respectively. A prime, (′), denotes spatial derivative with respect to x, and the over dot, (•), stands for time derivative. The beam’s torsional rigidity (GJ) is assumed to be very large compared with its warping rigidity (EΓ), and ends are free to warp, i.e., state of uniform torsion.As can be observed from Equations (1) and (2), the system is coupled by the end moment, Mzz,i. Exploiting the simple harmonic motion assumption to separate variables in space and time, and introducing the (sinusoidal) functions wix,t=Wisinωt and ϕix,t=θisinωt for lateral and torsional displacements, respectively, in expressions (1) and (2) leads to following form of the governing differential equations, in frequency domain:(3)EIiW″″i+Pi0Wi″+Mzz,iθi″−ρAiω2Wi=0,i=1−4,(4)GJiθi″+Pi0IP,iθi″+Mzz,iWi″−ρIP,iAiω2θi=0,i=1−4,where ω denotes the frequency and Wi and θi are the amplitudes of flexural and torsional displacements, respectively, of the beam segment i.


### 2.1. Free Mode Delamination Model

In contrast to the “constraint mode” model, the “free mode” delamination scheme considers no constraint or interaction between the upper and lowest layers in the defective zone, i.e., segments 2 and 3. The “free mode” model predicts the behaviour of the beam more accurately but it also exhibits some physically impossible modes of vibration. To solve the four pairs of coupled equations as one system numerically, the global end conditions at both extremities of the beam, as well as continuity requirements at the delamination tips, are to be enforced. The continuity conditions for lateral displacement, slope, torsional displacement, and axial force at the left side of the defective region (x=x1) are expressed as(5)W11=W21=W31,W11=W1x=x1,W21=W2x=x1,W31=W3x=x1,(6)W11′=W21′=W31′,W11′=W1′x=x1,W21′=W2′x=x1,W31′=W3′x=x1,(7)θ11=θ21=θ31,θ11=θ1x=x1,θ21=θ2x=x1,θ31=θ3x=x1,(8)P1=P2+P3.Similarly, continuity conditions for the right delamination tip are written as(9)W42=W22=W32,W42=W4x=x2,W22=W2x=x2,W32=W3x=x2,W42′=W22′=W32′,W42′=W4′x=x2,W22′=W2′x=x2,W32′=W3′x=x2,θ42=θ22=θ32,θ42=θ4x=x2,θ22=θ2x=x2,θ32=θ3x=x2,P4=P2+P3.Additionally, the requirement for the delamination tip faces to remain planar after deformation, at the left delamination tip, results in(10)u21−u31=H12W11′,where ui is the axial displacement of beam section i. Combining expression (10) with the similar expression from the right delamination tip leads to(11)Δu3−Δu2=H12W42′−W11′,where(12)Δu3=u32−u31,u32=u3x=x2,u31=u3x=x1,Δu2=u22−u21,u22=u2x=x2,u21=u2x=x1,W42′=W4′x=x2,W11′=W1′x=x1.It was originally assumed by Mujumdar and Suryanarayan [5] and was later shown in an exact way by Szekrényes [21] that for small deformations of a uniform beam, the axial displacement will behave according to the following expression (13):(13)uixi=Li−uixi=0=∫0LiFixiEiAixidxi=FiLiEiAi,where xi represents the local position along the length and Fi is the axial force induced in the segment i. Compressive force during vibration, Px,t, consists of applied static axial force, P0, and perturbed axial force, Ptx,t. The latter, preventing interlaminar slip, is induced in the delaminated layers during vibration. Knowing that the change in axial displacement for each segment is caused by induced perturbed component of axial force rather than its static component, substituting Equation (13) into (11) yields(14)P2tL2EA2−P3tL3EA3=H12W4′x=x2−W1′x=x1.Taking the perturbation of Equation (8) for continuity of axial force on the left delamination tip, one obtains(15)P2t+P3t=P1t.SinceP1=P0 and P0 is constant, P1t=0 which means: P2t+P3t=0. Substituting P3t=−P2t into Equation (14) results in(16)P3t=Λ∗W4′x=x2−W1′x=x1,where the coefficient Λ∗ is defined as(17)Λ∗=H12L2EA2EA3EA3+EA2.At the left and right delamination tips, respectively, continuity of bending moments leads to(18)M11=M21+M31−P2tH32+P3tH22,(19)M42=M22+M32−P2tH32+P3tH22,where Mij=Mix=xj. The internal bending moment, Mx, shear force, Sx, and torque, Tx, in each of the individual beam segments are written as(20)Mix=−EIiWi″,(21)Six=EIiWi‴+Mzz,iθi′+Pi0Wi′,(22)Tix=GJiθi′+Pi0IP,iθi′+Mzz,iWi′.Using Equations (20) and (18), the following continuity of bending moments is obtained(23)EI1W1″x=x1=EI2W2″x=x1+EI3W3″x=x1+ΛW4′x4=0−W1′x1=L1,where the parameter Λ is defined as(24)Λ=H124L2EA2EA3EA3+EA2.Likewise, to satisfy the continuity of shear forces at the left delamination tip, one should have(25)EI1W1‴+P10W1′+Mzz,1θ1′x=x1=EI2W2‴+P20W2′+Mzz,2θ2′x=x1+EI3W3‴+P30W3′+Mzz,3θ3′x=x1,and finally, the continuity of twisting moment at left delamination tip yields(26)GJ1θ1′+P10IP,1θ1′+Mzz,1W1′x=x1=GJ2θ2′+P20IP,2θ2′+Mzz,2W2′x=x1+GJ3θ3′+P30IP,3θ3′+Mzz,3W3′x=x1.Using Galerkin weighted residual formulation, the integral forms of the governing differential equations (3) and (4), also representing the virtual flexural (W¯f) and torsional (W¯t) works are written as(27)W¯f=∑i=14∫0LiδWiEIiW″″i+PiWi″+Mzz,iθi″−ρAiω2Widx=0,(28)W¯t=∑i=14∫0LiδθiGJiθi″+PiIP,iθi″+Mzz,iWi″−ρIP,iAiω2θidx=0,where δWi and δθi are weighting functions corresponding to flexural and torsional displacements, respectively, for beam segment i. Then, performing integrations by parts twice on Equation (27) and once on Equation (28) leads to(29)W¯f=∑i=14∫0LiEIiWi″δWi″−PiWi′δWi′+Mzz,iθi′δWi′+ρAiω2WiδWidx+∑i=14EIiWi‴+PiWi′+Mzz,iθi′δWi0Li−∑i=14EIiWi″δWi′0Li︸Δ=0,W¯t=∑i=14∫0LGiJiθi′δθi′+PiIP,iθ′δθ′+Mzz,iWi′δθi′−ρiIP,iAiω2θiδθidx−∑i=14GiJiθi′+PiIP,iθi′ʲMzz,iWi′δθi0Li︸Δ∗=0.Considering expressions (20)–(22), it can be shown that for all the system’s global classic end conditions (at x=0 and x=L), the corresponding relationships in boundary terms (Δ) and (Δ∗) will vanish. For example, for clamped-free boundary conditions, one has zero displacements W1=W1′=θ1=0 and zero virtual displacements δW1=δW1′=δθ1=0 at the clamped end (x=0) and null resultant shear force, S4x, bending moment, M4x, and twisting moment, T4x at the free end (x=L). This leads to(30)EI1W1‴+P10W1′+Mzz,1θ1′δW1x=0−EI1W1″δW1′x=0=0,EI4W4‴+P40W4′+Mzz,4θ4′δW1x=x4−EI4W4″δW4′x=x4=0,GJ1θ1′+P10IP,1θ1′+Mzz,1W1′δθ1x=0=0,GJ4θ4′+P40IP,4θ4′+Mzz,4W4′δθ4x=x4=0.The remaining terms in (Δ), corresponding to delamination edges (i.e., x=x1 and x=x2), can be resolved by applying the continuity conditions (23) and (25), with the following as a result:(31)∑i=14EIiWi‴+Pi0Wi′+Mzz,iθi′δWi0Li−∑i=14EIiWi″δWi′0Li=δW2x=x1EI1W1‴+P10W1′+Mzz,1θ1′x=x1−EI2W2‴+P20W2′+Mzz,2θ2′x=x1−EI3W3‴+P30W3′+Mzz,3θ3′x=x1︸∗−δW2′x=x1EI1W1″x=x1−EI2W2″x=x1−EI3W3″x=x1+δW2x=x2EI4W4‴+P40W4′+Mzz,4θ4′x=x2−EI2W2‴+P20W2′+Mzz,2θ2′x=x2−EI3W3‴+P30W3′+Mzz,3θ3′x=x2︸∗∗−δW2′x=x2EI4W4″x=x2−EI2W2″x=x2−EI3W3″x=x2.The terms (∗) and (∗∗) go to zero directly because of shear force continuity conditions. The following remaining terms, however, do not vanish:(32)∑i=14EIiWi‴+Pi0Wi′+Mzz,iθi′δWi0Li=δW2′x=x1−δW2′x=x2ΛW4′x4=0−W1′x1=L1.Likewise, the remaining terms in (Δ∗) can be resolved by applying the continuity conditions (7, 11):(33)∑i=14GJiθi′+Pi0IP,iθi′+Mzz,iWi′δθi0Li=δθ2x=x1GJ1θ1′+P10IP,1θ1′+Mzz,1W1′x=x1−GJ2θ2′+P20IP,2θ2′+Mzz,2W2′x=x1−GJ3θ3′+P30IP,3θ3′+Mzz,3W3′x=x1︸∗∗∗+δθ2x=x2GJ4θ4′+P40IP,4θ4′+Mzz,4W4′x=x2−GJ2θ2′+P20IP,2θ2′+Mzz,2W2′x=x2−GJ3θ3′+P30IP,3θ3′+Mzz,3W3′x=x2︸∗∗∗∗.The terms (∗∗∗) and (∗∗∗∗) vanish as a result of twisting moment continuity conditions for left delamination tip (Equation (26)) and similar continuity condition for right tip. With the end conditions and continuity requirements satisfied and knowing that the expressions (27) and (28) must satisfy the virtual work principle, i.e., W¯f+W¯t=0, the system discretization leads to(34)ΛδW2′x=x1−δW2′x=x2W4′x4=0−W1′x1=L1+∑i=14∑m=1n∫xmxm+1EIiWi″δWi″−Pi0Wi′δWi′+Mzz,iθi′δWi′+ρAiω2WiδWidx+∑i=14∑m=1n∫xmxm+1GiJiθi′δθi′+Pi0IP,iθi′δθi′+Mzz,iWi′δθi′−ρiIP,iAiω2θiδθidx=0,where n stands for the number of elements in beam segment i. Following the traditional Euler–Bernoulli FE development, cubic Hermite shape functions for flexural displacement and linear shape functions for torsion displacement are used to relate the displacement within each element (m) to their corresponding nodal displacements:(35)Wix=NxWp,δWix=NxδWp,θix=LxWp,δθix=LxδWp,where, Nx and Lx are row vectors of cubic shape functions and linear shape functions, respectively. Wp, δWp, θp, and δθp are vectors containing the nodal real and virtual displacements for each element (m), defined as(36)Wp=W1W1′W2W2′T,δWp=δW1δW1′δW2δW2′T,(37)θp=θ1θ2T,δθp=δθ1δθ2T.Substituting the above approximate displacements in (34),(38)δWn′ΛN′x=x2−N′x=x1N′x=x2−N′x=x1Wn′+∑i=14∑m=1nδWn∫xmxm+1EIiN″N″−Pi0N′N′+Mzz,iL′N′+ρAiω2NNdxWn+∑i=14∑m=1nδWn∫xmxm+1GiJiL′L′+Pi0IP,iL′L′+Mzz,iN′L′−ρiIP,iAiω2LLdxWn=0,where Wn′ and δWn′ are real and virtual nodal displacement vectors for the two nodes at the two ends of the delamination, while Wn and δWn are the elemental nodal vectors of real and virtual displacements (36) and (37), rewritten as(39)Wn=W1W1′θ1W2W2′θ2T,δWn=δW1δW1′δθ1δW2δW2′δθ2T.Assembly of terms in (38) and enforcing the system’s global boundary conditions leads to the following linear eigenvalue problem:(40)δWnKωWn=0,withKω=KG−ω2M,where KG=KB+KBT+KTB+KT+KD and KB and KT are the global bending and torsion stiffness matrices, respectively. KBT and KTB are global bending-torsion and torsion-bending stiffness matrices, respectively, resulting from coupled terms in bending and torsion equations. KD is delamination stiffness matrix generated from the term outside the integral in Equation (38), and M and KG represent the overall (global) mass and stiffness matrices, respectively. Finally the system’s natural frequencies and corresponding modes are the system’s eigenvalues and eigenvectors evaluated from eigenproblem (40), i.e., for arbitrary δWn, Kω=0. The above procedure is achieved through a code developed in Matlab®.


### 2.2. Constrained Mode Model

In the simplified “constrained mode” model, the delaminated layers are “constrained” to have the same transverse deformations. The delaminated beam configuration is analyzed as assembly of three beam segments I, II, and III, as shown in Figure1, where 0 ≤ x ≤ x1 represents segment I, x1 ≤ x ≤ x2 constitutes segment II, and x2 ≤ x ≤ L is segment III.For beam segment I, the governing equations are(41)EI1W″″I+P0WI″+MzzθI″−ρA1ω2WI=0,GJ1θI″+P0IP,1θI″+MzzWI″−ρIP,1A1ω2θI=0.For segment II,(42)EI2+I3W″″II+P0WII″+MzzθII″−ρA2+A3ω2WII=0,GJ2+J3θII″+P0IP,2+IP,3θII″+MzzWII″−ρIP,2A2+IP,3A3ω2θII=0,and for segment III,(43)EI4W″″III+P0WIII″+MzzθIII″−ρA4ω2WIII=0,GJ4θIII″+P0IP,4θIII″+MzzWIII″−ρIP,4A4ω2θIII=0.The continuity conditions for lateral displacement, slope, shear force, bending moment, and torsion torque at the left edge of delamination,x=x1, are(44)WI,1=WII,1,WI,1=WIx=x1,WII,1=WIIx=x1,(45)W′I,1=W′II,1,W′I,1=WI′x=x1,W′II,1=WII′x=x1,(46)θI,1=θII,1,θI,1=θIx=x1,θII,1=θIIx=x1,(47)EI1WI‴+P0WI′+MzzθI′x=x1=EI2+I3WII‴+P0WII′+MzzθII′x=x1,(48)EI1WI″x=x1=EI2+I3WII″x=x1−P2tH32+P3tH22,(49)GJ1θI′+P0IP,1θI′+MzzWI′x=x1=GJ2+J3θII′+P20IP,2+IP,3θII′+MzzWII′x=x1.The continuity conditions at the right delamination tip are also written in a similar way to expressions (44)–(49). The general solution steps, namely, the assembly of element equations, application of boundary conditions, forming the system’s linear eigenproblem, and extraction of the natural frequencies and modes for the “constrained mode” are identical to those stated above for the “free mode” model.


## 2.1. Free Mode Delamination Model

In contrast to the “constraint mode” model, the “free mode” delamination scheme considers no constraint or interaction between the upper and lowest layers in the defective zone, i.e., segments 2 and 3. The “free mode” model predicts the behaviour of the beam more accurately but it also exhibits some physically impossible modes of vibration. To solve the four pairs of coupled equations as one system numerically, the global end conditions at both extremities of the beam, as well as continuity requirements at the delamination tips, are to be enforced. The continuity conditions for lateral displacement, slope, torsional displacement, and axial force at the left side of the defective region (x=x1) are expressed as(5)W11=W21=W31,W11=W1x=x1,W21=W2x=x1,W31=W3x=x1,(6)W11′=W21′=W31′,W11′=W1′x=x1,W21′=W2′x=x1,W31′=W3′x=x1,(7)θ11=θ21=θ31,θ11=θ1x=x1,θ21=θ2x=x1,θ31=θ3x=x1,(8)P1=P2+P3.Similarly, continuity conditions for the right delamination tip are written as(9)W42=W22=W32,W42=W4x=x2,W22=W2x=x2,W32=W3x=x2,W42′=W22′=W32′,W42′=W4′x=x2,W22′=W2′x=x2,W32′=W3′x=x2,θ42=θ22=θ32,θ42=θ4x=x2,θ22=θ2x=x2,θ32=θ3x=x2,P4=P2+P3.Additionally, the requirement for the delamination tip faces to remain planar after deformation, at the left delamination tip, results in(10)u21−u31=H12W11′,where ui is the axial displacement of beam section i. Combining expression (10) with the similar expression from the right delamination tip leads to(11)Δu3−Δu2=H12W42′−W11′,where(12)Δu3=u32−u31,u32=u3x=x2,u31=u3x=x1,Δu2=u22−u21,u22=u2x=x2,u21=u2x=x1,W42′=W4′x=x2,W11′=W1′x=x1.It was originally assumed by Mujumdar and Suryanarayan [5] and was later shown in an exact way by Szekrényes [21] that for small deformations of a uniform beam, the axial displacement will behave according to the following expression (13):(13)uixi=Li−uixi=0=∫0LiFixiEiAixidxi=FiLiEiAi,where xi represents the local position along the length and Fi is the axial force induced in the segment i. Compressive force during vibration, Px,t, consists of applied static axial force, P0, and perturbed axial force, Ptx,t. The latter, preventing interlaminar slip, is induced in the delaminated layers during vibration. Knowing that the change in axial displacement for each segment is caused by induced perturbed component of axial force rather than its static component, substituting Equation (13) into (11) yields(14)P2tL2EA2−P3tL3EA3=H12W4′x=x2−W1′x=x1.Taking the perturbation of Equation (8) for continuity of axial force on the left delamination tip, one obtains(15)P2t+P3t=P1t.SinceP1=P0 and P0 is constant, P1t=0 which means: P2t+P3t=0. Substituting P3t=−P2t into Equation (14) results in(16)P3t=Λ∗W4′x=x2−W1′x=x1,where the coefficient Λ∗ is defined as(17)Λ∗=H12L2EA2EA3EA3+EA2.At the left and right delamination tips, respectively, continuity of bending moments leads to(18)M11=M21+M31−P2tH32+P3tH22,(19)M42=M22+M32−P2tH32+P3tH22,where Mij=Mix=xj. The internal bending moment, Mx, shear force, Sx, and torque, Tx, in each of the individual beam segments are written as(20)Mix=−EIiWi″,(21)Six=EIiWi‴+Mzz,iθi′+Pi0Wi′,(22)Tix=GJiθi′+Pi0IP,iθi′+Mzz,iWi′.Using Equations (20) and (18), the following continuity of bending moments is obtained(23)EI1W1″x=x1=EI2W2″x=x1+EI3W3″x=x1+ΛW4′x4=0−W1′x1=L1,where the parameter Λ is defined as(24)Λ=H124L2EA2EA3EA3+EA2.Likewise, to satisfy the continuity of shear forces at the left delamination tip, one should have(25)EI1W1‴+P10W1′+Mzz,1θ1′x=x1=EI2W2‴+P20W2′+Mzz,2θ2′x=x1+EI3W3‴+P30W3′+Mzz,3θ3′x=x1,and finally, the continuity of twisting moment at left delamination tip yields(26)GJ1θ1′+P10IP,1θ1′+Mzz,1W1′x=x1=GJ2θ2′+P20IP,2θ2′+Mzz,2W2′x=x1+GJ3θ3′+P30IP,3θ3′+Mzz,3W3′x=x1.Using Galerkin weighted residual formulation, the integral forms of the governing differential equations (3) and (4), also representing the virtual flexural (W¯f) and torsional (W¯t) works are written as(27)W¯f=∑i=14∫0LiδWiEIiW″″i+PiWi″+Mzz,iθi″−ρAiω2Widx=0,(28)W¯t=∑i=14∫0LiδθiGJiθi″+PiIP,iθi″+Mzz,iWi″−ρIP,iAiω2θidx=0,where δWi and δθi are weighting functions corresponding to flexural and torsional displacements, respectively, for beam segment i. Then, performing integrations by parts twice on Equation (27) and once on Equation (28) leads to(29)W¯f=∑i=14∫0LiEIiWi″δWi″−PiWi′δWi′+Mzz,iθi′δWi′+ρAiω2WiδWidx+∑i=14EIiWi‴+PiWi′+Mzz,iθi′δWi0Li−∑i=14EIiWi″δWi′0Li︸Δ=0,W¯t=∑i=14∫0LGiJiθi′δθi′+PiIP,iθ′δθ′+Mzz,iWi′δθi′−ρiIP,iAiω2θiδθidx−∑i=14GiJiθi′+PiIP,iθi′ʲMzz,iWi′δθi0Li︸Δ∗=0.Considering expressions (20)–(22), it can be shown that for all the system’s global classic end conditions (at x=0 and x=L), the corresponding relationships in boundary terms (Δ) and (Δ∗) will vanish. For example, for clamped-free boundary conditions, one has zero displacements W1=W1′=θ1=0 and zero virtual displacements δW1=δW1′=δθ1=0 at the clamped end (x=0) and null resultant shear force, S4x, bending moment, M4x, and twisting moment, T4x at the free end (x=L). This leads to(30)EI1W1‴+P10W1′+Mzz,1θ1′δW1x=0−EI1W1″δW1′x=0=0,EI4W4‴+P40W4′+Mzz,4θ4′δW1x=x4−EI4W4″δW4′x=x4=0,GJ1θ1′+P10IP,1θ1′+Mzz,1W1′δθ1x=0=0,GJ4θ4′+P40IP,4θ4′+Mzz,4W4′δθ4x=x4=0.The remaining terms in (Δ), corresponding to delamination edges (i.e., x=x1 and x=x2), can be resolved by applying the continuity conditions (23) and (25), with the following as a result:(31)∑i=14EIiWi‴+Pi0Wi′+Mzz,iθi′δWi0Li−∑i=14EIiWi″δWi′0Li=δW2x=x1EI1W1‴+P10W1′+Mzz,1θ1′x=x1−EI2W2‴+P20W2′+Mzz,2θ2′x=x1−EI3W3‴+P30W3′+Mzz,3θ3′x=x1︸∗−δW2′x=x1EI1W1″x=x1−EI2W2″x=x1−EI3W3″x=x1+δW2x=x2EI4W4‴+P40W4′+Mzz,4θ4′x=x2−EI2W2‴+P20W2′+Mzz,2θ2′x=x2−EI3W3‴+P30W3′+Mzz,3θ3′x=x2︸∗∗−δW2′x=x2EI4W4″x=x2−EI2W2″x=x2−EI3W3″x=x2.The terms (∗) and (∗∗) go to zero directly because of shear force continuity conditions. The following remaining terms, however, do not vanish:(32)∑i=14EIiWi‴+Pi0Wi′+Mzz,iθi′δWi0Li=δW2′x=x1−δW2′x=x2ΛW4′x4=0−W1′x1=L1.Likewise, the remaining terms in (Δ∗) can be resolved by applying the continuity conditions (7, 11):(33)∑i=14GJiθi′+Pi0IP,iθi′+Mzz,iWi′δθi0Li=δθ2x=x1GJ1θ1′+P10IP,1θ1′+Mzz,1W1′x=x1−GJ2θ2′+P20IP,2θ2′+Mzz,2W2′x=x1−GJ3θ3′+P30IP,3θ3′+Mzz,3W3′x=x1︸∗∗∗+δθ2x=x2GJ4θ4′+P40IP,4θ4′+Mzz,4W4′x=x2−GJ2θ2′+P20IP,2θ2′+Mzz,2W2′x=x2−GJ3θ3′+P30IP,3θ3′+Mzz,3W3′x=x2︸∗∗∗∗.The terms (∗∗∗) and (∗∗∗∗) vanish as a result of twisting moment continuity conditions for left delamination tip (Equation (26)) and similar continuity condition for right tip. With the end conditions and continuity requirements satisfied and knowing that the expressions (27) and (28) must satisfy the virtual work principle, i.e., W¯f+W¯t=0, the system discretization leads to(34)ΛδW2′x=x1−δW2′x=x2W4′x4=0−W1′x1=L1+∑i=14∑m=1n∫xmxm+1EIiWi″δWi″−Pi0Wi′δWi′+Mzz,iθi′δWi′+ρAiω2WiδWidx+∑i=14∑m=1n∫xmxm+1GiJiθi′δθi′+Pi0IP,iθi′δθi′+Mzz,iWi′δθi′−ρiIP,iAiω2θiδθidx=0,where n stands for the number of elements in beam segment i. Following the traditional Euler–Bernoulli FE development, cubic Hermite shape functions for flexural displacement and linear shape functions for torsion displacement are used to relate the displacement within each element (m) to their corresponding nodal displacements:(35)Wix=NxWp,δWix=NxδWp,θix=LxWp,δθix=LxδWp,where, Nx and Lx are row vectors of cubic shape functions and linear shape functions, respectively. Wp, δWp, θp, and δθp are vectors containing the nodal real and virtual displacements for each element (m), defined as(36)Wp=W1W1′W2W2′T,δWp=δW1δW1′δW2δW2′T,(37)θp=θ1θ2T,δθp=δθ1δθ2T.Substituting the above approximate displacements in (34),(38)δWn′ΛN′x=x2−N′x=x1N′x=x2−N′x=x1Wn′+∑i=14∑m=1nδWn∫xmxm+1EIiN″N″−Pi0N′N′+Mzz,iL′N′+ρAiω2NNdxWn+∑i=14∑m=1nδWn∫xmxm+1GiJiL′L′+Pi0IP,iL′L′+Mzz,iN′L′−ρiIP,iAiω2LLdxWn=0,where Wn′ and δWn′ are real and virtual nodal displacement vectors for the two nodes at the two ends of the delamination, while Wn and δWn are the elemental nodal vectors of real and virtual displacements (36) and (37), rewritten as(39)Wn=W1W1′θ1W2W2′θ2T,δWn=δW1δW1′δθ1δW2δW2′δθ2T.Assembly of terms in (38) and enforcing the system’s global boundary conditions leads to the following linear eigenvalue problem:(40)δWnKωWn=0,withKω=KG−ω2M,where KG=KB+KBT+KTB+KT+KD and KB and KT are the global bending and torsion stiffness matrices, respectively. KBT and KTB are global bending-torsion and torsion-bending stiffness matrices, respectively, resulting from coupled terms in bending and torsion equations. KD is delamination stiffness matrix generated from the term outside the integral in Equation (38), and M and KG represent the overall (global) mass and stiffness matrices, respectively. Finally the system’s natural frequencies and corresponding modes are the system’s eigenvalues and eigenvectors evaluated from eigenproblem (40), i.e., for arbitrary δWn, Kω=0. The above procedure is achieved through a code developed in Matlab®.


## 2.2. Constrained Mode Model

In the simplified “constrained mode” model, the delaminated layers are “constrained” to have the same transverse deformations. The delaminated beam configuration is analyzed as assembly of three beam segments I, II, and III, as shown in Figure1, where 0 ≤ x ≤ x1 represents segment I, x1 ≤ x ≤ x2 constitutes segment II, and x2 ≤ x ≤ L is segment III.For beam segment I, the governing equations are(41)EI1W″″I+P0WI″+MzzθI″−ρA1ω2WI=0,GJ1θI″+P0IP,1θI″+MzzWI″−ρIP,1A1ω2θI=0.For segment II,(42)EI2+I3W″″II+P0WII″+MzzθII″−ρA2+A3ω2WII=0,GJ2+J3θII″+P0IP,2+IP,3θII″+MzzWII″−ρIP,2A2+IP,3A3ω2θII=0,and for segment III,(43)EI4W″″III+P0WIII″+MzzθIII″−ρA4ω2WIII=0,GJ4θIII″+P0IP,4θIII″+MzzWIII″−ρIP,4A4ω2θIII=0.The continuity conditions for lateral displacement, slope, shear force, bending moment, and torsion torque at the left edge of delamination,x=x1, are(44)WI,1=WII,1,WI,1=WIx=x1,WII,1=WIIx=x1,(45)W′I,1=W′II,1,W′I,1=WI′x=x1,W′II,1=WII′x=x1,(46)θI,1=θII,1,θI,1=θIx=x1,θII,1=θIIx=x1,(47)EI1WI‴+P0WI′+MzzθI′x=x1=EI2+I3WII‴+P0WII′+MzzθII′x=x1,(48)EI1WI″x=x1=EI2+I3WII″x=x1−P2tH32+P3tH22,(49)GJ1θI′+P0IP,1θI′+MzzWI′x=x1=GJ2+J3θII′+P20IP,2+IP,3θII′+MzzWII′x=x1.The continuity conditions at the right delamination tip are also written in a similar way to expressions (44)–(49). The general solution steps, namely, the assembly of element equations, application of boundary conditions, forming the system’s linear eigenproblem, and extraction of the natural frequencies and modes for the “constrained mode” are identical to those stated above for the “free mode” model.


## 3. Numerical Results

In this section, numerical tests performed to confirm the predictability and accuracy of the proposed FE models for delaminated beams, are presented. A Matlab® code was developed to create the model, assembly, and the application of the system boundary conditions, and the resulting linear eigenproblem was solved using the “eig” function.In the first two examples, the natural frequencies of defective steel and graphite/epoxy beam configurations, in the absence of axial load and end moment (i.e.,P0=Mzz=0), are investigated to validate and verify the applicability of the formulation against existing data. The systems are assumed to have a central split, about the midsection (L1 = L4) and of various lengths up to 60% of the span (0 ≤ L2/L ≤ 0.6). The split is occurring symmetrically along the midplane of the beam and surrounded by intact beam segments, as has also been presented and studied in [4, 18]. The defective FEM models are then created and used to evaluate the natural frequencies and mode shapes of various prestressed delaminated beam configurations. As the benchmarks for comparison and validation purposes, the results from references [4, 5, 7, 13, 18] were used.


### 3.1. Validation of Presented Formulation

#### 3.1.1. Steel Beam

Consider a delaminated steel beam with elastic modulus ofE = 200 GPa, mass density of ρ = 7800 kg/m3, length of 8 m, and a rectangular cross-sectional area of width of 0.4 m and depth of 0.2 m. The defect is assumed to be a midplane (H2 = 0.5H1) delamination of variable length of 0 ≤ L2/L ≤ 0.6 (refer to Figure 1). A sample of convergence study is presented in Figure 2, where the analytical results presented by Wang et al. [4] are used as the benchmark to evaluate the performance of the presented formulation. The un-prestressed defective beam, in this case, is assumed to be clamped- clamped and isotropic, with P0 = 0, Mzz = 0, H2/H = 0.3, and L2/L = 0.4. To validate the present FE solution method, the first two nondimensional natural frequencies, λ2, of a defective clamped-clamped beam with a through-the-width delamination occurring symmetrically about the midsection (L1 = L4) on the midplane (H2 = H3) for various delamination lengths are compared with the analytical results reported by Wang et al. [4], as well as the DSM [18] and FEM data [19] by Erdelyi and Hashemi, where(50)λi4=ωi2ρA1EI1L4,i=1,2,….Figure 2
Convergence analysis for clamped-clamped isotropic beam, withP0 = 0, Mzz = 0, H2/H = 0.3, and L2/L = 0.4.In the formulations presented by Wang et al. [4] and Erdelyi and Hashemi [18, 19], no prestress effect and torsional vibration are investigated. Therefore, for comparison purposes, only bending equation is solved and applied axial force and end moment are set to zero; P0 = 0 and Mzz = 0. Erdelyi and Hashemi [19] presented FEM-based defective beam models and the frequency data calculated from six- and ten-element meshes of 2-node FEM elements. Tables 1 and 2 summarize the system’s first two natural frequencies obtained using the presented method, in comparison with those reported in the literature [4, 6, 17–19]. Referring to Tables 1 and 2 of reference [18], they showed exact match between the 1st frequency values obtained from both meshes. The 2nd frequency values obtained from the 6- and 10-element meshes, however, were found to be slightly different. Therefore, in what follows, frequency data obtained from the latter mesh are used for comparison.Table 1
The first nondimensional natural frequency (λ12) of an isotropic beam with a midplane delamination, clamped-clamped boundary conditions.
Delamination length (L2/L)
1st natural frequency
Present 10 elements
Wang et al. [4]
FEM [19]
DSM [18]
Della and Shu [17]
Layerwise FEM [6]
0.00
22.39
22.39
22.39
—
22.37
22.36
0.10
22.37
22.37
22.37
22.37
22.37
22.36
0.20
22.36
22.35
22.36
22.36
22.36
22.35
0.30
22.24
22.23
22.24
22.24
22.24
22.23
0.40
21.84
21.83
21.84
21.83
21.83
21.82
0.50
20.89
20.88
20.89
20.89
20.89
20.88
0.60
19.29
19.29
19.29
19.30
19.30
19.28Table 2
The second nondimensional natural frequency (λ12) of an isotropic beam with a midplane delamination, clamped-clamped boundary conditions.
Delamination length (L2/L)
2nd natural frequency
Present 10 elements
Wang et al. [4]
FEM [19]
DSM [18]
Della and Shu [17]
Layerwise FEM [6]
0.00
61.61
61.67
61.67
—
61.67
61.61
0.10
60.76
60.76
60.80
60.76
60.76
60.74
0.20
55.99
55.97
55.99
55.99
55.97
55.95
0.30
49.03
49.00
49.00
49.03
49.00
48.97
0.40
43.90
43.87
43.89
43.90
43.87
43.86
0.50
41.55
41.45
41.52
41.55
41.45
41.50
0.60
41.03
40.93
41.03
41.04
40.93
41.01As can be observed from Tables1 and 2, the natural frequencies calculated from the present method are in excellent agreement with the analytical results reported by Wang et al. [4], Della and Shu [16, 17], FEM [19], and layerwise FEM [6] data, with discrepancies less than 0.1%.


#### 3.1.2. Laminated ([0/90]2s) T300/934 Graphite/Epoxy Cantilever Beam

For the second validation case, consider an unstressed, defective, cantilevered, laminated beam ([0/90]2), made with T300/934 graphite/epoxy, containing a single, variable-length central delamination, located at the midplane, as also reported by Saravanos and Hopkins [13]. The dimensions of the beam are the same as the one investigated in the experimental study by Shen and Grady [7], i.e., 10 × 0.5 × 0.0435 inches. Since the presented model in this paper is based on isotropic materials, the classic lamination theory (CLT) [27] is used to find the properties of equivalent single layer beam for each of the four segments in Figure 1. In addition, since the beam is unstressed, the axial load and end moment are set equal to zero. As can be observed from Figure 3, the FEM results obtained from the present theory show excellent agreement with both analytical and experimental data reported by Saravanos and Hopkins [13] and Shen and Grady [7]. It is to be noted that the analytical and experimental frequency data are extracted from Figures 4 and 5, reported in references [7, 13], respectively.Figure 3
Effect of central midplane delamination on the fundamental natural frequency of an unstressed cantilevered [0/90]2s T300/934 graphite/epoxy beam, ──, present model; - - -, analytical [13]; and ○, experiment [7].Figure 4
Fundamental natural frequencyλ2 of homogeneous beam with clamped-clamped boundary condition, Mzz/Mb = 0.4 and central delamination located in midplane H2/H = 0.5, versus normalized axial force in “free mode” model for different values of L2/L.Figure 5
Fundamental natural frequencies (λ2) with respect to normalized applied end moment Mzz/Mb for both constrained and free modes with H2/H = 0.5, P0/Pcr = 0.4, and different values of L2/L.


### 3.2. Vibration Analysis of Prestressed Delaminated Beams

To further investigate the validity and practical applicability of the proposed formulation and the effect of axial load on the frequency response of the defective beam structures, following illustrative examples of prestressed symmetrically delaminated beams (Figure1), characterized by various delamination lengths, applied axial force and end moments are investigated. In each case, the change of normalized fundamental frequency, λ2, versus the normalized end moment, Mzz/Mb, compressive axial load, P0/Pcr, for “free mode” is investigated, where(51)P0Pcr=P04π2EI1L4,where Pcr and Mb, respectively, stand for the critical buckling load and buckling moment of an intact (with no delamination) beam. Figure 4 shows the system’s fundamental natural frequency, λ2, of a homogeneous beam with clamped-clamped boundary condition, Mzz/Mb=0.4, and central delamination located in midplane H2/H=0.5, versus normalized axial force in “free mode” model for different delamination sizes, L2/L. As can be observed, when the applied axial load, P0, approaches the system’s buckling load, the natural frequency, λ2, turns into zero, i.e., the structure has buckled. Moreover, as expected, the larger the delamination size, the smaller the Pcr and the fundamental frequency. This reconfirms the fact that structural stiffness decreases when the delamination size increased.Figure6 demonstrates the change in the first and second natural frequencies with respect to normalized buckling load P0/Pcr for Mzz/Mb = 0.4 and H2/H = 0.3. Results are shown for different delamination lengths of L2/L = 0.2, 0.3, 0.4, and 0.5. The figure reveals similar trends to those observed in Figure 2. Fundamental frequency results of “constrained mode” are compared to “free mode” data in Figure 7. As can be observed, the difference between predicted frequencies from these two modes increases as the applied axial load increases. But for both natural modes, the larger the compressive axial load, the smaller the frequency. It can be observed that as P0 increases, the natural frequency λ2 first decreases slowly, but by getting close to buckling load it decreases drastically and eventually goes to zero. As can also be seen from Figure 7, the “constrained mode” comparing to “free mode,” overpredicts the natural frequencies especially for compressive load values close to buckling load. Additionally, from the same trends observed in Figures 2–4, it can be inferred that compressive axial load causes a reduction in the structural stiffness leading, in turn, to lower natural frequencies. Similarly, increasing the delamination length results in lower stiffness and natural frequencies.Figure 6
Changes of the nondimensional first and second natural frequencies (λ2) with respect to normalized buckling load P0/Pcr for H2/H = 0.3, Mzz/Mb = 0.4, and different values of L2/L.Figure 7
Trend of change in fundamental natural frequencies (λ2) with respect to normalized buckling load P0/Pcr for both constrained and free modes with H2/H = 0.5, Mzz/Mb = 0.4, and different values of L2/L.Figure5 presents the normalized natural frequency, λ2, versus the normalized end moment, Mzz/Mb, for both constrained and free modes with H2/H = 0.5, P0/Pcr = 0.4, and different values of L2/L = 0.3, 0.4, and 0.5. Referring to Figure 5, the same trends as in the previous graphs (for axial force) is observed for variation of natural frequency versus buckling moment. The difference between predicted frequencies from these two modes increases as the applied end moment increases, and the larger the end moment, the smaller the frequency. It can also be observed that as Mzz increases, the natural frequency λ2 first decreases slowly, but by getting close to buckling load it decreases drastically and eventually goes to zero. In addition, the “constrained mode” comparing to “free mode,” overpredicts the natural frequencies. Once again, it can also be seen that as the end moment increases, the stiffness decreases, leading to lower natural frequencies.Figure8 depicts the mode shape associated with fundamental natural frequency of an intact beam in comparison with the first opening mode of vibration for a defective beam, characterized by a central delamination on the midplane, obtained using “free mode” model. This mode, however, cannot be detected using the constrained model. As also reported in the literature [19, 22], not all the mode shapes obtained from “free mode” model are admissible since unlike “constrained mode” no constraints are imposed between the two delaminated layers, which results in some physically impossible mode shapes. It is also worth noting that the experimentally observed delamination openings have been reported to be much fewer than those suggested by the free model, which would rather justify the use of the constrained model for further analysis [21].Figure 8
The first opening mode shape for a delaminated beam with a central delamination on the midplane using “free mode” model.Finally, it has been recently shown and experimentally verified by Szekrényes [22] that the opening is amplitude dependent. In other words, the delamination opening takes place only at certain critical amplitudes. The mode shapes are asymmetric and can be approximated by the superposition of the global shape of the entire beam and the local buckling eigenshape of the delaminated part based on a dynamic stability analysis [21]. Further investigation on this topic, however, is beyond the scope of this paper.


## 3.1. Validation of Presented Formulation

### 3.1.1. Steel Beam

Consider a delaminated steel beam with elastic modulus ofE = 200 GPa, mass density of ρ = 7800 kg/m3, length of 8 m, and a rectangular cross-sectional area of width of 0.4 m and depth of 0.2 m. The defect is assumed to be a midplane (H2 = 0.5H1) delamination of variable length of 0 ≤ L2/L ≤ 0.6 (refer to Figure 1). A sample of convergence study is presented in Figure 2, where the analytical results presented by Wang et al. [4] are used as the benchmark to evaluate the performance of the presented formulation. The un-prestressed defective beam, in this case, is assumed to be clamped- clamped and isotropic, with P0 = 0, Mzz = 0, H2/H = 0.3, and L2/L = 0.4. To validate the present FE solution method, the first two nondimensional natural frequencies, λ2, of a defective clamped-clamped beam with a through-the-width delamination occurring symmetrically about the midsection (L1 = L4) on the midplane (H2 = H3) for various delamination lengths are compared with the analytical results reported by Wang et al. [4], as well as the DSM [18] and FEM data [19] by Erdelyi and Hashemi, where(50)λi4=ωi2ρA1EI1L4,i=1,2,….Figure 2
Convergence analysis for clamped-clamped isotropic beam, withP0 = 0, Mzz = 0, H2/H = 0.3, and L2/L = 0.4.In the formulations presented by Wang et al. [4] and Erdelyi and Hashemi [18, 19], no prestress effect and torsional vibration are investigated. Therefore, for comparison purposes, only bending equation is solved and applied axial force and end moment are set to zero; P0 = 0 and Mzz = 0. Erdelyi and Hashemi [19] presented FEM-based defective beam models and the frequency data calculated from six- and ten-element meshes of 2-node FEM elements. Tables 1 and 2 summarize the system’s first two natural frequencies obtained using the presented method, in comparison with those reported in the literature [4, 6, 17–19]. Referring to Tables 1 and 2 of reference [18], they showed exact match between the 1st frequency values obtained from both meshes. The 2nd frequency values obtained from the 6- and 10-element meshes, however, were found to be slightly different. Therefore, in what follows, frequency data obtained from the latter mesh are used for comparison.Table 1
The first nondimensional natural frequency (λ12) of an isotropic beam with a midplane delamination, clamped-clamped boundary conditions.
Delamination length (L2/L)
1st natural frequency
Present 10 elements
Wang et al. [4]
FEM [19]
DSM [18]
Della and Shu [17]
Layerwise FEM [6]
0.00
22.39
22.39
22.39
—
22.37
22.36
0.10
22.37
22.37
22.37
22.37
22.37
22.36
0.20
22.36
22.35
22.36
22.36
22.36
22.35
0.30
22.24
22.23
22.24
22.24
22.24
22.23
0.40
21.84
21.83
21.84
21.83
21.83
21.82
0.50
20.89
20.88
20.89
20.89
20.89
20.88
0.60
19.29
19.29
19.29
19.30
19.30
19.28Table 2
The second nondimensional natural frequency (λ12) of an isotropic beam with a midplane delamination, clamped-clamped boundary conditions.
Delamination length (L2/L)
2nd natural frequency
Present 10 elements
Wang et al. [4]
FEM [19]
DSM [18]
Della and Shu [17]
Layerwise FEM [6]
0.00
61.61
61.67
61.67
—
61.67
61.61
0.10
60.76
60.76
60.80
60.76
60.76
60.74
0.20
55.99
55.97
55.99
55.99
55.97
55.95
0.30
49.03
49.00
49.00
49.03
49.00
48.97
0.40
43.90
43.87
43.89
43.90
43.87
43.86
0.50
41.55
41.45
41.52
41.55
41.45
41.50
0.60
41.03
40.93
41.03
41.04
40.93
41.01As can be observed from Tables1 and 2, the natural frequencies calculated from the present method are in excellent agreement with the analytical results reported by Wang et al. [4], Della and Shu [16, 17], FEM [19], and layerwise FEM [6] data, with discrepancies less than 0.1%.


### 3.1.2. Laminated ([0/90]2s) T300/934 Graphite/Epoxy Cantilever Beam

For the second validation case, consider an unstressed, defective, cantilevered, laminated beam ([0/90]2), made with T300/934 graphite/epoxy, containing a single, variable-length central delamination, located at the midplane, as also reported by Saravanos and Hopkins [13]. The dimensions of the beam are the same as the one investigated in the experimental study by Shen and Grady [7], i.e., 10 × 0.5 × 0.0435 inches. Since the presented model in this paper is based on isotropic materials, the classic lamination theory (CLT) [27] is used to find the properties of equivalent single layer beam for each of the four segments in Figure 1. In addition, since the beam is unstressed, the axial load and end moment are set equal to zero. As can be observed from Figure 3, the FEM results obtained from the present theory show excellent agreement with both analytical and experimental data reported by Saravanos and Hopkins [13] and Shen and Grady [7]. It is to be noted that the analytical and experimental frequency data are extracted from Figures 4 and 5, reported in references [7, 13], respectively.Figure 3
Effect of central midplane delamination on the fundamental natural frequency of an unstressed cantilevered [0/90]2s T300/934 graphite/epoxy beam, ──, present model; - - -, analytical [13]; and ○, experiment [7].Figure 4
Fundamental natural frequencyλ2 of homogeneous beam with clamped-clamped boundary condition, Mzz/Mb = 0.4 and central delamination located in midplane H2/H = 0.5, versus normalized axial force in “free mode” model for different values of L2/L.Figure 5
Fundamental natural frequencies (λ2) with respect to normalized applied end moment Mzz/Mb for both constrained and free modes with H2/H = 0.5, P0/Pcr = 0.4, and different values of L2/L.


## 3.1.1. Steel Beam

Consider a delaminated steel beam with elastic modulus ofE = 200 GPa, mass density of ρ = 7800 kg/m3, length of 8 m, and a rectangular cross-sectional area of width of 0.4 m and depth of 0.2 m. The defect is assumed to be a midplane (H2 = 0.5H1) delamination of variable length of 0 ≤ L2/L ≤ 0.6 (refer to Figure 1). A sample of convergence study is presented in Figure 2, where the analytical results presented by Wang et al. [4] are used as the benchmark to evaluate the performance of the presented formulation. The un-prestressed defective beam, in this case, is assumed to be clamped- clamped and isotropic, with P0 = 0, Mzz = 0, H2/H = 0.3, and L2/L = 0.4. To validate the present FE solution method, the first two nondimensional natural frequencies, λ2, of a defective clamped-clamped beam with a through-the-width delamination occurring symmetrically about the midsection (L1 = L4) on the midplane (H2 = H3) for various delamination lengths are compared with the analytical results reported by Wang et al. [4], as well as the DSM [18] and FEM data [19] by Erdelyi and Hashemi, where(50)λi4=ωi2ρA1EI1L4,i=1,2,….Figure 2
Convergence analysis for clamped-clamped isotropic beam, withP0 = 0, Mzz = 0, H2/H = 0.3, and L2/L = 0.4.In the formulations presented by Wang et al. [4] and Erdelyi and Hashemi [18, 19], no prestress effect and torsional vibration are investigated. Therefore, for comparison purposes, only bending equation is solved and applied axial force and end moment are set to zero; P0 = 0 and Mzz = 0. Erdelyi and Hashemi [19] presented FEM-based defective beam models and the frequency data calculated from six- and ten-element meshes of 2-node FEM elements. Tables 1 and 2 summarize the system’s first two natural frequencies obtained using the presented method, in comparison with those reported in the literature [4, 6, 17–19]. Referring to Tables 1 and 2 of reference [18], they showed exact match between the 1st frequency values obtained from both meshes. The 2nd frequency values obtained from the 6- and 10-element meshes, however, were found to be slightly different. Therefore, in what follows, frequency data obtained from the latter mesh are used for comparison.Table 1
The first nondimensional natural frequency (λ12) of an isotropic beam with a midplane delamination, clamped-clamped boundary conditions.
Delamination length (L2/L)
1st natural frequency
Present 10 elements
Wang et al. [4]
FEM [19]
DSM [18]
Della and Shu [17]
Layerwise FEM [6]
0.00
22.39
22.39
22.39
—
22.37
22.36
0.10
22.37
22.37
22.37
22.37
22.37
22.36
0.20
22.36
22.35
22.36
22.36
22.36
22.35
0.30
22.24
22.23
22.24
22.24
22.24
22.23
0.40
21.84
21.83
21.84
21.83
21.83
21.82
0.50
20.89
20.88
20.89
20.89
20.89
20.88
0.60
19.29
19.29
19.29
19.30
19.30
19.28Table 2
The second nondimensional natural frequency (λ12) of an isotropic beam with a midplane delamination, clamped-clamped boundary conditions.
Delamination length (L2/L)
2nd natural frequency
Present 10 elements
Wang et al. [4]
FEM [19]
DSM [18]
Della and Shu [17]
Layerwise FEM [6]
0.00
61.61
61.67
61.67
—
61.67
61.61
0.10
60.76
60.76
60.80
60.76
60.76
60.74
0.20
55.99
55.97
55.99
55.99
55.97
55.95
0.30
49.03
49.00
49.00
49.03
49.00
48.97
0.40
43.90
43.87
43.89
43.90
43.87
43.86
0.50
41.55
41.45
41.52
41.55
41.45
41.50
0.60
41.03
40.93
41.03
41.04
40.93
41.01As can be observed from Tables1 and 2, the natural frequencies calculated from the present method are in excellent agreement with the analytical results reported by Wang et al. [4], Della and Shu [16, 17], FEM [19], and layerwise FEM [6] data, with discrepancies less than 0.1%.


## 3.1.2. Laminated ([0/90]2s) T300/934 Graphite/Epoxy Cantilever Beam

For the second validation case, consider an unstressed, defective, cantilevered, laminated beam ([0/90]2), made with T300/934 graphite/epoxy, containing a single, variable-length central delamination, located at the midplane, as also reported by Saravanos and Hopkins [13]. The dimensions of the beam are the same as the one investigated in the experimental study by Shen and Grady [7], i.e., 10 × 0.5 × 0.0435 inches. Since the presented model in this paper is based on isotropic materials, the classic lamination theory (CLT) [27] is used to find the properties of equivalent single layer beam for each of the four segments in Figure 1. In addition, since the beam is unstressed, the axial load and end moment are set equal to zero. As can be observed from Figure 3, the FEM results obtained from the present theory show excellent agreement with both analytical and experimental data reported by Saravanos and Hopkins [13] and Shen and Grady [7]. It is to be noted that the analytical and experimental frequency data are extracted from Figures 4 and 5, reported in references [7, 13], respectively.Figure 3
Effect of central midplane delamination on the fundamental natural frequency of an unstressed cantilevered [0/90]2s T300/934 graphite/epoxy beam, ──, present model; - - -, analytical [13]; and ○, experiment [7].Figure 4
Fundamental natural frequencyλ2 of homogeneous beam with clamped-clamped boundary condition, Mzz/Mb = 0.4 and central delamination located in midplane H2/H = 0.5, versus normalized axial force in “free mode” model for different values of L2/L.Figure 5
Fundamental natural frequencies (λ2) with respect to normalized applied end moment Mzz/Mb for both constrained and free modes with H2/H = 0.5, P0/Pcr = 0.4, and different values of L2/L.


## 3.2. Vibration Analysis of Prestressed Delaminated Beams

To further investigate the validity and practical applicability of the proposed formulation and the effect of axial load on the frequency response of the defective beam structures, following illustrative examples of prestressed symmetrically delaminated beams (Figure1), characterized by various delamination lengths, applied axial force and end moments are investigated. In each case, the change of normalized fundamental frequency, λ2, versus the normalized end moment, Mzz/Mb, compressive axial load, P0/Pcr, for “free mode” is investigated, where(51)P0Pcr=P04π2EI1L4,where Pcr and Mb, respectively, stand for the critical buckling load and buckling moment of an intact (with no delamination) beam. Figure 4 shows the system’s fundamental natural frequency, λ2, of a homogeneous beam with clamped-clamped boundary condition, Mzz/Mb=0.4, and central delamination located in midplane H2/H=0.5, versus normalized axial force in “free mode” model for different delamination sizes, L2/L. As can be observed, when the applied axial load, P0, approaches the system’s buckling load, the natural frequency, λ2, turns into zero, i.e., the structure has buckled. Moreover, as expected, the larger the delamination size, the smaller the Pcr and the fundamental frequency. This reconfirms the fact that structural stiffness decreases when the delamination size increased.Figure6 demonstrates the change in the first and second natural frequencies with respect to normalized buckling load P0/Pcr for Mzz/Mb = 0.4 and H2/H = 0.3. Results are shown for different delamination lengths of L2/L = 0.2, 0.3, 0.4, and 0.5. The figure reveals similar trends to those observed in Figure 2. Fundamental frequency results of “constrained mode” are compared to “free mode” data in Figure 7. As can be observed, the difference between predicted frequencies from these two modes increases as the applied axial load increases. But for both natural modes, the larger the compressive axial load, the smaller the frequency. It can be observed that as P0 increases, the natural frequency λ2 first decreases slowly, but by getting close to buckling load it decreases drastically and eventually goes to zero. As can also be seen from Figure 7, the “constrained mode” comparing to “free mode,” overpredicts the natural frequencies especially for compressive load values close to buckling load. Additionally, from the same trends observed in Figures 2–4, it can be inferred that compressive axial load causes a reduction in the structural stiffness leading, in turn, to lower natural frequencies. Similarly, increasing the delamination length results in lower stiffness and natural frequencies.Figure 6
Changes of the nondimensional first and second natural frequencies (λ2) with respect to normalized buckling load P0/Pcr for H2/H = 0.3, Mzz/Mb = 0.4, and different values of L2/L.Figure 7
Trend of change in fundamental natural frequencies (λ2) with respect to normalized buckling load P0/Pcr for both constrained and free modes with H2/H = 0.5, Mzz/Mb = 0.4, and different values of L2/L.Figure5 presents the normalized natural frequency, λ2, versus the normalized end moment, Mzz/Mb, for both constrained and free modes with H2/H = 0.5, P0/Pcr = 0.4, and different values of L2/L = 0.3, 0.4, and 0.5. Referring to Figure 5, the same trends as in the previous graphs (for axial force) is observed for variation of natural frequency versus buckling moment. The difference between predicted frequencies from these two modes increases as the applied end moment increases, and the larger the end moment, the smaller the frequency. It can also be observed that as Mzz increases, the natural frequency λ2 first decreases slowly, but by getting close to buckling load it decreases drastically and eventually goes to zero. In addition, the “constrained mode” comparing to “free mode,” overpredicts the natural frequencies. Once again, it can also be seen that as the end moment increases, the stiffness decreases, leading to lower natural frequencies.Figure8 depicts the mode shape associated with fundamental natural frequency of an intact beam in comparison with the first opening mode of vibration for a defective beam, characterized by a central delamination on the midplane, obtained using “free mode” model. This mode, however, cannot be detected using the constrained model. As also reported in the literature [19, 22], not all the mode shapes obtained from “free mode” model are admissible since unlike “constrained mode” no constraints are imposed between the two delaminated layers, which results in some physically impossible mode shapes. It is also worth noting that the experimentally observed delamination openings have been reported to be much fewer than those suggested by the free model, which would rather justify the use of the constrained model for further analysis [21].Figure 8
The first opening mode shape for a delaminated beam with a central delamination on the midplane using “free mode” model.Finally, it has been recently shown and experimentally verified by Szekrényes [22] that the opening is amplitude dependent. In other words, the delamination opening takes place only at certain critical amplitudes. The mode shapes are asymmetric and can be approximated by the superposition of the global shape of the entire beam and the local buckling eigenshape of the delaminated part based on a dynamic stability analysis [21]. Further investigation on this topic, however, is beyond the scope of this paper.


## 4. Discussion and Concluding Remarks

A systematic FEM-based formulation and numerical solutions for the free vibration modeling and analysis of delaminated beams under axial compressive load and end moment were presented. The delaminated beam configurations are modeled and analyzed as the assemblage of four Euler–Bernoulli beams attached at the defect edges. Both the “free mode” and “constrained mode” models were created to carry out buckling and vibration analyses, and the results were validated against the analytical data and other numerical results available in the open literature. Parametric studies were also performed to investigate the effects of end moment and axial load on the natural frequencies and mode shapes of delaminated beams. The results exhibited a monotonic relation between the system’s natural frequency and both loads (i.e., axial compressive force and end moment). Both compressive axial load and end moment reduce the structural stiffness and result in lower predictions for natural frequencies. Similarly, increasing the length of delamination results in lower stiffness and natural frequencies. Based on the “free mode” model, the first and second mode buckling loads of the delaminated beams were also studied. In general, the “constrained mode,” when compared to “free mode,” was shown to overpredict the natural frequencies and specifically when compressive load approached to values close to buckling load. Finally, the difference between predicted frequencies from “free mode” and “constrained mode” delamination models increases as the applied axial load or end moment increases.



---

*Source: 1348970-2018-09-25.xml*