close

Вход

Забыли?

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

?

524

код для вставкиСкачать
INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING
Int. J. Numer. Meth. Engng. 42, 289–311 (1998)
A COMPUTATIONAL MODEL FOR ELASTO-VISCOPLASTIC
SOLIDS AT FINITE STRAIN WITH REFERENCE TO
THIN SHELL APPLICATIONS
∗
M. ROUAINIA AND D. PERIC
Department of Civil Engineering, University of Wales Swansea, Singleton Park, Swansea SA2 8PP, U.K.
ABSTRACT
This work extends a previously developed methodology for computational plasticity at nite strains that is
based on the exponential map and logarithmic stretches to the context of isotropic elasto-viscoplastic solids.
A particular form of the strain-energy function, given in terms of its principal values is employed. It is
noticeable that within the proposed framework, the small strain integration algorithms, and the corresponding
consistent tangent operators, automatically extend to the nite strain regime. Central to the eort of this
formulation is the derivation of the closed form of a tangent modulus obtained by linearization of incremental non-linear problem. This ensures asymptotically quadratic rates of convergence of the Newton–Raphson
procedure in the implicit nite element solution. To illustrate the performance of the presented formulation,
several numerical examples, involving failure by strain localization and nite deformations, are given. ? 1998
John Wiley & Sons, Ltd.
KEY WORDS: viscoplasticity; nite strain; localization; consistent linearization; nite elements
1. INTRODUCTION
Numerical modelling of nite strain plasticity problems has emerged in the last two decades or
so as a crucial part of the industrial design with a wide variety of applications involving complex
constitutive relations. From a purely theoretical standpoint, numerous technical papers have been
published, tackling a broad range of theoretical issues such as the notion of unstressed intermediate
conguration, material frame invariance, the role of the plastic spin, decomposition of strain rates
and related aspects, etc. (see References 1–4 among others). On the computational side, eort has
been directed toward the formulation of algorithms for integration of constitutive equations relying
mostly on the so-called operator-split methodology.5
The purpose of this paper is the formulation and numerical implementation of the constitutive
model for large deformations of elasto-viscoplastic solids, within the framework of the multiplicative decomposition.1 The elasto-viscoplastic constitutive equations are derived with reference to
the intermediate conguration. In addition, the present formulation is restricted to isotropic yield
criterion. The logarithmic stretches are employed as strain measures in conjunction with an exponential approximation of the viscoplastic ow rule, in the form rst proposed by Weber et al.6 As
∗
Correspondence to: D. Peric. Department of Civil Engineering, University of Wales Swansea, Singleton Park, Swansea
SA2 8PP, U.K. E-mail: D.Peric@Swansea.ac.uk
CCC 0029–5981/98/020289–23$17.50
? 1998 John Wiley & Sons, Ltd.
Received 26 November 1996
Revised 27 August 1997
290
M. ROUAINIA AND D. PERIC
an immediate consequence, the eects of nite strains appear only at kinematic level and many
properties of small strain algorithms are recovered. In particular, the standard return mapping
algorithm and the corresponding consistent tangent operators of the small strain situations can be
directly employed.
In this work, the eort is directed at the formulation of a complete computational model for nite
strain elasto-viscoplastic solids and the derivation of the associated consistent tangent modulus.
This is required by the Newton–Raphson solution of the associated incremental boundary value
problem. Here, the explicit expression of the tangent modulus is obtained by using the concept
of the consistent linearization of the weak form of the momentum balance equations. The key
result is the achievement of the asymptotically quadratic rate of convergence of the equilibrium
iterations.
Accounting for nite strain rate-independent elasto-plastic deformations in numerical simulations
often leads, as a result, to the emergence of strain localization phenomena. These phenomena,
which appear as a narrow region of intense deformations, frequently form a precursor to failure
mechanism. Numerical investigations of strain localization have shown that the solutions exhibit
strong mesh dependency. However, it is commonly known that the use of the rate-dependent
viscoplastic solids as regularization mechanism, preserves the ellipticity of the governing equations
and provides a mesh-objective solution.7 – 9
The remainder of the paper is organized as follows: Section 2 dispenses some preliminary
notions of kinematics, the essential assumptions and summarizes the constitutive equations based
on the Perzyna formulation.10 Issues pertaining to the algorithm for numerical integration of the
constitutive equations are addressed in Section 3 along with the small strain consistent tangent
operators. The closed form of the consistent tangent modulus, arising from the linearization of the
principle of virtual work, is derived in Section 4. An assessment of the eectiveness of the proposed
formulation is provided in Section 5. Typical problems, stretching of a perforated plate, pinched
hemispherical shell, Scordelis–Lo roof subjected to its own weight and plane strain stretching of
a thin sheet by a square punch, form the basis of this assessment. Finally, some conclusions are
drawn in Section 6.
2. FINITE STRAIN ELASTO-VISCOPLASTIC CONSTITUTIVE MODEL
2.1. Kinematical preliminaries
Denote by 0 ∈ R3 a reference conguration of the continuum body B with boundary @
0
and I ⊂ R+ the time interval of interest. At time t0 , each material point in 0 is labelled by
the co-ordinates X. The motion of B is described by e : 0 × I → R3 which maps, at the time
instant t, a material point onto x = e(X; t) of the current or deformed conguration. The associated
deformation gradient of the motion e is dened by:
F=
@e(X; t) @x
=
@X
@X
(1)
The kinematics of large elastoplasticity is based on the concept of multiplicative split of the
deformation gradient, F, into elastic and plastic parts:1
F = FeFp
? 1998 John Wiley & Sons, Ltd.
(2)
Int. J. Numer. Meth. Engng. 42, 289–311 (1998)
COMPUTATIONAL MODEL FOR ELASTO-VISCOPLASTIC SOLIDS
291
This choice introduces the concept of a local, unstressed or relaxed intermediate conguration B t .
One can obtain B t from the current conguration by a purely elastic unloading of a neighbourhood
of a material particle. Under the assumption of inelastic incompressibility, i.e. det[F p ] = 1, and since
det[F]¿0, it follows from (2) that det[F e ]¿0. Thus, from the polar decomposition theorem, F
and F e admit the unique representation
F = RU;
Fe = ReU e
(3)
where the right total and elastic stretch tensors, U and U e , respectively, are positive-denite
symmetric and the rotation tensors R and R e are proper orthogonal tensors. U and U e are related
to the right total and elastic Cauchy–Green deformation tensors, C and C e by
U 2 = C = FT F;
(U e ) 2 = C e = (F e )T F e
(4)
In relation to the constitutive model which is described in the next section, the Hencky logarithmic
strain E e , is used as the elastic Lagrangian strain measure
E e = ln [U e ]
(5)
Its conjugate rotated stress measure, which will be used in the description of the isotropic elastic
law, is denoted T := (R e )T c R e , where c denotes the Kirchho stress tensor. A detailed discussion
of stress conjugate to logarithmic strain adopted here can be found in Reference 11.
The time dierentiation of expression (1) and the use of the chain rule, yield the following
additive split form of the velocity gradient, L:
L = ḞF−1 = L e + L p
(6)
where L e and L p , the elastic and plastic velocity gradients, are, respectively, given by
L e = Ḟ e (F e )−1 ;
L p = F e Ḟ p (F p )−1 (F e )−1
(7)
p = Ḟ p (F p )−1 is used to derive the evolution of
Furthermore, the ‘modied’ plastic contribution L
p
F and, therefore, the intermediate conguration B t .
2.2. Isotropic elastic law
Following Reference 12, we assume the existence of a quadratic strain energy function W in
e
(i = 1; 2; 3). In fact, it is generally
the form of a scalar symmetric function of its stretches (i)
accepted that a suciently general constitutive model with possibly wide range of applications,
including metal plasticity, may be dened by employing the strain function W . This is given by
W (1e ; 2e ; 3e ) = [ln(1e ) 2 + ln(2e ) 2 + ln(3e ) 2 ] + 12 ln(J e ) 2
(8)
where and are positive material constants and (J e ) = 1e 2e 3e is the Jacobian. After applying
the standard procedure, the following hyperelastic constitutive equation is obtained:
T=
@W
= He : Ee
@E e
(9)
where H e = 2l + (I ⊗ I) denotes the fourth-order isotropic constant elastic tensor, with I and l
given in the component form as Iij = ij and Iijkl = 12 (ik jl + il jk ).
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 42, 289–311 (1998)
292
M. ROUAINIA AND D. PERIC
2.3. Viscoplastic constitutive equations
This section provides a brief review of the equations governing a material whose elastoviscoplastic model is based on a Von-Mises yield criterion and power law isotropic hardening.10
The model relies on the existence of the viscoplastic potential (T; T ) where (T ) : 0 × I → R is
the uniaxial yield stress conjugate to the equivalent viscoplastic strain E vp ∈ R. Since the attention
is restricted to the isotropic case and a single scalar hardening variable with linear hardening law,
the constitutive equations can be represented by
f(T; T ) = F(T) − T
p
F(T) = 3J2 (T)
(10)
(11)
T (E vp ) = 0 + H E vp
(12)
sym[Ḟ p (F p )−1 ] = h(T; T )i@T f(T; T )
(13)
E˙ vp = T˙ H −1 = h(T; T )i
(14)
where H is the standard isotropic linear hardening parameter and 0 is the initial uniaxial yield
stress. In equation (13), ∈ (0; ∞) is the uidity parameter, h·i is the ramp function dened as
hxi = (x + |x|)=2 and @T f(T; T ) = @f(T; T )=@T denotes the gradient of the yield function. Among
the various choices for the viscoplastic ow potential, (T; T ), the following form for isotropic
hardening is adopted:
N
F(T)
(T; T ) =
−1
(15)
T
where N is material parameter representing the rate sensitivity of the material. The evolution
problem described by (10)–(14) and (15), has a rm experimental basis and is widely accepted
as a description of rate-dependent deformations of solids. Additional formulative details of the
elasto-viscoplastic model used in the present study are given in Reference 13.
It is emphasized that the evolution of the plastic deformation gradient F p ( ow rule) is dened
by assuming the hypothesis of normality of the viscoplastic ow expressed by equation (13).
3. NUMERICAL INTEGRATION ALGORITHM
Assuming that an incremental displacement u is given, the problem consists of formulating
numerical procedures for updating the known variables {Fne ; Fnp ; Tn ; E vp
n } at time tn to obtain
p
e
; Fn+1
; Tn+1 ; E vp
}
at
time
t
.
For
this
purpose,
the
widely
used
algorithm based on oper{Fn+1
n+1
n+1
ator split methodology5 which involves two basic steps within a load increment, namely elastic
predictor and plastic corrector, is adopted. In the rst step, a purely elastic trial stress is computed;
if the yield condition is violated, the visco-plastic correction step based on the Newton–Raphson
procedure (in terms of (Tn+1 ; E vp
n+1 )) is performed on the xed intermediate conguration, such
that the nal stress is fully consistent with the constitutive integration algorithms. It should be
pointed out, however, that for a low rate-sensitive material (N ¿5) which causes the viscoplastic
evolution problem to be sti, these algorithms fail to converge for moderately large load steps.
In this case, a special numerical treatment is used consisting in taking the values (T̃n+1 ; E˜ vp
n+1 ) of
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 42, 289–311 (1998)
COMPUTATIONAL MODEL FOR ELASTO-VISCOPLASTIC SOLIDS
293
Box 1. Algorithm for integration of elasto-viscoplastic constitutive equations
(i) Initialization
(0)
= Fn
Fn+1
(ii) Update current conguration and deformation gradient
Fu = 1 + Grad u;
(i)
Fn+1 = Fu Fn+1
(iii) Evaluate elastic trial logarithmic strain
e trial
= Fn+1 (Fnp )−1
Fn+1
e trial
e trial T e trial
Cn+1
= (Fn+1
) Fn+1
e trial
e trial T e trial 1=2
Un+1
= [(Fn+1
) Fn+1 ]
e trial
e trial
e trial −1
Rn+1
= Fn+1
(Un+1
)
e trial
e trial
En+1
= ln[Un+1
]
(iv) Perform the small strain update : GOTO Box 2
(v) Update the rst Piola–Kirchho stress and intermediate conguration.
e
e trial
= Rn+1
Rn+1
e
e
Pn+1 = Rn+1
Tn+1 (Rn+1
)T (Fn+1 )−T
p
Fn+1
= exp[ t n+1 @T fn+1 ]Fnp
the rate-independent solution as the initial conditions to evaluate the viscoplastic response.13 This
procedure is additionally complemented by a line search method (see References 14 –16).
The governing equations for one load step of the Newton–Raphson iterative procedure for the
elasto-viscoplastic material are summarized in Boxes 1 and 2.
As a crucial point in the derivation of the algorithms, an exponential approximation is used
in the discretization of the ow rule in the plastic corrector step rst suggested by Weber and
Anand6 (see also References 17–19, 12).
p
Fn+1
= exp [t n+1 @T fn+1 ] Fnp
(16)
e trial
e trial T
Furthermore, if the material is assumed to be isotropic, with the consequence that Cn+1
= (Fn+1
)
e trial
Fn+1 and @T fn+1 commute, it can be easily shown that the above approximation (16) results in
e
e trial
En+1
= En+1
− t n+1 @T fn+1
? 1998 John Wiley & Sons, Ltd.
(17)
Int. J. Numer. Meth. Engng. 42, 289–311 (1998)
294
M. ROUAINIA AND D. PERIC
Box 2. Stress updating procedure-Small strain
(i) Initialize
k = 0;
(0)
E vp
n+1 = 0;
(0)
trial
Tn+1
= Tn+1
(ii) Check yield condition and evaluate residuals
(k)
(k) vp (k)
= f(Tn+1
; T (E n+1 ))
fn+1
e −1
(k)
(k) trial − Tn+1
)
(H ) : (Tn+1
@T fn+1
(k)
(k)
+ tn+1
Rn+1 =
(k)
vp
−1
(E vp
n+1 − E n )
(k)
(k)
IF k Rn+1
k ¡R and fn+1
¡f THEN
Set (•)n+1 = (•)(k)
n+1 and EXIT
ELSE
(iii) Compute algorithmic modulus
(k)
n+1
−1

=
vp
TT
n+1
E
Tn+1
E vp T
n+1
E vp E vp
n+1
−1

where
(k)
(k)
(k)
(k)
e −1
2
+ t(k)
TT
n+1 = (H )
n+1 @TT fn+1 + t@F n+1 @T fn+1 ⊗ @T fn+1
vp
(k)
E
Tn+1
= t@Evp (k)
n+1 @T fn+1
vp
(k)
E T
n+1
= t@F (k)
n+1 @T fn+1
vp vp
E E
n+1
= 1 − t@Evp (k)
n+1
(iv) Evaluate increment of stresses and internal variable
(k+1) −1
Tn+1
(k)
(k)
: Rn+1
vp (k+1) = − n+1
E n+1
(v) Update stresses and internal variable
(k+1) (k) (k+1) Tn+1
Tn+1
Tn+1
vp (k+1) =
vp (k) +
vp (k+1)
E n+1
E n+1
E n+1
Set
? 1998 John Wiley & Sons, Ltd.
k =k + 1
and GOTO (ii)
Int. J. Numer. Meth. Engng. 42, 289–311 (1998)
COMPUTATIONAL MODEL FOR ELASTO-VISCOPLASTIC SOLIDS
295
One of the most striking features of the adoption of the logarithmic strain measure along with the
exponential approximation for the viscoplastic ow rule, is the preservation of the standard return
algorithms of the small strain theories. It should be noted that, as consequence of the exponential
approximation, the incompressibility viscoplastic ow for the pressure-insensitive yield criterion is
also preserved.
3.1. Small strain tangent operator
From the computational standpoint, an important issue is the derivation of the small strain
tangent operator in a manner consistent with the constitutive integration algorithm of Box 2. Simo
and Taylor have shown, in their pioneering work20 concerning plasticity, that such an operator
preserves the important property of asymptotic quadratic convergence in the Newton–Raphson
iterative procedure. By contrast, the continuum modulus is known to destroy this important property
and frequently leads to a divergent iterative process, particularly when the load step size is larger.
After rather lengthy but straightforward algebraic manipulation, the following explicit expression
b of the present elasto-viscoplastic model, is obtained
for the small consistent tangent operator H,
(see Reference 13 for details):
h
i−1
b = @Tn+1 = TT − TEvp (Evp Evp )−1 Evp T
(18)
H
n+1
n+1
n+1
n+1
e
@En+1
vp
vp
vp vp
E T
TE
E E
where TT
n+1 , n+1 , n+1 and n+1
are given in Box 2.
4. LINEARIZED WEAK FORMULATION: CONSISTENT TANGENT MODULUS
The full Newton–Raphson iteration scheme is characterized by a convergence rate that is asymptotically quadratic. In order to preserve this feature, it is crucial that numerical models for the
description of general non-linear behaviour of solids are consistently linearized.
At time tn+1 , the weak form of the momentum balance, known as the principle of virtual work,
can be written, with respect to the initial conguration, as
Z
Z
G(en+1 ; W) =
(Pn+1 : ∇0 W − bn+1 ·W) d
0 −
tn+1 ·W d 0 = 0; ∀W ∈ V
(19)
0
@
0
where P = c·(F)−T is the rst Piola–Kirchho stress, fn+1 and tn+1 are, the body forces and the
surface tractions specied on the portion @
0 of the boundary, respectively. The symbol ∇0 denotes
the material gradient and V is the space of the admissible virtual displacements. The linearization
of (19) provides the basis of the Newton-type iterative schemes for solution of the momentum
balance equations. In this context, the directional derivative (see Reference 16) of the internal
virtual work plays an essential role. The virtual work of the internal forces
Z
Pn+1 : ∇0 W d
0
(20)
G int (en+1 ; W) =
0
has its derivative in the direction of an incremental displacement u given by the formula
Z
DG int (en+1 ; W)·u =
∇0 W : Dn+1 : ∇0 (u) d
0
(21)
0
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 42, 289–311 (1998)
296
M. ROUAINIA AND D. PERIC
where Dn+1 is the consistent tangent modulus given by
Dn+1 =
@Pn+1
@Fn+1
(22)
Remark 4.1. It should be observed that the nite strain constitutive model described in this
section is formulated in the material framework. This is in contrast to the model described by
de Souza Neto and Peric21 which is phrased in spatial co-ordinates. Although both models ultimately lead to identical results, the present model is more convenient for description of large
deformations of bodies which show certain preferred directions. Typical problems from this class
arise in modelling of nitely deformed shells.
4.1. Exact consistent tangent modulus
In order to obtain a step by step derivation of the exact expression of the consistent tangent
modulus, Dn+1 , it is convenient to rewrite the rst Piola–Kirchho stress Pn+1 in the alternative
form
Pn+1 = c·(Fn+1 )−T = Fn+1 ·Qn+1
(23)
e trial −1
e trial −T
Qn+1 = (Fnp )−1 (Un+1
) Tn+1 (Un+1
) (Fnp )−T
(24)
where Qn+1 is given by
The consistent tangent modulus in component form, follows straightforwardly as
Dijkl = ik Qlj + (Fn+1 )im
@(Qn+1 )mj
@(Fn+1 )kl
(25)
Since Qn+1 is a function of Tn+1 , the rotated stress tensor has to be derived with respect to the
total deformation gradient Fn+1 . This is accomplished by application of the chain rule as follows:
e trial
@(Cn+1
)
@(Tn+1 ) b
=H:L:
@(Fn+1 )
@(Fn+1 )
(26)
b is the small strain consistent tangent operator given by equation (18) and the fourth-order
where H
tensor L is dened as
L :=
e trial
@(En+1
)
e trial
@(Cn+1 )
(27)
e trial
. AddiUsing equation (5), the explicit form for L is the derivative of the tensor logarithm of Un+1
tionally, the remaining derivative term on the right-hand side of equation (26) follows easily from
e trial
e trial T e trial
:= Fn+1
·Fn+1 = (Fnp )−T ·
the algorithmic denition of the elastic trial Cauchy–Green tensor Cn+1
p
T
·Fn+1 ·(Fn )−1 , i.e.
Fn+1
e trial @(Cn+1
)
e trial
e trial
= (Fnp −1 )li (Fn+1
)kj + (Fn+1
)ki (Fnp −1 )lj
(28)
@(Fn+1 ) ijkl
To completely determine equation (25), it still remains to obtain an explicit formula for the
e trial
, with respect to the total deformation gradient
derivative of the right trial stretch tensor, Un+1
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 42, 289–311 (1998)
COMPUTATIONAL MODEL FOR ELASTO-VISCOPLASTIC SOLIDS
297
Fn+1 . Applying again the chain rule, yields
e trial −1
e trial
@(Un+1
@(Cn+1
)
)
=M:
@(Fn+1 )
@(Fn+1 )
(29)
in which the components of the right-hand side term are given in (28), and M is dened as
M :=
e trial −1
@(Un+1
)
e trial
@(Cn+1 )
(30)
e trial
e trial 1=2
By exploiting the fact that Un+1
= [Cn+1
] , the fourth-order tensor M is the derivation of
e trial
the inverse of the square root of Cn+1 . For more information regarding the computation of the
tensors L and M, as case of the derivative of general tensor function, the reader is referred to
References 21 and 22.
Finally, with all expressions involved in equation (25) at hand, the explicit formula of the
Cartesian components of the consistent tangent modulus follows as
b : L + K : M) : N]ijkl
Dijkl = ik Qlj + [(R : H
(31)
Here the fourth-order tensors R, K and N are conveniently dened as
e
e trial
e trial −1
Kijkl = [(Rn+1
)(Tn+1 )]il (Fnp −1 )jk + (Fn+1
)ik (Tn+1 )lm [(Fnp −1 )(Un+1
) ] jm
e
e trial −1
Rijkl = (Rn+1
)ik [(Fnp −1 )(Un+1
) ] jl
e trial
Nijkl = (Fnp −1 )li (Fn+1
)kj
+
(32)
e trial
(Fn+1
)ki (Fnp −1 )lj
Examination of the closed formula (31) clearly shows that the symmetric consistent tangent modulus D, contains both geometric and material nonlinearities contributions. The material related
b The remaining terms are
contributions appear only through the small strain consistent tangent H.
related exclusively to the kinematics of large strains.
5. NUMERICAL EXAMPLES
To demonstrate the performance of the elasto-viscoplastic model in conjunction with the numerical
scheme discussed earlier, numerical examples are presented in this section. All numerical simulations are performed by using an incremental solution based on the backward Euler time-stepping
scheme described in Sections 3 and 4 and the Newton–Raphson procedure. It is noted that the
convergence properties of the Newton–Raphson procedure are ensured by using the consistent tangent operator, given in Sections 3.1 and 4.1. The convergence of the incremental nite element
solution is established on the basis of the standard Euclidean norm of the out-of-balance forces.
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 42, 289–311 (1998)
298
M. ROUAINIA AND D. PERIC
Figure 1. Perforated plate subjected to stretching: geometry, nite element mesh and material properties
5.1. Stretching of a perforated plate
The problem consists of a two-dimensional thin plate (plane stress)† subjected to stretching. Due
to symmetry of the problem under consideration, only one quarter of the plate with appropriate
boundary conditions is modeled. The geometry and material properties are shown in Figure 1.
The half-length of the perforated plate is L0 = 18·0, the half-width is B0 = 10·0; the radius of
the hole is R0 = 5·0 and the thickness is t0 = 1·0. The material is characterised by an isotropic
hardening viscoplastic response. The following properties are assigned to the plate: elastic modulus
is E = 70, Poisson’s ratio is 0·2, yield stress is 0 = 0·243 and hardening parameter H = 0·2. The
nite element idealization consists of 613 three-noded constant strain membrane elements. Loading
is performed by controlling the displacement U2 at the top of the plate.
Figure 2 shows plots of the load versus edge displacement for dierent values of the rate sensitivity exponent N = 1; 10 and 100. The average deformation rate investigated ranges between
(V2 =L0 ) = 0·555 × 10−4 (1=s) and (V2 =L0 ) = 5·555 × 10−1 (1=s) and the rate being increased by a
factor of ten for every new simulation. After the peak load, the load-displacement curve drops
steeply. For high rate-sensitive material, N = 1, the plots highlight the eects of the deformation
rate on the increase of the maximum value of the force. For low rate-sensitive material, N = 100;
and for the whole range of the deformation rates, only a small dierence in the force–displacement
is observed. These plots also show that the force–displacement obtained using a very low rate
of deformation (V2 =L0 ) = 0·555 × 10−4 (1=s) is nearly equal to the rate-independent solution.
Figures 3 and 4 show evolution of the viscoplastic zone in the plate together with the corresponding typical deformed meshes at three dierent load levels. It can be noticed that the localization
†
The constitutive equations and associated numerical scheme are developed directly in the plane stress form with
33 = 32 = 31 = 0, as described in Reference 15
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 42, 289–311 (1998)
COMPUTATIONAL MODEL FOR ELASTO-VISCOPLASTIC SOLIDS
299
Figure 2. Stretching of a perforated plate: load–displacement curves for various rates of deformation and several values of
the rate sensitivity exponent
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 42, 289–311 (1998)
300
M. ROUAINIA AND D. PERIC
Figure 3. Stretching of a perforated plate: evolution of viscoplastic zone, plotted over initial conguration, at various stages
of loading: (a) U2 = 1·15; (b) U2 = 3·65; (c) U2 = 6·15
Figure 4. Stretching of a perforated plate: deformed nite element meshes for dierent load steps: (a) U2 = 1·15;
(b) U2 = 3·65; (c) U2 = 6·15
of plastication, which takes place around the hole, results in further narrowing of the weakest
section of the plate and forms a precursor to failure.
This simulation is successfully completed in 48, 82 and 216 load steps for N = 1; 10 and 100,
respectively. A standard check of the residual forces, based on the Euclidean norm, was carried
out with a convergence tolerance of 1 × 10−6 . The number of iterations required to meet this
tolerance was about 3–6 iterations for each load step. As expected, the convergence rate of the
Newton–Raphson method clearly exhibits the quadratic rates of asymptotic convergence whenever
the solution is within the radius of convergence. This fact is presented in tabular form for four
typical load steps and the rate sensitivity of material N = 100 in Table I. For the same load steps,
this fact is further elucidated graphically in Figure 5 where the residual norm kRk k of the iteration k
is plotted against the residual norm kRk+1 k of the subsequent iteration k +1. Observe that the slope
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 42, 289–311 (1998)
301
COMPUTATIONAL MODEL FOR ELASTO-VISCOPLASTIC SOLIDS
Table I. Perforated plate: residual norms for four typical load steps for N = 100 and average
rate of deformation V2 =L0 = 0·555 × 10−2 (1=s)
Iter
1
2
3
4
5
6
U2 = 0·15
U2 = 0·65
U2 = 1·15
U2 = 3·65
2·046703E−01
0·351815E−01
0·146249E−02
0·721279E−05
0·182781E−09
—
1·458925E−01
0·566755E−01
0·462543E−02
0·141538E−03
0·154090E−06
0·184024E−10
0·732901E−01
0·370454E−01
0·267615E−02
0·623028E−05
0·320327E−09
—
0·776389E−01
0·759442E−01
0·442407E−03
0·901602E−07
—
—
Figure 5. Stretching of a perforated plate: norms of the residual for four typical load steps corresponding to the rate
sensitivity exponent N = 10
1 : 2 shown in the gure corresponds to a quadratic convergence. As mentioned earlier, these results
were achieved using the closed-form consistent tangent operator given in Sections 3.1 and 4.1.
5.2. Pinched hemispherical shell
A pinched hemispherical shell loaded by two concentrated forces, 90◦ apart and characterized
by an isotropic hardening, is considered. Again, due to symmetry of the problem, only one quarter
of the hemispherical shell is considered in the nite element simulation. The geometry, the nite
element mesh, consisting of 128 triangular Morley thin shell elements, and the material properties
are shown in Figure 6. The radius of the shell is R = 10 and the thickness is t0 = 0·5. The material
is modelled with the following parameters : elastic modulus is E = 10, Poisson’s ratio is = 0·2,
yield stress is 0 = 0·02 and hardening parameter is assumed to be equal to H = 6·0. This example has been analysed by Simo and Kennedy23 by considering the stress resultant shell theory of
elasto-plasticity.
Figure 7 shows the sequence of deformed conguration corresponding to dierent values of
imposed concentrated force P = 0·03 × 10−01 , 0·10 × 10−01 , 0·20 × 10−01 and 0·30 × 10−01 . It is
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 42, 289–311 (1998)
302
M. ROUAINIA AND D. PERIC
Figure 6. Pinched hemispherical shell: geometry, nite element model and material parameters
Figure 7. Pinched hemispherical shell: deformed nite element meshes at various steps of loading: (a) P = 0·03 × 10−01 ;
(b) P = 0·10 × 10−01 ; (c) P = 0·20 × 10−01 ; (d) P = 0·30 × 10−01
worth noticing that the displacement of the node 289 (inward) in the nal conguration is almost
equal to the radius of the hemispherical shell.
Figure 8 shows the load–displacement curves for dierent uidity parameter ranging from
0·2 × 10−6 (1=s) to 0·2 × 10−2 (1=s). The range of the rate sensitivity exponent investigated is
N = 1; 10 and 100. Once again, the inuence of the material rate sensitivity is demonstrated by
increase of the force for the high rate-sensitive material represented by N = 1. Only a minor difference in the overall force–displacement curves is observed for the low rate-sensitive material
characterized by N = 100 over the whole range of the uidity parameter .
This simulation is successfully completed in 21 load steps. The number of iterations required
to meet the convergence tolerance was about 3– 4 iterations for each load step. Results of the
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 42, 289–311 (1998)
COMPUTATIONAL MODEL FOR ELASTO-VISCOPLASTIC SOLIDS
303
Figure 8. Pinched hemispherical shell: load–displacement curves for various values of the uidity parameter and several
values of the rate sensitivity exponent
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 42, 289–311 (1998)
304
M. ROUAINIA AND D. PERIC
Table II. Pinched hemispherical shell: residual norms for four typical load steps for N = 10
and uidity parameter = 0·2 × 10−3 (1=s)
Iter
P = 0·5 × 10−2
P = 0·10 × 10−1
P = 0·24 × 10−1
P = 0·30 × 10−1
1
2
3
4
0·410084E−01
0·389172E−02
0·138051E−03
0·532599E−06
0·221076E−01
0·142509E−02
0·873790E−03
0·251861E−05
0·281587E−01
0·104769E−02
0·611148E−03
0·842952E−06
0·191048E−01
0·120899E−02
0·249249E−03
0·219909E−06
Figure 9. Scordelis–Lo roof: geometry, material parameters and loading conditions
Euclidean norm of the global residual are shown in Table II and suggest that, again, a quadratic
rate of convergence is attained by the Newton–Raphson iterative procedure.
5.3. Elasto-viscoplastic Scordelis–Lo roof subjected to gravity loading
A Scordelis–Lo roof composed of elasto-viscoplastic material and subjected to gravity loading
is next considered. Due to symmetry considerations, only one-quarter of the roof, with appropriate boundary conditions, has been discretized. The geometry and the material properties for
this example are shown in Figure 9. The half-length of the roof is L = 7·6 m, radius is R = 7·6 m
and thickness is t0 = 0·076 m. The selected material properties of the roof are: elastic modulus is
E = 2·1 × 104 N=mm 2 , Poisson’s ration is = 0·2, yield stress is 0 = 4·2 N=mm 2 and hardening parameter H = 0·0. The gravity type loading is performed with the reference value of q0 = 4 kN=m 2 .
The deformations are restricted to be symmetric along the line X1 = 0 and X2 = 0. Accordingly,
the boundary conditions can be expressed as
U1 = 0
2 = 0
at X1 = 0
U2 = 0
1 = 0
at X2 = 0
U1 = U3 = 0
2 = 0
at X2 = L
where Ui are displacements in Xi directions and j are rotations about Xj -axis.
The gravitational force q=q0 applied to the roof is plotted against the vertical displacement of
point A in Figure 10, obtained with a mesh of 12 × 12 triangular Morley elements. Simulations are
carried out for rate sensitive material exponent N = 1 and for three dierent values of the uidity
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 42, 289–311 (1998)
COMPUTATIONAL MODEL FOR ELASTO-VISCOPLASTIC SOLIDS
305
Figure 10. Scordelis–Lo roof: Load-displacement curves for point A
parameter = 2·0 × 10−2 ; 2·0 × 10−3 and 2·0 × 10−4 . For comparison purposes, the underlying
rate-independent solution obtained by Peric et al.24 is also shown. From this example, the following
observations, already drawn from the previous tests, may be noted again: (a) the major inuence
of the rate of deformation, which is altered by varying the uidity parameter , is demonstrated by
the eect on the maximum value of the gravity force; (b) the descent of the force-displacement
curve observed in the rate-independent case is completely averted for the higher rate sensitive
material, N = 1.
For the rate-independent solution (see Figure 5.8 in Reference 24), a typical localization behaviour, as indicated by the sharp drop in the force-displacement curve, starts through the initiation
of plastic hinge around the X2 axis, followed by a second hinge in the X1 direction, leading to
failure. Not surprisingly, the failure mode is quite dierent in the rate-dependent case. Typical
deformed meshes of the roof in Figure 11 show that the appearance of localization and associated
failure is less pronounced with the increase of the rate sensitivity. This is evident from Figure 11,
where the formation of the second plastic hinge in X1 direction across the mid-section of the roof
is clearly delayed in comparison to the rate-independent results.
5.4. Plane strain stretching of a thin sheet by a square punch
A thin sheet stretched by a rectangular punch, with approximately plane strain conditions along
the centre section perpendicular to the longer side of the punch, is next considered. The analysis
is performed by employing a 3-D formulation, restricting the deformations to be symmetric along
the line X1 = 0 and imposing plane strain boundary conditions in the X2 direction.
The geometrical characteristics and material properties are shown in Figure 12. A mesh of
240 Morley thin shell nite elements is considered for the half of the blank, and 64 triangular
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 42, 289–311 (1998)
306
M. ROUAINIA AND D. PERIC
Figure 11. Scordelis–Lo roof: Deformed conguration of the shell at four dierent load steps: (a) q=q0 = 1·248; (b)
q=q0 = 1·688; (c) q=q0 = 2·008; (d) q=q0 = 2·408
Figure 12. Plane strain stretching of a thin sheet by a square punch: geometry and materials characteristics
facet elements are used to discretize both the surfaces of the die and the punch. The problem is
investigated for a coecient of friction f = 0·25 and uidity parameter = 0·002 s−1 .
Figure 13 shows plots of the punch force versus punch displacement for dierent rate sensitivity
exponent N = 1; 10 and 50. The loading program is altered by varying the punch speed between
Vp = 0·004 – 4 m=s, increasing the speed by a factor of ten for every new test. As expected, the major
inuence of the rate of loading, demonstrated by the eect on the maximum value of the punch
force, is observed for the highly rate-sensitive material represented by N = 1. No drop in the punch
force is observed in the simulation for the high rate sensitive material characterised by N = 1 and
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 42, 289–311 (1998)
COMPUTATIONAL MODEL FOR ELASTO-VISCOPLASTIC SOLIDS
307
Figure 13. Plane strain stretching of a thin sheet by a square punch: punch force versus punch displacement curves for
various punch speeds and several values of the rate sensitivity exponent
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 42, 289–311 (1998)
308
M. ROUAINIA AND D. PERIC
Figure 14. Plane strain stretching of a thin sheet by a square punch: Deformed nite element meshes at punch displacement
Dp = 4 mm, Dp = 8 mm, Dp = 12 mm, Dp = 16 mm and Dp = 20 mm
Figure 15. Plane strain stretching of a thin sheet by a square punch: distribution of the true strain plotted over the initial
conguration for punch speed Vp = 0·04 m=s: (a) top layer; (b) bottom layer
10, as seen in Figure 13. For low rate-sensitive material (N = 50), the force–displacement diagram
displays a drop after achieving a maximum value, and approaches the rate independent solution
for slow loading conditions (Vp = 0·4 × 10−2 m=s).
The deformed conguration for dierent values of the punch displacement are presented in
Figure 14. The distribution of the true strain in the X1 direction is shown in Figure 15 for the punch
speed Vp = 0·04 m=s and the rate sensitivity N = 1 and 10 for punch displacement Dp = 12 mm. For
comparison, the experimental results of Stoughton25 and the underlying rate-independent solution
are also presented. For the rate-independent solution, which is evidently in agreement with the
experimental results, the strain accumulates in a narrow band reaching high levels thus reecting
the tendency to localisation behaviour. However, the appearance of localisation is averted in the
higher rate-sensitive material characterized by N = 1.
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 42, 289–311 (1998)
309
COMPUTATIONAL MODEL FOR ELASTO-VISCOPLASTIC SOLIDS
Table III. Plane strain stretching of a thin sheet by a square punch: residual norms for four typical
load steps for N = 1 and uidity parameter = 0·2 × 10−2 (1=s)
Iter
1
2
3
4
5
Dp = 4 mm
Dp = 8 mm
Dp = 12 mm
Dp = 16 mm
0·711233E + 04
0·786184E + 03
0·443103E + 03
0·311620E + 01
0·457315E − 04
0·571463E + 04
0·870780E + 03
0·222233E + 03
0·700009E + 00
0·127989E − 05
0·640930E + 04
0·905270E + 03
0·136072E + 03
0·221860E + 00
0·231715E − 06
0·709472E + 04
0·619294E + 03
0·280641E + 03
0·110341E + 01
0·508533E − 05
Figure 16. Plane strain stretching of a thin sheet by a square punch: norms of the residual for ve typical load steps
In our numerical model, throughout the process of deformation, the contact region between
blank and punch has been restricted to a small area around the punch radius. This creates high
normal and frictional forces which suppress sliding of the sheet over the punch. This trend, which
is reported by Stoughton,25 is closely related to the bending eects which introduce an additional
restraint and highly inuence the overall deformation pattern.
It should be noted that, quadratic convergence is observed in a typical load step. As an illustration, values of the Euclidean norm of the global residual are reported in Table III and Figure 16
within several typical punch displacements. These results provide a clear indication of the quadratic
rates of asymptotic convergence near the solution point. The number of iterations needed to meet
the convergence tolerance of 1 × 10−3 is about 4 – 7 per increment.
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 42, 289–311 (1998)
310
M. ROUAINIA AND D. PERIC
6. CONCLUSIONS
Within the ecient hyperelastic-based multiplicative nite strain framework, a model for large
elasto-viscoplastic deformations has been proposed and applied to the large-scale problems involving failure by strain localization. A particularly simple algorithm for integration of the constitutive
equations, has been obtained by employing logarithmic strain as a kinematic variable and an exponential approximation of the viscoplastic ow rule. In this case the small strain format of the
stress updating procedure is preserved and the eects of nite strain appear only at the kinematic
level independently of the constitutive model.
On the basis of the consistent linearization of the incremental boundary value problem, the
explicit expression of the tangent modulus has been derived. It relies on the closed formula for
the derivative of a class of isotropic tensor functions. The performance of the described formulation
and the accuracy of the numerical results have been demonstrated for four representative problems.
It has been found that both, the material rate sensitivity and the strain rate have a pronounced
eect on the appearance of strain localisation and on the description of the failure mechanism.
The use of the consistent tangent operator given in Sections 3.1 and 4.1, has enabled quadratic
rates of convergence to be obtained in the equilibrium iterations. This increases the speed of the
computations by requiring only an average of three to six iterations for each load step in all the
numerical examples.
REFERENCES
1. E. H. Lee, ‘Elastic-plastic deformation at nite strains’, J. Appl. Mech., 36, 1– 6 (1969).
2. J. Mandel, ‘Equations constitutives et directeurs dans les milieux plastiques et viscoplastiques’, Int. J. Solids Struct.,
9, 725 – 740 (1973).
3. S. Nemat-Nasser, ‘On nite elasto-plasticity’, Int. J. Solids Struct., 18(10), 857 – 872 (1982).
4. J. C. Simo, ‘A framework for nite elastoplasticity based on maximum plastic dissipation and multiplicative
decomposition: Part I. Continuum formulation’, Comput. Meth. Appl. Mech. Engng., 66, 199 – 219 (1988).
5. J. C. Simo and T. J. R. Hughes, ‘General return mapping algorithms for rate independent plasticity’, in C. S. Desai
et al. (eds.), Constitutive Laws for Engineering Materials: Theory and Applications, Elsevier, Amsterdam, 1987,
pp. 221– 231.
6. G. Weber and L. Anand, ‘Finite deformation constitutive equations and a time integration procedure for isotropic,
hyperelastic-viscoplastic solid’, Comput. Meth. Appl. Mech. Engng., 79, 173 – 202 (1990).
7. B. Loret and J. H. Prevost, ‘Dynamic strain localisation in elasto(visco)-plastic solids, Part I: General formulation and
one dimensional examples’, Comput. Meth. Appl. Mech. Engng., 83, 247 – 273 (1990).
8. A. Needleman, ‘Material rate-dependence and mesh sensitivity in localisation problems’, Comput. Meth. Appl. Mech.
Engng., 67, 69 – 86 (1988).
9. L. J. Sluys, J. Block and R. de Borst, ‘Wave propagation and localisation in viscoplastic media’, in D. R. J. Owen
et al. (eds.), Proc. 3rd Int. Conf. on Computational Plasticity, Pineridge Press, Swansea, 1992, pp. 299 – 312.
10. P. Perzyna, ‘Thermodynamic theory of viscoplasticity’, in Advances in Applied Mechanics, Vol. 11, Academic Press,
New York, 1971, pp. 313 – 354.
11. A. Hoger, ‘The stress conjugate to logarithmic strain’, Int. J. Solids Struct., 23, 1645 –1656 (1987).
12. D. Peric, D. R. J. Owen and M. E. Honnor, ‘A model for nite strain elasto-plasticity based on logarithmic strains’,
Comput. Meth. Appl. Mech. Engng., 94, 35 – 61 (1992).
13. D. Peric, ‘On a class of constitutive equations in viscoplasticity: formulation and computational issues’, Int. J. Numer.
Meth. Engng., 36, 1365 –1393 (1993).
14. M. A. Criseld, Non-linear Finite Element Analysis of Solids and Structures. Vol. 1: Essentials, Wiley, Chichester,
1991.
15. M. Dutko, D. Peric and D. R. J. Owen, ‘Universal anisotropic yield criterion based on superquadratic functional
representation: Part I: algorithmic issues and accuracy analysis’, Comput. Meth. Appl. Mech. Engng., 109, 73 – 93
(1993).
16. J. C. Marsden and T. J. R. Hughes, Mathematical Foundations of Elasticity, Prentice-Hall, Englewood Clis, N.J.,
1983.
17. A. Cuitiño and M. Ortiz, ‘A material independent method for extending stress update algorithms from small strain
plasticity to nite plasticity with multiplicative kinematics’, Engng. Comput., 9, 437 – 451 (1992).
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 42, 289–311 (1998)
COMPUTATIONAL MODEL FOR ELASTO-VISCOPLASTIC SOLIDS
311
18. A. L. Eterovic and K. J. Bathe, ‘A hyperelastic based large strain elasto-plastic constitutive formulation with combined
isotropic-kinematic hardening using the logarithmic stress and strain measures’, Int. J. Numer. Meth. Engng., 30,
1099 –1114 (1990).
19. J. C. Simo, ‘Algorithms for static and dynamic multiplicative plasticity that preserve the classical return mapping
schemes of the innitesimal theory’, Comput. Meth. Appl. Mech. Engng., 99, 61–112 (1992).
20. J. C. Simo and R. L. Taylor, ‘Consistent tangent operators for rate-independent elastoplasticity’, Comput. Meth. Appl.
Mech. Engng., 48, 101–118 (1985).
21. E. A. de Souza Neto and D. Peric, ‘A computational framework for a class of fully coupled models for elastoplastic
damage at nite strains with reference to the linearization aspects’, Comput. Meth. Appl. Mech. Engng., 130, 179 –193
(1996).
22. D. E. Carlson and A. Hoger, ‘The derivative of a tensor-valued function of a tensor’, Quart. Appl. Math., 44(3),
409 – 423 (1986).
23. J. C. Simo and J. G. Kennedy, ‘On the stress resultant geometrically exact shell model. Part V. Nonlinear plasticity:
formulation and integration algorithm’, Comput. Meth. Appl. Mech. Engng., 99, 133 –177 (1992).
24. D. Peric and D. R. J. Owen, ‘The Morley thin shell nite element for large deformations problems: Simplicity
versus sophistication’, in N. Bicanic (ed.), Proc. Int. Conf. on Nonlinear Engineering Computations, Pineridge Press,
Swansea, 1991, pp. 121–142.
25. T. B. Stoughton, ‘Finite element modeling of 1008 AK sheet steel stretched over a rectangular punch with bending
eects’, in N. M. Wang et al. (eds.), Computer Modeling of Sheet Metal Forming Process, AIME, New York, 1985,
pp. 143 –159.
? 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 42, 289–311 (1998)
Документ
Категория
Без категории
Просмотров
2
Размер файла
336 Кб
Теги
524
1/--страниц
Пожаловаться на содержимое документа