Development of a Numerical Model to Simulate Laser-Shock Paint Stripping on Aluminum Substrates

: An explicit 3D Finite Element (FE) model was developed in the LS-Dyna code to simulate the laser shock paint stripping on aircraft aluminum substrates. The main objective of the model is to explain the physical mechanisms of the laser shock stripping process in terms of shock wave propagation, stress and strain evolution and stripping shape and size and to evaluate the effects of laser and material parameters on the stripping pattern. To simulate the behavior of aluminum, the Johnson–Cook plasticity model and the Gruneisen equation of state were applied. To simulate stripping, the cohesive zone modeling method was applied. The FE model was compared successfully against experiments in terms of back-face velocity proﬁles. The parameters considered in the study are the aluminum thickness, the epoxy paint thickness, the laser spot diameter, the fracture toughness of the aluminum/epoxy interface and the maximum applied pressure. In all cases, a circular solid or hollow stripping pattern was predicted, which agrees with the experimental ﬁndings. All parameters were found to affect the stripping pattern. The numerical results could be used for the design of selective laser shock stripping tests.


Introduction
In recent years, laser-induced shock waves were used in a wide range of engineering and materials science applications including the enhancement of fatigue strength [1] of metallic materials (laser peening) [2][3][4][5], the evaluation of adhesion strength (laser shock adhesion test) in adhesive joints [6][7][8] and the selective stripping of external coatings [9] in aerospace structures. Specifically, selective stripping is not a new process in aeronautical engineering, and it has a crucial role in the maintenance and restoration of aircraft during their lifetimes [9]. Usually, applied is a hazardous chemical stripping process with an important environmental impact followed by a plastic media or the eco-friendly walnut shell blasting, which damages the substrate material's surface. The main challenges in the selective laser stripping method are the high cost of the experimental set-up and the lack of standard experimental and numerical frameworks. This brings the need for the development of numerical tools to design process windows.
The laser shock selective stripping process can be proved to be a very useful ecofriendly tool, in comparison to media blasting and chemical stripping, that can be applied for the stripping of exterior primer coatings (top-coats) and for the stripping of structural primer coatings. Using this method, some parts of the aircraft that were stripped, without being damaged, can be recyclable and be used, after some proper post treatment, for any other application. At present, research on the laser shock selective stripping method aims to observe, understand and explain the physical phenomena that take place and upscale the method to an industrial level.
The development of an industrial laser stripping method aims to be a cost effective, eco-friendly method that restricts the damage in the coating area exclusively, through the optimization of the laser parameters. Because of the complicated physical phenomena that take place, any change in the laser parameters or material properties, such as inservice aging, affects the stripping result and thus, a comprehensive optimization of the method is needed. There are two ways to perform such an optimization: experimentally and numerically. Both can provide significant results through parametric studies, but in combination, they can provide a more detailed and complete description of the process. Using experimental tests, various stripping patterns can be obtained for different laser parameters and materials. Experiments can also provide important information that can be used for validation purposes by the numerical models. Using validated models, a detailed stress field in the material and the damage propagation for each timestep can be simulated. Numerical simulations have the potential to reduce the required number of tests for optimization but require appropriate input parameters that are sometimes difficult to be obtained.
Only a few works on the numerical simulation of the laser shock process were reported in the literature. Ivetic [2] numerically evaluated the effect of laser shock peening on thin sheets using a 3D explicit FE code developed in the Abaqus commercial software. Le Bras et al. [3] made an accurate pressure determination in the polymer confinement regime by coupling finite element simulations and experimental measurements. Peyre et al. [4], by finite element simulations, correctly reproduced the evolution of shock wave attenuation and elastic precursor amplitude and computed surface deformations induced by local impacts in Al-Cu-Li alloys. Ocana et al. [5] used the SHOCKLAS code to compute the through-the-thickness stresses in thin aluminum plates due to the laser shock processing treatment. Morales et al. [10] developed a simulation model (SHOCKLAS), dealing with the main aspects of laser shock processing in a coupled way. Using models, the authors simulated the hydrodynamic behavior of plasma expansion between the confinement layer and the material, using HELIOS (1D radiation-hydrodynamics Lagrangian fluid dynamic code). With this code, the plasma dynamics and the influence of the confining layer on the plasma pressure were studied. Bardy et al. [11] and Scius-Bertrand et al. [12], using ESTHER numerical code, describe the laser-matter interaction and thus predict the ablation pressure. Ecault et al. [13] developed a 3D FE model using the LS-DYNA FE code to study the laser shock and shear wave propagation induced by femtosecond laser irradiation in epoxy resins. The model was validated by experiments and used to understand the shock wave phenomenon. Bolis et al. [14], using HUGO and SHYLAC numerical codes, formed a comprehensive approach for the determination of mechanisms that take place in coating debonding and a quantitative evaluation of its fracture strength. Ecault et al. [6], in the first systematic effort, developed a reliable numerical model to be used for the simulation of the laser shock wave adhesion test of bonded composites. A calibration method was applied to set the model input parameters. The dynamic material parameters were identified from experimental results and validated through a complete campaign of laser shocks on various carbon fiber reinforced plastic (CFRP) materials (monolithic and bonded). The numerical results enabled the authors to understand the stress distribution within the composite assembly during wave propagation. Finally, Ecault et al. [6,7,13] used the models developed in [7] to perform an optimization of tunable pulse duration, double pulses and symmetrical laser shocks for adhesive bonding assessment. From the above literature overview, it becomes clear that there are no reported works on the simulation of laser shock stripping.
The aim of the present work is to develop a detailed 3D numerical model for the simulation of laser shock paint stripping on aluminum substrates. After being validated, it will be used to understand the basic physical mechanisms that take place between the aluminum and the epoxy materials during stripping, with thickness in the range of a few µm. Finally, a parametric study on the effects of laser and material parameters on the stripping pattern will be performed.

The Laser Shock Stripping Process
The laser shock stripping process is schematically described in 6 1. By focusing the high-power laser (1 J, 10 ns) onto a target, plasma creation occurs. With the help of plasma confinement (with water or solid), its duration becomes two times longer and four times higher [15]. The whole process can be distinguished into three parts. First, energy deposition from the laser pulse to the surface of the aluminum takes place and the creation of plasma through evaporation, ionization and expansion, begins.
The Gaia HP laser from THALES company (France) at Hephaistos, PIMM Lab, Paris, is used during the experiments. The Gaia HP laser is a Nd:YAG laser with 7.2 ns pulse duration, 14 Joules of energy (Gaussian pulse) and 2 Hz of repetition rate at 532 nm wavelength. Single beam is used, and energy calibration is done via a calorimeter (QE50LP-H-MB-QED, Gentec, Quebec, QC, Canada) before starting every experiment. For the conducted experiments, the spot size is kept constant as 4 mm through a lens with focal length of 198 mm. In addition, diffractive optical element (DOE) is used to obtain a top-hat shaped beam with an equal beam distribution along the spot size [16].
This high-pressure plasma-surface interaction initiates a shock wave that propagates to each material, according to its acoustic impedance, which is a function of density and the Hugoniot dimensionless coefficient S. The shock wave can be imagined as a propagating surface at which the displacement is continuous, and the mass density, temperature, stress and other variables are discontinuous. Interactions of the waves and reflections develop a complicated stress field composed of tensile and compressive stress faces that propagate inside the material. Stripping occurs due to the development of tensile stresses at the aluminum/epoxy interface ( Figure 1). μm. Finally, a parametric study on the effects of laser and material parameters on the stripping pattern will be performed.

The Laser Shock Stripping Process
The laser shock stripping process is schematically described in 6 1. By focusing the high-power laser (1 J, 10 ns) onto a target, plasma creation occurs. With the help of plasma confinement (with water or solid), its duration becomes two times longer and four times higher [15]. The whole process can be distinguished into three parts. First, energy deposition from the laser pulse to the surface of the aluminum takes place and the creation of plasma through evaporation, ionization and expansion, begins.
The Gaia HP laser from THALES company (France) at Hephaistos, PIMM Lab, Paris, is used during the experiments. The Gaia HP laser is a Nd:YAG laser with 7.2 ns pulse duration, 14 Joules of energy (Gaussian pulse) and 2 Hz of repetition rate at 532 nm wavelength. Single beam is used, and energy calibration is done via a calorimeter (QE50LP-H-MB-QED, Gentec, Quebec, QC, Canada) before starting every experiment. For the conducted experiments, the spot size is kept constant as 4 mm through a lens with focal length of 198 mm. In addition, diffractive optical element (DOE) is used to obtain a top-hat shaped beam with an equal beam distribution along the spot size [16].
This high-pressure plasma-surface interaction initiates a shock wave that propagates to each material, according to its acoustic impedance, which is a function of density and the Hugoniot dimensionless coefficient S. The shock wave can be imagined as a propagating surface at which the displacement is continuous, and the mass density, temperature, stress and other variables are discontinuous. Interactions of the waves and reflections develop a complicated stress field composed of tensile and compressive stress faces that propagate inside the material. Stripping occurs due to the development of tensile stresses at the aluminum/epoxy interface ( Figure 1). The shock wave propagation into the aluminum/epoxy system is described in the sketch in Figure 2. When the pressure from plasma is applied to the material, an elastic precursor shock (i) followed by an elastic-plastic compression shock (ii), are formed and propagate into the material [17]. The plastic-decompression shock (iii) and elastic-plastic decompression wave (iv) begin to propagate into the material after the pressure is removed. The desired interaction takes place between the elastic precursor (i) and the plastic decompression shock (iii). Usually in the literature, the precursor and the compression shock are referred as the "compression shock wave" and the decompression shock wave as the "release wave". The shock wave propagation into the aluminum/epoxy system is described in the sketch in Figure 2. When the pressure from plasma is applied to the material, an elastic precursor shock (i) followed by an elastic-plastic compression shock (ii), are formed and propagate into the material [17]. The plastic-decompression shock (iii) and elastic-plastic decompression wave (iv) begin to propagate into the material after the pressure is removed. The desired interaction takes place between the elastic precursor (i) and the plastic decompression shock (iii). Usually in the literature, the precursor and the compression shock are referred as the "compression shock wave" and the decompression shock wave as the "release wave".
In most cases of strong shocks, the reflections from the back free surface are irregular [18]. A schematic example of a regular reflection and a simple irregular reflection (Mach reflection) are presented in Figure 3. This small irregular area between reflected shock and reflection surface is called Mach stem and is characterized by high compression stress fields. As will be shown later in this paper, the Mach stem plays a significant role on the pattern of the developed stress field in the materials during the laser shock wave propagation. In most cases of strong shocks, the reflections from the back free surface are irregular [18]. A schematic example of a regular reflection and a simple irregular reflection (Mach reflection) are presented in Figure 3. This small irregular area between reflected shock and reflection surface is called Mach stem and is characterized by high compression stress fields. As will be shown later in this paper, the Mach stem plays a significant role on the pattern of the developed stress field in the materials during the laser shock wave propagation.

Specimen and Materials
The target specimen modeled consists of an aluminum AA 2024-T3 substrate, whose surface was treated by chemical etching, and a structural primer made from the epoxy CA7049 material by RESCOLL, Pessac, France. The primer's role is the adherence conservation between the substrate and the top-coat layers. A schematic of the specimen is given in Figure 1. The dimensions of the specimen are 80 mm × 125 mm. The thickness of the aluminum substrate is 0.97 mm and that of the epoxy primer is 0.025 mm.

Calculation of Pressure's Profile
The pressure profile that was modeled is shown in Figure 4 [9]. The maximum pressure value Pmax (GPa) in the profile was calculated by [19,20]   In most cases of strong shocks, the reflections from the back free surface are irregular [18]. A schematic example of a regular reflection and a simple irregular reflection (Mach reflection) are presented in Figure 3. This small irregular area between reflected shock and reflection surface is called Mach stem and is characterized by high compression stress fields. As will be shown later in this paper, the Mach stem plays a significant role on the pattern of the developed stress field in the materials during the laser shock wave propagation.

Specimen and Materials
The target specimen modeled consists of an aluminum AA 2024-T3 substrate, whose surface was treated by chemical etching, and a structural primer made from the epoxy CA7049 material by RESCOLL, Pessac, France. The primer's role is the adherence conservation between the substrate and the top-coat layers. A schematic of the specimen is given in Figure 1. The dimensions of the specimen are 80 mm × 125 mm. The thickness of the aluminum substrate is 0.97 mm and that of the epoxy primer is 0.025 mm.

Calculation of Pressure's Profile
The pressure profile that was modeled is shown in Figure 4 [9]. The maximum pressure value Pmax (GPa) in the profile was calculated by [19,20]

Specimen and Materials
The target specimen modeled consists of an aluminum AA 2024-T3 substrate, whose surface was treated by chemical etching, and a structural primer made from the epoxy CA7049 material by RESCOLL, Pessac, France. The primer's role is the adherence conservation between the substrate and the top-coat layers. A schematic of the specimen is given in Figure 1. The dimensions of the specimen are 80 mm × 125 mm. The thickness of the aluminum substrate is 0.97 mm and that of the epoxy primer is 0.025 mm.

Calculation of Pressure's Profile
The pressure profile that was modeled is shown in Figure 4 [9]. The maximum pressure value P max (GPa) in the profile was calculated by [19,20]  where, I 0 (GW/cm 3 ) is the laser's intensity, α is the part of the energy being used for the gas ionization, Z (g cm −2 /s −1 ) is the material's acoustic impedance, C 0 is the sound speed inside the material and u material's velocity.
where, I0 (GW/cm 3 ) is the laser's intensity, α is the part of the energy being used for the gas ionization, Z (g cm −2 /s −1 ) is the material's acoustic impedance, C0 is the sound speed inside the material and u material's velocity.

Constitutive Material Models
To simulate the high strain rate behavior of the aluminum due to the shock wave loading, the Johnson-Cook plasticity model [21] and the Gruneisen equation of state [21] were used. Thermal effects have not been simulated. The flow-stress expression of the material model is where ̄ is the effective plastic strain, * the normalized effective total strain-rate (VP.EQ.0), m the thermal softening and T* is the homologous temperature given by * = − − where − = internal energy (7) where Cp is aluminum's specific heat and V0 is the initial volume. Strain at fracture is calculated by where σ* is the ratio of pressure divided by the effective stress. The values of the parameters appearing in the above equations are listed in Table 1.

Constitutive Material Models
To simulate the high strain rate behavior of the aluminum due to the shock wave loading, the Johnson-Cook plasticity model [21] and the Gruneisen equation of state [21] were used. Thermal effects have not been simulated. The flow-stress expression of the material model is where ε p is the effective plastic strain, . ε * the normalized effective total strain-rate (VP.EQ.0), m the thermal softening and T* is the homologous temperature given by where where C p is aluminum's specific heat and V 0 is the initial volume. Strain at fracture is calculated by where σ* is the ratio of pressure divided by the effective stress. The values of the parameters appearing in the above equations are listed in Table 1.

Input Parameter Value
Input Damage Parameters

Stripping Simulation by Cohesive Zone Modeling
Failure of the materials' interface (stripping) was simulated using the cohesive zone modeling (CZM) method. The CZM method does not model any material but only the relation between the traction and separation forces and the critical energy release rates by the following equations where G IC is the critical energy release rate for mode I load, G IIC is the critical energy release rate for mode II load, T the peak traction force in the normal direction, S the peak traction force in the tangential direction, UND the ultimate displacement in the normal direction and UTD is the ultimate displacement in the tangential direction. The relation between traction and displacement (separation) was modeled using a bi-linear mixed-mode I + II traction-separation law schematically described in Figure 5. According to the law, the interface follows an elastic behavior until a specific stress value and after that there is a degradation of stiffness until final failure (complete debonding).

Stripping Simulation by Cohesive Zone Modeling
Failure of the materials' interface (stripping) was simulated using the cohesive zone modeling (CZM) method. The CZM method does not model any material but only the relation between the traction and separation forces and the critical energy release rates by the following equations where GIC is the critical energy release rate for mode I load, GIIC is the critical energy release rate for mode II load, T the peak traction force in the normal direction, S the peak traction force in the tangential direction, UND the ultimate displacement in the normal direction and UTD is the ultimate displacement in the tangential direction. The relation between traction and displacement (separation) was modeled using a bi-linear mixed-mode I + II traction-separation law schematically described in Figure 5. According to the law, the interface follows an elastic behavior until a specific stress value and after that there is a degradation of stiffness until final failure (complete debonding).  Due to the lack of data for the fracture toughness properties of the aluminum/epoxy interface, for which specific double-cantilever beam and end-notch flexure tests are required, properties for a CFRP/adhesive interface [23] were used. These properties are: G IC = 1018.52 J/m 2 and G IIC = 783.41 J/m 2 .

Finite Element Modeling
A 3D FE model of the aluminum/epoxy specimen was developed using the explicit FE software LS-DYNA. A mapped FE mesh consisting of different areas was applied. In the center of the specimen, a 4 mm diameter circular area, representing the laser spot, was meshed using elements with a size of 0.027 mm. Away from that area, elements with a size of 0.042 mm were used. The thickness of all elements is 0.005 mm. For the aluminum and epoxy part of the specimen, 3D solid elements with one integration point (ELFORM = 1) were used. The aluminum/epoxy interface (cohesive area) was simulated using zero thickness 8-noded, 4-point cohesive elements (ELFORM = 19). The FE meshes of the complete specimen and the individual parts are shown in Figure 6.

Finite Element Modeling
A 3D FE model of the aluminum/epoxy specimen was developed using the explicit FE software LS-DYNA. A mapped FE mesh consisting of different areas was applied. In the center of the specimen, a 4 mm diameter circular area, representing the laser spot, was meshed using elements with a size of 0.027 mm. Away from that area, elements with a size of 0.042 mm were used. The thickness of all elements is 0.005 mm. For the aluminum and epoxy part of the specimen, 3D solid elements with one integration point (ELFORM = 1) were used. The aluminum/epoxy interface (cohesive area) was simulated using zero thickness 8-noded, 4-point cohesive elements (ELFORM = 19). The FE meshes of the complete specimen and the individual parts are shown in Figure 6.  Table 2. The material model through which the CZM method were implemented in the LS-Dyna is the *MAT_138_COHESIVE_MIXED_MODE.
Using the FE model, second order objective stresses are computed and updated at every time-step. The timestep was set equal to 1 ns and the total analysis time was 2 μs. For the energy balance computation, the rigid-wall (or stone-wall) energy, the sliding in-  Table 2. The material model through which the CZM method were implemented in the LS-Dyna is the *MAT_138_COHESIVE_MIXED_MODE. Using the FE model, second order objective stresses are computed and updated at every time-step. The timestep was set equal to 1 ns and the total analysis time was 2 µs. For the energy balance computation, the rigid-wall (or stone-wall) energy, the sliding interface energy and the Rayleigh (or damping) energy dissipation were included. Furthermore, the hourglass energy, which refers to nonphysical distortions of the elements, is computed at the energy balance [21].

Parametric Study
The aim of the parametric study is to understand the influence on the stripping pattern (shape and size) of the aluminum thickness, the epoxy thickness, the spot diameter, the critical energy release rate values and the maximum applied pressure for a given laser intensity. Table 3 lists the variation range of each parameter. The thickness of the materials is a manufacturing parameter and can be modified. The laser spot diameter is a process parameter, which can be easily modified. The study of the deviation of G IC and G IIC is of increased importance due to fact that the actual values for the aluminum/epoxy interface are not available yet. Finally, the uncertainty in P max emanates from the methodology used to calculate P max from the laser intensity. The reference parameters used in the parametric study are listed in Table 3. Unless otherwise stated, the results presented hereafter will refer to the parameters of Table 3.

Model Validation
The numerical model was validated through the comparison of the computed backface velocity profiles (back face velocity vs. time) with experimental results obtained from [9]. In order to obtain material properties, velocity interferometer system for any reflector (VISAR) diagnostic is used which is based on the Michelson interferometer [25]. Since the target accelerates with effect of the applied shock wave, the Doppler shift occurs. Basically, with the VISAR diagnostic, change in Doppler shift of light which was reflected from the moving surface is measured. Moreover, the velocity of the target can be calculated, which can be linked to other material parameters [25]. In the frame of validation, a convergence study of the mesh density, especially through the thickness of the specimen, and the time-step was also performed. First, a comparison between experimental and numerical back face velocity graphs of pure aluminum 2024 for I = 3.09 GW/cm 2 , shown in Figure 7, were made. A good agreement between the first and second peak is observed. As it can be observed, there is a difference at the Hugoniot elastic limit (HEL) value between the two graphs. The numerical model underestimates this value, that is almost 100 m/s in comparison with the experimental, which is 150 m/s. This difference may occur due to the mechanical and thermodynamic properties of the aluminum, taken from previous works [2,22].      Figure 8 compares the predicted and experimental back face velocity profiles for three different laser intensities, namely, 2.32 GW/cm 2 (3072 MPa), 4.96 GW/cm 2 (3986 MPa) and 6.81 GW/cm 2 (4553 MPa). The experimental curves were created using VISAR measurements while the numerical curves were created by using the velocity values at the center node of the back free surface of the specimen, at the epoxy. In the analyses of Figure 8, stripping has not been simulated for the aluminum/epoxy specimen. As can be seen, there is a good agreement between the model and the tests in the sense that the numerical curve follows the experimental curve and there is a good convergence regarding the position and size of the velocity peaks. For all cases, unlike Figure 7, the HEL value is overestimated by the numerical models and is 200 m/s, while the experimental value is 150 m/s. This difference occurs due to lack of the exact mechanical properties of the epoxy primer and the interface between aluminum/epoxy. For the case of I = 2.32 GW/cm 2 (Figure 8a), the model predicts all velocity peaks but with an overestimation. For the cases of I = 4.96 GW/cm 2 ( Figure 8b) and I = 6.81 GW/cm 2 (Figure 8c), the model underestimates the initial secondary peak, predicts the first main peak very well and then it overestimates the forthcoming peaks. The three smaller peaks in Figure 8a-c are due to the main elastic-plastic shock, the elastic-plastic decompression shock and their reflections by the interface of aluminum/epoxy material. In Figure 8c, the small peak (5c) observed at the start of the experimental curve is a parasite signal of the laser pulse. Overall, it can be concluded that the model is capable of simulating the shock wave propagation inside the aluminum/epoxy specimen.

Stripping Simulation
For all cases, the predicted stripping pattern is circular, either solid or hollow, obviously due to the shape of the laser spot. This finding is in line with the experimental observations from [9]. Figure 9 depicts the predicted stripping pattern (gray-brown color) for the case of I = 1.75 GW/cm 2 and P max = 2400 MPa as well as the experimentally obtained stripping pattern. Both the tests and the preliminary numerical simulations have shown that the thickness of the ring in the stripping pattern and the transition from a hollow to a solid pattern and vice versa depend on the material and process parameters that are studied in the present work.
Aiming to understand the stripping process, a detailed stress analysis was performed [9]. The analysis was focused on the normal σ z stresses, which are responsible for the stripping. As the wave propagates, the compressive stress field tends to concentrate at the center, thus forming an outer ring of weaker compressive stresses. Next the compressive stress field at the outer ring turns into tensile, while at the center of the area the stresses remain compressive due to the development of a Mach stem by the irregular reflection of the shock front at the back free surface. When the Mach stem becomes weaker, a ring of tensile stress is developed at the center of the specimen. This stress propagation pattern at the surface of epoxy continues for as long as the shock propagates. to the mechanical and thermodynamic properties of the aluminum, taken from previous works [2,22].  Figure 8 compares the predicted and experimental back face velocity profiles for three different laser intensities, namely, 2.32 GW/cm 2 (3072 MPa), 4.96 GW/cm 2 (3986 MPa) and 6.81 GW/cm 2 (4553 MPa). The experimental curves were created using VISAR measurements while the numerical curves were created by using the velocity values at the center node of the back free surface of the specimen, at the epoxy. In the analyses of Figure 8, stripping has not been simulated for the aluminum/epoxy specimen. As can be seen, there is a good agreement between the model and the tests in the sense that the numerical curve follows the experimental curve and there is a good convergence regarding the position and size of the velocity peaks. For all cases, unlike Figure 7, the HEL value is overestimated by the numerical models and is 200 m/s, while the experimental value is 150 m/s. This difference occurs due to lack of the exact mechanical properties of the epoxy primer and the interface between aluminum/epoxy. For the case of I = 2.32 GW/cm 2 (Figure 8a), the model predicts all velocity peaks but with an overestimation. For the cases of I = 4.96 GW/cm 2 ( Figure 8b) and I = 6.81 GW/cm 2 (Figure 8c), the model underestimates the initial secondary peak, predicts the first main peak very well and then it overestimates the forthcoming peaks. The three smaller peaks in Figure 8a-c are due to the main elasticplastic shock, the elastic-plastic decompression shock and their reflections by the interface of aluminum/epoxy material. In Figure 8c, the small peak (5c) observed at the start of the experimental curve is a parasite signal of the laser pulse. Overall, it can be concluded that the model is capable of simulating the shock wave propagation inside the aluminum/epoxy specimen.

Stripping Simulation
For all cases, the predicted stripping pattern is circular, either solid or hollow, obviously due to the shape of the laser spot. This finding is in line with the experimental observations from [9]. Figure 9 depicts the predicted stripping pattern (gray-brown color) for the case of I = 1.75 GW/cm 2 and Pmax = 2400 MPa as well as the experimentally obtained stripping pattern. Both the tests and the preliminary numerical simulations have shown that the thickness of the ring in the stripping pattern and the transition from a hollow to a solid pattern and vice versa depend on the material and process parameters that are studied in the present work. Aiming to understand the stripping process, a detailed stress analysis was performed [9]. The analysis was focused on the normal σz stresses, which are responsible for the stripping. As the wave propagates, the compressive stress field tends to concentrate at the center, thus forming an outer ring of weaker compressive stresses. Next the compressive

Stripping Simulation
For all cases, the predicted stripping pattern is circular, either solid or hollow, obviously due to the shape of the laser spot. This finding is in line with the experimental observations from [9]. Figure 9 depicts the predicted stripping pattern (gray-brown color) for the case of I = 1.75 GW/cm 2 and Pmax = 2400 MPa as well as the experimentally obtained stripping pattern. Both the tests and the preliminary numerical simulations have shown that the thickness of the ring in the stripping pattern and the transition from a hollow to a solid pattern and vice versa depend on the material and process parameters that are studied in the present work. Aiming to understand the stripping process, a detailed stress analysis was performed [9]. The analysis was focused on the normal σz stresses, which are responsible for the stripping. As the wave propagates, the compressive stress field tends to concentrate at the center, thus forming an outer ring of weaker compressive stresses. Next the compressive 1b 2b 3b 4b 1c 2c 3c 4c 5c Figure 9. (a) A common stripping pattern predicted the model, (b) the experimental stripping pattern adapted from [9]. Figure 10 plots the evolution of σ z stress in the cohesive area for P max = 2400 MPa. In the plots, the stripped area is also shown with the gray-brown color (deleted elements). As shown, initially, the stripping has a four point axisymmetric shape (Figure 10a) and as it propagates the values of tensile stresses are decreasing (Figure 10b,c). As the stripping process approaches the end the compressive stress field covers the area around the stripped area and the tensile stresses increase (Figure 10d). A general observation is that the compressive stresses are acting as a stopping mechanism of the stripping. Figure 10 plots the evolution of σz stress in the cohesive area for max = 2400 MPa. In the plots, the stripped area is also shown with the gray-brown color (deleted elements). As shown, initially, the stripping has a four point axisymmetric shape (Figure 10a) and as it propagates the values of tensile stresses are decreasing (Figure 10b,c). As the stripping process approaches the end the compressive stress field covers the area around the stripped area and the tensile stresses increase (Figure 10d). A general observation is that the compressive stresses are acting as a stopping mechanism of the stripping.

Effect of Aluminum Thickness
With increasing aluminum thickness, a transition from the annular to the uniform stripping pattern takes place, as illustrated in Figure 11. This is due to the different shock propagation mechanism and the different distribution of the shock through the thickness stresses in the aluminum. The shock wave in a thinner specimen evolves faster and the Mach stem area dissolves quickly. Furthermore, the compressive field has lower values. Besides that, the highest tensile stress is closer to the epoxy and the tensile stress values

Effect of Aluminum Thickness
With increasing aluminum thickness, a transition from the annular to the uniform stripping pattern takes place, as illustrated in Figure 11. This is due to the different shock propagation mechanism and the different distribution of the shock through the thickness stresses in the aluminum. The shock wave in a thinner specimen evolves faster and the Mach stem area dissolves quickly. Furthermore, the compressive field has lower values. Besides that, the highest tensile stress is closer to the epoxy and the tensile stress values are lower. On the other hand, in a thicker specimen, the tensile stress field that develops after the Mach stem dissolution has a more uniform distribution, thus leading to more uniform stripping area. The differences of the shock wave propagation and the evolution of normal stresses between a thin and a thick aluminum substrate are schematically described by means of the comparison of the normal stress contours for the cases of t al = 0.7 mm and t al = 1.2 mm shown in Figure 12.
are lower. On the other hand, in a thicker specimen, the tensile stress field that develops after the Mach stem dissolution has a more uniform distribution, thus leading to more uniform stripping area. The differences of the shock wave propagation and the evolution of normal stresses between a thin and a thick aluminum substrate are schematically described by means of the comparison of the normal stress contours for the cases of tal = 0.7 mm and tal = 1.2 mm shown in Figure 12. are lower. On the other hand, in a thicker specimen, the tensile stress field that develops after the Mach stem dissolution has a more uniform distribution, thus leading to more uniform stripping area. The differences of the shock wave propagation and the evolution of normal stresses between a thin and a thick aluminum substrate are schematically described by means of the comparison of the normal stress contours for the cases of tal = 0.7 mm and tal = 1.2 mm shown in Figure 12.

Effect of Epoxy Thickness
With increasing epoxy thickness, a transition from a solid circular to an annular stripping pattern takes place ( Figure 13). For tep = 0.1 mm, a four-point axisymmetric damage

Effect of Epoxy Thickness
With increasing epoxy thickness, a transition from a solid circular to an annular stripping pattern takes place ( Figure 13). For t ep = 0.1 mm, a four-point axisymmetric damage area is predicted. This transformation of the stripping area is the opposite of the one that takes place with the increase of the aluminum thickness. For t ep = 125 µm and t ep = 150 µm, no stripping was predicted. This is attributed to the increased shock wave damping offered by the thicker epoxy and the lower effect of the Mach stem. For the latter values of t ep , an incomplete stripping is predicted.

Effect of Epoxy Thickness
With increasing epoxy thickness, a transition from a solid circular to an annular stripping pattern takes place ( Figure 13). For tep = 0.1 mm, a four-point axisymmetric damage area is predicted. This transformation of the stripping area is the opposite of the one that takes place with the increase of the aluminum thickness. For tep = 125 μm and tep = 150 μm, no stripping was predicted. This is attributed to the increased shock wave damping offered by the thicker epoxy and the lower effect of the Mach stem. For the latter values of tep, an incomplete stripping is predicted.  Figure 14 shows the simulated evolution of stripping for six different spot diameters dsp. For dsp = 2.5 mm and dsp = 3.0 mm, the stripping initiates as a solid circle and propagates with the same shape. This is due to the more uniform distribution of the laser's intensity. For dsp = 3.5 mm, the stripping initiates annularly, afterwards it also initiates at the center of the spot and ends as a solid circle. The same stands for dsp = 4.0 mm but with a different propagation rate. Finally, for dsp = 4.5 mm and dsp = 5.0 mm, the stripping initiates and propagates annularly; for these two cases, the stripping is incomplete.  Figure 14 shows the simulated evolution of stripping for six different spot diameters d sp . For d sp = 2.5 mm and d sp = 3.0 mm, the stripping initiates as a solid circle and propagates with the same shape. This is due to the more uniform distribution of the laser's intensity. For d sp = 3.5 mm, the stripping initiates annularly, afterwards it also initiates at the center of the spot and ends as a solid circle. The same stands for d sp = 4.0 mm but with a different propagation rate. Finally, for d sp = 4.5 mm and d sp = 5.0 mm, the stripping initiates and propagates annularly; for these two cases, the stripping is incomplete.  Figure 15 shows the stripping evolution and the final stripping patterns for the analyses with the different GIC and GIIC values described in Table 3. For the initial GIC, the stripping pattern starts as annular and turns into a solid. For GIC/1.25, a similar stripping process was predicted but with a larger final stripping area. For GIC/1.25, the stripping starts at the center of the loading spot, then annular stripping also starts, and both evolve into a  Figure 15 shows the stripping evolution and the final stripping patterns for the analyses with the different G IC and G IIC values described in Table 3. For the initial G IC , the stripping pattern starts as annular and turns into a solid. For G IC /1.25, a similar stripping process was predicted but with a larger final stripping area. For G IC /1.25, the stripping starts at the center of the loading spot, then annular stripping also starts, and both evolve into a final solid stripping area. The same stands for G IC /1.25 but with a higher evolution rate. For G IC /2.0, the stripping starts and evolves as a solid circle. The variation of the G IIC does not affect the stripping evolution and the final pattern as the stripping mechanism is mode-I dominated. Therefore, for the three cases with the different G IIC (bottom of Figure 15), the same stripping evolution as with the reference case was predicted.  Figure 15 shows the stripping evolution and the final stripping patterns for the analyses with the different GIC and GIIC values described in Table 3. For the initial GIC, the stripping pattern starts as annular and turns into a solid. For GIC/1.25, a similar stripping process was predicted but with a larger final stripping area. For GIC/1.25, the stripping starts at the center of the loading spot, then annular stripping also starts, and both evolve into a final solid stripping area. The same stands for GIC/1.25 but with a higher evolution rate. For GIC/2.0, the stripping starts and evolves as a solid circle. The variation of the GIIC does not affect the stripping evolution and the final pattern as the stripping mechanism is mode-I dominated. Therefore, for the three cases with the different GIIC (bottom of Figure  15), the same stripping evolution as with the reference case was predicted.

Effect of Maximum Applied Pressure
As there is an uncertainty in the correlation between the maximum applied pressure and the laser's intensity that applies in Equation (1), the effect of the peak pressure values on the stripping process were examined. Figure 16 depicts the final stripping patterns predicted for the different P values. As shown, by increasing P, the stripping pattern transforms from annular to solid and the stripping accumulates much quicker. As there is an uncertainty in the correlation between the maximum applied pressure and the laser's intensity that applies in Equation (1), the effect of the peak pressure values on the stripping process were examined. Figure 16 depicts the final stripping patterns predicted for the different P values. As shown, by increasing P, the stripping pattern transforms from annular to solid and the stripping accumulates much quicker.

Conclusions
In the present paper, a numerical model was developed to simulate the laser shock paint stripping on aluminum substrates. The main objectives of the model were to explain the physical mechanisms of the laser shock stripping process in terms of shock wave prop-

Conclusions
In the present paper, a numerical model was developed to simulate the laser shock paint stripping on aluminum substrates. The main objectives of the model were to explain the physical mechanisms of the laser shock stripping process in terms of shock wave propagation, stress and strain evolution and stripping shape and size and to evaluate the effects of laser and material parameters on the stripping pattern.
The main findings of the study are summarized as follows: • The model is capable of efficiently simulating the laser shock stripping process as indicated by the comparison of numerical and experimental results. • By increasing the aluminum thickness, a transition from the annular to the solid stripping pattern takes place. For values of aluminum thickness smaller than 1 mm, an incomplete stripping is predicted. • By increasing the epoxy thickness, a transition from a solid circular to an annular stripping pattern takes place. For 0.075 mm and 0.1 mm epoxy thickness, an incomplete stripping is predicted.

•
The laser spot's diameter significantly affects the stripping propagation and the final stripping pattern. For values of d sp larger than 4.0 mm, an incomplete stripping is predicted. • G IC affects the first stages of stripping evolution, while G IIC does not affect stripping evolution.

•
With increasing the maximum applied pressure, a transition from the annular to the solid stripping pattern takes place. For values smaller than 2500 MPa, an incomplete stripping (annular) is predicted.
From the above conclusions it is obvious that for certain combinations of the process and material parameters an incomplete stripping is predicted. If we try to give an initial rough guideline for the design of the tests, for the process parameters, we can say that the laser spot should be kept between 2.5 and 3.5 mm and the maximum applied pressure above 2550 MPa. Regarding the material parameters, the low thickness of the aluminum substrate and the epoxy paint makes stripping difficult and, in some cases, it leads to incomplete stripping. On the other hand, as expected, the lower the fracture toughness of the aluminum/epoxy interface the easier the stripping. The next step of the work is to design proper tests to characterize the fracture toughness of the aluminum/epoxy interface and after further experimental verification and validation to enable a realistic virtual testing of the stripping process.