close

Вход

Забыли?

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

?

Transient algorithms for heat transfer general developments and approaches for theoretically generating Nth-order time-accurate operators including practically useful second-order forms

код для вставкиСкачать
INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING
Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999)
STRESS-BASED FINITE ELEMENT ANALYSIS OF
PLANE PLASTICITY PROBLEMS
ZDZISLAW WI E
V CKOWSKI?;? , SUNG-KIE YOUN? AND BYUNG-SIK MOONД
Department of Mechanical Engineering, Korea Advanced Institute of Science and Technology,
373-1, Kusong-dong, Yusong-gu, Taejon 305-701, Korea
SUMMARY
A stress-based model of the nite element method is evolved for two-dimensional quasi-static plasticity
problems. The self-equilibrating elds of stresses are constructed by means of the Airy stress function,
which is approximated by three types of elements: the Bogner?Fox?Schmit rectangle, the Hsieh?Clough?
Tocher triangle and its reduced variant. Traction boundary conditions are imposed by the use of the Lagrange
multiplier method which gives the possibility of calculation of displacements for boundary points. The concept
of multi-point-constraints elements is applied in order to facilitate the application of this technique. The
iterative algorithm, analogous to the closest-point-projection method commonly used in the displacement-based
nite element model, is proposed for solving non-linear equations for each load increment. Two numerical
examples with stress- and displacement-controlled load are considered. The results are compared with those
obtained by the displacement model of FEM. Bounds for limit loads are obtained. Copyright ? 1999 John
Wiley & Sons, Ltd.
KEY WORDS:
stress function
nite element method; equilibrium model of FEM; plasticity; load limit bounds; variational inequalities;
INTRODUCTION
The stress-based nite element method utilizing the principle of complementary work or stationary
complementary energy is used rather rarely in the engineering practice. The possibility of application of this model of FEM is mentioned only in several nite element textbooks,1?4 written among
many others. It seems that the main reasons of the rare use of the model in comparison with the
displacement nite element method are two facts: the statically admissible elds of stresses are
more dicult to construct than the kinematically admissible elds of displacements (at least in case
of plane problems); the stress boundary conditions lead to linear equations when stress functions
are used, while?in the displacement formulation?nodal displacements can be prescribed usually.
Despite these diculties, the stress-based FE model is attractive since the use of it together with
? Correspondence to: Zdzislaw Wi e
Vckowski, Department of Mechanics of Materials, Technical University of L├od├z, Al.
Politechniki 6, 93-590 L├od├z, Poland. E-mail: zwi@kmm-lx.p.lodz.pl
? Post-doctoral fellow. Permanently Lecturer, Department of Mechanics of Materials, Technical University of L├
od├z, Al.
Politechniki 6, 93-590 L├od├z, Poland
? Associate Professor
Д Graduate student
Contract=grant sponsor: Korea Science and Engineering Foundation
CCC 0029?5981/99/101505?21$17.50
Copyright ? 1999 John Wiley & Sons, Ltd.
Received 19 January 1998
Revised 19 June 1998
1506
Z. WI E
VCKOWSKI, S.-K. YOUN AND B.-S. MOON
the displacement FE model provides the lower and upper bounds for the strain energy, load limit,
etc., and makes possible to estimate the error for the approximate solution in an easy way.
As the rst papers concerning the stress approach to FEM, those written by Fraeijs de Veubeke,5
and Fraeijs de Veubeke and Zienkiewicz6 should be cited, where the analogy between plane problems and the bending plate problem has been shown. Using this analogy, Morley7 and Elias8
applied the stress approach to the Kirchho theory of plates. Watwood and Hartz9 solved the
plane problems of elasticity. Fraeijs de Veubeke and Hogge10 considered the dual approach to
the two-dimensional problem of heat ow using the triangular element with linear shape functions
for components of the heat-ux vector. The equilibrium model has been used more frequently in
the analysis of torsion of a prismatic bar,1 also in the case of plastic analysis,11; 12 and stationary
seepage problems.13; 14
During the last two decades, some newer concepts were introduced. Taylor and Zienkiewicz15
used the penalty method to satisfy the equilibrium equations?this method can be regarded as a
mixed approach. Gallagher, Heinrich and Sarigul solved plane problems of elasticity using the
Airy stress function and the concept of blending interpolation technique in order to facilitate the
satisfaction of the boundary conditions.16; 17 Vallabhan and Azene18 also solved the plane problems,
but used Lagrange multipliers to impose the boundary conditions. Other recent works by Gallagher
and co-workers19; 20 should also be cited here. It should be noted that the stress boundary conditions
can also be satised by the use of the penalty method (see e.g. Reference 1) which does not
increase the number of unknown variables?this is a relatively simple but approximate technique.
The stress-based formulation of FEM was utilized in the theory of plasticity by Rybicki and
Schmit,21 Gallagher and Dhalla,22 and Azene.23 The plane problems are solved by the
above-mentioned authors?they used the elastic?plastic constitutive matrix to calculate the stress
increments in the iterative procedure applied for each load increment. Belytschko and Hodge24
used the stress approach for the limit analysis of the plane stress problem.
The application of complementary variational principles to elastodynamics and elastic buckling
was considered by Tabarrok, Sodhi and Simpson in References 25?27.
Using dual formulations of FEM, WiVeckowski found lower and upper bounds of eective moduli
for an elastic, brous composite material,28; 29 and analysed its eective constitutive relations in
the case of elastic?plastic behaviour.30 The application of the stress-based FEM to plane friction
problems was outlined by the same author in Reference 31.
In the present paper, the Airy stress function is employed to construct the self-equilibrating
stress elds. The function is approximated by the Bogner?Fox?Schmit rectangular32 and Hsieh?
Clough?Tocher triangular33 elements. The Lagrange multiplier technique is used to full the traction boundary conditions?the concept of multi-point-constraints elements is developed to make
this technique easier to implement and use. The proposed iterative procedure for solving the
non-linear system of nite element equations for each load increment is analogous to the closestpoint-projection method used commonly in the displacement formulation of FEM. It is shown that
the equilibrium model of FEM can be implemented relatively easily in the standard nite element
program.
SETTING OF THE PROBLEM
The problem of equilibrium of an isotropic elastic?plastic body subjected to plane stress or plane
strain state is considered. Although these two-dimensional problems are considered, for the sake
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999)
STRESS-BASED FINITE ELEMENT ANALYSIS OF PLANE PLASTICITY PROBLEMS
1507
of brevity, most of the mathematical formulas will be written for the general three-dimensional
case throughout the paper.
Let denote an open bounded domain in R3 with regular boundary @
which consists of two
parts u and such that u ? =@
; u ? = ?. Let I denote the time interval, I ?[0; T ]; T ┐0.
The solution to the considered quasi-static problem: the elds of displacements, ui (t), and
stresses, ij (t), satisfy the following relations:
(i) The strain?displacement relations
ij = 12 (ui; j + uj; i )
on (1)
where ij denotes the tensor of small deformations, ui the displacement vector.
(ii) The equilibrium equations
ij; j + bi = 0
on (2)
where ij is the stress tensor, and bi the vector of body forces.
(iii) The Prandtl?Reuss constitutive relations (see e.g. References 34 and 35)
?
?ij = ?eij + ?pij
?
?
?eij = Cijkl ?kl
on ?
?
p
?ij (ij ? ij )60 ? ? B
(3)
where the dot indicates the time derivative, ije and ijp are the elastic and plastic parts of the
strain tensor, respectively, and Cijkl denotes the tensor of elastic compliance,
1
(4)
[(1 + ) ? (1 + ) ]
E
where is equal to 0 or 1 for the plane stress problem or the plane strain problem, respectively,
C =
B = {: f() = q ? 0 60}
(5)
is the set of plastically admissible stresses, where 0 denotes the plastic limit in the uniaxial
tension test, and q is the stress invariant dened as follows:
q
q = 32 sij sij ; sij = ij + pij ; p = ? 13 ii
(the Huber?von Mises yield condition is considered).
(iv) The boundary conditions
ui = Ui
on
u;
ij nj = ti
on
where Ui and ti are given displacements and tractions on
the unit vector outwardly normal to the surface @
.
(v) The initial conditions are taken as
ui = 0; ij = 0
for t = 0
u
and
,
respectively, and ni is
(6)
where t is the time variable.
In the plane stress problem, the following conditions hold 3i ? 0, 3 ? 0, b3 ? 0, while in case
of the plane strain problem 3 ? 0, 3i ? 0, b3 ? 0. It is assumed above that the Latin indices vary
from 1 to 3 while the Greek ones from 1 to 2.
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999)
1508
Z. WI E
VCKOWSKI, S.-K. YOUN AND B.-S. MOON
DUAL VARIATIONAL FORMULATIONS OF THE PROBLEM
Let V and W denote the following spaces:
V = {v ? [H 1 (
)]3 : vi = U? i (t) on
u}
3
W = {v ? [BD(
)] : vi = U? i (t) on
u}
where H 1 (
) is the Sobolev space (see e.g. Reference 36, and BD(
) is the space of functions
of bounded deformation (see References 37 and 38). Let P denote the set
P = { ? [L2 (Y )]9 : ij = ji ; (x) ? B(x) a:e: on }
where L2 (
) is the space of square integrable functions (Reference 36), B is the set of plastically
admissible stresses dened in equation (5).
From the constitutive relations (3), the following inequality can be derived:
(?ij ? Cijkl ?kl )(ij ? ij )60
? ? B
After integrating this relation over domain , and using the strain?displacement relations (1) and
the denition of set P, we get the inequality
Z
(Cijkl ?kl ? u?i; j )(ij ? ij ) dx┐0 ? ? P
After multiplying the equilibrium equations (2) by (vi ? u?i ) and integrating the result over domain
, and then using the Green formula, we obtain the equation expressing the principle of virtual
work
Z
Z
Z
ij (vi; j ? u?i; j ) dx =
bi (vi ? u?i ) dx +
ti (vi ? u?i ) ds ?v ? V
The problem can be set variationally as follows: Find (; u?) : I ? P О W such that a.e. on I
b(;
? ? ) ? (e(u?); ? ) ┐ 0
? ? P
(; e(v ? u?)) = l(v ? u?)
(7)
?v ? V
(8)
(0) = 0
where the following notation is used:
Z
Cijkl ij kl dx;
b(; ) =
[e(v)]ij =
1
(vi; j + vj; i );
2
(9)
Z
(; ) =
ij ij dx
Z
l(v) =
Z
bi vi dx +
ti vi ds
(10)
The variational formulation (7)?(9) is the starting point for the formulation of the problem in
displacements, in stresses, and in mixed forms. In order to set the problem in stresses, we introduce
the following set of statically admissible elds of stresses
Y = { ? [L2 (
)]9 : ij = ji ; ij; j + bi = 0 on ; ij nj = ti on
Copyright ? 1999 John Wiley & Sons, Ltd.
}
Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999)
STRESS-BASED FINITE ELEMENT ANALYSIS OF PLANE PLASTICITY PROBLEMS
1509
The equation of virtual work (8) is satised identically for each ? Y , and the following equation
is true
Z
Z
U? i (ij ? ij )nj ds ?u? ? W; ? ? Y
u?i; j (ij ? ij ) dy =
(e(u?); ? ) ?
u
Then the considered problem can be formulated in the dual form as follows: Find : I ? K such
that a.e. on I
b(;
? ? ) ┐ L(U? ; ? )
? ? K
(0) = 0
where
(11)
(12)
Z
U? i ij nj ds;
L(U? ; ) =
K =Y ?P
u
The quasi-static variational formulation in stresses for the plasticity problem has been given, for
example, in References 34 and 38.
FINITE ELEMENT SOLUTION
Time discretization
Let us consider a partition of the time interval I into N subintervals as 0; t1 ; : : : ; tn ; : : : ; tN =
T (0Аt1 А и и и Аtn А и и и АtN ). The time derivative is approximated by the backward dierence
scheme
(t
? n) ?
=
n n ? n?1
?
tn
tn ? tn?1
(13)
where n ? (tn ).
Substitution of (13) into (11) gives the semi-discrete (with respect to time) formulation:
For n = 1; 2; : : : ; N , nd n ? K such that
n = n?1 + n
n
n
n
(14)
n
b( ; ? ) ┐ L(U ; ? ) ? ? K
0
=0
(15)
(16)
n
The solution to the inequality (15), , minimizes the functional of complementary energy
Z
Z
1
1
Cijkl ij kl dx ?
Uin ij nj ds
(17)
() = b(; ) ? L(U n ; ) ?
2
2 u
on the set K ? K(tn ).
Construction of self-equilibrating stress elds
According to the denition of set Y , the elds of stresses?elements of this set?must satisfy
the equilibrium equations inside region and on its boundary . In the two-dimensional problems
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999)
1510
Z. WI E
VCKOWSKI, S.-K. YOUN AND B.-S. MOON
considered in the paper, the equilibrium equations (2) are satised identically if the stresses are
expressed by the Airy stress function, F, as follows:
= e e F; + ?
(18)
where e and are the permutation and Kronecker symbols, respectively, and ? is the potential
of volumetric forces such that b = ? ?; : The problem of fullling the stress boundary conditions
is described in the Appendix.
Space discretization
Stresses being the second derivatives of the main variable, the Airy stress function, require the
function to be approximated with the help of elements of C 1 class. Several types of elements
have been used in previous studies. The Bogner?Fox?Schmit rectangular element (see References
32 and 3) has been used by Gallagher and Dhalla,22 Vallabhan and Azene,18 Gallagher et al.,16
Sarigul and Gallagher,17 and Wieckowski.28; 29; 31 The Bell triangular element (e.g. Reference 3)
with 18 degrees of freedom has been used by Azene,23 and Vallabhan and Azene.18 Rybicki
and Schmit21 have used the rectangular element with polynomials of the 5th degree as shape
functions. Gallagher et al.,16 and Sarigul and Gallagher,17 have used rectangular elements with the
so-called ?blending function interpolants? in order to utilize the second derivatives as degrees of
freedom at boundary nodes. Watwood and Hartz9 used triangular elements with linearly distributed
self-equilibrating stresses. Unlike the above-mentioned works, where nodal parameters related to
the Airy stress function are basic unknown variables, they introduced the so-called ?generalized
stresses?. Then they applied some concepts of ?building block? in order to obtain equilibrium
between the elements. One of these concepts leads to the same results as would be obtained by
using the Hsieh?Clough?Tocher macro-element (see References 33 and 3) utilized in the present
paper.
In the present paper, the Bogner?Fox?Schmit rectangular and Hsieh?Clough?Tocher triangular
elements are applied to approximate the Airy stress function.
The Bogner?Fox?Schmit element (BFS, see References 32 and 3), shown in Figure 1, has 4
degrees of freedom at each corner, [F @F=@x @F=@y @2 F=@x@y]T . The Airy function is approximated
by a polynomial of the order less than 4 with respect to both variables x and y
P
F=
cpq xp yq
(19)
p; q63
In order to apply the Gauss quadrature, 16 (4 О 4) integration points are introduced.
Figure 1. The Bogner?Fox?Schmit element
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999)
STRESS-BASED FINITE ELEMENT ANALYSIS OF PLANE PLASTICITY PROBLEMS
1511
Figure 2. The Hsieh?Clough?Tocher macro-element
Using equations (18) and (19), we obtain normal and tangential stresses distributed linearly and
quadratically, respectively, along the element edges.
The Hsieh?Clough?Tocher element originally introduced in Reference 33 (see also References
1?3) is a macro-element, subdivided into three triangles (see Figure 2) inside each of them the
Airy stress function is a complete cubic polynomial
?
on K1
? F1 = c1 + c2 x + c3 y + и и и + c9 xy2 + c10 y3
2
3
(20)
F = F2 = d1 + d2 x + d3 y + и и и + d9 xy + d10 y
on K2
?
F3 = e1 + e2 x + e3 y + и и и + e9 xy2 + e10 y3
on K3
The element has 6 nodes: 3 at the corners and 3 in the mid-sides. 12 degrees of freedom are dened:
three degrees of freedom, [F @F=@x @F=@y]T , at each corner node, and one, the normal derivative,
@F=@n, at each mid-side node. The 30 coecients existing in the polynomial representation of
function F (20) are calculated by solving the system of the following equations: the continuity
conditions for function F and its rst derivatives, written for corner nodes P1 ; P2 and P3 (18
equations), and the central point, C (6 equations); the continuity conditions for normal derivatives
of function F, written for external mid-side nodes P4 ; P5 and P6 , and internal mid-side nodes Q1 ; Q2
and Q3 (6 equations).
In the present paper, the reduced HCT element is also utilized (see e.g. References 1 and 3).
In this case, the normal derivative at each mid-side node is replaced by the average of the two
normal derivatives at the ends of the corresponding side?this leads to reduction of numbers of
nodes and degrees of freedom from 6 to 3, and from 12 to 9, respectively.
The above system of algebraic equations is established for real positions of nodes. The mapping
from the reference space to the real one, used for instance in case of isoparametric elements,
is not applied here, because in general it does not preserve angular quantities, which leads to
destruction of C 1 -continuity between elements resulting in the loss of equilibrium for tangential
stresses.
The Hammer numerical integration procedure (see e.g. Reference 39) with 3 points for each
sub-triangle is applied for both types of elements: the original (HCT6) and reduced (HCT3) ones.
The normal stresses are distributed linearly along the edges of both types of HCT element, while
the tangential stresses are linear for the HCT6 element and constant for the HCT3 element.
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999)
1512
Z. WI E
VCKOWSKI, S.-K. YOUN AND B.-S. MOON
Let us consider the following discrete set of statically admissible elds of stresses, Yh ? Y ,
Yh = { ? [L2 (
)]4 : ? PK (Ki ) ?Ki ? Th ;
= ; ; + b = 0 on ; n = t on
}
where PK (Ki ) is the space of polynomials dened on element Ki according to the structure of
shape functions of the considered elements, Th denotes the triangulation of region .
The discrete formulation of the considered problem can be written as follows: For n = 1; 2; : : : ; N ,
nd hn ? Kh such that
hn = hn?1 + hn
b(hn ; ? hn ) ┐ L(U n ; ? hn ) ? ? Kh
h0
=0
(21)
(22)
(23)
where Kh = Yh ? P.
Using matrix notation, we can write the system of nite element equations for the equilibrium
model of FEM in the form
Z
HT U dx ? F = 0
(24)
analogous to that obtained for the displacement model, where
?
?
@2 N1 @2 N2
и
и
и
?
? @y2
@y2
?
?
Z
2
2
?
?
@ N1 @ N2
T T
?
H n U ds; H = ?
F=
и и и?
?
2
2
@x
?
? @x
u
?
? @2 N
2
@ N2
1
?
иии
?
@x@y @x@y
? ?
? x ?
Ux
n
0 ny
U = y ; n = x
; U=
0 ny nx
Uy
? ?
xy
(25)
Equation (24) can be derived from (12) or (17) after representation of the Airy stress function by
nodal parameters a ? [a1 a2 и и и]T and shape functions N ? [N1 N2 и и и]T , F = Na.
Iterative solution
For each time increment, tn , the following iterative procedure is employed:
(i) Initialization of increments of deformations and stresses (at integration points)
ijn(0) = 0;
hn(0) = 0
(ii) For the ith iteration, calculation of the successive estimation of hn(i) ? Yh , solving the linear
problem
b(hn(i) ? hn(i?1) ; ? hn (i) ) = ?(n(i?1) ; ? hn(i) ) + L(U n ; ? hn(i) ) ? ? Yh
Copyright ? 1999 John Wiley & Sons, Ltd.
(26)
Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999)
STRESS-BASED FINITE ELEMENT ANALYSIS OF PLANE PLASTICITY PROBLEMS
1513
(iii) Calculation of increment of strains, n(i) , for each integration point, from the relations
ijn(i) = ijen(i) + ijpn(i)
n(i)
ijen(i) = Cijkl hkl
(
nij if ┐0
pn(i)
ij =
0
if 60
where
nij =
3 sij
2 q
= eijp(i) nij
n(i)
eijp(i) = ijpn(i?1) + Cijkl (hkl
? [P(hn(i) )]kl )
(27)
P is the operator of projection onto the yield surface.
(iv) Checking the accuracy of the solution. If the assumed accuracy criterion is satised the
calculation process can be stopped for the nth time increment, otherwise the next iteration
should be performed starting from point (ii).
The matrix of the linear system of nite element equations (26) is constant during the entire
calculation process.
The operator of stress projection onto the yield surface is dened as follows:
(ij ? [P()]ij )Cijkl (kl ? [P()]kl )6(ij ? ij )Cijkl (kl ? kl ) ? ? B
(28)
and has easy form in case of the plane strain problem
1
0
[P()]ij = kk ij + sij
3
q
while in the case of the plane stress problem the constrained minimization problem (28) has to
be solved numerically. Following the concept described in Reference 40, the Newton?Raphson
technique has been used to nd the minimum point of the augmented function
? 0 )
L(;
) = ( ? )C ( ? )?(q()
where is the Lagrange multiplier. In the trivial case when ij = 0, the relation [P()]ij = 0 should
be applied, in order to avoid the problem of non-uniqueness of projection onto the yield surface.
The procedure for calculation of the increment of plastic strains is illustrated in Figure 3.
The procedure is analogous to that commonly applied in the standard displacement nite element
method, known as the closest-point-projection method or backward Euler method (see e.g. References 41 and 42). The above procedure and the CPP method give the same results in the case of
boundary value problems, where stress elds are homogeneous.
Equation (26) has the following representation in matrix notation
Z
Z
T
n(i)
n(i?1)
H C H dx (a ? a
)= ?
HT ?n(i?1) dx + Fn
where terms of the matrix of elastic compliance, C, are dened in equation (4).
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999)
1514
Z. WI E
VCKOWSKI, S.-K. YOUN AND B.-S. MOON
Figure 3. Calculation of increment of plastic strains
Figure 4. Bending of rectangular plate, dimensions in centimeters
NUMERICAL EXAMPLES
Two numerical examples are considered: the bending of a rectangular plate, and the tension of a
perforated plate. The application of the three described elements: BFS, HCT6 and HCT3 are shown
in the case of the rst problem, while in the second one, only the triangular elements: HCT6 and
HCT3 are utilized. The results are compared with those obtained by the displacement model of
the nite element method with the help of 6-node isoparametric triangular elements, for which all
stress components are approximated by linear functions as in the case of the HCT6 element.
Bending of rectangular plate
The plane strain problem of a rectangular plate is considered. The plate is subjected to the
self-equilibrating loading according to Figure 4. Calculations have been made with the following
material data (steel): Young?s modulus E = 2 О 1011 Pa, Poisson?s ratio = 0 и3, plastic limit 0 =
2и5 О 108 Pa.
To stop the iteration process, the following convergence criterion has been checked:
|q(i) ? q(i?1) |
60и01
|q(i) |
R
where q(i) = HT U(i) dx (see notation in equation (25)) in the case of equilibrium model
of the nite element method, and q(i) is the vector of degrees of freedom calculated in the ith
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999)
1515
STRESS-BASED FINITE ELEMENT ANALYSIS OF PLANE PLASTICITY PROBLEMS
Table I. Lower and upper bounds of load limit
Element size (m)
0и05
0и025
0и0125
Element
type
pmax
(MPa)
Dierence
(per cent)
pmax
(MPa)
Dierence
(per cent)
pmax
(MPa)
Dierence
(per cent)
BFS
HCT6
HCT3
TRI6
102и4
106и8
100и4
117и4
12и777
9и029
14и480
?
109и2
108и4
106и6
114и0
4и211
4и912
6и491
?
111и6
111и2
110и0
112и6
0и888
1и243
2и309
?
iteration in the case of displacement-based nite element method, | и | the Euclidean length of a
vector.
The load increment has been set as p = 1 О 106 Pa, and reduced to p = 0и2 О 106 Pa for the
load level close to the limit value. Calculations have been made for three regular meshes with
elements having equal horizontal and vertical dimensions: 0и05; 0и025 and 0и0125 m. Because of
the symmetry, only the left part of the plate has been analysed.
The maximum values of load p, for which the required accuracy of the solution has been
achieved, are shown in Table I. The abbreviation TRI6 is used in the table for the 6-node triangular
element used in the displacement model of the nite element method. The values obtained by the
stress-based model are the lower bounds for the limit load, while those obtained by the standard
displacement model of FEM are the upper bounds. The table contains also the relative dierence
between the bounds, calculated as follows:
kin
st
? pmax
pmax
kin
pmax
kin
st
where pmax
and pmax
denote bounds obtained by the displacement and equilibrium models of FEM,
respectively. In the case of the nest mesh the dierences do not exceed 2и4 per cent, which is
satisfactory from the engineering viewpoint. The slight dierence between the results obtained for
BFS and HCT6 elements is observed. The HCT3 element provides a lower value for pmax due to
the fact that the tangential stresses are constant along its boundary, while in the case of HCT6
and BFS elements they are linear and quadratic functions, respectively.
The use of the method of Lagrange multiplier gives the possibility to calculate displacements
for the boundary where stresses are prescribed (see the Appendix). Assuming that points A and
B are xed in the vertical direction, the mean vertical displacement for the segment of boundary
L (see Figure 4), located near the symmetry axis, has been calculated, and compared with the
kinematically admissible displacement of the corresponding node. The comparison is shown in
Figure 5 in the form of load?deection diagrams. The relative dierences between kinematically
and statically admissible solutions are also shown in the gure for all elements used: BCT, HCT6
and HCT3. The dierences between displacements calculated by the two methods used are small
at the beginning of the load?deection paths, and grow rapidly when the load approaches the limit.
The BFS and HCT6 elements give similar results. The HCT3 element produces larger values of
displacements, but for the nest mesh, the results are quite satisfactory, as the dierences are
smaller than 1 per cent for the initial part of the path.
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999)
1516
Z. WI E
VCKOWSKI, S.-K. YOUN AND B.-S. MOON
Figure 5. Load?deection diagrams
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999)
STRESS-BASED FINITE ELEMENT ANALYSIS OF PLANE PLASTICITY PROBLEMS
1517
Figure 6. Plastic zones
Table II. Computation time in seconds (Silicon Graphics
MIPS R5000 Indy A2)
Element
type
BFS
HCT6
HCT3
TRI6
Element size (m)
0и05
0и025
0и0125
5и263
6и406
6и135
1и031
19и438
25и344
19и666
3и609
85и324
102и825
88и375
14и427
The comparison of plastic zones obtained with the help of the three self-equilibrating elements
with the kinematically admissible solution is shown in Figure 6 for three levels of load p. The
similarity between the results can be observed.
It is worth mentioning that the stress-based model appears to require approximately twice the
number of iterations and a larger number of integration points than that used in the displacement
model. The computation times for the self-equilibrating elements and the TRI6 element are compared in Table II for the rst 30 load increments. The comparison shows that the stress-based
model of FEM is more time consuming than the displacement-based one.
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999)
1518
Z. WI E
VCKOWSKI, S.-K. YOUN AND B.-S. MOON
Figure 7. Geometry of perforated plate, dimensions in millimeters
Figure 8. Exemplary meshes
Tension of perforated plate
The plane stress problem of the perforated plate shown in Figure 7 is considered. Unlike the
previous example, the load on the plate is displacement-controlled?it is assumed that the horizontal
displacements, U1 = ▒U , on the vertical sides of the plate are given. The considered material data
corresponds to aluminium: Young?s modulus E = 7 О 1010 Pa, Poisson?s ratio = 0и2, plastic limit
0 = 2и34 О 108 Pa.
Three triangular element meshes, shown in Figure 8, are used in the analysis. The values of
maximum element sizes for these meshes are as follows: (a) 3и86 mm, (b) 2и09 mm, (c) 1и12 mm.
As the plate is doubly symmetric, only the right upper quarter of its area is analysed.
The increment of given displacement is set as U = 2 О 10?6 m. The same accuracy criteria are
applied as in the previous example.
The results of calculations are shown in the form of P?U paths presented in Figure 9, where
P denotes the resultant force calculated for half height of the plate. As in the previous example,
the smaller values of forces, P, are obtained for the HCT3 element than for the HCT6 one. The
dierences between the load limits, calculated by the two nite element models, are not greater
than 1и5 per cent for the nest mesh.
The growth of plastic zones is presented in Figure 10, where their shapes are compared for
three values of displacement U . Similar plastic zone shapes are obtained for the three types of
elements used.
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999)
STRESS-BASED FINITE ELEMENT ANALYSIS OF PLANE PLASTICITY PROBLEMS
1519
Figure 9. P?U paths for perforated plate
CONCLUSIONS
The equilibrium model of the nite element method is applied to plane, quasi-static problems
of plasticity theory. The statically admissible elds of stresses are constructed by means of the
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999)
1520
Z. WI E
VCKOWSKI, S.-K. YOUN AND B.-S. MOON
Figure 10. Plastic zones for perforated plate
Airy stress function. Three kinds of elements are utilized to approximate the stress function: the
Bogner?Fox?Schmit rectangular element, and the Hsieh?Clough?Tocher triangular element in its
original and reduced versions.
The stress boundary conditions are fullled using the Lagrange multiplier method, allowing the
evaluation of the displacements on the boundaries with prescribed tractions. The concept of multipoint-constraints elements is proposed; this leads to an easier implementation of the method as
well as the easier use of a computer program.
The proposed iterative procedure, for solution of the non-linear system of equations for each load
increment, is analogous to the closest-point-projection method routinely utilized in the displacement
nite element method. The used procedure and the CPP method give the same results in case of
boundary value problems, where stress elds are homogeneous. The constant matrix of the set of
algebraic equations is applied in the procedure.
The equilibrium model of the nite element method provides a lower bound for the load limit
which is important from the design perspective. However, some eort is to be paid?this model
is more dicult to implement than its displacement counterpart, because of necessity of using
elements of C 1 class and multi-point-constraints imposing the appropriate boundary conditions.
The equilibrium model usually need the longer computation time than the displacement one when
the similar numbers of elements are applied to the specic problem in both the models of FEM.
The concepts, described in the present paper, somewhat soften these two drawbacks, as they make
it possible to implement the stress-based model in the standard nite element program without
special extensions.
The numerical examples show the convergence of the equilibrium model and good agreement
with the results obtained by the use of the standard displacement nite element method. Small
dierences between the lower and the upper bounds are obtained for load limits. Very good
agreement for displacements and shapes of plastic regions are also observed.
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999)
1521
STRESS-BASED FINITE ELEMENT ANALYSIS OF PLANE PLASTICITY PROBLEMS
Table III. Characteristics of elements
Distribution of stresses
Element
type
Normal
BFS
HCT6
HCT3
Linear
Linear
Linear
Tangential
Number of Lagrange
multipliers for
MPC element
Number of
d.o.f. in MPC
element
Quadratic
Linear
Constant
13
11
9
5
4
3
APPENDIX
Imposition of stress boundary conditions
Using notation shown in equation (25), we can express the boundary conditions in the form of
the matrix equation
nHa = t
(29)
where t = [tx ty ]T is the vector of given tractions. Equation (29) indicates that the multi-point
constraints (MPC) problem is to be solved unlike in the displacement formulation of the nite
element method, where, for simplest cases, the nodal displacements can be prescribed directly. In
the present paper, equation (29) is satised with the help of the method of Lagrange multiplier. This
technique leads to increasing the number of unknowns in the nal system of algebraic equations,
but enables calculation of displacements for points belonging to boundary , where tractions are
given. This follows from the fact that Lagrange multipliers are interpreted as displacements of
boundary points when used to impose equation (29).
Using (29), we obtain the same distribution of normal stresses and dierent distribution of
tangential stresses for the three elements utilized in the paper (see Table III). In the considered
nite element model, only the boundary conditions of the forms given in Table III can be satised
identically. In other cases, they should be approximated by functions given in the table.
In order to facilitate the implementation of the Lagrange multiplier technique as well as handling
the boundary conditions by a user of the computer program, the concept of MPC-elements is used.
The concept is described for HCT3 element. Using the local coordinate system (n; s) for the edge
of the element, shown in Figure 11, the normal and tangential stresses can be written as functions
of s, and are dependent on degrees of freedom dened at nodes P1 and P2 only
6
2
(2a1 + la3 ? 2a4 + la6 ) s ? 2 (3a1 + 2la3 ? 3a4 + la6 )
l3
l
1
T (s) = (a2 ? a5 )
l
N (s) =
(30)
where a1 ; : : : ; a6 are terms of the vector of degrees of freedom dened at nodes P1 and P2
aT = [F(P1 ) @F=@n(P1 ) @F=@s(P1 ) F(P2 ) @F=@n(P2 ) @F=@s(P2 )]T
If normal and tangential stresses are prescribed by means of constants A, B and C as follows:
N (s) = A s + B
T (s) = C
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999)
1522
Z. WI E
VCKOWSKI, S.-K. YOUN AND B.-S. MOON
Figure 11. MPC-element for HCT3 element
it follows from (30) that the boundary conditions are satised for arbitrary s if the following
system of equations is fullled:
? ?a ?
?
?
1?
0
6l?2 ?12l?3
0
6l?2
12l?3
?
? ? ?
?
? a2 ?
? ?A?
?
?
?1
?2
?1 ?
? ?6l?2
(31)
0
?4l
6l
0
?2l
.
? ? . ? = ?B?
?
. ?
?
C
?
?
?1
?1
? ?
0
l
0
0
?l
0
a6
or briey
Ga ? g = 0
(32)
The HCT3 element requires to introduce 3 Lagrange multipliers for each side located on in
order to satisfy equation (32). Let us dene MPC element corresponding to HCT3 element. Let
it consist of three nodes P1 , P2 , and the additional one, P3 , at which the auxiliary three degrees of
freedom are dened?3 Lagrange multipliers, associated with equations (31) (see Figure 11). Let
us dene the vector of degrees of freedom for this element, ae , and the following element square
matrix, ke , and the vector of right-hand-side terms, f e ,
?
?
F(P1 ) ?
?
?
?
?
?
?
?
?
?
?
?
@F
?
?
?
?
(P
)
1
?
? ?
?
?
?
@n
?
?
0?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
@F
0?
?
?
?
?
?
?
?
?
)
(P
1
?
?
?
?
?
?
?
?
@s
0
?
?
?
?
?
?
?
? F(P ) ?
? ?
?
?
?
?
2 ?
?
?0?
?
T
0G
e
e
e
; f = 0
; k =
a =
G 0
? @F
? ?
?
?
?
?
0?
?
?
?
(P2 ) ?
?
?
?
?
?
?
?
?
@n
?
?
?
A?
?
?
?
?
?
?
?
?
?
?
?
?
@F
?
?
?
?
B
?
?
?
?
(P
)
?
?
?
2
?
? ?
?
?
?
?
@s
C
?
?
?
?
?
?
?
?
1
?
?
?
?
?
2 ?
?
?
?
?
?
?
3
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999)
STRESS-BASED FINITE ELEMENT ANALYSIS OF PLANE PLASTICITY PROBLEMS
1523
Then treating matrices ke and f e of MPC element, during the assembly process for global matrices
K and F; in the same way as exibility matrices and right-hand-side vectors of triangular elements,
we impose the boundary conditions (32). Before assembling, matrix ke should be transformed to
the global co-ordinate system
kge = TT ke T
where
?
T
T= ? 0
0
0
T
0
?
0
0? ;
I
?
1
T = ? 0
0
?
0
ny ? ;
nx
0
nx
?ny
?
1
I = ?0
0
0
1
0
?
0
0?
1
The element matrices and vectors for MPC elements associated with HCT6 and BFS elements
can be dened in a similar way. The numbers of Lagrange multipliers and degrees of freedom of
all MPC elements used are summarized in Table III.
Lagrange multipliers N0 and T0 (N0 = 2 , T0 = 3 in case of HCT3 element) corresponding with
equations written for terms containing s0 in expressions for N (s) and T (s) can be interpreted as
Z l
Z l
uN ds and
uT ds
0
0
respectively, so the average normal and tangential displacements for each segment can be calculated
from the equations
av
=
uN
N0
l
and
uTav =
T0
l
The above interpretation of multiplier N0 is utilized in the rst numerical example.
The use of Lagrange multiplier technique means that zero terms appear on the diagonal of the
element matrix ke . Therefore, the pivoting or the proper ordering of degrees of freedom for the
global system of equations should be used. The easiest way is to group all additional nodes, where
Lagrange multipliers are degrees of freedom, at the end of the list of nodes, but this measure leads
to uneconomical growth of length of prole or half-width of the equations matrix. In the present
paper, a simple solution, consisting of two steps, is used. In the rst step, the well-known Cuthill?
McKee43 algorithm is applied to renumber only basic nodes related to variables dening the Airy
stress function. Then all numbers of the additional nodes are inserted to the list of all nodes in
such a way that the number of each additional node is set to n + 1, where n denotes the largest
number of a basic node belonging to the MPC element related to the considered additional node.
ACKNOWLEDGEMENTS
The rst author wishes to thank the Korea Science and Engineering Foundation for the nancial
support in the framework of the Post-doctoral Fellowship Program for Foreign Researchers.
REFERENCES
1. O. C. Zienkiewicz and R. L. Taylor, The Finite Element Method, McGraw-Hill, London, 1989.
2. R. H. Gallagher, Finite Element Fundamentals, Prentice-Hall, Englewood Clis, N.J., 1975.
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999)
1524
Z. WI E
VCKOWSKI, S.-K. YOUN AND B.-S. MOON
3. P. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, 1978.
4. M. Kleiber, Incremental Finite Element Modelling in Solid Mechanics, PWN=Ellis Horwood, Warszawa, 1989.
5. B. Fraeijs de Veubeke, ?Displacement and equilibrium models in the nite element method?, in O. C. Zienkiewicz and
G. S. Holister (eds.), Stress Analysis, Wiley, New York, 1965.
6. B. Fraeijs de Veubeke and O. C. Zienkiewicz, ?Strain-energy bounds in nite-element analysis by slab analogy?,
J. Strain Anal., 2, 265 ?271 (1967).
7. L. S. D. Morley, ?The triangular equilibrium element in the solution of plate bending problems?, Aero Quart., 19,
149?169 (1968).
8. Z. M. Elias, ?Duality in nite element method?, Proc. ASCE, J. Engng. Mech. Div., EM4, 94, 931?946 (1968).
9. V. B. Watwood and B. J. Hartz, ?An equilibrium stress eld model for nite element solution of two-dimensional
elastostatic problems?, Int. J. Solids Struct., 4, 857?873 (1968).
10. B. Fraeijs de Veubeke and M. Hogge, ?Dual analysis for heat conduction problems by nite elements?, Int. J. Numer.
Meth. Engng., 5, 65?82 (1972).
11. R. Glowinski, J. L. Lions and R. Tremolieres, Analyse Numerique des Inequations Variationelles, Dunod, Paris, 1976.
12. Z. Wi eVckowski, ?Dual nite element analysis for plasticity-friction torsion of composite bar?, Int. J. Numer. Meth.
Engng., 38, 1901?1916 (1995).
13. J. J. Connor and C. A. Brebbia, Finite Element Technique for Fluid Flow, Newnes-Butterworths, London, 1976.
14. A. R. S. Ponter, ?The application of dual minimum theorems to the nite element solution of potential problems with
special reference to seepage?, Int. J. Numer. Meth. Engng., 4, 85?93 (1972).
15. R. L. Taylor and O. C. Zienkiewicz, ?Complementary energy with penalty functions in nite element analysis?, in
R. Glowinski (ed.), Energy Methods in Finite Element Analysis, Wiley, New York, 1979.
16. R. H. Gallagher, J. C. Heinrich and N. Sarigul, ?Complementary energy revisited?, in Atluri, Gallagher and Zienkiewicz
(eds.), Hybrid and Mixed Finite Element Methods, Wiley, Chichester, 1983.
17. N. Sarigul and R. H. Gallagher, ?Assumed stress function nite element method: two-dimensional elasticity?, Int. J.
Numer. Meth. Engng., 28, 1577?1598 (1989).
18. C. V. G. Vallabhan and M. Azene, ?A nite element model for plane elasticity problems using the complementary
energy theorem?, Int. J. Numer. Meth. Engng., 18, 291?309 (1982).
19. R. H. Gallagher, ?Can and should the force method be revisited??, in S. N. Atluri, G. Yagawa and T. A. Cruse (eds.),
Proc. Int. Conf. Comp. Engng. Sci:; July 30?August 3; 1995; Hawaii; USA, Springer, Berlin, 1995, pp. 2?7.
20. C. Balakrishna, J. H. Kane and R. H. Gallagher, ?Assumed stress function nite element method: Towards unication?,
in S. N. Atluri, G. Yagawa and T. A. Cruse (eds.), Proc. Int. Conf. Comp. Engng. Sci.; July 30?August 3, 1995,
Hawaii, USA, Springer, Berlin, 1995, pp. 1578?1583.
21. E. F. Rybicki and L. A. Schmit, ?An incremental complementary energy method of nonlinear stress analysis?, AIAA
J., 8, 1805?1812 (1970).
22. R. H. Gallagher and A. K. Dhalla, ?Direct exibility nite element elastoplastic analysis?, Proc. Int. Conf. on Structural
Mechanics in Nucl. Reac. Tech., 1972.
23. M. Azene, ?A nite element complementary energy formulation for plane elastoplastic stress analysis?, Ph.D.
dissertation, Dept. of Civil Engineering, Texas Tech. Univ., Lubbock, Texas, May 1979.
24. T. Belytschko and P. G. Hodge, ?Plane stress limit analysis by nite elements?, Proc. ASCE, J. Engng. Mech. Div.,
EM4, 96, 931?944 (1970).
25. B. Tabarrok and D. Sodhi, ?On the generalization of stress function procedure for the dynamic analysis of plates?, Int.
J. Numer. Meth. Engng., 5, 532?542 (1972).
26. B. Tabarrok, ?Complementary variational principles in elastodynamics?, Comput. Struct., 19, 239?246 (1984).
27. B. Tabarrok and A. Simpson, ?An equilibrium nite element model for buckling analysis of plates?, Int. J. Numer.
Meth. Engng., 11, 1733?1751 (1977).
28. Z. Wi eVckowski, ?Duality in nite element method and its applications to some linear and non-linear problems of
mechanics of composite materials?, in Polish, Ph.D. dissertation, Institute of Civil Engineering, Technical University
of Lodz, Poland, February 1987.
29. Z. Wi eVckowski, ?Dual nite element methods in mechanics of composite materials?, J. Theoret. Appl. Mech., 33,
233?252 (1995).
30. Z. Wi eVckowski, ?Dual nite element methods in homogenisation for elastic-plastic brous composite material?, in
C. K. Choi, C. B. Yun and H. G. Kwak (eds.), Proc. 7th Int. Conf. on Computing in Civil and Building Engineering,
Seoul, Korea, 19?21 August 1997, pp. 1051?1056.
31. Z. Wi eVckowski, ?Dual nite element method in friction problems?, J. Theoret. Appl. Mech., 30, 535?543 (1992).
32. F. K. Bogner, R. L. Fox and L. A. Schmit, ?The generation of interelement compatible stiness and mass matrices by
the use of interpolation formulas?, Proc. Conf. Matrix Methods in Struct. Mech., October 1956, Wright Patterson
A. F. Base, Ohio, 1965.
33. R. W. Clough, J. L. Tocher, ?Element stiness matrices for analysis of plate bending?, Proc. Conf. Matrix Methods
in Struct. Mech., October 1956, Wright Patterson A. F. Base, Ohio, 1965.
34. G. Duvaut and J. L. Lions, Les inequations en mecanique et en physique, Dunod, Paris, 1972.
35. R. Hill, ?A comparative study of some variational principles in the theory of plasticity?, J. Appl. Mech., 17, 64?70
(1950).
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999)
STRESS-BASED FINITE ELEMENT ANALYSIS OF PLANE PLASTICITY PROBLEMS
1525
36. R. A. Adams, Sobolev Spaces, Academic Press, New York, 1975.
37. R. Temam and G. Strang, ?Functions of Bounded Deformation?, Arch. Rat. Mech. Anal., 75, 7?21 (1980).
38. C. Johnson, ?Mathematical and Numerical Analysis of Some Problems of Plasticity?, in Telega, J. J. (ed.), Functional
Analysis Methods in the Theory of Plasticity, Ossolineum, Wrocl aw, 1981, pp. 71?109 (in Polish).
39. G. Dhatt and G. Touzot, The Finite Element Method Displayed, Wiley, Chichester, 1984.
40. A. Samuelsson and M. Froier, ?Finite elements in plasticity?a variational inequality approach?, in J. R. Whiteman
(ed.), MAFELAP 1978 (III ), Academic Press, London, 1979.
41. J. C. Simo and R. L. Taylor, ?Consistent tangent operators for rate-independent elastoplasticity?, Comp. Meth. Appl.
Mech. Engng., 48, 101?118 (1985).
42. M. A. Criseld, Non-linear Finite Element Analysis of Solids and Structures, Volume 1: Essentials, Wiley, Chichester,
1991.
43. A. George and J. W. Liu, Computer Solution of Large Sparse Positive Denite System, Prentice-Hall, Englewood
Clis, N.J., 1981.
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 1505?1525 (1999)
Документ
Категория
Без категории
Просмотров
1
Размер файла
288 Кб
Теги
development, second, forma, practical, approach, generation, general, order, theoretical, times, transiente, heat, operatora, transfer, including, algorithms, accurate, useful, nth
1/--страниц
Пожаловаться на содержимое документа