close

Вход

Забыли?

вход по аккаунту

?

467

код для вставкиСкачать
INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING
Int. J. Numer. Meth. Engng. 44, 615—639 (1999)
A SYSTEMATIC CONSTRUCTION OF B-BAR FUNCTIONS
FOR LINEAR AND NON-LINEAR MIXED-ENHANCED
FINITE ELEMENTS FOR PLANE ELASTICITY PROBLEMS
R. PILTNER* AND R. L. TAYLOR
Department of Engineering Mechanics, ¼317.2 Nebraska Hall, ºniversity of Nebraska-¸incoln, NE 68588-0526, º.S.A.
Department of Civil Engineering, SEMM, ºniversity of California at Berkeley, Berkeley, CA 94270, º.S.A.
SUMMARY
In a previous paper a modified Hu—Washizu variational formulation has been used to derive an accurate
four node plane strain/stress finite element denoted QE2. For the mixed element QE2 two enhanced
strain terms are used and the assumed stresses satisfy the equilibrium equations a priori for the linear elastic
case. In this paper an alternative approach is discussed. The new formulation leads to the same accuracy for
linear elastic problems as the QE2 element; however it turns out to be more efficient in numerical
simulations, especially for large deformation problems. Using orthogonal stress and strain functions we
derive B functions which avoid numerical inversion of matrices. The B -strain matrix is sparse and has the
same structure as the strain matrix B obtained from a compatible displacement field. The implementation
of the derived mixed element is basically the same as the one for a compatible displacement element.
The only difference is that we have to compute a B -strain matrix instead of the standard B-matrix.
Accordingly, existing subroutines for a compatible displacement element can be easily changed to
obtain the mixed-enhanced finite element which yields a higher accuracy than the Q4 and QM6 elements.
Copyright 1999 John Wiley & Sons, Ltd.
KEY WORDS: mixed finite element; quadrilateral enhanced strain element
1. INTRODUCTION
Several methods have been developed to improve the performance of the standard four-node
compatible displacement element which yields poor results for problems with bending and, for
plane strain problems, at the nearly incompressible limit. Among the improvements we like to cite
the following methods:
(i) Selected reduced integration and the so-called B-bar method of Hughes.
(ii) The QM6-element of Taylor/Beresford/Wilson, who used four incompatible quadratic
displacement functions to calculate an approximated displacement gradient.
(iii) The method of ‘enhanced strains’ introduced by Simo and Rifai. Simo and Rifai also
demonstrated that the QM6 element can be viewed as an enhanced strain element with
four enhanced strain terms.
* Correspondence to: R. Piltner, Department of Engineering Mechanics, W317.2 Nebraska Hall, University of NebraskaLincoln, Lincoln, NE 68588-0526, U.S.A. E-mail: [email protected]
CCC 0029—5981/99/050615—25$17.50
Copyright 1999 John Wiley & Sons, Ltd.
Received 8 October 1997
Revised 7 May 1998
616
R. PILTNER AND R. L. TAYLOR
(iv) The hybrid stress elements of Pian and Sumihara and Yuan/Huang/Pian.
(v) The mixed finite element QE2 of Piltner and Taylor using a modified Hu—Washizu
variational formulation with bilinear displacement interpolations, seven strain and stress
terms in Cartesian co-ordinates, and two enhanced strain modes.
The enhanced concept became quite popular in recent years and has been used for both linear
and non-linear problems by several researchers, e.g. Simo/Armero, Simo/Armero/Taylor,
Andelfinger/Ramm, Crisfield/Moita/Jelenic/Lyons, Freischläger/Schweizerhof, Glaser/
Armero, Korelc/Wriggers, Nagtegaal/Fox, Roehl/Ramm, de Sausa Neto/Peric/
Huang/Owen, Wriggers/Reese.
In Reference 1, the stress functions for the QE2 element were chosen such that the homogeneous equilibrium equations are satisfied a priori. Therefore the linear version of the QE2 element
can also be considered as a Trefftz-type element (e.g. References 21—28). For the QE2 element
a matrix H has to be computed and inverted. The matrix H has the following block diagonal
structure:
H
0
(1)
0
H
where H is a 3;3 diagonal matrix. So with the assumed stresses of Reference 1 we are left at the
element level to invert the 4;4 symmetric sub-matrix H . In this paper we consider a method to
avoid numerical inversion (except the one needed for the static condensation of the unknowns
associated to the enhanced strain functions).
In addition to the improvement in efficiency of the mixed four-node finite element we address
the subject of non-physical instabilities of non-linear versions of enhanced finite elements
which have been observed by Wriggers and Reese and discussed by Crisfield et al. and De
Sousa Neto et al.
A mixed finite element with four enhanced strain terms has been tested in a series of numerical
examples. The element is denoted B -QE4. In several numerical tests for linear problems it is
shown that the elements B -QE4 and QE2 give identical results. The possibility of coding the
B -QE4 element like a displacement element and its excellent performance make it attractive for
non-linear applications.
H"
2. FINITE ELEMENT FORMULATION FOR THE LINEAR ELASTIC CASE
As in Reference 1 we use the following modified Hu—Washizu variational formulation:
4 2 e2Ee d»!4 u2f d»!1 u2T dS!4 r2(e!Du!eG) d»
%(u, eG, e, r)"
1
(2)
where eG is the enhanced strain field satisfying an orthogonality condition with respect to properly
chosen reference stresses. Without eG we have the classical Hu—Washizu formulation. Details
about the choice of enhanced strain terms are given in Reference 1. The motivation for using
a mixed variational formulation is the following: The good performance of the element QE2
showed that there can be an advantage of assuming stresses in Cartesian co-ordinates. Using both
assumed stresses and assumed strains we will be able to derive a formulation which allows a fast
implementation for both linear and non-linear problems. It will be shown that this implementation
of the mixed element will be very similar to an implementation of a pure displacement element.
Copyright 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 615—639 (1999)
SYSTEMATIC CONSTRUCTION OF B-BAR FUNCTIONS
617
Carrying out the variation in (2) we obtain the following equations:
4 dr2[e!Du!eG] d»"0
(3)
4 de2[Ee!r] d»"0
(4)
!
4 d(Du)2r d»"4 du2f d»#1 du2T dS
4 (deG)2r d»"0
(5)
(6)
The displacement, strain and stress fields are chosen in the following form:
u"Nq, e"Sa, r"Sb, eG"BGk
(7)
where N"N(m, g) is the matrix of compatible shape functions and q contains the nodal displacements of the finite element. The vectors a, b, k are strain, stress and enhanced strain parameters,
respectively. The strain field obtained from the compatible displacement field u can be written as
Du"Bq
(8)
where D is a linear differential operator matrix. The discretization of equations (3)—(6) with
arbitrary variations db, da, dq and dk leads to the following system of equations at the element
level:
0
!H L
LG
b
0
a
0
L2
H
2
0
0
0
q
LG2
0
0
0
k
!H
0
0
"
f
0
(9)
where
4 S2S d»,
H"
4 S2ES d»,
H "
2
4 S2B d»,
L"
4 S2BG d»
LG"
4 N2f d»#1 N2T dS
f "
(10)
From the first two equations of (9) we obtain the strain and stress parameters as
a"H\Lq#H\LGk
(11)
b"H\H a"H\H H\Lq#H\H H\LGk
2
2
2
(12)
Copyright 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 615—639 (1999)
618
R. PILTNER AND R. L. TAYLOR
Substitution of (11) and (12) into (9) gives us the following system of equations:
L2H\H H\L L2H\H H\LG
2
2
LG2H\H H\L LG2H\H H\LG
2
2
or
K C2
C Q
q
f
" k
0
q
f
" k
0
(13)
(14)
where
K"L2H\H H\L, C"LG2H\H H\L, Q"LG2H\H H\LG
2
2
2
(15)
Using a static condensation process with respect to the enhanced strain parameters k we finally
get the element stiffness matrix in the form
k"K!C2Q\C
(16)
In order to avoid calculating the matrices L, H, H , L, LG as well as performing the numerical
2
inversion of H and matrix multiplications in (15), we seek a more efficient way to compute the
matrices K, C and Q. The key for an efficient implementation of the mixed finite element is to
exploit the structures of the matrices H, L, LG for a proper choice of stress and strain functions
collected in the matrix S.
In order to achieve high accuracy with our finite element we choose the strain and stress
functions in Cartesian co-ordinates. This choice is motivated by the good numerical results of the
element QE2 and the Trefftz-type finite elements from References 21—28 for which stresses in
Cartesian co-ordinates have been chosen such that the equilibrium equations are satisfied a priori.
The disadvantage of the set of assumed stresses for element QE2 (see Equation (77) in Reference 1)
is that the matrix of assumed stresses is full and did not lead to a diagonal matrix H. In Reference
1 a block diagonal matrix with one 3;3 diagonal sub-matrix and a 4;4 fully populated
sub-matrix was constructed.
For the present element we do not require that the assumed stresses satisfy equilibrium at the
element level a priori for arbitrary stress parameters b, as was the case in Reference 1. Instead we
require only that the stresses can satisfy the equilibrium equations for proper relations between
the stress parameters b collected in the vector b. This implies that we now need more stress terms
H
than in Reference 1, where seven stress terms have been used.
For the mixed two-dimensional element we can assume the stress field in the form
p
S S S
0 0 0 0 0 0
VV
r" p " 0 0 0 S S S
0 0 0 [b]"Sb
WW
q
0 0 0 0 0 0 S S S
VW
(17)
where the S are linearly independent functions in Cartesian co-ordinates. For the strain field we
H
use the same matrix S as for the stress trial field. Unlike the assumed stresses for element QE2 (see
Copyright 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 615—639 (1999)
SYSTEMATIC CONSTRUCTION OF B-BAR FUNCTIONS
619
equation (77) in Reference 1) the assumed stress components of (17) are not coupled. The resulting
matrix H has the following block diagonal structure:
Hª
0
0
H" 0
Hª
0
0
0
Hª
(18)
If our linearly independent trial functions S ( j"1, 2, 3) do not satisfy the orthogonality
H
condition
4 SGSH d» O0
"0 for iOj
for i"j
(19)
we will get a fully populated 3;3 matrix Hª . Since the matrix H has to be inverted it is better to
orthogonalize a set of linearly independent functions so that Hª and therefore H become diagonal.
Possible choices for our initial (non-orthogonal) set of linearly independent trial functions S
H
are
S "1, S "xN ,
S "yN
(20)
and
S "1, S "m, S "g
(21)
where xN , yN are Cartesian co-ordinates with origin at the centre of the quadrilateral finite element.
The co-ordinates m, g are obtained as a linear combination of the Cartesian co-ordinates xN , yN from
the relationship
m
b
!a
xN
1
"
(22)
g
J !b
a
yN
The relationship between the global Cartesian co-ordinates x, y and the co-ordinates xN , yN , shifted
to the center of the element, is given as
xN "x!a ,
yN "y!b
(23)
where
a " (x #x #x #x ), b " (y #y #y #y )
The coefficients of the linear mapping (22) are given as
(24)
a " (!x #x #x !x ), b " (!y #y #y !y )
a " (!x !x #x #x ), b " (!y !y #y #y )
(25)
J "a b !a b
Finally, we can express the element coordinates xN , yN and m, g through the natural element
co-ordinates m, g via
xN "a m#a g#a mg
yN "b m#b g#b mg
Copyright 1999 John Wiley & Sons, Ltd.
(26)
Int. J. Numer. Meth. Engng. 44, 615—639 (1999)
620
R. PILTNER AND R. L. TAYLOR
and
J
m"m# mg
J
J
g "g# mg
J
(27)
where
a " (x !x #x !x ),
b " (y !y #y !y ),
(28)
and
J "a b !b a , J "a b !a b
(29)
Since the expressions (27) are simpler than (26) we choose [SK "1, SK "m, SK "g ] as initial set of
linearly independent trial functions to construct an orthogonal set of linearly independent
functions [S , S , S ].
This can be achieved by a Gram—Schmidt-orthogonalization process through the formula
SK S d»
H I
H\ 4
S "SK ! S
H
H
I
I S S d»
I I
4
(30)
The resulting set of orthogonal functions is
S "1
J
J
J
1
J
S "m! "m# mg! "
b xN !a yN ! (31)
3J
J
3J
J
3
J
J
J
1
J
S "g ! !f S "g# mg! !f S "
!b xN #a yN ! !f S
J
J
3J
3
3J
where
2J J
f"
(32)
3J!J#3J
It should be noted that the functions S and S used for the stress/strain approximations can be
expressed as polynomials in xN , yN , m, g or m, g-co-ordinates, whereas in a pure displacement
formulation the stresses and strains are rational functions in natural element co-ordinates for
a general quadrilateral element. Because the functions S can be expressed as a complete set of
H
linear functions in Cartesian co-ordinates, the special case of the seven parameter stress field of
the QE2 element satisfying the equilibrium equations a priori is included in the stress field (17).
With the choice of functions [S , S , S ] from equation (31) the entries of the matrix S for the
stress and strain fields are defined.
Copyright 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 615—639 (1999)
621
SYSTEMATIC CONSTRUCTION OF B-BAR FUNCTIONS
The compatible displacement field u is interpolated with standard shape functions according
u
1
u
u
" (1#m m)(1#g g) G " N G "Nq
(33)
G
G
G v
v
v
4
G
G
G
G
An admissible enhanced strain field with four terms can be obtained from the incompatible
displacement field
u"
uG"
N
0
N
0
0
N
0
N
j
j
j
j
(34)
where
N"1!m, N"1!g
(35)
Because the constant strain patch test cannot be satisfied when using the exact gradient of the
incompatible shape functions
*
y (m, g)
!y (m, g)
*x
1
E
K
N"
N"
H
H
*
J(m, g) !x (m, g)
x (m, g)
E
K
*y
*
*m
N"J\(m, g)
H
*
*g
*
*m
N
H
*
*g
(36)
the approximated gradient
*
*x
1
y (0, 0)
!y (0, 0)
E
K
N"
N"
H
H
J(m, g) !x (0, 0)
*
x (0, 0)
E
K
*y
*
*m
J
N" J\
H
J(m, g) *
*g
*
*m
N
H
*
*g
(37)
as proposed by Taylor/Beresford/Wilson in Reference 4 is used to get the following admissible
enhanced strain field:
eG
!b m
0
b g
0
VV
2
eG" eG "BGk"
0
a m
0
!a g
WW
J(m, g)
cG
a m
!b m !a g
b g
VW
j
j
j
j
(38)
where the determinant J of the Jacobian matrix can be expressed as
J(m, g)"J #J m#J g
(39)
With the assumed displacement, strain and stress fields, all integrals in the element formulation
can be calculated exactly with 2;2 Gauss-integration since the expression for the matrices H, H ,
2
L, and LG contain only polynomials in m and g. This is in contrast to the situation in the
Copyright 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 615—639 (1999)
622
R. PILTNER AND R. L. TAYLOR
displacement element formulation where we cannot get exact integral results for an arbitrary
quadrilateral element shape. This is due to the factor 1/J(m, g) in all integrals for the stiffness
matrix.
The diagonal matrix H is defined through the entries of the 3;3 matrix Hª . Through
analytical integration we get the diagonal coefficients of the matrix Hª as
hK "4J
Kh " J (3j!j#3)
4(3#2j!j#2j#2j j!j)
J
hK "
3(3j!j#3)
where j "J /J and j "J /J . The strain matrices B and BG have the structure
N
0
N
0
N
0
N
0
V
V
V
V
N
0
N
0
N
0
N
B" 0
W
W
W
W
N
N
N
N
N
N
N
N
W
V
W
V
W
V
W
V
and
*
*
N
0
N
0
*x *x BG"
0
*
N
*y 0
*
N
*y *
N
*y *
N
*x *
N
*y *
N
*x (40)
(41)
(42)
where */*x and */*y indicate that the derivatives of the incompatible shape functions are an
approximation given by equation (37).
It is important to notice that due to the sparse structure of S, B, and BG and the fact that these
matrices have repeated entries for their coefficients the matrices L and LG have only a few
non-vanishing different coefficients. The non-vanishing coefficients are denoted by f GH , f GH , f GI ,
V W V
f GI , where (i"1, 2, 3), ( j"1, 2, 3, 4), (k"1, 2) and given by
W
*
f GH" S (m, g) N (m, g) d»
V
G
*x H
4
4
*
S (m, g) N (m, g) d»
G
*y H
4
*
S (m, g) N (m, g) d»
G
*x I
f GH"
W
f GI "
V
(43)
*
4 SG (m, g) *y NI (m, g) d»
f GI "
W
Copyright 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 615—639 (1999)
SYSTEMATIC CONSTRUCTION OF B-BAR FUNCTIONS
623
There is no particular advantage to give analytical expressions for the results of the integrals in
equations (43). The exact values for f GH , f GH , f GI , and f GI can be obtained with 2;2 Gaussian
V W V
W integration. Having calculated the exact coefficients of the diagonal matrix H and the sparse
matrices L and LG we can calculate the strain parameters a from relationship (11). Now we can
write the assumed strain field in the form:
e"B q#B Gk
(44)
NM 0 NM 0 NM 0 NM 0
V
V
V
V
B " 0 NM 0 NM 0 NM 0 NM W
W
W
W
NM NM NM NM NM NM NM NM V
W
V
W
V
W
V
W
(45)
where
and
NM V
B G"
0
0
NM W
NM V
NM W
NM V
0
NM W
0
NM W
NM V
(46)
The matrices B , B G have the same sparse structure as the matrices B, BG. Taking into account that
S "1, the NM functions can be written as
f H f H
f H
NM H (m, g)" V # V S (m, g)# V S (m, g)
V
hK
hK
hK
f H f H
f H
NM H (m, g)" W # W S (m, g)# W S (m, g)
W
hK
hK
hK
f I
f I
f I
NM I
(m, g)" V# V S (m, g)# V S (m, g)
V
hK
hK
hK
f I
f I
f I
NM I (m, g)" W# W S (m, g)# W S (m, g)
W
hK
hK
hK
(47)
where ( j"1, 2, 3, 4) and (k"1, 2).
Note that the matrices B and B G contain polynomials in m, g whereas the matrices B and BG
contain rational functions of the natural element co-ordinates m and g. Constant stress and strain
modes can be represented exactly with the proposed formulation involving the matrices B and B G.
Instead of computing the matrices K, C and Q from equations (15), which involves time
consuming matrix multiplications of H, L, LG, H we can compute them faster by using the sparse
2
strain matrices B and B G:
4 B 2EB d»,
K"
Copyright 1999 John Wiley & Sons, Ltd.
4 B G2EB d»,
C"
4 B G2EB G d»
Q"
(48)
Int. J. Numer. Meth. Engng. 44, 615—639 (1999)
624
R. PILTNER AND R. L. TAYLOR
The external load vector involves compatible shape functions from equation (33) and is given by
4 N2f d»#1 N2T dS
f "
(49)
3. FORMULATION FOR LARGE DEFORMATION OF HYPERELASTIC MATERIALS
In this section we will first use tensor notation and later switch to matrix and vector notation. For
an overview on non-linear elasticity we refer to Beatty’s review articles. For the large
deformation case we will use the first (un-symmetric) Piola—Kirchhoff tensor P which is related to
the second (symmetric) Piola—Kirchhoff stress tensor S , the Kirchhoff stresses s and the Cauchy
stresses r through the relation
P"FS "sF\2"JrF\2
(50)
where F is the deformation gradient
*x
F" "I#Grad u
*X
(51)
In this notation, x is the position vector in the current configuration whereas X denotes the
position vector in the reference configuration. Using u as the displacement vector, the position
vectors are related through the equation
x"X#u
(52)
The displacement gradient with respect to the reference configuration is given for our twodimensional example by
*u *u
*X *½
Grad u"
(53)
*v *v
*X *½
The relations for the symmetric Kirchhoff stress tensor s and the Cauchy stress tensor r are given
by
s"PF2"FS F2
1
r" s
J
(54)
Here we want to consider the case of a hyperelastic material which is characterized through
a frame invariant stored energy function ¼(X, F). For objectivity the strain energy function must
satisfy the relation
W(X, QF)"¼(X, F)
(55)
for all proper orthogonal Q. For a hyperelastic material the un-symmetric Piola—Kirchhoff stress
tensor is given by
*¼
P"
*F
Copyright 1999 John Wiley & Sons, Ltd.
(56)
Int. J. Numer. Meth. Engng. 44, 615—639 (1999)
SYSTEMATIC CONSTRUCTION OF B-BAR FUNCTIONS
625
In the finite element implementation we will make use of the incremental relation
*P"A : *F
(57)
*P "A *F
G(
G(I* I*
(58)
or in index notation
The fourth order tensor A, which is called the first elasticity tensor, is given through
*¼
*¼
or A "
A"
G(I*
*F *F
*F*F
G( I*
(59)
Using an assumed displacement gradient G and an enhanced displacement gradient Grad u
in the reference configuration the modified Hu—Washizu formulation for large deformations can
be given as
4 Grad du : P d»#d% .
d%"
(60)
4 dP : [G!Grad u!Grad u ] d»
!
4 dG : *F >
*¼
#
F
I
!P d»
G
4 d[Grad u ] : P d»
#
where dF"Grad du and
*u
*u
*X
*½
uJ
uJ
G " 6 7 , Grad u "
vJ
vJ
*v
*v
6 7
*X
*½
P"
P
P
P
P
(61)
(62)
have been used. For a discretization of (60) we switch to matrix and vector notation. Using the
vectors
P2"[P , P , P , P ]
u2"[u , v , u , v ]
6 7 7 6
g2"[uJ , vJ , uJ , vJ ]
6 7 7 6
u2 "[u , v , u , v]
6
7
7
6
Copyright 1999 John Wiley & Sons, Ltd.
(63)
Int. J. Numer. Meth. Engng. 44, 615—639 (1999)
626
R. PILTNER AND R. L. TAYLOR
The variational formulation can be expressed as
4 du2P d»#d% .
d%"
(64)
4 dP2 [g2!
u!
u ] d»
!
4 dg [P (g)!P] d»
#
4 d
u2P d»
#
where P denotes the vector of Piola—Kirchhoff stresses obtained from *¼/*F by using the
assumed displacement gradient G . An increment for the Piola—Kirchhoff stress depending on an
increment *g can be calculated from
*P "A *g
2 (65)
An expression for the tangent moduli A for one of the chosen constitutive models is given in
2
Section 4.
For the increments of stresses and displacement gradients we assume the functions
*P"S*b, *g "S*a, (*u)"B*q, (*u )"B*k
(66)
S S S
0 0 0 0 0 0 0 0 0
0 0 0 S S S
0 0 0 0 0 0
S"
0 0 0 0 0 0 S S S
0 0 0
0 0 0 0 0 0 0 0 0 S S S
(67)
where
N
B"
6
0
N
7
0
0
N
7
0
N
6
N
6
0
N
7
0
0
N
7
0
N
6
N
6
0
N
7
0
0
N
7
0
N
6
N
6
0
N
7
0
0
N
7
0
(68)
N
6
B is the matrix for the enhanced displacement gradient. Unfortunately, for the large deformation case we cannot use a sparse matrix for the enhanced displacement gradient. Using a sparse
matrix B yields an element which is not frame invariant. The following B has been chosen for
Copyright 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 615—639 (1999)
SYSTEMATIC CONSTRUCTION OF B-BAR FUNCTIONS
627
the proposed enhanced element:
bm
!b b m !b b g
bg
am
!a a m !a a g
ag
1
B"
J(m, g) !a b m
a b m
a b g
!a b g
!a b m
a b m
a b g
!a b g
(69)
The above matrix can be obtained by using the following form for the enhanced displacement
gradient H:
J
F"H" J\H J\2
J (70)
where the Jacobian at the element origin is given by
a b
J " a b
and H is chosen as
(71)
j m j m
(72)
j g j g
Recently instabilities have been detected in a non-linear version of the QM6 element. The QM6
element proposed by Taylor/Beresford/Wilson in 1976 can be viewed as an enhanced strain
element with four enhanced strain terms. A non-linear version was proposed by Simo and
Armero. The non-physical instabilities in enhanced strain finite element simulations can occur in
a hyperelastic material under uniform compression, as initially observed by Wriggers and
Reese. The problem has been investigated recently by several researchers, e.g. Crisfield et al.,
de Souza Neto et al., Glaser and Armero.
In order to overcome the element instability problems, Glaser and Armero proposed to
exchange the original basis for the enhanced deformation gradient
H "
H "
j m j g
j m j g
(73)
by either expression (72) or by
j m
j m#j g
(74)
j m#j g
j g
Glaser and Armero considered the following two transformations for constructing the enhanced
part of the deformation gradient:
H "
Copyright 1999 John Wiley & Sons, Ltd.
J
(a) H" J H J\
J (75)
J
(b) H" J\2 H J\
J (76)
Int. J. Numer. Meth. Engng. 44, 615—639 (1999)
628
R. PILTNER AND R. L. TAYLOR
The enhanced deformation gradient in Reference 13 is chosen in the form
F"F H
(77)
as proposed by Simo/Armero and Taylor. Korelc and Wriggers proposed an enhanced
deformation gradient of the form
J
F"H" J\2 H J\
J (78)
where H is the same as in equation (72). However, Glaser and Armero pointed out that the use of
(78) yields a non-objective element. For the proposed element H is chosen as in equation (70).
Using (70) an objective element is obtained (see example 5.7).
After elimination of the strain parameters the assumed displacement gradient for an increment
takes the form
*g "B *q#B *k
(79)
where
B "
NM 6
0
NM 7
0
NM 6
0
0
NM 7
0
NM 7
0
NM 6
0
NM 7
0
NM 6
NM 6
0
NM 7
0
0
NM 7
0
NM 6
NM 6
0
NM 7
0
0
NM 7
0
(80)
NM 6
and
bM bM B " bM bM bM bM bM bM bM bM bM bM bM bM bM bM (81)
The functions bM (m, g) are defined as
GH
¸
¸
¸
bM (m, g)" GH # GH S (m, g)# GH S (m, g)
GH
hK
hK
hK
(82)
and
4 SI (m, g)bGH (m, g) d»
¸I "
GH
(83)
The functions b are the coefficients of the matrix B defined with equation (69).
GH
Because the assumed displacement gradient can be expressed in the form (79) involving the
sparse matrix B , the element tangent stiffness matrix can be computed in a very economic
manner. The element tangent stiffness matrix for the large deformation case is given as
k "K !C2 Q\C
2
2
2 2 2
Copyright 1999 John Wiley & Sons, Ltd.
(84)
Int. J. Numer. Meth. Engng. 44, 615—639 (1999)
SYSTEMATIC CONSTRUCTION OF B-BAR FUNCTIONS
629
where the matrices K , C , and Q can be computed in the following form:
2 2
2
4 B 2A2 B d»
K "
2
4 B A2 B d»
C "
2
(85)
4 B 2A2B d»
Q "
2
The residual vectors for step n are obtained by computing
4 B 2P d»,
f "
4 B 2PL d»
f "
(86)
4. CHOSEN CONSTITUTIVE MODELS
For the examples in the present paper we consider two constitutive models.
4.1. Model 1
The stored energy function for the first constitutive model is
K (I , J)" k(I !3)!k ln J# j(ln J)
¼"¼(I , I , I )"¼
where I , I , I are the invariants of the right Cauchy—Green deformation tensor
C"F2F
(87)
(88)
and
J"(I "(det C"det F
The second Piola—Kirchhoff stress tensor is obtained from
S "2
*¼ *I
*¼ *I
*¼ *I
*¼
#2
#2
"2
*C
*I *C
*I *C
*I *C
(89)
(90)
where
*I
*I
"I,
"(det C)C\"JC\"JF\F\2
*C
*C
(91)
For the above strain energy function ¼ we have
*¼ 1
*¼
*¼ *¼ *J *¼ 1
" k,
"0,
"
"
*I
2
*I
*I
*J *I
*J 2J
Copyright 1999 John Wiley & Sons, Ltd.
(92)
Int. J. Numer. Meth. Engng. 44, 615—639 (1999)
630
R. PILTNER AND R. L. TAYLOR
and
1
*¼ k
" #j ln J
J
J
*J
(93)
so that we get the constitutive relationship for S in the form
S "kI#[!k#jln J]F\F\2
(94)
For the other stresses we obtain
P"k[F!F\2]#jln JF\2, s"k[b!I]#jln JI,
k
j
r" [b!I]# ln JI
J
J
(95)
where
b"FF2
(96)
The first elasticity tensor A can be obtained conveniently by using components of the unsymmetric Piola—Kirchhoff stress tensor P, given from (56) by
1
1
P "kF # [!k#jln J]F , P "kF # [!k#jln J]F
J
J
1
1
P "kF ! [!k#jln J]F , P "kF ! [!k#jln J]F
J
J
(97)
Using relationship (59) we get
A
1
"k# [k#j!jln J]F F ,
J
1
A
"k# [k#j!jln J]F F
J
1
1
A
" [!k#jln J]# [k#j!jln J]F F
J
J
1
1
"! [k#j!jln J]F F
A
"! [k#j!jln J]F F , A
J
J
(98)
1
1
A
"! [k#j!jln J]F F , A
"! [k#j!jln J]F F
J
J
1
A
"k# [k#j!jln J]F F ,
J
1
A
"k# [k#j!jln J]F F
J
1
1
A
"! [!k#jln J]# [k#j!jln J]F F
J
J
In matrix notation the incremental constitutive equation for a point under consideration can be
written as
*P"A *F
2
Copyright 1999 John Wiley & Sons, Ltd.
(99)
Int. J. Numer. Meth. Engng. 44, 615—639 (1999)
SYSTEMATIC CONSTRUCTION OF B-BAR FUNCTIONS
631
where
A
A
A
A
A
A
A
A
A " 2
A
A
A
A
A
A
A
A
*P2"[P , P , P , P ]
(100)
*F2"
(*u2)"[*u , *v , *u , *v ]
6
7
7
6
It should be mentioned that although the symmetric matrix A is singular in the undeformed
2
configuration, the element stiffness matrix has the correct number of zero eigenvalues (i.e. three
zero eigenvalues associated to the rigid body modes).
4.2. Model 2
The second constitutive model used in the examples is taken from a paper by Knowles and
Sternberg. The strain energy function for this model is
¼(I , J)" k[I J\#2J!4]
(101)
Reference 31 gives a detailed discussion of this material model which is a special case of
a Blatz—Ko material.
5. NUMERICAL EXAMPLES
Several linear and non-linear problems have been selected to test the performance of the proposed
four-node element denoted as B -QE4. For all problems 2;2 Gaussian integration is used. The
element passes the patch test.
5.1. Beam bending
A beam modeled with five elements is subjected to two load cases (Figure 1). Plane stress
conditions are assumed in the model. The results of different elements for the maximum
displacement at point A and the normal stress p at point B are given in Table I. The elements
VV
used in the examples are: (i) the bilinear isoparametric displacement element Q4 [32, 33, 34].
(ii) the enhanced strain element QM6 of Taylor/Beresford/Wilson, (iii) the hybrid stress element
P—S of Pian/Sumihara, (iv) the enhanced mixed element QE2 by Piltner/Taylor, (v) the
proposed element B -QE4. The elements B -QE4 and QE2 give identical results.
5.2. Mesh distortion test for beam bending
In this test a beam under bending is analysed with only two plane stress elements (Figure 2).
The degree of distortion of the element is measured with the distortion parameter *. The material
parameters are E"1500 and l"0)25. From Table II we can see that the elements QE2 and
B -QE4 show the least sensitivity to mesh distortion even for every severe distortions. Elements
QE2 and B -QE4 provide identical results.
Copyright 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 615—639 (1999)
632
R. PILTNER AND R. L. TAYLOR
Figure 1. Finite element mesh for cantilever beam problem
Table I. Comparison of plane stress solutions obtained
with four node elements for cantilever beam problems
Case 1
Element
v
Q4
QM6
P-S
QE2
B -QE4
Exact
45)65
96)07
96)18
96)5
96)5
100
Case 2
p
VV
!1604
!2497
!3001
!3004
!3004
!3000
v
50)96
97)98
98)05
98)26
98)26
102)6
p
VV
!2146
!3235
!3899
!3906
!3906
!4050
Figure 2. Cantilever beam for the mesh distortion test
5.3. Cook’s membrane problem
The material parameters for the plane stress structure shown in Figure 3 are E"1 and l".
Displacement and stress results are listed in Table III.
Copyright 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 615—639 (1999)
SYSTEMATIC CONSTRUCTION OF B-BAR FUNCTIONS
633
Table II. Displacenment v of cantilever beam (Figure 2) for
different values of the mesh distortion parameter *
*
Q4
QM6
0
0)5
1
2
3
4
4.9
28)0
21)0
14)1
9)7
8)3
7)2
6)2
100)0
80)9
62)7
54)4
53)6
51)2
46)8
P-S
QE2
Displacement v
100)0
100)0
81)0
81)2
62)9
63)4
55)0
56)5
54)7
57)5
53)1
57)9
49)8
56)9
B -QE4
Exact
100)0
81)2
63)4
56)5
57)5
57)9
56)9
100
100
100
100
100
100
100
Figure 3. Cook’s membrane problem: plane stress structure with unit load uniformly distributed along right edge
(E"1, l")
5.4. Short beam under bending
A short beam, supported at the left end and subjected to a couple at the right edge, is modeled
with a mesh as shown in Figure 4. The purpose of this example is to test the influence of the mesh
pattern on the stress distribution in the finite element solution. In the QM6 element the stress has
the tendency to line up with the element edges of the element in the center of the mesh (Figure
5(a)). On the other hand, elements QE2 and B -QE4 preserve symmetries in the solution domain
although a very coarse irregular mesh is used (Figure 5(b)).
5.5. Single element test for a Neo-Hookean material model
Wriggers and Reese noticed that the enhanced element Q1/E4 can experience problems in
the analysis of Neo-Hookean hyperelastic materials under uniform compression. They suggested
the single element test shown in Figure 6. The material model (87) is used with the following
material parameters: j"40 000, k"80)2. At 30)4 per cent compression the element tangent
Copyright 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 615—639 (1999)
634
R. PILTNER AND R. L. TAYLOR
Table III. Results for the problem shown in
Figure 3
Element
N"2
N"4
N"16
Q4
QM6
P-S
QE2
B -QE4
Displacement
11)85
21)05
21)13
21)35
21)35
v at C
18)30
23)02
23)02
23)04
23)04
23)43
23)88
23)88
23)88
23)88
Q4
QM6
P-S
QE2
B -QE4
Maximum stress at A
0)1078
0)1814
0)1773
0)2225
0)1854
0)2241
0)1956
0)2261
0)1956
0)2261
0)2353
0)2364
0)2364
0)2364
0)2364
Figure 4. Finite element mesh for a short beam under bending
matrix of the Q1/E4 element has two negative eigenvalues. However, the element tangent matrix
should have only one negative eigenvalue for this deformation. At 39 per cent compression the
element Q1/E4 yields a non-physical instability. The system tangent matrix has one negative
eigenvalue at this deformation level. On the other hand the element B -QE4 keeps only one
negative eigenvalue for the element tangent matrix, up to almost 100 per cent compression, and
the system matrix has no negative eigenvalues.
5.6. Multi-element compression test for a Neo-Hookean material model
One-half of a block with unit sides is discretized with a 10;20 element mesh. The material
parameters for the Neo-Hookean model are chosen as in example 5.5. The element Q1/E4 yields
a non-physical instability at 30)6 per cent. However, the first physical instability should appear at
Copyright 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 615—639 (1999)
SYSTEMATIC CONSTRUCTION OF B-BAR FUNCTIONS
635
Figure 5. (a) Stress p obtained with element QM6; (b) Stress p obtained with element QE2 and B -QE4
VV
VV
Copyright 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 615—639 (1999)
636
R. PILTNER AND R. L. TAYLOR
Figure 6. Finite element under compression
Figure 7. (a) undeformed beam; (b) deformed beam and FE mesh
50)12 per cent compression. The element B -QE4 is able to detect this physical instability and the
following instabilities: At 50)12 per cent compression the system matrix based on the B -QE4
elements has one negative eigenvalue. At 50)26 per cent compression the system matrix has two
negative eigenvalues, and at 58)2 per cent compression the discretization with B -QE4 elements
shows three negative eigenvalues.
5.7. Objectivity test for large deformation
The beam shown in Figure 7(a) is fixed at the right end and subjected to a prescribed
displacement as illustrated in Figure 7(b). After calculating the displacements and stresses the
Copyright 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 615—639 (1999)
SYSTEMATIC CONSTRUCTION OF B-BAR FUNCTIONS
637
Figure 8. Finite element discretization for example 5)8
Figure 9. Material 2. Cauchy stress vs. axial stretch
structure is rotated and displacements and stresses are calculated for each rotation angle
(h"10°, 20°, . . . , 180°). For the B -QE4 element the results are independent of the rotation angle.
It should be noticed that the results are not frame independent if we use enhanced strains based
on equation (76) instead of (70).
5.8. Material 2 subjected to a homogenous plane deformation
One-half of a block of dimensions 1;1 is discretized with non-rectangular elements (Figure 8).
The block is subjected to plane strain uni-axial stress parallel to the x -axis. Material model 2 is
used with k"100. The analytical solution for p in this problem is given as a function of the
stretch j (Figure 9):
p "100(1!j\)
(102)
For every element in the mesh stresses can be calculated exactly.
Copyright 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 615—639 (1999)
638
R. PILTNER AND R. L. TAYLOR
In Figure 9, the solid line represents the analytical solution and the circles stand for stress
values obtained numerically for a chosen finite element in the mesh.
6. CONCLUDING REMARKS
For linear problems, the elements B -QE4 and QE2 give identical results and show an overall very
good behaviour. The advantage of element B -QE4 is that it can be coded very similar to
a displacement element and therefore runs faster than element QE2. The key for a fast linear and
non-linear element implementation is to use orthogonal stress and strain functions in the mixed
variational formulation. Due to the use of orthogonal functions matrix inversions for strain and
stress parameters at the element level can be avoided in a finite element program. In the current
implementation standard shape function derivatives N , N are exchanged by NM and NM
V
W
V
W
functions derived in the paper. For non-linear problems the displacement like implementation of
the mixed element B -QE4 has the obvious advantage that it is time saving compared to an
implementation of the mixed variational formulation with matrix inversions for strain and stress
parameters.
REFERENCES
1. R. Piltner and R. L. Taylor, ‘A quadrilateral finite element with two enhanced strain modes’, Int. J. Numer. Meth.
Engng., 38, 1783—1808 (1995).
2. T. J. R. Hughes, ‘Equivalence of finite elements for nearly incompressible elasticity’, J. Appl. Mech., 44, 181—193 (1977).
3. T. J. R. Hughes, ‘Generalizing of selective integration procedures to anisotropic and nonlinear media’, Int. J. Numer.
Meth. Engng., 15, 1413—1418 (1980).
4. R. L. Taylor, P. J. Beresford and E. L. Wilson, ‘A non-conforming element for stress analysis’, Int. J. Numer. Meth.
Engng., 10, 1211—1219 (1976).
5. J. C. Simo and M. S. Rifai, ‘A class of mixed assumed strain methods and the method of incompatible modes’, Int. J.
Numer. Meth. Engng., 29, 1595—1638 (1990).
6. T. H. H. Pian and K. Sumihara, ‘Rational approach for assumed stress elements’, Int. J. Numer. Meth. Engng., 20,
1685—1695 (1984).
7. K.-Y. Yuan, Y.-S. Huang and T. H. H. Pian, ‘New strategy for assumed stresses for 4-node hybrid stress membrane
element’, Int. J. Numer. Meth. Engng., 36, 1747—1763 (1993).
8. J. C. Simo and F. Armero, ‘Geometrically non-linear enhanced strain mixed methods and the method of incompatible
modes’, Int. J. Numer. Meth. Engng., 33, 1413—1449 (1992).
9. J. C. Simo, F. Armero and R. L. Taylor, ‘Improved versions of assumed enhanced strain tri-linear elements for 3D
finite deformation problems’, Comput. Methods Appl. Mech. Engng., 110, 359—386 (1993).
10. U. Andelfinger and E. Ramm, ‘EAS-elements for two-dimensional, three-dimensional, plate and shell structures and
their equivalence to HR-elements’, Int. J. Numer. Meth. Engng., 36, 1311—1337 (1993).
11. M. A. Crisfield, G. F. Moita, G. Jelenic and L. P. R. Lyons, ‘Enhanced lower-order element formulations for large
strains’, in D. R. J. Owen and E. Onate (eds.), Computational Plasticity—Fundamentals and Applications, Pineridge
Press, Swansea, 1995, pp. 293—320.
12. C. Freischläger and K. Schweizerhof, ‘On a systematic development of trilinear three-dimensional solid elements
based on Simo’s enhanced strain formulation’, Int. J. Solids Struct., 33, 2993—3017 (1996).
13. S. Glaser and F. Armero, ‘Recent developments in the formulation of assumed enhanced strain finite elements for
finite deformation problems’, Report No. ºCB/SEMM-95/13, University of California, Berkeley, 1995.
14. S. Glaser and F. Armero, ‘On the formulation of enhanced strain finite elements in finite deformations’, Engng.
Comput., 1996, submitted.
15. J. Korelc and P. Wriggers, ‘Consistent gradient formulation for a stable enhanced strain method for large deformations’, Engng. Comput., 13, 103—123 (1996).
16. P. Wriggers and J. Korelc, ‘On enhanced strain methods for small and finite deformations of solids’, Comput. Mech.,
18, 413—428 (1996).
17. J. C. Nagtegaal and D. D. Fox, ‘Using assumed enhanced strain elements for large compressive deformation’, Int. J.
Solids Struct., 33, 3151—3159 (1996).
18. D. Roehl and E. Ramm, ‘Large elasto-plastic finite element analysis of solids and shells with the enhanced assumed
strain concept’, Int. J. Solids Struct., 33, 3215—3237 (1996).
Copyright 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 615—639 (1999)
SYSTEMATIC CONSTRUCTION OF B-BAR FUNCTIONS
639
19. E. A. de Souza Neto, D. Peric, G. C. Huang and D. R. J. Owen, ‘Remarks on the stability of enhanced strain elements
in finite elasticity and elastoplasticity’, in D. R. J. Owen and E. Onate (eds.), Computational Plasticity—Fundamentals
and Applications, Pineridge Press, Swansea, 1995,
20. P. Wriggers and S. Reese, ‘A note on enhanced strain methods for large deformations’, Comput. Methods Appl. Mech.
Engng., 135, 201—209 (1996).
21. J. Jirousek and L. Guex, ‘The hybrid-Trefftz finite element model and its application to plate bending’, Int. J. Numer.
Meth. Engng., 23, 651—693 (1986).
22. J. Jirousek, ‘Hybrid-Trefftz plate bending elements with p-method capabilities’, Int. J. Numer. Meth. Engng., 24,
1367—1393 (1987).
23. J. Jirousek, A. Venkatesh, A. P. Zielinski and H. Rabemanantosa, ‘Comparative study of p-extensions based on
conventional assumed displacement and hybrid-Trefftz FE models’, Comput. Struct., 46, 261—278 (1993).
24. J. Jirousek, ‘Improvement of computational efficiency of the 9 DOF triangular hybrid-Trefftz plate bending element’,
Int. J. Numer. Meth. Engng., 23, 2167—2168 (1986).
25. A. Venkatesh and J. Jirousek, ‘An improved 9 DOF Hybrid-Trefftz triangular element of plate bending’, Engng.
Comput., 4, 207—222 (1987).
26. R. Piltner, ‘Recent developments in the Trefftz method for finite element and boundary element applications’, Adv.
Engng. Software, 24, 107—115 (1995).
27. R. Piltner, ‘A quadrilateral hybrid-Trefftz plate bending element for the inclusion of warping based on a threedimensional plate formulation’, Int. J. Numer. Meth. Engng., 33, 387—408 (1992).
28. A. P. Zielinski and O. C. Zienkiewicz, ‘Generalized finite element analysis with T-complete solution functions’, Int. J.
Numer. Meth. Engng., 21, 509—528 (1985).
29. M. F. Beatty, ‘Introduction to Nonlinear Elasticity’, in M. M. Carroll and M. A. Hayes (eds.), Nonlinear Effects in
Fluids and Solids, Plenum Press, New York/London, 1996, pp. 13—112.
30. M. F. Beatty, ‘Topics in finite elasticity: hyperelasticity of rubber, elastomers, and biological tissues—with examples’,
Appl. Mech. Rev., 40, 1699—1734 (1987).
31. J. K. Knowles and Eli Sternberg, ‘On the failure of ellipticity and the emergence of discontinuous deformation
gradients in plane finite elastostatics’, J. Elasticity, 8(4), 329—379 (1978).
32. I. C. Taig, ‘Structural analysis by the matrix displacement method’, English Electric Aviation Report No. S017, 1961.
33. T. J. R. Hughes, ¹he Finite Element Method, Prentice-Hall, Englewood Cliffs, N.J., 1987.
34. K.-J. Bathe and E. L. Wilson, Numerical Methods in Finite Element Analysis, Prentice-Hall, Englewood Cliffs, N.J.,
1976.
Copyright 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 615—639 (1999)
Документ
Категория
Без категории
Просмотров
2
Размер файла
304 Кб
Теги
467
1/--страниц
Пожаловаться на содержимое документа