## Abstract

This paper is a synthesis of earlier results supplemented by new results to define a comprehensive analysis of the growth rate of stress corrosion cracking (SCC). Two mechanisms, anodic dissolution (AD) and hydrogen embrittlement (HE), have been considered to calculate the SCC growth rate of AA 7050-T6 for a surface-breaking crack with blunt tip in an aqueous environment. The relative contributions of each mechanism and their mutual interactions have been quantitatively assessed. Results show that AD provides critical conditions for HE, which explains in part a stepwise propagation of the crack. Finally, the total crack growth rate due to the combined effects of AD and HE has been determined, and numerical results have been compared with experimental data, and a calculation of the crack growth rate for a practical configuration has been presented.

## 1. Introduction

Stress corrosion cracking (SCC) is an important form of fracture, because it can cause significant deterioration of the integrity of structures under mechanical loading conditions that would be quite safe without the presence of a corrosive environment. SCC is one of the most frequent failure mechanisms in the aeronautical, chemical, nuclear and offshore industry. In particular, high-strength 7XXX-series aluminium alloys present a continuing problem for the aeronautical industry. Significant efforts have been devoted to understanding SCC, but the prediction of the SCC rate is still a challenging subject because of the complex nature of the synergistic interactions among the environment, the materials and the stress that involve electrochemical and mechanical processes.

The system of interest in this paper is a surface-breaking crack in a structural aluminium alloy (AA 7050) in an aqueous environment, where anodic dissolution (AD) and hydrogen embrittlement (HE) are considered to be the principal mechanisms of damage [1–3]. This system is chosen because it widely occurs in the aeronautical industry and local chemistry measurement data are available [3]. To develop a comprehensive crack growth model, many factors need to be considered. Among the most important factors are: (i) the relationship between the bulk environment and the crack tip environment; (ii) the AD rate expected from the electrochemical kinetics for a given crack chemistry and the overpotential for hydrogen production; (iii) the relationship between the overpotential from hydrogen production and the internal hydrogen, or the hydrogen discharge current density and the hydrogen permeation current density; and (iv) stress-assisted and trapping-affected hydrogen diffusion ahead of the crack tip.

The link between the crack chemistry and the bulk chemistry has been established by the transport equation, which is important because different transport rates of each species can cause a crack chemistry which is quite different from the bulk environment [3–5]. Ions in an electrolyte move in response to an electric field (a process called *migration*) and also in response to concentration gradients (*diffusion*) and bulk fluid motion (*convection*) [6]. As SCC has no cyclic motion of the crack wake surface, as in fatigue [7,8], the convection is assumed to be negligible. The Nernst–Planck equation with the assumption of electroneutrality has been applied to evaluate the chemistry and electrochemistry at the crack tip [9–12]. At ambient temperature, concentrated solutions will most likely arise because of high solubility of metal cations [12]. To consider the interaction of ions [13,14] in a concentrated solution, the modified form of the equation that features an activity term in it has been used in this work [6,11]. The geometry of the crack has an important role in crack chemistry [7,8]. The geometrical effects can be captured by using the modified one-dimensional model proposed by Turnbull [12] or the two-dimensional model which has been used in this work. Chemical reactions have been treated by the widely used assumption of local equilibrium of chemical species [15]. Local equilibrium has been obtained with the well-established geochemical reaction code PHREEQC, which has the ability to solve thermodynamic equilibrium in a mixed aqueous/solid-phase system with homogeneous and heterogeneous reactions [16]. It is assumed that the crack solution is fully deaerated because of rapid depletion of *O*_{2} during the oxygen reduction reaction. On the surface near the crack mouth, deaeration affects crack chemistry and electrochemistry [12], which has, however, not been considered here. Two other factors, the formation of hydrogen bubbles and the precipitation of solute species, have been introduced to adjust the transport of species, and they jointly result in an effective diffusivity and mobility [17–19].

The boundary electrochemical kinetics is treated using the Tafel equation. The dependence of electrochemical kinetics on the crack chemistry has been obtained based on the bare metal kinetics experiment conducted by Ford *et al*. [20]. Gangloff *et al.* [21] have obtained the relationship between the overpotential for hydrogen production and the internal hydrogen for a Ni-based superalloy. Similarly, we have related the hydrogen flux at the crack tip to the hydrogen discharge current density based on the relation that the permeation current density is proportional to the square root of the discharge current density, as proposed by Iyer *et al.* [22].

The reported hydrogen diffusivity shows a large variability due to the effect of trapping at grain boundaries, dislocations and vacancies. According to Young & Scully [23], vacancies have the most influence on the effective diffusivity while grain boundaries have minimal influence in aluminium and its alloys. In this work, dislocations and vacancies have been considered in hydrogen transport analysis. The effect of local plasticity on the density of trap sites has also been considered [24]. To model the crack tip stress and strain fields, a blunted crack tip related to J2 plastic deformation [25–27] has been considered in the paper. It is well known that hydrostatic stress in front of the crack dilatates the lattice, which attracts hydrogen atoms. Hydrogen-enhanced decohesion has been considered as a mechanism for crack growth due to HE with the model proposed by Katz *et al.* [28]. Their model provides a reasonable prediction of the toughness degradation due to hydrogen localization. The total crack growth rate is determined based on the superposition model which combines Faraday's law for AD and a fracture mechanics model for HE which has been derived from the fracture criterion.

The purpose of this paper is to present a computational framework to determine the SCC growth rate. In future work, the results will be used as a critical building block for the prognosis of SCC.

## 2. Mathematical model of stress corrosion cracking

### (a) Coupling procedure

Before describing the mathematical models for each mechanism, it is of interest to draw a big picture which shows how the problem has been solved. As it is a coupled multi-physics and multi-mechanism problem, multiple software packages have been used. For the calculation of crack chemistry, PHREEQC (interfaced by IPhreeqc) and COMSOL have been used, and stress-assisted and trapping-affected hydrogen diffusion has been solved by ABAQUS and its user subroutines. To overcome the difficulties related to stiff equations in the discretized time domain of the numerical model, the two simultaneous and dependent processes that are multi-component transport and reaction are separated using the sequential non-iterative approach, also called the operator splitting or time-splitting approach [15].

Figure 1 illustrates the program structure for the coupling of COMSOL, ABAQUS and IPhreeqc using Matlab as the controlling application and interface for the three programs. The first column and the third column describe the calculation of crack chemistry and the crack growth rate due to AD. At the beginning of the simulation, the program initializes COMSOL and IPhreeqc computations and specifies boundary conditions. Precise concentration values for initial and boundary conditions in COMSOL have been determined by preliminary PHREEQC speciation calculations. After the initialization procedure, computations enter the iteration loop of the split operator, where COMSOL and IPhreeqc procedures alternate for the calculation of transport (TRANSPORT OF SPECIES) and reactions (BALANCE *c*_{i} AND GET pH). The second column indicates the procedure for the hydrogen transport and the crack growth rate due to HE. It is assumed that the crack is opened before exposure to the aqueous environment. Initial opening is simulated by fully coupled deformation and hydrogen transport analysis, and subsequent hydrogen transport analysis has been conducted under different boundary conditions based on the crack tip hydrogen flux *N*_{H} calculated from crack chemistry analysis. Further discussion will be provided in the model verification section.

### (b) Crack chemistry

#### (i) Governing equations

The modified form of the Nernst–Planck equation for moderately dilute electrolytes has been used to simulate mass transport inside a crack [6]. We have
*N*_{i} is the flux of species *i*, *D*_{i} is the diffusivity of species *i* and *c*_{i} is the concentration of species *i*. Also *γ*_{i} and *γ*_{n} are the activity coefficient of species *i* and a chosen ionic species *n*, respectively, *z*_{i} and *z*_{n} are the valences of species *i* and *n*, *u*_{i} is the mobility of species *i*, *F* is Faraday's constant and *ϕ*_{s} is the difference in the electrode potential, *E*, in the electrolytic solution (** ϕ_{s}**=

**(**

*E**x*=0)−

**(**

*E**x*)). The current density can be written as

In the case of an aqueous environment, the area near the tip is filled with gas bubbles generated from an electrochemical reduction reaction such as a hydrogen evolution reaction (HER), as observed by Singh *et al.* [29]. Thus, the crack solution can be viewed as a porous medium. To solve the mass transfer in a porous medium, averaging needs to be performed in a volume element within the electrode [6,19,30,31]. The averaged flux of species *i* can be expressed by using the volume fraction of solution, *ϵ*, with equation (2.1). Let *c*_{i} be the solution-phase concentration of species *i*, averaged over the pores. The superficial concentration is thus *ϵc*_{i}. Then
*D*_{i,eff} is the effective diffusion coefficient of species *i*, and *u*_{i,eff} is the effective mobility of species *i*. In equation (2.5), the volume fraction is inside the gradient operator, because it is the function of the position inside the crack [31]. The effective conductivity of the electrolyte is defined as
*κ*_{0} is the conductivity without bubbles, and *f* is the volume fraction of bubbles inside the electrolyte, *ϵ*=1−*f*. According to equations (2.6) and (2.7), the bulk values of diffusion coefficients and mobilities are multiplied in the paper by *ϵ*^{0.5} to yield *D*_{i,eff} and *u*_{i,eff} values corrected for tortuosity caused by hydrogen bubbles. The actual tortuosity and constriction of a crack are not considered in this work, but it is expected that they have a finite effect on mass transfer.

From the conservation of species and charge, we obtain
*n* coupled equations (equations (2.8)–(2.10)) has been solved for *c*_{i} and ** ϕ_{s}**.

To get the distribution of the volume fraction of gas bubbles, the crack is divided into several volume elements, and the fraction of hydrogen bubbles is the ratio of the volume of hydrogen bubbles and the volume element. It is assumed that the distribution of hydrogen bubbles follows the normal distribution function:
*x* is the coordinate along the crack, *μ*_{x} is the mean of the normal distribution function, which is adjusted to match the volume of hydrogen bubbles, and

#### (ii) Activity coefficient model and transport properties of electrolyte solutions

The calculation of the activity coefficient has been done by PHREEQC (figure 1 ‘BALANCE *c*_{i} AND GET pH IN SOLUTION’). As mentioned earlier, it is expected that the solution in the crack is relatively concentrated, so that the deviation from ideal may be significant. Specific ion interaction theory (SIT) is used to get an activity coefficient, which accounts for binary interactions of anions and cations with interaction parameters that are specific for individual ions in chloride or other media [33]. It is known that SIT is valid up to ionic strength of 3.5 mol kg^{−1}. The general SIT equation for a single ion, *i*, can be written as
*D*, the Debye–Hückel term, equals
*z*_{i} is the charge of ion *i*, and *m*_{j} is the molality of major electrolyte ion *j*, which is of opposite charge to ion *i*. Interaction parameters, *ϵ*_{ij}(*I*), refer to the interaction between ion *i* and major electrolyte ion *j* at ionic strength

Transport properties in the governing equations depend on composition and concentration. To take the dependence into account, they are updated in every loop (figure 1 ‘UPDATE TX PROPERTIES’). As the three properties, namely diffusivity, mobility and conductivity (or conductance), are interrelated, only conductance has been calculated in this work. The Turq–Blum–Bernard–Kunz (TBBK) conductance model [13] is used to get the composition and concentration-dependent conductance. In the TBBK model, the Fuoss–Onsager continuity equations were solved directly by a Green's function technique using the mean spherical approximation (MSA) pair distribution functions for the unrestricted primitive model (different ionic sizes). The single electrolyte solution consists of two types of free ions, cation (_{1}) and anion (_{2}), with total conductance given by
*i*th ion conductance. In the above equation, *δX*/*X* is the free ion relaxation force correction.

Equation (2.15) is valid for systems containing only a single cation and a single anion. Sharygin *et al.* [14] reviewed the literature on mixture models and found several different recommendations for essentially the same model although for different properties and with varying refinements. The consensus expression of the conductivity for a solution with a molar ionic strength *I* is
*M* and all negative species *X* with *Λ*_{MX} calculated at molar ionic strength *I*. The equivalent fractions of the species in the solution are given by *f*_{M}=*c*_{M}*z*_{M}/*c*_{eq} and *M*th cation (λ_{M}) is obtained by averaging over all anions that exist in the mixture, i.e.
*X*th anion, λ_{X}, is obtained by averaging over all cations, i.e.

The (simple) MSA theory has been found to be in agreement with experimental values for dissociated electrolyte solutions up to 2 mol l^{−1} [34]. The only unknown parameters are the diameters of the ions. The hydrated diameters of ions have been used in this work.

### (c) Hydrogen transport

Hydrogen is considered to reside at either normal interstitial lattice sites (NILSs) or trapping sites, e.g. grain boundaries, dislocations and vacancies. Young & Scully [23] observed the effect of trapping on diffusivity, in which the large influence of vacancies at relatively low concentrations is significant (e.g. 1×10^{−11} vacancies per lattice site lowers the diffusivity by almost a factor of 10 at 25°C), while large amounts of cold work (more than 10×10^{11} lines cm^{−2}) are needed for dislocation traps to appreciably affect the effective diffusivity. No evidence of hydrogen trapping was found at grain boundaries in either the pure aluminium or the Al–Zn–Mg alloy. According to first-principle calculations done by Lu & Kaxiras [35], a large amount of H atoms (up to 12) can be trapped at a single vacancy, which overcompensates the energy cost to form the vacancy defect. Therefore, current work only considers hydrogen trapping at dislocations and vacancies. According to Oriani's theory [36], the trapping site occupancies for dislocations *θ*^{(v)}_{T} are always in equilibrium with the lattice site occupancy,
*j* denotes *d* for dislocation or *v* for vacancy, *W*^{(j)}_{B} is the binding energy of hydrogen to the corresponding traps, and *R* is the universal gas constant equal to 8.31 J K^{−1} mol^{−1}. The hydrogen concentration in the NILSs is given by *c*_{HL}=*βN*_{L}*θ*_{L}, where *N*_{L}=*N*_{A}/*V*_{M} denotes the number of solvent atoms per unit volume, *V*_{M} is the molar volume of the host lattice, *N*_{A} is the Avogadro constant and *β* is the number of interstitial sites per solvent atom. The hydrogen concentration in trapping sites is calculated as *α*^{(j)} denotes the number of trapping sites per trap of type (*j*), and *N*^{(j)}_{T} is the corresponding trap density. The trap densities at dislocations, *N*^{(d)}_{T}, and vacancies, *V*_{H} is the partial molar volume of the hydrogen atoms in solid solution, *p*=Tr(*σ*)/3 is the hydrostatic stress, *ϵ*_{p} is the effective plastic strain and *D*_{eff} is an effective diffusion coefficient related to the NILS diffusion coefficient *D*_{H} by

### (d) Fracture criterion and crack growth rate model

In a purely brittle material, a crack will propagate when the stress intensity factor exceeds the intrinsic Griffith toughness, given by
*γ* is the required energy for the creation of new fracture surfaces, *E* is Young's modulus and *ν* is Poisson's ratio.

In systems where dislocations easily nucleate and propagate, the resulting plastic zone screens the applied mechanical load. When a dislocation nucleates, its stress field increases the critical value *K*_{Ic}. By superposition, Lin & Thomson [37] have obtained the stress intensity factor at the crack tip in the presence of a dislocation ‘cloud’. Based on their result, Huang & Gerberich [38] have proposed the relationship between *K*_{Ic} and *k*_{IG} which is given by
*σ*_{ys} is the yield stress, and *β*′ and *α*′′ are constants from a computer simulation of the dislocation structure around the crack tip.

To determine the reduction of the intrinsic fracture resistance due to H-decohesion, Katz *et al.* [28] coupled H localization with the dislocation model of the stressed crack tip to predict the threshold of *K*_{Ic} for hydrogen environment-assisted cracking as
*α* is a parameter related to the effectiveness of H decohesion. At the crack tip, hydrogen concentration lowers *K*_{Ic}, which facilitates the crack to propagate. The corresponding fracture criterion is

The crack growth rate due to HE is derived from the fracture criterion. A surface-breaking crack of depth *a*(*t*), under uniform distant tensile loading, *σ*, is considered. In linear elastic fracture mechanics, the stress intensity factor of a surface-breaking crack is known as
*x* at the crack mouth, the coordinate moving along with the crack tip may be expressed as *x*=*a*(*t*)+*ξ* and the time is defined as *τ*=*t*+*η*. A similar transformation has been used by Choi & Chudnovsky [39]. The crack grows when the stress intensity factor reaches the hydrogen-affected critical stress intensity factor, i.e.

The crack growth rate due to AD can be calculated from Faraday's law as
*j*, *i*^{(j)}_{a(tip)} is the anodic current density at the crack tip, *M*_{j} is the molar mass, *z*_{j} is the valence and *ρ*_{j} is the density.

## 3. Comparison of model predictions with experimental data

### (a) Electrolyte

At the start of the simulation, the solution contains 0.5 M Na_{2}CrO_{4}+0.05 M NaCl. The dissolution at the crack tip introduces alloying elements which undergo hydrolysis as follows:
_{3}, which is active if the solution becomes supersaturated with aluminium hydroxide [40]. In PHREEQC, this can be considered by adjusting the saturation index.

The reactions assumed to occur at the metal/solution interface are
*i*_{a} of Al^{3+}, 0.025*i*_{a} of Mg^{2+} and 0.025*i*_{a} of Zn^{2+}.

### (b) Electrochemical kinetics

Many protective anodic films are oxides and hydroxides whose dissolution depends upon the hydrogen ion concentration, and whose rate follows the Freundlich adsorption equation [42],
*k* is a constant, and *n* is related to the valency of the cations: *n*=0.5, 0.33, 0.25 for monovalent, divalent and trivalent ions, respectively. Equation (3.12) is used to get the crack wall current density. Due to passivity, AA 7050 shows a faster dissolution rate over the breakdown potential [43]. If the surface film is fully removed, the rate converges to the bare surface reaction rate [20]. In this paper, it is assumed that the tip is fully exposed to the environment, because the passive film is broken by plastic deformation and low pH. The data from Ford *et al*. [20] do not show much variation in the bare surface corrosion potential, *E*_{corr}, and the current density at corrosion potential, *i*_{corr}, with respect to pH. But the current density depends on pH due to repassivation, so that an equation similar to equation (3.12) is used to describe the dependence of the current density on the pH at the crack tip. The experimental data can be represented by
*E*_{corr} and *i*_{corr} are the corrosion potential and the corrosion current density on bare surface, respectively, *E*_{rp} and *i*_{rp} are on repassivating surface, and *β*_{a} and *β*_{rp} are the Tafel slopes for anodic reactions on bare surface and repassivating surface, respectively. *i*_{rp} is the function of pH as shown in equation (3.12). The electrode potential, *E*, can be calculated from
*ϕ*_{m} is the potential of AA 7050-T6, and *ϕ*_{m,film} is the potential over the film. The cathodic current density is
*i*_{cl} is the diffusion-limited hydrogen evolution current density which can be estimated from
*n* is the charge number of the cell reaction and *δ* is the diffusion layer thickness. The necessary parameters (shown in table 1) to construct the polarization curve have been obtained from the literature [1,20,43] where aluminium specimens are exposed to NaCl solution with various pH.

When the crack pops, it is quickly covered by a passive film, which causes an exponential decrease in current density. Therefore, it is more realistic to use transient reaction kinetics or time-averaged current density with the rate of production of bare metal [44]. According to the result obtained from this work, however, using bare surface kinetics is reasonable, because HE will only be activated when the crack tip chemistry reaches a critical condition. More importantly, the time interval between each pop is decreasing due to increasing hydrogen flux. As a result, the time to form the passive film is very limited.

### (c) Boundary condition

At the start of the simulation, oxygen is assumed to be depleted in the crack. The concentration outside the crack maintains the bulk concentration. Boundary reactions follow the electrochemical kinetics presented in the previous section. The dissolution of alloy elements has been obtained from the anodic current density, which is the sum of the contributions of each element. As the material is assumed to be homogeneous, the contribution ratio of each element is the same throughout the simulation. At the crack wall, the current density is three orders lower than at the the crack tip. The hydrogen discharge current density has also been obtained from the polarization curve. For stress-assisted hydrogen diffusion simulation, the flux boundary condition is used. The boundary conditions that are used in the simulations of HE are either the Dirichlet condition or a prescribed concentration boundary concentration [26]. Based on their study, Turnbull *et al.* [45] support the use of the flux boundary condition, especially in materials with high hydrogen diffusivity.

The rate constants for hydrogen adsorption, absorption, discharge (H^{+}+e^{−}+M→MH_{ads}) and recombination (2H_{ads}→H_{2}) can be obtained from a knowledge of the steady-state hydrogen permeation current, cathodic charging current, hydrogen overvoltage and hydrogen diffusivity [22]. When the steady-state permeation current *F* is Faraday constant and *μ* is the proportionality constant. In this work, *μ* is obtained from a calibration which matches the calculated crack growth rate with the experimental crack growth rate. In addition, the flux boundary conditions are chosen so that the resulting hydrogen concentration near the crack tip is within the experimentally measured concentration of hydrogen at the crack wake (0.007 to 0.30 wppm) [46].

### (d) Results

#### (i) Crack chemistry

The crack chemistry analysis has been worked out using a two-dimensional model with the geometry shown in figure 2*a*, where *a*=17 mm is the crack length, *h*_{m}=150 μm is the crack mouth opening displacement and *h*_{t}=2.4 μm is the crack tip opening displacement which corresponds to the applied stress intensity factor *K*_{I}=12.5 MPa. *h*_{t} has been obtained from transport–mechanical analysis using the geometry shown in figure 2*b*. A detailed explanation is given in §3d(ii). The governing equations are solved using COMSOL and PHREEQC, as explained in §2a. The two-dimensional model can effectively capture the dependence of crack chemistry on geometry, because the change in geometry changes the volume of the reaction vessels, which eventually modify the concentration of each species. The domain is discretized by quadratic triangular elements.

Figure 3*a* shows the changes in pH, [Cl^{−}] and *E*_{App}=−0.495 V_{SCE}. After an initial pH drop, a gradual decrease in pH has been observed, which is followed by a sharp drop when the chromate

The electrode potential within the local cell is one of the key parameters. In figure 3*b*, the crack tip potential shows a two-stage behaviour, the first one at 100 s is caused by the rapid drop of pH, which increases the current density sharply. The second drop, which shows a gradual decrease, is caused by the evolution of hydrogen bubbles and the precipitation of gibbsite (Al(OH)_{3}). This is consistent with experimentally observed data that a constriction formed by hydrogen bubbles and solid corrosion products results in a significant potential drop [47–49]. The calculated conductivity of an electrolyte solution ranges from 8.5 to 22 S m^{−1}. The conductivity measured from concentrated (but not saturated) crack solutions ranges from 8 to 16 S m^{−1} [50]. The difference is caused by the fact that a change in viscosity at higher concentrations entails changing ionic mobility, which is a manifestation of altered solvation structure around ions and short-range interactions between ions [51]. Due to the formation of hydrogen bubbles, the effective conductivity at the crack tip is as low as 2 S m^{−1}. Bubbles adhering to an electrode surface insulate a fraction of the surface that becomes inactive, and the actual current density is increased [52]. This means that the net flux of metal ions is the same as the one without bubbles. As a result of the decreased conductivity, the potential near the crack tip drops. It is also observed that the crack tip potential is a function of the fraction of hydrogen bubbles. As the fraction of coverage increases, the potential decreases. Figure 4*a* shows the distribution of crack potential, which matches well with the experimental result. The discrepancies are caused by the distribution of hydrogen bubbles and gibbsite.

Figure 4*b* shows the crack tip current density versus time. The anodic current density shows a sharp increase at 100 s, which corresponds to the drop in pH. Then it gradually decreases until 150 s. As the crack growth rate due to AD is proportional to the anodic current density, the shape of the curve will be the same for the crack growth rate by AD. Its contribution is diminishing with time. The cathodic current density keeps increasing until it reaches a steady-state condition. With the condition of low pH and potential calculated from the crack tip chemistry, the HER also occurs actively. As depicted in the cathodic polarization curve (equation (3.15)), significant amounts of cathodic current density can flow at the bare surface of the aluminium alloy. So from the result of the simulation, it is expected that the contribution of HE on the crack growth rate at lower pH cannot be ignored.

The limitation of the current AD simulation is that it does not consider the effect of crack propagation. Displacements during propagation affect the gas accumulation process and increase convective transport within the crack [53]. The distribution of hydrogen bubbles also requires further study. The transport of bubbles is affected by other phases, aqueous and solid, and surface roughness inside the crack. The nature of this process is probabilistic, therefore a probabilistic consideration will be introduced in future work.

#### (ii) Hydrogen transport to the region in front of the crack tip

Hydrogen transport analysis has been carried out using the small-scale yielding approach for mode I opening of the crack, as shown in figure 2*b*. The initial crack tip opening displacement is *h*^{0}_{t}, and the far-field radius of the model is taken as *et al*. [54]. A user subroutine, UMAT, has been used to calculate the effective plastic strain, the dilatational hydrogen-induced deformation and the hydrostatic stress. The material is assumed to obey the von Mises yield criterion with an associated flow rule. Another user subroutine, UMATHT, accepts these inputs to solve the governing equation (equation (2.22)). The 4-node plane-strain element, CPE4T, has been used. The number of elements was taken as 475, and the number of nodes as 516. The crack tip is modelled as a smooth notch to provide a good approximation of smooth blunting of a sharp crack at a higher deformation level.

The hydrogen transport analysis consists of multiple steps. It is assumed that the crack is opened before exposure to the aqueous environment. In step 1, hydrogen inside the computational domain is in equilibrium with atmospheric hydrogen. Coupled hydrogen transport–mechanical analysis under zero flux boundary conditions has been carried out. In step 2, the crack is filled with solution, which causes electrochemical reactions on the crack surfaces. Constant flux boundary conditions based on the crack chemistry have been considered during this step. After each step, except step 1, propagation of the crack is checked. If the crack has propagated, step 2 will be repeated with the step 1 result as the initial conditions. If the crack has not propagated, step 2 will also be repeated with the current result as the initial conditions until the fracture criterion is satisfied. The flowchart shown in figure 5*a* describes this process.

In step 1, a far-field boundary condition of the displacement, *u*_{i}, corresponding to an elastic *K*_{I}-field is applied incrementally over the outer layer of elements,
*E*=71.7 GPa, Poisson's ratio is *ν*=0.33 and *θ* is the angle of a point on the boundary. The target stress intensity factor is *h*_{t}, is three times greater than the initial value, *h*^{0}_{t}, requires blunting to three times the initial notch radius [55]. The uniaxial stress–strain law is given by the Ramberg–Osgood relationship of the form *ϵ*=*σ*/*E*+0.002×(*σ*/*σ*_{Y})^{n}, where the yield stress *σ*_{Y} is equal to 530 MPa, and the Ramber–Osgood parameter is *n*=10.

The lattice diffusion coefficient of hydrogen at 300 K is ^{−1} for dislocations and 68.6 kJ mol^{−1} for vacancies [23]. The interstitial hydrogen expands the lattice isotropically and its partial molar volume in solution is V_{H}=1.7×10^{−6} m^{3} mol^{−1}. The molar volume of aluminium is V_{M}=10.00 cm^{3} mol^{−1}. The uniform hydrogen concentration in the unstressed lattice, *c*_{0}=1.00×10^{3} atoms m^{−3} in equilibrium with the ambient environment, is used as the initial condition. The parameter *α* is taken to be equal to 1. The parameter, *β*, is set equal to 1, which corresponds to a maximum NILS concentration of hydrogen atoms per solvent lattice atom. The trap density at dislocations is expressed as a function of the dislocation density, *ρ*^{(d)}, and the lattice parameter, *a*=4.05×10^{−10} m: *et al.* [24], who have used the one-parameter approach of Kocks [56],
*b*_{0} is a constant, and *b*_{1} is the annihilation coefficient which depends on strain rate and temperature. However, *b*_{1} can be assumed to be constant for aluminium in the temperature range between 200 and 600 K. Equation (3.19) means that the dislocations are formed by the applied stress, *σ*, using *ρ*^{(d)}_{eq}=3.4×10^{16} m^{−2} is the equilibrium dislocation density and *ϵ*_{p}=0. The total vacancy concentration is the sum of the concentration in thermal equilibrium, *c*_{v,eq}, and the excess concentration, *c*^{ex}_{v}: *et al.* [57] has been used,
*k* is the Boltzmann constant equal to 1.38 J K^{−1}. Detemple *et al.* [24] have reported that as high as 10% of the total vacancy concentration is caused by plastic deformation, and they have proposed a model to predict the excess vacancy concentration. Due to the high binding energy at vacancy trap sites and the high excess vacancy concentration even with low strain rate (1×10^{−4} s^{−1}), the determination of the equilibrium concentration using Oriani's theory was impossible. In this work, only the vacancy concentration with the thermal equilibrium has been used, and the excess vacancy concentration will be considered in future work.

To use the criterion described in §2, it is required to get the parameters *β*′, *α*′′ and *α*. Gangloff [58] has used experimental data for high-strength steel and superalloys, and obtained a good fit for the two parameters. In this work, the parameters obtained from the experiment done by Jang *et al.* [59] using high-strength aluminium (AA 8090) have been used. Their experiment was designed to observe the effect of microstructure on the critical stress intensity factor. Their results show that the decrease in the threshold stress intensity factor is caused by the population of trap sites adjacent to the crack path at different heat treatment conditions. To obtain the total hydrogen concentration at the crack tip, computer experiments have been done at different charging times from 12 to 48 h. The result has been ploted in figure 5*b*. There are three unknowns, and only one set of experimental data is available. It has been assumed that *α*′ and *β*′′ are of the same order as for steel (*α*′′=0.0002 MPa m, *b*, *β*′ is determined as *α*′′ is 0.00025 MPa m, while *α* is obtained as *α* for steel. This is caused by the fact that we do not consider the excess vacancy concentration, which has high binding energy. In addition, pile-up of dislocations near obstacles such as grain boundaries or particles results in a heterogeneous distribution of trap sites. Therefore, 80 can be viewed as a correction factor for vacancy and dislocation concentrations.

For a different loading time, the normalized hydrostatic stress *σ*_{kk}/*σ*_{0} and the hydrogen concentration *N*^{(d)}_{T}/*N*_{L} are plotted, respectively, in figure 6*a*,*b* on the axis of symmetry (*θ*=0) ahead of the crack tip at the end of loading (*b*. Severe plastic deformation at the crack tip causes the formation of dislocation trap sites. The number of dislocation trap sites is saturated near the crack tip. As the hydrogen concentration is in equilibrium with the concentration in the NILSs, the location of the maximum concentration is not at the crack tip where the NILS concentration is lower.

In the steps other than step 1, the crack surfaces are exposed to an aqueous environment, which causes electrochemical reactions. In these steps, the flux boundary condition described in §3c is applied at the crack tip. Figure 7*a* shows *K*^{H}_{Ic}/*K*_{Ic} at the crack tip versus time. The critical value of the stress intensity factor *K*_{Ic} decreases with the increase in the crack tip hydrogen concentration, as expressed by equation (2.26). When the fracture work at the crack tip satisfies *b*). Therefore, the crack will stop eventually. At that point, an intermediate stress intensity factor can be expressed by *K*_{Ic} again, and crack propagation resumes. This provides a partial explanation of discontinuous crack propagation. During the process, the crack tip electrochemical kinetics are affected by the formation of a new surface and film. In this work, however, these effects have not been considered. Experimentally, the discontinuous crack growth behaviour is evidenced by crack arrest marks, stepwise changes in the direct current potential drop or strain gauge measurements [61].

The limitation of the current hydrogen diffusion analysis is that it does not consider the excess vacancy concentration due to plastic deformation. Parameters obtained from AA 8090 experimental data for hydrogen-affected fracture toughness may not be applicable because Al–Li alloys have higher solubility and different types of trap sites. It is also required to use more accurate boundary conditions, which should be justified either theoretically or empirically. In addition, Gangloff [62] has stated that the consideration of strain gradient plasticity is necessary.

### (e) Crack growth rate and comparison with experimental result

The crack growth rate due to AD and HE has been calculated. The contribution from AD is obtained by Faraday's law, and the crack growth rate for HE is obtained from fracture mechanics. The crack growth rate calculation shows that initially there is no contribution of HE because of low cathodic current density. After the crack tip, pH becomes low enough to activate the HER, and the contribution of HE to the crack growth rate can be observed. Figure 8*a* shows the calculated crack growth rate with respect to pH. The calculated (d*a*/d*t*)_{AD} well represents the crack growth rate for pH>3, and for pH lower than that (d*a*/d*t*)_{HE} presents a better fit. This indicates that the two different mechanisms act dominantly in different pH regions. The reason for the dominance of HE at lower pH is that the HER is more likely at lower pH and lower potential. For high hydrogen flux, the time to get the critical fracture condition is shorter than for low hydrogen flux. As the critical stress intensity factor *K*^{H}_{Ic} is the function of the crack tip hydrogen concentration, the numerator of the HE crack growth rate model (equation (2.32)), i.e. the time derivative of *K*^{H}_{Ic}, is higher for high hydrogen flux. The change in the gradient of *K*_{I}, because d*K*_{I}/d*a* is several orders lower than ∂*K*^{H}_{Ic}/∂*ξ*. This corresponds to the experimentally observed *K*_{I}-independent stage II crack growth.

Figure 8*b* shows the linear relationship between the cathodic current density, *i*_{c}, and the crack growth rate on a log–log plot. As it has been assumed that HER follows a coupled discharge–recombination mechanism, the steady-state permeation current density is proportional to the square root of the charging current density, as expressed in equation (3.17). Thus, the crack growth due to HE is a function of permeation current density

Under crack tip conditions of low pH and low potential, the cathodic current density is limited by diffusion in double layers, and the limiting current density is a function of pH,
*et al.* [27]. Thus, the total crack growth rate may be represented by a simple superposition law as

## 4. Model application to predict stress corrosion cracking growth rate

The model developed in §2 is now applied to a realistic surface-breaking crack. The crack geometry is almost the same as in figure 2*a* except for the crack length, which is *a*=1.7 mm, and the crack mouth opening displacement, which is *h*_{m}=17.2 μm. It is possible that this crack length will go undetected during manufacturing. Normally, the length of a crack detectable during manufacturing is longer than 1.9 mm using penetrant inspection only. It is assumed that the crack is exposed to an aqueous solution containing 0.1 M NaCl. The external loading is taken to result in a stress intensity factor of

Figure 3*a* shows the changes in pH, [Cl^{−}] at the crack tip at *E*_{App}=−0.860 V_{SCE} versus time. As there is no buffering species, the pH drops quickly as soon as the [Al^{3+}] is introduced to the solution. The precipitation of gibbsite is also important in pH control. Without precipitation, a higher pH ∼5 is calculated with the same conditions. In experiments, so-called ‘mud-crack’, which is caused by the precipitation of solute species, is observed [2]. The low pH and cathodic potential are favourable for the production of hydrogen bubbles, thus the fraction of hydrogen bubbles at the crack tip soon reaches the saturation point, *f*_{Max}. Thus, the effect of the maximum volume fraction of hydrogen bubbles at the crack tip is not as significant as in §3. The concentration of chloride is stable after reaching its highest point. The current density at the crack tip is shown in figure 9*b*. At more anodic potential (*E*_{App}=−0.495 V_{SCE}), the simulation yields a lower pH with higher potential at the crack tip, which causes higher anodic current density, but the cathodic current density is relatively constant. As a result, a faster crack growth rate has been calculated. At more cathodic potential (*E*_{App}=−1.200 V_{SCE}), the crack tip cannot reach critical conditions.

As shown in table 2, the crack growth rate varies with the applied potential. The steady-state growth rate is 1.83×10^{−8} mm s^{−1} (*E*_{App}=−1.200 V_{SCE}), and at a higher anodic applied potential (*E*_{App}=−0.860 V_{SCE}), the crack growth rate is 5.58×10^{−6} mm s^{−1}. If the crack is continuously exposed to electrolyte solution, the crack will reach an inspectable crack length (2.5 mm) in 240 days (*E*_{App}=−1.200 V_{SCE}). This result agrees with the data from Nguyen *et al.* [2], who have also observed a faster crack growth at more anodic applied potentials. Their observation of the deposition of corrosion product, as well as hydrogen bubbles, indicates that the gaseous and solid phases may have affected the transport of species.

## 5. Conclusion

A comprehensive SCC growth rate model has been developed. The model includes both AD and HE. For AD, the modified form of the Nernst–Planck equation with the effect of hydrogen bubbles and precipitates has been applied to simulate the transport of moderately dilute species. The coupling of PHREEQC through the sequential non-iterative approach makes the reactive transport simulation easier. For HE, a simple but elegant form of the crack growth rate equation has been derived, and the stress-assisted hydrogen diffusion with effective diffusivity has been used to obtain the distribution of hydrogen atoms which lowers cohesive energy. The most important finding from this work is that both AD and HE are important for crack growth rate prediction, and affect each other in the process of SCC.

The model shows good agreement with experimental data, as discussed in §3. A two-stage behaviour of potential drop has been observed. The crack tip pH is much lower than in bulk solution. Hydrogen bubbles readily form and accumulate within the crack during localized corrosion. These crack tip conditions promote hydrogen diffusion in the region in front of the crack tip, therefore both AD and HE contribute to the crack growth rate. The crack growth rate due to the combined effects of AD and HE has been determined, and numerical results for comparison with experimental data, as well as a calculation of the crack growth rate for a practical configuration, are presented.

The result also shows that the distribution of phases that impede the transport of electrolyte species is important for crack growth rate prediction. However, the calculation of the distribution of phases is a difficult task because of its probabilistic nature. The tortuosity of SCC makes it more complicated. Therefore, it is necessary to analyse this with a probabilistic technique, and this will be future work.

## Data accessibility

There are no data to report in this manuscript

## Authors' contributions

J.D.A. proposed the problem and supervised the research. D.L. formulated and worked out the details, and wrote the first draft. Y.H. supervised the research.

## Competing interests

All authors have declared that no competing interests exist

## Acknowledgements

The authors thank Professor Yun-Jae Kim at Korea University for providing the Abaqus user subroutines. They also wish to thank the anonymous reviewers for their constructive comments.

- Received September 17, 2014.
- Accepted April 17, 2015.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.