# Propagation of laser weak shock waves in a three-dimensional woven composite

Eduardo Cuenca, Mathieu Ducousso, Nicolas Cuvillier, Laurent Berthe, François Coulouvrat

► To cite this version:

Eduardo Cuenca, Mathieu Ducousso, Nicolas Cuvillier, Laurent Berthe, François Coulouvrat. Propagation of laser weak shock waves in a three-dimensional woven composite composite. Journal of the Acoustical Society of America, 2025, 157 (4), pp.2942-2955. 10.1121/10.0036496. hal-05056516

HAL Id: hal-05056516

https://hal.sorbonne-universite.fr/hal-05056516v1

Submitted on 5 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Propagation of laser weak shock waves in a three-dimensional woven composite

Eduardo Cuenca,<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> Mathieu Ducousso,<fnref n="1" />,<fnref n="3" /> Nicolas Cuvillier,<fnref n="1" /> Laurent Berthe,<fnref n="3" /> and François Coulouvrat<fnref n="2" />,ª)

## ABSTRACT:

The propagation of a laser-driven shock wave in an aeronautic composite material is investigated. The material, made of three-dimensional carbon fiber lattices embedded in an epoxy matrix, is heterogeneous and anisotropic, due to the intrinsic anisotropy of the carbon fibers and to the weaving process. The shock is generated by a 10 ns laser pulse, focused on the material surface. Its ablation results in an expanding plasma, which induces a shock wave in the material with a peak compression stress of a few GPa. Starting from an optical microscopy visualization of the weaving, a differentiation between resin and fibers and a segmentation of the fibers lead to an evaluation of the material's local elastic properties below the laser spot position. The shock propagation is simulated using a nonlinear source model, combined with a time domain finite difference discretization of the equations of linear elastodynamics with Lebedev's scheme adapted to the material anisotropy. The high frequency content of the signal, the material heterogeneity and anisotropy induce a complex propagation. A measurement campaign has been performed for several samples and repeated laser illuminations. Experimental data are statistically compared to the model outputs and discussed.

© 2025 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1121/10.0036496

(Received 19 July 2024; revised 1 April 2025; accepted 3 April 2025; published online 17 April 2025)

[Editor: Michael R Haberman]

Pages: 2942–2955

## I. INTRODUCTION

Reduction of carbon dioxide (CO₂) and other greenhouse gas emissions is a major challenge faced by the aeronautical transport sector. One of the most efficient ways for this alleviation is to replace metallic parts of the structure by composite materials. Safran, in its new LEAP (Leading Edge Aviation Propulsion) engine (manufactured by CFM International, a 50/50 joint venture between Safran and General Electric) uses fan blades made of organic composite materials, with only a titanium leading edge attached to each of them. The fan carter is also made of this same composite material. The total weight saving of the LEAP engine compared to the CFM56 is more than 400 kg per aircraft. This fan weight reduction allows to enlarge the engine diameter, increase the engine by-pass-ratio and reduce noise emissions. The LEAP engine is able to decrease the fuel consumption, and therefore its CO₂ emissions, by 15% compared to previous generation engines, while also reducing NOₓ emissions by 50% and satisfying Chapter 14's most stringent noise standards from the International Civil Aviation Organization. The considered material is made of carbon tows or yarns embedded in an epoxy matrix and woven in the three directions of space. Carbon tows themselves result from the packing and twisting of carbon

microfibers. Three-dimensional (3D) weaving pattern with tows waving up and down while joining together multiple plies, ensures mechanical strength through the material thickness (Chen et al., 2011; Mendoza et al., 2019). The composite material is therefore heterogeneous and doubly anisotropic, due to the intrinsic anisotropy of the carbon fibers and to the weaving process. A review of the mechanics of 3D woven material by (Wielhorski et al., 2022) outlines the richness and complexity of this issue. It distinguishes mostly predictive approaches attempting to predict the as-woven and as-manufactured states of woven fabrics starting from yarn properties and textile design, and descriptive approaches extracting the mesoscale geometry of actual preforms from x-ray micro-computed tomography (μ-CT) images. Both are said to be challenging and complementary.

To inspect bondage quality between the composite fan blade and its titanium leading edge, shock waves from a laser-generated plasma (Berthe et al., 1997) are an appealing approach because they induce sufficiently large mechanical traction forces at the bonding interface (Berthe et al., 2011; Ducousso et al., 2018). However, in case of complex materials such as woven composites, weak shock wave propagation in such complex media remains poorly understood. Homogenization (Carvelli and Poggi, 2001) is not adapted because the shock wavelength, on the order of several tens

<footnote n="1">¹Safran Tech, rue des Jeunes Bois, 78114 Magny les Hameaux, France</footnote>

<footnote n="2">²Sorbonne Université, CNRS, Institut Jean Le Rond d'Alembert, 4 place Jussieu, 75005 Paris, France</footnote>

<footnote n="3">³Laboratoire Procédés et Ingénierie en Mécanique et Matériaux, CNRS, Arts et Métiers Paris Tech, 151 boulevard de l'Hôpital, 75013 Paris, France</footnote>

<footnote>ª)Email: francois.coulouvrat@sorbonne-universite.fr</footnote>
of µm, is much smaller than the size of a carbon tow, on the order of several mm. Exact models for wave propagation in layered anisotropic media based on transfer matrix (Nayfeh, 1991; Rokhlin and Wang, 2002) are not applicable because of the fully 3D weaving process. A numerical approach appears more relevant, combining a mesoscale model describing the mechanical behavior of tows, and a direct simulation of wave propagation through the geometrical entanglement of tows embedded in the matrix material. Such an approach was achieved by Thierry et al. (2020) using wave finite elements, but for frequencies not higher than 2 × 10⁵ Hz, about hundred times smaller than those here considered. Shock wave propagation with characteristics similar to those here considered, was recently simulated by Tserpes and Kormpos (2024). It uses a detailed finite element (FE) model, taking into account in particular the variable orientation of tows. Comparison with a homogenized model points out considerable differences in the back-face velocity, with more internal wave reflections occurring in the FE model, and a significant effect of spot locations (two of them were considered) affecting the shock propagation path and the resulting back-face velocity.

Inspired by this, we also propose here a numerical approach for shock wave propagation at the mesoscale. Two types of 3D woven composites are considered, with a total of five samples of various thicknesses. For each sample, eight experimental laser shots and eight simulations are performed. This allows us to quantify the variability of the resulting back-face velocity signal when varying the location of the laser illumination spot, in terms of waveform, time of arrival and shock amplitude, with a quantitative comparison between simulations and experiments. For each of the 40 cases, local geometry is extracted from image analysis of cross-sections of the considered sample below the location of the shock loading, providing the local orientation of tows. Thus, following Wielhorski et al. (2022), we combine a predictive approach based on a numerical mechanical model, and a descriptive approach of the local geometry. However, to be tractable, we have to simplify the model. Compared to the few finite element approaches of the literature for shock wave propagation in 3D woven composites (Pankow et al., 2011; Ghosh and De, 2014; Tserpes and Kormpos, 2024), the present finite-difference time domain (FDTD) model is simpler to implement as it avoids the complex and lengthy process of meshing the intricate geometry of interfaces between the various material components (Su et al., 2023). Manageable computation times and the analysis of cross section images restrict our simulations to 2D cases (in-plane strain assumption). The laser shock source is described by a nonlinear approach (Berthe et al., 1997; Colombier et al., 2005) valid within a thin metallic layer above the sample, while the propagation is assumed linear. This hybrid approach has been proved successful for laser-generated shock wave propagation in homogeneous metals such as aluminum or titanium (Cuenca et al., 2020) and is here extended to composite materials.

In Sec. II, we will describe the composite, with a focus on the material heterogeneity and on the local orientation of tows.

In Sec. III, the experimental setup based on the generation and detection of shock-waves by laser-induced plasmas will be presented. The numerical approach will be introduced in Sec. IV, and comparisons between numerical simulations and observations will finally be discussed (Sec. V).

## II. DESCRIPTION OF THE MATERIAL

### A. Components and manufacturing process

The Tissue3D material is the 3D woven composite material designed, developed and produced by Safran Group, in partnership with the American company Albany International, to constitute the fan blades of its LEAP engine. A three-dimensional weaving process allows one to produce complex woven reinforcements with reduced weight compared to metallic materials, enhanced mechanical properties compared to classical composites, especially resistance to delamination and impact in the three directions of space, and adaptation to the desired shape of the product. It belongs to the class of carbon-fiber reinforced epoxy composites. Figure 1 [top figure, obtained from Mendoza Quispe (2019)], presents a µ-CT scan of the composite sample showing the fibers waving up and down while joining together multiples plies and therefore ensuring mechanical strength through the overall material thickness.

This 3D view is complemented (Fig. 1, bottom figure) by a photography of the cross section of the Tissue3D obtained with a microscope after cutting the sample with a microtome SECOTOM-50 cutting machine. The surface is then polished and photographied by a microscope Axio Observer (Zeiss) using RGB mode with a 2.5× zoom. Scales are not indicated for reasons of industrial

<figure><img src="image_2.png" /><figcaption>FIG. 1. Top: µ-CT Scan of a Tissue3D sample (Mendoza Quispe 2019). Bottom: photography of one cross section. The white structure corresponds to the epoxy resin. Carbon fiber reinforcement appears either in gray (in-plane fibers) or in black (out of plane fibers).</figcaption></figure>
confidentiality. Two different components can be distinguished. The darkest one is the carbon fibers, more precisely PolyAcrylonitrile (PAN) based fibers manufactured by Hexcel. The PAN fibers go through a complex manufacturing process involving high temperatures (above 1000 °C) and different atmospheres to be transformed into carbon microfibers (Khayyam et al., 2020). First is the precursor spinning: a PAN polymer dope, after being filtered, passes through the holes of a spinneret immersed in an aqueous sodium thiocyanate solution. This makes the fibers precipitate and then solidify with the form of a thread of diameter around 5 µm. After thermal stabilization, microfibers pass through heated rollers for drying. They are then packed and twisted in groups of a few thousands of them. The resulting thicker threads, called tows, therefore, have millimetric diameters. These tows are then intertwined between one another to conform the 3D reinforcement structure of the composite. We want to clarify for the reader that from now on these tows are going to be referred as “fibers.” The term microfibers will be used only for the elementary fibers of micrometric size that are part of the tows. The white material of the image is the epoxy resin, manufactured by Solvay. This class of reactive polymers, also called poly-epoxides, contains at least one epoxide group. It reacts with itself or with a co-reactant to create a cross-linking reaction known as curing, during which the mono-epoxy molecules interact with the amine reactants to create poly-epoxides (Tesoro, 1988). During the manufacturing process, the reinforcement made of 3D woven carbon fibers is woven as a raw form of the fan blade. Then, with a water jet cutter, this form is shaped and placed into a mold. The epoxy-amine resin is injected with high pressure into the mold. The curing process then takes place during a couple of hours at around 180 °C. The resin fills the carbon fiber reinforcement and gives integrity to the structure when cured. It fills not only the voids of the reinforcement, but also the space between the microfibers in the tows. Hence, the tows themselves can be considered as a heterogeneous and anisotropic material. The fiber volume fraction (FVF) which is defined as the percentage ratio between carbon fibers and epoxy resin in a considered volume, ranges between 80 and 20 and 90 and 10 within the tows. It is called the intratow FVF. The value of the overall FVF in the entire composite is different, but ranges between the same values. In Fig. 1(b), two different types of tows are visible. The first ones, visible as darker ellipses, indeed cross-the image plane mostly perpendicularly. The second ones, visible as gray thick threads, extend mostly within the image plane. Note that both “types” of carbon fibers are made of exactly the same material. Only their orientation is different. Fibers are woven following a layer-to-layer angle interlock. This means that even if the fibers go mainly along the horizontal axis x (for gray ones) or the transverse one z (for black ones), they also cross-two or more layers of fibers in the y axis. Regarding the weaving specifications, even though the sample geometry is 3D, an inspection of their cross section shows samples have a quasi-2D periodic horizontal pattern that is repeated all

along. In the following, two groups of materials will be examined, the second type having a quasi-period larger than the first one. Since the shock wave investigated in this manuscript will propagate mostly through the sample thickness in the y-direction, the disposition and orientation of the fibers that pack the sample will have a strong influence on the sample mechanical properties, and therefore on the wave propagation. Thus, due to its manufacturing process, each composite sample considered at the shock wave scale is going to be unique, strongly heterogeneous (because of the different materials used in the composite) and anisotropic (because of the manufacturing of the microfibers themselves, and because of the woven arrangement).

## B. Mechanical properties

The overall mechanical properties of Tisse3D are of course dependent on those of its individual constituents, epoxy resin, and carbon fibers. Mechanically speaking, the resin has an isotropic behavior, while a carbon fiber is recognized as a transversely isotropic material from a macroscopic perspective (Shirasu et al., 2020). Note that the typical wavelength is on the order of 50 to 100 µm, therefore much larger. Due to the manufacturing process of tows from individual microfibers, they are much stiffer in their longitudinal direction X than in the transverse one T. The objective of the weaving process is to make the Tisse3D material benefit from this high longitudinal stiffness of individual tows in all directions by having fibers orientated in all three directions. Mechanical engineering constants of both materials are reported in Table I. Values are, however, rounded for reasons of confidentiality, since they arise from data sheets and internal documents from Safran.

Assuming a linear elastic behavior locally for each phase, either transverse isotropic (for carbon fibers) or fully isotropic (for epoxy), the relations of linear elasticity between strain $\epsilon$ and stress $\sigma$ tensors are

$$ \sigma = C(x) \cdot \epsilon \Rightarrow \frac{\partial \sigma}{\partial t} = C \cdot \dot{\epsilon} \quad (1) $$

with position x indicating the local phase and, for fibers, its local orientation as will be detailed below. The strain tensor is defined as $\epsilon_{ij} = 1/2(\partial\delta_i/\partial x_j + \partial\delta_j/\partial x_i)$ with $\delta_k$ the material displacement in the kth direction. Dots indicate time derivative. Using Voigt notation ($T_1 = T_{XX}$, $T_2 = T_{YY}$, $T_3 = T_{ZZ}$, $T_4 = dT_{YZ}$, $T_5 = dT_{XZ}$, and $T_6 = dT_{XY}$ for both $T = \sigma$ stress and $T = \epsilon$ strain tensors with $d = 1$ for stress and $d = 2$ for strain), C can be written as a six-by-six

<table><caption>TABLE I. Mechanical properties of the carbon fiber and the epoxy resin. The letter E is for Young's modulus, G for shear modulus (unit: GPa), $\nu$ for Poisson's ratio and $\rho$ for density (unit: kg/m$^3$). Index T is for the transverse plane, and index x is for the fibers axial direction.</caption><thead><tr><th></th><th>E<sub>x</sub></th><th>E<sub>T</sub></th><th>G<sub>x</sub></th><th>ν<sub>x</sub></th><th>ν<sub>T</sub></th><th>ρ</th></tr></thead><tbody><tr><td>Carbon fiber</td><td>250</td><td>15</td><td>25</td><td>0.3</td><td>0.3</td><td>1780</td></tr><tr><td>Epoxy resin</td><td>5</td><td>5</td><td>1.9</td><td>0.316</td><td>0.316</td><td>1250</td></tr></tbody></table>
symmetric stiffness matrix. In the assumption of transverse isotropy (a particular case of orthotropic anisotropy) in the symmetry plane (Y, Z) with X the axial direction of fibers normal to this plane, the mechanical properties of the carbon fibers are described by five independent mechanical parameters reported in Table I with $E_X$ and $\nu_X$, respectively, the longitudinal Young's modulus and Poisson's ratio in the axial direction of fibers, $E_T$ and $\nu_T$, respectively, the Young's modulus and Poisson's ratio in the transverse symmetry plane and $G_X$ the shear modulus in the axial direction.

With these notations, the only non-zero stiffness coefficients are (Royer and Valier-Brasier, 2022)

$$C_{fib} = \begin{bmatrix} C_{11} & C_{12} & C_{12} & 0 & 0 & 0 \\ C_{12} & C_{22} & C_{23} & 0 & 0 & 0 \\ C_{12} & C_{23} & C_{22} & 0 & 0 & 0 \\ 0 & 0 & 0 & C_{44} = \frac{C_{22} - C_{23}}{2} & 0 & 0 \\ 0 & 0 & 0 & 0 & C_{55} & 0 \\ 0 & 0 & 0 & 0 & 0 & C_{55} \end{bmatrix} \quad (2)$$

with the following relations (Dowling, 1999; Ding et al., 2006) between the mechanical engineering parameters given in Table I and the stiffness matrix:

$$ \begin{aligned} C_{11} &= E_X + 4\nu_X^2 K_T, \\ C_{12} &= C_{13} = C_{21} = C_{31} = 2\nu_X K_T, \\ C_{22} &= C_{33} = K_T + G_T, \\ C_{23} &= C_{32} = K_T - G_T, \\ C_{44} &= G_T, \\ C_{55} &= C_{66} = G_X, \end{aligned} \quad (3) $$

where the transverse shear modulus $G_T$ and the transverse, plane-strain bulk modulus $K_T$ are given by

$$G_T = \frac{E_T}{2(1 + \nu_T)} \quad \text{and} \quad \frac{1}{K_T} = \frac{2(1 - \nu_T)}{E_T} - \frac{4\nu_X^2}{E_X}. \quad (4)$$

Note that these mechanical parameters are of the same order of magnitude as those used by Tserpes and Kormpos (2024).

With values of Table I, one has $K_T = 10.9$ GPa while the coefficient of axial rigidity $C_{11} \approx E_X$ is much larger than other ones. For epoxy (index $E$) assumed to be fully isotropic, the above equations are still valid with $\nu_T = \nu_X = \nu_E$, $G = G_T = G_X = G_E$, and $E_E = E_T = E_X = 2(1 + \nu)G_E$. In this case we recover the stiffness matrix of an isotropic material with $C_{44} = C_{55} = C_{66} = G_E = \mu_E$, $C_{12} = C_{13} = C_{23} = 2\nu_E G_E / (1 - 2\nu_E) = \lambda_E$, and $C_{11} = C_{22} = C_{33} = 2G_E (1 - \nu_E) / (1 - 2\nu_E) = \lambda_E + 2\mu_E$, where $\lambda_E$ and $\mu_E$ are Lamé's coefficients of epoxy.

Equation (1) assumes a linear elastic behavior which can be questionable for relatively strong shocks. However, laser intensity has been limited to levels similar to those investigated by Cuenca et al. (2020) in metals [aluminum (Al) and titanium]. The numerical approach in this reference

is similar to the present one and is hybrid, in the sense that the laser shock emission resulting from the illuminated surface ionization is strongly nonlinear, but the resulting bulk shock propagation is assumed to be linear. The extensive comparison in Cuenca et al. (2020) between numerical hybrid simulations and experiments with varying materials, sample thicknesses and laser intensities lead to a remarkable agreement, thus validating the hybrid approach and the assumption of bulk linear elasticity. Only the source emission is nonlinear, and numerically appears as a boundary condition for bulk, elastic propagation. In the present case of composite materials, this assumption will be all the more valid as the material heterogeneity will induced strong wave scattering and further reduce the shock amplitude.

### C. Image analysis and fibers orientation

To simulate the shock wave propagation, we need to separate numerically the fibers from the epoxy matrix and determine the local orientation of in-plane fibers. Indeed, due to the 3D weaving process, principal axes of symmetry of the fibers (X, Y, Z) generally do not correspond to propagation axes (x, y, z) shown in Fig. 1(b). Axis y is the main propagation direction oriented along the laser beam through the sample thickness, x is in the experimental cross section and parallel to the sample surfaces, and z in the third direction. As described in Cuenca et al. (2020) (see their Fig. 2), the laser source impulsion leads to a pressure pulse at the top of the sample of typical duration 100 ns with a half-width duration about 15 ns, and thus to a shock spectrum peaking in the 10 MHz frequency range but extending beyond 60 MHz. With stiffnesses on the order of tens of GPa in fibers, wavelengths corresponding to the peak of the spectrum are typically on the order of a tenth of a millimeter, quite smaller than the transverse size of the fibers (a few millimeters). The wavelengths associated with the half-width duration are even smaller. Therefore, a homogenization of the shock propagation is not possible because the problem is neither periodic nor at low frequency. To simulate it, it is therefore necessary to know the accurate geometry of fibers entanglement below each laser spot. Unfortunately, this is difficult to do repeatedly in 3D, and we had to limit in the present study to 2D images of material plane cross-sections [in the (xy) plane]. This limited 2D knowledge of the material geometry implies that simulation could be performed only by assuming an in-plane strain [in the same (xy) plane].

Cross-section images such as Fig. 1(b) are numerically inspected by means of FIJI (Fiji, 2025), an open source image processing package with many bundled plugins developed initially for the analysis of biomedical images (Abramoff et al., 2003; Schindelin et al., 2012). For the image analysis, we first create a MATLAB code that provides the various gray levels of each pixel within one image. Then we select the maximal gray range for the epoxy pixels, which appear as the lightest ones, and the minimal one for the perpendicular fibers which are the darkest ones. The intermediate values
correspond to the in-plane fibers. The different types of materials are, however, not perfectly separated: clouds of isolated light or dark points are spread over most of the image. To improve the material separation, we then use the IMAGEJ (Image Processing and Analysis in Java) plugin from FIJI to “clean” the material data by hand. To obtain the local orientation, we use within the OrientationJ plugin of FIJI, a tool called Vector Field that measures the local orientation of the fibers from an image. The procedure starts by transforming the image into a 16 bits image. We then select the local window (number of pixels for the orientation calculation) and the gradient of position, computed in our case with the recommended method of Cubic Spline interpolation. It provides the orientation in each pack of pixels determined by the “grid size” parameter, along with the position on the grid and the gradient in the x and y directions. Since we need a certain number of pixels to calculate the gradient in a certain zone, the dimensions of the calculated orientation matrix are reduced along with the grid size. For instance, if the original size of the image is 1000 × 100 pixels and the grid size is 10, the orientation matrix is a 100 × 10 matrix. To fit the orientation matrix to the original pixel size of the composite, we assume that every point within one grid size has the same orientation value, thus recovering the required 1000 × 100 orientation matrix. Orientation is then conserved only for in-plane fiber pixels. Out of plane fibers are assumed to be oriented exactly normal to the image plane (in this case, axis X is z). Figure 2 provides, in the case of the cross section shown in Fig. 1(b), the resulting colormap of the local orientation angle θ of the in-plane fibers with respect to the horizontal axis x in degrees, while epoxy appears as white and out of plane fibers as black. Green points identify areas with local orientation close to 0°, this means fibers with their axis X almost parallel to the sample surface. This null orientation is the average one, but locally fibers can strongly deviate from it, with values of θ reaching up to ± 80° (dark blue and bright yellow points): such fibers are locally orientated almost through the sample thickness as a result of the 3D weaving process.

## D. Stiffness coefficients and in-plane deformation

With the simplifying assumption of in-plane (x, y) strain, the Voigt notation for the in-plane strain rate vector

<figure><img src="image_3.png" /><figcaption>FIG. 2. Example of color map of the local orientation angle with respect to the horizontal direction, in degrees. Epoxy in white and fibers transverse to the plane in black.</figcaption></figure>

$\dot{\epsilon}^{2D}$ is $\dot{\epsilon}_{1}^{2D} = \partial u/\partial x$, $\dot{\epsilon}_{2}^{2D} = \partial v/\partial y$ and $\dot{\epsilon}_{3}^{2D} = \dot{\epsilon}_{6} = (\partial u/\partial y + \partial v/\partial x)$ with u (respectively, v) the material velocity in the x (respectively, y) direction. Similarly, one defines the stress vector $\sigma_{1}^{2D} = \sigma_{xx}$, $\sigma_{2}^{2D} = \sigma_{yy}$, and $\sigma_{3}^{2D} = \sigma_{xy}$.

For in-plane fibers, we can apply Eq. (1) in the local (X, Y) plane with axis X indicating fibers axis and corresponding material velocity noted U. In the perpendicular in-plane direction Y, material velocity is V. These directions are rotated by an angle +θ compared to the computation directions (x, y). One gets

$$ \frac{\partial \sigma_{XX}}{\partial t} = C_{11} \frac{\partial U}{\partial X} + C_{12} \frac{\partial V}{\partial Y} \\ = (E_X + 4\nu_X^2 K_T) \frac{\partial U}{\partial X} + 2\nu_X K_T \frac{\partial V}{\partial Y}, \qquad (5) $$

$$ \frac{\partial \sigma_{YY}}{\partial t} = C_{12} \frac{\partial U}{\partial X} + C_{22} \frac{\partial V}{\partial Y} \\ = 2\nu_X K_T \frac{\partial U}{\partial X} + (K_T + G_T) \frac{\partial V}{\partial Y}, \qquad (6) $$

$$ \frac{\partial \sigma_{XY}}{\partial t} = C_{55} \left( \frac{\partial U}{\partial Y} + \frac{\partial V}{\partial X} \right) = G_X \left( \frac{\partial U}{\partial Y} + \frac{\partial V}{\partial X} \right). \qquad (7) $$

It can be written

$$ \frac{\partial \sigma^{2D}}{\partial t} = C_{\text{fib}}^{2D} \cdot \dot{\epsilon}^{2D} = \begin{bmatrix} C_{11} & C_{12} & 0 \\ C_{12} & C_{22} & 0 \\ 0 & 0 & C_{55} \end{bmatrix} \cdot \dot{\epsilon}^{2D}. \qquad (8) $$

The equations in the referential (x, y) are deduced by introducing a rotation of angle −θ, thus transforming the rigidity matrix $C_{\text{fib}}^{2D}$ in the locally tilted one $C_{\text{loc}}$ (see details in Appendix A),

$$ C_{\text{loc}} = T(\theta) \cdot C_{\text{fib}}^{2D} \cdot D \cdot T(-\theta) \cdot D^{-1} \qquad (9) $$

with matrices T(θ) and D defined in Appendix A. In the computational frame of reference (Oxy), the constitutive equations take the following form:

$$ \frac{\partial \sigma_{xx}}{\partial t} = C_{11}^{\text{loc}} \frac{\partial u}{\partial x} + C_{12}^{\text{loc}} \frac{\partial v}{\partial y} + C_{13}^{\text{loc}} \left( \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \right), \qquad (10) $$

$$ \frac{\partial \sigma_{yy}}{\partial t} = C_{12}^{\text{loc}} \frac{\partial u}{\partial x} + C_{22}^{\text{loc}} \frac{\partial v}{\partial y} + C_{23}^{\text{loc}} \left( \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \right), \qquad (11) $$

$$ \frac{\partial \sigma_{xy}}{\partial t} = C_{13}^{\text{loc}} \frac{\partial u}{\partial x} + C_{23}^{\text{loc}} \frac{\partial v}{\partial y} + C_{33}^{\text{loc}} \left( \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \right). \qquad (12) $$

For example, an orientation θ = 10°, which is close to the mean value found in the Tissue3D composite, changes the stiffness matrix from

$$ C_{\text{fib}}^{2D} = C_{\text{loc}}(\theta = 0^{\circ}) = \begin{bmatrix} 254 & 6.53 & 0 \\ 6.53 & 16.7 & 0 \\ 0 & 0 & 25.0 \end{bmatrix} \text{MPa} \qquad (13) $$

to
$$C_{\text{loc}}(\theta = 10^{\circ}) = \begin{bmatrix} 242 & 11.1 & 32.9 \\ 11.1 & 19.2 & 7.63 \\ 32.9 & 7.63 & 29.6 \end{bmatrix} \text{MPa.} \quad (14)$$

At 10°, we observe that the elements C₁₃ and C₂₃ that were zero at 0° turn into positive values. Coefficient C₁₃ becomes the second largest in magnitude, one seventh of C₁₁ which is little modified. Coefficient C₁₂ has increased by about 70%. Indeed the large stiffness C₁₁ of the fibers in the axial direction is little affected by a small rotation, but it strongly influences other stiffness coefficients with index 1. At 45°, all coefficients are of the same order of magnitude, in the range 45–100 GPa.

For isotropic epoxy, we can apply Eq. (1) in plane (x = X, y = Y). We therefore have

$$C^{\text{loc}} = \begin{bmatrix} \lambda_E + 2\mu_E & \lambda_E & 0 \\ \lambda_E & \lambda_E + 2\mu_E & 0 \\ 0 & 0 & \mu_E \end{bmatrix}. \quad (15)$$

For out of plane fibers, we can apply Eq. (1) in the plane (x = Z, y = Y). We therefore have

$$C^{\text{loc}} = \begin{bmatrix} K_T + G_T & K_T - G_T & 0 \\ K_T - G_T & K_T + G_T & 0 \\ 0 & 0 & G_T \end{bmatrix}. \quad (16)$$

It has, as expected, the same isotropic form with transverse Lamé's coefficients $\lambda_T = K_T - G_T$ and $\mu_T = G_T$, the relation $K_T = \lambda_T + \mu_T$ being the usual one between bulk modulus and Lamé's coefficients in 2D space.

The constitutive equations Eqs. (10)–(12) with stiffness matrix C<sub>loc</sub> given by Eq. (9) for in-plane fibers, by Eq. (15) for epoxy and by Eq. (16) for out of plane fibers varying locally with material and local orientation (for in-plane fibers), are complemented by the equations of momentum balance

$$\rho \frac{\partial u}{\partial t} = \frac{\partial \sigma_{xx}}{\partial x} + \frac{\partial \sigma_{xy}}{\partial y}, \quad (17)$$

$$\rho \frac{\partial v}{\partial t} = \frac{\partial \sigma_{xy}}{\partial x} + \frac{\partial \sigma_{yy}}{\partial y}, \quad (18)$$

with values of density ρ also varying locally with material only.

To solve the equations of motion, it is necessary to introduce boundary conditions. Because they are strongly linked to the experimental setup, they are discussed in Sec. III describing this setup.

### III. EXPERIMENTAL SETUP

As sketched in Fig. 3(A), an intense Nd:YAG laser (Thales, France) emits a pulse at 532 nm with a repetition frequency of 2 Hz and an energy up to 7 J with a Gaussian temporal shape of duration 7 ns. The laser beam is focused

<figure><img src="image_4.png" /><figcaption>FIG. 3. (A) Scheme of the experimental setup. (B) Photography of the experimental setup with sample and DOE—note that here, the sample mounting system does not allow a VISAR measurement. (C) Photography of the rear, silver-coated surface of two samples. (D) Photography of illuminated Al-adhesive layer after 8 shots—dotted lines indicate the target shots alignment.</figcaption></figure>

by an optical lens of focal diameter 4.7 mm. A diffractive optical element (DOE), placed after the focusing lens [see Fig. 3(B)], leads to an almost ideal “top-hat” circular and spatially uniform intensity profile (standard deviation of only 0.023) (Cuenca et al., 2020). After surface cleaning with ethanol, the top surface of the Tisse 3D target is covered by a thin (70 µm) Aluminum adhesive tape. Any possible fold is removed thanks to a plastic card. The intense laser beam therefore ablates superficially this metallic layer over a thickness typically of 1.7–3 µm (Scius-Bertrand et al., 2021), hence sublimizing it into a dense plasma. Water, a transparent dielectric material, confines the expansion and therefore increases the surface pressure level imposed by the plasma expansion compared to the direct regime (Berthe et al., 1997; Bras et al., 2019). As a result of the recoil moment of the ablated momentum, a pulsed shockwave is generated and propagates in the bulk of the Al layer. It is then transmitted to the Tisse 3D material. Sublimation of the aluminum layer and emission of a shock wave are a well-controlled and quantifiable process (Bardy et al., 2016; Cuenca et al., 2020). The Al layer also protects the composite from the ablative laser impact, Fig. 3(D).

This experimental setup drives the boundary condition used in the numerical simulations. On the sample surfaces, the normal stress is assumed to be zero, except on the laser illuminated surface where the only non-zero stress component $\sigma_{yy} = p_{\text{load}}(t)H(x)$ describes the laser-induced pressure loading. Thanks to the use of a diffractive optical element (DOE), this one is assumed to be uniform over the laser spot surface $H(x) = 1$ and zero elsewhere $H(x) = 0$. To describe the process of laser light conversion into a mechanical stress, we rely on a 1D Lagrangian code, called Esther (Colombier et al., 2005; Bardy et al., 2016; Scius-Bertrand et al., 2021), describing both laser-matter interaction and
shock generation. This software simulates the deposition of laser energy and the evolution of target material state from solid to plasma while including radiative transfer and heat conduction. In the present work, Esther allows us to simulate the so-called ablation pressure $p_{load}(t)$ at the Al surface resulting from a laser illumination with a defined duration and intensity. One example of such temporal profile is shown in Fig. 2(a) of Cuenca et al. (2020). Using Esther simulations, it is possible to evaluate the rise time of the surface pressure, typically of several ns. Beyond the halfwidth peak around 15 ns, the pressure decreases slowly, leading to a pulse duration of typically 0.1 µs. Therefore, the frequency spectrum peaks around 10 MHz [see Fig. 3(a) in (Cuenca et al., 2020)]. The ablation pressure waveform, normalized by its peak value, little varies with the laser intensity. The peak overpressure [see Fig. 3(b) in (Cuenca et al., 2020)] follows a square root dependency with the intensity, in agreement with the analytical model proposed by (Fabbro et al., 1998)

$$ \max(p_{load}(t)) = 0.1T \left( \frac{\alpha Z I_0}{2\alpha + 3} \right)^{1/2}, \quad (19) $$

where the ablation pressure is given in kbar, $I_0$ (unit: GW/cm²) is the laser surface intensity, and $\alpha$ is the ratio of the thermal energy to the internal energy of the plasma, here chosen equal to 0.6 after calibration 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² s of, respectively, water and the illuminated material (Al). Indeed, a water droplet is used to confine the plasma expansion inside the material and enhance the efficiency of laser shock generation. Equation (21) is well known and has been confirmed extensively by Correa et al. (2015) or Cuenca et al. (2020). As the process of laser-induced material ionization is much more efficient and much better known in metals such as aluminum than in non-metallic materials, the illuminated surface of the composite 3D-woven material is covered by a thin aluminum adhesive tape. Its thickness (70 µm) is much larger than the ablation one (a few µm). We therefore assume the surface shock generation process takes place entirely in the aluminum tape (hence $Z_2$ is the impedance of aluminum), and that the resulting pressure profile $p_{load}(t)$ is fully transmitted to the composite. This transmission is described by a single transmission coefficient T that is left as a free parameter whose value is discussed at the end of this work.

For all other surfaces, a stress-free boundary condition is applied. For lateral boundaries, a non-reflective condition could have been expected, but they are sufficiently distant, and absorption of coherent signal due to scattering by the complex structure of the material is sufficiently high, so that this is not necessary.

The surface velocity of the opposite surface is measured by VISAR (velocity interferometer system for any reflector) laser interferometry (Barker and Hollenbach, 1972), a Michelson-type interferometer measuring with a ns time resolution the Doppler shift in the wavelength of a probe laser

(532 nm wavelength and 5 W power, from Coherent Company, Santa Clara, CA, USA), reflected by the sample opposite free surface and moving due to the shock wave successive arrivals. To obtain the required optical reflectivity, an ultra-thin silver layer (~100 nm) is deposited using an Electron Beam Evaporator PVD4-EB, Fig. 3(C). Since the thickness of the silver layer is much smaller than the shock wavelength (typically 100 µm), this one does not exert any traction in the coating: at the shock wavelength, the deposit layer is indistinguishable from the free surface where the mechanical stress vanishes. In Fig. 3(C) due to the very small thickness of the silver layer, which is not able to cover all the micro-holes present in the composite surface, the surface of the metallized samples is not completely smooth. This, however, does not create any problem as the VISAR system is calibrated before every shot.

The lateral dimensions of the examined samples are approximately 45 × 90 mm². Two types of 3D woven composites were illuminated, made of the same materials with the same manufacturing process, but type 2 has larger tows and a larger weaving period than type 1. Three samples of thicknesses 2, 3, and 4.7 mm of type 1 are illuminated, as well as two samples of thicknesses 3.7 and 7.5 mm of type 2. For each sample, two series of four laser shots each, parallel to one another and aligned along one of the main weaving axis, with shots and series separated by the distance of the quasi-periodic weaving pattern, are performed [see the eight surface laser impacts on Fig. 3(D)]. However, some shots turned out unsuccessful, with no shock clearly recorded by VISAR at the opposite face. These unsuccessful shots are most likely attributed to the optical detection. Indeed, VISAR relies on Doppler interferometry and thus requires an intense enough laser light reflected from the sample surface (Fleming and Broyles, 2003). In the case of intense scattering during propagation by the highly heterogeneous composite, the shock wavefront right at the VISAR laser position may be sufficiently deviated from the y axis, thus leading to a strongly tilted surface displacement with respect to the normal surface. This impairs efficient VISAR detection and likely explains most unsuccessful shots. The number of successful shots per sample is variable, ranging between 3 and 8 (Table II). The laser intensity is adjusted for each sample in order to obtain similar velocity peaks measured at opposite surface for all samples, in the range 100–200 m/s, so as to make the comparison easier and to compensate for increased amplitude losses for the thicker samples. We select the laser fluency listed in Table II by shooting on the sample a first time, measuring the

<table><caption>TABLE II. Parameters for the eight samples.</caption><thead><tr><th>Sample number</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th></tr></thead><tbody><tr><td>Thickness [mm]</td><td>2</td><td>3</td><td>3.7</td><td>4.7</td><td>7.5</td></tr><tr><td>Weaving pattern</td><td>1</td><td>1</td><td>2</td><td>1</td><td>2</td></tr><tr><td>Fluency [J/cm²]</td><td>7.3</td><td>13.5</td><td>16.21</td><td>16.43</td><td>20</td></tr><tr><td>Successful shots</td><td>3</td><td>8</td><td>8</td><td>5</td><td>5</td></tr><tr><td>Measured rise time [ns]</td><td>45.8</td><td>44.1</td><td>106.5</td><td>99.7</td><td>127.8</td></tr><tr><td>Simulated rise time [ns]</td><td>51.6</td><td>98.2</td><td>96.4</td><td>157</td><td>143.2</td></tr></tbody></table>
amplitude, and then adjusting the fluency to maintain the desired velocity amplitude. Note that a previous study (Cuenca et al., 2020) shows that such type of bulk wave propagation can be considered as linear in metals (aluminum and titanium), even though the process of shock emission is highly nonlinear. As in the present composite material sound is much more attenuated, we can expect this linear behavior to be even better satisfied here.

To examine how the structure of the composite affects the propagation of the shockwave, each sample, after having undergone a series of eight laser shots, is cut into three pieces using a SECOTOM-50 cutting machine. Both cuts are done along the two parallel lines that cross-four shots by their centers. Then the cross-sections are polished and photographed with a microscope Axio Observer Z1 of Zeiss, using RGB mode with a 2.5× zoom. Thanks to these images, we are able to correlate the VISAR measured velocity at the opposite free surface, and a 2D view of the composite structure exactly between the illumination and measurement points. This structure is used as the numerical geometrical input for the eight numerical simulations per sample, using the procedure of extraction of material and fibers orientation described previously, and the numerical method described in Sec. IV.

## IV. FDTD NUMERICAL METHOD

Among finite differences time domain (FDTD) methods, Virieux's scheme (Virieux, 1986) has proven to be highly efficient to simulate laser shock propagation in homogeneous and isotropic materials such as metals (Cuenca et al., 2020). The scheme also allows one to model propagation in the case of materials with more complex mechanical properties, such as tilted transversely isotropic elastic properties (Operto et al., 2009). However, for general anisotropic materials, Virieux's scheme cannot be applied straightforwardly because the grid staggering is not adapted anymore. For such cases, other finite difference schemes like the rotated staggered grid (Saenger et al., 2000) and Lebedev's schemes (Lebedev, 1964; Lisitsa and Vishnevskiy, 2010) have been developed. Like Virieux's scheme, Lebedev's scheme used here, relies on grid staggering both in time and in space. However, Lebedev's scheme locates all discrete velocities (noted generically W) on the same points of a first grid, and all stress components (noted generically S) on the same points of a second grid, staggered relatively to the velocity grid by half the spatial mesh $\Delta x/2$ in all directions. We can therefore note discrete spatial points ($x_i = i\Delta x/2, z_j = j\Delta x/2$) with ($i, j$) taking even values for velocity and odd ones for stress. Time grids are also staggered, with velocities computed at times $t_n = n\Delta t$ and stresses at times $t_{n+1/2} = (n + 1/2)\Delta t$. Spatial derivatives associated with one type of variable are computed by a fourth order, centered, finite difference

$$ \left(\frac{\partial W}{\partial x}\right)_{i,j}^n = D_x(W_{i,j}^n) + O(\Delta x^4), \qquad (20) $$

$$ \left(\frac{\partial S}{\partial x}\right)_{i,j}^{n+1/2} = D_x(S_{i,j}^{n+1/2}) + O(\Delta x^4), \qquad (21) $$

with

$$ D_x(A_{i,j}^n) = \frac{27(A_{i+1,j}^n - A_{i-1,j}^n) - (A_{i+3,j}^n - A_{i-3,j}^n)}{24\Delta x}, \qquad (22) $$

and with odd ($i, j$) values for Eq. (20) ($A_{i,j} = W_{i,j}$) and even ones for Eq. (21) ($A_{i,j} = S_{i,j}$). Derivatives in the y-direction are obtained by the same finite difference scheme with differences in j-index instead of i-index. Coefficients of material stiffness appearing in front of the spatial derivatives of velocity in the constitutive equations are evaluated at the corresponding mesh point with odd ($i, j$) values, while the material densities appearing in the equations of momentum balance are evaluated at mesh points with even ($i, j$) values. Time derivatives are evaluated by a second order, centered finite difference

$$ \left(\frac{\partial S}{\partial t}\right)_{i,j}^n = \frac{S_{i,j}^{n+1/2} - S_{i,j}^{n-1/2}}{\Delta t} + O(\Delta t^2), \qquad (23) $$

$$ \left(\frac{\partial W}{\partial t}\right)_{i,j}^{n+1/2} = \frac{W_{i,j}^{n+1} - W_{i,j}^{n}}{\Delta t} + O(\Delta t^2). \qquad (24) $$

The reader is referred to Fig. 2 in the work of Lisitsa and Vishnevskiy (2010) for a visualization of the grid-staggering. Thanks to this staggering, left-hand and right-hand sides of the constitutive equations (10)–(12) and equations of motion (17) and (18) are computed on the adequate time and space grid points. To ensure stability, time and spatial steps are linked by the value of the Courant-Friedrichs-Lewy (CFL) constant CFL = $c_{max}\Delta t/\Delta x = 0.33$ (computed on the highest sound velocity $c_{max} = \sqrt{C_{11}/\rho} \approx 13\,000\,\text{m/s}$ in the axial direction of fibers.). The mesh size $\Delta x$ is chosen equal to $8\,\mu\text{m}$. For a pressure signal of halfwidth 15 ns and therefore of frequency 66 MHz (about 7 times the signal peak frequency), propagating at the speed of shear waves in fibers $\sqrt{G_x/\rho} = 3750\,\text{m/s}$, the corresponding wavelength is about $57\,\mu\text{m}$, thus ensuring a discretization of seven points per wavelength, a sufficient value for a fourth-order scheme in space. This mesh size is four times the size of the pixel of the images, and choosing this instead of one pixel reduces by almost two orders of magnitude the computational time. Geometrically, it is much smaller than the thickness of the fibers, which is recalled to be millimetric. Hence, we expect the detailed geometry of the medium to be adequately discretized. The resulting time step $\Delta t$ equals about 0.2 ns. For the initial conditions, as laser is an impulsive source and because of the principle of causality, we assume the laser illumination starts only at initial time $t_n > \Delta t/2$, so that all velocity components ($u, v$) are set to zero at time step $n = 0$, and so are all stress components ($\sigma_{xx}, \sigma_{yy}, \sigma_{xy}$) at time step $n = 1/2$. For a simulation capturing three successive arrivals of the shock front at the rear surface of a 2-mm thick
sample, the computation time on a laptop with an i5-6300HQ CPU at 2.30 GHz is approximately 45 min.

For inner boundaries between various material components, FDTD simulations are known to introduce some artefacts, as regular mesh grid points do not match the physical interfaces. One way to deal with this issue is to introduce an immersed boundary treatment [see Zhang and LeVeque (1997) for wave equation, Lombard and Piraux (2004) for elastodynamics, or Gao et al. (2015) for application to Virieux's scheme in geophysics]. However, here the inner interfaces are extremely numerous and complex. For the sake of simplicity and numerical efficiency, we preferred to maintain Lebedev's algorithm unchanged and smooth the numerical outputs by spatial and temporal averaging (see Sec. V for more details). This choice is also coherent with the fact that VISAR experiments also induce some averaging, and because comparison of 2D simulations with 3D experiments can be done only statistically. In the homogeneous and isotropic case, Lebedev's algorithm has been validated by comparison with outputs of Virieux's scheme (Cuenca et al., 2020) applied to laser shock wave generated in aluminum, with perfect agreement between the two methods.

## V. RESULTS AND DISCUSSION

Analysis of experimental data and their comparison with simulations have to be performed from a statistical point of view because (i) the detailed arrangement of fibers and epoxy varies from sample to sample, and for a given sample, from shooting point to shooting point, (ii) 2D cuts used as geometrical inputs of simulations do not fully reflect the 3D actual local material structure beneath each shooting point. Therefore, for each sample thickness, we compare the VISAR measured velocities at the opposite face for the number of successful laser shots (reported in Table II and varying between 3 and 8), to the simulations performed using the eight local geometries extracted below the eight planed shots. For both simulations (in red) and measurement (in blue), we report on Fig. 4 the mean value of the temporal velocity waveform, and the 95% probability computed by assuming a Student's t-distribution (Arfken et al., 2013) as usual for samples of small size. To be able to perform statistics and to compare the data, the velocity waveforms have been normalized in time. The time chosen for normalization is the average time of arrival $\langle t_a \rangle$ of the first shock measured for a given sample, the time of arrival $t_a$ being defined as the time delay between the laser emission and the first time for which the measured (or simulated) velocity reaches 50% of its absolute peak. This also allows one to determine the first shock speed of sound $c_F = d/t_a$ for a given sample of thickness $d$, and the average one $\langle c_F \rangle = d/\langle t_a \rangle$. A comparison of the measured and simulated values of this shock speed is presented in Fig. 4(f). For confidentiality reasons again, all shock speeds are normalized by the average value $c_{ref}$ of all 29 measurements. For simulated signals, amplitudes have

been adjusted manually, also to allow a statistical comparison. The issue of amplitudes is discussed later on.

VISAR has a spatial resolution of 100 µm and a temporal one of 3 ns. For a comparison with VISAR results, and for a spatial grid size $\Delta x$ of about 8 µm, we therefore need to average the values of the simulated velocity over 13 spatial grid points, and over 15 temporal ones. However, because of (i) the “pixelization” of the various boundaries within the samples, (ii) the absence of particular treatment of these interfaces such as immersed boundary method, (iii) the noise in the evaluation of the local orientation of fibers, the outputs of the FDTD simulations are quite noisy, even after this averaging. To reach a noise level similar to the measured one, each simulation is smoothed spatially over 100 spatial grid points instead of 13 (but still over only 15 temporal ones).

Figure 4 shows that, despite the complex structure of the material, sharp wavefronts are still observable after the shock has travelled over one sample thickness, independent of its length. Multiple scattering results into an attenuation measurable as a finite shock rise time (time needed for the velocity to go from 10% to 90% of the peak amplitude). Values of rise times averaged for each sample are reported in Table II. The simulated rise time in general increases with sample thickness and is comparable to (for samples 2, 3, and 5) or larger than (for samples 3 and 4) the measured one. Larger simulated rise times may be due to the artificial numerical dispersion, and to noise induced by the lack of precise description of the numerous interfaces between various components. Also, clearly observable is the amplitude decay of the post-shock wave profile, which is well predicted by simulations. In particular, a second shock arrival (at around dimensionless time 3) is clearly visible for the two thinner samples. However, this arrival is not as sharp as the first one, contrarily to what is seen in homogeneous metals (Cuenca et al., 2020). As no intrinsic absorption has been introduced in the model, this means that structural attenuation simply by the scattering of waves at the different material heterogeneities explains most of the overall attenuation. This is an important result, pointing out that only linear elastic parameters of the Tisse3D material are needed to predict most of its behaviour under laser shock exposure. For other, thicker samples, the second arrival appears only as a weak oscillation. Their similar amplitudes when comparing simulations and observations again confirm the dominant role of scattering. For measurement, variability of measured waveforms from average ones keeps relatively small for samples 2, 3, and 4 but is more important for samples 1 and 5. This may be due to an insufficient number of successful shots for sample 1, and also to longer propagation distances for sample 5. Deviations of 2D simulations are within a ±10 m/s range in all cases. This suggests that a number of eight cases is necessary for a reasonable statistical convergence.

Figure 4 (bottom right subfigure) also compares the normalized measured and simulated first shock velocities $\langle c_F \rangle$ along with their 95% probability. Considering all samples
<figure><img src="image_5.png" /><figcaption>FIG. 4. Mean value (thick line) and 95% probability (colored areas) of VISAR measurement (in blue) and FDTD simulations (in red) of measured opposite face velocity of samples of thicknesses (from left to right and top to bottom) 2, 3, 3.7, 4, and 7 mm. For a given sample, time is normalized by the average of measured times arrival of first shock. Amplitudes of simulations are adjusted to fit measured ones. Right bottom figure: normalized velocity of the first shockwave versus sample thickness. Comparison between experiments (in blue) and simulations (in red), with mean value (circles) and 95% probability (error bars).</figcaption></figure>

and all 29 experimental values, the variations of $c_F$ are within an overall ± 4% around reference. For a given sample, these variations are within an overall ± 2% range. Simulated values are systematically lower than measured ones, with similar variability except for the thinnest sample. The too low values resulting from simulations can be explained by the 2D approximation, according to which the Tissue3D is made of three different types of constituents: epoxy, in-plane fibers and out of plane fibers. As the shock wave is clearly the first signal to arrive (there is no visible precursor), we can assume it propagates in each material

component along the direction y (the shortest source to receiver path) with the highest wave speed within this component. Otherwise, a signal would arrive before the shock, which is not what is observed. In epoxy (of average proportion 7%), the highest speed is the one of longitudinal waves. Within in-plane fibers (of average proportion 53%), the highest speed is the one of quasi-shear waves, polarized mainly in the x direction for small angles $\theta$ (see Appendix B). This speed increases with angle $\theta$, from 3748 m/s at 0° to 4348 m/s at 10° (approximately the mean value of fibers orientation). In out of plane fibers (of average proportion 40%), the highest
velocity is the speed of longitudinal, compression waves in the transverse section, equal to 3085 m/s (see Appendix B). However, the notion of out of plane fibers is somewhat artificial, resulting from the 2D simplifying approximation used for simulations: the orientation of out of plane fibers is in reality not perfectly perpendicular to the selected 2D plane, and taking into account this orientation would also increase the wave speed. This explains why the 2D simulations tend to underestimate the shock velocity. Regarding the variability of the shock velocity, it is similar between simulations and experiments, about ±2% for a given sample. Only for the thinnest 2 mm case do simulations show a significantly higher variability. This is due to one particular sample which has an epoxy spot much larger than the average in the xy plane.

Figure 5 displays the mean value and 95% probability of the measured peak velocity, compared to the simulated one. For this last one, we used for the amplitude of the pressure loading the value given by ESTHER simulations for aluminum. However, we accounted for the impedance mismatch between the aluminum adhesive tape and the composite layer. For this, we multiplied this amplitude by the transmission coefficient *T* between aluminum and composite material for a pure compression wave at normal incidence relative to the interface between the two materials. The impedance of the composite was evaluated by two methods giving very similar results. First, we relied on a standard mechanical homogenization (i.e., at low frequencies). Second, the impedance was determined as the product of the average density of the material times the average wave velocity $\langle c_F \rangle$ determined by the arrival time of the first shock, see Fig. 4(f). Both methods provided very similar values (differences of about 0.3%). Experimentally, we observe a huge variability of the peak velocity, which seems to be very sensitive to the detailed, local structure of the composite below the illumination point. This variability is maximum for the 3-mm-thick sample (eight successful

<figure><img src="image_6.png" /><figcaption>FIG. 5. Mean value (circles) and 95% probability (error bars) of VISAR measured (in blue) and FDTD simulated (in red) opposite face peak velocity of samples of thicknesses 2, 3, 3.7, 4, and 7 mm.</figcaption></figure>

shots) with a 150–250 m/s range, though this case showed the lowest variability of the shock wave velocity. The average values nevertheless tend to indicate an expected decrease in the peak value with thickness. These two trends (strong variability and decrease with sample thickness) are also observed for the numerical simulations, though the variability is much less pronounced except for the thinnest sample. Similar decrease with thickness once again indicates that scattering by the complex geometrical local structure of the material is the main source of signal attenuation. However, simulations tend to strongly underestimate the peak amplitude, except for the thinnest sample, which is an unexpected result. One would rather have expected a numerical overestimation, as the simulations do not take into account the impedance mismatches due to the adhesive layer of the aluminum tape. The influence of this one on the process of laser shock generation is investigated in detail in Ayad et al. (2023). Though this underestimation is still an open question, we have several possible explanations that would need further investigations to elucidate it. Experimentally, VISAR measurement assumes a rather homogeneous velocity field over the surface of the detection laser spot, a hypothesis that may not be satisfied here. Moreover, 3D-diffraction effects have been shown to be important for shock wave propagation in a homogeneous material (Cuenca et al., 2020), and shock wave resulting from emission by a piston-like pressure excitation are well-known in fluids to induce sharp peaks (Canney et al., 2008). A similar behavior here could not be captured by 2D simulations. The numerical method may overestimate the signal absorption due to a too rough discretization of the interfaces between the materials resulting from the cartesian geometry of the FDTD method. Also, the FDTD numerical scheme is known to induce some errors by not considering the effective interfaces between the various materials. Methods such as immersed boundaries (Lombard et al., 2008) are likely to provide more precise results, but with a large increase in complexity in such 3D cases with multiple interfaces. At last, the process of averaging the raw numerical results over the free surface is also likely to introduce artificial attenuation. Indeed, the raw signal without averaging shows a very noisy peak value over the opposite free surface, with peak values much larger than the average ones. Also, the numerical and experimental averaging processes may not be similar, the experimental one being unexplored.

## VI. CONCLUSION

The propagation of laser shock waves in 3D woven composite materials is explored experimentally, by measuring the velocity of the face opposite to the laser illumination, and numerically by a simple FDTD 2D approach allowing a statistical analysis for five samples of various thicknesses and weaving patterns. Indeed, both experiments and simulations clearly show a strong variability of the signal, depending on the details of the geometry of the composite between the illumination spot and the measurement
point: the number and shapes of the inner interfaces, and the proportion of various components (especially for thin samples) vary significantly. This induces more or less large changes in the temporal shape of the detected signal, in its velocity and its amplitude. This observation of the effect of spot location agrees with previous numerical simulations (Tserpes and Kormpos, 2024). Therefore, the comparison between experiments and simulations is performed statistically. The good agreement between the observed and simulated mean amplitude decay of the successive arrivals indicates that the scattering due to the strong heterogeneity of the material is the dominant effect explaining the signal attenuation. The first arrival nevertheless always remains quite sharp, so that a shock wave velocity can be measured. Simulations tend to slightly underestimate it, and this is explained by the 2D approximation. On the contrary, the peak amplitude shows a very strong variability. This is also observed numerically, but the mean value of the peak is largely underestimated by simulations. This may be due to the experimental uncertainties, the limits of the numerical model, and the difficulties of properly modeling the influence of the thin adhesive aluminum layer necessary for the process of shock generation. Further studies should explore this influence, carry out experiments on more samples and extend the numerical model to 3D. For this last point in particular, the material segmentation will be performed from tomographic volumes in a fully automated procedure, using deep learning method (Mendoza et al., 2021).

## ACKNOWLEDGMENTS

The silver deposition was performed by Loïc Becerra, research engineer at Institut des NanoSciences de Paris (INSP), mixed research unit (UMR) 7588 from Sorbonne Université and Centre National de la Recherche Scientifique (CNRS). Cross-section images were obtained by Nicolas Chauzy, technician at Safran Tech. E.C., as a Ph.D. student, benefited from an Industrial Convention of Formation by Research (CIFRE).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors state that they have no conflicts to disclose.

## DATA AVAILABILITY

The data that support the findings of this study are available from Safran Tech. Restrictions apply to the availability of these data, which were used under license for this study. Data are available from the authors (contact Mathieu Ducousso, mathieu.ducousso@safrangroup.com) upon reasonable request and with the permission of Safran Tech.

## APPENDIX A: TILTED TRANSVERSE ISOTROPY

To establish Eq. (11), we follow the calculations detailed by Duveneck and Bakker (2011) but in the

simplified in-plane strain case. Rotation from the local frame of reference (LFR) (X, Y) with X-axis oriented along fiber axis to the computational frame of reference (CFR) (x, y) is

$$ \begin{bmatrix} e_x \\ e_y \end{bmatrix} = \begin{bmatrix} \cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \end{bmatrix} \cdot \begin{bmatrix} e_X \\ e_Y \end{bmatrix} = R(\theta) \cdot \begin{bmatrix} e_X \\ e_Y \end{bmatrix} \quad (A1) $$

with $R(\theta)^T = R(\theta)^{-1} = R(-\theta)$. One notes $\sigma = [\sigma_{xx} \ \sigma_{xy} \ \sigma_{yy}]$ (respectively, $\epsilon$) the stress (respectively, the strain) 2-by-2 symmetric tensor in the CFR and $\sigma'$ (respectively, $\epsilon'$) the symmetric stress (respectively, the strain) 2-by-2 tensor in the LFR. One has (Royer and Valier-Brasier, 2022)

$$ \sigma = R(\theta) \cdot \sigma' \cdot R(\theta)^T \quad \text{and} \quad \epsilon = R(\theta) \cdot \epsilon' \cdot R(\theta)^T. \quad (A2) $$

Expanding this and using Voigt 2D notation with $\sigma^{2D}$ the 1-by-3 stress vector and $\epsilon^{2D}$ the 1-by-3 strain vector, one obtains

$$ \sigma^{2D} = T(\theta) \cdot (\sigma^{2D})' \quad \text{and} \quad \epsilon^{2D} = D \cdot T(\theta) \cdot D^{-1}(\epsilon^{2D})', \quad (A3) $$

where $T(\theta)$ is the coordinate transformation matrix

$$ T(\theta) = \begin{bmatrix} \cos^2\theta & \sin^2\theta & -2\sin\theta\cos\theta \\ \sin^2\theta & \cos^2\theta & 2\sin\theta\cos\theta \\ \sin\theta\cos\theta & -\sin\theta\cos\theta & \cos^2\theta - \sin^2\theta \end{bmatrix} \quad (A4) $$

and D is equal to the unit matrix except for the coefficient $D_{33} = 2$. This D matrix emanates from the coefficient $d = 2$ appearing in the definition of the strain vector in Voigt notation for nondiagonal terms. In the LFR one has

$$ (\sigma^{2D})' = C_{\text{fib}}^{2D} \cdot (\epsilon^{2D})'. \quad (A5) $$

Therefore, in the CFR, the local stiffness matrix $C_{\text{loc}}$ is deduced from the material stiffness matrix $C_{\text{fib}}^{2D}$ using Eqs. (A3) and (A5),

$$ C_{\text{loc}} = T(\theta) \cdot D \cdot C_{\text{fib}} \cdot T^{-1}(\theta) \cdot D^{-1}. \quad (A6) $$

It is easy to check from Eq. (A4) that $T^{-1}(\theta) = T(-\theta)$.

## APPENDIX B: WAVE VELOCITY IN FIBERS

We search for the velocity $c$ of a wave propagating through the fibers in the sample thickness direction $y$ depending on the orientation of fibers $\theta$ assumed to be constant, each velocity or stress field being simply proportional to $f(t - y/c)$. Replacing this form of solution in the equations of motion, we get

$$ (\rho c^2 - C_{33}^{\text{loc}})(\rho c^2 - C_{22}^{\text{loc}}) - (C_{23}^{\text{loc}})^2 = 0, \quad (B1) $$

which has two solutions, a slow one $c_S$ and a fast one $c_F$,
<figure><img src="image_7.png" /><figcaption>FIG. 6. Wave velocity in direction y in fibers depending on their orientation for in-plane strain. Blue line: $c_S$. Red line: $c_F$. Black dotted line: $c_{TF}$.</figcaption></figure>

$$ \sqrt{\frac{C_{33}^{loc} + C_{22}^{loc} \pm \sqrt{(C_{33}^{loc} - C_{22}^{loc})^2 + 4(C_{23}^{loc})^2}}{2\rho}} \quad (B2) $$

with – sign for the slow velocity and + for the fast one. These solutions are displayed in Fig. 6. One has $c_S(\theta = 0°) = \sqrt{(\lambda_T + 2\mu_T)/\rho} = 3058 \text{ m/s}$ the velocity of compression waves in the transverse section of fibers, $c_S(\theta = 90°) = c_F(\theta = 0°) = c_{TF} = \sqrt{G_X/\rho} = 3748 \text{ m/s}$ the velocity of shear waves polarized along the axis of fibers X (black dashed line), and $c_F(\theta = 90°) = \sqrt{(E_X + 4\nu_X^2 K_T)\rho} = 11944 \text{ m/s}$ the velocity of compression waves in the axis of fibers. This one is much higher because of the high axial stiffness of fibers. For moderate angles, the fastest wave is therefore a quasi-shear wave polarized mostly in the direction x along the axis of fibers, and whose celerity increases with inclination of fibers $\theta$. At 10°, the increase in celerity is 16%.

Abramoff, M., Magalhes, P., and Ram, S. (2003). "Image processing with ImageJ." Biophotonics Int. 11, 36–42.

Arfken, G. B., Weber, H. J., and Harris, F. E. (2013). "Probability and statistics," in Mathematical Methods for Physicists (Academic Press, New York), Chap. 23, pp. 1125–1179.

Ayad, M., Ünaldi, S., Scius-Bertrand, M., Le Bras, C., Fayolle, B., and Berthe, L. (2023). "Dynamical modeling of bi-layer aluminium adhesive tape for laser shock applications," Opt. Laser Technol. 163, 109366.

Bardy, S., Aubert, B., Berthe, L., Combis, P., Hébert, D., Lescoute, E., Rullier, J.-L., and Videau, L. (2016). "Numerical study of laser ablation on aluminum for shock-wave applications: Development of a suitable model by comparison with recent experiments," Opt. Eng. 56(1), 011014.

Barker, L. M., and Hollenbach, R. E. (1972). "Laser interferometer for measuring high velocities of any reflecting surface," J. Appl. Phys. 43(11), 4669–4675.

Berthe, L., Arrigoni, M., Boustie, M., Cuq-Lelandais, J. P., Broussillou, C., Fabre, G., Jeandin, M., Guipont, V., and Nivard, M. (2011). "State-of-the-art laser adhesion test (LASAT)," Nondestr. Test. Eval. 26(3-4), 303–317.

Berthe, L., Fabbro, R., Peyre, P., Tollier, L., and Bartnicki, E. (1997). "Shock waves from a water-confined laser-generated plasma," J. Appl. Phys. 82, 2826–2832.

Bras, C., Rondepierre, A., Seddik, R., Scius-Bertrand, M., Rouchausse, Y., Videau, L., Fayolle, B., Gervais, M., Morin, L., Valadon, S., Ecault, R.,

Furfari, D., and Berthe, L. (2019). "Laser shock peening: Toward the use of pliable solid polymers for confinement," Metals 9, 793.

Canney, M. S., Bailey, M. R., Crum, L. A., Khokhlova, V. A., and Sapozhnikov, O. A. (2008). "Acoustic characterization of high intensity focused ultrasound fields: A combined measurement and modeling approach," J. Acoust. Soc. Am. 124(4), 2406–2420.

Carvelli, V., and Poggi, C. (2001). "A homogenization procedure for the numerical analysis of woven fabric composites," Compos. Part A 32(10), 1425–1432.

Chen, X., Taylor, L. W., and Tsai, L.-J. (2011). "An overview on fabrication of three-dimensional woven textile preforms for composites," Text. Res. J. 81(9), 932–944.

Colombier, J.-P., Combis, P., Bonneau, F., Le Harzic, R., and Audouard, E. (2005). "Hydrodynamic simulations of metal ablation by femtosecond laser irradiation," Phys. Rev. B 71, 165406.

Correa, C., Peral, D., Porro, J. A., Díaz, M., Ruiz de Lara, L., García-Beltrán, A., and Ocaña, J. L. (2015). "Random-type scanning patterns in laser shock peening without absorbing coating in 2024-T351 Al alloy: A solution to reduce residual stress anisotropy," Opt. Laser Technol. 73, 179–187.

Cuenca, E., Ducousso, M., Rondepierre, A., Videau, L., Cuvillier, N., Berthe, L., and Coulouvrat, F. (2020). "Propagation of laser-generated shock waves in metals: 3D axisymmetric simulations compared to experiments," J. Appl. Phys. 128(24), 244903.

Ding, H., Chen, W., and Zhang, L. (2006). Elasticity of Transversely Isotropic Materials (Springer Science & Business Media, New York), pp. 1–27.

Dowling, N. E. (1999). Mechanical Behavior of Materials: Engineering Methods for Deformation, Fracture, and Fatigue, 4th ed. (Pearson Education, London, UK), pp. 215–224.

Ducousso, M., Bardy, S., Rouchausse, Y., Bergara, T., Jenson, F., Berthe, L., and Cuvillier, N. (2018). "Quantitative evaluation of the mechanical strength of titanium/composite bonding using laser-generated shock waves," Appl. Phys. Lett. 112(11), 111904.

Duveneck, E., and Bakker, P. M. (2011). "Stable P-wave modeling for reverse-time migration in tilted TI media," Geophysics 76(2), S65–S75.

Fabbro, R., Peyre, P., Berthe, L., and Scherpereel, X. (1998). "Physics and applications of laser-shock processing," J. Laser Appl. 10(6), 265–279.

Fiji (2025). "Fiji," https://fiji.sc/ (Last viewed April 10, 2025).

Fleming, K. J., and Broyles, T. A. (2003). "Shock analysis using the multi point velocimeter (VISAR)," Report No. SAND2003-3759, Sandia National Laboratories (SNL), Albuquerque, NM, and Livermore, CA.

Gao, L., Brossier, R., Pajot, B., Tago, J., and Virieux, J. (2015). "An immersed free-surface boundary treatment for seismic wave simulation," Geophysics 80, T193–T209.

Ghosh, R., and De, S. (2014). "Z-fiber influence on high speed penetration of 3D orthogonal woven fiber composites," Mech. Mater. 68, 147–163.

Khayyam, H., Jazar, R. N., Nunna, S., Golkarnarenji, G., Badii, K., Fakhrhoseini, S. M., Kumar, S., and Naebe, M. (2020). "PAN precursor fabrication, applications and thermal stabilization process in carbon fiber production: Experimental and mathematical modelling," Prog. Mater. Sci. 107, 100575.

Lebedev, V. I. (1964). "Difference analogues of orthogonal decompositions, basic differential operators and some boundary problems of mathematical physics. I," USSR Comput. Math. Phys. 4(3), 69–92.

Lisitsa, V., and Vishnevskiy, D. (2010). "Lebedev scheme for the numerical simulation of wave propagation in 3D anisotropic elasticity," Geophys. Prospect. 58(4), 619–635.

Lombard, B., and Piraux, J. (2004). "Numerical treatment of two-dimensional interfaces for acoustic and elastic waves," J. Comput. Phys. 195(1), 90–116.

Lombard, B., Piraux, J., Gélis, C., and Virieux, J. (2008). "Free and smooth boundaries in 2-D finite-difference schemes for transient elastic waves," Geophys. J. Int. 172(1), 252–261.

Mendoza, A., Schneider, J., Parra-Denis, E., Obert, E., and Roux, S. (2019). "Differentiating 3D textile composites: A novel field of application for digital volume correlation," Compos. Struct. 208, 735–743.

Mendoza, A., Trullo, R., and Wielhorski, Y. (2021). "Descriptive modeling of textiles using FE simulations and deep learning," Compos. Sci. Technol. 213, 108897.

Mendoza Quispe, A. (2019). "Image correlation for textile descriptors," Ph.D. thesis, Université Paris-Saclay, Saclay, https://theses.hal.science/tel-02343183/ (Last viewed April 10, 2025).
Nayfeh, A. H. (1991). "The general problem of elastic wave propagation in multilayered anisotropic media," *J. Acoust. Soc. Am.* **89**(4), 1521-1531.

Operto, S., Virieux, J., Ribodetti, A., and Anderson, J. (2009). "Finite-difference frequency-domain modeling of viscoacoustic wave propagation in 2D tilted transversely isotropic (TTI) media," *Geophysics* **74**, T75–T95.

Pankow, M., Waas, A., Yen, C.-F., and Ghiorse, S. (2011). "Shock loading of 3D woven composites: A validated finite element investigation," *Compos. Struct.* **93**(5), 1347–1362.

Rokhlin, S., and Wang, L. (2002). "Ultrasonic waves in layered anisotropic media: Characterization of multidirectional composites," *Int. J. Solids Struct.* **39**(21–22), 5529–5545.

Royer, D., and Valier-Brasier, T. (2022). *Elastic Waves in Solids 1 Propagation* (Wiley, New York), pp. 32–38, 266–270.

Saenger, E. H., Gold, N., and Shapiro, S. A. (2000). "Modeling the propagation of elastic waves using a modified finite-difference grid," *Wave Motion* **31**(1), 77–92.

Schindelin, J., Arganda-Carreras, I., Frise, E., Kaynig, V., Longair, M., Pietzsch, T., Preibisch, S., Rueden, C., Saalfeld, S., Schmid, B., Tinevez, J.-Y., White, D. J., Hartenstein, V., Eliceiri, K. W., Tomançak, P., and Cardona, A. (2012). "Fiji: An open-source platform for biological-image analysis," *Nat. Methods* **9**(7), 676–682.

Scius-Bertrand, M., Videau, L., Rondepierre, A., Lescoute, E., Rouchausse, Y., Kaufman, J., Rostohar, D., Brajer, J., and Berthe, L. (2021). "Laser induced plasma characterization in direct and water confined regimes:

New advances in experimental studies and numerical modelling,” *J. Phys. D: Appl. Phys.* **54**(5), 055204.

Shirasu, K., Goto, K., and Naito, K. (2020). “Microstructure-elastic property relationships in carbon fibers: A nanoindentation study,” *Compos. Part B 200*, 108342–108346.

Su, H., Si, X., Liu, Y., Xu, M-m., Huang, G.-Y., and Pan, J. (2023). "Experimental and numerical investigation of the behavior of three-dimensional orthogonal woven composite plates under high-velocity impact," *Mech. Adv. Mater. Struct.* **30**(16), 3293–3302.

Tesoro, G. (1988). *Epoxy Resins-Chemistry and Technology* (Wiley Online Library, New York).

Thierry, V., Mesnil, O., and Chronopoulos, D. (2020). “Experimental and numerical determination of the wave dispersion characteristics of complex 3D woven composites,” *Ultrasonics* **103**, 106068.

Tserpes, K., and Kormpos, P. (2024). "Detailed finite element models for the simulation of the laser shock wave response of 3D woven composites," *J. Compos. Sci.* **8**(3), 83.

Virieux, J. (1986). "P-SV wave propagation in heterogeneous media: Velocity-stress finite-difference method," *Geophysics* **51**, 889–901.

Wielhorski, Y., Mendoza, A., Rubino, M., and Roux, S. (2022). "Numerical modeling of 3D woven composite reinforcements: A review," *Compos. Part A 154*, 106729.

Zhang, C., and LeVeque, R. J. (1997). "The immersed interface method for acoustic wave equations with discontinuous coefficients," *Wave Motion* **25**(3), 237–263.