P Latil

L Orgéas

C Geindreau

P J J Dumont

S Rolland Du Roscoat

Towards the 3D in situ characterisation of deformation micro-mechanisms within a compressed bundle of fibres

Keywords: A -Polymer Matrix Composites, A -Fibre, A -Fabrics/textiles, D -X-ray microtomography, D -Rheology

The present study provides original experimental data concerning the evolution of the microstructure of a bundle of fibres during its deformation. For that purpose, a model saturated fibre bundle was processed and was subjected to a compression loading by using a specially designed micro-compression rheometer which was mounted on a synchrotron X-ray microtomograph. Thus, 3D images of the evolving fibrous microstructure could be obtained. Results first show that the compression induced both the bundle consolidation and liquid phase migration. Secondly, (i) the position, the orientation, the displacement and the deformation of each fibre together with (ii) the position and the evolution of each fibre-fibre contact were followed during the compression. The in situ tracking of these microstructure descriptors allows the consolidation micro-mechanisms to be analysed and provides useful information for theoretical or numerical models used to predict fibre bundle deformation during processing of fibre bundle reinforced polymer composites.

Introduction

In the past decades, the use of fibre-reinforced polymer composites in structural or functional applications has increased in various domains such as automotive or aerospace industries. Most of the fibrous reinforcements of these materials are made up of bundles containing few hundreds to few thousands continuous or discontinuous fibres. For this reason, these composite materials exhibit multiple scales: a macroscopic one, i.e. the part scale, a mesoscopic one, i.e. the fibre bundle scale, and a microscopic scale, i.e. the fibre scale. During their processing, e.g. sheet forming of long fibre reinforced thermoplatics, sheet forming of dry woven textiles for the Liquid Composite Moulding (LCM) processes, or compression molding of sheet moulding compounds, glass mat thermoplastics or hybrids, it is well known that the fibre bundles of these materials can be subjected to important displacements and deformations such as, for instance, orientation, consolidation, stretching, bending, twisting, shearing, bundle filamentation and breakage... [START_REF] Robitaille | Compaction of textile reinforcements for composites manufacturing. i: Review of experimental results[END_REF][START_REF] Buet-Gautier | Experimental analysis and modeling of biaxial mechanical behavior of woven composite reinforcements[END_REF][START_REF] Bickerton | The viscoelastic compression behavior of liquid composite molding preforms[END_REF][START_REF] Comas-Cardona | Unidirectional compression of fibre reinforcements. part 1: A non-linear elastic-plastic behaviour[END_REF][START_REF] Dumont | Processing, characterization and rheology of transparent concentrated fibre bundle suspensions[END_REF][START_REF] Badel | Simulation and tomography analysis of textile composite reinforcement deformation at the mesoscopic scale[END_REF][START_REF] Le | X-ray phase contrast microtomography for the analysis of the fibrous microstructure of SMC composites[END_REF][START_REF] Lomov | Full-field strain measurements in textile deformability studies[END_REF]. These bundle mesoscale deformation mechanisms are induced both by the complex relative motions and the deformation of fibres at the microscopic level within the bundles. They drastically affect the rheological behaviour of the composites during their forming by modifying, for instance, the permeability and the anisotropy of the fibrous networks [START_REF] Dumont | Processing, characterization and rheology of transparent concentrated fibre bundle suspensions[END_REF][START_REF] Lai | Model resin permeation of fibre reinforcements after shear deformation[END_REF][START_REF] Smith | The effect of shear deformation on the processing and mechanical properties of aligned reinforcements[END_REF][START_REF] Buntain | Compression flow permeability measurement: a continuous technique[END_REF][START_REF] Comas-Cardona | Unidirectional compression of fibre reinforcements. Part 2: A continuous permeability tensor measurement[END_REF][START_REF] Loix | Woven fabric permeability: from textile deformation to fluid flow mesoscale simulations[END_REF][START_REF] Loix | Flow of non-Newtonian liquid polymers through deformed composite reinforcements[END_REF]. They also affect final physical and mechanical properties of produced parts. For those reasons, mesoscale deformation mechanisms of fibre bundles have been deeply studied experimentally by using various observation and characterisation techniques (see above references). All these crucial experimental data are very useful and partially contribute to the (in)validation of some of the assumptions stated in meso-macro or micro-macro modelling approaches dedicated to this type of fibrous media [START_REF] Badel | Simulation and tomography analysis of textile composite reinforcement deformation at the mesoscopic scale[END_REF][START_REF] Toll | Packing mechanics of fiber reinforcements[END_REF][START_REF] Le Corre | Behavior of a net of fibers linked by viscous interactions: Theory and mechanical properties[END_REF][START_REF] Le Corre | Rheology of highly concentrated planar fiber supensions[END_REF][START_REF] Ausias | Direct simulation for concentrated fibre suspensions in transient and steady state shear flows[END_REF][START_REF] Lomov | Model of shear of woven fabric and parametric description of shear resistance of glass woven reinforcements[END_REF][START_REF] Alkhagen | Micromechanics of a compressed fiber mass[END_REF]. Among them, the X-ray microtomography technique which is now widely used in materials science [START_REF] Maire | On the application of x-ray microtomography in the field of material science[END_REF][START_REF] Rolland Du Roscoat | Estimation of microstructural properties from synchrotron X-ray microtomography and determination of the REV in paper materials[END_REF] is of particular interest since it can provide full 3D characterisa-tion of networks of fibres and their mesoscale deformation, even in composite materials with an optically opaque matrix. For instance, the orientation, the curvature, the compression and the breakage of discontinuous sized fibre bundles, occurring during the compression molding of an industrial SMC, were observed and quantified by using synchrotron X-ray microtomography with the phase contrast mode [START_REF] Le | X-ray phase contrast microtomography for the analysis of the fibrous microstructure of SMC composites[END_REF]. Similarly, by performing with dry woven fabrics in situ shear and biaxial tests inside a laboratory X-ray microtomograph, Badel et al [START_REF] Badel | Simulation and tomography analysis of textile composite reinforcement deformation at the mesoscopic scale[END_REF] have underlined the very large shape and surface changes the cross sections of continuous fibre bundles of deformed textile reinforcements could be subjected to.

By contrast and surprisingly, deformation mechanisms at the fibre scale, which are however responsible for the mesoscale bundle deformation, still remain not very well understood. What is the relative motion of fibres contained inside a bundle when it is subjected to a mesoscale consolidating/dilating and/or shear loading? Meanwhile, do fibres stretch, bend or twist? What is their orientation evolution? How many fibre-fibre contacts are there in the bundle?

Where are these contacts located? How do they evolve during the overall bundle deformation? What is their orientation? Can a fibre bundle be seen as a continuum at the mesoscale? Such type of crucial but still opened questions greatly hinder the development and the validation of the above meso-macro or micro-macro models. Discrete/finite element simulation performed at the fibre scale could be a relevant tool to answer them [START_REF] Le Corre | Rheology of highly concentrated planar fiber supensions[END_REF][START_REF] Ausias | Direct simulation for concentrated fibre suspensions in transient and steady state shear flows[END_REF][START_REF] Dumont | A numerical analysis of the evolution of bundle orientation in concentrated fibre-bundle suspensions[END_REF][START_REF] Switzer | Rheology of sheared flexible fiber suspensions via fiber-level simulations[END_REF][START_REF] Durville | Numerical simulation of entangled materials mechanical properties[END_REF]. However, this type of approaches must be previously validated and compared with experimental data collected at the fibre scale. The objective of this contribution is thus to propose an experimental methodology for the quantitative analysis of deformation mechanisms occurring at the fibre scale when a bundle of fibres is subjected to a mechanical loading. In this preliminary work, the methodology is built and illustrated by compressing a model bundle of fibres (section 2), but other mechanical loading conditions and other fibre bundles could be used in future. The compression experiment is performed with a micro-rheometrer that is mounted in a synchrotron X-Ray microtomograph, which allows in situ observations of the bundle microstructures to be performed during the test.

Therefrom, 3D relevant microstructure descriptors are quantified by using the 3D scanned images and image analysis subroutines (section 3). Thus, it is then possible to follow (section 4) and to discuss (section 5) the simultaneous evolutions of both the mesocopic mechanical response of the bundle and its associated microstructure descriptors.

Material and experimental procedure

In this preliminary study, the fibre bundle which was deformed was handmade. It consists of N = 69 rather straight and rather aligned fibres of length L = 11 mm and diameter D = 150 µm. These fibres were extracted from a fluorocarbon (PVDF) continuous elastic fishing wire (tensile Young's modulus E ≈ 2 GPa). The bundle was saturated with olive oil which behaves as a low viscosity Newtonian fluid at room temperature (shear viscosity µ ≈ 70 mPa s) and which is solid below 5 o C. This property is interesting since it allows the freezing of the bundle fibrous microstructure: this facilitates the bundle handling before mechanical testing which was performed at room temperature, i.e. when the oil is a fluid.

The fibre bundle was subjected to a simple compression loading at a low strain rate, i.e. ≈ 5 10 -3 s -1 , with a compression micro-rheometer which was purposely designed. As shown in the scheme displayed in figure 1, the bundle was put between an upper fixed plate and a lower one which was actuated vertically (with a piezo-motor not shown here). Before starting the test, the bundle was subjected to a very slight pre-compression, in order to ensure its upper and lower surfaces to be parallel with the rheometer plates. During this phase, no measurable force could be recorded.

The end of this pre-compression was arbitrarily chosen as the initial state. During the test, the upper and lower faces of the bundle were subjected to a compression force F 2 along the e 2 -direction which was recorded with a load cell (maximum capacity 500 N) together with the current height of the bundle l 2 . The other surfaces of the bundle remained stress free. The compression loading was interrupted at five compression strains ε 22 = ln(l 2 /l 20 ), l 20 being the initial height of the bundle. During these stops, slight decreases (≈ 10%) of the compression force were observed for ≈ 30 s [START_REF] Bickerton | The viscoelastic compression behavior of liquid composite molding preforms[END_REF], corresponding to a rearrangement of the fibrous microstructure and the fluid. Afterwards, the force remained constant and 3D images of the fibrous microstructure could be obtained. For that purpose, the micro-rheometer was installed in the X-ray microtomograph of the ID19 beamline at the European Synchrotron Radiation Facilities (Grenoble, France). At each compression stop, a fast scan of the bundle was achieved (overall scanning time ≈ 1 min) consisting in 1500 X-ray 2D projections obtained from the incremental rotation along e 2 of the micro-rheometer with respect to the synchrotron X-ray source (overall rotation = 180 o , beam energy = 25 keV). After proper reconstruction of the 2D projections, 3D maps of the X-ray absorption coefficient inside the bundle could be obtained (volumes of maximal size of 2048 × 2048 × 512 voxels, voxel size = 7.5 µm). The maps were denoised by using suitable filtering operations which faded their noise/signal ratio induced by the very short exposure time. Then, these five 3D maps could be cut in order to get cross sections of the micro-rheometer and the bundle (see figure 1).

They can also be used to characterise the bundle microstructure (see below).

Image analysis and data post-treatment

To analyse quantitatively both the mesoscopic mechanical response and the microstructure of the bundle of fibres, the five scanned volumes were subjected to several image analysis operations achieved with the softwares ImageJ and Matlab. These operations are described below and summarised in figure 2.

Volume of analysis -Mesoscale strains and stresses

The evolving volume of fibres V , (i) from which the mesoscopic mechanical response was estimated and (ii) within which the microstructure parameters of the bundle were extracted and analysed, was chosen in order to follow the mesoscopic deformation of the fibre bundle. After a first observation of the 3D images acquired during the test, it was found that the bundle flow was alike a lubricated compression, where the flow along the e 3 -direction was negligible compared to that in the (e 1 , e 2 ) plane. From this observation, V was chosen as a rectangular box (see figures 1 and 2(b)): its current height l 2 was given as the distance between the two plates of the micro-rheometer, its current width l 1 was set in order to contain all the fibres of the bundle and its length l 3 was set constant. As shown in figure 1(b), notice that due to the heterogeneity of fibres' lengths (induced by the manufacturing protocol), l 3 was such that it slightly and arbitrarily truncates (about 1 mm) the ends of some of the longest fibres which are not involved in the mechanical response of the bundle. It was then possible to estimate the mean strain the bundle was subjected to. For instance, the mean Hencki strain tensor ε could be estimated:

ε ≈ ε 11 e 1 ⊗ e 1 + ε 22 e 2 ⊗ e 2 with ε 11 ≈ ln l 1 l 10 and ε 22 ≈ ln l 2 l 20 (1)
together with its associated volume change ε v and deviatoric strain ε d :

ε v = ε : δ, ε d = 2 3 (ε - ε v 3 δ) : (ε - ε v 3 δ), (2)
where δ is the identity tensor. Furthermore, as the bundle was deformed at a very low strain rate with a very low viscosity saturating fluid, its compression could reasonably be considered as a drained compression during which interstitial fluid pressure effects are negligible: the monitored compression force was principally related to the microscale deformation of the fibres and their contacts. From this remark and from the estimation of the current volume of the bundle V , it is possible to estimate the mean Cauchy stress tensor of the bundle Σ as

Σ = Σ 22 e 2 ⊗ e 2 with Σ 22 ≈ F 2 l 1 l 3 . (3)

Detection / Description of fibres

Within the volume V , a standard thresholding operation was achieved in order to extract the fibre phase from the grey scale images (figure 2

choosing a volume V ′ centred in V and such that l ′ 1 = l 1 /3, l ′ 2 = l 2 and l ′ 3 = l 3 .
In order to detect the centrelines of all the fibres contained in V from the previous binarised images (see figure 2 Thirdly, the position of the centre of mass of each fibre cross section was calculated automatically and labellised (see figure 3(b-c)). Notice that due to the chosen resolution, the accuracy of their positions was about ±7.5 µm.

If this had negligible impact on the characterisation of fibre centrelines, this induced noticeable discrepancies on the detection of fibre-fibre contacts (see next subsection). Thirdly, the eleven detected centres of mass belonging to a given fibre i were used to form its centreline (this operation was done manually), here discretised by a piecewise linear 3D parametric curve x i (s i ), s i being the curvilinear abscissa of the considered fibre i (cf. figure 3(e)).

Therefrom, the fibrous medium was entirely reconstructed numerically (see

figure 2(d))
, by assuming that the fibre diameter D was constant. It was then possible to estimate numerically the direction of the lowest inertia axis pi of each fibre i. For the studied bundle, where the fibres are rather straight, pi can be considered as a relevant measure of the mean orientation of the fibre i.

Each centreline i was also associated with a local Frenet basis (t i , n i , b i ). This constitutes a possible way to characterise the local geometry of the centreline i (figure 3(f)). The tangent, normal and binormal unit vectors of this basis are respectively defined as (no summation on the indice i):

t i = dx i ds i , n i = 1 ||dt i /ds i || dt i ds i , b i = t i × n i . (4)
They are also linked by the Frenet-Serret formulae (no summation on i):

dt i ds i = κ i n i , dn i ds i = -κ i t i + τ i b i , db i ds i = -τ i n i , (5)
where κ i and τ i are the local geometrical curvature and torsion of the fibre centreline i at its curvilinear abscissa s i . In this work, a finite difference scheme was used to estimate the local Frenet bases, the geometrical curvatures and torsions from the discretised representation of the fibre cen-trelines. Lastly, let us remark that the geometrical torsion is not equivalent to the mechanical twist of the fibres. It is the rate of change of the osculating plane (t i , n i ) of the fibre centreline x i (s i ).

If it is different of zero, the centreline is not a 2D curve.

Results

Mesoscale mechanical behaviour

The mesoscale mechanical response of the bundle is summarised in figure 4.

This figure first shows that during its compression the fibre bundle is subjected to a pronounced consolidation. Indeed, its volumetric strain ε v is negative, its magnitude is as high as the deviatoric strain ε d and raises up to 0.45. This is accompanied with a drastic increase of the mean volume fraction of fibres Φ from Φ 0 = 0.29 to 0.45 and, consequently, with the migration of the fluid out of the bundle. Within the core of the bundle, i.e. in V ′ , a similar trend is observed, but the volume fraction of fibres is higher: during the test, Φ ′ goes from Φ ′ 0 =0.39 (showing the initial heterogeneity of the bundle) to 0.58. Finally, figure 4 also shows that the mean compression stress Σ 22 exhibits a significant hardening which is an increasing function of Φ.

In situ characterisation of fibres

(i) Orientation -Figure 5(a-b) gives the fibre orientation in the initial and the final steps. Two representations of the fibre orientation are used in this figure. The first one is the richest and uses unit spheres. On the surfaces of the spheres, spots represent the extremities of the N orientation vectors pi , their origin being located at the centre of the sphere. The second one, more compact, is the second order fibre orientation tensors A f [START_REF] Advani | The use of tensors to describe and predict fiber orientation in short fiber composites[END_REF], defined here in its discrete form as

A f = 1 N N i=1 pi ⊗ pi . (6)
As evident from figure 5(a), both representations show that in the initial state, fibres are mostly aligned along the e 3 direction, and that the disorientation is more pronounced in the plane (e 1 , e 3 ) than in the (e 1 , e 2 ) one. This figure also proves that these two trends are slightly accentuated during the compression.

(ii) Geometrical torsions and curvatures -The calculation of geometrical torsions τ i proved that the average geometrical torsion τ is negligible, i.e. below 10 -18 m -1 . This is not the case for the curvature. To illustrate this, the evolutions of the mean local radii of curvature Ri = 1/κ i with the position x 3 along the bundle have been plotted in figure 6(a) (κ i represents the average of the κ i estimated at x 3 for all the fibres in the bundle). This graph first shows that in the initial non-deformed state, fibres already exhibit finite curvature radii, with a mean value close to 50 mm. This is expected because the fluorocarbon fishing wire, from which the fibres were extracted, was initially rolled onto a 25 mm radius reel and thus presented an initial curved shape at rest, which mainly contributes to the initial curvatures observed in figure 6 [27]. The graph also points out that whatever the compression state, the curvature radii Ri are more or less constant and with values close to the average curvature radius R. ). However, in accordance with the previous point (see figure 5(c-d)), this graph proves that z′ is higher than z.

(iii) Orientation -The unit spheres displayed in figure 5(e-f), together with the second order contact orientation tensors A c defined as

A c = 1 N N j=1 qj ⊗ qj , (7)
first prove that contacts are mostly orientated in the (e 1 , e 2 ) plane: this is mainly due to the weak misalignment of fibres. Furthermore, these two orientation representations underline the preferential orientation of contacts parallel to the compression direction e 2 , whatever the considered compression state: for instance, A c 22 > A c 11 .

(iv) Contact geometry -Figure 7(c) gives the evolution of the dimensionless average contact distance d/ d0 with Φ: this distance slightly decreases during the first steps of compression, the decrease is sharper for the last two steps and reaches ≈ 4% at the end of the test. By the same time, the average contact length l/ l0 noticeably diminishes during the compression, up to 0.67 (see figure 7(d)), even if this trend seems to be reversed at the last compression step. It is interesting to notice that the error bars induced by the choice of h are smaller for these two descriptors than those obtained for z and z′ .

Discussion

The mesoscopic mechanical response of the studied bundle, which is characterised by marked consolidation and hardening effects (cf. figure 4), is a behaviour that is commonly observed during the compression of dry/saturated disordered/ordered fibrous media made up of contacting elastic fibres [START_REF] Robitaille | Compaction of textile reinforcements for composites manufacturing. i: Review of experimental results[END_REF][START_REF] Comas-Cardona | Unidirectional compression of fibre reinforcements. part 1: A non-linear elastic-plastic behaviour[END_REF][START_REF] Toll | Packing mechanics of fiber reinforcements[END_REF].

In particular, the experimental data plotted in figure 4(b) can be fitted by the following power-law type expression which is often proposed [START_REF] Toll | Packing mechanics of fiber reinforcements[END_REF]:

Σ 22 = kE(Φ m -Φ m 0 ). (8)
The two constitutive parameters k and m were here set to 54.8 and 14.25, respectively: the value of m is consistent with those reported in the literature for similar fibrous media, i.e. 7 ≤ m ≤ 16 [START_REF] Toll | Packing mechanics of fiber reinforcements[END_REF]. Beyond these phenomenological considerations, it is important to notice that trends reported in figure 4 are directly connected to the actual microstructure of the bundle and its evolution, which is ruled by micro-mechanisms arising at the fibre scale. By using standard micro-meso theories dedicated to discrete media (granular or fibrous), it would be possible, in principle, to account for these microscale effects and predict their mesoscale behaviour. For instance, the mesoscale stress tensor Σ of a fibrous medium containing N c fibre-fibre contacts inside a representative volume V f can be expressed as [START_REF] Toll | Packing mechanics of fiber reinforcements[END_REF][START_REF] Le Corre | Behavior of a net of fibers linked by viscous interactions: Theory and mechanical properties[END_REF][START_REF] Le Corre | Rheology of highly concentrated planar fiber supensions[END_REF][START_REF] Alkhagen | Micromechanics of a compressed fiber mass[END_REF]:

Σ = 1 V f Nc j=1 ξ j ⊗ f j , (9)
where f j is the contact force at the contact j, between a fibre i of centre of mass G i and a fibre k of centre of mass G k , and where ξ j = G i G k . In regard to the current experimental work, this expression brings up the following comments:

(i) The role of the fibrous microstructure on the mesocopic behaviour is obviously underlined from [START_REF] Lai | Model resin permeation of fibre reinforcements after shear deformation[END_REF]. In particular, such an expression shows that stress levels depend on the relative positions of contacting fibres, i.e. ξ j , which in turn depends on the fibre geometry (aspect ratio, cross section, waviness), content and orientation. Similarly, (9) also emphasises the role of the number of contacts N c in V f , or equivalently the coordination numbers Z in V f . The method proposed in this study is a possible "experimental way" to get estimates of these descriptors (e.g. figures 5-7). Nonetheless, this has to be achieved cautiously. For example, it was shown that the bundle displayed heterogeneous fibre content along its width and that this heterogeneity was preserved during the consolidation: this partially explains why the mean coordination number in V ′ , i.e. z′ , is higher than the one in V , i.e. z (cf. figure 7(a-b)). Edge effects at the lateral free surfaces x 1 = ±l 1 /2 could also contribute to this observation.

The method also provides experimental evidences to (in)validate microstructure models of the literature. Among them, one of the simplest is the soft core tube model [START_REF] Toll | Packing mechanics of fiber reinforcements[END_REF][START_REF] Le Corre | Rheology of highly concentrated planar fiber supensions[END_REF][START_REF] Toll | Note: On the tube model for fiber suspensions[END_REF][START_REF] Vassal | Upscaling the diffusion equations in particulate media made of highly conductive particles. Part II: Application to fibrous materials[END_REF]. By assuming that (a) fibres are uniformly distributed in the space and (b) can be seen as straight cylinders that can overlap, the tube model provides, for instance, an estimate of the average coordination number Z in an volume V f containing a sufficient number of fibres N f :

Z = 4φ   2 π r 1 N 2 f N f i=1 N f k=1 ||p i × pk || + 1 N 2 f N f i=1 N f k=1 |p i • pk | + 1   , (10)
This expression emphasises the role of the fibre aspect ratio r = L/D, orientation pi and content Φ on Z. It was used to estimate the average coordination numbers z and z′ , by setting respectively in it φ to Φ and Φ ′ , N f to N and N ′ , and by using the measured value of the orientation vectors pi in V and V ′ , respectively. If the experimental trend is more or less reproduced, figure 7(a-b) shows that the analytical tube model [START_REF] Smith | The effect of shear deformation on the processing and mechanical properties of aligned reinforcements[END_REF] largely overestimates the experimental coordination numbers. This is mainly due to the tested bundle which does not fulfil the hypothesis (a) of the tube model. Indeed, within V and V ′ , fibres are not homogeneously distributed spatially, in particular along the e 3 -direction. This trend is enhanced in V since noticeable variations of the fibre content are recorded (e.g. Φ < Φ ′ ). Besides, experimental estimations of z and z′ were achieved with a low number of fibres and by accounting for possible edge effects at fibre-plate contacts. This is further enhanced for z, since additional edge effects might occur at the free surface of the bundle, i.e. at (ii) The role of the deformation micro-mechanisms in ( 9) is also emphasised from the contact forces f j which contribute to the transfer, from fibre to fibre, of the mesoscopic compression loading. Firstly, contact forces induce the bending of fibres between consecutive contacts when fibres are not properly aligned: this explains the trends emphasised in figures 2(e) and 6. Secondly, while increasing the mesoscopic loading, load transfers from fibre to fibre also lead to deform contact zones and increase, in case of elastic fibres, the magnitude contact forces f j . These deformation mechanisms, together with the increase of the coordination numbers z or z′ (figure 7(a-b)), contribute to the trends reported in figure 4, i.e. the severe hardening of stress levels, the bundle consolidation and the drainage of the fluid. In the present study, it is interesting to notice that such mechanisms may be quite complex. For example, figure 7(c-d) shows that the contact distance d/ d0 and, more surprisingly, the contact length l/ l0 were decreasing during the compression. The last decrease could be due to the slight misalignment of fibres observed during the compression (figure 5(a-b)).

Conclusion

The preliminary results that have been obtained in this study are very encouraging. By combining a mechanical test on a fibrous medium, in situ 3D observations of its evolving microstructure and suitable image analysis subroutines, they prove that it was possible to characterise simultaneously meso but also micro deformation mechanisms in this fibrous medium. We have shown

with some examples that obtained results could provide useful information to understand the mechanics of fibrous media, but also to validate theoretical frameworks dedicated to the modelling of their microstructure and deformation. Future work will enter more deeply in this process, (i) by studying various types of fibrous media, (ii) by subjecting them to various types of mesoscale mechanical loading, (iii) by estimating quantitatively the mesoscale stress tensor ( 9).

  (a)), as illustrated in the binarised 3D views shown in figure 2(b). Hence, a first relevant microstructure descriptor, i.e. the volume fraction of fibres Φ, was simply estimated by the ratio [number of voxels of the fibrous phase] / [number of voxels contained in V ]. To analyse possible microstructure heterogenities, the volume fraction of fibre Φ ′ in the core of the bundle was also estimated similarly, by

  (c)), a procedure which is very efficient for almost parallel fibres with convex cross sections was used. Firstly, as depicted in figure3(d), eleven planar slices parallel to the (e 1 , e 2 ) plane and separated by a distance of 110 voxels were extracted from the scanned volumes. Secondly, on each 2D slice, adjoining fibre cross sections were separated by a standard 2D watershed algorithm without major difficulties, as shown in figure3(a-b).

3. 3

 3 figure 2(e)), (ii) its average contact distance dj , its average spatial position ȳj and (iii) its average orientation qj (see figure2(e)) (the indice "0" related to l j 0 and d j 0 in figure2(e) and in the following correspond to the initial state).

Figure 6 (

 6 b) lastly proves that the local bending of fibres is drastically increased during the compression, since R is decreased from approximately 50 mm to 20 mm at the end of the test. This can also be observed qualitatively by looking at, for instance, the two contacting fibres shown in figure2(e): these fibres are much more bent locally in the final compression state than in the initial one.4.3 In situ characterisation of contacts(i) Spatial distribution -As illustrated qualitatively in figure5(c) and (d), fibre-fibre contacts are preferentially located in the core of the bundle, i.e. within V ′ , where the fibre content is higher.(ii) Number -The number of contacts in V , N , substantially increases during the compression, from 61 to 239 (when h = D). To illustrate this, figure7(a)shows the evolution of the average number of contacts per fibre or equivalently the average coordination number z = 2n c /n f as a function of Φ, n c = N /V and n f = 4Φ/πD 2 L being the numbers of contacts and fibres per unit volume, respectively: z undergoes a non linear increase from ≈ 2.2 to ≈ 8.2 as Φ goes from 0.29 to 0.45. If the trend is preserved, noticeable shifts are observed by using h = D ± 7.5µm, as emphasised by the error bars. Similar conclusions can be drawn by considering the average coordination number z′ in V ′ (figure7(b)

x 1 =

 1 ±l 1 /2 (see above). By forgetting the assumption (a) of the tube model, i.e. by reconstructing numerically a bundle with N straight overlapping cylinders, the centres of mass and the orientations of which correspond to those measured in the tested bundle, figure7(a-b) shows that predictions of the tube model (without assumption (a)) are much closer to experimental values. However, due to the simplicity of assumption (b), i.e. soft core model plus no bending of fibres, the predicted coordination numbers remain overestimated.

Fig. 1 .

 1 Fig. 1. Scheme of the micro-rheometer and the model fibre bundle (reconstructed slices) which were put inside the ESRF X-ray microtomograph: vertical cross section in a plane parallel to (e 1 , e 2 ) (a), horizontal cross section in a plane parallel to (e 1 , e 3 ) (b).

Fig. 2 .

 2 Fig. 2. Illustration of the various image analysis operations used to estimate microstructure descriptors. The left pictures concern the initial non-deformed state, the right ones the last step of compression: reconstructed 3D grey scale volumes (a), volumes of fibres V after fibre segmentation (b), fibre centrelines (c), volumes of fibres V after reconstruction (d), detection of a fibre-fibre contact between two contacting fibres (e). In the left picture (e), two sections (A-A and B-B) are sketched in order to illustrate distances d j 0 and h used to detect fibre-fibre contacts (here with h = D).

Fig. 3 .

 3 Fig. 3. Method used to detect and characterise the fibre centrelines in V : (a) example of one of the 11 parallel slices shown in (d); (b) same slice after a 2D watershed; (c) same slice with the centres of mass of the cross sections; (d) generating the centreline of one fibre; (e) view of all centrelines; (f) local Frenet basis associated to one centreline.

Fig. 4 .

 4 Fig.[START_REF] Comas-Cardona | Unidirectional compression of fibre reinforcements. part 1: A non-linear elastic-plastic behaviour[END_REF]. Mechanical response at the mesoscale: evolution of the mean stress Σ 22 with the volume fraction of fibres Φ in V (outer graph), evolution of the mean deviatoric strain ε d as a function of the mean volumetric strain ε v . Around the graphs, 3D segmented volumes obtained at the five compression steps have been represented.

ΦFig. 6 .

 6 Fig. 6. Evolution (a) of the curvature radii Ri as a function of the position x 3 along the bundle for various compression states, (b) of the mean curvature radius R with the fibre content Φ.

Fig. 7 .

 7 Fig. 7. Evolution of the average coordination numbers z (a) and z′ (b) with the fibre contents Φ (a) and Φ ′ (b). Evolution of the dimensionless mean fibre-fibre distance d/ d0 ( d0 ≈ 145 ± 7.5 µm) (c) and the dimensionless mean contact length l/ l0 ( l0 ≈ 900 ± 7.5 µm) (d) with the fibre content Φ.

Acknowledgements -This work was performed within the ANR research program "3D discrete Analysis of micromechanisms of deformation in highly concentrated fibre suspensions" (ANAFIB, ANR-09-JCJC-0030-01) and the ESRF Long Term Project "Heterogeneous Fibrous Materials". The authors also thank the french region Rhône-Alpes for its financial support.


Denis Barthou
email: denis.barthou@labri.fr

Gilbert Grosdidier

Michael Kruse
email: michael.kruse@inria.fr

Olivier Pene
email: olivier.pene@th.u-psud.fr

Claude Tadonki
email: claude.tadonki@mines-paristech.fr

QIRAL: A High Level Language for Lattice QCD Code Generation

Quantum chromodynamics (QCD) is the theory of subnuclear physics, aiming at modeling the strong nuclear force, which is responsible for the interactions of nuclear particles. Lattice QCD (LQCD) is the corresponding discrete formulation, widely used for simulations. The computational demand for the LQCD is tremendous. It has played a role in the history of supercomputers, and has also helped defining their future. Designing efficient LQCD codes that scale well on large (probably hybrid) supercomputers requires to express many levels of parallelism, and then to explore different algorithmic solutions. While algorithmic exploration is the key for efficient parallel codes, the process is hampered by the necessary coding effort.

We present in this paper a domain-specific language, QIRAL, for a high level expression of parallel algorithms in LQCD. Parallelism is expressed through the mathematical structure of the sparse matrices defining the problem. We show that from these expressions and from algorithmic and preconditioning formulations, a parallel code can be automatically generated. This separates algorithms and mathematical formulations for LQCD (that belong to the field of physics) from the effective orchestration of parallelism, mainly related to compilation and optimization for parallel architectures.

Introduction

Quantum Chromodynamics (QCD) is the theory of strong subnuclear interactions [START_REF] Wilczek | What QCD Tells Us About Nature -and Why We Should Listen[END_REF]. Lattice QCD (LQCD) is the numerical approach to solve QCD equations. LQCD simulations are extremely demanding in terms of computing power, and require large parallel and distributed machines. At the heart of the simulation, there is an inversion problem:

Ax = b, (1)
where A is a large sparse (also sparse and implicit) matrix, called Dirac matrix, b is a known vector and x is the unknown. To model reality accurately, A should be of size around 2 32 × 2 32 elements at least. Current simulations on supercomputers handle sizes up to 2 24 × 2 24 elements. High performance LQCD codes on multicores, multinode and hybrid (using GPUs) architectures are quite complex to design: performance results from the interplay between the algorithms chosen to solve the inversion and the parallelism orchestration on the target architecture. Exploring the space of algorithms able to solve the inversion, such as mixed precision or communication-avoiding algorithms, aggressive preconditioning, deflation techniques or any combination of these is essential in order to reach higher levels of performance. Many handtuned, parallel libraries dedicated to LQCD, such as Chroma [START_REF] Edwards | The Chroma Software System for Lattice QCD[END_REF] or QUDA [START_REF] Clark | Solving lattice QCD systems of equations using mixed precision solvers on GPUs[END_REF] propose building blocks for these algorithms, so as to ease their development. Parallelism does not stem from the algorithm themselves but from the structure of the sparse matrices involved in the computation. However, there are a number of shortcomings to the library approach: First, any significant evolution of the hardware requires to tune the library to the new architecture. This can lead to change the grain of parallelism (as an adaptation to cache sizes for instance) or change the data layout. Then, sparse matrices for LQCD, in particular the matrix A, have a very regular structure and library functions take advantage of it. Combining different algorithms, such as preconditioners, leads to structural changes in matrices that are not supported by these libraries. These two limitations hinder considerably the time necessary to develop a parallel/distributed code. We propose in this paper a high-level Domain-Specific Language, QIRAL, to model LQCD problems and enable automatic parallel code generation. This novel approach makes possible the exploration and test of new algorithms for LQCD. Sparse matrices in LQCD can be structured using algebra operators on dense matrices and the key idea is to use this structure to express parallelism. This extends ideas proposed in SPIRAL [START_REF] Püschel | Encyclopedia of Parallel Computing[END_REF], the library generator for DSP algorithms to more complex codes and matrices. The problem formulation is first presented in Section 2, then the language is described in Section 3 and its implementation in Section 4.

Lattice QCD Description

The principle of LQCD is to discretize space-time and describe the theory on the resulting lattice. The lattice has to contain a very large number of sites since it has to be fine grained and describe large enough volumes. LQCD exists since 1974. Huge progress in hardware, software and algorithms have been achieved, but it is not enough to really sit on the parameters of nature. The light quarks, "u" and "d" are still described by quarks heavier than in nature while the heavy "b" quark is described as lighter. This implies systematic errors in our results. To break this limitation, several orders of magnitudes will be needed in the computing power and related resources. It will demand new hardware with several level of parallelism, and make coding more and more complicated. Our goal is to provide tools helping to face this new situation.

The heaviest part in the LQCD calculation is to generate a large Monte-Carlo sample of very large files (field configurations) according to an algorithm named "Hybrid Monte-Carlo". It is a Markovian process, every step of which takes several hours on the most powerful computers. The algorithm is complex but it spends most of its time in inverting very large linear systems which depend on the field configuration. Typically one will deal with matrices with billions of lines and columns. The second heaviest part in the calculation is to compute "quark propagators" which again boils down to solving large linear systems of the same type [START_REF] Lüscher | Local coherence and deflation of the low quark modes in lattice qcd[END_REF]. Therefore we will concentre on this task. The matrices involved in these computation represent 4D stencil computations, where each vertex is updated by the value of its 8 neighbours. Their structure is therefore very regular, statically known, and it is used in the following section.

A High-Level Domain-Specific Language

QIRAL is a high level language for the description of algorithms and for the definition of matrices, vectors and equations specific to LQCD. The objective of the algorithmic part is to define algorithms independently of the expression of sparse matrices used in LQCD. The

Dirac = I L⊗C⊗S + 2 * i * κ * µ * I L⊗C ⊗ γ 5 + κ * d∈D ((J -d L ⊗ I C ) * s∈L U (d)[s]) ⊗ (I S + γ[d]) + κ * d∈D ((J d L ⊗ I C ) * s∈L U (-d)[s]) ⊗ (I S -γ[d])
Figure 1: Definition of Dirac matrix on a Lattice L in QIRAL.

objective of the system of definitions and equations is to define properties and structure sparse matrices in order to be able to find parallelism. While the algorithmic part uses straightforward operational semantics, the equational part defines a rewriting system.

For the sake of simplicity, QIRAL is a subset of L A T E X, meaning that the QIRAL input can be either compiled into a pdf file, for rendering purposes, or compiled into executable codes. Algorithms and definitions correspond to different predefined L A T E X environments.

Variables, Types: Variables and constants used in QIRAL are vectors (denoted by the type V) of any length, matrices of any size (M), complex (C) or real numbers (R). Besides, counted loops are indexed by variables of type index, iterating over a domain (denoted indexset). Index variables are potentially multidimensional, and integers are a particular case of index value. Functions of any number of argument of these types can be also defined. The size of a vector is defined through the size of its index set. The size of an index set can be left undefined. If V1 is the vector and IS is an index set, V1[IS] denotes the subvector indexed only by IS.

Matrices and vectors cannot be manipulated or defined element-wise. Instead, matrices and vectors are built using either constant, predefined values such as identity I IS (for the identity on an index set IS), or operators such as +, -, * and the transposition, conjugate, direct sum and tensor product. The tensor product ⊗ and direct sum ⊕ are defined by:

A ⊗ B = [a ij B] ij , A ⊕ B = A 0 0 B .
These two operators define parallel operations and help defining the sparse matrices of LQCD.

Definition and Equations: All LQCD knowledge is given as a system of definition (for matrices and vectors) and a system of equations on these values. In particular, the structure of sparse matrices, explicitly defining where non-zero elements are located, is given in these definitions. This structure is then propagated in the algorithms. Algorithms, Preconditioners: Many algorithms proposed for solving Equation 1 are Krylov methods. Figure 2 presents two algorithms in QIRAL. The first is a variant of the conjugate gradient and is representative of iterative methods. The second one, Schur complement method, is a preconditioner: Computing x, solution of Equation 1 is achieved by computing two solutions to smaller problems. The keyword Match helps to define the algorithm as a rewriting rule: Whenever the statement in the Match condition is found, it can be rewritten by the algorithm. For the SCHUR preconditioner, this rewriting can be performed at the condition that the requirement is fulfilled: P 1 AP t 1 has to be invertible. Using identities defined through the equation system, the rewriting system proves here that indeed, the requirement is valid, when matrix A is Dirac matrix and P 1 is a projection matrix keeping only vertices of the lattice with even coordinates (even-odd preconditioning). QIRAL is expressive enough to represent different variants of con-

Algo 1: Conjugate Gradient, Normal Resolution [CGNR] Input : A ∈ M, b ∈ V, ǫ ∈ R Output : x ∈ V Match : x = A -1 * b Var : r, p, Ap, z ∈ V, α, β, nr , nz , n z1 ∈ R r = b ; z = A † * r; p = z ; x = 0 ; nz = (z | z) ; nr = (r | r) ; while (nr > ǫ) do Ap = A * p; α = nz /(Ap | Ap) ; x = x + α * p ; r = r -α * Ap ; z = A † * r ; n z1 = (z | z) ; β = n z1 /(nz ) ; p = z + β * p ; nz = n z1 ; nr = (r | r) ;
Algo 2: Definition of Schur complement method [SCHUR] 4 From QIRAL to parallel code

Input : A, P 1 , P 2 ∈ M, b ∈ V Output : x ∈ V Match : x = A -1 * b Var : v 1 , v 2 , x 1 , x 2 ∈ V, D
D 12 = P 1 * A * P t 2 ; v 1 = P 1 * b ; v 2 = P 2 * b ; x 2 = (D 22 -D 21 * D -1 11 * D 12 ) -1 * (v 2 -D 21 * D -1 11 * v 1 ) ; x 1 = D -1 11 * (v 1 -D 12 * x 2 ) ; x = P t 1 * x 1 + P t 2 * x 2 ;
The QIRAL compiler takes as input a file describing algorithms and equations and generates parallel code. The phases are described in 3 and presented in detail in the following. Rewriting System Generation: From its current L A T E X implementation, the QIRAL input is translated into a set of rewriting rules and equations, using Maude framework [START_REF] Clavel | The Maude System[END_REF]. Maude is a multi-purpose rewriting system, handling equational rewriting rules and reflection. More precisely, definitions are translated by a first phase, tex2maude into equations (or conditional equations) and algorithms into rewriting rules. The output file contains domain specific definitions (e.g. the Dirac operator from Figure 1), the desired program (in our case x = Dirac -1 b), and all algorithms and preconditioners we want to apply.

Definitions have a semantic similar to a rewriting rule: the left hand side is rewritten into the right hand side. To ensure the existence of a normal form, the system of equation has to be convergent and confluent. So far, our implementation does not use a tool such as an automatic Church-Rosser Checker [START_REF] Durán | A church-rosser checker tool for conditional order-sorted equational maude specifications[END_REF] . Algorithms are translated into conditional rules, applied is their prerequisite are checked. The left hand side of the rule correspond to the Match clause and the algorithm itself is the right hand side.

This rewriting system is merged with another one defining general algebraic properties, code generation rewriting rules and code optimizations. The main phase of the QIRAL compiler is therefore described through rewriting system, following works such as Stratego [START_REF] Visser | A core language for rewriting[END_REF]. QIRAL is a static strongly type language, vectors and matrices are defined by their index sets. Type checking is the first analysis achieved by this main phase.

Applying Algorithms and Simplifications:

The user provides the list of algorithms to compose and the initial program (here, we focused on the equation x = Dirac -1 b). Exploring different algorithms and sequences of algorithms boils down to change the list given to the QI-RAL compiler. Checking that algorithm requirements is automatically achieved by the rewriting system, using the equational theory provided through the system of definitions. The equational system, using both equations coming from LQCD definitions and algebraic properties involving the different operators, simplifies terms that are equal to zero.

The result is a program where the matrix A of the algorithms has been replaced either by the Dirac matrix, or by a matrix obtained through transformation by preconditioners.

Loop generation and parallelization: At this step, the initial statement has been replaced by the algorithm statements, directly using the Dirac matrix or a preconditioned version. Assignment statements are vector assignments: Values for all vertices of the lattice can be modified in one statement. Matrices are still described using tensor products and direct sums.

Loops are obtained by transforming all indexed sums, products into loops or sequences. For instance, the direct sum operator indexed by the lattice L in the definition of Dirac matrix, in Figure 1 is transformed into a parallel loop over all elements of the lattice. This loop is parallel, due to the meaning of ⊕. As the lattice is 4D, either 4 nested loops are created, or one single, linearized loop is created. The choice depends on a parameter in the QIRAL compiler.

Some usual compiler transformations are then used, computing dependences, fusioning loops, applying scalar promotion and other optimizations. At the end of this phase, an OpenMP code is produced, essentially by identifying for parallel loops the set of private variables. These transformations are driven by rewriting strategies, using reflection in Maude.

Matching Library Calls: The resulting code still uses some high level operators on dense matrices, such as tensor products on dense matrices, matrix-vector product. These operators are then replaced by library calls. For this step, it is sufficient to define for each library the expression it computes, as a rewriting rule. We defined LION, a set of hand-written library functions used for validation purposes.

Finally, post-processing phase, mostly syntactic rewriting, transforms this output into a C function. This function is compiled with a program skeleton and a runtime that initializes data, calls the generated code and stores the result. Custom functions can also be defined here

Conclusion

This paper has presented an overview of QIRAL, a high level language for automatic parallel code generation of Lattice QCD codes. The language is based on algorithmic specifications and on mathematical definitions of mathematical objects used in the computation.

The contribution of this short paper is to show that a high-level representation makes possible the automatic generation of complex LQCD parallel code. Parallelism directly stems from the structure of sparse, regular matrices used in LQCD. While the initial matrix represents a stencil computation, QIRAL is able to manipulate more complex structures obtained through preconditioning for instance, unlike Pochoir [START_REF] Tang | The pochoir stencil compiler[END_REF]. The approach, similar to the one proposed by Ashby et al. [START_REF] Ashby | Cross Component Optimisation in a High Level Category-Based Language[END_REF] enables the user to define new equations and domain-specific definitions. The QIRAL compiler is able to keep such information through transformations resulting from preconditioners or algorithms. The tensor product and direct sum operators are translated into parallel loops and lead to OpenMP code generation. This way, algorithmic exploration, the key for higher levels of performance, can be freed from the constraints and costs of parallel tuning. Besides generation of distributed codes with communications is within reach. For heterogeneous architectures, such as Cell or GPUs, further work for automatic data-layout optimization is required, in order to reach levels of performance of previous works (such as [START_REF] Vranas | Massively parallel quantum chromodynamics[END_REF][START_REF] Khaled | Implementing Wilson-Dirac operator on the cell broadband engine[END_REF][START_REF] Tadonki | An efficient CELL library for lattice quantum chromodynamics[END_REF] for the CELL BE and [START_REF] Clark | Solving lattice QCD systems of equations using mixed precision solvers on GPUs[END_REF] for the GPU).

Figure 1

 1 defines the sparse Dirac matrix. Declarations of index sets L, C, S, D are not shown here, as well as the declaration of the constant matrices used here.The direct sums indexed by the vertices of the lattice (L) define block-diagonal matrices. These diagonals are then shifted by d or -d columns by a permutation matrix J d L (predefined).

Figure 2 :

 2 Figure 2: Two algorithms. On the left: A variant of the conjugate gradient method; on the right: A preconditioner, the Schur complement method. jugate gradient, BiCGSTAB, methods with restart (such as GCR), and other methods proposed by physicists.

Figure 3 :

 3 Figure 3: Overview of QIRAL compilation chain

Figure 4 :

 4 Figure 4: Sample output C code produced for CGNR algorithm

Figure 5 :

 5 Figure 5: Convergence speed for a lattice of size 4 4 for different methods: Conjugate Gradient Normal Error (CGNE), Normal Resolution (CGNR), Biconjugate Gradient (BiCGSTAB), Modified conjugate gradient with preconditioning (MCR1).

  11 , D 12 , D 21 , D 22 ∈ M

	Require : invertible(P 1 * A * P t 1 ) D 21 = P 2 * A * P t 1 ; D 11 = P 1 * A * P t 1 ; D 22 = P 2 * A * P t 2 ;