- Email: [email protected]
Mathematical model for cathodic protection in a steel-saline water system
Accepted Manuscript Title: Mathematical Model for Cathodic Protection in a Steel Saline Water System Author: Khalid W. Hameed Aprael S. Yaro Anees A. Khadom PII: DOI: Reference:
S1658-3655(15)00076-X http://dx.doi.org/doi:10.1016/j.jtusci.2015.04.002 JTUSCI 185
To appear in: Received date: Revised date: Accepted date:
8-1-2015 30-3-2015 13-4-2015
Please cite this article as: K.W. Hameed, A.S. Yaro, A.A. Khadom, Mathematical Model for Cathodic Protection in a Steel - Saline Water System, Journal of Taibah University for Science (2015), http://dx.doi.org/10.1016/j.jtusci.2015.04.002 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Mathematical Model for Cathodic Protection in a Steel - Saline Water System
Khalid W. Hameed 1, Aprael S. Yaro2 and Anees A. Khadom*,3 Biochemical Engineering Department, Al – Khawarzim Engineering College, University of Baghdad, Aljadrea
ip t
1
71001, Baghdad, Iraq 2
Chemical Engineering Department, College of Engineering, University of Baghdad, Aljadrea 71001, Baghdad,
us
*
Chemical Engineering Department, College of Engineering, University of Daiyla, Baquba 32001, Daiyla, Iraq
Correspondence author: [email protected], mobile number: 00964 790 2305786
an
3
cr
Iraq
Abstract
The corrosion of steel tubes in sea water was controlled by cathodic protection. The impressed current technique
M
was used. The rate of reaction was evaluated as a function of the temperature, pH and solution velocity. In this technique, the polarization method was used to determine the protection potential and current. The rate of zinc
d
consumption, the protection potential, and the protection current are highly dependent on the variables of the study.
te
The boundary element technique was suitable for modelling corrosion problems. The average percentage of error among the experimental and theoretical data was 1.27%.
Ac ce p
Keywords: Corrosion; Electrochemical method; Metals; Polarization; Cathodic protection 1. Introduction:
Among the various corrosion control methods available, cathodic protection is commonly adopted to control the corrosion of steel. The cathodic protection system is aimed to shift the potential of steel to the least probable range for corrosion [1]. The first attempts at cathodic protection modelling were made using finite difference approaches, but this technique is limited to two-dimensional and axisymmetric problems. Recently, the finite element method has been used successfully by some researchers [2]. Unlike other applications of the finite element method, the object of a cathodic protection system is not the structure itself but the environment close to the structure. Therefore, it is required to construct a finite element mesh of the medium between the members of the structure and to extend the mesh along the distance to establish a realistic boundary. The boundary element method requires only the discretization of the anode and cathode surfaces; therefore, the numerical problem may be reduced 1
Page 1 of 13
in size, which permits better resolution and a reduction in computer time compared to other methods, particularly for complex geometries [3]. Cathodic protection involves the application of a direct current (DC) from an anode through the electrolyte to the surface to be protected. This is often thought of as "overcoming" the corrosion currents that
ip t
exist on the structure. Cathodic protection eliminates the potential differences between the anodes and cathodes on the corroding surface. A potential difference is then created between the cathodic protection anode and the structure such that the cathodic protection anode has a more negative potential than any point on the structure surface. Thus,
cr
the structure becomes the cathode of a new corrosion cell [4, 5]. There are two methods for applying cathodic
us
protection: sacrificial anode (galvanic) and impressed current. Each method depends on a number of economic and technical considerations. For every structure, there is a special cathodic protection system dependent on the
an
environment of the structure [6]. Current distribution in a cathodic protection system is dependent on several factors, the most important of which are the driving potential, anode and cathode geometry, spacing between the anode and cathode and the conductivity of the aqueous environment, which is favourable for a good distribution of current [2,
M
7]. Structures commonly protected are the exterior surfaces of pipelines, ships' hulls, jetties, foundation piling, steel sheet-piling and offshore platforms. Cathodic protection is also used on the interior surfaces of water-storage tanks
d
and water-circulating systems. However, because an external anode will seldom spread the protection over a
te
distance of more than two or three pipe-diameters, the method is not suitable for the protection of small-bore pipe work. [8]. There have been many laboratory studies on the corrosion of mild steel in saline water, but the application
Ac ce p
of mathematical models has been limited. This is of particular interest in developing a better scientific understanding of corrosion processes. Therefore, the present work considered the two types of cathodic protection that were studied in our previous work [9, 10] with the application of a mathematical model. This model was based on the boundary element method.
2. Mathematical Model Derivation
Corrosion engineers are interested in knowing the current and potential on the metal surfaces after two metals are electrically connected. The main objective was to provide a uniform potential distribution on the metal surface for the minimum possible power input. If the cathodic protection technique is developed with a homogeneous region Ω surrounded by a boundary Γ and with electrical conductivity k, the equation which relates the current with the potential is [11]:
2
Page 2 of 13
N ∂ci ∂E 2 I j = − F ∑ zi Di − F ∑ zi2ciui ∂x j ∂x j i =1 i =1 N
1
Defining the conductivity of the electrolyte [12]: N
k = F 2 ∑ zi2 ci ui
ip t
2
i =1
Equation 1 becomes:
I j = − F ∑ zi Di
3
us
i =1
∂ci ∂E −k ∂x j ∂x j
cr
N
The first term of equation 3 represents the portion of the current density sustained by the concentration gradient,
∂E ∂x j
Conservation of charge requires that:
∂x j
=
∂ ∂x j
− k ∂E = 0 ∂x j
5
d
∂I j
4
M
I j = −k
an
which is generally neglected in large scale simulations. Equation 3 is reduced to:
te
The conductivity is assumed to be constant for sea water; thus, equation 5 reduces to the Laplace equation for
Ac ce p
electrochemical potential [12, 13, 14, 15]:
k∇ 2 E = 0 Or:
6
∇2 E = 0
7
∂2 ∂2 ∂2 + + Where: ∇ = ∂x 2 ∂y 2 ∂z 2 2
The Laplace equation (Eq. 7) is a boundary integral formulation for the electric field that provides direct relationship between the potential and current. The solution of this equation can be obtained by applying Green's theorem, [16, 17] and this solution is readily available [17, 18]. To simplify the solution, many assumptions can be taken into account: the conductivity of sea water is assumed to be constant, each element has only one node and the shape
3
Page 3 of 13
function is constant and linear. The full details of these assumptions and the derivation of the boundary element model have been performed [8] and the following equation can be obtained:
where the nodal values
ε M e = ∑ ∑ I E * (ζ , P )φ dΓ m ζ e =1 m =1 m Γ∫ e
8
ip t
ε M c(P )E (P ) + ∑ ∑ E me ∫ I * (ζ , P )φ m dΓζ e =1 m =1 Γe
E me and I me are constants and can be brought outside the integrals. The principle of
cr
collocation means to locate the load point sequentially at all nodes of the discretization such that the domain variable at the load point E(P) coincides with the nodal value. Because linear and higher order polynomial shape functions
us
lead to nodes that belong to more than one element, it is worthwhile to introduce a global node numbering (k=1, ….,K) that does not depend on the element. If the load point is located on the first global node, the first equation of
an
the system reads:
Ĥ11
∫( φ) E * (ζ , P )dΓ + I (∫ φ) 1
1
2
E * (ζ , P1 )dΓ + ... + I K
Γ 2,e
G11
H1K K
E * (ζ , P1 )dΓ
9
Γ K ,e
G12
Ac ce p
The notation
∫( φ)
te
Γ 1, e
2
d
I1
H12
M
E1 ∫ φ1 I * (ζ , P1 )dΓ + c1 + E 2 ∫ φ 2 I * (ζ , P1 )dΓ + ... + E K ∫ φ K I * (ζ , P1 )dΓ = Γ ( 2 ,e ) Γ ( K ,e ) Γ (1,e )
G1K
∫( d) Γ means the sum of integrals contributed from these elements that contain the global node, where
Γ k ,e
Φk is the corresponding shape function. In matrix form, Eq. 9 is:
^ H 11
H 12
E1 E .... H 1K 2 = G 11 : E K
G12
I1 I .... G1K 2 : I K
10
where (^) denotes the element that contains the boundary term c(ζ). By collocation of the load point with nodes 2 to K, the additional equation of the system, Eq. 11 is obtained:
4
Page 4 of 13
H 12 ^
H 22 HK2
H 1K E G11 G12 .... G1K I1 1 .... H 2 K E 2 G21 G 22 .... G2 K I 2 = : : : : : ^ .... H KK E K G K 1 G K 2 .... G KK I K ....
11
ip t
H 11 H 21 : H K 1
And read in matrix notation:
HE = GI
cr
12
The diagonal elements of the matrices H and G contain singular integrals because the distance |ζ-P| vanishes at the
us
nodes. All other matrix elements contain regular integrals. 3. Experimental Data
an
In our previous work, the application of the cathodic protection technique was studied [9, 10]. Impressed current techniques were applied successfully. In this technique, the experimental work was conducted to determine
M
the potential and current density required in cathodic protection using polarization techniques for various conditions of temperature (0-45oC), rotating velocity (0 - 400 rpm) and pH (2-12) [9, 10]. The results of the above works were
d
used as the feed data of the mathematical model.
4. Results and discussion of the fundamental solution (E*, I*), example of the numerical solution and
te
comparison with the experimental values:
Ac ce p
For two dimensions and for simplicity, the load point is shifted to the origin, and constant elements with only one node located in the middle are chosen (i.e., M=1, Φ= 1 and c(P)=1/2). Laplace's equation of potential distribution is considered in the following example of a 2D rectangular domain with an aspect ratio 1:2, as depicted in Fig. 1. If Eq. 8 is written for the load point P, the following equation can be obtained: ε ε (ζ − P ) 1 1 ( ) E Pl + ∑ dΓ E e = ∑ − ln ζ − P dΓ I e 2 2 2π e =1 2π ζ − P e =1
13
l is the iteration of the load point, P. Four equations for four unknown boundary values are obtained if l takes the values 1 to 4 and Pl is located at the four nodes sequentially. The elements of the matrices are the integrals in Eq. 13 for different values of l and e. In matrix form:
5
Page 5 of 13
E1 E1 I1 E E I 1 2 + H 2 = G 2 E3 I 3 2 E3 E4 E 4 I 4
14
ip t
The boundary conditions were taken as a temperature of 22.5oC, rotating velocity of 0 rpm and pH 7. With these −
E 1 is the IR voltage between the cathode and anode ≈ 0, at node 2
cr
boundary conditions, at node 1 (element 1),
(element 2) the current density of the cathode ip = 61 µA/cm2, at node 3 (element 3), which is the wire connection
us
between the anode and cathode, the potential is -400 mV and at anode 4 (element 4), the potential of anode is equal to 0. Furthermore:
an
H14 = H 34 = H 32 = H12 G 14 = G 34 = G 32 = G 12
M
Spatial isotropy leads to:
17 18
Ac ce p
H 42 = H 24 and G 42 = G 24
te
And by virtue of symmetry:
H 31 = H 13 and G 31 = G 13
16
d
H 21 = H 41 = H 43 = H 23 G 21 = G 41 = G 43 = G 23
15
The main diagonal of matrix H vanishes:
H 11 = H 22 = H 33 = H 44 = 0
19
Because of symmetry:
G 33 = G 11
20
G 44 = G 22
21
Rewriting Eq. 14 in index notation and summation of the E terms leads to:
1 + H le E e = G le I e 2
22
Plugging the matrix elements and the known boundary data rearrangement leads to: 6
Page 6 of 13
E1 = 0 E2 = E 3 = −400 E 4 = 0
0.2695 − 0.0175 − 0.1119 − 0.0175 − 0.0210 0.3183 − 0.210 − 0.0420 − 0.1119 − 0.0175 − 0.2695 − 0.0175 − 0.0210 − 0.0420 − 0.0210 0.3183
I1 I = 61 2 I3 I4
ip t
0.0780 0.1621 0.9653 0.42 0.5000 0.1621 0.9653 0.5000
cr
0.5000 − 0.1621 0.9653 − 0.5000 0.0780 − 0.1621 0.9653 − 0.42
an
− 0.0780 − 0.1621 0 − 0.9653 − 0.4200 61 − 0.5000 − 0.1621 - 400 − 0.9653 − 0.5000 0
Solving the equations leads to:
Ac ce p
I 1 265 µA/cm 2 E 2 = − 816 mV I 3 − 265 µA/cm 2 2 I 4 − 58.3 µA/cm
d
0.5000 − 0.0175 − 0.9653 0.3183 − 0.0780 − 0.0175 − 0.9653 − 0.4200
I1 E 2 = I3 I4
M
− 0.1621 0.1119 0.0175 − 0.5000 0.0210 0.0420 − 0.1621 − 0.2695 0.0175 − 0.4200 0.0210 − 0.3183
te
− 0.2695 0.0210 0.1119 0.0210
us
After rewriting the equations with all unknown boundary data appearing on the left side.
Experimentally, E2=-805 mV; thus, the error is as 1.37 %. For protection, the summation of currents must be equal to zero [19]. Algebraic addition gives 61+265-58.3-265 = 2.7 µA/cm2 (too small). The results of the mathematical models are in good agreement with the experimental results. Table 1 shows the values of the experimental protection potential compared with the protection potential obtained from the boundary element model. Good agreement is obtained at all boundary conditions.
7
Page 7 of 13
Conclusions The boundary element technique was suitable for modelling corrosion problems when the surface has to be defined. The values of the potential and current density can be computed with high accuracy. Good agreement
ip t
between experimental and theoretical values was obtained. Acknowledgment
This work was supported by Baghdad University, Chemical Engineering Department, which is gratefully
cr
acknowledged.
[1].
us
References:
G.T. Parthiban, T. Parthiban, R. Ravi, V. Saraswathy, N. Palaniswamy, V. Sivan, Cathodic protection of
[2].
an
steel in concrete using magnesium alloy anode, Corros Sci. 50 (2008) 3329 – 3335. Montoya, R., Rendon, O., Genesca, J., 2009. Influence of conductivity on cathodic protection of reinforced alkali-activated slag mortar using the finite element method. Corros. Sci. 51, 2857–2862. M.H. Parsa, S.R. Allahkaram, A.H. Ghobadi, Simulation of cathodic protection potential distributions on
M
[3].
oil well casings, Jr. Petrol. Sci. Engin. 72 (2010) 215–219. [4].
O. Abootalebi, A. Kermanpur, M.R. Shishesaz, M.A. Golozar, Optimizing the electrode position in
d
sacrificial anode cathodic protection systems using boundary element method, Corros. Sci. 52 (2010) 678– 687.
J. Jezmar, Monitoring methods of cathodic protection of pipe lines, Jr. Corros. Measur. 2 (2002) 1-13.
[6].
E. Bardal, Corrosion and protection, Springer-Verlag London Limited, UK, 2004.
[7].
Douglas P. Riemer, Mark E. Orazem, A mathematical model for the cathodic protection of tank bottoms,
Ac ce p
te
[5].
Corros. Scie. 47 (2005) 849-868 [8].
K. W. Hameed, Cathodic protection of steel pipes in saline water, Ph.D. Thesis, College of Engineering, University of Baghdad, Baghdad, Iraq, 2006.
[9].
A. A. Khadom, K. W. Hameed, Prevention of Steel Corrosion by Cathodic Protection Techniques, International Journal of Chemical Technology 4 (2012) 17 – 30.
[10]. A. S. Yaro, K. W. Hameed, A. A. Khadom, Study for Prevention of Steel Corrosion by Sacrificial Anode Cathodic Protection, Theoretical Foundations of Chemical Engineering, 47 (2013) 266 – 273. [11]. J. Deconinck, Current distributions and electrode shape changes in electrochemical systems, 1st edition, Springer-Verlag, Berlin, 1992.
8
Page 8 of 13
[12]. M. Panayiotis, and C. W. Luiz, A BEM-based genetic algorithm for identification of polarization curves in cathodic protection systems, International journal for numerical method in engineering, 54 (2002) 159. [13]. R. A. Adey, Modeling of cathodic protection systems, vol. 12, WIT press, Southampton, UK, 2006.
ip t
[14]. J. C. Telles, , Santiago, J. A., A solution technique for cathodic protection with dynamic boundary conditions by the boundary element method, Advance in Engineering Software, 30 (1999) 663.
[15]. E. Chisholm, L. J. Gray, G. E. Giles, Solution of nonlinear polarization boundary conditions, Electronic
cr
Journal of Boundary Elements, 1 (2003) 418.
us
[16]. K.F. Riley, M.P. Hobson and S. J. Bence, Mathematical methods for physics and engineering, 3rd edition, Cambridge University Press, UK, 2006.
[17]. V.S. Ramanan1,M. Muthukumar2, S. Gnanasekaran, M.J. Venkataramana Reddy, B. Emmanuel, Green’s
an
functions for the Laplace equation in a 3-layer medium, boundary element integrals and their application to cathodic protection, Engin. Analy. with Boundary Elements 23 (1999) 777–786.
M
[18]. T. Milan, Laplace equation for boundary element method, 29 (2002) 10. [19]. R. A. Adey, P. Y. Hang, Computer simulation as an aid to corrosion control and reduction, National
Ac ce p
te
d
Association of Corrosion Engineering (NACE), conference 99, 1998.
9
Page 9 of 13
Symbols list:
δ( ζ −P ζ −P
) = Dirac delta function
= distance between the load point and the field point
ip t
ci = concentration of species i CO2 = concentration of oxygen
cr
D = diffusivity of dissolved oxygen Di = diffusion coefficient of species i
us
E = electrochemical potential Ea = activation energy
an
Ep = protection potential F = Faraday's constant
M
I* = fundamental current density Ij = component of the current density vector
kd = mass transfer coefficient
Ac ce p
ko = constant
te
k = rate constant of the reaction
d
J = mole flux of oxygen
N = number of species n = order of reaction
P = load point (source point) R = universal gas constant T = temperature
ui = mobility of species i zi = charge of species i
ζ = field point or influence point (receiver point) µ = seawater viscosity
10
Page 10 of 13
Figure captions:
→
ζ
→
) and load point P (marked by
P ).
Ac ce p
te
d
M
an
us
cr
common notation ζ for field point (marked by
ip t
Fig. 1 Calculation matrix elements H12 and G12 for the potential distribution in a rectangular domain. Note: the
11
Page 11 of 13
ip t
Tables:
cr
Table 1 Results of the protection potentials with different conditions X1
X2
X3
Ep, (Exp.) mV
Ep, (Theo.) mV
Absolute % error
1
9.5
85
4.11
-820
-825
0.61
2
9.5
85
9.89
-800
3
9.5
315
4.11
-825
4
9.5
315
9.89
-800
5
35.5
85
4.11
6
35.5
85
9.89
7
35.5
315
4.11
8
35.5
315
9.89
9
0.0
200
7.00
10
45.0
200
11
22.5
0
12
22.5
400
13
22.5
14
15-18
2.51
-833
0.97
-811
1.37
-825
-834
1.09
-800
-814
1.75
-830
-840
1.21
-805
-817
1.49
-805
-809
0.49
7.00
-810
-818
0.98
7.00
-805
-816
1.36
7.00
-810
-819
1.11
200
2.00
-850
-865
1.76
22.5
200
12.00
-790
-799
1.14
22.5
200
7.00
-805
-818
1.61
Ac ce p
te
d
an
-820
M
us
Run No.
o
X1 is the temperature range between 0 and 45 C, X2 is the rotating velocity range between 0 and 400 rpm, and X3 is the pH range between 2 and 12 [9, 10].
12
Page 12 of 13
ip t
Figures:
us
cr
Fig. 1
an
y
Direction of integration
Node
Γ3
2
M
Ē = -400 mV
Ī = 61 μA/cm2
d
Ē = 0 mV
Ac ce p
te
Node
Node →
→
n
ζ
→
Γ4
0
→
ζ−P
→
P
Γ1
Γ2
Node Ē ≈ 0 mV
x 1
13
Page 13 of 13
S1658-3655(15)00076-X http://dx.doi.org/doi:10.1016/j.jtusci.2015.04.002 JTUSCI 185
To appear in: Received date: Revised date: Accepted date:
8-1-2015 30-3-2015 13-4-2015
Please cite this article as: K.W. Hameed, A.S. Yaro, A.A. Khadom, Mathematical Model for Cathodic Protection in a Steel - Saline Water System, Journal of Taibah University for Science (2015), http://dx.doi.org/10.1016/j.jtusci.2015.04.002 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Mathematical Model for Cathodic Protection in a Steel - Saline Water System
Khalid W. Hameed 1, Aprael S. Yaro2 and Anees A. Khadom*,3 Biochemical Engineering Department, Al – Khawarzim Engineering College, University of Baghdad, Aljadrea
ip t
1
71001, Baghdad, Iraq 2
Chemical Engineering Department, College of Engineering, University of Baghdad, Aljadrea 71001, Baghdad,
us
*
Chemical Engineering Department, College of Engineering, University of Daiyla, Baquba 32001, Daiyla, Iraq
Correspondence author: [email protected], mobile number: 00964 790 2305786
an
3
cr
Iraq
Abstract
The corrosion of steel tubes in sea water was controlled by cathodic protection. The impressed current technique
M
was used. The rate of reaction was evaluated as a function of the temperature, pH and solution velocity. In this technique, the polarization method was used to determine the protection potential and current. The rate of zinc
d
consumption, the protection potential, and the protection current are highly dependent on the variables of the study.
te
The boundary element technique was suitable for modelling corrosion problems. The average percentage of error among the experimental and theoretical data was 1.27%.
Ac ce p
Keywords: Corrosion; Electrochemical method; Metals; Polarization; Cathodic protection 1. Introduction:
Among the various corrosion control methods available, cathodic protection is commonly adopted to control the corrosion of steel. The cathodic protection system is aimed to shift the potential of steel to the least probable range for corrosion [1]. The first attempts at cathodic protection modelling were made using finite difference approaches, but this technique is limited to two-dimensional and axisymmetric problems. Recently, the finite element method has been used successfully by some researchers [2]. Unlike other applications of the finite element method, the object of a cathodic protection system is not the structure itself but the environment close to the structure. Therefore, it is required to construct a finite element mesh of the medium between the members of the structure and to extend the mesh along the distance to establish a realistic boundary. The boundary element method requires only the discretization of the anode and cathode surfaces; therefore, the numerical problem may be reduced 1
Page 1 of 13
in size, which permits better resolution and a reduction in computer time compared to other methods, particularly for complex geometries [3]. Cathodic protection involves the application of a direct current (DC) from an anode through the electrolyte to the surface to be protected. This is often thought of as "overcoming" the corrosion currents that
ip t
exist on the structure. Cathodic protection eliminates the potential differences between the anodes and cathodes on the corroding surface. A potential difference is then created between the cathodic protection anode and the structure such that the cathodic protection anode has a more negative potential than any point on the structure surface. Thus,
cr
the structure becomes the cathode of a new corrosion cell [4, 5]. There are two methods for applying cathodic
us
protection: sacrificial anode (galvanic) and impressed current. Each method depends on a number of economic and technical considerations. For every structure, there is a special cathodic protection system dependent on the
an
environment of the structure [6]. Current distribution in a cathodic protection system is dependent on several factors, the most important of which are the driving potential, anode and cathode geometry, spacing between the anode and cathode and the conductivity of the aqueous environment, which is favourable for a good distribution of current [2,
M
7]. Structures commonly protected are the exterior surfaces of pipelines, ships' hulls, jetties, foundation piling, steel sheet-piling and offshore platforms. Cathodic protection is also used on the interior surfaces of water-storage tanks
d
and water-circulating systems. However, because an external anode will seldom spread the protection over a
te
distance of more than two or three pipe-diameters, the method is not suitable for the protection of small-bore pipe work. [8]. There have been many laboratory studies on the corrosion of mild steel in saline water, but the application
Ac ce p
of mathematical models has been limited. This is of particular interest in developing a better scientific understanding of corrosion processes. Therefore, the present work considered the two types of cathodic protection that were studied in our previous work [9, 10] with the application of a mathematical model. This model was based on the boundary element method.
2. Mathematical Model Derivation
Corrosion engineers are interested in knowing the current and potential on the metal surfaces after two metals are electrically connected. The main objective was to provide a uniform potential distribution on the metal surface for the minimum possible power input. If the cathodic protection technique is developed with a homogeneous region Ω surrounded by a boundary Γ and with electrical conductivity k, the equation which relates the current with the potential is [11]:
2
Page 2 of 13
N ∂ci ∂E 2 I j = − F ∑ zi Di − F ∑ zi2ciui ∂x j ∂x j i =1 i =1 N
1
Defining the conductivity of the electrolyte [12]: N
k = F 2 ∑ zi2 ci ui
ip t
2
i =1
Equation 1 becomes:
I j = − F ∑ zi Di
3
us
i =1
∂ci ∂E −k ∂x j ∂x j
cr
N
The first term of equation 3 represents the portion of the current density sustained by the concentration gradient,
∂E ∂x j
Conservation of charge requires that:
∂x j
=
∂ ∂x j
− k ∂E = 0 ∂x j
5
d
∂I j
4
M
I j = −k
an
which is generally neglected in large scale simulations. Equation 3 is reduced to:
te
The conductivity is assumed to be constant for sea water; thus, equation 5 reduces to the Laplace equation for
Ac ce p
electrochemical potential [12, 13, 14, 15]:
k∇ 2 E = 0 Or:
6
∇2 E = 0
7
∂2 ∂2 ∂2 + + Where: ∇ = ∂x 2 ∂y 2 ∂z 2 2
The Laplace equation (Eq. 7) is a boundary integral formulation for the electric field that provides direct relationship between the potential and current. The solution of this equation can be obtained by applying Green's theorem, [16, 17] and this solution is readily available [17, 18]. To simplify the solution, many assumptions can be taken into account: the conductivity of sea water is assumed to be constant, each element has only one node and the shape
3
Page 3 of 13
function is constant and linear. The full details of these assumptions and the derivation of the boundary element model have been performed [8] and the following equation can be obtained:
where the nodal values
ε M e = ∑ ∑ I E * (ζ , P )φ dΓ m ζ e =1 m =1 m Γ∫ e
8
ip t
ε M c(P )E (P ) + ∑ ∑ E me ∫ I * (ζ , P )φ m dΓζ e =1 m =1 Γe
E me and I me are constants and can be brought outside the integrals. The principle of
cr
collocation means to locate the load point sequentially at all nodes of the discretization such that the domain variable at the load point E(P) coincides with the nodal value. Because linear and higher order polynomial shape functions
us
lead to nodes that belong to more than one element, it is worthwhile to introduce a global node numbering (k=1, ….,K) that does not depend on the element. If the load point is located on the first global node, the first equation of
an
the system reads:
Ĥ11
∫( φ) E * (ζ , P )dΓ + I (∫ φ) 1
1
2
E * (ζ , P1 )dΓ + ... + I K
Γ 2,e
G11
H1K K
E * (ζ , P1 )dΓ
9
Γ K ,e
G12
Ac ce p
The notation
∫( φ)
te
Γ 1, e
2
d
I1
H12
M
E1 ∫ φ1 I * (ζ , P1 )dΓ + c1 + E 2 ∫ φ 2 I * (ζ , P1 )dΓ + ... + E K ∫ φ K I * (ζ , P1 )dΓ = Γ ( 2 ,e ) Γ ( K ,e ) Γ (1,e )
G1K
∫( d) Γ means the sum of integrals contributed from these elements that contain the global node, where
Γ k ,e
Φk is the corresponding shape function. In matrix form, Eq. 9 is:
^ H 11
H 12
E1 E .... H 1K 2 = G 11 : E K
G12
I1 I .... G1K 2 : I K
10
where (^) denotes the element that contains the boundary term c(ζ). By collocation of the load point with nodes 2 to K, the additional equation of the system, Eq. 11 is obtained:
4
Page 4 of 13
H 12 ^
H 22 HK2
H 1K E G11 G12 .... G1K I1 1 .... H 2 K E 2 G21 G 22 .... G2 K I 2 = : : : : : ^ .... H KK E K G K 1 G K 2 .... G KK I K ....
11
ip t
H 11 H 21 : H K 1
And read in matrix notation:
HE = GI
cr
12
The diagonal elements of the matrices H and G contain singular integrals because the distance |ζ-P| vanishes at the
us
nodes. All other matrix elements contain regular integrals. 3. Experimental Data
an
In our previous work, the application of the cathodic protection technique was studied [9, 10]. Impressed current techniques were applied successfully. In this technique, the experimental work was conducted to determine
M
the potential and current density required in cathodic protection using polarization techniques for various conditions of temperature (0-45oC), rotating velocity (0 - 400 rpm) and pH (2-12) [9, 10]. The results of the above works were
d
used as the feed data of the mathematical model.
4. Results and discussion of the fundamental solution (E*, I*), example of the numerical solution and
te
comparison with the experimental values:
Ac ce p
For two dimensions and for simplicity, the load point is shifted to the origin, and constant elements with only one node located in the middle are chosen (i.e., M=1, Φ= 1 and c(P)=1/2). Laplace's equation of potential distribution is considered in the following example of a 2D rectangular domain with an aspect ratio 1:2, as depicted in Fig. 1. If Eq. 8 is written for the load point P, the following equation can be obtained: ε ε (ζ − P ) 1 1 ( ) E Pl + ∑ dΓ E e = ∑ − ln ζ − P dΓ I e 2 2 2π e =1 2π ζ − P e =1
13
l is the iteration of the load point, P. Four equations for four unknown boundary values are obtained if l takes the values 1 to 4 and Pl is located at the four nodes sequentially. The elements of the matrices are the integrals in Eq. 13 for different values of l and e. In matrix form:
5
Page 5 of 13
E1 E1 I1 E E I 1 2 + H 2 = G 2 E3 I 3 2 E3 E4 E 4 I 4
14
ip t
The boundary conditions were taken as a temperature of 22.5oC, rotating velocity of 0 rpm and pH 7. With these −
E 1 is the IR voltage between the cathode and anode ≈ 0, at node 2
cr
boundary conditions, at node 1 (element 1),
(element 2) the current density of the cathode ip = 61 µA/cm2, at node 3 (element 3), which is the wire connection
us
between the anode and cathode, the potential is -400 mV and at anode 4 (element 4), the potential of anode is equal to 0. Furthermore:
an
H14 = H 34 = H 32 = H12 G 14 = G 34 = G 32 = G 12
M
Spatial isotropy leads to:
17 18
Ac ce p
H 42 = H 24 and G 42 = G 24
te
And by virtue of symmetry:
H 31 = H 13 and G 31 = G 13
16
d
H 21 = H 41 = H 43 = H 23 G 21 = G 41 = G 43 = G 23
15
The main diagonal of matrix H vanishes:
H 11 = H 22 = H 33 = H 44 = 0
19
Because of symmetry:
G 33 = G 11
20
G 44 = G 22
21
Rewriting Eq. 14 in index notation and summation of the E terms leads to:
1 + H le E e = G le I e 2
22
Plugging the matrix elements and the known boundary data rearrangement leads to: 6
Page 6 of 13
E1 = 0 E2 = E 3 = −400 E 4 = 0
0.2695 − 0.0175 − 0.1119 − 0.0175 − 0.0210 0.3183 − 0.210 − 0.0420 − 0.1119 − 0.0175 − 0.2695 − 0.0175 − 0.0210 − 0.0420 − 0.0210 0.3183
I1 I = 61 2 I3 I4
ip t
0.0780 0.1621 0.9653 0.42 0.5000 0.1621 0.9653 0.5000
cr
0.5000 − 0.1621 0.9653 − 0.5000 0.0780 − 0.1621 0.9653 − 0.42
an
− 0.0780 − 0.1621 0 − 0.9653 − 0.4200 61 − 0.5000 − 0.1621 - 400 − 0.9653 − 0.5000 0
Solving the equations leads to:
Ac ce p
I 1 265 µA/cm 2 E 2 = − 816 mV I 3 − 265 µA/cm 2 2 I 4 − 58.3 µA/cm
d
0.5000 − 0.0175 − 0.9653 0.3183 − 0.0780 − 0.0175 − 0.9653 − 0.4200
I1 E 2 = I3 I4
M
− 0.1621 0.1119 0.0175 − 0.5000 0.0210 0.0420 − 0.1621 − 0.2695 0.0175 − 0.4200 0.0210 − 0.3183
te
− 0.2695 0.0210 0.1119 0.0210
us
After rewriting the equations with all unknown boundary data appearing on the left side.
Experimentally, E2=-805 mV; thus, the error is as 1.37 %. For protection, the summation of currents must be equal to zero [19]. Algebraic addition gives 61+265-58.3-265 = 2.7 µA/cm2 (too small). The results of the mathematical models are in good agreement with the experimental results. Table 1 shows the values of the experimental protection potential compared with the protection potential obtained from the boundary element model. Good agreement is obtained at all boundary conditions.
7
Page 7 of 13
Conclusions The boundary element technique was suitable for modelling corrosion problems when the surface has to be defined. The values of the potential and current density can be computed with high accuracy. Good agreement
ip t
between experimental and theoretical values was obtained. Acknowledgment
This work was supported by Baghdad University, Chemical Engineering Department, which is gratefully
cr
acknowledged.
[1].
us
References:
G.T. Parthiban, T. Parthiban, R. Ravi, V. Saraswathy, N. Palaniswamy, V. Sivan, Cathodic protection of
[2].
an
steel in concrete using magnesium alloy anode, Corros Sci. 50 (2008) 3329 – 3335. Montoya, R., Rendon, O., Genesca, J., 2009. Influence of conductivity on cathodic protection of reinforced alkali-activated slag mortar using the finite element method. Corros. Sci. 51, 2857–2862. M.H. Parsa, S.R. Allahkaram, A.H. Ghobadi, Simulation of cathodic protection potential distributions on
M
[3].
oil well casings, Jr. Petrol. Sci. Engin. 72 (2010) 215–219. [4].
O. Abootalebi, A. Kermanpur, M.R. Shishesaz, M.A. Golozar, Optimizing the electrode position in
d
sacrificial anode cathodic protection systems using boundary element method, Corros. Sci. 52 (2010) 678– 687.
J. Jezmar, Monitoring methods of cathodic protection of pipe lines, Jr. Corros. Measur. 2 (2002) 1-13.
[6].
E. Bardal, Corrosion and protection, Springer-Verlag London Limited, UK, 2004.
[7].
Douglas P. Riemer, Mark E. Orazem, A mathematical model for the cathodic protection of tank bottoms,
Ac ce p
te
[5].
Corros. Scie. 47 (2005) 849-868 [8].
K. W. Hameed, Cathodic protection of steel pipes in saline water, Ph.D. Thesis, College of Engineering, University of Baghdad, Baghdad, Iraq, 2006.
[9].
A. A. Khadom, K. W. Hameed, Prevention of Steel Corrosion by Cathodic Protection Techniques, International Journal of Chemical Technology 4 (2012) 17 – 30.
[10]. A. S. Yaro, K. W. Hameed, A. A. Khadom, Study for Prevention of Steel Corrosion by Sacrificial Anode Cathodic Protection, Theoretical Foundations of Chemical Engineering, 47 (2013) 266 – 273. [11]. J. Deconinck, Current distributions and electrode shape changes in electrochemical systems, 1st edition, Springer-Verlag, Berlin, 1992.
8
Page 8 of 13
[12]. M. Panayiotis, and C. W. Luiz, A BEM-based genetic algorithm for identification of polarization curves in cathodic protection systems, International journal for numerical method in engineering, 54 (2002) 159. [13]. R. A. Adey, Modeling of cathodic protection systems, vol. 12, WIT press, Southampton, UK, 2006.
ip t
[14]. J. C. Telles, , Santiago, J. A., A solution technique for cathodic protection with dynamic boundary conditions by the boundary element method, Advance in Engineering Software, 30 (1999) 663.
[15]. E. Chisholm, L. J. Gray, G. E. Giles, Solution of nonlinear polarization boundary conditions, Electronic
cr
Journal of Boundary Elements, 1 (2003) 418.
us
[16]. K.F. Riley, M.P. Hobson and S. J. Bence, Mathematical methods for physics and engineering, 3rd edition, Cambridge University Press, UK, 2006.
[17]. V.S. Ramanan1,M. Muthukumar2, S. Gnanasekaran, M.J. Venkataramana Reddy, B. Emmanuel, Green’s
an
functions for the Laplace equation in a 3-layer medium, boundary element integrals and their application to cathodic protection, Engin. Analy. with Boundary Elements 23 (1999) 777–786.
M
[18]. T. Milan, Laplace equation for boundary element method, 29 (2002) 10. [19]. R. A. Adey, P. Y. Hang, Computer simulation as an aid to corrosion control and reduction, National
Ac ce p
te
d
Association of Corrosion Engineering (NACE), conference 99, 1998.
9
Page 9 of 13
Symbols list:
δ( ζ −P ζ −P
) = Dirac delta function
= distance between the load point and the field point
ip t
ci = concentration of species i CO2 = concentration of oxygen
cr
D = diffusivity of dissolved oxygen Di = diffusion coefficient of species i
us
E = electrochemical potential Ea = activation energy
an
Ep = protection potential F = Faraday's constant
M
I* = fundamental current density Ij = component of the current density vector
kd = mass transfer coefficient
Ac ce p
ko = constant
te
k = rate constant of the reaction
d
J = mole flux of oxygen
N = number of species n = order of reaction
P = load point (source point) R = universal gas constant T = temperature
ui = mobility of species i zi = charge of species i
ζ = field point or influence point (receiver point) µ = seawater viscosity
10
Page 10 of 13
Figure captions:
→
ζ
→
) and load point P (marked by
P ).
Ac ce p
te
d
M
an
us
cr
common notation ζ for field point (marked by
ip t
Fig. 1 Calculation matrix elements H12 and G12 for the potential distribution in a rectangular domain. Note: the
11
Page 11 of 13
ip t
Tables:
cr
Table 1 Results of the protection potentials with different conditions X1
X2
X3
Ep, (Exp.) mV
Ep, (Theo.) mV
Absolute % error
1
9.5
85
4.11
-820
-825
0.61
2
9.5
85
9.89
-800
3
9.5
315
4.11
-825
4
9.5
315
9.89
-800
5
35.5
85
4.11
6
35.5
85
9.89
7
35.5
315
4.11
8
35.5
315
9.89
9
0.0
200
7.00
10
45.0
200
11
22.5
0
12
22.5
400
13
22.5
14
15-18
2.51
-833
0.97
-811
1.37
-825
-834
1.09
-800
-814
1.75
-830
-840
1.21
-805
-817
1.49
-805
-809
0.49
7.00
-810
-818
0.98
7.00
-805
-816
1.36
7.00
-810
-819
1.11
200
2.00
-850
-865
1.76
22.5
200
12.00
-790
-799
1.14
22.5
200
7.00
-805
-818
1.61
Ac ce p
te
d
an
-820
M
us
Run No.
o
X1 is the temperature range between 0 and 45 C, X2 is the rotating velocity range between 0 and 400 rpm, and X3 is the pH range between 2 and 12 [9, 10].
12
Page 12 of 13
ip t
Figures:
us
cr
Fig. 1
an
y
Direction of integration
Node
Γ3
2
M
Ē = -400 mV
Ī = 61 μA/cm2
d
Ē = 0 mV
Ac ce p
te
Node
Node →
→
n
ζ
→
Γ4
0
→
ζ−P
→
P
Γ1
Γ2
Node Ē ≈ 0 mV
x 1
13
Page 13 of 13