Study on mechanism of differential concentration corrosion

The pipeline easily gets corroded in a seawater environment. The oxygen in the seawater is one of major parameters causing the corrosion. In practice, the corrosion due to the oxygen concentration difference, i.e. differential concentration corrosion (DCC), cannot be avoided. However, a one-dimensional DCC model cannot satisfactorily predict the corrosion because the oxygen distribution near the pipe wall is two-dimensional. In this regard, a two-dimensional DCC model was proposed in this study to numerically investigate the distribution of corrosion potential and current in the ionic conductive layer near the pipe wall as well as the overall corrosion current. The results show that DCC plays a significant role in determining the corrosion potential and current. Without considering DCC, a large corrosion potential and current exist at the location with high oxygen concentration near the pipe wall, whereas the occurrence of the low corrosion potential and current corresponds to the location with low oxygen concentration. However, as DCC is considered, at the location with high corrosion potential, cathodic polarization was produced and the corrosion rate decreases; at the location with low corrosion potential, anodic polarization was produced and the corrosion rate increases. In general, the corrosion potential can be homogenized in terms of DCC.

The differential concentration corrosion (DCC), which is due to the non-uniform distribution of oxygen concentration in an electrolyte 1 , is an important corrosion mechanism. For example, when a structure is immersed in seawater, the corrosion potential at its upper part is higher than that of the lower part because the oxygen concentration near the surface of seawater is higher than that in deep seawater. The potential difference generated between the upper and lower parts of the structure can lead to the transmission of electrons and subsequent corrosion on the structure.
Generally, DCC occurs only when the oxygen concentration difference reaches a certain level. In reality, it is not easy to measure the oxygen concentration throughout the solution. Therefore, DCC was always ignored in the engineering analysis. However, in a complex environment where the distribution of oxygen is highly uneven, the DCC should be considered to help explain the mechanism for some peculiar corrosion phenomena. For example, Matsumura 2 studied the failure of pipelines in Japanese Mihama Nuclear Power Plant. It was found that the outer elbows of these pipelines became thinner and thinner and were even damaged. It is difficult to explain the result by using traditional flow accelerated corrosion theory (FAC). Traditional theory 3 indicates that the bend in the pipeline has a large shear stress due to a fast fluid flow rate leading to a thinner boundary layer and a higher concentration of oxygen. Thus, the oxygen concentration gradient is larger, causing a higher mass transfer rate and higher electrochemical reaction rate. Correspondingly, the corrosion rate increases. Based on this traditional theory, the internal side of pipeline bend should be corroded firstly. However, the fact is contrary to the above estimation, the external elbow of the elbow pipe is first reduced and even damaged. The contradiction result only can be explained by the DCC mechanism.
In the past, this kind of research work [4][5][6][7][8][9][10][11] was carried out under laboratory conditions, and the oxygen concentration in the container was evenly distributed, which could not be applied to engineering. In practice, the distribution of oxygen is more complicated than that in the experimental conditions. For example, the fluid in the seawater pipeline is dissolved with a certain amount of oxygen, significant turbulence tends to cause a complex distribution of oxygen when the fluid flows in the pipeline. The oxygen distribution is very non-uniform even if the shape of a pipeline is not complex. The oxygen distribution can be non-uniform in both axial and circumferential directions. Laboratory conditions cannot meet the needs of engineering practice. In this regard, numerical simulation was applied to study the DCC under complex situations. For example, Lu et al. 12 proposed a model to predict the reducing-pipe flow accelerated corrosion. The concept of DCC was introduced into the traditional FAC model based on the difference in oxygen concentration at each end of a reducing pipe. The corrosion rate www.nature.com/scientificreports/ of the reducing pipe section was then calculated. Zhu et al. 13 predicted the corrosion rate of the loop pipeline elbow in a nuclear power plant, using a DCC model based on oxygen concentration difference between the inner bend and external elbow. There is a very thin ionic conductive layer between the inner bend and external elbow according to the model, forming an electronic conduction loop within the pipeline. An analytical method proposed by Song 14 was adopted to calculate the corrosion rate at several points along the conductive path between inner bend and the external elbow. The result indicates that the corrosion current considering DCC is higher than that without DCC by one order of magnitude. Therefore, the elbow can be corroded and damaged earlier.
The verified result obtained by Masumura 2 indicates that macro-battery corrosion caused by concentration difference must be considered when the distribution of oxygen is not uniform to a certain extent. Previous numerical studies on the corrosion of pipeline have almost exclusively used one-dimensional model assuming a uniform distribution of oxygen. However, such an approximate analysis is not sufficient. In practice, the distribution of oxygen is not uniform. Therefore, considering uneven distribution of the oxygen along the axial and circumferential directions in the pipeline, in this study, a two-dimensional DCC model was proposed to numerically investigate the distribution of corrosion potential and current in the ionic conductive layer near the pipe wall as well as the overall corrosion current. .

Differential concentration corrosion mechanism
The corrosion process can be explained using the following experiment. As illustrated in Fig. 1a, semi-permeable membrane is used for dividing a galvanic cell containing electrolytes into two parts. Two iron sheets with the same properties are submerged into both parts. Different amounts of oxygen are injected into both sides through a particular device to make different oxygen concentrations, i.e. c o 2,2 > c o 2,1 . The following electrochemical reactions can occur at both parts: The equilibrium potential, exchange current density and Tafel slope at both parts, i.e. E e,a , I 0,a , β a and E e,c , I 0,c , β c , are listed in Table 1, respectively. The coupling of reactions (1) and (2) leads to a mixed potential, E corr , and corrosion current, I corr . In general, E corr is far different from the equilibrium potential E e,a , E e,c , so the contraries of the two reactions can be ignored. In addition, the control step of the whole reaction is determined  www.nature.com/scientificreports/ by the discharge process if the concentration of oxygen appropriately increases. Therefore, the above coupling reaction can be described by the simplified Butler-Volmer formula: Consequently, E corr and I corr can be obtained by combining Eqs. (3) and (4) due to |I c | = I a = I corr : The electrochemical reactions on both parts of the galvanic cell can be described using Eqs. (5) and (6) in Fig. 1a. The equilibrium potential, E e,c , of electrode reaction (2) and the exchange current density, I 0,c , vary with different oxygen concentrations on both parts. Based on Nernest equation 15 , the equilibrium potential and exchange current density can be calculated as follows: where c o 2 ,0 is the oxygen concentration at inlet, corresponding to the equilibrium potential E e,c ,and its value is listed in Table 2. Equations (7)- (10) can been substituted into Eqs. (5) and (6), and then the corrosion potentials E 1 corr and E 2 corr as well as the corrosion current I 1 corr and I 2 corr on both parts can be obtained, respectively. This is the corrosion occurred on both parts when sampe 1 and 2 are not connected.
When parts 1 and 2 are connected using wires as illustrated in Fig. 1b, in the condition of E 2 corr > E 1 corr , the electrons flow through a wire from part 1-2 due to potential difference, and the ion directional movement also occurs in the solution. Thus, the external current is generated. The corrosion potential. leads to polarization, i.e. E 2 corr decreases and E 1 corr increases,which help produce different polarization current I 1 P and I 2 P . the change of corrosion current of parts 1 and 2 can be described as I ′ 1 corr = I 1 corr + I 1 P and I ′ 2

Two-dimensional DCC model
Calculation of oxygen distribution. The generative mechanism of the DCC was introduced using the aforementioned experimental devices. But in practice, the corrosion is more complex. The reduced pipeline,depicted in Fig. 2a, is used to help explain the corrosion mechanism in complex conditions. The oxgen distribution near the wall was simulated in ANSYS FLUENT 16.0 before carrying out DCC modelling,a multiphase flow model of mixure was applied. As shown in Fig. 2a, the inlet velocity and oxygen content, c o 2 ,0 , and outlet pressure were set. The values are given in Table 2. No sliding wall condition was applied. The turbulent flow in the pipeline was simulated by κ-ε model. Wall function method was used with the dimensionless distance, y + = 50 , which refers to the distance between the center point of the first layer element and the wall surface. Figure 2b,c show the computational mesh and oxygen distribution near the wall, respectively.
Due to the low density of oxygen, under the action of gravity, most of the oxygen is concentrated in the upper part of the pipeline (Fig. 2c). The results show that the oxygen concentration in the upper part is 3 to 10 times higher than that in the lower part. Therefore, to simplify the simulation condition, the lower part of the pipeline was neglected. Only upper part of the pipeline was used in the subsequent simulation while the www.nature.com/scientificreports/ boundary between the lower and upper part is regarded as the insulation boundary, as illustrated in Fig. 3a. The corresponding near wall boundary layer is given in Fig. 3b.
Numerical model of two-dimensional DCC. Viscous sublayer was used to understand DCC in the pipeline. The schematic of discretization of viscous sublayer is shown in Fig. 4a. Three zones (Zone I, II and III) can be identified. As shown in Fig. 5, a representative element (i, j) and its neighboring elements are taken to show potential and current flow between elements. If DCC is not considered, there only will be a natural corrosion potential E i,j corr and a natural corrosion current I i,j corr in each element, where (i, j) stands for arbitrary element number. But in fact, all elements are interconnected with one another, which inevitably generate current. Thus, the element corrosion potential ( E i,j corr ) can be polarized. The polarization potential is expressed as E i,j − E i,j corr and the polarization current is I i,j F (Fig. 6a,b). E i,j is the ultimate corrosion potential of element (i, j) after polarization. As the elements are connected to each other, the natural corrosion potential difference between them can drive a current, as illustrated in Figs. 4b, 5b, 6a and b. Taking element (i, j) as an example also, as DCC is not considered, the absolute values of anodic reaction current and cathodic current between the wall and element in the solution are equal, and no net current is generated. However, as DCC is considered, the current flow between elements leads to a polarization of the corrosion potential. The anode and cathode currents are no longer equal so that the net current, i.e. external current or polarization current, can be generated. At the same time, other elements around the element (i, j) will also have current flowing in or out, as illustrated in Figs. 4b, 5b, 6a and b. If all elements are considered as a circuit, the element (i, j) can be treated as a node in the circuit. When the current reaches a steady state, according to Kirchhoff 's Second Law 16 , a net flow of current through the element is zero. www.nature.com/scientificreports/ Derivation of discrete equations used in zone I and III. According to the balance of current at a steady state, an equation of the corrosion potential of the polarized element can be obtained. As an example, element (i, j) is surrounded by four neighboring elements. Since each element represents an electrolyte, according to the definition of corrosion potential, E i,j is the potential difference between the wall and the solution and −E i,j represents the potential difference between solution and wall. Therefore, the current flowing from element (i, j − 1) to (i, j) should be Similarly, the current from element (i, j) to element (i, j + 1) is Along the Y direction, the current flowing from element (i + 1, j) to element (i, j) is The current from element (i, j) to element (i − 1, j) is The Faraday current due to polarization is where R i,j P is the polarization resistance of element (i, j), Ω m 2 , more details please see "Discrete equations of elements at the boundary between zones I and II" section. S AA ′ B ′ B ≈ R�θ�x 1 , as illustrated in Fig. 6a. R I x1 ,R I x2 ,R I y2 ,R I y1 is the resistance between elements of the solution in X and Y direction , more details also can be found in "Discrete equations of elements at the boundary between zones I and II" section.
The polarization current is assumed to flow from wall to element. According to Kirchhoff 's Second Law, for each element (i, j), there exist Substitute (11)-(15) into Eqs. (16), (17): www.nature.com/scientificreports/  www.nature.com/scientificreports/ The elements at other boundaries can be treated similarly. Fig. 4b to obtain discrete equations similar to Eq. (17).

Derivation of discrete equations used in zone II. The Kirchhoff Second Law is applied to element (t, s) in Zone II as illustrated in
where R II x1 R II y1 R II x2 R II y2 are the resistances of element of Zone II in X and Y direcion. The detailed expression of them will be found in "Discrete equations of elements at the boundary between zones I and II" section.
It is easy to determine S A ′ EFB ′ according to the geometric relationship in Fig. 6a: where R x is the radius at coordinate x of reducing pipe section, its expression is: and k is r−R l 2 . Other definitions for r, R, l 1 , l 2 are shown in Fig. 3b.  Equation (19) is a standard form. When the element is located at the boundary, the standard form should be changed. For example, for element (m, n) shown in Fig. 4b, because an insulated boundary is above this element, I y in_flow is zero.As Kirchhoff 's Second Law is applied to this element, the equation becomes The elements at other boundaries in Zone II can be treated similarly.

Discrete equations of elements at the boundary between zones I and II.
As the elements (i, k) and (i, k + 1) are at the junction of a straight pipe and a variable-diameter pipe, as illustrated in Figs. 4b and 6c,the calculation of resistance is different from that in Zones I or II . Application of Kirchhoff 's Second Law to element (i, k) can For element (i, k + 1), the equation becomes: where R x12 is the solution resistance between element (i, k) and element (i, k + 1), R x1 , R y1 , R x2 , and R y2 can be expressed thoroughly in "Discrete equations of elements at the boundary between zones I and II" section. Calculation of element resistance in zone II. Because the width of the element in this zone decreases with increasing values of x, the calculation of resistance is different from that in Zones I and III. As illustrated in Fig. 6c, a typical element is highlighted with the red lines. The total resistance should be the integral of differential resistances between neighbouring elements,

Calculation of resistances. Calculation of element resistance in zones I and
Scientific Reports | (2020) 10:19236 | https://doi.org/10.1038/s41598-020-75166-7 www.nature.com/scientificreports/ The integral of Eq. (28) gives where x c 1 and x c 2 is the center coordinates for two adjacent elements. In a similar way, the Y-direction resistance can be calculated by an integral method. The resistance is The integral of Eq. (30) gives where x 1 , x 2 are the element start and end coordinates in X-direction.

Calculation of element resistance at the boundary between zones I and II.
For the elements at the junction of a straight pipe and a variable-diameter pipe, such as elements (i, k) and (i, k + 1), the calculation of resistance between two elements is The similar calculation can be done to obtain element resistance at the junction of zones II and III.
Calculation of plorization resitance. The polarization resistance between the pipe wall and the solution can be calculated according to the following formula: where β a , β c are the Tafel slope of iron and oxgen respectivly listed in Table 1. Figure 3c illustrates the oxygen concentration distribution near the pipeline wall. It can be seen that in the middle part of the upper wall, the concentration of oxygen is higher than that on both sides. In contrast, along the axial direction, near the outlet of the pipeline, the concentration of oxygen is the highest. www.nature.com/scientificreports/ The distribution of oxygen determines the corrosion of the pipeline wall. Figures 7a and 8a show the distribution of natural corrosion potential and current without considering DCC. The distribution of corrosion potential and current is closely in relation to the distribution of oxygen. The high natural corrosion potential and current occurs at location with high oxygen concentration, and vice versa.

Results and discussions
However, because the elements are actually connected, the difference in the natural corrosion potentials of the elements will inevitablely lead to a current flow between the elements. Natural corrosion potential is bound to polarize. Moreover, the absolute values of anodic and cathodic reaction current are no longer equal, causing polarization current. Figures 7b and 8b show the distribution of polarization potential and current, respectively. The final corrosion potential and current after polarization are illustrated in Figs. 7c and 8c, respectively.
The mechanism of concentration corrosion can be described in more details. Figures 9 and 10 illustrate the distribution of the above physical quantities, such as corrosion potential, current et al. at the selected representative rows and columns. Figure 9c shows the polarized potential of each column, indicating the different polarization extents. For example, in the first column, the degree of polarization is the highest and the polarization is anodic. Its means that the corrosion potential increases. Because the natural potential of the element is the lowest, the current mainly flows into the element, resulting in anodic polarization. Similarly, anodic polarization occurs in columns 48 and 99. In comparison to column 1, the higher natural corrosion potential in these two columns causes a lower degree of polarization.   www.nature.com/scientificreports/ For columns 222 and 247, cathodic polarization occurs in the middle part. Because the natural corrosion potential at these location are very high (Figs. 7a, 8b, and 11a). Therefore, the element current is mostly in the outflow state, and cathodic polarization is dominant.
The polarization of the elements in these columns has something in common, i.e. the polarization degree near the circumferential edges is higher than that at the middle part. This is because the concentration of oxygen is lower at the edge and the corresponding corrosion potential is also lower.
Driven by the polarization potential, the polarization current is generated, and its distribution is illustrated in Figs. 10b and 11b. In general, anodic polarization causes anodic current while cathodic polarization causes cathodic current and corrosion. The corrosion current is the algebraic sum of natural corrosion and polarization current, as illustrated in Figs. 8c and 10c. It can be seen that the polarization current is basically in the same order of magnitude as the natural corrosion current. Thus, the polarization current presents a significant influence on the final corrosion current distribution. It indicates that the concentration corrosion cannot be ignored in the corrosion analysis. Due to the concentration corrosion, at the location with high original corrosion current, the corrosion current tends to decrease, whereas the corrosion current can increase at the location with low original corrosion current. Additionally, if the solution resistance is not considered, the potentials of all elements will eventually tend to be uniform.

Conclusions
A two-dimensional DCC model was developed to predict the distribution of corrosion potential and current in the pipeline. The calculation results present a significant influence of concentration corrosion on the overall corrosion of seawater pipeline. The existence of concentration corrosion helps polarize the corrosion potential and subsequently causes the polarization current. The high natural corrosion potential area can cause cathodic polarization and cathodic current, leading to an increase in the corrosion rate. In contrast, the original low natural corrosion potential area can cause anodic polarization and anodic current so that the corrosion rate decreases. The corrosion potential tends to be homogenized due to the differential concentration corrosion. All these findings help clarify the corrosion mechanism in the seawater pipeline with the existence of the differential concentration.

Methods
Geometric modeling. Taking pipeline fluid as the research object. The geometry of the fluid is shown in Fig. 2b. The geometry is generated in software COMSOL and exported in format *.x_b.

Structured grid generation.
In order to esTablelish the DCC model successfully, the fluid must be divided into structured grids. This is done in ICEM,a model of Ansys software.The geometry of the fluid(in format *.x_b) was imported into ICEM and structured grid was generated using mapping method. At the same time, the boundary conditions are defined in this model. After meshing, it is saved into MSH format and imported into FLUENT.
Calculation and data export. In FLUENT, select the type of calculating model (multiphase flow of mixture, turbulence model κ-ε), select the fluid material (the primary phase is water, the secondary phase is oxygen), and then set various boundary conditionss (inlet velocity, oxygen concentration, etc.), and then calculate. After the calculation, UDF technology is used to output the element center coordinates and oxygen concentration close to the wall surface to the file for matlab processing.
Data processing and display in MATLAB. Use MATLAB script program to read the above data. Then calculate in the following order: