Copy Of Finite Element Formulation

Dr. B. S. SARMA

Overview  semi-discretization

-Variational Equations of Linear Solid Mechanics  Finite Element Discretization  Analytical Solution Procedures -Modal Superposition Method -Eigen value economizer Method -Free Response , Damped Eigen values  Forced Periodic Response  Transient Response -Frequency Response Method -Modal Decomposition Method

 With FEM, many dynamic and nonlinear numerical simulations

are now possible within the realm of engineering analysis and considerable effort has been directed towards the area of structural dynamics.  The most common approach to structural dynamics is the one in which the spatial and temporal solutions are performed independently. This is called SEMI-DISCRETIZATION.  The governing partial differential equations are first discretized in space yielding a coupled system of ordinary differential equations in time which are solved using various integration methods.  Finite Element Method is an example of SEMIDISCRETIZATION approach for a multi-dimensional linear continuum.



In the FINITE ELEMENT FORMULATION, one starts with the equation of a continuum and their variational forms and use the finite element method to construct the discrete equations of motion.



The continuum formulation can be extended to include structural mechanics theories such as Beams, Plates and Shells.



In the analysis of small deformations of a solid continuum, the governing equations are obtained from the MOMENTUM CONSERVATION LAW and they are referred to as EQUATIONS OF MOTION. In the absence of body forces they are written as

τij,j = ρüi

in Ω τij is the stress tensor

ρ is the density u is the displacement and Ω represents the continuum domain.

(1)

 The INFINITESIMAL STRAIN is related to displacements

by

εij

= ½ (ui,j + u j,i )

(2)

 The constitutive law relating stress to strain for elastic material

is given by

τij where Cijkl

= Cijkl εkl is the elastic tensor.

(3)

 The continuum has a boundary Γ which can be divided into

Γu where displacements are prescribed Γτ where stresses are prescribed.

These boundary conditions are written as ui = ui on Γu and τij nj = τi on Γτ *

*

and

(4)

 In addition we have initial conditions

u = u and ů = ů i(0)

0i

i(0)

0i

(5)

 Equations (1) (4) and (5) represent the strong form of the

initial/ boundary value problem of Elasto-Dynamics.  To facilitate the finite element discretization, the weak

form of the Variational statement corresponding to these equations is needed.



Let a test function w be defined such that it is a function of space only and is identically zero on the boundary Γu. Then the variational statement is as follows: Given the prescribed tractions τi* , initial displacements uoi and velocity ůoi , find ui (t) such that for all wi int

M (w,u) + F (w,u) = F

ext

Where M(w,u) = ∫Ω wi ρ üi dΩ int

F

ext

F

(w, τ)

(6) (7)

= ∫Ω wi,j τij dΩ

(8)

= ∫Ω wi τi* dΩ

(9)

Equations (7) to (9) correspond to the work done by the inertial, internal and external forces respectively. 

It is assumed that ui is chosen such that the essential boundary conditions are exactly satisfied.



For a linear elastic case equation (8) can alternatively be written by using the equation (3) and the symmetry of the stress tensor as int

F

= ∫Ωw(i,j) Cijkl u(k,l) dΩ

(10)

 



In Finite element semi-discretization, the domain Ω is subdivided into elements Ωe The displacement field in each element is approximated by ui(x,t) = NI (x) diI (t) (11) where the shape function NI (x) are functions of space alone and diI (t) include the time dependence. Introducing eqn.(11) in matrix form

u = N de

(12)

where de are the displacements local to element e. 

The discrete matrix form of the gradient u (i,j) can then be written as s u = B d e (13) where B contains the appropriate derivatives of Ni

Considering that the same shape functions are used for both ui and wi( as in the Galerkin approximation), the semi-discretization of eqn (6) yields int

δd (M

ext

+F -F ) = 0

T

(14)

where δd is a vector of nodal parameters resulting from the

interpolation. Global mass matrix M; internal force vector Fint ; external force vector F ext and nodal displacement vector d are assembled from element contributions where e

M

=AM

int

F F d

ext

e

e

int .e

e

ext.e

=A F

=A F

e e

= Ad

(15) (16) (17) (18)

 The element matrices are defined using Eqns.(11), (12) and (13)to be e

Me = int,e

F

ext,e

F

[ M ijIJ] = ∫Ωe δij Ni ρ Nj dΩ

= ∫Ωe Bt τ dΩ = {FiI

ext,e

}

(19) (20)

= > ∫Γeτ Ni τ*i dΓ

(21)

 Since Eqn. (14) must hold for any δd, the equation of motion becomes M

int

+F =

F

ext

(22)

 For a linear elastic material, the equation of motion can alternatively be written by using Eqn (10) as M

+ Kd =

ext

F

(23)

where K is the global stiffness matrix assembled from element e e contribution i .e, K = AK (24) and

K

= ∫Ωe Bt DB dΩ

(25)

 To include viscous damping in a consistent manner one can consider damping matrices of the form (Rayleigh’s Damping ) C = aM + bK



where the coefficients ‘a’ and ‘b’ control the amount of damping proportionality with mass and stiffness respectively ‘a’ and ‘b’ are determined experimentally. Now the damped equation of motion is given by M

+ C

+ Kd

ext

=F

(27)

 To account for mass proportional damping effect, Eqn (1) to be replaced by τij,j = ρüi + aρů in Ω (28) and for the stiffness proportional damping, Eqn (3) is to be augmented as τij = Cijkl (εkl + b ůck,l) (29)  Since M and C are derived from a variational form,they are referred to

CONSISTENT AND LUMPED MASS MODELS 









Consistent Mass Matrices obtained from Eqns.(15) and (19)are one form of mass matrices used in current FE S/W packages. In the context of Direct Time integration solutions Lumped Mass matrices are often considered. They are diagonal matrices with offdiagonal terms as zeros. Another alternative is the so called higher order matrices which are a combination of both consistent and lumped matrices. Implicit methods can operate on larger time steps but require the solution of equations at each time step. Consistent mass matrices have the same sparcity and bandwidth as that of the stiffness matrix. So the computational work required for the solution is the same for consistent or lumped matrices. Explicit time integration schemes require a diagonal mass matrix to be truly explicit. Therefore lumped mass models are used.

MASS LUMPING SCHEMES 

For Low order elements in continuous formulation the simplest procedure is to assign an equal fraction of total mass of each element to each node.

Me = [ MeijIJ ] = δij δIJ / nen ∫ e ρ dΩ Ω

where nen is the number of nodes per element. 

Another simple procedure is the row sum technique n en

e [ MeijIJ ] = Σ J=1 M ijIJ

The same can be written for continuum elements as

[ MeijIJ ] = δij ∫ e ρNI dΩ

n



Ω

en The above lumped mass scheme conserves mass in that

Σ MeijIJ = δij ρ Ωe

J=1







For higher order elements nodal wuadrature method is normally employed. Here the quadrature points coincides with the nodes resulting in a lumped mass. Labatto Quadrature rules such as Trapezoidal rule and Simpsons rule are used to diagonalize the mass amtrex. By definition of shape functions, Numerical quadrature finally becomes

[ MeijIJ ] = δij ρ (ξ1)J(ξ1)wi 



A draw back of nodal quadrature is that, in some cases, zero or negative masses are obtained. A Scheme that always results in positive lumped masses suggested is

[ MeijIJ ] = α δij ∫ ρNi2 dΩ Ω

Where α is so chosen that the total element mass is conserved.

ANALYTICAL SOLUTION PROCEDURES 







As a result of semi-discretization, many time – dependent problems can be reduced to a system of ordinary differential equations of the characteristic form given by ext M + C + Kd = F (30) Systems of linear ordinary differential equations with constant coefficients can always be solved in principle analytically without the introduction of additional approximations. Here we shall deal with ext a. Determination of free vibration response (i.e F (t)=0) ext b. Determination of periodic response (i.e (i.e F (t) is periodic) ext c. Determination of transient response (i.e (i.e F (t) is arbitrary) In the first two cases, initial conditions of the systems are of no importance and a general solution is sought.

MODAL ANALYSIS 

Regardless of the type of structure analyzed or semi-discretization method used, Eqn (30) represents a standard form of the discrete equation of motion for linear structural dynamics



The coupled system of ordinary differential equations possesses certain attributes such as M, C, K are symmetric and M as positive definite and K are semi-definite matrices.



Two Techniques to solve these equations, for the given initial conditions are * Direct Time Integration and Modal Superposition



In Direct Time Integration, the solution is marched through time step by step starting with the initial conditions and progressing to find the final solution.

 With Modal Superposition, the solution for all time is obtained simultaneously by a superposition of contributions from the natural modes of the system.



Modal Analysis is a procedure to transform the coupled system of ordinary differential equations described by Eqn (27) into a set of uncoupled equations representing the natural modes of the system.



Basic understanding of Modal Analysis is essential even if Direct Time Integration is used as a solution scheme.

THE MODAL SUPERPOSITION METHOD 

Consider the homogeneous, undamped equation of motion M + kd = 0 (31) Assume possible solutions to Eqn (31) of the form d (t) = Φ sin( ωt + θ) (32) where ω is the natural frequency and θ is a constant phase angle



Differentiating Eqn (32) approximately and substituting it into Eqn (31), we get (K Φ - ω2M Φ) sin ( ωt + θ) = 0 (33) i.e [K – ω2M] Φ = 0



(34)

This is an Eigen value problem and its non- trivial solution exists if Det |(K-ω2M)| = 0 (35)



When expanded out, Eqn (35) becomes an neqth degree algebraic equation called the characteristic equation. The neq roots of this equation are ω1,2ω22ω32 . . . ωn2 and these frequencies are generally placed in the order ω12 < ω22 <ω32 …. ωn2



For each ωr there exists a solution to Eqn (35) as Φr . The vector is called the mode shape for mode r



For convenience, Φ is normalized with respect to the mass matrix Φr

T

M Φr = 1

(36)

T

Pre multiplying Eqn(34) by Φr we have for the rth mode Φr

T

KΦr = ω

2 r

(37)



From the property of modal orthogonality, we have T

Φr M Φs = 0 for r ≠ s

(38)

T









Φr K Φs = 0 for r ≠ s (39) To prove the above statements, we can write, from Eqn (34) for any mode ωr2 MΦr = K Φr (39a) ωs2 MΦs = K Φs (39b) T

T

Pre-multiplying the first by Φs and the second by Φr and then subtracting, we get (Noting that M is symmetry) T 2 2 (ωr - ωs )(Φr M Φs) = 0 (39c) And if ωr ≠ ωs, the orthogonality condition for matrix M is proved. Immediately from this, orthogonality of the matrix K follows. Finally the modal matrix Φ defined as Φ = [Φ1, Φ2, Φ3, Φ4.. Φn eq] (40)



The coupled system of equations can now be transformed to a set of uncoupled equations of modal components,



Expanding the solution d(t) as the sum d(t) = Φ1Z1(t) + Φ2Z2(t) + … Φn eqZ n eq(t) = Φ Z(t)

(41)



where z(t) are the normal coordinates that incorporate the time dependence.



Substituting Eqn (41) into the Eqn (23) and pre-multiplying by T Φr , we have representing each mode of vibration T

T

T ext

Φr M Φ Z”(t) + Φr K Φ Z(t) = Φr F

(42)









Using orthogonality conditions Eqn (38,39) and definitions of Eqn(36,37) on Eqn (42), we obtain the uncoupled set of ordinary differential equations in time Z”(t) + ω r2 Zr(t) = Fr for r = 1,2,3….. n eqns

(43)

where

(44)

Fr = ΦrT Fext

In the modal superposition method, the global system of Eqns (23) is transformed into the set of single degree of freedom problems (43) by the normal coordination transformation. Each of the single degree of freedom problems is solved using a suitable transformation of initial conditions. Finally the solution in physical coordinates is found by the normal coordinate transformation which is a superposition of various modal equations.



If Rayleigh damping Eqn (26) is used, the same normal mode coordinate transformation may be applied to Eqn (27) to obtain a set of damped single degree of freedom equations of he form Zr”(t) + 2 ζr ωr Zr’(t) + ωr2 Zr(t) = fr , for r =1,2,3…n eq



(45)

The modal fraction of critical damping ζr is determined by the mass and stiffness proportionality coefficients a, b and ωr by the relation ζr = (a/ ωr + b ωr)/2

(46)

EIGEN VALUE ECONOMIZER METHOD: 

Computer effort of determining the Eigen values and Eigen vectors is larger by an order of magnitude than the solution of an equivalent static problem.



Fortunately, reasonably good Eigen values can be determined with fewer degrees of freedom than needed for static solutions.



If a rather fine sub-division is used in the dynamic analysis, a number of degrees of freedom can be eliminated and the mass and damping effects can be lumped at a reduced number ofnodal parameter. [similar to the sub-structure analysis used in static problems]



Let the nodal vector d be divided into two parts d= { d } (47) d Assume that the displacements ds depend in some unique way on the displacement dm. i.e ds is called master and dm is called slave variables. If dm = T ds (48) We have I d = [ T ] dm => T* dm (49) s

m

Where T is the matrix specifying the dependence. Now the dynamic equation of the whole system M + kd = 0 (50) can be reduced by applying the constraint on the deformation freedom implied in Eq (49)

Pre-multiplying Eqn (50) by T*T (principles of transformation) after substituting in Eqn (49), we get M*

m

+ K* dm = 0

Where K = T*T K T and M* = T*T M T Eqn (51) involves smaller number of variables. 

(51)

How to determine the relationship between the slave and the master deflections:- It is reasonable to assume that the general pattern of deformation will follow that which would be obtained by imposing displacements dm on an otherwise unloaded structure in static condition. - Thus partitioning with = = 0, the standard equation can be written asKss Ksm ds 0 Kd = [ Ksm Kmm ] {

dm }

= {

fm

}

(52)



As the slave loads are unloaded, we can get Kss ds + Ksm dm = 0 or,

ds = - Kss -1 Ksm dm i.e T = - Kss -1Ksm



(53) (54)

The best choice of the master nodes - The obvious choice is to eliminate first the nodes with little or no mass attached. - In an automatic procedure, the ratio of diagonal stiffness and mass terms be first calculated (kii / Mii) and the nodes with the largest values of this ratio be first eliminated . -One drawback of the condensation is that the mass matrix in the reduced problem is never diagonal .

FREE RESPONSE, DAMPED EIGENVALUES: For Free vibration conditions, the standard equation (27) is written as M

+c

+ K= 0

And substituting d= Φ e

(55) άt

The characteristic equation becomes 2

[(α M + α C + K) {Φ }] = 0 Where α and Φ in general are found to be complex

(56)

FORCED PERIODIC RESPONSE 

 



If the forcing term in the standard equation (30) is periodic or more generally if we can express it as αt f = fe (57) Then a general solution can be written as αt d = Φe (58) Where d is a complex value ( i.e α = α1 + i α2 ) - (30), we get Substituting the above equation in Eqn [α2 M + α c + k ]Φ => D Φ= - f (59) It is no longer an Eigen value problem. But it can be solved by inverting the matrix D i.e -1 Φ = -D f (60) Eqn(60) is precisely of the same form as that used for static problems. But, now it has to be determined in terms of complex quantities.



The computation can be arranged directly, noting that αt



e- = e- (cos- α2t + i sin α2t) f = f1 + if2 Φ = Φ1 + i Φ2 in which α1 , α 2 , f 1 , f 2 , Φ1, Φ2 are real quantities Inserting the above quantities in equation (59) Φ1 ((α 2 –α 2 )M + α C +K) (-2α α M -α c) 1

[ 

 



αt

2

(2α1α2 M + α2C)

1

1

2

1

((α12 – α22 )M + α1C + K )

(61)

f1

Φ2

] { } = -{ f2 } 62)

Eqn(62) form a system in which all quantities are real and from which the response to any periodic input can be obtained by the direct solution. The system is no longer positive definite (It is skew symmetric) With periodic Inputs , the solution after an initial transient is not sensitive to the initial conditions and this guessed solution represents the final solution. It is valid for problems of dynamic structural response

TRANSIENT RESPONSE BY ANALYTICAL PROCEDURES 

Steady state general solutions neither took care of any account of initial conditions nor of the non-periodic form of the forcing terms



The response of dynamic system taking these features requires either a full-time discretization or special analytical procedures.



Two broad possibilites include * FREQUENCY RESPONSE PROCEDURE * THE MODAL ANALYSIS PROCEDURE

FREQUENCY RESPONSE PROCEDURE 





The response of the system to any forcing terms of the general periodic type or in particular, to a periodic forcing function f = f eiωt was already discussed. As a completely arbitrary forcing function can be represented by a Fourier series m f(t) = f0 + Σ (am cos ωmt + bm sinm ωmt) i=1 or in the limit exactly as Fourier Integral, the response to such an input can be obtained by a synthesis of a curve representing the response of any quantity. e.g: Displacement at a particular point to all frequencies ranging from zero to infinity. In fact only a limited number of such forcing frequencies need to be considered and the results can be efficiently synthesized by fast Fourier transform techniques The technique of frequency response is adapted to problems where damping matrix C is of an arbitrary specified form.

MODAL DECOMPOSTION ANALYSIS 





This procedure is the most important and is widely used in practice. It provides insight into the behavior of the whole system, where strictly numerical procedures are used. Starting with general problem of equation(2 7) i.e M + C + Kd = F(t) (64) where f(t) is an arbitrary function of time. The general solution for the free response is of the form m αt

αit

d = Φ e i.e Σi=1Φi e (65) Where di are the Eigen values and Φi are Eigen vectors.  For Forced response, assume that the solution can be written in a linear combination of modes as d = Σ Φi yi = [ Φ1 Φ2 Φ3 … Φn] y

(66)

Where the scalar mode participation factor Y(t) is a function of





Such an assumed solution clearly gives the proportions of each mode occurring in the response. This assumption does not present any restriction since all the modes are linearly independent vectors Substituting Eqn (66) into equation(64) and pre-multiplying the result by Φit (i = 1,2,3,4.. n), then the resulting equations are simply a set of scalar independent equations.

miy”i + Ciy’ + ki yi = fi

(67)

mi = Φit M Φi => 1 (if modes normalized) Ci =

Φit C Φi ; Ki = Φit C Φi

and fi = Φit f(t).

As for true Eigen values Φi Φit MΦj = Φit C Φj = Φit K Φj = 0 (68) Each scalar equation of eqn (67) can be solved by usual



Assuming that the modes have been normalized so that mi = 1, the eqn (66) can be written in a standard form as of second degree: y” + 2ωci ’ yi’ + ωi2 yi = fi The general solution can be of the form -c ’ωi (t- τ) yi = ∫ fi e sin ωi(t-τ) dτ





(69)

(70)

Such integration can be carried out numerically and the modal response is obtained. In principle, superposition will result in the full transient response. Often a single calculation is carried out for each mode to determine the maximum responses and a suitable addition of these results is used as a standard procedure.

Thank You