Crack Extension Analysis of Atmospheric Stress Corrosion Based on Peridynamics

Abstract

Based on peridynamics, an atmospheric stress corrosion model was proposed. In this model, the role of hydrogen and stress in anodic-dissolution-dominated stress corrosion cracking was considered, and atmospheric corrosion was characterized by the change in liquid film thickness on the metal surface in the atmospheric environment. The near-field kinetic anodic dissolution model and the atmospheric corrosion model were coupled by varying the liquid film thickness. The thickness of the liquid film depended on factors such as the temperature, relative humidity, and hygroscopic salts. We validated the model using stress corrosion behavior from the literature for 304 stainless steel in a simulated atmospheric environment. The results of the model captured the crack expansion process. The obtained crack expansion direction and branching behavior agreed well with the experimental results in the literature.

1. Introduction

Atmospheric stress corrosion (ASCC) is the phenomenon of the stress corrosion cracking (SCC) of metals in the atmospheric environment and is one of the most prevalent and hazardous forms of metal corrosion. Metallic materials are subject to varying degrees of atmospheric corrosion from raw material inventory and component processing to assembly and service. Some common examples can be observed in land, water, and air transportation systems; electronic circuit boards; and urban and marine infrastructure. The degradation of metals due to atmospheric exposure has a significant impact on society [1]. According to recent reports, economic losses due to corrosion are estimated to be 3.4% of the global GDP, and losses due to atmospheric corrosion failure account for approximately more than 50% of losses in the overall corrosion field [2].

Stress corrosion cracking (SCC) is the process of rupture of a metallic material under the combined effect of stress and the corresponding environment. Stress corrosion cracking occurs in many common metallic materials in the corresponding environment, and in some harsh stress corrosion systems, stress corrosion in the material can even occur within a few minutes. The initial formation of SCC is microcracking, followed by the continuous expansion and fusion of microcracks to form macroscopic cracks. The environments in which SCC occurs are usually acidic, alkaline, and high-temperature water.

According to the corrosion mechanism, SCC is split into two types: anodic dissolution (AD) and hydrogen embrittlement (HE) [3]. Regarding the AD mechanism, under the combined influence of stress and the medium, the metal is dissolved by electrochemical reactions in the medium, resulting in crack creation and extension. The AD stress corrosion mechanism includes the slip dissolution mechanism, oxide film cracking mechanism, selective dissolution mechanism, and corrosion-promoted plastic deformation mechanism [4,5]. In the hydrogen embrittlement mechanism, hydrogen atoms diffuse into the metal matrix and become trapped by various hydrogen traps [6]. The metal matrix exposed to hydrogen becomes brittle due to reduced ductility, followed by cracking under stress. The hydrogen-enhanced debonding process, hydrogen-enhanced local plasticity, adsorption-induced dislocation emission, hydride production, cleavage, and other HE stress corrosion mechanisms are the most common [7,8,9,10].

During the atmospheric stress corrosion of metals, the variation of thin-film electrolyte thickness is a major factor in determining the corrosion rate and corrosion mechanism [11,12,13,14,15,16]. For example, the corrosion rates of iron [14,16] and steel [17,18,19,20] are significantly higher in uniform electrolyte thicknesses of 10–30μm than in monolithic electrolyte or immersion conditions [14,16,17,18,19,20], and the maximum corrosion rate of zinc occurs at electrolyte thicknesses of 8–10μm [21,22]. However, the thickness of thin-film electrolytes is affected by many factors, such as the temperature, relative humidity, wind, rainfall, and contaminants in the air and on the metal surface [23,24,25,26,27]. These environmental factors have an important influence on the thickness of the electrolyte, and this influence is both profound and extremely complex. The complexity is mainly due to the following reasons: these environmental factors are not constant and change dynamically, not only with time but also with spatial shifts; the effects of these environmental factors on the electrolyte thickness are not independent of each other but interact with each other. Furthermore, these environmental factors interact with each other and cause dynamic changes in each other.

ASCC has become an important research area for stainless steel because of the characteristics of sudden onset, large damage, and being difficult to monitor. Simulating the crack sprouting and expansion in materials plays an important role in the reliability analysis and health monitoring of engineering structures. Most of the existing numerical models for simulating crack emergence and expansion in materials are based on the finite element method (FEM). In the conventional finite element method, the displacement field discontinuity generated by the crack expansion is corrected by redefining the object. That is, the crack is defined as a boundary and the mesh needs to be re-defined after each crack extension. It is also necessary to provide mathematical expressions for the kinetic relationships that control crack growth for the simulated fracture method to describe how the crack expands under local conditions. Considering the above factors, it is difficult to solve the multi-crack problem with multiple interactions and a complex expansion mode using the traditional finite element method.

The extended finite element method (XFEM) [28] is one of the ways to solve these challenges by introducing additional nodal degrees of freedom and local strengthening functions on top of the standard finite elements. These reinforcement functions contain discontinuous displacement fields and can simulate the growth of cracks and fractures without the need to divide the mesh. However, this approach requires the reintroduction of external criteria. Meanwhile, the XFEM model performs poorly in terms of crack sprouting and fracture. Molecular dynamics simulation (MDS) [29] uses interatomic forces to simulate crack sprouting and extension, and is a detailed and realistic method for predicting material fracture. However, because of computational capacity constraints, the time and size scales that can be studied in simulations differ greatly from those required in real-world engineering applications.

Peridynamics (PD) was proposed by Silling et al. [30], who constructed new elastodynamic equations based on space integral equations and the idea of nonlocal action. The theory has unique advantages when dealing with damage and fracture problems by solving spatial integral equations to simulate mechanical behavior. It also has a strong ability to capture complex phenomena. Examples include crack nucleation, crack spreading, multi-crack joining, and crack halting in both isotropic and non-isotropic materials [31]. In PD theory, cracks can naturally sprout and expand freely without imposing a restriction criterion on crack expansion [32]. PD theory has been widely used in the simulation of crack initiation and extension in a variety of materials. PD corrosion models were shown to be a reliable and flexible tool for simulating corrosion damage and were used for damage evolution in pitting and intergranular corrosion in a variety of alloys [33,34,35,36].

Siavash Jafarzadeh et al. [36] coupled a PD corrosion model with a PD mechanics model to simulate the mechanochemical effects of stress on the anodic dissolution rate. Dennj De Meo et al. [31] proposed a numerical multi-physics field near-field kinetic framework for simulating the stress corrosion cracking of adsorbed hydrogen based on the AIDEC mechanism. Ziguang Chen et al. [37] proposed a PD model of stress corrosion that was dominated by an anodic dissolution mechanism. Xiu Ran et al. [38] simulated the hydrogen embrittlement of AISI 4340 steel. However, their simulation and experimental models were conducted by submerging the samples in an electrolyte solution and did not consider the stress corrosion cracking of metallic materials in an atmospheric environment, i.e., ASCC. Meanwhile, in their PD models, hydrogen embrittlement or anodic dissolution was singled out as the only mechanism for the occurrence of the stress corrosion of metals. However, according to relevant studies, anodic dissolution and hydrogen embrittlement do not exist in isolation when SCC occurs. The effects of anodic dissolution and hydrogen atoms are simultaneous and interact during stress corrosion [39,40,41,42,43,44], and the dominant mechanism of stress corrosion may also change when the environmental medium is changed [45,46,47]. Katona, R M et al. [48] experimentally confirmed the importance of hydrogen precipitation reactions during atmospheric corrosion.

There were two innovations in this study: (1) the application of PD in ASCC and (2) the introduction of a new PD model that organically combines the synergistic effects of anodic dissolution and hydrogen atoms in stress corrosion for the simulation of the ASCC of 304 stainless steel.

2. Materials and Methods

A peridynamics atmospheric stress corrosion model was proposed for the first time. In this model, the role of hydrogen and stress in anodic-dissolution-dominated stress corrosion cracking was considered, and atmospheric corrosion was characterized by the change in liquid film thickness on the metal surface in the atmospheric environment. The near-field kinetic anodic dissolution model and the atmospheric corrosion model were coupled by varying the liquid film thickness. We validated the model using stress corrosion behavior in the literature for 304 stainless steel in a simulated atmospheric environment.

2.1. Atmospheric Stress Corrosion Coupling Model with Peridynamics

2.1.1. Peridynamics Fracture Model

In the peridynamics theory proposed by Silling et al. [30], it is assumed that the material is composed of material points that are uniformly distributed in space. At any time, each material point in the material represents the position of a material particle, and each particle has its particle volume and mass density. Particles move and deform when they are affected by a force and velocity. Peridynamics involves the study of the motion of objects due to the interaction between particles in objects. It is assumed that the interaction force between a particle $x_i$ and a particle beyond a certain range (near field range δ) is negligible, and the particles within its near field range are called family $H_{x_i}$. PD theory studies the interaction between a particle and all particles within its near-field range.

The existence of interaction between the mass point $x_i$ and all the mass points in the family $H_{x_i}$ produces a force density vector that acts on the mass point $x_i$, which can be regarded as a force that is applied to the mass point $x_i$ by the mass point $x_j$. Similarly, the mass point $x_j$ is acted upon by all the other mass points in the family $H_{x_j}$. In the bonded peridynamics theory, a bond describes the interaction between two material locations, and a broken bond means that there is no longer an interaction force between the mass points. The mass points act in pairs, and the relative position state changes of two interacting material points when they are subjected to external forces are shown in Figure 1: $x_i$ and $x_j$ denote the initial positions of two material points, and $x_i^{\prime }$ and $x_j^{\prime }$ denote the positions of two mass points after a certain time t of external forces.

The fundamental equations of motion for near-field dynamics are constructed based on Newton’s second law. The mechanical parameters and characteristics of the materials are incorporated into the fundamental equations of motion through the instanton equations to achieve numerical simulations of damage to different materials. ${\mathit{u}}_i$ denotes the displacement of matter point $x_i$, where its unit is $m$, while ${\mathit{u}}_j$ denotes the displacement of matter point $x_j$, $m$. Let $\mathbf{\eta }$ and $\mathbf{\xi }$ represent the relative displacement and relative position of the material point ${\mathit{x}}_i$ and material point ${\mathit{x}}_j$, respectively, where the expressions are $\mathbf{\xi }={\mathit{x}}_j-{\mathit{x}}_i$ and $\mathbf{\eta }=\mathit{u}\left({\mathit{x}}_j,t\right)-\mathit{u}\left({\mathit{x}}_i,t\right)=\left({\mathit{x}}_j^{\prime }-{\mathit{x}}_j\right)-\left({\mathit{x}}_i^{\prime }-{\mathit{x}}_i\right)$, $m$. Within the theoretical framework of bonded peridynamics, the equation of motion of the material point $x_i$ in time $t$ can be formulated as [49]

$${\mathsf{\rho }}_i\left(x_i\right){\ddot{\mathit{u}}}_{ }{}_i\left(x_i,t\right)=\underset{H_{x_i}}{\overset{}{{\displaystyle \int }}}\mathit{f}\left(\mathbf{\eta },\mathbf{\xi },t\right)d{V}_j+{\mathit{b}}_i\left(x_i,t\right),\forall x_j\in H$$

(1)

where ${\mathsf{\rho }}_i$ denotes the mass density of matter point $x_i$, $kg/m^3$; ${\mathit{b}}_i$ denotes the external load in terms of the body force density (force per unit volume of matter), $N/m^3$; $H_{x_i}$ is the near-field range of matter point $x_i$; $V_j$ is the volume of matter point $x_j$, $m^3$; ${\mathsf{\rho }}_i\left(x_i\right){\ddot{\mathit{u}}}_{ }{}_i\left(x_i,t\right)$ represents the combined force that affects matter point $x_i$; and ${\int }_{H_{x_i}}\mathit{f}\left(\mathbf{\eta },\mathbf{\xi },t\right)d{V}_j$ represents the sum of the bond forces acting with all matter points in the near-field region $H_{x_i}$ of matter point $x_i$.

The force density vector $\mathit{f}\left(\mathbf{\eta },\mathbf{\xi },t\right)$ represents the inter-mass point-to-point force response function, i.e., the intrinsic force function, which is defined as the force vector per unit volume squared applied to the mass point $x_i$ by the mass point $x_j$, $N/m^6$. For austenitic steels, which undergo quasi-dissociative brittle fracture in a hydrogen environment, the initial microscopic elastic–brittle (PMB) model can be used, in which the relationship between the force density vector $\mathit{f}$ and the elongation $S$ is a combination of being linear and nonlinear, as shown in Figure 2.

Combining the results of Silling et al. [50] and Jianfeng Wang et al. [51], the present constitutive equation can be expressed as

$$\mathit{f} \left(\mathbf{\eta },\mathbf{\xi }\right)=\left\{\begin{array}{c}\left[C_{1}Sexp\left(-\frac{S-S_{ec}}{S_{oc}}\right)-C_2\right]\frac{\mathbf{\eta }+\mathbf{\xi }}{\left|\mathbf{\eta }+\mathbf{\xi }\right|}, S_c \le S<S_{ec}\\ \left(C_{1}S-C_2\right)\frac{\mathbf{\eta }+\mathbf{\xi }}{\left|\mathbf{\eta }+\mathbf{\xi }\right|}, S_{ec}\le S\le S_{et}\hfill \\ \left[C_{1}Sexp\left(-\frac{S-S_{et}}{S_{ot}}\right)-C_2\right]\frac{\mathbf{\eta }+\mathbf{\xi }}{\left|\mathbf{\eta }+\mathbf{\xi }\right|}, S_{et}<S\le S_t\\ 0, others\hfill \end{array}\right.$$

(2)

where $S_{et}$ denotes the linear elastic elongation of the bond in tension and $S_{ec}$ denotes the linear elastic elongation of the bond in compression, corresponding to the elastic deformation part of the material; $S_{ot}$ and $S_{oc}$ are the critical elongation of the bond in tension and compression, respectively, corresponding to the maximum value of the near-field force; $S_t$ and $S_c$ are the fracture elongation of the bond in tension and compression, respectively, corresponding to the fracture of the bond; S is the elongation of the bond, which can be interpreted as strain in classical continuum medium theory; $C_1$ and $C_2$ are peridynamic material parameters, $N/m^6$; $S$, $C_1$, and $C_2$ are calculated as shown below:

$$S=\frac{\left|\mathbf{\eta }+\mathbf{\xi }\right|-\left|\mathbf{\xi }\right|}{\left|\mathbf{\xi }\right|}$$

(3)

$$C_1=C=\frac{6E}{\pi {\delta }^4\left(1-2v\right)}$$

(4)

$$C_2=C\alpha$$

(5)

where $E$ is Young’s modulus, $v$ is Poisson’s ratio, the bond constant $C_1$ is equivalent to the “modulus of elasticity” of the material in elastodynamics, and $\alpha$ is the coefficient of thermal expansion of the material. The failure criterion of PD includes the critical bond elongation criterion and the critical bond energy criterion. The critical elongation $S_0$ is the extreme value of the bond deformation and can be solved based on the assumption that the fracture is completely separated and there is no other form of energy dissipation near the crack tip. It is calculated as

$$S_0=\sqrt{\frac{10G_0}{\pi c{\delta }^5}}$$

(6)

$G_0$ is the critical energy release rate, $N/m^{ }$, and is calculated as follows:

$$G_0=\frac{K_{IC}^2}{E/\left(1-{\nu }^2\right)}$$

(7)

Bond fracture occurs when $S\ge S_0{}_{}$. Any material point is impacted by some bonds within its near-field range. If the bonds acting on the material point in the region of the near field break, the damage increases. The damage of a material point $x$ at moment $t$ can be expressed as

$$\phi =\frac{n_d}{n}$$

(8)

where $n_d$ denotes the number of bonds that were broken in the near-field region of the material point and $n$ denotes the total number of bonds in the near-field region of the material point.

In many practical applications, cracks are present at the beginning or distributed in several parts of the structure; therefore, we use the PD method to prefabricate these initial cracks. This is also a common method used in studies related to near-field dynamics.

2.1.2. Atmospheric Corrosion Model

The atmospheric corrosion of metals is a phenomenon that is prevalent in nature; it is very different from the corrosion of metals in their native solution, with the main difference being the thickness of the electrolyte liquid film on the surface of the metal. According to the thickness of the liquid film, atmospheric corrosion is divided into three categories:

(1)

Dry atmospheric corrosion: There is no liquid film coating on the metal surface. It is mostly generated by pure chemical reactions and its corrosion efficiency and destruction are often modest.

(2)

Humid atmospheric corrosion: There is a 10 nm–1 $\mu m$ liquid film coating on the metal surface. This thin liquid film is generated on a metal surface via capillary action, adsorption, or chemical condensation. Atmospheric oxygen travels considerably more readily through the water film on the metal surface than through the liquid layer when totally submerged.

(3)

Wet atmospheric corrosion: There is a 1 $\mu m$–1 mm liquid film coating on the metal surface. It is broadly compatible with the electrochemical corrosion mechanism in the body solution.

The main research objects of atmospheric corrosion are humid and wet atmospheric corrosion. They fundamentally belong to the same category as corrosion in the native solution: electrochemical corrosion. During the atmospheric corrosion of metals, the variation of the film electrolyte thickness ($\lambda$) is the main determinant of the corrosion rate. The atmospheric corrosion rate increases with the liquid film thickness and then decreases. When the liquid film thickness is ${\lambda }_0$, the corrosion rate reaches a maximum and then gradually becomes smaller. This is mainly due to the liquid film being too thick, preventing the oxygen in the air from diffusing through the electrolyte solution to the metal surface where the anodic reaction occurs. At this time, the oxygen reduction reaction at the cathode begins to be blocked and the reduction reaction of hydrogen begins to occur.

The anodic reaction of atmospheric corrosion is

$$M⟶M^{n+}+n{e}^{-}$$

(9)

When the electrolyte layer is very thin, the cathodic process is dominated by the oxygen reduction reaction, and the hydrogen precipitation reaction can be disregarded. The cathodic process is

$$O_2+2H_{2}O+4e^{-}⟶4O{H}^{-}$$

(10)

For a higher electrolyte thickness of ${\lambda }_0$, the oxygen reduction reaction is blocked and hydrogen precipitation begins to occur in the cathodic process:

$$2H^{+}+2e^{-}⟶H_2$$

(11)

Van den Steen, N et al. [14] proposed a dynamic electrolyte film corrosion model based on the continuous medium assumption, i.e., the assumption that the liquid film on the metal surface is continuous. The model estimates the dynamics of the liquid film thickness with environmental changes by providing an initial liquid film thickness ${\lambda }_B$; it is based on a fundamental condensation/evaporation model of standard heat and mass transfer, it is driven by the known evolution of temperature and relative humidity with time, and it takes into account the effect of hygroscopic salts. The model uses the heat and mass transfer equation shown below:

Heat transfer equation:

$$q_{conv}=h\left(T_s-T_{\infty }\right)$$

(12)

Mass transfer equation:

$$\dot{m}=k_m\left(w_{\infty }-w\right)$$

(13)

Linking the heat and mass transfer coefficients through Lewis’ constant:

$$\frac{h}{k_m}=\left(C_{p,air}+w{C}_{p,water}\right)L{e}^{2/3}$$

(14)

where $q_{conv}$ is the heat flux at the surface/air interface, $W/m^{ }{}^2$; $h$ is the mass transfer coefficient; $T_s$ is the temperature of the film surface, $K$; $T_{\infty }$ is the temperature of the air, $℃$; $\dot{m}$ is the mass flux of water toward or away from the film surface; $k_m$ is the mass transfer coefficient, $kg/\left(m^2\cdot s\right)$; $w_{\infty }$ and $w$ are the absolute humidity and absolute saturation humidity of the air at the film/air interface, respectively, $kg/m^3$ ; $C_{p,air}$ is the heat capacity of air, $J/\left(kg\cdot K\right)$; $C_{p,water}$ is the heat capacity of water, $J/\left(kg\cdot K\right)$; and $Le$ is the Lewis constant. When $w_{\infty }-w>0$, the film thickness will increase; otherwise, the film will become thinner. By using the mass flux $\dot{m}$, the metal surface area $S(m^2$), the mass density $\rho$, and the substrate thickness ${\lambda }_B$, the transient average liquid film thickness of the metal surface can be estimated from the following equation:

$$\lambda ={\lambda }_B-\frac{\dot{m}}{S\rho }$$

(15)

Under the hygroscopic effect of hygroscopic salts or mixtures, the model describes the effect of the hygroscopic effect on the liquid film thickness by correcting the partial pressure of water due to the activity of water. The absolute humidity of air at the film/air interface $w_{\infty }$ and the absolute saturation humidity $w$ at the film/air interface in the model are calculated as shown below:

$$w_{\infty }=\frac{0.622p_{\infty }^{s}R{H}_{}}{p-p_{\infty }^{s}R{H}_{}}$$

(16)

$$w=\frac{0.622R{H}_{\infty }}{p-a_w{p}^s}$$

(17)

where $p$ indicates the atmospheric pressure, $Pa$; $p_{\infty }^s$ indicates the vapor pressure of air, $Pa; p^s$ is the saturation vapor pressure of air, $Pa;$ $a_w$ is the activity of water ; $R{H}_{}$ is the relative humidity of air; and $R{H}_{\infty }$ is the equilibrium relative humidity.

$$p_{\infty }^s=610.78exp\left(\frac{17.269T_{\infty }^{\prime }}{T_{\infty }^{\prime }+238.3}\right)$$

(18)

$$p^s=610.78exp\left(\frac{17.269T^{\prime }}{T^{\prime }+238.3}\right)$$

(19)

where $T_{\infty }^{\prime }$ is the air temperature expressed in degrees Celsius and $T^{\prime }$ is the surface temperature expressed in degrees Celsius.

In order to simplify the numerical model of atmospheric corrosion, according to the variation of the electrolyte thickness $\lambda$, the following assumptions are made qualitatively in this study: no anodic dissolution type of stress corrosion occurred in the dry atmospheric corrosion and the liquid film thickness ${\lambda }_0$ was taken as the critical point in tidal atmospheric corrosion and wet atmospheric corrosion. When $\lambda <{\lambda }_0$, there was only the reduction reaction of oxygen in the cathodic process, which involved ignoring the hydrogen precipitation reaction, and thus, ignoring the influence of hydrogen atoms on the anodic dissolution type of SCC. When $\lambda \ge {\lambda }_0$, the hydrogen precipitation reaction was introduced in the cathodic process, and the effect of hydrogen atoms on the anodic dissolved SCC was considered. Referring to [17], ${\lambda }_0$ was taken as 30 μm.

2.1.3. Peridynamics Anodic Dissolution Model

In the research model of PD, corrosion is considered a process of dissolution diffusion. F. Bobaru et al. [52,53] studied heat and mass transfer in discontinuous objects with evolving capabilities and proposed the following PD equation for mass transfer:

$$\frac{\partial C\left(x_i,t\right)}{\partial t}={{\displaystyle \int }}_{H_{x_i}}^{}k\left(x_i,x_j\right)\frac{C\left(x_j,t\right)-C\left(x_i,t\right)}{{||x_j-x_i||}^2}d{V}_{x_j}$$

(20)

Equation (20) simulates mass transfer by anodic dissolution at the metal interface and diffusion of metal ions in the electrolyte, where $C\left(x,t\right)$ is the concentration of plasmas $x_i$ at time $t$, $mol/m^3$. The concentration-dependent damage model describes the relationship between the damage and concentration change during corrosion; $k\left(x_i,x_j\right)$ is the diffusion coefficient of the bond between plasmas $x_i$ and $x_j$, which is defined as a function of mechanical damage $d\left(x,t\right)$; and $k$ is calculated as follows:

$$k\left(x_i,x_j\right)=\left\{\begin{array}{c}{K}_L\left(D\right), d\left(x_i,t\right)=1\cup d\left(x_j,t\right)=1\hfill \\ 0 , d\left(x_i,t\right)<1\cup d\left(x_j,t\right)<1\hfill \\ K_{diss}, \left[d\left(x_i,t\right)ord\left(x_j,t\right)\right]<1\cup \left[d\left(x_i,t\right)ord\left(x_j,t\right)\right]=1\hfill \end{array}\right.$$

(21)

Equation (21) indicates that if both masses belong to the solid phase, no diffusion occurs, and if one point is in the solid phase and the other in the liquid phase, the bond is an interfacial bond that carries the anodic dissolution microflux for which $K=K_{diss}$ and is called microdiffusivity. $K_L\left(D\right)$ is the microdiffusivity of the liquid, which can be calculated based on the diffusivity $D$ of the electrolyte:

$$K_L\left(D\right)=\frac{9D}{2\pi {\delta }^3}$$

(22)

The damage function $d\left(x,t\right)$ in Equation (21) is

$$d\left(x,t\right)=\left\{\begin{array}{c}1 , C\left(x_i,t\right)\le C_{sat}\\ \frac{C_{solid}-C\left(x_i,t\right)}{C_{solid}-C_{sat}}, C_{sat}<C\left(x_i,t\right)<C_{solid}\\ 0 , C\left(x_i,t\right)=C_{solid}\end{array}\right.$$

(23)

Equation (23) indicates that when the concentration of metal ions in the solid undergoes dissolution below the saturation concentration of metal ions in the electrolyte solution $C_{sat}$, the solid becomes a liquid, i.e., the damage is 1, and the damage to the intact solid is 0. The damage to the region can be determined from the concentration of metal ions in Equation (23), and this damage is expressed by removing the mechanical bonds in the corroded region.

Hydrogen is able to enter the specimen and promote anodic dissolution, and there is a synergistic effect of hydrogen and stress in promoting corrosion. Under the combined action of hydrogen and external stress, the metal anodic dissolution current density is

$$i_{\left(\sigma -H\right)}=k_{\left(H\right)}{k}_{\left(\sigma \right)}{k}_{\left(\sigma ,H\right)}{i}_0$$

(24)

Under external stress alone, the metal anode dissolution current density is

$$i_{\left(\sigma \right)}=k_{\left(\sigma \right)}{i}_0$$

(25)

In atmospheric corrosion, Equations (24) and (25) can be combined and rewritten as

$$i=\left\{\begin{array}{c}{k}_{\left(H\right)}{k}_{\left(\sigma \right)}{k}_{\left(\sigma ,H\right)}{i}_0, {\lambda }_{}<{\lambda }_0\\ k_{\left(\sigma \right)}{i}_{0,} \lambda \ge {\lambda }_0\end{array}\right.$$

(26)

where $i_{\left(\sigma -H\right)}$ is the dissolution current (proportional to the dissolution rate) generated by adding a constant load after hydrogen charging,$A$; $i_0$ is the corrosion current when no hydrogen is loaded, $A$; $k_{\left(H\right)}$ and $k_{\left(\sigma \right)}$ are the factors that increase the dissolution when hydrogen and stress are each present alone, respectively, and $k_{\left(\sigma ,H\right)}$ is the hydrogen and stress interaction factor. The related equations are shown below:

$$k_{\left(H\right)}=exp\left[\frac{∆U_H-T∆S_H}{RT}\right]$$

(27)

$$k_{\left(\sigma \right)}=exp\left[\frac{W{\sigma }^2}{2E\rho RT}+\frac{\eta W}{\left(1+n\right)P^n\rho RT}\left(n{\sigma }^{n+1}+{\sigma }_y^{n+1}\right)+\frac{∆\phi }{RT}\right]$$

(28)

$$k_{\left(\sigma ,H\right)}=exp\left[\frac{{\overline{V}}_H{\sigma }_h}{RT}\right]$$

(29)

Since the dissolution current at the anode is proportional to the dissolution rate, Equation (26) is rewritten as

$$k_{diss}^s\left(x_L,x_S\right)=\left\{\begin{array}{c} k_{diss}^u{k}_{\left(\sigma \right)}, {\lambda }_B<{\lambda }_0\\ k_{diss}^u{k}_{\left(H\right)}{k}_{\left(\sigma \right)}{k}_{\left(\sigma ,H\right)}, {\lambda }_B\ge {\lambda }_0\end{array}\right.$$

(30)

where $x_L$ and $x_S$ denote the liquid and solid ends of the interfacial bond, respectively; $k_{diss}^s$ is the microsolubility of an elastically stressed material; and $k_{diss}^u$ is the microsolubility of an inelastically stressed material.

2.1.4. Atmospheric Corrosion–Fracture Coupling Model with Peridynamics

In this part, the mechano-chemical coupled kinetic model for near-field kinetic ASCC is introduced. In the previous sections, the PD fracture model, the air corrosion process, and the anodic dissolution model were covered briefly. In the fracture model, the nodes carry information about the displacement ($u$) and mechanical strain damage index ($d_s$); in the atmospheric corrosion mechanism, they carry information about the liquid film thickness ($\lambda$) and corrosion mechanism ($c_i$); and in the anodic dissolution model, they carry information about the concentration ($C$) and corrosion damage index ($d_c$). Figure 3 illustrates the links between them.

Starting from the top of the figure, environmental changes lead to changes in the electrolyte film thickness, where the film thickness is updated according to Equation (15) and the corrosion mechanism is determined according to Equation (26). Starting from the corrosion evolution and moving in a clockwise direction: electrochemical corrosion leads to a phase change in the corrosion front of the metal according to Equation (23); then, the corrosion damage index is updated, the displacement and bond strain are updated according to Equation (2) under the current load, and finally, the microsolubility is updated according to Equation (30). In the counterclockwise direction: alterations in the corrosion process’s progression lead to the distribution of metal displacement under a mechanical force, and when the bond stretch exceeds the critical bond strain, mechanical damage occurs, and according to the corrosion damage connection in Equation (21), the cumulative mechanical damage accelerates corrosion.

We numerically simulated atmospheric stress corrosion cracking using the algorithm shown in Figure 4. It refers to the algorithm reported by Siavash Jafarzadeh et al. [36]. There were two main differences: the hydrogen atoms were involved in the anodic dissolution type of SCC and atmospheric stress corrosion cracking was formulated through the thickness variation of the electrolyte liquid film.

Solution steps: The region in the defined material was discretized with a uniform grid and the relevant parameters were entered for modeling. The initial concentration C, the initial liquid film base thickness ${\lambda }_b$, the relevant atmospheric parameters were set, and the initial damage d and the initial liquid film thickness $\lambda$ were calculated. Equation (1) was solved using the mechanics solver step and $k_{diss}^s$ was found using Equation (30). For the AD diffusion solver: the diffusion of metal ions at each time step was solved using the near-field kinetic diffusion model (Equation (20)), and the concentration C and the damage index d of the metal ions were updated. For the mechanics solver: when the metal ion diffusion step cycle was completed, the displacement of the material point was solved using Equation (1) and the damage d at the material point was updated.

2.2. Simulation of Atmospheric Stress Corrosion Fracturing in 304 Stainless Steel

In this section, we discuss the 3D simulations of the atmospheric corrosion of 304 stainless steel under controlled laboratory conditions using a coupled near-field kinetic atmospheric corrosion–fracture model in which artificial seawater was used, where stress corrosion of 304 stainless steel in an acidic chloride ion atmosphere is widely believed to be dominated by anodic dissolution, with the participation of hydrogen atoms in promoted anodic dissolution corrosion of the metal [54,55].

2.2.1. Experiment

To verify whether the coupled near-field kinetic atmospheric corrosion–fracture model proposed above could capture the crack extension of the atmospheric corrosion of 304 stainless steel, we compared the experimental results with the numerically simulated crack extension results.

Masami Mayuzumi et al. [56] placed 304 stainless steel in a laboratory at constant temperature (353 K) with adjustable humidity and dripped 30μLto 50μL of $MgC{l}_2$ artificial seawater on the surface of the specimen, applied a constant load during the experiment, and investigated the morphology of crack extension of the specimen at 35% relative humidity and 75% relative humidity (Figure 5).

Figure 5 shows the scanning electron micrographs of corrosion cracks on the specimen surface. The main purpose of this study was to investigate the crack expansion in atmospheric stress corrosion, and it was reasonable to prefabricate the cracks in the simulation experiments in order to simplify the model and improve the computational efficiency.

The chemical composition of 304 stainless steel is shown in

Table 1. Chemical composition of the specimen.

C	Si	Mn	P	S	Ni	Cr	Fe
0.06%	0.61%	0.95%	0.028%	0.011%	8.12%	18.10%	Bal.

.

2.2.2. Calculation Model Setup

To simplify the model and improve the computational efficiency, we set the numerical simulation model as a rectangular plate of $150\times 4\times 100{\mathrm{mm}}^3$. As shown in Figure 6, the initial crack length was 25 $\mathrm{mm}$ and a constant load was applied at the top and bottom edges of the plate. A layer of nodes to the right of the initial crack was used as the initial corrosion area, and the metal ion concentration in the initial corrosion area was set to zero and maintained at zero throughout the corrosion process. Zero-flux conditions were applied to the rest of the boundaries. The region describing the material was discretized using a uniform grid with nodes situated in the center of each cell, which was equivalent to employing Gaussian quadrature with a single point in the vicinity.

The region defining the material was discretized with a uniform grid with the nodes located in the middle of each cell, which was equivalent to using a single-point Gaussian quadrature in the neighborhood.

The computational model consisted of 60,000 material points of 1$mm$ size, with eight particle layers at the bottom and top edges as loading layers, with a grid spacing $\Delta x=1mm$ and a near-field radius $\delta =3\mathrm{mm}$ (corresponding to m = 3), which was uniformly discretized in the $x$, $y$, and $z$ directions with $\Delta x$. The Young’s modulus of 304 stainless steel $E=194020\mathrm{MPa}$ and the Poisson’s ratio $v=0.3$; we chose $C_{solid}=143\mathrm{M}$, $C_{sat}=5.1\mathrm{M}$, and ${\lambda }_0=30\mathsf{\mu }\mathrm{m}$, which are typical values for stainless steel, and applied $C_{sat}$ at the location of the crack tip and along the newly created crack node. We assumed that the solubility of the corrosion front was $0.02{\mathsf{\mu }\mathrm{m}}^2/\mathrm{s}$, and to improve the modeling rate, by referring to the calculation in [57], we used $k_H=1.07$, $k_{\sigma }=1.001$, and $k_{H,\sigma }=1.08$. Typically, the plastic deformation of steel must be accounted for in the fracture model. However, as the plasticity of stainless steel decreases in corrosive environments, the ductile fracture transforms into a quasi-brittle fracture, necessitating the use of the intrinsic force model of brittle material with linear and nonlinear mechanical behavior in the solution process.

In the process of stress corrosion, the time for corrosion diffusion is much longer than the time for crack expansion; therefore, the time t for corrosion diffusion calculation in the computational model should be much larger than the time $T$ for the mechanical solution step. Moreover, the given initial liquid film thickness ${\lambda }_B$ has a great influence on the results of the numerical simulation. Therefore, we conducted experiments to determine the appropriate time interval and initial liquid film thickness in the following section.

3. Results and Discussion

In order to study the effect of the corrosion diffusion time $t$ with mechanical loading time $T$ and the initial liquid film thickness on the numerical simulation of the stress corrosion process, we performed the simulation solution based on the data model with RH = 35% relative humidity under the condition that all the boundaries were not set to zero flux. We first studied the effect of the corrosion diffusion time t and mechanical loading time $T$. Four sets of simulations were done: T = 30t (40t, 50t, 60t) when t = 30 (40, 50, 60), i.e., the corrosion diffusion solver ran 30 (40, 50, 60) time steps after the fracture solver ran once. The results are shown in Figure 7. The red-brown part of the energy sketch (left) indicates the mechanical damage release energy. The red-brown part of the corrosion diffusion domain diagram (right) corresponds to the region where the metal ion concentration in the metal matrix became zero.

Figure 7 shows that in the case of a coupled solver with the same running time, the longer the corrosion diffusion solver ran, the greater the corroded region. The direction of fracture expansion was toward the region of corrosion diffusion in this investigation.

Taking ${\lambda }_B=0.5\mathsf{\mu }m \left(5\mathsf{\mu }m,50\mathsf{\mu }m,500\mathsf{\mu }m\right)$, $T=40t$, setting the initial corrosion region metal ion concentration to zero, and applying zero flux to the rest of the boundary, we did four sets of simulation experiments to determine the effect of the initial electrolyte film thickness, and the simulation results are shown in Figure 8. As shown in Figure 8: As ${\lambda }_B$ becomes larger, the area of corrosion diffusion also becomes larger, and is accompanied by the appearance of branching cracks. When ${\lambda }_B$ to increase to a certain extent, the same time step under the area of corrosion diffusion in turn becomes smaller. This is in line with the effect of electrolyte film thickness on corrosion rate.

Lu et al. [58] demonstrated that oxygen plays a key role in the nucleation and extension of stress corrosion cracks. For IGSCC cracks, oxidation first occurs before the SCC crack tip. The crack then expands by connecting with microcracks within the grain interface oxide. For TGSCC cracks, the nucleation and growth of microcracks in the oxide along the deformation zone is an important mechanism for TGSCC crack extension. This was consistent with the partially discontinuous cracks in Figure 8b,d,f,h in the simulation experiments, further justifying the model.

Through simulation experiments, we found that the thickness of the electrolyte liquid film not only affected the rate of stress corrosion to a great extent but also affected the sprouting and evolution of branching cracks during the stress corrosion process. Through the above two simulation experiments, we found that the selection of the appropriate time interval and the initial liquid film thickness had an important influence on the simulation solution.

Based on the above simulation experiments, the metal ion concentration in the initial corrosion region was set to zero and zero flux was applied to the boundary of the remaining prefabricated cracks. We used $T=40t$ and ${\lambda }_B=0.5\mathsf{\mu }m$ and numerically simulated the atmospheric stress corrosion of 304 stainless steel in a laboratory environment with a relative humidity RH of 35%. We used $T=40t$ and ${\lambda }_B=5\mathsf{\mu }m$ and numerically simulated the atmospheric stress corrosion of 304 stainless steel in a laboratory environment with a relative humidity RH of 75%. The results are shown in Figure 9 and Figure 10. We found that the results of the numerical simulation of the model and the experimental results displayed good agreement. Both were dominated by transverse cracks and had some similarities in the direction of branching cracks.

Nyrkova et al. [5], through the study of stress corrosion of X70 pipeline steel in a near-neutral environment, found that the stress corrosion process can proceed both according to the anodic dissolution mechanism and according to the hydrogen embrittlement mechanism. 304 stainless steel samples due to local anodic dissolution on the surface formed the initial corrosion areas, these initial corrosion areas are the original defects of crack generation. In these defect areas, will cause a large stress concentration. At the same time, the hydrogen atoms generated during the anodic dissolution process will be enriched in the stress concentration area, and the hydrogen atoms will diffuse through the defects to the metal interior under the effect of concentration gradient and stress gradient.

Chiari et al. [59], through a positron annihilation lifetime spectroscopy study on 304 stainless steel, confirmed that hydrogen diffuses into a deeper region after the application of tensile stress. Simultaneously, vacant hydrogen complexes are generated in this region and develop into vacant clusters. This study corroborated the role of hydrogen and stress in anodic-dissolution-dominated SCC.

Haruna et al. [60] studied the stress corrosion of iron-containing chloride ion rust layers that occurred in humidity-controlled air. The results of the study showed that the relative humidity affected the rate of hydrogen uptake by the rust layer of iron. In the region of 42–74% RH, the hydrogen uptake rate increased with increasing RH. The pH of the rust layer was weakly acidic in the range of 4.2–4.3 during the corrosion in the tested RH range. Therefore, the higher the relative humidity within the range of 42–74%, the more hydrogen the metal absorbed during the stress corrosion. The greater the stress gradient and concentration gradient (hydrogen diffusion in metal mainly occurs via stress diffusion), the greater the diffusion range of hydrogen atoms in the metal. Under tensile stress, the direction of the crack extension remained essentially the same as the direction of diffusion of hydrogen atoms. This was highly similar to our simulation tests.

Chengshuang Zhou et al. [61], in their study of X80 pipeline steel hydrogen embrittlement, found that when hydrogen embrittlement occurs in this material, the hydrogen originated from internal hydrogen and surface-absorbed hydrogen. When the material was charged with hydrogen, hydrogen atoms were first adsorbed on the metal surface, and then hydrogen penetrated the metal and became internal hydrogen. When the hydrogen charging was stopped, the absorbed hydrogen on the surface disappeared and the internal hydrogen atoms diffused outward. We set the boundary conditions with zero concentration of hydrogen atoms outside the model, i.e., only hydrogen atoms existed inside the model. In this case, the diffusion direction of hydrogen should be toward the right outer boundary of the model under the stress gradient and concentration gradient. In Figure 9 and Figure 10, the results of our numerical simulations were consistent with the hydrogen atom diffusion law of Chengshuang Zhou et al. [61], which showed that our model was reasonable.

4. Conclusions

In this study, a coupled model of atmospheric stress corrosion was proposed, which was based on the near-field kinetic theory, near-field kinetic anodic dissolution theory, and atmospheric corrosion theory. In the near-field kinetic anodic dissolution theory, the factors of hydrogen, stress, and the hydrogen and stress synergy are added to express the role of hydrogen and stress in stress corrosion that is dominated by anodic dissolution. It should be noted that in our model, the role of hydrogen atoms was to promote the anodic dissolution of the metal rather than to produce hydrogen embrittlement. In the theory of atmospheric corrosion, the influence of the thickness of the electrolyte liquid film on the rate and mechanism of atmospheric corrosion is the main consideration. The variation in the liquid film thickness was calculated using the main factors that affect the variation of liquid film thickness. Then, a coupled near-field kinetic anodic dissolution and atmospheric corrosion model was proposed. This coupled model enabled the numerical simulation of stress corrosion of metals that occur in an atmospheric environment.

We used the atmospheric stress corrosion of 304 stainless steel as a control experiment. In the simulation experiments, it was found that ASCC cracks started to sprout from the prefabricated cracks and expanded toward the boundary of the sample under the combined effect of stress and hydrogen atoms, and the relative humidity had an important effect on the rate and direction of crack expansion. Within a certain relative humidity range (35–75%), the increase in relative humidity accelerated the corrosion rate. The results of the numerical simulations were in good agreement with the experimental results.

Because our model simplified the role of hydrogen and stress in the anodic dissolution process, it did not consider the influence of the material microstructure, and the calculation of the liquid film thickness was relatively simple in the atmospheric corrosion model. Therefore, the numerical simulation results had some errors relative to the experimental results.