Abstract
A dynamic stiffness element for flexural vibration analysis of delaminated multilayer beams is developed and subsequently used to investigate the natural frequencies and modes of two-layer beam configurations. Using the Euler-Bernoulli bending beam theory, the governing differential equations are exploited and representative, frequency-dependent, field variables are chosen based on the closed form solution to these equations. The boundary conditions are then imposed to formulate the dynamic stiffness matrix (DSM), which relates harmonically varying loads to harmonically varying displacements at the beam ends. The bending vibration of an illustrative example problem, characterized by delamination zone of variable length, is investigated. Two computer codes, based on the conventional Finite Element Method (FEM) and the analytical solutions reported in the literature, are also developed and used for comparison. The intact and defective beam natural frequencies and modes obtained from the proposed DSM method are presented along with the FEM and analytical results and those available in the literature.
1. Introduction
Layered structures have seen greatly increased use in civil, shipbuilding, mechanical, and aerospace structural applications in recent decades, 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. The replacement of traditionally metallic structural components with laminated composites has resulted in new and unique design challenges. Metallic structures exhibit mainly isotropic material properties and failure modes. By contrast, composite materials are anisotropic, which can result in more complex failure modes. Delamination is a common failure mode in layered structures. It may arise from loss of adhesion between two layers of the structure, from interlaminar stresses arising from geometric or material discontinuities, or from mechanical loadings. The presence of delamination may significantly reduce the stiffness and strength of the structures. A reduction in the stiffness will affect the vibration characteristics of the structures, such as the natural frequencies and mode shapes. Changes in the natural frequency, as a direct result of the reduction of stiffness, may lead to resonance if the reduced frequency is close to an excitation frequency.
The dynamic modeling of flexible delaminated multi-layer beams has been a topic of interest for many researchers. With the increase in use of laminated composite structures, the requirement for accurate delamination models has also grown. The earliest delamination models, formulated in the 1980s [1], dealt with the vibration of two-layer sandwich beams, where layers were governed by the Euler-Bernoulli bending beam theory. The upper and lower intact portions of the delaminated segment were assumed to vibrate freely—independent of each other; as a result, this model is known as “free mode” delamination. It was later discovered that the free mode underpredicted natural frequencies for off-midplane delaminations due to unrestricted penetration of the beams into each other. This was accounted for in 1988 [2] by constraining the transverse displacements of the top and bottom beams to be equal. The resulting model, known as the “constrained mode” delamination model, predicts vibration behaviour much more accurately for off-midplane delamination. However, in modeling terms, the constrained mode is simply a limiting case of the free mode delamination model, from which it can be derived. Thus, in the present study, the free mode delamination model will be investigated, and the constrained mode delamination model can be derived in the same manner as presented. It is also worth noting that, in general, a laminated composite beam may have orthotropic layerwise material properties, resulting in displacement coupling behaviour. The model used in this study assumes an isotropic material and is not readily applicable to fibre-reinforced laminated composite beams, as there would, in general, be a torsional and/or extensional response coupled with flexural vibration [3–5]. While the model used in this study assumes an isotropic materials (i.e., no coupled response), the proposed technique could be extended—keeping the delamination continuity conditions—to include more realistic composite response.
The accuracy of dynamic analysis, and forced response calculation, of a flexible structure depends greatly on the reliability of the modal analysis method used and the resulting natural frequencies and modes. There are various numerical, semianalytical, and analytical methods to extract the natural frequencies of a system. The conventional Finite Element Method (FEM) has a long well-established history and is the most commonly used method for modal analysis. The FEM is a general systematic approach to formulate the constant element mass and stiffness matrices for a given system and is easily adaptable to complex systems containing variations in geometry or loading. Nonuniform geometry, for example, is often modeled as a stepped, piecewise-uniform configuration. The conventional FEM formulation, based on polynomial shape functions, leads to constant mass and stiffness matrices and results in a linear eigenvalue problem from which the natural frequencies and modes of the system can be readily extracted. Lee [6], amongst others, investigated the free vibration analysis of laminated beams with delamination using a conventional FEM. Based on layerwise theory, equations of motion were derived from Hamilton’s principle, a Finite Element Method (FEM) was developed to formulate the problem, and the effects of location, size, and number of delaminations on vibration frequencies of delaminated beams were investigated [6]. During the last decade, sandwich and composite elements have been made available in some commercial software packages and are used to analyze the vibration of composite structures.
Alternatively, one can use semianalytical formulations, such as the so-called Dynamic Finite Element (DFE) method [7, 8], to carry out structural modal analysis. The hybrid DFE formulation results in a more accurate prediction method than traditional and FEM modeling techniques, allowing for a reduced mesh size. The main principle of the DFE is the Weighted Residual Integral formulation, providing a general systematic modeling procedure. The word Dynamic in DFE acronym refers to the frequency-dependent shape functions used to express the displacements, which in turn lead to the stiffness matrix of the system. The DFE technique follows the same typical procedure as the FEM by formulating the element equations discretized to a local domain, where element stiffness matrices are constructed and then assembled into a single global matrix. The application of the DFE to free vibration analysis of a delaminated 2-layer beam has been reported in an earlier work by the authors [9].
Analytical methods, namely the dynamic stiffness matrix (DSM), have also been used for the vibrational analysis of isotropic [10, 11], sandwich [12–14], and composite structural elements [15] and beamstructures [15, 16]. The DSM approach makes use of the general closed-form solution to the governing differential equations of motion of the system to formulate a frequency-dependent stiffness matrix. The DSM describes the free vibration of the system and exhibits both inertia and stiffness properties of the syetem. Based on this exact member theory, the DSM produces exact results for simple structural elements, such as uniform beams, using only one element [10, 11]. Banerjee and his coworkers [10–16] have developed a number of DSM formulations for various beam configurations, where the root-finding technique proposed by Wittrick-Williams (W-W) [17] was exploited to determine the eigenvalues of the system. The DSM has also been used by Wang et al. [18] to simulate a cracked beam. Wang [19] also investigated the effects of a through-thickness crack on the free vibration modes, aeroelastic flutter, and divergence of a composite wing. Borneman et al. [20] presented explicit expressions of a DSM for the coupled composite beams, exhibiting both material and geometric couplings. These expressions were consequently used to develop a cracked DSM formulation, and the free vibration of doubly coupled cracked composite beams was investigated. Given these considerations, the DSM method for a single beam can be modified to accurately model delaminated multilayer beams. A DSM-based preliminary analysis of a two-layer split beam has laso been presented in an earlier work by the authors [21].
The aim of this paper is to present a DSM formulation for the free vibration analysis of a delaminated two-layer beam, using the free mode delamination model. The delamination is represented by two intact beam segments; one for each of the top and bottom sections of the delamination. The delaminated region is bounded on either side by intact, full-height beams. The beams transverse displacements are governed by the Euler-Bernoulli slender beam bending theory. Shear deformation and rotary inertia, commonly associated with Timoshenko beam theory, are neglected. For harmonic oscillation, the governing equations are developed and used as the basis for the DSM development. Continuities of forces, moments, displacements, and slopes at the delamination tips are enforced, leading to the DSM of the system. Assembly of element DSM matrices and the application of boundary conditions results in the nonlinear eigenvalue problem of the defective system. In addition, two computer codes, based on the conventional Finite Element Method (FEM) and the analytical solutions reported in the literature [7, 8], taking into account the same continuity conditions, are also developed and used as a benchmark for comparison. The FEM model exploits cubic Hermite interpolation functions of approximation to express the flexural displacement functions, that is, field variables and weighting functions [22, 23]. Both DSM and FEM models are used to compute the natural frequencies of an illustrative example problem characterized by delamination zone of variable length. The frequency values are then compared with those from the literature. Certain modal characteristics of the system are also discussed.
2. Mathematical Model
Figure 1 shows the general coordinate system and notation for a delaminated beam, with total length , intact beam segment lengths and , delamination length , and total height . This model incorporates a general delamination, which can include laminated composites or bilayered isotropic materials, with different material and geometric properties above and below the delamination plane. Thus, the top layer has thickness , Young’s modulus , density , cross-sectional area , and second moment of area . The bottom layer has corresponding properties, with subscript 3. The delamination tips occur at stations and , and torsion, shear deformation, axial (warping effects and axial deformation), and out-of-plane delamination are ignored. Following this notation, the general equation of motion for the th Euler-Bernoulli beam in free vibration is written as [8, 9]: For harmonic oscillations, the transverse displacements can be described in the frequency domain by using the transformation where is the circular frequency of excitation of the system, is the amplitude of the displacement , and subscript “” represents the beam segment number. By backsubstituting (2) into (1), the equations of motion reduce to The general solution to the 4th-order, homogeneous differential equation (3) can be written in the following form: which represents the bending displacement of beam segment “,” is the beam segment length, and stands for nondimensional frequency of oscillation, defined as: Coefficients , , , and are evaluated to satisfy the displacement continuity requirements of the beam segments and the system boundary conditions. As also observed and reported by several researchers [8, 9], the inclusion of delamination into the beam model results in a coupling between axial and transverse motion of the delaminated beam segments. This is primarily due to the continuity requirements imposed on the delaminated beam endpoints at the delamination tips. Since the delamination tip cross-sections are assumed to remain planar after deformation, the ends of the top and bottom beams must have the same relative axial location after deformation, preventing interlaminar slip. Since the midplanes (assumed to be the neutral axes of the beam segments) in the delaminated segments are located at a distance from the midplanes of the intact segments, they will not have the same axial deformation unless some internal axial force is imposed. This imposed axial force is fully derived and discussed in [2], however, the final result will be briefly presented here for completeness.
Consider a delamination tip after deformation. According to the numbering scheme in Figure 1, and since no external axial load is applied, the top and bottom beam segments must have equal and opposite internal axial forces, that is, , applied to prevent interlaminar slip (Figure 2), where is the slope of the th beam segment, where “prime” represents the differentiation with respect to the beam longitudinal axis, , and the parameter is defined as which can be further simplified if the cross-sectional shape is known. With explicit expressions (6) and (7) for the internal axial force, continuity conditions for bending moment can be derived as follows:
at stations , continuity of bending moments leads to Using expression (8) and the previous conditions (6) and (7) for internal axial force, and noting that from beam theory, bending moments and shear forces in beam segment “” are related to displacements, , through , and , respectively, one can write for , where
Likewise, using (9) for , a similar relationship can be derived. Two boundary conditions at the intact beam ends, continuity of displacements and slopes, at the delamination tips results in 12 equations. Along with an additional four equations resulting from the continuity of bending moments and shear forces at the delamination tips, the total 16 equations can be used to solve for the 16 unknowns, for each beam, , as appearing in (4). This solution method, based on finding the coefficient matrix of the system, herein refered to as the “Coefficients Method (CM),” has been used to predict the vibration behavior of different systems of varying complexity (see, e.g., [7, 8]). However, it remains a relatively problem-specific solution technique. Thus, in what follows, this technique is reformulated into an equivalent, yet more readily and more conveniently applicable DSM formulation.
Through continuity conditions, a coupling relationship can be found within the delamination region to reduce the total number of unknowns from eight (, , for the top and bottom beams within the delaminated region) to four. Of particular interest are the continuity conditions for displacement and slope at the delamination tips, from which a coupling between the coefficients for the top beam and the bottom one can be derived. Stemming from the requirement that the displacement and slope of each beam, at the delamination tips, must be equal, the transverse displacements of beam segments 2 and 3 can be linked through the following relationship: Similarly, using the continuity of shear forces and moments across the delamination tips, and the beam theory relationships, the shear force and bending moment response at the delamination tips can be represented as a function of both sets of coefficients, written as: where the matrices are functions of the problem geometry and continuity conditions, related to the coefficients using beam theory relationships and (4). Using the relationships given in expressions (12) and (13), the force vector can be written as a function of one set of coefficients only (i.e., or ). Here, the choice was made to have the top beam’s coefficients, , as the reference parameters. Consequently, from expressions (12) and (13) one can write: furthermore, from (4), the end displacements and slopes can be related to coefficient vector, , through the following expression: Finally, using expressions (14) and (15) leads to where is the frequency-dependent, dynamic stiffness matrix of the system. The standard assembly process similar to FEM leads to the nonlinear eigenvalue problem of the system: where is the overall (global) dynamic stiffness matrix, and represents the vector of defrees of freedom of the system. The solution of the problem consists of finding the eigenvalue, , and corresponding eigenvector, , that satisfy (17) and the boundary conditions imposed using, for example, the penalty method [22]. Powerful algorithms exist for solving a linear eigenvalue problem (i.e., system’s natural frequencies), resulting from discrete or lumped mass models. In the case of the nonlinear eigenproblem (17), involving frequency-dependent dynamic stiffness matrices arising from the DFE or DSM formulations, one can use the Wittrick-Williams (W-W) root-finding technique [17] to determine the eigenvalues of the system. The W-W algorithm is a simple method of calculating the number of natural frequencies of a system that are below a given trial frequency value. The method exploits the bisection method and the Sturm sequence properties of the dynamic stiffness matrix to converge on any particular natural frequency of the system, to any desired accuracy. This allows one to solve for any specific frequency number without having to solve for all previous frequencies, which is the requirement of some linear eigenvalue solvers. Consequently, the corresponding modes can be evaluated [10–17, 24].
3. Numerical Tests
Numerical checks were performed to confirm the predictability, accuracy, and practical applicability of the proposed DSM method. DSM and FEM formulations, as well as the Coefficient Method (CM), were programmed in Matlab codes. To solve the nonlinear eigenproblem (17) resulting from DSM formulation, a determinant search method was used; the nondimensional frequency was swept, searching a particular frequency, , which would make the determinant of the global dynamic stiffness matrix zero, , whose corresponding eigenvector, , represented the degrees of freedom of the mode shape associated with the natural frequency. The linear eigenvalue problem resulting from the conventional FEM formulation, was solved using Matlab “eig” function. The use of the nondimensional frequency (5) in the calculations removed material dependencies from the system, provided that the material was isotropic, or at least orthotropic with principal axes aligned with the Cartesian coordinate system in Figure 1.
In what follows, an illustrative example of fixed-fixed, homogeneous, 2-layer delaminated beam is examined. The natural frequencies of the system with a central split, about the midsection (), of various lengths up to 60% of the span , occurring symmetrically along the midplane of the beam and surrounded by intact beam segments, are considered. This split beam configuration has also been presented and studied in [1, 2, 8, 9]. The DSM was used to compute the natural frequencies and mode shapes of various delamination cases. As the benchmarks for comparison, the results from references [1, 8] —and reference [2] for the constrained mode—were used to validate the solution method presented here. As also suggested in [1], the first few frequencies were computed for a delamination length of 0.0002 and showed negligible discrepancies from those of a solid intact beam to check for numerical instability when the split length becomes extremely small. The effect of the longitudinal motion of the upper and lower parts of the split region on the frequencies, examined in [1], has been neglected here for this class of example problems.
Table 1 summarizes the DSM results for the first two natural frequencies of the system. The DSM results are compared to those presented by Wang et al. [1] and Della and Shu [8] and Erdelyi and Hashemi [9]. The DSM model incorporates a total of only three elements; one intact element on each end of the delamination representing the undamaged beam segments (1 and 4) obtained using the methods outlined in [25, 26], and one fully delaminated element (Figure 1). The DSM natural frequencies are in excellent agreement with those reported in references [1, 7], with a maximum difference of 0.24% (see the last row in Table 1), which could be simply attributed to the extensive use of matrix manipulations in both DSM and CM methods (the results obtained from the CM-based code developed by the authors [23] are not reported here, as they were found to be in perfect match with tabulated frequency data and graphs reported in the literature). It is also worth noting that the authors found a slight dissimilarity between the 2nd natural frequency values (61.67, 60.76, 55.97, 49.00, 43.87, 41.45, 40.93, resp.) reported in Table 1 of [1] and the same data appearing in 5th column of Table 1 of [8], cited and reported from [1] (see the 5th column of Table 1). Conventional FEM natural frequencies obtained based on layerwise theory, as reported by Lee [6], are also presented for comparison (last two columns in Table 1). Excellent agreement was found between the DSM and the FEM results.
A split beam FEM, exploiting cubic Hermite [22] interpolation functions, was also developed [23]. The weighted residual method is applied on the differential equations (3), governing the free vibration of 2-layer delaminated beams. The residual was made orthogonal to a virtual displacement over the domain of the element, and two integrations by parts were carried out to reduce the continuity requirements of displacement functions. The principle of virtual work was used to determine the element system equations. As presented earlier, the differential stretching of the top and bottom layers should be present to keep the delamination faces planar after deformation (i.e., no interlaminar slip at the delamination faces). The FEM formulation results in an additional stiffness term not present if interlaminar slip were included. This “delamination stiffness” has the effect of stiffening the system at the delamination tips (for more information on the split beam FEM, the reader may refer to [23]). Table 2 summarizes the first two natural frequencies obtained using the developed (cubic Hermite-type) Finite Element Model (FEM), with 6- and 10-element discretizations of midplane delaminated region (60% of span). The intact beam segments were modeled using single-beam elements. As can be seen from Table 2, the FEM frequencies exhibit a convergence towards the DSM results, as the number of elements is increased. Conventional FEM natural frequencies reported by Lee [6] are also presented for reference.
Figure 3 shows the first two natural modes of the 2-layered beam, with 60% of span midplane delamination, compared with those of an intact configuration. It is worth noting that the conventional FEM-based models are characterized by constant mass and stiffness matrices of limited number of total degrees of freedom (DOF), that is, number of nodes times number of DOF per node. Accordingly, the natural modes obtained from the conventional FEM model—being the eigenmodes of the governing linear eigenvalue problem—have the same dimension as the total degrees of freedom of the FEM model. Unlike the conventional FEM (e.g., 4-DOF Hermite beam element), the DSM and frequency-dependent FEM matrices are formulated based on continuous element assumptions, which introduces infinite number of degrees of freedom within each element (see, e.g., [11–17, 25, 26]). Therefore, through the use of these techniques, additional modes of vibration can be found. These modes are the result of the denominator of the global stiffness matrix going to zero, and correspondingly the determinant of the global stiffness matrix approaching infinity, . Also known as the poles of a system, they can represent real physical mode shapes, describing the structure vibrating at zero nodal displacements [18, 25] outside of the delaminated region. Through simplification, it was found that the denominator of the stiffness matrix, in this case, has the following form: While the mode shapes of the poles were not analytically important in this analysis, their corresponding natural frequencies are important when using more advanced root solving techniques [18, 25]. Zero-nodal-displacement modes have also been observed and reported in the literature for other structural configurations (see, e.g., [18, 24–26]). There are also certain frequencies captured through the system modal analysis whose mode shapes, while mathematically possible, do not represent physically admissible displacements. These modes—for example a second mode in the case of present study—are simply the result of the free model assumptions [5]. They correspond to interpenetration of the beams, as illustrated in Figure 4, and would not be present in a constrained mode analysis. As seen in Figure 4, the vibration of the top and bottom delaminated beams would be inadmissible due to non-linear phenomena such as contact, which cannot be modeled in the frequency domain. Similar inadmissible partial and complete interpenetration modes have also been reported in the literature [27]. In addition to real natural modes of vibration, poles and inadmissible interpenetration modes examined above, under small vibration amplitudes a split layered beam may exhibit a mode at a frequency corresponding to a delamination-opening mode. Figure 5 shows the first opening mode for a delaminated beam with top beam thickness equal to 40% the height of the intact beam, 60% of span, off-midplane delamination, obtained using 3-element DSM and FEM models (FEM nodes visualized). Similar opening modes have also been reported in the literature (see, e.g., [7, 8]).
(a)
(b)
4. Conclusion
Based on the “exact” dynamic stiffness matrix (DSM) formulation, a new element for the free vibration analysis of a delaminated layered beam has been developed using the free mode delamination model. The DSM element exploits the closed form solution to the governing equation of the system and is “exact” within the limitations of the theory. For homogeneous beams with a central, midplane delamination, a 6-element model of the delaminated system provides excellent agreement with those models presented in the literature. A conventional finite element model was also briefly discussed. System natural modes, pole behaviour and opening modes for both midplane and off-midplane delaminations were also examined and illustrated.
Acknowledgments
The authors wish to acknowledge the support provided by NSERC, Ontario Graduate Scholarship (OGS), and Ryerson University, as well as the reviewers for their helpful comments.