Compurers dr Srructures Vol. Printed in Great Britain.
17, No. 4, PP. 555-561,
0045-7~9/83/10055~7$03.00/0 Q 1983 Pergamon PressLtd.
1983
A FINITE ELEMENT PROCEDURE TO TIME DEPENDENT CONTACT ANALYSIS E. ZOLTI Interatom-Intemationale
Atomreaktorbau
GmbH Postfach, West Germany
D-5060
Berg&h
Gladbach
1,
(Received 7 July 1982; received for publication 21 October 1982)
Abstract-The subject of this paper is a finite element approach to solve contact problems involving contact-time dependence and interactions and frictional conditions. The procedure find8 application in elastic or inelastic analysis of structural systems with two or more parts separated by small gaps and subjected to cyclic static loadings, as encountered, e.g. in nuclear reactor core components. The procedure is based on the load increment method and the introduction of orthotropic gap elements. Their stiffnesses both in the directions normal and tangential to the opposing surfaces assume an infinitesimal or a finite value according to the open or close configuration and to the frictional conditions. At every time step the locations of the contacts are at first identified through the non-positive values of the Jakobian in the gap elements. Then the new global stiffness is determined with the contribution of the stiff elements at the contact regions limited to the present time step to account for the preceding free relative motion. Two elementary illustrative analyses show the validity of the procedure in determining the time points of the gap closure and opening and in accounting for contact interactions. Finally an application to a real contact problem in the reflector subassembly of a nuclear reactor core is given.
1. INTRODUCTION When two or more members of a structural system are separated by small gaps and subjected to variable external loads, they can come into contact and disconnect repetitively. The gap closure could result in a mechanical interaction of primary importance for the structural response of the system, as high peak stresses and strains could arise at the contact locations and the associated large fatigue and creep damage could shorten remarkably the structural lifetime of the component in comparison to the case without contacts. Such a kind of problems are encountered for instance in nuclear reactor core components and typically in fuel pins, where the extent of the pelletcladding gap influences considerably both the thermal and the mechanical loads of the cladding. Several works on the solution of contact problems have been published (e.g.[l+). However, they are mostly limited to the elastic material behavior, or to frictionless contacts or to continuous and plane contact surface. Attempting a step to a more realistic treatment of the reactor core situations, the concern was therefore focussed on the development of a procedure based on already available inelastic finite element program and in which the central question was the identification of the load histogram time points of the gap closures and openings and of their location, accounting also for possible interactions in case of multiple contacts. The determination of the strains and stresses in the closed configuration, i.e. once the contact has been identified and the actual displacements have been obtained ensues straightforward, as it is a matter of a usual finite element calculation. The finite element program in which the procedure is incorporated is IAFETIN[S], which using 6 node isoparametric triangular elements permits thermal and inelastic mechanical analyses of 2 dimensional structures. Since the code IAFETIN has no dynamic capability
and dynamic impact situations are not of primary importance for nuclear reactor core components, the procedure is at present limited to static problems. The procedure is briefly outlined in the next section.
2. ANALYSISPROCEDURE The main points of the IAFETIN-solution of contact problems consist in the recognition of possible material overlapping in the structural members considered free to deform at the opposing surfaces under the action of the external load, and, if the case, in the change of the global stiffness to determine the real displacements. More in detail, the solution is obtained in the following way (Fig. 1); (a) At the finite element schematization of the geometry also the gap regions where a contact is expected are provided with finite elements. The materials of these gap elements can be real materials for the thermal analysis, e.g. sodium, gas, but represent fictive orthotropic materials for the mechanical analysis. For the E-Moduli E, and En in the direction tangential and normal to the opposing surfaces two groups of values are specified in input to be used in the open or closed configuration. Frictional conditions can be simulated through adequate E,, as, e.g. E, = lo- lo N/mm’ or E, = E, for the limit cases of frictionless sliding or contact without sliding respectively. (b) In order to better handle the time dependence of the problem with adequate time steps the incremental procedure normally used for inelastic analysis is initiated, even if the structure is supposed to behave only elastically (c) For every time increment A$ a first value of the nodal displacements (6 1; corresponding to the present total loads cr), and to the contact-free configuration is 5.55
4
Fig. 1. Flow of the time step computations
I
total loads fj - stiffness matrix k
-
in the IAFETIN
procedure
to contact analysis.
A finite element procedure to time dependent contact analysis
a)
t, = 0
bl
t
= tj_1 _
Cl
t
=
551
I tj
Fig. 2. Time changes of the gap elements.
calculated through the solution of the algebraic system
[Kl’{61j= cr>, where [a is the global stiffness matrix obtained with the infinitesimal E-Modulus of the fictive gap materials, i.e. practically with a zero-contribution of the gap element stiffness. Thanks to the very small (but not zero) values of the gap elements properties, no difficulty arises during the inversion of the matrix[K]‘. (d) On the basis of the displacements {S}, the nodal point coordinates are updated and then the compatibility of every gap element is checked through the Jakobian J at the integration points. The non-positive value of J indicates the element has been destroyed, because the boundaries of the deformed gap element cross and denote the occuring of the contact at the present time interval Atj. Figure 2(a-c) illustrate schematically the reduction of gap from S, at t,, = 0 to sj_ , at tj_ , and the overlapping at the time t!. (e) If no overlapping has been identified (1 > 0 everywhere) the calculation continues as usual (see point (g)). Otherwise, every element, where at least in one integration point J $0, is considered to be entirely in contact. The E-Moduli of the orthotropic material in this element are switched to the second values specified for the closed configuration and the displacement calculation is rerun with a new stiffness matrix[K]. However, the simultaneous use of the total load (r}, and of the new stiffness matrix[K] in the (1) would mean to consider the new “stiff’ contact elements as acting since the time t = 0. In reality until the end of the previous time step, i.e. the time t,_ ,, the gap elements practically have had zero stiffness and so have not
hindered the relative motion of the opposing surfaces (corresponding to the gap reduction S, - sj_ , in Fig. 2). In order to limit the effects of the additional stiffnesses of the detected contact elements to the present time step A$ an incremental approach is used, based on the old displacements {Sz _ , (at the time fj_ ,). Only the increments of the displacements {AcS}~ m the present time step A$ are determined at first with the new [K] through the solution of the system
(2) where now {Afjj is the increment of the load vector corresponding to the increment of the external loads and of the initial strains in A$. These increments A6), are then added to the old total displacements a}, _ , to get the total displacements (~3)~at the time tj
{a},={a},-1+ {Ad}, In the first time step of course the two approaches (total vs incremental) coincide. (f) A new compatibility check through the Jakobian values is performed for all the gap elements, permitting new contact regions due to the changed global stiffness to be detected. If the case, the calculations of the point (e) are repeated until nowhere element overlapping results. (g) Once the final values of the displacements are available the computation continues, as usual, to the determination of the total strains, the stresses, the increments and the total values of the inelastic strains. (h) During the inelastic analysis iterations in the
558
E.
analysis with temperature-dispiacemeut interactions, the temperature dist~butions not yet available for the next time points are determined using new thermal properties for the “stiff’ elements at the contact regions, similarly to the mechanical analysis.
E=16,5i0' N/mm' 1= 0.3
al
Kric:hgesetr:a (p &,61940' e1.986.1
Test example The following elementary test analysis, ap preachable also with simple hand calculations, shows the capability of the IAFETIN-procedure in the determination of the contact time dependence. The bar of Fig. 3(a) is axially loaded and after an elongation equal to the small gap s comes into contact with the rigid obstacle H. The axial load varies with time as indicated in the load histogramm OABCDE of Fig. 3(b). The gap closes at the times ti and ts after elongation due to elastic or creep strains respectively and opens at the time t, during unloading.
?I_
-H
%LTr
-I-
Contact interactions In case of multiple contacts gap closure or opening can be due not only to external loads but also to contacts at other locations. Such a situation of contact interaction occurs, e.g. in the illustrative example of the cantilever beam of Fig. 4 whose deflections are limited to the small gaps between the end of the branches A, B, C, D and the rigid obstacle. Although the thermal expansion would suffice for some gaps to close the actual open or closed configuration depends on the bending deflections associated with the contacts at the other locations. In the IAFETIN-approach this problem is handled with a set of further iterations in the time step, in which material 0verIapping at some locations is interpreted as definitive local gap closure only after the deformations due to the sure contacts have occurred. The locations of the sure contacts and of those ones to be iteratively confirmed are specified in input by the user on the basis of the residual gap and load increment at the present time step with help of restart runs if necessary. In the ilIustrative example of Fig. 4 above the IAFETIN-results show that irrespective of the thermal expansion the gap closes only at the location A at the time r’ (Fig. 4b) and A and D at the time r” (Fig. 4c).
df
5+--+-f
102
!*
203
Fig. 3. Simple test calculation of the gap closure and opening.
time step At, the procedure of the load and displaeement increment (see point (e)) is applied everytime at least one contact is present. Thereby the inelastic strain increments of the previous iteration are used to build the nodal point forces due to initial strains and there is no need to invert again the global stiIIness matrix of the previous iteration as long as the contact gap elements do not change. (i) After the inelastic calculations of a time increment A$ have been concluded, the next step At,, , is commenced. Before that, however, in case of
3. ANALYSIS EXAMPLE
As example of application to a real contact problem the case of a reflector subassembly of a fast breeder reactor core is now discussed. (See schema of Fig. 5). Due to the neutron flux and temperature gradients non-uniform swelling and thermal deformations result, causing the refIe&or block and the outer hexagonal tube to come into contact after a certain operation time. The analysis aims at the determination of the zones and of the time points of the contacts. Since the stresses are irrevelant in this analysis (swelling depends only on neutron fluenee and temperature), the reflector subassembly could be schematized with a 2-~mensional model, which responses with the same deformations and therefore reproduces the same contact situation. In this model the inner block and the outer tube are substituted with an inner and an outer (wrapper) beam respectively, having the same bending stiffness (Fig. 6).
A finite element procedure to time dependent contact analysis
559
a) t I 0
b) t = t’
ALh = aL,AT s
s,
Ab = abAT
*
%
A4 ii a&Al
+
SC
Ab = aL,,Abf *
Fig. 4. Illustrative example of contact interaction problem.
t
!
2
6-
-t
-.L *
5
_i-
z Fig. 5. Geometry schema, loading and boundary conditions of the reflector subassembly.
b
a) Location A
location A mid
1:
i,. _I_.__,_~ PI
t=o
-
b)
t = 14175 h
Location B
Fig. 6. Two dimensional schematization of the reflector subassembly for the IAFETIN contact analysis
The E-moduh of the material 1 replacing the tube steel I .498 1 and of material 2 replacing the block steel 1.4948are therefore obtained from the conditions
and
where & and JI are the bending moments of inertia of the outer and of the inner beams equivalent to the outer tube and to the inner block respectively. Figure 5 shows also the boundary and the loading conditions, i.e. the distribution in x and z directions of the temperature and of the neutron flux, In the finite element discretization a fine mesh has been chosen in the regions where the contacts were expected (A and 3 in Fig. 6). For the fictive orthotropic materials of the gap elements (mat. 3) the following E-Modulus values were specified: E, = E, = 10-'N/mm2 in the open
Fig, 7. Plotter of the initial and deformed geometry at the two contact locations.
configuration and E,;= Emar,Z and Ei= 10m9N/mm2 after gap closure with E: corresponding to the assumption of frictionless sliding at the contact focations, The computation showed that after a first gap closure at the lower end of the block (zone B of Fig. 6) an additional contact at the region A of Fig. 6, i.e. near the mid core, occurs, resulting in a further bending loading of the outer tube. This twofoid contact ~onfiguratjon which is actually of concern from the structural point of view resulted to begin the 578th and the 590th day of the nearly 1I years full power operation. Fig. 7(a, b) show the plotter of deformed geometries at the 2 contact locations.
A finite element procedure to time dependent contact analysis 4. SUMMARY
AND CONCLUSIONS
The approach presented permits the detailed analysis of the structural response of systems with two or more parts separated by small gaps and subjected to cyclic static loadings. In particular the local peak stresses and strains can be accurately determined and used, e.g. for a realistic evaluation of the fatigue and creep damage and therefore of the structural lifetime. The degree of the solution accuracy can be established by the user choosing adequately the finite element mesh at the contact regions and the time steps preceding gap closure or opening. However, the solution accuracy is directly related to the analysis cost. In fact every time a gap closes or opens a new global stiffness matrix must be built and inverted in “contact’‘-iterations, and the corresponding computer time adds to the time for the inversion anyway required at the beginning of every time step. Thus the total computer time can be considerable and easily achieve exorbitant values, especially if large structures and strong inelastic material behaviour are involved. The engineering judgement will help in the cases of such complex problems to find the optimal compromise between the accuracy and the economy of the solution.
561
Acknowledgemenrs-This work was performed in the frame of the R&D activities of the German fast breeder program, financially supported by the BMFT in Bonn. The author wishes to acknowledge the useful discussions and suggestions provided by Dr. H. Tiibbe of the department of core component design and development at INTERATOM.
REFERENCES
1. J. T. Stadter and R. 0. Weiss, Analysis of contact through finite element gaps. Cornput. Structures 10, 867-873 (1979). 2. F. F. Mahmoud et al. A direct automated procedure for frictionless contact problems. Inr. J. Num. Mefh. Engng 18, 245-257 (1982). 3. N. Okamoto and M. Nakazawa, Finite element contact analysis with various frictional conditions. Inr. J. Num. Mefh. Engng 14, 337-357 (1979). 4. B. Fredriksson, Finite element solution of surface nonlinearities in structural mechanics with special emphasis to contact and fracture mechanics problems. Compur. Structures 6, 281-290 (1976). 5. E. Zolti, Localized thermostructural analysis of LMFBR core components through the computer code IAFETIN. 61h SMiRT Co& Paris, Aug. 1981.