Machine Translated by Google

Democratic and People's Republic of Algeria

Ministry of Higher Education and Scientific Research.
Abderrahmane Mira University – Bejaia
Faculty of Technology
Department of Civil Engineering

Course
Finite Element Method

Prepared and presented by

Abdelghani SEGHIR
Doctor of Science from A. Mira University, Béjaia, Algeria Doctor of
Civil Engineering from Paris-Est University, Marne-la-Vallée, France

2005-2014
Machine Translated by Google

Contents
chapter 1 Presentation of MATLAB .................................................. ............................. 5

1.1 Introduction................................................. .................................................. ...... 5

1.2 Introduction of matrices........................................................ ................................... 7

1.3 Operations on matrices .............................................. ................................... 8

1.4 Operations on vectors.............................................. .................................. 10
1.5 Functions on matrices.............................................. ................................... 11

1.6 Access to matrix elements............................................ .......................... 12

1.7 Loops in MATLAB............................................ ................................. 13

1.8 Logic tests ............................................. ............................................. 14

1.9 Function files............................................. ............................................. 15

1.10 Application 1............................................ .................................................. .... 16

1.11 Exercise................................................ .................................................. ........... 17

1.12 Duty - TP........................................................ .................................................. .......... 17

chapter 2 General presentation of the Finite Element Method............................................ 1

2.1 Introduction................................................. .................................................. ...... 1

2.2 History of the method .......................................... ............................................. 2

2.3 The main points of the method........................................................ ............................. 3

2.4 Variational formulation ................................................ .................................. 5

2.4.1 Strong form................................................ .................................................. 5

2.4.2 Weak form ................................................ ................................................. 7

2.5 Discretization of the domain........................................................ ................................... 8

2.6 Approximation on the element............................................ ................................... 9

2.6.1 Polynomial approximation and nodal approximation................................ 9

chapter 3 Differential equations and boundary problems........................................................ 11

3.1 First order differential equation ......................................... ............ 11
3.2 Steps to follow ............................................... .................................................. .11
Machine Translated by Google

3.2.1 Formulation of the problem............................................. ........................... 11

3.2.2 Discretization of the domain........................................................ .......................... 12

3.2.3 Discretization and interpolation on the element ............................................. .. 13

3.2.4 Properties of shape functions ............................................. ................ 14
3.2.5 Lagrange type element......................................................... ........................... 15
3.2.6 Elementary matrices................................................ ................................. 16

3.2.7 Noticed ................................................. .................................................. 19

3.3 Step 4: Assembly ............................................. ........................................... 19

3.4 Step 5a: Resolution – Application of CALs ......................................... .......... 23

3.5 Step 5b: Resolution – Calculation of the solution ......................................... ........... 24

3.6 Study of convergence............................................ ....................................... 24

3.7 Assignment No. 2 ............................................. .................................................. ......... 27

3.8 2nd order differential equations ............................................ ................... 27

3.9 MATLAB program........................................................ ....................................... 29

3.10 Comparison with step-by-step integration........................................ ................ 31

3.11 Duty No. 3 .................................................. .................................................. .... 33

3.12 General program for 2nd order equations ......................................... 35

3.13 Duty No. 4 .................................................. .................................................. .... 45

chapter 4 Bar Element............................................ .................................................. ... 47

4.1 Governing equation................................................ ............................................. 47
4.2 Wording of the element ............................................. ................................... 48

4.3 Example of a recessed rod........................................................ ................................. 50

4.4 Plane lattice structures .............................................. ................................... 51

4.5 Example of two bars............................................ ....................................... 56

4.6 Assembly techniques for bar elements............................................ .58

4.7 Finite element program for the bar element............................................ ..... 59

4.7.1 Data entry............................................... ............................................. 59

4.7.2 Calculation of stiffness and mass matrices............................................ ..... 61

3
Machine Translated by Google

4.7.3 Application of support conditions............................................. ............. 64

4.7.4 Discrete system resolution............................................. .................... 66
4.7.5 Main program................................................ ................................. 67
4.8 TP No. 3 ............................................. .................................................. ............... 70

chapter 5 Beam Element .......................................... .................................................. .71

5.1 Governing equation................................................ ............................................. 71
5.2 Wording of the element ............................................. ................................... 75

5.2.1 Variational formulation ................................................ ......................... 75

5.2.2 Discretization.............................................. ............................................. 76

5.2.3 Elementary matrices................................................ ................................. 78

5.3 Example of a console beam ......................................... .............................. 81

5.3.1 Case of concentrated loading at the free end............................................ 81
5.3.2 Case of distributed loading ......................................... ......................... 82

5.3.3 MATLAB program .................................................. ......................... 85

5.3.4 SAP model .................................................. .............................................. 86

5.4 Duty No. 6 .................................................. .................................................. .... 88

4
Machine Translated by Google

chapter 1

MATLAB Overview

1.1 Introduction
MATLAB is an interactive software based on matrix calculation (MATrix LABoratory).
It is used in scientific calculations and engineering problems because it allows
complex numerical problems to be solved in less time required by common
programming languages, thanks to a multitude of built-in functions and several tool
programs. tested and grouped according to use in files called toolboxes.

When starting MATLAB on a PC, the interface looks like this:

Figure 1.2: MATALB main window (version 7.0)

Machine Translated by Google

MEF course – Chapter 1: Presentation of MATLAB A. Seghir 2005-2014

The main MATLAB window contains two secondary windows that can be moved or
closed.
On the right, the commands window allows you to enter commands line by line and
display the results. The >> symbol indicates waiting for a command.
On the left, the Workspace, Current Directory and sometimes Command History windows are nested in
panes .

- Workspace allows you to display the variables used with their sizes.

- Current Directory displays the current path or directory with all

files and subdirectories.

- Command History displays commands that have been entered.

Declarations and commands can also be entered in the form of a script

in a text file with ".m" extension . MATLAB is equipped with a text editor for entering
script files. The edit prog1 command opens the editor and loads the prog1.m file if it
exists, otherwise the editor opens with an empty file.
The following figure

Execution of the script (the commands one after one) is done using the Save and
Run button or with the debug/Save and Run menu or, simply, by pressing the F5
function key. The execution results are displayed in the execution window.
commands in the same manner as if commands are entered in this window.

Notes :
A command line can contain multiple statements separated by commas or

6
Machine Translated by Google

MEF course – Chapter 1: Presentation of MATLAB A. Seghir 2005-2014

with semicolons. The result of a statement followed by a semicolon will not be displayed.

1.2 Introduction of matrices

MATLAB essentially works with matrices. The components can be real, complex or
symbolic. Scalars are interpreted as matrices
1ÿ1!!

There are several ways to introduce a matrix into MATLAB:

1- Entry explicitly by list of elements
>> A = [ 1.1 1.2
2.1 2.2 ]

>> A = [ 1.1, 1.2 ; 2.1, 2.2]

The columns of the matrix are separated by commas (,) or spaces and the rows by
semicolons (;) or line breaks (RETURN).
2- Generated by an internal function (built-in function)
>>A = eye(3) % identity matrix
A=
1 0 0
0 1 0
0 0 1

>>B = ones(3) % square matrix containing 1s

B=

1 1 1
1 1 1
1 1 1

>>C = zeros(2,3) % zero rectangular matrix

C=

0 0 0
0 0 0

>>D = diag([1:3]) % matrix whose diagonal varies by 1

at 3
D=

1 0 0
0 2 0

7
Machine Translated by Google

MEF course – Chapter 1: Presentation of MATLAB A. Seghir 2005-2014

0 0 3
>> A = rand(2,3) % random matrix, a 2nd call of
rand

A= % returns another matrix 2ÿ3

0.7271 0.8385 0.3704

0.3093 0.5681 0.7027

3- Input from a script file with the editor

4- Loaded from a data file or imported from an application.
>>load matrix.txt

>>A = load('matrix.txt')

The matrix.txt file contains numerical values arranged in a table.

The first command loads the values of the matrix and assigns the matrix name to the
result, but the second allows you to choose a name for the loaded matrix.
The structure of the matrix.txt file is:

157

542

421

Spacings between components may be irregular. Commas or semicolons can be

used to separate columns and rows. An incomplete line causes a runtime error.

1.3 Operations on matrices

The following matrix operations are available in MATLAB
ÿ Addition C=A+B
ÿ

Substraction C=A-B
* C=A*B
Multiplication
^
Power C = A^2 or C = A * A
'
Transposed C = A' or C = transpose(A)
\ left division x = A\b

/ right division x = b/A

These matrix operations also apply to scalars. An error is caused in the case where
the ranks of the matrices are incompatible.

8
Machine Translated by Google

MEF course – Chapter 1: Presentation of MATLAB A. Seghir 2005-2014

Note 1
Matrix division involves solving systems of linear equations. If A is an invertible
square matrix:

x = A\b is the solution of the system A*x = b, with x and b are column vectors

x = b/A is the solution to the system x*A = b, with x and b now being row vectors

>> A = [3 2; 4 7]
A=
3 2
4 7
>> b = [2;3]
b=
2
3
>> x = A\b
x=
0.6154
0.0769
>>A* x
years =

2
3
y = b'
y=
2 3
>> z = y/A
z=
0.1538 0.3846
>>z * HAS

years =

2 3

Note 2
If a point is added to the left of the operators * ^ \ and / they will be applied to the
elements of matrices with index matching:
>>A*A % equivalent to A ^2
years =

17 20

9
Machine Translated by Google

MEF course – Chapter 1: Presentation of MATLAB A. Seghir 2005-2014

40 57
>>A.*A % equivalent to A.^2
years =

9 4
16 49

1.4 Vector operations

Vectors are considered rectangular matrices of one row or one column. All matrix
operations apply to vectors but care must be taken to match the sizes. Operators
preceded by a period are particularly useful in vector calculations and in graphical
representations.
Examples of vector manipulation:

>> a = [2 1 3] % line vector

a=
2 1 3

>> b = [1;3;2] % column vector

b=
1
3
2

>> a *b % scalar product

years =

11

>> a .* b' % product of components (note b')

years =

2 3 6

>> x = 0:0.2:1 % fixed step vector

x=
0 0.2000 0.4000 0.6000 0.8000 1.0000

ÿ
>> y = x .^ 2 .* cos(x) % y = x cos(x) (note the points)
y=

10
Machine Translated by Google

MEF course – Chapter 1: Presentation of MATLAB A. Seghir 2005-2014

0 0.0392 0.1474 0.2971 0.4459 0.5403

>> plot(x,y) % displays a window of the graph y(x)

1.5 Functions on matrices

MATLAB contains several functions that apply to matrices. Those that may interest us
currently are the following, to find out more type help or demo on the command line.

debt determining

>> d = det(A)

inv reverse

>> B = inv(A) (I = B * A)

eig eigenvalues and vectors

>> [Vects,Vals] = eig(A)

chol Cholesky decomposition

>> C = chol(A) (A = C'*C

read LU decomposition
>> [L,U] = read(A) A = L *U

qr QR decomposition
>> [Q,R] = qr(A) A = Q *R

svd decomposition into singular values

>> [U,S,V] = svd(A) A = U*S*V'

standard standard

Several definitions (type help norm)

Scalar functions such as cos, sin, exp, log…etc. apply to matrix elements

Example of calculating eigenvalues and vectors:

>> A = [5 2 1; 2 7 3; 1 3 8] ;

>> [vects,vals] = eig(A)

11
Machine Translated by Google

MEF course – Chapter 1: Presentation of MATLAB A. Seghir 2005-2014

vectors =
0.7040 0.6349 0.3182
-0.6521 0.4005 0.6437
0.2812 -0.6607 0.6959
values =
3.5470 0 0
0 5.2209 0
0 0 11.2322
>> vects' * A * vects
years =

3.5470 -0.0000 -0.0000

-0.0000 5.2209 -0.0000
-0.0000 -0.0000 11.2322

1.6 Access to matrix elements

Access to the elements of a matrix is done using their row and column indices
placed, after the name of the matrix, in parentheses: A(i,j) is the element of row i
and column j .
It is possible to obtain a submatrix by specifying a list of indices and two points in
place of an index or a list of indices are interpreted as an entire row or column:

>> A = [1.1 1.2 1.3 1.4

2.1 2.2 2.3 2.4
3.1 3.2 3.3 3.4
4.1 4.2 4.3 4.4]

>> A(2,2:4)
years =

2.2000 2.3000 2.4000

>> A(3,:)
years =

3.1000 3.2000 3.3000 3.4000

>> A(:,3)
years =

12
Machine Translated by Google

MEF course – Chapter 1: Presentation of MATLAB A. Seghir 2005-2014

1.3000
2.3000
3.3000
4.3000
>> t = [1,4]
t=
1 4

>> A(t,t) % This writing is useful in MEF

years =

1.1000 1.4000
4.1000 4.4000

1.7 Loops in MATLAB

The structure of loops in MATLAB does not differ from that of other programming languages. There are two types of

loops: for and while.

The for loop has the following structure:

for value = initial_value: not: final_value

statement 1
instructions 2
. . . . . . .
end

The while loop has the following structure:

initializations
while logical expression
statement 1
instructions 2
. . . . . . .
end

The for loop varies a variable from an initial value to a final value with a constant increment step which is taken equal

to one if omitted. The instructions inside the loop are executed at each value of the variable. The while loop

executes as long as a logical expression is true; it often requires initialization before being launched and can
run indefinitely if the controlling logical expression always remains true.

13
Machine Translated by Google

MEF course – Chapter 1: Presentation of MATLAB A. Seghir 2005-2014

Both loops can be broken by the break command. We generally insert break with an if condition test .

The following two loops give the same vector: f = w = 2 4 6 8 10

for i = 1:5 % for loop with unit step

f(i) = 2*i; % increment of i is automatic
end;

i = 1; % initialization of i
while i <= 5 % loop as long as i is less than or equal to 5
w(i) = 2*i;
i = i+1; % the increment is important otherwise the
end; % loop is infinite

1.8 Logic tests

Before the curls
Logic tests in MATLAB are done using the keywords if, else and end. The relational operators used in logical instructions

are:
< lower <= less than or equal
> greater == equal >= greater than or equal
~= not equal (different)

and the logical operators are:

& and example: c = a & b or c = and(a,b)
| or example: c = a | b or c = or(a,b)
~ no example: c =~a or else c = not(a)

Example of if test
if a < b
m=a
else
m=b
end

14
Machine Translated by Google

MEF course – Chapter 1: Presentation of MATLAB A. Seghir 2005-2014

1.9 Function files

Script files allow commands to be written and saved to a disk file. Function files are
more important, they provide extensibility to MATLAB. You can create a function
specific to a given problem and use it in the same way as the other functions offered
by MATLAB. A function file has the following structure:

function [res1, res2, ...] = functionname(var1, var2, ...)

% put here the comment lines to display

% in the command window following a help request

% with help functionname

instruction

instruction

the function keyword is essential. The input variables are var1, var2 ... and the results are res1, res2, ... As many variables and results as

necessary can be used. The file must be saved under the name functionname.m and may contain more than one function.

A return keyword can be added at the end of each function.

function R = rot3dz(theta)

%R = rot3dz(theta)
% returns a 3d rotation matrix

% of an angle around theta of the z axis % theta in degrees

measured from the x axis

theta = theta * pi/180; % transformation into radian

c = cos(theta); % cos and sin for greater speed

s = sin(theta);
R = [ cs 0
-sc 0

001

];
return

To have a rotation of 30° around z, we can call the function as follows:

>> r = rot3dz(30)
r=

0.8660 0.5000 0
-0.5000 0.8660 0

15
Machine Translated by Google

MEF course – Chapter 1: Presentation of MATLAB A. Seghir 2005-2014

0 0 1.0000

1.10 Application 1
1) Open MATLAB

2) Create two matrices A and B of 3ÿ3 in the command window

3) Create a column vector a and a row vector b of length 3

4) Calculate the transposes of matrices A, B and vectors a, b.

5) Calculate A + B and A – B.

6) Calculate A*B and A.* B, comment

7) Calculate the vectors x = A \ b and r = A * x.

8) Calculate the vectors x = a / A and r = x * A.

9) Calculate the matrices X = A / B and X = A \ B, check A*X, X*A, B*X and X*B

10) Calculate the matrices X = A ./ B and X = A .\ B, check A*X, X*A, B*X and X*B

11) Calculate the matrices R = A ^ 2, R = A .^2 and R = A .^B, comment

12) Calculate (A'\a')' and R = (b'\A')', comment

13) Calculate inv(A)*b and a*inv(A) , comment

14) Calculate a*b , b*a, b'*a, a'*b, b'*a', comment

15) Calculate a.*b , b.*a, b'.*a, a'.*b, b'.* a', comment

16) Create three functions for the three 3D rotations around the x, y and z axes.
Calculate a rotation matrix R composed of a rotation of 30° around x , 60°
around y and 45° around z. Then calculate Ar = R'*A*R.

17) Compare trace(A), trace(Ar), and eig(A), eig(Ar). Comment

18) Create a batch file using MATLAB editor and redo 2-10.

19) Create a text file containing the elements of a 3ÿ3 matrix and load the
matrix using a command line.

16
Machine Translated by Google

MEF course – Chapter 1: Presentation of MATLAB A. Seghir 2005-2014

1.11 Exercise

Create matrix A using the diag function and a vector b with a fixed step
A=
1 1 0
1 2 2
0 2 3
b=
1
2
3

Solve the system of equations A x ÿ b using:

1) The division operator,

2) The LU decomposition,

3) QR decomposition.

Use help for diag, LU and QR functions . Example help read.

1.12 Duty - TP

Write MATLAB scripts that allow you to:

1) Calculate the strain tensor as a function of the stress tensor in

using the relation ÿ ÿ (1 ÿ ÿ)/E ÿ ÿ ÿ/E s I.

2) Calculate the stress tensor as a function of the strain tensor using the
relation ÿ ÿ 2 ÿ ÿ + ÿ e I .

3) Calculate the relationships ÿ ÿ D ÿ and ÿ ÿ C ÿ using vector writing and

ÿÿ
elasticity matrix D, check that C ÿD _

4) Calculate the principal constraints and directions

5) Check that a rotation of 60° around y gives the main directions of the
stress tensor ÿ taken as an example below

6) Check the stress invariants and the strain invariants

17
Machine Translated by Google

MEF course – Chapter 1: Presentation of MATLAB A. Seghir 2005-2014

7) Calculate the stress vector and the strain vector on a normal n.

8) Calculate the normal stress and the unit expansion along the normal n.

9) Propose a structure of several function files for this type of problem

ÿ 7 0 ÿ

33 ÿ
ÿ ÿ

Take the stress tensor as an example: ÿÿ

ÿ
0 25 0 ÿ
MPa ,
ÿ
ÿ

330 13 ÿ

ÿ ÿ

The elastic modulus: E ÿ 200,000 MPa, Young's modulus: ÿ ÿ 0.25 and

1 1 1 T
The normal vector: n ÿ ÿ /ÿ3 /ÿ3 /ÿ3 ÿ

Reminder :

ÿ ÿ ÿ ÿ 2ÿ ÿ ÿ 0 0 0 ÿ

ÿ ÿ ÿ 2ÿ ÿ 0 0 0
ÿ

ÿ
ÿ ÿ ÿ2 0 0 0 ÿ

D ÿ
ÿ ÿ
;
2ÿ 0 0
ÿ ÿ

Sym 2ÿ 0 ÿ

ÿ ÿ

ÿ 2ÿ ÿ

ÿ 1 ÿÿÿÿ 0 0 0 ÿ
ÿ
1 0 0 0
ÿ

ÿÿ
ÿ ÿ

1 ÿ
1 0 0 0 ÿ

VS ÿ
ÿ ÿ

E 1ÿÿ 0 0
ÿ ÿ

Sym 1ÿÿ 0 ÿ

ÿ ÿ

1ÿÿ
ÿ ÿ

I1 ÿ ÿii ÿ ÿI ÿ ÿII ÿ ÿIII ;

I2 ÿ ½(ÿii ÿjj ÿ ÿij ÿji) ÿ ÿIÿII ÿ ÿIIÿIII ÿ ÿIIIÿI ;

I3 ÿ det(ÿ) ÿ ÿIÿIIÿIII.

ÿ ÿ ÿE / (1ÿ2ÿ)(1ÿÿ) ; 2ÿ ÿ E/(1ÿÿ)

18
Machine Translated by Google

chapter 2

General presentation of the

Finite Element Method

2.1 Introduction
To analyze a natural phenomenon in general or an engineering problem in particular,
we are often required to develop a mathematical model that can describe the problem
in question as reliably as possible.
The development of a mathematical model is generally based on a few basic
postulates and several simplifying hypotheses to arrive at governing equations which
are often differential equations to which boundary conditions are added. Example,
the theory of elasticity is based on the fundamental postula of the existence of the
stress vector and the general equations of isotropic linear elasticity are obtained with
the assumptions of small deformations, homogeneity and isotropy of materials as
well as the linearity of the relationships linking stresses and deformations.

The analytical resolution of differential equations sometimes poses insurmountable

difficulties, and an exact solution that describes the problem studied is not always
easy to find. The use of physical models and experimental simulation to find a
solution similar to the solution sought can prove costly in terms of time and resources.

With advances in the field of computing and increasingly greater computer

performance, it has become possible to solve very complex systems of differential
equations. Several digital resolution techniques
Machine Translated by Google

MEF Course – Chapter 2: Presentation of the Method A. Seghir 2005-2014

have been thus developed and successfully applied to provide satisfactory solutions to a
wide variety of engineering problems.
The finite element method is one of the most powerful numerical techniques.
One of the major advantages of this method is the fact that it offers the possibility of
developing a program to solve, with few modifications, several types of problems. In
particular, any complex form of a geometric domain where a problem is well posed with all
boundary conditions, can be easily treated by the finite element method.

This method consists of dividing the physical domain to be treated into several subdomains
called finite elements with non-infinitesimal dimensions. The sought solution is replaced in
each element by an approximation with simple polynomials and the domain can then be
reconstituted with the assembly or summation of all the
elements.

2.2 History of the method

The finite element method is the result of two areas of research:
mathematics and engineering sciences. [26]
Mathematics: Tools that go back to the weighted residuals of Gauss (1775), Galerkin (1915)
and Biezenokoch (1923), as well as the variational methods of Rayleigh (1870) and Ritz
(1909).
Engineering sciences: Whose contribution began in the 1940s with Hrenikoff (1941), Henry
(1943) and Newmark (1949) who touched continuous structures for the first time, making an
approximation on small portions in a continuous problem of a long bar. Hence the basic idea
of finite elements.

Argyris (1955), Turner (1956), Glough (1956) and Martin (1956) made a direct analogy by
adopting a simplified behavior for small portions: they represent a two-dimensional elastic
continuous medium by an assembly of triangular panels, on which the displacements are
assumed to vary linearly as for each bar or beam of the discrete system: each panel is
described by a

2
Machine Translated by Google

MEF Course – Chapter 2: Presentation of the Method A. Seghir 2005-2014

rigidity matrix and the assembly gave the overall rigidity of the continuous medium. Hence
the birth of finite elements with “panels” as their name.
Argyris and Kelsy (1960) use the notion of energy in the analysis of structures and use
mathematical methods (weighted residues, variational principles, etc.).
The term "finite element" was used for the first time by Glough (1960), and from then on, it
there is rapid development of the method.
In the sixties; Zienkiwicz (1965), De Arante (1968), Oliviera (1968), Green
(1969), Tones (1969), Lay (1969), Storne (1969), and Finlayson (1975) reformulated the
method from energetic and variational considerations in the general form of weighted
residuals, hence the mathematical model of the MEF.
In 1969 FEM was recognized as a general tool for solving PDE, and used to solve nonlinear
and nonstationary problems in several domains.
In fluid mechanics, the solution of the incompressible Navier Stokes equations by finite
elements using the velocity–pressure formulation began in the 1970s.

2.3 The main points of the method

In this paragraph, we will try to present in a simplified manner, the steps for applying the
finite element method and the tools necessary for its implementation.
The resolution of a physical problem by finite elements roughly follows the following steps
(Fig. 1.1):

Physical Problem

Domain discretization Mathematical Model

geometric

Choice of interpolation functions Variational formulation

Writing elementary matrices

Assembly of global matrices,

Application of boundary
conditions and System resolution
global
3
Figure 1.1: General steps of the finite element method
Machine Translated by Google

MEF Course – Chapter 2: Presentation of the Method A. Seghir 2005-2014

Step 1 : Formulation of the governing equations and boundary conditions.

The majority of engineering problems are described by partial differential equations
associated with boundary conditions defined on a domain and its contour. The application
of the MEF requires a rewriting of these equations in integral form. Weak wording is often
used to include conditions to
boundaries.

Step 2 : Division of the domain into subdomains.

This step consists of discretizing the domain into elements and calculating the connectivities
of each as well as the coordinates of its nodes. It thus constitutes the phase of preparation
of geometric data.

Step 3 : Approximation on an element.

In each element the variable such as displacement, pressure, temperature, is approximated
by a simple linear, polynomial or other function. The degree of the interpolation polynomial
is related to the number of nodes of the element. The nodal approximation is appropriate.
It is in this step that the construction of the elementary matrices takes place.

Step 4: Assembly and application of boundary conditions.

All the properties of the element (mass, rigidity, etc.) must be assembled in order to form
the algebraic system for the nodal values of the physical variables. It is at this level that
we use the connectivities calculated in step 2 to construct the global matrices from the
elementary matrices.

Step 5: Resolution of the overall system:

The overall system can be linear or non-linear. It defines either an equilibrium problem
which concerns a stationary or static case or a problem of critical values where it is necessary

4
Machine Translated by Google

MEF Course – Chapter 2: Presentation of the Method A. Seghir 2005-2014

determine the system's eigenvalues and vectors which generally correspond to the
eigenfrequencies and modes of a physical system.
A propagation problem which concerns the transient (non-stationary) case in which it
is necessary to determine the variations over time of physical variables and the
propagation of an initial value. Step-by-step integration methods are the most common
such as central finite difference method, Newmark method, Wilson method.
These methods must be associated with iteration techniques to deal with the non-
linear case. The most famous is Newton Raphson's method.

2.4 Variational formulation

Currently, the principle of virtual work is well known and very widespread. It is often
formulated in terms of the equality of work done by external and internal forces during
any virtual movement. This concept is essential for solving partial differential
equations. Indeed, the displacements are replaced by an arbitrary function continuous
on the domain and the equation is rewritten in integral form

2.4.1 Strong form

A classic problem of differential equations governing a physical system is stated as

follows:

Find a function u ÿ V; V space of functions, such that:

A(u) = 0 | ÿ ; B(u) = 0 | ÿ (2.1)

where A(u) is the set of governing equations defined on the domain ÿ and B(u) is the
set of boundary conditions that the functions u must satisfy on the contour ÿ
(fig.3.2). The function u can be a scalar such as temperature, pressure, etc. or a
vector such as displacement, speed, etc.
ÿ

ÿ B(u) ÿ 0

A(u) ÿ 0

Figure 2.2: geometric domain and contour

5
Machine Translated by Google

MEF Course – Chapter 2: Presentation of the Method A. Seghir 2005-2014

The variational problem associated with system (2.1) is written by taking the integral of the system of
governing equations weighted by weight functions, the statement becomes:

Find u ÿ V such that:

ÿ: w V ÿ ÿÿ w A ud
()ÿÿ0 (2.2)

This equation is called the strong integral form of the differential equation (or system of differential
equations). It is analogous to the expression of virtual works. In fact the solution of (2.2) has even more
scope, we can affirm that if it is satisfied for any weight function w, then the differential equation (2.1) is
satisfied at every point of the
domain ÿ.

(i) Example
We consider the following second order differential equation:

of2
Tou ( ) ÿ
ÿ 10ÿÿx ; 0ÿxÿ1 (2.3)
dx 2

Defined in the one-dimensional domain ÿ ÿ ÿ0, Lÿ with boundary conditions:

u( xÿ0 )ÿ0 and u( xÿ1 )ÿ0 (2.4)

In this case B(u) is the set of values imposed on the two edges of the domain. In one dimension, ÿ reduces
to two points. The integral form associated with the equation A(u)
is written:

2
of 1 of 2 1
w (1) xd
2 ÿÿÿÿ 0; ÿ w dx wÿ (x
ÿ 1) dx ÿ

(2.5)
ÿÿ 0 dx 2 0
dx

With the form of the terms to be integrated, we see that the search for approximate solutions for the
unknown function u requires the use of polynomials or interpolation functions differentiable at least twice.
In addition, the boundary conditions must be checked separately since they do not intervene in the integral
equation above, hence the introduction of the weak integral form.

6
Machine Translated by Google

MEF Course – Chapter 2: Presentation of the Method A. Seghir 2005-2014

2.4.2 Weak form

To satisfy the boundary conditions we have two ways of proceeding; either by the choice of the
weighting function, or check that:

ÿÿ w Bud
()ÿÿ0 (2.6)

In practice, it is possible to integrate (2.2) by part and replace it with:

ÿÿ C(w) D(u) dÿÿ E(w)

ÿÿ F(u) dÿ ÿ 0 (2.7)

The operators C, D, E and F contain derivatives of lower order, hence a wider choice of u
approximation functions.

This equation is the weak formulation of the differential equation, it forms the basis of the finite
element approximation.

Note: To obtain such integral formulations, we have two methods at our disposal: the first is the
method of weighted residues known as the Galerkin method, the second is the determination of
variational functionals which we seek to make stationary.

In practice, it is not always easy to find a functional one. The first process is more used; it
consists of choosing wÿÿÿÿu Dirac function (a disturbance of the desired function) and using the
nodal approximation for the discretization. w
is also called weight function hence the word: “weighted”.

(i) Example
To obtain the weak variational form of the previous example (equation 2.5) we integrate the first
term in parts.

1 dw of ÿ du ÿ 1 1
ÿ dx wÿ
ÿ ÿ 0 w (x 1) dx (2.8)
ÿ ÿ

0 dx dx ÿÿ
dx ÿÿ
0

We now see that the boundary conditions, particularly on the derivatives, are explicitly taken into
account in the equation. The values imposed on the function itself will be processed when
solving discrete systems.

7
Machine Translated by Google

MEF Course – Chapter 2: Presentation of the Method A. Seghir 2005-2014

2.5 Domain discretization

The finite element method is a subdomain approximation method, so before any application
it is necessary to divide the domain to be studied into elements. Each element is defined
geometrically by a well-determined number of nodes which generally constitute its vertices.
(Fig. 2.3)

Border
of the domain

Element

Border node
Internal node

Figure 2.3: Domain discretization – triangular elements

The geometric discretization must respect the following rules:

20) A node of an element must not be interior to one side of another of the same
type. (Fig. 2.4a)

21) No two-dimensional element must be flat, avoid angles close to 180°

or 0°. (Fig. 2.4b)

22) Two distinct elements can only have points in common located within their common
boundaries; recovery is excluded. (Fig. 2.4c)

23) The set of all elements must constitute a domain as close as

possible of the given domain; holes between elements are excluded. (Fig. 2.4d)

(has)
(d)
(b)
(vs)

8
Figure 2.4: Discretization rules
Machine Translated by Google

MEF Course – Chapter 2: Presentation of the Method A. Seghir 2005-2014

The result of the discretization process must contain two essential data which are the
coordinates of the nodes and the connectivities of the elements. We must number all
nodes and elements in order to have global matrices with a small bandwidth, for this, the
numbering is done according to the smallest width of the domain.

2.6 Approximation on the element

After defining the element, we can replace the exact function with an approximate one.
We often use polynomials or functions that are easy to implement on a computer.

2.6.1 Polynomial approximation and nodal approximation

The approximate function is expressed, in the one-dimensional case, by:

not

i
u axi (2.9)
ÿÿ iÿ 0

Which can be written in the following matrix form:

ÿ has 0 ÿ

2 ÿ ÿ

u ÿÿ 1 xx ÿ
ÿ ÿ has 0 ÿ
ÿÿ P x( )ÿAÿÿ (2.10)
ÿ ÿ

ÿ
ÿ ÿ

This form of approximation is called polynomial interpolation. If we express the function on

all the nodes we obtain for each node i with coordinate xi :

ÿÿ )
u i P(x i
ÿ ÿÿ
HAS
ÿÿ Pij a j (2.11)
j

Or for all nodes:

ÿ ÿ ÿP1ja { }j ÿ

ÿu ÿ .....
ÿ ÿ

ÿ
ÿÿ {..} ÿ
ÿ Un = Pn an (2.12)
ÿ
ÿÿ F nja { }j ÿ

ÿ ÿ

9
Machine Translated by Google

MEF Course – Chapter 2: Presentation of the Method A. Seghir 2005-2014

One : represents the values at the nodes of the function.

Pn : values of polynomials at nodes with coordinates xi.
an : generalized variables which are the factors of the polynomials.

The disadvantage of this interpolation lies in the use of the parameters ai as a base variable,
additional calculations are necessary to calculate the desired function u. In order to avoid these
calculations, we can put the values of the function u at the nodes as base variables by proceeding
as follows:

From equation (3.12), we can draw the ans as a function of ones and replace them in equation
(3.10). Which give :
ÿ

1
u ÿÿ P (x)ÿ PU ÿ ÿ N x ÿ
not not
() A (2.13)

This is the form most used because its variables are the values of the function at the nodes, the
resolution directly gives these values.

This type of approximation is called nodal interpolation, the Ni functions are called shape
functions, they depend on the type of element used for the geometric discretization.

10
Machine Translated by Google

chapter 3

Differential equations and

boundary problems

Evolution problem and boundary problem

3.1 First order differential equation

3.2 Steps to follow

In the case of differential equations, the problem is formulated directly in mathematical
form and amounts to determining an unknown function u defined in a domain ÿ
and governed by a differential equation with boundary conditions

3.2.1 Problem formulation

We take as an example the following ordinary differential equation of order one:

of
ÿ 2x(u ÿ1) ÿ 0 with u(0) ÿ 0
dx

In this problem, ÿ ÿ ÿ0, 2ÿ is a domain of dimension 1; its border at two ÿÿ is reduced

points: 0 and 2.
2
The exact solution to this equation is: u 1 e x
ÿ

ÿÿ

e
Machine Translated by Google

MEF course – Chapter 3: Differential equations A. Seghir 2005-2014

3.2.2 Domain discretization

The domain ÿ is divided into n segments (called elements) of size 1/n. Each element contains two
nodes on which the function u is interpolated.

1 Element 2 Element 3 Element 4 Element 5

1 2 3 4
xÿ0 x ÿ 0.5 x ÿ1 x ÿ 1.5 xÿ2

The division of the domain ÿ into several elements is called meshing. We use two tables to describe
the mesh: element connectivity table and node coordinate table. For an example of three elements
we obtain the two tables as follows:

Connectivity table Coordinate table

Element Node 1 (start) Node 2 (end) Node Coordinate ( x )

1 1 2 1 0.0
2 2 3 2 0.5
3 3 4 3 1.0
4 4 5 4 1.5
5 2.0

The MATLAB command lines which allow you to obtain the two tables are:
n = 4; x = % number of elements

0:2/n:2; % node coordinates

for i = 1:n % start of loop on elements

t(i,:) = [i, 1+i]; % connectivity of each element

end; % end of loop

or, in more compact form:

x = [0:2/n:2]'

t([1:n],:) = [[1:n]',[2:n+1]']

Noticed

12
Machine Translated by Google

MEF course – Chapter 3: Differential equations A. Seghir 2005-2014

For a higher number of elements, simply increase the value of nx(i) gives the coordinate of node i ; .
x(2) displays: 0.5

t(i,:) gives the two nodes of element i ; t(2,:) displays: 2 3

t(:,i) gives column i of table t ; t(:,2) displays: 2

3
4
5

3.2.3 Discretization and interpolation on the element

We can interpolate the function u sought in an element by a polynomial. The order of the
polynomial conditions the precision of the approximate solution. For an element with two nodes
we can take:

u ÿ aÿ ÿ aÿ x (3.)

u ÿ aÿ ÿ aÿ x u ÿ uÿ

u ÿ uÿ

Node 1 Node 2
x ÿ xÿ x ÿ xÿ

or in vector form:

ÿ has 0 ÿ
_ ÿÿu 1 x ÿ ÿ ÿ ÿ u ÿ p an
has
ÿ 1ÿ

not

with p row vector containing the monomials x

and a column vector containing the factors of the polynomial.

This approximation of the unknown function u is called polynomial interpolation, it is a function of

aÿ and aÿ which are coefficients without physical value. To use the values of u at the nodes we
look for an interpolation according to uÿ and uÿ.
The polynomial interpolation at the nodes is written:

13
Machine Translated by Google

MEF course – Chapter 3: Differential equations A. Seghir 2005-2014

ÿ u1 ÿ ÿ 1 x1 ÿ ÿ ÿ has 0 ÿ

ÿ ÿ
ÿ
ÿ ÿ
ÿ one ÿ Pn year
u2 x2
_ ÿÿ
has
1
ÿ ÿ ÿÿ1 ÿ

The inverse of this system of equations gives the parameters an.

ÿ has 0 ÿ 1 ÿ xxu 1 ÿ

ÿ
ÿÿ 2 1ÿ ÿ ÿ
an ÿ Pn Aÿ ÿ ÿ
ÿ
ÿ ÿ
has
1
xx 2 1 ÿ ÿ

1 1 u2
ÿ ÿ ÿÿ ÿÿ ÿ

By replacing the an we can now approach the function u by:

1 ÿ x2 _
ÿ

xu 1ÿ ÿ1ÿ ÿ xx 2 ÿ

xx 1
ÿ

ÿ u1 ÿ
u 1 ÿÿ ÿx _ ÿ ÿ
ÿÿ ÿ ÿ ÿ
xx 2 1
ÿ ÿ

1 1 u2 xx 2 1 ÿ

xx 2 1 ÿ

u2
ÿÿ ÿÿ ÿ ÿ ÿ

u ÿ N One

with N is a row vector containing functions of x called shape functions.

This interpolation is called nodal interpolation since it depends on the values at the nodes of the
unknown function u.

3.2.4 Properties of shape functions

It is interesting to note the following properties for functions of form N :

xx 2 ÿ

1 xx 1ÿ xx 1ÿ

1 xx 2ÿ
Nx()
1
ÿ ÿ
; Nx()
2
ÿ ÿ

xx 2 1 ÿ

xx 2ÿ xx 2 1 ÿ

xx 1ÿ
ÿÿÿ0 ÿÿÿ0

N 1(x) ÿ N (x)
2 ÿ1 ÿxÿ x ; x
ÿ
1 2
ÿ

1) They take the unity value at nodes with the same index and the zero value at
other nodes.

2) Their sum is equal to unity over the entire interval of the element:

1
N1

N2

x1 x2

14
Machine Translated by Google

MEF course – Chapter 3: Differential equations A. Seghir 2005-2014

Noticed
The two shape functions can be written as Jacobi polynomials:

xx ÿ

N i (x) ÿ
ÿ xx ÿ
j

ji ÿ ij

3.2.5 Lagrange type element

Thus elements for which the shape functions are Jacobi polynomials are called
elements of Jacobi type. We can directly construct, without going through polynomial
interpolation, shape functions for higher order elements (more than two nodes).

For example the shape functions of the element with three nodes are:

( xxxx
ÿ

2 )( ÿ

3 )
Nx()
1
ÿ
;
( xxxx 2 )(
ÿ ÿ

1 13 )

( xxxx
ÿ

1 )( ÿ

3)
Nx ()
2
ÿ
;
( xxxx 1 )(
ÿ ÿ

2 23 )

( xxxx )
ÿ ÿ

1 2 )
Nx () ÿ
3
( xxxx ( )( 3 2
ÿ ÿ

2 1 )

We can easily check in the figure below the properties at the nodes of the three
functions as well as their sum over the element.

15
Machine Translated by Google

MEF course – Chapter 3: Differential equations A. Seghir 2005-2014

ÿ1ÿ xx ÿ1 ÿ0ÿ xx ÿ1 ÿ0ÿ xx ÿ1

N x 1( ) ÿ
ÿ
0 xx ÿ2 ; N x 2( ) ÿ
ÿ
1 xx ÿ2 ; N x 3( ) ÿ
ÿ
0 xx ÿ2

xx ÿ3 xx ÿ3 xx ÿ3
ÿÿ0 ÿÿ0 ÿÿ1

The MATLAB script which allows you to calculate and plot these three functions is as follows:
x1 = 2;
x2 = 5;
x3 = 8;
x = x1:0.1:x3;
N1 = (x-x2).*(x-x3)/(x1-x2)/(x1-x3); N2 = (x-x1).*(x-x3)/(x2-
x1)/(x2-x3); N3 = (x-x1).*(x-x2)/(x3-x1)/(x3-x2);

S = N1+N2+N3;
plot(x,N1,x,N2,x,N3,x,S)

Note that it is preferable to use a more general-purpose function file.

3.2.6 Elementary matrices

The calculation of elementary matrices involves rewriting the problem in integral form:

ÿ of
ÿÿ ÿ u ÿ 2 xu
(
ÿ

1)ÿÿdÿ 0
ÿÿ
dx ÿÿ

with ÿu is a weighting function taken equal to a disturbance of the unknown function u.

The domain ÿ includes the interval 0 to 2, dÿ ÿ dx and with nodal interpolation we have:

of dN
ÿ
One ; and ÿu ÿ N ÿUn
dx dx

Since only the functions N depend on x and the disturbances only affect the values of u.

16
Machine Translated by Google

MEF course – Chapter 3: Differential equations A. Seghir 2005-2014

For convenience we write: ÿu ÿ (ÿ Nn TT

)N
The integral from 0 to 2 can be replaced by the sum of the integrals from xi to xi+1 (or: the integral
over ÿ is the sum of the integrals over ÿe, with ÿe is the domain of each element)

2 xi ÿ 1 x2

ÿÿ ÿ ÿÿ ÿ

elements e
ÿ
ÿ ÿÿ ÿ ÿ
0
ÿ

in 1,
ÿ
xi
ÿ

elements
x1

The integral form of the differential equation then becomes for each element:

x2 of ÿ x2 x2
ÿ
ÿ ÿu dx ÿ ÿ ÿ
ÿ uxu
2 dx ÿ ÿ

ÿ ÿ ÿ
ÿ u dx ÿ
0
x1
ÿÿ
dx ÿÿ
x1 x1

x2
ÿ TT dN ÿ x2
TT
x2
TT
ÿ ÿA not
U dx ÿ not
ÿ ÿ ÿ
2 dx
ÿ UN x NUnot not
ÿ ÿ

ÿ ÿ
A
ÿ2 xdx
not
ÿ ÿ
0
x1
ÿÿ
dx ÿ x1 x1

This discretized writing is valid for all types of elements. In the particular case of a linear element with
two nodes, it is written as follows:

x2 ÿ NOT ÿ ÿu ÿ
1 1
ÿÿ u1 ÿ ÿu
2 ÿ ÿ ÿ
ÿ ÿ dN dN dx
1 2 ÿ ÿ ÿ
x1
ÿ
NOT u
2ÿ ÿ 2ÿ

x2 ÿ NOT1 ÿ ÿu1 ÿ
ÿÿ u1 ÿ ÿu 2 ÿ ÿ ÿ
2 x NN
ÿ dx
1 2 ÿ ÿ ÿ
ÿ

x1
ÿ
NOT
2ÿ ÿ
u2 ÿ

x2 ÿ NOT ÿ
1
ÿÿ u1 ÿ ÿu ÿ ÿ ÿ
2 xdx =0
2
x1 ÿ
NOT
2ÿ

T
We see that it is possible to simplify (ÿa) since it is not zero and returns every time
term. By the way, it's the other terms that must be zero! we see it clearly with this example:

ÿv1 ÿ ÿw1ÿ
ÿ ÿ uu1 2 ÿÿ ÿ ÿ ÿ ÿ ÿ ÿ uu
ÿ1 2 ÿ ÿ
ÿ
0
ÿ
v2 ÿ ÿ
w2 ÿ

gives: ÿu1 v1 ÿ ÿu2 v2 ÿ ÿu1 w1 ÿ ÿu2 w2 ÿ 0 ; either : v1 ÿ w1 ÿ 0

17
Machine Translated by Google

MEF course – Chapter 3: Differential equations A. Seghir 2005-2014

ÿ v1 ÿ ÿ w1 ÿ

And v2 ÿ w2 ÿ 0 ; Or : ÿ ÿ
ÿ ÿ ÿ
ÿ
0 .
v2 w2
ÿ ÿ ÿ ÿ

Finally the discretized integral equation is put into matrix form:

ÿ x2
_
ÿ #1 ÿ x2
_
ÿ #1 ÿ ÿÿ u1 ÿ

ÿ ÿ ÿ
ÿ d N d N dx ÿÿ ÿ ÿ ÿ ÿ 2 xÿ NN 1dx ÿ ÿ ÿ2 ÿ ÿ ÿ
ÿ

1 2
x1 NOT x1 NOT u2
ÿÿÿ ÿ 2 ÿ ÿ 2 ÿ ÿ

x 2 ÿ Nÿÿ
1
ÿ _ ÿ ÿ 2 x dx = 0 ÿ
x1 No. 2

ÿ x2 ÿ N d NN d N ÿ ÿ x2 ÿ 2N
N xN 2 xN 2 ÿ ÿ ÿ dx ÿÿÿ ÿ ÿ u1 ÿ ÿ2x 2 x N 1
_
1 1 1 2 ÿ dx
_
1 1 1
ÿ ÿÿÿ ÿ ÿÿ
ÿ = ÿ ÿ ÿ dx ÿ
N d NN d N 2 2 N xN2 N
2 xN 2 2 2 u2
ÿÿÿ
x1
2 1
x1
ÿÿ 1 ÿ
x1
ÿ ÿ2xN _ 2

Which is a system of linear equations Ke Ue ÿ Fe ; with Ke and Fe are called elementary matrix and vector of the

equation system. In the case of the present differential equation K is the sum of two matrices: Ke ÿ Keÿ ÿ Keÿ such

that:

x2
_
ÿ N 1dN N dN 1ÿ
1 2ÿ x2
_ N2 xN
ÿ2 N 1xN 1 ÿ1 ÿ 2ÿ
Ke 1 ÿ
ÿ dx ÿ ;
ÿ
Ke 2 ÿ
ÿ dx ÿ
x1 N 2dN N
1 dN 2ÿ 2
x1 N xN
2 2N 1xN 2 2 2
ÿ

By replacing the shape functions and their derivatives by their respective expressions we obtain:

xx 2 1 xx 2 1
ÿ ÿ

ÿÿ ÿ

x2
xxxx 2 1 2 xxxx 2 2 1
ÿ ÿ ÿ ÿ

_
1 1
Ke
1
ÿ
ÿ ÿ ÿ dx ; ÿ ÿ ÿ
xx 1 1 xx 2 1
ÿ ÿ

x1

xxxx 1 1 2 xxxx 1 2 1
ÿ ÿ ÿ ÿ

ÿÿ 2 2

xx 2 xx 2 xx 2 xx 1
ÿ ÿ ÿ ÿ

ÿÿ
ÿ
2 x 2 x
x2
_
xx 1 2 xx 1 2 xx 1 2 xx 2 1
ÿ ÿ ÿ ÿ

Ke 2
ÿ ÿ ÿ ÿ dx ÿ
x1 xx 1 xx 2 xx 1 xx 1
ÿ ÿ ÿ ÿ

2x _ 2x _
xx 2 1 xx 1 2 xx 2 1 xx 2 1
ÿ ÿ ÿ ÿ

ÿÿÿÿÿÿ
ÿÿ

Or after integration of the components of the two matrices:

1 ÿ
ÿ

11
ÿ
Ke1 ÿ ;
2 ÿ
ÿ

ÿÿ1
1ÿ

18
Machine Translated by Google

MEF course – Chapter 3: Differential equations A. Seghir 2005-2014

xxxxxx
2
ÿ

1
ÿ3 12
ÿ
2 1
ÿ ÿ
Ke 2 ÿ

6 ÿ xxxx 3
ÿ ÿ

ÿÿ
112 2 ÿ

The vector F is given by:

x2 ÿ 2xxxx ÿ

)(
ÿ

)ÿ
Fe ÿ
ÿ (ÿ
212
ÿ dx
ÿ2
xxxx
( 1) ( )ÿ
x1
ÿ ÿ

21

either :

xxxx
2
ÿ

1 ÿ
2 12
ÿ ÿ
Fe ÿ
ÿ ÿ
3 ÿ
xx1 2 ÿ 2 ÿ

3.2.7 Note
It is possible to use MATLAB to analytically integrate elementary matrices. The script can
be as follows:
syms x x1 x2 real symbolic % declaration of variables

N = [(x-x2)/(x1-x2) (x-x1)/(x2-x1)] % shape functions

dN = simple(diff(N,x)) % derived from functions of
shape

Ke1 = simple( int(N' * dN , x, x1, x2) ) % matrix Ke1

Ke2 = simple( int(N' * 2*x * N , x, x1, x2) ) % matrix Ke2
Fe = simple( int(N' * 2*x , x, x1, x2) ) % Fe vector

3.3 Step 4: Assembly

The calculation of elementary matrices makes it possible to obtain the following systems
of elementary equations for the four elements:
Element 1 : xÿ ÿ 0; x2 ÿ ½ ;

1 1 ÿ 1 1 1
ÿ 0 ÿ 3 ÿ0 2ÿ 0 ÿ ÿ
ÿ

2 2 2 2
Ke1
(1) ÿ
; Ke 2(1) ÿ
;
1 1 ÿ 21 1
6 ÿ 0 03
ÿ ÿ ÿ ÿ
ÿ

ÿÿ
ÿ 2 2ÿ 2ÿ

19
Machine Translated by Google

MEF course – Chapter 3: Differential equations A. Seghir 2005-2014

11 13 ÿÿ
ÿ ÿ

(1) (1) 1 ÿ 0.4583 0.5417 ÿ 0.4583 ÿ

(1) Ke Keÿ Ke ÿ1
2
ÿ

ÿÿ ÿÿ ÿ ;
24 ÿ

11 15 ÿ ÿ ÿ

0.6250 ÿ

12
ÿ

ÿ ÿ ÿÿ 2ÿ 0 0 12 ÿ 1 ÿ 1 ÿ ÿ ÿ 0.0833
(1) Fe ÿ
ÿ
ÿ
ÿÿ ÿÿ
ÿ
ÿ ÿ
3 0 2ÿ ÿ 12 2 0.1667 ÿ
12 ÿ ÿ

u1
ÿ

(1) (1) (1) 1 ÿ ÿ ÿ 11 13 ÿ 1 ÿ 1 ÿ

Ke Ue ÿ Fe ; ÿÿ ÿÿ ÿ
ÿ
ÿÿ ÿÿ
24 ÿ

11 15 u2 12 2
ÿÿ ÿ

We denote by U the values of the function u at the five nodes, the values of the
elementary vector Ueÿÿÿ of element 1 correspond to the components uÿ and uÿ global
vector U, those of Ueÿÿÿ of the element 2 correspond to the second and third of the
global vector U, those of Ueÿÿÿ to the third and fourth and those of Ueÿ4ÿ to the fourth
and fifth.
Element 2 : x1 ÿ ½; x2 ÿ 1;

(2) 1
ÿ

15 ÿÿ ÿ ÿ ÿ 7
ÿ

ÿ 0.2917 0.6250 ÿ ÿ 1 ÿÿ4 ÿ ÿ 0.3333

; Fe (2)
ÿ
ÿ ÿ
Ke ÿ
ÿÿÿ ÿ ÿ
ÿÿ 12 5 0.4167
24 ÿÿ
ÿ

9 19 ÿ

0.3750 0.7917 ÿ ÿ

(2) (2) (2) 1 ÿ

ÿ

7 15 u2 ÿ 1 4
ÿÿÿ
Ke Ue ÿ Fe ; ÿ ÿ
ÿ
ÿ 12 ÿ
ÿ 24 ÿ

ÿÿÿ9 u 5ÿ
ÿ 19 ÿ ÿ 3 ÿ

Element 3 : x1 ÿ 1; x2 ÿ ÿ /ÿ ;

1 ÿ
ÿ

3 17 ÿ ÿ ÿÿ
ÿ

0.1250 0.7083 1 ÿ 7 ÿ 0.5833

ÿÿÿÿ
(3)
Ke ÿ

ÿ7
; Fe (3) ÿ
ÿ ÿ ÿ 12
ÿ
0.6667 ÿ
ÿ 24 ÿ

23 ÿ ÿ ÿ

ÿ ÿ 0.2917 0.9583 ÿ ÿ8ÿ

ÿ

1 ÿ
ÿ

3 17 ÿÿ u3 ÿ 1 7
ÿÿÿ
(3) (3) (3)
Ke Ue ÿ Fe ; ÿ
ÿ
ÿ 12 ÿ
ÿ 24 ÿ

ÿ ÿ 7 23 u4 8ÿ
ÿ ÿÿ ÿ

ÿ
Element 4 : x1 ÿ /ÿ ; x2 ÿ 2 ;

(4) 1 ÿ 1 19 ÿ ÿ ÿÿ 0.0417 0.7917 1 ÿÿÿÿ10 ÿ 0.8333

ÿÿÿÿ
Ke ÿ

ÿ5
; Fe (4) ÿ
12 ÿ 11 ÿ
ÿ
0.9167 ÿ
ÿ 24
ÿ
ÿ

27 ÿ ÿ ÿ

ÿ ÿ 0.2083 1.1250 ÿ

20
Machine Translated by Google

MEF course – Chapter 3: Differential equations A. Seghir 2005-2014

1 ÿ 1 19 ÿ ÿ ÿ ÿ u 4 ÿ 1 ÿ ÿ 10
(4) (4) (4)
Ke Ue ÿ Fe ; ÿÿ ÿÿ ÿ
ÿ
ÿÿ ÿ
24 ÿ

5 27 u5
ÿ
12 11 ÿ

By rewriting the elementary systems according to all the components of U we

gets:

ÿ
ÿ

11 13 0 0 0 0 ÿ ÿ 11 15 0 0 ÿ ÿ u 1 ÿÿ ÿÿÿÿ2
1
00 ÿÿÿÿÿ
ÿ

ÿÿÿÿÿÿ0000000ÿÿ u2 ÿ0

1 00 1 12 ÿ ÿ 0
ÿÿÿÿÿ

Element 1 : 000ÿÿÿÿ0000ÿÿÿÿ 0 u3 ÿÿ
ÿ
0ÿ

ÿÿÿÿ ÿÿ
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ u4
24 ÿÿ
ÿÿÿÿ ÿ u5
ÿÿÿ ÿÿÿÿÿ

ÿ ÿÿ ÿÿ ÿ0ÿ0009019
0 ÿ0ÿ0ÿÿÿÿÿ ÿ u1 ÿÿ ÿ ÿ 0ÿ ÿÿ ÿÿ ÿ ÿ

ÿÿÿÿÿÿÿÿ00000ÿ 4
0 7 15 0 0
ÿÿÿÿÿÿÿ
ÿ

u2
1 ÿ 1 ÿ0ÿ

Element 2 : ÿ

u3 ÿÿ
ÿ
5
24 12
00000 u4 0
u5
ÿÿÿÿÿ

ÿ ÿÿ ÿÿ ÿ0ÿ000003017 ÿ ÿ0ÿÿÿÿÿ ÿ u1 ÿÿ ÿÿ0

ÿÿÿÿ

ÿ 0ÿ 0ÿ 0ÿ 0
ÿ ÿ0 ÿ ÿ 0 0 0 0 0 ÿ u2
ÿ1ÿÿ7ÿ
0
1 ÿ ÿÿÿÿÿÿ
0ÿ

Element 3 : ÿ

u3 ÿÿ
ÿ

24 12
0 0 7 23 0 ÿ

u4 8
u5
ÿÿÿÿÿ

ÿ ÿÿ ÿÿ ÿ0ÿ0ÿ000000ÿ0ÿ0ÿ ÿÿ ÿ u1 ÿÿ ÿ ÿ 0 ÿÿ 1ÿ ÿ ÿ ÿ

ÿ 0ÿ 00 00 00 01 19 ÿ ÿ ÿ ÿ ÿ u2 0
ÿ ÿ ÿ ÿ 10 ÿ ÿ

1 ÿ ÿ 0 0 0 5 27 ÿ ÿ ÿ ÿ ÿ 11 ÿ

Element 4 : u3 ÿÿ
ÿ
0
24 12
u4

u5
ÿÿÿÿÿ

Now taking the sum (ÿ sum of the integrals), the global system is finally written:

21
Machine Translated by Google

MEF course – Chapter 3: Differential equations A. Seghir 2005-2014

11 13 0 0 0 ÿ ÿu 1 ÿ 1
ÿ

ÿ ÿ ÿ
ÿ ÿ ÿ ÿ ÿ ÿ
ÿ

1115715 ÿ

0 0 u 2 4ÿ
2
1 1
ÿ ÿ ÿ

ÿ KU ÿ F
ÿ ÿÿ ÿÿ ÿÿ

0 9193 170 ÿ u 3 ÿ ÿ 5 7ÿ
ÿ ÿ
ÿ
ÿ ÿ
ÿ
24 12
ÿ ÿ

ÿ u 8 1ÿ0
ÿ ÿ ÿ ÿ

0 0 ÿ

723119
ÿ ÿ
4
ÿ ÿ ÿ ÿ

ÿ
0 0 0 ÿ

527 ÿ u 11
ÿ ÿ ÿ 5ÿ ÿ ÿ
ÿ ÿ ÿ ÿ

This operation is called assembly. In practice, we do not proceed in this way for reasons
of saving memory and calculation time but we assemble the elementary matrices using
the connectivities of the elements. The global matrix K is first initialized to the zero matrix,
then at each elementary matrix construction, we locate with a table where it must be
added to the global matrix. In the case of one degree of freedom per node (this case) the
location table corresponds to the connectivity table. Example ; for element 1 the location
table is t ÿ ÿ1 2ÿ we therefore add the matrix Keÿÿÿ to K(1,1) , K(1,2) , K(2,1) and K(2 ,2);
and for element 2 the location table is t ÿ ÿ2 3ÿ we add the matrix Keÿ2ÿ

to K(2,2), K(2,3), K(3,2) and K(3,3). And so on.

The MATLAB script that allows you to do this type of assembly is simple:
K = zeros(n+1); % matrix initialization

global (n+1 nodes)

F = zeros(n+1,1); % vector initialization

global (note ,1)

for i = 1:n % loop over elements

t = [i i+1]; % location table

each element
x1 = x(i); x2 = x(i+1); % coordinates of the two nodes

of the element
[Ke,Fe] = MatElt2Nd(x1,x2); % Function to calculate the
mast. and vectors. elementals
K(t,t) = K(t,t) + Ke; % Die assembly

global
F(t) = F(t) + Fe; % Assembly of the global vector
end; % end of loop

22
Machine Translated by Google

MEF course – Chapter 3: Differential equations A. Seghir 2005-2014

3.4 Step 5a: Resolution – Application of CALs

Before solving the system it is necessary to apply the boundary conditions of the function u.
At node 1 with coordinate x ÿ 0, u ÿ 0; which results in uÿ ÿ 0. This induces the reduction of
the total number of equations to solve. The overall system then becomes:

ÿ 15 7 15 ÿ 0 0 ÿ ÿu 2 ÿ 2 4ÿ
ÿ

ÿ ÿ

9 19 3 17 0 u3 5 7ÿ
ÿ ÿ ÿ ÿ ÿ

1 1
ÿ ÿ

ÿ ÿ ÿ ÿ

ÿ ÿ ÿ
ÿ ÿ ÿ ÿ
24 ÿ 0 ÿ

7 23 1 19 ÿ ÿ
u4 12 8 10
ÿ
ÿ ÿ ÿ ÿ

ÿ ÿ

0 0 ÿ

5 27 ÿ
u5 ÿ ÿ
11 ÿ

ÿ ÿ ÿ ÿ ÿ ÿ

The row and column with index 1 which correspond to the value u1 have been deleted from
the matrix K and the vector F. If u(0) ÿ a ÿ 0 (uÿ ÿ a) then we subtract from F the product of
the
time
1
column of K by a.
The MATLAB command which allows you to delete a row or column from a matrix consists of
assigning it a zero (empty) vector:
A(i, :) = [] % deletes line i from the
matrix A
A(:, i) = [] % deletes column i
V(i) =[] % removes component i from
vector V

To apply the boundary condition uÿ ÿ 0, simply write:

K(1,:) = []; % deletion of 1st line

of K
K(:,1) = []; % deletion of the 1st
column of K
F(1,:) = []; % deletion of the 1st
component of F

Note: for u(0) ÿ a, we insert, before the deletions, the command: F = F -

K(:,1)*a.

23
Machine Translated by Google

MEF course – Chapter 3: Differential equations A. Seghir 2005-2014

3.5 Step 5b: Resolution – Calculation of the solution

The solution can now be calculated by one of the linear system resolution algorithms.
Note that in this case the matrix K is tridiagonal, but with MATLAB it is enough to
use the left division U = K \ F. The results are given in the following table and figure:

Coord. Solution Solution Error

x Exact MEF (%)

0.5000 0.2212 0.2488 12.4555

1.0000 0.6321 0.6673 5.5705

1.5000 0.8946 0.9154 2.3225

2.0000 0.9817 0.9843 0.2694

3.6 Study of convergence

Any study that uses an approximate numerical solution must be subject to a
convergence examination. In finite element studies, convergence can be achieved
by increasing the degree of the interpolation polynomials (number of nodes per
element), by using a fairly large number of elements with refinement of the mesh or
by seeking a good adaptation of the mesh to the problem addressed. The figure
after gives the percentage of error relative to the point x ÿ 0.5 as a function of the
number n of elements used. We see that the decrease in the relative error is a function of

24
Machine Translated by Google

MEF course – Chapter 3: Differential equations A. Seghir 2005-2014

1/ n and tends to cancel out in a monotonous way. A number of elements n ÿ 30

type
is enough to have an acceptable error.

The following MATLAB program combines all the previous steps.

function [U, Ue, Err, x] = ExampleEquaDiff(n)
%

% du/dx - 2*x*(1-u) = 0
% u(0) = 0
%

%------------------------------------------------- ---

% Solution with interpolation functions at 2Nds

%------------------------------------------------- ---

x = [0:2/n:2]'; % vector of coordinates

K = zeros(n+1) ; % initializations

F = zeros(n+1,1);

for i = 1:n % loop over elements

j = i+1;
t = [ij];
x1 = x(i);

x2 = x(j);
[Ke, Fe] = MatElt2Nd(x1,x2);
K(t,t) = K(t,t) + Ke;
F(t) = F(t) + Fe;
end;

25
Machine Translated by Google

MEF course – Chapter 3: Differential equations A. Seghir 2005-2014

K(1,:) = []; % application of CALs

K(:,1) = [];

F(1,:) = [];

U = K\F; % resolution

U =[0;U]; % incorporation of value

initial for plot

Ue = 1-exp(-x.^2); Err = 100*(U- % calculation of the exact solution

Ue)./Ue; % relative error calculation

plot(x,U,x,Ue) % plots both solutions

return

%--------------------------------------------

% Calculation of the Ke matrix and the Fe vector

%--------------------------------------------

function [Ke, Fe] = MatElt2Nd(x1,x2) % declaration of the function

Ke1 = 1/2*[ -1 1 % matrix Ke1

-1 1 ] ;

Ke2 = (x2-x1)/6* [ 3*x1+x2 x1+x2 % matrix Ke2

x1+x2 x1+3*x2];

Ke = Ke1 + Ke2; % matrix Ke

Fe = (x2-x1)/3 * [2*x1+x2 ; x1+2*x2]; % Fe vector

return

The command: [U, Ue, Err, x] = ExampleEquaDiff(30) allows you to have both the
finite element solution U, the exact solution Ue, the relative error Err and the
coordinates x with 30 elements.
Note :
Note that this program practically allows you to solve all differential equations of
order 1 with the boundary condition u(0) ÿ 0. Only the function MatElt2Nd(x1,x2)
requires a change for elementary matrices.

26
Machine Translated by Google

MEF course – Chapter 3: Differential equations A. Seghir 2005-2014

3.7 Duty No. 2

1) Modify the program ExampleEquaDiff(n) to solve the same equation using quadratic form functions
(three-node element)

2) Write a MATLAB function which allows you to have the expressions of the functions
1D shape .

3) Write a MATLAB finite element program which allows you to solve in the
of u(0) ÿ1
domain [0, 3] the equation: dx ÿ (2x ÿ1) u ÿ 0 with . Compare with
2
xx
the exact solution:
ÿ

eue ÿ
and examine the convergence at the point x ÿ 0.5.

3.8 2nd order differential equations

As an example of a 2nd order equation , we will treat in the domain ÿ ÿ [0, 1], the equation:

2
of of
ÿ uxx
69ÿÿÿ (1 ) with the boundary conditions: u(0) ÿ 0 And
2
dx dx
ÿdx ÿ
xÿ 1
ÿ
0.

3x 3
(8 ÿ ÿ xe 12)x
The exact solution is: u xx)(3 4) ÿ ÿ e 5 4 1
(1
ÿ

ÿ ÿ

e 127

The domain ÿ will be subdivided (mesh) into n linear elements with 2 nodes. The strong integral form of
this equation is written:

1 of 2 1 du 1 1
ÿ ÿu 2
dx 6 ÿ ÿÿ u dx 9ÿÿuuÿdx ÿ ÿ u x(1 x) dx
ÿ ÿ

0 dx 0 dx 0 0

When dealing with differential equations of second order and higher, we are faced with two difficulties.
Interpolation of derivatives and satisfaction of boundary conditions. In fact, the greater the order of the
derivatives, the higher the degrees of the polynomials or interpolation functions to use must be. A linear
element with two nodes cannot therefore be used for this type of equation. Moreover, the strong integral
form does not reveal the condition
ÿdx ÿ
xÿ 1
ÿ
0.

27
Machine Translated by Google

MEF course – Chapter 3: Differential equations A. Seghir 2005-2014

These two difficulties led to the use of the weak integral form which can be obtained, in this case
with integration by part of the first term:
1
1 of2 1 ÿ of of of ÿ ÿ
ÿ ÿu dx ÿÿ
ÿ dx ufrom ÿ ÿ
0 dx 2 0 dx dx dx ÿÿ ÿÿ
0

We immediately see the major advantage that this expression offers. It reduces the order of
the derivatives (hence the weak term) and allows the boundary condition to be taken into account
ÿdx ÿ
xÿ 1
ÿ
0.

Moreover since u (0) ÿ 0, the term ÿu also vanishes at x ÿ 0 (there cannot be any disturbances in
known or zero values).

Thus the integral form is written taking into account the boundary conditions as follows:

1 ÿ of of 1 1 1
ÿ

ÿ dx 6 uÿ ÿ ÿ from dxÿ9 ÿ ÿ x) dx
ÿ uu dx ÿ ÿ u x(1 ÿ

0 dx dx 0 dx 0 0

Its discretization gives for the integrals on the elements:

x2 T x2 x2 x2
dN dN T dN T T
ÿ

ÿ 6ÿ ÿ NOT ÿ ÿ dx N xxdx
(1 dx9 NN ÿ
ÿ ÿ

)
x1 dx dxdx x1 dx x1 x1

With this writing it is possible to use a linear element; the elementary matrix Ke
is the sum of three matrices: Ke ÿ Keÿ ÿ Ke2 ÿ Ke3 , with :

x2 ÿ1 xx
12
ÿ

)ÿ 1 ÿ
ÿ

1 1ÿ
Ke
1
ÿÿ
ÿ (ÿ ÿ
ÿ 1 ( xx 1 2 ÿ

) 1 ( xx dx 2 1 ÿ ÿ ) ; Ke 1
ÿ

x1
ÿ1 (
xx
21
ÿ

) ÿ
xx2 1 ÿ
ÿ

ÿ 1 ÿ

1ÿ
ÿ

x2 xxxx ÿ ÿ

ÿ
ÿ

33 ÿ
ÿ 2 )( 12 )ÿ
Ke 2 ÿ
6ÿ (ÿ ÿ
ÿ 1 ( xx 1 2 ÿ

) 1 ( xxdx 21
ÿÿ ) ; Ke2
ÿ

33
ÿ ÿ

x1
ÿ

xxxx ÿ ÿ

ÿ ÿ
ÿ( 1) ( 21 ) ÿ

x2 ÿ xxxx
ÿ

)(
ÿ

)ÿ
2
Ke 3 ÿ
9ÿ (ÿ
12
ÿ
ÿ ÿ (xxxx
12
2 )(
ÿ

) ( ) xxxxdx1 ) (
ÿ

2
ÿÿ
1 ;
x1 xxxx
1) ( )ÿ
ÿ ÿ

ÿ 21

xx2 1 ( ÿ ÿ 6 3 ÿ
Ke 3 ÿ

2 ÿ3 6 ÿ
ÿ ÿ

28
Machine Translated by Google

MEF course – Chapter 3: Differential equations A. Seghir 2005-2014

and the elementary vector Fe is:

x2 ÿ xxxx ÿ

)(
ÿ

)ÿ
2 12
Fe ÿ
ÿ (ÿ ÿ xx(1)
dx
ÿ

;
x1
ÿ(
xxxx
1) ( )ÿ
ÿ ÿ

21

xx 2 1 ÿ

ÿ (4 3
ÿÿ

xxx 2 12)
1
ÿ ÿ (2 xx )
22 ÿ
Fe ÿ
ÿ ÿ
12 ÿ (4 3
ÿÿ

x2
_ 2)xxxx ÿÿ
1 2 1 1 (2 ) ÿ

3.9 MATLAB program

Note the modifications introduced in the previous program to process this equation.

function [U, Ue] = EquaDiff2(n)

%--------------------------------------
% d²u/dx² + 6 du/dx + 9 u = x(1-x)
% u(0) = 0 du(1)= 0
%--------------------------------------
x = [0:1/n:1]'; % modification of 1 terminal

integration
K = zeros(n+1) ;
F = zeros(n+1,1);
for i = 1:n

j = i+1;
t = [ij];
x1 = x(i);
x2 = x(j);
[Ke,Fe] = MatElt2Nd(x1,x2);
K(t,t) = K(t,t) + Ke;
F(t) = F(t) + Fe;
end;
K(1,:) = [];
K(:,1) = [];
F(1) = [];

U = K\F;
U = [0.0;U];

t = 0:0.01:1; Ue =
1/27*(1-t).*(3*t-4)+1/54*exp(-3*t).*(8+t*exp(3)-12*t);

29
Machine Translated by Google

MEF course – Chapter 3: Differential equations A. Seghir 2005-2014

plot(x,U,'-.',t,Ue)
return

%--------------------------------------------

% Calculation of the Ke matrix and the Fe vector

%--------------------------------------------

function [Ke,Fe] = MatElt2Nd(x1,x2)

Ke1 = 1/(x2-x1)*[ -1 1 % changes do not

touch

1 -1 ] ; % essentially that the

matrices

Ke2 = [ -3 3 % elementary

-3 3 ] ;

Ke3 = (x2-x1)* [ 3 3/2

3/2 3];

Ke = Ke1 + Ke2 + Ke3;

Fe = (x2-x1)/12* [(4-3*x1-2*x2)*x1+(2-x2)*x2;

(4-3*x2-2*x1)*x2+(2-x1)*x1 ];
return

A call to this program with the EquaDiff2(15) command gives the figure
next :

30
Machine Translated by Google

MEF course – Chapter 3: Differential equations A. Seghir 2005-2014

3.10 Comparison with step-by-step integration

Compared to step-by-step integration methods which deal with initial value problems
(or evolution problems), the finite element method deals with boundary value problems.
This results in the boundary conditions or initial conditions required by the two methods.

For comparison we will treat the differential equation

of2
0ÿuÿxÿ
dx2

of
with the boundary conditions: u(0) ÿ 0 And
(4ÿ) ÿ 0 for the MEF
dx
and the initial conditions:

of of 2
u(0) ÿ 0 , (0) ÿ 0 And 2
ÿ
(0) 0 dx for step by step integration
dx

The exact solution being ue ÿ x ÿ sin x. The MATLAB command that provides this
solution is:
>> ue = simple(dsolve('D2u + u -t','u(0)=0','Du(4*pi)=0'))

The MATLAB function for step-by-step integration of one-degree-of-freedom systems

with the linear acceleration method is used to solve the equation numerically. A
detailed account of this method can be found in the work on structural dynamics by
Penzien and Clough.
function [u,v,a] = SdofTH(m,c,k,acc,dt)
% [u,v,a] = SdofTH(m,c,k,acc,dt)
% u, v, a: history of the response in displacements u,
% speeds v and accelerations a of a 1DOF system
% m, c, k: mass, damping and stiffness
% acc : acceleration at the base in vector
% dt : no time
%

% A. Seghir, 08/19/04

u(1) = 0.0;
v(1) = 0.0;
a(1) = - acc(1);

31
Machine Translated by Google

MEF course – Chapter 3: Differential equations A. Seghir 2005-2014

n = length(acc);
K = k + 6/dt^2 * m + 3/dt * vs;
for i = 1:n-1
A = 6/dt * v(i) + 3.0 * a(i);
B = 3.0 * v(i) + dt/2 * a(i);
P=m * ( acc(i) - acc(i+1) + A ) + B * vs;
du = P/K;
dv = 3/dt * du - 3.0 * v(i) - dt/2 * a(i);
u(i+1) = u(i) + du;
v(i+1) = v(i) + dv;
P = -c * v(i+1) - k * u(i+1) - m * acc(i+1);
a(i+1) = P/m;
end
return

The MATLAB commands that allow you to plot the two solutions are as follows:
clear
clc
n = 20 % number of integration steps
or finite elements
dt = 4*pi/n; t = % no integration time
0:dt:4*pi; % discretization of time
acc = -t; % inertia force Fi = - m *

acc
up = sdofth(1,0,1,acc,dt); step-by-step % function call
integration

[um,xm] = EquadiffMEF(n); % finite element solution

x = 0:0.01:4*pi; % x coordinates to plot the
exact solution
ue= x-sin(x); % exact solution
plot(t,up,'.',xm,um,'o', x,ue) % plots the three solutions

Note : The EquadiffMEF(n) function has the same structure as EquaDiff2(n). The
difference lies only in the calculation of the elementary matrices and vectors.
Also note the output variables!
The plot of the three solutions with a number of integration steps equal to the number
of elements shows in the two figures below that the step-by-step integration is a little
closer to the exact solution when the number of elements is small and both

32
Machine Translated by Google

MEF course – Chapter 3: Differential equations A. Seghir 2005-2014

solutions coincide almost perfectly with the exact solution when the number of elements is
large.
The difference in the case of a small number of elements is due to the fact that the element
used is linear, therefore the interpolation polynomial is of the first order. On the other hand,
step-by-step integration uses an approximation with a cubic polynomial since the method is
linearly accelerated.

10 elements and 10 integration steps 20 elements and 20 integration steps

3.11 Duty No. 3

1) Modify the program EquaDiff2(n) to obtain EquadiffMEF(n). Calculate the elementary

matrices and vectors with MATLAB and change the function MatElt2Nd(x1,x2) to
calculate the new matrices and the new
vector.

2) Introduce the SdofTH function into MATLAB and retrace the two figures
previous
3) Use the function FF1D(n) below to calculate the shape functions of a quadratic
element (element with three nodes) and calculate the expressions of

33
Machine Translated by Google

MEF course – Chapter 3: Differential equations A. Seghir 2005-2014

2
of
elementary matrices and vectors of the previous equation: 0ÿuÿxÿ with
dx2
the same CALs.
function [N,dN] = FF1D(n)
%------------------------------------------------- ------------------
-

% [N,dN] = FFDim1(n)

% N: Jacobi-type one-dimensional element shape function

% dN: derivatives of shape functions
% n: number of nodes of the element

% the function returns the symbolic expression of the functions of

shape

% A. Seghir, 12/09/04
%------------------------------------------------- ------------------
-

x = sym('x','real'); % character x

for i=1:n % loop to add the

characters 1 2 3 ...

c = char(48+i); % attention to n > 9

xn(i) = sym(['x',c],'real');
end;

for i=1:n % loop for N(i)

N(i) = sym(1); for j=1:n % initialization to 1

% loop for product

if j~=i

N(i) = N(i)*(x-xn(j))/(xn(i)-xn(j));
end
end
end

N = expand(N); N = simple(N); % rearrangement of the terms of

NOT

dN = simple(diff(N,x)); % calculation of derivatives

return

34
Machine Translated by Google

MEF course – Chapter 3: Differential equations A. Seghir 2005-2014

3.12 General program for the equations of 2nd

order
The structure of the MATLAB programs produced so far does not allow all second-order differential
equations to be processed. More or less slight modifications are necessary for each case, they mainly
affect the elementary matrices
and boundary conditions.

It is possible to think about putting the function for calculating elementary matrices in another file and
creating another function for processing boundary conditions to keep the rest of the program in a function
which will not be modified. But, in addition to the fact that the type of element remains fixed, the expressions
of the elementary matrices themselves always pose a problem since they must be calculated by integration
over the element.
However, in certain cases, integration can be difficult or almost impossible; digital integration is a solution
that is perfectly suitable, particularly if it is associated with an element with fixed terminals called a parent
element or reference element.

Furthermore, the coefficients of the equation can also be functions of the variable x. To implement a
program that treats very diverse equations without modification, we must first characterize the variables
which are data of the problem.

A boundary value problem defined in a domain ÿ ÿ [x0, xÿÿ and governed by a 2nd order linear differential
equation is written as follows:

2
of of
2
ÿ cx ( ) ÿ dx kxufx
() ÿ
() , with u(x 0) ÿ u 0 And fromx(ÿ)ofÿ
ÿ
dx dx

If the second order term of the equation is multiplied by a function mÿxÿ then it is better to divide the whole
equation by this term otherwise the variational formulation
weak will involve one more integral term containing the derivative of m(x) with respect to
x.

The parameters of the equation are the functions k, c and f as well as the values u0 and duÿ, those of the

domain ÿ (of geometry) are x0 and xÿ and those specific to the method are the type of the element
characterized by the number of nodes it contains nne and the number

35
Machine Translated by Google

MEF course – Chapter 3: Differential equations A. Seghir 2005-2014

of net total elements . It is best to group all these parameters into a data structure
called Eq like Equation.
The finite element program that we want to obtain can be called with data preparation
in the following way:

function TstMefEqDiff;
%--------------------------------
% Test of the MefEquaDiff program
%

% A. Seghir, 14/10/04
%--------------------------------

Eq.c = @coef_c; the % the coefficients of

equation point
Eq.k = @coef_k ; % on function names
implemented Eq.f
= @func_f ; % lower

Eq.nne = 3; NetEq = 20; % nbr nodes per elt and nbr.

elt. total

Eq.xi = 0; Eq.xf = pi; % integration terminals and CAL

Eq.u0 = 1; Eq.Du = 1;

[U,x] = MefEquaDiff(Eq); % call to program

t = Eq.xi:0.01:Eq.xf; Ue = % exact solution

1/3*cos(t)+2/3*cos(2*t)+1/2*sin(2*t);

plot(x,U,'.',t,Ue); % plot of both solutions

return

function c = coef_c(x); c = 0; return; % declaration of coefficients function k = coef_k(x); k = 4; return; % of

equation
function f = func_f(x); f = cos(x); return;

We can clearly see that with this type of script, it is easy to process very diverse equations; it suffices to give the

expression of the coefficients in the functions coef_k, coef_c and func_f, the limits of the domain as well as the

boundary conditions can be

36
Machine Translated by Google

MEF course – Chapter 3: Differential equations A. Seghir 2005-2014

changed using the values of Eq.xi, Eq.xf, Eq.u0 and Eq.Du. In addition we now have the choice of the type of elements

and their number with the variables Eq.nne and

Eq.net.

Noticed
The exact solution of the differential equation treated as an example can be obtained
with these commands:
>>ue = dsolve('D2u + 4* u = cos(t)','u(0)=1','Du(pi)=1');
>>ue = simple(ue)
eu =
1/3*cos(t)+2/3*cos(2*t)+1/2*sin(2*t)

To achieve such flexibility, the program must be fairly comprehensive. Its general structure
is:
Start of the program
Geometry (coordinates and connectivity)
Initialization of matrices
Start of loop on elements

Function call to calculate ÿ dNT dN dx

Function call to calculate ÿ NT c(x) dN dx
Function call to calculate ÿ NT k(x) N dx
Function call to calculate ÿ NT f(x) dx
Assembling the K matrix
Assembling vector F
End of loop on elements

Application of boundary conditions

Resolution
End of program

The MefEquaDiff function file is as follows:

function [U,p] = MefEquaDiff(Eq)

%------------------------------------------------- ------------
% [U,p] = MefEquaDiff(Eq)
% U: vector containing the values of the solution at the nodes

% p: node coordinates
% Eq: structure containing the parameters of the equation

37
Machine Translated by Google

MEF course – Chapter 3: Differential equations A. Seghir 2005-2014

% A. Seghir, 16/10/04
% ------------------------------------------------- -----------

[t,p] = mesh1d(Eq.net,Eq.nne,Eq.xi,Eq.xf); % connectivity and

contact details
nnt = length(p); % total number of nodes

K = zeros(nnt); % matrix initialization

K
F = zeros(nnt,1); % initialization of vector F

for ie = 1:Eq.net % loop over elements

te = t(ie,:); % element connectivity

current

X = p(te)'; % element coordinates

current

Me = -intdNTdN(X); % calculation of the integral ÿ dNT

dN dx

Ce = intNTcxdN(Eq.c,X); % calculation of the integral ÿ NT c

dN dx

Ke = intNTkxN (Eq.k,X); % calculation of the integral ÿ NT

k(x) N dx

Fe = intNTfx (Eq.f,X); % calculation of the integral ÿ NT

f(x)dx

K(te,te) = K(te,te) + Me + Ce + Ke; % die assembly

K
F(te) = F(te) + Fe; % assembly of vector F
end; % end of loop

F = F -K(:,1)*Eq.u0; % application of CALs

F(nnt) = F(nnt) - Eq.Du;
K(1,:) = [];
K(:,1) = [];
F(1) = [];
U = K\F; % resolution

U = [Eq.u0;U];
return

The first thing we notice in this version of the program is that the domain is meshed
with the mesh1d function which allows us to calculate both the coordinates of the
nodes and the connectivities of the elements. This way of doing things is particularly
useful for cases of two or three dimensions, in particular when we want to use a file of

38
Machine Translated by Google

MEF course – Chapter 3: Differential equations A. Seghir 2005-2014

mesh result of an automatic meshing program. For this simple case, we propose this
function as follows:
function [t,p] = mesh1d(net,nne,xi,xf)
%------------------------------------------------- ----

% t,p] = mesh1d(net,nne)
%t : table of connectivity of line elements
%p : table of node coordinates
% net: total number of elements
% nne: number of nodes per element
% xi,xf: limits of the integration domain
%

% A. Seghir, 14/10/04
%------------------------------------------------- ----
t(1,:) = 1:nne;
for i = 2:net
t(i,:) = t(i-1,:)+nne-1;
end;
n = (nne-1)*(net);
dx = (xf-xi)/n;
p = xi:dx:xf;
return

The second remark is more important, it concerns the calculation of elementary

matrices. The MatElt2N function used in previous versions is now replaced by four
functions, three of which: intdNTdN, intNTcxdN, intNTkxN return the elementary
matrices after integration and the fourth intNTfx calculates the elementary force
vector. The input arguments are a vector _ _ _ _ _ _ ÿxÿ. Since the size of the vector
X is not specified, the functions support any type of line element.

It should be remembered that it is in these functions that the integration of matrices

takes place, and this is precisely what limited previous versions, so it is here that we
expect modifications. Indeed, the exact integration of the components of the matrices
is replaced by a weighted sum of the values of these components at specific points
according to the Gaussian numerical integration method. This method, which will be
explained later, allows us to evaluate the integral of a function using a sum, i.e.:

39
Machine Translated by Google

MEF course – Chapter 3: Differential equations A. Seghir 2005-2014

b not

ÿ fx(dx
) wfx i ( )i
has
ÿÿ ÿi 1

With f any function of x, n is the number of points xi in which the function is evaluated and
wi are weighting values. The coordinates xi as well as the weights wi
are calculated and the precision of this numerical integration depends on the number of
points n, the result is exact if the function f is a polynomial of order 2n. We give below the
points and weights for the integration on the interval ÿÿ1,ÿ1ÿ:

n=3
x -0.34641016151378 0 0.34641016151378
w 0.55555555555556 0.88888888888889 0.55555555555556
n=4
x -0.861136311594053 0.861136311594053 -0.339981043584856
0.339981043584856
w 0.347854845137454 0.347854845137454 0.652145154862546
0.652145154862546
n=6
x -0.932469514203152 0.932469514203152 -0.661209386466265
0.661209386466265
w 0.171324492379170 0.171324492379170 0.360761573048139
0.360761573048139
x -0.238619186083197 0.238619186083197
w 0.467913934572691 0.467913934572691

The MATLAB function implemented for the program is based on a number of points n ÿ 6,
which allows us to use elements with up to 11 nodes practically without loss of precision.

function [xsi,w] = wpgL6p

%------------------------------------------------- ------------------
---

% [xsi,w] = wpgL6p
% xsi, w: points and Gaussian weights between -1 and +1 in line vectors

% number of points: 6 points

% A.Seghir, 03/08/04

40
Machine Translated by Google

MEF course – Chapter 3: Differential equations A. Seghir 2005-2014

%------------------------------------------------- ------------------
---

xsi = [ -0.932469514203152 , 0.932469514203152 , ...

-0.661209386466265 , 0.661209386466265 , ...
-0.238619186083197 , 0.238619186083197 , ...
];

w = [0.171324492379170 , 0.171324492379170 , ...

0.360761573048139 , 0.360761573048139 , ...
0.467913934572691 , 0.467913934572691 , ...
];
return

With this integration method, the exact expressions of the elementary matrices are no longer necessary,
especially since we do not know, a priori, the functions to be integrated; the coefficients of the equation
can be arbitrary.

The last remaining point consists of making the transition between the integration terminals of
the real element and the terminals ÿ1 and ÿ1 of the reference element. The transformation

geometric :

xÿNÿ

with ÿ the coordinate in the reference element, offers the best solution because it uses the N shape
functions also defined in the reference element. We can easily verify that if ÿ ÿ ÿ 1, x ÿ x1, x2 for a linear
element and x ÿ x1, x3 for a quadratic element, … etc.

With this geometric transformation, the derivatives of the shape functions in the two elements are linked by
the relation:

dN dN dx
ÿ
;
dÿ dx d ÿ

dN ÿ J d dN
Or in more general form ÿ dx with J becomes a matrix (Jacobian) in the

case of two or three dimensions: J ÿ ÿ dx

N ()x ÿ ddN x d .
ÿ dÿ not not

1
Example of a linear element: NOT ÿ ÿ1ÿ
2
ÿ, 1ÿ ÿ ÿ ; 1dN ÿ_ ÿ2ÿ1 , 1ÿ ; And

Jÿxÿx1_( _ _ _
2 2 1)

41
Machine Translated by Google

MEF course – Chapter 3: Differential equations A. Seghir 2005-2014

The limits of the integrals are thus transformed = ÿ1:

x2 ÿ 1
ÿ ÿ
ÿ dx f NJ d
( ) ÿ ( ( )) fx ÿ
x1 ÿ

Finally, with these two major modifications which concern the geometric transformation and
numerical integration, the four functions of the elementary matrices and vectors will have
the following structure:

Start of function

Gaussian integration points and weights

Initialization of the elementary matrix (vector)
Start of loop on Gauss points
Evaluation of shape functions at the Gaussian point
Calculation of J = dN * Xn

Calculation of dN/dx = Jÿÿ dN

Sum of products dNT.dN, dNT.N, NT.N with multiplication by
weights and J
End of loop on Gauss points
End of function

The source code for the four functions is as follows:

function Me = intdNTdN(Xn)
%----------------------------------

% integration of: dN' * dN

%----------------------------------

[xsi, w] = wpgL4p;
npg = length(xsi);
nne = length(Xn);
Me = zeros(nne);
for i =1:npg
[N, dN] = FF1D(xsi(i), nne);
J = dN * Xn;
dN = dN * inv(J);
Me = Me + det(J) * w(i) * (dN' * dN);
end
return

42
Machine Translated by Google

MEF course – Chapter 3: Differential equations A. Seghir 2005-2014

function Ce = intNTcxdN(cx,Xn)
%--------------------------------------

% integration of: N' * c(x) * dN

%--------------------------------------

[xsi, w] = wpgL4p;
npg = length(xsi);
nne = length(Xn);
This = zeros(nne);
for i =1:npg
[N, dN] = FF1D(xsi(i), nne);
x=N* Xn;
J = dN * Xn;
dN = dN * inv(J);
c = feval(cx,x);
Ce = Ce + det(J) * w(i) * vs * (N' * dN);
end
return

function Ke = intNTkxN(kx,Xn)
%----------------------------------

% integration of: N' * k(x) * N

%----------------------------------

[xsi, w] = wpgL4p;
npg = length(xsi);
nne = length(Xn);
Ke = zeros(nne);
for i =1:npg
[N, dN] = FF1D(xsi(i), nne);
x=N* Xn;
J = dN * Xn;
k = feval(kx,x);
Ke = Ke + det(J) * w(i) * k * (N' * N);
end
return

function Fe = intNTfx(fx,Xn)
%------------------------------

% integration of: N' * f(x)

%------------------------------

[xsi, w] = wpgL4p;
npg = length(xsi);

43
Machine Translated by Google

MEF course – Chapter 3: Differential equations A. Seghir 2005-2014

nne = length(Xn);
Fe = zeros(nne,1);

for i =1:npg
[N, dN] = FF1D(xsi(i), nne);
x=N* Xn;
J = dN * Xn;

f = feval(fx,x);
Fe = Fe + det(J) * w(i) * f * N';
End
return

Shape functions as well as their derivatives for reference elements of order n

are :
function [N,dN] = FF1D(x,n)
%------------------------------------------------

% [N,dN] = FF1D(x,n)

% N: shape functions of 1D line elements

% dN: derivatives of shape functions

% x: function evaluation point

% n: number of nodes of the element
%

%HAS. Seghir 14/10/04

%------------------------------------------------

if isa(x,'sym') symbolic % determines if x is

one = sym(1);
zero = sym(0);
else

one = 1.0;
zero = 0.0;

if (x < -1) | (x > 1); error('wrong x coordinate'); end;

end;

xn =-1:2/(n-1):1;

for i=1:n

p(i) = one;
q(i) = one;
s(i) = zero;

44
Machine Translated by Google

MEF course – Chapter 3: Differential equations A. Seghir 2005-2014

for j=1:n if j~=i

p(i) = p(i) * (x - xn(j));

q(i) = q(i) * (xn(i) - xn(j));
r = one;
for k =1:n

if (k~=i & k~=j); r = r*(x - xn(k)); end;

end;
s(i) = s(i) + r;
end;
end;
end;
N = p./q;
dN = s./q;
return

3.13 Duty No. 4

2 2
1) Repeat the equation from homework No. 3, question 3: du dx ÿ u ÿ x ÿ 0 with this

program. Modify the functions coef_k, coef_c and func_f, as well as the domain limits and boundary conditions.

Carry out a study of convergence in type and number of elements.

2) Solve the equation with the MefEquaDiff program:

2 2
of 1 of 1ÿ x
ÿ ÿ

2
dx x dx x

of
with the boundary conditions: u(0) ÿ1 And (4ÿ) ÿ ÿ2
dx

digital integration
function TstIntNum
clear
clc
xn =[-2 , 2];
[xsi,w, npg] = wpgL6p;

45
Machine Translated by Google

MEF course – Chapter 3: Differential equations A. Seghir 2005-2014

sn = 0;
for i=1:npg
[N,dN] = FF2Nd(xsi(i));
x = xn*N;
dx = xn*dN;
f = func(x);
sn = sn + dx* w(i) * f;
end;
x = sym('x','real'); se =
int(func(x),xn(1),xn(2))
int(func(x))
sn
se = double(se)
err = sn-se

function f = func(x)
f = x.*exp(-x);
return

function [N,dN] = FF2Nd(xsi)

N = 0.5*[1-xsi; 1+xsi];
dN = 0.5*[-1; 1];
return;

46
Machine Translated by Google

chapter 4

Bar Element

4.1 Governing equation

The bar element is used in assemblies of bars or rods working in tension or
compression. They are mainly found in metal frames and in trellis systems.

To formulate this element, we consider a bar of section A of length L subjected to a

traction P(x) varying from P0 to PL.

ÿ
P0 PL P P+ Pdx
ÿx _

x
dx

Fig. 2.1: Elementary balance

An infinitesimal portion of length dx located at the coordinate x along the bar is in

dynamic equilibrium under the following system of forces:

ÿF ÿ mÿ

ÿ ÿ Pÿ ÿ 2u
P ÿ ÿ ÿ ÿ( 2dxÿ PA
ÿ t dx ) (2.1)
ÿÿ
x ÿÿ
Machine Translated by Google

MEF Course – Chapter 4: Bar Element A. Seghir 2005-2014

In this expression u designates the longitudinal displacement, x is the coordinate and t is the time.

If we designate by E the modulus of elasticity of the material with which the bar is made, Hooke's law
gives the axial stress as a function of the longitudinal deformation:

P
ÿ ÿ ÿxÿ x E (2.2)
HAS

The deformation is linked to the displacement by the derivative with respect to x :

ÿu
ÿÿx (2.3)
ÿx

By substituting (2.3) into (2.2) we obtain the expression of the traction P as a function of the displacement
u:

ÿu
PEA ÿ ÿ
(2.4)
x

Hence the differential equation of the balance of the element dx which is obtained by replacing (2.4) in
(2.1):

ÿÿu( 2u
EA ÿ ÿ) ÿ A (2.5)
ÿx ÿx 2
ÿt

The boundary conditions of this equation depend on the embedding or displacements imposed for u,
and on the loading at the nodes for the derivatives of u (Equation 2.4).

4.2 Element wording

Taking ÿu as a weighting function, the strong variational formulation of equation (2.5) is written:

L ÿ ÿu 2u

ÿ 0ÿ
u ( EA ÿ )] ÿ ÿu ÿA 2 ÿ t dx ÿ
0 (2.6)
[ ÿx ÿx

The weak formulation is written by taking the integration by parts of the first term:

48
Machine Translated by Google

MEF Course – Chapter 4: Bar Element A. Seghir 2005-2014

L
L ÿu ÿu L 2u ÿ ÿ dx ÿ ÿ u u
EA ÿ dx ÿÿÿ ÿ u2Aÿ t AE ÿ ÿ 0 (2.7)
ÿ 0ÿ ÿx ÿx 0 ÿx ÿ ÿ ÿ ÿ
0

The last term is only the difference in the forces applied to the ends of the bar: ÿu|
xÿL PL ÿ ÿu|xÿ0 P0.
For the discretization of equation (2.7) take linear form functions with x1 ÿ 0 and x2 ÿ L.
The expressions of these functions and their derivatives are:

x L ÿÿ ÿ

x ÿ

01
N x( ) 0 ÿ0 L xx ÿÿ ÿ ÿ ; 1 dN
x ( ) ÿ ÿ ÿ1 1ÿ (2.8)
ÿ

L L ÿ

L L

So by replacing for the different operators we obtain:

ÿ ÿL (x)/W ÿ ÿ ÿ ÿ ÿÿ 1/ L ÿ
ÿ uÿÿ ÿ u0 ÿ uL ÿ ÿ
ÿ uu ÿÿ ÿ u0 ÿ uÿ Lÿ ÿ ÿ
(2.9a)
/xL ÿ
ÿx 1/ L ÿ

ÿÿ ÿ
ÿ uÿÿ ÿ u0 1 ÿÿ ÿ ÿ ÿ uL ÿ 0 ÿ u0 ÿ ÿÿ ÿ ÿ ÿ ÿ ÿ u0 uL ÿ ÿ (2.9b)
x 0ÿ x Lÿ
1

ÿÿ

ÿu ÿ u0 ÿ ÿ 2u 1 ÿ u0 ÿ

LL ÿÿ ÿ ÿ
1/ 1/ ÿ ÿ 2
tL L xx ÿ ÿ ÿ ÿ ÿ ÿ ÿÿ
ÿ (2.9c)
ÿx ÿ uL ÿ ÿ uL ÿ

Or :
ÿÿ

u represents the second derivative with respect to the travel time; that's to say
acceleration.

With these expressions the weak integral equation (2.7) becomes after simplification of
ÿ ÿu0 ÿuL ÿ :

L ÿÿ 1/ ÿ u0
ÿ ÿ
EA ÿ ÿ ÿ ÿÿ 1/LL
1/ dx ÿ ÿ
ÿ
0 1/ L
Lÿÿ _
u
ÿ L
ÿÿ (2.10)
L 1 ÿLÿxÿ ÿ 1 ÿLÿxÿ ÿ ÿ u0 ÿ ÿÿ P0ÿ
ÿ ÿ ÿA ÿ
dx ÿ ÿÿ
ÿ
ÿ
ÿ
0 L ÿ x ÿ
L ÿ x u P
ÿ ÿ L ÿ ÿÿ ÿ ÿ L ÿ

ÿÿ

or in matrix form: Ke Ue ÿ MeUe ÿ Fe , with :

49
Machine Translated by Google

MEF Course – Chapter 4: Bar Element A. Seghir 2005-2014

L ÿÿ 1/ ÿ
K
e
ÿ
ÿ ÿ
EA ÿ ÿ ÿ 1/ 1/
LL dx
0
Lÿÿ _
1/ L ÿ

(2.12)
1 L ÿ 1 ÿ

1 EA ÿ
1 ÿ

1ÿ
ÿ
ÿ EA ÿ dx
ÿÿÿ
L2
ÿ
0
ÿÿ
ÿ

1 1 Lÿ ÿ

1 1ÿ ÿ

L 1 ÿLÿxÿ ÿ1ÿ
M ÿ
ÿ HAS L xxdx ÿ ÿ ÿ ÿ ÿ
e
0 L
ÿÿ
x L
(2.13)
Lx 21
2
ÿA L () Lxÿ
ÿ

ÿ ÿ ÿ ÿ ÿ ÿ ÿ ÿ ( ) xA.L.
ÿ

ÿ
ÿ dx
L2 0
ÿ (Lxx )
ÿ

2x _ 6 ÿÿ12ÿ

ÿ ÿ 0ÿ ÿ ÿ ÿÿ1ÿÿÿÿ ÿÿ 0
ÿ
F ÿ
P ÿ
ÿÿ
P0 ÿ (2.14)
e L
1 0 Pÿÿ _
P
L ÿ

The elementary matrices Ke and Me are called respectively rigidity matrix and mass matrix. Fe is the
loading vector of the element. These expressions are valid if the section A, the density ÿ and the Young's
modulus E do not vary along the entire bar.
In the case where these quantities vary, it is possible to subdivide the bar into several elements and take
average values for each element.

The expression for the mass matrix as obtained in (2.13) is called coherent or distributed mass. It is
possible to concentrate the mass of the element at its ends.
Half of the total mass of the element is attributed to each of the two nodes, i.e.:

Me vs
ÿ
ÿ ÿ10
ÿAL
(2.15)
2 ÿÿÿ01ÿ

4.3 Example of a recessed rod

A rod of length L ÿ 3m and section A ÿ 25cmÿ is embedded at one of its ends
and subjected to a tension P ÿ 250KN at the other end. The density and modulus of

Young of the rod material are: ÿ ÿ7800 Kg/ mÿ and E ÿ 210,000 MPa.

250KN _

Fig. 2.2: Recessed bar

50
Machine Translated by Google

MEF Course – Chapter 4: Bar Element A. Seghir 2005-2014

The mass and rigidity matrices as well as the loading vector are:

ÿ 1 ÿ

1ÿ ÿ1 0 ÿ ÿ 0 ÿ
M 29.25
ÿ N/ m ;
ÿ
K ÿÿ
175 106
ÿ
kg ; Fe ÿ
4 ÿ
NOT
ÿ
ÿ0 1
ÿ ÿ

25 10
1 1ÿ
ÿ

ÿ
ÿ

ÿ ÿ ÿ

Taking into account the condition u(0) ÿ 0. The resolution gives a displacement u(L) ÿ
1.4mm

and a period T ÿ0.0026 sec.

4.4 Plane lattice structures

Lattice structures are made up of assemblies of bars linked by joints in such a way that the external
loading is taken up only by axial forces in the bars. Figure 2.3 shows an example of a truss system
composed of an assembly of 13 bars and subjected to a loading of two forces.

F1 F2

Fig. 2.3: Lattice system

The bars making up this system undergo two movements at their ends; a horizontal component and
another vertical. However, only the axial displacement of the bar gives rise to the axial force. Thus,
the elementary matrices of a bar

51
Machine Translated by Google

MEF Course – Chapter 4: Bar Element A. Seghir 2005-2014

two-dimensional (fig.2.4) become 4ÿ4 matrices since the vector of elementary

displacements becomes:
T
Ue ÿ ÿ u1 v1 u2 v2 ÿ (2.16)

v1 v2

u1 EA u2

Fig. 2.4: Two-dimensional bar element

All the components of the stiffness matrix associated with the component v of the
displacement are zero. Those of the mass matrix, on the other hand, are not all. The
effects of inertia also exist when the ends of the bar move in the orthogonal direction,
only there is no coupling between the two components u and v of the displacement.
The two matrices are written as follows:

ÿ1010ÿÿÿ1
ÿ

ÿÿ
ÿ ÿÿ ÿ2 10 01 20 0ÿ ÿÿ ÿÿ ÿ
0000 ÿÿ
AE
ÿ AL 0201
Ke ÿ
; Me ÿ
(2.17)
L ÿ

010ÿÿ 6

ÿÿÿÿ
0000 0102

The concentrated mass matrix is written in the same way as in (2.15) by assigning half
of the mass of the element for all the components of the displacement:

ÿÿ ÿ1000ÿ
0100 ÿ
ÿ AL
Me ÿ
(2.18)
2 ÿ0010ÿÿ

ÿÿÿÿ
0001

52
Machine Translated by Google

MEF Course – Chapter 4: Bar Element A. Seghir 2005-2014

The expressions of the elementary matrices as they are written by equations (2.17) and (2.18) are
valid in a system of axes coincident with the longitudinal axis of the bar element. In practical cases,
the bar can be inclined relative to the reference axes. It is then necessary to rotate the axes to
return to the axes of the bar. v2
uÿ2
vÿ2

u2

y
yÿ ÿ

xÿ
v1
vÿ1 uÿ1

u1
oÿ x
o

Fig. 2.5: Inclined bar element

Figure 2.5 shows a bar inclined at an angle ÿ with respect to the horizontal axis of the coordinate
system (oxy). We note (u,v) the components of the displacement in this system and (uÿ,vÿ) those
of the displacement in the reference frame (oÿxÿyÿ) linked to the bar. We also note R the reference
rotation matrix which makes it possible to make the transition from the system (oxy) to the new
system (oÿxÿyÿ). We can thus write the following relationships:

ÿ
ÿ
u1ÿ ÿ cs ÿÿÿÿ u1ÿ ÿ ÿu2ÿ ÿ cs ÿÿÿÿ u2 ÿ

ÿ ÿ
ÿ
ÿ ÿ
; ÿ ÿ ÿ
ÿ
ÿ ÿ
ÿ
v1 ÿ

sc v1 ÿ
v2
ÿ ÿÿ
ÿ

sc ÿ
v2
ÿ
ÿ ÿ ÿÿ ÿ ÿ

with c ÿ cos(ÿ) and s ÿ sin(ÿ).

Or with the totality of the vectors of elementary displacements Ue and Uÿe :

53
Machine Translated by Google

MEF Course – Chapter 4: Bar Element A. Seghir 2005-2014

ÿ
ÿ
u1ÿ ÿÿ
cs 0ÿÿÿ0 ÿÿ0 u1 ÿ

0ÿÿ
ÿ
v1
ÿ

sc ÿÿÿÿÿÿ v1
ÿ ÿ
ÿ ÿÿ ÿÿ ; Eu ÿ RUe (2.19)
ÿ
u2 00 cs u2

ÿ
ÿÿÿÿÿ v 2 ÿÿÿÿ ÿÿÿÿ
00 ÿ

sc v2
ÿÿÿÿ

It now remains to link the loading to the displacements in the frame (oxy), we have the
Fe Ke Ue . The relationship between the two vectors
ÿ ÿ ÿ
relationship in the frame linked to the ÿ

bar forces in the two frames is:

, or : Fe ÿ RT Fÿe ; with: Fe ÿ ÿ Fx1 Fy1 Fx2 Fy2 ÿ

ÿ
T
Fe ÿ R Fe (2.20)

from which it results:

Fe ÿ RT KÿeeUÿ ÿ R Kÿ
T
RU
e e ÿ Ke Ue
(2.21)

We thus obtain the relationship between the expressions of the rigidity matrix in the two
benchmarks:

T ÿ

KReKe R
ÿ
(2.22)

Either :

cs ÿ

ÿ00ÿÿ ÿÿ ÿ1010ÿÿÿ ÿ
ÿ

cs 00ÿÿÿ ÿ

ÿ ÿ ÿ
sc 00 AE 0000 ÿ

sc 00
Ke ÿ

ÿÿÿÿ00 vs
ÿ

s L ÿ

ÿÿ1010ÿÿ 00 cs
ÿÿ
ÿÿ
00 sc ÿÿ ÿÿÿÿ
0000 00 ÿ

sc ÿ

2
ÿ tbsp 2 tbsp c ÿÿ

cs ÿÿ

ÿ 2 2
AE ÿ ÿ cs sec
ÿÿ

cs sec
ÿ ÿ
2
L ÿ

tbsp
ÿ

2 tbsp c cs

2 2
ÿ

cs ÿ

sec cs sec
ÿÿ ÿÿÿÿÿ

the stiffness matrix, the coherent or concentrated mass matrix remains unchanged at the
reference rotation. Indeed, from the point of view of physics, the mass is independent of the
orientation of the reference frame.

54
Machine Translated by Google

MEF Course – Chapter 4: Bar Element A. Seghir 2005-2014

The MATLAB script which allows verifications by calculation of elementary matrices

is:

%--------------------------------------------
% Elementary matrices of the bar element
%

% A. Seghir, 08/10/04
%--------------------------------------------
clear % clears all variables
in global memory
clc % clears the window
orders
syms sc real symbolic % declares variables

R = [ cs 0 0 % rotation matrix
-sc 0 0
0 0 cs
0 0 -sc
]

K = [1 0 -1 0 % stiffness matrix in the

bar mark
0000 % the A*E/L factor is
intentionally omitted

-1 0 1 0
0000
]

M=[20 10 % mass matrix in the reference

linked to the bar
0 2 01 % ale factor rho*A*L/6 is
also omitted
1 0 20
0 1 0 2
]

Ke = R' * K * R % stiffness matrix in the

global benchmark (oxy)
Me = R' * M * R % mass matrix in the reference
global (oxy)

55
Machine Translated by Google

MEF Course – Chapter 4: Bar Element A. Seghir 2005-2014

4.5 Example of two bars

Figure 2.6 shows a structure composed of two bars embedded at their upper ends and
assembled at their lower ends with a ball joint from which a weight P ÿ 3400 KN is suspended.
y

The bars have the same characteristics: 2

1ÿ ÿ

Length: L ÿ 1 m; section: A ÿ 25 cmÿ

1 2
3 x
and Young’s modulus E ÿ 210,000 MPa.

The embedding angle is: ÿ ÿ 36.87°. P

Which gives for bar No. 1: c ÿ ÿ 0.8 and s ÿ 0.6 Fig. 2.5: Two-bar trellis

And the first elementary system: Ke1 Ue1 ÿ Fe1

1 e1
ÿÿ
0.6 4 0.4 8
ÿ ÿ

ÿ ÿ 0.6 4 0.4 8 ÿ ÿ e ÿ ÿF ÿ ÿ

u1 1x _

ÿÿ e1 1st _

6
ÿ

0.4 8 0.3 6 0.4 8 0.3 6 ÿ v 1

ÿ

F1
y1
525 1 0
ÿ

ÿ
ÿ
ÿÿ
1
ÿÿ
e ÿÿ
ÿ

0.6 4 0.4 8 0.6 4 0.4 8 ÿ ÿ ÿ ÿ

ÿ

and Fÿ
ÿ 2_ 2x _

e1 1st _

0.4 8 0.3 6 0.4 8

ÿ ÿ

0.3 6 v2 F2
y
ÿÿÿÿ ÿÿÿÿ ÿÿ ÿÿÿÿ

By involving the displacement vector of all the nodes (from 1 to 3) the system becomes:

ÿÿÿ
0.6 4 ÿ

0.4 8 0 0 0.6 4 0.4 8 ÿ

u1 ÿÿ
ÿ F 1st _
ÿÿ

x1
ÿ ÿ
e1
ÿ

0.4 8 ÿ ÿ ÿ 0.3 6 0 0 0.4 8 0.3 6 ÿ ÿ ÿ ÿ ÿ ÿ

v1 Fÿÿÿÿ

ÿ 1y

0 000 0 0 u2 0
6
525 1 0
ÿ

ÿÿ
ÿ
ÿ ÿÿ

0 000 0 v2 0
0ÿ
ÿ

0.6 4 0.4 8 0 0 0.6 4 0.4 8 ÿ

u3 F 1st _

x2

0.4 8 ÿ

ÿ ÿ ÿ ÿ ÿ 0.3 6 0 0 0.4 8 0.3 6 ÿ ÿ

ÿ

v3 F e1
ÿÿÿÿÿÿ ÿÿÿÿÿ
ÿ 2y
ÿÿÿ
ÿÿÿÿÿÿÿ

For bar No. 2: c ÿ 0.8 and s ÿ 0.6, the system with all the components of the displacement
vector is:

56
Machine Translated by Google

MEF Course – Chapter 4: Bar Element A. Seghir 2005-2014

ÿ00ÿ 0 0 0 ÿÿÿÿ0
u1 ÿÿ ÿÿ
0 ÿÿ

0ÿ
00ÿÿ 0 0 0 ÿÿ v1 0
ÿ 0.4
e2
0 0 0.6 4 0.4 8 ÿ ÿ ÿ 0 0 0.6 ÿ

0.6 4 ÿ

ÿÿÿÿÿÿ 8 u2 F 1x _
6
525 1 0
ÿ

0.4 8 ÿÿ
ÿ
ÿÿ e2 ÿÿ

0 0 0.4 8 0.3 6 0.4 8 ÿ ÿ

0.3 6 v2 F
ÿ 1y
e2
4 0.4 8 ÿ ÿ ÿ ÿ

0.6 4 ÿÿÿ u3 ÿÿÿÿ Fÿ

2x _

e2
0 0 0.4 8 0.3 6 0.4 8
ÿ ÿ

0.3 6 v3 F
ÿÿÿÿÿ ÿÿ y 2 ÿÿÿÿÿÿÿ

The assembly of the two elementary systems gives the overall system:

ÿÿÿ
0.6 4 0.4 8 ÿ

0 0 ÿ

0.6 4 0.4 8 ÿ ÿ 0.3 ÿ ÿ u 1 ÿÿ ÿÿ

0 ÿÿ

6
ÿ

0.4 8 0.3 6 0 0 0.4 8 ÿ

ÿÿÿ v1 0
ÿ
0 0 0.6 4 0.4 8 ÿ

0.6 4 0.4 8 ÿ ÿ 0.3 ÿ

u2 0
6
525 1 0
ÿ

6 ÿÿ
ÿ
ÿÿ ÿÿ

0 0 0.4 8 0.3 6 0.4 8 ÿ ÿ

ÿÿÿ v2 0
ÿ0
ÿ

0.6 4 0.4 8 ÿ

0.6 4 0.4 8 ÿ

1.2 8 ÿ

ÿ u3 0
ÿ 6
ÿÿÿÿÿÿ
0.4 8 ÿ

0.3 6 0.4 8 ÿ ÿ

0.3 6 0 0.7 2 ÿ ÿ v3
ÿÿÿÿÿ ÿÿÿÿÿ
ÿÿ
3.4 1 0 ÿÿÿÿÿ

The application of boundary conditions (embedding: u1 ÿ v1 ÿ u2 ÿ v2 ÿ 0) reduces this

system to two unknowns:

ÿ 1.28 0 u3 ÿ ÿ 0 ÿ
6
525 10
ÿ

ÿ ÿ
ÿ
ÿ ÿ
6
ÿÿ
ÿ ÿ ÿ 0 0.72 ÿ ÿ v 3 ÿ ÿ
ÿÿ
3.4 10 ÿ

The solution gives: u3 ÿ 0; v3 ÿ ÿ 8.99 mm which is a downward shift.

To calculate the reactions, it is sufficient to perform the matrix products of the elementary
systems with the displacements now known and the internal forces unknown; we obtain
for bar N°1:
e1
ÿÿ
0.6 4 0.4 8 ÿ ÿ

0.6 4 0.4 8 ÿ ÿ u1 ÿ
0 ÿ
ÿ
ÿ

2265.4 8
ÿ 1
6
ÿ

0.4 8 0.3 6 0.4 8 ÿ

ÿ 0.3 6 ÿ ÿ ÿ e ÿ
0 ÿ ÿ ÿ 1699.1
v1
525 1 0
ÿ ÿ ÿ KN ÿ
ÿ e1 ÿÿ ÿ 1ÿ
ÿ

0.6 4 0.4 8 0.6 4 0.4 8 ÿ

u2 ÿ
0 2265.4 8
ÿ
1 ÿ

3
0.4 8 ÿ

0.3 6 0.4 8 ÿ

0.3 6 ÿ ÿ ÿ e ÿÿÿ
8.9 9 1 0
ÿ

ÿ 1699.1 1 ÿ
ÿÿÿÿ v2 ÿÿÿÿ
ÿÿÿÿÿ

The reactions correspond to the first two components of the Fe vector and they are equal
to the opposite of the forces applied on the other end and which correspond to the third
and fourth components.

57
Machine Translated by Google

MEF Course – Chapter 4: Bar Element A. Seghir 2005-2014

The axial force in the bar is obtained with the same reference rotation defined above:

Fx ÿ c Fe1 ÿ s Fe2 = 2831.85 KN .

The calculations give the same value for the second bar.

4.6 Assembly techniques for bar elements

It should be noted when assembling the two elementary systems that components 5 and 6 of the global
matrix K are the sum of components 3 and 4 of the elementary matrices. We therefore see that the
components which correspond to node 3 occupy positions 2ÿ3 ÿ 6 and 2ÿ3ÿ1 ÿ 5 in the global matrix.

Generalizes this assembly by calculating a location table from the

connectivity table as follows:

L([2*i-1 2*i]) =[2*t(i)-1 2*t(i)].

The following MATLAB function file returns a location table in the case of

two degrees of freedom per node:

function L = Locate(t)
% L = Locate(t)
% t: element connectivity table
% L: location table

nne = length(t);
for i= 1:nne

L([2*i-1 2*i]) = [2*t(i)-1 2*t(i)];

end
return

In the case of any number of degrees of freedom, we propose the following function which is more general:

function L = Loc(t,n)
%-------------------------------------------------

% L = Loc(t,n)
% t: element connectivity table

58
Machine Translated by Google

MEF Course – Chapter 4: Bar Element A. Seghir 2005-2014

% L: location table
% n: number of degrees of freedom per node
%

% A. Seghir, 08/24/04 modified 10/23/04

%------------------------------------------------- -
e = eye(n);
L = kron(n*t,e(n,:));
for i=1:n-1
L = L + kron(n*ti,e(ni,:));
end
return

4.7 Finite element program for the bar element

The implementation of a finite element program for bar assembly structures first
requires good organization of the data file. The program must have a structure that
allows it to perform three main tasks: reading data, calculating the solution and
displaying the results.

4.7.1 Data entry

The part of the program intended for reading data must be able to construct variables
for element geometry, material characteristics and external loading.

- The geometry concerns the coordinates of the nodes p, the connectivities of the
elements t and their straight sections A as well as the embedding conditions
which are grouped into a vector e.
- The characteristics of the material are the modulus of elasticity E and the mass
rho volume .

- The external loading is put in the form of a column vector F.

This data can be entered in a text file but for reasons of simplicity offered by MATLAB
functions, or use a function file.
In the case of the previous example of two bars (§2.5), this file is the following:

function [t,p,e,A,E,rho,F] = expl2barres;

59
Machine Translated by Google

MEF Course – Chapter 4: Bar Element A. Seghir 2005-2014

t = [1 3 % connectivity of element 1

2 3 % connectivity of element 2

];

p = [ 0.0 0.6 % coordinates of node 1

1.6 0.6 % coordinates of node 2

0.8 0.0 % coordinates of node 3

];

e = [ 1; 2 % embedding node 1

3; 4 % embedding node 2

] ;

A = [1 1] * 25e-4; % sections of bar elements

E = [1 1] * 210e9; % modulus of elasticity of the elements

rho=[1 1] * 2800; % density 2800 kg/m3

F = [ 0; 0 % node 1 not loaded (embedding)

0; 0 % node 2 not loaded (embedding)
0; -3.4e6 % node 3 loaded in vertical direction
];

If we want to display the structure for visual verification, we use the following
commands:

>> [t,p,b,A,E,rho,F] = expl2bars;

>> plotmesh(t,p,1,1,'b');

The plotmesh function displays a figure of the mesh, it therefore allows, with plotting
the structure entered, a visual verification of the elements and nodes. Its script is:

function plotmesh(T,P, NodeLabels, EltLabels,color)

% plotmesh(T, P, NodeLabels, EltLabels, color)

%

% T = Connectivity table (fem.mesh.t)

% P = Coordinate table (fem.mesh.p)

% NodeLabels: true to display elements

% EltLabels: true to display nodes: line color

% color

% A. Seghir, 02/08/04 modified 07/08/04 Net = size(T,1);

60
Machine Translated by Google

MEF Course – Chapter 4: Bar Element A. Seghir 2005-2014

for ie = 1:Net

XY = P(T(ie,:),:);
X = [[XY(:,1)]; XY(1,1)];
Y = [[XY(:,2)]; XY(1,2)];
line(X,Y,'color',color)
if(EltLabels)

x = mean(XY(:,1));

y = mean(XY(:,2));
text(x,y,num2str(ie),'color','b')
end

end

if(NodeLabels)

Np = size(P,1);
for i=1:Np
text(P(i,1),P(i,2),num2str(i),'color','m')
end
end

4.7.2 Calculation of stiffness and mass matrices

Once the geometry is entered, we can calculate the assembled stiffness and mass matrices. We create a
function truss2dKM for both matrices K and M at the same time.
The input variables common to both matrices are: the connectivity table of the

elements t, the table of coordinates of the nodes t and the sections of the elements A. We add the elastic
modules E for the stiffness matrix K and the densities rho for the mass matrix M. We take different
characteristics for each element (A, E and rho of the vectors) to be able to calculate assembly structures
for bars of several types of material and different sections.

The script for the truss2dKM function is:

function [K,M] = truss2dM(t,p,A,E,rho)

% [K,M] = truss2dKM(t,p,A, rho)

% K: assembled stiffness matrix

% M: assembled mass matrix
% t: element connectivity table

% p: table of node coordinates

% A: sections of elements (bars)
% E: elastic moduli in vector of elements

61
Machine Translated by Google

MEF Course – Chapter 4: Bar Element A. Seghir 2005-2014

% rho: density

% A. Seghir, 08/07/04, modified: 10/26/04

net = size(t,1); nnt = 2*size(p,1); % total number of items

% total number of nodes

K = sparse(nnt,nnt); % initialization of matrices K and M

M = sparse(nnt,nnt); % by zero sparse matrices

for i = 1:net % start of loop on elements

ti = t(i,:); % connectivity tables and

Li = locate(ti); % location of current element

Ke = truss2dKe(p(ti,:),A(i),E(i)); % stiffness matrices

Me = truss2dMe(p(ti,:),A(i),rho(i)); % and elementary mass

K(Li,Li) = K(Li,Li) + Ke; % K matrix assembly
M(Li,Li) = M(Li,Li) + Me; % assembly of matrix M
end % end of loop

return % end of function

It should be noted that we initialized the matrices K and M with sparse matrices . The advantage offered by
sparse matrices is that only non-zero elements are saved in memory. MATLAB is equipped with a data
structure and several functions to support all matrix operations on this type of matrices and in a way that is
practically imperceptible to the user.

The construction of the two matrices K and M is not very different from what was presented for the differential
equations. We have introduced here a localization table Li which we calculate with the localise function
(§2.6) and which we use instead of the connectivity table t for the assembly of K and M.

For the calculation of the elementary matrices Ke and Me we practically only need to enter the components.
The length and direction cosines of the element are calculated from the coordinates of its nodes:

L ÿ (x ÿ 2x ) ÿ1( 2y ÿ y 2) 1
2
; cos(ÿ) ÿ (x2 ÿ x1 )/ L ; sin(ÿ) ÿ ( y2 ÿ y1 )/ L (2.24)

The function files of the two elementary matrices are:

function ke = truss2DKe(XY,A,E)

62
Machine Translated by Google

MEF Course – Chapter 4: Bar Element A. Seghir 2005-2014

% ke = truss2dKe(XY,A,E)
% Calculation of the elementary matrix for a bar element with two nodes (x1,y1)
% (x2,y2)
% A: section of the bar

% E: modulus of elasticity of the material

% XY: coordinates of the nodes XY = [x1,y1; x2,y2]

% A. Seghir, 06/08/04

[c,s,L] = EltLen(XY); % length and orientation of the element

cc = c*c; % cos(angle)^2
cs = c*s % cos(angle)*sin(angle)
ss = s*s; % sin(angle)^2

ke = (A*E/L) * [ cc cs -cc -cs % elementary matrix Ke

cs ss -cs -ss

-cc -cs cc cs

-cs -ss cs ss

];
return

function Me = truss2dMe(XY,A,rho)

% Me = truss2dMe(XY,A,rho)
% Calculation of the elementary mass matrix Me for a bar element with two nodes (x1,y1) (x2,y2)
%

% A: section of the element

% rho: density of the element % XY: coordinates of the nodes XY =

[x1,y1; x2,y2]

% A. Seghir, 07/08/04

L = EltLen(XY);

Me = 0.5*rho*A*L*eye(4);
return

For the distributed mass replace the line Me = 0.5*rho*A*L*eye(4); by :

Me =(rho*A*L/6) * [ 2 0 10

0 2 01

1 0 20

0 1 02

63
Machine Translated by Google

MEF Course – Chapter 4: Bar Element A. Seghir 2005-2014

];

The mass matrix can be calculated by considering a distributed mass or a concentrated

mass. To change, you must put the first expression in a comment line and remove the %
from the second expression. If necessary, it is preferable to introduce an option op variable
in the function header: truss2dMe(XY,A,rho,op) and then choose with an if test what to
return as a result. It should still be noted that in most cases the result is approximately the
same and it is for its diagonal matrix configuration that the concentrated mass is preferred.

Finally, we give the source code of the EltLen function :

function [ds,c,s] = EltLen(XY)

% [ds,c,s] = Getcosine(XY)
% Calculate element length, cos and sin for matrix
% rotation of a bar element with two nodes (x1,y1) (x2,y2)
% ds: element length
% c: cosine of the angle between the element axis and the x axis
% s: sine of the angle
% XY: coordinates of the nodes XY = [x1,y1; x2,y2]
%

% A. Seghir, 06/08/04

dx = XY(2,1) - XY(1,1);
dy = XY(2,2) - XY(1,2);
ds = sqrt(dx^2 + dy^2);
c = dx/ds;
s = dy/ds;

return

4.7.3 Application of support conditions

The support conditions are the set of zero displacements at the support levels.
The vector e in the data file is used to specify the degrees of freedom to block. To apply
this condition we eliminate the rows and columns of the matrices

64
Machine Translated by Google

MEF Course – Chapter 4: Bar Element A. Seghir 2005-2014

K and M as well as the vector F since the loading is entered for all the nodes including
the support ones.
Care should be taken when applying the MATLAB commands A(i,:)=[] and A(:,i)=[]

which allow you to delete rows and columns from a matrix. With each column or row
deleted the size of the matrix is reduced and the indices are shifted to the left in the
case of a column and upwards in the case of a row.
Thus, before proceeding to delete a list of degrees of freedom, it must first be sorted
in an ascending or descending order and then delete the rows and columns by going
through it from the largest to the smallest index. This way we don't have to worry about
offsets in the list values.
The MATLAB command which allows you to order a list or a vector L is: sort(L).

The function which allows the support conditions to be applied is called DelDOFs and
expects as input arguments a variable A which can be a square matrix or a column
vector, and a list L of the degrees of freedom to be deleted. Its script is as follows:
function A = DelDOFs(A,L)
% A = DelDOFs(A,L)
%

% A: overall matrix after assembly

% L: list of degrees of freedom to eliminate
%

% A. Seghir, 01/08/04 modified 08/08/04

L = sort(L); % sorting of the list of DDLs

n = length(L); % length of L

if (size(A,2) == 1) % case of a vector

for i = n:-1:1 % browse L from most
big index
A(L(i)) =[]; % removal of component
associated with the current DDL
end
else % case of a matrix
for i = n:-1:1
A(L(i),:) =[]; % deletion of line L(i)
A(:,L(i)) =[]; % column deletion
L(i)
end

65
Machine Translated by Google

MEF Course – Chapter 4: Bar Element A. Seghir 2005-2014

end
return

4.7.4 Resolution of the discrete system

Once the boundary conditions are applied, all that remains is to solve the discrete system.
We have two types of analysis:
1) Static analysis which consists of determining the displacements and internal forces
in the bars under the effect of static loading
2) Dynamic analysis which consists of determining the dynamic behavior of the
structure and its response to a response spectrum or an accelerogram.
In the present case we will limit ourselves to the static analysis and the simple calculation
of the eigenmodes. The dynamic response will be the subject of a separate chapter. The
displacement solution is obtained with the command U = K\F and the axial forces in the
bars (internal forces) are obtained with the elementary systems Fe ÿ Ke U. For each
element it is necessary to recalculate the stiffness matrix K, extract the elementary
displacements from the solution U and carry out the product Ke Ue and finally apply a
reference rotation to return to the local reference of the bar. The vector Ue must also
include the zero values of the movements blocked at the supports (see the example).
Finally, the reactions and the loading are the sum of the elementary forces of all the
elements linked to the supports. The MATLAB function that allows you to perform these operations is:
function [F,R] = TrussForces(t,p,A,E,U)

% [F,R] = TrussForces(t,p,A,E,U)
% F: axial forces in the bars

% R: forces at the nodes = reactions in the case of supports

% t: element connectivity table

% p: table of node coordinates

% A: sections of elements

% E: moduli of elasticity

% U: solution on the go

% A. Seghir, 08/07/04 modified: 08/27/04

net = size(t,1); % total number of items

nnt = 2*size(p,1); % total number of nodes

R = zeros(2*nnt,1); % Forces at the nodes

for ie = 1:net % loop over elements

66
Machine Translated by Google

MEF Course – Chapter 4: Bar Element A. Seghir 2005-2014

L = locate(t(ie,:)); % location table

ke = truss2dKe(p(t(ie,:),:),A(ie),E(ie)); % matrix
elementary

ue = U(L); % node movements

fe = ke*ue; % elemental forces in

(oxy)
[Len,c,s] = EltLen(p(t(ie,:),:)); % direction cosines
F(ie,:) = -(c*fe(1)+s*fe(2)); % rotation at local reference
R(L) = R(L) + fe; % all elements related to
node

end % end of loop

return

We can calculate the eigenmodes of the structure with the eigs function which
ÿ
system ÿK ÿ ÿ resolves M ÿ ÿ ÿ 0. It takes as input argument the two matrices K and M and
returns two matrices as a result: one orthogonal for the eigenvectors and the other
diagonal for the eigenvalues. The eigenmodes are therefore the columns of the
eigenvector matrix and the eigenperiods are given by: T ÿ 2ÿ/ÿ.

4.7.5 Main program

Now all the functions necessary for the implementation of a program for calculating
lattice structures have been presented. All that remains is to preferably place the
assemblies in a function file like this:
function [U,P,R, T,phi]=trussfem(ffd)
% [U,P,T,phi] = trussfem(ffd)
% resolution of two-dimensional bar assembly systems

% U: solution in nodal displacements

%P : axial forces in the bars

% R: reactions to support
%T : specific periods of the structure
% phi: eigenmodes of the structure

%ffd: problem data function file

% A. Seghir, 08/06/04 modified 10/27/04

[t,p,e,A,E,rho,F] = feval(str2func(ffd)); % appeals to

function file

% containing data

plotmesh(t,p,1,1,'b'); % displays structure

loaded

67
Machine Translated by Google

MEF Course – Chapter 4: Bar Element A. Seghir 2005-2014

[K,M] = truss2dKM(t,p,A,E,rho); % calculation of matrices K and

M
K = DelDOFs(K,e); % Application of
terms
M = DelDOFs(M,e); % at limits
F = DelDOFs(F,e);
U = K\F; % static solution % addition of
U = AddDOFs(U,e); node DDLs
recessed
[P,R] = TrussForces(t,p,A,E,U); % forces in the bars and
the knots
[phi, omega2] = eigs(K,M); omega = % clean modes
sqrt(diag(omega2)); T = 2*pi./sort(omega); % pulsations
% clean periods per
Ascending

% display results
net = size(t,1); nnt = size(p,1);
disp('Calculation results of the 'bar assembly' structure)
disp([' data file: ',ffd])
disp([' Number of bars: ',num2str(net)])
disp([' Number of nodes: ',num2str(nnt)])

disp(sprintf('\n Movements to nodes:'))

disp(sprintf(' Node\t\t Ux\t\t\t\t Uy '))
for i=1:nt
disp(sprintf(' %d\t\t\t%+5.4f\t\t\t%+5.4f',i,U(locates(i))))
end

disp(sprintf('\n Efforts in the bars:'))

disp(sprintf(' Bar \t\t P'))
for i=1:net
disp(sprintf(' %d\t\t\t%+1.4E',i,P(i)))
end

disp(sprintf('\n Reactions to presses:'))

disp(sprintf(' Press\t\t Rx\t\t\t\t\t Ry '))
for i=1:nt
L = locate(i);
if find(e==L(1))

68
Machine Translated by Google

MEF Course – Chapter 4: Bar Element A. Seghir 2005-2014

disp(sprintf(' %d\t\t\t%+1.4E\t\t%+1.4E',i,R(L)))
end
end

disp(sprintf('\n Structure's own periods:'))

disp(sprintf(' mode\t\tT'))
for i=1:size(T)
disp(sprintf(' %d\t\t\t%5.4f',i,T(i)))
end
shapeplot(U,p,t,10) % draws the structure
distorted

Notice how the data file is called by this function. The ffd variable

must be of character type (string), str2func allows you to construct a function address from a name and feval evaluates it. To

use trussfem, you must therefore prepare a data file similar to expl2barres.m and use the command:

>> [U,P,T] = trussfem('expl2barres');

Executing it displays the following lines as results:

Bar connection structure calculation results
data file: expl2barres
Number of bars: 2
Number of nodes: 3

Movements at nodes:
Node Ux Uy
1 +0.0000 +0.0000
2 +0.0000 +0.0000
3 +0.0000 -0.0090

Efforts in the bars:

Rod P
1 +2.8333E+006
2 +2.8333E+006

Reactions to support:
Supports Rx Ry
1 -2.2667E+006 +1.7000E+006
2 +2.2667E+006 +1.7000E+006

69
Machine Translated by Google

MEF Course – Chapter 4: Bar Element A. Seghir 2005-2014

Specific periods of the structure:

Fashion T
1 0.0011
2 0.0014

4.8 TP No. 3
Model and calculate the structure below using the trussfem program and
SAP2000 software . All the bars of the structure have a cross section A ÿ
6cmÿ6cm and are made of a material with a density ÿ ÿ 7800 Kg/mÿ and a
modulus of elasticity E ÿ 210000 MPa. Distances are given in meters and loads
in KN . Give the displacements of the nodes, the reactions to the supports and
the internal forces in the bars.

15KN _ 15KN _ 15KN _

10KN _
4 7 6 1 8
1 1
4 6 8 1
0
3 2
7
3 5 5 9

2
1

1 1.0 1.0 1.0

Fig. 2.6: Structure to calculate in

TP No. 3

70
Machine Translated by Google

chapter 5

Beam Element

5.1 Governing equation

The beam element is used to take on, in addition to the axial force like the bar element, a loading
perpendicular to its axis. Beams are found in many civil engineering and mechanical construction structures.
The most common cases are porticos constituting residential buildings, bridges, etc.

We consider beams to be slender pieces (in reinforced concrete or steel), which have a very large dimension
compared to the other two and which generally work in
bending.

The formulation of the beam element can be obtained based on the theory of material strength; we consider
a beam of section A and length L subjected to a loading q(x) varying along its longitudinal axis as shown in
figure 3.1 below:

y
q(x)

Fig. 3.1: loaded beam

Under the effect of loading, the beam bends and moves vertically by a displacement v(x). It is assumed that
after this deformation, the straight sections remain flat and
Machine Translated by Google

MEF Course – Chapter 5: Beam Element A. Seghir 2005-2014

perpendicular to the mean line; they therefore undergo a small rotation of angle ÿ in the
(oxy) plane. Consider an element dx of the beam delimited by two neighboring sections,
one straight and the other inclined.

y
u ÿ

not

y
ÿ
v(x)

ÿ
ÿ vx( )
dx vx( ) ÿ dx
ÿx
v(x)
x

Fig. 3.2: deformation of a section

The rotation of the deformed section is the tangent of the curved mean line:

v
ÿÿ (3.1)
ÿ ÿx

Because of the rotation, the points of the section undergo a horizontal displacement u
varying linearly from the lower fiber to the upper fiber. At a point in the section this
displacement is worth:

v
u ÿ ÿ ÿ ÿ ÿÿ xÿ yy (3.2)

where y denotes the distance from the mean line (center of the section)
Within the framework of the hypothesis of small deformations, the axial deformation along
x along the section is:

72
Machine Translated by Google

MEF Course – Chapter 5: Beam Element A. Seghir 2005-2014

u
ÿÿÿ (3.3)
x
ÿx

If we designate by E the modulus of elasticity of the material of the beam, Hooke's law gives the
distribution of stresses along the section:

ÿu ÿ 2v _
ÿ ÿx ÿ ÿ E x
E ÿÿ
Ey2ÿx (3.4)
ÿx

The moment created by these constraints must balance the bending moment M created by the external
loading:

ÿ 2v _ 2
ÿ 2v _
2
ÿ ÿE 0. ÿ.
M yÿds ÿÿ ÿ y ds ME ÿ ÿ
ÿ y ds (3.5)
s 2x _ s 2x _ s

with s designates the area of the cross section. We set Iz the moment of inertia with respect to
2
the z axis perpendicular to the plane (xy) I yÿ ds
ÿ , the expression of the moment becomes:
s

ÿ 2v _
M EI ÿÿ (3.6)
2x _

Now consider the static equilibrium of an element dx.

q
M
MÿdM
x

T TÿdT
dx

Fig. 3.3: Static equilibrium

The sum of the moments relative to its center of gravity gives:

73
Machine Translated by Google

MEF Course – Chapter 5: Beam Element A. Seghir 2005-2014

dx dx
ÿ ÿTM
ÿÿ2 ( T dT ) ( ÿ M ÿ dM )ÿ 0
2

dx
Tdx dT
ÿ ÿ dM ÿ 0
2

By neglecting the second order terms, we obtain the relationship between the shear force and
the bending moment:

dM ÿ ÿ 2v
T ÿÿÿÿ
( EI ) (3.7)
dx ÿx ÿx 2

Balance of vertical forces for positive loading in the direction of the y axis (figure 3.3):

T ÿ dT ÿ q dxÿT ÿ 0 (3.8)

gives the relationship between the loading q and the shear force T and relates the loading to
the displacement v by:

2 2
dT ÿ ÿ v
q ÿÿÿ

2 (
EI
2)
(3.9)
dx x ÿ ÿx

This equation translates the static balance of the beam. In the case of a dynamic movement, a
term reflecting the inertia forces must be added to equation (3.8):

2
ÿ u
F im A dx ÿ ÿ ÿ ÿ 2
(3.10)
ÿt

with ÿ is the density of the material, A the section of the beam and t represents the time.

Equations (3.8) and (3.9) become:

T ÿ dT ÿ q dxÿT ÿ Fi (3.11)

ÿ 2v ÿ ÿ 2v
ÿA ÿ ( EI ) ÿ
qx( ) (3.12)
ÿ tx2 2 ÿ ÿx 2

Equation (3.11) is the Euler-Bernoulli equation for beam bending. The displacement v is a
function of the coordinate x along the axis of the beam and the time t.

74
Machine Translated by Google

MEF Course – Chapter 5: Beam Element A. Seghir 2005-2014

5.2 Element wording

5.2.1 Variational formulation

By denoting the length of the beam by L and taking ÿv the weight function, the strong
variational formulation associated with equation (3.12) is written:

L 2v _
L ÿ2 ÿ 2v _ L
ÿ ÿ ÿÿ v A 2 ÿ dx ÿ ÿÿ v ( IS ) dx ÿ ÿ ÿ vqx (dx
) (3.13)
0 t 0 ÿ 2x _ ÿ 2x _ 0

The weak integral form is obtained with two integrations by parts of the second term.
L
L 2 L
ÿ ÿ 2v _ ÿÿ v ÿ ÿ 2v _ ÿ ÿ 2v _
v ( IS ) dx ÿÿ ÿ ( IS ÿ ) dx ÿÿÿÿ v ( IS (3.14a)
ÿ 0ÿ ÿ 2x _ ÿ 2x _ 0 ÿ xx ÿ ÿ 2x _ ÿx ÿ 2x _
ÿ )ÿ 0ÿ

L
L ÿÿ v 2
ÿ ÿ 2v _ L ÿÿ v ÿ 2v _ ÿ ÿÿ ÿ ÿ v ÿ 2v _ ÿ

ÿ ( IS ) dx
ÿÿ
ÿ IS dx ÿ IS (3.14b)
0 ÿ ÿ xx ÿ 2x _ 0 ÿ 2x _ ÿ 2x _ ÿx _ ÿ 2x _
ÿÿ
0

Taking into account expressions (3.6) and (3.7), the second terms represent the difference
in the loadings in forces (T0 and TL) and in moments (M0 and ML) applied to the ends of the
beam. Furthermore, we can replace the derivative of the disturbances of the displacements
by a disturbance of the rotations (equation 3.1):ÿÿv ÿx ÿ ÿÿ .

We write the boundary conditions as follows:

L
ÿ ÿÿ ÿ v ÿ 2v _
IS ÿ ÿ ÿ ÿÿ ÿ M0 ÿ ÿÿ M
xÿ 0 x LL
ÿ

ÿÿ x ÿ 2x _
0

(3.15a)
L
ÿ ÿÿ( 2v _
ÿv _ EI ÿ ÿ ÿ ÿ v T v T x LL 0 0
2ÿx ÿ) ÿ
xÿ
ÿx
ÿÿ ÿÿ 0

(3.15b)

Now substituting (3.15) into (3.14) and the result into (3.13) we obtain the expression for the
weak variational form:

75
Machine Translated by Google

MEF Course – Chapter 5: Beam Element A. Seghir 2005-2014

L ÿ 2v L ÿ 2ÿ v ÿ 2v
ÿ ÿ go
ÿ dx ÿ ÿ EI dx M
ÿ ÿÿ ÿ
0 ÿt2 0 ÿx 2 ÿx 2 x L.L.
ÿ

(3.16)
L
ÿ ÿ ÿ ÿ ÿÿ v T v TM vqx dx ÿ ÿÿ ()
x L.L.
ÿ
x 00ÿ
x 0
ÿ 0
0

5.2.2 Discretization
For the discretization of this equation we consider an element with two nodes: a node at each end
of the beam. The presence of derivatives of order two requires the use of quadratic polynomials or
more. Furthermore, we see that the expression of the boundary conditions involves rotation at the
ends, it is therefore more interesting to take two degrees of freedom per node in order to ensure at
the same time the continuity of the displacements and their derivatives which are the rotations. The
number of degrees of freedom thus reaches four and the interpolation polynomial must be cubic
(four constants).
The vector of elementary displacements and rotations is therefore written as follows:

T
U ÿÿ v1 ÿ1 v2 ÿ2 ÿ
not
(3.17)

v1 (F1) v2 (F2)

ÿ1 (Mÿ) ÿ2
(M2)

1 2

x1 ÿ 0 x2 ÿ L
Fig. 3.4: Beam element with two nodes

The displacements and rotations along the beam are approximated by:

2 3
v(x) ÿ a 0ÿ ax 1ÿ ax ÿ2 ax 3 (3.18a)

2
ÿ(x) ÿ a ÿ1ax ÿ2 ax2 3 3 (3.18b)

The evaluation of these polynomials at the nodes gives:

v(0) ÿ a 0ÿ v 1
; ÿ(0) ÿ a 1ÿ ÿ 1 (3.19a)

76
Machine Translated by Google

MEF Course – Chapter 5: Beam Element A. Seghir 2005-2014

2 3 2
v(L) ÿ v 2ÿ v ÿ1 ÿ L 1ÿ a L2ÿ a L 3 ; ÿ(L) ÿ ÿ ÿ 2ÿ ÿ a 1L ÿ 2a L2 3 3 (3.19b)

Solving (3.19b) for a2 and a3 gives:

3 1 2 1
ÿ
( vv2 1 ) ÿ ÿ1 (2
ÿ ÿ2 ) ; ( vv1 2 (3.20)
ÿ

has
2 has
3
ÿ ÿ

) ÿ ( ÿ1 ÿ ÿ2 )
2 3 2
L L L L

Replace these parameters in equation (3.18) and after arrangement of the terms we
write the nodal interpolation of the displacements in the form:

v(x) ÿ N 1(x) v ÿ1 N (x)

2 ÿ ÿ N
1 (x) 3v ÿ N2 (x) ÿ4 2 (3.21)

The Ni shape functions are called Hermite polynomials, their expressions are:
2 3 3
3x 2x 2 2x x
N x1 (ÿ)ÿ1ÿ ; N x (ÿ )x ÿ 2 ÿ ;
2
L L
3
L L2

x2 2x 3
x
3
x
2
N x3 3ÿ ( ) ÿ

; N 4x( ) ÿ ÿ (3.22)
L2 L3 L
2
L

We can verify that the sum of the shape functions associated with the displacements is
equal to unity: N1ÿN3 ÿ1. They also take values equal to one at the nodes that
correspond to them and zero values at the opposite nodes. This remark is not valid for
the functions associated with rotations since the rotations are themselves derivatives
of the displacements. The MATLAB program which allows you to calculate and plot the
Hermite form functions is as follows:
clear, clc
syms x L real % declare real symbolic x and L

P = inline('[1 xx^2 x^3]') % 3rd degree polynomial dP = inline(diff(P(x)))

% derivative of polynomial
Pn = [P(0); dP(0); P(L); dP(L) ] % evaluation at knots

N = inline(( P(x) * inv(Pn))) dN = inline((dP(x) * % shape functions and

inv(Pn))) % their derivatives

% graphs of N and dN
t = 0:0.01:1;
subplot(2,1,1), plot(t, N(1,t')), title(' Functions of form N1 N2 N3
N4')
subplot(2,1,2), plot(t,dN(1,t')), title(' Derivatives dN1 dN2 dN3 dN4 ')

77
Machine Translated by Google

MEF Course – Chapter 5: Beam Element A. Seghir 2005-2014

The figure below shows the curves of the shape functions and their derivatives for a
beam element of unit length (L ÿ 1):

Fig. 3.5: Shape functions of a beam element and their derivatives

5.2.3 Elementary matrices

We now replace the displacement v in the variational form (3.16) by its approximation
(3.21), we obtain for the disturbances:
T
TT T
dN
ÿv ÿ ÿAn N ; ÿÿ ÿ ÿU with ÿU Tÿ ÿ ÿv1 ÿÿ1 ÿv1 ÿÿ2 ÿ (3.23)
dx
not not

The derivatives become:

ÿ 2v 2
dN 2
ÿ v of
2
not
ÿÿ

2
ÿ

2 A ; 2
ÿ
NOT
2
ÿ
NAKEDnot (3.24)
ÿx dx ÿt dt

78
Machine Translated by Google

MEF Course – Chapter 5: Beam Element A. Seghir 2005-2014

2
dN 1
with : ( x) ÿÿÿ 12 6 , 6 6 12 x4 L 2Wx
ÿ

, LL x Lx L
ÿ

, 6 ÿ

2 2
ÿ
(3.25)
dx 2
L3

Taking into account the values of the functions N and dN at the nodes (Fig. 3.5), the conditions at
limits are written:

T T T
ÿ ÿ ÿ v T v TU x LL 0
ÿ
x 0
ÿ
ÿÿÿ not (0010 ÿ TL ÿÿÿ 1000 T0 )
T T
(3.26a)
ÿ ÿ ÿ UT
ÿ not 0 0 TL 0ÿ

ÿÿ Mx ÿ ÿÿ MU 0 0 ÿ ÿ ÿ ( 0T 0
01 ÿ T M0100ÿÿÿ T
M )
LL
ÿ
x ÿ not L 0
(3.26b)
ÿ Tÿ 0
0 ÿ ÿUMM ÿ T
not 0 L

These conditions correspond to the external loading applied to the nodes, that is to say
at the junctions between elements such as post-beam junctions for example. The elementary
vector Fn of the forces and moments concentrated at the nodes is written as follows:

T
Fc ÿ ÿ ÿ T0 ÿ M0 TL ML ÿ (3.26c)

The loading Fe distributed on the beam element corresponds to the second term of the variational
form (3.16):
L L L
T T T
( ) dx UNÿ qx
ÿ ÿvqx ÿ dx ÿ not
() ; Fe ÿNÿqx dx ( ) (3.27)
0 0 0

This loading is a function of the distribution q(x). In simple cases we can integrate and give the
explicit expression of Fe, otherwise, it is possible to subdivide the beam (or the element) into
several small elements so as to be able to replace the distributed loading q by two equivalent
concentrated forces at the extremities.

We give below the expression of Fe for a case of trapezoidal loading varying from q0 to qL.
Uniform and triangular loadings are only particular cases of trapezoidal loading.

x
q(x) q (qL q )ÿ ÿ 0 ;
ÿ

0 (3.28a)
L

79
Machine Translated by Google

MEF Course – Chapter 5: Beam Element A. Seghir 2005-2014

L T
F eÿ ÿ (9 ÿ 2 1 ) , (2 ÿ 3 ) , q(20L1qL
ÿ 9q ) , qL 0 qL q 0 ÿ (3 ÿ 2q) ÿ
L qL 0 (3.28b)
60

There now remain two integral terms in equation (2.16) to discretize:

2 2 2 T 2
L ÿÿ v ÿ v L
T dN dN T
ÿ 2
EI 2
dx U 2
EI 2
U dx (3.29a)
0 ÿx ÿx ÿ ÿ ÿ0 not

dx dx
not

L ÿ 2v L
TT
ÿÿ

ÿ ÿ go
ÿ dx UNNU dx
ÿ ÿ ÿ0
ÿ (3.29b)
ÿt2
not not
0

T
The elimination of ÿA of the discrete integral equation gives the two mass matrices
Me and Ke rigidity :

L d 2N T
d 2N
Ke EI dx (3.29)
ÿ ÿ0 dx 2 dx 2
L
T
ÿ ÿ dx
Me NAN ÿ (3.30)
0

If the beam is made of the same homogeneous material and of the same section (E, I and ÿ are
constant), then the explicit expressions of the elementary matrices can be obtained and are
written as follows:

ÿ 12 6 L ÿ

12 6 Lÿ ÿ 156 22 L 54 13 ÿ

Lÿ
2
EI LL 6 LL
2 2 22LL 13 L 3L 2
ÿ ÿ

4
ÿ ÿ

26 4
ÿ AL
ÿ ÿ

Ke
ÿ ÿ ÿ

; Me ÿ ÿ ÿ

L3 ÿ
ÿÿ

12 6 L 12 6 ÿ

L ÿ 420 ÿ 54 13 L 156 22 ÿ

L ÿ

2 2
LL 6 LL
ÿ ÿ ÿ ÿ

262 ÿ

4 ÿ

2 13LL
3 ÿ ÿ

22 L 4L ÿ
ÿ ÿ ÿ

(4.31)

The same remark concerning the mass matrix of the bar element can be made for the beam
element. We can associate half of the total mass of the element with the degrees of translation
of each node. We only take translations because rotations do not produce inertial forces. The
concentrated mass matrix is therefore written as
follows:

80
Machine Translated by Google

MEF Course – Chapter 5: Beam Element A. Seghir 2005-2014

ÿ 0ÿ 0
1000ÿÿ
AL
ÿÿÿÿÿ0010ÿÿ (4.32)
Me ÿ

2 ÿ0000ÿ

It should be noted that the two concentrated and distributed matrices give the same total mass
with the sum of the components associated with the degrees of freedom of translation: M(1, 1) ÿ
M(1, 3) ÿ M(3, 1) ÿ M(3, 3) ÿ ÿAL.

5.3 Example of a console beam

5.3.1 Case of concentrated loading at the free end

A reinforced concrete beam with density ÿ ÿ 2.5 tonnes/ mÿ and modulus of elasticity E ÿ 32000
MPa is embedded at one of its ends and loaded at the other end with a concentrated load Fy ÿ
ÿ90 KN and a moment M ÿ 60 KN m

The geometric characteristics of the beam are:

Fy
- length L ÿ 150 cm
h M
- width b ÿ 30cm b
L
- height h ÿ 40 cm Fig. 3.6: Console loaded at the end
ÿ
The inertia of the beam is: I ÿ 3 ÿ 4 / 12 = 16 cmÿ ;

The stiffness matrix is written:

ÿ
86 ÿ

86ÿ66 ÿ
6
3.2 16 3 ÿ 3ÿ ÿ

K e ÿ MN / m
3 1.5 ÿÿ
ÿÿ

86 ÿ86ÿ

ÿ
ÿ
6366ÿ ÿ

2ÿÿ

3 T
The force vector according to equation (3.26) is: F ÿ10 ÿ 0 0 ÿ 90 60ÿ

The elementary system is written:

81
Machine Translated by Google

MEF Course – Chapter 5: Beam Element A. Seghir 2005-2014

ÿÿ
86 ÿ86ÿÿÿ v1
ÿ

ÿ ÿ
0
ÿ
3.2 1ÿ6 3ÿ 6663ÿÿÿ86ÿÿÿÿ ÿ1 ÿ

0
6 ÿ
K ÿ ÿ
10 ÿ
ÿ ÿ
ÿ
3 ÿÿ ÿ

1.5 2 ÿÿ

86 v2 ÿ ÿ

ÿ ÿ ÿ 90000 ÿ

ÿÿ
6366 ÿ

ÿ2
ÿÿÿÿ ÿÿÿÿÿ
60000 ÿÿ

The embedding conditions (v1 ÿ 0; ÿ1 ÿ 0) reduce the system to two equations:

ÿ86ÿÿÿ
ÿ

v2 ÿ ÿÿ ÿÿ90000
ÿ 60000
ÿ
7ÿ
K ÿ
2.27556 10 ÿ

6 ÿ
ÿ
ÿ
ÿ

6ÿÿÿ
ÿÿ 2 ÿ

The solution gives v2 ÿ ÿ0.6592 mm ; ÿ2 ÿ ÿ0.2197ÿ 10ÿÿ ; a negative downward displacement

and a negative rotation in the opposite direction to the trigonometric direction (Fig. 3.4)
The basic efforts are:

ÿÿ
86 ÿ86ÿ ÿ
ÿ

0 ÿ
F1 ÿ ÿ 90ÿ ÿ

ÿÿ ÿ
6663ÿÿÿ86ÿÿÿ6
ÿ

0 ÿ

3
M1 75
K ÿ
7 2.27556 1 0 ÿ

3 ÿ 10 ÿ
ÿ ÿÿ
ÿ
310 ÿ ÿ
ÿÿ

86 66ÿÿ ÿ ÿ

ÿ ÿ ÿ 0.6592 F2_ ÿ

90ÿ
ÿ

ÿ ÿ

ÿ ÿ 0.2197 M 60ÿ
ÿÿÿÿ ÿ ÿÿÿÿÿ 2 ÿÿÿÿ ÿÿÿÿÿ

The first two components are the supporting reactions:

Ry = 90 KN ; Mz ÿ 90 ÿ 1.5 ÿ 60 ÿ 75 KN m.
The moment is positive because of the sign convention adopted in Figure 3.4: At node 1, the
positive moment stretches the upper fiber to rotate in the trigonometric direction. The last two
components correspond to the nodal forces applied to the free end.

It should be noted that the solution obtained, whether for the arrow or the reactions, is an exact
3
solution. This is due to the fact that the arrow is a function in x of form are and functions
polynomials of the third degree. It follows that the Hermite polynomials make it possible to
obtain the exact solution.

5.3.2 Case of distributed loading

82
Machine Translated by Google

MEF Course – Chapter 5: Beam Element A. Seghir 2005-2014

We now consider the same previous load distributed uniformly along the beam, i.e.:

q ÿ 90/1.5 = 60 KN / ÿ. ; q ÿ 60 KN / mÿ

and we keep the same moment at the free end.

M ÿ 60 KNm

Fig. 3.7: Distributed loading

The force vector due to loading q, is written, according to equation (3.28),

(with q0 ÿ q1 ÿ q ÿ ÿ60 KN / mÿ):

3 T
Fq ÿ10 ÿ ÿ 45 ÿ11.25 ÿ 45 11.25ÿ (check the sum of Fi )

3 T
The vector concentrated forces from the moment is: FC ÿ10 ÿ 0 0 0 60ÿ
T
The total force vector is therefore: F ÿ10 ÿ ÿ3 45 ÿ11.25 ÿ 45 71.25ÿ

The rigidity matrix being the same, the solution of the reduced system gives:

v2 ÿ 0.5768 mm and ÿ2 ÿ 1.0986 ÿ 10ÿÿ .

The elementary reaction vector:

T
R ÿ KU ÿ10 ÿ 453 ÿ3.75 ÿ 45 71.25ÿ

These efforts are represented in Figure

3.8 opposite.

They are a part of the reactions F2 ÿ ÿ450 KN

F1 ÿ 450 KN
which correspond to the displacement
and rotation of the end M2 ÿ71.25 KN m
Mÿ ÿ ÿ3.75 KN m

free. The presence of the force q

distributed over the element gives rise
Fig. 3.8: Efforts at the nodes
to forces and

nodal moments (vector F) to which the

second part of the reactions a

corresponds.

83
Machine Translated by Google

MEF Course – Chapter 5: Beam Element A. Seghir 2005-2014

The total reaction therefore corresponds to the sum of the two types of reactions:

Ry ÿ Rÿ ÿ (ÿFqÿ) ÿ 90 KN. ; Mz ÿ R2 ÿ (ÿFqÿ) ÿ 7.5 KN m.

At the free end the sum of the two reactions cancel each other out in force and their sum equals the
applied moment:

M ÿ Rÿ ÿ (ÿFqÿ) ÿ 60 KN m.

Generally speaking, the moment acting on a cutting section at a distance x

of the embedding is:
Ry
ÿ
q (x)
M(x) ÿ ÿMz ÿ Ry x ÿ ½ qx M(x)
Mz
ÿ
ÿ ÿ7.5 ÿ 90 x ÿ 30 x ÿKN mÿ
x
at distance x ÿ 0.5 m : M(0.5) ÿ 30.0 KN m
Fig. 3.9: The moment in the element
at distance x ÿ 1.0 m : M(1.0) ÿ 52.5 KN m

at the distance x ÿ 1.5 m : M(1.5) ÿ 60.00 KN m.

3
This moment cancels at the position xÿÿÿ 2 0.08579
2
Note : In the case of a trapezoidal load (expression 3.28a), it is necessary to replace by the sum of
ÿ
½ tbsp the moment of a uniform load q0 and that of a load

triangular (qL ÿ q0) xÿL or:

2
x3
Mq qx 1 2 0
ÿ
( ÿ qL
ÿ q 0)
6L

The concentrated and distributed mass matrix of the beam are:

ÿ 225 0 ÿ 00 ÿ ÿ 624 132 216 ÿ

ÿ

78 ÿ
00 00 132 36 78 27
ÿ
ÿ

450 ÿ

Mc
ÿ ÿ ÿ

; Mr. ÿ ÿ ÿ

ÿ 0 0 225 0 ÿ 1680 ÿ 216 78 624 132 ÿ

ÿ ÿ

00 00
ÿ ÿ

ÿ ÿ ÿ
ÿ

78 27 132
ÿÿ

36 ÿ

The natural periods of the beams are given by T ÿ 2ÿ/ÿ with ÿ solution of
ÿ
the equation: det(K ÿ ÿ M) ÿ 0.

For a concentrated mass matrix: T1 ÿ 0.01397 s ; T2 ÿ 0.0 s

84
Machine Translated by Google

MEF Course – Chapter 5: Beam Element A. Seghir 2005-2014

For a coherent mass matrix: T1 ÿ 0.00968 s ; T2 ÿ 0.00098 s

5.3.3 MATLAB program

If we want to do these calculations with MATLAB, we can provide a script similar to this one
this :

%------------------------------------------------- ------------------

% console.m

% script to solve the problem of a console beam

%------------------------------------------------- ------------------

clear, clc % clears variables and window

L = 1.5; % Console length

b = 0.3; % Section width
h = 0.4; % Section Height

E = 3.2e10; % Elasticity module

rho = 2.5e03; % Volumic mass

A = b*h; % Section

I = b*h^3/12; % Inertia

q = -60000; % Distributed load

P = -90000; % Concentrated load
Mz=60000; % Concentrated moment

Fq = q*L/60*[30 5*L 30 -5*L]'; % Force vector associated with q

Fm = [0 0 0 Mz]'; % Force vector associated with Mz
Fp = [0 0 P 0 ]'; % Force vector associated with P

F1 = Fm + Fp; % First loading case

F2 = Fm + Fq; % Second loading case

% Stiffness matrix
K = E*I/L^3 * [ 12 6*L -12 6*L
6*L 4*L^2 -6*L 2*L^2
-12 -6*L 12 -6*L

6*L 2*L^2 -6*L 4*L^2];

% Concentrated mass matrix

Mc = rho*A*L/2 * diag([1 0 1 0]);

85
Machine Translated by Google

MEF Course – Chapter 5: Beam Element A. Seghir 2005-2014

% Coherent mass matrix

Mr = rho*A*L/420 * [ 156 22*L 54 -13*L
22*L 4*L^2 13*L -3*L^2
54 13*L 156 -22*L
-13*L -3*L^2 -22*L 4*L^2];

U1 = K([3 4],[3 4])\F1([3 4]) % solution for the 1st case

U2 = K([3 4],[3 4])\F2([3 4]) % solution for the 2nd case
R1 = K*[0;0;U1] % Reactions for the 1st case
R2 = K*[0;0;U2]-Fq % Reactions for the 1st case

syms x real; % Moment versus x

Mzx = inline(-R2(2) + R2(1) * x + 0.5* q *x.^2)
Mzx([0 0.5 1 1.5]')

% clean periods case of a concentrated mass

Tc = 2*pi*eig(K([3 4],[3 4]),Mc([3 4],[3 4])).^-0.5

% proper periods case of a distributed mass

Tr = 2*pi*eig(K([3 4],[3 4]),Mr([3 4],[3 4])).^-0.5

5.3.4 SAP model

The SAP model of this console beam can be produced in the following way:
1) Create a Frame element between the points (X0.0,Y0.0) and (X1.5,Y0.0). See the steps
1 to 6 described for example 2.5. Embed the node at point (0,0)

2) In the material definition window, put 3.2e10 for Elasticity Modulus,

2500 for Mass per Unite Volume and 0 in all other input fields.

3) For the section, you must put 0.3 for Depth (t3) and 0.4 for Width (t2) since the model is
planar and the height of the beam is in the Y direction instead of the Z direction.

4) Define the three static load cases FP, FQ and FM in the Define Static Load Case Names window ; choose the Live type

and enter 0 for Self Weight Multiplier.

5) Once the loads have been defined, you can make the FM+FQ and FM+FP combinations
with the Define/Load Combinations menu. Click on the Add New Combo button. In the
Load Combination Data window , enter FQM for both Load Combination Name and
Title fields, leave ADD for Load Combination Type. In the

86
Machine Translated by Google

MEF Course – Chapter 5: Beam Element A. Seghir 2005-2014

Define Combination pane , choose FQ Load Case for Case Name and keep 1 as Scale
Factor then press the Add button. Add FM in the same way and validate with OK. Now
define the FPM combination in the same way .
Use the Delete button to erase previous loads or the Modify button for modified ones.
Close both windows by validating in both cases with OK.
6) Select the element to assign the distributed load to it with the Assign/Frame Static Loads/
Point and Uniform menu. Choose FQ for Load Case Name, Forces and GlobalY for
Load Type and Direction and enter the value -60000 in the Uniform Load field.

7) Now select the node at the free end to assign it the load P
with the Assign/Joint Static Loads/Forces menu . Choose FQ for Load Case
Name and enter -90000 for Force GlobalY. Select the node again to assign the
moment M to it. Now choose FM for Load Case Name and enter 60000 for Moment
GlobalZZ. The moment is only visible in the drawing window if it is in 3D mode or in
the XZ plane, change view if necessary.
8) In the Analysis Options window, check only UY and RZ for Available DOF's, also
activate the Dynamic Analysis option and press Set Dynamic Parameters to set
Number of Modes to 1.
9) Start the analysis. Before closing the calculation progress window, we can read:
Period = 0.014265, the second vibration mode is deleted by the SAP since it has
zero periods.
10) The results can be displayed according to the defined combinations or according to the
charges alone.

Note :
Note the definition of loads and their combinations in steps 4 to 6 and compare with the
console.m script. In fact, this remark is generalized for all data declarations, the SAP2000
graphical interface creates a data file which is saved before analyzing the model to be
loaded and read by the software modules. A version in ASCII format can be created with
the File/Export/SAP2000.S2K menu. It is interesting to do this operation and call the

87
Machine Translated by Google

MEF Course – Chapter 5: Beam Element A. Seghir 2005-2014

file for example console.s2k to load it with a text editor like WordPad or NotePad in
order to study it.

5.4 Duty No. 6

Find the expressions of the stiffness and mass matrices of the beam element
(equations: 4.31)
Force F in the case of a triangular distributed load ÿ and

88