Integrated Numerical-Experimental Assessment of the E ﬀ ect of the AZ31B Anisotropic Behaviour in Extended-Surface Treatments by Laser Shock Processing

: In recent years, an increasing interest in designing magnesium biomedical implants has been presented due to its biocompatibility, and great e ﬀ ort has been employed in characterizing it experimentally. However, its complex anisotropic behaviour, which is observed in rolled alloys, leads to a lack of reliable numerical simulation results concerning residual stress predictions. In this paper, a new model is proposed to focus on anisotropic material hardening behaviour in Mg base (in particular AZ31B as a representative alloy) materials, in which the particular stress cycle involved in Laser Shock Processing (LSP) treatments is considered. Numerical predictions in high extended coverage areas obtained by means of the implemented model are presented, showing that the realistic material’s complex anisotropic behaviour can be appropriately computed and—much more importantly—it shows a particular non-conventional behaviour regarding extended areas processing strategies.


Introduction
Since 1930, an increasing interest in magnesium alloys has been documented in the aerospace industry [1,2] as a light-weight material. It is attractive for the biomedical industry in implant design, since it is degradable, with outstanding mechanical properties in comparison with polymers [3][4][5]. This has led to a great number of studies regarding mechanical properties enhancement [6][7][8], corrosion resistance [5,[7][8][9][10][11] and microstructure modifications [12][13][14], among others. Single Mg crystal responses to shock loading in various conditions and orientations is also documented, including the study of the spall fracture [15], confirming the main mechanisms of hexagonal-closed-packed (hcp) crystals subject to uniaxial shock compression. In addition, the shock loading response of a single Mg crystal is compared with the one observed in a Mg alloy. The spall fracture results, from reference [15], were complemented in a recent publication with additional experiments at very high strain rates (up to 10 −7 ·s −1 ) [16].
Magnesium alloys are presented in two different forms: The first ones are called casting alloys, which have defects such as inclusions or pores. The second ones are called wrought alloys, whose mechanical properties are enhanced with respect to casting alloys. This group includes rolled and extruded alloys. This paper is centered on the study of a rolled Mg AZ31B alloy, in which the grain is deformed along the rolling direction leading to a relevant asymmetry in its mechanical behaviour [17].
This anisotropic behaviour is motivated by the fact that different deformation modes are active depending on the loading path [18].
In general terms, a great effort has been employed in the experimental characterization of physical properties of magnesium alloys, but there is still a lack of simulation procedures regarding the accurate prediction of the residual stresses field. This may be motivated by the complex anisotropy and the atypical hardening behaviour of rolled Mg AZ31B alloy, in which great differences are observed in stress-strain curves in the normal direction (ND), rolling direction (RD) and transverse direction (TD). This anisotropic hardening has been studied by Tucker in Mg AZ31B alloy for different strain rates up to 4300 s −1 , in which the effect of the specimen temperature in results is also documented [19,20]. Jäger et al. [21] developed quasi static tensile behaviour experiments of hot rolled Mg AZ31 alloys up to 673 K. Biaxial tests were performed at room temperature to study the mechanical formability and elongation to failure [22]. Most of the researchers oriented these studies to understand the high-complexity physics involved in the different slip modes presented, twinning and detwinning in alternative paths, dislocations reorientation and the influence of grain size and texture in the mechanical properties in different loading orientations.
Regarding the developed mechanical constitutive models, several approaches to characterize stress-strain curves have been proposed in different situations [23,24]. Koh used Voce-Swift type isotropic hardening law to model the tensile stress-strain curves in RD and TD directions with temperatures ranging from 373 K to 573 K, with the purpose of evaluating the increased formability in hot rolled specimens [23]. Several authors used a visco-plastic self-consistent model to characterize Mg AZ31 alloys in a large variety of conditions [24]. However, a great number of constants need to be calibrated to fix these models. Thus, the identification of the particular stress cycles, temperatures and strain rates involved in LSP is necessary in order to develop a material model valid for residual stress determination.
In Laser Shock Processing (LSP), a high intensity pulsed laser beam irradiates the material's surface, forcing a sudden vaporization of the irradiated metallic target. The generated plasma expansion develops in a high intensity shockwave that propagates through the material leading to high strain rate deformations. The stresses in a material subject to LSP are characterized by alternative stress states in different orientations in the deviatoric plane. A hydrostatic pressure is simultaneously applied, which increases as the deviatoric stress does [25,26]. The hydrostatic pressure in metallic materials is usually neglected. However, several authors included an equation of state (Mie-Grüneisen) in their models [27,28], which correlates the applied pressure with the material's density. Nevertheless, the material's compressibility in the range of pressures involved in the considered alloy (from 0 to 5 GPa approximately) is properly defined by the Hooke's Law inherent to the elasto-plastic behaviour. Therefore, the hypothesis of neglecting the hydrostatic pressure is assumed and LSP will be considered as a process in which alternative stresses in different directions are applied.
Experimental evidence shows that anisotropic hardening behaviour occurs when alternative cycles are involved: In general, once a material is deformed by a compressive state, further lower amplitude tensile states result in plastic straining (Bauschinger effect). This is usually modelled by a kinematic hardening rule. Although the explicit consideration of the yield surface displacement within a saturated isotropic hardening rule is normally used to model aluminum alloys subject to LSP treatments in order to prevent the prediction of an undefined expansion of the yield surface [26], for the case of the considered alloy, reasonably good agreement is still possible for low-density treatments using a cyclicity-independent (i.e., conventional Johnson-Cook) hardening model.
In this paper, Hill's yield surface is used to properly model material hardening in both direct and reverse deformation paths. The reasonably good agreement between numerical residual stress predictions and experimental results for the highly anisotropic material considered (AZ31B alloy) confirms the appropriateness of the explicit consideration of anisotropic deformation behaviour in generic highly deformable materials.

Material Description
The material used in the present study is a 5 mm thickness plate of Mg AZ31B alloy (3 wt.% Al, 1 wt.% Zn, balance Mg) provided by Magnesium Elektron (Magnesium Elecktron Ltd., Manchester, UK). The Young modulus, E, is equal to 45 GPa. The Poisson ratio is 0.35. The anisotropic stress-strain curves along ND, RD and TD at quasi-static and dynamic conditions were adapted from Tucker's results [19], in which the stress-strain curves are represented for different strain rates. Considering that each of the stress-strain points in the original publication are not provided, a graph digitalization software was used to obtain each of the experimental results (Figure 1), which is necessary in order to develop and calibrate a material model.

Material Description.
The material used in the present study is a 5 mm thickness plate of Mg AZ31B alloy (3 wt.% Al, 1 wt.% Zn, balance Mg) provided by Magnesium Elektron (Magnesium Elecktron Ltd., Manchester, UK). The Young modulus, E, is equal to 45 GPa. The Poisson ratio is 0.35. The anisotropic stressstrain curves along ND, RD and TD at quasi-static and dynamic conditions were adapted from Tucker's results [19], in which the stress-strain curves are represented for different strain rates. Considering that each of the stress-strain points in the original publication are not provided, a graph digitalization software was used to obtain each of the experimental results (Figure 1), which is necessary in order to develop and calibrate a material model.

Anisotropic Hardening Formulation based on Hill's Yield Surface to Model Alternative Loading Paths Presented in LSP.
Classic plasticity theory assumes that the yield condition in stress space is governed by the second deviatoric stress invariant, . This implies that the elastoplastic surface in stress space is considered as independent of the hydrostatic pressure and the third deviatoric stress invariant, . The second invariant, , is defined as a function of a generic stress state as follows: The most widely implemented yield criterion based on 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 = and = = 0, where is the yield stress of the material and , and are the three principal stresses. Thus, the yield surface results in a coaxial cylinder with the hydrostatic axis ( = = ) in the principal stresses space:

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: 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 = σ y and σ 2 = σ 3 = 0, where σ 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 (σ 1 = σ 2 = σ 3 ) in the principal stresses space: 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, σ yi (i = 1 to (3), expressed as a function of a generic yield stress, σ y , where R ii are the first three potentials (Equation (5)). The other three potentials, R ij (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].
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.

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 great 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 effect 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 effects 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 affected depth for a steel alloy subject to LSP is estimated at 30 µm [32], while the plastically affected depth typically reaches more than 1 mm. As a result, thermal effects 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 effect 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, σ ND0 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 different strain rates, and ∆ε p is the equivalent plastic strain.
Considering that experimental results are obtained up to 4300 s −1 , the function F . ε r is set as constant for further strain rates, since there are no experimental data beyond 4300 s −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. Table 1. Calibrated constants for the analytic predictions of compressive stress-strain curves (ND).

Parameter Value
Metals 2020, 10, 195 6 of 18 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, ∆ , 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. 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, ∆ε 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.
where σ RDe and σ TDe are the quasi static yield stresses along RD and TD respectively. σ RDd and σ TDd are the dynamic yield stresses. ∆ε p c is the critical equivalent plastic strain beyond which the high amplitude Voce function is activated.

Parameter
Value Metals 2020, 10, 195 7 of 18 Where and are the quasi static yield stresses along RD and TD respectively. and are the dynamic yield stresses. is the critical equivalent plastic strain beyond which the high amplitude Voce function is activated.     The calibrated stress-strain curves following the presented procedure are plotted together in Figure 5. The reference compressive yield stress, σ y , is set equal to σ D0 . Thus, potentials, R ij , are obtained dividing the corresponding yield stress in ND, TD and RD by the reference yield stress, σ D0 . The determination of Hill's coefficients is remarkably conditioned by the process, according to reference [36].  The calibrated stress-strain curves following the presented procedure are plotted together in Figure 5. The reference compressive yield stress, , is set equal to . Thus, potentials, , are obtained dividing the corresponding yield stress in ND, TD and RD by the reference yield stress, . The determination of Hill's coefficients is remarkably conditioned by the process, according to reference [36]. 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

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][39][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/cm 2 . 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.
Metals 2020, 10,195 9 of 18 by a detailed calculation of plasma pressure and thermal flux (HELIOS code) as described in previous publications by the authors [38][39][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, = 28.5 GW/cm 2 . 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. 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. 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. 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.  The general form of the spatial distribution is defined by Equation (13), where P(t) is the temporal pressure profile, a φ , R c 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.
Metals 2020, 10,195 10 of 18 The general form of the spatial distribution is defined by Equation (13), where ( ) is the temporal pressure profile, ɸ , and are constants to be calibrated, and 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.

Experimental Characterization of the Residual Stresses for Different Input Parameters.
The residual stresses have been obtained experimentally by means of the hole drilling method (ASTM E837-13 standard) [43], in which the material is drilled and the resulting deformation as the material is removed is measured by a precision strain gage up to 1 mm thickness (HBM K-RY61-1.5/120R-3 prewired). The specimens subject to LSP treatment were obtained from a 5 mm thickness rolled Mg AZ31B plate. Four different treatment strategies were implemented experimentally, covering a treated squared area of 900 mm 2 (Figure 9 shows a schematic representation): 1) Setting a laser spot diameter of 1.5 mm, with a distance between successive pulses of 0.66 mm, developing in an Equivalent Overlapping Density, EOD [44], of 225 pp/cm 2 , and an Equivalent Local Overlapping Factor, ELOF [44], of 4. The peening direction (PD) is set coincident with the transverse direction, TD.

Experimental Characterization of the Residual Stresses for Different Input Parameters
The residual stresses have been obtained experimentally by means of the hole drilling method (ASTM E837-13 standard) [43], in which the material is drilled and the resulting deformation as the material is removed is measured by a precision strain gage up to 1 mm thickness (HBM K-RY61-1.5/120R-3 prewired). The specimens subject to LSP treatment were obtained from a 5 mm thickness rolled Mg AZ31B plate. Four different treatment strategies were implemented experimentally, covering a treated squared area of 900 mm 2 (Figure 9 shows a schematic representation): (1) Setting a laser spot diameter of 1.5 mm, with a distance between successive pulses of 0.66 mm, developing in an Equivalent Overlapping Density, EOD [44], of 225 pp/cm 2 , and an Equivalent Local Overlapping Factor, ELOF [44], of 4. The peening direction (PD) is set coincident with the transverse direction, TD.
(2) Setting a laser spot diameter of 1.5 mm, with a distance between successive pulses of 0.66 mm, developing in an Equivalent Overlapping Density, EOD, of 225 pp/cm 2 , and an Equivalent Local Overlapping Factor, ELOF, of 4. The peening direction (PD) is set coincident with the rolling direction, RD.
(3) Setting a laser spot diameter of 1.5 mm, with a distance between successive pulses of 0.50 mm, developing in an Equivalent Overlapping Density, EOD, of 400 pp/cm 2 , and an Equivalent Local Overlapping Factor, ELOF, of 7. The peening direction (PD) is set coincident with the transverse direction, TD.
(4) Setting a laser spot diameter of 1.5 mm, with a distance between successive pulses of 0.50 mm, developing in an Equivalent Overlapping Density, EOD of 400 pp/cm 2 , and an Equivalent Local Overlapping Factor, ELOF, of 7. The peening direction (PD) is set coincident with the rolling direction, RD.  The EOD corresponds to the number of applied pulses per unit of area, while the ELOF represents the average number of pulses at a given location. All the studied configurations can be considered as low-density treatments according to their respective EOD's and ELOF's. This is consistent with the proposed material model, in which the cyclic plasticity modelling is neglected. Higher density treatment-modelling would require further research in the proper calibration of stress-strain asymmetry. The EOD corresponds to the number of applied pulses per unit of area, while the ELOF represents the average number of pulses at a given location. All the studied configurations can be considered as low-density treatments according to their respective EOD's and ELOF's. This is consistent with the proposed material model, in which the cyclic plasticity modelling is neglected. Higher density treatment-modelling would require further research in the proper calibration of stress-strain asymmetry.
The experimental results from the material's surface up to 1 mm for configurations (1), (2), (3) and (4) are represented in Figures 10 and 11. In these figures, the in-depth minimum, S min , and maximum in plane, S max are represented with their corresponding uncertainties, calculated with the aid of Hdrill software. Several conclusions can be extracted from the experimental results: (1) For all the experimental measurements developed, the minimum stress, S min , is almost coincident with the stress along the PD, while the maximum in-plane stress, S max , is similar to the stress along a perpendicular direction to the ND and PD This result is independent of the relative orientation between the PD and the RD.
(2) While significant differences are observed between treatment (1) and the others, great similarities are presented between configurations (2), (3) and (4). This suggests that the residual stresses tend to reach a saturation value for relatively low-density treatments, which is as expected considering the relatively low yield stress of the present alloy in comparison with most of metallic materials.
(3) Setting the peening direction (PD) coincident with the rolling direction (RD) leads to greater peak compressive residual stresses for low-density treatments ( Figure 10).
Metals 2020, 10,195 12 of 18 considering the relatively low yield stress of the present alloy in comparison with most of metallic materials.
3) Setting the peening direction (PD) coincident with the rolling direction (RD) leads to greater peak compressive residual stresses for low-density treatments ( Figure 10).

Realistic Modelling Results for Extended Surface High-Coverage LSP Treatments.
The advances in computational resources in the last decade have motivated an increase in the number of publications regarding numerical predictions of extended LSP treatments in steel, aluminum and titanium alloys, with different loading strategies. However, there is still a lack of

Realistic Modelling Results for Extended Surface High-Coverage LSP Treatments.
The advances in computational resources in the last decade have motivated an increase in the number of publications regarding numerical predictions of extended LSP treatments in steel, aluminum and titanium alloys, with different loading strategies. However, there is still a lack of

Realistic Modelling Results for Extended Surface High-Coverage LSP Treatments
The advances in computational resources in the last decade have motivated an increase in the number of publications regarding numerical predictions of extended LSP treatments in steel, aluminum and titanium alloys, with different loading strategies. However, there is still a lack of results concerning magnesium alloys which may be caused by the significant anisotropy presented. In this section, the corresponding numerical predictions of four particular treatments are obtained by means of two different models: (1) Considering the anisotropic model presented in Section 2.2, in which the stress-strain behaviour along ND, RD and TD directions is fully characterized. (2) Considering a purely isotropic model, in which both the plastic loading and unloading are modelled by the compressive stress-strain curve in ND direction. The purpose of comparing these two models is to show the impact in numerical predictions of the explicit consideration of anisotropy, since conventional LSP models consider a unique stress-strain curve. Results confirm that representative differences are predicted between the anisotropic model and the isotropic one. Better agreement with experimental results is obtained by means of the anisotropic model.
In realistic high-coverage extended treatments, such as the ones presented in this study, the peening direction motivates anisotropy in the residual stress predictions even in isotropic materials. Concretely, higher compressive residual stresses are usually predicted numerically and obtained experimentally along the PD [45]. This effect is also presented in the considered anisotropic alloy, in which results show that the minimum principal stress direction is coincident with the PD. The relative orientation between the PD and RD motivates variations in S max and S min which cannot be predicted by purely isotropic hardening models.
The numerical predictions have been extracted and averaged from a representative squared area of 2 mm 2 up to 1 mm depth, which represents a similar volume than the one removed by means of the hole drilling method. Figures 12-15 show a comparison between analytical predictions and experimental results for strategies (1), (2), (3) and (4) respectively. Table 4 shows a comparison between experimental results and the analytical predictions obtained by means of the anisotropic model. (d max ) i is the depth at which the maximum peak compressive residual stress is presented. min(S min ) i represents the peak compressive residual stress. i = isot and i = anisot corresponds to the analytical predictions of isotropic and anisotropic models respectively, and i = exp to the experimental results.
Metals 2020, 10, 195 13 of 18 In realistic high-coverage extended treatments, such as the ones presented in this study, the peening direction motivates anisotropy in the residual stress predictions even in isotropic materials. Concretely, higher compressive residual stresses are usually predicted numerically and obtained experimentally along the PD [45]. This effect is also presented in the considered anisotropic alloy, in which results show that the minimum principal stress direction is coincident with the PD. The relative orientation between the PD and RD motivates variations in and which cannot be predicted by purely isotropic hardening models.
The numerical predictions have been extracted and averaged from a representative squared area of 2 mm 2 up to 1 mm depth, which represents a similar volume than the one removed by means of the hole drilling method. Figures 12-15 show a comparison between analytical predictions and experimental results for strategies 1), 2), 3) and 4) respectively. Table 4 shows a comparison between experimental results and the analytical predictions obtained by means of the anisotropic model. ( ) is the depth at which the maximum peak compressive residual stress is presented. ( ) represents the peak compressive residual stress. = and = corresponds to the analytical predictions of isotropic and anisotropic models respectively, and = to the experimental results.        Regarding the principal stress directions, the PD sets the minimum stress curve, S min . The maximum in plane stress, S max , is perpendicular to the ND and the PD. This is predicted correctly by both models. However, the effect of the relative orientation between the PD and the RD in the stress curves, which is significant between the low-density treatments (strategies (1) and (2)), can only be simulated by the anisotropic model.
Concerning the comparison between experimental results and the analytical predictions, the anisotropic model provides a reasonably accurate calculation the effect of the relative orientation between the PD and RD for strategies (1) and (2). A comparison between Figures 12 and 13 shows that setting the PD coincident to the RD leads to greater peak compressive residual stresses. In addition, the depth at which this peak value is achieved also increases.
The anisotropic model provides, in general, more accurate predictions than the isotropic one. In particular, the isotropic model predicts high compressive residual stresses at material's surface, which is not consistent with experimental results. In treatments where the EOD is set equal to 400 pp/cm 2 (strategies (3) and (4)), experimental evidence shows that the residual stresses tend to saturate. Reasonably good agreement is obtained by means of the anisotropic model for strategy (3) while an overestimation of the compressive residual stresses is predicted for configuration (4). This is not surprising since the model is conceived for low-density treatments, which are the most suitable ones for the present alloy.

Discussion
This paper provides a complete methodology to obtain residual stress predictions in anisotropic materials subject to low-density LSP treatments. In summary: Firstly, the material's stress-strain curves in ND, RD and TD directions were properly modelled both at low and high strain rates. Then, the spatial-temporal pressure pulse distribution was calibrated with the aid of experimentally determined deformation profiles. Finally, the average in-depth residual stresses were obtained for four different treatment strategies, showing the effect of the treatment density and the relative orientation between the PD and RD in results. A reasonably good agreement is presented between the experimental results and the analytical predictions obtained by means of the proposed anisotropic model.
Experimental results suggest that the PD sets the minimum principal stress direction, which implies that the anisotropy in residual stress distribution is mainly motivated by the treatment strategy itself. In fact, the minimum stress, S min , is obtained along the PD for the four treatment strategies experimented. This is properly simulated by both the isotropic and the anisotropic models. Nevertheless, the relative orientation between the PD and the RD motivates differences in the minimum stress curve, S min , and the maximum stress, S max , in low-density treatments. This effect cannot be predicted by the isotropic model, since a unique stress-strain curve is defined.
Regarding the comparison between the numerical predictions and the experimental results, better agreement is obtained by means of the implemented anisotropic model. This is observed specially at material's surface, where the isotropic model predicts high compressive residual stresses that are not presented in experimental results. The residual stress distribution from 200 µm to 700 µm is predicted with great accuracy for configurations (1), (2) and (3). The greatest differences between the analytical predictions of the anisotropic model and the experimental results are presented near to the material's surface, where numerical results tend to predict higher tensile residual stresses than the ones observed experimentally. This may be caused by the evaporation of the material's surface during the application of successive pulses, leading to removal of already deformed material. This effect is not considered in LSP simulations [25][26][27]32,34,45].
The presented procedure is particularized for Mg AZ31B alloy. However, it is applicable to any material of interest subject to relatively low-density LSP treatments. In addition, the effect of the explicit consideration of mechanical anisotropy may differ depending on the selected material, since it is motivated by variations in the stress-strain curves along different loading paths.

Conclusions
This paper provides a complete guideline to model anisotropic behaviour of metallic materials subject to LSP treatments. Experimental evidence, confirmed by the corresponding numerical predictions, leads to the identification of the material's stress-strain anisotropy as a key factor to appropriately compute residual stress predictions in low-density treatments. The main conclusions are as follows: (1) The relevance of including anisotropic behaviour in material modelling in residual stress predictions has been demonstrated in the particular case of Mg AZ31B alloy. Conventional isotropic models are not able to predict accurately the in-depth residual stress distribution. In addition, incorrect principal stress directions are obtained by means of purely isotropic models.
(2) The experimental results show that the anisotropy in residual stresses after treatment are motivated mainly by the treatment strategy. The minimum stress, S min , is similar to the stress along the PD, while the maximum in plane residual stress, S max , is similar than the one presented along a perpendicular direction to the ND and the PD. Both the implemented isotropic and anisotropic models predict properly this effect. Nevertheless, better agreement between numerical predictions and experimental results is achieved by the proposed anisotropic model.
(3) S min and S max curves are influenced by the relative orientation between the PD and the RD. Setting the PD coincident to the RD leads to greater compressive residual stresses for low-density treatments (EOD = 225 pp/cm 2 ). This effect cannot be computed properly by isotropic models. Therefore, the anisotropic hardening model is required for this purpose.
Regarding further developments, several issues may be considered: (1) Modelling the material loss during target irradiation. The application of successive pulses leads to the evaporation of already deformed surfaces, which may modify residual stresses especially at material's surface.
(2) Analyzing the effect of the anisotropic modelling in residual stress predictions with a large variety of different input parameters, including additional relative orientations between the PD and material's anisotropic directions.