Advanced Finite Element: Professor: Carlos Felippa
Advanced Finite Element: Professor: Carlos Felippa
Advanced Finite Element: Professor: Carlos Felippa
=
on
o
An alternative version of above equation that displays more physical meaning is the force-
balance form:
p() = (, A)
In this form p denotes the configuration-dependent internal forces resisted by the structure
whereas f are the control-dependent external or applied forces, which may also be configuration
dependent.
The residual r is either r p f or r f p, the two versions being equivalent except for sign.
The force equilibrium equations r 0 or f p express the fact that the total potential energy is
stationary with respect to variations of the state vector when the structure is in static
equilibrium.
7
If a total potential energy exists, the decomposition becomes:
=
o
o
=
o
o
Where U and P are the internal and external energy components, respectively, of n = -
2-2-2-Stiffness and Control Matrix:
Varying the vector r with respect to the components of u while keeping fixed yields the
Jacobian matrix K:
K
]
=
or
o
]
This is called the tangent stiffness matrix in structural mechanics applications. The inverse of K,
if it exists, is denoted byF = K
-1
; a notation suggested by the name flexibility matrix used in
linear structural analysis for the reciprocal of the stiffness. If r derives from a potential, both K
and F are symmetric matrices.
Varying the negative of r with respect to while keeping u fixed yields the control matrix as
below:
]
= -
or
oA
]
2-2-3-Staging:
Multiple control parameters are quite common in real-life nonlinear problems. They are the
analog of multiple load conditions in linear problems. But in the linear world, multiple load
conditions can be processed independently because any load combination is readily handled by
superposition. In nonlinear problems, however, control parameters are not varied independently.
It means that we cannot superpose the effect of each load case on the structure.
As an example, consider a bridge which has different load cases such as dead, live, earthquake,
wind, temperature changes and foundation settlement. There are six control parameters
corresponding to each of load cases.
A = {A
1
, A
2
, .. A
6
]
In reference configuration the is the vector of zeros.
A
c]
= {,,,,,]
Suppose the is 1 for dead load and 4 for earthquake:
8
A
1
= A
d
= {A,,,,,] A
q
= {,,, A,,]
In order to get from A
c]
to A
d
we do the first staging (A
1
). Then in order to get from A
d
to
A
q
we do the following stage:
A
2
= {A,,, A,,]
The other load cases are involved in the staging as the above procedure.
2-2-4-One parameter residual:
From this section on we reduce the number of stages to one, because the problem that we are
planning to solve has only one load case. Hence, the load control parameter vector has only one
component, so we only have ={}. The residual vector becomes as follows:
r(, z) =
The corresponding residual-derivative equations are:
r = K - qz
=
or
o
q = -
or
oz
Where K is the tangent stiffness matrix and q is the incremental load vector.
= K
-1
q
This vector is called the incremental velocity vector and is an important component of all
methods.
2-2-5-Separable Residuals and Proportional Loading
The force-balance equivalent for a one-parameter residual equation is:
() = (, z)
If the right hand side, which represents the external force vector, does not depend on the state
parameters u, the system is called separable. So:
() = (z)
Furthermore, if f is linear in the loading is said to be proportional. Obviously q f/ is then
a constant vector.
9
2-3-Critical points (Limit and Bifurcation point):
It is assumed that the structural system is conservative and consequently K is symmetric.
Response points at which K becomes singular are of great interest in the applications because of
their intimate connection to structural stability. These are called critical points, and also
nonregular or singular points. At these points the velocity vector v is not uniquely determined by
q. This means that the structural behavior cannot be controlled by the parameter .
It is convenient to distinguish the following types of critical points:
1. Isolated limit points, at which the tangent to the equilibrium path is unique but normal to
the axis so V becomes infinitely large.
2. Multiple limit points, at which there the tangent lies in the null space of K and is not
unique but still normal to the axis.
3. Isolated bifurcation points also called branch points or branching points, from which two
equilibrium path branches emanate and so there is no unique tangent. The rank deficiency
of K is one.
4. Multiple bifurcation points, from which more than two equilibrium path branches
emanate. The rank deficiency of K is two or greater.
To classify critical points we proceed as follows. Let z be a null right eigenvector of K at a
critical point, that is:
K =
Since K is assumed symmetric
1
K = ; that is, z is also a left null eigenvector. The parametric
differential equation of the equilibrium path is u = z
u
= _
_
_
Although these points generally do not have physical meaning, they can cause special problems
in path following solution procedures because of turnback effects.
In my project they are some combinations of loads which cause turning point in the response
curve.
2-5-kinematic descriptions of geometrically nonlinear finite element
Three kinematic descriptions of geometrically nonlinear finite element analysis are in current use
in programs that solve nonlinear structural problems. They can be distinguished by the choice of
reference configuration.
1. Total Lagrangian description (TL). The reference configuration is seldom or never
changed: often it is kept equal to the base configuration throughout the analysis. Strains
and stresses are measured with respect to this configuration.
2. Updated Lagrangian description (UL). The last target configuration, once reached,
becomes the next reference configuration. Strains and stresses are redefined as soon as
the reference configuration is updated.
12
3. Corotational description (CR). The reference configuration is split. Strains and stresses
are measured from the corotated configuration whereas the base configuration is
maintained as reference for measuring rigid body motions.
In this project the Total Lagrangian Formulation is used as the kinematic description. Hence
our focus of discussion would be TL description only.
2-5-1-Total Lagrangian description (TL)
As explained before, in TL formulation the reference configuration is kept unchanged along
the solution procedure. In this case the Green Lagrange strain tensor and Piola Kirchhoff
stress tensor are used to measure the strains and stresses respectively. It should be noted that
the stress and strains are measured in the reference configuration.
2-5-1-1-Strain measure:
In our problem the Green Lagrange strain tensor has 3 components as the following:
XX =
o
X
oX
+
_
o
X
oX
]
2
+ _
o
oX
]
2
=
o
o
+
_
o
o
]
2
+ _
o
X
o
]
2
X =
_
o
X
o
+
o
oX
] +
_
o
X
oX
o
X
o
+
o
o
o
oX
_ = X
The following figure shows the nonlinear problem in Total Lagrangian description.
The geometrically nonlinear problem in a Lagrangian kinematic description Coordinate systems,
reference and current configurations, and displacements for many (but not all) problems.
13
2-5-1-2-Stress measure:
As long as the strains remain infinitesimal, in spite of having finite deformations, the
hooks constitutive relation for elastic materials is valid and can be used in the
formulations. The constitutive equation is as below:
s
= s
0
+
]
]
Where
and s
0
are stresses in the reference configuration (also called prestresses) and
]
are constant
elastic moduli with
]
=
]
.
]
=
( - v)( + v)
_
- v v
v - v
- v
_ foi plane stiain pioblem in B
]
=
( - v
2
)
_
v
v
- v
_ foi plane stiess pioblem in B
2-5-1-3-Strain Energy Density:
The strain energy density U in the current configuration measured per unit volume of the
reference configuration is as following:
U =
]+
s
It can be observed that the strain energy density is quadratic in the Green-Lagrange strains.
The strain energy in the current configuration is obtained by integrating this energy density
over the reference configuration:
= _ U
0
in B
Above expression forms the basis for deriving finite elements based on the Total Lagrangian
description.
14
2-6-Solution algorithm
Before the mid-1970s, geometrically nonlinear structural problems were usually treated with
purely incremental methods under load control. These methods have the disadvantage of causing
computed solutions to drift away from the equilibrium path. This drift error is step-size
dependent and often accumulates during the analysis, thus requiring a very fine step-size for
accurate analysis. This accumulated drift error can only be practically assessed by re-running the
problem with different step-sizes. Furthermore, traversal of critical points is difficult or even
impossible if one enforces load control, as discussed below.
These shortcomings motivated the development of incremental/ iterative methods where the
increment was followed by equilibrium-correcting iterations that bring the solution back to the
equilibrium path. Introduction of a corrector has the advantage that the drift error is eliminated
and thus, as long as the iterative phase converges, the computed equilibrium path is independent
of the increment step-size.
Geometrically non-linear structures usually reach a maximum load level. At that state the
structure is unable to withstand further load increases until a significant change in geometry
occurs. A load control strategy may be able to detect a limit point but cannot generally traverse
it. Traversal is often desirable to assess whether the structure has residual load carrying
capabilities after what might be a localized instability.
A number of methods of traversing the equilibrium path beyond limit points have been described
in the literature. One can mention the method of artificial springs, methods based on controlling
the load increment with the current stiffness parameter and suppressing iterations around limit
points and the displacement control techniques first introduced by Argyris.
During the past 20 years important improvements have been made by allowing loads and
displacements to be simultaneously varied in each incremental step. The most practically
important instances of these strategies are the hyper plane displacement control developed by
Simons, Bergan and Nygard and the arc-length methods originally proposed by Riks and
Wempner and later refined by Bathe et al, Batoz and Dhatt, Crisfield, Ramm and Riks.
Unfortunately, none of these algorithms are best for all problems. However, the arc-length type
algorithms are generally considered as the most versatile algorithms in terms of the range of
problems they can solve.
15
Among all the different solution algorithms that are available for solving the nonlinear problems,
the Arc length method is selected which includes the corrector step as well.
2-6-1-Conventional Arc length Algorithm:
Arc length method can be summarized in the following figure:
Nested hierarchy in nonlinear solution methods:
stages, increments and iterations (predictor and corrector)
The equilibrium path is a curve in the N +1 dimensional space spanned by the loading parameter
and the N degrees of freedom vector v. An incremental displacement vector v with matching
load increment can be written as the augmented displacement vector A
:
A
= _
A
Az
_
Where prefix is used for the incremental step (or predictor step) of the arc length algorithm
For the equilibrium iterations the increment in the augmented space is written as:
o
= _
o
oz
_
Where prefix denotes corrector changes. The following sections describe the predictor step and
corrector iterations in further detail.
Predictor step:
Advancing from a converged solution, (z,
`
) gives the incremental load solution wq0
from solving (
`
)
q0
= q(z). The incremental displacement step is set to:
A
= Az
q0
= Az _
q0
_
By requiring the step in the load-displacement space to have prescribed length s, the
vector length definition:
16
As =
A
1
A
Gives:
Az = _
As
+
q0
1
q0
The sign for can be determined by assuming that the equilibrium path is smooth, and
that in an advancing solution process the present predictor step must form a positive
vector product with the previous predictor step.
Corrector step:
Among different corrector steps the normal plane iterations is selected for the present
work.
Normal plane iterations:
The normal plane method was developed by Riks and Wempner and is often referred to
as the Riks-Wempner method. Here the name normal plane iterations is adopted due to
its more descriptive character. The corrector steps should reach back to the equilibrium
path as quickly as possible. Because the predictor step direction wq0 is a good
approximation of the tangent direction of the equilibrium path, in the absence of further
information one may assume that the shortest distance to the equilibrium path lies in a
direction orthogonal to the predictor step, as depicted in following picture:
Predictor and corrector step for the normal plane method.
17
Because the normal plane to the predictor step is a hyper plane, subsequent corrector
iteration are forced to lie on that surface. This is the rationale behind the normal plane
iteration method. The iteration displacement vector can be written as:
o
= _
o
oz
_
Note that the symbol is here used to denote a corrector iteration increment, whereas is
used for a predictor step increment. To force ov
)W
r
= i(v
, i) and K(v
)W
q
= q(i) respectively.
The combined solution is:
o
+ oz
q
= _
W
r
_ + oz _
W
q
_
Scaling the load increment oz so that o
+
q0
1
q
Which satisfies the constraint
q0
1
o
=
The following table summarizes the Arc length type algorithm:
18
Predictor step with arc length As:
Set = K(v
)
-1
q(i)
For n= 1 to number of increments
Solve K(v
)W
q0
= q(i) with respect to W
q0
Set = _ +
q0
1
q0
If
q0
1
0
> Az =
As
]
Else Az = -
As
]
Set = AzW
q0
Update z z + Az and
`
`
+ AzW
q
Corrector iterations:
For i=1 to max iteration number for correction
Set W
q
= K(v
)
-1
q(i)
Set W
r
= -K(v
)
-1
i(v
, z)
Set oz = -
w
q0
T
w
r
1+w
q0
T
w
q
Set z = z + oz
Set v
= v
+ (
r+
oz
q)
Until [r(v
= _ X
c
n
=1
= _
c
n
=1
= _
c
n
=1
x
= _ x
c
n
=1
= _ y
c
n
=1
= _
c
n
=1
(X
) and (x
x
= _
x
c
n
=1
= _
c
n
=1
u
x
=x-X u
y
=y-Y
20
3-3-Shape functions:
The following bilinear shape functions are used for interpolation:
1
c
=
1
4
( - )( - p)
2
c
=
1
4
( + )( - p)
3
c
=
1
4
( + )( + p)
4
c
=
1
4
( - )( + p)
3-4-Partial derivative computations:
Partial derivatives of shape functions with respect to the Cartesian coordinates X and Y are
required for the strain and stress calculations. Because shape functions are not directly functions
of X and Y but of the natural coordinates and , the determination of Cartesian partial
derivatives is not trivial. The derivative calculation procedure is presented below for the case of
an arbitrary Iso-parametric quadrilateral element with n nodes.
3-5-Jacobian:
In quadrilateral element derivations we will need the Jacobian of two-dimensional
transformations that connect the differentials of {x, y} to those of {, } and vice-versa. Using
the chain rule:
[ =
(X,)
({,q)
= _
X
{
{
X
q
q
_ = _
[
11
[
12
[
21
[
22
_ , [
-1
=
({,q)
(X,)
= _
{
X
q
X
{
_ =
1
ct(])
_
[
22
-[
12
-[
21
[
11
_
3-6-Shape function derivative:
The shape functions of a quadrilateral element are expressed in terms of the quadrilateral
coordinates and . The derivatives with respect to X and Y are given by the chain rule:
N
i
c
X
=
N
i
c
{
{
X
+
N
i
c
q
q
X
N
i
c
=
N
i
c
{
{
+
N
i
c
q
q
21
3-7-Derivatives of Cartesian coordinates w.r.t. natural coordinates:
oX
o
= (X
n
=1
o
c
o
)
oX
op
= (X
n
=1
o
c
op
)
o
o
= (
n
=1
o
c
o
)
o
op
= (
n
=1
o
c
op
)
3-8-Displacement derivative:
The displacement derivatives are used in the nonlinear formulations, so it is worth deriving the
derivatives of displacements with respect to initial configuration as below:
o
X
c
oX
= (
x
n
=1
o
c
o
)
o
oX
+ (
x
n
=1
o
c
op
)
op
oX
o
X
c
o
= (
x
n
=1
o
c
o
)
o
o
+ (
x
n
=1
o
c
op
)
op
o
o
c
oX
= (
n
=1
o
c
o
)
o
oX
+ (
n
=1
o
c
op
)
op
oX
o
c
o
= (
n
=1
o
c
o
)
o
o
+ (
n
=1
o
c
op
)
op
o
3-9-Strain derivatives with respect to state variables:
oXX
o
X
=
o
2
X
oXo
X
+
o
2
X
oXo
X
o
X
oX
+
o
2
oXo
X
o
oX
oXX
o
=
o
2
X
oXo
+
o
2
X
oXo
o
X
oX
+
o
2
oXo
oX
o
o
=
o
2
oo
+
o
2
oo
o
+
o
2
X
oo
o
X
o
o
o
X
=
o
2
oo
X
+
o
2
oo
X
o
o
+
o
2
X
oo
X
o
X
o
oX
o
X
= .
o
2
X
oo
X
+
o
2
oXo
X
+.
o
2
X
oXo
X
o
X
o
+
o
2
X
oo
X
o
X
oX
+
o
2
oo
X
o
oX
+
o
2
oXo
X
o
o
oX
o
= .
o
2
X
oo
+
o
2
oXo
+.
o
2
X
oXo
o
X
o
+
o
2
X
oo
o
X
oX
+
o
2
oo
oX
+
o
2
oXo
o
All above derivatives are used to measure the GL strain tensor for each element which is used to
determine the internal strain energy as discussed before.
22
The gradient of the strain energy for each element gives the residual and tangent stiffness matrix.
The resulting r and k of each element is assembled into the r and K of the whole structure. The
process will be explained in details in the further sections.
3-10-Selective Reduced Integration: (SRI method)
In Quadrilateral element if the element is very thin, the element becomes stiff in inplane bending,
which is called the shear locking phenomenon. In order to get rid of the shear locking, the
selective reduced integration is used.
In this method the E matrix is split into a volumetric and a shear portion. Thus, the internal
energy is also split into a volumetric and a shear part. 1x1 node Gauss rule is applied on the
shear portion and 2x2 node Gauss rule is applied on the volumetric part when we integrate the
energy over the entire volume of the element.
The E tensor for the plane stress problem is split as below:
]
=
( + v)
_
_ +
v
- v
2
_
_ =
+
x
23
4-Explanation of the FEM code (assumptions and abilities
of the code)
The code consists of three major parts which are:
Preprocessing (Reading input data) :
1. Getting the material properties:
Elastic modulus of material
Poissons ratio
2. Prestress value in 3 directions: It is a vector with initial stress in reference
configuration in
3 Cartesian directions. Dimension is 3x1
3. Number of elements
4. Number of nodes
5. Number of unrestrained degrees of freedom
6. Nodal coordinates: The dimension of this matrix is (number of node x 2)
7. Applied loads on nodes: The dimension of this matrix is (the number of nodes that
have forces on x 3): The first component of each row is the node number and the
other two components of each row correspond to the loads in X and Y directions
respectively.
8. Connectivity between the nodes: that shows that every element is located between
which four nodes. The nodes assignment is counterclockwise. The dimension of
connectivity matrix is (number of element x 4).
9. ID matrix: that assigns an equation number to each degree of freedom of all the
nodes. For example if node one is restrained in X direction and free in Y
direction, the first column of the ID Matrix is |
0
1
]. The dimension of ID is ( 2 x
number of nodes)
10. Getting the data required for iteration procedure:
As : the Arc length
maxitei: the maximum number of iteration in corrector step
Conv_eps1: the convergence criterion for force and displacement residuals
Conv_eps1: the convergence criterion for energy
Inc_num: the number of increments
24
z : initial value for z at the beginning of the iterations
: initial value of u vector at the beginning of the iterations. The length of
u is equal to the number of free degrees of freedom.
Solution:
1. Setting up the load vector: that reduces the degrees of freedom which are
restrained. The length of this vector is equal to the number of free degrees of
freedom or the number of unknowns.
2. Setting up the LM matrix: In this matrix every column corresponds to one
element. The rows of each column give the equation numbers of the corner nodes
in X and Y directions (as set up in ID before). The nodes are assigned in
counterclockwise direction. The dimension of LM matrix is (8 x number of
elements)
3. Determining the reference residual vector and tangent stiffness matrix for the
whole structure considering u and z equal to zero.
Calling Global_r_k function: this function computes the r and K for the
whole structure. The followings are the adopted procedures in this
function to compute the global r and K :
o Loop over the number of elements
Setting up the u vector for each element
Setting up the nodal Cartesian coordinates vector for each
element
Calling r_and_k function: this function calculates the r and k
for each element
o Applying Selective Reduced Integration to measure
the internal energy over the volume of the element.
Calling one_node_Gauss and two_node
_Gauss functions to do all the required
calculations in the element level to measure the
derivatives of displacement, strains, internal
energy and in the end the r
e
a and k
e
.
Establishing the global r vector and global K tensor
25
Assembling the local r and k to the global ones.
o end the loop
Subtracting the z - q from the r_global to get to the assembled
global r vector for the whole system.
4. Starting the iterations as mentioned in the table for arc length type algorithm.
In this part the r and K will be updated due to changes in u vector and i
Post-processing:
1-storing all the incremented z and displacements in some vectors
2-plotting the equilibrium path for each degree of freedom (in z and displacement space)
The following flow chart summarizes different parts of the code in a hierarchical fashion:
26
27
28
5-Verification of the code:
5-1-verification examples:
In order to verify the code, a slender column is modeled with 6 quadrilateral elements and the
buckling load of the system is determined. The column behaves as a clamped bar since it is
restrained at the end against rotation and displacement. The size of each element is .25x3 meter.
So the dimension of the column in XY plane is 9.0 by 0.5 which looks very slender. The
Poisons ratio and the E modulus of the structure are 0.2 and 1.0 KN/m
2
respectively
in all the
calculations. The shape of the structure is shown below:
Geometry of the column
The buckling load of the above-shown column is determined by applying a horizontal load on
nodes 4,8,12 and a very small vertical load on node 12 to perturb the tip of the column.
The following table shows the nodal forces on nodes:
Node 4 8 12
load
X direction P1=-0.25 P2=-0.5 P1=-0.25
Y direction 0 0 P0=0.00001
The resultant response curves for nodes 8 (DOFs 12, 13) and 11(DOFs 17, 18) in both X & Y
directions are plotted below:
29
-0.025 -0.02 -0.015 -0.01 -0.005 0
0
1
2
x 10
-4
state variable u12
l
o
a
d
c
o
n
t
r
o
l
p
a
r
a
m
e
t
e
r
l
a
m
b
d
a
response curve
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
0
1
2
x 10
-4
state variable u13
l
o
a
d
c
o
n
t
r
o
l
p
a
r
a
m
e
t
e
r
l
a
m
b
d
a
response curve
0
0.2
0.4
0.6
0.8
-0.03
-0.02
-0.01
0
0
1
2
3
x 10
-4
state variable u13
state variable u12
l
a
m
b
d
a
30
-0.03 -0.025 -0.02 -0.015 -0.01 -0.005 0
0
1
2
x 10
-4
state variable u17
l
o
a
d
c
o
n
t
r
o
l
p
a
r
a
m
e
t
e
r
l
a
m
b
d
a
response curve
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35
0
1
2
x 10
-4
state variable u18
l
o
a
d
c
o
n
t
r
o
l
p
a
r
a
m
e
t
e
r
l
a
m
b
d
a
response curve
0
0.1
0.2
0.3
0.4
-0.03
-0.02
-0.01
0
0
1
2
3
x 10
-4
state variable u13
state variable u12
l
a
m
b
d
a
31
As can be seen the buckling load in which the response curve has an abrupt change in the
direction is . -
-4
KN. Now we calculate the analytical buckling load for this clamped
column to compare the results:
c
=
n
2
I
2
= . -
-4
K
Where:
E=1 KN/m
2
L=9.0 I=1/12x1x.5
3
=0.0104 m
4
The reasons of the difference between the analytical and numerical solutions can be the effect of
poisons ratio and the finite element approximation. On the other hand, in analytical solution we
neglect the shear effect. Considering these facts, the results of FEM code are acceptable.
5-2-some specific load combinations on the column:
Moreover, Different nodal forces are applied on the structure and the corresponding response
curves are plotted in (z , u) space for nodes 4,8,12.
The following pictures are the response curves for different values of forces on nodes 4,8,12 in X
& Y directions.
Node 4 8 12
load
X direction P1=-0.25 P2=-0.5 P1=-0.25
Y direction 0 0 P0=0.0
-25 -20 -15 -10 -5 0
-0.1
0
0.1
0.2
0.3
0.4
0.5
0.6
state variable u12
l
o
a
d
c
o
n
t
r
o
l
p
a
r
a
m
e
t
e
r
l
a
m
b
d
a
response curve
32
-6 -5 -4 -3 -2 -1 0 1
x 10
-12
-0.1
0
0.1
0.2
0.3
0.4
0.5
0.6
state variable u13
l
o
a
d
c
o
n
t
r
o
l
p
a
r
a
m
e
t
e
r
l
a
m
b
d
a
response curve
-16 -14 -12 -10 -8 -6 -4 -2 0
-0.1
0
0.1
0.2
0.3
0.4
0.5
0.6
state variable u17
l
o
a
d
c
o
n
t
r
o
l
p
a
r
a
m
e
t
e
r
l
a
m
b
d
a
response curve
-0.04 -0.03 -0.02 -0.01 0 0.01 0.02 0.03
-0.1
0
0.1
0.2
0.3
0.4
0.5
0.6
state variable u18
l
o
a
d
c
o
n
t
r
o
l
p
a
r
a
m
e
t
e
r
l
a
m
b
d
a
response curve
33
Node 4 8 12
load
X direction P1=0.25 P2=0.5 P1=0.25
Y direction 0 0 P0=0.0
0 5 10 15 20 25
0
2
4
6
8
10
12
state variable u12
l
o
a
d
c
o
n
t
r
o
l
p
a
r
a
m
e
t
e
r
l
a
m
b
d
a
response curve
-0.5 0 0.5 1 1.5 2 2.5 3 3.5 4
x 10
-15
0
2
4
6
8
10
12
state variable u13
l
o
a
d
c
o
n
t
r
o
l
p
a
r
a
m
e
t
e
r
l
a
m
b
d
a
response curve
34
0 5 10 15
0
2
4
6
8
10
12
state variable u17
l
o
a
d
c
o
n
t
r
o
l
p
a
r
a
m
e
t
e
r
l
a
m
b
d
a
response curve
-0.35 -0.3 -0.25 -0.2 -0.15 -0.1 -0.05 0
0
2
4
6
8
10
12
state variable u18
l
o
a
d
c
o
n
t
r
o
l
p
a
r
a
m
e
t
e
r
l
a
m
b
d
a
response curve
35
Node 4 8 12
load
X direction P1=0.0 P2=0.0 P1=0.0
Y direction 0 0 P0=-1.0
-10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
state variable u12
l
o
a
d
c
o
n
t
r
o
l
p
a
r
a
m
e
t
e
r
l
a
m
b
d
a
response curve
-20 -18 -16 -14 -12 -10 -8 -6 -4 -2 0
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
state variable u13
l
o
a
d
c
o
n
t
r
o
l
p
a
r
a
m
e
t
e
r
l
a
m
b
d
a
response curve
36
-6 -5 -4 -3 -2 -1 0 1
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
state variable u17
l
o
a
d
c
o
n
t
r
o
l
p
a
r
a
m
e
t
e
r
l
a
m
b
d
a
response curve
-14 -12 -10 -8 -6 -4 -2 0
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
state variable u18
l
o
a
d
c
o
n
t
r
o
l
p
a
r
a
m
e
t
e
r
l
a
m
b
d
a
response curve
37
Node 4 8 12
load
X direction P1=-1.0 P2=-2.0 P1=-1.0
Y direction 0 0 P0=1.0
-7 -6 -5 -4 -3 -2 -1 0
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
x 10
-4
state variable u12
l
o
a
d
c
o
n
t
r
o
l
p
a
r
a
m
e
t
e
r
l
a
m
b
d
a
response curve
0 1 2 3 4 5 6 7 8
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
x 10
-4
state variable u13
l
o
a
d
c
o
n
t
r
o
l
p
a
r
a
m
e
t
e
r
l
a
m
b
d
a
response curve
38
0
2
4
6
8
-8
-6
-4
-2
0
0
1
2
x 10
-4
state variable u13
state variable u12
l
a
m
b
d
a
-3 -2.5 -2 -1.5 -1 -0.5 0
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
x 10
-4
state variable u17
l
o
a
d
c
o
n
t
r
o
l
p
a
r
a
m
e
t
e
r
l
a
m
b
d
a
response curve
-3 -2.5 -2 -1.5 -1 -0.5 0
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
x 10
-4
state variable u17
l
o
a
d
c
o
n
t
r
o
l
p
a
r
a
m
e
t
e
r
l
a
m
b
d
a
response curve
39
Node 4 8 12
load
X direction P1=1.0 P2=2.0 P1=1.0
Y direction 0 0 P0=1.0
0
1
2
3
4
-3
-2
-1
0
0
1
2
x 10
-4
state variable u13
state variable u12
l
a
m
b
d
a
-0.5 0 0.5 1 1.5 2 2.5
0
0.01
0.02
0.03
0.04
0.05
0.06
state variable u12
l
o
a
d
c
o
n
t
r
o
l
p
a
r
a
m
e
t
e
r
l
a
m
b
d
a
response curve
40
0 0.5 1 1.5 2 2.5 3
0
0.01
0.02
0.03
0.04
0.05
0.06
state variable u13
l
o
a
d
c
o
n
t
r
o
l
p
a
r
a
m
e
t
e
r
l
a
m
b
d
a
response curve
-0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8
0
0.01
0.02
0.03
0.04
0.05
0.06
state variable u17
l
o
a
d
c
o
n
t
r
o
l
p
a
r
a
m
e
t
e
r
l
a
m
b
d
a
response curve
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
0
0.01
0.02
0.03
0.04
0.05
0.06
state variable u18
l
o
a
d
c
o
n
t
r
o
l
p
a
r
a
m
e
t
e
r
l
a
m
b
d
a
response curve
41
As we can see in the responses, by changing the load combinations on different nodes we can
capture limit points, buckling point (bifurcations) and the turning points. The slender column is a
good exercise for capturing all the important points that may occur in nonlinear analysis.
The arc length type solution algorithm with hyperplane method for corrective stage gives
reasonable results for this problem.
5-2- Assumption of the code:
The following are the assumptions made in developing the nonlinear quadrilateral Iso-P element:
1. The material is elastic with finite deformation but small strains. Hence the hook
constitutive equation is still valid for this element. As we know this code is not good for
modeling the elastic membranes with large strains such as rubbers or polymers.
2. The loads are assumed to be conservative. The code is not capable of solving systems
with non-conservative loads such as wind or hydrostatic loads which change direction
along with the deformation of the structure.
3. The material is assumed to be compressible. This code is not suitable to model the
material with Poissons ratio close to 0.5.
4. The structure has 2D plane stress behavior. For plane strain problems we should add the
corresponding E tensor for plane strain behavior to the code.
5. Only one stage is defined in this code. The only force that can be applied is the nodal
forces.
5-3- Directions for further research:
The following issues could be investigated and added to improve the code:
Extending the element to 3D
Adding drilling degrees of freedom to each node. This has been done by
Haugen in 1994 but the kinematic description was Corotational which
assumed a linearized strain in the formulation. Since the resulting element
will no longer be compatible due to the drilling degrees of freedom, the
large strain and finite rotations must be investigated carefully.
42
Changing the constitutive relation in such a way that the element can
model the finite deformations with large strains such as rubbers and
polymers.
Adding plasticity to the code.( material nonlinearity)
6-Appendix A:
The code Matlab file