Int. J. Pres. Ves. & Piping 58 (1994) 41-50
On the unloading solution by viscoplastic B EM of axisymmetric pressure vessels with three material hardening models Liu Yong & Hong Qichao Department of Mechanical Engineering, Zhejiang Institute of Technology, Hangzhou 310014, People's Republic of China (Received 10 February 1993; accepted 1 April 1993)
This paper presents a viscoplastic boundary element formulation and numerical implementation of the nonlinear unloading problem of axisymmetric pressure vessels with new material hardening models for which the reverse yield stress does not lie between the kinematic hardening yield stress and the isotopic hardening yield stress of present nonlinear theory. Three nonlinear unloading stress-strain models of high strength low alloy steel with both the strain hardening and Bauschinger effect factor are discussed. The modified Aitken method is used for the unloading iteration process. Comparisons of this paper's results with the data of FEM and experiment are given.
initial yield stress, .~/is mixed-hardening material parameter; ~7/= 1 corresponds to the isotopic hardening model; AT/=0 corresponds to the kinematic hardening model; 37/= - 1 corresponds to the soft model. This mixed-hardening material model requires that the reverse stress lines between the kinematic hardening yield stress and the isotopic hardening yield stress. This gives many limitations to their use in engineering, as for actual mixed material models, the unloading reverse yield stress does not lie between the yield stress of the kinematic hardening model and the isotopic hardening model (see Fig. 1). Thus, three hardening material models with both the strain hardening and Bauschinger effect have been proposed in this paper for BEM so that the BEM can be used for more and more applications in engineering. This paper will present a viscoplastic boundary element formulation and numerical implementation of the axisymmetric unloading problem. Three nonlinear unloading stress-strain material models are presented and the modified Aitken iteration method is used in this paper. A comparison of the results of BEM, FEM and experiment is given finally.
1 INTRODUCTION In previous work there have been many researches an axisymmetric problems using BEM, but almost all the papers (see Refs 1-7) concerning the axisymmetric BEM only consider the loading problem or the elastic unloading problem. This is far from true for many materials, i.e. for the high strength steel with strain hardening and Bauschinger effect. 8-1° In general, the ideas for the mixed hardening material model in FEM are as follows. The loading yield function is F = ½(s,j - aq)(Siy - a~) - loy(~P, 1~/1)
where Oy = as + f AT/day, a~ = a°(E p, A'/, ~P) As At is constant
oy( p, M) = as + MEpf
d~ p
where E p is the hardening coefficient, os is the Int. J. Pres. Ves. & Piping 0308-0161/94/$07.00 (~ 1994 Elsevier Science Limited. 41
42
Liu Yong, Hong Qichao -
-
. . . . mE
0
tractions and boundary displacements are prescribed, eqn (3) can be written as
bilinear model actual model O"
{A~ = f + D~# '' & = - A ' i + f' + D'b"
in which i is u n k n o w n quantities, and i', i" are the prescribed values of tractions and b o u n d a r y displacements. Substituting the first equation of (4) into the second equation can p r o d u c e
j .-4" ~ -
-~. ~
i= Gb"+fi 6e = K~"+ m
kinematic m o d e l ~ ' - isotopic m o d e l
Fig. 1. The actual stress-strain material model.
2 BASIC VISCOPLASTIC A X I S Y M M E T R I C BEM E X P R E S S I O N S Assuming L b o u n d a r y elements and M internal cells, the discretized version of the axisymmetric B E M equation can be written as follows: T M u q• N bp d r p} ,)
n= l
-
ejkiN p dff2 OF}m) m
~
|
n-i
(5)
where G = A - t D b,
cijtij = 2 ~ ~
(4)
K=D I-A'A
IDb,
li=A-lf lii=f'-A'A-lf
The incremental form of eqn (5) is A X = G A o ° + An Ate = KAo" + A m
(6)
Equation (6) is the f u n d a m e n t a l expression for iteration. It can be seen that eqn (6) represents a set of ordinary differential equations for stresses at related boundary nodes and internal points which can be solved by 'initial stress' m e t h o d s as used in FEM.
(1) 3 BASIC E X P R E S S I O N S OF VISCOPLASTIC S T R E S S - S T R A I N RELATIONS
m
n
eq~N p d Q 0~}"~ rn = I
(2)
.1
where the N" and N I a r e interpolation functions adopted for b o u n d a r y elements and internal cells respectively. T h e Ui*, p*, ej]t and eqk~ are the Kelvin fundamental solutions shown in Ref. 13 and in Refs 11 and 14. Applying eqns (1) and (2) to all b o u n d a r y nodes and internal points, the following matrix relationship is generated Hfi = G[~ + D"b" /¢ = G'[~ H'ti + Olo ° -
(3)
where matrices H, G, H ' , G' are the same as those obtained for pure elastic B E M analysis. D ° and D j stand for the initial stress integral matrix. It is k n o w n that when a sufficient n u m b e r of
In the usual m a n n e r for problems it is assumed that can be separated into elastic (evr) c o m p o n e n t s , so that the be expressed as:
nonlinear continua the total strain (~) (Ee) and viscoplastic total strain rate can
/~ = / L +/~vp
(7)
The total stress rate d e p e n d s on the elastic strain rate according to
where D is the elasticity matrix. T h e onset of viscoplastic behaviour is governed by a scalar yield condition of the form F ( o , evp) - E, = 0
(9)
in which F,, is the uniaxial yield stress which may itself be a function of strain hardening and Bauschinger effect parameters. It is assumed that viscoplastic flow occurs for values of F > E, only.
43
Unloading solution of axisymmetric pressure vessels Substituting eqn (7) into eqn (8), we obtain:
where
=D~,
~O=D~.~p
It is necessary to choose a specific law defining the viscoplastic strains. The simplest option is one in which the viscoplastic strain rate depends only on the current stresses, so that /~vp= f ( o )
y(~(F))~--~
(12)
in which Q=Q(a, evp, m , f ) is a 'plastic' potential and ), is a fluidity parameter controlling the plastic flow rate. The term ~ ( x ) is a positive monotonic increasing function for x > 0 and the notation ( ) implies =~(x) (~(x)) 0
{(~(x))
~ - ' = ~p + H " A ~
(11)
This relationship can be generalized to include strain hardening and temperature dependence, and the influence of state d e p e n d e n t variables, such as damage parameters for rupture theories, can also be considered. One explicit form of (11) which has wide applicability is offered by the following viscoplastic flow rule kvp =
For 0 = 0 we obtain the Euler time integration scheme which is also referred to as 'fully explicit' since the strain increment is completely determined from conditions existing at time t,. For 0 = 1-0 we obtain the 'fully implicit' scheme with the strain increment being determined from the strain rate corresponding to the end of the time interval. To define ~"+~ vvp in (16) we can use a limited Taylor series expansion and write
for for
x>0 x~<0
(13)
If we restrict ourselves to the associated viscoplastic situation, in which case F = Q , expression (12) reduces to OF kvp= 7 ( ~ ( F ) ) O-o
(14)
Different choices have been r e c o m m e n d e d for the function ~ . The two most c o m m o n versions are • (F) = e MICF-F')/F''I- 1
(15a)
• (F) = ( F - 17°1 u\ \ Fo / '
(15b)
where rl°
=
AEv~p= i~pA/. + ~AI¢'
(16c)
c" = OAt~H"
(16d)
where
4 THE THREE
MATERIAL
(16a)
MODELS
As in Fig. 1, the actual material chosen for investigation is high strength low alloy steel which shows behaviour different from the kinematic and the isotopic hardening cases and requires the mixed material model of strain hardening and Bauschinger effect factor. The experiment utilizes uniaxial tension-compression specimens for the loading and unloading processes. From a large n u m b e r of experiments, it has been shown that the stress-strain curve during loading can be replaced with sufficient accuracy by a bilinear elastic-plastic model, that is
Ee
Ae~p n+l ] " = At,[(1 - 0)/~p -{- 0evp
H°(°")
and A o n is the stress change occurring in the time interval At, = t , + l - t , . Thus eqn (16) can be written as
and
in which M and N are arbitrary prescribed constants. With the strain rate law expressed by (14) we can define a strain increment A/~v"poccurring in a time internal At, = t , + ~ - t , using an implicit time stepping scheme, as
(16b)
F = oy =
for os q- EPE p for
Ory< Os try >I as
(17)
where E is Young's modulus, E p is loading hardening constant, os is initial yield stress which corresponds to O. 1% offset. For the unloading portion, choosing a new coordinate system O* (o*, e*) with the origin at the point before unloading, the unloading relation for strain and stress is described by three kinds of model.
44
Liu Yong, Hong Qichao
(i) Bilinear elastic-plastic unloading model
F*
,=~E*e* = CrY
(a~* + EP*e p*
for for
mE
/
Cry Cr*
0"
y v
where a* = F,* = [1
+f(ePo.)](crs + EPeo~.)
eo°. is the loading plastic strain just before unloading occurs. E* = (0.9-1.0)E, E °* = m*E*/(1-m*), m* = m * ( e ~ . ) is the unloading strain hardening factor, f=f(e°o .) is the Bauschinger effect factor. The expressions for m* and f are obtained also by a large n u m b e r of experimental data, as follows. 1,~
Fig. 3. Trilinear model.
The stress-strain relationship is given as
F*
f(eP,) = 0.969918209 - 2-52133293e~ + 2.80426104e~. m*(e~.) = 0.534296201 - 0.238271643e~,
* ~- a y
for
Oy* < O*l
ok + A(EP*)"
for
as~ ~< cry* < O*z
m'E*
=
a~2 -~ -1 -- m* ep* for
(19)
Expression (19) shows a decrease in f and m* with an increasing a m o u n t of tensile plastic strain. For complicated stresses in 2D or 3D problems, the plastic strain e p, e w, eov. should be replaced by the effective plastic strains Go, G0*,
E'e*
Cry* ~
O*s2
(20) where A and ce are material constants d e t e r m i n e d by experiment, a* = [1
+ f(e~.)][Os + EPe°o.],
a~2 - (1.1021.15)o*, A = (0"2 -
(ii) Linear-Ramberg-Osgood-linear elastoplastic unloading model As in Fig. 1, the unloading bilinear m o d e l at yield point is not good e n o u g h for agreeing with actual material curve as shown in unloading. So it is necessary to use the R a m b e r g - O s g o o d function for describing a small area near the yielding point, as in Fig. 2.
The expressions for m* and f are the same as in eqn (19).
(iii) Tri-linear elastic-plastic unloading model This model is shown in Fig. 3. Its expression is written as
F* = a* =
¢r mE
0"
I E* e*
for
ay < a*l
ok +
for
a*, ~
for
O'y ~ Ors~2
, , L °*2-~miE*'ep*l -m~
Oy
(21) where o*~ = [1 +f@S*)][crs +
o
in
Fig. 2. Linear-Ramberg-Osgood-linear model.
0"~2 :
r - p * p* Os~l + ~--JI e l , ~ :
EPe°o.], 1"0
The expression for m* and f are the same as in eqn (19). Er* is the material hardening conatant given by experiment.
Unloading solution of axisymmetricpressure vessels 5 THE NUMERICAL IMPLEMENTATION OF THE ITERATION PROCESS IN UNLOADING Suppose the unloading process occurs at coordinate origin 0"(0", e*); the definitions of residual stresses and strains are I o'i, =
- o,j
S~=SiI--Sij
(22)
L E i j - 8ij -- Fit
Equation (6) together with eqns (16)-(21) can represent the unloading iteration process. The convergence condition is expressed in
[Atm~ ~vp(m)/At,~ ~vp(1)]<-Tolerl This shows that as t = Y, Atm ~
oo the
45
solutions are obtained. Hence the result is just the elastoplastic solution (see Ref. 15). The residual stresses convergence condition is: IIAo°'(m) - Aa°'(m - 1)l[= <~ Toler2 IIAo°'(m)ll~
(24)
where Tolerl and Toler2 are convergence tolerances. In the actual iterative process, for accelerating the convergence, it is necessary to used the modified Aitken method which has been used successfully in FEM. The acceleration is done in every two steps of iteration and its process is given as Aa°'(m + 1) = Ao°'(m) + &(m)
(23)
steady state
x (Aa°'(m + 1) - AcC"(m))
(25)
where &(m) is the acceleration factor,
&(m) =
{
1.
m=0,2,4,... I2Ao°'(m) - Ao°'(m + 1) - Ao°'(m - 1)]T[Ao°'(m) -- Ao°'(m - 1)] m = 1, 3, 5, • • • [ 2 A ~ ' ( m ) - Ao°'(m + 1) - A~e"(m - 1)]a[2Ao°'(m) -- A~C"(m + 1) - Ao°'(m - 1)] (26)
The unloading calculation steps of axisymmetric B E M are as follows: (i) Calculate the complete elastic unloading residual stresses (0/)). (ii) Check to see if the O~ satisfies the reverse yielding condition P
O~>oy
(v) Do the iteration loops, suppose the m loops now are: (a) calculate the Age'(m) in terms of eqn
(6) o' A~i~(m), Arvp by (b) calculate the AO(m), ' eqns (8)-(16) (c) calculate the increment of Aoijo,(m) at every Gauss point and nodes, i.e.
(27)
where f3
t
t
A<}'(m) -
¢
if yes, let R' = O[,/4 if no, stop the program and output the elastic unloading residual stresses and strains. (iii) Store the elastic results of reverse load Pmax/R', let the increment of the reverse loading step be 1(1_1 AY' = -N
AoJ(m -
(d) check to see the steady state condition and convergence condition: IAt,,~ ~(m)/At,~
~(1) <~Tolerl
and l i A r ' ( m ) - A ~ ' ( m - 1)11~
)Pmax
where N is the unloading step prescribed. (iv) Do the step loops of increment of unloading and let the last step of viscoplastic stress be the initial stress.
1)
IIA~'(m)ll~
~< Toler2
if yes, o~'= o~0'+ Ao~'(m), go to (vi) if no, using the Aitken accelerating technique as in eqns (24) and (25), then, go to (v). (vi) go to (iv), do the next unloading step till
Liu Yong, Hong Qichao
46
the unloading is zero, i.e. the reverse loading Ay' = p where Ay/)=Pm,x/R'. (vii) output the final results. For the above numerical implementation, one B E M program with loading and unloading processes has been developed for three kinds of material models.
L
....
, Pb
(a) Thick walled cylinder (BEM) Boundary elements 20 Internal elements 12 Here, the elements of FEMof cylinder are 20
6 THE RESULTS OF S T E A D Y STATE SOLUTIONS OF VISCOPLASTIC BEM The materials used in this paper are PCrNiMo and PCrNi3MoV, their stress and strain relations are based on a uniaxial tension-compression experiment, as in Figs 1-3. The main material constants are shown in Table 1 in which the explicit time integration scheme and Von Mises yield criterion are used. The viscoplastic parameters and associated constants are taken as
L2
(b) Cone and Cylinder vessel (FEM) Elements are 45
I Section Ill
I Section II
[ Section I
v = 0.001 1/day, At0 = 0-01, 0 = 0 (explicit form) Equation (15b) is used and N = 1. The samples chosen in this paper are thick walled cylinders and thick-walled cone combined cylinder shaped vessels, as shown in Fig. 4. The discretization of B E M and F E M for the samples is shown in Fig. 4; there are two group of samples in the calculation and experiment, one group is the thick-walled cylinders, another group is the vessels of cone combined cylinders.
(i) Group one: thick walled cylinders The sizes and internal pressures of the samples are: No. 1 Pm,x = 555"0 MPa, p~,/p~ = 1.80 No. 2 Pm,x = 607"0 MPa, Pt,/P, = 1.96 No. 3 Pma× 533-0 MPa, Pb/P, = 1.97 No. 4 Pm,x = 671"3 MPa, Ph/Pa = 2"4 L = 5 0 m m , Pa = 15 m m , Pb/P~ is the thick wall ratio, where No. 1-No. 3 samples are of steel PCrNiMo, No. 4 sample is of steel PCrNi3MoV. The experimental condition is open ended for cylinders. The results are shown in Figs 5-8 and Table 2. =
,
I
(c)
Cone and Cylinder vessel (BEM) Boundary elements 16 Internal elements 12
Fig. 4. The samples and meshes.
Figures 5-7 show the curves of pressure versus strain at the outside surface of the tube with BEM, F E M and experiment. Figure 8 shows the comparison of three unloading material models. Table 2 shows the hoop residual stresses at the outside surface with experimental results and three material models. From Figs 5 - 8 and Table 2, we can see that the results of B E M , F E M and experiment shows good agreement. The error of the bi-linear model is bigger, but this model is simple and easy to use in engineering. The error of the linear-R.O.linear model which approaches the actual material model is the smallest.
Table 1. Mechanical properties of steels PCrNi3MoV and PCrNiMo Material PCrNiMo PCrNi3MoV
E(MPa)
v
o,(MPa)
EP(MPa)
o<
EP*(MPa)
E~°*(MPa)
2-0 x 105 2.06 x 10s
0.3 0-3
850-0 1080-0
1380.0 1250-0
0.3 0-57
6100.0 1-86 x 10s
1-2 x 105 1.9 x 105
Unloading solution of axisymmetric pressure vessels
-" -" -= BEM (linear-R.O.-linear model) m- -m---m FEM f o o o Exp.
500
-" -" =- BEM
500
o
47
(bi-linear
model)
o o Exp.
/
400 375 300 fl_ v 0..
250 n
200
125 100 0
t/
i
0
400
.!
i
I
800
I
1200 1600 EOb (~te)
I
I
2000
2400
I
L
0,1/t 0
Fig. 5. p v e r s u s eob f o r N o . 1.
.¢
I 500
I I 1000 1500 e Qb (~e)
I 2000
I 2500
Fig. 7. p v e r s u s eob o f N o . 3.
600 " •
•
" BEM (Tri-linear model) • Exp.
z/,
500
600 -
500
i
-
# E
_
~ h
400
400
0_
: : -_ bi-linear model BEM (2) . . . . . linear-R.O.-linear model BEM
~ (3) • / 7 °
• • Tri-linear model BEM o o Exp. results
///~
300 a_
3o0
200
200
/
/
/
100 100
0
500
Fig. 6. p
1000 1500 e Qb(~te) versus
eoh
2000
2500
f o r No. 2.
0
400
800 e Qb(pe)
1200
1600
Fig. 8. p v e r s u s Eob o f N o . 4.
(ii) Group two: cone combined cylinder vessels The sizes of this t y p e o f vessel are L = 115 m m , L , = 40 m m , L2 = 75 m m p~ = 30 m m , Pa, = 15 m m , P,2 = 11-5 m m , p m , x = 764.9 MPa
There are six samples used in the computation, their material is PCrNi3MoV. The viscoplastic parameters and associated constants are as in the above example (Group (i)). The results are shown in Figs 9-12. The material hardening model chosen for investigation is a bi-linear
Liu Yong, Hong Qichao
48
Table 2. The hoop residual stresses at outer surface unit (MPa); test results compared with results from analyses Sample No.
Bi-linear model Tri-linear model Linear-R.O.-linear model Exp. results ~
--x IOMPa 0 120
2
3
4
Average errors compared with exp. results '6
117 112
96 93
64 58
40 37
21% 13%
110 101
90 81
55 50
36 33
9.7% --
FEM
- 10
1
p = 22.5mm
Cone part : 2 0 - 5 5 m m
BEM 20
30
40
50
60
70_
80
80
--q
40
J
OZ
"[pz
--z(mm)
0
Gp -40
-80
Fig. 9. The loading end stresses of No. 4 sample along the axis at radius p = 22.5 mm.
x IOMPa 10
x
20
30
40
50
60
70
80
6O
40
20
O' z z(mm)
0 O'p
-20 i -40
l ---
FEM
- -
BEM
Fig. 10. T h e u n l o a d i n g residual stresses of No. 4 along t h e axis at radius p = 22-5 m m .
Unloading solution o f axisymmetric pressure vessels
49
model and because of space limitations we can only give here the results of No. 4 sample. Figures 9-12 show the steady state stresses of No. 4 sample along the radal direction and axial direction both at the loading end and after unloading computed by BEM and FEM. From these results we can see that both the BEM and FEM have good agreement.
x IOMPa 100
80
60
40 a b c
a : section b : section c : section
2O 1.5
0
15
7 CONCLUSIONS
I! III
20
25
30 mm)
a
-20
-40
b
z
c
/
/
. . . . . FEM BEM
-60
-80
Fig. 11. The loading end stresses o v, o~) at three sections.
x IOMPa c' 0 6O a
b
c
40
20
0
-20
-40
The nonlinear unloading problem is a very complicated problem which is of great interest in actual engineering. Since the actual unloading material model is different from the kinematic and isotopic hardening model, its expression is more complicated than the loading material model, so the present paper first obtained the three actual unloading material models which consider both the strain hardening and Bauschinget effect factors. Then this paper presents an axisymmetric nonlinear BEM implementation with three actual material modes for the unloading problem. The modified Aitken method is used for accelerating the iteration process. From the above results and discussion, it can be seen that the results of BEM, FEM and experiment are in good agreement for the unloading problem which shows the effectiveness of the paper's method. The linear-R.O.-linear unloading material model is closest in its approach to the experimental results; the trilinear model is also good for approaching with experimental results. The bilinear unloading model produces the biggest errors of the three models, but it is simple and easy to use in engineering.
ACKNOWLEDGEMENT
r /
-60
/
This research project was supported by the National Natural Science Foundation of China (No. 19202018).
/ ~
~M
/
-80 t
-10011.5
ab :: section section II[ c : section II| ll5
I 20
[ I 25
(ram) 30
Fig. 12. The unloading residual stresses o'p, try')at three sections.
REFERENCES 1. Kermamidis, T., A numerical solution for axially symmetrical elasticity problems, Int. J. Solids Struct., 11 (1975) 493-500. 2. Cruse, T. A., Snow, D. W, & Wilson, R. B. Numerical
50
3.
4. 5.
6.
7.
8. 9.
Liu Yong, Hong Qichao solutions in axisymmetric elasticity. Comp. Struct., 7 (1977) 445-51. Cathie, D. N. & Banerjee, P. K., Boundary element methods for axisymmetric plasticity. In Innovative Numerical Methods for the Applied Engineering Science, ed. R. P. Shaw et al. University of Virginia Press, 1980. Sarihan, V. Mukherjee, S., Axisymmetric viscoplastic deformation by boundary element method'. Int. J. Solids Struct., 18 (1982) 1113-28. Henry, D. P. & Banerjee, P. K., A substructured boundary element analysis of thermo-plastic axisymmetric bodies. J. Engng Mech. ASCE, 113 (12) (1987) 937-56. Henry, D. P., Pape, D. A. & Banerjee, P. K., New axisymmetric BEM formulation for body forces using particular integrals. J. Engng Mech. ASCE, 113 (1987) 671-88. Henry, D. P. & Banerjee, P. K., A variable stiffness type boundary element formulation for axisymmetric elastoplastic media. Int. J. Num. Meth. Engng, 26 (1988) 1005-27. Chen, P. C. T., The Bauschinger and hardening effect on residual stress in an autofrettaged thick walled cylinder. J. Pres. Ves. Technical., 108 (1986) 108-12. Milligan, R. V. et al., The Bauschinger effect in a high strength steel. J. Basic Engng, $8 (1966) 480-8.
10. Liu, Yong & Shen, Naijie, Residual stress analysis for actual material model of autofrettaged tube by nonlinear BEM. Int. J. Pres. Ves. & Piping, 48 (1) (1991) 21-36. 11. Liu Yong, PhD thesis, East China Institute of Technology, People's Republic of China, 1990. 12. Liu, Y. & Shen, N. J., The residual stress and strain solutions of autofrettaged pressure vessels with a cone and cylinder connection. J. Strain Anal. 27 (1) (1992) 7-14. 13. Brebbia, C. A., et al., Boundary Element Techniques: Theory and Applications in Engineering. Springer° Verlag, Berlin, 1983. 14. Telles, J. C. F., A boundary element formulation for axisymmetric plasticity. Proc. 5th Int. Conf. on BEM, Pergamon, 1983. 15. Liu, Yong & Shen, N. J., A BEM implementation for elastoplastic analysis by using the steady state solutions of viscoplastic problems. The computational mechanics and its applications in engineering. Proceedings of 3rd Computational Mechanics Conference of East China. China University of Science & Technology Ltd, 1992, pp. 41-6. 16. Feng, G. B. Masters thesis, East China Institute of Technology, People's Republic of China, 1990.