INTRODUCTION TO THE FINITE ELEMENT METHOD Evgeny Barkanov ... Free changes of are called the variations. The mathematical condition of the minimum of ...

2 downloads 114 Views 1MB Size

Riga, 2001

Preface

Today the finite element method (FEM) is considered as one of the well established and convenient technique for the computer solution of complex problems in different fields of engineering: civil engineering, mechanical engineering, nuclear engineering, biomedical engineering, hydrodynamics, heat conduction, geo-mechanics, etc. From other side, FEM can be examined as a powerful tool for the approximate solution of differential equations describing different physical processes. The success of FEM is based largely on the basic finite element procedures used: the formulation of the problem in variational form, the finite element dicretization of this formulation and the effective solution of the resulting finite element equations. These basic steps are the same whichever problem is considered and together with the use of the digital computer present a quite natural approach to engineering analysis. The objective of this course is to present briefly each of the above aspects of the finite element analysis and thus to provide a basis for the understanding of the complete solution process. According to three basic areas in which knowledge is required, the course is divided into three parts. The first part of the course comprises the formulation of FEM and the numerical procedures used to evaluate the element matrices and the matrices of the complete element assemblage. In the second part, methods for the efficient solution of the finite element equilibrium equations in static and dynamic analyses will be discussed. In the third part of the course, some modelling aspects and general features of some Finite Element Programs (ANSYS, NISA, LS-DYNA) will be briefly examined. To acquaint more closely with the finite element method, some excellent books, like [1-4], can be used.

Evgeny Barkanov Riga, 2001

2

Contents

PREFACE……………………………………………………………………………….…2 PART I

THE FINITE ELEMENT METHOD……………………………………..…5

Chapter 1 Introduction………………………………………………………………...…5 1.1 Historical background………………………………………………………………...5 1.2 Comparison of FEM with other methods………………………………………….....5 1.3 Problem statement on the example of “shaft under tensile load”…………………….6 1.4 Variational formulation of the problem……………………………………………....9 1.5 Ritz method………………………………………………………………………….10 1.6 Solution of differential equation (analytical solution)………………………………12 1.7 FEM……………………………………………………………………………...….13 Chapter 2 Finite element of bending beam…………………………………………….20 Chapter 3 Quadrilateral finite element under plane stress………………………...…23 PART II SOLUTION OF FINITE ELEMENT EQUILIBRIUM EQUATIONS….30 Chapter 4 Solution of equilibrium equations in static analysis……………………….30 4.1 Introduction…………...……………………………………………………………..30 4.2 Gaussian elimination method……………………………………………………..…31 4.3 Generalisation of Gauss method………………………………………………….…31 4.4 Simple vector iterations…………………………………………………….……….33 4.5 Introduction to nonlinear analyses……………………………………………….….34 4.6 Convergence criteria………………………………………………………………...37 Chapter 5 Solution of eigenproblems…………………………………………………...39 5.1 Introduction…………………………………………………………………………39 5.2 Transformation methods…………………………………………………………….40 5.3 Jacobi method……………………………………………………………………….41 5.4 Vector iteration methods…………………………………………………………….42 5.5 Subspace iteration method…………………………………………………………..43 Chapter 6 Solution of equilibrium equations in dynamic analysis…………………...45 6.1 Introduction………………………………………………………………………….45 6.2 Direct integration methods…………………………………………………………..45 6.3 The Newmark method………………………………………………………………46 6.4 Mode superposition………………………………………………………………….47 6.5 Change of basis to modal generalised displacements……………………………….48 6.6 Analysis with damping neglected…………………………………………………...49 6.7 Analysis with damping included…………………………………………………….50

3

PART III EMPLOYMENT OF THE FINITE ELEMENT METHOD……………...53 Chapter 7 Some modelling considerations…………………………………………..…53 7.1 Introduction………………………………………………………………………….53 7.2 Type of elements…………………………………………………………………….53 7.3 Size of elements……………………………………………………………………..55 7.4 Location of nodes…………………………………………………………………....56 7.5 Number of elements…………………………………………………………………56 7.6 Simplifications afforded by the physical configuration of the body………………..58 7.7 Finite representation of infinite body………………………………………………..58 7.8 Node numbering scheme……………………………………………………………59 7.9 Automatic mesh generation…………………………………………………………59 Chapter 8 Finite element program packages…………………………………………..60 8.1 Introduction………………………………………………………………………….60 8.2 Build the model……………………………………………………………………...60 8.3 Apply loads and obtain the solution………………………………………………...61 8.4 Review the results…………………………………………………………………...62 LITERATURE……………………………………………………………………………63 APPENDIX

A typical ANSYS static analysis……………………………………...64

4

PART I

THE FINITE ELEMENT METHOD

Chapter 1

Introduction

1.1

Historical background

In 1909 Ritz developed an effective method [5] for the approximate solution of problems in the mechanics of deformable solids. It includes an approximation of energy functional by the known functions with unknown coefficients. Minimisation of functional in relation to each unknown leads to the system of equations from which the unknown coefficients may be determined. One from the main restrictions in the Ritz method is that functions used should satisfy to the boundary conditions of the problem. In 1943 Courant considerably increased possibilities of the Ritz method by introduction of the special linear functions defined over triangular regions and applied the method for the solution of torsion problems [6]. As unknowns, the values of functions in the node points of triangular regions were chosen. Thus, the main restriction of the Ritz functions – a satisfaction to the boundary conditions was eliminated. The Ritz method together with the Courant modification is similar with FEM proposed independently by Clough many years later introducing for the first time in 1960 the term “finite element” in the paper “The finite element method in plane stress analysis” [7]. The main reason of wide spreading of FEM in 1960 is the possibility to use computers for the big volume of computations required by FEM. However, Courant did not have such possibility in 1943. An important contribution was brought into FEM development by the papers of Argyris [8], Turner [9], Martin [9], Hrennikov [10] and many others. The first book on FEM, which can be examined as textbook, was published in 1967 by Zienkiewicz and Cheung [11] and called “The finite element method in structural and continuum mechanics”. This book presents the broad interpretation of the method and its applicability to any general field problems. Although the method has been extensively used previously in the field of structural mechanics, it has been successfully applied now for the solution of several other types of engineering problems like heat conduction, fluid dynamics, electric and magnetic fields, and others. 1.2

Comparison of FEM with other methods

The common methods available for the solution of general field problems, like elasticity, fluid flow, heat transfer problems, etc., can be classified as presented in Fig. 1.1. Below FEM will be compared with analytical solution of differential equation and Ritz method considering the shaft under tensile load (Fig. 1.2).

5

Methods

Analytical

Numerical

Approximate

Exact (e.g. separation of variables and Laplace transformation methods)

Numerical solution

FEM

(e.g. Rayleigh-Ritz and Galerkin methods)

Numerical integration

Finite differences

Fig. 1.1 Classification of common methods. 1.3

Problem statement on the example of “shaft under tensile load”

The main task of the course “Strength of Materials” is determination of dimensions of a shaft cross section under known external loads. Applying the general plan for the solution of problems in the field of mechanics of deformable solids, tree group of equations should be written: 1) equilibrium equations (statics) The equilibrium equation for the separate element with the length dx has the following form

∑X

=0

or

− σF + (σ + dσ ) F + qdx = 0

After some transformations we have dσ F +q=0 dx Taking into account that σ = εE = d 2u dx 2

du E we obtain the static equilibrium equation dx

EF + q = 0

2) geometric equations ε=

du dx

3) physical equations σ = εE From this system of equations it is possible to determine all necessary values.

6

σ dx

dx

q q

F

σ+dσ x

x z

y

Fig. 1.2 Shaft under tensile load. Another approach for the solution of the problem examined exists also. This is utilisation of the principle of “minimum of the potential energy” which means: a system is in the state of equilibrium only in the case when it potential energy is minimal. Correctness of this principle may be observed on the following simple examples: - a ball is in the state of equilibrium only in the lower point of surface (Fig. 1.3), - a water on the rough surface takes the equilibrium state in the lower position, - a student tries to take examination with the minimum expenditures of labour. From the condition that the potential energy takes the minimum, it is possible to determine the unknown values. The general algorithm of solution in this case is following: 1) an expression for the potential energy of elastic system under external loads is written, 2) conditions of minimum of the potential energy are written, 3) unknown values are determined from the condition of minimum, 4) a strength problem is solved.

Π

P

R

P Π min ∆

R

Fig. 1.3 Principle of “minimum of the potential energy”.

7

∆

Complete potential energy of the deformable system consists from the strain energy U stored in the system and energy W lost by the external forces (Fig. 1.4). That is why the work of the external forces W is negative value Π = U −W Since the tension of a shaft is examined, U is the potential energy of tension. Then for the tension we have Π = U −W =

1 1 P∆ − P∆ = − P∆ 2 2

The force loses energy − P∆ , but the shaft acquires the tensile energy 1 2 ( P∆ ) . The second part goes on overcoming the friction forces, internal heat, changes into kinetic energy, etc. After removal of load, the system can gives back only the energy equal to the potential energy of tension U=

1 P∆ 2

P P

Π -

U +

∆ ∆

W -P ∆ H

P

H-∆

P(H - ∆ ) - PH = -P∆ = W final initial

Fig. 1.4 Energy balance. Π

δu(x)

u(x) Π(u+δu)

Π(u)

δΠ

Fig. 1.5 Variational formulation.

8

1.4

Variational formulation of the problem

A numerical value of the potential energy of tension Π = U − W is dependent from the function u( x ) to be used. Because Π is a functional, since a functional is a value dependent from the choice of function. This can be explained by the help of Fig. 1.5. In the lower point, an infinitesimal change of the function u( x ) equalled to δu( x ) will not give an increase of the functional δΠ . In the point of minimum: δΠ = 0 . Free changes of δu , δΠ are called the variations. The mathematical condition of the minimum of potential energy can be written as δΠ = 0 . How it can be seen, variation in the case of functional investigation has the same meaning as differential in the case of function investigation. Let’s investigate the functional Π of a tensile shaft under distributed load q 2

l 1l du Π ( u( x )) = U − W = ∫ EF dx − ∫ qudx 2 0 dx 0 2

l

l 2 2 l 2 2 2 σ F ε E F 1l 1l N2 du 2 dx = ∫ EFε dx = ∫ EF dx dx = ∫ dx = ∫ U =∫ 2 EF 2 EF 2 EF 20 2 0 dx 0 0 0

Let’s determine the variation δΠ as a difference of two values – the potential energy with and without increment δu δΠ = Π ( u + δu ) − Π ( u ) =

l 1l du du EF 2 δ dx − ∫ qδudx = 2 0∫ dx dx 0

l l du l l d 2 u l du dδu = ∫ EF dx − ∫ qδudx =EF δu o − ∫ δu 2 dx − ∫ qδudx = 0 dx dx dx dx 0 0 0 0

du ; s = δu dx Integration by parts: t= l

l

l

∫ s' tdx = st 0 − ∫ st' dx

0

0

The potential energy will has the minimum value, if δΠ = 0 , or by other words, if all items equal to zero in the last expression. Boundary conditions for our problem are: 1) x = 0 : u = 0 du 2) x = l : =0 dx At all length l: δu ≠ 0 . Applying these boundary conditions to the variation of the functional Π , we obtain l

− EF ∫ δu 0

l l d 2u dx − ∫ qδudx = − ∫ δu EF + q dx = 0 dx 2 dx 2 0 0

d 2u

9

d 2u

+ q = 0 . Moreover, this condition presents the dx 2 static equilibrium equation. Expressions obtained show that the potential energy of system has the minimum, if: 1) the equilibrium equations will be realised 2) the boundary conditions will be realised The second boundary conditions, so called as natural boundary conditions for the functional Π , since they are obtained from the minimum of functional, realise automatically. But it is necessary to satisfy without fail to the first boundary conditions. Otherwise, these conditions are not taken into account anywhere. These boundary conditions are called principal. In the case of beam bending: - natural conditions are forces, - principal conditions are displacements. The problem of determination of u( x ) can be solved by two ways: 1) by solution of the differential equation, 2) by minimising the functional Π . Solving the problem by the approximate methods using computers, the second way is more suitable. This equation can be solved if EF

1.5

Ritz method

By the Ritz method it is possible to determine an approximate min Π . An unknown function of displacements u( x ) is found in the form u( x ) = ∑ a k ϕ k ( x ) k

where ak are coefficients to be determined, ϕ k ( x ) are coordinate functions given so that they satisfy to the principal boundary conditions. By insertion u( x ) into functional Π and then to integrate, it is possible the problem of the functional minimisation to come to the problem of determination the function minimum Π = Π ( a k ) from unknowns ak . To minimise the function of the potential energy obtained, it is necessary to equate to zero the derivatives on a k ∂Π ∂Π ∂Π = 0 , ... = 0, = 0, ∂a3 ∂a 2 ∂a1 After this operation, the system of algebraic equations is obtained and solved to find the unknowns ak . In the Ritz method, the choice of function u( x ) closed enough to a truth is a complicated problem requiring a good idea of the result expected.

10

Example 1.1 1) u = ax The principal boundary conditions: x = 0 , u = 0 . Π [u( x )] =

2

l 1l 1 l2 du 2 − = − EF dx qudx EFa l qa ∫ 2 0∫ 2 2 dx 0

q = const ;

du =a dx

∂Π ql 2 ql 2 ql = EFal − = 0 , hence a = = ∂a 2 2 EFl 2 EF ql du ql ql u= x; ε = = ; σ = εE = 2 EF dx 2 EF 2F u( l ) =

2 ql 2 l ql ; u = 2 EF 2 4 EF

2) u = a1 x + a2 x 2 The principal boundary conditions: x = 0 , u = 0 . Π [u( x )] = ... ∂Π ∂Π = 0; =0 ∂a1 ∂a 2 Then the system consisting from 2 linear algebraic equations is solved and unknown coefficients a1 and a 2 are determined. After this example, it is possible to write the general algorithm of the Ritz method: Presentation of Π . Determination of boundary conditions. Approximation for all construction. Integration of Π for all construction. Determination of the minimum of Π . Solution of the system of linear algebraic equations. 6) Calculation of displacements. 7) Calculation of stresses. For a beam with few areas, it is more easily to guess the deflection function for separate areas than for the whole of beam. Moreover, the function for one area will be a more simple than for the whole of beam. The idea of division of the investigated object is used in FEM. In the Ritz method, the accuracy can be increased choosing more terms of the approximated function, but in FEM – increasing the quantity of finite elements. To simplify the problem solution by computer, the finite elements and approximated functions are chosen the same. 1) 2) 3) 4) 5)

11

1.6

Solution of differential equation (analytical solution)

Analytical solution means determination of the displacement function u( x ) from the equilibrium equation. Example 1.2 EF

d 2u dx 2

+q=0

To solve the problem of u( x ) determination it is necessary to satisfy to the following boundary conditions: 1) x = 0 : u = 0 du 2) x = l : =0 dx d 2u dx 2

=−

q EF

After integration we obtain du q qx =− dx + C1 = − + C1 ∫ dx EF EF u=−

qx 2 + C1 x + C 2 2 EF

The constants are determined using the boundary conditions: qx 2 + C1 x + C 2 = 0 ; C 2 = 0 2 EF ql ql 2) − + C1 = 0 ; C1 = EF EF 1) −

Then we have qx 2 qlx u=− + 2 EF EF u( l ) = −

ql 2 ql 2 ql 2 + = 2 EF EF 2 EF

ql 2 ql 2 3ql 2 l u = − + = 8 EF 2 EF 8EF 2 du qx ql σ = εE = E=− + dx F F du qx ql ε= =− + dx EF EF

12

1.7

FEM

FEM was treated previously as a generalisation of the displacement method for shaft systems. For a computation of beams, plates, shells, etc. by FEM, a construction is presented in a view of element assembly. It is assumed that they are connected in a finite number of nodal points. Then it is considered that the nodal displacements determine the field of displacements of each finite element. That gives the possibility to use the principle of virtual displacements to write the equilibrium equations of element assembly so, as made for a calculation of shaft systems. Let’s have a look the finite element of tensile shaft (Fig. 1.6). The displacement function can be chosen in the following form u = C 0 + C1 x Using boundary conditions for a single finite element in his local coordinate system, we have u1 = C 0 u 2 = C0 + C1l Now our purpose to express coefficients through the nodal displacements of the finite element C 0 = u1 u − u1 C1 = 2 l Then the displacement function for a single finite element can be written in the following form u −u u = u1 + 2 1 x = u1 1 − l N1 ( x ) = 1 −

x x + u 2 = u1 N1( x ) + u 2 N 2 ( x ) l l

x x , N 2 ( x ) = - nodal functions. l l

u1 1 (i) q

l

L

q

2 (i+1) u2 x x

Fig. 1.6 Finite element of tensile shaft.

13

The potential energy of the finite element can be expressed as follows 2

l 1l du = ∫ EF dx − ∫ qudx = 2 0 dx 0

Π

e

=

l l l 1 EF u12 ∫ ( N1' ( x ))2 dx + 2u1u 2 ∫ N1' ( x )N 2 ' dx + u 22 ∫ ( N 2 ' ( x ))2 dx − 2 0 0 0 l

l

0

0

− ( qu1 ∫ N1( x )dx + qu 2 ∫ N 2 ( x )dx ) =

{

}

1 K11u12 + 2 K12 u1u 2 + K 22 u 22 − {u1 F1 + u 2 F2 } 2

du = u1 N1' ( x ) + u 2 N 2 ' ( x ) dx 1 N1' ( x ) = − l 1 N 2' ( x ) = l

l

K11 = EF ∫ ( N1' ( x ))2 dx = EF 0 l

1 l

2

EF l

l=

EF 11 K12 = EF ∫ N1' ( x )N 2 ' ( x )dx = EF − l = − l ll 0 l

K 22 = EF ∫ ( N 2 ' ( x ))2 dx = EF 0

1 l

2

l=

EF l

K ij - elements of the stiffness matrix.

l

l

F1 = q ∫ N1( x )dx = q ∫ 1 − 0 0 l

l

l 2 x x = q l − l = ql dx = q x − l 2l 2 2 0

x x2 F2 = q ∫ N 2 ( x )dx = q ∫ dx = q l 2l 0 0

l

= 0

ql 2

Fi - nodal forces. Now the potential energy of the finite element presents the function of his nodal displacements Π e = Π e ( u1 , u 2 )

14

Let’s rewrite the potential energy of finite element in the matrix form K 1 Π e = {u1 u 2 } 11 2 K 21

K12 u1 F1 1 eT e e eT e − {u1 u 2 } = d K d − d F K 22 u 2 F 2 2 1x2 2x2 2x1

1x2 2x1

u de = 1 u 2 ql F Fe = 1 = 2 F2 ql 2 K12 EF 1 − 1 K K e = 11 = l − 1 1 K 21 K 22 Then it is possible to determine the potential energy of structure consisting from the separate energies of the finite elements Π =

N

1

∑ Π e = 2 d T Kd − d T F

l =1

where K =

N

∑Ke

is the stiffness matrix of construction as a sum of stiffness matrices of

e =1

separate finite elements, d is the vector of nodal unknowns of construction, F is the vector of given external nodal forces. By this way, the potential energy of structure is expressed in a view of function dependent on unknown nodal displacements d . The condition of the functional minimum turns into condition the function minimum Π ( d ) ∂Π =0 ∂u i Solution of this system is unknown displacements d = K −1 F It is necessary to note that solving the present system of equations, it is necessary to take into account conditions of structure supports, that is to say the principal boundary conditions. After determination of the nodal displacements d , the internal forces and stresses σ are computed. Then they are used for a valuation of the structure’s strength.

15

Example 1.3 u1 = 0; F1 =

1

I

L

q

ql 2

u 2 ; F2 = F( 2I ) + F(1II) =

II 2

3 L = l - length of FE 2

u3 ; F3 = F( II2 ) =

ql ql + = ql 2 2

ql 2

x

Fig. 1.7 Finite element model of shaft under tensile load. The potential energy of shaft under tensile load can be expressed as follows 1 I 2 I I 2 Π = Π I + Π II = K11 u1 + 2 K12 u1u 2 + K 22 u 2 − u1 F1I + u 2 F2I + 2 1 II 2 II II 2 + K11 u 2 + 2 K12 u 2 u 3 + K 22 u 3 − u 2 F1II + u 3 F2II = 2 1 I II 2 II II 2 = K 22 + K11 u 2 + 2 K12 u 2 u 3 + K 22 u 3 − (u 2 F2 + u 3 F3 ) = 2 1 1 = K 22 u 22 + 2 K 23u 2 u 3 + K 33u 32 − (u 2 F2 + u 3 F3 ) = d T Kd − d T F 2 2

(

( {( {

)(

)(

)

)

)

}

}

1x2 2x2 2x1 1x2 2x1

I II K 22 = K 22 + K11 II K 23 = K12

- in the global coordinate system.

II K 33 = K 22

K K = 22 K 32 K=

EF l

K 23 ; K 33

2 − 1 − 1 1 ;

u d = 2; u 3 ql F = ql 2

F F = 2 F3

16

Now it is necessary to determine the nodal displacements of the structure using the principle of “minimum of the potential energy”.

Π(u2,u3)

min u2,u3

Fig. 1.8 Possible solution. These unknowns are determined from the following system of linear algebraic equations ∂Π ∂u = 0 2 ∂Π =0 ∂u 3 ∂Π ∂u = K 22 u 2 + K 23u 3 − F2 = 0 2 ∂Π = K 23u 2 + K 33u 3 − F3 = 0 ∂u 3 3ql 2 u = 2 2 EF 2 2 u = ql 3 EF

since l =

L , then we have 2

3qL2 u 2 = 8 EF 2 u = qL 3 2 EF

Stresses can be calculated by the following way u( x ) = u1 N1( x ) + u 2 N 2 ( x ) ε=

du = u1 N1′ ( x ) + u 2 N 2′ ( x ) dx

17

σ I = Eε = Eu1 N1' ( x ) + Eu 2 N 2 ' ( x ) = E

2 3qL2 3 qL = L 8EF 4 F

=0

N2( x ) =

x ; l

1 2 N 2' ( x ) = = l L

2 2 qL2 3 qL qL 1 qL 2 3qL σ II = EN1' ( x )u 2 + EN 2 ' ( x )u 3 = E − +E =− + = L 2 EF 4 F F 4 F L 8 EF

N1 ( x ) = 1 −

x ; l

1 2 N1' ( x ) = − = − l L

Example 1.1, 1.2, 1.3 Let’s compare the analytical solution with the solution obtained by FEM and Ritz methods.

ql F

σ

u

analytical solution

Ritz method FEM 3 ql 4F

analytical solution 2

ql 4 EF

1 ql 4F

2

3ql 8 EF

ql 2F

ql 2 2 EF

Fig. 1.9 Analytical, FEM and Ritz solutions.

18

FEM

Ritz method

After this example, it is possible to write the general algorithm of FEM: 1) Presentation of Π . 2) Determination of boundary conditions: δΠ = 0 . 3) Approximation for the finite element. 4) K e - integration (analytical or numerical). 5) Finite element meshing. 6) Building of K .

By computer

7) Determination of the minimum of Π . Solution: Kd = F, d = K −1F , where K is symmetrical and banded. 8) Output of displacements. 9) Computation of stresses. 10) Output of stresses. FEM Accuracy of FEM: ∆FEM approximate ≤ ∆exact , Π approximate ≥ Π exact .

19

Chapter 2

Finite element of bending beam

The functional Π of bending beam loaded by the concentrated forces Pi , bending moments M j and distributed load q k can be written in the following form dw j L 1L 2 Π = U − W = ∫ EJ ( w′′ ) dx − ∑ Pi wi − ∑ M j − qwdx 20 dx 0∫ i j

(1)

Then it is necessary to describe the boundary conditions. In our case, the principal boundary conditions are w,

dw dx

and the natural boundary conditions w′′ ( M = − EJw′′) , w′′′ (Q = − EJw′′′)

w

M

P q

x

dw dx

dx L w

wi ' wi i

w i l

w1 '

'

w2 ' w2

w1

i+1

w i+1

i+1

i 1

i+1

2 l

Fig. 2.1 Finite element of bending beam.

20

x

Besides, we have additional conditions - two principal boundary conditions which should be realised at each end of the finite element. These are conditions of joining of two neighbouring elements ′ = w(′ i +1 ) w( i ) = w( i +1 ) , w(i) To satisfy these four conditions, let’s choose the polynom with four coefficients as coordinate function w( x ) = a 0 + a1 x + a 2 x 2 + a3 x 3 In such view the coordinate function w( x ) does not satisfy to the boundary conditions yet. Therefore, let’s change it so, that coefficients a 0 , a1 , a 2 , a3 were expressed through unknowns in the nodal points of element ends - w1 , w1′ , w2 , w2′ , where 1 and 2 are the numbers of nodal points w1 = a0 (when x = 0 ) and etc. Then the system of equations is solved in relation to a0 , a1 , a 2 , a3 . Substituting these expressions into coordinate function and introducing the nodal functions N1( x ), N 2 ( x ), N 3 ( x ), N 4 ( x ) , we obtain w( x ) = w1 N1 ( x ) + w1′ N 2 ( x ) + w2 N 3 ( x ) + w2′ N 4 ( x )

N1 ( x ) = 1 − 3

x2 l2

N2( x ) = x − 2 x2

N3( x ) = 3

l2

N4( x ) = −

+2

(2)

x3 l3

x 2 x3 + l l2 −2

x3 l3

x 2 x3 + l l2

In such view the coordinate function w(x) satisfies to the principal boundary conditions. Then we substitute the expression (2) in (1) and obtain after integration the potential energy of the finite element Π e = Π e ( w1 , w1′ , w2 , w2′ ) Πe=

1 eT e e d K d − d eT F e 2 1x4 4x4 4x1 1x4 4x1

21

where K e is the stiffness matrix of the finite element of bending beam, d e is the vector of nodal unknowns of the finite element, F e is the vector of given external nodal efforts, when the external load is presented by the forces and moments in the nodal points.

12 6l 4l 2 EJ e K = l3

ql 2 − 12 6l F1 ql 2 w1 − 6l 2l 2 e M 1 12 e w1′ , d = , F = = ql 12 6l F2 w2 M 2 2 w2′ 4l 2 ql 2 − 12

Now it is possible to determine the potential energy of structure consisting from the separate energies of the finite elements Π = ∑Π e N

The complete potential energy is a function of unknowns – displacements and angles of rotations in the nodal points. To obtain the minimum of the potential energy, as in the Ritz method, we take derivatives on unknowns, equate to zero and obtain the system of algebraic equations for determination of unknown values. Assuming that a beam consists from one finite element ( Π = Π e ), the condition of minimum can be written as ∂Π e ∂Π e ∂Π e ∂Π e = 0, = 0, = 0, =0 ∂w1 ∂w1′ ∂w2 ∂w2′

22

Chapter 3

Quadrilateral finite element under plane stress

Since general relations of plane strain and plain stress differ only by the elastic constants, a solution of the plane problem in the theory of elasticity we examine on the base of plane stress. For the calculation of plates loaded in their plane, the functional of complete potential energy for the plane stress is written in the following form: Π = U −W =

1 ( σ x ε x + σ y ε y + τ xy γ xy )dΩ − ∫ ( p x u + p y v )dL 2 Ω∫ L

(1)

where σ x , σ y , τ xy are the normal and tangential stresses, ε x , ε y ,γ xy are the linear and angle strains, u , v are the linear displacements of the points on the middle plane of plate in relation to axes x and y, p x , p y are the vector components of external loading in relation to axes x and y, dΩ , dL are infinitely small element of two-dimensional area and outline. For the plane problem in the theory of elasticity we have σ x px u d = , F = , σ = σ y , ε = v py τ xy

ε x ε y γ xy

(2)

For the isotropic material, the general relations of the plane stress can be presented as E σ x 1 − υ 2 υE σ = σ y = 2 τ 1 − υ xy 0 ∂ εx ∂x ε = ε y = 0 γ ∂ xy ∂y

υE 1−υ E

2

1−υ 2 0

εx 0 ε y , E γ xy 2( 1 + υ ) 0

0 ∂ u ∂y v ∂ ∂x

(3)

(4)

or in the matrix form:

23

σ = Eε , ε = Dd

(5)

where E is the matrix of elasticity, D is the matrix of differentiation. Now the functional of complete potential energy of the plate loaded in it plane can be written in the compact form: Π =

1 T ε σ dΩ − ∫ F T d dL ∫ 2Ω L

(6)

For the building of stiffness matrix, it is necessary to set the displacement approximation for the finite element area and to connect it with the degrees of freedom. For an existence of the functional of complete potential energy, the approximation functions of displacements should contain terms not lower than first order. The linear polynomial from two variables contains three terms. To connect four nodes of the quadrilateral finite element with necessary quantity of constant coefficients of the approximation functions, the following form of these polynomials are taken u( x , y ) = a1 + a 2 x + a3 y + a 4 xy

(7)

v( x , y ) = a5 + a 6 x + a7 y + a8 xy This model corresponds to the linear distribution of displacements along coordinate axes. The number of linearly independent coefficients is twice more than the number of finite element nodes. On this reason, for each node it is possible to give two degrees of freedom. Thus, the finite element has eight degrees of freedom (Fig. 3.1). The vectors of nodal displacements and nodal reactions have the following form: R x1 u1 R v y1 1 Rx2 u 2 R v y2 2 d= , R= R u 3 x3 R y3 v3 u R 4 x4 v 4 R y 4 v1

v2 u 2

1

2

4

3

v4

v3

x

b

u1

(8)

u4

u3

a y

Fig. 3.1 Quadrilateral finite element.

24

The stiffness matrix K with dimension 8x8 connects these vectors by the following way R = Kd

(9)

Let’s express the linearly independent constant coefficients of approximation functions by the nodal displacements. For this purpose coordinates ( x1 , y1 ) for the first node are substituted to the expression (7) and we have u1 = a1 + a 2 x1 + a3 y1 + a 4 x1 y1 The same operation is repeated for the second and other nodes. After that we have the system of four linear algebraic equations in relation to the constant coefficients ai (i = 1,2,3,4) : u1 1 u 1 2 = u 3 1 u 4 1

x1 x2 x3 x4

y1 y2 y3 y4

x1 y1 a1 x 2 y 2 a 2 x3 y 3 a 3 x 4 y 4 a 4

or in the compact form: d = Ca

(10)

Solving the system (10), the constant coefficients ai (i = 1,2,3,4) are determined a = C −1d

(11)

Then approximation functions can be written in the form: u( x , y ) =

4

∑ di Ni

(12)

i =1

v( x , y ) =

8

∑ di Ni

i =5

where d i (i = 1,2 ,...,8) are the degrees of freedom of the finite element, N i (i = 1,2 ,...,8) are the nodal functions. As it is seen, only coefficients for the function u( x , y ) were determined, since the coefficients for the function v( x , y ) have the same form. In the detailed form, the functions can be expressed as u( x , y ) =

1 1 1 1 ( a − x )( b − y )u1 + x( b − y )u 2 + xyu3 + ( a − x ) yu 4 ab ab ab ab

v( x , y ) =

1 1 1 1 ( a − x )( b − y )v1 + x( b − y )v 2 + xyv3 + ( a − x ) yv4 ab ab ab ab

(13)

From the principle of possible displacements we have ( k ij )r = ∫ σ j ε i dΩ r

(14)

Ω

where σ j are the stresses on the finite element area from displacement d j = 1 , ε i are the strains on the finite element area from displacement d i = 1 . If degrees of freedom have the

25

physical meaning of displacements, ( k ij )r , as an element of the stiffness matrix, is an effort arising along i-th degrees of freedom from j-th unit displacement under condition that all others ( i ≠ j ) degrees of freedom d i = 0 . Thus, we can obtain now the stiffness coefficients using expression (14) and taking into account the plane stress state and expression (5) ab

k ij = h ∫ ∫ (Eε i )T ε j dxdy

(15)

00

where h is the thickness of plate, ε i ( i = 1,2,...,8 ) is the strain vector (4) on the finite element area in the case, when the node displacement with number i is equal to unit but all other displacements are zero, ε j ( j = 1,2 ,...,8 ) is the strain vector (4) on the finite element area in the case, when the node displacement with number j is equal to unit but all other displacements are zero. As an example, let’s express the stiffness matrix element k13 = k u1u 2 presenting the reaction arising in the node 1 along axis x from the unit displacement of the node 2 in the same direction. The numbering of degrees of freedom is given in relation of their recording in the column (8). At the beginning we build the strain vector ε1 = ε u1 corresponding to the deformation state on the finite element area from the unit displacement u1 , when all other nodal displacements are zero. In this case the vector of approximate functions is formed from the expression (13) taking into account u1 = 1 and u 2 = u3 = u 4 = v1 = v 2 = v3 = v 4 = 0 : u( x , y ) =

1 ( a − x )( b − y ) = N1 ab

(16)

v( x , y ) = 0 The strain vector ε1 is formed in relation with the expression (4): ∂ εx ∂x ε1 = ε y = 0 γ ∂ xy ∂y

1 0 − ab ( b − y ) ∂ u( x , y ) 0 = ∂y v( x , y )1 1 − ( a − x ) ∂ ab ∂x

(17)

By the same way, the strain vector ε 3 = ε u 2 is build. This vector corresponds to the deformation state on the finite element area from the unit displacement u 2 , when all other nodal displacements are zero. 1 ab ( b − y ) ε3 = 0 − 1 x ab

(18)

26

Substituting the vectors (17) and (18) into expression (15), we have ab

E 1 b− y E 1 x k13 = h ∫ ∫ − (b − y ) ( a − x ) dxdy = + 2 ab ab 2( 1 + υ ) ab ab 0 0 1−υ Eh ab 3 1 − υ a 3b , =− − 2 6a 2 b 2 1 − υ 2 3a 2 b 2 or after introduction of m = b a : k13 = −

Eh m 1 − υ − 1 − υ 2 3 12m

To obtain remaining elements of the stiffness matrix of quadrilateral finite element under plane stress, the same procedure should be applied and we have k11 =

Eh 1 + υ Eh m 1 − υ Eh 1 3υ , k14 = + , k12 = − + , 6m 1−υ 2 8 8 1−υ 2 8 1−υ 2 3

k15 =

Eh m 1 − υ Eh 1 + υ Eh m 1 − υ , k16 = − − − , k17 = − , 8 8 1−υ 2 1 − υ 2 6 12m 1−υ 2 6

k18 =

Eh 1 3υ Eh 1 1 − υ Eh 1 3υ + m , k 23 = − , k 22 = − , 2 8 2 8 6 1−υ 1−υ 2 8 8 1 − υ 3m

k 24 =

Eh 1 1 − υ Eh 1 + υ Eh 1 1 − υ − m , k 25 = − m , − , k 26 = − 2 6m 2 6 8 1−υ 1−υ 1 − υ 2 6m 12

k 27 =

Eh 1 3υ Eh 1 1 − υ Eh m 1 − υ − − , k 28 = − + m , k 33 = + , 6m 1−υ 2 8 8 1 − υ 2 3m 12 1−υ 2 3

k 34 =

Eh 1 + υ Eh − , k 35 = 8 1−υ 2 1−υ 2

k 37 =

Eh 1 + υ Eh m 1 − υ Eh 1 1 − υ , k 44 = + m , − − , k 38 = 2 6 12m 2 8 6 1 − υ 2 3m 1−υ 1−υ

k 45 =

Eh 1 + υ Eh 1 3υ Eh 1 1 − υ , + m , k 47 = − , k 46 = − 2 8 2 8 1−υ 1 − υ 3m 12 1−υ 2 8

k 48 =

Eh 1 + υ Eh 1 1 − υ Eh m 1 − υ , − − m , k 55 = + , k 56 = 6m 1 − υ 2 6m 12 1−υ 2 8 1−υ 2 3

k 57 =

Eh m 1 − υ Eh 1 3υ Eh 1 1 − υ − + , k 66 = − + , k 58 = + m , 6 1−υ 2 8 8 1 − υ 2 3m 1 − υ 2 3 12m

k 67 =

Eh 1 3υ Eh 1 1 − υ Eh m 1 − υ − m , k 77 = − , k 68 = + , 2 8 2 8 6 6m 1−υ 1 − υ 6m 1−υ 2 3

Eh 1 3υ m 1−υ − , k 36 = − + , 6 6m 1−υ 2 8 8

27

k 78 =

Eh 1 + υ Eh − , k88 = 8 1−υ 2 1−υ 2

1 1−υ + m 6 3m

Based on expressions (5), stresses on the finite element area are determined using known nodes displacements of the system. The strain vector with three components is expressed by the vector of approximate functions by the following way: σ = Eε = EDd

(19)

where E is the matrix of elasticity, D is the matrix of differentiation, d is the vector of approximate functions consisting of two components u( x , y ) and v( x , y ) . The approximate functions (13) are written in the matrix form: 0 ( a − x )( b − y ) 0 ( a − x )( b − y ) x( b − y ) 0 0 x( b − y ) u( x , y ) 1 = xy 0 v( x , y ) ab 0 xy ( a − x )y 0 0 ( a − x ) y

T

u1 v 1 u 2 v 2 u 3 v3 u 4 v 4

(20)

After substituting expressions (3), (4) and (20) into expression (19), we have 1−υ υ( y − b ) (x−a) ( y −b) 2 1−υ υ( x − a ) (x−a) ( y −b) 2 1−υ − ( y − b ) − υ( y − b ) − x 2 1−υ σ x − − − − υ x x ( y b ) E 2 = σ y 2 1−υ τ ( 1 − υ )ab y x υy xy 2 1−υ υx x y 2 1−υ −y − υy − ( x − a ) 2 1−υ − y − υ( x − a ) − ( x − a ) 2

T

u1 v 1 u 2 v 2 u 3 v3 u 4 v 4

(21)

where x and y are coordinates of the points on the finite element area. As it is seen from the expression (21), stresses on the finite element area are the linear functions of coordinates. In the centre of gravity of the finite element (x = a 2 , y = b 2) , the stress vector has the following form:

28

b −2 a − υ 2 b 2 a σ x − υ E 2 σ y = 2 b τ ( 1 − υ )ab 2 xy a υ 2 −b 2 a υ 2

b −υ 2 a − 2 b υ 2 a − 2 b υ 2 a 2 b −υ 2 a 2

1−υ − a 4 1−υ − b 4 1−υ − a 4 1−υ b 4 1−υ a 4 1−υ b 4 1−υ a 4 1−υ − b 4

T

29

u1 v 1 u 2 v 2 u 3 v3 u 4 v 4

PART II

SOLUTION OF FINITE ELEMENT EQUILIBRIUM EQUATIONS

Chapter 4

Solution of equilibrium equations in static analysis

4.1

Introduction

When FEM is used for a solution of static problems, we get a set of simultaneous linear equations, which can be stated in the form KX = F where K is the stiffness matrix of a structure, X is the displacement vector and F is the load vector. This equation can be expressed in the scalar form as k x + k x + K k x = f 1n n 1 11 1 12 2 k 21 x1 + k 22 x 2 + K k 2n x n = f 2 M k x + k x + K + k x = f n2 2 nn n n n1 1 where the coefficients k ij and the constants f i are given. The problem is to find the values of xi , if they exist. In the matrix form k11 k K = 21 n×n M k n1

K k1n K k 2n , k n 2 K k nn k12 k 22

x1 x X = 2, n×1 M x n

f1 f F = 2 n×1 M f n

It is necessary to note that in the finite element analysis, the order of the matrix K is very large. The methods available for solving of the systems of linear equations can be divided into two types: direct or iterative. Direct methods are those, which, in the absence of round-off and other errors, will yield the exact solution in a finite number of elementary arithmetic operations. In practice, because a computer works with a finite word length, sometimes the direct methods do not give good solutions. Indeed the errors arising from round-off and truncation may lead to extremely poor or even useless results. The fundamental method used for direct solutions is Gaussian elimination, but even within this class there are a variety of choices of methods, which vary in computational efficiency and accuracy. Iterative methods are those, which start with an initial approximation, and which by applying a suitably chosen algorithm, lead to successively better approximations. When the process converges, we can expect to get a good approximate solution. The accuracy and

30

the rate of convergence of iterative methods vary with the algorithm chosen. The main advantages of iterative methods are the simplicity and uniformity of the operations to be performed, which make them well suited for use on computers and their relative in sensitivity to the growth of round-off errors. Matrices associated with linear systems in the finite element analysis are classified as sparse and have very few nonzero elements. Fortunately, in most finite element applications, the matrices involved are positive definited, symmetric and banded. Hence solution techniques, which take advantage of the special character of such systems of equations, have also been developed. 4.2

Gaussian elimination method

The basic objective of this method is to transform the given system into an equivalent triangular system, whose solution can be more easily obtained. We shall consider the following system of three equations to illustrate the process x1 2 x1 4 x 1

−

+ 3 x3 + x3 − x3

x2

+ 3x2 + 2 x2

= 10 = 15 = 6

To eliminate the terms x1 from equations (2) and (3), we multiply equation (1) by -2 and 4 and add respectively to equation (2) and equation (3) living the first equation unchanged. We will have then x1

−

x2 5x2 6 x2

+

3 x3

− 5 x3 − 13 x3

=

10

= −5 = − 34

To eliminate the term x2 from equation (3), we multiply equation (2) by −

6 and add to 5

equation (3). We will have the triangular system x1

−

x2 5 x2

+ 3 x3 − 5 x3 − 7 x3

=

10

= −5 = − 28

This triangular system can be solved now by the back substitution. From equation (3) we find x3 = 4 . Substituting this value for x3 into equation (2) and solving for x2 , we obtain x 2 = 3 . Finally, knowing x3 and x2 , we can solve equation (1) for x1 , obtaining x1 = 1 . 4.3

Generalisation of Gauss method Let’s the given system of equations be written as

31

k11 x1 k 21 x 2 k 31 x3 0 0

+ k12 x 2 + k 22 x 2 + k 32 x 2 + k 42 x 2 + 0

+ + + + +

k13 x3 k 23 x3 k 33 x3 k 43 x3 k 53 x3

+ + + + +

0 k 24 x 4 k 34 x 4 k 44 x 4 k 54 x 4

+ 0 + 0 + k 35 x5 + k 45 x5 + k 55 x5

= = = = =

f1 f2 f3 f4 f5

The wide of matrix band is 3. After elimination of x1 , we have k11 0 0 0 0

k12

k13

0

(1) k 22

(1 ) k 23

(1 ) k 24

(1) k 42

(1 ) k 43

(1 ) k 44

(1) k 32

0

(1 ) k 33 (1 ) k 53

(1 ) k 34 (1 ) k 54

0 0 (1 ) k 35 ; (1 ) k 45 (1 ) k 55

f (1) 1 f 2( 1 ) (1) f3 (1) f4 f (1) 5

where new coefficients are expressed by the following way k (1) k 22 = k 22 − k 21 12 k11 k (1 ) = k 23 − k 21 13 k 23 k11 k (1) k 32 = k 32 − k 31 12 k11 k (1) k 33 = k 33 − k 31 13 k11 f f 2( 1 ) = f 2 − k 21 1 k11 f f 3( 1 ) = f 3 − k 31 1 k11 The upper index (1) has been used to denote the first elimination. The general relation of an arbitrary coefficient after first elimination has the following form k ij( 1 ) = k ij − k i1

k1 j k11

, i, j > 1

To elimination with number n corresponds the following general relation k ( n −1 ) (n) ( n −1 ) ( n −1 ) nj − k in , k ij = k ij ( n −1 ) k nn

i, j > n

Analogous formulas are obtained for a vector of load

32

( n −1 ) (n) ( n −1 ) ( n −1 ) f n fi , = fi − k in ( n −1 ) k nn

i>n

The symmetry in coefficients after elimination operation maintains. Thus the matrix decomposition may be carried out using only coefficients situated on the main diagonal and above it. Therefore, it is not necessary to store the full matrix. Another positive feature of the Gaussian elimination, which can be used, is that matrix elements situated out of band do not influence on the elimination process. They equal to zero. Hence, it is not necessary to store them. That gives the possibility to store the global stiffness matrix as a rectangular array with wide equal to the wide of the matrix band. After applying the above procedure ( n − 1 ) times, the original system of equations reduces to the following single equation ( n −1 )

k nn

( n −1 )

xn = f n

from which we can obtain xn =

−1 f n( n ) ( n −1 )

k nn

The values of the remaining unknowns can be found by the back substitution. Note:

If zero or negative diagonal element occurs in the Gauss elimination the structure is not stable or, by other words, the global stiffness matrix is not positive definite.

The Gauss elimination scheme falls under the category of direct methods. This category includes also Choleski method (a direct method for solving a linear system which makes use of the fact that any square matrix can be expressed as the product of an upper and lower triangular matrices), the Givens factorization (rotation matrices are used to reduce the global stiffness matrix into upper triangular form), the Householder factorization (reflection matrices are used). 4.4

Simple vector iterations

The power and inverse iteration methods are the methods not used widely now, but they should be examined, since help to understand more complex modern algorithms. Algorithm of the power method: The unit vector X 0 is chosen. Then, for k = 1,2,3,... 1) To form Fk = KX k −1 F 2) To set the rate of X k = k Fk 3)

To subject X k on the convergence test.

33

The inverse iteration method is a power method applied to K −1 . It is not necessary to make an inverse matrix K , instead of that we change Fk = KX k −1 in the power method on 1’)

To solve in relation of Fk the following equation KFk = X k −1 .

In the class of iterative, the Gauss-Seidel method is well known. The conjugate gradient and Newton’s methods are other iterative methods based on the principle of unconstrained minimisation of a function. It is to be noted that the indirect methods are less popular than the direct methods in solving large systems of linear equations. 4.5

Introduction to nonlinear analyses

In the linear analysis we assumed that displacements of the finite element assemblage are infinitesimally small and that a material is linearly elastic. In addition we also assumed that a nature of boundary conditions remains unchanged during application of loads on the finite element assemblage. Figure 4.1 gives a classification that is used very conveniently in practical nonlinear analysis because this classification considers separately material nonlinear effects and kinematic nonlinear effects. The basic problem in a general nonlinear analysis is to find the state of equilibrium of a body corresponding to the applied load. Assuming that the externally applied loads are described as a function of time, the equilibrium conditions of a system of finite elements representing the body under consideration can be expressed as t

F−t R = 0

(1)

where the vector t F stores the externally applied nodal loads and t R is the vector of nodal point forces that are equivalent to the element stresses. This relation must express the equilibrium of the system in the current deformed geometry taking due account of all nonlinearities. Considering the solution of the nonlinear response, we recognise that the equilibrium relation examined must be satisfied throughout the complete history of load application, i.e. the time variable may take on any value from zero to the maximum time of interest. The basic approach in the incremental step-by-step solution is to assume that the solution for the discrete time t is known, and that the solution for the discrete time t + ∆t is required, where ∆t is a suitably chosen time increment. Hence, considering the previous relation at time t + ∆t we have t + ∆t

F − t + ∆t R = 0

(2)

Since the solution is known at time t , we can write t + ∆t

R=t R + R

(3)

where R is the increment in nodal point forces corresponding to the increment in element displacements and stresses from time t to time t + ∆t . This vector can be approximated using a tangent stiffness matrix t K which corresponds to the geometric and material conditions at time t

34

σ

∆

P/2

σ ,ε

L

E

P/2 L

1

σ = P/ A ε =σ / E

ε

ε < 0.04

∆ =ε L

(a) Linear elastic (infinitesimal displacements).

∆

σY

P/2 L

1

σ = P/ A σ y σ −σY + ε= P/2 E ET

L

P/ A

σ

ET

E 1

ε

ε < 0.04

(b) Materially-nonlinear-only (infinitesimal displacements, but nonlinear stress-strain relation).

∆'

y'

x' y

ε'

L

L

ε ' < 0.04 ∆' = ε ' L

x

(c) Large displacements and large rotations but small strains. Linear or nonlinear material behaviour. Fig. 4.1 Classification of analyses.

35

(d) Large displacements, large rotations and large strains. Linear or nonlinear material behaviour.

P/2

P/2

∆

(e) Change in boundary condition at displacement ∆. Fig. 4.1 (continuation) R ≅ t KX

(4)

where X is a vector of incremental nodal point displacements. Substituting two last expressions (3) and (4) into previous, we obtain t

KX = t + ∆t F − t R

(5)

and solving for X , we can calculate an approximation to the displacements at time t + ∆t t + ∆t

X= t X + X

(6)

The exact displacements at time t + ∆t are those that correspond to the applied loads t + ∆t

F . We calculate in equation (6) only an approximation to these displacements because equation (4) was used. Having evaluated an approximation to the displacements corresponding to time t + ∆t , we can now solve for an approximation to the stresses and corresponding nodal point forces at time t + ∆t , and could then proceed to the next time increment calculations. However, because of the assumption in equation (4) such a solution may be subject to very significant errors and, depending on the time or load step sizes used, may indeed be

36

unstable. In practice, it is therefore frequently necessary to iterate until the solution of equation (2) is obtained to sufficient accuracy. A widely used iteration procedure is the modified Newton iteration, in which we solve for i = 1,2,3,... ∆F (i −1 ) = t + ∆t F − t + ∆t R (i −1 ) t

(7)

∆K (i) = ∆F (i −1 )

t + ∆t

(8)

X(i) = t + ∆t X(i −1 ) + ∆X(i)

(9)

with the initial conditions t + ∆t

X( 0 ) = t X ;

t + ∆t

R( 0 ) =t R

These equations were obtained by linearizing the response of the finite element system about the conditions at time t . In each iteration we calculate in equation (7) an out-ofbalance load vector, which yields an increment in displacements obtained in equation (8), and we continue the iteration until the out-of-balance load vector ∆F (i −1 ) or the displacement increments ∆X(i) are sufficiently small. The most frequently used iteration schemes for the solution of nonlinear finite element equations are some forms of Newton-Raphson iteration, when the equation (7), (8), (9) are a special case. In the Newton-Raphson iteration, in general the major computational cost per iteration lies in the calculation and factorisation of the tangent stiffness matrix. As an alternative to forms of Newton iteration, a class of methods known as matrix update methods or quasi-Newton methods has been developed for iteration on nonlinear systems of equations. These methods involve updating the stiffness matrix to provide a secant approximation to the matrix from iteration ( i − 1 ) to i .

4.6

Convergence criteria

If a solution strategy based on iterative methods is to effective, realistic criteria should be used for the termination of the iteration. At the end of each iteration, the solution obtained should be checked to see a convergence. If the convergence tolerances are too loose, inaccurate results are obtained, and if the tolerances are too tight, much computational effort is spent to obtain needless accuracy. Since we are seeking the displacement corresponding to time t + ∆t , it is natural to require that the displacements at the end of each iteration be within a certain tolerance of the true displacement solution. Hence, a realistic convergence criterion is X(i) t +t

X

2 ≤e D 2 1/ 2

n 2 where e D is a displacement convergence tolerance, X 2 = ∑ x k k =1 norm. The vector t + ∆t X is not known and must be approximated.

37

is the Eucledean

A second convergence criterion is obtained by measuring the out-of-balance load vector. For example, we may require that the norm of the out-of-balance load vector be within a preset tolerance e F of the original load increment t +t

F − t + t R (i)

t +t

t

F− R

2 ≤e F

2

In the case, when both the displacements and the forces are near their equilibrium values, a third convergence criterion may be useful, in which the iteration (i.e. the amount of work done by the out-of-balance loads on the displacement increments) is compared to the initial internal energy increment. Convergence is assumed to be reached when, with e E a preset energy tolerance, X(i)T ( t + t F − t + t R (i −1 ) ) X( 1 )T ( t + t F − t R )

≤ eE

38

Chapter 5

5.1

Solution of eigenproblems

Introduction

The forced vibration equation after the finite element discretization of structure can be expressed as follows && + KX = F MX

(1)

where M and K are the mass and stiffness matrices of structure; F is the external load && are the displacement and acceleration vectors. In the free vibration vector; X and X analysis, the external load vector is zero and the displacements are harmonic as X = Xe iωt

(2)

After substitution of Equation (2) into Equation (1), we obtain

[K − ω M]X = 0 2

or KX = λMX

(3)

where X represents the amplitudes of the displacements X called the mode shape or eigenvector, ω denotes the natural frequency of vibration and ω 2 = λ is called eigenvalue. Equation (3) is the generalised eigenvalue problem. It has a nonzero solution

[

]

for X , when the determinant of the following matrix K − ω 2 M is zero, i.e. K − ω 2M = 0

(4)

The structure with n degrees of freedom has n natural frequencies. It was assumed that the rigid body degrees of freedom were eliminated in Equation (1). If rigid body degrees of freedom are not eliminated in deriving the matrices K and M , some of the natural frequencies ω would be zero. In such case, for a general three-dimensional structure, there will be six rigid body degrees of freedom and hence six zero frequencies. In most of the numerical methods used for the solution of Equation (3), the generalised eigenvalue problem is first converted into the form of a standard eigenvalue problem, which can be stated as HX = λX or

[H − λI ]X = 0

(5)

Multiplying Equation (3) by M −1 , we obtain Equation (5), where H = M −1K

(6)

However, in this form the matrix H is in general nonsymmetric, although M and K are both symmetric. Since a symmetric matrix is desirable from the points of view of storage

39

and computer time, we adopt the following procedure to derive a standard eigenvalue problem with symmetric H matrix. Assuming that M is symmetric and positive definite, we use Choleski decomposition and express M as M = U T U , where U is the upper triangular matrix. By substitution of M into Equation (3), we obtain KX = λU T UX i.e. U T −1KX = λUX i.e. U T −1KU −1UX = λUX

(7)

Defining a new vector Y as Y = UX , Equation (7) can be written as a standard eigenvalue problem as

[H − λI ]Y = 0 where the matrix H is now symmetric and is given by H = U T −1KU −1 Then we apply the inverse transformation to obtain the desired eigenvectors X i = U −1Yi corresponding to the eigenvalues λi . Two general types of methods, namely, transformation methods and iterative methods, are available for solving eigenvalue problems. The transformation methods such as Jacobi, Givens and Householder schemes are preferable, when all the eigenvalues and eigenvectors are required and the dimension of the eigenvalue problem is small. The iterative methods such as the power method, subspace iteration and Lanczos methods are preferable, when few eigenvalues and eigenvectors are required only and the eigenvalue problem has a large dimension. 5.2

Transformation methods Transformation methods employ the basic properties of eigenvectors in the matrix

X, XT KX = Λ

(8)

XT MX = I

(9)

Since the matrix X , of order n × n , which diagonalizes K and M in the way given in Equations (8) and (9) is unique, we can try to construct it by iteration. The basic scheme is to reduce K and M into diagonal form using successive pre- and post-multiplication by matrices PkT and Pk , respectively, where k = 1,2,... Specifically, if we define K 1 = K and M1 = M , we form

40

K 2 = P1T K 1P1 K 3 = P2T K 2 P2 M T K k +1 = Pk K k Pk M

M 2 = P1T M1P1 M 3 = P2T M 2 P2 M T M k +1 = Pk M k Pk M

(10)

(11)

where the matrices Pk are selected to bring K k and M k closer to diagonal form. Then K k +1 → Λ and M k +1 → I as k → ∞ and l is examined as the last iteration, X = P1P2 ...Pl . In practice, it is not necessary that M k +1 converges to I and K k +1 to diagonal form. Namely, if K k +1 → diag ( K r ) and M k +1 → diag ( M r ) as k → ∞ then with l indicating the last iteration K ( l +1 ) Λ = diag r M ( l +1 ) r 1 X = P1P2 ...Pl diag M ( l +1 ) r

The basic idea described above is used in Jacobi and Householder-QR methods applied effectively in the finite element analysis. 5.3

Jacobi method

The basic Jacobi method has been developed over a century ago for the solution of standard eigenproblems. A major advantage of the procedure is its simplicity and stability. Considering the standard eigenproblem KX = λX , the k iteration step defined in Equation (10) reduces to K k +1 = PkT K k Pk

(12)

where Pk is an orthogonal matrix, i.e. Equation (11) gives PkT Pk = I In Jacobi solution the matrix Pk is a rotation matrix, which is selected in such way that an off-diagonal element in K k becomes zero. If element (i, j ) is to be reduced to zero, the corresponding orthogonal matrix Pk is

41

i − column

j − column

M M 1 O M M 1 M M − sin θ L L L i − row cos θ O Pk = 1 L L L j − row sin θ cos θ 1 O 1 where θ is selected from the condition that element (i, j ) in K k +1 be zero. Denoting element (i, j ) in K k by kij( k ) , we use tan 2θ =

θ =

2k ij( k ) k ii( k ) − k (jjk )

π 4

for k ii( k ) ≠ k (jjk ) (k )

for k ii

(k )

= k jj

Since K k is symmetric for all k , the upper (or lower) triangular part of the matrix, including its diagonal elements, is used. It is necessary to note that although the transformation in Equation (12) reduces an off-diagonal element in K k to zero, this element will again become nonzero during the transformations that follow. Therefore for the design of an actual algorithm, we have to decide which element to reduce to zero. One choice is to always zero the largest offdiagonal element in K k . However, the search for the largest element is time consuming and it may be preferable to simply carry out the Jacobi transformations systematically, row-by-row or column-by-column, which is known as the cyclic Jacobi procedure. The disadvantage of this procedure is that the element may already be nearly zero and a rotation is still applied. A procedure that has been used very effectively is the threshold Jacobi method, in which the off-diagonal elements are tested consequently, namely rowby-row (or column-by-column), and a rotation is only applied if the element is larger than the threshold. 5.4

Vector iteration methods In the vector iteration methods the basic relation is K X = λM X

(13)

42

If assume a vector for X , say X1 , and assume a value for λ , say λ = 1 , we can evaluate then the right-hand side of Equation (13), i.e. we may calculate R1 = ( 1 )MX1 Since X1 is an arbitrarily assumed vector, we do not have, in general, that KX1 = R1 . Instead, we have the following equation KX 2 = R1 , X 2 ≠ X1 We may assume that X 2 may be a better approximation to an eigenvector than was X1 . By repeating the cycle we obtain an increasingly better approximation to an eigenvector. The procedure described above is the basic of inverse iteration. We will see that other vector iteration techniques work in a similar way. Specifically, in forward iteration, in the first step we evaluate R1 = KX1 and then obtain the improved approximation X 2 to the eigenvector solving MX 2 = R1 . 5.5

Subspace iteration method

This method is very effective in finding the first few eigenvalues and the corresponding eigenvectors of large eigenvalue problems. The various steps of this method are given below briefly. Algorithm: 1)

Start with q initial iteration vectors X1 , X 2 ,..., X q , q > p , where p is the number of eigenvalues to be calculated. Bathe and Wilson suggested a value of q = min( 2 p , p + 8 ) for good convergence. Define the initial modal matrix X 0 as X 0 = X1X 2 ...X q

[

]

and set the iteration number as k = 0 . 2)

Use the following subspace iteration procedure to generate an improve modal matrix X k +1 : ~ a) Find X k +1 from the relation ~ KX k +1 = MX k b) Compute ~ ~ K k +1 = XTk +1KX k +1 ~ ~ M k +1 = XTk +1MX k +1 c)

d)

Solve for the eigenvalues and eigenvectors of the reduced system K k +1Q k +1 = M k +1Q k +1Λ k +1 and obtain Λ k +1 and Q k +1 Find an improved approximation to the eigenvectors of the original system as

43

~ X k + 1 = X k + 1 Q k +1 Note: (1) It is assumed that the iteration vectors converging to the exact eigenvectors X1exact , X exact ,... are stored as the 2 (2)

(3)

columns of the matrix X k +1 . It is assumed that the vectors in X 0 are not orthogonal to one of the required eigenvectors.

If λ(i k ) and λ(i k +1 ) denote the approximations to the i eigenvalue in the iterations k − 1 and k respectively, we assume convergence of the process whenever the following criteria is satisfied: ( k +1 )

λi

(k )

− λi

( k +1 ) λi

≤ ε , i = 1,2,..., p

where ε ≅ 10 −6 .

44

Chapter 6

6.1

Solution of equilibrium equations in dynamic analysis

Introduction The dynamic equation of motion of a structure can be written as && + CX & + KX = F MX

(1)

where M , C and K are the mass, damping and stiffness matrices of structure, F is the & and X && are the displacement, velocity and acceleration vectors external load vector, X , X of a finite element assemblage. It should be noted that Equation (1) is derived from considerations of static at time t, i.e. Equation (1) may be written as FI (t) + FD (t) + FE (t) = F(t)

(2)

&& ; F (t) are the damping forces, F (t) = CX &; where FI (t) are the inertia forces, FI (t) = MX D D FE (t) are the elastic forces, FE (t) = KX , all of them are time-dependent. Mathematically, Equation (1) represents a system of linear differential equations of second order and, in principle, the solution of the equations can be obtained by standard procedures for the solution of differential equations with constant coefficients. However, these procedures can be very expensive if the order of the matrices is large. Therefore few effective methods have been elaborated to apply them in practical finite element analysis. These methods are divided into direct integration and mode superposition. Although these two techniques may at first sight appear to be quite different, in fact, they are closely related, and the choice for one method or the other is determined only by their numerical effectiveness. We will consider only solution of the linear equilibrium equations (Equation (1)). 6.2

Direct integration methods

In direct integration the equations in (1) are integrated using a numerical step-by-step procedure, the term “direct” meaning that prior to the numerical integration, no transformation of the equations into a different form is carried out. In essence, direct numerical integration is based on two ideas. First, instead of trying to satisfy Equations (1) at any time t , it is assumed to satisfy Equations (1) only at discrete time intervals ∆t apart. Therefore it appears that all solution techniques employed in static analysis can probably also be used effectively in direct integration. The second idea on which a direct integration method is based is that a variation of displacements, velocities and accelerations within each time interval ∆t is assumed. This assumption determines the accuracy, stability and cost of the solution procedure.

45

In the solution the time span under consideration T is subdivided into n equal time T intervals ∆t (i.e. ∆t = ) and the integration scheme employed establishes an n approximate solution at times 0 , ∆t ,2∆t ,3∆t ,...,t ,t + ∆t ,...,T . Since an algorithm calculates the solution at the previous times considered, we derive the algorithms by assuming that the solutions at times 0,∆t ,2∆t ,3∆t ,...,t are known and that the solution at time t + ∆t is required next. Let’s have a look now a more effective and widely used in the general purpose finite element programs (including ANSYS), the Newmark method. 6.3

The Newmark method

The Newmark integration scheme can be examined as an extension of the linear acceleration method. The following assumptions are used

[

]

t + ∆t & t & && + δ t + ∆t X && ∆t X= X + (1 − δ ) t X t + ∆t

(3)

& ∆t + 1 − α t X && + α t + ∆t X && ∆t 2 X= t X+ t X 2

(4)

where α and δ are parameters that can be determined to obtain integration accuracy and stability. The linear acceleration method: a linear variation of acceleration from time t to time t + ∆t is assumed. When δ = 1 2 and α = 1 6 , relations (3) and (4) correspond to the linear acceleration method. Newmark originally proposed, as an unconditionally stable scheme, the constantaverage-acceleration method (also called trapezoidal rule), in which case δ = 1 2 and α = 1 4. In addition to (3) and (4), for solution of the displacements, velocities and accelerations at time t + ∆t , the equilibrium Equations (1) at time t + ∆t are also considered && + C t + ∆t X & + K t + ∆t X= t + ∆t F M t + ∆t X

(5)

(

1 t && t + ∆t && X+ X 2

t && X

)

t + ∆t &&

X

t + ∆t

t

Fig. 6.1 Newmark’s constant-average-acceleration scheme.

46

&& in terms of t + ∆t X and then substituting for t + ∆t X && into (3), Solving from (4) for t + ∆t X && and t + ∆t X & , each in terms of the unknown displacements we obtain equations for t + ∆t X t + ∆t

& and t + ∆t X && are substituted into (5) to solve X only. These two relations for t + ∆t X && and t + ∆t X & can also be calculated. for t + ∆t X , after which, using (3) and (4), t + ∆t X Step-by-step solution using Newmark integration method: A. Initial calculations: 1. Form stiffness matrix K , mass matrix M and damping matrix C . & and 0 X && . 2. Initialise 0 X, 0 X 3.

4.

Select time step size ∆t , parameters α and δ , and calculate integration constants: δ ≥ 0.50 ; α ≥ 0.25( 0.5 + δ ) 2 1 1 1 δ ; a2 = ; a3 = −1 ; a0 = ; a1 = α∆t α∆t 2α α∆t 2 ∆t δ δ a 4 = − 1 ; a5 = − 2 ; a6 = ∆t (1 − δ ) ; a7 = δ∆t 2 2 α Form effective stiffness matrix K : K = K + a0 M + a1C .

5. Triangularize K : K = LDLT B. For each time step: 1. Calculate effective loads at time t + ∆t : t + ∆t & +a tX && ) + C(a t X + a t X & +a tX && ) F = t + ∆t F + M(a0 t X + a 2 t X 3 1 4 5 2. Solve for displacements at time t + ∆t : LDK T t + ∆t X= t + ∆t F 3.Calculate accelerations and velocities at time t + ∆t : t + ∆t && & − a tX && X = a t + ∆t X− t X − a t X 0

(

)

t + ∆t & t & && + a X= X + a6 t X 7

6.4

2 t + ∆t &&

3

X

Mode superposition

The numbers of operations required in the direct integration are directly proportional to the number of time steps used in the analysis. Therefore, in general, the use of direct integration can be expected to be effective, when the response for a relatively short duration is required. However, if the integration must be carried out for many time steps, it may be more effective to first transform the equilibrium equations (Equation (1)) into a form in which the step-by-step solution is less costly. In particular, since the number of operations required is directly proportional to the half-bandwidth mk of the stiffness matrix, a reduction in mk would decrease proportionally the cost of the step-by-step solution.

47

6.5

Change of basis to modal generalised displacements

We propose to transform the equilibrium equations into a more effective form for direct integration by using the following transformation on the finite element nodal point displacements X X(t) = PU(t)

(6)

where P is a square matrix and U(t) is a time-dependent vector of order n . The transformation matrix P is still unknown and will have to be determined. The components of U are refered to as generalised displacements. Substituting (6) into (1) and premultiplying by P T , we obtain && (t) + C U & (t) + KU(t) = F (t) MU

(7)

M = P T MP, C = P T CP, K = P T KP, F = P T F

(8)

The objective of the transformation is to obtain new system stiffness, mass and damping matrices, K, M and C , which have a smaller bandwidth than the original system matrices, and the transformation matrix P should be selected accordingly. In theory, there can be many different transformation matrices P , which would reduce the bandwidth of the system matrices. However, in practice, an effective transformation matrix is established using the displacement solutions of the free vibration equilibrium equations with damping neglected, && + KX = 0 MX

(9)

By substitution X = Φ e iωt , Equation (9) becomes KΦ = ω 2MΦ

(10)

where Φ is an eigenvector and ω 2 is an eigenvalue, and Equation (10) is called the generalised eigenvalue problem with unknowns Φ and ω 2 . The eigenproblem (10) yealds

(

)(

)

(

)

the n eigensolutions ω12 , Φ1 , ω 22 , Φ 2 , …, ω n2 ,Φ n , where the eigenvectors are Morthonormalized, i.e. = 1, i = j ΦiT M Φi = 0 , i ≠ j

(11)

0 ≤ ω12 ≤ ω 22 ≤ ... ≤ ω n2

(12)

The vector Φi is called the i –mode shape vector and ω i is the corresponding frequency of vibration. Defining a matrix Φ , whose columns are the eigenvectors Φi and a diagonal matrix Ω 2 , which stores the eigenvalues ω i2 on its diagonal, i.e.

48

Φ = [Φ1 , Φ 2 ,..., Φ n ] ,

ω12 ω 22 2 Ω = O ω n2

(13)

we can write the n solutions to (10) as KΦ = MΦΩ 2

(14)

Since the eigenvectors are M–orthogonal, we have ΦT KΦ = Ω 2 ,

ΦT MΦ = I

(15)

It is now apparent that the matrix Φ would be a suitable transformation matrix P in (6). Using X(t) = ΦU(t)

(16)

we obtain equilibrium equations that correspond to the modal generalised displacements && (t) + ΦT CΦU & (t) + Ω 2 U(t) = Φ T F(t) U

(17)

The initial conditions on U(t) are obtained using (16) and the M–orthonormality of Φ , i.e. at time 0 we have 0

U = ΦT M 0 X ,

0 & & U = ΦT M 0 X

(18)

The equations in (17) show that if a damping matrix is not included in the analysis, the finite element equilibrium equations are decoupled, when using in the transformation matrix P the free vibration mode shapes of the finite element system. Since the derivation of the damping matrix can in many cases not be carried out explicitly, but the damping effects can only be included approximately, it is reasonable to use a damping matrix that includes all required effects, but at the same time allows an effective solution of the equilibrium equations. 6.6

Analysis with damping neglected

If velocity-dependent damping effects are not included in the analysis, each equation presents the equilibrium equation of a single degree of freedom system with unit mass and stiffness ω i2 . In summary, the response analysis by mode superposition requires first, the solution of the eigenvalues and eigenvectors of the problem, then the solution of the decoupled equilibrium equations and, finally, the superposition of the response in each eigenvector. The choice of whether to use direct integration or mode superposition will be decided by considerations of effectiveness only. The essence of a mode superposition solution of a dynamic response is that frequently only a small fraction of the total number of decoupled equations need be considered, in order to obtain a good approximate solution to the actual response of the system. However, the finite element mesh should be chosen such that all important exact frequencies and vibration mode shapes are well approximated. In this case, the mode superposition procedure can be much more effective than direct integration.

49

For earthquake loading, in some cases only the 10 lowest modes need be considered, although the order of the system n may be larger than 1000. On the other hand, for blast or shock loading, many more modes need generally be included, and number of eigenmodes required p may be as large as 2 n /3. Finally, in vibration excitation analysis, only a few intermediate frequencies may be excited, such as all frequencies between the lower and upper frequency limits ω l and ω u , respectively. Considering the problem of selecting the number of modes to be included in the mode superposition analysis, it should always be kept in mind that an approximate solution to the dynamic equilibrium equations is sought. In summary, assuming that the decoupled equations have been solved accurately, the errors in a mode superposition analysis using p < n are due to the fact that not enough modes have been used, whereas the errors in a direct integration analysis arise because too large a time step is employed. 6.7

Analysis with damping included

Considering the analysis of system in which damping effects can not be neglected, we still would like to deal with decoupled equilibrium equations, merely to be able to use essentially the same computational procedure whether damping effects are included or neglected. In general, the damping matrix C can not be constructed from element damping matrices, such as the mass and stiffness matrices of the element assemblage, and its purpose is to approximate the overall energy dissipation during the system response. The mode superposition analysis is particularly effective if it can be assumed that damping is proportional, in which case ΦiT C Φ j = 2ω iξ i δ ij

(19)

where ξ i is a modal damping parameter and δ ij is the Kronecker delta (δ ij = 1 for i = j , δ ij = 0 for i ≠ j ). Therefore, using Equation (19), it is assumed that the eigenvectors Φi , i = 1,2 ,..., n are also C–orthogonal and the equations in (17) reduce to n equations of the form u&&i ( t ) + 2ω i ξ i u& i ( t ) + ω i2 u i ( t ) = f i ( t )

(20)

where f i ( t ) and the initial conditions on u i ( t ) are defined as f i ( t ) = ΦiT F( t ) ,

i = 1,2,..., n

(21)

ui t = 0 = ΦiT M 0 X & u& = ΦT M 0 X i t =0

(22)

i

We note that Equation (20) is the equilibrium equation governing motion of the single degree of freedom system. In considering the implications of using (19) to take account of damping effects, the following observations are made. Firstly, the assumption in (19) means that the total damping in the structure is the sum of individual damping in each mode. The damping in one mode could be observed, for example, by improving initial conditions corresponding

50

to that mode only (i.e. 0 X = Φi for mode i ) and measuring the amplitude decay during the free damped vibration. A second observation is that in the numerical solution we do not calculate the damping matrix C , but only the stiffness and mass matrices K and M . However, assume that it would be numerically more effective to use the direct stepby-step integration and that the realistic damping ratios ξ i , i = 1,2 ,..., p are known. In that case, it is necessary to evaluate the matrix C explicitly, which when substituted into (19) yields the established damping ratios ξ i . If p = 2 , Rayleigh damping can be assumed, which is of the form C = αM + βK

(23)

where α and β are constants to be determined from two given damping ratios that correspond to two unequal frequencies of vibration by the following way ΦiT ( αM + βK )Φi = 2ω iξ i

(24)

α + βω i2 = 2ω iξ i

In actual analysis it may well be that the damping ratios are known for many more than two frequencies. In that case two average values, say ξ1 and ξ 2 , are used to evaluate α and β . If more than only two damping ratios are used to establish C , a more complicated damping matrix may be suggested. Assume that the p damping ratios ξ i , i = 1,2 ,..., p are given to define C . Then a damping matrix that satisfies the relation in (19) is obtained using the Caughey series, C=M

∑ a k [M −1K ]

p −1

k

(25)

k =0

where the coefficients a k , k = 1,2,..., p are calculated from the p simultaneous equations ξi =

1 a0 + a1ω i + a 2ω i3 + ... + a p −1ω i2 p −3 2 ωi

(26)

We should note that with p = 2 , Equation (25) reduces to Rayleigh damping, as presented in (23). An important observation is that if p > 2 , the damping matrix C in (25) is, in general, a full matrix that considerably increases the cost of analysis. Therefore Rayleigh damping is assumed in most practical analyses using direct integration. A disadvantage of Reyleigh damping is that the higher modes are considerably more damped than the lower modes, for which the Rayleigh constants have been selected. In the case of nonproportional damping (analysis of structures with widely varying material properties), it may be resonable to assign in the construction of the damping matrix different Rayleigh coefficient α and β to different parts of the structure, which results into a damping matrix that does not satisfy the relation in (19). Another case of nonproportional damping is encountered, when concentrated dampers corresponding to specific degrees of freedom (e.g. at the support points of a structure) are specified. The solution of the finite element system equilibrium equatoins with nonproportional damping can be obtained using the direct integration algorithms without modifications, because the property of the damping matrix did not enter the derivation of the solution

51

procedures. In the mode superposition method for the case of nonproportional damping, the equilibrium equations in the basis of mode shape vectors are no longer decoupled. An exact mathematical formulation may be presented in an alternative analysis procedure, where the decoupling of the finite element equilibrium equations is achieved by solving a quadratic eigenproblem, in which case complex frequencies and vibration mode shapes are calculated.

52

PART III

Chapter 7

7.1

EMPLOYMENT OF THE FINITE ELEMENT METHOD

Some modelling considerations

Introduction

An establishment of appropriate finite element model for an actual practical problem depends to a large degree on the following factors: understanding of the physical problem including a qualitative knowledge of the structural response to be predicted, knowledge of the basic principles of mechanics and good understanding of the finite element procedures available for analysis. Discretization of the domain into finite elements is the first step in the finite element method. This is equivalent to replacing the domain having an infinite number of degrees of freedom by a system having finite number of degrees of freedom. The shape, size, number and configuration of elements have to be chosen carefully so that the original body or domain is simulated as closely as possible without increasing the computational effort needed for the solution. The various considerations taken in the discretization process are: • type of elements • size of elements • location of nodes • number of elements • simplifications afforded by the physical configuration of the body • finite representation of infinite bodies • node numbering scheme • automatic node generation After meshing of the body it is necessary to add the material properties, external loads, and apply the boundary conditions. Before start of the problem, only parameters of the calculation regime should be added to the input file. 7.2

Type of elements

Often the type of elements to be used will be evident from the physical problem itself and geometry of the body. Let’s consider briefly various types of finite elements, which are subject to certain static and kinematic assumptions. •

Truss and beam elements Truss and beam elements are very widely used in structural engineering (Fig. 7.1, 7.2) to model for example building frames and bridges.

53

FE y Across section A-A: σ x is uniform, all other stress components are zero. A

x

A

P1

P2

Fig. 7.1 Uniaxial stress condition: frame under concentrated loads.

FE

y

Across section A-A:

P

M

A

σ x ≠ 0, σ y ≠ 0, τ xy ≠ 0

x

A

Fig. 7.2 Beam under concentrated load and moment. y q

FE x

structural 2-D solid

triangle

rectangle

y M

q

x

σ x , σ y , τ xy are uniform across the thickness. All other stress components are zero.

P

Fig. 7.3 Plane stress conditions: membrane and beam under in-plane actions.

54

z

y

y z x

u ≠ 0, v ≠ 0, w = 0 ε z = γ xz = γ yz = 0 σz ≠ 0

1

x Fig. 7.4 Plane strain conditions: long dam subjected to water pressure and long retaining wall subjected to soil pressure. •

Plane stress and plane strain elements Plane stress elements are employed to model membranes, the in-plane action of beams and plates and so on (Fig. 7.3). In each of these cases a two-dimensional stress situation exists in x-y plane with the stresses σ z ,τ xz ,τ yz equal to zero. Plane strain elements are used to represent a slice (of unit thickness) of a structure in which the strain components ε z , γ xz , γ yz are zero. This situation arises in the analysis of long dam, retaining wall and so on (Fig. 7.4). •

Plate and shell elements The basic proposition in plate and shell analyses is that the structure (Fig. 7.5) is thin in one dimension and therefore the following assumptions can be made: 1) The stress through the thickness ( σ z = 0 ) of the plate/shell is zero. 2) Material particles that are originally on a straight line perpendicular to the mid-surface of the plate/shell remain on a straight line during deformations. In the Kirchhoff theory, shear deformations are neglected and the straight line remains perpendicular to the mid-surface during deformations. In the Mindlin theory, shear deformations are included and therefore the line originally normal to the mid-surface does in general not remain perpendicular to the mid-surface during the deformations. In certain problems, the given body can not be represented as an assemblage of only one type of elements. In such cases, we may have to use two or more types of elements for idealisation. Examples of this would be the analysis of an aircraft wing or the analysis of a car body. 7.3

Size of elements

Size of elements influences the convergence of the solution directly and hence it has to be chosen with care. If the size of elements is small, the final solution is expected to be more accurate. However, we have to remember that the use of elements of smaller size will also mean more computational time. Sometimes, we may have to use elements of different sizes in the same body. For example, in the case of stress analysis of a plate with a hole, elements of different sizes have to be used. Size of elements has to be very small near the

55

z

mid-surface

FE

y x h

z

z

mid-surface

x

y

x

y

h

σz =0 All other stress components are nonzero.

Fig. 7.5 Plate and shell structures. hole (where stress concentration is expected) compared to far away places. In general, whenever steep gradients of the field variable are expected, we have to use a finer mesh in those regions. Another characteristic related to the size of elements, which affects the finite element solution, is the aspect ratio of elements. The aspect ratio describes the shape of element in the assemblage of elements. For two-dimensional elements, the aspect ratio is taken as the ratio of the largest dimension of the element to the smallest dimension. Elements with an aspect ratio of nearly unity generally yield best results. 7.4

Location of nodes

If the body has no abrupt changes in geometry, material properties and external conditions (like load, temperature, etc.), the body can be divided into equal subdivisions and hence the spacing of nodes can be uniform. On the other hand, if there are any discontinuities in the problem, nodes have to be introduced obviously at these discontinuities, as shown in Figure 7.5, where (a) and (b) - discontinuity in loading, (c) discontinuity in geometry, (d) - discontinuity in material properties, (e) - discontinuity in material. 7.5

Number of elements

The number of elements to be chosen for idealisation is related to the accuracy desired, size of elements and number of degrees of freedom involved (dimension of problem). An increase in the number of elements generally means more accurate results (Fig. 7.7). However, the use of large number of elements involves large number of degrees of freedom and we may not be able to store the resulting matrices in the available computer memory.

56

q P

aluminium steel node

nodal line

a) Concentrated load on a beam.

q1

d) Discontinuity in material properties.

q2

nodes

q

q

node

b) Abrupt change in the distributed load.

e) A cracked plate under loading.

node c) Abrupt change in the cross section of a beam. Fig. 7.6 Location of nodes at discontinuities.

exact solution no significant improvement beyond N 0

FEM

N0

number of elements

Fig. 7.7 Convergence in results.

57

y (v)

x(u)

v =0

v=0 u=0

Fig. 7.8 A plate with a hole with symmetric geometry and loading. 7.6

Simplifications afforded by the physical configuration of the body

If configuration of the body, as well as the external conditions, is symmetric, we may consider only half or quarter of the body for finite element idealisation (Fig. 7.8). The symmetry conditions, however, have to be incorporated in the solution procedure. 7.7

Finite representation of infinite bodies

In some cases, like in the case of analysis of foundations and semi-infinite bodies, the boundaries are not clearly defined. Fortunately it is not necessary to idealise the infinite body. So, in the case of analysis of foundation, the effect of loading decreases gradually with increasing distance from the point of loading and we can consider only the continuum in which the loading is expected to have significant effect. In this case, the boundary conditions for this finite body have to be incorporated in the solution. In the present example, the semi-infinite soil has been simulated considering only a finite portion of the soil (Fig. 7.9). In some applications, the determination of size of the finite domain may pose a problem. In such cases, one can use infinite elements for the modelling.

P

Footing

P

u=v=0 or v=0

y (v )

Semi-infinite soil

u=0

u=0 x (u )

Fig. 7.9 A foundation under concentrated load.

58

7.8

Node numbering scheme

Since most of the matrices involved in the finite element analysis are symmetric and banded, the required computer storage can be considerably reduced storing only the elements involved in the half bandwidth instead of storing the whole matrix. The bandwidth of the assemblage matrix depends on the node numbering scheme and the number of degrees of freedom considered per node. If we can minimise the bandwidth, the storage requirements, as well as solution time, can also be minimised. Since the number of degrees of freedom per node is generally fixed for any given type of problem, the bandwidth can be minimised using a proper node numbering scheme. As an example, consider a three-bay frame with rigid joints, 20 storeys high (Fig. 7.10). A shorter bandwidth (B) can be obtained numbering the nodes across the shortest dimension of the body. 7.9

Automatic mesh generation

For large systems, the procedure of node numbering becomes nearly impossible. Hence automatic mesh generation algorithms capable of discretizing any geometry into an efficient finite element mesh without user intervention are applied. Most of the automatic bandwidth renumbering schemes permits an arbitrary numbering scheme initially. Then the nodes are renumbered through an algorithm to reduce the bandwidth of the system equations. After the system equations are solved, the node numbers are often converted back into the original ones.

1

2

3

4

1

21 41 61

5

6

7

8

2

22 42 62

3

23 43 63

20

40 60 80

77

78 79

80

B=15 a) along the shorter dimension

B=63 b) along the longer dimension

Fig. 7.10 Node numbering scheme.

59

Chapter 8

8.1

Finite element program packages

Introduction

The general applicability of the finite element method makes it a powerful and universal tool for a wide range of problems. Hence a number of computer program packages have been developed for the solution of a variety of structural and solid mechanics problems. Among more widely used packages are ANSYS, NASTRAN, ADINA, LS-DYNA, MARC, SAP, COSMOS, ABAQUS, NISA. Each finite element program package consists from three parts: • programs for preparation and control of the initial data, • programs for solution of the finite element problem, • programs for processing of the results. Let’s consider now some general features of a more widely applied finite element program - ANSYS. The ANSYS program is a computer program for the finite element analysis and design. The ANSYS program is a general-purpose program, meaning that you can use it for almost any type of finite element analysis in virtually and industry - automobiles, aerospace, railways, machinery, electronics, sporting goods, power generation, power transmission and biomechanics, to mention just a few. “General purpose” also refers to the fact that the program can be used in all disciplines of engineering - structural, mechanical, electrical, electromagnetic, electronic, thermal, fluid and biomedical. The ANSYS program is also used as an educational tool at universities. ANSYS software is available on many types of computers including PC and workstations. Several operating systems are supported. The procedure for a typical ANSYS analysis can be divided into three distinct steps: • build the model, • apply loads and obtain the solution, • review the results. 8.2

Build the model

“Build the model” is probably the most time-consuming portion of the analysis. In this step, you specify the job-name and analysis title and then use pre-processor (PREP 7) to define the element types, element real constants, material properties, and the model geometry. The ANSYS element library contains over 80 different element types. Each element type is identified by unique number and prefix that identifies the element category:

60

BEAM4, SOLID96, PIPE16, etc. The following categories are available: BEAM, COMBINation, CONTACT, FLUID, HYPERelastic, INFINite, LINK, MASS, MATRIX, PIPE, PLANE, SHELL, SOLID, SOURCe, SURFace, USER, and VISCOelastic (or viscoplastic). The element type determines, among other things, the degree-of-freedom set (which implies the discipline - structural, thermal, magnetic, electric, fluid, or coupled-field), the characteristic shape of the element (line, quadrilateral, brick, etc.), and whether the element lies in 2-D space or 3-D space. BEAM4, for example, has 6 structural degrees-of-freedom (UX, UY, UZ, ROTX, ROTY, ROTZ), is a line element and can be modelled in 3-D space. Element real constants are properties that are specific to a given element type, such as cross-sectional properties of a beam element. For example, real constants for BEAM3, the 2-D beam element, are area (AREA), moment of inertia (IZZ), height (HEIGHT), shear deflection constant (SHEARZ), initial strain (ISTRN), and added mass per unit length (ADDMAS). Material properties are required for most element types. Depending on the application, material properties may be linear, nonlinear, and/or anisotropic. The main objective of the step “Creating the model geometry” is to generate a finite element model - nodes and elements - that adequately describes the model geometry. There are two methods to create the finite element model: solid modelling and direct generation. With “solid modelling”, you describe the geometric boundaries of your model and then instruct the ANSYS program to automatically mesh the geometry with nodes and elements. You can control the size and shape of the elements that the program creates. With “direct generation”, you “manually” define the location of each node and the connectivity of each element. Several convenience operations, such as copying patterns of existing nodes and elements, symmetry reflection, etc. are available. 8.3

Apply loads and obtain the solution

In this step, you use SOLUTION menu to define the analysis type and analysis options, apply loads, specify load step options, and initiate the finite element solution. The analysis type is chosen based on the loading conditions and the response you wish to calculate. For example, if natural frequencies and mode shapes to be calculated, you would choose a modal analysis. The following analysis types are available in the ANSYS program: static, transient, harmonic, modal, spectrum, buckling, and substructuring. Not all analysis types are valid for all disciplines. Modal analysis, for example, is not valid for a thermal model. Analysis options allow you to customise the analysis type. The word “loads” as used in the ANSYS program includes boundary conditions as well as other externally and internally applied loads. Loads in the ANSYS program are divided into six categories: • DOF constraints, • forces, • surface loads, • body loads, • inertia loads, • coupled-field loads. Most of these loads can be applied either on the solid model (keypoints, lines and areas) or the finite element model (nodes and elements).

61

Load step options are options that can be changed from load step to load step, such as number of substeps, time at the end of a load step, and output controls. A load step is simply a configuration of loads for which you obtain a solution. In a structural analysis, for example, you may apply wind loads in one load step and gravity in a second load step. Load steps are also useful in dividing a transient load history curve into several segments. Substeps are incremental steps taken within a load step. They are mainly used for accuracy and convergence purposes in transient and nonlinear analyses. Substeps are also known as time steps - steps taken over a period of time. After SOLVE command, the ANSYS program takes model and loading information from the database and calculates the results. Results are written to the results file (Jobname.RST, Jobname.RTH, or Jobname.RMG) and also to the database. The difference is that only one set of results can reside in the database at one time, whereas all sets of results (for all substeps) can be written to the results file. 8.4

Review the results

Once the solution has been calculated, you can use the ANSYS postprocessors to review the results. Two postprocessors are available: POST 1 and POST 26. POST 1, the general postprocessor, is used to review results at one substep (time step) over the entire model. You can obtain contour displays, deformed shapes, and tabular listings to review and interpret the results of the analysis. Many other capabilities are available in POST 1, including error estimation, load case combinations, calculations among results data, and path operations. POST 26, the time history postprocessor, is used to review results at specific points in the model over all time steps. You can obtain graph plots of results data versus time (or frequency) and tabular listings. Other POST 26 capabilities include arithmetic calculations, and complex algebra.

62

Literature

1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11.

Bathe K.-J. and Wilson E. L. Numerical Methods in Finite Element Analysis. – Prentice-Hall, Inc., 1976. Sigerlind L. J. Applied Finite Element Analysis. – John Wiley and Sons, Inc., 1976. Варвак П. М., Бузун И. М., Городецкий А. С., Пискунов В. Г., Толокнов Ю. Н. Метод конечных элементов. – Головное издательство издательского объединения «Вища школа»: Киев, 1981. Rao S. S. The Finite Element Method in Engineering. – Pergamon Press, 1989. Ritz W. Über eine Neue Metode zur Lösung gewisser Variationsprobleme der Matematischen Physik // J. Reine Angew. Math., 1909, Vol. 135, P. 1-61. Courant R. Variational methods for the solution of problems of equilibrium and vibrations // Bulletin of the American Mathematical Society, 1943, Vol. 49, P. 1-23. Clough R. W. The finite element method in plane stress analysis. // Proc. American Society of Civil Engineers (2nd Conference on Electronic Computation, Pitsburg, Pennsylvania), 1960, Vol. 23, P. 345-378. Argyris J. H. Energy theorems and structural analysis // Aircraft Engineering, 1954, Vol. 26, Part 1 (Oct. – Nov.), 1955, Vol. 27, Part 2 (Feb. – May). Turner M. J., Clough R. W., Martin H. C. and Topp L. J. Stiffness and deflection analysis of complex structures // Journal of Aeronautical Science, 1956, Vol. 23, No. 9, P. 805-824. Hrennikov A. Solution of problems in elasticity by the frame work method // Journal of Applied Mechanics, 1941, Vol. 8, P. 169-175. Zienkiewicz O. C. and Cheung Y. K. The Finite Element Method in Structural and Continuum Mechanics. – McGraw-Hill: London, 1967.

63

APPENDIX

A typical ANSYS static analysis

The goal of this example is to model the problem as shown in Figure A1. This is a 3D plate model where ANSYS general shell elements will be used to predict the displacement and stress behaviour of the plate subjected to concentrated loads from one side. The structure is clamped at the left end so that no translations or rotations are allowed there.

(2,135,0)

2,135,1)

1i

(2,135,3) (2,135,4)

2i n 1i

n 2i

45

y(α)

n

n

(0,0,0)

(1,135,1)

(1,135,3)

0

(0,0,4)

F 10

F

in

F

(10,0,0)

F

4i

h=0.1 in F=0.25 lbs

(10,0,4)

Fig. A1 3D plate model. Material properties: a) Isotropic. b) Young`s modulus, E=30e6 psi. c) Poisson`s ratio, υ=0.3.

64

n

X(R)

GEOMETRIC MODELLING STEP 1:

To build geometry of the model more easily, we change the coordinate system - from the default decart to cylindrical.

ANSYS Utility Menu WorkPlane > Change Active CS to > Global Cylindrical STEP 2:

Since the geometry of the model is 3D, we change the display from top view to isonometric view.

ANSYS Utility Menu Plot Ctrls > Pan,Zoom,Rotate… > Iso STEP 3:

Create ten keypoints allocated in corners of the plate by opening the window, defining the keypoint number and coordinate values, where X coordinate is radius, Y coordinate is angle, Z coordinate is height.

ANSYS Main Menu Preprocessor > -Modelling-Create > Keypoints > In Active CS… Keypoint number X,Y,Z Location in active CS Apply

1 0

0

0

Keypoint number X,Y,Z Location in active CS Apply

2 0

0

4

etc. OK STEP 4:

Create thirteen lines by joining two keypoints. You have to specify start and end of the line by picking with mouse to keypoints. Keypoints must be joined as shown in Figure A2.

ANSYS Main Menu Preprocessor > -Modelling-Create > Lines > Straight Line STEP 5:

Model four new areas by picking with mouse corresponding four lines as shown in Figure A2.

ANSYS Main Menu Preprocessor > -Modelling-Create > -Areas-Arbitrary > By Lines

65

10 9 L12 L11 A4 L13 8 5 L6 L7 L9 1 L10 7 A2 A3 L5 L8 L1 6

2

L4 A1 L2 4 L3 3

Fig. A2 Keypoints, lines and areas of the geometric model.

DESCRIPTION OF FINITE ELEMENTS STEP 6:

For the plate model we choose four-node shell element.

Ansys Main Menu Preprocessor > Element Type > Add/Edit/Delete…> Add…> Structural Shell Elastic 4 node 63 Ok Close STEP 7:

Plate thickness must be also definite.

Ansys Main Menu Preprocessor > Real Constants... > Add... > OK Shell thickness at node I, J, K, L = 0.1 Ok Close

66

MATERIAL MODELLING STEP 8:

Add material data for the isotropic material: Young`s modulus E=30e6 psi and Poisson`s ratio υ=0.3.

Ansys Main Menu Preprocessor > Material Props > Constant-Isotropic... OK Young’s modulus Poisson’s ratio (major) OK

EX=30e6 PRXY=0.3

FINITE ELEMENTS MESHING STEP 9:

For smooth mesh we choose all lines divided into three elements excluding the division of two long sides divided into eight finite elements. After choosing size of the mesh, we allow computer to perform mesh of all marked areas. The finite element mesh is presented in Figure A3.

Ansys Main Menu Preprocessor > -Meshing-Size Cntrls > Picked Lines No. of element divisions = 3 No. of element divisions = 8 OK Ansys Main Menu Preprocessor > -Meshing-Mesh > -Areas-Free OK

APPLICATION OF LOADS AND BOUNDARY CONDITIONS STEP 10: Mark the right-side edge nodes and input the value of applied force -0.25 lbs. The applied load is shown in Figure A3. Ansys Main Menu Solution > -Loads-Apply > -Structural-Force/Moment > On Nodes Direction of force/mom Apply As Value OK

FY Constant value -0.25

67

STEP 11: Apply boundary conditions to the left-side edge nodes defining no deformation and rotation for all nodes. The boundary conditions are shown in Figure A3. Ansys Main Menu Solution > -Loads-Apply > -Structural-Displacement > On Nodes DOF’s to be constrained OK

All DOF

SOLUTION STEP 12: Analysis type is static. Now solution can be started. When the window “Solution is done” appears, solution is completed. Ansys Main Menu Solution > Analysis type-New Analysis... > STATIC

OK Ansys Main Menu Solution > -Solve-Current LS OK

Fig. A3 Finite element mesh, applied load and boundary conditions.

68

ANALYSIS OF RESULTS STEP 13: After static analysis you can plot a deformation shape of the plate (Fig. A4). Ansys Main Menu General Postproc > Plot Results > Deformed Shape... Def + undeformed STEP 14: Besides deformation, stress state of the plate (Fig. A5) can be calculated by von Mises theory. Ansys Main Menu General Postproc > Plot Results > -Contour Plot-Nodal Solu... Item to be contoured OK

Stress

von Mises SEQV

Fig. A4 Deformation shape of the plate.

69

Fig. A5 Stress state of the plate.

70