Beam Dynamic Analysis by FEM Roberto Bernetti ii Contents I Intoduction to Math Software Package 1 1 Introduction to Math Language 3 1.1 Basic Matrix Manipulation . . . . . . . . . . . . . . . . 3 1.1.1 Matrix Definition . . . . . . . . . . . . . . . . . 4 1.1.2 The Range Operator . . . . . . . . . . . . . . . 6 1.1.3 Arithmetic Operators . . . . . . . . . . . . . . . 7 1.1.4 Matrix Operator . . . . . . . . . . . . . . . . . 8 1.2 Program Flow Control . . . . . . . . . . . . . . . . . . 11 1.3 Graphics Functions . . . . . . . . . . . . . . . . . . . . 13 1.4 File Scripting and Function File . . . . . . . . . . . . . 16 II Beam Dynamics 21 2 Dynamic Analysis of Beam 23 2.1 Kinematics of Beam . . . . . . . . . . . . . . . . . . . 23 2.1.1 Euler-Bernoulli Theory . . . . . . . . . . . . . . 24 2.2 Dynamics of Beam . . . . . . . . . . . . . . . . . . . . 28 2.2.1 Least Action Principle . . . . . . . . . . . . . . 28 2.2.2 Finite Element . . . . . . . . . . . . . . . . . . 34 2.2.3 Hermite Polynomials . . . . . . . . . . . . . . . 37 2.3 Numerical Exercise . . . . . . . . . . . . . . . . . . . . 38 2.3.1 Natural Frequncy . . . . . . . . . . . . . . . . . 40 2.3.2 Mode Shapes . . . . . . . . . . . . . . . . . . . 40 3 Mixed Finite Elements 49 3.1 Governing Equation . . . . . . . . . . . . . . . . . . . . 49 iii iv CONTENTS 3.1.1 Dynamic Equilibrium . . . . . . . . . . . . . . . 50 3.1.2 Thick Beam . . . . . . . . . . . . . . . . . . . . 52 3.2 Weighted Residual . . . . . . . . . . . . . . . . . . . . 54 3.2.1 Weak Formulation . . . . . . . . . . . . . . . . 55 3.2.2 Shape Functions . . . . . . . . . . . . . . . . . 59 3.2.3 Displacement Only Method . . . . . . . . . . . 62 3.3 Full Time Integration . . . . . . . . . . . . . . . . . . . 65 3.4 Numerical Exercise . . . . . . . . . . . . . . . . . . . . 68 3.4.1 Eigen Modes . . . . . . . . . . . . . . . . . . . 68 3.4.2 Time Integration . . . . . . . . . . . . . . . . . 70 4 Spatial Frame Applications 77 4.1 The column behaviour . . . . . . . . . . . . . . . . . . 77 4.2 Matrix Rotation . . . . . . . . . . . . . . . . . . . . . . 80 4.2.1 Spatial Behaviour . . . . . . . . . . . . . . . . . 81 4.2.2 Local to Global System . . . . . . . . . . . . . 84 4.3 Numerical Exercise . . . . . . . . . . . . . . . . . . . . 90 Preface The use of Finite Element Method (FEM hereon) together with the computational power of personal computers has widespread the abil- ity of engineers to solve real world problems, on the other hand com- puters and FEM modelling will help to increase the speed and the level of comprehension of the physic laws underlying problem. The aim of this book is to present both the theory and the numerical to allow the reader to understand all the formal and programming as- pect to perform a personal analysis on the influence of the various parameter that can influence the structure response. The book focus on the use of the Mixed Finite Element formulation that allow for an higher level of accuracy respect to the standard formulation, this due to the easy of use of higher polynomial degrees in the approximation function. The book is intended for graduate student all ready familiar with basic notion of dynamic analysis and Finite Element Method, here they can test and improve the ability to model and understand beam structure behavior in dynamics. The first part of the book give an introduction to the of high level programming language octave and freemat . In the second part both stationary and variational approach are presented in the derivation of FEM. In chapter 2 the standard, station- ary, approach to the FEM is presented together some octave/freemat functions to analyze the free vibrations problem. In chapter 3 the Mixed formulation is presented using a weighted residual approach that start from the dynamic equilibrium equation instead than from a stationary principle, in this way the reader will experience both procedure to deduce governing numerical equation. v vi CONTENTS The script by which it is possible to reproduce the plottings presents in the book are available at the following URL: http://www.robertobernetti.com/bin/view/Numerical In the recent years the FEM has been laid on a strong mathe- matical basis. The emphasis of the text is not on mathematical and numerical aspect anyway the necessary theory will be presented to both better understand the potential of the method and to interpret the results obtained. Part I Intoduction to Math Software Package 1 Chapter 1 Introduction to Math Language The language described in the present chapter (for easy hereon called MatLang) was firstly defined in the MATLAB(tm) program, then has been implemented in other mathematical packages like FreeMat and Octave. FreeMat and Octave are distribuited under the GNU Public License terms (GPL), detailed information about their distribution, development and initiation can be found on the following web sites: Matlab http://www.mathworks.com Octave http://www.octave.org FreeMat http://freemat.sourceforge.net In the following section a short introduction to MatLang will be pre- sented, enough to code the example showed in the book. 1.1 Basic Matrix Manipulation MatLang is an interpreted language, the output of the requested com- mand is assegned to a default variable ans. The following code is an example of the simplest operation: > > 2.1*1.23e11 ans = 2.5830e+11 3 4 1.1. Basic Matrix Manipulation > > MatLang can work as a simple calulator for you. Every numerical object in MatLang is considered a matrix, a scalar is a matrix of dimension 1x1. We introduce the first built-in function: size(x) return the dimension of the object x an example follow: > > size(ans) ans = 1 1 > > A wide number of built-in functions are available which size is an example of. The call is obtained by the function name followed by a variable numbers of arguments enclosed by parenthesys. Each func- tion returns an array object like size returns a two elements row vec- tor: > > sin(pi/2) ans = 1 > > The dimensions of the returned object depends on the function it- self and on the arguments dimension, this behaviour will be addressed later. We can use variable name for the numerical objcts, we will always get an answer from the program but this time the result of the numer- ical operation will be assigned to the variable insted to the built-in ans: > > x=2*3 x = 6 > > 1.1.1 Matrix Definition Matrix can be directly input by the user or can be definied using suited MatLang function. 1.1. Basic Matrix Manipulation 5 A matrix can be defined by input of its elements which shall be supplied surrounded by brackets, if all elements are on the same line rows shall be separated by a simi-colon. Examples: > > x=[ 1 3 ; -3 1] x = 1 3 -3 1 > > otherwise each row shall be terminated by the return key: > > x=[ 1 2 > -2 1 ] x = 1 2 -2 1 > > > Can be defined a row vector: > > x = [1 3 5 7 ] x = 1 3 5 7 > > or a column vector > > x = [1 ; 3; 5; 7 ] x = 1 3 5 7 > > 6 1.1. Basic Matrix Manipulation Some functions are available to directly build matrixes, a short list follows: zeros(m,n) creates a matrix which elements are all 0 with m rows and n colums ones(m,n) creates a matrix which elements are all 1 with m rows and n colums eye(m) create an identity matrix of size m A=diag(x,j) insert the x vector into the j position above (j > 0) or below (j < 0) the principal diagonal (j = 0 or absent) in the matrix A. If x is a matrix the operationis the extraction of the diagonal and its assigment to the vector y 1.1.2 The Range Operator One useful feature is the easy definition of equispaced set of number. The semicolon operator or range operator : defines the elements of an array or addresses a range of indexes in an array, the structure of the command is: start:inc:end in whichstart is the starting value,inc is the in- crement andend is the ending value Examples: This define a row vector containing the first 6 even num- bers: > > x=2:2:12 x = 2 4 6 8 10 12 > > while here it is possible to take the 2nd, 4th and 6th ele- ments of the vector: > > x(2:2:6) 1.1. Basic Matrix Manipulation 7 ans = 4 8 12 > > as it is possible to assign only a part of an arrray: > > x(2:2:6)=1:3 x = 0 1 0 2 0 3 > > The semicolon operator can be otherwise used as a place- holder to assign/reference an intere row/column: > > a(2,:)=2:2:8 a = 1 1 1 1 2 4 6 8 1 1 1 1 1 1 1 1 > > 1.1.3 Arithmetic Operators All the usual arithmetic operator are available: x+y is the matrix sum x.+y is the ement-by-element matrix sum, the matrix shall have the same dimension x*y is the matrix product and the columns of x shall be equal the rows of y x.*y is the ement-by-element matrix product x./y is the ement-by-element matrix divition 1./y is the ement-by-element inverse 8 1.1. Basic Matrix Manipulation 1.1.4 Matrix Operator Every data object is an array so it is very straightforward to use operator for matrix manipulation. Here are listed the more common: " " the transpose operator det(x) returns the determinant of the matrix x rank(x) returns the number of linerly indepentent rows/columns of x cond(x) returns the condition number of the matrix x inv(x) returns the inverse matrix of x eig(x) returns the eigenvalues and eigenvectors of matrix x; lu(x) returns the LU decomposition of the matrix x Example: The transpose of an array is very simple: > > x ans = 2 4 6 8 10 12 > > We can define a matrix: > >A=[ 1 -2 3 > -2 3 2 > 3 2 1] 1.1. Basic Matrix Manipulation 9 A = 1 -2 3 -2 3 2 3 2 1 > > and to calculate the determinant: > >det(A) ans = -56 > > It is possible now to verify if the matrix is invertible using the rank function: > > rank(A) ans = 3 > > calculate the condition number > >cond(A) ans = 1.3054 > > the inverse operator can be explicitly called: > >inv(A) ans = 0.017857 -0.142857 0.232143 -0.142857 0.142857 0.142857 0.232143 0.142857 0.017857 > > or by using the power operator: > >A^-1 10 1.1. Basic Matrix Manipulation ans = 0.017857 -0.142857 0.232143 -0.142857 0.142857 0.142857 0.232143 0.142857 0.017857 > > here a different calling syntax is presented, in fact if the call to the function is not in an assigment statment only the first object is loaded into the ans variable and showed in the output, otherwise if the call is in an assigment stat- ment the returned objects can be referenced using differ- ent variable names in the left hand side of the statement: > >[ V L ]=eig(A) V = -6.4464e-01 -7.0711e-01 -2.9057e-01 -4.1093e-01 -1.8489e-15 9.1167e-01 6.4464e-01 -7.0711e-01 2.9057e-01 L = -3.27492 0.00000 0.00000 0.00000 4.00000 0.00000 0.00000 0.00000 4.27492 > > which returns: V the eigenvectors in the columns of the matrix A L the eigenvalues in the diagonal element of the matrix A The same condition apply to the lu function, It is now possible to use the the LU decomposition to solve the linear system asociated to A: 1.2. Program Flow Control 11 > >[ L U ]=lu(A) L = 0.33333 -0.61538 1.00000 -0.66667 1.00000 0.00000 1.00000 0.00000 0.00000 U = 3.00000 2.00000 1.00000 0.00000 4.33333 2.66667 0.00000 0.00000 4.30769 > > The returned objects of the lu function are according to the calling sysntax: L the row-permuted lower triangular matrix U the uppper triangular marix 1.2 Program Flow Control To allow program flow control various type of logical operator are defined: == equality test operator ~= not equals operator (in octave equivalent also != <> ) <= less equal operator >= great equal operator < less then operator > great then operator & logical AND opeator ~ logical NOT operator 12 1.2. Program Flow Control | logical OR operator The comparison and logical opeators can be used in any logical espres- sion. Such expression can be used very effectivly as index parameter. Example > > x=2*pi/10:2*pi/10:2*pi; > > y=sin(x); > > y(pians = -5.8779e-01 -9.5106e-01 -9.5106e-01 -5.8779e-01 -2.4492e-16 > > The Logical expressions are used to control program flow in the lan- guage structures listed below: for loop control cycle with the following syntax forivar=range statements end the command setstatementsare executed for each value of the variableivar in the setrange, while loop conditional execution with syntax whiletest_condition statements end the command setstatementsare executed only if and till the conditiontest_conditionhas value TRUE before to enter the while loop if conditional statement with the following syntax: iftest_condition statements_1 1.3. Graphics Functions 13 else statements_2 end if thetest_conditionis TRUE (it is not necessary it is a logical variable but only that it is a number different from zero FALSE ) the set of command statements_1 will be executed otherwise the set statements_2 will be executed; the else clause can be omitted as can be substituted by an elseif statemennts folowed by a further condition. switch selective execution of statements with the following syntax switchexpression case value_1 statements_1 case value_2 statements_2 · · · · · · · · · · · · otherwise statements_d end according to the numerical or string values ofexpression the corresponding set of statements_iwill be executed, the clause otherwise enclose thestatements_dexecuted in the default case. 1.3 Graphics Functions Graphical data rendering is a key tool in numerical algorithm de- velopment, MatLang has very powerful functions to easy plotting 2 and 3 dimensional data. In the view of the following discussion the main graphic functions are listed (not all the functions here listed are in octave version 2.1.7 because they are in development phase; the 14 1.3. Graphics Functions same rendering can be obtained using gnuplot command from inside octave): plot(x,y,fmt,...) plot function graphs thex vector data versus theyvector data using the formatfmtto choose the line type and color as follow: - Solid line style : Dotted line style -. Dot-Dash-Dot-Dash line style Dashed line style r Color Red g Color Green b Color Blue k Color Black c Color Cyan m Color Magenta y Color Yellow variable numbers of data sets x,y,fmt are accepted contour(x,y,z) generates a contour plot (2D graph) of the spatial data stored inx,y,z matrixes, the matrixx andy con- tain respectivly the x and y coordinate of the point whose hight is stored inz surf(x,y,z) genearates a surface plot using data stored in thex,y,z matrixes, the matrixx andy contain respectivly the x and y coordinate of the point whose hight is stored inz subplot(m,n,act_win) divedes the plot window in mxn subregions (mrows andncolumns) and activate for subseequent plot commands the subregionacti_win axes(a) defines the range of x and y axes using the four element vector a = [xmin xmax ymin ymax] 1.3. Graphics Functions 15 title(s) insert thes string as graph title xlabel(s) insert thes string as label at x axes ylabel(s) insert thes string as label at y axes Examples the following example will draw the graph of an hermite polynomial as reported in figure 1.1 with axes labels and continuos blue line, the graph is then saved in png format using the command print. > > x=0:0.01:1; > > y=x.*(1-x).^2; > > plot(x,y, -b ) > > xlabel( x ); > > ylabel( y ); > > grid > > print( plotHermite.png ) the folowing example show the code to represent a spatial Hemite polynomial as a 3D surface, see figure 1.2: > > x=0:0.1:1; > > y=0:0.1:1; > > z=(1-3*x.^2+2*x.^3) *(1-3*y.^2+2*y.^3); > > surf(x,y,z) > > view(3) > > print( HermiteSurf.png ) > > In the following example the same polynomial surface will be showed using a contouring plot getting the drawing in figure : > > x=0:0.01:1; 16 1.4. File Scripting and Function File Figure 1.1: Graph of an Hermite polynimial > > y=0:0.01:1; > > z=(1-3*x.^2+2*x.^3) *(1-3*y.^2+2*x.^3); > > contour(x,y,z) here a little bit more sophisticated use of the grafics with handler, colorbar and contour line labels: > > clabel(h,[1:-0.2:0]) > > h=contour(x,y,z); > > colorbar > > clabel(h,[1:-0.2:0]) > > print( HermiteContLabel.png ) 1.4 File Scripting and Function File The sequence of commands to performe a difinite task can be stored in a file, called script-file, and then loaded into the MatLang environment to execute the task. The file name shall follow the convention of using 1.4. File Scripting and Function File 17 Figure 1.2: Surface representing an Hermite polynomial Figure 1.3: Contour plot of spatial Hermite polynomial 18 1.4. File Scripting and Function File Figure 1.4: Contour plot with lateral color bar scale and line labels the m extension like in interfer.m, the load of the file and its execution into the environment is done writing the name without the extension as in: > > interfer The script files, called m-files, are used to add more functionalities to the language by the definition of new function, in this case the name of the file and the name of the function stored into it should be the same. The structure of an m-file is described in the following example (the name of the m-file shall be function_name.m ): %% %% descriptive comment on the function task and arguments %% %% Copyright 1998 Roberto Bernetti %% %% This is free software; you can redistribute it and/or modify %% it under the terms of the GNU General Public License as published function [arg_out1 arg_out2 ... arg_outn]=function_name(arg_in1,arg_in2 ... arg_inn) ....... statements ....... end The starting comments will be used as help string, infact when the help command is issued in the following form 1.4. File Scripting and Function File 19 > > help function_name the output is the comment string set. The first excutable state- ment in the file shall be thefunctionstatement, its prototype is function [arg_out1 arg_out2 ..arg_outn]=function_name(arg_in1,arg_in2, ... arg_inn) where arg_outi are the output argument arg_ini are the input arguments It is possible to allow for a variable number of input arguments, in this case the alst or the only argument of the function shall bevarargin as in function [a b] = fun(x,y,varargin) then in function body you refer to the arrayvararginwhose number of elements arenargin. Here a simple example of a function creating the coefficient of the Newmark integration method: function [a,b,c]=tethaCoe(alpha,delta,tetha,dt) %%--------------------------------------------------------------- %% Purpose: %% Create coefficient for the Equilibrium collacation method %% theta=1 => Newmark method %% theta=1.4 => Wilson method %% %% Synopsis: %% [a,b,c]=tethaCoe(alpha,delta,tetha,dt) %% %% Variable Description: %% a - vector of "a" coefficent %% b - vector of "b" coefficent %% c - vector of "c" coefficent %% alpha - dispalacements difference coefficient %% delta - velocity difference coefficient %% tetha - Collocation parameter %% dt - time integration step %% %% Copyright 1998 Roberto Bernetti %% %% This is free software; you can redistribute it and/or modify 20 1.4. File Scripting and Function File %% %%--------------------------------------------------------------- a=[ 1/(alpha*(tetha*dt)^2) , ... 1/(alpha*tetha*dt) , ... (1/(2*alpha) -1) ]; b=[ delta/(alpha*tetha*dt) , ... (delta/alpha - 1) , ... (delta/alpha - 2)*tetha*dt/2 ]; c=[ a(1)/tetha , ... a(2)/tetha , ... (1/(2*alpha*tetha) - 1) ]; end Part II Beam Dynamics 21 Chapter 2 Dynamic Analysis of Beam In a rational approach to Mechanics the laws of dynamics can be usu- ally deduced from stationary conditions of some functionals, one of these is the Least Action Principle [Landau 1982]. From another point of view the Newton s Dynamic Laws incorporate the stationary conditions for appropriate functionals. This twofold approach is ap- plied in the Finite Element Method, where it is possible to get the discretized equations either from the stationarity condition of some functionals or directly from Newton s Dynamic Laws. In this chap- ter the stationary approach is used to deduce the Finite Element approximation of the dynamic behaviour of a beam disregarding its shear deformation. In the first section the kinematic relations, linking the displacements of the beam to the internal strain and stress, are described. In the second one the application of the Least Action prin- ciple leads to the approximated equation ruling the beam behaviour and in the last section some illustrative examples are presented. 2.1 Kinematics of Beam Two different approach has been generally used to describe the rela- tion between the displacements and the strain in a beam: " thin beam condition (Euler-Bernoulli theory [Calcote 1969]) " thick beam condition (Mindlin-Reissner theory [Auricchio 1994] ) 23 24 2.1. Kinematics of Beam z y x L Figure 2.1: Classical beam analysis reference frame The principal difference between the two theory is the different as- sumption about the shear deformation of the beam cross section. The thick beam model will be used in the next chapter presenting the vari- ational approach to the Finite Element Method. The reference system used in the following is showed in figure 2.1 For the sake of simplicity this section will focus on a 2-Dimensional example considering deflection of the beam only in the x - z plane. . 2.1.1 Euler-Bernoulli Theory The Euler-Bernoulli model of the beam behaviour is based on some assumptions that allow semplification in the description of the struc- ture deformation, and can be summarized as follow [Calcote 1969]: a) points laying on a normal to the beam center line will still lay on the normal to the deformed beam center line b) the elongation of a line drawn along the thickness, normal to the beam center line can be neglected c) line normal to the beam centerline before the deformation remain normal to the center line after the deformation 2.1. Kinematics of Beam 25 z ¸ Deflection line w x Figure 2.2: Physical quantities in the beam deflection From an energetic point of view hypothesis a), b) and c) state that only bending strain energy will be stored inside the beam as will be explained by the folllowing results. Under the action of a load normal to the center line of the beam, the structure deformation can be sketched as in figure 2.2. According to hypotesis b) the displacement in the z direction can be described by the displacement of the center line which is function of the x coordinate only: w = w(x) Following assumption a) at points with positive z coordinate the dis- placements in the x direction are negative and proportional to the rotation ¸y around the y axis of the normal (positive clockwise) and can be expressed by: u = z¸y w = w (2.1) According to hypothesis c) and considering the axis z = 0 as the iniital configuration the rotation results: "w ¸y = - (2.2) "x 26 2.1. Kinematics of Beam for the sake of brevity we will hereon use the following convention: "w wx = "x Inserting equation 2.2 into 2.1 the u displacement can be expressed as a function of w: u = -z · wx (2.3) Equation 2.2 or 2.3 is based on the assumption that shear stiffness of the cross section is infinite. This statement should be considered in a relative way respect to the bending stiffness. The correct assertion should be: the bending stiffness is several order of magnitude less then the shear stiffness. If, for any reason, the shear stifness decreases to a value comparable to the bending one the relationship represented by equation 2.2 is no more valid. To finally express the internal action and the deformation as func- tion of beam displacement now the stress shall be take into account. Considering small strain the only non-zero deformation tensor com- ponent is: µx = u (2.4) ,x Using the assumption b) and equation 2.3: µx = -z · wxx (2.5) inserting equation 2.5 into the Hook s law the stress inside the beam is: Ãx = E · z · ¸y,x (2.6) where: E is the Young modulus From the stress equation 2.6 it is possible to get the internal bend- ing moment: b h 2 2 Mx = E · z2 · ¸y,xdz b h - 2 2 in which: 2.1. Kinematics of Beam 27 z y x M dx Figure 2.3: Infinitesimal Beam deformation b width of the cross section beam h height of the cross section beam Recalling the definition of the moment of inertia of the cross section it is possible to set: b h 2 2 J = z2 · dz b h - 2 2 leading to the relation connecting cross section rotation to the bending moment: Mx = EJ · ¸y,x (2.7) According to the Euler assumptions inserting equation 2.2 into 2.7 results in the differential relation: Mx = -EJ · wxx (2.8) It is worth recall that the minus sign is due to the convention used for addressing the positive rotation in the right-handed reference system as depicted in figure 2.3. 28 2.2. Dynamics of Beam 2.2 Dynamics of Beam Kinematic relations is now used to state the energy principle then a stationary method is used for the deduction of the equations that rule the dynamic behaviour of an elastic beam. 2.2.1 Least Action Principle The principle of "Least Action states that the integral of the differ- ence of Kinetic and Potential energy between two instant t1 and t2 shall be minimum, in mathematical form:
t2 S(w) = (K(‡ (x, t)) - U(w(x, t))) dt (2.9) t1 where: w(x, t) displacements of the centroid of the beam cross section at x position at t time in the z direction K(‡ ) Kinetic energy of the beam U(w) Potential energy of the beam For the case under examination the expressions of the Kinetic and Potential Energy can be deduced considering the bending energy only according to the Euler-Bernoulli beam theory described in section 2.1.1 [Landau 1979]:
L 1 K(‡ ) = m · ‡ 2dx 2 0 2
L 1 "2w U(w) = E · J dx (2.10) 2 "x2 0 in which the reference system is showed in figure 2.1 and L is the beam length m is the mass per unit length of the beam 2.2. Dynamics of Beam 29 E is the Young modulus J is the second order momentum of the cross section dw ‡ = is the total time derivative dt The Action integral is therefore:
2 t2 L 1 1 "2w S(w) = m · ‡ 2 - E · J dx · dt (2.11) 2 2 "x2 t1 0 The configuration, described by the function w, will be the solution of the mechanical problem (for a correct statement of the problem we need the initial and boundary condition that we will deduce further on) if the Least Action Principle is satisfied, hence S(w) d" S(u) in which: u is in a Ä…-neighbor of w defined as follow: u = w + Ä… · ´w Ä… " R ´w satisfies the following conditions: ´w(t1) = 0 ´w(t2) = 0 (2.12) The functions w that describe the admissible deformation state of the beam are said admissible. A more rigorous definition of the set into which searching for the admissible function w and ´w will be discussed later. Now, the calculation of the first order differential of the functional in equation 2.11, as a function of Ä… define the variation of the functional S(w) [Smirnov 1980]: "S(u) ´S(Ä…) = · dÄ… "Ä… 30 2.2. Dynamics of Beam The stationary condition implies that ´S(0) = 0 hence: "S(u) = 0 (2.13) "Ä… Using the derivation under the integral sign the variation results:
t2 L ´S(0) = [m · ‡ ´‡ - E · J · wxx´wxx] dx · dt (2.14) t1 0 Using the Green theorem it is possible to rewrite the right side mem- ber of equation 2.14:
t2 L t2 L m‡ ´‡ · dx · dt = - m ´w · dx · dt t1 0 t1 0
t2 L EJ · wxx · dx · dt = EJ · wxx · ´wx|L - EJ · wxxx´w|L + 0 0 t1 0
t2 L + EJ · wxxxx · dx · dt (2.15) t1 0 In the first equation the condition expressed by relation 2.12 has been used. The variation of the Action functional is then modified as follow:
t2 L ´S = - [m + EJ · wxxxx] ´w · dx · dt + t1 0 + EJ · wxx · ´wx|L - EJ · wxxx´w|L (2.16) 0 0 The stationarity condition is satisfied if the right member is equal to zero for any choice of the variation ´w, which means each addend must be zero, leading to the following conditions: m + EJ · wxxxx = 0 EJ · wxx · ´wx|L = 0 0 EJ · wxxx´w|L = 0 (2.17) 0 The first of equations 2.17 is a Partial Differential Equation (PDE), the second and the third define the boundary conditions. In fact considering the second of equations 2.17 the left member might be zero 2.2. Dynamics of Beam 31 either if the first spatial derivative of the variation on the boundary is zero ´wx(0) = ´wx(L) = 0 or if the bending moment of the solution is zero: EJ · wxx(0) = EJ · wxx(L) = 0 the same consideration can be applied to the third of 2.17. According to the previous considerations the following classification is used: Essential Boundary Conditions boundary conditions that define the characteristic features of the set into which the solution will be searched (every function of the set satisfies these conditions) Natural Boundary Condition conditions that are naturally satis- fied by the function that verified the stationary condition hence makes the integral in equation 2.16 zero Now it is possible to define the mechanical problem of the dynamic oscillations of a beam fully restrained at the left end and free at the right one: define the set H2 of all integrable functions that satisfies the fol- lowing condition: &! = {x|0 d" x d" L}
t2 L [m · ‡ ´‡ - E · J · wxx´wxx] dx · dt = 0 "´w " H (2.19) t1 0 Remark 2.2.1 More restrictions have to be applied to the admissible function w in order the functional, in the form expressed by equa- tion 2.11, has a physical meaning. In fact admissible functions with 32 2.2. Dynamics of Beam first derivative discontinuos on a numerable subset &!n " &! hence w " C0(&!n), can still belong to the set H2(&!). From a physical point of view this feature means that dislocations are present inside the beam and the second addend in equation 2.11 can no more rep- resent the strain energy stored inside the beam without adding local energy contribution. For this reason, if no dislocations are present, C1(&!) continuity have to be required for the admissible functions. Moreover if we are searching solution in a set Hm and the ba- sis functions belong to Cm-1 simplier convergence condition can be deduced [Carey 1983]. According to equations 2.16 e 2.17 the statement in equation 2.19 is equivalent to the solution of the following PDE: m + EJ · wxxxx = 0 wx(0) = w(0) = 0 EJ · wxx(L) = EJ · wxxx(L) = 0 (2.20) Interesting consideration can be drawn if we assume that the time dependence is of harmonic type: w = w(x) · T (t) (2.21) this assumption is commonly used in linear free or forced vibration analysis and it is supported by the fact that under some non restric- tive conditions any function can be decomposed in an infinite sum of harmonic functions. According to this hypothesis the time derivative results:\ ´‡ = i · ´w w = i · w where: is considered the circular frequency Then equation 2.10 can be rewritten as follow:
t2 L t2 L -2 m · w´w · dx · dt - E · J · wxx´wxxdx · dt = 0 t1 0 t1 0 2.2. Dynamics of Beam 33 according to the definition of equation 2.10 the previous relation takes on the following form: - 2´K - ´U = 0 (2.22) leading to the following relation: ´U 2 = - (2.23) ´K Now considering the following functional: U(w) ²(w) = (2.24) K(w) the stationary condition results: ´U · K - U · ´K = 0 K2(w) equivalent to the numerator equal to zero the following relation apply: ´U U = = ²stat (2.25) ´K K Putting together equations 2.25 and 2.23 results: -2 = ²stat(w) Å» In other words the function w " H, which makes the quotient in Å» equation 2.24 stationary, verifies the stationary condition expressed by equation 2.22 too. In fact according to equation 2.25: ´U(w) = ²stat´K(w) Å» Å» hence: ²stat´K(w) - ´U(w) = 0 Å» Å» which shows that the function w satisfies the stationary condition of Å» the functional in equation 2.11. This results in a deduction from the Least Action Principle of the Rayleigh Principle [Krasnov 1984] also known in mathematics as the Courant min-max theorem [Courant 1943]. 34 2.2. Dynamics of Beam È i Ć Ć i-1 i i-1 i i+1 Figure 2.4: Shape function Ćj and basis function Èi 2.2.2 Finite Element The Finite Element Method (FEM) is a systematic way to find ap- proximation of the function w and its variation ´w. In this section the focus is on w(x) the stationary part of the function in equation 2.21. It is worth mentioning that there are finite element formulations taking into account the time dependent function too. The main steps of the approximation method are: 1. subdividing the x domain into n non-overlapping subdomain called exactly Finite Elements (FE) Ei = {x|xi-1 d" x d" xi} 2. each boundary point xi is called node 3. on each FE define the shape functions Ćj, one for each node that assumes unitary value on node j and zero on the others 4. define a set of n + 1 linearly independent functions Èi " H2 obtained patching togher the shape functions of the elements Ek with the node i in common; the functions Èi have finite support extended on all the elements attached to the corresponding node i. 2.2. Dynamics of Beam 35 The functions Èi form a finite basis for a subset Hn ‚" H2, in which the minimization of the functional in equation 2.26 will be searched. In fact any function in the subdomain Hn can be represented as: n
u = Ä…iÈi(x) i=1 in which Ä…i are unknown parameters. Hereafter the following conven- tions are used: n Ä…iÈi = Ä…iÈi repeated indexes are summed over the definition of i=1 the repeated index "2 Èi = Èi,xx the comma in the index represents differentiation "x2 According to remark 2.2.1 some restriction have to be applied to the functions Èi. In fact shape functions Ćj belonging to C1 over the element j can lead to basis functions Èi that have discontinuos first derivative on the inter-element boundary, resulting in: Èi " C0 For this reason, if no dislocations are present, C1 continuity have to be required for the basis function Èi at the node xi, in this case the element is said comforming element ensuring a stable and convergent behaviour decreasing the size of the FE. Some Finite Elements have been developed that do not respect requirment expressed by remark 2.2.1, mainly for 3D applications, they are defined non-conforming elements and require special attention to assure convergence, on the other hand their implementation is simpler respect to the conforming ones [Carey 1983]. The Action corresponding to the function u is calculated using equation 2.11:
t2 L 2 1 1 Sn(u) = m · Ä…iÈi · j - E · J (Ä…iÈi,xx)2 dx · dt 2 2 t1 0 (2.26) 36 2.2. Dynamics of Beam The function in subset Hn making S minimum must have the Ä…i parameters satisfying the following conditions:
t2 L 2 1 " 1 m · Ä…iÈi · j - E · J (Ä…iÈi,xxT )2 dx · dt = 0 "Ä…j t1 0 2 2 j = 1, 2, · · · n After some mathematics and using definition 2.2.1:
t2 L -2Ä…i m · ÈiÈjdx T · dt + t1 0
t2 L - Ä…i EJ · Èi,xxÈj,xxdx T · dt = 0 t1 0 j = 1, 2, · · · n (2.27) Due to the condition T (t) = 0 the integrand in equation 2.27 needs
to be zero in order to be satisfied, leading to the following:
L L 2Ä…i m · ÈiÈjdx + Ä…i EJ · Èi,xxÈj,xxdx = 0 0 0 j = 1, 2, · · · n (2.28) Using the following position:
L Mji = m · ÈiÈjdx mass matrix 0
L Kji = EJ · Èi,xxÈj,xxdx stiffness matrix 0 it is possible to write down the well known eigenvalue equation in vector notation: K · u + 2M · u = 0 · K (2.29) where: K = [Kij] the stiffness matrix M = [Mij] the mass matrix u = [Ä…i] the vector of the unknown parameters The free vibration of beam is approximated by the eigenvalue prob- lem represented by equation 2.29, frequency and eigenfunction in the continuos model are substituited by eigenvalues and eigenvectors. 2.2. Dynamics of Beam 37 z Ć Ć i-1 3 i 1 Ć i i-1 4 x Ć i 2 Figure 2.5: Basis functions relevant to node i 2.2.3 Hermite Polynomials Now it is possible to build the shape function and the resulting basis functions according to the requirments of remark 2.2.1 and definition 2.18. The simplest way is using third order Hermite polynomial, hence for the generic i element the shape function results [Paz 1985]: x2 x3 Ći 1 = 1 - 3 + 2 2 3 li li 2 x Ći 2 = -x 1 - li x2 x3 Ći 3 = 3 - 2 2 3 li li
x2 x Ći 4 = - - 1 (2.30) li li The unknown parameters now assume a clear meaning, they are the nodal displacements and rotations at the nodes. Patching together the shape functions of two adiacent Finite Elements we forms the basis function (see figure 2.4) 38 2.3. Numerical Exercise 2.3 Numerical Exercise Some example are now presented to show the main features of the FEM applied to the dynamics of an Euler-Bernoulli beam. The fol- lowing data are used: L = 1 beam length b = 0.1 cross section width h = 0.1 cross section height E = 1 Young modulus Á = 1 mass density Two dfferent boundary conditions have been used, both ends free and left end fixed and right one free. The script that has been used is exebeam.m it has the following stages: 1. calculating the local stiffness and mass matrix 2. assembling of the local matrix into the global matrix 3. applying homogeneus boundary condition 4. calulating eigenvalues and associated eigenvector 5. calculating internal action The start section sets the geometrical and mechanical data of the beam and the dimension of the global stiffness and mass matrix, then a loop starts on each element to calculate the element stiffness and mass matrix using the function BeamLocM, its prototype is the follow- ing: [K,M,ML]=BeamLocM(L,E,J,rho,A) For the mass matrix two output are available: Mloc the consistent mas matrix calculated using equation 2.28 2.3. Numerical Exercise 39 MlocL the lumped mass matrix with only the diagonal elements, corresponding to the translational DOF, different from zero Then the local matrixes are assembled by the function assmat its prototype is: [kglo]=assmat(kglo,kloc,index) the assembling process is ruled by the vector index which contains at its j position the values of the global DOF corresponding to the j - th local DOF For the Fixed-Free ends case the boundary conditions at the left end have to be applied by eliminating the rows and columns corre- sponding to the DOFs at the restrained node. In this case this has been accomplished sending to the eig function only the submatrix of active DOFs as in the command [V2L,D2L]=eig(Kglo(3:DOF,3:DOF),MgloL(3:DOF,3:DOF)); Obtained the eigenvectors the natural frequency and mode shapes are calculated according to the resulta of [Paz 1985] in the function TheFreq which prototype is: [F,u,M]=TheFreq(ibou,L,E,J,rho,A,x) The next step is the calculation of the internal action for the FE model using function IntAction, its prototype is: IntAction(u,L,E,J,rho,A,n) The internal action is calucalted using the matrix product between the stifness local matrix and the local DOFs solution vector. In this way the moment at the left end is positive if clockwise (x - z plane figure 2.2) while the internal action is positive according to the curvature if counter-clockwise. For this reason there is a sign change in the moment calculation. 40 2.3. Numerical Exercise 1st-mode Exact FEM 4 Ele. FEM 16 Ele. Consistent Mass 0.1027920 0.1029037 0.1027925 Lumped Mass 0.1027920 0.0868712 0.1015660 Table 2.1: Exact and approximated frequency for Free-Free ends beam 1st-mode Exact FEM 4 Ele. FEM 16 Ele. Consistent Mass 0.01615393 0.0161545 0.01615400 Lumped Mass 0.01615453 0.0868712 0.01612510 Table 2.2: Exact and approximated frequency for Fixed-Free ends beam 2.3.1 Natural Frequncy The results of the calulation are reported in table 2.1 for the Free- Free ends and in table 2.2 for the Fixed-Free ends. The decreasing of the element size is the better way to verify the model correctness and accuracy. It is possible in fact to verify the stationary condition of the natural frequncy as stated by equation 2.23, with the following consideration: " using the consistent mass approach, which means the Kinetic energy is calculated in a way consistent with the displacement field used for the strain energy, the exact natural frequency is a minimum for the admissible function set, reaching with few element an acceptable level of accuracy " using the lumped mass matrix calculation the natural frequancy is a maximum for the Rayleigh quotientover the admissible func- tion set. It is left to reader modify the script to try with different boundary conditions ad test different natural frequency orders. 2.3.2 Mode Shapes Ineresting consideration can be drawn from the analysis of the shape mode. First it is worth noting that the current examples allow to 2.3. Numerical Exercise 41 analize the meaning of the Natual Boundary Condition. In fact as defined in section 2.2.1 the Natural boundary conditions are not ap- plied locally but they are satisfied by the convergence process. As showed by figure 2.6 the moment is different from zero at the free ends of the beam. As the convergence increase using more elements ("x 0) the values of the bending moment at the free ends converge to zero as depicted in figure 2.10. So the boundary condition is natu- rally satisfied by the numericla model without applying any restraint to the trial function used. Another consideration steam directly from the definition of the subset Hn, the moment is discontinuos at the FE nodes as clearly showed by figures 2.6 and 2.8. The basis functions are continuos till the first derivative. This lack of continuity is balanced by the convergence process which ensure for the decrease in the jump at the node as described by the figure 2.10 and 2.11. The amoount of this jump is often used to assess the level of accuracy of the model used in real life applications. Different behaviour is showed by the use of lumped mass matrix, in this case the bending moment is continuos through the node too. This difference is showed by the results of the same calculation done using the commercial Finite Element code ANSY Stm in figure 2.12. 42 2.3. Numerical Exercise (a) (b) Figure 2.6: Mode shape for displacements (a) and Moment (b) for first mode Free-Free ends (4 elements) 2.3. Numerical Exercise 43 (a) (b) Figure 2.7: Mode shape for displacements (a) and Moment (b) for first mode Free-Free ends (4 elements) with Lumped Mass Matrix 44 2.3. Numerical Exercise (a) (b) Figure 2.8: Mode shape for displacements and Moment for first mode Free-Free ends (4 elements) 2.3. Numerical Exercise 45 (a) (b) Figure 2.9: Mode shape for displacements (a) and Moment (b) for first mode Free-Free ends (4 elements) with Lumped Mass Matrix 46 2.3. Numerical Exercise (a) (b) Figure 2.10: Mode shape for displacements (a) and Moment (b) for first mode Free-Free ends (16 elements) 2.3. Numerical Exercise 47 (a) (b) Figure 2.11: Mode shape for displacements (a) and Moment (b) for first mode Free-Free ends (16 elements) 48 2.3. Numerical Exercise (a) (b) Figure 2.12: Bending moment for a Free-Free ends with consistent (a) and lumped (b) mass matrix using ANSYS Chapter 3 Mixed Finite Elements In chapter 2 the Finite Element Method has been introduced us- ing the stationary approach and using the beam displacements and roatations as unknown parameters (refer is done to section 2.2.2). In present chapter the variational approach is presented together with the inclusion of shear effect in the cross sectin deformation. Moreover in the numerical approximation presented in chapter 2 the primary variables, displacements and rotations, are more accurate then the derived one, bending moment and shear, because the last are usu- ally obtained by derivative of the primary ones. From an engineering point of view the stress inside the material is of main interest for this reason the family of Mixed Finite Elements has been developed to increase the accuracy of the more interesting variables. This goal is obtained using as primary variable the stresses instead of rotation re- sulting in a lower degrees of derivative in the numerical scheme. Here the Mixed Finite Elements formulation is presented applying it ot the beam dynamics assuming as primary variable the displacements and the bending moment, instead of rotations. 3.1 Governing Equation In the variational approach the starting point is the Partial Differ- ential Equations which will be satisfaied using a weak formulation. Considering the beam structure in the mixed formulation the mo- 49 50 3.1. Governing Equation z V x M* q(x) Figure 3.1: Force diagram sketch ment is assumed as primary variable as the displacement. The final governing equation should be slightly modified to take this change into account, and in the following this formulation is presented. First the equilibrium equation of a rectilinear beam will be formulated then the compatibility equation connecting deformation and internal action is added to the system. 3.1.1 Dynamic Equilibrium The application of the second Newton law of dynamics to the infineti- samal beam stretch depicted in figure 3.1 leads to the equilibrium differential equation. The rotational equilibrium applied to the beam stretch dx of figure 3.1 reads: "Mx -Vxzdx - Mx + Mx + dx = 0 "x after simplification: "Mx Vxz = (3.1) "x 3.1. Governing Equation 51 and the translation in the z direction: "Vxz "2w -Vxz + Vxz + dx + q(x) = Ac "x "t2 after simplification: "Vxz "2w + q(x) = Ac (3.2) "x "t2 in which: w is the transverse beam displacement in the z direction Vxz is the beam shear q(x) is the applied load per unit length Ac = ÁA is the mass per unit length of the beam Combining equations 3.1 and 3.2 it is possible to obtain: "2Mx "2w + q(x) = Ac (3.3) "x2 "t2 Assuming no shear deformation and the ¸y expression presented in section 2.2 results: "2w ¸y,x = - "x2 in which coma stand for derivation, and considering the relationship connecting bending moment and beam cross section rotation as de- duced in section 2.1 at equation 2.8 the final system of governing equations results the following one:
"2w Mx + E = 0 "x2 (3.4) "2M + q(x) - Ac "2w = 0 "x2 "t2 The boundary condition shall be added to system 3.4 to uniquely determine the solution. Usually two different sets of conditions are defined: 52 3.1. Governing Equation " essential boundary conditions that are applied to the un- known function in the form w = bw(t) or M = bm(t) at one or both the beam ends, where bw(t)and bm(t) are known function; " natural boundary conditions that are applied on the share force and on the rotation in the form Vxz = bv(t)or ¸y,x = b¸(t) at one or both the beam ends, where bv(t) and b¸(t) are known function; 3.1.2 Thick Beam The needs to consider more general deformation conditions lead to the Timoshenko-Mindlin-Reissner theory [Auricchio 1994] which drops the condition c) of section 2.1.1, but holding a) and b) ones, assuming that points laying on a straight line before the deformation remain on a straight line. According to the previous assumption the points laying along the normal to the reference surface before the deformation lay no more on the normal to the deformed surface. For this reason the relations, expressed by equations 2.2, have to be written in a more general form. According to the sign convention represented in figure 3.2 it is possible to define: ¸y rotation of a line normal to the reference surface around the y axis (positive clockwise in the sheet plane) Å‚xz angular shear deformation in the x - z plane (positive clockwise) w x component of the normal to the deformed beam centre ,x line (positive counterclockwise) and drawn the following relationships: ¸y = Å‚xz - w (3.5) ,x Using the rotation ¸y it is possible to express the displacements using the equation 2.1: u = z¸y w = w (3.6) 3.1. Governing Equation 53 z n x Å‚xz Deformed line x ¸y (a) didascalia tetay Figure 3.2: Shear and bending deformation relation, arrows point to positive rotation If the shear effect has to be taken into account the first equation of 3.4 should be modified considering the following hypothesis: " the shear deformation is assumed to be uniform in the cross section of the beam and is called Å‚xz As a consequence the relation between the shear stress Äxz and the internal action Vxz assumes the following form: h h 2 2 Vxz = b Äxzdz = µ · b GÅ‚xzdz h h 2 2 in which: b is the cross section base µ is the cross section shape factor taking into account for the differ- ence between real shear deformation variation along the thick- ness and the constant assumption G is the shear modulus in the x - z plane of the beam 54 3.2. Weighted Residual Using the assumption of constant shear deformation the final relation is: Vxz = µG · A · Å‚xz (3.7) where: A is the cross sectional area Inserting equation 3.1 into 3.7 reads: "M = µG · A · Å‚xz (3.8) "x Now deriving respect to x the shear deformation definition: Å‚xz,x = ¸y,x + w (3.9) ,xx performing the substitution of equation 2.7 and 3.8 it is possible to write: 1 "2M M = + w ,xx µG · A "x2 EJ leading to the system of governing equation: Å„Å‚ Mx 1 "2Mx "2w - + = 0 òÅ‚ EJ µG "x2 "x2 (3.10) ół "2Mx + q(x) - Ac "2w = 0 "x2 "t2 in which the shear deformation effect has been taken into account. The boundary conditions are applied in the same manner as in section 3.1.1. 3.2 Weighted Residual The numerical solution of the system of equations 3.4 can be ob- tained via the usual Finite Element technique where the transverse displacement function w and the bending moment Mx are assumed as unknowns in both formulation with and without the shear effect. For a coherent approximation theory it is necessary to state a little bit more formally the previous problem as follow: 3.2. Weighted Residual 55 find Mx(x), w(x) " C2 on the interval [0, L] satisfying the bound- ary conditions and the following differential equation: Å„Å‚ Mx 1 "2M "2w - + = 0 ôÅ‚ EJ µGA "x2 "x2 ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ "2M ôÅ‚ + q(x) - Ac "2w = 0 ôÅ‚ "x2 "t2 òÅ‚ (3.11) ôÅ‚ w(0, t) = 0 ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ ¸y(0, t) = b¸(0, t) ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ Mx(L, t) = 0 ôÅ‚ ół Vxz(L, t) = bv(t) The previous statement defines the problem with the classical formu- lation. As consequence it is not always possible to find a solution for the problem due to the restriction imposed to the degree of continuity of the solution. 3.2.1 Weak Formulation Finite Element Method, instead of a classical formulation of the prob- lem, starts usually from less restrictive set of possible solution to find then in this new set the approximated solution, this is usually referred as weak formulation . The unknown function can be approximated in a broader set then the classical one C2. This approximated set, called H1, is the set of all functions with at least the first derivative square integrable. The next step is to define a subset H1 ‚" H1 con- h structed using the base function Ni(x). In the following development the subscript i connects the base function to the node of the elements used to discretize the domain by the following: Ni(x) is different from zero only on the elements that are attached to or contain the node i The unknown function can then be expressed by: Mh(x) = MiNi(x) wh(x) = wiNi(x) i = 1, 2, ..n (3.12) in which the repeating indeces summation convention has been used and: Mh approximated function of the bending moment 56 3.2. Weighted Residual wh approximated transverse displacement Mi value of the bending moment at node i wi value of the transverse displacement at node i n total number of nodes The above relations if substituted in system 3.10 give rise to the fol- lowing residues: Å„Å‚ Mh 1 "2Mh "2wh ôÅ‚ - + = R1(x, t) òÅ‚ EJ µGA "x2 "x2 ôÅ‚ ół "2Mh + q(x) - Ac "2wh = R2(x, t) "x2 "t2 Now the Finite Element Method looks for a solution of the problem in the subset H1 weighting the previous residues with trial functions h equal to the Nj(x). This approach was introduced by Galerkin and named Galerkin method after him. The aim is to represent the physical problem described by the system of equations 3.11 using a different class of solution functions and it is formulated as follow: 1 find Mh(x), wh(x) " Hh on the interval [0, L] satisfying the bound- h ary conditions wh(0, t) = 0, ¸y (0, t) = 0 and the following system of 1 equations for each vh " Hh: Å„Å‚
L L Mh 1 "2Mh vhdx - vhdx+ ôÅ‚ ôÅ‚ 0 EJ 0 µGA "x2 ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ h ôÅ‚ +wi 0L "2wh vhdx + (¸y - b¸(t))vh = 0 ôÅ‚ ôÅ‚ "x2 0 òÅ‚ (3.13) ôÅ‚ ôÅ‚ ôÅ‚ L L ôÅ‚ "2Mh ôÅ‚ vhdx + q(x)vhdx+ ôÅ‚ 0 "x2 0 ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ L
ół h
-Ac 0L "2wh vhdx - (Vxz - bv(t))vh = 0 "t2 in which: (a) the repeated indexes summation convention is used 3.2. Weighted Residual 57 (b) the boundary condition, expressed through the functions b¸(t) and b¸(t), for the shear force (last term in the left hand member of the second equation) and for the rotation (last term in the left hand member of the first equation), are applied in a weak form called Natural Boundary Con- ditions (c) the boundary conditions for the displacement and the mo- ment in a strong form, Essential Boundary Conditions. As clearly stated by system 3.13 the Natural Boundary Conditions will be enforced by the solution of the linear system of equations whilst as it results from definition 3.13 Essential Boundary Conditions will be enforced by the appropriately choice of the nodal values. The relation expressed by equation 3.13 is satisfied if it is true for each base function, leading to: Å„Å‚
Mi L Mi L "2Ni NiNjdx - Njdx+ ôÅ‚ ôÅ‚ EJ 0 µGA 0 "x2 ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚
L Mi 0L "2Ni Njdx + q(x)Njdx+ "x2 0 ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ L
ôÅ‚ ôÅ‚ h ôÅ‚ -Ac d2wi L NiNjdx - (Vxz - bv(t))Nj = 0 ôÅ‚ dt2 0 ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ ół j = 1, ..n where use of the equation 3.12 has been done. Applying integration by parts it is possible to highlight the presence of the boundary con- ditions: Å„Å‚ L
Mi L Mi L "Ni "Nj Mi "Ni ôÅ‚ NiNjdx + dx - Nj 0 + ôÅ‚ EJ 0 µGA 0 "x "x µGA "x ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ L h
ôÅ‚ L "Ni "Nj ôÅ‚ ôÅ‚ -w dx + wi "Ni Nj 0 - (¸y - b¸(t))Nj 0 = 0 ôÅ‚ 0 "x "x "x òÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ L L
ôÅ‚ ôÅ‚ ôÅ‚ -Mi 0L "Ni "Nj dx + Mi "Ni Nj 0 + q(x)Njdx+ ôÅ‚ "x "x "x 0 ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ L ôÅ‚ ół h -Ac d2wi L NiNjdx - (Vxz - bv(t))Nj = 0 dt2 0 58 3.2. Weighted Residual Referring to equation 3.9 and 3.1 it is possible to write: Mi "Ni "Ni h ¸y = - wi µG "x "x "Ni h Vxz = Mi "x Except for the boundary node the base function The approximated system of equations finally is: Å„Å‚ Mi L Mi L "Ni "Nj NiNjdx + dx - wi 0L "Ni "Nj dx ôÅ‚ ôÅ‚ EI 0 µG 0 "x "x "x "x ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ L ôÅ‚ ôÅ‚ Mi h i ôÅ‚
ôÅ‚ ôÅ‚ ôÅ‚ -Mi 0L "Ni "Nj dx + q(x)Njdx + -Ac d2wi L NiNjdx ôÅ‚ "x "x dt2 0 ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ L h L ół + Mi "Ni Nj 0 - (Vxz - bv(t))Nj = 0 "x and Å„Å‚ Mi L Mi L "Ni "Nj NiNjdx + dx - wi 0L "Ni "Nj dx+ ôÅ‚ ôÅ‚ EJ 0 µGA 0 "x "x "x "x ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ L ôÅ‚ h h h ôÅ‚ + ¸y Nj 0 - (¸y - b¸(t))Nj 0 - ¸y Nj = 0 ôÅ‚ ôÅ‚ òÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ L ôÅ‚ ôÅ‚ -Mi 0L "Ni "Nj dx + q(x)Njdx - Ac d2wi L NiNjdx+ ôÅ‚ "x "x 0 dt2 0 ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ L h L h ół h + VxzNj - (Vxz - bv(t))Nj - VxzNj 0 = 0 Considering now the boundary conditions the shape function is zero, for the displacement at x = 0, with Nj(0, t) = 0 " j and for the moment at x = L with Nj(L, t) = 0 " j, canceling out the relevant 3.2. Weighted Residual 59 terms. Finally the system of equations reads: Å„Å‚
Mi L Mi L "Ni "Nj ôÅ‚ NiNjdx + dx - wi 0L "Ni "Nj dx+ ôÅ‚ EJ 0 µGA 0 "x "x "x "x ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ + b¸(t)Nj|0 = 0 ôÅ‚ ôÅ‚ òÅ‚ (3.14) ôÅ‚ ôÅ‚
ôÅ‚ L ôÅ‚ ôÅ‚ -Mi 0L "Ni "Nj dx + q(x)Njdx - Ac d2wi L NiNjdx+ ôÅ‚ "x "x 0 dt2 0 ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ ół + bv(t)Nj|L = 0 The previous equation give rise to a system of 2n linear equations the whom solution determines the nodal values of the displacement and moment. 3.2.2 Shape Functions It is useful to write in matrix form equation 3.14 : Kb · m + Ks · m - Kn · w = -¸ (3.15) d2 -Kn · m - Mi w = -Q - V (3.16) dt2 in which:
L 1 (Kb)ij = NiNjdx Bending compliance matrix EJ 0
L Qj = q(x)Njdx Distribute load Vector 0 mi = Mi Vector of the bending moment nodal values wi = wi Vector of the transverse displacements nodal values ¸ Vector of the boundary condition applied to the end beam rotations 60 3.2. Weighted Residual V Vector of the boundary condition applied to the end beam shear force The matrix elements can be obtained once the function Ni of the Finite Elements has been chosen. Here it is presented, as an exam- ple, the case of polynomial function of 2th degrees with the following expressions:
L x - (x - L) 2 N1(x) = L2 2 x (x - L) L N2(x) = L - L 2 2
L x x - 2
N3(x) = L L L - 2 in which: l is the length of the finite element The first derivative of the shape functions are:
L 2 x - 2 (x - L) d 2 N1(x) = + dx L2 L2 d 4 (x - L) 4x N2(x) = - - dx L2 L2
L 2 x - 2x d 2 N3(x) = + dx L2 L2 Now it is straitforward to obtain the matrix: ëÅ‚ öÅ‚ 2L L L - 15 15 30 1 L 8L L íÅ‚ Å‚Å‚ Kb = 15 15 15 EI L L 2L - 30 15 15 ëÅ‚ öÅ‚ 7 8 1 - 3L 3L 3L 1 8 16 8 íÅ‚ Å‚Å‚ Ks = - - 3L 3L 3L µG 1 8 7 - 3L 3L 3L ëÅ‚ öÅ‚ 7 8 1 - 3L 3L 3L 8 16 8 íÅ‚ Å‚Å‚ Kn = - - 3L 3L 3L 1 8 7 - 3L 3L 3L 3.2. Weighted Residual 61 ëÅ‚ öÅ‚ 2L L L - 15 15 30 L Å‚Å‚ Mi = Ac íÅ‚ 15 8L L 15 15 L L 2L - ëÅ‚ öÅ‚30 15 15 L 6 2L íÅ‚ Å‚Å‚ Q = 3 L 6 The two vector ¸ and V have all null element but the first and last ones associated to the boundary values. The matrix equations 3.15 and 3.16 can be combined to use only one matrix equation by defining the following coefficient matrixes:
Kb + Ks -Kn K = -Kn 0
0 0 M = 0 Mi
-¸ L = -Q - V and for the unknown vector:
m u = w leading to the following usual matrix equation for the discrete system: d2 K · u - M u = L (3.17) dt2 The form of equation 3.17 is very useful for static problems where the equation reduces to: K · u = L and the boundary conditions can be applied in the usual form with the same subroutines used for any others Finite Element program. In the dynamic case a different formulation has to be used due to the singularity of the mass matrix M. 62 3.2. Weighted Residual 3.2.3 Displacement Only Method The systems of equations 3.15 and 3.16 can be solved one at a time, deducing the moment from the first system as a function of the dis- placement and then substituting in the second system, as hereafter described (with no lack in generality an for the sake of simplicity no shear effect is considered): m = Kb-1 · (Kn · w - ¸) d2 -Kn · Kb-1 · (Kn · w - ¸) - Mi w = -Q - V (3.18) dt2 The solution of system of equations 3.18 requires that the matrix Kb is non singular as it is assured by the linearly independent shape functions. The second matrix equations of 3.18 now has only dis- placement degrees of freedom and the matrix Mi is non singular too for the same reason stated above. Usual methods for structural time integration can now be used considering the following stiffness matrix: and for The governing system of equations can now be written in the following form: d2 - K · w - Mi · w = L (3.19) dt2 in which: K = Kn · Kb-1 · Kn the stiffness matrix L = -Q - V - (Kn · Kb-1 · ¸) the load vector Some attention shall be placed in considering the boundary condi- tions, that allows to make the stiffness matrix non-singular. This topic will be discussed in the following section. 3.2.3.1 End Rotations In this case the boundary conditions are expressed imposing the ro- tations at both ends in the form: ¸y(0) = b¸1(t) ¸y(L) = b¸2(t) 3.2. Weighted Residual 63 This means that the moments at the same ends are unknown. The ap- plication of system 3.18 is straightforward: the vector ¸ is completely known, all elements zero but the first and the last ones: îÅ‚ Å‚Å‚ b¸1(t) ïÅ‚ śł 0 ïÅ‚ śł ¸ = (3.20) ïÅ‚ śł . . ðÅ‚ ûÅ‚ . -b¸2(t) The nodal moment values are expressed by: m = Kb-1 · (Kn · w - ¸) (3.21) Substituting in the second of the equation 3.18 the final form is ob- tained: d2 - K · w - Mi w = -Q - V - (Kn · Kb-1 · ¸) (3.22) dt2 For the application of displacement boundary conditions the proce- dure is outlined hereafter in section3.2.3.3. 3.2.3.2 Applying End Moments More attention has to be payed in the case of end moments constrain that can be expressed as: Mx(0) = bm1(t) Mx(L) = bm2(t) (3.23) Now the end rotations are unknown and the boundary condition, expressed by equations 3.23, have to be added to system 3.26 to solve the problem. Different solution techniques have to be used depending on the type of analysis, they will be described in the next sections. Static Analysis In this case equations 3.23 have to be added to system 3.19 to obtain the solution in term of displacements and end rotations. The two boundary conditions can be introduced into the 64 3.2. Weighted Residual first and last equations of the first matrix equation of system3.18, they read: (Kb-1 · Kn)1j · wj - Kb-1¸1 = bm1(t) j = 1, 2..n 11 (Kb-1 · Kn)nj · wj - Kb-1¸n = -bm2(t) j = 1, 2..n(3.24) nn The two equations in 3.24 can now be added to the system 3.22 as the initial and final rows along with the two additional columns, deduced from the matrix equation 3.18: (Kn · Kb-1)j1 · ¸1 j = 1, 2..n (Kn · Kb-1)jn · ¸n j = 1, 2..n The resulting stiffness matrix reads: îÅ‚ Å‚Å‚ -Kb-1 (Kb-1 · Kn)1i -Kb-1 11 1n ðÅ‚ K = (Kn · Kb-1)j1 -(Kn · Kb-1 · Kn)ji (Kn · Kb-1)jn ûÅ‚ -Kb-1 (Kb-1 · Kn)ni -Kb-1 n1 nn (3.25) The system can now be written, referring to equation 3.18 (the minus sign of the stiffness matrix is put inside the stiffness matrix definition) considering the static conditions: îÅ‚ Å‚Å‚ îÅ‚ Å‚Å‚ ¸1 0 ðÅ‚ ûÅ‚ ðÅ‚ ûÅ‚ K · w = L (3.26) ¸n 0 in which: L = -Q - V is the load vector 3.2.3.3 End Displacements/Forces The solution of both matrix equations 3.22 and 3.26 require the ap- plication of displacements boundary condition that shall be applied following the standard Finite Element procedures: 1. if the displacement is imposed the relevant system row is taken out from the solution procedure and the unknown values of the shear force will be determined after the time integration 3.3. Full Time Integration 65 2. if the shear force is imposed the known values is set on the right hand side of equation 3.22 When the displacements are calculated the equation 3.21 allows for the determination of the nodal values of the moment. 3.3 Full Time Integration In the case of dynamic analysis the general discrete equation 3.19 can be used with some slight modifications to the stiffness and mass ma- trix. The stiffness matrix is exactly the same determined in equation 3.25, for the mass matrix the contribution of the end rotations is nil so the modified matrix is easily obtained: îÅ‚ Å‚Å‚ 0 0 0 ðÅ‚ ûÅ‚ Mw = 0 M i 0 0 0 0 The final equation reads: îÅ‚ Å‚Å‚ îÅ‚ Å‚Å‚ îÅ‚ Å‚Å‚ ¸1 0 d2 ¸1 ðÅ‚ ûÅ‚ ûÅ‚ ðÅ‚ ûÅ‚ - K · w - Mw ðÅ‚ w = -Q(t) - V (t) (3.27) dt2 ¸n ¸n 0 The numerical solution of equation 3.27 will be calculated using the standard time integration techniques for structural dynamics, the Wilson-¸method is the one used in this case. Here only the main results for the development of the method will be addressed and for a deeply discussion the reader is reminded to [Hilber 1976]. The method is based on the assumption that the acceleration varies linearly inside a time step. The time step is comprised be- tween the two temporal instants tk and tk + ¸"t. The acceleration is expressed by the following equations: ć k+¸ = (1 - ¸)ć k + ¸Ä‡ k+1 (3.28) in which: ć second time derivative of the unknown vector of equation 3.27 66 3.3. Full Time Integration tk+1 = tk + "t discrete time Integration of acceleration in the interval [tk, tk +¸"t] gives the veloc- ity and the displacement as polynomial expression of ¸"t [Hilber 1976, Paz 1985]: ‹ k+¸ = ‹ k + ¸"t[(1 - ´)ć k + ´Ä‡ k+¸] (3.29) 1 xk+¸ = xk + ¸"tvk + ¸2"t2[( - Ä…)ć k + ąć k+¸] (3.30) 2 Substitution of the expression for the kinematics quantities in matrix equation 3.27 reads: K · xk+¸ + C‹ k+¸ + Mw · ć k+¸ = F k+¸ in which: ëÅ‚ öÅ‚ 0 íÅ‚ Å‚Å‚ F = Qk+¸ + V k+¸ k+¸ 0 Expressing discrete acceleration and velocity as a function of discrete displacement through equations 3.28, 3.29 and 3.30 leads, after re- arraying of parameters, to the final expression: [a0Mw + b0C + K] · xk+¸ = F + k+¸ +Mw [a0xk + a1‹ k + a2ć k] + +C [b0xk + b1‹ k + b2ć k] (3.31) in which: 1 a0 = Ä…¸2"t2 1 a1 = Ä…¸"t 1 a2 = - 1 2Ä… ´ b0 = Ä…¸"t ´ b1 = - 1 Ä… 3.3. Full Time Integration 67 ´ b2 = (2Ä… - 1)¸"t The equation 3.31 can be used to determine the unknown vector at time instant tk+¸ once the whole kinematic states are known at time instant tk. Then using equations 3.28 it is possible to calculate the acceleration at time tk+1 = tk + "t: ak+1 = c0(xk+¸ - xk) - c1‹ k - c2ć k in which: a0 c0 = ¸ a1 c1 = ¸
1 c2 = - 1 2Ä…¸ Then using equations 3.29 and 3.30 with ¸ = 1 it is possible to cal- culate the displacements and the velocities at time tk+1 = tk + "t. If ¸ = 1 the algorithm is the classical Newmark. The following worth characteristics of the algorithm can be listed : Stability for unconditional stability, in choosing the time step "t, the following inequalities must hold [Hilber 1976]: 1 ´ e" 2 ´ Ä… e" 2 2Ä… ¸ e" 1 - 2Ä… Accuracy for every values of ¸ the order of accuracy is two if [Hilber 1976]: 1 ´ = 2 Dissipation No numerical damping for [Hilber 1976]: ´ = 0.5 Ä… = 0.25 68 3.4. Numerical Exercise and the better values for optimal dissipation characteris- tics, for both low and high frequency is [Hilber 1976]: (´ + 0.5)2 Ä… = 4 Damping The numerical solution results to be under-damped if the following inequalities hold: 1 ´ e" 2 (´ + 0.5)2 Ä… e" 4 The better behaviour in terms of frequency induced dissi- pation for unconditionally stable algorithm is for (damp- ing increasing with frequency): ¸ H" 1.4 Example of application of the integration algorithm will be showed in the next sections. 3.4 Numerical Exercise A simple model, for which analytical solution is available, has been chosen to test the performance of the Mixed Finite Element formu- lation. On that model, eigen mode extraction and time integration analysis have been run comparing the result of the Mixed Finite Ele- ments (heron called MFE), the Classical Finite Elements (heron called CFE) and if available analytical solution. The script used for the test are eig-solve.m for mode/frequency extraction and dyn_solve.m for the full time integration. 3.4.1 Eigen Modes For the sake of simplicity, without losing in generality, a uniform rect- angular cross section beam with free end translations and restrained 3.4. Numerical Exercise 69 length 10 m N Young Modulus 2.1E+11 m2 Kg density 7250 m3 width 0.1m height .25m Table 3.1: Geometrical and Mechanical data Mode nth Exact MFE ANSYS 1 0.0 0.2247e-05 0.2008E-06 2 6.1011 6.1042 6.0992 3 24.4045 24.588 24.449 4 54.9101 56.739 55.677 Table 3.2: Eigen Frequency results end rotations is considered. For the MFE equation 3.22 has been solved, for the CFE the solution has been obtained using the general purpose Finite Element Program ANSYS and reference for the ana- lytical solution can be found in [Paz 1985]. The following geometrical and mechanical data has been used: The beam has been discretized using a Finite Element size equal to 2.5m resulting in 5 node, a 2nd degree polynomial shape function for the MFE, and the usual 3rd degree polynomial shape function for the CFE. Eigen frequency results, for the various methods, are reported in table 3.2. They show a good agreement between both numerical methods and the exact solution considering that only 9 degree of freedom are present for the MFE and 10 for the CFE. A maximum relative error, on the third mode, of 3.3% for the MFE formulation and 1.37% for the CFE. The higher value for the MFE can be easily explained with the lower degree of the approximating polynomial for the displacement unknowns used respect to the CFE. This kind of behaviour is completely reversed considering the stresses that is in this case the bending moment value. The approximating polynomial is of order two in the MFE formulation and the shape func- tion is continuous, instead for the CFE the approximating polynomial is of order one and the shape function is discontinuous crossing the 70 3.4. Numerical Exercise 4e+06 4e+06 Mixed FE Mixed FE Exact Sol. Exact Sol. Classic FE Classic FE 3e+06 3e+06 2e+06 2e+06 1e+06 1e+06 0 0 -1e+06 -1e+06 -2e+06 -2e+06 -3e+06 -3e+06 -4e+06 -4e+06 0 2 4 6 8 10 0 2 4 6 8 10 Figure 3.3: Bending moment for the 1st mode. Space coordinate (m) on abscissa and Moment on ordinate (Nm) element border. As it is clearly showed by figures 3.3, 3.4 and 3.5 the MFE formulation gives higher level of approximation both in terms of continuity and precision. In the CFE formulation the discontinuity on the nodes and the accuracy of the values decrease quickly with the increasing of the mode number. In the light of the previous results it is possible to draw the conclusion that with a lower degree of approx- imation in the displacements there is no sensible lack in the frequency accuracy calculation while it is possible to gain a strong increase in the accuracy of the bending moment calculation. This increase is much more impressive if using the same degree of approximation of the CFE for the shape functions, ie third degree polynomials. Numerical and exact solutions are very close in this case as figure 3.6. This behaviour can still be recognized in the bending moment shape mode for the 4th mode as described in figure 3.7. 3.4.2 Time Integration For the test of the time integration performances of the MFE method the case of the structure presented in section 3.4.1, with simply sup- 3.4. Numerical Exercise 71 2e+07 Mixed FE Exact Sol. Classic FE 1.5e+07 1e+07 5e+06 0 -5e+06 -1e+07 -1.5e+07 -2e+07 0 2 4 6 8 10 Figure 3.4: Bending moment for the 2nd mode. Space coordinate (m) on abscissa and Moment on ordinate (Nm) Mixed FE Exact Sol. Classic FE 4e+07 2e+07 0 -2e+07 -4e+07 0 2 4 6 8 10 Figure 3.5: Bending moment for the 3rd mode. Space coordinate (m) on abscissa and Moment on ordinate (Nm) 72 3.4. Numerical Exercise Mixed FE Exact Sol. 4e+07 2e+07 0 -2e+07 -4e+07 0 2 4 6 8 10 Figure 3.6: Bending moment for the 3rd mode with four node Mixed Element. Space coordinate (m) on abscissa and Moment on ordinate (Nm) Mixed FE Exact Sol. 4e+07 2e+07 0 -2e+07 -4e+07 0 2 4 6 8 10 Figure 3.7: Bending moment for the 4th mode with four node Mixed Element. Space coordinate (m) on abscissa and Moment on ordinate (Nm) 3.4. Numerical Exercise 73 ported ends, one constant point load suddenly applied at midspan, at time t = 0, has been chosen with the same discretization used for the modal analysis. For this particular load configuration an analytical solution [Paz 1985] is available for the beam deflection in the following form:
" 2P0 1 nÄ„ nÄ„ z(x, t) = sin (1 - cosÉnt)sin x AcL É2 2 L n=1 and for the bending moment:
" 2Ä„2P0 n2 nÄ„ nÄ„ m(x, t) = EI sin (1 - cosÉnt)sin x AcL3 É2 2 L n=1 in which Én is the nth eigen-value. The numerical results (obtained with the octave script dyn_solve.m revision 1.13 with iload set to 1) has been obtained, both for the MFE and CFE method, using the classical Newmark method, the only available in the ANSYS software. Neither numerical nor physical damping has been used in both model. Results are showed in terms of displacement and bending moment at the loaded point in figure. The bending moment of the first node of the last element is showed in figure. The time history of the numerical and exact solution are very close, almost identical for the displacement and with some error for the bending moment. For a better understanding of the accuracy of the numerical solutions and for a comparison between the MFE and the CFE the RMS of the FFT of the time history, that describe the energy content of the time series signal, has been calculated and reported in table 3.3. The CFE data are on an element basis so for each node the right and left element data have been considered. Remembering that the numerical solution is always stiffer that the exact one as it is expected the exact solution has the lower energy level. The MFE numerical solution shows a close RMS value to the exact one respect to the CFE method, calculated with the ANSYS software. 74 3.4. Numerical Exercise 5e-08 MFE Analytical 0 -5e-08 -1e-07 -1.5e-07 -2e-07 -2.5e-07 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Figure 3.8: Displacement of the loaded point. Time on the abscissa (sec.) and displacement on the ordinate (m) MFE Analytical 0 -0.05 -0.1 -0.15 -0.2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Figure 3.9: Bending moment at the loaded point. Time on the ab- scissa (sec.) and moment on the ordinate (Nm) 3.4. Numerical Exercise 75 MFE Analytical 0 -0.05 -0.1 -0.15 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Figure 3.10: Bending moment at the first node of the last element. Time on the abscissa (sec.) and moment on the ordinate (Nm) node Exact MFE CFE right CFE left 5 2.9352 2.9566 NA 2.9652 4 2.3512 2.3519 2.3614 2.3600 3 1.6534 1.6528 1.6622 1.6602 2 0.86640 0.86668 0.87646 0.87268 1 0.000000 0.000000 0.027169 NA Table 3.3: RMS comparison of the Fast Fourier Transform for the bending moment time series data 76 3.4. Numerical Exercise Chapter 4 Spatial Frame Applications Application of the Displacement Only Method in spatial frame re- quires that rotation of element matrix and axial stress component should be accounted for. 4.1 The column behaviour The effect of the axial force has to be considered in spatial beam structures, in this the equilibrium equation is (reference is made to figure 4.1): "2u " Ac dx = P - P "t2 in which: u is the axial displacement P is the axial force " "P P = P + dx "x hereafter for other symbols reference is made to section 2.1.1 and 3.1.1. After some simplification and with the compatibility equation the system that rules the rod behaviour is: "u P - EA = 0 "x 77 78 4.1. The column behaviour z y P* P x dx Figure 4.1: Beam axial equilibrium sketch "P "2u - + Ac = 0 (4.1) "x "t2 Using the same weak formulation and finite element approximation of the section 3.2.1 the equation 4.1 is transformed as:
L L "Ni NjNidx · Pi - EA Nj dx · ui = 0 "x 0 0
L L "Ni d2ui h - Nj dx · Pi + Ac NjNidx · + (P - bp(t))Nj 0 = 0 "x dt2 0 in which: Pi value of the axial force at node i ui value of the axial displacement at node i Using integration by parts in the equilibrium equation, the final form of the system is:
L L "Ni NjNidx · Pi - EA Nj dx · ui = 0 "x 0 0
L L "Nj d2ui + · Nidx · Pi + Ac NjNidx · - bp(t)Nj|L = 0 0 "x dt2 0 0 With the following position: 4.1. The column behaviour 79
L (Ki)ji = NjNidx 0
L "Nj (Ku)ij = (Kp)ji = · Nidx 0 "x K · p - EA · Ku · u = 0 Kp · p + AcKi · ü = 0 in which : p is the vector of the axial force nodal value u is the vector of the axial displacement nodal values Now it is possible to condensate the axial force in the matrix form: p = EA · Ki-1 · Ku · u Kp · Ki-1 · Ku · u + AcKi · ü = 0 (4.2) where the stiffness matrix is: Kax = Kp · Ki-1 · Ku Considering the same shape function of section 3.2.2 the matrix rela- tionship for a three node element is: ëÅ‚ öÅ‚ ëÅ‚ öÅ‚ ëÅ‚ öÅ‚ ëÅ‚ öÅ‚ 2 1 -1 P1 -1 2 -1 u1 2 2 3 6 L 2 íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚ 1 8 1 P2 - EA -2 0 u2 = 0 3 3 15 1 -1 1 2 P3 -2 1 u3 2 ëÅ‚ öÅ‚ ëÅ‚ öÅ‚ ëÅ‚6 3 2 öÅ‚ ëÅ‚ öÅ‚ -1 2 -1 P1 2 1 -1 ü1 2 3 6 2 1 2 íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚ -2 0 P2 + AcL 1 8 1 ü2 = 0 3 3 15 1 -2 1 P3 -1 1 2 ü3 6 3 2 2 Eliminating the axial force from the equilibrium equation it is possible to get: ëÅ‚ öÅ‚ ëÅ‚ öÅ‚ 7 -8 1 u1 3 3 3 EA íÅ‚ Å‚Å‚ íÅ‚ -8 16 -8 u2 Å‚Å‚ + 3 3 3 L 1 -8 7 u3 3 3 3 ëÅ‚ öÅ‚ ëÅ‚ öÅ‚ 2 1 -1 ü1 2 1 íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚ +AcL 1 8 1 ü2 = 0 15 -1 1 2 ü3 2 With similar calculations it is possible to deduce the stiffness matrix for shape function of different polynomial degrees. 80 4.2. Matrix Rotation w 3 u 3 M 3 w 2 M 2 u 2 ² w 1 u 1 z G M 1 x G Figure 4.2: DoF sketch for 2D spatial frame 4.2 Matrix Rotation In 2D geometry it is necessary to consider for each node two transla- tional Degrees of Freedom, one in the x direction and the other in the z direction, as showed in the figure 4.2, where a general configuration for an element is sketched. The stiffness matrix in the global refer- ence system (see figure 4.2) has to be obtained by its rotation from the local reference system (see figure 4.3). The local stiffness matrix is obtained coupling the column behaviour, described in section 4.1, and the beam behaviour described in section 3.2.1 represented by the equation 3.14. For the sake of comprehension of the following analysis it is useful to rewrite equation 3.14 with both end rotations prescribed 4.2. Matrix Rotation 81 w w w 2 1 3 u u 1 3 u 2 M M M 2 1 3 zL xL Figure 4.3: Local reference system for the Element (the first equation results multiplied for EJ respect to equation 3.14): Å„Å‚
L MiEJ "Ni "Nj ôÅ‚ Mi 0L NiNjdx + dx - EJ · wi 0L "Ni "Nj dx+ ôÅ‚ µGA 0 "x "x "x "x ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ òÅ‚ - EJ · b¸(t)Nj|L + EJ · b¸(t)Nj|0 = 0 ôÅ‚ L ôÅ‚ ôÅ‚ -Mi 0L "Ni "Nj dx + q(x)Njdx - Ac d2wi L NiNjdx+ ôÅ‚ "x "x 0 dt2 0 ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚ ół + bv(t)Nj|L = 0 Along with the stiffness matrix the load vector has to undergone the same process of transformation too, taking into account that the load vector has a contribution from the end rotations. This effect is quan- tified by (referring to equations 3.14, 3.15 and 3.20): q¸ = -(Kn · Kb-1 · ¸) It has to be considered that the vector ¸ has only two elements differ- ent from zero, the end rotations (again reference is made to equation 3.20) and that the inversion of the matrix Kb has to be referred to the structure assembled matrix. 4.2.1 Spatial Behaviour To express the stiffness matrix in an form suitable for the rotation procedure it is worthwhile to re-arraying together both equations, 4.2 82 4.2. Matrix Rotation and 3.18, in the compatibility equation:
Ki 0 p · + 0 Kb m L
-EA · Ku 0 u 0 = (4.3) 0 -EI · Kn w EI · ¸ L L and in the equilibrium equation:
Kp 0 p · + 0 -Kn m L
d2 Mi 0 u 0 · = 0 Mi w -Q - V dt2 L L Obtaining p and m from the first equation: -1 p Ki 0 = · m 0 Kb L
EA · Ku 0 u 0 · · + 0 EI · Kn w -EI · ¸ L L it is possible to express the equilibrium equation in the displacement nodal value unknowns: -1 Kp 0 Ki 0 · · 0 Kn 0 Kb
EA · Ku 0 u · · + 0 EI · Kn w L
d2 Mi 0 u 0 + · = + 0 Mi w Q + V dt2 L L
Kp 0 + · 0 Kn -1 Ki 0 0 · 0 Kb EI · ¸ L 4.2. Matrix Rotation 83 Re-arraying the nodal displacement vector as follow (for an n nodes elements): îÅ‚ Å‚Å‚ u1 ïÅ‚ śł w1 ïÅ‚ śł ïÅ‚ śł . . UL = ïÅ‚ śł . ïÅ‚ śł ðÅ‚ ûÅ‚ un wn The equilibrium equation results: d2 KL · UL + ML · UL = Q + KR · G · Q¸ (4.4) dt2 in which (for the sake of simplicity a three node element has been chosen): îÅ‚ Å‚Å‚ 0 0
ïÅ‚ śł EI 0 ¸1 ïÅ‚ śł Q¸ = G = ïÅ‚ śł . . . . ¸n ðÅ‚ ûÅ‚ . . 0 -EI (4.5) îÅ‚ Å‚Å‚ (Kax)11 0 (Kax)12 0 (Kax)13 0 ïÅ‚ 0 (Kfl)11 0 (Kfl)12 0 (Kfl)13 śł ïÅ‚ śł ïÅ‚ śł (Kax)21 0 (Kax)22 0 (Kax)23 0 ïÅ‚ śł KL = ïÅ‚ 0 (Kfl)21 0 (Kfl)22 0 (Kfl)23 śł ïÅ‚ śł ðÅ‚ ûÅ‚ (Kax)12 0 (Kax)12 0 (Kax)33 0 0 (Kfl)31 0 (Kfl)32 0 (Kfl)33 The same arrangement has to be used for the mass matrix MLand for the node rotation coefficient matrix. In the analysis of frame structures, respect to the analysis developed for supported beams in section 3.2.3.1, the minus sign of the left node rotation and EI have to be moved to the rotation coefficient matrix, KR, due to the 84 4.2. Matrix Rotation assembling process of the rotation equilibrium equation of the frame joint. The rotation coefficient matrix, KR, results: îÅ‚ (Kp · Ki-1)11 0 (Kp · Ki-1)12 ïÅ‚ 0 (Kn · Kb-1)11 0 ïÅ‚ ïÅ‚ (Kp · Ki-1)21 0 (Kp · Ki-1)22 ïÅ‚ KR = ... ïÅ‚ 0 (Kn · Kb-1)21 0 ïÅ‚ ðÅ‚ (Kp · Ki-1)12 0 (Kp · Ki-1)12 0 (Kn · Kb-1)31 0 Å‚Å‚ 0 (Kp · Ki-1)13 0 (Kn · Kb-1)12 0 (Kn · Kb-1)13 śł śł śł 0 (Kp · Ki-1)23 0 śł (Kn · Kb-1)22 0 (Kn · Kb-1)23 śł śł ûÅ‚ 0 (Kp · Ki-1)33 0 (Kn · Kb-1)32 0 (Kn · Kb-1)33 Ki and Kb must result from the assembling process of the elements belonging to the stretch under analysis. In the previous matrix defi- nitions the following position has been made: Kax = Kp · (EA · Ki-1) · Ku Kfl = Kn · (EI · Kb-1)Kn Equation 4.4 is the building block for the assembling process, neces- sary for the solution of 2-D frame structures. 4.2.2 Local to Global System Considering the transformation between the local and global reference systems, the relationship connecting the local displacements compo- nents of the j node to global ones reads: ULj = R · UGj in which the rotation matrix is: îÅ‚ Å‚Å‚ cos xLxG cos xLĆyG cos xLĆzG Ć ðÅ‚ xG yG zG ûÅ‚ R = cos yLĆ cos yLĆ cos yLĆ cos zLĆ cos zLwG cos zLÄ™ G xG 4.2. Matrix Rotation 85 Applying the transformation to the vector of displacements of the finite element gets (considering for simplicity a three node element): îÅ‚ Å‚Å‚ îÅ‚ Å‚Å‚ îÅ‚ Å‚Å‚ UL1 R 0 0 uG1 ðÅ‚ ûÅ‚ ðÅ‚ ûÅ‚ ðÅ‚ UL2 = 0 R 0 UG2 ûÅ‚ UL3 0 0 R UG3 In a more concise form: UL = T · UG (4.6) in which: T is the transformation matrix Thus for a plane-structural problem the transformation matrix is: îÅ‚ Å‚Å‚ cos² sin² 0 0 0 0 0 0 0 ïÅ‚ śł -sin² cos² 0 0 0 0 0 0 0 ïÅ‚ śł ïÅ‚ śł 0 0 1 0 0 0 0 0 0 ïÅ‚ śł ïÅ‚ śł 0 0 0 cos² sin² 0 0 0 0 ïÅ‚ śł ïÅ‚ śł T = 0 0 0 -sin² cos² 0 0 0 0 ïÅ‚ śł ïÅ‚ śł 0 0 0 0 0 1 0 0 0 ïÅ‚ śł ïÅ‚ śł 0 0 0 0 0 0 cos² sin² 0 ïÅ‚ śł ðÅ‚ ûÅ‚ 0 0 0 0 0 0 -sin² cos² 0 0 0 0 0 0 0 0 0 1 If a 2D element is considered as in figure 4.2:
uLj ULj local vector of node j displacements wLj
uGj UGj global vector of node j displacements wGj Paying attention to the 2D case and considering ² the angle between the xL and the xG axis the transformation equation becomes:
cos² sin² ULj = UGj -sin² cos² 86 4.2. Matrix Rotation Using the orthogonality feature of the rotation matrix results:
cos² -sin² UGj = ULj sin² cos² Applying the transformation from global to local of the vector of displacements for the whole finite element gets: îÅ‚ Å‚Å‚ îÅ‚ Å‚Å‚ îÅ‚ Å‚Å‚ uL1 cos² sin² 0 0 0 0 uG1 ïÅ‚ śł ïÅ‚ śł ïÅ‚ śł wL1 -sin² cos² 0 0 0 0 wG1 ïÅ‚ śł ïÅ‚ śł ïÅ‚ śł ïÅ‚ śł ïÅ‚ śł ïÅ‚ śł uL2 0 0 cos² sin² 0 0 uG2 ïÅ‚ śł ïÅ‚ śł ïÅ‚ śł = ïÅ‚ śł ïÅ‚ śł ïÅ‚ śł wL2 0 0 -sin² cos² 0 0 wG2 ïÅ‚ śł ïÅ‚ śł ïÅ‚ śł ðÅ‚ ûÅ‚ ðÅ‚ ûÅ‚ ðÅ‚ uL3 0 0 0 0 cos² sin² uG3 ûÅ‚ wL3 0 0 0 0 -sin² cos² wG3 here îÅ‚ Å‚Å‚ cos² sin² 0 0 0 0 ïÅ‚ śł -sin² cos² 0 0 0 0 ïÅ‚ śł ïÅ‚ śł 0 0 cos² sin² 0 0 ïÅ‚ śł T = (4.7) ïÅ‚ śł 0 0 -sin² cos² 0 0 ïÅ‚ śł ðÅ‚ ûÅ‚ 0 0 0 0 cos² sin² 0 0 0 0 -sin² cos² Then applying the transformation relationship to the stiffness matrix results: T KG = T · KL · T (4.8) The same transformation of equation 4.8 is applied to the mass ma- trix: T MiG = T · MiL · T (4.9) The frame rotation, mathematically described in equations 4.8 and 4.9, has to be applied to the element characteristic quantities before the assembling process will be performed. To this aim bring back equation 4.4 and applying to it the rules of equation 4.6 it reads: d2 KL · T · UG + ML · T · UG = T · QG + KR · G · Q¸ dt2 4.2. Matrix Rotation 87 T then multiplying both sides for T : d2 T T · KL · T · UG + T · ML · T · UG = T · T · QG + dt2 T +T · KR · G · Q¸ T and finally using the definition of equation 4.8 and the relation T · T = I: d2 T KG · UG + MG · UG = QG + T · KR · G · Q¸ (4.10) dt2 Here in a 2D condition the end rotation vector undergoes an identity transformation during reference frame rotation. The equation 4.10 allows for using a standard element stiffness, mass matrix and load vector assembling process only for the end ro- tations contribution some attention has to be payed. 4.2.2.1 End Rotations Considering a linear structure with constant cross section, the matrix T product T ·KR is identical for two adjacent elements of the structure. Moreover the vector of end rotations has the following form (where the adjacent elements are labeled i and j with the common node l): îÅ‚ Å‚Å‚ îÅ‚ Å‚Å‚ 0 0 0 0
ïÅ‚ śł ïÅ‚ śł EI 0 EI 0 ¸h ¸l ïÅ‚ śł ïÅ‚ śł G·Q¸j = G·Q¸i = ïÅ‚ śł ïÅ‚ śł . . . . . . . . ¸l ¸m ðÅ‚ ûÅ‚ ðÅ‚ ûÅ‚ . . . . 0 -EI 0 -EI (4.11) as a consequence the assembling process will cancel out the contri- bution of the element common end rotations but the case of element with different cross section and T joint configurations. In the latter T cases the matrix product T ·KR is different for two adjacent elements although the common node rotation is identical, as a consequence the joint (or discontinuity node) rotations will still appear in the right hand member of equation 4.10, unknown at the solution stage. For this reason some other conditions have to be introduced to account for the node rotations. 88 4.2. Matrix Rotation The rotational equilibrium of the joint is the condition that makes the system 4.10 solvable. To this aim, firstly it is necessary to account for the global compo- nents of the node displacements, in a manner similar to the one used for equation 4.4, introducing the following position: îÅ‚ Å‚Å‚ p1 ïÅ‚ śł 1 ïÅ‚ śł ïÅ‚ śł . . s = ïÅ‚ śł . ïÅ‚ śł ðÅ‚ ûÅ‚ pn n to re-array the coefficient matrix of equation 4.3 : KS · sL + KD · UL = -G · Q¸ (4.12) Considering a three node element the matrix KS and KD are as follow: îÅ‚ Å‚Å‚ (Ki)11 0 (Ki)12 0 (Ki)13 0 ïÅ‚ 0 (Kb)11 0 (Kb)12 0 (Kb)13 śł ïÅ‚ śł ïÅ‚ śł (Ki)21 0 (Ki)22 0 (Ki)23 0 ïÅ‚ śł KS = ïÅ‚ 0 (Kb)21 0 (Kb)22 0 (Kb)23 śł ïÅ‚ śł ðÅ‚ ûÅ‚ (Ki)12 0 (Ki)12 0 (Ki)33 0 0 (Kb)31 0 (Kb)32 0 (Kb)33 îÅ‚ EA(Ku)11 0 EA(Ku)12 ïÅ‚ 0 EI · (Kn)11 0 ïÅ‚ ïÅ‚ EA · (Ku)21 0 EA(Ku)22 ïÅ‚ KD = ... ïÅ‚ 0 EI · (Kn)21 0 ïÅ‚ ðÅ‚ EA · (Ku)12 0 EA(Ku)12 0 EI · (Kn)31 0 Å‚Å‚ 0 EA · (Ku)13 0 EI · (Kn)12 0 EI · (Kn)13 śł śł śł 0 EA · (Ku)23 0 śł EI · (Kn)22 0 EI · (Kn)23 śł śł ûÅ‚ 0 EA · (Ku)33 0 EI · (Kn)32 0 EI · (Kn)33 4.2. Matrix Rotation 89 Jm Jn ¸Jm ¸Jn ¸l ¸m ¸k Jl Jk Figure 4.4: Example sketch of a frame structure It is important to highlight that the previous matrix KS and KD must be assembled for each node belonging to the stretch under anal- ysis to get the complete form 4.12. Equation 4.12 expresses the in- ternal action, Mx and P , as a function of the deformation (U and ¸) of the element. Considering figure 4.4 the bending moment exerted on joint Jn is a function of the vertical displacements and rotation of the orizontal stretch between Jn and Jm. To calculate the moment on joint Jn the deformation of each element on the orizontal stretch have to be accounted for. Then assembling matrixes KS and KD and recalling equation 4.11 only joint rotations, ¸J and ¸J , give rise to a n m contribution different from zero. The resulting matrix equation reads (here the matrix KS and KD resutlt from the assembling performed on the element of the stretch from Jn to Jm): sL = -K-1 · (KD · T · UG + G · Q¸) (4.13) S in which the global compoments of the displacements have been used. From relation4.13 the row correspong to the moments acting on the joints of the structure has to be cut off (for the sake of semplicity the subscript L , referring to the local component, has been dropped and substituted by the vector component number and by the label of the joint which the moment acts on and the stretch which belongs to) : (sJ )s = -(KS)-1 · (KD)jmTml · (UG)l + m 2j -(KS)-1 · (Q¸)1 + (KS)-1 · (Q¸)2 22 2k (sJ )s = -(KS)-1 · (KD)jmTml · (UG)l + m kj 90 4.3. Numerical Exercise -(KS)-1 · (Q¸)1 + (KS)-1 · (Q¸)2 k2 kk j, m, l = 1, 2, ......k (4.14) in which k is the number of the nodes in the stretch s . On the right end of the element the internal positive bending moment laied in the negative side of global reference system, for this reason the second equation have to be multiplyed by -1: (sJ )s = -(KS)-1 · (KD)jmTml · (UG)l + m 2j -(KS)-1 · (Q¸)1 + (KS)-1 · (Q¸)2 22 2k (sJ )s = +(KS)-1 · (KD)jmTml · (UG)l + m kj +(KS)-1 · (Q¸)1 - (KS)-1 · (Q¸)2 k2 kk Getting equation 4.14 for each stretch that converge to the joint it is possible to write down the rotational equilibrium condition (referring is made to the example figure 4.4): (sJ )1 + (sJ )2 = 0 (4.15) n n Adding, for each joint , the equilibrium equation 4.15 to the sys- tem of equations 4.10 the solution of the structural problem can be obtained using the usual methods for inversion of the system of al- geabric equations. When the displacements and joint rotations of the structure has been obtained the moment should be calculated using the first equation of the system 4.14. 4.3 Numerical Exercise The MBEM (Mixed Beam Element Method) has been tested versus the standard Finite Element Methods using its implementation in the general purpose package ANSYS. A simple frame geometry has been used with unitary dimension as represented in figure 4.5, with a sud- denly applied unitary force at time t = 0. The frame characteristics are listed in table 4.1. The element size used was "h = 0.3333 with a total of 10 node for both model and 3 element for the MBEM and 9 for the ANSYS model. The calculation for the MBEM is done by the script dyn3Dframe.m. 4.3. Numerical Exercise 91 F 1 1 Figure 4.5: Sketch of the Frame geometry E 4.6 · 109 Á 1.6 · 103 b 0.1 h .7 · 10-2 Table 4.1: Mechanical and geometric characteristic of the cross section of the frame 92 4.3. Numerical Exercise The results are reported, in figure 4.6, in terms of horizontal dis- placements and bending moments of the upper left node. Both figure show how the MBEM gets better results above all in term of bending moment that for the ANSYS model is discontinuos at the node, two values have to be reported. The dynamic displacements oscillates about the value results from the statical application of the load. Such value is 0.019015 the average value of the dynamic analysis is 0.018998 with a percent error of 0.09% for the MBEM while for the ANSYS model the average value is 0.018635 with a percent error of 1.998%. This kind of result is more stressed for the bending moment which, in the ANSYS model, is discontinuos while in the MBEM is continuos and with an error of 0.773% on the average value. 4.3. Numerical Exercise 93 0.045 Mixed Beam Element ANSYS Solution 0.04 0.035 0.03 0.025 0.02 0.015 0.01 0.005 0 0 0.5 1 1.5 2 2.5 3 (a) Displacement 1.2 Mixed Beam Element ANSYS Solution ANSYS Solution 1 0.8 0.6 0.4 0.2 0 0 0.5 1 1.5 2 2.5 3 (b) Bending Moment Figure 4.6: Time history of the ANSYS and Mixed-FE model solu- tions at the upper-left node 94 4.3. Numerical Exercise Bibliography [Courant 1943] R. Courant, "Variational methods for the solution of problems of equilibrium and vibrations," Bull. Amer. Math. Soc., 49 (1943), pp. 1-23. [Courant 67] R. Courant, K. Friedrichs, and H. Lewy, "On the partial difference equations of mathematical physics," IBM J. Res. Develop., 11 (1967), pp. 215-234. [Calcote 1969] Lee. R. Calcote, "The Analysis of Laminated Com- posite Plate" Van Nostrand Reinhold Company 1969 [Hilber 1976] Hans M. Hilber, Analysis and design of numerical inte- gration method in structural dynamics , College of engi- neering University of California, Report No. EErc 76-29 November 1976 [Landau 1979] L. D. Landau, E. M. Lifsits, Fisica teorica/Teoria dell elasticitÄ… , Editori Riuniti, Roma 1979 [Smirnov 1980] V. I. Smirnov, Corso di Matematica Superiore , Vol. IV parte prima, Editori Riuniti 1980 [Landau 1982] L. D. Landau, E. M. Lifsits, Fisica teor- ica/Meccanica , Editori Riuniti, Roma 1982 [Carey 1983] G. F. Carey, J. T. Oden, Finite Elements: A second course, Prentice-Hall 1983 [Krasnov 1984] M.L. Krasnov, G. I. Makarenko, A.I. Kiselev, Calcolo delle variazioni , Edizioni MIR Mosca 1984 95 96 BIBLIOGRAPHY [Paz 1985] M. Paz, Structural Dynamics Theory & Computation , Van Nostrand Reinhold New York, 1985 [Jonson 1987] C. Jonson, Numerical solution of differential equations by the finite element method , Studentlitteratur, Lund 1987 [Auricchio 1994] F. Auricchio, R.L. Taylor, A shear deformable plate element with an exact thin limit , Computer Methods in Applied Mechanics and Engineering, vol 118,pg 393- 412, 1994 [Kwon 1997] Y. W. Kwon, H. Bang, The Finite Element Method using Matlab , CRC Press 1997 Index admissible function, 29, 31 integrable functions, 30 arithmetic operator, 7 internal action, 37 inverse matrix, 8 bending moment, 26 bending stiffness, 26 Kinetic energy, 28 center line, 24 language, 3 comforming element, 34 laws of dynamics, 23 condition number, 8 Least Action, 27 Least Action Principle, 23, 28, 32 decomposition, 8 logical opeators, 11 deformation tensor, 26 determinant, 7 MATLAB, 3 Matrix, 4 eigenvalue equation, 35 Mindlin-Reissner theory, 23 eigenvalues, 8 moment of inertia, 27 eigenvectors, 8 Essential Boundary conditions, 30 Natural Boundary Conditions, 30 Euler-Bernoulli theory, 23 non-conforming elements, 34 normal, 25 finite basis, 33 numerical, 4 Finite Element Method, 23 FreeMat, 3 Octave, 3 function, 4 Potential energy, 27 functionals, 23 program flow, 11 GPL, 3 range operator, 6 harmonic functions, 31 Rayleigh Principle, 32 Hermite polynomial, 36 reference system, 24 homogeneus boundary, 37 Hook s law, 26 shape functions, 33 97 98 INDEX shear deformation, 23 shear stiffness, 26 small strain, 26 stationary conditions, 23 strain energy, 31 thick beam, 23 thin beam, 23 transpose, 7 unknown parameters, 34, 36 variation, 33