Nonlinear Vibrations of 3 D Laminated Composite Beams

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 consideringGreen’s strain tensor.The equation ofmotion 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.


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][2][3][4][5][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 axialbending 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 -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 threedimensional 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.

Mathematical Model
A laminated composite beam is considered (Figure 1); 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 .The length, the width, and the height of the beam are denoted by , , and ℎ, respectively.

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-bendingtorsional 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 (, , , ) and the transverse displacements V(, , , ) and (, , , ), which are functions of space coordinates and time, are expressed by functions of the longitudinal coordinate  in the following way: where the subscript "0" represents the reference line, which is the line of the twist centres of the cross section, (  ,   ) is the position of the centroid of the cross section,   is the rotation of the cross section about the longitudinal axis ,   and   denote rotations of the cross section about  and  axes, respectively, and (, ) is the warping function.

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: where Ω  is the cross section of the th layer and Γ  is its contour and  is the number of layers.C 55 and C 66 are the shear modulus of the material of the th layer from the reduced stiffness matrix.  and   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−1 66 = C−1 55 = 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: where   ,  = 1, . . ., , are the weighting functions.

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 The stress-strain relations of laminate , rotated on angle  about the  axis, are expressed as where the coefficients of the stiffness matrix of ( 5) are given to be [23]: +  66 (sin () 4 + cos () 4 ) .
Taking into account that the beam is slender, the following assumption for the stresses can be introduced: Considering assumption (7) and the relation ( 5), the reduced stress-strain relation for beam is obtained: The coefficients of the reduced stiffness matrix are given in the appendix.

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

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 the -version finite element method.The -FEM has several advantages over the ℎ-FEM; the most important is that it requires far fewer degrees of freedom than the ℎ-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,  0 (, ), V 0 (, ),  0 (, ),   (, ),   (, ), and   (, ), are expressed by shape functions and generalized coordinates in a local coordinate system: where q() is the vector of generalized coordinates: and N() is the matrix of shape functions: , and N  z ()  are, in this order, the row vectors of longitudinal, transverse along , transverse along , torsional, rotational about , and rotational about  shape functions and  ∈ [−1 ⋅ ⋅ ⋅ 1] is the local coordinate,  = 2/.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: where  V ,  in , and  E are the virtual work of internal, inertia, and external forces due to a virtual displacement d: and d is the vector of displacement components.The variations of the work of internal, inertia, and external forces are defined as where  represents the virtual strains, d is the acceleration of a point of the beam,  is the density of the beam, and F 0 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: where M is the mass matrix obtained from the inertia forces and K(q()) is the stiffness matrix obtained from the internal forces; it depends on the vector of generalized coordinates q(), 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 2  1 is used, where   1 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 arclength continuation method.

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 axis .The
The bending natural frequencies in  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  0 is coupled with the longitudinal displacement  0 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 V 0 and the torsion   exists also in linear analysis.Three-dimensional finite element software Elmer [28] is used to further validate the derived beam model.The natural frequencies, including the transverse modes in  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;  = 0.0254 m, ℎ = 0.0254 m, and the length is  = 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.
The 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 Figure 2. The amplitudes are normalized, for better interpretation of the mode shapes and the differences between them.
There are no differences in the orientations of the layers for Case 2 and Case 3; the layers are in pairs, oriented in   ( 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() = A cos(/),  = −/2 ⋅ ⋅ ⋅ /2.The results of transverse displacement  0 , for different amplitudes of the static force, are presented in Tables 6 and 7.The results are in agreement with the ones from [31].

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 The torsional rigidity of the cross section of Case 1 (Figure 3(a)) is obtained to be 242.24Nm 2 , while the  torsional rigidity of the second cross section, which is Case 2 (Figure 3(b)), is 225.17Nm 2 .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.49rad/s and for the beam of Case 2 4876.05rad/s (Tables 3 and 4).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.
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 As a result, a system of algebraic nonlinear equations is obtained, with unknowns, the coefficients of the harmonics, q  ,  = 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, V 0 and  0 , 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 V 0 , 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  = 2000 rad/s, are presented in Figure 5.
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 Table 8.It can be seen that, for the same nondimensional frequency of vibration, the amplitudes of the first harmonics are different.

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 dimensions  = 0.0254 m, ℎ = 0.0254 m, and  = 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,   =   = 10000 cos(1200), and the moment is   = 10 cos(1200).
The transient and the steady-state time responses, of both transverse displacements and torsion, are presented in Figures 6, 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  0 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.
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 Figure 9. Generally, a model with geometrical nonlinear terms results into smaller amplitudes .Black: Case 2 and blue: Case 3 (torsion exists only for Case 3). 1 : shape of first harmonic function of V 0 ,  3 : shape of third harmonic of V 0 , Θ 1 : shape of first harmonic of   , and Θ 3 : shape of third harmonic of   . 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)).

Conclusion
A model for 3D laminated composite beams was presented.The equation of motion was derived by the principle of virtual work and the -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 the -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.  are zero.These coefficients are expressed in terms of the coefficients defined in equation ( 6

Figure 4 :
Figure 4: Main branches of first bending modes of free vibration of composite beams. 1 : amplitude of the first harmonic of the transverse displacement  0 ,  1 : amplitude of the first harmonic of the transverse displacement V 0 , red: Case 1, black: Case 2, and blue: Case 3.

Figure 5 :
Figure 5: Shapes of vibration of (a) transverse displacement V 0 and (b) torsion   , for frequency of vibration  = 2000 rad/s, from bifurcation diagram shown on Figure 4(b).Black: Case 2 and blue: Case 3 (torsion exists only for Case 3). 1 : shape of first harmonic function of V 0 ,  3 : shape of third harmonic of V 0 , Θ 1 : shape of first harmonic of   , and Θ 3 : shape of third harmonic of   .

Figure 9 :
Figure 9: Time responses of (a) dimensionless transverse displacement  0 and (b) torsion   , due to external force   =   = 10000 cos(1200) and   = 10 cos(1200), applied on the middle of the beam.Line: nonlinear model and dash: linear model.
the same directions.The difference comes from the vertical order of the layers.The first bending frequency in  plane is 1496.34rad/s for Case 2, and for Case 3 it is significantly smaller-1096.80rad/s.The associated mode shapes of both cases are purely in  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  plane.On the other hand, the first frequency of bending in  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.15rad/s for Case 3.

Table 8 :
Comparison of the nondimensional amplitudes of the first harmonics ( 1 /ℎ) of the transverse displacement  0 , presented in the bifurcation diagrams, Figure3(a).