- Email: [email protected]
Using the Gaussian function to simulate constant potential anodes in multiobjective optimization of cathodic protection systems
Engineering Analysis with Boundary Elements 73 (2016) 35–41
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
Using the Gaussian function to simulate constant potential anodes in multiobjective optimization of cathodic protection systems W.J. Santos, J.A.F. Santiago n, J.C.F. Telles Department of Civil Engineering, COPPE/UFRJ, Caixa Postal 68506, CEP21941-972, Rio de Janeiro, RJ, Brazil
art ic l e i nf o
a b s t r a c t
Article history: Received 22 March 2016 Received in revised form 7 July 2016 Accepted 31 August 2016
The purpose of this work is to numerically find the optimum location of constant potential anodes to ensure complete structure surface protection using a cathodic protection technique. The existence of sacrificial anodes is originally introduced through the boundary conditions of the corresponding boundary value problem (BVP). However, if constant potential galvanic regions are introduced through its boundaries, then finding their optimal location is not an easy task due to the necessity of redefining boundary geometric nodes and the arrangement of virtual sources for the standard method of fundamental solutions (MFS) formulation. Therefore, in this work, the galvanic anodes are introduced as source terms using a Gaussian function. Hence, the boundary remains the same for different anode positions. The optimization process includes the identification of the following parameters characterizing the Gaussian function: the optimum coordinates of the centre of the anode, a factor that involves the inherent potential of the electrode and a proportionality factor for the electrode diameter. The MFS methodology coupled with a genetic algorithm presented good results for this multiobjective optimization procedure. This fact can be seen in the several results of applications that are discussed in this paper, considering numerical simulations in finite regions in R2. & 2016 Elsevier Ltd. All rights reserved.
Keywords: Genetic algorithms Method of fundamental solutions Dual reciprocity
1. Introduction Cathodic protection (CP) is a technique to reduce the corrosion rate of the metal surface by making it the cathode of an electrochemical cell. In practice, the goal is to provide a uniform potential distribution on the metal surfaces, limited by a minimum potential value to guarantee protection from the corrosion using an external power source (impressed current CP) or by utilizing a sacrificial anode [1,2]. In this work, it is considered the case of sacrificial anode system, where the galvanic relationship between a sacrificial anode material, such as zinc or magnesium, and the structure is used to supply the required CP current. The most commonly used methods for modelling cathodic protection systems are finite element method (FEM) and boundary element method (BEM). However, the method of fundamental solutions (MFS) is a technique which can also be applied to CP problems (see references [3,4]). Just like BEM, MFS is applicable when a fundamental solution of the differential equation in question is known, with the advantage of not requiring any integration procedure or specific treatment for the singularities of n
Corresponding author. E-mail addresses: [email protected] (W.J. Santos), [email protected] (J.A.F. Santiago), [email protected] (J.C.F. Telles). http://dx.doi.org/10.1016/j.enganabound.2016.08.014 0955-7997/& 2016 Elsevier Ltd. All rights reserved.
the fundamental solution. The BEM and the MFS are the most appropriate techniques to solve problems involving galvanic corrosion and CP systems, mainly to solve large problems and considering homogeneous conductive medium. These methods require only the representation of anodes and cathodes surfaces, which leads to better resolution and reduction in computer run time when compared to FEM. Mathematical simulations of cathodic protections systems using FEM can be seen in [5–7], where compact support functions were used to simulate constant potential electrodes in electrochemical process. Due to its accuracy and simplicity of mesh generation, the BEM is usually used for numerical simulations of sacrificial anode CP systems. Miyasaka et al. [8], for example, evaluated the computational accuracy of BEM to estimate the galvanic corrosion and CP in an actual field. Abootalebi et al. [9] determined the optimum location of zinc anode electrode using BEM in 2D. In addition, the influence of anode length and paint defect on corrosion current density and potential distributions of sacrificial anode CP system were investigated. Several different applications of BEM to study CP systems have been reported in the literature, including reference to practical analyses performed by offshore oil companies [10–12]. The BEM implementation includes a Newton–Raphson solution algorithm to accommodate possible nonlinear boundary conditions [13]. Coupled with the numerical method to solve the Laplace equation, optimization algorithms can be used to
36
W.J. Santos et al. / Engineering Analysis with Boundary Elements 73 (2016) 35–41
determine the optimum location and the corresponding current intensity values of the anodes in order to satisfy a protection criterion. Hence, the minimization of an objective function using, for example, genetic algorithms (GAs) and a penalty method for handling constraints can be adopted. This type of optimization can and has also been successfully performed using BEM [14,15]. Kupradze and Aleksidze [16] first proposed the basic ideas for the formulation of the MFS. In order to construct the solution, the MFS uses only a superposition of fundamental solutions associated to the problem, with singular points (virtual sources) located outside the domain. The accuracy of the MFS numerical solution depends on the radius of such a circle or on the distance from the virtual sources over the geometrically similar boundary contour to the problem boundary, especially due to possible ill-conditioning and/or rank-deficiency of the algebraic system of equations formed. The MFS has successfully been applied for solving several problems. For example, Costa et al. [17] developed numerical frequency domain formulations to simulate the 2D acoustic wave propagation in the vicinity of an underwater configuration which combines two sub-regions using the MFS. More recently, Fontes Jr. et al. [18] applied a regularized method of fundamental solutions coupled with the numerical Green's function procedure to solve linear elastic fracture mechanics (LEFM) problems. Problems with nonlinear boundary conditions solved by MFS can be treated as nonlinear least squares problems [19]. Therefore, to determine the coefficients of the linear superposition of the fundamental solutions and the positions of the virtual sources, a nonlinear least squares algorithm is found necessary. In the present work, the minimization of the nonlinear functional is done using the MINPACK [20] routine LMDIF, which is a modified version of the Levenberg–Marquardt algorithm [21]. The Levenberg–Marquardt method has been highly recommended when Jacobian is rankdeficiency or nearly so. Santos et al. [3] were the first to use standard MFS successfully in the numerical simulations of CP systems. In the paper cited, the authors proposed a GA with the MFS to simulate cathodic protection systems with nonlinear boundary conditions. The adopted GA was used to minimize a nonlinear error function, whose design variables were the coefficients of the linear superposition of fundamental solutions and the positions of the virtual sources, randomly distributed outside the problem domain. In Santos et al. [4] was presented a formulation using a GA and the MFS to determine the optimum location and the optimum current intensity of the anodes inserted in the electrolyte leading to a practical optimized design procedure. The results presented in this paper included a comparison with a direct boundary element (BEM) solution procedure. The purpose of the present paper is to use a genetic algorithm (GA) with the MFS for optimizing the anode position in a galvanic cathodic protection system. Furthermore, the constant potential circular anodes are here mathematically represented by a Gaussian function. Thus, the GA will be used to search the Gaussian function parameters: coordinates of the centre of the anode, a factor that involves the inherent potential of the anode and a proportionality factor for the anode diameter. The main advantage of considering constant potential regions as the source term is that the boundary conditions remain the same, no matter where the anodes are located. The circular anodes are considered and, therefore, a good approximation for the source term is a function with circular compact support like the Gaussian function. In the MFS case, it is not necessary for the arrangement of the virtual sources inside the galvanic anodes (cavities), when they are introduced as source terms. A particular solution to the inhomogeneous differential equation resulting will be evaluated using the dual reciprocity method (DRM). In the DRM, the source term is approximated by a finite series of radial basis functions with an approximation to
particular solution calculated analytically from source [22]. The corresponding homogeneous solution is found by MFS. This text is organized as follows: in Section 2 is presented a standard electrochemical potential problem, i.e., the Laplace equation with boundary conditions given by constant potential (anode) and polarization curve (cathode). The decision of introducing the Gaussian function to model constant potential galvanic regions is also discussed in Section 2. The standard MFS formulation with the dual reciprocity method for solving the Poisson equation with nonlinear boundary conditions is discussed in Section 3. The multiobjective optimization problem and the genetic algorithm are showed in Section 4. This algorithm is coupled with MFS to guarantee protection from the corrosion by utilizing a sacrificial anode in the numerical examples solved in Section 5. The paper ends with some discussions found in Section 6.
2. Boundary value problem (BVP) For a homogeneous and isotropic electrolyte (domain Ω) system of conductivity k, as illustrated in Fig. 1, the electrochemical potential problem studied obeys the Laplace equation:
k∇2ϕ(x) = 0,
x ∈ Ω,
(1)
subjected to the following boundary conditions:
i(x) = F (ϕ),
ϕ(x) = ϕ0,
x ∈ Γ1,
x ∈ Γ2,
(2)
(3)
where Γ ≡ Γ1 ∪ Γ2 is the boundary of Ω, ϕ0 is the fixed potential of anode, i(x) is the current density in the outward normal direction n and F (ϕ) is a nonlinear function of ϕ given here by the following polarization curve [23]:
i = F (ϕ) = e
ϕ + 693.91 β1
−1 ϕ + 521.6 ⎤ ϕ + 707.57 ⎡1 − β3 − ⎢ + e β2 ⎥ − e , ⎣ i1 ⎦
(4)
with ϕ and i having units mV and μA/cm2, respectively, and β1, β2, β3 and i1 are given constant parameters: β1 = 24 mV ,
Fig. 1. The original potential problem.
W.J. Santos et al. / Engineering Analysis with Boundary Elements 73 (2016) 35–41
β2 = 23.47 mV , β3 = 55 mV e i1 = 86.06 μA/cm2 . The conductivity of the electrolyte is equal to k = 0.0479 Ω−1 cm−1 and the critical value of the electrochemical potential is ϕc = − 850 mV (vs. SCE). The boundary value problem (1)–(3) can be modified in order to simplify the optimization of positioning of anodes for cathodic protection. For this, it had been decided to consider the constant potential condition ϕ0 for Γ2 on the right-hand side of Poisson equation and not in the boundary condition. Thus, the corresponding BVP can be written in the following form:
k∇2ϕ(x) = b(x), x ∈ Ω,
(5)
with boundary condition given by Eq. (2). The term b presented in Eq. (5) is a function representing the anode as external source. Here, the source term b(x) is mathematically represented by a Gaussian function, defined in two-dimensional space (x,y) by
b(x, y) = Ae−nn(x − x0)
2
− nn(y − y0 )2
,
(6)
where A is height of the Gaussian function, the coordinates x0 and y0 are the centre of this function and nn is a proportionality factor that controls the width of the Gaussian. The volume (V ) under the Gaussian function is given by
V=
Aπ , ∫Ω b(x, y)dΩ = nn
The integral in Eq. (7) is 1 if A = and, for this case, the Gaussian function is a probability density function of a normal distribution that can be written as
nn −nn(x − x0)2 − nn(y − y0 )2 b(x, y) = e . π
(8)
In Eq. (8), if the parameter nn tends to infinity, then the function b(x, y ) tends to Dirac δ “function”. Treating the anode as Gaussian function, the term b(x) is equal to expression (8) multiplied by density of the source P (x0, y0 ). Therefore:
b(x) = b(x, y) = P (x 0 , y0 )
nn −nn(x − x0)2 − nn(y − y0 )2 e . π
(9)
The conservation of current between anode and cathode in the original problem (see Fig. 1) is automatically obtained by satisfying the condition
∫Γ i(x)dΓ = − ∫Γ i(x)dΓ . 1
2
(10)
On the other hand, using Eqs. (2) and (5), the corresponding BVP satisfies the following condition:
∫Γ i(x)dΓ = − P(x0, y0 ), 1
(11)
which ensures that there is no loss of current. To impose equivalence between the problems, i.e., between Eqs. (10) and (11), it is necessary to guarantee that
P (x 0 , y0 ) =
∫Γ i(x)dΓ . 2
ϕ(x) = ϕp(x) + ϕh(x),
(13)
where ϕp(x) satisfies the inhomogeneous equation (5) but does not necessarily satisfy the nonlinear boundary condition presented in Eq. (2) and ϕh(x) satisfies the Laplace equation and the corresponding boundary condition for the homogeneous problem, which in the nonlinear case, is given by
ih(x) = F (ϕp + ϕh) − ip(x),
x ∈ Γ1,
(14)
ϕp is a particular solution and ip(x) is its normal derivative. ∂ϕ Once ϕp and p are known, ϕh(x) can be determined by standard ∂n where
MFS. The determination of the particular solution depends on the type of source b(x) in (5). In the simple cases, as point sources for example, ϕp(x) can be evaluated analytically. For general cases, the particular solutions can be determined by numerical techniques [24,25]. Here, a particular solution will be derived using the context of the dual reciprocity method using a positive definite compactly supported radial basis function (CS-RBF) of Wendland [26]. Thus, the function b(x) is approximated by a finite series of radial basis functions ψ (r ) with an approximation ϕ^ to ϕ calcup
k∇2ϕ^p(x) =
N+L
N+L
∑ ψj(x)aj = ∑ ψj( j=1
x − x j )aj , (15)
j=1
where N are the boundary nodes, L are the internal nodes and the unknown coefficients aj 's are determined solving the system of equations N+L
∑ ψj(x)aj = b^(x),
(16)
j=1
^ with the point values of the known function b denoted for b . For practical implementation, the Wendland radial basis function of compact support is scaled with scale β as follows:
⎛r⎞ ψ ( x − x j ) = ψ ⎜ ⎟. ⎝ β⎠
(17)
The system (16) is symmetric and can be dense or sparse depending on the support parameter β. The accuracy and efficiency of the Wendland CS-RBF depends on the scale of the support and determining the scale of support can be uncertain. For this reason, the singular value decomposition (SVD) procedure is used to provide acceptable solutions to possible ill-conditioned system of equations generated by the Wendland CS-RBF. After determining the coefficients aj 's using (16) an approximation to a particular solution can be written as
ϕ^p(x) =
N+L
∑ Ψj(x)aj,
(18)
j=1
(12)
Eq. (12) shows that the determination of P depends not only on the potential at which the anode is set, but also on the geometry of the domain and the location of coordinates x0 and y0. The diameter of the anode is related to the parameter nn in Eq. (9).
3. Method of fundamental solutions for poisson's equation In order to solve the corresponding BVP, Eqs. (2) and (5), using MFS, it is necessary to eliminate the inhomogeneous term in (5) by the use of a particular solution. Let
p
lated analytically from ψ [27]. Mathematically, the following differential equation is solved:
(7) nn π
37
where
Ψj is the solution of
2
k∇ Ψj(x) = ψj(x).
(19)
The solution of Eq. (19) can easily be achieved considering the radial part of Laplacian operator and by analytical integration of Eq. (19). In this paper, the following Wendland radial basis function is chosen: 6⎛ ⎛ ⎞ ⎛ ⎞ r r⎞ r2 r ψ ⎜⎜ ⎟⎟ = ⎜ 1 − ⎟ ⎜⎜ 35 2 + 18 + 3⎟⎟. β⎠ ⎝ β β ⎝ β⎠ ⎝ ⎠
An explicit analytical representation of
(20) r Ψ(β)
corresponding to
38
W.J. Santos et al. / Engineering Analysis with Boundary Elements 73 (2016) 35–41
Wendland's CS-RBF presented in Eq. (20) can be derived by straightforward integration [28]:
4. Multiobjective optimization using GA
r Ψ( ) β =
1 k ⎧ 7r10 64r 9 105r 8 64r 7 35r 6 7r 4 3r 2 ⎪ − + − + − + , 2 8 7 6 5 4 4 ⎪ 20β 4β 27β 16β 7β 6β ⎪ ⎨ r ≤ β, ⎪ ⎪ 3517β 2 β2 ⎛ r ⎞ ⎪ 15120 + 6 ln⎜ β ⎟, r > β. ⎝ ⎠ ⎩
(21)
After obtaining Ψj by Eq. (19), it is possible to determine a particular solution using Eq. (18) and, consequently, the approximate solution of the associated homogeneous problem by MFS. The method of fundamental solutions establishes that the approximate solution can be constructed by a summation of similar problem solution given by the following superposition: nvs
ϕh(x) =
∑ G(x,
xvs j )cj,
j=1
(22)
with nvs being the total number of virtual sources, xvs j represents the position of the virtual sources and the coefficients cj 's are the unknown source intensities. The function G(x, xvs j ) is a fundamental solution of Laplace's equation given by
G(x,
xvs j )
1 1 = ln , 2π k R
(23)
where R is the Euclidean distance between point xvs and the field point x . ∂G Similarly, defining H = k ∂n , the homogeneous solution for the current density (ih) is given as nvs
ih(x) =
user-specified maximum number of function evaluations is reached.
∑ H (x, xvsj )cj. j=1
(24)
The idea of MFS is to determine the coefficients cj 's by imposing satisfaction of the boundary conditions at certain boundary nodes. The boundary condition given by Eq. (14) generates a nonlinear least squares problem in which the design variables are the coefficients cj 's and the positions of the virtual sources. In least squares problems, the objective function f has the following special form:
1 f (c, xvs) = 2
nbn
∑ ϵ2j (c, xvs), j=1
(25)
where nbn is the number of boundary nodes, c is a vector containing the coefficients cj 's and each ϵj is a smooth function referred to as a residual given by
ϵj = igj − F (ϕgj ),
j = 1, …, nbn .
(26)
Considering the arrangement of the virtual sources on a circular contour, it is only necessary to search for the radius ρ of such a circle. Thus, design variables of Eq. (25) are the coefficients cj 's and a radius ρ. In this paper, the minimization of the functional of Eq. (25) is carried out using the MINPACK routine LMDIF, which is a modified version of the Levenberg–Marquardt algorithm. In LMDIF, the Jacobian is evaluated internally by finite differences. LMDIF terminates when either a user-specified tolerance is achieved or the
In this work, there are two optimization goals: (1) to pre-establish the protection criterion of metallic structures and (2) to guarantee the equivalence between the original BVP (1)–(3) and the corresponding BVP, Eqs. (2) and (5). In other words, it is necessary to search the parameters of Gaussian function in Eq. (9) in order to satisfy this multiobjective optimization. The aim of CP systems is to reduce the potential over the structure surface below a critical value (ϕc ). By optimizing the sacrificial anode (Γ2) location, the maximum corrosion current over the structure surface will be minimized in such a way as to satisfy, as uniformly as possible, close to the critical potential value: ϕ ≤ ϕc . The optimum anode location can be achieved by minimizing the proposed objective function:
z1(x 0 , y0 ) =
1 nbn
nbn
2
∑ ⎡⎣ ϕi − ϕc⎤⎦
+v
i=1
1 * nbn
nbn
∑ kl2(x 0, y0 ) , i=1
(27)
where z1 calculates the root mean square error (RMSE) between the electrochemical potential at each boundary node and the cri* is the number of boundary nodes tical potential. Furthermore, nbn that do not satisfy the protection criterion, the constant v is a penalty number and function kl is equal to:
kl(x 0 , y0 ) = (ϕi − ϕc )u(ϕi − ϕc ),
(28)
where u is the unit step function. Typical values of v for the penalty method are within the range 102–105 [29]. In addition, it is also necessary to search the optimum solution for the Gaussian function parameters nn and P (x0, y0 ). The diameter of the anode is related to the parameter nn and the height of the Gaussian function P is related to the constant potential condition, but also on the geometry of the domain and on the location of coordinates (x0, y0 ). The boundary nodes in Γ2 are internal nodes in corresponding BVP and the potential value in these nodes must be the anode prescribed potential ϕ0. Thus, the second objective function z2 is modelled in order that the internal potential values corresponding to Γ2 obey the Dirichlet condition of the original BVP (1)–(3). For this, it was decided to minimize the following objective function:
z2(nn, P ) =
1 nint
nint
2
∑ ⎡⎣ ϕinti − ϕ0⎤⎦
,
(29)
i=1
where nint is the number of internal nodes that represent the boundary of the anode in the original problem and ϕint are the potential values at these internal nodes solved for each pair (nn,P). In this optimization problem there are two objectives functions to be optimized: z1 and z2. There are some classical methods for multiobjective optimization [30]. Here, a method of objective weighting is applied, where multiple objective functions zi are combined into one overall objective function Z:
Z (x 0 , y0 , nn, P ) = w1z1 + w2z2,
(30) 2 ∑i = 1 wi
where the weights wi ∈ [0, 1] and = 1. In this paper, the minimization of Eq. (30) is achieved using a typical genetic algorithm (GA) with binary representation [31]. However, some characteristics have been included as the twopoint crossover, the elitism and the probabilities of mutation and crossover can vary linearly over the generations [32]. The weighting parameters w1 and w2 are chosen through the results of
W.J. Santos et al. / Engineering Analysis with Boundary Elements 73 (2016) 35–41
the multiobjective optimization process, which are obtained by using genetic algorithms (GAs). Different weight vectors provide different solutions.
5. Numerical results 5.1. Example 1 For the purpose of testing the Gaussian function methodology, the problem of a square metal structure protected by one sacrificial anode, shown in Fig. 1, was initially solved. First, it is considered the BVP (1)–(3) making ϕ0 = − 1280 mV , where this value was chosen with the goal of providing a potential distribution over the metal surface below the critical potential (ϕ ≤ ϕc = − 850 mV). Second, the Gaussian function is used to model the sacrificial anode originally represented by boundary Γ2 in Fig. 1. After solving original BVP (1)–(3), Eq. (12) can be used to determine the Gaussian function parameter P. Only in this example, the parameter nn is estimated of empirical form, i.e., some values for nn are tested in order to ϕint ≈ − 1280 mV = ϕ0 . The function that represents the cathodic polarization curve is given by Eq. (4) The BVP (1)–(3) is solved using a standard MFS and the LMDIF routine. The boundary is represented with 156 boundary nodes in Γ1 and 14 boundary nodes in Γ2 have been adopted. The initial solution to the coefficients cj 's is considered equal to (10.0, − 10.0, … , 10.0, − 10.0). The virtual sources are created using two circles: the first circle that is associated to the boundary Γ1 has 154 virtual sources and it is located outside the domain and the cavity, whereas the second circle with 14 virtual sources is associated to the boundary Γ2 and therefore will be within the cavity. The initial radii chosen were ρ1 = 70 cm for the first virtual circle and ρ2 = 2.5 cm for the second virtual circle. The centre of the virtual circles is (50 cm, 50 cm). The optimum radii determined by routine LMDIF after 11,220 function evaluations (NFEV) are ρ1* = 145.464 cm and ρ2* = 0.992 cm . Fig. 2 presents the potential distribution in the electrolyte using MFS for the original BVP (1)–(3). Now, the corresponding BVP is solved using the Gauss function, Eq. (9). The parameter P (50, 50) was calculated using Eq. (12) and a simple Gauss quadrature rules with four integration points. The flux integral in Γ2 obtained was P (50, 50) = − 47, 040.226 μA and after empirical tests it was chosen nn¼ 0.14. The calculations were performed for 156 boundary nodes (Γ1) and 155 virtual sources with initial radius ρ1 = 70 cm . The value for support parameter β is equal to 40.0 cm . The optimum radius determined by routine LMDIF with NFEV ¼9024 is ρ1⁎ = 81.922 cm . Fig. 3 shows the potential distribution obtained on the metal surface considering original BVP and the corresponding BVP. The potential distribution in the electrolyte considering Gauss function like sacrificial anode is presented in Fig. 4. The value found for the ϕint was −1279.1 mV( ≈ − 1280 mV). As seen in Figs. 2–4, satisfactory coincidence of results between the original problem and the corresponding BVP using Gauss function has been obtained.
Fig. 2. Potential in the electrolyte for the original BVP.
Fig. 3. Potential distribution on the metal.
5.2. Example 2 In order to assess the efficiency of the multiobjective optimization, example 1 is solved again considering this time the sacrificial anode initially randomly distributed within the electrolyte. The purpose of this optimization is to search the following design variables: the optimum location of the sacrificial anode (x0, y0 ), the parameter nn and the height of the Gaussian function P. This purpose will be reached minimizing the objective function Z by
Fig. 4. Potential in the electrolyte for the corresponding BVP.
39
40
W.J. Santos et al. / Engineering Analysis with Boundary Elements 73 (2016) 35–41
Fig. 5. Potential distribution on the metal (corresponding BVP).
genetic algorithm. The parameters x0, y0, nn and P have ranges equal to [10.0 cm, 90.0 cm], [10.0 cm, 90.0 cm], [0.01, 0.2] and [ − 60, 000.0 μA, − 30, 000.0 μA], respectively. In addition, the following GA values have been used: population size =30, maximum number generations =100, required precision =0.01, initial crossover probability =0.6, final crossover probability =0.5, initial mutation probability =0.035, final mutation probability =0.055 and the penalty number v = 103. The weights of the multiple objective function Z are w1 ¼0.5 and w2 ¼ 0.5. The optimum location determined by GA is (49.498 cm, 50.101 cm) and the optimum parameters for nn and P are 0.159 and −47, 197.896 μA , respectively. As can be seen in Figs. 5 and 6, the surface of the structure is properly protected after optimization (ϕ ≤ ϕc = − 850 mV) and the anode is quite close to the centre of the square as expected. 5.3. Example 3 In this example a two-dimensional galvanic field formed by a rectangular low carbon steel structure and two aluminium alloy anodes in electrolyte with conductivity k = 0.0479 Ω−1 cm−1 is analysed. The boundary was represented with 134 points
uniformly distributed and 133 virtual sources. Galvanic anodes were treated as having a constant effective potential based on the material of which they were constructed. This method is strictly valid at low current densities. The fixed potential for the aluminium alloy is ϕ0 = − 1055.0 mV (vs. SCE). The polarization curve is given by Eq. (4) and the dimensions of the structure are 100 cm × 50 cm . The starting position of each anode have coordinates x0 and y0 with ranges equal to [10.0 cm, 90.0 cm] and [10.0 cm, 40.0 cm], respectively. The ranges for the Gaussian functions parameters nn and P are assumed as [0.01, 0.2] and [ − 30, 000.0 μA, − 10, 000.0 μA], respectively. Now, the goal is to provide a potential distribution on the metal within the range −950 mV ≤ ϕ ≤ − 850 mV (pre-established protection criterion) and to guarantee the equivalence between the original BVP and the corresponding BVP using the objective function z2, Eq. (29). The entrance parameter values for the GA are the same presented in the previous example. However, the weights of the multiple objective function Z are w1 ¼ 0.3 and w2 ¼0.7. In this example, two sacrificial anodes are considered for the cathodic protection. For this reason, a potential distribution on the metal within the pre-established protection criterion can be reached more easily when compared with the previous examples. The objective function z1 establishes the protection criterion of metallic structures, whereas the objective function z2 guarantees the equivalence between the standard boundary value problem (BVP) and the corresponding BVP. Hence, in order to improve the results, it was necessary to increase the value of weighting parameter of the objective function z2. Fig. 7 shows the potential distribution obtained on the metal surface after multiobjective optimization, where the boundary node 1 has coordinates (0, 0). This figure presents the results considering the original BVP and the corresponding BVP in order to analyse the pre-established protection criterion and the equivalence between the problems. The optimum location of the anodes determined by GA after 91 generations is shown in Fig. 8 (corresponding BVP) and Fig. 9 (original BVP), where the similarity of the results in the electrolyte can also be seen.
6. Concluding remarks The results here presented showed that the proposed strategy is effective in the optimization of the anode position in a galvanic cathodic protection system using GA and MFS. Furthermore, the constant potential circular anodes were successfully represented by a Gaussian function whose parameters were reached using a multiobjective optimization methodology. The dual reciprocity method (DRM) was implemented using a
Fig. 6. Potential in the electrolyte after GA optimization (corresponding BVP).
Fig. 7. Potential distribution on the metal.
W.J. Santos et al. / Engineering Analysis with Boundary Elements 73 (2016) 35–41
41
References
Fig. 8. Potential in the electrolyte after GA optimization (corresponding BVP).
Fig. 9. Potential in the electrolyte after GA optimization (original BVP).
positive definite compactly supported radial basis function (CSRBF) of Wendland. The Wendland CS-RBF guaranteed a good approximation for the source term (Gauss function) and the solvability of the linear system used to solve a particular solution of the Poisson equation. The examples presented show satisfactory equivalence between the original problem (sacrificial anode CP system) and the corresponding BVP, where the sacrificial anodes were introduced as source terms using the probability density function (Gaussian). Thus, it was not necessary to redefine the boundary geometric nodes and the arrangement of virtual sources for the standard method of fundamental solutions (MFS) formulation in the process of optimal positioning of sacrificial anodes. The adopted MINPACK routine LMDIF showed satisfaction to solve the nonlinear least squares problem that is a result of the nonlinear boundary conditions in the context of the MFS formulation. However, a large number of function evaluations (NFEV) were required for the convergence of the minimization process. This is due to the large number of design variables of the nonlinear least squares problem. In general, increasing the number of design variables slowed down the convergence of the minimization process.
Acknowledgements The research was sponsored by CNPq-Brazil.
[1] Fontana MG, Greene ND. Corrosion engineering.New York: McGraw-Hill; 1967. [2] Roberge PR. Handbook of corrosion engineering.New York: McGraw-Hill; 1999. [3] Santos WJ, Santiago JAF, Telles JCF. An application of genetic algorithms and the method of fundamental solutions to simulate cathodic protection systems. Comput Model Eng Sci 2012;87:23–40. [4] Santos WJ, Santiago JAF, Telles JCF. Optimal positioning of anodes and virtual sources in the design of cathodic protection systems using the method of fundamental solutions. Eng Anal Bound Elem 2014;46:67–74. [5] Montoya R, Rendon O, Genesca J. Mathematical simulation of cathodic protection system by finite element method. Mater Corros 2005;56:404–11. [6] Montoya R, Aperador W, Bastidas DM. Influence of conductivity on cathodic protection of reinforced alkali-activated slag mortar using the finite element method. Corros Sci 2009;51:2857–62. [7] Montoya R, Galvan JC, Genesca J. Using the right side of Poisson's equation to save on numerical calculations in FEM simulation of electrochemical systems. Corros Sci 2011;53:1806–12. [8] Miyasaka M, Hashimoto K, Aoki S. A boundary element analysis on galvanic corrosion problems - computational accuracy on galvanic fields with screen plates. Corros Sci 1990;30:299–311. [9] Abootalebi O, Kermanpur A, Shishesaz MR, Golozar MA. Optimizing the electrode position in sacrificial anode cathodic protection systems using boundary element method. Corros Sci 2010;52:678–87. [10] Telles JCF, Mansur WJ, Wrobel LC, Marinho MG. Numerical simulation of a cathodically protected semisubmersible platform using PROCAT system. Corrosion 1990;46:513–8. [11] Santiago JAF, Telles JCF. On boundary elements for simulation of cathodic protection systems with dynamic polarization curves. Int J Numer Methods Eng 1997;40:2611–22. [12] BEASY User Guide Beasy. Computational mechanics. Ashurst, Southampton, UK: BEASY Ltd.; 2000. [13] Azevedo JPS, Wrobel LC. Nonlinear heat conduction in composite bodies: a boundary element formulation. Int J Numer Methods Eng 1988;26:19–38. [14] Aoki S, Amaya K. Optimization of cathodic protection system by BEM. Eng Anal Bound Elem 1997;19:147–56. [15] Wrobel LC, Miltiadou P. Genetic algorithms for inverse cathodic protection problems. Eng Anal Bound Elem 2004;28:267–77. [16] Kupradze VD, Aleksidze MA. Aproximate method of solving certain boundaryvalue problems. Soobshch Akad Nauk Gruz SSR 1963;30:529–36. [17] Costa EGA, Godinho L, Santiago JAF, Pereira A, Dors C. Efficient numerical models for the prediction of acoustic wave propagation in the vicinity of a wedge coastal region. Eng Anal Bound Elem 2011;35:855–67. [18] Fontes Jr. EF, Santiago JAF, Telles JCF. On a regularized method of fundamental solutions coupled with the numerical Green's function procedure to solve embedded crack problems. Eng Anal Bound Elem 2013;37:1–7. [19] Karageorghis A, Fairweather G. The method of fundamental solutions for the solution of nonlinear plane potential problems. IMA J Numer Anal 1989;9:231–42. [20] Garbow BS, Hillstrom KE, Mor JJ. MINPACK project. Argonne National Laboratory, Argonne, Illinois, USA; 1980. [21] More L. Levenberg–Marquardt algorithm: implementation and theory. In: Proceedings of the biennial conference held at Dundee, June 28–July 1. Argonne, IL, USA: Argonne National Laboratory; 1977. [22] Partridge P, Brebbia CA, Wrobel LC. The dual reciprocity boundary element method.Southampton: Computational Mechanics Publications; 1992. [23] Yan JF, Pakalapati SNR, Nguyen TV, White RE. Mathematical modelling of cathodic protection using the boundary element method with nonlinear polarisation curves. J Electrochem Soc 1992;139:1932–6. [24] Atkinson KE. The numerical evaluation of particular solution for Poisson's equation. IMA J Numer Anal 1985;5:319–38. [25] Golberg MA. The numerical evaluation of particular solutions in the BEM - a review. Bound Elem Commun 1995;6:99–106. [26] Wendland H. Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree. Adv Comput Math 1995;4:389–96. [27] Golberg MA. The method of fundamental solutions for Poisson's equation. Eng Anal Bound Elem 1995;16:205–13. [28] Chen CS, Golberg MA, Schaback RA. Recent developments of the dual reciprocity method using compactly supported radial basis functions. In: Transformation of domain effects to the boundary. Southampton, Boston: WIT Press; 2003. p. 183–225. [29] Fletcher R. Practical methods of optimisation.Chichester, UK: John Wiley & Sons; 1987. [30] Srinivas N, Deb K. Multiobjective optimization using nondominated sorting in genetic algorithms. Evol Comput 1994;2:221–48. [31] Michalewicz Z. Genetic algorithms þ data structures¼ evolution programs. Springer-Verlag, New York: Springer-Verlag, Berlin Heidelberg; 1996. [32] Mitchell M. An introduction to genetic algorithms.Cambridge, MA: The MIT Press; 1997.
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
Using the Gaussian function to simulate constant potential anodes in multiobjective optimization of cathodic protection systems W.J. Santos, J.A.F. Santiago n, J.C.F. Telles Department of Civil Engineering, COPPE/UFRJ, Caixa Postal 68506, CEP21941-972, Rio de Janeiro, RJ, Brazil
art ic l e i nf o
a b s t r a c t
Article history: Received 22 March 2016 Received in revised form 7 July 2016 Accepted 31 August 2016
The purpose of this work is to numerically find the optimum location of constant potential anodes to ensure complete structure surface protection using a cathodic protection technique. The existence of sacrificial anodes is originally introduced through the boundary conditions of the corresponding boundary value problem (BVP). However, if constant potential galvanic regions are introduced through its boundaries, then finding their optimal location is not an easy task due to the necessity of redefining boundary geometric nodes and the arrangement of virtual sources for the standard method of fundamental solutions (MFS) formulation. Therefore, in this work, the galvanic anodes are introduced as source terms using a Gaussian function. Hence, the boundary remains the same for different anode positions. The optimization process includes the identification of the following parameters characterizing the Gaussian function: the optimum coordinates of the centre of the anode, a factor that involves the inherent potential of the electrode and a proportionality factor for the electrode diameter. The MFS methodology coupled with a genetic algorithm presented good results for this multiobjective optimization procedure. This fact can be seen in the several results of applications that are discussed in this paper, considering numerical simulations in finite regions in R2. & 2016 Elsevier Ltd. All rights reserved.
Keywords: Genetic algorithms Method of fundamental solutions Dual reciprocity
1. Introduction Cathodic protection (CP) is a technique to reduce the corrosion rate of the metal surface by making it the cathode of an electrochemical cell. In practice, the goal is to provide a uniform potential distribution on the metal surfaces, limited by a minimum potential value to guarantee protection from the corrosion using an external power source (impressed current CP) or by utilizing a sacrificial anode [1,2]. In this work, it is considered the case of sacrificial anode system, where the galvanic relationship between a sacrificial anode material, such as zinc or magnesium, and the structure is used to supply the required CP current. The most commonly used methods for modelling cathodic protection systems are finite element method (FEM) and boundary element method (BEM). However, the method of fundamental solutions (MFS) is a technique which can also be applied to CP problems (see references [3,4]). Just like BEM, MFS is applicable when a fundamental solution of the differential equation in question is known, with the advantage of not requiring any integration procedure or specific treatment for the singularities of n
Corresponding author. E-mail addresses: [email protected] (W.J. Santos), [email protected] (J.A.F. Santiago), [email protected] (J.C.F. Telles). http://dx.doi.org/10.1016/j.enganabound.2016.08.014 0955-7997/& 2016 Elsevier Ltd. All rights reserved.
the fundamental solution. The BEM and the MFS are the most appropriate techniques to solve problems involving galvanic corrosion and CP systems, mainly to solve large problems and considering homogeneous conductive medium. These methods require only the representation of anodes and cathodes surfaces, which leads to better resolution and reduction in computer run time when compared to FEM. Mathematical simulations of cathodic protections systems using FEM can be seen in [5–7], where compact support functions were used to simulate constant potential electrodes in electrochemical process. Due to its accuracy and simplicity of mesh generation, the BEM is usually used for numerical simulations of sacrificial anode CP systems. Miyasaka et al. [8], for example, evaluated the computational accuracy of BEM to estimate the galvanic corrosion and CP in an actual field. Abootalebi et al. [9] determined the optimum location of zinc anode electrode using BEM in 2D. In addition, the influence of anode length and paint defect on corrosion current density and potential distributions of sacrificial anode CP system were investigated. Several different applications of BEM to study CP systems have been reported in the literature, including reference to practical analyses performed by offshore oil companies [10–12]. The BEM implementation includes a Newton–Raphson solution algorithm to accommodate possible nonlinear boundary conditions [13]. Coupled with the numerical method to solve the Laplace equation, optimization algorithms can be used to
36
W.J. Santos et al. / Engineering Analysis with Boundary Elements 73 (2016) 35–41
determine the optimum location and the corresponding current intensity values of the anodes in order to satisfy a protection criterion. Hence, the minimization of an objective function using, for example, genetic algorithms (GAs) and a penalty method for handling constraints can be adopted. This type of optimization can and has also been successfully performed using BEM [14,15]. Kupradze and Aleksidze [16] first proposed the basic ideas for the formulation of the MFS. In order to construct the solution, the MFS uses only a superposition of fundamental solutions associated to the problem, with singular points (virtual sources) located outside the domain. The accuracy of the MFS numerical solution depends on the radius of such a circle or on the distance from the virtual sources over the geometrically similar boundary contour to the problem boundary, especially due to possible ill-conditioning and/or rank-deficiency of the algebraic system of equations formed. The MFS has successfully been applied for solving several problems. For example, Costa et al. [17] developed numerical frequency domain formulations to simulate the 2D acoustic wave propagation in the vicinity of an underwater configuration which combines two sub-regions using the MFS. More recently, Fontes Jr. et al. [18] applied a regularized method of fundamental solutions coupled with the numerical Green's function procedure to solve linear elastic fracture mechanics (LEFM) problems. Problems with nonlinear boundary conditions solved by MFS can be treated as nonlinear least squares problems [19]. Therefore, to determine the coefficients of the linear superposition of the fundamental solutions and the positions of the virtual sources, a nonlinear least squares algorithm is found necessary. In the present work, the minimization of the nonlinear functional is done using the MINPACK [20] routine LMDIF, which is a modified version of the Levenberg–Marquardt algorithm [21]. The Levenberg–Marquardt method has been highly recommended when Jacobian is rankdeficiency or nearly so. Santos et al. [3] were the first to use standard MFS successfully in the numerical simulations of CP systems. In the paper cited, the authors proposed a GA with the MFS to simulate cathodic protection systems with nonlinear boundary conditions. The adopted GA was used to minimize a nonlinear error function, whose design variables were the coefficients of the linear superposition of fundamental solutions and the positions of the virtual sources, randomly distributed outside the problem domain. In Santos et al. [4] was presented a formulation using a GA and the MFS to determine the optimum location and the optimum current intensity of the anodes inserted in the electrolyte leading to a practical optimized design procedure. The results presented in this paper included a comparison with a direct boundary element (BEM) solution procedure. The purpose of the present paper is to use a genetic algorithm (GA) with the MFS for optimizing the anode position in a galvanic cathodic protection system. Furthermore, the constant potential circular anodes are here mathematically represented by a Gaussian function. Thus, the GA will be used to search the Gaussian function parameters: coordinates of the centre of the anode, a factor that involves the inherent potential of the anode and a proportionality factor for the anode diameter. The main advantage of considering constant potential regions as the source term is that the boundary conditions remain the same, no matter where the anodes are located. The circular anodes are considered and, therefore, a good approximation for the source term is a function with circular compact support like the Gaussian function. In the MFS case, it is not necessary for the arrangement of the virtual sources inside the galvanic anodes (cavities), when they are introduced as source terms. A particular solution to the inhomogeneous differential equation resulting will be evaluated using the dual reciprocity method (DRM). In the DRM, the source term is approximated by a finite series of radial basis functions with an approximation to
particular solution calculated analytically from source [22]. The corresponding homogeneous solution is found by MFS. This text is organized as follows: in Section 2 is presented a standard electrochemical potential problem, i.e., the Laplace equation with boundary conditions given by constant potential (anode) and polarization curve (cathode). The decision of introducing the Gaussian function to model constant potential galvanic regions is also discussed in Section 2. The standard MFS formulation with the dual reciprocity method for solving the Poisson equation with nonlinear boundary conditions is discussed in Section 3. The multiobjective optimization problem and the genetic algorithm are showed in Section 4. This algorithm is coupled with MFS to guarantee protection from the corrosion by utilizing a sacrificial anode in the numerical examples solved in Section 5. The paper ends with some discussions found in Section 6.
2. Boundary value problem (BVP) For a homogeneous and isotropic electrolyte (domain Ω) system of conductivity k, as illustrated in Fig. 1, the electrochemical potential problem studied obeys the Laplace equation:
k∇2ϕ(x) = 0,
x ∈ Ω,
(1)
subjected to the following boundary conditions:
i(x) = F (ϕ),
ϕ(x) = ϕ0,
x ∈ Γ1,
x ∈ Γ2,
(2)
(3)
where Γ ≡ Γ1 ∪ Γ2 is the boundary of Ω, ϕ0 is the fixed potential of anode, i(x) is the current density in the outward normal direction n and F (ϕ) is a nonlinear function of ϕ given here by the following polarization curve [23]:
i = F (ϕ) = e
ϕ + 693.91 β1
−1 ϕ + 521.6 ⎤ ϕ + 707.57 ⎡1 − β3 − ⎢ + e β2 ⎥ − e , ⎣ i1 ⎦
(4)
with ϕ and i having units mV and μA/cm2, respectively, and β1, β2, β3 and i1 are given constant parameters: β1 = 24 mV ,
Fig. 1. The original potential problem.
W.J. Santos et al. / Engineering Analysis with Boundary Elements 73 (2016) 35–41
β2 = 23.47 mV , β3 = 55 mV e i1 = 86.06 μA/cm2 . The conductivity of the electrolyte is equal to k = 0.0479 Ω−1 cm−1 and the critical value of the electrochemical potential is ϕc = − 850 mV (vs. SCE). The boundary value problem (1)–(3) can be modified in order to simplify the optimization of positioning of anodes for cathodic protection. For this, it had been decided to consider the constant potential condition ϕ0 for Γ2 on the right-hand side of Poisson equation and not in the boundary condition. Thus, the corresponding BVP can be written in the following form:
k∇2ϕ(x) = b(x), x ∈ Ω,
(5)
with boundary condition given by Eq. (2). The term b presented in Eq. (5) is a function representing the anode as external source. Here, the source term b(x) is mathematically represented by a Gaussian function, defined in two-dimensional space (x,y) by
b(x, y) = Ae−nn(x − x0)
2
− nn(y − y0 )2
,
(6)
where A is height of the Gaussian function, the coordinates x0 and y0 are the centre of this function and nn is a proportionality factor that controls the width of the Gaussian. The volume (V ) under the Gaussian function is given by
V=
Aπ , ∫Ω b(x, y)dΩ = nn
The integral in Eq. (7) is 1 if A = and, for this case, the Gaussian function is a probability density function of a normal distribution that can be written as
nn −nn(x − x0)2 − nn(y − y0 )2 b(x, y) = e . π
(8)
In Eq. (8), if the parameter nn tends to infinity, then the function b(x, y ) tends to Dirac δ “function”. Treating the anode as Gaussian function, the term b(x) is equal to expression (8) multiplied by density of the source P (x0, y0 ). Therefore:
b(x) = b(x, y) = P (x 0 , y0 )
nn −nn(x − x0)2 − nn(y − y0 )2 e . π
(9)
The conservation of current between anode and cathode in the original problem (see Fig. 1) is automatically obtained by satisfying the condition
∫Γ i(x)dΓ = − ∫Γ i(x)dΓ . 1
2
(10)
On the other hand, using Eqs. (2) and (5), the corresponding BVP satisfies the following condition:
∫Γ i(x)dΓ = − P(x0, y0 ), 1
(11)
which ensures that there is no loss of current. To impose equivalence between the problems, i.e., between Eqs. (10) and (11), it is necessary to guarantee that
P (x 0 , y0 ) =
∫Γ i(x)dΓ . 2
ϕ(x) = ϕp(x) + ϕh(x),
(13)
where ϕp(x) satisfies the inhomogeneous equation (5) but does not necessarily satisfy the nonlinear boundary condition presented in Eq. (2) and ϕh(x) satisfies the Laplace equation and the corresponding boundary condition for the homogeneous problem, which in the nonlinear case, is given by
ih(x) = F (ϕp + ϕh) − ip(x),
x ∈ Γ1,
(14)
ϕp is a particular solution and ip(x) is its normal derivative. ∂ϕ Once ϕp and p are known, ϕh(x) can be determined by standard ∂n where
MFS. The determination of the particular solution depends on the type of source b(x) in (5). In the simple cases, as point sources for example, ϕp(x) can be evaluated analytically. For general cases, the particular solutions can be determined by numerical techniques [24,25]. Here, a particular solution will be derived using the context of the dual reciprocity method using a positive definite compactly supported radial basis function (CS-RBF) of Wendland [26]. Thus, the function b(x) is approximated by a finite series of radial basis functions ψ (r ) with an approximation ϕ^ to ϕ calcup
k∇2ϕ^p(x) =
N+L
N+L
∑ ψj(x)aj = ∑ ψj( j=1
x − x j )aj , (15)
j=1
where N are the boundary nodes, L are the internal nodes and the unknown coefficients aj 's are determined solving the system of equations N+L
∑ ψj(x)aj = b^(x),
(16)
j=1
^ with the point values of the known function b denoted for b . For practical implementation, the Wendland radial basis function of compact support is scaled with scale β as follows:
⎛r⎞ ψ ( x − x j ) = ψ ⎜ ⎟. ⎝ β⎠
(17)
The system (16) is symmetric and can be dense or sparse depending on the support parameter β. The accuracy and efficiency of the Wendland CS-RBF depends on the scale of the support and determining the scale of support can be uncertain. For this reason, the singular value decomposition (SVD) procedure is used to provide acceptable solutions to possible ill-conditioned system of equations generated by the Wendland CS-RBF. After determining the coefficients aj 's using (16) an approximation to a particular solution can be written as
ϕ^p(x) =
N+L
∑ Ψj(x)aj,
(18)
j=1
(12)
Eq. (12) shows that the determination of P depends not only on the potential at which the anode is set, but also on the geometry of the domain and the location of coordinates x0 and y0. The diameter of the anode is related to the parameter nn in Eq. (9).
3. Method of fundamental solutions for poisson's equation In order to solve the corresponding BVP, Eqs. (2) and (5), using MFS, it is necessary to eliminate the inhomogeneous term in (5) by the use of a particular solution. Let
p
lated analytically from ψ [27]. Mathematically, the following differential equation is solved:
(7) nn π
37
where
Ψj is the solution of
2
k∇ Ψj(x) = ψj(x).
(19)
The solution of Eq. (19) can easily be achieved considering the radial part of Laplacian operator and by analytical integration of Eq. (19). In this paper, the following Wendland radial basis function is chosen: 6⎛ ⎛ ⎞ ⎛ ⎞ r r⎞ r2 r ψ ⎜⎜ ⎟⎟ = ⎜ 1 − ⎟ ⎜⎜ 35 2 + 18 + 3⎟⎟. β⎠ ⎝ β β ⎝ β⎠ ⎝ ⎠
An explicit analytical representation of
(20) r Ψ(β)
corresponding to
38
W.J. Santos et al. / Engineering Analysis with Boundary Elements 73 (2016) 35–41
Wendland's CS-RBF presented in Eq. (20) can be derived by straightforward integration [28]:
4. Multiobjective optimization using GA
r Ψ( ) β =
1 k ⎧ 7r10 64r 9 105r 8 64r 7 35r 6 7r 4 3r 2 ⎪ − + − + − + , 2 8 7 6 5 4 4 ⎪ 20β 4β 27β 16β 7β 6β ⎪ ⎨ r ≤ β, ⎪ ⎪ 3517β 2 β2 ⎛ r ⎞ ⎪ 15120 + 6 ln⎜ β ⎟, r > β. ⎝ ⎠ ⎩
(21)
After obtaining Ψj by Eq. (19), it is possible to determine a particular solution using Eq. (18) and, consequently, the approximate solution of the associated homogeneous problem by MFS. The method of fundamental solutions establishes that the approximate solution can be constructed by a summation of similar problem solution given by the following superposition: nvs
ϕh(x) =
∑ G(x,
xvs j )cj,
j=1
(22)
with nvs being the total number of virtual sources, xvs j represents the position of the virtual sources and the coefficients cj 's are the unknown source intensities. The function G(x, xvs j ) is a fundamental solution of Laplace's equation given by
G(x,
xvs j )
1 1 = ln , 2π k R
(23)
where R is the Euclidean distance between point xvs and the field point x . ∂G Similarly, defining H = k ∂n , the homogeneous solution for the current density (ih) is given as nvs
ih(x) =
user-specified maximum number of function evaluations is reached.
∑ H (x, xvsj )cj. j=1
(24)
The idea of MFS is to determine the coefficients cj 's by imposing satisfaction of the boundary conditions at certain boundary nodes. The boundary condition given by Eq. (14) generates a nonlinear least squares problem in which the design variables are the coefficients cj 's and the positions of the virtual sources. In least squares problems, the objective function f has the following special form:
1 f (c, xvs) = 2
nbn
∑ ϵ2j (c, xvs), j=1
(25)
where nbn is the number of boundary nodes, c is a vector containing the coefficients cj 's and each ϵj is a smooth function referred to as a residual given by
ϵj = igj − F (ϕgj ),
j = 1, …, nbn .
(26)
Considering the arrangement of the virtual sources on a circular contour, it is only necessary to search for the radius ρ of such a circle. Thus, design variables of Eq. (25) are the coefficients cj 's and a radius ρ. In this paper, the minimization of the functional of Eq. (25) is carried out using the MINPACK routine LMDIF, which is a modified version of the Levenberg–Marquardt algorithm. In LMDIF, the Jacobian is evaluated internally by finite differences. LMDIF terminates when either a user-specified tolerance is achieved or the
In this work, there are two optimization goals: (1) to pre-establish the protection criterion of metallic structures and (2) to guarantee the equivalence between the original BVP (1)–(3) and the corresponding BVP, Eqs. (2) and (5). In other words, it is necessary to search the parameters of Gaussian function in Eq. (9) in order to satisfy this multiobjective optimization. The aim of CP systems is to reduce the potential over the structure surface below a critical value (ϕc ). By optimizing the sacrificial anode (Γ2) location, the maximum corrosion current over the structure surface will be minimized in such a way as to satisfy, as uniformly as possible, close to the critical potential value: ϕ ≤ ϕc . The optimum anode location can be achieved by minimizing the proposed objective function:
z1(x 0 , y0 ) =
1 nbn
nbn
2
∑ ⎡⎣ ϕi − ϕc⎤⎦
+v
i=1
1 * nbn
nbn
∑ kl2(x 0, y0 ) , i=1
(27)
where z1 calculates the root mean square error (RMSE) between the electrochemical potential at each boundary node and the cri* is the number of boundary nodes tical potential. Furthermore, nbn that do not satisfy the protection criterion, the constant v is a penalty number and function kl is equal to:
kl(x 0 , y0 ) = (ϕi − ϕc )u(ϕi − ϕc ),
(28)
where u is the unit step function. Typical values of v for the penalty method are within the range 102–105 [29]. In addition, it is also necessary to search the optimum solution for the Gaussian function parameters nn and P (x0, y0 ). The diameter of the anode is related to the parameter nn and the height of the Gaussian function P is related to the constant potential condition, but also on the geometry of the domain and on the location of coordinates (x0, y0 ). The boundary nodes in Γ2 are internal nodes in corresponding BVP and the potential value in these nodes must be the anode prescribed potential ϕ0. Thus, the second objective function z2 is modelled in order that the internal potential values corresponding to Γ2 obey the Dirichlet condition of the original BVP (1)–(3). For this, it was decided to minimize the following objective function:
z2(nn, P ) =
1 nint
nint
2
∑ ⎡⎣ ϕinti − ϕ0⎤⎦
,
(29)
i=1
where nint is the number of internal nodes that represent the boundary of the anode in the original problem and ϕint are the potential values at these internal nodes solved for each pair (nn,P). In this optimization problem there are two objectives functions to be optimized: z1 and z2. There are some classical methods for multiobjective optimization [30]. Here, a method of objective weighting is applied, where multiple objective functions zi are combined into one overall objective function Z:
Z (x 0 , y0 , nn, P ) = w1z1 + w2z2,
(30) 2 ∑i = 1 wi
where the weights wi ∈ [0, 1] and = 1. In this paper, the minimization of Eq. (30) is achieved using a typical genetic algorithm (GA) with binary representation [31]. However, some characteristics have been included as the twopoint crossover, the elitism and the probabilities of mutation and crossover can vary linearly over the generations [32]. The weighting parameters w1 and w2 are chosen through the results of
W.J. Santos et al. / Engineering Analysis with Boundary Elements 73 (2016) 35–41
the multiobjective optimization process, which are obtained by using genetic algorithms (GAs). Different weight vectors provide different solutions.
5. Numerical results 5.1. Example 1 For the purpose of testing the Gaussian function methodology, the problem of a square metal structure protected by one sacrificial anode, shown in Fig. 1, was initially solved. First, it is considered the BVP (1)–(3) making ϕ0 = − 1280 mV , where this value was chosen with the goal of providing a potential distribution over the metal surface below the critical potential (ϕ ≤ ϕc = − 850 mV). Second, the Gaussian function is used to model the sacrificial anode originally represented by boundary Γ2 in Fig. 1. After solving original BVP (1)–(3), Eq. (12) can be used to determine the Gaussian function parameter P. Only in this example, the parameter nn is estimated of empirical form, i.e., some values for nn are tested in order to ϕint ≈ − 1280 mV = ϕ0 . The function that represents the cathodic polarization curve is given by Eq. (4) The BVP (1)–(3) is solved using a standard MFS and the LMDIF routine. The boundary is represented with 156 boundary nodes in Γ1 and 14 boundary nodes in Γ2 have been adopted. The initial solution to the coefficients cj 's is considered equal to (10.0, − 10.0, … , 10.0, − 10.0). The virtual sources are created using two circles: the first circle that is associated to the boundary Γ1 has 154 virtual sources and it is located outside the domain and the cavity, whereas the second circle with 14 virtual sources is associated to the boundary Γ2 and therefore will be within the cavity. The initial radii chosen were ρ1 = 70 cm for the first virtual circle and ρ2 = 2.5 cm for the second virtual circle. The centre of the virtual circles is (50 cm, 50 cm). The optimum radii determined by routine LMDIF after 11,220 function evaluations (NFEV) are ρ1* = 145.464 cm and ρ2* = 0.992 cm . Fig. 2 presents the potential distribution in the electrolyte using MFS for the original BVP (1)–(3). Now, the corresponding BVP is solved using the Gauss function, Eq. (9). The parameter P (50, 50) was calculated using Eq. (12) and a simple Gauss quadrature rules with four integration points. The flux integral in Γ2 obtained was P (50, 50) = − 47, 040.226 μA and after empirical tests it was chosen nn¼ 0.14. The calculations were performed for 156 boundary nodes (Γ1) and 155 virtual sources with initial radius ρ1 = 70 cm . The value for support parameter β is equal to 40.0 cm . The optimum radius determined by routine LMDIF with NFEV ¼9024 is ρ1⁎ = 81.922 cm . Fig. 3 shows the potential distribution obtained on the metal surface considering original BVP and the corresponding BVP. The potential distribution in the electrolyte considering Gauss function like sacrificial anode is presented in Fig. 4. The value found for the ϕint was −1279.1 mV( ≈ − 1280 mV). As seen in Figs. 2–4, satisfactory coincidence of results between the original problem and the corresponding BVP using Gauss function has been obtained.
Fig. 2. Potential in the electrolyte for the original BVP.
Fig. 3. Potential distribution on the metal.
5.2. Example 2 In order to assess the efficiency of the multiobjective optimization, example 1 is solved again considering this time the sacrificial anode initially randomly distributed within the electrolyte. The purpose of this optimization is to search the following design variables: the optimum location of the sacrificial anode (x0, y0 ), the parameter nn and the height of the Gaussian function P. This purpose will be reached minimizing the objective function Z by
Fig. 4. Potential in the electrolyte for the corresponding BVP.
39
40
W.J. Santos et al. / Engineering Analysis with Boundary Elements 73 (2016) 35–41
Fig. 5. Potential distribution on the metal (corresponding BVP).
genetic algorithm. The parameters x0, y0, nn and P have ranges equal to [10.0 cm, 90.0 cm], [10.0 cm, 90.0 cm], [0.01, 0.2] and [ − 60, 000.0 μA, − 30, 000.0 μA], respectively. In addition, the following GA values have been used: population size =30, maximum number generations =100, required precision =0.01, initial crossover probability =0.6, final crossover probability =0.5, initial mutation probability =0.035, final mutation probability =0.055 and the penalty number v = 103. The weights of the multiple objective function Z are w1 ¼0.5 and w2 ¼ 0.5. The optimum location determined by GA is (49.498 cm, 50.101 cm) and the optimum parameters for nn and P are 0.159 and −47, 197.896 μA , respectively. As can be seen in Figs. 5 and 6, the surface of the structure is properly protected after optimization (ϕ ≤ ϕc = − 850 mV) and the anode is quite close to the centre of the square as expected. 5.3. Example 3 In this example a two-dimensional galvanic field formed by a rectangular low carbon steel structure and two aluminium alloy anodes in electrolyte with conductivity k = 0.0479 Ω−1 cm−1 is analysed. The boundary was represented with 134 points
uniformly distributed and 133 virtual sources. Galvanic anodes were treated as having a constant effective potential based on the material of which they were constructed. This method is strictly valid at low current densities. The fixed potential for the aluminium alloy is ϕ0 = − 1055.0 mV (vs. SCE). The polarization curve is given by Eq. (4) and the dimensions of the structure are 100 cm × 50 cm . The starting position of each anode have coordinates x0 and y0 with ranges equal to [10.0 cm, 90.0 cm] and [10.0 cm, 40.0 cm], respectively. The ranges for the Gaussian functions parameters nn and P are assumed as [0.01, 0.2] and [ − 30, 000.0 μA, − 10, 000.0 μA], respectively. Now, the goal is to provide a potential distribution on the metal within the range −950 mV ≤ ϕ ≤ − 850 mV (pre-established protection criterion) and to guarantee the equivalence between the original BVP and the corresponding BVP using the objective function z2, Eq. (29). The entrance parameter values for the GA are the same presented in the previous example. However, the weights of the multiple objective function Z are w1 ¼ 0.3 and w2 ¼0.7. In this example, two sacrificial anodes are considered for the cathodic protection. For this reason, a potential distribution on the metal within the pre-established protection criterion can be reached more easily when compared with the previous examples. The objective function z1 establishes the protection criterion of metallic structures, whereas the objective function z2 guarantees the equivalence between the standard boundary value problem (BVP) and the corresponding BVP. Hence, in order to improve the results, it was necessary to increase the value of weighting parameter of the objective function z2. Fig. 7 shows the potential distribution obtained on the metal surface after multiobjective optimization, where the boundary node 1 has coordinates (0, 0). This figure presents the results considering the original BVP and the corresponding BVP in order to analyse the pre-established protection criterion and the equivalence between the problems. The optimum location of the anodes determined by GA after 91 generations is shown in Fig. 8 (corresponding BVP) and Fig. 9 (original BVP), where the similarity of the results in the electrolyte can also be seen.
6. Concluding remarks The results here presented showed that the proposed strategy is effective in the optimization of the anode position in a galvanic cathodic protection system using GA and MFS. Furthermore, the constant potential circular anodes were successfully represented by a Gaussian function whose parameters were reached using a multiobjective optimization methodology. The dual reciprocity method (DRM) was implemented using a
Fig. 6. Potential in the electrolyte after GA optimization (corresponding BVP).
Fig. 7. Potential distribution on the metal.
W.J. Santos et al. / Engineering Analysis with Boundary Elements 73 (2016) 35–41
41
References
Fig. 8. Potential in the electrolyte after GA optimization (corresponding BVP).
Fig. 9. Potential in the electrolyte after GA optimization (original BVP).
positive definite compactly supported radial basis function (CSRBF) of Wendland. The Wendland CS-RBF guaranteed a good approximation for the source term (Gauss function) and the solvability of the linear system used to solve a particular solution of the Poisson equation. The examples presented show satisfactory equivalence between the original problem (sacrificial anode CP system) and the corresponding BVP, where the sacrificial anodes were introduced as source terms using the probability density function (Gaussian). Thus, it was not necessary to redefine the boundary geometric nodes and the arrangement of virtual sources for the standard method of fundamental solutions (MFS) formulation in the process of optimal positioning of sacrificial anodes. The adopted MINPACK routine LMDIF showed satisfaction to solve the nonlinear least squares problem that is a result of the nonlinear boundary conditions in the context of the MFS formulation. However, a large number of function evaluations (NFEV) were required for the convergence of the minimization process. This is due to the large number of design variables of the nonlinear least squares problem. In general, increasing the number of design variables slowed down the convergence of the minimization process.
Acknowledgements The research was sponsored by CNPq-Brazil.
[1] Fontana MG, Greene ND. Corrosion engineering.New York: McGraw-Hill; 1967. [2] Roberge PR. Handbook of corrosion engineering.New York: McGraw-Hill; 1999. [3] Santos WJ, Santiago JAF, Telles JCF. An application of genetic algorithms and the method of fundamental solutions to simulate cathodic protection systems. Comput Model Eng Sci 2012;87:23–40. [4] Santos WJ, Santiago JAF, Telles JCF. Optimal positioning of anodes and virtual sources in the design of cathodic protection systems using the method of fundamental solutions. Eng Anal Bound Elem 2014;46:67–74. [5] Montoya R, Rendon O, Genesca J. Mathematical simulation of cathodic protection system by finite element method. Mater Corros 2005;56:404–11. [6] Montoya R, Aperador W, Bastidas DM. Influence of conductivity on cathodic protection of reinforced alkali-activated slag mortar using the finite element method. Corros Sci 2009;51:2857–62. [7] Montoya R, Galvan JC, Genesca J. Using the right side of Poisson's equation to save on numerical calculations in FEM simulation of electrochemical systems. Corros Sci 2011;53:1806–12. [8] Miyasaka M, Hashimoto K, Aoki S. A boundary element analysis on galvanic corrosion problems - computational accuracy on galvanic fields with screen plates. Corros Sci 1990;30:299–311. [9] Abootalebi O, Kermanpur A, Shishesaz MR, Golozar MA. Optimizing the electrode position in sacrificial anode cathodic protection systems using boundary element method. Corros Sci 2010;52:678–87. [10] Telles JCF, Mansur WJ, Wrobel LC, Marinho MG. Numerical simulation of a cathodically protected semisubmersible platform using PROCAT system. Corrosion 1990;46:513–8. [11] Santiago JAF, Telles JCF. On boundary elements for simulation of cathodic protection systems with dynamic polarization curves. Int J Numer Methods Eng 1997;40:2611–22. [12] BEASY User Guide Beasy. Computational mechanics. Ashurst, Southampton, UK: BEASY Ltd.; 2000. [13] Azevedo JPS, Wrobel LC. Nonlinear heat conduction in composite bodies: a boundary element formulation. Int J Numer Methods Eng 1988;26:19–38. [14] Aoki S, Amaya K. Optimization of cathodic protection system by BEM. Eng Anal Bound Elem 1997;19:147–56. [15] Wrobel LC, Miltiadou P. Genetic algorithms for inverse cathodic protection problems. Eng Anal Bound Elem 2004;28:267–77. [16] Kupradze VD, Aleksidze MA. Aproximate method of solving certain boundaryvalue problems. Soobshch Akad Nauk Gruz SSR 1963;30:529–36. [17] Costa EGA, Godinho L, Santiago JAF, Pereira A, Dors C. Efficient numerical models for the prediction of acoustic wave propagation in the vicinity of a wedge coastal region. Eng Anal Bound Elem 2011;35:855–67. [18] Fontes Jr. EF, Santiago JAF, Telles JCF. On a regularized method of fundamental solutions coupled with the numerical Green's function procedure to solve embedded crack problems. Eng Anal Bound Elem 2013;37:1–7. [19] Karageorghis A, Fairweather G. The method of fundamental solutions for the solution of nonlinear plane potential problems. IMA J Numer Anal 1989;9:231–42. [20] Garbow BS, Hillstrom KE, Mor JJ. MINPACK project. Argonne National Laboratory, Argonne, Illinois, USA; 1980. [21] More L. Levenberg–Marquardt algorithm: implementation and theory. In: Proceedings of the biennial conference held at Dundee, June 28–July 1. Argonne, IL, USA: Argonne National Laboratory; 1977. [22] Partridge P, Brebbia CA, Wrobel LC. The dual reciprocity boundary element method.Southampton: Computational Mechanics Publications; 1992. [23] Yan JF, Pakalapati SNR, Nguyen TV, White RE. Mathematical modelling of cathodic protection using the boundary element method with nonlinear polarisation curves. J Electrochem Soc 1992;139:1932–6. [24] Atkinson KE. The numerical evaluation of particular solution for Poisson's equation. IMA J Numer Anal 1985;5:319–38. [25] Golberg MA. The numerical evaluation of particular solutions in the BEM - a review. Bound Elem Commun 1995;6:99–106. [26] Wendland H. Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree. Adv Comput Math 1995;4:389–96. [27] Golberg MA. The method of fundamental solutions for Poisson's equation. Eng Anal Bound Elem 1995;16:205–13. [28] Chen CS, Golberg MA, Schaback RA. Recent developments of the dual reciprocity method using compactly supported radial basis functions. In: Transformation of domain effects to the boundary. Southampton, Boston: WIT Press; 2003. p. 183–225. [29] Fletcher R. Practical methods of optimisation.Chichester, UK: John Wiley & Sons; 1987. [30] Srinivas N, Deb K. Multiobjective optimization using nondominated sorting in genetic algorithms. Evol Comput 1994;2:221–48. [31] Michalewicz Z. Genetic algorithms þ data structures¼ evolution programs. Springer-Verlag, New York: Springer-Verlag, Berlin Heidelberg; 1996. [32] Mitchell M. An introduction to genetic algorithms.Cambridge, MA: The MIT Press; 1997.