Eduardo Cuenca

Mathieu Ducousso

Alexandre Rondepierre

Laurent Videau

Nicolas Cuvillier

Laurent Berthe

François Coulouvrat

Propagation of laser-generated shock waves in metals: 3D axisymmetric simulations compared to experiments

I. INTRODUCTION

Sound generation by electromagnetic wave is well known since G. Bell first experiments in 1880 using rapidlyinterrupted beam of sunlight [START_REF] Bell | The production of sound by radiant energy[END_REF] . At first, this discovery lead to none application, but following the laser discovery by T. Maiman in 1960, acoustic generation could be obtained using laser pulses [START_REF] White | Generation of elastic waves by transient surface heating[END_REF] . Critical parameters of sound generation are the laser pulse duration and its energy.

When laser energy is under the ablation threshold, sound generation is in the thermoelastic regime. In that case, the most common way of sound generation results from the sudden thermal expansion due to the heating of a small volume near the surface of the illuminated material, controlled by the laser diameter and the optical penetration. The frequency content of the generated ultrasound is determined on first order by the laser pulse duration. Frequencies up to THz can be generated using picoseconds or shorter pulses laser pulses. Such field of laser ultrasonics, named picosecond ultrasonics, found application at the nanometer scale from fundamental physics [START_REF] Greener | Highfrequency elastic coupling at the interface of van der waals nanolayers imaged by picosecond ultrasonics[END_REF] and 4 , to material characterization [START_REF] Wang | Imaging grain microstructure in a model ceramic energy material with optically generated coherent acoustic phonons[END_REF] or biological evaluation [START_REF] Ducousso | Evaluation of mechanical properties of fixed bone cells with sub-micrometer thickness by picosecond ultrasonics[END_REF] . In picosecond ultrasonics, acoustic wavelength is comparable to optical penetration of the illuminated material. a) Electronic mail: eduardo.cuenca@safrangroup.com Thus the acoustic source must be considered over a volume under the surface. This phenomena induces, in the case of the thermoelastic regime, the generation of only longitudinal waves. Using nanosecond laser pulses allows to generate acoustic frequencies in the MHz range, hence wavelengths of a few tens of micron. Laser ultrasound lead to a variety of applications in Non Destructive Testing (NDT), from surface [START_REF] Shan | Surface-breaking fatigue crack detection using laser ultrasound[END_REF] to bulk [START_REF] Yan | Characterization of progressive fatigue damage in solid plates by laser ultrasonic monitoring of zero-group-velocity lamb modes[END_REF] or bounding [START_REF] Hodé | Nondestructive evaluation of structural adhesive bonding using the attenuation of zero-group-velocity lamb modes[END_REF] evaluation. In this configuration, the acoustic generation (in opaque materials with respect to the laser wavelength) can be simply considered as surface generation. Resulting generated acoustic waves are both longitudinal and transverse of comparable magnitudes, but the later predominates at the epicenter. The direction of propagation of the generated waves shows a significant angle, of few tens of degrees, with respect to the normal of the surface [START_REF] Scruby | EnglishLaser ultrasonics : techniques and applications[END_REF] .

As the laser energy is increased and exceeds the ablation threshold of the illuminated material [START_REF] Tsibidis | The influence of dynamical change of optical properties on the thermomechanical response and damage threshold of noble metals under femtosecond laser irradiation[END_REF] , a high pressure plasma (in the range of GPa) is created. In this case, the acoustic generation results from the resulting plasma expansion away from the illuminated surface. For metals, such ablation threshold is classically of a few MW /cm 2 , but laser intensities up to several GW /cm 2 can be used for shock wave emission [START_REF] Berthe | Shock waves from a water-confined laser-generated plasma[END_REF] . Transition from thermoelastic to ablative regimes does not significantly affect the acoustic spectrum, with a first order dependence with laser pulse duration. However, nature and directivity of radiated waves are significantly changed. Indeed, while both longitudinal and transverse waves are observed in oblique directions in the thermoelastic regime(for ns laser pulses), mostly a longitudinal wave in the direction normal to the surface is radiated in the ablative regime [START_REF] Scruby | EnglishLaser ultrasonics : techniques and applications[END_REF] . The wave amplitude is directly correlated to the laser intensity, from a few MPa for intensities just above the ablation threshold, to more than 150 GPa for laser intensities up to 10 3 GW /cm 2 13 . The induced acoustic wave goes from a smooth linear elastic wave, to a nonlinear elastic wave, to a non-linear elastoplastic one, possibly inducing material phase-changes.

The correct modelling of the induced transient phenomena obviously depends on the wave amplitude. It ranges from analytical expressions for the lowest amplitudes [START_REF] Mi | Parametric studies of laser generated ultrasonic signals in ablative regime: Time and frequency domains[END_REF] to ab initio calculation for the more intense ones [START_REF] Mundy | Ultrafast transformation of graphite to diamond: An ab initio study of graphite under shock compression[END_REF] . Between these extrema, GPa laser-generated shock waves are classically described using a hydrodynamic numerical solver, including adapted equations of states relating energy and stress into the shock [START_REF] Ramis | Multi-fs -a computer code for laser-plasma interaction in the femtosecond regime[END_REF] . However, such simulations are very demanding in terms of numerical resources, and present several restrictions. Thus, despite promising applications of laser-generated shock waves for some engineering applications 171819 , simulations in 2D or 3D dimensions remain a major lock, particularly for applications on thick and/or complex structures such as aeronautical composite materials. Moreover, the equations of state related to the laser shock process are not well established since they depend, for instance, on the properties of the material, the induced pressure and the duration of the loading. This regime of shock wave, up to a few GPa, is difficult to simulate, and experimental verifications are systematically required. However, such shock waves remain weak and correspond to an acoustical Mach number M of around 10 -2 only. Thus, neglecting the surface residual stress induced by the laser pulse, linear acoustical modelling should be efficient to simulate the propagation of the wave in the bulk of the material while the acoustic wave generation at the illuminated material surface has to be considered as a fundamentally nonlinear process. This hybrid approach, separating wave nonlinear generation and linear propagation, is the purpose of the present work. This work discusses the validity of using linear elastodynamic simulations in the case of shock wave propagation with laser intensities up to the air breakdown threshold, i.e. for simulating simple table-top experiments without requirements to vacuum condition. Water [START_REF] Berthe | Shock waves from a water-confined laser-generated plasma[END_REF] or transparent solid materials [START_REF] Bras | Laser shock peening: Toward the use of pliable solid polymers for confinement[END_REF] can be used to increase the amplitude of the shock wave with respect to air condition. Stresses up to several GPa are reached, without significant modification of the directivity diagram of the wave generation : a longitudinal shock wave which propagates through the normal direction of the surface towards the bulk of the material. However, a diffractive optical component is used in experiments to obtain a laser top-hat profile on the material surface. We will demonstrate that this discontinuity at the edges of the laser illumination creates, by diffraction, shear waves of significant amplitude. We will show such diffraction process can be described correctly only using a 3D approach. We will also demonstrate the validity of our proposed hybrid simulation approach to simulate laser-generated shock waves (of magnitude up to a few GPa) in both aluminum and titanium materials from systematic comparison to experiments for several laser diameters, sample thicknesses and laser energies. Excellent agreement will be systematically obtained. Such work will pave the way towards simple and efficient simulations of laser-generated shock wave applications, for NDT applications such as evaluation of mechanical strength of structural bonding, [START_REF] Ducousso | Quantitative evaluation of the mechanical strength of titanium/composite bonding using laser-generated shock waves[END_REF] or delamination generation in composite materials [START_REF] Ghrib | Generation of controlled delaminations in composites using symmetrical laser shock configuration[END_REF] .

The second part of the article introduces laser-generated shock waves principles and the laser facility used for the experimental verification. The third section is devoted to numerical simulations. First, we use the 1D hydrodynamic code named Esther to quantify (in amplitude and in time profile) the generated surface pressure from laser illumination. This will play the role of an acoustic source in the hybrid approach. Such pressure is then propagated into an isotropic sample using a 3D model of linear elastodynamics, assuming asymmetry for the sake of simplicity. Virieux numerical scheme [START_REF] Virieux | P-sv wave propagation in heterogeneous media: Velocity-stress finite-difference method[END_REF] is chosen, and its numerical parameters (grid size 0.005 mm, compared to the wavelength 0.1 mm for a 15 ns period, Courant-Friedrichs-Lewy condition [START_REF] Courant | On the partial difference equations of mathematical physics[END_REF] , CFL = 0.2) are determined to minimize numerical dispersion and error. In a fourth part, a detailed discussion of the induced waves from a top-hat laser illumination is proposed, outlining the importance of diffraction to explain all observed arrivals on the opposite face. An extensive comparison between simulations and experiments is finally performed for several experimental parameters (laser intensity, laser diameter, material, sample thickness), demonstrating the validity of the present hybrid approach to simulate wave phenomena associated to laser-generated shock waves in such regime. The analysis is completed by a discussion about the possible nonlinear mechanisms that could limit the model validity.

II. GENERAL PRINCIPLES AND EXPERIMENTAL SETUP

A. Principle of shockwave laser generation Fig. 1 presents a general view of experimental configuration. An intense laser is focused over the surface of the target material. It ablates its surface, over a thickness typically of 1.7-3 µm 23 , changing superficially the material into a dense plasma (Fig. 1). A transparent dielectric material, such as water or polymer, is generally used to confine the expansion and therefore increase the surface pressure level induced by the plasma expansion compared to the direct regime 1220 . As a result of the recoil moment of the ablated momentum, a pulsed shockwave is generated in the bulk of the target. It propagates through the target and emerges at its rear free surface. Using diagnostic like VISAR (Velocity Interferometer System for Any Reflector) (see subsection 2.3), measured velocity provides a signature of laser shock travel inside the material (typical rear free surface velocity temporal profiles will be shown later).

FIG. 1. Scheme of experimental set-up for laser shock generation and rear free surface velocity measurements. The focused laser, confined by water, illuminates the upper face of the target plate. In the opposite face, an optical detection, VISAR, records the free surface velocity.

B. Hephaïstos laser facility

The Laser Hephaïstos Facility (PIMM Laboratory for Processes and Engineering in Materials and Mechanics, ENSAM ParisTech and CNRS) is used to generate laser shots [START_REF] Bras | Laser shock peening: Toward the use of pliable solid polymers for confinement[END_REF] . Two Gaïa HP lasers from THALES (Elancourt, France) emits two Nd:YAG synchronized or delayed lasers pulses, both at 532 nm, with a repetition frequency of 2 Hz. Through a Gaussian temporal shape of 7 ns of duration (Fig. 2(a)), it is able to deliver an energy up to 14 J when both laser pulses are superposed. Each laser beam can be focused through an optical lens (Fig. 1) controlling the focal diameter as desired, here between 3 and 5 mm, as is represented in Fig. 2(b). This image has been recorded using a CCD (Charge-Coupled Device) camera located at the focus of the laser. Pixel resolution is 5.5 µm. A Diffractive Optical Element (DOE), placed after the focusing lens, was used to obtain the top-hat spatial profile. The shape of the focused laser appears almost perfectly circular with a homogeneous intensity profile measured with a standard deviation of only 0.023.

C. Rear free surface velocity measurement

The velocity of the plate opposite free surface is measured by a non-contact optical diagnostic tool, so called VISAR [START_REF] Barker | Laser interferometer for measuring high velocities of any reflecting surface[END_REF] . This Michelson-type interferometer allows to measure the Doppler shift in the wavelength of a probe laser of 532 nm wavelength and 5W power (Coherent Company, Santa Clara, CA, USA) reflected on the target free surface moving due to the shockwave. Such device has already been used for laser shockwave detection, both in direct [START_REF] Tollier | Study of the laser-driven spallation process by the velocity interferometer system for any reflector interferometry technique. i. laser-shock characterization[END_REF] and water-confined [START_REF] Berthe | Shock waves from a water-confined laser-generated plasma[END_REF] regimes. This kind of detection yields a time resolution of the opposite surface velocity of 1 ns.

III. NUMERICAL SIMULATION A. Shockwave generation

To get plasma pressure profile, a 1D Lagrangian code describing both laser-matter interaction and shock propagation is used (so called Esther [START_REF] Bardy | Numerical study of laser ablation on aluminum for shock-wave applications: development of a suitable model by comparison with recent experiments[END_REF] ). It simulates the deposition of laser energy and the evolution of target material state from solid to plasma. The code also includes radiative transfer and heat conduction. Mechanical models reproduce elastoplasticity of material [START_REF] Steinberg | A constitutive model for metals applicable at high-strain rate[END_REF] up to damaging and fracturation [START_REF] Johnson | Dynamic fracture and spallation in ductile solids[END_REF] . In this work, Esther allows to simulate the pressure at the surface, also called ablation pressure, resulting from a laser illumination of a defined duration and intensity. One example of temporal profile is shown in Fig. 2(a) for an intensity of 26.7 J/cm 2 in comparison with laser beam time profile. From that simulation, it is possible to evaluate the rise time and the Full Width at Half Maximum (FWHM) of the surface pressure, here observed at 8 and 15 ns, respectively. Beyond the peak the pressure decreases slowly, leading to a pulse duration, T pulse , of typically 0.1 µs and therefore a frequency spectrum up to 10 MHz.

Esther ablation pressure time profiles for laser intensity ranging from 6.7 J/cm 2 to 40 J/cm 2 in aluminum are shown on Fig. 3(a), with an inset in the right showing the shock rise and peak during the first 40 ns. The waveform, normalized by its peak value, little varies with laser intensity, except a small widening : the shock arrives slightly sooner and the pressure decays a little bit slower for higher intensities. Peak overpressure, displayed on Fig. 3(b), follows a square root dependency with intensity in agreement with the analytical model 29 and [START_REF] Fabbro | Physical study of laser-produced plasma in confined geometry[END_REF] :

P = 0.10 α 2α + 3 1/2 Z 1/2 I 1/2 0 ; ( 1
)
where P is the peak pressure given in [kbar], α is a physical parameter, here chosen equal to 0.6, calibrated by experiments. Impedance Z is given by 2/Z = 1/Z 1 + 1/Z 2 where Z 1 and Z 2 are the impedances in [g/cm 2 s] of water and metal.

The laser surface intensity is

I 0 in [GW /cm 2 ]
. This analytical model is well known and confirmed by other works such as [START_REF] Correa | Random-type scanning patterns in laser shock peening without absorbing coating in 2024-t351 al alloy: A solution to reduce residual stress anisotropy[END_REF] . From such data we are able to generate a shockwave source term. Spatially, we impose a uniform normal stress over the laser spot surface, exponentially matched to zero to nonablated area, according to the spatial intensity profile shown on Fig. 2. The temporal waveform is determined by the Esther simulation. This source term is applied as a boundary condition. It, obviously, induces severe gradients of stress, both temporally, in the form of a sharp shockwave, and spatially, with a sharp pressure decrease at the edges of the laser focal spot. Longitudinal (compression) and transverse (shear) waves will thus be generated, respectively from the central and the outer part of the illuminated zone. The longitudinal release wave created just behind the compression wave has been deeply studied and understood by the literature 123220 . However, focus is generally on first longitudinal arrival, occurring much before the shear wave arrival which has drawn little attention. Thermo-elastic regime, related to low laser energies, induces shear waves of same magnitude as the longitudinal waves. For the ablation regime considered here, the shear waves due to the ablation process are frequently assumed almost negligible, at least for the few existing analytical solutions [START_REF] Miller | The field and radiation impedance of mechanical radiators on the free surface of a semi-infinite isotropic solid[END_REF] . However, a realistic source of finite size has not been considered yet to our knowledge. With this in view, contribution of such shear waves is integrated here over 2π radians revolution, in order to properly describe laser shock Hephaïstos experiments.

B. Numerical model

To model the shock propagation, we rely on the Virieux Finite-Difference Time-Domain (FDTD) method [START_REF] Virieux | P-sv wave propagation in heterogeneous media: Velocity-stress finite-difference method[END_REF] . Such numerical scheme has been purposely designed to discretize the elastodynamics equations by finite differences, on a grid staggered both in time and space (see Appendix B). The advantages of this approach are: (1) stability whatever Poisson's ratio, (2) small numerical dispersion and anisotropy and (3) easiness to implement sources and boundary conditions. Staggered grid in time, between velocities on one side and stresses on the other side, induces a globally second order precision [START_REF] Virieux | Dynamic faulting studied by a finite difference method[END_REF] . Fourth order spatial finite differences, following [START_REF] Levander | Fourth-order finite-difference p-sv seismograms[END_REF] , are here chosen to minimize numerical dispersion. Both two and three dimensional versions of the algorithm, the 3D one with axial symmetry, are implemented (see Appendix A).

For modelling semi-infinite plates, the chosen boundary conditions are: (1) imposed normal stress according to Fig. 2(a The temporal and spatial steps must be chosen small enough regarding the impulsive behavior of the input shockwave. In our case, the shock duration at half amplitude is 15 ns which corresponds to a longitudinal wavelength of about 0.1 mm for both aluminum and titanium. Therefore spatial mesh, ∆x, has to be a fraction of this value. We explore the range one fifth (0.02 mm) [START_REF] Levander | Fourth-order finite-difference p-sv seismograms[END_REF] , one tenth (0.01 mm) [START_REF] Virieux | P-sv wave propagation in heterogeneous media: Velocity-stress finite-difference method[END_REF] , one twentieth (0.005 mm) and one fortieth (0.0025 mm) to check numerical convergence illustrated by Fig. 4 (b), and zooms in Figs. 4 (c) and (d). Time step, ∆t, has to satisfy the CFL stability condition, given by [START_REF] Levander | Fourth-order finite-difference p-sv seismograms[END_REF] :

∆t ≤ 0.606 ∆x c L ; (2)
Both maximum value and half of CFL constant have been tested on Fig. 4 (a).

The materials investigated in this article are aluminum and titanium. Both material are anisotropic at the grain scale thanks to their crystal structure, Face-Centered Cubic (FCC) in the case of aluminum and Hexagonal Closed Packed (HCP) in the case of the titanium. However, thanks to the random orientation and shape of grains, they can nevertheless be treated at a large scale, as a quasi-isotropic materials [START_REF] Brenner | Elastic anisotropy and yield surface estimates of polycrystals[END_REF] . In a first approach, using a purely elastic code, we have to check that the laser excitation does not reach the elastic limit, recalled in Table 1, along with the elastic properties of the two metals [START_REF] Committee | ASM Handbook: Properties and Selection: Nonferrous Alloys and Special-Purpose Materials (ASM Handbook)[END_REF] . Fig. 4 (a) presents influence of CFL for ∆x=0.005 mm and Fig. 4 (b),(c) and (d) show the convergence with mesh refinement for CFL = 0.303. In all cases titanium has been illuminated by a 26.7 J/cm 2 laser pulse with a focal spot diameter of 3 mm. The displayed quantity is the normal rear free surface velocity at the opposite face (2 mm thickness) at the center of the focal spot in the 3D axisymmetric case. The observed signals are successively the first longitudinal shock (about t = 0.33 µs), followed by the first release wave (t = 0.4 µs). This pattern is repeated every 0.66 µs for each echo. In between a smooth shear wave appears with first arrival at about t = 0.78 µs in agreement with shear wave velocity (c T =3100 m/s). At t = 1.3 µs we observe an arrival linked to longitudinal to shear conversion at the free surfaces. Large values of CFL (equal to maximum value 0.606 for stability, or three fourth of it) lead to dispersion mostly before the shock. Lower values (0.303 or 0.152) induce dispersion after the shock, all the more intense as CFL is small. Hence, the value CFL= 0.303 appears a good compromise with no visible dispersion of the shock and small dispersion after it. It is chosen for the following simulations. Strong dispersion is visible after the shocks for the two largest grid meshes (∆x = 0.02 and 0.01 mm ) and increases as expected with time. Smaller grids lead to a moderate dispersion. To compromise with the increasing computational time (see Table 2), ∆x=0.005 is finally selected. Note on the various zooms that both the release and the shear arrivals are of lower frequency content than the longitudinal shock, and are therefore much less affected by numerical dispersion.

Density [kg/m 3 ] Young modulus [GPa] Poisson's ratio HEL [GPa] c L [m/s] c T [m/s]
In Fig. 5, a 2D simulation is superposed to a VISAR measurement in the same conditions as previously described, expect the laser energy, now 6.7 J/cm 2 . The data show the same succession of arrivals (first shock, first release wave, first shear wave, second shock, etc. . . ). Compared to axisymmetric simulations of Fig. 5(b), it is obvious that 2D simulations produce only tiny release and shear waves (dotted ellipses). Indeed these waves emanate from the circular edge of the focal spot and arrive in phase at the measurement point so that they cumulate. This phenomenon cannot be described using the 2D approach where only two points represent the edge of the focal illumination. In the following only 3D axisymmetric simulations will be considered.

IV. ACOUSTICAL DESCRIPTION OF THE GENERATED WAVES

We now propose to describe the different waves generated during a laser-driven plasma expansion. The optical VISAR detection placed on the back face of the sample at the epicenter of the shock generation, can only give access to the back surface velocity. Such experimental detection is presented as a blue line in Fig. 6(b), compared to simulations in red. Data are obtained illuminating a 5 mm thick titanium plate with a 3 mm diameter focal spot while delivering an energy of 41 J/cm 2 . This 'thick' configuration allows a clear separation between the different waves. From such experimental signal, 3 different acoustic contributions can be isolated in a 3 µs time-window. First, two longitudinal pressure shockwaves, around 0.8 µs and 2.5 µs, are easily identifiable. Second, just after each longitudinal shock wave, a release wave (with negative velocities) is observed. Finally, at 1.6 µs, a third arrival is observed, to be later identified as a shear wave. The agreement of the arrival times of the different waves on Fig. 6(b) is remarkable. The two release waves are somewhat overestimated by the simulation, as is also the second shock. It could be due to some underestimation of the source or to absorption effects neglected in the simulation. The slight delay of the simulated second shock may also arise from neglected non-linearities. The shear wave turns out to be well described. Fig. 6(a) presents in color levels the computed velocity along the thickness of the plate at the center of it as a function of time. Shock waves appear in blue straight lines, as coming from the center of the laser spot. Release and shear waves are due to diffraction at the edges of the spot, therefore appear later on in time (time necessary to propagate from the edge to the center), initially as a curve in the time-space diagram. Shear wave is initially positive and then negative but after some propagation the negative part dominates. Mode conversions occurring at the reflection on the two free surfaces are visible, and labeled. For example: PS for a pressure shock wave reflecting as a shear wave or SP for a shear wave reflecting as a pressure shock wave. PSP is for both SP and PS waves reflection on the illuminated surface.

The origin of such waves can be better understood thanks to the time resolved simulations. Fig. 7 presents the same case but with 2D velocity maps at 4 successive times. The most recognizable signal in Fig. 7 is the pressure wave (P) generated by the plasma expansion, appearing as a quasi-plane waveform, with a geometric extent equal to the spot diameter at the beginning of the propagation, Fig. 7(a). Due to diffraction, a release longitudinal wave (R) is also generated all around the edge of the focal spot as result of the huge discontinuity of the laser illumination. This wave advances from the edge of the focal spot with the same c L velocity of the longitudinal pressure wave, thus creating a toroidal wave front tangential to the edges of the P wave, Fig. 7(b). A shear wave is also radiated from the discontinuity at the spot edge, but propagates slower. After some propagation time (Fig. 7(c)), the inner part of the release wave reaches the center, converging into a negative velocity release wave. Fig. 7(d) shows the similar but retarded convergence of the shear wave, also creating a negative velocity field.

The simulations also provide the stress value at each point and instant as shown in Fig. 8(a) for the normal stress component σ zz along the axis. The overall shape of the figure is very similar to the velocity field (Fig. 6(a)) except the change of sign of all waves after each reflection on the free stress surfaces. The time evolution of the same stress component at three different depths (Fig. 8b) illustrates this alternation of compressions and tractions. In particular, the red position outlines a constructive interference between the reflected shear shock, now a traction, and the direct shear wave, also a traction. This induces an increase of the traction of about 20% compared to other positions.

V. PARAMETRIC COMPARISON BETWEEN SIMULATIONS AND EXPERIMENTS

A. Sample analysis

Once the numerical model has been established and all phenomena have been explained, it is now possible to perform an extensive parametric survey by varying key parameters of the process and systematically comparing simulation outputs with Hephaïstos experimental VISAR data. To our knowledge such a large-scale parametric comparison and experimental validation has never been carried out. In a first stage, we vary the sample parameters, e.g. the structural material and the sample thickness. We compare the same configuration (6,7 J/cm 2 laser with 3 mm focal diameter for samples 1 and 2 mm thick) for both aluminum and titanium in Fig. 9. It is remarkable that, as both materials have very similar acoustic properties c L and c T , the measured velocity profiles are very similar. Agreement between simulations and experiments is remarkable, with each arrival (labeled as in section 4) finely predicted. Over the considered time window of 1.7 µs, the number of arrivals is much larger for "thin" (1 mm) than for "thick" (2 mm) samples: arrivals up to 9P/9R are detected in the first case, only arrivals 5P/5R in the second one. Only the direct shear wave arrivals (S) is visible but three mode conversions are detected (2PS, 4PS, 6PS) in the thinner case, one in the thicker case. Also clearly visible for both thicknesses is the time delay between the first compression shock P and the first release wave R. This delay is due to the longer time needed for the release wave to propagate from the edge of the focal spot, compared to the direct propagation time of the shock P. As a counterpart of this time separation, the re-lease wave is less intense in the 1 mm case. For other arrivals, the propagation time is dominated by the number of crossed thicknesses, and the influence of the initially longer path for the release wave gets negligible. Therefore 3R or 5R wave arrives just after 3P or 5P waves. This effect is also visible for the first arrivals in case of thicker samples (see 5 mm case, Fig. 6). Differences between simulations and measurement appear only for the amplitudes. For titanium, the agreement remains pretty good. The simulated shock amplitude is slightly underestimated, maybe due to some underestimation of the source term. However, when propagating back and forth, the shock amplitude does not change significantly. Absorption does not play a major role, may be due to the high elasticity limit of this material. On the contrary, for aluminum, the shock amplitude visibly decays, a phenomenon that cannot be captured by our simulations. As arrival times are always very consistent, this is most likely explained by absorption: nonlinearities would be accompanied by a time delay (shock wave would arrive earlier and earlier) that is not observed experimentally. Absorption also affects release and shear waves as is observed, while nonlinearities would affect mostly high amplitude P-shocks.

B. Effect of pressure level and diameter of the laser source

Regarding the laser illumination, two parameters were changed: the laser energy density and the diameter of the laser spot. For the energy, three configurations are presented, with 2 mm thick titanium plates illuminated by a laser spot of 3 mm of diameter. The laser energy density changes from 6.7 J/cm 2 in Fig. 9 (c), to 26.7 J/cm 2 in Fig. 10 (a) and to 41 J/cm 2 in Fig. 10 (b). It is remarkable that for all these three configurations, the computed amplitudes of the first arrival of the compression shock (P), the release wave (R) and the shear wave (S) are all in excellent agreement with measured ones. This clearly shows that the source computation by Esther provides an excellent input, and that absorption and nonlinear effects are not important over such a short propagation distance (2 mm). The shape of the most intense simulated shock is however slightly narrower than the recorded one. Second (3P) and third (5P) arrivals having propagated respectively over 6 mm and 10 mm, arrive at perfectly simulated times, but tend to be overestimated, especially at 26.7 J/cm 2 . This indicates that absorption is likely to play some role over such propagation distances, as confirmed by Fig. 6(b) in a 5 mm thick plate, where again the second arrival has significantly decayed after 15 mm of propagation. Nonlinear effects could also play some role, see following section 6 for discussion. This good agreement for titanium also confirms the adequacy of the purely elastic scheme for the simulations. For aluminum (not shown here), its lower elastic limit induces much larger deviations from the simulations in the cases of high intensity experiences. This again is discussed in last section 6. Another important parameter to consider is the size of the laser spot. Knowing that the release and shear waves em-anate from diffraction at the edge of the laser spot, the source size directly influences the arrival of these waves. For a titanium plate of 2 mm thickness illuminated by a 6.7 J/cm 2 laser, one can compare three spot sizes : 3 mm (Fig. 9(c)), 4 mm (Fig. 11(a)) and 5 mm (Fig. 11(b)). Increasing the focal size increases the propagation distance for the first release wave, which is all the more delayed compared to the pressure shock. Similarly the first shear wave arrives later. For the largest spot (5 mm), it arrives almost simultaneously with the next pressure shock / release wave (3P/3R) which may explain the amplitude decrease of the second shock.

VI. MODEL LIMITATIONS AND PLASTICITY

For sake of simplicity, we considered here a linear approach for modelling the propagation in an ideal elastic medium. This seems contradictory with the shock generation process, which is intrinsically a nonlinear one (as taken into account by Esther approach). Indeed, it is possible to identify three physical mechanisms that may contradict our approach : nonlinear elasticity, elastoplasticity and absorption. For weak shock waves, nonlinear elasticity assumes a reversible material behavior, but with a nonlinear relation between stress and strain. In an isotropic medium and for compression waves (by far the most intense ones in all contemplated cases), nonlinear elasticity would appear at first order as an additional quadratic term to the stress-strain relation, measured by a nonlinear dimensionless parameter β . The resulting equation is of inviscid Burgers' type and can thus lead to shock generation or modification. Nonlinear effects can be quantified through the distance L = 1/(β κM) necessary for significant non-linear phenomena, with κ = 1/(c L T pulse ) the wavenumber, and M = v/c L the Mach number. For aluminum, β ranges between 3,54 and 8,41 depending on the alloy, we take here 5,5 [START_REF] Jacob | Acoustic nonlinearity parameter measurements in solids using the collinear mixing of elastic waves[END_REF] . Taking the maximum value of the VISAR measured free surface velocity (v = 200m/s), divided by 2 because of doubling due to reflection at the free surface, yields a value L = 6.7 mm for a longitudinal celerity c L = 6100 m/s. This indicates that for millimetric distances here considered, nonlinear elasticity would not be fully negligible and is expected to be really significant for multiple cross-overs of the sample thickness. To quantify this, we use a numerical solver of Burgers' equation [START_REF] Coulouvrat | A quasi-analytical shock solution for general nonlinear progressive waves[END_REF] that allows us to propagate compression waves in 1D in a fully non-linear way. The result (not presented in this work) shows that velocity waveforms do not agree at all with the experience data: nonlinear elasticity leads to a wave pulse that arrives much too early and that widens a lot. Nonlinear elasticity is therefore not relevant in the present case. We thus have to simply solve the usual elastodynamics equations with linear Hooke's law.

The second model limitation deals with assumption of elasticity itself. Regarding the tensile yield strength Y 0 of both materials used in this work (503 MPa for aluminum and 880 MPa for titanium) and the stress produced by the shock within the structure (up to 4 GPa), it might be considered that the elastic assumption is not satisfied at all. However, yield stress is established for quasi-static tensile testing. In the case of shockwave propagation, the Hugoniot Elastic Limit (HEL) should be considered. The HEL is defined as the critical shock pressure at which a solid yields under the uniaxial strain of a plane shock wave [START_REF] Kanel | Elastic-plastic response of solids under shock-wave loading[END_REF] . Both limit values are related to one another by HEL = Y 0 (1ν)/(1 -2ν) with ν Poisson's coefficient. The HEL of the used materials are therefore much higher, 2.8 GPa for titanium [START_REF] Hopkins | Hugoniot and shear strength of titanium 6-4 under shock loading[END_REF] and 0.65 GPa for aluminum [START_REF] Peyre | Electromagnetic gauge study of laser-induced shock waves in aluminium alloys[END_REF] . The pertinence of the HEL is verified in Figs. 12(a) (titanium) and (b) (aluminum) showing the VISAR measured velocity of first arrivals on opposite face at different laser intensities. For aluminum, Fig. 12(b) shows a clear deviation from a sharp shockwave beyond around 80 m/s (for laser intensities above 12.87 J/cm 2 ), outlining the beginning of plastic effects. To estimate the corresponding stress σ , we just have again to divide the measured velocity by a factor 2 (because of doubling due to free surface reflection) and then recall that v = σ /ρc L for a plane, compression wave. In this case the observed value of 80 m/s corresponds indeed to a HEL of 0.65 GPa. In the case of titanium (Fig. 12(a)) the appearance of plastic effects is visible around 230 m/s (only for an intensity of 51.64 J/cm 2 ) which, following the same steps, provides an HEL of 3.08 GPa. These considerations allow us to keep the linear, elastic approach for low laser energy simulations in the case of the aluminum, but up to 3 GPa for the titanium.

The last physical mechanisms not considered in our simplified model is absorption, that can be either linear (viscoelasticity) or nonlinear (through shock formation in the case of nonlinear elasticity, or due to plastic deformations). To quantify it, we present on Figs. 12(c) and (d) the VISAR measured velocities on the opposite free surface of the three first arrivals, normalized by the peak of the first P-arrival. In the case of aluminum (Fig. 12(d)), the loss of amplitude increases with source intensity. This clearly indicates again an amplitudedependant, nonlinear phenomenon, most likely due to elastoplastic deformation arising quite soon. Moreover, when comparing arrivals 3P and 5P, one can see that for lowest intensities 6.65 and 12.87 J/cm 2 , the waveform keeps quite similar with a sharp shock structure clearly maintained, while for higher intensities, the waveform decays much more and progressively loses its typical shock profile. This again is an indication of nonlinear phenomena. The change in waveform and the visible delay in arrival time with increasing amplitude is completely opposite to what would be observed in the case of nonlinear elasticity, thus confirming this mechanism does not play any role here. For titanium (Fig. 12(c)), we observe some little absorption in the low amplitude case, approximately doubled when comparing respectively the ratios of 3P/P and 5P/P peak values. This is an indication of mostly linear viscoelastic absorption. For higher laser intensities, the decay between 3P and P peak arrivals is a bit more significant, but much less when comparing 5P to 3P peak values. The likely explanation is that, during the first sample cross-over (P-3P) some tiny material places may undergo a stress beyond the HEL limit and therefore plastic deformation and dissipation. During the second cross-over on the contrary (3P-5P), the stress does not reach the HEL limit and attenuation keeps small and linear. Another indication of limited elastoplastic effects is that the waveform hardly changes, and arrival times keep almost identical for all intensities.

VII. CONCLUSIONS

Laser shock generation and propagation in the ablation regime is a quite complex process that involves multiphysics and is governed by numerous parameters. In order to better understand and optimize this process, numerical simulations appears an efficient approach. However, full 3D simulations including all multiphysics with deposition of laser energy, solid to plasma phase change, material elastoplasticity up to damaging and fracturation, radiative transfer, heat conduction, wave and shock propagation, diffraction by source edges and material viscoelastic absorption is beyond current numerical capacities. Therefore we propose here a much less demanding approach, combining a fully nonlinear, multiphysical approach at 1D to describe the source term, localized in the tiny region of material plasma expansion, and a 3D (here with geometrical axisymmetry) linear elastomechanical wave propagation over the bulk of the target material. The present article demonstrates the validity of this simplified approach. Numerical convergence tests allowed us to determine the optimal numerical parameters in order to minimize the numerical dispersion during the propagation phase. Diffraction by the edges of the laser focal spot turns out essential to quantify the release and shear waves. A 2D approach strongly underestimates these ones, as in phase diffracted waves would emanate from two points instead of a full circle in 3D. This 3D approach, here simplified for axisymmetric geometries, has been extensively validated by systematic comparison with data obtained with Hephaïstos laser facility, varying target material, sample thickness, diameter of laser spot and energy density. In all cases, an excellent agreement is obtained between experiments and simulations for the overall shape of the recorded velocity at the opposite face, and times of arrivals of various waves. All of them are easily explained by linear elastic wave propagation, with compression shocks emanating from the laser plasma expansion, and release and shear waves due to wave diffraction at the sharp edges of the spot. Mode conversions from compression to shear waves (or vice versa) result from reflexions at the stress-free surfaces. Internal stress evaluation indicate possible constructive interferences, leading to increased traction within the material. Good agreement for amplitudes of the first arrivals also confirms Esther simulations provides an excellent source term for the propagation part of the code. Agreement is not so good after multiple reflexions, indicating that material absorption cannot be neglected for thick samples or long travel paths with multiple reflexions. The observed absorption is, at least in part, due to an elastoplastic behavior when stress goes beyond the Hugoniot Elastic Limit. For titanium, plasticity does not perturb much the prediction of waveforms and arrival times, at least for the considered intensities. On the contrary, for aluminum, elastic domain is limited to the lowest values of laser intensities. Nonlinear elasticity seems to play no role at all, even though the laser source is highly nonlinear and emits a sharp shock. Future works will include absorption and elastoplastic mechanisms in the wave propagation part of the algorithm, consider structural assemblies including aeronautical composite materials, and optimize traction within the bulk of the material for index n denotes discretized time t n = n∆t, while i (resp. j) is for discretized spatial variable x i = i∆x (resp. y = j∆y), ∆x, ∆y and ∆t denoting the grid steps for the space and time variables. Note the spatial discretization at point i + 1/2 is computed by a centered finite difference through a staggered grid, involving points i + 1/2 ± 1/2 and i + 1/2 ± 3/2. Thus different variables are computed on different nodes, staggered by half of a grid step, according to the arrangement of Fig. B1(a) imposed by elastodynamics equations: if horizontal velocities are allocated at the cell vertices, then vertical velocities lie at the cell centers and stress components are placed in the middle of the edges. Material parameters are also to be allocated accordingly, density at the same nodes as velocity components, stiffness components c 11 , c 13 and c 33 at the same nodes as σ r and c 55 at the same nodes as σ rz . Similarly, staggering is also used for time dependency, using the leap-frog algorithm with standard centered finite difference at second order:

∂ M ∂t

FIG. 2 .

 2 FIG. 2. (a) Left: Esther simulation of the time evolution of the pressure generated by a laser shock of 26.7 J/cm 2 exposition in Aluminum with water as confinement. Inset right in red the first instants of the pressure wave. Inset left in blue, temporal profile of the laser intensity. (b) Right: Over the black background, CCD image of the 3 mm laser spot. In blue, the measured spatial intensity profile in arbitrary units, in red the profile used for the simulations.

FIG. 3 .

 3 FIG. 3. (a) Normalized Esther simulation of the time evolution of the pressure generated by a laser shock for different laser energy density in aluminum. Zoom in the first 40 ns. (b) Ablation pressure in aluminum as function of laser energy density for Esther simulations (red crosses) compared to an analytical model 30 (blue line).

  ) for time waveform and Fig. 2(b) for spatial distribution on the illuminated side, (2) stress free condition on the opposite face and (3) PML (Perfectly Matched Layer) conditions for the two lateral faces 363738 .

  FIG. 4. Computed velocity at the center of the opposite face for a 2 mm thick titanium plate illuminated by a 26.7 J/cm 2 laser pulse with a focal spot diameter of 3 mm. (a) Comparison between four CFL conditions. (b) Comparison between four grid size with zoom on the first longitudinal arrival in (c) and on the second one and first shear arrival in (d).

FIG. 5 .

 5 FIG. 5. (a) Comparison between numerical simulation (red continuous line) with a 2D approach and VISAR measurement of the opposite free surface velocity (blue dashed line). Case of a 2 mm thick titanium plate illuminated by a 6.7 J/cm 2 laser pulse with a focal spot diameter of 3 mm. The dashed ovals indicate the release and shear waves. (b) Same as Fig. 5(a) but with a 3D axisymmetric simulation.

FIG. 6 .

 6 FIG. 6. Case of a laser shock (41 J/cm 2 ) of 3 mm focal spot diameter in a 5 mm titanium plate. P, R and S note the pressure, release and shear waves respectively. a) Time-space representation of the computed velocity along the thickness at the central position. The colorbar measures the velocity [m/s] but it is saturated in order to present in the same image all the waves involved in the phenomena. b) Opposite face velocity: comparison between simulations (red line) and VISAR measurement (blue line).

FIG. 8 .

 8 FIG. 8. (a) Same as Fig. 7(a) but for σ zz stress component. (b) stress component σ zz time waveform at three depths identified by their colour on Fig.8(a)

FIG. 10 .

 10 FIG. 10. Comparison of experimental and numerical results in a 2 mm thickness titanium plate illuminated by a laser of 3 mm diameter and of different energies per surface unit. (a) 26.7J/cm 2 (b) 41 J/cm 2

FIG. 11 .

 11 FIG. 11. Comparison of experimental and numerical results of both 2 mm thickness titanium, 6.7 J/cm 2 laser shock for different spot sizes. (a) 4 mm diameter spot (b) 5 mm diameter spot

FIG. 12 .

 12 FIG. 12. VISAR measurements of the opposite free surface velocity for four different laser intensities (plate thickness : 2 mm -focal spot: 3 mm). a) titanium -zoom on the first arrival, b) aluminumzoom on the first arrival, c) titanium -normalization by the peak of the first arrival and d) aluminum -normalization by the peak of the first arrival.

  and stress components are computed at interleaved time steps, as schematized by Fig. B1(b).

FIG. 13 .

 13 FIG. 13. (a) Spatial staggering of Virieux scheme: relative positions of nodes where velocity components, stress components and material properties are evaluated. (b) Staggering of leapfrog temporal scheme.

TABLE I

 I

	Aluminum 2810	71.7	0.33	0.65 6140	3100
	Titanium 4430	113	0.34	2.8	6060	3120

. Mechanical properties of the materials used in the laser experience. (See Appendix A for further explanation)

bonding inspection.

VIII. DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

IX. SUPPLEMENTARY MATERIAL

The reader can find in supplementary material two files. First, the 3D axisymmetric code used in this publication. The code consist of a main code, where the user can introduce the initial condition, the equations are developed and can see the result of the simulation, and four accessory functions. Further information is to be found in the zip file in order to run the code. Please cite this article in case of use of it. Second, the non-linear simulation cited in section VI and used to verify the model limitation.

Appendix A: Equations of elastodynamics

This Appendix recalls the equations of linear elastodynamics used in the present paper, in both the 2D (with parameter m = 0) and 3D axisymmetric (with parameter m = 1) cases for an isotropic material of density ρ. Variable t denotes time, while z is the spatial variable in the plate thickness (z=0 is the illuminated interface) and r measures distance from the laser optical axis. The corresponding velocities are respectively u z and u r . The stress tensor components are noted for the sake of simplicity σ z = σ zz , σ r = σ rr and, in the 3D case σ θ = σ θ θ . Equations (A1) and (A2) are the momentum equations relating the medium material velocity to the stress divergence. Eqs.(A3) to (A6) are the time derivative of the linear Hooke's constitutive equations.

The components of the stiffness matrix C are related to the Lamé parameters λ and µ by:

(A7)

with Lamé parameters linked to Young modulus E and Poisson coefficient ν by We here use Virieux numerical scheme [START_REF] Virieux | P-sv wave propagation in heterogeneous media: Velocity-stress finite-difference method[END_REF] , generalized to fourth order for spatial discretization [START_REF] Bardy | Numerical study of laser ablation on aluminum for shock-wave applications: development of a suitable model by comparison with recent experiments[END_REF] , to solve the first order hyperbolic system of equations of elastodynamics recalled in Appendix A. Spatial partial derivative in the x direction is approximated by