Cyril Sola

Bruno Castanié

Laurent Michel

Frédéric Lachaud

Arnaud Delabie

Emmanuel Mermoz

On the role of kinking in the bearing failure of composite laminates

Keywords:

Fibre kinking is one of the main failure modes of composite laminates under compression loading. In this paper, the role of kinking in the failure of quasi isotropic composites subjected to a bearing load is inves tigated. High resolution CT scans show that kinking is largely involved in the events leading to laminate collapse, notably by triggering other damage modes such as delamination. Kink bands develop extremely progressively, leading to the formation of a wide localization zone (or FPZ, failure process zone). Such behaviour calls for a non local modelling approach. Local damage models would lead to overly conserva tive sizing. A simple model, based on Hashin failure criteria and non local effective stresses is confronted to experiments, and its limits are highlighted. It will be shown that proper modelling of the bearing fail ure requires the characteristic behaviour of kink bands to be taken into account.

Introduction

Bearing failure has been widely studied, both experimentally and numerically [1]. However, the complex nature of this failure mode still makes it a challenging test case for damage models [2]. Most experimental studies have focused on the bearing beha viour of pinned or bolted composite joints, with typical bearing strength (Eq. (1)) ranging from 400 MPa to 1000 MPa.

r bearing F t Á d ð1Þ
Thus, the boundary conditions pertaining to each joint configu ration appear to have considerable influence on the bearing strength. The most influential ones are the lateral restraint due to bolt preloading [3], the bolt hole clearance [4], or secondary bending for single lap shear joints. The damage mechanisms involved in the bearing failure of bolted joints have been studied by Xiao and Ishikawa [5], who showed that kinking (often referred to as fibre microbuckling) plays an important role in the failure process, as did Seike et al. [6]. Camanho et al. [7] showed that the progressiveness of the bearing failure could be explained by the accumulation of subcritical damage. They emphasised the role of delamination, suggesting the use of three dimensional failure criteria, which they did in a later paper [8]. However, cohesive zone models do not seem to be an absolute requirement [9]. Despite the many models proposed in the literature, which are often based on traditional failure criteria [8,10 12], the sizing of composite joints is still mainly based on allowable bearing strengths obtained for different joint configurations.

Hence, there still is a need for improved failure models built upon a sound physical basis. Mutual feedback should be estab lished between the numerical models and the experiments to look beyond the load displacement curve. To do this, it would be desir able to use very simple experiments, such as pinned contact [13] or half hole pinned bearing [14 17], in order to reduce the number of unknown or uncertain variables and to focus on the damage mech anisms. A more general knowledge of the bearing failure would also make it possible to diverge from the traditional bolt/rivet mechanical assemblies in specific industrial applications.

In this paper, bearing failure and the role of kinking are further investigated using a half hole pinned bearing test and modern experimental techniques such as computed tomography (CT) and digital image correlation (DIC), for both woven and UD fabric lam inates, with the same epoxy matrix. An enhanced engineering fail ure model is then introduced, exhibiting large similarities with most of the existing failure theories dedicated to bearing failure simulation. The enhancement mainly stems from the use of an embedded non local approach, which prevents an excessive impact of subcritical damage on the failure load. However, the aim of this paper is less to propose a new failure theory than to point out the deficiencies of such theories, and to suggest other The damage variables are then replaced by their regularized counterparts.

These damage variables are then condensed into two damage variables, related to fibre damage or to matrix damage. They will be used to translate the effect of damage on the compliances. No specific damage variable is used to model shear behaviour. This simplified approach is similar to that proposed by Gutkin and Pinho [START_REF] Gutkin | Proceedings of the 18th international conference on composite materials, ICCM18[END_REF].

d f maxðd ft ; d fc ; d N f Þ ð 8Þ
d m maxðd mt ; d mc ; d N m Þ ð 9Þ
The superscript N refers to the previous converged iteration. This definition of damage ensures that damage can only increase or remain stable once it has initiated.

The damaged compliance matrix then reads:

e S 1 E 1 ð1 d f Þ m 12 E 1 m 13 E 1 0 0 0 m 12 E 1 1 E 2 ð1 dmÞ m 23 E 2 0 0 0 m 13 E 1 m 23 E 2 1 E 3 ð1 dmÞ 0 0 0 0 0 0 1 G 12 ð1 dmÞ 0 0 0 0 0 0 1 G 23 ð1 dmÞ 0 0 0 0 0 0 1 G 13 ð1
Its inverse, the stiffness matrix, is also the tangent operator, since damage evolution is not related to stresses or strains.

In the proposed model, stresses in the failure criteria are replaced by effective stresses after the initiation of damage:

r ij rij ð11Þ
where r ij

r ij 1 d f if i j 1 and r ij r ij
1 dm if i; j 2 ½2; 3. Non local enhancement. The non local approach is used to eval uate the failure criteria in a non local way, much as in the FCV method proposed Hochard and coworkers [START_REF] Hochard | A ply scale non-local fibre rupture criterion for CFRP woven ply laminated structures[END_REF]. The FCV approach is built upon an external subroutine because a user material model is not able to handle several integration points simultaneously (the routine loops through every integration point in a sequential way). In the current model, the procedure can be directly embedded in the finite element code by using the non local capabilities of Sam cef,2 which makes it possible to establish a ''dialogue" between neighbouring Gauss points. Using the built in implicit gradient model, it is possible to compute non local variables in a user mate rial routine by simply enhancing the stress vector and the tangent matrix. In contrast with the FCV, the implicit gradient enhanced approach is not a simple averaging, but a weighted average similar to the Gaussian averaging of traditional integral non local approaches [START_REF] Peerlings | A critical comparison of nonlocal and gradient-enhanced softening continua[END_REF]. In fact, the implicit gradient enhanced model is a diffusion equation, relating a local variable f to its non local counter part f and involving an internal length c. It is associated with an equation enforcing zero flux at the boundaries of the material domain considered, X, which ensures that a homogeneous field remains unchanged.

f c 2 Df f in X rf :ñ 0 o n @X 8 < : ð12Þ
Samcef allows only one internal length to be used. Thus, the internal length cannot be fine tuned to each damage mode, which would be preferable [START_REF] Germain | Composite layered materials: anisotropic nonlocal damage models[END_REF].

The fact that damage is instantaneous, denoting perfectly brittle behaviour, might seem to contradict the idea that damage should be non local to simulate the quasi brittle nature of composite materials. In fact, in this model, the non local approach, applied to the failure criteria only, is used as a simple way of emulating the effect of the progressive failure without actually needing a complete set of unconventional material parameters, which are tricky to identify and most often not available to the engineer at the design stage.

In the non local approach, effective stresses r are replaced by the effective non local stresses r when evaluating failure criteria.

Non local averaging can be applied to any component of the stress tensor. Alternatively, non local strains could be computed. The constitutive relation between stresses and strains remains local, which preserves the original elastic behaviour.

Results

A typical coupon was modelled using hexahedral elements, with one layer of elements per ply in the refined zone of interest. The pin was assumed to be perfectly rigid, which improved calcu lation times while providing results very close to a flexible/flexible contact. The friction coefficient was taken to be equal to 0.1. This value was selected on the basis of an analysis of the DIC results focusing on the shapes of the different strain fields (shear and transverse strains in particular). It has no direct physical meaning, and should merely be seen as an average of the local friction coef ficients, which would vary depending on the local angle between the fibres and the slip directions.

Local model. The local model predicts a peak strength that is much lower than the actual one. Since damage almost instanta neously increases to 0.99, the laminate has no further possibility to absorb (by deforming elastically) or dissipate (by failing pro gressively) strain energy, and it thus fails immediately.

It should be noted that in the model proposed by Camanho and Matthews [8], which is based on the same failure criteria as the current model, damage is deemed to have much less influence on the components of the damaged stiffness tensor. For example, the fibre stiffness would only decrease by 86% upon failure initia tion under compression, whereas in the present model, the dam aged stiffness is virtually reduced to 0 when damage is total. While, regarding stiffness, Camanho's approach is probably closer to reality, as compressive damage appears unlikely to lead to a complete stiffness loss, its major drawback is that it enables con siderable energy storage in the laminate, since non completely damaged elements would continue to load elastically. Whether this energy storage is real or artificial remains to be debated, espe cially when one considers the kink band broadening phenomenon, which supposedly takes place at constant stress. What is certain is that such an approach, although physically motivated, greatly delays the numerical peak strength, which can even be non existent (Fig. 14).

Non local model. The non local model delays initiation depend ing on the internal length value, effectively limiting the impact of the first subcritical damage on the failure load by simply ignoring it. It was found that, to obtain a numerical peak strength close to the experimental value, the model should be fully non local, mean ing that all the failure criteria should be evaluated using non local effective stresses, with a very large internal length of 1.5 mm.

FEA package commercialised by Siemens.


A Diet

M Biancheri-Astier

Y Le Bihan

C Conessa

K Gbafa

F Alves

M Grzeskowiak

M Benamara

G Lissorgues

A Pozzebon

A switched reader complementary-loops structure for detecting LF RFID tagged pebbles

Keywords: LF RFID, complementary loops, mutual coupling

This paper presents a reader switched complementary loops structure for pebbles RFID tracking. These so-called "smart pebbles" are tagged with LF RFID glasstags and deployed on beaches, in the context of coastal morpho-dynamics survey. At different periods, the detection of these pebbles over the beach has to be done in the minimum of time. The goal of the project LARGE is to provide a RFID reader antenna structure for reducing the time of the glasstags detection. Technical challenges are that LF glasstags detection is very sensitive to its orientation. As the prototype has to adapt to an unknown orientation and depth of the RFID tags, two key ideas are exploited in this work: (i) the use of complementary loops for reducing the orientation sensitivity and (ii) the switching in loop sizes for adapting to optimal depth of tag detection. A prototype is build and tested with a commercial RFID reader and conclusions of our measurements show that detection performances can be increased with this multiple loops and complementary loops structure.

I. THE CONTEXT OF PROJECT "LARGE"

RFID technology is widely used for tracking objects and or animals in various environments, including applications in an outdoor context [START_REF] Paret | Identification Radio-Fréquences et Cartes à puce sans contact[END_REF][2][3] [START_REF] Grzeskowiak | Pebbles Tracking Thanks to RFID LF Multiloops Inductively Coupled Reader[END_REF]. The operating frequency between the RFID reader/interrogator and the transponders (tags), is the discriminant factor due to multiple tradeoffs mainly concerning the detection range, data rate transfer and sensitivity to the environment (water, metal…). The RFID chosen band of frequency consequently determines the electromagnetic physical link to be in near-field (magnetic coupling at LF 125/134 kHz, or HF 13,56 MHz) or alternatively in far-field (backscattering of radiating waves, usually in UHF 433/860 MHz and sometimes at 2,4 GHz).

In this paper, we focus on tracking pebbles, on the beach, as a part of the project called LARGE (LocAlisation with Rfid of tagGed pEbbles). The key idea is to insert on the beach some tagged pebbles and to recover them at different time slots in order to build a behavioral model of their displacements due to the coastal morpho-dynamics. As salt water and sand are present, the choice of LF frequency band is preferred. LF RFID tags are dedicated technology for this type of study in that it can be efficiently used to track the movements of sediments on coarse grained beaches, for both the under and outside sections of a beach. As a consequence, LF RFID tags can be embedded inside autochthonous pebbles turning them into tracers that can be positioned directly on the beach. These tracers can be then detected, identified and possibly recovered to study their movements and variations in terms of size and shape. These data are fundamental to evaluate the effects of coastal erosion on a beach and to plan the countermeasures to be deployed [START_REF] Benelli | An RFID-based toolbox for the study of under-and outside-water movement of pebbles on coarse-grained beaches[END_REF]. LF is the only frequency band that can be employed in such a scenario in that is the only one that is able to provide a Long Range communication channel even under sea water, allowing the detection of Smart Pebbles also in the submerged portion of a beach. Tags on pebbles has a minimum volume and the dedicated structure is a glasstag type which is typically made of a 1-1,5 cm long ferrite rod, with diameter in the range of 1-2 mm with enameled copper wire soldered to a RFID chip. As resumed in Figure 1, The RFID tag in LF enable to store and read data in presence of salt water and sand, and consequently to date and localize precisely the movement of the pebble (the use of such data for building the morpho-dynamics behavioral model is not the subject of the paper). In LF, most of the tags does not have battery and need being powered by magnetic coupling before communicating thanks to the load-modulation (see Figure 1). Moreover, the penalty of the glasstag is it high sensitivity to the orientation because the ferrite rod is equivalent to a magnetic dipole with a very small effective surface due to the rod diameter value. As the pebbles can be positioned randomly, the RFID reader loop antenna is supposed to detect them whatever their orientation. This is technologically challenging because the tag need to be powered with enough magnetic field, fruitfully distributed over the equivalent magnetic dipole surface, before being able to be detected. To improve the tag powering ability, high quality factor coils can be used, but the quality factor also reduce the bandwidth of the system and can penalize the correct interpretation of the data, preventing the detection. Consequently, a tradeoff is present on the quality factor of the transfer function between the tag and the reader, which is based on the impedance value seen by the reader, including the mutual inductance [START_REF] Grzeskowiak | Pebbles Tracking Thanks to RFID LF Multiloops Inductively Coupled Reader[END_REF][5] [6][7]. This mutual inductance can be defined by the Neumann formula, seen in Figure 2, where M is the value, in Henry, between two closed paths Γ1 and Γ2 where r12 is the distance between the two elementary vectors dl1 and dl2. The mutual inductance is exclusively based on geometrical considerations: the position and orientation of the paths (loops in classical examples). In this Figure 2, we suppose that the reader loop area is higher in size than the tag loop area, to the contrary of case 1, because of the surface of detection, which determine the total time spend to control entire beach surface. This drives to a very week value of the mutual inductance in a case of single loop configuration (the case 2). As expressed in Figure 2 (case 3 and equation) it is possible to separate the reader loop path defined by Γ1 into loop sub-paths (multi-loop or series loops) Γ1 a and Γ1 b which will drive to the addition of the two partial mutual inductances parts, one of which being potentially better suited for mutual coupling with a given geometry, for example those of the RFID glasstag.

Additionally, the tag orientation implies to study the problem of mutual coupling maximization for the two possible orientation mode of the glasstag equivalent surface, which are the horizontal and vertical ones regardless to the plane defined by the reader loop surface [START_REF] Grzeskowiak | Pebbles Tracking Thanks to RFID LF Multiloops Inductively Coupled Reader[END_REF][5] [START_REF] Imura | Maximizing air gap and Efficiency of Magnetic Resonant Coupling for Wireless Power Transfer using Equivalent Circuit and Neumann formula[END_REF]. The ability to detect in both horizontal and vertical modes is more fruitful in the case of pebbles detection than to increase the detection range in only one mode.

The idea developed in that paper is to improve the structure of the reader loop by modifying the mutual coupling in the case of a glasstag. The multi-loop principle is explored in that way and combined with complementary loops principle [START_REF] Grzeskowiak | Pebbles Tracking Thanks to RFID LF Multiloops Inductively Coupled Reader[END_REF][5], for different size.

II. STRUCTURE OF THE RFID READER LOOP

The problem of detecting the tagged pebbles on the beach can be simplified by an equivalent schematic shown in Figure 3, where the beach is represented as a rectangular area.

A : Wide control area:

-Fast (meas. Campaign) -Low coupling with tags whatever the depth low recovery rate In Figure 3, three scenarii are proposed in function of the surface of control of the RFID system. If the surface is too wide (case A), the detection can be penalized by the weak magnetic coupling between the reader loop and the tag coil. If this surface is too small (case C), the detection can suffer from the two high relative distance between the tag equivalent area and the reader loop area, decreasing the magnetic coupling. Also, the time spend to perform the control of the entire beach is directly proportional to the inverse of the surface area. Optimal case B suppose that the depth of the pebbles is evaluated and help in designing the area size of the reader loop. As it is difficult to evaluate precisely a reference depth, an interval of 4 cm to 40 cm is supposed, implying us to find a solution involving the combination of different reader loop sizes. Also the relative orientation of the glasstag can be overcome thanks to the complementary-loops structures for which the possibility of mutual coupling in both horizontal and vertical modes are higher than those of the single and multiple "sub-loops in-phase" ones.

A. Complementary-loops structure

The orientation of the tag brings a challenge on the detection performances. The complementary loops structure can help in the reduction of sensitivity due to the mode of detection. To illustrate this assumption, simulations of the mutual induction generated between a small circular tag and rectangular reader loops are performed with MATLAB. The calculus consists in generating the magnetic induction distribution in the space and summing the vectorial components orthonormal to a given delimited surface, following the definition of the induced magnetic flux. In Figure 4, the simulation setup is shown and we identified three scenarii with the same rectangular reference surface of 10x20 cm². The surface (a) corresponds to the classical single turn rectangular loop. Surface (b) and (c) corresponds to the use of sub-loops of 10x10 cm² for, respectively, in-phase loops and complementary loops (i.e. out-of-phase) structures. In order to quantify the ability of such structures for detection the tag at a fixed distance (4 cm herein), we defined the factor S as the sum of the mutual inductance absolute value for a given displacement of tag. This S factor is evaluated for quantifying, at a given distance of the reader loop, the amount of mutual induction available for potential position of the tags. Results of mutual coupling evaluation and S factor are reported in Figure 5 and Figure 6. The results reported in Figure 5 and Figure 6 shows that the case (c) corresponds to the best case in regards of the S factor. In horizontal mode, the case (b) does not improves the case (a) because the distribution of the magnetic field are destructive at the middle of the surface, preventing a potential mutual coupling. To the contrary, the case (c) corresponds to a constructive addition of the magnetic field generated by the two complementary loops, in horizontal and vertical modes. Additionally, the vertical mode does not suffer from the used nul of magnetic coupling at the middle of the rectangular surface, as seen for case (b) in Figure 6. The values of S factor, increasing for the complementary loops structure, are due to the fact that extrema of the "M" values are higher for the case (c) not only at specific position but over the surface. Nevertheless these results are only concerning a fixed distance between the tag and the reader loops, it can be noticed the important benefit of using complementary loops for reducing the sensitivity to the tag orientation.

B. Switching the reader loops of different sizes

It was shown in the previous part that complementary loops structures are fruitful, at a given distance of detection. As was studied in [START_REF] Grzeskowiak | Pebbles Tracking Thanks to RFID LF Multiloops Inductively Coupled Reader[END_REF] [START_REF] Diet | Improving LF Reader Antenna volume of detection for RFID token tag thanks to Identical Coaxial Loops (ICLs) and in/out-of phase multiple-loops structures[END_REF], the optimal size of a rectangular reader loop depends on the tag size and the distance fixed for the detection, that is to say the depth in our context. As the depth can be from 4 to 40 cm, we decide to increase the surface of the reader loop and to superpose different reader loops whose structures are made of complementary loops, as it is represented in Figure 7. The three structures chosen are all included in a reference surface of 40x40 cm (as it is the maximum depth) and we fixed the classical single turn 40x40 cm² loop and its sub-division by a factor 2 and 4, respectively, the 4x(20x20) cm² C-loop (complementary loop) and the 16x(10x10 cm²) C-loop. The magnetic field simulations, with CST studio, of these structures are seen in Figure 8. The simulations results in Figure 8 are obtained for an arbitrary current of 1A and the color ramp is clamped at the same level for the three simulations. As it can be seen, the magnetic field magnitude is concentrated near the edges, and the curvature radius increase with the inverse of the sub-loop area. This potentially creates some zeros of detection for the larger loops, which can be compensated by using the sub-loops under the hypothesis of a not too deep tag. The prototype of such a structure involves turning the wire of the loop in order to respect the phase of the current from one sub-loop to another. The next part describes how the prototype is realized for combining switched reader loops and sub-loops.

III. PROTOTYPING AND TESTS

To benefit from the sub-loop structure detection performances (in terms of sensitivity to the tag orientation), the complementary loops are realized line by line and column by column following the method depicted in Figure 9. This method has the advantage of making a line of "n" sub-loops for which the electrical port is on the same side of the surface, which highly simplify the practical realization of the prototype and avoid some penalizing extr-length of wire (loss…). In this figure, the equivalent electrical schematic for the case of two and "n" sub-loops is draw in the idea to evaluate the analytical approximation of the inductance seen by the RFID reader. The formula given considers the mutual coupling of sub-loops only between two adjacent loops and neglect the other possible coupling. The result of this hypothesis is that the equivalent inductance seen by the reader for a line (or a column) of complementary sub-loops is proportional to the number of subloops.

As each reader loop, and each of the lines and columns of the complementary loops structure is switched, the evaluation of the inductance is mandatory to be adapted to the RFID reader. The prototype is realized with Litz wire (to reduce losses), and the inductance is adapted to 26 μH, which is the value given in the datasheet of the TIRIS 134 kHz RFID reader. Photography of the prototype realized over a wood structure is seen in Figure 10, in which the RFID reader and the different mechanical switches can be seen.

L eq(N) Results of detection are shown in Figure 11 for the two mode of detection, horizontal and vertical, and for three cases: (i) the single turn coil (40x40) and, in the idea to reduce the amount of data, only one line and one column of the two others structures (ii) and (iii). This selection of a line and a column can be clearly seen in the figure .   As was predicted in simulation with CST in the previous part, the wide loop presents some zeros of detection for the vertical mode, see the figure at top right and middle right. Also, the horizontal mode benefit from the small C-loops because the detection range near the edges of the reference surface is higher these structures than for the single turn. By mixing the results of horizontal and vertical modes, the ability to detect arbitrarily orientated glasstag is possible at several centimeters, taking into account the fact that the line and column can be additionally swept.

I 1 I 2 V 1 V 2 L 1 L 2 M L eq (2)     M L j M L L j I V V L eq        2 2 ) 2 (
V 1 V 2 L 1 L 2 M L eq (N) V 4 L 3 L 4 M V 3 M           I I I with 2Mjω 2M L jω n 2M 1 - n - L n jω I V V (n) L p k 1 1 n even p p n odd k k eq           
The prototype uses a BANANA PI board for driving the switches, enabling an automated sweep in the range of some seconds, fruitful for a practical scanning of a given surface of the beach in the measurement campaign. IV. CONCLUSION

In this work, we presented a switched structure of single and complementary loops, whose interest is to improve the detection performances of RFID tags detection. Practical results showed that complementary loops enables to detect the RFID tags both in horizontal and vertical modes, as was predicted by simulation. Additionally, the use of switched structures increases the detection performances of the RFID tags at different depths, which are not known in the context of the project. Practical sweeping of all the switches, in the range of several seconds, are done with a simple BANANA PI in the measurement. To conclude, RFID tag detection can be improved by using this structure, because the switched configuration enables the detection at different optimal depths, at different sub-areas, and for the two possible modes (Hm and Vm). This increases the possibility of detection over the area.

Figure 1 :

 1 Figure 1 : Principle of pebble tagged in RFID with LF glasstag.

Figure 2 :

 2 Figure 2 : definition of the Mutual coupling and possible vectorial decomposition of the path in a multi(sub)-loop case

BFigure 3 :

 3 Figure 3 : Equivalent problem of RFID tagged pebbles detection and potential scenarii A, B and C.

Figure 4 :

 4 Figure 4 : Simulation of the complementary-loops structure in the case of an equivalent circular surface (radius 1mm) and 3 reader loop cases. Definition of the S factor.

Figure 5 :

 5 Figure 5 : Mutual induction, in horizontal mode, for the three rectangular loops cases in Figure 4, and evaluation of the S factor.

Figure 6 :

 6 Figure 6 : Mutual induction, in vertical mode, for the three rectangular loops cases in Figure 4, and evaluation of the S factor.

Figure 7 :

 7 Figure 7 : Superposition of the different reader loops

Figure 8 :

 8 Figure 8 : Magnetic field magnitude, with CST studio, for the single turn 40x40 cm² loop, the 4x(20x20) cm² C-loop and the 16x(10x10 cm²) C-loop.



  What is L eq (N) for a line made of N complementary loops, with:-Lx = Ly -M < 0 -Mutual Coupling only with the next loop, other neglected

Figure 9 :Figure 10 :

 910 Figure 9 : method for turning a line of "n" complementary loops, and equivalent electrical schematic for determining the inductance value.i1 i2

Figure 11 :

 11 Figure 11 : Detection performances, in cm, of a RFID glasstag in horizontal (left) and vertical (right) mode, for the three structures: the single turn 40x40 cm² loop (top), the 4x(20x20) cm² C-loop (middle) and the 16x(10x10 cm²) C-loop (bottom)

ACKNOWLEDGMENT

The authors want to thanks Joachim Fert for the valuable work concerning the final prototyping.


J.-E Franzkowiak
email: jean-eloi.franzkowiak@ensam.eu

P Mercier

G Prudhomme

And L Berthe

Multiple light scattering in metallic ejecta produced under intense shockwave compression

A roughened metallic plate, subjected to intense shock wave compression, gives rise to an expanding ejecta particle cloud. Photonic Doppler velocimetry (PDV), a fiber-based heterodyne velocimeter, is often used to track ejecta velocities in dynamic compression experiments and on nanosecond time scales. Shortly after shock breakout at the metal-vacuum interface, a particular feature observed in many experiments in the velocity spectrograms is what appear to be slow-moving ejecta, below the free-surface velocity. Using Doppler Monte Carlo simulations incorporating the transport of polarization in the ejecta, we show that this feature is likely to be explained by the multiple scattering of light, rather than by possible collisions among particles, slowing down the ejecta. As the cloud expands in a vacuum, the contribution of multiple scattering decreases due to the limited field of view of the pigtailed collimator used to probe the ejecta, showing that the whole geometry of the system must be taken into account in the calculations to interpret and predict PDV measurements.

INTRODUCTION

As a shockwave reaches the free surface of a roughened metal, high-speed ejecta particles are produced from a limiting case (Atwood number A ≃ -1) of the Richtmyer-Meshkov instability phenomena [START_REF] Dimonte | Richtmyer-Meshkov instability in the turbulent regime[END_REF] and freely expand in a vacuum, constituting an expanding cloud of droplets. To measure the areal mass of ejecta, Asay foil [START_REF] Asay | Ejection of material from shocked surfaces[END_REF], x-ray [START_REF] Monfared | Experimental observations on the links between surface perturbation parameters and shock-induced mass ejection[END_REF][START_REF] Zellner | Effects of shock-breakout pressure on ejection of micron-scale material from shocked tin surfaces[END_REF], and piezoeletric pins [START_REF] Zellner | Exploring Richtmyer-Meshkov instability phenomena and ejecta cloud physics[END_REF] have been commonly used in dynamic compression experiments. Photonic Doppler velocimetry (PDV), a fiber-based heterodyne velocimeter [START_REF] Strand | Compact system for high-speed velocimetry using heterodyne techniques[END_REF][START_REF] Moro | Defining parametric dependencies for the correct interpretation of speckle dynamics in photon Doppler velocimetry[END_REF][START_REF] Warren | Simulating translation-induced laser speckle dynamics in photon Doppler velocimetry[END_REF], is employed to measure particle velocities in their direction of expansion and is a promising diagnostic for a non-invasive characterization of ejecta particles. This effort investigates the multiple scattering of light in the ejecta and its influence on time-velocity spectrograms measured by PDV. What appear to be slow-moving ejecta detected by PDV shortly after breakout (below the free-surface) are explained by collisions among ejecta particles in Ref. [START_REF] La Lone | Spall Strength and Ejecta Production of Gold Under Explosively Driven Shock Wave Compression[END_REF]. As already presented in Ref. [START_REF] Andriyash | Optoheterodyne Doppler measurements of the ballistic expansion of the products of the shock wave-induced surface destruction: experiment and theory[END_REF] in the case of a micro-spall particle cloud (particle sizes in the range of 10-100 μm and velocity gradients in the range of 300 m/s), we suggest an explanation to this observation in the case of a micro-jetting particle cloud (particle sizes in the range of 0.5-10 μm and velocity gradients in the range of 1000-2000 m/s), based on the multiple scattering of light in the ejecta. The solution is obtained using a Doppler Monte Carlo (MC) approach, taking into account the transport of polarization in the ejecta cloud. Furthermore, we link the contribution of multiple scattering to the geometry of the PDV pigtailed collimator, justifying why this particular signature decreases with time, and showing that the whole geometry of the system must be taken into account to interpret precisely PDV measurements in metallic ejecta experiments.

METHODS

In a PDV setup, a collimated and coherent light beam at λ 1.55 μm and frequency ω 0 is used to illuminate an expanding metallic particle cloud. Doppler-shifted scattered fields are collected in the pigtailed collimator in the backscattering direction and mixed with a reference beam (local oscillator) with frequency ω LO ω 0 and intensity P LO , using a heterodyne configuration, to generate an intensity modulation I t in the photodetector plane. The square modulus of the short-term Fourier transform (STFT) calculated on I t is the optical power spectral density (PSD), so-called time-velocity PDV spectrograms in the shock community. To emphasize how the multiple scattering of light impacts PDV results, we perform Doppler MC simulations. Doppler MC is a powerful tool, often employed in tissue optics and laser Doppler flowmetry [START_REF] Ramella-Roman | Three Monte Carlo programs of polarized light transport into scattering media: part I[END_REF][START_REF] Jentink | Monte Carlo simulations of laser Doppler blood flow measurements in tissue[END_REF][START_REF] De Mul | Laser Doppler velocimetry and Monte Carlo simulations on models for blood perfusion in tissue[END_REF], to simulate complex systems. The Doppler spectrum, originating from the interaction of a coherent light beam with moving scatterers, can be recovered using such statistical methods. We employed this method in our study due to the high complexity of the ejecta cloud: the polydispersity must be taken into account, as well as the spatial evolution of volume particle densities, from the head of the cloud to the location of the free surface.

A. Description of the Cloud

A particle cloud evolving freely in vacuum along z is defined by its particle size and velocity distributions (resp., f d and f v ), between diameters d min and d max , and between free-surface velocity V fs and V max . f v V d V is the approximate number of particles between velocities V and V dV , and f d d p d d p the number of particles between diameters d p and d p d d p . The integrated areal mass-velocity function of ejecta M V (kg∕m 2 ) between V and V max is defined by

M V ≃ ρ p hV p i S Z V max V f v V dV ; (1)
with ρ p (kg∕m 3 ) the material density and S (m 2 ) the surface of ejection. The total areal mass of ejecta is M s M V fs . hV p V i (m 3 ) is an average particle volume determined from the normalized size distribution f d :

hd p i 3 Z d max d min d 3 p f d d p dd p 6 π hV p i: (2)
We have assumed that particle sizes and velocities are uncorrelated. This may not be exactly true, and a dependence on velocity may be introduced to account for correlations between particle sizes and velocities, d max increasing with decreasing velocities [START_REF] Sorenson | Ejecta particle size measurements from the break-up of micro-jets in vacuum and helium gas using ultraviolet in-line Fraunhofer holography[END_REF]. However, this dependence is very small, and we assume uncorrelated sizes and velocities will not change the results presented in this paper.

A pigtailed collimator (lens diameter Φ p ) is used to illuminate ejecta particles in their direction of propagation z. Scattered fields are collected in the same aperture in the backscattering direction. The solid angle Ω within which photons are collected is limited by the fiber core diameter Φ c and the focal length f 0 of the probe through the relationship Ω πΦ 2 c ∕4f 02 π tan α 2 . Below a critical distance d c Φ p ∕2 tan α, the optical coupling efficiency of the scattered fields in the system is constant. Shortly after breakout, when liquid metal sheets have broken up into micrometer-sized particles, the cloud of length L is discretized into N slabs of width δz L∕N , to each of which is associated a velocity (in decreasing order) V i ∈ V max ; V fs and a scattering mean free path 1∕μ i (m) defined, ∀ i ∈ 1; N and using the independent scattering approximation, by

μ i 3hQ ext iδM V i 2ρ p δzhd p i ; ( 3
)
with δM V i M V i1 -M V i (kg∕m 2 )
and

hQ ext i Z d min d max f d d p Q ext d p dd p ; (4)
the average extinction efficiency determined using Mie theory of light scattering (spherical particles). Q ext d p is the extinction efficiency for a single spherical particle of diameter d p . Equation (3) is valid if particle volume fractions do not exceed ∼1%, otherwise correlations between scatterers will imply a nonlinear dependence of μ i with particle densities [START_REF] Garca-Valenzuela | Insights into the dependent-scattering contributions to the extinction coefficient in highly scattering suspensions[END_REF].

B. Doppler Monte Carlo Simulations

In the preceding section, the optical properties of the ejecta particle cloud have been defined and MC simulations can be performed to simulate the interaction. A photon is launched with initial weight W 1 1, direction 0; 0; 1 u 1;x ; u 1;y ; u 1;z u 1 , position x 1 ; y 1 ; 0, and null Doppler shift V d (m/s) (Fig. 1), x 1 and y 1 being sampled using a Gaussian distribution of waist ω. A pseudo-random number η between 0 and 1 is sampled and j initialized to 1. Aslog η remains smaller than δz P j i1 μ i , j is incremented by one and the distance l s traveled before a first scattering by δz. When the inequality is violated,

l s j -1δz - 1 μ j log η δz X j-1 i1 μ i ! ; (5)
a scattering occurs on a particle of diameter d p;1 and velocity

V 1 -V 1 u z .
Assuming a power-law dependence [START_REF] Sorenson | Ejecta particle size distributions for shock loaded Sn and Al metals[END_REF] for particle sizes f d ∝ d -α c p , d p;1 is directly sampled using the inverse cumulative distribution function and the random deviate η:

d p;1 d 1-α c max V -d 1-α c min η d 1-α c min 1 1-α c : (6)
To get the direction of the photon u 2 after a first scattering, the scattering properties of the particle must be defined. The single scattering matrix is expressed as

Mθ 0 B B @ aθ bθ 0 0 bθ aθ 0 0 0 0 d θ -eθ 0 0 eθ d θ 1 C C A ; (7)
with

1 2 jS 1 j 2 jS 2 j 2 ; bθ 1 2 jS 2 j 2 -jS 1 j 2 ; d θ 1 2 S 2 S 1 S 1 S 2 ; eθ i 2 S 2 S 1 -S 1 S 2 :
S 1 and S 2 define the complex amplitude scattering matrix calculated using Mie theory. Note that matrices Mθ are calculated for particle diameters ranging from d min to d max before the MC simulation starts. To incorporate the transport of a polarized light beam in the ejecta cloud, we follow the Euler MC method described in Ref. [START_REF] Ramella-Roman | Three Monte Carlo programs of polarized light transport into scattering media: part I[END_REF]. For an incident Stokes vector S 0 I 0 ; Q 0 ; U 0 ; V 0 , the scattering phase function p s θ; ϕ is defined by

p s θ; ϕ aθI 0 bθQ 0 cos 2ϕ U 0 sin 2ϕ: (8)
At each photon-particle interaction site, the polar and azimuthal angles θ and ϕ are sampled using a rejection method [Eq. ( 8)]. The new direction of the photon is given by

u 2 R θ 1 R ϕ 1 u 1 , R θ 1 being
the rotation matrix of angle θ 1 in 3D space around the y axis and R ϕ 1 the rotation matrix of angle ϕ 1 around the z axis. The weight W 1 is decreased using an implicit scheme to account for part of the radiation being absorbed by the metal:

W 2 W 1 1 -Q abs d p;1 ∕Q ext d p;1 . The Doppler shift in velocity V d (m/s) is also incremented according to V d V d 0.5u 2 -u 1 • V 1 .
The incident Stokes vector S 0 is transformed into a scattered Stokes vector S upon scattering. As described in Ref. [START_REF] Ramella-Roman | Three Monte Carlo programs of polarized light transport into scattering media: part I[END_REF], using the 4 × 4 matrix defined by

Q ϕ 0 B B @ 1 0 0 0 0 cos 2ϕ sin 2ϕ 0 0 -sin 2ϕ cos 2ϕ 0 0 0 0 1 1 C C A ; (9)
we get

S Mθ 1 Q ϕ 1 S 0 ; (10)
defined in the local Euler reference frame of the photon. From this position, a new distance l s is sampled to determine the location of the next scattering event, etc. The free surface is treated as a perfect reflector moving at velocity V fs -V fs u z .

To speed up the calculation, a variance reduction technique is used: as scattering events remain within the field of view (FOV) (see Fig. 1) of the PDV probe (ω), we store the probability that, after each scattering, the photon is collected in the backscattering direction. The photon being scattered n e times in the cloud within the FOV, for k ∈ 1; n e , the probability of detection, corrected from an exponential attenuation factor, is

P k Y k l1 1 - Q abs;l Q ext;l exp - X j k i1 μ i δz p k π -θ k : (11)
j k is the index of the layer in which each scattering occurs, and p k is the probability of backscattering for the k-th scattering on a particle of diameter d p;k . θ k is the angle between -u z and u k .

If the photon propagates towards the probe (decreasing z), p k π -θ k is replaced with p k θ k in Eq. ( 11). To P k is associated a given Doppler shift V d;k (m/s), accumulated by the photon over the scattering path (Doppler frequency ω k ):

V d;k λ 4π ω k -u z u k • V k 2 1 2 X k-1 l1 u l1 -u l • V l ; (12)
and a given Stokes vector S k , calculated after k scattering events. At each interaction site, we determine the backscattered Stokes vector. A last transformation, which rotates S k to the meridian plane of the probe, must be computed to get the final Stokes vector S k;f in the reference plane (Oxy), shown in Fig. 1. This rotation is computed using the local Euler basis (see Ref. [START_REF] Ramella-Roman | Three Monte Carlo programs of polarized light transport into scattering media: part I[END_REF]).

The backscattered complex field amplitude Et can be described by a sum over all different paths k followed by the photons in the ejecta:

Et ∝ X k ffiffiffiffiffi P k p expiω 0 ω k t ϕ k e k;f ; (13)
where e k;f is a unitary vector defining the polarization state of each backscattered photon associated to the Stokes vector S kf , P k is given by Eq. ( 11), and ϕ k and ω k are, respectively, the phase and the Doppler frequency associated to each path (ϕ k are uniformly distributed between -π and π). The intensity of the local oscillator P LO being much larger than the intensity scattered by the cloud, we can neglect the homodyne interference term (mutual interferences between the scattered waves), as well as the components oscillating at optical frequencies 2ω 0 ω k (DC terms). The beating intensity I t measured by the photodetector is proportional to (heterodyne interference):

It ∝ 2 ffiffiffiffiffiffiffiffi P LO p X k ffiffiffiffiffi P k p expiω k t ϕ k -ϕ LO e k;f • e LO c:c:; (14)
where "c:c:" stands for "complex conjugate." e LO is the unit norm vector defining the polarization state of the local oscillator with constant phase ϕ LO . Following the reasoning presented in Ref. [START_REF] De Mul | Laser Doppler velocimetry and Monte Carlo simulations on models for blood perfusion in tissue[END_REF] (Appendix A), we define N the total number of frequency intervals in the spectrum, and we assume that all photon paths corresponding to the same frequency (number of paths ≫1) can be merged into a single Doppler frequency component ω j , as done in Ref. [START_REF] Andriyash | Optoheterodyne Doppler measurements of the ballistic expansion of the products of the shock wave-induced surface destruction: experiment and theory[END_REF]. Using the δ-Kronecker notation, Eq. ( 14) simplifies to

I t ∝ 2 ffiffiffiffiffiffiffiffi P LO p X N j1 X k ffiffiffiffiffi P k p expiω k t ϕ k -ϕ LO × δω k -ω j e k;f • e LO c:c: I t ∝ 2 ffiffiffiffiffiffiffiffi P LO p X N j1
ffiffiffiffi ffi P 0 j q expiω j t ϕ j c:c::

P 0 j and ϕ j (uniformly distributed between -π and π) depend on time t. The expectation value of P 0 j is given by hP 0 j i

X k P k C k δω k -ω j ; (16)
with C k ∈ 0; 1 the interference contrast between e k;f and e LO :

C k ke k;f • e LO k 2 : (17)
e k;f and e LO represent two unique points on the surface of the Poincaré sphere. If e k e LO , C k 1. If there are no preferential backscattered polarizations, a polarization diversity receiver may be used to mitigate speckle fluctuations in the PDV spectrum [START_REF] Goodman | Speckle Phenomena in Optics: Theory and Applications[END_REF]. The values of C k for two different polarization states are uniquely determined by the distance between these two points on the surface of the Poincaré sphere. After a few scatterings, the polarization states of the backscattered photons will be randomized, and C k will decrease to an average value of 1/2. Since all particles are moving in the cloud, a Doppler shift is calculated for each scattering event. A typical PDV probe (Φ c ≃ 8 μm, f 0 ≃ 2.5 mm) has an aperture angle α ≃ 0.1°, and it will be assumed that only backscattered photons contribute to the Doppler spectrum (cos α ≃ 1, an angle of 0.1°having a negligible impact on the Doppler measurement). A photon exiting two times the FOV is terminated in the calculation. The Fourier transform of the autocorrelation of the intensity signal I t ∝ EE can be calculated to determine the PSD (see Ref. [START_REF] De Mul | Laser Doppler velocimetry and Monte Carlo simulations on models for blood perfusion in tissue[END_REF]). Finally, the authors show in Ref. [START_REF] De Mul | Laser Doppler velocimetry and Monte Carlo simulations on models for blood perfusion in tissue[END_REF] that the heterodyne Doppler spectrum at a given frequency ω j is proportional to the optical energy detected in the corresponding domain. Using this result, and n ph photons being propagated in the ejecta in the simulation, the Doppler velocity spectrum Φ mc V is deduced from Eqs. ( 11), ( 12), ( 16) and ( 17), ∀ V :

Φ mc V X n ph i1 X n e i k1 P k C k δV -V d;k ! n ph →∞ hΦV i: ( 18
)
hΦV i is the expectation value of a typical noisy PDV spectrum. As the number of photons n ph propagated in the MC simulation increases, the fluctuations σ 2 Φ V hΦ mc V 2 i -hΦ mc V i 2 observed in the Doppler spectrum are reduced.

C. PDV Velocity Spectrograms

Independently, PDV diagnostic is often used to get a Doppler spectrum ΦV on a given ejecta cloud, and experimental results may be directly compared to MC results [Eq. [START_REF] Ramella-Roman | Three Monte Carlo programs of polarized light transport into scattering media: part I[END_REF]]. A noise-corrupted spectrum is recovered, due to shot noise from the photodetector and electronic and digitization noise. For a random particle arrangement, it was shown in Ref. [START_REF] Andriyash | Optoheterodyne Doppler measurements of the ballistic expansion of the products of the shock wave-induced surface destruction: experiment and theory[END_REF] that the probability density function of the Doppler spectrum ΦV measured by PDV, originating from the Gaussian fluctuations of the fields, follows a fully developed speckle statistics described by

PΦV 1 hΦV i hB Φ i exp - ΦV hΦV i hB Φ i ; ( 19
)
where hB Φ i is an additive background noise in the spectrum. The decorrelation time of intensity is directly the width of the Fourier window T (s) used to process the data, and the velocity decorrelation in the spectrum is δV λ∕2T . For a given experiment, the properties of the cloud are unknown, and it may be possible to calculate which Doppler spectrum hΦV i is the most likely to explain the observation, using Eqs. ( 11), ( 12), ( 17), and (18).

EXAMPLE

A single example is presented in this paper in Fig. 2. We demonstrate the relevance of performing Doppler MC simulations to illustrate particular features observed in PDV velocity spectrograms in shock compression ejecta experiments, and we show why the whole geometry of the system must be taken into account to interpret PDV data.

A. Properties of the Cloud

The areal mass-velocity distribution of ejecta M V is supposed to be

M V M s exp -β V V fs -1 : (20)
This simple exponential model may be a good approximation in the case of small surface perturbations [START_REF] Durand | Mass-velocity and size-velocity distributions of ejecta cloud from shock-loaded tin surface using atomistic simulations[END_REF][START_REF] Dimonte | Ejecta source model based on the nonlinear Richtmyer-Meshkov instability[END_REF] and for Taylor wave shocked tin [START_REF] Monfared | Experimental observations on the links between surface perturbation parameters and shock-induced mass ejection[END_REF]. M s 0.5 mg∕cm 2 , β 8.5, V fs 2000 m∕s, and V max 3600 m∕s. The ejected material is tin, whose optical index at λ 1.55 μm is n 3.31 8.67i [START_REF] Palik | Handbook of Optical Constants of Solids[END_REF] and bulk density ρ p 7.3 × 10 3 kg∕m 3 . The size distribution is f d ∝ d -5.2 p [START_REF] Sorenson | Ejecta particle size distributions for shock loaded Sn and Al metals[END_REF] with assumed minimum and maximum diameters d min 1 μm and d max 10 μm. Mie scattering phase functions used to propagate the photons in the medium are discretized between 0 and π on 1000 points. For each scattering event, the particle diameter is sampled using Eq. ( 6). The cloud is discretized in N 1600 domains (V max -V fs ) corresponding to a velocity step 1 m/s between neighboring particle slabs. The FOV of the PDV probe (collimated beam) is ω 250 μm. The initial length of the cloud is L 1 mm. At t 0, the shock releases at the metal-vacuum interface, giving rise to an expanding cloud of droplets. The polarization states of the incoming beam and of the local oscillator are linearly polarized along the Fig. 2. Noisy Doppler velocity spectrum of an ejecta cloud expanding in vacuum, reconstructed using Eqs. ( 11) and ( 12) with Φ mc ≃ hΦV i. Time step δt 625 ns. hB Φ i -10 dB.

x axis (see Fig. 1): S 0 I 0 ; Q 0 ; U 0 ; V 0 1; 1; 0; 0. A total number of 32 MC simulations is performed at 32 different times to reconstruct the time-velocity spectrograms near breakout. Between successive simulations, L is incremented by 1 mm (expansion of the cloud), which corresponds to a time step δt 625 ns. The position of the PDV probe is set arbitrarily below the critical distance d c , and simulations are performed before particles impact the probe. 15000 photons are propagated in the ejecta for each simulation (one slice of the Doppler spectrum shown in Fig. 2) and the boundaries of the cloud along x and y are set to 1.5 mm, much larger than the FOV. Figure 2 shows the velocity spectrograms. At each time of motion of the particle cloud, the average spectrum hΦV i is approximated by Φ mc V , determined using Doppler MC calculations. Finally, we replace hΦV i with Φ mc V in Eq. ( 19) to construct and represent the PDV spectrum.

B. Interpretation

A non-negligible PSD is detected below the free-surface velocity near breakout. As the cloud expands, this particular behavior progressively disappears. This feature is very likely to be explained by the multiple scattering of light in the ejecta, as also established in Ref. [START_REF] Andriyash | Optoheterodyne Doppler measurements of the ballistic expansion of the products of the shock wave-induced surface destruction: experiment and theory[END_REF], rather than by collisions among particles (two particles impacting would potentially drag before being swept up by the surface). Our results may explain the experimental observation presented in Ref. [START_REF] La Lone | Spall Strength and Ejecta Production of Gold Under Explosively Driven Shock Wave Compression[END_REF] on gold ejecta production (micro-jetting) under explosively driven shock wave compression, and may be supported by the experimental observations presented in Ref. [START_REF] Buttler | Second shock ejecta measurements with an explosively driven two-shockwave drive[END_REF] on tin plate ejecta experiments: as the areal mass of ejecta increases with shock breakout pressure, the amplitude below the free-surface velocity appears in the PDV spectrum. At later times, there is as much multiple scattering in the ejecta as just after shock breakout. However, high-order scattering events occur outside the FOV (Fig. 1), and photons cannot be collected due to the very limited aperture angle of the probe (≃0.1°).

Figure 3 shows the different contributions to the Doppler spectrum. As the areal mass is small (0.5 mg∕cm 2 ), a semitransparent cloud is observed. Unscattered photons can go through the cloud before being reflected by the free surface (red curve). The single scattering (SS) approximation is invalid below V ≃ 2800 m∕s, even for a very moderate areal mass. The number of scatterings depicted in Fig. 3 does not exceed 9 (for visibility), and above ∼15th order in our calculation, the contribution to the Doppler spectrum is negligible due a predominant absorption (the photon's weight approaches 0). A strong amplitude in the spectrum at V V fs is obtained for scattering orders between 1 and ∼5. Since particle sizes d p ∼ λ, photons are preferentially scattered in the forward direction and, because they do not acquire any Doppler shift (u i1 -u i ∼ 0), they do not contribute to the creation of the time-velocity spectrogram. Consequently, even after few scatterings, the largest contribution to the Doppler shift accumulated by the photon may come from the reflection on the free surface. This explanation may also justify why, in many ejecta experiments, a dark region between the ejecta cloud and the free-surface velocity is observed on PDV time-velocity spectrograms [START_REF] Buttler | Second shock ejecta measurements with an explosively driven two-shockwave drive[END_REF]. Figure 4 illustrates how the geometrical properties of the probe influence the contribution of multiple scattering to the Doppler spectrum. At a given time of motion of the cloud, this contribution increases with the beam diameter 2ω, showing that the whole geometry of the system must be taken into account to interpret precisely PDV data.

Figure 5 shows the distribution of the backscattered polarization states on the surface of the Poincaré sphere, after three (top) and five (bottom) scattering events undergone by the photons inside the FOV. As expected, we observe the progressive randomization of the backscattered polarizations (for an incoming linearly polarized light along the x axis). As the beam goes deeper in the ejecta and the scattering order increases, the backscattered light gets progressively depolarized. Fig. 3. Doppler velocity spectrum Φ mc at time t 3.12 μs (corresponding to the fifth slice shown in Fig. 2). SS is limited to the high-velocity region, and the non-zero power spectral density (PSD) detected below the free surface is justified by high-order scattering.

DISCUSSION

The only work dealing with multiple light scattering in metallic ejecta clouds is presented in Ref. [START_REF] Andriyash | Optoheterodyne Doppler measurements of the ballistic expansion of the products of the shock wave-induced surface destruction: experiment and theory[END_REF]. The authors apply a twoflux Kubelka-Munk solution [START_REF] Ishimaru | Wave Propagation and Scattering in Random Media[END_REF] of the radiative transfer equation (RTE) to estimate the Doppler spectrum and explain experimental PDV data. A two-flux theory supposes that the incident light is directly transformed into a diffused flux as it enters the ejecta. This may not be a realistic assumption in the case of a semi-transparent cloud, since the fastest region is very dilute, and the beam cannot be immediately diffused. Furthermore, a two-flux theory cannot be used for collimated light [START_REF] Ishimaru | Wave Propagation and Scattering in Random Media[END_REF]. They use a scattering coefficient σ, which is the same for the backward and forward propagating diffused fluxes. This may suppress the strong angle dependence of scattering and induce an overestimation of the single scattering contribution to the Doppler spectrum, since the probability of backscattering is actually much smaller than the probability of forward scattering for d p ∼ λ. Since a clear transition between single and multiple scattering regimes appears as the beam penetrates deeper into the ejecta, MC simulations or a complete solution of the time-dependent RTE are necessary to solve the problem (the polydispersity may be introduced with difficulty in the RTE). At later times, high-order multiple scattering may have a negligible impact on the Doppler signature and may be neglected to speed up the calculation and facilitate a possible inverse approach of determining ejecta properties using noninvasive PDV diagnostic. If the areal mass of ejecta is known, MC simulations can be performed to find out which particle size distribution is the most likely to explain the data. When a larger areal mass of particles is ejected (typically 5 to ∼15 mg∕cm 2 ), light scattering in the ejecta may be dependent, i.e., correlations between particles in the medium must be taken into account. In this regime, the relationship between the extinction coefficient of light and the particle volume fractions is non-linear (as opposed to independent scattering). Preliminary work on the interaction between a near infrared beam and a dense distribution (up to volume fractions 20%) of micrometer-sized metallic particles shows that the second-order Keller approximation [START_REF] Ishimaru | Attenuation constant of a coherent field in a dense distribution of particles[END_REF] to the extinction coefficient may well describe this dependent scattering regime. Such macroscopic variables may be integrated in the MC calculation to study the interaction of a PDV beam with a dense particle distribution, a situation that may be encountered when the ejecta are dragging in a gas.

CONCLUSION

In conclusion, this article addresses the multiple scattering of light in an expanding cloud. Using Doppler MC calculations incorporating the transport of polarization, we explain particular features observed in PDV Doppler spectra measured in shock compression ejecta experiments, and we show that the whole geometry of the system studied must be taken into account to interpret precisely PDV data. Extensive MC simulations are being carried out to improve our current understanding of ejecta physics probed by PDV and to reveal possibly hidden information encoded in the spectrum. For particle sizes and PDV wavelength on the same order of magnitude, multiple light scattering in the ejecta must be taken into account, even for moderate areal mass densities.

Fig. 1 .

 1 Fig. 1. Photon path in the ejecta at two different instants. Velocities are colinear to the z axis. At late time (right), the impact of multiple scattering on the measurement is reduced due to the limited FOV of the probe lens.

  aθ

Fig. 4 .

 4 Fig. 4. Doppler velocity spectra Φ mc calculated for a 1-mm-thick tin cloud (using 10000 photons) for three different beam waists ω. The properties of the cloud are described in Section A.

Fig. 5 .

 5 Fig. 5. Distribution of the backscattered polarization states (blue) on the surface of the Poincaré sphere after three (top) and five (bottom) scatterings undergone by the photons inside the FOV; red dot point: incoming polarization; RC: right-handed circular polarization; LC: left-handed circular polarization; H: horizontally polarized light (along x); V: vertically polarized light (along y); 45°: linearly polarized light at 45°with respect to x (see Fig. 1). The properties of the cloud are defined in Section. A.


Iréne Maillet
email: irene.maillet@intradef.gouv.fr

Laurent Michel
email: laurent.michel@isae.fr

Frédéric Souric

Yves Gourinat

Mode II fatigue delamination growth characterization of a carbon/epoxy laminate at high frequency under vibration loading

Keywords:

A fatigue crack propagation test under mode II has been set up for unidirectional composite laminate characterization at high frequency. The specimen is clamped on a shaker and loaded on its first bending resonance mode. This dynamical test allows an accurate delam ination extension monitoring based on the resonance frequency shift. Dynamical tests are conducted at 260 Hz and 400 Hz and their results are compared to those obtained with 10 Hz 3 ENF tests. Fatigue crack propagation results are presented taking into considera tion the loading ratio difference between both test configurations. All results are collapsed in a unique propagation rate curve. imen compliance. Fatigue behavior of carbon/epoxy laminate in pure mode II up to 10 6 cycles has been widely studied

Introduction

During their service life aeronautical structures are submitted to vibrations due to turbulent aerodynamic flow. These fluid structure interactions induce not only sonic acoustic fatigue spectra with low amplitude loads at high frequency range but also buffeting or flutter type fatigue with high loads at lower frequency range [START_REF] Dhainaut | Nonlinear random response of panels in an elevated thermal-acoustic environment[END_REF][START_REF] Kouchakzadeh | Panel flutter analysis of general laminated composite plates[END_REF]. At high frequency a very high num ber of cycles can be reached. Today certification of composite structures under damage tolerance is generally based on the concept of no growth damage [START_REF] Amc- | Proof of structure -fatigue and damage tolerance in AMC 20-29 composite aircraft structure[END_REF]. So whatever the load level and the frequency range, it must be proved that no delamina tion defect propagates under vibration loading. This raises the questions of frequency effects on crack propagation and of methods enabling the investigation of crack propagation under very high cycle fatigue (VHCF).

Delamination will propagate under one of the three following pure modes or a combination of them: the opening mode (Mode I), the sliding mode (Mode II) and the tearing mode (Mode III). Tests for pure mode I and II are largely used to characterize composite material under fatigue loading [4 6]. For the most critical mode, mode I, ASTM has defined the Double Cantilever Beam (DCB) test as a standard [START_REF]Standard test method for mode I fatigue delamination growth onset of unidirectional fiber-reinforced polymer matrix composite[END_REF]. Currently no standard has been officially released for mode II fatigue test. Several test configurations can be used for mode II: the three point bending end notch flexure (ENF), the four point bending end notch flexure (4ENF) and the end loaded split (ELS) [START_REF] Blackman | Mode II fracture testing of composites: a new look at an old problem[END_REF]. Mode II fatigue tests on fiber reinforced polymers (FRP) are commonly carried out at a frequency up to 10 Hz, under controlled loading ratio and with a test stop criterion related to spec [9,4,6,10 12]. At this frequency of 10 Hz, more than 11 days are needed to reach 10 6 cycles. In these conditions the explo ration of high cycle fatigue ð10 6 10 9 Þ is greatly time consuming. One solution to reduce the test time is to increase the test frequency. However, this question doesn't seem to be addressed in the literature for delamination propagation under mode II.

Frequency effects on fatigue life have been widely investigated in the literature on polymers and composite materials [ 13 15]. For composites most of these studies concern fatigue tests performed at a rather low frequency range from 0.01 Hz to 30 Hz on tension/compression specimens [16,14,17 19]. They show a large dependence of the fatigue behavior of FRP on loading frequency. This is due to the inherent viscoelasticity of organic matrix. Indeed there are two notable fre quency effects which oppose each other: creep and hysteretic heating [START_REF] Hahn | The effect of loading parameters on fatigue of composite laminates: part IV information systems[END_REF] and it reduces fatigue resistance with a decrease in frequency. At higher frequencies, hysteretic behavior leads to a self heating which may turn into a thermal failure and result in a reduction of fatigue life of the component. Very significant temperature rises have been monitored for the highest frequency tests (10 30 Hz). Xiao [START_REF] Xiao | Modeling of load frequency effect on fatigue life of thermoplastic composites[END_REF] studied angle ply APC2 laminates and computed a temperature increase of up to 135 °C. Barron [17] also reported temperature rises up to 90 °C for a carbon fiber/epoxy composite at 20 Hz. This is particularly true for highly hysteretic or viscoelastic materials at the test temperature and also for disoriented lay ups towards the loading direction [START_REF] Barron | Frequency effects on the fatigue behaviour on carbon fibre reinforced polymer laminates[END_REF][START_REF] Xiao | Modeling of load frequency effect on fatigue life of thermoplastic composites[END_REF][START_REF] Minak | On the determination of the fatigue life of laminated graphite-epoxy composite by means of temperature measurement[END_REF]. Rotem [START_REF] Rotem | Effect of the frequency of the alternating load on fatigue delamination of laminates[END_REF] tested a quasi isotropic graphite/epoxy laminate under reverse loading. He observed that between 2.8 Hz and 10 Hz, when the frequency increases, the heat generation increases too. High temperatures lead to earlier crack initiation in fatigue. Temperature increase due to self heating might also be used to monitor the fatigue life of laminated composites [START_REF] Minak | On the determination of the fatigue life of laminated graphite-epoxy composite by means of temperature measurement[END_REF]. However it is also reported that at low stress levels this self heating is reduced and disables any thermal failure. Indeed stress level is one of the key factors [START_REF] Herztberg | Fatigue of engineering plastics[END_REF][START_REF] Kharrazi | Frequency-dependent fatigue damage accumulation in fiber-reinforced plastics[END_REF][START_REF] Minak | On the determination of the fatigue life of laminated graphite-epoxy composite by means of temperature measurement[END_REF] for frequency effects. Some experimentations at much higher frequencies were already achieved in order to accelerate fatigue tests and to reach a very large number of cycles [ 21 23]. For a frequency of 100 Hz, several authors [START_REF] Michel | Fatigue strength of carbon fibre composites up to the gigacycle regime (gygacycle-composites)[END_REF][START_REF] Hosoi | High-cycle fatigue characteristics of quasi-isotropic CFRP laminates over 10 8 cycles (initiation and propagation of delamination considering interaction with transverse cracks)[END_REF] reported limited effects on fatigue life of composite laminates. Wang et al. [START_REF] Wang | Fatigue crack growth of bonded composite repairs in gigacycle regime[END_REF] have even developed an experiment at 20 kHz to perform fatigue tests on epoxy bonded joints. These experiments were possible because specific specimen dimensions and the limited strain amplitudes have restricted any inconvenient self heating effects.

Similarly the geometry of notched specimen used for crack propagation is likely to reduce self heating which is largely dependent upon the capacity of samples to dissipate heat. It is related to the thermal transfer characteristics of the material, to extension of the heated zone and to the specimen surface area to volume ratio [START_REF] Herztberg | Fatigue of engineering plastics[END_REF]. Indeed for notched specimen, the zone of high stresses is localized and the rest of the specimen is able to dissipate enough heat. This limits any large temper ature rise and any thermal failure.

Frequency effects on fatigue crack propagation in polymers are largely dependent upon the polymer itself [13,24 26]. Yuen et al. [START_REF] Yuen | The effects of loading frequency, tensile overload and compressive underload on the fatigue crack propagation behaviour of polymethyl methacrylate[END_REF] reported no effect of frequency at room temperature for the PMMA. Merah et al. [START_REF] Merah | Modeling the combined effects of temperature and frequency on fatigue crack growth of chlorinated polyvinyl chloride (CPCV)[END_REF] achieved tests on Chlorinated PolyVinyl Chloride (CPVC) to find the influence of temperature and frequency between 0.1 Hz and 10 Hz. The frequency sensitivity increases with an increase of the fixed test temperature. It is linked by the authors to a major crazing fatigue mechanism. In a study performed on several polymers, Ramsteiner [START_REF] Ramsteiner | Fatigue crack growth in polymers[END_REF] displayed da=dN f ðMKÞ curves for 1 and 10 Hz which are almost identical. However the author reported that crack propagation is faster at 10 Hz than at 1 Hz. It is difficult to get a clear comprehension of frequency effects from these studies.

There are very few studies concerning frequency effects on crack propagation in composites. Nakai et al. [START_REF] Nakai | Effects of loading frequency and environment on delamination fatigue crack growth of CFRP[END_REF] performed fatigue crack propagation tests on CF/epoxy laminates at different frequencies and environment. The authors reported an effect of cycle dependance for mode I test and of time dependance for mode II test. Studies on temperature effects are more common. For example Sjorgen et al. [START_REF] Sjogren | Effects of temperature on delamination growth in a carbon/epoxy composite under fatigue loading[END_REF] investigated the temperature effect on the delamination growth of carbon/epoxy composite in mode II. Fatigue propagation rates at 100 °C at same energy release rate (ERR) are faster than at room temper ature and a decrease of 30% of the critical value is observed in static loading.

In the field of composite delamination propagation, researchers limited their investigation to frequencies lower than or equal to 10 Hz. So the first aim of this work is to develop a simple test approach to perform fatigue crack propagation (FCP) under mode II at higher frequencies. As performed before for mode I [START_REF] Maillet | A new test methodology based on structural resonance for mode I fatigue delamination growth in unidirectional composite[END_REF], fatigue test of cracked specimens is performed under vibration loading at resonance. An ELS type test specimen is loaded to its first bending resonance mode. The loading ratio is then R 1. When delamination propagates, resonance frequency decreases [START_REF] Saravanos | Effects of delaminations on the damped dynamic characteristics of composite laminates: analysis and experiments[END_REF]. So the resonance frequency is used to monitor the crack propagation. Dynamic tests are achieved for two frequencies 260 Hz and 400 Hz. Results of these tests are compared with more classical ones which are performed with the ENF (End Notched Flexure) configuration on a servo hydraulic machine at 10 Hz and R 0:1. Effects of frequency on the behavior of a carbon fiber/epoxy composite are then investigated and discussed together with the possibility to accelerate FCP tests.

Material and specimens

The material used in this study was the Hexply Ò M21 resin and T700 Carbon fiber unidirectional prepreg supplied by Hexcel Ó . Beam specimens of 180 Â 25 Â 5 mm 3 were cut out from unidirectional ½0 10 sym plates cured in a thermal press. The curing procedure followed a standard cycle with a homogenization stage at 150 °C and a curing one at 180 °C as recom mended by the supplier. A 25 lm thickness Teflon film was inserted at mid plane of laminates during stacking operations to create an artificial starter crack (Fig. 1). All specimens were then pre cracked to get at least a 2 mm crack extension before fatigue tests. Mechanical properties of the material were also measured and main results are listed in Table 1.

Low frequency fatigue tests

Test configuration

Low frequency fatigue tests have been performed with an ENF configuration on a servo hydraulic fatigue machine with a 5 kN load cell capacity, at a controlled displacement ratio of R 0:1 and at a frequency of 10 Hz. The dimensions of the ENF test device used in this study are given in Fig. 2. All specimens had an initial crack length ða 0 Þ of 15 mm before pre cracking. dynamical fatigue at resonance should largely accelerate test durations. However, further tests have to be completed and performed for other load ratios and other materials in order to support this eventuality.

  . At low frequencies, creep behavior is predominant

	Nomenclature
	DG	G max G min
	DG eq	equivalent delta G
	da dN c	crack propagation rate Hojo formula material constant
	m lt m tt	Poisson ratio transverse Poisson ratio
	q	density
	u	phase between input acceleration and specimen displacement
	a	current delamination length
	a 0	initial delamination length
	a i	specimen delamination size
	b	specimen width
	C	compliance
	C 0	compliance law parameter
	d imp	imposed displacement
	dep numerical numerical model displacement
	dep specimen
		specimen test displacement
	E f	flexural Young modulus
	E numerical numerical model bending modulus
	E specimen fitted specimen bending modulus
	E t	transversal Young modulus
	ERR	energy release rate
	f f	final frequency of the specimen
	f numerical numerical model frequency
	f sweep	specimen resonance frequency obtained by frequency sweep
	G	ERR
	G c	critical ERR
	G final	dynamical test ERR
	G lt	shear modulus
	G max	maximum ERR
	G numerical numerical model ERR
	K	stress intensity factor
	k	Paris law parameter
	k D	DG eq Paris law parameter
	L	distance between ENF rollers
	L 0	free specimen length
	L e	specimen length
	m	compliance law parameter
	m	G N curve parameter
	m D	DG eq N curve parameter
	N	number of cycle
	n	Paris law parameter
	N c	critical number of cycles
	n D	DG eq Paris law parameter
	P	load
	p	GN curve parameter
	p D	DG eq N curve parameter
	R	fatigue loading rate