close

Вход

Забыли?

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

?

152

код для вставкиСкачать
INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING
Int. J. Numer. Meth. Engng. 44, 155–176 (1999)
AN IMPROVED ITERATIVE METHOD FOR
LARGE STRAIN VISCOPLASTIC PROBLEMS
H.-L. CAO∗ AND M. POTIER-FERRY
Laboratoire de Physique et MÃecanique des MatÃeriaux, UMR CNRS 7554, UniversitÃe de Metz,
Ile du Saulcy, 57045 Metz Cedex, France
ABSTRACT
In this work, some techniques are proposed to improve the usual Newton–Raphson Method (NRM) used
in the numerical analysis of large strain viscoplastic problems. These techniques, based on the ÿrst-order
perturbation technique, allow to deÿne an adaptive step strategy and to improve the trial solution (guessed
solution) for the ÿrst iteration at each step. To assess the eciency of the proposed techniques, a number
of numerical examples are presented: necking of a circular bar, plane strain bending, and an axisymmetric
hydrostatic bulging process. Compared with the classical Newton–Raphson method, the proposed scheme
requires less matrix inversions, less time steps and less CPU time. For a typical axisymmetric hydrostatic
bulging process (606 degrees of freedom), the proposed method needs 5 times less matrix inversions, 1·7
times less time steps, and 4 times less CPU time. Copyright ? 1999 John Wiley & Sons, Ltd.
KEY WORDS:
Newton–Raphson method; perturbation technique; viscoplasticity; adaptive step
1. INTRODUCTION
In the numerical simulation of metal forming processes, the viscoplastic constitutive laws have
been widely employed in analysing large strain processes by the ÿnite element method.1–5 Because
of material and geometrical non-linearities, some iterative methods such as the Newton–Raphson
Method (NRM) are generally used. This kind of methods has two major diculties: the diculty
to choose an optimal step size and the diculty to estimate a good trial solution for the ÿrst
iteration at each step.
For the optimal step size, some ÿnite element codes (ABAQUS, etc) adopt a so-called automatic
step size, which is limited according to some criteria (the maximum admissible strain increment,
for example). This technique permits the avoidance of very large steps when non-linearities come
essentially from geometric variations or from plastic strain increments. It does not, however, allow
to choose an adequate step size when the non-linearity is due to the variation of strain rate. For
the trial solution, Chung and Wagoner4 and many other authors1–3 have examined the choice
of dierent trial solutions. The major idea is to utilize the solution of the previous step or the
solution extrapolated from that of several previous steps. As these techniques are based only on
∗
Correspondence to: H.-L. Cao, LPMM, URA CNRS 1215, UniversitÃe de Metz, Ile de Saulcy, 57045 Metz Cedex, France.
E-mail: cao@lpmm.univ-metz.fr
CCC 0029–5981/99/020155–22$17.50
Copyright ? 1999 John Wiley & Sons, Ltd.
Received 17 October 1997
Revised 20 April 1998
156
H.-L. CAO AND M. POTIER-FERRY
the solutions of the previous steps, it is dicult to predict a good trial solution particularly when
the strain rate varies signicantly.
The perturbation technique, coupled with the discretization method, has been applied to numerical
analysis of geometrically non-linear problems.6 – 8 Some signicant advances on the perturbation
technique, called incremental perturbation method or Asymptotic Numerical Method (ANM), have
been achieved by Yokoo et al.,9 Damil and Potier-Ferry,10 as well as by Cochelin et al.11 The
same method has been extended to the large strain viscoplastic problems by Potier-Ferry et al.12
Our objective in this paper is to use the same perturbation technique to improve the classical
Newton–Raphson Method. The improvements concern particularly the two principal diculties
cited above.
Recently, a study carried out by Mole et al.13 employed a rst-order perturbation technique
to solve large strain viscoplastic problems, although the terminology ‘perturbation technique‘ has
not been mentioned in their paper. Their idea is to estimate a solution based on the acceleration (the rst time derivative of velocity) at each step without any iteration. Unfortunately, the
given numerical examples were restricted to strain rate sensitive material. The investigation in
the present work follows the same principle as Mole et al.13 but diers from theirs in four
aspects. Firstly, the proposed improving techniques are deduced in the framework of the perturbation technique and can be extended straightforward to higher orders. Secondly, the proposed
computational procedure deals not only with strain rate sensitive materials but also applicable to
work hardening materials. Thirdly, the predicted solution, based on the rst time derivative of
the velocity, is just taken as the trial solution for the rst iteration at each step. This ensures
some accuracy and numerical stability. Finally, on the basis of both strain increment and strain
rate increment, an adaptive step strategy is considered to optimize the computational procedure.
In addition, the proposed improvements were readily inserted into an existing viscoplastic nite
element code.
The paper is arranged in the following way. In Section 2, we will give the mechanical model
used in this study. In Section 3, a brief review on the classical NRM will be given. In Section 4,
we will present the perturbation technique in the framework of large strain viscoplasticity. In Section 5, some techniques will be provided to dene an adaptive step size and to improve the trial
solution for the rst iteration at each step. A new computational procedure will be detailed at the
end of this section. To assess the eciency of the proposed techniques and check their validity, a
number of numerical examples are examined in Section 6 such as necking of a circular bar, plane
strain bending, and axisymmetric hydrostatic bulging processes. The computational results obtained
by the proposed scheme are compared with that of the classical NRM using equal or adaptive
steps, with the numerical solution of the ANM, also with a published computational solution14 as
well as with the available experimental results.15 Finally, in Appendix I, an usual computational
procedure used in the classical NRM will be given to facilitate the comparison with respect to
the proposed method. For the same reason, The ANM will be briey described in Appendix II.
Appendix III gives in detail the description to deduce the rst time derivative formulation in large
strain viscoplasticity.
2. THE MECHANICAL MODEL TO BE SOLVED
This section describes the governing equations of large strain viscoplastic problems and some
notations used in the paper.
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 155–176 (1999)
IMPROVED ITERATIVE METHOD FOR LARGE STRAIN VISCOPLASTIC PROBLEMS
157
2.1. Equilibrium equations and boundary conditions
By neglecting the inertial and volumic forces, the equilibrium of a solid body occupying a
deformed region with boundary @
can be expressed as
div b = 0
in ;
V = Vd
on @
v ;
b · n = td
on @
(1)
where @
v ∪ @
= @
and @
v ∩ @
= ∅. In (1), V and b denote the velocity and stress elds over
. V d represents the prescribed velocity over the kinematic boundary @
v and t d the prescribed
traction over @
with as a load parameter, which depends on time in a process. For the present
work, is taken as time without loss of generality. Finally, n gives the external normal vector to
the surface @
.
2.2. The viscoplastic law
In the present work, the elastic strain is neglected and a viscoplastic law is assumed:
b = b 0 − pI
with b 0 = cD 0
tr D = 0
(2)
in which b denotes the stress tensor, p the hydrostatic pressure, c the material consistency coecient, and D the strain rate tensor which is the symmetric part of the velocity gradient tensor. The
prime 0 and the superscript t imply a deviatoric tensor and a transposed tensor, respectively. The
symbol tr(·) gives tr D = Dii , which denes the volume strain rate. The evolution of c is generally
governed by
r
Z t
2 0 0
n m−1
2
D : D = 0 +
with D =
D ds
(3)
c = 3 k( + ) D
3
t0
Here k, , n and m are material constants, D is the equivalent strain rate, is the cumulated
plastic strain and 0 denes the initial value of the plastic equivalent strain at the beginning of
time step t0 : The symbol (:) implies the contracted product of two tensors D : D = Dij Dij .
As in most numerical codes, the material incompressibility is treated by using a penality method.
In this case, (2)3 can be replaced by
tr D = −
p
g
(4)
with g a very large constant.
To avoid the division by a very small value when the strain rate D is near to zero, equations
(3) are regularized as
r
Z t
2 0 0
n m−1
2
DA ds
(5)
D : D + Á 2 ; = 0 +
c = 3 k( + ) DA ; DA =
3
t0
in which Á is a small numerical parameter. If Á is zero, equations (5) reduce to (3). In most
viscoplastic nite element codes, a tolerance is introduced to prevent the division by a very small
value. These two techniques have not much essential dierence. From now on, equations (4) and
(5) will be employed in place of (2)3 and (3).
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 155–176 (1999)
158
H.-L. CAO AND M. POTIER-FERRY
2.3. Principle of virtual velocities
In practice, a spatial discretization technique such as the nite element method is commonly used
to obtain approximate solutions to the above equations. To do that, it is necessary to reformulate
equations (1) and (2) into a weak form—the principle of virtual velocities. They are thus rewritten
as
Z
Z
(b 0 : D 0 − p tr D) d
=
t d · V d
@
Z p
p d
= 0
(6)
tr D +
g
for all weight functions p and V which satisfy the homogenous boundary conditions on @
v . In
(6), represents boundary surfaces. The unknowns are the velocity V and the pressure p elds.
The integration is carried out on the deformed conguration estimated at the end of the time
step. At a given time step, the Newton–Raphson method uses the linearized form of (6) to iterate
until its residual norm is admissible at the end of the time step.
In order to compute the time derivatives of the velocity and the pressure at each time step, the
above equations should be expressed on a known reference conguration (the conguration at the
beginning of the time step for instance), designated as 0 with boundaries @
0 , in the form
Z
Z
Z
p
p d
= 0
(7)
T : L d
=
t d · V d ;
J tr D +
g
0
@
0
0
with T being the rst Piola–Kirchho stress tensor, and t d being the traction vector,
T = J bF−t
or
T = J (b 0 − pI)F−t ;
t d = J bF−t N
(8)
Here, F is the deformation gradient tensor, L the velocity gradient tensor with respect to the
reference conguration, J the determinant of F and N indicates the external normal vector to the
surface @
0 . In the cases where the stress is prescribed b = b d instead of t d on the surface, t d
varies as a function of F (in an hydrostatic bulging process, for instance).
As in most viscoplastic nite element codes, the process here is described using an updated
Lagrangian scheme. The co-ordinate vector x of the deformed conguration is updated at the end
of each time step as
Z t
V ds
(9)
x = x0 +
t0
where x0 is the co-ordinate vector at the beginning of time step t0 : This new updated conguration
is taken as the reference conguration for the subsequent time step.
3. REVIEW ON THE CLASSICAL NEWTON–RAPHSON METHOD
The Newton–Raphson method is well known. A general review on the method can be found in
Reference 16. Here, we make precise only the NRM for large strain viscoplastic problems using the
updated Lagrangian scheme, which allows us to understand more clearly the proposed improving
techniques and the facilities of its insertion into an existing viscoplastic nite element code. An
usual computational procedure used in the classical NRM is given in Appendix I.
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 155–176 (1999)
IMPROVED ITERATIVE METHOD FOR LARGE STRAIN VISCOPLASTIC PROBLEMS
159
Generally, the classical NRM uses the linearized form of (6) at each iteration.1 This form can
be expressed as
Z
[(M c : dD 0 ) : D 0 − dp tr D] d
= R ext − Rint
Z dp
p d
= Rp
tr dD +
(10)
g
with
nDA t
2c
D0 ⊗ D0
(m − 1·0) +
M = cI I +
( + )
3DA2
Z
1
dD = (dl + dlt ); Rint =
(b 0 : D 0 − p tr D) d
2
Z
Z p
p d
tr D +
t d · V d ; Rp = −
R ext =
g
@
c
where M c is the eective viscoplastic modulus tensor (symmetric). dl, dD 0 , and dp are the
corrections between two successive iterations to the eld velocity gradient l, to the eld D 0 and to
the eld p, respectively. t represents the time increment. The integration and its integral function
in (10) are evaluated by using the newly obtained solution of the preceding iteration. Finally, the
0
.
symbols and ⊗ imply: (I I)ijkl = ik jl and (D 0 ⊗ D 0 )ijkl = Dij0 Dkl
By using the nite element discretization, we obtain a set of equations from (10), schematically
as
[K c ]{dU } = {R}
(11)
where {dU } includes the correction of velocity {dV } and the correction of pressure {dp}. [K c ]
stands for the tangent stiness matrix (symmetric) and {R} denotes the residual force array. This
is the system to be solved at each iteration in the classical NRM.
4. PERTURBATION TECHNIQUE
Let us represent the weak form (7) by the symbol (U; ) to simplify the discussion. Here U
stands for both the unknown velocity and the unknown pressure, and denotes the load parameter.
Since U and are functions of time, (U; ) can also be considered as a function of time.
Supposing that at time t0 we have U(0) and (0) which satisfy the weak form (7):
(U(0) ; (0) ) = 0
(12)
we can take the Taylor expansion of U, and (U; ) in the vicinity of t0 for any time varying
from t0 to t:
U() = U(0) + U; t |t0 ( − t0 ) +
() = (0) + ; t |t0 ( − t0 ) +
(U(); ()) = 0
Copyright ? 1999 John Wiley & Sons, Ltd.
1
2! U; tt |t0 (
1
2! ; tt |t0 (
− t0 ) 2 + · · ·
− t0 ) 2 + · · ·
(13)
Int. J. Numer. Meth. Engng. 44, 155–176 (1999)
160
with
H.-L. CAO AND M. POTIER-FERRY
(U(); ()) = (U(0) ;
(0) ) + ; t |t0 ( − t0 ) +
1
2! ; tt |t0 (
− t0 ) 2 + · · ·
where the symbol (•); t |t0 denotes the rst time derivative taken at time t0 and (•); tt |t0 is the
second time derivative. From now on, we will use the symbol (•)(1) to designate (•); t |t0 and (•)(2)
for (•); tt |t0 : As relation (13)3 should be valid whatever is the value , we obtain, by taking into
account (12):
(1) = 0;
(2) = 0; : : :
(14)
Futhermore, (14) can be reformulated as
(1) = ; U |t0 U(1) + ; |t0 (1) = 0
(2) = ; U |t0 U(2) + ; |t0 (2) + ; UU |t0 U(1) U(1)
+ 2; U |t0 U(1) (1) + ; |t0 (1) (1) = 0
:::
(15)
From (15)2 , we observe that the second time derivatives U(2) and (2) have the same operator as
the rst time derivatives U(1) and (1) : The remaining terms involve only terms dependent on lower
time derivatives (the rst time derivatives in (15)2 , for instance). For the higher time derivatives
U(3) , (3) ; : : : ; the same characteristic can be observed. It is this feature which is exploited in the
ANM12; 18 because it needs only one matrix inversion to estimate any higher time derivatives.
Neglecting the higher time derivatives than the one in (13) and taking into account (15), we
have
U() = U(0) + U(1) ( − t0 ); () = (0) + (1) ( − t0 )
; U |t0 U(1) + ; |t0 (1) = 0
(16)
Of course, we can keep the higher derivatives until n order in the computation. In this case,
the method is called as asymptotic numerical method10 or incremental perturbation method.9 The
interested reader can nd more details for large strain viscoplastic problems in Reference 12.
A short description is given in Appendix II. In the present work, we follow the same idea as
Mole et al.13 by using only the rst time derivative to improve the classical Newton–Raphson
method.
Introducing the rst time derivative of (7) at time t0 into (16)3 we obtain
Z
0∗
ext
[(M : D(1)
) : D 0 − p(1) tr D] d
= R(1)
− Rv ∗
0
Z 1
∗
(17)
tr D(1) + p(1) p d
= Rp ∗
g
0
with
D∗
(1)ij
ext
R(1)
@V(1) j
1 @V(1)i
2c(m − 1) 0
=
+
D ⊗ D0
; M = cI I +
2 @xj
@xi
3DA2
Z
d
d
=
((1) t d + t(1)
) · V d ; t(1)
= 0; if t d is prescribed
@
0
d
= b d (tr (l)I − lt )N;
t(1)
Copyright ? 1999 John Wiley & Sons, Ltd.
if b = b d is prescribed
Int. J. Numer. Meth. Engng. 44, 155–176 (1999)
IMPROVED ITERATIVE METHOD FOR LARGE STRAIN VISCOPLASTIC PROBLEMS
R
v∗
Z
=
0
161
{b[tr(l)I − lt ] : l + [(Rc)D 0 − M : RD ] : D 0 } d
(Rc) = 23 kn( + 0 ) (n−1) DA ; RD = 12 (lc + lct ); lc = ll
Z p
p∗
t
R = −
tr(l)I + l : [tr(l)I − l ] p d
g
0
in which the integration and all the integral functions are estimated on the conguration at time t0 :
From this set of equations we obtain a linear system of unknown elds V(1) and p(1) once V(0) and
p(0) are known. Some details about the derivation of equations (17) are presented in Appendix III.
Comparing system (10) with (17), it can be noted that they have the same structure with their
respective unknowns (dV, dp in (10) and V(1) , p(1) in (17)). If system (10) is constructed by
using the last convergent solution at time t0 , the only dierence between the operators of their
unknowns concerns the eective viscoplastic modulus M c and M: This dierence comes from the
work hardening of the material. If n = 0, M is identical to M c , which has been dealt with by
Mole et al.13 If t = 0, M is also identical to M c : This last point is taken into account in the
ANM to nd out the solution at time t0 and the corresponding time derivatives. In general, M c
and M are dierent. However, M can be approximated by M c for most applications. In fact, the
term involving work hardening in M c can be rewritten as
−1
+ 0
nDA t
=n 1 +
+
DA t
The value of this term is small with regard to |m − 1| for most metallic materials where n is
around 0·2 and m is less than 0·5. Particularly, equivalent strain increment is usually controlled
under 5 per cent in numerical simulations. If the term involving work hardening is neglected, it
can decrease the convergent rate in the NRM. On the contrary, the substitution of M with M c in
(17) gives a good approximation of V(1) and p(1) without constructing and inversing the stiness
matrix once more. This will be shown by the numerical examples.
After nite element discretization, system (17) becomes
[K]{U(1) } = {R∗ }
(18)
where {U(1) } involves the rst time derivatives of velocity {V(1) } and pressure {p(1) } at time t0 :
[K] stands for the stiness matrix, which will become [K c ] if M in (17) is replaced by M c :
ext
, Rv ∗ and Rp ∗ : This residual
Finally, {R∗ } denotes the residual force array coming from R(1)
force can be calculated in a similar way as {R} in (11).
5. IMPROVED NEWTON–RAPHSON METHOD
Three points will be discussed in this section, the adaptive step size, the trial solution for the rst
iteration at each step and the consequent computational procedure.
5.1. Adaptive step size
As indicated by Chung and Wagoner,4 an adaptive step size is generally better than an optimal
equal one in terms of computational eort. This is why the adaptive step size is implemented
in some commercial codes (ABAQUS, etc.). The problem for this strategy concerns essentially
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 155–176 (1999)
162
H.-L. CAO AND M. POTIER-FERRY
their chosen criteria by which the step is controlled. Normally, the usual criteria are the maximum
admissible strain increment and the maximum incremental rotation. These limitations allow to
avoid divergence if the non-linearities originate from the variation of geometry or that of the
equivalent strain. However, when the non-linearity comes from the variation of strain rates, the
previous limitations will lose their eciency.
As demonstrated in Section 4, the rst time derivatives of velocity and pressure may be computed
easily once the velocity and pressure elds are known. We can use these data to control the step
size apart from the classical limitations on strain increment and incremental rotation. Based on
both the limitations of velocity variation and that of the velocity, the adaptive step size becomes
more ecient even if the strain rate varies signicantly.
In the realized nite element code, the following criteria are chosen: limitation on the increment
of velocity gradient |l| and limitation on the relative variation of velocity gradient with respect
to its velocity gradient |l(1) |=|l|. Of course, other limitations can also be implemented into the
code, limitations based on strain increment and the variation of strain rate, for example. All these
will not inuence the adaptive step size considerably.
Compared with the method proposed by Mole et al.,13 the introduction of the present adaptive
step size is more satisfactory, which can be observed from the numerical examples presented in
Section 6.
5.2. Trial solution for the rst iteration at each step
If V(0) and p(0) at time t0 are obtained, we can use the last inversed stiness matrix [K c ]−1
instead of [K]−1 to compute V(1) and p(1) for system (18). The next time step size t can be
chosen according to the criteria discussed above. After the choice of step size, the trial solution
for the rst iteration at time t0 + t can be predicted as follows:
V = V(0) + V(1) t;
p = p(0) + p(1) t
(19)
As a consequence, the trial conguration x and the trial equivalent strain at the time t0 + t
can be estimated as
x = x0 + tV(0) + 12 t 2 V(1) ;
= 0 + tDA(0) + 12 t 2 DA(1)
(20)
To keep the consistancy with (20), the conguration x and equivalent strain at the end of
each iteration can be updated as
x = x0 + 12 t(V(0) + V);
= 0 + 12 t(DA(0) + DA )
(21)
where V is the newly obtained solution at the time t0 + t, and DA is the equivalent strain rate
corresponding to V and x:
Furthermore, the viscoplastic modulus M c are modied to take into account the updated scheme
(21), which is rewritten as
nDA t
2c
D0 ⊗ D0
(m
−
1)
+
(22)
M c = cI I +
2( + )
3DA2
It is this modulus which is used to construct the stiness matrix [K c ] in the improved Newton–
Raphson method.
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 155–176 (1999)
IMPROVED ITERATIVE METHOD FOR LARGE STRAIN VISCOPLASTIC PROBLEMS
163
5.3. Computational procedure
With the consideration of the two techniques above, the new computational procedure is conceived as follows for the cases where the boundary conditions are supposed to be unchanged. For
the rst time step, the classical NRM is employed to nd the solution. Some measures proposed
by Wang1 and others2; 4; 5 can be used to increase the convergence rate for this rst time step. In
this study, we focus on the improved techniques only for the numerical solutions after the rst
step.
(0) Given a prescribed time step t, the initial time t0 , the initial reference conguration, the
initial equivalent plastic strain eld and the trial solution, nd out the solution {U } at
time t (1) = t0 + t by using the classical NRM.
Loop of the deformation process for l = 2 to lmax
(1) Construct the residual force array {R}∗ of system (18).
(2) Compute {U(1) } by using the last inversed stiness matrix [K c ]−1 at the time step t (l−1) :
(3) Choose the time step size t a according to some criteria (the two criteria described above,
for example).
(4) With the chosen step size, estimate the trial solution by (19), the trial conguration and
the trial equivalent strain by (20) at time t (l−1) + t a :
Loop of the equilibrium at each step for n = 1 to n max
(5) Estimate the residual forces {R} of (11).
(6) Check if the norm of the residual forces {R} is less than the admissible tolerance. If yes,
pass to (12); if not, continue.
(7) Form the tangent stiness matrix [K c ] of (11) with consideration of the modied viscoplastic modulus M c in (22).
(8) Solve system (11) for the correction elds {dU }.
(9) Compute the new solution
{U } (n) = {U } (n−1) + ÿ{dU }
(23)
where ÿ is a relaxation coecient varying from 0 to 1.
(10) Update the reference conguration and the equivalent strain eld according to (21) with
the determined time step size t a :
(11) end of the equilibrium loop.
(12) Update the time:
t (l) = t (l−1) + t a
(13) Verify if the time is superior or equal to the nal time of the process. If yes, pass to
(15); If not, continue.
(14) End of the deformation process loop.
(15) End of the process.
In the above procedure lmax denotes the maximum number of time steps and n max is the maximum
number of iterations.
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 155–176 (1999)
164
H.-L. CAO AND M. POTIER-FERRY
As compared to the classical NRM with adaptive step size presented in Section 3, the proposed
computational procedure has three major features. Firstly, the adaptive step size is determined
according to the criteria based not only on the current velocity and pressure but also on their
rst time derivatives. Secondly, the trial solution for the rst iteration at each step is calculated
according to (19) in place of the solution at the preceding time step. Finally, the eective viscoplastic modulus is modied to take into account the updated scheme expressed in (22). Between
the proposed computational procedure and that described in Section 3, the dierent points are
concerned with the steps about the choice of the trial solution and that of the step size. The rest
of the procedure is the same. This characteristic allows one to insert easily the present scheme in
an existing viscoplastic nite element code.
The scheme described by Mole et al.,13 which is restricted to strain rate sensitive material
without work hardening, adopts the equal step sizes and considers the trial solution by (19) as
the real solution for the time t0 + t: Their scheme needs only one iteration per step after the
rst step of the deformation process. To ensure a good solution, however, the step size should
be suciently small. If not, the obtained solution may become completely erroneous. Comparing
their scheme with the proposed one, the latter takes into account the work hardening eect and
adopts an improved adaptive step strategy. In addition, the former can be considered as a particular case of the latter. It is sucient to relax all limitations on the step size and put n max equal
to 1.
In comparison with the ANM (a short description is given in Appendix II), the present scheme
needs only the rst time derivatives while the ANM requires generally until 10th or higher time
derivatives. The proposed scheme thus reduces the memory storage. Since the present scheme
needs only the rst time derivatives, it is not necessary to regularize strongly the used formulations
(presenting a sharp variation) as recommended in Reference 17 for the ANM. Moreover, it is easier
to implement the present scheme into an existing code than the ANM.
6. NUMERICAL EXAMPLES
In this section, several numerical examples are presented to assess the performance of the improving
techniques proposed above. For each example, ve solution procedures are employed: the classical
NRM with equal time steps, the classical NRM with adaptive time steps, the one-iteration scheme
using equal time steps (it becomes the scheme proposed by the Mole et al.13 in the case of strain
rate sensitive material), the ANM until the 15th time derivative, and the present proposed scheme.
The numerical results obtained by these ve dierent procedures are compared in terms of their
computational requirements (the number of time steps and that of matrix inversions) and their
nal kinematic variables (velocity and equivalent strain). Three kinds of examples are given here.
The rst is concerned with the necking of a circular bar, which has been considered by Mole
et al.13 The second is plane strain bending under hydrostatic pressure. The nal example deals
with an axisymmetric hydrostatic bulging process, simulated by Kim and Yang14 and realized
experimentally by Iseki et al.15
In the simulations, a 4-node isoparametric element with bilinear velocity interpolation and constant pressure is used. A relative residual norm of 10−2 is generally adopted. For the simulations except that using the ANM, ÿ is taken equal to 0·75 (described in (23)) to facilitate
the comparison. As information, all the simulations have been carried out on a HP9000=K200
server.
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 155–176 (1999)
IMPROVED ITERATIVE METHOD FOR LARGE STRAIN VISCOPLASTIC PROBLEMS
165
Figure 1. (a) The geometry of the circular bar (mm); (b) the initial mesh; (c) the deformed mesh corresponding to the
nal vertical displacement of the upper head 7·1 mm
6.1. Necking of a circular bar
In this example, the geometry, the initial nite element mesh, and the deformed mesh for a half
of a cross-section of the specimen are shown in Figure 1. The specimen is assumed to be pulled
by an uniform velocity of 1000·0 mm=s at the upper head and is xed in the radial direction. One
hundred elements are used in the computation.
The rst simulation is carried out with k = 1·0 Pa s−0·1 , m = 0·1, and n = 0 as the material
constants, which have been considered by Mole et al.13 For Á, it is taken equal to 10−5 =s.
The rst step size in the simulation is set to 10−4 s. This step is solved by using the classical
NRM. After the rst step, we apply the ve computational procedures for the rest of the deformation process. In Table I are registered the equivalent strain at the centre of the specimen, the radial
velocity of the external node of the central section, and the computational requirements except the
rst step. All the results correspond to the nal vertical displacement of the upper head 7·1 mm.
The reference results are obtained by using the classical NRM with equal time steps t = 10−5 s.
From Table I, we observe that the NRM with adaptive steps needs less matrix inversions and
time steps than the NRM with equal steps to reach the same accuracy in terms of : This remark
conrms the observation given by Chung and Wagoner.4 The most ecient case of Mole’s scheme
and that of the proposed scheme require generally less matrix inversions (2 times less) than the
optimal NRM. Moreover, the present scheme obtains a very good solution (inferior to 1·0 per cent
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 155–176 (1999)
166
H.-L. CAO AND M. POTIER-FERRY
Table I. The equivalent strain at the centre, the radial velocity of the external
node at the central section, and the computational requirements except the rst
step (time steps, matrix inversions simplied as MI), obtained when the tensile process is stopped at 7·1 mm vertical displacement of the upper head. Here
Vref = −985·7 mm=s and ref = 1·4180
Steps
MI
V=Vref
=ref
NRM
equal step
t (10−3 s)
0·1
0·25
0·5
71
28
14
138
56
28
1·0149
1·0016
1·0058
1·0045
1·0133
1·0333
NRM adaptive
step |l|max and
tmax = 0·75 (10−3 s)
0·05
0·1
43
21
55
37
1·0321
1·0610
0·9922
0·9925
Mole’s scheme
t (10−3 s)
0·1
0·25
0·5
71
29
15
71
29
15
1·0000
1·0007
1·0021
0·9950
0·9941
0·9937
Present scheme
|l|max and
(|l(1) |=|l|)max
0·01
0·1
0·1
0·1
0·15
0·1
223
223
1·0000
0·9950
23
23
1·0001
0·9945
15
15
1·0005
0·9937
10
10
0·9991
0·9946
ANM (time derivatives)
15
dierence) with respect to the reference. Finally, the ANM, as expected, gives an accurate solution
with less matrix inversions and time steps with respect to all other schemes.
Another simulation is carried out with a material having both, work hardening eect and strain
rate sensitivity. The material constants are: k = 1·0 Pa s−0·01 , m = 0·01, n = 0·20, and = 0·001.
In Table II are registered the same kind of data as in Table I. Some similar remarks can be
observed as in the preceding example. On the contrary, The one-iteration scheme, whenever the
matrix [K c ] or [K] in (18) is used, gives similar results. It generally requires less matrix inversions
and time steps than the optimal NRM. However, the obtained velocity is completely erroneous
with respect to the reference in the cases of t = 0·125 × 10−3 s and t = 0·25 × 10−3 s. This
gives us some sign of numerical instability of the scheme when the time step is a little bigger.
6.2. Plane strain bending
In this example, a metallic blank of 24 mm in width and 0·314 mm in thickness is subjected to
a time-dependent hydrostatic pressure. The material constants are taken as: k = 180·0 MPa s−0·2 ,
m = 0·0, n = 0·29, = 0·000769, and Á = 10−5 =s.
Because of the symmetry, only a half of the blank is used in the computation. Concerning the
boundary conditions, the blank is clamped at its border and xed at the centre in the horizontal
direction. On the upper surface, a hydrostatic pressure is applied with a rate of 10−3 MPa=s. The
blank is discretized by means of two elements across the thickness and 20 elements in horizontal
direction. In Figure 2 the deformed mesh is displayed corresponding to t = 1353·3 s after the start
of the deformation process.
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 155–176 (1999)
167
IMPROVED ITERATIVE METHOD FOR LARGE STRAIN VISCOPLASTIC PROBLEMS
Table II. The equivalent strain at the centre, the radial velocity of the external node at the central section, and the computational requirements except the
rst step (time steps, matrix inversions simplied as MI), obtained when the tensile process is stopped at 7·1 mm vertical displacement of the upper head. Here
Vref = −1134·1 mm=s and ref = 1·3480
Steps
MI
V=Vref
=ref
NRM
equal step
t (10−3 s)
0·1
0·125
0·25
70
56
100
1·0411
104
1·0503
No convergence
0·9874
0·9887
NRM adaptive
step |l|max and
tmax = 0·25 (10−3 s)
0·025
0·04
0·05
79
49
113
1·0259
84
1·0361
No convergence
0·9874
0·9867
One-iteration
scheme with [K c ]
t (10−3 s)
0·1
0·125
0·25
0·35
70
56
28
70
0·9986
56
0·7555
28
0·4646
Erroneous results
0·9941
0·9941
0·9503
One-iteration
scheme with [K]
t (10−3 s)
0·125
0·25
0·35
56
28
56
0·7417
28
0·3154
Erroneous results
0·9941
0·9573
Present scheme
|l|max and
(|l(1) |=|l|)max
0·01
0·1
0·05
0·1
0·075
0·1
205
206
0·9987
0·9941
42
46
0·9993
0·9934
ANM (time derivatives)
15
No convergence
39
40
0·9989
0·9941
Figure 2. The deformed mesh corresponding to time t = 1353·3 s in the plane strain bending process
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 155–176 (1999)
168
H.-L. CAO AND M. POTIER-FERRY
Figure 3. The variation of vertical velocity at the centre as a function of time
Table III. The equivalent strain at the centre and the computational requirements (time
steps, matrix inversions simplied as MI), obtained when the plane strain bending
process is stopped at t = 1353·3 s. Here ref = 0·2562
Steps
NRM
equal step
t (s)
1·0
5·0
10·0
12·5
MI
=ref
1354
272
137
3865
0·2564
1230
0·2569
852
0·2568
No convergence
1·0008
1·0027
1·0023
NRM adaptive
step |l|max and
tmax = 25·0 s
0·01
0·02
0·03
0·04
225
112
74
1097
0·2566
656
0·2569
796
0·2571
No convergence
1·0016
1·0027
1·0035
One-iteration
scheme with [K c ]
t (s)
1·25
2·5
5·0
7·5
1059
530
266
1161
0·2559
632
0·2560
368
0·2561
Erroneous results
0·9988
0·9992
0·9996
One-iteration
scheme with [K]
t (s)
0·5
1·0
1·25
2646
2768
0·2559
Erroneous results
Erroneous results
0·9987
Present scheme
|l|max and
(|l(1) |=|l|)max
0·01
0·1
0·02
0·1
ANM (time derivatives)
15
275
340
0·2561
0·9996
162
308
0·2562
1·0000
243
311
0·2561
0·9996
For the rst time step, t is set to 2·5 s to be able to attain rapidly the convergence (three
iterations). With this rst time step size, the one-iteration scheme has not obtained the right solution
for t superior to 0·5 s. In order to be able to increase the time step size, a bigger time step size
t = 30·0 s is used at the rst step for the simulations using the one-iteration scheme. The reason
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 155–176 (1999)
IMPROVED ITERATIVE METHOD FOR LARGE STRAIN VISCOPLASTIC PROBLEMS
169
Figure 4. The deformed mesh corresponding to time t = 2090·0 s in the axisymmetric hydrostatic bulging process
for this problem can be explained by the sharp time variation of the central vertical velocity at
the beginning of the process shown in Figure 3. From this gure, it is observed that it is very
dicult to estimate the real solution without any iteration near to t = 15·0 s if the step is too big.
Table III shows the equivalent strain at the centre and the computational requirements, corresponding to the time t = 1353·3 s. It is recalled that the simulations using the one-iteration scheme
adopts a time increment of 30·0 s for the rst step instead of 2·5 s as used in other simulations. The
reference results are obtained by using the classical NRM with equal steps of t = 5·0 × 10·0−2 s.
From this table, we observe the same phenomena as in the preceding examples. We emphasize here
that the one-iteration scheme using [K c ] shows much better performance than that of using [K]:
Moreover, the present scheme illustrates well its eciency with respect to others even compared
to the ANM.
6.3. Axisymmetric hydrostatic bulging
For this example, a metallic blank of 24 mm in diameter and 0·314 mm in thickness is subjected
to a time-dependent hydrostatic pressure. The material constants are the same as in the preceding
example of the plane strain bending. This example has been simulated by Kim and Yang14 using
the classical NRM and experimentally realized by Iseki et al.15
In the simulations, the similar boundary conditions and the same mesh topology are applied
as for the plane strain bending. In Figure 4 the deformed mesh is displayed corresponding to
t = 2090 s after the start of the deformation process. For the rst time step, t is set to 2·5 s.
Table IV shows the equivalent strain at the centre and the computational requirements, corresponding to t = 2090·0 s. The reference results are obtained by using the classical NRM with equal
steps of t = 0·1 s. We note here that the most ecient case of the proposed scheme requires 2
times less matrix inversions and almost 3 times less time steps than the optimal case of the oneiteration scheme. With respect to the best of the NRM, the gain is 1·5 times in terms of matrix
inversions. Regarding the ANM, it still stays the most satisfactory.
Figure 5 depicts the variation of the time steps during the deformation process, obtained by the
proposed scheme with |l|max = 0·01 and the maximum of (|l(1) |=|l|) equals to 0·1. It is remarked
that the step varies signicantly during the process. In Figure 6 the variation of polar thickness
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 155–176 (1999)
170
H.-L. CAO AND M. POTIER-FERRY
Table IV. The equivalent strain at the centre and the computational requirements (time
steps, matrix inversions simplied as MI), obtained when the axisymmetric hydrostatic
bulging process is stopped at t = 2090·0 s. Here ref = 0·5278
Steps
NRM
equal step
t (s)
6·25
12·5
25·0
37·5
336
169
85
NRM adaptive
step |l|max and
tmax = 25·0 s
0·01
0·02
0·03
163
90
85
One-iteration
scheme with [K c ]
t (s)
1·0
2·5
5·0
7·5
One-iteration
scheme with [K]
t (s)
5·0
6·25
7·5
Present scheme
|l|max and
(|l(1) |=|l|)max
0·01
0·1
0·015
0·1
0·02
0·1
ANM (time derivatives)
15
MI
1313
0·5258
776
0·5235
734
0·5161
No convergence
729
566
788
=ref
0·9962
0·9919
0·9778
0·5237
0·5181
0·5159
0·9922
0·9816
0·9775
2092
838
420
2093
0·5279
839
0·5279
421
0·5278
Erroneous results
1·0081
1·0002
1·0000
420
421
0·5277
Erroneous results
Erroneous results
0·9998
206
247
0·5278
1·0000
157
224
0·5278
1·0000
137
421
0·5265
0·9975
96
97
0·5278
1·0000
Figure 5. The variation of time step size as a function of time
strain versus hydrostatic pressure is shown, corresponding to, respectively, the present analysis,
the viscoplastic nite element analysis using the NRM14 and the experimental data published in
Reference 15. It is seen that the present numerical results agree well with the experimental data15
and the published computational results.14
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 155–176 (1999)
171
IMPROVED ITERATIVE METHOD FOR LARGE STRAIN VISCOPLASTIC PROBLEMS
Figure 6. The variation of polar thickness strain as a function of the hydrostatic pressure
Table V. The equivalent strain at the centre and the computational requirements (time steps, matrix
inversions simplied as MI and CPU time), obtained when the axisymmetric hydrostatic bulging process is stopped at t = 2000·0 s. Here the blank is discretized by two elements across the thickness
and 100 elements in the radial direction
Steps
NRM
equal step
t (s)
2·0
3·0
6·25
NRM adaptive
step |l|max and tmax = 25·0 s
0·0025
0·01
One-iteration
scheme with [K c ]
t (s)
Present scheme
|l|max and (|l(1) |=|l|)max
ANM (time derivatives)
MI
CPU time (h)
1000
3055
0·4191
28·19
No convergence at t = 1931·5 s
No convergence at t = 1558·75 s
695
2288
0·4196
21·56
No convergence at t = 1617·8 s
1·0
1·5
2·5
1999
2001
0·4193
22·79
Erroneous results from t = 1990·0 s
Erroneous results from t = 1895·0 s
0·01
0·1
598
617
0·4193
6·81
365
652
0·4193
8·42
15
In order to check if the proposed scheme still keeps its eciency for problems where the number
of degrees of freedom is relatively important, the preceding axisymmetric hydrostatic bulging
process is simulated with a ner mesh (two elements across the thickness and 100 elements in
the radial direction). In Table V the computational requirements (matrix inversions, time steps and
CPU time) are registered, corresponding to t = 2000·0 s. It is observed clearly that the proposed
scheme shows its superiority with respect to any other scheme in terms of matrix inversions, time
steps and also CPU time. In addition, the one-iteration scheme also gives good results with respect
to the classical NRM using equal steps, 1·5 times less in terms of matrix inversions and 1·24 times
less in terms of CPU time. This result does not agree well with the observation given by Mole
et al.13 where Mole’s scheme required 3 times more CPU time than the classical NRM.
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 155–176 (1999)
172
H.-L. CAO AND M. POTIER-FERRY
7. CONCLUSIONS
In the present work, an improved Newton–Raphson scheme is presented for numerical analysis
of large strain viscoplastic problems. The improvements, based on the perturbation technique, are
concerned with the adaptive step strategy and with the rst trial solution at each step.
From all the performed numerical simulations, some remarks can be made here. (1) The proposed
scheme generally obtains an accurate solution with less computational requirements with respect
to the classical NRM and also with respect to the one-iteration scheme (Mole’s scheme in the
cases of strain rate sensitive materials). For a typical axisymmetric hydrostatic bulging process,
the proposed scheme requires 5 times less matrix inversions, 1·7 times less time steps, and 4 times
less CPU time with respect to the classical NRM using equal steps; (2) it is always easier for
users to prescribe some limitations with physical signication, which are independent of the treated
problems, than to look for the relative optimal time step used in the NRM by doing numerous
tests; (3) although the ANM seems in general better than the proposed scheme in this work, it
needs also more memory storage and more modications to be inserted in an existing viscoplastic
nite element code.
It is useful to note here that the proposed improvements are given in the framework of large
strain viscoplastic problems but the basic idea is amenable to other non-linear problems. Additionally, as a restriction of the proposed improving techniques, a direct solver should be used to
solve the linear system formed by nite element discretization. If an iterative solver is used, the
proposed scheme will lose its advantages.
APPENDIX I. AN USUAL COMPUTATIONAL PROCEDURE USED
IN THE CLASSICAL NRM
In a viscoplastic nite element code, the computational procedure can be summarized as follows
if the boundary conditions are supposed to be unchanged.
Known parameters are: the prescribed time step t, the initial reference conguration and the
initial equivalent plastic strain eld.
Loop of the deformation process for l = 1 to lmax
(0) The trial solution is given and t (0) = t0
Loop of the equilibrium at each step for n = 1 to n max
(1) Estimate the residual forces.
(2) Check if the norm of the residual forces {R} in (11) is less than the admissible tolerance.
If yes, pass to (9); if not, continue.
(3) Construct the tangent stiness matrix [K c ].
(4) Solve the equation system (11) to obtain the correction elds {dU }.
(5) Compute the new solution:
{U } (n) = {U } (n−1) + ÿ{dU }
where ÿ is a relaxation coecient varying from 0 to 1.
(6) Choose the time step size t c according to the used criteria (the maximum admissible
equivalent strain increment, for example).
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 155–176 (1999)
IMPROVED ITERATIVE METHOD FOR LARGE STRAIN VISCOPLASTIC PROBLEMS
173
(7) Update the reference conguration and the equivalent strain eld:
{X } (n) = {X } (n−1) + t c {V } n
{} (n) = {} (n−1) + t c {D} n
(8) End of the equilibrium loop.
(9) Update the time: t (l) = t (l−1) + t c :
(10) Verify if the time is superior or equal to the nal time of the process. If yes, pass to
(12); if not, continue.
(11) End of the deformation process loop.
(12) End of the process.
For step (0), the trial solution can be taken as the solution of the previous step except for the
rst step, where the trial solution may be the solution of the linear Newtonian problem (m = 1·0
and n = 0·0), or that of an elastic problem. In step (5), the value ÿ has an important inuence
on the convergence rate. Wang1 has proposed some variation of ÿ for several typical forming
processes. If step (6) is removed from the above procedure, we obtain the usual procedure used
in the NRM with equal steps.
APPENDIX II. AN ASYMPTOTIC NUMERICAL METHOD FOR
LARGE STRAIN VISCOPLASTIC PROBLEMS
The ANM, based on the perturbation technique presented in Section 4, consists in expanding
the unknowns (velocity V and pressure p) into power series (13) and in solving the obtained
linear equations for each time derivative in a recursive manner. As all the linear systems at a
time step t0 have the same stiness matrix, only one matrix inversion is necessary to estimate
each time derivative. With the time derivatives available, it is trivial to compute the velocity and
pressure for any time according to their respective power series. As the power series have their
own convergence radius, the estimated solution can be accepted when the time increment t from
t0 is inside the radius. A bisection method18 is used to determine this convergence radius or the
time interval during which the relative residual norm and the volumic strain rate are inferior to the
prescribed tolerances. Since the time step is estimated posteriorly, there is no divergence problem
as it is often encountered in the NRM. With the valid time increment t, it is easy to compute the
velocity and pressure as well as other variables at time t0 + t: Following the updated Lagrangian
scheme, the conguration and any used variables at time t0 + t are considered as a new starting
point on the equilibrium path. As the solution estimated by the power series at time t0 + t is not
always the ‘just’ solution after the updated conguration and the change of the starting point for
the expansion, an iterative Newton–Raphson loop is introduced to improve the solution but without
any time increment. After this, the same procedure is repeated to compute the time derivatives at
time t0 + t and so on.
The computational procedure for the ANM is briey described as follows:
(0) Given a prescribed time step t, the initial time t0 , the initial reference conguration,
the initial equivalent plastic strain eld and the trial solution, nd out the solution {U }
at time t (1) = t0 + t by using the classical NRM.
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 155–176 (1999)
174
H.-L. CAO AND M. POTIER-FERRY
Loop of the deformation process for l = 2 to lmax
Loop of the adjustment at time t (l−1) for n = 1 to n max without time increment
(1) Estimate the residual forces {R} of (11) using the last obtained solution.
(2) Check if the norm of the residual forces {R} is less than the admissible tolerance. If yes,
pass to (7); if not, continue. Here, the criteria are stringent in terms of relative residual
norm and volumic strain rate.
(3) Form the tangent stiness matrix [K] of system (18).
(4) Solve system (11) for the correction elds {dU }.
(5) Compute the new solution
{U } (n) = {U } (n−1) + ÿ{dU }
where ÿ is a relaxation coecient varying from 0 to 1.
(6) End of the adjustment loop.
Loop of the time derivatives for r = 1 to rmax
(7) Form the residual force array such as {R}∗ of system (18), but corresponding to the time
derivative at order r:
(8) Compute {U(r) } by using the last inversed stiness matrix [K]−1 .
(9) End of the time derivative loop.
(10) Determine the time interval t a , during which the criteria in terms of relative residual
norm and volumic strain rate are satised by using a bisection method.18 The criteria are
less stringent than that in (2).
(11) Estimate the solution of {p} and {V } according to their respective power series at time
t (l−1) + t a : Other variables can be evaluated according to the corresponding relations
regarding to {p} and {V }:
(12) Update the time:
t (l) = t (l−1) + t a
(13) Update the reference conguration and the corresponding variables at time t (l) .
(14) Verify if the time is superior or equal to the nal time of the process. If yes, pass to
(16); If not, continue.
(15) End of the deformation process loop.
(16) End of the process.
APPENDIX III. THE SET OF EQUATIONS FOR
THE FIRST TIME DERIVATIVE FIELDS V(1) AND p(1)
Equations (17) with unknown elds V(1) and p(1) can be deduced by putting equal to zero the
rst time derivative of equations (7) at time t0 , which comes from the condition described by
equation (16)3 .
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 155–176 (1999)
IMPROVED ITERATIVE METHOD FOR LARGE STRAIN VISCOPLASTIC PROBLEMS
175
From (7), we can express readly their rst time derivatives at time t0 and identify them into
zero as follows:
Z
Z
d
T(1) : L d
=
( (1) t d + t(1)
) · V d
Z
0
@
0
p(1)
p
+ J tr D(1) +
p d
= 0
J(1) tr D +
g
g
0
(24)
where the symbol (•)(1) denotes the rst time derivative of the corresponding variable at time t0 .
d
and the known elds
The main task is to establish the relations between T(1) , tr D(1) , J(1) , t(1)
V and p at time t0 as well as their rst time derivative elds V(1) and p(1) at the same time.
From the relation D = 12 (ḞF−1 + F−t Ḟ t ), we can deduce its rst time derivative:
∗ −R
D(1) = D(1)
D
with RD = 12 (lc + lct )
(25)
in which D∗(1) and lc are dened in (17).
With (25) available, it becomes trivial to have the rst time derivative of tr D and that of D 0 ,
∗ − tr R ; D 0 = D 0 ∗ − R 0
(26)
tr D(1) = tr D(1)
(1)
D
D
(1)
From the relations b 0 = cD 0 and c = 23 k( + ) n D m−1 , we can express the rst time derivative
of b 0 as
2c(m − 1·0) 0 0
0
0
= c(1) D 0 + cD(1)
; c(1) = (Rc) +
D : D(1)
(27)
b(1)
2
3DA
where (Rc) is also dened in (17).
Similarly, the rst time derivative T(1) can be obtained by using the relations T = J bF−t and
T = J (b 0 − pI)F−t ,
0
− p(1) I)F−t
T(1) = J(1) bF−t − J bl t F−t + J (b(1)
(28)
By introducing successively (27), (26), and (25) into (28) and taking into account the relation
J(1) = J tr(l), (28) can be rewritten as
0∗
− p(1) I]F−t + J {b[tr(l)I − l t ] + [(Rc)D 0 − M : RD0 ]}F−t
(29)
T(1) = J [M : D(1)
In this relation, all terms in the right-hand side can be calculated directly from V(1) , p(1) , V and p
at time t0 .
Concerning (1) , it can be considered as a prescribed function due to the given load parameter .
In the case where is identied as time, (1) = 1.
d
, it depends on the treated boundary conditions. If t d is prescribed,
For t(1)
d
=0
t(1)
(30)
d
If b = b d is prescribed, the rst time derivative t(1)
can be deduced from the relation t d = J b d F−t N,
d
= J b d (tr (l)I − l t )F−t
t(1)
(31)
d
with (30) or (31), D(1) with (25), J(1) with
Finally, by substituting in (24) T(1) with (29), t(1)
J tr(l), and taking into consideration J = 1, F = I for the updated Lagrangian scheme, we arrive
at equations (17) after putting together all terms related to V(1) and p(1) into the left-hand side
and other terms in the right-hand side.
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 155–176 (1999)
176
H.-L. CAO AND M. POTIER-FERRY
ACKNOWLEDGEMENTS
The authors are grateful to Professor L. Toth for his kind help to improve the text.
REFERENCES
1. N. M. Wang, ‘A rigid-plastic rate-sensitive nite element method for modelling sheet metal forming processes’,
in J. F. T. Pittman, O. C. Zienkiewicz, R. D. Wood and J. M. Alexander (eds.), Numerical Analysis of Forming
Processes, Wiley, London, 1984, pp. 117–164.
2. S. Kobayashi, ‘Thermoviscoplastic analysis of metal forming problems by the nite element method’, in J. F. T. Pittman,
O. C. Zienkiewicz, R. D. Wood and J. M. Alexander (eds.), Numerical Analysis of Forming Processes, Wiley, London,
1984, pp. 45 – 69.
3. O. C. Zienkiewicz, ‘Flow formulation for numerical solution of forming processes’, in J. F. T. Pittman,
O. C. Zienkiewicz, R. D. Wood and J. M. Alexander (eds.), Numerical Analysis of Forming Processes, Wiley,
London, 1984, pp. 1– 44.
4. K. Chung and R. H. Wagoner, ‘Numerical improvement of viscoplastic, non-linear, nite element analysis’, Int.
J. Mech. Sci., 29, 45 –59 (1987).
5. C. Bohatier and J. L. Chenot, ‘Finite element formulation for non-steady-state viscoplastic deformation’, Int. J. Numer.
Meth. Engng., 21, 1697–1708 (1985).
6. A. C. Walker, ‘A nonlinear nite element analysis of shallow circular arches’, Int. J. Solids Struct., 5, 97–107 (1969).
7. J. Connor and N. Morin, ‘Perturbation techniques in the analysis of geometrically nonlinear shells’, Proc. 1970 IUTAM
Conf. High Speed Computing of Elastic Structures, Universite de Liege, 1971, pp. 683 –705.
8. L. Azrar, B. Cochelin, N. Damil and M. Potier-Ferry, ‘An asymptotic numerical method to compute the post-buckling
behaviour of elastic plates and shells’, Int. J. Numer. Meth. Engng., 36, 1251–1277 (1993).
9. Y. Yokoo, T. Nakamura and K. Uetani, ‘The incremental perturbation method for large displacement analysis of
elastic–plastic structures’, Int. J. Numer. Meth. Engng., 10, 503 –525 (1976).
10. N. Damil and M. Potier-Ferry, ‘A new method to compute perturbed bifurcation: Application to buckling of imperfect
elastic structures’, Int. J. Engng. Sci., 28, 704 –719 (1990).
11. B. Cochelin, N. Damil and M. Potier-Ferry, ‘The asymptotic numerical method, an ecient perturbation technique for
non-linear structural mechanics’, Revue Europeenne des Elements Finis, 3, 281–297 (1994).
12. M. Potier-Ferry, H. L. Cao, J. Descamps and N. Damil, ‘An asymptotic numerical method for numerical analysis of
large deformation viscoplastic problems’, Comput. Struct.; submitted.
13. N. Mole, J. L. Chenot and L. Fourment, ‘A velocity based approach including acceleration to the nite element
computation of viscoplastic problems’, Int. J. Numer. Meth. Engng., 39, 3439 –3451 (1996).
14. Y. J. Kim and D. Y. Yang, ‘A rigid-plastic nite element formulation considering the eect of geometric change and
its application to hydrostatic bulging’, Int. J. Mech. Sci., 27, 453 – 463 (1985).
15. H. Iseki, T. Jimma and T. Murota, ‘Finite element method of analysis of the hydrostatic bulging of a sheet metal’,
Bull. JSME, 20, 285 –291 (1977).
16. G. Powell and J. Simons, ‘Improved iteration strategy for non-linear structures’, Int. J. Numer. Meth. Engng., 17,
1455 –1467 (1996).
17. M. Potier-Ferry, N. Damil, B. Braikat, J. Descamps, J. M. Cadou, H. L. Cao and A. Elhage Hussein, ‘Traitement
des fortes non-linearites par la methode asymptotique numerique’, Compte Rendus de l’Academie des Sciences, 324,
171–177 (1997).
18. J. Descamps, H. L. Cao and M. Potier-Ferry, ‘An asymptotic numerical method to solve large strain viscoplastic
problems’. Proc. 5th int. Conf. on Computational Plasticity. Models Software and Applications, Barcelona, 1997,
pp. 393 – 400.
Copyright ? 1999 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 44, 155–176 (1999)
Документ
Категория
Без категории
Просмотров
2
Размер файла
167 Кб
Теги
152
1/--страниц
Пожаловаться на содержимое документа