Mathematical Problems in Engineering

Mathematical Problems in Engineering / 2018 / Article

Research Article | Open Access

Volume 2018 |Article ID 2724781 |

L. Hanten, G. Giunta, S. Belouettar, V. Salnikov, "Free Vibration Analysis of Fibre-Metal Laminated Beams via Hierarchical One-Dimensional Models", Mathematical Problems in Engineering, vol. 2018, Article ID 2724781, 12 pages, 2018.

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

Academic Editor: Anders Eriksson
Received30 Jan 2018
Accepted31 May 2018
Published26 Jun 2018


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 order , 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.

1. 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 to 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 ) 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 to mm) are often made of aluminium alloys such as , , or . 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 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 [2125] and extended to beams [2631]. 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 () 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 . Cross section geometry and reference system are reported in Figure 1. The displacement field is

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 strain-displacement 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 , , and are the following differential matrix operators:

and is the unit matrix. Under the hypothesis of linear elastic behaviour of the material, the generalised Hooke law holds:

in which is the material stiffness matrix. According to (2) and (3), it reads

Matrices , , , and in (8) are

where terms are the material stiffness coefficients:


and are Young’s moduli, Poisson’s ratios, and the shear moduli, respectively.

3. 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

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].

4. 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, represents the internal work, and 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:

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 operates on and the differential matrix operator operates on .

Equation (19) in a compact vectorial form reads

where is the differential linear stiffness matrix nucleus, that is, the representative term of the governing equations. is the differential linear matrix nucleus of the boundary conditions. The components of the differential linear stiffness matrix are

The generic term is a cross section moment:

As far as the boundary conditions are concerned, the components of are

4.2. Virtual Variation of the Inertial Work

The virtual variation of the inertial work is

Double dots stand for a second derivative with respect to time (). Accounting for (12), (24) becomes


and being Kronecker’s delta, then, the components of the inertial matrix are

Thus, the compact vectorial form of (25) is

4.3. Governing Equations’ and Boundary Conditions’ Fundamental Nuclei

Finally, the explicit form of the fundamental nucleus of the governing equations is

The boundary conditions are

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).

5. Closed Form Analytical Solution

The resulting boundary value problem is solved via a Navier-type 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 is the global mass matrix, is the global stiffness matrix, and 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

with being the length, m the width, and m the height of the beam; see Figure 1. A slenderness ratio as high as (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 GPa, Poisson’s ratio equal to , and density equal to kg/ (subscript “” represents the metallic material). Modes with a half-wave () are studied and the natural frequency is put into the following dimensionless form:

[GPa] [GPa] [GPa] [GPa] [GPa] [GPa] []

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 1600

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 ” where the subscript “” 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 and in the second one . 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 , the DOFs are then 843. For a fixed approximation order , the total DOFs of the proposed solutions are

In the case of the highest considered expansion order (), they are .

Bending (on and plane), torsional, shear (on and plane), and axial modes are investigated. For slender beams (), 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 , 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 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 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 (). Classical theories do not accurately predict the different modes. For the ARALL beam, low-order models () 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 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 . 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.

Mode Mode Mode Mode

FEM 3.5597 6.0478 1.7718 6.6721
FEM 3.5698 6.0650 1.7743 6.6911
3.5700 6.0652 1.7753 6.6913
3.5700 6.0654 1.7758 6.6914
3.5700 6.0654 1.7762 6.6915
3.5701 6.0655 1.7763 6.6915
3.5702 6.0655 1.7763 6.6916
3.5703 6.0659 1.7784 6.6916
3.5703 6.0659 1.7784 6.6920
3.5704 6.0667 1.7802 6.6921
3.5705 6.0670 1.7802 6.6930
3.5705 6.0679 1.9716 6.6935
3.5708 6.0687 1.9716 6.6945
3.5702 6.0662 6.6922
3.5707 6.0690 - 6.6923

Bending mode in the plane .
Bending mode in the plane .
Torsional mode.
Axial mode.
Mode not provided by the theory.

Mode Mode Mode Mode

FEM 3.4870 5.7544 1.9753 6.3546
FEM 3.4871 5.7546 1.9753 6.3547
3.4908 5.7623 1.9781 6.3618
3.4927 5.7631 1.9790 6.3639
3.4946 5.7701 1.9809 6.3711
3.4958 5.7724 1.9810 6.3711
3.4958 5.7724 1.9810 6.3738
3.5026 5.7736 1.9831 6.3738
3.5026 5.7736 1.9831 6.3750
3.5216 5.7905 1.9947 6.3750
3.5216 5.7906 1.9947 6.3947
3.5518 5.8719 2.1225 6.3950
3.5519 5.8724 2.1225 6.4773
3.4771 5.7497 6.3422
3.4775 5.7515 - 6.3422

Bending mode in the plane .
Bending mode in the plane .
Torsional mode.
Axial mode.
Mode not provided by the theory.

Mode Mode Mode Mode Mode Mode

FEM 3.4597 5.6766 1.7724 6.5591 1.8280 2.1088
FEM 3.4696 5.6929 1.7750 6.5762 1.8280 2.1088
3.4698 5.6932 1.7760 6.5770 1.8309 2.1140
3.4700 5.6936 1.7765 6.5776 1.8318 2.1141
3.4702 5.6939 1.7768 6.5778 1.8328 2.1142
3.4702 5.6939 1.7770 6.5782 1.8328 2.1142
3.4713 5.6956 1.7770 6.5783 1.8329 2.1146
3.4714 5.6958 1.7789 6.5813 1.8329 2.1147
3.4744 5.7024 1.7789 6.5816 1.8344 2.1162
3.4744 5.7028 1.7806 6.5939 1.8345 2.1165
3.4800 5.7281 1.7809 6.5943 1.8428 2.1235
3.4800 5.7286 1.9716 6.6467 1.8435 2.1246
3.5019 5.7936 1.9716 6.6474 2.0889 2.3280
3.5003 5.7919 6.6922 2.0910 2.3309
3.5537 6.0441 - 6.6923 - -

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.

Mode Mode Mode Mode Mode Mode

3.4293 5.5176 1.9753 6.2849 2.2640 2.3267
3.4294 5.5179 1.9753 6.2850 2.2640 2.3267
3.4330 5.5248 1.9781 6.2926 2.2721 2.3275
3.4348 5.5256 1.9790 6.2949 2.2750 2.3278
3.4368 5.5322 1.9809 6.3024 2.2802 2.3290
3.4379 5.5343 1.9810 6.3031 2.2802 2.3291
3.4381 5.5343 1.9810 6.3059 2.2806 2.3293
3.4444 5.5353 1.9831 6.3060 2.2808 2.3294
3.4445 5.5358 1.9831 6.3072 2.2928 2.3313
3.4623 5.5512 1.9948 6.3090 2.2937 2.3319
3.4624 5.5581 1.9950 6.3291 2.3314 2.3612
3.4922 5.6316 2.1225 6.3601 2.3318 2.3651
3.4970 5.6706 2.1226 6.4463 2.4103 2.6280
3.4233 5.5598 6.3422 2.4115 2.6248
3.4613 5.7280 - 6.3422 - -

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.

Mode Mode Mode Mode Mode Mode

1.2835 1.9524 3.5486 1.2492 1.9610 2.3446
1.2872 1.9581 3.5539 1.2577 1.9619 2.3445
1.2872 1.9582 3.5556 1.2521 1.9641 2.3506
1.2875 1.9586 3.5566 1.2525 1.9649 2.3513
1.2877 1.9588 3.5574 1.2526 1.9664 2.3516
1.2877 1.9588 3.5576 1.2528 1.9665 2.3518
1.2889 1.9603 3.5577 1.2529 1.9666 2.3532
1.2889 1.9603 3.5609 1.2547 1.9667 2.3541
1.2924 1.9667 3.5610 1.2547 1.9669 2.3603
1.2924 1.9668 3.5636 1.2628 1.9671 2.3625
1.2993 1.9924 3.5656 1.2628 1.9710 2.3923
1.2994 1.9927 3.9430 1.3012 1.9737 2.4005
1.3278 2.0619 3.9432 1.3014 2.1988 2.6006
1.3263 2.0618 1.3360 2.2074 2.6191
1.4014 2.3883 - 1.3385 - -

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.

Mode Mode Mode Mode Mode Mode

1.3095 1.9821 3.9505 1.2167 2.3598 2.5200
1.3097 1.9824 3.9507 1.2167 2.3597 2.5199
1.3109 1.9844 3.9563 1.2185 2.3682 2.5218
1.3115 1.9847 3.9580 1.2190 2.3715 2.5223
1.3124 1.9870 3.9620 1.2207 2.3766 2.5248
1.3128 1.9876 3.9621 1.2213 2.3768 2.5252
1.3130 1.9877 3.9622 1.2219 2.3773 2.5256
1.3151 1.9879 3.9664 1.2220 2.3783 2.5258
1.3152 1.9886 3.9666 1.2222 2.3922 2.5291
1.3211 1.9935 3.9899 1.2234 2.3954 2.5313
1.3212 2.0015 3.9914 1.2276 2.4385 2.5851
1.3328 2.0247 4.2453 1.2509 2.4405 2.6001
1.3394 2.0701 4.2456 1.2707 2.5135 2.8777
1.3111 2.0364 1.2684 2.5186 2.8665
1.3654 2.2635 - 1.2684 - -

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.

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 length-to-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 . 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.

Mode Mode Mode Mode Mode Mode

4.9660 6.1883 1.8199 7.2386 1.9987 2.1544
4.9661 6.1884 1.8200 7.2387 1.9987 2.1544
4.9849 6.2109 1.8265 7.2647 2.0067 2.1619
4.9850 6.2111 1.8265 7.2648 2.0068 2.1619
4.9872 6.2127 1.8268 7.2651 2.0070 2.1624
4.9875 6.2132 1.8294 7.2680 2.0071 2.1627
4.9917 6.2188 1.8294 7.2687 2.0110 2.1640
4.9917 6.2195 1.8342 7.2791 2.0110 2.1644
4.9923 6.2244 1.8342 7.2801 2.0229 2.1654
4.9925 6.2245 1.8373 7.2882 2.0229 2.1656
5.0006 6.2262 1.8375 7.2884 2.0306 2.1670
5.0008 6.2264 2.0073 7.2933 2.0316 2.1683
5.0833 6.2899 2.0073 7.2936 2.2013 2.3693
5.0806 6.2878 7.3153 2.2035 2.3727
5.2443 6.6069 - 7.3154 - -

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.

Mode Mode Mode Mode Mode Mode

4.8378 5.9934 2.0625 6.8880 2.3883 2.4585
4.8380 5.9937 2.0625 6.8881 2.3883 2.4585
4.8466 6.0038 2.0683 6.9005 2.4019 2.4594
4.8498 6.0115 2.0691 6.9101 2.4042 2.4603
4.8529 6.0136 2.0701 6.9130 2.4079 2.4605
4.8571 6.0213 2.0771 6.9131 2.4081 2.4607
4.8582 6.0221 2.0771 6.9213 2.4336 2.4642