# Nonlinear Vibrations of 3D Laminated Composite Beams

Authors: S. Stoykov; S. Margenov
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

A model for 3D laminated composite beams, that is, beams that can vibrate in space and experience longitudinal and torsional deformations, is derived. The model is based on Timoshenko’s theory for bending and assumes that, under torsion, the cross section rotates as a rigid body but can deform longitudinally due to warping. The warping function, which is essential for correct torsional deformations, is computed preliminarily by the finite element method. Geometrical nonlinearity is taken into account by considering Green’s strain tensor. The equation of motion is derived by the principle of virtual work and discretized by thep-version finite element method. The laminates are assumed to be of orthotropic materials. The influence of the angle of orientation of the laminates on the natural frequencies and on the nonlinear modes of vibration is presented. It is shown that, due to asymmetric laminates, there exist bending-longitudinal and bending-torsional coupling in linear analysis. Dynamic responses in time domain are presented and couplings between transverse displacements and torsion are investigated.


---

## Body

## 1. Introduction

The use of composite materials in the engineering applications has been significantly increased in the last decades. Composite materials can be adjusted to meet the design requirements of structures, such as stiffness and strength, by varying the orientations of the layers. The ability to change the dynamical properties of structures is one of the most important advantages of composite materials over ordinary materials.Laminated composite beams have a variety of applications in industry. Beam models can be used to model, for example, helicopter blades, wind turbine blades, aircraft wings, and so forth. Often these structures vibrate in space due to a variety of aerodynamic forces. Furthermore, due to asymmetric position of the laminates, bending-torsional and bending-longitudinal coupling appears, even in linear models. Therefore, accurate models of 3D laminated composite beams are essential for engineers and researchers.A number of analytical and numerical methods for free vibration problems of laminated beams have been proposed in the literature, for example, [1–6]. Goyal and Kapania [7] studied the response of asymmetrically laminated composite beams including torsional and warping effects, as well as shear deformation. Jun et al. [8] used dynamic finite element method to investigate the natural frequencies and mode shapes of generally laminated composite beams. Banerjee [9, 10] derived exact expressions for the frequency equation and mode shapes of composite Timoshenko beams. Coupling between bending and torsional modes, shear deformation, and rotary inertia were taken into account. Machado and Cortínez [11] developed a geometrically nonlinear model for thin-walled composite beams and investigated the postbuckling behaviour of simply supported beams. Sapountzakis [12, 13] used the boundary element method to compute the warping function of composite bars. Damages of the composite beams, such as delamination and interlayer slips, were modelled numerically in [14, 15].Improvements of the classical Bernoulli-Euler and Timoshenko beam theories have been developed in order to better approximate the cross-sectional behaviour due to axial-bending coupling of the different layers. Such improvements are known as equivalent single layer (ESL) and discrete layer theories (DLT) and, as a result, the unknowns in the equation of motion are increased [16]. In ESL theories the number of unknowns is independent of the number of layers, but they lead to poor approximation of the shear stresses. In DLT the number of unknowns is assumed layer by layer, which leads to sufficiently accurate results. The DLT are more accurate than the ESL theories but are computationally expensive. An important subclass of DSL theories are the zigzag theories, which assume a zigzag distribution of the longitudinal displacement, and the number of kinematic variables is independent of the number of layers [17, 18]. The usage of zigzag theories is essential for beams with considerably different material properties of the layers.In this work, 3D beam model for laminated composite beams is derived considering geometrical type of nonlinearity. The model is based on Timoshenko’s theory for bending and assumes that, under torsion, the cross section rotates as a rigid body and deforms longitudinally due to warping [19, 20]. The warping function is obtained preliminarily by solving a partial differential equation with Neumann boundary conditions. The equation of motion of the beam is derived by the principle of virtual work and discretized by the p-version finite element method. The model is validated by comparing the bending natural frequencies with available results from the literature, considering different orientations of the laminates and different boundary conditions. Further to validate the 3D beam model, a comparison with three-dimensional finite elements is presented. It is shown that accurate results can be obtained by the proposed model with few degrees of freedom. The nonlinear part of the model is validated by comparing the static displacement, of the nonlinear model, with available results. Nonlinear modes of vibration, for three different orientations of the laminates, are compared. Finally, the dynamic steady-state responses of the nonlinear model are presented.


## 2. Mathematical Model

A laminated composite beam is considered (Figure1); the layers are assumed to be of orthotropic materials. The principal material axes can be oriented in an arbitrary angle with respect to the longitudinal coordinate x. The length, the width, and the height of the beam are denoted by l, b, and h, respectively.Figure 1
Geometry of laminated composite beam.


### 2.1. Displacement Field

The displacement field is based on Timoshenko’s theory for bending [21] and it assumes that the cross section rotates as a rigid body in its own plane and deforms longitudinally due to warping [22]. This displacement field is suitable for capturing the bending-bending-torsional coupling due to geometrical nonlinearity [19] or due to asymmetric cross sections [20]. Zigzag functions are not included in the current model, because they are meaningless when the laminates are with similar material properties, particularly when the only difference between the laminates is the angle of orientation [23, 24]. The longitudinal u(x,y,z,t) and the transverse displacements v(x,y,z,t) and w(x,y,z,t), which are functions of space coordinates and time, are expressed by functions of the longitudinal coordinate x in the following way:
(1)u(x,y,z,t)=u0(x,t)-(y-yc)ϕz(x,t)iiiiiiiiiiiiiiiiiiii+(z-zc)ϕy(x,t)+ψ(y,z)∂θx∂x(x,t),v(x,y,z,t)=v0(x,t)+ycos⁡(θx(x,t))iiiiiiiiiiiiiiiiiiii-y-zsin(θx(x,t)),w(x,y,z,t)=w0(x,t)+ysin(θx(x,t))iiiiiiiiiiiiiiiiiiiii+zcos⁡(θx(x,t))-z,
where the subscript “0” represents the reference line, which is the line of the twist centres of the cross section, (yc,zc) is the position of the centroid of the cross section, θx is the rotation of the cross section about the longitudinal axis x, ϕy and ϕz denote rotations of the cross section about y and z axes, respectively, and ψ(y,z) is the warping function.The functions on the reference lineu0(x,t), v0(x,t), w0(x,t), θx(x,t), ϕy(x,t), and ϕz(x,t) are expressed by shape functions and generalized coordinates in a local coordinate system.


### 2.2. Warping Function

The warping function, which not only defines the longitudinal deformation of the beam due to torsion, but also defines the torsional rigidity of the beam, is computed numerically by the finite element method. The warping function is defined by a partial differential equation with Neumann boundary conditions. The equation is derived from the equilibrium equations of elasticity [22], where the stresses due to torsion are computed from the displacement field (1). It is assumed that there are no external forces acting on the lateral surface of the beam in the axial direction:
(2)C~66k∂2ψ∂y2+C~55k∂2ψ∂z2=0inΩk,(C~66k-C~66k-1)(∂ψ∂y-z)ny+(C~55k-C~55k-1)(∂ψ∂z+y)nz=0iiiiiiiiiiiiiiiiiiiiiiiiion⁡Γk,k=1,…,M,
where Ωk is the cross section of the kth layer and Γk is its contour and M is the number of layers. C~55k and C~66k are the shear modulus of the material of the kth layer from the reduced stiffness matrix. ny and nz are the components of the normal vector to the boundary; it points in the outer direction for the outer boundaries and it is considered to point in the lower numbered layer (Figure 1) for the boundaries between different layers. It is considered that C~66k-1=C~55k-1=0, for the outer boundaries.The boundary conditions between different layers are derived by employing that the traction vectors on the interfaces, separating the layers, are equal in magnitude and opposite in direction.Equation (2) is written in weak form and then it is solved by the finite element method:
(3)∫Ωk(C~66k∂υj∂y∂ψ∂y+C~55k∂υj∂z∂ψ∂z)dΩk=∫Γkυj((C~66k-C~66k-1)zny-(C~55k-C~55k-1)ynz)dΓk,
where υj, j=1,…,n, are the weighting functions.


### 2.3. Constitutive Equations

The laminates are assumed to be of orthotropic materials. The stress-strain relations of orthotropic material, expressed in the principal coordinates of the material, are given by(4){σxσyσzτyzτxzτxy}=[C11C12C13000C12C22C23000C13C23C33000000C44000000C55000000C55]{εxεyεzγyzγxzγxy}.
The stress-strain relations of laminate k, rotated on angle α about the z axis, are expressed as
(5){σxσyσzτyzτxzτxy}k=[C-11kC-12kC-13k00C-16kC-12kC-22kC-23k00C-26kC-13kC-23kC-33k00C-36k000C-44kC-45k0000C-45kC-55k0C-16kC-26kC-36k00C-66k]{εxεyεzγyzγxzγxy}k,
where the coefficients of the stiffness matrix of (5) are given to be [23]:
(6)C-11k=C11cos⁡⁡(α)4+2C12sin⁡(α)2cos⁡⁡(α)2iiiiiiiii+4C66sin⁡(α)2cos⁡⁡(α)2+C22sin⁡(α)4C-12k=(C11+C22-4C66)cos⁡⁡(α)2sin⁡(α)2iiiiiiiii+C12(sin⁡(α)4+cos⁡⁡(α)4)C-13k=C13cos⁡⁡(α)2+C23sin⁡(α)2C-16k=(C11-C12-2C66)cos⁡⁡(α)3sin⁡(α)iiiiiiiii+(C12-C22+2C66)sin⁡(α)3cos⁡⁡(α)C-22k=C11sin⁡(α)4+C22cos⁡⁡(α)4iiiiiiiii+2(C12+2C66)sin⁡(α)2cos⁡⁡(α)2C-23k=C13sin⁡(α)2+C23cos⁡⁡(α)2C-26k=(C11-C12-2C66)cos⁡⁡(α)sin⁡(α)3iiiiiiiii+(C12-C22+2C66)sin⁡(α)cos⁡⁡(α)3C-33k=C33C-36k=(C13-C23)cos⁡⁡(α)sin⁡(α)C-44k=C44cos⁡⁡(α)2+C55sin⁡(α)2C-45k=(C55-C44)cos⁡⁡(α)sin⁡(α)C-55k=C44sin⁡(α)2+C55cos⁡⁡(α)2C-66k=(C11-2C12+C22-2C66)cos⁡⁡(α)2sin⁡(α)2iiiiiiiii+C66(sin⁡(α)4+cos⁡⁡(α)4).
Taking into account that the beam is slender, the following assumption for the stresses can be introduced:
(7)σy=σz=τyz=0.
Considering assumption (7) and the relation (5), the reduced stress-strain relation for beam is obtained:
(8){σxτxzτxy}k=[C~11k0C~16k0C~55k0C~16k0C~66k]{εxγxzγxy}k.
The coefficients of the reduced stiffness matrix are given in the appendix.


### 2.4. Strain Expressions

Geometrically nonlinear deformation is considered and the axial and shear strains are derived from Green’s strain tensor [25], where the longitudinal terms of second order are neglected:
(9)iεx=∂u∂x+12(∂v∂x)2+12(∂w∂x)2,γxz=∂u∂z+∂w∂x+∂v∂x∂v∂z+∂w∂x∂w∂z,γxy=∂u∂y+∂v∂x+∂v∂x∂v∂y+∂w∂x∂w∂y.
A square root of the shear correction λ factor is applied solely to the bending terms of the shear strains. This formulation is preferred here, instead of applying the shear correction factor in Hooke’s law, because of the torsional terms which appear in the shear strains, for which a shear correction factor is not necessary. With this formulation, the expressions from the virtual work of internal forces related to torsion will not contain the shear correction factor, while the ones related to bending will be multiplied by the shear correction factor [19]. In the numerical results, a shear correction factor equal to 5/6 is used. Substituting the displacement field (1) into the expressions for the strains (9) and assuming the shear correction factor λ, the strains become
(10)εx=∂u0∂x-(y-yc)∂ϕz∂x+(z-zc)∂ϕy∂xiiiiiii+ψ∂2θx∂x2+12(∂v0∂x)2+12z2(∂θx∂x)2+12(∂w0∂x)2+12y2(∂θx∂x)2-z∂v0∂x∂θx∂x+y∂w0∂x∂θx∂x,γxz=(∂ψ∂z+y)∂θx∂x+λ∂w0∂x+λϕy-λ∂v0∂xθx,γxy=(∂ψ∂y-z)∂θx∂x+λ∂v0∂x-λϕz+λ∂w0∂xθx.


### 2.5. Equation of Motion

The equation of motion is derived by the principle of virtual work and it is discretized into a system of ordinary differential equation by thep-version finite element method. The p-FEM has several advantages over the h-FEM; the most important is that it requires far fewer degrees of freedom than the h-FEM [26]. One element is used for the laminated composite beam and improvement of the accuracy is achieved by increasing the number of the shape functions. Thus, the functions on the reference line, u0(x,t), v0(x,t), w0(x,t), θx(x,t), ϕy(x,t), and ϕz(x,t), are expressed by shape functions and generalized coordinates in a local coordinate system:
(11)d0={u0(ξ,t)v0(ξ,t)w0(ξ,t)θx(ξ,t)ϕy(ξ,t)ϕz(ξ,t)}=N(ξ)q(t),
where q(t) is the vector of generalized coordinates:
(12)q(t)={qu(t)qv(t)qw(t)qθx(t)qϕy(t)qϕz(t)}
and N(ξ) is the matrix of shape functions:(13)N(ξ)=[Nu(ξ)T000000Nv(ξ)T000000Nw(ξ)T000000Nθx(ξ)T000000Nϕy(ξ)T000000Nϕz(ξ)T].Nu(ξ)T, Nv(ξ)T, Nw(ξ)T, Nθx(ξ)T, Nϕy(ξ)T, and Nϕz(ξ)T are, in this order, the row vectors of longitudinal, transverse along y, transverse along z, torsional, rotational about y, and rotational about z shape functions and ξ∈[-1⋯1] is the local coordinate, ξ=2x/l.The shape functions have to satisfy the geometric boundary conditions. The sets of shape functions, used in previous works for isotropic beams, are implemented here [19, 26]. They are suitable for beams with clamped-clamped boundary conditions. If one wants to apply simply supported or free boundary conditions, additional shape functions, such as Hermite cubics, have to be included.The equation of motion is derived by the principle of virtual work:(14)δWV+δWin+δWE=0,
where WV, Win, and WE are the virtual work of internal, inertia, and external forces due to a virtual displacement δd:
(15)δd={δuδvδw}
and d is the vector of displacement components. The variations of the work of internal, inertia, and external forces are defined as
(16)δWV=-∫VδεTσdV,δWin=-∫VρδdTd¨dV,δWE=∫Vδd0TF0dV,
where δε represents the virtual strains, d¨ is the acceleration of a point of the beam, ρ is the density of the beam, and F0 represents external forces applied on the reference line.Replacing (16) into (14) and introducing mass proportional and frequency dependent damping, the following system of second order ordinary differential equations is obtained:
(17)Mq¨(t)+βωMq˙(t)+K(q(t))q(t)=Fcos⁡⁡(ωt),
where M is the mass matrix obtained from the inertia forces and K(q(t)) is the stiffness matrix obtained from the internal forces; it depends on the vector of generalized coordinates q(t), but it also has constant terms. F represents vector of generalized external forces; a harmonic excitation with frequency ω is applied. β is damping factor, which can be estimated from experiments. In the following work, β=0.01ωl12 is used, where ωl1 is the fundamental frequency. The equation of motion (17) is solved in time domain by Newmark’s method [27] and the resulting nonlinear algebraic system is solved by Newton-Raphson’s method [27], by using numerical computation of the Jacobian. Results of nonlinear modes of vibration are also presented; for that purpose the vector of generalized coordinates is expressed in Fourier series and (17) is transformed into a system of nonlinear algebraic equations, which is solved by the arc-length continuation method.


## 2.1. Displacement Field

The displacement field is based on Timoshenko’s theory for bending [21] and it assumes that the cross section rotates as a rigid body in its own plane and deforms longitudinally due to warping [22]. This displacement field is suitable for capturing the bending-bending-torsional coupling due to geometrical nonlinearity [19] or due to asymmetric cross sections [20]. Zigzag functions are not included in the current model, because they are meaningless when the laminates are with similar material properties, particularly when the only difference between the laminates is the angle of orientation [23, 24]. The longitudinal u(x,y,z,t) and the transverse displacements v(x,y,z,t) and w(x,y,z,t), which are functions of space coordinates and time, are expressed by functions of the longitudinal coordinate x in the following way:
(1)u(x,y,z,t)=u0(x,t)-(y-yc)ϕz(x,t)iiiiiiiiiiiiiiiiiiii+(z-zc)ϕy(x,t)+ψ(y,z)∂θx∂x(x,t),v(x,y,z,t)=v0(x,t)+ycos⁡(θx(x,t))iiiiiiiiiiiiiiiiiiii-y-zsin(θx(x,t)),w(x,y,z,t)=w0(x,t)+ysin(θx(x,t))iiiiiiiiiiiiiiiiiiiii+zcos⁡(θx(x,t))-z,
where the subscript “0” represents the reference line, which is the line of the twist centres of the cross section, (yc,zc) is the position of the centroid of the cross section, θx is the rotation of the cross section about the longitudinal axis x, ϕy and ϕz denote rotations of the cross section about y and z axes, respectively, and ψ(y,z) is the warping function.The functions on the reference lineu0(x,t), v0(x,t), w0(x,t), θx(x,t), ϕy(x,t), and ϕz(x,t) are expressed by shape functions and generalized coordinates in a local coordinate system.


## 2.2. Warping Function

The warping function, which not only defines the longitudinal deformation of the beam due to torsion, but also defines the torsional rigidity of the beam, is computed numerically by the finite element method. The warping function is defined by a partial differential equation with Neumann boundary conditions. The equation is derived from the equilibrium equations of elasticity [22], where the stresses due to torsion are computed from the displacement field (1). It is assumed that there are no external forces acting on the lateral surface of the beam in the axial direction:
(2)C~66k∂2ψ∂y2+C~55k∂2ψ∂z2=0inΩk,(C~66k-C~66k-1)(∂ψ∂y-z)ny+(C~55k-C~55k-1)(∂ψ∂z+y)nz=0iiiiiiiiiiiiiiiiiiiiiiiiion⁡Γk,k=1,…,M,
where Ωk is the cross section of the kth layer and Γk is its contour and M is the number of layers. C~55k and C~66k are the shear modulus of the material of the kth layer from the reduced stiffness matrix. ny and nz are the components of the normal vector to the boundary; it points in the outer direction for the outer boundaries and it is considered to point in the lower numbered layer (Figure 1) for the boundaries between different layers. It is considered that C~66k-1=C~55k-1=0, for the outer boundaries.The boundary conditions between different layers are derived by employing that the traction vectors on the interfaces, separating the layers, are equal in magnitude and opposite in direction.Equation (2) is written in weak form and then it is solved by the finite element method:
(3)∫Ωk(C~66k∂υj∂y∂ψ∂y+C~55k∂υj∂z∂ψ∂z)dΩk=∫Γkυj((C~66k-C~66k-1)zny-(C~55k-C~55k-1)ynz)dΓk,
where υj, j=1,…,n, are the weighting functions.


## 2.3. Constitutive Equations

The laminates are assumed to be of orthotropic materials. The stress-strain relations of orthotropic material, expressed in the principal coordinates of the material, are given by(4){σxσyσzτyzτxzτxy}=[C11C12C13000C12C22C23000C13C23C33000000C44000000C55000000C55]{εxεyεzγyzγxzγxy}.
The stress-strain relations of laminate k, rotated on angle α about the z axis, are expressed as
(5){σxσyσzτyzτxzτxy}k=[C-11kC-12kC-13k00C-16kC-12kC-22kC-23k00C-26kC-13kC-23kC-33k00C-36k000C-44kC-45k0000C-45kC-55k0C-16kC-26kC-36k00C-66k]{εxεyεzγyzγxzγxy}k,
where the coefficients of the stiffness matrix of (5) are given to be [23]:
(6)C-11k=C11cos⁡⁡(α)4+2C12sin⁡(α)2cos⁡⁡(α)2iiiiiiiii+4C66sin⁡(α)2cos⁡⁡(α)2+C22sin⁡(α)4C-12k=(C11+C22-4C66)cos⁡⁡(α)2sin⁡(α)2iiiiiiiii+C12(sin⁡(α)4+cos⁡⁡(α)4)C-13k=C13cos⁡⁡(α)2+C23sin⁡(α)2C-16k=(C11-C12-2C66)cos⁡⁡(α)3sin⁡(α)iiiiiiiii+(C12-C22+2C66)sin⁡(α)3cos⁡⁡(α)C-22k=C11sin⁡(α)4+C22cos⁡⁡(α)4iiiiiiiii+2(C12+2C66)sin⁡(α)2cos⁡⁡(α)2C-23k=C13sin⁡(α)2+C23cos⁡⁡(α)2C-26k=(C11-C12-2C66)cos⁡⁡(α)sin⁡(α)3iiiiiiiii+(C12-C22+2C66)sin⁡(α)cos⁡⁡(α)3C-33k=C33C-36k=(C13-C23)cos⁡⁡(α)sin⁡(α)C-44k=C44cos⁡⁡(α)2+C55sin⁡(α)2C-45k=(C55-C44)cos⁡⁡(α)sin⁡(α)C-55k=C44sin⁡(α)2+C55cos⁡⁡(α)2C-66k=(C11-2C12+C22-2C66)cos⁡⁡(α)2sin⁡(α)2iiiiiiiii+C66(sin⁡(α)4+cos⁡⁡(α)4).
Taking into account that the beam is slender, the following assumption for the stresses can be introduced:
(7)σy=σz=τyz=0.
Considering assumption (7) and the relation (5), the reduced stress-strain relation for beam is obtained:
(8){σxτxzτxy}k=[C~11k0C~16k0C~55k0C~16k0C~66k]{εxγxzγxy}k.
The coefficients of the reduced stiffness matrix are given in the appendix.


## 2.4. Strain Expressions

Geometrically nonlinear deformation is considered and the axial and shear strains are derived from Green’s strain tensor [25], where the longitudinal terms of second order are neglected:
(9)iεx=∂u∂x+12(∂v∂x)2+12(∂w∂x)2,γxz=∂u∂z+∂w∂x+∂v∂x∂v∂z+∂w∂x∂w∂z,γxy=∂u∂y+∂v∂x+∂v∂x∂v∂y+∂w∂x∂w∂y.
A square root of the shear correction λ factor is applied solely to the bending terms of the shear strains. This formulation is preferred here, instead of applying the shear correction factor in Hooke’s law, because of the torsional terms which appear in the shear strains, for which a shear correction factor is not necessary. With this formulation, the expressions from the virtual work of internal forces related to torsion will not contain the shear correction factor, while the ones related to bending will be multiplied by the shear correction factor [19]. In the numerical results, a shear correction factor equal to 5/6 is used. Substituting the displacement field (1) into the expressions for the strains (9) and assuming the shear correction factor λ, the strains become
(10)εx=∂u0∂x-(y-yc)∂ϕz∂x+(z-zc)∂ϕy∂xiiiiiii+ψ∂2θx∂x2+12(∂v0∂x)2+12z2(∂θx∂x)2+12(∂w0∂x)2+12y2(∂θx∂x)2-z∂v0∂x∂θx∂x+y∂w0∂x∂θx∂x,γxz=(∂ψ∂z+y)∂θx∂x+λ∂w0∂x+λϕy-λ∂v0∂xθx,γxy=(∂ψ∂y-z)∂θx∂x+λ∂v0∂x-λϕz+λ∂w0∂xθx.


## 2.5. Equation of Motion

The equation of motion is derived by the principle of virtual work and it is discretized into a system of ordinary differential equation by thep-version finite element method. The p-FEM has several advantages over the h-FEM; the most important is that it requires far fewer degrees of freedom than the h-FEM [26]. One element is used for the laminated composite beam and improvement of the accuracy is achieved by increasing the number of the shape functions. Thus, the functions on the reference line, u0(x,t), v0(x,t), w0(x,t), θx(x,t), ϕy(x,t), and ϕz(x,t), are expressed by shape functions and generalized coordinates in a local coordinate system:
(11)d0={u0(ξ,t)v0(ξ,t)w0(ξ,t)θx(ξ,t)ϕy(ξ,t)ϕz(ξ,t)}=N(ξ)q(t),
where q(t) is the vector of generalized coordinates:
(12)q(t)={qu(t)qv(t)qw(t)qθx(t)qϕy(t)qϕz(t)}
and N(ξ) is the matrix of shape functions:(13)N(ξ)=[Nu(ξ)T000000Nv(ξ)T000000Nw(ξ)T000000Nθx(ξ)T000000Nϕy(ξ)T000000Nϕz(ξ)T].Nu(ξ)T, Nv(ξ)T, Nw(ξ)T, Nθx(ξ)T, Nϕy(ξ)T, and Nϕz(ξ)T are, in this order, the row vectors of longitudinal, transverse along y, transverse along z, torsional, rotational about y, and rotational about z shape functions and ξ∈[-1⋯1] is the local coordinate, ξ=2x/l.The shape functions have to satisfy the geometric boundary conditions. The sets of shape functions, used in previous works for isotropic beams, are implemented here [19, 26]. They are suitable for beams with clamped-clamped boundary conditions. If one wants to apply simply supported or free boundary conditions, additional shape functions, such as Hermite cubics, have to be included.The equation of motion is derived by the principle of virtual work:(14)δWV+δWin+δWE=0,
where WV, Win, and WE are the virtual work of internal, inertia, and external forces due to a virtual displacement δd:
(15)δd={δuδvδw}
and d is the vector of displacement components. The variations of the work of internal, inertia, and external forces are defined as
(16)δWV=-∫VδεTσdV,δWin=-∫VρδdTd¨dV,δWE=∫Vδd0TF0dV,
where δε represents the virtual strains, d¨ is the acceleration of a point of the beam, ρ is the density of the beam, and F0 represents external forces applied on the reference line.Replacing (16) into (14) and introducing mass proportional and frequency dependent damping, the following system of second order ordinary differential equations is obtained:
(17)Mq¨(t)+βωMq˙(t)+K(q(t))q(t)=Fcos⁡⁡(ωt),
where M is the mass matrix obtained from the inertia forces and K(q(t)) is the stiffness matrix obtained from the internal forces; it depends on the vector of generalized coordinates q(t), but it also has constant terms. F represents vector of generalized external forces; a harmonic excitation with frequency ω is applied. β is damping factor, which can be estimated from experiments. In the following work, β=0.01ωl12 is used, where ωl1 is the fundamental frequency. The equation of motion (17) is solved in time domain by Newmark’s method [27] and the resulting nonlinear algebraic system is solved by Newton-Raphson’s method [27], by using numerical computation of the Jacobian. Results of nonlinear modes of vibration are also presented; for that purpose the vector of generalized coordinates is expressed in Fourier series and (17) is transformed into a system of nonlinear algebraic equations, which is solved by the arc-length continuation method.


## 3. Results and Discussion

### 3.1. Free Vibration

The derived model is validated by comparing the natural frequencies of laminated composite beam with results available in the literature, as well as with results obtained by three-dimensional finite element analysis. 10 shape functions are used for each displacement component; that is, the resulting system consists of 60 degrees of freedom (DOF). All layers are assumed to be of equal material, oriented in different angles about the transverse axisz. The assumed material is graphite-epoxy (AS4/3501-6), with the following material properties:
(18)E1=144.8GPa,E2=9.65GPa,G12=G13=4.14GPa,G23=3.45GPa,ν12=ν13=0.3,ν23=0.5,ρ=1389kg/m3.
First, a beam with width b=0.0254 m, height h=0.0254 m, and length l=0.381 m, composed of four layers with asymmetric cross ply [0°/90°/0°/90°], is considered. The bending natural frequencies in xz plane, for different boundary conditions, are presented and compared in Tables 1 and 2. It is noted that, for the case of asymmetric cross ply [0°/90°/0°/90°], the transverse displacement w0 is coupled with the longitudinal displacement u0 even in linear analysis. Furthermore, also due to asymmetry of the laminates, the twist centre and the centroid do not coincide. Hence, additional coupling between the transverse displacement v0 and the torsion θx exists also in linear analysis.Table 1
Natural frequencies (Hz) of clamped-clamped beam of asymmetric cross ply [0°/90°/0°/90°].
Mode
Clamped-Clamped
Present
Reference [8]
Reference [32]
Reference [33]
1
1053.9
1054.4
1062.2
1051.8
2
2596.9
2509.2
2612.0
2505.4
3
4292.1
4281.5
4293.7
4265.1
4
6314.5
6215.8
6309.0
6150.3
5
8271.0
8239.9
8245.7
8062.5Table 2
Natural frequencies (Hz) of clamped simply supported and simply supported simply supported beams of asymmetric cross ply [0°/90°/0°/90°].
Mode
Clamped Simply Supported
Simply Supported Simply Supported
Present
Reference [8]
Reference [32]
Present
Reference [8]
Reference [33]
1
793.3
784.9
818.5
557.9
558.9
557.8
2
2244.2
2223.9
2271.1
1904.7
1907.7
1892.6
3
4062.8
4038.6
4077.9
3782.6
3787.2
3730.9
4
6060.0
6031.5
6071.0
5819.4
5824.8
5669.5
5
8131.6
8104.8
8132.9
7974.9
7980.5
7672.0Three-dimensional finite element software Elmer [28] is used to further validate the derived beam model. The natural frequencies, including the transverse modes in xy plane and the torsional modes, are compared in Tables 3, 4, and 5, for clamped-clamped beams with three different orientations of the laminates, denoted by Case  1 [0°/0°/0°/0°], Case  2 [0°/90°/90°/0°], and Case  3 [0°/90°/0°/90°]. The width and the height are the same, as in the previous example; b=0.0254 m, h=0.0254 m, and the length is l=1.0 m. Fine mesh of quadratic tetrahedrons is generated with Gmsh [29], for the three-dimensional finite element analysis. The mesh consists of 222 246 elements, which results in about 1 million degrees of freedom. The resulting large-scale system is solved on parallel processors, the library MUMPS [30] (multifrontal massively parallel solver), in which a parallel direct sparse solver is used. It was verified that, by reducing the size of the elements, the results obtained by Elmer are converged.Table 3
Natural frequencies (rad/s) of clamped-clamped laminated composite beam with symmetric cross ply [0°/0°/0°/0°], validation with Elmer.
Mode
Case  1 [0°/0°/0°/0°]
Beam 60 DOF
Elmer1 million DOF
Mode shape
Difference %
1
1589.83
1589.05
w
0
0.049
2
1589.83
1589.06
v
0
0.049
3
4198.73
4136.98
w
0
1.493
4
4198.73
4137.00
v
0
1.492
5
5057.49
4999.45
θ
x
1.161
6
7618.51
7596.04
w
0
0.296
7
7618.51
7596.11
v
0
0.295
8
10102.52
10001.40
θ
x
1.011
9
11833.45
11698.23
w
0
1.156
10
11833.45
11698.39
v
0
1.155Table 4
Natural frequencies (rad/s) of clamped-clamped laminated composite beam with symmetric cross ply [0°/90°/90°/0°], validation with Elmer.
Mode
Case  2 [0°/90°/90°/0°]
Beam60 DOF
Elmer1 million DOF
Mode shape
Difference %
1
1188.70
1184.76
v
0
0.333
2
1496.90
1496.34
w
0
0.037
3
3196.79
3146.51
v
0
1.598
4
3958.58
3900.35
w
0
1.493
5
4876.05
4821.07
θ
x
1.141
6
5979.48
5898.63
v
0
1.371
7
7198.00
7169.94
w
0
0.391
8
9498.68
9269.03
v
0
2.478
9
9740.08
9644.57
θ
x
0.990
10
11197.12
11053.57
w
0
1.299Table 5
Natural frequencies (rad/s) of clamped-clamped laminated composite beam with asymmetric cross ply [0°/90°/0°/90°], validation with Elmer.
Mode
Case  3 [0°/90°/0°/90°]
Beam60 DOF
Elmer1 million DOF
Mode shape
Difference %
1
1102.01
1096.80
(w0,u0)
0.476
2
1188.69
1186.15
(v0,θx)
0.215
3
2967.67
2912.99
(w0,u0)
1.877
4
3196.75
3153.58
(v0,θx)
1.369
5
4932.52
4874.94
θ
x
1.181
6
5564.63
5461.35
(w0,u0)
1.891
7
5979.39
5919.25
(v0,θx)
1.016
8
8856.42
8583.29
(w0,u0)
3.182
9
9498.00
9313.68
(v0,θx)
1.979
10
9851.14
9751.74
θ
x
1.019The mode shapes of the first four natural frequencies, for the beam with asymmetric cross ply [0°/90°/0°/90°], which introduces bending-longitudinal and bending-torsional coupling in linear analysis, are shown in Figure2. The amplitudes are normalized, for better interpretation of the mode shapes and the differences between them.First four mode shapes (with normalized amplitudes) of laminated asymmetric cross ply [0°/90°/90°/0°] with clamped-clamped boundary conditions,l=1.0 m. (a) First mode; —w0, - - - -u0, (b) second mode: —v0, - - - -θx, (c) third mode: —w0, - - - -u0, and (d) fourth mode: —v0, - - - -θx.
(a)
(b)
(c)
(d)There are no differences in the orientations of the layers for Case  2 and Case  3; the layers are in pairs, oriented in the same directions. The difference comes from the vertical order of the layers. The first bending frequency inxz plane is 1496.34 rad/s for Case  2, and for Case  3 it is significantly smaller—1096.80 rad/s. The associated mode shapes of both cases are purely in xz plane. The bending-longitudinal coupling, which appears due to asymmetric cross ply in Case  3, is the reason for the significant change of the fundamental frequency of bending in xz plane. On the other hand, the first frequency of bending in xy plane couples with torsion for layers of Case  3, but in both cases the frequencies remain close: 1184.76 rad/s for Case  2 and 1186.15 rad/s for Case  3.


### 3.2. Static Deformations

In the previous section only the linear part of the model was validated. Here, the nonlinear model is validated by comparing the static displacements of composite beams with results from the literature [31]. The beam is with dimensions as follows: h=0.2794 m, b=0.2794 m, and l=6.35 m, with the following material properties:
(19)E1=241.5GPa,E2=E3=18.89GPa,G12=G13=5.18GPa,G23=3.45GPa,ν12=ν13=0.24,ν23=0.25,ρ=2015kg/m3.
Two different beams are considered, with symmetrical and with nonsymmetrical layers about the middle line. The first beam is with three layers at orientation [90°/0°/90°], with thicknesses 0.25 h, 0.5 h, and 0.25 h. The second beam is with two layers at orientation [90°/0°], with thicknesses 0.5 h and 0.5 h. The beam is assumed to be clamped-clamped with static force: F(x)=Acos⁡⁡(πx/l), x=-l/2⋯l/2. The results of transverse displacement w0, for different amplitudes of the static force, are presented in Tables 6 and 7. The results are in agreement with the ones from [31].Table 6
Static deformations of the transverse displacementw0, layers [90°/0°/90°].
Force(106 N)
1.5
2
2.5
3.5
4.5
8
10
15
Current model
−0.0628
−0.0792
−0.0935
−0.1176
−0.1374
−0.1883
−0.2101
−0.2533
Reference [31]
−0.0625
−0.0788
−0.0932
−0.1173
−0.1372
−0.1882
−0.2102
−0.2535
Difference %
0.44
0.44
0.34
0.25
0.18
0.05
0.03
0.09Table 7
Static deformations of the transverse displacementw0, layers [90°/0°].
Force(106 N)
5
7
10
12
13
14
15
17
Current model
−0.1339
−0.1641
−0.1993
−0.2187
−0.2275
−0.2358
−0.2436
−0.2583
Reference [31]
−0.1341
−0.1645
−0.1999
−0.2193
−0.2281
−0.2365
−0.2444
−0.2591
Difference %
0.17
0.23
0.29
0.28
0.28
0.28
0.31
0.30


### 3.3. Warping Function

The warping function, which is computed numerically by solving (2), is presented in Figure 3 for two different orientations of the layers, in order to demonstrate the differences in the longitudinal deformation due to warping. The torsional rigidity also changes with the orientation of the layers. It is defined to be
(20)GIt=∑k=1M∫Ωk(C~55k(∂ψ∂z+y)2+C~66k(∂ψ∂y-z)2)dΩk.
The torsional rigidity of the cross section of Case  1 (Figure 3(a)) is obtained to be 242.24 Nm2, while the torsional rigidity of the second cross section, which is Case  2 (Figure 3(b)), is 225.17 Nm2. These differences of the torsional rigidity are due to different shear moduli of the layers, which are result of the angle orientation and due to the warping distribution. The differences of the torsional rigidities result in different torsional natural frequencies. For the beam of Case  1, the fundamental torsional frequency is 5057.49 rad/s and for the beam of Case  2 4876.05 rad/s (Tables 3 and 4).Warping deformation for beam with rectangular cross section. (a) Layers orientation [0°/0°/0°/0°] and (b) layers orientation [0°/90°/90°/0°].
(a)
(b)The numerical solution of the warping function, together with numerical computation of the cross-sectional properties, will allow one to model composite beams with arbitrary cross sections.


### 3.4. Nonlinear Free Vibration

The main branches of the nonlinear free vibrations are presented in this section for the three cases of laminate orientations, discussed in Section3.1. The beam is with the same dimensions as in Section 3.1, b=0.0254 m, h=0.0254 m, and l=1.0 m, and the orientations of the laminates are denoted by Case  1: [0°/0°/0°/0°], Case  2: [0°/90°/90°/0°], and Case  3: [0°/90°/0°/90°]. The same orthotropic material graphite-epoxy (AS4/3501-6) is assumed.For the nonlinear free vibration analysis, damping and external force are not considered, the vector of generalized coordinates is expressed in Fourier series, and the harmonic balance method is applied. Harmonics up to third order are considered, so the vector of generalized coordinates is expressed as(21)q(t)=12q0+∑k=13qicos⁡⁡(iωt).
As a result, a system of algebraic nonlinear equations is obtained, with unknowns, the coefficients of the harmonics, qi, i=0,…,3, and the frequency of vibration, ω. This system is solved in the frequency domain by the arc-length continuation method.The first bending linear modes and the corresponding linear frequencies, of each transverse displacement,v0 and w0, are used to start the continuation procedure. It is noted that these linear modes correspond to the first two natural frequencies, presented in Tables 3–5. The main branches of the three cases are presented in Figure 4. In addition to the comparison of the natural frequencies, given in Tables 3–5, Figure 4 shows the influence of the orientation and the position of the layers on the nonlinear normal modes. On Figure 4(b) the main branches of Cases  2 and 3 almost coincide. This is a consequence from the linear frequencies of the first bending mode of transverse displacement v0, which are close in both cases (Tables 4 and 5). The nonlinear normal modes also remain close. In Case  3, where the layers are asymmetrical, the torsion is different from zero. Hence, the mode shapes are different, even though the amplitudes of the transverse displacement are close and the linear natural frequencies are almost equal. The shapes of vibration, for frequency ω=2000rad/s, are presented in Figure 5.Main branches of first bending modes of free vibration of composite beams.W1: amplitude of the first harmonic of the transverse displacement w0, V1: amplitude of the first harmonic of the transverse displacement v0, red: Case  1, black: Case  2, and blue: Case  3.
(a)
(b)Shapes of vibration of (a) transverse displacementv0 and (b) torsion θx, for frequency of vibration ω=2000rad/s, from bifurcation diagram shown on Figure 4(b). Black: Case  2 and blue: Case  3 (torsion exists only for Case  3). V1: shape of first harmonic function of v0, V3: shape of third harmonic of v0, Θ1: shape of first harmonic of θx, and Θ3: shape of third harmonic of θx.
(a)
(b)Further to present the influence of the orientation and the position of the layers on the nonlinear modes of vibration, the nondimensional amplitudes of vibration of the three cases are compared. The results are presented in Table8. It can be seen that, for the same nondimensional frequency of vibration, the amplitudes of the first harmonics are different.Table 8
Comparison of the nondimensional amplitudes of the first harmonics (W1/h) of the transverse displacement w0, presented in the bifurcation diagrams, Figure 3(a).
ω
/
ω
l
Case  1
Case  2
Case  3
1.1
0.607
0.783
0.569
1.2
0.883
1.137
0.827
1.3
1.110
1.430
1.039
1.4
1.314
1.693
1.231
1.5
1.506
1.940
1.409
1.6
1.690
2.177
1.580
1.7
1.867
2.405
1.746
1.8
2.039
2.627
1.908
1.9
2.205
2.841
2.066
2
2.364
3.046
2.220


### 3.5. Dynamic Response

The dynamics responses of laminated beams due to harmonic excitations are investigated in this section, with particular attention to the influence of the orientation of the laminates on the amplitudes of vibration. Considering the above investigation on the natural frequencies, beam with dimensionsb=0.0254 m, h=0.0254 m, and l=1.0 m is considered and composed of four layers of graphite-epoxy (AS4/3501-6), with orientations defined as in Cases  1 to 3.In order to excite all displacements, a combined harmonic force is applied, in the middle of the beam, in both transverse directions plus a moment. The excitation frequency is 1200 rad/s; it is close to the fundamental natural frequency of beam with orientation of layers defined as Case  2. The excitation frequency is also close to the first and second natural frequencies of beam with layers orientation from Case  3. The transverse forces are chosen to be with equal amplitudes,Fy=Fz=10000cos⁡⁡(1200t), and the moment is Mx=10cos⁡⁡(1200t).The transient and the steady-state time responses, of both transverse displacements and torsion, are presented in Figures6, 7, and 8. The phase plot of the torsion is also given in Figure 8. Even though the external transverse forces are equal in amplitudes, the transverse responses become different just by changing the orientation of the layers. Furthermore, the position of the laminates, in the cases of equally oriented laminates such as Cases  2 and 3, is also important. The amplitude of transverse displacement w0 becomes bigger when orientation of Case  3 is used in comparison to Case  2. The torsional response also undergoes changes, due to the orientation of the laminates. The torsional response changes not only its amplitude of vibration, but also the harmonics involved in the periodic steady response, as can be seen from Figure 8. The appearance of the constant term in Fourier expansion, as well as higher harmonics, is obvious for the torsional response from Case  2.Transient time responses of dimensionless transverse displacements: (a)v0 and (b) w0 and (c) torsion due to external force Fy=Fz=10000cos⁡⁡(1200t) and Mx=10cos⁡⁡(1200t), applied on the middle of the beam. Red: Case  1, black: Case  2, and blue: Case  3.
(a)
(b)
(c)Time responses of dimensionless transverse displacements: (a)v0 and (b) w0 due to external force Fy=Fz=10000cos⁡⁡(1200t) and Mx=10cos⁡⁡(1200t), applied on the middle of the beam. Red: Case  1, black: Case  2, and blue: Case  3.
(a)
(b)(a) Time responses and (b) phase plots of torsionθx due to external force Fy=Fz=10000cos⁡⁡(1200t) and Mx=10cos⁡⁡(1200t), applied on the middle of the beam. Red: Case  1, black: Case  2, and blue: Case  3.
(a)
(b)Finally, to point out the influence of the geometrically nonlinear terms, on the dynamic response of the laminated beam, results obtained with and without geometrical nonlinearity are compared in Figure9. Generally, a model with geometrical nonlinear terms results into smaller amplitudes of the displacements than the linear model. Nevertheless, the absolute amplitude of the torsional response for Case  2 remains similar for the linear and the nonlinear models (Figure 9(b)).Time responses of (a) dimensionless transverse displacementw0 and (b) torsion θx, due to external force Fy=Fz=10000cos⁡⁡(1200t) and Mx=10cos⁡⁡(1200t), applied on the middle of the beam. Line: nonlinear model and dash: linear model.
(a)
(b)


## 3.1. Free Vibration

The derived model is validated by comparing the natural frequencies of laminated composite beam with results available in the literature, as well as with results obtained by three-dimensional finite element analysis. 10 shape functions are used for each displacement component; that is, the resulting system consists of 60 degrees of freedom (DOF). All layers are assumed to be of equal material, oriented in different angles about the transverse axisz. The assumed material is graphite-epoxy (AS4/3501-6), with the following material properties:
(18)E1=144.8GPa,E2=9.65GPa,G12=G13=4.14GPa,G23=3.45GPa,ν12=ν13=0.3,ν23=0.5,ρ=1389kg/m3.
First, a beam with width b=0.0254 m, height h=0.0254 m, and length l=0.381 m, composed of four layers with asymmetric cross ply [0°/90°/0°/90°], is considered. The bending natural frequencies in xz plane, for different boundary conditions, are presented and compared in Tables 1 and 2. It is noted that, for the case of asymmetric cross ply [0°/90°/0°/90°], the transverse displacement w0 is coupled with the longitudinal displacement u0 even in linear analysis. Furthermore, also due to asymmetry of the laminates, the twist centre and the centroid do not coincide. Hence, additional coupling between the transverse displacement v0 and the torsion θx exists also in linear analysis.Table 1
Natural frequencies (Hz) of clamped-clamped beam of asymmetric cross ply [0°/90°/0°/90°].
Mode
Clamped-Clamped
Present
Reference [8]
Reference [32]
Reference [33]
1
1053.9
1054.4
1062.2
1051.8
2
2596.9
2509.2
2612.0
2505.4
3
4292.1
4281.5
4293.7
4265.1
4
6314.5
6215.8
6309.0
6150.3
5
8271.0
8239.9
8245.7
8062.5Table 2
Natural frequencies (Hz) of clamped simply supported and simply supported simply supported beams of asymmetric cross ply [0°/90°/0°/90°].
Mode
Clamped Simply Supported
Simply Supported Simply Supported
Present
Reference [8]
Reference [32]
Present
Reference [8]
Reference [33]
1
793.3
784.9
818.5
557.9
558.9
557.8
2
2244.2
2223.9
2271.1
1904.7
1907.7
1892.6
3
4062.8
4038.6
4077.9
3782.6
3787.2
3730.9
4
6060.0
6031.5
6071.0
5819.4
5824.8
5669.5
5
8131.6
8104.8
8132.9
7974.9
7980.5
7672.0Three-dimensional finite element software Elmer [28] is used to further validate the derived beam model. The natural frequencies, including the transverse modes in xy plane and the torsional modes, are compared in Tables 3, 4, and 5, for clamped-clamped beams with three different orientations of the laminates, denoted by Case  1 [0°/0°/0°/0°], Case  2 [0°/90°/90°/0°], and Case  3 [0°/90°/0°/90°]. The width and the height are the same, as in the previous example; b=0.0254 m, h=0.0254 m, and the length is l=1.0 m. Fine mesh of quadratic tetrahedrons is generated with Gmsh [29], for the three-dimensional finite element analysis. The mesh consists of 222 246 elements, which results in about 1 million degrees of freedom. The resulting large-scale system is solved on parallel processors, the library MUMPS [30] (multifrontal massively parallel solver), in which a parallel direct sparse solver is used. It was verified that, by reducing the size of the elements, the results obtained by Elmer are converged.Table 3
Natural frequencies (rad/s) of clamped-clamped laminated composite beam with symmetric cross ply [0°/0°/0°/0°], validation with Elmer.
Mode
Case  1 [0°/0°/0°/0°]
Beam 60 DOF
Elmer1 million DOF
Mode shape
Difference %
1
1589.83
1589.05
w
0
0.049
2
1589.83
1589.06
v
0
0.049
3
4198.73
4136.98
w
0
1.493
4
4198.73
4137.00
v
0
1.492
5
5057.49
4999.45
θ
x
1.161
6
7618.51
7596.04
w
0
0.296
7
7618.51
7596.11
v
0
0.295
8
10102.52
10001.40
θ
x
1.011
9
11833.45
11698.23
w
0
1.156
10
11833.45
11698.39
v
0
1.155Table 4
Natural frequencies (rad/s) of clamped-clamped laminated composite beam with symmetric cross ply [0°/90°/90°/0°], validation with Elmer.
Mode
Case  2 [0°/90°/90°/0°]
Beam60 DOF
Elmer1 million DOF
Mode shape
Difference %
1
1188.70
1184.76
v
0
0.333
2
1496.90
1496.34
w
0
0.037
3
3196.79
3146.51
v
0
1.598
4
3958.58
3900.35
w
0
1.493
5
4876.05
4821.07
θ
x
1.141
6
5979.48
5898.63
v
0
1.371
7
7198.00
7169.94
w
0
0.391
8
9498.68
9269.03
v
0
2.478
9
9740.08
9644.57
θ
x
0.990
10
11197.12
11053.57
w
0
1.299Table 5
Natural frequencies (rad/s) of clamped-clamped laminated composite beam with asymmetric cross ply [0°/90°/0°/90°], validation with Elmer.
Mode
Case  3 [0°/90°/0°/90°]
Beam60 DOF
Elmer1 million DOF
Mode shape
Difference %
1
1102.01
1096.80
(w0,u0)
0.476
2
1188.69
1186.15
(v0,θx)
0.215
3
2967.67
2912.99
(w0,u0)
1.877
4
3196.75
3153.58
(v0,θx)
1.369
5
4932.52
4874.94
θ
x
1.181
6
5564.63
5461.35
(w0,u0)
1.891
7
5979.39
5919.25
(v0,θx)
1.016
8
8856.42
8583.29
(w0,u0)
3.182
9
9498.00
9313.68
(v0,θx)
1.979
10
9851.14
9751.74
θ
x
1.019The mode shapes of the first four natural frequencies, for the beam with asymmetric cross ply [0°/90°/0°/90°], which introduces bending-longitudinal and bending-torsional coupling in linear analysis, are shown in Figure2. The amplitudes are normalized, for better interpretation of the mode shapes and the differences between them.First four mode shapes (with normalized amplitudes) of laminated asymmetric cross ply [0°/90°/90°/0°] with clamped-clamped boundary conditions,l=1.0 m. (a) First mode; —w0, - - - -u0, (b) second mode: —v0, - - - -θx, (c) third mode: —w0, - - - -u0, and (d) fourth mode: —v0, - - - -θx.
(a)
(b)
(c)
(d)There are no differences in the orientations of the layers for Case  2 and Case  3; the layers are in pairs, oriented in the same directions. The difference comes from the vertical order of the layers. The first bending frequency inxz plane is 1496.34 rad/s for Case  2, and for Case  3 it is significantly smaller—1096.80 rad/s. The associated mode shapes of both cases are purely in xz plane. The bending-longitudinal coupling, which appears due to asymmetric cross ply in Case  3, is the reason for the significant change of the fundamental frequency of bending in xz plane. On the other hand, the first frequency of bending in xy plane couples with torsion for layers of Case  3, but in both cases the frequencies remain close: 1184.76 rad/s for Case  2 and 1186.15 rad/s for Case  3.


## 3.2. Static Deformations

In the previous section only the linear part of the model was validated. Here, the nonlinear model is validated by comparing the static displacements of composite beams with results from the literature [31]. The beam is with dimensions as follows: h=0.2794 m, b=0.2794 m, and l=6.35 m, with the following material properties:
(19)E1=241.5GPa,E2=E3=18.89GPa,G12=G13=5.18GPa,G23=3.45GPa,ν12=ν13=0.24,ν23=0.25,ρ=2015kg/m3.
Two different beams are considered, with symmetrical and with nonsymmetrical layers about the middle line. The first beam is with three layers at orientation [90°/0°/90°], with thicknesses 0.25 h, 0.5 h, and 0.25 h. The second beam is with two layers at orientation [90°/0°], with thicknesses 0.5 h and 0.5 h. The beam is assumed to be clamped-clamped with static force: F(x)=Acos⁡⁡(πx/l), x=-l/2⋯l/2. The results of transverse displacement w0, for different amplitudes of the static force, are presented in Tables 6 and 7. The results are in agreement with the ones from [31].Table 6
Static deformations of the transverse displacementw0, layers [90°/0°/90°].
Force(106 N)
1.5
2
2.5
3.5
4.5
8
10
15
Current model
−0.0628
−0.0792
−0.0935
−0.1176
−0.1374
−0.1883
−0.2101
−0.2533
Reference [31]
−0.0625
−0.0788
−0.0932
−0.1173
−0.1372
−0.1882
−0.2102
−0.2535
Difference %
0.44
0.44
0.34
0.25
0.18
0.05
0.03
0.09Table 7
Static deformations of the transverse displacementw0, layers [90°/0°].
Force(106 N)
5
7
10
12
13
14
15
17
Current model
−0.1339
−0.1641
−0.1993
−0.2187
−0.2275
−0.2358
−0.2436
−0.2583
Reference [31]
−0.1341
−0.1645
−0.1999
−0.2193
−0.2281
−0.2365
−0.2444
−0.2591
Difference %
0.17
0.23
0.29
0.28
0.28
0.28
0.31
0.30


## 3.3. Warping Function

The warping function, which is computed numerically by solving (2), is presented in Figure 3 for two different orientations of the layers, in order to demonstrate the differences in the longitudinal deformation due to warping. The torsional rigidity also changes with the orientation of the layers. It is defined to be
(20)GIt=∑k=1M∫Ωk(C~55k(∂ψ∂z+y)2+C~66k(∂ψ∂y-z)2)dΩk.
The torsional rigidity of the cross section of Case  1 (Figure 3(a)) is obtained to be 242.24 Nm2, while the torsional rigidity of the second cross section, which is Case  2 (Figure 3(b)), is 225.17 Nm2. These differences of the torsional rigidity are due to different shear moduli of the layers, which are result of the angle orientation and due to the warping distribution. The differences of the torsional rigidities result in different torsional natural frequencies. For the beam of Case  1, the fundamental torsional frequency is 5057.49 rad/s and for the beam of Case  2 4876.05 rad/s (Tables 3 and 4).Warping deformation for beam with rectangular cross section. (a) Layers orientation [0°/0°/0°/0°] and (b) layers orientation [0°/90°/90°/0°].
(a)
(b)The numerical solution of the warping function, together with numerical computation of the cross-sectional properties, will allow one to model composite beams with arbitrary cross sections.


## 3.4. Nonlinear Free Vibration

The main branches of the nonlinear free vibrations are presented in this section for the three cases of laminate orientations, discussed in Section3.1. The beam is with the same dimensions as in Section 3.1, b=0.0254 m, h=0.0254 m, and l=1.0 m, and the orientations of the laminates are denoted by Case  1: [0°/0°/0°/0°], Case  2: [0°/90°/90°/0°], and Case  3: [0°/90°/0°/90°]. The same orthotropic material graphite-epoxy (AS4/3501-6) is assumed.For the nonlinear free vibration analysis, damping and external force are not considered, the vector of generalized coordinates is expressed in Fourier series, and the harmonic balance method is applied. Harmonics up to third order are considered, so the vector of generalized coordinates is expressed as(21)q(t)=12q0+∑k=13qicos⁡⁡(iωt).
As a result, a system of algebraic nonlinear equations is obtained, with unknowns, the coefficients of the harmonics, qi, i=0,…,3, and the frequency of vibration, ω. This system is solved in the frequency domain by the arc-length continuation method.The first bending linear modes and the corresponding linear frequencies, of each transverse displacement,v0 and w0, are used to start the continuation procedure. It is noted that these linear modes correspond to the first two natural frequencies, presented in Tables 3–5. The main branches of the three cases are presented in Figure 4. In addition to the comparison of the natural frequencies, given in Tables 3–5, Figure 4 shows the influence of the orientation and the position of the layers on the nonlinear normal modes. On Figure 4(b) the main branches of Cases  2 and 3 almost coincide. This is a consequence from the linear frequencies of the first bending mode of transverse displacement v0, which are close in both cases (Tables 4 and 5). The nonlinear normal modes also remain close. In Case  3, where the layers are asymmetrical, the torsion is different from zero. Hence, the mode shapes are different, even though the amplitudes of the transverse displacement are close and the linear natural frequencies are almost equal. The shapes of vibration, for frequency ω=2000rad/s, are presented in Figure 5.Main branches of first bending modes of free vibration of composite beams.W1: amplitude of the first harmonic of the transverse displacement w0, V1: amplitude of the first harmonic of the transverse displacement v0, red: Case  1, black: Case  2, and blue: Case  3.
(a)
(b)Shapes of vibration of (a) transverse displacementv0 and (b) torsion θx, for frequency of vibration ω=2000rad/s, from bifurcation diagram shown on Figure 4(b). Black: Case  2 and blue: Case  3 (torsion exists only for Case  3). V1: shape of first harmonic function of v0, V3: shape of third harmonic of v0, Θ1: shape of first harmonic of θx, and Θ3: shape of third harmonic of θx.
(a)
(b)Further to present the influence of the orientation and the position of the layers on the nonlinear modes of vibration, the nondimensional amplitudes of vibration of the three cases are compared. The results are presented in Table8. It can be seen that, for the same nondimensional frequency of vibration, the amplitudes of the first harmonics are different.Table 8
Comparison of the nondimensional amplitudes of the first harmonics (W1/h) of the transverse displacement w0, presented in the bifurcation diagrams, Figure 3(a).
ω
/
ω
l
Case  1
Case  2
Case  3
1.1
0.607
0.783
0.569
1.2
0.883
1.137
0.827
1.3
1.110
1.430
1.039
1.4
1.314
1.693
1.231
1.5
1.506
1.940
1.409
1.6
1.690
2.177
1.580
1.7
1.867
2.405
1.746
1.8
2.039
2.627
1.908
1.9
2.205
2.841
2.066
2
2.364
3.046
2.220


## 3.5. Dynamic Response

The dynamics responses of laminated beams due to harmonic excitations are investigated in this section, with particular attention to the influence of the orientation of the laminates on the amplitudes of vibration. Considering the above investigation on the natural frequencies, beam with dimensionsb=0.0254 m, h=0.0254 m, and l=1.0 m is considered and composed of four layers of graphite-epoxy (AS4/3501-6), with orientations defined as in Cases  1 to 3.In order to excite all displacements, a combined harmonic force is applied, in the middle of the beam, in both transverse directions plus a moment. The excitation frequency is 1200 rad/s; it is close to the fundamental natural frequency of beam with orientation of layers defined as Case  2. The excitation frequency is also close to the first and second natural frequencies of beam with layers orientation from Case  3. The transverse forces are chosen to be with equal amplitudes,Fy=Fz=10000cos⁡⁡(1200t), and the moment is Mx=10cos⁡⁡(1200t).The transient and the steady-state time responses, of both transverse displacements and torsion, are presented in Figures6, 7, and 8. The phase plot of the torsion is also given in Figure 8. Even though the external transverse forces are equal in amplitudes, the transverse responses become different just by changing the orientation of the layers. Furthermore, the position of the laminates, in the cases of equally oriented laminates such as Cases  2 and 3, is also important. The amplitude of transverse displacement w0 becomes bigger when orientation of Case  3 is used in comparison to Case  2. The torsional response also undergoes changes, due to the orientation of the laminates. The torsional response changes not only its amplitude of vibration, but also the harmonics involved in the periodic steady response, as can be seen from Figure 8. The appearance of the constant term in Fourier expansion, as well as higher harmonics, is obvious for the torsional response from Case  2.Transient time responses of dimensionless transverse displacements: (a)v0 and (b) w0 and (c) torsion due to external force Fy=Fz=10000cos⁡⁡(1200t) and Mx=10cos⁡⁡(1200t), applied on the middle of the beam. Red: Case  1, black: Case  2, and blue: Case  3.
(a)
(b)
(c)Time responses of dimensionless transverse displacements: (a)v0 and (b) w0 due to external force Fy=Fz=10000cos⁡⁡(1200t) and Mx=10cos⁡⁡(1200t), applied on the middle of the beam. Red: Case  1, black: Case  2, and blue: Case  3.
(a)
(b)(a) Time responses and (b) phase plots of torsionθx due to external force Fy=Fz=10000cos⁡⁡(1200t) and Mx=10cos⁡⁡(1200t), applied on the middle of the beam. Red: Case  1, black: Case  2, and blue: Case  3.
(a)
(b)Finally, to point out the influence of the geometrically nonlinear terms, on the dynamic response of the laminated beam, results obtained with and without geometrical nonlinearity are compared in Figure9. Generally, a model with geometrical nonlinear terms results into smaller amplitudes of the displacements than the linear model. Nevertheless, the absolute amplitude of the torsional response for Case  2 remains similar for the linear and the nonlinear models (Figure 9(b)).Time responses of (a) dimensionless transverse displacementw0 and (b) torsion θx, due to external force Fy=Fz=10000cos⁡⁡(1200t) and Mx=10cos⁡⁡(1200t), applied on the middle of the beam. Line: nonlinear model and dash: linear model.
(a)
(b)


## 4. Conclusion

A model for 3D laminated composite beams was presented. The equation of motion was derived by the principle of virtual work and thep-FEM. The model assumes Timoshenko’s theory for bending and considers that under torsion the cross section rotates as a rigid body and deforms longitudinally due to warping. The warping function was included in the model; it was obtained preliminarily by the finite element method. The inclusion of the warping function is essential, for correct torsional and bending-torsional modes, because it influences the torsional rigidity of the beam.The beam model was validated with available results from the literature and, with equivalent beam structure, discretized by three-dimensional finite elements. It was shown that the beam model, discretized by thep-FEM, gives results in agreement with the large-scale model but with far fewer degrees of freedom. The couplings between different displacement components, in the case of asymmetrical laminates, were pointed out.Nonlinear free vibrations of composite beams were compared for different orientations and positions of the layers. It was shown that changes of the orientation or the positions of the layers can lead to significant changes in the nonlinear modes of vibration.Finally, the dynamic steady-state responses of the nonlinear model were presented. It was shown that the amplitudes of vibration change not only with the orientation of the laminates, but also with the vertical order of the laminates. The torsional response undergoes changes in its Fourier expansion.



---

*Source: 892782-2014-04-24.xml*