INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING Int. J. Numer. Meth. Engng. 44, 1823–1841 (1999) NUMERICAL ANALYSIS OF LOCALIZATION USING A VISCOPLASTIC REGULARIZATION: INFLUENCE OF STOCHASTIC MATERIAL DEFECTS M. A. GUTIERREZ AND R. DE BORST ∗ Koiter Institute Delft=Department of Civil Engineering; Delft University of Technology; P.O. Box 5048; 2600 GA Delft; The Netherlands SUMMARY 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 1. INTRODUCTION 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: R.deBorst@ct.Tudelft.nl CCC 0029–5981/99/121823–19$17.50 Copyright ? 1999 John Wiley & Sons, Ltd. Received 24 November 1997 Revised 2 July 1998 1824 M. A. GUTIERREZ AND R. DE BORST 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. 2. THE FINITE ELEMENT RELIABILITY METHOD FOR VISCOPLASTIC SOLIDS 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 | (1) (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) INFLUENCE OF STOCHASTIC MATERIAL DEFECTS 1825 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 as Vi = V(xic ) (2) 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 analysis. 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 @ 1 b · n = b(t) in @ 2 (3) 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 1 ḃ = De U̇ − (b − bp ) 1 ˙ = − ( − p ) (4) 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. (5) Int. J. Numer. Meth. Engng. 44, 1823–1841 (1999) 1826 M. A. GUTIERREZ AND R. DE BORST otherwise the material is in the elastic regime, ḃ = De U̇ (6) The underlying rate-independent elasto-plastic problem (denoted by p above) follows the usual relation ! @f (7) ḃp = De U̇ − ˙p @bp p when the yield function vanishes f(bp ; p ) = 0 (8) ḟ(bp ; p ) = 0 (9) 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 b(t) (10) 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 (11) 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 (12) 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) INFLUENCE OF STOCHASTIC MATERIAL DEFECTS 1827 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 Z f c (c ) dc (13) 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))) (14) Since Y follows an uncorrelated standard normal distribution, the probability of failure can be expressed as Z (y) dy (15) Pf = g(y)¡0 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 (16) 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) 1828 M. A. GUTIERREZ AND R. DE BORST 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 ≈ (−) (17) 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 (18) 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) )) (19) 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 T 1 (20) ∇g(k) y(k) − g(k) ∇g(k) y(k+1) = (k) 2 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. (21) Int. J. Numer. Meth. Engng. 44, 1823–1841 (1999) INFLUENCE OF STOCHASTIC MATERIAL DEFECTS 1829 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 3. COMPUTATION OF THE MECHANICAL TRANSFORMATION The mechanical transformation dened in Section 2.2 must be computed by discretizing u in nite elements, P ai Ni (22) u= i 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 t+t a = t a + a t+t = 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 t fint (a; v) = t fd (24) is obtained, where the nodal internal force vector is dened as Z t fint = BTt b d (25) 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) 1830 M. A. GUTIERREZ AND R. DE BORST through the constitutive equations. Using equations (23) the evolution of the equilibrium path can be written incrementally as Z BT b d = fd (26) 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 (27) 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 Z BT b d = fd (28) aT a = l2 which system is solved with the Newton–Raphson method a (k+1) = a (k) − K 2a T (k) −fd r(a; )(k) (29) 0 where K is the consistent tangent stiness matrix, which will be discussed in the next section, and the residual r is dened as Z T B b d − f d (30) 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) INFLUENCE OF STOCHASTIC MATERIAL DEFECTS 1831 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 (31) where A = I+ f(b; ) = p 3p De P + p ) ( t p 3J2 − () (32) () = 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 0 0 2 (33) 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 − where H= H(@f=@b)(@f=@b)T H h + (@f=@b)T H(@f=@b) @2 f I + p De 2 @b (34) −1 De (35) Equations (4) are integrated with an Euler-backward algorithm as well, resulting in the explicit formulae b = 1 (De U + t=bp − t=(t b − t bp )) 1 + t= 1 = (t=p − t=(t − t p )) 1 + t= (36) The consistent tangential relation between b and U has been derived by Wang et al.,3 Dc = De + t=Dcp 1 + t= (37) 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) 1832 M. A. GUTIERREZ AND R. DE BORST 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 = @ @vi (38) 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 Z BT @i b d − @i f̂ = 0 (39) T 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 + @b @i U @U (40) where the term @b=@U is the consistent tangent operator (37). Substituting (40) into (39), the system Z Z @b T B @i b|U d + BT @i U d − @i f̂ = 0 @U (41) T 2a @i a = 0 results, which, using @i U = BT @i a, can be recast as Z # " T K −f̂ − B @i b|U d @i a = (k) @i 0 2aT 0 (42) 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. (43) Int. J. Numer. Meth. Engng. 44, 1823–1841 (1999) INFLUENCE OF STOCHASTIC MATERIAL DEFECTS 1833 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 = 1 (@i De U + De @i U + t=@i bp − t=(@i t b − @i t bp )) 1 + t= (44) @i = 1 (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. NUMERICAL SIMULATIONS 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 ) (45) 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) 1834 M. A. GUTIERREZ AND R. DE BORST 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 ) (46) 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 ) (47) where VM is the expectation of V and F is any matrix such that FT F = (Cov[V])−1 (48) 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) INFLUENCE OF STOCHASTIC MATERIAL DEFECTS 1835 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) 1836 M. A. GUTIERREZ AND R. DE BORST 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 ) (49) 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) INFLUENCE OF STOCHASTIC MATERIAL DEFECTS 1837 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) 1838 M. A. GUTIERREZ AND R. DE BORST 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) INFLUENCE OF STOCHASTIC MATERIAL DEFECTS 1839 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) 1840 M. A. GUTIERREZ AND R. DE BORST 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. 5. CONCLUSIONS 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. REFERENCES 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) INFLUENCE OF STOCHASTIC MATERIAL DEFECTS 1841 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, 1984. 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 (1978). 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)

1/--страниц