*2.2. Methodology*

2.2.1. Anisotropic Hardening Formulation based on Hill's Yield Surface to Model Alternative Loading Paths Presented in LSP

Classic *J*2 plasticity theory assumes that the yield condition in stress space is governed by the second deviatoric stress invariant, *J*2. This implies that the elastoplastic surface in stress space is considered as independent of the hydrostatic pressure and the third deviatoric stress invariant, *J*3. The second invariant, *J*2, is defined as a function of a generic stress state as follows:

$$J\_2 = \frac{1}{6} \left[ (\sigma\_{22} - \sigma\_{33})^2 + (\sigma\_{11} - \sigma\_{33})^2 + (\sigma\_{11} - \sigma\_{22})^2 \right] + (\sigma\_{23})^2 + (\sigma\_{31})^2 + (\sigma\_{12})^2 \tag{1}$$

The most widely implemented yield criterion based on *J*2 theory used to model metallic materials is the von Mises criterion. It assumes that the material yields when the distortion energy per unit of volume reaches a limit. This limit is defined as the distortion energy achieved in a uniaxial tensile test, in which the yield condition is presented when σ1 = <sup>σ</sup>*y* and σ2 = σ3 = 0, where <sup>σ</sup>*y* is the yield stress

of the material and σ1, σ2 and σ3 are the three principal stresses. Thus, the yield surface results in a coaxial cylinder with the hydrostatic axis (<sup>σ</sup>1 = σ2 = <sup>σ</sup>3) in the principal stresses space:

$$\left(\sigma\_{\mathcal{Y}}\right)^{2} = \frac{1}{2}\left[\left(\sigma\_{22} - \sigma\_{33}\right)^{2} + \left(\sigma\_{11} - \sigma\_{33}\right)^{2} + \left(\sigma\_{11} - \sigma\_{22}\right)^{2}\right] \\ + 3\left(\left(\sigma\_{23}\right)^{2} + \left(\sigma\_{31}\right)^{2} + \left(\sigma\_{12}\right)^{2}\right) \\ \tag{2}$$

$$\left(\sigma\_{\mathcal{Y}}\right)^{2} = \frac{1}{2} \left[ \left(\sigma\_{2} - \sigma\_{3}\right)^{2} + \left(\sigma\_{1} - \sigma\_{3}\right)^{2} + \left(\sigma\_{1} - \sigma\_{2}\right)^{2} \right] \tag{3}$$

In 1950, Hill extended the von Mises yield criterion to define anisotropic yield surfaces (Equation (4)), in which the parameters involved *F*, *G*, *H*, *L*, *M* and *N* can be defined as functions of the plastic strain. Thus, the material's elasto-plastic anisotropic behaviour can then be fully modelled. However, it is easier for researchers to define the yield stresses in three different directions, <sup>σ</sup>*y<sup>i</sup>* (*i* = 1 to (3), expressed as a function of a generic yield stress, <sup>σ</sup>*y*, where *Rii* are the first three potentials (Equation (5)). The other three potentials, *Rij* (*i*, *j* = 1 to 3, *i* - *j*, *i* < *j*), affect the material's behavior when shear stresses are presented. These potentials can be correlated with Hill's original parameters [29].

$$\left(\sigma\_y\right)^2 = \frac{1}{2} \left[ F(\sigma\_{22} - \sigma\_{33})^2 + G(\sigma\_{11} - \sigma\_{33})^2 + H(\sigma\_{11} - \sigma\_{22})^2 \right] \\ + 3 \left( L(\sigma\_{23})^2 + M(\sigma\_{31})^2 + N(\sigma\_{12})^2 \right) \tag{4}$$

$$
\sigma\_{y\dot{\imath}} = \mathcal{R}\_{\dot{\imath}\dot{\imath}} \sigma\_y \tag{5}
$$

In LSP, the plasma expansion develops in a high intensity shockwave that propagates through the material, which experiences a phenomenon called uniaxial strain [26], in which the material is subject to alternative loading paths (plastic loading and unloading). In the case of materials of highly anisotropic behaviour, as the considered case of AZ31B alloy, the proper modelling of these alternative loading paths requires an extended yield surface definition, in which the material's compressive behaviour along ND, RD and TD is properly modelled. Since the yield stress of the present alloy is relatively small, low-density treatments are expected to be the suitable ones to enhance its mechanical properties. Thus, the yield surface displacement, usually modelled by means of a kinematic hardening rule, will not be considered.

In the present paper, the plastic loading is modelled by means of the compressive stress–strain curve along the normal direction (ND), which is coincident with the direction of the applied pressure. Therefore, a biaxial compressive state in rolling and transverse directions (RD and TD) is considered during plastic unloading. Intuition suggests the use of the tensile stress–strain curve in the normal direction (ND) during plastic unloading. However, the effect of anisotropy in RD and TD would be neglected as a consequence, which is not realistic since different behaviour is expected along these directions.

### 2.2.2. Model Calibration Based on the Anisotropic Stress-Strain Curves of Mg AZ31B Alloy

The anisotropic stress–strain behaviour of Mg AZ31B alloy has been documented by many authors [30,31]. Several of them have focused their studies on warm temperatures (373 K to 573 K) since magnesium formability is relatively low at room temperature. Others have characterized the stress–strain curves in the ND subjected to uniaxial compressive states [20]. Concerning the cyclic behaviour, several authors studied the hysteresis loops in a wide variety of conditions: temperature, strain amplitude, directions and strain rates [17,30].

A grea<sup>t</sup> percentage of the stress–strain characterization works are centered on the RD and TD. The ND direction is harder to study due to its short length. In consequence, a lack of results is presented especially in the tensile direction. In spite of these difficulties, quasi static ND tensile stress–strain curves have been characterized in cyclic behaviour studies, in which the specimens were obtained from a 50 mm thickness rolled plate [18,30]. However, the dynamic behaviour at high strain rates is not characterized. In this section, an anisotropic model is calibrated based on the experimental dynamic characterization along ND, RD and TD.

The particularities of LSP processes need to be considered in order to select the most suitable experimental results to define the material model. Results need to include high strain rate dependence, while the e ffect of temperature and high compressive stress triaxialities could be ignored. In LSP, the explicit consideration of asymmetry in stress–strain curves may be considered for high-density treatments [26]. Due to the relatively small yield stress of AZ31B alloy, low-density treatments are expected to suitable and, consequently, the explicit consideration of asymmetry is neglected in the present study.

Regarding thermal e ffects at the material's surface, although the generated plasma usually reaches high temperatures, it is maintained during a very short period of time. Thus, the material is not expected to experience large temperature variations and room temperature is used for model calibration. For instance, the thermally a ffected depth for a steel alloy subject to LSP is estimated at 30 μm [32], while the plastically a ffected depth typically reaches more than 1 mm. As a result, thermal e ffects are usually neglected in LSP processes and the hypothesis of considering it as a purely mechanical process is generally adopted [32].

The experimental results that best fit with the above LSP characteristics definition are the ones presented by Tucker et al. [19]. These results include the quasi static and high strain rate behaviour in ND (only compression), RD (both tensile and compressive) and TD (both tensile and compressive) directions.

The plastic loading will be modelled by means of the compressive stress–strain curve in the ND. Therefore, compressive stress–strain curves in RD and TD characterize the biaxial stress state during plastic unloading. The yield stress of the present alloy is relatively small compared with most of the metallic materials subject to LPS treatments. As a consequence, low-density treatments are expected to be the suitable ones. This implies that the e ffect of the cyclic plasticity, modelled by means of a kinematic hardening rule, could be ignored. Thus, the yield surface in the present model will experience an anisotropic expansion in the stress space when the material is plastically deformed.

Regarding the stress–strain curve calibration along ND, a Voce type hardening law [33] has been used, including strain rate dependence (Equation (6)), where σ*ND* is the compressive yield stress in the ND, σ*ND*0 is the first compressive yield stress, *Q* is the yield stress saturated increment, . ε0 is a reference strain rate, *b* and *b*0 are constants to be calibrated, *F* . ε*r* is a strain rate function that models the saturation equivalent plastic strain at di fferent strain rates, and <sup>Δ</sup><sup>ε</sup>*p* is the equivalent plastic strain.

Considering that experimental results are obtained up to 4300 s<sup>−</sup>1, the function *F* . ε*r* is set as constant for further strain rates, since there are no experimental data beyond 4300 s<sup>−</sup>1. It is documented by some authors that extrapolating experimental results for higher strain rates usually leads to overestimations in the compressive residual stress predictions [34]. Table 1 presents the calibrated constants obtained for the best fit option. Figure 2 shows a comparison between the experimental results and the numerical predictions.

ε*r*

=

⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩

.

ε1

2 .

ε0

> ε*r* > .

ε*r*

ε1

$$
\sigma\_{\rm ND}(\Delta \varepsilon\_p) = \sigma\_{\rm ND0} + Q \left( \dot{\varepsilon}\_r \right) \left( 1 - e^{-\left( F(\dot{\varepsilon}\_r)b + b\_0 \right) \left( \Delta \varepsilon\_p \right)} \right) \tag{6}
$$

$$
F(\dot{\varepsilon}\_r) = \begin{cases}
1 & \dot{\varepsilon}\_r < \dot{\varepsilon}\_0 \\
& \dot{\varepsilon}\_1 + 1 \ \dot{\varepsilon}\_0 \le \dot{\varepsilon}\_r \le \dot{\varepsilon}\_1 \\
& & \dot{\varepsilon}\_1
\end{cases} \tag{7}
$$

ε1

$$^{113}$$

**Table 1.** Calibrated constants for the analytic

**Parameter Value** σ*ND*0 (MPa) 178 *Q* (MPa) 125 *b*0 (–) 19 *b* (–) 18 . ε0 (s−1) 0.001 . ε1 (s−1) 4300

predictions

 of

compressive

stress–strain

 curves (ND).

**Figure 2.** Calibrated model vs. experimental stress–strain curves along ND.

The material experiences softening at a certain plastic strain (approximately when the strain is equal to 0.07 at high strain rates), which is not possible to model with a purely Voce expression (monotonically increasing function). Koh used a Swift–Voce equation to predict softening [23]. However, in LSP, the equivalent plastic strain, <sup>Δ</sup><sup>ε</sup>*p*, usually exceeds the material's ultimate tensile stress with no failure. This is due to the high compressive stress triaxialities, η, involved in LSP treatments, ensuring the complete protection of the material. In 2004, Bao and Wierbzicky [35] postulated a cut-off value for the compressive stress triaxiality beyond which failure will never occur (η < −1/3). This is the particular situation presented in LSP treatments. Therefore, a saturated hardening curve is considered to model the stress-strain scenarios beyond the ultimate tensile stress. Thus, Voce's equation itself predicts correctly this particular situation. It is important to point out that the isotropic models, which are widely accepted in LSP simulations (i.e., Johnson–Cook model) predict an unbounded expansion of the yield stress when high-density treatments are applied. This leads to overestimated compressive residual stress predictions, which are not observed experimentally. Thus, the use of a Voce equation (with saturated yield stress) may represent more accurate predictions than classic isotropic hardening models widely implemented in LSP simulations [26].

Regarding the constitutive equations to model the plastic unload, both RD and TD compressive stress–strain curves are considered. These curves present an atypical behaviour at an approximately equivalent plastic strain of 0.059, beyond which a new deformation mechanism is activated and yield stress increases abruptly. In order to model this behaviour, the superposition of two Voce equations will be used: The first one is characterized by a low amplitude. The second one, with higher amplitude, is activated once the equivalent plastic strain reaches the critical value. Equations (8)–(11) represent the stress–strain behaviour in RD and TD, both for quasi static and dynamic conditions. The dynamic behaviour from quasi static conditions up to the highest calibrated strain rates is obtained by linear interpolation. The dynamic behaviour beyond the experimental data is set as constant. The calibrated parameters for the best fit option are shown in Table 2. Figures 3 and 4 show a comparison between the experimental results and the corresponding analytical predictions along RD and TD respectively.

$$\sigma\_{\rm RD}(\Delta \varepsilon\_p) = \min \left( \sigma\_{\rm D0} + Q\_D \left( 1 - e^{-b\_D(\Delta \varepsilon\_p)} \right) + \partial \left( \Delta \varepsilon\_p \right) Q\_{\rm RD2} \left( 1 - e^{-b\_{\rm Dk}(\Delta \varepsilon\_p)} \right), \sigma\_{\rm RD\max} \right) \tag{8}$$

$$\sigma\tau\text{Dc}\{\Delta\varepsilon\_{p}\} = \min\{\sigma\Box 0 + Q\Box\left(1 - e^{-b\_{\Box}\left(\Delta\varepsilon\_{p}\right)}\right) + \bar{\sigma}\left(\Delta\varepsilon\_{p}\right)Q\tau\Box\left(1 - e^{-b\_{\Box\varepsilon}\left(\Delta\varepsilon\_{p}\right)}\right), \text{ }\sigma\tau\text{Dmax}\right\} \tag{9}$$

$$
\sigma\_{\rm RDd} \left( \Delta \varepsilon\_p \right) = \sigma\_{\rm D0} + Q\_D \left( 1 - e^{-b\_D \left( \Delta \varepsilon\_p \right)} \right) + \partial \left( \Delta \varepsilon\_p \right) \cdot Q\_{\rm RDd} \left( 1 - e^{-b\_{\rm RDd} \left( \Delta \varepsilon\_p \right)} \right) \tag{10}
$$

$$
\sigma\_{\rm TDM}(\Delta \varepsilon\_p) = \sigma\_{\rm D0} + Q\_{\rm D} \left( 1 - e^{-b\_{\rm D}(\Delta \varepsilon\_p)} \right) + \partial \left( \Delta \varepsilon\_p \right) \cdot Q\_{\rm TDd} \left( 1 - e^{-b\_{\rm TDd}(\Delta \varepsilon\_p)} \right) \tag{11}
$$

$$
\partial \begin{pmatrix} \Delta \varepsilon\_{\mathcal{P}} \end{pmatrix} = \begin{pmatrix} 0, & \Delta \varepsilon\_{\mathcal{P}} < \begin{pmatrix} \Delta \varepsilon\_{\mathcal{P}} \end{pmatrix}\_{\mathcal{L}} \\\ 1, & \Delta \varepsilon\_{\mathcal{P}} \ge \begin{pmatrix} \Delta \varepsilon\_{\mathcal{P}} \end{pmatrix}\_{\mathcal{L}} \end{pmatrix} \tag{12}
$$

where σ*RDe* and σ*TDe* are the quasi static yield stresses along RD and TD respectively. σ*RDd* and σ*TDd* are the dynamic yield stresses. <sup>Δ</sup><sup>ε</sup>*pc* is the critical equivalent plastic strain beyond which the high amplitude Voce function is activated.

**Table 2.** Calibrated constants for the analytic predictions of stress–strain curves.


**Figure 3.** Calibrated model vs. experimental stress–strain curves along RD.

**Figure 4.** Calibrated model vs. experimental stress–strain curves along TD.

The calibrated stress–strain curves following the presented procedure are plotted together in Figure 5. The reference compressive yield stress, <sup>σ</sup>*y*, is set equal to σ*D*0. Thus, potentials, *Rij*, are obtained dividing the corresponding yield stress in ND, TD and RD by the reference yield stress, σ*D*0. The determination of Hill's coefficients is remarkably conditioned by the process, according to reference [36].

**Figure 5.** Calibrated model vs. adapted experimental stress–strain curves for ND, RD and TD.

2.2.3. Experimental Determination of the Spatial Pressure Pulse Profiles

In LSP, the high intensity laser beam applied to the surface forces a sudden vaporization generating a plasma, whose expansion develops in pressures about 5 GPa on the material's irradiated area. This

pressure generates a shockwave that propagates through the material, leading to plastic deformations and a final residual stress state. Thus, both the characterization of the plasma pressure and the simulation of the shockwave propagation through the material are key issues to model properly LSP treatments.

Regarding the plasma characterization, several authors adopt a simplified expression to estimate the maximum pressure as a function of laser intensity, in which the hypothesis of adiabatic expansion is considered [37]. However, in the present study, the plasma characterization is provided by a detailed calculation of plasma pressure and thermal flux (HELIOS code) as described in previous publications by the authors [38–40].

The laser device used in the present study is a Q-switched Nd:YAG (Spectra-Physics, Santa Clara, USA) operating at 10 Hz and providing 9 ns FWHM, 2 J per pulse with near Gaussian spatial profile. The spot diameter was set to 1.5 mm. This leads to a peak laser intensity, *I* = 28.5 GW/cm2. The material's reflectivity was calculated with the aid of Wu and Shin model [41,42], which is set as an input parameter to the HELIOS code. The in-time pressure profile is represented in Figure 6.

**Figure 6.** In-time pressure evolution calculated by means of HELIOS code.

The deformation profiles after the application of 1 to 5 concentric pulses in rolled Mg AZ31B alloy were obtained in the Laser Center UPM by means of a confocal microscope LEICA DCM 3D (Leica Microsystems GmbH, Wetzlar, Germany). In Figure 7, a characteristic map obtained for the application of a single pulse is presented.

**Figure 7.** Confocal image of the deformed surface after the application of the first pulse.

The general form of the spatial distribution is defined by Equation (13), where *P*(*t*) is the temporal pressure profile, *<sup>a</sup>*φ, *Rc* and *H* are constants to be calibrated, and *r* represents the distance from a generic surface point to the center of the laser spot. In Figure 8, a comparison between the experimental softened profiles and its corresponding numerical predictions for the best fit option is plotted. Table 3 presents these calibrated constants.

$$P(r,t) = P(t)e^{-a\_{\phi} \left(\frac{r}{R\_c}\right)^{2H}} \tag{13}$$

**Figure 8.** Experimental deformation profiles vs. numerical predictions.

**Table 3.** Calibrated constants for spatial pressure definition.

