Dissertation
Dissertation
By
DOUGLAS P. RIEMER
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.
ii
TABLE OF CONTENTS
Page
ACKNOWLEDGEMENTS . . . . . . . . . . . . . . . . . . . . . . . . . . . ii
LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x
ABSTRACT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxi
CHAPTERS
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
iii
3 SOLUTION OF GOVERNING PDE . . . . . . . . . . . . . . . . . . . . 29
3.1 Boundary Element Method Development . . . . . . . . . . . . . . . 30
3.2 Symmetric Galerkin Boundary Element Method . . . . . . . . . . . 32
3.3 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.3.1 Bare Steel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.3.2 Coated Steel . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.3.3 Anodes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.4 Infinite Domains . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.5 Half Spaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.6 Layers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.7 Development of Finite Element Method for Internal Domain . . . . 41
3.7.1 Pipe Shell Elements . . . . . . . . . . . . . . . . . . . . . . . . 43
3.7.2 Applying Elements to the FEM . . . . . . . . . . . . . . . . . 44
3.8 Bonds and Resistors . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.9 Coupling BEM to FEM . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.10 Current Flow in the Pipe . . . . . . . . . . . . . . . . . . . . . . . . . 48
3.11 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
5 MESH GENERATION . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
5.1 Basic Mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
5.1.1 Discretization and Shape Functions . . . . . . . . . . . . . . 76
5.1.2 Continuity of the Shape Functions . . . . . . . . . . . . . . . 78
5.1.3 Quadratic Iso-Parametric Rectangles . . . . . . . . . . . . . . 81
5.1.4 Quadratic Triangles . . . . . . . . . . . . . . . . . . . . . . . . 83
5.2 Long Pipes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
5.3 Holidays . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
5.3.1 Round Holidays . . . . . . . . . . . . . . . . . . . . . . . . . 86
iv
5.3.2 Rectangular Holidays . . . . . . . . . . . . . . . . . . . . . . 91
5.4 Pipe Crossings and Shadowing . . . . . . . . . . . . . . . . . . . . . 92
5.5 Soil Type Divisions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
5.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
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
v
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
APPENDICES
vi
B OpenGL CODE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 205
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
REFERENCES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 252
vii
LIST OF TABLES
Table page
9.2 Normalized benchmark result data for many types of computers. . 131
9.5 Conditions used for all hard drive performance tests . . . . . . . . . 135
11.3 Results of two coupon model. The small coupon used a polariza-
tion curve with ilim,O2 = 1.0 µA/cm2 (curve 1). The large coupon
used ilim,O2 that was an order of magnitude smaller (curve 2). . . . . 160
viii
11.4 Results of two coupon model. The large coupon used a polariza-
tion curve with ilim,O2 = 1.0 µA/cm2 (curve 1). The small coupon
used a ilim,O2 that was an order of magnitude lower(curve 2). . . . . 160
ix
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-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 dΓi and dΓ j . . . . . . . . . . . . . . . . . . . . . . . . . . 34
x
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
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
√ 2 discretized with 8 linear triangular elements. The area is
2 2r for the 8 triangles which is a 10.% error. . . . . . . . . . . . . 77
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
xi
5-9 Current density distribution on a circular holiday. Parameter J is a
function of the average current density iavg , the inverse Tafel slope
αc F
RT
, the holiday radius ro , and the soil resistivity κ. . . . . . . . . . . 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-15 Mesh at a crossing of two pipes. Both pipes have the highest den-
sity of elements at the crossing point. . . . . . . . . . . . . . . . . . . 93
5-17 Completed mesh with pipes passing through soil type division. . . 94
xii
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
8-1 A screen shot and movie (click on image above) that shows some of
the interactive graphics techniques needed for three dimensional
geometry. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 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 Microsoftr 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
xiii
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-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
10-1 Pipe that exhibits many coating holidays. Pipe was coated in the
late 60’s and dug up in 1996. . . . . . . . . . . . . . . . . . . . . . . . 137
xiv
10-6 Distant view of a under-protected holiday. Lower right portion of
pipe with holiday is also exhibiting stray current. . . . . . . . . . . 143
xv
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 ξ = 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 Ω 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-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-6 Radial potential and current density distribution when the oxygen
concentration
p 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
xvi
12-7 Potential distribution with distributed anodes and uniform oxygen
concentration profile. . . . . . . . . . . . . . . . . . . . . . . . . . . . 180
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
xviii
KEY TO SYMBOLS
i = Current density
βi = Tafel slope
F = Faraday’s constant
R = Gas constant
T = Temperature
k = Kinetic parameter
κ = Conductivity
ci = Concentration of species i
t = Time
zi = Charge on species i
ui = Mobility of species i
I = Ionic strength
xix
qi = Kirchhoff transformation variable for species i
Q = Quasi-potential
Γ = Boundary of a domain
G = Green’s function
Ψ = defined as V − Φ
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
By
Douglas P. Riemer
December 2000
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-
A butane line explosion, which cost the lives of two persons, was caused by exter-
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.
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
xxi
mathematical model that can describe correctly the complex interactions seen in
soil, steel, and connections between pipes and anodes. The model accounts for
gen evolution reactions. It also accounts for the nonuniform current and potential
Current and potential distributions are calculated around the circumference and
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
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
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
rosion mitigation strategies has in general been successful. Most pipeline failures
loss of $75 million and 4 lives.∗ Additional failures, later deemed to have been
∗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.
1
2
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
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
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
3
(a)
(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
4
pipeline safety. The joint bill5 and amendments6 propose new research that in-
pipelines.
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,
constructed during the construction boom of the 1960s and 1970s are aging, and
property damage and loss of life. The design of cathodic protection systems is
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.
a more desirable one on another surface. In the case of buried metallic structures,
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.
6
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
High resistance coatings are used on most pipelines to reduce the amount
and scrapes which expose bare steel cannot be avoided during installation. Its
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
the boundary conditions are constant on any given side. However, the boundary
conditions for pipelines under cathodic protection are not uniform or constant on
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
n · ∇Φ = c at r = R) since
use the uniform boundary condition (Φ = c at z = 0 or ~
7
the potential and current density will be different for a discrete holiday than for
that relate the net current density of the chemical reactions on the pipe surface to
Some assumptions about the knowledge of the reader are made in the fol-
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
trode. They proposed that the electrochemical reactions that take place on the
structures under cathodic protection in soil form alkaline environments.3 The cor-
R.D.S.
Fe + 2OH− −−−→ Fe(OH)2 + 2e− (1-2)
and
Fe(OH)2
Fe2+ + 2OH− (1-3)
8
where R. D. S. stands for Rate-Determining Step. The sum of reactions (1-2) and
Fe
Fe2+ + 2e− (1-4)
which is an anodic reaction when the reaction proceeds from left to right.
Fe(OH)2
Fe2+ + 2OH− (1-5)
and
R.D.S.
Fe(OH)2 + 2e− −−−→ Fe + 2OH− (1-6)
VF
ia = ka,Fe aOH− e[ RT ] (1-7)
where ka,FE is the kinetic rate constant, aOH− is the activity of OH− ions, V is de-
the Gas Constant and T is temperature. The activity of the OH− ions is linear
−VF
ic = kc,Fe aFe2+ aOH− e[ RT ] (1-8)
where aFe2+ is the activity of Fe2+ ions, and the net current density i is given by
i = ia + ic (1-9)
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
RT κc,Fe aFe2+
Vo = ln (1-10)
F κa,Fe
where Vo is defined as the equilibrium potential. Then the exchange current den-
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
VF −VF
i = io e[ RT ] − e[ RT ] (1-13)
or
VF
i = io 2 sinh (1-14)
RT
The number of unknowns in equation (1-13) can be reduced by changing the
2.303F
βa = (1-15)
RT
2.303F
βc = (1-16)
RT
and then moving all pre-exponential terms into the exponent to get
V − Ea
i = 10[ βa ] − 10[ V−βcEc ] (1-17)
10
For freely corroding systems, conservation of charge requires that the anodic
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
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, 18, 19
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-
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
The total current is given by the sum of equations of the form of (1-17). There
Equation (1-22) can often be simplified by observing that some terms can be ig-
from this equation by setting it equal to zero and solving for V = Vcorr . This is the
potential is taken with a reference electrode just outside the surface of the metal.
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
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
where ilim,O2 is the mass transfer limited current density for oxygen reduction.
This means that the portion of the current due to oxygen reduction cannot exceed
since it is assumed that there is always sufficient water next to the pipe for that
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-
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
1.6 Coatings
ings are applied to the steel surface before burial. These coatings are usually
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
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.
14
Φinner Φinner
Riemer and Orazem have modified equation (1-23) to take into account the
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
Φ − Φin
i= (1-24)
ρδ
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, ρ is the resistivity of the coating
and δ 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
15
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
V −Φin − EO
!−1 −(V −Φin − EH )
A (Φ − Φin ) V −Φin − EFe
1 βO
2
βH
2
= 10 βFe
− − 10 2 − 10 2
Apore ρfilm δfilm (1 − αblock ) ilim,O2
(1-26)
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
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 − αblock ) 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
100
10-2
10-3
10-4
10-5
10-6
10-10
-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.
where the half of interest is the corrosion (oxidation) of the anode material. Both
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
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-
because the driving force for the reduction reaction is higher than that for steel
The three parameters of the model are generally known for all types of gal-
Impressed current anodes do not obtain their driving force from the activity
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
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
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
and
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
i = iO2 10(V−Φ−4Vrectifier −EO2 )/βO2 − 1 (1-30)
where 4Vrectifier 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 βO2 is the tafel slope for the oxygen
evolution reaction.
19
Table 1.2: Parameters for the oxygen and chlorine evolution reactions.25
reaction Equil Potential (E) Tafel Slope (β )
O2 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
Parameters for reactions (1-28) and (1-29) are given in Table 1.2 except for iO2
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
model has been developed under the constraints of the kinetics are transport as-
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
Furthermore, the appendices assume the reader is well instructed in the C and
There are several good sources for C and especially C++ programing informa-
20
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.
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
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
Results of the model and their implications for design of cathodic protection
The intervening chapters discuss how the model was verified, and how per-
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.
tions for the soil domain which is governed by Laplace’s equation. The equation
is given by
i = −κ∇Φ (2-1)
where i is the current density vector, κ is the conductivity of the domain and Φ
21
22
bulk yields
∇·i =0 (2-2)
∇2 Φ = 0 (2-3)
a dilute electrolyte, the flux for any species can be written as the sum of three
where v is the fluid velocity, Di is the diffusion coefficient for species i, zi is the
equation
Di = RTui (2-6)
The current density i, that results from equations (2-4) and (2-5) can be expressed
as
i = F ∑ zi Ni (2-7)
i
neutrality
∑ zi ci = 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 ∑ zi ci =0 (2-9)
i
The reaction terms also sum to zero because of charge conservation in a homoge-
neous reaction
F ∑ zi Ri = 0 (2-10)
i
κ = F2 ∑ z2i ui ci (2-13)
i
equation
∇2 Φ = 0 (2-14)
24
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
1
2∑
I= z2i ci (2-15)
i
κ = F2 ∑ z2i ui ci (2-16)
i
F2
RT ∑
κ= z2i Di ci (2-17)
i
The diffusivity of a typical species is10−5 (cm2 /s). The conductivity is then written
as
1 F2 D
κ= =
ρ RT i ∑ z2i ci (2-18)
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)
2F2 Dρ1.5
Using typical values of T = 290K, and D, and ρ taking values as above, a value
The resistivity of pipe steel is 9.6 x 10−6 Ωcm. For an 18 inch diameter pipe,
the resistance per linear foot of pipe is about 2.3 x 10−6 Ωft . If 10 Amps must travel
5 miles on this pipe, the potential drop in the steel must be 600 mV which is a
drop in the pipe steel made by Esteban et al. cannot be made for long pipes.33 The
The flow of current through the pipe steel, anodes and connecting wires is
∇ · (κ∇V) = 0 (2-21)
where V is the departure of the potential of the metal from an uniform value
For instance, the wire connecting the pipe the anode will be made of a different
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
which is generally valid, does not apply within a diffusion layer (about 50 µm
in thickness41 ) near the pipe and anode where electrochemical reactions generate
tion can be neglected, the system is at steady-state, and that concentration and
Since ci is a function of Φ,
dci
∇c i = ∇Φ (2-24)
dΦ
thus
dci
Ji = − Di − zi ui Fci ∇Φ (2-25)
dΦ
Under the assumption that the diffusion coefficient is constant,
Ji = f i (Φ)∇Φ (2-26)
where
dci zi F
f i (Φ) = − Di − Di − ci (2-27)
dΦ RT
27
A new variable, obtained through the Kirchhoff transformation,43 has the prop-
erties
and
Z Φ
qi = − f i (Φ)dΦ (2-29)
0
Ji = −∇qi (2-30)
The current density is given by the sum of contributions from the flux of each of
i = F ∑ zi Ji (2-31)
i
i = − F ∑ z i ∇q i (2-32)
i
Q = F ∑ zi qi (2-33)
i
such that
i = −∇ Q = F ∑ zi Ji (2-34)
i
Laplace’s equation
∇2 Q = 0 (2-35)
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
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
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.
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
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-
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 first task is to find a method to solve the differential equation given ar-
the boundaries to a form for which a solution exists. Moulton first described the
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
29
30
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
tion.54 Aoki and also Telles provided the first practical demonstrations of utilizing
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
soil conductivity changes for the case of a single pipe with no angular variations
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
31
or weighted residual of the governing differential equation. If the PDE takes the
form
∇2 u − f = 0 (3-1)
where f is a forcing term that is a function of position only, then one would write
where Ω is the domain and w is any weighting function. In the case of Laplace’s
This equation holds for all weighting functions w. Equation (3-3) is integrated
with Γ 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
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
∇2 Gi, j + δi = 0 (3-6)
32
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
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 Φi to the
boundary. A Cauchy principle value is introduced in the integral on the left side
non-symmetric matrices. Even though the diagonals of the matrices have the
33
largest magnitude, they are not diagonally dominant. This means that the fast
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
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,
φ, of the source element, i, and integrate the resulting equations over the entire
boundary
Z Z Z Z Z
φi Φi dΓi + Φ ~ n · ∇Φ) dΓ j dΓi
φi n · ∇Gi, j dΓ j dΓi = φi Gi, j (~ (3-11)
Γ Γ Γ Γ Γ
and
Z Z Z
φi κ (~ni · ∇Φi ) dΓi + Φκ ~
φi ni · ∇ ~
n j · ∇Gi, j dΓ j dΓi =
Γ Γ Z Γ Z
~ni · ∇Gi, j κ ~n j · ∇Φ dΓ j dΓi
φi (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
∫dΓi
∫dΓj
Figure 3-1: Symmetric integration. Lines between the two circles represent the
Green’s function relation between the two integration points on the surfaces dΓi
and dΓ j .
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
can only be used when the boundary conditions are constant. The symmetry falls
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
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
of the buried structure where bare steel is exposed to the soil (coating holiday or
uncoated pipe).
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
3.3.3 Anodes
The boundary conditions for the anodes must be an essential type 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
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-
Equations (3-13) and (3-14) are then used to form a well-posed problem definition
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
a second boundary, Γ̄, 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
Z Z
c i Φi + Φ ~
n · ∇Gi, j dΓ + Φ ~
n · ∇Gi, j dΓ̄ =
Γ Z Γ̄ Z
n · ∇Φ) dΓ +
Gi, j (~ n · ∇Φ) dΓ̄
Gi, j (~ (3-15)
Γ Γ̄
If the radius of the new boundary is taken to infinity the limit of the integrals
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
π 2π Φ 2
Z Z
r sin φ dθ dφ = 4π Φ
lim − (3-18)
r→∞ 0 0 r2
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
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
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 x0 , one may
write66
which implies that the two source intensities σ (x), and σ (x0 ) are equal and have
the same sign. The sign is the same because the outward normal vectors have
38
1 1
Gi, j = + (3-20)
4π r(xi , x j ) 4π r(xi , x0j )
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
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
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
Re f lections+1
1
Gi, j = ∑ 4π r(xi , x j,k )
(3-21)
k
where the index k > 1 refers to the reflection about some plane k. For k = 1, x j,k is
the field point on the real object. If the soil surface is the only reflection used, the
Reflection
Air
r'
r
Source
Soil
r''
Rock Reflection
Figure 3-2: Reflections of Green’s function to account for boundaries with zero
normal current.
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
is seen in Figure 3-4. For these calculations, the anode was placed far from the
40
Figure 3-3: Rectangular enclosure created with a primary Green’s function and 6
reflections.
Total Current (Amps)
0.21
0.19
0.17
0.15
0 50 100 150 200
Depth of zero conductivity (ft)
Figure 3-4: Total cathodic current driven to a pipe with an underlying non-
conducting region.
41
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
a
I = Ir→∞ − √ (3-22)
r
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
The data used in Figure 3-4 was then then fit to the equation
a
i = √ + ir=∞ (3-23)
r
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.
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
42
0.25
0.15
0
0 0.2 0.4 0.6 0.8 1
1/sqrt(depth)
−0√.064
Figure 3-5: Fit of the data from Figure 3-4 to i = r
+ 0.2118.
equation (2-21)
Z
w∇ · ~
~κ · ∇VdΩ = 0 (3-24)
Ω
where w is any weighting function. For pipe steel, κ is just a constant times an
identity matrix
1 0 0
~~κ = κ 0 1 0
(3-25)
0 0 1
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-
where ~
n is the outward normal vector from the boundary of the domain. The do-
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-
n
Z Z Z ZZ
∑ (e)
κ∇φ(e)
i · ∇φ(e)
j dx dy dz V j = − κφ(e)
i (~n · ∇V)(e)
i ds (3-27)
j=1
where κ is the scalar component of the property tensor, φi is a shape function and
the sum goes from 1 to the total of all the shape functions of all the 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
where ξ and η define the outside surface of the pipe and ζ is the outward normal.
cause the scale of the problem is many orders of magnitude greater in the ξ 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.
44
These functions are obtained through the product of two Lagrange interpolating
polynomials of the same order, one in ξ and one in η . The result is a family of el-
ements with square parent elements. They include the four-node linear element,
mation
dx dy dz = |J| dξ dη dζ (3-28)
where xk is one of the cartesian coordinates, |J| is the determinate of the Jacobian,
and Ω(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
45
shell elements used for pipes, the integral is only performed numerically over
the surface dξ dη 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
" ZZ #
∂φi ∂φ j
q ZZ
∑ h ∑
Ω(e) k ∂ xk
κ·
∂ xk
|J|dξ dη = − wκ (~n · ∇V) ds
Γ(e)
(3-31)
i
with the thickness of the steel, h, a parameter. The normal vector has the same
rule, but starting from the derivatives in cartesian coordinates, one writes
∂φ ∂ x ∂φ ∂ y ∂φ ∂ z ∂φ
+ + ∂ z ∂ξ
∂ξ
∂ x ∂ξ ∂ y ∂ξ
∂φ ∂ x ∂φ ∂ y ∂φ ∂ z = ∂φ (3-32)
∂ x ∂η + ∂ y ∂η + ∂ z ∂η ∂η
∂φ ∂ x + ∂φ ∂ y + ∂φ ∂ z
∂φ
∂ x ∂ζ ∂ y ∂ζ ∂ z ∂ζ ∂ζ
If J− 1
i, j is the element from the ith row and the jth column of the inverse of the
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 ξ and η directions, a nine-
point Gauss rule (3x3) applied in both directions will give an exact (to machine
precision) result.
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
47
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
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
1 −1
R R
Ke = (3-36)
−1 1
R R
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
48
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
Mho/m. Using the variables for potential in the pipe, V, and potential in the soil,
κsoil
n · ∇V) =
(~ n · ∇Φ)
(~ (3-38)
κsteel
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
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 determine the direction of current flowing through a bond. The method for
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
9
∂φk
iξ (ξ, η ) = −κ ∑ Vk (3-40)
k=1
∂ξ
and
9
∂φk
iη (ξ, η ) = −κ ∑ Vk (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
∂ VI
ix, I = −κ = (J −1 )11 iξ (ξ I , η I ) + (J −1 )12 iη (ξ I , η I ) (3-42)
∂x
∂ VI
i y, I = −κ = (J −1 )21 iξ (ξ I , η I ) + (J −1 )22 iη (ξ I , η I ) (3-43)
∂x
∂ VI
iz, I = −κ = (J −1 )31 iξ (ξ I , η I ) + (J −1 )32 iη (ξ I , η I ) (3-44)
∂x
for the z component, where (ξ I , η I ) are the optimal current sampling points and
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 (ξ, η ) = c1 + c2 ξ + c3 η (3-45)
where the constants are found through a least squares regression. Therefore, a
where ~ioptimal 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
IV
ex = ∑ [c1 + c2 ξk + c3 ηk − 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-
∂ ex
=0 (3-48)
∂cj
which yields
IV
∑ [c1 + c2 ξk + c3 ηk − ik ] = 0 (3-49)
k= I
for j = 1
IV
∑ [c1 + c2 ξk + c3 ηk − ik ]ξ = 0 (3-50)
k= I
for j = 2, and
IV
∑ [c1 + c2 ξk + c3 ηk − ik ]η = 0 (3-51)
k= I
51
or in matrix notation as
and [Q] is
1 1 1 1
[Q] = 1
−√3
1
√ √1 √1 (3-56)
− 3 3 3
1
√ √1 1
√ √1
− 3 3 − 3 3
The 9 nodal values of the flux are then found from equation (3-45) and written
in matrix form as
i1 1 −1 −1
i2 1 0 −1
i
3 1 1 −1
1 −1 0
i4
c1
i = (3-58)
1 0 0 c2
5
i6
1 1 0 c3
i
7 1 −1 1
i8 1 0 1
i9 1 1 1
Upon substituting the values of c j from equation (3-54), the transformation matrix
becomes
√ √
1+2 3
1 1 1−2 3
√√ √ √
1− 3 1+ 3
1+ 3 1− 3
√ √
1 1−2 3 1+2 3 1
√ √ √ √
1+ 3 1+ 3 1− 3 1− 3
1
[TR] = 1 1 1 1 (3-59)
4
√ √ √ √
1− 3 1− 3 1+ 3 1+ 3
√ √
1 1+2 3 1−2 3 1
√ √ √ √
1− 3 1+ 3 1− 3 1+ 3
√ √
1−2 3 1 1 1+2 3
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.
and
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
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.
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
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.
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
55
56
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-
mate functions on each individual element. Example functions may be 0th 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-
summed together form the complete integral over the entire boundary
Z Z
ci φi + ∑ n · ∇Gi, j dΓ j = ∑ Gi, j (~
Φ ~ n · ∇Φ) dΓ j
(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.
value, the notation used here will denote the matrix resulting from the left hand
n · ∇Φ)
HΦ = G (~ (4-2)
boundary.
57
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
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-
Then the gradients of potential are known to be zero everywhere which leads to
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.
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
58
may be obtained at any point within the domain by use of equation (3-8). In this
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
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
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-
singularity never appears if the interior point chosen is sufficiently far from the
boundary.
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
59
small compared to the size of a pipe or holiday. Therefore, it can be assumed that
mately the spot where V is evaluated and the same nodal locations can be used
in both methods.
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 − Φ,
4.2 Self-Equilibration
ment formulation as developed in Chapter 3 will not show this feature without
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
system56, 69
Z
∑ Γi
κ~nˆ · ∇ΦdΓi = 0 (4-7)
i
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
values that are placed in this column come from equation (3-17) and are all 4π or
1 depending on where the 4π 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
that differ by a constant. Therefore, at least one element in the system must have
unique solution.
Riemer and Orazem describe in detail the development needed to model in-
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
and the column matrix q defined in the usual way. When there are two or more
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
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.
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
62
current distributions.24 This meant a great many triangles were needed where
sharp profile of the solution in these regions. Therefore, only small sections of
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
special methods to enforce them in the boundary element method. Aoki et al.
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
tained. One starts by writing a residual after the columns corresponding to the
n · ∇Φ. Therefore equation (4-2) does not hold. For exposed steel in a soil
tion ~
equation (1-23). The Taylor series expansion truncated at the linear term gives
Φ(k)
i ∆Φ(k)
i
0= H
··· +
······
n · ∇Φ)(k)
(~ i ∆(~
n · ∇Φ)(k)
i
Φ(k)
i
0
−G
· · · · · · · · ·
+
···············
(4-12)
∂ f i Φ(k)
− f i Φ(k) i − ∂Φ
i
∆(~
n · ∇Φ)(k)
i
The residual error, RResidual , is obtained by subtracting equation (4-11) from equa-
tion (4-12)
..
0 . 0 ∆Φ(k)
i
RResidual = − H + G · · · · · · · · ·· · · · · · · · · ········· (4-13)
.. ∂ f i Φ(k)
∆ (~
n · ∇Φ)(k)
i
0 . − ∂Φ i
64
h i
The matrix containing −∂ f i Φ(k)
i /∂ Φ is an m × n diagonal matrix which is
matrix, the resultant that is added to the H matrix is square. The new values of Φ
Φ(k
i
+1)
= Φ(k)
i + ∆Φi
(k)
n · ∇Φ)(k
(~ i
+1)
= (~n · ∇Φ)(k)
i + ∆(~
n · ∇Φ)(k)
i (4-14)
The new values can be substituted back into equation (4-11) to obtain a new resid-
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
A subroutine takes the column matrix containing the Newton-Raphson step and
n · ∇Φ) and λ is a
where x is the column matrix of unknown variables (Φ and ~
4-1. The initial calculation is done with the full Newton step (λ = 1). If the sum
of squares or the residual error, g, obtained from equation (4-15) is reduced, the
selected using the same Newton step and a new residual and objective function
65
Calulate
guess x Objective
Function
Calulate
Newton Step
Get λ
Calulate
Objective
Function
No
Is λ too No
Objective
small? smaller?
Yes
Yes
Cannot Stopping
No
converge criteria met?
Yes
Done
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
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 ZZ
fi = ∑ − Γ w`,i κsoil (~n · ∇Φ) ds (4-17)
`=1 (e)
where f i 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
n · ∇Φ is either a known
surfaces of all structures in the model. Since the value of ~
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
nce ZZ nen
fi = ∑ − Γ w`,i κsoil ∑ φ`,k (s) (~n · ∇Φ)`,k ds (4-18)
`=1 (e) k=1
of the normal electric field within the element n, and φ(s)k is the shape function
whose value is one at node k. Using the property of integrals that an integral of a
67
nce nen ZZ
fi = ∑ ∑ − Γ w`,i κsoil φ`,k (s) (~n · ∇Φ)`,k ds (4-19)
`=1 k=1 (e)
or
n · ∇Φ]
f i = F̂i [~ (4-21)
These sub-matrices are then assembled with the right hand column matrix from
n · ∇Φ]
K [V] = F̂ [~ (4-22)
where F̂ is a matrix formed from the assembly of equation (4-21) such that the
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
n · ∇Φ]
calculated using the BEM. The result was the values of the load vector, [~
the FEM as shown in equation (4-22). The advantage of the technique is that
68
tion can be used. The drawback, is that the iterative application of the nonlinear
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
The set of variables that appear in the combined soil domain metal domain
problem are potentials outside the surface of all pipes or tanks (Φ), unknown
Φ) 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
The coupled methods are placed into a global matrix system and one instance
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
n · ∇Φ), and for all other structures it is the potential (Φ). The reordered
density (~
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.
If the boundary conditions are constant, the unknown variables can be ob-
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
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
Φp
n · ∇Φa
~
f (V, Φ p )
[A] Vp − [B] = R (4-27)
n · ∇Φa )
f (V, ~
Va
Φ∞
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
∂ Fi
Ji, j = (4-29)
∂ x j x`6= j
71
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
tried including the minimal residual techniques such as the modified Powell
gate gradient method. These techniques would converge sometimes if the guess
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
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 −Φ− EO
!−1 V −Φ− EH
V −Φ− EFe
1 βO
2 − βH
2
i = 10 βFe
− + 10 2 10 2 (4-30)
ilim,O2
In equation (4-30), the difference V − Φ appears many times. A simple transfor-
Ψ= V−Φ (4-31)
72
where Ψ represents the driving force for the electrochemical kinetics. The vari-
V =Ψ+Φ (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
for coated metal there are two simultaneous equations in two unknowns that
− (Φ − Φin )
i= (4-34)
ρfilm δfilm
and
!−1
Apore Ψ+Φ−Φin − EFe 1
Ψ+Φ−Φin − EO
2 (
− Ψ+Φ−Φin − EH
2 )
βO βH
i= 10 β Fe
− − 10 2 − 10 2
A (1 − αblock ) ilim,O2
(4-35)
The set of equations for the finite element method must also be rewritten as
K [Ψ + Φ] = f
(4-38)
73
where f and K do not change. The boundary element method is not changed.
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
Ψ+Φ=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,
Once the equations are set up, all unknown variables are placed on the left-hand-
side, which are Φ and Ψ for cathodes (pipes, tank bottoms, ship hulls, etc.) and
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
equation within the soil and it accounts for the non-zero voltage drop in long
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-
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)
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
Elements use a shape function to approximate the solution along the geom-
75
76
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
it must not only be able to describe the solution accurately, but must also describe
The basic mesh must account for the current and potential distribution both
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
curvi-linear system, the element is called the parent element since all the elements
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
77
√
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-
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-
78
freedom are needed with linear or constant elements to achieve the same accu-
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
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
be square integrable (H 0 )† for the boundary element method or that their deriva-
tives be continuous (C1 ) or square integrable (H 1 ), just about any piece-wise con-
1 1
0.8
0.6
u η , y0
u η , y1 0.4
u η , y1
u η , y0
0.2
1 x2
0
1 x2
0.2
1 x2
0.4
1 x2
0.6
0.8
1
1 1 0.8 0.6 0.4 0.2 0 0.2 0.4 0.6 0.8 1
1 u η , x0 , u η , x1 , u η , x1 , u η , x0 , x , x , x , x 1
Figure 5-3: Four quadratic elements with eight degrees of freedom representing
a circle.
2
80
0.5
50
40
30
0 20
10
0
0
10
20
tinuous function could be used as the basis functions. However, the finite ele-
ment method used within the steel of the pipe must use square integrable (H 0 )
shape functions.74 It is also necessary to have the solutions between the soil and
steel match up. Therefore the same element needs to be used in both numerical
methods. Lagrange polynomial basis functions are very well developed in the
literature,67, 74 and are often needed for boundary element methods.63 They work
well for FEM, so they are used when rectangular shapes are needed such as the
walls of the pipe. In the case of this code, quadratic Lagrangian iso-parametric
triangles are used for this purpose. They exhibit C0 continuity with the Lagragian
type elements. The triangles can then be used to increase or decrease the local
81
Figure 5-5: Use of triangles to increase the degrees of freedom in the center region
while maintaining C0 continuity.67
can also be used to describe difficult shapes such as coating holidays, ends of
pipes, and tank bottoms. Figure 5-5 shows how triangular elements can be used
to increase the density of elements in a small region. Small regions occur any
of the mesh that constitutes the pipes and anodes. The parent element contains 9
nodes and nine shapes functions. The geometry and solution is written in terms
9
u(η, ξ ) = ∑ φi (η, ξ )ui (5-1)
i=1
where φi is the shape function and η and ξ are homogeneous coordinates in the
parent element that take values from -1 to 1, u is the quantity of interest such as
82
7 8 9
4 5 6
η
1 2 3
ξ
Figure 5-6: Quadratic iso-parametric rectangular parent element. Numbering of
nodes and shape functions given in figure.
the x coordinate or the solution to the differential equation, and ui is the value of
the quantity of interest at node i. The nodes are numbered as in Figure 5-6. The
1
φ1 = η (η − 1)ξ (ξ − 1) (5-2)
4
1
φ2 = (1 − η 2 )ξ (ξ − 1) (5-3)
2
1
φ3 = η (1 + η )ξ (ξ − 1) (5-4)
4
1
φ4 = η (η − 1)(1 − ξ 2 ) (5-5)
2
φ5 = (1 − η 2 )(1 − ξ 2 ) (5-6)
83
η
6
1 4 2
ξ
Figure 5-7: Quadratic iso-parametric triangular parent element. Numbering of
nodes and shape functions given in figure.
1
φ6 = η (1 + η )(1 − ξ 2 ) (5-7)
2
1
φ7 = η (η − 1)ξ (1 + ξ ) (5-8)
4
1
φ8 = (1 − η 2 )ξ (1 + ξ ) (5-9)
2
and
1
φ9 = η (1 + η )ξ (1 + ξ ) (5-10)
4
Quadratic triangles used in this code have 6 nodes per triangle and are C0
continuous with the rectangular elements also used. Figure 5-7 shows the num-
bering of the nodes in the parent element. The shape functions for each node are:
φ1 = [1 − (ξ + η )] [1 − 2(ξ + η )] (5-11)
84
φ2 = ξ (2ξ − 1) (5-12)
φ3 = η (2η − 1) (5-13)
φ4 = 4ξ [1 − (ξ + η )] (5-14)
φ5 = 4ξη (5-15)
and
φ6 = 4η [1 − (ξ + η )] (5-16)
Long pipes must be meshed with a minimum number of nodes, yet account
for the nature of the solution. If elements of an aspect ratio of one are used to
create the mesh for the entire pipe, only very short sections of pipe could be
fore, the aspect ratio of any element cannot be restricted to something near one.
In fact, the aspect ratios need to be as large as 5000 to 1 (length to width) to model
the complete potential and current distributions for long pipes with the limited
resources available on a PC (October, 2000 - 1GB RAM). However, there are se-
vere numerical difficulties encountered with elements of these aspect ratios. The
High aspect-ratio elements cannot be used at the ends of the pipe or where
there is a change in the nature of the boundary condition i.e., a change in the
85
parameters for the boundary condition. It is also not wise to have extremely
to gradually increase the aspect ratio until a maximum value is reached. At that
point, elements of the same size are used until a point that a change in boundary
conditions or the end of the pipe is approached. Then, the aspect ratio must
be gradually decreased until a low aspect ratio element is placed at the point
of change. The trick is have the meshing algorithm place the low aspect ratio
element at precisely the point of the change. The algorithm creates the mesh
from both end simultaneously, meeting in the middle with a high aspect-ratio
element. The effect is to make the mesh between any two boundary condition
changes symmetric about the center of the segment of pipe being meshed. Figure
5.3 Holidays
As stated in chapter 1, coating flaws that expose steel to the soil environment
(holidays) are a serious problem. Therefore, the model must be able to account
for the higher rates of reaction that take place on the exposed section of steel. This
86
is done through the kinetic expressions used as the boundary conditions. Since
any individual element must have the same boundary condition expression, the
mesh must provide entire elements that define the size and location of a holiday.
Also, when using second order polynomials to describe the solution over any
individual element, more than one element is needed on the holiday to account
properly for the nonuniform current distribution. The sample calculations pre-
sented in Figure 5-9 illustrate the need for refined meshing of coating holidays.
For even small values of J = −αc Fro iavg / RT κ, a single second order polynomial
cannot adequately represent the solution. Therefore, the diameter of the holiday
needs at least two quadratic elements for small values of J and four or more for
larger values.
Holidays are meshed is two possible ways. If the holiday is round, it is easiest
to use a constrained Dalauny mesh generator.75, 76 If the holiday is not round such
For the this model, the method that was used to generate a Dalauny mesh
comes from Jonathan Richard Shewchuck.77 The code is publically available with
line graph (PSLG). It then creates a mesh of linear or quadratic triangles that
author of the code provided a simple interface to the generator that utilized the
87
1.8
1.6
1.4
1.2
1
i/i avg
0.8
0.6
J=1.92 J=-αcFroiavg/RTκ
0.4 J=3.30
J=4.88
0.2
J=0.35
0
0 0.2 0.4 0.6 0.8 1
r/ro
Figure 5-9: Current density distribution on a circular holiday. Parameter J is a
function of the average current density iavg , the inverse Tafel slope αRT
cF
, the holiday
radius ro , and the soil resistivity κ.
88
Figure 5-10: A mesh that includes two round holidays. This is the output from
a standard mesh generator.77 Further processing is needed to map the mesh to a
pipe and adjust the holidays so they are round.
same command line options as the stand-alone code. The PSLG input is created
as the boundary to the region being divided into elements. This is usually a
rectangle, a circle, or a half circle which are trivial to create and a list of holes
(holidays) that must be included. The output is then the mesh created by the
library.
The output of the generator is not ideally suited for for a pipe. The code gen-
erates a mesh that is on a 2-dimensional plane as in Figure 5-10. The mesh must
then be post-processed in order to have all the desired properties. The holidays
are dealt with first. The output along the edge of the holiday is a piece-wise set of
straight lines. Since the elements are quadratic in nature, a better job can be done
by moving the center node of each line segment onto the circle used to define the
holiday. As seen in Figure 5-3, the quadratic nature of each side of the element
Figure 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.
Figure 5-11 shows the locations of the center nodes on each side of an ele-
ment by subdividing each element up into four smaller elements using the center
nodes. Segments of the holiday boundary are seen to bee straight lines consist-
ing of two segments. After the nodes at the center of the lines are moved to their
proper positions, the circle is defined as best that can be done with a fixed number
of quadratic line segments see Figure 5-12. The movement can be done so long
as the angle at any corner of the triangle is greater than 0 degrees. In fact, it is de-
sirable to have the angles of all the corners of the elements be greater than about
15 degrees.67 This avoids numerical difficulties with large aspect ratio triangles
when blind Gaussian integration is used. The problem is greatly alleviated using
Figure 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.
Figure 5-13: A gouge at the bottom of a pipe that was meshed using rectangular
elements.
91
Figure 5-14: Closeup of one edge of a rectangular holiday where two triangles
have been used to increase the element density of the holiday.
the gouge, three elements are needed. Along the length of the gouge, a minimum
rectangular, it is best to use the rectangular elements to create the mesh. In this
case triangular elements are used to increase the element density at the holiday as
described in section 5.1.2. Two triangles are placed at either end of the holiday to
increase the number of elements across the holiday from one to three. See Figure
There are a number of problems placing holidays on the pipe. For both mesh-
ing practices described, there can be some distortion of the pipe surface due to
wrapping the planar elements to the pipe. It is most often seen as a ripple in the
surface near the edge of the holiday. Great effort has been made to minimize the
92
unwanted effect, but if too few nodes are specified to mesh the circumference of
When two pipes cross, attention must be paid to the nature of the mesh gen-
erated. Both pipes must account for, in their solution, the effect of the close ap-
proach of the other pipe. Pipeline operating companies can attest to the fact that
the pipes will greatly interfere with each other and that special precautions are
needed.
In order to account for the close approach of another pipe in the solution,
the mesh must be adjusted. There must be sufficient nodes located at the point of
closest approach, but too many nodes is a waste of computer time. For simplicity,
the same technique used for a change in coating properties is used. At the point
of crossing, a coating property change in placed to force the mesh generator to in-
crease the element density at that point, see Figure 5-15. With the addition of the
extra elements, the model has the flexibility in the solution to closely approximate
The boundary element method must have the boundary of changes in bulk
mains, the boundary that is to be divided into elements extends to infinity. Mesh
generators cannot create a mesh that extends to infinity with standard elements.
Therefore, a mesh is created that extends sufficiently far into the domain such
93
Figure 5-15: Mesh at a crossing of two pipes. Both pipes have the highest density
of elements at the crossing point.
that most of the variations in potential and current density have died out. At the
end of this mesh, a special element known as an infinite element is placed. The
Soil divisions may need to have holes in the mesh to allow pipes the pass
through. Since the location and number of pipes is arbitrary, the same mesh
generator used to create holidays is used.77 Sample output of the mesh generator
Pipes are added to the mesh such that the nodes of the pipe at the point of
crossing the soil type division, match exactly the nodes at the edges of the holes
Figure 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.
Figure 5-17: Completed mesh with pipes passing through soil type division.
95
5.6 Conclusions
The generation of the boundary element and finite element mesh is a critical
component for the solution of the differential equations describing cathodic pro-
tection of pipelines. The solution obtained form the BEM and FEM can only be
trusted if the underlying mesh is of sufficient quality such that the shape func-
tions can closely approximate the true solution. The next two chapter discuss
how errors using the type of mesh described in this chapter can be reduced and
There are a great many factors that affect the accuracy of a solution obtained
by the Boundary Element Method. These include the basis (shape) functions used
to approximate the solution on the finite element, the size and location of an
element, and the numerical method used to evaluate the integrals seen in the
Errors in the boundary element method come from several sources. They can
be discretization, integration,80, 81, 82, 83, 84, 85, 86 shape function choice, and matrix
the solution to the PDE on the element for which it is applied. Choosing a shape
function - element size pair that cannot represent the solution to a given error
level can greatly affect the solution throughout the domain. This is because the
boundary element method applied to any given point in the domain or on the
boundary is a function of all the points on the boundary, including the point
96
97
Integration presents one of the most difficult and time consuming aspects of
the Boundary Element Method. There have been many papers devoted to eas-
ing the difficulty and improving the accuracy of the results.80, 81, 82, 83, 84, 85, 86, 87, 88, 89
The reason for this is that the integration kernel (in 3D) contains a weak singu-
larity (1/r) and a strong singularity (1/r2 ) which must be evaluated in a Cauchy
ply used a basic Gaussian quadrature to evaluate the singular integrals.82 After
recognizing that the results of the integrals were the major source of error, sev-
eral techniques have been utilized to remove the singularity. For simple shape
may be used to evaluate the weakly singular integrals.81, 82 However, the strongly
singular integrals cannot be evaluated easily. All of the strongly singular inte-
grals are located on the diagonal of the H matrix (see section 4.1.1 for details of
this matrix). For these integrals there are two choices if the physics of the prob-
lem warrant. The first and easiest is to make use of a unit displacement or a
not elasto-plastic in nature (mechanics) or does not store charge in the domain
98
(steady-state), then the solution to the normal flux at the boundary is known to
HI = 0 (6-1)
where I is a column matrix all of whose elements are 1’s. Then the unknown
strongly singular integral, of which there is only 1 per row, can be solved explic-
The second technique is to solve the Cauchy principle value explicitly. There are
techniques for doing this for the direct boundary element method described by
83, 84
Guiggianni and for the Symmetric Galerkin Boundary Element Method by
90
the same author.
There are several ways to improve the accuracy of a given integral over an
arbitrary element. The easiest is to increase the order of the Gaussian quadrature.
This technique though, can never integrate a weakly singular integral since it
tions significantly. Cools et al.91 have made available general integration source
rules as “Cubature.”
Singular integrals can benefit from a change of the coordinate system to polar
coordinates. Polar coordinates reduce the order of the singularity by one. For
99
the 1/r dependence of the G integral, the singularity is removed as well as all
dependence on 1/r as one approaches the singularity. For the 1/r2 dependence
The final choice made for this project was to use a specially modified form
dard routines like IMSL could not be used.92 The routines that are used are based
take into account the high aspect ratios. The modified subroutine takes any arbi-
trary triangle and seeks to subdivide it such that the subdivisions are performed
in the area with the largest error contribution to the total error for that element. It
also tries to extrapolate the subdivisions with a nonlinear algorithm to reduce the
a value for the integral over a certain specified geometric shape and an estimate
of the error in the integration. This is accomplished by actually having two inte-
gration rules of different order, where the higher order formula is a super-set of
the other. This means that for a 13-15 rule, 13 points form the lower order rule
and the addition of only two more points to the rule gives the 15 point formula.
It is computationally cheap because the first 13 points of the 15 point rule are
identical to the 13 point rule. Adapting the integration routine for high aspect
100
Section 1 Section 2
Figure 6-1: Sub-element with a high aspect ratio subdivided by the adaptive in-
tegration routine to reduce the aspect ratio of the subsections.
Section 1 Section 2
Section 3 Section 4
Figure 6-2: Sub-element with a desired aspect ratio subdivided by the adaptive
integration routine.
ratio elements is very straight forward. After finding the subsection that has the
largest error contribution, the aspect ratio of the real sub-element is calculated.
If the ratio is greater than 2 or less than 0.5, the element is divided into two sub-
sections such that the aspect ratio of the sub-sections is reduced (or enlarged) to
bring it closer to the desired range. If the aspect ratio is already within the desired
A new technique has been introduced by Karafiat that could reduce the time
the vertices of a small number of sub-triangles to optimize the effect of the quadra-
ture. However, there is some ambiguity in the cost function that places these
nodes, in that a local minimum can be found without obtaining the global mini-
mum. Also, the technique was only described for a case where the singularity lies
in the same transformed plane as the element itself. For 3D models of pipelines,
the singularity may be located anywhere (i.e. above the center of the element) in
Singular integrals are handled by first insuring that the point of singularity
important for accuracy that the subsections that actually contain the singularity
have an aspect ratio of 1.0. For instance, suppose a quadratic Lagrangian element
is to be integrated, and the singularity is at the lower left corner. The element is
then broken into 3 separate integration domains. One that contains the singular-
ity and two the are “nearly singular” because of the proximity of the singularity.
the order of the singularity, while the other two subsections are integrated using
tions, Gaussian elimination is a very poor choice on computers. The reason for
102
element. That error is added to every other element in the column through the
row operation used to eliminate elements in the matrix A. This is repeated during
back-substitution which carries the error back up through the column to elements
tion used in the code is a standard subroutine from LAPACK named dgesv and
its dependencies.96 This subroutine makes extensive use of BLAS (Basic Linear
Algebra Subroutines).97 The source code for these libraries can be found on the
be linked together to perform almost any type of linear algebra function. There
is an added benefit in that the BLAS routines have been highly optimized for al-
most all computing platforms and are freely available. The optimized versions
can have a performance impact of 2 to 16 over just compiling the basic source
code. The libraries for Intelr processors running a Win32r operating system can
2000), Intel had not posted the libraries for the Linux/FreeBSD environments yet
After the code has been written to solve a differential equation, it is important
to benchmark it against known solutions. The code used herein solves Laplace’s
equation in a half space. There are two analytic solutions to Laplaces equation
for a half space that have been used to verify different portion of the code. The
first is due to Newman where a disk electrode is placed on the plane defining a
half-space and the primary current distribution is found. This solution provides a
in the solution that causes the current density at the edge of the disk to approach
infinity. Most numerical solutions have a great deal of difficulty obtaining this
solution.
The second analytic solution is due to Kasper.8 This solution is directly rele-
vant to pipes because it solves Laplace’s equation for two parallel infinitely long
cylinders. This solution shows the importance of the solution around the cir-
cumference of the pipe and can give an indication of the “fineness” of the mesh
A third related solution is from Dwight.3 His equation is used to find the re-
103
104
i 0.5
=p (7-1)
iavg 1 + (r/ro )2
where ro is the radius of the disk, i is the current density at a distance r from
the center of the disk and iavg is the average current density for the entire disk.
The disk is placed on the plane that defines the half-space. The counter-electrode
must be at a great distance so that it does not interfere with the current distri-
bution. It also should have a greater surface area than the electrode so that the
The first comparison to the analytic solution given in equation (7-1) is for a
model setup that is nearly identical to that used in obtaining equation (7-1). The
model setup is shown in Figure 7-1. The model used a disk of 15.24 meters diam-
eter. The anode was placed 500 meters from the disk. Hydrogen evolution Tafel
kinetics were used for the boundary condition. At a sufficient current level, the
For hydrogen evolution, the numerical calculation approached the primary dis-
tribution when the average current density was about 0.75µA/cm2 . The results of
the calculation and the analytic solution are plotted together in Figure 7-2, where
the solid line represents the analytic solution. The results show excellent agree-
105
10
i/i avg
0.1
0.001 0.01 0.1 1
1-r/ro
ment with the analytic solution over all but the edge of the disk. The difficulty at
the edge is that a numerical method cannot approximate infinity very well. When
the solution is placed on a log-log plot, a straight line should be seen as 1 − r/ro
approaches zero. A good approximation to the straight line is made when the the
mesh is sufficiently fine at the edge. This fact must be taken into account when
disk shaped holiday on a pipe also responds as theory predicts. A holiday was
placed on a 10 mile long pipe of 18 inches diameter. A portion of the pipe with the
holiday is shown in Figure 7-3. A single anode was placed 150 meters from the
pipe. The anode output was adjusted to obtain certain values of the parameter J
107
1.1
i/iavg (center)
0.9
0.8
0.7
0.6
0.5
0 2 4 6 8 10
J
2.303 ro iavg
J= (7-2)
βc κ
where βc is the cathodic tafel slope for the hydrogen evolution reaction, ro is the
radius of the holiday, iavg is the average current density on the holiday, and κ is
half-plane. The results of the calculation are shown in Figure 7-4. The mesh for
the holiday used two elements in the radial direction from the center of the hol-
iday and eight elements in the angular direction with a total of 15 quadratic tri-
angles (see Figure 10.12(b) for an example of the mesh). The results demonstrate
cylinders that are not concentric. The solution assumes a constant potential dis-
tributions around each cylinder and is given in terms of the ratio of the maximum
Kasper also assumed that the cylinders were infinitely long and therefore
dimensional model, an aspect ratio for the cylinders of 1000 to 1 was used. The
ability of the 3 dimensional model to handle such an aspect ratio and provide an
A plot of the current distribution for the case of the cylinders being placed
with a separation of 3 time the diameter is given in Figure 7-6. The ratio of the
developed as an anode resistance formula. This means that one calculates a re-
sistance with this formula by using the resistivity of the surrounding soil and the
surface area of the anode/pipe. After the resistances of a pipe - anode pair is
obtained through the use of the basic form of Ohm’s Law, V = I R, where V is the
potential shift of the pipe from the anode and R the the series resistance of the
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0 1 2 3 4 5
Distance between cylinders in terms of the diameter
1
0.95
0.9
0.85
0.8
0.75
0.7
0.65
0.6
0.55
0.5
0 0.1 0.2 0.3 0.4 0.5
Angular Distance (Normalized)
pipe and anode are remote from each other. When the pipe and anode are close
to each other (5+ diameters8 ), the distributions have a significant effect on the
actual current flowing between the anode - pipe pair. When the anode is outside
this limit, the current driven to the pipe is higher than that predicted by Dwight’s
The solution seen in Figure 7-7 was obtained by setting the potential distribu-
tion to be constant on both the anode and pipe and then calculating the current
density distribution and integrate to get the total current. The constant potential
distribution was used to remove the resistance associated with the kinetics from
the solution. These are some of the same assumptions made in deriving Dwight’s
equation3 . The others are that the pipe and anode are effectively at an infinite dis-
tance from each other. This ensures the current is strictly ohmic controlled. As
111
1 100
90
58 ft Pipe/Anode
Dwight’s Equation 80
Total Current (Amps)
Error
70
60
0.1 50
40
30
20
10
0.01 0
1 10 100 1000 10000
Separation (pipe diameters)
Figure 7-7: Error in Dwight’s equation as a function of the separation of two cylin-
ders. 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.
112
can be seen in Figure 7-7, the model approaches Dwight’s solution as the pipes
Figure 8-1: A screen shot and movie (click on image above) that shows some of
the interactive graphics techniques needed for three dimensional geometry.
The interface to the model has been designed as a Microsoft Windowsr ap-
terface ”100 and “The Windows Interface Guidelines for Software Design”. 101
The guidelines describe how the user should interact with the data using the
mouse and keyboard. However, these two books do not provide information in
113
114
Figure 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 Microsoftr Win-
dows.
interaction with 3-dimensional data. For this, other sources are needed and are
The interface starts with some common elements that most programs use and
with which most users would be familiar. First there is the main window as
shown in Figure 8-2. The window contains both a standard and application spe-
cific tool bar, a menu, and at the bottom, a status bar. The standard toolbar sup-
ports common functions such as opening and saving models. The top title bar
has the standard button to control the size of the window. The main part of the
window is where the models and their interface elements will be displayed.
A model is displayed in the main window when either a new one is created,
or one is opened. The model windows is divided into two parts (see Figure 8-3).
The left part has dual purpose. The first is calculation control and feedback, the
other is for controls related the the right section. The right-hand section is for
Figure 8-3: A screen shot of the basic window with an empty model. The window
has two parts: the first is for control and feedback, the other is for data display in
3 dimensions.
A solution of a PDE is only useful if one can examine it and derive useful
conclusion from it. One of the simplest ways to examine a solution is to plot
it versus position in the domain. This works well for a 1-D problem, but what
solution of a 2D problem on the z axis for a rectangular domain, but these types of
plots are difficult to use and precise information cannot be obtained from them.
aries. The boundaries are then colored using a color scale that corresponds to the
solution. Vectors can be overlaid on the surface to indicate flux directions (see
116
Figure 8-1). This technique requires the ability to interact with the data since not
Two dimensional problems are straight forward when solutions are obtained
by the boundary element method. Since the techniques solves for the unknown
essential or natural boundary conditions, a plot of the solutions along the bound-
The most meaningful plots of the solutions are either one along the length
(axial) of the pipe at a particular angle, or a plot along the circumference at some
point along the length. Extracting the necessary data involves interpolating be-
The first case, which is shown in Figure 8-4, is relatively simple. Here, one
simply determines the desired angular position to display and then interpolate
between the nodes where the solution is known using the shape functions of the
The second case is significantly more complicated. Again, one should extract
the data for the plot such that the angle from the the top of the pipe does not
change as one moves along the length extracting the data. However, as seen in
Figure 8-5, the mesh may not be regular along the surface of the pipe. Also, the
density of nodes changes in the region of the holidays, meaning simple book-
keeping that can be used to make a plot from the data in Figure 8-4 cannot be
done.
117
Figure 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.
118
Figure 8-5: The 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.
119
plot data with several extra bookkeeping structures that must be created at the
time the mesh is generated. When rectangular holidays are used, the mesh is
generated using successive circles of nodes. If the center point, and angle of each
circle of nodes is saved along with the number of nodes that make up each circle,
then once the angular position on the pipe where the plot data is to extracted
is calculated, the plot data can be extracted one node circle at a time. At each
node circle, a calculation must be done to find the element the plot line passes
through and then the homogeneous coordinates for that point. The homogeneous
coordinates can then be used with the same element to get the desired solution
8.4 3D Visualization
and perspective.102 If the video card is fast enough, it also can be made to be
interactive. This means that a user of the model can use a mouse to change the
viewing angle and distance from the object of interest in a true interactive fashion.
The code has been written using interactive techniques described by Hall and
describe how to take the viewer into account when designing the interface. For
instance, using the mouse to control how objects move within the display is crit-
ical. The mouse is only capable of a 2-dimensional motion, but the objects must
this software, depth control. By selecting the mode of motion, movement of the
Hall and Forsyth also supplied a toolkit to learn the techniques.104 The code in
the toolkit was not used for the development of this software. Only the OpenGL
API (application programing interface) was used.105, 106 This meant that the code
would be portable across platforms that support the OpenGL API. However,
code for other parts of the graphical user interface are not portable. Therefore
the code was separated into two distinct sections: portable code which contains
the geometry engine and most of the OpenGL drawing primitives, and non-
portable code which contains the interface elements such as windows, menus,
OpenGL can draw only primitive objects made up of planar triangles, poly-
gons (subdivided into planar triangles), lines and points. The boundary element
method and finite element method also uses polygons to represent the geome-
try so OpenGL is ideally suited to the task. Additionally, OpenGL can shade the
polygons, lines, and points using garuad shading based on colors specified at
the vertices of the primitive. Garaud shading is a determining the color at any
point in a triangle by linear interpolation of the colors at the vertices using the
Figure 8-6: Shaded mesh of boundary element calculation grid. OpenGL allows
instant visual inspection of the mesh for correctness.
Drawing the objects as they are created (meshed) provides the user with an
excellent method for feedback on its correctness. Then, if an object has been
misplaced, the user knows immediately and can fix the problem before time is
OpenGL uses linear interpolation schemes and is only fast when linear poly-
gons are used. This means that the polygon lies in a single plane in cartesian
coordinates. To alleviate this problem somewhat, the grid used for the calcula-
tion can be divided up into smaller sub-units over which the linear interpolation
can be applied. For instance, a nine-node Lagrangian element can be divided into
four four-node polygons each of which has a linear interpolation done, and are
Figure 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 graphics elements.
Figure 8-6 shows a view of a pipe with each quadratic Lagrange element di-
vided up into four linear four-sided polygons. This method does not provide a
true estimate of the paraboloid that describes the curvature of each element.
A better way to visualize the true variation of each element in both geometry
and values of potential and current when restricted to piece-wise linear 2D sur-
faces is to break each real element into 16 linear surfaces. Such a subdivision is
shown in Figure 8-7. Here, one can see how well a small number of quadratic
Most of the code for the program was written as a set of C++ classes. Follow-
ing the philosophy of C++, objects should be self contained black boxes when
complete. One of the necessary functions of all geometry objects is that they can
draw themselves into the OpenGL rendering engine. Therefore the base class
handles drawing tasks for all derived geometry objects. The single argument
a wireframe and if shaded, which color set to use (potential, current, corrosion
current, etc.) This bases class maintains the list of nodes and elements used to
Before Render is called, the OpenGL engine is told how all polygons sent to
glPolygonMode(GL_FRONT_AND_BACK,GL_LINE);
Other commands can tell the engine to fill the polygons with a color, or to inter-
polate the color at every point within the polygon by the colors at the vertices.
Setting up the OpenGL engine is platform dependent and one should refer to
the appropriate literature for their platform. Because of this, the setup for this
software is contained in the interface code. Appendix B has the setup necessary
of ways. The first is to use an object that the field engineer will be familiar with,
a Cu/CuSO4 reference electrode. The electrode can then be moved over the sur-
face of a pipe to get a reading of the desired quantity. Unlike a real reference
electrode, a software version can output more than just a potential. The reference
electrode used in model provides on-potential, off-potential, net current, and cor-
rosion current.
CHAPTER 9
COMPUTATIONAL PERFORMANCE
Computer performance for the boundary element method is a very active area
of research.107 Since most of the computational time is spent evaluating the inte-
grals in equation (3-9) on page 32, it is the best place to look for performance
improvements. One of the best ways to improve the speed after algorithms
cheaper to hook several single cpu machines together with a network than to buy
a multi-CPU machine and that such arrangements are readily available in most
There are several ways to break down the Boundary Element Method into
The smallest tasks that can be executed in parallel are the calculation of the
one point until all the points have been calculated. For this method to be effi-
cient, an extremely high speed bus, cache and memory must be present or all of
the CPUs may have to wait for one CPU to put its result into memory to allow
calculation of the final value of the integral can be calculated. This method is very
efficient on vector type computers(Cray C-90, Y-MP) and on shared memory pro-
cessor arrays (SGI Origin 2000, Intel Paragon). Each processor caches its gauss
125
126
point data and instructions and simply receives the physical element where the
point will be evaluated. This allows for efficient use of CPU power but requires
The next level up in the size of a parallel task is for each CPU to perform an
entire single integration. In this scenario, each cpu receives a set of data necessary
to complete a single integral. This integral forms a single entry in one of the
matrices in equation (4-2). This method also works well in small processor arrays
and for small distributed memory computers. The downside of this method is
that the ratio of time spent communicating between processors to the time spent
in calculation higher.
now a matter of determining when the communication takes place. One can as-
sign a CPU a block of point-element pairs and send all the necessary data at the
pairs immediately after each set of integrals have been completed and results
returned.
These types of computers are made up of a single large block of physical mem-
ory with a large number of processors having access to the memory. The SGI
127
Origin 2000 and Intel Paragon are examples of these types of computers. They
have the potential of being the fastest of all types of supercomputers. Generally,
dard networking system. Instead of being able to directly address another proces-
Another problem, is that the memory is also distributed, so any data needed by
the a processor must also be sent along with the calculation commands. This
means that the level at which the parallelism is made must be at a coarser level
The point at which the parallelism is applied depends greatly on the overhead
ers making up the parallel machine. It is important that the overhead of commu-
nicating is less than the time required to perform the calculation or the gain of
having multiple CPUs is lost.108 Many people have solved this problem by sim-
ply assigning each machine a large block of the matrices that appear in equation
(4-2). The method assumes that all the CPUs have the same performance capabil-
ities and that no other processes are running. It also does not allow for any error.
If for some reason one of the computers crashes, shuts down or loses contact with
There are several codes available to handle the data and command communi-
These include Parallel Virtual Machine (PVM), and Message Passing Interface
(MPI). For the work done here a new approach has been taken. The new method
makes use of a faster and lower level interface for inter-computer communica-
tion called Remote Procedure Call (RPC). RPC allows a computer program to call
C contains details for writing RPC code to separate functions onto several com-
puters.
cate data to each other. The communication can be greatly simplified if takes a
a star configuration with the master computer at the center of the star.
In a star configuration, the master computer takes the input from the user and
constructs the mesh for the problem. Boundary conditions are supplied by the
user. The master computer then splits up the calculation according to the parallel
method being used and sends tasks out to each computer in the star. When the
outer point computers complete the assigned task, they return the results to the
central computer. The central computer compiles the various result to form the
complete solution. Code for the communication on both the client and server
1200
400
0
0 400 800 1200
Expected Value (MFLOPS)
This is done by installing the remote module on the stand-alone machine and
telling the main program to use its own remote module. A simple test gave the
results in Table 9.1. The reason for the difference is fairly simple. The amount
of data needed by the small integration module is very small and can usually
130
fit within the CPU cache. This is referred to as good locality of data. The main
program also has to keep track of where in the main matrices the data returned
used to determine the efficiency of process and the computational power applied
puter was timed using a random selection of integrals to establish what its maxi-
mum stand-alone performance would be. The results have been normalized and
are shown in Table 9.2. These results have been checked by using them to pre-
dict the time required to solve a real problem and have been found to be accurate
to within 5%. By normalizing the result with respect to the master computer,
any differences that may have occurred because of the selection of integrals is
removed. Since the total FLOPS used to complete the calculation are known, the
tual model of a coated pipeline with coating holidays and a 5 anode ground bed.
The total number of floating point operations needed to evaluate all of the inte-
Table 9.2: Normalized benchmark result data for many types of computers.
Machine Processor OS and RAM Result
1 Pentium II 300Mhz NT 4.0, 384MB 1.0
2 Alpha 21164, 300Mhz NT Server 4.0, 128MB 0.976
3 Pentium II 300Mhz NT 4.0, 128MB 0.923
4 Pentium, 133Mhz NT 4.0, 32MB 0.298
5 Pentium 84Mhz NT 4.0, 32MB 0.088
6 Pentium 200Mhz NT 4.0, 32MB 0.439
7 Pentium 225Mhz NT 4.0, 32MB 0.494
8 Pentium 200Mhz Windows 98, 32MB 0.411
9 Pentium 133Mhz NT 4.0, 32MB 0.296
10 Pentium 100Mhz NT 4.0, 32MB 0.225
11 Pentium 100Mhz NT 4.0, 32MB 0.297
12 Dual Pentium 200 NT 4.0, 128MB 0.806
13 Pentium 133Mhz Window 95, 48MB 0.292
14 AMD K6 300Mhz NT 4.0, 64MB 0.419
15 Pentium 83Mhz NT 4.0, 40MB 0.100
16 i486 66hz NT 4.0, 32MB 0.030
17 Pentium 120Mhz NT 4.0, 48MB 0.296
18 Pentium 133Mhz NT 4.0, 64MB 0.296
19 Pentium Pro 200 NT 4.0, 128MB 0.553
20 Pentium 200Mhz NT 4.0 Server, 64MB 0.450
21 Pentium II 233Mhz NT 4.0, 64MB 0.747
22 Pentium II 233Mhz NT 4.0, 128MB 0.773
23 Pentium 200Mhz NT 4.0, 32MB 0.407
24 Pentium 200Mhz NT 4.0, 64MB 0.442
25 Pentium II 233Mhz NT 4.0, 64MB 0.762
132
A timing mechanism was built into the software to check the wall-clock time
before and after the calculation was done. Using the time and the number of float-
ing point operations needed, an average through-put for the distributed com-
The procedure used to build a performance scaling line was to run the model
first on the master only to get the first point in the line. Then machines were
added to the distributed computer, 1 at a time, with the model being run once
as each computer was added. To eliminate any effect that the order of adding
computers may have on the shape of the line, several runs were made of adding
all of the available computers using a different random order of machine addition
Often, with distributed machines, the actual node in the machine may be used
for some job by another user. This effect can be seen in the center section of rows
in Table 9.3. Those machines that are under heavy use by other users contribute
very little to the distributed computer. On the other hand, if there is very little
external activity, the performance scales very well. An abbreviated table of data
is shown in Table 9.4. Figure 9-1 shows the pooling of several runs of the type
shown in table 9.3. In this case, linear scaling through all of the machines is
seen up to about 800 mega-FLOPS which corresponds to 12.5 times the speed
of the master computer. Scaling is still good through the tested range of 1100
It is apparent that this method will scale up to the greatest amount of free
CPU time available on all of the machines within a distributed computer. This
Loading and saving data from persistent storage can be a source of slow per-
degrade by over a factor of 100. It is also advantageous to have only one of the
matrices (H or G) in RAM at any time and to have both matrices only in persis-
tent storage when forming the A matrix. This can increase the size of the solvable
√
problem by a factor of 2.
is best to store these matrices in the FORTRAN style column major order. This is
opposed to C/C++, where the natural form is row major. Some experiments were
done using a small set of matrices and a large set. The first experiment involved
135
Table 9.5: Conditions used for all hard drive performance tests
Computer Pentium Pro 200
RAM 128 MB
Hard Drive 2GB 5400RPM Ultra SCSI
Small H matrix size 1.2 MB
Small G matrix size 2.8 MB
Large H matrix size 13.7 MB
Large G matrix size 30.7 MB
where the matrices H and G are stored in persistent storage. Table 9.6 shows the
results of the time required to calculate the residual. The conditions for all hard
The next test was to build the matrix A as described in section 4.5 on page
62. The results are shown in Table 9.7. Here it is apparent that a substantial per-
The performance improvement get better as the matrices increase in size until the
optimum block size for data transfer from the hard drive is reached. The value
varies from drive to drive but is usually between 16K and 64K. 64K corresponds
for their processors. Intelr has built a set of well tuned libraries as made them
freely usable in compiled form.98 These libraries have the added benefit in that
they can utilize more than one processor on the machine where the code is run.
Intel’s data shows linear scaling of CPU efficiency on up to four processors. The
code has been tested for this work on single and dual processor machines and
Intel also supplies optimized math functions such as ln, exp, sin and cos. A
factor of 2 improvement over the stdlibc routines that accompany all compilers
can be seen.
CHAPTER 10
HOLIDAY EFFECTS ON SINGLE PIPES AND PIPELINE NETWORKS
The goal of this research was to be able to predict the effects seen in current
and potential distributions in pipeline networks. One contributing item was coat-
ing holidays. Coating holidays are a small sections of the pipe where the coating
has somehow been removed from the pipe and bare steel is exposed to the soil
Figure 10-1: Pipe that exhibits many coating holidays. Pipe was coated in the late
60’s and dug up in 1996.2
109 110
The early work by Kennelly et al. and Orazem et al. on modeling the
influence of coating holidays employed the Finite Element Method. While their
results showed that the majority of current from the CP system went to the hol-
iday, there was a significant error in the closure of current (3% to > 12%) when
modeling holidays.109 Thus, their work demonstrated the weakness of using the
classical finite element method where the flux is calculated from potential gradi-
137
138
role of coating holidays.33, 111, 112 The geometry considered was a single pipeline
protected by sacrificial ribbon anodes such as those used for the Trans-Alaska
the current near the coating-holiday edge. Their work supported the conclusion
that it is much more difficult to drive current to a holiday than to a uniformly poor
coating. The results of the model for a pipe protected by ribbon anodes showed
that the variation of potential around the pipe at a holiday can be as large as 550-
650mV. This value is the difference between the most cathodic potential on the
coated pipe and the most anodic potential at the holiday. Numerical difficulties
limited the application of the model to short sections of pipe (<40 ft.).
The first effect of holidays to investigate is to see how the current and po-
Orazem et al. 110 first demonstrated the effects. Figure 10-2 shows a large holiday
system operating at -10 volts from the connection to the pipe. At the edge of the
Cu/CuSO4 which is marginally protected. The holiday has the greatest potential
at the center, -863 mV Cu/CuSO4 . It might be expected that the coated region
close to the edge of the holiday is in the worst shape. However, the corrosion rate
is what is most important. Examining the corrosion rate at the center of the holi-
139
Figure 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 cri-
terion while the holiday does. Even so, the corrosion rate at the holiday is much
higher than at any of the coated sections.
day reveals that it is corroding at a rate that is 3 orders of magnitude higher than
that under the coating. The reason for the low rate under the coating is that the
coating causes the rate of oxygen transport to be decreased and it adheres to the
the transport of electrolytes through coating, and how water uptake affects pore
size, adherence properties and transport.27 They have shown that electrochemical
reactions (corrosion, oxygen reduction, etc.) are taking place at the steel surface
still observed with larger variations in the off-potential of the coating. A plot of
the angular potential and current distributions around the pipe at the holiday
is given in Figure 10-3. In this case the off-potential of the coating exhibits a 50
140
-1000 100
Holiday Location
-1010
10
µA/cm )
Potential (mV Cu/CuSO4)
2
-1020
Current Density (µ
1
-1030
0.1
-1040
0.01
-1050
-1060 0.001
0 60 120 180 240 300 360
Position (degrees)
Figure 10-3: Current density and potential distribution around the circumference
of a pipe through a holiday. Zero degrees is the top of the pipe. The blue line is
the current density.
current at the edge of the holiday. The reason for the draw down in current is
that a galvanic couple forms between the coated steel next to the holiday and the
exposed steel of the holiday. The coated pipe is at a more negative potential than
the bare steel because, among other things, the coating blocks oxygen transport to
the underlying steel thereby causing a shift in the corrosion potential. The effect
shift the potential of bare steel a given amount than coated steel.
To investigate situations involving two pipes, a model was set up with two
10-mile-long pipes running parallel to each other in the same right-of-way. Both
pipes were protected by a single remote groundbed placed at the 2.5 mile point.
141
-550 10
Holiday Location
µA/cm )
2
-650 0.1
Current Density (µ
-700 0.01
-750 0.001
-800 0.0001
-850 0.00001
0 60 120 180 240 300 360
Position (degrees)
Figure 10-4: Current density and potential distribution around the circumference
of a pipe through an under-protected holiday. Zero degrees is the top of the pipe.
The blue line is the current density.
The coating of the flawed pipe was assumed to be aged. The adjacent pipe was
Figure 10-4 shows the potential and current density distribution around the
circumference of the pipe at the point of the holiday. The holiday is under-
protected at the center and is corroding at a rate of about 0.5 mils per year. The
current density near the holiday is an order of magnitude smaller than that seen
The adjacent pipe shows some effects of the neighboring coating holiday. Fig-
ure 10-5 shows the off-potential and current density distributions for the well-
coated pipe. The effect of the holiday is apparent but it is not to the extent of
that seen for the pipe with the holiday. The effect is called shadowing. The cusp
seen at the distributions is a numerical effect due to the joint between elements
occurring exactly at the point of closest approach. The effect can be removed by
142
-835 0.0002
µA/cm )
0.00019
2
-840
Current Density (µ
0.00018
0.00017
-845
0.00016
-850 0.00015
0 60 120 180 240 300 360
Position (Degrees)
Figure 10-5: Current density and potential distribution around the circumference
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.
changing the number of elements used to create the mesh around the circumfer-
ence of the pipe such that the center point of an element is at the point of closest
approach.
20 feet away from the holiday, the potential distribution on the well-coated
pipe is much closer to uniform and all points around the circumference are pro-
tected. Even though the section of pipe closest to the holiday has fallen below the
tems was developed Section 4.3. A pipeline network was setup as shown in Fig-
ure 10-6. For this system, two pipes were placed adjacent to each other with one
143
pipe leaving the right-of-way at a 45 degree angle 1 mile from the starting point
of the model. The diverging pipe was given two different coatings of the form of
equations (1-24) and (1-25). One was for the section that ran parallel to the second
pipe with parameters to emulate an “aged” coating. The other section consisting
of the part diverging away from the second pipe was given a new coating similar
to the original. The second pipe was considered to have a new coating except
for a holiday at the bend in the pipe. The two pipes have separate CP systems
consisting of distributed groundbeds that are relatively close to each other. This
scenario can be thought of as an old oil pipeline and a new gas pipeline sharing a
portion of a single right-of-way. The model was used to examine the point where
A small portion of the pipes being modeled is shown in Figure 10-7. The shad-
ing represents the calculated potential distribution around the pipe. The stray
current associated with the two CP systems caused a portion of well-coated pipe
144
seen more clearly in Figure 10-8 where the current distribution around the well-
potential distribution is given in Figure 10-9 for a section of the pipe exhibiting
stray current effects and for a section of the pipe farthest from the section of the
pipe that exhibits stray current effects. In the region where stray current effects
are evident, the entire pipe falls below the -850 mV Cu/CuSO4 protection cri-
terion. The point closest to the neighboring poorly-coated pipe is anodic, and
the section opposite the neighboring pipe is at only slightly cathodic potentials.
These results can be compared to the current distribution given in Figure 10-8
Even in the section of the pipe that was not driven to anodic currents by the
interference between CP systems, the pipe is not adequately protected in the re-
gion that is closest to the neighboring pipe. For this system, rectifier settings
145
-1
0 90 180 270 360
Position (degrees)
Figure 10-8: Current Density around pipe in area exhibiting stray current in Fig-
ure 10-7. Positive current densities are anodic.
-500
Potential (mV Cu/CuSO4)
-600
-700
Point farthest from
neighboring pipe Point closest to
-800 neighboring pipe
-850 mV Cu/CuSO4
-900
-1000
0 90 180 270 360
Position (degrees)
could be found that would cause the entire system to meet or exceed the -850
was over-protected. Satisfactory protection for this case was achieved only when
the pipes were bonded such that stray current effects were eliminated. Again, the
smooth distributions presented in Figure 10-8 and Figure 10-9 reflect the accuracy
of the meshing and the adaptive integration schemes used in these calculations.
Soil surface on-potential distributions were calculated for the pipes exhibiting
stray current described in section 10.3. If there is going to be much to see in an on-
potential survey of the pipe, the results of section 10.3 should show it. Figures 10-
10 and 10-11 show the on-potential distributions. Figure 10-11 has been split so it
is apparent where the pipes are located with respect to the soil surface calculation.
On-potential surveys near the anodes tend to only show the anodes. This is made
147
Figure 10-11: Distant 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.
apparent in the calculations. The interference of the anodes masks the behaviour
of the pipes and even the stray current situation. The only indications of the pipe
under-protected large coating holiday next to a small holiday would adversely ef-
fect the small holiday. To test the idea, a model was setup with two holidays very
close to each other (see Figure 10-12). The large holiday exposed 17.7 cm2 bare
steel. The small holiday exposed 12.6 mm2 of steel. The holidays were placed
with 2.91 cm space between them. This was the closest that two holidays could
be placed using the mesh generator. Both holidays were assumed to have the
148
Figure 10-12: Setup for calculation of a pinhole holiday next to a large holiday: a)
the dimensions of the holidays and b) the automatically generated mesh for the
holidays.
same polarization curve (equation (1-23)) and parameters. The CP system was
set such that the large holiday was under-protected. Figure 10-13 shows the off-
potential distribution of the two holidays. From the color scale, it is apparent that
all but the edges of the large holiday are under-protected. The small holiday re-
Therefore, the situations that should concern pipeline operators are large coat-
ing holidays. These are the most difficult to protect so monitoring systems need
to indicate when a problem will occur on a large holiday versus a small one.
Chapter 11 discusses one monitoring technique and how the size of holidays af-
Coupons are a small piece of pipeline grade steel placed in the same soil envi-
ronment as the pipeline. The coupon is connected to the CP system through a test
station so the polarization level of the coupon can be monitored without affect-
ing the CP provided to the pipe. It is thought that a coupon represents a holiday
of similar size. Mathematical models were used to explore the conditions un-
tive assessment of cathodic protection was obtained when the area of the coupon
was larger than the area of the largest expected holiday. The performance of the
11.1 Introduction
can be difficult to ascertain. Normally, all that can be obtained is the off-potential
of the pipeline itself which is said to represent some kind of a weighted average
of the true off-potential of all the holidays. When one cannot remove all sources
cating factors include direct bonded galvanic anodes, telluric currents, and stray
currents.
150
151
large holiday if there are many smaller holidays also on the pipe which sum
to a similar area as the larger holiday. One solution to this problem is to use a
structures.113, 114, 115 The coupon can be connected to the pipe or CP system through
a test station in such a way that all sources of current to the coupon may be turned
off without disrupting the current to the pipeline. It has been recommended
ing protected.3, 114 However, there has been little theoretical work done to show
Stears et al.114 and Lawson and Thompson115 showed that the off-potentials
measured for coupons varied with location relative to the pipeline. This result
current density for reduction of oxygen has a critical influence on the corrosion
tion will influence the coupon off-potential. Coupon configuration, location, and
size can also influence coupon readings. Lawson and Thompson115 observed that
small coupons tended to be polarized more than large coupons. This observation
is consistent with the calculations presented by Orazem et al.33, 109, 110, 111 Lawson
and Thompson115 also reported that in some cases, the coupon off-potential was
significantly more negative than the off-potential of the pipe, suggesting that the
exposed surface area was greater on the pipe than on the coupon. The differ-
ence between coupon and pipe off-potentials was influenced by the amount of
152
fied. The difference between coupon and pipe off-potential was generally small,
but, in some cases, was on the order of 200 mV. The biggest differences were seen
for cases where the coupon reading of -780 mV indicated that a 200 mV potential
shift criterion was met but where the pipe off-potential of -580 mV suggested that
influence the coupon readings relative to pipe condition. This work takes a first-
coupon and the off-potential of a holiday through the solution of the governing
114, 115
differential equations. Recently published experimental work provides a
basis for comparison to the models. Their work indicates that there should be a
the setting of a rectifier and the change in the off-potential of a pipe with coating
It should be noted that, because the calculations provide detailed results for
individual holidays, the modeling approach presented here should yield greater
insight into the relationship among coupons and holidays than can be obtained
Parameter Values
ilim,O2 0.1, 1.0, 3.16, 10.0 µA/cm2
EFe -526 mV(CSE)
βFe 59 mV/decade
E O2 -104 mV(CSE)
β O2 59 mV/decade
EH2 -955 mV(CSE)
βH2 118 mV/decade
each coupon and holiday. A coupon with a new bright finish, for example, would
rent density for oxygen reduction; whereas, to account for the presence of scale
code can also handle many holidays on the same pipe, each sized differently and
The polarization curve used for bare metal surfaces (coupons and holidays)
was equation (1-23). The parameter values for equation (1-23) are given in Table
(11.1).
The coating was modeled using a set of equations given in Chapter 1. The
values in Table 11.1 were used for the polarization behavior of the steel beneath
the coating and the addition parameters listed in Table 11.2 determined from ex-
Table 11.2: Coating parameters used in model. Note: These parameters were
determined after water uptake of the coating.116
Parameter Values
A/ Apore 1000
ρ 5 × 109 Ohm-cm
δ 20 mil
αblock 0.99 (99% effective)
11.3 Approach
The pipeline used for the model was 6.44 km in length and had a diameter
of 457 mm. The pipe was protected by a single impressed current ground bed at
one end of the pipe section. The holiday-free sections of the pipe were assumed
to have a uniform 0.5 mm (20 mil) FBE coating. Coating holidays were assumed
to expose bare steel. The anode used for this model was an impressed-current
ground bed, sized at 76 mm OD and 15.25 m length. Its current output could
be adjusted until the desired reading of the off potential for the coupon was ob-
tained.
Two coupon designs were used for the simple cases. The first was a band
of bare steel of length 6.35 cm on a 5.08 cm diameter pipe which follows the
design described by Stears et al.114 The second was a 1.27 cm rod with an exposed
tip of 2.5 cm which follows the design described by Lawson and Thompson.115
Schematic representations of the coupons are given in Figures 11.1(a) and 11.1(b)
respectively.
Two sets of calculations were performed. One set involved a single holiday
with a single coupon. The coupon was placed 6 meters from the holiday with 1.2
m top cover over the exposed portion. For each of the coupon types, four ratios
155
(a) (b)
Figure 11-1: Schematic representations of the two cylindrical coupons used in the
simulations: a) large coupon114 with a surface area of ∼ 50 cm2 ; and b) smaller
coupon115 with a surface area of 9.5 cm2 .
of the holiday area to the coupon area were used. For each of these cases, four
different polarization curves were used. Finally the current output of the anode
was adjusted iteratively such that the average off-potential of the coupon was -
750, -850, or -950 mV CSE. The calculation matrix therefore included 96 different
cases.
The more complicated scenario involved the same pipe and anode, but the
pipe was assumed to have five coating holidays of different sizes. Three of the of
them were placed at the top of the pipe (0o ), and the remaining two were placed at
the bottom (180o ). The bottom of the pipe was assumed to rest in water-saturated
soil while the top was assumed to be in dry soil. To simulate the influence of
water level, the holidays at the bottom of the pipe were assigned a polarization
156
curve with a ilim,O2 that was an order of magnitude lower than that used for the
The design of the large coupon was changed for the more complicated sce-
nario. From the results of the first simulations, it was apparent that there may
nary exploration of the role of coupon design, the large coupon was redesigned
to match closely the design of the small coupon, Figure 11.1(b), but with the same
The difference between the measured potential of the coupon and the most
positive potential on the holiday is presented in Figure 11-2 for the large coupon.
The corresponding results for the small coupon are given in Figure 11-3.
potential of the coupon and the most positive potential seen on the holiday em-
phasizes the behavior of the system where the corrosion rate is the highest. The
sign on the y-axis indicates whether the holiday is at a more cathodic potential
(negative sign in front) or at a more anodic potential than the coupon. Both sets
As seen in Figures 11-2 and 11-3, the coupon does not provide an exact mir-
a change in the rectifier output results in a similar shift in the off-potential of the
coupon and the holiday. The reliability of the coupon off-potential as an indi-
1000 1000
ΦHoliday-ΦCoupon (mV)
ΦHoliday-ΦCoupon (mV)
100 100
10 0.1 µA/cm2 10
1.0 µA/cm2 0.1 µA/cm2
3.16 µA/cm2 1.0 µA/cm2
10.0 µA/cm2 3.16 µA/cm2
1 10.0 µA/cm2
1
-1000 -950 -900 -850 -800 -750 -700 -1000 -950 -900 -850 -800 -750 -700
Potential of Coupon (mV CSE) Potential of Coupon (mV CSE)
1000 1000
0.1 µA/cm2
−(ΦHoliday-ΦCoupon )(mV)
1.0 µA/cm2
−(ΦHoliday-ΦCoupon )(mV)
3.16 µA/cm2
10.0 µA/cm2
100 100
10 10 0.1 µA/cm2
1.0 µA/cm2
3.16 µA/cm2
2
10.0 µA/cm
1 1
-1000 -950 -900 -850 -800 -750 -700 -1000 -950 -900 -850 -800 -750 -700
Potential of Coupon (mV CSE) Potential of Coupon (mV CSE)
Figure 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 holiday is at a potential cathodic to the coupon.
The lack of a sign indicates the opposite. Panels a) through d) correspond to
different values of ξ = Aholiday / Acoupon .
158
ΦHoliday-ΦCoupon (mv)
2
10 µA/cm
100 100
10 0.1 µA/cm
2 10
2
1.0 µA/cm
3.16 µA/cm2
10 µA/cm2
1 1
-1000 -950 -900 -850 -800 -750 -700 -1000 -950 -900 -850 -800 -750 -700
Potential of Coupon (mV CSE) Potential of Coupon (mV CSE)
1000 1000
0.1 µA/cm2 0.1 µA/cm2
1.0 µA/cm2
−(ΦHoliday-ΦCoupon )(mV)
1.0 µA/cm2
ΦHoliday-ΦCoupon (mV)
10 µA/cm2 100
2
100 10 µA/cm
10 10
1 1
-1000 -950 -900 -850 -800 -750 -700 -1000 -950 -900 -850 -800 -750 -700
Potential of Coupon (mV CSE) Potential of Coupon (mV CSE)
Figure 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 holiday is at a potential cathodic to the coupon.
The lack of a sign indicates the opposite. Panels a) through d) correspond to
different values of ξ = Aholiday / Acoupon .
159
the off-potential of the holiday and the off-potential of the coupon. The sign of
the differences given in Figures 11-2 and 11-3 determines whether the coupon
provides a conservative assessment of the state of the holiday. For coupons that
are larger than the biggest holiday, the differences were negative, and the hol-
iday was more protected than the coupon. For coupons that are smaller than
the biggest holiday, the differences were positive, and the coupon was more pro-
parameter. While no specific critical value of ξ was found that gave consistently
protection applied to the holidays. To compare the effects the geometry has on
the performance of the coupon, similar ratios of the holiday area to the coupon
area can be used. By comparing Figure 11.2(a) and Figure 11.3(b), one can see the
design of the smaller coupon, Figure 11.1(b), results in a slightly closer match to
the off potential of the holiday at high values of ilim,O2 . As seen in Figure 11.3(a),
which represents a holiday with a surface area equal to the surface area of the
large coupon, the small coupon seriously over-predicted the level of cathodic
protection at the holiday. In this case, at high values of ilim,O2 , the holiday po-
tential was very close to the corrosion potential for conditions where the small
coupon had an off-potential showing a 280 mV shift from the corrosion potential.
size. If the coupons did behave this way, then all the plots in Figures (11-2) and
(11-3) would have their lines at a 0 difference between the holiday and coupon at
160
Table 11.3: Results of two coupon model. The small coupon used a polarization
curve with ilim,O2 = 1.0 µA/cm2 (curve 1). The large coupon used ilim,O2 that was
an order of magnitude smaller (curve 2).
Object Size Location Pol. Curve Pot. (CSE) Avg. Curr. Den.
Holiday 1 7.81 cm2 0o 1 -839 mV 1.11 µA/cm2
Holiday 2 12.9 cm2 0o 1 -834 mV 1.10 µA/cm2
Holiday 3 130.6 cm2 0o 1 -787 mV 1.05 µA/cm2
Holiday 4 64.5 cm2 180o 2 -847 mV 0.228 µA/cm2
Holiday 5 6.45 cm2 180o 2 -855 mV 0.244 µA/cm2
Small Coupon 9.36 cm2 122 cm 1 -850 mV 1.13 µA/cm2
Large Coupon 101 cm2 167 cm 2 -854 mV 0.241 µA/cm2
Table 11.4: Results of two coupon model. The large coupon used a polarization
curve with ilim,O2 = 1.0 µA/cm2 (curve 1). The small coupon used a ilim,O2 that
was an order of magnitude lower(curve 2).
Object Size Location Pol. Curve Pot. (CSE) Avg. Curr. Den.
Holiday 1 7.81 cm2 0o 1 -867 mV 1.19 µA/cm2
Holiday 2 12.9 cm2 0o 1 -861 mV 1.18 µA/cm2
Holiday 3 130.6 cm2 0o 1 -815 mV 1.08 µA/cm2
Holiday 4 64.5 cm2 180o 2 -873 mV 0.312 µA/cm2
Holiday 5 6.45 cm2 180o 2 -883 mV 0.350 µA/cm2
Small Coupon 9.36 cm2 167 cm 2 -884 mV 0.356 µA/cm2
Large Coupon 101 cm2 122 cm 1 -850 mV 1.13 µA/cm2
all potentials. The results show that the coupon does not yield the off-potential
the difference between the off-potentials of the coupon and the holiday and the
mass transfer limited current density for oxygen reduction (ilim,O2 ). The ratio of
the surface areas then indicates whether holiday will be cathodic or anodic to the
coupon.
The results of the simulations for the multi-holiday, two-coupon system are
summarized in Tables 11.3 and 11.4. The results shown in Tables 11.3 and 11.4
161
simulate the use of coupons for a pipeline with a distribution of holiday sizes and
with the water table at the center-line of the pipe. Thus, the bottom of the pipe
is in wet soil, and the top of the pipe is dry. The least conservative assessment
of cathodic protection was obtained with the large coupon placed at the bottom
of the pipe (Tables 11.3). In this case, both coupons met the protection criterion
of −850 mV(CSE), yet only the smallest of the five holidays met the protection
with the large coupon placed at the top of the pipe (Tables 11.4). Again, both
coupons met the protection criterion of −850 mV(CSE). Only the largest of the
five holidays did not meet the protection criterion, and this holiday was only 35
actions on the coupon was the same as that on bare steel exposed by coating
holidays, the most important parameters were ξ = Aholiday / Acoupon and iO2 , the
current density for reduction of oxygen. All the calculated results could be placed
yield off-potentials more positive than exposed steel are indicated by black sym-
bols, cases where coupons yield off-potentials more negative than exposed steel
are indicated by open symbols, and cases where coupons yield off-potentials less
than 15 mV more negative than exposed steel are indicated by gray symbols.
Calculations to the left of the line at Aholiday / Acoupon = 0.8 were considered
10
2
iO2, µA/cm
Conservative Non-Conservative
0.1
0.01 0.1 1 10 100
Aholiday/Acoupon
Figure 11-4: Coupon performance diagram for coupons placed in the same soil
environment as the pipeline coating defects that expose bare steel. The soil re-
sistivity was 10,000 Ω cm. Black symbols correspond to calculations where a
protected coupon corresponded to a protected 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 sym-
bols correspond to cases where the coupon reading is not conservative, but the
error was less than 15 mV.
163
will yield a conservative assessment of the pipe condition so long as each coating
holiday that exposes bare steel exposes no more than 80 percent of the coupon
area. The dashed line indicates the approximate boundary for which the degree
to which the coupon was not conservative was smaller than 15 mV.
increase in soil resistivity from 10,000 Ω cm to 100,000 Ωcm moved the position
Figure 11.5 did not provide a reliable indication of coupon performance for
cases where the polarization curve for the coupon corresponded to a less ag-
gressive environment than that experienced by the bare steel exposed by coating
that of the steel gave consistently non-conservative results, even for coupons with
areas 13 times the area of exposed steel. (Aholiday / Acoupon < 0.077).
It should be noted that the calculations presented here were for cases where
the only coating defects on the pipe exposed bare steel. The present work does
not address the use of coupons to assess the possible role of corrosion under intact
disbonded coatings.
11.6 Conclusion
coupons should provide a conservative value such that the largest expected hol-
iday is protected when the coupon is protected. The present work shows that
the area of the coupon should be larger than the area of the largest expected
ent sizes are to be used, the largest coupons should be placed in regions where
to the polarization curve indicates that coupons should be placed to sample all
posing bare steel. The coupon performance diagram addressed both polarization
Storage tanks are located at the terminus of most pipelines. They may be used
to store crude oil unloaded from ships, furnace fuel waiting for the fall demand,
or fuel depot reserves. Most cities have a fuel depot where gasoline, diesel, and
other types of petroleum products are stored before distribution to filling stations.
The storage facility allows pipelines to run continuously and efficiently, instead
a tank farm. They are built of steel and rest on the soil. The tank bottom is then
subject to the same materials problems as the buried pipelines. As the geometry
of a tank is much different than pipelines, tanks have different problems with CP
design than do pipes. Since tanks sit on the ground and are operated at atmo-
spheric pressure, the steel is thinner (0.25”) and, therefore, low corrosion rates
can be catastrophic.117
of underlying cache basins to prevent any leakage from entering the aquifer. The
liner of the cache basin is usually made of a very high dielectric membrane which
does not allow passage of electrical current. The CP systems for these types
of tanks are designed to place the anodes between the liner and the steel tank
bottom.118
165
166
Concentration profiles for oxygen under a storage tank are cyclic in nature.
The reason is that when the tank is empty, the center rises above the soil creating
a void. The void fills will air and because the tank is no longer in contact with
the soil, none of the oxygen in the soil at the center of the tank is consumed
through the oxygen reduction reaction (second term of equation (1-23)). When
the tank comes into contact with the soil again after filling, the mass-transfer-
limited current density for oxygen reduction is the same at the center as at the
periphery.
The model has been adapted to handle tank bottoms. It can also handle the
insulating cache basin if it is assumed that no current can pass through the liner
The boundary element mesh for the tank bottom was created using isopara-
The triangles are used to increase the degrees of freedom as the mesh proceeds
radially from the center. If the Lagrange elements where used exclusively, either
the center of the tank would be over-meshed, or the outside edge would have
a too-course mesh. The triangles are used to double the degrees of freedom by
placing a triangle at every other element on the last ring of elements that would
double the diameter from the previous spot where the number of degrees of free-
dom were doubled. Therefore, starting from the center, rings that use triangles
are 1 (the center which is all triangles), 2, 4, 8, 16, etc (see Figure 12-1). The den-
sity of elements is then set by specifying the density at the center and the total
167
Figure 12-1: High resolution mesh used for the boundary element mesh.
number of rings. The mesh generator takes care of the rest. The last ring is di-
vided into two rings to help handle the possible edge effects. The ring at the very
1 2
edge is 3
the size of the rest of the rings while the next to last is 3
the size.
Individual rings of elements were assumed to have the same parameters for
curve parameters does not mean that an individual ring will have the same po-
tential everywhere since different points in the rings can be at different points
mean that the concentration of oxygen in the soil near the ring was assumed to be
constant for that ring. The implicit assumption is that oxygen is consumed and
tion having an axisymmetric profile is valid if all the oxygen that diffuses to any
and that the diffusion constant for oxygen (DO2 ) does not change in any part of
Oxygen can get to the steel of a tank bottom with or without an underlying
cache basin through the soil around the edges of the tank. If it is assumed that the
oxygen is consumed according to equation (1-23), and that the concentration does
not vary with depth, then the consumption can be obtained through a balance
equation
∂ c O2
= ∇ · NO2 − RO2 (12-1)
∂t
where the flux of oxygen, N, is given by
where DO2 is the diffusion coefficient for oxygen transport through the soil. The
since virtually all of the oxygen that gets to the steel surface is reduced. This can
be written as
R O2 = kc (12-3)
169
where k is the proportionality constant which accounts for the thin diffusion layer
∂ cO2
= − DO2 ∇2 cO2 − kc (12-4)
∂t
Under assumptions of a steady state and that oxygen concentration does not vary
d2 c D dc
D 2+ − kc = 0 (12-5)
dr r dr
Upon moving the constants into the last term and multiplying through by r2
r !2
d2 c dc k
r2 2 + r − r c=0 (12-6)
dr dr D
which resembles a Bessel’s equation. To get the equation into the standard form
d2 c dc
ξ2 2
+ ξ − ξ2c = 0 (12-10)
dξ dξ
170
whose solution is a modified Bessel function. The boundary conditions for the
at r = R, c = c∞ (12-11)
and
at r = 0, c is finite (12-12)
and
at ξ = 0, c is finite (12-14)
The solution to equation (12-10) is the modified Bessel function of the first kind.
c∞
c(ξ ) = q I0 (ξ ) (12-17)
I0 (R Dk )
The solution is presented as a function of radial position in Figure 12-2 for several
p
values of the scaling parameter R k/ D.
171
1.00
R k
D =1
0.75
R k
D = 1.5
cO /cO ,∞
2
0.50 R k
D
=2
2
0.25
R k
D =5
0.00
0.00 0.25 0.50 0.75 1.00
ξ /R k
D
12.3 Implementation
which was solved using the boundary element method (BEM). The boundary
condition for the bare steel was given by equation (1-23). The anodes were also
condition at the steel surface of the tank bottom. The parameter, i O2 in equa-
tion (1-23) was adjusted according to equation (12-17). However, the value of
boundary element mesh was assigned a single value of the parameter (see Figure
12-1 for the element layout) which remained constant in that ring.
a different value at each ring introduced an error at the edge of the ring in the
current density. The boundary element method compensates for the error at the
center of the element. Therefore the point at which to sample the solution is
critical for this type of simulation. The best points to sample are the Gauss points
There were three anode arrangements used corresponding to the three stan-
dard anode arrangements used in industry.120 The first was a remote ground-bed,
ence of the tank and finally, ribbon anodes placed directly under the tank-bottom.
The remote ground-bed cannot be used when there is an insulating cache basin
12.4 Results
A single tank was modeled using the three anode arrangements described
above. The tank was 100 feet diameter with bare steel exposed to the soil. In
all cases the anode output was adjusted to 14.2 Amps. If the current density is
uniform for the entire tank-bottom, a protection level of 1.9 µA/cm2 would be
The possibility of an insulating cache under the electrolyte where the CP sys-
tem was installed was also modeled for the anode arrangements that can work
with it. The model could account for the soil surface or tank bottom and an op-
The results are presented in a series of cases based upon the three possible
Smyrl and Newman provide a criterion for determining the maximum size
∆Φo κ
ro = (12-18)
0.363iavg
where ∆Φo is the maximum over-polarization of the steel from the point of min-
takes place at -1200 mV CSE, then ∆Φo is equal to 350 mV. For an average current
density of 1.9 µA/cm2 and 10,000 Ωcm soil, the maximum radius is 51 cm. If the
174
soil resistivity is only 1000 Ωcm, then a tank of 5.1 m could be protected. A soil
The current distribution on the tank bottom is not necessarily uniform but is
of the tank should have a much higher current density than the center if Tafel ki-
netics apply (hydrogen evolution).122 If the center of the tank is poorly polarized,
the net current density will be below the mass-transfer-limited current density for
Therefore, it is not expected that equation (12-18) would be of much value for CP
design. The polarization curve used in the model for bare steel can account for
The cathodic protection system was modeled using a deep well remote ground-
bed at 250 feet below the tank. The groundbed output was 14.2 Amps. It was
assumed that the soil resistivity was constant throughout the soil domain and
equal to 10,000 Ωcm. Two cases were calculated, one for the uniform oxygen con-
centration when the tank is filled, and one with an oxygen concentration profile
centration of oxygen in the soil to simulate the conditions that exist when the
tank has just been filled and the oxygen beneath the tank is not yet depleted. The
anode output was adjusted to 14.2 Amps. Equation (12-18) would predict the dif-
ference in potential between the outer edge and the center of the tank to be -10.5
volts if there is a uniform current distribution. Figure 12-3 shows the off-potential
175
distribution on the tank bottom. The center of the tank is at -605 mV Cu/CuSO4
-100 mV polarization shift criteria.124 The outside edge of the tank is polarized to
-1189 mV Cu/CuSO4 , which is significantly more positive than the -10.5 volts pre-
dicted by equation (12-18). The potential variation is smaller because the current
density distribution was nonuniform, due to the role of corrosion and hydrogen-
as a function of the radial position from the center of the tank. The current density
176
-500 100
-700
10
-800
i/iavg
-900
1
-1000
-1100
-1200 0.1
0 0.2 0.4 0.6 0.8 1
r/r0
Figure 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 distribution would be.
(see equation (7-1) on page 104). The value of J for the secondary distribution,
described by Newman,122 was equal to 536 (see equation (7-2) in section 7.1.2).
Thus, any small change in the current density at the center of the tank used to
increase the protection there, will result in significantly elevated current densities
at the edge of the tank. The more nonuniform current density distribution results
distribution under the tank was implemented by changing the value of the mass-
the solution to the diffusion equation for oxygen given in equation (12-17). The
profile for the oxygen distribution under the tank was taken using the scaling pa-
p
rameter, R k/ D, set to 2. The anode size and location was not changed for this
177
simulation. The anode output was 14.2 Amps for the same potential setting as the
first simulation. The off-potential distribution is shown in Figure 12-5. This time,
the off-potential was much closer to being uniform. The polarization level at the
center was also significantly improved to -906 mV Cu/CuSO4 . The corrosion po-
tential at the center, where it is most negative, was -550 mV Cu/CuSO4 . Therefore
the center was polarized by over 350 mV. As one proceeds radially from the cen-
ter, the polarization level only increases. These results mean that the entire tank
was over-protected and money has been wasted by running the anode at too high
of a level.
178
-500 100
-700
10
-800
i/iavg
-900
1
-1000
-1100
-1200 0.1
0 0.2 0.4 0.6 0.8 1
r/r0
Figure 12-6: Radial potential and current density distribution when thepoxygen
concentration follows the profile given in Figure 12-2 corresponding 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.
Figure 12-6 shows both the potential and current distribution for the nonuni-
the primary distribution and the total current was the same as before. The dif-
ference was in the potential distribution, especially at the center of the tank. A
substantial reduction in current could be safely made while still maintaining pro-
the current level would have to be increased again until the steady-state oxygen
where an arbitrary number of anodes are arranged in a circle around the outside
edge of the tank. The anodes are close to the surface so it may be possible to use
179
the following examples, 8 anodes were placed around the circumference of the
the same 100 ft diameter tank used in the previous 2 examples. The total output
of the anodes was set to match closely that of the previous two examples (14.2
Amps). The anodes were placed such that they were two feet from the edge of
Uniform Oxygen Distribution Figure 12-7 shows the potential distribution for
the tank when the oxygen concentration distribution is uniform. The center of the
tank was polarized only -16 mV from the corrosion potential as opposed to the
remote groundbed where the center was polarized by -75 mV at the same total
current. The corrosion current density was 0.57 µA/cm2 which is about half that
A plot of the current density is shown in Figure 12-8. Since the distributions
were no longer axisymmetric, the position that is selected for a plot of the distri-
butions is important. The distribution is periodic however, with one period being
the angular distance between two anodes. There are two possibilities for the plot,
select a line from the center to a point closest to any one anode, or select a line
from the center to a point between two anodes. For either case, where the pri-
mary distribution shows a value of i/iavg equal to 0.5 at the center of the tank, the
value for the distributed anodes is 0.25. Near the outside edge, there is a much
sharper turn in the distribution as it goes toward infinity even though the value
of J is the same as for the case with the remote groundbed (J = 536). The corro-
sion current at the center of the tank is more than an order of magnitude higher
180
Figure 12-7: Potential distribution with distributed anodes and uniform oxygen
concentration profile. The distribution is must worse than the primary distribu-
tion. The center is only polarized -16 mV
-500 10
Off-Potential (mV Cu/CuSO4)
-600
-700
-800
i/iavg
1
-900
-1000
-1100
-1200 0.1
0 0.2 0.4 0.6 0.8 1
r/ro
Figure 12-8: Current density distribution on tank bottom with distributed anodes
and uniform oxygen concentration profile. The distribution was taken on a line
that went from the center of the tank to a point on the edge exactly between two
anodes.
181
than that of the remote groundbed. In this case icorr,Fe is equal to 0.57 µA/cm2
compared to 0.0275 µA/cm2 for the remote ground bed protected tank.
The reason for the poor performance of the distributed anodes is geometry.
The edge of the tank is significantly closer to the anodes than the center. There-
fore, is is much easier to drive current to the edge of the tank than to the center.
When a remote anode is used, all points on the tank are equally distant from the
anode.
the nonuniform oxygen concentration distribution, the potential and current dis-
tributions obtained showed better protection at the center of the tank with the
was total current driven to the steel. The corrosion current is 3 orders of magni-
tude lower than when the oxygen concentration was uniformly high. The poten-
The potential and current density distribution are plotted in Figure 12-10. It is
interesting to note that the current density at the center of the tank was the same
as that for the uniformly high oxygen profile. The additional current needed to
reduce the extra oxygen for the case with the high oxygen profile comes from the
corrosion of the steel tank. Because the oxygen concentration was lower in the
present case, the same current density was better able to protect the metal.
beds cannot be used. Distributed anodes can be used with an underlying mem-
brane because they can be placed within the soil between the membrane and the
182
Figure 12-9: Potential distribution with distributed anodes and nonuniform oxy-
gen concentration profile. The distribution is much worse than the primary dis-
tribution.
-600 10
Off-Potential (mV Cu/CuSO4)
-700
-800
i/iavg
-900 1
-1000
-1100
-1200 0.1
0 0.2 0.4 0.6 0.8 1
r/ro
Figure 12-10: Potential and current density distribution on a tank bottom with
distributed anodes and nonuniform oxygen concentration profile. The distribu-
tion was taken on a line went from the center of the tank to a point on the edge
exactly between two anodes.
183
tank. As might be expected the situation is worse than if the insulating layer did
not exist. The reason to expect poorer performance is that there are significantly
fewer current paths from the anode to the center of the tank due to the narrow
layer of soil in which the current must flow. There should also then be an increase
Figure 12-11 shows the potential distribution for a tank protected in such a
manner. The center of the tank was polarized only -13 mV from the corrosion
potential which is about 25% less than the case without the insulating cache basin.
Therefore, intuition has served well. The current density at the center of the tank
the oxygen at the center of the tank has been consumed, the polarization im-
proves significantly, but not to a required level. The center of the tank was only
polarized by -80 mV. The remote groundbed with no cache basin was able to po-
larize the center of the tank by -350 mV. The potential distribution is shown in
Figure 12-12. Again, the current density distribution was nearly identical for the
very close to the tank-bottom. The ribbon configuration can yield a more uniform
current distribution. Two anode arrangements were used for the evaluation of
185
Figure 12-13: Tank bottom with 5 and 11 parallel impressed current anodes. The
gray grid lines are on 10 ft intervals. The tank with 5 anodes uses 4 inch diameter
anodes and the one with eleven uses 1 inch diameter ribbons. The mesh for the
tank with 11 ribbons is more refined.
ribbons. The first used 5 parallel impressed current ribbons of 4 inch diameter,
the second used 11 1 inch anodes in a similar configuration. The setup for the
The mesh for the tank that had 11 anodes was made to be finer than the one
which had 5 anodes. A finer mesh was used to provide a reasonable level of
detail with the more closely spaced anodes, since there would not even be an axi-
periodic distribution seen with the distributed anodes. The anodes were assumed
to be cylindrical and were discretized in both the axial and angular directions.
The angular discretization could then account for the nonuniform current and
Figure 12-14: Potential distribution for tank bottoms with 5 and 11 parallel im-
pressed current anodes with secondary containment. The gray grid lines are on
10 ft intervals. Peak corrosion rate for tank with 5 anodes is an order of magni-
tude higher than that of the tank with 11 anodes. Total current driven to tank is
14.2 Amps for both cases.
Uniform Oxygen Distribution with Insulating Cache Basin For both anode
design cases, the total current output to the tank bottom was adjusted to 14.2
Amps. The potential distributions for the two systems are shown in Figure 12-
14. The peak corrosion current for the tank with 5 anodes was 1.6 µA/cm2 . This
large value, higher than could be expected from oxygen reduction alone, was due
to galvanic coupling between the highly polarized metal that is very close to an
anode and the poorly polarized metal between the anodes. Poorly polarized re-
gions are shown in red in Figure 12.14(a) and highly polarized regions are colored
in blue. The same color scheme and scale are used in Figure 12.14(b). Here the
potential distribution was much better, but not yet optimal, even though the an-
ode surface area was half that of the example with 5 anodes. The anode-to-anode
187
distance was cut in half with the 11 anode system. The result is that 25% smaller
driving force was needed to drive the same amount of current to the tank bottom,
Uniform Oxygen Distribution without Insulating Cache Basin When the in-
sulating cache basin was removed, some unexpected phenomena were observed
while adjusting the impressed current setting in order to drive 14.2 Amps to the
tank surface. It was estimated that roughly the same driving force would be
required to achieve the same current levels as was used for the cases using the
insulating cache. It was found, instead, that a substantial increase in driving force
was needed to obtain the same total current. To understand why there is such a
large disparity between the two systems, even though the anodes are identical, it
helps to look at the current density distribution along the length of the anodes.
The axial current density distribution for the long center anode for both the 5 and
11 anode cases is shown in Figure 12-15. In both cases, when there is no electrical
insulator under the the anodes, most of the current comes from either end. The
The current density and corrosion current distribution for the worst case of
the 5 ribbons is shown in Figure 12-16 The largest corrosion current is next to the
long center ribbon at either end. There is evidence of galvanic coupling about
two-thirds of the distance from the center of the tank to the edge. The corro-
sion current density is about 2 µA/cm2 which is double the rate due to general
corrosion.
188
10000
10000
Without Insulating Cache Basin
µA/cm )
Without Insulating Cache Basin
2
µA/cm )
2
Current Density (µ
1000
Current Density (µ
100
100
10
10
0 20 40 60 80 100 0 20 40 60 80 100
Position (ft) Position (ft)
Figure 12-15: Axial current density distributions on anodes for both the 5 and
11 anode cases. For each case, the current density distribution is plotted for the
situations of having an insulating cache basin and only having conductive soil
under the tank.
Figure 12-16: Potential and corrosion current distribution for tank bottoms with
5 parallel impressed current anodes with no secondary containment. The gray
grid lines are on 10 ft intervals. Total current driven to tank is 14.2 Amps. There
are significantly more anodic sections in this case with 5 anodes than for the case
with an insulating cache basin.
189
Figure 12-17: Potential distribution for tank bottoms with 5 and 11 parallel im-
pressed current anodes and secondary containment with nonuniform oxygen dis-
tribution. The gray grid lines are on 10 ft intervals. Total current driven to tank
is 14.2 Amps for both cases. Center of tank is more polarized than in the case
shown in Figure 12-14
a nonuniform oxygen concentration profile was modeled for both the 5 and 11
anode CP systems. The anodes were again adjusted such that the total output
was 14.2 Amps. The potential distributions are shown in Figure 12-17. As ex-
pected, the center of the tank was polarized to a more negative potential then
that shown in Figure 12-14. The current distribution was nearly identical to that
case involved removing the restriction of the underlying insulating cache. The
190
Figure 12-18: Potential distribution for tank bottoms with 5 and 11 parallel im-
pressed current anodes and no secondary containment with nonuniform oxygen
distribution. The gray grid lines are on 10 ft intervals. Total current driven to
tank is 14.2 Amps for both cases.
potential distribution is shown in Figure 12-18. In this case, the center of the tank
was not as well polarized as in the case where the cache basin is present. Again,
the anodes. The current density was much lower at the centers of the anodes than
it was with the cache membrane present. To compensate, the current density was
higher at the ends of the anodes which, then, promoted a larger current density
at the edge of the tank. The current was better distributed than when remote or
Use of larger numbers of ribbons also served to improve the current distribu-
tion for all cases where ribbons were used. If the spacing was too large, there was
a possibility of galvanic coupling between the sections of the tank that were more
191
polarized directly above the anodes and sections in between the anodes. For the
100-ft diameter tank bottom used in the examples, it appears that a 10-ft spacing
12.5 Conclusions
Even under the best of situations, it is very difficult to provide cathodic pro-
around the circumference of the tank are used for protection, the best that can be
hoped for is a current distribution that approaches the primary as given in equa-
tion (7-1). In the case of distributed anodes around the perimeter, the situation
is much worse than for the remote ground-bed. Because of the proximity of the
anodes to the outside edge of the tank, most of the current goes to polarizing the
The only solution that can work is to place anodes directly under very close
to the steel they are to protect. Design of such systems is required. If the spacing
between anodes is too great, anodic regions can form on the steel between the
anodes. The goal should also be to provide the largest amount of anode surface
area possible to reduce to driving force needed to supply a sufficient total current
to the tank. With a smaller driving force, there is significantly less likelihood of
Cyclic oxygen concentration gradients along the tank bottom due to periodic
emptying and filling of the tank can cause periodic over- and under-protection.
The location of the reference electrode is very dependent on the design of the
around the circumference of the tank, is is most important to monitor the center
of the tank. Other design require modeling of the design to determine optimal
locations. As shown in Figure 12-14 the optimal location is not at the center of
the tank but near the longest ribbon about 3/4 of the distance from the center to
When a nonuniform oxygen distribution occurs, the parallel ribbon type an-
ode system cannot be adjusted at the center of the tank. A high level of polariza-
tion occurs at the center because of the lack of controlability. A new design for
the anode layout has been proposed by Koszewski.118 The arrangement consists
In the end, it was quite by accident that the best means of properly protecting a
tank bottom was found. By restricting the current paths, the lining plays an active
role in cathodic protection design and forces a more uniform current distribution
that all new tanks have secondary containment was designed to contain spills,
in Chapter 2, was that the ionic strength of the electrolyte within the domain
was constant. Many times, however, it may be desirable to relax that restriction.
Greens function solutions for these types of problems are very difficult to ob-
tain. There have been some solutions reported61 for the case when the conductiv-
ity is a linearly varying function in only one of the three dimensions. However,
side of the region with the linear variation. If one does not do this, then one end
of the region will have a negative conductivity which has no physical meaning
while the extreme other end will have near superconducting characteristics. If a
Green’s function cannot be used, the only alternative for the boundary element
method is to divide the domain into regions. At the point of division between
Often, such as at river crossings, the change in soil conductivity is abrupt and
can be modeled as a plane. Figure 13-1 shows how the division must be modeled
using the boundary element method. The point of the change in soil properties is
193
194
Figure 13-1: Plane meshed for a soil type change. Three pipes are shown passing
through a point where the soil conductivity changes.
195
two domains of constant conductivity. The two domains must then be coupled
together through the common boundary of the mesh at the soil type division.
If steady state is assumed, then there can be now accumulation of charge on any
either side of the division. The reason for allowing separate values of potential on
either side of the division is due to a phenomenon known as the liquid junction
tween the anions and cations. The change in potential can be calculated through
the Henderson formula if the concentrations and diffusion coefficients for all the
charged species are known. The continuity condition at the soil division becomes
Φ I = Φ I I + ΦLJ (13-2)
where ΦLJ is the liquid junction potential. Careful attention must be paid to the
13.2 Implementation
Soil divisions are implemented within the boundary element method only.
Laplace’s equation is written for each domain which is now bounded by the soil
196
surface (see section 3.5), and one or two soil type division half-planes. Then a set
of algebraic equations from equation (4-1) can be written for for the first domain
H j , j H j ,i Φ j G j , j G j ,i ~
n · ∇Φ j
= (13-3)
Hi, j Hi,i Φi Gi, j Gi,i ~n · ∇Φi
where i represents nodes and elements on the soil division interface and j is for
all other nodes and elements within the domain. The second domain would be
written as
Hi,i Hi,k Φi Gi,i Gi,k ~n · ∇Φi
= (13-4)
Hk,i Hk,k Φk Gk,i Gk,k ~n · ∇Φk
where i again represents nodes and elements on the same boundary as in equa-
tion (13-3) and k represents all nodes and elements not on the domain boundary
and within domain two. Applying the continuity conditions from equations (13-
where the subscripts 1 and 2 refer the the domain the matrix elements come from
and ΦLJ,1−2 comes from the liquid junction potential between domains 1 and 2
calculated using the Henderson formula. If there are n domains in the problem
197
then n − 1 extra columns are added to the right-hand matrix to account for liquid
junction potentials.
The equations for the conduction of current through the pipe steel, bonds and
anodes is solved exactly the same as described in section 3.7. The only caveat
is that there are no entries in the finite element matrices for the soil domain divi-
sions since they are not part of the steel domain where the FEM is applied. Again,
Many domains can be strung together using the above method. Figure 13-
2 shows how the mesh for a system with two soil divisions and therefore three
different soil types. This model setup may be appropriate for a situation where
198
Figure 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.
two pipes are passing through a river bottom where the saturated section of soil
as in equation (13-6).
4.5 and 4.9, following the procedure outlined in Figure 4-1. Sparse techniques can
199
Sparse techniques cannot usually be used with mixed boundary condition prob-
lems under the formulation shown in equation (13-7). After all of the unknown
boundary conditions are moved to the left-hand-side, many elements in the the 0
submatrices above the K submatrices are filled with G submatrices. When non-
linear boundary conditions are used, many of the remaining terms with zeros in
the upper triangle above the K submatrices are filled with terms from the partial
14.1 Conclusions
Buried pipelines carry most of the United States’ natural gas and petroleum
products from wells to refineries and to final consumption points. All of the
pipes are subject to corrosion and therefore have cathodic protection systems to
protect them. To reduce the amount of current required to protect the pipes,
high-resistivity coatings have been applied to the pipes. Coating defects expose
portions of bare steel to the soil environment, and cathodic protection is conse-
quired because the standard design equations cannot account for the behavior
of the cathodic protection system when coating defects are present. The design
of cathodic protection systems has been further complicated as demand has in-
creased because new pipes have been placed next to old pipes in the same right-
of-ways.
today. New methods for treating the coating have been developed to predict
more accurately the protection levels of coated pipes with holidays. Complete
axial and angular current density and potential distributions can now be calcu-
200
201
Interactions between pipes with coating holidays have been shown in Chapter
pipes.
It has been shown that the model can easily be extended to geometries other
than pipes through the investigation of cathodic protection of tank bottoms. Through
The largest weakness of the model is the lack of good polarization data for
both the bare and coated steel. Polarization parameters are strong functions of
steel composition and pre-treatments, soil environment, time, and history. Ken
be modified to integrate the time evolution of the system. Because the bound-
ary conditions are time-dependent, additional derivations are needed. One at-
tempt has been made by Santiago and Telles.126, 127 They incorporated the time-
Carson has provided data for a similar formulation of the polarization curves for
The model can be easily extended to other applications. For example, ship
hulls can be treated in a similar manner as pipes. The major modification needed
is to provide a means to mesh the geometry of the hull and propellers, and to
202
provide polarization curves for the metals used in their construction. Wells could
also be treated if a technique is available to account for the multiple soil layers
the well passes though. One method may be to use the techniques described in
Chapter 13. The domain divisions normal to the soil surface would be replaced
Cathodic protection of the interior of pipes, storage tanks, reactor vessels, heat
exchangers can also be treated with some simple changes to the formulation. In
all cases, high-quality polarization data would be needed for the specific envi-
ronments.
APPENDIX A
PARTIAL DERIVATIVES OF THE BOUNDARY CONDITIONS
When solving the system of nonlinear equations that result from the bound-
ary element method, the linearization process requires partial derivatives of the
boundary conditions for each of the independent variables. Using the equation
(1-26) is really the sum of two equation eliminating i. The original two equations
are
Apore
V−Φ −E
in Fe
i = 10 βFe
A
!−1
1
V −Φin − EO
2 (
− V −Φin − EH
2 )
βO βH
− − 10 2 − 10 2 (A-4)
(1 − αblock ) ilim,O2
203
204
and an IR drop
Φ − Φin
i= (A-5)
ρfilm δfilm
Equations (A-4) and (A-5) are solved simultaneously to get the current density i
∂i
(A-6)
∂Φ
and
∂i
(A-7)
∂V
are obtained numerically.
∂i i(Φ + h, V) − i(Φ + h, V)
= (A-8)
∂Φ 2h
where h is a small number. The size of h must be carefully chosen such that
i(Φ + h, V) − i(Φ + h, V) is not equal to zero because h is smaller than the noise
due to finite machine precision and i(Φ + h, V) − i(Φ + h, V) does not straddle
large variations in the function. The value of h is selected based on the mag-
nitude of the variable it is modifying, i.e., Φ. The value of the machine , the
smallest number added to 1 such that the result is different than 1, is used to find
h according to
√
h = max(Φ, 1.0) (A-9)
ing commands are platform independent, so code that utilizes OpenGL can easily
be ported across different platforms. The only caveat is that the OpenGL device
setup is not platform independent. The reason is that windows are created in dif-
NT,‡ Iris,§ Solaris,¶ AIX,k etc. Therefore, a small subset of code written by the
The OpenGL device on the Windows NT platform must be created and at-
tached to a window before any drawing commands are issued. The device is
created during the window creation subroutine. As soon as the operating system
has created a generic window, a Pixelformat is specified for the OpenGL device
based of the capabilities of the window. Once the the Pixelformat has been suc-
cessfully chosen, the rendering device is created and attached to the window
205
206
text may be active at any given time, a command is given to the new context to
make it the current one. The command must be issued every time a drawing cy-
cle begins. Source listing B.1 contains the code used to set up the OpenGL device
{
36 // failed to set the pixel format
AfxMessageBox("Failed to set the pixel
38 format. Be sure that OpenGL is installed");
return -1;
40 }
60 // white background
glClearColor(1.0f,1.0f,1.0f,1.0f);
62 // color properties of materials will be used
glEnable(GL_COLOR_MATERIAL );
64
68 glEdgeFlag(TRUE);
// polygons are numbered counterclockwise
70 glFrontFace(GL_CCW);
// dont draw faces that face away from the viewer
72 glCullFace(GL_BACK);
glClearDepth( 1.0 );
74 return 0;
}
Line 4-47 are used to determine the number of colors OpenGL will use. The
The drawing cycle begins with a message from the operating system or from
the interactive code written by the programmer that says something within the
window has changed and it needs to be redrawn. The first step is to make the
windows OpenGL device the current one so that drawing commands are not
sent to the wrong window. Next all the platform independent OpenGL drawing
code is called. Finally SwapBuffers is called to place all the newly rendered
objects into the window. The code below in source listing B.2 demonstrates the
procedure.
Source Listing B.2: Draw command called by operating system whenever a win-
dow needs updating.
void CObjView::OnDraw(CDC* pDC)
2 {
// get the size of the window
4 CRect rect;
GetClientRect(&rect);
6
26 DrawRecFace(&rect);
lighting setup and texture setup. First, the perspective transform matrix is cre-
ated and model view transform matrix are created. Then polygon rendering
properties are set. Finally all the objects are asked to draw themselves into the
OpenGL device.
Source Listing B.3: Platform independent code to setup to model view and pro-
jection transformation matrices.
void CObjView::DrawRecFace(CRect *prect)
2 {
// and the view port to the window size
4 glViewport(0,0,prect->right,prect->bottom);
// clear the color and the depth
6 glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
glEnable(GL_DEPTH_TEST);
8
glMatrixMode(GL_MODELVIEW);
20
glLoadIdentity();
22 // Position / translation / scale
// move the eye back from the screen
24 glTranslated(0.0,0.0,m_dzTrans);
// rotate the drawing
26 glRotated(m_dxRotation, 1.0, 0.0, 0.0);
glRotated(m_dyRotation, 0.0, 1.0, 0.0);
28 // now change to the z=depth coordinate system
glMultMatrixd(m_dMatrixCorrd);
30 // spot to lookat (world x, y, z corrd)
glTranslated(m_dxTrans,m_dyTrans,m_dDepthTrans);
32 glScaled(m_dScaling,m_dScaling,m_dScaling);
The RenderScene command cycles through all of the objects in the model,
asking each to draw itself. It also draws reference objects such as the soil surface
An example piece of C++ code is given below in listing B.4 that draws a set
rendering.
Source Listing B.4: Low level OpenGL code to draw a wireframe rendering of a
single object.
void CBaseMesh::Render(PARAMS *pParas) {
2 int ii, jj;
// start drawing 4 node faces
4 glBegin(GL_QUADS);
// use black wires for a wire frame
6 glColor3f(0.0f,0.0f,0.0f);
for(ii=0; ii<m_nNumQuadFaces; ii++)
8 {
for(jj=0; jj<4; jj++)
10 {
glVertex3d(m_pNodes[m_pQuadFaces[ii].nPoints[jj]].x,
12 m_pNodes[m_pQuadFaces[ii].nPoints[jj]].y,
m_pNodes[m_pQuadFaces[ii].nPoints[jj]].z);
14 }
}
16 glEnd();
As can be seen in listing B.4 lines 11 and 24, the actual boundary element or
The mouse is used for much of the interaction with the model. All three
buttons and the wheel are used. The middle mouse button is used strictly for
OpenGL interactions. It used used to turn and rotate the model, zoom in and
out, translate from one position to another, and change the depth the viewer is
looking at. Much of the technique used to develop the interaction mechanisms
tions to the techniques have been made to make make them consistent with the
OpenGL API and to allow interactions specific to pipeline networks and buried
structures.
Objects are moved around within the OpenGL window using the mouse. The
sent from the operating system. There are several different messages depending
on the type of mouse action. For each message, there is a corresponding function
OnRButtonDblClk. The three functions used for changing the view of an OpenGL
are needed since the display should only move while the middle button is held
1. The user clicks and holds the middle mouse-button (sometimes a wheel)
2. The user moves the mouse while holding the middle mouse button down.
The display updates in real time as the mouse is moved. The number of
A movie of the process is given in Figure B-1. The movie shows one of the pos-
sible modes of movement known as rotation and orbiting. It also shows how
The code used to accomplish the mouse control of movement utilizes several
variables. The first is to set a flag to indicate whether or not the middle mouse
Source Listing B.5: Code to respond to the middle mouse button being depressed.
void CObjView::OnMButtonDown(UINT nFlags, CPoint point)
2 {
// set flag that the mouse is depressed
4 m_bMButtonDown=TRUE;
// keep mouse with this program even
6 // if the mouse moves outside its window
SetCapture( );
8 // save the point where the mouse is first clicked
m_ptMDownPos = point;
10 }
The function SetCapture() allows the mouse to function within the current
program even if the user drags the mouse outside the window for the program.
214
Figure B-1: A movie demonstrating how the user can use the middle mouse but-
ton 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.
215
The work of changing the appearance of the view window is done in the function
OnMouseMove.
Source Listing B.6: Code to respond to mouse movement within the OpenGL
window.
void CObjView::OnMouseMove(UINT nFlags, CPoint point)
2 {
CRect rect;
4 GetClientRect(&rect);
// check if a button is depressed
6 if(m_bMButtonDown)
{
8 CConView *pView=GetConView();
// check which type of movement has been desired
10 if(m_nMouseMode==0)
{
12 // Zoom mode
m_dzTrans-=fabs(m_dzTrans)*(m_ptMDownPos.y-
14 point.y)/(0.5*rect.bottom);
if(m_dzTrans>0.0)
16 m_dzTrans=-m_dzTrans;
}
18 else if(m_nMouseMode==1)
{
20 // Orbit and rotate mode
m_dyRotation -= (m_ptMDownPos.x - point.x)/3.0;
22 m_dxRotation -= (m_ptMDownPos.y - point.y)/3.0;
//put limits on the rotation angles
24 // or loop back to zero when 360 is exceeded
if(m_dyRotation>360.0)
26 m_dyRotation-=360.0;
if(m_dyRotation<-360.0)
28 m_dyRotation+=360.0;
if(m_dxRotation>90.0)
30 m_dxRotation=90.0;
if(m_dxRotation<-90.0)
32 m_dxRotation=-90.0;
}
34 else if(m_nMouseMode==2)
{
36 // translate across the soil surface
// going to use both bottom and right member
38 // first use x movement on screen to left -right
// movement of the 3D soil surface
216
40 m_dxTrans+=((fabs(m_dzTrans)+1.)*sin((m_dyRotation)*
PI/180.0)*(m_ptMDownPos.x - point.x)/rect.right);
42 m_dyTrans-=((fabs(m_dzTrans)+1.)*cos((m_dyRotation)*
PI/180.0)*(m_ptMDownPos.x - point.x)/rect.bottom);
44
In the OnMouseMove function, the first thing done is a check for which mouse
button has been depressed. For the case of the middle mouse button, there
are 4 different modes of movement which are specified through the variable
variable depending on the movement mode. Looking at the code for the OpenGL
perspective and model view matrix transformations, it can be quickly seen how
217
modifying each of the variables with the mouse subroutines translates into some
Source Listing B.7: Code to change view of objects due to mouse control.
void CObjView::DrawRecFace(CRect *prect) {
2 ...
// set the projection view
4 glMatrixMode(GL_PROJECTION);
glLoadIdentity();
6 // set clip to one tenth the eye distance and 5 times
// the largest extent
8 gluPerspective(45.0,aspect,fabs(m_dzTrans)*0.1,5.0*
pDoc->GetDataSize());
10
There are only two rotation modes used (lines 19 and 20 in listing B.7) be-
cause a third is redundant. Line 22 of the same listing is very important because
OpenGL uses a left-handed coordinate system. All of the data is created using a
system change.
218
Context menus are displayed within the OpenGL window when a using right-
clicks on an object. The type of object determines the menu displayed. Since a
3D view can have several objects appear at the same place because one is directly
in back of another, there has to be a mechanism to select the closest object. That
method is called picking and is described in detail in the OpenGL 1.2 Program-
ming Guide.106 It involves rendering the seen through a special buffer to find the
closest object to the viewer that is at the point the person clicked on. Once an
The picking code first creates a memory buffer to render into to find the object
being selected (lines 4 and 19 of listing B.8) and tells the OpenGL engine to use the
selection rendermode (line 22). It then sets up the OpenGL model and projection
matrices to be identical to the ones used to display the objects (lines 28-37). After
the objects have been rendered into the special buffer, the number of objects hit
is returned (line 45). Then the code finds the hit with the smallest z depth (not
depth in the soil but distance from the view straight into the computer screen).
Source Listing B.8: Code to find what object was clicked over.
int CObjView::ProcessMouseHit(int xPos, int yPos, BOOL bIsMenu) {
2 CObjDoc* pDoc = GetDocument();
GLuint *buffer;
4 buffer=new GLuint[512];
// save the sky variable (don’t want
6 // hits to the sky faces)
BOOL bUseSky=this->DisplaySky();
8 m_bDisplaySky=FALSE;
GLint hits;
10 // clear the error buffer
int nError=glGetError();
219
12
glSelectBuffer(512,buffer);
20 // tell the OpenGL device we are checking for
// mouse clicks on a face
22 glRenderMode(GL_SELECT);
// Set up the OpenGL Device
24 glInitNames();
glPushName(0);
26 glEnable(GL_CULL_FACE);
glEnable( GL_DEPTH_TEST );
28 glMatrixMode(GL_PROJECTION);
30 glLoadIdentity();
gluPickMatrix((GLdouble)(xPos), (GLdouble)(m_viewport[3] -yPos),
32 1, 1, m_viewport);
34 glMultMatrixd(m_projMatrix);
glMatrixMode(GL_MODELVIEW);
36 glPushMatrix();
glLoadMatrixd(m_modelMatrix);
38
// Start rendering...
40 pDoc->RenderScene(this, 1);
glPopMatrix();
42 glFlush();
96 delete buffer;
// tell everyone that no faces were hit
98 return -1;
}
221
vide the user with data about the object under the point clicked. The mouse then
can provide any type of data calculated on any surface. The soil does not re-
strict it. The model currently outputs potential, net current density, off-potential,
is placed.
APPENDIX C
CODE STRUCTURE
The source code for the model was written as a Microsoft Windows NT ap-
plication utilizing the OpenGL Application Programming Interface (API) for in-
teractive display of the pipeline networks. The code uses a simple C++ structure
to help manage code development. It may be best to explain the structure from
the interface down through the data structures and then into the calculation en-
gine. The code is divided into three major components: i) the interface ii) the
C.1 Interface
chitecture. Each of the two parts is implemented as a C++ class. A C++ class is
a piece of self contained code that can perform a specific task, manage a data set,
C++ class is far beyond the scope of this work. The reader is referred to some
of the classic texts for C++ programming.38, 129 Bruce Eckel has a very good ex-
planation of C++ objects and classes in Chapter 1 of his book “Thinking in C++,
Volume 1”.39
class. The division of responsibilities is made under the premise that one class
222
223
Figure C-1: 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.
manages data (the Document) while the other class controls how the data appears
The model adds a second View type class because of the nature of the data
and its display. Figure C-1 shows the main window as it is divided into two
section. The left half is a control window meaning it contains buttons, edit boxes
and other controls. The right half contains the OpenGL 3D display of the pipes.
Each half requires its own class to manage user input a feedback. Both classes are
connected to the same data, so there is only one class for data management.
The names for the classes follow their purpose. The control console is called
CConView, which is short for Class Console View Window. The OpenGL
window class is call CObjView which is short for Class Object View Window.
The document class is called CObjDoc, which is short for Class Object Doc-
ument.
224
There are also several support classes needed. The model is setup to allow
of all of the open models within the main program. The main program is built
as a class. Its purpose is to establish how all of the classes are linked together to
form a single coherent model. The last support class is the main frame window.
This widow contains the standard menu, toll button bars and status bars. It takes
messages when a user clicks a tool bar button or selects a menu item, passes
the message to the document manager which then forwards the message to the
top level model that is open within the main program. Details on this type of
structure may be found in a number books130, 131 and in the documentation for
The document class is used to manage all aspects of the the data. It owns
all of the major data structures such as the boundary element and finite element
meshes, boundary condition database and current and potential distribution so-
lutions. It handles such tasks as saving and restoring the data to disk (persis-
tence), creating new data such as pipes (allocation) and destruction of data (clos-
It takes a single argument of an CArchive object. The object contains all the
methods needed to read or write data to a disk. Therefore, data owned by the
class is written by calling functions from the CArchive object with the data as
225
parameters. However, there are very few variables within the Document class
that are saved that way. Instead, the document Serialize function asks the
data structures it owns to save themselves passing the CArchive object to the
Source Listing C.1: Code to save and restore a model to persistant storage (disk).
void CObjDoc::Serialize(CArchive& ar)
2 {
if (ar.IsStoring())
4 {
// sentinel (file version and format)
6 ar<<m_dFileVersion;
}
8 else ... //read the file
{
10 // sentinel (file version and format)
ar>>m_dFileVersion;
12 // check version
if(!CheckVersion())
14 {
throw(CFileException());
16 return;
}
18 }
// save all the pipes
20 m_pPipeList->Serialize(ar);
// color scales
22 m_pclPotRange->Serialize(ar);
// Boundary Condition database
24 m_pBCDataBase->Serialize(ar);
// render data
26 ((CChildFrame *)GetParentFrame())->Serialize(ar);
// calculation engine setup
28 m_pBemSpec->Serialize(ar);
// soilzone database
30 m_pSoilZoneDB->Serialize(ar);
// CP zone list
32 m_pCPZoneNameList->Serialize(ar);
// bond list
34 m_pBondList->Serialize(ar);
226
// soil-surface sections
36 m_pSoilSurfaceSectionList->Serialize(ar);
}
Of course, everything saved is binary data. The data structures are so large
and complex, there is no point to using a human readable format. Also, a binary
The document also passes commands from the View to appropriate data struc-
tures, or creates new structures as needed. For instance, the view will want all
of the pipes to draw themselves into the client area of the window. It does this
by asking the document. A function is called in the document by the view where
the document then calls the appropriate data structures to draw themselves.
Source Listing C.2: Render scene code from the document. Code asks data struc-
tures to draw themselves.
void CObjDoc::RenderScene(CObjView *pView, int nMode) {
2
int ii;
4 int nStart=0;
if(pView->DisplaySky())
6 {
// draw the sky
8 GLint polymode[2];
glGetIntegerv(GL_POLYGON_MODE, polymode);
10 glPolygonMode(GL_FRONT , GL_FILL);
m_pSky->Render();
12 glPolygonMode(GL_FRONT , polymode[0]);
}
14 if(pView->RenderAxis()&& nMode==0)
m_pAxis->Render();
16
if(pView->ShowBonds())
26 {
// draw the bonds
28 pos=m_pBondList.GetHeadPosition();
while(pos!=NULL)
30 m_pBondList.GetNext(pos)->RenderBond(nMode);
}
32 // draw the soil surface
pos=m_pSoilSurfaceSection.GetHeadPosition();
34 while(pos!=NULL)
m_pSoilSurfaceSection.GetAt(pos)->Render(pView);
36 }
The Document provides support for to the code that starts the calculation
The view class controls how the user interacts with the data. It is responsible
for everything that happens or appears with the client area of the main window.
The client area is the part of the window that is not the title bar, borders or resiz-
ing handles. It is also responsible for most of the menu items and toolbar buttons.
The view class must display dialog boxes as they are needed in response to
user actions. It also handles how the objects appear to the user. The techniques
used to display and interact with the objects in the model is described in Ap-
pendix B.
The second view class is used to support the console window (see Figure C-
2). It provides code to respond to all of the controls in the console window and
provide the feedback information. The upper half of the window is for the cal-
culation engine, the lower half is for the OpenGL window and virtual reference
electrode.
228
Figure C-2: Console window. This window provides model control, calculation
control and feedback.
229
the user clicking the “RunModel” button. The process begins by checking for
necessary minimums in the model such as a pipe and at least one anode. It then
checks to see if the model has ever been saved. If it has not yet been saved,
the user must save it before the model can be run. The code then creates a data
structure object to pass to the calculation engine. The object contains pointers
to the pipe and anode data, the boundary condition database, and a pointer to
functions needed to supply the user with where to model is in the progress of the
calculation. After that, the file that contains the calculation engine is loaded into
has been successfully loaded and is the correct version for the machine’s CPU,
again far beyond the scope of this dissertation, so the reader is referred to some
Source Listing C.3: Code to spin off a thread to run the calculation. Two functions
are needed.
void CConView::OnConRunModel() {
2 CObjView *pView=GetRenderView();
CObjDoc* pDoc = pView->GetDocument();
4 // make sure the CP systems have been organized
pDoc->OrganizeCPSystems();
6 // make sure there are some objects
if(!pDoc->CheckCalcMins())
8 {
return;
10 }
// force a save if the document does not have a path yet
12 if((pDoc->GetPathName( )).GetLength()==0 )
{
14 // force a save
SendMessage(WM_COMMAND, ID_FILE_SAVE, 0);
230
16 }
m_ctrlNRError.SetWindowText("not set");
18 m_ctrlFEMError.SetWindowText("not set");
// start the elapsed timer
20 m_unStartTick=GetTickCount();
// set auto update interval to 0.5 seconds
22 SetTimer(0x00000001, 500, NULL);
m_ctrlRunModel.EnableWindow(FALSE);
24 m_ctrlPauseRun.EnableWindow(TRUE);
m_ctrlCancelRun.EnableWindow(TRUE);
26 // fill in the nodes and elements in the bemspec
pDoc->FillBemSpecForCalc();
28 BEMSTARTSTRUCT *pbstart;
pbstart=new BEMSTARTSTRUCT[1];
30 pbstart->pSpecs=pDoc->GetBEMSpec();
pbstart->pSpecs->SetParenthWnd(m_hWnd);
32 pbstart->pSpecs->SetBemStopFlag(FALSE);
pbstart->pSpecs->SetBemPauseFlag(FALSE);
34 pbstart->pSpecs->SetBondList(&(pDoc->m_pBondList));
pbstart->pSpecs->SetSoilSurfaceList(
36 pDoc->GetSoilSurfaceSecList());
pbstart->pWnd=m_hWnd;
38
54 ResetControls();
delete pbstart;
56 return;
}
58 // check the cpu
CPUCHECK checkcpu=(CPUCHECK)GetProcAddress(m_hInstEng,
60 "fnIsCPUCorrect");
231
if(checkcpu==NULL)
62 {
AfxMessageBox("Version of calculation engine wrong,
64 or cannot determine the cpu in this machine.
Contact your software vendor.");
66 delete pbstart;
ResetControls();
68 return;
}
70 if(checkcpu()==FALSE)
{
72 AfxMessageBox("The calculation engine installed will
not work with the cpu on this machine.
74 Contact your software vendor.");
delete pbstart;
76 ResetControls();
return;
78 }
pbstart->m_hInstEng=m_hInstEng;
80 if(m_hInstEng!=NULL)
{
82 AfxBeginThread(MyControllingFunction, pbstart,
THREAD_PRIORITY_BELOW_NORMAL);
84 pDoc->CanEdit(FALSE);
}
86
}
88
106
// cleanup
108 ((BEMSTARTSTRUCT *)pParam)->pSpecs->DeleteElementPointers();
The calculation engine is divided into several sections. The first is the calcula-
tion of the coefficient matrices corresponding the boundary element method and
the finite element method. The matrices are kept in separate files because of the
sparse nature of the coupled system (see section 4.8). The actual integration is
performed by a separate section of code described in section C.3. After the ma-
trices have been assembled, the boundary conditions are applied and a solution
obtained. With the solution in hand, post-solution calculations are done. Cur-
rently, only on-potential distributions on the soil surface within user specified
The coefficients of the boundary element matrices (equation (4-2) are the in-
tegrals given in equation (4-1). The two types of integrals are placed into their
respective matrices as each set of integrals is completed. A set consists of the in-
tegrals resulting from a collocation point - element pair. If the element has nine
nodes, then there are 18 integrals in the set, half of which belong to the H matrix
233
of equation (4-2), and the other half go into the G matrix. For the six node tri-
angles, 12 integrals are calculated with the same division between matrices. The
data is aligned within the G matrix such that all of the results for a given element
are next to each other if the storage format is row major. Therefore, the transpose
of the matrix must be taken after all of the integrals are completed. The results of
the integration belonging to the G matrix are fed directly to the harddrive as each
integral set is completed. The results for the H matrix are retained in memory.
The final format of the files on disk must be column major. This has significant
the boundary conditions and solution of the set of nonlinear equations) as dis-
cussed in section 9.7. Column major means that columns are written sequentially
to disk instead of rows. Row major is the default method for C/C++ for both
matrices in memory and on disk. Fortran always uses column major and owes
The integration makes use of remote processing which means that many com-
puters can contribute to the completion of the integrals in parallel. Section C.3
calculated on the local machine use the same integration technique as used for
The finite element matrices are assembled directly by the calculation engine
without the use of any remote processing. The code has been written for parallel
computation on the local machine so all processors are utilized. The reason for
234
local computation only is that the finite element calculation is essentially a two
The code is written as a C++ class that either reads the data from disk, or per-
forms the calculation. The outer loop of the assembly code is over the elements
and calls a function to create the element stiffness matrix. The outermost loop of
the element stiffness matrix code is over the Gauss points to greatly reduce the
overhead of calculating the Jacobian matrix. The inner loop is over each of the
entries in the element stiffness matrix. Since the element stiffness matrix is sym-
metric, only the upper triangle needs to be calculated. However, the complete
matrix must be used when integrated with the BEM because the BEM is not sym-
metric. Even if it was, the solution technique for the nonlinear set of equations
Once the boundary element and finite element matrices are complete, the
boundary conditions can be applied and the system of nonlinear equations can
be solved. It starts by reordering the variables in equation (4-39) such that all of
the unknowns are on the left-hand-side. The result is equation (4-41). The so-
lution continues with an initial guess of the unknown boundary conditions and
then a calculation of the known boundary conditions from the guessed condi-
tions. The first residual vector is then calculated and a sum of squares error is
made to estimate how far off the initial guess is. If the guess is not close enough
minimize the norm of the residual. The code for the line search is given in code
listing C.5.
44 Info=0;
char trans=’N’;
46 int strlength=1;
/*****************************************************/
48 // solve Ax=b
#ifdef _DEBUG
50 int rVal;
Vector<Subscript> ipiv(nLim);
52 rVal=LU_factor(m_mA,ipiv, (UPDATEF)UPDATE);
if(rVal==-1)
54 {
AfxMessageBox("Singular Matrix, Cannot solve.");
56 m_pBemSpecs->SetBemStopFlag();
m_mA.newsize(0,0);
58 return;
}
60 LU_solve(m_mA,ipiv,p);
#else
62 Vector<Fortran_integer> index(nLim);
DGESV(&N, &NRHS, &m_mA(1,1), &LDA, &index(1),
64 &p(1), &LDB, &Info, (void *)UPDATE);
if(Info>0)
66 {
AfxMessageBox("Singular Matrix, exiting engine.");
68 m_pBemSpecs->SetBemStopFlag();
m_mA.newsize(0,0);
70 return;
}
72 pBEM->UpdateStatus(89, 4);
#endif
238
74
/*****************************************************/
76 // update the user interface
PostMessage(m_pBemSpecs->GetParentWnd(), WM_COMMAND,
78 32822, 10);
// do the line search
80 // lnsrch tries to find lambda <=1.0 such
// that the norm of the residual is reduced
82 // lambda is the fraction of the full Newton step
lnsrch(pBEM, nLim, m_vR, xold, fold, g, p, m_vX,
84 &f, stpmax, &Info);
// get NRerror
98 if(pBEM->UseAttenuation())
{
100 int ONE=1;
m_dNRError=cblas_dnrm2(m_nNumNodes, &m_vR(1), ONE);
102 m_dFEMError=cblas_dnrm2(NPZ, &m_vR(m_nNumNodes+1), ONE);
PostMessage(m_pBemSpecs->GetParentWnd(), WM_COMMAND,
104 32802, (LPARAM)(&m_dNRError));
PostMessage(m_pBemSpecs->GetParentWnd(), WM_COMMAND,
106 32841, (LPARAM)(&m_dFEMError));
}
108 else
{
110 int ONE=1;
m_dNRError=cblas_dnrm2(nLim, &m_vR(1), ONE);
112 PostMessage(m_pBemSpecs->GetParentWnd(), WM_COMMAND,
32802, (LPARAM)(&m_dNRError));
114 }
if(test<xTol)
116 break; // break the while loop
118 test=0.0;
239
Soil surface on and off-potentials are calculated using equation (3-8). Again,
at the soil surface. This time, the source points are not on a boundary, so there
are no singular integrals. The results of the integration are again placed in a pair
The remote integration module consists of two parts: a client and a server.
The server is a piece of code the actually performs the integration. The client
between the client and server is handled through Remote Procedure Calls (RPC).
RPC functions by taking and ordinary function call and passing it over a network
function is actually executed. The results of the function are then passed back to
data structures that are to be passed to the function in terms of basic data types
such as double, int, or char. It also contains a globally unique identifier to insure
that no-one anywhere else in the world has made a definition the same as this
one. Finally the code contains the definitions of the functions that are to be called.
Each parameter passed to the function contains a label as to whether the param-
eter is an input, output or both. This is done to reduce the amount of data that
must be transported across the network when the function is called and when it
18
24 long nUseHighSpeed;
}INTPARAMETERS;
26
DOMAINRESULTS IntegrateDomainPoints(
42 [in] handle_t hBind,
[in] INTDATA data);
44
RESULTS IntegrateElement(
46 [in] handle_t Binding,
[in] INTDATA theData);
48
long GetNumProcessors(
50 [in] handle_t Binding);
};
There is an extra function added to tell the client program how many proces-
sors are available to do calculations. Then the client program can create as many
CORBA interface has also been written for the model code.
242
DCE RPC has been built with cross-platform capability in mind. Cross-platform
means the client and server do not have to be on the same operating system or
hardware type. Therefore, a standard must be created for the transport of data
between the client and server so that when the data is reassembled, it is identical
to the original that was sent. If the client has Little Endian hardware (least sig-
nificant byte first) and the server is a Big Endian machine (most significant byte
first), there must be mechanisms to correctly assemble the data. This is where
data marshaling in needed. Data marshaling takes the data being sent from one
machine to another and translates it into a standard byte stream format. Then the
bytes can be reasembled into the original data using the appropriate standard for
The code to marshal the data is created by an IDL compiler. The compiler
reads the IDL file and creates a client stub that exposes the functions to the client
program, marshals the data and sends the function call and data to the remote
machine. The compiler also creates the server stub that calls the function on the
remote machine. The server stub also takes the data sent to it and formats it into
TCP or UDP which are protocols in the TCP/IP family of network protocols. In
The server portion of the code is a stand-alone program that runs on every
computer that will be used in the parallel computation. Since it would be very
inconvenient to start each program manually, the code is written as a system ser-
vice and network listener. That way, the code is always ready to do a calculation
but does not use any computer time when not in use.
The code is loaded into memory whenever the computer boots and does not
shut down again until the computer is shutdown. It also listens on a TCP port for
requests to calculate an integral. Since the computer can possibly receive requests
from several computers simultaneously, or the server computer may have more
than one CPU, the code is multi-threaded. This means there can be more than
When the code is loaded into memory, a startup routine is called which then
task at startup is
1. Select a protocol for network transport. This is the language used to com-
municate data between the client program and the server program.
2. Register the existence of the listener with the RPC endpoint mapper.
Then the listener is registered with the RPC endpoint mapper with calls to the
if (status)
10 goto cleanup;
// register the endpoints
12 status = RpcEpRegister(CPInt_v4_0_s_ifspec,
pBindingVector,
14 NULL,
NULL);
245
Finally, the code tells the operating system that all is well and the program is
listening for requests as shown in listing C.9. There is one more step in the code,
setting the priority class for the service in line 6 of listing C.9. It is set to idle so
that the code does not interfere with other processes the user of the machine may
be using. In essence, the code will only use what would otherwise be unused
CPU cycles.
Source Listing C.9: Code to start listening for integration calculation requests.
ReportStatusToSCMgr(
2 SERVICE_RUNNING, // service state
NO_ERROR, // exit code
4 0);
// set the priority
6 SetPriorityClass(GetCurrentProcess(),
IDLE_PRIORITY_CLASS);
8 //start listening for calls
status = RpcServerListen(cMinCalls,
10 cMaxCalls,
fDontWait );
for the code to respond to requests, extra threads are created that actually run the
functions.
the calculation engine of the main program. It consists of two classes that work
together to complete all of the integration with the least amount of local CPU
effort. The first of the two classes (CDomainData) will have only one instance
during the life of the integration process, while there will be one instance of the
246
second class (CRpcIntegrate) for every CPU on every machine involved with
the integration.
The integration process starts by creating a list of “jobs” to be done. Each job
over which the integrals are to be computed are triangles or rectangles. Each
job is stored in a last-in-first-out (LIFO) stack saved to disk. This task is carried
out by the first class described above (CDomainData). Once the data has been
sorted and placed as discrete jobs in a stack, the main program spawns off a
new thread for each CPU. Each thread creates an instance of CRpcIntegrate
with the parameters of the machine name that the thread should send integration
Source Listing C.10: Code to spin off integration threads for each CPU on each
machine.
void CLinearSystem::AssembleMatrices(CBEM *pBEM) {
2 ...
// start a thread for each machine
4 POSITION pos=pDB->GetHeadPos();
while(pos!=NULL)
6 {
CComputerRecord *pRec=pDB->GetNext(pos);
8 CString str=pRec->GetName();
pData->bUseHighSpeed=FALSE;
24 pData->bWantSing=FALSE;
DWORD dwThreadID;
26 pThread=CreateThread(NULL, 0,
(LPTHREAD_START_ROUTINE)IntegrationControl,
28 pData, CREATE_SUSPENDED, &dwThreadID);
// set thread priority high
30 SetThreadPriority(pThread, THREAD_PRIORITY_ABOVE_NORMAL);
ResumeThread(pThread);
32 pThreadList.AddTail(pThread);
}
34 }
BOOL done=FALSE;
36 // now monitor
unsigned int count=0;
38 while(!done)
{
40 count++;
Sleep(250); // sleep a quarter second
42 BOOL bStillAlive=FALSE;
DWORD bStatus;
44 POSITION pos=pThreadList.GetHeadPosition();
while(pos!=NULL)
46 {
HANDLE pT=pThreadList.GetNext(pos);
48 GetExitCodeThread(pT , &bStatus);
if(bStatus==STILL_ACTIVE )
50 bStillAlive=TRUE;
}
52 if(!bStillAlive)
done=TRUE;
54 }
}
Listing C.10 shows the code that creates the threads that send integration jobs
to remote machines. Line 26 contains the code that starts a new thread in a sus-
pended state. This allows the main thread to set the thread priority before its
execution is resumed. The high priority level given to the service threads may
seem as though they should cause poor execution of other tasks in the program.
248
That is not true, however, because most the the time the thread is blocked waiting
is to talk to. It then goes into a loop consisting of three tasks: 1) get a job from
CDomainData, 2) make a remote procedure call, 3) If the call succeeds, write the
results to the matrices; otherwise write the job to the error stack for some other
thread to complete.
32 LeaveCriticalSection(&CRPCIntThread::m_csDataLock);
}
34 RpcExcept(1) {
ulCode = RpcExceptionCode();
36 fContinue = 0;
}
38 RpcEndExcept
if(status!=RPC_S_OK)
40 return -1;
BOOL bDone=FALSE;
64 BOOL bGotSing;
INTDATA intData; // data to be sent
66 JOBSTACK job; // job to do
CBemSpecs *pSpecs;
68 CSoilZoneDataBase *pSZDB;
EnterCriticalSection(&CRPCIntThread::m_csFileLock);
70 {
pSpecs=m_pThreadData->GetBemSpecs();
72 pSZDB=pSpecs->GetSoilZoneDB();
}
74 LeaveCriticalSection(&CRPCIntThread::m_csFileLock);
CSoilZoneRecord *pRec;
76 int ii;
250
while(!bDone)
78 {
EnterCriticalSection(&CRPCIntThread::m_csDataLock);
80 {
if(!m_pThreadData->GetJob(job, &bGotSing,
82 m_bWantSing))
{
84 LeaveCriticalSection(
&CRPCIntThread::m_csDataLock);
86 bDone=TRUE;
break;
88 }
}
90 LeaveCriticalSection(&CRPCIntThread::m_csDataLock);
BOOL bResultsGood=TRUE;
92 RESULTS theResults; // the returned data
// make the RPC call
94 RpcTryExcept
{
96 theResults = IntegrateElement(hCP3D, intData);
}
98 RpcExcept(1) {
ulCode = RpcExceptionCode();
100 fContinue = 0;
bResultsGood=FALSE;
102 HandleRPCError(pSpecs, job, intData, hCP3D);
EnterCriticalSection(&CRPCIntThread::m_csDataLock);
104 {
// check for quit message
106 if(m_pThreadData->ShouldStopCalc())
{
108 LeaveCriticalSection(
&CRPCIntThread::m_csDataLock);
110 return 0;
}
112 }
LeaveCriticalSection(&CRPCIntThread::m_csDataLock);
114 }
RpcEndExcept
116 if(bResultsGood)
{
118 EnterCriticalSection(&CRPCIntThread::m_csDataLock);
{
120 m_pThreadData->WriteFinishedWork(job,
&theResults);
251
122 }
LeaveCriticalSection(&CRPCIntThread::m_csDataLock);
124 }
}
126 RpcBindingFree(&hCP3D);
128 return 0;
}
cause the technique supplies a relatively small packet of information that can be
2. NTSB, Pipline Accident Summary Report: Pipeline Rupture, Liquid Butane Re-
lease, and Fire, Lively, Texas August 24 1996, Technical report, National Trans-
portation Safety Board, Washington, D.C. (1998).
3. J. Morgan, Cathodic Protection, 2nd edition (Houston, TX: NACE, Interna-
tional, 1993).
8. C. Kasper, “The Theory of the Potential and the Technical Practice of Elec-
trodeposition: IV. The Flow between and to Circular Cylinders.” Transactions
of the Electrochemical Society, 78 (1940) 147–161.
252
253
25. D. A. Jones, Principles and Prevention of Corrosion (Upper Saddle River, NJ:
Prentice Hall, 1996).
254
28. I. Margarit and O. Mattos, “About Coatings and Cathodic Protection: Pos-
sibilities of Impedance as Monitoring Technique,” in Electrochemical Methods
in Corrosion Research VI, Parts 1 and 2 (Zurich-Eutikon: Transtec Publications,
1998) 279–292.
29. I. Thompson and D. Campbell, “Interpreting Nyquist Responses From De-
fective Coatings on Steel Substrates,” Corrosion Science, 36 (1994) 187–198.
38. H. M. Deitel and P. J. Deitel, C++: How to Program, 3rd edition (Upper Sadle
River, NJ: Prentice Hall, 2000).
39. B. Eckel, Thinking in C++, Volume 1, 2nd edition (Upper Saddle River, NJ:
Prentice Hall, 2000).
40. A. W. Peabody, Control of Pipeline Corrosion (Houston, TX: NACE, 1978).
41. K. E. Jeffers, Electrochemical Impedance Spectroscopy for the Characterization of
Corrosion and Cathodic Protection of Buried Pipelines, Master’s thesis, Univer-
sity of Florida, Gainesville, Florida (1999).
42. D. R. Baker, M. W. Verbrugge, and J. Newman, “A Transformation for the
Treatment of Diffusion and Migration. Application to Stationary Disk and
Hemispherical Electrodes,” Journal of Electroanalytic Chemistry, 314 (1991)
23–44.
43. D. R. Baker, “Reducing Nonlinear Systems of Transport Equations to
Laplace’s Equation,” SIAM Journal of Applied Math, 53 (1993) 419–439.
44. B. Pillay and J. Newman, “Modeling Diffusion and Migration in Dilute Elec-
trochemical Systems Using the Quasipotential Transformation,” Journal of
the Electrochemical Society, 140 (1993) 414–420.
45. M. W. Verbrugge, D. R. Baker, and J. Newman, “Dependent-Variable Trans-
formation for the Treatment of Diffusion, Migration, and Homogeneous Re-
actions,” Journal of the Electrochemical Society, 140 (1993) 2530.
46. K. Allahar and M. Orazem, “Modeling Disbonded Coatings using the
Quasipotential Transformation,” in Passivity and Localized Corrosion, R. G.
Kelly, B. MacDougall, M. Seo, and H. Takahashi, editors, volume PV-99-27
(Pennington, N.J.: Electrochemical Society, Inc., 1999) 609.
47. H. S. Carslaw and J. C. Jaeger, Conduction of Heat in Solids (New York: Oxford
University Press, 1959).
48. H. F. Moulton, “Current Flow in Rectangular Conductors,” Proceedings of the
London Mathematical Society, ser. 2 (1905) 104–110.
49. F. Bowman, Introduction to Elliptic Functions with Applications (New York:
John Wiley and Sons, 1953).
50. M. E. Orazem and J. Newman, “Primary Current Distribution and Resis-
tance of a Slotted-Electrode Cell,” Journal of the Electrochemical Society, 131
(1984) 2857–2861.
51. M. E. Orazem, “Calculation of the Electrical Resistance of a Compact Ten-
sion Specimen for Crack-Propagation Measurements,” Journal of the Electro-
chemical Society, 132 (1985) 2071–2076.
256
52. M. E. Orazem and W. Ruch, “An Improved Analysis of the Potential Drop
Method for Measuring Crack Lengths in Compact Tension Specimens,” In-
ternational Journal of Fracture, 31 (1986) 245–258.
87. W. Sun and N. Zamani, “An Adaptive h-r Boundary Element Algorithm for
the Laplace Equation,” Intenational Journal for Numerical Methods in Engineer-
ing, 33 (1992) 537–552.
88. N. Kamiya and M. Koide, “Adaptive Boundary Element for Multiple Subre-
qions,” Computational Mechanics, 12 (1993) 69–80.
259
105. R. S. Wright and M. Sweet, OpenGL SuperBible, 2nd edition (Waite Group,
1999).
106. M. Woo, J. Neider, T. David, D. Shriner, and T. Davis, OpenGL(r) 1.2 Program-
ming Guide, Third Edition: The Official Guide to Learning OpenGL, Version 1.2,
3rd edition (Addison-Wesley, 1999).
131. J. Prosise, Programming Windows With MFC, 2nd edition (Microsoft Press,
1999).
132. “Visual C++ Enterprise Edition version 6.0,” Software and Documentation
(1998).
BIOGRAPHICAL SKETCH
the author started attendance at the University of Florida, under the direction
Following graduation the author will take up residence in New Jersey where he
263