Document Text (Pages 11-20) Back to Document

2D Nonlinear Finite Element Analysis of Masonry Bedjoint

by Osvald, Andrej

Page 11


1.2.2 Solution

Finite element method uses algebraic equations in matrix form. We will use a
simple example for explaining the solution procedure.

Fig. 1.2: Finite element of rectangular shape Fig. 1.3: Common nodes

In Fig. 1.2 we see a finite element of rectangular shape. Writing an equation for
each of the nodes describing a displacement of a node, as a function of its co-nodes,
gives us 8 equations from the whole element. It forms a matrix of 8 equations.
Continuing according fundamental laws of mechanics we relate displacements to
stresses. From stresses we obtain strain energy and from this we derive potential
energy and finally from minimum potential energy we obtain a pair of system
equations for the complete element. Every element can be described by the stiffness

� �




The matrix operates with a vector displacement d, which is a displacement d of the
whole element. K is called a stiffness matrix. Displacement and stiffness are related
to load F. Stiffness matrix of the element with 8 nodes can be written:



















Carrying out this process for every element in mesh we get a stiffness matrix for each
element. Next step is to combine all individual matrices together to obtain a stiffness
of a whole system. Any two neighboured elements have got nodes in common (Fig.
1.3). Individual matrices can be combined by a merging technique this process is
called assembly (Fig. 1.4). Using standard procedure we eliminate parts of the matrix.
Rows of the matrix represent a set of simultaneous equations. We solve a first
equation and substitute to remaining ones. At the end, when last matrix is added, we
are left with a solution of a single node its displacement. We are getting backwards,
applying the acquired information, in to the rest of equations, as a key of a system,

Page 12


until every displacement of node is obtained. From these results corresponding
stresses can be calculated.

Fig. 1.4: Merged area

1.2.1 Post-processing

Postprocessor software uses methods for sorting, printing and plotting
solutions of finite element method. There are many possible ways to present the
results: sort element stresses according intensity, check equilibrium, plot deformed
structure, animate dynamic model behaviour. As usual, every computer output
needs to be checked by an engineering judgment, if the solution is physically

1.3 Linearity

Linear analysis is based on linear stress-strain relationship, where Hook’s law
is applied. This method can be used only in elastic part of stress-strain diagram.
Results from different loads are compatible, because of linear superposition validity.
We use linear analysis method mainly because the solutions are easier to find,
computational requirements are small and solutions can be superposed on each
In many cases, linear analysis is not suitable and nonlinear analysis is
necessary when:
designing high performance components
trying to establish the cause of failure
simulating true material behaviour
better understanding of physical phenomena is needed

Page 13


1.3.1 Non-Linear Analysis

The first type is a geometric nonlinear analysis. It is used in cases with large
strains, displacements, and rotations and other. A good example is an airplane wing,
on which large load affects. After every load step a new geometry of a component is
redefined by adding nodal displacements. This guarantees that the component
geometry will have true adapted (deformed) shape for next load step.
The second type is a material nonlinear analysis. Material nonlinearities occur
when the stress-strain or force-displacement law is not linear, or when material
characteristics change with the applied loads. In this case, total load on component is
applied in small steps. In each step, nonlinear stress-strain relationship is
approximated by a linear one. Because our model is controled by a displacement it is
advantageous that the displacement values can be modified after each load step
depending on demands.

1.4 Applications

The Finite element method is very universal tool has got a world-wide usage
in many mechanical engineering disciplines (for example in biomechanical,
aeronautical and automotive industries). The main aim is the design and
development of products. An optimized design means the cost reduction by
minimizing a weight and material consumption. The FEM results allow detailed
visualization of stresses and displacements distribution (also shows weak parts of
model, places with inappropriate bend or twist...). This method is also very practical,
there is no need to physically test the product every time a part or material
characteristic is changed (experimental analysis is expensive and we do not create a
vain models).
Perfect examples of FEM usage are crash-tests (Fig. 1.5). FEM allows designing
a prototype, testing it virtually, redesigning and manufacturing it. The
implementation of FEM has significantly decreased the time to take products from
concept to the production line. In conclusion, benefits of FEM involve increased
accuracy, improved design and better understanding into critical design parameters,
virtual modelling, fewer physical prototypes, faster and less expensive design cycle
and increased productivity.

Fig. 1.5: The FEM application in crash-tests

Page 14


1.5 Limitations of FEM

Even trough, FEM is very universal and effective, and can enable designers to
obtain information about the behaviour of various structures with almost arbitrary
loads, the achieved results must be carefully examined before they can be used.
The most considerable limitation of the finite element method is that the
accuracy of the obtained solution is mostly a function of the mesh resolution. A
problem around point loads and supports is obvious. Large forces on a very small
area mean almost infinite stresses. These places with highly concentrated stress must
be carefully analysed with the use of an adequately refined mesh.
Reaching a solution with finite element method often requires a significant
amount of computer’s and user’s time and it is important not to underestimate user’s
knowledge and experiences. Present packages can solve many sophisticated
problems; there is a strong will to get results without doing the hard work of
thinking. Modern finite element packages are powerful tools that are expandable to
mechanical design and analyses. They become so approachable that a common user
has to avoid making easy mistakes.

1.6 Convergence criteria

The finite element method provides a numerical solution to an overall
problem. It may be expected that the solution converges to the exact solution under
certain circumstances. It can be shown that the displacement formulation of the
method leads to the actual stiffness of the structure. Hence the sequence of
progressively finer meshes is expected to convergence to the exact solution if
presumed element displacement fields satisfy certain criteria.

The displacement field within an element must be continuous. It does not
yield a discontinuous value of function but more like a smooth variation of
function, which do not involve openings, overlap or jumps. This condition can
be satisfied by using a polynomials for the displacement model, for example:

2 3

W C1 C2x C3x C4x ... (1.3)

The displacement polynomial should include a constant term, representing a
rigid body displacement, which should occur at any unrestrained component
when subjected to external load. It also should contain linear terms, which on
differentiation give constant strain terms. Constant strain is the logical
condition as the element size reduces to a point in the limit.
Compatibility of the displacement and its derivatives, up to required order,

must be satisfied across inter-element boundaries. Otherwise, the
displacement solution may result in separated or overlapped inter-element

Page 15


Beside the convergence and compatibility requirements, one of the important
considerations in choosing proper terms is that the element should have no
preferred direction. It means that the displacement shapes will not change
with a change in local coordinate system. This property is known as geometric

1.7 Nonlinear solution methods

It is not possible to solve the nonlinear equation system in one single step as it
could be done for linear equation system. From this reason the solution must be
applying iterative algorithms.
When we use the loading as a number of increments it is called an
incremental-iterative solution method. These methods are also known as path
following algorithms, because the structural behaviour is examined for every step of
the load history and then for the complete load path.
In can be shown in equation 1.4, that total load vector F depends on n scalar
increase factors, the load factors λi. The single load vectors Fi represent the nodal
force values caused by random load collectives of point, line and surface.

F � �1F1 � �2F2 ... � �nF



To simplify following considerations we can restrict to only one parameter
load system. It means, all load factors will be constant and only one will be
incremented. Simplification the load factor only one is remaining – λ and F0 is a
reference load.


� �F



Algorithms that allow us to follow an arbitrary nonlinear path in the loaddisplacement-space
are called following algorithms. Two of them are introduced

1.7.1 Newton-Raphson method

The Newton-Raphson method is one of the most frequently used iteration
scheme. It is based on the incremental solution concept. For a constant current load
level the iterative corrections of displacement are obtained from the out of balance
forces. From this reason the tangential stiffness matrix is used. The two following
attributes describe this method:

It is pure force controlled method. The load level is constant and algorithm is
trying to iteratively find an equilibrium situation for that considered load

Page 16


While we use the exact tangential stiffness matrix in the iteration the algorithm
converges in a quadratic form. It is obvious that the error approximation
decreases with quadratic order from one to the next iteration.

Suppose that initial displacements d0 are known, the first guess of nodal
displacement for load F is obtained by solving a set of linear algebraic equations


KT (0)



F (1.6)



K d )


0) T




is a tangent stiffness matrix calculated for initial displacements.

As the displacements d1 are most probably not accurate, the equilibrium equation 1.8
(matrix representation of a set of nonlinear algebraic equations for unknown nodal
displacements d) is not satisfied

K d ) F



and that means there are unbalanced (or residual) nodal forces



K d ) F



By calculating new tangential stiffness matrix


K d )


1) T




and solving new set of algebraic linear equations

T (

d r

1) 1 1


we will get an improved solution



d1 � �d1

k2 K d ) F 0


If the procedure is iterated until the sufficiently accurate solution is
obtained. The iteration process is shown in Fig. 1.6.

Page 17


Fig. 1.6: The Newton-Raphson method

1.7.2 Modified Newton-Raphson method

The modified Newton-Raphson method is kind of simpler way how to reach
equation equilibrium. The method is based on using the same tangent stiffness
matrix for all iterations. The stiffness is determined once at the beginning of the
iterative process and is used unchanged until the equilibrium is reached. The two
following attributes describe this method:

It is pure force controlled method. The load level is constant and algorithm is

trying to iteratively find equilibrium for that considered load level.
Unlike the standard Newton-Raphson method, the algorithm demonstrates

slower convergence, but there is no need to re-establish the stiffness matrix
after every iteration. On the other hand, number of iterations is usually larger.

There are some occasions, when the tangent stiffness must be modified,
especially when calculation is not able to converge due to zero stiffness of some

Page 18


Fig. 1.7: The Modified Newton-Raphson method

Page 19


2 Theory of elasticity

2.1.1 Elasticity

Almost all structural materials exhibit an elastic response at a certain range of
deformation. For example, if some external forces produce deformation of a structure
and they do not reach the elastic limit, the deformation disappears and no residual
deformation remains.

2.1.2 Stress

In the most simplest case of a prismatic bar loaded by a tensile force uniformly
distributed over the ends (Fig. 2.1), we can assume the internal forces will be also
uniformly distributed over any cross section n-n’.

Fig. 2.1: Prismatic bar

Then the intensity of this distribution (the stress) can be obtained by dividing the
total tensile force F by the area A of the cross section n-n’.
In the general case, when stress is not distributed uniformly over cross section
(Fig. 2.2), to obtain the magnitude of stress appears on small area δA, cut out of the
section area n-n’ at any point 0, we assume that the forces appear across this
elemental area can be reduced to a resultant δF.

Fig. 2.2: General object

Page 20


The limiting value δF/δA ratio presents the magnitude of the stress acting on the
cross section n-n’ at the point 0. Direction of the stress is the limiting direction of the
resultant δF.
In the general case when the stress direction is related to the area δA on which
it acts, the stress is composed of normal and shearing stresses acting in the plane of
the area δA.

2.1.1 Notation for forces and stresses

There are two kinds of external forces which may act on bodies. The first are
forces distributed over the surface of the body (for example the pressure of one body
to another or hydrostatic pressure). The second are forces distributed through the
volume of a body (for example gravitational forces, magnetic forces or if the body is
subjected to acceleration, inertial forces).
We use to mark normal stress with letter σ and shearing stress with letter τ. To
differ on which plane the stresses are acting it is used to add subscripts to these
letters. If we imagine a small cubic element at a point 0 (Fig. 2.3) with sides parallel to
the coordinate axes, the stress acting on the side of this element takes as positive
direction with the axis orientation.

Fig. 2.2: Infinitesimal cube subjected to stresses

The shearing stress has got similar notation except it consists of a two
subscripts. The first one indicates the direction of the normal to the plane and the
second indicates the direction of the component of the stress.

© 2009 All Rights Reserved.