1

Chapter 3: BASIC ELEMENTS Section 3.1: 3 1: Preliminaries (review of solid mechanics)

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

Outline
Most structural analysis y FE codes are displacement p based In this chapter we discuss interpolation methods and elements l t based b d on di displacement l ti interpolations t l ti Stiffness matrix formulations will be presented Shortcomings and restrictions of the elements due to the interpolations used will be discussed We will review the governing equations (for solids elastic bodies) to help us understand the solution methods and accuracy

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

Review of Solid Mechanics

The analysis y of any y solid elastic body y has to define and develop the following quantities and/or relations
Stress St Strain (strain-displacement relations) Constitutive Properties (Stress-Strain (Stress Strain relations) Compatibility Equilibrium Equations Boundary Conditions

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

Stress
Stresses are distributed internal forces that result from externally applied forces
F1 B1 F2 F1 F

F3 F5 F4 F5

Note: There are two types of forces: Surface forces that act on an area of external surface and body force that acts on the volume of the body
R.T. Haftka EML5526 Finite Element Analysis University of Florida

Stress/Force acting on a surface

A force acting on a surface can be resolved into two components: One tangential to the surface (shear force) and the other normal to the surface

Force is a 1st order tensor (vector) Stress is a 2nd order tensor. Why is it a tensor?

Needs two vectors to specify it it. One is the vector of the internal force and the other is the outward normal of the defining g area

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

Representation of stress at a point

In 3 3-D D space the stress at a point is denoted by the stresses acting on three mutually perpendicular planes at that point Often we use a simple infinitesimal rectangular solid and indicate the stresses on the faces of that solid Stress terms have two indices.
The first indicating the plane on which it acts The second indicating the direction in which it acts

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

Stress at a point
State of stress at a point in 3-D Cartesian Coordinates
yy
Matrix representation of state of stress y x

yz zy

yx

xy xx

zx xz zz

xx xy zx xy yy yz xz yz zz

Sign Conventions: Normal stresses are positive when acting outward from a surface (tension) Shear stresses are positive when they act in the +ve direction on a positive face and ve direction on a ve face
R.T. Haftka EML5526 Finite Element Analysis University of Florida

Stress on an arbitrary plane (2-D)

We often need to enforce stress boundary conditions on surfaces that are not always rectangular Let cos() = l and sin()= m If length of side BC = A, then length of sides OC=Al OC Al and OB= OB Am If we write the force equilibrium q in X and Y-directions, we have
C Y N

xy

/2
O

XA = x ( Al ) + xy ( Am)
Which simplifies to

xy y y

X = x l + xy m Y = y m + xy l
EML5526 Finite Element Analysis University of Florida

R.T. Haftka

Stress on an arbitrary plane (3-D) (3 D)

The forces (per unit area) in X, Y and Z-directions on an arbitrary bit plane l defined d fi d by b its it normal l vector t N is i given by

X = l x + m xy + n zx

y B N

Y = l xy + m y + n yz Z = l xz + m yz + n z
Where l, m and n are the direction cosines of the normal vector of the x arbitrary plane

O C z

l = cos Nx, m = cos Ny, n = cos Nz

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

10

Stress transformations in 3D
Stress transformation in 3D space can be defined using the directions cosines as follows. follows

x ' = l12 x + m12 y + n12 z + 2m1l1 xy + 2m1n1 yz + 2n1l1 zx

2 2 y ' = l22 x + m2 y + n2 z + 2m2l2 xz + 2m2 n2 yz + 2n2l2 zx 2 2 z ' = l32 x + m3 y + n3 z + 2m3l3 xy + 2m3n3 yz + 2n3l3 zx

x ' y ' = l1l2 x + m1m2 y + n1n2 z + (l1m2 m1l2 ) xy

+ (m1n2 n1m2 ) yz + (n1l2 l1n2 ) zx

y ' z ' = l2l3 x + m2 m3 y + n2 n3 z + (l2 m3 m2l3 ) xy y

+ (m2 n3 n2 m3 ) yz + (n2l3 l2 n3 ) zx

x ' z ' = l3l1 x + m3 m1 y + n3n1 z + (l3m1 m3l1 ) xy

+ (m3 n1 n3m1 ) yz + (n3l1 l3 n1 ) zx
R.T. Haftka EML5526 Finite Element Analysis University of Florida

11

Direction cosines in 3D
The direction cosines l.m and n between the new coordinate axes x, , y and z and the original coordinate system x, y and z are defined as follows

x x y' z l1 l2 l3

y m1 m2 m2

z n1 n2 n3

Since the transformation is orthogonal, the direction cosines must satisfy the following properties

l12 + m12 + n12 = 1

2 l12 + l2 + l32 = 1

l1l2 + m1m2 + n1n2 = 0 l1m1 + l2 m2 + l3m3 = 0

Where l1 =cos xx, Where, m2= cos yy

R.T. Haftka EML5526 Finite Element Analysis University of Florida

12

Equilibrium Equations (2-D) (2 D)

The 2-D stresses are shown on a volume given of length dx and dy in Xand Y-directions and unit thickness in the Z Zdirection

y yx + d yx
xy
dy
FX

x
Summation of forces in X-direction

x + d x

dx

FY

yx
=0
University of Florida

( x + d x )dy y x dy y + ( xy + d xy )dx xy dx + FX
R.T. Haftka EML5526 Finite Element Analysis

13

Equilibrium of forces in X-direction X direction

( x + d x )dy x dy + ( xy + d xy )dx xy dx + FX dxdy = 0 ( x + d x )dy d x dy d + ( xy + d xy )dx d xy dx d + FX dxdy d d =0
d x =
d xy
1 1 2 2

x dx x xy = dy y y

d xy d x dxdy + dxdy + FX dxdy = 0 dx dy

x xy x + y + FX
R.T. Haftka

dxdy = 0
University of Florida

EML5526 Finite Element Analysis

14

2-D 2 D Equilibrium Equations

The force equilibrium provide the relations shown below referred to as differential equation of equilibrium

x xy + +X=0 x y xy x + y y +Y =0

Establishing moment equilibrium by the same method will provide the condition for symmetry p y y of the stress tensor

xy y =y yx
R.T. Haftka EML5526 Finite Element Analysis University of Florida

15

Strain
Why do we need the strain measures? Will displacement not suffice? Strain better q quantifies the deformation of the body and eliminates rigid body motion/ rotation y general g terms is a measure of Strain in very relative deformation
Relative to what?
Undeformed body : Lagrangian strain Deformed body: Eulerian strain

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

16

Strain-Displacement Strain Displacement Relations

For uniaxial (1-D) case: l l
y
P Q
dy

l = l

P dx Q

P, P

du dv

u
P

P=P(x+u,y+v,z+w)

u =ui + vj + wk P P(x y z) P=P(x,y,z

R.T. Haftka EML5526 Finite Element Analysis

University of Florida

17

Strain 2-D 2 D : Normal strain

y
P Q
dy

P P P,

du dv

du =
dv =

P dx Q

u dx x
v dx x

x =
2

dx dx dx
2 2

Q 1 u u v x = 1+ 2 + + x x x 2
2

u v 2 dx + dx + dx dx x x x = dx

x = 1+ 2

u u v + + 1 x x x
2 2

2 2 u 1 u v x = + + x 2 x x

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

18

STRAIN: Shear Strain

O O
dy dx dx

dx

u dx Q x

B v dx x Q

O
P

dx

dy

Q P

1
dy

= 1 + 2

x v P dy y A

u dy y

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

19

Shear strain
With a bit of trigonometry (see for example, Allen and Haisler, Introduction to Aerospace Structural Analysis, p.60)
v u u v u u + + + x y x x y y

xy y =

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

20

For small displacements

The normal and shear strains are expressed as

x =

u x

y =

v y

xy =
1 v

v u + x y

Engineering Strains

+ xy = x y 2

Tensor strains

Shear strain definitions are of two forms. The above form is referred to as engineering strains. The alternate form referred as tensorial strains have a factor of applied to engineering strains. To apply coordinate transformations need the tensor form.

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

21

Compatibility
Deformation must be such that the pieces fit together without any gaps or overlap. overlap Why is this an issue? In 2-D we require only 2 displacements u, and v to describe deformation, but have three strain quantities x,y, and xy. This implies only two of the three strain terms are independent.

x =

u x

y =

v y

xy =

v u + x y

2 2 x y 2 2 u 2 v ( xy ) = 2 + 2 y 2 + x 2 xy y x x y

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

22

Stress-Strain Stress Strain Relations

The stress-strain relations in solid mechanics is often referred to as the Hookes Hooke s Law Hookes law of proportionality stated as extension is proportional to the force refers to the axial extension of a bar under an axial force This can be extended to 3-D stress/strain state referred to as the G Generalized li d H Hookes k L Law relates l t the th components t of f the th 3-D 3 D stress t state to 3-D strains as follows.

x c11 c y 21 c31 z = xy c41 yz c51 zx c61

R.T. Haftka

c12 c22 c32 c42 c52 c62

c13 c23 c33 c43 c53 c63

c14 c24 c34 c44 c54 c64

c15 c25 c35 c45 c55 c65

c16 x c26 y c36 z c46 xy c56 yz c66 zx

University of Florida

EML5526 Finite Element Analysis

23

Generalized Hookes Law

In the most general form the generalized Hookes Law requires 36 constants to relate the terms of a 3-D Stress state to its corresponding 3-D 3 D strain state for an elastic material However, from symmetry of the strain energy terms, it can be shown that cij = cji This reduces the number of unknown constants to 21

x c11 c y 12 c13 z = xy c14 yz c15 zx c16

R.T. Haftka

c12 c22 c23 c24 c25 c26

c13 c23 c33 c34 c35 c36

c14 c24 c34 c44 c45 c46

c15 c25 c35 c45 c55 c56

c16 x c26 y c36 z c46 xy c56 yz c66 zx

University of Florida

EML5526 Finite Element Analysis

24

Hookes Hooke s Law: Orthotropy and Isotropy

If we assume the x, y and x coordinates provide the planes of symmetry we can further reduce the number of constants to 9 9.
x c11 c y 12 z c13 = xy 0 yz 0 0 zx c12 c22 c23 0 0 0 c13 c23 c33 0 0 0 0 0 0 c44 0 0 0 0 0 0 c55 0 0 x 0 y 0 z 0 xy 0 yz c66 zx

This corresponds to an fully orthotropic material

Isotropy assumes that there is no directional variation on property. Using this argument we can obtain

c11 = c 22 = c33 c12 = c13 = c 23 c 44 = c55 = c66

R.T. Haftka EML5526 Finite Element Analysis University of Florida

25

Hookes law : Engineering Elastic Constants

The two engineering elastic constants used to relate stress to strain for isotropic materials are the Elastic modulus, E and the Poissons ration . For uniaxial loading, strain in the loading direction obtained from Hookes law, states

x =

x
E

Transverse to loading direction y = z = x =

x
E

The relation between the shear stress component and its corresponding shear-strain shear strain component is called the modulus of rigidity or modulus of elasticity in shear and is denoted by the letter G. E G= = G= 2(1 + )
R.T. Haftka EML5526 Finite Element Analysis University of Florida

26

Generalized Hookes Law from Isotropic M t i l Materials

The generalized Hookes law expressed in engineering elastic constants
1 x ( y + z ) E 1 y = y ( x + z ) E 1 z = z ( x + y ) E

x =

] ]

xy = yz zx

2(1 + ) xy G E yz 2(1 + ) yz = = G E 2(1 + ) = zx = zx G E

xy

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

27

Plane state of stress

There are a large class of problems for which the stresses normal lt to th the plane l of f the th solid lid are absent b t or negligibly li ibl small. If we assume that the stresses are restricted to the xy plane, then

z = xz = yz = 0
This simplifies p the stress strain relationship p to the form as shown below.
1 x y E 1 y = y x E ( x + y ) z = E

x =

x =

E x + y 1 2 E y = y + x 1 2

] ]

xy = xy

xy

G = G xy

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

28

Plane strain
Strains that deform the body normal to the reference plane are absent or are negligible

z = xz = yz = 0
This indicates that the stress normal to the p plane of strain is dependent on the stresses in the plane of the 1 strain z = z ( x + y ) = 0 z = ( x + y ) E Substituting z into other strain expressions we obtain

x =

1 + (1 ) x y E 1 + (1 ) y x y = E

] ]

E (1 ) x + y (1 + )(1 2 ) E (1 ) y + x y = (1 + )(1 2 ) E z = + (1 + )(1 2 ) x y

x =

] ]

xy =

xy
G

xy = G xy

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

29

Conversion from plane strain to plane stress and d vice-versa i

The solution obtained for the stress and strains in plane stress and plane strain states are qualitatively similar. To use a p plane strain solution for a p plane stress or vice versa, we simply interchange the appropriate constants as shown below For plane stress the expressions in E, For plane stress the expressions in E*,

E* * E= , = 1 * 1 *

or E * = E

1 + 2 (1 + )2

and =

1 +

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

30

Solution of 3D Elasticity Problems

3 6 6 ---- x , y , z , xy , yz , zx 15 Equilibrium equations Strain-displacement equations Hookes law equations

x , y , z , xy , yz , zx
u,v, w

3 Equilibrium equations p y 3 Compatibility 6 Hookes law equations ----12 x , y , z , xy , yz , zx x , y , z , xy , yz , zx

3 Equilibrium equations in displacements ----u , v, w 3 3 3 Equilibrium equations Compatibility equations expressed in stress

3 6

---- x , y , z , xy , yz , zx 9

Equilibrium equations Strain displacement equations

u , v, w

---- x , y , z , xy , yz , zx 6

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

31

Principle of superposition
Effects of several forces acting together are equal to the combined effect of the forces acting separately separately. This is valid only when
The stresses and displacements are directly proportional to the load The geometry and loading of the deformed object does not differ significantly from the undeformed configuration g
R.T. Haftka EML5526 Finite Element Analysis University of Florida

32

Energy Principles
Strain Energy Density: When an elastic body is under the action of external forces, the body deforms and work ork is done b by these forces forces. The work ork done b by the forces is stored internally by the body and is called the strain energy. Let us consider the unit element of volume dxdydz with only the normal stress x acting on it it. The work done done, or work stored in the element is
z D C D C
x = x

x
A

B B A

B x

u d u + dx dydz x x x =0 =
x = x x =0

x = x
x =0

d (u )dydz
x

u dxdydz x

y
R.T. Haftka EML5526 Finite Element Analysis University of Florida

33

Strain Energy
Using Hookes law
x = x x =0

u = x = x x E x
1 x d x dxdydz = dxdydz E 2 E
2

Work done =

For shear stresses, stresses it can be similarly shown that the work done is

The h strain energy stored d in an element l dxdydz under d a general three dimensional stress system is calculated as

12 dxdydz 2G

dU =
R.T. Haftka

1 ( x x + y y + z z + xy xy + yz yz + zx zx )dxdydz 2
EML5526 Finite Element Analysis University of Florida

34

Strain Energy Density

The strain energy gy density y refers to strain energy gy per p unit volume

dU 0 =

1 ( x x + y y + z z + xy xy + yz yz + zx zx ) 2
dU 0 =

Using principal stresses and strains, this can be expressed as

1 ( 11 + 2 2 + 3 3 ) 2 1 2 dU 0 = I1 2(1 + ) I 2 2E I1 = 1 + 2 + 3

I 2 = 1 2 + 2 3 + 1 3

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

35

3 2 Interpolation Functions 3.2

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

36

Introduction to Interpolation Functions

Interpolation is to devise a continuous function that satisfies prescribed conditions at a finite number of points points. In FEM, the points are the nodes of the elements & the prescribed conditions are the nodal values of the field variable Polynomials P l i l are th the usual l choice h i f for FEM

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

37

Polynomial interpolation
The polynomial function (x) is used to interpolate a field variable based on its values at n n-points points

( x) = ai xi

or

X = 1 x x

i =0

= X {a}

... x n and

{a} = a0

a1

a2

... a n

The number of terms in the polynomial is chosen to match the number of given quantities at the nodes. With one quantity per node node, we calculate ais s using the n n-equations equations resulting from the expressions for i at each of the n-known points

( x j ) = ai x ij
i =0
R.T. Haftka

{e } = [A ]{a}
EML5526 Finite Element Analysis

{a} = [A]1 {e }
University of Florida

38

Shape Functions or Basis Functions

Traditional interpolation takes the following steps 1. Choose a interpolation function 2. Evaluate interpolation function at known points 3. Solve equations to determine unknown constants

= X{a}

{e } = [A ]{a}

{a} = [A]1 {e }

= X{a}

In FEM we are more interested in writing in terms of the nodal values

= X{a}

{e } = [A ]{a}
= X [ A]
1

{a} = [A]1 {e }
= N {e }
N = X [ A]
University of Florida

{e }

R.T. Haftka

EML5526 Finite Element Analysis

39

Degree of Continuity
In FEM field quantities are interpolated in piecewise fashion over each element This implies that is continuous and smooth within the element However, may not be smooth between elements m An interpolation function with C continuity provides a continuous variation of the function and up to the mm derivatives at the nodes
For example in a 1-D interpolation of f(x) C continuity indicates that f is continuous at the nodes and f f,x is not continuous. continuous 0 If the displacement u(x) is C then displacements are continuous between elements, but the strains are not (bar elements)
0

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

40

Degree of Continuity
Function 1 is C0 continuous while 2 is C1 continuous
d dx

(x)
1 2

d1 dx

C0

d 2 dx

C1 x1 x2 x x0

C0 x1

C1 x2 x

x0

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

41

Example: Deriving a 1D linear interpolation shape h f function ti

N {e } each interpolation From = function is zero at all the dofs except one. This can allow us to derive interpolation functions one at a time. For F linear li i interpolation t l ti b between t x1 and d x 2, N1(x1)=1, N1(x2)=0, N1=a1x+a2. So obviously, b i l N1=1-(x-x 1 ( 1)/L, )/L L L=x2-x1.

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

42

C0 Interpolation -1D Element

1

x1
1

x2
1

x3

N1 =

(x2 x )(x3 x ) (x2 x1 )(x3 x1 ) (x1 x )(x3 x ) (x1 x2 )(x3 x2 ) (x1 x )(x2 x ) (x1 x3 )(x2 x3 )

N1 =
1

x2 x x 2 x1 x x1 x 2 x1
x1

N2 =

x2

x3
1

N2 =

1 x1 1

2 = N 1 2 x2 2

N3 =

x1

x2

x3

1
x1 1

2
x2 2

3
x3 3

1 = N 2 3

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

43

Lagrange Interpolation Formula

Shape functions shown for the C0 interpolations are special forms of the Lagrangian interpolation functions
f ( x) = N k f k
k =1 n

Nk =

( x1 x )( x2 x ) ...[ xk x ]... ( xn x ) ( x1 xk )( x2 xk ) ...[ xk xk ]... ( xn xk )

In above expressions for Nk the terms in square brackets are omitted

R.T. Haftka EML5526 Finite Element Analysis University of Florida

44

Properties of C0 Interpolation Shape F Functions ti All shape functions Ni, along with function are polynomials of the same degree For any shape function Ni, Ni = 1 at node i ( (x=x d zero at t all ll other th nodes d j j, ( (xj xi) i) and C0 shape functions sum to one

N
k =1
R.T. Haftka

=1
University of Florida

EML5526 Finite Element Analysis

45

C1 Interpolation p
Also called Hermitian interpolation (Hermite polynomials) Use the ordinate and slope information at the nodes to interpolate
,x2 ,x1
C0 interpolation curve

C1 interpolation curve

1 x

2 x
University of Florida

R.T. Haftka

EML5526 Finite Element Analysis

46

Hermitian interpolation used for beam elements l t

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

47

2-D 2 D and 3-D 3 D Interpolation

The 2-D and 3-D shape p functions follow the same procedure as for 1-D We now have to start with shape functions that have two or more independent terms terms. For example a linear interpolation in 2-D from 3 nodes will require an interpolation function

f ( x, y ) = 1 x

y a1

a2

a3

If there are two or more components (e (e.g., g u, v and w displacements) then the same interpolation function is used for all components
R.T. Haftka EML5526 Finite Element Analysis University of Florida

48

Principle of Virtual Work

The principle of virtual work states that at equilibrium q the strain energy gy change g due to a small virtual displacement is equal to the work done by the forces in moving through the virtual displacement. A virtual displacement is a small imaginary change in configuration that is also a admissible displacement An admissible displacement satisfies kinematic boundary y conditions Note: Neither loads nor stresses are altered by the virtual displacement.
R.T. Haftka EML5526 Finite Element Analysis University of Florida

49

Principle of Virtual Work

The principle of virtual work can be written as follows T T T { } { }dV = {u} {F }dV + {u} {}dS The same can be obtained by the Principle of Stationary Potential Energy The total potential energy of a system is given by

= U W = U + V = 0
U is strain energy, W is work done, or V is potential of the forces

= U W

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

50

Element and load derivation

{u} = Interpolation {u} = [ N ]{d} u v w Strain displacement { } = [ B ]{d} [ B ] = [ ][ N ] Vi t l {u} = {d} [ N ] and { } = {d} [ B ] Virtual Constitutive law { } = [ E ]{ } Altogether
T T T T T T

{d} ( [ B ] [ E ][ B ] dV {d} [ B ] [ E ]{ 0 } dV + [ B ] { 0 } dV T T [ N ] {F} dV [ N ] { } dS ) = 0

T T T T
EML5526 Finite Element Analysis

R.T. Haftka

University of Florida

51

Stiffness matrix and load vector

Equations of equilibrium [ k ]{d} = {re } Element stiffness matrix

[ k ] = [ B ] [ E ][ B ] dV
T

Element El tl load d vector t

{re} =[ B] [ E]{0} dV[ B] {0} dV+ [ N] {F} dV+ [ N] {} dS ) Loads due to initial strain, initial stress, y forces and surface tractions body
T T T T
R.T. Haftka EML5526 Finite Element Analysis University of Florida

Plane Problems: Constitutive Equations

Constitutive C tit ti equations ti for f a linearly li l elastic l ti and d isotropic material in plane stress (i.e., z=xz=yz=0):

52

Initial thermal strains

x 0 = y 0 = T , xy0 = 0

In matrix form, form

where

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

53

Plane Problems: Strain-Displacement Relations

.

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

Plane Problems: Displacement Field Interpolated

54

From the previous two equations,

where B is the strain-displacement strain displacement matrix.
R.T. Haftka EML5526 Finite Element Analysis University of Florida

Constant Strain Triangle (CST)

55

The node numbers sequence must go counter clockwise p field so strains are constant! Linear displacement

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

Constant Strain Triangular Element

CST ELEMENT

Decompose two-dimensional domain by a set of triangles. Each triangular element is composed by three corner nodes. Each element shares its edge g and two corner nodes with an adjacent j element Counter-clockwise or clockwise node numbering Each node has two DOFs: u and v displacements interpolation using the shape functions and nodal displacements. Displacement is linear because three nodal data are available available. Stress & strain are constant. y v3 3 v1 1 u1 v2 2 u2 x
56

u3

CST ELEMENT cont.

Displacement Interpolation
Since two-coordinates are perpendicular, u(x,y) and v(x,y) are separated. u(x,y) needs d t to be b i interpolated t l t di in t terms of f u1, u2, and d u3, and v(x,y) in terms of v1, v2, and v3. interpolation function must be a three term polynomial in x and y. Since we must have rigid body displacements and constant strain terms in the interpolation function, the displacement interpolation must be of the form
u( x, y ) = 1 + 2 x + 3 y v ( x, y ) = 1 + 2 x + 3 y

The goal is how to calculate unknown coefficients i and i, i = 1, 2, 3, in terms of nodal displacements.

u( x, y ) = N1( x, y )u1 + N2 ( x, y )u2 + N3 ( x, y )u3

57

CST ELEMENT cont.

Displacement Interpolation
x-displacement: Evaluate displacement at each node
u( x1, y1 ) u1 = 1 + 2 x1 + 3 y1 u( x , y ) u = + x + y 2 2 2 1 2 2 3 2 u( x3 , y 3 ) u3 = 1 + 2 x3 + 3 y 3

In matrix notation
u1 u2 = u 3 1 x1 1 x2 1 x3 1 y1 y 2 2 y 3 3

Is the coefficient matrix singular?

58

CST ELEMENT cont.

Displacement Interpolation
1 2 = 3 1 x1 1 x2 1 x3
1 u1 f1 f2 y1 1 b1 b2 y 2 u2 = 2 A c1 c2 y 3 u3

u1 f3 b3 u2 c3 u3

where f1 = x2 y 3 x3 y 2 , b1 = y 2 y 3 , c1 = x3 x2 f = x y x y , b = y y , c = x x 2 3 1 1 3 2 3 1 2 1 3 f3 = x1y 2 x2 y1, b3 = y1 y 2 , c3 = x2 x1 Area:

1 x1 1 A = det 1 x2 2 1 x3 y1 y2 y3

59

CST ELEMENT cont.

1 (f u + f u + f3u3 ) 2A 1 1 2 2 1 2 = (b u + b2u2 + b3u3 ) 2A 1 1 1 3 = (c u + c2u2 + c3u3 ) 2A 1 1 1 =

Insert to the interpolation equation

u( x, y ) = 1 + 2 x + 3 y 1 [ (f u + f u + f3u3 ) + (b1u1 + b2u2 + b3u3 )x + (c1u1 + c2u2 + c3u3 )y ] 2A 1 1 2 2 1 (f + b1x + c1y )u1 = N1(x,y) (x y) 2A 1 1 + (f + b2 x + c2 y )u2 N2(x,y) 2A 2 1 + (f + b3 x + c3 y )u3 2A 3 N (x,y) =
3

60

CST ELEMENT cont.

Displacement Interpolation
A similar procedure can be applied for y-displacement v(x, y).
u( x, y ) = [N1 N2 u1 N3 ] u2 u 3 v1 N3 ] v 2 v3

v ( x, y ) = [N1 N2

1 N1( x, y ) = (f + b1x + c1y ) 2A 1 N ( x, y ) = 1 (f + b x + c y ) 2 2 2 2A 2 1 N3 ( x, y ) = (f + b3 x + c3 y ) 2A 3 Shape Function

N1, N2, and N3 are linear functions of x- and y-coordinates. Interpolated displacement changes linearly along the each coordinate direction.

61

CST ELEMENT cont.

Displacement Interpolation
Matrix Notation
u = {u} v N1 0 N2 0 N1 0 0 N2 N3 0 u1 v 1 0 u 2 N3 v2 u3 v3

{u( x, y )} = [N( x, y )]{q} [N]: 26 matrix, {q}: 61 vector. For a given point (x,y) within element, calculate [N] and multiply it with {q} to evaluate displacement at the point (x,y).

62

CST ELEMENT cont.

Strain Interpolation
differentiating the displacement in x- and y-directions. differentiating shape function [N] because {q} is constant.
xx yy
3 u = = N x y u ( , ) i i x x i =1 3 v = = ( , ) N x y v i i y y i =1

N xi ui = i =1 N yi v i = i =1
3

i ui 2A i =1 3 i vi 2A i =1

b c

xy

u v + = y x

ci b ui + i v i 2A i =1 i =1 2 A

63

CST ELEMENT cont.

Strain Interpolation
u x b1 0 b2 v 1 0 c1 0 { } = = y 2 A c1 b1 c2 u v + y x u1 v 1 0 u2 c3 [B]{q} v 2 b3 u 3 v 3

0 c2 b2

b3 0 c3

[B] matrix t i is i a constant t t matrix t i and d depends d d only l on th the coordinates di t of f the three nodes of the triangular element. the strains will be constant over a given element

64

65

Constant Strain Triangle (CST): Stiffness Matrix Strain-displacement relation, =Bd

A is the area of the triangle and xij=xi- xj. (textbook has results for a coordinate system with x aligned with side 1-2

From the general formula

k = BT EBtA

where t: element thickness (constant)

R.T. Haftka EML5526 Finite Element Analysis University of Florida

CST ELEMENT cont.

Strain Energy: U (e ) = h {}T [C]{} dA(e )
2
A

h (e ) T (e ) {q } [B]T 63 [C]33 [B]36 dA {q } 2

A

1 (e ) T (e ) {q } [k ]66 {q( e ) } 2

Element Stiffness Matrix: [k( e ) ] = hA[B]T [C][B] Different from the truss and beam elements elements, transformation matrix [T] is not required in the two-dimensional element because [k] is constructed in the global coordinates.

The strain energy of the entire solid is simply the sum of the element strain energies
U=

U (e ) =
e=1

NE

NE assembly 1 1 (e ) T (e ) (e ) U = {Q s }T [K s ]{Q s } { q } [ k ]{ q } 2 2 e=1

66

BEAM BENDING EXAMPLE

-F 2 1 4 6 5 5m 7 8 9 F 10 1m

along l th the x-axis i and d linear along y-axis Exact Solution: xx = 60 MPa Max deflection vmax = 0.0075 0 0075 m

xx is constant

Max v = 0.0018

xx
x
67

Linear Strain Triangle (LST)

68

The element has six nodes and 12 dof.

R.T. Haftka EML5526 Finite Element Analysis University of Florida

69

Linear Strain Triangle (LST)

The quadratic displacement field in terms of generalized coordinates:

The linear strain field:

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

Bilinear Quadrilateral (Q4): CQUAD4 in NASTRAN

Q4 element has four
nodes and eight dof. dof Can be quadrilateral; but for now rectangle g .

70

Displacement field:
So, u and v are bilinear in x and y. Because of form, sides id are stiffer tiff th than di diagonals-artificial l tifi i l anisotropy! i t !
R.T. Haftka EML5526 Finite Element Analysis University of Florida

71

Q4: The strain fields

Strain field:

Observation 1: x f(x) Q4 cannot exactly model the beam where x x

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

Q4: Behavior in Pure Bending of a Beam

Observation 2: When 4 0, 0 x varies linearly in y desirable characteristic for a beam in pure bending because normal strain varies linearly y along g the depth coordinate. But xy0 is undesirable because there is no shear strain.

72

Fig. Fig (a) is the correct deformation in pure bending while (b)
is the deformation of Q4 (sides remain straight). Physical y interpretation: p applied pp moment is resisted by ya spurious shear stress as well as flexural (normal) stresses.
R.T. Haftka EML5526 Finite Element Analysis University of Florida

73

Q4: Interpolation functions

Easy E to t obtain bt i interpolation i t l ti f functions ti

where matrix N is 2x8 and the shape functions are

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

74

Q4: The Shape (Interpolation) Functions

N1=1, N2=N3=N4=0 at node 1, x=-a, y=-b, so u= N1 u1= u1
at that node. Similarly Ni=1 while all other Ns are zero at node i i.

Fig 3.6-2

See Eqn. 3.6-6 for strain-displacement matrix (= Nd=Bd). All in all all, Q4 converges properly with mesh refinement and works better than CST in most problems.
R.T. Haftka EML5526 Finite Element Analysis University of Florida

75

Coarse mesh results

Q4 element is over-stiff in bending. For the following
problem, deflections and flexural stresses are smaller than the exact values and the shear stresses are g greatly y in error:

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

BEAM BENDING PROBLEM cont.

Sxx Plot
Max v = 0.0051

Stress is constant along the x x-axis axis (pure bending) linear through the height of the beam Deflection is much higher than CST element element. In fact fact, CST element is too stiff. However, stress is inaccurate.
76

BEAM BENDING PROBLEM cont.

Caution:
In numerical integration, we did not calculate stress at node points. Instead, we calculate stress at integration points. Lets L t calculate l l t stress t at t th the b bottom tt surface f for f 3(1,1) 4(0,1) element 1 in the beam bending problem. Nodal Coordinates:1(0,0), 2(1,0), 3(1,1), 4(0,1) Nodal Displacements: u = [0, 0.0002022, -0.0002022, 0] [0, , 0.0002022, , 0.0002022, , 0] ] v=[ 2(1 0) 2(1,0) 1(0 0) 1(0,0) Shape functions and derivatives

N1 = ( x 1)( y 1) N 2 = x( y 1) N 3 = xy N 4 = ( x 1) y

N1 / x = ( y 1) N 2 / x = ( y 1) N 3 / x = y N 4 / x = y

N1 / y = ( x 1) N 2 / y = x N 3 / y = x N 4 / y = ( x 1)
77

BEAM BENDING PROBLEM cont.

At bottom surface, surface y = 0
N1 / x = 1 N 2 / x = 1 N 3 / x = 0 N 4 / x = 0
x
4 3

N1 / y = x 1 N 2 / y = x N 3 / y = x N 4 / y = ( x 1)
1 2

Strain 4 N I xx = u I = 1 0.0002022
I =1 4

u = [0, 0.0002022, -0.0002022, 0] v = [0, 0.0002022, 0.0002022, 0]

yy =
4

N I vI = 0.0002022 x + 0.0002022 x = 0 I =1 y

xy =

N I N vI + I u I = 0.0002022 0.0004044 x y I =1 x

Stress:

{ } = [C]{ } = {4.44, 1.33, 1.55} 107

78

RECTANGULAR ELEMENT
y-normal y normal stress and shear stress are supposed to be zero zero.

yy Plot
xx is a linear function of y alone yy if a linear function of x alone xy is a linear function of x and y

xy Plot

N1 / x = ( y 1)

N1 / y = ( x 1)

N 2 / x = ( y 1) N 2 / y = x N 3 / x = y N 3 / y = x N 4 / x = y N 4 / y = ( x 1) )
79

xx

4 N NI = uI yy = I v I I =1 y I =1 x

RECTANGULAR ELEMENT
Discussions
Cant represent constant shear force problem because xx must be a linear function of x. Even E if xx can represent t linear li strain t i i in y-direction, di ti the th rectangular t l element cant represent pure bending problem accurately. Spurious shear strain makes the element too stiff.

u = 1 + 2 x + 3 y + 4 xy v = 1 + 2 x + 3 y + 4 xy

Exact

xx = 2 + 4 y yy = 3 + 4 x xy = ( 3 + 2 ) + 4 x + 4 y
4 0

Rectangular element

80

RECTANGULAR ELEMENT
Two-Layer Two Layer Model
xx = 3.48107 vmax = 0.0053

81

BEAM BENDING PROBLEM cont.

Distorted Element
Max v = 0.004

As element is distorted distorted, the solution is not accurate any more more.

82

BEAM BENDING PROBLEM cont.

Constant Shear Force Problem
Max v = 0.0035

Sxx is supposed to change linearly along x x-axis. axis But But, the element is unable to represent linear change of stress along x-axis. Why? Exact solution: v = 0.005 m and xx = 6e7 Pa.
83

BEAM BENDING PROBLEM cont.

Higher-Order Higher Order Element?
8-Node Rectangular Element
4 8 7 3 6

u ( x, y ) = a0 + a1 x + a2 y + a3 x + a4 xy
2

+ a5 y 2 + a6 x 2 y + a7 xy 2
Strain

u ( x, y ) = a1 + 2a3 x + a4 y + 2a6 xy + a7 y 2 x
Can this element accurately represent pure bending and constant shear force problem?

84

BEAM BENDING PROBLEM cont.

8-Node 8 Node Rectangular Elements

Tip Displacement = 0.0075 m, Exact!

85

BEAM BENDING PROBLEM cont.

If the stress at the bottom surface is calculated, calculated it will be the exact stress value.

Sxx

Syy

86

Q6: Additional Bending Shape Functions

Q4: Artificial shear deformation under pure bending Additional shape functions to solve the issue
u = Ni (s, t )ui + (1 s 2 )a1 + (1 t 2 )a2 v = Ni (s, t )v i + (1 s 2 )a3 + (1 t 2 )a4
i =1 i =1 4 4

Bubble modes

Strain xx can vary linear along x-dir. Shear strain xy can vanish for pure bending Nodeless DOFs, a1, a2, a3, and a4, are condensed in the element level ( (total 12 DOFs) )
87

88

Modeling Bending with the Q6 Element

K dd Kad fd Kda d = Kaa a fa
1 { a } = K aa { fa K ad d} 1 1 [Kdd KdaK aaK ad ]{d} = { fd K daK aa fa }

Modeling the previous bending problem with Q6 elements gives the following stresses:

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

Quadratic Quadrilateral (Q8)

89

3.7-1

4 corner nodes and 4 side nodes and 16 nodal dof.

R.T. Haftka EML5526 Finite Element Analysis University of Florida

Quadratic Quadrilateral (Q8): Displacement

The displacement field, which is quadratic in x and y:

90

Two types of shape functions (=x/a, =y/b) :

The edges x=a deform into a parabola (i.e., quadratic displacement in y) (same for y=b)

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

Quadratic Quadrilateral (Q8): Strains

The strain field:

91

Strains have linear and quadratic terms. Hence, Q8 can

represent many strain states exactly. For example, p states of constant strain, bending g strain, etc.
R.T. Haftka EML5526 Finite Element Analysis University of Florida

92

3.11 Nodal loads

Consistent (work-equivalent) (work equivalent) loads
W = d Tre = u T FdV + u TdS

Mechanical loads: concentrated loads surface traction loads, traction, body forces.

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

93

Example: Beam under uniform loads

Normal forces are obvious. For moments
M 11 = 2 x 2 x3 qN N 21dx d = q( x + 2 )d dx L L 0 0 1 2 1 qL2 = qL ( + ) = 2 3 4 12
2 L L

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

Work equivalent (consistent) normal loads

Normal surface traction on a side of a plane element whose sides remain straight (q is force/length):

94

R.T. Haftka EML5526 Finite Element Analysis University of Florida

95

Loads on Quadratic sides

*

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

96

Distributed Shear Traction

Shear traction on a side of a plane element whose sides remain straight (q is force/length):

In (b), a Q4 element and two LSTs share the top midnode so that the nodal loads from Q4 and the right LST are combined.
R.T. Haftka EML5526 Finite Element Analysis University of Florida

Uniform Body Force

W Work-equivalent k i l t nodal d lf forces corresponding di t to weight i ht as a body force:

97

LST has no vertex loads and vertex loads of Q8 are upwards! The resultant in all cases is W, the weight of the element.
R.T. Haftka EML5526 Finite Element Analysis University of Florida

Since the p previous p plane elements have translational dof only, y

no moment can be applied to their nodes. The connection (a) of a beam and plane elements cannot t transmit it a moment t and d the th beam b can f freely l rotate. t t (Singular (Si l K!)

Connecting beam and plane elements

98

. the beam is extended. Rotational dof at A, B and C are In (b) associated with the beam elements only only. A plane element with drilling dof would also work but is not recommended.
R.T. Haftka EML5526 Finite Element Analysis University of Florida

99

Elements with Drilling DOF

Drilling dof: rotational dof about axis normal to the plane.
A CST with these added to each node has 9 dof. This dof allows twisting and bending rotations of shells under some loads to be represented. See Section 3.10

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

10 0

Stress calculation
Combining g Hookes law with strain-displacement p equation q
= E ( Bd 0 ) e.g. x E y = 2 1 xy 0 1 T 1 0 Bd T 0 0 0 (1 ) / 2

Stresses are most accurate inside elements

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

Improving stresses at nodes and b boundaries d i

One common technique q is averaging, g g, but

10 1

There are interpolation and extrapolation techniques that we will study y later
R.T. Haftka EML5526 Finite Element Analysis University of Florida

10 2

Examples of poor meshing

Do not create unnecessary y discontinuities!

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida

Example of discontinuities in unaveraged d stresses t

10 3

R.T. Haftka

EML5526 Finite Element Analysis

University of Florida