Antoinette Maarawi

Zoé Anxionnaz-Minvielle

Pierre Coste

Nathalie Di

Miceli Raimondi

Michel Cabassud

Numerical Investigation of the Heat Transfer in a Meandering Millichannel by a coupled 1D-Flow-Channel / 3D-Surrounding-Solid Model

Keywords: heat exchanger/reactor, modelling, heat transfer, exothermic reaction, process intensification, corrugated channel

Heat exchanger/reactors with corrugated channels are used to implement especially exothermal chemical reactions. Heat conduction effects in their solid material are not negligible and induce in particular heat flows between adjacent channels. In this context, this work proposes a model of a platetype heat exchanger/reactor with a 2D-meandering millichannel, taking into account these conduction effects by coupling a 1D approach for the channel flow, with correlations for the friction factor and the Nusselt number, to a 3D approach for the solid part. This model is a compromise between the 1D classical reactor models that don't account for 3D conduction effects in the solid and the 3D simulations of the channel flow with the surrounding solid which are computationnally expensive or unachievable for a large reactor. In the tested case, the 1D/3D and 3D CFD results are in good agreement (5 min vs 4 h computational time), which suggests valid model assumptions and internal model consistency.

Simulated and experimental conversion rates of a fast highly exothermic reaction are then compared, showing the relevancy of the model. Further investigations are performed to highlight the effect of conjugate heat transfer in millistructured heat exchangers and to quantify its limitations versus the 1D/3D model.

INTRODUCTION

Process intensification has been defined in several ways in the literature [START_REF] Rivas | Process intensification education contributes to sustainable development goals. Part 1[END_REF]- [START_REF] Van Gerven | Structure, Energy, Synergy, Time -The Fundamentals of Process Intensification[END_REF]. However, all the definitions focus on the multidisciplinary nature of this engineering approach with the goal of developing innovative technologies that offer processes that are more environmentally friendly and energy efficient, as well as processes that are more productive and operate safely. Process intensification may lead to design multifunctional and compact units such as compact heat exchangers/reactors. The combination of two unit operations in the same miniaturised apparatus (heat exchange and reaction) enables process intensification. Different kinds of compact heat exchangers/reactors (HEX reactors) exist and are summarized in the review of Anxionnaz et al. [START_REF] Anxionnaz | Heat exchanger/reactors (HEX reactors): Concepts, technologies: State-of-the-art[END_REF]. They are promising for the implementation of chemical syntheses, especially fast and exothermic ones, by offering efficient heat and mass transfer and temperature control. However, carrying out chemical reactions requires a sufficient residence time for the chemistry. Fluid velocities are though moderate, resulting in laminar to transitional flow regime conditions in miniaturised reactors. Therefore, in order to maintain a plug flow and intensify heat and mass transfer, generating vorticity in the flow regions is a key success factor. This can be achieved using several methods that exist in the literature, such as jets [START_REF] Nastase | Impinging jets -a short review on strategies for heat transfer enhancement[END_REF] and turbulence promoters or vorticity generators for generating embedded vortices [START_REF] Ghanem | Energy efficiency in process industry -High-efficiency vortex (HEV) multifunctional heat exchanger[END_REF], [START_REF] Li | A flapping vortex generator for heat transfer enhancement in a rectangular airside fin[END_REF], and wall curvature [START_REF] Anxionnaz-Minvielle | Influence of the meandering channel geometry on the thermo-hydraulic performances of an intensified heat exchanger/reactor[END_REF]- [START_REF] Schwolow | Design and application of a millistructured heat exchanger reactor for an energy-efficient process[END_REF] for inducing secondary flows (Dean vortices). This last technique leads to implementing corrugated geometries in compact HEX reactors. Plate HEX reactors with wavy channels exist in the literature and are commonly used due to their compactness and modularity regarding the number of plates. Examples for such devices are the LFR and AFR from Corning [START_REF] Wu | Hydrodynamic Study of Single-and Two-Phase Flow in an Advanced-Flow Reactor[END_REF], the Marbond from Chart Martson [START_REF] Phillips | Development of a novel integrated chemical reactor-heat exchanger[END_REF], the Flow Plate from Lonza and Ehrfeld Mikrotechnik [START_REF] Mielke | Micro-reactor mixing unit interspacing for fast liquid-liquid reactions leading to a generalized scale-up methodology[END_REF], the OPR from Alfa Laval Vicarb [START_REF] Prat | Performance Evaluation of a Novel Concept "Open Plate Reactor" Applied to Highly Exothermic Reactions[END_REF], the MR500 from 3M Technical Ceramics [START_REF] Schwolow | Model-based scale-up and reactor design for solvent-free synthesis of an ionic liquid in a millistructured flow reactor[END_REF], the ART plate reactor PR37 from Ehrfeld Mikrotechnik GmbH [START_REF] Rave | Characterization of a Modular, Scalable Millistructured Plate Reactor[END_REF] and the DeanHex from the Laboratoire de Génie Chimique (LGC, France) and the Commissariat à l'Energie Atomique et aux Energies Alternatives (CEA, France) [START_REF] Anxionnaz-Minvielle | Influence of the meandering channel geometry on the thermo-hydraulic performances of an intensified heat exchanger/reactor[END_REF]. Their corrugated geometry allows the generation of secondary motions (i.e. Dean vortices), that are responsible for flow mixing and heat and mass transfer enhancement. Therefore, these micro/millistructured heat exchanger/reactors are preferably used to conduct highly exothermic reactions that cannot be carried out in a batch reactor without high risk of thermal runaway.

The present work focuses on the "DeanHex" plate HEX reactor that consists of a stack of process plates (where the reaction takes place) and utility plates (for cooling or heating), where millichannels are engraved by laser machining (cf. Figure 1). Experiments conducted on the "DeanHex" reactor showed a shift to a lower Reynolds number (equal to 200) of the transition between laminar and turbulent flow compared to straight channels thanks to the corrugation of the reaction channel [START_REF] Anxionnaz-Minvielle | Influence of the meandering channel geometry on the thermo-hydraulic performances of an intensified heat exchanger/reactor[END_REF], [START_REF] Théron | Characterization of the performances of an innovative heat-exchanger/reactor[END_REF]. Heat and mass transfers are thus intensified and plug flow behaviours are obtained while working under flow conditions that usually lead to laminar flow regime.

The number of plates is dependent on the application, in order to provide the required residence time to complete the chemistry. The geometry of the reaction channel consists of periodic zigzag units and is represented in Figure 2. In this work, a square process channel with a hydraulic diameter dh of 2 mm, a radius of curvature of the bends R of 1.5 mm, a 7 mm long straight section L between two bends and an angle θ between two straight sections equal to 90 °, is considered. Even if the use of micro/millistructured heat exchanger/reactors is an interesting alternative to "batch" systems, the number of parameters to consider for their design (hydrodynamics, heat and mass transfer, etc.) makes it difficult to scale them up to industrial scale. In this context, numerical simulation is a useful tool to study the system and analyse the influence of channel dimensions on the flow during the scale-up process.

Since the flow in micro-and millichannels is of the plug flow type, classical 1D models based on the same assumptions as those used for modelling conventional continuous reactors can be developed [START_REF] Nauman | Chemical Reactor Design, Optimization, and Scaleup[END_REF].

However, these models do not generally consider all the heat conduction effects in the solid material which are not negligible in micro-and millistructured heat exchangers/reactors (high solid/fluid volume ratio) and induce in particular heat flows between adjacent channels. Indeed, in these intensified devices (cf. Figure 1), the "process" channel meanders through the "process" plate to give the sufficient fluid residence time. The conductive materials used and the proximity of the adjacent channels because of compactness mean that there are heat flows between these channels, not only between the "process" channel and the "utility" channel. In order to increase the compactness of heat exchangers/reactors, the distance between adjacent channels is decreased, reducing the thermal resistance due to the walls. Heat is not only transported from one fluid perpendicularly through the wall to the other fluid, but mainly in all directions (three-dimensional heat flow). This represents the so-called "conjugate heat transfer" and must be taken into account, especially in micro/millistructured compact devices. Boundary conditions such as a constant temperature on the channel wall or a constant heat flux on the wall no longer apply and the heat transfer in the solid-fluid system must be calculated simultaneously with the flow field [START_REF] Rave | Modeling and experimental investigation of the spatial heat transfer in a plate reactor with meandering millichannels[END_REF], [START_REF] Dai | Impact of tortuous geometry on laminar flow heat transfer in microchannels[END_REF]. To our knowledge, there is only one 1D modelling work that has addressed these conduction effects, in which a model for the PR37 plate ART heat exchanger/reactor designed by Ehrfeld Mikrotechnik GmbH was developed [START_REF] Rave | Modeling and experimental investigation of the spatial heat transfer in a plate reactor with meandering millichannels[END_REF]. This model consists in calculating the heat flow between a given position in the process channel and all sections of the utility channel using a thermal resistance network. The only adjustable parameters of the model are the parameters of the Nusselt correlations used to calculate the heat transfer resistances in the channels. According to the authors, this model allows to simulate the heat transfer in the reactor with few computational resources and with high accuracy over a wide range of operating conditions.

Furthermore, throughout the years, three-dimensional (3D) CFD simulations have been conducted for corrugated micro-and millichannels, for a better representation of the local physico-chemical phenomena. However, Table 1 shows that the simulated domain, in several heat transfer numerical studies from the literature, carried out for other corrugated channels than the DeanHex, consists of "only" one to a few periodic units to minimize the computational time and resources. Nevertheless, the number of mesh elements required to conduct the simulations is still high, making the simulations long and complex. In fact, the complex geometries of the channels (especially the presence of bends) with a high L/dh ratio and the large number of units which may not be all periodic (cf. Figure 1) would make 3D CFD simulations highly demanding in staff and computational resources. For the "DeanHex" reactor, 3D CFD simulations of 3 to 11 periodic units (depending on the case) of the heat transfer in the square zigzag millichannels has been carried out in the recent study of Shi et al. [START_REF] Shi | Numerical study of heat transfer in square millimetric zigzag channels in the laminar flow regime[END_REF].

The simulations were conducted in stationary laminar flow conditions, using the finite volume CFD code, ANSYS CFX 16. Several channel geometries (Ls=2-12 mm) with a total length of 0.1 m were modelled. In terms of meshing, respectively 4.3 million (at Re=224) and 7.6 million (at Re=560) mesh elements were needed for the wavy channels consisting of 8 and 4 periodic units.

Consequently, for the purposes of rapid simulations of case studies and of building a general model for the millistructured heat exchanger/reactor (DeanHex) that is able to predict the performance of this device in different applications and at scales up to the industrial one, this work deals with an 1D/3D numerical appoach which aims at:

 in comparison with the 1D classical approaches of reactor modelling: taking into account the thermal conduction effects in the solid material which are not negligible in the DeanHex (high solid volume to fluid volume ratio) and induce, in particular, heat transfer between adjacent channels; and providing a finer representation of the heat transfer process, with local heat transfer coefficients from wall to fluid used in the simulations;

 in comparison with the 3D CFD simulations: providing a less resource demanding simulation tool for millichannel geometry with a high L/dh ratio.

The full description of this approach is detailed in a previous work [START_REF] Maarawi | 1D/3D Numerical Approach for Modelling a Milli-Structured Heat Exchanger/Reactor[END_REF]. In the present work, the validation of the model is addressed in two parts: first with 3D CFD simulations and then with previous experimental works conducted in the "DeanHex" zigzag millichannel with a fast and exothermic chemical reaction. In each part, a comparison of results is made in order to discuss the relevance and the accuracy of the proposed model.

NUMERICAL APPROACH 2.1 Principle

The 1D/3D numerical approach [START_REF] Maarawi | 1D/3D Numerical Approach for Modelling a Milli-Structured Heat Exchanger/Reactor[END_REF] used in the present work consists in representing the meandering flow channel by a 1D model coupled to a 3D model for the surrounding solid in which the channel is embedded. The 1D model considers all process variables but with reduced calculation time compared to 3D approaches, which are generally focused on local effects. In the 1D model, calculated variables are averaged over the cross-sections of the channel and discretised along the channel; heat is exchanged in all normal directions to its walls with the surrounding solid. Correlations, related to thermal performance and friction laws for example, that have been obtained during experimental works are integrated in this 1D model. Correlations may also be obtained with a validated 3D CFD approach on a limited number of periodic units. The 1D model is coupled to a 3D model for the surrounding solid in order to take into account the thermal conduction effects in the heat transfer process simulation.

Consequently, this approach not only avoids meshing the channel cross-section with a complex 3D mesh, but also provides a simulation method capable of predicting the performance of the heat exchanger/reactor (conversion rate, selectivity, etc...). It intends to bridge the gap between the simulation at small scale and the simulation at the reactor scale, at a reasonable CPU cost. The numerical resolution of this 1D/3D model is performed using COMSOL Multiphysics V5.5.

Computational domain

The three-dimensional zigzag channel is represented by a one-dimensional curve "C" following the axial direction of the flow in the channel. The channel is embedded in a solid matrix to form the reactor; thus in the COMSOL model, "C" is surrounded by the 3D solid "S". However, doing so, the geometry of the solid is not well represented which will lead to the application of non-realistic boundary conditions in the "S" domain, in particular in terms of heat flux. The solid material volume to fluid volume ratio should be conserved in order to properly represent the heat flux in the equipement. Consequently, the proposed solution is to define, around the curve "C", the three-dimensional channel as a highly radially conductive fictive solid "S2" [START_REF] Maarawi | 1D/3D Numerical Approach for Modelling a Milli-Structured Heat Exchanger/Reactor[END_REF], as shown in Figure 3. The idea is to transmit in each channel cross section the temperature from "C" to the walls of "S2", consistently with the cross section averaging of the 1D fluid approach.

The computational domain thus consists of a one-dimensional zigzag channel "C" of total developed length equal to L, where the process fluid flows; a three-dimensional zigzag solid domain "S2" of total developed length equal to L; and a three-dimensional solid bloc "S" surrounding "S2", which represents the solid material in the heat exchanger/reactor, with a thickness t, width w and length l. A sensitivity study was carried out in order to consistently determine the material properties of "S2" [START_REF] Maarawi | 1D/3D Numerical Approach for Modelling a Milli-Structured Heat Exchanger/Reactor[END_REF]. In order to be able to set the radial thermal conductivity, a curvilinear coordinate system was defined for "S2". The axial (resp. radial) components of the "S2" thermal conductivity are then set along this curvilinear coordinate system to zero (resp. to 10 5 W‧m -1 ‧K -1 ) [START_REF] Maarawi | 1D/3D Numerical Approach for Modelling a Milli-Structured Heat Exchanger/Reactor[END_REF].

Model equations

The 1D simulation work is based on the resolution of classical continuity and momentum equations for steady state and an incompressible fluid flowing inside the channel [START_REF] Maarawi | 1D/3D Numerical Approach for Modelling a Milli-Structured Heat Exchanger/Reactor[END_REF]:

ρ u • ∇u = -∇p -Ʌ ρ 2d , u|u| (1)
∇ • (S ρ u) = 0 ( 2
)
where u is the cross section averaged velocity (m‧s -1 ), ρ the fluid density (kg‧m -3 ), p pressure (Pa), Ʌ (dimensionless) the Darcy friction factor, d the mean hydraulic channel diameter (m) and S the channel cross section area available for flow (m 2 ). The subscript p refers to the process side.

The Darcy friction factor Ʌ is calculated from Théron et al. [START_REF] Théron | Characterization of the performances of an innovative heat-exchanger/reactor[END_REF]'s correlation, established for the zigzag millichannel, in function of the Reynolds number Re as shown in equation ( 3). The Reynolds number in the channel is calculated using equation ( 4), where μ is the fluid viscosity (Pa‧s). (3)

Re = ρ • u • d μ (4)
Let et = (et,x et,y et,z) be the unit tangent vector to the channel axis (Figure 4). The 1D simulation work solves for the tangential velocity up that is defined as follows:

u = u • e (5)
In addition, the solved energy equation for an incompressible fluid flowing in the channel is the following:

ρ S C , u • ∇T = ∇ • S λ ∇T + Ʌ ρ S 2d , |u| + Q (6)
Where Tp is the process fluid temperature (K), λ is the fluid thermal conductivity (W‧m -1 ‧K -1 ), C is the heat capacity at constant pressure (J‧kg -1 ‧K -1 ) and Qwall represents external heat exchange through the channel wall (W‧m -1 ). This radial heat transfer from the surroundings into the channel is given by:

Q = h P (T -T ) ( 7
)
where P is the wetted perimeter of the channel (m), h is the convective heat transfer coefficient (W•m - 2 •K -1 ) inside the channel and TS the external temperature outside the channel (K).

h is deduced from the Nusselt number correlation (8) established for the zigzag millichannel of the DeanHex heat exchanger/reactor [START_REF] Anxionnaz | RAPIC project: Toward competitive heat-exchanger/reactors[END_REF] for 550 < Rep < 8623. Pr is the Prandtl number in the process channel, calculated using equation [START_REF] Rave | Characterization of a Modular, Scalable Millistructured Plate Reactor[END_REF].

Nu = h • d , λ = 0.16 • Re . • Pr . ( 8
)
Pr = C • μ λ (9)
Ts in equation ( 7) is given as the temperature field computed in the 3D surrounding solid. In this way, heat transfer coupling is done between 1D channel flow and the 3D side as a line source. In other words, the temperature coupling between the channel and the surrounding domain is implemented as a line heat source in the 3D domain. The source strength is proportional to the temperature difference [START_REF] Li | A flapping vortex generator for heat transfer enhancement in a rectangular airside fin[END_REF] between the channel fluid and the surrounding domain.

In this work, in the case of implementing a chemical reaction that takes place in a homogeneous phase, the source term Q , which is the heat released (W•m -1 ) by the reaction, is added to the energy balance equation ( 6) and a mass transfer equation ( 10) is added to the above-mentioned equations resolved by the 1D model, where C and R are the concentration (mol•m -3 ) and the rate of consumption (mol

•m -3 •s -1 )
of the limiting reactant, respectively.

u • ∇C = R (10)
In addition, a heat balance equation is resolved for the 3D surrounding solid, taking into account the conduction in the solid material as well as the convective heat transfer through the solid S walls [START_REF] Maarawi | 1D/3D Numerical Approach for Modelling a Milli-Structured Heat Exchanger/Reactor[END_REF]:

∇ • (-λ ∇T ) = q (11)
where λ is the thermal conductivity (W‧m -1 ‧K -1 ) of the 3D solid, T the solid temperature (K) and q the convective heat flux (W•m -2 ) at the solid S walls [START_REF] Maarawi | 1D/3D Numerical Approach for Modelling a Milli-Structured Heat Exchanger/Reactor[END_REF].

The validation of the proposed approach is done by comparing the obtained 1D/3D results, firstly, to 3D

CFD calculations, secondly, to experimental results. Finally, in order to show the importance of the 3D solid representation, a sensitivity study about the influence of the solid material thermal conductivity is conducted, as well as a comparison of 1D/3D results to 1D calculations.

3D CFD VS 1D/3D SIMULATIONS

Straight channel

In order to validate the proposed 1D/3D approach in the simple case of a straight channel, a comparison is made with a full 3D simulation using ANSYS Fluent CFD code.

Figure 5: Schematic of the simulated geometry using the 1D/3D approach for the straight channel configuration.

Computational domain and meshing

The computational domain consists of a square cross-section (2×2 mm 2 ) 1 meter long straight channel, embedded in a parallelepipedic (w=12 mm, t=4 mm and l=1 m) solid matrix. The latter is made of structural steel with a thermal conductivity of 44 W‧m -1 ‧K -1 . The simulated geometry using the 1D/3D approach is shown in Figure 5, where the one-dimensional straight channel is added (which is no more than the pipe axis).

For the 3D CFD simulations, the cross-sectional mesh is made of cells with a typical size of about 3×10 - 4 m in the center of the channel which get finer (five cell layers) close to the walls in order to ensure the most accurate resolution at the wall (cf. Figure 6a). The surrounding solid is meshed three times coarser than the center of the fluid. The 3D CFD mesh is thus made of 1 030 000 cells, while in the 1D/3D calculations, it consists of only 60 000 cells (cf. Figure 6,b). Grid independence studies, for both simulations, were conducted to ensure that the solution, in terms of thermal profile, is not influenced by the size of the grid.

Equations and boundary conditions

For both simulations, at the channel inlet, constant velocity and temperature conditions are assumed (equal to 0.069 m‧s -1 and 20 °C respectively). At the channel outlet, an average pressure of 1 atm is set.

The fluid is heated through its passage in the channel due to a constant temperature boundary condition (equal to 60 °C) on the solid external walls.

In the 3D CFD simulations, a no-slip condition is applied to the channel walls and a steady laminar flow is considered. Continuity, momentum (Navier-Stokes) and energy equations for an incompressible fluid are solved. In the 1D/3D model, the following Nusselt correlation [START_REF] Nellis | Heat transfer[END_REF] (equation ( 12)) is integrated into the 1D model, where x is the axial coordinate of the channel, in order to take into account the fact that the flow is not established close to the inlet. Note that Nu = 2.98 corresponds to a fully developped laminar flow in channels of square cross-sections, with a uniform temperature boundary condition on the four external channel walls [START_REF] Shah | Rectangular Ducts[END_REF].

Nu = 2.98 + 0.049 + 0.020 Pr • Re • Pr • d , x . 1 + 0.065 • Re • Pr • d , x . (12)

Results

The computed temperature field in the 1D/3D calculation corresponds to average temperature values in the channel cross sections, along the channel axis (x-axis). In the 3D CFD calculation, the average temperature at axial location x is calculated, during post-processing, from the computed velocity v (m•s - 1 ) and temperature T (K) fields, as follows:

T = ∬ T • v • da ∬ v • da (13)
The average temperature profiles for both calculations are represented in Figure 7. It can be observed that the 1D/3D model describes well the heat transfer in the channel. The maximum relative error in temperature between the 1D/3D and 3D CFD simulations is less than 7%, which is acceptable. This difference may result from the use of the above-mentioned Nu correlation [START_REF] Phillips | Development of a novel integrated chemical reactor-heat exchanger[END_REF] in the 1D/3D simulation because the condition of a constant temperature at the walls of the channel is not verified. Instead, it is the temperature at the external walls of the solid S that is fixed in the simulation. Therefore, the Nusselt correlation may introduce an inaccuracy in this simulated case.

Zigzag Channel

To go further, more complexity is added to the model. A corrugated geometry is now considered. It consists of the 2×2 mm 2 zigzag channel used by Théron et al. [START_REF] Théron | Characterization of the performances of an innovative heat-exchanger/reactor[END_REF]. The 1D/3D simulated thermal profiles are compared to 3D CFD results obtained using the finite volume CFD code ANSYS Fluent 2020 R1.

Computational domain and meshing

For the comparison with 3D CFD results, only a part of the 2×2 mm 2 zigzag process channel is modelled since the full 3D modelling of the zigzag channel is very expensive in computational resources. The process channel length L is reduced to about 0.1 m, which represents roughly the first row of the original process channel shown in Figure 1, a.

The computational domain consists of a zigzag channel that starts and ends at half a straight section (cf.

Figure 8) where the fluid flows in the 3D CFD simulations and the one-dimensional zigzag channel is added (which is no more than the three-dimensional zigzag channel axis) in the 1D/3D simulations; and a zigzag solid matrix (t=w=6 mm and l=0.1 m), surrounding the three-dimensional channel, representing the solid material of the process and closing plates (the utility plates are not modelled). The solid material is the stainless steel with a thermal conductivity of 16.3 W‧m -1 ‧K -1 . For the 3D CFD simulations, in order to provide an efficient and accurate resolution, a refined meshing is constructed for the channel, as shown in Figure 9, a. In the channel cross-section where the fluid flows, the mesh size is 10 -4 m at the center and 2×10 -5 m near the wall, with five layers. The solid matrix is meshed five times coarser than the center of the fluid. The mesh consists of about six million cells.

For the 1D/3D mesh, a free tetrahedral meshing is applied to the 3D solid domains where only thermal conduction takes place. There is no need for a finely constructed meshing (cf. Figure 9,b). The mesh consists of 205 681 cells. Grid independence studies were conducted to ensure that the solution, in terms of thermal profile, is not influenced by the size of the grid.

Equations and boundary conditions

For the 3D CFD calculations, the SST k-ꞷ turbulence model is used. In fact, in the simulated cases, the Reynolds number varies between 700 and 2250. Numerical studies in zigzag geometries similar to the one studied in the present work show that the flows can become unsteady at Reynolds number around 500 [START_REF] Shi | Numerical study of heat transfer in square millimetric zigzag channels in the laminar flow regime[END_REF], [START_REF] Zheng | Transient laminar heat transfer simulations in periodic zigzag channels[END_REF]. All simulations are at steady state. Continuity, momentum (Navier-Stokes) and energy equations for an incompressible fluid are solved.

In a previous work [START_REF] Maarawi | 1D/3D Numerical Approach for Modelling a Milli-Structured Heat Exchanger/Reactor[END_REF], the consistency of the 1D/3D numerical approach for the simulation of heat transfer in the DeanHex was shown by comparison with experimental data from [START_REF] Théron | Characterization of the performances of an innovative heat-exchanger/reactor[END_REF]. In the present study, the same conditions are used to compare the two types of simulations (1D/3D vs 3D CFD).These experimental conditions are presented in Table 2. Qp and Qu (kg‧h -1 ) are the mass flowrates of the process and utility fluids, respectively. Note that the high utility flowrate guarantees that the heat transfer process is not limited by the local heat transfer coefficient on the utility side (the thermal resistance on the utility side is 3 to 6 times lower than that on the process side). At the channel inlet, a constant temperature Tp,in (°C) is set. At the outlet, an atmospheric pressure condition is set.

Since the utility channels are not modelled, a heat flux boundary condition "qS" ( 14) is imposed on the upper and lower external walls of "S" to represent the cooling side. "qS" is equal to zero at the remaining sides of "S". Figure 10, b illustrates this representation in the 1D/3D model. The heat flux boundary condition "qs" is expressed as follows:

q = h • (T -T ) • (14)
where the subscript u refers to the utility side. T is the utility fluid temperature (K). It is assumed uniform and constant along the channels. This hypothesis is confirmed by the inlet and outlet temperature measurements on the utility side in [START_REF] Théron | Characterization of the performances of an innovative heat-exchanger/reactor[END_REF] (relative difference less than 4%) where the utility fluid flowrate is 8 to 23 times higher than the process fluid flowrate.

The convective heat transfer coefficient h (W•m -2 •K -1 ) in the zigzag 2×2 mm 2 utility channels, shown in Figure 1, b, is deduced from the following Nusselt correlation [START_REF] Shi | Numerical study of heat transfer in square millimetric zigzag channels in the laminar flow regime[END_REF]:

Nu = h • d , λ = 0.2 • (Re . + 8.9) • Pr . • d , L , . (15)
where Re and Pr are the Reynolds and Prandtl numbers for the utility side. d , and L , are, respectively, the hydraulic diameter (m) and the straight length (m) between the bends in the utility zigzag channels. The physicochemical properties of the utility fluid (water), such as its thermal conductivity λ (W•m -1 •K -1 ), are estimated at T (K).

However, when representing the utility channels of heat transfer surface area equal to A (m 2 ) by a heat flux boundary condition on the wall of "S" of surface area equal to A (m 2 ), a correction factor equal to should be applied, in order to take into consideration the difference in heat transfer areas [START_REF] Maarawi | 1D/3D Numerical Approach for Modelling a Milli-Structured Heat Exchanger/Reactor[END_REF].

Utility

Results

For the 3D CFD simulations, the average temperature in the channel cross section is calculated from the computed temperature and velocity fields using equation ( 13). The average temperature profiles for both calculations are presented in Figure 11 in function of the channel length. A maximum relative error in temperature of less than 2% highlights a very good agreement between 1D/3D and 3D CFD simulations. Once the results consistency between 1D/3D and 3D CFD is shown, it is interesting to point out the following observation that highlights the CPU related interest of the 1D/3D model. A 0.1 m long channel requires 6 million cells to be modelled using 3D CFD approach leading to up to 4 hours lasting simulations, whereas only 205 681 cells are used in the 1D/3D simulations that converge in less than 5 min (CPU used: Intel(R) Xeon(R) Silver 4214 CPU @ 2.20GHz 2.19 GHz (2 processors)). Assuming that the heat transfer correlations are known, the proposed model in this paper allows a significant reduction of computational resources compared to full 3D simulations, which is of interest when it comes to modelling the full-scale heat/exchanger reactor, with a channel total developed length of 6.6 m.

Therefore, based on the comparison between 1D/3D and 3D CFD simulations for the straight and zigzag channels, the performance of the 1D/3D model proves to be acceptable. In the next section, experimental results from Théron et al. [START_REF] Théron | Characterization of the performances of an innovative heat-exchanger/reactor[END_REF] are used as a supplementary validation of the proposed approach.

EXPERIMENTS VS 1D/3D SIMULATIONS 4.1 Description of the heat exchanger/reactor used in the experiments

The experimental work of Théron et al. [START_REF] Théron | Characterization of the performances of an innovative heat-exchanger/reactor[END_REF] is about a DeanHex reactor at pilot-scale. This HEX reactor is made of stainless steel (thermal conductivity equal to 16.3 W‧m -1 ‧K -1 ) and consists of a stack of three process plates (P1, P2 and P3) and four utility plates (U1, U2, U3 and U4), which configuration is represented in Figure 12. 2 mm square cross-section channels, with periodic zigzag units, were engraved on both plates by laser machining (Figure 1). The process and utility plates are sandwiched between 2 mm thick metal plates (closing plates) (Figure 10,a). In the work of Théron et al. [START_REF] Théron | Characterization of the performances of an innovative heat-exchanger/reactor[END_REF], experiments were performed with the implementation of a highly exothermic oxidation reaction. The latter is added in the 1D/3D model and simulations of the heat transfer with chemical reaction are conducted in order to make comparisons with the experimental results and to validate the proposed model.

Utility Inlet Utility Outlet

Process Inlet

Process Outlet U1 U2 U3 U4 P1 P2 P3

Computational domain and meshing

The first process plate of 2.2 m long zigzag channel can be modelled as shown in Figure 13 [START_REF] Maarawi | 1D/3D Numerical Approach for Modelling a Milli-Structured Heat Exchanger/Reactor[END_REF], where the configuration "C-S2-S" is established. However, in order to have the same residence time of the reactants into the reactor as in the experiments, the process channel of the whole reactor (6.6 m of total developed length) should be considered. Therefore, three geometries like the one illustrated in Figure 13 are placed in series to have the 6.6 m long channel for the reaction. The free tetrahedral mesh consists of 1 068 720 cells for the 6.6 m long zigzag channel.

Equations and boundary conditions

For the 1D/3D simulations, the experimental conditions listed in Table 3 are employed. Tp,in is applied at the inlet of "C" and an atmospheric pressure at its outlet. The cooling side is represented by a heat flux boundary condition "qS" ( 14) on the upper and lower walls of the solid S, as described previously.

Tu in ( 14) is equal to the average of the inlet and outlet temperatures of the utility fluid obtained during each experiments and listed in Table 3. Adiabatic boundary conditions are assumed on the other sides of the solid S.

Reaction kinetics

The considered reaction is the oxidation reaction of sodium thiosulfate (Na S O ) by hydrogen peroxyde (H O ) which is written:

2Na S O + 4H O → Na S O + Na SO + 4H O
For this first order reaction with respect to both reactants, the reaction rate "r" is expressed as a function of the reactants' concentrations, as follows:

r = k • C • C (16)
where the kinetic constant k is assumed to be governed by an Arrhenius law, as follows:

k = k • e ℛ• (17)
The kinetic parameters of equation ( 17), obtained by Grau et al. [START_REF] Grau | Batch and semibatch reactor performance for an exothermic reaction[END_REF] (k0=8.13 × 10 8 m 3 ‧mol -1 ‧s -1 , Ea/ℛ=9156 K) are used in the simulations, since they best fit the experimental results of Théron et al. [START_REF] Théron | Characterization of the performances of an innovative heat-exchanger/reactor[END_REF]. The reaction takes place in a homogeneous liquid phase, therefore the mixture fluid properties and total process flowrates are considered in the 1D model. The consumption rate of the limiting reactant R in equation ( 10) is R = -2r. The 1D model solves the mass balance equation for C . Thus, C is replaced by the following expression that is deduced from the mass balance in the reactor, taking into consideration the stoechiometric coefficients of the reaction:

C = C , -2 • C , -C (18) where C , and C
, are the concentrations of sodium thiosulfate and hydrogen peroxide (mol.m -3 ) at the reactor inlet.

This irriversible and fast reaction is highly exothermic with a heat of reaction ∆H =-586.2 kJ•mol -1 of Na S O . The heat generated (W•m -1 ) during the reaction, integrated in the 1D model, is thus equal to:

Q = R • ∆H • S (19)
Where Sp (m 2 ) is the cross section of the channel.

Results

Since in the experiments the high utility flowrate guarantees that the heat transfer process is not limited by the local heat transfer coefficient on the utility side (the thermal resistance on the utility side is 3 to 6 times lower than that on the process side), the comparison of the experimental results with the 1D/3D simulation results will show the effect of the proposed approach on the process side. Therefore, referring to Théron et al. [START_REF] Théron | Characterization of the performances of an innovative heat-exchanger/reactor[END_REF]'s experiments for the validation of the 1D/3D approach is suitable.

The simulated and experimental results are compared in terms of outlet process temperature and conversion rate, as shown in Table 4. The maximum relative difference in outlet process temperature (

) is less than 8%, which is acceptable.

Since hydrogen peroxide is in excess during the reaction, the conversion rate obtained in the simulations is calculated according to the concentration loss of sodium thiosulfate. Théron et al. [START_REF] Théron | Characterization of the performances of an innovative heat-exchanger/reactor[END_REF] have calculated the conversion rate using two different methods:

 Reactor: based on a thermal balance between the inlets and outlets of the process and utility streams at steady state;

 Dewar: based on measuring the adiabatic temperature rise in the Dewar vessel from the sampling until reaching the equilibrium temperature.

Both methods induce some errors in the calculation of the rate of conversion of the sodium thiosulfate.

The differences in temperature, used in the thermal balance, between the inlet and outlet of the utility fluid (cf. Table 3) may fall in the experimental error of the temperature measurements and thus the conversion rate calculated thanks to the first method may be inaccurate. This may justify the differences in conversion rate obtained in the simulations compared to the experimental "reactor" conversions, especially in experiment n°1 where Tu,out-Tu,in=0.2°C. Regardless this experiment, the relative difference in conversion ( ) is less than 12%.

Furthermore, the sampling time in the Dewar method may not be negligible in comparison with the residence time of the reactants in the reactor, which will lead to an overestimation of the conversion rate at the reactor outlet, as seen in Table 4 when comparing simulated to experimental "Dewar" conversions.

Nevertheless, the relative difference in conversion (

| | ) is less than 5%.
Overall, in terms of outlet process temperature and conversion rates, the simulations seem to succeed in reproducing the experimental results. Therefore, it could be deduced that the 1D/3D model with the radially infinite thermal conductivity imposed for "S2" proposed in this work is valid for the DeanHex reactor considering these experimental data.

N° Exp. It is also interesting to visualize the temperature distribution in the solid material of the reactor plate, which is one of the assets of the proposed model in particular when comparing to 1D classical models.

This is done for a cut plane that corresponds to an orthogonal 2D coordinate system embedded in the 3D geometry. This plan passes through the middle of the fictive solid "S2". The temperature distribution obtained when simulating experiment n°4 (Table 3) is illustrated in Figure 14. The solid temperature varies with the process fluid temperature which is plotted in Figure 15 (SS plot).

The reactants enter the channel at 20.7°C and are heated by the utility side (Tu=50.1°C) and by the heat generated by the occurring reaction. The process temperature increases gradually and reaches 60°C.

Then, it decreases towards Tu due to the cooling effect of the utility side. This temperature variation is observed in Figure 14 from the solid perspective. Furthermore, temperature gradients are observed between adjacent channels. Far from the process channel, the solid temperature is close to the utility temperature.

In the next section, the actual influence of the conjugate heat transfer on the thermal performances of the heat exchanger/reactor is discussed.

INFLUENCE OF CONJUGATE HEAT TRANSFER

Sensitivity study about the solid material conductivity

In order to better understand the influence of conjugate heat transfer in the case of millistructured heat exchanger/reactors, a sensitivity study about the influence of the solid material thermal conductivity has been conducted. Therefore, a negligible thermal conductivity (0.01 W‧m -1 ‧K -1 ) is tested to simulate the case without heat transfer between channels and three other typical values of material conductivities for an intensified heat exchanger: glass (1 W‧m -1 ‧K -1 ), stainless steel (16.3 W‧m -1 ‧K -1 ) and silicon carbide (170 W‧m -1 ‧K -1 ).

In the previously described simulated case (section 4), the reactor solid material is stainless steel. Hence, the obtained results when simulating the experiment n°4 (Table 3) are used in the sensitivity study for the stainless steel case. Considering the same conditions, new simulations are done by only changing the thermal conductivity of the solid material to 0.01; 1 and 170 W‧m -1 ‧K -1 . Four temperature and conversion profiles are shown in Figure 15 to demonstrate the influence of solid material conductivities on the thermal behaviour of the reactor. Simulated fluid temperatures are plotted as lines over the channel length and reaction conversions are shown as symbols. When increasing the conductivity of the solid, the temperature inside the channel is better controlled. The heat generated by the chemical reaction is dissipated in the solid and the temperature peak is reduced. Whereas in the case of an insulator, the process temperature is out-of-control and continues to increase as shown in Figure 15. In addition, the reaction kinetics depend on the process temperature (equation ( 17)), thus the effect of the solid conductivity is also observed on the obtained conversions. This highlights the effects of heat conduction in the solid material on the heat exchanger/reactor performances. In addition, in terms of heat flux distribution in the solid material, Figure 16 shows the y-axis component of the conductive heat flux for the SiC case on the cut plane that passes through the middle of the fictive solid "S2". This y-axis component of the heat flux illustrates the heat flow between adjacent channels. This heat conduction flux between adjacent channels is indeed clearly observed, especially between channel rows 1-2 and 2-3. This flow magnitude may reach 0.5×10 5 to 2×10 5 W•m -2 and hence, it is not negligible.

Comparison with the 1D classical model results

A temperature profile is calculated assuming a one-dimensional heat flow between utility and process channels and is then compared to the temperature profile of the 1D/3D model.

Basic heat transfer without chemical reaction is considered. Therefore, on the one hand, the first process plate of the "DeanHex" is simulated using the 1D/3D model by applying the conditions listed in Table 2 (Qp=6.9 L•h -1 ). On the other hand, based on a 1D classical plug-flow model, a temperature profile is calculated from the energy balance over an element dl (m) of the process channel length (equation ( 20)).

dT dl = 4 • 1 h + 1 h + e λ • T -T Q • C , (20)
Where e (m) is the solid thickness between process and utility plates.

The obtained thermal profiles are shown in Figure 17. Differences between both temperature profiles are observed for each tested solid material. This shows that the heat is not only transported from the process fluid perpendicularly through the wall to the utility fluid as it is the case with 1D calculations. Axial heat conduction effects exist and are taken into account in the 1D/3D model. When decreasing the conductivity of the solid material to nearly 0, these effects are less significant and thus the 1D/3D and 1D temperature profiles match better. This result shows once again the effect of the conduction in the solid material of the heat exchanger/reactor and justifies the idea of this work that is to take into account the axial heat conduction in the channel walls.

CONCLUSIONS

This paper deals with a new 1D/3D numerical approach proposed in order to model the DeanHex reactor at a reasonable CPU cost despite the complexity of the channel geometry and its high L/dh ratio. The model is based on the coupling of a 1D representation of the flow channel with a 3D surrounding solid.

In this way, the solid to fluid volume ratio in the reactor is conserved and the heat transfer process is well simulated, especially when implementing extremely high exothermic chemical reactions. Firstly, the heat exchange without chemical reaction has been investigated in the DeanHex reactor using the 1D/3D model. The consequent simulated thermal profiles have been firstly compared to 3D CFD simulations (less than 2% relative difference and significant computational time reduction (5min vs 4h)).

Secondly, a highly exothermic reaction was implemented and the simulated results were compared to experimental ones in terms of outlet temperature (less than 8% relative difference) and conversion rates (less than 12% relative difference). The proposed model provides a precise computation of the conversion rate from the concentration loss of reactants, which is a crucial parameter for the reactor.

Based on these comparisons, the 1D/3D approach has proved to be promising and valid for the DeanHex reactor. Furthermore, the effect of conjuguate heat transfer is investigated by comparing 1D/3D results to 1D calculations and by conducting a sensitivity study about the thermal conductivity of the solid material. Both investigations showed the presence of conjugate heat transfer in heat flows between adjacent channels and their influence on the millistructured heat exchanger/reactor performances. These results highlight the importance of taking into account the conduction effects in the solid material of the reactor in the numerical model. Moreover, assuming the availability of heat transfer and friction laws, either from experiments or from dedicated validated 3D CFD calculations on a restricted number of periodic units, the 1D/3D model could be used for other corrugated channels, usually with high L/dh and solid volume to fluid volume ratios, exchanging heat in all directions with the surrounding medium.

Other reactions that occur in a homogeneous phase can also be modelled using the approach described in this paper. The next step is thus to extend the use of this model for two-phase flows (heterogeneous reactions and mass transfer limitations) as well as in dynamic mode to evaluate the effects of operating condition deviations on reaction control regarding the material thermal capacitance. At the end, the aim is to build a generalized scaling-up model for the purpose of predicting, at industrial scale, the performance of the heat exchanger/reactor in terms of conversion rate and energy efficiency.

Figure 1 :

 1 Figure 1: DeanHex heat exchanger/reactor pilot: (a) process plate and (b) utility plate [16].

Figure 2 :

 2 Figure 2: Schematic of the zigzag channel with geometric parameters [17].

Figure 3 :

 3 Figure 3: Illustration of the simulated geometry using the 1D/3D approach.

Ʌ = 24 . 3 •

 243 Re . for 20 < Re < 200 Ʌ = 6 • Re . for 200 < Re < 2250

Figure 4 :

 4 Figure 4: Unit tangent vector to the channel axis.

Figure 6 :

 6 Figure 6: Illustration of the mesh in the cross-section of the a) 3D CFD and b) 1D/3D simulated geometry.

Figure 7 :

 7 Figure 7: Comparison between 1D/3D and 3D CFD simulations in terms of mean temperature profile of the fluid for a straight channel.

Figure 8 :

 8 Figure 8: The modelled portion of the channel in the 1D/3D and 3D CFD simulations.

Figure 9 :

 9 Figure 9: Illustration of the mesh for the process channel and the surrounding solid used in the a) 3D CFD simulations and b) 1D/3D simulations.

Figure 10 :

 10 Figure 10: Schematic of: (a) the plates' configuration in the reactor with a zoom-in on the channel inlet and (b) its representation in the 1D/3D model.

Figure 11 :

 11 Figure 11: Evolution of the 1D/3D and 3D CFD simulated process fluid temperature along the process channel for different process flowrates.

Figure 12 :

 12 Figure 12: Process and utility plates configuration in Théron et al. [16]'s pilot.

Figure 13 :

 13 Figure 13: Illustrations from the 1D/3D model of: (a) the first process plate of Théron et al.[22]'s pilot with the inlet and outlet of the flow and (b) a zoom-in on the process inlet where the C-S2-S chain is showed [17].

Figure 14 :

 14 Figure 14: Temperature distribution in the solid material corresponding to experiment n°4.

Figure 15 :

 15 Figure 15: Thermal profile (lines) and reaction conversion (symbols) in the channel for different plate solid materials.

Figure 16 :

 16 Figure 16: Y-axis component of the conductive heat flux distribution in the solid material corresponding to experiment n°4.

Figure 17 :

 17 Figure 17: 1D/3D vs 1D process temperature profiles along the channel length for different solid materials.

NOMENCLATURE S cross-sectional area (m 2 )A heat exchange area (m 2 )Λ 1 ) 2 ) 1 ) 3 )

 221213 Cp specific heat of the fluid (J‧kg -1 ‧K -1 ) dh hydraulic diameter (m) t thickness of solid matrix (m) l solid matrix length (m) w solid matrix width (m) e solid thickness between process and utility plates (m) Darcy friction factor h heat transfer coefficient (W‧m -2 ‧K -1 ) L channel developed length (m) Ls straight section length (m) Q volumetric flow rate (m 3 ‧s -Nu local Nusselt number Pr Prandtl number q heat flux (W‧m -Rc radius of curvature (m) Re Reynolds number T temperature (K) v velocity vector (m‧s -1 ) u cross-section mean fluid velocity (m‧s -1 ) up tangential velocity (m‧s -1 ) Qr heat generated by the reaction (W‧m -1 ) C concentration of the reactant (mol‧m -3 ) r reaction rate (mol.m -3 .s -1 ) R consumption rate of the reactant (mol.m -3 .s -Ea activation energy (J‧mol -1 ) k reaction rate constant (m3‧mol -1 ‧s -1 ) k0 specific rate constant (m 3 ‧mol -1 ‧s -1 ) ℛ universal gas constant (J‧mol -1 ‧K -1 ) da element of cross-sectional area (m 2 ) x coordinate over the channel length (m) y coordinate over the channel width (m) z coordinate over the channel depth (m) et unit tangent vector Abreviations Na2S2O3 sodium thiosulfate H2O2 hydrogen peroxide Na2SO4 sodium sulfate H2O water 1D one-dimensional 2D two-dimensional 3D three-dimensional CFD computational fluid dynamics Greek symbols ϴ angle between two straight sections (°) λ thermal conductivity (W‧m -1 ‧K -1 ) μ dynamic viscosity (Pa‧s) ρ fluid density (kg‧m -ΔHr heat of reaction (kJ‧mol -1 ) Subscripts p process side u utility side s solid matrix in channel inlet out channel outlet Na2S2O3 relative to the reactant sodium thiosulfate H2O2 relative to the reactant hydrogen peroxide simu relative to simulation exp relative to experiment

Table 1 :

 1 Simulated domains and number of mesh elements used in 3D CFD simulations of corrugated channels.

	References	Year of the simulation	Number of simulated periodic units	Total number of computational volume elements
	Facão and Oliveira[25]	2005	1	833 000
	Rosaguti et al.[22]	2007	1	210 000
	Sui et al. [26]	2012	1	136 653
	Zheng et al.[23]	2013	10	10 000 000
	Zheng et al.[24]	2013	14	13 000 000
	Karale et al.[27]	2013	1	450 000
	Abed et al.[28]	2015	10	2 340 000
	Premkumar D et al.[21]	2016	3	2 030 188

Table 2 :

 2 Operating condition applied for the simulations.

		side		Process side
	Qu (kg‧h -1 )	Tu (°C)	Qp (kg‧h -1 )	Rep	Tp,in (°C)
			2.4	554
			4	937
	152	15.6			76
			5.5	1303
			6.9	1641

Table 3 :

 3 Experimental conditions for simulating the heat exchange with chemical reaction.

Table 4 :

 4 Comparison between simulated and experimental results for the reaction of sodium thiosulfate oxidation by hydrogen peroxide.

		Théron et al. [16]		Simulation (this work)
		Tp,out,exp (°C)	Conversion Reactor Dewar	Tp,out,simu (°C)	Conversion	Residence time (s)
	1	43.9	60	59	40.6	73	6.9
	2	41.4	82	94	40.2	92	19.3
	3	43.4	88	91	40.7	86	13.8
	4	51.0	93	100	50.3	96	13.8

ACKNOWLEDGEMENTS

This research did not receive any specific grant from funding agencies in the public, commercial, or notfor-profit sectors.


Abdul Rochim
email: abdoul.rochim@univ-nantes.fr.

Didier Marot

Luc Sibille

Fateh Bendahmane

The Effect of Fine Content and History of Hydraulic Loading on the Characterization of Suffusion Sensibility of Cohesionless Soil

Keywords: Cohesionless soils, erodimeter, suffusion, water seepage power

Suffusion involves a selective erosion of fine particles in a soil by the water seepage in the matrix formed by coarser particles. This phenomenon can induce a change in particle size distribution, porosity and hydraulic conductivity of the material. With the objective to characterize the soil sensibility and propose suffusion sensibility classification, a series of one-dimensional downward seepage flow tests was realized with an erodimeter. Tests were performed under controlled hydraulic gradients in binary mixtures of sandy gravel soils. Two controlling parameters were used: 1) fine content and 2) hydraulic loading history. An analysis based on energy induced by the seepage flow is proposed to characterize the hydraulic loading and the cumulative eroded dry mass to characterize the soil response. The results show that for a given soil different hydraulic loading histories result in different soil sensibilities. Soils with the lowest fine contents tend to require larger flow energy for the development of erosion. These results demonstrate that this approach is effective to characterize suffusion sensibility for cohesionless soils.

Introduction

During the cycle of seasons, hydraulic structures made of soil may be in contact with variable interstitial flows, possibly generating the detachment and the transport of some constituent particles of the structure or their foundations. Different categories can be distinguished in this problem generally is called internal erosion: concentrated leak erosion, suffusion, contact erosion and backward erosion [START_REF] Fell | The state of the art of assesing the likelihood of internal erosion of embankment dams, water retaining structures and their foundations, Internal Erosion of Dams and their Foundations[END_REF]. Suffusion corresponds to selective erosion of fine particles by seepage flow leaving behind the matrix of coarser particles. However, only soil sensibility classifications for interface erosion were proposed so far [START_REF] Hanson | Erodibility of cohesive streambeds in the loess area of the midwestern U.S.A[END_REF][START_REF] Wan | Investigation of rate of erosion of soils in embankment dams[END_REF][START_REF] Marot | Energy based method for providing soil surface erodibility rankings[END_REF], and classification of soils with respect to suffusion sensibility is still missing.

Literature Review

Suffusion is a complex matter as it depends on geometric, mechanical and hydraulic condition. The size of fine fraction should be smaller than the size of constrictions, in addition the volume of fine particles should be less than the volume of void formed by the coarse particles. This can result in criteria on the shape of the grading curves [START_REF] Lafleur | Filtration of broadly graded cohesionless soils[END_REF][START_REF] Wan | Assessing the Potential of Internal Instability and Suffusion in Embankment Dams and Their Foundations[END_REF]. Hydraulic conditions should be studied to investigate the onset of suffusion even if the transport of particles is geometrically possible. The hydraulic loading applied on a soil is often described by the hydraulic gradient, the pore velocity or hydraulic shear stress. Recently a new method has been proposed in [START_REF] Marot | Energy based method for providing soil surface erodibility rankings[END_REF] to describe the hydraulic loading from the energy dissipated by the fluid seepage.

The stream power related to water seepage, and called here "erosion power", Perosion can be expressed by: ' '

where Jw is the unit weight of water, 'P = PA -PB, the pressure drop; 'Z = ZA -ZB, the altitude change for a one dimensional flow between an inlet section A and an outlet B; and Q the fluid flow rate. The erosion energy, Perosion, is defined as the time integration over the test duration of Perosion.

When the erodibility, or erosion sensibility, is estimated with respect to the hydraulic shear stress (assumed to be representative of the hydraulic loading), an erodibility coefficient indicates the erosion rate for a unit excess of shear stress through an erosion function [START_REF] Hanson | Channel erosion study of two compacted soils[END_REF][START_REF] Wan | Investigation of rate of erosion of soils in embankment dams[END_REF]: [START_REF] Hanson | Erodibility of cohesive streambeds in the loess area of the midwestern U.S.A[END_REF] where is the soil erosion rate; kd the erodibility coefficient; W the hydraulic shear stress at the soil-water interface; and Wc the critical shear stress at onset of erosion.

It is worth stressing that Equation 2 has been developed only for interface erosion. However in case of suffusion, the expression of the hydraulic shear stress within a soil defined in [START_REF] Reddi | Comparison of internal and surface erosion using flow pump test on a sand-kaolinite mixture[END_REF] can be used. For a vertical downward flow the shear stress is given by:

P (3)
where 'h is the hydraulic head drop, k the hydraulic conductivity, P the viscosity, and n the porosity.

Downward Seepage Test

Tested Gradations and Specimens

Three grain size distributions corresponding to sandy gravel soils were tested (Figure 1). Properties of gradations are summarized in Table 1. The grain size distributions are assessed Table 1 with respect to several gradation-based criteria [START_REF] Wan | Assessing the Potential of Internal Instability and Suffusion in Embankment Dams and Their Foundations[END_REF][START_REF] Kezdi | Soil physics selected topics[END_REF][START_REF] Kenney | Internal stability of granular filters[END_REF][START_REF] Li | Comparison of two criteria for internal stability of granular soil[END_REF]. The definitions of H/F and d15/d85 are given in [START_REF] Kezdi | Soil physics selected topics[END_REF] and [START_REF] Kenney | Internal stability of granular filters[END_REF] respectively. Cu: uniformity coefficient; Gr: gap ratio (see [START_REF] Chang | Extended internal stability criteria for soils under seepage[END_REF]); P: percentage of particle smaller than 0.063 mm U: internally unstable gradation, S: internally stable gradation [START_REF] Wan | Assessing the Potential of Internal Instability and Suffusion in Embankment Dams and Their Foundations[END_REF][START_REF] Kezdi | Soil physics selected topics[END_REF][START_REF] Kenney | Internal stability of granular filters[END_REF][START_REF] Li | Comparison of two criteria for internal stability of granular soil[END_REF] are gradation criteria: Wan and Fell, Kezdi, Kenney and Lau, Li and Fannin respectively

Erodimeter

An erodimeter as shown in Figure 2 was used to characterize the soil sensibility. It consists of an erosion cell, a demineralised water supply system, a soil collection system, and a water collection system. The sample is saturated by flushing water in the upward direction, and erosion test is then performed by injecting through the sample in the downward direction from the injection tank. The funnel-shaped draining system is connected to effluent tank by a glass pipe. The effluent tank is equipped with an overflow outlet in order to control the downstream hydraulic head and a rotating sampling system containing 8 beakers for catching eroded particles. Overflow water is continuously weighed to determine flow rate. A series of tests was conducted in two steps: specimen preparation (production of the specimen, and saturation) and downward seepage test. The sand grains and gravel are then first mixed with a fixed moisture content. The specimens are then prepared using a single layer semi-static compaction technique. The mixture is placed in a mould of 50 mm diameter and 50 mm height and subsequently compressed under the action of two pistons until the initial fixed dry density is reached. In the erosion cell device, the specimen is placed on a 4 mm pore opening grid and wrapped within membrane. This pore opening allowed the migration of all particles of sand. For the saturation phase carbon dioxide is first injected in the specimen for duration of 5 minutes, before

Results and Discussion

Erosion rate and Hydraulic conductivity

Erosion rate and hydraulic conductivity can be interpreted as a soil response to hydraulic loadings. Figure 3 shows the evolution of hydraulic conductivity of tested specimen during test and Figure 4 depicts the variation of erosion rate versus hydraulic shear stress. The arrow signs in both figures highlight the time when hydraulic conductivity reaches constant and significant eroded mass is captured. The evolution of hydraulic conductivity present a typical trend: a decrease is firstly observed before it progressively increases and finally reaches a constant value from the arrow sign as shown in Figure 3. We assume this decrease of the hydraulic conductivity is attributed to some detached particles that were transported under the imposed water seepage and then filtered further within the soil itself. This filtration leads to partial clogging corresponding the first decreasing of the hydraulic conductivity. Hydraulic conductivity increases only latter, for much larger hydraulic gradients able to detach again these filtrated particles.

From the variation of hydraulic conductivity of specimens subjected to a multistage hydraulic gradients (A-90a, A-90b, B-90a, B-97a, B-97b, R2-90a, R2-90b) three predominant processes can be identified: filtration, process of erosion and finally steady state (represented by constant value of hydraulic conductivity). However such successive phases do not seem to occur when a single stage hydraulic gradient equal to 4. Considering test B-90a, by increasing progressively the applied hydraulic gradient from 1.3 to 4 the hydraulic conductivity increases by a factor of 20. Whereas in the case of test B-90c, performed with a single-stage hydraulic gradient, even by applying hydraulic gradient of 4, the hydraulic conductivity continuously decreases. Thus for the need to the classification of suffusion sensibility, the method by increasing hydraulic gradients is preferred.

Soil sensibility

In conformity with methods proposed for interface erosion classification, a first approach to define a suffusion sensibility classification can consist to investigate the variation of the erosion rate with the hydraulic shear stress. It is worth stressing that the erosion function proposed by [START_REF] Wan | Investigation of rate of erosion of soils in embankment dams[END_REF] and [START_REF] Hanson | Channel erosion study of two compacted soils[END_REF] has been developed only for interface erosion. Erosion rate (kg/s/m 2 )

Times (min) Hydraulic conductivity (m/s)

A-90a A-90b A-90c B-90a B-90c B-97a B-97b R2-90a R2-90b
A-90a A-90b A-90c B-90a B-90c B-97a B-97b R2-90a R2-90b
Erosion coefficient, kd

A-90c hydraulic gradient, i = 1.3 hydraulic gradient, i = 4 B-90c
The hydraulic shear stress within the soil is computed from Equation 3 and the erosion rate coefficient kd is assumed to be representative of the erosion sensibility. In the necessity to evaluate the classification of suffusion sensibility of tested specimens a first approach proposed by [START_REF] Wan | Investigation of rate of erosion of soils in embankment dams[END_REF] for interface erosion can be used. This approach is based on the determination of kd coefficient. In order to apply such approach for suffusion process, the value of kd is defined from linear approximation of the plot of the erosion rate with respect to the critical hydraulic shear stress. The critical hydraulic shear stress is determined at the point where the hydraulic conductivity presents a minimum or where the hydraulic conductivity starts to inflect. For instance, for specimen B-90a this point corresponds to a hydraulic gradient equal to 1.3 as shown in Figure 3 and4. The difficulty to determine kd coefficient is encountered for a single-stage hydraulic gradient (A-90c and B-90c) as shown in Figure 4. Thus these specimens are not taken into account for proposition of suffusion sensibility classification.

The approach based on the erosion energy, as introduced in Section 2, and corresponding to the energy expanded by the water to seep through the soil, is now considered. Figure 5 shows the relationship between cumulative expanded energy and cumulative eroded mass. The arrow signs indicate the point from which the hydraulic conductivity is stabilized.

Classification of suffusion sensibility

The proposed classification is built from the approach based on the erosion energy (Figure 5). This classification is divided into five classes from resistant to highly erodible as presented in Figure 6. The points highlighted by arrows correspond to the end of development of suffusion (i.e. time from which hydraulic conductivity is stabilized) in order to classify the soil sensibility. It is worth noting that the points before and after the arrows do not represent the soil sensibility. The eroded mass in the beginning may be attributed to loss mass during saturation phase and during filtration process. The after-arrow points also do not represent the soil sensibility since they

1,0E-02 1,0E-01 1,0E+00 1,0E+01 1,0E+02 1,

Effect of hydraulic loading history and fine contents

Figure 5 shows that soils A, B and R2 imposed to hydraulic gradient a are less erodible than those imposed to hydraulic gradient b. It can be also shown in Figure 6, given B-90 and R2-97 different hydraulic loading histories give different classifications. Thus hydraulic loading history can be an important parameter to soil sensibility. With respect to fine content, it seems that a specimen with a low fine content soil, specimens A tends to require a larger energy for the onset of erosion than a specimen with the same initial density but a higher fine content, specimens B. Indeed, since the same initial global density is obtained for specimens A and B, but with a larger content of coarse particles A, density proper to the coarse fraction should be higher in A than in B. Therefore constriction sizes of the coarse skeleton may be smaller for specimen A, limiting the possibility of erosion of the fine fraction. This result is in good agreement with the test results presented in [START_REF] Ke | Strength reduction of cohesionless soil due to internal erosion induced by onedimensional upward seepage flow[END_REF]. Given three different fine contents: 16.7%, 20%, and 25%, it is demonstrated that the samples with the lowest fine content (16.7%) required a larger critical hydraulic gradients for the onset of internal erosion for relative density of 0.2 and 0.6 respectively. These results show that the energy-based approach should be effective to distinguish suffusion sensibility for cohesionless soils.

Conclusion

In this study, the effect of history of hydraulic loading and fine content on the characterization of suffusion sensibility of cohesionless soil was investigated through a series of tests using erodimeter. When a soil is susceptible to erosion, the soil sensibility has to be characterized by imposing seepage flow in suffusion tests. The test should be performed by progressively increasing the applied hydraulic gradient and it should be carried on until the hydraulic conductivity reaches a steady state. With the objective to characterize independently the hydraulic loading and the induced erosion, the cumulative eroded dry mass and the cumulative energy expanded by the seepage flow, Eerosion are computed. Finally with the aim to classify suffusion sensibility, the first time hydraulic conductivity reaches constant value can be determined as the time to suffusion sensibility classification of a soil on the diagram of the cumulative eroded dry mass versus the cumulative expanded energy. For given soils, it demonstrates different soil sensibility classification for different hydraulic loading history and fine content. The gradient of erosion coefficient for hydraulic gradients with increment 1 (hydraulic gradient b) shows more erodible than ones with increment 0.1 (hydraulic gradient a). This result shows how hydraulic loading is an important parameter for characterization of suffusion.

Figure 1 .

 1 Figure 1. The particle size distribution of tested specimens

  water. The whole saturation phase required approximately one night. Finally, the specimen is subjected to a hydraulic flow in a downward direction using demineralized water.A series of nine tests with two different hydraulic loadings was performed: multistage (a and b) and a single-stage (c) hydraulic gradients. The first (a) hydraulic loading consisted of increasing the hydraulic head by steps from 0.1 to 2, then by steps of 0.5 between 2 and 4 and by steps of 1 beyond; whereas steps were equal to 1 for the second kind of hydraulic loading (b). For each step the hydraulic gradient was kept constant during 10 minutes. However c corresponds to a single-stage hydraulic gradient equal to 4 m/m.

Figure 2 .

 2 Figure 2. Schematic diagram of the erodimeter

Figure 3 .

 3 Figure 3. Time series of hydraulic conductivity

Figure 4 .

 4 Figure 4. Erosion rate versus hydraulic shear stress

Filtration

Figure 5 .

 5 Figure 5. Cumulative eroded dry mass in terms of the cumulative energy expanded by the water to seep through the soil.

  steady state. If the series of final points is taken into account, the soil sensibility becomes less erodible.

Figure 6 .

 6 Figure 6. Classification of suffusion sensibility based on the energy expanded by the water seepage in the soil

Table 1 .

 1 Tested gradations properties and assessment of soil internal stability by recent criteria

	Tested	Cu	Gr	P	(H/F)min	d (H/F=min) d15/d85		Assessment method
	gradation			(%)		(mm)		[9]	[10]	[11]	[6]
	A	17.06	2.14	1.227	0.038	0.400	8.761	U	U	U	U
	B	19.52	2.14	1.533	0.035	0.400	8.741	U	U	U	S
	R2	24.46	WG	1.200	0.165	0.212	9.653	U	U	U	S

Acknowledgement

The authors thank The Indonesian Directorate General of Higher Education (DIKTI), Sultan Agung Islamic University Indonesia, GeM Institute -University of Nantes-ECN-CNRS for providing financial support.