A Peridynamic Model for Crevice Corrosion Damage

A new peridynamic (PD) model for crevice corrosion damage is introduced and cross-validated with experimental results reported in the literature. The model defines a simple metal ion concentration-dependent corrosion rate along the metal-electrolyte interface. Crevice problems have ratios of gap-to-length of 1/100 or lower. To increase the computational efficiency when having a domain with an extreme aspect ratio, the PD formulation with spherical horizons is modified to accommodate arbitrary-shape horizons. The model is validated against experimental results on immersed bolted washers. The PD simulations predict the observed corrosion kinetics, and the deepest corrosion trenches form at a critical distance from the mouth of the crevice. The location of the most severe corrosion attack does not have to be specified as an input. Instead, it results directly from solving the problem with PD. The evolution of crevice corrosion damage is autonomous in PD, and only the geometry, initial, and boundary conditions control the evolution of the corrosion process.


Introduction
Crevice corrosion is a type of localized corrosion that occurs in locations where the metallic surface is exposed to a confined, stagnant electrolyte in a "crevice," while the rest of the metallic surface is in contact with the bulk electrolyte [1].Restricted flow in the crevice slows down the transport of chemical species and leads to local acidification (pH drop), triggering a self-accelerating anodic dissolution of the metal surface in the crevice [1].Crevice corrosion damage is considered a significant problem in many industries.Joints, fasteners, and most types of contacts in ships, aircrafts, infrastructures, or any other structures in offshore and marine environments, are highly susceptible to crevice corrosion attack [2,3].For example, in bridges' tendons, crevice corrosion occurs between strands and the grout surrounding them and even between twisted wires in a strand, and may contribute to catastrophic failure of bridges [4].
Computational models for corrosion damage phenomena, if predictive, are of significant interest since they provide a tool to virtually investigate the potential damage caused by environmental factors [5].Computational models output chemical speciation, the evolution of various electrolyte properties such as pH, conductivity, and, more importantly, potential and current density profiles, which ultimately can determine the corrosion rate.Models for crevice corrosion are of two classes: 1) the first class uses a stationary domain, where the governing equations (usually mass transport and electrostatics) are solved within a fixed time-invariant domain (of the gap only) [6][7][8][9][10][11][12][13][14][15][16][17][18].The second class of crevice corrosion models considers evolving geometry, where the anodic dissolution changes the shape of the crevice in time as the corrosion progresses [19][20][21][22][23].The first class of models is computationally more efficient because they do not deal with an evolving domain.However, they have obvious limitations because, at best, they can only provide rough estimates of the corrosion profile based on the computed current densities over the original crevice domain (see for example [24,25]).Since geometrical changes of crevices influence the transport phenomenon and the electric potential distribution on the anodic surface, the second class is more realistic in simulating the corrosion damage front profile.This profile is important because, under mechanical loadings, cracks can initiate from the deep trenches carved by crevice corrosion attack [4].
Traditionally, corrosion problems with evolving domains have been described by Partial Differential Equations (PDEs)-based models.Some of such models have regarded corrosion as a moving boundary problem.For example, the Finite Element Method has been used to solve the PDEs in the domain configuration at each time step and the level set [21] or Arbitrary Lagrangian-Eulerian (ALE) [20] methods are used to update the domain boundary, leading to a new domain configuration for the next time step.Models with moving boundaries face serious challenges in their discretization as they need to change the domain and meshes at each time step.More discussion on the limitation of these models are provided in [5,26].Some other PDE-based models consider corrosion propagation as a moving interface problem in a domain consisting of both the liquid phase (electrolyte) and a solid phase (metal) (see, e.g., [19]).However, because of continuous changes on the solid-liquid interface, it is hard to predict those changes, forcing such PDE-based models into enforcing boundary conditions on the moving interface that do not always reflect reality.Combinations of numerical methods and sometimes ad-hoc techniques have been adopted to simulate moving boundaries/interfaces [5].For example, in one study, the Finite Volume Method (FVM) was used for PDEs inside the domain, and a Voxel method was used to solve the moving interface problem [19].
More recently, a new class of models has emerged that solves the governing equations defined over a two-phase, electrolyte-solid domain, and predict the evolution of the corrosion front caused by anodic dissolution more efficiently This approach eliminates the need to explicitly track the corrosion front, simplifies numerical complexity, and improves a model's applicability to more complex situations.One such class of models is Cellular Automata (CA) [27,28].However, given their discrete nature and the heuristic rules for "cell" transformation they implement, they are difficult to calibrate and less applicable for quantitative predictions [5].While CAs may offer results that replicate certain qualitative aspects of an observed phenomenon, validations against experimental results are almost non-existent.Peridynamic (PD) models of corrosion [29][30][31][32] are another class of models that does not require tracking the moving boundary explicitly.PD is a nonlocal approach that replaces spatial derivatives with integrals in its formulations.This change allows PD models to naturally capture autonomous emergence and evolution of discontinuities, moving boundaries, and critical features in modeling corrosion damage.PD models have been shown to be remarkably accurate in modeling fracture [33,34], corrosion-induced fracture [35], and stress corrosion cracking [36,37].
Phase-field corrosion models [22,38] are PDE-based models of corrosion that approximate the material discontinuity at the corrosion front with a smooth transition function over a small length scale, so that spatial derivatives can exist.Phase-field models have shown promise in some corrosion problems, but challenges persists.For example, unrealistically thick cracks/damages develop when simulating corrosion-induced fracture and damage [39,40].
To the best of our knowledge, no crevice corrosion model has yet produced results that have been quantitatively cross-validated with published experimental results in terms of the damage propagation in time.In this study, we introduce a simple and predictive peridynamic model for simulating crevice corrosion damage.The model is based on simplifying the complex phenomenon using a metal ion concentration-dependent parameter.In addition, in order to be able to efficiently handle problems with high aspect ratio geometries, a modified version of the peridynamics formulation (with non-circular horizon regions) is presented, allowing domain discretizations that mimic the given geometry extreme aspect ratio.Crevices with micrometer-sized gaps and centimeter-sized lengths in fasteners are examples of such geometries.The model is validated against published experimental images on the progression of crevice corrosion damage in bolted washers.

Peridynamic corrosion models
In this section, we briefly review the basics of peridynamic (PD) theory, the general formulation of PD corrosion models, and the employed discretization method.

Peridynamics
Peridynamics theory is a nonlocal extension of continuum mechanics [41,42].In this theory, each material point interacts with other material points that are located within its finite size neighborhood.For a point with the position vector , the finite size neighborhood (  ), which is usually taken to be a sphere in 3D (or a disk in 2D, a line segment in 1D) centered at  with the radius  called horizon size.In Section 4 we present a formulation for horizon sizes of arbitrary shapes, including extremely elongated ones, useful in treating problem with extreme aspect ratios.Other points inside   are called the "family" of  and are denoted by the position vector  ̂.Fig. 1, schematically shows a PD body, a generic point , its horizon, and family nodes.Schematic of a peridynamic body , and the nonlocal interactions between a generic material point and its family (from [32]).
The term PD bond refers to objects that carry the nonlocal interactions between two family points.There are different types of bonds, depending on the type of interaction.For example, in a mechanical problem, PD mechanical bonds transmit force densities between points, while for diffusion problems, PD diffusion bonds carry heat/mass between family points.

Peridynamics modeling of corrosion damage
The PD corrosion damage model was originally introduced in [29] and later modified in [32].In this part, we briefly layout the PD corrosion damage formulation based on the modified version in [32].

Formulation of peridynamic corrosion damage
PD corrosion model is based on a damage-dependent nonlocal diffusion equation that governs the mass transport in a two-phase domain consisting of both metal (solid) and electrolyte (liquid) phases.The governing equations are as follows:   Schematics of a PD corrosion domain with focus on the solid-liquid interface and the three different types of PD bonds: solid-solid, liquid-liquid, and solid-liquid (interfacial) bonds [32].
According to Eq.( 2), for the liquid-liquid bonds (bonds connecting two liquid points)  =  diff . diff is called the micro-diffusivity and is a function of the classical diffusivity of the electrolyte () and the horizon size (): This parameter allows modelling of diffusion of metal ions in the electrolyte.Note that we use the italic  for diffusivity and the non-italic D to denote the dimensions.Eq. ( 2) requires solid-solid bonds (bonds connecting two solid points) to have  = 0, implying that no mass transfer in the metal phase is considered.For interfacial bonds that connect solid and liquid points and cross the solid-liquid interface (,  ̂, ) =  diss , where  diss denotes the micro-dissolvability, a parameter that controls the dissolution rate and is calibrated to the current density ().We use a numerical calibration procedure to find the correlation between  diss and .This process is described in Section 2.2.3.
Different corrosion types can be modeled via defining the appropriate formula for current density according to the particular anodic dissolution kinetics of the desired corrosion type.For example, pitting [30,32,43], intergranular [31], stress-dependent [36,37], and galvanic corrosion [26], are respectively modeled by defining  as a function of concentration, various solid phases, deformation, and electric potential.In Section 3, we discuss how to define  diss for crevice corrosion.
The damage index that defines phases in Eq. ( 2), is determined from Eq. ( 3).To use this definition of , one needs to consider the mechanical bonds between solid points.Eq. ( 3) and Eq. ( 4) state that  at a point  is equal to the number of intact mechanical bonds, divided by the total number of mechanical bonds connected to .In order to model the corrosion-induced damage in time, as the anodic dissolution takes place by the mass transfer from solid to liquid, one breaks mechanical bonds accordingly.To this aim, the concentration-dependent damage (CDD) model in Eq. ( 5) is used, to give an expected damage index   proportional to the concentration drop at the solid points near the interface (solid points that are connected to liquid points via interfacial dissolution bonds) [26].In Eq. ( 5),  solid is the molar concentration of pristine metal, equal to the molar mass divided by mass density. sat is the saturation concentration: the maximum possible concentration of metal ions in the electrolyte.Once   is found by Eq. ( 5), one then uses a stochastic procedure to break the mechanical bonds accordingly.This is achieved by randomly assigning (,  ̂, ) = 0 for the bonds of , such that  in Eq. ( 3) is approximately equal to   from Eq. ( 5).One then updates the solid to liquid phase-change by updating  using Eq.(3) and Eq. ( 4).The stochastic bond-breaking procedure was introduced in [29], and is briefly reviewed in Section 2.2.2, where the discretization method is presented.
It is noteworthy to say that this model results in a corrosion front with a -thick graded damaged layer in the solid phase, because of the bonds that have one leg in the solid and the other in the electrolyte.This solid layer is referred to as the dissolving solid as opposed to the intact solid being the rest of the solid phase (see Fig. 2).The intact solid domain does not participate in the computations.The dissolving solid region has also a graded metal concentration between  solid and  sat , and is similar to the partially damaged/dissolved subsurface layer experimentally observed in different corroding metals [44][45][46][47][48].This layer is weaker than the bulk metal (intact solid) and is a potential site for the initiation of cracks under mechanical loading.PD corrosion models naturally capture this corrosion-induced embrittlement [37,46].

Discretization
For the numerical discretization of the model presented in the previous section, we employ a quadraturebased meshfree method [49].We first discretize the domain with a uniform grid (see Fig. 3): Fig. 3. Discretization of space using uniform grid spacing (from [36]).
Let  be the total number of nodes/grids.The PD integral can be approximated by mid-point (one-point Gaussian) quadrature as follows: where   denotes the position vector for node , and Δ  is the volume (area in 2D) of the node   that is covered by the horizon of   .While for most family nodes, Δ  = Δ 2 (in 2D), there are nodes near the horizon edge whose volumes are not fully covered by the horizon of   .We use Eq. ( 8) to approximate Δ  in 2D [50]: We use a first-order ODE solver for integrating in time.In previous studies, explicit Forward Euler has been used for this purpose.However, in explicit schemes, the size of time steps is restricted by stability conditions.In problems were the diffusion in the electrolyte is important, the large diffusion coefficient of the electrolyte restricts the time step to very small values, and therefore, computations for relatively long corrosion times would be very costly with explicit time marching.In this study, we use the implicit backward Euler for time integration which is stable for any time step size [51]: The subscripts  andq refer to the nodal coordinates   and   and the superscripts  and  + 1 refer to the current and next time steps (  and  +1 =  + Δ) respectively.In Eq. ( 9) we use  at the time step   , not  +1 , which means that the phase-change process is explicit, while the transport is solved implicitly.
At each time step, Eq. ( 9) updates the concentration field, and   +1 is computed for each node from Eq.
(5).Then, the following stochastic algorithm is used for breaking mechanical bonds accordingly [29]: After bond breaking,   +1 is updated from: which is the discrete version of Eq. ( 3).  +1 is then used to identify the nodal phase at the next time step.

Numerical calibration of micro-dissolvability
As mentioned in the previous section, micro-dissolvability is numerically calibrated to the anodic current density.The following relationship is assumed between  diss and  [32]: where  trial is the current density obtained from a trial PD simulation assuming an activation-controlled uniform corrosion with  diss =  diss trial , a trial micro-dissolvability. Activation controlled condition here is modeled by setting the (, ) = 0 everywhere in the electrolyte region [29].This eliminates the dependency of the dissolution rate on the transport in the electrolyte and  diss trial becomes the only parameter that controls the corrosion rate.Note that the choice of horizon size and spatial discretization size in the trial simulation should be the same as the ones in the main simulation.For more details on trial simulation please see [29,30].The current density from a trial simulation can be computed using Faraday's second law: where  is the charge number,  is the Faraday's constant,  is the corroding area,  is the corrosion time,  is the total number of nodes in the domain, and Δ  is the nodal volume at   .The numerator in Eq. (12) gives the total mass loss due to anodic dissolution from the initial time until time .Division by  gives the dissolution mass flux which can be translated into current density by multiplying with .
Once Eq.( 11) is used to define  diss in Eq.( 2) in terms of the current density , any type of corrosion can be modeled using the particular formula for the current density specific to the local anodic dissolution kinetics for that particular corrosion type.This approach was used in PD corrosion models for predict pitting [30,32,43], intergranular [31], galvanic [26], and stress-dependent [36,37] corrosion, with the only difference in these models being the different  diss (from different current densities) formulas.In the present study, we will establish a relationship for  diss for crevice corrosion damage, based on the underlying electro-chemo-physics.

Peridynamic crevice corrosion model
In this section, we introduce the new PD model for crevice corrosion damage.To this aim, based on the underlying electro-chemo-physics of crevice corrosion, we construct a formula for  diss that effectively reproduces the kinetics of anodic dissolution inside crevices.
Crevice corrosion is known to be driven by the local environmental changes inside the crevice [1].The steps involved in the dominant mechanism are: 1) restricted flow of electrolyte in the crevice results in accumulation of dissolved positively charged metal ions, produced by anodic dissolution (even with rates as low as the passive current density), or by micro-galvanic corrosion induced by impurities and inclusions on the metal surface inside the crevice; 2) electro-neutrality causes migration of chloride (or other aggressive anions) from the bulk electrolyte into the crevice; 3) as a result, hydrolysis reaction increases, and pH drops; 4) local acidification in the crevice increases the anodic dissolution rate which produces more positively charged ions at a faster rate.These steps are then repeated.Fig. 4(a) shows these four steps in the crevice corrosion mechanism.This self-accelerating dissolution process is restricted by saturation of the electrolyte and salt layer formation that may occur due to the slow mass transfer in the crevice (diffusion path is long and narrow).In locations along the crevice where the electrolyte is saturated (usually near the closed end), anodic dissolution is controlled by the mass transfer and therefore follows a diffusion-controlled regime [1].
The self-accelerating dissolution process and its restriction by saturation of the electrolyte, can be simplified in the following statement: the local corrosion rate increases as the concentration of dissolved metal ion increases, up to saturation of electrolyte.
where  is an increasing function of .
Using Eqs. ( 11) and ( 13), and knowing that the PD corrosion model provides the evolution of metal ion concentration in the electrolyte, we define  diss to be a metal ion concentration-dependent quantity: where  S and  L are respectively the solid and the liquid ends of an interfacial bond.For any specific corrosion system, i.e. metal and environment,  sat and the function () need to be determined. sat is a quantity that can be found in the literature.() however, is a new concept and no standard methods for obtaining it exist.Inspired by [25], the approach that we use in this work to obtain () is the following: 1) Find the relationship between anodic current density and pH at a given potential: (pH), using polarization curves measured at different pH values.2) Substitute the pH in the relationship, with mathematical models and/or empirical equations that calculate the pH in terms of concentration of metal ions: pH().
This two-step process is the same as the one plotted in Fig. 4(b).In the example in Section 5, we show how perform this procedure for a specific case.

Modified PD formulation for using discretization grids with extreme aspect ratios
The extreme aspect ratio of the geometry in crevices (long and narrow) present a significant computational challenge: if the same spacing is used in a uniform domain discretization, the computational cost may be too large; therefore, a discretization that matches the geometry aspect ratio (large spacing along the long direction, and small spacing along the short direction) is desirable.Using such a grid with the standard PD formulation (with spherical horizon) would not work, since we may leave covering nodes only in the dense direction, and no nodes in the coarse direction.Note also that to reduce grid dependency and have an acceptable accuracy in the quadrature used for computing the PD integral, grid spacing should not be larger than ¼ or, at most, 1/3 of the horizon size [52].On the other hand, it is well understood that the horizon size should be smaller than the smallest geometrical feature of the domain to prevent undesired/unrealistic nonlocal effects [53].For crevice corrosion, using the spherical horizon, this means that  has to be several times smaller than the gap size.
To resolve this issue, we introduce the PD formulation for non-spherical horizons (or non-circular in 2D) so that it can work with grids that have extremely different grid densities in different directions.
Note that, for noncircular horizon, (,  ′ , ) in Eq. ( 1) cannot be computed from Eq. ( 2), because that relationship is obtained by a calibration process that assumed a spherical/circular horizon.We propose the following modification to Eq. ( 2) for noncircular horizons:   (,  ̂, ) = {  L  (,   ) , (, ) = 1 and ( ̂, ) = 1 0 , (, ) < 1 and ( ̂, ) < 1  diss  , (, ) = 1 xor ( ̂, ) = 1 (15) where   is the nonlocality range along the bond direction that makes an angle  with the -axis (polar angle).In 3D, the direction is given by the spherical coordinates, and in 1D,   has only two values along the positive and the negative coordinate directions. L (,   ) is computed from Eq. ( 16) which is obtained by replacing  in Eq. ( 6) with   .This makes the micro-diffusivity a direction dependent quantity, in order to match a given direction-independent diffusivity constant D: We now show that the new formulation is consistent with classical isotropic diffusion for the 2D case.Similar proofs can be carried out for the 1D and 3D cases as well.
In [54], the calibration of  in PD diffusion with spherical horizons is carried out by finding  such that the classical flux ( classic ) and the peridynamics flux ( Peri ) are equal for a linear concentration (constant flux) profile.Below, we follow similar steps to show that Eq. ( 16) is a valid calibration of  for PD diffusion equation with arbitrary (non-spherical) horizon.
For homogeneous isotropic diffusion the classical flux is given by: where  denotes the gradient operator, and   and   are the unit vectors in  and  Cartesian directions.
According to [54], PD flux can be defined by: where ℋ  + is the part of the horizon with ( ̂, ) > (, ), and   ̂ is the unit vector along the bond  ̂− .
We consider a linear profile for (, ), and an arbitrary-shape horizon for the PD integral (see Fig. 5).Fig. 5.
A linear concentration profile (a constant flux) and its projections into  and  coordinates; an elliptical PD horizon as an example for demonstration of a non-spherical horizon; and the direction-dependent nonlocality range (  ).
By projecting the linear concentration gradient into the Cartesian coordinates as shown in the Fig. 5, the PD flux can be expressed as: We substitute , with the formula given by Eq.( 16): Given that (, ) profile is linear, one can write: As shown above, using the proposed direction dependent micro-diffusivity in Eq. ( 16) recovers the classical flux for isotropic diffusion in 2D.This formulation works for any horizon shape.The elliptical horizon in Fig. 5, is just an example for demonstration.
Note that the non-spherical horizon in this study does not lead to anisotropic behavior as it may in other studies (e.g., [55]).The reason is that  L  here is calibrated for an isotropic classical model.From a physical point of view, the direction dependent micro-diffusivity in Eq. ( 16) has a lower value in directions where   is larger, and a higher value in directions where   is smaller.This leads to a balanced transport in all directions (at the continuum level), and isotropy is maintained.Now we derive the direction-dependent micro-dissolvability  diss  in Eq. ( 15) for a non-spherical horizon.Assume PD diffusion with a spherical horizon of radius  ref .Given Eqs.( 6) and ( 16 .This implies that we can write: where  diss is a micro-dissolvability calibrated for a spherical horizon with the radius  ref (see Section 2.2.3 for calibration of  diss ).
Note that Eq. ( 6), and consequently Eq. ( 16) and Eq. ( 22), are derived for the particular kernel used in the PD integral in Eq. ( 1).If the PD diffusion equation uses another kernel (e.g.see [52,56]), one needs to modify all derivations accordingly.
For the crevice corrosion simulation presented in the next section, we use an elliptical horizon with the long axis aligned with the crevice length direction, and the short axis along the crevice gap direction.This will allow us to choose coarse grid spacing in the crevice length direction and dense grid spacing in the crevice gap direction.

Model Validation
Here we validate the PD crevice corrosion model described in Section 3 against the experimental results shown in [57].

Brief description of experimental setup
In an experiment reported in [57], two washers were held together using a nut and bolt fastener (see Fig. 6(a)).The washers, bolt, and nut were all made of Nickel alloy 625.The bolted washers were immersed in ASTM artificial ocean water at room temperature, and potentiostatic tests were carried out at 200 mV (vs SCE).The experimental study focused on the crevice corrosion propagation between the washers [57].
Fig. 6.Schematics of the 3D actual geometry (a); a zoom-in for the crevice cross-section in (b); and the 2D domain used in the simulation in (c).

Model input data and the calibration of concentration-dependent 𝒌 𝐝𝐢𝐬𝐬
The input data for our PD model of crevice corrosion are discussed next.Molar concentration of the pristine metal can be approximated by dividing the alloy's mass density over its molar mass:  solid = 140 M. Saturation concentration of the alloy in the electrolyte used here is reported to be  sat = 5.6 M [25].The diffusion coefficient of metal ions in the electrolyte used in the experiment is  = 720 μm 2 /s [57].To find the concentration-dependent  diss formula for this metal-electrolyte system, one first needs to determine the anodic current density as a function of metal ions' concentration.To this aim, we follow the procedure described in Section 3. From published polarization curves carried out at different pH levels [25], we read the current density values associated with different pH values at E=200 mV (SCE).Using Matlab's curve fitting toolbox, we fit an exponential function to these data points.Fig. 7 shows the experimental data points read from [25] and the fitted function.Data points are collected from polarization curves given in [25].
The fitted function gives  (A/cm 2 ) as a function of pH: In the same study [25], pH is provided in terms of metal ions concentration from electrochemical and phenomenological relationships between species.Fitting a function to such data points provides pH as a function of (M) (see Fig. 8): Fig. 8.
As discussed in the procedure explain in Section 2.2.3, we calibrate  diss to  by using a trial simulation of a uniform corrosion that assumes a trial micro-dissolvability.The trial simulation using a spherical horizon with  = 4 μm and uniform grid spacing Δ = Δ = 1 μm with  diss trial = 0.04 μm −1 results in  trial = 2.22 × 10 4 A/m 2 .Using Eq. ( 14) gives: diss ( S ,  L , ) = { (1.8 × 10 −6 ) × 10 1.151 exp{−723[7.405exp(−0.08036( L ,))−6.498]}−3.647( L , ) <  sat 0 ( L , ) ≥  sat (26) Note that in this example we used the pH-dependent polarization curves and the pH-concentration relationships because they were available to us from literature.One can use (or propose) any other experimental/analytical method that provides a reasonable relationship between anodic current density and metal ions concentration to construct the concentration-dependent current density.For example, if data on corrosion rate in terms of chloride concentration is available, one can use the principle of electroneutrality to approximate the corresponding metal ions concentration (at a given chloride concentration) and find ().

Computational model setup
We use a 2D peridynamic model of transverse cross-section of the washers to simulate the crevice corrosion in the system described in [57] (see Fig. 6).The crevice between the washers is measured to be a wedge-shape of length 1.27 cm and with a gap size of 10 μm at the closed end and 50 μm at the mouth.We use the symmetry of the geometry (washers are identical) to define the domain as one-half of the system: one washer and a wedge crevice on the top with the length 1.27 cm and half of the original gap size (5 μm at the end and 25 μm at mouth).Fig. 6 shows how the 2D domain is chosen from the actual 3D geometry, and Fig. 9 shows the initial and boundary conditions used in the simulation.
Fig. 9. Geometry of the 2D computational model, with the initial and boundary conditions used.
For the boundary conditions, we set  = 0 for a -thick layer at the crevice mouth to represent the connection to the bulk dilute electrolyte.The thickness is required because of the special way that nonlocal boundary conditions are defined in PD problems.Details on nonlocal boundary conditions are available in [58,59].The rest of the boundaries (including the symmetry line) are free boundaries and the no-flux conditions is naturally satisfied.
The initial condition  =  solid is imposed over the washer region and  = 0 over the electrolyte region, except for a small area of length 500 μm at the crevice closed end where we impose  = 0.99 sat .The nearly saturated electrolyte near the end causes a local spike in current density according to Eq. ( 26) near the closed end.This end-condition acts as corrosion initiation at the tip of the crevice which could occur from microgalvanic dissolution of or around metallic inclusions on the surface, or due to passive film breakdown caused by other reasons.The length 500 μm is simply selected because it was approximately the minimum amount that could kick off the self-acceleration mechanism.The 500 μm was found by trying several different lengths from 100 and 1000 μm.Using lengths smaller than 500 μm, caused the metal-ion concentration to diffuse out quickly before they could accumulate enough to result in an anodic current density high enough to sustain the self-accelerating cycle illustrated in Fig. 4.
As noticed from the dimensions shown in Fig. 9, the crevice length is two orders of magnitude longer than the gap.We use the PD formulation from Section 4 with an elliptical horizon so that grid spacing along the length can be selected to be 25 times larger than the spacing in the gap direction.We choose the elliptical horizon with  0 = 100 μm and  /2 = 4 μm (subscripts are the values for the polar angle ) along the major and minor axes, respectively, and set grid spacings Δ = 25μm and Δ = 1 μm to discretize the domain.A measure of grid density inside the horizon for PD with spherical horizon and uniform grid spacing is the m-factor defined by /Δ .Note that this choice of horizon and grid spacing results in an m-factor of 4 in both directions:  0 Δ ⁄ =  Note that   in the Backward Euler scheme (Eq.( 9)), is obtained from Eq. ( 26) with  +1 .This leads to a nonlinear system of equations in terms of  +1 .We use a modified Polak-Ribiere nonlinear Conjugate Gradient method [60] to solve the nonlinear system at each time step.The total simulated corrosion time is 72 hours, and the time step is Δ = 5 s.
For the computer simulation, we coded the model into an in-house Fortran 90 program with OpenACC enabled GPU parallel computation.The simulation was performed on a Linux cluster with one Intel Xeon Gold 6248 processor (2.50 GHz, and 27.5 MB Cache) and one Tesla V100 GPU.The simulation took about 30 hours to finish.

Simulation results and discussion
Fig. 10 shows the PD simulation results in several snapshots in comparison with the corresponding experimental observations (at same times and using the same geometrical scales).The experimental graphs (left column) show the corrosion profiles at different times obtained by Optical Profilometry, scanning the depth along a radial line on the washer surface [57].The PD simulation results (right column) show the evolution of damage and of the metal ion concentration in the crevice.The PD snapshots use the same length scales as the ones employed for the experimental results from [52].The window frame shown on the top simulation snapshot represents the corresponding window used for the experimental results shown in the left column.The PD simulation is also provided in Video 1 (see Supplementary Materials).
As observed, the complex evolution of corrosion damage is predicted very well by the PD model.Damage starts at the closed end and moves towards the mouth as time passes, affecting only a superficial layer of material.Progression of the active site toward the mouth stops after about 30 hours.The dissolution then localizes at a critical distance from the mouth, being controlled, autonomously, by the diffusion conditions near the crevice mouth, where dilute electrolyte enforces conditions that are wellapproximated by the boundary condition we imposed in the model at that end (see Fig. 9).This stagnation of the active site leads to deep carving into the washer near the crevice mouth.This observation is consistent with other studies on crevice corrosion in Nickel Alloy 625 washers [61][62][63].
The PD model presented here helps us explain the underlying mechanism in crevice corrosion.According to the simulation results, the accumulated ions at the closed end locally increase the current density, which produces more dissolved metal ions; the closed end quickly saturates locally, due to the slow diffusion rate along the almost one-dimensional (narrow and long) path towards the crevice mouth; the current density peak starts to travel to towards the mouth, along with the location in the electrolyte between the saturated and the dilute regions, where the concentration is high enough (to cause a large current) to lead to dissolution but smaller than  sat ; as the solution saturates, the solid does not passivate but saturation induces diffusion controlled corrosion with a significantly lower diffusion/dissolution rate; corrosion damage slows down in the vertical direction, into the washer; as the current peak keeps moves towards the mouth, it reaches a location at a critical distance from the mouth where the diffusion rate is high enough to prevent further saturation (as the shorter distance to the mouth results in a higher diffusion flux), and the dissolution rate and diffusion rate become balanced; when this process reaches semi-steady state transport, the peak current density location with high  <  sat stops from translating to the left, and stabilizes , causing a deep attack in that particular region.Note that all this complex behavior is obtained Fig. 10.
Comparison of experimental results (left column, from [57]; note that each panel is from distinct washers coming from distinct experiments) and PD simulation results for crevice corrosion.The time and length scales for simulation and experiments are identical.The colors in the metal region (see Fig. 9) show the evolution of the nodal damage index, while those in the electrolyte region (see Fig. 9) indicate the metal ion concentration in the crevice.
autonomously by the PD model that only uses a simple concentration-dependent  diss .The PD crevice corrosion model introduced here is the first computational model to validate experimental results on crevice corrosion damage evolution with such details.
As seen from the results in Fig. 10, the deep trenches carved in the crevice corrosion process can serve as initiation points for cracks in SCC.Given the easy and accuracy with which PD can model fracture, a uniquely valuable advantage of PD corrosion models is that one can now easily simulate SCC by simply coupling the corrosion damage model to a PD fracture model as done in, for example, [37], where pit-tocrack transition in a turbine steel was accurately predicted using a 3D PD model.
We note that the model introduced is not only capable of qualitative match with experiments, but also quantitative.As seen from Fig. 10, the evolution of the corrosion damage front and the depth of the local attack found by the PD model are very similar to those from experiments.Please note that the panels showing the experimental data are showing the corrosion front for different washers used in distinct tests.Some small differences between the simulation results and the experiments are in terms of surface roughness at the main damage site.This can be attributed to the stricter boundary condition (fixed location where  = 0 is set) imposed in the model at the crevice mouth than likely exists in reality, and to the presence of microstructural heterogeneities in the metal that were not considered in the current PD model but could be added by using explicit representation of grain and grain boundaries (e.g.[31]).To reduce the likely high computational cost of such a model, an alternative would be to incorporate microstructural influences using the ideas from the intermediate homogenization (IH) PD modeling [64,65].Another small difference between the simulation and the experiment is the location of the deep attack (critical distance).Note that the experiments in [57], in addition to the variability between different washers used in distinct tests, show a considerable degree of variability for this critical distance along different radial directions of the same washer in the same experiment (see Figure .70 in [57]).There are several reasons for this variability, including imperfection on the shape of the washers, the pressure between them, slight asymmetries, etc.On the other hand, the model used here has many simplifications and assumptions, including that a two-dimensional approximation was used to simulate the actual 3D crevice problem.While a 3D PD simulation can be attempted, an axisymmetric PD corrosion formulation could offer similar results at a fraction of the cost.This is planned for the future.

Conclusions
A new peridynamic (PD) model for crevice corrosion damage was introduced and validated against experimental results from the published literature.We simplified the self-accelerating anodic dissolution kinetics in crevice corrosion to a metal-ion concentration dependent current density relationship.This relationship defines the local micro-dissolvability for interfacial PD transport bonds that carry anodic dissolution micro-fluxes.To be able to compute efficiently problems defined over domains with extreme aspect ratios, discretizations with similar aspect ratios are desired, but they were not possible in the standard PD formulation.To solve this problem, we presented a generalized version of the PD formulation that allows horizons with arbitrary shapes, which allows discretizations with highly different grid spacings in different directions.This plays a crucial role in efficiently simulating crevice corrosion, since crevices are often long but very narrow (with differences of two-three orders of magnitude between their length and gap).The model was validated against an experiment from literature on crevice corrosion between two washers of nickel alloy 625.We found the concentration-dependent current density from experimental polarization curves and analytical/empirical relationships.A PD model with the simple concentration-dependent dissolution formulation is able to predict the kinetics of anodic dissolution and damage evolution inside the crevice in great detail.

Fig. 1 .
Fig. 1.Schematic of a peridynamic body , and the nonlocal interactions between a generic material point and its family (from[32]).

( 5 )
In Eq. (1), (, ) denotes the concentration of the dissolving species (here metal atoms/ions) at point  and time .The integrand is the mass flow density which is the molar amount that a unit volume at point  receives from the unit volume at a family point  ̂ in one second.(,  ̂, ) is a constant that is determined from Eq.(2), based on the phases of the points  and  ̂ (solid or liquid) at time .In this model, phases are represented by their damage index: the scalar (, ) ∈ [0,1].If (, ) = 1 at a point  and time , then  is in the liquid phase, and if 0 ≤ (, ) < 1 , then  is the solid phase.

Fig. 2 .
Fig. 2.Schematics of a PD corrosion domain with focus on the solid-liquid interface and the three different types of PD bonds: solid-solid, liquid-liquid, and solid-liquid (interfacial) bonds[32].

Fig. 4 .
Fig. 4. Self-accelerating anodic dissolution mechanism in crevice corrosion: (a) the 4-step cycle, where [M + ] is the concentration of positively charged metal ions, i.e. , and [Cl − ] denotes the concentration of chloride; (b) 3-step equivalent cycle with a concentration dependent pH, and a pH-dependent current density; (c) a 2-step simplified version with a concentration-dependent current density.

Fig. 7 .
Fig. 7. Anodic current density in logarithmic scale in terms of pH value at 200 mV (SCE).
For each node   compute the probability of bond breaking   +1 =