FE Analysis of Laser Shock Peening on STS304 and the Effect of Static Damping on the Solution

: Laser shock peening creates compressive residual stress on the surface of the material, reducing stress corrosion cracking and increasing fatigue life. FE simulation of laser shock peening is an effective way to determine the mechanical effects on the material. In conventional FE simulations of laser shock peening, explicit analysis is used while pressure loads are applied and switched into implicit analysis to dissipate kinetic energy. In this study, static damping was adopted to dissipate kinetic energy without conversion into implicit analysis. Simulation of a single laser shock and multiple shocks was performed, and deformation and minimum principal stress were compared to evaluate the static damping effect. The history of the internal and kinetic energy were analyzed to compare the stabilization time depending on the damping value. Laser shock peening experiments were also performed on stainless steel 304 material. The residual stress of the specimen was measured by the hole drilling method and it was compared to the FE simulation result. The residual stress from the experiment and the simulation results showed similar distributions in the depth direction. Anisotropic residual stress distribution due to the laser path was observed in both results.


Introduction
The peening process generates a compressive residual stress field on the surface of a material. Shot peening and ultrasonic peening not only largely deform the material surface but also generate the compressive residual stresses within a few micrometers of depth. On the other hand, laser shock peening (LSP) is known to improve the fatigue life by 4 to 5 times compared to other peening processes by forming compressive residual stress at the depth of several millimeters from the material surface. The LSP generates compressive residual stress by a plasma shock induced by the reaction between a high-energy laser pulse and the material [1]. When the plasma pressure is applied on the material surface, it deforms plastically and pushes near the material. As an interaction, the elastic region pushes back to the plastically deformed zone, and the compressive residual stress remains as a result. The compressive residual stress and the grain refinement generated by laser shock peening are known to increase the stress corrosion cracking resistance and fatigue life of the alloy [2,3]. However, the LSP process is mainly applied to high-added-value products because of the high process cost compared to other peening processes [4]. For example, mechanical parts subjected to fatigue loads, such as turbine blades, or subjected to corrosive environments, such as welded parts in nuclear reactors, are the main applications for LSP [5].
Since FEM analysis of the LSP was firstly introduced by Braisted and Brockman [6] in 1999, several studies have been reported based on ABAQUS [7][8][9][10][11][12] and ANSYS LS-DYNA [13][14][15]. The LSP process is usually performed within an instant. Plasma pressure loading takes less than 200 ns and the material plastically deforms. Due to this process characteristic, FE simulation of the LSP is analyzed by explicit methods. The explicit analysis is considered in the case of a non-linear problem, but it requires a much longer calculation time than implicit analysis due to the small time step size. During the pressure loading period, the FE model is solved by an explicit method. To de-excite the model after pressure loading, the explicit analysis result is exported to implicit analysis and it repeats until final loading, as shown in Figure 1 [16]. Peyre et al. [8] performed LSP FE simulation of 12% Cr stainless steel based on the strain-rate-dependent stress-strain relationship of the Johnson and Cook model. Two-step analysis was repeated, where the first step was loading and the second was calculating the deformation and residual stress field. However, the procedure needed 1 ms for stabilization; as a result, it took one hour (wall-clock) per shot. Amarchinta et al. [9] performed LSP FE simulation of Ti-6Al-4V material by using the Johnson-Cook and Zerilli-Armstrong stress-strain models. Depending on the laser pulse energy, residual stress was evaluated according to depth. As a simulation method, the pressure loads were applied through explicit analysis and residual stress was analyzed by implicit analysis. As shown in Figure 1, Hasser et al. [17] proposed SEATD analysis, which does not need to switch into implicit analysis but applies time-dependent damping after the loading period to stabilize the model. This method was performed six-times faster in terms of the total calculation time than the 2N + 1 analysis method that was proposed by Brockman et al. [18]. Bahmare et al. [15] performed LS-DYNAbased LSP simulation with the Zerilli and Armstrong model for Ti-6Al-2Sn-4Zr-2Mo material. Only explicit analysis with damping was used for stabilization, but the damping value was not described in detail. Fameso et al. [19] used explicit analysis only and they employed time-stepped damping during the loading step and viscous damping during the no-loading step. Residual stress derived from the simulation and experiment showed a close correlation in terms of distribution, and total CPU time was reduced to one third compared to conventional explicit-implicit methods.  Compared to the LS-DYNA, the Autodyn requires the optional license to transfer the explicit analysis result to an implicit solver. In this work, to settle the model without converting the analysis from explicit to implicit, static damping is employed. Similar to Fameso et al.'s study [19], the static damping was applied during both the loading and relaxation period. The effect of the static damping on the FE simulation was evaluated from the single LSP and overlapped multiple LSP simulation results. Comparison of the residual stress of the LSP-treated material and LSP simulation was discussed.

Conservation Equations for Explicit Analysis
The partial difference equations (PDE) expressed as mass, momentum and energy conservation need to be solved in explicit analysis. The FE model, boundary condition, initial condition and loading condition are also calculated. Each conservation equation can be expressed as [20] follows. The ρ is density; V is volume; m is mass;ẍ,ÿ andz are the x, y and z direction acceleration; b is the momentum term; σ ij is the stress tensor;ė is the energy rate; andε ij is the strain rate tensor.
Mass conservation equation: Momentum conservation equation: Energy conservation equation:

Static Damping
The ANSYS Autodyn solver was used to perform transient analysis and it presents a guideline for a proper static damping value [21,22]. The procedure introduces a damping force that is proportional to the nodal speed and aims to critically damp the lowest vibration mode of a static system. When the damping option is used, the nodal velocity is calculated every iteration in the solver as Equation (4): where v is the nodal velocity, R d is a static damping value, n is the time of the current iteration of the solver and ∆t is a time step during the solving process.
Here, T is a period of the lowest mode of vibration. The values of ∆t and T can be obtained by after first solving the FE model without the damping option. These values should be estimated by the user and the T should not be underestimated. The values that were used in this study will be discussed later.

Stress-Strain Model
Every material has its own stress-strain curve. These curves are mostly obtained by tensile tests. However, the LSP process deforms the material surface in a short time, and, in the case of a high strain rate (strain per time) process, the stress-strain curves change depending on the strain rate. There are several models of the high strain rate stress-strain relationship that are well known, such as as Johnson-Cook (JC) [23], Zerilli-Armstrong (ZA) [24] and the Khan-Huang-Liang (KHL) [25] equation.
The equation (6) is the JC model, where A, B, C, n and m are the material-dependent coefficients that are determined by the experiment. A is the yield strength; B is the hardening modulus; C is the strain rate dependency coefficient; ε is the effective plastic strain; ε is the effective plastic strain rate;ε 0 is the reference strain rate; n is the work hardening exponential; m is the term of thermal softening [26]. The thermal term of the JC model is neglected because the LSP is regarding a non-thermal process [11,14,27]. The JC model was used in this study for FE simulation.

Pressure Load Model
Several pressure loading models have been proposed in previous studies. The pressure load is expressed as a function of time, a simple triangular load model [13] and a step load model [14]. The pressure load profile applied in this work is shown in Figure 2. This normalized pressure load model was proposed in several studies [9,28]. The peak pressure can be calculated from Equation (7). Here, P peak (GPa) is the peak plasma pressure; α is the correction factor of the internal energy and the thermal energy (usually α ≈ 0.25 ); I (GW · cm −2 ) is the laser pulse power density, where P pulse is the laser pulse power and A peen is the laser spot area at the working point; Z (g · cm −2 s −1 ) is the combined acoustic impedance of the confinement layer Z water (≈0.15 × 10 6 ) of water and the target material Z target (≈4.5 × 10 6 ) of stainless steel [8,12,29,30]. 2 From Table 1 and Equations (7) and (8), the pulse power density of the laser beam I was calculated to be 5.94 GW · cm −2 and the peak pressure P peak was calculated to be 3.5 GPa. Pressure load sharply increased to peak pressure (3.5 GPa) in 3 ns, maintained for 3 ns, dropped to 75% of peak pressure (2.6 GPa) in 3 ns, decreased to 30% of the peak pressure (1.0 GPa) in 81 ns and remained until 80 ns and then dropped to zero. Additionally, the 2D pressure load distribution from the beam center to the radial direction is also considered. Depending on the laser beam pulse profile, Gaussian or flat-top models have been suggested [12]. Sun et al. [31] proposed a mixed shape of Gaussian and flat-top according to measurements from the beam profiler. In this work, a flat-top spatial energy laser beam was used for the LSP experiment and for simulation, and the pressure distribution in radial direction was assumed to be uniform.

Overlapped Multiple LSP
LSP is usually performed by overlapping each single laser pulse. In this work, the laser beam was overlapped by 50% in the horizontal and vertical direction, as shown in Figure 3. Through these repeated shots, the dotted box area appears repeatedly. To reduce the computational time, nine shots of LSP pressure loads were applied and the cross-section of the repeated area was evaluated.

LSP Experiment Setup
In order to verify the FE simulation results, the LSP experiment was performed with the conditions shown in Table 1.
Continuum Powerlite Furie was used as a laser system, the laser pulse energy was measured by Gentec-EO Maestro console with QE65LP-H-MB-QED-D0 energy detector and the beam diameter was checked by a burn paper pattern with low-level pulse energy at the working distance. An actual LSP experimental setup is shown in Figure 4. The beamfocusing head focuses a 20 mm diameter input beam to increase the energy density, and the laser spot size can be adjusted depending on the working distance. The robot made the specimen move and water was supplied near the laser irradiation point to form the 1 mm thickness of the confinement layer. The specimen dimensions were 70 × 50 × 10 mm and the surface was ground condition.

LSP FE Simulation Setup
ANSYS Autodyn was used for LSP FE simulation and for the explicit dynamic analysis. The coordinate system, pressure loading and boundary condition faces are shown in Figure 5. The dimension of the geometry was 20 × 20 × 5 mm, the x-y plane was parallel to the top plane of the geometry, and the z-axis was parallel to the thickness direction. The nine yellow-colored circles indicate the faces where the pressure loads were to be applied according to the LSP process model setting and the 5 faces (bottom and four lateral faces) shown in red were set to the non-reflecting boundary condition. This allowed propagating traveling waves in the outward direction of the geometry so that the energy did not reflect back into the FE geometry [13,14].
The edges of the circles where pressure loads were applied were evenly divided into 72 and the element size was set to 0.15 mm for fine mesh generation, as shown in Figure 5. For the region were pressure loads were not applied, the multi-zone meshing method was employed to generate a coarse mesh with a size of 0.25 mm. The material properties are listed in Table 2 [32]. To apply overlapped pressure loads on the surface, the surfaces were divided into 40 regions. The pressure loads were applied in sequence as shown in Figure 6. For example, in the case of the yellow-colored surface No. 16 in Figure 6, in total, four pressure loads needed to be applied. If the laser pulse period was 10 µs, the loads were applied 1st (t = 0), 2nd (t = 10 µs), 5th (t = 40 µs) and 6th (t = 50 µs). However, each number in Figure 6 is given for convenience and is independent of the order of the loading sequence.
The static damping value was calculated from Equation (5) by solving the FE model without static damping. The time step size ∆t was 1.4 × 10 −10 from the solver output and the T was estimated to be 6 µs from the energy summary in Figure 7. In this study, R d values of 0, 4.6 × 10 −5 and others were chosen to compare how the static damping affected the result of the deformation and minimum principal stress distribution. Each pressure loading period was set to 10 µs.   Figure 7 shows the each energy variation according to time. The total energy of the pressure pulse load was converted into the internal energy, kinetic energy and hourglass energy. Hourglass energy cannot exist in practice, but for calculation stability, it exists when solving the problem of deformation with high velocities [13]. In order to estimate the stabilization time depending on the static damping value R d , the derivative of internal energy with time was plotted. In the case of the static damping R d = 0, it converged to zero with vibration. The stabilization time for the no-damping model was assumed to be 10 µs, and for multiple LSP simulation, the period of each single shot was also set to 10 µs. This is because, at 10 µs on the kinetic energy graph, the remaining kinetic energy is relatively small, compared to the kinetic energy converted from the pressure loading instance. If the remaining kinetic energy affects the multiple LSP simulation, the solution of the no-damping model can diverge or the result may be significantly different from others. In the case of R d = 4.6 × 10 −5 , the variation in internal energy was observed before 6 µs but it later converged to zero. For R d = 50 × 10 −5 , it only took 2 µs to become stabilized without variation in the internal energy. Hasser et al. [17] also reported that the kinetic energy decreases quickly as damping increases. The model without damping showed the largest internal energy. The reason for this is that the applied pressure load was mostly converted into internal and kinetic energy without energy loss due to static damping.

Transient Analysis
Since the LSP is achieved within a few microseconds, the LSP process can be visualized by performing transient analysis. Figure 8 shows the z-directional deformation of the nodamping model at every 0.2 µs. A pressure load was delivered to the geometry within 0.17 µs, as shown in Figure 2. However, deformation due to the pressure load was continued even afterwards. After pressure loading, the maximum amount of deformation appeared at t = 0.5 µs, and the geometry was deformed continuously in the downward direction. After 1.1 µs, plastic deformation was completed and vibration was gradually reduced.

Deformation Distribution
To evaluate the z-directional deformation depending on the R d values, the crosssections of the middle of the geometry were plotted as in Figure 9. In the case of a large damping value, relatively small deformation was observed. In addition, z-directional deformation values were sampled at the top of the surface. As shown in Figure 10, the center of the LSP was slightly less deformed than the surrounding area, and the surface next to the region where the pressure load was applied showed small deformation in the opposite direction.

Minimum Principal Stress Distribution
When overlapped and complicated multiple LSP are simulated, a scalar invariant quantity such as minimum principal stress is the best practice to derive meaningful information from simulation [18]. The minimum principal stress distribution of the cross-sections was captured as shown in Figure 11. To compare the damping effect, the minimum principal stress value was sampled in the depth direction from the center point of the LSP, and it is plotted in Figure 12. Large compressive stress was formed on the surface in the case of large damping values but was relatively shallow. In contrast, deep compressive stress was achieved in the cases in which no or small damping is used. Exceptionally, from Figure 11d, only a small amount of compressive stress remained near the surface. This is due to the large damping, which significantly reduced the energy converted from the pressure load.

Discussion on Single LSP FE Analysis
As shown by the energy summary in Figure 7, when a large static damping value was used, the kinetic energy was quickly dissipated and, as a result, a relatively small amount of deformation was observed. It was also expected that kinetic energy would be concentrated on the surface not being transferred in the depth direction during a short period of time, and relatively large minimum principal stress was intensively formed near the surface. Table 3 and Figure 13 show the static damping effect on the single LSP simulation result. For ease of comparison, all values were normalized based on the result of R d = 0. As expected, when the damping value was increased in the simulation, the model could quickly stabilize but the difference in the deformation and minimum principal stress result became larger compared to the in the case of R d = 0. With R d = 4.6 × 10 −5 , the simulation time could be reduced to 60% within a 10% difference from the no-damping result. This means that the static damping value R d = 4.6 × 10 −5 evaluated from Equation (5) was reasonable. Appropriate static damping can be effective when performing multiple LSP simulation because it can reduce the calculation time for every individual shot.

Deformation Distribution
Z-directional deformation values were sampled in the same manner as single LSP analysis. As shown in Figures 14 and 15, faces No. 10, 16, 25 and 32 were the most deformed because of the four-times overlapped pressure loads. Similar to the single LSP result, as the damping value was increased, the amount of deformation was decreased.

Minimum Principal Stress Distribution
The minimum principal stress distribution of the multiple LSP FE simulation is shown in Figures 16 and 17. In the case of a large value of static damping, even though the pressure loads were overlapped, minimum principal stress was not formed deeply but was formed intensively near the surface of the geometry. In the multiple LSP FE simulation, the effect of the static damping was similar to the single LSP simulation. As the pressure loads were overlapped, the amount of deformation and the minimum principal stress were increased. In other words, the laser pulses should be overlapped to maximize the LSP effect.
From the deformation and the stress analysis results, the no-damping simulation was not diverged and the result was not significantly different from the others. This means that even though there exists a small kinetic energy in the model, it does not significantly effect the solution of the LSP simulation.

Residual Stress Distribution of LSP Material
To verify the simulation results, the residual stress of the actual LSP-treated specimen was measured by the hole drilling method as shown in Figure 18. The main purpose of the LSP experiment was to determine whether the compressive residual stress was formed up to a depth of 1 mm. Accordingly, the residual stress was measured up to 1 mm in the depth direction at the center of the LSP test specimen. The measurement conditions are listed in Table 4. Table 4. Hole drilling test parameters for residual stress measurement.

Material
Analysis Method

Hole Step
Step Method Hole Depth (mm)

Analysis
Step To derive the residual stress, normal stress σ x and σ y were analyzed. Figure 19 shows the residual stress distribution on the solution geometry surface. The residual stress evaluated from the simulation was not evenly distributed. On the other hand, in-plane residual stresses from the hole drilling method within a certain region are assumed to be constant. Keller et al. [33] averaged the nodal stresses in the laser shock peened area from the top surface to the depth direction. Similarly, the residual stress in the dotted area shown in Figure 3 was averaged at every 0.125 mm in the depth direction for comparison with the hole drilling result.

Analysis
The residual stress(σ x and σ y ) of the specimen and the simulation according to depth are plotted in Figure 20. The figure shows that the distribution of compressive residual stress in both the FE simulation and the experiment were similar to each other, in the case of R d = 0 and R d = 4.6 × 10 −5 . The residual stress distribution expressed as σ x and σ y in the depth direction was not significantly different from the minimum principal stress. Comparison of experiment and simulation result The residual stress distribution σ x and σ y showed a small anisotropic relationship, as shown in Figure 20. Kallien et al. [34] reported that, depending on the overlap and sequence of laser shots, a non-equibiaxial residual stress distribution can be observed in AA2024-T3 material. The mechanism of anisotropic stress generation in LSP of ferrite-pearlite steel was reported by Hirano et al. [35]. Depending on the laser pulse irradiation sequence and the pattern, residual stress distribution can be isotropic or anisotropic. Hirano et al. introduced the coverage ratio term, which is expressed as Equation (9).
C is the coverage ratio, N is the total number of laser pulses irradiated, D is the spot diameter and S is the LSP processed area. When C > 2, LSP stress tends to be anisotropic. In other words, the laser pulses are densely irradiated within the unit area, and the difference between σ x and σ y becomes larger. A similar relationship was observed in stainless steel 304 [36]. In this study, the coverage ratio C could be calculated to be 1.77 from the LSP parameters, and it indicates that there was a small non-equibiaxial residual stress distribution.

Conclusions
In this study, LSP FE simulation was performed and both deformation and the residual stress were evaluated. The simulation result was also compared with the LSP experiment specimen to determine whether the proposed simulation process is appropriate. The study can be summarized as follows: • The single and multi-shot explicit dynamic LSP FE analysis was conducted with the proposed static damping model by ANSYS Autodyn. • An implicit analysis that involved stabilizing the geometry can be skipped by adding the static damping value. • The static damping quickly settles the model after shock loading, but with large damping, residual stress is formed near the surface of the geometry, albeit not deeply, and relatively small plastic deformation can be obtained. • A comparison of residual stress measurement results on the LSP specimen with compressed residual stress obtained from the simulation showed a similar tendency. • The residual stress of the FE simulation and that measured by hole drilling showed a small anisotropic relationship due to the LSP process parameters. • By using the proper static damping value in LSP FE simulation, the calculation time can be reduced and there is no significant effect on the simulation result.
Compared to other peening methods, LSP is a high-cost material surfacing process. To evaluate the actual residual stress distribution, additional destructive tests such as hole drilling are also required. Therefore, when quick verification of the residual stress distribution is needed according to the depth formed by the LSP, the proposed FE simulation model can be appropriate.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to an ongoing study.