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



код для вставкиСкачать
Int. J. Numer. Meth. Engng. 44, 1823–1841 (1999)
Koiter Institute Delft=Department of Civil Engineering; Delft University of Technology; P.O. Box 5048;
2600 GA Delft; The Netherlands
The nite element reliability method is applied to evaluate the eects of stochastic imperfections on the
localization behaviour of elasto-viscoplastic softening solids. Material properties as the Young’s modulus, the
softening modulus and the initial yield stress are considered to be random elds. Failure criteria are dened
for the peak load and the dissipated energy. Likely congurations of imperfections which lead to failure
are obtained and the nature and relative importance of the corresponding localization patterns are analyzed.
Copyright ? 1999 John Wiley & Sons, Ltd.
KEY WORDS: localization; viscoplasticity; stochastic imperfections; nite element reliability method
Mode-II localization in softening solids can be described by means of inelastic material models.
In particular, the inclusion of rate dependence in the constitutive description by means of elastoviscoplastic models is of great interest, since these models do not suer from the mathematical
ill-conditioning, which is observed when classical rate-independent elasto-plastic models are combined with constitutive laws which exhibit a softening branch. This ill-conditioning manifests itself
as the emergence of innite inelastic strains within a zone of measure zero while no energy is dissipated. Finite element approximations of such a process are therefore meaningless, since they try to
accommodate the inelastic strains in the smallest zone possible, e.g. a row of elements. The energy
dissipation is consequently pathologically dependent on the discretization size. Elasto-viscoplastic
models introduce strain-rate eects in the evolution of the stress–strain path by means of a uiditylike parameter. The mathematical description of the problem then remains well-posed and nite
element simulations converge to a unique, nite energy dissipation upon mesh renement.1 Different viscoplastic formulations have been proposed by Perzyna,2 Duvaut and Lions3 and Wang
et al.4 However, a proper mathematical description of the problem is not sucient to understand
the behaviour of a solid in which localization takes place. The presence of material defects can
be of great importance for the response of a viscoplastic solid.5 Indeed, an imperfection is often
Correspondence to: R. de Borst, Koiter Institute Delft=Department of Civil Engineering, Delft University of Technology,
P.O. Box 5048, 2600 GA Delft, The Netherlands. E-mail:
CCC 0029–5981/99/121823–19$17.50
Copyright ? 1999 John Wiley & Sons, Ltd.
Received 24 November 1997
Revised 2 July 1998
necessary to trigger localization in numerical analysis. However, such practice is by no means a
computational ‘trick’. Localized deformations in laboratory specimens are inuenced by imperfections to a considerable extent, which can be found in the material properties and the boundary
conditions.6 When dierent imperfection patterns lead to dierent equilibrium paths, a meaningful
simulation can only be carried out if a stochastic description of the imperfections is considered.
Research in this subject, making use of Monte Carlo simulations, has been presented by Carmeliet
and de Borst7 for concrete fracture. Unfortunately, Monte Carlo simulations are expensive from a
computational point of view and, while they give accurate results, they do not provide any insight
in the nature of the most probable localization patterns.
In this study a dierent strategy is used. Material parameters of a viscoplastic solid like the initial
yield stress and the elastic and softening moduli are considered as random elds and coupled to
the peak load or the energy dissipation through a non-linear nite element algorithm. The nite
element reliability method8 is then used to obtain realizations of the material random elds that
furnish a local maximum of the probability density of failure. This is achieved by means of a
gradient-based optimization algorithm. When the local maxima have been found, the marginal
probability of dierent modes and the global probability of failure can be approximated.
The basic idea of the nite element reliability method is to approximate the statistics of the mechanical response through a transformation of this response into the space of random material
properties. Since these properties are known as a random eld, they must be discretized into a
vector of random variables. The realization of this vector which represents a local maximum of its
probability density and which violates the safety criteria, is found through an optimization algorithm, that requires the evaluation of the response gradient with respect to the material properties.
2.1. Discretization of the material random elds
Let V be a homogeneous random eld describing a material property in the physical domain .
Homogeneity means in this context that the statistical parameters of the eld are invariant to translations on the domain. Therefore, V is completely dened by a probability distribution function
and a correlation function (x1 ; x2 ) which gives the correlation coecient between two dierent
points x1 and x2 and is exclusively dependent on the distance |x1 − x2 |. Common expressions of
in structural mechanics are
|x1 − x2 |
max 0; 1 −
|x1 − x2 |
(x1 ; x2 ) = exp −
|x − x |2
 exp − 1 2 2
which are known as triangular, exponential and Gaussian correlation respectively. The parameter ‘
accounts for the scaling of the correlation decay as the distance between the two points increases
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 1823–1841 (1999)
and is dened as correlation length. The random eld V must be discretized into a vector of
random variables V, which is characterized by a joint probability density function fV . Several
methods exist for this purpose.9 Without loss of generality, the mid-point method has been adopted
in this study, in which the domain is discretized into a nite number of subdomains i . The random
eld within i is considered to be constant and equal to a random variable Vi , which is dened
Vi = V(xic )
where xic is the central point of i . The marginal statistics of Vi are obviously those of V and the
correlation coecient between Vi and Vj is given by (xic ; xjc ). In order to facilitate the computer
implementation, it is usual to select the subdomains i coinciding with patches of the nite element
discretization of the mechanical problem, which exhibits to a certain extent similarities with the
usual practice of placing a defect in a patch of nite elements to trigger localization in deterministic
The vector V of random variables can be converted into a vector of uncorrelated, standard
normally distributed variables Y by means of a multivariate transformation.10 In the case that V
is already normally distributed, this transformation reduces to an ane transformation. Converting
V into Y is not strictly necessary; however, the uncorrelated standard normal space is easier to
operate with, because the corresponding probability density function is monotonically decreasing
as the distance to the origin increases.
2.2. Mechanical transformation: from random properties to random response
The relation between the (discretized) random eld of material properties and the response in
quasi-static conditions, neglecting body forces, is governed by the boundary value problem
∇·b = 0
in u = u(t)
in @
b · n = b(t)
in @
at each instant t, where b is the stress tensor, u is the displacement eld, @
1 ∪ @
2 = @
, n is
are the prescribed boundary displacements and
the outward normal vector to @
and u(t)
and b(t)
loading. The relation between the stress tensor, the displacements eld and the material properties
is formalized by means of a Duvaut–Lions elasto-viscoplastic material model. For such material,
the constitutive equations in rate form are written as
ḃ = De U̇ − (b − bp )
˙ = − ( − p )
where U is the strain tensor, is the hardening (softening) parameter, De is the linear elastic
tangent operator, is a uidity-like parameter and the subindex p signies that these quantities
belong to a rate-independent elasto-plastic problem. Equations (4) are valid when the yield function
is non-negative
f(b; )¿0
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 1823–1841 (1999)
otherwise the material is in the elastic regime,
ḃ = De U̇
The underlying rate-independent elasto-plastic problem (denoted by p above) follows the usual
@f (7)
ḃp = De U̇ − ˙p
@bp p
when the yield function vanishes
f(bp ; p ) = 0
ḟ(bp ; p ) = 0
and Prager’s consistency condition
is complied with. The dependence of these equations on the material parameters can be found in
the elastic operator De and the yield function f.
When the stress does not exhibit a monotonic increasing evolution law, i.e. when softening takes
place, it is possible that the boundary value problem (3) has no real solution at each t for arbitrary
that furnishes a solution of (3) for all
prescribed boundary forces b(t).
Each pair of elds (u; b)
t is dened as equilibrium path of the solid . If we restrict ourselves to the case of proportional
loading, that is
= (t)bd
where is a scalar valued multiplier and bd is a xed, design load, we can represent the equilibrium
path corresponding to bd as (u; ). Functionals can be dened of the equilibrium path to represent
characteristics of practical interest. For instance,
c = max t
represents the maximal load multiplier, which is typically attained when ˙ = 0. The critical load
is then dened as
bc = c bd
Since c is coupled to the random material properties, it can be considered as a random variable,
so that c is a realization of c . A random measure of the capacity of the structure to resist the
design load bd is then given by
Z(c ) = c − 1
According to this equation, positive realizations of Z represent a safe structure whereas negative
realizations reect failure. In a similar fashion, energy functionals can be dened to study the
degree of brittleness.
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 1823–1841 (1999)
Figure 1. Probability content of the set z(c )¡0
2.3. Statistics of the response
The properties of function (12) enable us to dene the probability of failure in a straightforward
manner as
f c (c ) dc
Pf =
z(c )¡0
This is schematically depicted in Figure 1. However, the probability distribution function of V is
available rather than that of c . Making use of the mechanical transformation dened in Section 2.2
and the multivariate transformation introduced in Section 2.1, the state of the structure can be
formally recast as
G(Y) = Z(c (V(Y)))
Since Y follows an uncorrelated standard normal distribution, the probability of failure can be
expressed as
(y) dy
Pf =
where is the uncorrelated standard normal probability density function. Evaluation of the integral
in (15) can be dicult for arbitrary shapes of the surface g(y) = 0. The usual strategy in the
reliability method is then to expand this surface in a Taylor series around a point y∗ ,
g(y∗ ) + ∇g(y∗ )(y − y∗ ) + 12 (y − y∗ )T ∇(2) g(y∗ )(y − y∗ ) + · · · = 0
in order to approximate integral (15). In (16) the superscript T implies a transpose. If development
(16) is truncated after the linear terms, the approximation is called First-Order Reliability Method
(FORM) whereas if second-order terms are included, the approximation is known as Second-Order
Reliability Method (SORM). In this study only rst-order approximations are considered. Due to
the spherical symmetric properties of , the major contribution to the probability integral (15)
comes from the closest points to the origin. It is then meaningful to choose y∗ as the most
central point of g(y) = 0. The distance of this point to the origin is called reliability index and
usually denoted as . In accordance with this notation, the most central point is called -point and
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 1823–1841 (1999)
Figure 2. Approximation of the probability of failure by means of FORM
represented by y . It can readily be proved that follows a standard normal distribution. When
FORM is used, the probability of failure is consequently approximated by the expression
Pf ≈ (−)
where is the one-dimensional standard normal cumulative distribution function. This approximation is illustrated in Figure 2.
Finding the -point is the main task in this approximation procedure. This is formulated as a
constrained optimization problem
minimize kyk
subject to g(y) = 0
Several algorithms have been proposed to solve this problem.11 They usually consist of a recursive
sequence of points in the form
y(k+1) = F(y(k) ; z(y(k) ); ∇g(y(k) ))
where k accounts for the iteration number and F represents a function and=or a procedure (e.g.
a line search). Sequence (19) is iterated until convergence is attained upon a criterion based on
incremental norms. In this study, the HL-RF method12; 13 has been adopted, for which (19) is
written as
∇g(k) y(k) − g(k) ∇g(k)
y(k+1) =
k∇g k
According to this algorithm, the -point is obtained as a scaling of the gradient of the constraint g.
An accurate computation of this gradient is therefore of capital importance for the stable and
meaningful convergence of (20). By means of the chain rule, ∇g is expressed as
∇g = ∇c ∇v
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 1823–1841 (1999)
when denition (12) of Z is used. The computation of c and its gradient by means of the nite
element method will be described in Section 3.
The constrained optimization problem (18) can exhibit several local minima. This problem
becomes especially signicant when softening behaviour is considered and localization takes place.
Each local minimum of (18) can lead to a dierent localization pattern. Unfortunately, when (19)
has converged to a local minimum there is no way to force convergence to another minimum.
Dierent solutions are obtained if (19) is started from dierent points. In simple problems, the
most adequate starting points can be chosen so as to keep or break symmetries. Algorithms for
automatic generation of all likely failure mechanisms are only available for frame structures.14 An
extension of these algorithms to continua is a topic of current research. When the most signicant
-points have been found, the global probability of failure can be bounded.15
The mechanical transformation dened in Section 2.2 must be computed by discretizing u in nite
ai Ni
where a is dened as the nodal displacement vector and N represents the shape functions. The
discretized version of (3) with the load parameter as dened in (10) is then integrated in time
with an Euler-forward scheme according to the expressions
a = t a + a
= t + (23)
where the notation t = (t) has been used for the dependence of any object on time. Instead, the
material equations (4)–(9) are solved at each integration point with an Euler-backward scheme.
The gradient of the mechanical transformation is obtained from a careful dierentiation of the
global as well as the local algorithms. A study on the existence, behaviour and computation of
these derivatives in softening solids can be found in Reference 16.
3.1. Computation of the equilibrium path
3.1.1. Global equilibrium. Elaborating a weak form of equation (3) for discretization (2) of
the material properties and discretization (22) of the displacements eld, the non-linear system of
algebraic equations
fint (a; v) = t fd
is obtained, where the nodal internal force vector is dened as
fint =
BTt b d
B is the strain-nodal displacement matrix, fd accounts for the nodal design forces and b represents
the stress tensor, which is related to the nodal displacement vector a and the material properties v
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 1823–1841 (1999)
through the constitutive equations. Using equations (23) the evolution of the equilibrium path can
be written incrementally as
BT b d
= fd
assuming that equilibrium is complied with at time t. A load control procedure, in which is supplied by the analyst, cannot be used because this increment can be incompatible with the
equilibrium path (see Section 2.2). Instead, a selective arc-length control procedure is used.17 In
the selective arc-length method, both the nodal incremental displacements and external forces are
parametrized by the length of a projection of the nodal incremental displacement vector. This is
formally done by introducing a constraint equation
aT a = l2
with a projection matrix to select the degrees-of-freedom which contribute to the increment of
arc-length l, which, in a context of quasi-static loading, can be considered congruent with the
time increment t. However, such a restriction is not necessary and in general it can be considered
that l is any function of t. Combining (26) and (27) we obtain
BT b d
= fd
aT a = l2
which system is solved with the Newton–Raphson method
 −
T (k)
 r(a; )(k)
where K is the consistent tangent stiness matrix, which will be discussed in the next section, and
the residual r is dened as
r(a; ) =  
aT a − l2
From this procedure and equations (23) a time-discretized equilibrium path is obtained. However,
a point such that ˙ = 0, which is associated with c , does not generally exist in the context of
spatially discretized systems with softening behaviour. In this study it has been considered that
the peak load multiplier c can be interpolated between t−t , t and t+t when t is larger than
both t−t and t+t .
3.1.2. Integration of the constitutive relation. When the yield criterion (5) is satised at a
generic integration point of the nite element discretization, equations (7) – (9) for the underlying
rate-independent elasto-plastic problem must be integrated in time to obtain a time discretized
increment of stress b. This is done with an Euler-backward return-mapping algorithm.18
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 1823–1841 (1999)
Without loss of generality, a Von Mises material in plane-stress conditions with a linear hardening
(softening) rule is considered,
(I − A)t bp + De U = Abp
f(t bp + bp ; t p + p ) = 0
A = I+
f(b; ) =
De P
+ p )
t p
3J2 − ()
= 0 + h
J2 = 12 bT Pb
In these equations, h is the hardening (softening) modulus, is the yield stress and
2=3 −1=3 0
2=3 0 
P =  −1=3
Equations (31) are solved with any suitable iterative procedure for a given increment of strain
U. The consistent tangential relation between bp and U is given by the matrix18
Dcp = H −
H(@f=@b)(@f=@b)T H
h + (@f=@b)T H(@f=@b)
@2 f
I + p De 2
Equations (4) are integrated with an Euler-backward algorithm as well, resulting in the explicit
b =
(De U + t=bp − t=(t b − t bp ))
1 + t=
(t=p − t=(t − t p ))
1 + t=
The consistent tangential relation between b and U has been derived by Wang et al.,3
Dc =
De + t=Dcp
1 + t=
This consistent tangent operator is assembled in the stiness matrix in (29) in order to ensure a
quadratic rate of convergence of the Newton–Raphson procedure. Moreover, the accurate evaluation
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 1823–1841 (1999)
of the response gradient requires the use of the consistent tangent operator (see Section 3.2). A
closed form of equations (36) is given by Simo et al.,19 but will not be used here because it does
not provide a consistent tangent operator.
3.2. Computation of the gradient of the equilibrium path
A general methodology has been presented in Gutierrez and de Borst16 to evaluate the gradient of
the equilibrium path in softening solids. The particular case for a Duvaut–Lions elasto-viscoplastic
solid will be given here. Adopting the notation
@i =
for the partial derivative with respect to a single coecient of the discretization of the material
property, dierentiation of equation (28) with respect to vi leads to
BT @i b d
− @i f̂ = 0
2a @i a = 0
The derivative of the stress increment b with respect to vi can be decomposed into two contributions
@i b = @i b|U +
@i U
where the term @b=@U is the consistent tangent operator (37). Substituting (40) into (39), the
B @i b|U d
+ BT
@i U d
− @i f̂ = 0
2a @i a = 0
results, which, using @i U = BT @i a, can be recast as
 Z
− B @i b|U d
@i a
@i 0
2aT 0
This is a system of linear simultaneous equations with a matrix identical to that in equation (29)
at a converged state.
The computation of the right-hand side of (42) is explained next. Dierentiation of equations
(31) with respect to vi yields
@i At+t bp + A@i bp = (I − A)@i t bp + @i De U + De @i U
3t+t bpT P(@ti bp + @i bp )
t p + p )
= @i (
2 (
t+t p )
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 1823–1841 (1999)
where the terms @i t bp and @i t p are known from the previous time increment. Developing terms in
(43) for constant U, a linear system of four simultaneous equations is obtained which is solved
for @i bp |U and @i p |U . Dierentiation of (36) then gives
@i b =
(@i De U + De @i U + t=@i bp − t=(@i t b − @i t bp ))
1 + t=
@i =
(t=@i p −t=(@i t − @i t p ))
1 + t=
For constant U, equations (44) then yield @i b|U and @i |U which furnish the right-hand
side of (41). Solving (41) @i a and @i are obtained. The derivatives of the strain increment
are obtained from @i U = B@i a and those of bp , p , b and are solved from (43) or
(44). Finally, the derivatives at t + t are updated as in (23).
4.1. Shear layer problem
To investigate the most probable width of the localization zone appearing under pure shear
conditions, the shear layer of Figure 3 has been analyzed. The shear layer is an innitely long strip
in the x-direction with a stress jump at the boundaries. Material imperfections depend exclusively
on the y-coordinate. A scheme of the nite element model and the discretization of the random eld
in the layer, which has a height A = 100 mm and a thickness of 1 mm, is depicted in Figure 3.
The problem is basically one-dimensional and could have been analyzed with line elements as
well. Use of two-dimensional elements requires the introduction of linear constraint equations in
the horizontal displacements. The displacements in the vertical direction are prevented, the bottom
is xed and the loading is controlled by a selective arc-length procedure such that the horizontal
velocity u̇x = 1 mm=s at the top. The Poisson’s ratio is taken to be zero and the Young’s modulus
E, the hardening (softening) modulus h and the initial yield stress 0 are normally distributed
random elds according to the parameters
E ∝ N (100 000 N=mm2 ; 10 000 N=mm2 )
H ∝ N (−4 000 N=mm2 ; 400 N=mm2 )
0 ∝ N (1 000 N=mm2 ; 100 N=mm2 )
in which N denotes the normal distribution with the mean value and the standard deviation as
parameters, and an exponential correlation function (cf. equation (1)). Using a Gaussian distribution
for these random elds implies that a very small probability is accepted that E and 0 attain
a negative value. This is not meaningful from a physical point of view. Since the probability
content of the sets of negative realizations is negligible we shall retain the Gaussian distribution
for simplicity in the multivariate transformation (see Section 2.1), which reduces to an ane
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 1823–1841 (1999)
Figure 3. Innitely long shear layer: static scheme, applied loading and nite element mesh
transformation in such case. Representing the vector V as
V = (E1 ; : : : ; En ) t (H1 ; : : : ; Hn ) t ( 01 ; : : : ; 0n )
where n is the number of subdomains dened as in Section 2.1, the multivariate transformation to
a vector Y of uncorrelated standard normal variables becomes
Y = F(V − VM )
where VM is the expectation of V and F is any matrix such that
FT F = (Cov[V])−1
The matrix F is usually obtained by inverting the Cholesky factor or the square root of Cov[V].
Calculations have been carried out with three dierent values for the correlation length, namely
‘ = 10, 20 and 40 mm. A generic limit state is dened for a peak load of 1000 N. Two dierent
values have been used for the uidity parameter, = 0·015 and 0·030 s. The layer and the random
eld have been discretized into 100 six-noded triangular nite elements with a three-point Gauss–
Legendre integration quadrature and 50 random subdomains respectively. Algorithm (20) has been
used starting from y(0) = 0. A -point has been obtained after about 10 –15 iterations, depending on
the choice for ‘ and . Due to the symmetry of the problem in the y-direction and the symmetry
of the initial realization of the imperfections y(0) , the obtained collapse modes are symmetric as
well and exhibit the same properties, i.e., they consist of a shear band with a width that is basically
governed by the uidity parameter . The correlation length has been found to inuence the width
of the shear band only at the onset of localization. When the shear band has developed fully, its
width is exclusively controlled by the uidity parameter.
The development of the equivalent plastic strain until ≈ 70 per cent residual loading for
= 0·015 and = 0·030 s is represented in Figures 4 and 5 for dierent values of the correlation
length ‘, respectively. Figure 6 shows the total displacement eld corresponding to the collapse
mechanisms obtained with = 0·015 and = 0·030 s at ≈ 70 per cent residual loading, with the
correlation length ‘ = 20 mm.
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 1823–1841 (1999)
Figure 4. Evolution of the equivalent plastic strain at the -point until ≈ 70 per cent residual load for dierent values
of the correlation length ‘, with = 0·015 s
4.2. Biaxial tension test
The localization behaviour has also been analyzed for the biaxial tension test of Figure 7. The
specimen has a size of 64 × 64 mm2 and a thickness of 5 mm. The domain is discretized into
a structured mesh of 16 × 16 eight-noded plane-stress nite elements with a four-point Gauss–
Legendre integration quadrature. While coaxiality of the load has been enforced, linear constraint
equations are used at the top and the bottom to ensure that they can rotate freely. The loading is
controlled by a selective arc-length procedure such that the distance between both ends increases
with a velocity of 1 mm=s.
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 1823–1841 (1999)
Figure 5. Evolution of the equivalent plastic strain at the -point until ≈ 70 per cent residual load for dierent values
of the correlation length ‘, with = 0·030 s
The Young’s modulus, the hardening modulus and the initial yield stress are normally distributed
according to
E ∝ N (20 000 N=mm2 ; 2000 N=mm2 )
H ∝ N (−500 N=mm2 ; 50 N=mm2 )
0 ∝ N (100 N=mm2 ; 10 N=mm2 )
and an exponential correlation function. Each random eld is discretized into 256 random variables
corresponding with the nite element mesh. The Poisson’s ratio has been taken = 0·2. The design
load is 28 000 N.
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 1823–1841 (1999)
Figure 6. Total displacement pattern at ≈ 70 per cent residual load for ‘ = 20 mm. Left: = 0·015 s. Right: = 0·030 s
Figure 7. Biaxial tension test: static scheme and applied loading
4.2.1. Analysis of probable localization patterning: inuence of the uidity parameter. A common experimental observation in several materials is that biaxial tests on rectangular specimens
usually fail in an asymmetric shear band mode. However, sometimes two cross-diagonal, conjugate
shear bands are observed. When numerical simulations are carried out the ultimate failure mode
depends on the imperfection pattern which is used to trigger localization. In this sense, symmetrically and asymmetrically placed imperfections lead to symmetric and asymmetric localization
patterns, respectively. To investigate which are the most probable localization patterns when a
viscoplastic material model is used, the presented algorithms have been run for the uidity parameters = 0·0015 and = 0·0075 s, with a correlation length ‘ = 10 mm. Algorithm (20) has been
started from dierent initial imperfection congurations y(0) , which basically consist of a reduction
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 1823–1841 (1999)
Figure 8. Initial, symmetric and asymmetric imperfection patterns for the biaxial tension test
Figure 9. Imperfection patterns at -point for symmetric and asymmetric modes
of the initial yield stress with an amount of 10 per cent. The chosen initial imperfection patterns
are depicted in Figure 8 and the obtained -points for the yield stress after running algorithm (20)
are shown in Figure 9. For the small uidity ( = 0·0015 s) it is found that the reliability index of
two cross-diagonal shear bands (Figure 10) is = 4·0646 and the reliability index of a single diagonal shear band (Figure 11) is = 3·9399, which correspond to marginal probabilities of failure
Pf = 2·407 × 10−5 (symmetric mode) and Pf = 8·165 × 10−5 (union of both asymmetric modes),
respectively. According to this result, an asymmetric localization pattern is almost four times as
probable as a symmetric pattern. This is in agreement with the experimental observations. When
the uidity is increased to = 0·0075 s, convergence is found to be a single, symmetric -point,
starting algorithm (20) from either the symmetric or the asymmetric conguration. The reliability
index is = 4·5935 which represents a probability Pf = 2·186 × 10−6 and the corresponding localization pattern is a necking mode (Figure 12). This, again, agrees with the experimental results,
according to which ductile materials subject to tension exhibit a symmetric, necking failure mode.
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 1823–1841 (1999)
Figure 10. Localization pattern at symmetric -point for ≈ 70 per cent residual load and = 0·0015 s. Left: total displacements (scaling factor 25). Right: equivalent plastic strain Figure 11. Localization pattern at asymmetric -point for ≈ 70 per cent residual load and = 0·0015 s. Left: total displacements (scaling factor 25). Right: equivalent plastic strain 4.2.2. Analysis of probable localization patterning: inuence of the correlation length. In
Section 4.1 it has been shown that the correlation length has little inuence in the width of
the shear bands. To study the consequences of a variation of the correlation length on the shear
band patterning, computations have been carried out for the uidity parameter = 0·0015 s and the
correlation length ‘ = 100 mm (originally: ‘ = 10 mm). Starting algorithm (20) from both symmetric and asymmetric congurations, the reliability indices = 2·0210 and = 2·2664 are obtained
that correspond to probabilities of failure Pf = 2·167 × 10−2 and Pf = 2·346 × 10−2 (for the union
of two homologous asymmetric modes). Surprisingly, these probabilities are in the same order of
magnitude while the width of the shear bands remains unaltered with respect to the case ‘ = 10 mm.
This increased marginal probability of a symmetric mode compared to an asymmetric mode can
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 1823–1841 (1999)
Figure 12. Localization pattern at -point for ≈ 70 per cent residual load and = 0·0075 s. Left: total displacements (scaling factor 25). Right: equivalent plastic strain be explained by the smoothness of the imperfections eld corresponding to a large correlation
length. Asymmetric modes require sharper defects than symmetric modes. The probability of accommodating a sharp defect is lower than that of accommodating a smooth defect for a large
correlation length. This property has also relevant eects for the relation between the uidity and
the reliability index. For ‘ = 100 mm, the marginal probability of an asymmetric mode has also
been studied for the uidities = 0·00075 s and = 0·0030 s. The obtained reliability indices are
= 2·8419 and = 1·9453 which correspond to marginal probabilities of failure Pf = 2·247 × 10−3
and 2·591 × 10−2 , respectively. A larger uidity thus leads to a smoother imperfection pattern at
the -point, which for relatively large values of the correlation length ‘ corresponds to a larger
probability of failure.
The nite element reliability method has been used in combination with a Duvaut–Lions elastoviscoplastic model to evaluate the eect of stochastic material imperfections on the localization
behaviour. It has been found that the width of the localization zone is governed by the uidity parameter, whereas the correlation length is of importance only at the onset of localization.
Numerical simulations of a biaxial tension test have shown that asymmetric localization modes
are more likely to occur than symmetric modes when the correlation length is small, whereas
the opposite eect is found when the correlation length is large. Also, ductile materials exhibit
a larger probability of failure than brittle materials when a large correlation length is considered.
The observed interaction between the uidity parameter, the correlation length and the reliability
index provides a valuable insight in the mechanical behaviour of viscoplastic solids.
1. R. de Borst, L. J. Sluys, H.-B. Muhlhaus and J. Pamin, ‘Fundamental issues in nite element analysis of localization
of deformation’, Engng. Comput., 10(2), 99 –121 (1993).
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 1823–1841 (1999)
2. P. Perzyna, ‘Fundamental problems in viscoplasticity’, Recent Advances in Applied Mechanics, Vol. 9, Academic
Press, New York, 1966, pp. 243 –377.
3. G. Duvaut and J. L. Lions, Les inequations en Mecanique et en Physique, Dunod, Paris, 1972.
4. W. M. Wang, L. J. Sluys and R. de Borst, ‘Viscoplasticity for instabilities due to strain softening and strain-rate
softening’, Int. J. Numer. Meth. Engng., 40, 3839 –3864 (1997).
5. W. M. Wang, L. J. Sluys and R. de Borst, ‘Interaction between material length scale and imperfection size for
localisation phenomena in viscoplastic media’, Eur. J. Mech. A=Solids, 15(3), 447– 464 (1996).
6. J. Desrues, La localisation de la deformation dans les materiaux granulaires, Dissertation, USMG-INPG Grenoble,
7. J. Carmeliet and R. de Borst, ‘Stochastic approaches for damage evolution in standard and non-standard continua’,
Int. J. Solids Struct., 32, 1149 –1160 (1995).
8. A. der Kiureghian and J. B. Ke, ‘The stochastic nite element method in structural reliability’, Probab. Engng. Mech.,
3, 83 –91 (1988).
9. C. C. Li and A. der Kiureghian, ‘Optimal discretization of random elds’, J. Eng. Mech. ASCE, 119(6),
1136 –1154 (1993).
10. M. Rosenblatt, ‘Remarks on a multivariate transformation’, Ann. Math. Statist., 23(3), 470 – 472 (1952).
11. P. L. Liu and A. der Kiureghian, ‘Optimization algorithms for structural reliability’, Struct. Safety, 9, 161–177 (1991).
12. A. M. Hasofer and N. C. Lind, ‘Exact and invariant second-moment code format’, J. Eng. Mech. ASCE, 110(1),
111–121 (1974).
13. R. Rackwitz and B. Fiessler, ‘Structural reliability under combined load sequences’, Comput. Struct., 9, 489 – 494
14. Y. Murotsu, S. Shao and S. T. Quek, ‘Some studies on automatic generation of structural failure modes’, in
A. der Kiureghian and P. Thoft-Christensen (eds.), Reliability and Optimization of Structural Systems ’90; Proc.
3rd IFIP WG 7.5 Confer., Springer, Berlin, 1991.
15. O. Ditlevsen and H. O. Madsen, Structural Reliability Methods, Wiley, Chichester, 1996.
16. M. A. Gutierrez and R. de Borst, ‘Studies in material parameter sensitivity of softening solids’, Comp. Meth. Appl.
Mech. Engng., 162, 337–350 (1998).
17. R. de Borst, ‘Computation of post-bifurcation and post-failure behaviour of strain-softening solids’, Comput. Struct.,
25, 211–224 (1987).
18. R. de Borst and P. H. Feenstra, ‘Studies in anisotropic plasticity with reference to the Hill criterion’, Int. J. Numer.
Meth. Engng., 29, 315 –336 (1990).
19. J. C. Simo, J. G. Kennedy and S. Govindjee, ‘Non-smooth multisurface plasticity and viscoplasticity. Loading=unloading
conditions and numerical algorithms’, Int. J. Numer. Meth. Engng., 26, 2161–2186 (1988).
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 1823–1841 (1999)
Без категории
Размер файла
226 Кб
Пожаловаться на содержимое документа