Abstract
The dynamic analysis of flexible delaminated layered beams is revisited. Exploiting Boolean vectors, a novel assembly scheme is developed which can be used to enforce the continuity requirements at the edges of delamination region, leading to a delamination stiffness term. The proposed assembly technique can be used to form various beam configurations with through-width delaminations, irrespective of the formulation used to model each beam segment. The proposed assembly system and the Galerkin Finite Element Method (FEM) formulation are subsequently used to investigate the natural frequencies and modes of 2- and 3-layer beam configurations. Using the Euler-Bernoulli bending beam theory and free mode delamination, the governing differential equations are exploited and two beam finite elements are developed. The free bending vibration of three illustrative example problems, characterized by delamination zones of variable length, is investigated. The intact and defective beam natural frequencies and modes obtained from the proposed assembly/FEM beam formulations are presented along with the 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. Delamination, a common failure mode in layered structures, 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, in turn, will affect the vibration characteristics of the structures. 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.
Several experimental methods exist to predict the onset, size, and growth of delamination as a failure mode in composite materials. Using acoustic emission (AE) sensors, different levels of amplitude signals emitted by the materials can be monitored [1]. The different amplitudes correspond to loading types and are assigned damage mechanisms. Using this technique, continuous monitoring of damage is possible experimentally. Another research [2] has shown that acoustic emission is a viable and effective tool for identifying damage and distinguishing damage types in self-reinforced polyethylene composites. More recently, further research has been done using neural networks and unsupervised learning techniques applied to the data set of acoustic emission signals [3]. The signals were successfully used to classify the AE patterns caused by different damage mechanisms in carbon-reinforced composites (delamination and matrix cracking).
The time-domain stability of vibrating delaminated systems has also been an area of study. Particularly, the mechanisms that cause the different laminates to separate from each other are typically not defined in most theoretical applications. They are described as physically inadmissible mode shapes whose existence in the frequency domain is a product of delamination tip boundary conditions. The study of the real phenomenon, however, has shown that time-dependent normal forces in the delaminated segments do not influence the global free vibration frequencies but may contribute to localized buckling [4, 5]. Instability and critical dynamic forces can be predicted, allowing for study of the onset of delamination opening.
The vibration modelling and analysis of delaminated multilayer beams has been a topic of interest for many researchers. The earliest delamination models formulated in the 1980s [6] addressed the vibration of two-layer beams, where each layer was modelled using Euler-Bernoulli bending beam theory. The upper and lower portions of the delaminated segment were assumed to vibrate independent of each other, that is, “free mode” delamination. The free mode, however, underpredicts natural frequencies for off-midplane delaminations due to unrestricted penetration of the beams into each other. In 1988 [7], Mujumdar and Suryanarayan proposed the “constrained mode” delamination model, which assumes equal transverse displacements for the top and bottom beams and the rigid connector assumption. The rigid connector assumption states that, for the beam models presented, the delamination faces, which are planar and normal to the neutral axis of the undeformed beam, remain planar and normal to the neutral axis of the deformed beam. This assumption produces a set of kinematic and force continuity conditions at the delamination tips. In a recent work by Szekrényes [8], an extensive literature survey of the research works related to the vibrations of delaminated elements was presented and, based on coupled flexural-longitudinal vibration model, 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 that delamination buckling can take place if the normal force is compressive in one of the half periods of the vibration and reaches a critical value [8].
The constrained mode delamination model, predicting vibration behaviour much more accurately for off-midplane delamination, is in fact simply a limiting case of the free mode delamination model. However, opening delamination modes, that is, where the layers separate from each other, commonly seen in experimental analysis [9–11], cannot be captured using the constrained model. Therefore, in the present study, the free mode delamination model will be investigated and the constrained mode delamination model can be derived in a similar manner.
The accuracy of dynamic/forced response analysis 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 predict the natural frequencies and mode shapes of such a system. Several exact solution methods exist for well-defined systems, such as delaminated isotropic beams with constant geometric and material properties. Single [9, 12, 13], multiple [14], and various overlapping and enveloped delamination conditions in space and on various elastic media, such as Pasternak soil [15], have been studied using analytical solution methods. Some work has also been done on delaminated sandwich structures [16], albeit with some mathematical simplification. These solution methods generally use the same procedure as Mujumdar and Suryanarayan [7] to formulate the kinematic continuity conditions across the delamination tips. The power of this type of formulation lies in the ability to be applied to any number of different system configurations. However, a potential drawback to this procedure is that the system equation must be reformulated after any configuration change, potentially limiting its applicability.
The conventional Finite Element Method (FEM) has a long, well-established history and with the advent of digital computers it is commonly used for structural analysis. The FEM is a general and systematic approach to formulate the element matrices for a given system and is easily adaptable to complex systems, such as nonuniform geometry, often modeled as a stepped, piecewise-uniform configuration. Exploiting polynomial interpolation (shape) functions, the FEM leads to constant element mass and stiffness matrices and ultimately a linear eigenvalue problem from which the natural frequencies and modes of the system can be readily extracted. The FEM method for a single beam can be modified to accurately model delaminated multilayer beams. Among others, Lee [17–19] used the layerwise FEM theory to investigate the free vibration of delaminated beams. In the recent years, layered, sandwich, and composite elements have been integrated in certain commercial software and are used to analyze the vibration of composite structures. However, modelling a delaminated configuration in commercial software packages such as ANSYS is not straightforward and 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) [20] to enforce the displacement and slope continuity at the edges of delamination region [21].
Semianalytical formulations, such as the so-called Dynamic Finite Element (DFE) method [22], have also been developed 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 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 the preliminary free vibration analysis of a delaminated 2-layer beam has been reported in an earlier work by the authors [23].
Analytical methods, namely, the Dynamic Stiffness Matrix (DSM), have also been used for the vibrational analysis of isotropic, sandwich, and composite structural elements and beam-structures. The DSM approach exploits the general, closed-form solution to the governing differential equations of motion of the system to formulate a frequency-dependent stiffness matrix. The DSM produces exact results for simple structural elements, such as uniform beams, and Banerjee and his coworkers [24, 25] have developed a number of DSM formulations for various beam configurations. The DSM method for a single beam can be modified to accurately model delaminated multilayer beams. A DSM-based analysis of a two-layer split beam has also been presented in earlier works by the authors [26, 27].
The aim of this paper is to present an FEM formulation for the linear, 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, and shear deformation and rotary inertia are neglected. Continuity conditions for forces, moments, displacements, and slopes at the delamination tips are enforced through a novel Boolean vector assembly scheme, leading to the integral FEM model of the system. In fact, through the presented method, one obtains specific matrices for an intact, fully delaminated, and delaminated elements attached to an intact segment from left/right. Therefore, a direct assembly method can be directly used to form various multiple-delaminated beam configurations, without the need to manually create the model and to use a constraint element (e.g., ANSYS element MPC184 [20]), to enforce the displacement and slope continuity at the edges of delamination region. Thence, the direct assembly of element matrices and the application of system’s global boundary conditions results in the linear eigenvalue problem of the defective system. In addition, two MATLAB-based computer codes, based on the Dynamic Stiffness Matrix (DSM) method [26, 27] and the analytical solutions reported in the literature [28–31], are developed and used as a benchmark for comparison. Two 2-noded and 3-noded beam FEM elements are presented, where cubic Hermite and quartic interpolation functions of approximation, respectively, are used to express the flexural displacement functions, that is, field variables and weighting functions [32]. The FEM models are used to compute the natural frequencies of an illustrative defective beam example, characterized by a single delamination zone of variable length. The frequency values are then compared with DSM data and those from the literature. Certain modal characteristics of the system are also discussed. It is worth noting that while the model used in this study assumes isotropic materials, further research is underway to extend it to sandwich [33] and fibre-reinforced laminated composite beams, characterized by an extensional response coupled with flexural/torsional and coupled bending-torsion vibration [34, 35].
2. Mathematical Model
Figure 1 shows the general coordinate system and notation for a single-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 [26–28, 31]For harmonic oscillations, the transverse displacements can be written aswhere is the amplitude of the displacement , subscript “” represents the beam segment number, and is the circular frequency of excitation of the system. Back-substituting (2) into (1), the equations of motion reduce toThe 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 asCoefficients , , , and are evaluated to satisfy the displacement continuity requirements of the beam segments and the system boundary conditions. As also reported by several researchers [8, 26–28, 31], the inclusion of delamination into the beam model leads to coupled axial-transverse motion of the delaminated beam portions, primarily associated with the continuity requirements imposed at the delamination endpoints. In order for the delamination tip cross sections 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. The midplanes, that is, the neutral axes of the beam segments, in the delaminated region are located at a distance from the midplanes of the intact segments. Hence, they will not have the same axial deformation unless some internal axial force is imposed. As mentioned earlier in this paper, based on coupled flexural-longitudinal vibration model [8], the equality of axial forces in the top and bottom beams has been recently derived and shown in an exact way and the continuity of the effective bending moments was related to the equilibrium equations. However, in what follows, this imposed axial force is briefly presented for completeness, following the method derived and discussed in [7].
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 (see Figure 2). Additionally, the requirement that the delamination tip faces remain planar after deformation results in, at the left delamination tip,where is the axial displacement of beam section and from the kinematic continuity conditions. If this is combined with the same formulation from the right delamination tip,The assumption is made by Mujumdar and Suryanarayan [7] and by other researchers (see, e.g., [27, 31]), that the axial displacement will behave according to the following, for small deformations and material and geometric properties which remain constant along the length of the beam:where is the axial stiffness of beam section . Substituting this into (7) yieldsUsing the continuity of axial forces across the delamination tip, ,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 asExpression (11) can be further simplified if the cross-sectional shape is known, and continuity of bending moments at the left and right delamination tips, respectively, leads to the following equations:Exploiting (12a) and conditions ((10), (11)) for internal axial force, and noting that bending moments and shear forces in beam segment “” are related to displacements, , through , and , respectively, it can be shown that, for continuity of bending moments,where the coefficient is defined asTo satisfy the continuity of shear forces about the left delamination tip,
Likewise, using (12b) for , relevant relationships similar to (13) and (15) can be derived for the left delamination tip. Two boundary conditions at the intact beam ends, continuity of displacements and slopes, at the delamination tips result in 12 equations. Along with 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 appear in Expression (4). In other words, the 16 equations can be solved simultaneously, using a root finding algorithm to find the natural frequencies and mode shapes of the system. Thus, an analytical solution can be produced for each set of system’s imposed boundary conditions (see, e.g., [28–31]). This solution method, based on finding the coefficient matrix of the system, herein referred to as the “Coefficient Method (CM),” has been used to predict the vibration behavior of different systems of varying complexity (see, e.g., [28–31]). However, one of the main drawbacks of this method is that it remains a relatively problem-specific solution technique, which should be reformulated every time the system’s global boundary conditions change. In an earlier publication by the authors [27], an equivalent, yet more conveniently applicable Dynamic Stiffness Matrix (DSM) formulation was presented. The DSM, however, lacks generality and is limited to uniform beam configurations. Therefore, in what follows a general FEM-based model is presented, which can be readily extended to more complex cases, with variable geometric and material parameters.
2.1. Finite Element Method (FEM) Formulation
The finite element approach used here is based on the Galerkin method of weighted residuals. The equations of motion for each beam are used as the basis of this solution method. Simple harmonic motion is again assumed, and the equations of motion, according to the Euler-Bernoulli beam theory, take the following form:where is the actual transverse displacement of beam , and the same nondimensionalization used in (5) has been applied. An approximate transverse displacement is introduced in place of the actual displacement, such that . This results in the following residual equation:where is the residual of the approximate equation. Following the Galerkin method of weighted residuals, the residual above is weighted by a virtual displacement and the integral is set to zero across the domain of the system. Since the system is composed of four distinct beam sections occupying their own subset of the domain, the following is representative of the Galerkin method applied to the delaminated system:withwhere are the shape functions of the beam elements, which will be defined later. Since the virtual displacement is applied to the entire domain, and the four different beam sections occupy unique subdomains, . In order to produce the force and displacement continuity terms, a set of integrations by parts is performed on the above, resulting in the following weak form:The terms in , above, represent the boundary and continuity conditions imposed on the system. Using Euler-Bernoulli beam theory, the shear force and bending moment at any point are defined based on the transverse displacement asFor the endpoints of beam sections 1 and 4, the following is true for free vibration:where is the external virtual work caused by applied external forces on the system, causing virtual displacements. For the free vibration of this system, the total external work is . The remaining terms in above can be resolved by applying the necessary continuity conditions for displacements (e.g., , slopes: , and shear forces: ), leading toThe terms and in (23), as well as the external work term, go to zero directly as a result of the shear force continuity conditions. However, the remaining terms do not vanish, since the continuity of bending moments ((13), (14)) contains an additional implicit bending-axial coupling term, such thatWith the boundary and continuity conditions satisfied, the system can be discretized into elements, which will each be approximated using their own basis functions, from which FE shape functions can be found. The system can be discretized as follows, using the result of (24):where “” is the number of elements in beam section .
It is worth noting that following the same above-described procedure, expressions (13) and (14), respectively, for the case of double-delaminated three-layer beam configurations can be written aswhere subscripts stand for the intact beam (left and right extremity) segments and the three delaminated layers are represented by . The necessary displacements, slopes, and shear force continuity conditions at the delamination tips, respectively, are then written as
2.2. 2-Node Beam Element (#1)
Following the traditional Euler-Bernoulli finite element development, Hermite cubic polynomials [32] were used as the basis functions of approximation for each beam, such that, for a two-node, 2-degree-of-freedom per node beam element, that is, transverse displacement and slope defined at each node (Figure 3),where is a column vector of unknown constant coefficients. The following represents the vector of nodal displacements used in further FE development:Thuswhere is a row vector of shape functions, which describe the displacements at any point along the domain of the element in terms of the nodal displacements and slopes at the endpoints of the element domain, . Additionally, the shape functions may also be used to approximate the virtual displacements, . With the shape functions fully defined, they may be substituted for the approximate displacements in (25)Frequency-dependent and non-frequency-dependent terms above can be gathered to form the following eigenvalue problem, common to structural vibration analysis with FEM, with a modification caused by the presence of the delamination:ifwhere is the structural stiffness matrix formed by assembling the associated beam elements, as per (31), is the delamination stiffness matrix, from the term appearing outside the integral expression in (31), and is the structural mass matrix. From this formulation, the simplest solution methods involve eigensolutions.
2.3. 3-Node Beam Element (#2)
Exploiting the same concept of a polynomial interpolation function, a 3-node, 5-DOF element is also developed (see Figure 4). Making use of a higher-order polynomial interpolation functions increases the accuracy of the solution. Whereas for the 2-node beam element (Figure 3) a 3rd-order polynomial was required, for a higher-order interpolation of 4th-order one requires the addition of another single degree of freedom to the system. This was accomplished by adding a midpoint node with only one degree of freedom (lateral displacement) to the beam model used previously. This third node, while increasing the mesh fineness, allows for a greater solution accuracy and possibly faster convergence, which will be investigated. The 3-node beam element (Figure 4) was developed in the same way as the 2-node beam element, except using the following interpolation function:Consequently, the degrees of freedom for the system were modified as discussed. The addition of the midpoint node and its associated lateral degree of freedom are compensated for by using the following degrees of freedom:Thus
3. Numerical Tests
Numerical checks were performed to confirm the predictability, accuracy, and practical applicability of the proposed delaminated FEM models. Both the FEM and Coefficient Method (CM) [31] formulations were programmed in MATLAB codes. The linear eigenvalue problem resulting from the conventional FEM formulations 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, two illustrative examples of fixed-fixed, homogeneous, 2-layer, and 3-layer delaminated beams are examined. In the first example, 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 [6, 7, 27–31]. The split FEM and DSM models were created and used to compute the natural frequencies and mode shapes of various delamination cases. As the benchmarks for comparison and validation purposes, the results from [6, 7, 27] for the constrained mode, as well as an alternative formulation from [6] were used. As also suggested in [6], the first two frequencies were computed for a delamination length of , to check for numerical instability when the split length becomes extremely small. This case showed negligible discrepancies from those of a solid intact beam. The effect of the longitudinal motion of the upper and lower parts of the split region on the frequencies, examined in [6], has been neglected here for this class of example problems. As discussed earlier in this paper, the differential stretching of the top and bottom layers was 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, called “delamination stiffness,” which has the effect of stiffening the system at the delamination tips.
Tables 1 and 2 summarize the first two natural frequencies obtained using the beam element #1 (4-DOF) Finite Element Method (FEM), with 6- and 10-element discretizations of midplane delaminated region (up to 60% of span). The intact beam segments were modeled using single beam elements. As can be seen from Tables 1 and 2, the FEM formulation produced excellent agreement with analytical results and the layerwise FEM theory [17] and exhibits a convergence towards the DSM results, as the number of elements is increased. While the FEM result discrepancies are generally low, even for a coarse mesh size, it is observed to be even lower for the first natural mode than the second one. This is consistent with traditional FEM theory, where more elements are required to guarantee accurate solutions for higher mode numbers.
Figure 5 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 matrices are formulated based on continuous element assumptions, which introduces infinite number of degrees of freedom within each element (see, e.g., [24, 25]). Therefore, as also reported in [27], using DSM technique, 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 [36] outside of the delaminated region. Zero-nodal-displacement modes have also been observed and reported in the literature for other structural configurations (see, e.g., [22, 36]). 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 (λ = 31.0) in the case of present study, are simply the result of the free model assumptions [27]. They correspond to interpenetration of the beams, as illustrated in Figure 6, and would not be present in a constrained mode analysis. Similar inadmissible partial and complete interpenetration modes have also been reported in the literature [37]. In addition, as also reported in the literature (see, e.g., [27–31]), under small vibration amplitudes a split layered beam may exhibit a mode at a frequency corresponding to a delamination-opening mode. The first opening mode for a delaminated beam with top beam thickness equal to 40% of the height of the intact beam, 60% of span, off-midplane delamination, obtained using a 17-element FEM model, is depicted in Figure 7.
(a)
(b)
In order to further assess the accuracy of the proposed FEM method, in what follows two different 3-layer beam configurations with central dual delaminations are analyzed (refer to Figure 8).
(a)
(b)
The FEM frequency results, obtained using both 2- and 3-node elements for both models above, are validated against various analytical solutions obtained from a MATLAB code written based on the “exact/analytical” (CM) method [31], those extracted by interpolation from a graph by Della and Shu [31], and the frequency data obtained from a DSM code (see an earlier work by the authors [27]), as well as those gathered from the literature. For the relatively simple models presented in this section, the discretization (Figure 9) was used to mesh the domains.
In the first delamination model tested (see Figure 8) the top and centre delaminated beams each have a height of 30% of the intact beam height. In addition, the delamination length, a, was varied as a percentage of the total beam length from 20% to 50%. The delamination is central, meaning that the left and right intact segments have equal lengths. The frequency results for these delaminated clamped-clamped beam configurations are presented in Table 3.
As it can be observed from Table 3, both 2- and 3-node beam elements perform well with respect to both the analytical solution and those taken from the literature. As expected, light deviations (0.26% for 2-node mode 2, with respect to the exact/analytical solution) are present for larger delamination sizes and for higher modes of vibration. For higher modes, the 3-node beam tends to perform slightly better than the 2-node beam. This difference is expected to increase with an increase in mode number. However, for the first mode and for small delamination sizes, the 2- and 3-node beam elements perform similarly and differences between the results for each element were negligible. This could be used to justify the use of a 2-node beam element in such situations to save on processor requirements and solution times.
Since identifying nonphysical or physically inadmissible behaviour can be essential in improving a given model, the system’s natural modes are also an important consideration in free vibration modeling. Therefore, the first two mode shapes were compiled in Figure 10, generated from the system’s eigenvectors results, the shape functions derived and presented earlier. Since the variance between mode shapes for the 2- and 3-node beam elements was minimal, only the mode shape for the 3-node beam will be plotted. Note that, in Figure 10, the element end nodes are represented by circular markers, and the midpoint nodes are represented by × markers.
In the second delamination model tested (see Figure 8), the total thickness is the same as in Model 1; however, each beam segment’s thickness is assumed to be one-third (1/3) of the total thickness. As before, the results reported from [31] were converted from chart to numerical form. The final results of this second analysis are detailed in Table 4.
As it can be observed from Table 4, for smaller delamination sizes, the 2-node and 3-node beam elements performed with similar accuracy, and both converged to a reasonable accuracy (largest deviation approximately 0.16% from the exact solution). However, for higher modes and larger delamination sizes, the 3-node beam exhibits higher accuracy and better convergence characteristics. It is also worth noting that the FEM results presented in this section were obtained with only slight modifications to the single delamination technique presented earlier, whereas the analytical solution had to be completely modified. This lends further credence to the advantages of presented FEM-based formulations, since good correlation with analytical results can be achieved with less development overhead.
The system’s first two natural modes evaluated in the same manner as in Model 1 are presented in Figure 11. Furthermore, some mode shapes emerged from the analysis, which involved physically inadmissible mode shapes. In this case, the inadmissibility comes from the interpenetration of different beam layers with each other (see Figure 12). That is, one beam segment would vibrate laterally in one direction and another beam segment, occupying the same axial domain, would vibrate laterally in the opposite direction. In the present study, the physically inadmissible mode shapes were found to occur when the difference in flexural stiffness of the beams is nonzero and worsens with increasing difference, seen for both exact and FEM solutions. It was observed that if the difference in beam stiffness between the three delaminated beams was sufficiently large, these modes would appear to be slight interpenetrations.
(a)
(b)
Delaminations not centrally located along the beam (i.e., ) were observed to generate mode shapes deviating from standard beam mode shapes, and opening of the beam segments can be seen (see Figure 13). The delaminated beam segment with the smallest bending stiffness deforms the largest and will thus be forced into the other, stiffer beam segments. Although not physically admissible, there is no provision within the free model to prevent this behaviour. This phenomenon was observed for all doubly delaminated beam models with disparate delaminated beam heights (), as well as for any noncentral delamination ().
(a)
(b)
Finally, it is also worth mentioning that, by proving the existence of the parametric excitation in delaminated beams [4], it was recently shown and experimentally verified that the opening is amplitude dependent; that is, the delamination opening takes place only at a certain critical amplitude. In reality, 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 [8]. Furthermore, the experimentally observed delamination opening was reported to be significantly less than those calculated by the free model, which would rather justify the use of the constrained model for further analysis [8].
4. Conclusion
Based on the conventional Finite Element Method (FEM) formulation and employing Boolean vectors, a novel assembly scheme for the free vibration modelling and analysis of a delaminated layered beams was developed. Based on the free mode delamination model, 2- and 3-node beam elements equipped with cubic and quartic shape functions, respectively, were used to extract the flexural natural frequencies of single- and double-delaminated beam configurations. While the results showed very good agreement with those obtained from analytical and Dynamic Stiffness Matrix (DSM) models presented in the literature, the use of 3-node beam elements with quartic shape functions, as opposed to the more common Hermite cubic shape functions found in 2-node beams, did not produce significantly better results for the same mesh size. System’s natural modes and opening modes for both midplane and off-midplane delaminations were also examined and illustrated. The constrained model can be applied to prevent the interpenetration modes, at the expense of increased system stiffness [7, 31].
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Authors’ Contribution
The authors declare that this paper presents the results of a research conducted by the first author, under the supervision of the second (corresponding) author.
Acknowledgments
The authors wish to acknowledge the support provided by The Natural Sciences and Engineering Research Council of Canada (NSERC), Ontario Graduate Scholarship (OGS), and Ryerson University.