Finite element modeling of melt pool dynamics in laser powder bed fusion of 316L stainless steel

In laser powder bed fusion (LPBF), the stability of melt pool dynamics determines the overall quality of a manufactured component. In this work, a numerical model of the LPBF process was developed in order to study and fully understand the behavior of the melt pool dynamics. The numerical model takes into account most of the manufacturing parameters, thermophysical properties, an enhanced thermal conductivity approach, and a volumetric heat source in order to precisely simulate LPBF. This research assumes that the energy emitted by the laser interacts with the metal powder with an absorptivity gradient through the layer thickness in order to calculate the thermal history of the process and the evolution of the melt pool dimensions. The obtained results determined that melt pool dimensions follow a thermal pattern, which is caused by the laser scanning strategy of the LPBF process. A new effective width criterion was proposed in the present research in order to accurately relate both calculated and measured dimensions of the melt pool, reducing the relative error of the model and obtaining data scattering with a standard deviation of ± 7.21 µm and a relative error of 2.92%.


Introduction
Within the additive manufacturing (AM) processes known as powder bed fusion (PBF) technologies, one of the most studied laser based techniques is commonly known as selective laser melting (SLM). SLM supplies a thin layer (20 µm-1 mm) of metal powder on a plate or substrate and afterwards, a high-density power laser melts selective regions of metal powder according to data from a CAD model. Once the laser scanning has been completed, a new powder bed is deposited on the surface of the consolidated metal and the process repeats until the part is completely manufactured [1].
Understanding the different consolidation mechanisms that occur during the manufacturing process is of great relevance. For instance, in metals, both melting and solidification phenomena are fundamental mechanisms behind the SLM process, which consolidates metal particles in order to manufacture a fully dense and functional part. Most of the common defects associated with SLM are the development of porosity, non-consolidated particles, high level of residual stress, and not connected layers [2]. Due to the rapid melting and solidification of the powder bed generated by a mobile heat source (laser), different regions of the workpiece experience repeated cycles of melting/solidification which directly affect its microstructure and final properties [3][4][5][6]. Besides, as a result of the nature of the additive process, experimental measurements of the temperature in the workpiece are only possible on easy-access surfaces, not in internal regions [5]. Transient thermal fields and histories are previous requirements to establish an understanding of thermal and physical phenomena that alter the metallurgical quality of the final product, such as cooling velocity, solidification parameters, microstructure, residual stress, and geometric changes (mainly shrinkage) [7]. For such reasons, numerical simulations are an alternative to determine the temperature evolution during SLM, as well as to achieve a more thorough understanding of the process parameters and defects regarding thermal gradients and parts quality.
Some of the challenges that still need to be overcome when modeling and simulating SLM are (1) accurate modeling of the heat source and thermal radiation absorptance of the material, (2) consideration of the heat loss during the process, (3) anisotropy of the process, and (4) the solidification behavior and phase transformation of the material [8,9]. To date, numerous models have been proposed to simulate the SLM process, from techniques such as finite element (FE) to computational fluid dynamics (CFD). One of the main factors in determining the type of simulation technique is the domain size. For example, some authors used CFD method and metal powder discretization which it makes a more robust simulation [10][11][12][13][14], but these studies were limited to domains so small that they were only a few hundred micrometers in single-track simulations. Besides, adding a mass transport term to energy equations and powder discretization may result in complexity and a high computational load. Additionally, the previous studies do not show the evolution of the molten pool geometry while the efficiency a finite element study enables the use of larger domains [2] and multi-track multi-layer simulations of SLM [15].
Some authors have identified and studied the behavior of the melt pool geometry. Weirather et al. [16] carried out a numerical simulation model of SLM based on the smoothed-particle hydrodynamics (SPH) meshless computational method and observed a transient behavior in the geometry of melt pool through the laser track. Bruna-Rosso et al. [2] detected a reduced melt pool size and a lack of fusion of metal powder in the first tracks of both FE simulation and experimental validation. The thermal simulations results from Song et al. [17] demonstrated that the scanning strategies take a great influence on temperature fields and hardly affect the melt pool size. A multi-layer single-track FE model was proposed by Liu et al. [7] to study the subsequent thermal cycling SLM, and the results showed that as the number of layers increases, the maximum temperature, dimensions, and liquid lifetime of the molten pool increase as well. A three-dimensional FE model was developed by Li et al. [18] and found that the width, depth, and length of the molten pool region are inversely proportional to the laser power and scan speed.
Despite efforts mentioned in previous studies, the characterization of melt pool dynamics was limited to observations and comments, and there is no approach to mathematically characterize the melt pool dynamics and geometries during SLM. Therefore, a comprehensive study of the melt pool geometries during SLM has not been studied yet, neither the description of direct thermal phenomena that affects melt pool dimensions, which are an important indicator to detect critical thermal zones. The aim of this work is to evaluate the temperature behavior over time and the dimensional temperature distribution during SLM to describe the dynamics of the melt pool dimensions under the influence of the scanning strategy. This study presents a non-linear and transient three-dimensional FE model to simulate the SLM process of SS316L and to predict the thermal response of the system. The multi-layer model considered most of the processing parameters involved in SLM, such as laser power and velocity, powder bed thickness, and inherent properties to the material which include absorptance, particle size, thermal conductivity, heat capacity, and enthalpy (in liquid, solid and powder state).
The laser model is considered as a volumetric heat source, which has a dynamic radiation absorptance depending on the thickness of the powder bed. The model describes a rectangular strategy to consolidate the metal powder layer-bylayer with the purpose of studying the behavior of molten pool dimensions from the point of view of dynamic systems. Finally, an anisotropic enhanced thermal conductivity methodology was introduced for temperatures above the melting point due to convection produced inside the molten pool (Marangoni effect).

Thermal modeling
In order to model SLM, the energy conservation and heat conduction equations are considered, resulting in the thermal equilibrium equation for tridimensional heat transfer in a material with anisotropic thermal properties as described in Eq. (1) [9].
where ρ is the material density (kg/m 3 ), C p is the isotropic heat capacity (J/kg-K), T is temperature (K), t is time, Q is the heat generated by the heat transfer medium in a unit volume (W/m 3 ), ∇ is a divergence operator and r is relative reference or spatial coordinate in the system volume domain, and k is the material's thermal conductivity (W/m K) in function of temperature.
The inert atmosphere is usually kept at 20 °C, while the metal substrate is preheating at 200 °C as during the actual SLM process. Previous values are considered as initial condition in the simulation according to Eq. (2).
(2) T(x, y, z)| t=0 = 20 • C at gas domain 200 • C at metal domain In SLM, the metal powder is consolidated by laser thermal energy and forms the molten pool. Therefore, the three main ways of heat transfer (conduction, convection, and radiation) are described according to Eq. (3).
where ϵ is the thermal emissivity radiation, h conv is the thermal convection coefficient, and σ St is Stefan-Boltzmann constant.
Due to material phase changes and thermophysical properties in SLM, a non-linear implementation is necessary [5,9,19]. Moreover, a special approach is presented for those material properties such as density, thermal conductivity, and enthalpy, which depend on whether the state of SS316L is as powder, solid or liquid. The thermophysical properties of the present alloy are obtained from Mills [20]. However, the powder bed is a mixture composed of metal particles and voids among them occupied by inert gas. Therefore, powder bed density pwd can be described as a rule of mixtures between metal particles s and argon g according to Eq. (4) [21].
where ε is the powder bed porosity. The thermophysical properties of SS316L are depicted in Fig. 1. The density behavior of both powder and solid state can be observed in Fig. 1a, where a spike is observed in the transition region. Figure 1b shows the thermal conductivity for metal powder and solid state, and Fig. 1c presents the heat capacity and enthalpy of 316L. The thermal conductivity of the metal powder was determined according to the equations presented by Sih and Barlow [22,23].
In order to compute the powder thermal conductivity k p , the morphology of particles is assumed to be spherical and behave according to Eq. (5), where k g is defined as the thermal conductivity of inert gas surrounding the particles (4) pwd = g + (1 − ) s Fig. 1 Thermophysical properties of 316L steel alloy. a Metal powder and solid densities. b Metal powder and solid thermal conductivities. c Heat capacity and enthalpy taken from Lemmon and Jacobsen [24], ε is the powder bed porosity, and k s is the thermal conductivity of solid 316L steel. The emission of the total radiation amongst particles inside of the powder bed is calculated by Eq. (6), which considers the fraction of pores on the surface A H (Eq. 7), ϵ H is the emissivity of vacancies over the powder surface (Eq. 8), and k r defines the heat transfer accounted for radiation among particles of diameter D p according to Eq. (9), where σ defines Stefan-Boltzmann constant.
Simulations of heat transfer by convection require calculations of mass velocity fields, which are difficult and computationally demanding. However, there are evidences that show this accelerated heat transfer is not taken into account in the molten pool which may result in inaccurate temperature records and cooling velocity [5,25].
It is possible to consider two important effects due to thermal dependence of thermophysical properties developed inside the melt pool. The first effect is the thermal capillarity convection phenomenon or so-called Marangoni [26,27], and the second effect is the buoyancy by density change. The convection by thermo-capillarity occurs when a superficial tension gradient is developed inside the flow, as the surface tension of the liquid metal is higher compared to other liquids (due to a combination of dispersion force and metallic bonds), the metal tends to flow [28].
Due to the nature of the FE model, the mass transfer was neglected. However, some researchers [25,29,30] have used an enhanced thermal conductivity approach to include the heat transfer by convection in the molten pool, which can manipulate the thermal distribution, and the (5) effects may be similar to those induced by molten pool convection [25]. Hence, the energy equation was determined by the introduction of an enhanced thermal conductivity factor k ′′ ii by the Eqs. (10) and (11).
where k is defined as the isotropic thermal conductivity related to the temperature, the subscript ii is the spatial coordinate, and F ii is the anisotropic enhanced factor for a given spatial coordinate described by Eq. (10). Based on previous work [31], it was concluded that the most suitable enhanced factor has the setup (F xx , F yy , F zz ) = (1, 5, 10). Moreover, the molten pool surface and top layer are assumed to be flat and the simulation does not consider heat and mass losses due to vaporization.
During SLM, temperature changes due to the heat flux are attributed to sensitive heat, while the release or absorption of thermal energy at constant temperature is due to phase changes or latent heat of the material. Both thermal phenomena are introduced by enthalpy H [32], in volume terms, defined by Eq. (12).
Numerous approaches have been proposed to model the laser as a mobile heat source [21,26,29,[33][34][35][36] for SLM SS316L. However, considering the laser as a surface heat source is acceptable if the power value is less than 100 W [32,37] and this was confirmed on previous work [31]. Otherwise, modeling the laser with a volumetric heat source is suitable for high power and speed values [15,21,[38][39][40][41]. The analysis achieved by Gusarov et al. [42,43] resulted in a deduction of a dynamic absorptivity model, which consists in considering a volumetric heat source along with two factors. The first factor is the absorptivity, which takes place when the energy is absorbed from the powder bed surface. The second factor is the reflected radiation, which enters interparticle voids (pores) that are placed beneath the powder bed surface.
The absorptivity A is specified as the net density rate of absorbed radiation energy flux by the powder bed to the laser power density, and its behavior is based on the powder bed thickness [42]. Therefore, the laser radiation travels through the powder bed pores to an effective depth of several particle diameters due to multiple diffuse reflection. Consequently, laser energy is not only supplied on the powder bed surface but through the powder bed thickness (z direction) according to Eq. (13) where β is defined as an extinction coefficient, λ is defined as the effective penetration thickness of the laser energy (which is equivalent to the thickness of the layers), and ρ r describes the metal powder's hemispherical reflectivity (when in dense state). Willy et al. [44] performed optical measurements for different materials which used a SLM 280 HL machine equipped with a 1.07-μm fiber laser wavelength; the reflectance value for SS316L for this wavelength is 0.682. Introducing ρ r determines the presumption that the reflection is diffuse and non-specular. Moreover, the porosity of the SS316L metal powder was measured, which is shown in Table 1. The variables D, a, λ, and β, are determined using Eqs. (14) to (17). Table 1 provides the parameters that were used to determine the powder bed dynamic absorptivity, in which the thickness is characteristic data of the process. The particle diameter D p is assumed according to the value of Table 1 due to the spreading characteristics of metallic powder (30 μm layer thickness) during SLM; the powder roller (coater) probably selects particles smaller than 30 μm.
The behavior of the dynamic absorptivity distributed in the depth of the powder bed is represented in Fig. 2a. Gusarov et al. [42] proposed that the absorbed energy density by the metal powder bed caused by the laser irradiance is assumed as a volumetric heat source based on dynamic absorptivity (Eq. 18), where the dA/dz component is the absorptivity exchange rate of the powder bed through the thickness, see Fig. 2b.
The Q 0 component is the surface power density with Gaussian distribution which is expressed in Eq. (19), in which the laser power is represented as P, the normal distribution factor as f (f = 1), the scan speed as v, and the laser radius as r b . The x-direction is considered parallel to the movement of the heat source, whereas the y-direction is considered as normal to the laser path [45].  Fig. 2 a Absorptivity A through the powder bed thickness (z) and b absorptivity exchange rate through the powder bed thickness (z) [31] 3 Numerical simulation A three-dimensional thermal model was developed to simulate a rectangular multi-layer scanning strategy in a SLM process using a FE simulation. The geometry used for this study is shown in Fig. 3. The schematic indicates the three main domains. The first one represents the base plate, also known as substrate, whose main function is to support the powder bed and the consolidated part in the SLM process. Another domain represents the metal powder, and the last one is the solidified metal by the laser action. The thermophysical properties for these three domains are described in Fig. 1 where the substrate and consolidated metal is considered as solid state. The selected dimensions prevent undesirable effects due to the boundary conditions of non-heat transfer on the lateral faces and guarantee a suitable thermal interaction between the solid parts and powder. The FE model was implemented in ANSYS MECHANICAL APDL ® with a high refinement mesh at the interaction area between laser/powder (central region) in which the element size was 10 µm 3 and the coarser mesh was up 196 × 94.3 × 33.3 µm 3 . Both mesh regions can be observed in Fig. 4, and the entire geometric model was defined with a SOLID70 element type characterized by its 3D thermal conduction capacity. SOLID70 elements consist of 8 nodes and one degree of freedom (temperature). Because of the additive nature of the process, the Element Birth and Death function is used to activate layers during the SLM simulation.
During the SLM process, the substrate was preheated to about 200 °C. Therefore, the simulation substrate domain has a fixed temperature of 200 °C (473 K) during the process. The bottom and lateral faces were considered adiabatic. The top surface of powder/metal is in contact with inert gas (argon) where radiation and convection cooling occur, so the boundary condition is implemented according to Eq. (3) with a gas temperature in the chamber of 20 °C and the convective heat transfer coefficient h = 80 Wm −2 K −1 was taken from [46].
The Jacobi conjugate gradient (JCG) iterative solver was selected to determine the transient thermal phenomena which is efficient for heat transfer field problems [47]. The thermal elements have a single thermal degree of freedom per node and the matrices are very well conditioned and its initial JCG implementation uses a sparse matrix file formed from the frontal assembly process [48].
The laser scanning strategy is zigzag according to Fig. 5; six paths per layer and the laser/material interaction zone is 500 × 620 μm 2 . Due to the discrete nature of the model, the laser movement has been discretized; it advances each 20 µm in the indicated direction with an irradiance time of 86.96 µs (calculated from the scanning speed). Figure 5 contains a series of red circles (laser spot), and each one is spaced 20 µm away (point distance). Each path consists of 22 spots, 132 spots per layer, and 660 spots in total (5 layers). The process parameters were established according to recommendations by the manufacturer. A customized subroutine code was implemented in ANSYS Parametric Design Language (APDL) to perform the movement of a volumetric heat source as rectangular scanning strategy. Furthermore, a secondary subroutine was developed to determine whether the powder temperature overcomes the melting point, and if so, turns it into liquid metal and subsequently solidifies when the laser spot moves away. Table 2 presents the processing parameters used in FE simulation, and the general flowchart of the simulation procedure is illustrated in Fig. 6, which explains the main steps involved in thermal analysis.

Experiments
A SLM Solutions 280 machine was used in this study, which is equipped with an IPG fiber laser system of 400 W nominal power, a wavelength of 1065 nm, and a spot diameter of 70 µm. The selected metallic powder was SS316L, and test samples were manufactured with 10 × 10 × 10 mm as dimensions. The test samples were manufactured in an inert atmosphere using argon on the AM system. The specimens were traditionally cut in the perpendicular plane to the thickness. Subsequently, the samples were prepared with standard metallographic procedures and an electrochemical attack was performed to reveal the microstructure. The etchant used was composed by a mix with a proportion of 10 g of oxalic acid and a hundred milliliters of distilled water; then, a voltage of 6 V (direct current) was applied for 60 s to expose the microstructure. A Zeiss optical microscope AX10 was used to examine the microstructure. The SLM processing parameters are enlisted in Table 2. Experimental measurements of the obtained melt pools of the test samples were performed to estimate the temperature profile. The obtained microstructural characteristics of the test samples can be observed in Fig. 7, where melt pool measurements are illustrated. The obtained morphology of the melt pool determines that its width and depth can be directly related with those regions that were exposed to temperatures above the melting point. The measured width (y-axis) and depth (z-axis) can be further examined in contrast with the simulation results. The obtained melt pool can be considered as a series of ovoids, which are commonly generated by the laser beam. Measurements molten pools were a total of 234 for both width and depth.

Results and discussion
The implementation of a volumetric heat source proposed by Gusarov et al. [42,43] resulted in the normal distribution of surface power Q 0 , which is presented in Fig. 8. For this condition, the laser power and diameter are 400 W and 70 µm, accordingly. Based on the rate of change of dynamic absorptivity dA/dz shown in Fig. 2 and power Q 0 , it is possible to calculate the power distribution through the powder bed thickness. The absorbed power distribution based on Gusarov model is observed in Fig. 9 for 10 µm (element size) increments in z-direction. Once the temperature histories are obtained with thermal effects from previous layers and based on the study of Luo and Zhao [41], the last layer of the powder bed was examined. Figure 10 shows the thermal distribution of the powder bed surface at spots 539 and 550, respectively. As observed, thermal histories developed during the simulation process did not reach the boiling point of 316L steel (3200 K).
For the measurement of the simulated melt pool, nodes with temperatures higher than melting point (≥ 1678 K) were selected and located in the laser irradiance region. This can be seen in Fig. 11, where a considerable melt pool depth is observed. Using this methodology, all three dimensions of the melt pool (width, length, and depth) were calculated in the last layer (132 molten pools). In previous work [31], a superficial heat source was studied to model the laser with constant absorptivity coefficient with several values of thermal conductivity anisotropic enhancement factors by design of experiments (DoF) to minimize the error of the calculated results and experimental data. It was concluded that although enhancement factors were optimized, the simulated molten pool geometries were far away from experimental ones when the model laser was a superficial heat source. For this reason, the investigation devised another laser model and the volumetric heat source with dynamic absorptance through bed powder thickness was employed. Table 3 compares both laser models, the mean, and standard deviation of the experimental and simulated data highlighting that difference is so close to the volumetric heat source with dynamic absorptance. Table 4 only compares the volumetric heat source and experimental data by the relative error when all simulated pool is considered.
Thermal histories were analyzed in those nodes located in critical zones, where a thermal accumulation caused by the scanning strategy (zigzag) was observed. Figure 12 shows the critical regions using red circles with a nomenclature from A to E. Thermal histories through time can be seen in Fig. 13 during the processing of last layer, where the second digit corresponds to layer number.
From Fig. 13 can be deduced that while the laser consolidates particles in the last layer, a thermal influence exists in nodes located in previous layers, obtaining a melt pool depth of up to 4 layers. The last layer can be seen as an example, in which thermal histories from A1 to A4 surpass the melting line. Therefore, the portion above the melting temperature indicates that the metal is in a liquid phase and the projection

Dynamic melt pool and effective width criterion
Following the dynamic phenomena discussed above, the change in melt pool dimensions throughout the manufacturing process was further investigated. The behavior of the melt pool width through the laser path in the last layer was analyzed as shown in Fig. 14, in which the horizontal axis corresponds to the absolute position of the laser (see Fig. 3) and vertical axis corresponds to the melt pool width. The first scan path initially showed an increase and transient behavior of the melt pool width until a permanent response was reached (~ 130 µm) due to the cooling of the powder bed. Subsequently, once the powder bed was thermally affected by the first scan path, all subsequent laser paths started with a maximum melt pool width due to the laser scanning strategy, but the behavior tends to decrease until reaching a steady state (~ 152 µm). This behavior was also reported by Weirather et al. [16] who recorded the changes of the melt pool length, depth, and width numerically determined by simulating the distance covered by the laser. Figure 15 shows the melt pool width evolution for the 1st and nth laser scan. This curve represents an average curve from the 2nd to 6th laser scan. It can be observed in Fig. 15 that an exponential regression was carried out in both curves to consider the thermal system as a transient and steady-state response. The equations shown in Fig. 15 are composed by two terms: the first term in the  [16] proposed the following net melting rate in relation to the laser path shown in Eq. (20), which can be defined in terms of volume.
where (dV/dt) f defines the volume of particles transforming from solid to liquid in a certain amount of time, whereas (dV/dt) s describes the volume rate associated with the solidifying liquid. If dV/dt > 0, the melt pool size increases. Contrarily, dV/dt < 0 relates to a shrinking size, whereas dV/dt = 0 points out the melt pool's steady state. Figure 16 shows the net melting rate supporting the sense of transient response at the beginning of the scan in both cases, although the nth scan is obtained a higher rate due to the thermal influence caused by the previous scan (i.e., nth − 1 scan).
Based on the melt pool width response in steady state discussed above and the fact that melt pool measurements were carried out in central areas of the samples, a new evaluation between numerical and experimental melt pool dimensions was proposed, in which the simulated widths were considered at t > 4τ. Table 5 shows an average width and standard deviation for all considered melt pools and those at time t > 4τ in nth scan and last layer. It was observed that for a criterion of melt pool width in steady state, a standard deviation of ± 7.21 µm and relative error (20)  where T f is defined as the alloy's melting point, T 0 represents the temperature of the preheated substrate, Q abs is the absorbed power by the alloy, k indicates the material's thermal conductivity, V is the laser scan speed, and α is the material's thermal diffusivity. The terms r and ξ refer to the melt pool size and shape, as shown in Fig. 17, in which ξ represents the gap from the heat source that is determined in the same trajectory in which the laser beam travels, and r = (ξ 2 + y 2 ) 0.5 indicates the gap to the melt pool edge from the center of the heat source. The latter, where the highest melt pool width is located, r = r* in accordance to Fig. 17.
In order to simplify the previous expressions, constants M and N are used and are defined next: Substituting in the Rosenthal equation and after mathematical deductions [49] yields However, Tang et al. [49] assumed that the term −ln(r * M) ∕r * M ≈ 0 , while, for purposes of this study, a value of 0.29 was estimated according to thermal properties and some process parameters summarized in Table 6.
Additionally, a new relationship between r* and φ* and the half-width was found.
Solving Eqs. (24) to (26) an estimated melt pool width of W = 150.93 µm is determined, which was graphed in Fig. 13 as a red and horizontal dashed line. There is a relationship between the estimated melt pool W and steady-state response w(t) in nth scan where this latter converges approximately in W and shows the permanent regime of w(t).   Table 6 Thermal properties of 316L steel at melting temperature and some process parameters

Conclusions
This research analyzed the modeling of a thermal process to calculate temperature histories and the behavior of a SS316L manufactured by SLM. A non-linear and transient FE model was applied to predict the thermal response of the material in each of its states (powder, liquid, and solid). The model included the process parameters that are involved in SLM as input data. For a deeper analysis, temperature histories in critical zones were monitored due to thermal accumulation caused by the selected scanning strategy (zigzag). This analysis was carried out in each layer in which there is no record of temperature above the boiling point. For the 4th and 5th layers, those critic zones showed several cycles of fusion/solidification, affecting up to four previous layers and lifetime of approximately ~ 1 ms of molten material.
A dynamic thermal response was identified and characterized, specifically on the melt pool dimensions. This dynamic behavior on the melt pool width was divided into two types of laser scans: the first and nth scans which behave differently and were separately analyzed by exponential regressions. Both curves showed a transient and steady-state melt pool width response whose response regimes were characterized numerically, and a relationship was found between the steady-state response w(t) of nth scan and the estimated molten pool width by the expressions derived from the Rosenthal equation.
From previous work, a new comparison between simulated and experimental melt pool dimensions was evaluated, in which the simulated widths were considered as long as the response was in steady state, in nth scan, and neglecting molten pool widths in the first scan. This comparison raised a new effective width criterion in order to more accurately relate calculated and measured dimensions of the melt pool. Funding The authors received financial support provided by the Institutional Fund for Regional Scientific, Technological and Innovation Development (FORDECyT) of the National Council of Science and Technology (CONACyT) of Mexico, through the project "Aeronautical Strengthening in Mexico's North East Region." Availability of data and material The data presented in this study are available upon request to the corresponding author. All the data is presented within the article.

Code availability
The code used to obtain the FEM results in this study is available upon request to the corresponding author.