Phahath Thammavong

Dr Marie Debacq
email: marie.debacq@cnam.fr

Stéphane Vitu

Marc Dupoizat

Experimental apparatus for studying heat transfer in externally heated rotary kilns

Keywords: Rotary kiln, rotary cylinder, lifter, flight, gas/solid reactor, heat transfer, heat transfer coefficient, pilot plant, experimentation

Rotary kilns are widely used in several branches of the chemical industry. In order to control the temperature of the solid and the gas flowing through the kiln, it is important to understand the heat exchange phenomena that occur. This article describes the design and construction of a novel experimental device to study heat exchange in rotary kilns. The device, which comprises a rotary kiln equipped with an external electrical heating system, enables the study of the influence of various parameters such as the solid flow rate, the kiln inclination angle, the rotational speed, or the presence of lifters on heat exchange and in particular on the heat exchange coefficient between the solid and the wall.

Preliminary experimental results concerning the influence of the solid flow rate and the rotational speed on the solid/wall heat exchange coefficient are here presented.

Introduction and background to the problem

For over a century [1], rotary kilns have been widely used in the inorganic chemistry industry. They are a key element in the production processes of cement, lime and pigments (titanium dioxide calcination) and in extractive metallurgy for the reduction of ore [START_REF] Boateng | Rotary Kilns[END_REF]. Rotary kilns are well suited for the pyrolysis of solid waste [START_REF] Descoins | [END_REF][START_REF] Descoins | [END_REF], and are also used to convert uranium fluoride into uranium dioxide for the manufacture of nuclear fuel.

Natural uranium cannot be processed directly because its content in isotope 235, the only fissionable isotope, is too low. To enrich the mineral, the uranium is therefore first converted into uranium hexafluoride (UF6). Enrichment in isotope 235 is carried out in the gas phase, either by gaseous diffusion or by centrifugation. The depleted uranium is then converted into a chemically stable and insoluble oxide, U3O8, and stored. The enriched uranium is converted into solid UO2 used to produce the solid fuel for nuclear power plants [5]. Both processes are conducted in conversion kilns, which comprise a hydrolysis reactor and a rotary kiln fitted with lifters. Rotary conversion kilns are externally heated by resistors located close to the outside wall of the rotating cylinder.

Various studies have been carried out on conversion kilns in an attempt to simulate their hydrodynamic, chemical and thermal behavior [5]. In this type of kiln, hydrodynamic, thermal and chemical modelling processes are closely linked, as it is clear that the temperature profiles of the solid and the gas in the kiln directly influence the chemical reaction kinetics that occur. In order to correctly model the operation of the kiln it is essential to be able to calculate precisely the temperature of the solid and of the gas at any point in the rotary kiln. Understanding the heat exchanges taking place in rotary kilns is therefore crucial.

Heat transfer between the solid charge, the counter-current gas and the wall of the kiln can take place by convection, conduction, and radiation. To calculate the heat balance, it is necessary to determine the various exchange surfaces as well as the emissivities, thermal conductivities, and heat transfer coefficients. Several studies [5][6][7] have shown that for conversion kilns the dominant heat transfer in a cross-section of the kiln is transfer between the solid and the wall, followed by solidto-gas transfer, and then by gas-to-wall transfer. Calculation of the temperature profile of the powder flowing through a rotary kiln therefore strongly depends on the solid/wall heat transfer coefficient (hsw) that is used.

Various correlations are reported in the literature for estimating the heat transfer coefficient hsw. Table 1 gives the hsw values calculated from some of these correlations [START_REF] Wachters | Proc. of the 3rd European Symposium on Chemical Reaction Engineering[END_REF][START_REF] Dumont | [END_REF][10][11][12] for the dimensions and typical operating conditions of an industrial kiln used for uranium conversion (diameter a bit less than 1 m and L / D ratio greater than 10). Most of the correlations mentioned in Table 1 have been obtained thanks to data measured in industrial kilns and thus can be used for the kiln chosen as an example. Other correlations have also been published [13][14][15][16], mainly based on penetration theory. It can be clearly seen in Table 1 that the heat transfer coefficient values obtained by using these correlations are highly scattered. Moreover, as experimental measurements of heat transfer in rotary kilns are scarce, the correlations were determined or validated from a limited amount of experimental data [10,[15][16][17].

The aim of the present study is therefore to design and construct an adaptable experimental setup to determine primarily the heat transfer coefficient between the wall of a rotary kiln and the solid flowing through the kiln. A key point in our approach is that the operating mode (and the main technical characteristics) of the experimental setup must be similar to those of the rotary kilns used for the conversion of uranium. In particular, our kiln is externally heated and equipped with lifters so as to investigate the influence of the latter on heat transfer.

The article is organised as follows: after a brief summary of experimental rotary kilns reported in the literature, the design of our experimental apparatus is described, focusing particularly on the heating mode chosen and on the instrumentation. The preliminary results obtained with this pilot-scale kiln are also presented.

Previous experimental equipments

As mentioned above, experimental studies of heat transfer in rotary kilns are relatively scarce. The most relevant ones involved large-scale kilns equipped with internal heating systems such as burners, or alternatively used air that had been preheated before being introduced into the rotating cylinder containing the solid. Other studies concerned kilns externally heated by steam or resistors.

There are currently three main types of experimental setups in use.

Internal heating using hot gas

Tscheng and Watkinson [10] constructed an experimental rotary kiln with a length (L) of 2.44 m and an internal diameter (D) of 0.19 m, giving an L / D ratio of 12.8. In this type of setup, the kiln was continuously fed with solid by a screw feeder. The kiln was also fed with preheated air flowing counter-current to the solid. As the water vapor and CO2 content of air is very low, gas to wall and gas to solids radiation could be neglected. These experiments were mainly conducted with sand (particle diameter ca. 0.7 mm).

Thermocouples were inserted in the bed, in the gas, and close to the rotating wall. The thermocouples were fixed to the rotating cylinder and positioned such that temperature measurements were made at the same axial locations. It was assumed that the gas phase and solid phase were in plug flow and at uniform temperatures at each axial position.

With this setup, the authors quantified the heat transfer coefficients from gas to wall (hgw), from gas to solid (hgs), and from solid to wall (hsw) by studying the influence of several parameters such as the inclination angle of the kiln, its rotational speed, the mass flow rate of the gas, the solid feed rate, and particle size.

Internal heating with a burner

Most industrial rotary kilns are internally heated by a burner. This was the heating mode adopted by Brimacombe and Watkinson [17][18][19] for the construction of a pilot rotary kiln 5.5 m in length and with an inside diameter of 0.406 m (giving an L / D ratio of 13.6). The burner was fed by a mixture of natural gas and air, and here again the solid used was sand with a mean particle size of 0.6 mm.

The kiln was equipped with 52 thermocouples for temperature measurements of the gas, solid, and wall. The thermocouples were fixed to the kiln wall and thus rotate with it.

As in the above study the solids temperatures and inner wall temperatures were very close together, the heat flow rate between the wall and the solid was neglected in order to simplify the data processing. The authors' main focus was on estimating the heat flux between the gas and the wall, and between the gas and the solid.

Dhanjal et al. [20] studied radial heat exchange in a bed of varying particle-size distributions inside a short kiln (L/D = 1.5) internally heated by a burner, and developed a model based on the active layer.

External heating

Few experimental studies have been undertaken on externally heated rotary kilns. To investigate heat transfer in this type of kiln, Lehmberg et al. [16] used a horizontal copper cylinder 0.6 m in length and 0.25 m in diameter (L / D ratio = 2.4). After filling the cylinder with solid, it was closed, set into rotation and heated to 200°C until homogenisation. The heating was then stopped and the copper wall of the kiln was cooled by spraying water on it. A constant wall temperature was maintained by spraying and the solid particles were cooled. During the cooling process, heat transfer was studied using thermocouples to register the temperatures of the wall and the particle bulk. Experiments were conducted with various sand particle diameters in order to study the impact of particle size on heat transfer. Suzuki et al. [21] measured the solid/wall heat transfer coefficient during the combustion of waste in an externally-heated batch-type rotary-kiln pyrolyzer (L = 1.5 m, D = 1 m; L / D = 1.5).

Their pyrolytic kiln was externally heated using the combustion exhaust of natural gas.

Temperature measurements taken at different locations within the thick cylinder wall showed that the heat exchange coefficients for waste composed of paper, wood, polyethylene, polypropylene, polystyrene and grain, when burnt at temperatures between 450 and 650°C, was about 20 to 50 W m -2 K -1 . They observed a slight increase in the coefficient with the number of lifters, their height, and the rotational speed of the cylinder [22].

Wes et al. [15] used an industrial-scale reactor (L = 9 m, D = 0.6 m; L / D = 15) to study heat transfer in rotary kilns. The cylinder was charged with potato starch and the kiln heated by steam flowing around the rotating cylinder. The cylinder temperature was recorded by a thermocouple and was taken to be constant along the cylinder length. The gas and solid temperatures were measured by several thermocouples positioned along the cylinder.

To estimate the heat transfer coefficient between the wall and the bulk solid (hsw), no air was passed through the cylinder during the experiments. The heat transfer coefficient between the wall and the gas (hgw) was measured in an empty kiln.

As the above summary makes clear, there is a lack of experimental studies of heat transfer in rotating cylinders. Not only have studies often been performed in specific conditions of heating mode or particle movement, but the results obtained are seldom comparable or applicable to different operating conditions.

We therefore decided to design a new experimental rotary kiln that is both versatile and fully instrumented so as to study heat transfer in conditions comparable to those encountered in the industrial rotary kilns used for the conversion of uranium.

Experimental device

General description

As can be seen on Figures 1 and2, the experimental device comprises a rotating cylinder in Inconel, an external heating system, a solid feed system and a tank for the recovery of the solid after it has passed through the kiln. The rotating cylinder has a length of 1.95 m, an outer diameter of 0.114 m and an inner diameter of 0.101 m. The L / D ratio is therefore greater than 19. The inside wall of the cylinder is smooth.

The device is equipped with a hopper for storage of the solid. The rotating cylinder is fed with solid by means of a loading chute and a feed screw located beneath the hopper. The feed screw and the cylinder are rotated by two variable frequency electric motors. The solid exiting the cylinder is collected in a tank. A control box placed at the lower end of the cylinder monitors the inflow of gas into the cylinder during the heating phase. The gas flow rate can be varied from 0 to

7 Nm 3 h -1 .
The outer wall of the rotating cylinder is heated by a furnace equipped with heating resistors. The furnace comprises two identical and independent heating zones (see Fig. 1), each measuring 40 cm and placed side by side. The heating resistors are insulated, the whole forming the heating device, which is placed around the rotating cylinder. The maximum power that can be delivered by the heating elements is 10 kW, with each zone delivering a power of 5 kW.

Variable parameters

The apparatus was designed to be adaptable so as to enable the study of the influence of many different parameters on heat transfer. Table 2 lists the variable parameters of the setup.

The angle of kiln inclination can be manually adjusted between 0° and 5°; this corresponds to industrial kilns, in which the tilt angle does not exceed 5°.

A wide range of rotational speeds is possible by means of the variable frequency motor. This range can be further extended by replacing one of the driving pinions of the cylinder. The rotational speed can be varied from 0.5 to 12 rpm. Note that the rotational speed of industrial rotary kilns is always well below 10 rpm. Obviously, it is difficult to compare an industrial kiln and a laboratory kiln operating at the same rotational speed since there is a large difference in solid flow rate.

However, we chose the widest possible range of rotational speed in order to avoid to be restricted in the study of this parameter influence.

The rotational speed of the feed screw is also adjustable within the range of 0.5 to 12 rpm.

Given the properties of the sand used in the experiments, the solid flow rate can be varied between 1.5 and 7.5 kg h -1 .

Dams can be placed at the end of the rotating cylinder, making it possible to select the height of the solid bed at the cylinder exit between 0 and 33.5 mm, and to vary the filling ratio of the kiln irrespective of the solid flow rate and the cylinder rotational speed.

With the heating power installed, the temperature of the solid passing through the kiln can reach 1000°C.

Since rotary kilns for the conversion of uranium are fitted with lifters, one of the goals of our study was to quantify the possible influence of the presence of lifters on heat transfer. The setup was therefore designed to operate with and without lifters. Removable lifters, made of the same material as the rotating cylinder, can be fixed to the inner wall of the cylinder and rotate with it in order to mix the solid. Figure 3 shows the lifters fitted to the rotating cylinder in the presence of solid.

Temperature measurements

Particular attention was paid to the instrumentation of the kiln, since in order to calculate the heat balance, it is essential to know the exact temperatures of the wall, the gas and the solid at any point in the kiln. To achieve this, the device was equipped with 36 thermocouples spaced in 4 series along the length of the cylinder (positions 1, 2(1), 2(2) et 2(3) on figure 1). 4 series of 4 K-type thermocouples were placed on the outside of the cylinder, as close as possible to the rotating wall. To insert the thermocouples, small holes were drilled in the heating unit. The tip of the thermocouples was located in the space between the heating resistors and the outer wall of the cylinder (see Fig. 4a and Fig. 5). The thermocouples are 3 mm in diameter and their tips are sheathed in ceramic to protect them from the radiation emitted by the heating resistors.

In order to validate the measurements recorded by the external thermocouples, the temperature of the outside wall of the cylinder was measured with a pyrometer via a hole in the heating unit.

Lastly, to measure temperatures up to 1300°C, 4 series of 5 N-type thermocouples 1.5 mm in diameter were placed inside the rotating cylinder on a rod made of the same material as the cylinder (see Figure 4b). These 5 thermocouples measure the temperature of the gas (3 thermocouples per series) and of the solid (2 thermocouples per series) (see Fig. 5). The rod is removable and is introduced along the kiln axis; its external diameter is 21 mm. The rod and internal thermocouples remain fixed during rotation of the kiln. To limit the mass and surface of metal inside the kiln, the dimensions of the rod were as small as possible. The axial positions of the internal and the external thermocouples were identical. Figure 5 shows the arrangement of the two series of thermocouples for a given cross-section of the kiln.

As the kiln can operate with and without lifters, two rods were constructed: one for the experiments without lifters, the other for the experiments with lifters. The main difference between the two rods is the distance between the tip of the thermocouples and the inside wall of the cylinder.

In the experiments conducted with lifters, the thermocouples need to be placed at a certain distance from the wall so as to avoid any contact between the rotating lifters and the fixed thermocouples.

It is important to note that only two internal thermocouples are in contact with the solid in the kiln, and that the rods are hollow to enable the passage of the extension wires connecting the thermocouples to the data acquisition system. The operating parameters of the device and the temperatures were recorded at regular intervals and then processed.

The angle of repose formed by the solid in the rotary kiln depends, among other things, on the shape and particle size of the solid. In order to adapt the technique to different types of solid, the angle of the rod supporting the internal thermocouples can be adjusted so as to maintain a right angle between the thermocouples inserted into the bed and the bed surface. The heating unit can also be adjusted slightly around the cylinder to ensure that the internal and external thermocouples are always opposite each other.

Materials

The physical properties of the sand used for the experiments are listed in Table 3. The mean particle size was 0.54 mm, the bed density was 1420 kg m -3 and the static angle of repose approximately 39 degrees. The static angle of repose of the sand has been measured using the fixed cone method: the powder flows out of a funnel and forms a heap from which the static angle of repose can be estimated.

Thermal conductivity, specific heat and sand emissivity were not measured; the values given in table 3 are those reported by [START_REF] Descoins | [END_REF]10,12].

Preliminary tests and results

Axial and transverse motion of sand

Before undertaking the heat study experiments, it was important to check the correct hydrodynamic operation of the pilot kiln and to quantify the mean residence time. We therefore conducted repeated measurements of the residence time distribution under the nominal operating conditions of the setup, to ascertain that the flow in the pilot was a plug flow with low axial dispersion. The pilot was therefore operated cold for at least 5 hours with a solid feed rate of 2.5 kg h -1 , a rotational speed of 3 rpm, a 23.5 mm dam and an inclination angle of 2°. As it could then be considered that from a hydrodynamic point of view the steady-state regime had been reached, we promptly introduced 12 g of red sand grains into the inlet of the rotating cylinder (just downstream the feed screw outlet). We thus made a pulse tracer input signal. For residence time distribution experiments, it is necessary that the tracer properties be the same as those of the normal charge entering the system. That is why coloured sand was used as tracer in our experiments.

After the impulse tracer input, samples of 3.4 g of sand were collected at the cylinder outlet every 30 s. To describe the concentration-time function of tracer in the exit stream of the cylinder, the red grains were sorted out from the samples and weighed with a precision balance, enabling the residence time distribution function shown on Figure 6 to be determined. Indeed, at each collected sample corresponds a time (ti) and a mass of tracer (mi). The RTD function shown on Figure 6 was obtained using the equation

      dt t m t m t E function RTD 0 i i     (1)
The results obtained demonstrate first that the trial can be correctly reproduced, and second that peak broadening is negligible. The mean residence time in these operating conditions is 46.5 min and the variance is 2.2 min 2 (6 10 -4 s 2 ). The Peclet number is therefore extremely high (around 2000), as is usual for rotary kilns with a high L / D ratio [23][24][25]. The flow regime in our pilot-kiln can thus be considered to be a plug flow with low axial dispersion. The results further show that five hours is ample time to establish the regime before undertaking RTD measurements.

Concerning the transverse bed motion, visual observations allowed us to conclude that the bed is in rolling mode in the investigated rotational speed and kiln inclination range. Indeed, the bed shows a constant angle of inclination and a flat planar surface. There is no slipping of the bed against the wall. Moreover, the Froude number was calculated under the nominal operating conditions (above-mentioned) of our experimental device. The result obtained is Fr = 1.6 10 -4 . This value is in the range of Froude numbers corresponding to the rolling mode according to Mellmann [26]. Even if the Froude number is an insufficient criterion to characterize bed behaviour, this result confirms our visual observations.

Validation of the steady state

The study detailed above enabled us to quantify the time required to establish the hydrodynamic regime. With a near-perfect plug flow, it was thus ascertained that twice the mean residence time is more than sufficient to attain steady-state from the hydrodynamic point of view.

The next step was to determine the heat regime. The pilot was therefore set in operation with a solid feed rate of 2.5 kg h -1 , a rotational speed of 3 rpm -1 , a dam of 23.5 mm and an inclination angle of 2°. A further specification was to apply a temperature of 600°C to zone 2. Figure 7 shows the temperature readings of certain thermocouples as a function of time (cf. figure legends on Figures 1 and5). It can be seen on Figure 7a that after 5 hrs the temperatures appear to level off.

Figure 7b shows the profiles obtained after 5 hrs in greater detail, confirming that steady state is reached in five hours. The power box displays a perfectly constant value of power consumed (to within 10W, out of 0.83 kW), thus confirming that when the temperatures are stable, steady state has been attained from the thermal point of view.

In all our experiments to determine the value of the solid/wall heat transfer coefficient, the pilot-kiln was therefore operated for 5 to 5.5 h and the time-averaged temperatures were then calculated during the half-hour following this period in order to represent the steady state operation of the setup.

Heat balance

Measurements were made in steady state without lifters, with a continuous flow of sand and no circulation of air (the gas inlet and outlet of the rotating cylinder were closed).

The measurement points 2(1) and 2(3) are considered respectively to be the entry and exit points of our system where the heat balance was determined; the measurement point 2 [START_REF] Boateng | Rotary Kilns[END_REF] corresponds to the inside of the system. The heat flow required to heat the solid moving through the system is therefore:

  ) 1 ( 2 ) 3 ( 2 s s Ts Ts Cp m    • (2)     4 ) 2 ( 2 4 ) 2 ( 2 sg g w s Tw Tp A F 1             (3)
The net radiative flux between the solid and the gas is:

  4 ) 2 ( 2 4 ) 2 ( 2 sg g s Tg Ts A         (4)
The heat fluxes exchanged by modes other than radiation between the solid and the wall, and between the solid and the gas, are respectively:

  ) 2 ( 2 ) 2 ( 2 sw sw Tw Ts A h    (5) and 
 ) 2 ( 2 ) 2 ( 2 sg sg Tg Ts A h    (6)
The mass flow of the sand s m

• was imposed and checked at the outlet of the cylinder. The heat capacity of the sand Cps is given in Table 3. The solid/gas heat transfer coefficient hsg was estimated to be 3 W m -2 K -1 using a correlation that is applicable to natural convection above a horizontal solid plane for air in the temperature range measured by our pilot-kiln [START_REF] Perry | Perry's chemical engineers' handbook[END_REF]; the solid/wall heat transfer coefficient hsw is the magnitude to be determined. The heat exchange areas between solid and gas Asg, and between solid and wall Asw, were calculated as a function of the filling degree of the cylinder; the filling degree was experimentally determined by completely emptying the cylinder after stopping the sand feed. As the mass and the bulk density of the sand in the cylinder are known (see Table 3), the volume occupied by the solid can be calculated, and was found to be:

L 2 sin R 2      (7)
R is the cylinder radius,  the filling angle and L the total length of the rotating cylinder (the error due to the shape of the powder bed is minimal as the bed is almost flat). Using the value of the filling angle thus obtained, the exchange surfaces for the heat balance can be calculated:

balance sg L 2 sin R 2 A            (8) and balance sw
L R A     (9)
Lbalance is the length of the system on which the heat balance was measured, i.e. the distance between the measurement points 2(1) and 2(3), or 200 mm. The Stefan-Boltzmann constant was taken to be 5.67 10 -8 W m -2 K -1 . The emissivity of the cylinder (Allegheny Ludlum 800HT, nickel-base alloy) is 0.9 [manufacturer's data]. The emissivity of the sand is given in Table 3. The air is mainly composed of symmetrical diatomic gases and can therefore be considered to be transparent.

Given the ambient humidity in the laboratory, the emissivity of the air is below 0.01 [START_REF] Welty | Fundamentals of Momentum, Heat and Mass Transfer[END_REF], meaning that the radiative exchange term with the gas can be neglected. The correction factor for our geometry [START_REF] Vriendt | La transmission de la chaleur -Introduction au rayonnement thermique, 3 e Ed[END_REF] is:

wg sg w w s s A A ε ε 1 1 ε ε 1 1 F       with   balance wg L R Γ π 2 A     (10)
All the temperatures were obtained from the time-averaged temperatures measured by the thermocouples during the half-hour following at least 5 to 5.5 h of steady-state regime. The temperature of the solid, Ts, was calculated by averaging the temperature registered by the thermocouple inserted in the sand close to the wall (Tsw on Fig. 5) and the thermocouple inserted in the sand slightly above (Tsi) when the difference between the two was less than 10°C. Any difference in excess of 10°C means that the thermocouple measuring Tsi does not in fact penetrate the sand because the bed is too shallow; in this case we took Ts = Tsw. The temperature of the gas, Tg, was measured by the thermocouple plunged in the gas in close proximity to the solid (Tg on Fig. 5). The wall temperature, Tw, was calculated by averaging the temperatures registered by the thermocouples placed around the wall (they are close to one another, as can be seen in Table 4).

The heat balance in steady state can therefore be written:

              0 Tw Ts Tw Ts F 1 A Tw Ts A h Tg Ts A h Ts Ts Cp m 4 ) 2 ( 2 4 ) 2 ( 2 g s 4 ) 2 ( 2 4 ) 2 ( 2 g w s sg ) 2 ( 2 ) 2 ( 2 sw sw ) 2 ( 2 ) 2 ( 2 sg sg ) 1 ( 2 ) 3 ( 2 s s                                • (11)
Neglecting the radiative exchange term with the gas as explained above, the solid to wall heat transfer coefficient can be written:

          ) 2 ( 2 ) 2 ( 2 sw 4 ) 2 ( 2 4 ) 2 ( 2 g w s sg ) 2 ( 2 ) 2 ( 2 sg sg ) 1 ( 2 ) 3 ( 2 s s sw Ts Tw A Tw Ts F ε 1 ε ε σ A Tg Ts A h Ts Ts Cp m h                    • (12)
The values used for calculation in the nominal operating case are given in Table 4; the repeated trials demonstrate the accuracy of our results. The results obtained by varying independently the mass flow rate s m

• and the rotational speed of the cylinder N are given in Table 5. It can be seen that the heat transfer coefficient hsw varies only slightly as a function of N in the range investigated, whereas Figure 8 exhibits a linear correlation between hsw and s m

• for a constant rotational speed of 3 rpm.

These results may seem surprising since several previous studies reported a significant effect of the rotational speed on the wall to solid heat transfer coefficient. However, in our study, the investigated rotational speed range (from 1 to 3 rpm) is a bit narrow and we are planning to extend the investigated range soon.

FIGURESFigure 1 .Figure 2 .

 12 FIGURES

Figure 3 .

 3 Figure 3. Photograph of the experimental rotary kiln equipped with lifters.

Figure 4 .

 4 Figure 4. a) External temperature measurement. Thermocouples are placed between the electrical resistance and the external cylinder wall. b) Internal temperature measurement. At a given section, five thermocouples are used to measure the solid and gas temperatures at different radial positions in the rotary kiln.

Figure 5 .

 5 Figure 5. Thermocouple positions in a cross section of the rotary kiln. 4 external thermocouples to measure the wall temperature. 3 internal thermocouples to measure the gas temperature in the kiln. 2 internal thermocouples to measure the bed temperature. This measurement device is repeated four times along the kiln in the heating zone. The positions of the four devices are specified in Fig 1.(1, 2(1), 2(2) and 2(3)).

Figure 6 .Figure 7 .

 67 Figure 6. Residence Time Distribution repeated twice in the same conditions ( s m • = 2.5 kg h -1 , N = 3 rpm, dam = 23.5 mm and  = 2°).

Figure 8 .

 8 Figure 8. Wall to sand heat transfer coefficient as a function of sand mass flowrate at a constant rotational speed of 3 rpm.

  With regard to the influence of the solid flow rate, a comparison of our results with previous literature results is very tricky. Actually, the operating conditions used in this study are very far from those used in the few experimental previous works.Finally, let us stress that this work has been conducted because in the open literature, experimental measurements concerning the wall to solid heat transfer coefficient in our conditions (i.e. external heating, thin wall without layer of refractory) are almost nonexistent.

	SYMBOLS USED
	A	[m 2 ]	area
	Cp	[J kg -1 K -1 ]	specific heat
	D	[m]	kiln diameter
	Fr	[-]	Froude number
	hsg	[W m -2 K -1 ] gas to solid heat transfer coefficient
	hgw	[W m -2 K -1 ] gas to wall heat transfer coefficient
	hsw	[W m -2 K -1 ] wall to solid heat transfer coefficient
	L	[m]	kiln length
	Lbalance [m]	length of heat balance system
	m	[g]	mass of tracer in a sample
	• m	[kg s -1 ]	mass flowrate
	N	[rpm]	rotational speed
	R	[m]	inside radius of cylinder
	T	[K]	temperature
	t	[min]	time
	Greek letters
		[rad]	kiln inclination angle
		[-]	emissivity
		[m]	dam height
		[W m -2 K -1 ] Stefan-Boltzmann constant
		[rad]	filling angle
	Subscripts
	2(1)	first measurement position in the second heating zone
	2(2)	second measurement position in the second heating zone
	2(3)	third measurement position in the second heating zone
	g	gas
	s	solid (sand)
	w	wall

Table 1 .

 1 TABLES Wall to solid heat transfer coefficient obtained using different correlations.

	Correlation	hsw [W m -2 K -1 ] (order of magnitude)
	Wachters & Kramers [8]	44
	Dumont & Belanger [9]	30
	Tscheng & Watkinson [10]	185
	Imber & Paschkis [11]	0.1
	Lybaert [12]	0.08

Table 2 .

 2 Variable parameters of the experimental device.

	Parameter	Range
	Kiln inclination	0 to 5 degrees
	Kiln rotation speed	0.5 to 12 rpm
	Screw rotation speed	0.5 to 12 rpm
	Solid feed rate (with sand)	1.5 to 7.5 kg h -1
	Bed Height	0 to 33.5 mm
	Internal kiln temperature	Up to 1000°C
	Lifters	Removable

Table 3 .

 3 Physical properties of the solid used in this study.

	Material	Sand
	Mean particle diameter	0.54 mm
	Bulk density	1420 kg m -3
	Static angle of repose	39 degrees
	Thermal conductivity [10]	0.268 W m -1 K -1
	Specific heat [12]	835 J kg -1 K -1
	Emissivity [3]	0.76

Table 4 .

 4 Calculation of wall to sand heat transfer coefficient ( s m • = 2.5 kg h -1 , N = 3 rpm,  = 23.5 mm and  = 2°).

		rotational speed (N)	rpm	3	3	3
		sand mass flowrate (m s ) 	kg h -1	2.48	2.48	2.44
		power	kW	0.83	0.83	0.83
		filling degree	%	9	9	9
		Tsw	°C	480	481	481
	2(3) 2(1) 2(2)	Tsi Ts 2(1) (mean) Tsw -Tsi Tw1 Tw2 Tw3 Tw 2(2) (mean) max -min Tg' Tg = Tg 2(2) Tg' -Tg Tsw Tsi Ts 2(2) (mean) Tsw -Tsi Tsw Tsi Ts 2(3) (Tsw)	°C °C °C °C °C °C °C °C °C °C °C °C °C °C °C °C °C °C	478 479 2 582 582 585 583 3 543 519 23 550 549 549 1 566 540 566	479 480 2 583 583 586 584 3 545 523 21 550 549 550 1 566 542 566	480 481 1 583 583 584 583 1 537 514 23 550 550 550 0 566 544 566
		Tsw -Tsi	°C	25	24	22
			rad	1.57	1.57	1.57
		A sg	m²	0.0142	0.0142	0.0142
		A wg	m²	0.0476	0.0476	0.0476
		A sw	m²	0.0158	0.0158	0.0158
		m s .Cp s .(Ts 2(3) -Ts 2(1) ) 	W	49.89	49.52	48.07
		h sg .A sg .(Ts 2(2) -Tg 2(2) )	W	1.27	1.13	1.53
		radiation	W	-48.43	-48.27	-47.25
		A sw .(Ts 2(2) -Tw 2(2) )	m² K	-0.540	-0.538	-0.527
		wall to sand heat transfer coefficient (h sw ) W m -2 K -1	5.1	4.4	4.5

Table 5 .

 5 Wall to sand heat transfer coefficient as a function of sand mass flowrate ( s m  = 23.5 mm and  = 2°. .48 2.44 2.88 3.18 3.18 3.35 3.59 2.47 2.49 2.50

	•
	) and



Acknowledgements

This work was supported by a grant from Areva NC. We gratefully acknowledge Pr. Bontemps, GRETh, for his kind help in this work.

Conclusions

We have designed and built an adaptable experimental device to determine the heat transfer coefficient between the wall and the solid flowing through a rotary kiln. The kiln is externally heated by an electrical heating unit comprising two independent zones to allow both the heating and the cooling of the solid charge to be studied; the rotating cylinder can also be fitted with lifters. The following parameters can be independently varied: rotational speed of the cylinder, solid flow rate, kiln inclination, bed height at the cylinder exit, and the temperatures of the two heating zones.

Lastly, the device can be adapted to the specific characteristics of the powder used (in particular, the angle of repose).

The experimental procedure and the calculations developed in this study have enabled the determination of the solid/wall heat transfer coefficient under the trial conditions, thanks to temperature measurements and various other parameters. We have established that, in the nominal operating conditions of our system, and using sand as the solid charge, the solid/wall heat transfer coefficients obtained are of a comparable order of magnitude to those reported in the literature, under more or less similar operating conditions. Among the variable parameters investigated, we have found that the solid/wall heat transfer coefficient increases linearly with the solid flow rate, but that it shows negligible variation with the rotational speed of the cylinder. Future work will focus on the influence of the other variable parameters of the setup, namely the kiln inclination, the bed height at the cylinder exit, and the heating temperature. We also plan to examine the case of cooling, the effect of particle size, and the presence of lifters


Léo Morin
email: leo.morin@ensam.eu

Jean-Claude Michel

Jean-Baptiste Leblond

A Gurson-type layer model for ductile porous solids with isotropic and kinematic hardening

Keywords: Porous ductile solids, Isotropic hardening, Kinematic hardening, Sequential limitanalysis

 for isotropic hardening to mixed isotropic/kinematic hardening. The model is finally assessed through comparison of its predictions with the results of some micromechanical finite element simulations of the same cell. First, the numerical and theoretical overall yield loci are compared for given distributions of isotropic and kinematic pre-hardening. Then the predictions of the model are investigated in evolution problems in which both isotropic and kinematic hardening parameters vary in time. A very good agreement between model predictions and numerical results is found in both cases.

A C C E P T E D M A N U S C R I P T 1 Introduction

The failure of metals, whose impact on the integrity of engineering structures need not be stressed, is one of the most challenging problems faced by the scientific and industrial communities. Indeed its analysis and modelling are complex tasks because it is a multiscale problem. Various mechanisms, at the microscale (e.g. structures of dislocations or grain boundaries) and mesoscale (e.g. hard precipitates or voids), can induce damage leading ultimately to macroscopic cracks. In particular, a difficult but essential task consists in providing predictive micromechanically-based models that permit to account for both monotonic and cyclic loadings.

In the case of ductile materials considered in this paper, failure essentially takes place in three steps (see e.g. [START_REF] Benzerga | Ductile fracture by void growth to coalescence[END_REF]; [START_REF] Pineau | Failure of metals I: Brittle and ductile fracture[END_REF]; [START_REF] Benzerga | Ductile failure modeling[END_REF] for recent reviews of the topic): (i) the nucleation of voids, (ii) their growth, change of shape and rotation, and finally (iii) their coalescence leading to final failure.

The modelling of these mechanisms started with the pioneering work of [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I-Yield criteria and flow rules for porous ductile media[END_REF], who combined homogenization and limit-analysis of a hollow sphere made of a rigidideal-plastic isotropic material to derive a model of ductile materials incorporating void growth. This model has been extended in various directions to account for phenomena not included in its original version, notably void shape effects [START_REF] Gologanu | Approximate models for ductile metals containing non-spherical voids-Case of axisymmetric prolate ellipsoidal cavities[END_REF][START_REF] Gologanu | Approximate models for ductile metals containing nonspherical voids-Case of axisymmetric oblate ellipsoidal cavities[END_REF][START_REF] Gologanu | Recent extensions of Gurson's model for porous ductile metals[END_REF][START_REF] Garajeu | A micromechanical approach of damage in viscoplastic materials by evolution in size, shape and distribution of voids[END_REF]Madou andLeblond, 2012a,b, 2013;Madou et al., 2013), plastic anisotropy [START_REF] Benzerga | Plastic potentials for anisotropic porous solids[END_REF][START_REF] Monchiet | Macroscopic yield criteria for plastic anisotropic materials containing spheroidal voids[END_REF][START_REF] Keralavarma | A constitutive model for plastically anisotropic solids with non-spherical voids[END_REF]Morin et al., 2015c) and void coalescence (Thomason, 1985;Tekoglu et al., 2012;[START_REF] Benzerga | Effective yield criterion accounting for microvoid coalescence[END_REF]Morin et al., 2015b). It has met, in both its original and improved forms, considerable success in the reproduction of experimental tests of failure of ductile materials under monotonic loading conditions.

The failure of ductile metals under cyclic loadings is less well understood and mastered.

Experiments (Schmidt et al., 1991;[START_REF] Kobayashi | The cyclic J-R curve and upper-limit characteristic of fatigue-crack growth in 2-1/2 Cr-Mo steel[END_REF] have shown that the strain to fracture is considerably lower, for a given load, if it is reached under cyclic conditions rather than monotonically. This reduction of ductility is commonly attributed to an effect of gradual increase of the mean porosity (volume fraction of voids) during each cycle

A C C E P T E D M A N U S C R I P T

termed the ratcheting of the porosity. This phenomenon was first evidenced in micromechanical finite element simulations performed by [START_REF] Gilles | Analysis of cyclic effects on ductile tearing strength by a local approach of fracture, in: Advances in Fracture/Damage Models A C C E P T E D M A N U S C R I P T for the Analysis of Engineering Problems[END_REF] under conditions of constant overall triaxiality (in absolute value), and later confirmed by several authors [START_REF] Devaux | On continued void growth in ductile metals subjected to cyclic loadings[END_REF][START_REF] Besson | An extension of the Green and Gurson models to kinematic hardening[END_REF][START_REF] Brocks | Damage models for cyclic plasticity[END_REF][START_REF] Rabold | Cell model simulation of void growth in nodular cast iron under cyclic loading[END_REF]Steglich et al., 2005;[START_REF] Mbiakop | On void shape effects of periodic elasto-plastic materials subjected to cyclic loading[END_REF][START_REF] Lacroix | Numerical study and theoretical modelling of void growth in porous ductile materials subjected to cyclic loadings[END_REF]. As explained by [START_REF] Lacroix | Numerical study and theoretical modelling of void growth in porous ductile materials subjected to cyclic loadings[END_REF], the ratcheting of the porosity is fundamentally tied to two features of the material behaviour, namely strain hardening and elasticity. The effect of elasticity, although important in the context of the ratcheting of the porosity under cyclic loadings, will not be considered in this paper, its study and modelling being postponed to a future paper. We shall thus focus exclusively on the effect of strain hardening, as a first step toward a complete modelling of ductile rupture under cyclic loading conditions. [START_REF] Devaux | On continued void growth in ductile metals subjected to cyclic loadings[END_REF], among other things, showed that [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I-Yield criteria and flow rules for porous ductile media[END_REF]'s classical model does not predict the effect of ratcheting of the porosity under cyclic loadings, but a stabilization of the evolution of the porosity right from the first semi-cycle. The explanation of this shortcoming of [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I-Yield criteria and flow rules for porous ductile media[END_REF]'s model lies in the crude modelling of strain hardening within this model, and more specifically in the fact that the same "average yield stress of the matrix" appears in both the "square" and the "cosh" terms of the yield function.

A few theoretical works have tried to improve the modelling of strain hardening effects within [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I-Yield criteria and flow rules for porous ductile media[END_REF]'s model:

• For isotropic hardening, [START_REF] Leblond | An improved Gurson-type model for hardenable ductile metals[END_REF] proposed a heuristic extension of [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I-Yield criteria and flow rules for porous ductile media[END_REF]'s yield function involving distinct "average yield stresses of the matrix" in the "square" and "cosh" terms. Their approach was based on some approximate analytical solution to the problem of a hollow sphere made of some rigid-hardenable material and subjected to some arbitrary loading. Their model, which accounts for the heterogeneity of isotropic hardening within the matrix, has notably permitted to qualitatively reproduce the ratcheting of the porosity observed in micromechanical finite element simulations of porous cells subjected to cyclic loadings under conditions of constant overall triaxiality (in absolute value) [START_REF] Leblond | An improved Gurson-type model for hardenable ductile metals[END_REF]. From a quantitative point of view, however, the comparison was not fully satisfactory. This obviously arose from the fact that the model was not a priori designed for cyclic loadings, since it was based on an assumption of positively proportional straining which is inadequate for

A C C E P T E D M A N U S C R I P T
such loadings. This hypothesis was relaxed by [START_REF] Lacroix | Numerical study and theoretical modelling of void growth in porous ductile materials subjected to cyclic loadings[END_REF], at the expense of introduction of radial discretization of an underlying mesoscopic hollow sphere and calculation and storage of the hardening parameters in each of the spherical layers thus defined. This resulted in a much improved agreement of model predictions and results of micromechanical finite element simulations.

• For kinematic hardening, [START_REF] Mear | Influence of yield surface curvature on flow localization in dilatant plasticity[END_REF] introduced a macroscopic backstress in [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I-Yield criteria and flow rules for porous ductile media[END_REF]'s model. However, this backstress was not linked to some heterogeneously distributed microscopic counterpart in the matrix, and its evolution equation was complex: it was adjusted in such a way that the isotropic (Gurson) and kinematic models yielded identical predictions for proportional loadings. The aim of this kinematic model was not to deal with cyclic loadings but rather to evidence the impact of the local curvature of the yield locus upon macroscopic strain localization. This aspect was examined in detail by [START_REF] Becker | Effect of yield surface curvature on necking and failure in porous plastic solids[END_REF].

The importance of the effect of strain hardening upon ductile failure -and especially the ratcheting of the porosity under cyclic loadings -acts as a strong incentive to develop models accounting better for the heterogeneous distribution of hardening in the plastic matrix. In particular, the incorporation of kinematic hardening appears to be necessary in order to deal with cyclic plasticity [START_REF] Chaboche | On some modifications of kinematic hardening to improve the description of ratcheting effects[END_REF]. The development of a Gurson-type, micromechanically-based model accounting for both isotropic and kinematic hardening, with possibly complex evolution laws for the parameters governing the latter type of hardening, seems of great interest to reproduce the failure of ductile metals under cyclic loadings.

The aim of this paper is to derive such a model. It is organized as follows:

• Section 2 presents the basic ingredients of the theoretical approach and notably the "sequential limit-analysis" approach, which extends the methods and results of classical limit-analysis to materials exhibiting strain hardening.

• Section 3 is devoted to the derivation of a Gurson-like "layer model" extending that developed by [START_REF] Lacroix | Numerical study and theoretical modelling of void growth in porous ductile materials subjected to cyclic loadings[END_REF], accounting for both isotropic and kinematic hardening.

• Section 4 finally compares the predictions of the theoretical model to the results of some new micromechanical finite element simulations.

A C C E P T E D M A N U S C R I P T 2 Position of the problem

In the major part of the paper, no hypothesis is made whatsoever on the magnitude of the displacements and strains, and a large displacement -large strain formulation is used; d represents the local Eulerian strain rate and σ the local Cauchy stress tensor.

Geometry, material, admissibility conditions of velocity fields

In order to derive the overall constitutive law of the porous medium, we consider, following [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I-Yield criteria and flow rules for porous ductile media[END_REF], a spherical "elementary cell" Ω containing a concentric spherical void ω.

The porosity (void volume fraction) is defined by

f = vol(ω) vol(Ω) = a 3 b 3 , ( 1
)
where a is the void's radius and b the cell's external radius. The spherical coordinates and associated local orthonormal basis are denoted r, θ, ϕ and (e r , e θ , e ϕ ).

The material is assumed to be rigid-plastic (no elasticity) and exhibit a mixed, isotropic and kinematic hardening; it is thus supposed to obey the following criterion:

φ(σ(x)) = (σ(x) -α(x)) 2 eq -σ2 (x) ≤ 0, ∀x ∈ Ω -ω, (2)
where σ is the current yield stress and (σα) 2 eq is defined by

(σ -α) 2 eq = 3 2 (σ -α) : (σ -α). (3)
In this expression, σ = σ -1 3 (tr σ) I (where I is the second-order unit tensor) is the deviator of σ, and α is a traceless backstress tensor due to kinematic hardening. The Prandtl-Reuss flow rule associated to the criterion via the normality property reads

d = λ ∂φ ∂σ (σ) = 3 λ(σ -α), (4)
where λ ≥ 0 is the plastic multiplier. The evolution equations of the yield stress σ and

A C C E P T E D M A N U S C R I P T
the backstress α characterizing kinematic hardening will be presented in due time.

The spherical cell is subjected to conditions of homogeneous boundary strain rate:

v(x) = D • x , ∀x ∈ ∂Ω, (5)
where v denotes the local velocity, x the present position-vector and D some overall strain rate tensor.

The velocity must verify the property of incompressibility imposed by the absence of elasticity and the plastic flow rule:

tr d(x) = div v(x) = 0 , ∀x ∈ Ω -ω. (6)

Principles of limit-analysis

Limit-analysis combined with Hill-Mandel homogenization is a convenient framework to derive constitutive equations for porous ductile solids. It permits to effectively operate the scale transition by evidencing the effects of microstructural features at the macroscopic scale.

Classical limit-analysis is limited to rigid-ideal-plastic materials within a small displacement -small strain (linearized) framework. Under such assumptions the macroscopic yield locus can be determined using the classical upper-bound theorem (see e.g. [START_REF] Benzerga | Ductile fracture by void growth to coalescence[END_REF]). The macroscopic stress and strain rate tensors Σ and D being defined as the volume averages of their microscopic counterparts σ and d, the fundamental inequality enunciated by the theorem,

Σ : D ≤ Π(D), (7)
leads to the parametric equation of the yield locus

Σ = ∂Π ∂D (D), (8)
A C C E P T E D M A N U S C R I P T
where D is arbitrary and independent of Σ in inequality ( 7), but tied to it through the macroscopic constitutive law in equation ( 8). The macroscopic plastic potential Π(D) in equations ( 7) and ( 8) is defined by:

Π(D) = inf v∈K(D) (1 -f ) π(d) Ω-ω , (9)
where the notation . Ω-ω stands for volume averaging over the sound volume Ωω.

In this definition the set K(D) consists of those velocity fields v which are kinematically admissible with D and verify the property of incompressibility, and π(d) is the microscopic plastic potential defined for any traceless d by the formula

π(d) = sup σ * ∈C σ * : d, ( 10
)
where C is the microscopic convex domain of reversibility.

Sequential limit-analysis (Yang, 1993;[START_REF] Leu | Analytical and numerical investigation of strain-hardening viscoplastic thick-walled cylinders under internal pressure by using sequential limit analysis[END_REF] heuristically extends the methods and results of classical limit-analysis by incorporating the effects of strain hardening and geometric changes. The idea is, still disregarding elasticity, to consider a hardenable material as the sequence of different, successive rigid-ideal plastic materials. At a given instant, a hardenable material without elasticity behaves, if the hardening and the geometry are considered as fixed, like a rigid-ideal plastic material with some pre-hardening modifying its yield criterion and flow rule. An instantaneous limit-load can thus be determined using the classical limit-analysis theorems. In order to account for changes of the strain hardening and geometry, the local hardening parameters and present positions are then updated approximately using the trial velocity field used in the limit-analysis, integrated in a small time step.

A Gurson-type model accounting for isotropic and kinematic hardening

In this section we will use sequential limit-analysis to derive a Gurson-type model incorporating both isotropic and kinematic hardening. Hardening will be introduced locally in the criterion as a fixed pre-hardening, and the resulting instantaneous loads promoting plastic flow of the cell will be evaluated. The hardening parameters in the matrix will then be considered to evolve a posteriori according to the trial velocity field adopted in

A C C E P T E D M A N U S C R I P T
the limit-analysis.

Macroscopic plastic potential

In order to derive a Gurson-type model accounting for both isotropic and kinematic hardening, we first need to evaluate the macroscopic plastic potential for an initial arbitrary distribution of the quantities characterizing hardening, that is σ(x) and α(x).

Velocity fields

To approximately calculate the plastic potential, we consider [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I-Yield criteria and flow rules for porous ductile media[END_REF]'s trial incom-

pressible velocity field 1 v(x) = b 3 r 2 D m e r + D • x, (11)
where D m =1 3 tr D denotes the mean overall strain rate and D = D-D m I the deviatoric overall strain rate tensor. The associated microscopic strain rate reads

d(x) = D + b 3 r 3 D m (-2e r ⊗ e r + e θ ⊗ e θ + e ϕ ⊗ e ϕ ). (12)

Expression of the approximate plastic potential

Using equations ( 2) and ( 10), the value of the microscopic plastic potential π(d) is calculated to be, for any traceless d,

π(d) = σd eq + α : d, (13)
where the equivalent strain rate d eq is defined by

d eq = 2 3 d : d. ( 14
) A C C E P T E D M A N U S C R I P T
The approximate macroscopic plastic potential thus reads, adopting the trial velocity field defined by equation ( 11) and the associated strain rate d given by equation ( 12):

Π(D) = 1 vol(Ω) Ω-ω π(d) dΩ = 1 vol(Ω) Ω-ω (σd eq + α : d) dΩ = Π iso (D) + Π kine (D), (15)
where the "isotropic" and "kinematic" contributions Π iso (D) and Π kine (D) to Π(D) are

given by

                 Π iso (D) = 1 vol(Ω) Ω-ω σd eq dΩ Π kine (D) = 1 vol(Ω) Ω-ω α : d dΩ. ( 16
)
Some further approximations are needed in order to evaluate these integrals analytically.

One necessary assumption will be that the local hardening parameters are distributed in a certain, specific way within the matrix. The cell is thus considered to be composed of a finite number N of phases distributed in concentric spheres of radii2 (see Figure 1)

a = r 1 < ... < r i < ... < r N +1 = b. ( 17
)
The phase contained within the interval [r i , r i+1 ] is denoted P i . The distribution of the quantity σ characterizing isotropic hardening will be considered as homogeneous in each phase P i . The distribution of the quantity α characterizing kinematic hardening will also be considered as "homogeneous" in each P i , in a more elaborate sense specified below.

It should be noted that this kind of approach is basically similar, in the context of ductile rupture, to that of [START_REF] Herve | N-layered inclusion-based micromechanical modelling[END_REF], in the context of Eshelby's inclusion problem.

A C C E P T E D M A N U S C R I P T a b r i r i+1
Fig. 1. Hollow sphere: definition of some geometric parameters.

Calculation of the isotropic contribution Π iso

The expression ( 16) 1 of the isotropic part Π iso (D) of the macroscopic plastic potential reads

Π iso (D) = 1 vol(Ω) b a 4πr 2 σ d eq (r) S(r) dr, ( 18
)
where the symbol . S(r) denotes an average value over the sphere S(r) of radius r.

In order to evaluate analytically the integral in equation ( 18), we use [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I-Yield criteria and flow rules for porous ductile media[END_REF]'s classical approximation detailed by [START_REF] Benzerga | Ductile fracture by void growth to coalescence[END_REF] and [START_REF] Leblond | Gurson's criterion and its derivation revisited[END_REF]; the potential is evaluated as:

Π iso (D) 1 vol(Ω) b a 4πr 2 σ d 2 eq (r) S(r) dr = 1 vol(Ω) b a 4πr 2 σ D 2 eq + 4b 6 r 6 D 2 m dr. ( 19
)
In order to calculate the potential Π iso (D), we introduce the following approximation on the spatial distribution of the yield limit σ:

A 1 : In each phase P i , the yield limit σ = σi is considered as uniform.

The expression (19) of the isotropic part Π iso (D) of the macroscopic plastic potential can thus be simplified into:

Π iso (D) N i=1 Π iso i (D), ( 20
)
where the "partial" isotropic potential Π iso i (D) in phase P i is given by

Π iso i (D) = 3σ i b 3 r i+1 r i r 2 D 2 eq + 4b 6 r 6 D 2 m dr, ( 21
) A C C E P T E D M A N U S C R I P T
or equivalently upon use of the change of variable u = b 3 /r 3 , by

Π iso i (D) = σi b 3 /r 3 i b 3 /r 3 i+1 D 2 eq + 4D 2 m u 2 du u 2 . ( 22
)
Define now the "local volume fraction" f i as

f i = r i b 3 . ( 23
)
The expression of the partial isotropic potential then becomes

Π iso i (D) = σi 1/f i 1/f i+1 D 2 eq + 4D 2 m u 2 du u 2 = σi   - D 2 eq u 2 + 4D 2 m + 2D m ln   2D m u D eq + 1 + 4D 2 m u 2 D 2 eq     u=1/f i u=1/f i+1 . ( 24
)

Calculation of the kinematic contribution Π kine

In order to calculate the kinematic part Π kine (D) of the potential, we introduce the following approximation on the spatial distribution of the backstress α:

A 2 : In phase P i , the backstress α is of the form:

α = α i = A i 1 + A i 2 (-2e r ⊗ e r + e θ ⊗ e θ + e ϕ ⊗ e ϕ ), (25)
where A i 1 is a second-order traceless tensor and A i 2 a scalar, both uniform within P i . (Note that the term A i 2 (-2e r ⊗ e r + e θ ⊗ e θ + e ϕ ⊗ e ϕ ) is not uniform since the vectors e r , e θ , e ϕ vary within P i ). It should be noted that the form considered in equation ( 25) has been chosen in order to be compatible with the microscopic strain rate given by equation ( 12). This is a reasonable choice since the rate of the backstress depends on the microscopic strain rate, regardless of the model considered for kinematic hardening.

The expression ( 16) 2 of the kinematic potential Π kine (D) then becomes

Π kine (D) = 1 vol(Ω) b a 4πr 2 α : d S(r) dr = 1 vol(Ω) N i=1 r i+1 r i 4πr 2 α i : d S(r) dr. ( 26
) A C C E P T E D M A N U S C R I P T
Let us first study the quantity α i : d S(r) :

α i : d S(r) = A i 1 : D + A i 2 b 3 r 3 D m b : b S(r) where b = -2e r ⊗ e r + e θ ⊗ e θ + e ϕ ⊗ e ϕ . (27)
In this expression the property b S(r) = 0 has been used [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I-Yield criteria and flow rules for porous ductile media[END_REF]. One then gets upon calculation of b : b S(r) :

α i : d S(r) = A i 1 : D + 6A i 2 b 3 r 3 D m . ( 28
)
It follows that

Π kine (D) = N i=1 3 b 3 r i+1 r i A i 1 : D r 2 dr + N i=1 r i+1 r i 18A i 2 D m dr r = A 1 : D + 3A 2 D m , (29)
where

             A 1 = N i=1 A i 1 (f i+1 -f i ) A 2 = N i=1 2A i 2 ln f i+1 f i .
(30)

Macroscopic yield criterion

The macroscopic yield criterion is given by the parametric equation

Σ = ∂(Π iso + Π kine ) ∂D (D) = Σ iso + Σ kine , (31)
where the "isotropic" and "kinematic" contributions Σ iso and Σ kine to the stress Σ are defined by

               Σ iso = ∂Π iso ∂D (D) Σ kine = ∂Π kine ∂D (D). (32) A C C E P T E D M A N U S C R I P T Isotropic contribution.
Since the isotropic contribution to the plastic potential Π iso (D) depends only on D m and D eq , the isotropic stress Σ iso reads

Σ iso = ∂Π iso ∂D m ∂D m ∂D + ∂Π iso ∂D eq ∂D eq ∂D = 1 3 ∂Π iso ∂D m I + ∂Π iso ∂D eq 2D 3D eq = Σ iso m I + Σ iso eq 2D 3D eq , (33)
where the mean and equivalent isotropic contributions Σ iso m and Σ iso eq to the stress are defined by

                 Σ iso m = 1 3 ∂Π iso ∂D m = N i=1 Σ iso m,i Σ iso eq = ∂Π iso ∂D eq = N i=1 Σ iso eq,i (34)
with

                       Σ iso m,i = 1 3 ∂Π iso i ∂D m = 2 3 σi ln 2ξu + 4ξ 2 u 2 + 1 u=1/f i u=1/f i+1 Σ iso eq,i = ∂Π iso i ∂D eq = σi   -4ξ 2 + 1 u 2   u=1/f i u=1/f i+1 , ξ = D m D eq . (35)
Equations ( 34) and ( 35) provide a parametric representation of the quantities Σ iso m and Σ iso eq , ξ acting as a parameter.

Kinematic contribution. The calculation of the kinematic stress Σ kine is straightforward and yields

Σ kine = ∂Π kine ∂D = A 1 + A 2 I. ( 36
)
Summary. By equations ( 31), ( 33) and ( 36), the macroscopic stress inducing plastic flow is given by

Σ -(A 1 + A 2 I) = Σ iso m I + Σ iso eq 2D 3D eq . ( 37
)
Separating the mean and deviatoric parts in this expression and calculating the "von Mises norm" of the latter part, one obtains the macroscopic yield locus in the parametric

A C C E P T E D M A N U S C R I P T form              Σ m -A 2 = Σ iso m (ξ) (Σ -A 1 ) eq = Σ iso eq (ξ), (38)
where Σ iso m , Σ iso eq , A 1 and A 2 are given by equations ( 34), ( 35) and (30).

Flow rule

Since the property of normality of the flow rule is preserved in the homogenization procedure (see [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I-Yield criteria and flow rules for porous ductile media[END_REF]), it reads

D = Λ ∂Φ ∂Σ (Σ) , Λ            = 0 if Φ (Σ) < 0 ≥ 0 if Φ (Σ) = 0, (39)
where D denotes the plastic Eulerian strain rate, Λ the plastic multiplier and Φ(Σ) the macroscopic yield function.

It is worth noting that writing the flow rule (39) explicitly is more difficult than it seems at first sight, since the yield locus is defined by the parametric equations ( 38), the corresponding yield function Φ having no explicit expression. Morin et al. (2015a) have analyzed this problem for a parametric criterion of the form (38) and shown that the flow rule may be rewritten, with a suitable re-definition of the plastic multiplier being now denoted Λ, in the following parametric form:

D = Λ - 1 3 dΣ iso eq dξ (ξ)I + dΣ iso m dξ (ξ) 3 2 Σ -A 1 Σ iso eq (ξ) , ( 40
)
where Σ is the deviator of Σ.

A C C E P T E D M

A N U S C R I P T

Evolution of hardening parameters

Isotropic hardening. Let g 1 denote the function providing the local yield limit σ of the matrix as a function of the local cumulated plastic strain p:

σ = g 1 (p). ( 41
)
We introduce the following approximation on the evolution of isotropic hardening:

A 3 : The "mean yield limit" σi in phase P i is then supposed to depend on some average value pi of the cumulated plastic strain in this phase, through the formula

σi = g 1 (p i ). ( 42
)
The cumulated plastic strain pi is taken on the mid-surface of the spherical shell occupied by phase P i , that is at the radial position ri = 1 2 (r i + r i+1 ); its expression reads

pi = t 0 di eq (τ )dτ, ( 43
)
where di eq is the average equivalent plastic strain rate at the position ri associated to Gurson's trial velocity field:

di eq = d eq S(r i ) ≈ d 2 eq S(r i ) = D 2 eq + 4 b 6 r6 i D 2 m , (44)
where Gurson's classical approximation has again been used.

In particular, for a power-law isotropic hardening, the function g 1 (p) is given by

g 1 (p) = σ0 + hp m , ( 45
)
where σ0 is the initial yield stress, h the isotropic hardening parameter and m the hardening exponent.

A C C E P T E D M

A N U S C R I P T Kinematic hardening. Let g 2 denote the function providing the local rate α of the backstress α as a function of this backstress itself, the strain rate d, the cumulated plastic strain p and the equivalent plastic strain rate d eq : α = g 2 (α, d, p, d eq ). ( 46)

Note that the function g 2 must be homogeneous of degree 1 with respect to the strain rate d since only rate-independent behaviors are considered in the limit-analysis theory.

We introduce the following approximation on the evolution of kinematic hardening:

A 4 : The rate αi of the backstress α i in phase P i is then supposed to depend locally on this backstress itself, the plastic strain rate d i and some average values pi , di eq of the cumulated plastic strain and equivalent plastic strain rate in this phase through the formula

αi = g 2 (α i , d i , pi , di eq ). ( 47
)
We again take the strain rate d i on the mid-surface of the spherical shell defining phase

P i : d i (θ, ϕ) = d(r i , θ, ϕ) ∀(θ, ϕ), (48)
and the average values pi , di eq of the cumulated plastic strain and equivalent plastic strain rate are defined by equations ( 43) and (44), respectively. Equation ( 47) is based on the implicit assumption that the form postulated for αi is compatible with that postulated for α i , equation ( 25). This is not the case for all possible functions g 2 . However the forms ( 25) and ( 47) are compatible for [START_REF] Armstrong | A mathematical representation of the multiaxial Bauschinger effect[END_REF]'s quite (though not fully) general kinematic hardening law,

g 2 (α, d, p, d eq ) = c(p)d -k(p)αd eq , (49)
where c(p) is the "kinematic hardening slope" and k(p) the "strain recovery parameter" (which may both depend on p). For such a law combination of equations ( 25) and ( 47)

A C C E P T E D M A N U S C R I P T
leads to the following evolution equations for the parameters A i 1 and A i 2 :

                 Ȧi 1 = c(p i )D -k(p i )A i 1 D 2 eq + 4 b 6 r6 i D 2 m Ȧi 2 = c(p i ) b 3 r3 i D m -k(p i )A i 2 D 2 eq + 4 b 6 r6 i D 2 m , (50)
where equations ( 12) and ( 44) have been used.

Evolution of geometrical parameters

The evolution of the porosity is classically given by ḟ

= (1 -f ) tr D. (51)
The proposed approximate evolution equation of the internal radii reads

ṙi = ri f i D m . (52)
Note that this equation accounts for the geometry changes due to the sole hydrostatic part of the loading.

3.6 A simplified version of the model and its link with [START_REF] Leblond | An improved Gurson-type model for hardenable ductile metals[END_REF]'s model 3.6.1 [START_REF] Leblond | An improved Gurson-type model for hardenable ductile metals[END_REF]'s model for isotropic hardening [START_REF] Leblond | An improved Gurson-type model for hardenable ductile metals[END_REF] proposed to replace [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I-Yield criteria and flow rules for porous ductile media[END_REF]'s heuristic approach to isotropic hardening effects in porous ductile materials (briefly recalled in Appendix A) with some micromechanical approach based on some approximate analytical solution to the problem of a hollow sphere made of some rigid-hardenable material and subjected to some arbitrary loading. The criterion they obtained was a heuristic extension of that of [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I-Yield criteria and flow rules for porous ductile media[END_REF] involving distinct "average yield stresses of the matrix" in the "square" and "cosh" terms

A C C E P T E D M A N U S C R I P T of the yield function: Σ 2 eq Σ 2 1 + 2f cosh 3 2 Σ m Σ 2 -1 -f 2 = 0, (53)
where Σ 1 and Σ 2 are macroscopic internal variables connected to the spatial distribution of the local yield stress σ. These variables govern the yielding in purely deviatoric and purely hydrostatic loadings, respectively:

• For a purely deviatoric loading (Σ m = 0), the value Σ LPD eq of the overall yield stress is given by

Σ LPD eq = (1 -f )Σ 1 . (54)
• For a purely hydrostatic loading (Σ eq = 0), the value Σ LPD m of the overall yield stress is given by

Σ LPD m = - 2 3 Σ 2 ln f. (55)
In the original version of the model [START_REF] Leblond | An improved Gurson-type model for hardenable ductile metals[END_REF], the expressions of Σ 1 and Σ 2 were obtained in two steps, using the approximate solution of the hollow hardenable sphere problem referred to above. First, Σ 1 and Σ 2 were expressed as functions of the distribution of the local yield stress, using estimates of the overall yield stresses of the hollow sphere under purely deviatoric and purely hydrostatic loadings. Second, the distribution of the local yield stress resulting from the previous mechanical history was obtained in an analytical form under the assumption of positively proportional straining. Although the model did provide a better description of isotropic strain hardening effects in plastic porous materials than that of Gurson recalled in Appendix A, its predictions for cyclic loadings were found to only qualitatively capture the ratcheting of the porosity observed in numerical micromechanical simulations [START_REF] Devaux | On continued void growth in ductile metals subjected to cyclic loadings[END_REF], the quantitative agreement remaining mediocre. The explanation was, of course, the inadequacy of the hypothesis of positively proportional straining in the case of cyclic loadings.

In an improved version of the model [START_REF] Lacroix | Numerical study and theoretical modelling of void growth in porous ductile materials subjected to cyclic loadings[END_REF], the hypothesis of positively proportional straining was dropped. In the absence of any specific hypothesis on the evolution of strain in time, the analytical calculation of the distribution of the local yield stress resulting from the previous mechanical history was no longer possible and numerical integration became necessary. The spherical cell was therefore discretized radiallyexactly like in the layer model developed here which extends the idea so as to incorporate

A C C E P T E D M

A N U S C R I P T kinematic hardening. Using our notations, the expressions of Σ 1 and Σ 2 were given by [START_REF] Lacroix | Numerical study and theoretical modelling of void growth in porous ductile materials subjected to cyclic loadings[END_REF] as

                 Σ 1 = 1 1 -f N i=1 σi (f i+1 -f i ) Σ 2 = - 1 ln f N i=1 σi ln f i+1 f i .
(56)

A simplified version of the model

It is worth noting that in the case of purely isotropic hardening, the limit-loads given by the parametric criterion (38) coincide exactly with the predictions of Leblond et al.

(1995)'s model, for both purely hydrostatic and purely deviatoric loadings. Indeed:

• Hydrostatic loading. In this case, the parameter ξ goes to infinity. The limit-load then reads

Σ iso m (ξ → ∞) = 2 3 N i=1 σi ln f i+1 f i = Σ LPD m . (57)
• Deviatoric loading. In this case, the parameter ξ is equal to zero. The limit-load then reads

Σ iso eq (ξ = 0) = N i=1 σi (f i+1 -f i ) = Σ LPD eq . (58)
A natural idea is then to replace the parametric form given by equations ( 38) for Σ iso m (ξ) and Σ iso eq (ξ) by the simpler explicit yield function proposed by [START_REF] Leblond | An improved Gurson-type model for hardenable ductile metals[END_REF]. We thus propose a simplified version of the model involving the explicit criterion

(Σ iso eq ) 2 Σ 2 1 + 2f cosh 3 2 Σ iso m Σ 2 -1 -f 2 = 0, ( 59
)
or equivalently by equations ( 38):

(Σ -A 1 ) 2 eq Σ 2 1 + 2f cosh 3 2 Σ m -A 2 Σ 2 -1 -f 2 = 0, ( 60
)
where Σ 1 , Σ 2 , A 1 and A 2 are given by equations ( 56) and (30), respectively.

It is worth noting that this simplified form coincides with the full parametric form in the case of a homogeneous distribution of the isotropic hardening parameter σ but a

A C C E P T E D M A N U S C R I P T
heterogeneous distribution of the kinematic hardening parameter α; indeed in this case the parametric expressions (38) of Σ iso m (ξ) and Σ iso eq (ξ) exactly yield, upon elimination of ξ, the criterion (59), with Σ 1 = Σ 2 = σ.

Assessment of the model

The model will now be assessed by comparing its predictions to the results of some micromechanical finite element simulations of the elementary cell considered in its derivation.

Yield surfaces with initial pre-hardening

Description of the simulations

First, in order to study the macroscopic criterion, we consider all internal parameters (geometry and hardening) as fixed. We thus solve the limit-analysis problem for a given, fixed pre-hardening (resulting from some given prestraining), without any geometry update, using the finite element method (FEM).

The calculations are performed with a homemade code using 2D axisymmetric meshes subjected to conditions of homogeneous boundary strain. Figure 2 shows the mesh used for a void volume fraction f = 0.01. Eight-node quadratic elements subintegrated with 2×2 Gauss points are used. The mesh contains 3,540 elements and 10,859 nodes (21,718 degrees of freedom). This discretization is adequate for the numerical calculations envisaged, further mesh refinement making no appreciable difference to the results. Axisymmetric loadings are considered: Σ 11 = Σ 22 = 0, Σ 33 = 0, and Σ ij = 0 otherwise. Two Lode angles (denoted θ L ) are considered: θ L = 0 corresponding to Σ 33 -Σ 11 > 0 and θ L = π corresponding to Σ 33 -Σ 11 < 0. The simulations are performed by solving an elasticplastic evolution problem, the limit-load being considered as reached when the overall stress components no longer evolve [START_REF] Michel | Effective properties of composite materials with periodic microstructure: a computational approach[END_REF].

A C C E P T E D M

A N U S C R I P T Fig. 2. Spherical mesh used in the finite element calculations (f = 0.01).

Isotropic pre-hardening

We consider two cases with isotropic pre-hardening (and no kinematic pre-hardening: α = 0 everywhere in the matrix):

• Case 1. Hardening is assumed to be more important near the cavity; the yield stress is supposed to vary linearly with r from the value σ(r = a) = 1.5 σ0 to the value σ(r = b) = 0.5 σ0 .

• Case 2. Hardening is assumed to be more important near the cell's boundary; the yield stress is supposed to vary linearly with r from the value σ(r = a) = 0.5 σ0 to the value σ(r = b) = 1.5 σ0 . Note that this case is probably unrealistic in practice since hardening will likely concentrate near the voids's boundary; it is however interesting to consider it in order to study the robustness of the model. (1977)'s model does not account for the different values of the macroscopic yield stresses in the "square" and "cosh" terms of the yield criterion (resulting from the heterogeneous distribution of the local yield stress within the matrix), two values of Gurson's unique "overall yield stress" Σ (see Appendix A) are envisaged, Σ = Σ 1 and Σ = Σ 2 , these quantities being given by equation ( 56). Some comments are in order here: • In the two cases considered, the model developed is in very good agreement with the finite element results for all the triaxialities considered; in particular, the hydrostatic point is perfectly reproduced.

A C C E P T E D M A N U S C R I P T 0 0.5 1 1.5 2 2.5 3 3.5 Σ m /σ 0 -1 -0.5 0 0.5 1 (Σ 33 -Σ 11 )/σ 0 Present model Simplified model Gurson, Σ = Σ 1 Gurson, Σ = Σ 2 FEM (a) 0 0.5 1 1.5 2 2.5 3 3.5 Σ m /σ 0 -1 -0.5 0 0.5 1 (Σ 33 -Σ 11 )/σ 0 Present model Simplified model Gurson, Σ = Σ 1 Gurson, Σ = Σ 2 FEM (b)
• In both cases, the full and simplified models yield very similar results: the simplified model with an explicit form of the yield criterion thus seems to be a viable alternative to the more complex one.

• The comparison with [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I-Yield criteria and flow rules for porous ductile media[END_REF]'s model highlights the influence of the heterogeneous distribution of pre-hardening on the strength of the porous material: this model involving a single overall yield limit Σ can capture either the hydrostatic point or the deviatoric point, depending on the choice made for this overall yield limit, but not both.4

• Finally, it is worth noting that the finite element results are slightly sensitive to the Lode angle θ L (small dissymmetry of the numerical results with respect to the horizontal axis), in contrast to the model developed that does not account for a Lode angle dependency;

this is due to [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I-Yield criteria and flow rules for porous ductile media[END_REF]'s approximation (leading from equation ( 18) to equation ( 19) above) that erases the Lode angle dependency in the macroscopic plastic potential [START_REF] Cazacu | On the combined effect of pressure and third invariant on yielding of porous solids with von Mises matrix[END_REF][START_REF] Leblond | Gurson's criterion and its derivation revisited[END_REF].

A C C E P T E D M

A N U S C R I P T 4.1.3 Kinematic pre-hardening

We now consider three cases with kinematic pre-hardening (and no isotropic pre-hardening: σ = σ0 = Cst. everywhere). In all cases the local backstress α is of the form

α = α 1 (r) + α 2 (r)(-2e r ⊗ e r + e θ ⊗ e θ + e ϕ ⊗ e ϕ ), (61)
where the traceless tensor α 1 (r) and the scalar α 2 (r) depend only on r.

• Case 1. Pre-hardening is assumed to affect the sole deviatoric stress; the non-zero components of α 1 are: α 1(11) = α 1(22) = -α 1(33) /2 = -σ 0 /6 and α 2 is nil.

• Case 2. Pre-hardening is assumed to affect the sole hydrostatic stress; α 2 is supposed to vary linearly from the value α 2 (r = a) = -σ 0 /6 to the value α 2 (r = b) = 0, and α 1 is nil.

• Case 3. Pre-hardening is assumed to affect both the hydrostatic and deviatoric stresses; everywhere in the matrix),5 and the finite element results, for a porosity f = 0.01. (In this case, the simplified model exactly coincides with the more complete one since there is no isotropic pre-hardening, so its predictions need not be shown).

In all three cases considered, the model developed is in very good agreement with the finite element results for all the triaxialities considered; in particular, the model reproduces the "translatory motion" of the yield surface perfectly, in contrast to [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I-Yield criteria and flow rules for porous ductile media[END_REF]'s model without kinematic hardening. As expected from the theoretical form of the criterion, the term A 1 affects only the position of the yield surface in the "deviatoric direction" while the term A 2 affects only its position in the "hydrostatic direction".

Evolution problems

Description of the simulations

We now wish to compare the model's predictions to the results of micromechanical finite element simulations of evolution problems. The hardening is now supposed to evolve locally within the matrix but geometry changes are still disregarded. The motivation for this choice is to comply with the model's hypothesis of spherical voids, in order to assess the sole evolution of the hardening; the modification of the geometry would induce some void shape effects that are disregarded by the model and thus would introduce some discrepancies.

A C C E P T E D M A N U S C R I P T

Again, the simulations are performed using 2D axisymmetric meshes subjected to conditions of homogeneous boundary strain (see Figure 2). The macroscopic stress tensor is constrained to be of the form

Σ = C(t)Σ 0 , Σ 0 =            Σ 0 11 0 0 0 Σ 0 11 0 0 0 Σ 0 33            . (62)
In this expression, C(t) is a time-varying scalar and Σ 0 is a constant tensor whose components Σ 0 11 and Σ 0 33 are given by

Σ 0 11 = cosβ - 1 3 sinβ , Σ 0 33 = cosβ + 2 3 sinβ, ( 63
)
where the angle β depends on the stress triaxiality T :

β = argtan 1 T . (64)
Only the value θ L = 0 of the Lode angle is considered, corresponding to Σ 33 -Σ 11 > 0.

The simulations are performed incrementally, for proportional loading paths -implying a constant stress triaxiality [START_REF] Michel | Effective properties of composite materials with periodic microstructure: a computational approach[END_REF]. Finally the macroscopic strain tensor is denoted E.

We investigate the cases of isotropic hardening and linear kinematic hardening. Isotropic hardening is supposed to follow the power-law defined by equation ( 45), and kinematic hardening the law defined by equation ( 49) with c = Cst. and k = 0 (no strain recovery).

The values of the material parameters considered are as follows: Young's modulus, E = 10, 000 GPa, Poisson's ratio, ν = 0.25; isotropic hardening parameters, σ0 = 100 MPa, h = 400 MPa and m = 0.2; kinematic hardening slope, c = 400 MPa. Note that the value chosen for Young's modulus is extremely high so that elasticity is negligible in the simulations, in agreement with the hypothesis made in the derivation of the model. The distribution of the local equivalent von Mises stress is represented in Figure 6 for the present model and the finite element simulations. The present model globally accurately reproduces the heterogeneous distribution of hardening. In particular, the finite element results reveal that the hypothesis made in the model of an essentially radial variation of the yield limit σ is acceptable. In the case of a pure hydrostatic loading, T = +∞, this approximation becomes exact. In the cases of low or high triaxiality, T = 1/3 and 2.333, the maps show that the hypothesis of essentially radial dependence is acceptable in a large domain far from the void's boundary, and that the stress level predicted is correct.

The hypothesis is no longer verified near the void's boundary, but the angular average of the finite element values is qualitatively well predicted by the model.

Kinematic hardening

Figure 7 compares the strain-stress curve predicted by the present model with N = 30 phases to the results obtained by the finite element method.

In this case, the predictions of the model coincide almost perfectly with the finite element results, emphasizing that the "macroscopic backstresses" A 1 and A 2 are sufficient to capture the effect of the heterogeneous distribution of kinematic hardening. Again it is worth noting that the geometry is not allowed to evolve in the simulations, which ensures that void shape effects are disregarded. It should be noted that the very good results observed for yield surfaces as well as evolution problems have been achieved with a reasonable number of layers (N = 10 for isotropic hardening and N = 30 for kinematic hardening). These discretizations were found to be sufficient to accurately describe both the macroscopic results and the distributions of microscopic hardening, further refinement making no appreciable difference to the results.

It should be noted too that the model developed, with the number of layers considered, does not require a significantly larger CPU time than Gurson's original model. This is due to the fact that the integration over the layers represents only a small number of operations in the local projection algorithm, representing itself a modest part of the global iterative algorithm. However, the new model requires more memory than Gurson's original model in order to store the internal variables in the layers.

Conclusion

The aim of this paper was to develop a model for ductile porous materials accounting for both isotropic and kinematic hardening. This model differs (although it is inspired) from that of [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I-Yield criteria and flow rules for porous ductile media[END_REF] including only isotropic hardening, and also from its extension to kinematic hardening due to [START_REF] Mear | Influence of yield surface curvature on flow localization in dilatant plasticity[END_REF], in that it is no longer based on a purely heuristic and macroscopic approach, but on some detailed analysis of the effect of the heterogeneous distribution of microscopic hardening parameters near the voids.

An approximate yield criterion was derived by performing a "sequential limit-analysis" of a hollow sphere made of a rigid-hardenable matrix. To approximately account for the heterogeneity of hardening, the cell was discretized into a finite number of spherically distributed phases in which the quantities characterizing hardening were considered as The model was assessed numerically using micromechanical finite element simulations.

First, overall yield loci were investigated for both isotropic and kinematic pre-hardening: a very good agreement was observed between the numerical results and the predictions of the model. Then these predictions were assessed on evolution problems: they were found to be in very good agreement with the finite element results, with regard to both the overall stress-strain curves and the distributions of isotropic and kinematic hardening parameters.

Future developments of the work will include:

• Numerical micromechanical simulations for cyclic loadings including geometry changes, with some comparisons with the predictions of the present model. Such studies are desirable in order to understand the effect of mixed isotropic/kinematic hardening on the ratcheting of the porosity mentioned in the Introduction.

• Finite element implementation of the model and applications to numerical studies of the ductile rupture of actual test specimens subjected to cyclic loadings.

• Derivation of a model including void shape effects, in order to investigate cyclic ductile rupture under conditions of low stress triaxiality.

Figure 3

 3 Figure 3 compares the yield surfaces associated to the theoretical model of Section 3.2 and its simplified version of Section 3.6 with N = 10 phases 3 , Gurson (1977)'s model without pre-hardening and the finite element results, for a porosity f = 0.01. Since Gurson

Fig. 3 .

 3 Fig. 3. Yield surfaces for isotropic pre-hardening: theoretical model given by equations (38) (Present model), simplified model given by equations (60) (Simplified model), Gurson's model (Gurson) and finite element results (FEM). (a) Case 1, (b) Case 2 (see text).

  the nonzero components of α 1 are α 1(11) = α 1(22) = -α 1(33) /2 = -σ 0 /4, and α 2 is supposed to vary linearly from the value α 2 (r = a) = -σ 0 /4 to the value α 2 (r = b) = 0.

Figure 4

 4 Figure 4 compares the yield surfaces associated to the theoretical model developed in Section 3.2 with N = 30 phases, Gurson (1977)'s model without pre-hardening (with σ = σ0

Fig. 4 .

 4 Fig. 4. Yield surfaces for kinematic pre-hardening: theoretical model given by equations (38) (Present model), Gurson's model without hardening (Gurson) and finite element results (FEM). (a) Case 1, (b) Case 2, (c) Case 3 (see text).

Fig. 5 .

 5 Fig. 5. Stress-strain curve in the case of isotropic hardening: predictions of the present model (Present model), Gurson (1977)'s model with a heuristic modelling of isotropic hardening (Gurson) and finite element results (FEM).

Fig. 7 .

 7 Fig. 7. Stress-strain curve in the case of kinematic hardening: predictions of the present model (Present model) and finite element results (FEM).

Fig. 8 .

 8 Fig. 8. Distribution of the local equivalent backstress α eq at E : Σ 0 = 0.05 in the case of kinematic hardening. (a) T = 1/3, (b) T = 2.333, (c) T = +∞. Top: finite element results, bottom: predictions of the present model.

  homogeneous. The macroscopic yield locus was characterized by an overall criterion expressed in a parametric form, wherein the heterogeneous local hardening parameters were accounted for through macroscopic variables. A simplified version of the model involving A C C E P T E D M A N U S C R I P T an explicit form of the overall criterion was then proposed. This version reduces to Lacroix et al. (2016)'s model in the case of purely isotropic hardening.

The component proportional to D m of this velocity field reproduces the exact solution for a hollow sphere made of plastic material (with or without hardening) and subjected to a hydrostatic loading.

The thickness of the phases is not necessarily uniform and may be chosen arbitrarily.

In all the simulations, the thickness of the phases has been chosen uniform for simplicity.

The reason why[START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I-Yield criteria and flow rules for porous ductile media[END_REF]'s model fails so utterly to reproduce the numerical results is that the variation of σ considered here in the matrix is very important; for a more moderate variation the results would have been more acceptable.

[START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I-Yield criteria and flow rules for porous ductile media[END_REF]'s criterion, which does not account for kinematic hardening, is represented with the sole purpose of evidencing the influence of pre-hardening.

Acknowledgments

Fruitful discussions with P. Suquet and D. Kondo are gratefully acknowledged.

A C C E P T E D M A N U S C R I P T

The response is investigated for three stress states: (i) pure hydrostatic loading, T = +∞;

(ii) high triaxiality, T = 2.333; (iii) low triaxiality, T = 1/3. We are interested in:

(1) Macroscopic quantities. We study the evolution of the stress amplitude C(t) versus the "measure of deformation" E : Σ 0 .

(2) Microscopic quantities -isotropic hardening. For this type of hardening we study the distribution of the local equivalent von Mises stress σ eq (at the end of the simulation for which E : Σ 0 = 0.05). This variable reduces to the local yield stress σ in the model which assumes that the matrix is plastic everywhere; but this is not true in the finite element calculations since the material is not necessarily entirely plastic.

(3) Microscopic quantities -kinematic hardening. For that type of hardening we study the distribution of the local equivalent backstress α eq = 3 2 α : α (at the end of the simulation for which E : Σ 0 = 0.05). The present model's predictions are globally very close to the numerical results in the three cases considered. In the intermediate case (T = 2.333) the model slightly overestimate the finite element results. It is worth noting that the very good agreement observed at very low triaxiality (T = 1/3) is due in part to the fact that the geometry is not allowed to evolve in the simulations, which ensures that the void remains spherical. [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I-Yield criteria and flow rules for porous ductile media[END_REF]'s predictions are in good agreement with the numerical results for low and moderate triaxialities but become poor at very high triaxiality. This highlights the detrimental effect of Gurson's hypothesis that the overall yield stresses under purely hydrostatic and purely deviatoric loadings may be related to a single average yield stress of the matrix. is acceptable; the comments made above for isotropic hardening also apply to kinematic hardening. Note, however, that for kinematic hardening, unlike for isotropic hardening, formula (25) predicts some dependence of α upon the spherical angle θ, which is confirmed by the finite element results, see case (b).

Isotropic hardening

A C C E P T E D M A N U S C R I P T

A Appendix: Gurson's heuristic approach to isotropic hardening effects in porous ductile materials [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I-Yield criteria and flow rules for porous ductile media[END_REF]'s original model was obtained for a rigid-ideal plastic material obeying von Mises's criterion and the associated Prandtl-Reuss flow rule, from an approximate limit-analysis of a hollow sphere made of such a material and subjected to conditions of homogeneous boundary strain rate. The macroscopic yield criterion obtained read

where σ0 is the yield stress, uniform in the matrix.

To introduce isotropic hardening effects into the model, instead of pursuing his micromechanical analysis (as was later done by [START_REF] Leblond | An improved Gurson-type model for hardenable ductile metals[END_REF]), Gurson dropped homogenization and limit-analysis and adopted a heuristic approach. This approach consisted in retaining the analytic form (A.1) of the criterion, replacing σ0 with some "average yield stress" Σ of the matrix given by: Σ = g 1 (P ); (A.2) in this equation g 1 (p) denotes the function providing the local yield limit σ as a function of the local cumulated plastic strain p, like in Subsection 3.4 above, and P represents some "average equivalent cumulated strain" in the heterogeneous, porous material. The evolution of P was assumed to be governed by the following equation:

(

which expressed the heuristic assumption of equality of the plastic dissipations in the real, heterogeneous porous material and in a fictitious "equivalent" homogeneous material with equivalent cumulated strain P and yield stress Σ.