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



код для вставкиСкачать
39,619-633 (1996)
National Laborarory for Structural Analysis of Industrial Equipment, Dalian University of Technology,
Dalian 116023, P.R. China
A finite element method for gradient elasto-plastic continuum in which the yield strength of strain
hardening/softeningmaterials not only depends on the effective plastic strain but also on its Laplacian is
presented. The consistent integration algorithm to update the stress and the internal state variable at
integration points and the consistent compliance matrix for the gradient plasticity are formulated in the
non-local sense.
The methodology to derive the finite element formulation for the gradient plasticity at large strains
presented in this paper is applicable to general finite element analysis; the formulation in the context of the
two-dimensional four-noded mixed finite element with one integration point and mean von Mises yield
criterion is particularly derived.
Numerical examples are tested to demonstrate the capability and performance of the present finite
element method at large strain in solving for the strain localization problem.
KEY WORDS gradient plasticity; strain localization; mixed strain element method; large strain; consistent algorithm
Localization of deformation into narrow bands of intense straining caused by strain softening is
a characteristic feature of plastic deformation. It has been experimentally observed in many
engineering materials such as metals, polymers, concrete and geological materials. Localization
phenomena are often associated with a significant reduction of the load-carrying capacity of the
structures; hence, the onset of localization is naturally considered as the point of threshold to
predict rupture of engineering structures. Over the last decade considerable effort has been
devoted to obtain a comprehensive understanding of the problem and to describe the behaviour
quantitatively. Numerous attempts to simulate the behaviour using classical plastic continuum
theories have been unsatisfactory.The principal defect of the classical plastic continuum models is
that they do not incorporate an internal length scale or higher-order continuum structure and,
therefore, the governing equations under quasi-static loading lose ellipticity. Consequently, as
strain softening behaviour is incorporated into a computational model the ill-posedness of the
problem will lead to pathologically mesh-dependent results. Furthermore, the energy dissipated
at strain softening is incorrectly predicted to be zero, and the finite element solutions converge to
incorrect, physically meaningless ones as the element mesh is refined.'
*Permanent address: M.S.M.Department, University of Liege, Quai Banning 6,4000,Liege, Belgium
CCC 0029-5981 /96/040619- 15
0 1996 by John Wiley & Sons, Ltd.
Received 22 August 1994
Revised 29 January 1995
To remedy the situation generalized continuum theories, in which internal length scales or
high-order continuum structures are introduced as regularization mechanisms for loss of ellipticity, should be adopted. Among those are gradient continuum theories,'-' Cosserat continuum
theory" and rate-dependent continua.'' In the present work we will utilize the gradient continuum approach. Many efforts have been devoted to developing the models using the gradient
continuum theories. Bazant presented the non-local gradient formulation with a Laplacian of
strain2 and applied the formulationin solving a one-dimensionallocalization p r ~ b l e mAifantis4*
show that the introduction of the second gradient of deformation into the expression of the flow
stress for plastic materials can preserve ellipticity of the governing equations. Lasry and Beiytschko6 introduced additional higher-order terms in the strain displacement relation for the
strain softeningportions of the domains for the proposed non-local formulations. Muhlhaus and
Aifantis' proposed their non-local gradient model on the basis of a variational principle for
gradient plasticity, in which higher-order spatial gradients of the effective plastic strain are
included in the yield condition.
Recently De Borst and his c o - w o r k e r ~presented
the formulations and algorithms for the
gradient plastic continuum in a finite element context and in these the yield strength depends not
only on the effective plastic strain but also on its Laplacian. It makes the presented formulations
and algorithms able to solve two- or three-dimensional boundary value problems. As the
Laplacian of the effective plastic strain is involved in the yield function and the increment of the
plastic multiplier cannot be obtained at a local (integration point) level, in contrast to the
conventional approach in computational plasticity, the plastic multiplier is taken as an independent global variable. The consistency condition of the yield function is enforced only in a weak
form; in addition, it is not satisfied at each iteration but only at the end of a load step.
With the use of the gradient plastic model including Laplacian of the effective plastic strain in
the yield condition, a finite element method for gradient elasto-plastic continuum is presented in
this work. The consistent integration algorithm to update the stress and the internal state variable
at integration points and the consistent compliance matrix for the gradient plasticity are
formulated in the non-local sense. Even though the yield condition at each integration point is no
longer local, it is satisfied in a pointwise fashion and at each iteration of a load step. The plastic
multipliers are internal state variables locally defined at each integration point. The Laplacian of
the effective plastic strain at an integration point is evaluated from the derivatives of the
polynomial, which is interpolated by using the values of the plastic multipliers at neighbouring
integration points. As the plastic multipliers are not taken as independent global unknowns there
is no requirement on C'-continuous shape functions for the interpolation of the plastic multiplier.
The methodology to derive the finite element formulations for the gradient plasticity at large
strains presented in this paper is applicable to general finite element analysis. Even the formulations in this paper are particularly derived in the context of the two-dimensional four-noded
mixed finite element with one integration point and von Mises yield criterion," the methodology
can be applied to general cases with the use of displacement-based finite elements, more general
yield conditions and non-local softening law such as that for quasi-brittle materials exhibiting
microcracking and internal friction.
We start with the rate constitutive equation for elasto-plastic materials:
where 8, %, 8, are total, elastic and plastic strain rates, respectively. D, is the matrix of elastic
moduli. The rate of plastic strain vector can be expressed as
where the rate of plastic multiplier
associated plasticity is defined by
= in
is a non-negative scalar and the flow vector n for the
n = af/aa
Here the yield functionf for the gradient plasticity is non-local and can be assumed in the form
f =f(a, E P , V%P)
The assumed effectiveplastic strain Ep can be calculated by the integration of incremental effective
plastic strain AEP:
where AEp is defined by
As the von Mises yield criterion is considered, the flow vector n takes the form
with s being the deviatoric stress vector and 6 the effective stress defined by
5 = (3STS)l’Z
The non-local von Mises yield surface can simply be written as
f (a,E P , VZEP) = @(a) - ay(EP,VZEP) = 0
Here ay is the current yield stress and is assumed as piecewise linear function of EP as well as
a linear function of the Laplacian of the effective plastic strain:
ay = uyo
+ h,(EP)P
- cV%P
where ayois the initial yield strength, h, is the hardening/softening modulus and is a function of
Ep, c is a non-local material parameter defined by
c = af/avZEp
According to equations (2) and (6)-(8) it is given that
. .
1=g p
The variation of the yield function with time can then be written as
f = 8 - h,X
Consider a load step from time t to time t
+ At. Integration of equation (1) with the time gives
- AADent+Ar
Gf+At =a;+~t
where the trial elastic stress at time t
+ At is defined by
= De(&r+Af-
The flow vector at the point eE in the stress space is
where sEis the deviatoric component of 8 and
5 E = (3SETSE)1/2
It can be verified for plane strain, axisymmetric and three-dimensional solid that
nf+At = nt+Ar
Pre-multiplying by nT both sides of equation (14) and using equation (18) gives
4 + a t = 8,"+~t KEAJ
where the effective elastic stiffness (Gis the shear modulus) is given by
KE= nTD,n = 3G
Substitution of the variation 8 evaluated from equation (19) into equation (13) results in
- (KE + hp)X + cV2A
It is observed that, unlike in classical (local) computational plasticity, incremental plastic
multiplier cannot be determined at a local (integration point) level.
Consider a typical time step [t, t + At] for an arbitrary material point i in plastic yielding. The
yield condition at the current iteration k can be written in an iterative manner as
5 . k =fi,k-
SA = A , k -
- ( K E+ hp)id(AAi) + CV26(AAi) = O
As shown in equation (22) the yield condition at an integration point is related to not only the
local plastic variable but also to that in the vicinity of the integration point.
To evaluate V2 A&, we consider the neighbouring points around the integration point i. Let
N i be the number of the neighbouring integration points, including the point i itself. As the plastic
multipliers at these N iintegration points are used to compute V2AJi, we are able to express
V2 A l i in terms of Adj ( jE N i ) in the form
To determine the coefficients gij, a complete second-order polynomial in two co-ordinates is
assumed to represent the function of plastic multiplier around point i, i.e.
A 1 = aTv
a = [al a2 a3 a4 a5 aJ
v = [ 1 x y x2 x y y2-J‘
The coefficient vector a is obtained by minimization of
Error =
(AA - Ahj)’
j= 1
with AA = [A& . . . A&JT denoting the values of Al, at the neighbouring integration points of
the integration point i, and we get the result of vector a in the form
Ga = M A A
C in equation (30) stands for Cy:,.
Matrix M is defined by
M = [vI
. . . vj . . . vNi]
with vj being the vector v at jth neighbouring integration point.
Hence, the coefficientsa can be expressed in terms of the values of the AA of the neighbouring
integration points
a = G - l M A A = G-l
j= 1
Though V2 AA can be evaluated, with equations (25), (26) and (32), as
V2 AI = 2(a4
+ aa)
to express V2A& explicitly in terms of AAj (j E Ni) in the form (24)we should determine coefficient
gij given by
gij = g'vj
gT = 2(4th row of G-' + 6th row of G-')
It is remarked that the coefficients gij only depend on the co-ordinates of the neighbouring
integration points. In case of small displacements, they are computed once for all and used at each
time step to evaluate the Laplacian by (24). However, in case of large strains, they should be
updated when geometry changes are significant. The choice of the neighbouring integration
points is simply based on the mesh topology. Generally, the integration points of the neighbouring elements are used. However, if these points are not sufficient to allow the computation of G(i.e. if G is singular or ill-conditioned),indirect neighbours (i.e. the neighbours of the neighbours)
must be considered.
It is worthwhile to point out that the polynomial approximation assumed for the plastic
multiplier function in the neighbourhood of the integration point i is used only as a means to
determine approximately the Laplacian of the plastic multiplier at the point, and does not mean
that the plastic multiplier itself is approximated in a least-squares sense.
Substitution of expression (24) into equation (22) gives
for each integration point i (i = from 1 o N,, N , is the total number of the integration points in the
finite element mesh).
To ensure quadratic convergence of the global Newton iterative procedure, a consistent
tangent matrix for the gradient-dependent von Mises elasto-plasticity has to be formulated. For
classical continuum models we can write
for each integration point i (i = 1 to Nt),where iri and iiare the stress and the strain rates at the
local point i, Dep,is the local consistent tangent matrix. It is noted that the consistent tangent
matrix at an integration point for the gradient plasticity is no longer local, the stress rate 6i at
point i depends not only on the internal state variable at the local point i, but also on those at the
neighbouring integration points. Consequently, in contrast to the form (37), the consistent
tangent constitutive relation of the point i will be given in the form
j= 1
where N, stands for the number of integration points within the elasto-plastic zone, in which
point i lies. However, it can be numerically observed that only those Depvrelating to integration
pointj, which are among the neighbouring ones of the point i, give the essential contributions to
the constitutive relation (38). Hence, it is accurate enough to assume that
Based on this fundamental consideration, let us consider the non-local consistent tangent
matrices Dep, (j = 1,N i ) at integration point z, with the use of the mixed four-node element.”
Let L(i) and N ( i ) be the set of the neighbouring integration points and the number of these
points for integration point i. Define Li(k) = L ( i ) n L ( k )and let Ni(k) be the number of points in
the intersection Li(k). The variation of equation (19),together with the consistency condition
f k = 0 and equation (13), for each point k E L(i) at time t -k At gives
j = l,jeL.i(k)
Solution of the equation set (40) gwes the rates of the effective plastic strains for each point
k E L(i)
N (i)
j= 1
Substitution of equation (41) into the rate form of the equation (19) for point i gives
In the strain driven finite element analysis the rate of the trial elastic stress +Af is entirely
determined by the strain rate. A constitutive relation between hEand the strain rate can be built
up. To do this, we first review the mixed finite element12 used in the present work briefly. The
strain of the mixed element is chosen in the form
+ h,$
and correspondingly, the stress takes the form
u = 6 + h*,UX
where a overbar denotes a constant field, the second term on the right-hand side in each of
equations (44)and (45) varies with co-ordinates. As the strain vector for the plane strain analysis
with the form E = [ex. eyy E,, cxyITis assumed, we have
where A is the area of the finite element.
The constant fields (6,E) are uncoupled with the anti-hourglass mode (a",sX).The stress m and
the strain E can be further decomposed into deviatoric and volumetric parts by
s =S
+ h,,a'
+ h,.Er
= E - = f-T'f
e =8
withm= [l 1 1 0IT.
Since the anti-hourglass mode is introduced into the element, one integration point scheme is
used for two-dimensional four-noded element while neither zero-energy modes nor locking
phenomena for shear or incompressible materials arises. The effective stress defined in equation
(8) is replaced by the mean value of it over an element as
cF2 = - / A ($sTs)dA
Correspondingly,CE redefined as the mean value of it over the element and its derivative can be
given by
= - JA
Define the ratio
The rate of the deviatoric stress can be divided into a local and a non-local part as
si = sf:+ 5;'
After some classical deductions the integration of the virtual work done by the stress rate
6 over the element i gives
The first and the second terms on the right-hand side of the equation represent the local and
non-local tangent stiffness, with the local consistent compliance matrices
C k i = A&P
- B3SiST)
CieSi= - 2jI,~,(af)~H~
C:evi= - 2P3HiaT$
and the non-local consistent compliance matrices
The coefficients in equations (62) can be given by
and the matrix P for the plane strain and the axisymmetric solid problems is given by
To analyse geometrically non-linear problems with finite strains and/or finite rotations by
utilizing the mixed element, a co-rotational formulation is adopted. For each element a local
co-ordinate system is embedded at the centre (the integration point) of the element and rotates
with it. The element displacements,strains and stresses are evaluated with reference to the local
co-ordinates at the current time.”
Consider a load step from time t to t At. The algorithm described in the proceeding
equations can be summarized as follows:
(1) Compute the element stiffness matrices according to the local and the non-local consistent
tangent matrices given by (62) and (63) for each element.
(2) Solve for the global displacement increments and then obtain the total displacements at
time t + At referred to the global co-ordinates.
(3) Determine the orientation of local axes at time t + At for each element according to the
total displacements at time t At, and the element displacementsat the last converged load
step referred to the local co-ordinate system of the element at time t. The incremental
strains at time t + At referred to the local axes at time t + At is then obtained.12
(4) Check the yield function for each integration point according to equation (23), i.e. check if
the point is in elasticity, or on the current yield surface or outside the current yield surface.
Then 6(AAi)for each integration point i is calculated non-locally according to equation (36).
Since this equation is linear in the unknowns and is diagonal dominated, Gauss-Seidel
iterations are used at this level.
( 5 ) Accumulate the effective plastic strain, its Laplacian and update the stress vector for each
integration point.
(6) Check global convergence criterion. If convergence does not occur go to (1).
We consider a rectangular bar with initial length Lo = 5 and initial width Bo, meshed by N x M
four-noded rectangular mixed elements, shown in Figure l(a). The bar fixed at the bottom and
loaded by an increasing prescribed horizontal displacement at the top, is analysed as a plane
strain problem. In order to have a pure shear problem, the vertical displacements of all the nodes
are prevented and the horizontal displacements of nodes pertaining to the same horizontal line
are imposed to be equal. It is noted that at least N = 3 elements over the width Bo of the bar are
required. Indeed, with N = 1 or 2 and only one integration point per mixed element, it would be
impossible to approximate the plastic multiplier by (25).
To demonstrate that pathological mesh dependence has been overcome by using the present
finite element method for gradient plasticity, we consider first an example. Four test cases with
different dimensions in x co-ordinate and element meshes are specified Case 1: Bo = 3-0, M = 5,
N = 3; Case 2 Bo = 1-75, M = 15, N = 5; Case 3: Bo = 1-0,M = 25, N = 5; Case 4 Bo = 0.7,
M = 35, N = 5. The values of material elastic properties used for the entire bar are E = 210000
and p = 0.0, The initial yield strength ayo= 240 except for the elements which are between
y = - 05 and y = 0.5. The initial yield strength for these elements has been weakened to the
value of cy0= 230. As the effectivestress at a material point reaches oy0,a softening modulus
h, = - loo00 is used for all elements, including the weakened elements. A value of the non-local
material parameter c = 2000,which is related to the internal length I by8
is used for the four test cases.
Figure l(a). A two-dimensional bar with length L and width B subjected to a shear stress
Numerical results show that the thickness of the localization zone and the effective plastic
strain distribution along y axis over the zone rapidly converge to a unique solution, which is
independent of the finite element mesh and is only related to the softening modulus h, and the
non-local material parameter c. This is illustrated in Figures 1(b) and 1(c). Figure 1(b) illustrates
the convergence of the shear stress-displacement (at the nodes .on the top face) curves. The
solution for 3 x 5 element mesh is softer, however, the solutions for the other three cases with
5 x 15, 5 x 25, 5 x 35 elements converge to a unique solution. Particularly, the curves for 5 x 25
and 5 x 35 element meshes almost overlap. Figure 1(c) illustrates the effective plastic strain
at the top
distributions along y axis for a prescribed horizontal displacement Au = 8-56 x
face of the bar. Because of symmetry only the distribution on the half of the bar is illustrated. It is
observed that the effective plastic strain profiles for 5 x 15, 5 x 25 and 5 x 35 element meshes
almost coincide,though the profile for 3 x 5 element mesh deviates somewhat from the converged
solution. It is remarked that the localization zones gradually grow up with the increasing
prescribed displacement in the test cases. The thicknesses of the converged localization zones for
the latter three test cases agree very well with the analytical solution w determined by8
w = 2 z J X
The softer behaviour of the coarse mesh is easily understood considering the strain distribution
(Figure l(c)). Indeed, the plastic strain at y = 0 (the integration point of the central element) is
much higher than that with the finer meshes: 0.037 instead of 0.026 approximately. From (10)and
Displacement (E-2)
Figure 1 (b). Shear stress curves with increasing displacement Au of the top face: (0)3 x 5 elements,(A) 5 x 15 elements,
( x ) 5 x 25 elements, (0)5 x 35 elements
Coordinate in Y aixs
Figure l(c). Effectiveplastic strain distributionalong Y axis of the bar subjected to a pure shear loading with a weakened
part between Y = 0.5 and Y = 0.5 for different discretizations: (0)3 x 5 elements, (A) 5 x 15 elements, ( x ) 5 x 25
elements, (0)
5 x 35 elements
Figure 2(a). Shear stress curves with increasing displacement Au of the top face for different v d u a of the non-local
material plastic parameter c: ( X ) c = 450, (A) c = 275, (0)
c = 100
with the chosen values of h, and c, the yield stress is much lower for the coarse mesh (despiteof the
attenuating effect of the non-local term), which explains the lower shear stress obtained.
To show the dependence of the thickness of localization zone on the non-local material
parameter c relating to the internal length, we consider a second example. Three test cases, with
the same dimensions, element mesh and the initial yield strength as those used for the test case 2 in
the first example, are executed. To display the performance of the present method at extremely
large strain, we use the material elastic properties E = 100o0,p = 0.0,and the softening modulus
h, = - 2500 for the entire bar. The non-local material parameters used for the three cases are:
c = 450, c = 275, c = 100, respectively. The shear stress-displacement curves are given in
Figure 2(a). Figure 2(b)illustrates the thicknesses of the localization zones and the effectiveplastic
strain distributions within the zones for the three test cases. It is observed that using larger value
of the internal length (i.e. the non-local material parameter c) results in the wider localization zone
and the lower peak value of the effective plastic strain. It can be explained by the fact that more
material points in the wider localization zone contribute their residual load-carrying capacity to
resist the further plastic deformation due to the gradient (non-local) effect.
A finite element method for the gradient plasticity, in which the yield function depends not only
on the local internal state variable, but also its spatial gradient, is proposed. The major novel
character of the present approach is that the plastic multiplier in plasticity at each integration
point is still treated as the internal state variable, determined by the yield criterion at the local
point in the non-local sense. Hence the non-local yield condition and the constitutive equation
Figure 2(b). Effective plastic strain distribution along Y axis of the bar subjected to a pure shear loadingwith a weakened
part between Y = - 0 5 and Y = 0.5 for different values of the non-local material plastic parameter c: ( x ) c = 450, (A)
c = 275, (0)c = 100
can be fulfilled in a pointwise fashion at each iteration of a load step, even though as a non-local
quantity, the spatial gradient of the local internal state variable is introduced into the yield
condition and it is approximately determined by the local plastic multipliers at neighbouring
integration points surrounding each individual integration point in a least-squares sense.
The authors are pleased to acknowledge the support of this work by the National Science
Foundation of China and the Education Committee of China. The second author would also like
to thank the FNRS (Belgium National Science Foundation) for its support during his sabbatical
leave at Dalian University of Technology.
1. Z. P.Bazant, T. Belytschko and T. P.Chang, ‘Continuum theory for strain softening’, J . Eng. Mech. ASCE, 110,
1666-1692 (1984).
2. Z. P.Bazant, ‘Imbricate continuum and its variational derivation’, J. Eng. Mech. ASCE, 110,1693-1712 (1984).
3. Z. P. Bazant and L. Cedolin, Stability of Struczures: Elastic, Inelastic, Fracture and Damage Theories, Oxford
University Press, New York, 1991.
4. E. C. Aifantis, ‘On the microstructural origin of certain inelastic models’, Trans. ASME. J . Eng. Mat. Tech. 106,
326-330 (1986).
5. E. C. Aifantis, ‘The physics of plastic deformation’, Int. J. Plasticity, 3,211-247 (1988).
6. D.Lasry and T. Belytschko, ‘Localization limiters in transient problems’, Int. J. Solids Struct., 24, 581-597 (1988).
7. H.B. Muhlhaus and E. C. Aifantis, ‘A variational principle for gradient plasticity’, I#. J. Solids S m r . , 29,845-857
R. de Borst and H. B. Muhlhaus, ‘Gradient-dependent plasticity: formulation and algorithmic aspects’, Int. j . numer.
methods eng., 35, 521-539 (1992).
9. L. J. Sluys, Borst and H. B. Muhlhaus, ‘Wave propagation, localization and dispersion in a gradient-dependent
medium’, Int. J . Solids Srruct., 30, 1153-1171 (1993).
10. H. B. Muhlhaus, ‘Application of Cosserat theory in numerical solutions of limit load problems’, Ing.-Archiv, 59,
124-137 (1989).
11. A. Needleman, ‘Material rate dependence and mesh sensitivity in localization problems’, Comp. Methods Appl. Mech.
Eng., 67, 69-86 (1988).
12. Ph. Jetteur and S. Cescotto, ‘A mixed finite element for the analysis of large inelastic strains’, 1nt.j. numer. methods
eng., 31,229-239 (1991).
13. J. C. S i o and R. L. Taylor, ‘Consistent tangent operations for rate independent plasticity’, Comp. Methods Appl.
Mech. Eng., 48, 101-118 (1985).
14. 0.C . Zienkiewicz and R. L. Taylor, The Finite Element Method, Vol. 1, 4th edn, Mcgraw-Hill, New York, 1989.
15. 0. C. Zienkiewicz and R. L. Taylor, The Finite Element Method, Vol. 2,4th edn, Mcgraw-Hill, New York, 1991.
Без категории
Размер файла
651 Кб
Пожаловаться на содержимое документа