Finite Element Based Physical Chemical Modeling of Corrosion in Magnesium Alloys

Abstract

Magnesium alloys have found widespread applications in diverse fields such as aerospace, automotive, bio-medical and electronics industries due to its relatively high strength-to-weight ratio. However, stress corrosion cracking of these alloys severely restricts their applications in several novel technologies. Hence, it will be useful to identify the corrosion mechanics of magnesium alloys under external stresses as it can provide further insights on design of these alloys for critical applications. In the present study, the corrosion mechanics of a commonly used magnesium alloy, AZ31, is studied using finite element simulation with a modified constitutive material damage model. The data obtained from the finite element modeling were further used to formulate a mathematical model using computational intelligence algorithm. Sensitivity and parametric analysis of the derived model further corroborated the mechanical response of the alloy in line with the corrosion physics. The proposed approach is anticipated to be useful for materials engineers for optimizing the design criteria for magnesium alloys catered for high temperature applications.

1. Introduction

Research in magnesium (Mg) alloys has found widespread interest in recent years due to their attractive strength-to-weight ratio, lightweight and high stiffness. These extraordinary properties make Mg alloys particularly useful for applications requiring high performance and energy savings such as structural engineering, automotive and biomedical applications [1,2,3]. It should also be noted that the weak corrosion resistance of Mg alloys restricts their applications in several novel fields [4,5,6,7]. Most of the existing studies in Mg alloys are focused on investigating the mechanical strength based on specific applications such as structural, mechanical or bio-medical engineering. Hence, there is a need to study the mechanics of Mg alloys under corrosive conditions. It is also well known that Mg alloys are more susceptible to stress corrosion cracking (SCC). Hence, an effective deployment of Mg alloys in new and emerging fields will require an understanding of the behavior of corroded Mg alloys under mechanical and thermal loading conditions.

Most of the literature studies on Mg alloys are based on their mechanical properties for specific applications. Argade et al. [8] studied the mechanism of AZ31 Mg alloy under slow strain rate testing technique. They found that the corrosion of AZ31 is predominantly influenced by its microstructure which results in a favorable texture for increased adsorption of hydrogen. Su et al. [9] deployed a controlled phosphate deposition technique for improving the surface biocompatibility of AZ60 Mg alloy. They found that the coating process is influenced by temperature and pH, which affects the corrosion resistance of AZ60 alloy. Choudary et al. [10] investigated the SCC mechanism of various aluminum free Mg alloys in a simulated human fluid environment. It was found that the SCC of Mg alloys in the simulated condition were attributed to their electrochemical and microstructural characteristics. Gaur et al. [11] developed a silane based biodegradable coating for improving corrosion resistance of Mg alloys in physiological environment. They found that formation of a stable and uniform hydroxide layer on alloys surface facilitated adhesion of silane coating. Wei et al. [12] investigated the corrosion behavior of plasma formed coatings on Mg alloys with respect to the aluminum composition in the alloy mix. It was found that the difference in porosity level of coatings is caused by different discharging behaviors in Mg alloys.

In addition to experimental studies, some simulation studies have also been conducted to investigate the corrosion mechanics of Mg alloys. Grogan et al. [13] developed a physical corrosion model based on finite element (FE) technique for predicting mechanical strength of corroded metal stents. The developed model was able to capture the change in surface due to corroding environment. Duddu et al. [14] developed an extended FE model for modeling galvanic corrosion in magnesium alloys. It was demonstrated that the proposed approach can be used to model crevice geometry without remeshing. Shang et al. [15] deployed molecular dynamic simulations to study the adsorption behavior of composite coatings on Mg alloys. The presented approach provided an effective application based surface modification of Mg alloys. In addition to FE modeling, analytical models based on computational intelligence techniques have also been widely used for predicting corrosion. For instance, Narimani et al. [16] deployed artificial neural network for predicting the effects of chemical composition and corrosion characteristics in magnesium alloy pipelines. Similarly, Zadeh Shirazi and Mohammadi [17] used a hybrid neural network model for predicting corrosion rates in underwater pipelines.

It is clear from the literature studies that experiments and simulation have been deployed to study corrosive characteristics in Mg alloys, with simulation being more preferred mode of investigation. This is because simulation techniques provide a viable alternative in saving cost related to time and materials. The current state-of-the-art in scientific modeling of corrosion is directed at two fronts, viz. mechanical and chemical model (Figure 1). These models are capable of generating accurate solutions in predicting corrosion characteristics with minimal cost and high rapidity. In addition to these science based models, process modeling based on computational intelligence technique has also been conducted, as can be seen from the previous literature studies. The advantages of process modeling include formation of an explicit mathematical relationship which can help in understanding hidden relationship between various input parameters in corrosion process. It can be seen that the science based models and process models have their own unique advantages. As a result, it will be useful to integrate the models which can combine the advantages of accuracy of FE models with model formation ability of process modeling.

Hence, the main objective of the present research proposes an integrated science based model and process model for predicting corrosion characteristics of Mg alloy (Figure 2). The Mg alloy considered in this study is the AZ31 alloy used in aircraft components, biodegradable metallic implants, cell phone and laptop cases, etc. In the present study, a corroded AZ31 alloy with cracks and voids typically resulting from an environment that leads to SCC (e.g., hydrogen embrittlement) is modeled using FE simulation. The corrosive film factor is included in the study by considering the Pilling–Bedworth ratio (P–B ratio), which is defined as the ratio of the volume of corrosive metal oxide to the volume of the metal. As can be seen in Figure 2, the study focuses on understanding the mechanical strength of AZ31 alloy with respect to the input parameters viz. the P–B ratio, strain rate and temperature. As can be seen in Figure 2, the FE model is used to predict the mechanical strength of AZ31 alloy in corrosive environment with respect to the input parameters, viz. the P–B ratio, strain rate and temperature. The process modeling based on genetic programming (GP) is the used to obtain the mathematical relationship between the mechanical strength based on the input parameters. Scientific analysis of the formulated model is also undertaken based on actual corrosion science which will provide vital design inputs for maximizing the performance of AZ31 alloys in corrosive environment. The various approaches deployed in investigating the corrosion characteristics in AZ31 alloys and the proposed research objective is depicted in Figure 1 and Figure 2, respectively.

2. Finite Element Modeling of Corroded AZ31 Alloy

2.1. Geometry Description and Boundary Conditions

The FE modeling of AZ31 alloy is performed using ABAQUS/Explicit Version 6.14 software (Dassault Systemes, Vélizy-Villacoublay, France, 2014) with fully coupled thermal stress analysis. The thermal analysis enables modeling the strength characteristics of AZ31 alloy at elevated temperatures. The model geometry constructed in ABAQUS consists of AZ31 alloy of dimensions 450 mm × 20 mm × 10 mm. The most predominant byproduct of corrosion in Mg alloys is magnesium oxide (MgO). This corrosive product is included in the geometry by adding a film of MgO on top of AZ31 alloy. In the ABAQUS/CAE model, one end of the model geometry is fixed with an ENCASTRE boundary condition that arrests all rotational and translational degrees of freedom. Translational degree of freedom along one axis is provided to the other end of the model geometry for tensile loading condition. The description of AZ31 model geometry with corrosive film along with the applied boundary condition is depicted in Figure 3. The material properties of AZ31 alloy and the corrosive oxide film used in the FE modeling is illustrated in

Table 1. Properties of alloy and film used in FEM simulation [18,19]. FEM: Finite Element Modelling.

Parameter	AZ31 Alloy	MgO Film
Thermal Conductivity (k)	96 W/m-K	42 W/m-K
Density (ρ)	1770 kg/m3	3580 kg/m3
Young’s modulus (E)	44.8 GPa	249 GPa
Poisson’s ratio (υ)	0.35	0.18
Specific heat (Cp)	1000 J/kg/K	877 J/kg/K

. The defined model geometry is then meshed into discrete eight-node (CPE8RT) type elements, following which boundary conditions are applied to create mechanical loading of the AZ31 part.

2.2. Constitutive Modeling of Voids in AZ31 under Corrosion

AZ31, similar to most magnesium alloys, is susceptible to SCC. SCC results in the formation of mild cracks and voids on the surface of the material. These cracks can nucleate and grow further under application of external loads which results in catastrophic failure of the material. This phenomenon in which the nucleation and growth of cracks and voids that lead to ductile failure of metal is modeled in FE simulation using the Gurson–Tvergaard–Needleman (GTN) model [20]. Hence, in this work, the AZ31 alloy subjected to mechanical loading is modeled using the GTN damage model with integrated Yld2000 yield criterion [21,22]. The model is programmed in ABAQUS solver using the user-defined subroutine module VUMAT. This approach can provide a better prediction of crack evolution and failure of AZ31 alloys in corrosive environment. The yield function of GTN damage model is described mathematically as [23],

$$\mathsf{\Phi }=\frac{q^2}{{\mathsf{\sigma }}_y^2}+2q_1{f}^{*}\cos h\left(\frac{-3q_{2}p}{2{\mathsf{\sigma }}_y}\right)-\left(1+q_3{f}^{*2}\right)=0$$

(1)

where σ_y is the equivalent stress; p is hydrostatic stress; q is the effective stress; and q₁, q₂, and q₃ are fitting parameters to improve prediction of void interaction effects in the Gurson damage equation.

The reduction in mechanical strength of the material due to void coalescence is defined by the damage parameter $f^{*}$, expressed as a function of void volume fraction f. For $f^{*}=0$, the GTN yield condition becomes the normal yield condition and $f^{*}=1$ indicates that the material has failed. It is defined as,

$$f^{*}=\{\begin{array}{cc}f& f\le f_{\mathrm{c}}\\ f_{\mathrm{c}}+\mathsf{\delta }\left(f-f_{\mathrm{c}}\right)& f_{\mathrm{c}}<f\le \\ {\overline{f}}_F& f\ge f_{\mathrm{F}}\end{array}{f}_{\mathrm{F}}$$

(2)

where,

$$\mathsf{\delta }=\frac{{\overline{f}}_{\mathrm{F}}-f_{\mathrm{c}}}{f_{\mathrm{F}}-f_{\mathrm{c}}}, {\overline{f}}_{\mathrm{F}}=\frac{q_1+\sqrt{q_1^2-q_3}}{q_3}$$

(3)

where $f_{\mathrm{c}}$ is critical void volume fraction at which void coalescence first occurs, and $f_{\mathrm{F}}$ is the void volume fraction at which the material fails due to crack propagation.

The rate of increase in total void volume fraction $\dot{f}$ is due to growth of existing voids as well as nucleation of new voids, $\dot{f_{\mathrm{new}}}$. This can be described as,

$$\dot{f}=\dot{f_{\mathrm{growth}}}+\dot{f_{\mathrm{new}}}$$

(4)

The growth rate of existing voids is proportional to the hydrostatic component of plastic strain rate (${\dot{\mathsf{\epsilon }}}_{kk}^p$), given by,

$$f_{\mathrm{growth}}=\left(1-f\right){\dot{\mathsf{\epsilon }}}_{kk}^p$$

(5)

The nucleation rate of new voids is generally defined by strain-controlled nucleation rule which assumes that the voids undergo nucleation as second phase particles which is normally distributed for the total particle population. This is described as,

$$\dot{f_{\mathrm{new}}}=\frac{f_N}{S_N\sqrt{2\mathsf{\pi }}}\exp \left[-\frac{1}{2}\left(\frac{{\overline{\mathsf{\epsilon }}}^p-{\mathsf{\epsilon }}_N}{S_N}\right)\right]{\overline{\mathsf{\epsilon }}}_m^{pl}$$

(6)

where $f_N$ represents the volume fraction of void-nucleating particles, and ${\mathsf{\epsilon }}_N$ and $S_N$ are the average and standard deviation of the strains at which particles nucleate voids, respectively.

The steps involved in the modified GTN damage model using VUMAT subroutine is illustrated in the form of a flowchart [22] in Figure 4.

The modified GTN model is used to predict the strength characteristics of AZ31 alloy by varying the P–B ratio (x₁), temperature (x₂) and strain rate (x₃). The range of tested input parameters tested are given in

Table 2. Input parameters for tensile loading simulation of AZ31 alloy.

Process Input Variable	Values	Units
P–B ratio	0.2, 0.3, 0.4, 0.5	No units
Temperature	300, 350, 400, 450	K
Strain rate	0.0001, 0.001, 0.01	s−1

. The simulation computes the engineering stress and engineering strain, from which the true stress and true strain are calculated using classical mechanics relations.

2.3. Validation of Finite Element Model

The stress–strain plot of the AZ31 alloy under uni-axial tensile loading is shown in Figure 5. For comparison of FE results with actual data, the experimental tensile loading characteristics of AZ31 obtained from Ref. [24] are also included in the figure. It can be seen from this plot that the FE model was able to predict the stress characteristics of AZ31 with an accuracy of 92%. It is noted that the model predicted a relatively high yield stress compared to experiments. This can be attributed to several factors such as the material defects, temperature fluctuations, machine settings in tensile tester, etc. Furthermore, the FE geometry of AZ31 also factored the corrosive film component which may have caused certain variations in the yield stress prediction. This clearly demonstrates the ability of the FE simulation for generating required data sets that can be subsequently utilized by optimization algorithms.

The tensile loading simulation performed by ABAQUS/CAE software is depicted in the form of screen shots in Figure 6. It can be seen that the stress distribution in the AZ31 alloy as well as the corrosive film depends markedly on the simulation temperature due to the development of thermal stresses in the material. A total of 5000 steps were used to complete the tensile loading simulation of the AZ31 material in the CAE model.

3. Analytical Modeling of Stress Corrosion Characteristics of AZ31 Alloy

The FE model generates the strength data of AZ31 alloy by systematically varying the input parameters related to P–B ratio, temperature and strain rate. A total of 48 data sets are then used to obtain non-linear analytical model using GP technique. In GP cluster, the training data are used to train the GP model and the validity of obtained model is tested on testing data. Kennard-and-Stone algorithm is integrated within the GP algorithm to ensure randomness of data selection for training and testing sets. The GP cluster is then used to derive mathematical model that describes the relationship between the true strength of the alloy (y) as a function of input variables, viz. P–B ratio (x₁), temperature (x₂) and strain rate (x₃). The descriptive statistics of the data sets generated from the FE model is shown in

Table 3. Descriptive statistics of true stress data from FE simulation.

Parameter	P–B Ratio, x1 (No unit)	Temperature, x2 (K)	Strain Rate, x3 (s−1)	True Stress, y1 (MPa)
Mean	0.35	375	0.034	137.43
Median	0.35	375	0.001	121.92
Standard deviation	0.11	56.49	0.047	49.30
Kurtosis	−1.38	−1.38	−1.53	−0.49
Skewness	−1.842 × 10−15	0	0.730	0.75
Minimum	0.2	300	10−4	69.82
Maximum	0.5	450	0.1	235

.

The accuracy and effectiveness of the evolutionary algorithm framework depends on parameter settings, which control the generation of different sizes of mathematical models. A trial and error approach is used by the authors to determine the optimum parameter settings in the present study. The optimal settings determined are listed as follows: the population size is 400, the number of iterations is 20, the number of maximum genesis 6, the probability for crossover is 0.80, and the number of generations is 200. The authors deployed the GP algorithm from the GPTIPS toolbox in MATLAB framework developed by Searson et al. [25]. The best GP model (Equation (7)) for the true stress of AZ31 alloy (Y_AZ31) is selected based on the minimum mean absolute percentage error value among all the runs.

$$Y_{\mathrm{AZ}31}\left(\mathrm{MPa}\right)= 68.18+\left(1036.91\times x_1\right)+\left[0.71\times \tan \left\{\left(x_2\right)\times \left(\log \left(\log \left(x_3\right)\right)\right)\right\}\right]-\left[920.58\times \left\{\log \left(\cos \left(\log \left(x_2\right)\right)\right)\right\}\right]-\left[1130.62\times \left\{\tan \left\{\left(\tan \left(x_1\right)\right)\right\}\right]+\right[0.53\times \left\{\tan \right(\tan \left(\log \left(\log \left(x_3\right)\right)\right))\}]+[0.49\times \{x_2\times (\tan \left(\tan \left(x_1\right)\right))\}]$$

(7)

where $x_1,x_2,x_3$ are values of P–B ratio, temperature (K) and strain rate (s⁻¹), respectively. $Y_{\mathrm{AZ}31}$ is the value of true stress derived from the model based on the values of $x_1,x_2,x_3$.

3.1. Performance Evaluation of the Analytical Model

The evaluation of the FEM-GP model performance is conducted by analyzing the statistical metrics of the model as defined by the root mean square error (RMSE), the coefficient of determination (R²), and the mean absolute percentage error (MAPE). These metrics provides an indicator of model performance as they measure the deviations of the predicted from actual values. The definition of these model metrics represented by mathematical equations is available in Ref. [26]. The performance and validation analysis of the FEM-GP model is shown in Figure 7. It can be seen from the figure that the FEM-GP model is able to predict the true stress of AZ31 with close fit and higher accuracy. Hence, it can be concluded that the FEM-GP model can be used as an effective approach for generalizing the true stress of AZ31 alloy under corrosion.

3.2. Model Analysis Based on Corrosion Physics

Corrosion of material is a complex process in which several factors exhibit considerable influence on the material behavior. Hence, any analytical model describing the material strength of the component under corrosion must be able to capture the actual physics involved during corrosion. For this, parametric analysis is conducted which shows the variation in the true stress based on the considered input factors. Figure 8 shows the effect of each of the considered individual factors on the true stress of AZ31 alloy under corrosion. It can be seen that as the P–B ratio increases the material strength degrades rapidly. The degradation is more pronounced when P–B ratio exceeds 0.3. Based on corrosion physics, it can be stated that the corrosive film does not protect the alloy from degradation. Similar trend can be seen for the case of temperature in which a corrosive AZ31 alloy exhibits lower material strength at higher working temperatures. This indicates that AZ31 alloy under corrosive conditions is not suitable for high temperature applications. However, the strain rate seems to have no effect on the material strength of AZ31 alloy which agrees well with previous experimental study [24]. Hence, at ambient temperatures, the effect of loading (i.e., strain rate) does not have any negative or positive impact on the strength of the AZ31 alloy. This corroborates with the fact that material failure under SCC occurs not because of external load itself, but due to the crack growth and nucleation of voids that results from residual stresses, loading, etc. Hence, from this analysis it can be concluded that the derived analytical model was able to confirm well with the observed physics in corrosion process.

A sensitivity analysis was also conducted to determine the relative importance of each of the input parameters on the considered output. The sensitivity analysis is defined by partial derivative of strength with respect to temperature. While the partial derivative is computed, the other inputs such as the P–B ratio and strain rate are kept constant at their mean values. The effect of the two inputs P–B ratio and strain rate, which are at mean values, is still there on the strength. Therefore, the strength depends on the inputs in a complicated way. The relative importance of the considered input parameters on the mechanical strength of AZ31 alloy under corrosion is depicted in

Table 4. Sensitivity analysis of true stress of AZ31 alloy.

Input Parameter	Percentage Contribution to True Stress (%)
P–B ratio	36
Temperature (K)	62
Strain rate (s−1)	2

. It can be seen that the most important amongst those studied under corrosion is the temperature, followed by P–B ratio and strain rate. This analysis corroborates with actual corrosion mechanisms in magnesium alloy system, as they are more susceptible to SCC at elevated temperatures. As AZ31 is primarily used for high temperature applications, controlling the system temperature is crucial to avoid tensile strength degradation under corrosion, which can lead to material failure.

4. Conclusions

The present work introduced an integrated FE based computational intelligence approach for modeling corrosion mechanics of AZ31 alloy under mechanical loading. FE modeling of AZ31 alloy with a modified constitutive damage model was used to predict the true stress responses of AZ31 alloy with varying corrosive film concentration, temperature and strain rate. The data obtained from the model were used to derive a genetic programming model, which could predict the strength of the alloy with an accuracy of 92% and confirmed the data trend with actual corrosion physics. Sensitivity and parametric analysis of the genetic programming model showed that the mechanical strength of the AZ31 alloy highly dependent on system temperature and does not show much variation with regards to applied strain rate. The proposed model can be used to determine optimal design parameters for AZ31 alloy under corrosive environments without the need of experiments. This can help in reducing costs associated with material fabrication and testing, which can also lead to sustainable materials design.