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

Authors: L. Hanten; G. Giunta; S. Belouettar; V. Salnikov
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

This paper presents a free vibration analysis of beams made 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 orderN, 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. Three- and 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 three-dimensional elements showing that the proposed approach is accurate yet computationally effective.


---

## Body

## 1. Introduction

In Fibre-Metal Laminates (FMLs), Figure1, 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).Figure 1
FML beam geometry and dimensions.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 solid-like 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–25] and extended to beams [26–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 three-dimensional FEM solutions obtained by commercial code Ansys.


## 2. Preliminaries

In presenting the proposed model, the methodology in Giunta et al. [32] has been followed. A beam is a structure whose axial extension (l) 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: y- and z-axes are two orthogonal directions lying on Ω. The x coordinate is coincident to the axis of the beam and is bounded such that 0≤x≤l. Cross section geometry and reference system are reported in Figure 1. The displacement field is(1)ux,y,z=uxx,y,zuyx,y,zuzx,y,zTin whichux, uy, and uz are the displacement components along x-, y-, and z-axes. Superscript ‘T’ represents the transposition operator. The strain, ε, and stress, σ, tensors in Voigt notation are grouped into two vectors: εn, σn lying on planes orthogonal to Ω where(2)εn=εxxεxyεxzσn=σxxσxyσxzandεp, σp lying on the cross section:(3)εp=εyyεzzεyzσp=σyyσzzσyzIn the case of a linear analysis, the following strain-displacement geometrical relations hold:(4)εn=ux,xux,y+uy,xux,z+uz,xεp=uy,yuz,zuy,z+uz,ySubscripts ‘x’, ‘y’, and ‘z’, when preceded by comma, represent derivation versus the corresponding spatial coordinate. A compact vectorial notation can be adopted for (4):(5)εnT=Dnx+DnΩuεpT=DpuwhereDnΩ, Dnx, and Dp are the following differential matrix operators:(6)DnΩ=000∂∂y00∂∂z00Dnx=I∂∂x,Dp=0∂∂y000∂∂z0∂∂z∂∂y,andI is the unit matrix. Under the hypothesis of linear elastic behaviour of the material, the generalised Hooke law holds:(7)σ=Cεin whichC is the material stiffness matrix. According to (2) and (3), it reads(8)σn=Cnnεn+Cnpεpσp=Cpnεn+CppεpMatricesCnn, Cnp, Cpn, and Cpp in (8) are(9)Cnn=C11000C66000C55Cnp=CpnT=C12C130000000Cpp=C22C230C23C33000C44where termsCij are the material stiffness coefficients:(10)C11=1-ν23ν32E11δ,C12=ν12+ν32ν13E22δC13=ν13+ν12ν23E33δC22=1-ν13ν31E22δC23=ν23+ν21ν13E22δC33=1-ν12ν21E33δC44=G23C55=G31C66=G12where(11)δ=1-ν12ν21-ν23ν32-ν31ν13-2ν21ν32ν13andEij are Young’s moduli, νij Poisson’s ratios, and Gij the shear moduli, respectively.


## 3. Hierarchical Beam Approximation

The displacement field is, a priori, assumed over the cross section in the following way:(12)ux,y,z=Fτy,zuτxwithτ=1,2,…,NuN u represents the number of unknowns while Fτ 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 Fτ. Nu and Fτ as functions of the expansion order N 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 N-order displacement field is(13)ux=ux1+yux2+zux3+⋯+yNuxN2+N+2/2+zNuxN+1N+2/2uy=uy1+yuy2+zuy3+⋯+yNuyN2+N+2/2+zNuyN+1N+2/2uz=uz1+yuz2+zuz3+⋯+yNuzN2+N+2/2+zNuzN+1N+2/2Table 1
Maclaurin’s polynomials terms via Pascal’s triangle [13].
N N u F τ 0 1 F 1 = 1 1 3 F 2 = y   F 3 = z 2 6 F 4 = y 2   F 5 = y z   F 6 = z 2 3 10 F 7 = y 3   F 8 = y 2 z   F 9 = y z 2   F 10 = z 3 … … … N N + 1 N + 2 2 F N 2 + N + 2 / 2 = y N   F N 2 + N + 4 / 2 = y N - 1 z   …   F N N + 3 / 2 = y z N - 1     F N + 1 N + 2 / 2 = z NAs far as the first-order approximation order is concerned, the kinematic field is(14)ux=ux1+ux2y+ux3zuy=uy1+uy2y+uy3zuz=uz1+uz2y+uz3zClassical models, such as Timoshenko beam theory (TBT),(15)ux=ux1+ux2y+ux3zuy=uy1uz=uz1,and Euler-Bernoulli beam theory (EBT),(16)ux=ux1-uy1,xy-uz1,xzuy=uy1uz=uz1,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].


## 4. Governing Equations

The strong form of the governing differential equations and the boundary conditions is obtained via the Principle of Virtual Displacements (PVD):(17)δLint+δLine=0whereδ stands for a virtual variation, Lint represents the internal work, and Line stands for the inertial work.


### 4.1. 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:(18)δLint=∫VδεTσdV=∫VδεnTσn+δεpTσpdVwhereV 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 (19)δLint=∫lδuτT∫ΩDnΩTFτCnnDnx+DnΩFs+CnpDpFsdΩusdx-∫lδuτT∫ΩDnxFτCnnDnx+DnΩFs+CnpDpFsdΩusdx+∫lδuτT∫ΩFτDpTCpnDnx+DnΩFs+CppDpFsdΩusdx+δuτT∫ΩFτCnnDnx+DnΩFs+CnpDpFsdΩusx=0x=lwhere the differential matrix operatorDnxT operates on us and the differential matrix operator DnΩT operates on Fτ.Equation (19) in a compact vectorial form reads(20)δLint=∫lδuτTKτsusdx+δuτTΠτsusx=0x=lwhereKτs is the differential linear stiffness matrix nucleus, that is, the representative term of the governing equations. Πτs is the differential linear matrix nucleus of the boundary conditions. The components of the differential linear stiffness matrix K¯τs are (21)K¯xxτs=Jτ,ys,y66+Jτ,zs,z55-Jτs11∂2∂x2,K¯xyτs=Jτ,ys66-sJτs,y12∂∂x,K¯xzτs=Jτ,zs55-Jτs,z13∂∂xK¯yyτs=Jτ,ys,y22+Jτ,zs,z44-Jτs66∂2∂x2,K¯yxτs=Jτ,ys12-Jτs,y66∂∂x,K¯yzτs=Jτ,ys,z23+Jτ,zs,y44K¯zzτs=Jτ,ys,y44+Jτ,zs,z33-Jτs55∂2∂x2,K¯zxτs=Jτ,zs13-Jτs,z55∂∂x,K¯zyτs=Jτ,zs,y23+Jτ,ys,z44The generic termJτ,ϕs,ξgh is a cross section moment:(22)Jτ,ϕs,ξgh=∫ΩCghFτ,ϕFs,ξdΩAs far as the boundary conditions are concerned, the components ofΠ¯τs are (23)Π¯xxτs=Jτs11∂∂x,Π¯xyτs=Jτs,y12,Π¯xzτs=Jτs,z13Π¯yyτs=Jτs66∂∂x,Π¯yxτs=Jτs,y66,Π¯yzτs=0Π¯zzτs=Jτs55∂∂x,Π¯zxτs=Jτs,z55,Π¯zyτs=0


### 4.2. Virtual Variation of the Inertial Work

The virtual variation of the inertial work is(24)δLine=∫Vρδuu¨dV=∫l∫Ωδρuu¨dΩdxDouble dots stand for a second derivative with respect to time (t). Accounting for (12), (24) becomes(25)δLine=∫lδuτ∫ΩρFτFsdΩu¨sdxPosing(26)Iτs=∫ΩFτρFsdΩandδij being Kronecker’s delta, then, the components of the inertial matrix Mτs are(27)Mijτs=δijIτswithi,j=x,y,zThus, the compact vectorial form of (25) is(28)δLine=∫lδuτMτsu¨sdx


### 4.3. Governing Equations’ and Boundary Conditions’ Fundamental Nuclei

Finally, the explicit form of the fundamental nucleus of the governing equations is(29)δuxτ:-Jτs11uxs,xx+Jτ,ys,y66+Jτ,zs,z55uxs+-Jτs,y12+Jτ,ys66uys,x+-Jτs,z13+Jτ,zs55uzs,x+Iτsu¨xs=0δuyτ:+-Jτs,y66+Jτ,ys12uxs,x-Jτs66uys,xx+Jτ,ys,y22+Jτ,zs,z44uys+Jτ,ys,z23+Jτ,zs,y44uzs+Iτsu¨ys=0δuzτ:-Jτs,z55+Jτ,zs13uxs,x+Jτ,ys,z44+Jτ,zs,y23uys-Jτs55uzs,xx+Jτ,ys,y44+Jτ,zs,z33uzs+Iτsu¨zs=0The boundary conditions are(30)δuxτJτs11uxs,x+Jτs,y12uys+Jτs,z13uzsx=0x=l=0δuyτJτs,y66uxs+Jτs66uys,xx=0x=l=0δuzτJτs,z55uxs+Jτs55uzs,xx=0x=l=0where 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).


## 4.1. 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:(18)δLint=∫VδεTσdV=∫VδεnTσn+δεpTσpdVwhereV 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 (19)δLint=∫lδuτT∫ΩDnΩTFτCnnDnx+DnΩFs+CnpDpFsdΩusdx-∫lδuτT∫ΩDnxFτCnnDnx+DnΩFs+CnpDpFsdΩusdx+∫lδuτT∫ΩFτDpTCpnDnx+DnΩFs+CppDpFsdΩusdx+δuτT∫ΩFτCnnDnx+DnΩFs+CnpDpFsdΩusx=0x=lwhere the differential matrix operatorDnxT operates on us and the differential matrix operator DnΩT operates on Fτ.Equation (19) in a compact vectorial form reads(20)δLint=∫lδuτTKτsusdx+δuτTΠτsusx=0x=lwhereKτs is the differential linear stiffness matrix nucleus, that is, the representative term of the governing equations. Πτs is the differential linear matrix nucleus of the boundary conditions. The components of the differential linear stiffness matrix K¯τs are (21)K¯xxτs=Jτ,ys,y66+Jτ,zs,z55-Jτs11∂2∂x2,K¯xyτs=Jτ,ys66-sJτs,y12∂∂x,K¯xzτs=Jτ,zs55-Jτs,z13∂∂xK¯yyτs=Jτ,ys,y22+Jτ,zs,z44-Jτs66∂2∂x2,K¯yxτs=Jτ,ys12-Jτs,y66∂∂x,K¯yzτs=Jτ,ys,z23+Jτ,zs,y44K¯zzτs=Jτ,ys,y44+Jτ,zs,z33-Jτs55∂2∂x2,K¯zxτs=Jτ,zs13-Jτs,z55∂∂x,K¯zyτs=Jτ,zs,y23+Jτ,ys,z44The generic termJτ,ϕs,ξgh is a cross section moment:(22)Jτ,ϕs,ξgh=∫ΩCghFτ,ϕFs,ξdΩAs far as the boundary conditions are concerned, the components ofΠ¯τs are (23)Π¯xxτs=Jτs11∂∂x,Π¯xyτs=Jτs,y12,Π¯xzτs=Jτs,z13Π¯yyτs=Jτs66∂∂x,Π¯yxτs=Jτs,y66,Π¯yzτs=0Π¯zzτs=Jτs55∂∂x,Π¯zxτs=Jτs,z55,Π¯zyτs=0


## 4.2. Virtual Variation of the Inertial Work

The virtual variation of the inertial work is(24)δLine=∫Vρδuu¨dV=∫l∫Ωδρuu¨dΩdxDouble dots stand for a second derivative with respect to time (t). Accounting for (12), (24) becomes(25)δLine=∫lδuτ∫ΩρFτFsdΩu¨sdxPosing(26)Iτs=∫ΩFτρFsdΩandδij being Kronecker’s delta, then, the components of the inertial matrix Mτs are(27)Mijτs=δijIτswithi,j=x,y,zThus, the compact vectorial form of (25) is(28)δLine=∫lδuτMτsu¨sdx


## 4.3. Governing Equations’ and Boundary Conditions’ Fundamental Nuclei

Finally, the explicit form of the fundamental nucleus of the governing equations is(29)δuxτ:-Jτs11uxs,xx+Jτ,ys,y66+Jτ,zs,z55uxs+-Jτs,y12+Jτ,ys66uys,x+-Jτs,z13+Jτ,zs55uzs,x+Iτsu¨xs=0δuyτ:+-Jτs,y66+Jτ,ys12uxs,x-Jτs66uys,xx+Jτ,ys,y22+Jτ,zs,z44uys+Jτ,ys,z23+Jτ,zs,y44uzs+Iτsu¨ys=0δuzτ:-Jτs,z55+Jτ,zs13uxs,x+Jτ,ys,z44+Jτ,zs,y23uys-Jτs55uzs,xx+Jτ,ys,y44+Jτ,zs,z33uzs+Iτsu¨zs=0The boundary conditions are(30)δuxτJτs11uxs,x+Jτs,y12uys+Jτs,z13uzsx=0x=l=0δuyτJτs,y66uxs+Jτs66uys,xx=0x=l=0δuzτJτs,z55uxs+Jτs55uzs,xx=0x=l=0where 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).


## 5. Closed Form Analytical Solution

The resulting boundary value problem is solved via a Navier-type solution. The following displacement field is adopted:(31)uxsx=Uxscos⁡αxeiωtuysx=Uyssin⁡αxeiωtuzsx=Uzssin⁡αxeiωtwhereα is(32)α=mπlm represents the half-wave number along the beam axis, ω is the angular frequency, and Uis:i=x,y,z are the maximal amplitudes of the displacement components.The displacement field in (31) satisfies boundary conditions (30) since(33)ux,x0=ux,xl=0uy0=uyl=0uz0=uzl=0By substituting (31) into linear algebraic system (29), the following generalised eigenvalue problem is obtained:(34)δUxτ:α2Jτs11+Jτ,ys,y66+Jτ,ys,z56+Jτ,zs,y56+Jτ,zs,z55Uxs+α-Jτs,y12-Jτs,z14+Jτ,ys66+Jτ,zs56Uys+α-Jτs,y14-Jτs,z13+Jτ,ys56+Jτ,zs55Uzs=ω2IτsUxsδUyτ:α+Jτs,y66+Jτs,z56-Jτ,ys12-Jτ,zs14Uxs+α2Jτs66+Jτ,ys,y22+Jτ,ys,z24+Jτ,zs,y24+Jτ,zs,z44Uys+α2Jτs56+Jτ,ys,y24+Jτ,ys,z23+Jτ,zs,y44+Jτ,zs,z34Uzs=ω2IτsUysδUzτ:+Jτs,y56+Jτs,z55-Jτ,ys14-Jτ,zs13Uxs+α2Jτs56+Jτ,ys,y24+Jτ,ys,z44+Jτ,zs,y23+Jτ,zs,z34Uys+α2Jτs55+Jτ,ys,y44+Jτ,ys,z34+Jτ,zs,y34+Jτ,zs,z33Uzs=ω2IτsUzsIn a compact vectorial form, (29) reduces to(35)K-ω2Mq=0whereM is the global mass matrix, K is the global stiffness matrix, and q is the amplitude of the displacement field.


## 6. Results and Discussion

Beams with a square cross section are considered. The beam is bounded such that(36)0≤x≤l,-b2≤y≤b2,-h2≤z≤h2withl being the length, b=1 m the width, and h=1 m the height of the beam; see Figure 1. A slenderness ratio l/h 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 “c” is used to address the composite material. Mechanical properties of aluminium are: Young’s modulus Em equal to 71.7 GPa, Poisson’s ratio νm equal to 0.33, and density ρm equal to 2.80·103 kg/m3 (subscript “m” represents the metallic material). Modes with a half-wave (m=1) are studied and the natural frequency is put into the following dimensionless form:(37)ω¯=l2hρmEmωTable 2
Properties of FRP layers [14].
FRP E c x E c y E c z ν c x y ν c x z ν c y z G c x y G c x z G c y z ρ c [GPa] [GPa] [GPa] [GPa] [GPa] [GPa] [kg/m3] Aramid/Epoxy 67.0 4.7 4.7 0.45 0.34 0.34 2.0 2.0 1.6 1440 Carbon/Epoxy 60.8 58.2 58.2 0.07 0.07 0.40 4.5 4.5 5.0 1600In 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 “FEM3Dn” where the subscript “n” stands for the number of elements along the y- or, 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 n=10 and in the second one n=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 n are (38)DOFsn=33n+1n+1Forn=20, the DOFs are then 3′843. For a fixed approximation order N, the total DOFs of the proposed solutions are(39)DOFsN=3N+1N+22In the case of the highest considered expansion order (N=19), they are 630.Bending (onxz and xy plane), torsional, shear (on xz and xy plane), and axial modes are investigated. For slender beams (l/h=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 xz plane, the second one (“Mode 2”) corresponds to a bending mode in the xy 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 xz plane, and, finally, the sixth one (“Mode 6”) is a shear mode in the xy 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 three-dimensional FEM solution of less than 1% is required. The shear mode in the xz 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 (l/h=5). Classical theories do not accurately predict the different modes. For the ARALL beam, low-order models (N≤4) will accurately predict the bending frequencies in the xz plane and the torsional and the shear frequencies in the xz plane. Higher-order theories are required for the other modes. A fifth-order theory predicts accurately the bending and shear modes in the xy plane, the error being less than 1% compared to the reference solution. For the CARALL beam, a fourth-order 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 xy plane, a fifth-order is sufficient for the same accuracy. For the shear mode in the xz plane, a seventh-order model is needed.Table 3
Three-layer ARALL beam,l/h=100.
ω¯ Mode1a Mode2b Mode3c Mode4d × 1 0 2 × 1 0 2 × 1 × 1 FEM3D20 3.5597 6.0478 1.7718 6.6721 FEM3D10 3.5698 6.0650 1.7743 6.6911 N=19 3.5700 6.0652 1.7753 6.6913 N=14 3.5700 6.0654 1.7758 6.6914 N=10 3.5700 6.0654 1.7762 6.6915 N=9 3.5701 6.0655 1.7763 6.6915 N=8 3.5702 6.0655 1.7763 6.6916 N=7 3.5703 6.0659 1.7784 6.6916 N=6 3.5703 6.0659 1.7784 6.6920 N=5 3.5704 6.0667 1.7802 6.6921 N=4 3.5705 6.0670 1.7802 6.6930 N=3 3.5705 6.0679 1.9716 6.6935 N=2 3.5708 6.0687 1.9716 6.6945 TBT 3.5702 6.0662 - e 6.6922 EBT 3.5707 6.0690 - 6.6923 aBending mode in the plane xz. bBending mode in the plane xy. cTorsional mode. dAxial mode. eMode not provided by the theory.Table 4
Three-layer CARALL beam,l/h=100.
ω¯ Mode1a Mode2b Mode3c Mode4d × 1 0 2 × 1 0 2 × 1 × 1 FEM3D20 3.4870 5.7544 1.9753 6.3546 FEM3D10 3.4871 5.7546 1.9753 6.3547 N=19 3.4908 5.7623 1.9781 6.3618 N=14 3.4927 5.7631 1.9790 6.3639 N=10 3.4946 5.7701 1.9809 6.3711 N=9 3.4958 5.7724 1.9810 6.3711 N=8 3.4958 5.7724 1.9810 6.3738 N=7 3.5026 5.7736 1.9831 6.3738 N=6 3.5026 5.7736 1.9831 6.3750 N=5 3.5216 5.7905 1.9947 6.3750 N=4 3.5216 5.7906 1.9947 6.3947 N=3 3.5518 5.8719 2.1225 6.3950 N=2 3.5519 5.8724 2.1225 6.4773 TBT 3.4771 5.7497 - e 6.3422 EBT 3.4775 5.7515 - 6.3422 aBending mode in the plane xz. bBending mode in the plane xy. cTorsional mode. dAxial mode. eMode not provided by the theory.Table 5
Three-layer ARALL beam,l/h=10.
ω¯ Mode1a Mode2b Mode3c Mode4d Mode5e Mode6f × 1 × 1 × 1 0 - 1 × 1 0 - 1 × 1 0 - 2 × 1 0 - 2 FEM3D20 3.4597 5.6766 1.7724 6.5591 1.8280 2.1088 FEM3D10 3.4696 5.6929 1.7750 6.5762 1.8280 2.1088 N=19 3.4698 5.6932 1.7760 6.5770 1.8309 2.1140 N=14 3.4700 5.6936 1.7765 6.5776 1.8318 2.1141 N=10 3.4702 5.6939 1.7768 6.5778 1.8328 2.1142 N=9 3.4702 5.6939 1.7770 6.5782 1.8328 2.1142 N=8 3.4713 5.6956 1.7770 6.5783 1.8329 2.1146 N=7 3.4714 5.6958 1.7789 6.5813 1.8329 2.1147 N=6 3.4744 5.7024 1.7789 6.5816 1.8344 2.1162 N=5 3.4744 5.7028 1.7806 6.5939 1.8345 2.1165 N=4 3.4800 5.7281 1.7809 6.5943 1.8428 2.1235 N=3 3.4800 5.7286 1.9716 6.6467 1.8435 2.1246 N=2 3.5019 5.7936 1.9716 6.6474 2.0889 2.3280 TBT 3.5003 5.7919 - g 6.6922 2.0910 2.3309 EBT 3.5537 6.0441 - 6.6923 - - aBending mode in the plane xz. bBending mode in the plane xy. cTorsional mode. dAxial mode. eShear mode in the plane xz. fShear mode in the plane xy. gMode not provided by the theory.Table 6
Three-layer CARALL beam,l/h=10.
ω¯ Mode1a Mode2b Mode3c Mode4d Mode5e Mode6f × 1 × 1 × 1 0 - 1 × 1 0 - 1 × 1 0 - 2 × 1 0 - 2 FEM3D20 3.4293 5.5176 1.9753 6.2849 2.2640 2.3267 FEM3D10 3.4294 5.5179 1.9753 6.2850 2.2640 2.3267 N=19 3.4330 5.5248 1.9781 6.2926 2.2721 2.3275 N=14 3.4348 5.5256 1.9790 6.2949 2.2750 2.3278 N=10 3.4368 5.5322 1.9809 6.3024 2.2802 2.3290 N=9 3.4379 5.5343 1.9810 6.3031 2.2802 2.3291 N=8 3.4381 5.5343 1.9810 6.3059 2.2806 2.3293 N=7 3.4444 5.5353 1.9831 6.3060 2.2808 2.3294 N=6 3.4445 5.5358 1.9831 6.3072 2.2928 2.3313 N=5 3.4623 5.5512 1.9948 6.3090 2.2937 2.3319 N=4 3.4624 5.5581 1.9950 6.3291 2.3314 2.3612 N=3 3.4922 5.6316 2.1225 6.3601 2.3318 2.3651 N=2 3.4970 5.6706 2.1226 6.4463 2.4103 2.6280 TBT 3.4233 5.5598 - g 6.3422 2.4115 2.6248 EBT 3.4613 5.7280 - 6.3422 - - aBending mode in the plane xz. bBending mode in the plane xy. cTorsional mode. dAxial mode. eShear mode in the plane xz. fShear mode in the plane xy. gMode not provided by the theory.Table 7
Three-layer ARALL beam,l/h=5.
ω¯ Mode1a Mode2b Mode3c Mode4d Mode5e Mode6f × 1 0 - 1 × 1 0 - 1 × 1 0 - 1 × 1 0 - 2 × 1 0 - 2 × 1 0 - 2 FEM3D20 1.2835 1.9524 3.5486 1.2492 1.9610 2.3446 FEM3D10 1.2872 1.9581 3.5539 1.2577 1.9619 2.3445 N=19 1.2872 1.9582 3.5556 1.2521 1.9641 2.3506 N=14 1.2875 1.9586 3.5566 1.2525 1.9649 2.3513 N=10 1.2877 1.9588 3.5574 1.2526 1.9664 2.3516 N=9 1.2877 1.9588 3.5576 1.2528 1.9665 2.3518 N=8 1.2889 1.9603 3.5577 1.2529 1.9666 2.3532 N=7 1.2889 1.9603 3.5609 1.2547 1.9667 2.3541 N=6 1.2924 1.9667 3.5610 1.2547 1.9669 2.3603 N=5 1.2924 1.9668 3.5636 1.2628 1.9671 2.3625 N=4 1.2993 1.9924 3.5656 1.2628 1.9710 2.3923 N=3 1.2994 1.9927 3.9430 1.3012 1.9737 2.4005 N=2 1.3278 2.0619 3.9432 1.3014 2.1988 2.6006 TBT 1.3263 2.0618 - g 1.3360 2.2074 2.6191 EBT 1.4014 2.3883 - 1.3385 - - aBending mode in the plane xz. bBending mode in the plane xy. cTorsional mode. dAxial mode. eShear mode in the plane xz. fShear mode in the plane xy. gMode not provided by the theory.Table 8
Three-layer CARALL beam,l/h=5.
ω¯ Mode1a Mode2b Mode3c Mode4d Mode5e Mode6f × 1 0 - 1 × 1 0 - 1 × 1 0 - 1 × 1 0 - 2 × 1 0 - 2 × 1 0 - 2 FEM3D20 1.3095 1.9821 3.9505 1.2167 2.3598 2.5200 FEM3D10 1.3097 1.9824 3.9507 1.2167 2.3597 2.5199 N=19 1.3109 1.9844 3.9563 1.2185 2.3682 2.5218 N=14 1.3115 1.9847 3.9580 1.2190 2.3715 2.5223 N=10 1.3124 1.9870 3.9620 1.2207 2.3766 2.5248 N=9 1.3128 1.9876 3.9621 1.2213 2.3768 2.5252 N=8 1.3130 1.9877 3.9622 1.2219 2.3773 2.5256 N=7 1.3151 1.9879 3.9664 1.2220 2.3783 2.5258 N=6 1.3152 1.9886 3.9666 1.2222 2.3922 2.5291 N=5 1.3211 1.9935 3.9899 1.2234 2.3954 2.5313 N=4 1.3212 2.0015 3.9914 1.2276 2.4385 2.5851 N=3 1.3328 2.0247 4.2453 1.2509 2.4405 2.6001 N=2 1.3394 2.0701 4.2456 1.2707 2.5135 2.8777 TBT 1.3111 2.0364 - g 1.2684 2.5186 2.8665 EBT 1.3654 2.2635 - 1.2684 - - aBending mode in the plane xz. bBending mode in the plane xy. cTorsional mode. dAxial mode. eShear mode in the plane xz. fShear mode in the plane xy. gMode not provided by the theory.Figures2 and 3 present the axial and the shear mode in the xy 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.Figure 2
Free vibration axial mode for a three-layer ARALL beam,l/h=5.Figure 3
Free vibration shear mode in thexy plane for a three-layer ARALL beam, l/h=5.Finally, results for five-layer FML beams with a length-to-height ratio equal to ten are reported in Tables9 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 xy plane. Fourth-order and second-order models predict accurately the torsional and the axial frequencies, respectively. For the shear mode in the xz 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.Table 9
Five-layer ARALL beam,l/h=10.
ω¯ Mode1a Mode2b Mode3c Mode4d Mode5e Mode6f × 1 × 1 × 1 0 - 1 × 1 0 - 1 × 1 0 - 2 × 1 0 - 2 FEM3D20 4.9660 6.1883 1.8199 7.2386 1.9987 2.1544 FEM3D10 4.9661 6.1884 1.8200 7.2387 1.9987 2.1544 N=19 4.9849 6.2109 1.8265 7.2647 2.0067 2.1619 N=14 4.9850 6.2111 1.8265 7.2648 2.0068 2.1619 N=10 4.9872 6.2127 1.8268 7.2651 2.0070 2.1624 N=9 4.9875 6.2132 1.8294 7.2680 2.0071 2.1627 N=8 4.9917 6.2188 1.8294 7.2687 2.0110 2.1640 N=7 4.9917 6.2195 1.8342 7.2791 2.0110 2.1644 N=6 4.9923 6.2244 1.8342 7.2801 2.0229 2.1654 N=5 4.9925 6.2245 1.8373 7.2882 2.0229 2.1656 N=4 5.0006 6.2262 1.8375 7.2884 2.0306 2.1670 N=3 5.0008 6.2264 2.0073 7.2933 2.0316 2.1683 N=2 5.0833 6.2899 2.0073 7.2936 2.2013 2.3693 TBT 5.0806 6.2878 - g 7.3153 2.2035 2.3727 EBT 5.2443 6.6069 - 7.3154 - - aBending mode in the plane xz. bBending mode in the plane xy. cTorsional mode. dAxial mode. eShear mode in the plane xz. fShear mode in the plane xy. gMode not provided by the theory.Table 10
Five-layer CARALL beam,l/h=10.
ω¯ Mode1a Mode2b Mode3c Mode4d Mode5e Mode6f × 1 × 1 × 1 0 - 1 × 1 0 - 1 × 1 0 - 2 × 1 0 - 2 FEM3D20 4.8378 5.9934 2.0625 6.8880 2.3883 2.4585 FEM3D10 4.8380 5.9937 2.0625 6.8881 2.3883 2.4585 N=19 4.8466 6.0038 2.0683 6.9005 2.4019 2.4594 N=14 4.8498 6.0115 2.0691 6.9101 2.4042 2.4603 N=10 4.8529 6.0136 2.0701 6.9130 2.4079 2.4605 N=9 4.8571 6.0213 2.0771 6.9131 2.4081 2.4607 N=8 4.8582 6.0221 2.0771 6.9213 2.4336 2.4642 N=7 4.8675 6.0470 2.0912 6.9246 2.4340 2.4651 N=6 4.8718 6.0489 2.0912 6.9535 2.4658 2.4739 N=5 4.8784 6.0603 2.0950 6.9619 2.4662 2.4744 N=4 4.8905 6.0604 2.0952 6.9750 2.4668 2.4764 N=3 4.9722 6.0871 2.2658 6.9757 2.4687 2.4780 N=2 5.0105 6.1277 2.2659 7.0072 2.5805 2.7291 TBT 4.8775 6.0271 - g 6.8956 2.5783 2.7260 EBT 4.9813 6.2279 - 6.8957 - - aBending mode in the plane xz. bBending mode in the plane xy. cTorsional mode. dAxial mode. eShear mode in the plane xz. fShear mode in the plane xy. gMode not provided by the theory.


## 7. 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 allows obtaining results as accurate as desired through an appropriate choice of the approximation order without presenting all the complexities that higher-order models present in their formulation.



---

*Source: 2724781-2018-06-26.xml*