On the foundations of anisotropic interior beam theories

Composites Part B 87 (2016) 299e310

Contents lists available at ScienceDirect

Composites Part B journal homepage: www.elsevier.com/locate/compositesb

On the foundations of anisotropic interior beam theories Anssi T. Karttunen*, Raimo von Hertzen Department of Applied Mechanics, School of Engineering, Aalto University, FI-00076 Aalto, Finland

a r t i c l e i n f o

a b s t r a c t

Article history: Received 28 May 2015 Received in revised form 24 October 2015 Accepted 28 October 2015 Available online 6 November 2015

This study has two main objectives. First, we use the Airy stress function to derive an exact general interior solution for an anisotropic two-dimensional (2D) plane beam. Second, we cast the solution into the conventional form of 1D beam theories to clarify some basic concepts related to anisotropic interior beams. The derived general solution provides the exact third-order interior kinematic description for the plane beam and includes the Levinson/Reddyekinematics as a special case. By applying the Clapeyron's theorem, we show that the stresses acting as surface tractions on the lateral end surfaces of the interior beam need to be taken into account in all energy-based considerations related to the interior beam in order to avoid artificial end effects. Exact 1D interior beam equations are formed from the general 2D solution. Finally, we develop an exact interior beam finite element based on the general solution. With full anisotropic coupling, the stiffness matrix of the element becomes initially asymmetric due to the interior nature of the plane beam. By redefining the generalized nodal axial forces of the element, the stiffness matrix takes a symmetric form. © 2015 Elsevier Ltd. All rights reserved.

Keywords: B. Anisotropy B. Elasticity C. Analytical modelling C. Finite element analysis (FEA)

1. Introduction Efficient use of anisotropic composite materials in mechanical design requires thorough understanding of anisotropic elasticity and accurate tools of analysis. Motivated by this, we study here a two-dimensional (2D) plane beam with anisotropic coupling effects within the framework of 2D linear elasticity. There are two well-known complex variable formulations for 2D linearly elastic anisotropic plane problems, the Lekhnitskii and Stroh formalisms [1e3]. When it comes to 2D interior plane beam problems, where the end effects are neglected by virtue of the Saint Venant's principle, a number of more straightforward polynomialbased stress function approaches can be found in the literature, e.g. Refs. [4e6]. It is common for these classical polynomial approaches that a solution is generated only for one problem at a time. In this study, we first provide a more versatile method for 2D beams which is based on a general interior solution derived using the Airy stress function. Our approach is different from the generalization of Silverman's method [4] by Ding et al. [7] in the way that, rather than constructing stress functions for each problem separately, the focus is on solving the cross-sectional force and moment resultants along a beam for each case in the same way as in conventional beam theories, while the core of the used stress function is always the same. * Corresponding author. Tel.: þ358 (0) 405272713. E-mail addresses: anssi.karttunen@aalto.fi (A.T. Karttunen), raimo.von.hertzen@ aalto.fi (R. von Hertzen). http://dx.doi.org/10.1016/j.compositesb.2015.10.026 1359-8368/© 2015 Elsevier Ltd. All rights reserved.

There is a number of ingenious theories meant specifically for anisotropic beams, e.g. Refs. [8e13]. For a historical survey on the topic, see the book by Hodges [14]. The mentioned anisotropic beam theories are typically applicable to a wider variety of practical problems than the linearly elastic beam with a rectangular cross-section studied in this paper. However, it is commonplace in the referenced treatments to employ engineering assumptions and/or to rely on a quite heavy theoretical machinery. Therefore, it is often difficult to see the underlying structure of the developments in a fully explicit (and assumption-free) form in order to state something fundamental on anisotropic beams in general. Thus, for the remainder of this study, we view the derived 2D solution for the anisotropic plane beam as a conventional one-dimensional (1D) beam theory in order to clarify certain concepts related to anisotropic interior beams. As the starting point for the 1D considerations, the general solution provides the exact third-order interior kinematic description for the plane beam. By “third-order kinematics” we mean that the displacement components defined at the central axis of the beam are expanded by third-order polynomials throughout the height of the beam. Many approximate beam and plate theories, which are used also in association with anisotropic composite materials, are based on similar, but assumed, displacement fields. For reviews on assumed third-order kinematics, see the works by Jemielita [15] and Reddy [16,17]. The exact interior kinematic description derived from the general solution in this paper includes, as a special case, the displacement

300

A.T. Karttunen, R. von Hertzen / Composites Part B 87 (2016) 299e310

field which is used to derive the widely known Levinson and ReddyeBickford beam theories [18e20]. The long-standing belief in the literature is that the Levinson beam theory is “variationally inconsistent”, that is to say, it cannot be derived using the principle of virtual displacements. In contrast to this, we showed in a very recent study that the Levinson theory is actually variationally consistent [21]. The variational formulation was carried out by taking into account the fact that the stresses of the beam act as surface tractions on the lateral end surfaces of the Levinson beam. On the other hand, it was shown that the boundary layer behavior of the ReddyeBickford beam is artificial. In the present study, the methodology utilized for the isotropic Levinson theory in Ref. [21] is further elucidated within the 2D interior framework of anisotropic elasticity. The rest of the paper is organized in the following way. In Section 2, we formulate the anisotropic plane beam problem to which the solution is then given in terms of a stress function. The strains are calculated from the stresses under plane stress conditions and the exact 2D interior displacement field is obtained by integrating the strains. Calculation examples are presented. In Section 3, three kinematic variables defined at the central axis of the plane beam are formed from the 2D interior displacement field. Using these variables, the third-order kinematic description for the beam is given. Clapeyron's theorem is employed to facilitate energy-based considerations. Variational and vectorial approaches for interior theories are discussed and finally 1D beam equations for the 2D plane beam are obtained by a direct method. In Section 4, an exact interior beam finite element is developed both by a force-based method and from the total potential energy of the anisotropic plane beam. Conclusions are presented in Section 5.

2.1. Problem formulation A 2D linearly elastic homogeneous anisotropic plane beam under a uniform pressure p is shown in Fig. 1. We have chosen the uniform load as a representative load for our developments. The beam has a rectangular cross-section of constant thickness t and the length and height of the beam are L and h, respectively. The load resultants N, M and Q stand for the axial force, bending moment and shear force, respectively, and act at an arbitrary cross-section of the beam. These cross-sectional load resultants are calculated from

NðxÞ ¼ t

Zh=2 sx dy ; MðxÞ ¼ t

h=2

Zh=2 sx ydy ; Q ðxÞ ¼ t

h=2

sy ðx; h=2Þ ¼ p ;

sy ðx; h=2Þ ¼ 0 ;

txy ðx; ±h=2Þ ¼ 0:

sx ¼

v2 J ; vy2

sy ¼

v2 J ; vx2

v2 J ; txy ¼  vxvy

(3)

which satisfy the two-dimensional equilibrium equations. To ensure compatibility, it is required that the stress function satisfies the governing equation [23].

s22

v4 J v4 J v4 J v4 J  2s26 3 þ 2ðs12 þ s66 Þ 2 2  2s16 4 vx vx vy vx vy vxvy3

(4)

where sij are the elastic compliances, see Eqs. (15)e(17). The solution to the interior plane beam problem is obtained by finding a solution of Eq. (4) that satisfies the stress boundary conditions (2) of the beam.

2.2. Interior stresses By starting from the general polynomial of the fifth degree and adapting a solution procedure outlined by Barber [22, Chap. 5] to solve the polynomial coefficients, the stress function satisfying the stress boundary conditions (2) can be found as

txy dy ;

 Jðx; yÞ ¼ c1 y2 þ c2 y3  c3

h=2

(2)

The boundary conditions are satisfied in a strong (pointwise) sense on the upper and lower surfaces, but at the beam ends the tractions are specified only through the load resultants and, thus, the boundary conditions are imposed only in a weak sense [22]. The replacement of the true boundary conditions at the beam ends by the statically equivalent weak boundary conditions (load resultants) means that the exponentially decaying end effects of the anisotropic plane beam are neglected by virtue of the Saint Venant's principle and only the interior solution of the beam is under consideration. The interior solution represents actually a beam section with fully-developed interior stresses which has been cut off from a complete beam far enough from the real lateral boundaries at which the true boundary conditions could be set. Using the Airy stress function J(x,y), the stresses of the plane beam are obtained from the equations

v4 J þ s11 4 ¼ 0; vy

2. General interior approach for a plane beam

Zh=2

which can be used to impose force and moment boundary conditions at x ¼ ±L/2. The boundary conditions on the upper and lower surfaces of the plane beam are

(1)

3 2 1 s16 4 h xy  xy3  y 4 2 s11

 þ Jq ;

(5)

where c1, c2 and c3 are constant coefficients and the part that depends on the nature of the applied load is

# " s216 5 q 2 2 Jq ¼ 8 2 y  5x ðh  yÞðh þ 2yÞ 240I s11 h   i q 2s12 y5 þ 5s16 h2  2y2 xy2 þ s66 y5 ;  120Is11

(6)

where q ¼ pt is the uniform line load and I ¼ th3/12 is the second moment of the cross-sectional area. The stresses are calculated from Eq. (3), after which the load resultants are obtained from Eq. (1) and can be written as Fig. 1. 2D anisotropic plane beam under a constant uniform pressure. The load resultants act at an arbitrary cross-section of the beam. The positive directions of the load resultants are defined for later use in calculation examples and finite element developments.

NðxÞ ¼ 2Ac1 þ 6Ic3

s16 ; s11

(7)

A.T. Karttunen, R. von Hertzen / Composites Part B 87 (2016) 299e310

MðxÞ¼6Iðc2 þc3 xÞþ

 i q h 20s211 x2 þh2 4s216 2s11 s12 þs11 s66 ; 2 40s11 (8)

Q ðxÞ ¼ 6Ic3 þ qx;

(9)

where A ¼ ht is the area of the cross-section. The load resultants satisfy the global equilibrium equations

vNðxÞ ¼0; vx

vMðxÞ ¼ Q ðxÞ ; vx

vQ ðxÞ ¼ q: vx

(10)

s16 ½Q ðxÞ  qx: s11

vUx vUy þ ; vy vx

(20)

where Ux(x,y) and Uy(x,y) are the displacements in the directions of x and y, respectively. The displacements are integrated from Eqs. (18) and (19), and the resulting arbitrary integration functions are solved by substituting the calculated displacements into (20). This standard procedure is described, for example, in Ref. [23]. As the result, we get for the general 2D interior displacements the expressions q

Note that due to the anisotropic coupling term s16 in Eq. (7), the stretching of the beam is not uncoupled from the shearing. Furthermore, using Eqs. (7) and (9), we can write

NðxÞ ¼ 2Ac1 þ

gxy ¼

301

(11)

Ux ¼ 2c1 ðs11 x þ s16 yÞ þ 3c2 ð2s11 x þ s16 yÞy þ Cy þ D1 þ Ux " !  2  s2 h þ y2 x þ c3 3s11 x2 y þ 2 16  s12 y3 þ 3s16 s11 4   # 3 þ s66 h2  y2 y ; 4

(21) By solving Eqs. (7)e(9) for c1, c2 and c3, the stresses (3) calculated using the stress function (5) are obtained in the form

sx ¼

  NðxÞ MðxÞy Q ðxÞ s16 y2 þ  1  12 2 A I A s11 h !  s216 2s12 þ s66 3qy qy3  ; þ 4 2 10A s11 6I s11

(12)

 q  3 h þ 3h2 y  4y3 ; sy ¼  24I

(13)

  2 2 s16 h  4y2 : qy txy ¼ Q ðxÞ þ 3 s11 8I

(14)

Indeed, when compared to an isotropic case, a marked effect due to the anisotropy is that the shear force Q(x) appears in the axial stress (12). The above interior stress distributions are valid for any plane beam with a rectangular cross-section. To obtain stress distributions for other loads, one needs to update only the load-dependent part (6) of the stress function (5).

2.3. Interior strains and displacements Under plane stress, that is, for a beam with a narrow rectangular cross-section, the stressestrain relations read

εx ¼ s11 sx þ s12 sy þ s16 txy ;

(15)

εy ¼ s12 sx þ s22 sy þ s26 txy ;

(16)

gxy ¼ s16 sx þ s26 sy þ s66 txy ;

(17)

where εx , εy and gxy are the axial normal strain, transverse normal strain and transverse shear strain, respectively. To obtain the interior strains, we substitute Eqs. (12)e(14) into Eqs. (15)e(17). The strain-displacement relations are

  q Uy ¼ 2c1 s12 y þ 3c2 s12 y2  s11 x2  Cx þ D2 þ Uy     s 3 þ c3 3s12 xy2  s11 x3 þ 2s12 16 y3 þ s26 h2  y2 y ; 4 s11 (22) where Uxq and Uyq depend on the particular load at hand and are given in Appendix A. The integration constants C, D1 and D2 constitute the degrees of freedom of the anisotropic plane beam as a rigid body. The constant C relates to a small clockwise rotation about the origin and D1 and D2 correspond to translations in the directions of x and y, respectively. Note that as soon as the crosssectional load resultants are available, the coefficients c1, c2 and c3 can be calculated from Eqs. (7)e(9). 2.4. Calculation examples In this Section, we study four calculation examples, which are presented in Fig. 2. We note that the interior stresses are of main interest here. In the calculation of the displacements, we follow the customary practice, that is, the displacement conditions to solve the constants C, D1 and D2 from Eqs. (21) and (22) are chosen so as to imitate the actual conditions at the ends of the beams. To impose true boundary conditions, we would need to model the boundary layer behavior using PapkovicheFadle eigenfunctions [22]. 2.4.1. End-loaded cantilever beam First we consider a cantilever beam subjected to an end-load P at x ¼ L/2 as shown in Fig. 2(a). The structure is statically determinate and the load resultants using the positive directions of Fig. 1 are

PðL þ 2xÞ ; Q ¼ P; N ¼ 0; M ¼  2 from which it follows by Eqs. (7)e(9) that

c1 ¼

vUx ; εx ¼ vx

(18)

vUy ; εy ¼ vy

(19)

P s16 ; 2A s11

PL ; c2 ¼  12I

P c3 ¼  : 6I

By substituting c1, c2 and c3 into the displacements (21) and (22) and choosing the conditions Ux ðL=2; 0Þ ¼ Uy ðL=2; 0Þ ¼ vUx =vyðL=2; 0Þ ¼ 0, the rigid body constants become

302

A.T. Karttunen, R. von Hertzen / Composites Part B 87 (2016) 299e310

Fig. 3. The relative difference between the Airy stresses and those from a 2D FE cantilever model at y ¼ 0 with the elastic compliances s11 ¼ 4.49,1012, s22 ¼ 8.13,1012, s12 ¼ 9.48,1015, s16 ¼ 5.84,1013, s26 ¼ 2.13,1013 and s66 ¼ 1.31,1011 m2/N. The relative difference is calculated from (sx,FE  sx,Airy)/sx,Airy  100.

By substituting c1, c2 and c3 into the displacements (21) and (22), and choosing the conditions Ux(L/2,0) ¼ Uy(±L/2,0) ¼ 0, we find the rigid body constants to be Fig. 2. Calculation examples. (a) Cantilever beam subjected to an end-load. (b) Simplysupported beam under a uniform load. (c) Cantilever beam under a uniform load. (d) Propped cantilever beam under a uniform load.

 PLs16 PL  2 2 5L s11  2h2 s216 þ 3h2 s11 s66 ; ; D2 ¼ 48Is11 4A   P 9L2 s211  2h2 s216 þ 3h2 s11 s66 : C¼ 24Is11 D1 ¼

The strains are obtained from Eqs. (18)e(20) and the stresses by inverting Eqs. (15)e(17), or directly from Eqs. (12)e(14). The axial normal and transverse shear stresses are

  i P h2 4y2 P h 6s11 ðLþ2xÞys16 h2 12y2 ; txy ¼ : sx ¼ 12Is11 8I

Fig. 3 shows the relative difference at the centerline (y ¼ 0) between the above stresses and those calculated from a 2D plane stress FE model with the boundary conditions Ux(L/2,y) ≡ Uy(L/ 2,y) ≡ 0 using the parameter values L ¼ 1.6 m, h ¼ 0.2 m, t ¼ 0.04 m, P ¼ 1 kN. The shear stress reflects the “rule of thumb” used in isotropic cases, that is, the interior solution is a good approximation when the axial distance to an end of the beam is at least equal to the height of the beam. However, the axial normal stress given by the 2D FE model starts to deviate from the interior solution well before this limit. It is known that the boundary layer can be thicker in anisotropic than in isotropic cases [24]. 2.4.2. Simply-supported beam As another example, let us consider a simply-supported beam under a constant uniform load q according to the configuration presented in Fig. 2(b). The axial force, bending moment and shear force along the beam are given by

NðxÞ ¼ 0 ;

  2 x L2  ; MðxÞ ¼ q 2 8

We obtain from Eqs. (7)e(9)

c1 ¼ c3 ¼ 0 ; c2 ¼

Q ðxÞ ¼ qx:

h i q h2 s11 ð2s12 þ s66 Þ  4h2 s216  5L2 s211 2 240Is11

C ¼ 0;

D1 ¼

qL ð4hs12  Ls16 Þ; 16A

i qL2 h 25L2 s211  16h2 s216 þ 6h2 s11 ð3s12 þ 4s66 Þ : D2 ¼  1920Is11 The axial normal and transverse shear stresses are

   i qy h 15s211 L2  4x2 þ 4s216 3h2  20y2 sx ¼  2 120Is11 h     q s12 40y3  6h2 y þ 10s16 h2  12y2 x  120Is11   i þ s66 20y2  3h2 y ;

txy ¼

  q ð3s11 x þ 2s16 yÞ h2  4y2 : 24Is11

Fig. 4 shows how the above shear stress behaves near the center of the beam. The parameter values are the same as in the previous example with the exceptions of P ¼ 0 and q ¼ 40 kN/m. Results from a 2D plane stress FE model are also presented. The shear stress distribution is cubic due to anisotropy over the height of the beam at the center but starts to turn quickly into a parabola outside the cross-section x ¼ 0.

2.4.3. Cantilever beam under uniform load As a further example, we study a cantilever beam under a uniform load, see Fig. 2(c). It is easy to see that this case can also be solved by the superposition of the two cases presented above by choosing P ¼ qL/2 for the end-loaded cantilever. By superposing the stresses and calculating the load resultants (1), we obtain for the uniformly-loaded cantilever

N¼0;

M¼

qðL þ 2xÞ2 ; 8

Q¼

qðL þ 2xÞ : 2

An instructive discussion on the use of superposition in 2D plane beam problems can be found in Ref. [25].

A.T. Karttunen, R. von Hertzen / Composites Part B 87 (2016) 299e310

fðxÞ ≡

303

 vUx 3  ðx;0Þ ¼ 2c1 s16 þ 6c2 s11 x þ c3 4x2 s11 þ h2 s66 þ C 4 vy q h 3 2 þ 4x s11 þ 2h2 s216 x  h2 s11 ð3s12 x 24Is11 i þ hs26 Þ ; (25)

respectively. Using these kinematic central axis variables, we can present the displacements (21) and (22) of the plane stress state in the form

 vuy s16 4y3  fþ vx s11 3h2 ! 2 3 vf y s12 s16   2 vx 6 s11 s11

Ux ¼ ux þ yf  þ Fig. 4. Shear stress distributions near the center of the uniformly-loaded simplysupported beam.

2.4.4. Propped cantilever under uniform load As the final example, we examine the propped cantilever under a uniform load shown in Fig. 2(d). The load resultants can be written as 2

N¼0;

M¼

qðL þ 2xÞ FðL þ 2xÞ ;  2 8

Q¼

qðL þ 2xÞ  F; 2

y2 s16 2 s11



v2 f þ UxA ; vx2

(26)

  s12 vux y vfx þ s11 vx 2 vx   i v2 f  y h 2 8y  3h2 s12 s16 þ 3h2  4y2 s11 s26 þ UyA ; þ 2 vx2 24s11

Uy ¼ uy þ y

(27) where the functions UxA and UyA depend on the applied load and are

where F is the unknown vertical reaction force at the propped end. The coefficients c1, c2 and c3 are obtained from Eqs. (7)e(9) after which C, D1, D2 and F are calculated using the conditions Ux(L/ 2,0) ¼ Uy(L/2,0) ¼ vUx/vy(L/2,0) ¼ Uy(L/2,0) ¼ 0. The reaction force is found as

 q 30L3 s211  12h2 Ls216  h2 s11 ð9Ls12 þ 10hs26  18Ls66 Þ

  ; F¼ 10 8L2 s211 þ h2 3s11 s66  2s216 after which the calculation of the stresses is straightforward. In conclusion, the general interior solution enables the solution of many 2D plane beam problems in a unified and effortless manner.

given in Appendix A. The displacements (26) and (27) represent the exact third-order interior kinematics of a linearly elastic homogeneous anisotropic beam under a uniform load that has a narrow rectangular cross-section. In the isotropic case, the elastic compliances become

s16 ¼ s26 ¼ 0 ; s11 ¼ s22 ¼ 1=E ; s12 ¼ n=E ; s66 ¼ 1=G; (28) where E and G are the Young's modulus and shear modulus, respectively, and n is the Poisson ratio. Assuming isotropy, the displacements (26) and (27) simplify to

Ux ¼ ux þ yf  3. 2D solution in 1D form 3.1. Third-order kinematics of the beam

Uy ¼ uy  ny

The general 2D interior solution studied in Section 2 can be presented in terms of three kinematic variables derived from the displacements (21) and (22) at the central axis of the beam so that the variables depend only on x. We obtain for the axial displacement and the transverse deflection of the central axis, and for the clockwise positive rotation of the cross-section at the central axis the expressions

3 qx ðs x ux ðxÞ ≡ Ux ðx; 0Þ ¼ 2c1 s11 x þ h2 c3 s16 x þ D1 þ 4 4A 16  2hs12 Þ;

vux vx

(23)

uy ðxÞ ≡ Uy ðx; 0Þ ¼ 3c2 s11 x2  c3 s11 x3  Cx þ D2 i qx2 h 2 2 2s11 x þ 4h2 s216  3h2 s11 ðs12 þ s66 Þ ;  48Is11 (24)

  vuy 4y3 ny3 v2 f þ f þ ; vx 6 vx2 3h2

vux ny2 vf þ UyI ;  vx 2 vx

(29)

(30)

where

UyI ¼

  i qy h ð2h þ 3yÞ n2  1 h2 þ 2y3 ð1 þ 2nÞ : 48EI

(31)

We see that the Poisson effect (lateral contraction/expansion) is fully included in the displacements (29) and (30). To derive the 2D displacement field in terms of the central axis variables for isotropic material behavior, one may also consider the Marguerre function for the case of plane stress and the Love strain function for plane strain instead of the Airy stress function, see, e.g. Ref. [26]. If we neglect the load term and the Poisson effect in Eqs. (29) and (30), the 2D displacement field is exactly of the same form as in the Levinson and ReddyeBickford beam theories [18e20]. These theories are formulated by assuming that the central axis variables ux(x), uy(x) and f(x) do not initially take the polynomial forms (23)e(25), but are arbitrary, sufficiently smooth functions. Next, to discuss interior theories in general and the well-known

304

A.T. Karttunen, R. von Hertzen / Composites Part B 87 (2016) 299e310

discrepancies between the Levinson and ReddyeBickford theories, we turn to the energetic aspects of the plane beam.

3.2. Applying Clapeyron's theorem The strain energy of the 2D anisotropic plane beam and the external work due to the uniform load are given by

U¼

1 2

Z



sx εx þ sy εy þ txy gxy dV ;

V

(32)

ZL=2 Wq ¼ 

qUy ðx; h=2Þdx;

for the Levinson interior beam theory was carried out by accounting for the external virtual work due to the stresses on the end surfaces of the interior beam [21]. The resulting beam equations were exactly the same as those in the case of the vectorial derivation by Levinson [18]. The variational formulation for the Levinson beam theory given in Ref. [21] did not actually include the axial displacement ux(x). With this in mind, let us complement that formulation by the following consideration, which illustrates how consistent equilibrium equations for interior theories can be obtained when employing the principle of virtual displacements. Example e variational treatment of a rod with Poisson effect. We study an isotropic interior rod in pure tension so that M(x) ≡ Q(x) ≡ uy(x) ≡ f(x) ≡ 0 and q ¼ 0. The displacements (29) and (30) are

L=2

respectively. While the contributions (32) to the total potential energy of the beam are apparent, the following consideration may not be that obvious at first. We recall from the problem definition in Section 2.1 that the interior solution part represents merely a beam section with fully-developed interior stresses which has been cut off from a complete beam so that the boundary layer is not modeled. Thus, we obtain for the work due to the stresses on the lateral end surfaces of the interior beam

Zh=2 Ws ¼ t

sx ðL=2; yÞUx ðL=2; yÞdy  t h=2

Uy ¼ ny

sx ðL=2; yÞUx ðL=2; yÞdy þ t

vux ; vx

(35)

where the central axis variable according to Eqs. (23) and (28) is

ux ¼ Ux ðx; 0Þ ¼

2c1 x þ D1 : E

(36)

In the following, we want to carry out a variational formulation which leads to the governing differential equation for the axial rod. Our aim is to consider under what conditions the variational procedure leads to the same result as the polynomial Airy stress function method. We assume next that ux(x) is a sufficiently smooth function, but otherwise arbitrary. We take as the starting point of the variational formulation the displacement field (35) of the rod

Uy ðx; yÞ ¼ nyu0x ;

Ux ðx; yÞ ¼ ux ;

Zh=2 h=2

Ux ¼ ux ;

(37)

where comma denotes differentiation with respect to x. The axial normal and transverse shear stresses are

(33)

Zh=2

sx ¼ Eεx ¼ Eu0x ;

txy ðL=2; yÞUy ðL=2; yÞdy  t h=2

txy ¼ Ggxy ¼ Gnyu00x :

(38)

The internal virtual work is

Zh=2

ZL=2 Z

txy ðL=2; yÞUy ðL=2; yÞdy:

dU ¼

h=2

L=2

¼

2U  Wq  Ws ¼ 0:

(34)

The above calculation shows that in static equilibrium the strain energy of the beam is equal to one-half of the work done by the surface tractions if they were to move (slowly) through their respective displacements. This is exactly according to the Clapeyron's theorem (see, e.g. Ref. [23]). The recent variational formulation

sx dεx þ txy dgxy dAdx ¼

ZL=2

 Ndu0x  nSdu00x dx

L=2

A

ZL=2 By substituting the polynomial expressions according to the general solution presented in Section 2 for sx, εx , sy, εy , txy, gxy, Ux and Uy into (32) and (33), and by handling the laborious calculations by the aid of Mathematica (see the online Supplementary File), we find that



L=2 ðN 0 dux þ nS0 dux Þdx þ Ndux  nSdu0x þ nS0 dux L=2 ;

L=2

(39) where the axial force and the higher-order load resultant are given by

Z N¼ A

sx dA ¼ EAu0x ;

Z S¼

nh2 GAu0x ; txy ydA ¼  12

(40)

A

respectively. The external virtual work due to the interior stresses on the lateral end surfaces calculated using (33), (37) and (38) results in

A.T. Karttunen, R. von Hertzen / Composites Part B 87 (2016) 299e310

L=2 dWs ¼ Ndux  nSdu0x L=2 :

(41)

εy ¼

By applying the principle of virtual displacements, dU ¼ dWs, the equilibrium equation and interior boundary conditions are found as

N 0 þ nS00 ¼ 0

and

L=2

½nS0 dux L=2 ¼ 0:

nS ð±L=2Þ ¼ 0:

(43)

n2 h2 0 0 00 00 000 GAu000 f ≡ nS0 ¼  x / f ¼ N ¼ EAux / f ¼ EAux : 12 (44) It follows that f00 ¼ b2f / f ¼ a1ebx þ a2ebx, where b2 ¼ 24(1 þ n)/ (nh)2. With the interior boundary conditions (43), that is, f(±L/ 2) ¼ 0, we obtain a1 ¼ a2 ¼ 0 / f ≡ 0. Therefore, the equilibrium equation in (42) simplifies to

  vuy 4y2 f þ vx h2 " !# s16 vf 4 vux 1 s16 v2 f þy 2 þ þ y þ gAxy ; 2 s11 vx2 s11 vx h vx 1

(48)

where the load-dependent parts εAx , εAy and gAxy are given in Appendix A. Using the strains (46)e(48) and by inverting the plane stress constitutive relations (15)e(17), we can calculate the load resultants (1) in terms of the central axis variables and substitute them into the equilibrium equations (10) to obtain 1D beam equations. For full anisotropic coupling this procedure leads to unwieldy expressions. For the sake of brevity, we consider here only the isotropic case, which can be easily compared to existing 1D beam theories. We apply Eq. (28) to the strains (46)e(48) and to the constitutive relations (15)e(17), after which the calculation of the load resultants (1) gives

N ¼ EA

vux qhn ;  2 vx

(49)

M ¼ EI

vf qh2  ð2 þ 5nÞ; vx 40

(50)

Q¼

  vuy 2 GA f þ : 3 vx

(51)

Substitution of the load resultants into the equilibrium equations (10) leads to

v2 u x ¼ 0; vx2

(52)

  vuy v2 f 2 GA f þ ¼ 0;  vx vx2 3

(53)

EA

EI

2 vf v2 uy GA þ 3 vx vx2

! ¼ q:

(54)

By uncoupling Eqs. (53) and (54), we obtain alternatively

3.3. 1D beam equations Instead of a variational formulation similar to the previous example, we may derive the 1D beam equations directly by using the global equilibrium equations (10). In terms of the polynomial central axis variables (23)e(25), the strains (18)e(20) calculated using the displacements (21) and (22) of the plane stress state can be written as

vux vf y2 s16 v2 f þ þy þ εAx ; vx vx 2 s11 vx2

gxy ¼

(45)

which may be compared to the first of (10). The solution to (45) is given by ux in (36), that is, by the Airy elasticity solution. Note that if the function f above would not vanish, the differential equation for ux would be of fourth order and its solution would also contain exponential terms. We see that the higher-order load resultant S has no physical meaning in the interior context. We conclude that with the inclusion of the lateral contraction and expansion (the Poisson effect) in the rod model, it becomes apparent that the virtual work due to the interior stresses on the lateral end surfaces of the interior rod is a crucial part of the variational formulation in order to avoid spurious exponential end effects and violation of the axial equilibrium of the rod. Similar conclusions were arrived at in the case of the Levinson beam theory, for details and further discussion, see Ref. [21].

εx ¼

(47) 

Next we use (40) and (42) and write

N 0 ¼ EAu00x ¼ 0;

  s12 vux vf þy vx s11 vx   i v2 f  1 h 2 h  4y2 s11 s26  h2  8y2 s12 s16 þ 2 þ εAy ; vx2 8s11

(42)

Note that the second term on the left-hand side of the equilibrium equation is due to the Poisson effect and that it violates the axial equilibrium condition if it is nonzero. In addition, according to the interior beam definition, the virtual displacement dux is free in the interior region, including the interior end surfaces. This means that the interior boundary conditions simplify to

0

305

(46)

EI

v3 f ¼ q; vx3

(55)

EI

v4 uy ¼ q: vx4

(56)

In summary, the general solution to Eqs. (52)e(54) is given by Eqs. (23)e(25) and (28), and the solution is expanded into the whole interior region of the 2D plane beam through Eqs. (29)e(31). These equations constitute an alternative representation of the elasticity solution presented in Section 2 in the isotropic case. We can see, for example, that the obtained beam equations (55) and (56) are

306

A.T. Karttunen, R. von Hertzen / Composites Part B 87 (2016) 299e310

exactly the same as those in the Timoshenko beam theory with a constant uniform distributed load. Eqs. (51) and (54) are the same as in the static Levinson beam theory. 1D beam presentations may be obtained for more complicated loads by updating the loaddependent part (6) of the stress function (5). The provided exact 1D interior presentation enables a more thorough study of the pros and cons of approximate interior beam theories, such as those of Levinson and Timoshenko, on the level of governing differential equations instead of simplistic (numerical) comparisons between specific 2D elasticity and 1D beam solutions. 4. Exact anisotropic beam finite element

where

4.1. Elasticity solution in terms of FE-degrees of freedom The general interior displacements (21) and (22) can be used as the basis for the derivation of an exact 1D beam finite element. Similarly to the 1D beam equations given in the previous section, the finite element equations constitute an alternative presentation of the exact 2D interior elasticity solution presented in Section 2. Fig. 5 presents the setting according to which the element is developed. Both nodes in Fig. 5 have three degrees of freedom. For nodes i ¼ 1,2, we have axial displacements ux,i, transverse displacements uy,i and rotations fi. By the aid of the central axis variables (23)e(25), we obtain for nodes 1 and 2 six equations

ux;1 ¼ ux ðL=2Þ ; ux;2 ¼ ux ðL=2Þ; uy;1 ¼ uy ðL=2Þ ; uy;2 ¼ uy ðL=2Þ; f1 ¼ fðL=2Þ ; f2 ¼ fðL=2Þ:

By substituting Eqs. (58) and (59) into Eqs. (21) and (22) we can write the 2D displacements in the form

Ux ðx; yÞ ¼ Nx u þ UxF ;

(61)

Uy ðx; yÞ ¼ Ny u þ UyF ;

(62)

u ¼ ux;1

 3h2 s16  2 uy;1  uy;2 þ Lðf1 þ f2 Þ 4DL

2L2 s11 þ 3h2 s66  ux;1  ux;2 2DL i q h 2 2L s11 s12 þ 3h2 ðs12 s66  s16 s26 Þ ; þ 4Dt h  i f1  f2 q L2 s211 þ h2 2s216  3s11 s12 ;  c2 ¼ 2 6Ls11 144Is11

uy;1

f1

ux;2

uy;2

f2

T

(63)

is the displacement vector. The shape functions Nx and Ny are given in Appendix B, as well as the load-dependent functions UxF and UyF . Therefore, once the nodal displacements are known, the 2D displacement field can be calculated by substituting them into Eqs. (61) and (62), after which the calculation of the 2D interior strains and stresses is straightforward.

4.2. Finite element equations To obtain the finite element equations, we substitute Eq. (58) into Eqs. (7)e(9) to calculate the load resultants at nodes i ¼ 1,2, with the notion that the positive directions are taken to be according to Fig. 1 so that

 2s  4s16  ux;1  ux;2  11 2 uy;1  uy;2 þ Lðf1 þ f2 Þ DL DL 2q þ ðs11 s26  s12 s16 Þ; Dt (58)

  i L2 s2 3 uy;1  uy;2 h 2 2 11 L s11 þ h2 s11 s66  s216 þ ðf1 þ f2 Þ C¼ DL 2D 

(60)

where



c3 ¼

  D ¼ 2L2 s211 þ 3h2 s11 s66  s216 :

(57)

We can solve the six unknowns c1, c2, c3, C, D1 and D2 from Eq. (57). This results in the following relations

c1 ¼

Fig. 5. Set-up according to which the exact interior beam finite element is developed.

qL2 s11 Ls11 s16  ux;1  ux;2  ðs11 s26  s12 s16 Þ; D 2Dt

ux;1 þ ux;2 qL2 s16  ; 2 16A  1  D2 ¼ 4 uy;1 þ uy;2 þ Lðf1  f2 Þ 8  i qL2 h 2 2  L s11 þ h2 6s11 s66  4s216 ; 384Is11

D1 ¼

N1 ¼ NðL=2Þ ; N2 ¼ NðL=2Þ ; Q1 ¼ Q ðL=2Þ ; Q2 ¼ Q ðL=2Þ ; M1 ¼ MðL=2Þ ; M2 ¼ MðL=2Þ:

The conventional presentation for the 1D beam element is obtained by writing Eq. (64) in matrix form. Before doing so, we also derive the FE equations from the total potential energy

P ¼ U  Wq  Ws :

(59)

(64)

(65)

The stresses on the end surfaces in Eq. (33) are written as given by Eqs. (12)e(14), where the load resultants are expressed as nodal forces according to Eq. (64). Then, by calculating (32) and (33) and by applying the principle of minimum total potential energy

A.T. Karttunen, R. von Hertzen / Composites Part B 87 (2016) 299e310

vP ¼0; vux;i

vP ¼0; vuy;i

vP ¼0 vfi

ði ¼ 1; 2Þ

stiffness matrix (68) becomes symmetric. It is also noteworthy that with full anisotropic coupling, every nodal load is associated with every nodal degree of freedom in the end surface work Ws. For example, N1 is not conjugate only to ux,1 in the conventional manner, but also to all the other degrees of freedom. If we were to neglect the interior nature of the plane beam and derive, in the absence of the distributed load, the stiffness matrix in the conventional way on the basis of the strain energy U, we would obtain a symmetric, but incorrect, stiffness matrix. In the isotropic case, the asymmetry of the stiffness matrix (68) vanishes along with the coupling term s16 and the equations simplify to

(66)

we obtain the finite element equilibrium equations. The forcebased method and the total potential energy approach result in the same equations, which can be written in the form

Ku ¼ f þ q:

307

(67)

In explicit form, with a ¼ 12I/D, we have

2

s2 A þ 16 aLs11 Ls11

6 6 6 6 6 2s 6  16 6 L 6 6 6 6 s16 6 6 K ¼ a6 ! 6 6 s216 A 6 þ 6 aLs11 Ls11 6 6 6 6 2s16 6 6 L 6 6 4 s16

f ¼ fN1

Q1

M1

N2

s2 A  þ 16 aLs11 Ls11

!

s  16 L

s  16 2

2s11 L

s11

2s16 L

2s  11 L

s11

D þ 6L2 s211 12Ls11

s16

s11

s16 L

s16 2

s2 A þ 16 aLs11 Ls11

s  16 L

2s  11 L

s11

2s  16 L

2s11 L

s11

D  6L2 s211  12Ls11

s16

s11

Q2

M2 gT ;

s16 L

3

7 7 7 7 7 7 s11 7 7 7 2 2 7 D  6L s11 7  7 12Ls11 7 7 7 7 s16 7 7 2 7 7 7 7 7 s11 7 7 2 2 7 D þ 6L s11 5 12Ls11

(68)

(69)  EA 1 L 1

9 8 2hL2 s211 s12  2h3 s12 s216  h3 s11 ðs16 s26  3s12 s66 Þ > > > > > > > > > > 2Ds11 > > > > > > > > > > 2 3 2 2 > > h s ð2hs  3Ls Þ þ 2L s  h s ð2hs  3Ls Þ > > 16 12 16 11 26 66 11 > > > > > > > > 2D > > > > # " > >   > > 2 2 2 3 3 > > > > > >  1 2h s16 þ 3h ð3s12  s66 Þ  60h Ls11 s26 þ 10L2 1 þ 6h s12 s16 > > > > > > 2 D DL s11 = < 120 s11 q¼q > > 2 2 3 2 3 > > ðs Þ 2hL s s  2h s s  h s s  3s s > > 12 12 11 16 26 12 66 11 16 > >  > > > > > > 2s11 D > > > > > > > > > > h2 s ð2hs þ 3Ls Þ  2L3 s2  h2 s ð2hs þ 3Ls Þ > > > > 16 12 16 11 26 66 11 > > > > > > > > 2D > > > " # > > >   > > > > 1 2h2 s2 2 3 3 > > 3h ð3s  s Þ 60h Ls s 6h s s > > 12 66 11 26 12 16 2 16 > > 1  þ þ þ 10L ; : 2 120 s D DL s11 11

s  16 2

1 1



ux;1 ux;2



 ¼

N1 N2

 þ

qnh 2



 1 ; 1

9 38 6L 12 6L > > uy;1 > > = 2 2 7< ð4 þ FÞL 6L ð2  FÞL 7 f1 5 6 12 6L u > > > : y;2 > ; f2 ð2  FÞL2 6L ð4 þ FÞL2 9 8 L > > > > > > 8 9 > > 2 2 > > > > 10L  3h ð2 þ 5nÞ Q > > > > 1 > > > >  < = q< = M1 60 ¼ þ ; Q > > > > 2> L > > : 2> ; > > > > M2 > > 2 2 > > > > > 10L  3h ð2 þ 5nÞ > ; : 60

(71)

2

(70)

The remarkable feature of Eqs. (67) is the asymmetry of the stiffness matrix (68). The reciprocal theorem (e.g. Ref. [22]) requires the stiffness matrix to be symmetric. In light of this, the asymmetry of the stiffness matrix can be easily remedied by redefining the generalized axial forces in Eq. (64) as N1 ≡ 2N(L/2) and N2 ≡ 2N(L/2). In other words, if we multiply the first and fourth row of the linear system of equations (67) by a factor of two, the

12 6 6L EI 6 ð1 þ FÞL3 4 12 6L

(72)

where F ¼ 3h2(1 þ n)/L2. Eqs. (71) and (72) can be written concisely as

Kr ur ¼ f r þ qr ;

(73)

308

A.T. Karttunen, R. von Hertzen / Composites Part B 87 (2016) 299e310

Kb ub ¼ f b þ qb ;

(74)

where Kr and Kb are the rod and beam element stiffness matrices, respectively, and ur and ub are the rod and beam nodal displacement vectors, respectively. The force vectors are fr and qr for the rod element, and fb and qb for the beam element, respectively. Note that the stiffness matrix Kb is equal to the stiffness matrix of the Levinson beam element and the Timoshenko beam element when the Timoshenko shear coefficient has the value of 2/3 [21,27]. The presented finite element formalism for the plane beam extends the applicability of the 2D elasticity solution.

(A.2) The load-dependent terms in Eqs. (26) and (27) are

UxA 5. Conclusions In this paper, a general interior elasticity solution for a 2D linearly elastic anisotropic plane beam under a uniform load was derived. The exact third-order interior kinematic description for the beam was obtained from the solution. It was shown that in an energy-based formulation, which is founded exclusively on the interior kinematics of the beam, one has to take into account the fact that the interior stresses do work on the lateral end surfaces of the interior beam. A number of higher-order beam theories can be found in the literature that are based exclusively on interior kinematics. However, these higher-order constructions are incomplete as interior beam theories because they lack the aforementioned work due to the stresses at the beam ends. From a different point of view, the higher-order beam theories attempt to model both the interior part and the boundary layer of a beam simultaneously with the same (interior) kinematic description. In search of an accurate beam theory that would better serve this end, the kinematic description should also include information on the true boundary layer behavior. Adapting the methodology presented in this paper, a kinematic description for the boundary layer may be constructed using PapkovicheFadle eigenfunctions in association with the Airy stress function [22].

    q h 3s12 h2 þ 4y2 x2 þ s22 2y3  2h3  3h2 y y 48I i qy s h     i 16 s26 h2  2y2 y  2s12 h2  4y2 x  2s11 x4 þ 24I s11  q s66  2 q h 3h s11 x2  2s12 y4 þ 4s12 s216 y4 þ 48I s11 24Is211     i  2s11 s212 y4 þ h2 s216 x2 þ s211 s26 3h2  4y2 xy :

q

Uy ¼

"  s3   qy2 s12 s16  15h2 þ 8hy  24y2  16 h2  4y2 ¼ 24Is11 6 s11 #   s s   11 26 2 2 2 2 9h þ 4hy  6y ; þ s16 s66 h  3y  3 (A.3)

UyA

"  s s   s212  3 qy 2h þ 3h2 y  4y3  11 22 2h3 þ 3h2 y  2y3 ¼ 24Is11 2 2 #    s12 s216  2  y h  4y2 þ s16 s26 y h2  2y2  s12 s66 y3 : s11 (A.4)

The load-dependent terms in Eqs. (46)e(48) are given by

εAx ¼

εAy

i qy3 h 2 2s16  s11 ðs12 þ s66 Þ ; 6Is11

(A.5)

"     q ¼ s2 h3 þ 3h2 y  8y3  s11 s22 h3 þ 3h2 y  4y3 24Is11 12 #    s12 s216  2 2 y h  8y2 þ 2s16 s26 y h2  4y2  4s12 s66 y3 ; s11 (A.6) "

Appendix A. Load-dependent expressions

gAxy ¼

The load-dependent terms in Eqs. (21) and (22) read

q

Ux ¼

   q  3 3 q 3 4 2 3 2 2 s s þ h x2 x y þ s y  s s xy þ 12y 66 16 11 16 11 48I 6Is211    q s16 h  2 þ y s12 h  4y2 y þ 2s16 h2 þ 4y2 x 24I s11   i  q h  3 s12 h þ 3h2 y þ 4y3 x þ s66 h2  3y2 y  24I   i þ s26 h3 þ 3h2 y  2y3 y ; (A.1)

   s3  qy 3h2 þ 4hy  8y2  2 16 h2  8y2 s s 24Is11 12 16 s11 #      s11 s26 3h2 þ 4hy  4y2 þ 2s16 s66 h2  6y2 : (A.7)

Appendix B. Shape functions and load functions The shape functions in Eqs. (61) and (62) are In addition, in Eqs. (61) and (62) we have the load-dependent terms

A.T. Karttunen, R. von Hertzen / Composites Part B 87 (2016) 299e310

309

9 8     > > 2 3 3 3 2 2 2 2 > > > > 2L ðL  2xÞs þ 16y s  3ðL  2xÞs  h s s ðL  2xÞs þ 8y ðyðs þ s Þ  3xs Þ 2ðL þ 2xÞys  s 3h 16 66 11 16 16 12 66 16 > > 11 16 11 > > > > > > > > > > 2DLs 11 > > > > > >  > >   2 > > 2 2 2 2 > > > > y 3 L  4x  8y s  4ys ð3xs  yðs þ s ÞÞ s 11 16 12 66 > > 11 16 > > > > > > > > DL > > > > > >        > >  > > 3 2 3 2 2 2 2 2 2 2 > > > y LðL þ 6xÞðL  2xÞs11  3h ys16 þ s11 s16 3h ðL  2xÞ  8Ly s16 þ 3h ys66 þ s11 4Ly s12 þ 2LðL  6xÞys16 þ 4Ly  3h ðL  2xÞ s66 > > > > > > > > > = < 2DLs 11 T ; Nx ¼     > > > > > > 2L2 ðL þ 2xÞs311  16y3 s316 þ 3ðL þ 2xÞs211 2ðL  2xÞys16 þ h2 s66  s11 s16 3h2 ðL þ 2xÞs16  8y2 ðyðs12 þ s66 Þ  3xs16 Þ > > > > > > > > > > 2DLs11 > > > > > > > >  > >   > > 2 2 2 2 2 > > y 3 L  4x  8y s  4ys ð3xs  yðs þ s ÞÞ s > > 11 16 12 66 > > 11 16 > > > >  > > > > DL > > > > >    >     > >  > > 3 2 3 2 2 2 2 2 2 2 > > > y LðL  6xÞðL þ 2xÞs11 þ 3h ys16 þ s11 s16 3h ðL þ 2xÞ  8Ly s16  3h ys66 þ s11 4Ly s12  2LðL þ 6xÞys16 þ 4Ly  3h ðL þ 2xÞ s66 > > > > > > > > > ; : 2DLs11 (B.1)

9 8 > >         > > > > > > 8y3 s12 s216  s211 2L2 ys12  x L2  4x2 s16 þ ys11 3h2  4y2 s16 s26 þ 3s12 4xys16  h2 s66 > > > > > > > > > > > > DLs > > 11 > > > >         > > > > > > 2 2 2 3 2 2 2 3 2 > > 2 3h s  s 24xy s 2ðL  2xÞ ðL þ xÞs þ s y  8y  3h ðL  2xÞs s þ 6h y  8y  3h ðL  2xÞs > 16 12 16 11 12 26 66 > 11 > > > > > > > > > > 2DL > > > > > >       > >  > > 3 2 2 2 2 2 2 > > 2 2 2 s 12xys s  L x s þ y 3h  8y s  ys þ 3h  4y 4x > > 12 16 11 12 26 > > 11 L  4x þ 4y ðs12 =s11 Þ > > > > þ = < 8L 2D T Ny ¼ :         > > > > > > 8y3 s12 s216  s211 2L2 ys12  x L2  4x2 s16 þ ys11 3h2  4y2 s16 s26 þ 3s12 4xys16  h2 s66 > > > > > >  > > > > DLs11 > > > > > > >       >  > > > > 2 2 2 3 2 2 2 3 2 > > 2 3h s þ s 24xy s 2ðL þ 2xÞ ðL  xÞs  s y  8y þ 3h ðL þ 2xÞs s þ 6h y  8y þ 3h ðL þ 2xÞs > 16 12 16 11 12 26 66 > 11 > > > > > > > > > > 2DL > > > > > >       > >  > > 3 2 2 2 2 2 2 > > 2 2 2 > > s 12xys s  L x s þ y 3h  8y s  ys þ 3h  4y 4x 12 16 11 12 26 > > 11 L  4x þ 4y ðs Þ =s > > 12 11 > > þ  > > > > ; : 8L 2D (B.2)

UxF ¼

n        q 24h2 y2 h2  4y2 s516  4h2 y2 s11 s316 15h2 þ 16hy  24y2 s12 þ 48xys16 þ 6 2h2  7y2 s66 2 192DIs11  h           16L2 x L2  4x2 ys511 þ s211 s16 32h3 y3 s212 þ 3h4 L2  4x2  4h2 L2 þ 36x2 y2 þ 64L2 y4 s216 þ 24h2 y2 h2  3y2 s266      i h  þ 8h2 y2 s16 9h2 þ 8hy  6y2 s26 þ 36xys66 þ 4h2 y2 s12 24xðy  hÞs16 þ 15h2 þ 8hy  24y2 s66  2s411 32L2 xy3 s12             i   þ L2 h2 L2  4x2 þ 4 L2  12x2 y2 s16 þ 4y 3h3 L2  4x2 þ 6h2 L2 y  4L2 y3 s26 þ x 3h2 L2  4x2 þ 8L2 y2 s66 h            þ s311 8xy 3h2 L2  4x2 þ 16L2 y2 s216  8h2 y2 s66 9h2 þ 4hy  6y2 s26 þ 12xys66 þ s16 96h3 xy2 s26 þ 3h4 4x2  L2         io ; þ 4h2 L2 þ 36x2 y2  48L2 y4 s66  8ys12 3h3 4x2  L2  5h2 L2 y þ 8L2 y3 s16 þ 4h2 y2 ðhs26 þ 3xs66 Þ (B.3)

310

UyF ¼

A.T. Karttunen, R. von Hertzen / Composites Part B 87 (2016) 299e310

   2  48h2 y2 h2  4y2 s12 s416  2 L3  4Lx2 s511

q 384DIs211 h     þ s411 16 L2 y 12x2  L2 ys12 þ 2y3  2h3         3h2 y s22 þ x h3 L2  4x2 þ 6h2 L2 y  8L2 y3 s26   i h   3h2 5L4  24L2 x2 þ 16x4 s66 þ 4h2 s11 s216 2y2 12y2      9h2  16hy s212 þ 3s16 h2 4x2  L2 s16        4y2 h2  2y2 s26 þ 12ys12 2x h2  4y2 s16  i h   þ y 5y2  h2 s66 þ 2s211 12h2 ys212 2h3 þ 3h2 y       4y3 s66  8hxys16 þ 3h2 s16 8y2 h2  2y2 s26 s66      þ s16 4y 2h3 þ 3h2 y  2y3 s22 þ 8xy 4y2  3h2 s26       þ 4ys12 y 16L2 y2  h2 36x2 þ 5h2 L2  4x2 s66      þ L2 s216  6h2 y3 s266  12h2 s16 h h2  2y2 s26 þ x h2  i h   þ s311 16L2 y 2h3 þ 3h2 y  4y3 s212  4y2 s66    þ h2 11L4  56L2 x2 þ 48x4 s216 þ 32L2 y2 h2       2y2 s16 s26 þ 2h2 8h 3h2 y  4y3 s226 þ 12y 2y3       2h3  3h2 y s22 þ 2x 3h2  4y2 s26 s66  9h2 L2        4x2 s266 þ 8s12 2x h3 4x2  L2  4h2 L2 y      þ 16L2 y3 s16 þ y2 24h3 xs26 þ 3h2 12x2  L2  i :  4L2 y2 s66 (B.4)

Appendix A. Supplementary data Supplementary data related to this article can be found at http:// dx.doi.org/10.1016/j.compositesb.2015.10.026.

References [1] Lekhnitskii SG. Theory of elasticity of an anisotropic body. San Francisco: HoldeneDay; 1963. [2] Stroh AN. Dislocations and cracks in anisotropic elasticity. Philos Mag 1958;3(30):625e46. [3] Stroh AN. Steadyestate problems in anisotropic elasticity. J Math Phys 1962;41(2):77e103. [4] Silverman IK. Orthotropic beams under polynomial loads. J Eng Mech Div ASCE 1964;90(5):293e319. [5] Hashin Z. Plane anisotropic beams. J App Mech 1967;34(2):257e63. [6] Ding HJ, Huang DJ, Wang HM. Analytical solution for fixed-fixed anisotropic beam subjected to uniform load. Appl Math Mech 2006;27(10):1305e10. [7] Ding HJ, Huang DJ, Chen WQ. Elasticity solutions for plane anisotropic functionally graded beams. Int J Solids Struct 2007;44(1):176e96. [8] Giavotto V, Borri M, Mantegazza P, Ghiringhelli G, Carmaschi V, Maffioli GC, et al. Anisotropic beam theory and applications. Comput Struct 1983;16(1): 403e13. [9] Bauchau OA. A beam theory for anisotropic materials. J Appl Mech 1985;52(2): 416e22. [10] Berdichevsky V, Armanios E, Badir A. Theory of anisotropic thin-walled closed-cross-section beams. Compos Eng 1992;2(5):411e32. [11] Murakami H, Reissner E, Yamakawa J. Anisotropic beam theories with shear deformation. J Appl Mech 1996;63(3):660e8. [12] Yu W, Hodges DH, Volovoi VV, Fuchs ED. A generalized Vlasov theory for composite beams. Thin Wall Struct 2005;43(9):1493e511. [13] Yu W, Hodges DH, Ho JC. Variational asymptotic beam sectional analysisean updated version. Int J Eng Sci 2012;59:40e64. [14] Hodges DH. Nonlinear composite beam theory. Reston: US: American Institute of Aeronautics and Astronautics; 2006. [15] Jemielita G. On kinematical assumptions of refined theories of plates e a survey. J Appl Mech 1990;57(4):1088e91. [16] Reddy JN. A general nonlinear 3rdeorder theory of plates with moderate thickness. Int J Nonlin Mech 1990;25(6):677e86. [17] Reddy JN. Mechanics of laminated composite plates and shells: theory and analysis. 2nd ed. Boca Raton, Florida: CRC Press; 2003. [18] Levinson M. A new rectangular beam theory. J Sound Vib 1981;74(1): 81e7. [19] Bickford W. A consistent higher order beam theory. Dev Theor Appl Mech 1982;11:137e50. [20] Reddy JN. A simple highereorder theory for laminated composite plates. J Appl Mech 1984;51(4):745e52. [21] Karttunen AT, von Hertzen R. Variational formulation of the static Levinson beam theory. Mech Res Commun 2015;66:15e9. [22] Barber JR. Elasticity. 3rd ed. New York: Springer; 2010. [23] Sadd MH. Elasticity e theory, applications and numerics. 3rd ed. Oxford: Academic Press; 2014. [24] Horgan CO. Recent developments concerning SainteVenant's principle: an update. Appl Mech Rev 1989;42(11):295e303. [25] Doghri I. Mechanics of deformable solids. Berlin: Springer; 2000. [26] Soutas-Little RW. Elasticity. New York: Dover; 1973. [27] Kosmatka JB. An improved twoenode finite element for stability and natural frequencies of axialeloaded Timoshenko beams. Comput Struct 1995;57(1): 141e9.