A new peridynamic (PD) model for crevice corrosion damage is introduc ed 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 bo lted 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 inp ut. 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. Keywords: crevice corr osion; peridynamics; corrosion; modeling; nonlocal mass transfer; nickel alloy 1. 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 l ocal acidification (pH drop), triggering a self -accelerating anodic dissolution of the metal surface in the crevice [1]. Crevice corrosion damage is considered a significant pr oblem in many industries. Joints, fasteners, and most type s 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 bridg es [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 corrosio n 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-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-23]. The first class of models is computationally more efficient *Corresponding Author Email Address : fbobaru2@unl.edu 2 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] ). Sin ce 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, und er 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 do main 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 d iscussion 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 liq uid 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 nume rical 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, an d 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 elimi nates 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 discre te 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 ce rtain qualitative aspects of an observed phenomenon , validations against experimental results are almost non -existent. Peridynamic (PD) models of corrosion [29-32] are another class of models that do es not require track ing the moving boundary explicitl y. PD is a nonlocal approach that replace s spatial derivatives with inte grals 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 i n 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 c orrosion 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 persist s. For example, unrealistically thick cracks/damage s develop when simulating corrosion -induced fracture and damage [39, 40] . To the best of our knowledge, no crevi ce 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 c revice 3 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. 2. Peridynamic corrosion models In this section, we briefly review the basi cs of peridynamic (PD) theory, the general formulation of PD corrosion models, and the employed discretization method. 2.1. 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 fi nite 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 deno ted by the position vector 𝒙̂. Fig. 1 , schematically shows a PD body, a generic point 𝒙, its horizon , and family nodes. Fig. 1. 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. 2.2. Peridynamics modeling of corrosion damage The PD corrosion damage model was originally introduced in [29] and la ter modified in [32]. In this part, we briefly layout the PD corrosion damage formulation based on the modified version in [32]. 4 2.2.1. 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 cons isting of both metal (solid) and electrolyte (liquid) phases. The governing equations are as follows: 𝜕𝐶(𝒙,𝑡) 𝜕𝑡=∫𝑘(𝒙,𝒙̂,𝑡)𝐶(𝒙̂,𝑡)−𝐶(𝒙,𝑡) |𝒙̂−𝒙|2𝐻𝑥d𝑉𝑥̂ (1) 𝑘(𝒙,𝒙̂,𝑡)= {𝑘diff(𝐷,𝛿),if 𝑑(𝒙)=1 and 𝑑(𝒙̂)=1 0 ,if 𝑑(𝒙)<1 and 𝑑(𝒙̂)<1 𝑘diss(𝑖) ,if 𝑑(𝒙)=1 xor 𝑑(𝒙̂)=1 (2) 𝑑(𝒙,𝑡)=1−∫𝜇(𝒙,𝒙̂,𝑡)𝐻𝑥d𝑉𝑥̂ ∫d𝑉𝑥̂𝐻𝑥 (3) 𝜇(𝒙,𝒙̂,𝑡)= {1,if there is an intact mechanical bond between 𝒙 and 𝒙̂ at time 𝑡 0,if there is no mechanical bond between 𝒙 and 𝒙̂ at time 𝑡 (4) 𝑑𝑐(𝒙,𝑡)={1 ,𝐶(𝒙,𝑡)≤𝐶sat 𝐶solid−𝐶(𝒙,𝑡) 𝐶solid−𝐶sat ,𝐶sat<𝐶(𝒙,𝑡)≤𝐶solid (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 th e 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. 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 ( 𝛿): 5 𝑘diff(𝐷,𝛿)= { 𝐷 𝛿 ,for 1D 4𝐷 𝜋𝛿2 ,for 2D 9𝐷 𝜋𝛿3 ,for 3D (6) 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 numerica l calibration procedure to find the correlation between 𝑘diss and 𝑖. This process is described in Section 2.2.3 . Different corrosion types can be m odeled 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 b y 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 bo nds, 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 mec hanical 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 t hat 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-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] . 6 2.2.2. Discretization For the numerical discretization of the model presented in the previous section, we emp loy a quadrature - based 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: d𝐶 d𝑡(𝒙𝑝,𝑡)=∑𝑘(𝒙𝑝,𝒙𝑞,𝑡)𝐶(𝒙𝑞,𝑡)−𝐶(𝒙𝑝,𝑡) |𝒙𝑞−𝒙𝑝|2Δ𝑉𝑝𝑞 𝒙𝑗∈𝐻𝒙𝑖; 𝑝=1,2,…,𝑁 (7) 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]: Δ𝑉𝑝𝑞= { Δ𝑥2 ,if |𝒙𝑞−𝒙𝑝|≤(𝛿−Δ𝑥 2) [(𝛿+Δ𝑥 2)−|𝒙𝑞−𝒙𝑝| Δ𝑥]Δ𝑥2,if (𝛿−Δ𝑥 2)<|𝒙𝑞−𝒙𝑝|≤(𝛿+Δ𝑥 2) 0 ,if (𝛿+Δ𝑥 2)<|𝒙𝑞−𝒙𝑝| (8) 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 u se the implicit backward Euler for time integration which is stable for any time step size [51]: 𝐶𝑝𝑛+1−𝐶𝑝𝑛 Δ𝑡=∑𝑘𝑝𝑞(𝑑𝑝𝑛,𝑑𝑞𝑛)𝐶𝑞𝑛+1−𝐶𝑝𝑛+1 |𝒙𝑞−𝒙𝑝|2Δ𝑉𝑝𝑞 𝒙𝑞∈𝐻𝒙𝑝; 𝑝=1,2,…,𝑁 (9) 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 7 𝑡𝑛, 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]: ▪ For each node 𝒙𝑖 compute the probability of bond breaking 𝑃𝑝𝑛+1=1 1−𝑑𝑐𝑛+1(𝐶𝑝𝑛−𝐶𝑝𝑛+1 𝐶solid−𝐶sat); ▪ For each intact mechanical bond connected to 𝒙𝑖, generate a random number from a uniform distribution in (0,1); - If this random number is smaller than 𝑃𝑝𝑛+1, then break this mechanical bond: 𝜇𝑝𝑞𝑛+1=0; After bond breaking, 𝑑𝑖𝑛+1 is updated from: 𝑑𝑝𝑛+1=1−∑𝜇𝑝𝑞𝑛+1𝒙𝑞∈𝐻𝒙𝑝 ∑1 𝒙𝑞∈𝐻𝒙𝑝; 𝑝=1,2,…,𝑁 (10) which is the discrete version of Eq. (3). 𝑑𝑝𝑛+1 is then used to identify the nodal phase at the next time step. 2.2.3. Numerical calibration of micro -dissolvability As mentioned in the previous section, micro -dissolvability is numerica lly calibrated to the anodic current density. The following relationship is assumed between 𝑘diss and 𝑖 [32]: 𝑘diss= (𝑘disstrial 𝑖trial)𝑖 (11) where 𝑖trial is the current density obtained from a trial PD simulation assuming an activation -controlled uniform corrosion with 𝑘diss=𝑘disstrial, a trial micro -dissolvability. Activation controlled condition here is modeled by setting the 𝐶(𝒙,𝑡)=0 everywhere in th e electrolyte region [29]. This eliminates the dependency of the dissolution rate on the transport in the electrolyte and 𝑘disstrial 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: 𝑖trial= 𝑧𝐹∑[𝐶(𝒙𝑗,0)−𝐶(𝒙𝑗,𝑡)]𝑁 𝑗=1 Δ𝑉𝑗 𝐴𝑡 (12) where 𝑧 is the charge number, 𝐹 is the Faraday’s constant, 𝐴 is the corroding area, 𝑡 is the corrosion t ime, 𝑁 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 multip lying 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 8 present study, we will establish a relationship for 𝑘diss for crevice corrosion damage, based on the underlying electro -chemo -physics. 3. 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 insid e 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 crevic e; 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 transf er 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 . Fig. 4(c) shows this simplified interpretation. Consequently, one can describe the local anodic current den sity 𝑖(𝒙,𝑡) as a function of the local metal ion concentration: 𝑖(𝒙,𝑡)= {𝑓[𝐶(𝒙,𝑡)]𝐶(𝒙,𝑡)<𝐶sat 0𝐶(𝒙,𝑡)≥𝐶sat (13) 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 : 𝑘diss(𝒙S,𝒙L,𝑡)= {(𝑘disstrial 𝑖trial)𝑓[𝐶(𝒙L,𝑡)] 𝐶(𝒙L,𝑡)<𝐶sat 0 𝐶(𝒙L,𝑡)≥𝐶sat (14) 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(𝐶). 9 This two -step process is the same as the one plotted in Fig. 4(b). In th e example in Section 5, we show how perform this procedure for a specific case. 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 concen tration dependent pH, and a pH -dependent current density; (c) a 2 -step simplified version with a concentration -dependent current density. 4. 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 no des 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 no nlocal 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: 10 𝑘𝜃(𝒙,𝒙̂,𝑡)= {𝑘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 directi ons. 𝑘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: 𝑘L(𝐷,𝛿𝜃)= { 𝐷 𝛿𝜃 ,1D 4𝐷 𝜋𝛿𝜃2 ,2D 9𝐷 𝜋𝛿𝜃3 ,3D (16) 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 (const ant 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: 𝒒classic= −𝐷𝛁𝐶=−𝐷(𝜕𝐶 𝜕𝑥𝒆𝒊+𝜕𝐶 𝜕𝑦𝒆𝒋) (17) where 𝛁 denotes the gradient operator, and 𝒆𝒊 and 𝒆𝒋 are the unit vectors in 𝑥 and 𝑦 Cartesian directions. According to [54], PD flux can be defined by: 𝒒peri=−∫𝑘𝐶(𝒙̂,𝑡)−𝐶(𝒙,𝑡) |𝒙̂−𝒙|𝒆𝒙̂ d𝑉𝒙̂ ℋ𝑥+ (18) 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 ). 11 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: 𝒒peri=−𝒆𝒊∫∫𝑘𝐶(𝒙̂,𝑡)−𝐶(𝒙,𝑡) |𝑥̂−𝑥|cos2𝜃 𝑟d𝑟d𝜃𝛿𝜃 0𝜋 2 −𝜋 2 −𝒆𝒋∫∫𝑘𝐶(𝒙̂,𝑡)−𝐶(𝒙,𝑡) |𝑦̂−𝑦|sin2𝜃 𝑟d𝑟d𝜃𝛿𝜃 0𝜋 0 (19) We substitute 𝑘, with the formula given by Eq. (16): 𝒒peri=−𝒆𝒊∫∫4𝐷 𝜋𝛿𝜃2𝐶(𝒙̂,𝑡)−𝐶(𝒙,𝑡) |𝑥̂−𝑥|cos2𝜃 𝑟d𝑟d𝜃𝛿𝜃 0𝜋 2 −𝜋 2 −𝒆𝒋∫∫4𝐷 𝜋𝛿𝜃2𝐶(𝒙̂,𝑡)−𝐶(𝒙,𝑡) |𝑦̂−𝑦|sin2𝜃 𝑟d𝑟d𝜃𝛿𝜃 0𝜋 0 (20) Given that 𝐶(𝒙,𝑡) profile is linear, one can write: 𝒒peri=−𝐷𝜕𝐶 𝜕𝑥𝒆𝒊∫∫4 𝜋𝛿𝜃2cos2𝜃 𝑟d𝑟d𝜃𝛿𝜃 0𝜋 2 −𝜋 2−𝐷𝜕𝐶 𝜕𝑦𝒆𝒋∫∫4 𝜋𝛿𝜃2sin2𝜃 𝑟d𝑟d𝜃𝛿𝜃 0𝜋 0 =−𝐷(𝜕𝐶 𝜕𝑥𝒆𝒊+𝜕𝐶 𝜕𝑦𝒆𝒋)=𝒒classic (21) 12 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 deri ve 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), one c an write 𝑘L𝜃= 𝑘L(𝛿ref 𝛿𝜃)𝑛 where 𝑛=1,2,3 for 1D, 2D ,and 3D , respectively . This means that one can use a 𝑘 that is calibrated for the spherical horizon with 𝛿=𝛿ref, in a model with non -spherical horizon, by multiplying the flow density in Eq. (1) by (𝛿ref 𝛿𝜃)𝑛 . This implies that we can write: 𝑘diss𝜃= 𝑘diss(𝛿ref 𝛿𝜃)𝑛 (22) 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. 5. Model Validation Here we validate the PD crevice corrosion model described in Section 3 against the experimental results shown in [57]. 5.1. 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]. 13 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). 5.2. 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 allo y’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 μm2/s [57]. To find the concentration -depen dent 𝑘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. Fig. 7. Anodic current density in logarithmic scale in terms of pH value at 200 mV (SCE). Data points are collected from polarization curves given in [25]. The fitted function gives 𝑖 (A/cm2) as a function of pH: log(𝑖)=1.151exp(−723pH)−7.647 (23) 14 In the same study [25], pH is provided in terms of metal ions concentration from electrochemical and pheno menological relationships between species. Fitting a function to such data points provides pH as a function of 𝐶(M) (see Fig. 8 ): Fig. 8. pH in terms of dissolved metal concentration from [25], and the fitted function pH=7.405exp(−0.08036 𝐶)−6.498 (24) Where 𝐶 is in mol/m3.Substituting pH in Eq. (23) with the function in Eq. (24), we find: 𝑖(𝐶)=101.151exp{−723[7.405exp(−0.08036 𝐶)−6.498]}−3.647 (25) where 𝑖 is in A/m2. 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 𝑘disstrial= 0.04 μm−1 results in 𝑖trial=2.22×104 A/m2. Using Eq. (14) gives: 𝑘diss(𝒙S,𝒙L,𝑡) = {(1.8×10−6)×101.151exp{−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 i n 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 𝑖(𝐶). 15 5.3. Computational model setup We use a 2D peridynamic m odel 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 nearl y 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 co rrosion 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 length s smaller than 500 μm, cause d the metal -ion concentration to diffuse out quickly b efore they c ould 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 directi on. 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 de nsity inside the horizon for PD with spherical horizon and 16 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Δ𝑥⁄=𝛿𝜋 2Δ𝑦⁄=4. 𝛿ref and Δ𝑥 in the calibration simulation are 4 μm and 1 μm respectively. 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 simulat ed 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 . 5.4. 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 prof iles 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 employe d for the experimental results from [ 52]. The wind ow frame shown on the top simulation snapshot represent s 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 well - approximated 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 wa shers [61-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 dens ity 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 p articular region. Note that all this complex behavior is obtained 17 Fig. 10. Comparison of experimental results (left column, from [57]; note that each panel is from distinct washers coming from distinct experiment s) 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. 18 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 corrosi on 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 -to- crack 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 experim ental 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 stri cter 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 cou ld be added by using explicit representation of grain and grain boundaries (e.g. [31]). To reduce the likely high computational cost of such a mode l, 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 dir ections 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 fo r the future. 6. 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 corrosi on 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 defin ed 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 horiz ons 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 t wo-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 nick el alloy 625. We found the concentration -dependent current density from experimental p olarization 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. 19 7. Acknowled gements The work of S.J., J.Z. and F.B. was supported by the National Science Foundation under CMMI award No. 1953346, and by a Nebraska System Science award from the Nebraska Research Initiative. The computations have been performed utilizing the Holland Computing Center of the University of Nebraska, which receives support from the Nebraska Research Initiative. 8.