Free Vibration Analysis of Fibre-Metal Laminated Beams via Hierarchical One-Dimensional Models

This paper presents a free vibration analysis of beamsmade of fibre-metal laminated beans. Due to its attractive properties, this class of composites has gained more and more importance in the aeronautic field. Several higher-order displacements-based theories as well as classical models (Euler-Bernoulli’s and Timoshenko’s ones) are derived, assuming Carrera’s Unified Formulation by a priori approximating the displacement field over the cross section in a compact form.The governing differential equations and the boundary conditions are derived in a general form that corresponds to a generic term in the displacement field approximation. The resulting fundamental term, named “nucleus”, does not depend upon the approximation order N, which is a free parameter of the formulation. A Navier-type, closed form solution is used. Simply supported beams are, therefore, investigated. Slender and short beams are considered. Threeand five-layer beams are studied. Bending, shear, torsional, and axial modes and frequencies are presented. Results are assessed for three-dimensional FEM solutions obtained by a commercial finite element code using threedimensional elements showing that the proposed approach is accurate yet computationally effective.


Introduction
In Fibre-Metal Laminates (FMLs), Figure 1, metallic and Fibre-Reinforced Polymers (FRPs) are synergistically merged such that improved properties with respect to the individual constituents are obtained.Thanks to the presence of the metallic layer, FMLs present higher fracture toughness, yield stress, and compressive strength when compared to conventional FRPs [1].Resistance to impact and environmental solicitations is improved by placing a metallic layer at top and the bottom of the structure.FMLs also offer several advantages upon monolithic metallic laminates such as high strength-and stiffness-to-weight ratios and fibre-bridging mechanism against crack growth in the metallic layers.In the case of heating due to fire, monolithic aluminium skin layers last about 20 to 30 seconds while the metallic layers in FMLs do not melt away and maintain a residual load carrying capability because the fibres may act as firewall [2].Although FMLs per kilogram can cost five to ten times more than aluminium laminates, the cost of a finished product can be very close, even disregarding the positive effect of FMLs on direct operative costs (such as maintenance).Saving in weight of the structure (about 20%) can be obtained [3].FMLs idea has been conceived in the late 70s at the Delft University of Technology when trying to lower the cost of full composites while maintaining the weight saving they comported.A breakthrough was obtained by inserting reinforcing fibres to the adhesive of already available bonded metallic sheets.In 1982, aramid reinforced aluminium laminates (ARALL) were commercially available.Unidirectional glass fibre embedded reinforced FMLs (GLARE) followed [3].Besides aramid or glass fibres, FRPs made of unidirectional carbon (CARALL), vinylon (Virall), and both carbon and Kevlar or jute natural fibres in polymeric matrices have been presented.Woven fibre-reinforced epoxy plies have been also considered.The metallic layers (whose thickness range is about 0.2 to 0.5 mm) are often made of aluminium alloys such as 7075, 7475, or 2024.Other possible solutions account for magnesium or titanium.Mechanical or adhesive bonding is, then, used to connect the metallic and FRPs layers.FMLs have been conceived to meet the need of reduced weight and high performances in aeronautics.Fokker was one of the first to investigate the use of FMLs as panels for the wing of the F-27.Following this example, several companies (e.g., Airbus, Boing, and Embraer) started to use FMLs [4].For instance, FMLs are used in the fuselage skins of the Airbus A380.Titanium-carbon FRPs solutions have been conceived for high temperature applications in supersonic transport and space applications.Over the years, different applications have been considered on FMLs tubes for chemical and nuclear industry, electronics packaging and fuel tanks, and blast resistant containers as well as sport goods (e.g., bike frames, tennis rackets).Botelho et al. [4] presented a review of FMLs material properties, whereas Wu and Yang [5] thoroughly reviewed GLARE properties.Van Rooijen et al. [6] investigated the influence of constituents properties on FMLs mechanical properties by means of a rule of mixture in terms of the metal volume fraction.Cortes and Cantwell [7] studied numerically and experimentally the tensile properties of a lightweight FML based on a titanium alloy and carbon fibre-reinforced FMLs.Sinke [8] presented a review of several GLARE structural solutions as well as manufacturing.FMLs offer a wide variety of choices due to the high number of involved parameters and possible configurations [9].Due to the presence of FRPs layers, both the structure and the laminates have to be accounted for in the design process.Alderliesten [10] suggested that a top-down design approach (from structural requirements down to individual constituents) should be used.Sinke [8] outlined the need for multiscale modelling and design.Sawant and Muliana [11,12] used this approach for the thermoviscoelastic behaviour of FMLs.
As laminated structures, advanced models developed for laminates can be used [15] provided they are extended to account for the hybrid behaviour deriving from the metallic layers [8].As outlined in Alderliesten [10], FMLs mechanical response is generally predicted by the classical laminate theory [16].Teply et al. [17] used a linear layerwise plate element for the analysis of ARALL laminates.A geometrical and material nonlinear solid-like shell element was formulated by Hashagen et al. [18].Hagenbeek [19] proposed a solidlike shell element for the thermomechanical problem and results were assessed for experimental test.Hausmann et al. [20] presented some analytical and numerical solutions for determining the thermal residual stresses due to processing.
A free vibration analysis of FML beams via hierarchical one-dimensional models is addressed in this paper.Models are derived via Carreras Unified Formulation that was previously derived for plates and shells [21][22][23][24][25] and extended to beams [26][27][28][29][30][31].Through a concise notation for the displacement field, the governing differential equations and the corresponding boundary conditions are reduced to a "fundamental nucleus" that does not depend upon the approximation order.This latter can be assumed as a formulation free parameter.Displacement-based theories that account for nonclassical effects, such as transverse shear and cross section in-and out-of-plane warping, can be formulated.Governing differential equations are solved via a Navier's closed form solution.Slender as well as short beams are investigated.Three-and five-layer FML beams are studied.The presented models are validated towards threedimensional FEM solutions obtained by commercial code Ansys.

Preliminaries
In presenting the proposed model, the methodology in Giunta et al. [32] has been followed.A beam is a structure whose axial extension () is predominant if compared to any other dimension orthogonal to it.To identify the cross section (Ω), it suffices to intersect the beam with planes that are orthogonal to its axis.A Cartesian reference system is adopted: and -axes are two orthogonal directions lying on Ω.The  coordinate is coincident to the axis of the beam and is bounded such that 0 ≤  ≤ .Cross section geometry and reference system are reported in Figure 1.The displacement field is u (, , ) = {  (, , )   (, , )   (, , )}  (1) in which   ,   , and   are the displacement components along -, -, and -axes.Superscript '' represents the transposition operator.The strain, , and stress, , tensors in Voigt notation are grouped into two vectors:   ,   lying on planes orthogonal to Ω where and   ,   lying on the cross section: In the case of a linear analysis, the following straindisplacement geometrical relations hold: Subscripts '' , '' , and '' , when preceded by comma, represent derivation versus the corresponding spatial coordinate.A compact vectorial notation can be adopted for (4): where D Ω , D  , and D  are the following differential matrix operators: ] and I is the unit matrix.Under the hypothesis of linear elastic behaviour of the material, the generalised Hooke law holds: in which C is the material stiffness matrix.According to (2) and (3), it reads Matrices C  , C  , C  , and C  in (8) are ] where terms   are the material stiffness coefficients: where and   are Young's moduli, ]  Poisson's ratios, and   the shear moduli, respectively.

Hierarchical Beam Approximation
The displacement field is, a priori, assumed over the cross section in the following way: represents the number of unknowns while   is a generic expansion function over the cross section.The compact expression is based on Einstein's notation: a repeated index indicates summation.This generic kinematic field can be used to formulate several displacement-based theories.Thanks to this notation, problem's governing differential equations and boundary conditions can be derived in terms of a single "fundamental nucleus".In this paper, Maclaurin's polynomials are assumed as approximation functions   .  and   as functions of the expansion order  can be obtained via Pascal's triangle as shown in Table 1; see Giunta et al. [13].The actual governing differential equations and boundary conditions due to a fixed approximation order and polynomials type are obtained straightforwardly via summation of the nucleus corresponding to each term of the expansion.According to the previous choice of polynomial function, the generic -order displacement field is Table 1: Maclaurin's polynomials terms via Pascal's triangle [13].
As far as the first-order approximation order is concerned, the kinematic field is Classical models, such as Timoshenko beam theory (TBT), and Euler-Bernoulli beam theory (EBT), are directly derived from the first-order approximation model.Higher-order models yield a more detailed description of the shear mechanics (no shear correction coefficient is required), of the in-and out-of-section deformations, of the coupling of the spatial directions due to Poisson's effect, and of the torsional mechanics than classical models do.EBT theory neglects them all, since it was formulated to describe the bending mechanics.TBT model accounts for constant shear stress and strain components.In the case of classical models, to contrast the phenomenon known in literature as Poisson's locking, reduced Hook's law should be used; i.e., the material stiffness coefficients should be corrected.For more details see Giunta et al. [28] and Carrera and Brischetto [33].

Governing Equations
The strong form of the governing differential equations and the boundary conditions is obtained via the Principle of Virtual Displacements (PVD): where  stands for a virtual variation,  int represents the internal work, and  ine stands for the inertial work.

Virtual Variation of the Strain Energy.
According to the grouping of the stress and strain components in (2) and ( 3), the virtual variation of the strain energy is considered as sum of two contributes: where  represents the volume of the beam.Using the geometrical relations, (5), the material constitutive equations, (8), and the unified hierarchical approximation of the displacements, (12), and integrating by parts, (18) reads where the differential matrix operator D   operates on u  and the differential matrix operator D  Ω operates on   .Equation (19) in a compact vectorial form reads where K  is the differential linear stiffness matrix nucleus, that is, the representative term of the governing equations.
The generic term  ℎ  (,)  (,) is a cross section moment: Double dots stand for a second derivative with respect to time ().Accounting for ( 12), (24) becomes Posing and   being Kronecker's delta, then, the components of the inertial matrix M  are Thus, the compact vectorial form of ( 25) is where either displacement at one beam end is prescribed (the virtual variation is then nil) or the terms between brackets are nil (corresponding to a zero higher-order resultant coherent with a simply supported configuration).

Closed Form Analytical Solution
The resulting boundary value problem is solved via a Naviertype solution.The following displacement field is adopted: where  is represents the half-wave number along the beam axis,  is the angular frequency, and {  :  = , , } are the maximal amplitudes of the displacement components.The displacement field in (31) satisfies boundary conditions (30) since By substituting (31) into linear algebraic system (29), the following generalised eigenvalue problem is obtained: In a compact vectorial form, (29) reduces to where M is the global mass matrix, K is the global stiffness matrix, and q is the amplitude of the displacement field.

Results and Discussion
Beams with a square cross section are considered.The beam is bounded such that with  being the length,  = 1 m the width, and ℎ = 1 m the height of the beam; see Figure 1.A slenderness ratio /ℎ as high as 100 (slender beams) and as low as five (thick beams) is accounted for.ARALL and CARALL FMLs are considered.Three-and five-layer beams are investigated.Table 2 presents the properties of the FRP layers (see Ravishankar et al. [14]).Subscript "" is used to address the composite material.Mechanical properties of aluminium are: Young's modulus   equal to 71.7 GPa, Poisson's ratio ]  equal to 0.33, and density   equal to 2.80 ⋅ 10 3 kg/m 3 (subscript "" represents the metallic material).Modes with a half-wave ( = 1) are studied and the natural frequency is put into the following dimensionless form: In order to validate the models, the results are compared to the three-dimensional finite element solution obtained through the commercial code Ansys.Reference solutions are addressed within the paper as "FEM 3D  " where the subscript "" stands for the number of elements along the yor, equivalently, z-axis (the axial to cross section length ratio being equal to ten).The quadratic twenty-node solid element "SOLID186" is used.In order to show the convergence of the reference three-dimensional FEM solution (up to four significant digits for all the considered results), two meshes are considered: The first has  = 10 and in the second one  = 20.As far as the computational costs are concerned, the degrees of freedom (DOFs) of the three-dimensional FEM solution over a beam cross section as function of  are For  = 20, the DOFs are then 3  843.For a fixed approximation order , the total DOFs of the proposed solutions are In the case of the highest considered expansion order ( = 19), they are 630.
Bending (on  and  plane), torsional, shear (on  and  plane), and axial modes are investigated.For slender beams (/ℎ = 100), the shear modes are much higher than the other considered modes (after all classical models hypotheses are satisfied in this case) and they are not reported.In all the tables, the frequencies are ordered according to the order of appearance in the reference three-dimensional FEM solution: the first mode ("Mode 1") is a bending mode in the  plane, the second one ("Mode 2") corresponds to a bending mode in the  plane, the third one ("Mode 3") is a torsional mode, the forth one ("Mode 4") is an axial mode, the fifth one ("Mode 5") is a shear mode in the  plane, and, finally, the sixth one ("Mode 6") is a shear mode in the  plane.Higher-order beam model respects this order, whereas this is not always the case for classical and some low-order model.The mode swapping can be easily recognised by the fact that for a given model the frequencies are not in an ascending order.It should be noted that for two types of considered FML beams the order of apparition of the different modes is the same.Tables 3 and 4 present the dimensionless frequencies in the case of a three-layer slender beam.Classical theories accurately predict the bending and axial frequencies.A fourth-order model is required to provide good results for the torsional mode, the maximal difference being about 1%, at worst.The case of a length-to-height ratio equal to ten is presented in Tables 5 and 6.For the ARALL beam, a third-order model accurately predicts the bending and the shear frequencies.A fourth-order theory is required for the torsional and axial mode, the error being less than 1% compared with the reference solution.Classical theories yield good results for the bending and axial frequencies of CARALL beams.A fourth-order model is required to provide the torsional and shear frequencies when a difference with the reference threedimensional FEM solution of less than 1% is required.The shear mode in the  plane for the CARALL beam is the only exception.For this one, a seventh-order model is needed.Tables 7 and 8 show the frequencies for thick beams (/ℎ = 5).Classical theories do not accurately predict the different modes.For the ARALL beam, low-order models ( ≤ 4) will accurately predict the bending frequencies in the  plane and the torsional and the shear frequencies in the  plane.Higher-order theories are required for the other modes.A fifth-order theory predicts accurately the bending and shear modes in the  plane, the error being less than 1% compared to the reference solution.For the CARALL beam, a fourthorder approximation is required for the bending and the axial frequencies; the error being less than 1%.For the torsional mode and the shear mode in the  plane, a fifth-order is sufficient for the same accuracy.For the shear mode in the  plane, a seventh-order model is needed.
Figures 2 and 3 present the axial and the shear mode in the  plane of three-layer ARALL FML beams with  length-to-height ratio equal to five.Those for CARALL are very similar and are not reported for the sake of brevity.From these figures, the difference in the material properties of the aluminium and the fibres is clearly visible since the cross section is no longer rigid.
Finally, results for five-layer FML beams with a lengthto-height ratio equal to ten are reported in Tables 9 and  10.Classical theories do not accurately predict the different modes.For an ARALL beam, a third-order theory is needed for the bending modes and the shear mode in the  plane.Fourth-order and second-order models predict accurately the torsional and the axial frequencies, respectively.For the shear mode in the  plane, a seventh-order theory is required for a percentage error less than 1%.For the CARALL beam, classical theories are accurate only for the axial mode, higher-order models being required for the other modes.
For the bending modes, fifth-and sixth-order models are needed.An eighth-order theory yields good results for the torsional mode, whereas shear modes call for a ninth-order theory.

Conclusions
Several models for the free vibration analysis of FML beam structures have been derived via Carrera's Unified Formulation.Via this approach, higher-order theories that  account for shear deformations, in-and out-of-plane warping can be derived straightforwardly since the displacement field is formulated in a general manner regardless of the approximation order over the beam cross section.Classical models, such as Euler-Bernoulli and Timoshenko, are obtained as particular cases.A closed form, Navier-type solution has been used.Natural frequencies associated with bending, torsional, axial, and shear modes have been investigated.Slender and thick, three-and five-layer FML beams have been considered.Three-dimensional FEM solutions obtained via the commercial code Ansys have been considered for validating the proposed models.Classical models are accurate only in the case of bending frequencies of slender beams.As soon as short beams are considered, higher approximation orders are required to describe accurately the free vibrations problem.It has been shown that the proposed formulation
Bending mode in the plane . Bending mode in the plane . Torsional mode. Axial mode. Mode not provided by the theory.
Bending mode in the plane . Bending mode in the plane . Torsional mode. Axial mode. Mode not provided by the theory.

Table 5 :
Three-layer ARALL beam, /ℎ = 10.Bending mode in the plane . Bending mode in the plane . Torsional mode. Axial mode. Shear mode in the plane . Shear mode in the plane . Mode not provided by the theory.
Bending mode in the plane .Bending mode in the plane .Torsional mode.Axial mode.Shear mode in the plane .Shear mode in the plane .Mode not provided by the theory.

Table 8 :
Three-layer CARALL beam, /ℎ = 5. Bending mode in the plane . Bending mode in the plane . Torsional mode. Axial mode. Shear mode in the plane . Shear mode in the plane . Mode not provided by the theory.

Table 9 :
Five-layer ARALL beam, /ℎ = 10.Bending mode in the plane .Bending mode in the plane .Torsional mode. Axial mode. Shear mode in the plane . Shear mode in the plane . Mode not provided by the theory.