text:
Gérald Portemont

Julien Berthe
email: julien.berthe@onera.fr

Alain Deudon

François-Xavier Irisarri

Static and dynamic bearing failure of carbon/epoxy composite joints

Keywords: CFRP, Dynamic testing, Bearing failure, Bolted joints

Mechanical fastening is a common method used to join composite materials in aeronautical industry. Various studies have been performed dedicated to the behaviour of composite bolted joints under quasi-static loadings, but only few studies deal with the dynamic behaviour (crash or impacts). The aim of this work is to study the loading rate influence on the bearing response of a carbon/epoxy laminate loaded by a pin. For that purpose, a double shear test fixture has been specially designed to measure the global behaviour and the local response around the pin. Infrared thermography and Digital Image Correlation techniques have been used to detect, map and characterize dissipative phenomena evolution. The tests have been performed on a servohydraulic jack with a loading rate ranging from 10 -4 m/s to 1 m/s. An increase of the peak bearing load of more than 20% is observed with the loading rate increase. A decrease of the load plateau of more than 60% is obtained. Simultaneous measurements of thermal and kinematic fields in this work give access to the evolution of the damage-related dissipative

Introduction

Over the last decades, Carbon Fibre Reinforced Polymers (CFRPs) have been increasingly used in transportation industry, especially in aeronautics, in order to reduce the weight of vehicles. CFRPs represent more than 50% of the mass of the structure of the most recent Airbus and Boeing airliners.

Thus, joining composite components to other composite or metallic parts is of utmost importance in the design process of such structures. Mechanical fastening is the most widespread joining method used in aeronautics for composite materials, due to the need of disassembly for maintenance tasks and the lack of reliable non-destructive control methods for bonded joints.

Mechanical fastening of composite laminates has been thoroughly examined in the literature since the 70's. Review articles on the topic can be found in [START_REF] Camanho | Stress analysis and strength prediction of mechanically fastened joints in FRP: a review[END_REF] or [START_REF] Thoppul | Mechanics of mechanically fastened joints in polymermatrix composite structures A review[END_REF]. First experimental studies of the composite bolted-joint behavior focused on the characterization of the main macroscopic failure modes and influencing parameters on the joint behavior (see for instance [START_REF] Collings | On the bearing strengths of CFRP laminates[END_REF][START_REF] Collings | Bearing deflection behaviour of a loaded hole in CFRP[END_REF] or [START_REF] Lj | Bolted Joints In Graphite-Epoxy Composites[END_REF]). In-plane failure of the joint mostly occurs by tension, shear or bearing failure, or combined failures modes. With laminates, tension and shear are catastrophic failure modes whereas bearing failure is progressive.

Bearing failure develops for high width-to-hole-diameter and end-distanceto-hole-diameter ratios. Bearing strength is mostly influenced by the lay-up, the hole-diameter-to-laminate-thickness ratio, the bolt-hole clearance and the fastening technology. In particular, bolt lateral clamping has been shown to significantly increase the laminate quasi-static bearing strength, up to 100%, by preventing out-of-plane damage effects under the washers [START_REF] Park | Bearing failure analysis of mechanically fastened joints in composite laminates[END_REF][START_REF] Kelly | Bearing strength of carbon fibre/epoxy laminates: effects of bolt-hole clearance[END_REF][START_REF] Khashaba | Effect of washer size and tightening torque on the performance of bolted joints in composite structures[END_REF][START_REF] Irisarri | Experiments and modeling of clamping effects on the bearing strength of mechanically fastened joint in CFRP laminates[END_REF]. In complex junctions involving more than a single row of fasteners, the joint failure is also affected by the unequal load distribution between the fasteners [START_REF] Tate | Preliminary investigation of the loads carried by individual bolts in bolted joints[END_REF][START_REF] Ross | An Electrical Computer for the Solution of Shear-Lag and Bolted-Joint Problems[END_REF][START_REF] Mccarthy | A simple method for determining the effects of bolthole clearance on load distribution in single-column multi-bolt composite joints[END_REF][START_REF] Ekh | Finite element modeling and optimization of load transfer in multi-fastener joints using structural elements[END_REF][START_REF] Gray | A global bolted joint model for finite element analysis of load distributions in multi-bolt composite joints[END_REF][START_REF] Lecomte | An analytical model for the prediction of load distribution in multi-bolt composite joints including hole-location errors[END_REF]. The load distribution between fasteners is itself influenced by the stiffness of the single-bolt elementary joint in the loading direction. Consequence of the unequal load distribution in a line of fasteners is that most elementary joints are submitted to a combination of bypass and bearing loadings that triggers a competition between tensile failure mode and bearing failure mode [START_REF] Crews | Bearing-bypass loading on bolted composite joints[END_REF]. Fractographic analyses of the bearing failure in laminates show that bearing is due to the accumulation of local damages on the loaded side of the hole. In particular, bearing failure is driven by both fiber-kinking in the load-sustaining plies parallel to the loading direction and delamination between the plies [START_REF] Kelly | Bearing strength of carbon fibre/epoxy laminates: effects of bolt-hole clearance[END_REF][START_REF] Wu | Modeling bearing failure initiation in pin-contact of composite laminates[END_REF][START_REF] Ireman | On damage development in mechanically fastened composite laminates[END_REF][START_REF] Xiao | Bearing strength and failure behavior of bolted composite joints (part I: Experimental investigation)[END_REF][START_REF] Irisarri | Progressive damage and failure of mechanically fastened joints in CFRP laminates Part I: Refined Finite Element modelling of single-fastener joints[END_REF][START_REF] Sola | On the role of kinking in the bearing failure of composite laminates[END_REF].

The aforementioned studies, deal with quasi-static testing and modeling of the composite joints. However, energy absorption capability of the bolted joint is an important factor for structural crashworthiness. The energy dissipated by one bolted joint is small, but considering the huge number of joints, the sum is no longer negligible at the scale of a complete airliner. Few studies have been published about the influence of the loading rate on the global behavior of the joint. Ger et al. [START_REF] Ger | Dynamic tensile strength of composite laminate joints fastened mechanically[END_REF] compared quasi-static and dynamic response (loading rates of 3 to 5 m/s) of CFRP and hybrid carbon and kevlar fiber reinforced plastics for various joint configurations. The authors concluded that energy absorption as well as ultimate load decreased with increasing loading rates. However, Li et al. [START_REF] Li | Static and dynamic behaviour of composite riveted joints in tension[END_REF] pointed out that the inertia effect had not been considered to calculate the dynamic failure load in [START_REF] Ger | Dynamic tensile strength of composite laminate joints fastened mechanically[END_REF] which could lead to unreliable conclusions. Several CFRP joint configurations with two rivets were tested at loading rates from 0.5 mm/min to 8 m/s in [START_REF] Li | Static and dynamic behaviour of composite riveted joints in tension[END_REF]. The authors concluded that increasing loading rate had no significant influence on the joint stiffness and strength, but observed an increase of the total energy absorption of the single-lap joint, mostly explained by a change of failure mode from net-tension to bearing for the highest loading rates. Coherent results were obtained by Heimbs et al. [START_REF] Heimbs | Static and dynamic failure behaviour of bolted joints in carbon fibre composites[END_REF] who tested CFRP singlelap countersunk bolted joints up to 10 m/s. Egan et al. [START_REF] Egan | Static and high-rate loading of single and multi-bolt carbonepoxy aircraft fuselage joints[END_REF] investigated single-lap joint behavior of CFRP bolted joints at quasi-static and dynamic (5 m/s and 10 m/s) loading speeds. It was shown that energy absorption of the joint can be maximized by delaying the transition from progressive bearing failure to ultimate failure. Increasing the loading rate tends to favor ultimate failure by fastener pull-through over fastener failure, thus allowing more extensive bearing damage and energy absorption. Pearce et al. [START_REF] Pearce | Experimental Investigation of Dynamically Loaded Bolted Joints in Carbon Fibre Composite Structures[END_REF] studied the response of single fastener single-lap countersunk bolted joints with in-plane loading (bearing) or normal loading (pull-through). Loading rates vary from 0.1 m/s to 10 m/s. Only minor loading rate dependency was observed in pull-through. For the specimen loaded in bearing, a transition was observed from fastener failure at low loading rates to pull-through at 1 m/s and above. The energy absorption capability of bearing mode failure was specifically investigated by Heimbs and Bergmann [START_REF] Heimbs | Bearing Mode Absorber On the Energy Absorption Capability of Pulling a Bolt through a Composite or Sandwich Plate[END_REF]. A test set-up was devised to continuously pull a bolt through the laminated plate in a configuration somewhat similar to a double-lap pinned joint. Energy absorption was shown to decrease significantly with increasing loading rates in CFRP laminates.

The literature review clearly exhibits that studies related to the loading rate dependencies are classically performed on joints close to the technological solutions. In such cases, the understanding of the rate dependency is complex because of the interaction between the rivet and the composite plate which leads to a competition between bearing of the composite plate and pullthrough of the rivet. The present work is focused on the study of the loading rate dependency of the bearing response of composite plates. In section 2, a specific experimental set-up on a servo-hydraulic jack has been devised on the principle of double shear lap test in order to perform global (load and displacement) and local measurements (digital image correlation and infrared thermography close to the pin) of the bearing response of the composite plate. Experimental results are presented and discussed in Section 3. Tests are performed at different loading rates from 10 -4 m/s to 1 m/s. Global and local measurement are analysed in order to evaluate the influence of the loading rate on the bearing response.

Material and methods

Preparation of the specimens

The composite plates used in this study have been manufactured from Hexply M21/35%/268/T700GC (Hexcel France) unidirectional carbon fibre prepreg plies based on a third generation epoxy resin system. The laminated plates have been cured under specific cure cycle in a hydraulic press, after hand lay-up. A generic 16-plies quasi-isotropic stacking sequence [(45/90/-45/0) 2 ] s is used in this study (thickness 4.02mm). The specimens have been obtained after two machining operations (cutting and milling with appropriate tools) with the following manufacturing tolerance for the holes: 6.35 +0.02 -0.00 . The characterisation of the bearing response of the CFRP laminate under quasi-static loading is normalised by the ASTM D5961 procedure [START_REF]Standard test method for bearing response of polymer matrix composite laminates[END_REF]. The specimen geometry and dimensions used in the present work are shown in Figure 1. Because of the holder geometry of the testing machine, two notches have to be machined in the clamped area of the specimen. With this specimen design the value of the geometrical ratios influencing the behaviour of the assembly are: W/D=5.6 and E/D=4.7 with D the hole diameter, W the width of the sample and E the distance between the hole and the end of the specimen. Preliminary quasi-static tests have been performed to validate the bearing failure of the selected geometry.

Experimental set-up

The present work is focused on the analysis of the loading rate effect on the bearing failure of the laminate. According to the ASTM D5961 standard [START_REF]Standard test method for bearing response of polymer matrix composite laminates[END_REF], the bearing response of the laminate can be studied using either a singlelap configuration or a double-lap test. However, in the case of the single-lap test, the bending effect due to the excentricity of the load path generates out-of-plane effects that can lead to fastener pull-through (see for example [START_REF] Egan | Static and high-rate loading of single and multi-bolt carbonepoxy aircraft fuselage joints[END_REF]). To prevent this unwanted failure mode, a double-lap configuration has been chosen in the following.

A specific double shear test fixture for hydraulic jack has been developed for this study. As it can be seen in the exploded diagram of the Figure 2, the test fixture brings several improvements with respect to the standard setup. Metallic plates (numbered 2 and 3) with high elastic limit (Ramax 1400 MPa) are used to transfer the load to a metallic pin. Observation windows have been machined into the plates to enable observation of the specimen surface in front of the pin in the loading direction. The distance between the two plates is controlled by the two spacers (numbered 5 and 6) which ensure sufficient space for the specimen with a clearance of approximatively 0.2 mm (pre-stress is not considered in this study). The spacers also increase the overall stiffness of the set-up and prevent out-of-plane bending of the plates. The two loading plates are linked to the hydraulic jack through the part number 1 in Figure 2. To avoid relative displacement between the two loading plates during the loading, grooves have been machined into the plates.

Finally, a metallic pin with a specific thermal treatment (cementation) is used to load the specimen. The test fixture has been used with a hydraulic jack to perform displacement controled tests with a loading rate varying from 5 mm.min -1 to 1 m.s -1 .

Global and local measurements

Global and local measurements have been performed to characterise the bearing response of the composite plate under dynamic loading. The experimental set-up is shown in Figure 3. A ±200 kN piezoelectric load cell (Kistler) is used to measure the applied load on the specimen. The load cell is localised between the lower holder and the hydraulic jack support. The pin displacement is measured by means of different systems depending of the loading rate. Between 5 mm.min -1 and 0.1 m.s -1 , a Keyence laser sensor is used to monitor the displacement of the pin end. For higher loading speeds, a Rudoph optical extensometer is used. The measurement is based on the tracking of a black-and-white target with a resolution of 10 -4 of the measuring range, which leads here to a resolution of 2 µm. To measure the strain in the loading direction, a strain gauge (TML YFLA-2) is glued on each specimen at a distance of 25 mm of the hole boundary, between the hole and the fixed jaw. The strain gauge is conditioned with a VISHAY 2310 conditioner with 75 kHz cutting frequency.

Using the observation windows (see Figure 2), displacement and thermal field measurements have been performed on the surfaces of the specimen in the purpose of detecting in the vicinity of the hole the damage phenomena related to the progressive bearing failure of the laminate. The kinematic field are measured using 2D Digital Image Correlation (DIC) technique. The images are captured by a Photron high speed camera (Photron SA-X). In this study, two configurations of the camera have been used: 1024x1024 pixels with a frequency of 10 kHz for loading speeds above 0.1 m.s -1 and 512x504 pixels with a frequency of 40 kHz for higher loading speeds. The black and white speckle pattern applied on the specimen is lightened with two high power LED projectors (150 W). That choice allows to reduce the exposure time in order to limit blurred area in the images. Additionally, that light technology does not influence the thermal field measurement performed on the other side of the sample [START_REF] Crump | Measuring visco-elastic behaviour of composites at high strain rates using optical techniques[END_REF]. To obtain the displacement field, images are divided into small areas called zone of interests (ZOI) and characterised by their grey level distribution. This characteristic distribution is used by the algorithm to find the ZOI in the deformed image and to build the displacement field. The computation of the displacement fields and strain fields are performed with the ARAMIS 6.1 software (GOM). For each loading speed, a ZOI of 15 x 15 pixels and a distance between the centres of two successive ZOIs of 13 x 13 pixels have been employed for DIC computations.

An infrared camera is used to measure the evolution of the thermal field in the hole area during dynamic loading. The camera used is a CEDIP JADE3 LWIR with 320x240 MCT detectors sensitive to the long infrared wavelength (between 8 and 12 µm) equipped with a 50 mm LW objective.

To obtain a quantitative temperature measurement the camera has to be calibrated. The calibration is performed by putting a blackbody in front of the camera. Images are taken for different surface temperature of the blackbody and a global response law of the matrix is built to link the digital level measured to the temperature. In this study, the objective is to measure small temperature gradients at the surface of the specimen over few pixels of the IR camera. Indeed, the first bearing damages are localised in a small area near the loaded boundary of the hole. Some authors have shown that in such a case an alternative calibration method is required to increase the accuracy of the temperature measurement [START_REF] Poncelet | htrognits intrinsques et structurales des essais d'auto-chauffement et de fatigue grand nombre de cycles[END_REF][START_REF] Bodelot | Experimental setup for fully coupled kinematic and thermal measurements at the microstructure scale of an aisi 316l steel[END_REF]. In this alternative method the calibration law of the matrix is replaced by a calibration law for each pixel. In this study, the pixel-by-pixel calibration procedure has been performed in the To synchronise all measurement systems and to record data a 1-MHz data acquisition system has been used (NICOLET).

Results

To study in detail the loading rate dependency, tests have been performed for five different speeds: 8.3•10 -5 m.s -1 , 1•10 -3 m.s -1 , 8.3•10 -3 m.s -1 , 0.1 m.s -1 and 1 m.s -1 . As shown in Table 1, for each speed at least three tests have been performed to check repeatability. In the following, the experimental results are analysed in three phases. First, the consistency of the different load and displacements measurements systems employed is checked using a quasi-static test. Second, the influence of the loading rates on the global bearing response of the specimens is analysed. Finally, a local analysis is performed thanks to digital image correlation and infrared thermography.

The objective of this analysis is to gain insight into the influence of the loading rate on the material behaviour near the hole.

Consistency of the load and displacement measurement systems

As it is shown in Figure 4, for each tests the applied load (blue curve), the pin displacement (red curve) and the strain measured by the strain gauge (green curve) are recorded. Concerning the pin displacement, two different measuring system are used depending of the loading rate. A comparison of the displacement obtained by the two means have been performed on quasistatic test (8.3•10 -5 m.s -1 ) and a difference of 10 µm is obtained for the value of the displacement at the bearing peak, which corresponds to a variation of 2.5% of the bearing peak value.

The Figure 4 also shows that the strain curve and the load curve have the same shape. Apparently, the strain measurement is not disturb by the perforation as it is performed at a sufficient distance of the hole. Therefore, the applied load can also be calculated with the gauge measurement knowing the elastic properties of the laminate, which is interesting for high loading rate when the load measurement can be disturbed by the load cell resonance frequency.

Global bearing response of the laminates

Experimental results for the three tests performed at 1 m.s -1 plotted in Figure 5 illustrate the repeatability of the tests even for dynamic loading rates. The applied load is plotted on the y-axis and the pin displacement is plotted on the x-axis. Based on this figure, the bearing response of the laminate plate can be divided in several phases:

1. at the begin of the test, some displacements can be observed without applied load on the specimen, this phase corresponds to the decrease of the existing clearance in the experimental set-up, 2. next the assembly exhibits a linear behaviour, characterised by a rigidity denoted K e in the sequel, followed by the begin of the non linear behaviour of the assembly, 3. this non linear behaviour leads to the bearing peak for a an applied load denoted F peak (between 15 kN and 16 kN in Figure 5) and a pin displacement denoted δ peak (between 0.4 mm and 0.5 mm in Figure 5) followed by a brutal decrease of the strength of the assembly, 4. finally a load plateau appears after a stabilisation phase (for a displacement larger than 1.5 mm in Figure 5), the load plateau amplitude is denoted F plat in the following and is evaluated as the mean value between 5 mm and 10 mm.

The evolution of the bearing response with respect to the loading rate is plotted in Figure 6. For sake of clarity, the displacement is plotted until 9 mm as no major evolution is observed after such a displacement, even if all the signals have been recorded until the end of the test. The Figure 6 clearly exhibits a rate dependency of the bearing response for this laminate. The rate dependency is also observable in Table 1 and Figure ??, in which the results obtained on the different tests for the previously defined parameters (K e , δ peak and F plat ) are summarised. Firstly, an increase of the peak bearing load is observed with the loading rate increase, from an average value of 12.7 kN for 8.3•10 -5 m.s -1 tests to 15.5 kN for 1 m.s -1 tests. The pin displacement at the bearing peak also increases from 0.4 mm to 0.5 mm. Secondly, the magnitude of the load plateau is decreasing with the loading rate increase, from 10.3 kN for 8.3•10 -5 m.s -1 tests to 6.3 kN for 1 m.s -1 tests (see Table 1).

Contrary to the aforementioned parameters, the assembly linear behaviour does not appear to be loading rate dependent, this is an expected results for carbon/epoxy quasi-isotropic plates as the stiffness of the 0 o plies is not rate dependent. On the other hand, the strain rate dependency of the compressive strength of the 0 o plies could be responsible of the increase of the bearing peak, as the bearing peak is due to the kinking of the fibres in the 0 o plies [START_REF] Camanho | Failure mechanisms in bolted cfrp[END_REF]. The Infrared Thermography (IT) analysis is performed until the bearing peak at which the 45 o surface ply disturbs the measurement by being unsticked of the rest of the laminate. As the computer based triggering of the JADE III infrared camera is not as efficient as the direct triggering of visible cameras, the infrared films have been readjusted using the bearing peak as stopping point. Indeed, this brutal phenomenon with intense dissipative activity is easily detectable on infrared images. Given the properties of the infrared camera, measurements have only been performed for loading speeds lower to 0.1 m.s -1 . For the tests performed, the infrared camera configurations are summarised in Table 2.

In the following, the analysis of the thermography images is based on the three zones in the pin area described in Figure 8. The difference between the two radius of each zone is equal to the half of the radius of the pin. Here, damage-related small variations of the temperature are expected (see [START_REF] Fruehmann | Thermoelastic stress and damage analysis using transient loading[END_REF]), consequently the change of temperature ∆T is preferred to the absolute temperature value to perform the analysis. The temperature change is mainly due to the applied loading through thermoelasticity and the apparition of dissipative phenomenon through damage. In order to detect damage appearance, the change of temperature ∆T is computed using a sliding window with a length corresponding to 1% of the total number of images for the tests:

∆T (i) = T (i) -T (i - 1 100 N tot ) (1)
here N tot is the number of images for the studied test. As the load increase is almost linear until the bearing peak (see Figure 4), the change of temperature due to thermoelasticity is almost constant for all the ∆T computed. Given the length of the sliding window, its magnitude is limited.

In Figure 9, the evolution of the mean temperature change with respect to the applied load is plotted in zone 1 for the A-01 test. As expected, thanks to the sliding window processing, the temperature change ∆T does not evolve until the first damage event appearance. For an applied load of about 6.7

kN, the first significant ∆T variation highlights the appearance of a matrix crack in the 45 o surface ply long before the bearing peak. Shortly after the appearance of the crack matrix, another ∆T variation can be observed, with higher magnitude and duration. The two cartographies taken for an applied load of about 8.7 kN and 9.6 kN and illustrating this temperature change peak, illustrate a more diffuse dissipative phenomenon. This explains the higher magnitude observed since the mean value is more impacted by a diffuse phenomenon. It is assumed here that this dissipative phenomenon can be related to delaminations occurring near the loaded boundary of the hole. As a matter of facts, such local damages have already been reported in the litterature for quasi-static tests (see for instance Camanho et al. [START_REF] Camanho | Failure mechanisms in bolted cfrp[END_REF],

using Xradiography and sectioning).

In Figure 10, a comparison of the temperature change mean values in the different zones previously defined is performed on the same low speed tests.

These curves highlights the propagation of the dissipative phenomenon to distant regions well before the bearing peak. The magnitude of the peak decreases as the zone is farther from the pin. This can be explain in different ways. Firstly, as the mean value is considered, the magnitude of the peak is influenced by the distribution of the dissipative phenomenon on the con-sidered area. Secondly, the intensity of the dissipative phenomenon could be lower in an area far from the pin.

Digital Image Correlation analysis for a quasi-static test

In this study, Digital image correlation is used to analyse strain fields (ε xx , ε yy and ε xy ) ahead of the pin. In Figure 11, the longitudinal ε xx (with the x-axis corresponding to the loading direction), transversal ε yy and shear ε xy strain fields are plotted just before the bearing peak. Some discontinuities can be observed on these strain fields particularly in the area near to the pin for ε yy and ε xy . The shape of these discontinuities is similar to the matrix cracking in the 45 o ply, previously observed with IT measurement.

An efficient way to detect crack with DIC measurement is to analysed the residual field [START_REF] Roux | Digital image correlation and fracture: an advanced technique for estimating stress intensity factors of 2d and 3d cracks[END_REF]. As Aramis software is used in this study, the authors do not have access to the residual field. In this case, a visual analysis of the strain field can be performed in order to evaluate the appearance time of the crack. However, this analysis is user dependent and will lead to an important uncertainty. Consequently, a similar approach to the one applied on the thermal field in Section 3.3.1 is proposed in the following to analyse the strain fields.

The strain mean values are computed in the same three areas previously defined in Figure 8 for IT analysis. The variation of the mean values of the longitudinal strain ε xx in the three areas is plotted in Figure 12 with respect to the applied load. The figure clearly exhibits that the deformation is mostly localised in the first area. It can also be noticed that the evolution of ε xx in the first zone is initially linear and becomes more and more non-linear with the load increase. The transition between the linear and non-linear evolution occurs to appear between 6 kN and 8 kN. This order of magnitude can be linked to the increase of temperature previously observed in IT analysis.

Finally, the field of strain variation ∆ε during the test with the same sliding window as for the IT analysis can be studied. In Figure 13, the variation of the shear strain field ∆ε xy is plotted before and after a crack appearance in the surface ply. These figures clearly reveal the crack and only the crack. For this loading speed, the matrix crack is observed for a load applied on the specimen of 77.1% of the bearing peak load (see Table 3).

This value is higher than the one obtained for IT analysis (54.5% of the bearing peak load, see Table 2). This can be explained first by the fact that the measurements are performed on opposite faces which can lead to a delay. The delay can also by explained by the fact that due to the use of ZOI to compute displacement and strain field, DIC cannot capture the crack appearance close to the pin. The crack is captured by DIC when its propagation is sufficient, contrary to the infrared thermography which can measure dissipative phenomenon close to the pin.

Scattering of the test results

Even if the scattering is relatively low for the global measurement results, higher scattering has been obtained on local measurement results. For example, the comparison of the results obtained with IR analysis for two tests performed at 8.3•10 -3 m.s -1 is plotted in Figure 14 for zone 1. For the second test, the matrix crack appears for lower load level. Nevertheless, the observed scenario remains the same: the crack appearance in the surface 45 o ply is followed by a more diffuse dissipative phenomenon ahead of the pin. Such scattering can also be observed for DIC analysis at the same loading rate as shown in Figure 15a. This scattering could be due to differences in the way to apply the load to the composite plate. For example, a small coaxiality default can lead to a difference in the applied load to the surface ply. Finally, at higher loading rate, a smaller scattering is observed as it can be noticed in Figure 15b. The scattering decrease with the increase of the loading rate has also been observed for the results obtained with IT analysis.

Loading rate dependency

Concerning the influence of the loading rate, the evolution of the temperature change mean value is plotted in Figure 16 for a loading speed of 1•10 -3 m.s -1 and in Figure 17 for the tests performed at 8.3•10 -3 m.s -1 . Given the properties of the infrared camera, measurements have only been performed

for loading speeds lower to 0.1 m.s -1 . As previously, the appearance of a matrix crack in the 45 o surface ply is observed. The scenario is similar to one previously described, the cracking is followed by a diffuse dissipative activity ahead of the pin. As observed in Figure 17, for the higher loading rate multiple cracks are observed contrary to lower loading rate tests. The first conclusion regarding the influence of the loading speed is that the load level leading to matrix cracking in the surface ply is rate dependant (see Table 2).

Rate dependency is also observed in DIC results. In Figure 18, the evolution of the mean longitudinal strain field in zone 1 is plotted for different loading rate. The sliding window processing can only be efficiently applied on the 4 lower speeds (not enough points to obtain relevant results for the 1 m/s test). A delay in the appearance of the transition between linear and non linear evolution is observed with the increase of the loading rate. This results is in good agreement with the delay in the appearance of the matrix crack in the surface ply.

Finally, the analysis can also be performed regarding different areas ahead of the pin. For the 1•10 -3 m.s -1 , this analysis is performed in Figure 16. It appears that the area in which dissipative phenomenon happens is smaller than for the low speed tests. A limited dissipative activity is noticed in the second area of analysis and almost no activity in the third area. Again, these results are confirmed by tests performed at 8.3•10 -3 m.s -1 . The evolutions of the mean value of ∆T in the second and third area are respectively limited and unnoticeable. This leads to the second conclusion of this IT analysis which is that the area of dissipative activity ahead of the pin is less extended with the increase of the speed. The higher the speed, the shorter the damage area is.

Conclusions

In this paper, the loading rate dependency of the bearing response of quasi-isotropic T700GC/M21 composite laminates has been studied. For that purpose a specific double-lap experimental set-up has been designed in order to perform global and local measurements of the bearing response of the composite plate. Tests have been performed at different loading rates from 5 mm/min up to 1 m/s with a servo-hydraulic jack.

Regarding the global response of the laminate, the experimental results

show a significant influence of the loading rate. The peak bearing load increases by about 20% with respect to the quasi-static response when the loading rate increases from 8.3•10 -5 m/s to 1 m/s. On the same range of loading rate, the load plateau falls by about 40%. Conversely, the stiffness of the joints in its phase of linear behaviour does not seem to be influenced by the loading rate.

Concerning local measurements, DIC and IT exhibit local non linear behaviour, with for example matrix cracking in the surface ply, which is not observed on global measurement. The local non linearity observed with these two techniques is loading rate dependent. Firstly, a delay in the appearance of the non linear phenomenon is observed with the loading rate increase.

Secondly, the area of dissipative activity ahead of the pin is less extended with the increase of the loading rate. The higher the loading rate, the shorter the damage area is.

13 Map of the shear strain field variation ∆ε xy plotted just before the matrix cracking in the surface ply on the left hand side

  environment of the subsequent experimental investigation. The interval of temperature used for the calibration phase ranges from 18 o C to 40 o C by steps of 0.5 o C. For each temperature value 100 images of the SR800 blackbody (CI Systems) have been taken and the mean of these images has been used to generate the calibration law (6 order polynomial law) of each pixel using MATLAB. In this study, after calibration, the difference between the mean temperature measure by the camera and the real blackbody temperature was lower than 0.015 o C between 18 o C and 30 o C and lower than 0.02 o C between 30 o C and 40 o C.

3. 3 .

 3 Local bearing response of the laminates 3.3.1. Infrared Thermography analysis for a quasi-static test

F

  = 11750N and just after on the right hand side F = 11778N for a test speed of 8.3•10 -5 m.s -1 . . . . . . . . . . . . . . . . 36 14 Comparison of the evolution of the temperature change mean value for two tests performed at 8.3•10 -5 m.s -1 in the first zone close to the pin . . . . . . . . . . . . . . . . . . . . . . . 37 15 Evolution of the mean of the strain ε xx in zone1 for all tests performed for the loading rate of 8.3•10 -5 m.s -1 (left hand side) and for the loading rate of 0.1 m.s -1 (right hand side) . 38 16 Evolution of ∆T mean value for a test speed of 1•10 -3 m.s -1 in different zone close to the pin . . . . . . . . . . . . . . . . . 38 17 Evolution of ∆T mean value in the zone 1 (see Figure 8) for the three tests performed at 8.3•10 -3 m.s -1 . . . . . . . . . . . 39 18 Evolution of the mean of the strain ε xx for all speed in zone1 . 40 a-8.

Figure 1 Figure 2 :Figure 3 :

 123 Figure 1: Specimen dimensions

Figure 5 :

 5 Figure 5: Evolution of the different measured values for the three tests performed at 1 m.s -1

Figure 6 : 1 0

 61 Figure 6: Evolution of the bearing response with respect to the loading rates

Figure 7 :

 7 Figure 7: Test results evolution of the global bearing parameters with respect to the loading rate

Figure 8 :Figure 9 :Figure 10 :

 8910 Figure 8: Description of the three processing zones for the IT and DIC analyse

Figure 11 :

 11 Figure 11: Strain fields observed just before the bearing peak F = 12950N for a test speed of 8.3•10 -5 m.s -1

Figure 12 :

 12 Figure 12: Evolution of the mean of the strain ε xx for a test speed of 8.3•10 -5 m.s -1 in different zone close to the pin

Figure 13 : 5 m

 135 Figure 13: Map of the shear strain field variation ∆ε xy plotted just before the matrix cracking in the surface ply on the left hand side F = 11750N and just after on the right hand side F = 11778N for a test speed of 8.3•10 -5 m.s -1

Figure 15 :

 15 Figure 15: Evolution of the mean of the strain ε xx in zone1 for all tests performed for the loading rate of 8.3•10 -5 m.s -1 (left hand side) and for the loading rate of 0.1 m.s -1 (right hand side)

Figure 16 :Figure 17 :

 1617 Figure 16: Evolution of ∆T mean value for a test speed of 1•10 -3 m.s -1 in different zone close to the pin

Figure 18 :

 18 Figure 18: Evolution of the mean of the strain ε xx for all speed in zone1

Table 1 :

 1 3•10 -5 m.s -1 Test results evolution of the global bearing parameters with respect

	Test	Ke	F peak	δ peak	F plat
	Number	(kN.mm -1 )	(kN)	(mm)	(kN)
	A-01	39.6	12.4	0.41	10.3
	A-02	37.2	12.8	0.46	9.8
	A-03	38.0	12.9	0.44	10.2
	mean	38.3	12.7	0.44	10.1
	standard deviation	1.2	0.3	0.03	0.3
		b-1•10 -3 m.s -1
	A-04	40.9	13.4	0.4	9.9
	A-05	39.9	13.7	0.4	9.9
	A-06	42.2	13.8	0.4	10.2
	mean	41.0	13.6	0.4	10.0
	standard deviation	1.2	0.2	0.0	0.2
		c-8.3•10 -3 m.s -1
	A-07	39.1	14.6	0.42	8
	A-08	40.9	14.4	0.41	8.2
	A-09	37.8	14.3	0.40	7.7
	mean	39.3	14.4	0.41	8.0
	standard deviation	1.6	0.2	0.01	0.3
		d-0.1 m.s -1
	A-10	39.9	15.7	0.61	6.6
	A-11	43.2	15.6	0.63	7.4
	A-12	44.0	15.0	0.54	7.3
	mean	42.4	15.4	0.59	7.1
	standard deviation	2.2	0.4	0.05	0.5
		e-1 m.s -1
	A-13	38.4	15.4	0.51	6.3
	A-14	39.2	15.9	0.47	6.2
	A-15	36.6	15.2	0.47	6.3
	mean	38.1	15.5	0.48	6.3
	standard deviation	1.3	0.4	0.03	0.1
	to the loading rate

Acknowledgements

The authors gratefully acknowledge funding of VulcompII project by the French National Research Agency (ANR) and of the PRF TSD project by ONERA.

List of Figures


P Kerleau

Y Swesi

V Meille

Dr Isabelle Pitault
email: isabelle.pitault@lgpc.cpe.fr

F Heurtaux

Total catalytic oxidation of a side-product for an autothermal restoring hydrogen process

Keywords: Combustion, Dehydrogenation, Heat exchanger reactor, Coupling, Coated reactor

The catalytic total oxidation of toluene is used to bring the required calories to restore hydrogen from methylcyclohexane. Two coupled reactions, both catalyzed by Pt/alumina, are thus considered in one autothermal heat-exchanger-reactor (HER): endothermic dehydrogenation of methylcyclohexane and catalytic combustion of less than 10% of produced toluene (6% for the chemical reaction, about 4% for heat losses). Both reactions are carried out at 350°C. Methylcyclohexane and toluene light-off curves have been performed in one elementary module. Methylcyclohexane combustion has been eliminated due to its safety drawbacks, although the MCH combustion provides best stability to autothermal HER.

Introduction

To feed an on-board fuel cell, hydrogen can either be produced in situ by a chemical reaction or chemically or physically stored. The here-studied process consists of releasing hydrogen from naphthenic molecules, hydrogen storing compounds [1][2][3]. The aromatic compounds, the by-products of the dehydrogenation reaction, can be recycled ex-situ in a hydrogenation plant to restore the hydrogen-rich reactants [4]. The molecules are used as hydrogen carriers. In the global processing system, the hydrogen and aromatics mixture is first separated by condensation at room temperature. Aromatics are afterwards stored in a tank and hydrogen (containing 2%vol. Aromatics for MCH/Toluene couple) is then purified using adsorption on polymers. After purification, aromatics content in hydrogen being lower than 1ppm, hydrogen can safety feed a FC [5].

The chosen hydrogen-rich reactant is methylcyclohexane (MCH) because of its high hydrogen storage capacity (6 wt-%) and the lower toxicity of the produced toluene (TOL) in comparison with benzene. As dehydrogenation of MCH is highly endothermic and occurs above 300°C, heat can be provided by burning toluene with both reactions performed in a catalytic heat exchanger.

The heat exchanger reactors (HER) consist of coupling both reactions (endothermic and exothermic) and recovering the heat produced by the exothermic reaction to heat the endothermic reaction. In this study, to reach highest efficiency in the HER, the oxidation reaction has been performed on coated stainless steel plate walls that are in contact with the dehydrogenation reaction. This configuration corresponds to recuperative mode in Kolios's terminology [6; 7]. The main advantage of using stainless steel plates is the high thermal conductivity, which reduces the thermal gradient within the walls and the required transfer area. Moreover, the use of structured plates increases the turbulence, which also reduces choking compared with the use of tubes. Nevertheless, high thermal conductivity should also be a drawback, it increases heat losses and the energy necessary to ignite the combustion. One of the first catalytic HERs with coated walls was described by Hunter and McGuire [8] in a patent. They presented the use of coated catalyst to improve heat transfer and decrease combustible consumption in a domestic boiler. The flame absence inducing a low production of gas pollutants is presented as an advantage. Their main difficulty was the distribution of heat flows.

The goal of this study is to show the feasibility of coupling a fast exothermic reaction such as combustion with a slow endothermic reaction such as dehydrogenation, focusing on combustion reactors. Firstly, the choice of combustible based on light-off curves and the use of structured plate for combustion are discussed. Secondly, the reactions of dehydrogenation and combustion are independently studied to quantify heat losses and give critical estimates of flows and concentrations suitable to be coupled in an HER. Finally, an improvement in the design of the HER is shown.

Experimental

Toluene (TOL) is used as received from Acros. It is pumped by a HPLC pump (Gilson 305) and vaporized via a T-junction through a capillary tube perpendicular to the air flow. The TOL-air mixture is preheated up to 180°C and introduced into the combustion reactor. The range of operating conditions is 1-3Nlair.min -1 , 1-2.2% TOL, 350-450°C. The combustion effluents are cooled to room temperature and water is condensed in a Peltier condenser (M&C Cooler ECP 100) before analysis in ABB detectors (O 2 : ABB Magnos 106 (Magnetomechanical detector; range 0-100 vol%); CO 2 : ABB Uras 14 (IR absorption; range 0-100 vol %); CO: two ABB Uras 14 (IR absorption; ranges 0-1000 ppm and 0-100 vol%). Each detector is equipped with internal calibration allowing a precision greater than 0.5%.

At the inlet of the dehydrogenation reactor, methylcyclohexane (>99%, Acros) is pumped by a HPLC ISMATEC pump (Reglo CPF/FMI 005), vaporised and preheated (250°C). The evaporator is close to the reactor inlet to avoid heat loss. MCH enters the reactor with a small amount of hydrogen to avoid rapid deactivation of the catalyst. For activation and regeneration purposes, air and N 2 can also be fed. The range of operating conditions is 2-11 mlMCH.min -1 , 350-400°C, 1.5 bars. The dehydrogenation products are condensed by two condensers in series. A smaller third condenser is used to sample products before GC analysis. This configuration is used to protect the Hydrogen detector. MCH conversions are measured using a GC equipped with PONA column.

The same nature of catalyst is used for combustion and dehydrogenation reactions. It was prepared by wet impregnation of gamma-alumina (SASOL, SBa-200) with a toluenic solution of platinum acetylacetonate allowing Pt/Al 2 O 3 catalyst to be obtained. The powder obtained is used in the dehydrogenation reactors. Concerning the combustion reactors, the coated walls were obtained by dip-coating in a suspension of the catalyst in diluted nitric acid, followed by drying and calcination steps.

One single combustion reactor (A) and two different heat-exchanger reactors (B and C) (HER) were designed in the lab and tested in the experimental set-up presented in Figure 1.

All reactors and modules are made of stainless steel. The reactors can be heated up to 600°C by two copper heating plates positioned above and below the stacks.

Single combustion module (A)

The single combustion module (5 cm long * 5cm wide * 0.5 cm thick) is coated with 0.12 g of 2.1wt-% Pt/alumina (Figure 2). A Ni-Cr foam is used at the gas inlet to increase the pressure drop and evenly distribute the reactants.

First generation heat-exchanger reactor (B)

It consists of a long and thin dehydrogenation reactor (15 cm long * 5 cm wide * 0.4 cm thick) coupled to 6 combustion reactors (5 cm long * 5 cm wide * 0.5 cm thick) (Figure 3).

The dehydrogenation reactor is filled with 22 g of 0.88wt-% Pt/alumina, whereas each combustion reactor contains 0.12 g of 2.1wt-% Pt/alumina, coated on the heat-exchanging walls. The temperatures are measured with thermocouples located in a space between the combustion modules and the dehydrogenation reactor, one of them being inserted in the powder bed as Figure 4.

Second generation heat-exchanger reactor (C)

The second generation is more compact than the first [9]. Its dimensions are 5.4 cm long, 5.4 cm wide and 2.5 cm high. The HER consists of 6 plates, 4 of them are 4.5 mm thick and 2 of them (at the top and the bottom of HER) are 3 mm thick (Figure 5). The four central plates in the HER centre each form half the dehydrogenation reactor on one side and half the combustion reactor on the other side. The 2 end plates comprise only the combustion reactor halves. Each dehydrogenation section consists of one channel 45 cm long (4mm wide, 3mm high), containing a mean of 4.6 g catalyst powder (200 µm diameter particles of 0.5 wt-%Pt/Al 2 O 3 ) and coiled in a 5.4 cm*5.4 cm plate. Each combustion section consists of a 5cm*5cm*3mm chamber (0.7 g of 2.1 wt-%Pt/Al 2 O 3 ) containing pillars and wash-coated on its two faces. The HER is equipped with ten thermocouples. Four of them measure inlet and outlet gas temperatures. These temperatures are necessary to calculate the heat balance and loss. The others (T_1 to T_6) are inserted in the walls between the dehydrogenation and combustion chambers (Figure 6). They give indications of temperature profiles along and across the reactors. The arrows show the flow direction for MCH in the dehydrogenation reactors. Combustion chambers can be fed in co-or counter-current mode with respect to MCH.

For HER B and C, the stacks are first electrically heated up to 360-375°C, under N 2 or air flow through the dehydrogenation and combustion reactors respectively. Then, MCH is fed to the dehydrogenation reactors. Finally, TOL is injected into the combustion reactors.

Results

First generation reactor -1 single combustion module (configuration A)

In a first part, MCH and TOL combustions were studied in one of the six modules of HER B.

The programmed temperature curves (Figure 7) show that the ignition and extinction temperatures of TOL are greater than for MCH. The MCH extinction curve could not be obtained with the operating conditions in our set-up, heat losses being lower than heat produced by combustion. Although the TOL combustion was unreliable below about 300°C, toluene was chosen due to the risk of flashback for MCH.

The temperature reached during toluene combustion evidently depends on the air flow rate and the toluene content (Figure 8). Moreover, a temperature gradient between the inlet and the outlet of each combustion module is observed. The localization of the hot spot also varies from the entrance of the module at low toluene fraction to the exit of the module at higher toluene fraction, with a non-linear propagation front, despite the use of the metallic foam as a distributor. Despite the high thermal conductivity of stainless steel and the thickness of the walls, the thermal gradient remains relatively high, about 10 to 15°C. From these results, the six modules were laid out so as to lightly stagger the 6 hot spots on the both sides of HER B.

In this case, the axial thermal flow provided by the modules to the dehydrogenation reactor is more homogeneous.

First generation reactor -3 top combustion modules (configuration B)

The global heat-exchanger reactor was validated step by step. Firstly, three combustion modules running together were studied. The temperatures observed at the wall of the different modules are reported for various flow rates and toluene fraction. The difference between the top and the bottom modules is due to the fact that only the top modules are fed with air and toluene. What is more surprising is that the temperature decreases from Ta6 to Ta2, showing that the distribution of air and toluene is not the same in each module (Figure 9). To overcome this drawback, the pressure drop in each module was measured and a restriction was added at the entrance to each module to give all modules the same pressure drop. These measures were, however, not sufficient.

First generation reactor -dehydrogenation (configuration B)

The dehydrogenation reaction was also studied alone. Temperature profiles in the dehydrogenation reactor (central part of the heat-exchanger reactor) (Figure 10) show a large temperature gradient between walls and bed centre. Moreover, the greater the MCH flow is, the greater the temperature gradient is. This gradient is entirely due to low thermal conductivity of the catalytic bed. On the other hand, the temperature gradients for Ta3 and Tb2 are quite similar and show the catalytic bed is symmetrical, the thermal flows induced by the dehydrogenation reaction should be equal above and below the dehydrogenation reactor.

First generation reactor -coupling (configuration B)

Despite the disappointments encountered during the study of each reaction, the coupling of both reactions was attempted (Figure 11). The dehydrogenation is started first, with electrical heating. After 70 minutes, the combustion (in the 6 modules) is started. The temperature rise in the dehydrogenation section is now due only to the calories produced by the combustion modules. The temperature controller reduces the electrical heating to zero. The coupling works, but a very large temperature gradient (50°C) is observed between the centre of the catalyst bed and the wall at the exit of the reactor. Moreover, the observed temperature gradients in the wall between the reactor entrance and exit using the three modules increase as the dehydrogenation progresses. In this situation, energy produced by combustion at the dehydrogenation reactor entrance provided insufficient heat to the reactor. Whilst, energy produced by combustion at the dehydrogenation reactor exit provided an excess of heat to the reactor.

After the all experiments, a single module was disassembled at the end of the combustion experiments. This allowed the observation that the quality of the coating remained good, but some corrosion of the stainless steel (Figure 12) and decomposition of the Ni-Cr foam was also observed. The Ni-Cr was stressed during the temperature change and turned into powder.

The presence of black marks suggests homogeneous combustion in the foam.

The second generation reactor was designed with these observations in consideration.

Second generation reactor -coupling (configuration C)

The internal structures and arrangement of the plates have been designed in order to promote heat transfer between dehydrogenation and combustion catalysts while maintaining ideal plug flow in each reactor. The foam used to homogenize the gas flow was removed and the structures of the combustion plates were designed to break the flow and increase the turbulence. To avoid an axial thermal gradient due to poor gas distribution in the tubes between several modules, a stack of plates was designed with the same entrance box and the same exit box for all combustion plates. The same configuration was chosen for the dehydrogenation chamber entrance and exit.

The main observed behaviours in the HER C (Figure 13) are the following. Temperatures are nearly homogeneous in the HER under electrical heating with no reaction. With electrical heating and MCH feed, the significant temperature undershoot close to the inlet illustrates the endothermicity of the reaction. The temperature measured at the same distance from the inlet, but at different depths in the stack (T_1&T_2; T_3&T_4) are close to each other but significantly different. The temperatures measured with the thermocouples furthest from the heating plates are evidently the coldest. MCH conversions obtained during the periods (2) of Figure 13 are about 85%.

In the second heat-exchanger reactor configuration, the reactors being perfectly symmetrical, two different modes of operation were studied: counter-current and co-current. Autothermal conditions are successfully achieved both in co-current and counter-current flow modes. In both cases, the use of central combustion module removes the temperature difference between T_1&T_2 or T_3&T_4. Figure 13 clearly shows the advantage of co-current operation. Lower temperature excursions along the reactor and higher MCH conversions are achieved in cocurrent flow mode than in counter-current flow mode. In co-current flow mode, the coupling occurs at a constant temperature of 330°C. But, in all cases, T_1 (resp. T_3) is much closer to T_2 (resp. T_4) because heat supply (combustion) is uniformly spread over the stack depth.

The MCH conversions achieved in co-current flow are 78% against 59% in counter-current flow. In all cases, TOL combustion is total.

Comparison between HER B and HER C:

The enthalpies of MCH dehydrogenation and Toluene combustion being 204kJ/mol and 3771kJ/mol respectively, only 6% of the produced toluene is necessary to heat the reactor with no heat loss. For HER B and HER C, 10.0% and 9.7% of the produced toluene was burned with a MCH flow of about 3.9 g/min and the same hydrogen production. Heat losses are therefore of the same order of magnitude for both laboratory reactors. Moreover, the higher the MCH flows are, the lower the toluene proportions to compensate for heat loss are (for example, for both the HER, for 5g/min MCH flow, 7.4% of the produced toluene is burned and for 2g/min MCH flow, 13.9% of the produced toluene is burned). Both reactors have the same efficiency in terms of hydrogen production. However, the HER B is less compact than the HER C. The masses of dehydrogenation catalysts in the reactors are very different: 22g of 0.88wt-%Pt (i.e. 193mg of Pt) and 9g of 0.5wt-%Pt (i.e. 45mg of Pt) in HER B and HER C respectively. The masses of combustion catalysts are 0.7g and 2.4g of the same 2.1wt-%Pt in HER B and HER C respectively. In all combustion modules, the reaction is hardly limited by mass and heat transfer, the mass of Pt only affects the reaction rate during the ignition step. In conclusion, if the compactness is an effectiveness criterion, the HER C is better than HER B.

Conclusions -perspectives

The coupling of one endothermic reaction and one exothermic one in two autothermal heatexchanging reactors was carried out with success.

Present work has concerned the study of different designs of the structured heat-exchanger reactor. It highlighted the problems of corrosion linked to the coating procedure on stainless steel walls. This drawback is what stopped the future process development. Methods to protect walls have to be developed. Nevertheless, both the reactors have dramatically enhanced the efficiency of the dehydrogenation reaction, as the same conversion can not be reached at the same WHSV in classical reactors.

Figure 1 .

 1 Figure 1. Experimental set-up to carry out the coupling of toluene combustion and methylcyclohexane dehydrogenation.

Figure 2 .

 2 Figure 2. Coated catalyst in the module A structure.

Figure 3 .

 3 Figure 3. Autothermal reactor (HER B) coupling MCH dehydrogenation (blue central section) and TOL combustion (red sections).

Figure 4 .

 4 Figure 4. Diagram of thermocouple locations in the autothermal reactor (HER B).

Figure 5 .

 5 Figure 5. Second generation heat-exchanger reactor (HER C). From left to right: stack of plates, combustion zone, dehydrogenation zone.

Figure 6 .

 6 Figure 6. Thermocouple positions in the HER C. Left: top view. Right: side view. Pink rectangles: Combustion chambers. Yellow rectangles: Dehydrogenation Beds. Gray lines: Boundaries of the 6 plates.

Figure 7 .

 7 Figure 7. Extinction and ignition curves of combustion of toluene and methylcyclohexane in

Figure 8 .Figure 9 .

 89 Figure 8. Temperature profiles in one combustion module at different toluene fractions (from left to right: 0.61%, 1.21%, 1.51%, 1.80%, 2.09%, 2.11%). D air =0.3 l STP .min -1 , regulation temperature = 350°C

Figure 10 .

 10 Figure 10. Temperature profiles in the reactor for MCH dehydrogenation, heat is provided by electric power. Between 0 and 20 min, D MCH =5 ml.min -1 ; between 20 and 56 min, D MCH =7.8 ml.min -1 ; After 56 min, D MCH =10.4 ml.min -1 .

Figure 11 .

 11 Figure 11. Dehydrogenation and toluene combustion coupling. Preheating with electric power and nitrogen flows. At 30 min, start of MCH feeding (D MCH =3.9 ml.min -1 ). At 68 min, start of toluene and air feeding (D air =2.0 l STP .min -1 ; y TOL =1.97%).

Figure 12 .

 12 Figure 12. Single combustion module after run.

Figure 13 .

 13 Figure 13. Dehydrogenation and toluene combustion coupling. (1) Preheating with electric power and nitrogen flows at 390°C. (2) MCH feeding (D MCH =2.3 ml.min -1 ). (3) Toluene and air feeding (D air =3 l STP .min -1 ; y TOL =0.25%) in counter-current mode. (4) Toluene and air feeding (D air =3 l STP .min -1 ; y TOL =0.25%) in co-current mode.

Acknowledgements

Renault Direction des Technologies Automobile Avancées is warmly acknowledged for the financial support of the present study.


Gratien Bonvin
email: gratien.bonvin@mines-paristech.fr

Sophie Demassey
email: sophie.demassey@mines-paristech.fr

Welington De Oliveira
email: welington.oliveira@mines-paristech.fr

Robust design of pumping stations in water distribution networks

Keywords: Pressurized Water Network Design, Stabilized Benders' Decomposition, Mixed Integer Nonlinear Programming

Restricted to gravity-fed networks, most water network design models minimize investment costs under a static peak water demand scenario. In networks equipped with pumping stations, design models should also account for operation costs incurred by the pump energy consumption that depends on dynamic demand and tariff. Evaluating the lifetime operation costs amounts to solve a large-scale non-convex combinatorial optimization problem for each considered design. In this paper, we address the pressurized water network design problem with a joint optimization of the pump investment and operation costs through a stabilized Benders' decomposition. To reduce the complexity of the operational subproblem, we decompose the scheduling horizon in representative days, and relax the discrete and non-convex components of the hydraulic model. We also evaluate the design robustness on stress-day scenarios and derive feasibility cuts using a dominance argument. Experiments on a typical rural branched water distribution network with one year of historical data show the accuracy of our approximations and the significant savings expected from the optimal pump resizing.

Introduction

While the lifetime of pipes and water tanks usually reaches 100 years, the mean lifetime of a pump is closer to 20 years. Operators of water networks must then periodically proceed to the rehabilitation of pumping stations with the characteristics of the other network assets already fixed. The problem is complex because, besides the strategical level -which pump combination to install? -it requires to investigate the operational level -how to operate the installed pumps? -to evaluate and minimize the lifetime costs over the set of pump combinations. Moreover, evaluating the minimum operation costs brings into play, together, dynamic water demand and energy tariff profiles, discrete pump scheduling decisions, non-convex hydraulic laws, and uncertain long-term forecasts. Actually, solving the operational subproblem, known as pump scheduling problem, deterministically on a daily horizon is already considered challenging (see e.g. [START_REF] D'ambrosio | Mathematical programming techniques in water network optimization[END_REF]).

Optimization methods dedicated to the design of water distribution networks have been proposed for more than four decades. A common approach is to combine an evolutionary algorithm with a hydraulic simulator (see e.g. [START_REF] Marques | Multi-objective optimization of water distribution systems based on a real options approach[END_REF][START_REF] Murphy | Optimum design and operation of pumped water distribution systems[END_REF] and references therein). The approach deals accurately with the short-term dynamic of the network operation, but it is inherently a heuristic and provides no performance guarantee. Mathematical programming approaches handle the hydraulic explicitly as non-linear constraints to address the pipe layout design problem [START_REF] Alperovits | Design of optimal water distribution systems[END_REF][START_REF] Bragalli | On the optimal design of water distribution networks: a practical MINLP approach[END_REF][START_REF] Raghunathan | Global optimization of of nonlinear network design[END_REF]. They provide guarantee, but neglect the dynamic by evaluating the feasibility of operating the network on a static worst-case water demand scenario.

In this paper, we address the pumping station design problem precisely, with an optimization approach combining the strengths of both previous strategies. While the overall approach is generic, the operational level is presented in the context of the FRD network, a branched network that is described in details in [START_REF] Bonvin | A convex mathematical program for pump scheduling in a class of branched water networks[END_REF]. The approach is built on a decomposable formulation of the problem including the operation scheduling subproblem as a large-scale non-convex mixedinteger non-linear program with uncertain data. We apply a stabilized Benders' decomposition [START_REF] Van Ackooij | Inexact stabilized benders' decomposition approaches with application to chance-constrained problems with finite support[END_REF] to solve the problem and propose different approximations for the subproblem to reduce its complexity while maintaining some performance guarantees. First, we decompose the multi-year scheduling horizon in a restricted set of seasonal representative days. This common assumption in long-term planning models fits well with the daily periodicity of drinking water distribution, and allows to separate the scheduling subproblem in independent daily subproblems (once a pump configuration is given). As optimizing iteratively each subproblem would remain too time consuming, we propose to relax the integrality constraints and to convexify the hydraulic constraints. This results in convex continuous non-linear programs providing under-estimates of the operation costs and dual information to derive Benders cuts. Second, we handle the long-term uncertainties by enforcing the robustness of the solution on hypothetical stress days characterized, e.g., by a high water demand and the outage of one pump. By disregarding optimality but forcing feasibility on stress days we can, for the class of networks in consideration, aggregate identical pumps and, therefore, reduce the size and complexity of the subproblem. Furthermore, we exhibit a dominance relation between pump combinations that allows to generate multiple feasibility cuts from one infeasible solution. Our experiments on the FRD network show the accuracy of our approximations: the impact of the horizon decomposition is negligible while the continuous convex relaxation induces a deviation of the optimum lower than 4.3%. Finally, the annual savings (in terms of investment and operation costs) expected from the rehabilitation are estimated to 32%.

The paper is organized as follows: Section 2 presents a formulation of the problem. Section 3 describes the stabilized Benders' decomposition and Section 4 the management of infeasibility under the dominance argument. Computational results are provided in Section 5 and some conclusions in Section 6.

Optimal design of a pumping station

This section defines our design model, by first describing how the network operates in our application case. This simpler operation model detailed in [START_REF] Bonvin | A convex mathematical program for pump scheduling in a class of branched water networks[END_REF] is a specialization of standard models to the network characteristics [START_REF] Burgschweiger | Optimization models for operative planning in drinking water networks[END_REF].

Operation of a branched network

The water distribution network FRD [START_REF] Bonvin | A convex mathematical program for pump scheduling in a class of branched water networks[END_REF] can be represented as a directed graph G = (J, L). The water flows from a source r ∈ J to elevated tanks j ∈ J T ⊆ J through, successively: fixed speed pumps k ∈ K ⊆ L set in parallel at the pumping station s ∈ J, directed pipes l ∈ L P ⊆ L connected by junctions j ∈ J J ⊆ J, and pressure reducing valves l ∈ L V ⊆ L. The dynamic state of the system over a given time horizon T is driven by the water demand D ∈ R J T ×T + at the tanks and is governed by complex hydraulic laws of conservation of flow and pressure through the network. A standard model is defined as

P K T = {(x, q, h) ∈ {0, 1} K×T × R L×T + × R J×T + | (1) -(9)},
with x the on/off state of the pumps, q the flow through the arcs, and h the head (sum of pressure and elevation in m) at the nodes, and:

ij∈L qijt = ji∈L qjit, t ∈ T , j ∈ JJ (1)
ij∈L qijt = Sj(hjt -hjt-1) + Djt, t ∈ T , j ∈ JT (2)
hj0 = hjT = H 0 j , j ∈ JT

H min j ≤ hjt ≤ H max j , t ∈ T , j ∈ JT (4) hit -hjt ≥ 0, t ∈ T , ij ∈ LV (3)
α min ≤ hst -hrt ≤ α max , t ∈ T (5)
Q min κ k x kt ≤ q kt ≤ Q max κ k x kt , t ∈ T , k ∈ K (7) hit -hjt ≥ Φij(qijt), t ∈ T , ij ∈ LP (8) hst -hrt ≤ Ψκ k t(qkt) + M (1 -x kt ), t ∈ T , k ∈ K. (6)
In this model, the time horizon is discretized T = {1, . . . , T } with a resolution of typically 1 hr or 2 hrs in which the system is assumed to operate in steady state.

Flows q lt and demands D jt are given in volume (in m 3 ) for the duration of a time step t ∈ T . Constraints (1) and ( 2) enforce the conservation of flow at junctions and tanks. In (2), a tank j ∈ J T is assumed to be a vertical cylinder of area S j (in m 2 ). Constraints (3) fix the initial and final volumes and Constraints (4) the security limits with 0 ≤ H min j ≤ H 0 j ≤ H max j (in m 3 ). Constraints (5) model the pressure reduction at the valves. Constraints (6) enforce bounds 0 < α min ≤ α max on the head increase to limit leaks for instance. Let κ k denote the class of pump k ∈ K. Constraints (7) limit the pump operation range, given 0

< Q min κ ≤ Q max κ (in m 3
), and bind flow values and pump activation states.

Constraints (8) models the head loss due to friction in pipes. For each pipe ij ∈ L p , the head loss-flow coupling function Φ ij can be accurately approximated by a quadratic function 9) synchronize (given a large enough M value) the head increase of each active pumps. Function Ψ κt depends on the manufacture characteristics κ and on the ageing t of the pump. It can be accurately fitted from operating points as a quadratic function Ψ κt (q) = α κt -β κt q 2 , with α κt ≥ 0 and β κt ≥ 0. We highlight that the head-flow coupling constraints ( 8) and ( 9) are actually equalities in the original -thus non-convex -formulation of the pump scheduling problem, but it is shown in [START_REF] Bonvin | A convex mathematical program for pump scheduling in a class of branched water networks[END_REF] that the optimality gap of this relaxation is small.

Φ ij (q) = A ij q + B ij q 2 with A ij ≥ 0, B ij ≥ 0. Constraints (
Finally, the financial cost of operation plan (x, q, h) is mainly incurred by the purchase of the electricity consumed by pumping, namely:

C K T (x, q, h) = t∈T C t k∈K Γ κ k t (x kt , q kt ), (10)
with C t ≥ 0 the actualized electricity price on period t and Γ κt the power consumption function for each active pump of class κ and ageing t defined as a linear fit Γ κt (x kt , q kt ) = λ κt x kt + µ κt q kt .

Robust design and relaxed operation

In the considered water network design problem, only the set K of pumps must be sized in a way to satisfy the future water demand D and to minimize the global cost -sum of the actualized investment and operation costs -over a life span T of typically 20 years. The number of pumps in K is limited by the capacity N ∈ N * of the pumping station. Each pump is selected from a given set of candidate classes κ ∈ K and acquired new at time t = 0 at a fixed investment cost I κ ≥ 0. We assume that the maximal efficiency and the head increase at constant flow for all pumps decrease from 1% each year, according to the empirical ageing model of [START_REF] Bunn | Ageing pump efficiency: The hidden cost thief?[END_REF]. Under this hypothesis, we assume in the definition of Φ ij that α κt = (1 -vt 100 )α κ0 and β κt = (1 -vt 100 )β κ0 for any time t ∈ T with v t the duration from time 0 to t in years. Also, because the power consumption can alternatively be formulated as the product of the flow and the head increase divided by the efficiency, we assume that it does not change in time, i.e. λ κt = λ κ and µ κt = µ κ for t ∈ T in the definition of Γ κt .

Assuming that the water demand (D t ) t∈T and the electricity price (C t ) t∈T are given, the optimal design problem is to find a set K of at most N pumps of classes in K and an operation plan (x, q, h) ∈ P K T of minimum investment and operation costs. As the life span extends into years while the operation model requires a time resolution in minutes or hours, this model is not practical due to its complexity, its dimension and the stochasticity of the water demand. To address these issues, we formulate a robust variant of the problem where the time horizon is decomposed into a given set of regular (resp. peak) days d ∈ D R (resp. d ∈ D P ). Each day d ∈ D = D R ∪ D P is characterized by water demand and electricity tariff profiles over the daily time horizon T d which are representative for a number L d of days over the life span T . Thanks to Constraints (3) (which forces the tank's water level at the end of the day to be equal to that of the beginning of the day) the subproblem of determining the minimum cost operation plan (x, q, h) ∈ P K T associated to a pump combination K is decomposed into independent daily subproblems: min

(x,q,h)∈P K T C K T (x, q, h) = d∈D L d min (x d ,q d ,h d )∈P K T d C K T d (x, q, h). (11)
Because optimizing on P K T d remains challenging, we propose to relax some operational constraints which are not structural for the long term horizon. We consider a different relaxation R K T d whether d is a regular or peak day. The purpose of peak days is to enforce the robustness of the solutions by simulating stress operation conditions, including a high demand. We check the feasibility of the pump combinations under these extreme conditions but, if feasible, we neglect the associated operation costs as these conditions are considered as exceptional. More specifically, we consider the sum in [START_REF] Perry | Perrys chemical engineers handbook 7th ed[END_REF] only over the index set of regular days D R , and we use the peak days to check feasibility (of the pump configuration), as detailed in Section 3 below. In such stress days we also ignore the minimum flow constraints in [START_REF] D'ambrosio | Mathematical programming techniques in water network optimization[END_REF] by setting Q min κ = 0. Under these assumptions, all the pumps in a class can be assumed to operate equally and the size of model R K T d can be reduced accordingly, by aggregating the variables per pump class and by relaxing the integrality constraints on the aggregated state variables for each class κ ∈ K such that α κ0 ≥ α max .

To allow a quick computation of the operation costs of a pump combination over the regular days, as well as the use of dual information in the overall solution algorithm described in the next section, we define R K T d as the continuous relaxation of P K T d for any regular day d ∈ D R . I.e., the feasible sets

P K T d in (11) are replaced with R K T d .

Pump investment variables, constraints and costs

As the considered network has only one pumping station, a pump combination is uniquely determined by the number of pumps of each class it contains. An alternative representation, with more symmetries but that is better suited to Benders cut generation, is given by y κn ∈ {0, 1} for κ ∈ K and n ∈ {1, . . . , N } such that y κn = 1 if and only if the combination has at least n pumps of class κ. For a fixed ordering of the pumps, each combination corresponds to exactly one binary vector y such that

Let N κ (y) = N n=1 y κn be the number of pumps of class κ in the combination y. The investment cost of acquiring N κ (y) pumps of class κ is I κ N κ (y), with given I κ ≥ 0. Having defined all the elements related to the variables of investment and operation, we are now in position to state the considered formulation for the robust design optimization problem:

     min y,x,q,h κ∈K IκNκ(y) + d∈D R L d C K T d (x d , q d , h d ) s.t. y binary satisfying (12) and x d ≤ y ∀d ∈ DR ∪ DP (x d , q d , h d ) ∈ R K T d ∀d ∈ DR ∪ DP . ( 13
)
The constraint x d ≤ y means that a pump is operational at day d if the configuration y accounts its installation.

Benders' decomposition

Problem (13) has a decomposable structure: we can separate the investment variables y from the operational ones (x, q, h) to rewrite (13) as min y f (y) s.t. y binary satisfying [START_REF] Raghunathan | Global optimization of of nonlinear network design[END_REF] and c(y) = 0 .

In this formulation, the objective function f is

f (y) = κ∈K I κ N κ (y) + d∈D R L d      min x d ,q d ,h d C K T d (x d , q d , h d ) s.t. (x d , q d , h d ) ∈ R K T d ∀d ∈ D R x d ≤ y ∀d ∈ D R . (15)
Peak days are handled by the constraint function c:

c(y) = 0 if there exist (x d , q d , h d ) ∈ R K T d s.t. x d ≤ y ∀d ∈ D R ∪ D P ∞ otherwise.
In the process of solving (14), both f and c are approximated by cuts. Benders cuts. Note that function f is convex but nonsmooth. Therefore, optimality cuts for ( 14) is nothing but linearization of f , computed by making use of subgradients. Let ỹ be a given combination of pumps. Then a subgradient s of f at ỹ is a vector of same dimension as ỹ constructed from the dual variables associated to the constraints x d ≤ y, for all d ∈ D R , as described in [1, Lemma 1]. The optimality cut is thus given by f (ỹ) + s, • -ỹ and satisfies (due to convexity of f ) f (ỹ) + s, y -ỹ ≤ f (y) for all y. Feasibility cuts. Suppose that ỹ is infeasible for problem (14). In order to exclude such a point from the set of candidate solutions for (14) we consider the following linear constraint κ∈K y κNκ(ỹ)+1 ≥ 1, with convention y κN +1 = 0.

Stabilized cutting-plane algorithm

Let be an iteration counter and z (resp. z ) be the lower (resp. upper) bound available at iteration for the optimal value of (14). We initialize z 0 = 0 and z 0 = ∞ in the stabilized Benders' decomposition of [START_REF] Van Ackooij | Inexact stabilized benders' decomposition approaches with application to chance-constrained problems with finite support[END_REF] to generate a sequence {y } of be the trial points for (14). We split such a sequence according to feasibility of its elements: O = {ι ∈ {0, . . . , } : c(y ι ) = 0} and F = {0, . . . , }\O .

Note that F gathers (up to iteration ) the points that have been proved infeasible for problem (14). By setting = 0, starting with a vector y 0 and defining the incumbent point ŷ0 = y 0 , our variant of [START_REF] Van Ackooij | Inexact stabilized benders' decomposition approaches with application to chance-constrained problems with finite support[END_REF]Alg. 4] defines trial points by iteratively solving the following MILP:

y +1 ∈          arg min y 1 2 1 -ŷ , y s.t.
y binary satisfying ( 12)

f (y ι ) + s ι , y -y ι ≤ (z + z )/2 ∀ι ∈ O κ∈K y κNκ(y ι )+1 ≥ 1 ∀ι ∈ F . ( 16
)
If this master program is infeasible, then z +1 = (z + z )/2 is a valid lower bound for ( 14). Otherwise, we set z +1 = z and check feasibility of z +1 : if it is infeasible we set F +1 = F ∪ { + 1} (a new feasibility cut is added). On the other hand, if c(z +1 ) = 0 we then update

O +1 = O ∪ { + 1} (a new Benders cut is added), compute f (z +1
) and a subgradient s +1 , and set z +1 = min{z , f (y +1 )}. If z > f (y +1 ) then the algorithm updates the incumbent point: ŷ +1 = y , otherwise it remains unchanged. In any case, the algorithm updates = + 1 and repeats this process until (z -z )/z ≤ tol, for a given tolerance tol > 0. We refer to [START_REF] Van Ackooij | Inexact stabilized benders' decomposition approaches with application to chance-constrained problems with finite support[END_REF] for further details on the algorithm.

Infeasibility, dominance and fault tolerance

In this section, we show how to accelerate the convergence of the solution algorithm by using a dominance relation on the set of combinations to generate more than one feasibility cut at a time, including in a preprocessing step. To simplify the presentation, we fairly assume that a combination which is infeasible for a regular day is also infeasible for a peak day. We finally integrate to the definition of peak days a model of robustness to pump failures.

Dominance

By definition, the maximal flow a pump combination y can deliver for a given head increase α = h s -h r corresponds to activate all N κ (y) pumps of each class κ ∈ K with α κ0 > α. It is thus equal to Q y (α) = κ∈K N κ (y) max(0, ακ0-α βκ ). We say that a combination y dominates y , and we note y y, if Q y ≤ Q y on the allowed operation range [α min , α max ]. The following proposition shows that all combinations dominated by an infeasible combination are infeasible too. Proposition 1. For any combinations y and y and peak day d ∈ D P , if y y and y is infeasible for (14), so is y .

Proof. For any combination y, function Q y is strictly decreasing with the head increase, then

A y = {(q, α) ∈ R 2 | 0 ≤ q ≤ Q y (α min ), α min ≤ α ≤ Q -1
y (q)} identifies the set of pairs of total flow and head increase which can be operated by y at any time t (because we relaxed lower bound Q min κ on peak days). Then for any point in

R d (y ) = {(x , q , h ) ∈ R K T d s.t.
x d ≤ y ∀d ∈ D P }, which is built on a sequence of T d elements of A y ⊆ A y , there exists another point R d (y).

Generation of feasibility cuts

Suppose that, in the cutting-plane algorithm, a candidate combination y +1 is proved to be infeasible on some peak day. Proposition 1 implies that several feasibility cuts can be generated at the same iteration instead of a single one: it amounts to search for combinations that are dominated by y +1 . Dominance clearly defines a partial order relation on the set of combinations. The dominance of y by y can be checked by proving that Q y -Q y ≥ 0 at each point of nondifferentiability {α min , α max } ∪ {α κ | κ ∈ K}. However, the class of dominated combinations can reasonably be computed only for small values of |K| or N . We propose to compute it in a heuristic way by observing that, if p κ denotes the power of a pump κ at its maximal efficiency then the reference power P (y) = κ∈K N κ (y)p κ of y is likely greater than the reference power of the combinations it dominates. In our implementation, we first compute and sort the list C of combinations of at most N pumps by increasing order of reference power. When y is proved to be infeasible, we iterate on the combinations with lower reference power, by decreasing order, and check if they are dominated. Each dominated combination is then used to enlarge the set F in (16) indexing feasibility cuts.

As a preprocessing step for the algorithm, to accelerate its convergence, we also propose to initialize F 0 by greedily evaluating the list C. We iteratively pick up a combination y ∈ C of median value P (y) and check the feasibility of the continuous relaxation of R K T d on peak days d ∈ D P . If y is infeasible, then we updated F 0 as previously. Otherwise, we remove all y from C with P (y ) ≥ P (y) s.t. y y as they are likely feasible too.

Robustness to a pump outage

Finally, we propose to add another degree of robustness to the design problem by simulating the inoperability of one pump on peak days (due to possible failure or maintenance). To that purpose, we redefine that a combination y is feasible for (14) if all variations of y defined by setting off one pump (of any type) on peak days are feasible. This new assumption only impacts (and speeds up) the algorithms described above to initialize and update the index set F of feasibility cuts.

Numerical assessments on a real-life instance

In this section, we consider the rural drinking water distribution network FRD investigated in [START_REF] Bonvin | A convex mathematical program for pump scheduling in a class of branched water networks[END_REF]. The resulting design problem (13) were constructed with data from the historical demand profiles of the year 2013 (see [START_REF] Bonvin | A convex mathematical program for pump scheduling in a class of branched water networks[END_REF] for more details). We assume that all the 6 existing pumps should be replaced by at most N = 6 new ones selected in the KSB manufacturer catalog [START_REF] Ksb | Multitec : High-pressure pumps in ring-section design -booklet with performance curves[END_REF]. We have preselected a set K of 19 classes of pumps compatible with the network allowed range of pressure The unique peak day is built from the day in 2013 with the highest instantaneous demand and by initializing the tanks at their minimum level H min . We fix the time resolution (duration of the time steps) to 2 hours for regular days and to 4 hours for peak days.

The computations were performed on a Xeon E5-2650V4 2.2GHz with 254 GB RAM. The algorithms were implemented in Python (with tol = 10 -3 ), and the master and slave programs solved with Gurobi 7.0.2.

Numerical results

The preprocessing algorithm identified 1,902 infeasible combinations of at most 5 pumps (among 42,504) in 275 seconds, 83% of the time being to compare combinations. This resulted in an initial set F 0 of 20,757 infeasible combinations of at most 6 pumps, which represent 96% of F at the last iteration of the stabilized cutting-plane algorithm. This suggests that the preprocessing algorithm separates feasible and infeasible combinations with good accuracy.

The cutting-plane algorithm ran in 732 seconds and 37 iterations (28 infeasible candidates). In average, the time of one iteration can be roughly decomposed in 12 seconds for solving the master program, 6 seconds for the twelve regular slave programs (which could be reduced by parallelization), and 3 seconds for the peak slave program and to compute the dominated infeasible combinations.

To evaluate the quality of our approximations, we computed the non-relaxed operation cost of the optimal combination y * when optimizing on the original operation scheduling model P K T (using the heuristic in [START_REF] Bonvin | A convex mathematical program for pump scheduling in a class of branched water networks[END_REF]) on the 12 regular days of 2013 and on the 365 days of year 2013. By applying the discount rate, we get the operation costs over the 20 years: 460,388 euros in the former case and 460,247 in the latter case. This negligible deviation suggests that the chosen number of representative days is enough. By adding the investment costs on the 12-day case, we get a lifetime cost of 598,748 euros. Compared to the final lower bound z = 570, 898, it gives an optimality gap of 4.3% which confirms that the continuous relaxation R d offers a good estimate of the operation cost.

Finally, we compared the operation of network FRD before and after resizing the pumping station according to the optimal combination y * . First, the number of pumps is reduced from 6 to 5 pumps and the average reference power from 127 to 66 kW. The purchase cost of the new combination is then half the present value of the installed combination. Second, the operation costs for year 2013, when estimated with P K T , are reduced by 24.3%, a reduction that is mainly driven by a better usage of the pumps (+18.9% of the average efficiency) which proves the better adequation to the demand.

Conclusion

In this paper, we tackled the water network design problem in the context of networks equipped with pumping stations, where the operation costs must be estimated dynamically in addition to the investment costs. We formulated different approximations to estimate the operation costs to be embedded in a stabilized Benders' decomposition approach. We also ensured a certain robustness of the solutions to stress operation conditions and derived dominance arguments to accelerate the cutting-plane algorithm. Experimented on a realistic instance, the approximations turned to be accurate and the algorithm fast. While some features of our implementation are specific to the FRD network, and more generally to a class of branched network defined in [START_REF] Bonvin | A convex mathematical program for pump scheduling in a class of branched water networks[END_REF], the method can be generalized to a variety of water networks after identifying an accurate continuous relaxation for the operation problem. The potential higher complexity and larger optimality gap should however be evaluated in practice.

y

  κn ≤ N, y κn ≤ y κn-1 , for all κ ∈ K and n = 2 . . . , N.

6 0

 6 [α min = 91, α max = 140]. For each class κ ∈ K, we interpolated curves Ψ κ and Γ κ from the catalog data and estimated the investment cost I κ = 4 400 pκ 74.following ([11], Table 9-50). Considering a planning horizon of T =20 years, we built a set D R of 12 regular days, each representative of a week day or a week-end in six 2-months periods (January-February, March-April,...) of any year. For each regular day d ∈ D R , we define L d = L 0 d 20 l=1 (1 + τ ) 1-l with L 0 d the number of days represented by d in 2013 and τ = 5% the discount rate. Demand and tariff profiles D d and C d for a regular day d are computed in average over the L 0 d represented days of 2013.