¨ KAISERLSLAUTERN TECHNISCHE UNIVERSIT AT D EPARTMENT OF M ATHEMATICS

MASTER’S THESIS

Finite element methods for geometrically linear curved beams

by: Tabitha Wangari Kinyanjui

Supervisors: Dr - ing. Joachim Linn (ITWM) Dr. Joseph M. Maubach (TU/e)

Kaiserslautern, Germany November 22, 2007

The Faculty members appointed to examine the thesis by Tabitha W. Kinyanjui find it satisfactory and recommend that it be accepted.

TU-Kaiserslautern

TU-Eindhoven

Date:

Date:

Contents 1 Introduction

8

2 Straight Beam Theory

2.1

2.2

11

2.0.1

Geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11

2.0.2

Stresses . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11

Euler - Bernoulli Beam Theory . . . . . . . . . . . . . . . . . . . . .

12

2.1.1

Total potential energy . . . . . . . . . . . . . . . . . . . . . .

13

Timoshenko Beam Theory . . . . . . . . . . . . . . . . . . . . . . .

14

2.2.1

15

Total potential energy . . . . . . . . . . . . . . . . . . . . . .

3 Finite Element Method 3.0.2 3.1

3.2

16

Formulation of the Finite Element Method . . . . . . . . . .

17

Bubnov-Galerkin’s Approximation Method . . . . . . . . . . . . . .

17

3.1.1

The Stiffness Matrix and the Load vector . . . . . . . . . . .

19

Assembly of global stiffness matrix and force vector . . . . . . . . .

21

4 Linear Curved Beams

22

4.1

The Equilibrium Equations . . . . . . . . . . . . . . . . . . . . . . .

24

4.2

Normal stresses in Curved Beams . . . . . . . . . . . . . . . . . . .

25

4.3

Deflection of a Curved Beam . . . . . . . . . . . . . . . . . . . . . .

28

4.4

Shear Deformation . . . . . . . . . . . . . . . . . . . . . . . . . . .

30

4.5

Finite Element Implementation . . . . . . . . . . . . . . . . . . . . .

32

4.5.1

34

Stiffness Matrix and Load Vector . . . . . . . . . . . . . . .

1

Contents 4.5.2

Timoshenko beam stiffness matrix . . . . . . . . . . . . . . .

34

4.5.3

Euler-bernoulli beam stiffness matrix . . . . . . . . . . . . .

36

4.6

Numerical Integration . . . . . . . . . . . . . . . . . . . . . . . . . .

39

4.7

Displacement Transformation Matrix . . . . . . . . . . . . . . . . . .

39

4.8

Common Boundary Conditions . . . . . . . . . . . . . . . . . . . . .

42

4.8.1 4.9

Numerical results of beams under different boundary conditions 44

Shear and Membrane Locking . . . . . . . . . . . . . . . . . . . . .

49

4.10 Numerical Test Examples . . . . . . . . . . . . . . . . . . . . . . . .

50

4.10.1 Example 1: An arc for various subtended angles. . . . . . . .

51

4.10.2 Castigliano’s Energy Method . . . . . . . . . . . . . . . . . .

53

4.10.3 Example 2 :A quarter-Circular Cantilever Beam . . . . . . . .

54

4.10.4 Results for Timoshenko curved Beam . . . . . . . . . . . . .

57

4.10.5 Results for the Euler-Bernoulli curved beam . . . . . . . . . .

58

4.11 Points to Note: Timoshenko curved beam strains definition . . . . . .

59

4.12 Summary and Conclusion . . . . . . . . . . . . . . . . . . . . . . . .

62

Appendix

64

.1

Matlab Programmes . . . . . . . . . . . . . . . . . . . . . . . . . . .

64

.1.1

Common programmes and Parameters for both beam theories

64

.1.2

Curved Beam - Euler-Bernoulli beam theory . . . . . . . . .

66

.1.3

Curved Beam - Timoshenko Beam Theory . . . . . . . . . .

71

Declaration

79

2

List of Figures 2.1

Geometry of a simple Beam . . . . . . . . . . . . . . . . . . . . . .

11

2.2

Deformation of cross-section of an Euler - Bernoulli beam . . . . . .

13

2.3

Deformation of the cross-section of a Timoshenko beam . . . . . . .

15

4.1

Bending stress in curved beams . . . . . . . . . . . . . . . . . . . . .

23

4.2

Forces and displacements on an infinitesimal curved beam element . .

24

4.3

Geometry of a beam’s element vertical displacement of magnitude y .

25

4.4

(a.)Tangential displacement (b.)Radial displacement(c.)cross-sectional displacements and rotations of the cross-sectional plane around the yaxis (d)About the z-axis. . . . . . . . . . . . . . . . . . . . . . . . .

4.5

29

Geometry of a curved beam element with three degrees of freedom at each node. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

33

4.6

curved beam element related to local and global displacements . . . .

40

4.7

straight beam element displacement . . . . . . . . . . . . . . . . . .

40

4.8

Relation of curved beam element in local and global coordinates system 41

4.9

Different beam supports and their reactions . . . . . . . . . . . . . .

44

4.10 Tangential displacement and Radial displacements for Euler -Bernoulli and Timoshenko beams fixed at both ends . . . . . . . . . . . . . . .

45

4.11 Cross sectional Rotation for Euler -Bernoulli and Timoshenko beams fixed at both ends . . . . . . . . . . . . . . . . . . . . . . . . . . . .

46

4.12 Tangential displacement and Radial displacements for cantilever Euler -Bernoulli and Timoshenko beams subjected to external load of 500 .

46

4.13 Cross sectional Rotation for cantilever beam Euler -Bernoulli and Timoshenko beams,subjected to external load of 500

3

. . . . . . . . . . .

47

List of figures 4.14 Tangential Displacement and Rotational Displacement for increased external load from 100 - 500 of a beam clamped at both ends . . . . .

47

4.15 Cross sectional Rotation for increased external load from 100 - 500 of a beam clamped at both ends . . . . . . . . . . . . . . . . . . . . . .

48

4.16 Displacement of an Euler- Bernoulli Beam and a Timoshenko Beam that are both simply supported at both ends,subjected to external load from 100 - 500. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

48

4.17 Cross sectional Rotation of an Euler- Bernoulli Beam and a Timoshenko Beam that are both simply supported at both ends, subjected to external load from 100 - 500. . . . . . . . . . . . . . . . . . . . .

49

4.18 Effects of reduced membrane energy term integration for Euler-Bernoulli beam, and reduced shear and membrane energy terms integration for Timoshenko beam . . . . . . . . . . . . . . . . . . . . . . . . . . . .

50

4.19 An arc fixed on one end and supported on rollers and subjected to point load on the other end . . . . . . . . . . . . . . . . . . . . . . . . . .

51

4.20 Radial displacement (a)With reduced integration, (b) Without reduced integration for an arc fixed on one end and supported on rollers and subjected to point load on the other end . . . . . . . . . . . . . . . .

52

4.21 A quarter-circular cantilever beam subjected to radial point load Q, at the free end. (Lee, Sin [22])

. . . . . . . . . . . . . . . . . . . . . .

54

4.22 Free body diagram of a quarter-circular cantilever beam subjected to radial point load Q, and we introduce axial dummy force Θ and a dummy moment - m at the free end. . . . . . . . . . . . . . . . . . .

55

4.23 Maximum Absolute Error for Tangential, Radial displacements and Cross sectional Rotation when the conventional strains of a Timoshenko curved beam are used . . . . . . . . . . . . . . . . . . . . . . . . . .

4

61

List of Tables 4.1

Parameters used for a curved beam subtended by angle π . . . . . . .

4.2

Comparison of the finite element radial displacement with reduced in-

44

tegration for the Timoshenko beam and Euler-Bernoulli beam with, the method by Saffari and Tabatabaei [2] 4.3

. . . . . . . . . . . . . . . . .

53

Comparison of the two different finite element solution models loaded with tip radial force, of a quarter circular cantilever beam with Castigliano’s energy solutions for Timoshenko beam. . . . . . . . . . . .

4.4

57

Comparison of the two different finite element solution models loaded with tip radial force, of a quarter circular cantilever beam with Castigliano’s energy solutions for Euler-Bernoulli Beam. . . . . . . . . .

4.5

58

Comparison of the two different strain definition for fem solution models loaded with tip radial force, of a quarter circular cantilever beam with Castigliano’s energy solutions for Timoshenko Beam. . . . . . .

5

60

6

Preface I would like to express my sincere thanks and appreciations to my supervisor Dr - ing. Joachim Linn for his great effort and a lot of time spent in introducing me into this area of study and guiding me through. Without his undoubted advice and guidance, it would have been impossible to get through this work. I would like to thank my lecturers at TU - Kaiserslautern and TU - Eindhoven. They kept us very busy, we barely noticed days passing by, but they tireless helped and assisted us all the way. I would like to thank the consortium for industrial mathematics for giving me this chance to participate in the Erasmus Mundus programme. Indeed it is a privilege to be a participant and the knowledge, and skills acquired during my two years of study are without doubt enormous. Many thanks to my friends in Eindhoven, and in Kaiserslautern. I wish to mention especially Changgih, Henry and Vincent. Thanks for keeping the long hours in the discussion groups. Lastly, i thank my parents and my siblings, for their tireless and relentless love, continuous support, and countless sacrifices they’ve had to make on my behalf. Thank you for believing in me. Tabitha W. Kinyanjui November 22, 2007 Kaiserslautern

Chapter 1

Introduction Beams are the most common structural components found in most mechanical structures. They can be defined as thin long structures that are capable of carrying loads in flexure, and are characterized by (a) The longitudinal dimension, which defines the axial direction of the beam is considerably larger than the transverse dimensions, the dimension defining the directions normal to the axial directions. (b) If the beam is restricted to deform in plane, then it resists primarily loading applied in one plane and has a cross section that is symmetric with respect to that plane, hence it can be analyzed with two-dimensional idealizations [1]. A variety of practical engineering problems for example some of the structures such as tyres, pipelines, rings, arc-like structures, just to mention a few however are modeled by curved finite elements. A lot can be learned by studying the performance of simple (essentially one-dimensional) models so as to gain insight into the more general behavior of the overall structures. The circular ring serves as a simple yet realistic model in addition for analyzing the effects of symmetry in structural systems elements. Since the deformations are symmetric, a closed ring can be regarded as an example of a curved beam with restrained sections or ends. From conditions of symmetry, the distribution of stress in one quadrant is know to be the same as in another. Thus one quadrant may be considered to be the curved beam in which the behavior of the ring under load, at any section is to be found. Furthermore, curved beams are known to be more efficient in transfer of loads than straight beams because the transfer is affected by bending, shear, and membrane action (Saffari, Tabatabaei [2]). Therefore the main objective in this thesis is 8

(a) Give a detailed derivation of equations of a curved beam in three dimesions, with reference to the quoted literature, having first introduced a straight beam model as a starting ground. (b) To investigate the behavior of a plane curved beam, resulting from loading applied to the beam with cross-sections symmetric about that plane. (c) The behavior of the beam under the assumption that there is a linear variation of displacement over the cross-section, the Euler-Bernoulli beam theory. This implies that the cross-sections which are normal to the centroid axis of the beam before deformation remain plane, undeformed and perpendicular to the deformed axis.This model accounts for bending moment effects on stresses and deformations. This implies that transverse shear forces are recovered from equilibrium but their effect on beam deformations is neglected. (d) Account for shear deformation by constraining the normal to remain straight, but not necessarily normal, to the centroidal axis after shearing and bending; Timoshenko beam theory. (e) Implement the finite element method to model the numerical analysis of the behavior of the displacement of curved beams. (f) To check how accurate the finite element we are going to use here, we compare the finite element approximation solution of the analytical solution, and this comparison is carried out with already existing finite element results in literature which have been proven to be very accurate in their approximation. The examples we are going to use are the typical examples of the shear and membrane locking. They are frequently dealt with to check the capability of analysis of the elements developed and to see if the locking phenomena is property overcome. Assumptions The main assumptions are that (a) The beam deformation is geometrically linear. Therefore the displacement gradients are small compared to unity1 , implying that the difference between the deformed and the undeformed beam is very small. 1

∂U k=: ε ¿ 1 ∂X where U is the displacement vector, and X is position of material point p in a reference configuration with ||G || =k

respect to a given origin O.

9

(b) The beam is prismatic with constant circular cross section and constant radius of curvature. (c) We are neglecting changes in dimensions of the cross sections as the beam is displaced.

10

Chapter 2

Straight Beam Theory To set ground for curved beams, we first give an overview of the deformation of a straight beam with a square cross section, while adopting definitions according to [3].

2.0.1 Geometry Y

q(x) Y

X

h Z

L

b Cross − Section

Figure 2.1: Geometry of a simple Beam (a) The beam is long and thin, implying that L >> b and h. (b) Loading is in the y− direction to obtain the bending moment. (c) Constant cross section. (d) The beam does not experience torsion.

2.0.2 Stresses (a) No load applied in the z- direction implying the stresses are zero in that direction. (b) The only significant stress are σxx and τxy . 11

2.1 Euler - Bernoulli Beam Theory We define v the axial deflection and u as the transverse deflection of midplane. After deformation, (a) Plane sections remain plane and perpendicular to the midplane. For this to be true, we require that (i). The beam should be bent only with bending moments, therefore there are no shear on transverse planes, as we will explain below; (ii). The applied loads are such that no twisting occurs. (b) In addition to this, we use the assumptions that the axis of the rod is inextensible. (c) The Euler-Bernoulli beam theory assumes that the internal energy of a beam is entirely due to bending strains and stresses, therefore the predominant stresses and strains are σxx and εxx . Therefore

εyy =

∂u =0 ∂y

εxy =

1 2

µ

∂u ∂v + ∂x ∂y

(2.1) ¶ =0

(2.2)

Thus from (2.1) u is not a function of y, thus u = u(x)

(2.3)

This implies that every point of the cross section at a given location along the x−axis has the same vertical displacement.(Timoshenko, Goodier [4]) Therefore

εxy = 0 ⇒

du ∂v ∂u =− =− ∂y ∂x dx

(2.4)

Thus v = −y

du = −yθ dx

(2.5)

implying that the axial displacement vary linearly with the y distance from the neutral axis. We also note that the equation (2.5) implies that the cross-sections’ rotational (bending) angle directly depend on the gradient of deformation, i.e

θ=

du dx

(2.6) 12

which along with (2.3) confirms the assumption (b) above that any cross-sectional plane stays a plane and normal to the deformed longitudinal axis, though it displaces and rotates a bit from its original position.

y,u

deformed state

θ

Normal to the reference beam Normal to deformed beam axis, which is also the direction of the deformed crosssections. du/dx

u

undeformed state x, v

Figure 2.2: Deformation of cross-section of an Euler - Bernoulli beam

2.1.1 Total potential energy For a given loaded elastic body under geometrical constrains or boundary conditions, the potential energy of the deformed body assumes a stationary value, and it attains an absolute minimum when the displacement of the body are those of equilibrium configuration. From (2.5), the axial strain distribution is

εxx =

∂v ∂ 2u d2u = −y 2 = −y 2 ∂x ∂x dx

(2.7)

Using the Hooke’s law, (Timoshenko,Goodier, [4]), the constitutive equation is thus

σxx = E εxx = −Ey

d2u dx2

(2.8)

the stress distribution in terms of displacement field, where E is the Young’s modulus of elasticity. The resulting bending moment M is defined as Z

d2u M = −yσxx dx = E 2 dx A

Z A

y2 dA = EI κ for κ =

Where I is the area moment of inertia

R

Ay

2 dA

d2u dx2

(2.9)

a property of the beam that is used to

predict its resistance to bending and deflection, and κ is the curvature of the beam. Since the internal energy for the Euler-Bernoulli beam accounts only for the bending 13

moment deformation, the total potential energy of the beam, subjected to force q acting in the y− direction is 1 2

Wi

=

We

= −

Z V

σxx εxx dV =

Z L 0

Z L 0

µ EI

d2u dx2

¶2 dx

(2.10)

qudx

1 = Wi +We = 2

⇒Π

1 2

Z L 0

µ EI

d2u dx2

¶2 −

Z L 0

qudx

Where Wi and We are internal and external virtual work respectively and L is the length of the beam.(Oden, [5])

2.2 Timoshenko Beam Theory As cited by Felippa [1] the Timoshenko beam corrects the classical beam theory with first-order shear deformation effects. Therefore the effect of the shear stresses on the deformation are taken into account. After deformation, (a) Plane sections are assumed to remain plane and rotate about the same neutral axis, with the assumption under (a) - (ii) in section 2.1 also holding for this case too. (b) Unlike the euler-Bernoulli beam, the cross sections do not remain perpendicular to the deformed longitudinal axis. Therefore

εxx =

∂v dθ = −y ∂x dx

1 εxy = 2

µ

∂v ∂u + ∂x ∂x

(2.11) ¶

µ ¶ 1 1 du = −θ + = ϒ 2 dx 2

(2.12)

Thus the deviation from the normal direction which is the difference between the normal to the longitudinal axis and the plane section rotation is the shear angle ϒ. It consists of the contributions from the gradient of the deflection (the rotation of the longitudinal axis), and the cross sectional rotation bending angle θ . This is further illustrated in figure (2.3). For Timoshenko beam the cross section deformation is the sum of two contributions that is one is due to bending, and the other is the shear deformation. 14

γ − Shear angle y, u

du/dx Normal to the reference beam θ

Normal to defromed beam axis du / dx

Direction of deformed cross section

u

x, v

Figure 2.3: Deformation of the cross-section of a Timoshenko beam

2.2.1 Total potential energy From equations (2.11) and (2.12), the axial and shear strains for a Timoshenko beam are

ε

= −y

dθ dx

ϒ = −θ +

(2.13) du dx

(2.14)

Therefore Z

Π

= =

Z

L 1 (σxx εxx + σxy εxy ) dV − qudx 2 V 0 ¶ ¶ µ µ Z Z Z L EI L d θ 2 κ GA L du 2 dx + dx − qudx −θ + 2 0 dx 2 0 dx 0

(2.15)

(2.16) where κ is the shear correction factor, G the shear modulus,A is the cross sectional area, and q is the external load applied to the beam, as has been explained by (Kwon, Bang [6]).

15

Chapter 3

Finite Element Method In this chapter we introduce the finite element numerical methods, that will be adopted for the problems to follow. Most numerical methods involve an approximation to an unknown function U by a new function U , which is a linear combination of the basis functions

n

Uˆ = ∑ Ui Wi i=1

where U are the unknown coefficients, Wi are the chosen weight functions. While classical methods employ global basis functions of various frequencies, approximation of this form will work adequately for simple geometries but yield poor approximations for more complex geometries, as global basis functions cannot reflect properties of the solution that are induced by the geometry. Instead of employing global basis functions, the finite element method approximates a solution for problems which are described by partial differential equations or formulated as functional minimization, based on the concept of discretization of the solution domain into sub-domains called finite elements. As summarized by (Nikishkov, [7]), this can be achieved with the following general steps. (a) Discretize the domain, that is to divide the solution region into finite elements. The solution domain is divided into several simpler finite elements, where each element has a simple geometry, so appropriate assumed solutions can easily be written for the element. (b) Selection of interpolation functions to describe solution over an element is the next step. 16

If we use a polynomial function to interpolate U within elements from nodal values of U , then for U to approach the exact values as the finite element mesh is repeatedly refined and if α is the highest derivative of U , then (i) Within each element, the assumed field for U must contain a complete polynomial of degree α or higher, (ii) Across the elements’ boundaries, there must be a continuity of order

α − 1 for U and its derivatives. (c) Establish the matrix equation for the finite element which relates the nodal values of the unknown function to other parameters.Different approaches used are for example the most the variational approach and the Galerkin method. (d) To find the global equation system for the whole solution, all element equations must be assembled. This involves the combination local element equations for all elements used for discretization. Element connectivities are used for the assembly process. Before solution, boundary conditions (which are not accounted in element equations) should be imposed. (e) Solving the global equation system. Since the finite element global equation system is typically sparse, symmetric and positive definite, direct and iterative methods can be used for solution. The nodal values of the sought function are produced as a result of the solution,since approximating functions are determined in terms of nodal values of as physical field which is sought.

3.0.2 Formulation of the Finite Element Method As we have mentioned in (c) above, several approaches can be used to transform the physical formulation of the problem to its finite element discrete state. If the physical formulation of the problem is described by a differential equation then the most popular method of its finite element formulation is the Galerkin method. If the physical problem can be formulated as minimization of a functional then variational formulation of the finite element equations is usually used.

3.1 Bubnov-Galerkin’s Approximation Method The formulation given here is standard, and can be found in many finite element resources. The construct of the finite element given below is in reference mainly from (Hughes, [8]); Also we refereed to approach used by (Nikishkov, [7]); ( 17

Strang et al., [9]); With a general example, let the problem to be solved in one dimensional be defined by the equation Lu

=

f in Ω

(3.1)

u|Γ

= uρ on Γ1

up|α p

= uβ on Γ2

1

Γ2

for L a given operator acting on functions u ∈ L p (Ω), such that the functions satisfy in some sense the boundary conditions. The fundamental question is then: To match such a space of functions u with a class of inhomogeneous terms f , in such a way that to each f , there corresponds one and only one solution u. Let the space

δ = {u | u ∈ H s , u|Γ = uρ }

(3.2)

1

be defined as a function space for trail solutions u. Define a second collection of functions know as test functions such that we require that they vanish at boundaries corresponding to the essential boundary condition uρ . That is V = {ω | ω ∈ H s , ω|Γ = 0}

(3.3)

L : H s → (H s )∗ , and f ∈ (H s )∗ the dual space of H s .

(3.4)

1

and

Then the weak formulation of the problem is to find u ∈ δ , such that for all

ω ∈ V, a(u, ω ) = ω T L u = (ω , f )

(3.5)

Where a : δ × V → R is a bilinear operator linear in its second argument. The

∂ upα p for any α with p α p ≤ s exist in the weak sense, and derivatives of small order needn’t exist, like in the classical sense. The Galerkin’s method with piecewise polynomial subspaces δ h and V h is the finite element method. Let the approximations of sets δ and V , denoted by δ h and V h respectively, be such that

δ h ⊂ in δ (i.e.; if u ∈ δ h , then uh ∈ δ )

(3.6a)

V h ⊂ in V (i.e.; if ω ∈ V h , then ωh ∈ V )

(3.6b)

18

Consequently from (3.1) uh|Γ = uξ

and

1

ω|hΓ = 0

(3.6c)

1

since this is an essential boundary condition then it must be satisfied by every function ω h Implying that the functions ω h , should be zero at that boundary. Given the collection V h , then to each member ν h ∈ V h , we construct a function uh ∈ δ h , by uh = ν h + uhρ

(3.7)

where uρh = uρ . Therefore δ h is all functions of the form (3.6a). The problem is now defined as: Given f , uρ , and uβ , find uh = ν h + uρh , where ν h ∈ V such that for all ω h ∈ V h ,

a(ω h , ν h ) = (ω h , f ) + ω|hΓ uβ − a(ω h , uρh ) 2

(3.8)

uρh and V h have to be explicitly defined. Note that the approximate solution is required to fulfil the essential boundary conditions, but not the natural boundary condition, which doesn’t prevent the approximations form converging in the H s norm, to a solution u, which does satisfy up|α p = uβ . Γ2

3.1.1 The Stiffness Matrix and the Load vector The Galerkin method leads to a coupled system of linear algebraic equations. Every test function ν h is determined by its nodal parameters, which are the unknowns ci of the discrete problem. Each of these nodal parameters is the value at a given node, of either the function or its derivatives. Let V h consist of all linear combinations of given functions which are denoted by φi : = Ω → R, for i = 1, 2, ..., n. Then if ω h ∈ V h , then there exists constants ci such that n

ω h (x) = ∑ ci φi (x)

(3.9)

i=1

for ci = wh (xi ) for ci = wh (xi ), the nodal values for wh (x), ( 1 if i = j φ j (xi ) = 0 if i 6= j i, j = 1, ..., n 19

and (φi )|Γ = 0, 1

i = 1, 2, ..., n

Thus V h is a linear space of dimensions n with the basis {φi }ni=1 . To define members of δ , we need to specify uξ . Defining another shape function as φn+1 : = Ω → R, with the property that / V h , then (φn+1 )|Γ = 1, φn+1 ∈ 1

uhρ

= uφn+1

and thus uρh |

= uρ

Γ1

For uh ∈ δ h , and constanst di , uh

i = 1, 2, ..., n

= ν h + uρ

(3.10)

n

∑ di φi + uρ φn+1

=

i=1

Substituting (3.9) and (3.10), into (3.8), then à ! à ! " a

n

n

i=1

j=1

∑ c i φi , ∑ d j φ j

n

=

∑ ci φi , f

i=1

n

+

∑ ci φi|Γ1

#

à uβ − a

i=1

n

!

∑ ci φi , uρ φn+1

i=1

(3.11) which simplifies to n

∑ ci Di = 0

(3.12)

i=1

For n

Di =

∑ a (φi , φ j ) d j − (φi , f ) − φi|Γ1 uβ + a (φi , φn+1 ) uρ

(3.13)

j=1

Since the solution hold ∀ ω h ∈ V h , and ∀ {ci }ni , i = 1, ..., n, arbitrary in (3.12), then n

∑ a (φi , φ j ) d j = (φi , f ) + φi|Γ1 uβ − a (φi , φn+1 ) uρ

(3.14)

j=1

for each Di , i = 1, .., n. Since everything is know besides d j , then (3.14) constitute a system of n equations, in n unknowns. Let Ki j = a (φi , φ j )

(3.15a) 20

Fi = (φi , f ) + φi| uβ − a (φi , φn+1 ) uρ

(3.15b)

Γ1

Then (3.14) while adopting a matrix notation can compactly be restated as solving for d in the algebraic equations defined as Kd = F

(3.16)

3.2 Assembly of global stiffness matrix and force vector The elements stiffness must now be assembled into the beam’s stiffness, implying now that to add values of element stiffnesses in global axes, we need to adjust the numbering into a numbering system of the beam as a whole. This particular task is accomplished by the location matrix.The location matrix (if we denote it by loc) dimensions are the number of element nodes by the number of elements. Thus give a particular degree of freedom, and an element number, for example i and j respectively, then the value returned by the location matrix array is the corresponding global equation number A˜ such that   j       j+1   j+2 A˜ = loc(i, j) =   .  ..       j+n −1 el

if i = 1 if i = 2 if i = 3

for j = 1, 2, ..., nel

if i = nel

Therefore the contribution of an element to the beams’ global stiffness is obtained by adding to the beams’ global stiffness matrix term in row loc(i), and (e)

column loc( j), the term Kloc(i) , Kloc( j) , the element stiffness term Ki, j . This applies to the force vector too, where the representation of the global force vector is now (e)

Floc(i) = Floc(i) + Floc(i)

(3.17)

In the next two subsections, we are going to discuss a plane two-dimensional, linear Euler-Bernoulli beam finite element, and Timoshenko beam element. From the definitions of the potential energy for each beam,we derive the displacement formulation on the basis of the galerkin method already described.

21

Chapter 4

Linear Curved Beams In our discussion in the previous chapters, the deformation of the beam subjected to a load, has been restricted to beams with the assumption that the longitudinal elements have the same length. The xy- longitudinal plane has been assumed to be the plane of symmetry with the load applied in this plane. This has restricted the beam theory formulated above to initially straight beams of constant cross section. Although considerable deviations from this restriction can be tolerated in real problems, when the initial curvature of the beams becomes significant, the linear variations of strain over the cross section is no longer valid, even though the assumption of plane cross sections remaining plane is valid. The cross section of part of an initially curved beam is described by the figure (4.1), with xy plane the plane of symmetry. The radii R references the location of the centroid of the cross sectional area; Rn references the location of the neutral axis; and r references some arbitrary point p of area element dA on the cross section. For curved beams, the longitudinal deformation of any element will be proportional to the distance of the element from the neutral surface, implying that the total deformations are proportional to the distance from the neutral axis. This implies then that the linear variation of the strain over the cross section is no longer valid as the beam’s elements are not of equal length, even though the assumptions of the plane cross section remaining plane after deformation is valid, as we assumed for the straight beam case. This citation can be found at university of Washington’s department of mechanical engineering website [13]. Since the stress during deformation is proportional to the strain of the beam, then

22

Stress distribution Centroidal Axis

A’ A’ Rn − r M p

Rn Neutral Axis

A

R

p

r

A

∆θ

∆θ

R = Radius of curvature Rn = Radius of Neutral axis r = Radius of general fiber in the beam M = Bending moments ∆θ = Cross sectional rotations

Figure 4.1: Bending stress in curved beams the elastic stress of a curved beam is not promotional to the distance of the element from the neutral axis. For the same reason, the neutral axis in a curved beam does not pass through the centroid section. We define the geometry of the curved beam with reference to (Oden, [5]). Consider a beam, having a constant cross section and a constant radius of curvature in the plane of bending. The geometry of the curved beam is defined by establishing two orthogonal coordinates systems: one fixed coordinate system (x, y, z), with its origin located on a convenient cross-section, and another curvilinear system (s, y, z) in which s is tangent to the curve, y is a radial coordinate directed toward the center of the geometric axis, and z is directed normal to the plane of the beam. The assumptions taken are that, (a) The beam material is homogeneous, isotropic and linearly elastic. (b) The beam axis always lies in the plane of bending. (c) The Radius of curvature is constant and the cross section is constant and has symmetry about the plane of curvature, that is xy-plane. 23

(d) External loads act in the plane of curvature, therefore the plane of bending and the plane of curvature are the same. (e) Plane sections remain plane, and the cross section does not deform. We note that these assumptions are almost the same as the ones used for straight beams only that the beam is initially curved in this case.

∆s

s

My Vy y R Ms

V y + ∆ Vy

M s + ∆ Ms

Ns

Ns + ∆ M s My +∆ M y

z

Mz + ∆ Mz Vz + ∆ V z s

Vz Mz

Figure 4.2: Forces and displacements on an infinitesimal curved beam element

4.1 The Equilibrium Equations With reference to figure (4.2), the equilibrium equations at a point p in the deformed beam element have been derived and found to be dNs Vy − + q1 ds R dVy Ns + + q2 ds R dVz + q3 ds

= 0

(4.1)

= 0 = 0

For the equilibrium of the moments we get dMs My − + m1 ds R dMy Ms + −Vz + m2 ds R dMz +Vy + m3 ds

= 0 = 0 = 0 24

(4.2)

Where Ns ,Vy ,Vz are Normal force along s, and shear forces in in y and z direction, while My , Mz are bending moments around y and z axis, and Ms is the twisting moment around the longitudinal axis, and {qi , mi ,

i = 1, ..., 3} are external

loads and moments applied in the tangential, radial direction and z direction respectively. Detailed derivation are derived by (Oden, [5]), (Blake [14], pg. 575 ) and ( Lebeck, Knowlton, [15]) derives the equilibrium equations in polar coordinates.

4.2 Normal stresses in Curved Beams We follow the derivation given by (Oden, [5], pg. 64). Following similar derivation is (Wallerstein, [16], pg. 99). To be able to find the displacement of the beam, we find the stress-strain relations,assuming that the material is elastic and obeys the Hooke’s law. First the derivation of the pure bending of the curved beam is presented, implying that we use the assumption that plane sections normal to the beam’s axis before deformation remain plane and normal to this axis after deformation. With this assumption,the components of displacement in a direction normal to a cross section must satisfy the equation of a plane u = cˆ0 + cˆ1 y + cˆ2 z

(4.3)

where cˆ0 , cˆ1 , cˆ2 are constants. Let the geometry of element between the cross-sections of the rod be as in the diagram below.

∆s ∆s−y

y R

Figure 4.3: Geometry of a beam’s element vertical displacement of magnitude y Materials located at a radial distance y from the axis are of different length ∆sy . Thus taking limit as ∆s approaches zeros, 1 ds = dsy 1 − y/R 25

From strain-displacement formula,

∂v ∂ sy

εs =

(4.4)

Substituting this to (4.3), then

εs

ds ds ds + c1 y + c2 z dsy dsy dsy ds (c0 + c1 y + c2 z) dsy

= c0 =

for d cˆ0 d cˆ1 d cˆ2 , c1 = , c2 = ds ds ds

c0 =

Implying that 1 (c0 + c1 y + c2 z) 1 − Ry

εs =

(4.5)

Thus the normal stress is

σs = E εs =

E (c0 + c1 y + c2 z) 1 − Ry

(4.6)

To determine c0 , c1 , c2 then: Z

Ns

=

Mz

=

Z

1 σs dA = c0 E y dA + c1 E 1 − A A R

Z

Z

A

σs ydA = c0 E

Z

My

=

Z

Z

A

σs zdA = c0 E

y y dA + c2 E 1 − A R

y y dA + c1 E A 1− R

z y dA + c1 E A 1− R

Z

Z

Z

z dA y (4.7) 1 − A R

y2 y dA + c2 E A 1− R

yz y dA + c2 E A 1− R

Z

Z

yz y dA A 1− R

z2 y dA A 1− R

The inertia terms, also derived by (A.O.Lebeck, J.S.Knowlton, [15]) can be defined as Z

Iy Iyz

z2 y dA A 1− R Z yz = y dA A 1− R =

(4.8)

Z

Iz

=

y2 y dA A 1− R

We note here when the dimensions of the cross section are small compared to the radius of curvature of the longitudinal axis, that is y > 1. Further, there is deterioration of radial deflection convergence, as the values are not converging towards unity but away from unity. Since we have used reduced integration, there is generally an improvement of convergence, as seen in figures (4.18), and (4.20) but as from the results in table

57

(4.3), even with reduced integration for the choice of the interpolation functions used, the reduced integration approach is not trouble free.

4.10.5 Results for the Euler-Bernoulli curved beam If we again compare the analytical solution with the numerical solution for the example give in figure(4.21) for an Euler-Bernoulli beam, the analytical solution is without the shearing term, that is vc

=

uc

=

duc /dx

=

QR3 QR − 2EI 2EA π QR3 π QR + 4EI 4EA 2 QR EI

(4.73a) (4.73b) (4.73c)

Comparing to the results to an already existing model by Saffari and Tabatabeai [2] , we get the results in (4.4) for the displacement of the quarter circular EulerBernoulli beam fixed at one end and subjected to a radial point force at the other free end. We note that first of all the tangential displacement in this case is not as R/b

1 Element - Saffari and Tabatabaei

1 Element - Present Element

v f /vc

u f /uc

θ f /θc

v f /vc

u f /uc

(du f /ds)/(duc /ds)

4

1.0133

1.0000

1.0000

1.5040

0.6906

0.5284

10

1.0053

1.0000

1.0000

1.5367

0.7406

0.5682

20

1.0013

1.0000

1.0000

1.5451

0.7582

0.5791

50

1.0002

1.0000

1.0000

1.5497

0.7690

0.5851

100

1.0000

1.0000

1.0000

1.5511

0.7726

0.5870

200

1.0000

1.0000

1.0000

1.5519

0.7745

0.5879

500

1.0000

1.0000

1.0000

1.5523

0.7756

0.5885

1000

1.0000

1.0000

1.0000

1.5524

0.7759

0.5887

Table 4.4: Comparison of the two different finite element solution models loaded with tip radial force, of a quarter circular cantilever beam with Castigliano’s energy solutions for Euler-Bernoulli Beam. accurate as for the Timoshenko beam, for Euler-bernoulli beam subject to bending moments only. The beam cannot represent the zero radial shear strains assumption by our use of low order interpolation functions, hence we used reduced integration, but still the results are not very accurate. Considerable improvement in radial approximation is noted, where we used cubic hermitian polynomials to 58

approximate the displacements u, while the cross sectional rotation approximation is not better than for the Timoshenko beam.

4.11 Points to Note: Timoshenko curved beam strains definition We have explained how the equations governing the deformation of a curved Euler-Bernoulli beam are derived, based on the derivation given by (Oden, [5]). To consider the shearing effect, hence the Timoshenko beam term, we have used the stress tensor σsy defined from the same source to formulate the differential equation governing the shear deformation. These strains are the conventional strains definitions used in the deformation of a curved shearable beam, for example (Prathap and Ramesh, [26] ), and (Zienkiewicz et al. [30]) use this definition of strain. Based on minimum potential energy of a deformed beam, we have calculated the displacement of a loaded beam under boundary constrains. While the energy strains of the curved Euler-Bernoulli beam are correct, the representation of the Timoshenko curved beam energy strains, based on the definition we have used here are not according to studies by (Day, Potts, [31]) entirely correct. The strain definition used here is

dv u − ds R

ε

=

χ

= −

γ

= −θ +

(4.74a)

d2u dv − ds2 Rds

(4.74b)

du ds

(4.74c)

We have used rotation due to bending θ as a degree of freedom, where the curvature changes due to bending is defined therefore as

χ

= −

dθ dv − ds Rds

(4.75)

Day and Potts, [31], while considering the rigid-body displacements, have showed that θ = du/ds is not equal to the physical rotation of the cross section for a 59

shearable beam and hence is not a suitable parameter to be used as a nodal degree of freedom, whose continuity is to be imposed. We can explain this in other way as: we know that for a straight beam, the rotation is a function depending only on the vertical displacement u, but when the beam is initially curved, the tangential displacement will also now affect the cross sectional rotations, therefore the total rotation θˆ at a point becomes du/ds + v/R. Therefore the appropriate choice of nodal degree of freedom, instead of θ , would be θ + v/R, hence we can define a single variable for example θˆ = θ + u/R so that now the strains for the curved timoshenko beam will instead be defined as

ε

=

dv u − ds R

(4.76a)

χ

=

d θˆ ds

(4.76b)

γ

du v = −θˆ + + ds R

(4.76c)

This definition is used in many problems relating to shearable curved beams, and is not a new concept. If we now compare this definition of strains with the conventional strains we have previously used for the example two, we get the following results in table (4.5) . R/b

1 Element - Days and Pott

1 Element - Present Element

v f /vc

u f /uc

θ f /θc

v f /vc

u f /uc

θ f /θ c

4

1.2386

0.6110

1.1985

1.2386

0.6110

0.7117

10

1.1914

0.5945

1.1985

1.1914

0.5945

0.7306

20

1.1775

0.5885

1.1985

1.1775

0.5885

0.7361

50

1.1696

0.5847

1.1985

1.1696

0.5847

0.7393

100

1.1670

0.5835

1.1985

1.1670

0.5835

0.7403

200

1.1657

0.5828

1.1985

1.1657

0.5828

0.7408

500

1.1649

0.5825

1.1985

1.1649

0.5825

0.7411

1000

1.1647

0.5823

1.1985

1.1647

0.5823

0.7412

Table 4.5: Comparison of the two different strain definition for fem solution models loaded with tip radial force, of a quarter circular cantilever beam with Castigliano’s energy solutions for Timoshenko Beam.

60

We can note that there is a difference in values approximated for the displacements when using the two definitions of strains especially in cross sectional rotations. We still do not achieve convergence to the exact value of unity, even with this definition of strain, as we still use low-order linear basis to interpolate the displacements. When we use the conventional definition of strains for a shearable curved beam, then we use an inappropriate parameter θ = θˆ − u/R as a nodal degree of freedom. According to (Day, Potts, [31]), this parameter when chosen as the rotation degree of freedom is dependent on the geometry and orientation of the elements, while the rotation, θˆ , is uniquely defined at all points. Therefore for normalized displacements approximation for strains definition in (4.76), with the conventional strain definition (4.74) for the quarter-circular cantilever beam in figure (4.21), we get the following error in figure (4.23)in our approximation Strain displacement errors for a Timoshenko curved beam 0.5 0.45

cross sectional rotation error

Maximum Absolute Error

0.4 0.35 0.3 0.25 0.2 0.15 0.1

Tangential Displacement 8 Error x 10

Radial Displacement 8 Error x 10

0.05 0 0.5

1 1.5 2 2.5 3 3.5 Tangential, Radial displacements and Cross Sectional rotation errors

Figure 4.23: Maximum Absolute Error for Tangential, Radial displacements and Cross sectional Rotation when the conventional strains of a Timoshenko curved beam are used The rotation due to bending, as a degree of freedom for a timoshenko beam results into incorrect results, and this is magnified in the cross sectional rotations, as the tangential and radial displacements are approximated with a very small and negligible error, while the approximations for cross sectional rotations have a large error.

61

4.12 Summary and Conclusion We have investigated the definition and consequently the difference between the deformation of a beam subjected to external loading for an Euler-Bernoulli beam, and a beam under Timoshenko beam theory in chapter 2. Introduction to finite elements is covered in chapter three, while the curved beam is introduced in chapter4. In order to investigate how a ring deformers when subjected to external loads and moments, and from the definitions given herein, cited from references quoted, we have noted that the deformation of a curved beam, unlike the straight beam is subject to nonlinear variation of strain distribution even though the cross section of the beam remains plane after deformation. Further more, the equations governing the deformation of the curved beam are coupled. This makes the displacement of a curved beam to be a bit more difficult to tackle than the straight beam, but there is an advantage to that, as curved beams are efficient in load transfer due to inclusion of bending, shearing and membrane action in the transfer. We have implemented the equations of the curved beam to study their displacement when subjected to external loads in finite element method, and have assumed linear basis functions, for Timoshenko beam, while for Euler-Bernoulli beam, have assumed linear basis functions for tangential displacements, and cubic hermitian for radial displacement. To remedy shear and membrane locking, we have used 1− point gaussian integration for the shear strain and membrane strain. The numerical results we have gotten, we have compared with already existing different beam elements that gives very accurate approximation, and for the second example that is used is compared with the analytical solution, For the choice of the basis functions used, the results are not an accurate approximation of the exact solution. From our results, (a) Although there is notable improvement in the convergence when we used reduced integration, (4.18) the finite element solution does not converge to the analytical solution for our test examples with only use of one element, instead it converges to a value lower that the expected. We choose to use one element because the elements we are using for comparison, already existing in literature have been proven to be very accurate with just the use of one element, hence we use this as a comparison bench for the curved

62

beam element we have used here. (b) The use of reduce integration comes along with its disadvantages. When used, reduced integration admits deformation modes which are able to occur without any changes in strain energies (hence the name zero - energy modes).This results because the nodal displacements when using reduced integration are now disassociated from the gauss points, therefore there is no stiffness associated with this mechanism and these can lead to considerably inaccurate solutions. (c) Errors especially in rotation of the cross section arises from the use of the incorrect definitions of a Timoshenko beam strains, and we have with reference to [31] given the proper definition of strains. (d) Therefore though they are the most elementary and simple finite elements to implement, the finite elements implemented here for a curved beam under both Timoshenko and Euler-Bernoulli beam theory are not the most appropriate elements to use when implementing in practical problems, as they do not give very accurate results, and are prone to severe shear and membrane locking.

63

Appendix .1 Matlab Programmes .1.1 Common programmes and Parameters for both beam theories %------------------------------------------------------------------clear; clear figure clc nel

= 30;

%Number of elements

nnel

= 2;

%Number of nodes per element

ndof

= 3;

%Number of degrees of freedom per node

nnode = (nnel-1)*nel+1; %Total number of nodes edof

= nnel*ndof;

% Total number of degrees of freedom per element

sdof

= nnode*ndof; % Total number of degrees of freedom in all elements

%-------------parameters-------------------------------------------E =

27.6*10^6;

%Young’s Modulus of Elasticity (Stainless steel)

G = 3.8*10^6;

% Shear Modulus of Elasticity (Stainless steel)

b = 1;

% Breadth of the beam

A = pi*(b/2)^2

%Cross sectional Area (Circular cross section)

I = pi*(b/2)^4/4;%Area moment of inertia of cross section tleng = 40;

%Beam span length

h = tleng/nel;

%Mesh Size

R = 15;

%Radius

theta = pi;

%Subtended Angle of the beam

cf = 5/6;

%Shear correction factor

Q

%External load

= 500;

u_m=Q; v_m=Q; w_m=Q; %External Tangential, Radial and Moment respectively. %----ESSENTIAL--BOUNDARY--CONDITIONS-------------------------------------%Example boundary for a clamped beam on both ends, which means that %the first three and the last three dofs values is zero. bcdof = [1,2,3,nel*3+1,nel*3+2,nel*3+3]; bcval = [0, 0, 0,0,0,0]; 64

%-------------LOCATION MATRIX--------------------------------------------function[locmatrix] = eldof(iel,nnel,ndof,nel) %compute dofs associated with each element %locmatrix = matrix system of dof vector associated with element iel %iel

= element number whose global dofs are to be determined

edof = nnel* ndof; count = (iel -1)* (nnel -1)*ndof; for i = 1:edof locmatrix(i) = count + i; %gives location of an element at global level end %----------APPLY

ESSENTIAL BCS TO THE MATRIX EQN---Kd = F --------------

function[kk,ff] = febc(kk,ff,bcdof,bcval,nel) n

= length(bcdof);

sdof

= size(kk);

for i = 1:n c

= bcdof(i); % gives the dofs where bcs are to be applied.

for j = 1:sdof kk(c,j) = 0; %NOTE end kk(c,c) = 1; ff(c)

= bcval(i);

end NOTE: Removing the equations corresponding to zero boundary conditions is very convenient for small matrices,because it reduces down on the number of equations to solve. But when we were implementing in Matlab the equations have to be rearranged. With rearrangement it was taking a bit more time than solving the equations. To apply boundary conditions without rearranging the equations we set to zero rows and columns corresponding to prescribed zero displacements as well as the in the force vector too. We then place ones on the diagonal to maintain non-singularity. %----------GAUSS INTEGRATION --POINTS---(Adopted from Kwon,Bang [6],pp 178;)------function [point1,weight1]= gaussint(ngl,h) %

ngl

- number of integration points

65

%

point1

- vector containing integration points

%

weight1 - vector containing weighting coefficients

%------------------------------------------------------------------point1=zeros(ngl,1); weight1=zeros(ngl,1); %find corresponding integration points and weights if ngl==1

% 1-point quadrature rule

point1(1)=0.0; weight1(1)=2.0; elseif ngl==2

% 2-point quadrature rule

point1(1)=-0.577350269189626; point1(2)=-point1(1); weight1(1)=1.0; weight1(2)=weight1(1); elseif ngl==3

% 3-point quadrature rule

point1(1)=-0.774596669241483; point1(2)=0.0; point1(3)=-point1(1); weight1(1)=0.555555555555556; weight1(2)=0.888888888888889; weight1(3)=weight1(1); elseif ngl==4

% 4-point quadrature rule

point1(1)=-0.861136311594053; point1(2)=-0.339981043584856; point1(3)=-point1(2); point1(4)=-point1(1); weight1(1)=0.347854845137454; weight1(2)=0.652145154862546; weight1(3)=weight1(2); weight1(4)=weight1(1); end

.1.2 Curved Beam - Euler-Bernoulli beam theory %---------INITIALIZE EMPTY STIFFNESS, LOAD and LOCATION MATRICES---------ff =zeros(sdof,1); kk = zeros(sdof,sdof);

66

locmatrix = zeros(nnel*ndof,1); %set up location matrix, to renumber the % dofs when doing assembly into global % system. for iel = 1:nel locmatrix

%Loop over the elements

= eldof(iel,nnel,ndof);

%Element stiffness matrix------------------------------kc = eulerelemetmatrix(E,I,G,A,cf,R,h,theta,ndof); %Element load vector----------------------------------fs = eulerelementvector(u_m,v_m,h,w_m,R,theta,ndof); %Global load vector-----------------------------------ff = eulerglobalvector(ff,fs,locmatrix,nel); %Global stiffness matrix-----------------------------kk = eulerglobalmatrix(kk,kc,locmatrix); end [kk,ff] = febc(kk,ff,bcdof,bcval); %Implement Essential Boundary conditions sol = kk\ff; % Solve for displacement %-------------DETERMINATION OF ELEMENT STIFFNESS MATRIX ---------------function[kc] = eulerelemetmatrix(E,I,G,A,cf,R,h,theta,ndof) %-------------Rotation matrix------------------------------------------rot = [cos(theta) sin(theta) 0 0 0 0;... -sin(theta) cos(theta) 0 0 0 0;... 0 0 1 0 0 0;... 0 0 0 cos(theta) sin(theta) 0;... 0 0 0 -sin(theta) cos(theta) 0;... 0 0 0 0 0 1]; [axial ]

= euleraxial(R,h,ndof); %Axial Strain Stiffness Matrix

%Contribution [bend,Bc ] = eulerbend(R,h,ndof); %Bending Strain Stiffness Matrix %Contribution ktt = (E*A/2)*axial + (E*I/2)*bend; kc = rot’*ktt*rot; %--------------AXIAL- STRAIN--MATRIX----------------------------function[axial ] = euleraxial(R,h,ndof) ngl =

1;

% 1- point integration

[point1,weight1] = gaussint(ngl,h);

% extract integration points as weights

axial = zeros(ndof*2,ndof*2); for int=1:ngl x

= point1(int);

67

wt = weight1(int); %---------Linear---Interpolation functions---------------------s_1

= 1/2*(1 - x); s_2

= 1/2*(1 + x);

st_1 = -1/2; st_2 = 1/2; % 1st derivative of the shape functions %---------------cubic--Hermitian Interpolation functions--------c_1 = 1/4*(2- 3*x + x^3); c_2 = 1/4*(1 - x - x^2 + x^3); c_3=1/4*(2+3*x - x^3); c_4 = 1/4*(-1-x + x^2 + x^3); ------First derivative---------------------------------------c_u = [(-3/4+3/4*x^2)*2/h ,-1/4-1/2*x+3/4*x^2, (3/4-3/4*x^2)*2/h, -1/4+1/2*x+3/4*x^2]; ------Second derivative-------------------------------------c_uu = [(3/2*x)*4/h^2, (-1/2+3/2*x)*(2/h), (-3/2*x )*4/h^2, 1/2+3/2*x*(2/h)]; c_u1 = 6*x/h^2; c_u2 = (-1 + 3*x)/h; c_u3 = -6*x/h^2; c_u4 = (1 + 3*x)/h; %-------Axial strain Matrix-----------------------------------------a_ns = [ 2*st_1/h, -c_1/R, -c_2/R, 2*st_2/h, -c_3/R, -c_4/R ]; A = a_ns’* a_ns; for i = 1: ndof*2 for j = 1:ndof*2 axial(i,j) = axial(i,j) + A(i,j)*wt*h/2; end end end %--------------BENDING- STRAIN--MATRIX--------------------------function[bend,Bc ] = eulerbend(R,h,ndof) ngl = 4 ;

% 4 - point integration

[point1,weight1]=gaussint(ngl,h); bend = zeros(ndof*2,ndof*2); Bc = zeros(ndof*2,ndof*2); for int=1:ngl x=point1(int); wt=weight1(int); %---------Linear---Interpolation functions---------------------s_1

= 1/2*(1 - x);

s_2

= 1/2*(1 + x);

---1st derivative of the shape functions----------------

68

st_1 = -1/2; st_2 = 1/2; %---------------cubic--Hermitian Interpolation functions--------c_1 = 1/4*(2- 3*x + x^3); c_2 = 1/4*(1 - x - x^2 + x^3); c_3=1/4*(2+3*x - x^3); c_4 = 1/4*(-1-x + x^2 + x^3); ------First derivative---------------------------------------c_u = [(-3/4+3/4*x^2)*2/h ,-1/4-1/2*x+3/4*x^2, (3/4-3/4*x^2)*2/h, -1/4+1/2*x+3/4*x^2]; ------Second derivative-------------------------------------c_uu = [(3/2*x)*4/h^2, (-1/2+3/2*x)*(2/h), (-3/2*x )*4/h^2, 1/2+3/2*x*(2/h)]; c_u1 = 6*x/h^2; c_u2 = (-1 + 3*x)/h; c_u3 = -6*x/h^2; c_u4 = (1 + 3*x)/h;; %-------Bending strain matrix ---------------------------------b_ns = [-st_1*2/(h*R), -c_u1, -c_u2, -st_2*2/(h*R), -c_u3, -c_u4]; B =

b_ns’* b_ns;

for i = 1: ndof*2 for j = 1:ndof*2 bend(i,j) = bend(i,j) + B(i,j)*h/2*wt; Bc(i,j) = Bc(i,j) + P_e(i,j)*h/2*wt; % Stiffness Matrix %contribution for Robin type boundary condition, end end end %--------------GLOBAL ELEMENT MATRIX-------------------------------function[kk] = eulerglobalmatrix(kk,kc,locmatrix) N = length(locmatrix); % # of dofs per element for i = 1:N % gives element row position ii = locmatrix(i) ;%global row for j = 1:N % element column position jj = locmatrix(j); % global column kk(ii,jj) = kk(ii,jj)+kc(i,j); %assemble global matrix end end %-------ELEMENT---LOAD----VECTOR-------------------------------------

69

function[fs] = eulerelementvector(u_m,v_m,h,w_m,R,theta,proc,ndof) rot = [cos(theta) sin(theta) 0 0 0 0;... -sin(theta) cos(theta) 0 0 0 0;... 0 0 1 0 0 0;... 0 0 0 cos(theta) sin(theta) 0;... 0 0 0 -sin(theta) cos(theta) 0;... 0 0 0 0 0 1]; Ax = eulergaussvtwo(R,h,ndof); Bn = eulergaussvfour(R,h,ndof); ftt = [u_m*Ax + v_m*Bn]; fs

= rot*ftt;

%---------------TANGENTIAL DISPLACEMENT FUNCTION CONTRIBUTION---------function[Ax] = eulergaussvtwo(R,h,ndof) ngl = 2; %two - point integration [point1,weight1]= gaussint(ngl,h); Ax = zeros(ndof*2,1); for int = 1:ngl x

= point1(int);

wt = weight1(int); %---------Liner shape functions-------------------------------------------s_1

= 1/2*(1 - x);

s_2

= 1/2*(1 + x);

st_1 = -1/2; st_2 = 1/2; %--------Load vector calculation-----------------------------------------a = [s_1,

0,

0 , s_2, 0, 0]’;

for i = 1: ndof*2 Ax(i) = Ax(i) + a(i)*wt*h/2; end end %---------------RADIAL DISPLACEMENT FUNCTION CONTRIBUTION----------function[Bn] = eulergaussvfour(R,h,ndof) ngl = 4; [point1,weight1]= gaussint(ngl,h); Bn = zeros(ndof*2,1); for int = 1:ngl x

= point1(int);

70

wt = weight1(int); %-----------cubic hermitian interpolation functions-----------------c_1 = 1/4*(2-3*x+x^3); c_2 = 1/4*(1-x-x^2+x^3); c_3 = 1/4*(2 +3*x-x^3); c_4 = 1/4*(-1-x + x^2 + x^3); %------------Load vector calculation----------------------------------b = [0, c_1, c_2, 0, c_3, c_4 ]’; for i = 1: ndof*2 Bn(i) = Bn(i) + b(i)*wt*h/2; end end %--------------GLOBAL ELEMENT MATRIX---------------------------------function[ff] = eulerglobalvector(ff,fs,locmatrix,nel) N = length(locmatrix); for i = 1:N

%

j = locmatrix(i); ff(j)= ff(j)+ fs(i);

gives local vector row position % global position %assembled matrix

end

.1.3 Curved Beam - Timoshenko Beam Theory %---------INITIALIZE EMPTY STIFFNESS, LOAD and LOCATION MATRICES--------------ff = zeros(sdof,1); kk = zeros(sdof,sdof); locmatrix = zeros(nnel*ndof,1); %set up location matrix, to renumber the % dofs when doing assembly into global % system. for iel = 1:nel locmatrix

%Loop over the elements

= eldof(iel,nnel,ndof);

%Element stiffness matrix---------------------------Kel = shearelementmatrix(E,I,A,G,cf,R,h,theta,ndof,tleng); %Element load vector----------------------------------Fel = shearelementvector(u_m,v_m,w_m,h,R,theta,ndof,tleng,cf,E,A,I,G); %Global load matrix---------------------------------F

= shearglobalvector(F,Fel,locmatrix,nel);

%Global stiffness matrix------------------------------

71

KK

= shearglobalmatrix(KK,Kel,locmatrix);

end [kk,ff] = febc(kk,ff,bcdof,bcval); %Implement Essential Boundary conditions sol = kk\ff; % Solve for displacement %-------------DETERMINATION OF ELEMENT STIFFNESS MATRIX ---------------function[Kel]=shearelementmatrix(E,I,A,G,cf,R,h,theta,ndof,tleng) %-------------Rotation matrix------------------------------------------rot = [cos(theta) sin(theta) 0 0 0 0;... -sin(theta) cos(theta) 0 0 0 0;... 0 0 1 0 0 0;... 0 0 0 cos(theta) sin(theta) 0;... 0 0 0 -sin(theta) cos(theta) 0;... 0 0 0 0 0 1]; %----Energy strains --stiffness matrix contribution --------------------[axial, shear] = axialnshear(R,h,ndof); %Axial Strain Stiffness Matrix %Contribution [bend,Bc]

= bendenergy(R,h,ndof); %Bending Strain Stiffness Matrix

%Contribution kel = (E*A/2)*axial + (E*I/2)*bend + (cf*G*A/2)*shear; kc

= rot’*ktt*rot;

Kel = rot’*kel*rot; %--------------AXIAL and SHEAR STRAIN MATRICES---------------------function[axial, shear] = axialnshear(R,h,ndof) ngl=

1;

% 1- point integration

[point1,weight1]=gaussint(ngl,h); axial = zeros(ndof*2,ndof*2); shear = zeros(ndof*2,ndof*2); for int=1:ngl x=point1(int); wt=weight1(int); %---------Linear---Interpolation functions---------------------s_1

= 1/2*(1 - x);

s_2

= 1/2*(1 + x);

-------1st derivative of the shape functions------------------st_1 = -1/2; st_2 = 1/2;

72

%-------Strain Matrix-----------------------------------------a_x = [2/h*st_1,

-s_1/R,

0 ,

2/h*st_2,

-s_2/R,

0];

s_h = [0, 2/h*st_1, -s_1, 0, 2/h*st_2, -s_2 ]; A_e =

a_x’* a_x;

S_e =

s_h’* s_h; for i = 1: ndof*2

for i = 1:

ndof*2

for j = 1: ndof*2 axial(i,j) = axial(i,j) + A_e(i,j)*h/2*wt; shear(i,j) = shear(i,j) + S_e(i,j)*h/2*wt; end end end %--------------BENDING STRAIN MATRIX--------------------------function[bend,Bc ] = ngl = 2 ;

bendenergy(R,h,ndof)

% 2 - point integration

[point1,weight1]=gaussint(ngl,h); bend = zeros(ndof*2,ndof*2); Bc

= zeros(ndof*2,ndof*2);

for int=1:ngl x=point1(int); wt=weight1(int); %---------Linear---Interpolation functions---------------------s_1

= 1/2*(1 - x);

s_2

= 1/2*(1 + x);

% -------1st derivative of the shape functions----------------st_1 = -1/2; st_2 = 1/2; %----strain matrix ---------------------------------------a_x = [2/h*st_1,

-s_1/R,

0 ,

2/h*st_2,

-s_2/R,

0];

s_h = [0, 2/h*st_1, -s_1, 0, 2/h*st_2, -s_2 ]; p = [N_1,0,0,N_2,0,0]; % Boundary contribution to the stiffness matrix A_e =

a_x’* a_x;

S_e =

s_h’* s_h;

P_e = 1/R*p’*p; for i = 1: ndof*2 for j = 1:ndof*2 bend(i,j) = bend(i,j) + B_e(i,j)*h/2*wt;

73

Bc(i,j) = Bc(i,j) + P_e(i,j)*h/2*wt; % Stiffness Matrix %contribution for Robin type natural boundary condition, %example when the moment is zero, then du/dx - v = 0; end end end %--------------GLOBAL ELEMENT MATRIX-------------------------------function[kk] = eulerglobalmatrix(kk,kc,locmatrix) N = length(locmatrix); for i = 1:N % gives element row position ii = locmatrix(i) ;%global row for j = 1:N % element column position jj = locmatrix(j); % global column kk(ii,jj) = kk(ii,jj)+kc(i,j); %assemble global matrix end end %-------ELEMENT---LOAD----VECTOR------------------------------------function[Fel] = shearelementvector(u_m,v_m,w_m,h,R,theta,meth,ndof,tleng,cf,E,A,I,G) rot = [cos(theta) sin(theta) 0 0 0 0;... -sin(theta) cos(theta) 0 0 0 0;... 0 0 1 0 0 0;... 0 0 0 cos(theta) sin(theta) 0;... 0 0 0 -sin(theta) cos(theta) 0;... 0 0 0 0 0 1]; [Ax,Bn,Mn] = sheargaussvector(R,h,ndof); fst

= [u_m*Ax + v_m* Bn + w_m*Mn];

Fel

= rot*fst;

%-------Gauss--Numerical integration --------------------------------function[Ax,Bn,Mn] = sheargaussvector(R,h,ndof) ngl = 2; [point1,weight1]= gaussint(ngl,h); Ax = zeros(ndof*2,1); Bn = zeros(ndof*2,1); Mn = zeros(ndof*2,1); for int=1:ngl x=point1(int); wt=weight1(int);

74

% extract integration points as weights

%-------Linear Basis Functions----------------------------------------s_1

= 1/2*(1 - x);

s_2

= 1/2*(1 + x);

st_1 = -1/2; st_2 = 1/2; %------------Load vector calculation----------------------------------ax = [s_1,

0,

0 , s_2, 0, 0]’;

bn = [0, s_1 , 0, 0, s_2, 0]’; sh = [0 , 0 , s_1, 0, 0, s_2]’; for i = 1: ndof*2 Ax(i) = Ax(i) + ax(i)*wt*h/2; Bn(i) = Bn(i) + bn(i)*wt*h/2; Mn(i) = Mn(i) + sh(i)*wt*h/2; end end %--------------GLOBAL FORCE VECTOR-----------------------------------function[ff] = shearglobalvector(F,Fel,locmatrix,nel) N = length(locmatrix); % # of dofs per element for i = 1:N

%

j = locmatrix(i); F(j)= F(j)+ Fel(i);

gives local vector row position % global position %assembled matrix

end

75

Bibliography [1] Carlos Felippa, http://www.colorado.edu/engineering/CAS/courses. d/NFEM.d/NFEM.Ch09.d/NFEM.Ch09.pdf [2] H. Saffari and R. Tabatabaei:“A Finite Circular Arch Element Based on Trigonometric Shape Functions”Mathematical Problems in Engineering. Vol. 2007, Article ID 78507 [3] P.Lagace,ocw.mit.edu/.../16-20Structural-MechanicsFall2002/ 3A499B9B-094C-4923-BFC9-9F5CB06E3E2C/0/unit13.pdf [4] S.Timoshenko and J.N, Goodier, Theory of Elasticity, McGraw-Hill Co.,Tokyo [5] J.T.Oden, Mechanics of Elastic Structures, MacGraw-Hill Co.1967 [6] Young W. Kwon, Hyochoong Bang, The Finite Element Method Using Matlab, CRC Press NY.1997 [7] G.Nikishkov, www.u-aizu.ac.jp/∼niki/feminstr/introfem/introfem. html-12k [8] Thomas J.R.Hughes, The Finite Element Method, Linear and Dynamic Finite Element Analysis, Dover Publ.,NY.2000 [9] Gilbert Strang, George J.Fix,An analysis of the Finite Element Method, Prenctice Hall, New Jersey [10] Hinton Owen, An introduction to Finite Element Computations, Pineridge Press Ltd., Swansea, 1979 [11] Gilbert W. Stewart, Afternotes on Numerical Analysis, SIAM, 1996 [12] L.D.Landau ,E.M Lifshitz, Theory of Elasticity, Butterworth-Heinemann, 1986 [13] University of washington, Mechanics of Materials Laboratory http:// courses.washington.edu/me354a/Curved%20Beams.pdf 76

[14] Alexander Blake, Handbook of Mechanics, Materials, and Structures, WileyIEEE, 1985 [15] A.O.Lebeck and J.S.Knowlton, “A finite element for three-dimensional deformation of a circular ring” International Journal for Numerical Methods in Engineering. Vol. 21, 421 - 435, 1985 [16] David. V. Wallerstein, A Variational Approach to Structural Analysis , John Wiley and Sons, 2002 [17] Z. Friedman, J. B. Kosmatka,“An accurate two-node finite elelemnt for shear deformable curved beams” International Journal for Numerical Methods in Engineering. Vol.41, 1998. 473-498. [18] Seok-Soon Lee, “Development of a new curved beam element with shear effect”, Engineering Computations Vol.13 No. 2/3/4, pp. 9 - 25, 1996. [19] Henryk Stolarski and Ted Belytschko, “Shear and Membrane Locking in Curved C(0) Elements”, Computer Methods in Applied Mechanics and Engineering, Vol. 41, pg 279 - 296 , 1983 [20] A.F. Saleeb and T.Y. Chang,“ On the hybri-mixed formulation of C(0) curved beam elements ”, Computer Methods in Applied Mechanics and Engineering, Vol. 60, 95 - 121, 1987. [21] A. E. H. Love, A Treatise on Mathematical Theory of elasticity, Dover Publ.,New York. [22] Fumio Kikuchi, “On the validity of the finite element analysis of circular arches represented by an assemblage of beam elements” Computer Methods in Applied Mechanics and Engineering,Vol 5. pg 253-276, 1975. [23] David J. Just, “Circularly curved beams under plane loads”, Journal of the Strucutral Division, Proceedings of the American Society of Engineers, Vol. 108, No.ST8, August 1982. [24] P.Lee and H.Sin,“Locking-free curved beam element based on curvature”, International Journal for Numerical Methods in Engineering vol. 37, no. 6, pp. 9891007, 1994. [25] P. Raveendranath, G. Singh, and B. Pradhan “A two-noded locking-free shear flexible curved beam element ” International Journal for Numerical Methods in Engineering, vol. 44, no. 2, pp. 265 - 280, 1999. 77

[26] G. Prathap and G. R. Bhashyam,“ Reduced integration and the shear flexible beam element”, International Journal for Numerical Methods in Engineering, Vol. 18, pp. 195-210, 1982 [27] Irbmeusz Kreja and Zbignifw Cywi´nski, “Is reduced Integration just a numerical trick” Computers & Structures.Vol. 29, No. 3, pp. 491-496, 1988 [28] A. Krishnan and Y. J. Suresh, “A simple cubic linear element for static and free vibration analysis of curved beams”, Computers & Structures, vol. 68, no. 5, pp. 473489, 1998. [29] RoyMech,http://www.roymech.co.uk/Useful Tables/Beams/ Beam energy methods.html [30] Zienkiewicz et al.“A simple and efficient element for axisymmetric shells ”, International Journal for Numerical Methods in Engineering, Vol. 11, 10 , pp.1545 - 1558, 2005 [31] R. A. Day and D. M. Potts, “Curved Mindlin beam and axi-symmetric shell elements - A new approach” International Journal for Numerical Methods in Engineering, Vol. 30, 1263 -1274 ,1990

78

Declaration of authorship I hereby declare

• that I have written this thesis without any help from others and without the use of documents and aids other than those stated in the bibliography, • that I have mentioned all used sources and that I have cited them correctly.

Signature:

Date:

79