Evaluation of strength of component-laminates in strip-based mechanisms

Composite Structures 100 (2013) 1–16

Contents lists available at SciVerse ScienceDirect

Composite Structures journal homepage: www.elsevier.com/locate/compstruct

Evaluation of strength of component-laminates in strip-based mechanisms q Hemaraju Pollayi a,⇑, Dineshkumar Harursampath b, Wenbin Yu a a

Solid Mechanics and Structures Laboratory (SMS Lab), Department of Mechanical and Aerospace Engineering, Utah State University, UMC 4130, Logan, UT 84322-4130, USA Nonlinear Multifunctional Composites Analysis and Design Laboratory (NMCAD Lab), Department of Aerospace Engineering, Indian Institute of Science, Bangalore 560012, Karnataka, India b

a r t i c l e

i n f o

Article history: Available online 11 January 2013 Keywords: Load-carrying capacity Four-bar mechanism Failure criterion Reserve factor Revolute joints VAM

a b s t r a c t This paper deals with the evaluation of the component-laminate load-carrying capacity, i.e., to calculate the loads that cause the failure of the individual layers and the component-laminate as a whole in fourbar mechanism. The component-laminate load-carrying capacity is evaluated using the Tsai–Wu–Hahn failure criterion for various lay-ups. The reserve factor of each ply in the component-laminate is calculated by using the maximum resultant force and the maximum resultant moment occurring at different time steps at the joints of the mechanism. Here, all component bars of the mechanism are made of fiber reinforced laminates and have thin rectangular cross-sections. They could, in general, be pre-twisted and/ or possess initial curvature, either by design or by defect. They are linked to each other by means of revolute joints. We restrict ourselves to linear materials with small strains within each elastic body (strip-like beam). Each component of the mechanism is modeled as a beam based on geometrically non-linear 3-D elasticity theory. The component problems are thus split into 2-D analyses of reference beam cross-sections and non-linear 1-D analyses along the three beam reference curves. For the thin rectangular cross-sections considered here, the 2-D cross-sectional nonlinearity is also overwhelming. This can be perceived from the fact that such sections constitute a limiting case between thin-walled open and closed sections, thus inviting the non-linear phenomena observed in both. The strong elastic couplings of anisotropic composite laminates complicate the model further. However, a powerful mathematical tool called the Variational Asymptotic Method (VAM) not only enables such a dimensional reduction, but also provides asymptotically correct analytical solutions to the non-linear cross-sectional analysis. Such closed-form solutions are used here in conjunction with numerical techniques for the rest of the problem to predict more quickly and accurately than would otherwise be possible. Local 3-D stress, strain and displacement fields for representative sections in the component-bars are recovered, based on the stress resultants from the 1-D global beam analysis. A numerical example is presented which illustrates the failure of each component-laminate and the mechanism as a whole. Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction This paper is concerned with the evaluation of the componentlaminate load-carrying capacity, i.e., to calculate the loads that cause the failure of the individual layers and the componentlaminate as a whole in four-bar mechanism. The only restriction in this analysis is that the strains within each elastic body (strip-like beam) remain small. This work, in other words, does not deal with materials exhibiting non-linear constitutive laws at the 3-D level.

q Presented as Paper AIAA 2012-1388 at the 53rd AIAA/ASME/ASCE/AHS/ACS Structures, Structural Dynamics, and Materials Conference, Honolulu, Hawaii, USA, April 23–26, 2012. ⇑ Corresponding author. Address: Department of Mechanical and Aerospace Engineering, Utah State University, Logan, UT 84322-4130, USA. Tel.: +1 435 757 6204; fax: +1 435 797 2417. E-mail addresses: [email protected] (H. Pollayi), dinesh@ aero.iisc.ernet.in (D. Harursampath), [email protected] (W. Yu).

0263-8223/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compstruct.2012.12.030

Failure criteria are used to assess the possibility of failure of a material. A number of phenomenological failure criteria for anisotropic materials have been developed which essentially provide a practical basis for the design. These phenomenological criteria can be roughly divided into two categories: Stress-based criteria, such as the maximum stress criterion, and Strain-based criteria, such as the maximum strain criterion. These two criteria predict physically different failure envelopes. The strain-based criteria for anisotropic materials have not received the same degree of attention as that of the stress-based criteria because of the difficulties in measuring strains to failure. Amongst the stress-based criteria, one of the most rationally developed interactive failure criteria is the stress tensor polynomial criterion of Tsai and Wu [1] which allows for the interaction of anisotropic stresses on failure. The most commonly used stress-based failure criteria i.e., the Tsai and Hahn [2] failure criterion is used in the present work for failure analysis. To solve for the complete flexible, non-linear multi-body beam system dynamics, which is a 3-D non-linear problem, one has to

2

H. Pollayi et al. / Composite Structures 100 (2013) 1–16

split the component beam problems into two separate problems each in order to solve them in simultaneously efficient and accurate manner. In the work of Berdichevsky [3], the 3-D elasticity representation of a beam was shown to give rise to two separate problems: a linear 2-D problem over the beam cross-section, which provides a set of elastic constants (elements of linear stiffness matrix) and a set of ‘‘recovering relations’’ for 3-D displacements, strains, and stresses; and a non-linear 1-D problem along the beam reference curve. The process of splitting the 3-D non-linear elasticity problem into two separate problems by taking advantage of the small parameters inherent to the problem is called dimensional reduction. A powerful mathematical tool called the Variational Asymptotic Method (VAM) has been used to formulate the complete problem for dimensionally reducible structures. Before applying the VAM to these dimensionally reducible structures, the kinematics of the problem is formulated using the procedure outlined by Danielson and Hodges [4]. Based on this approach, Cesnik and Hodges [5] solved linear cross-sectional problems for materials having generally anisotropic and inhomogeneous properties. Popescu et al. [6] studied the obliqueness effects in asymptotic cross-sectional analysis by treating the effect exactly for a prismatic beam and regarded the effect as a small parameter for the initially twisted and curved beam. Rajagopal and Hodges [7] proposed a beam theory for the in-plane deformation of an initially curved composite strip and this theory serves as a good validation tool for VABS. Later it was generalized by Yu [8] and still maintaining a cross-sectional analysis tool called Variational Asymptotical Beam Sectional Analysis (VABS) which is a 2D Finite Element Code written in Fortran 90. VABS is used to model composite beams-like structure with arbitrary geometry made with arbitrary material to obtain the cross-sectional properties such as neutral axes, centroid, center of gravity, shear center, classical beam model, Timoshenko beam model, and Vlasov beam model and recover the 3-D strain and stress fields distribution throughout the cross section. VABS [9] is obtaining reputation as the state-of-the-art cross-sectional analysis program with an increasing user community throughout the world. In the present analysis, VABS is used for comparing linear stiffness coefficients with the results from the analytical expressions for non-linear stiffness coefficients of each component-laminate in the four-bar mechanism. Cross-sectional analyses are usually linear, but there are at least a couple of exceptions. Popescu and Hodges [10] have introduced geometric nonlinearity at cross-sectional level by discretizing the cross-section of a strip-like composite beam into finite-elements. They mainly studied the trapeze effect, a non-linear effect, which is typically included in the analyses of rotating structures such as helicopter rotor blades, propellers, and turbomachinary blades because of the presence of large centrifugal forces. It leads to a change in the torsional rigidity that varies with axial force. The trapeze effect is caused by non-linear extension-twist coupling in beams undergoing large axial forces and is caused by the presence of certain non-linear terms in the strain field because of moderate local rotation [11]. They concluded that as the cross-section departs from a strip configuration, the trapeze effect becomes less and less important compared to the overall torsional rigidity. Harursampath et al. [12] studied non-linear cross-sectional problems such as the bending of thin-walled, hollow, circular tubes, which leads to a non-linear, moment–curvature relation. Harursampath [13] has developed closed-form analytical solutions for fiber reinforced laminated beams with thin rectangular crosssections. These strips could, in general, be pre-twisted and/or possess initial curvature. He compared the results with linear crosssectional models and margin between the linear and non-linear results are very large for certain problems. In the present study, we

are using these closed-form analytical solutions for the cross-sectional properties and warping of such beams. Beams can be defined as elastic bodies each of whose volume is that spanned by a cross-section translating along a smooth reference curve for distances much larger than the largest cross-sectional dimension. The equations of motion resulting from the modeling of multi-body beam systems present distinguishing features: they are stiff, non-linear, differential–algebraic equations. The ‘‘stiffness’’ of the equations stems not only from the presence of high frequencies in the component beams, but also from the infinite frequencies associated with the kinematic constraints. Bauchau and Kang [14] presented a multi-body formulation for helicopter nonlinear dynamic analysis and validated it finding a good correlation between analytical solutions and multi-body formulations. Bauchau [15] developed two computational schemes, energy preserving and energy decaying schemes for the dynamic analysis of flexible, non-linear multi-body systems and simulated realistic multi-body systems to illustrate the efficiency and accuracy of these two schemes. Bauchau and Hodges [16] studied the behavior of non-linear multi-body systems involving elastic members made of laminated, anisotropic composite materials and found some interesting elastic couplings in the system. Bauchau’s commendable and comprehensive work on this area is based entirely on numerical techniques, while the current work is partly analytical. An integrated approach for efficient and accurate analysis of composite structures by connecting high-fidelity models of composite beams with a versatile multi-body dynamic environment is introduced by Yu [17]. Pollayi and Harursampath [18] demonstrated the importance of geometrically non-linear cross-sectional analysis of certain composite beam-based four-bar mechanisms in predicting system dynamic characteristics and analysed two-different ways of arranging the beams in the mechanism, named as width-wise and thickness-wise connections, and simulated using the available commercial software (I-DEAS + NASTRAN + ADAMS) to see the practical applicability by arranging the beams in the mechanism(s) in distinctly different orientations w.r.t. each other. 1.1. Contribution of this paper The complete formulation is developed by incorporating geometric nonlinearity both at the beam cross-sectional level and the beam reference curve level. With the model, tools and procedure in place, the authors evaluated the component-laminate load-carrying capacity in the present four-bar mechanism problem using the Tsai–Wu–Hahn failure criterion. The complete analysis methodology can be viewed as a four-step procedure: First, the cross-sectional properties of each initially twisted and curved bar of the mechanism is determined analytically based on an asymptotic procedure, starting from Classical Laminated Shell Theory (CLST) and taking advantage of its thin strip geometry. Second, the dynamic response of the non-linear, flexible four-bar mechanism is simulated by treating each bar as a 1-D beam reference curve, discretized using mixed finite elements, and employing energy-preserving and-decaying time integration schemes for unconditional stability. Third, local 3-D deformations and stresses in the entire system are recovered, based on the 1-D responses predicted in the previous step and the analytical recovery relations obtained in the first step. Finally, the laminate load-carrying capacity is evaluated based on the 3-D stress states predicted in step 3. 2. Deformed state representation: beam configuration and base vectors All component beams of the flexible, multi-body system are made of generally anisotropic fiber reinforced composite laminates

3

H. Pollayi et al. / Composite Structures 100 (2013) 1–16

and have thin rectangular cross-sections. They could, in general, be pre-twisted and/or possess initial curvature either by design or by manufacturing and/or handling defects. They are linked to each other by means of joints whose behavior constitutes end constraints on the beams. Based on geometrically non-linear 3-D elasticity theory, each component of the mechanism is modeled as a beam undergoing large deflections, rotations, and non-linear warping of cross-sectional size and shape. The component problems are thus split into non-linear, yet analytical, 2-D analyses of reference beam cross-sections and non-linear FEM-based 1-D analyses along the beam reference curves. The beam(s) can be idealized as a reference curve(s) and a typical reference cross-section(s), as shown in Fig. 1. It is convenient to introduce a reference frame A, in which are fixed dextral, mutually perpendicular, unit vectors Ai for i = 1, 2, 3. The frame A is an absolute frame as far as deformation is concerned, in that the orientation of the local undeformed beam cross-section in A is a function only of x1 and not of time t. The motion of A in an inertial frame I is, however, supposed to be known for all time. This assumption is easily relaxed for applications to flexible multi-body dynamics; it is made only for the sake of simplicity. Let x1 denote arc-length along a curved refernce curve r within an undeformed, but initially curved and twisted beam. Let xa denote lengths along straight lines that are orthogonal to each other and to the reference P curve r within a cross-section ðx1 Þ. Here a point on the undeformed beam reference curve r is located relative to a point fixed in frame A by the position vector r(x1). At each point along r define a frame b in which are fixed orthogonal unit vectors bi for i = 1, 2, 3 such that b2(x1) and b3(x1) are tangent to the coordinate curves x2 and x3 at r and b1 is tangent to r. The frame b has an orientation that is fixed in A for any fixed value of x1 but varies along the beam if the beam is initially curved or twisted. Notice that n = x2b2 + x3b3 = xaba is the position vector of an arbitrary point within a particular cross-section relative to the point in that cross-section where r interects it. A particle of the beam is then located from a fixed point in space by the position vector ^rðx1 ; x2 ; x3 Þ, given by

^r ¼ r þ xa ba ¼ r þ n

ð1Þ

The locus of material points along r has now assumed a different curved line denoted by R. Let s denote arc-length along R. The locus of points belonging to the initially planar refernce cross-section of the undeformed beam has undergone a rigid-body translation and rotation, as well as a warping displacement. The rigid body transla  tion is along the vector u x1 , the position vector from the point on the undeformed beam reference curve at x1 ¼ x1 to the point on the

  deformed beam reference curve at s ¼ s x1 . At each point along R introduce the frame B in which are fixed orthogonal unit vectors   Bi(x1) for i = 1, 2, 3, with B1 x1 normal   to the deformed beam reference cross-sectional plane and Ba x1 lying in this plane. Note that B1 = B2  B3 is not necessarily tangent to R unless one adopts the Euler–Bernoulli hypothesis: that the reference cross-section remains normal to R when the beam is deformed. In addition to the rigid-body motion of the set of particles making up the reference cross-section of the undeformed beam, this initially plane set of particles warps. The frame B is chosen so that the portion of the displacement relegated to the warping is small, so that the deformed beam reference cross-sectional plane is the plane that is closest to those material points of the deformed beam that make up the reference cross-section of the undeformed beam at x1 of the undeformed beam. In order to represent the deformed state mathematically, we must first deal with the rotation. Rotation from bi to Bi is accomplished by pre-dot multiplication with an orthogonal tensor that we call the global rotation tensor CBb:

Bi ¼ CBb :bi ¼ C Bb ij bj

ð2Þ

Introduce R = r + u, where u = uBiBi = uAiAi = uibi is the displacement of a point on the reference curve, one can represent the position of a particle in the deformed beam that had position ^r in the undeb 1 ; x2 ; x3 Þ, given by formed beam as Rðx

^ ¼ R þ CBb :ðn þ wÞ ¼ r þ u þ x2 B2 þ x3 B3 þ wi Bi R

ð3Þ

where w = wibi is a vector that represent the small warping displacement field such that wi = wi(x1, x2, x3). The crucial step in the cross-sectional analysis is the expression of the 3-D strain field of the beam in terms of 1-D generalized strains whose definitions agree with those from the Cosserat theory, a director theory. The vector-dyadic definitions of force and moment generalized strains can be expressed as:

c ¼ CbB :R0  r0 j ¼ CbB :K  k

ð4Þ

where ()’ denotes the derivative w.r.t. x1 (for vectors, the derivative is taken in (A) and K is the curvature vector of the deformed beam, defined to be analogous to angular velocity, so that

e ¼ ðCBA Þ0 :CAB K

ð5Þ

Note that the derivative is a partial derivative when time enters the problem. When K is expressed in the B basis, K1 is the twist per unit length and K2 and K3 are curvatures of the deformed beam. Similarly, k is the curvature vector for the undeformed beam

~ ¼ ðCbA Þ0 :CAb k

ð6Þ

and when k is expressed in the b basis, k1 is the twist per unit length and k2 and k3 are curvatures of the undeformed beam. The fact that b1 is tangent to r allows the simplification that

~ ¼e r0b þ kr 1 b

ð7Þ T

where e1 = [1 0 0] so that one of the simplest forms of c becomes

~  e1 c ¼ C Bb ðe1 þ u0 þ kuÞ

ð8Þ

Now

j¼K k

ð9Þ

where

~ bB e ¼ ðC Bb Þ0 C bB þ C Bb kC K

Fig. 1. Schematic of the beam deformation.

ð10Þ

In both the force and moment strains, the initial curvature vector measures in the b basis, contained in the column matrix k, are known.

4

H. Pollayi et al. / Composite Structures 100 (2013) 1–16

3. Initially curved and twisted strip-like beam cross-sectional analysis (2-D, non-linear) The closed-form analytical analysis is restricted to thin rectangular cross-sections for generally anisotropic composites. Fig. 2 shows the sketch of an initially twisted and curved strip configuration, coordinate system and cross-section. The purpose of the non-linear analysis of the beam cross-section is to determine the elements of the non-linear 1-D stiffness matrix, SNL, and the recovering relations in closed-form. Recovering relations, for example, provide the relationship between the 3-D strain tensor components, Cij(i, j = 1, 2, 3), and the generalized 1-D strain measures, ew. In the present notation, the 1-D constitutive law can be written as

f H ¼ SNL eH

h i 1 T 1 k ‘ ½S‘ ‘ þ T‘ ½S‘n n þ Tn ½Sn n þ 2 T‘ Sc1 ‘ ‘ 2 2 2 2 h i  c k T þ 2 T‘ Sc2 ‘ ‘ þ k2 ‘ S‘n n 2

ð12Þ

where k2 is initial curvature in flatwise direction of strip, and the linear and non-linear 1-D strain vectors, ‘ and n, respectively, are defined as follows:

‘ ¼ ½c11 j 1 j 2 j 3 T n ¼ ½j 21 j 22 j 2 c11 j 2 j 3 j 2 j 1 T

ð14Þ

The cross-sectional stiffness matrix is obtained by expressing the resultant forces on the beam cross-section as

f H D½F 1 F 2 F 3 M1 M2 M 3 T ¼

 T @U 1D 1 @U 1D 1 @U 1D @U 1D @U 1D @U 1D @ c11 2 @ c12 2 @ c13 @ j1 @ j2 @ j3 ð15Þ

This is the most general form of the cross-sectional analysis for class S and class T beams defined by Hodges [11]. The coefficients of the non-linear stiffness matrix, SNL, obtained from the non-linear analysis of the beam cross-section is defined in the Appendix. All the 1D stiffness matrices and their corrections are defined in the work of Harursampath [13].

ð11Þ

where fw = [F1 F2 F3 M1 M2 M3]T with F1 as the axial force, F2 and F3 as transverse shear forces, M1 as the twisting moment, M2 and M3 as the bending moments; ew = [c11 2c12 2c13 j1 j2 j3]T with c11 as the axial stretching measure of the beam, 2c12 and 2c13 as transverse shear measures, j1 as the elastic twist per unit length, j2 and j3 as elastic bending components of the curvature. The final expression for the strain energy per unit length (or 1-D strain energy density) of the curved and twisted beam component is

U 1D ¼

ðÞ ¼ ðÞj2c12 ¼2c13 ¼0

ð13Þ

11 ; j  1; j  2 and j  3 relate to their unHere, the barred quantities c barred counterparts as

4. Initially curved and twisted strip-like beam analysis (1-D, non-linear) The dynamic response of non-linear, flexible multi-body systems is simulated within the framework of energy-preserving and energy-decaying time integration schemes that provide unconditional stability for non-linear systems. The kinematic description of bodies and joints in their undeformed and deformed configurations adopted here follows Bauchau [15] and makes use of three orthogonal triads. Fig. 3 shows the beam in the undeformed and deformed configurations. First, an inertial triad, S I , is used as a global reference for the system with unit vectors, ~ ı1 , ~ ı2 , and ~ ı3 . A second triad, S 0 , is attached to the body and defines its orientation in the reference configuration with unit vectors, ~ e01 ; ~ e02 , and ~ e03 . Finally, a third triad, S  , defines the orientation of the body in its deformed configuration with unit vectors, ~ e1 , ~ e2 , ~ measured in S I and S  are deand ~ e3 . The components of vector ðÞ noted by () and ()⁄, respectively. For completeness, we outline the relevant portions of Bauchau [15]. The kinetic and strain energies of the beam are

K¼

1 2

Z

‘

tT M t dx1 ; U ¼

0

Z

‘

11 ; j  1; j  2; j  3 Þdx1 U 1D ðc

ð16Þ

0

respectively, where ‘is the length of the beam; x1 is the curvilinear coordinate along the reference curve; M⁄ and tw are the sectional inertia matrix and velocity vector, respectively. The equations of motion of the beam can be derived from the Hamilton’s principle as follows:

Fig. 2. Sketch of (a) initially twisted strip configuration and coordinate system and (b) its cross-section.

Fig. 3. Beam in the undeformed and deformed configurations.

H. Pollayi et al. / Composite Structures 100 (2013) 1–16

Z ti

tf

Z

5

‘

ðdtT M t  dU 1D þ dW a Þdx1 dt ¼ 0

ð17Þ

0

where dt⁄T are the virtual variations in sectional velocities, dU 1D is the variation of 1-D strain energy density, and dW a is the virtual work done by the externally applied forces. The final form of the equations of motion of the beam is as follows:   0 ~_ ~00 þ u ~ 0 ÞRR0 f  ¼ q ðRR0 p Þ þ U½uRR 0 p  ðRR0 f Þ  U½u

ð18Þ

where u0 the displacement vector from a global reference for the system to body attached reference configuration and u the displacement vector from body attached reference configuration to deformed configuration and the sectional momenta and elastic forces are defined as p⁄ = M⁄t⁄ and as before f⁄ = SNLe⁄ respectively; and q are the external forces. An energy-preserving discretization of these equations of motion, Eq. (18), is performed. The inertial and elastic forces are discretized to yield the following discretized equations of motion

Rf R0 pf  Ri R0 pi

Dt

þU

  0 pf þ pi  ~i ~f  u u Ra R0  Rb R0 f m Dt 2

~ 00 þ u ~ 0m ÞRb R0 f m  U½u ¼ qm

ð19Þ

where ()m = [()f + ()i]/2; and all the rotation operators and the discretization of finite rotations are defined in Appendix A and B of Bauchau [15]. For the case of revolute joint elements, the two kinematic constraints, C1 ¼ 0 and C2 ¼ 0, are defined in the deformed configuration utilizing the condition of no relative displacements and the third constraint, C3 ¼ 0, obtained from the definition of the relative rotation / between the two bodies where Ci ði ¼ 1; 2; 3Þ are defined as follows: l C1 ¼ g 31 ¼ ekT 3 e1 ¼ 0 l C2 ¼ g 32 ¼ ekT 3 e2 ¼ 0

C3 ¼ g 11 sin / þ g 12 cos / ¼ 0

ð20Þ

l kT l where g 11 ¼ ekT 1 e1 and g 12 ¼ e1 e2 . Here (e1,e2,e3) form orthonormal bases in the deformed configurations for two bodies k and l. The above formulation is explained in detail in Pollayi and Harursampath [18].

5. Recovery analysis (3-D, approximate asymptotically) To recover any or all 3-D results at a specific section, one needs to provide information obtained both from the 1-D global multibeam system analysis along with warping and strain recovery relations, a by-product of 2-D non-linear cross-sectional analysis. Using closed-form analytical expressions for 3-D warping and stress field derived by Harursampath [13], one can recover the in-plane and out-of-plane warping of the cross-section, the total displacement field, the 3-D strain and stress fields. The above mentioned complete analysis of non-linear, flexible multi-body beam systems is depicted as a flow-chart in Fig. 4. 6. Failure criteria and strength of component-laminates Failure criteria associated with a lamina can be divided into two categories: [1] Failure criteria NOT directly associated with failure modes such as Tsai–Hill, Tsai–Wu, Modified Tsai–Wu, and all polynomial, tensorial, or parametric criteria. [2] Failure criteria associated with failure modes such as Maximum strain or stress, Hashin and Rotem, Yamada and Sun, Hashin, Hart-Smith, Puck, and Kriging.

Fig. 4. Flow-chart of variational asymptotic analysis of thin-composite beam systems.

6.1. Polynomial criterion of failure tensor A simple form-invariant expression is extensively used for the failure tensor polynomial function. The theoretical predictions may or may not seriously effect because of the rank of the tensors and the number of terms appearing in the respective equation. The following is the failure function by considering tensor polynomial coefficients up to the fourth rank

f ðrÞ ¼ r:H:r þ h:r  1 ¼ 0

ð21Þ

With an appropriate definition of the Cartesian components of the fourth-and second-rank tensors H and h, respectively, many anisotropic failure criteria can be written in the form of above relation. The stress differential effect second-rank tensor h, as well as the diagonal components of failure fourth-rank tensor H. H and h are uniquely determined for the criteria assuming the general form of the above equation. Differences arising only in the definition of off-diagonal components Hij(i – j) express the various phenomenological conjectures of these criteria.

Hij ¼ 0ði – jÞ for Tsai  Wu ðTWÞ Criterion 1 1 Hij ¼  ðHii Hjj Þ2 for Tsai  Hahn ðTHÞ Criterion 2 1 Hij ¼ ðHkk  Hii  Hjj Þði; j; k 6 3; i – j – kÞ for Elliptic Paraboloid 2 Failure Surface ðEPFSÞ Criterion

6

H. Pollayi et al. / Composite Structures 100 (2013) 1–16

The quadratic stress tensor polynomial failure criterion for anisotropic materials can be represented, using the usual contracted matrix notation, as

F i ri þ F ij ri rj ¼ 1ði; j ¼ 1; 2; . . . ; 6Þ where Fi and Fij are strength tensors of the second and fourth rank, respectively. Fi, Fij, and the recovered 3-D stress tensor, ri, are given below.

F Ti ¼ ½F 1 F 2 F 3 F 4 F 5 F 6  2

F 11

6 6 6 6 F ij ¼ 6 6 6 6 4

F 13

F 14

F 15

F 16

F 22

F 23

F 24

F 25

F 33

F 34

F 35

F 44

F 45

F 26 7 7 7 F 36 7 7 F 46 7 7 7 F 56 5

F 55

ð23Þ

ð24Þ

6.2. Expanding Normal Components of Failure Tensor

F1 ¼ F2 ¼ F3 ¼

rT1 rC1 1

rT1 1

rT2 1

rT3

  

1

rC1 1

rC2 1

rC3

ð29Þ

6.5. Non-dimensionalization of Tsai–Wu–Hahn criterion

pffiffiffiffiffiffiffi

rTi ¼ ½r1 r2 r3 r4 r5 r6  ¼ ½r11 r22 r33 r23 r31 r12 

; F 22 ¼

1 1 F 36 ¼  ðF 33 F 66 Þ2 2

The equation for non-dimensionalization of Tsai–Wu–Hahn failure criterion is derived by using the following definitions:

3

F 12

sym

1

1 1 F 35 ¼  ðF 33 F 55 Þ2 ; 2 1 1 F 46 ¼  ðF 44 F 66 Þ2 2

ð22Þ

F 66

F 11 ¼

1 1 F 34 ¼  ðF 33 F 44 Þ2 ; 2 1 1 F 45 ¼  ðF 44 F 55 Þ2 ; 2 1 1 F 56 ¼  ðF 55 F 66 Þ2 2

1

rT2 rC2

; F 33 ¼

1

ð25Þ

rT3 rC3

¼ ðrC1  rT1 ÞF 11

F1 rf1 D F 11 r1 ; f F 1 D pffiffiffiffiffiffiffi

F 11 F2 f F 2 D pffiffiffiffiffiffiffi F 22 pffiffiffiffiffiffiffi F 3 rf3 D F 33 r3 ; f F 3 D pffiffiffiffiffiffiffi F 33 pffiffiffiffiffiffiffi F4 g r F 4 D pffiffiffiffiffiffiffi F 44 r23 ; f 23 D F 44 pffiffiffiffiffiffiffi F 5 g r F 5 D pffiffiffiffiffiffiffi F 55 r13 ; f 13 D F 55 pffiffiffiffiffiffiffi F6 g F 66 r12 ; f r F 6 D pffiffiffiffiffiffiffi 12 D F 66 pffiffiffiffiffiffiffi rf2 D F 22 r2 ;

ð30Þ

Non-dimensionalize and substitute in non-dimensional form of Tsai–Wu–Hahn failure criterion implies a quadratic equation, quadratic in reserve factor, R, and is as follows:

AR2 þ BR  1 ¼ 0

¼ ðrC2  rT2 ÞF 22

ð26Þ

ð31Þ

where 2 2 2 f1 2 þ r f2 2 þ r f3 2 þ r g g g f1 r f1 r f1 r f2 þ r f3 þ r g A¼r 23 þ r 13 þ r 12  ½ r 23

¼ ðrC3  rT3 ÞF 33

where rTi is Tension (T) failure stress in i-direction (i = 1, 2, 3) and rCi is Compression (C) failure stress in i-direction (i = 1, 2, 3).

f1 r f1 r f2 r f2 r f2 r f2 r f3 r g g f3 þ r g g g g þr 13 þ r 12 þ r 23 þ r 13 þ r 12 þ r 23 f f g g g g g g g g þ r3 r13 þ r3 r12 þ r23 r13 þ r23 r12 þ r13 r12  fg fg f1 þ f f2 þ f f3 þ f g B¼f F1 r F2 r F3 r F4 r 23 þ F 5 r 13 þ F 6 r 12

6.3. Expanding shear components of failure tensor 6.6. The procedure

F 44 ¼ F4 ¼ F5 ¼ F6 ¼

1

;

rSþ4 rS4 1

rSþ4 1

rSþ5 1

rSþ6

  

1

r

S 4

1

rS5 1

rS6

F 55 ¼ ¼ ¼ ¼



1

rSþ5 rS5

;



F 66 ¼

1

ð27Þ

rSþ6 rS6

rS4  rSþ4 F 44









rS5  rSþ5 F 55

ð28Þ

rS6  rSþ6 F 66

where rSþ is Shear (S) strengths + ve in i-plane (i = 4, 5, 6) and rSi is i Shear (S) strengths-ve in i-plane (i = 4, 5, 6).

1 1 F 13 ¼  ðF 11 F 33 Þ2 ; 2 1 1 F 16 ¼  ðF 11 F 66 Þ2 2 1 1 F 24 ¼  ðF 22 F 44 Þ2 ; 2

1. [F M]T = The stress resultants obtained from 1-D global beam(s) analysis as distributions w.r.t. time and all three beam reference curves0 coordinates. 2. Find F rmax and M rmax (Available @ different times). 3. Input [F M]T twice (for F rmax and M rmax ) to cross-sectional analysis code based on current work by keeping Recovery Flag = True to obtain recovered 3-D strain and hence stress results. 4. Determine stresses in the principal material directions using the above recovered 3-D stresses by applying the following stress transformation relations:

8 9 ðnÞ > > > < rL > =

6.4. Expanding F ij ði – jÞ

1 1 F 12 ¼  ðF 11 F 22 Þ2 ; 2 1 1 F 15 ¼  ðF 11 F 55 Þ2 ; 2 1 1 F 23 ¼  ðF 22 F 33 Þ2 ; 2 1 1 F 26 ¼  ðF 22 F 66 Þ2 2

The following is the step-wise procedure used to predict the first-ply-failure and the failure of present mechanism as a whole.

1 1 F 14 ¼  ðF 11 F 44 Þ2 2

F 25

1 1 ¼  ðF 22 F 55 Þ2 2

> > :

rðnÞ T > > ; rðnÞ LT

2

2

cos2 aðnÞ

6 ¼4

sin aðnÞ

2

cos2 aðnÞ

sin aðnÞ sina cosa ðnÞ

ðnÞ

sinaðnÞ cosaðnÞ

9 38 > < r11 > = 7 sin2aðnÞ 5 r22 > > : r12 ; cos2aðnÞ sin2aðnÞ

ð32Þ 8 ðnÞ 9 > r > > < LT 0 > = > > :

2

cosaðnÞ

sinaðnÞ

6 ¼ 4 sinaðnÞ cosaðnÞ rðnÞ TT 0 > > ; 0 0 rðnÞ T0

9 38 0 < > r13 > = 7 0 5 r23 > > : r33 ; 1

ð33Þ

7

H. Pollayi et al. / Composite Structures 100 (2013) 1–16

Here a(n) denotes the ply-angle of nth ply. 5. Determine strains in the principal material directions are determined using the above recovered 3-D strains by applying the following strain transformation relations:

8 9 ðnÞ > > > < L > = > > :

2

2

cos2 aðnÞ

6 2 ¼4 ðnÞ sin aðnÞ T > > ðnÞ ðnÞ ; sina cosaðnÞ LT

sin aðnÞ cos2 aðnÞ sinaðnÞ cosaðnÞ

9 38 > < 11 > = 7 sin2aðnÞ 5 22 > > : 12 ; cos2aðnÞ sin2aðnÞ

ð34Þ 8 ðnÞ 9 >  > > < LT 0 > = > > :

2

cosaðnÞ

6 ¼ 4 sinaðnÞ ðnÞ TT > > ðnÞ ; 0 T 0

0

sinaðnÞ cosaðnÞ 0

9 38 0 < > 13 > = 7 0 5 23 > > : 33 ; 1

Fig. 6. The four-bar mechanism problem.

ð35Þ

Here a(n) denotes the ply-angle of nth ply. 6. Calculate reserve factor (R) in each ply by applying Tsai–Wu–Hahn 3-D failure criterion using above calculated stresses in the principal material directionsi:e: h iT ðnÞ ðnÞ ðnÞ ðnÞ ðnÞ ðnÞ . rL rT rLT rLT 0 rTT 0 rT 0 7. Then RLaminate = min. of {R1, R2, R3, . . . , Rn} for each of the three bars in the mechanism. Fig. 5 shows the complete algorithm for finding the failure stress level in the component-laminate of the present problem, the four-bar mechanism. 7. Numerical results and discussions The numerical example adopted from Bauchau [15] deals with a four-bar mechanism problem depicted in Fig. 6. Bar 1, (AB) is of length ‘1 = 0.12 m and is connected to the ground at point A by means of a revolute joint. Bar 2, (BC) is of length ‘2 = 0.24 m and is connected to bar 1 at point B with a revolute joint. Finally, bar 3, (CD) is of length ‘3 = 0.12 m and is connected to bar 2 and the ground at points C and D, respectively, by means of two revolute joints. At time, t = 0, in the undeformed configuration, the bars of this mechanism intersect each other at 90° angles. However, the axes of rotation of all the revolute joints are at non-zero angles with respect to the normal to the plane of the mechanism to simulate initial defects in the mechanism. Such initial defects shall be defined in the form of initial twists (k1) and curvatures (k2) of all

the three bars. A torque is applied on bar 1 at point A so as to enforce a constant angular velocity X = p rad/s. For elastic bars, complete motion becomes possible, but generates large internal forces. The laminated composite strip of bar-1 is similar to that considered by Armanios et al. [19] and was fabricated from the ICI Fiberite T300/954-3 Graphite/Cyanate material system, whereas the laminated composite strips of bar-2 and bar-3 are similar to those considered by Kosmatka [20] and were made from the T300/5208 Graphite/Epoxy material system. In this work, the stacking sequence [a2/(a  90°)4/a2/  a2/(90°  a)4/  a2]T is used for bar-1 and bar-3 whereas for bar-2 it is [a/(a  90°)2/a/  a/ (90°  a)2/  a]T where a = 20°. Each strip has a uniformly distributed pre-twist of 5° about the longitudinal axis and the initial flatwise curvature causing an angle between its ends is 2°. The physical characteristics of the three bars are tabulated in Table 1. 0 0 0 The indices 1, 2, and 3 suffixing E s, G s, and m s denote the principal material axes. The specific type of antisymmetrical layup used for all the bars ensures hygrothermal stability of the laminates (see Winckler [21]) and eliminates the initial warping that results from the curing stresses while exhibiting extension-twist coupling under an applied mechanical load. Thus all the three bars are antisymmetric laminates and have the following properties:

A16 ¼ A26 ¼ 0;

B11 ¼ B22 ¼ B66 ¼ B12 ¼ 0;

D16 ¼ D26 ¼ 0

ð36Þ

For comparing linear stiffness coefficients from FE code VABS with the results from the analytical expressions for non-linear stiffness coefficients, the cross-sections of bar-1 and bar-3 are divided into 160 eight-node isoparametric finite elements with 10 elements each along the width direction and 16 elements along the thickness direction. Similarly, bar-2 is divided into 80 eight-node isoparamet-

Fig. 5. Algorithm for finding the failure stress level in component-laminate of four-bar mechanism.

8

H. Pollayi et al. / Composite Structures 100 (2013) 1–16

Table 1 Material properties and dimensions of all three bars in the mechanism. Property (dimension)

Bar-1

Bar-2

Bar-3

EL (GPa) ET ¼ ET 0 (GPa) GLT ¼ GLT 0 (GPa) GTT 0 (GPa) mLT ¼ mLT 0

135.6 9.9 4.2 3.3 0.3 0.5 0.12 0.012 0.0012 0.075 1520.6 0.7275 0.291

132.2 10.75 5.65 5.65 0.239 0.4 0.24 0.012 0.0006 0.075 1530.8 0.3638 0.1455

132.2 10.75 5.65 5.65 0.239 0.4 0.12 0.012 0.0012 0.075 1530.8 0.7275 0.291

mTT 0 ‘ (m) b (m) h (m) Ply thickness (mm) q (kg/m3) k1 (rad/m) k2 (rad/m)

Table 3 Stiffness coefficients of all three bars in the mechanism: 1-Extension, 4-Torsion, 5Flatwise Bending, 6-Edgewise Bending. The bold-faced values in brackets are from the linear cross-sectional analysis using FE code VABS with refined Timoshenko-like beam model. Stiffness

Present work

6 SNL 11 ð10 NÞ

SNL 14 ðN:mÞ

SNL 15 ðN:mÞ SNL 16 ðN:mÞ

2 SNL 45 ðN:m Þ

2 SNL 46 ðN:m Þ

2 SNL 55 ðN:m Þ

2 SNL 56 ðN:m Þ

2 SNL 66 ðN:m Þ

0.4053

Bar-3 0.8107

(0.7723)

(0.3996)

(0.7848)

158.804

36.5353

146.141

+9.7351j1 (151.5718)

+4.864j1 (38.1372)

+9.7281j1 (143.8132)

0

0

0

(0.05956)

(0.0003808)

(0.04814)

0

0

0

(0.0)

(0.0)

(0.0)

0.09468

0.0127

0.1016

+9.7351c11 +0.005513j1 +0.000008575j2 +0.0003154j21

+4.864c11 +0.00115j1 +0.000001853j2 +0.0001576j21

+9.7281c11 +0.005057j1 +0.000007412j2 +0.0003152j21

+0.00002949j22 (0.0877)

+0.00001271j22 (0.0127)

+0.00002549j22 (0.09564)

0.000006248

0.0000006752

0.000005402

+0.000008575j1 0.00004294j2 0.0002658j3 +0.00005898j1j2 (0.0002549)

+0.000001853j1 0.000009282j2 0.0002574j3 +0.00002549j1j2 (0.00003071)

+0.000007412j1 0.00003713j2 0.0002574j3 +0.00005099j1j2 (0.0002318)

0.0002658j2

0.0002574j2

0.0002574j2

(0.0)

(0.0)

(0.0)

0.07402

0.0102

0.08161

0.00004294j1 +0.00001627j2 +0.00002949j21

0.000009282j1 +0.000003043j2 +0.00001275j21

 0.00003713j1 +0.00001217j2 +0.00002549j21

+0.00005592j22

+0.00002091j22

+0.00004183j22

0.0009134j23 (0.1069)

0.001769j23 (0.01371)

0.0008844j23 (0.1082)

0.0002658j1

0.0002574j1

0.0002574j1

0.001827j2j3 (0.0)

0.003538j2j3 (0.0)

0.001769j2j3 (0.0)

9.735

4.864

9.728

0.0009134j22 (8.5317)

0.001769j22 (4.6761)

0.0008844j22 (8.9096)

3

2

0

0

0.5

1

1.5

2

2.5

Time [s]

Table 2 The sectional mass properties of all three bars in the mechanism. Property

Bar-1

Bar-2

Bar-3

l (kg/m)

0.218966E01 0.262760E08 0.262760E06

0.110218E01 0.330653E09 0.132261E06

0.220435E01 0.264522E08 0.264522E06

i2 (kg m2/m) i3 (kg m2/m)

Bar-2

0.8113

4

Rotations [rad]

ric finite elements with 10 elements along the width direction and eight elements along the thickness direction. The sectional mass properties are given in Table 2 for all three bars. Here l is the mass per unit length, i2 is the mass moment of inertia about the direction along the width, and i3 is the mass moment of inertia about the direction along the thickness. The stiffness matrices obtained from the linear cross-sectional analysis of bars 1, 2, and 3 are given in Table 3 in bold-face and kept inside brackets. The non-linear crosssectional stiffness matrix has strong extension-twist coupling, as expected. The beam-1, beam-2, and beam-3 reference curves are divided into 4-equal elements each for 1-D non-linear analysis. For this four-bar mechanism problem, the automated time step size procedure is used. The desired local error level for the new time step size is set to be ^e ¼ 1:0  1006 and a new time step size is evaluated at each time step. This problem is simulated for a total of 8.0 s using the energydecaying scheme. If the four revolute joints had their axes of rotation orthogonal to the plane of the mechanism, the response of the system would be purely planar, and bars 1 and 3 would rotate at constant angular velocities around points A and D, respectively. The initial twist and/or curvature of bars causes a markedly different response. Fig. 7 shows the time history of the relative rotations at points A and D, as well as the absolute rotation of the mid-point of bar 2. Bar 1 rotates at a constant angular velocity under the effect of the applied torque, but bar 3 now oscillates back and forth, never completing an entire turn. When the direction of rotation of bar 3 reverses, bar 2 undergoes large rotations, instead of mere translation. In both linear and present non-linear analyses, the behavior of the mechanism is very similar and one can observe smooth time histories of rotations of the system. The authors have simulated the four-bar mechanism(s) using the available commercial software (I-DEAS + NASTRAN + ADAMS) to see the practical applicability by arranging the beams in the mechanism(s) in distinctly different orientations. Fig. 8 shows the noticeable snapshots of the deformation shapes of the four-bar mechanism problem about 1000 frames in 1 sec for three complete rotations of beam-1 about point ‘A’. It was observed that beam-1 makes a full rotation about point ‘A’ and beam-3 oscillates back and forth, never completing an entire turn. Only beam-2 undergoes large deformations and large rotations. Noticeable deformations of beam-2 start

2 SNL 44 ðN:m Þ

Bar-1

Linear : Relative Rotation@A Linear : Relative Rotation@D Linear : Absolute Rotation of Beam−2 Non−Linear : Relative Rotation@A Non−Linear : Relative Rotation@D Non−Linear : Absolute Rotation of Beam−2

Fig. 7. Time history of rotations of the system.

3

3.5

H. Pollayi et al. / Composite Structures 100 (2013) 1–16

at about t = 0.098 s for the first rotation. The maximum deformation of beam-2 occurred about t = 0.243 s in the first rotation. Reserve factor can be defined as a common factor by which all of the stresses applied on a lamina need to be multiplied by for the lamina to just fail. The ply in the component-laminate resists the combination of stresses without failure if R > 1. The reserve factor of each ply in the component-laminate is calculated by using the

9

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi maximum resultant force, F r ¼ F 21 þ F 22 þ F 23 and the maximum qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi resultant moment, Mr ¼ M 21 þ M 22 þ M 23 instead of choosing either maximum of {F1, F2, F3} or maximum of {M1, M2, M3}. The positive root for R shows the failure stress level with loads in the same direction as the initial values and the negative root shows the failure stress level with all loads reversed. The component-lam-

Fig. 8. Snap-shots of the deformation shapes of the four-bar mechanism at various time instants.

10

H. Pollayi et al. / Composite Structures 100 (2013) 1–16

Table 4 Laminate strength properties. Laminate

Strength parameters (MPa)

T300/954-3 T300/5208

rT1

rC1

rT2

rC2

rSþ6

1314 1500

1220 1500

43 40

168 246

48 68

inate strength properties are given in Table 4 and these values are used to evaluate the component-laminate load-carrying capacity in the present four-bar mechanism problem. The Reserve Factor, R, of

component-laminates in the present four-bar mechanism are given in Table 5 as point-wise. In this work, the component-laminate load-carrying capacity is evaluated based on maximum Fr and maximum Mr using the Tsai– Wu–Hahn failure criterion for various lay-ups of the four-bar mechanism. The non-linear through-the-thickness stress and strain distributions are reported at different points. Tsai–Wu–Hahn failure criterion is used to predict first ply failure and mechanism as a whole. Fig. 9(i) shows the strain distribution through-the-laminate-thickness at ‘A’ of bar-1 (AB) for maximum resultant force, Fr = 38.03 N which occurs at t = 1.51 s. Fig. 9(ii) shows the stress distribution through-the-laminate-thickness at ‘A’ of bar-1 (AB)

Table 5 Reserve factor (R) of component-laminates of the present four-bar mechanism. Point

Frmax

A (AB) B (AB) B (BC) C (BC) C (CD) D (CD)

38.03 N (1.51 s) 38.08 N (1.51 s) 38.08 N (1.51 s) 38.29 N (1.51 s) 38.16 N (1.51 s) 38.21 N (1.51 s)

R (Clock)

R (Anti-Clock)

85.24

44.28

771.22

1474.13

400.24

215.62

63.44

36.57

172.81

211.8

124.65

146.43

Mrmax

R (Anti-Clock)

14.86

8.89

224.42

325.04

35.06

14.46

28.89

51.15

150.39

188.9

63.29

51.02

1.2e+06

2.5e−05

1e+06

2e−05

800000

1.5e−05

600000

1e−05

400000

5e−06

200000

0

0

−5e−06

−200000

−1e−05

−400000

−1.5e−05

−600000

−2e−05

−800000

−2.5e−05 −0.0006

−0.0004

−0.0002

0

0.0002

0.0004

−1e+06 −0.0006

0.0006

−0.0004

−0.0002

0

(i)

0.0002

0.0004

(ii)

Fig. 9. Stress–strain distribution through-the-laminate-thickness (@A of AB) for max. Resultant force, Fr = 38.03 N @ t = 1.51 s.

1200

0 Clockwise

1000

−100

Reserve Factor

Reserve Factor

Strain distribution

R (Clock)

0.47 N m (0.005 s) 0.03 N m (0.0005 s) 0.1005 N m (0.00056 s) 0.083 N m (0.0038 s) 0.041 N m (1.49 s) 0.122 N m (0.0002 s)

800 600 400

−200 −300 −400

200

−500

0

−600

Anti−Clockwise

1

2

3

4

5

6

7

8

9 10 11 12 13 14 15 16

1

2

3

4

5

6

7

8

9 10 11 12 13 14 15 16

Ply No.

Ply No.

(i)

(ii)

Fig. 10. Reserve factor (R) in each ply (RLaminate = +85.24 (clockwise), 44.28 (anti-clockwise)).

0.0006

11

Strain distribution

H. Pollayi et al. / Composite Structures 100 (2013) 1–16

0.00015

8e+06

1e−04

6e+06

5e−05

4e+06

0

2e+06

−5e−05

0

−0.0001

−2e+06

−0.00015 −0.0006

−0.0004

−0.0002

0

0.0002

0.0004

−4e+06 −0.0006

0.0006

−0.0004

−0.0002

0

0.0002

0.0004

0.0006

(ii)

(i)

Fig. 11. Stress–strain distribution through-the-laminate-thickness (@A of AB) for max. Resultant moment, Mr = 0.47 N m @ t = 0.005 s.

160

0 Clockwise

140

−20

Reserve Factor

Reserve Factor

120 100 80 60

−40 −60 −80

40 −100

20

Anti−Clockwise

0

1

2

3

4

5

6

7

8

9 10 11 12 13 14 15 16

−120

1

2

3

4

5

6

7

8

9 10 11 12 13 14 15 16

Ply No.

Ply No.

(i)

(ii)

Fig. 12. Reserve factor (R) in each ply (RLaminate = +14.86 (clockwise), 8.89 (anti-clockwise)).

200000

3e−06

150000

Strain distribution

2e−06

100000

1e−06

50000 0 0 −1e−06

−50000

−2e−06

−100000

−3e−06

−150000

−4e−06 −0.0003

−0.0002

−0.0001

0

0.0001

0.0002

0.0003

(i)

−200000 −0.0003

−0.0002

−0.0001

0

0.0001

0.0002

0.0003

(ii)

Fig. 13. Stress–strain distribution through-the-laminate-thickness (@B of BC) for max. Resultant force, Fr = 38.08 N @ t = 1.51 s.

for maximum resultant force, Fr=38.03 N which occurs at t = 1.51 s. The distributions are significantly affected by cross-sectional nonlinearity. Here it can be observed that the layer orientation influences the stress distribution of the component-laminate and it is also very important to determine the location of highly stressed regions in the component-laminate and the stress values in all the layers of the system. Stresses rL and rLT reach their maximum values at the bottom of the laminate of bar-1 (AB) whereas stress rT is maximum at 0.4 mm from the bottom of the same laminate. The

highest stress, as expected, is rL due to bending in this laminate. Fig. 10(i) shows the reserve factor, R, in each ply at ‘A’ of the laminate ‘AB’ for the applied clockwise constant angular velocity, X = p rad/s. In this case, ply-1 of bar-1 is predicted to fail first for maximum Fr and its R = +85.24. For maximum Fr and with reversal of the applied constant angular velocity, ply-3 of bar-1 is predicted to fail first and its R = 44.28 as observed from Fig. 10(ii). Fig. 11(i) shows the strain distribution through-the-laminatethickness at ‘A’ of bar-1 (AB) for the case of maximum resultant

12

H. Pollayi et al. / Composite Structures 100 (2013) 1–16

3000

0 Clockwise

−500

Reserve Factor

Reserve Factor

2500 2000 1500 1000

−1000 −1500 −2000

500

−2500

0

−3000

Anti−Clockwise

1

2

3

4

5

6

7

8

1

2

3

4

5

Ply No.

Ply No.

(i)

(ii)

6

7

8

Fig. 14. Reserve factor (R) in each ply (RLaminate = +400.24 (clockwise), 215.62 (anti-clockwise)).

1e−04

4e+06

8e−05

3e+06

Strain distribution

6e−05 2e+06

4e−05 2e−05

1e+06

0 0

−2e−05 −4e−05

−1e+06

−6e−05

−2e+06

−8e−05 −0.0001 −0.0003

−0.0002

−0.0001

0.0001

0

0.0002

−3e+06

−0.0003

0.0003

−0.0002

−0.0001

0

(i)

0.0001

0.0002

0.0003

(ii)

Fig. 15. Stress–strain distribution through-the-laminate-thickness (@B of BC) for max. Resultant moment, Mr = 0.1005 N m @ t = 0.00056 s.

200

0

180

Clockwise

−20

140

Reserve Factor

Reserve Factor

160 120 100 80 60 40

−40 −60 −80 −100

20 0

Anti−Clockwise

1

2

3

4

5

6

7

8

−120

1

2

3

4

5

Ply No.

Ply No.

(i)

(ii)

6

7

8

Fig. 16. Reserve factor (R) in each ply (RLaminate = + 35.06 (clockwise), 14.46 (anti-clockwise)).

moment, Mr = 0.47 N m occuring at t = 0.005 s. Fig. 11(ii) shows the corresponding stress distribution through-the-laminate-thickness at ‘A’ of bar-1 (AB). In this case also, the stress and strain distributions are significantly affected by cross-sectional non-linearities. Stresses rL and rLT reach their maximum values at the bottom of the laminate of bar-1 (AB) whereas stress rT is maximum at 0.4 mm from the bottom of the laminate. The highest stress, as expected, is rL due to bending in this laminate. Fig. 12(i) shows the reserve factor, R, in each ply at ‘A’ of the laminate ‘AB’ for the applied clockwise constant angular velocity, X = p rad/s. In this case

also, ply-1 of bar-1 is predicted to fail first for maximum Mr and its R = + 14.86. For maximum Mr and with the reversal of the applied constant angular velocity, ply-3 of bar-1 is predicted to fail first and its R = 8.89 as observed from Fig. 12(ii). For bar-2 (BC) at ‘B’, the strain and stress distributions throughthe-laminate-thickness for maximum resultant force, Fr = 38.08 N occurs at t = 1.51 s are shown in Fig. 13(i) and (ii), respectively. The strain component, L is constant and its value is equal to 1.0E06 through-the-thickness of the laminate. The T value is increasing from 2.0E06 and converging L at the top of the

13

Strain distribution

H. Pollayi et al. / Composite Structures 100 (2013) 1–16

4e−06

400000

3e−06

300000

2e−06

200000

1e−06

100000

0 0 −1e−06 −100000

−2e−06

−200000

−3e−06 −4e−06 −0.0006

−0.0004

−0.0002

0

0.0002

0.0004

−300000 −0.0006

0.0006

−0.0004

−0.0002

0

(i)

0.0002

0.0004

0.0006

(ii)

Fig. 17. Stress–strain distribution through-the-laminate-thickness (@D of CD) for max. Resultant force, Fr = 38.21 N @ t = 1.51 s.

14000

0 Clockwise

−1000

12000

Reserve Factor

Reserve Factor

−2000 10000 8000 6000 4000

−4000 −5000 −6000 −7000 −8000

2000 0

−3000

−9000 1

2

3

4

5

6

7

8

9 10 11 12 13 14 15 16

Anti−Clockwise

−10000

1

2

3

4

5

6

7

8

9 10 11 12 13 14 15 16

Ply No.

Ply No.

(i)

(ii)

Fig. 18. Reserve factor (R) in each ply (RLaminate = + 124.65 (clockwise), 146.43 (anti-clockwise)).

Strain distribution

laminate. At the same ‘B’, the reserve factor, R, in each ply of the laminate ‘BC’ for the applied clockwise and with the reversal of the applied constant angular velocity are shown in Fig. 14(i) and (ii), respectively. At the same point, ‘B’, for maximum resultant moment, Mr = 0.1005 N m occurs at t = 0.00056 s, the strain and stress distributions are shown in Fig. 15(i) and (ii), respectively. The Reserve Factor, R, in each ply at ‘B’ of the laminate ‘BC’ for the applied clockwise and with the reversal of the applied constant angular velocity are shown in Fig. 16(i) and (ii), respectively. In both these cases, ply-7 of bar-2 is predicted to fail first for maximum Fr and Mr and its Reserve Factors R = +400.24 and R = +35.06, respectively.

For maximum Fr and Mr with the reversal of the applied constant angular velocity, ply-2 of bar-2 is predicted to fail first in both cases and its R = 215.62 and R = 14.46, respectively. At ‘D’ of bar-3 (CD), for maximum resultant force, Fr = 38.21 N occurs at t = 1.51 s and for maximum resultant moment, Mr = 0.122 N m occurs at t = 0.0002 s. At this point, the strain– stress distributions are shown in Fig. 17(i) and (ii) for maximum Fr = 38.21 N. The strain component T is symmetric about the center of the laminate which reaches its maximum value 4.0E06 at the bottom and minimum value 4.0E06 at the top. Also, LT is symmetric about the center of the laminate which reaches its max-

1e−05

800000

8e−06

600000

6e−06

400000

4e−06 200000

2e−06 0

0

−2e−06

−200000

−4e−06

−400000

−6e−06 −600000

−8e−06 −1e−05 −0.0006

−0.0004

−0.0002

0

(i)

0.0002

0.0004

0.0006

−800000 −0.0006

−0.0004

−0.0002

0

0.0002

0.0004

(ii)

Fig. 19. Stress–strain distribution through-the-laminate-thickness (@D of CD) for max. Resultant moment, Mr = 0.122 N m @ t = 0.0002 s.

0.0006

14

H. Pollayi et al. / Composite Structures 100 (2013) 1–16

0

2500 Clockwise

−2000

Reserve Factor

Reserve Factor

2000 1500 1000 500

−4000 −6000 −8000 −10000 −12000 Anti−Clockwise

0

1

2

3

4

5

6

7

8

9 10 11 12 13 14 15 16

−14000

1

2

3

4

5

6

7

8

9 10 11 12 13 14 15 16

Ply No.

Ply No.

(i)

(ii)

Fig. 20. Reserve factor (R) in each ply (RLaminate = + 63.29 (clockwise), 51.02 (anti-clockwise)).

imum value 2.0E06 at the top and minimum value 2.0E06 at the bottom. The stress component rT reaches its maximum value 320,000 N/m2 at about h = 0.00045 m and minimum value 280,000 N/m2 at about h = 0.00045 m. At the same ‘D’, the reserve factor, R, in each ply of the laminate ‘CD’ for the applied clockwise and with the reversal of the applied constant angular velocity are shown in Fig. 18(i) and (ii), respectively. In this case, ply-3 of bar-3 is predicted to fail first for maximum Fr and its R = +124.65. For maximum Fr and with reversal of the applied constant angular velocity, ply-14 of bar-3 is predicted to fail first and its R = 146.43. Fig. 19(i) and (ii) shows the strain–stress distributions for maximum resultant moment, Mr = 0.122 N m. It is observed that the strain components L and LT follows the same pattern from bottom to top of the laminate and also the stress components rL and rT follows the same pattern from bottom to top of the laminate. Fig. 20(i) and (ii) shows the reserve factor, R, in each ply at ‘D’ of the laminate ‘CD’ for the applied constant angular velocity, X = p rad/s in clockwise and anti-clockwise. In this case, ply-14 of bar-3 is predicted to fail first for maximum Mr and its R = +163.29. For maximum Mr and with reversal of the applied constant angular velocity, ply-3 of bar-3 is predicted to fail first and its R = 51.02. Fr and Mr are evaluated at the other joints and along the reference curves of all the component-laminates and tabulated in the Table 5. It is observed that the mechanism as a whole fails due to maximum resultant moment, Mr = 0.47 N m which occurs at t = 0.005 s. Work that is underway will be reported in later papers.

Appendix A. Non-linear stiffness for initially curved and pretwisted strip(s)

8. Conclusion

SNL 15 ¼

The component-laminate load-carrying capacity is evaluated using the Tsai–Wu–Hahn failure criterion for various layups of the three bars in the four-bar mechanism. The reserve factor, R of each ply in the component-laminate is calculated by using the maximum resultant force, Frmax, and the maximum resultant moment, Mrmax occurring at different time steps along the reference curves of all the component-bars and at the joints of the mechanism. The non-linear through-the-thickness stress and strain distributions are reported at different points. This illustrates the failure of each component-laminate and the mechanism as a whole. Tsai–Wu–Hahn failure criterion is used to predict both the first-ply-failure and the failure of the mechanism as a whole. Also, it is observed that the layer orientation influences the stress distribution in the component-laminates and it is also very important to determine the location of highly stressed regions in the component-laminate amongst all the layers of each component of the system. Work under progress will be reported in later papers.

The following are the analytical expressions for the coefficients of the 1-D non-linear stiffness matrix, SNL.

 @F 1 @ @U 1D ¼ @ c11 @ c11 @ c11  1 4 2 4 4 bA11 720D222 þ j22 b B212 þ k2 b B212 þ 2k2 b B212 j2 =D222 ¼ 720  @F 1 @ @U 1D ¼0 SNL Þ ¼ 12 ¼ @ð2c12 Þ @ð2c12 @ c11  @F 1 @ @U 1D SNL ¼0 ¼ 13 ¼ @ð2c13 Þ @ð2c13 Þ @ c11

SNL 11 ¼

 @F 1 @ @U 1D ¼ @ j1 @ j1 @ c11  1 2 2 ¼ b 120960B16 D222  5040b A11 k1 D222  5040b A11 j1 D222 60480

SNL 14 ¼

4

4

4

168b A11 B12 k1 j2 D22  168j1 b A11 B12 j2 D22 þ 168j22 b A11 B12 D26 6

4

2 6

þj22 b A11 B212 k1  168k2 b A11 B12 k1 D22 þ 4k2 b B212 k1 A11 2 4

4

4

þ168k2 b B12 A11 D26  168k2 b B12 A11 j1 D22 þ 8k2 b B212 j2 k1 A11 4 þ336k2 b B12 j2 A11 D26 =D222  @F 1 @ @U 1D ¼ @ j2 @ j2 @ c11  1 4 4 ¼ b 30240B11 D222  84b A11 B12 k1 j1 D22  42b A11 B12 j21 D22 30240 4

4

4

126b A11 B12 D12 j22  84b A11 B212 j2 c11 þ 168b A11 B12 j2 j1 D26 6

þb A11 B212 j2 j 4 þ4k2 b B212

2 4 1 k1  42k2 b A11 B12 D12

4  84k2 b B212 A11 11

4

SNL 22



j1 k1 A11 þ 168k2 b B12 j1 A11 D26  168k2 b B12 A11 D12 j2 =D222

 @F 1 @ @U 1D ¼0 ¼ @ j3 @ j3 @ c11  @F 2 @ @U 1D ¼ ¼ ¼0 @ c11 @ c11 @ð2c12 Þ  @F 2 @ @U 1D ¼ ¼ ¼ 0 ðBut we use SNL 22 ! 1 @ð2c12 Þ @ð2c12 Þ @ð2c12 Þ

SNL 16 ¼ SNL 21

c

4

in numerical calculationsÞ  @F 2 @ @U 1D ¼0 SNL ¼ ¼ 23 @ð2c13 Þ @ð2c13 Þ @ð2c12 Þ

15

H. Pollayi et al. / Composite Structures 100 (2013) 1–16





@F 2 @ @U 1D ¼ ¼0 @ j1 @ j1 @ð2c12 Þ  @F 2 @ @U 1D ¼0 ¼ ¼ @ j2 @ j2 @ð2c12 Þ  @F 2 @ @U 1D ¼ ¼ ¼0 @ j3 @ j3 @ð2c12 Þ  @F 3 @ @U 1D ¼ ¼ ¼0 @ c11 @ c11 @ð2c13 Þ  @F 3 @ @U 1D ¼ ¼ ¼0 @ð2c12 Þ @ð2c12 Þ @ð2c13 Þ  @F 3 @ @U 1D ¼ ¼ ¼ 0 ðBut we use SNL 33 ! 1 @ð2c13 Þ @ð2c13 Þ @ð2c13 Þ

SNL 24 ¼ SNL 25 SNL 26 SNL 31 SNL 32 SNL 33

 @M2 @ @U 1D ¼ @ c11 @ c11 @ j2  1 4 4 ¼ b 30240B11 D222  84b A11 B12 k1 j1 D22  42b A11 B12 j21 D22 30240

SNL 51 ¼

4

6

þb A11 B212 j2 j 4 þ4k2 b B212

SNL 53

4

4

168b A11 B12 k1 j2 D22  168j1 b A11 B12 j2 D22 þ 168j 4



j1 k1 A11 þ 168k2 b B12 j1 A11 D26  168k2 b B12 A11 D12 j2 =D222

 @M 2 @ @U 1D ¼0 ¼ @ð2c12 Þ @ð2c12 Þ @ j2  @M 2 @ @U 1D ¼0 ¼ ¼ @ð2c13 Þ @ð2c13 Þ @ j2

2 4

4

6

4

6

4

4

20160b A11 D12 k1 j1 D22  30240b A11 B12 D12 j2 c11

2 6

6

4

þ360b A11 D12 j2 j1 B12 k1 þ 20160c11 b A11 B12 j1 D26 6

4

þ120c11 b A11 B212 j1 k1  10080j21 b A11 D12 D22 6

8

2

8

2

240j21 b A11 B12 D26 k1 þ 40j21 b A211 k1 D22 þ 9j21 b A11 k1 B212

 @M1 @ @U 1D ¼0 ¼ ¼ @ð2c12 Þ @ð2c12 Þ @ j1  @M1 @ @U 1D ¼0 ¼ ¼ @ð2c13 Þ @ð2c13 Þ @ j1

@U 1D @ j1



4

þ120j23 b A211 D22  20160j21 b A11 D226  3628800D11 D222

2 4 2 b A11 B12 D26

4



c

4

5040k2 b A11 D212  5040c211 b A11 B212  360j23 b A11 B212

þ168k2 b B12 A11 D26  168k2 b B12 A11 j1 D22 4 4 þ8k2 b B212 j2 k1 A11 þ 336k2 b B12 j2 A11 D26 =D222

@M 1 @ SNL 44 ¼ @ j1 ¼ @ j1

4  84k2 b B212 A11 11

 @M 2 @U 1D @ SNL 55 ¼ @ j2 ¼ @ j2 @ j2  4 4 1 b 60480b A11 D12 j2 j1 D26  30240b A11 D212 j22 ¼  3628800

þj22 b A11 B212 k1  168k2 b A11 B12 k1 D22 þ 4k2 b B212 k1 A11

SNL 43

4

2 4 1 k1  42k2 b A11 B12 D12

SNL 52 ¼

 @M 1 @ @U 1D ¼ @ c11 @ c11 @ j1  1 2 2 ¼ b 120960B16 D222  5040b A11 k1 D222  5040b A11 j1 D222 60480

SNL 42

4

126b A11 B12 D12 j22  84b A11 B212 j2 c11 þ 168b A11 B12 j2 j1 D26

SNL 41 ¼

2 4



@M 1 @ @U 1D ¼ @ j3 @ j3 @ j1  1 7 ¼ k2 b A11 3B212 þ A11 D22 j2 =D222 30240

SNL 46 ¼

in numerical calculationsÞ  @F 3 @ @U 1D ¼0 SNL ¼ 34 ¼ @ j1 @ j1 @ð2c13 Þ  @F 3 @ @U 1D SNL ¼ ¼0 35 ¼ @ j2 @ j2 @ð2c13 Þ  @F 3 @ @U 1D ¼0 ¼ SNL 36 ¼ @ j3 @ j3 @ð2c13 Þ

6



6

4

þ960k2 b D12 j1 k1 A11 B12 þ 40320k2 b D12 j1 A11 D26 4 4 30240k2 b D212 A11 j2  20160k2 b D12 A11 B12 c11 =D222

 4 4 2 6 1 b 10080k2 b A11 B12 c11 D22 þ 40320j2 b D22 A11 k1 D26 14515200D66 D222  960k2 b k1 A11 D26 B12 ¼  3628800

2

4

4

2 4

2

2

4

2

þ1209600b B16 k1 D222  45360b A11 k1 D222  302400c11 b A11 D222  68040b A11 j21 D222 þ 1814400b 4

j1 D222 B16  302400j2 b2 D222 B11

j1 D222 A11 k1  720j2 b6 D22 A11 B12 k21 6 4 4 6 8 2 2160b A11 j2 j1 D22 B12 k1 þ 60480b A11 j2 j1 D22 D26  10080j22 b A11 D12 D22  240j22 b A11 B12 D26 k1 þ 40j22 b A211 k1 D22 20160j22 b A11 D226  20160k2 b A11 D226  10080b A11 D12 j2 c11 D22  136080b 8

2

6 2

4

4

2 8 2

6

þ9j22 b A11 k1 B212  720k2 b k1 D22 A11 B12 þ 40320k2 b k1 D22 A11 D26 þ 40k2 b k1 A211 D22  1080k2 b þ60480k2 b

4

@M 1 @ SNL 45 ¼ @ j2 ¼ @ j2

4

4

2 2 A11 D26

j1 D22 A11 D26  10080k2 b j2 A11 D12 D22  40320k2 b j



@U 1D @ j1



4

1 bð20160b ¼  1814400

4

6

 1920k2 b

4

þ15120j

2 6 1 b A11 D22 B12 k1

 540j

6



j2 k1 A11 D26 B12 þ 80k2 b8 j2 k21 A211 D22 =D222

j1 D22 A11 k1 D26  151200b2 B11 k1 D222  151200b2 j1 D222 B11 þ 3628800D16 D222 4

j1 D22 A11 B12 k21  5040j1 b4 A11 D22 B12 c11 4 þ 90b A11 D12 j22 B12 k1 þ 15120b A11 D12 j22 D26 þ 10080c11 b A11 B12 j2 D26

10080b A11 D12 k1 j2 D22  5040b A11 B12 k1 c11 D22  5040k2 b A11 D12 k1 D22  360b 2 4 1 b A11 D22 D26

j1 D22 k1 A11 B12

6

4

6

4

4

6

8

2

þ60c11 b A11 B212 j2 k1  20160b A11 j2 j1 D226  10080b A11 j2 j1 D12 D22  240b A11 j2 j1 B12 D26 k1 þ 40b A211 j2 j1 k1 D22 8

2

2 6

2 4

6

4

þ9b A11 j2 j1 k1 B212 þ 120k2 b D12 k1 A11 B12 þ 5040k2 b D12 A11 D26 þ 480k2 b D12 j2 k1 A11 B12 þ 20160k2 b D12 j2 A11 D26 4

4

4

þ240k2 b B212 c11 k1 A11 þ 10080k2 b B12 c11 A11 D26  5040k2 b 960k2 b

6

j1 A11 D12 D22  20160k2 b4 j1 A11 D226

j1 k1 A11 D26 B12 þ 40k2 b8 j1 k21 A211 D22  180k2 b6 A11 j3 B212 þ 60k2 b6 A211 j3 D22 Þ=D222

16

H. Pollayi et al. / Composite Structures 100 (2013) 1–16

@M2 @ SNL 54 ¼ @ j1 ¼ @ j1



@U 1D @ j2



 4 2 2 1 b 20160b j1 D22 A11 k1 D26  151200b B11 k1 D222 151200b j1 D222 B11 þ 3628800D16 D222 ¼  1814400

4

4

4

10080b A11 D12 k1 j2 D22  5040b A11 B12 k1 c11 D22  5040k2 b A11 D12 k1 D22  360b 4

6

6

6

j1 D22 A11 B12 k21  5040j1 b4 A11 D22 B12 c11

4

4

þ15120j21 b A11 D22 D26  540j21 b A11 D22 B12 k1 þ 90b A11 D12 j22 B12 k1 þ 15120b A11 D12 j22 D26 þ 10080c11 b A11 B12 j2 D26 6

4

4

6

8

2

þ60c11 b A11 B212 j2 k1  20160b A11 j2 j1 D226  10080b A11 j2 j1 D12 D22  240b A11 j2 j1 B12 D26 k1 þ 40b A211 j2 j1 k1 D22 8

2

2 6

2 4

6

4

þ9b A11 j2 j1 k1 B212 þ 120k2 b D12 k1 A11 B12 þ 5040k2 b D12 A11 D26 þ 480k2 b D12 j2 k1 A11 B12 þ 20160k2 b D12 j2 A11 D26 4

4

j1 A11 D12 D22  20160k2 b4 j1 A11 D226 6  180k2 b A11 j3 B212 þ 60k2 b A211 j3 D22 =D222

þ240k2 b B212 c11 k1 A11 þ 10080k2 b B12 c11 A11 D26  5040k2 b 6

960k2 b

8

2 2 1 k1 A11 D22

j1 k1 A11 D26 B12 þ 40k2 b j

4

6

 @M 2 @ @U 1D ¼ @ j3 @ j3 @ j2 1 5 2 2 ¼ b A11 ð84B12 j3 D22  6b j2 j3 B212 þ 2b j2 j3 A11 D22 30240

SNL 56 ¼

 3k2 b

2

j1 B212 þ k2 b2 j1 A11 D22 Þ=D222

 @M3 @ @U 1D ¼0 ¼ @ c11 @ c11 @ j3  @M 3 @ @U 1D ¼0 ¼ ¼ @ð2c12 Þ @ð2c12 Þ @ j3  @M 3 @ @U 1D ¼0 ¼ ¼ @ð2c13 Þ @ð2c13 Þ @ j3

SNL 61 ¼ SNL 62 SNL 63

 @M 3 @ @U 1D ¼ @ j1 @ j1 @ j3  1 7 k2 b A11 3B212 þ A11 D22 j2 =D222 ¼ 30240  @M 3 @ @U 1D ¼ ¼ @ j2 @ j2 @ j3  1 5 2 2 b A11 84B12 j3 D22  6b j2 j3 B212 þ 2b j2 j3 A11 D22 ¼ 30240

SNL 64 ¼

SNL 65

3k2 b

SNL 66 ¼

2

j1 B212 þ k2 b2 j1 A11 D22 =D222

 @M 3 @ @U 1D ¼ @ j3 @ j3 @ j3  1 3 2 4 b A11 2520D222 þ 84b B12 j2 D22  3j22 b B212 ¼ 30240 4

2

2 4

2 4

þj22 b A11 D22 þ 84k2 b B12 D22  3k2 b B212 þ k2 b A11 D22 =D222 References [1] Tsai SW, Wu EM. A general theory of strength for anisotropic materials. Journal of Composite Materials 1971;5:58–80. [2] Tsai SW, Hahn HT. Introduction to Composite Materials. Technomic Publishing Company, Lancaster, Pennsylvania.1980.

[3] Berdichevsky VL. Variational-asymptotic method of constructing the nonlinear shell theory. Theory of shells. Amsterdam: North-Holland Publishing Company; 1980. pp. 137–161. [4] Danielson DA, Hodges DH. Nonlinear beam kinematics by decomposition of the rotation tensor. Journal of Applied Mechanics 1987;54:258–62. [5] Cesnik CES, Hodges DH. VABS: a new concept for composite rotor blade crosssectional modeling. Journal of the American Helicopter Society 1997;42:27–38. [6] Popescu B, Hodges DH, Cesnik CES. Obliqueness effects in asymptotic crosssectional analysis of composite beams. Computer and Structures 2000;76:533–43. [7] Rajagopal A, Hodges DH. Analytical beam theory for the in-plane deformation of a composite strip with in-plane curvature. Composite Structures 2012;94:3793–8. [8] Yu W. Variational asymptotic modeling of composite dimensionally reducible structures, PhD thesis. Aerospace Engineering, Georgia Institute of Technology, Georgia, USA; 2002. [9] Yu W, Ho JC, Hodges DH. Variational asymptotic beam sectional analysis – an updated version. International Journal of Engineering Science 2012;59:40–64. [10] Popescu B, Hodges DH. Asymptotic treatment of the trapeze effect in finite element cross-sectional analysis of composite beams. International Journal of Non-Linear Mechanics 1999;34:709–21. [11] Hodges DH. Nonlinear composite beam theory, vol. 213 of Progress in Astronautics and Aeronautics. American Institute of Aeronautics and Astronautics, Inc., 1801 Alexander Bell Drive, Reston, Virginia; 2006. p. 20191–4344. [12] Hodges DH, Harursampath D, Volovoi VV, Cesnik CES. Non-classical effects in non-linear analysis of pretwisted anisotropic strips. International Journal of Non-Linear Mechanics 1999;34:259–77. [13] Harursampath D. Non-classical non-linear effects in thin-walled composite beams, PhD thesis. Aerospace Engineering, Georgia Institute of Technology, Georgia, USA; 1998. [14] Bauchau OA, Kang NK. A multi-body formulation for helicopter structural dynamic analysis. Journal of the American Helicopter Society 1993;38:3–14. [15] Bauchau OA. Computational schemes for flexible, nonlinear multi-body systems. Multibody System Dynamics 1998;2:169–225. [16] Bauchau OA, Hodges DH. Analysis of nonlinear multi-body systems with elastic couplings. Multibody System Dynamics 1999;4:168–88. [17] Yu W. Efficient high-fidelity simulation of multibody systems with composite dimensionally reducible components. Journal of the American Helicopter Society 2007;52:49–57. [18] Pollayi H, Harursampath D. Geometrically non-linear dynamics of composite four-bar mechanisms. International Journal of Non-Linear Mechanics 2012;47:837–50. [19] Armanios EA, Makeev A, Hooke D. Finite displacement analysis of laminated composite strips with extension-twist coupling. Journal of Aerospace Engineering 1996;9:80–91. [20] Kosmatka JB. Extension-bend-twist coupling behavior of nonhomogeneous anisotropic beams with initial twist. AIAA Journal 1992;30:519–27. [21] Winckler SI. Hygrothermally curvature stable laminates with tension-torsion coupling. Journal of the American Helicopter Society 1985;31:56–8.