MODELING CATHODIC PROTECTION FOR PIPELINE NETWORKS

By

DOUGLAS P. RIEMER

A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL

OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT

OF THE REQUIREMENTS FOR THE DEGREE OF

DOCTOR OF PHILOSOPHY

UNIVERSITY OF FLORIDA

2000

ACKNOWLEDGEMENTS

This work was jointly supported by the Pipeline Research Council Interna-

tional and the Gas Research Institute through Contract Number PR-101-9512.

TABLE OF CONTENTS

Page

ACKNOWLEDGEMENTS ................... ........ ii

LIST OF TABLES ................................viii

LIST OF FIGURES ................... ............ x

KEY TO SYM BOLS ...............................xix

A BSTRACT . . . . . . . . . . . . . . . . . .xxi

CHAPTERS

1 PIPELINES AND CATHODIC PROTECTION ............... 1

1.1 Introduction .. .............................. 1

1.2 Electrode Kinetics (Butler-Volmer) ................ . 7

1.3 Mixed Potential Kinetics ........................ 10

1.4 Application to Corrosion in an Electrolyte ....... .... 11

1.5 Passivation and Film Formation ........... ...... ... 12

1.6 Coatings ........... ............... ....... 13

1.7 Anode Polarization ................... ........ 16

1.7.1 Galvanic Anodes ................... ..... 16

1.7.2 Impressed Current Anodes ........ ......... 17

1.8 Conclusions ................... ......... 19

2 MATHEMATICAL DEVELOPMENT ..................... 21

2.1 Soil Domain ..... ........................... 21

2.2 Justification of Dilute Approximation ........ ....... 24

2.3 Internal Domain ...... ............ ........ 25

2.4 Quasipotential Transformation ........ .......... . 26

2.5 Conclusions .......... ......... .......... 28

3 SOLUTION OF GOVERNING PDE ................. .

3.1 Boundary Element Method

3.2 Symmetric Galerkin Bound

3.3 Boundary Conditions ..

ary Element Met

3.3.1 Bare Steel ..............

3.3.2 Coated Steel ..............

3.3.3 Anodes .................

Infinite Domains ..............

Half Spaces ...................

Layers . . . . . . . . . . .

Development of Finite Element Method for

3.7.1 Pipe Shell Elements ..........

3.7.2 Applying Elements to the FEM ...

Bonds and Resistors ..............

Coupling BEM to FEM ............

Current Flow in the Pipe ...........

Conclusions ................

:hod

Internal

Domain

4 APPLYING KINETICS TO THE BEM AND FEM .. .........

4.1 Pipe Discretization ............................

4.1.1 Discretization of the Boundary Element Method .......

4.1.2 Discretization of the Finite Element Domain .........

4.2 Self-Equilibration .............................

4.3 M multiple CP Systems ...........................

4.4 Coating H holidays .............................

4.5 Nonlinear Boundary Conditions without Attenuation in the Pipe

Steel

. . . .. ... . . . ...... . . .. . .. 62

4.6 Nonlinear Boundary Conditions with Attenuation in the Pipe Steel

4.7 Solution by Successive Substitutions .. ...............

4.8 Solution of Combined System by Newton-Raphson Technique . .

4.9 Variable Transformation to Stabilize Convergence .. ........

4.10 Conclusions .. . .. ... ... ... ... .. ... ... ... ..

5 MESH GENERATION .. . ......................

5.1 Basic M esh . . . . . . . . . . . . . . . . .

5.1.1 Discretization and Shape Functions .. ...........

5.1.2 Continuity of the Shape Functions .. .............

5.1.3 Quadratic Iso-Parametric Rectangles .. ............

5.1.4 Quadratic Triangles .. . ....................

5.2 Long Pipes . . . . . . . . . . . . . . . . .

5.3 H holidays . . . . . . . . . . . . . . . . .

5.3.1 Round Holidays .. . .....................

3.4

3.5

3.6

3.7

3.8

3.9

3.10

3.11

np~rpl nnm Pnf

5.3.2 Rectangular Holidays

Pipe Crossings and Shadowing

Soil Type Divisions .......

Conclusions .. .........

6 SOLUTION ACCURACY .......................... 96

6.1 Sources of Error ............... . . . . ...... 96

6.2 Integration Accuracy .................. ....... 97

6.3 Integration Techniques .................. ....... 98

6.3.1 Adaptive Integration ................... . . 99

6.3.2 Singular Integrals .......... . . . . . . ......101

6.4 Matrix Inversions .......... ................... 101

7 MODEL VERIFICATION .......................... 103

7.1 Primary Current Distribution to a Disk ................ 104

7.1.1 Disk in a Half Space ......... .............. 104

7.1.2 Disk-Shaped Coating Holiday ................ 106

7.2 Comparison to the Solution of Kasper ................ 108

7.3 Dwight's Equation .. . .... .................. 108

8 SOLUTION VISUALIZATION AND INTERFACE DESIGN ....... 113

8.1 Interface Design .............................. 113

8.2 Data Visualization ........... .................. 115

8.3 2D Solution Visualization ......... ............... 116

8.4 3D Visualization ......... ................... 119

8.4.1 Using OpenGL ..........................119

8.4.2 Coding OpenGL for BEM and FEM Meshes .......... 123

8.5 Visual Feedback .............................. 124

9 COMPUTATIONAL PERFORMANCE . . . . . . . ... 125

Parrallel Computing Methods .....

Processor Arrays .. ..........

Distributed Computation Methods .

Communication Between Computers.

Improving Performance by Separating

Benchmarks and Scaling ........

Hard Drive Performance ........

High Performance Libraries ......

Calculation Code

10 HOLIDAY EFFECTS ........................... 137

10.1 Single Pipes .......... ..... ...... ......... 138

10.2 Two Pipes with a Coating Flaw ... . . . . .... 140

10.3 Pipes Exhibiting Stray Current ... . . . . .... 142

10.4 Soil Surface On-Potentials .......... ....... ........ 146

10.5 Adjacent Holidays .................. ....... . 147

11 COUPONS ................... ............ ..150

11.1 Introduction ....... ........ ....... ....... 150

11.2 Boundary Condition Parameters ..... . . . . . ... 153

11.3 Approach ......................... ...... 154

11.4 Results and Discussion .................. .. .... .. 156

11.5 Coupon Performance Diagram ... . . . . . . 161

11.6 Conclusion ................ . . . . . . ... 163

12 TANK BOTTOMS ............................. 165

12.1 Tank Bottom Mesh ................... . . . .. 166

12.2 Model for oxygen consumption ... . . . . .... 168

12.3 Implementation ............. . . . . . . .. 172

12.4 Results ....... .............. 173

12.4.1 Remote Ground-bed ..... . . . . . . ....173

12.4.2 Distributed Anodes. .... . . . . . . .... 178

12.4.3 Ribbon Anodes ......... ................. 184

12.5 Conclusions .......... ......... .......... 191

13 MULTIPLE SOIL DOMAINS .................. .... 193

13.1 Continuity Equations ........ ........ ......... 193

13.2 Implementation ................... .......... 195

13.3 Steel Resistance in Multiple Soil Domains ........... . 197

14 CONCLUSIONS AND FUTURE WORK . . . . . . ... 200

14.1 Conclusions ................... ............ 200

14.2 Future Work ................................ 201

APPENDICES

A PARTIAL DERIVATIVES OF THE BOUNDARY CONDITIONS .. .. 203

B OpenGL CODE ......................

B.1 The OpenGL Device ................... ....... 205

B.2 The Drawing Cycle ................... ........ 208

B.3 Platform Independent Code ..... . . . . . . ....209

B.4 Using the Mouse for Interaction ... . . . . . ..212

B.4.1 Orienting an Object ................. .. . .. 212

B.4.2 Selecting Objects ............... . . . .. 218

C CODESTRUCTURE ............................222

C.1 Interface ....... ...... ... .... ... ......... 222

C.1.1 Data Management ................ ...... . 224

C.1.2 Data Interaction ......... . . . . . . ......227

C.2 Calculation Engine ....... . . . . . . . . ... 232

C.2.1 Assembly of Boundary Element Matrices . . . . ... 232

C.2.2 Assembly of Finite Element Matrices . . . . . ... 233

C.2.3 Solution of Nonlinear System ..... . . . . . ..... 234

C.2.4 Soil Surface On-Potentials ..... . . . . . . 239

C.3 Remote Integration .............. . . ... 239

C.3.1 Interface Definition ................. .. . . 240

C.3.2 Data Marshaling .................. .... . 242

C.3.3 Server Side. ............ . . . . . ... 243

C.3.4 Client Side ............ . . . . . ... 245

REFERENCES ................... ............. 252

BIOGRAPHICAL SKETCH... ............... ....... 263

. . . . 205

LIST OF TABLES

Table page

1.1 Parameters for common galvanic anodes. . . . . . .. 17

1.2 Parameters for the oxygen and chlorine evolution reactions. . . 19

5.1 Error of Constant, Linear, and Quadratic elements in representing

a circle whose circumference is 7r and has unit diameter. ...... .78

9.1 Comparison of location of integration source code. . . . ... 129

9.2 Normalized benchmark result data for many types of computers. 131

9.3 An example run of the scaling of a heterogeneous distributed com-

puter .......... ........ ........ ........ 133

9.4 An second example run of the scaling of a heterogeneous distributed

computer. This run shows better scaling because the load from

other users was less. .................. ...... .. 134

9.5 Conditions used for all hard drive performance tests . . . ... 135

9.6 Implementing Vector operations in the Residual Vector Calculation 135

9.7 Implementing Vector operations in building the A matrix . . 135

11.1 Polarization parameters used in the model. . . . . . ... 153

11.2 Coating parameters used in model. Note: These parameters were

determined after water uptake of the coating. . . . . . . 154

11.3 Results of two coupon model. The small coupon used a polariza-

tion curve with lim,02 1.0 /A/cm2 (curve 1). The large coupon

used lim,02 that was an order of magnitude smaller (curve 2). . . 160

11.4 Results of two coupon model. The large coupon used a polariza-

tion curve with ilim,2 =1.0 p/A/cm2 (curve 1). The small coupon

used a ilim,2 that was an order of magnitude lower(curve 2). . . 160

LIST OF FIGURES

Figure page

1-1 Post explosion flare-off of leaking gas from a burst pipeline. The

burst was later attributed to corrosion weakened walls of the pipe. 2

1-2 Failed section of the butane pipe. a) an overview showing the size

of the ruptured seam; b) a closeup of the failed section of pipe. The

burst seam in the pipe followed a line of severe corrosion. . . 3

1-3 Pipe and anode system used for cathodic protection of buried pipelines.

The anode is the smaller vertical cylinder on the right. The red

cylinder in-line with the connecting wire is a resistor that is usu-

ally inserted in the circuit to control the amount of current driven

to the pipe in order to limit hydrogen evolution. . . . . . 5

1-4 Pouraix diagram for iron. ................... ...... 13

1-5 Diagram of coating model. ................ ...... . 14

1-6 Typical polarization curves associated with bare and coated steel. 16

3-1 Symmetric integration. Lines between the two circles represent the

Green's function relation between the two integration points on

the surfaces dFi and d . . . ........ . . ... 34

3-2 Reflections of Green's function to account for boundaries with zero

norm al current. ............. . . . . . 39

3-3 Rectangular enclosure created with a primary Green's function and

6 reflections. .............. . . . . . .... 40

3-4 Total cathodic current driven to a pipe with an underlying non-

conducting region. .................. .......... 40

3-5 Fit of the data from Figure 3-4 to i = 64 + 0.2118. . . . . 42

+ .11 .. .. 4

3-6 Diagram of the shell element used to calculate the potential drop

within the pipe steel. The element is assumed to have no variation

in the ( direction. The curvilinear coordinate system is displayed

on top of the element. ................ . . . 43

3-7 Black lines connecting pipes and anode are bonds. The red cylin-

ders represent resistors. They are modeled as 1-D finite elements

that connect to the bond connection points. . . . . . ... 47

3-8 Normalized current vectors drawn on top of the pipe. Three com-

ponents of the vectors are calculated in cartesian coordinates even

though the flow of current is only in two directions. The fact that

the vectors are always tangent to the pipe indicates that the coor-

dinate transformation was done correctly. . . . . . ..... 53

4-1 Newton-Raphson line search flow diagram. . . . . . ... 65

5-1 Type of mesh needed for model. Three pipes are shown passing

through a point where the soil conductivity changes. . . ... 75

5-2 Circle discretized with 8 linear triangular elements. The area is

2 2r2 for the 8 triangles which is a 10.% error. . . . . ... 77

5-3 Four quadratic elements with eight degrees of freedom represent-

ing a circle .... . ................ ...... 79

5-4 Continuity of the basis functions. Two elements are shown. Basis

functions that share a node with a value of one must be continuous. 80

5-5 Use of triangles to increase the degrees of freedom in the center

region while maintaining CO continuity. . . . . . . . 81

5-6 Quadratic iso-parametric rectangular parent element. Numbering

of nodes and shape functions given in figure. . . . . ... 82

5-7 Quadratic iso-parametric triangular parent element. Numbering

of nodes and shape functions given in figure. . . . . ... 83

5-8 Decrease of element aspect ratio at a boundary condition change. 85

5-9 Current density distribution on a circular holiday. Parameter J is a

function of the average current density iavg, the inverse Tafel slope

, the holiday radius ro, and the soil resistivity K. . . . . . 87

5-10 A mesh that includes two round holidays. This is the output from

a standard mesh generator. Further processing is needed to map

the mesh to a pipe and adjust the holidays so they are round. . . 88

5-11 The same mesh as shown in Figure 5-10. Here the elements for the

holiday have been removed. The elements have been divided into

four triangles to illustrate the locations of the center nodes on each

side of the triangles. .................. ...... .. 89

5-12 The same mesh as shown in Figures 5-10 and 5-11. The nodes on a

holiday boundary that did not lie on the circle defining the holiday

have been moved for the large holiday. The task has not yet been

done for the small holiday. ................... . 90

5-13 A gouge at the bottom of a pipe that was meshed using rectangular

elements. .................. .............. .. 90

5-14 Closeup of one edge of a rectangular holiday where two triangles

have been used to increase the element density of the holiday. . 91

5-15 Mesh at a crossing of two pipes. Both pipes have the highest den-

sity of elements at the crossing point. . . . . . . . 93

5-16 Output of the Delauny mesh generator when creating a division

between two soil domains. The holes for the pipes need post-

processing before the mesh can be used. . . . . . . .. 94

5-17 Completed mesh with pipes passing through soil type division. . 94

6-1 Sub-element with a high aspect ratio subdivided by the adaptive

integration routine to reduce the aspect ratio of the subsections. . 100

6-2 Sub-element with a desired aspect ratio subdivided by the adap-

tive integration routine. ................ ...... .. 100

7-1 Primary current distribution on a disk electrode in a half-plane.

Model solution at center of disk is within 2.0% of the analytic so-

lution of i/iavg = 0.5. This view is looking along the z axis toward

the disk. .......... ........ ........ ........ 105

7-2 Comparison of primary current distribution on a disk electrode in

a half-plane. Red line without circles is the analytic solution, blue

line is the numerical solution from the disk shown in Figure 7-1. 105

7-3 Disk holiday on a pipe. Shading corresponds to the off-potential.

Note the draw-down of potential in the coated region around the

holiday. ......... ........................ 106

7-4 Comparison of current density at the center of a disk electrode in

a half-plane (solid line) to that of a disk shaped holiday on a pipe

(circles). ...... .......... ...... .. ........107

7-5 Solution to the current distribution on a cylinder given by Kasper.

Plot represents the ratio of the maximum to minimum current den-

sity. At infinite separation, the ratio becomes 1. . . . . ... 109

7-6 Distribution of current around half the circumference when the

cylinders are placed at 3.0 diameters separation. The ratio of the

max current to the min is 0.72. Kasper's solution is 0.78 ....... .110

7-7 Error in Dwight's equation as a function of the separation of two

cylinders. Solution given by the 3 dimensional model approaches

Dwight's solution at large separation. The difference between Dwight's

solution and the model is due to end effects. . . . . . .... 111

8-1 A screen shot and movie (click on image above) that shows some of

the interactive graphics techniques needed for three dimensional

geom etry. ............ ... ...... .. ..... 113

8-2 A screen shot of the basic window used for the interface to the

model. The window contains elements familiar to anyone who

uses Microsoft@ Windows. ............... . ... .. 114

8-3 A screen shot of the basic window with an empty model. The win-

dow has two parts: the first is for control and feedback, the other

is for data display in 3 dimensions. . . . . . . ..... 115

8-4 Red line represents the points of the pipe surface that are extracted

to create a X-Y plot. X values are the how far from the starting

end of the pipe the value of the potential, current, off-potential, or

corrosion current is taken. ................... ..... .117

8-5 The red line represents the points of the pipe surface that are ex-

tracted to create a X-Y plot. X values are the how far from the start-

ing end of the pipe the value of the potential, current, off-potential,

or corrosion current is taken. ................. . .118

8-6 Shaded mesh of boundary element calculation grid. OpenGL al-

lows instant visual inspection of the mesh for correctness. ..... .121

8-7 Same view as figure 8-6 on page 121, with four times the number of

elements. This has been done to better approximate the quadratic

nature of the real elements using linear interpolation for the graph-

ics elements ........ ........ ............. 122

9-1 Plot of performance vs expected performance for a distributed com-

puter ................... ....... ........ 129

10-1 Pipe that exhibits many coating holidays. Pipe was coated in the

late 60's and dug up in 1996. ................... . 137

10-2 Off-potential distribution for a large holiday in a good quality coat-

ing. The edge of the coating next to the holiday does not meet the

-850 mV criterion while the holiday does. Even so, the corrosion

rate at the holiday is much higher than at any of the coated sections. 139

10-3 Current density and potential distribution around the circumfer-

ence of a pipe through a holiday. Zero degrees is the top of the

pipe. The blue line is the current density. . . . . . . 140

10-4 Current density and potential distribution around the circumfer-

ence of a pipe through an under-protected holiday. Zero degrees

is the top of the pipe. The blue line is the current density. . . 141

10-5 Current density and potential distribution around the circumfer-

ence of a pipe at the spot closest to the under-protected holiday on

the neighboring pipe. Zero degrees is the top of the pipe. The blue

line is the current density. .................. ... .142

10-6 Distant view of a under-protected holiday. Lower right portion of

pipe with holiday is also exhibiting stray current. . . . ... 143

10-7 Closeup of an under-protected holiday. . . . . . . . 144

10-8 Current Density around pipe in area exhibiting stray current in

Figure 10-7. Positive current densities are anodic. . . . ... 145

10-9 Potential distribution around pipe at center of section exhibiting

stray current and at end of section where all parts of the circumfer-

ence of the pipe are cathodic. ................. . .145

10-10Distant view of soil surface on-potentials over an under-protected

holiday. ................... ........ ........146

10-11Distant view of soil surface on-potentials with a portion of the soil

surface removed to see the pipe. The soil surface on-potential only

indicates the presence of the anodes. The holiday is not seen. . . 147

10-12Setup for calculation of a pinhole holiday next to a large holiday. 148

10-130ff-potential distribution of a pinhole holiday next to a large hol-

iday. The large holiday is under-protected by standard criteria

while the pinhole is well protected. . . . . . . ...... 149

11-1 Schematic representations of the two cylindrical coupons used in

the simulations: a) large coupon with a surface area of 50 cm2;

and b) smaller coupon with a surface area of 9.5 cm2 . . . ... 155

11-2 Large coupon Difference between the most positive potential of

the holiday and the coupon at three different rectifier outputs that

result in the coupon potential being measured at -750, -850 and

-950 mV (CSE). The negative sign on y-axis label indicates the hol-

iday is at a potential cathodic to the coupon. The lack of a sign

indicates the opposite. Panels a) through d) correspond to differ-

ent values of J = Aholiday/Acoupon. ... . . . . . . .. 157

11-3 Small coupon Difference between the most positive potential of

the holiday and the coupon at three different rectifier outputs that

result in the coupon potential being measured at -750, -850 and

-950 mV (CSE). The negative sign on y-axis label indicates the hol-

iday is at a potential cathodic to the coupon. The lack of a sign

indicates the opposite. Panels a) through d) correspond to differ-

ent values of J = Aholiday/Acoupon. ... . . . . . . .. 158

11-4 Coupon performance diagram for coupons placed in the same soil

environment as the pipeline coating defects that expose bare steel.

The soil resistivity was 10,000 Q cm. Black symbols correspond

to calculations where a protected coupon corresponded to a pro-

tected coating holiday. Open symbols correspond to calculations

where the coupon could indicate adequate protection for pipes

with coating holidays that are seriously underprotected. The gray

symbols correspond to cases where the coupon reading is not con-

servative, but the error was less than 15 mV. . . . . . . 162

12-1 High resolution mesh used for the boundary element mesh. . . 167

12-2 Oxygen concentration profile within the soil at several values of

the scaling parameter R k/D. ..... . . . ... 171

12-3 Potential distribution when the oxygen concentration is uniform

at its largest value. This condition occurs when the tank has been

recently filled after being empty. The grid spacing is 10 ft square. 175

12-4 Radial potential and current density distribution when the oxygen

concentration is uniform and at c,. Potential is given by the red

line. The black line indicates where a uniform current density dis-

tribution would be. ............. . . . . . . .......176

12-5 Potential distribution when the oxygen concentration follows the

profile given in Figure 12-2 corresponding to R k/D =2.0. The

grid spacing is 10 ft square. .............. . . . 177

12-6 Radial potential and current density distribution when the oxygen

concentration follows the profile given in Figure 12-2 correspond-

ing to R k/D = 2.0. Potential is given by the red line. The black

line indicates where a uniform current density distribution would

be.. . . . ....... . . . . . . . . . . . .. 178

12-7 Potential distribution with distributed anodes and uniform oxygen

concentration profile. ................ . . . . 180

12-8 Current density distribution on tank bottom with distributed an-

odes and uniform oxygen concentration. . . . . . . . 180

12-9 Potential distribution with distributed anodes and nonuniform oxy-

gen concentration profile. ............... . . ... . 182

12-10Potential and current density distribution on a tank bottom with

distributed anodes and nonuniform oxygen concentration. . . 182

12-11Potential distribution with distributed anodes, underlying insulat-

ing layer at 2 ft and uniform oxygen concentration profile. . . 183

12-12Potential distribution with distributed anodes, underlying insulat-

ing layer at 2 ft and nonuniform oxygen concentration profile. . 184

12-13Tank bottoms with 5 and 11 parallel impressed current anodes. . 185

12-14Potential distribution for tank bottoms with 5 and 11 parallel im-

pressed current anodes with secondary containment. . . ... 186

12-15Axial current density distributions on anodes for both the 5 and 11

anode cases. ........ ........ .............. 188

12-16Potential and corrosion current distribution for tank bottoms with

5 parallel impressed current anodes with no secondary containment. 188

12-17Potential distribution for tank bottoms with 5 and 11 parallel im-

pressed current anodes and secondary containment with nonuni-

form oxygen distribution. ................ . ... . 189

12-18Potential distribution for tank bottoms with 5 and 11 parallel im-

pressed current anodes with nonuniform oxygen distribution. . 190

13-1 Plane meshed for a soil type change. Three pipes are shown pass-

ing through a point where the soil conductivity changes. . . 194

13-2 Three domain model setup utilizing two soil division meshes. The

soil surface grid is drawn using 1 mile squares. There are two pipes

of 24 in diameter in a 50 ft right-of-way. . . . . . .... 198

xvii

B-1 A movie demonstrating how the user can use the middle mouse

button the change the OpenGL view of two pipelines. This is one

of the modes of movement. The movie is accessed by clicking on

the above image. .................. ........ 214

C-1 Interface with two section:One for feedback and control, one for

3D data display. ............. . . . . . .... 223

C-2 Console Window ......... . . . . . . 228

xviii

KEY TO SYMBOLS

D = Potential and potential distribution

V = Voltage distribution in the pipe steel and anode material

i = Current density

Ei = Equilibrium potential for species i

ji = Tafel slope

F = Faraday's constant

R = Gas constant

T = Temperature

k = Kinetic parameter

K = Conductivity

ci = Concentration of species i

t = Time

Ji = Flux vector of species due to fields

Ni = Net flux vector of species i

Ri = Rate of generation of species i

Di = Diffusion coefficient for species i

zi = Charge on species i

ui = Mobility of species i

I = Ionic strength

qi = Kirchhoff transformation variable for species i

Q = Quasi-potential

Q = Domain for which a differential equation applies

F = Boundary of a domain

n = Outward normal vector of boundary

w = Arbitrary weighting function

G = Green's function

J = Source point in Green's function

x = Field point in Green's function

6 = Dirac delta function

r = Euclidian distance between two points

G = Matrix of coefficients from the right hand integral of 3-9

H = Matrix of coefficients from the left hand integral of 3-9

A = Matrix of columns from the H and G matrices

Y = defined as V <

n = Outward normal vector of boundary

K = Matrix resulting from Finite Element Method

n = total number of nodes in a problem

ne = total number of elements in a problem

nen = number of nodes in an element

nce = number of elements having a given node in common

XX

Abstract of Thesis Presented to the Graduate School

of the University of Florida in Partial Fulfillment of the

Requirements for the Degree of Doctor of Philosophy

MODELING CATHODIC PROTECTION FOR PIPELINE NETWORKS

By

Douglas P. Riemer

December 2000

Chair: Dr. Mark E. Orazem

Major Department: Chemical Engineering

An extensive network of buried steel pipelines is used throughout the world

to transport natural gas and liquid petroleum products. In the United States

alone, there are over 1.3 million miles of pipe used to transport natural gas. Fail-

ure to prevent external corrosion of pipelines can have disastrous consequences.

A butane line explosion, which cost the lives of two persons, was caused by exter-

nal corrosion attributed to inadequate cathodic protection. Proper cathodic pro-

tection is particularly important today because the pipelines constructed during

the construction boom of the 1960s and 1970s are aging, and the encroachment of

suburban housing construction increases the potential for property damage and

loss of life.

Standard cathodic protection design equations, based on the assumptions

that pipes are isolated and have a uniform bare metal surface, are not valid for

modern complex pipeline networks. The object of this work was to develop a

mathematical model that can describe correctly the complex interactions seen in

cathodic protection of complex pipeline networks. The governing differential

equation is solved in 3-dimensions for the conductance of current through the

soil, steel, and connections between pipes and anodes. The model accounts for

the kinetics of the corrosion, mass-transfer-limited oxygen reduction, and hydro-

gen evolution reactions. It also accounts for the nonuniform current and potential

distributions due to pipe-pipe interactions, coating flaws, and anode locations.

Current and potential distributions are calculated around the circumference and

along the length of the pipes.

This dissertation describes the solution techniques needed to solve the sys-

tem of equations that couple nonlinear reaction kinetics, passage of current through

soil, and passage of current through pipeline steel. Sample calculations will be

presented to show the complicated nature of the interactions between adjacent

pipelines.

xxii

CHAPTER 1

PIPELINES AND CATHODIC PROTECTION

1.1 Introduction

Since the beginning of the 20th century, petroleum products and natural gas

have been transported over long distances by buried steel pipelines. There are

over 1.28 million miles of buried steel main-line pipe for the transport of natural

gas alone. 1 There are, in addition, 170,000 miles of pipeline for transport of crude

oil and refined products.2

These pipes are generally in soil environments which contain oxygen in an

electrolyte and therefore are subject to corrosion. By the late 1920s, the number of

leaks had begun to increase alarmingly, and, by the early 1930s, all major pipeline

owners were providing some measure of corrosion mitigation to their pipelines,

including application of coatings and cathodic protection.3 The application of cor-

rosion mitigation strategies has in general been successful. Most pipeline failures

are now attributed to mechanical damage. Nevertheless, over a six-year period,

(1994-1999), 408 separate incidents were attributed to corrosion, resulting in a net

loss of $75 million and 4 lives.* Additional failures, later deemed to have been

caused by corrosion, are listed under "other." It is therefore difficult to assess,

Data compiled from the U.S. Office of Pipeline Safety statistics for 1994-1999 for transmission

and distribution of liquid petroleum products and natural gas.

Figure 1-1: Post explosion flare-off of leaking gas from a burst pipeline. The burst

was later attributed to corrosion weakened walls of the pipe.2

from government statistics, the net loss of property and life that can be attributed

to failure to provide proper corrosion protection.

For example, a pipeline that ran through a residential neighborhood in a small

Texas valley burst releasing butane in 1996.2 Two persons tried to escape by

driving away, but the ignition system of their vehicle ignited the vapor cloud,

killing them. The post explosion results of the effects of unmitigated corrosion

are shown in Figure 1-1.

Investigation of the section of pipe showed that the burst was caused by weak-

ening of the steel through corrosion from the outside of the pipe Figures 1.2(a)

and 1.2(b). It was determined in the investigation that the cathodic protection

systems were not properly maintained to U.S. DOT standards. The pipeline op-

erating company was fined (January 2000) 30 million dollars for this and other

non-fatal leaks.4

There are other examples of catastrophic failure due to corrosion, but the

larger problem is often small leaks due to coating failure and insufficient cathodic

protection. There are currently 7 bills before congress (October 2000) that address

(b)

Figure 1-2: Failed section of the butane pipe. a) an overview showing the size of

the ruptured seam; b) a closeup of the failed section of pipe. The burst seam in

the pipe followed a line of severe corrosion.2

pipeline safety. The joint bill5 and amendments6 propose new research that in-

cludes a call for better understanding of corrosion and cathodic protection of

pipelines.

Typically such catastrophic failures are prevented by aggressive monitoring of

the cathodic protection systems coupled with inspection of suspect regions of the

pipeline. Yet, wall loss due to corrosion is sometimes evident even on pipes that,

according to standard CP monitoring techniques, are nominally well protected.

There is clearly a need to improve our understanding of cathodic protection. The

issue of cathodic protection is particularly relevant today because the pipelines

constructed during the construction boom of the 1960s and 1970s are aging, and

the encroachment of suburban housing construction increases the potential of

property damage and loss of life. The design of cathodic protection systems is

complicated by the current practice of placing an increasing number of pipes in a

right-of-way, and by the increasing number of crossings due to more pipes being

laid. The combination of these factors greatly increases the potential for stray

current effects.

Cathodic protection is a means in which an undesirable reaction is replaced by

a more desirable one on another surface. In the case of buried metallic structures,

the undesirable reaction is the dissolution of the metal.

M -- M"+ + ne (1-1)

Cathodic protection works by having a second surface within the conductive soil

or electrolyte at a more negative potential than the pipe steel would have. This

surface is then connected to the pipe with a wire as shown in Figure 1-3.

5

Figure 1-3: Pipe and anode system used for cathodic protection of buried

pipelines. The anode is the smaller vertical cylinder on the right. The red cylinder

in-line with the connecting wire is a resistor that is usually inserted in the circuit

to control the amount of current driven to the pipe in order to limit hydrogen

evolution.

Since the anode and pipe are electrically linked, the more active (negative)

surface will tend to favor a reduction reaction, while an oxidation reaction will

be favored on the more positive surface. With careful selection of anode metals,

the rate of oxidation of the pipe steel can be reduced to a level which will result

in a much longer service life.

High resistance coatings are used on most pipelines to reduce the amount

of current necessary to protect the pipeline. However, the introduction of nicks

and scrapes which expose bare steel cannot be avoided during installation. Its

is therefore absolutely necessary to use cathodic protection (CP) on coated pipes

since experience has shown that not using CP results in accelerated failure.7 Even

though CP is necessary for coated pipes, the cost savings when coatings are used

makes the practice desirable.

A better understanding of the behavior of a coated pipeline under cathodic

protection can be obtained by solving a mathematical model of the governing

differential equations. A model is usually set up in a rectangular domain where

the boundary conditions are constant on any given side. However, the boundary

conditions for pipelines under cathodic protection are not uniform or constant on

any plane (i.e., ( = OV at z = Ift or n Vc = 1/A/cm2 at r = R). There are three

reasons for this. The first is that the pipes are cylindrical. Laplace's equation for

right circular parallel cylinders can be solved analytically if the boundary condi-

tion is constant on the surface of the pipe.8 However, this method of solution can-

not be used for the problem treated here because the anodes do not run parallel

to the pipe (see Figure 1-3). Second, coating failures can be discrete. One cannot

use the uniform boundary condition (( = c at z = 0 or n Vc = c at r = R) since

the potential and current density will be different for a discrete holiday than for

coated steel. The boundary conditions are nonlinear mathematical expressions

that relate the net current density of the chemical reactions on the pipe surface to

the potential driving force.

1.2 Electrode Kinetics (Butler-Volmer)

Some assumptions about the knowledge of the reader are made in the fol-

lowing sections. It is assumed he has a basic knowledge of electrochemistry.

Good sources of the necessary concepts can be found in Bockris and Reddy9,10

and in Bard and Faulkner.11 The best treatment of the mathematics describing the

physics is in Newman.12

The first step taken in developing a mathematical model for cathodic protec-

tion of steel is to understand the role of the electrochemical reactions that take

place on the steel and anode surfaces.

Bockris et al. 13 have developed a pH-dependent mechanism for an iron elec-

trode. They proposed that the electrochemical reactions that take place on the

surface of an electrode may be developed in the following way:

Iron corrodes according to two mechanisms depending on the pH. Buried

structures under cathodic protection in soil form alkaline environments.3 The cor-

rosion step for high pH values is

R.D.S.

Fe + 20H- R- Fe(OH)2 + 2e (1-2)

and

Fe(OH)2 = Fe2 + 20H

(1-3)

where R. D. S. stands for Rate-Determining Step. The sum of reactions (1-2) and

(1-3) results in the net reaction

Fe ~ Fe2+ + 2e (1-4)

which is an anodic reaction when the reaction proceeds from left to right.

For iron deposition the reactions are written as

Fe(OH)2 = Fe2 + 20H (1-5)

and

R.D.S.

Fe(OH)2 + 2e R Fe + 20H (1-6)

which are cathodic reactions when proceeding from left to right.

Newman provides a development for reactions with multiple single electron

steps in Chapter 8 of his text.12 The anodic current density as a function of pH

and potential was given by Bockris et al. to be13

ia = ka,FeaoH e[ (1-7)

where ka,FE is the kinetic rate constant, aOH- is the activity of OH ions, V is de-

fined as Ometal Osolution, F is Faraday's constant (96487 coulombs /equiv), R is

the Gas Constant and T is temperature. The activity of the OH ions is linear

with respect to pH.

The cathodic current was given by

ic = kc,FeaFe2+aOH- e[r (1-8)

where aFe2+ is the activity of Fe2+ ions, and the net current density i is given by

i = ia + ic

(1-9)

At equilibrium, where only equations (1-7) and (1-8) are taking place, the equi-

librium potential Vo can be found by equating equations (1-7) and (1-8) and solv-

ing for V

WV =RT nKc,FeaFe2+,e (1-10)

F Ka,Fe

where Vo is defined as the equilibrium potential. Then the exchange current den-

sity can be defined to be

r VoF

io = -a,FeaOHe- ] (1-11)

using the anodic term. The cathodic term would have the same current density

at Vo. Substituting in the expression for Vo into equation (1-11) obtains the proper

expression

1/4 3/4 3/4 (1-12)

Lo = Ka,Fe,FeFe2+aOH (1-12)

Equation (1-9 can be rewritten in terms of the exchange current density as

i= io (e[ e[ F]) (1-13)

or

i = io2 sinh (1-14)

(RT

The number of unknowns in equation (1-13) can be reduced by changing the

base from e to 10 and defining Tafel slopes as

2.303F

a (1-15)

RT

2.303F

/3c = (1-16)

RT

and then moving all pre-exponential terms into the exponent to get

(1-17)

i= 10[Ei 10[E ]

where Ea is fit to experimental data. Using large-scale experiments, ARCO per-

sonnel found values of Ea to be around -0.56 V referenced to Cu/CuSO4.14

1.3 Mixed Potential Kinetics

For freely corroding systems, conservation of charge requires that the anodic

corrosion reaction, e.g.

Fe Fe2+ + 2e (1-18)

be balanced by a cathodic reaction such as

02 + 2H20 + 4e -, 40H (1-19)

which is called the oxygen-reduction reaction. The expression for the rate of

the oxygen reaction can also be assumed to be Tafelian in nature. Thus, a two-

parameter expression can be obtained in the same way as was done for the iron

reactions. When an active metal anode (Zn or Mg) is connected to a pipe, the po-

tential shifts sufficiently negative and a second cathodic reduction reaction takes

place

2H20 + 2e -, H2 + 20H (1-20)

Reaction (1-20) is commonly called the hydrogen evolution reaction. The mecha-

nism of the hydrogen evolution reaction is currently a topic of research.15' 16,17,1819

The hydrogen evolution reaction follows Tafel kinetics so an equation of the same

type as the cathodic reaction for iron is used. It also can have the number of un-

knowns reduced to two following equation (1-17)

i VEaH2] 1] Ec,H2

i = 10 aH2 10.L cH2 J (1-21)

The corrosion potential Vcorr is the potential difference between a reference

electrode and a freely corroding surface. The corrosion potential differs from

equilibrium potential described by equation (1-10) because more than one reac-

tion takes place. The corrosion potential is also a function of the oxygen con-

tent and transport characteristics for oxygen within the electrolyte of the sys-

tem. When the level of oxygen is low, a more negative corrosion potential is seen

which can be associated with a lower corrosion rate.

The total current is given by the sum of equations of the form of (1-17). There

is one for each reaction, and the sum can be written as

i = 10L 10 ,e +

SUV-Ea,02 V c,02

10L _- 10L +

10 ,H2 1 L H (1-22)

Equation (1-22) can often be simplified by observing that some terms can be ig-

nored in the potential range of interest. A corrosion potential may be calculated

from this equation by setting it equal to zero and solving for V = Vcorr. This is the

same potential that is measured under a zero current experimental condition. A

potential is taken with a reference electrode just outside the surface of the metal.

1.4 Application to Corrosion in an Electrolyte

If there are two different interconnected metals contained in a continuous elec-

trolyte, one metal will corrode at a higher rate than it would alone while the other

metal will corrode at a lower rate. If the kinetics are known, they can be used as

boundary conditions in a numerical method. Equations of the form of (1-22) are

used to describe the kinetics at the boundary. Often, one of the reactions will

have a mass-transfer limitation. This is because the reactant must come from the

electrolyte. One common form is one that was published by Nisoncioglu20 21,22

and adapted for steady state soil systems by Yan et al. 23 and Kennelly et al. 24 For

bare steel it takes the form

VEFe V 02 2

V--e V- (~E 1 -(V- EH2)

i = 10 re 10 "02 10 "H2 (1-23)

lim,02 /

where ihm,02 is the mass transfer limited current density for oxygen reduction.

This means that the portion of the current due to oxygen reduction cannot exceed

the value ilm,02. There is no such limitation for the hydrogen-evolution-reaction,

since it is assumed that there is always sufficient water next to the pipe for that

reaction to proceed under kinetic limitation.

1.5 Passivation and Film Formation

The production of OH by oxygen reduction and hydrogen evolution reac-

tions (1-19) and (1-20), respectively, increases the value of pH at the pipe surface.

The pH on bare steel has been reported to be on the order of 10-12.3 This is be-

cause the by product of corrosion is the formation of hydroxide ions. At this pH

level, Fe(OH)2 becomes insoluble and forms a film on the surface of the metal.

Such films can block the transport of oxygen and cause the corrosion potential to

change.25 One may use a Pourbaix diagram to determine what films are stable (if

any) at a given pH (see Figure 1-4).

06

|6) eS asA Fe(OH)3

0 Fe2+

(Fe2-02 )=1 ( 1 4) (9) (HFe0 1)= 10-'

la4)

-0.6 (5)Co

(Fe+) = 10-' (5) 10

Fe(OH), {2) (3)

HFeO;2-- 10

-1.4 Fe

1.8 ----

-2 0 2 4 6 8 10 12 14 16

pH

Figure 1-4: Pouraix diagram for iron.25

1.6 Coatings

To reduce the amount of current required to provide cathodic protection, coat-

ings are applied to the steel surface before burial. These coatings are usually

polymers but can sometimes be cement-based products.

When a coating is applied to a buried pipeline, cathodic protection must be

used, or the pipeline will experience perforation due to corrosion in a far shorter

time than would be experienced if the pipe were left bare. Isolated faults in the

coating (holidays) can form a galvanic couple with the coated parts, accelerating

the corrosion rate of the holiday.

It is then necessary to model the behaviour of coatings and the steel under the

coating when under cathodic protection. It is known that some current will pass

through some coatings after they have absorbed water.26,27,28,29,30 Corfias et al.

have shown that the pore structure expands as the coating absorbs water.27 They

also shown that the conductivity greatly increases with immersion time.

Coating as a solid transport Coating with pores that

inhibiter with diminished allow transport of ions and

diffusion coefficients 02

(Douter (Douter

inner __

Figure 1-5: Diagram of coating model.

Riemer and Orazem have modified equation (1-23) to take into account the

existence of an electrically resistive diffusion barrier.31'32 This barrier may be an

insoluble oxide film or a polymeric coating applied before the pipe was buried.

Figure 1-5 shows two possible behaviours of a conducting coating model. The

left side of the figure shows a uniform diffusion/electrical barrier. The right side

shows a barrier with small micropores through which an electrolyte can flow and

form an electrical connection.

The Potential drop through the film or coating can be expressed as 33

i= n (1-24)

p6

where ( is the potential in the electrolyte next to the coating, (in is the potential at

the underside of the coating just above the steel, p is the resistivity of the coating

and 6 is the thickness of the coating. The current density can also be written

in terms of the electrochemical reactions that take place under the coating if the

coating is viewed as either a porous structure or as a uniform diffusion barrier

E EO2 _V-in-EH2)

Apore EFe 1 E 1

-- = 10 -^e 10 02 10_H2

A (1 Oblock) ilim,02 2

(1-25)

A

where is the effective surface area available for reactions, and block is the

reduction to the transport of oxygen through the barrier. It is also assumed that

the coating has absorbed enough water that the hydrogen evolution reaction is

not mass-transfer limited. Equations (1-24) and (1-25) are solved simultaneously

to get a value of the current density. The current density is eliminated to give

A ((D- (in) V-in-EFe ( 1 V-j-E (Vin-EH2)

= 10 Fe 10 10 2

AporePfilmfilm (1 block) lim,02 10 H2

(1-26)

which is solved by a Newton-Ralphson method to get (in. The current density

can subsequently be readily calculated from equation (1-24).

Recent work by NOVA Gas34 and CC Technologies,29'35,36 showed that the

coatings on steel pipe form a diffusion barrier when placed in aqueous environ-

ments and that the coating absorbs water. They have shown that it is possible to

polarize slightly the steel under a disbonded coating.

Examples of polarization curves associated with equations (1-23) and (1-25)

are given in Figure 1-6. The corrosion potential can be calculated from equation

(1-25) since there is no current through the coating. If (1 Oblock) is smaller than

the pore ratio, then the metal will be more active than bare steel. If they are the

same, then the coated pipe and bare steel will have the same corrosion potential.

If (1 block) is larger than the pore ratio, the metal under the coating will be more

noble than bare steel and form a galvanic couple that enhances the corrosion rate

of the holiday with no CP.

100

c4- 10-1 -

0 10-2

S10-4

) High Oxygen Blocking

10-8 Low Oxygen Blocking

3 10 Bare Steel

10-9 -

10-1 I I I I

-1.50 -1.25 -1.00 -0.75 -0.50 -0.25

Potential (V CSE)

Figure 1-6: Typical polarization curves associated with bare and coated steel.

1.7 Anode Polarization

Anodes undergo similar electrochemical processes as the pipes. The corro-

sion term can again be expressed as one half of an oxidation-reduction reaction

where the half of interest is the corrosion (oxidation) of the anode material. Both

galvanic and impressed current anodes are used.

1.7.1 Galvanic Anodes

A galvanic anode is a metal that is more active than the metal that needs to be

protected. For carbon steel, zinc and magnesium are used in soil environments

and aluminum is used in sea water. Aluminum does not have a large enough

driving force to be effective in soil where the resistivity is much higher.

A simple model for a galvanic anode would be:

i = io2 (10(V- coErr)/ anode- 1 (1-27)

Table 1.1: Parameters for common galvanic anodes.7

Anode Type Eeql (mV CSE) 3 (mV/decade) io2 (pA/cm2)

Al -1.0 60 1.0

Zn -1.1 60 1.0

Mg (Standard) -1.5 60 1.0

Mg (High Potential) -1.75 60 1.0

where io2 is the mass-transfer-limited current density for oxygen reduction, Ecorr

is the free corrosion (equilibrium) potential of the anode and / is the Tafel slope

for the anode corrosion reaction. Hydrogen evolution was ignored for galvanic

anodes because it only makes a small contribution to the net current at operating

potentials.

The reduction reaction term (the negative contribution to equation (1-27)) was

assumed to be constant because the operating conditions are such that the re-

duction is always mass-transfer-limited. The reaction is mass-transfer-limited

because the driving force for the reduction reaction is higher than that for steel

which when freely corroding, is at the mass-transfer-limitation.

The three parameters of the model are generally known for all types of gal-

vanic anodes. Some typical values are given in Table 1.1.

1.7.2 Impressed Current Anodes

Impressed current anodes do not obtain their driving force from the activity

of one metal relative to another. Instead, an outside current source such as a

rectifier provides the activity. This means that the current densities obtainable on

impressed current systems are usually an order of magnitude higher than can be

obtained with galvanic anodes.

Impressed current systems consist of a dimensionally stable anode connected

to the positive terminal of a direct-current (DC) rectifier. The negative terminal

of the rectifier is connected to the pipe. The rectifier can then be used to push the

potential of the anode to any desired potential that is more negative than the pipe.

Impressed current system can protect a much longer section of pipe because they

can provide a larger potential driving force for current flow in the soil. They also

can be used to greater effect in highly resistive soils.

Since the anodes do not react, the usual reaction of metal dissolution, equa-

tion (1-4), is not appropriate. The most likely reactions are water oxidation and

chloride oxidation. They proceed respectively as follows

2H20 02 + 4H+ + 4e (1-28)

and

2C1 C12 + 2e (1-29)

The model for an impressed current anode is similar to that for galvanic an-

odes. An additional term is added to the exponent to account for the potential

setting of the rectifier, i.e.,

i = i2 10(Vo(_ rectifierEO2)/02 1) (1-30)

where A Vrectifier is the voltage added by the rectifier, V is the voltage of the anode,

( is the voltage just outside the surface of the anode, Eo2 is the equilibrium po-

tential for the oxygen evolution reaction, and /3o2 is the tafel slope for the oxygen

evolution reaction.

Table 1.2: Parameters for the oxygen and chlorine evolution reactions.25

reaction Equil Potential (E) Tafel Slope (3)

02 evolution -172 mV (CSE) 100 mV/decade

Cl2 evolution 50 mV (CSE) 100 mV/decade

If chloride ions are present, reaction (1-29) can occur which is the chlorine

evolution reaction, the reaction occurs in environments such as salt mashes and

estuaries where there is salt in the soil.

Parameters for reactions (1-28) and (1-29) are given in Table 1.2 except for io2

which must be determined experimentally.

1.8 Conclusions

This chapter has presented the problems facing engineers combating corro-

sion of buried pipelines. Cathodic protection is used to slow the rate of corrosion

but is not well understood for complicated pipeline networks. To further the

understanding of the operation of cathodic protection systems, a mathematical

model has been developed under the constraints of the kinetics are transport as-

sociated with corrosion in a soil system.

The model is developed in the next three chapters under the assumption that

the reader is familiar with the principles of the calculus of variations, vectors and

tensors, and transport in electrolytes. Chapter 11 of Newman's text provides a

good introduction to transport in electrolytes.12 Aris provides an excellent treatise

on vectors and tensors.37

Furthermore, the appendices assume the reader is well instructed in the C and

C++ programing languages and Microsoft Windows program development.

There are several good sources for C and especially C++ programing informa-

tion. One standard textbook is by Deitel and Deitel.38 One of the best books to

convey the concepts and power of C++ is by Bruce Eckel.39 There are a myriad

of books that concern themselves with how to write Microsoft Windows based

applications. Almost any bookstore will have an entire aisle devoted to them.

The mathematical model for cathodic protection of pipeline networks is re-

ported in the next chapters. Chapter 2 shows the derivations for the equations

that describe the flow of current through the soil and through the pipe steel.

Chapter 3 describes the numerical techniques that were used to solve the differ-

ential equations. Finally, Chapter 4 goes into the details of applying the nonlinear

mathematical expressions for the kinetics of pipes under cathodic protection as

boundary conditions for the differential equations. Chapter 4 also describes how

addition equations and unknowns are added to the numerical method to describe

several more aspects of the physics that are not accounted for in the differential

equations derived in Chapter 2, such as multiple independent cathodic protec-

tion systems within the same environment.

Results of the model and their implications for design of cathodic protection

systems are reported in Chapters 10 to 12. Specifically, Chapter 10 shows how

holidays affect neighboring pipes and other holidays in the vicinity.

The intervening chapters discuss how the model was verified, and how per-

formance and accuracy were obtained.

CHAPTER 2

MATHEMATICAL DEVELOPMENT

The model for cathodic protection must account for the flow of current in the

soil, in the pipes, and in the circuitry. Until recently, most models of cathodic

protection of pipelines assumed that the pipe steel was an equipotential surface.

While it may be true for small electrodes with high conductivities and low cur-

rent densities, long pipes exhibit a non-negligible potential difference along the

steel.3 7,40

Therefore, there are two separate domains for the flow of current. The first is

the soil domain up to the surfaces of the pipes and anodes. The boundary condi-

tions for the soil domain are the kinetics of the corrosion reactions as described

in section 1.2. The second domain is the internal pipe metal, anode metal and

connecting wires for the return path of the cathodic protection current.

2.1 Soil Domain

The corrosion problem discussed in Chapter 1 serves as the boundary condi-

tions for the soil domain which is governed by Laplace's equation. The equation

can be derived by imposing conservation of charge on Ohm's Law. Ohm's Law

is given by

i = XV7( (2-1)

where i is the current density vector, K is the conductivity of the domain and (

is the potential distribution through the domain. Conservation of charge in the

bulk yields

V i = 0 (2-2)

The result in terms of potential and constant conductivity is

V20 = 0 (2-3)

which is known as Laplace's equation.

The solution can also be obtained by performing a material balance over a

small volume element

Oci

= (V Ni) + R, (2-4)

Ot

where ci is the concentration of species i, Ni is the net flux vector for species i

and Ri is the rate of generation of species i due to homogeneous reactions. In

a dilute electrolyte, the flux for any species can be written as the sum of three

contributions: convection, diffusion and migration:

Ni = vc DiVci ziuiFciV( (2-5)

where v is the fluid velocity, Di is the diffusion coefficient for species i, zi is the

charge on species i, ui is the mobility, F is Faraday's constant, and ( is the po-

tential. The mobility is related to the diffusion coefficient by the Nerst-Einstein

equation

Di = RTui (2-6)

The current density i, that results from equations (2-4) and (2-5) can be expressed

as

i = F ziNi (2-7)

which is the charge flux within the continuum.

A derivation of the governing differential equation can be obtained by making

two assumptions concerning the electrolyte. Under the assumption of electro-

neutrality

zici = 0 (2-8)

i

and that heterogeneous reactions only take place at the boundaries, then one can

multiply equation (2-4) by Fzi and sum over all species.12 The final assumption is

that each of the species in the domain has a uniform concentration distribution.

This causes the diffusion term to vanish and the convection term to sum to zero

by equation (2-8)

V.( (F zici)) =0 (2-9)

The reaction terms also sum to zero because of charge conservation in a homoge-

neous reaction

F ziRi =0 (2-10)

i

The term involving the time derivative also sums to zero

SF zici = 0 (2-11)

This leaves a time independent equation

V (F z 2ci=VO 0 (2-12)

where the conductivity K

K = F2Y 2uici (2-13)

is a constant due to uniform concentration. Equation (2-12) then reduces to Laplace's

equation

V20 = 0

(2-14)

2.2 Justification of Dilute Approximation

The electrolyte must be dilute for equations (2-5) and (2-6) to hold. For soil

environments, the assumption that the soil is dilute can be justified by calcu-

lating the ionic strength. A soil with a low resistivity will be in the range of

p = 100Ocm. The ionic strength is defined to be

1

I= 2 zc (2-15)

The conductivity of an electrolyte is given in terms of mobility by

= F2 Z 2UiCi (2-16)

and in terms of diffusivities by

F2

= RT z2Dici (2-17)

The diffusivity of a typical species is105(cm2 /s). The conductivity is then written

as

1 FD 2ci (2-18)

p RT i

Substituting for the sum in equation (2-18), with equation (2-15), and accounting

for the porosity of the soil an expression for the ionic strength is obtained:

RT

I =- (2-19)

2F2Dp1.5

Using typical values of T = 290K, and D, and p taking values as above, a value

for the ionic strength is:

mol

I = 0.013 1 (2-20)

Since this value is closer to a maximum value in soil environments, the dilute

solution approximations used in equation (2-5) holds.

2.3 Internal Domain

The resistivity of pipe steel is 9.6 x 10 6Qcm. For an 18 inch diameter pipe,

the resistance per linear foot of pipe is about 2.3 x 10 -6. If 10 Amps must travel

5 miles on this pipe, the potential drop in the steel must be 600 mV which is a

significant potential drop. Therefore the assumption that there is no potential

drop in the pipe steel made by Esteban et al. cannot be made for long pipes.33 The

amount of pipe modeled by Esteban was approximately 40 ft.

The flow of current through the pipe steel, anodes and connecting wires is

also governed by the 3D Laplace equation

V (VV) = 0 (2-21)

where V is the departure of the potential of the metal from an uniform value

and K is the material conductivity. K does not necessarily have to be uniform.

For instance, the wire connecting the pipe the anode will be made of a different

material i.e., copper.

Equation (2-21) can be derived following the approach presented in section 2.1

if one treats electrons as a species.12 On a pipe surface, the analysis must account

for injection of current through bonds and through the pipe wall itself. If the

wires can be assumed to be well insulated, the potential drop though a wire can

be obtained by Ohm's law

AV = IR = IpL (2-22)

where R is the resistance of the wire, p is the electrical resistivity, L is the length

of wire, and A is the cross-sectional area.

2.4 Quasipotential Transformation

As described in section 2.1, use of Laplace's equation is predicated on the

assumption that concentrations of ionic species are uniform. This assumption,

which is generally valid, does not apply within a diffusion layer (about 50 pm

in thickness41) near the pipe and anode where electrochemical reactions generate

large surface concentrations of ionic species. The relaxation of the assumption of

a uniform concentration is possible within the context of the boundary element

formulation by use of the quasipotential formulation.42 43, 44, 45

The basic assumptions of the quasipotential transformation are that convec-

tion can be neglected, the system is at steady-state, and that concentration and

potential can be written as single valued functions of a common variable, qi.

Equation (2-5) therefore appears as the Nernst-Planck equation

i = -DiVci ziuiFci (2-23)

Since ci is a function of (,

dci

Vc = V (2-24)

d

thus

Ji = -Didi ziuiFci VDO (2-25)

Under the assumption that the diffusion coefficient is constant,

Ji = fi()VQ (2-26)

where

fi(_) = -Di ,Di zci (2-27)

( d(D RT

A new variable, obtained through the Kirchhoff transformation,43 has the prop-

erties

-Vqi = f;i()VD (2-28)

and

q j = f- I i()dD (2-29)

which allows the flux for species i to be written as

Ji = -Vqi (2-30)

The current density is given by the sum of contributions from the flux of each of

the charged species, i.e.,

i = F zili (2-31)

Equation (2-31) can be written in terms of the transformation variable as

i = -FFziVqi (2-32)

The quasipotential is defined to be

Q = F ziqi (2-33)

such that

i = -VQ = F ziJ, (2-34)

i

Under the assumption that charge is conserved, V i = 0, and Q is seen to satisfy

Laplace's equation

V2Q = 0 (2-35)

The quasipotential could be used to correct for local concentration gradients

by patching an inner solution of Laplace's equation for Q to an outer solution of

Laplace's equation for

the diffusion potential for the estimated concentration gradients, and this correc-

tion is typically less than 40 mV. The importance of the transformation is not so

much in the correction it gives to the calculated potential but the insight it leads to

changes in the local electrolyte chemistry caused by electrochemical reactions.46

The pH at the surface of the pipe was calculated by Allahar and Orazem to reach

values as high as 12 when the solution far from the pipe had a pH value of 7.46

Such dramatic changes in pH have significant influence on deposition of corro-

sion products or calcareous films.

2.5 Conclusions

Through the developments in sections 2.1 and 2.4, a complete model for cur-

rent transport through the soil for pipes under cathodic protection is obtained.

The quasipotential transformation accounts for the necessary effects of concen-

tration gradients near the surface of the pipes and anodes. However, it is desired

to not have to break the domain into two parts, one where convection is impor-

tant and one where it is not. Through experiments by Jeffers, the thickness of the

region that is governed by the quasipotential transformation is only a few mi-

crons. Therefore, because the parameters of the expressions for the kinetics are fit

to experimental data and because they contains expressions to account for trans-

port effects, the region governed by the quasipotential transformation is lumped

into the boundary condition.

CHAPTER 3

SOLUTION OF GOVERNING PDE

Laplace's equation has been solved for many boundary conditions and do-

mains.47 In the case of this research it is desired to be able to solve this equation for

arbitrary arrangements of pipes and anodes within the domain. It is also desired

to introduce the nonlinear boundary conditions that arise from the chemistry at

the boundaries as discussed in Chapter 1.

The first task is to find a method to solve the differential equation given ar-

bitrary boundary geometry. One way is to use complex transformations to map

the boundaries to a form for which a solution exists. Moulton first described the

technique for conduction of current through a solid.48 Bowman supplies details

for analytic solutions using the Schwaz-Christoffel transformation.49 Orazem and

Newman provide a semi-analytical implementation of the Schwarz-Christoffel

transformation for the more complicated structure of a slotted electrode.50 Orazem

used the technique for a compact tension fracture specimen51 which was later

further refined.52 Diem et al. take the semi-analytical technique a step further by

allowing for insulating surfaces to form an arbitrary angle with the electrode.53

Unfortunately, there are not enough of these transformations to cover all the

possible configurations of pipes within a domain. Also, the soil domain is not

bounded. The other way to solve the equation is to use a completely numerical

technique. Of the available techniques, the boundary element method is very

attractive. It can provide the same level of quality as the Schwarz-Christoffel

transformations while not being restricted to any geometry. The method only

solves the governing equation on the boundaries. This is ideal for corrosion prob-

lems since all the activity takes place at the boundaries. Brebbia first applied the

boundary element method for potential problems governed by Laplace's equa-

tion.54 Aoki and also Telles provided the first practical demonstrations of utilizing

the boundary element method with simple nonlinear boundary conditions.55'56

Zamani and Chuang demonstrated optimization of cathodic current through ad-

justment of anode location.5

The pipe steel domain is solved using the finite element method. Brichau

first demonstrated the technique of coupling a finite element solution for pipe

steel to a boundary element solution for the soil.58 He also demonstrated stray

current effects from electric railroad interference utilizing the same solution form-

ulation.59 However, Brichau's method was limited in that it assumed that the

potential and current distributions on the pipes and anodes were axisymmetric

allowing only axial variations. Since it is desired to have a solution for the current

and potential distributions both around the circumference and along the length

of the pipe, Brichau's method must be modified. Aoki presented a technique

similar to Brichau's that included optimization of anode locations and several

soil conductivity changes for the case of a single pipe with no angular variations

in potential and current distributions.60' 61

3.1 Boundary Element Method Development

The boundary element method can be derived from the same technique used

to obtain the classical finite element method. One starts by writing a variational

or weighted residual of the governing differential equation. If the PDE takes the

form

V2u f = 0 (3-1)

where f is a forcing term that is a function of position only, then one would write

the weighted residual as

j (V2u f) dQ = 0, Vw (3-2)

where Q is the domain and w is any weighting function. In the case of Laplace's

equation, f = 0 and u = (; thus

SwV27Id = 0,Vw (3-3)

This equation holds for all weighting functions w. Equation (3-3) is integrated

analytically by using the divergence theorem

SVwVcdQ I w( V7) dF = 0, Vw (3-4)

with F being the boundary of the domain. Equation (3-4) is classical weak form

of the finite element method for Laplace's equation. At this stage, the boundary

element method development departs from the finite element method by using

a second application of the divergence theorem. The highest order derivative is

thereby moved to the weighting function, i.e.,

/ OV2wdQ I w(i V)dF d + (D (i. Vw) dF = 0, Vw (3-5)

At this point the weighting function function needs to be specified to show why

the second application of the divergence theorem was done. If one observes that

the solution to the equation

V2Gi,j + i = 0

(3-6)

is the Greens function for Laplace's equation, one can simplify equation (3-5) by

picking the weighting function to be the Green's function.62 The first integral in

equation (3-5) becomes

/ D (-) df = -Qi (3-7)

Substitution of equation(3-7) into equation (3-5) yields a simpler equation valid

within the domain

iD + (. V Gij) dF = Gi, (n VD) dF (3-8)

where the highest order derivative has been removed and all the integrals are

along the boundaries only. Equation (3-8) is still exact in so far as the Green's

function is known and is valid for determining the values of the potential at any

interior point in the domain given that the potential and current density distri-

butions on the boundary are known. This implies that, for Laplace's equation,

values of interior points are fully specified by integrals along the boundary only.

The last step in deriving the Boundary Element method is to take ci to the

boundary. A Cauchy principle value is introduced in the integral on the left side

of equation (3-8). It is usually represented by a constant appearing in front of the

first term in equation (3-8)

Ci(;+ I D (Wn VGi,j) dF = Gj (W V)D) dF (3-9)

For a smooth surface at the point i, the constant, ci, is equal to 7r.63

3.2 Symmetric Galerkin Boundary Element Method

The direct Boundary Element Method as described above results in dense,

non-symmetric matrices. Even though the diagonals of the matrices have the

largest magnitude, they are not diagonally dominant. This means that the fast

iterative techniques of solving Ax = b cannot be used. A form of the boundary

element method that is symmetric has been developed.64 Symmetric formula-

tions are possible because both the Greens function for Laplace's equation and

its normal derivatives are symmetric. The method can be derived as follows:65

start with the standard boundary element formulation (equation (3-8)) and take

the normal derivative of it with respect to the source point (point i)

di V4i + j D ((i V (, VGi,j)) dF = j (i VGi,j) (y V(D) dF (3-10)

which is known as the hyper-singular boundary integral formulation. It is called

hyper-singular because the highest order singularity is greater than two. The

next step is to multiply through equations (3-8) and (3-10) by the shape function,

0, of the source element, i, and integrate the resulting equations over the entire

boundary

j icidFi+ I Of Ij D (-. VGi,) dFjdFi= / I Gij (. VD) drjdFi (3-11)

and

jK (d VcD,) dF, + I ji Ij (n, V (t VGi,j)) dFjdF=

J i j (di VGi,j) j ( V1) dFjdFri (3-12)

The additional outer integration now results in a symmetric matrix on the left

side. Equation (3-11) is written on Dirichlet surfaces and equation (3-12) is written

for Neumann surfaces. The symmetry of the left-hand-side matrix is only seen

when all the unknown variables have been moved to the left hand side. The

integration can be more easily seen in Figure 3-1 on page 34.

Figure 3-1: Symmetric integration. Lines between the two circles represent the

Green's function relation between the two integration points on the surfaces dFi

and dFj .

The symmetric positive definite system of equations resulting from this method

can now utilize Cholesky decomposition for their solution. This method is signif-

icantly faster than LU decomposition which is needed for the standard boundary

element method. However, the symmetric Galerkin boundary element method

can only be used when the boundary conditions are constant. The symmetry falls

apart for nonlinear boundary conditions as demonstrated in Chapter 4.

3.3 Boundary Conditions

The form of the boundary conditions for steel pipe under cathodic protection

have been developed in Chapter 1. They take three possible forms one for bare

steel, one for coated steel and one for anodes.

3.3.1 Bare Steel

For the static problem, the boundary condition relates the potential difference

between the steel and the soil just outside the outer Helmholtz plane, V (, to

the current density, i = nR V7 at the soil pipe interface. The expression for the

corrosion kinetics developed in Chapter 1, equation (1-23), is used for sections

of the buried structure where bare steel is exposed to the soil (coating holiday or

uncoated pipe).

3.3.2 Coated Steel

Coated steel has a transport barrier in the form of the coating. The coating is

also bonded to the steel and therefore many of the reaction sites are blocked. The

equation developed in Chapter 1 also relates the net current density at the pipe

- soil interface to the voltage difference between the pipe steel at the interface

and the potential in the soil just above the interface (past the outer Helmholtz

plane). The expression, equation (1-26), also introduces an inner potential which

is defined to be the potential at the coating steel interface. This is the potential

that is used to drive the kinetics.

3.3.3 Anodes

The boundary conditions for the anodes must be an essential type boundary

condition in order to have a well-posed problem. An essential type boundary

condition is one imposed on the potential at the boundary. A condition imposed

on the normal flux at the boundary would be referred to as a natural boundary

condition. However, the equations describe the reactions on the anode are also

nonlinear functions of the potential. Equation (1-27) can be quickly solved for the

potential

V ( = Eeql + anode log(i/i02 + 1)

(3-13)

which can then serve as an essential boundary condition. Many cathodic pro-

tection systems may not have a galvanic anode but an impressed current anode.

The equation for the current density for impressed current anodes, equation (1-

30) can also be rewritten in where i is the independent variable

V D = Eo2 + Vrectifier + 02 log (i/ieql + 1) (3-14)

Equations (3-13) and (3-14) are then used to form a well-posed problem definition

for the solution the governing differential equations.

3.4 Infinite Domains

Everything done to this point has been done under the assumption that the

boundary encloses the domain. In many situations, the domain lies outside the

boundary. In order to get a solution for this situation, it is necessary to introduce

a second boundary, T, placed around the surface of interest and centered on the

source point on that surface. Adding the enclosing boundary to equation (3-9),

one obtains

Ci;i + c (. VGi,j) dF +c I D (i. VGi,j) d =

SGi,j (nd V(D) dF + I Gi,j (n VcD) dT (3-15)

If the radius of the new boundary is taken to infinity the limit of the integrals

over the external boundary are taken to infinity

lim ( W C V Gij) dT Gi, (d V7D) dT (3-16)

When integrating across this boundary the integral will be

D ( VGi,y) dFO = Hi,, = 47r0 (3-17)

Jl00

where the minus sign indicates the direction of the normal vector. The other

integral vanishes as the limit is taken. The value of integral in (3-17) can be more

easily seen if it is transformed to spherical coordinates

lim (r2 sin ) dOd = 47r (3-18)

where r is also the normal vector since the surface is a sphere centered at a source

point. ( at infinity is often assumed to be zero; thus satisfying the zero radiation

condition at infinity exactly.

3.5 Half Spaces

A half space is simply an infinite domain split by a plane. The half space is the

space lying on one side of the plane. If either the Dirichlet or Neumann condition

vanishes at the plane, the Green's function presents an interesting opportunity

to satisfy that condition exactly with a very small additional computation when

evaluating the kernels of the integrals. This is done by making use of the Green's

functions reflection properties.63 64, 65, 66 In the case of buried pipelines, the Neu-

mann condition vanishes at the plane boundary, that is there is no current flowing

out of the soil into the air, and none is flowing from the air into the soil.

If one starts with the boundary condition on the plane being zero normal cur-

rent and places a source at x and its reflection about the plane at x', one may

write66

a (x) ( VG (x, )) + a (x') (n. VG (x', )) = 0 (3-19)

which implies that the two source intensities o(x), and a(x') are equal and have

the same sign. The sign is the same because the outward normal vectors have

opposite signs. The final form of the Green's function is

1 1

Gi = (+ (3-20)

Gi 47rr(xi, xi) 47r(xi, x) )

where x. is the reflected source point.

3.6 Layers

Layers are created using the same types of Green's function reflections as de-

scribed above. The only restriction being that one of the two boundary condi-

tions at the interface between the two layers must be zero (either ( = 0 V or

n V( = OiA/cm2).66 In the context of cathodic protection of buried pipelines,

the only boundary condition that has physical meaning is a zero normal current

condition since there is no easy way to have an arbitrary plane within the elec-

trolyte that has a potential equal to zero. Including a plane with a zero Neuman

(natural) condition implies that there is one region in which current may flow

that is bounded by regions of zero conductivity whose boundaries are defined by

the Green's function reflections. An example would be an underlying rock layer

which has zero conductivity.

The resulting functions can be obtained by using equations of the form of (3-

19). The Green's function in three dimensions would then be of the form

Reflections+1l 1

Gij = 4 (3-21)

k 47Fr(xi,Xj ,k)

where the index k > 1 refers to the reflection about some plane k. For k = 1, Xj,k is

the field point on the real object. If the soil surface is the only reflection used, the

result is the same as equation (3-20).

Reflection

Figure 3-2: Reflections of Green's function to account for boundaries with zero

normal current.

An example of a pipeline in a halfspace with an underlying rock layer is

shown in Figure 3-2 on page 39. There are two reflections used. The first is to

account for the zero normal current at the air soil interface. This is represented

by the dashed pipe in the air. The second reflection is to account for the zero

normal current at the rock soil interface and is represented by the lower dashed

pipe.

A model was run where the non-conducting layer was placed at a couple of

inches below the pipe and slowly moved away from the pipe. The total cathodic

current that the pipe received from the anode was calculated with no other pa-

rameters changed except for the depth of the non-conducting layer. A plot of the

total current is given in Figure 3-4.

The influence of the nonconducting layer on the current delivered to a pipe

is seen in Figure 3-4. For these calculations, the anode was placed far from the

Air

Soil

Rock

Reflection

I

Figure 3-3: Rectangular enclosure

reflections.

o 0.21

E

5 0.19

- 0.17

0

I-

0.15

created with a primary Green's function and 6

0 50 100 150

Depth of zero conductivity (ft)

200

Figure 3-4: Total cathodic current driven to a pipe with an underlying non-

conducting region.

pipe such that the current distribution around the circumference of the pipe was

affected only by the screening of the insulating ledge. The total current tended

toward that predicted by Dwight's formula as the distance between the ledge and

the pipe increased. As seen in Figure 3-5, the current increased with the distance

according to

I = Ir- (3-22)

Vr

where a is a fitted parameter that depends on the geometry and soil resistance

and Ir,, is the value obtained from Dwight's formula which does not account

for the presence of a nonconducting layer.

The data used in Figure 3-4 was then then fit to the equation

i = + ir=o (3-23)

Vr

where a is a fit parameter that contains the area of the pipe, the soil resistivity,

the applied voltage of the anode, etc. After calculating a linear fit, the data can be

plotted again along with the line representing the linear fit.

3.7 Development of Finite Element Method for Internal Domain

The finite element method has been selected for the domain of the pipe steel,

copper connection wires, and internal anode material. It is ideal for completely

bounded domains where the material properties change from one location to an-

other, i.e., current flow from the pipe to the copper wire connecting the pipe to

the anode.

The development of the finite element method for the pipe steel domain starts

by writing a weighted residual for the strong form of Laplace's equation given in

0.25

0.2

S0.15

S0.1 i = -0.064/sqrt(r) + 0.2118

0 R2 = 0.9789

0.05

0

F-

0 I I I I

0 0.2 0.4 0.6 0.8 1

1/sqrt(depth)

Figure 3-5: Fit of the data from Figure 3-4 to i 0.064 +0.2118.

equation (2-21)

SwV - VVdQ = 0 (3-24)

where w is any weighting function. For pipe steel, is just a constant times an

identity matrix

100

= K 0 1 0 (3-25)

001

with step changes in when the type of metal changes. After substituting the

material properties into equation (3-24) and integrating by parts using the diver-

gence theorem, the following Weak Form is obtained

J n Vw. VVdxdydz= -f #v(i. VV)ds (3-26)

where n is the outward normal vector from the boundary of the domain. The do-

main is divided into elements using piecewise continuous polynomial isopara-

Figure 3-6: Diagram of the shell element used to calculate the potential drop

within the pipe steel. The element is assumed to have no variation in the ( direc-

tion. The curvilinear coordinate system is displayed on top of the element.

metric shape functions. The shape functions approximate the geometry and so-

lution, V, over the elements. The discritized form of equation (3-26) is

V ) [/ I e) e)dx dy dz] V = (f. VV) e) ds (3-27)

where K is the scalar component of the property tensor, Oi is a shape function and

the sum goes from 1 to the total of all the shape functions of all the elements.

3.7.1 Pipe Shell Elements

A special type of thin shell elements is introduced here and shown in Figure

3-6. They are specifically designed for potential problems on shells where the

absolute value of the material property (K in equation (3-27)) is large.

The elements are defined in orthonormal curvilinear coordinates, J, r, and (

where J and r define the outside surface of the pipe and ( is the outward normal.

( is obtained through the cross product J x r.

Variations of the potential in the ( direction are assumed to be negligible be-

cause the scale of the problem is many orders of magnitude greater in the J and

directions. Variations of the potential parallel to the surface are allowed to vary in

a piece-wise continuous way using bi-quadratic shape functions for the elements.

These functions are obtained through the product of two Lagrange interpolating

polynomials of the same order, one in J and one in q. The result is a family of el-

ements with square parent elements. They include the four-node linear element,

the nine-node quadratic element, the 16 node cubic element, etc.

3.7.2 Applying Elements to the FEM

The integral (3-27) must be transformed to the curvilinear coordinate system

of the parent elements. A differential volume, dx dy dz can be transformed to the

curvilinear system by the determinate of the Jacobian of the coordinate transfor-

mation

dxdydz = J d drld( (3-28)

where J is the Jacobian.

The Jacobian is composed of the partial derivatives of the coordinates x, y and

z with respect to each of the curvilinear coordinates

Ox ay Oz

ap ap ap

j= x ay az (3-29)

a7 a7- a5

ax ay Oz

a ac ac -

Integral (3-27) transformed to the curvilinear coordinate system is

ifi fe i |xlJdj dridCd = wK(W. V V)ds (3-30)

k JXk 9Xk J = d #F7e)

where Xk is one of the cartesian coordinates, IJI is the determinate of the Jacobian,

s is the surface of the element n is the outward normal vector from the surface,

and Q(e) is the domain of the parent element. The limits of integration in the

parent element are from -1 to +1 for all three of the coordinates. For the special

shell elements used for pipes, the integral is only performed numerically over

the surface dd drj which is then scaled by the physical thickness of the shell h, the

result of the integral over (. This requires the Jacobian to be modified to a surface

Jacobian. In this case it is simply the square root of the magnitude of the Jacobian.

Rewriting equation (3-30), the two dimensional integral over the surfaces of the

pipes is obtained

Sh 0i dd =_0 wn (. V V) ds (3-31)

with the thickness of the steel, h, a parameter. The normal vector has the same

direction as the one from the boundary element method.

Partial derivatives in equation (3-31) with respect to cartesian coordinates

need to be expressed in terms of the curvilinear coordinates. Using the chain

rule, but starting from the derivatives in cartesian coordinates, one writes

Ox 9^- + 9y 9, 9z 9^ 9

a 7x + 9y 7, 9z = aq a (3-32)

9Ox a a Oy ai Oz ( a9

or, in terms of the Jacobian

a 99 = 190 (3-33)

Since 0 is a function of J, q, and (, the partial derivatives of b with respect to the

cartesian coordinates can be found by inverting the Jacobian.

ao = J-1 >g (3-34)

7Y 717

If J1 is the element from the ith row and the jth column of the inverse of the

Jacobian, then the needed partial derivatives can be expanded to

= ( )11 + (j1)12 +(1)13

= (J1)21 + (j1)22 (J1)23

= (J)31i + (J )32 (J1)33 (3-35)

The assumption that the thickness of the shell is small with respect to the ra-

dius and length of the shell means that there is negligible variation in the poten-

tial in the ( direction. Therefore, all terms in equation (3-35) that involve partial

derivatives with respect to ( can be assumed to be equal to zero. All of the terms

in equation (3-31) can be evaluated numerically. Since the shape functions for

the geometry and solution vary quadratically in the J and qr directions, a nine-

point Gauss rule (3x3) applied in both directions will give an exact (to machine

precision) result.

3.8 Bonds and Resistors

Bonds and resistors are used to electrically tie two pipes or a pipe and an an-

ode together. They are implemented as a linear 1-D element that goes between the

connection node on one pipe to the connection node on the second pipe. There-

fore, no extra nodes are introduced. The material property is set by accounting

for the real length and gauge of the wire that ties the pipes together. If a resistor

is specified within the wire, it is added to the total resistance of the bond. An

illustration of a bond is given in Figure 3-7. The red cylinders represent a resistor

Figure 3-7: Black lines connecting pipes and anode are bonds. The red cylinders

represent resistors. They are modeled as 1-D finite elements that connect to the

bond connection points.

that is sometimes added in series with the bond. If one is not used, its value of

resistance is set to zero.

For a 1-D element, a one-dimensional version of Laplace's is used. The devel-

opment for the linear 1-D element is done in the same way as in section 3.7. The

result of the integration over the bond is the total conductance of the bond or

Ke = R (3-36)

-1 1

where R is the total resistance of the bond.

3.9 Coupling BEM to FEM

The finite element domain is coupled to the boundary element domain at the

interface between the two domains. Ohm's Law holds within each domain as

-~4~5~1r+

stated in section 2.1. Therefore, at any arbitrary surface that forms a boundary

between two domains it can be shown that the flux on either side of the boundary

is related by the material property.

K1nii* VI1 = K2ni* V(I2 (3-37)

which is simply a balance on charge at the pipe/soil interface, where K is the

material property. For potential problems, K is the conductivity of the material in

Mho/m. Using the variables for potential in the pipe, V, and potential in the soil,

(, the interface condition is written as

( V V) = 'soil (n VID) (3-38)

Steel

where the quantity n VV is used in the finite element load integral.

Equation (3-38) can be inserted into equation (3-31) to obtain

h /.- I1ddrlq] _f-e oil (. Vc) ds (3-39)

The right side of equation (3-39) links the soil domain to the steel domain through

the current density generated by the kinetics of the corrosion reaction at the steel

surface given by equations (1-23) and (1-26).

3.10 Current Flow in the Pipe

In addition to using false-color surfaces to indicate potential or current density

normal to the surface, Figure 3-7 includes vectors laid on the surface of the pipe

which represent the direction and relative magnitude of the current flowing in the

pipe steel. Calculation of the current density within the pipeline steel can be used

to isolate regions of inadequate protection, to identify the electrical midpoint, and

to determine the direction of current flowing through a bond. The method for

calculating the current within the pipes is presented in this section.

The value of the two components of the current in the parent element at each

of the four optimal Gauss points is found by differentiating the shape functions

which describe the potential in the steel, i.e.,

S )= Vk (3-40)

k= 1

and

i, = Vk-0k (3-41)

k=1

Since the current vector is needed in cartesian coordinates, but the currents are

only known in the curvilinear system as above, the Jacobian of the coordinate

transformation must be used to get the three partial derivatives in cartesian co-

ordinates (see equations (3-32) and (3-34)). The values of the current are found at

each of the optimal sampling points through

DVI

ix,i = -Ox = (l I, I) = (U 1)12 I, I) (3-42)

for the x component,

DVI

iy,I (" = = (121 I, rI) + (J 1)22i/( i, r I) (3-43)

for the y component, and

iz,i = i = (J 1)3li(I, rI) + (J 1)32ij(I, II) (3-44)

for the z component, where (Ji, rI) are the optimal current sampling points and

i and i, are defined in equations (3-40) and (3-41), respectively.

The values of the flux at the nodes are found by fitting three planes to the

three components of the four optimal values of the current by the equation

ix(, rl) = c1 + C2 + C3r] (3-45)

where the constants are found through a least squares regression. Therefore, a

transformation matrix is needed to perform the regression

Inodes = [TR]ioptimal (3-46)

where optimal is a 4-by-3 matrix and inodes is a 9-by-3 matrix, hence [TR] is a 9-by-4

regression matrix. To find [TR], the error of fitting a plane to the current at each

node is minimized through the sum of squares

IV

ex = [c1 + C2k C3k ik]2 (3-47)

k=I

where ex is the error in the x component of i. Two more equations of the form

of equation (3-47) are written for the y and z components of i. The error is mini-

mized for each component by setting

Oex

S= 0 (3-48)

Ocj

which yields

IV

[C1 + C2k +C3rk ik] = 0 (3-49)

k=I

for j = 1

IV

[C1 + C2& + C3k k] = 0 (3-50)

k=I

for j = 2, and

IV

[Cl + C2k +C3rk -k] = 0 (3-51)

k=I

for j = 3. Equations (3-49) through (3-51) can be written in matrix form as

k=IV y k=IVr

Lk=I k Lk=I Ik

k- IV k k ( k-IV

k= IV k=IV

k=I k Lk=I r lkk

1 1

C2 I

\C3 / I

1 1 1

rII III I zv

7 il /III TIV

or in matrix notation as

[P]{c} = [Q]}optimal

Then c for each component of i can be written as

{C} = [P] [Q]{}optimal

For the optimal Gauss points, { }, [P] is

4 0 0

[P]= 0 0 0

0 0

and [Q] is

[Q]=

1 1 1 1

1 1 1 1

3 3 31 1 3

1 1 1 1

3 3 -3 3

so that [P] [Q] is

1

1

[P] I[Q] = 4 -

3 f

1 1 1

yk=IV 1

k=I I

yk-IV k

k IV

L k IV rk

im

iii

(3-52)

(3-53)

(3-54)

(3-55)

(3-56)

(3-57)

1

The 9 nodal values of the flux are then found from equation (3-45) and written

in matrix form as

ii 1 -1 -1

i2 1 0 -1

i3 1 1 -1

i4 1 -1 0 c

i = 1 0 0 C2 (3-58)

i6 1 1 0 kc3

i7 1 -1 1

is 1 0 1

ig L1 1 1

Upon substituting the values of cj from equation (3-54), the transformation matrix

becomes

1+23 1 1 1-2-3

1+ \3 1+ 3 1-

1 1-2 1+2 3 1

1+ 3 1+3 1- 1-3

1

[TR] = 1 1 1 1 (3-59)

4

1-V3 1 3 1+~ 1+ 3

1 1+2-3 1-2-3 1

1-3 1+3 1-V 1+ 3

1-2-3 1 1 1+23,f

Thus, the three components of the flux are found through

{ix}nodes = [TR]{ix}optimal

(3-60)

Figure 3-8: Normalized current vectors drawn on top of the pipe. Three compo-

nents of the vectors are calculated in cartesian coordinates even though the flow

of current is only in two directions. The fact that the vectors are always tangent

to the pipe indicates that the coordinate transformation was done correctly.

{iy}nodes = [TR]iy optimal

{iz}nodes = [TR]{z}optimal

(3-61)

(3-62)

The production code normalizes the current vectors by the largest current vector

when they are drawn to the screen. An example of the vectors that are calculated

with the above method is shown in Figure 3-8. The concepts presented in this

section represent an extension of the methods described by Burnett.67

and

3.11 Conclusions

This chapter has presented the numerical method for solution of the govern-

ing differential equations within the soil domain and within the steel domain.

Special elements have been introduced to simplify the solution for the potential

distribution within the pipe steel. Techniques have also been introduced to deal

with the soil surface and underlying rock ledge without the need to discretize

those boundaries. The reason is that both boundaries have a normal current den-

sity equal to zero, which can be explicitly handled by the Green's function for

Laplace's equation.

The boundary conditions for the differential equations and the coupling be-

tween the steel and soil domains have been briefly mentioned for completeness.

Details of the application of the boundary conditions to the numerical methods

and coupling of the two domains are contained in the next chapter.

CHAPTER 4

APPLYING KINETICS TO THE BEM AND FEM

A numerical method has been developed in Chapter 3 for solution of the gov-

erning differential equations for pipelines under cathodic protection. Now, the

soil domain boundaries must be divide into elements and the pipe steel domain

must be divided into finite volumes. The volumes will use the elements described

in section 3.7.1, so that the loads (current density) at the boundary with the soil

domain will match what is calculated in the soil domain.

While the boundary conditions have been mentioned in section 3.3, it is not

been explained how they are applied to the numerical method. This chapter cov-

ers the application of the boundary condition to the numerical methods in detail.

4.1 Pipe Discretization

New pipelines are often placed very close to old pipes because there are a lim-

ited number of right-of-way corridors. The need to place pipelines in close prox-

imity means there is a strong possibility of pipes interfering with the protective

current being supplied to other pipes. To account adequately for any variations

in the current and potential distributions, the pipes need to be discritized in both

the angular and axial directions. This allows any type of distributions to be ap-

proximated by the basis functions. It also makes the matrices generated by the

BEM very large.

4.1.1 Discretization of the Boundary Element Method

Equation (3-9) may now be applied to a surface that has its boundaries broken

up into individual finite elements (see Chapter 5 for details of the mesh genera-

tion). The solution to the differential equation can be represented by approxi-

mate functions on each individual element. Example functions may be Oth order

or constant value, 1st order or linear, etc. The simplest case is the constant ele-

ment. The value of ( is assumed to be constant across each element. Then one

can write an equation of the form of (3-9) for each degree of freedom that makes

up the surface and where the subscript i refers to the element number. The inte-

grals are broken up into N sub-intervals corresponding to N elements that when

summed together form the complete integral over the entire boundary

cii +, I D (nd VGi,j) dFj = I j Gi,j (n VD) dFj (4-1)

j J j J

When constant element are used, the unknown element potentials and normal

current densities can be brought outside the integrals. Since a well posed problem

has one half of the boundary conditions specified, the result is N equations with

N unknowns.

Since it is desired to be able to change the boundary condition type as well as

value, the notation used here will denote the matrix resulting from the left hand

integral as H and the right side as G

HQ = G (n. VD) (4-2)

where n. Vc is the normal current density divided by the conductivity at the

boundary.

The diagonal elements of the matrix H are often unknown at this point be-

cause linear or higher order elements have been used and the constant ci in equa-

tion (4-1) is know only as a Cauchy principle value. A theorem by Gibbs can be

used to easily find the values of the diagonal.68 The theorem states that for any

domain governed by Laplace's equation, if the potential is the same everywhere

throughout the domain and that the gradient of the potential at any boundary

of the domain including infinity if it exists is equal to zero, then the gradient of

the potential is equal to zero everywhere within the domain. The theorem is im-

plemented by setting a uniform potential equal to one throughout the domain.

Then the gradients of potential are known to be zero everywhere which leads to

the matrix equation

H1,1H1,2 1 0 (4

(4-3)

H2,1H2,2 1 0

where the diagonal terms are unknown. The diagonal terms are now specified

by the negative sum of the off-diagonal terms of each row

Hi,i = Hi, (4-4)

j7i

where j goes from one to the number of terms in a row.

When the known boundary conditions are applied, columns are swapped be-

tween the two matrices such that all the known variables are on the right side.

If the resulting matrix is multiplied against the known boundary conditions, a

standard linear system is obtained

Ax = b (4-5)

The solution to this system of equations yields the necessary unknown boundary

conditions. Following the solution on the boundary, the value of the potential

may be obtained at any point within the domain by use of equation (3-8). In this

case ci is the unknown potential inside the domain, and (, and (n. Vc) are the

solutions and known boundary conditions found previously.

If one wishes to obtain current densities within the domain, a little more work

is needed. The gradient of equation (3-8) is taken to get the components of the

current density vector within the domain

VDi + f (DV ( *. VGi,j) dF = I ( VG, ( VD) dF (4-6)

One immediately notices that the nabla operator has been taken inside the inte-

grals. This can be done because the limits of the integral are not affected. The

first integral in equation (4-6) does not vanish. The reason for this is that the nor-

mal derivative of the Green's function has been taken at the boundary and so one

does not obtain the Laplacian of the Green's function.

Equation (4-6) can be split into 3 separate equations, one for each component

of the current density vector. This equation is commonly referred to as the hyper-

singular boundary element method because it contains a 1/r3 singularity. This

singularity never appears if the interior point chosen is sufficiently far from the

boundary.

4.1.2 Discretization of the Finite Element Domain

The finite element domain is discretized using the same surface mesh that is

used in the boundary element method. This means there is a one-to-one corre-

spondence between the degrees of freedom for the boundary element method

in the soil domain and the finite element method in the pipe/anode interior do-

main. It also assumes that the diffusion length at the pipe surface is negligibly

small compared to the size of a pipe or holiday. Therefore, it can be assumed that

the point at which ( is evaluated in the Boundary Element Method, is approxi-

mately the spot where V is evaluated and the same nodal locations can be used

in both methods.

The match-up is necessary to simplify the evaluations of the boundary condi-

tions, equations (1-23) and (1-26). The boundary conditions use both the solution

from the finite element method in the steel and the solution from the boundary

element method in the soil. The difference between the two solutions, V (,

provides the driving force for the reactions.

4.2 Self-Equilibration

Cathodic protection systems by nature, do not loose or gain electrons from

their surroundings. In other words, charge is conserved. The boundary ele-

ment formulation as developed in Chapter 3 will not show this feature without

some modification. Instead of having a specified potential at infinity (most often

set to zero) which serves as an infinite source or sink for charge, an unknown

value of potential at infinity is used such that no current enters or leaves through

that boundary. To implement this condition, an extra equation is added to the

system56,69

Y Kf VdFi = 0 (4-7)

Equation (4-7) simply states that there is no current lost to or gained from infinity.

The left hand side is added as a new row at the bottom of the G matrix of equation

(4-2), while the matrix H in (4-2) receives a row of zeros. A column is added

to the H matrix which corresponds to the unknown potential at infinity. The

values that are placed in this column come from equation (3-17) and are all 47r or

1 depending on where the 47 from the Green's function is placed. This would

make the H matrix singular because there is still a row of zeros in it. However,

that would only be the case if Neumann type boundary conditions are specified

everywhere. The Neumann problem results in an infinite number of solutions

that differ by a constant. Therefore, at least one element in the system must have

a Dirichlet boundary condition to make the H matrix nonsingular and result in a

unique solution.

4.3 Multiple CP Systems

Riemer and Orazem describe in detail the development needed to model in-

teractions among CP systems.32 To allow stray current to occur between separate

CP systems, additional rows of the type in equation (4-7) can be added.70 For

each separate CP system, one extra column in the H matrix is also added. The

appearance of the matrices for 2 CP systems in the same domain will be

G1,1 G1,2 G1,3 G1,4

G G2,1 G2,2 G2,3 G2,4 (4-8)

G= (4-8)

0 0 A3 A4

A1 A2 0 0

H1,1 H1,2 0 1

H H2,1 H2,2 1 0

H= (4-9)

0 0 00

0 0 00

with the column matrix u given by

U1

U2

u = (4-10)

Uoo,System2

Uoo,Systeml

and the column matrix q defined in the usual way. When there are two or more

CP systems within a domain, a Dirichlet condition on at least one element in each

CP system must be specified to prevent the H matrix from being singular. The

added equations and unknown potentials at infinity are sufficient to allow several

CP system to interact with each other while enforcing that the total current on

each CP system sum to zero.

4.4 Coating Holidays

A coating holiday is a small region of pipe where the coating has somehow

been removed. This can happen during handling, back fill, or digging in the

region of the pipe. In this way, bare steel is exposed to the corrosive soil envi-

ronment. Without cathodic protection, such exposed areas corrode very rapidly.

The design equations of Morgan are incapable of handling nonuniform poten-

tial and current distribution in the angular direction. Therefore another solution

method is needed. Early attempts at solving this problem used finite elements.24

While the method used had a 15% error in the sum of all current being zero, it did

show that the majority of current entered the pipe through the holiday. Bound-

ary element techniques fared much better at correctly calculating the current.33

Kennelly et al. used constant triangular elements to represent the potential and

current distributions.24 This meant a great many triangles were needed where

the discontinuities in the boundary conditions occurred in order to represent the

sharp profile of the solution in these regions. Therefore, only small sections of

pipe (about 20 feet) could be modeled successfully when holidays were to be

included.

To include holidays in the boundary element model, one must have a group

of elements whose outside boundary form the edge of the holiday. Since it has

been stated that the discontinuity of the solution at the holiday/coating boundary

results in a very sharp curve in the solution, there must be sufficient elements in

the neighborhood of the boundary condition change such that the solution can

be adequately represented by their basis functions.

4.5 Nonlinear Boundary Conditions without Attenuation in the Pipe Steel

The boundary conditions described in Chapter 1 are nonlinear and require

special methods to enforce them in the boundary element method. Aoki et al.

describe a method to address nonlinear boundary conditions.71

First, it is assumed that V is zero everywhere, in all of the pipes and anodes.

This assumption is justified if the conductivity of the steel of the pipes is several

orders of magnitude greater than that of the soil environment and the pipes are

short. The pipe length is constrained by the requirement that there be no ap-

preciable potential drop from one end of the pipe to the other as compared to

the smallest potential drop through the soil between any two points between an

anode and pipe.

Then, the procedure is to do a Taylor series expansion of equation (4-2) around

the known boundary conditions. A Newton-Raphson iterative procedure is ob-

tained. One starts by writing a residual after the columns corresponding to the

unknown boundary conditions have been moved into the matrix H

RResidual = H ......... G ......... (4-11)

S-W. V o = f (()

The residual results because only a nonlinear relationship is know between (

and n. Vc so ( must be guessed in order to obtain the known boundary condi-

tion n. Vc. Therefore equation (4-2) does not hold. For exposed steel in a soil

environment under cathodic protection, that relationship may be of the form of

equation (1-23). The Taylor series expansion truncated at the linear term gives

Y (k) r A Ik)

O=H ... + ......

(n0. VcI)) A(n. VcI) k)

G ......... + ............... (4-12)

SA(n.VJ))

The residual error, RResidual, is obtained by subtracting equation (4-11) from equa-

tion (4-12)

0 : 0 AO)

RResidual = -H + G .................. ......... (4-13)

I a (u V(D))k)

The matrix containing L-fi f (c(k)) /9Oc] is an m x n diagonal matrix which is

obtained from the polarization curves (see Appendix A). Since G is an n x m

matrix, the resultant that is added to the H matrix is square. The new values of (

and n. Vc are obtained from:

(k+1) = (k) + A(k)

(ft V 4))k+l) Y 1())k) k) (4-14)

(nI-V =(n Iv +A(W-V^) (4-14)

The new values can be substituted back into equation (4-11) to obtain a new resid-

ual. The process is repeated until the residual is sufficiently small.

It should be noted that the above method will often fail if the initial guess is

not close enough to the solution. Therefore a line search method is needed to

insure convergence.72 The method works by defining an objective function that

must be reduced at each iteration of the technique

g:= R2+...+R (4-15)

A subroutine takes the column matrix containing the Newton-Raphson step and

updates the unknown variables according the the function

x = Xod + AAx (4-16)

where x is the column matrix of unknown variables (( and n V() and A is a

number with a value 0 < A < 1. The method is illustrated schematically in Figure

4-1. The initial calculation is done with the full Newton step (A = 1). If the sum

of squares or the residual error, g, obtained from equation (4-15) is reduced, the

updated values of x are accepted. If g is not reduced, a smaller value of A is

selected using the same Newton step and a new residual and objective function

Cannot

converge

Yes

Done

Figure 4-1: Newton-Raphson line search flow diagram.

calculated. This is repeated until the objective function is reduced from that of

the previous iteration. Then the boundary conditions are recalculated and a new

Jacobian is calculated and the method is repeated until the residual, g is close

enough to zero to stop. At any point in the method,if a value of A cannot be

found such that g is reduced, the method fails.

4.6 Nonlinear Boundary Conditions with Attenuation in the Pipe Steel

The external load integral from the finite element method can be written for

each node as a sum of contributions from each element that has that node in

common

nce

fi = I f ,isoil (n V4) ds (4-17)

= 1 (e)

where fi is the load for node i and nce is the number of contributing elements

to the load at node i. This integral must be evaluated across all nodes on all

surfaces of all structures in the model. Since the value of n- Vc is either a known

or unknown boundary condition from the boundary element method, and the

kinetics of the electrochemistry, it must be factored out of the integral. Using the

shape functions that describe the values of n Vc across each element in terms

of the nodal values n. Vc and the weights, one can rewrite equation (4-17) as

nce nen

fi = I '., Ksoi I ,k(s) (W" V7D),k ds (4-18)

e=1 (e) k= 1

where nen is the number of nodes in the element, (n. VI))k is the nodal value

of the normal electric field within the element n, and O(s)k is the shape function

whose value is one at node k. Using the property of integrals that an integral of a

sum is the sum of integrals, equation (4-18) can be rewritten as

nee nen

fi = 11- Ie, ,-...,. ., (s)(n. V))t,k ds (4-19)

= k 1 k=1 (e)

Equation (4-19) can then be written in matrix form as

(n. VID)1,1

fi= [..- weiKs e,k(s)ds ... (4-20)

L J (e) I

(if- V<)ne,nen

or

fi = i [n. VID] (4-21)

These sub-matrices are then assembled with the right hand column matrix from

equation (4-2) to form the global finite element system

K[V] = F[. VcI] (4-22)

where F is a matrix formed from the assembly of equation (4-21) such that the

column matrix corresponding to the Boundary Element Method is formed on the

right hand side. F is not symmetric.

4.7 Solution by Successive Substitutions

The first attempt to couple the FEM to the BEM, was to use a successive sub-

stitutions approach. The idea was to assume an initial value of the voltage in

the pipe steel to be zero everywhere. Then the current density distribution was

calculated using the BEM. The result was the values of the load vector, [n. V]c

needed in equation (4-22). Then, updated values of V were determined using

the FEM as shown in equation (4-22). The advantage of the technique is that

the sparse symmetry of the matrix K is preserved so that a Cholesky decomposi-

tion can be used. The drawback, is that the iterative application of the nonlinear

polarization curves must be repeated for every update of V. Unfortunately, the

simplistic approach of successive substitution failed to converge for all but the

simplest case of a single short well-coated pipe. Also, the approach took a great

deal of time to converge because of the multiple applications of the nonlinear

boundary conditions within the BEM.

4.8 Solution of Combined System by Newton-Raphson Technique

The set of variables that appear in the combined soil domain metal domain

problem are potentials outside the surface of all pipes or tanks (0), unknown

normal electric field on anodes (n. VO), and unknown potential difference (V -

O) for coated and bare protected structures. Introduction of the attenuation code

from equation (4-22) to the problem introduces n more equations and n more

unknowns where n is the number of nodes used to describe the boundary mesh

for the problem. The total number of algebraic equations to be solved is now

2n + k, where k is the number of separate CP systems.

The coupled methods are placed into a global matrix system and one instance

of equation (4-7) is added for self-equilibriation of each CP system

H,,, H,,p 0 0 -47r p Gp,p Ga,p

H,,a Ha, O 0 -47 (D, Gp,a Ga,a

0 0 K, 0 0 V, = 0 0 (4-23)

0 0 0 Ka 0 Va 0 a

0 0 0 0 0 oo Ap Aa

where the subscripts a and p refer to anodes and pipe respectively, double sub-

scripts are from the boundary element method and refer to where the collocation

point and field point lie respectively, and A is the area of each section.

The columns of equation (4-23) are sorted such that all of the unknown vari-

ables are on the left hand side. For anodes, the unknown variable is the current

density (nf V)), and for all other structures it is the potential (0). The reordered

system then looks like

Hp, -Ga,p 0 0 -47r Op Gp,p -Ha,p

H,a G,,, 0 0 -47T VVD, Gp,, -Ha,a

0 0 Kp 0 0 Vp = P 0

0 -fa 0 Ka 0 Va 0 0

0 -Aa 0 0 0 Ioo Ap 0

(4-24)

where the column matrix on the left hand side is the set of unknown variables

and the column matrix on the right side is the set of known boundary conditions.

Equation (4-24) can be rewritten as

[A] V, =[B] (4-25)

Va

where A and B correspond to the matrices in equation (4-24).

If the boundary conditions are constant, the unknown variables can be ob-

tained by a simple matrix inversion. For problems in electrochemistry, the set of

known boundary conditions are given as a set of nonlinear functions of the un-

known variables. Therefore, the solution must be obtained using a technique for

coupled sets of nonlinear algebraic equations. This is done by rewriting equation

(4-25) as a function that is equal to zero

Op 0

n Voa I 0

n. V, = f(V,
[A] V -[B] = 0 (4-26)

Va L J 0

oa = f(v, W. Vo.) 0

(" 0

Since one of the column matrices contains unknown variables, a guess for the

values must be made. Since the guess will not be the correct solution to equation

(4-26), the right hand column vector will not be equal to zero, but a residual

n .VOfl

f( V-Op) [ 1

[A] V -[B] =f( R (4-27)

Va

Rewriting equation (4-26) as a general set of equations of an unknown vector x

F(x) = 0 (4-28)

a Jacobian J, of the set of equations can be written where an element in the Jaco-

bian is given by

i F, = O j (4-29)

Dixj Oxj

Then the Newton-Raphson iteration can be performed as in section 4.5 on

page 62 to obtain the solution vector x such that equation (4-28) holds. Again the

line search is necessary for the the method to work. However, in most non-trivial

cases, even the line searches fail.

In order to get the method to converge to a solution, many techniques were

tried including the minimal residual techniques such as the modified Powell

method,3 Levenberg-Marquort techniques and other methods such as the conju-

gate gradient method. These techniques would converge sometimes if the guess

was nearly the correct value.

4.9 Variable Transformation to Stabilize Convergence

The system of equation that result from (, V, i as the set of variables for the

solution has very poor convergence characteristics. Great effort was required

to bring the system to convergence. Standard line searches and even minimal

residual techniques failed to converge.

A simple variable transformation has been developed that enables good con-

vergence properties using the simple line search method. The variable transfor-

mation needed can be seen easily from the equations used to describe the kinetics

V- Fe ( 1 V-( E 1 V-EH2

i = 10 OFe 1+10 2 10 H2 (4-30)

\1im,02

In equation (4-30), the difference V ( appears many times. A simple transfor-

mation can be written

= V- (

(4-31)

where T represents the driving force for the electrochemical kinetics. The vari-

able V is then chosen as the dependent variable and rewritten as

V = Y + D (4-32)

All of the equations used for the kinetics, bare metal (1-23), coated metal (1-26),

galvanic anodes (1-27), and impressed current anodes (1-30) can be rewritten in

terms of the new variable. For bare metal, one obtains

EFe 1 EO 1 Y EH2

i = 10 Fe + 10 2 ) 10 H2 (4-33)

1lim,02 I

for coated metal there are two simultaneous equations in two unknowns that

describe the current density

i = -( n) (4-34)

Pfilm film

and

Apoe -1 in-Fe in-EO2(Y inEH2)

pore- 10 Pe - b) 10 2 10 H2

A (1 Oblock) lim,02 10

(4-35)

for galvanic anodes

Y = Eanode + 3anode log(i/icorr + 1) (4-36)

and finally for impressed current anodes

S= A Vrectifier + E02 + 102 log(i/icorr + 1) (4-37)

where A Vrectifier is the known setting for the rectifier.

The set of equations for the finite element method must also be rewritten as

K [ + ] = [f]

(4-38)

where f and K do not change. The boundary element method is not changed.

The set of equations that is solved becomes:

0 0 Hpp H,p -47r p Gp,p Ga,p

0 0 H,,, H,,, -45r Y, G,,, G,

0 0 Hp,a Ha,a -47 T, Gp,a Ga,a

W VVOP

K, 0 K, 0 0 O, = f 0 (4-39)

0 K, 0 K, 0 D, 0 F L

O 0 0 0 0 oo Ap Aa

At this point the essential boundary condition for V (T + O) is applied. There is

one degree of freedom in the problem for this potential so one node on one pipe

is selected to have a 0 potential. This is done be replacing the FEM equation for

the single selected node in the system given by equation (4-39) with the equation

Y + = 0 (4-40)

The FEM equation (lower half of the system) is replaced because that is where the

degree of freedom in the reference potential is. If there are multiple CP systems,

then there must be an equal number of essential boundary conditions on T + O.

Once the equations are set up, all unknown variables are placed on the left-hand-

side, which are ( and T for cathodes (pipes, tank bottoms, ship hulls, etc.) and

n Vc and ( for anodes. The reordered system of equations is

0 -G,,p Hp,p Ha,p -47Tr Y Gp,p 0

0 -G,,, H,,a H,,, -47r I V 4, G,,, 0

K, 0 K, 0 0 o, = Fp 0

0 F 0 K, 0 D, -K,

0 -A, 0 0 0 ~o A, 0

(4-41)

The above system of equations is then solved using the same technique as de-

scribed in section 4.5. This system has excellent convergence characteristics and

uses about the same number of iterations as the method without attenuation.

4.10 Conclusions

A method has been developed to model cathodic protection of complicated

pipeline networks. It utilizes the boundary element method to solve Laplac's

equation within the soil and it accounts for the non-zero voltage drop in long

pipes using a finite element solution.

A new variable transformation has been described to stabilize the conver-

gence of the application of the nonlinear boundary conditions.

There are still some issues to be dealt with in the development of a working

model. The creation of an optimized mesh is treated in the next chapter. Verifica-

tion of the solution techniques is done in Chapter 7.

CHAPTER 5

MESH GENERATION

Figure 5-1: Type of mesh needed for model. Three pipes are shown passing

through a point where the soil conductivity changes.

The Boundary Element Method (BEM) and the Finite Element Method (FEM)

rely on a mesh in order to solve the differential equations. A mesh is a set of

lines, triangles, rectangles or other shapes that when pieced together, form the

boundaries of the geometry of interest. The finite element method also must have

the domain enclosed by the boundaries divided up into these shapes. The BEM

needs a mesh on the boundaries only. These shapes are called elements hence the

name of the numerical methods.

Elements use a shape function to approximate the solution along the geom-

etry described by the element. The number of elements should be sufficient to

allow the shape functions to approximate the solution to the governing differen-

tial equation.

One of the most difficult problems to address in computer code is how to dis-

cretize an object such that the solution obtained from the BEM and FEM contains

a desired minimum level of information. For a discretization method to succeed,

it must not only be able to describe the solution accurately, but must also describe

the geometry sufficiently well so that the solution can be obtained.

5.1 Basic Mesh

The basic mesh must account for the current and potential distribution both

axially and around the circumference of the pipe. It is accomplished be divid-

ing the pipe surface (boundary to the soil domain) into 2-dimensional elements

that are piece-wise continuous. The elements are not planar. They are wrapped

onto the pipe geometry like wrapping paper to a gift. Therefore the elements

are 2-dimensional only in a curvi-linear coordinate system. A coordinate trans-

formation must be derived to map the 2-dimensional curvi-linear system to the

3-dimensional cartesian system where the differential equation is solved. In the

curvi-linear system, the element is called the parent element since all the elements

have the same size and shape in that system.

5.1.1 Discretization and Shape Functions

The family of elements that use the same basis (shape) functions to describe

the geometry over the surface of the element as well as the solution to the dif-

ferential equation over the element are called iso-parametric elements.67,74 These

Figure 5-2: Circle discretized with 8 linear triangular elements. The area is 2/2r2

for the 8 triangles which is a 10.% error.

elements have been used in the code that has been developed. The choice of el-

ement shape has been dictated by the geometry being represented. The element

that best fits the sides of a cylinder is a quadrilateral element while triangles are

best used on the ends. The only remaining issue then is to choose how many

nodes and therefore the order of the basis polynomials used to describe the vari-

ation of the geometry and solution over the individual elements.

Basis functions considered for the elements were: 1) constant, 2) linear, 3)

quadratic, and 4) cubic. Of these, the curved functions, quadratic and higher,

make the most sense because of the curvature of the geometry. They also make

sense for representing the solution since they can more easily describe the curve

heading into the singularity in the current density that occurs at junctions be-

tween different boundary conditions where the angle of the geometry at the junc-

tion is not 90 degrees. Table 5.1 shows the capability of 3 types of elements to rep-

Table 5.1: Error of Constant, Linear, and Quadratic elements in representing a

circle whose circumference is 7r and has unit diameter.

Nodes Constant %err Linear %err Quadratic %err

6 3.464 10.3% 3.0 4.51% 3.094 1.51%

8 3.314 5.48% 3.062 2.55% 3.125 0.53%

10 3.249 3.43% 3.090 1.64% 3.134 0.23%

12 3.215 2.35% 3.106 1.14% 3.138 0.12%

16 3.183 1.31% 3.121 0.64% 3.140 0.038%

24 3.160 0.58% 3.133 0.29% 3.1414 0.0077%

resent a circle as a boundary to an area.32 At 12 degrees of freedom, the quadratic

element is better by an order of magnitude in the error. In fact, 32-36 degrees of

freedom are needed with linear or constant elements to achieve the same accu-

racy representing a circle as do 12 degrees of freedom using quadratic elements.

Figure 5-3 shows eight linear triangles used to represent a circle. The error in

the surface area is 10.%. That error may not be the most important however. In

the end, it is how well the mesh can approximate the solution to the differential

equation that is most important.

5.1.2 Continuity of the Shape Functions

It is known from examination of Laplace's equation that the potential distri-

bution along any boundary is at least piece-wise continuous. Therefore, the basis

functions must also exhibit this type of continuity. These elements are commonly

referred to as CO continuous elements. Since it is not required for the elements to

be square integrable (Ho)t for the boundary element method or that their deriva-

tives be continuous (C1) or square integrable (H1), just about any piece-wise con-

tThese terms (Co, Ho, etc.) refer to Sobolev spaces.74