A composite beam element with through the thickness capabilities based on global-local superposition

Accepted Manuscript A composite beam element with through the thickness capabilities based on global-local superposition André Lima, Alfredo Faria PII: DOI: Reference:

S0263-8223(17)32968-9 https://doi.org/10.1016/j.compstruct.2017.11.051 COST 9117

To appear in:

Composite Structures

Received Date: Revised Date: Accepted Date:

12 September 2017 13 November 2017 20 November 2017

Please cite this article as: Lima, A., Faria, A., A composite beam element with through the thickness capabilities based on global-local superposition, Composite Structures (2017), doi: https://doi.org/10.1016/j.compstruct. 2017.11.051

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

A composite beam element with through the thickness capabilities based on global-local superposition André Limaa,*, Alfredo Fariaa a

Department of Mechanical Engineering, Instituto Tecnológico de Aeronáutica (ITA), Praça Marechal Eduardo Gomes, 50, São José dos Campos, SP 12228-900, Brazil E-mail adresses: [email protected] (A. Lima), [email protected] (A. Faria). Abstract: A beam element is proposed that captures through the thickness effects in composite laminated beams, namely, transverse shear and normal stresses and strains. Stress continuity along the thickness is inherently enforced leading to a system of algebraic equations that is solved in the element level, permitting independence between the number of layers and the number of degrees of freedom, with all of them possessing a clear physical significance. Global-local superposition is performed in the thickness direction, where a cubic global displacement field, that guarantees imposition of the boundary conditions at the top and bottom surfaces of the beam, is combined with a layerwise linear local displacement distribution that assures zig-zag behavior of the stresses and displacements. The element behavior for different length-to-thickness ratios is assessed and compared to the analytical elasticity solution, as well as a commercial finite element alternative. Key-words: composite beam; through the thickness effects; zig-zag; global-local superposition. Funding: This work was funded by the Brazilian agencies CAPES grant 1642839; CNPq grant 300886/2013-6.

Introduction The application of multilayered structures has increased in many fields of engineering, accompanied by a burst in the number of theories and models to represent their behavior. The main reason driving that burst was the inability of classical theories of beams, plates and shells [1–4] to reproduce piecewise continuous displacement fields and transverse shear stress fields in the thickness direction [5]. Those theories also neglect transverse normal stresses, that are fundamental in failure analysis of multilayered structures (e.g., delamination). This piecewise form of displacement fields and transverse stresses is often described in the open literature as zig-zag (ZZ), and the theories that describe this behavior, as well as interlaminar continuity (IC), are called Zig-Zag theories [6]. Those two effects are summarized by the acronym C0z requirements, meaning that displacements and transverse stresses must be C0 continuous functions in the thickness direction. *

Corresponding author.

1

A historical review of ZZ theories was made by Carrera [6] who identified Lekhnitskii [8] as the first contributor to the theory of multilayered structures. Despite its great potential, Carrera observed that Lekhnitskii's contribution was systematically forgotten, except for the work by Ren [7,8]. Moreover, in his original work, Lekhnitskii did not fully complied with the Koiter’s recommendation that reads: “A refinement of Love’s first approximation theory is indeed meaningless, in general, unless the effects of transverse shear and normal stresses are taken into account at the same time” [9]. Koiter’s recommendation must not be overseen. From a mathematical point of view, for a bi-dimensional (xz) problem, neglecting transverse shear cripples the equilibrium equation in the x-direction, and neglecting transverse normal stresses cripples the equilibrium equation in the z-direction. Complete and comprehensive reviews on laminated composite plate and beam models were written by Sayyad and Ghugal [10,11] with a massive number of references of the different models. Nevertheless, in both works it was highlighted that theories considering the effect of transverse shear strain and transverse normal strain and stress are rarely found in the literature and need more attention. Also, thermo-mechanical analyses would greatly benefit from the consideration of transverse stresses and strains since there is usually substantial mismatch between the in-plane and transverse coefficients of thermal expansion in composite laminates. In an overall review of laminate theories based on displacement hypothesis Liu and Li [5] present the limitations and physical consequences of adopting different hypothesis and assumptions. If displacements are assumed smooth and continuous through the thickness, double-valued interlaminar stresses are predicted in the laminate interfaces, usually overestimating the stresses, even for higher-order theories, leading to over-dimensioned structures. In addition, this leads to inaccurate prediction of interlaminar failure and delamination, which are of fundamental importance to composite structures. In trying to circumvent some of the limitations just mentioned, Li and Liu [12] proposed the use of models based on a global-local superposition of displacement fields, enabling C0 continuous displacement and stress fields in the thickness direction. Zhen et al. [13] applied the concept of global-local superposition to build finite elements including transverse normal strain for linear elastic composite plates, and later Zhen and Li [14] included thermo-mechanical behavior. In these works, a superposition of cubic global and local fields was adopted for in-plane displacements and special assumptions needed to be made to ensure C0z continuity and still have enough equations to solve the problem, but the transverse normal displacement is assumed constant through the thickness [13] or it varies exclusively with temperature [14]. From a practical point of view, through the thickness stresses play decisive roles in many engineering applications. Composite laminates fabricated through curing cycles are prone to 2

fabrication induced residual stresses and severe distortions [15,16]. Such effects can only be analytically described by finite elements that capture through the thickness effects, typically requiring the use of very large models built on solid elements [17]. Moreover, problems involving contact require application of forces on the bottom and top surfaces of laminates that may come to touch other surfaces. These forces inevitably end up inducing large transverse stresses and strains. The element formulation developed in the present work comes to fulfill Koiter’s recommendation, to meet the suggestions proposed by Sayyad and Ghugal and to provide new capabilities to commercial finite element codes that may have new elements added to their suite of options. The formulated element includes a number of innovative features: piecewise continuous displacements, transverse normal strain/stress, transverse shear strain/stress and nonhomogeneous boundary conditions. The proposed element assumes a linear local displacement field within each layer, for both in-plane and transverse displacements, and continuity is inherently enforced without restraining the solution. This assumption is more intuitive and proves to be sufficient to build an efficient element, but still, the element is restricted to small displacement gradients, a fundamental hypothesis of linear elasticity. However, extensions to nonlinear application are straightforward since the degrees of freedom involved consist of either translations or rotations. The element behavior is assessed for different length-to-thickness ratios and compared to the analytical elasticity solution by Pagano [18] for a metallic sandwich of isotropic materials and for a [0°/90°/0°] composite laminate. Additionally, a cantilever beam problem with a tearing load at the free tip, whose analytical solution is not available, is solved and a comparison to a commercial 2D element mesh is established, revealing the suitability of the element for practical applications in relevant composite structures.

Constitutive relations Figure 1 illustrates a general kth lamina. Axis z is the through the thickness axis. Axis 1 is aligned with the fibre direction. Axis 2 is perpendicular to axis 1 and lies on the xy plane. Layer k has its orientation angle θk with respect to the structural coordinate system xyz. It is assumed that there are three planes of material symmetry (planes 12, 13 and 23). Hence, the constitutive flexibility matrix expressed in the material coordinate system as a function of engineering constants is

3

k k  1/ E1k −ν 21 / E2k −ν31 / E3k 0 0 0   k k  k k k 1/ E2 −ν32 / E3 0 0 0  −ν12 / E1 k −ν k / Ek −ν 23 / E2k 1/ E3k 0 0 0  Sk =  13 1  k 0 0 1/ G23 0 0   0  0 0 0 0 1/ G13k 0    0 0 0 0 1/ G12k   0

(1)

k k k k where E1k, E2k, E3k are Young moduli, ν12k , ν21 , ν13k , ν31 , ν23 , ν32k are Poisson ratios and G23 , G13k , G12k

are shear moduli. As required, Sk = (Sk)T. The constitutive stiffness matrix can be obtained by inversion of the constitutive flexibility matrix leading to C11k C12k C13k 0 0 0   k  k k C12 C22 C23 0 0 0  Ck Ck Ck 0 0 0  Ck = (Sk )−1 =  13 23 33 k   0 0 0 C44 0 0   0 0 0 0 Ck 0  55   k  0 0 0 0 0 C66

(2)

k k where, evidently, the components C11k , C12k , C13k , C22 , C23k , C33k , C44 , C55k and C66k can be computed

using the engineering constants present in Sk. The constitutive stiffness matrix in the structural xyz system is obtained as a function of Ck and the orientation angle θk. To shorten notation, we define s = sinθk and c = cosθk. Thus,

C11k  k C12 C13k k C = 0 0  k C16

C12k C13k

0

0

C22k C23k C23k C33k

0 0

0 0

0 0

0 0

C26k C36k

C44k C45k C45k C55k 0

0

C16k   C26k  C36k   0 0  C66k 

(3)

where the transformed C k matrix components are:

4

k 2 C13k = C13k c2 + C23 s

k 4 k C11k = C11k c4 + C22 s + 2s2c2 (C12k + 2C66 ) 2 2

4

4

2 2

k 2 C23k = C13k s2 + C23 c

C = (C + C )s c + C (c + s ) − 4s c C k 12

k 11

k 4 11

k 22

k 4 22

k 12

2 2

k 66

k C36k = (C13k − C23 )cs

C = C s + C c + 2s c (C + 2C ) k 22

k 12

k 66

k k C16k = (C11k − C12k )c3s + (C12k − C22 )cs3 − 2cs(c2 − s2 )C66 k k C26k = (C11k − C12k )cs3 + (C12k − C22 )c3s + 2cs(c2 − s2 )C66 k k C66k = (C11k + C22 − 2C12k )c2s2 + C66 (c2 − s2 )2

(4)

k C33k = C33 k 55

k 2 44

k 2 55

C =C s +C c

k k C45k = (C55 − C44 )cs k 2 k 2 C44k = C44 c + C55 s

Figure 2 depicts a composite beam structural axes (xyz) and the stacking sequence of N layers. Any given layer k has the bottom surface positioned at zk−1 and the top surface positioned at zk.

In the xyz coordinate system each layer has the following constitutive relations: σ xk  C11k  k  k σ y  C12  k = k σ z  C13 τ xyk  C16k  

k C12k C13k C16k ε x   k  τ yzk  C44k C45k γ yzk  C22k C23k C26k ε y  ,    k = k k k  k  C23k C33k C36k ε z  τ xz  C45 C55 γ xz   k C26k C36k C66k γ xy 

(5)

Particularization of Eq. (5) for the 2D plane must be done with care. A plane stress state in the xz plane is assumed such that σyk = τxyk = τyzk = 0. When these zero stresses are substituted in Eq. (5) they lead to the reduced material constitutive matrix, split in two parts

σ xk  Q11k Q13k ε xk  , τ xzk = Q55k γ xzk  k = k k  k  σ ε Q Q  z   13 33  z 

(6)

where the components of the reduced material constitutive matrix are: −1

Q11k Q13k  C11k C13k  C12k C16k  C22k C26k  C12k C23k  (C45k ) 2 k k = − , Q = C −  k          55 55 k k k k k k k k k C44k Q13 Q33  C13 C33  C23 C36  C26 C66  C16 C36 

(7)

From Eq. (6) it can be inferred that the transverse shear strain is linked to the transverse shear stress through Q55k . For a plane strain state in the xz plane εyk = γxyk = γyzk = 0. When these zero strains are substituted in Eq. (5) they lead to the reduced material constitutive matrix split in two parts:

σ xk  C11k C13k ε xk  , τ xzk = C55k γ xzk  k =  k k  k  σ z  C13 C33 ε z 

(8) 5

Kinematic relations Layerwise linear local distributions of displacements are assumed for both the axial displacement u and the transverse displacement w as shown in Fig. 3. Notice that the bottom (u0, w0) and top (uN, wN) displacements associated with the piecewise linear distributions are zero. Superimposed to the piecewise linear distributions are cubic functions defined for u and w based on four degrees of freedom for each: ub, φb, ut, φt for u and wb, ψ b, wt, ψ t for w. The subscribed b and t refer to the bottom and top surface of the laminate, respectively. The ub, ut, wb and wt, denote the translational displacement degrees of freedom and the φb, φt, ψ b, and ψ t denote the rotational degrees of freedom for the top and bottom surfaces. One can define a nondimensional through the thickness coordinate ζk for each layer k whose domain is −1 ≤ ζk ≤ 1. Two linear interpolation functions can then be defined as L0 (ζ k ) =

1− ζ k 1+ ζ k , L1(ζ k ) = 2 2

(9)

where ζk = 2z/hk and hk is the kth-layer thickness. The global cubic functions adopted are the Hermite cubics given by 1 h H0 (ζ ) = (ζ 3 − 3ζ + 2) , H1(ζ ) = (ζ 3 − ζ 2 − ζ +1) 4 8 1 h H2 (ζ ) = (−ζ 3 + 3ζ + 2) , H3 (ζ ) = (ζ 3 + ζ 2 − ζ −1) 4 8

(10)

where h is the laminate total thickness and ζ = 2z/h is the laminate nondimensional through the thickness coordinate whose domain is −1 ≤ ζ ≤ 1. Within layer k the axial and transverse displacements are therefore interpolated as: uk = uk−1L0 + uk L1 + ub H0 + φb H1 + ut H2 + φt H3

(11)

wk = wk−1L0 + wk L1 + wb H0 +ψbH1 + wt H2 +ψt H3

(12)

The rotational degrees of freedom φb and φt are now eliminated using the nonhomogeneous boundary conditions. Assuming that prescribed transverse shear stresses τxzb and τxzt at the bottom and top surfaces are known, Eqs. (6) and (9)-(12) can be used to yield

6

u  h1



 uN−1  + wt ,x + φt   hN 

τ xzb = Q551  1 + wb,x +φb  , τ xzt = Q55N  − 

(13)

Hence, φb and φt can be expressed as:

φb =

τ xzb 1 55

Q

− wb,x −

τ u u1 , φt = xztN − wt ,x + N−1 h1 Q55 hN

(14)

Similarly, the rotational degrees of freedom ψ b and ψ t can also be eliminated. Assuming that prescribed transverse normal stresses σzb and σzt at the bottom and top surfaces are known Eqs. (6) and (9)-(12) can again be used to yield:

 w1   w  +ψ b  , σ zt = Q13N ut ,x + Q33N  − N −1 +ψ t   h1   hN 

σ zb = Q131 ub, x + Q331 

(15)

Hence, ψ b and ψ t can be expressed as

ψb =

σ zb Q131 1 33

Q

σ zt

−

1 33

Q

ub,x −

w1 σ zb w = 1 −ν bub, x − 1 h1 Q33 h1

(16)

Q13N

w σ w ψ t = N − N ut,x + N −1 = ztN −ν t ut,x + N −1 hN Q33 hN Q33 Q33

where ν b = Q131 / Q331 and ν t = Q13N / Q33N . Equation (14) can be substituted in Eq. (11) resulting in τ τ u u  uk = uk−1L0 + uk L1 + ub H0 − wb,x H1 + ut H2 − wt ,x H3 +  xzb1 − 1 H1 +  xztN + N−1 H3  Q55 h1   Q55 hN 

(17)

Equation (16) can be substituted in Eq. (12) resulting in σ σ w w  wk = wk−1 L0 + wk L1 + wb H0 −ν b ub,x H1 + wt H 2 −ν t ut , x H3 +  zb1 − 1 H1 +  ztN + N −1 H3  Q33 h1   Q33 hN 

(18)

In order to compute strains the derivatives of u k and wk with respect to x and z are required. Those can be evaluated using Eqs. (17) and (18) to yield: u,kx = uk−1, x L0 + uk ,x L1 + ub,x H0 − wb,xxH1 + ut , x H2 − wt ,xxH3 − 7

u1,x h1

H1 +

uN −1,x hN

H3 +

τ xzb,x 1 55

Q

H1 +

τ xzt,x Q55N

H3

(19)

w,kx = wk −1,x L0 + wk ,x L1 + wb,x H0 −ν b ub, xx H1 + wt ,x H2 −ν t ut , xx H 3 −

σ zb,x Q331

H1 +

σ zt, x Q33N

w1,x h1

H1 +

wN −1,x hN

H3 +

(20)

H3

u,kz = uk−1L0,z + uk L1,z + ub H0,z − wb,x H1,z + ut H2,z − wt ,x H3,z −

τ τ u1 u H1,z + N−1 H3,z + xzb1 H1,z + xztN H3,z h1 hN Q55 Q55

(21)

w,kz = wk −1 L0, z + wk L1, z + wb H 0, z −ν b ub, x H1, z + wt H 2, z −ν t ut , x H 3, z −

(22)

w1 w σ σ H1, z + N −1 H 3, z + zb1 H1, z + ztN H3, z h1 hN Q33 Q33

The element formulation The degrees of freedom φb, φt, ψ b and ψ t have been explicitly eliminated from the formulation as indicated by Eqs. (14) and (16). Thus, N + 1 degrees of freedom related to the axial displacement (ub, u 1, ..., u N−1, ut) and N + 1 degrees of freedom related to the transverse displacement (wb, w1, ..., wN−1, wt) remain. However, most of these degrees of freedom can be computed as functions of only u b, u t, wb, and wt, as this section demonstrates. In order to achieve this reduction in the number of degrees of freedom the stress continuity along the thickness will be imposed. Stresses σzk and τxzk computed in Eqs. (6) must be continuous, provided there is no delamination. When written in terms of displacements σzk and τxzk are given by:

σ zk = Q13k ε xk + Q33k ε zk = Q13k u,kx + Q33k w,kz

(23)

τ xzk = Q55k γ xzk = Q55k (u,kz + w,kx )

(24)

Continuity of the normal transverse stress σz in the interface between layers k and k + 1 requires that

σzk(ζk = 1) = σzk+1(ζk+1 = −1). Therefore, using Eqs. (19), (22) and (23) one can write  w − wk −1  Q33k  k + wb H0k, z +ψ b H1k ,z + wt H2k ,z +ψ t H3k, z  +  hk  k Q13 (uk ,x + ub,x H0k +φb,x H1k + ut ,x H2k +φt , x H3k ) =  w −w  Q  k +1 k + wb H0k,z +ψ b H1k ,z + wt H2k ,z +ψ t H3k, z  +  hk +1  k +1 Q13 (uk ,x + ub,x H0k +φb, x H1k + ut , x H2k +φt , x H3k )

(25)

k +1 33

where H0k, H1k, H2k, H3k are the Hermite cubics evaluated at the interface between layers k and k + 1 where z = zk. Substitution of Eqs. (14) and (16) into (25) leads to

8

w  w − Qˆ 33k wk −1 + (Qˆ 33k + Qˆ 33k +1 )wk − Qˆ 33k +1wk +1 + (Q33k +1 − Q33k ) 1 H1k,z − N−1 H3k, z  = hN  h1  τ xzt, x τ xzb, x  σ  σ (Q13k +1 − Q13k ) 1 H1k + N H3k  + (Q33k +1 − Q33k ) zb1 H1k ,z + ztN H3k, z  + Q55 Q33  Q55   Q33  k +1 k (Q13 − Q13 )(ub,x H0k + ut ,x H2k − wb,xx H1k − wt ,xx H3k ) +

(26)

(Q33k +1 − Q33k )(wb H0k ,z + wt H2k ,z −ν bub,x H1k ,z −ν t ut ,x H3k ,z ) + u1,x uN −1, x   (Q13k +1 − Q13k )uk, x − H1k + H3k  h1 hN   where Qˆ33k = Q33k / hk . For a laminate with N layers Eq. (26) must be enforced N − 1 times for k = 1, 2, ..., N − 1. Recalling that w0 = wN = 0 these N − 1 equations can be cast in a matrix form to yield

ˆ w = w w + w w +u u +u u + w w + w w +σ σ +σ σ + Q u Q 33 b b t t b, x b, x t, x t,x b, xx b, xx t , xx t , xx zb zb zt zt 13 , x

(27)

w = {w1 w2 K wN−2 wN −1}T

(28)

where

(29)

u = {u1 u2 K uN−2 uN−1}

T

2 33

1 33

1 33

σzb = {(Q − Q )H11,z / Q

3 33

2 33

1 33

(Q − Q )H12,z / Q

N−1 33

K (Q − Q N 33

1 T 33

)H1( N−1),z / Q }

(30)

σzt = {(Q332 − Q331 )H31,z / Q33N (Q333 − Q332 )H32,z / Q33N K (Q33N − Q33N−1)H3( N−1),z / Q33N }T 2 33

1 33

3 33

2 33

N−1 33

2 33

1 33

3 33

2 33

N−1 33

wb = {(Q − Q )H01,z (Q − Q )H02,z K (Q wt = {(Q − Q )H21,z (Q − Q )H22,z K (Q

N−2 33

−Q

N−2 33

−Q

)H0( N−1),z}

(32)

N−1 33

T

)H2( N−1),z}

(33)

)H2( N−2),z (Q − Q N 33

(31)

N−1 33

)H0( N−2),z (Q − Q N 33

T

wb,xx = {(Q131 − Q132 )H11 (Q132 − Q133 )H12 K (Q13N−2 − Q13N−1)H1( N−2) (Q13N−1 − Q13N )H1( N−1)}T

(34)

wt,xx = {(Q131 − Q132 )H31 (Q132 − Q133 )H32 K (Q13N−2 − Q13N−1)H3( N−2) (Q13N−1 − Q13N )H3( N−1)}T

(35)

2 13

1 13

2 33

1 33

3 13

2 13

N −1 13

ub,x = {(Q − Q )H01 (Q − Q )H02 K (Q − Q 3 33

N 13

2 33

)H0( N−1)} − T

{(Q − Q )ν b H11,z (Q − Q )νb H12,z K (Q − Q33N−1)νb H1(N−1),z }T N 33

ut ,x = {(Q132 − Q131 )H21 (Q133 − Q132 )H22 K (Q13N − Q13N−1 )H2( N−1)}T − {(Q332 − Q331 )νt H31,z (Q333 − Q332 )νt H32,z K (Q33N − Q33N−1 )ν t H3( N−1),z }T

 ˆ 1 ˆ 2 (Q332 − Q331 ) 2 H11,z − Qˆ33 0 Q33 + Q33 + h 1  3 2  ˆ 2 (Q33 − Q33) 2 3 3 − Q + H12,z Qˆ33 + Qˆ33 − Qˆ33 33  h1 ˆ = Q 33 (Q334 − Q333 ) 3 3 4  H13,z − Qˆ33 Qˆ33 + Qˆ33  h1  M M M  N N−1 ( Q − Q ) 33 33  H1( N−1),z 0 0 h1 

9

L L L O L

 (Q332 − Q331 ) H31,z  hN  (Q333 − Q332 )  − H32,z  hN  4 3 (Q − Q )  − 33 33 H33,z  hN  M  N N−1 (Q − Q ) Qˆ33N−1 + Qˆ33N − 33 33 H3( N−1),z  hN 

(36) (37)

−

(38)

 2 H11  1   0 0  (Q13 − Q13)1−  h1    (Q133 − Q132 ) − H12 Q133 − Q132 0  h 1  Q13 =  (Q134 − Q133 ) − H13 0 Q134 − Q133  h1  M M M  N N−1 ( Q − Q ) − 13 13 H 0 0 1( N−1)  h1 

L L L O L

 (Q132 − Q131 ) H31  hN   (Q133 − Q132 ) H32  hN  (Q134 − Q133 )  H33  hN  M  H3( N−1)  N N−1   (Q13 − Q13 )1+ hN  

(39)

Similarly, continuity of the transverse shear stress τxz in the interface between layers k and k + 1 requires that τxzk(ζk = 1) = τxzk+1(ζk+1 = −1). Therefore, using Eqs. (20), (21) and (24) one can write:  u −u  Q55k  k k −1 + ubH0k ,z +φb H1k ,z + ut H2k,z +φt H3k ,z + wk ,x + wb,x H0k +ψb,x H1k + wt,x H2k +ψt,x H3k  =  hk   u −u  Q55k+1 k+1 k + ub H0k,z + φbH1k ,z + ut H2k,z + φt H3k ,z + wk,x + wb,x H0k +ψb,x H1k + wt,x H2k +ψt,x H3k   hk+1 

(40)

Substitution of Eqs. (14) and (16) into (40) leads to u  u k k k +1 k +1 − Qˆ 55 uk −1 + (Qˆ 55 + Qˆ 55 )uk − Qˆ 55 uk +1 + (Q55k +1 − Q55k ) 1 H1k, z − N −1 H3k, z  = hN  h1  σ zt,x τ   σ zb,x  τ (Q55k +1 −Q55k ) xzb1 H1k, z + xztN H3k, z  + (Q55k +1 − Q55k ) 1 H1k + N H3k  + Q55 Q33  Q55   Q33  k +1 k (Q55 −Q55 )(ub H0k ,z + ut H2k, z − wb,x H1k ,z − wt,x H3k ,z ) +

(41)

(Q55k +1 −Q55k )(wb,x H0k + wt,x H2k −ν bub, xx H1k −ν t ut, xx H3k ) + w1, x wN−1, x   (Q55k +1 −Q55k ) wk ,x − H1k + H3k  h1 hN   k = Q55k / hk . For a laminate with N layers Eq. (41) must be enforced N − 1 times for k = 1, 2, where Qˆ55

..., N − 1. Recalling that u0 = u N = 0 these N − 1 equations can be cast in a matrix form to yield ˆ u = u u + u u + w w + w w + u u + u u +τ τ +τ τ + Q w Q 55 b b t t b,x b,x t ,x t ,x b,xx b,xx t ,xx t ,xx xzb xzb xzt xzt 55 ,x

(42)

τxzb = {(Q552 − Q551 )H11,z / Q551 (Q553 − Q552 )H12,z / Q551 K (Q55N − Q55N−1)H1( N−1),z / Q551 }T

(43)

τxzt = {(Q552 − Q551 )H31,z / Q55N (Q553 − Q552 )H32,z / Q55N K (Q55N − Q55N−1)H3( N−1),z / Q55N }T

(44)

where

10

ub = {(Q552 − Q551 )H01,z (Q553 − Q552 )H02,z K (Q55N−1 − Q55N−2 )H0( N−2),z (Q55N − Q33N−1)H0( N−1),z}T 2 55

1 55

3 55

2 55

N−1 55

ut = {(Q − Q )H21,z (Q − Q )H22,z K (Q 1 55

2 55

2 55

3 55

N−2 55

−Q

N−1 55

ub,xx = {(Q − Q )νb H11 (Q − Q )νbH12 K (Q

N−1 55

)H2(N−2),z (Q − Q N 55

(45)

)H2(N−1),z }

(46)

T

T

− Q )νb H1( N−1)}

(47)

ut ,xx = {(Q551 − Q552 )νt H31 (Q552 − Q553 )νt H32 K (Q55N−1 − Q55N )νt H3(N−1)}T

(48)

N 55

wb,x = {(Q552 − Q551 )(H01 − H11,z ) (Q553 − Q552 )(H02 − H12,z ) K (Q55N − Q55N−1)(H0(N−1) − H1( N−1),z )}T 2 55

1 55

3 55

2 55

N−1 55

wt ,x = {(Q − Q )(H02 − H13,z ) (Q − Q )(H22 − H32,z ) K (Q − Q N 55

 ˆ 1 ˆ 2 (Q − Q ) 2 H11,z − Qˆ55 0 Q55 + Q55 + h 1  3 2  ˆ 2 (Q55 − Q55) 2 3 3 − Q + H12,z Qˆ55 + Qˆ55 − Qˆ55 55  h1 ˆ = Q 55 (Q554 − Q553 ) 3 3 4  H13,z − Qˆ55 Qˆ55 + Qˆ55  h1  M M M  N N−1 ( Q − Q ) 55 55  H1( N−1),z 0 0 h1  2 55

1 55

 2 H11  1   0 0  (Q55 − Q55)1−  h1    (Q553 − Q552 ) − H12 Q553 − Q552 0  h 1  Q55 =  (Q554 − Q553 ) − H13 0 Q554 − Q553  h1  M M M  N N−1 ( Q − Q ) − 55 55 H 0 0 1( N−1)  h1 

L L L O L

L L L O L

)(H2( N−1) − H3( N−1),z )}

T

 (Q − Q ) H31,z  hN  (Q553 − Q552 )  − H32,z  hN  (Q554 − Q553 )  − H33,z  hN  M  N N−1 (Q − Q ) Qˆ55N−1 + Qˆ55N − 55 55 H3( N−1),z  hN  −

2 55

(49) (50)

1 55

 (Q552 − Q551 ) H31  hN   (Q553 − Q552 ) H32  hN  (Q554 − Q553 )  H33  hN  M  H3( N−1)  N N−1   (Q55 − Q55 )1+ hN  

(51)

(52)

Notice that material discontinuities between layers will reflect directly in σzT, σzb, σzt, τxzb, τxzt, ub, ˆ , Q , Q ˆ and Q . Matrices Q ˆ and Q ˆ ut, wb, wt, ub,x, ut,x, wb,x, wt,x, ub,xx, ut,xx, wb,xx, wt,xx, Q 33 13 55 55 33 55 have special forms: they are tridiagonal "N"-shaped, meaning that they are tridiagonal all over, except that their first and last columns are full. Therefore, efficient procedures can be easily devised for their inversions. Equations (27) and (42) may be approximately solved to obtain u and w as functions of wb, wt, ub, ut and their derivatives with respect to x. Once u and w are available it becomes a simple task to compute strains. Since a finite element shall be proposed, the axial displacement u will be assumed to vary at most parabolically with x such that u,xxx = 0 and ub,xxx = ut,xxx = 0, whereas w will be assumed to vary at most cubically with x such that w,xxxx = 0 and wb,xxxx = wt,xxxx = 0. Also, stresses

σzb, σzt, τxzb, τxzt are assumed to be constant in the element length; therefore, the derivatives with respect to x of the stresses at boundary conditions σzb, σzt, τxzb, τxzt are also assumed to be zero. Hence, Eq. (27) can be differentiated twice with respect to x leading to:

11

ˆ w = w w + w w +u u +u u + w w Q 33 , xx b b, xx t t , xx b, x b, xxx t , x t , xxx b, xx b, xxxx + wt , xx wt , xxxx + σ zbσ zb,xx + σ ztσ zt, xx + Q13u,xxx = wb wb, xx + wt wt ,xx

(53)

Equation (42) is differentiated once with respect to x and Eq. (53) is used leading to: ˆ −1 (u u + u u + w w + w w + Q w ) = u,x = Q 55 b b,x t t ,x b, x b, xx t , x t ,xx 55 , xx ˆ −1[u u + u u + w w + w w + Q Q ˆ −1 Q 55 b b,x t t ,x b, x b, xx t , x t ,xx 55 33(wb wb,xx + wt wt , xx )]

(54)

Equation (54) is now substituted in Eq. (27) to obtain w: ˆ −1 (w w + w w + u u + u u + w w + w w + σ σ + σ σ ) + w =Q 33 b b t t b, x b, x t , x t, x b, xx b, xx t , xx t , xx zb zb zt zt ˆ −1Q Q ˆ −1 ˆ −1 Q 33 13 55[ubub, x + ut ut , x + wb, x wb, xx + wt , x wt , xx + Q55Q33 (wb wb, xx + wt wt , xx )]

(55)

Since w is available in Eq. (55) it is possible to compute w,x as: ˆ −1(w w + w w + u u + u u + w w + w w ) + w,x = Q 33 b b,x t t ,x b,x b,xx t , x t , xx b,xx b,xxx t , xx t , xxx ˆ −1Q Q ˆ −1 ˆ −1 Q 33 13 55[ubub, xx + ut ut ,xx + wb,x wb, xxx + wt ,x wt ,xxx + Q55Q33(wb wb,xxx + wt wt , xxx)]

(56)

Substitution of Eq. (56) into Eq. (42) gives

ˆ −1[u u + u u + w w + w w + u u + u u + τ τ + τ τ + u=Q 55 b b t t b,x b,x t ,x t ,x b, xx b, xx t ,xx t ,xx xzb zb xzt zt ˆ −1(w w + w w + u u + u u + w w + w w ) + Q55Q 33 b b,x t t ,x b, x b, xx t ,x t ,xx b,xx b,xxx t , xx t , xxx

(57)

ˆ −1Q Q ˆ −1 ˆ −1 Q55Q 33 13 55[ubub,xx + ut ut ,xx + wb,x wb,xxx + wt, x wt , xxx + Q55Q33(wbwb, xxx + wt wt, xxx)] Examination of Eqs. (54)-(57) shows that the computation of u, w, u,x and w,x always involve multiplication of a vector by either wb, wt, ub, ut or their derivatives with respect to x. Tables 1 and 2 schematically present these multiplications. The following convention is now proposed: the vector located in row i and column j of either Tab. 1 or 2 is denoted sj,i. Hence, u, w, u,x and w,x can be written as: u = su,1 + su,ubub + su,utut + su,ubxxub,xx + su,utxxut ,xx + su,wbxwb,x + su,wtxwt ,x + su,wbxxxwb,xxx + su,wtxxxwt ,xxx

(58)

w = sw,1 + sw,ubxub,x + sw,utxut ,x + sw,wbwb + sw,wt wt + sw,wbxxwb,xx + sw,wtxxwt ,xx

(59)

u,x = sux,ubxub,x + sux,utxut,x + sux,wbxxwb,xx + sux,wtxxwt,xx

(60) 12

w,x = swx,ubxxub,xx + swx,utxxut ,xx + swx,wbxwb,x + swx,wtxwt ,x + swx,wbxxxwb,xxx + swx,wtxxxwt ,xxx

(61)

Equations (58)-(61) provide a quick and easy way of computing any of the values of u1, u2, ..., uN−1, w1, w2, ..., wN−1 or their first derivatives with respect to x. For instance, if uk,x (k = 1, 2, ..., N − 1) is required, then simply access Eq. (60) and retrieve the corresponding entry: uk,x = sux,ubx,k ub,x + sux,utx,k ut,x + sux,wbxx,k wb,xx + sux,wtxx,k wt,xx, where sux,ubx,k, sux,utx,k, sux,wbxx,k, sux,wtxx,k are the corresponding entries in vectors sux,ubx, sux,utx, sux,wbxx, sux,wtxx respectively. A 6-noded element is proposed as illustrated in Fig. 4. The parabolic interpolation functions for ub and ut in x are adopted to improve the element accuracy under transverse shear, but it calls for two nodes positioned in the middle of the bottom and top surfaces as shown in Fig. 4. In addition, in order to improve the element accuracy under bending, cubic interpolations for wb and wt are used. Therefore, the interpolation scheme is given by: ub (ξ ) = ub1N0 (ξ ) + ub2 N1(ξ ) + ub3N2 (ξ ) ut (ξ ) = ut1N01(ξ ) + ut 2 N1 (ξ ) + uN2 (ξ )

(62)

wb (ξ ) = wb1P0 (ξ ) +θ yb1P1(ξ ) + wb2 P2 (ξ ) +θ yb2 P3 (ξ ) wt (ξ ) = wt1P0 (ξ ) +θ yt1P1(ξ ) + wt 2 P2 (ξ ) +θyt2 P3 (ξ )

The interpolation functions present in Eq. (62) are

N0 (ξ ) =

ξ (ξ −1) 2

, N1(ξ) =

ξ (ξ +1) 2

, N2 (ξ ) = 1− ξ 2

l 1 P0 (ξ) = (ξ 3 − 3ξ + 2) , P1(ξ ) = e (ξ 3 − ξ 2 − ξ +1) 4 8 le 3 2 1 3 P2 (ξ ) = (−ξ + 3ξ + 2) , P3 (ξ ) = (ξ + ξ −ξ −1) 4 8

(63)

where le is the element length and ξ ∈ [−1, 1] is the nondimensional axial coordinate. The element vector of degrees of freedom is:

qe ={ub1 wb1 θyb1 ut1 wt1 θyt1 ub2 wb2 θyb2 ut 2 wt 2 θyt2 ub3 ut3}T

(64)

The strains within layer k (εxk, εzk, γxzk) can be computed with the aid of Eqs. (19)-(22). It is clear that the strains depend on u1, u2, ..., uN−1, w1, w2, ..., wN−1 and their derivatives with respect to x.

13

However, these quantities are computed from Eqs. (58)-(61). Therefore, substitution of Eqs. (58)(61) in Eqs. (19)-(22) results in:

ε xk = (L0sux,ubx,k−1 + L1sux,ubx,k + H0 + sux,ubx,N −1H3 / hN − sux,ubx,1H1 / h1)ub,x + (L0sux,utx,k−1 + L1sux,utx,k + H2 + sux,utx,N −1H3 / hN − sux,utx,1H1 / h1)ut,x + (L0sux,wbxx,k −1 + L1sux,wbxx,k − H1 + sux,wbxx,N −1H3 / hN − sux,wbxx,1H1 / h1)wb,xx +

(65)

(L0sux,wtxx,k −1 + L1sux,wtxx,k − H3 + sux,wtxx, N−1H3 / hN − sux,wtxx,1H1 / h1)wt ,xx

ε zk = (L0,z sw,1,k−1 + L1,z sw,1,k + sw,1,N−1H3,z / hN − sw,1,1H1,z / h1) + (L0,z sw,ubx,k−1 + L1,z sw,ubx,k + sw,ubx,N−1H3,z / hN − sw,ubx,1H1,z / h1 −νbH1,z )ub,x + (L0,z sw,utx,k−1 + L1,z sw,utx,k + sw,utx,N−1H3,z / hN − sw,utx,1H1,z / h1 −ν t H3,z )ut,x + (L0,z sw,wb,k−1 + L1,z sw,wb,k + sw,wb,N−1H3,z / hN − sw,wb,1H1,z / h1 + H0,z )wb + (L0,z sw,wt,k−1 + L1,z sw,wt,k + sw,wt,N−1H3,z / hN − sw,wt,1H1,z / h1 + H2,z )wt +

(66)

(L0,z sw,wbxx,k−1 + L1,z sw,wbxx,k + sw,wbxx,N−1H3,z / hN − sw,wbxx,1H1,z / h1)wb,xx + (L0,z sw,wtxx,k−1 + L1,z sw,wtxx,k + sw,wtxx,N−1H3,z / hN − sw,wtxx,1H1,z / h1)wt,xx + (σ zb + σ1zT )H1,z / Q331 + (σ zt +σ zTN )H3,z / Q33N

γ xzk = (L0,z su,1,k−1 + L1,z su,1,k + su,1,N−1H3,z / hN − su,1,1H1,z / h1) (L0,z su,ub,k−1 + L1,z su,ub,k + su,ub,N−1H3,z / hN − su,ub,1H1,z / h1 + H0,z )ub + (L0,z su,ut,k−1 + L1,z su,ut,k + su,ut,N−1H3,z / hN − su,ut,1H1,z / h1 + H2,z )ut + (L0,z su,ubxx,k−1 + L1,z su,ubxx,k + su,ubxx,N−1H3,z / hN − su,ubxx,1H1,z / h1)ub,xx + (L0,z su,utxx,k−1 + L1,z su,utxx,k + su,utxx,N−1H3,z / hN − su,utxx,1H1,z / h1)ut,xx + (L0,z su,wbx,k−1 + L1,z su,wbx,k + su,wbx,N−1H3,z / hN − su,wbx,1H1,z / h1 − H1,z )wb,x + (L0,z su,wtx,k−1 + L1,z su,wtx,k + su,wtx,N−1H3,z / hN − su,wtx,1H1,z / h1 − H3,z )wt,x + (L0,z su,wbxxx,k−1 + L1,z su,wbxxx,k + su,wbxxx,N−1H3,z / hN − su,wbxxx,1H1,z / h1)wb,xxx + (L0,z su,wtxxx,k−1 + L1,z su,wtxxx,k + su,wtxxx,N−1H3,z / hN − su,wtxxx,1H1,z / h1)wt,xxx + (L0swx,ubxx,k−1 + L1swx,ubxx,k + swx,ubxx,N−1H3 / hN − swx,ubxx,1H1 / h1 −νb H1)ub,xx + (L0swx,utxx,k−1 + L1swx,utxx,k + swx,utxx,N−1H3 / hN − swx,utxx,1H1 / h1 −νt H3 )ut,xx + (L0swx,wbx,k−1 + L1swx,wbx,k + swx,wbx,N−1H3 / hN − swx,wbx,1H1 / h1 + H0 )wb,x + (L0swx,wtx,k−1 + L1swx,wtx,k + swx,wtx,N−1H3 / hN − swx,wtx,1H1 / h1 + H2 )wt,x + (L0swx,wbxxx,k−1 + L1swx,wbxxx,k + swx,wbxxx,N−1H3 / hN − swx,wbxxx,1H1 / h1)wb,xxx + (L0swx,wtxxx,k−1 + L1swx,wtxxx,k + swx,wtxxx,N−1H3 / hN − swx,wtxxx,1H1 / h1)wt,xxx + 1 55

(67)

τ xzbH1,z / Q +τ xztH3,z / Q

N 55

Equations (65)-(67) can be cast in a more compact form if one defines matrices Bk and BSk , whose terms are given in Eqs. (68)-(70).

14

b1k,1 = (L0sux,ubx,k−1 + L1sux,ubx,k + H0 + sux,ubx,N−1H3 / hN − sux,ubx,1H1 / h1)N0,x b1k,2 = (L0sux,wbxx,k−1 + L1sux,wbxx,k − H1 + sux,wbxx,N−1H3 / hN − sux,wbxx,1H1 / h1)P0,xx b1k,3 = (L0sux,wbxx,k−1 + L1sux,wbxx,k − H1 + sux,wbxx,N−1H3 / hN − sux,wbxx,1H1 / h1)P1,xx b1k,4 = (L0sux,utx,k−1 + L1sux,utx,k + H2 + sux,utx,N−1H3 / hN − sux,utx,1H1 / h1)N0,x b1k,5 = (L0sux,wtxx,k−1 + L1sux,wtxx,k − H3 + sux,wtxx,N−1H3 / hN − sux,wtxx,1H1 / h1)P0,xx b1k,6 = (L0sux,wtxx,k−1 + L1sux,wtxx,k − H3 + sux,wtxx,N−1H3 / hN − sux,wtxx,1H1 / h1)P1,xx b1k,7 = (L0sux,ubx,k−1 + L1sux,ubx,k + H0 + sux,ubx,N−1H3 / hN − sux,ubx,1H1 / h1)N1,x b1k,8 = (L0sux,wbxx,k−1 + L1sux,wbxx,k − H1 + sux,wbxx,N−1H3 / hN − sux,wbxx,1H1 / h1)P2,xx b1k,9 = (L0sux,wbxx,k−1 + L1sux,wbxx,k − H1 + sux,wbxx,N−1H3 / hN − sux,wbxx,1H1 / h1)P3,xx b1k,10 = (L0sux,utx,k−1 + L1sux,utx,k + H2 + sux,utx,N−1H3 / hN − sux,utx,1H1 / h1)N1,x b1k,11 = (L0sux,wtxx,k−1 + L1sux,wtxx,k − H3 + sux,wtxx,N−1H3 / hN − sux,wtxx,1H1 / h1)P2,xx b1k,12 = (L0sux,wtxx,k−1 + L1sux,wtxx,k − H3 + sux,wtxx,N−1H3 / hN − sux,wtxx,1H1 / h1)P3,xx b1k,13 = (L0sux,ubx,k−1 + L1sux,ubx,k + H0 + sux,ubx,N−1H3 / hN − sux,ubx,1H1 / h1)N2,x b1k,14 = (L0sux,utx,k−1 + L1sux,utx,k + H2 + sux,utx,N−1H3 / hN − sux,utx,1H1 / h1)N2,x

15

(68)

b2k,1 = ( L0, z sw ,ubx, k −1 + L1, z sw ,ubx, k + sw ,ubx, N −1H 3, z / hN − sw ,ubx,1H1, z / h1 − ν b H1, z ) N 0, x b2k, 2 = ( L0, z sw , wb , k −1 + L1, z sw , wb , k + sw , wb , N −1H 3, z / hN − sw , wb ,1H1, z / h1 + H 0, z ) P0 + ( L0, z sw , wbxx , k −1 + L1, z sw , wbxx , k + sw , wbxx , N −1H 3, z / hN − sw , wbxx ,1H1, z / h1 ) P0, xx b2k,3 = ( L0, z sw , wb , k −1 + L1, z sw , wb , k + sw , wb , N −1H 3, z / hN − sw , wb ,1H1, z / h1 + H 0, z ) P1 + ( L0, z sw , wbxx , k −1 + L1, z sw , wbxx , k + sw , wbxx , N −1H 3, z / hN − sw , wbxx ,1H1, z / h1 ) P1, xx k 2, 4

= ( L0, z sw ,utx, k −1 + L1, z sw ,utx, k + sw ,utx , N −1H 3, z / hN − sw ,utx,1H1, z / h1 − ν t H 3, z ) N 0, x

k 2, 5

= ( L0, z sw , wt , k −1 + L1, z sw , wt , k + sw , wt , N −1H 3, z / hN − sw , wt ,1H1, z / h1 + H 2, z ) P0 +

b b

( L0, z sw , wtxx , k −1 + L1, z sw , wtxx , k + sw , wtxx , N −1H 3, z / hN − sw , wtxx ,1H1, z / h1 ) P0, xx b2k, 6 = ( L0, z sw , wt , k −1 + L1, z sw , wt , k + sw, wt , N −1H 3, z / hN − sw , wt ,1H1, z / h1 + H 2, z ) P1 + ( L0, z sw , wtxx , k −1 + L1, z sw , wtxx , k + sw , wtxx , N −1H 3, z / hN − sw , wtxx ,1H1, z / h1 ) P1, xx k 2, 7

= ( L0, z sw ,ubx, k −1 + L1, z sw ,ubx, k + sw ,ubx, N −1H 3, z / hN − sw ,ubx,1H1, z / h1 − ν b H1, z ) N1, x

k 2, 8

= ( L0, z sw , wb , k −1 + L1, z sw , wb , k + sw , wb , N −1H 3, z / hN − sw , wb ,1H1, z / h1 + H 0, z ) P2 +

b b

( L0, z sw , wbxx , k −1 + L1, z sw , wbxx , k + sw , wbxx , N −1H 3, z / hN − sw , wbxx ,1H1, z / h1 ) P2, xx b2k,9 = ( L0, z sw , wb , k −1 + L1, z sw , wb , k + sw , wb , N −1H 3, z / hN − sw , wb ,1H1, z / h1 + H 0, z ) P3 + ( L0, z sw , wbxx , k −1 + L1, z sw , wbxx , k + sw , wbxx , N −1H 3, z / hN − sw , wbxx ,1H1, z / h1 ) P3, xx k 2,10

b

= ( L0, z sw ,utx ,k −1 + L1, z sw ,utx, k + sw ,utx, N −1H 3, z / hN − sw ,utx,1H1, z / h1 − ν t H 3, z ) N1, x

b2k,11 = ( L0, z sw , wt , k −1 + L1, z sw , wt , k + sw , wt , N −1H 3, z / hN − sw , wt ,1H1, z / h1 + H 2, z ) P2 + ( L0, z sw , wtxx , k −1 + L1, z sw , wtxx , k + sw , wtxx , N −1H 3, z / hN − sw , wtxx ,1H1, z / h1 ) P2, xx b2k,12 = ( L0, z sw , wt , k −1 + L1, z sw , wt , k + sw, wt , N −1H 3, z / hN − sw , wt ,1H1, z / h1 + H 2, z ) P3 + ( L0, z sw , wtxx , k −1 + L1, z sw , wtxx , k + sw , wtxx , N −1H 3, z / hN − sw , wtxx ,1H1, z / h1 ) P3, xx b2k,13 = ( L0, z sw ,ubx, k −1 + L1, z sw ,ubx, k + sw ,ubx, N −1H 3, z / hN − sw ,ubx,1H1, z / h1 −ν b H1, z ) N 2, x b2k,14 = ( L0, z sw ,utx ,k −1 + L1, z sw ,utx, k + sw ,utx, N −1H 3, z / hN − sw ,utx,1H1, z / h1 − ν t H 3, z ) N 2, x + bs1k,1 = (L0,z su,ub,k−1 + L1,z su,ub,k + su,ub,N−1H3,z / hN − su,ub,1H1,z / h1 + H0,z ) N0 (L0,z su,ubxx,k−1 + L1,z su,ubxx,k + su,ubxx,N−1H3,z / hN − su,ubxx,1H1,z / h1 )N0,xx (L0 swx,ubxx,k−1 + L1swx,ubxx,k + swx,ubxx,N−1H3 / hN − swx,ubxx,1H1 / h1 −νb H1) N0,xx bs = (L0,z su,wbx,k−1 + L1,z su,wbx,k + su,wbx,N−1H3,z / hN − su,wbx,1H1,z / h1 − H1,z )P0,x + k 1, 2

(L0 swx,wbx,k−1 + L1swx,wbx,k + swx,wbx,N−1H3 / hN − swx,wbx,1H1 / h1 + H0 )P0,x + (L0,z su,wbxxx,k−1 + L1,z su,wbxxx,k + su,wbxxx,N−1H3,z / hN − su,wbxxx,1H1,z / h1 )P0,xxx + (L0 swx,wbxxx,k−1 + L1swx,wbxxx,k + swx,wbxxx,N−1H3 / hN − swx,wbxxx,1H1 / h1)P0,xxx bs1k,3 = (L0,z su,wbx,k−1 + L1,z su,wbx,k + su,wbx,N−1H3,z / hN − su,wbx,1H1,z / h1 − H1,z )P1,x + (L0 swx,wbx,k−1 + L1swx,wbx,k + swx,wbx,N−1H3 / hN − swx,wbx,1H1 / h1 + H0 )P1,x + (L0,z su,wbxxx,k−1 + L1,z su,wbxxx,k + su,wbxxx,N−1H3,z / hN − su,wbxxx,1H1,z / h1 )P1,xxx + (L0 swx,wbxxx,k−1 + L1swx,wbxxx,k + swx,wbxxx,N−1H3 / hN − swx,wbxxx,1H1 / h1)P1,xxx

16

(69)

bs1k,4 = (L0,z su,ut,k−1 + L1,z su,ut,k + su,ut,N−1H3,z / hN − su,ut,1H1,z / h1 + H2,z )N0 (L0,z su,utxx,k−1 + L1,z su,utxx,k + su,utxx,N−1H3,z / hN − su,utxx,1H1,z / h1)N0,xx (L0swx,utxx,k−1 + L1swx,utxx,k + swx,utxx,N−1H3 / hN − swx,utxx,1H1 / h1 −νt H3 )N0,xx bs1k,5 = (L0,z su,wtx,k−1 + L1,z su,wtx,k + su,wtx,N−1H3,z / hN − su,wtx,1H1,z / h1 − H3,z )P0,x + (L0swx,wtx,k−1 + L1swx,wtx,k + swx,wtx,N−1H3 / hN − swx,wtx,1H1 / h1 + H2 )P0,x + (L0,z su,wtxxx,k−1 + L1,z su,wtxxx,k + su,wtxxx,N−1H3,z / hN − su,wtxxx,1H1,z / h1)P0,xxx + (L0swx,wtxxx,k−1 + L1swx,wtxxx,k + swx,wtxxx,N−1H3 / hN − swx,wtxxx,1H1 / h1)P0,xxx bs = (L0,z su,wtx,k−1 + L1,z su,wtx,k + su,wtx,N−1H3,z / hN − su,wtx,1H1,z / h1 − H3,z )P1,x + k 1,6

(L0swx,wtx,k−1 + L1swx,wtx,k + swx,wtx,N−1H3 / hN − swx,wtx,1H1 / h1 + H2 )P1,x + (L0,z su,wtxxx,k−1 + L1,z su,wtxxx,k + su,wtxxx,N−1H3,z / hN − su,wtxxx,1H1,z / h1)P1,xxx + (L0swx,wtxxx,k−1 + L1swx,wtxxx,k + swx,wtxxx,N−1H3 / hN − swx,wtxxx,1H1 / h1)P1,xxx bs1k,7 = (L0,z su,ub,k−1 + L1,z su,ub,k + su,ub,N−1H3,z / hN − su,ub,1H1,z / h1 + H0,z )N1 (L0,z su,ubxx,k−1 + L1,z su,ubxx,k + su,ubxx,N−1H3,z / hN − su,ubxx,1H1,z / h1)N1,xx (L0swx,ubxx,k−1 + L1swx,ubxx,k + swx,ubxx,N−1H3 / hN − swx,ubxx,1H1 / h1 −νb H1)N1,xx bs1k,8 = (L0,z su,wbx,k−1 + L1,z su,wbx,k + su,wbx,N−1H3,z / hN − su,wbx,1H1,z / h1 − H1,z )P2,x + (L0swx,wbx,k−1 + L1swx,wbx,k + swx,wbx,N−1H3 / hN − swx,wbx,1H1 / h1 + H0 )P2,x + (L0,z su,wbxxx,k−1 + L1,z su,wbxxx,k + su,wbxxx,N−1H3,z / hN − su,wbxxx,1H1,z / h1)P2,xxx + (L0swx,wbxxx,k−1 + L1swx,wbxxx,k + swx,wbxxx,N−1H3 / hN − swx,wbxxx,1H1 / h1)P2,xxx bs = (L0,z su,wbx,k−1 + L1,z su,wbx,k + su,wbx,N−1H3,z / hN − su,wbx,1H1,z / h1 − H1,z )P3,x + k 1,9

(L0swx,wbx,k−1 + L1swx,wbx,k + swx,wbx,N−1H3 / hN − swx,wbx,1H1 / h1 + H0 )P3,x + (L0,z su,wbxxx,k−1 + L1,z su,wbxxx,k + su,wbxxx,N−1H3,z / hN − su,wbxxx,1H1,z / h1)P3,xxx + (L0swx,wbxxx,k−1 + L1swx,wbxxx,k + swx,wbxxx,N−1H3 / hN − swx,wbxxx,1H1 / h1)P3,xxx bs1k,10 = (L0,z su,ut,k−1 + L1,z su,ut,k + su,ut,N−1H3,z / hN − su,ut,1H1,z / h1 + H2,z )N1 (L0,z su,utxx,k−1 + L1,z su,utxx,k + su,utxx,N−1H3,z / hN − su,utxx,1H1,z / h1)N1,xx (L0swx,utxx,k−1 + L1swx,utxx,k + swx,utxx,N−1H3 / hN − swx,utxx,1H1 / h1 −νt H3 )N1,xx bs1k,11 = (L0,z su,wtx,k−1 + L1,z su,wtx,k + su,wtx,N−1H3,z / hN − su,wtx,1H1,z / h1 − H3,z )P2,x + (L0swx,wtx,k−1 + L1swx,wtx,k + swx,wtx,N−1H3 / hN − swx,wtx,1H1 / h1 + H2 )P2,x + (L0,z su,wtxxx,k−1 + L1,z su,wtxxx,k + su,wtxxx,N−1H3,z / hN − su,wtxxx,1H1,z / h1)P2,xxx + (L0swx,wtxxx,k−1 + L1swx,wtxxx,k + swx,wtxxx,N−1H3 / hN − swx,wtxxx,1H1 / h1)P2,xxx bs = (L0,z su,wtx,k−1 + L1,z su,wtx,k + su,wtx,N−1H3,z / hN − su,wtx,1H1,z / h1 − H3,z )P3,x + k 1,12

(L0swx,wtx,k−1 + L1swx,wtx,k + swx,wtx,N−1H3 / hN − swx,wtx,1H1 / h1 + H2 )P3,x + (L0,z su,wtxxx,k−1 + L1,z su,wtxxx,k + su,wtxxx,N−1H3,z / hN − su,wtxxx,1H1,z / h1)P3,xxx + (L0swx,wtxxx,k−1 + L1swx,wtxxx,k + swx,wtxxx,N−1H3 / hN − swx,wtxxx,1H1 / h1)P3,xxx

17

bs1k,13 = (L0,z su,ub,k−1 + L1,z su,ub,k + su,ub,N−1H3,z / hN − su,ub,1H1,z / h1 + H0,z ) N2 (L0,z su,ubxx,k−1 + L1,z su,ubxx,k + su,ubxx,N−1H3,z / hN − su,ubxx,1H1,z / h1)N2,xx (L0swx,ubxx,k−1 + L1swx,ubxx,k + swx,ubxx,N−1H3 / hN − swx,ubxx,1H1 / h1 −ν bH1) N2,xx

(70)

bs1k,14 = (L0,z su,ut,k−1 + L1,z su,ut,k + su,ut,N−1H3,z / hN − su,ut,1H1,z / h1 + H2,z )N2 (L0,z su,utxx,k−1 + L1,z su,utxx,k + su,utxx,N−1H3,z / hN − su,utxx,1H1,z / h1)N2,xx (L0swx,utxx,k−1 + L1swx,utxx,k + swx,utxx,N−1H3 / hN − swx,utxx,1H1 / h1 −ν t H3 )N2,xx

The strains at layer k are therefore:

0 ε xk  k    k  = B qe +  + ε z  L0,z sw,1,k −1 + L1,z sw,1,k + sw,1,N −1H3,z / hN − sw,1,1H1,z / h1  0    1 N H1,z (σ zb ) / Q33 + H3, z (σ zt ) / Q33 

γ xzk = BkS qe + L0,z su,1,k−1 + L1,z su,1,k + su,1,N−1H3,z / hN − su,1,1H1,z / h1 + H1,z

τ xzb 1 55

Q

+ H3,z

(71)

τ xzt Q55N

(72)

The virtual work of the stresses within element e is:

 k T Q11k Q13k  k  (B )  B + (BkS )T Q55k BkS qe dxdz  ∫ k k   k =1 zk −1 x =0  Q13 Q33   z l k k N k e 0 Q Q   δqTe ∑ ∫ ∫ (Bk )T  11k 13k  dxdz + k =1 zk −1 x =0 Q13 Q33 L0,z sw,1,k −1 + L1,z sw,1,k + sw,1,N −1 H3,z / hN − sw,1,1 H1, z / h1  N

zk

δUe = δqTe ∑ ∫

N

zk

δqTe ∑ ∫

le

le

∫ (B )

k T S

Q55k (L0,z su,1,k −1 + L1,z su,1,k + su,1,N −1 H3,z / hN − su,1,1 H1,z / h1 )dxdz

(73)

k =1 zk −1 x =0

T e

δq

k 0  k   k T Q11 Q13  ( B ) ∑  k 1 N dxdz + ∫ ∫ k  k =1 zk −1 x =0 Q13 Q33 H1, z (σ zb ) / Q33 + H3,z (σ zt ) / Q33  N

zk

le

N

zk

le

δqTe ∑ ∫

∫ (B )

k =1 zk −1 x =0

k T S

 τ τ  Q55k  H1,z xzb1 + H3, z xztN dxdz Q55 Q55  

The first integral in Eq. (73) contains the element stiffness matrix and the other integrals are consequences of stresses prescribed on the top and bottom surfaces.

Numerical procedures Initially, notice that the element topology, as presented in Fig. 4, is unnecessary - i.e. the actual coordinates of nodes 1 to 6 must not be informed. All that is required is that the x coordinate of either node 1 or 2, node 3 or 4, and node 5 or 6 are provided, and that the thicknesses of the layers 18

composing the laminate are given. This remark is important because it allows for the use of standard 1D mesh generators that geometrically represent a beam element as a straight line segment. The numerical procedures involved are straightforward, except for special pre-processing steps that are detailed in this section. In the pre-processing phase the vectors presented in Tabs. 1 and 2, and subsequently used in Eqs. (58)-(61), are computed and stored. This must be done only once for each beam element in the ˆ and Q ˆ is required. However, explicit mesh. Notice that, apparently, the inversion of both Q 33 55 ˆ −1u is to be computed, then the system inversion is not necessary. For instance, if su,ub = Q 55 b ˆ s = u is solved. Thus, matrix Q ˆ can be LU decomposed in-place only once and then LUsu,ub Q 55 u,ub b 55 ˆ matrix. = ub is solved. The same applies to the Q 33 Once the required vectors of Tabs. 1 and 2 are computed a numerical integration procedure starts to evaluate the element stiffness matrix and load vectors presented in Eq. (73). The highest degree of polynomials defined in the x-direction is 3 and the highest degree of polynomials defined in the zdirection is also 3. Thus, in the evaluation of the element arrays in Eq. (73) polynomials of highest degree 6 will appear, calling for 4 integration points along x and z. Thus, the term related to the shear stress in Eq. (73), (BkS )T Q55k BkS , can be fully integrated. Implementations that are more efficient can be achieved if the integrations are carried out analytically. According to Eqs. (68)-(70), the element matrices and vectors involve 14 degrees of freedom. However, the two mid-side nodes with the two degrees of freedom u b3 and ut3 can be eliminated in the element level. The usual procedure for this elimination is presented in Eq. (74)

Kii Kid qi  fi  −1 T −1 KT K q  = f  ⇒ (Kii − Kid KddKid )qi = fi − Kid Kddf d dd  d   d   id

(74)

where qi = {ub1 wb1 θyb1 ut1 wt1 θyt1 ub2 wb2 θyb2 ut2 wt2 θyt2}T and qd = {ub3 ut3}T. Matrix Kdd in Eq. (74) is 2×2 and can be analytically inverted. Once matrix K = K ii − K id K −dd1 K Tid and vector

f = fi − Kid K−dd1 f d are computed, the usual assembly procedure continues to store element arrays in the global arrays. After the assembly and solution of the global systems of algebraic equations, Kqi = f, the nodal displacements vector qi is available. Recall that the mid-node degrees of freedom

19

were eliminated through Eq. (74). Those degrees of freedom are needed to evaluate strains and stresses; therefore, they must be recomputed using qd = (Kdd)−1(fd − KidTqi).

Numerical examples In this section, the new finite element capabilities are assessed. First the element solution was compared to the analytical elasticity solution by Pagano [18] for a problem of a composite laminate beam under cylindrical bending. The analyzed problem, illustrated in Fig. 5, consists of a threelayered composite of length L and thickness h under a uniform distributed load q. In his original paper, Pagano solves the problem for a sinusoidal load given by:

q = q0 sin(px) , p =

nπ L

(75)

However, the constant load distribution can be expressed in the form of a Fourier series as:

q=

4

1  nπ  sin  π n=1,3,5,... n  L  ∞

∑

(76)

It is important to use a large number of terms in the Fourier expansion in order to have a good approximation for q. Therefore, in the solution tested, n in Eq. (76) varies from 1 to 2001. Also, Pagano enforces σx(0, y) = σy(L, y) = 0, but in the finite element method (FEM) this boundary condition must not be enforced, since it is inherently imposed by the method. Therefore, based on St. Venant's principle, in order to avoid effects of the boundary conditions, the strains, stresses and displacements were evaluated at x = 3L/4 and compared to the new finite element solution. Two different material configurations were evaluated: Beam (1): a metallic sandwich with the outer layers in aluminum (al) with a steel (st) core. Est = 210 GPa, νst = 0.3, Gst = Est /2(1 + νst), Eal = 70 GPa, νal = 0.3, Gal = Eal/2(1 + νal). Beam (2): 0°/90°/0° composite laminate [13]. E1 = 172.5 GPa, E2 = E3 = 6.9 GPa, G12 = G13 = 3.45 GPa, G23 = 1.38 GPa, ν12 = ν13 = ν23 = 0.25. For a unidirectional lamina oriented with 0°, directions (1, 2, 3) are aligned with (x, y, z) directions. The beams were analyzed with different S = L/h ratios to assess the behavior of thin (S = 20) and thick (S = 4) beams. Figures 6 and 7 present the results of the normalized strains ((a) to (f)), stresses ((g) to (l)) and displacements ((m) and (n)) along the normalized thickness. All quantities were 20

normalized by their maximum absolute value, computed based on the elasticity solution (ELST). As for the normalized thickness, the values for z are −0.5 ≤ z ≤ 0.5.

No special treatments were performed to recover any of the quantities - i.e. once nodal displacements are obtained, strains are computed directly from Eqs. (71)-(72), whereas stresses are computed from Eqs. (23)-(24). Good agreement for εx, σx and u are observed, for both the metallic sandwich and composite laminate. It is noted in Figs. 6e, 6f, 6k, 6l, 7e, 7f, 7k and 7l that τxz and γxz are affected by the length-tothickness ratio S, deviating from ELST solution as S decreases (thicker beams). Nevertheless, the transverse shear stresses predicted by the new finite element are above those predicted by the elasticity solution for this problem. Furthermore, S = 20 is already a relatively thick composite laminate, revealing an accurate approximation for usual composite structures, with no need of integration of the shear strain and stress. If thicker structures are to be modeled, the integration of the equilibrium equation in Eq. (77) can be performed to recover more precise values, since σx results are very accurate.

z

τ xz = −∫σx,xdz

(77)

z0

For σz, it is seen in Figs. 6i, 6j, 7i and 7j that the present FEM solution significantly deviates from the elasticity solution for both thin and thick beams. Integration of the equilibrium equation in the zdirection in Eq. (78) produces much better results, represented by the curves FEM-I.

z

σz = −∫τ xz,xdz

(78)

z0

Therefore, although the computation of a shear correction by integration is not essential for many common laminate configurations (higher S ratios), transverse normal stress integration through the thickness is recommended as prescribed by Eq. (78). Concerning εz, Figs. 7c and 7d reveal that, for usual composites, mechanical properties vary significantly between layers (the Young modulus in the longitudinal direction E1 is 25 times higher in the outer layers than in the central layer), the present finite element solution significantly deviates

21

from the elasticity solution. A more accurate approximation to εz can be evaluated by Eq. (79), using the very accurate εx and the integrated σz obtained from Eq. (78).

ε zk =

σ zk − Q13kε xk

(79)

Q33k

The behavior of the FEM solution with respect to the S ratio is presented in Fig. 8a for the metallic sandwich and Fig. 8b for the 0°/90°/0° laminate. In those figures, the maximum normalized deflection is given by

wmax =

100E*h3w(L / 2,0) q0L4

(80)

where E* is Est for the metallic sandwich and E2 for the composite laminate. It is important to state that the form of the Pagano solution for isotropic and transversely isotropic materials differs from the one for orthotropic materials (lamina). The results in Fig. 8 reinforce that for thinner beams the finite element solution converges to the analytical solution. And most importantly, the present finite element formulation follows the same trend as the analytical solution, a behavior not shared by the classical laminate theory [18]. The last discussed example is a cantilever beam with a tearing load in its free tip, Fig. 9. This simple example does not possess an analytical solution; therefore, the results where compared to a finite element model in MSC Nastran using 2D QUAD8 elements. A mesh with 1800 QUAD8 elements was used for the Nastran model (3 along the thickness of each layer and 200 along the beam length), and 32 equal length beam elements were used with the present formulation with the 0°/90°/0° laminate mechanical properties. In Figs. 10 and 11 it is seen that the w and σz fields predicted along the thickness by the presently formulated element are in accordance with the QUAD8 elements used in Nastran. This reveals a great computational advantage for the present beam element formulation.

Conclusions The present work presents a new finite element beam formulation that satisfies Koiter’s recommendations, considering effects of transverse shear and normal stresses, and C0z requirements 22

by a global-local superposition displacement field. The element is developed with nonhomogeneous boundary conditions and the number of degrees of freedom is kept independent of the number of layers, with all of them having a clear physical significance. The element comes to meet the needs highlighted by recent reviews in the scientific literature concerning beam and plate theories for composite structures. The present element delivers results that follow the same tendency of those obtained from the elasticity solution. The element behavior is accurate for usual length-to-thickness ratios in composite structures, dismissing the use of integration of the equilibrium equation to recover transverse shear stress and strain. However, the transverse normal stress and strain may significantly differ from those analytically obtained, even for thin structures, suggesting that a more safe procedure is to evaluate them by post-process integration of the equilibrium equation. The element also proved to be an efficient alternative to commercial finite element solution and other existent multilayered theories to model a problem with a tearing load, a loading condition with no analytical solution available. Straightforward expansion to laminated plate elements is envisioned. In this case, perhaps the only drawback would be the necessity to use 8-noded serendipity-like elements, because the elimination of the in-plane displacements in the element level, as done in Eq. (74), is no longer possible. The alternative would be to maintain the mid-side nodes in the mesh, what is not very practical since bookkeeping would be cumbersome (different nodes with different number of degrees of freedom). However, a 4-noded bilinear element combined with the numerical integrations put forward in Eqs. (77) and (78) might also deliver good results. The newly formulated element is applicable to situations when top (τxzt σzt) and bottom (τxzb σzb) surface stresses play a decisive role. For instance, normal contact forces/stresses in contact problems modeled by beam or plate elements could be precisely taken into account. Delamination might also be predicted, as long as the element formulation is enhanced to allow for the possibility of crack openings in between layers. Thermal effects can also be easily included in the formulation and will be the core of future developments, to address the need highlighted by the scientific community. With such capabilities, the element will be applicable to study very important phenomena, as fabrication induced distortions in composite laminates, fibre-metal laminates and sandwich beams.

23

REFERENCES [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18]

Timoshenko SP. X. On the transverse vibrations of bars of uniform cross-section. Philos Mag 1922;43:125–31. doi:10.1080/14786442208633855. Kirchhoff G. Ueber das Gleichgewicht und die Bewegung eines unendlich dünnen elastischen Stabes. J Für Die Reine Und Angew Math 1859:285–313. doi:10.1515/crll.1859.56.285. Reissner E. On transverse bending of plates, including the effect of transverse shear deformation. Int J Solids Struct 1975;11:569–73. doi:10.1016/0020-7683(75)90030-X. Mindlin RD. Influence of rotatory inertia and shear in flexural motions of isotropic elastic plates. ASME J Appl Mech 1951;18:1031–6. Liu D, Li X. An Overall View of Laminate Theories Based on Displacement Hypothesis. J Compos Mater 1996;30:1539–61. doi:10.1177/002199839603001402. Carrera E. Historical review of Zig-Zag theories for multilayered plates and shells. Appl Mech Rev 2003;56:287. doi:10.1115/1.1557614. Ren JG. A new theory of laminated plates. Compos Sci Technol 1986:225–39. Ren JG. Bending theory of laminated plates. Compos Sci Technol 1986:225–48. Koiter WT. A consistent first approximations in the generaltheory of thin elastic shells. Proc Symp. Theory Thin Elastic Shells, 1959, p. 12–23. Sayyad AS, Ghugal YM. On the free vibration analysis of laminated composite and sandwich plates: A review of recent literature with some numerical results. Compos Struct 2015;129:177–201. doi:10.1016/j.compstruct.2015.04.007. Sayyad AS, Ghugal YM. Bending, buckling and free vibration of laminated composite and sandwich beams: A critical review of literature. Compos Struct 2017;171:486–504. doi:10.1016/j.compstruct.2017.03.053. Li X, Liu D. Generalized laminate theories based on double superposition hypothesis. Int J Numer Methods … 1997;40:1197–212. doi:10.1002/(SICI)10970207(19970415)40:7<1197::AID-NME109>3.0.CO;2-B. Zhen W, Lo SH, Sze KY, Wanji C. A higher order finite element including transverse normal strain for linear elastic composite plates with general lamination configurations. Finite Elem Anal Des 2012;48:1346–57. doi:10.1016/j.finel.2011.08.003. Zhen W, Li T. C0-type global-local higher-order theory including transverse normal thermal strain for laminated composite plates under thermal loading. Compos Struct 2013;101:157– 67. doi:10.1016/j.compstruct.2013.02.002. Parlevliet PP, Bersee HEN, Beukers A. Residual stresses in thermoplastic composites-A study of the literature-Part I: Formation of residual stresses. Compos Part A Appl Sci Manuf 2007;38:651–65. doi:10.1016/j.compositesa.2006.07.002. Sarrazin H, Kim B, Ahn S-H, Springer GS. Effects of processing Temperature and Layup on Springback. J Compos Mater 1995;29:1278–94. doi:10.1177/002199839502901001. Nyman T, Rudlund M, Bohlin J, Berg R, Engineering V, Sky C. Shape Distortion Analysis of a Complex Shaped Wing 2015:19–24. Pagano NJ. Exact Solutions for Composite Laminates in Cylindrical Bending. J Compos Mater 1969;3:398–411. doi:10.1177/002199836900300304.

VITAE André Schwanz de Lima BS in Mechanical Engineering (UFES, 2014), sandwich graduation at École des Mines d’Alès (2012-2013), MSc in Aeronautical and Mechanical Engineering (ITA, 2016). PhD student of the department of mechanical engineering of ITA, Solid Mechanics and Structures field. Alfredo Rocha de Faria BS and MSc in Mechanical Engineering (ITA, 1993, 1995), PhD in Aerospace Science & 24

Engineering (University of Toronto, 2000). Associate Professor in the department of mechanical engineering of ITA. Responsible for teaching undergraduate and graduate courses. CNPq fellowship holder (PQ 1D). Publications: 30 papers in refereed international journals, 44 papers in refereed conferences, and 1 book chapter. Expert in aerospace engineering on topics such as design of aerospace structures, finite elements, structural optimization and composite materials.

FIGURE CAPTIONS Fig. 1. Lamina with structural axes (xyz), material axes (123) and orientation angle θk. Fig. 2. Stack-up sequence in a laminate. Fig. 3. Axial and transverse displacement distributions. Fig. 4. Proposed beam finite element and possible boundary conditions. Fig. 5. Simply supported beam under uniform load. Fig. 6. Elasticity and finite element results for the metallic sandwich beam. Fig. 7. Elasticity and finite element results for the 0°/90°/0° composite laminate. Fig. 8. Normalized maximum deflection by S ratio (a) metallic sandwich beam, (b) 0°/90°/0°

laminate. Fig. 9. Cantilever beam with tearing load. Fig. 10. w displacement field for a cantilever beam under tearing load by the (a) present formulation, (b) Nastran solution. Fig. 11. σz stress field for a cantilever beam under tearing load by the (a) present formulation, (b) Nastran solution.

TABLES Table 1 Computation of u and u,x u − 1 ˆ (τ τ + τ τ ) 1 Q 55 xzb xzb xzt xzt −1 ˆ ub Q u

u,x −

−

55 b

ut

ˆ −1u Q 55 t

−

ub,x

−

ut,x

−

ˆ −1u Q 55 b ˆ −1u Q

ub,xx

ˆ −1[u + Q Q ˆ −1 ˆ −1 Q 55 b,xx 55 33 (ub,x + Q13Q55ub )] ˆ −1[u + Q Q ˆ −1(u + Q Q ˆ −1u )] Q

− − − −

wt,x

− − −1 ˆ ˆ −1w ) Q55(wb,x + Q55Q 33 b −1 ˆ ˆ Q (w + Q Q−1w )

wb,xx

−

wt,xx

−

ˆ −1(w + Q Q ˆ −1 Q 55 b,x 55 33wb ) ˆ −1(w + Q Q ˆ −1w ) Q

ut,xx wb wt wb,x

55

t ,xx

55

55

33

t ,x

55 t

t,x

55

13

33

−

55 t

−

t

55

25

t ,x

55

33

t

wb,xxx wt,xxx

ˆ −1Q Q ˆ −1 ˆ −1 ˆ −1 Q 55 55 33[wb, xx + Q13Q55(wb,x + Q55Q33wb )] ˆ −1Q Q ˆ −1[w + Q Q ˆ −1(w + Q Q ˆ −1w )] Q 55

55

33

t ,xx

13

55

t,x

55

33

− −

t

Table 2 Computation of w and w,x w −1 ˆ 1 Q33 (σ zbσ zb + σ ztσ zt ) ub − ut − −1 ˆ ˆ −1u ) ub,x Q33(ub,x + Q13Q 55 b − 1 ˆ (u + Q Q ˆ −1u ) ut,x Q 33

t,x

13

w,x −

− − − −

55 t

ˆ −1(u + Q Q ˆ −1 Q 33 b,x 13 55ub ) ˆ −1(u + Q Q ˆ −1u ) Q

ub,xx

−

ut,xx

−

wb

−

wt

ˆ w Q b −1 ˆ Q w

wb,x

−

wt,x

−

ˆ −1w Q 33 b ˆ Q−1w

wb,xx

−

wt,xx

ˆ [w + Q Q ˆ (w + Q Q ˆ w )] Q b, xx 13 b,x 55 b − 1 − 1 − 1 ˆ [w + Q Q ˆ (w + Q Q ˆ w )] Q

wb,xxx

−

wt,xxx

−

ˆ −1[w + Q Q ˆ −1 ˆ −1 Q 33 b, xx 13 55(wb, x + Q55Q33wb )] ˆ −1[w + Q Q ˆ −1(w + Q Q ˆ −1w )] Q

33

t,x

−1 33

33

−1 33

33

13

55

33

−1 33

t,x

55 t

−

t

−1 55

t ,xx

13

55

33

t

−

t

33

t ,xx

26

13

55

t,x

55

33

t

0.5

0.5

ELST FEM

ELST FEM

S=4

0.25

S=20

0.25

z

z 0

0

-0.25

-0.25

-0.5

-1

-0.75 -0.5 -0.25

0

0.25

εx

0.5

0.75

-0.5

1

-1

-0.75 -0.5 -0.25

(a)

0

0.25

εx

0.5

0.75

1

(b) 0.5

0.5

ELST FEM

ELST FEM

S=4

0.25

S=20

0.25

z

z 0

0

-0.25

-0.25

-0.5

-1

-0.75 -0.5 -0.25

0

εz

0.25

0.5

0.75

-0.5

1

-1

-0.75 -0.5 -0.25

(c)

0

0.25

εz

0.5

0.75

1

(d) 0.5

0.5

ELST FEM

ELST FEM 0.25

0.25

z

z S=4

0

-0.25

-0.25

-0.5

S=20

0

0

0.25

0.5

γxz

0.75

1

1.25

-0.5

0

0.25

0.5

γxz

0.75

1

1.25

(f)

(e) 0.5

0.5

ELST FEM

ELST FEM

S=4

0.25

S=20

0.25

z

z 0

0

-0.25

-0.25

-0.5

-0.5 -1

-0.75 -0.5 -0.25

0

σx

(g)

0.25

0.5

0.75

1

-1

-0.75 -0.5 -0.25

0

σx

(h)

0.25

0.5

0.75

1

0.5

0.5 ELST FEM FEM-I

S=4

0.25

S=20

0.25

ELST FEM FEM-I

z

z 0

0

-0.25

-0.25

-0.5

0

0.25

0.5

0.75

σz

-0.5

1

0

0.25

0.5

0.75

σz

1

(j)

(i) 0.5

0.5 ELST FEM

ELST FEM

S=4

0.25

S=20

0.25

z

z 0

0

-0.25

-0.25

-0.5

0

0.25

0.5

0.75

τxz (k)

1

-0.5

1.25

0

0.25

0.5

0.75

τxz (l)

1

1.25

0.5

0.5

ELST FEM

ELST FEM

S=4

0.25

S=20

0.25

z

z 0

0

-0.25

-0.25

-0.5

-1

-0.75 -0.5 -0.25

0

u

(m)

0.25

0.5

0.75

1

-0.5

-1

-0.75 -0.5 -0.25

0

u

(n)

0.25

0.5

0.75

1

0.5

0.5

ELST FEM

ELST FEM

S=4

0.25

S=20

0.25

z

z 0

0

-0.25

-0.25

-0.5

-0.5 -1

-0.75 -0.5 -0.25

0

εx

0.25

0.5

0.75

1

-1

-0.75 -0.5 -0.25

0

εx

0.25

0.5

0.75

1

(b)

(a)

0.5

0.5

ELST FEM

ELST FEM

S=4

0.25

S=20

0.25

z

z

0

0

-0.25

-0.25

-0.5

-0.5 -1

-0.75 -0.5 -0.25

0

εz

0.25

0.5

0.75

1

-1

-0.75 -0.5 -0.25

0

εz

0.25

0.5

0.75

1

(d)

(c)

0.5

0.5

ELST FEM

ELST FEM

S=4

0.25

S=20

0.25

z

z

0

0

-0.25

-0.25

-0.5

-0.5

0

0.25

0.5

γxz (e)

0.75

1

1.25

0

0.25

0.5

γxz (f)

0.75

1

1.25

0.5

0.5

ELST FEM

ELST FEM

S=4

0.25

S=20

0.25

z

z 0

0

-0.25

-0.25

-0.5

-1

-0.75 -0.5 -0.25

0

0.25 0.5 0.75

σx

-0.5

1

-1

-0.75 -0.5 -0.25

0

σx

0.25 0.5 0.75

1

(h)

(g) 0.5

0.5 ELST FEM FEM-I

S=4

0.25

S=20

0.25

ELST FEM FEM-I

z

z

0

0

-0.25

-0.25

-0.5

-0.5 0

0.25

0.5

0.75

σz

1

0

0.25

0.5

0.75

σz

1

(j)

(i) 0.5

0.5

ELST FEM

ELST FEM

S=4

0.25

S=20

0.25

z

z 0

0

-0.25

-0.25

-0.5

0

0.25

0.5

τxz

0.75

1

1.25

-0.5

0

0.25

0.5

(k)

τxz

0.75

1

1.25

(l)

0.5

0.5 ELST FEM

ELST FEM

S=4

0.25

S=20

0.25

z

z 0

0

-0.25

-0.25

-0.5

-1

-0.75 -0.5 -0.25

0

u

0.25

(m)

0.5

0.75

1

-0.5

-1

-0.75 -0.5 -0.25

0

u

0.25

(n)

0.5

0.75

1

2

4 ELST FEM

ELST FEM

1.75

3

__

__

wmax

wmax

1.5

2

1.25

1

1

0

10

20

S (a)

30

40

0

0

10

20

S (b)

30

40