R Gantois

A A⇤

Cantarel

G Dusserre

J.-N Félices

F Schmidt

NUMERICAL SIMULATION OF RESIN TRANSFER MOLDING USING BEM AND LEVEL SET METHOD

Keywords: Resin Transfer Molding (RTM), Level Set, Boundary Element Method (BEM)

Resin Transfer Molding is widely used to produce fiber-reinforced materials. In the process, the resin enters a close mold containing the dry fiber preform. For mold designer, numerical simulation is a useful tool to optimize the mold filling, in particular to identify the best positions of the ports and the vents. An issue in mold filling simulation is the front tracking, because the shape of the resin front changes during the flow. In particular, topological changes can appear resulting from internal obstacles dividing the front or multi-injection. A previous approach [1] using the Boundary Element Method (BEM) in a moving mesh framework shows the capability of the method to compute accuratlely the front propagation at low CPU time. The present paper describes a method developed to handle complex shapes, using BEM together with a Level Set approach. Numerical results in two dimensions are presented, assuming a Newtonian non-reactive fluid, and an homogeneous and not-deformable reinforcement. The resin flow in the fibrous reinforcement is modeled using Darcy's law and mass conservation. The resulting equation reduces to Laplace's equation considering an isotropic equivalent mold. Laplaces equation is solved at each time step using a constant Boundary Element Method to compute the normal velocity at the flow front. It is extended to the fixed grid and next used to feed a Level Set solver computing the signed distance to the front. Our model includes a boundary element mesher and a Narrow Band method to speed up CPU time. The numerical model is compared with an analytical solution, a FEM/VOF-based simulation and experimental measurements for more realistic cases involving multiple injection ports and internal obstacles.

INTRODUCTION

Liquid Composite Molding (LCM) processes are widely used in industry. Among them, Resin Transfer Molding (RTM) is one of the most popular. It consists in injecting the liquid resin in the dry preform held in position in a mold. Recent aerospace programms focus on forming structural parts using LCM. As mechanical performances strongly depend on filling conditions, resin flow prediction is important in mold design. In particular gates' locations are adjusted so that the resin impregnates correctly the entire preform. In that task, numerical simulation can be usefull. As the resin flows, a tracking technique is employed to follow the moving front. It can be performed using moving mesh methods [START_REF] Schmidt | Numerical simulation of resin transfer molding using linear boundary element method[END_REF][START_REF] Um | A study on the mold filling process in resin transfer molding[END_REF][START_REF] García | Meshless methods with application to liquid composite molding simulation[END_REF], Volume Of Fluid methods [START_REF] Trochu | Numerical analysis of the resin transfer molding process by the finite element method[END_REF][START_REF] Hirt | Volume of fluid (vof) method for the dynamics of free boundaries[END_REF] or Level Set methods [START_REF] Sethian | Level Set Methods and Fast Marching Methods Evolving Interfaces in Computational Geometry[END_REF][START_REF] Soukane | Application of the level set method to the simulation of resin transfer molding[END_REF]. The main advantage of the last one is that the front is accurately captured in an Eulerian ⇤ Corresponding author: Postal address: Ecole des Mines Albi, Campus Jarlard, F-81013 Albi, France. Phone : (+33) 5.63.49.30.00. Email address : renaud.gantois@mines-albi.fr framework. The present paper focus on a technique combining a Boundary Element Method(BEM) [START_REF] Schmidt | Numerical simulation of resin transfer molding using linear boundary element method[END_REF][START_REF] Um | A study on the mold filling process in resin transfer molding[END_REF][START_REF] Soukane | Application of the level set method to the simulation of resin transfer molding[END_REF][START_REF] Yoo | Numerical simulation of the resin transfer mold filling process using the boundary element method[END_REF] and a Level Set method.The first part of this paper considers the governing equations. The second part describes the implemented model. The last part covers some applications.

GOVERNING EQUATIONS

PRELIMINARY TRANSFORMATIONS

Our model assumes that the liquid resin is a Newtonian non-reactive fluid, flowing in isothermal conditions through an homogeneous and not deformable fibrous reinforcement. In the impregnated area ⌦, the macroscopic resin motion is governed by the modified Darcy's law [START_REF] Simacek | Role of acceleration forces in numerical simulation of mold filling processes in fibrous porous media[END_REF] and incompressibility equation

< ! v >= [K] µ ! rp 0 (1) ! r. < ! v >= 0 (2) 8 > > < > > :
where < ! v > is the macroscopic velocity, [K] the permeability tensor, µ the liquid resin viscosity and p 0 the acting pressure. That pressure given by p 0 = p + ⇢gz includes a gravity term, where ⇢ is the resin specific mass and g gravity. Combining the previous equations and transforming the coordinates in the isotropic equivalent domain ⌦ e of permeability K e , we obtain

! v = K e µ✏ ! rp 0 (3)
4p 0 = 0 (4) 8 > < > :
where ! v is the resin velocity, computed from the macroscopic velocity using ✏ (medium porosity), and K e is the equivalent isotropic permeability. The isotropic equivalent transformation [START_REF] Schmidt | Numerical simulation of resin transfer molding using linear boundary element method[END_REF][START_REF] Adams | Radial penetration of viscous liquid into a planar anisotropic porous medium[END_REF] is

" x e y e # = p K e " 1 p K 1 0 0 1 p K 2 #" x y # (5)
where K 1 and K 2 are the principal permeabilities of the reinforcement. The isotropic equivalent permeability, is given as

K e = p K 1 K 2 .

BOUNDARY CONDITIONS

Let us consider e the isotropic equivalent front bounding ⌦ e . Boundary conditions are assigned according to the location of the point M 2 e under consideration. Dirichlet conditions (imposed pressure) are prescribed on the gate ep 0 and on the free edge

ep f p 0 = ( p 0 + ⇢gz 0 8M 2 ep 0 p f + ⇢gz f 8M 2 ep f (6)
and Neuman conditions (imposed normal pressure gradient) are prescribed on the mold wall eq 0 for a nonpenetration condition

! rp 0 . ! n = 0 8M 2 eq 0 (7)
where ! n is the unit outwards vector at point M .

NUMERICAL METHOD

OUTLINE

Our programm is implemented using Matlab environement. At the beginning, pre-processing imports a standard mesh file, performs the isotropic equivalent transformation, assigns material and processing data, and locates the injection ports. Next, the filling routine based on a Level Set formulation advances the front. At the end, the post-processing plots the results once the real domain is recovered.

The filling stage is divided into a finite number of quasi steady states. At each step time, governing equations are solved using a constant Boundary Element Method. More details can be found in section 3.2. The front is advanced by feeding a Level Set solver with the BEM-computed velocities. How it is done is covered in section 3.3. It is repeated until the mold is completely filled.

BOUNDARY ELEMENT METHOD

For clarity, the subscript e referring to the isotropic equivalent domain is omitted by the next. As mentioned earlier, BEM confines the calculation on the front which is a close curve bounding the calculation domain ⌦. A known value of pressure p or normal pressure gradient q is prescribed on the edges p or q . Laplace's equation is multiplied by the Green function p ? . Integrated twice by parts over the calculation domain and using Green's theorem leads to Somigliana's equation [START_REF] Schmidt | Numerical simulation of resin transfer molding using linear boundary element method[END_REF][START_REF] Brebbia | Boundary elements: an introductory course[END_REF][START_REF] Mathey | Optimisation numérique du refroidissement des moules dinjection de thermoplastiques basée sur la simulation des transferts thermiques par la méthode des éléments frontières[END_REF]

c i p i + ˆ pq ? d = ˆ p ? qd with c i = ✓ 2⇡ (8)
where p i is the value of the pressure at a point M i on the boundary, q ? the pressure gradient associated with p ? and ✓ the internal angle of the corner in radians. For a twodimensional domain, p ? and q ? are given as

p ? = 1 2⇡ ln ✓ 1 r ◆ and q ? = 1 2⇡ ! r . ! n r 2 (9)
where r is the distance from the point M i of application of the Dirac delta function to any point under consideration.

Meshing the boundary into N constant boundary elements and applying Equation 8 leads to

c i p i + N X j=1 p j ˆ j q ? d = N X j=1 q j ˆ j p ? d (10)
Rewritten in a matrix form using H and G, Equation 10 is then transformed into

N X j=1 p j H ij = N X j=1 q j G ij (11)
where

H ij = 1 2 ij +
´ j q ? d and G ij =

´ j p ? d are N 2 matrix. Finally, the previous equation is reordered to take the form of a linear system AX = F , where X is a vector of N unknowns.

LEVEL SET METHOD

The moving front is captured using the signed-distance function , defined so that the zero level set corresponds to the interface [START_REF] Sethian | Level Set Methods and Fast Marching Methods Evolving Interfaces in Computational Geometry[END_REF][START_REF] Sethian | Level set methods for fluid interfaces[END_REF]

(t) = (x, y) 2 R 2 / (x, y, t) = 0 (12)
For each point under consideration, distance is signed negative if located in the impregnated area and positive otherwise. We consider two meshes : a fixed grid made of unstructured triangle elements and a moving mesh made of beam elements. Level Set equations [START_REF] Sethian | Level Set Methods and Fast Marching Methods Evolving Interfaces in Computational Geometry[END_REF][START_REF] Sethian | Level set methods for fluid interfaces[END_REF] govern the evolution of the signed-distance function. To combine with BEM, we use a formulation involving a propagating interface with a velocity in its normal direction, given as follows

@ @t + F | ! r | = 0 (13) (x, y, t = 0) = 0 (14) 8 > < > :
where F is the extended normal velocity, built to coincide with the velocity of the front (it is extrapolated outside whithout physical meaning). The filling stage is initialized by choosing 0 as the signed-distance to the inlet gates.

Next, values for each points are updated using an Euler scheme

(x, y, t + t) = (x, y, t) F | ! r | t ( 15
)
where t is the time step, adjusted to match its upper limit (CFL conditions).

For the sake of efficiency, is not updated on the entire grid, but only on few nodes around the front in a "narrow band" [START_REF] Sethian | Level Set Methods and Fast Marching Methods Evolving Interfaces in Computational Geometry[END_REF]. Next, front is rebuilt using a boundary element mesher, based on interpolating the zero level set on the grid.

The method directly handles topological changes, involving merging or dividing fronts, but does not ensure that the resin remains inside the mold. Contact with mold is implemented using a fixed level set describing the mold walls. It acts by correcting as follows

= max( , mold ) ( 16
)
where mold is the signed-distance to the mold.

APPLICATIONS 4.1 ISOTROPIC RADIAL INJECTION

We consider here the standard case of an isotropic radial injection, for a 0.25 m 2 square plate. The resin enters the part from the center and flows through the preform describing circular patterns. Four vents located on the corners ensure the complete filling. Material and processing parameters are summerized in Tables 1 and2. In the first case, gravity is neglected for a comparison with an analytical solution [START_REF] Schmidt | Numerical simulation of resin transfer molding using linear boundary element method[END_REF][START_REF] Adams | Radial penetration of viscous liquid into a planar anisotropic porous medium[END_REF]. It is given as

Table 1: Material data K[m 2 ] µ[P a.s 1 ] ⇢[kg.m 3 ] ✏[ ] 1.10 9 0.1 1150 0.5
p 0 [P a] p f [P a] r 0 [m] g[m.s 2 ] 2.
 r f r 0 2 .  2 ln ✓ r f r 0 ◆ 1 + 1 = 4K (p 0 p f ) t ✏µr 2 0 ( 17
)
where r f and r 0 are the radii of the moving front and the inlet gate, p 0 and p f the inlet and outlet pressures. The mold is meshed using 2500 triangle elements. CPU time is around 20 s on a 2.26 GHz / 1.93 Go of RAM laptop. The predicted filling time is 98 s. Our model is assessed by fitting the numerical fronts using the analytically predicted circles. Figure 1 shows the comparison between our numerical results (left) and the analytical results (right). Analytical results (in bold lines) are overlaid on our results, for 10.7 s and 63.3 s elapsed time. Using L1 norm leads to the relative error of 2.1%, which shows the accuracy of our numerical model. The second case takes into account gravity effects. On Figure 2 we compare our results (left) with a FEM/VOFbased code (PAM-RTM TM ) (right) using the same mesh.

As expected, the flow is slightly deported in the lower part of the plate. The predicted filling time is 101 s for our simulation and 89 s for PAM-RTM TM , which shows a good qualitative accordance. The difference is due to the contact implementation which is different in Level Set (see Equation 16) and VOF. Further developements will improve the contact accuracy by refining the mesh on the mold wall.

ANISOTROPIC INJECTION INVOLVING COMPLEX SHAPES

We compared our numerical results with experimental data, performing an infusion experiment involving a more realistic case. It consisted in impregnating an anisotropic knitted glass preform (1x1 rib knit fabric manufactured by Textile Aero Tarn) using a canola oil. The reinforcement (one layer) was placed under a transparent flexi-ble bag and filled using two inlet gates and two vaccum lines (on the upper and lower sides). The dimensions of the mold are 0.3 m per 0.35 m. Some internal obstacles were placed to simulate contacts with mold walls. Material data (see reference [START_REF] Dusserre | Effect of deformation on knitted glass preform in-plane permeability[END_REF] for measurement procedure) and processing parameters are given in Tables 3 and4. CPU time is around 200 s for a model involving 3693 elements. Figure 3 shows the comparison between our results (at the top) and experimental data (at the bottom) at different times.The agreement is fair at any time. In particular, fronts merging (at 3 s) and dividing on the obstacles (at 7 s) are accurately predicted.

CONCLUSION

We developed a software to predict a two-dimensional resin impregnation for both isotropic and anisotropic cases. Results were compared with an analytical solution, a FEM/VOF-based simulation and experimental data, with a fair agreement. The implemented model includes a gravity term, but it can easily be modified to take into account other body forces. Further developements will include a 3D approach (technique remains unchanged), and a high permeability layer to simulate infusion process.

Figure 1 :

 1 Figure 1: Comparison with analytical solution

Figure 2 :

 2 Figure 2: Comparison with FEM/VOF (PAM-RTM TM )

.10 3 Figure 3 :

 33 Figure 3: Comparison with infusion experiment

Table 2 :

 2 Processing data

Table 3 :

 3 Material data

	K 1 [m 2	]	K 2 [m 2	]	µ[P a.s 1	]	✏
	1.50.10 9	7.75.10 10	0.067		0.705

Table 4 :

 4 Processing data

ACKNOWLEDGEMENT

This work was supported by DAHER-Socata and CRCC. The authors are grateful to B. Cosson (ENSM Douai) and M. Bordival (ICA-Albi) for their contributions.


Denis Frath

Quyen Van

Frédéric Nguyen

Pascal G P Lafolet

Jean-Christophe Martin

Lacroix

Electrografted monolayer based on a naphthalene diimide-ruthenium terpyridine complex dyad: efficient creation of large-area molecular junctions with high current densities

were recently reported to act as rectifiers even though the rectification ratio (RR) remains modest. High RR were obtained using liquid mercury as the top electrode in a porphyrin/fullerene bilayer device 13 or using self-assembled alkylferrocene monolayers, contacted by chemisorption to a bottom electrode and physisorption to a top electrode. 14,15 These devices are among the best molecular rectifiers to date but their fabrication techniques suffer from the intrinsic stability limit of the SAM-based layer and the eGaIn contact. More recently, ferrocene, has been combined with a polychloro-triphenylmethyl unit to control the RR by changing the electronic nature of the component. 16 Different strategies have been used to construct the active layer, including Langmuir-Blodgett techniques, 17,18 layer-bylayer assembly 11,12 and electrochemical deposition. Diazonium electroreduction 19 is one of the best techniques, as the electrogenerated layers are compact, strongly adherent to the metal surface and pinhole-free. Consequently it is possible to evaporate the top electrode directly onto the molecular layer, which makes the creation of devices based on this approach very attractive.

20,21

Electrografting from organometallic diazonium derivatives has been reported [22][23][24] and was used to develop memories [START_REF] Cabaret | 14th IEEE International Conference on Nanotechnology[END_REF] or MJ with robust bipolar light emission and charge transport properties.

26

Naphthalene diimide (NDI) derivatives are commonly used electron acceptors in many organic electronic devices..

27

Recently they have been used to create various type of MJ that show moderate rectification properties. 28,29 Combining both Ru complexes and NDI results in donor-acceptor systems with multiple redox states. 28 Such systems could be of interest for memories or rectifiers. In order to investigate these aspects we have developed a new graftable dyad for the creation of large-area molecular junctions.

The molecular structure was designed with an NDI electronaccepting unit based on covalently bound to a ruthenium bisterpyridine complex (Scheme 1). The bridging phenyl ring is expected to decouple the two subunits. Steric hindrance between the hydrogen of this phenyl ring and the NDI oxygen atoms induces a large twist with respect to the NDI plane that should limit delocalization. 31 A hexyl chain was introduced to enhance solubility. The structure also features a phenylamine terminal group in order to enable in-situ preparation of a diazonium cation and subsequent grafting onto various substrates. The first synthetic step consists of the preparation of the alkylated precursor 1 by condensation of 1,4,5,8naphthalene-tetracarboxylic dianhydride and hexylamine. A previously reported procedure was adapted to introduce (4aminophenyl)-2,2':6',2"-terpyridine (tpyPhNH 2 ) and to generate NDI 2. 30 The resulting ligand was then complexed with the (4-aminophenyl)-2,2':6',2"-terpyridine ruthenium trichloride (Ru[tpyPhNH 2 ]Cl 3 ) precursor using procedures already reported. 32 Silica gel chromatography followed by counter-ion exchange was used to afford pure dyad 3 (NMR and HRMS are provided in Supporting Information files).

The electrochemical behavior of a solution of 3 (0.5 mM in MeCN + 0.1M TBAPF 6 ) was studied by cyclic voltammetry (CV) on a glassy carbon (GC) electrode (Figure 1a and Table S1 in SI)

Scheme 1. Synthetic pathway towards dyad 3:

(i) C6H13NH2, DMF, reflux, 40% (ii) tpyPhNH2, isoquinoline cat., m-cresol, 200°C, 74% (iii) Ru[tpyPhNH2]Cl3, Et3N, EtOH, Reflux, 18%
The curve displays two waves in the positive potential range, corresponding to the irreversible aniline oxidation at E pa = 1.08 V/SCE and the quasi-reversible oxidation of the Ru metal center at E 1/2 = 1.32 V/SCE. In the negative region, there are four partially reversible one-electron waves at E 1/2 = -0.52 V/SCE and E 1/2 = -0.95 V/SCE, corresponding to the reduction of the NDI fragment, and at E 1/2 = -1.30 V/SCE and E 1/2 = -1.49 V/SCE, corresponding to the reduction of the terpyridine ligands. Compared to the reference compound Ru(tpy) 2 PhNH 2 , there was no significant effect of the electron-withdrawing NDI fragment on oxidation of the ruthenium complex. Moreover, the complex only slightly affects the reduction waves of NDI (Table S1 and Figure S1 andS2). These results suggest that conjugation through the central phenyl group connecting the donor and acceptor parts of dyad 3 is negligible.

The diazonium cation derived from dyad 3 was generated insitu by addition of 15 eq. of tert-butyl nitrite (experimental details in Supporting Information). Ten successive cycles of the electrochemical reduction between 0 and -0.5 V/SCE were run on the resulting solution (Figure S3). In the first reduction cycle the irreversible peak at ca. -0.25 V/SCE corresponds to diazonium reduction. In the following cycles, the current intensity decreases, which suggests that dyad 3 is being grafted onto the surface. The electrochemical behavior of the modified electrode was studied by CV (Figure 1b and Table S1).

Reversible oxidation of the ruthenium occurs at E 1/2 = 1.32 V/SCE. The region of negative potential displays four reversible waves corresponding to the reduction of the NDI fragment at E 1/2 = -0.50 V/SCE and E 1/2 = -0.93 V/SCE, and for the reduction of the terpyridine ligands at E 1/2 = -1.24 V/SCE and E 1/2 = -1.53V/SCE. The additional peaks at E pa = 1.01 V/SCE and E pc = -0.72 V/SCE could be explained by a charge trapping phenomenon. 33 Similar behavior of electrografted layers of Rucomplexes has been reported. 22 After a few scans in oxidation or in reduction, these peaks decrease and the film remains stable upon further cycling (Figure S4). Integration of the peaks gives an apparent concentration of ca. 4.8x10 -10 mol.cm -2

. Taking bulkiness into account, it is consistent with values below 20x10 -10 mol.cm -2 usually observed for close-packed monolayers of aromatic core such as nitrobenzene.

34 To further characterize the dyad 3 layer obtained by electrografting, a gold electrode was modified by running ten successive cycles between 0 and -0.5 V/SCE. X-Ray photoelectron spectroscopy (XPS) was used to characterize the film (Figure 2a and Table S2). The attenuation of the Au4f signal of the modified substrate compared to a bare gold electrode indicated that the film is 2.7 nm-thick. This value is consistent with the molecular length obtained by molecular modeling (2.9 nm) and excludes the formation of multilayer. A characteristic peak at 281.4 eV can be assigned to Ru(II) 3d 5/2 . The corresponding Ru(II) 3d 3/2 peak at 285.4 eV is masked by the C1s signal (Figure 2b). On top of the main C1s peak at 285.1 eV, several other signals can be attributed to terpyridine ligand CN bonds at 286.3 eV, NDI imide CO bonds at 288.7 eV and a π-π* satellite due to aromatic cores at 292.4 eV. The single narrow N1s peak at 400.1 eV in the high-resolution spectra rules out the possibility of azo-bridged grafting (Figure 2c).

23

Experimental atomic ratios for C, N and Ru reproduce fairly well the expected composition of the layer (Figure 2d). The P/Ru ratio is approximately half the theoretical value. This can be rationalized by the fact that the first reduction of the grafted NDI occurs at E 1/2 = -0.5V/SCE. Since this potential is reached during grafting, the radical-anion of the NDI is formed and compensates one of the two charges on the ruthenium complex. It is reasonable to think that a PF 6 counter-ion is consequently removed from the film. Gold microelectrodes (ca. 20 μm x 2 cm and 50 nm-thick on Si/SiO 2 ) were functionalized with dyad 3 using the same electrografting conditions described above. The presence of the molecular layer was confirmed by AFM (Figure 2e). Surface topography reveals a roughness of ca. 0.8 nm, of the same order of magnitude as that observed for bare electrodes (ca. 0.6 nm). This shows that the molecular layer is quite densely packed and homogeneous. The thickness of the molecular layer is measured using the cross-sections of the modified and unmodified gold steps. A histogram of the height data is fitted by two separate Gaussian functions, with the height expressed as the difference between the centers of the two functions and the uncertainty given by the quadrature addition of the two bestfit σ values (Figure 2f). It gives an average value of 2.9 nm, close to the length of the molecule and confirming that the film is a monolayer (Figure 2g). Next, the modified microelectrodes were used to create MJs using photo-lithography step (see SI files). The top contact electrode was obtained by direct evaporation of Ti (2 nm) and Au (30 nm) onto the 2.9 nm-thick organometallic layer. Largearea molecular junctions were fabricated; current density vs. applied voltage (J-V) curves are presented in Figure 3. A high yield in operating junctions is reached (80%, 24 out of 30 junctions fabricated) showing that the monolayer is robust enough to enable MJ creation despite the fact that the top electrode is directly evaporated. This yield can be compared with that found when SAM-based monolayers are used: only 1.2% of the devices "work". 35 Another interesting feature is that very high current densities were reached: 9 A/cm 2 was observed when a bias of -2.5 V was applied. This has to be compared with current densities below 1mA/cm 2 in MJ incorporating alkane chains bearing redox moities, which are limited by the large attenuation factor of the chains.

15,16

Moreover, J-V curves feature an asymmetric profile with a higher current density when the bottom electrode is negatively polarized at a moderate potential of -2.5 V. In other words, the preferred electron flow is from the donor to the acceptor moities, i.e. anti-Aviram rectification. Similar behavior was observed with several systems 3,4 while opposite preferred current flows have been reported for bilayer MJ.

29

The average rectification ratio observed at -2.5V for the 24 operating junctions is 300 (see SI file for an histogram) and can be compared to recent results ranging from 350 to 2800 at 2.5V 6 or close to 900 in few junctions 36 with EGaIn top contact.

Note also that a value 10 5 was recently obtained but associated with low current density due to alkane moities. The exact mechanism behind the direction of the observed rectification is beyond the scope of this study and will be the subject of a more detailed article. Let us just recall that rectification behavior can be reversed when the effective tunneling distance between the electrode is changed by using redox moities incorporated at various distances from the electrodes, 38 or when the relative strength of the coupling between the molecular layer and the top or bottom electrode is changed. 39 Further investigation is needed to understand the rectifying properties of the present MJ, and we anticipate that they may be tuned and controled by changing the initial redox state of the molecular junction and the amount of counter-ion incorporated in the device. To conclude, a donor-acceptor dyad has been designed for the creation of large-area molecular junctions. This system can be conveniently grafted by diazonium electroreduction to form well controlled monolayers on various substrates. Films on glassy carbon characterized by cyclic voltammetry showed five stable redox states. Modified gold electrodes were obtained and fully characterized by XPS and AFM. The robustness of dyad 3 monolayer on gold substrates makes it possible to create large-area molecular junctions in high yield despite direct metal deposition on the monolayer. Thanks to the very low thickness of the device and the molecular structure of dyad 3, the current densities through these MJ are several orders of magnitude higher than those observed for other monolayer-based devices. Moreover, strong rectification behavior was observed with prefered electron transport from the donor to the acceptor moities. Detailed mechanisms of this asymmetric charge transport are being investigated in order to better understand this phenomenon and to optimize the rectification ratio of these very promising monolayer-based electronic devices.

Figure 1 .

 1 Figure 1. Cyclic voltammetry of dyad 3 (a) in MeCN (0.5 mM), TBAPF6 0.1M on GC at 0.1 V/s. (b) grafted layer on GC in MeCN, TBAPF6 0.1M at 0.1 V/s

Figure 2 .

 2 Figure 2. Characterization of Dyad 3 functionalized gold electrodes: (a) XPS survey (b) C1s-Ru3d high-resolution XPS spectra (c) N1s high-resolution XPS spectra (d) atomic ratio of the layer determined by XPS (e) AFM topography of the surface 5x5 μm 2 (top) and cross section (bottom) of the modified nanostructured electrode corresponding to the thickness measurement, (f) AFM statistical analysis of the thickness, fitted with gaussian (g) schematic representation of the monolayer.

37

 37

Figure 3 .

 3 Figure 3. a) schematic illustration of the junction. b) JV and semilogarithmic JV curve for Au/dyad 3 monolayer/Ti/Au MJ with standard and log-standard deviations (7 junctions have been tested).

We thank the the Centre National de la Recherche Scientifique, and Agence Nationale de la Recherche (ANR-15-CE09 0001-01) for their financial support. We also warmly thank B. Lemarque and P. Decorse for XPS analyses and discussions.

Notes and references


A Lefebvre-Lepot
email: aline.lefebvre@cmap.polytechnique.fr

Numerical simulation of suspensions: lubrication correction, including fluid correction

Keywords: Suspensions, Numerical simulations, Stokes, Lubrication, Simulation numérique de suspensions : prise en compte des forces de lubrification

, where the singular field is tabulated. Here, we propose a new method, based on an asymptotic expansion of the singular field. This method allows to catch the effect of the lubrication on the whole velocity and pressure flows and can deal with any forms of particles. We focus on a toy problem in two dimensions to present the method. Numerical results are given, using a finite element discretization.

Résumé. Lors de simulations numériques d'écoulements de particules dans un fluide de Stokes, se pose inévitablement la question de la gestion des interactions entre particules proches. En effet, quand deux solides se rapprochent, les champs de vitesse et de pression sont singuliers et il devient difficile de les approcher numériquement. Or, la bonne prise en compte de ces interactions est primordiale, tant d'un point de vue physique que numérique. Par ailleurs, les expérimentateurs ont besoin de résultats de plus en plus précis, prenant en particulier en compte l'effet de ces interactions sur la totalité du champ fluide. La méthode que nous proposons ici consiste à résoudre le problème fluide/particules en le décomposant en deux sous problèmes : un problème singulier (quand les distances inter-particulaires tendent vers zéro) et un problème régulier. La partie singulière du champs est supposée connue et le problème initial se ramène alors à résoudre un problème régulier. Une première approche avait été proposée dans [START_REF] Lefebvre-Lepot | An accurate method to include lubrication forces in numerical simulations of dense Stokesian suspensions[END_REF], dans laquelle le champs singulier était tabulé. Nous proposons ici une nouvelle méthode, basée sur un développement asymptotique du champs singulier. Cette méthode permet de prendre en compte les effets de la lubrification sur les champs de vitesse et de pression dans tout le domaine et permettent la prise en compte de particules de formes quelconques. Afin de présenter la méthode, on se concentre ici sur un problème jouet en deux dimensions d'espace. Nous présentons des résultats numériques basés sur une discrétisation éléments finis.

Mots clés. Suspensions, simulations numériques, Stokes, lubrification

Introduction

Dense suspensions are present in various domains coming from industry (nuclear waste reprocessing, concretes, reinforced plastics) to natural environment (silting up), biology (blood tests) or even sanitary concerns (wastewater treatment). This large area of applications gave rise to a great amount of research. However, the flow properties of these systems are still partly misunderstood. The applications we have in mind such as studying migration or understanding the interphase stress contribution to the total fluid stress involve dense systems, close to the maximal packing and for which the grains are most of the time in contact (or quasi-contact). It consists in a challenging problem, for which no result is available (neither experimental, nor theoretical, nor numerical).

From a computational point of view, it leads to suspension simulations with density between 30% and 50%, containing up to about one million of suspended particles. Experimental results show that the key points to make a breakthrough in the domain is to compute the velocity and pressure fields in the whole fluid domain for a dense suspension, to model carefully the multi-body lubrication and contact interactions and their feedback on the flow. An evaluation using numerical simulations is not possible with classical methods since they do not calculate accurately the fluid pressure and the resulting flow when the particles are close to each other. As a consequence, there is a great need of designing new mathematical models, stable and efficient numerical methods, together with HPC programming.

The objectives are the following : calculate accurately the whole pressure and velocity fields even for coarse meshes and small distances, deal with any form of particles (poly-disperse suspensions of spheres, fibres, anisotropic particles...), take advantage of the existing fast solvers for Stokes equations and obtain a parallelizable method. For the sake of clarity, we consider in this manuscript a toy problem involving 2 particles in 2 dimensions. Each of the methods presented can be easily extended to multi-particle suspensions in 3 dimensions with general relative displacements.

In section 2, we present the problem and the lubrication singularity. Then, we explain how the relative movement of close particles can be decomposed in the sum of two movements : a regular one and a singular one. We explain how this decomposition has been used since the Stokesian Dynamics was proposed in 1988 [START_REF] Durlofsky | Dynamic simulation of hydrodynamically interacting particles[END_REF] to compute precisely the total force exerted on the particles. Then, we present a recent work [START_REF] Lefebvre-Lepot | An accurate method to include lubrication forces in numerical simulations of dense Stokesian suspensions[END_REF] where the decomposition was used to compute precisely the whole fluid and pressure fields. In both of these methods, tabulations are used to compute accurately the singular forces or fields, which prevents the use of the methods for non mono-disperse suspensions.

In section 3, we present a new method to decompose the velocity and pressure fields using an explicit expansion of the singular fields, which is available for any form of particles. We explain how the decomposition can be implemented using finite elements and finally, we present numerical results and compare the precision of the new method with a classical direct method.

2. Numerical simulation of suspensions : a singular problem.

2.1.

A toy problem. We consider (see Fig. 1) Ω = [0, 1] × [0, 1] ⊂ R 2 . We denote by (B i ) i=1,2 two rigid spheres in Ω. B = ∪ i B i is the whole rigid domain. We suppose that Ω \ B is filled with a Newtonian fluid governed by the Stokes equations. The two spherical particles are aligned along the horizontal direction : if X i ∈ R 2 denotes the center of particle i, we have X 1 = (0.5 -d/2 -r 1 , 0.5) and X 2 = (0.5 + d/2 + r 2 , 0.5) where r i is the radius of particle i and d stands for the distance between the two particles.

Ω \ B X 1 × B 1 X 2 × B 2 r 1 r 2 d V 1 V 2 Figure 1. Notations.
We consider homogeneous Dirichlet boundary conditions on ∂Ω. We suppose that the velocities V 1 , V 2 of the particles are given and are oriented along the horizontal direction. The problem we consider consists in computing the corresponding pressure field, velocity field and forces :

(1)

                               Find (u, p) ∈ H 1 (Ω \ B) × L 2 (Ω \ B) such that -u + ∇p = 0 in Ω \ B ∇ • u = 0 in Ω \ B u = 0 on ∂Ω u = V i on ∂B i , i = 1, 2
and compute the corresponding forces exerted on the particles :

F i = ∂Bi σ(u, p) • n, i = 1, 2.
where σ is the stress tensor and n the outward normal to the boundary of the particles.

2.2.

A singular problem due to lubrication. It is well known that the main difficulty in designing numerical methods for approximating (1) arises when the configuration involves close particles. Indeed, when the two particles are almost in contact, the possible discrepancy between the velocities on both sides of the thin gap between them leads to large variations of the velocity fields. This phenomena combines with the incompressibility constraint to generate lubrication forces : large hydrodynamic forces with densities located in small areas. More precisely, let us consider the toy problem for two spheres of radius 1 with opposite velocities V 1,2 = ±e x (see Figure 2, extracted from [START_REF] Lefebvre-Lepot | An accurate method to include lubrication forces in numerical simulations of dense Stokesian suspensions[END_REF] where the force distribution was numerically computed using a decomposition in spherical harmonics). The magnitude of the force density exerted on ∂B 1 is displayed on the figure for different values of d. We can observe that the force densities concentrate on spherical caps with radii of order of √ d. For such large localized densities, if one wants to compute the hydrodynamic forces exerted on the particles, a large number of degrees of freedom is required in order to capture the relevant phenomenon. Therefore, simulations based on volume discretization (such as finite element or finite volume methods) will require very fine meshes to take these forces into account. For spectral discretizations (such as multipole methods) this corresponds to a truncation order L satisfying L 1/ √ d. Such a precision cannot be achieved in simulations of 3-dimensional dense suspension and a numerical strategy has to be implemented to take the lubrication effects into account.

2.3.

A singular/regular decomposition of the movement. All the numerical methods to include lubrication into simulations are based on two key ingredients. First, the movement of two close particules can be decomposed in the sum of two movements : a regular one plus a singular one. Then, once the decomposition is achieved, the behavior of the singular part of the solution can be investigated in a theoretical framework and the regular movement can be easily approximated using numerical solvers.

More precisely, using the linearity of Stokes equations, problem (1) can be decomposed in two similar problems with modified boundary conditions on the particles (see Fig. 3). In the first regular problem, the two particles undergo a global movement of translation. The solution of this problem do not contain singularity in the sense that the norm of the solution is bounded when the distance goes to zero. In the second singular problem, they move one towards the other with the same velocity. The key idea to estimate the singular effects due to the singular boundary conditions is to note that the asymptotic behavior of the solution is well known, provided we consider that the two particles are embedded in an unbounded fluid.

Such a strategy has been proposed by Durlofsky, Brady and Bossis in 1988 to take lubrication into effects in Stokesian Dynamics simulations [START_REF] Durlofsky | Dynamic simulation of hydrodynamically interacting particles[END_REF] and has been widely used since then. Recently, their method has also been used to include lubrication effects in direct solvers [START_REF] Gallier | A fictitious domain approach for the simulation of dense suspensions[END_REF]. Their founding remark is the following : supposing that the particles are embedded in an infinite fluid and have the same radius, the solution to the problem with opposite velocities V sing 1,2 = ±e x only depends on the distance between the particles. Then, in that case, the forces exerted on the two particles can be expressed as a (vector-valued) function of d : d → F sing i,R 2 (d). This function can be accurately approximated for any distance using its well-know asymptotic expansion when d goes to zero [START_REF] Cox | The motion of suspended particles almost in contact[END_REF], together with an accurate tabulation for greater distances. As a consequence, since the asymptotic behavior when the distance goes to zero of the initial problem in Ω is the same as the one in R 2 , the forces exerted on the particles can be computed as follows :

B 1 B 2 d V 1 = V 1 ex V 2 = V 2 ex = V 1 +V 2 2 × regular movement B 1 B 2 d ex ex + V 1 -V 2 2 × singular movement B 1 B 2 d ex -ex
F i = F sing i,R 2 (d) + F reg i , i = 1, 2
where the remaining part of the force F reg i is regular (i.e. is bounded when the distance goes to zero) and can be computed using any numerical solver with a small number of degrees of freedom.

This singular/regular decomposition provides an accurate approximation of the forces exerted on the particles. However, using this method, there is no way to obtain velocity and pressure fields including lubrication effects. Here, the fluid solver is used to compute the regular part of the forces and no correction is achieved on the flow to provide an approximation of the fields (u, p) solution to [START_REF] Cox | The motion of suspended particles almost in contact[END_REF]. Moreover, a precise tabulation of the forces is possible, provided the corresponding function only depends on the distance between the particles. In the case it depends on two parameters or more, this step is too demanding in CPU time, even if achieved offline. This prevents the use of this correction method in simulations of suspensions made of poly-disperse spherical particles (the singular part of the forces depends both on the distance and the ratio between the radius of the particles) or suspensions of non spherical particles (the singular part of the forces depends on the distance, the form and the relative position of the particles).

A singular/regular decomposition of the flow

3.1. Description of the method. We proposed in [START_REF] Lefebvre-Lepot | An accurate method to include lubrication forces in numerical simulations of dense Stokesian suspensions[END_REF] a new method based on a singular/regular decomposition of the whole flow. We write [START_REF] Durlofsky | Dynamic simulation of hydrodynamically interacting particles[END_REF] u = u reg + u sing , p = p reg + p sing .

The singular field is supposed to catch the singularity of the solution when the distance goes to zero in such a way that the remaining regular field is bounded independently of the distance in H 1 (Ω \ B) × L 2 (Ω \ B). Let us for the moment suppose that such a singular field can be computed (this point will be discussed later). Then, to solve the initial problem (1) for (u, p), it is enough to compute (u reg , p reg ), solution to the following regular problem :

(3)

                                 Find (u reg , p reg ) ∈ H 1 (Ω \ B) × L 2 (Ω \ B) such that -u reg + ∇p reg = u sing -∇p sing in Ω \ B ∇ • u reg = -∇ • u sing in Ω \ B u reg = -u sing on ∂Ω u reg = V i -u sing on ∂B i , i = 1, 2
and compute the corresponding forces exerted on the particles :

Freg i = ∂Bi σ(u reg , p reg ) • n, i = 1, 2.
3.2. Choice for the singular field. If (u reg , p reg ) is computed as a solution to problem (3), the field (u, p) given by ( 2) is solution the initial problem (1), whatever the singular field is. The remaining question is to manage to design (u sing , p sing ) such that the regular problem (3) can be numerically approximated with a number of degrees of freedom not depending on the distance between the particles. This can be achieved using, as for the Stokesian Dynamics correction, the singular part of the movement described in Fig. 3. Note that the decomposition has to be achieved as soon as the relative movement of the two particules is non zero i.e. as soon as V 1 is different from V 2 . In [START_REF] Lefebvre-Lepot | An accurate method to include lubrication forces in numerical simulations of dense Stokesian suspensions[END_REF], the singular field (u sing , p sing ) is defined as

(u sing , p sing ) = V 2 -V 1 2 (u sing R 2 , p sing R 2 ) where (u sing R 2 , p sing R 2 )
is the solution to the fluid-particle problem in R 2 with the boundary conditions V sing 1,2 = ±e x . Again, for particles with the same radius, this singular field only depends on the distance between the particles and in [START_REF] Lefebvre-Lepot | An accurate method to include lubrication forces in numerical simulations of dense Stokesian suspensions[END_REF] we propose to tabulate its coefficients in a well-chosen basis. This leads to a precise method to compute the solution to problem [START_REF] Cox | The motion of suspended particles almost in contact[END_REF]. Here, the lubrication is taken into account in the computed flow (u, p). However, since the singularity is computed using tabulations, the method is limited to mono-disperse suspensions of spherical particles. Moreover, the singular field is not truncated and therefore, the boundary conditions for the regular problem has to be modified in the whole boundary of Ω \ B, which will induce heavy computations when extending the method to dense suspensions.

Here, we propose a new choice for the singular field, based on an asymptotic expansion of the flow when the distance goes to zero. The computations to obtain such an expansion are similar to the ones achieved in [START_REF] Cox | The motion of suspended particles almost in contact[END_REF] to obtain the asymptotic behavior of the total lubrication force. More recently, in [START_REF] Hillairet | Justification of lubrication approximation : an application to fluid/solid interactions[END_REF], the authors provide explicit formulas for the expansion and prove that the remaining part of the flow is regular in the sense that its norm is bounded in

H 1 0 (Ω \ B) × L 2 (Ω \ B)
independently of the distance between the particles. We propose to use

(4) (u sing , p sing ) = V 2 -V 1 2 (u sing exp , p sing exp )
where (u sing exp , p sing exp ) is a truncation of the explicit asymptotic expansion for the fluidparticle problem in R 2 with the boundary conditions V sing 1,2 = ±e x . An example of such a field is given in Fig 4.

Then, solving the corresponding regular problem (3) and computing the solution as the sum of the regular and singular fields(2), the lubrication is taken into account in the computed flow (u, p). The new singular field is explicit and only depends (explicitly) on the local geometry of the particles near the (almost) contact point between them. Then, the method can be implemented for poly-disperse suspensions of spherical particles and for various forms of particles. Moreover, the singular field is localized in space, which avoids too many computations when modifying the boundary conditions to solve the regular problem (3).

Numerical implementation.

We propose to solve (3) using a finite element method. Let V h ⊂ H 1 (Ω \ B) and M h ⊂ L 2 (Ω \ B) be a stable pair of finite element spaces for the Stokes problem. The singular field is given by the truncated asymptotic expansion (4). We want to compute the approximation of the regular field (u reg h , p reg h ) solution to

(5)

                 Find (u reg h , p reg h ) ∈ V h × M h such that -u reg h + ∇p reg h = u sing -∇p sing in Ω \ B ∇ • u reg h = -∇ • u sing in Ω \ B u reg h = 0 on ∂Ω u reg h = V i -u sing on ∂B i , i = 1, 2,
where the right-hand side is computed as precisely as possible. Then, from (2) ,the discrete field approximating (u, p) solution to problem (1) is

u h = u sing + u reg h ∈ u sing + V h p h = p sing + p reg h ∈ p sing + M h
Non-homogeneous Stokes problems as (3) or ( 5) have a unique solution if and only if a compatibility condition is satisfied by the right-hand side of the problem, expressing the fact that one may have

Ω\ B ∇ • u reg = ∂Ω\ B u reg • n and Ω\ B ∇ • u reg h = ∂Ω\ B u reg h • n.
Concerning the continuous problem (3), the compatibility condition is verified since

Ω\ B ∇ • u sing = ∂Ω\ B u sing • n.
It is not true anymore for the discrete problem [START_REF] Hillairet | Justification of lubrication approximation : an application to fluid/solid interactions[END_REF]. Indeed, inside Ω, the divergence equation gives

Ω\ B ∇ • u reg h = Ω\ B ∇ • u sing .
Concerning the boundary we have

u reg h | ∂Ω\ B = V i -Proj u sing on ∂B i
where Proj is the projection onto the space made of the traces on ∂Ω \ B of the functions of V h . Then,

Ω\ B ∇ • u reg h = Ω\ B ∇ • u sing = ∂Ω\ B Proj u sing • n = ∂Ω\ B u reg h • n.
As a consequence, problem (5) may not have solutions.

To circumvent this problem, we chose ε > 0 and we introduce a modified discrete regular problem :

(6)                    Find (u reg h,ε , p reg h,ε ) ∈ V h × M h such that -u reg h,ε + ∇p reg h,ε = u sing -∇p sing in Ω \ B ∇ • u reg h,ε + εp reg h,ε = -∇ • u sing in Ω \ B u reg h,ε = 0 on ∂Ω u reg h,ε = V i -u sing on ∂B i , i = 1, 2.
This problem has a unique solution. Testing the divergence equation with v h ≡ 1 ∈ V h we see that

ε Ω\ B p = - Ω\ B ∇ • u sing + ∂Ω\ B Proj u sing • n = 0.
The classical unknown constant in the pressure field will be chosen by the solver so as to compensate the fact that the compatibility condition is not satisfied from the numerical point of view.

The constant ε has to be chosen sufficiently small in order that the error introduced in the divergence equation is smaller than the one due to the finite element discretization. On the other hand, the finite element discretization being fixed, letting ε go to zero in the numerical tests shall induce numerical errors.

Numerical results

The numerical tests have been implemented with freefem++ [3]. We chose V h = P 1 b and M h = P 1 . The configuration is the one given in Fig. 1 with r 1 = 0.07, We plot in Fig. 5 the numerical solutions obtained for d = r 2 /10. The solution is given by the singular/regular decomposition with ε = 10 -4 and using a quadrature formula of order 10 to compute the right-hand side of the regular problem. The direct solver Crout is used to solve the linear system. The reference is computed using the direct finite element method for (1), the same direct linear system solver, and a very fine mesh (until convergence). The constant in the pressure field is computed such that p = 0. We can observe that, even for a very coarse mesh, the solution is well approximated.

r 2 = 0.1, V 1 = e x , V 2 = -e x .
To be more precise, we plot on Fig. 6 the relative errors in norm L 2 (Ω \ B) for the pressure and H 1 (Ω \ B) for the velocity. The convergence of the direct method and the decomposition method are compared for different distances. We can observe that the classical order of convergence for the P 1 b ×P 1 approximation of Stokes equations is recovered. The decomposition method is equivalent to the direct method for d = r 2 /2 but becomes more precise than the direct method when the distance between the particles decreases. These numerical tests also confirm that the precision is less dependent on the distance when the decomposition method is used.

Conclusion

We presented in this manuscript a new method to include lubrication in numerical simulations of particles. The lubrication effects are transmitted to the whole velocity and pressure flows. The singularity is taken into account by solving a numerical regular problem which is a Stokes problem with a modified right-hand side. This allows the use of existing fast solvers to compute its solution. It can deal with poly-disperse suspensions and with any forms of particles since the explicit asymptotic expansion only depends on the local geometry of the particles.

This method can straightforwardly be extended to the three dimensional case with any number of particles. In this manuscript we tested the singularity induced by normal relative translation of the two particles. The general relative movement of two particles can be decomposed as a regular movement (a global rigid movement) plus a singular movement including normal and tangential relative translations as well as relative rotations. All these singular relative movements can be taken into account in the simulations since the asymptotic expansion of the corresponding singular fields are available.

We believe this method can be useful for intensive simulations of dense threedimensional suspensions. Indeed, in that case, the lubrication effects are predominant but only coarse meshes can be used. Moreover, the size of the problem imposes to use iterative solvers. Using the decomposition method, we end with a regular numerical problem and then, less iterations of the solvers shall be needed to obtain a given precision. Moreover, the computation of the right-hand side can be parallelized since the contributions of each pairs of particles are independent.

Figure 2 .

 2 Figure 2. Pair of close spheres of radius 1. Notations (left). Magnitude (middle) and normalized magnitude (right) of the force density near the (almost) contact point on ∂B 1 for various distances : d = 2.10 -2 (• • •), d = 4.10 -2 ( ), d = 8.10 -2 (+ + +) and d = 16.10 -2 (• • •).

Figure 3 .

 3 Figure 3. Decomposition of the movement into the sum of two movements : a regular movement and a singular movement.

Figure 4 .

 4 Figure 4. Truncated asymptotic singular field. r 1 = 0.07, r 2 = 0.1, d = r 2 /10. p sing exp (top) and u sing exp , horizontal component (bottom left) and vertical component (bottom right).

  Three distances have been tested : d = r 2 /2, d = r 2 /5 and d = r 2 /10.

Figure 5 .

 5 Figure 5. Numerical solution. Reference (left) and numerical solution (right). Zoom of the mesh in the gap between the particles (top), pressure (line 2), horizontal velocity (line 3) and vertical velocity (line 4)

Figure 6 .

 6 Figure 6. Relative errors versus h for the direct method (• • • ) and the decomposition method (-). The tests are run for different values of the distance : d = r 2 /2 (+), d = r 2 /5 (•) and d = r 2 /10 ( ). Slope 2 is given on the top plot (black line). Slope 1 is given on the middle and bottom plot (black line).

Aknowledgements

The author would like to thank warmly Pr. François Alouges (CMAP-Ecole Polytechnique) for helpful discussions about this work.

This work was partially supported by Investissements d'Avenir of LabEx PALM (ANR-10-LABX-0039-PALM) and Labex LMH (ANR-11-IDEX-003-02 and ANR-11-LabX-0056-LMH).

Références


Ronan-Alexandre Cherrueau
email: ronan-alexandre.cherrueau@inria.fr

Matthieu Simonin
email: matthieu.simonin@inria.fr

Alexandre Van Kempen
email: alexandre.van-kempen@inria.fr

EnosStack: A LAMP-like stack for the experimenter

Keywords: Repeatability, Reproducibility, Application deployment, Performance, Grid5000, Chameleon

Reproducibility and repeatability dramatically increase the value of scientific experiments, but remain two challenging goals for the experimenters. Similar to the LAMP stack that considerably eased the web developers life, in this paper, we advocate the need of an analogous software stack to help the experimenters making reproducible research. We propose the EnosStack, an open source software stack especially designed for reproducible scientific experiments. EnosStack enables to easily describe experimental workflows meant to be re-used, while abstracting the underlying infrastructure running them. Being able to switch experiments from a local to a real testbed deployment greatly lower code development and validation time. We describe the abstractions that have driven its design, before presenting a real experiment we deployed on Grid'5000 to illustrate its usefulness. We also provide all the experiment code, data and results to the community.

I. INTRODUCTION

In early 2000, the LAMP stack was considered as the Graal for web developers. The stack was composed of Linux as Operating System, Apache as Web server, MySQL as database backend and PHP as scripting language. The stack was packaged all together in the major linux distributions increasing its adoption. It was opinionated but decoupled clearly the different parts needed to build a web application. The application logic later moved to dedicated frameworks (e.g Django, Ruby On Rails) which further abstract the database backend through the use of ORMs (Object-Relationnal Mapping) As a consequence the efficiency of the developer increased due to easy switching between a local development mode and production deployment.

In this paper, we argue that experimenters also require their LAMP stack to assist them in making reproducible research. Indeed, while reproducible research is a key component of the scientific method, it is far from being the norm, especially in the context of computational science [?], [?]. Even if experimenters are committed to automate their experiment code, some tasks remain tedious [?]. First, maintaining a code to keep up with the development of the upstream software code can be time consuming (e.g OpenStack is released every six months). Second, sharing the code with others requires to follow coding and packaging gold standards [?]. The latter helps experimenting in a repeatable way, the former addresses the problem of revisiting results in different time frames.

Leveraging concepts introduced in frameworks such as OMF [?], we propose the EnosStack, an open source software stack, as a solution to such challenges by following the model of the LAMP stack. Parallely the authors of Popper [?] paved the way towards a set of common conventions that would ensure an experiment to be easily re-executed and validated. EnosStack is a pragmatic framework in which such conventions can be implemented.

The software stack is composed of a group of open source software that are typically installed together to ease the orchestration of reproducible experiments. It includes Python, Ansible, Docker and the EnosLib, a library we specially developed for this work. Large distributed applications were the initial targets of the EnosStack, but it appears flexible enough to target all kind of applications.

In order to introduce the concepts of EnosStack while easing the description of its implementation, we present a red-thread example that we will rely on throughout the paper to illustrate more concretely various concepts. This example consists in evaluating the scalability of a message-orientedmiddleware (MOM, e.g. RabbitMQ1 ) stressed by multiple agents communicating through it. To perform this experiment on a real testbed (e.g Chameleon [?], Grid'5000 [?]), an experimenter would need to configure all the agents as well as the MOM and then deploy some code that benchmarks it for a given number of agents. This process would have to be repeated with an increasing number of agents until reaching the limit of the MOM. During these experiments, various metrics would have to be monitored and backed up before the resulting data analysis. As shown in the next sections, a such seemingly simple experiment leads to numerous design challenges when having to be deployed on real testbeds.

More specifically, the contributions of this paper are as follows:

• We identify key concepts and good practices to deploy experimental workflows. • We introduce the EnosStack, a novel LAMP-like stack for the experimenter. • We present a use case experiment that evaluates the scalability of a RabbitMQ bus, deployed on Grid'5000 thanks to the EnosStack. The remainder of the paper is structured as follows: section II discusses high level requirements, section III introduce concrete technologies for EnosStack. The section IV presents the most important concepts of the EnosStack that help in Acquire some resources Deploy some code

Run the experiments

Retrieve the results

Analyze the results

Destroy

Fig. 1. Experimental workflow fulfilling the previous requirements. Finally, section V will present the evaluation of the real deployment of our red thread example on Grid'5000 before concluding in Section VI.

II. HIGH-LEVEL REQUIREMENTS

We envision the process of creating an experiment in four steps: Designing, validating, scaling and sharing. Each step is in general challenging [?] thus a stack-for-experimenter must ease each one of them. In this section, we outline some highlevel requirements that emerge from the above steps.

(i) Designing the experiment: A typical experiment can be expressed as a workflow of tasks as in Fig. 1. This workflow starts by [START_REF]Docker images for the use case experiments[END_REF] getting some resources from a specific testbed, (2) deploying the application itself and all the third-party software needed (e.g instrumentation), (3) running the wanted workload (4) retrieving the results, (5) destroy everything before iterating on all the previous steps for other configurations, or parameter set. Once the required parameter space of the experiment has been explored [START_REF] Balouek | Adding virtualization capabilities to the Grid'5000 testbed[END_REF] the results analysis can finally be performed. The enactment of the above workflow must be fully automated and repeatable. As a result, being able to easily describe the workflow and to execute it without sacrificing fast iterations on the code is a strong requirement for the EnosStack.

(ii) Validating the experiment: The typical experimenting process can be split in two phases: first, the development phase and then, the production phase. In the development phase an experimenter needs to iterate over the code of the experiment until reaching the desired level of automation and correctness. In the production phase, the experimenter will actually run the code of the experimentation on the targeted platform (e.g Grid'5000, Chameleon). Most probably, experimenters will need to go back and forth between the two phases while refining their experimentation code (bug fixing, scalability issue, etc.). Production platforms are often shared between users and rarely offer a friendly development environment. Thus, skipping the first phase leads an experimenter to waste precious resources on the production platform and be less efficient while iterating on the code. As a consequence, being able to transparently switch between environments is another strong requirement for the EnosStack.

(iii) Scaling the experiment: After validating the workflow and its correct execution, experimenters likely need to scale it. This is typically achieved by deploying more processes of the application under study. To run those extra processes one can choose (a) to claim more compute resources to the underlying infrastructure to run them or (b) to increase the density of the processes on the servers. In both cases, an experimenter need a simple way to specify it and an effortless way to enforce it.

(iv) Sharing the experiment: By sharing, we mean sharing the whole experiment or part of it. The former provides reusability. It is catalyzed, for instance, when an easy path to install (e.g package) is offered or a friendly user experience (e.g command line interface) is provided. The latter refers to the modularity of a software project. In any cases, following any high-level language guidelines of software development leads to these properties. As a consequence, an experimenter who writes experiment as code needs to be given the same development experience as in a regular software project.

III. THE ENOSSTACK LAYERS

The EnosStack's cornerstone is provided by the EnosLib, a library we specially developed for this work. Leveraging the EnosLib, we built the EnosStack by including well known open source software namely Python, Ansible and Docker. We describe hereafter these four components in more details.

1) Python: Python is a versatile language, that comes with one of the most mature package libraries. It is particularly adapted to scientific experiments thanks to powerful libraries tailored for data science such as numPy and sciPy [?]. Moreover, it neatly interfaces with Ansible thus making it the language of choice for the EnosStack.

2) Ansible: Ansible is an orchestration tool that automates software provisioning, application deployment and configuration management. We detail hereafter some of its concepts we will use throughout the paper. In Ansible, a task is an action that is remotely controlled on a machine. Tasks are launched on groups of nodes using a fork-join model. The sequence of tasks to run is described in a playbook. Ultimately, playbooks can be packaged alongside their required files and variables to form an Ansible role. Finally, an Inventory file allocate hosts into different groups to logically map physical nodes to Ansible roles.

Including Ansible in the EnosStack has been driven by the fact that, contrary to other similar tools, Ansible only requires an SSH connection as it follows a push model and is thus agentless. This avoids the need of any software installation on the targeted machines, making it very lightweight to setup. Additionally it's a robust software backed by a large community of users. Its role mechanism is particularly interesting as it provides a good level of modularity in the code produced or re-used from a third-party user.

3) Docker: Docker is an open-source project designed to ease the deployment and execution of applications by using containers. Containers are like packages that include all the libraries and other dependencies required by a given application. Being self-sufficient, containers offer an inherent portability, while considerably easing the lifecycle of containerized applications. These are properties particularly suitable in our context of reproducible experiments as explained in [?]. Moreover they provide a unified abstraction for the life-cycle management of programs. This eases the description of the deployment of multi-containers applications.

4) Enoslib: The Enoslib [?] is an open source library we developed that enables a developer to describe and run an experimentation workflow on different infrastructures. It abstracts the underlying infrastructure to a developer, such that only few changes are required to migrate an experiment from one testbed to another testbed. It also allows the experimenter to define tasks and combined them to form the desired experimental workflow in a way that allows incremental code development. Enoslib takes its roots in EnOS [?] but goes beyond as it gets agnostic to the application under study.

IV. ENOSSTACK ABSTRACTIONS

We now present the main abstractions the EnosStack rely on. In order to fulfill the requirements explained in the previous section, we carefully designed the EnosStack to abstract concepts such as resources, services or tasks as explained hereafter.

A. Resource Model

In the EnosStack, the concept of providers is the key that allows experimenter to switch from one execution context to another. EnosStack is shipped with three providers: Vagrant, Grid'5000, OpenStack. The former is mainly used in the development phase while the others will be preferably used in the production phase. Note that an implementation for the Chameleon platform is also possible and is inherited from the OpenStack provider. EnosStack models two types of resources as described in the next sections.

a) Machines: In EnosStack, machines are compute resources and can be grouped. A group is given a count and a set of roles. In the current implementation, machines in one group share the same hardware characteristics. A group is an abstract description and it is the provider responsability to concretize this description by : (1) claiming count machines from the underlying infrastructure (2) distributing those machines in the roles. The roles of the machines will be used at deployment time to differentiate the configuration to apply on each machine.

b) Networks: In EnosStack, networks provide connectivity between the different machines. Networks are assigned to machines, allowing to model non trivial topologies. Depending on the underlying provider used, the implementation may differ : bridge networks on Vagrant, vlans on Grid'5000 or private networks on OpenStack. Similar to a group of machines, a network can be given a set of roles. For example, m a c h i n e s :

f l a v o r : l a r g e c o u n t : 3 r o l e s :

mom, t e l e g r a f i n f l u x d b n e t w o r k s :

c o n t r o l n e t w o r k i n t e r n a l n e t w o r k f l a v o r : t i n y c o u n t : 10 r o l e s :

mom-a g e n t s t e l e g r a f n e t w o r k s :

- c) Example: Fig. 2 gives the declaration of the resources for the same experiment. Fig. 2(a) will start virtual machines and get two networks from the local VirtualBox hypervisor. Fig. 2(b) will deploy bare-metal servers on the Grid'5000 platform and get one isolated VLAN network (the second network is the default provided by Grid'5000). The description of machines is slightly different : size refers to a predefined amount of vCPU and Memory in the Vagrant case whereas cluster name is used in the Grid'5000 case. Also, network description for Grid'5000 allows to specify inhomogeneous network card configuration. More concretely here, specifying a count greater than one deploys several agents of the MOM in a clustered mode where the internal traffic (e.g replication traffic between RabbitMQ brokers) is isolated in a VLAN. Note also that the experimenter can give the two roles to the same network and traffic will not be isolated. This flexibility is allowed by the service abstraction presented in the next part.

B. Service model

During the deployment of an experiment, software will be installed and configured. Those software can be libraries but also running daemon processes that can optionally require a persistent storage to keep their states. In EnosStack, we call service an abstraction of one or several software that serves the same purpose regarding the application to deploy. A service has inputs (e.g configurations options for the encapsulated programs) and actions that transition the state of the service. This service abstraction is beneficial for both the extensibility of an experiment (e.g considering new scenarios) and re-usability (e.g using some building blocks in another experimental workflow).

Going back to the red thread example, a monitoring infrastructure is required to collect performance metrics. Since it's likely to be re-used in another context we model this infrastructure as a service. It is composed of monitoring agents (Telegraf2 ), metrics collector (InfluxDB3 ) and a graphical frontend (Grafana4 ). Note that the monitoring service can be seen as a composite service composed of three sub-services: telegraf, influxdb and grafana. To ease the management of their life-cycle, we chose to containerize all the corresponding components. We also consider the following actions and semantics: deploy installs and configures all the programs according to the inputs, backup fetches the persistent storage of the metrics collector (the docker volume attached to the collector container) to keep all the metrics generated. Finally, destroy removes all the corresponding containers as well as their associated volumes.

At deployment time, EnosStack allows to map a service to a role and thus to all the machines in that role. More concretely, according to the description given in Fig. 2, the telegraf subservice will be deployed on thirteen machines. At this point increasing the number of machines associated with the telegraf role will scale automatically the number of running monitoring agents. Additionally changing the monitoring agent software to an equivalent one is as simple as switching the monitoring agent service implementation. In conclusion, following the above service model allows to benefit from a certain degree of modularity and scalability.

C. Task model

A typical experiment follows a workflow of tasks as outlined in the section II. For instance, the experimenter of the redthread example first claims resources on a testbed. Then, she deploys the MOM and monitoring services. Next, she stresses the MOM with a given number of agents, and backups results of the experiment. Finally, she destroys the experiment to start a new one in a fresh environment. The new experiment will follow a nearly similar workflow of tasks (only the stress task will vary by increasing the number of agents) which emphasis the need for fast iterations on a task to adapt it to new experiments.

In the EnosStack, a task is a python function. By itself, a python function already provides mechanisms for task's adaptation. This is the case of function parameters that adapt the task when they are given a value. But, such mechanism are not sufficient to provide fast iteration on an experiment. To get fast iteration, the experimenter needs facilities to easily stop, then adapt, and finally restart the experiment at the moment of a specific task. Indeed, during the development phase, an experimenter often has to exit the experiment execution, e.g., to fix it, before returning to it. Unfortunately, this is not possible in a normal python program. If the experimenter stop the experimentation at a certain task, then she looses the execution state of that task (hidden by the python runtime environment) and cannot return directly to it latter. To get back to that task, she has to rerun the experimentation from the beginning. A rerun that is time consuming when the experimentation has heavy tasks such as deploying a specific Operating System on the testbed resources.

The EnosStack offers facilites to reifies the state of an experiment so that every task of the workflow is filled by a data structure that represents the execution state of that experiment. Executing a task modifies that data structure, that is then passed to the next task. The created data structure thus can be accessed by the experimenter, instead of being hidden in the python runtime environment. Thanks to that mechanism, an experimenter can easily recall a stopped task with its execution state, and thus restart the experimentation where she left it. For instance, the red-thread example exposes five tasks (see II). The first two one are in charge of getting resources from a testbed and deploying the application. When EnosStack executes the first task, it fills it with the resource model (see Fig. 2) and gets in return the concrete resources list. EnosStack then passes that resources list to the second task and gets in return, e.g., the IP addresses of deployed services. EnosStack stores the data structure in the filesystem and thus, make it accessible to the experimenter. So, if the experimenter hits ctrl+c in the middle of the second task, she can then recall that second task with its execution state (i.e., the concrete resources list) stored on the filesystem. Doing so proceeds the experimentation and save the need to restart the experimentation from the beginning. Also note that an experimenter can easily extend the EnosStack with other programming languages (e.g., C, Go, Ruby) by accessing and modifying the reified state.

To fast iterate on experiments, the EnosStack also favors idempotent tasks. In computing, an operation is idempotent if it can be applied many times without changing its result. For instance, calling twice the deployment task of the MOM service will do nothing the second time, since the MOM has been already deployed the first time. The result is a deterministic run, even when tasks imply side effects. Getting a deterministic run helps during the development phase of an experimentation. It lets the experimenter develop interactively on its experimentation until it works as expected. This, once again, speeds up the development phase.

Facilities provided by the EnosStack are idempotent. For instance, providers, that get testbed resources (see IV-A), are idempotent. Thus, calling a provider twice will do nothing the second time. More generally, the EnosStack relies on Ansible modules to favor idempotent tasks. An ansible module is an idempotent action that can be performed locally or remotely (e.g., copy a configuration file, start a container, execute a specific process, etc.). All together, Ansible modules form a Domain Specific Language (DSL) for scripting tasks that are idempotent by design.

V. MESSAGE ORIENTED MIDDLEWARE EVALUATION

Benchmarking a software component is a common experiment that can provide a deeper understanding about the component behavior, its scalability, its reliability, etc. For this example, we chose to evaluate the scalability of a messageoriented-middleware, namely RabbitMQ, to illustrate how a real experiment could be deployed thanks to the EnosStack. We describe this experiment in more details in the next section before presenting a related demo scenario and sharing the experiment material.

A. Experimentation

In order to illustrate the usefulness of the EnosStack, we present in this section some experiments we conducted in the context of the OpenStack performance working group. The experiments presented hereafter are part of a much broader test plan 5 that aims to evaluate the OpenStack message bus in the context of a massively distributed cloud architecture. We first detail some background about these experiments, and explain how we implemented them before presenting briefly some results.

1) Context: OpenStack is a cloud operating system, composed of many services that has to communicate between them. One way these services can communicate is through a Remote Procedure Call (RPC) system, built on top of a message bus. The OpenStack services relies on Olso Messaging, a library that provides a messaging API which supports RPC over a number of different messaging transports. In other words, Olso Messaging abstracts the underlying messaging middleware technology, allowing various messaging buses such as RabbitMQ, Qpid or ZeroMQ to be deployed. Being able to compare the performance of these different buses for different configurations would thus be highly valuable for the community. However, deploying such experiments in a repeatable way and at scale is a tedious task and can quickly become a time-consuming endeavor.

2) Implementation: We only describe here the experiments and results related to RabbitMQ as it is one of the most popular open source message bus. In order to stress the RabbitMQ bus, clients periodically sent RPC request messages to servers over the queue. To generate the RPC requests, we relied on ombt (Oslo Messaging Benchmarking Tool), an open source tool able to create RPC agents (clients or servers), and measure the latency and throughput of RPC transactions for each agent. In addition to these two metrics, we also collected numerous others, such as the CPU, memory, I/Os, etc. for every agent. Those metrics are collected using the monitoring service described in section IV-B. More widely, all the services were deployed in docker containers to ease the experiment management. Indeed, backuping the experiment consists in retrieving all the metrics and all the statistics from the RPC agents. This allows very deep post-mortem exploration or even in real-time. Destroying the experiment consists in removing all the running docker container as well as their attached states. Practically this means that the workflow (Fig. 1) can be repeated without the need to start fresh resources each time. This is particularly interesting in our cases because iterating over more than 50 sets of parameters were required. Thanks to the EnosStack, we were able to orchestrate in a fully automated way different executions of the experiment. 3) Results: We deployed these experiments on the Grid'5000 testbed on a cluster of 74 nodes, equipped with Intel Xeon E5 2.4Ghz, 128GB of memory and 10Gb connections. We emphasize that, leveraging the concept of providers in the EnOSLib, we first validated this workflow locally on Vagrant, before effectively deploying it on Grid'5000 on thousands of agents. This enabled us to greatly reduce the development time, while significantly easing the debugging process. We ran 50 experiments where we gradually increased the number of RPC clients from 250 to 2000 with a number of RPC servers equals to 250, 500, 750 and 1000. One node was entirely dedicated to host the RabbitMQ message bus, while another one was fully dedicated to orchestrate all the experiments. Each of the 72 remaining nodes was thus hosting either multiple RPC clients or servers agents. Each of the clients sent a total of 10000 RPC requests, waiting 0.1 seconds between each request to avoid overloading the servers. The average latency and throughput of these requests are respectively plotted in Fig. 3 and Fig. 4 depending on the number of clients and servers.

Results clearly indicate that both the latency and throughput linearly increase with the number of clients, up to 1000 clients. With more than 1000 clients, the throughput reaches a limit, while the latency increases dramatically due to RabbitMQ limitations. Correlated with the internal queue size of the bus and indirect measures like CPU and Memory we deduced that the bus was buffering a lot of messages when reaching this scale. It appeared, after further investigations that RPC servers blocked while handling requests leading to be considered as slow consumers from the bus perspective. In addition, we observed that the lower the number of RPC servers, the smaller the latency and the higher throughput. We suspect the bus delivery time to be dependent on the number of consumers a queue has, but it has to be confirmed in the near future.

These first experiments pave the way to larger scale benchmarks, deployed on different messaging bus technologies. Comparing these results to more distributed message bus for example is one of our close future work. Extending these experiments to other scenarios such as including network failures, or different communication patterns would provide developers a much better understanding of the messageoriented-middleware they aim to deploy.

B. Proposed Demonstration

We propose a live demonstration of the EnosStack. This demo will be based on the previous experimentation and will show how an experimenter can first validate its workflow in a local setup and then migrate it to an experimental platform (either Grid'5000 or Chameleon). We'll show how a list of parameters can be explored (each requiring a run of the workflow) while keeping all the wanted informations. Ultimately we'll give some hints on how the experimenter can create an environment allowing to analyse the results.

C. Reproducibility

As we put forward the importance of reproducible and repeatable experiments throughout this paper, we provide all the code and data that served to perform the experiments presented in section V. More specifically, the interested reader will find the followings:

VI. CONCLUSION

In this paper, we presented the EnosStack, an open source software stack we specially designed to assist experimenters. Deploying real experiments that are repeatable and reproducible remains challenging, and we argued that experimenters should benefit from a LAMP-like stack to help them in running their experiments. We provided precise requirements for such a stack, and detailed how the abstractions we conceived fulfilled

Fig. 2 .

 2 Fig. 2. Resource descriptions for Vagrant (left) and Grid'5000 (right)these roles can be used at the deployment time to segregate the network traffic.c) Example: Fig.2gives the declaration of the resources for the same experiment. Fig.2(a) will start virtual machines and get two networks from the local VirtualBox hypervisor. Fig.2(b) will deploy bare-metal servers on the Grid'5000 platform and get one isolated VLAN network (the second network is the default provided by Grid'5000). The description of machines is slightly different : size refers to a predefined amount of vCPU and Memory in the Vagrant case whereas cluster name is used in the Grid'5000 case. Also, network description for Grid'5000 allows to specify inhomogeneous network card configuration. More concretely here, specifying a count greater than one deploys several agents of the MOM in a clustered mode where the internal traffic (e.g replication traffic between RabbitMQ brokers) is isolated in a VLAN. Note also that the experimenter can give the two roles to the same network and traffic will not be isolated. This flexibility is allowed by the service abstraction presented in the next part.

Fig. 4 .

 4 Fig. 3. Average Latency

•

  The EnosLib source code on GitHub [?]. This code is released under the GPL3.0 • The python code to deploy the experiments, as well as the default configuration files, also on GitHub [?] under GPLv3.0. • All the docker images used to deploy services [?]. • The raw data results in JSON format, with all the parameters used [?]. • An online pre-written Jupyter notebook to extract meaningful results [?]. • The set of Gnuplot scripts and data files to plot the Figures in section V [?].

https://www.rabbitmq.com/

https://www.influxdata.com/time-series-platform/telegraf/

https://www.influxdata.com/time-series-platform/influxdb/

https://grafana.com/

https://docs.openstack.org/performance-docs/latest/test plans/massively distribute rpc/plan.html

them. Finally we presented a real experiment deployed on Grid'5000 and showed how the EnosStack greatly eased its deployment.

This work has been supported by the Discovery Inria project lab and the joint lab between Orange and Inria (CRE amqp evaluation). Experiments presented in this paper were carried out using the Grid'5000 testbed, supported by a scientific interest group hosted by Inria and including CNRS, RENATER and several Universities as well as other organizations (see https://www.grid5000.fr).