A beam finite element for analysis of composite beams with the inclusion of bend-twist coupling

Accepted Manuscript A Beam Finite Element for Analysis of Composite Beams with the Inclusion of Bend-Twist Coupling Pavel Babuska, Richard Wiebe, Michael R. Motley PII: DOI: Reference:

S0263-8223(17)32056-1 https://doi.org/10.1016/j.compstruct.2018.01.036 COST 9276

To appear in:

Composite Structures

Received Date: Revised Date: Accepted Date:

3 July 2017 20 November 2017 9 January 2018

Please cite this article as: Babuska, P., Wiebe, R., Motley, M.R., A Beam Finite Element for Analysis of Composite Beams with the Inclusion of Bend-Twist Coupling, Composite Structures (2018), doi: https://doi.org/10.1016/ j.compstruct.2018.01.036

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

A Beam Finite Element for Analysis of Composite Beams with the Inclusion of Bend-Twist Coupling Pavel Babuskaa , Richard Wiebea , Michael R. Motleya a

Department of Civil & Environmental Engineering, University of Washington, Seattle, WA

Abstract Material-induced bend-twist coupling in laminated composite beams has seen applications in engineered structures for decades, including airplane wings and turbine blades. Numerical studies and analytical formulations of the dynamics of bend-twist coupled laminated beams and plates have been investigated in recent years, yet can be cumbersome to implement quickly and efficiently. In early stages of design, employing a stiffness method approach to predict the load-deformation relationship and structural natural frequencies can be more efficient than developing a shell finite element model for each design iteration. A weak-form approach to the development of an accurate bend-twist coupled composite laminate beam element is presented herein. Comparisons are made between the stiffness matrix terms using the presented method and a shell finite element model of an idealized beam; the proposed method shows good agreement for a suite of beams with varying degrees of bend-twist coupling. The method is then extended to the calculation of natural frequencies by combining the new stiffness matrix with a corresponding consistent mass Email address: [email protected] (Pavel Babuska) Preprint submitted to Elsevier

January 10, 2018

matrix to formulate the eigenvalue problem and shows very good agreement with both analytical and shell element solutions. Keywords: element stiffness matrix, stiffness method, weak form, bend-twist coupling, composite turbine blades 1. Introduction The classical stiffness matrix method is a formulation of the equilibrium equations of ‘stick-type’ structures (e.g. beams, frames, trusses). While it may be formulated using variational principles with the finite element method (FEM), the resulting equations require no reference to general FEM theory. For beam-type elements, the deformations between nodes can usually be solved exactly, allowing for the formulation of the stiffness matrix directly (usually called the direct stiffness method) without reliance on more typical variational approaches of FEM. These advantages have made the stiffness matrix method the standard approach of structural analysis [1]. The application of the stiffness method to the analysis of laminated composite structures has been limited, however, as the anisotropic laminate response is often more complex than that of an isotropic, homogeneous material. The structural response of steel beam elements, for instance, can be solved using the Euler-Bernoulli beam equation, whereas this is more difficult for anisotropic composite beam elements. In the design phase of such laminates, where structural geometry or laminate schedule may be unknown, it may be useful to employ a stiffness matrix approach rather than build a new finite element shell model for each iteration. In [2], the authors derive the nonlinear equations of motion for elastic

2

coupled bending and torsion of homogeneous, isotropic beams. The equations are valid for beams whose coupling is a result of non-coincident mass and elastic axes, as well as nonuniform mass and stiffness section properties. The solutions of the presented equations of motion could be achieved via Galerkin’s method, however, composite beams with inherent bend-twist coupling due to material layup are not discussed. The authors of [3] investigated the effects of unbalanced laminate-induced bend-twist coupling on the vibration of cantilever composite plates. Also presented were the requisite terms to be included in a Rayleigh-Ritz modal analysis to accurately predict the first three vibrational modes. The Rayleigh-Ritz mode shapes were assumed as trigonometric sine and cosine functions and were not solved exactly, but rather presented in the context of the minimum set of contributing mode shapes to accurately capture the first three natural frequencies of a vibrating cantilevered plate with flexural-torsional coupling. No static investigation was performed, and the interaction between the kinematics and constitutive behavior observed in this paper was not discussed. In [4], a Vlasov-type linear theory solution for warping composite open-section beams is presented and compared to the results obtained by experiments on autoclave-fabricated beams. General solutions are presented for symmetric I-beams under bending and torsional loads, however, the laminate coupling terms (D16 , D26 of the ABD matrix) for beams whose flanges are unbalanced are not included in the stiffness term contribution. The correlation between experimental results and predicted solutions was satisfactory for extension-twist coupled beams but notably less satisfactory for bend-twist coupled beams. Bend-twist coupling as related to aeroelasticity and divergence instability was explored in [5]

3

utilizing a thin-walled anistropic composite beam model over classical platebeam or solid-beam models. Extensive investigation into dynamic stiffness matrices of bend-twist coupled structures was conducted by Banerjee, beginning with [6]. In [6], a dynamic stiffness matrix is developed for general bend-twist coupled Euler beams (plane sections remain plane), but Banerjee also extended the development of dynamic stiffness matrices for Timoshenko beams (shear-deformable), axially loaded beams, cross section warping, and various boundary conditions in succeeding works. The governing equations for each of the above cases and boundary conditions must be exactly solved to generate the dynamic element stiffness matrix. Additionally, the results are derived in the frequency domain to calculate the structure’s amplitudefrequency curve, but do not consider static response. A general theory for the development of exact stiffness matrices of shear-deformable thin-walled beams with closed and open cross sections that may contain coupled effects is derived in [7]. The work in [7] addresses coupling induced by offset shear and elastic centers that arise from open cross sections but does not deal with material coupling. Work in [8] follows a general method for the derivation of an element stiffness matrix for a bend-twist coupled laminate beam based on the weak form, but does not include an important coupling term and is inaccurate for highly-coupled laminates where the twisting response due to bending is large. Furthermore, the coupling stiffness parameters in the stiffness matrix are not formulated or presented analytically in [8], but rather generated by a multi-physics software developed by the National Renewable Energy Laboratory (NREL) called FAST which is used to simulate the coupled dynamic response of wind turbines. The bending, torsional, and

4

coupling stiffness parameters for flat, laminated plates are well known (see [9]), but have not been implemented in a stiffness method approach. In the current work, an accurate and efficient stiffness matrix formulation of static and dynamic response of bend-twist coupled composite laminate beams is developed. The work presented here introduces new, analyticallyderived shape functions which provide insight into the fundamental mechanics of composite beams which include bend-twist coupling. The resulting stiffness matrix allows for exact static force-deformation calculations using a single element and conceptually simple (although requiring additional discretization) natural frequency calculations. The goal was not to develop a method that is more accurate in natural frequency prediction than existing methods, but rather present a method which retains accuracy alongside its simplicity. As such, the proposed method is compared to an ABAQUS commercial FEM solution as well as one of the existing analytical methods as validation, but a broad comparison against all current methods was not conducted. 2. Derivation 2.1. Constitutive Law For thin, flat, composite laminated beams that have coincident shear, flexural, and elastic centers and are loaded at the midline, the bend-twist coupling is a result of the laminate layup schedule. A weak-form approach with assumptions from classical laminated plate theory (CLPT) can be taken to develop an element stiffness matrix for these coupled laminates if the shape functions are appropriately defined for the coupled structural response. Con5

sider a laminated beam with the coordinate system shown in Figure 1, where α corresponds to the fiber angle of each ply. Nodal coordinate notation follows the right-hand rule, consistent with the approach presented in [1]. This paper will focus on stacked single-orientation laminates where the fiber angle remains the same for each ply as shown schematically in Figure 1. This type of laminate results in the highest degree of bend-twist coupling and serves as a useful baseline for exploring the development of this method. However, the formulation presented herein is easily applicable to all symmetric laminates with any amount of bend-twist coupling, including quasi-isotropic and general unbalanced laminates. The properties of composite laminate plates are conventionally expressed in terms of an ABD matrix which relates three in-plane forces and three bending moments with three in-plane strains and three curvatures at the mid-plane. The A, B, and D terms refer to the corresponding 3x3 submatrices that make up the 6x6 ABD matrix. Any laminate which has a non-zero D16 term in its ABD matrix can benefit from 𝐹𝑧1 , 𝑤𝑧1

𝑀𝑦1 , θ𝑦1 𝐹𝑦1 , 𝑤𝑦1

𝑀𝑧1 , θ𝑧1

𝒏𝒐𝒅𝒆 𝟏

𝑀𝑥1 , θ𝑥1

𝐹𝑥1 , 𝑤𝑥1

𝛼

𝑡

𝑧

𝐹𝑧2 , 𝑤𝑧2 𝑀𝑧2 , θ𝑧2

𝑀𝑦2 , θ𝑦2

𝑙

𝐹𝑦2 , 𝑤𝑦2 𝑦

𝑏

𝑥

𝒏𝒐𝒅𝒆 𝟐

𝑀𝑥2 , θ𝑥2 𝐹𝑥2 , 𝑤𝑥2

Figure 1: Fiber angle offset and coordinate system for unidirectional laminate

6

the formulation presented in this work (see [10] for composite mechanics, ABD derivation, and CLPT). Because this work focuses on beam elements rather than plates, the derivation instead begins with a general three-dimensional beam theory approach which simplifies the results as in [11]. Equilibrium and kinematics are established in [11] which result in the load-strain relation of F = C ,

(1)

where C is the constitutive relation matrix, F is a vector containing the inplane forces, and  is a vector containing mid-plane strains. For a bend-twist coupled beam, the previous equation can be written as      Fx Ex A 0 0 0 w0      x      Mx   0 GJ K 0   θx0   =  ,      My   0 K Ex Iy 0   θy0       Mz 0 0 0 Ex Iz θz0

(2)

where Ex A represents the axial rigidity, GJ represents the St. Venant torsional rigidity, Ex Iy represents the weak axis bending rigidity, Ex Iz represents the strong axis bending rigidity and K represents the bending-torsion coupling rigidity. The derivatives in the strain-vector, , are with respect to the element length, or ‘x’, direction (as per Figure 1). This constitutive relation differs from that of an uncoupled structural element whose C matrix has zero-valued off-diagonal elements. The constitutive rigidities can be

7

represented by the following expressions per [9] or [12]: Ex A = Ex b t,   2 D26 GJ = 4 b D66 − , D22   D12 D26 K = 2 b D16 − , D22   2 D12 Ex Iy = b D11 − , D22 b3 t Ex Iz = Ex , 12 where Ex is the approximate x-direction elastic modulus, the Dij terms are components of the laminate ABD matrix, and b is the beam or plate width (Figure 1). The Dij terms relate bending moments and plate curvature. The dominant coupling parameter in K is D16 , which relates My and θx0 , or the relationship between beam bending moment and corresponding twisting rate. For uncoupled laminates, this value is zero; it also tends to be close to zero for laminates with small amounts of bend-twist coupling. 2.2. Weak-Form Formulation A 3D beam element can now be defined with 2 nodes and 6 degrees of freedom per node. The displacement vector takes the form h iT d = wx1 , wy1 , wz1 , θx1 , θy1 , θz1 , wx2 , wy2 , wz2 , θx2 , θy2 , θz2 ,

(3)

where w and θ correspond to the nodal displacements and rotations for the directions denoted in the subscript. Strains from Equation (1) can be expressed via the strain-displacement relation as  = B d, 8

(4)

where B is a matrix which contains functions relating strains and nodal displacements; these relations will be discussed in the following subsection. The conventional weak-form derivation based on potential energy minimization results in the following expression relating loads and displacements via an element stiffness matrix. This can be expressed as Z L  T P= B C B dx d,

(5)

0

or more simply P = K d,

(6)

where P is the vector of nodal forces and K is the resulting element stiffness matrix. 2.3. Strain-Displacement Matrix Derivation To this point, this method is consistent with the derivation of an uncoupled beam element. For a bend-twist coupled beam element, however, the strain-displacement relation matrix, B, is more complicated. Uncoupled beam relations denote twisting strains only due to applied torques at the nodes. As a point of comparison, in [8], while bend-twist coupling is included in the constitutive law (i.e. the K term is present as in (2)), the strain-displacement relationships are derived from classic isotropic beam elements. This approach leads to inaccuracy for strongly coupled laminates where twisting strain is strongly influenced by transverse displacements and beam rotations, the importance of which cannot be ignored in the straindisplacement matrix. The effects of out-of-plane nodal beam displacements and nodal beam rotations, wz and θy , respectively, on beam twist must also be accounted for to completely establish strain-displacement relations. 9

In order to fully populate the strain-displacement matrix for this type of laminate, the shape functions corresponding to unit element displacements and rotations must be known. Equations (7-12) are the assumed forms of the six deformation fields (three translations, and three rotations). Using the standard discretization approach, each deformation field is approximated using the nodal degrees of freedom and their respective shape functions. Outside of the standard linear shape functions for axial deformations (NA ) and the four cubic bending (NB ) and rotation (NR ) shape functions wellknown to beam theory, shape functions for beam twist due to transverse nodal displacements and beam twist due to nodal rotations must be derived. Midspan twist induced by transverse displacements and bending rotations at the nodes can be observed in experiments and FEM simulations, indicating that new shape functions must be added. The requirement for new shape functions can be understood by examining the displacement equations for each degree of freedom, noting that the θx equation contains additional expressions to fully quantify the twisting response. wx = NA1 wx1 + NA2 wx2

(7)

wy = NB1 wy1 + NR1 θz1 + NB2 wy2 + NR2 θz2

(8)

wz = NB1 wz1 + NR1 θy1 + NB2 wz2 + NR2 θy2

(9)

θx = NT 1 θx1 + Nθx z1 wz1 + Nθx θy1 θy1 + (10) NT 2 θx2 + Nθx z2 wz2 + Nθx θy2 θy2 θy = Nθy θy1 θy1 + Nθy z1 wz1 + Nθy θy2 θy2 + Nθy z2 wz2

(11)

θz = Nθz θz1 θz1 + Nθz y1 wy1 + Nθz θz2 θz2 + Nθz y2 wy2

(12)

10

The boxed expressions in Equation (10) identify the additional shape functions and nodal displacements that contribute to the θx twist equation. The naming convention on the shape functions is as follows. Nθx z1 refers to the contribution to θx (twist) by a unit displacement, wz , at node 1. Following suit, Nθx θy2 indicates the contribution to θx (twist) by a unit rotation, θy , at node 2. All shape functions are functions of x through the element length, yet will be denoted simply by N for brevity. These additional shape functions can be solved for explicitly in the following fashion. To solve for the θx relationship due to an imposed nodal displacement, wz , we have (from equilibrium or Equation (2)), Mx = GJ θx0 + K θy0 ,

(13)

My = K θx0 + Ex Iy θy0 .

(14)

With Euler-Bernoulli kinematics, wz0 = −θy , Equations (13) and (14) become Mx = GJ θx0 − K wz00 ,

(15)

My = K θx0 − Ex Iy wz00 .

(16)

As there are no applied loads in-span, My is a linear function of position and thus My00 = 0. Equation (16) can be differentiated twice to give 0 = K θx000 − Ex Iy wz0000 .

(17)

Likewise, the torque is constant within the element, meaning that the derivative of Equation (13) is also equal to zero. Combining this and Equation (17) yields two uncoupled differential equations,   K2 Ex Iy − wz0000 = 0, GJ 11

(18)



K2 Ex Iy − GJ



θy000 = 0,

(19)

which reduce to the standard governing beam equations wz0000 = 0,

(20)

θy000 = 0.

(21)

The solution form for transverse displacement, wz , with the nodal deformations representing the boundary conditions yields wz (x) = wz1 N1 (x) + θy1 N2 (x) + wz2 N3 (x) + θy2 N4 (x).

(22)

In this example derivation, imposing only a displacement at node 1 results in a simplified displacement equation, wz (x) = wz1 N1 (x),

(23)

with N1 (x), in this case, equal to the standard beam bending shape function, N1 (x) = 1 − 3 ( Lx )2 + 2 ( Lx )3 . Equation (23) can be combined with the derivative of Equation (13) (because Mx0 = 0) to give GJ θx00 − K wz1 N1000 (x) = 0.

(24)

Evaluating N1000 (x), twice-integrating θx00 , and applying zero θx boundary conditions gives the resulting shape function for θx due to a displacement of wz1 θx =


6K x2 − Lx wz1 = Nθx z1 wz1 . 3 GJ L

(25)

A similar approach can be taken to yield the shape function for θx due to a rotation, θy (noting that the boundary conditions must be altered accordingly), of θx =


3K Lx − x2 θy1 = Nθx θy1 θy1 . 2 GJ L 12

(26)

The additional bend-twist shape functions come directly from the governing differential equations in the same way as the traditional isotropic EulerBernoulli shape functions. One difference, however, is that they include material properties and thus are not purely kinematic relationships. The magnitude of the deformation in the extra relations is conditional on the relationship between coupling and torsional rigidities, not on kinematics alone. The shape functions, N , representing the element displacement field for each degree of freedom are listed in Table 1. Having established the appropriate shape function relations for coupled laminated beams, the complete strain-displacement matrix, B, can now be written using the first derivatives with respect to x of the appropriate shape functions, and the strains can now be more completely expressed as i = Bij (x) dj =

d Nij (x) dj , dx

where 1

NA1

1−

NA2

x L

x L

0.5 0 0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

1 0.5 0

1

NT 1

1−

NT 2

x L

NB1

1−3

x2 L2

NB2

x2 L2

x L

0.5 0

1 0.5 0

1

+2

x3 L3

0.5 0

1

3

−2

x3 L3

0.5 0

13

(27)

0.2

2

NR1

L ( Lx − 2 Lx 2 +

NR2

L ( Lx 3 −

x3 ) L3

0.1 0 0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

0

3

x2 ) L2

-0.1 -0.2

0

Nθx z1

6K GJ L3

(x2 − Lx)

-1 -2

2

Nθx z2

6K − GJ (x2 − Lx) L3

1 0

1

Nθx θy1

3K GJ L2

(Lx − x2 )

0.5 0

1

Nθx θy2

3K − GJ (Lx − x2 ) L2

Nθy z1

6 Lx2 − 6 Lx 3

Nθy z2

−6 Lx2 + 6 Lx 3

Nθy θy1

1 − 4 Lx + 3 Lx 2

Nθy θy2

3 Lx 2 − 2 Lx

Nθz y1

−6 Lx2 + 6 Lx 3

Nθz y2

6 Lx2 − 6 Lx 3

Nθz θz1

1 − 4 Lx + 3 Lx 2

Nθz θz2

3 Lx 2 − 2 Lx

0.5 0

2

2

1 0

0

2

-1 -2

1

2

0.5 0

1

2

0.5 0

0

2

-1 -2

2

2

1 0

1

2

0.5 0

1

2

0.5 0

Table 1: Shape functions table

14

The complete, coupled B matrix that is used in the analysis presented here is expressed as the following:   0 0 N 0 0 0 0 0 NA2 0 0 0 0 0  A1    0 Nθ0 x z1 NT0 1 Nθ0 x θy1 0 0 0 Nθ0 x z2 NT0 2 Nθ0 x θy2 0   0  , B=  0 0 0  0  0 N 0 N 0 0 0 N 0 N 0 θ θ y y2 θ z θ θ θ z y y y 1 y1 2   0 0 0 0 0 Nθz y1 0 0 0 Nθz θz1 0 Nθz y2 0 0 0 Nθz θz2 where the boxed terms are the strain-displacement relations that arise due to the newly-derived shape functions. In a conventional, uncoupled, approach, these boxed terms would be zero. Note that Equation (2) only requires the four strains, wx0 , θx0 , θy0 , θz0 , which are represented by the respective rows in B. 15

2.4. Element Stiffness Matrix Calculation The integration from Equation (5) can now be performed with the constitutive relation from Equation (1) to yield the 12x12 element stiffness matrix for the coupled beam. The boxed Ki expressions in the stiffness matrix are expressed in Equations (28-31) and are a result of the inclusion of the extra strain-displacement relations in B. Without the extra relations, the integration would simply yield the standard uncoupled stiffness matrix expressions of the form

12 Ex Iy 6 Ex Iy 4 Ex Iy , L2 , L , L3

and

2 Ex Iy . L

The new stiffness expressions are

the following: 12 (−K 2 + GJ Ex Iy ) , GJ L3 6 (K 2 − GJ Ex Iy ) K2 = , GJ L2 −3 K 2 + 4 GJ Ex Iy K3 = , GJ L −3 K 2 + 2 GJ Ex Iy K4 = . GJ L

K1 =

16

(28) (29) (30) (31)

The fully coupled, symmetric element stiffness matrix, K, used in the analysis presented here is expressed as the following:  Ex A L

0

17

  12 Ex Iz  0 L3    0 0    0 0    0 0   6 Ex Iz  0 L2  K= − ELx A 0    0 − 12 EL3x Iz    0 0    0 0    0 0  6 Ex Iz 0 L2

0

0

0

0

− ELx A

0

0

− 12 EL3x Iz

0

0

0

0

0

0

0

0

0

6 Ex Iz L2

K1

0

K2

0

0

0

−K1

0

K2

0

GJ L

K L

0

0

0

0

− GJ L

−K L

K2

K L

K3

0

0

0

−K2

−K L

K4

0

0

0

4 Ex Iz L

0

− 6 ELx2 Iz

0

0

0

0

0

0

0

Ex A L

0

0

0

0

0

0

0

− 6 ELx2 Iz

0

12 Ex Iz L3

0

0

0

−K1

0

−K2

0

0

0

K1

0

−K2

0

− GJ L

−K L

0

0

0

0

GJ L

K L

K2

−K L

K4

0

0

0

−K2

K L

K3

0

0

0

2 Ex Iz L

0

− 6 ELx2 Iz

0

0

0

0



     0    0    0   2 Ex Iz   L .  0   6 Ex Iz  − L2    0    0    0   6 Ex Iz L2

4 Ex Iz L

3. Validation 3.1. Idealized Beam Comparison The performance of this stiffness method approach to the analysis of bending-torsion coupled beams was evaluated by comparison to a linear elastic ABAQUS/Standard [13] FEM of an idealized test beam. The idealized beam is 2.54 meters long, 0.0254 meters wide, and 0.001048 meters thick. The beam has high length-to-width and width-to-thickness ratios to ensure that the Euler-Bernoulli assumptions are valid as well as to minimize the effects of localized through-width or through-thickness deformation concentrations that arise due to the imposition of displacements at individual nodes. The material used in the analyses was Hexcel IM7/8552 unidirectional carbon fiber composite prepreg. Material properties for the IM7 are as follows: E1 = 164 GPa, E2 = 8.96 GPa, G12 = 4.93 GPa, and ν12 = 0.316. A total of 1696 four-node, reduced integration shell elements (S4R) were used in the FEM mesh. The beam ends were kinematically constrained to a single point to represent the respective points on which boundary conditions could be applied rather than on nodes contained within the mesh, as this would impart unwanted localized stresses and deformation. A unit displacement at node 2 was imposed in each degree of freedom and the respective reaction forces at node 1 were recorded by the ABAQUS model. Evaluation of the stiffness matrix was done by direct comparison of the top-right quadrant to the finite element model reaction forces. Because unit magnitude boundary condition displacements were imposed on node 2 in the FE model, the top-right quadrant represents the reaction forces at node 1 to displacements at node 2. 18

3.2. ABAQUS vs Stiffness Matrix Values Comparison The shell FEA simulations were conducted for fiber angles ranging from 0◦ to 90◦ and each set of reaction forces was recorded. In the top-right quadrant, the five terms that govern the beam bending and torsional response are the following: −K1 , K2 , K4 , − GJ , and − K . Both coupled and uncouL L pled analytical stiffness matrix expressions are plotted for all fiber angles and compared to values obtained from the commercial FEM solution. The comparisons can be seen in Figures 2 - 6. It should be noted that the − GJ L and − K terms are identical for the coupled and uncoupled cases. This is L because the additional shape functions for the fully-coupled element do not contribute to those respective stiffnesses. Note that the uncoupled analytical solution greatly over-predicts stiffness for fiber angles ranging from near 0◦ to near 50◦ . This is because it doesn’t

0.3

Abaqus Model Uncoupled Coupled

K1 [N=m]

0.25

0.2

0.15

0.1

0.05

0 0

10

20

30

40

50

60

70

80

90

,

Figure 2: K1 , translational stiffness for an imposed transverse nodal displacement

19

0.4

Abaqus Model Uncoupled Coupled

0.35

K2 [N=m]

0.3 0.25 0.2 0.15 0.1 0.05 0 0

10

20

30

40

50

60

70

80

90

,

Figure 3: K2 , translational stiffness for an imposed nodal rotation

0.35

Abaqus Model Uncoupled Coupled

0.3

K4 [N=m]

0.25 0.2 0.15 0.1 0.05 0 -0.05 0

10

20

30

40

50

60

70

80

90

,

Figure 4: K4 , rotational stiffness for an imposed nodal rotation

20

0.07

Abaqus Model Coupled

0.06

[N=m]

0.05

GJ L

0.04

0.03

0.02

0.01 0

10

20

30

40

50

60

70

80

90

,

Figure 5:

GJ L ,

torsional stiffness for an imposed nodal twist

0.07

Abaqus Model Coupled

0.06

K L

[N=m]

0.05 0.04 0.03 0.02 0.01 0 -0.01 0

10

20

30

40

50

60

70

80

90

,

Figure 6:

K L,

rotational stiffness for an imposed nodal twist

21

allow for the element twisting, which leads to artificial stiffening through constraint. The coupled solution is in fact exact as all the shape functions are solved explicitly from the governing equations. The deviation from the FEM solutions is possibly a result of how the displacements were applied or how the reaction force was recovered in the finite element model. 4. Sensitivity to the Coupling Parameter, K It is important to the engineer to understand the role that bend-twist coupling plays in the structure and when the inclusion of the coupled relations are necessary to the analysis. Observing how the uncoupled and coupled stiffness terms vary in Equations (28-31), the difference is seen to be conditional on the ratio of

K2 . Ex Iy GJ

The sensitivity ratio can be calculated directly

from CLPT as all of the parameters are functions of D terms from the ABD matrix as the following: 

Λ= 

D12 D26 D22

2

D16 −  D2 D11 − D12 D66 − 22

2 D26 D22

.

(32)

This ratio, Λ, serves as a nondimensional parameter that dictates how sensitive the laminate is to the coupling terms in the stiffness matrix. As shown in Figure 7, the ratio is conditional only on fiber angle and is significant for a large suite of unidirectional laminate fiber angles. For quasi-isotropic laminates or generally arbitrary layups, the sensitivity ratio can be calculated as a metric to aid in the design and analysis of coupled laminates. A ratio near zero suggests that uncoupled strain-displacement relations should be sufficiently accurate in analysis. However, the simplicity of the exact formulation leaves little need to avoid the completely coupled stiffness matrix. 22

0.7 0.6 0.5

$

0.4 0.3 0.2 0.1 0 0

10

20

30

40

50

60

70

80

90

,

Figure 7: Sensitivity of laminate fiber angle to the difference between coupled and uncoupled stiffness matrix terms

5. Applications One application already investigated is the use of the method for static analyses. The structural engineering community is very familiar with stiffness matrix assembly, making this method immediately available for their use. By utilizing discretization, the method is readily extendible for use with in-span loads. Additionally, through the use of rotation matrices, it is suitable for 3D frame analysis. The analysis could also be adapted into many existing nonlinear co-rotational methods for large deformations. This section will focus instead on dynamic analysis, as the stiffness matrix is immediately suitable for this application. The usefulness of the element stiffness matrix can be extended into the structural dynamics regime by formulating an element mass matrix.

23

5.1. Element Mass Matrix Mass distribution for use in the stiffness matrix method can be established in a number of ways, including direct mass lumping or variational mass lumping. Following the variational mass lumping approach establishes the kinetic energy as, T =

 i 1 h  ˙2 ρ A wx + w˙y2 + w˙z2 + I0 θ˙x2 , 2

(33)

where the dotted terms are the respective degree of freedom velocities, ρ is the material density, A is the cross-sectional area, and I0 = Iy +Iz is the polar moment of inertia for the twisting body by the perpendicular axis theorem. The mass matrix can then be represented by the following integral, Z L M = ρA (NA )T NA dx 0 Z L +ρ A (NB )T NB dx 0 Z L +ρ I0 (NT )T NT dx,

(34)

0

where NA , NB , and NT are shape function matrices containing the relevant shape functions for axial, bending, and torsion terms, respectively. For a consistent mass matrix, shape functions for the mass matrix are the same as those for the stiffness matrix. However, only the shape functions corresponding to displacement equations for wx , wy , wz , and θx in Equations (7-10) are needed due to negligible θy and θz rotational inertia. Thus, NA would contain the shape functions present in the first row of the shape function matrix in (35). Similarly, NB corresponds to the shape functions included in the second and third rows of (35), while NT contains shape functions present 24

in the fourth row of (35). The resulting shape function matrix is 4x12 in dimension and results in a 12x12 consistent mass matrix.

25

The shape function matrix, N, used to derive the consistent mass matrix, M, is expressed as the following:   N 0 0 0 0 0 NA2 0 0 0 0 0  A1     0 NB1 0 0 0 NR1 0 NB2 0 0 0 NR2  , N= (35)    0 0 NB1 0 NR1 0 0 0 NB2 0 NR2 0    0 0 Nθx z1 NT 1 Nθx θy1 0 0 0 Nθx z2 NT 2 Nθx θy2 0 where the boxed terms are the newly-derived shape functions. In a conventional, uncoupled, mass matrix derivation, these positions in the matrix would be zero. The new expressions for coupled mass terms are as follows: 26 M1 =

13 A L ρ 6 I0 K 2 ρ + , 35 5 GJ 2 L

11 A L2 ρ 3 I0 K 2 ρ + , 210 5 GJ 2 11 A L2 ρ 3 I0 K 2 ρ M3 = − + , 210 5 GJ 2 A L3 ρ 3 I0 K 2 L ρ M4 = + , 105 10 GJ 2 9 A L ρ 6 I0 K 2 ρ M5 = − , 70 5 GJ 2 L 13 A L2 ρ 3 I0 K 2 ρ M6 = − − , 420 5 GJ 2 M2 =

(36) (37) (38) (39) (40) (41)

13 A L2 ρ 3 I0 K 2 ρ M7 = − , 420 5 GJ 2 A L3 ρ 3 I0 K 2 L ρ M8 = − − . 140 10 GJ 2

(42) (43)

The entire mass matrix can be presented as the following, where the boxed terms represent the terms

27

expressed by Equations (36-43).  ALρ 0 0 0 0 0  3  13 A L ρ 11 A L2 ρ  0 0 0 0 35 210   I0 K ρ  0 0 M1 − 2 GJ M2 0   Kρ I0 L ρ  0 0 − I20 GJ − I04KGJL ρ 0 3    0 0 M2 − I04KGJL ρ M4 0  2  11 A L ρ A L3 ρ  0 0 0 0 210 105 M= ALρ  6 0 0 0 0 0   9ALρ 13 A L2 ρ  0 0 0 0 70 420   I0 K ρ  0 0 M5 M7 0 2 GJ   Kρ I0 L ρ  0 0 − I20 GJ − I04KGJL ρ 0 6   I0 K L ρ  0 0 M6 M8 0 4 GJ  2 L3 ρ 0 − 13 A420L ρ 0 0 0 − A140

ALρ 6

0

0

9ALρ 70

0

0

0

0

0

M5

Kρ − I20 GJ

M6

0

0

I0 K ρ 2 GJ

I0 L ρ 6

I0 K L ρ 4 GJ

0

0

M7

I0 K L ρ 4 GJ

M8

420

0

0

0

ALρ 3

0

0

0

0

0

13 A L ρ 35

0

0

0

0

0

M1

I0 K ρ 2 GJ

M3

0

0

I0 K ρ 2 GJ

I0 L ρ 3

I0 K L ρ 4 GJ

0

0

M3

I0 K L ρ 4 GJ

M4

0

0

0

0

0

13 A L2

− 11 A210L

0

ρ

2

ρ

0

0

0

 

2  − 13 A420L ρ  

0

     0    0  3 AL ρ  − 140     0  11 A L2 ρ  − 210     0    0    0  A L3 ρ 105

5.2. Natural Frequency Evaluation Having established element mass and stiffness matrices, the system natural frequencies can be solved as a multiple degree of freedom linear dynamics eigenvalue problem. Solving the eigenvalue problem in Equation (44) yields the structural circular frequencies, ωn . |Kel − ωn2 Mel | = 0

(44)

  Kel − ωn2 Mel φn = 0

(45)

The mode shapes are the corresponding eigenvectors, φn , that satisfy Equation (45). To date, bend-twist coupling has only been applied in cantilevered settings such as wings or turbine blades, thus a cantilever beam subject was implemented to evaluate the performance of the proposed stiffness and mass matrices in calculating structural natural frequencies. The stiffness matrix method model was compared to an ABAQUS/Standard finite element model, as well as to the analytical solution to the coupled equations of motion as presented by Banerjee in [6] and adapted for bend-twist composites by Kramer in [12]. For the remainder of this work, the term “coupled analytical solution” denotes the solution presented in [12]. Natural frequencies in the ABAQUS FEM were calculated via linear perturbation, utilizing the Lanczos method for the eigensolution, standard to the software program. The coupled analytical solution was implemented primarily in Mathematica due to its powerful abilities of symbolic computation. The stiffness matrix method of this paper, referred to hereafter as “SMM,” was evaluated for 1 element, 2 elements, 5 elements, and 20 elements in MATLAB. 28

The natural frequencies were non-dimensionalized as per [9] and [12] to allow for cross comparison between geometries and material properties. The non-dimensionalization is done for ωn as r Ωn = ωn L2

ρt , D0

where Ωn is the non-dimensionalized natural frequency and D0 =

E1 t3 12 (1−ν12 ν21 )

corresponds to the reference bending stiffness at α = 0◦ . This enables the results presented in this work to be used as a baseline performance metric for others to check their own implementations. In order to further demonstrate the need to include coupling effects in the analysis of these structures, the comparison between beam first natural frequencies found via the commercial FEM solution, the coupled analytical solution, and the uncoupled analytical solution is shown in Figure 8 for cantilever boundary conditions. The uncoupled analytical solution incorporates the effective bending rigidity, Ex Iy , but has no coupling in the equations q of motion, thus resulting in the expression ωn = αn2 ρEAx LIy4 , where αn are the constants relating to the mode of vibration for an isotropic, homogeneous cantilever, 1.875, 4.694, 7.885, etc, as presented in [14]. The uncoupled analytical solution clearly overestimates the first natural frequency as it artificially constrains twisting in the bending modes, and vice-versa, thereby increasing stiffness and frequency. Comparing the coupled analytical solution to the commercial FEM solution in Figure 8 shows reasonable agreement, though the ABAQUS FEM slightly overpredicted Ω values at larger fiber angles compared to the coupled analytical solution. This was also the case in [12]. The coupled analytical solution exactly solves the coupled equations of motion and thus could likely be assumed to be more accurate. 29

5 FEM Coupled Analytical Uncoupled Analytical

4

+

3 2 1 0

0

10

20

30

, Figure 8: First natural frequency comparisons for FEM, uncoupled, and coupled analytical models of cantilever

To appropriately compare the SMM model, it is important to observe the convergence with the increase in number of elements. More elements give the SMM model greater accuracy which converges to the coupled analytical solution. The results of first bending natural frequency for the SMM model of various element meshes is presented in Figure 9. Convergence is similar for higher bending modes and torsional modes. Results comparing the 20-element SMM model and the coupled analytical solution are presented in Figure 10 where very good agreement is shown for first and second bending natural frequencies as well as the first twisting natural frequency. ABAQUS model mode shapes for a 10◦ fiber angle beam are presented alongside the figures to give clarity to the reader about how the coupling affects the deformation state. A secondary comparison of

30

1 element 2 elements 5 elements 20 elements

4

+

3 2 1 0

0

10

20

30

, Figure 9: Mesh refinement resulting in 1st bending natural frequency convergence

mode shapes, generated by the SMM model for 0◦ and 10◦ fiber angle beams, is presented in Figure 11 illustrating the difference between the uncoupled and coupled vibrational modes and that the coupled mode matches the results from the ABAQUS solution. Extensive literature exists regarding mode shapes of bend-twist-coupled beams and plates both in air and in water. A more thorough presentation of coupled mode shapes can be found in [12]. Convergence of the stiffness matrix method values to the coupled analytical solution is shown to be within 2% with as few as five elements and within 0.1% with as few as 20 elements. A tabular-form comparison of the nondimensionalized natural frequencies between analysis methods is presented in Table 2 as well. Note that the 20-element SMM agrees very well with the coupled analytical solution for the presented modes while the ABAQUS solution tends to predict slightly higher frequencies.

31

SMM 20 Element Coupled Analytical

4

+

3 2 1 0

0

10

20

30

,

(b) 10◦ beam 1st bending mode shape

(a) 1st bending frequency 25 SMM 20 Element Coupled Analytical

20

+

15 10 5 0

0

10

20

30

,

(d) 10◦ beam 2nd bending mode shape

(c) 2nd bending frequency

SMM 20 Element Coupled Analytical

40

+

30

20

10 0

10

20

30

,

(f) 10◦ beam 1st twisting mode shape

(e) 1st torsional frequency

Figure 10: Natural frequency comparisons for SMM 20 element and coupled analytical models presented alongside representative bend-twist coupled ABAQUS model mode shapes

32

0.04

0.04

0.02

0.02

0

0

-0.02

-0.02

-0.04

-0.04

-0.06 0.1

-0.06 0.1 0.2

0.2

0.15

0

0.15

0

0.1

0.1

0.05 -0.1

0.05 -0.1

0

(a) 0◦ 1st bending

0

(b) 10◦ 1st bending

0.04

0.04

0.02

0.02

0

0

-0.02

-0.02

-0.04

-0.04

-0.06 0.1

-0.06 0.1 0.2

0.2

0.15

0

0.15

0

0.1

0.1

0.05 -0.1

0.05 -0.1

0

(c) 0◦ 2nd bending

0

(d) 10◦ 2nd bending

0.04

0.04

0.02

0.02

0

0

-0.02

-0.02

-0.04

-0.04

-0.06 0.1

-0.06 0.1 0.2

0.2

0.15

0

0.15

0

0.1

0.1

0.05 -0.1

0.05 -0.1

0

(e) 0◦ 1st twisting

0

(f) 10◦ 1st twisting

Figure 11: Comparison of mode shapes for coupled and uncoupled beams

33

Ω 1st flexural

2nd flexural

1st twisting

ABAQUS

3.51

21.97

20.01

SMM (20 elements)

3.50

22.02

18.34

Coupled Analytical

3.51

21.97

18.34

ABAQUS

2.68

15.10

27.19

10 SMM (20 elements)

2.45

14.95

25.69

Coupled Analytical

2.45

15.10

25.33

ABAQUS

1.80

9.70

34.61

20 SMM (20 elements)

1.56

9.67

34.63

Coupled Analytical

1.56

9.70

34.25

ABAQUS

1.33

7.30

35.41

30 SMM (20 elements)

1.17

7.30

35.11

Coupled Analytical

1.17

7.30

35.25

α

0

Method

Table 2: Comparison of nondimensionalized natural frequencies between analysis methods for cantilever beams of varying fiber angle and level of coupling

34

6. Conclusions In this work, the theory and derivation of an element stiffness matrix for bend-twist coupled composite laminated beams was presented. The accuracy of the stiffness matrix terms was compared with those generated by an ABAQUS finite element model of an idealized beam geometry. The comparison was made by imposing unit displacements in each degree of freedom at one node and recording the reaction forces from the other, which showed very good agreement for each degree of freedom. Specifically of interest was the comparison between stiffness terms derived from an uncoupled beam versus a coupled beam. The coupled stiffnesses are reduced because the bend-twist coupled beams twist between the boundaries regardless of the degree of freedom displaced, especially for laminates with high degrees of coupling. The sensitivity of structural analysis results to the inclusion of full-coupling was presented in the form of a non-dimensional sensitivity ratio. This ratio denotes the degree of difference between stiffness values derived from uncoupled and coupled formulations. The sensitivity ratio can be calculated for any laminate and the authors recommend the inclusion of bend-twist coupling in the analysis if the ratio is not near zero. This method should be implemented for any non-zero unidirectional or non-cross-ply laminate as it is more complete and accurate for all laminates with a non-zero D16 term in their ABD matrix. Also presented are all requisite shape functions for the generation of a consistent element mass matrix for a bend-twist coupled laminate beam. For static analyses, only one element is required, with the exception of cases that include in-span loading. With the element mass and stiffness matri35

ces, the common eigenvalue problem can be solved to yield the structural natural frequencies. Non-dimensional natural frequencies derived from the stiffness matrix method, finite element analysis, and closed-form analytical solutions for both coupled and uncoupled laminates are compared showing near identical agreement between the stiffness matrix method and the analytical solution for a coupled system. Finite element solutions tend to slightly overpredict natural frequency, while the analytical solution from the uncoupled equations of motion is shown to be drastically different for non-zero fiber angle laminates. The primary impacts of this paper have been two-fold. First, a beam element which is coupled both constitutively and kinematically in mass and stiffness at the element scale was developed in order to wholly characterize the bend-twist coupled deformation response of such a laminate. Second, the application of the beam element into a stiffness matrix method approach of analysis was presented to illustrate the simplicity and accuracy of the formulation when compared against a commercial finite element software and a closed-form continuum analytical solution. Calculating the structure’s response by stiffness matrix method can be more versatile than solving the coupled equations of motion exactly while being simpler than developing a shell or continuum finite element model. Stiffness matrix methods are also wellknown and preferred by many analysts and can be a more useful approach than analytical methods for coupled systems which become complicated very quickly. Lastly, while the main focus of the work has been on development and validation of the stiffness matrix method approach of analysis for bend-twist

36

coupling, some important effects of this coupling are also discussed. Primarily, the effect of the coupling on the resulting stiffness is explored through a non-dimensional parameter in Figure 7. The coupling is shown to lead to a significant reduction in the beam stiffness. Second, the mode shapes shown in Figures 10 and 11 highlight the coupling effects in the deformation of the vibrating beam. In the 0◦ case, the bending and twisting are entirely decoupled while, in the 10◦ case, there is significant twist and deflection in the midline of the beam. Of course, as observed in other studies, this coupling affects not only the mode shapes, but also the vibration frequencies to a large degree, thus the correct predictions of natural frequencies are of great importance. 7. Acknowledgements The authors would like to acknowledge Tyler van Iderstein for his continued support in discussing and evaluating shortcomings during the development of this work. The authors would also like to acknowledge Ramona Barber for her help in reviewing and editing the content of this work. References [1] R. Hibbeler, Structural Analysis, 8th Edition, Pearson Prentice Hall, Upper Saddle River, New Jersey 07458, 2012. [2] D. Hodges, E. Dowell, Nonlinear equations of motion for the elastic bending and torsion of twisted nonuniform rotor blades, techreport, NASA (Dec. 1974).

37

[3] D. Jensen, E. Crawley, J. Dugundji, Vibration of cantilevered graphite/epoxy plates with bending-torsion coupling, Journal of Reinforced Plastics and Composites 1 (1982) 254–269. [4] R. Chandra, I. Chopra, Experimental and theoretical analysis of composite i-beams with elastic couplings, AIAA Journal 29 (12) (1991) 2197–2206. [5] L. Librescu, O. Song, On the static aeroelastic tailoring of composite aircraft swept sings modelled as thin-walled beam structures, Journal of Composites Engineering 2 (5-7) (1992) 497–512. [6] J. Banerjee, Coupled bending-torsional dynamic stiffness matrix for beam elements, International Journal for Numerical Methods in Engineering 28 (6) (1989) 1283–1298. [7] K. Moon-Young, K. N. II, Y. Hee-Taek, Exact dynamic and static stiffness matrices of shear deformable thin-walled beam-columns, Journal of Sound and Vibration 267 (2003) 29–55. [8] R. D. F. Lopez, A 3D finite beam element for the modelling of composite wind turbine wings, Masters thesis, Royal Institute of Technology (KTH) (2013). [9] T. A. Weisshaar, B. L. Foist, Vibration tailoring of advanced composite lifting surfaces, Journal of Aircraft 22 (2) (1985) 141–147. [10] J. Vinson, R. Sierakowski, The Behavior of Structures Composed of Composite Materials, Martinus Nijhoff, 1987. 38

[11] L. Andersen, S. R. Nielsen, Elastic beams in three dimensions, DCE lecture notes no. 23, Aalborg University Lecture Notes (August 2008). [12] M. R. Kramer, Z. Liu, Y. L. Young, Free vibration of cantilevered composite plates in air and in water, Journal of Composite Structures 95 (2013) 254–263. [13] SIMULIA by Dassault Systemes, ABAQUS Version 6.14 Documentation (2014). [14] S. S. Rao, Vibration of Continuous Systems, John Wiley & Sons, Inc., 2007.

39