A Continuum Based Solid Shell Element Based on EAS and ANS

Master's Degree Thesis ISRN: BTH-AMT-EX--2015/D11--SE

Supervisors: Maciej Wysocki and Mohammad Rouhi, SICOMP AB Sharon Kao-Walter, BTH

.

Department of Mechanical Engineering Blekinge Institute of Technology

Karlskrona, Sweden 2015

Waleed Ahmad Mirza

A
Continuum
Based
Solid

Shell
Element
Based
on
EAS

and
ANS

Waleed
Ahmad
Mirza

Department
of
Mechanical
Engineering Blekinge
Institute
of
Technology

Karlskrona,
Sweden. 2015

Thesis
submitted
for
completion
of
Master
of
Science
in
Mechanical
Engineering
with
emphasis
on
Structural
Mechanics
at
the
Department
of
Mechanical
Engineering,
Blekinge
Institute
of
Technology,
Karlskrona,
Sweden.

Abstract

This
work
is
a
stepping
stone
towards
developing
higher
order
shell
element
for
simulating
composite
manufacturing
procedure.
In
this
study,
a
continuum
approach
suitable
for
combined
material
and
geometrically
nonlinear
analysis
for
an
eight
node
solid
shell
element
SS8
is
explained.
The
formulation
of
SS8
comprises
two
ingredients
to
alleviate
undesirable
locking
effects:
1)
Assumed
Natural
Strain
concept,
which
has
proven
to
alleviate
the
curvature
thickness
and
transverse
shear
locking
problems.
2)
Enhanced
Assumed
Strain,
which
adds
enhanced
degrees
of
freedom
to
improve
the
in-plane
response
of
the
element
and
the
curvature
thickness
locking
problem.
This
formulation
has
been
extended
to
represent
geometric
and
material
non-linearity
using
Total
Lagrangian
approach.
Finally,
finite
strain
formulation
has
been
verified
by
numerical
examples.
Results
when
compared
to
continuum
shell
element
in
ABAQUS
show
a
reasonable
agreement
with
a
relative
error
of
less
than
2%.

Keywords

Acknowledgement

This research work is carried out in Swerea SICOMP AB and Blekinge Insti-tute of Technology, Karlskrona, Sweden, under the supervision of Dr. Maciej Wysocki, Dr. Mohammad Rouhi and Prof. Sharon Kao-Walter.

I am grateful to Dr. Maciej Wysocki who provided me the opportunity to con-duct this research. Moreover my sincere appreciation goes to Dr. Mohammad Rouhi for helping me with technical difficulties and to all my colleagues at SICOMP for providing me a conducive atmosphere to carry out this research and above all for their valuable support and advice. Last but not least, my father, mother and girlfriend, who owes my deepest gratitude and love for being an endless source of motivation for me.

Gothenburg, May 2015

Contents

1 Notation  1.1 Variables . . .  1.2 Abbreviations . . .  1.3 Superscripts . . .  1.4 Operators . . .  2 Introduction  2.1 Motivation . . .  2.2 Solid Shell Elements – State of the art review . . .  2.3 Outline of thesis and proposed methodology . . . 

3 Theory of Finite Strain Analysis 

3.1 Total Lagrangian Formulation . . .  3.1.1 Integrating EAS and TL Formulation . . . 

4 SS8 Formulation- Small Strain kinematics 

4.1 Linear Element Formulation . . .  4.2 Thickness Strain . . .  4.3 Transverse Shear Strain . . .  4.4 In-plane Strain Response . . .  4.5 Transformation to Orthonormal Coordinate System . . .  4.6 Enhanced Strain . . .  4.7 Element Stiffness Matrix . . . 

5 SS8 Formulation- Finite Strain kinematics 

5.4 Algorithm . . . 

6 Numerical Analysis 

6.1 Benchmarking . . .  6.2 Cantilever beam bending . . .  6.2.1 Results . . .  6.2.2 Discussion . . .  6.3 Analysis
of
thin
Eeam
(
L
_t =30,
300
and
3000)
.
.
.
.
.
.
.
.
.
.
.
.

6.3.1 Conclusion . . .  6.4 Cook’s Membrane Problem . . .  6.4.1 Results . . .  6.4.2 Discussion . . .  6.5 Patch Test . . .  6.5.1 Result . . .  6.5.2 Discussion . . .  6.6 Distortion test . . .  6.6.1 Results . . .  6.6.2 Discussion . . . 

7 Conclusion and Recommendations 

7.1 Conclusion . . .  7.2 Learning outcomes . . .  7.3 Future Work . . .   8 References 9 Appendix 

Notation

Variables

Benh BNL BOMH B b Cijrs c D E t+t o E eij F F

Enhanced
strain
displacement
tensor
Non-linear
strain
displacement
tensor
Non-linear
strain
displacement
matrix
Linear
strain
displacement
tensor

Left
green
Cauchy
tensor Component
of
Fourth
order
material
tensor

Right
green
Cauchy
stress Energy
dissipation Elasticity
modulus Green
Lagrangian
strain
tensor
at
time
t
+
t

fenh
Internal
force

due
to
enhanced
strain
matrix
and
stress
state

Internal stress state tensor

Shear modulus

ith column of Jacobian matrix

ith column of spatial basis vector

Fourth order Identity Tensor

Determinant of deformation gradient

Linear stiffness tensor

Non-linear stiffness tensor

Spatial velocity gradient

Tensor of order 3x24 containing shape functions

Spatial normal to a surface

Shape function at node I fJ O U
_G Gi
gi
I₄
J KL
KNL
l
N
n NI

N_Xk Derivative of Shape function at node k with respect to spatial coordinate
X

P Stress power

R External load

S Second Piola-Kirchhoff stress

Sv Second Piola-Kirchhoff stress in vector formulation

t True traction vector

u Nodal displacement

ui,j Displacement component i with respect to component j

Vo Undeformed volume

v Poisson ratio

X, Y, Z Global coordinates( material coordinates)

x, y, z Deformed coordinates

ξ, η, ζ Natural coordinates

α Enhanced degree of freedom

ηij Non-linear strain tensor

εij Components of Green Lagrangian strain tensor

Ψ Free energy

τ Cauchy stress

σxx Longitudinal stress along x-axis

σyy Longitudinal stress along y-axis

σxy In plan sheer stress

Assumed Natural Deviatoric Strain

Assumed Natural Strain

Abbreviations

ANDES

ANS

SS8 Eight Node Solid Shell Element Total Lagrangian Updated Lagrangian TL UL

Superscripts

c Orthonormal Coordinate system

h Discrete

n Natural coordinate system

Operators

δ Variation operator

But
these
advantages
come
at
cost
of
many
undesirable
phenomenas,
popu-larly
referred
as
locking
effects
[1].
Locking
occurs
when
a
shell
element
is
modelled
like
a
solid
brick
element
using
displacement
interpolation
which
tends
to
‘lock
out’
realistic
displacement
of
element
response
by
activating
extraneous
strains
that
require
much
higher
energy
input
than
strains
of
the
realistic
mode.
These
locking
behaviours
prevent
solid
elements
to
be
used
for
shell
like
structures.

The
goal
of
this
project
is
to
formulate
an
eight
node
solid
shell
element

Introduction

 They
are
computationally
effective
and
reliable
in
terms
of
capturing
in plane
and
transverse
response
compared
to
brick
elements.
For
instance,
in application
such
as
sheet
metal
deforming
membrane
stretching,
bending
and shearing
are
very
dominant.
In
such
scenario,
solid
shell
elements
can
be successfully
implemented
to
capture
such
intricate
responses.

 Unlike
planer
shell
elements,
using
solid
shell
elements
boundaries
of
three dimensional
structures
can
be
modelled
without
introducing
any
kinematic assumptions.

review

Solid
shell
elements
are
quite
similar
to
brick
elements
in
nodal
formulation
and
exhibit
only
translational
degree
of
freedom.
Along
with
several
advan-tages,
these
elements
come
with
a
number
of
disadvantages
resulting
from
smaller
span
of
thickness
dimension
compared
to
the
lateral
dimensions
[1].
Solid
shell
concept
was
developed
to
overcome
well
known
degenerate
concept.
Using
nodes
at
upper
and
lower
surface
and
displacement
degree
of
freedom,
stresses
and
strains
in
longitudinal,
transverse
and
thickness
direction
are
cal-based
on
continuum
approach
using
isotropic
hyper-elastic
constitutive
ma-terial
model.
This
formulation
is
coupled
with
techniques
such
as
Assumed
Natural
Strain
(ANS)
and
Enhanced
Assumed
Strain
(EAS)
aimed
at
re-moving
locking
effects.
First
a
solid
shell
element
is
formulated
for
small
strain,
linear
kinematics
and
later
extended
to
finite
strain,
nonlinear
mate-rial
kinematics
using
Lagrangian
approach.
At
the
end
case
studies
have
been
proposed
to
verify
the
formulation.

2.1
Motivation

This
work
is
part
of
a
bigger
project
aimed
at
developing
a
simulation
tool
for
modelling
composite
manufacturing
process
of
a
wide
range
of
popular
infusion
techniques
using
higher
order
solid
shell
elements.
Higher
order
shell
elements
have
extensive
application
in
fields
like
porous
media
theory
where
pressure
field,
a
derivative
of
displacement,
is
approximated
with
linear
variation
[].
As
a
consequence
of
linear
pressure
variation,
the
displacement
is
rendered
quadratic.
This
quadratic
displacement
trend
can
be
very
well
captured
with
a
20
node
solid
shell
element.
Moreover,
in
other
applications
of
structural
mechanics,
higher
order
shell
elements
are
better
capable
of
capturing
curved
geometries.

The
current
project
is
a
stepping
stone
towards
developing
higher
order
solid
shell
elements.
The
project
involves
formulating
a
3D-shell
element
compris-ing
8
nodes
based
on
hyper-elastic
material
model.
The
approach
developed
for
formulating
eight
node
solid
shell
element
will
be
used
for
modelling
higher
order
solid
shell
element
in
the
follow
up
project.
Such
higher
order
element
will
be
implemented
in
simulating
composite
manufacturing
procedures
like
vacuum
infusion.

culated accurately. Along with several advantages, these elements come with a number of disadvantages resulting from smaller span of thickness dimension compared to the lateral dimensions. In literature several techniques such as ANDES, ANS and EAS [1], [4] have been applied to overcome deficiencies such as extra undesirable energy modes and locking phenomenas. A review on the locking behaviour is as follow.

1- Transverse shear locking / trapezoidal locking is the inability of an element to exhibit zero shear strain when subjected to pure bending. This defect is owed to the formulation of standard strain displacement matrix us-ing displacement interpolation, as a result of which shear strain terms arise in strain displacement matrix because of in plane strain terms. This idea can be illustrated in the following example.

Figure

2.1:

Deformed
and
un-deformed
state
of
element
under
pure
bending.

Consider
the
4
node
quadrilateral
element
as
illustrated
in
figure
2.1.
The
element
is
subjected
to
pure
bending
which
is
supposed
to
render
zero
shear
strain.
Let’s
assume
the
displacement
vector
in
the
current
deformation
is
as
follows:

(u₁, v₁, u₂, v₂, u₃, v₃, u₄, v₄) = (1, 0,−1, 0, 1, 0, −1, 0) (2.1) These nodal displacements will trigger shear strain owing to the following strain displacement formulation.

⎧ ⎪ ⎨ ⎪ ⎩ εxx εyy εxy ⎫ ⎪ ⎬ ⎪ ⎭= 1/4

−(1−y) 0 (1−y) 0 (1+y) 0 −(1+y) 0 0 −(1−x) 0 −(1+x) 0 (1+x) 0 (1−x) −(1−x) −(1−y) −(1+x) (1−y) (1+x) (1+y) (1−x) −(1+y)

This
problem
can
be
minimized
by
using
reduced
integration.
Shear
strains
are
evaluated
at
x=0
and
y=0,
which
results
in
zero
shear
strains
as
can
be
seen
in
equation
(2.2).

2- Volumetric
locking
arise
in
problems
comprising
incompressible
or
nearly compressible
constitutive
material
models
where
1oisson
ratio
is
equal
to
0.5. This
value
renders
the
material
matrix
equal
to
infinity
as
shown
in
equation (2.3-2.6)
leading
to
a
very
high
value
of
stress.
There
are
several
methods
to avoid
this
behaviour
one
of
which
is
reduced
integration
or
using
constraints such
as
given
in
equation
(2.6).

(2.3) C = E (1 + v)(1− 2v) ⎡ ⎢ ⎣ 1− v v 0 v 1− v 0 0 0 1− 2v ⎤ ⎥ ⎦ (2.4)

Equation (2.4) becomes as follows at v=0.5.

C =∞ ⎡ ⎢ ⎣ 0.5 0.5 0 0.5 0.5 0 0 0 0 ⎤ ⎥ ⎦ (2.5) exx=−eyy (2.6) ⎡ ⎢ ⎣ σxx σyy σxy ⎤ ⎥ ⎦= E (1 + v)(1− 2v) ⎡ ⎢ ⎣ 1− v v 0 v 1− v 0 0 0 1− 2v ⎤ ⎥ ⎦ ⎧ ⎪ ⎨ ⎪ ⎩ εxx εyy 2εxy ⎫ ⎪ ⎬ ⎪ ⎭

The
constraint
limits
the
value
of
stress
on
each
integration
points
by
can-celling
the
denominator
term
of
(1-2v)
in
equation
(2.4).

(2.7)

 Curvature/Trapezoidal
locking
occurs
when
element
edges
in
thickness direction
are
not
perpendicular
to
the
mid
plane.
This
type
of
situation
arises when
curved
geometries
are
modelled
with
solid
shell
element.
This
defect can
be
reduced
by
using
Assumed
Natural
Strain
concept>@

 Membrane
locking
happens
when
the
element
is
subjected
to
in-plane longitudinal
or
transverse
(shear)
loads
and
the
low
order
shape
functions
are not
capable
of
modelling
the
physical
behaviour
of
the
element.
ANDES
and ANS
approach
can
used
to
alleviate
this
behaviour
[1].

Assumed
Natural
Strain
(ANS)
concept
was
first
introduced
in
1978
by
Park
and
Stanley
[]
for
doubly
curved
thin
shell.
This
technique
is
effective
in
alleviating
shear
locking.
A
similar
technique
called
Mixed
Interpolation
Tensoral
Components
(MITC)
was
developed
by
Bathe
et
al.
[4].
The
As-sumed
Natural
Deviatoric
Strain
concept
(ANDES)
presented
by
Felippa
and
Militello,
represents
a
combination
of
free
formulation
of
Bergan
[]
<>
and
has
been
extensively
used
by
researchers
for
alleviating
membrane
strain.
AN-DES
and
MITC,
though
very
effective,
has
not
been
used
much
in
the
past.
EAS
technique
originated
from
variational
framework
presented
by
Simo
and
Rifai
[15]
which
ultimately
evolved
to
EAS
variant.
EAS
is
mainly
used
to
counter
Poisson
thickness
locking.
All
locking
alleviating
techniques
have
been
extensively
used
in
both
small
and
finite
strain
analysis.

2.3

Outline of thesis and proposed

method-ology

To accomplish the goals of this project, following methodology has been fol-lowed.

1- Literature review to comprehend theory and equations constituting finite strain kinematics and Enhanced Assumed Strain concept (Chapter 3).

2- SS8 formulation for small strain linear kinematic using Assumed Natural Strain (ANS) and Enhanced Assumed Strain (EAS) concepts (Chapter 4).

3- Finite strain formulation of SS8 element using Total Lagrangian approach (Chapter 5).

4- Implementation and verification of formulated element in the proposed case studies (Chapter 6).

o Sij

Theory
of
Finite
Strain

Analysis

This
chapter
details
theory
and
constitutive
equations
behind
the
adopted
fi-nite
strain
methodology.
There
are
two
popular
methodologies
to
approach
a
finite
strain
problem:
Total
Lagrangian
(TL)
and
Updated
Lagrangian
(UL)
formulation.
According
to
Bathe
et
al.
[4],
in
Total
Lagrange
formulation
(TL),
the
reference
configuration
is
the
undeformed
or
material
configuration
as
opposed
to
Updated
Lagrangian
(UL)
formulation
where
the
reference
con-figuration
is
set
as
the
current
concon-figuration
from
the
last
converged
incre- ment.
In
TL
formulation,
integrals
are
solved
over
the
undeformed
configura-tion
and
Second
Piola-Kirchhoff
stress
and
Green
Lagrangian
strain
are
used
as
stress
and
strain
measures.
Whereas
in
UL,
formulation
integrals
are
solved
over
current
configuration
and
if
the
increments
are
small,
Cauchy
stress
and
Rate
of
Deformational
tensor
are
employed
as
stress
strain
measures.
The
downside
of
Updated
Lagrangian
approach
includes
more
computation
since
reference
state,
volume
and
stress
orientation
are
updated
in
every
incre-ment.
But
unlike
Total
Lagrangian
formulation,
Updated
Lagrangian
does
not
contain
any
initial
displacement
effect
in
strain
measure.
In
this
study,
Total
Lagrangian
formulation
has
been
followed
coupled
with
EAS
approach
to
model
ILQLWHstrain
response.

.1
Total
Lagrangian
Formulation

Consider
energetically
conjugate
pair
of
stress
and
strain
at
time
t
+
t
re-ferred
to
reference
state
at
time
t=0
denoted
as
t+t
and
t+to εij
.
Principle



Vo

t+t

o Sijδt+t₀ εijdVo=t+tR (3.1)

In an incremental approach solution at time t is known (for example t₀Sij, t

0ui,j etc.). Therefore stresses and strains are decomposed as follows: t+t o Sij
= t P
Sij
+P
Sij (3.2) t+t o εij
= t P
εij
+Pεij (3.3)

In
incremental
procedure,
the
stress
tP
Sij
and
strain
tP
εij
stated
at
time
t
are

known
while
increments
i.e.
P
εij
,
P
Sij
are
unknown.
Therefore,
these
unknown

increments
are
required
to
be
estimated
at
any
given
time
step.
Defining
the
Green
Lagrangian
Strain
Tensor
at
time
t
and
t
+
t
in
following
equations. t Pεij
= 1 2( t

Pui,j
+tP
uj,i
+tP
uk,itPuk,j
) (3.4)

and t+t P εij = 1 2( t P+tui,j+ t P+tuj,i+ t P+tuk,i tP+t uk,j) (3.5)

Now substituting equation (3.4) and (3.5) in equation (3.3) yields the follow-ing. Pεij
= 1 2(0ui,j
+0
uj,i
+ t P
uk,i
tP
uk,j
)
+
1 2uk,iPuk,j (3.6)

Now,
linear
strain
increment
Peij
and
non-linear
strain
increment
Pηij
can
be

defined
as
follows:

Variation in incremental Green Lagrangian Strain yields the following.

δoεij
=
δPηij
+
δPeij (3.10)

Now equation of virtual work becomes,

 Vo oSijδoεijdVo+  Vo t PSij
δ0ηij
dVo
=
t+tR
−  Vo t

PSij
δPeij
dVo









(3.11)

The above mentioned equation is a non-linear function of the unknown dis-placement increment. Therefore, equation (3.11) is linearized to obtain the following form.

t

oK



U
=
t+tR
−to
f (3.12)

Equation (3.12) is an approximation of equation (3.11) obtained by neglecting all higher order terms in displacement increment. A detailed description of linearization procedure can be found in [36]. Here final form is mentioned.

 V_ooCijrsoersδoeijdVo+  V_o t oSijδoηijdVo =t+t R−  V_o t

PSij
δPεij
dVo



(3.13)

3.1.1

Integrating
EAS
and
TL
Formulation

The
EAS
formulation
was
introduced
in
1990
by
J.C.
Simo
et
al.
[15].
En-hanced
Assumed
Strain
Approach
is
capable
of
treating
volumetric,
thickness
and
shear
locking.
This
formulation
enhances
the
in
plane
strain
and
thick-ness
strain
response
by
adding
extra
degree
of
freedom
in
strain
displacement
matrix.
In
this
section
we
will
layout
the
mathematical
formulation
governing
this
approach.
In
EAS,
Linear
strain
displacement
matrix
B
is
enhanced
by
adding
few
extra
columns,
denoted
by
B_enh.

Bnew=B Benh



(3.18)

Matrix
Benh
is
formulated
by
meeting
the
following
conditions
[2]

Benh

should
be
linearly
independent
from
B
as
given
in
following
equation.
Ignoring
the
condition
will
render
matrix
singularity
and
will
give
a
non-unique
solution.

Benh∩ B = φ (3.19)

2- Second
Piola-Kirchoff
Ttress
and
enhanced
strain
displacement
matrix should
be
orthonormal
i.e.



Vo

B_enhT SvdVo= 0 (3.20)

Since Second Piola-Kirchoff stress is constant, therefore equation (3.20) is written as:



Vo

BenhdVo = 0 (3.21)

Now the new form of equation (3.17) is as follow.

 VoB T_CBdV o+_VoBnlgT S9BnlgdVo _VoBTCBenhdVo  VoB T enhCBdVo  VoB T enhCBenhdVo   u α  =  Rt+t 0  −   VoB

Sv
dVo  V_oBenhSvdVo  (3.22)

Where α are the extra degree of freedom added as a result of enhancing the strain displacement matrix. In short form equation (3.22) can be written as:

_t oKL+toKNL GT G A   u α  =  Rt+t 0  − _t ofint t ofenh  (3.23)

Lastly, these extra degree of freedoms are condensed out using static conden-sation method. Finally displacement degrees of freedom are left as given in the following equation:

(toKL+toKNL− GTA−1G)u =t+to R−tofint+ GTA−1tofenh (3.24) t

o

In
linear
formulation
terms
such
as
toKNL
and
fint
will
be
ignored

and
Gollowing
equation
JT
MFGU

(KL− GTA−1G)u = R− GTA−1fenh (3.25)

SS8
Formulation-
Small

Strain
kinematics

In
chapter
3,
continuum
mechanics
behind
Total
Lagrangian
method
and
Enhanced
Assumed
Strain
(EAS)
method
were
discussed.
In
this
chapter,
finite
element
implementation
of
these
concepts
is
explained.
Techniques
such
as
Enhanced
Assumed
Strain
and
Assumed
Natural
Strain
will
be
used
to
model
small
and
finite
strain
response
of
SS8
solid
shell
element.
First,
moving
on
from
equation
(3.25)
small
strain
response
of
SS8
element
is
modelled.
Later
using
TL
formulation
(as
discussed
in
chapter
3),
this
formulation
will
be
extended
to
finite
strain,
non-linear
material
analysis.
MATLAB
code
based
on
this
formulation
will
be
verified
with
case
studies
in
chapter
6.
Small
strain
formulation
of
SS8
element
proceeds
as
follow.

.1
Linear
Element
Formulation

Figure

4.1:

Geometry
of
SS8
Solid
Shell
Element
[28]. X =  X Y Z T (Xi, i = 1, 2, 3) (4.1) ξ =  ξ η ζ T (ξi, i = 1, 2, 3) (4.2) x =  x y z T (xi, i = 1, 2, 3) (4.3) J =G₁ G₂ G₃= ⎡ ⎢ ⎣ ∂X ∂ξ ∂X∂η ∂X∂ζ ∂Y ∂ξ ∂Y∂η ∂Y∂ζ ∂Z ∂ξ ∂Z∂η ∂Z∂ζ ⎤ ⎥ ⎦ (4.4)

Table

4.1:

Coordinate
System.

Coordinate system Base vectors Coordinates Environment

Orthonormal r₁ , r₂ , r₃ x,y,z Local

Natural ξ₁, ξ₂, ξ₃ ξ, ζ, η Local

Orthonormal Coordinate System is constructed at element centre ξ = 0, η = 0, ζ = 0 by the following system of equations.

P₁ = ⎡ ⎢ ⎣ ∂X ∂ξ ∂Y ∂ξ ∂Z ∂ξ ⎤ ⎥ ⎦ (4.5) P₃ = ⎡ ⎢ ⎣ ∂X ∂ξ ∂Y ∂ξ ∂Z ∂ξ ⎤ ⎥ ⎦× ⎡ ⎢ ⎣ ∂X ∂η ∂Y ∂η ∂Z ∂η ⎤ ⎥ ⎦ (4.6) P₂ = P₁× P₃ (4.7) ri = 1 PiPi, i = 1, 2, 3 (4.8)

Orthonormal coordinate system is used to compute parameters such as Second Piola-Kirchhoff stress (in nonlinear analysis) and Green Lagrangian strain (in both linear and nonlinear analysis). One advantage of constructing such sys-tem is that in nonlinear analysis local stress and strain components between iterations can be added without any transformation. If these parameters were in global form then they would have to be rotated between the increments. Lastly, the natural coordinate system is used to build strain interpolation and enhanced strain interpolation matrices.

In order to estimate coordinates of points across the element, isoparametric formulation is used. Isoparametric formulation relates global coordinates of points within the element to nodal coordinates using shape functions.

Xh= 8  I=1 N₁(ξ, η, ζ)XI (4.9) Where, NI(ξ1, ξ2, ξ3) = 1 8(1 + ξ1ξ)(1 + η1η)(1 + ζ1ζ) (4.10) Superscript h is used to denote finite element approximation and XI

XI=XI, YI, ZI (4.11)

A similar isoparametric formulation is also used to interpolate displacements in terms of shape function.

Uh =

8



I=1

N₁(ξ, η, ζ)UI (4.12)

Similar to XI, UI represents displacement of nodal point I of SS8 element,

shown as follow:

UI = [UxI, UyI, UzI] (4.13)

Furthermore, N tensor is defined as follows:

N = [N₁I, N₂I, N₃I, N₄I, N₅I, N₆I, N₇I, N₈I] (4.14)

Where
I
is
3x3
identity
matrix
and
N
is
a
tensor
of
order
3x24.
Once
coor-dinate
systems
and
kinematic
parameters
have
been
defined,
strain
displace-ment
matrix
of
SS8
eledisplace-ment
will
be
estimated.

4.2
Thickness
Strain

Thickness
strain
response
is
formulated
using
Green
Lagrangian
strain
coupled
with
Assumed
Natural
Strain
approach.

Ec
=
E_ijn
Gi

⊗

Gj (4.15) Where E_ijn represents the entries of Green Lagrangian Strain tensor in natural coordinate system represented by equation (4.16) and Gi represents columns of inverse Jacobian matrix.

En= (E₁₁n, E₂₂n, En₃₃, 2E₁₂n, 2E₁₃n, 2E₃₂n)T (4.16)

In general E_ijn is represented as,

Note that the nonlinear term ∂U_∂ξ

i

∂U

∂ξ_j in equation (4.17) is ignored in linear

analysis. Now in order to calculate thickness strain, i, j= 3 will be inserted to the equation and equation (4.9) and (4.12) will be substituted in equation (4.17). Final form of thickness Green Lagrangian Strain proceeds as below.

E₃₃n = GT₃N₃U (4.18)

Where N₃ represents the derivative of tensor N shape function vector with respect to ζ. Now the derivative of E₃₃n with respect to U yields out of the plane strain displacement matrix B₃₃.

B₃₃= GT₃N₃ (4.19)

Now
using
Assumed
Natural
Strain
method,
B33
is
interpolated
at
four

col-location
points
(shown
in
Figure
4.2)
and
the
four
values
are
interpolated
linearly
to
obtain
the
assumed
strain
field. B₃₃ANS = 4  L=1 1 4(1 + ξLξ)(1 + η1η)(G L 3)TN3L (4.20)

Where
G₃L,N₃L
are
calculated
at
points
A1,
A2,
A3,
A4.
Schwarze
et
al.
[]
showed
that
constructing
element
with
this
formulation
alleviates
thickness curvature
locking.

(4.21)

Figure

4.3:

Natural
coordinates
of
collocation
points
used
in
transverse
strain
interpolation
[35].

4.3
Transverse
Shear
Strain

Transverse
shear
strains
are
also
formulated
in
similar
lines
as
thickness strain
E₃₃
n using
Green
Lagrangian
Strain
tensor
and
ANS
concepts.
Here
only
final
form
of
derivation
is
presented.  B₁₃ANS B₂₃ANS  = 1 2  (1− η)B₁₃B + (1− η)B₁₃D (1− ξ)B₂₃A + (1− ξ)BC₂₃  Cardoso
et
al.
[29]
have
shown
that
this
approach
for
constructing
the
trans-verse
shear
strain
field
is
effective
in
alleviating
transverse
shear
locking
in
the
element
behaviour.
It
is
to
be
noted
that
in
[1],
Mustafa
has
used
full
integration
by
utilizing
4
collocation
points
for
each
transverse
shear
strain.
However,
similar
to
the
current
work,
Quak
[27]
also
used
two
collocation
points
for
enhancing
the
computational
efficiency
of
the
algorithm.

4.4
In-plane
Strain
Response

⎛ ⎜ ⎝ B₁₁ B₂₂ B₁₂ ⎞ ⎟ ⎠= ⎛ ⎜ ⎝ GT₁N₁ GT₂N₂ GT₂N₁+ GT₁N₂ ⎞ ⎟ ⎠ (4.22)

4.5

Transformation to Orthonormal

Coor-dinate
System

Until
now
strain
displacement
matrix
is
estimated
in
natural
coordinate
sys-tem.
Now
these
strain
components
are
transformed
into
orthonormal
coordi-nate
system.
Following
series
of
steps
are
taken
to
undergo
this
transforma-tion.

Ec = F En (4.23)

Where
F
is
the
transformation
matrix
converting
Green
Lagrangian
strain
from
natural
coordinates
to
orthonormal
coordinates
[2].

F
=

(4.24)

As
seen,
F
matrix
is
written
in
terms
of
components
of
T
matrix.
Note
that
F
is
formulated
with
respect
to
the
components
of
strain
mentioned
in
order
given
by
equation
(4.16).
T
matrix
is
evaluated
with
the
following
equation
[2].

T = [r₁ r₂ r₃TJ−1T] (4.25)

matrix
in
orthonormal
coordinate
system
is
also
arranged
in
order
of
strain
vector
given
in
equation
(4.16).

4.6
Enhanced
Strain

Now
in
order
to
counter
membrane,
shear
and
thickness
locking
SKnomenas,
columns
of
strain
displacement
matrix
are
enhanced.
Strains
are
enhanced
ac-cording
to
equation
(3.22).
The
assumed
natural
thickness
strain
calculated
in
equation
(4.20)
is
increased
by
one
extra
degree
of
freedom
to
alleviate
thickness
locking.
In
an
attempt
to
improve
the
in-plane
response
of
element
[1],
five
extra
degrees
of
freedom
are
used
(represented
in
first
two
and
fourth
row
of
enhanced
strain
displacement
matrix)
in
equation
(4.26).

Benh(ξ, η, ζ) = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ξ 0 0 0 ξη 0 0 η 0 0 −ξη 0 0 0 0 0 0 ζ 0 0 ξ η ξ2− η2 0 0 0 0 0 0 0 0 0 0 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (4.26) Benh(ξ, η, ζ) = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ξ ζ 0 0 ξη ξζ 0 0 0 0 η ζ 0 0 ξη ηζ ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (4.27)

M1= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ξ 0 0 0 0 0 0 0 0 0 0 η 0 0 0 0 0 0 0 0 0 0 0 ζ 0 0 0 0 0 0 0 0 0 ξ η 0 0 0 0 ξζ 0 0 0 0 0 ξ ζ 0 0 0 0 0 0 0 0 0 0 η ζ 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (4.28) M₂= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 0 0 0 ξη ξζ 0 0 0 0 0 0 0 0 0 0 ξη ηζ 0 0 0 0 0 0 0 0 0 0 ξζ ηζ 0 0 0 0 0 0 0 0 0 0 ξη ηζ 0 0 0 0 0 0 0 0 0 0 ξη ξζ 0 0 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (4.29) M₃= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 0 0 ξηζ 0 0 0 0 0 0 0 0 0 0 ξηζ 0 0 0 0 ηζ 0 0 0 0 0 ξηζ 0 0 0 0 ξη 0 0 0 0 0 ξηζ 0 0 0 0 ξζ 0 0 0 0 0 ξηζ 0 0 0 0 ηζ 0 0 0 0 0 ξηζ ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (4.30)

Now combining M1 M2 and M3 together in a matrix as follow:

Benh(ξ, η, ζ) =



M₁ M₂ M₃ (4.31) In SS8 formulation, Benhmatrix used in equation (4.31) is implemented. The

enhanced strain displacement matrix is converted to orthonormal coordinate system using the following.

(4.32) Benhc =
FPBenh

Where
Fo
is
the
value
of
transformation
matrix
at
ξ
=
0,
η
=
0,
ζ
=
0.

4.7


Element
Stiffness
Matrix

(5.1)

(5.2)

Where
xI
represent
the
coordinate
of
node
I
and
is
a
matrix
of
3x8
order.
Spatial
coordinates
are
now
used
to
define
spatial
basis
vectors
(also
known
as
columns
of
spatial
Jacobian
vector)
as
follow.

SS8
Formulation-
Finite

Strain
kinematics

In
this
chapter,
kinematic
formulation
of
SS8
solid
shell
element
is
ex-tended
to
finite
strain,
material
non-linear
analysis.
As
explained
in
chapter
3,
Total
Lagrangian
together
with
EAS
concept
will
be
used
in
formulating
finite
strain
analysis.
In
last
section
of
this
chapter,
an
incremental
procedure
for
implementation
of
SS8
element
in
FEM
is
introduced.

5.1
Geometric
and
Kinematic
Formulation

Relation between Gi and gi is as follows:

gi = Gi+ ∂U

∂ξi (5.3)

(5.4) Equation
(5.3)
can
also
be
written
as:

gi
=
FGi

Where
F
is
the
deformation
gradient
not
the
transformation
matrix
F
given
in
chapter
4.

5.1.1

Strain
and
Enhanced
Strain
Parametrization

Green
Lagrangian
Strain
tensor
components
are
formulated
in
local
coordi-nate
system
as
follow.

t PE

D

=
t

PEij (5.5)

(5.6)

Strain displacement matrix of equation (5.6) is as follow.

Bij = GTi Nj+ GTjNi+ UTNiTNj+ UTNjTNi (5.7)

For the sake of simplicity from now onward time increments will not be men-tioned in superscripts. Equation (5.7) can be written in simplified form by exploiting equation (5.4).

Bij = giTNj+ gjTNi (5.8)
_G

i

_⊗

_G

j

O

(5.9)

Enhanced
Strain
displacement
matrix
is
formulated
in
exactly
the
same
way
as
mentioned
in
the
last
chapter.
Transformations
are
carried
out
in
accor-dance
with
equation
(4.23)
and
equation
(4.24)
to
obtain
strain
displacement
matrices
Bc
in
orthonormal
coordinate
system.

5.1.2

Internal
load
vectors

The
internal
force
vectors
are
also
calculated
in
the
local
orthonormal
coor-dinate
system
as
follow.

fint=  Vo (Bc)TSvdVo (5.10) fenh=  V_o(B c enh)TSvdVo (5.11) KL=  Vo (Bc)TCBcdVo (5.12)

Stress stiffening matrix is computed with series of equations and transforma-tions given below.

B = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ g₁TN1 g₂TN₂ 4 l=1 1 4(1 + ξ1Lξ1)(1 + ξL2ξ2)(g3L)TN3L g₂TN₁+ g₁TN₂ 0.5[(1− n)B₁₃B + (1 + n)B₁₃D] 0.5[(1− ξ)B₂₃A + (1 + ξ)B₂₃C] ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ Where
Sv

is
Second
Piola-Kirchhoff
stress
in
vector
formation
containing
en-tries
of
S
tensor.

5.2



Stiffness
Computation

As
stated
earlier
total
stiffness
of
SS8
element
consists
of
linear
stiffness
KL

which
comes
from
the
strain
displacement
matrices
and
KNL
which
comes

KNL=  Vo (Bnlg)TS9BnlgdVo (5.13) Where, S₉ = ⎡ ⎢ ⎣ S 0 0 0 S 0 0 0 S ⎤ ⎥ ⎦ (5.14)

Note that S is a 3x3 tensor of Second Piola-Kirchhoff stress.

 Bnlg_k  = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ N_Xk 0 0 N_Yk 0 0 N_Zk 0 0 0 N_Xk 0 0 N_Yk 0 0 N_Zk 0 0 0 N_Xk 0 0 N_Yk 0 0 N_Zk ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ k = 1· · · 8 (5.15)

All entries denoted by B_knlg are assembled in a 9x24 matrix as follow.

Bnlg_{ = B}nlg

1 , B2nlg, B3nlg, B4nlg, B5nlg, B6nlg, B7nlg, B8nlg (5.16)

In
next
section,
equations
used
to
estimate
Second
Piola-Kirchhoff
stress
ten-sor
and
fourth
order
material
tensor
C
are
discussed.

5.3
Constitutive
Equations
of

Hyperelas-ticity

It
is
assumed
that
the
free
energy
density
Ψ(per
unit
volume
of
the
refer-ence
configuration)
acts
as
potential
function
for
stress
due
to
total
(elastic)
deformation.
In
order
to
assure
objectivity
of
the
response
because
of
super
imposed
rigid
body
motions,
free
energy
Ψ
is
related
to
the
deformation
via
right
Cauchy–Green
deformation
tensor
c
as
follow
[30]:

The
constitutive
relations
for
hyper
elastic
response
are
obtained
for
zero
en-ergy
dissipation
D
at
constant
temperature.
This
is
expressed
as:

D = P − Ψ = 0 (5.18)

Where P is the stress power. Now with these fundamental assumptions, Neo-Hookean model based on isochoric-volumetric split is employed to estimate Cauchy stress tensor.

τ = τiso+ τvol (5.19)

Where,

τiso = GJ−23_(b−I1

31) (5.20)

Where invariant I₁ is defined as I₁ =trb = 1:b

τvol= KJ (J − 1)1 (5.21) Now Second Piola-Kirchhoff stress is calculated by push forward transforma-tion equatransforma-tion as follows.

5.4
Algorithm

(5.26) u(k+1)
=
0

1- Computation
for
internal
and
external
forces
is
as
follows: fJOU(k+1)
=  V_o[B c (uk+1)]TS(uk+1)dVo (5.27) R =  VoN T
_ρ ob
dVP
+  ∂Vo NTtdA (5.28)

2- Calculate out of balance force.

Re= [G(k+1)]TA− (5.29) Gk+1
=  Vo (B c enh)TCBc(uk+1)dVo (5.30) K =  Vo [Bc(uk+1)]T
C
Bc(uL_)dVo
+  V_o (Bnlg) T
_S 9BnlgdVo (5.31) fFOI
=  Vo (B_enhc )T

4

WdVo (5.32) A =  Vo (Bcenh)TCBenhc dVo (5.33)

4- Assembling and solving these equations gives following results:

K_T(k+1) = ηelm_ e=1 K(k+1)− [G]TA−1G (5.34) e
1fFOI
+
R
−
fJOU(k+1) If
R5=
Ae η =1 elm_R e
≤
tolerance
−→
END

3- Compute
components
of
tangential
stiffness
matrix.

Algorithm
to
implement
the
proposed
formulation
in
a
finite
element
program
is
as
follows:

Numerical
Analysis

In
this
chapter,
a
number
of
case
studies
will
be
solved
using
the
formulated
MATLAB
code
and
results
will
be
benchmarked
against
results
obtained
from
8
node
continuum
shell
element
in
ABAQUS.
The
purpose
of
benchmarking
is
to
verify
the
formulated
MATLAB
code
and
response
of
SS8
element
under
given
natural
and
essential
boundary
conditions.
For
this
purpose,
several
case
studies
have
been
chosen
from
the
literature.
A
summary
of
these
case
studies
is
mentioned
in
table
6.1.

Table

6.1:

Description
of
proposed
case
studies

Case
Studies. Challenge/Purpose.

Cantilever
beam
under
tip
loading. To
verify
transverse
response
of
SS8
element. Analysis
of
thin
beam. To
verify
SS8
element
is
free
_{of
shear
locking.} Cook’s
membrane
problem
[1],
[35]. 1- To
verify
in
plane
re-_{sponse
of
SS8
element.}
2-To
verify
SS8
element
is
free of
membrane
locking.

Distortion
test
[35]. SS8
element’s
response
sensitivity
to
element
distortion.

Figure

6.1:

Material
Model
in
ABAQUS.

6.2
Cantilever
beam
bending

This
case
study
involves
a
three
dimensional
cantilever
plate
subjected
to
tip
loading.
Consider
a
simply
supported
cantilever
beam
subjected
to
static
load
of
F = 1E5
N
applied
vertically
at
its
tip,
as
shown
in
Fig.
6.2.
Let
the
length
of
the
plate
be /
=
0.3 m,
its
cross
section
a
rectangle
of
height
H
and
thickness
t,
thus
with
area
A
=
H.t
(L
>>
A
and
H/t
>
5).
The
beam
is
made
of
an
isotropic,
hyper-elastic
material
with

E=200.E9
1P

and
Poisson
ratio
v
=
0.3
and
its
shear
modulus
by
G
_{= 2(v}E +1).

6.1
Benchmarking

Figure

6.3:

Comparison
of
results
from
finite
strain
models.

2- Secondly,
values
of
in-plane
strain
from
MATLAB
and
ABAQUS
are compared
(Table
6.2).
Here
only
the
in-plane
strains
are
compared
because ABAQUS
is
only
equipped
to
give
in-plane
strain
parameters,
whereas
MAT-LAB
code
is
also
capable
to
give
transverse
and
out
of
the
plane
strains
in

Figure

6.2:

Cantilever
beam
case
Ttudy

6.2.1

Results

In
this
section,
results
of
finite
strain
analysis
for
given
beam
bending
problem
are
displayed.

Strain Parameter Results from ABAQUS Results from MATLAB Exx 9.09E-2 8.98E-2 Eyy 1.603E-2 1.58E-2 Exy 7.278E-3 7.211E-3 output.
Deformed
and
undeformed
states
of
beam
obtained
from
both
models
are
illustrated
in
fig.
6.4.

Table

6.2:

Comparison
of
strain
values
of
finite
strain
model
from
ABAQUS
and
MATLAB.

Figure

6.4:

Deformed
and
undeformed
beam
geometry
obtained
from
ABAQUS
and
MATLAB

6.2.2

Discussion

The
vertical
tip
displacement
values
obtained
from
ABAQUS
and
MATLAB
are
0.26338
m
and
0.261
m
respectively.
The
relative
percentage
difference
error
is
0.9%.
When
strain
values
are
compared
as
shown
in
table
6.2,
a
relative
percentage
error
of
0.92%
is
obtained.
This
study
proves
that
SS8
element
exhibits
an
accurate
transverse
response
when
compared
with
8
node
continuum
shell
element
of
ABAQUS.





"OBMZTJT
PG
UIJO
CFBN
-U



BOE


increasing the stiffness of the beam. In this study, SS8 element is used to analyse beam with _{T hickness}Length = 30, 300 and 3000 and compared results with an
8
node
brick
element
in
ABAQUS
to
prove
that
EAS
modified
element
is
free
of
shear
locking.
To
compensate
the
reduction
of
beam’s
stiffness
in
this
analysis,
tip
load
is
reduced
by
1000
when
thickness
is
decreased
by
a
factor
of
10. Results
of
this
comparison
are
outlined
in
table
6.3.

Table

6.3:

Analysis
of
thin
beam.

Length to thickness ratio. wmax from contin-uum shell element [m]. wmax from SS8 [m]. wmax from brick ele-ment [m]. Percentage error be-tween column 3 and 4. 30(thick beam)

2.14E-1 2.123E-1 2.047E-1 3.58%

300 2.242E-1 2.198E-1 1.783E-1 20.47%

3000(thin beam)

2.26E-1 2.19E-1 1.274E-1 43.46%

6.3.1

Conclusion

It
can
be
seen
that
brick
elements
show
a
sizeable
increase
in
beam
stiffness
in
thin
beam
case
as
opposed
to
SS8
element,
whose
bending
stiffness
does
not
change
much
with
increasing
L/t
ratio.
For
brick
element
the
value
of
stiffness
increases
from
3.58%
to
43.46%
as
L/t
ratio
is
increased.
The
reason
of
this
increase
in
stiffness
is
that
8
node
brick
elements
does
not
involve
modifications
such
as
ANS
and
EAS
in
its
formulation.
Hence
in
case
of
thin
structures,
they
cannot
counter
locking
phenomenas.
Whereas,
formulation
of
SS8
element
involves
modifications
that
counter
locking.
Hence,
this
analysis
indicates
that
SS8
element
is
free
of
shear
locking.

6.4
Cook’s
Membrane
Problem

Figure

6.5:

Cook’s
membrane
study.

6.4.1

Results

For
the
given
shear
load
and
geometric
properties,
analysis
has
been
carried
out
for
finite
strain
models
of
ABAQUS
and
MATLAB.
A
maximum
tip
displacement
is
measured
to
be
1.504
units
and
1.491
units
obtained
from
ABAQUS
and
MATLAB
respectively.
A
graph
representing
mesh
conver-gence
for
the
given
problem
is
presented
in
fig.
6.6.
Behavior
of
SS8
element
for
values
of
Poisson
ratio
0.4,
0.495
and
0.4995
are
plotted
and
displayed
in
fig.
6.7.
Deformed
state
of
the
skewed
membrane
obtained
from
ABAQUS
and
MATLAB
are
shown
in
the
figure
6.8.

Figure

6.6:

Mesh
Convergence
test
for
Cook’s
membrane
problem
at
1oisson
ratio=1/3.

Figure

6.8:

Results
of
Cook’s
membrane
problem
carried
out
at
very
fine
mesh.

u1= (x + y 2)× 10 −3_{, u} 2 = (y +x 2)× 10 −3_{, u} 3 = 0.0 (6.1)

Mesh description is shown in fig. 6.9.

Figure

6.9:

Element
mesh
for
patch
test
[35].

6.4.2

Discussion

This
reasonable
agreement
in
value
of
maximum
vertical
displacement
shows
the
membrane
response
is
well
captured
by
SS8
element
and
there
is
no
membrane
locking.
It
can
be
seen
in
GJg.
6.7
that
SS8
element
gives
an
accurate
membrane
response
until
v
=
0.495.
But
for
0.4995,
SS8
element
fails
to
converge.
Hence
SS8
element
cannot
model
response
of
material
with
Poisson
ratio
greater
than
0.495.
Moreover,
mesh
convergence
test
shows
that
8
node
continuum
shell
element
converges
faster
than
SS8
element.

6.5
Patch
Test

6.5.1

Result

As
a
result
of
subjected
boundary
conditions,
the
plate
undergoes
deformation
and
nodal
displacements
outlined
in
table
6.4
are
obtained.

Table

6.4:

Nodal
degrees
of
freedom
in
patch
test. Nodal

dis-placements

u v u v u v u v

Analytical solution

5E-5 4E-5 1.95E-4 1.2E-4 2.E-4 1.6E-4 1.2E-4 1.2E-4 Solution obtained from SS8 elements 4.99E-5 3.9E-5 1.94E-4 1.19E-1 1.99E-4 1.599E-4 1.199E-4 1.2E-4 Moreover,
a
constant
stress
state
is
achieved
in
all
five
elements
BT
GPMMPX. S₁₁= S₂₂= 1499, S₁₂= 399.8 (6.2)

6.5.2

Discussion

Nodal displacements show a reasonable agreement with analytical formulation [1]. Moreover, constant stress DFURVV DOO HOHPHQWV shows that SS8 element has passed the patch test. This should be seen as no surprise since in chapter 3 it was assumed that volume integration of enhanced strain displacement matrix is zero, which is necessary condition for passing the patch test.

6.6
Distortion
test

)LJXUH0HVK'LVWRUWLRQDQDO\VLV

Figure

6.10:

Illustration

of

mesh

distortion

study.

6.6.1

Results

6.6.2

Discussion

1- Beam subjected to end loading: In this study cantilever beam is sub-jected to a tip load and its response in terms of vertical displacement and in-plane strain is compared with ABAQUS. Comparison shows a small rel-ative error of less than 1%, verifying accurate transverse response of SS8 element.

2- Analysis of Thin Beam: This study is carried out to assess SS8 ele-ment’s response under shear locking. Results of SS8 elements are compared with 8 node hexahedral brick element in ABAQUS. It is observed that brick

Conclusion
and

Recommendations

In
this
chapter
the
results
of
this
research
project
will
be
summarized
and
potential
follow
up
avenues
of
research
will
be
explored.

7.1
Conclusion

element shows a stiff behaviour as length to thickness ratio is increased. When length to thickness ratio is increased from 30 to 3000, stiffness of brick ele-ment increases from 3.58% to 43.46%. On contrary to brick eleele-ment, SS8 element retains its stiffness with increase in length to thickness ratio. This study shows that that SS8 element doesn’t show any extraneous locking in thin structures when subjected to transverse loading.

 Cook’s
Membrane
Problem:
This
study
constitutes
a
skewed
beam subjected
to
an
in-plane
shear
load.
Accurate
membrane
response
of
the skewed
beam
estimated
by
SS8
element
proves
that
SS8
element
is
free
of membrane
locking.
Whereas,
1oisson
locking
test
revealed
that
SS8
element does
not
converge
for
materials
with
1oisson
ratio
greater
than
0.495,
thus failing
the
1oisson
locking
test.
Mesh
convergence
of
Cook’s
membrane
study shows
that
SS8
element
exhibits
slow
convergence
when
compared
with
8
node continuum
shell
element
in
ABAQUS

In above mentioned studies, error ranging from 0.5%-1.8% are observed. These errors can be attributed to the difference of element model used in ABAQUS and SS8 formulation, explained as follow.

1- In ABAQUS, reduced integration has been carried out coupled with hour glass mode remedies in order to remove shear locking from the continuum shell element. While in SS8 element, ANS and EAS approach has been car-ried out for this purpose.

2- In ABAQUS, the continuum shell element is capable to exhibit 6 degrees of freedom (3 translational DOF and 3 rotational DOF) per node while SS8

4- Patch Test: Patch test has been carried out to check quality of

formu-lated SS8 element. Patch check is used to verify if given element is able to fulfil the conditions of compatibility, convergence and completeness in finite element analysis. In the study, it is shown that a plate with irregular mesh of SS8 element exhibits constant stress field under given displacement boundary conditions, which is a necessary condition to pass the patch test.

5- Mesh Distortion Test: In this study a cantilever beam is subjected to

7.3
Future
Work

For
future
work,
following
recommendations
are
proposed.

1- The
methodology
proposed
for
developing
SS8
element
can
be
used
to develop
a
higher
order
20
node
element,
as
the
one
shown
in
fig.7.1.
It
is
learnt from
the
literature
review
that
by
far
such
element
has
not
been
developed mainly
because
of
lack
of
its
application.
Mostly
in
structural
analysis
using more
elements
with
less
number
of
nodes
is
favoured
than
using
less
higher order
elements.
In
opinion
of
our
research
team,
20
node
thick
shell
element
is very
feasible
for
simulating
composite
manufacturing
procedures
like
vacuum infusion
where
linear
pressure
field
[34]
is
only
possible
if
displacement
field is
quadratic
or
higher
order.
Developing
such
an
element
will
take
more computation
but
essentially
the
same
methodology.

is
only
capable
of
exhibiting
3
degrees
of
freedom
per
node.
This
indicates
a
difference
in
respective
formulations
of
SS8
and
continuum
shell
element.

3- The
nature
of
model
used
in
ABAQUS
for
8
node
continuum
shell
element could
not
be
explored
indepth
because
of
the
lack
of
a
detailed
description
in ABAQUS
theory
manual
[33].

7.2
Learning
outcomes

In
this
project
following
learning
outcomes
and
goals
have
been
achieved. 1- An in-house code consisting of 8 node solid shell element is developed, which is capable to capture finite strain and material non-linearity.

2- A great deal of learning in fields of non-linear continuum mechanics and finite element analysis is achieved.

3- In-depth experience of structured programming in MATLAB.

Numerical
Example Challenge

Stretching
of
a
cylinder
with
free
ends.

Geometric
non-linear
be-haviour
of
the
element
against
bending
and
mem-brane
modes.

Bending
patch
test. Ability
of
the
element
to
re-produce
a
constant
bending. Twisted
beam. Performance
of
wrapped

structures
when
curvature
locking
is
dominant.

Compression
of
Clamped
square
plate.

Studying
element’s
response
under
buckling.

Figure

7.1:

20
node
solid
shell
element

2-
In
this
study,
SS8
element
is
only
verified
to
be
free
of
shear
and
membrane
locking.
Membrane
and
shear
locking
are
most
dominant
of
all
lockings
while
less
dominant
lockings
such
as
curvature,
thickness
and
curvature
lockings
are
yet
to
be
verified.
Some
case
studies
that
can
further
verify
SS8
element,
are
proposed
in
table
7.1.

References

1. Mohammad Reza Mostafa, (2011), A geometric nonlinear solid-shell el-ement based on ANDES, ANS and EAS concepts. Unpublished doctoral dissertation, University of Colorado, Boulder.

2. Voyiadjis, George Z., and Dimitrios Karamanlidis, eds., (2013) Advances in the Theory of Plates and Shells. Elsevier.

3. Arnold, Douglas N., Alexandre L. Madureira, and Sheng Zhang, (2002), On the range of applicability of the Reissner–Mindlin and Kirchho–Love plate bending models, Journal of elasticity and the physical science of solids, Vol. 67, No. 3, Pg.171-185.

4. Bucalem, M. L., and K. J. Bathe, (1997), Finite element analysis of shell structures. Archives of Computational Methods in Engineering Vol.4, No.1, Pg.3-61.

5. Wall, Wolfgang A., Michael Gee, and Ekkehard Ramm, (2000), The chal-lenge of a three–dimensional shell formulation the conditioning problem, In Proceedings of ECCM, Vol. 99.

6. Alves de Sousa, R. J., R. M. Natal Jorge, R. A. Fontes Valente, and J. M. A. César de Sá., (2003), A new volumetric and shear locking-free 3D enhanced strain element, Engineering Computations, Vol.20, No. 7, Pg. 896-925.

7. De Sousa, RJ Alves, J. W. Yoon, R. P. R. Cardoso, RA Fontes Valente, and J. J. Grácio. (2007), On the use of a reduced enhanced solid-shell (RESS) ele-ment for sheet forming simulations, International Journal of Plasticity, Vol.23, No. 3, Pg. 490-515.

el-ements, International Journal for Numerical Methods in Engineering, Vol.36, No. 8, Pg.1311-1337.

9. Park, K. C., and G. M. Stanley, (1986), A curved Co shell element based on assumed natural-coordinate strains, Journal of applied Mechanics, Vol.53, No. 2, Pg.278-290.

10. Cardoso, Rui PR, and Jeong Whan Yoon, (2005), One point quadrature shell element with through-thickness stretch, Computer Methods in Applied Mechanics and Engineering, Vol.194, No. 9, Pg.1161-1199.

11. Cardoso, Rui PR, Jeong Whan Yoon, and Robertt A. Fontes Valente,(2006), A new approach to reduce membrane and transverse shear locking for one-point quadrature shell elements: linear formulation, International Journal for Numerical Methods in Engineering, Vol.66, No. 2, Pg.214-249.

12. Cardoso, Rui PR, Jeong-Whan Yoon, José J. Grácio, Frédéric Barlat, and José MA César de Sá, (2002), Development of a one point quadrature shell element for nonlinear applications with contact and anisotropy, Computer methods in applied mechanics and engineering, Vol.191, No. 45, Pg.5177-5206.

13. Cardoso, Rui PR, and Jeong-Whan Yoon, (2007), One point quadrature shell elements: a study on convergence and patch tests, Computational Me-chanics, Vol. 40, No. 5, Pg.871-883.

14. Cardoso, Rui PR, Jeong Whan Yoon, Made Mahardika, S. Choudhry, R. J. Alves de Sousa, and R. A. Fontes Valente, (2008), Enhanced assumed strain (EAS) and assumed natural strain (ANS) methods for one-point quadrature solid-shell elements, International Journal for Numerical Methods in Engi-neering, Vol. 75, No. 2, Pg.156-187.

15. Simo, Juan C., and M. S. Rifai, (1990), A class of mixed assumed strain methods and the method of incompatible modes, International Journal for Numerical Methods in Engineering, Vol.29, No. 8, Pg.1595-1638.

17. M. Schwarze and S. Reese. (2009), A reduced integration solid-shell finite element based on the EAS and the ANS concept-Geometrically linear prob-lems, International Journal for Numerical Methods in Engineering, Vol. 80, No. 10, Pg. 1322-1355.

18. F. Abed-Meraim and A. Combescure (2002), SHB8PS-a new adaptative, assumed-strain continuum mechanics shell element for impact analysis, Com-puters and Structures, Vol. 80, Pg. 791-803.

19. Farid Abed-Meraim and Alain Combescure (2009), An improved assumed strain solid-shell element formulation with physical stabilization for geomet-ric non-linear applications and elastic-plastic stability analysis, International Journal for Numerical Methods in Engineering, Vol. 80, No. 13, Pg. 1640-1686.

20. Coultate, John K., Colin HJ Fox, Stewart McWilliam, and Alan R. Malvern, (2008), Application of optimal and robust design methods to a MEMS accelerometer, Sensors and Actuators A: Physical, Vol.142, No. 1, Pg.88-96.

21. Criseld, M-A, (1981), A fast incremental/iterative solution procedure that handles “snap-through”, Computers Structures Vol.13, No. 1, Pg.55-62.

22. Criseld, M. A., and G. F. Moita, (1980), A COROTATIONAL FORMU-LATION FOR 2D CONTINUA INCLUDING INCOMPATIBLE MODES, International Journal for Numerical Methods in Engineering, Vol.39, No. 15, Pg. 2619-2633.

23. P. G. Bergan, (1980), Finite elements based on energy orthogonal func-tions, International Journal for Numerical Methods in Engineering, Vol.15, No.10, Pg.1541-1555.

24. Bergan, P. G., and M. K. Nygård, (1984), Finite elements with increased freedom in choosing shape functions, International Journal for Numerical Methods in Engineering, Vol.20, No. 4, Pg. 643-663.

26. De Sá, César, M. A. José, Renato M. Natal Jorge, Robertt A. Fontes Valente Pedro M.Almeida Areias, (2002), Development of shear locking free shell elements using an enhanced assumed strain formulation, International Journal for Numerical Methods in Engineering, Vol.53, No. 7, Pg.1721-1750.

27. Quak, Wouter, (2007), A solid-shell element for use in sheet deformation processes and the EAS method.

28. Schwarze, Marco, and Stefanie Reese, (2009), A reduced integration solid-shell finite element based on the EAS and the ANS concept—Geometrically linear problems, International Journal for Numerical Methods in Engineering, Vol.80, No. 10, Pg.1322-1355.

29. Cardoso, Rui PR, Jeong Whan Yoon, Made Mahardika, S. Choudhry, R. J. Alves de Sousa, and R. A. Fontes Valente, (2008), Enhanced assumed strain (EAS) and assumed natural strain (ANS) methods for onepoint quadrature solidshell elements, International Journal for Numerical Methods in Engineer-ing, Vol. 75, No. 2, Pg. 156-187.

30. Ragnar Larsen, (1997), MULTIPLICATIVE FINITE STRAIN HYPER ELASTO—PLASTICITY — Basic Theory and its Relation to Numerical Methodology.

31. Bathe, K. J., Ramm, E., Wilson, E. L. (1975), Finite element formulations for large deformation dynamic analysis, International Journal for Numerical Methods in Engineering, Vol.9, No. 2, Pg.353-386.

32. RD Cook, DS. Malkus, and ME. Plesha. Concepts and applications of finite element analysis. John Wiley and Sons, New York, 3rd edition, 1989.

33. Systèmes, D. (2012), ABAQUS 6.12 Theory manual, Dassault Systèmes Simulia Corp., Providence, Rhode Island.

34. Mohammad S. Rouhi (2007), Infusion Modeling Using Two Phase Porous Media Theory, Independent Master’s Thesis in Solid and Fluid Mechanics, Chalmers University of Technology, Göteborg, Sweden.

36. Bathe, K. J. (2006), Finite element procedures, Klaus-Jurgen Bathe, Prentice Hall, Pearson Education, Inc. ISBN: 978-0-9790049-0-2.

37. Klinkel, S., Wagner, W. (1997), A geometrical non-linear brick element based on the EAS-method, International Journal for Numerical Methods in Engineering, Vol. 40, No. 24, Pg. 4529-4545.

Appendix

The related codes in our routines regarding Newton iteration read as below:

%

%===============ʰʰʰʰʰ====ʰʰʰ= MAIN FEM ANALYSIS PROCEDUR======================================

%

% w Nodal displacements. Let w-i^a be ith displacement

component

% at jth node.

% dw Correction to nodal displacements. Let w_i^a be ith displacement

% Component at jth node. Then dofs (w_1^1, w_2^1, w_3^1, w_1^2,

w_2^2,

% K Global stiffness matrix. Stored as [K_1111 K_1112 K_1121 K_1122… % K_1211 K_1212 K_1221 K_1222… % K_2111 K_2112 K_2121 K_2122…] % % F % % R % b

Force vector. Currently only includes contribution from tractions acting on element faces (i.e. body forces are neglected)

Volume contribution to residual RHS of equation system

%

clear all; clc;

% Run preproc routine

preproc; % dw = zeros (nnode*ndof,1);dw=sparse(dw); w = zeros (nnode*ndof,1);w=sparse(w); K= zeros (nnode*ndof,nnode*ndof);K=sparse(K); R= zeros (nnode*ndof,1);R=sparse(R); b= zeros (nnode*ndof,1);b=sparse(b);

%

% Here we specify how the Newton Raphson iteration should run % Load is applied in nsteps increments;

% tol is the tolerance used in checking Newton-Raphson convergence % maxit is the max no. Newton-Raphson iterations

% relax is the relaxation factor (Set to 1 unless big convergence problems) %

Nsteps = 25;

tol = 0.0001;

maxit = 200;

relax = 1.;

for step=1 : nsteps

loadfactor = step/nsteps; err1 = 1 . ; err2 = 1 . ; nit = 0; disp(‘%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%’) disp(‘%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%’) disp(‘%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%’)

disp (‘%%%%%%%%%%%%%%%%%%% Step Load %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%’)

disp ([step loadfactor])

disp(‘%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%’)

disp(‘%%%%ʩʩʩʩ%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%’)

disp(‘%%%%%%%%%%%%%%ʩ%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%’)

while (err2>tol) & (nit<maxit)

nit = nit + 1;

K = globalstiffness (nrecord, ndof, nnode, cords, …

nelem, maxnodes, nelnodes, connect, materialprops, w);

%Calculate global traction vector

F = globaltraction (nol2, py, ivno, nnode);

%Calculate global residual vector

R = globalresidual (nrecord, ndof, nnode, cords, …

nelem, maxnodes, nelnodes, connect, materialprops, w);

%out of balance vector

for n=1 : nnode*ndof

b(n) = loadfactor*F(n) – R(n) ;

end

for n = 1 : nfix

rw = ndof* (fixnodes (1,n)-1) + fixnodes (2,n) ;

for cl=1 : ndof*nnode ;

K (rw,cl) = 0 ;

end

% Solve for the correction

dw = (K\b) ; % _{Check Convergence} err1 = 0 ; err2 = 0 ; wnorm = 0 ; for n=1 : ndof*nnode w(n) = w(n) + relax*dw (n) ; wnorm = wnorm + w(n)*w(n) ; err1 = err1+dw (n)*dw (n) ; err2 = err2+b(n)*b(n); end

% Store traction and displacement for plotting later

end

% Run post processing routine

School of Engineering, Department of Mechanical Engineering Blekinge Institute of Technology

SE-371 79 Karlskrona, SWEDEN

Telephone: E-mail: