- Email: [email protected]
Application of finite element method to ring compression test
Nuclear Engineering and Design 44 (1977) 75-85 © North-Holland Publishing Company
75
APPLICATION OF FINITE ELEMENT METHOD TO RING COMPRESSION TEST Rayji YAMADA, Motoe SUZUKI and Yasuo HARAYAMA Japan Atomic Energy Research Institute, Tokai Research Establishment, Tokai-mura, Naka-gun, lbaraki.ken, Japan Received 3 April 1977
The relation between deflection and load during the ring compression test used to study mechanical properties of a fuel cladding tube is examined by a computer program which deals with the elasto-plastic deformation by the two-dimensional finite element method. Under Hill's criterion, in the early stage of ring deformation the results calculated by the plane strain condition, and in the late stage the results calculated by the plane stress condition, agree well with the experimental ones. It is revealed that plastic deformation occurs even at the initial stage of the test, and crack locations are suggested. The relation between the ring compression test and the uniaxial tension test are also discussed.
1. Introduction The four experimental ways to examine the mechanical properties of zircaloy cladding tubes of a light water reactor's fuel pins during pre- and post-irradiation are: (1) the tensile test [1,2] ; (2) the ring tension test [1] ;(3) the burst test [1,2]; and (4) the ring compression test. The ring compression test has merits of easy fabrication of specimens and simple testing procedures; these are particularly important points for the post-irradiation test, because activated tubes to be cut for the test are difficult to handle in hot laboratories. The ring compression test is considered to be suitable for the evaluation of strength and elongation ' where tension and compression exist, because stresses become tensile or compressive, depending on the location of the tubing in one specimen, and it is also used to evaluate the brittleness of tubes caused by oxidation in steam in loss-of-coolant accidents. That is to say, the maximum deflection of an oxidized zircaloy tube [3], or the number of cracks at a given deflection [4], or whether the ring deformation is plastic or not at the occurrence of a crack [5], are measured to detect the ductility of the tube. As the stress and strain distributions of a tube at each deformation stage are not assessed by the ring compression test, we intend to use the computer
simulation to understand the phenomena that gives a qualitative explanation of the behaviour of a deformed specimen, and emphasizes the relation between the experimental results and the mechanical properties. As the incremental theory gives good correspondence between theory and experiments [6], the finite element method is a powerful tool to solve the elasto-plastic problems [7,8]. We have developed the new program to simulate the ring compression test based on the two-dimensional elasto-plastic computer programs described in refs. [9,10], see figs. 1 and 2. We experimented uniaxial tension test to obtain the work-hardening coefficients used for the program's input data, and the ring compression test to check the program's fitness. 2. Experimental The tube subjected to the ring compression test was cold-rolled and stress-relieved zircaloy-4. The specimens, taken from commercially fabricated tubes, had 10.75 mm o.d., 9.30 mm i.d., and 15 mm length. The main compositions except Zr (weight percent) were Sn 1.46, Fe 0.22, Cr 0.10, O 0.128, C 0.008 and H 0.002. The specimens were annealed at 1000°C for 10 min and thereafter quenched in water. The following four test conditions were chosen: the crosshead speeds of testing machine (Instron-type) were
76
R. Yamada et al. / Application o f finite element method to ring compression test
set at 20 mm/rrlin or 5 mm/min, depending on whether molybdenum disulfide (MoS2) paste was used or not. As there was little difference between them, friction and relaxation were neglected. The load was measured by a load cell using resistive wire. The amount of deflection d and flattening value f are defined as follows, d = D - h and f = (D - h)/D, where D is the initial fabricated outer diameter and h is the distance between the two crossheads. The experimental results for the ring compression test are shown in figs. 3 - 6 by the solid lines (curves 1). The load per unit length of tubing is plotted against the amount of deflection. An additional experiment was carried out in order to measure the uniaxial tensile properties of the workhardening coefficients. The specimens were plates of 1.0 mm thickness and 20 mm gauge length. Their heat treatment was the same as for the foregoing ring compression test. The test was performed by the lnstrontype tension machine at the crosshead speed of 5 mm/min. The stress and plastic strain (o - ep) curve obtained was divided into five regions in order to be fitted by the exponential function, namely o = X ( Y + ep) n . In each region the fitting parameters were calculated by the least squares method. The yield stress of the specimen was ~35.5 kg/mm 2 .
F(Gdetp - HdeaP)2+ G(Hdeap - FdeP) 2 +H(Fderp - GderP)2 X
(FG + G H + H F ) 2
+ d')'taP2 + d'r~P~ + d"/rp211/2 2L
Zircaloy cladding tube is well known as an anisotropic material. For an anisotropic material that possesses three mutually orthogonal planes of symmetry, Hill's criterion is used as the yield criterion. If its radial, tangential and axial directions are the principal axes for the anisotropy in the cylindrical coordinates (r,t,a), Hill's criterion is given in ref. [11] as °=
//2
3 (F+G+H)
[F(°t -- O'a)2 + G(°a -- or)2
-2N-'-]
{ao)
=
[De
I -
---
(1)
where ~-is the equivalent yield stress; suffixes r, t and represent radial, tangential and axial directions, respectively. F, G, H, L, M and N are anisotropy parameters. The equivalent plastic strain increment d~-p related to ~-is defined in the following way: d~-p = V 2 ( F + G + H) 3
(2)
---
-
-= [Dp] {2~e),
(3)
where [De] is the elastic stress-strain matrix and W = do/dep is the work-hardening coefficient dependent on plastic strains. If we express {Ae} and (Ao} by the displacement matrix [B] and [Dp], they become (Ae} = [B] CA/i} and (Ao) = [Dp] (Ae}, respectively. The incremental form of the finite element equations of equilibrium, neglecting body forces and initial strains, are expressed [7] as
=f[Bl T [Dpl [BId V. (aS}
= [KI {a63,
(4)
where the integral is taken over the whole volume of each element. As the radial and tangential directions in a given element, which must be chosen for the principal axes of the anisotropy, vary with continuous deformation of the ring, the finite element equations of equilibrium require to be represented initially in the local coordinates (r,t,a). We transform them to the global coordinate (x,y, z) and sum them over all dements. Let the angular transformation matrix be [T], then the relation between local and global coordinates becomes {A6~) = [T] {A~g}
+ H(ar - 002 + 2(Lr2a + g'ra2r + Nr2t)] 1/2,
"
We assume that Hill's criterion is applicable to the zircaloy tube [12]. The relation between stress increments (Ao} and strain increments {Ae} is represented by [8]
{A/) 3. Calculation
2M
and {AFt} = [T] CA/g),
(5)
where suffixes £ and g represent local and global coordinates, respectively. From eqs. (4) and (5), we obtain {Z~g) = I T ] - l [K] [TI ( a 6 g ) .
(6)
The summation of eq. (6) over the whole elements gives the total euqations of equilibrium. The matrix {ASg} gives the stress increments {Ao~} as {Ao~) = [Dp] [B] [T] {A6g} .
(7)
77
R. Yamada et al. / Application of finite element method to ring compression test Geometry of two dimensional
i-1+ t i ]
t i
inite element
Yc- Yc
Z~Ycl
) t i i
Load
ty~ "~ rr-Ye~ y~-1_y i - - - ~ l z--,-dy~ '= - AtY~b,(i=m)
ZircaLoytube dimension inner diam. 9.30 mrn outer diam. 10.75 mm
~
i sfep number t: temporary value Y: y- coordinate F : nodal force in the
direcfionof y axis
Ay?= at..~yjrn, J=o,b,c FOrced disPlQcement of the node b and c
e>O
x
fes
~,u.,j
,
Fig. 1. Schematic diagram for calculation by the finite element method. or, at and oa involved in eq. (1) are obtained by the summation of eq. (7). We assume two-dimensional plane conditions for the calculation of eq. (6); namely, the plane stress condition (Oa = 0) and the plane strain condition (% = 0). Triangular elements are adopted for finite elements whose grids and fixed conditions are shown in fig. 1. When we apply the finite element approach to the ring compression test, the most important problem is how to handle the contact between the tube and the crosshead of the compression tester. From eq. (4) the incremental method can obviously reduce the nonlinear problem, such as the plastic flow, to the linear problem. The simulation of the ring compression test using the linearity is as follows: eq. (6) is solved to obtain the temporary nodal forces and displacements by applying a certain negative force displacement increment to the constrained nodal points that take a maximum y-value. If y-values of some free nodal points exceed the constrained one at a certain deformation stage, one proportional coefficient (ix) is calculated to equalize them, and then the temporary values are multiplied by it in order to obtain the true values. These free nodal points are dealt with constrained ones from the next stage. If forces of the constrained nodal points become positive at a certain stage, another proportional coefficient (~]) is calculated in order to make the forces zero, and then these constrained nodal points are converted into the free points from the next stage. The computer block dia-
I
dF=t , (i=n) ~' . _ _ ia t I1 Ayj =~'. Ayj, j=a,b,c
I
Free the node c
Fig. 2. Block diagam to flatten nodal points on the tube surface. gram for three nodal points (Yc >Yb > Ya at initial stage) are given in fig. 2. The proportional coefficients a and/3 can exclude the recalculation of the contact problem in that example.
4. Results and discussion 4.1. The effects o f anisotropy on the results o f simulation
The results by which we simulate the ring compression test are strongly dependent on the number of nodal points and triangular elements. These results approach the experimental ones and increase with the number of elements and nodal points. The dashed curves 2 in figs. 3 and 4 show the results of the simulation adopting the van Mises yield criterion under plane strain and plane stress conditions, respectively. It is seen that these results become less than the experimental ones in the initial plastic deformation region. The elastic solution of the ring compression is given analytically [13] as
12e/(bE),
(83
78
R. Yamada et aL /Application o r finite element method to ring compression test
I
I
'
E
[
I
Plane Strain
~:, ~} '~'
Loadl,"
•
-
1") Exp.
Exp [sofropic a,
Plane strain ",c v c-
20 E m
~) H=O.5F,G:O5F. N : 2 OF
J
H=O5F,G=I.0F.N=20F 4 , H:O.5F,G= 15F, N = 2 OF
2O
4-
!
E
lO
A,
/
,
,
~
,,j)
10
d
i
0
2
4
6
Deflection
J
Plane stress
I
'
]
Plane stress
] sotropic
(mm)
range, but deviate at the advanced one, whereas the calculated with smaller G value (i.e. G = 0.5 F) approaches the experimental ones. These results suggest that anisotropy parameters should take different values depending on the degree of plastic deformation. Maki and Ooyama have measured the plastic deformation and reported that the anisotropy of zircaloy-2 changes with the progress of plastic deformation [14]. As they did not reveal values of the anisotropy parameters, we estimated values from their experimental data and simulated using Hill's criterion with the estimated values shown in figs. 3 and 4 (curves 3). The estimated anisotropy parameters do not appear to give good results in the initial plastic region under
i
~, Exp
:2)
6
Fig. 5. Dependence of calculated results on anisotropy parameters under the plane strain condition.
F
E
4 Deflection
where P is the load in the deflected direction, b and T are the length and the thickness of the tube, respectively, and R is the mean radius of the tube. The result obtained by eq. (8) is indicated in figs. 3 and 4 by a dotted line. It is seen in the elastic region that the finite element method results agree well with the analytical elastic solution. The effects of the anisotropy parameters on the results of simulation are examined by varying the values of F, G, H a n d N divided by F. These results, shown in figs. 5 and 6, indicate that if G becomes larger than H and F (i.e. G = 1.5 F) under the plane strain condition, the calculated results agree with the experinaental ones in the initial plastic deformation
A
2
(mm)
Fig. 3. Calculated results under the plane strain condition. (a) Anisotropy parameters estimated from the data of ref. [14]; (b) anisotropy parameters listed in table 1.
I
0
8
:~
'
Exp.
I
]
,~ H:O.5F.G=O.5F.N=2OF
E
(~
H = O.5F, G= 10F, N=2.0F
4)
H = 0 5 F . G= 1.SF, N = 2 0 F
,~4;, b 20
:E D ; L&,®
cc-r
__m t-
E
~ lO
!}
1o
123
3
_J ,
I
2
,
I
4 Deflection
i
I
i
6 ( mm )
Fig. 4. Calculated results under the plane stress condition.
i
0
I
i
2
4 Deflection
I
J
I
6 (mm)
-Fig. 6. Dependence of calculated results on anisotropy parameters under the plane stress condition.
R. Yamada et al. /Application of finite element method to ring compression test
Table 1 Anisotropy parameters used in the present work
~p < 1 × 10 - 5 1 X 10 -5 ~
H
G
N
0.5F 0.5F 0.5F 0.5F 0.5F
1.5F 1.0F 0.5F OAF 0.3F
2.0F 2.0F 2.0F 2.0F 2.0F
b o t h the plane strain and stress conditions. From the comparison of anisotropy parameters' dependence in figs. 5 and 6, we chose the values listed in table 1 where the G value becomes larger than the one estimated from Maki and Ooyama's data in the range
79
below 10 - 3 of % . It is understood from the results in figs. 3 and 4 (curves 4) that in the early stage of ring compression the calculated values under the plane strain condition and in the advanced stage, the calculated values under the plane stress condition, agree with the experimental ones. The reasons for these conditions are considered in the following section. 4.2. The stress and strain distributions
The stress and strain distributions in the flattened tube during the deformation in figs. 3 and 4 (curves 4) are shown in figs. 7 - 1 0 , A and B in the figures indicate the distributions under the plane strain condition whereas C and D indicate those under the plane stress condition. The plastic region at each deflection is dis-
A
Fig. 7, Contour map of er distribution. The unit for st/ain is a fraction. A,B,C and D represent the contour maps at each deflection indicated by the same characters as in figs. 3 and 4. These characters are commonly used in figs. 7-11.
R. Yamada et aL /Application o f finite element method to ring compression test
80
2
~
4
D
"
2
:
~
?
r~ ~
B
CONTOUR 1 ? 2 6 3 4 4 3 S ] 6 O ? ] 8 3 9 4 10 G 11 ?
V~LUE SO00 E + 1 OOOO E + 1 SO00 E + 1 O000 E + 1 SO0@ E + 1 O000 SO00 E + 1 @000 E + 1 SO00 E + 1 0000 E + ] SO00 E + 1
SiOM~-T GEOMETRY
CONTOUR SoO E - 1
INE ~__j
Fig. 8. Contour map of o t distribution. The unit for stress is kg/mm 2.
palyed by Y in fig. 11. It shows that the spread of the plastic area coincides with the increase in the amount of deflection. It is to be noted that the regions inside and outside of the tube near the x and y coordinates are already plastic even at small amount of deflection where the load is linearly dependent on deflection, that is, point A of the deflection. The radial stresses are nearly zero except the compressive stress which occurs at the region of contact with the crosshead, and the small tensile and compressive stresses which occur at the outside and inside of the neutral surface, respectively, at the severe deflection. The er distribution in fig. 7 is remarkably large compared with the or distribution because of the influence of large ot. As the o a in the initial deformation take small values that do not exceed the yield
stress (~35.5 kg/mm2), the tensile and compressive stresses are considered to be cancelled out macroscopic ally in order that the plane strain condition can be fulfilled in that range. In the advanced deformation, however, the increases exceed the yield stress and the tube becomes plastically deformed accompanied by the stress relaxation that is thought to satisfy the plane stress condition. If a criterion for the occurrence of crack is defined as more than 10% tensile strain, then cracks will occur along the radial direction from the outer surface of tube (et > 0.1) around the x coordinate, or along the tangential (Ca > 0.1) or axial (er > 0.1 ) direction from the inner surface around it. The experimental evidence from the ring compression test shows that cracks occurred along the radial direction at the outer
81
R. Yamada et aL / Application or finite element method to ring compression test
C
~
~
~
l
B
@ .
"
\//~ ,1
Fig. 9. Contour map of et distribution. surface of the tube along the x coordinate and were also found to stretch along the axial direction. These were thought to be formed by the accumulation of many small cracks along the tangential direction at the inner surface. It appears that more than 10% tensile strain can be a criterion for the occurrence of cracks. It is worth noticing that the maximum tensile strain which results in the occurrence of a crack is taken by the tangential strain at the outer surface of the tube in the x coordinate of fig. 1. The variance of the tangential strain of the finite element at that surface against the flattening value f i s shown in fig. 12. The computed strain takes a larger value than the experimental result of Ueki et al. [15] at large elonga-
tion. The reason is thought to be that the surface elongation measured by strain gauge gives mean values averaged by the larger area to be detected than the area of the finite element (i.e. the length of the outer surface element is about 0.1 mm).
4. 3. Estimation of yieM stress and maximum elongation The ring compression test cannot directly give the yield stress and maximum uniaxial elongation up to fracture obtained by the uniaxial tension test. Because it requires more difficult specimen fabrication procedures compared with the ring compression test described in the introduction, it is useful if the ring
82
R. Yamada et aL / Application o f finite element method to ring compression test
~St~'TOU£ L J~,lF = . 0 £ 1 ........ Fig. 10. Contour map of oa and ea distributions. compression test is used to estimate the yield and maximum elongation. We investigated the possibility in order to estimate them by using the perfect elastoplastic model (i.e. non-work-hardening in the plastic range) for simplicity. The amounts of loads, where the linearity between deflection and load disappears, are estimated by our program by varying the inputs of yield stresses. The results are shown in fig. 13 where loads are normalized to correspond with eq. (8). It may be said that the mean tensile stress during plastic deformation in the uniaxial tension test can be found by the normalized load at the appearance of nonlinearity in the ring compression test from fig. 13 because it shows fairly good correspondence between the estimation and the foregoing experimental results that give ~2000 kg/mm for the normalized load and
~50 kg/mm 2 for the average work-hardening stress. If the fracture energy is estimated by fO-dgv, the maximum elongation is considered to be calculated as follows. As the deflection at the fracture is known by experiment, the amount of f~d~-p up to it at the fractured element is calculated by our program. It should be noted that there is no difference in choosing an element through the calculation of fod~'p if the stresses and strains are converted basically into the equivalent stress ~, and plastic strain fp. For the fractured element, we chose the outer surface element at the x-coordinate. We estimated the flattening value at the fracture and the maximum equivalent plastic strain and the maximum tangential strain in each yield stress of the perfect elasto-plastic model, by assuming the fractured deformation stage to possess the equiva-
R. Yamada et aL /Application of finite element method to ring compression test
g g g~g~#~gYyg'~9'~y~'~t~F~yyY~yyY
Y~,y* Yy rY rr yY Yr yy y yY 't ~"#
,,-z
~#~Y ~.Y~
~
-
~ . _
~ L £ . ,
-
r.
-
~.
"~7
-
83
' > e
' ,'
~ \\
<~7(.. .
,
•
.
~.,-.~:.
:~"~'Y
f
".f
.
,
, ....,
".".'
....
'~: : : :
A ""
"~N Y
".\\
\\ \
\\
,\ \\
\
\
\
7
OEONETRY
Fig. 11. Yield region of flattened tube.
I
,50
ko c-
~,0
././/o,
-e
,oo
'
I
~ mb-r
/°
I
J°/°'_
/°
'-'
'
5O E ,
20 Flottening vQlue (%}
~0
Fig. 12. The tangential strain at the outer surface of tube vs. flattening value,
Z
[ 40 Yield
i
[
i
60 stress
80 Kg/mm 2
Fig. 13, Normalized load vs. yield stress in the perfect elastoplastic model.
R. Yamada et al. / Application o f finite element method to ring compression test
84
I
I
'
I o
~ao-
80~
~ 60-
60 m
(5) The possibility of estimating the mean stress and maximum elongation for a uniaxial tension test is suggested by the simulation of the ring compression test using the perfect elasto-plastic model.
t~
m
Acknowledgement r~
g~o-
~E
20-
20 5
1
I
i
The authors express their thanks to Mr. A. Morishima and Dr. S. Kawasaki for their encouragement and useful discussion, and also to Mr. M. Ono, CRC Co. Ltd., for programming of contour maps.
I
0 50 70 Flattening value at fracture (%) F'ig. 14. The m a x i m u m t a n g e n t i a l s t r a i n a n d e q u i v a l e n t plastic strain o f the finite e l e m e n t at the f r a c t u r e o f t u b e in the perfect elastic-plastic m o d e l .
lent fracture energy of the former fod~p. This estimation is shown in fig. 14. We think that the maximum tangential strain at the outer surface of the tube, which may have close connection to the maximum elongation of the uniaxial tension test, is ably obtained by fig. 14, knowing experimentally the flattening value of the tube at the occurrence of crack.
5. Conclusion
Nomenclature b d D [Del [Dpl E
f
length of tube displacement matrix d e f l e c t i o n of t u b e initial f a b r i c a t e d d i a m e t e r elastic s t r e s s - s t r a i n m a t r i x plastic s t r e s s - s t r a i n m a t r i x Young's modulus = flattening value o f t u b e
F,G,H,
L,M,N h P r, t, a R
x, y, z
We attempted to simulate qualitatively the ring compression test by developing a computer program based on the two-dimensional elasto-plastic finite element method. The results obtained are as follows. (1) The isotropic yield criterion (i.e. the yon Mises yield criterion) gives a little difference, whereas use of the plastic strain dependence from Hill's anisotropy parameters gives a good fit with experiment. (2) In the early stage of ring deformation the results calculated under the plane strain condition and in the late stage the results under the plane stress condition agree with experiment. It is considered thai the axial stress strongly influences these results. (3) A crack can be expected to occur in a location where there is more than 10% tensile strain. (4) A small portion of tube is plastically deformed even when linearity exists between load and deflection.
= = = = = = =
x,r,n T
[T] W
(} LJ
= a n i s o t r o p y p a r a m e t e r s of Hill's y i e l d c r i t e r i o n = d i s t a n c e b e t w e e n the t w o c r o s s h e a d s of the compression testing machine = load w e i g h t e d to t u b e = radial, t a n g e n t i a l , a x i a l d i r e c t i o n s in the cylindrical c o o r d i n a t e s = m e a n r a d i u s of t u b e = Cartesian coordinates = f i t t i n g p a r a m e t e r s in the e x p o n e n t i a l f u n c t i o n = t h i c k n e s s of t u b e = angular t r a n s f o r m a t i o n m a t r i x = work-hardening coefficients = c o l u m n vector = r o w vector
Greek symbols = n o r m a l stress and s t r a i n of ith c o m p o n e n t oi, ei = shear stress a n d strain rii, 7# = e q u i v a l e n t stress and strain = nodal displacement increments = n o d a l force i n c r e m e n t s = stress i n c r e m e n t s = strain i n c r e m e n t s = proportional coefficients Subscripts and superscripts e = elastic p = plastic = local c o o r d i n a t e s g = global coordinates.
R. Yamada et al. / Application o f finite element method to ring compression test
References [1] D.O. Pickman, Nucl. Eng. Des. 21 (1972) 212. [2] A.A. Bauer, L.M. Lowry and J.S. Perrin, BMI-NUREG1948 (1976). [3] S. Kawasaki, T. Furuta and M. Hashimoto, JAERI-M 6181 (1975) (in Japanese). [4] D.O. Hobson and P.L. Rittenhouse, Oak Ridge Report, ORNL4758 (1972). [5 ] D.H. Bradhurst and P.M. Heuer, J. Nucl. Mater. 55 (1975) 311. [6] R. Hill, The Mathematical Theory of Plasticity (Oxford University Press, 1950) p. 38. [7] Y. Yamada, N. Yoshimura and T. Sakurai, Int. J. Mech. Sci. 10 (1968) 343.
85
[8] R.H. GaUagher, Y. Yamada and J.T. Oden, Recent Advances in Matrix Method of Structural Analysis and Design (The University of Alabama Press, 1971) p. 283. [9] O.C. Zienkiewicz, Finite Element Method in Engineering Science (McGraw-Hill Book Co., 1971) p. 435. [10] Y. Yamada, Plasticity and Viscoelasticity (Baihukan, (1972) p. 173 (in Japanese). [11] R. HiU, ref. [6], p. 317. [12] Y. Miyamoto, Y. Komatsu, T. Kakuma and N. Nagai, J. Nucl. Mater. 61 (1976) 53. [13] S. Timoshenko, Strength of Materials, Part 1 (D. Van Nostrand Co., 1955) p. 381. [14] H. Maki and M. Ooyama, J. Nucl. Sci. Technol. 13 (1976) 43. [15] I. Ueki, K; Fukai and H. Ichikawa, J. At. Energy Soc. Jpn. 18 (1976) 591 (in Japanese).
75
APPLICATION OF FINITE ELEMENT METHOD TO RING COMPRESSION TEST Rayji YAMADA, Motoe SUZUKI and Yasuo HARAYAMA Japan Atomic Energy Research Institute, Tokai Research Establishment, Tokai-mura, Naka-gun, lbaraki.ken, Japan Received 3 April 1977
The relation between deflection and load during the ring compression test used to study mechanical properties of a fuel cladding tube is examined by a computer program which deals with the elasto-plastic deformation by the two-dimensional finite element method. Under Hill's criterion, in the early stage of ring deformation the results calculated by the plane strain condition, and in the late stage the results calculated by the plane stress condition, agree well with the experimental ones. It is revealed that plastic deformation occurs even at the initial stage of the test, and crack locations are suggested. The relation between the ring compression test and the uniaxial tension test are also discussed.
1. Introduction The four experimental ways to examine the mechanical properties of zircaloy cladding tubes of a light water reactor's fuel pins during pre- and post-irradiation are: (1) the tensile test [1,2] ; (2) the ring tension test [1] ;(3) the burst test [1,2]; and (4) the ring compression test. The ring compression test has merits of easy fabrication of specimens and simple testing procedures; these are particularly important points for the post-irradiation test, because activated tubes to be cut for the test are difficult to handle in hot laboratories. The ring compression test is considered to be suitable for the evaluation of strength and elongation ' where tension and compression exist, because stresses become tensile or compressive, depending on the location of the tubing in one specimen, and it is also used to evaluate the brittleness of tubes caused by oxidation in steam in loss-of-coolant accidents. That is to say, the maximum deflection of an oxidized zircaloy tube [3], or the number of cracks at a given deflection [4], or whether the ring deformation is plastic or not at the occurrence of a crack [5], are measured to detect the ductility of the tube. As the stress and strain distributions of a tube at each deformation stage are not assessed by the ring compression test, we intend to use the computer
simulation to understand the phenomena that gives a qualitative explanation of the behaviour of a deformed specimen, and emphasizes the relation between the experimental results and the mechanical properties. As the incremental theory gives good correspondence between theory and experiments [6], the finite element method is a powerful tool to solve the elasto-plastic problems [7,8]. We have developed the new program to simulate the ring compression test based on the two-dimensional elasto-plastic computer programs described in refs. [9,10], see figs. 1 and 2. We experimented uniaxial tension test to obtain the work-hardening coefficients used for the program's input data, and the ring compression test to check the program's fitness. 2. Experimental The tube subjected to the ring compression test was cold-rolled and stress-relieved zircaloy-4. The specimens, taken from commercially fabricated tubes, had 10.75 mm o.d., 9.30 mm i.d., and 15 mm length. The main compositions except Zr (weight percent) were Sn 1.46, Fe 0.22, Cr 0.10, O 0.128, C 0.008 and H 0.002. The specimens were annealed at 1000°C for 10 min and thereafter quenched in water. The following four test conditions were chosen: the crosshead speeds of testing machine (Instron-type) were
76
R. Yamada et al. / Application o f finite element method to ring compression test
set at 20 mm/rrlin or 5 mm/min, depending on whether molybdenum disulfide (MoS2) paste was used or not. As there was little difference between them, friction and relaxation were neglected. The load was measured by a load cell using resistive wire. The amount of deflection d and flattening value f are defined as follows, d = D - h and f = (D - h)/D, where D is the initial fabricated outer diameter and h is the distance between the two crossheads. The experimental results for the ring compression test are shown in figs. 3 - 6 by the solid lines (curves 1). The load per unit length of tubing is plotted against the amount of deflection. An additional experiment was carried out in order to measure the uniaxial tensile properties of the workhardening coefficients. The specimens were plates of 1.0 mm thickness and 20 mm gauge length. Their heat treatment was the same as for the foregoing ring compression test. The test was performed by the lnstrontype tension machine at the crosshead speed of 5 mm/min. The stress and plastic strain (o - ep) curve obtained was divided into five regions in order to be fitted by the exponential function, namely o = X ( Y + ep) n . In each region the fitting parameters were calculated by the least squares method. The yield stress of the specimen was ~35.5 kg/mm 2 .
F(Gdetp - HdeaP)2+ G(Hdeap - FdeP) 2 +H(Fderp - GderP)2 X
(FG + G H + H F ) 2
+ d')'taP2 + d'r~P~ + d"/rp211/2 2L
Zircaloy cladding tube is well known as an anisotropic material. For an anisotropic material that possesses three mutually orthogonal planes of symmetry, Hill's criterion is used as the yield criterion. If its radial, tangential and axial directions are the principal axes for the anisotropy in the cylindrical coordinates (r,t,a), Hill's criterion is given in ref. [11] as °=
//2
3 (F+G+H)
[F(°t -- O'a)2 + G(°a -- or)2
-2N-'-]
{ao)
=
[De
I -
---
(1)
where ~-is the equivalent yield stress; suffixes r, t and represent radial, tangential and axial directions, respectively. F, G, H, L, M and N are anisotropy parameters. The equivalent plastic strain increment d~-p related to ~-is defined in the following way: d~-p = V 2 ( F + G + H) 3
(2)
---
-
-= [Dp] {2~e),
(3)
where [De] is the elastic stress-strain matrix and W = do/dep is the work-hardening coefficient dependent on plastic strains. If we express {Ae} and (Ao} by the displacement matrix [B] and [Dp], they become (Ae} = [B] CA/i} and (Ao) = [Dp] (Ae}, respectively. The incremental form of the finite element equations of equilibrium, neglecting body forces and initial strains, are expressed [7] as
=f[Bl T [Dpl [BId V. (aS}
= [KI {a63,
(4)
where the integral is taken over the whole volume of each element. As the radial and tangential directions in a given element, which must be chosen for the principal axes of the anisotropy, vary with continuous deformation of the ring, the finite element equations of equilibrium require to be represented initially in the local coordinates (r,t,a). We transform them to the global coordinate (x,y, z) and sum them over all dements. Let the angular transformation matrix be [T], then the relation between local and global coordinates becomes {A6~) = [T] {A~g}
+ H(ar - 002 + 2(Lr2a + g'ra2r + Nr2t)] 1/2,
"
We assume that Hill's criterion is applicable to the zircaloy tube [12]. The relation between stress increments (Ao} and strain increments {Ae} is represented by [8]
{A/) 3. Calculation
2M
and {AFt} = [T] CA/g),
(5)
where suffixes £ and g represent local and global coordinates, respectively. From eqs. (4) and (5), we obtain {Z~g) = I T ] - l [K] [TI ( a 6 g ) .
(6)
The summation of eq. (6) over the whole elements gives the total euqations of equilibrium. The matrix {ASg} gives the stress increments {Ao~} as {Ao~) = [Dp] [B] [T] {A6g} .
(7)
77
R. Yamada et al. / Application of finite element method to ring compression test Geometry of two dimensional
i-1+ t i ]
t i
inite element
Yc- Yc
Z~Ycl
) t i i
Load
ty~ "~ rr-Ye~ y~-1_y i - - - ~ l z--,-dy~ '= - AtY~b,(i=m)
ZircaLoytube dimension inner diam. 9.30 mrn outer diam. 10.75 mm
~
i sfep number t: temporary value Y: y- coordinate F : nodal force in the
direcfionof y axis
Ay?= at..~yjrn, J=o,b,c FOrced disPlQcement of the node b and c
e>O
x
fes
~,u.,j
,
Fig. 1. Schematic diagram for calculation by the finite element method. or, at and oa involved in eq. (1) are obtained by the summation of eq. (7). We assume two-dimensional plane conditions for the calculation of eq. (6); namely, the plane stress condition (Oa = 0) and the plane strain condition (% = 0). Triangular elements are adopted for finite elements whose grids and fixed conditions are shown in fig. 1. When we apply the finite element approach to the ring compression test, the most important problem is how to handle the contact between the tube and the crosshead of the compression tester. From eq. (4) the incremental method can obviously reduce the nonlinear problem, such as the plastic flow, to the linear problem. The simulation of the ring compression test using the linearity is as follows: eq. (6) is solved to obtain the temporary nodal forces and displacements by applying a certain negative force displacement increment to the constrained nodal points that take a maximum y-value. If y-values of some free nodal points exceed the constrained one at a certain deformation stage, one proportional coefficient (ix) is calculated to equalize them, and then the temporary values are multiplied by it in order to obtain the true values. These free nodal points are dealt with constrained ones from the next stage. If forces of the constrained nodal points become positive at a certain stage, another proportional coefficient (~]) is calculated in order to make the forces zero, and then these constrained nodal points are converted into the free points from the next stage. The computer block dia-
I
dF=t , (i=n) ~' . _ _ ia t I1 Ayj =~'. Ayj, j=a,b,c
I
Free the node c
Fig. 2. Block diagam to flatten nodal points on the tube surface. gram for three nodal points (Yc >Yb > Ya at initial stage) are given in fig. 2. The proportional coefficients a and/3 can exclude the recalculation of the contact problem in that example.
4. Results and discussion 4.1. The effects o f anisotropy on the results o f simulation
The results by which we simulate the ring compression test are strongly dependent on the number of nodal points and triangular elements. These results approach the experimental ones and increase with the number of elements and nodal points. The dashed curves 2 in figs. 3 and 4 show the results of the simulation adopting the van Mises yield criterion under plane strain and plane stress conditions, respectively. It is seen that these results become less than the experimental ones in the initial plastic deformation region. The elastic solution of the ring compression is given analytically [13] as
12e/(bE),
(83
78
R. Yamada et aL /Application o r finite element method to ring compression test
I
I
'
E
[
I
Plane Strain
~:, ~} '~'
Loadl,"
•
-
1") Exp.
Exp [sofropic a,
Plane strain ",c v c-
20 E m
~) H=O.5F,G:O5F. N : 2 OF
J
H=O5F,G=I.0F.N=20F 4 , H:O.5F,G= 15F, N = 2 OF
2O
4-
!
E
lO
A,
/
,
,
~
,,j)
10
d
i
0
2
4
6
Deflection
J
Plane stress
I
'
]
Plane stress
] sotropic
(mm)
range, but deviate at the advanced one, whereas the calculated with smaller G value (i.e. G = 0.5 F) approaches the experimental ones. These results suggest that anisotropy parameters should take different values depending on the degree of plastic deformation. Maki and Ooyama have measured the plastic deformation and reported that the anisotropy of zircaloy-2 changes with the progress of plastic deformation [14]. As they did not reveal values of the anisotropy parameters, we estimated values from their experimental data and simulated using Hill's criterion with the estimated values shown in figs. 3 and 4 (curves 3). The estimated anisotropy parameters do not appear to give good results in the initial plastic region under
i
~, Exp
:2)
6
Fig. 5. Dependence of calculated results on anisotropy parameters under the plane strain condition.
F
E
4 Deflection
where P is the load in the deflected direction, b and T are the length and the thickness of the tube, respectively, and R is the mean radius of the tube. The result obtained by eq. (8) is indicated in figs. 3 and 4 by a dotted line. It is seen in the elastic region that the finite element method results agree well with the analytical elastic solution. The effects of the anisotropy parameters on the results of simulation are examined by varying the values of F, G, H a n d N divided by F. These results, shown in figs. 5 and 6, indicate that if G becomes larger than H and F (i.e. G = 1.5 F) under the plane strain condition, the calculated results agree with the experinaental ones in the initial plastic deformation
A
2
(mm)
Fig. 3. Calculated results under the plane strain condition. (a) Anisotropy parameters estimated from the data of ref. [14]; (b) anisotropy parameters listed in table 1.
I
0
8
:~
'
Exp.
I
]
,~ H:O.5F.G=O.5F.N=2OF
E
(~
H = O.5F, G= 10F, N=2.0F
4)
H = 0 5 F . G= 1.SF, N = 2 0 F
,~4;, b 20
:E D ; L&,®
cc-r
__m t-
E
~ lO
!}
1o
123
3
_J ,
I
2
,
I
4 Deflection
i
I
i
6 ( mm )
Fig. 4. Calculated results under the plane stress condition.
i
0
I
i
2
4 Deflection
I
J
I
6 (mm)
-Fig. 6. Dependence of calculated results on anisotropy parameters under the plane stress condition.
R. Yamada et al. /Application of finite element method to ring compression test
Table 1 Anisotropy parameters used in the present work
~p < 1 × 10 - 5 1 X 10 -5 ~
H
G
N
0.5F 0.5F 0.5F 0.5F 0.5F
1.5F 1.0F 0.5F OAF 0.3F
2.0F 2.0F 2.0F 2.0F 2.0F
b o t h the plane strain and stress conditions. From the comparison of anisotropy parameters' dependence in figs. 5 and 6, we chose the values listed in table 1 where the G value becomes larger than the one estimated from Maki and Ooyama's data in the range
79
below 10 - 3 of % . It is understood from the results in figs. 3 and 4 (curves 4) that in the early stage of ring compression the calculated values under the plane strain condition and in the advanced stage, the calculated values under the plane stress condition, agree with the experimental ones. The reasons for these conditions are considered in the following section. 4.2. The stress and strain distributions
The stress and strain distributions in the flattened tube during the deformation in figs. 3 and 4 (curves 4) are shown in figs. 7 - 1 0 , A and B in the figures indicate the distributions under the plane strain condition whereas C and D indicate those under the plane stress condition. The plastic region at each deflection is dis-
A
Fig. 7, Contour map of er distribution. The unit for st/ain is a fraction. A,B,C and D represent the contour maps at each deflection indicated by the same characters as in figs. 3 and 4. These characters are commonly used in figs. 7-11.
R. Yamada et aL /Application o f finite element method to ring compression test
80
2
~
4
D
"
2
:
~
?
r~ ~
B
CONTOUR 1 ? 2 6 3 4 4 3 S ] 6 O ? ] 8 3 9 4 10 G 11 ?
V~LUE SO00 E + 1 OOOO E + 1 SO00 E + 1 O000 E + 1 SO0@ E + 1 O000 SO00 E + 1 @000 E + 1 SO00 E + 1 0000 E + ] SO00 E + 1
SiOM~-T GEOMETRY
CONTOUR SoO E - 1
INE ~__j
Fig. 8. Contour map of o t distribution. The unit for stress is kg/mm 2.
palyed by Y in fig. 11. It shows that the spread of the plastic area coincides with the increase in the amount of deflection. It is to be noted that the regions inside and outside of the tube near the x and y coordinates are already plastic even at small amount of deflection where the load is linearly dependent on deflection, that is, point A of the deflection. The radial stresses are nearly zero except the compressive stress which occurs at the region of contact with the crosshead, and the small tensile and compressive stresses which occur at the outside and inside of the neutral surface, respectively, at the severe deflection. The er distribution in fig. 7 is remarkably large compared with the or distribution because of the influence of large ot. As the o a in the initial deformation take small values that do not exceed the yield
stress (~35.5 kg/mm2), the tensile and compressive stresses are considered to be cancelled out macroscopic ally in order that the plane strain condition can be fulfilled in that range. In the advanced deformation, however, the increases exceed the yield stress and the tube becomes plastically deformed accompanied by the stress relaxation that is thought to satisfy the plane stress condition. If a criterion for the occurrence of crack is defined as more than 10% tensile strain, then cracks will occur along the radial direction from the outer surface of tube (et > 0.1) around the x coordinate, or along the tangential (Ca > 0.1) or axial (er > 0.1 ) direction from the inner surface around it. The experimental evidence from the ring compression test shows that cracks occurred along the radial direction at the outer
81
R. Yamada et aL / Application or finite element method to ring compression test
C
~
~
~
l
B
@ .
"
\//~ ,1
Fig. 9. Contour map of et distribution. surface of the tube along the x coordinate and were also found to stretch along the axial direction. These were thought to be formed by the accumulation of many small cracks along the tangential direction at the inner surface. It appears that more than 10% tensile strain can be a criterion for the occurrence of cracks. It is worth noticing that the maximum tensile strain which results in the occurrence of a crack is taken by the tangential strain at the outer surface of the tube in the x coordinate of fig. 1. The variance of the tangential strain of the finite element at that surface against the flattening value f i s shown in fig. 12. The computed strain takes a larger value than the experimental result of Ueki et al. [15] at large elonga-
tion. The reason is thought to be that the surface elongation measured by strain gauge gives mean values averaged by the larger area to be detected than the area of the finite element (i.e. the length of the outer surface element is about 0.1 mm).
4. 3. Estimation of yieM stress and maximum elongation The ring compression test cannot directly give the yield stress and maximum uniaxial elongation up to fracture obtained by the uniaxial tension test. Because it requires more difficult specimen fabrication procedures compared with the ring compression test described in the introduction, it is useful if the ring
82
R. Yamada et aL / Application o f finite element method to ring compression test
~St~'TOU£ L J~,lF = . 0 £ 1 ........ Fig. 10. Contour map of oa and ea distributions. compression test is used to estimate the yield and maximum elongation. We investigated the possibility in order to estimate them by using the perfect elastoplastic model (i.e. non-work-hardening in the plastic range) for simplicity. The amounts of loads, where the linearity between deflection and load disappears, are estimated by our program by varying the inputs of yield stresses. The results are shown in fig. 13 where loads are normalized to correspond with eq. (8). It may be said that the mean tensile stress during plastic deformation in the uniaxial tension test can be found by the normalized load at the appearance of nonlinearity in the ring compression test from fig. 13 because it shows fairly good correspondence between the estimation and the foregoing experimental results that give ~2000 kg/mm for the normalized load and
~50 kg/mm 2 for the average work-hardening stress. If the fracture energy is estimated by fO-dgv, the maximum elongation is considered to be calculated as follows. As the deflection at the fracture is known by experiment, the amount of f~d~-p up to it at the fractured element is calculated by our program. It should be noted that there is no difference in choosing an element through the calculation of fod~'p if the stresses and strains are converted basically into the equivalent stress ~, and plastic strain fp. For the fractured element, we chose the outer surface element at the x-coordinate. We estimated the flattening value at the fracture and the maximum equivalent plastic strain and the maximum tangential strain in each yield stress of the perfect elasto-plastic model, by assuming the fractured deformation stage to possess the equiva-
R. Yamada et aL /Application of finite element method to ring compression test
g g g~g~#~gYyg'~9'~y~'~t~F~yyY~yyY
Y~,y* Yy rY rr yY Yr yy y yY 't ~"#
,,-z
~#~Y ~.Y~
~
-
~ . _
~ L £ . ,
-
r.
-
~.
"~7
-
83
' > e
' ,'
~ \\
<~7(.. .
,
•
.
~.,-.~:.
:~"~'Y
f
".f
.
,
, ....,
".".'
....
'~: : : :
A ""
"~N Y
".\\
\\ \
\\
,\ \\
\
\
\
7
OEONETRY
Fig. 11. Yield region of flattened tube.
I
,50
ko c-
~,0
././/o,
-e
,oo
'
I
~ mb-r
/°
I
J°/°'_
/°
'-'
'
5O E ,
20 Flottening vQlue (%}
~0
Fig. 12. The tangential strain at the outer surface of tube vs. flattening value,
Z
[ 40 Yield
i
[
i
60 stress
80 Kg/mm 2
Fig. 13, Normalized load vs. yield stress in the perfect elastoplastic model.
R. Yamada et al. / Application o f finite element method to ring compression test
84
I
I
'
I o
~ao-
80~
~ 60-
60 m
(5) The possibility of estimating the mean stress and maximum elongation for a uniaxial tension test is suggested by the simulation of the ring compression test using the perfect elasto-plastic model.
t~
m
Acknowledgement r~
g~o-
~E
20-
20 5
1
I
i
The authors express their thanks to Mr. A. Morishima and Dr. S. Kawasaki for their encouragement and useful discussion, and also to Mr. M. Ono, CRC Co. Ltd., for programming of contour maps.
I
0 50 70 Flattening value at fracture (%) F'ig. 14. The m a x i m u m t a n g e n t i a l s t r a i n a n d e q u i v a l e n t plastic strain o f the finite e l e m e n t at the f r a c t u r e o f t u b e in the perfect elastic-plastic m o d e l .
lent fracture energy of the former fod~p. This estimation is shown in fig. 14. We think that the maximum tangential strain at the outer surface of the tube, which may have close connection to the maximum elongation of the uniaxial tension test, is ably obtained by fig. 14, knowing experimentally the flattening value of the tube at the occurrence of crack.
5. Conclusion
Nomenclature b d D [Del [Dpl E
f
length of tube displacement matrix d e f l e c t i o n of t u b e initial f a b r i c a t e d d i a m e t e r elastic s t r e s s - s t r a i n m a t r i x plastic s t r e s s - s t r a i n m a t r i x Young's modulus = flattening value o f t u b e
F,G,H,
L,M,N h P r, t, a R
x, y, z
We attempted to simulate qualitatively the ring compression test by developing a computer program based on the two-dimensional elasto-plastic finite element method. The results obtained are as follows. (1) The isotropic yield criterion (i.e. the yon Mises yield criterion) gives a little difference, whereas use of the plastic strain dependence from Hill's anisotropy parameters gives a good fit with experiment. (2) In the early stage of ring deformation the results calculated under the plane strain condition and in the late stage the results under the plane stress condition agree with experiment. It is considered thai the axial stress strongly influences these results. (3) A crack can be expected to occur in a location where there is more than 10% tensile strain. (4) A small portion of tube is plastically deformed even when linearity exists between load and deflection.
= = = = = = =
x,r,n T
[T] W
(} LJ
= a n i s o t r o p y p a r a m e t e r s of Hill's y i e l d c r i t e r i o n = d i s t a n c e b e t w e e n the t w o c r o s s h e a d s of the compression testing machine = load w e i g h t e d to t u b e = radial, t a n g e n t i a l , a x i a l d i r e c t i o n s in the cylindrical c o o r d i n a t e s = m e a n r a d i u s of t u b e = Cartesian coordinates = f i t t i n g p a r a m e t e r s in the e x p o n e n t i a l f u n c t i o n = t h i c k n e s s of t u b e = angular t r a n s f o r m a t i o n m a t r i x = work-hardening coefficients = c o l u m n vector = r o w vector
Greek symbols = n o r m a l stress and s t r a i n of ith c o m p o n e n t oi, ei = shear stress a n d strain rii, 7# = e q u i v a l e n t stress and strain = nodal displacement increments = n o d a l force i n c r e m e n t s = stress i n c r e m e n t s = strain i n c r e m e n t s = proportional coefficients Subscripts and superscripts e = elastic p = plastic = local c o o r d i n a t e s g = global coordinates.
R. Yamada et al. / Application o f finite element method to ring compression test
References [1] D.O. Pickman, Nucl. Eng. Des. 21 (1972) 212. [2] A.A. Bauer, L.M. Lowry and J.S. Perrin, BMI-NUREG1948 (1976). [3] S. Kawasaki, T. Furuta and M. Hashimoto, JAERI-M 6181 (1975) (in Japanese). [4] D.O. Hobson and P.L. Rittenhouse, Oak Ridge Report, ORNL4758 (1972). [5 ] D.H. Bradhurst and P.M. Heuer, J. Nucl. Mater. 55 (1975) 311. [6] R. Hill, The Mathematical Theory of Plasticity (Oxford University Press, 1950) p. 38. [7] Y. Yamada, N. Yoshimura and T. Sakurai, Int. J. Mech. Sci. 10 (1968) 343.
85
[8] R.H. GaUagher, Y. Yamada and J.T. Oden, Recent Advances in Matrix Method of Structural Analysis and Design (The University of Alabama Press, 1971) p. 283. [9] O.C. Zienkiewicz, Finite Element Method in Engineering Science (McGraw-Hill Book Co., 1971) p. 435. [10] Y. Yamada, Plasticity and Viscoelasticity (Baihukan, (1972) p. 173 (in Japanese). [11] R. HiU, ref. [6], p. 317. [12] Y. Miyamoto, Y. Komatsu, T. Kakuma and N. Nagai, J. Nucl. Mater. 61 (1976) 53. [13] S. Timoshenko, Strength of Materials, Part 1 (D. Van Nostrand Co., 1955) p. 381. [14] H. Maki and M. Ooyama, J. Nucl. Sci. Technol. 13 (1976) 43. [15] I. Ueki, K; Fukai and H. Ichikawa, J. At. Energy Soc. Jpn. 18 (1976) 591 (in Japanese).