Fem Stress Analysis

Finite Element Method (FEM) Of Stress Analysis Finite Element Method Engineering Design always involved complicated shapes that are often composed of several different materials. Often the complexity of a part made it impossible to solve using classical method (also known as analytical method). Finite Element Method is the most popular numerical technique of solving these complex problems. The application of FEM ranges from stress analysis to heat transfer and fluid flow. For the purpose of this lecture, we will discuss application of FEM on Stress Analysis. Introduction If a load W, is applied at a point on the surface of a body, the resultant displacement / deformation on the other point of the surface is given by the function: ∆ = f(W) If systems of loads W1 , W2 ,ÉÉWn is applied to a body at points 1,2,Én, assuming that the system is linear and the resulting deformations and strains are small, the Principle of Superposition can be applied to determine the deformation ∆i = fi,1(W1) + fi,2(W2) + É+ fi,n(Wn) where ∆i is the deflection at point i measured in the direction of Wi. Combining all the equations and place it into matrix form, we have:

[∆] = [F] [W]

where ∆1

f11 f12 f13 L f1n

∆2 [∆ 1 ] = M M ∆n

[F ] =

W1

f 21 f 22 f 23 L f 2n

[W ] = W2

M M M f n1 f n 2 f n 3 L f nn

[∆] = [F] [W]

M Wn

eqn 1

Manipulating the equation 1 to solve in terms of load, [F]-1 [∆] = [F]-1 [F] [W] [F]-1 [∆] = [I] [W] [W] = [F]-1 [∆] [W] = [K] [∆]É. eqn 2

If equation 1 is used in the analysis, the method is known as the displacement or flexibility method. If equation 2 is used in the analysis, the method is known as the force or stiffness method. Both flexibility and stiffness matrices are symmetrical and same dimension n*n.

Finite Element Method

In Finite Element Method, the analysis begins with an approximation of the region of interest and its subdivision into numbers of meshes. These meshes are associated with nodes and element that becomes finite element. These elements can be triangular or quadrilateral element or even a complex shape element. Finite Element Analysis generates a series of nodes that is connected to each other by an element. A triangular element has 3 nodes while a quadrilateral has 4 nodes that are inter-connected. Every node in the element is represented by a mathematical equation with accommodating boundary condition. The principle of the method: Example:

a rectangular is subjected to a normal stress along part of itÕs longer edge and it is required to determine the stress distribution within it.

Load, W

The procedure is to split the rectangle into a series of element, preferably all of the same type, which are physically separated from each other except at the nodal point.

The continuum has now been given the from of a series of connected elements and can be regarded as structure to which the rules of structural analysis can be applied.

Stiffness Matrix of a Triangular Element When a triangular element is stressed in its plane, there are six possible components of displacement, two at each of the nodal points.

V2 U2 V3 U3 V1 U1

Y(v) x3 2 x2 3 y2

y3

1 X(u)

From the above displacement, Internal deformation can take place within an element (u , v). Both u and v are assumed to vary linearly in the x and y directions respectively. Based on the triangular element shown above, the internal deformation of each node in the triangular element is found to be: ui = α1 + xiα2 + yiα3 vi = α4 + xiα5 + yiα6 where u and v is the displacement in the x and y-direction respectively. By assembling the internal deformation of these 3 nodes (for triangular element) into matrix form:

u1 u2 u3 v1 v2 v3

=

1 1

x1 x2

y1 y 21

0 0

0 0

0 0

1 0

x3 0

y 31 0

0 1

0 x1

0 y1

0 0

0 0

0 0

1 1

x2 x3

y2 y3

α1 α2 α3 α4 α5 α6

Based on the local coordinate system, node 1 is located on the origin of the reference coordinate system, this concludes that:

x1 = y1 = 0 which leads to u1 = α 1 v1 = α 4

in simplified matrix form,

u1 u2 u3 v1

=

v2 v3

1 1

0 x2

0 y 21

0 0

0 0

0 0

1 0

x3 0

y 31 0

0 1

0 0

0 0

0 0

0 0

0 0

1 1

x2 x3

y2 y3

α1 α2 α3 α4 α5 α6

Further simplification leads to:

U V

=

[A] [α ]

eqn 3

where u and v are the components of the displacement related to the local frame α1É.αn are the components of the displacement related to some general system of coordinates. It will be called the generalized coordinate displacement components.

Internal strain Strain is defined by amount of elongation / stretch when subjected to a tensile load. It is shown by change in the length over the original length. The strain in the x direction, y direction and xy plane is shown:

εx =

δu δx

εy =

δv δy

γ xy =

u v

=

δu δv + δy δx

1 x y 0 0 0 0 0 0 1 x y

substituting these strain in each node, we have:

α1 M αn

εx εy

=

γ xy

0 1 0 0 0 0

α1

0 0 0 0 0 1 0 0 1 0 1 0

M αn

which is in the form of [ε] = [B] [α].............eqn 4 The evaluation of Internal Stresses The elastic relationship between stress and strains components for the plane stress is :

σx =

(ε

E 1− v2

x

+ vε y )

E (ε y + vε x ) 1− v2 E (1 − v ) τ xy = γ xy 2 1− v2

σy=

(

)

in matrix form:

σx σy τ xy

which is in the from of:

E = 1− v2

1 v v 1

0 0 1− v 0 0 2

εx εy γ ny

[σ] = [D] [ε] ..................eqn 5 but: [ε] = [B] [α]

∴ [σ] = [D] [B] [α].............eqn 6 Internal work δWi Internal work done onto the system is defined by : Internal work = internal strain * internal stress * volume σ δV ε t but ε=Bα

Virtual strain

ε=B α ε t = α t Bt ε t = [ I ] Bt δWi = [ I ] Bt σ δV

but σ = Dε = DBα

Unit virtual displacement in the direction of generalized coordinates

δWi = [ I ] Bt D B α ( T δx δy) δWi = T Bt D B α δx δy Wi = T Bt D B α ∫ ∫ δx δy Wi = T Bt D B α ( Area of element) Wi = T Bt D B ( Area of element) α

W = k ‡ i

where :

k = T B t D B ( Area of element)

External work A distance where the force is applied defines external work done onto the system. External work = displacement * applied force Assume that the resultant stress acting at the nodal points of the element in the direction of the generalized coordinate displacement components are:

S1 [S ] =

S2 M S6

We = α

t

S

We = S

From CastiglianoÕs theorem, the external work done on the system must equal with the internal work

Therefore, Wi = We k α = S ∴ S = k α

Force in the direction of generalized coordinate displacement

Generalized coordinate displacement

k is the stiffness of the element expressed in terms of the general coordinate displacement vector α

The relationship between the stiffness expressed in terms of the generalized and local coordinates

Total energy stored in the system is the same no matter what coordinate axes are used.

(

)

(

)

1 1 U W + U W + U W + ...... = α S + α S + α S + ...... 1 1 2 2 3 3 2 2 3 3 2 2 1 1 Ut W = αt S But: U V

[α ]

[A] [α ]

=

=

U

[A]−1

W = [ A −1 ]t k α

[α ] = [A]

U V

−1

t

V

W = [ A −1 ]t k [ A]−1

∴ u

t

W =

W = k

U V

t

U V

t

[ A −1 ]t S

U V

where k = [ A −1 ]t k [ A] −1

For example

As an illustration of the finite element method, let us consider the simple case of a thin square plate subjected in plane stress to two unit point loads applied at two corners, as shown in the figure:

A

B 1

Fixed location D

C 1

1,2

3,4 b

Triangular elements with 2 elements

a 7,8

5,6

If at this stage we ignore any boundary conditions then every nodal point can be assume to have two degree of freedom, i.e. it is capable of displacement in both the horizontal and vertical directions At each nodal point assign an odd number to the horizontal displacement and an even number to represent the vertical displacement.

7 1

1 2

7 1 0.625 − 0.125 1 2 − 0.125 0.125

5 3

8 4

2 5

− 0.5 0

0.375 − 0.25 − 0.125 − 0.125 0 0.125

0 0.5 − 0.25 5 3 − 0.5 E k = a (1 − v 2 ) 8 4 0.375 − 0.125 − 0.25 0.625 0 0.25 − 0.5 2 5 − 0.25 0 − 0.125 6 6 − 0.125 0.125

Local coordinate Global coordinate

2

(1,2)

With local coordinates a (7,8) 1

(5,6) 3

6 6

0.25 − 0.5

0 − 0.125

0.5 0

0 0.125

1 1 1 1 5 2

5 2 0.5 0

0 0.125

3 3

2 4

− 0.5 − 0.125

6 5

0 0.125

4 6 0.25 0

− 0.25 − 0.125

3 3 − 0.5 − 0.125 0.625 − 0.125 − 0.25 0.375 E k = b (1 − v 2 ) 2 4 0 0.125 − 0.125 0.125 0 − 0.125 0 − 0.25 0 0.5 − 0.5 6 5 0.25 0.625 4 6 − 0.25 − 0.125 0.375 − 0.125 − 0.5

1

3 (3,4)

(1,2) With local coordinates

b

2 For element a u 1 element a represent the displacement 7; u 2 element a represent the displacement 1; u 3 element a represent the displacement 5; u 4 element a represent the displacement 8; u 5 element a represent the displacement 2; u 6 element a represent the displacement 6;

(5,6)

For element b u 1 element a represent the displacement 1; u 2 element a represent the displacement 5; u 3 element a represent the displacement 3; u 4 element a represent the displacement 2; u 5 element a represent the displacement 6; u 6 element a represent the displacement 4;

The overall stiffness matrix by adding the matrix a and b according to global coordinate:

0.625 0

k=

E

− 0.5 − 0.25

0  1- v2 √ ↵ 0.375 − 0.125 − 0.125

0 − 0.5 − 0.25 0.625 − 0.125 − 0.125 0.375 0.375 − 0.125 − 0.125 0.625 0.625 − 0.125 − 0.125 0.375 0

0.375 0 0.25 − 0.5

− 0.125 − 0.125 − 0.25 − 0.5 0 0 0

0

0.375 0 − 0.25 − 0.5

0.625

0

0

0.625

− 0.5 − 0.25

− 0.125 − 0.125

− 0.125 − 0.125 − 0.25 − 0.5 0 0 0

0

− 0.5 − 0.125

− 0.25 0.375

0.625

0.375

0.375

0.625

Boundary conditions:

From the above figure, point a and D are fixed, i.e., the boundary condition u 1, u 2, u 7, and u 8 is 0. This leads to zero displacement in row and columns 1, 2, 7 and 8. In order to make the matrix calculation easier, rows and columns that has a zero displacement is removed from the matrix. Thus the new stiffness matrix:

k=

E

0.625

0.375

− 0.125 − 0.25

0.375

0.625

− 0.125

− 0.5

0.625 0

0 0.625

 1 - v 2 √ − 0.125 − 0.125 ↵ − 0.25 − 0.5

From the relationship: W = K∆ Where

-1 W=

u3

0 -1 0

and

∆=

u4 u5 u6

this leads to 4 simultaneous equations:

0.625 u 3 + 0.375 u 4 - 0.125 u 5 - 0.25 u 6 = - 1 0.375 u 3 + 0.625 u 4 - 0.125 u 5 - 0.5 u 6 = 0 - 0.125 u 3 - 0.125 u 4 + 0.625 u 5 = -1 - 0.25 u 3 -

0.5 u 4

+ 0.25 u 6 = - 1

Solving the equation simultaneously:

u 3 = -2.710 u 4 = 1.032 u 5 = -1.935 u 6 = -0.258

Validity if the results: It is obvious that the results obtained from the example are meaningless. In the perfect material, such as assumed, the horizontal displacement of point b and c must both be inwards and equal in magnitude. Similarity, the vertical displacement of B and C must both be outwards and equal in magnitude. The reason is that the grid chosen was not sufficiently refined and also not symmetrical.