MS5019 FEM

2.1. Definition of the Stiffness Matrix

Stiffness matrix, k is a matrix such that f = kd where k relates local - coordinate ( x, y, z ) nodal displacement to local force f of a single element.

The bold notation denotes a matrix or a vector, and the symbol denotes quantities referred to a local - coordinate system set up to be convenient for the element. For a continuous medium or structure with n element, a stiffness matrix K relates global - coordinate ( x, y, z ) nodal displacement d to global force F of the whole medium or structure.
MS5019 FEM 2

2.2. Stiffness Matrix for a Spring Element

Consider a linear spring element shown in Figure 2-1. References point 1 and 2 (nodes) are located at the ends of the element. We want to develop a relationship between nodal force and nodal displacement for a spring element. This relationship will be the stiffness matrix as follows.

f1x k11 = f 2 x k 21

k12 d1x k 22 d 2 x
k

(2.2.1)

where the element kij of the k matrix in Eq. (2.1) are to be determined.
1
f1x , d1x

2
f2x , d2x

Figure 2-1 Linear spring element with positive nodal displacement and force conventions
MS5019 FEM 3

Step 1 Select Element Type

Consider a linear spring subjected to resulting nodal forces T directed a long the spring axial direction x as shown in Figure 2-2, so as to be in equilibrium. The spring is represented by labeling nodes at each end. The original distance between nodes before deformation is denoted by L.

T
d2x

d1x
1

L
k
2

Figure 2-2 Linear spring subjected to tensile forces

MS5019 FEM 4

Step 2 Select a Displacement Functions A displacement function is assumed. Here a linear displacement varation along the axis of the spring is assumed because a linear function with specified endpoints has a unique path. Therefore

u = a1 + a2 x

(2.2.2)

In general, the total number of coefficients a is equal to the total number of degrees of freedom (dof) associated with the element. Here the total number of dof is two an axial displacement at each nodes of the element. In matrix form Eq. (2.2.2) becomes
a u = [1 x ] 1 a 2
MS5019 FEM

( 2.2.3)
5

We now want to express as a function of the nodal displacement. We achieve this by evaluating at each node and solving for a1 and a2 from Eq. (2.2.2) as follows

u (0) = d1x = a1 u ( L) = d 2 x = a2 L + d1x or, solving Eq. (2.5) for a2 , d d a2 = 2 x 1 x L

d d u = 2 x 1x x + d1x L
MS5019 FEM

(2.2.4) (2.2.5)

(2.2.6)

On substituting Eqs. (2.2.4) and (2.2.6) into Eq. (2.2.2), we have

(2.2.7)

In matrix form, we express Eq. (2.2.7) as

x u = 1 L or Here u = [N1 x d1x L d 2 x

( 2.2.8)

d N 2 ] 1x d 2 x x x N1 = 1 and N 2 = L L

( 2.2.9) (2.2.10)

are called the shape function because the Nis express the shape of the assumed displacement function over the domain of the element when the ith element degree of freedom has unit value and all other dof are zero.
MS5019 FEM 7

In the case, N1 dan N2 are linear functions that have the properties that N1 = 1 at node 1 and N1 = 0 at node 2, whereas N2 = 1 at node 2 and N2 = 0 at node 1. Also, N1 + N2 = 1 for any axial coordinate along the bar. In addition, the Nis are often called interpolation function because we are interpolating to find the value of a function between given nodal value. The interpolation function may be different from the actual function except at the endpoints or nodes where the interpolation function and actual function must be equal to specified nodal values.

MS5019 FEM

Step 3 Define the Strain-displacement and Stress-strain Relationships The tensile forces T produce a total elongation (deformation) of the spring. For the linear spring, T and are related through Hookes law by.

T = k (2.2.11) where, because is the deformation of the spring, we have = u ( L ) u ( 0) (2.2.12) Making use of Eq. (2.2.7), Eq. (2.2.11) becomes = d d (2.2.13)
2x 1x

MS5019 FEM

Step 4 Derive the Element Matrix and Equations

By the sign convention for nodal forces, we have f = T and f = T (2.2.14)

1x 2x

Using Eqs. (2.2.11), (2.2.13), and (2.2.14), we have T = f1x = k (d 2 x d1x ) (2.2.15) T = f 2 x = k (d 2 x d1x ) or rewriting Eq. (2.2.15), we obtain f = k (d d )
1x 1x 2x

f 2 x = k (d 2 x d1x )

(2.2.16)

MS5019 FEM

10

Now expressing Eq. (2.2.16) in a single matrix equations yields f1x k k d1x (2.2.17) = f 2 x k k d 2 x This relationship holds for the spring along the x axis. For our basic definition of a stiffness matrix and use of Eq. (2.2.1) applied to Eq. (2.2.17), we obtain k k= k k k (2.2.18)

as the stiffness matrix for a linear spring element. We observe that k is a symmetric square matrix.

MS5019 FEM

11

Step 5 Assemble the Element Equations to Obtain the Total/ Global Equations and Introduce Boundary Equations The total stiffness matrix and total force vector are assembled using nodal force equilibrium equations, force/deformation and compatibility equations from Section 2.2, and the direct stiffness method described in Section 2.4. This step applies for structures composed of more than one element such that

K = [K ] = k ( e )
e =1

and

F = {F } = f ( e )
e =1

(2.2.19)

where k and f are now element stiffness and force matrices expressed in a global reference frame.
The sign used in this context means that all element matrices must be assembled properly according to the direct stiffness method described in Section 2.4.
MS5019 FEM 12

Step 6 Solve for the Nodal Displacements The displacements are then determined by imposing boundary conditions and solving a system of equations, F = K d, simultaneously.

Step 7 Solve for the Element Forces The element forces are determined by back-substitution, applied to each element, into equation similar to Eq. (2.2.16).

Step 8 Interpret the results Finally, the element forces resulted can be analysed to check if all element forces are below its allowable values.
MS5019 FEM 13

2.3. Example of a Spring Assemblage

We will consider the specific example of the two-spring assemblage shown in Figure 2-3. Here we fix node 1, and apply axial forces at node 2 and node 3. The stiffness of spring elements 1 and 2 are k1 and k2 respectively. The x axis is the global axis of the assemblage. In this case, the local axis of each element coincides with the global axis of the assemblage.
1 1 2

x
F2 x

k1

F3 x
Figure 2-3 Two-spring assemblage

MS5019 FEM

14

Using Eq. (2.17) the element stiffness matrix for each element cam be expressed for element 1
f1x k1 = f 3 x k1 and for element 2, f3 x k2 = f 2 x k 2 k1 d1x k1 d 3 x k 2 d 3 x k 2 d 2 x (2.3.1)

(2.3.2)

Furthermore, element 1 and 2 must remain connected at common node 3 throught the displacement. This is called the continuity or compatibility requirement. The compatibility requirement yields (superscript refers to the element number.)

d 3(1) = d 3( 2 ) = d 3 x x x
MS5019 FEM

(2.3.3)
15

Based on the sign convention for element nodal forces given in Figure 2-1, we can write nodal equilibrium equations at node 3, 2, and 1 as
1 F3 x = f 3(x ) + f 3(x2 )

(2.3.4) (2.3.5) (2.3.6)

F2 x = f 2(x2 ) F1x = f1(x1)

Where F1x results from the reaction at the fixed support. To further clarify the resulting Eqs. (2.3.4) (2.3.6), free-body diagrams of each element and node (using the established sign conventions for element nodal forces) are shown in Figure 2-4.
1
( f1x1)

(1 f 3x )

( f 3x2 )

f 2(x2)

F1x
MS5019 FEM

F3 x
Figure 2-4 Nodal forces consistent with element forces sign convention

F2 x
16

Application of Newtons third law for each node and element gives
F3 x = k1d1x + k1d 3 x + k 2 d 3 x k 2 d 2 x F2 x = k 2 d 3 x + k 2 d 2 x F1x = k1d1x k1d 3 x ( 2 .3 .7 )

In matrix form, Eqs. (2.3.7) are expressed by

F3 x k1 + k 2 F2 x = k 2 F k 1 1x k2 k2 0 k1 d 3 x 0 d 2 x k1 d1x (2.3.8)

Or rearranging, Eqs. (2.3.8) in numerically increasing order of the nodal dof, we have
F1x k1 F2 x = 0 F k 3x 1
MS5019 FEM

0 k2 k2

k1 d1x k 2 d 2 x k1 + k 2 d 3 x

(2.3.9)

17

Eq. (2.3.9) is now written as the single matrix equation

F = Kd

(2.3.10) d1x d 2 x is called d 3x

F1x where F2 x is called the global nodal force vector , F 3x the global nodal displacement vector , and 0 k1 k1 0 k2 K= k2 k1 k 2 k1 + k 2 is called the total or global stiffness matrix.

(2.3.11)

MS5019 FEM

18

2.4.Assembling the Total Stiffmess Matrix by Superposition (Direct Stiffness Method)

This method is based on proper superposition of the individual element matrices making up a structure. Referring to the-spring assemblage of Section 2.3, the element stiffness matrices are given in Eqs. (2.3.1) and (2.3.2) as d1x d3 x d3x d2x k1 k2 k k (2.4.1) k (1) = 1 k ( 2) = 2 k1 k1 k 2 k 2 Here the d ix ' s written above the column in the k ' s indicate the degrees of freedom associated with each element.

MS5019 FEM

19

To superpose the element matrices, all of them must be expended to the order of the total structure (spring assemblage) stiffness matrix so that each element matrix is associated with all the dof of the structure. To expand each element stiffness matrix to the order of total stiffness matrix, we simply add rows and columns of zero for the those displacement not associated with that particular element. The expanded form of each element equation can be expressed as

for element 1
(1) (1) 1 0 1 d1x f1x ( 1 k1 0 0 0 d 21x) = f 2(x ) 1 0 1 d 3(1) f 3(x1) x

(2.4.2)

MS5019 FEM

20

for element 1
( 2) (2) 0 0 0 d1x f1x ( (2.4.3) k 2 0 1 1 d 22 ) = f 2(x2 ) x ( 2) (2) 0 1 1 d 3 x f 3 x Now considering force equilibrium at each node results in

f1(x1) 0 F1x ( 2) 0 + f 2 x = F2 x f (1) f ( 2 ) F 3x 3x 3x Using Eqs. (2.4.2) and (2.4.3) in Eq. (2.4.4), we obtain
(1) (2) 0 0 0 d1x F1x 1 0 1 d1x ( ( k1 0 0 0 d 21x) + k 2 0 1 1 d 22 ) = F2 x x d (1) 0 1 1 d 3( 2 ) F3 x 1 0 1 3 x x
MS5019 FEM

(2.4.4)

(2.4.5)
21

Simplifying Eq. (2.4.5) results in k1 d1x F1x 0 k1 0 k 2 d 2 x = F2 x k2 k1 k 2 k1 + k 2 d 3 x F3 x

(2.4.6)

Here the superscripts indicating the element numbers associated with 1 the nodal displacements have been dropped because d1(x) is really d1x , ( 2 d 22 ) is really d 2 x , and, by Eq. (2.3.3) d 3(1) = d 3( x ) = d 3 x , the node 3 x x displacement of the total assemblage. Equation (2.4.6), obtained through superposition, is identical to Eq. (2.3.9).
The method of directly assembling individual element stiffness matrices to form the total structure stiffness matrix and the the total set of stiffness equations is called direct stiffness method. It is the most important step in the FEM.
MS5019 FEM 22

2.5.Boundary Conditions
We must specify boundary (or support) conditions for structure models such as the spring assemblage of Figure 2-4, or K will be singular; that is the determinant of K will be zero and, therefore, its inverse will not exist. Without specifying adequate kinematic constraints or support conditions, the structure will be free to move as a rigid body. Boundary conditions (BC) are of two general type: 1. Homogeneous BC the most common occur at locations that are completely prevented from movement, 2. Non- homogeneous BC occur where finite non-zero values of displacement are specified, such as the settlement of a support.

MS5019 FEM

23

In general, specified support conditions are treated mathematically by partitioning the global equilibrium equations as follows:

Simplifying Eq. (2.4.5) results in K11 K12 d1 F1 (2.5.1) K = 21 K 22 d 2 F2 where we let d1 be the unconstrained or free displacement and d 2 be the the specified displacements nodes. From. (2.5.1), we have K11d1 = F1 K12d 2 (2.5.2) and F2 = K 21d1 + K 22d 2 (2.5.3) where F1 are the known nodal forces and F2 are the unknown nodal forces at the specified displacement nodes. F2 is found from Eq. (2.5.3) after de termining d1 from Eq. (2.5.2). In the Eq. (2.5.2), we assume that K11 is no longer singular, thus allowing for the determination of d1.
MS5019 FEM 24

Homogenoue BC
To illustrate the two general types of BC, let us consider Eq. (2.4.6), derived for the spring assemblage of Figure 2-4. We will first consider the case of homogenous BC. Hence, all BC are such that the displacements are zero at certain nodes.Here we have d1x = 0 because node 1 is fixed. Therefore, Eq. (2.4.6) can be written as. k1 0 F1x 0 k1 0 k2 d2 x = F2 x (2.5.4) k2 k1 k2 k1 + k2 d3x F3x Equation(2.5.4),writtenin expandedfrombecomes k1 (0) + (0)d2 x k1d3x = F1x 0(0) + k2d2 x k2d3x = F2 x (2.5.5) k1 (0) k2d2 x + (k1 + k2 )d3x = F3x
MS5019 FEM 25

Consider the second and third of Eqs. (2.5.5), written in matrix form, we have k 2 d 2 x F2 x k2 (2.5.6) k k + k d = F 2 1 2 3x 3x
We have now effectively partitioned off the first column and row of K and the first row of d and F to arrive at Eq. (2.5.6). For homogenous BC, Eq. (2.5.6) could have been obtained directly by deleting the row and column of Eq. (2.5.4) corresponding to the zerodisplacement degree of freedom. Here row 1 and column 1 are deleted because d1x = 0. However, F1x is not necessary zero and must be determined as follows.

MS5019 FEM

26

After solving Eq. (2.5.6) for d 2 x and d 3 x , we have 1 1 1 + k 2 F2 x k 2 k1 k1 F2 x d 2 x k 2 (2.5.7) = = 1 F3 x d 3 x k 2 k1 + k 2 F3 x 1 k1 k1 Using Eq. (2.5.7) in the first of Eqs. (2.5.5), we obtain the reaction F1x as F1x = k1d 3 x (2.5.8) We can express the unknown nodal force at node 1 (the reaction force) in terms of applied nodal forces F2 x and F3 x by using Eq. (2.5.7) for d 3 x sub stituted into Eq. (2.5.8). The results is F1x = F2 x F3 x (2.5.9)
1

MS5019 FEM

27

Note for homogenoue BC. For all homogenous BC, we can delete the rows and culumns corresponding to the zero-displacement dof from the original set of equatons and then solve for the unknown displacements. This procedure is useful for hand calculations. (More practical, computer-assisted scheme for solving the system of simultaneous equations will be discussed in another chapter.) Students are encouraged to review the numerical method to solve the system of simultaneous equations, i.e. Gauss or Gauss-Seidel method or other.

MS5019 FEM

28

Non-homogeneous BC
Next, we consider the case of non-homogenous BC. Hence, some of the specified displacements are non-zero. For simplicitys sake, let d1x = , where is a known displacement, in the Eq. (2.4.6). We now have

0 k1 F1x k1 0 (2.5.10) k2 k 2 d 2 x = F2 x d F k1 k 2 k1 + k 2 3 x 3 x Equation (2.5.10) written in expanded form becomes k1 + 0d 2 x k1d 3 x = F1x 0 + k 2 d 2 x k 2 d 3 x = F2 x k1 k 2 d 2 x + (k1 + k 2 )d 3 x = F3 x
MS5019 FEM 29

(2.5.11)

Considering the second and third of Eqs. (2.5.11) because they have known right - side nodal forces F2 x and F3 x 0 + k 2 d 2 x k 2 d 3 x = F2 x (2.5.12) k1 k 2 d 2 x + (k1 + k 2 )d 3 x = F3 x Transforming the known to the right side of Eqs. (2.5.12) results in k 2 d 2 x k 2 d 3 x = F2 x (2.5.13) k 2 d 2 x + (k1 + k 2 )d 3 x = k1 + F3 x Rewritting Eqs. (2.5.13) in matrix form, we have k2 k 2 k 2 d 2 x F2 x = k1 + k 2 d 3 x k1 + F3 x (2.5.14)

MS5019 FEM

30

Note for homogenoue BC. When dealing with non-homogeneous BC, one cannot initially delete row 1 and column 1 of Eq. (2.5.1)), corresponding to the nonhomogeneous BC, as indicated by the resulting Eq. (2.5.14). Had we done so, the k1 term in Eq. (2.5.14) would have been neglected, resulting in an error in the solution for the displacement. For non-homogeneous BC, we must, in general, transform the terms associated with the known displacement to the right-side force matrix before solving for the unknown nodal displacements. This was illustrated by transforming the k1 term of the second row of Eqs. (2.5.12) to the right-side of the second row of Eqs. (2.5.13). Finally, we could solve for the displacement in Eq. (2.5.14) in a manner similar to that used to solve Eq. (2.5.6).
MS5019 FEM 31

Some properties of the stiffness matrix (that are also applicable to the generalization of the finite element method). 1. K is symetric, as is each of the element stiffness matrices. 2. K is singular and thus no inverse exists until sufficient BC are imposed. 3. The main diagonal terms of K are always positive. Otherwise, a positive nodal force Fi could produce a negative displacement di a behavior contrary to the physical behavior of any actual structure.

Example 2.1

MS5019 FEM

32

2.6.Potential Energy Approach

One of the alternative methods often used to derive the element equations and the stiffness matrix for an element is based on the principle of minimum potential energy (POMPE). Adantages of the POMPE:
1. More general than previous method (equilibrium equations of nodal and element). 2. More adaptable for the determination of element equations for complicated elements (those with large numbers of dof) such as plane stress/strain element, the axisymetric stress element, the bending plate element, and the 3D solid stress element). 3. Only applicable for elastic materials (principle of virtual work is applicable for any material behavior). 4. Has lack of physical insight.
MS5019 FEM 33

The total potential energy (TPE), p of a structure is expressed in term of displacements. In the FE formulation, these will generally be nodal displacements such that p = p (d1, d2, d3, ..., dn). When p is minimized with respect to these displacements, equilibrium equations result. For the spring element, we will show that the same nodal equilibrium equations result as previously derived in Section 2.2. The POMPE is stated as follows: Of all the displacements that satisfy the given boundary conditions of a structure, those that safisfy the equations of equilibrium are distinguishable by the stationary value of the potential energy. If the stationary value is a minimum, the equilibrium state is stable. To explain this principle, we must first explain the concepts of PE and of a stationary value of a function.
MS5019 FEM 34

Total potential energy (PTE) is defined as the sum of the internal strain energy U and the potential energy of the external force ; thus is
Strain energy PE of external force

p =U +

(2.6.1)

Strain energy is the capacity of internal forces (or stresses) to do work through deformation (strains) in the structure. Potential energy (PE) of the external force , is the capacity of external forces such body force, surface traction forces, and applied nodal forces to do work through deformation of the structure.

MS5019 FEM

35

Recall that a linear spring has force related to deformation by F = kx, where k is the spring constant and x is the deformation of the spring (Figure 2-5). F

k
k x
Figure 2-5 Force-deformation curve for linear spring

F
x

The differential internal work (or strain energy) dU in the spring is the internal force multiplied by the change in displacement through which the force moves, given by.

dU = F dx
MS5019 FEM

(2.6.2)
36

Now we express F as (2.6.3) F = kx Using Eq. (2.6.3) in Eq. (2.6.2), the differential strain energy becomes dU = kx dx (2.6.4) The total strain energy is then given by U = kx dx
0 x

( 2.6.5)

Upon explicit integration of Eq. (2.6.5), we obtain U = 1 kx 2 (2.6.6) 2 Using Eq. (2.6.3) in Eq. (2.6.6), we have U = 1 (kx) x = 1 Fx ( 2.6.7) 2 2 Equation (2.6.7) indicates that strain energy is the area under force - deformation curve.
MS5019 FEM 37

The PE of external force, being opposite in sign from the external work expression because the PE of external force is lost when the work is done by the external force, is given by
= Fx (2.6.8) Therefore substituting Eqs. (2.6.6) and (2.6.8) into (2.6.1), the TPE becomes

p = 1 kx 2 Fx 2

(2.6.9)

To apply the POMPE that is, to minimize of p we take the variation of p defined in general as

p =
MS5019 FEM

p d1

d1 +

p d 2

d 2 + L +

p d n

d n

(2.6.10)
38

The principles states that equilibrium exists when the di define a structure state such that p = 0 for arbitrary admissible variations di from the equilibrium state. An admissible variation is one in which the displacement field still satisfies the BC and inter-element continuity.

To satisfy p = 0 (for d i 0), all coefficients associated with the d i must be zero independently. Thus p p = 0 (i = 1,2,3, L, n) or =0 (2.6.11) d i {d } where n equations must be solved for the n values of di that define static equilibrium state of the structure.

MS5019 FEM

39

Equation (2.6.11) shows that for our purposes throughout this text, we can interpret the variation of p (p) as a compact notation equivalent to differentiation of p with respect to the unknown nodal displacements for which p is expressed. We now derive the spring element equations and stiffness matrix using the POMPE by analyzing a linear single-dof spring as shown in Figure 2-6.

f1x

f2x

Figure 2-6 Linear spring subjected to nodal forces

MS5019 FEM

40

Using Eq. (2.6.9), the TPE becomes 2 p = 1 k d 2 x d1x f1x d1x f2 x d 2 x 2

( (

(2.6.12)

where d 2 x d1x is the deformation of the spring in Eq. (2.6.9). Simplifying Eq. (2.6.12), we obtain p = 1 k d 22x 2d 2 x d1x + d12x f1x d1x f2 x d 2 x (2.6.13) 2

Minimization of p with respect to each nodal displacement requires that p 1 = 2 k 2d 2 x + 2d1x f1x = 0 d

1x

p 1 = 2 k 2d 2 x 2d1x f 2 x = 0 d
2x
MS5019 FEM 41

(2.6.14)

Simplifying Eq. (2.6.14), we have k (d + d ) = f

2x 1x 1x (2.6.15) k (d 2 x d1x ) = f 2 x In matrix form, we express Eq. (2.6.15) as k k d1x f1x (2.6.16) k k = d 2 x f 2 x Since f = k d, we have the stiffness matrix for the spring element obtained from Eq. (2.6.16) : k k (2.6.17) k= k k As expected, Eq. (2.6.17) is identical to the stiffness matrix obtained in Section 2.2, Eq. (2.2.18).
42

MS5019 FEM

We have developed the FE spring element equations by minimizing the TPE with respect to the nodal displacements. Now, we show that the TPE of an entire structure (here an assemblage of spring elements) can be minimized with respect to each nodal dof and this minimization results in the same FE equations used for the solution as those obtained by the direct stiffness method explained in Section 2.4.

Example 2.5

MS5019 FEM

43

Reference: 1. Logan, D.L., 1992, A First Course in the Finite Element Method, PWS-KENT Publishing Co., Boston. 2. Imbert, J.F.,1984, Analyse des Structures par Elements Finis, 2nd Ed., Cepadues. 3. Zienkiewics, O.C., 1977, The Finite Eelement Method, 3rd ed., McGraw-Hill, London.

MS5019 FEM

44