Response of fast-reactor core subassemblies to pressure transients

Nuclear Engineering and Design 68 ( 198I) 153-174 North-Holland Publishing Company

153

RESPONSE OF FAST-REACTOR CORE SUBASSEMBLIES TO PRESSURE TRANSIENTS J. DONEA, S. GIULIANI and J.P. H A L L E U X Applied Mechanics Division, Joint Research Centre of the Commission of the European Communities, Ispra Establishment, Italy

P. GRANET, P. H A M O N a n d P. PERMEZEL Commissariat h rEnergie A tomique, Dbpartement des Rbacteurs h Neutrons Rapides, C.E.N. Cadarache, B.P. 1, F-13115 St. Paul-lez-Duranee, France

Received 24 August 1981

An Arbitrary Lagrangian-Eulerian (ALE) finite element method for analyzing the transient, nonlinear fluid-structure interaction in fast-reactor core subassemblies is presented. The method combines the basic attributes of the finite element technique - namely, ease in modeling complex geometries and in mixing fluid elements with structural elements - and the flexibility in moving the fluid mesh offered by the ALE description. The result of this combination is a very versatile modeling technique which permits accommodation of large fluid and structure displacements and logically simple, but accurate fluid-structure coupling. Numerical results are presented to illustrate the effectivenessof the proposed method. These include the response of a single hexagonal duct to internal, static or dynamic, pressure loading for which numerical predictions are compared to experimental data, and applications to clustered hexcans.

1. Introduction The ability to predict the nonlinear response of a cluster of core subassemblies exposed to transient dynamic loading is of considerable importance in the safety analysis of liquid metal fast breeder reactors. The numerical analysis of problems of this type requires very versatile modeling techniques capable of treating complex fluid-structure systems undergoing large displacements. The purpose of the present paper is to present an Arbitrary Lagrangian-Eulerian (ALE) finite element method for analyzing the transient, nonlinear fluid-structure interaction in fast-reactor core subassemblies. The method combines the basic attributes of the finite element technique - namely, ease in modeling complex geometries and mixing fluid elements with structural elements - and the flexibility in moving the fluid mesh offered by the ALE kinematical description. The result of this combination is a very versatile modeling technique which permits accommodation of large fluid and structure displacements and a logically simple, but accurate fluid-structure coupling. The paper is organized as follows. In section 2, we briefly recall the basic conservation equations for mass, momentum and energy in an ALE reference frame and

derive the associated weak variational forms which form the basis for spatial discretization using finite element modeling. Section 3 is devoted to the presentation of the semi-discrete equations governing the flow of an inviscid, compressible fluid; 4-noded isoparameteric elements are employed to model the hydrodynamic domain. Numerical time integration of the semi-discrete equations is discussed in section 4. An algorithm which automatically prescribes the displacements of the fluid mesh is described in section 5. Such an automatic algorithm enables us to avoid excessive distortions of the computational grid and yet allows the fluid boundaries, including deforming structures, to be tracked with the necessary accuracy. Section6 is concerned with the modeling of the structural part of the fluid-structure problem and includes a short description of finite element models. Fluid-structure coupling is discussed in section 7 with the emphasis being placed on the ease and accuracy with which the coupling may be achieved. when the fluid motion is described in the ALE formulation. Special attention is given to problems which arise when the interface between an inviscid fluid and a structure possesses sharp corners. Section 8 presents

0 0 2 9 - 5 4 9 3 / 8 1 / 0 0 0 0 - 0 0 0 0 / $ 0 2 . 7 5 © 1981 N o r t h - H o l l a n d

154

J. l)onea et al. , Respon,se q[ filst reactor core .~u/~assemhlies

material modeling including the effect of fast transient loading. Section9 contains purely structural problems related to code validation exercises and featuring both dynamic and quasi-static tests. Finally, representative calculations related to the simulation of explosive accidents in various configurations of fast-reactor core subassemblies are presented in section 10 as an illustration of the effectiveness of the proposed modeling procedures.

0 0 a t ( p e J ) - = J ~ x I ( p e ( , , ) -- v,))

+J(ovjb,_

The arbitrary Lagrangian-Eulerian formulation was originally developed in finite difference formats, see e.g. ref. [1], in an effort aiming at combining to a certain extent the best features of both the Lagrangian and the Eulerian approaches. ALE finite element methods have been reported more recently (see ref. [2] for a survey), but have already been shown [3,4] to have considerable potential for application to a wide variety of fluidstructure problems encountered in reactor safety analysis. The ALE formulation has no basic dependence on particles and treats the computational mesh as a reference frame which may be moving with an arbitrary velocity w in the laboratory system. When w = v, where v is the particle velocity, the classical Lagrangian viewpoint is obtained; when w = 0 the reference frame is fixed in space and this corresponds to the classical Eulerian viewpoint. When w v~ v #= 0 the reference frame moves in space at a velocity w which is different both from the particle velocity v and from 0 and such a reference frame is called arbitrary Lagrangian-Eulerian. 2.1. Basic conservation laws in the A L E description

As shown in refs. [2,4], the ALE form of any differential law may be obtained from the corresponding law expressed in Eulerian description. In particular, the ALE differential forms of the basic conservation laws for mass, momentum and energy are obtained as follows for the case of an inviscid, compressible fluid: a

o

a t (p J ) = J~-~xj (#(w/ - v , ) )

~7 (p'~'J)

=j~

(mass),

(2.1)

(total energy). (2.3)

0

0

(piJ) :J~,x, (pi(w, ~,)) -

ave -JP ~x,

2. ALE kinematical description of the fluid domain

a (pvj))

(internal energy).

(2.4)

In these expressions, p is the fluid density, p the pressure, v the fluid velocity, w the reference frame arbitrary velocity, b a body force, and e the total specific energy defined by e = ½ v 2 + i,

(2.5)

where i is the specific internal energy. J represent the Jacobian: IOx, I J=[~a/ ,

(2.6)

for the transformation between the mixed coordinates x (identifying the position vector of any point in the ALE reference frame) and the material coordinates a. The Jacobian J provides a mathematical link between the current volume element d V (function of the mixed coordinates) in the reference frame and the associated volume element dV0 (function of the material coordinates) in the initial configuration: d V = J d V 0.

(2.7)

One may show that the time rate of change of the Jacobian determinant is given by aJ =J at

--

v.

,,,.

(2.8)

2.2. Weak variational forms

To provide the basis for the derivation of spatially discrete models of the ALE conservation equations using finite element modeling, we shall now develop weak variational forms of eqs. (2.1)-(2.4). 2.2.1. Mass and energy equations In view of relation (2.8) the ALE differential form (2.1) of the mass equation may be rewritten as

a~, (p,), (wj - v,)) (momentum), (2.2)

ap _ ( , , - v ) . a-7-

vo-

o v-,~.

(2.9)

J. Donea et al. / Response of fast-reactor ('ore subassemblies Multiplying this equation by an arbitrary weighting function ~ and integrating over a control volume V(t), we obtain the following weak form of the ALE mass conservation equation:

fv (t)

dt

at

f

JV(t)

8v,o=-dV=f & , , o ( ~ - , , ) . v v , dV Ot

JV(t)

+f

dV= Jv(t) f

w-v). v

oi dV (2.11)

2.2.2. Momentum equation In order to isolate the acceleration term in the ALE differential form (2.2) of momentum conservation, we rewrite this equation as

pJ~-+v,~(pJ)

=J ~-~j(Pvi(wj°oj))

and make use of the mass equation (2.1) to reduce (2.12) to the simple form

~v, =p(,~-v).vv; +pb, Ox," Op p~-f

(2.15)

3. ALE finite element model of the fluid domain

The Galerkin formulation of the method of weighted residuals is employed to generate spatially discrete models of the ALE conservation equations. The modeling of the fluid domain will be based on 4-noded isoparametric quadrilateral elements. The element geometry, the fluid velocity v and the arbitrary grid velocity w are interpolated over a typical finite element by bilinear shape functons. The fluid density, specific internal energy and pressure are assumed to be elementwise uniform. Under this hypothesis, the weighting function q, in the weak forms (2.10) and (2.11) of the mass and internal energy equations may be taken as unity and the local gradients of p and pi are zero. It follows that conservation of mass and internal energy is expressed at the element level by the following integral statements:

~t fv~(opdV=~s~(t) p ( w - v ) ' n d S

(mass), (3.1)

(2.13)

The required weak variational form of momentum conservation is obtained by operating on eq. (2.13) and on the associated stress boundary relation

r,. = -,,~ 8~.,p,

p~O-~-(6vi)dV

Svft) °xi

+ ~s(t)SviTi dS.

(2.10)

- fv(t)4~(pi+p) v . v d V

8v,ob,dV+ f

Jv(t)

In a similar fashion, the weak form of the internal energy equation (2.4) is obtained in the form

fv it)

desired weak form of momentum conservation:

dv=fv,, (w-v)'vodV -- fv(t)~p ~7 "vdV.

155

(2.14)

where T, represents the components of prescribed boundary loads (force per unit area), nj are the direction cosines of the unit outward normal to the boundary, 6j; is the Kronecker delta and p the fluid pressure. To this aim, we multiply eq. (2.13) by an arbitrary admissible variation 6v~ of the fluid velocity and integrate over a control volume V(t). Similarly, we multiply eq. (2.14) by 8v~ and integrate along the bounding surface S(t). We then add the two integral equations and obtain, after integration by parts of the pressure term, the

- f

p V • v dV

J V "(t)

(internal energy)

(3 2) According to the Galerkin method, the velocity variation in the weak form (2.15) of momentum conservation is defined as

8 v, = N, (x) 6v,;,

(3.3)

where the N t's are the bilinear velocity shape functions, Independent variations of all the velocity degrees of freedom in the fluid mesh are introduced in turn into the variational form (2.15) and by invoking the arbitrariness of such variations, the discrete analog to the momentum equation is obtained The semi-discrete equation governing a typical velocity node I is found in

J. Donea et al. / Response oJ fast-reactor core sut, assembhe~

156

the form

+~fv'Nt(pb'-~x~)dV-~ where ~

sL, NtT'dS'

(3.4)

vanish in the Lagrangian formulation, i.e. when , ..... v. It is this reduction that leads one to first solve the Lagrangian conservation equations in the first step of the calculational cycle, phase I, and then to add the transport contributions as a separate step in the last phase, phase III. The middle step, making up phase II, adds an implicit pressure calculation that permits solutions to be obtained at all flow speeds.

indicates a summation over all the elements

e

to which node I belongs, and the repetition of subscript J implies a summation over the four velocity nodes in element e. The governing momentum equations for the whole fluid mesh form, as shown by eq. (3.4), a set of ordinary differential equations in time which may be written in condensed form as [ M ] { d } = {Ft} + {Fb} + (Fp) + {/~},

4.1. Phase I. Explicit Lagrangian calculation This step calculates the Lagrangian nodal velocities resulting from the previous cycle's pressure and body forces. From eq. (3.5) with w = v in eq. (3.4) and using lumped masses, we obtain

(3.5)

where [M] is the global mass matrix; (t~} the total nodal acceleration vector; {Ft) global nodal loads induced by transport of momentum; (Fb) global nodal loads due to body forces pb; {Fp} the nodal loads induced by the fluid pressure p; and (if} accounts for the externally applied loads T. The integrals involved in the semi-discrete conservation equations (3.1), (3.2) and (3.4), (3.5) are evaluated either analytically (surface integrals) or numerically (volume integrals) using a 2 × 2 Gaussian quadrature rule. Full details concerning the evaluation of the above integrals may be found in refs. [2,4] where the stabilizing effect of the introduction of some selective 'upwinding' into the various transport terms is also discussed.

where the tilde indicates an intermediate time for this phase of cycle n + 1, and the superscript n indicates previous cycle values. The specific internal energy i is advanced to an intermediate time from the integral equation (3.2) with

i = i N- - - ~ ) p

V

dV.

(4.2)

In this expression, M e is the element mass and a mean velocity is used, as suggested by Pracht [6], in order to conserve the total energy exactly in the finite element formulation.

4.2. Phase II." Implicit Lagrangian calculation 4. Numerical time integration Since the partial differential equations governing unsteady inviscid, compressible flows are hyperbolic, time-centered schemes are in principle required to ensure numerical stability of an explicit time integration procedure. However, if an upwind discretization of the transport terms is employed, a numerical diffusion is introduced into the semi-discrete conservation equations and stability may be achieved on using a simple first-order explicit time integration. An interesting computational scheme based on firstorder explicit time integration and valid for both compressible and incompressible flows has been proposed by Hirt et al. [1] and Stein et al. [5]. In deriving such a scheme, one notes that the transport terms in the semidiscrete conservation equations (3.1), (3.2) and (3.4)

The explicit Lagrangian phase detailed above yields tilde values for the fluid velocities and specific internal energies, while the corresponding densities and pressures have yet to be calculated before the calculation enters the third (transport) phase to be described later. If a purely explicit procedure is desired, one may obtain the updated densities and pressures immediately from the equation of conservation of mass and the equation of state of the fluid, viz, pL(l + A t V "v L) = p " ,

(4.3)

pL =f(pL iL),

(4.4)

where we have taken v L = t3, i L = /-and used a timeadvanced discrete form of eq. (3. I) with w = v. In cases (especially those involving low Mach numbers) where this procedure is found to be too noisy, the

J. Donea et al. / Response of fast-reactor core subassemblies following implicit phase may be employed. The object is to solve eqs. (4.3) and (4.4) simultaneously with ( v L) = {v"} + A t { v ( p e ) } M" and

iL=i n

At t"

-~-~Jvp

L

V

(Vn~V L - - )

(4.5)

dV.

(4.6)

Notice that advanced values (superscript L) of all the fluid variables have been used in the Lagrangian conservation equations (except in eq. (4.6) which employs an average velocity value), but that the masses and integration volumes are at time level n, a choice made for computational convenience. Eqs. (4.3)-(4.6) are nonlinear relations for the four quantities OL, pL, 19e and iL, and we solve them using a pseudo-Newton scheme (see, e.g., Pracht and Brackbill [7]). In this scheme, a quantity H is defined which is reduced iteratively to zero for each element:

H=p-f(p,i).

157

evaluated. In phase III we now add the contributions from these terms to the time-advanced level L values. We assume that the mesh velocities w have in some way been specified; how this may be automatically accomplished is described in section 5. The change in volume of the element from time level n to n + l is: V "+ f - V n • AtfsWin idS.

(4.9)

From eq. (3.1) the updated element mass is

M~'+' : M e n 4-At/sePL(Wi--v L) n i d S

,(4.10)

and the new density results from

p"+' = M e + ' / V "+'.

(4.11)

From eq. (3.2), the contribution from the transport terms to the specific internal energy yields

(4.7)

Choosing some arbitrary perturbation ( A p ) e from the pressure at level n in element e, one calculates in tum (Av) e, (Ai) e, (A f ) e and ( A H ) e and so numerically evaluates the derivative (d H / d p ) e for each element. The iteration proper is then entered, with starting values p n On, 6, i. H is calculated from eqs. (4.4) and (4.7) and the new estimated time-advanced pressure determined from

i , + l - M;+, M~' i L + MAt ; +' /S pLiL(w,--V~) n, dS(4.12) and the transport of momentum contribution to the velocity components at node I, from eq. (3.4), is

~i,! --t~),l +

NIoL(wi--V i ) ~ x i d V .

Finally, the new pressure is calculated from the equation of state p " + ' = f ( p " + l , in+l).

where ~0 is a relaxation factor with value close to unity. If the change in p is below a pre-determined tolerance level for every element, then the iteration is complete and the calculation proceeds to phase III. Otherwise, the new pressure values are used to calculate in turn the velocities, densities and internal energies from eqs. (4.5), (4.3) and (4.6) respectively and the whole procedure is repeated. Since the basic objective of the above implicit Lagrangian phase is to obtain new velocities that have been accelerated by pressure forces at the advanced time level, this phase eliminates the usual numerical stability condition that limits sound waves to travel no further than one element per time step.

4.3. Phase III: Rezone or convective-flux calculation Phases I and II were assumed Lagrangian, that is, the contributions from the transport terms were not

(4.13)

(4.14)

This completes the steps contained in a cycle of the time integration procedure.

5. Automatic rezoning The freedom in moving the fluid mesh offered by the arbitrary Lagrangian-Eulerian formulation is very attractive. However, it can be overshadowed by the burden of specifying grid velocities well suited to the particular problem. As a consequence, the practical implementation of the ALE description, and in particular that of the transport phase of the calculation cycle, requires that an automatic mesh displacement prescription algorithm be supplied. The purpose of this algorithm is to provide the ALE computer code with a capability for automatic and continuous rezoning that conserves, as far as possible, the regularity of the computational mesh

15g

J. Donea et al. / Response of fast-reactor core suhassemhlie.~

and yet allows the fluid boundaries including deforming structures to be tracked with the accuracy characteristic of Lagrangian methods. Rezoning techniques are almost entirely based on heuristic developments and their acceptance rests on their performances in practical applications. The particular technique that will be described here has been implemented into EURDYN-1M, an ALE finite element program for transient dynamic fluid-structure interaction [8], and successful applications have been made to problems in reactor safety analysis. To achieve the desired continuous and automatic rezone, the components of the grid velocity at a typical ALE node I are computed at each time step by the following relationship: t+,at 1 ~ t -~--Wj, t wl'i N J

0.1 1 t 1 8t - 8t + $ 7 ~ ~ L , . , ~ ~ "'i '" j

Ltls

'

(5.1)

where N indicates the number of nodes connected to node I via element sides and diagonals, LIj is the current distance between node I and connected node J, and 8 represents total nodal displacements. The above formula was arrived at by trial and error on the following grounds. The basic idea was to give node I a current grid velocity equal to the mean of the grid velocities of the neighbouring nodes at the previous time step. However, this was occasionally found insufficient to prevent an excessive squeezing of some elements, especially in the vicinity of fluid-fluid interfaces, and led us to introduce a corrective term which enhances the grid velocity when the distance between adjacent nodes tends to become too short. The components of the grid velocity computed by relation (5.1) are then limited by the following empirical inequalities which are dictated by numerical stability and accuracy conditions: I - - y < ~ wI''<<- 1 + y .

(5.2)

[l,i

In this expression, v/, ~ are the components of the fluid velocity and y is a numerical coefficient (0 ~
they are constrained to move with the fluid by setting wl, , = vz. r At this point, the assignment of grid velocities is completed and the nodal displacements may be updated by explicit time integration. It is important to note that the presence in the ALE conservation equations of transport terms accounting for arbitrary displacements of the mesh nodes automatically guarantees conservation of mass, momentum and energy during the continuous rezoning procedure. This represents a significant advantage over the rezoning techniques employed to remap the mesh at fixed times in Lagrangian programs.

6. D e s c r i p t i o n

of structural e l e m e n t s

In the analysis of complicated coupled fluidstructure problems as those pertaining to reactor safety investigations, the prime goal in the search for a solution scheme is to define a numerical procedure ensuring natural continuity between and great flexibility inside the various physical simulations. A numerical concept providing these features is the Finite Element Method (FEM) which permits realistic calculations as well in the 3-D field as in the analysis of thin structural members. The co-rotational finite element method is particularly suited to problems involving large displacements and rotations but in which the strains remain moderate. The basic equations for a FEM solution of nonlinear dynamic problems by means of co-rotational coordinates are well known and were first derived by Belytschko and Hsieh [9]. The co-rotational technique consists in expressing element stresses and strains in non-deformable local systems of axes that rotate and translate with the elements. An orthogonal transformation matrix I T ] e is defined for each element e and relates local vectorial quantities to a global system of axes. In an element e, deformation displacements d" in convected axes are then related to the total displacement d in global axes by (~e = [ T ] e ((~ _ ( ~ ) ,

(6.1)

where (~ represents an approximate rigid body rotation and translation of the element. If a small strain hypothesis is used, it is possible to show [9] that the strain tensor ~iejmeasured in convected coordinates is related to the deforming displacements by the classical linear strain-displacement operator [B]: i, =- [B]

(6.2)

J. Donea et al. / Response of fast-reactor core subassemblies

propagation type problems prompts the choice of an explicit time integration. Indeed, though an implicit integrator would guarantee stability for time increments of any size, precision would only be achieved if the increments are sufficiently small or if an iteration procedure is used, resulting in high computer costs. In the case of an iterative procedure also the risk of nonconvergence exists. Maximum efficiency of an explicit scheme is achieved if it is possible to put the mass matrix [M] into a diagonal form. Eqs. (6.5) are then decoupled and no system solver is needed. Moreover, the combination of an explicit time integrator with a lumped mass matrix offers two advantages [10]: (i) the use of a diagonal mass matrix decreases the frequency response and compensates for the increases due to the explicit integrator, and (ii) the reduction in frequency response yields a less restrictive stability limit. Velocities and displacements are obtained by the following undamped integration scheme (central difference):

The stress tensor a~ is now derived using an adequate constitutive law and a set of internal forces ~e conjugated to the deforming displacements ~ can be obtained by means of the principle of virtual work:

ie=fv[B] o,j dV, T.e

(6.3)

where V~ is the element volume. The internal forces are transformed to global axes: (6.4)

I ~ = [ T ] ~ j e,

to permit assembly into the global internal force vector F which can be written as F = [K]~ where [K] is the stiffness matrix. Inertial forces being included in a D'Alembert sense, the system of equations to solve reads: [M] g=F-

F~xt,

(6.5)

where [M] is a mass matrix, ~ acceleration components (the dot - stands for total time derivation) and F~xt external forces. The evolution equations (6.5) are best integrated in time by a finite difference method. The presence of strong nonlinearities as are often encountered in wave

+~(~"

"lf

j

0 Fig. I. Euler-Bernouilli rectilinear beam element.

+ ~"+ ~),

"'+' =~" +At(~" +~n),

Y

r~

159

g

X

(6.6) (6.7)

160

.l. D o m , a et al. / Resiponse of..[ast-react°r core suhassemhhe.~

where n stands for time t and n + 1 for t + At, At being a uniform time increment. The Euler-Bernoulli beam theory is used in a corotational coordinates formulation (fig. I). The transverse displacement is assumed to be cubic in .~, the longitudinal displacement linear in 2. 0 is assumed to approximate the rigid body rotation of the element (small strain hypothesis). If ¢p~ and 4~2 are the total rotations obtained from the equilibrium equations, the deforming rotations (rotations relative to the .~ axis) are

~, =¢k, --0,

~2 =';b2 --0.

(6.8)

Writing the membrane strain as

I(l+ L)

(6.9)

where u = u~2 - u~ I,

v = Uy2 - u~.I

( u = displacements),

£ = x2 -- xl,

.V=Y2 --.vl

( x , y = coordinates),

and using normalized coordinates ranging from - 1 to 1, ~ = (2s - I)/l, ~ = 2f,/h, one obtains ~v ~ - ~ m - - h T / [ ( 3 ~

-

At<-lmin(l,~),

I)d~l -}- ( 3 ~ - t - 1 ) @ 2 ] .

(6.10)

(6.19)

where c is the sound wave speed in an unconstrained bar. For elements with an ideal aspect ratio l = 3 ~ , a practical condition can be derived:

A t = a - -/rain ,

2 ( u X + v y ) + U 2 nk V 2

~"~ =

In order to obtain stable solutions it is necessary to provide an estimate for the stability limit. If a central difference time integrator is used together with an elastic beam element assuming linear longitudinal and cubic lateral displacements, the following limitation is found [10]:

C

a ~< 1,

(6.20)

where lm~n is the shortest element length. Practical experience seems to indicate that elastoplastic problems do not exhibit poorer stability properties than purely elastic problems. As a general empirical rule it is advised to use a time step that is shorter than one half the critical value (a ~<0.45). It should also be noted that excessive refinement of the time step at constant mesh discretization degrades the quality of the response. As a matter of fact, the accuracy is maximized by operating as close to the stability limit as possible [101.

The internal nodal forces are obtained according to a virtual work principle [9]:

/,~= £,

-

7. Fluid-structure coupling

h xf.

dl~ drl,

h = ~ f a ~ d~ d~,

(6.12)

h 2 1"

nl I = -- -~- jvf (3t~ -- l) a~ dl~ d~,

h2 /'n 2 --

fie

-

-

8

fn(3~+ 1) 0~d~ dn,

3h 2 4l fne.

3h 2 r f2r = ~ - J ~ ° ~

(6.11)

de d,, d~ d~.

(6.13) (6.14) (6.15) (6.16)

All integrations are performed using a Gaussian rule with two points along the length and five through the thickness of the element. Lumped masses are obtained through geometrical considerations:

M t = M 2 = ohl,/2,

(6.17)

I, = 12 = 0h13/24.

(6.18)

Since finite element algorithms have been employed to analyze both the hydrodynamic and the structural parts of the coupled fluid-structure problem, the equations of motion in fluid and solid regions are all expressed in terms of nodal loads and are exactly of the same form. The equations of motion governing the coupled system may thus be solved simultaneously and a strong coupling between fluid and solid responses is achieved. However, suitable conditions must be prescribed along fluid-solid interfaces to allow relative sliding of the fluid and solid and it is the purpose of this section to show that fluid-structure coupling may be achieved in a very simple and elegant m a n n e r if the fluid is treated in the ALE formulation.

7.1. Coupling procedure In order to illustrate the coupling procedure along fluid-solid interfaces, let us consider a structural member embedded in a fluid which is allowed to slide along

J. Donea et al. / Response of fast-reactor core subassemblies

161

similar conditions at Lagrangian fluid nodes, as' "boundary conditions" for the automatic fluid mesh displacement algorithm described in section 5, Conditions (b) and (c) are enforced through the nodal loads. The condition of common normal velocity implies that l) 3n = /3In •

0 FLUID NODES

Y

S = STRUCTLIRAL ELEMENT

SO that v 3n = V nI = v ~ - -4v. - n

both faces of the structure. As shown in fig. 2, we place two nodes at each point of the interface: one fluid node and one structural node. Since the fluid is treated in the A L E formulation, the movement of the fluid mesh may be chosen completely independent of the movement of the fluid itself. In particular, we may constrain the fluid nodes to remain contiguous to the structural nodes, so that all nodes on the sliding interface remain permanently aligned. It is clear that such a permanent alignment of nodes along the interface greatly facilitates the flow of information between fluid and structural domains and permits fluid-structure coupling to be effected in the simplest and most elegant manner. The following conditions are prescribed at each point of the interface between an inviscid fluid and a deforming structure: (a) The grid velocity w of the fluid coincides with the velocity v s of the solid. (b) The normal velocity v~ of the fluid coincides with the normal velocity vg of the solid. (c) The tangential velocity vtF of the fluid is unconstrained. In order to specify the above conditions, a local coordinate system (t, n) is set up for each point of the interface (fig. 2), so that t is the tangent to the sliding interface and n is 90 ° clockwise from t. Whenever a corner occurs in the interface, t is the average of the two tangent directions, as explained in section 7.2. Referring to the node numbering in fig. 2, condition (a) implies that w4 = v 2,

(7 3)

/-9 I

x-

and

(7.2)

and since the solid element is based on the Kirchhoff hypothesis, we further have the equality , _-- / 2n2 ,

F = FLUID ELEMENT

Fig. 2. Sliding interfaces in ALE formulation.

w3 ----vI

v , ~-- 1 9 2n ,

and

STRUCTURAL NODES

(7.1)

where w indicates the grid velocity of the fluid and v the velocity of the solid. These conditions are imposed directly at each time step and serve, together with

(7.4)

Hence,the four nodes must have a common normal velocity which is obtained as follows. Let VF be the velocity of a fluid node at the interface and v s the velocity of the adjacent solid node. Since the two nodes must have a c o m m o n normal velocity at all times, we prescibe that d

[(v v -Vs)'n]

= 0

(7.5)

or

dn ( r F -- V s ) ' n + (VF -- VS) "-d-i = 0.

(7.6)

Now dn dt

n t+At --n

t

At

(7.7)

and since (v F -Vs),n t =0,

(7.8)

eq. (7.6) reduces to . .v .s v. ~. --

At

v? (7.9)

Thus, the condition of common normal velocity implies different normal accelerations at the fluid and solid nodes if the normal to the interface changes direction with time. To evaluate these normal accelerations we formulate -n = a n

VF

n

q-a F ,

vs'" = a" + a~,

(7.10)

where a n is a c o m m o n normal acceleration and a~, a S are individual complementary accelerations. The common normal acceleration a n is obtained

J. Donea el aL / Rev~onse of Jast-reactor core s'uhassemhfies

162 from the equation of motion

7.2. Interfaces with sharp corners

( M s + MF) a" = f " ,

(7.11)

where, referring to fig. 2, M s is the mass contributed by the structural elements adjacent to nodes 1 and 2, while M v = M w + MF2 represents the mass contributed by the fluid elements on either side of the interface The normal load f " is obtained by transformation into the (t, n) system of the assembled internal nodal loads ( f ~, f " ) contributed by the relevant fluid and solid elements: f n = ( f s ~ +ff~, + f ¢ 2 ) c o s a + ( f s ~' +f{:, +f'{2) s i n a ,

(7.12)

where a is the angle between x and n directions, as shown in fig. 2. In view of relations (7.10) and (7.9), the complementary accelerations a~ and a~ are linked by n aF

n __

v~ - v~ _ A~

--as

(7.13)

At

and a second relationship between these accelerations may be obtained from the requirement of equilibrium which reads

MFa ~ + Msa ~ = 0 .

(7.14)

It follows that " -

aF

m__Ms M F + M s An;

a~ -

-Mr

M v + Ms

A'.

(7.15)

This completes the computation of normal accelerations which ensure a c o m m o n normal velocity at the fluid and solid nodes on the interface. The tangential velocities at the fluid and solid nodes are unconstrained and result from the equations of motion

Mrlt~ =f~l, MF26~ = # 2 ,

Before concluding this section on fluid-structure coupling in the A L E description, we want to draw the attention of the reader to particular problems that arise when the interface between an inviscid fluid and a structure possesses sharp corners [16]. In ideal fluid flows, the velocity vector can have only one value and it is therefore impossible for a streamline to change direction abruptly. However, sharp corners in inviscid fluid flows are points where these rules are broken. The situation is sketched in fig. 3 which shows a corner in a wall, the wall turning away from the flow. At the corner, the streamline which coincides with the wall changes direction abruptly and the velocity vector has two directions. At this singularity, the velocity of the ideal fluid is infinite. In these conditions, it is easily understood that special precautions must be taken to achieve accuracy in the numerical prediction of the flow near a sharp corner. As mentioned before, whenever a corner occurs in an interface, the velocity vector at the corner is assumed to be parallel to an average of the two tangent directions. As can be seen from fig. 4, the consequence of this choice is that a spurious transport of fluid takes place across the intersecting parts of the walls, the fluid flows out on portion L I of the interface and flows in on portion L2. These spurious fluxes must indeed be accounted for when expressing the conservation statements of mass, m o m e n t u m and energy over the fluid finite elements meeting at the corner and this may prevent from achieving conservation in the fluid mesh if the amount of fluid lost on portion L I of the interface is not exactly compensated by the amount gained on portion L 2. Fortunately, exact compensation may always be achieved by a suitable choice of the angle fl defining the flow direction at the corner. Referring to fig. 4 and recalling that we are using linear fluid elements, the following equality must hold for exact compensation of the spurious fluxes: L i v I = L2v2,

(7.17)

(7.16)

where the tangential loads are again obtained by transformation into the (t, n) system of the assembled internal nodal forces (f.~,f.v). It may thus be concluded that the problem of fluidstructure coupling becomes very simple, since fluid nodes at the structure can remain attached to it thoughout the calculation without imposing any restriction on the fluid tangential velocity.

Fig. 3. Illustration of singularity at a corner in the flow of an inviscid fluid.

163

J. Donea et aL / Response of fast-reactor core subassemblies

FLUID

8.1.2. Plastic flow During plastic flow, the total strain increments are supposed to be linearly separable into elastic and plastic (incompressible) parts:

.

(8.3)

d , , j : d C + d,~,,,

so that the classical Hooke's law gives doij = Dij .... dee,,....

Fig. 4. Treatment of flow around sharp corners.

(8.4)

where Diy .... is the elastic matrix depending only on Young's modulus (E) and Poisson's ratio (v): Dijmn : 0

and since

except:

v t = v sin(a t - / 3 ) ;

v2=vsin(B-a2)

,

Do.,,,( i ----m =/:j = n ) : 2#,

%m,,(i:j=m:n) : X + 2 J , ,

(7.18)

(8.5)

D, jm,,( i :j=#= m : n ) : •,

the optimal value of angle/3 is found to be given by

with the Lam~ constants k and # defined by

L t sin a t + L 2 sin a 2 tan/3 = Li cos a t + L 2 cos a z '

(7.19)

Eu )~ = (1 + v)(1 -- 2 v ) '

E /~= G - 2(1 + v ) '

(8.6)

which reduces to /3: (at +a2)/2

(7.20)

in the particular case where L~ = L 2. Notice also that angle/3 in (7.19) is given by the slope of the line joining the end points of segments L~ and L 2 of the interface.

G or/~ being the shear modulus. Relations (8.3) imply that the strains remain small. The plastic rule relates the incompressible plastic strain increments ( d ~ ) to the current stress state (normality rule): dcPj = d y

8. Material modeling

30

O(~ij,

(8.7)

d ' / b e i n g a parameter of proportionality depending on strain increments.

8.1. Classical elasto-plastic material modeling

8.1.1. Yield criterion Most materials exhibit a reversible behaviour (often linear) up to a certain loading level above which metallic materials, for instance, do experience permanent strains (plasticity). The limit between elastic and elasto-plastic behaviour is represented in the stress space (a,j) by a yield surface:

f(ai)):O--ffy

:0,

(8.1)

where Oy is the elastic or yield limit in a tension test and 6 the equivalent stress. Using the Von Mises criterion the equivalent stress reads

O = [(olt + 022 + 033) 2 - 3(o t,o22 + 022033 -I- o t to33 - 0122 - 0223

8.1.3. Hardening rule The hardening rule explicitly defines the loading function describing the evolution of the yield surface during plastic loading:

F(o,,, 'r,) = o.

(8.8)

Two hardening rules are considered hereafter: isotropic and kinematic strain hardening. Comparative studies of various hardening rules [12,13] tend to indicate that agreement with experiments is generally less satisfactory for kinematic hardening than for isotropic hardening and do identify the so-called sublayer model as being the most adequate tool to take the Bauschinger effect into account. This latter model is also described. 8.2. "'Static" material behaviour

2 " "1 I / 2 --

o,3)1

.

(8.2)

8.2.1. lsotropic hardening First a hardening rule conserving isotropy is consid-

J. Donea et a/. / Response q/last-reactor core .sut)assembhe~

164

ered: the yield surface grows uniformly and keeps its original shape,

translation of the yield surface (initially equal to zero) given by (Prager's kinematic hardening [14])

F ( o , j , (P,j) = 6 - % ( , P ) = 0.

dc~ =-CdcP.

(8.9)

A given elasto-plastic stress state then always corresponds to the same yield surface independently of the deformation history. In a uniaxial stress situation this means that when reversed loading is applied (compression after a tension test), plastic deformation begins only when the compressive stress becomes equal to the last tension stress that produced plastic deformation. From an experimental point of view the isotropic model is realistic only as long as no reversed loading is applied: in this case the Bauschinger effect (apparent reduction of the yield point if loading is reversed after some plastic deformation has occurred) has to be taken into account. This phenomenon destroys the initial isotropy of the material. For isotropic hardening, the following incremental stress-strain relationship can be obtained [11]:

do,,=

D,j ....

9G 2 o,?(Ep+36) (%j-}3''°~a)

1 X ( % , , - ~3,,,,o,,)1 de .....

(8.10)

(8.t3)

However, difficulties are encountered for the quantitative definition of the hardening coefficient C. An incremental stress-strain relationship can be established for a constant linear plastic strain hardening Ep (eq. (8.11 )): 9G 2

D,i .....

do,,=

o,2( Ep + 3G)

(S,i--%,)

x (&,,, - ,~,,,,,)[ d,,,,,,.

(8.14)

!

8.2.3. Sublayer model In this model [15] the material is supposed to consist of superimposed "layers" all subjected to the same state of strain but with different elasto-perfectly plastic stress-strain curves (no hardening within a sublayer). For a material with a multilinear tension curve (,l segments), the sublayers (i) are defined as n elasticperfectly plastic materials and the stress is obtained as a weighted sum of the sublayer stresses: 11

o = ~ tio i

(8.15)

i=1

where do~) and de,) are stress and total strain increments, respectively, Du,,, the elastic matrix (8.5), oy the current yield limit (8.9) corresponding to the uniformly expanding yield surface, ojj the current total stress state, G the shear modulus (8.6) and E v the plastic modulus: Ep = E E t / ( E - - E , ) ,

(8.11)

E and E t being respectively the elastic and plastic tangent moduli. Numerical implementation of eq. (8.10) for finite increments includes a trial elastic step; if at the end of this step the yield condition is violated a corrective term is applied for the part of strain increments contributing to plastic deformation (determined by linear scaling).

E i - E,+ i t i --

In order to take into account the Bauschinger effect, a natural idea consists in replacing the isotropically expanding yield surface, discussed before, by a rigidly translating yield surface in deviatoric stress space (Si/):

) -% =0,

(8.12)

where % is the initial yield stress and a~ represents the

El

,

l<-i<~nwithE,,~l=O,

(8.16)

where Ej is the elastic Young's modulus of the original material. These relations derive from the requirement that continuous loading should reproduce the original stress-strain curve. Note that if the n th plastic slope E, is equal to zero, the n th layer vanishes as can be expected. All sublayers (i) behave elastically with a Young's modulus E, up to a yield stress given by i--I Oy i --

I _ OYi

8.2.2. Kinematic hardening

F=/az(Si,-%i)(S,.,--aii

with

2 tjolj J=' ~

Z

;

l
'

~ OyI

----oy.

(8.17)

tj

jet

When reversed loading is applied, it is obvious from the structure of the first layer that the Bauschinger effect is taken into account. An elastic-perfectly plastic material, i.e. one layer, can be modelled according to relation (8.10).

165

J. Donea et al. / Response of fast-reactor core subassemblies 8.3. Dynamic material behaviour If attention is restricted to small strain, the strain rate dependent mechanical behaviour of stainless steel can be idealized in the following way: - the elastic behaviour is strain rate independent; - the yield criterion depends on strain rate; - the hardening rule at constant strain can be strain rate dependent. It is further supposed that the normality law (8.7) still applies. Multiaxial strain rate can tentatively be defined by a simple scalar quantity:

d~ ~---- d~-'

(8,18)

It is then easily seen that the average strain rate during one step can be evaluated with reference to the purely elastic calculation (trial step) as

9.2. Hexcan material behaviour

The duct responses are largely dependent on the steel mechanical properties, which can vary in a wide range according to the operating conditions, the metallurgical history of the metal and the manufacturing of the duct. Only selected cases will be considered. 9.2.1. Static properties The unirradiated CW type 316 stainless steel of PHENIX hexcans has been well characterized in uniaxial tensile tests performed at "static" strain rate (3 × 10 4 s - i ) . Five mm gage length specimens were directly cut from standard ducts by an electrospark device and were set on a hydraulic machine applying a

350 Q_

--

-~----

EAt

/1

I

'

(8.19)

3OO

-",,,,. where #(o~/) is the equivalent stress at the end of the purely elastic step (trial), o(%/ - "-] ) the previous equivalent stress and At the time step. Relation (8.19) is general in the sense that it can obviously be used in both elastic and elasto-plastic steps.

[

(a)

J: 250

>

I

I

20O

[

I

l

l

0 10 Distance (ram)

I

~

30

12C,

9. Single hexcan response to internal pressure loading 9.1. Introduction

The structural part of the EURDYN-IM computer code [8] was checked on a series of out-of-pile experiments performed at the ZESFIR experimental area of C E N / C A D A R A C H E on small-length subassembly ducts (PHENIX hexcans) subjected to quasi-static and dynamic internal pressures at room temperature. The ~ 30 cm long duct sections were axially constrained in order to obtain plane strain conditions in the mid-sections. Tests were performed on ducts of two different 316 SS types - annealed and 20% cold-worked - with well-known static mechanical properties. The experimental equipment includes measurements of pressures, strains and displacements (only in the quasi-static cases). Hexcan responses under dynamic loading obtained with an explosive source were also compared to calculations performed with a strain rate dependent plasticity model based on the dynamic properties of AISI-316 stainless steel.

A

g_

~ ~

_..---- 2 3

Cold worked

Annealed

(b) i

I

J

I

4O Strain (%)

t

I 6O

Fig. 5. Material properties of 20% CW and annealed type 316 stainless steel, at room temperature: (a) Hardness test results, performed in the midwidth, of the ducts under a 100 N load, as a function of the position of measurements. (b) Engineering static stress-strain tensile test curves of CW and annealed steel, represented until the ultimate tensile stress. The CW steel is characterized by the 3 curves corresponding to the regions indicated in the insert.

J. Donea eta/. / Response of fast-reactor core subassemblies

166

tensile load at constant cross head speed until failure. As it was revealed by hardness diamond tests (see fig, 5a) that the drawing process during the manufacture of the cold-worked sections introduces a circumferential gradient of mechanical properties, it was decided to characterize the duct with 3 tensile test specimens belonging to the midflat, to the corner and to an intermediate zone of the flat. Annealed steel duct sections were prepared subjecting the CW hexcans to heat treatment in order to increase the ductility. As it was checked that the material in this case is quite uniform, only one tensile test was sufficient to characterize the annealed duct. Typical stress-strain curves describing the coldworked and annealed ducts at ambient temperature are presented in fig. 5b. The difference between the ultimate tensile strengths of the midflat and the corner reaches 15% and one can note the low strength of the annealed steel (see also ref. [17]). The 316 SS at 20°C is further defined by density Young's modulus Poisson's ratio

7890 k g / m 3, 194 GPa, 0.216.

9.2.2. Dynarnic properties Under dynamic loading, high strain rates can be reached depending on the slope of the pulse, on the steel and on the analyzed zone. Current strain rates are of the order of one thousand s i in the inner corners of the duct (region where stresses and strains are maximum) and several hundreds in the outer midflat (region of maximum displacements) as illustrated in fig. 6. The static uniaxial tensile test programme presented above has been extended to the dynamic regime up to strain rates of 1000 s i (work in progress). Unfortunately, these experimental data are still unavailable; standard properties of 316 SS obtained at J R C ISPRA [18] in the range up to 500 s ] were used. The experimentally observed strain rate hardening is taken into account by the following formula:

o(d, ¢ ) / 0 ( 0 , c) = 1 + ( ~ / c ) " ,

(9.1)

where o ( t , ( ) and o(0, c) are the dynamic and static stresses, respectively; c and n are two adjustable parameters, generally dependent on the strain c. Fig. 7 gives a representation of the dynamic to static stress ratio as a function of strain rate for selected strains. The c and n parameters, determined by least square fits, are found to be well represented by the following laws: c = exp[(6.4 + 8 6 . 2 c ) / ( 1 + 6.5~)],

800r

n = 0.23

--Q~

(9.2)

These data will be used as input in the strain rate dependent model (see section 9.4.2) in order to analyze the dynamic experiments.

' Cornel" inner s u r f a c e

60¢

(no strain dependence observed).

Midflat outer surface

9.3. Hexcan hydrostatic pressurization ~-- 40(;

Two static tests have been performed on annealed and CW materials. The main goal is to check the structural part of the computer code by comparison with simple experiments where reliable values of pressures, strains and displacements can be obtained without strain rate effects being taken into account.

E u~20 O

9.3.1. Annealed steel duct

I

0,5

T i m e (ms)

1

1,5

Fig. 6. Typical strain rates obtained in a 20% CW steel duct section submitted to a dynamic internal pressurization of 22 MPa (see section 9.4.3). The results are obtained with the EURDYN code.

A 35 cm long hexcan section is welded between two heavy and rigid flanges bolted together with large crosssectional area bolts reducing the end motions (fig. 8) to a minimm. As was pointed out in a previous publication [ 19], if the section is sufficiently long, such an apparatus permits to realize plane strain conditions in the vicinity of the duct mid-section (this was confirmed by axial strain measurements which were found to be an order of magnitude below the circumferential strains).

J. Donea et al. / Response of fast-reactor core subassemblies

2 -

167

Strain (%)

o

==

•

5

[]

lo

•

20

o

.u_ E >.

c~

n 10-5

I 10.3

l

I

i

10"1

1 101

i

I 103

Strain-rate (s -l)

Fig. 7. Strain rate effect for selected strain values, relative to type 316 stainless steel, fitted with formula (9.1) to JRC Ispra tensile test data [I 8].

Fig, 8. Experimental arrangements used in static (left-hand picture) and dynamic (right-hand picture) single hexcan experiments.

168

J. Donea et al. / Response of fi~st-reactor core subassemblies

Measurements of internal pressures, circumferential outer midflat strains and midflat and corner displacements were done by piezoelectric transducers, strain gages and inductance transducers, respectively. These data were doubled (pressures and displacements) or tripled (strains) to provide sufficient redundancy and to reveal possible differences due to material and geometry asymmetry. The natural period of the hexcan sections subjected to an internal pulse is of 0.25 ms. Using a loading time much larger than this period, it is possible to render negligible the dynamic effects and to analyze the static deformations with an explicit code (while the use of an explicit computer code is not recommended for this type of problems, it serves here for validation purposes). 1/12 of the duct was modeled with 9 beam elements

Test

/ EUNDYN ~ / . / 3.73 ~ calculation

e = 3.43

....

e

A

E .~

E

! 5

o

, , . ¢ ~ . _ .

/

/

;/

I

0

I

I

10 Internal pressure (MPa)

20

Fig. 10. Comparison of experimental and calculated strains versus internal hydrostatic pressure for an annealed 316 steel hexcan,

Peak pressure (MPa) = •

/'

j

24

z~ 18

o

12

,

6

30 Test ....

e = 3.43

EURDYN

$

._¢ 20

I (a)

I

lO3

10 2

lO4

10 S

Pulse rate (MPa/s)

10

20

~

,

"~ --

g

Test

//,

3 . 3 103 } ! Pulse rate 104 ?

o

///,1/, ///R'I

E (a) Midflat

E E C

o:

to

//

i5

10

//

o ~ i 5

t(b) l 10

15

20

25

Internal pressure (MPa)

Fig. 9. Pulse rate influence on annealed 316 steel hexcan responses to hydrostatic pressurization: (a) Kinetic energy of the duct as a function of the pulse rate for different final pressures. (b) Comparison of duct responses for different pulse rates. The shaded area corresponds to strain-gage dispersion.

{b) Corner

-to

i 0

I

I

I

10 20 Internal pressure (MPa)

i 30

Fig. 11. Comparison of (a) midflat and (b) corner displacements with EURDYN calculations as a function of internal hydrostatic pressure (annealed steel hexcan).

J. Donea et al. / Response of fast-reactor core subassemblies and plane strain conditions were assumed. The hexcan dimensions used in the analyses correspond to mean measured values (3.6 and 124 mm for the width and the external flat to flat distance, respectively). Several calculations have been made with different loading slopes ranging from 3.3 X 103 M P a / s (20 MPa in 6 ms) to 6 X 104 MPa/s. In fig. 9a the different values of the kinetic energy are compared; it is obvious that dynamic effects are negligible for slopes less than 3.3 x l03 MPa/s. The influence of the loading rise time on the midflat displacement is also shown (fig. 9b) as a function of the internal pressure. As was pointed out in an earlier publication [19], stress reversals occurred in the inner midflat and outer corner surfaces and for pressures above ~ 15 MPa, the hexcan deforms into a cylindrical shape. This means that neither isotropic nor kinematic strain hardening models are able to give correct responses for a given pressure level: in this case, it is needed to use the sublayer strain hardening model available in the EURDYN- 1M code. Figs. 10 and I1 show the comparison between EURDYN computations and experimental values for the mid flat outer strain and midflat and comer displacements. The comparison is made until the failure of the last strain gage (-- 22 MPa) or until the duct failure (occurring in a comer) at ~ 26 MPa. As the thickness dispersion of the hexcan varied between 3.43 and 3.73 ram, two calculations corresponding to these values were performed. The agreement is quite reasonable; the small remaining discrepancies could be attributed to non-uniform geometry and variable material properties a n d / o r limitations of small strain theory (at 20 MPa, the inner comer strain is of the order of 30%) and plasticity model. 9.3.2. Cold-worked duct section The calculation of another test using a 15% CW .type 316 stainless steel duct was also made. In this test, a 1 m long section was free to move axially: plane strain conditions are only approximately satisfied. The nonstandard 15% CW steel properties were derived from 20% CW 3 material data of fig. 5 using correlations between hardness and stresses (corrections of ~ 7%). The only measurements concern the flat to flat increase as a function of the applied pressure, which is compared in fig. 12 with 2 EURDYN calculations (performed with 2.5 X 103 M P a / s rise time) using a 3 and 1 material (from midflat) description, respectively. It can be seen that at l0 MPa, the cold-working gradient influence is as large as 60%. The calculation with the implicit INCA code [20] is also shown to be in good agreement with the

169

5r I|

o ....

Test INCA EURDYN 3 materials

4~-!

.....

EURDYN Midflat Properties

g

/

/

/

;'1-

//

'rk O

~

I

i I I 5 10 Internal pressure (MPa)

I 15

Fig. 12. Comparison of computed and experimental midflat displacement of a 15% CW 316 steel duct section subjected to an internal hydrostatic pressurization.

test results. This test, performed up to 15 MPa (without failure), was interpreted with the isotropic strain hardening model of the 2codes because the limited ductility of the steel did not permit complete stress reversals in the locations mentioned above. 9.4. Hexcan dynamic pressurization The aim of the EURDYN code is to simulate accidental conditions occurring in the fast dynamic loading regime. Therefore, two duct pressurization tests were performed with different materials using an explosive source which permits to produce pulses with l0 to 30 MPa peaks reached in ~ ! ms rise time. These pulses are not strictly representative of subassembly accidents but are used as a means of validating the codes. 9.4.1. Gas generator technique and apparatus Axisymmetric dynamic pressure loadings on ~ 30 cm duct sections are obtained by a low density explosive detonating in a perforated steel chamber. The highpressure gaseous products are controlled by the crosssectional area of the chamber holes and the axial venting of the apparatus. This permits the adjustment of the rise time of the amplitude and of the decay rate of the pulse. The measurements consist of pressures, which are

170

J, Donea et aL / Response of Jast-reactor core subassemhlies

found to be roughly uniform in the hexcan, and of outer midflat strains in the mid-section. To insure plane strain conditions, the same arrangement as in the static case was used (see fig. 8). 9,4.2. Strain rate dependent model Until now, there is no reliable and simple model related to steel constitutive equations able to take into account strain rate effects. In particular, it is not known how to move from one tensile test curve to another when the strain rate is rapidly varying. The simple formulation used in the E U R D Y N code consists in defining an equivalent strain increment at each step time, for every integration point and to calculate an equivalent strain rate; this permits to obtain the equivalent stress using the previously defined law (see section 9.2.2). In addition, the following assumptions are made: (i) If the current strain rate is larger than the previous one, loading occurs instantaneously on the elastic line and eventually on the plastic curve if the calculated equivalent stress is larger than the one obtained by the above law. (ii) If the current strain rate is smaller than the previous one, 2 cases are possible: - the strain decreases and the unloading occurs, elastically; - the strain increases and it is admitted that the stress point remains on the frontier of the plasticity surface which is contracting.

This procedure was used only in relation with the isotropic strain hardening model. 9.4.3. 20% Cold-worked hexcan test The experimental response of the hexcan to a 22 M P a pulse is compared in fig. 13 with two E U R D Y N calculations without strain rate effects (i.e. using static steel properties) and with dynamic properties and the strain rate model described above. The final strain discrepancy is reduced from ~ 35% to -- 8% when strain rate effects are taken into account. 9.4.4. Annealed duct dynamic loading In this test the section is submitted to a 25 MPa pulse in 1.5 ms. One can show that above l ms, stress reversals being present, it is not possible to obtain correct results. Restricting ourselves to the range less than I ms, experimental and computationl results using static and dynamic steel properties are compared (fig. 14). One can note the slight difference between the test results and the E U R D Y N calculation performed without strain rate effects. This tends to indicate that the strain rate dependence of the annealed steel is smaller than that of 20% CW steel: this is confirmed by preliminary results obtained in the uniaxial dynamic tensile test programme mentioned above. The remaining differences observed between these

•

~ c

4

°

f

....

//

!

'

/ i/Y ~/1

2,'

0.5

Fig. 13. C o m p a r i s o n

/e~%

ms

~

Test Static properties

....

Strain-rate model

1.0 Time (ms)

of experimental

E

5

O

1.5

/

// .

l 0,2

and

EURDYN

strain

result versus time, with and without strain rate model, for a 20% CW 316 hexcan subjected to the internal pulse shown in the insert.

---c----Test . . . . Static properties

/ /"

i

I

0.4

n

.

.

.

1

,t,o,o.,amo0e, to

0.6 Time (ms)

,

I

0.8

J

I

1.0

Fig. 14. Comparison of experimental and calculated midflat strain, as a function of time, for an annealed 316 steel hexcan. The loading pulse is indicated in the insert.

J. Donea et al. / Response of fast-reactor core subassemblies

171

Table 1 Material data for clustered hexcan problems Material

Initial density Po (kg/m3 )

Water

103

Constitutive law

p=p~[l

O'I4(I-V)] +0"28E/V-V-

(MPa)

Pl = 105A(0.1483 + 2.086A - 1-398A2)

A

=

t

2.086(1- V ) - 1 + ~ [ 2 . 0 8 6 ( 1 - V ) - 1 1 2 + 0 . 8 2 9 3 ( 1 - V) 2 - -

[ 0.0,

V=Oo/P; Stainless steel

7.9 × 103

V~I.0 .

.

.

.

.

.

when V= 1.0 E=energy/unit reference volume (MPa-m3/m3);

Poisson's ratio: v= I/3 Trilinear stress-strain law (isotropic hardening): E¢ -~2.0X 10 s MPa; Ep~=6.808× 103 MPa; Oylield=4.5X 102 MPa;

E 0 = 0.0

-2 -_ 1.g05X 103 MPa Ep

2 -O'yield - - 5 . 3 X 10 2 M P a

/ A

Fig. 17. Deformed configuration of simplified hexcan model at time t=0.6 ms.

Fig. 15. Simplified model of clustered hexcans.

2 tests and the E U R D Y N results will probably be reduced when more representative dynamic material data will be available and with some improvement in the development of the strain rate dependent model.

10. Application to clustered hexcans

10.1 Simplified model

Fig. 16. Finite element mesh for simplified hexcan model.

The first model is a simplified model, shown in fig. 15, which represents only the accident (pressurized) subassembly, the channel of fluid, and the adjacent subassembly. The loading consists of a triangular pressure pulse with rise and decay times of 0.25 ms and a

J. Donea

172

et

al. / Response of filst-reactor core suhasse, zhlze,s

z,.O

a

b

L.o

~0

3.O

2.0

2.0

N 'o

c

tO

OJO 0.0

o.,

i

6'

i

oi~

i

o'~

o'.~

l

l

0.0

o.~

0.0

Time (msec)

i

~l

i

012

r'-~,

, - ,

0.3

,

0.~ Time

,

,

0.5 (msec)

,

0.6

d

C

4.0

3.0

3.0

2.0 c

10 m

0.0 0,0

I

I

0.1

(3.2

~

Iw

I --

0,3

I - - - - -I % . , , ~l

O.t, Time

0,5 (reset)

I

0.0

O.G

O0

i

l

0,1

i

0.2

i

i

0.3

i

i

i

0.4 Time

O.G (msec)

l

i

0.6

Fig. 18. Time histories of longitudinal strains at the innersurface of the adjacent subassembly mid-faces.

peak pressure of 300 bar. The pressure is applied uniformly on the inside wall of the accident subassembly. The wall material is stainless steel and it is modeled by an elastic-plastic law in a state of plane strain. Isotropic hardening is included, but no strain rate effects are considered. The fluid between the hexcans is water; 8.0 3.0~

%

6.C

•

~.0

I" LO 2.0

.

0.0

o.I

0.2

O.3

13.4

0.5 Tir~ (mse¢)

0.6

'0.0

Fig. 19. Pressure and impulse time histories in the channel element on the horizontal symmetry line.

it is modeled as inviscid, compressible with an equation of state in the form p = f ( p , i). A pressure cut-off was included so that no tensile pressures are permitted. The numerical values used for the material data are given in table 1 *. The fluid nodes in contact with the structure are constrained to displace with the adjacent solid nodes, so that nodes along the interfaces remain permanently aligned. The intemal fluid nodes are treated in the ALE formulation and are moved automatically by the program as explained in section 5. The results for this problem are displayed in figs. 17-19. Fig. 17 shows the deformed configuration of the fluid-structure system at time t -- 0.6 ms and it is noted that the A L E description permits to keep the fluid mesh perfectly regular. The time histories of the strains at the inner surface of the adjacent subassembly mid-faces are represented in fig. 18. The pressure history in the channel element on the horizontal symmetry line and the corresponding impulse are given in fig. 19. * The finite element model for this problem is shown in fig. 16.

173

J. Donea et al. / Response of fast-reactor core subassemblies

L , LI

Fig. 21. Deformed configuration of complete model of clustered hexcans at time t= 1.0 ms.

Fig. 20. Complete model of clustered hexcans.

a

4.0 3.0

b 3.0[

:[~'-'

i

9 2.C~ E ".~ 1.0

0.0

0.0

d,G O.0

O.Z

0.4

0,6

0.8 (msec)

Time

-1.0 0.0

1.0

I

0%

'

OIj*

I

016

I

0~8 Time (msec)

I

1

1,0

d

C

/

c

3.(

3.C

z.c

2.c

c: "~ I.C

1.C

m

0,(

-"%

I

I 0.2

I

I

I

0.6

I

I

0.8 Time (msec)

I

|

1.0

-"°o.o

I

o12

I

oi,

I

I

o16 ~me

I 0.8

I

I 1.0

(mse¢)

Fig. 22. Time histories of longitudinal strains at the innersurface of the adjacent subassembly mid-faces. 10.2. Complete model

The complete model, shown in fig. 20, represents the pressurized subassembly, the channel of fluid, and the adjacent subassembly filled with water. The loading and material data are the same as for the simplified model.

Fluid-structure coupling is again treated in the ALE formulation with fluid and structural nodes on the interfaces permanently aligned. The automatic rezone procedure is employed to displace the nodes in the bulk of the hydrodynamic subdomains. The results for this problem are displayed in figs.

174

~

J. Donea et al. / Responxe of./d~t reactor core suhassembhe.v

2.5

1.0

2.0

0,8 x

x

0.6

i

=

~o

3.~

0,5

3.2

0.0 0.0

0,2

OA

0.6 Time

0.8 {rnsec)

1.0

E

3,0

Fig. 23. Pressure and impulse time histories in the channel element on the horizontal symmetry line.

21-23. Fig. 21 shows the d e f o r m e d configuration of the f l u i d - s t r u c t u r e system at time t = 1.0 ms. The time histories of the strains at the inner surface of the adjacent subassembly mid-faces are represented in fig. 22. The pressure history in the c h a n n e l element o n the horizontal s y m m e t r y line a n d the corresponding impulse are s h o w n in fig. 23.

References

[1] C.W. Hirt, A,A. Amsden and J.L. Cook, An arbitrary Lagrangian-Eulerian computing method for all flow speeds, J. Comput. Physics 14 (1974) 227. [2] J. Donea, Arbitrary Lagrangian-Eulerian finite element methods, chapter in Transient Analysis, eds. T. Belytschko and T.J.R. Hughes, to be published. (North-Holland, Amsterdam). [3] T. Belytschko and J.M. Kennedy, Computer models for subassembly simulation, Nucl. Engrg. Des. 49 (1978) 1738. [4] J. Donea, Finite element analysis of transient dynamic fluid-structure interaction, in: Advanced Structural Dynamics, ed. J. Donea (Applied Science Publishers, 1980) pp. 255-290. [5] L.R. Stein, R.A. Gentry and C.W. Hirt, Computer Methods in Applied Mechanics and Engineering 11 (1977) pp. 57-74.

[6] W.E. Pracht, J. Comput. Phys. 17 (1975) pp. 132-159. {7] W.E. Pracht and J.U. Brackbill, Report LA-6342, Los Alamos Scientific Laboratory (August 1976). [8] J. Donea, P. Fasoli-Stella, S. Giuliani, J.P. Halleux and A.V. Jones, The computer code EURDYN-IM for transient dynamic fluid-structure interaction, Part. I: Finite element modeling, Report EUR 6751, Commission of the European Communities (1980). [9] T. Belytschko and B.J. Hsieh, Nonlinear transient finite element analysis with convected coordinates, Internat. J. Numer. Meths. Engrg. 7 (1973) pp. 225-271. [10] S.W. Key, Transient response by time integration: Review of implicit and explicit operators, in: Advanced Structural Dynamics, ed. J. Donea (Applied Science, 19801. [11] F. Frey, Le calcul elasto-plastique des structures par ta methode des el6ments finis, Rapport No. 33, LMMSC, Universit6 de Liege (1973). [12] H. Armen, Plasticity in general purpose software, Workshop on Inelastic Constitutive Equations for Metals, Rensselaer Polytechnic Institute (April 1975). [13] J.A. Stricklin, B. Hunsaker and D.K. Vaughan, A coinparison of current work hardening models used in the analysis of plastic deformation, Workshop on Inelastic Constitutive Equation for Metals, Rensselaer Polytechnic Institute (April 1975). [14] W. Prager, A new method of analysing stress and strains in work-hardening plastic solids, J. Appl. Mech. 23 (1956) 493. [15] P. Duwez, On the plasticity of crystals, Phys. Rev. 457 (1935) 494. [16] J. Donea and S. Giuliani, Fluid-structure interaction in problems with interfaces involving sharp corners, Paper BI/3, 6th SMIRT Conf., Paris, August t7-21, 1981. [17] P. Granet et al., Virgin and irradiated hexagonal subassembly duct crushing test and analysis, Paper E2/6, 6th SMIRT Conf., Paris. August 17-21, 1981. [ 18] C. Albertini et al., Effects of irradiation on the mechanical properties of austenitic stainless steels under dynamic loading, Paper presented at the 9th Int. Syrup. on Effects of Radiation on Structural Materials, Richland, July 1978. [19] T.J. Marciniak et al., LMFBR subassembly response to simulated load pressure loadings, Nucl. Engrg. Des. 38 (1976) 1. [21] A. Combescure, Programme INCA, NT EMT 79/36.