10th HSTAM International Congress on Mechanics Chania, Crete, Greece, 25 – 27 May, 2013

1

MODAL SHAPE CONTROL OF SMART COMPOSITE BEAMS USING PIEZOELECTRIC ACTUATORS

Georgia A. Foutsitzi1, Evangelos P. Hadjigeorgiou2, Christos G. Gogos1 and Georgios E. Stavroulakis3

1Technological Educational Institution of Epirus, Department of Finance and Auditing,

P.O. BOX 169, GR-48100 Preveza, Greece e-mail: gfoutsi@teiep.gr, web page: http://preveza.teiep.gr/faculty/gfoutsi/

e-mail: cgogos@teiep.gr, web page: http://preveza.teiep.gr/faculty/chgogos/

2University of Ioannina Department of Materials Science and Engineering,

GR-45110 Ioannina, Greece e-mail: ehadjig@cc.uoi.gr

3Technical University of Crete,

Dept. of Production Engineering and Management Institute of Computational Mechanics and Optimization

University Campus, Kounoupidiana, GR-73100 Chania, Greece e-mail: gestavr@dpem.tuc.gr http://users.isc.tuc.gr/~gestavroulakis/

Keywords: smart structures, piezoelectric laminated element, modal control.

Abstract. A new finite element formulation for laminated structures with bonded piezoelectric patches that can also handle coupled analysis is developed. The developed FE model can handle smart patches embedded in asymmetrically stacked laminated composite beams. These developed elements are used to perform dynamic analysis and modal shape control to a cantilever smart beam. Numerical results are presented to show the optimal locations and optimal values of electrical voltages in the actuators to achieve modal shape. 1 INTRODUCTION

Smart laminated composite beams with bonded/embedded piezoelectric layers as sensors and/or actuators offer excellent potentials to be used in advanced structural applications including aerospace and automotive industries. Excellent sensing and actuating capabilities of piezoelectric materials made them the most practical smart materials to integrate with laminated structures. Analysis of smart laminated structures requires an accurate and efficient model with capability to take into account the strong inhomogeneities through the thickness of the structure. A vast number of analytical and computational models for smart piezoelectric structures using various theories have been reported in the review papers by Sunar and Rao [1], Saravanos and Heyliger [2]. Khdeir et al. [3] proposed analytical models and solutions for the free vibration of cross-ply laminated beams with extension piezoelectric actuators. Apart from the analytical modeling of smart laminated beams, finite element modeling techniques have also been developed [4], [5].

The recent advances in smart composite structures have prompted interest in structural shape control [6]. Static shape control is focused on finding the optimal values of some parameters, like location and applied voltages of actuators for achieving the structural shape [7], [8], [9].

In this work, the super-convergent FE approach is used to derive beam elements with embedded piezoelectric patches for both uncoupled and coupled material models. The element uses higher order interpolating polynomials that are derived by solving the static part of the governing equations of motion and hence gives an exact elemental stiffness matrix. Each node has three degrees of freedom (dof), which include extension, bending and rotation. First-order theory is used for modeling transverse shear deformation. The formulated element is an extension of the two-node Timoshenko beam element of Ref. [10] to incorporate the extension degree of freedom and to handle smart piezoelectric patches embedded in asymmetrically stacked laminated composite beams. The developed FE model is validated by comparing the free vibration response to those available in literature. A cantilever cross-ply laminated beam is used for modal shape control applications. The optimal values for the locations of the piezo-actuators are determined and optimal voltages for shape control are obtained for cantilever beams by using a genetic optimization procedure developed in [9].

Georgia A. Foutsitzi, Evangelos P. Hadjigeorgiou, Christos G. Gogos and Georgios E. Stavroulakis

2

2 GOVERNING DIFFERENTIAL EQUATIONS OF A LAMINATED BEAM

Consider a laminated beam as shown in Figure 1. The laminate is formed from two or more layers bonded together to act as a single layer material with piezoelectric sensor and actuator layers. The bond between two layers is assumed to be perfect. The classical formulation of laminated materials is followed and complemented with electromechanical coupling terms. It is assumed that the beam centroidal and material axes coincide with the x-axis. The piezoelectric layers have poling direction along z-axis and the electric field is applied through the thickness direction. The displacement field equations for the beam using first order shear deformation at any point through the thickness are presented as

( ) ( ) ( ) ( ) ( ) ( )0 0, , , , , , , 0, , , , ,,x y y zu x y z t z x t u x y z u x y z t w x tu x t θ= − = = (1) where t denotes time; u0 and w0 denote the axial and transverse displacements of the beam mid–plane, respectively and θy is the rotation of the beam cross section about the positive y-axis. The electric field E in the pth piezoelectric layer is expressed as

{ } { } ( )1, , 0 0 ,T

p p p p px y zp p p

E E E E B x thφ φ φ⎡ ⎤⎡ ⎤= = − = −⎣ ⎦ ⎢ ⎥⎣ ⎦ (2)

where ph is the thickness and ( ),p x tφ is the electrical voltage of the pth piezoelectric layer, respectively. For a one-dimensional beam where the width in the y-direction is stress free and by using the plane stress assumption, the general 3D constitutive equations can be reduced to:

σ x

τ xzDz

⎧

⎨⎪⎪

⎩⎪⎪

⎫

⎬⎪⎪

⎭⎪⎪k

=

!Q11 0

0 !Q55

00

!e31 0 0

⎡

⎣

⎢⎢⎢⎢⎢

⎤

⎦

⎥⎥⎥⎥⎥k

ε xγ xz0

⎧

⎨⎪⎪

⎩⎪⎪

⎫

⎬⎪⎪

⎭⎪⎪

−!e310

− !ξ33

⎧

⎨⎪⎪

⎩⎪⎪

⎫

⎬⎪⎪

⎭⎪⎪k

Ezk (3)

where ,x xzσ τ is the normal and shear stress, respectively, ,x xzε γ is the normal and shear strain respectively,

zD is the electric displacement, zE is the electric field and  the  coefficients   !Q11, !Q55, !e31  and  !ξ33  are given in

[9]. The total strain energy U and the kinetic energy T are calculated using the expressions

[ ]0

12

L

x x xz xz z zA

U E D dAdxσ ε τ γ= + −∫ ∫ T = 12

ρ !ux2 + !uz

2⎡⎣ ⎤⎦dAA∫ dx

0

L

∫ (4)

where dA is the area of cross-section of the beam.

Applying Hamilton’s principle, we get the following four governing equations corresponding to 4 degrees of freedom

δu0 : I0!!u0 − I1!!θ y − A11u0,xx + B11θ y ,xx − Ae31p Ez ,x

p

p∑ = 0

δw0 : I0 !!w0 − A55 w0,xx −θ y ,x( ) = 0δθ y : I2 !!θ y − I1!!u0 + B11u0,xx − D11θ y ,xx + A55 θ y −w0,x( )+ Be31

p Ez ,xp

p∑ = 0

δEzp : −Ae31

p u0,x + Be31p θ y ,x − Aξ33

p Ezp = 0

(5)

The associated forced boundary conditions are

Georgia A. Foutsitzi, Evangelos P. Hadjigeorgiou, Christos G. Gogos and Georgios E. Stavroulakis

3

31

31

11 0, 11 ,

11 0, 11 ,

p px y x e z

p

p px y x e z

p

N A u B A E

M B u D B E

θ

θ

= − +

= − + −

∑

∑

( )55 0,x yV A w θ= − (6)

where N is the axial force, M is the moment and V is the shear force. The cross sectional properties

3111 55 11 11, , , , peA A B D A and

31

peB are given by

A11,B11,D11⎡⎣ ⎤⎦ = !Q11 1, z, z2⎡⎣ ⎤⎦

A∫ dA = b !Qij 1, z, z

2⎡⎣ ⎤⎦dzzk

zk+1

∫k=1

N

∑ A55 = k !Q55z2

A∫ dA = kb !Q55z

2 dzzk

zk+1

∫k=1

N

∑

[ ] 20 1 2, , 1, ,

A

I I I z z dAρ ⎡ ⎤= ⎣ ⎦∫ Aξ33p = !ξ33

p∫ dA Ae31

p ,Be31p⎡

⎣⎤⎦ = !e31 1, z⎡⎣ ⎤⎦

p∫ dA (7)

where k=5/6 is the shear coefficient. In the above expression, the integration over p means that this operation is performed only to pth piezoelectric layer.

3 SUPER-CONVERGENT FINITE ELEMENT FORMULATION The governing equations (5) and the associated boundary conditions (6) will be used for stiffness and electromechanical coupling matrix formulation. This is done by assuming appropriate polynomials for the displacement field based on the order of the static part of the differential equations and substituting these back into the governing equations. In this process certain constants can be eliminated and, at the same time, certain constants become dependent on material and sectional properties. The details of the formulation are given in the next section.

3.1 Uncoupled Formulation For finite element beam formulation based on uncoupled constitutive model, the electric field is considered as a known variable, hence ,

pz xE =0. Hence, the governing equations and the associated boundary conditions become

δu0 : I0!!u0 − I1!!θ y − A11u0,xx + B11θ y ,xx = 0

δw0 : I0 !!w0 − A55 w0,xx −θ y ,x( ) = 0δθ y : I2 !!θ y − I1!!u0 + B11u0,xx − D11θ y ,xx + A55 θ y −w0,x( ) = 0

(8)

and

( )

31

31

11 0, 11 ,

11 0, 11 ,

55 0,

p px y x e z

p

p px y x e z

p

x y

N A u B A E

M B u D B E

V A w

θ

θ

θ

= − +

= − + −

= −

∑

∑ (9)

It is to be noted that the electric field will appear only in the force boundary conditions.  The finite element formulation begins by assuming the interpolating functions of appropriate order for the six degrees of freedom ( 0 0, , yu w θ at the two nodes). Looking at the governing equations (8), we see that the axial displacement u0 and the rotation about the y-axis θy require quadratic polynomial, while the transverse displacement w0 requires cubic polynomials. Hence the interpolating polynomials for the three mechanical degrees of freedom can be assumed as

2 2 3 20 1 2 3 0 4 5 6 7 8 9 10, , yu c c x c x w c c x c x c x c c x c xθ= + + = + + + = + + (10)

where jc , j=1,2,...,10, are the ten unknown coefficients to be determined by the six nodal degrees of freedom. Hence, there are only six independent constants and four dependent constants, which can be obtained by

Georgia A. Foutsitzi, Evangelos P. Hadjigeorgiou, Christos G. Gogos and Georgios E. Stavroulakis

4

substituting Equations (10) into the static part of Equations (8). The remaining six independent constants are calculated using the boundary conditions (9) at the beam ends x=0, L. The resulting explicit form of the three shape functions are given by

{ } [ ]{ } [ ][ ] [ ] { }0 0, ,T T

y u wu w N d N N N dθθ ⎡ ⎤= = ⎣ ⎦ (11)

where { } { }1 1 1 2 2 20 0 0 0, , , , ,

T

y yd u w u wθ θ= is the array of nodal displacements and

[ ]

( ) { }

( ){ }

( ) { }

( ){ }

2

2

2

2

16131

6131

Tu

v L

vN

v L

v

ξβ ξ ξ

β ξ ξ

ξβ ξ ξ

β ξ ξ

−⎡ ⎤⎢ ⎥⎢ ⎥−

+⎢ ⎥⎢ ⎥⎢ ⎥−⎢ ⎥+⎢ ⎥=⎢ ⎥⎢ ⎥⎢ ⎥− −

+⎢ ⎥⎢ ⎥⎢ ⎥−⎢ ⎥+⎣ ⎦

, [ ]

( ) ( ){ }

( )

( ){ }

( )

3 2

3 2

3 2

3 2

01 2 3 11

2 11 2 2

01 2 31

11 2 2

Tw

vv

Lv

N

v

Lv

ξ ξ νξ

ν νξ ξ ξ

ξ ξ νξ

ν νξ ξ ξ

⎡ ⎤⎢ ⎥⎢ ⎥− − + +

+⎢ ⎥⎢ ⎥

⎧ ⎫⎛ ⎞ ⎛ ⎞⎢ ⎥− + + +⎨ ⎬⎜ ⎟ ⎜ ⎟⎢ ⎥+ ⎝ ⎠ ⎝ ⎠⎩ ⎭⎢ ⎥=⎢ ⎥⎢ ⎥⎢ ⎥− − −⎢ ⎥+⎢ ⎥

⎧ ⎫⎢ ⎥⎛ ⎞ ⎛ ⎞− − −⎨ ⎬⎜ ⎟ ⎜ ⎟⎢ ⎥+ ⎝ ⎠ ⎝ ⎠⎩ ⎭⎣ ⎦

,

[ ]

( ) { }

( ) ( ) ( ){ }

( ) { }

( ) ( ){ }

2

2

2

2

06

11 3 4 11

06

11 3 21

T

v L

vN

v L

v

θ

ξ ξ

ξ ν ξ ν

ξ ξ

ξ ν ξ

⎡ ⎤⎢ ⎥⎢ ⎥−

+⎢ ⎥⎢ ⎥⎢ ⎥− + + +⎢ ⎥+⎢ ⎥=⎢ ⎥⎢ ⎥⎢ ⎥− −

+⎢ ⎥⎢ ⎥⎢ ⎥− −⎢ ⎥+⎣ ⎦

(12)

where 11

11

BA

β = , 11 11 11 112

11 55

12 A D B BvL A A

⎛ ⎞−= ⎜ ⎟

⎝ ⎠ and x

Lξ = .

The stiffness matrix is obtained using forced boundary conditions (9), where ,pz xE terms are used to get the

equivalent nodal force due to applied voltage. Next, the consistent element mass matrix [M] is derived using the expression:

[ ] [ ] [ ]1 20

1 00 1 0

0

LnlayerT

kk A

zM N N dAdx

z zρ

=

−⎡ ⎤⎢ ⎥= ⎢ ⎥⎢ ⎥−⎣ ⎦

∑ ∫ ∫ (13)

It is noted that the formulated element can be used for actuator applications. The formulated element is expected to have super convergent property as it uses exact solution to the governing equation as its interpolation function. Also, good accuracy in dynamic analysis can be expected from this element using smaller system sizes. This is because the stiffness of the structure is exactly represented even though the inertial distribution is approximate. The performance of the formulated element is examined for free vibration and modal control problems.

3.1 Coupled Formulation

For finite element beam formulation based on coupled constitutive model, electric field terms pzE will be

considered as an unknown variable. Differentiating Equation (5)4 with respect to x and substituting to the static parts of Equations (5)1-(5)3 and to the boundary conditions (6), we get

Georgia A. Foutsitzi, Evangelos P. Hadjigeorgiou, Christos G. Gogos and Georgios E. Stavroulakis

5

δu0 : −A11* u0,xx + B11

*θ y ,xx = 0

δw0 : A55 w0,xx −θ y ,x( ) = 0δθ y : B11

* u0,xx − D11*θ y ,xx + A55 θ y −w0,x( ) = 0

(14)

N = A11* u0,x − B11

*θ y ,xM = −B11

* u0,x + D11*θ y ,x

V = A55 w0,x −θ y( ) (15)

where

31 31

33

*11 11

p pe ep

A AA A

Aξ= −∑ , 31 31

33

*11 11

p pe ep

p

A BB B

Aξ= −∑ , 31 31

33

*11 11

p pe ep

p

B BD D

Aξ= −∑

The formulation of this element is proceed in a similar manner as the uncoupled case and the stiffness and mass matrices and the electric induced load vectors take the same form as uncoupled case with A11;B11 and D11 replaced by * *

11 11,A B and *11D .

After calculating the mechanical displacement, the electric field can be obtained using Equation-(5)4. It is noted that this element can be used for sensing applications.

4 MODAL SHAPE CONTROL

The concept of this work is similar to the quasi-static shape control problem. The main difference in modal shape control is that the structural shape is defined based on the i-th mode shape (eigenvector) Vi. and scalar factor G as d

iw =G*Vi.

4.1 Problem Formulation

The shape control problem considered in this section focuses on piezoelectric actuator design optimization in terms of finding optimum values for applied voltages and actuator positions. When the location of piezoelectric actuators α and the applied voltage vector φ are considered as design variables, the quasi-static shape control problem can be defined, in the context of optimization formulation, as follows: Find (φ ,α) to minimize

( )( )21

,r

di i

if w wφ α

=

= −∑ (16)

Subject to

min maxiφ φ φ≤ ≤ , 0 1ia or= , iia number of actuators=∑ (17)

where diw is the desired nodal displacement value and r is the number of concerned displacements. In this work,

the actuator position is modeled using a Boolean type discrete variable for each element and the electric potential of the actuator using a bounded continuous variable. The Mixed Integer Problem that arises is highly nonlinear and is solved using a modified genetic algorithm procedure developed in [9].

5 NUMERICAL EXAMPES

5.1 Validation Example

In order to validate the present finite element model, a composite beam of the following dimensions is considered [3]: L = 0.06 m and h = 0.006 m. It is assumed that the beam is made of 6 layers: 2 piezoelectric actuators layers and 4 composite material layers (PZT5H / composite / composite / composite / composite / PZT5H). The thickness of each piezoelectric layer is tp = 0.001 and the thickness of each composite layer is tc = 0.001. The ply stacking sequence of the beam is (0°/ 0°/ 90°/ 0°/ 90°/0°). In Figure 1 the location of piezoelectric actuator is d1 = 0.02 m. The piezoelectric patches are bonded at local position on the beam surface and have length Lp = 0.02m. .

Georgia A. Foutsitzi, Evangelos P. Hadjigeorgiou, Christos G. Gogos and Georgios E. Stavroulakis

6

There is no mechanical load in structure. The composite has the following material constants for each layer: E1=132 GPa, E2=10.8 GPa, G12 =5.65 GPa, G13= 5.65GPa, G23= 3.38 GPa, v12=0.24, ρc= 1540 kg/m3. The properties of the piezoelectric actuators include the elastic constants: c11=c22=127.2GPa, c33=117.44 GPa, c12= 80.2 GPa, c13=c23=84.67 GPa, c44=c55=23 GPa, c66=23.5 GPa, the density ρ=7500 kg/ m3 and piezoelectric constants e13= -6.55 C/ m2, e15= 17 C/ m2. The beam is divided evenly into 30 finite elements. The left end of the beam has been clamped.

The dynamic analysis of the beam is performed and the first four circular frequencies are listed in table 1. An excellent agreement is obtained when compared with the results obtained by [3], where exact solutions were obtained by using the state-space approach. Figure 2 depicts the first four vibration modes of the beam. The shapes are similar with that obtained by [3].

Mode No Analytical Results [3] Present First 976.7 976.8

Second 6098.9 6100.7 Third 15076.0 15096.5 Fourth 29450.2 29454.3

Table 1.The first four natural frequencies (Hz) of the smart laminated beam of Fig.1

 Number

of Elements

Mode

1 2 3 4

1 -41.77 -84.25 -800.00 -649.10 2 0 0 0 0 3 -20.24 0 0 0 4 23.27 15.92 0 0 5 0 0 0 799.99 6 -113.43 0 0 0 7 0 0 502.13 0 8 0 55.87 218.54 0 9 0 0 366.33 0

10 0 0 749.69 0 11 0 0 0 799.97 12 0 0 0 0 13 0 40.63 0 -669.47 14 0 0 0 0 15 0 0 173.19 0 16 0 65.52 0 -783.69 17 0 0 -799.83 0 18 0 0 0 -799.97 19 0 0 0 565.61 20 -3.08 0 -799.96 0 21 -7.94 34.02 0 0 22 -62.31 0 -800.00 799.95 23 62.47 0 0 800.00 24 0 14.49 0 0 25 0 110.61 0 0 26 0 -91.03 0 0 27 -9.89 0 0 0 28 0 0 0 0 29 0 0 0 0 30 4.26 44.39 -799.98 799.83

Fitness 2.02e-16 2.58e-16 2.09e-13 1.38e-12

Table 2. Optimal location and voltages of actuators within the 30 finite element mesh

Georgia A. Foutsitzi, Evangelos P. Hadjigeorgiou, Christos G. Gogos and Georgios E. Stavroulakis

7

5.2 Modal Shape Control

A beam with similar material and geometric properties as described in section 5.1, is considered to calculate the optimal location and applied voltage of actuators in order to modify its shape. Next, all the thirty elements are candidates for locations of 10 piezoelectric actuator patches of length equal to FE mesh. In addition, the lower as well as the upper piezoelectric patches acts as actuators. The left end of the beam has been clamped. The lower limit of the voltage is set to be -800 V and the upper limit is set to be 800 V (limit imposed due to depoling of actuators). The desired shape is given by d

iw =G*Vi, where Vi is the i-th mode shape (eigenvector) and scalar factor of modal shape is G=10-5 for the mode shape 1 and 2 and is G=10-4 for the mode shape 3 and 4. Table 2 shows the optimal solutions for placement of the actuators and the corresponding optimal voltages for the first four mode shapes. The genetic algorithms were run using the following parameters: Generations=500, Population=100, EliteCount=4.

Figure 2 shows the desired modal shapes and the calculated shape in each mode. The results in Table 1 shows that it is easy to achieve the modal shape of modes 1 and 3 with low applied voltages, whilst modal shapes of modes 2 and 4 require very high energy to achieve. It is obvious that more actuators are required to achieve the fourth mode since the desired structural shape is more complex.

Figure 1. The geometry of a beam with piezoelectric patches

0 0.01 0.02 0.03 0.04 0.05 0.060

1

2

3

4

5x 10-7 First Mode

0 0.01 0.02 0.03 0.04 0.05 0.06-2

-1.5

-1

-0.5

0

0.5

1x 10-7 Second Mode

0 0.01 0.02 0.03 0.04 0.05 0.06-5

0

5

10x 10-7 Third Mode

0 0.01 0.02 0.03 0.04 0.05 0.06-8

-6

-4

-2

0

2

4

6

8x 10-7 Fourth Mode

desired shapeoptimal

Georgia A. Foutsitzi, Evangelos P. Hadjigeorgiou, Christos G. Gogos and Georgios E. Stavroulakis

8

Figure 2. Desired and calculated shapes at different modes of the free vibration

6 CONCLUSIONS

A finite element formulation considering electromechanical effects has been developed for axial-flexural-shear coupled deformation in asymmetrically stacked laminated composite beams with piezoelectric patches. The developed elements are used to perform dynamic analysis and modal shape control to a cantilever smart laminated beam. Numerical simulation of composite beam structures with piezoelectric actuators is conducted and presented. The optimization problem of finding the optimal location and voltages of actuators is done by using the hybrid genetic algorithm developed in authors’ previous work [9]. The proposed model can be developed for vibration energy harvesting applications as a future work. Acknowledgements: This research has been co-financed by the European Union (European Social Fund – ESF) and Greek national funds through the Operational Program "Education and Lifelong Learning" of the National Strategic Reference Framework (NSRF) - Research Funding Program: ARCHIMEDES III. Investing in knowledge society through the European Social Fund. The authors gratefully acknowledge this support.

REFERENCES

[1] Sunar, M. and Rao, S. (1999) “Recent advances in sensing and control of flexible structures via piezoelectric materials technology”, Appl. Mech. Rev. 52, pp. 1–16.

[2] Saravanos, D. and Heyliger, P.R. (1999) “Mechanics and computational models for laminated piezoelectric beams, plates, and shells” Appl. Mech.Rev. 52, pp. 305–320.

[3] Khdeir, A.; Darraj, E. and Aldraihem, O. J. (2012) “Free Vibration of Cross Ply Laminated Beams with Multiple Distributed Piezoelectric Actuators” J of Mechanics 28 (1), pp. 217-227.

[4] Aldraihem, O. J., Wetherhold, R. C. and Singh, T. (1997) “Distributed control of laminated beams : Timoshenko Vs. Euler-Bernoulli Theory,” J. of Intelligent Materials Systems and Structures 8, pp. 149–157.

[5] Benjeddou,A., Trindade, M. A. and Ohayon, R. (1999) “New shear actuated smart structure beam finite element,” AIAA J., vol. 37, pp. 378–383, 1999.

[6] Irschik, H. (2002) “A review on static and dynamic shape control of structures using piezoelectric actuation” Comput. Mech. 26, pp. 115–128.

[7] Hadjigeorgiou E. P., Stavroulakis G. E., Massalas C. V. (2006) “Shape control and damage identification of beams using piezoelectric actuation and genetic optimization” Int J Eng Sci 44, 7, pp. 409-421

[8] Tong D., Williams R. L. and Agrawal S. K. (1998) “Optimal Shape Control of Composite Thin Plates with Piezoelectric Actuators” Journal of Intelligent Material Systems and Structures 9(6), pp. 458-467.

[9] Foutsitzi, G., Gogos, C., Hadigeorgiou, E and Stavroulakis, G. (2012) “Actuator Location and Voltages Optimization for Shape Control of Smart Beams Using Genetic Algorithms” Submitted.

[10] Friedman, Z. and Kosmatka, J.B. (1993) “An improved two-node Timoshenko beam finite element” Comput. Struct. 47, pp. 473–481.