UFDC Home  myUFDC Home  Help 



Full Text  
MODEL FOR INTERPRETATION OF PIPELINE SURVEY DATA By CHENCHEN QIU A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2003 ACKNOWLEDGMENTS I would like to express my sincerest gratitude and appreciation to Dr. Mark E. Orazem for his support, inspiration and guidance in directing this study. His hard work and attention to detail were an excellent example during my studies. He always helped me to learn more and taught me the importance of thinking creatively. Many of the advances achieved during my graduate study would not have been possible without his support. I would also like to gratefully acknowl edge Dr. Loc Vu Quoc for teaching me the boundary element method at the early stage of this work. Moreover, I wish to thank Dr. Anthony J. Ladd, Dr. Oscar D. Crisalle, and Dr. Darryl Butt for their time, useful discussions and guidance as members of the supervisory committee. I would like to acknowledge the fi nancial support of the Pipeline Research Council, International and Gas Research Institute. My colleagues, Douglas P. Riemer, Kerry Allahar, Nelliann PerezGarcia, Pavan Shukla, who have contributed to this research by their valuable discussions and friendship, are also gratefully acknowledged. Finally, I would like to thank my husband, Kuide Qin, for his unselfish support and understanding throughout my doctorate work, which would not have been possible without him. I would also like to thank my parents and my sister for their encouragement and inspira tion. Even though they were half a planet away from me, they have always been with me in my heart. TABLE OF CONTENTS page ACKNOWLEDGMENTS ............................ LIST OF TABLES ................................ LIST OF FIGURES ............................... ABSTRACT .................. N O TATIO N . . . . . . . . . CHAPTER 1 INTRODUCTION ............................. 2 PIPELINE CORROSION, PROTECTION AND MEASUREMENTS . 2.1 Introduction . . . . . . . . 2.2 Basic Concepts of Corrosion ....................... 2.3 Corrosion Protection ........................... 2.3.1 C oating . . . . . . . . 2.3.2 Cathodic Protection ....................... 2.4 Cathodic Protection Criteria ....................... 2.4.1 850mV Potential Criterion ................... 2.4.2 100mV Potential Criterion ................... 2.5 Pipeline Measurements .......................... 2.5.1 Potential Survey ......................... 2.5.2 Line Current Survey ....................... 2.5.3 Other Survey Techniques ................... 2.6 Conclusions . . . . . . . . 3 APPLICATION OF TWODIMENSIONAL FORWARD AND MODELS TO CORROSION ................... 3.1 Introduction . .. 3.2 Thin Plate Method . . 3.2.1 Introduction . . INVERSE . . 16 . . 16 . . 17 . . 17 . . . . . . . . . . . . 3.2.2 Evaluation Process .................. ....... 18 3.2.3 Summary ....... ...... ... ....... 20 3.3 TwoDimensional Boundary Element Method ............. 20 3.3.1 Introduction .......... ......... ...... 20 3.3.2 Evaluation Process .................. ....... 22 3.3.3 TwoDimensional Problem .... . . . 22 3.4 Inverse Analysis in Two Dimensions . . . 24 3.4.1 Objective Function .................. ....... 24 3.4.2 Regression Method Analysis . . . 25 3.4.3 Downhill Simplex Method .................. ..25 3.4.4 Accuracy of the Parameters ........ ......... 26 3.4.5 Regression Results .................. ..... 28 3.5 Conclusions .. ........................... 29 4 DEVELOPMENT OF THREEDIMENSIONAL FORWARD AND INVERSE MODELS FOR PIPELINE WITH POTENTIAL SURVEY DATA ONLY .30 4.1 Introduction .............. ............... 30 4.2 Construction of the Forward Model .............. 33 4.2.1 Cathodic Protection System .......... ....... 33 4.2.2 Governing Equation and Boundary Conditions ...... 34 4.2.3 Theoretical development ......... ....... .. 35 4.3 ThreeDimensional Boundary Element Method ............ 38 4.3.1 Infinity Domain ........... .............. 38 4.3.2 HalfInfinity Domain .......... ....... ...... 41 4.3.3 Boundary Discretization ....... .... .......... 44 4.3.4 Coordinates Definition and Transformation ......... 45 4.3.5 Discretization of Boundary Element Method ........ 46 4.3.6 Row Sum Elimination ........... ...... ...... 49 4.3.7 SelfEquilibrium ........... ............. 50 4.4 Forward Model ................... ........ 52 4.4.1 Constant Steel Potential Assumption .............. 52 4.4.2 Simulation Results ........... ...... ....... 54 4.5 Inverse Model ............................... 54 4.5.1 Objective Function ........... ...... ....... 55 4.5.2 Analysis of Regression Methods ........ ...... 55 4.5.3 Simulated Annealing Method ....... ........ 56 4.5.4 Simulation Results and Discussions .............. 57 4.5.5 Inverse Strategies ........... ............. 62 4.6 Conclusions .............................. 70 5 DEVELOPMENT OF THREEDIMENSIONAL FORWARD AND INVERSE MODELS FOR PIPELINE WITH BOTH POTENTIAL AND CURRENT DENSITY SURVEY DATA ......................... 72 5.1 Forward Model .............................. 72 5.1.1 Introduction ........ ......... ......... 72 5.1.2 Pipeline with Varying Steel Potential .............. 72 5.1.3 Simulation Results ................... ....... .. 82 5.1.4 Analysis of the Simulation Results . . ..... 86 5.1.5 Validation with CP3D ............... . 87 5.2 Inverse Model ........... .... ......... 91 5.2.1 Introduction ........... .... ....... 91 5.2.2 Regression to NoiseFree Data . . . 93 5.2.3 Regression to Noisy Data ... . . . 99 5.3 Conclusions ............. . . ... 102 6 APPLICATION OF OFFPOTENTIAL DATA TO THE INVERSE MODEL 104 6.1 Forward ModelCalculation of the Soil Surface OffPotential . 104 6.1.1 Physical Process .......... ....... ....... 104 6.1.2 Assumptions ........... ............... 105 6.1.3 Model Implementation ... . . .105 6.1.4 Simulation Results and Analysis . . . ... 106 6.2 Inverse ModelUsing Three Kinds of Heterogeneous Data Sets 107 6.2.1 Data without Noise .......... .......... ... 107 6.2.2 Data with Added Noise ................ 109 6.3 Conclusions .......... .......... .......... .. 112 7 SUMMARY AND FUTURE WORK . . . .. ... 113 7.1 Summary ........ ........... ...... ........ 113 7.2 Future Work ................................ 114 APPENDIX A BOUNDARY INTEGRAL ......................... 115 A.1 Integral along One Object .......... ....... ........ 115 A.1.1 Integral over Cylindrical Elements ......... ...... 115 A.1.2 Integral over Circle Elements ............. 116 A.2 Integral between Different Objects ............... 117 A.2.1 Integral over Cylindrical Elements ......... ...... 118 A.2.2 Integral over Circle Elements ............. 119 B CODESTRUCTURE ................... ....... 121 B.1 Code Structure of the Forward Model ............... 121 B.2 Code Structure of the Inverse Model ................. 121 C INTERFACE ... ............... .. ......... 124 C.1 Input Windows .............................. 124 C.2 Output Windows ................... ......... 124 REFERENCES .................................. 127 BIOGRAPHICAL SKETCH ........................... 134 LIST OF TABLES Table page 3.1 Regression result for twodimensional inverse model described by Aoki et al. The underlined symbols represent the parameters esti mated by regression. .................. ....... .. 28 4.1 Parameter values obtained using the threedimensional inverse model developed in the present work for a 10m pipe segment with two coating defects. .................. ........... .. 60 4.2 Test case parameters with five coating defects on the pipe used to demonstrate the method for determination of the number of statis tically significant parameters (see Figure 415). The intact coating resistivity had a value of 5.0 x 107 / m. . . . 62 4.3 Regression results from the threedimensional inverse model for a 100m pipe segment with three coating defects. The sequential re gression procedure was used to identify the number of defects that were statistically significant. The intact coating resistivity had a value of 5.0 x 107 / Qm. ................. ..... .. 64 5.1 Regression results using either equation (519) or (518) for homoge neous synthetic data without added noise or equation (520) for het erogeneous synthetic data without added noise. The initial values for each defect was Xk = 50 m, pk = 3.5 x 107 Qm, and ok = 0.92 m. 95 5.2 Regression results using equation (524) for heterogeneous synthetic data with added noise. The initial values for each defect was Xk = 50 m, pk = 3.5 x 107 Qm, and ak = 0.92 m. . . .... 101 6.1 Regression results from the data without added noise obtained by using the general form of the objective function (525), (526) and (527). The initial values for each defect was Xk = 50 m, pk = 3.5 x 107 Qm, and ok = 0.92 m. ................ ..... .. 108 6.2 Regression results from the noiseadded data by using general form of the objective function (525), (526) and (527). The initial values for each defect was Xk = 50 m, pk = 3.5 x 107 Qm, and ok = 0.92 m. 111 LIST OF FIGURES Figure page 21 Polarization curve of the cathodic protection. .............. 8 22 Sacrificial electrode CP ............... ........ 9 23 Impressed current CP ............... ........ 9 24 Potential survey method. ............. .. ..... 12 25 Pipetosoil potential survey method and importance of the halfcell location. ................. ...... ......... 13 26 Line current survey. .............. . . .... 14 31 Finite difference method. ................ ........ 19 32 Comparisons of the forward and inverse results of TPA method. .. 20 33 Schematic representation of the twodimensional problem. . 22 34 Reproducing results for the twodimensional forward model. . 23 35 Schematic representation of the twodimensional inverse model. .. 24 36 Schematic representation of the downhill simplex method. . 26 41 Schematic illustration of the pipe segment and anode used to test the inverse model. ............. . . ...... 33 42 Linear relationship between potential drop and current density over the pipe coating .................. ............ 36 43 Schematic illustration of the pipe surface resistivity model. . 38 44 Interior problem. .................. .......... .. 39 45 Multi connected region of interior problem. . . ... 40 46 Schematic illustration of mirror reflection technique. . ... 43 47 The Fundamental solution to the halfinfinity domain satisfying the Neuman b.c's. ............... . . ..... 43 48 Pipe discretision and collocation points. . . . 44 49 Coordinate rotation .... ........... ....... .. 46 410 False color image of the onpotential on the soil surface that was generated by the forward model corresponding to Figure 41. . 55 411 Flow chart for the inverse model calculations. . . .... 58 412 Grid showing the location of 303 surface onpotentials calculated using the threedimensionalforward model developed in the present work. The grid shown is for a 10m pipe segment. A scaled version of the grid was used for a 100m pipe segment. . . ... 59 413 The regression objective function as a function of the number of evaluations for a pipe coating with one defect region. The simu lated annealing method was used for this regression. . ... 60 414 Comparison of the set and fitted results for pipe coating with two coating defects. .................. ........... .. 61 415 The regression statistic as a function of the number of coating de fects assumed for the model. The minimum in this value is used to identify the maximum number of coating defects that can be justi fied on statistical grounds. ................... ...... 63 416 Comparison between the input values and regression results for noise = 0.1 mV. ................ ..... .. .. .. .. 66 417 Comparison between the input values and regression results for noise = 1.0 mV. ........... . . ...... 67 418 Comparison between the input values and regression results for noise = 2.0 mV. ........... . . ....... 69 51 Cathodic protection system with variant potential along the pipe and anode. ................. . . .... 73 52 Relation of potential V and current density. . . ... 74 53 Axial direction current flow along the pipeline. . . .. 75 54 Axial direction current flow along the anode. . . .... 76 55 Potential and current density along the anode. . . ... 83 56 Simulated axial distributions along the pipeline. . . ... 85 57 Potential of pipe calculated by using CP3D .. . . ... 88 58 Potential of pipe steel calculated by using CP3D . . .... 89 59 Axial direction current density of pipe calculated by using CP3D 90 510 Radial direction current density of anode calculated by using CP3D 90 511 Schematic illustration of the coupling of a soilsurfacelevel poten tial survey method to an aerial magnetometer current survey to as sess the condition of a buried pipeline. . . . ..... 92 512 Synthetic pipe current and surface potential data corresponding to Table 5.1. Dashed lines indicate the location of coating defects. 93 513 The X2/ statistic corresponding to Table 5.1 as a function of the number of coating defects assumed in the regression. . ... 96 514 The X2/ statistic corresponding to regression of equation (521) (see Table 5.1) as a function of the number of coating defects as sumed in the regression. ................. ..... .. 98 515 Synthetic data with added noise. The line represents the noisefree values and the symbols represent synthetic data with added noise. a) surface onpotential value with = 0.4 mV; b) axial current den sity value with proportional noise corresponding to 2 percent of the value ................. ........ .......... 100 61 Switch off CP power source to determine instant offpotential. 105 62 Onpotential and offpotential on soil surface. . . ... 106 63 Typical closeinterval potential graph. . . . 107 64 Decide the number of defects using noise free soil surface on and offpotential and current density data. . . . 109 65 Set value and noiseadded offpotential on soil surface data. . 110 66 Decide the number of defects using noiseadded on and offpotential and current density data. .................. ....... 111 A1 Integral between different objects. .. . . . 118 B1 Code structure of the forward model. ..... . . . 122 B2 Code structure of the inverse model. .... . . . 123 C1 Window for pipe properties. ................... .. 125 C2 Window for anode properties ............. . 125 C3 Window for coating properties ...... . . . 126 C4 Window for figures .................. . .. 126 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy MODEL FOR INTERPRETATION OF PIPELINE SURVEY DATA By Chenchen Qiu December 2003 Chair: Mark E. Orazem Major Department: Chemical Engineering Pipeline corrosion annually costs the U.S. seven billions of dollars. Ninety percent of the cost is the capital, operation and maintenance. It is desirable to optimize this expenditure while ensuring that the integrity of the pipeline, which represents a valuable asset, is maintained. Regular inspection of pipelines using a variety of survey techniques is routinely conducted, but interpretation of the results is confounded by stochastic and systematic errors. The objective of the present work was to investigate use of inverse model that could interpret survey data in the context of the physics of the system. An inverse analysis model was developed which provides a mathematical framework for interpretation of survey data in the presence of random noise. A boundaryelement forward model was coupled with a weighted nonlinear regres sion algorithm to obtain pipe surface properties from two types of survey data: soilsurface potentials and local values of current flowing through the pipe. The forward model accounted for the passage of current through a threedimensional homogeneous medium and yielded soil surface potentials for given pipe/anode configurations and pipe coating properties. The number of regressed parameters was reduced by using a function for coating resistivity that allowed specification of coating defects. A weighted simulateannealing nonlinear regression algorithm facilitated analysis of noisy data. A method to determine the appropriate number of fitted parameters was developed. The model was demonstrated for synthetic data generated for a section of a coated underground pipeline electrically connected to a vertical sacrificial anode. The success of the regression was sensitive to the relative weighting applied in the objective function to the respective types of data. A generalized weighted and scaled objective function was proposed. NOTATION Potential and potential distribution, V 4 Potential and potential distribution refer to the reference electrode, defined as V V Pk Resistance reduction of the coating defect k, Q.m Jk The half width of the coating defect k, m if Normal vector of boundary Cj Concentration of species j, mol/cm3 Di Diffusion coefficient of species j, cm2/s F Faraday's constant, 96487C/equiv i Current density, A/m2 Uj Mobility of species j, cm2 mol/J s V Voltage distribution of pipe steel or anode metal, V v Bulk velocity, cm/s Xk Center point of the coating defects k, m zi Charge number of species j p Resistivity of the pipe coating, Q.m CHAPTER 1 INTRODUCTION In the U.S., there are 779,000 km (484,000 mi) of gas and liquid transmission pipelines. Seven billion of dollars is spent due to pipeline corrosion every year. 90% of the cost is the capital, operation and maintenance. It is desirable to op timize this expenditure while ensuring that the integrity of the pipeline, which represents a valuable asset, is maintained. Regular inspection of pipelines using a variety of survey techniques is routinely conducted, but interpretation of the results is confounded by stochastic and systematic errors. The goal of this work was to develop a computer program that can be used to extract the condition of the pipe taking into account a wide variety of field data, such as potential and line current survey data. In addition, the intention of this preliminary effort was to establish a proof of concept which can motivate future development. The pipeline protection techniques that are used most often are cathodic pro tection and coating. The normal use of computer programs for modeling of ca thodic protection is to assess whether a pipeline with an assumed coating condi tion can be protected by a given cathodic protection design. In this strategy, called a forward solution, all physical properties of the pipe, anodes, ground beds, and the soil are assumed, and the corresponding distribution of current and potential on the pipe is calculated. Overprotected and underprotected sections of the pipeline can be identified in this way. Such programs can also be used to calculate the po tential and the current density at other locations such as at the soil surface or at 2 the specified locations where coupons may be buried. This strategy lends itself well to answering "what if" questions. For example, "what will happen if discrete coating holidays expose five percent of the pipe surface?", or "what will happen if a new pipe is constructed in the rightofway?" The forward strategy is not appropriate, however, for the interpretation of field data to assess whether a given pipe is protected. The concept behind the approach developed in this work is to solve the inverse problem in which the properties of the pipe coating are inferred from measurements of the current and the potential distributions. This approach would allow interpretation of field data in a manner that would take into account the physical laws that constrain the flow of electrical current from anode to pipe. In addition, the development of the inverse model can not proceed alone. Only when the forward model is well developed can the inverse model be able to inter pret multiple kinds of survey data. The basic concepts of the corrosion, corrosion protection methods, and pipeline survey methods are presented in Chapter 2. Verification of published methods for model twodimensional systems are pre sented in Chapter 3. This chapter includes the thinplate method and the bound ary element method to solve the twodimensional Laplace's equation. The two dimensional inverse model is also established in this chapter. This work involved the construction of the objective function, the formalism of regression strategies, and regression methods analysis. The development of the forward and the inverse models for twodimensional systems provided insights into the methods needed to address the inverse prob lem for pipelines. However, the twodimensional forward model cannot describe precisely the geometry of the pipeline cathodic protection system. In addition, the 3 more sophisticated threedimensional forward models such as CP3D,1 PROCAT,2 and OKAPPI3 are not appropriate, since they require too much computation work for the inverse model. A simplified version of a forward is presented in Chapter 4. Both the forward model and inverse model are applied to soil surface potential sur vey data. For the forward part, a threedimensional boundary element method with cylindrical elements was developed in this work. The simulated annealing method was shown to be ideally suitable for the threedimensional inverse model. In addition, several strategies were explored to assess the confidence level of the inverse model results. The forward model presented in Chapter 4 provided the foundation for an ex tended model, presented in Chapter 5, which accounted for the potential drop along the pipe steel. Both the potential on soil surface and the current density along the pipeline could be obtained from the forward model. An inverse model was created to apply these data. A new general form of the objective function was established by adding different weighting factors for each data point and by including the number of the data points. Chapter 6 verifies the general form of the objective function built in Chapter 5. Offpotential on the soil surface was simulated in the forward model. Therefore, three kinds of data sets, i.e., onpotential, offpotential and current density data, were used for the inverse model. Summary of the study and future work are presented in Chapter 7. CHAPTER 2 PIPELINE CORROSION, PROTECTION AND MEASUREMENTS 2.1 Introduction The corrosion of metallic structures has significant impact on the U.S. economy in many fields. In 1975, BattelleNBS's benchmark study estimated that the cost of corrosion was $70 billion per year, 4.2 percent of the nation's gross national product (GNP). If the effective and presently available corrosion technology can be applied, $10 billion of this cost could be avoided. From 1999 to 2001, the Federal Highway Administration (FHWA) indicated that the total direct cost of corrosion was about $276 billion per year, 3.1 percent of the U.S. gross domestic product (GDP). The indirect costs is unpredictable .4 Pipelinerelated corrosion costs approximately $7.0 billion annually. There are 528,000 km (328,000 mi) of natural gas transmission pipelines, 119,000 km (74,000 mi) of crude oil transmission pipelines and 132,000 km (82,000 mi) of hazardous liquid transmission pipelines in the U.S.4 In the past few years, a number of gas and liquid pipeline failures have drawn people's attention to the pipeline safety. In order to preserve the asset of the pipeline and to avoid those failures which may jeopardize public safety, result in product loss, or cause property and environmental damage, some new regula tions were established. Regular pipeline inspections, such as hydrostatic testing, direct assessment, and inline inspection (ILI), are required. Furthermore, corro sion prediction models need to be developed in order to determine and prioritize the most effective corrosion preventive strategies. Development of an understanding of cathodic protection of buried structures such as pipelines requires an understanding of the electrochemical reactions asso ciated with corrosion. This section provides an overview of the principles of corro sion, methods for cathodic protection, and the criteria used to assess the effective ness of corrosionprevention strategies. Moreover, the pipeline survey methods are included. 2.2 Basic Concepts of Corrosion Corrosion takes place in response to the tendency to reduce the overall free en ergy of a system. For example, when metal is in a dilute aerated neutral electrolyte atmosphere, moisture films, which contain oxygen and water, usually cover the metal surface. Some atoms in the metal surface tend to give out electrons and be come ions in the moisture layer. In this way, they are in a lower energy state in solution than when they are in the lattice of the solid metal. At the same time the metal surface builds up a huge excess negative charge. These electrons move out from the metal and attach themselves to protons (H+ ions) and molecules oxygen (02). The process repeats itself on various parts of the surface. Then the metal dissolves away as ions.5 The primary constituent of pipelinegrade steels is iron. Therefore the overall corrosion reaction can be written in terms of dissolution of iron, e.g., 2Fe + 02 + 2H20  2Fe(OH)2 (21) Reaction (21) can be considered as the result of two halfcell reactions. One is Fe Fe2+ + 2e (22) This is an oxidation reaction with an increase of oxidation state for iron from 0 to 2+. It is called an anodic reaction. The other is 2H20 + 02 + 4e 40H (23) This reaction is a reduction reaction with a decrease of oxidation state for 02 to OH It is defined to be a cathodic reaction. The change in free energy associated with reaction (21) can be expressed in terms of a cell potential as AG = nFAE (24) where n is the number of electrons exchanged in the reaction, F is Faraday's con stant, 96487 coulombs/equivalent,6 and E is the electrochemical potential. Ac cording to the Nernst equation, obtained by neglecting both activity coefficients and liquidjunction potential,7'8 RT RT AE = AE ln(H(a~')) = AEo ln(H(c~')) (25) nEE where EO is the equilibrium potential, and ai is the activity of the species i, si is the stoichiometric coefficient of species i, and ci is the concentration of the species i. For reaction (22), the potential in Nernst equation form is RT E, = Eoe2+Fe 2 n ([Fe2+]) (26) Likewise, for reaction (23), the potential in Nernst equation form is RT In [OH ]4 (27) Ec = E2/OH n P2 where a means anodic reaction and c refers to the cathodic reaction. To the whole reaction, where n, = 2 and nc = 4, AG = 2AGa + AGc = 2naFEa ncFEc = 4F(E, + Ec) (28) For the two halfcell reactions, E, and Ec are positive. Hence, the free energy of the overall reaction must be negative, which means the corrosion reaction occurs spontaneously. 2.3 Corrosion Protection To control the corrosion situations, several methods have been proved effec tive, which include coating and cathodic protection. 2.3.1 Coating Electrically insulating materials can be use to cover the pipe surface and thereby isolate the pipe metal from contact with the surrounding electrolyte. In this way, a high electrical resistance in the anodecathode circuit is added. Theoretically, since no significant corrosion current flows from the anode to the cathode, no cor rosion would occur. Many kinds of materials have been applied as coatings, such as enamels, tapes, and plastic coating. Recent improvements in pipeline coating materials will also reduce the risk of a corrosionrelated failure. However, a recent survey of major pipeline companies indicated that about 30% of the primary loss of the pipeline protection was due to coating deterioration.4 Therefore, coating alone can not provide full protection for the pipeline. Practically, effective pipeline corrosion control comprises use of good coatings along with cathodic protection as a secondary defense. 2.3.2 Cathodic Protection In 1824, Sir Humphrey Davy successfully protected copper against corrosion from seawater by using iron anodes. It was the first application of cathodic pro tection (CP) and, at that time, it had no theoretical foundation. From that begin ning, CP has been applied to marine and underground structures, water storage tanks, pipelines, oil platform supports, reinforcing steel and many other facilities exposed to a corrosive environment. By now, the theory of CP is well established. Natural gas and oil companies have already been using CP for economic, as well as safety reasons, since the 1930's. The Pipeline Safety Act of 1972 made the appli 0.4 Fe Fe2+ + 2e 0.6 0.8 H20 +2e H2+20H S1.0 {.. O_02+H20+4e >40H 1.2 0.01 0.1 1 10 100 Current Density, mA/ft2 Figure 21: Polarization curve of the cathodic protection. cation of CP on pipelines transporting hazardous material mandatory for safety concerns. Cathodic protection results from cathodic polarization of a corroding metal surface to reduce the corrosion rate. In this system, the corrosion potential reduces the rate of the halfcell reaction (22) with an excess of electrons, which drives the equilibrium from right to left. The excess of electrons also increases the rate of oxygen reduction and OH production by reaction (23) in a similar manner during cathodic polarization. Cathodic protection reduces the corrosion rate of a metallic structure by reducing its corrosion potential, bringing the metal closer to an immune state, which is represented by the black curve in Figure 21. Two methods are usually used to achieve this goal.6 There include sacrificial anode and impressed current. Sacrificial Anode. One method to provide CP is to connect a sacrificial metal (with a higher natural electromotive force) through a metallic conductor or a wire to the structure intended to be protected (as Figure 22). Magnesium is a common sacrificial or galvanic anode. This type of galvanic cathodic protection relies on the Ground level Earth Current environment *Pipe 4 Mg anode Figure 22: Sacrificial electrode CP. i  rrenti r Current F f  SCurrent I I Figure 23: Impressed current CP. natural electrical potential between the two metals to cause the cathodic protection current to flow. Since the driving voltage is limited to the very small potential difference existing between the metals, and the current output is relatively low, this type of cathodic protection is normally associated with very small or very well coated structures. ImpressedCurrent. The second common method to provide CP involves use of impressedcurrent. It relies on an external direct current source such as a rectifier or battery (as Figure 23). An anode material, placed in the electrolyte with the protected structure, is made more positive than the structure by connecting both the anode and the structure to the direct current supply. Any conductive materials can be utilized as an impressed current anode, but since corrosion takes place at the anode, materials with very low consumption rates are most desirable. 2.4 Cathodic Protection Criteria As the cathodic protection is widely used, the reliable criteria need to be devel oped to direct the cathodic protection systems to an optimum level. The criteria are based on the potential differences between the pipeline and its environment. Two criteria will be introduced in this section, 850mV and 100mV potential crete ria. 2.4.1 850mV Potential Criterion The criteria for protecting steel pipe, RP1069, was established by NACE in 1969. It was revised later in 1972, 1983 and again in 1996. It defines a negative po tential at least 850mV with respect to a saturated copper/copper sulfate reference electrode(CSE) when CP is applied. This criterion is mostly used since the poten tial measurements with current applied requires minimum equipment, personnel and takes less time to obtain in the field survey. Normally, the desired potential range of the cathodic protection for a pipeline is between 850mV and 1200mV measured with respect to a Cu/CuSO4 refer ence electrode located at the surface of the steel or coating defects.9 If excessive amounts of cathodic protection are applied, the direct reduction of water becomes thermodynamically possible, e.g., 2H20 + 2e , H2 + 20H (29) Hydrogen is evolved upon the metal, and, at large cathodic overpotentials, a small fraction of hydrogen can enter the metal, which makes the metal brittle. This undesirable process is called "hydrogen embrittlement". If the pipeline has a coating, the developing hydrogen gas bubble can exert tremendous pressure. When the pressure is created under the torn coating, a stripping action exposes the metal and results in rapid deterioration. Obviously, if smaller amounts of CP is applied (resulting in potential less negative than 850mV), the structure will not be fully protected. In addition, the 850mV criterion does not directly address the polarization. It is valid only when the effect of IR drop has been considered and eliminated. 2.4.2 100mV Potential Criterion The 100mV polarization criterion was first proposed by S.P. Ewing in 1951.10 The criterion means that the change in potential polarization to protect pipeline is always less than 100mV. To utilize this criterion, the IR drop effect at measurement area needs to be estimated. Normally, coupons are used to get the IR drop values by measuring the potentials when the coupon is "on" and "off", and comparing the two potentials curves to establish the potential differences between the two curves, which is the IR drop value. Recent research in the UK has shown that the applied current density could be reduced by 25% after three month of using the 100mV criterion and by 65% after one year, without losing any protection.11 2.5 Pipeline Measurements Cathodic protection design should be based on actual data. Monitoring meth ods are necessary to investigate the corrosive conditions and to evaluate the de gree of the applied CP. The techniques used include measurement of potential along the ground above the pipe, measurement of current flowing through the pipes, and use of auxiliary electrodes, such as a reference electrode and coupon. 2.5.1 Potential Survey Construction of a potential survey involves measuring the potential difference between the pipe and a reference electrode at the ground surface level. The refer ence electrode is moved to different locations to sample many positions at ground level near the pipe. In this way, a mapping of potential is created which can be Voltmeter CuS04 Electrode Figure 24: Potential survey method. used to identify regions where the pipe is grossly underprotected. The survey involves use of a voltmeter and a reference electrode which is typically based on the Cu/CuSO4 reaction. A schematic representation of the survey method is pre sented in Figure 24.12 Potential survey measures the potentials between buried pipeline and environment.9 This survey needs an appropriate voltmeter with its negative terminal connected to the pipeline and the positive terminal connected to a reference electrode. The procedure involves moving the reference electrode at 3 to 10 ft intervals down the full length of the pipe and using the trailing cable to determine distance and to connect to the pipeline. A vast quantity of data is gathered by a data logger. From the survey readings, the worst corrosion takes place where the readings are the highest (less negative value) and little or no corrosion should take place where the potential is more negative. For this measurement, the location of the reference electrode (halfcell) is im portant.13 The most desirable place is directly above the pipeline. If the reference halfcell is placed 2 3 meters to the side of the pipeline, the measured potential will be a group of scattered data as the second curve shown in Figure 25. In such a case, it is difficult to make a judgement. To avoid this case, measurements on both sides of the pipe may be required to make sure that the over the pipe potentials I:1 I5 Reference electrode placed 2~3 meter to side of pipeline +  Figure 25: Pipetosoil potential survey method and importance of the halfcell location. are taken as close to the pipes as possible. An instant offpotential (IRfree) survey9'14 is used to eliminate the potential drop through the soil in the pipe to soil potential measurement by interrupting the flow of protective current to the pipeline. This technique can be used to determine the effectiveness of the cathodic protection system. It is based on the principle that the IR effect in the potential measurement decays almost instantaneously, while the pipetosoil interface polarization decays relatively slowly, thus allowing the correct pipetosoil potential to be measured free of the IR error. 2.5.2 Line Current Survey The line current survey technique is used to measure the electrical current flow ing on the pipe. If corrosion is taking place on a pipeline, current will flow to the line at some points and flow out of the line at others. For small local cells, the current path may be too short to detect. For large cells, the current may follow the pipe for hundreds or thousands of feet. It is these long line currents that can be detected in a line current survey. 350 300  250 200 150 100 50 IT 0 0 1000 2000 3000 4000 Pipeline Length /feet Figure 26: Line current survey. Because the pipe itself has some resistance to the flow of electric current, there will be a voltage drop in the pipe if current is flowing through the pipe steel. Normally, the voltage drops are very small. Thus the span of the two test points is often fixed at 100 feet to increase the span resistance. Knowing the span resistance of the pipe being surveyed, the voltage drops may be converted to equivalent current by the application of Ohm's Law: E(potential) I(current) = (l (210) R(resistance) The values of current together with the direction of flow then may be plotted as function of line length.9 From Figure 26, it is noted that in one area current flows from both directions toward a particular point on the line. This point must be a place of current discharge. 2.5.3 Other Survey Techniques The airborne cathodic monitoring systems (ACMS)15 detect upset conditions on the pipelines protected by impressed current. Sensitive and filtered magnetic field coils are installed on a helicopter, which continuously measure the magnetic field generated by the ripple from an alternating current source. The principle is based on B = ( (211) 2 R where B is the magnetic field generated by the impressed current I, R is the dis tance from the center of the pipe to the sensor installed on the helicopter, and the constant 0o/27 is the permeability of the medium. Thus, if the magnetic field can be measured along with the distance, then the current can be calculated. In addition, use of an insulated coupon is a variation of the instantoff potential to soil method. The coupon usually is buried close to the protected structures to have the same characteristic of the chemical environment and electrical field as those at the structure surface. The exposed coupon surface simulates a coating defect. Cathodic protection current to the coupon can be interrupted without any effect on the protected structures. In this way, the difficulties of the instantoff potentialtosoil method can be eliminated. 2.6 Conclusions A brief introduction of the corrosion situations was given in this chapter. The basic concepts of corrosion were explained. The pipeline protection techniques, es pecially the cathodic protection technique, have been introduced. In addition, the commonly used pipeline survey methods, such as potential survey, line current survey, airborne cathodic monitoring systems and coupons, are also summarized. CHAPTER 3 APPLICATION OF TWODIMENSIONAL FORWARD AND INVERSE MODELS TO CORROSION 3.1 Introduction Given the governing equation and the boundary conditions to estimate the po tential and current density distribution on the soil surface by solving the govern ing equation directly is known as the forward model. In contrast, in an inverse prob lem, usually a model or a governing equation and measurements of some vari ables are given, such as the potential and current density measurements. Bound ary conditions and the rest of the variables may not be known explicitly. The purpose of the inverse model is to identify the unknown parameters using the in formation from the measurements. It is a type of illconditioned problem in which the solution is extremely sensitive to the measured data. In this chapter, the forward and inverse models applied in the literature are re viewed and reproduced. Two methods were studied, the thin plate approximation method and the boundary element method. The experiences of developing and evaluating these models are summarized. 3.2 Thin Plate Method 3.2.1 Introduction Inglese16,17 solved the inverse problem for two electrode separated by a small gap, a, as following: V2 =0 in 0Q; (x) = 00(x) x [0, 1], y = 0; (31) On(x) + 7(x)9(x) = f(x) x e [0,1], y = a; n (0, y) = O(l, y)= 0 y [0,a] where 7 is the corrosion coefficient, which represents the corrosion rate, 0 is the potential, and ,n is the outward normal flux and a is the thickness of the do main. The goal was to find 7 from the potential 00 on the accessible boundary x [0, 1],y = 0. The author introduced the thin plate approximation (TPA) method. The idea is to perform an expansion of the solution 7 in terms of powers of the thickness of the domain Q. Owing to the assumption that the thickness a is much smaller than 1 (a < 1), higher order terms of the expansion can be ignored. Inglese deduced a formula of 7 as function of ,xx(x, 0), which is the second derivative of 0 with respect to x. #xx(x, 0)+ So(x) rTPA(x) = (X 0) + 0 (32) ((x,0) In order to prove the accuracy of equation (32), the direct problem was solved for a prescribed corrosion rate give by 7(x) = exp(100(x 0.25)2) (33) First of all, for the forward part, O(x, 0) is calculated according to the given 7(x). After that, ,xx, Oo can be obtained without problems. For the inverse part, the corrosion rate 7 was obtained from equation (32) and it can be compared with the value of 7(x) from equation (33). 3.2.2 Evaluation Process The domain of Q in x direction is [0, 1] and in y direction is [0, a]. To simplify the problem, the y direction is normalized to [0, 1]. Thus, the decomposed form of equations (31) is: a29xx + yy = 0 in Q; x(O, Y)= 0 x(,y) = 0 (34) Oy(X, O) = a2 Oy(x, 1) + a27(x)(x, 1) = 0 ay(x) = exp[100(x 0.25)2] To solve equation (34), the Laplace's equation with Robin boundary condi tions, two methods were tested. One was the separation of variables method and the other was the finite difference method. However, the separation variable method failed for this kind of problem because a contradiction occurs between the result and boundary condition Oy(x, 0) = a2 0. To apply the finite difference method as shown in Figure 31, we set o,,j (xi, yj) xi = ih (35) yj = jh Ax = Ay i1 ,j ij i+1 ,j i,j1 Figure 31: Finite difference method. Therefore, Oy,y Ox(1, y) Oy(x, 0) ' (i+l,j  01,j 1,; j i,i i, 1 20ij, + i i,j) 20i,i + i,j+1) Equation (34) becomes: i,j1 + a20i+1i, 2(a2 + 1)oi,j + a20i1, + i,j+l = 0 01,] 01,i = 0 N+ ,] PN1, = 0 i,i 0i,1 = 2ha2 0 2[ i,N1 + Oi,N 1] = a2/(ih)O(i, N) = a(ih)Oi,N i = 0, 1, ... N; I a(x) = exp[100(x 0.25)2] the solution to equation (37) yields Ox,x(x, 0), O(x) and O(x, 0). Due to the norma ization in y direction, correspondingly, equation (32) will be, aoxx(x, 0) + O(x)  (  S (x,0) The results are shown in Figure 32. Where the solid curve represents 7(x) (36) (37) 8) 0,1,.. N; 0, 1, ...N; 0, 1, .. N; al 08 calculatedgamma exactgamma 06 04 02 0 .' ... .. .' . 00 01 02 03 04 05 06 07 08 09 10 Figure 32: Comparisons of the forward and inverse results of TPA method. exp[100(x 0.25)2] and the points represent 7(x) calculated from the 0 and its derivatives according to equation (38). 3.2.3 Summary This method is only useful for twodimension problem of a thin strip and for threedimension problem of a thin plate since the instability increases dramatically with thickness a. Therefore, this method is inapplicable to the buried underground pipeline under cathodic protection system. In addition, the corrosion parameter 7 provides an inadequate description of the corrosion problem. These limitations show that TPA method is not suitable for this application. 3.3 TwoDimensional Boundary Element Method 3.3.1 Introduction Aoki et al.,18s19 solved a twodimensional inverse problem. The nonlinear rela tionship between the potential and current density was assumed to be known on the boundary of cathode as equation (39) and anode as equation (310). q =qo{exp( exp (0C } (39) 3PC ( (c ) I q = qoa{exp (O ) exp } (310) Both the 0potential and qcurrent density were unknown. Aoki applied a Newton Raphson iterative procedure and used singular value decomposition (SVD) method to get 0 and q. Where qoa, qoc, P,, aL, pa ca and ac are parameters. Except for the parameter to represents the connecting position of two electrodes, four param eters are used for each electrode. Thus, these parameters need to be fitted in the inverse process. Kenji and Aoki20'21 simplified equations (39) and (310) by using an inverse hyperbolic sine function. The polarization curves for anode and cathode turn out to be, f(q) = a1 sinhl q + P/ (311) f(q) = a2 sinh1 q + 2 (312) where a,, pi, a2, and 32 are the parameters. The number of parameters was re duced a factor of 2. The authors used a conjugate gradient method to minimize the objective function (313) k g(ai, ;i, xo) = (j j)2 (i = 1, 2) (313) j1 where (j and )j represent measured potential and calculated potential respec tively. xo is the connecting point of two electrodes. However, the results of the fit ted parameters and their standard deviations were not given. In order to improve the accuracy of the results, the authors checked the eigen values of the Hessian matrix of the objective function at the convergence point and found that three of the five eigen values were approaching zero. In other words, the five parameters of the objective function actually occupied about a twodimensional subdomain q=O or 4 =const q=O q=O Cathode Anode q=f(f) q=fa(4) ) guess ) guess Figure 33: Schematic representation of the twodimensional problem. instead of a fivedimensional domain. This is the reason that the results cannot be close to the exact solutions and why it is difficult to give the accuracy of the param eters. In addition, an a priori technique was used in their work. A priori means the probable values of some parameters that need to be estimated are known. In such cases, one may want to perform a fit that takes this advance information properly into account, neither completely freezing a parameter at a predetermined value nor completely leaving it to be determined by the data set. 3.3.2 Evaluation Process The procedure of evaluating the inverse model described by the literature took place in two stages. At the first stage, a forward model for the twodimensional system was created, which was similar to the work of Aoki et al. 's.l1820 Based on the forward model, the second stage involved testing the regression approaches and comparing them. 3.3.3 TwoDimensional Problem The twodimensional problem studied here is in a rectangular domain as shown in Figure 33, which is by described Aoki et al. .18,20 The lower horizontal bound ary consists of a cathode and an anode. The cathode may be considered to play the role of the pipeline, and the anode plays the role of a sacrificial anode. The upper 084 ** * Electrode 2 IElectrode Potential on Boundary Electrode Surface E Boundary > 088 .. I ! 1, t Potential at c 0  C 092 Ground Surface " 0 96 S*2. 2 1 00 I I I 00 02 04 06 08 10 00 02 04 06 08 10 Position / m Position Im (a) (b) Figure 34: Reproducing results for the twodimensional forward model: a) Poten tials at the boundary of domain; b) Current density along the electrode surface. horizontal boundary can be considered as the soil surface. Neither the potential nor current density q are known on the electrode boundary, but the nonlinear re lationship between them is known. The nonlinear relationship between potential Sand current density q is the polarization curve of the electrodes, which are given as equations (39) and (310). The parameters are qoc = qoa = 1Am 2, 0c = 0.845V, a = 0.985V, ac = 0.001V, a = 0.025V, )c = 0.05V and a = 0.05V. The parame ters a, 3 and qo are polarization parameters to be obtained by the inverse technique described in section 3.4. The subscripts c and a represent the cathode and anode, respectively. Potential within the interior of the system is governed by Laplace's equation. We followed Aoki's method, applying a NewtonRaphson iterative procedure for the forward problem. The potential and current density distributions from thefor ward model are shown in Figure 34. As seen in Figure 3.4(a), Aokiet al. showed that the variation of potential on the upper surface is much less than that seen on the electrode boundary.1820 It is evident that a large distance between the elec trodes and the soil surface will blur the measurement data on the soil surface. X measured potential j i=O in=O Cathode Anode in=fc() in=fa) Polarization curve is unknown Figure 35: Schematic representation of the twodimensional inverse model. 3.4 Inverse Analysis in Two Dimensions The purpose of this portion of effort was to apply an inverse boundary ele ment method to develop an efficient approach for identification of the polariza tion curve in a cathodic protection system. Different regression methods will be explored in this section. 3.4.1 Objective Function An objective function, which describes the difference between the measured potential and calculated potential on the soil surface, was constructed as following k g(ai, A;, xo)= (j __ j)2 (i = 1, 2) (314) j1 where (j and Oj represent measured and calculated soil surface potentials, respec tively. For the twodimensional inverse analysis, Aoki et al.18'20 described it in Figure 35. The parameters of the polarization curves were not known. What was known is the measured potential on the soil surface. The goal of the inverse analy sis is to identify the parameters of the polarization curves that would provide the minimum value for the objective function. 3.4.2 Regression Method Analysis We tested several regression methods, including the Le';ildergMA lirquart, conju gate gradient and quasiNewton methods. These methods, distinguished by the need to evaluate both the function value and the function gradient, were rejected for the pipeline inverse model. For the present problem, evaluation of the function gradi ent can only use finitedifference method, which requires extra calculations of the forward model. It inevitably adds significant work to the computation. Besides, the most serious deficiency was that these methods were found to be extremely sensitive to the value of an initial guess and to have a tendency to find local rather than global minima. The downhill simplex method requires evaluation only of the objective function itself, i.e., the derivative of the objective function is not needed. As the results of this work, it is shown that the downhill simplex method is the most robust one, independent of the measurement position and the initial guesses. It was successful for our twodimensional inverse analysis and will be further de veloped for the threedimensional problems. 3.4.3 Downhill Simplex Method The downhill simplex method was first created by Nelder and Mead.22 Only the objective function needs to be evaluated, no such need to its derivative. An objective function, which has N number of adjustable parameters, is evaluated for N + 1 points. The N + 1 points are generated by defining a starting point Po, then other points are obtained by using Pi = Po + Aiei (315) where ei are N unit vectors, Ai are constants set according to the problem's charac teristic length scale at different directions. The geometrical figure including N + 1 vertices, all the interconnecting line segments and polygonal faces is called a 'sim Reflection high val reflection and expansion high value contraction 4 high value / multiple contraction high value Figure 36: Schematic representation of the downhill simplex method. plex'. According to the function values on N + 1 points, by using steps, such as re flection, reflection and expansion, contraction and multiple contraction, as shown in Figure 36, one can find a new point to continue the function evaluation until the objective function encounters a minimum. 3.4.4 Accuracy of the Parameters One of the important issues of the regression process is how to determine the accuracy of the parameters. For a x2 merit function N 9i iYi2 x2(ak) = Y (316) i=1 i ak represents the parameters. gi and yi are the model and the measurement values respectively. ,6 is the standard deviation for the measurement data point yi. In this section, we are going to discuss the uncertainties of the parameters ak. The variance and the covariance of parameter ak are given as i Dai dak 0a2 ( ak j,k = 12,... M (317) ak,aj J o i )i If = k, Oa or 2 a j)2 (318) aj I Oyi 2ak is a M x M matrix. We can set 2ak = [C] (319) The gradient of X2 with respect to the parameters ak will be zero at the X2 mini mum. Thus, from equation (316), we get OX2 N [i y k = 1,2,..., M. (320) Dak i=1 I \i OakJ Taking an additional partial derivative gives 02 X2 N 09i i d. yi1 (321) =2 (9i Yi) (321) aijOak i=1 aj Oak (ajak 2 Note that the components of equation (321) depend both on the first derivatives and on the second derivatives of the basis functions with respect to their param eters. Inclusion of the second derivative term can in fact be destabilizing if the model fits badly. From this point on, we always define ajk as 1 02 2 i y 1 (322) 2 OaiOak i=1 aj k _ Comparing equation (319) with equation (322), we get [C] = [a1 (323) The diagonal elements of the matrix [a] 1 are the uncertainties of the parameters. These methods will be used to estimate the uncertainties of the parameters in the next section. Table 3.1: Regression result for twodimensional inverse model described by Aoki et al. The underlined symbols represent the parameters estimated by regression. Initial Parameters Final Parameters g Standard Deviations al = 0.2, al = 0.1, 1.32603 x 10 12 a = 1.594 x 10 6 1 = 0.7, 1 = 0.7, a2 = 1.0, a2 = 1.0, 32 =0.0, 32 =0.0, Connection = 0.5. Connection = 0.5. al = 0.1, al = 0.1, 1.32603 x 1012 aU, = 1.46 x 10 6 A3 = 0.1, 1 = 0.7, a2 = 1.0, a2 = 1.0, 32 =0.0, 32 =0.0, Connection = 0.5. Connection = 0.5. al = 0.2, al = 0.1576, 7.25157 x 104 oa = 3.16 x 10+ 1 = 0.7, i = 0.7, a2 = 1.3, a2 = 1.2112, 02 = 1.13 x 10+3 /32 =0.0, /32 = 0.0, Connection = 0.5. Connection = 0.5. al = 0.2, al = 0.164, 7.35439 x 103 aa = 1.18 x 10+2 1 = 0.5, /1 = 0.6632, p, = 3.16 x 10+1 a2 = 1.3, a2 = 1.0708, 02 = 1.26 x 10+3 /32 = 0.1, /32 = 0.497, 32 = 1.0 x 10+2 Connection = 0.5. Connection = 0.5. Uxo,,,,cion = 4.35 x 101 3.4.5 Regression Results The results of the regression for the twodimensional problem are presented in Table 3.1. The underlined symbols represent the parameters estimated by regres sion. It is clear that the number of independent parameters that can be obtained in a statistically significant manner is limited. When four parameters are fixed, the fifth could be obtained with high confidence. When two and more parameters were determined from the regression, the confidence interval became too large. This analysis provided two important insights. The first was that the number of statistically significant parameters that could be obtained by regression was limited by the amount and quality of the data. The second was that the weighted regression procedure proposed in this work could provide an indication of when the number of statistically significant parameters was exceeded. 3.5 Conclusions Different regression techniques, such as conjugate gradient method, quesiNewton method and downhill simplex method, were compared. As the results of the stud ies in this chapter, the downhill simplex approach was shown to be the most robust method and it is independent of the initial guesses. It will be applied further for the threedimension problem. However, the results from all the regression procedures were insufficiently re liable. One possible reason may be the insufficiency of data, since only one series of potential data has been used. Future work will incorporate not only the on potential data, but also the offpotential and the current survey data. The complex objective function to be minimized will be: 2(a am) (i yi)2' + n ( yi)2 L = i onpotential i i off potential S (y[ i yi)2] (324) =1 current where al,...,am are the parameters to be fitted, ci is the standard deviation of a measured value, and ni, n2, n3 are the number of data points. CHAPTER 4 DEVELOPMENT OF THREEDIMENSIONAL FORWARD AND INVERSE MODELS FOR PIPELINE WITH POTENTIAL SURVEY DATA ONLY 4.1 Introduction From the previous chapter, it is known that pipeline surveys are usually used to determine the adequacy of cathodic protection (CP) for the underground pipelines. However, these survey data may contain significant scatter or noise. The limita tions of measurement techniques, instruments, and analysis methods make the data interpretation difficult. Therefore, the objective of this chapter is to develop a mathematical framework for interpretation of survey data in the presence of ran dom noise. There are two types of approaches used in modeling cathodic protection. A forward model yields the distribution of current and potential for a given system geometry and for given physical properties of pipes, anodes, ground beds, and soil. An inverse model yields system properties such as pipe coating resistivity given values for current and/or potential distributions. The history of the development of analytic and numerical forward models for cathodic protection of pipelines encompasses more than fifty years. Waber et al.23 derived an analytic model in the form of a Fourier series, which can only be used for simple geometries and boundary conditions. Pierson et al.24 developed a se ries of analytic equations to account for attenuation for coated pipelines which extend the usual resistance formulas such as Dwight's equation.9 26 Doig et al.27 utilized the finite difference method to simulate the galvanic corrosion with com plex polarization. In the 1980's, Brebbia28 developed the revolutionary boundary element method (BEM). The BEM has since been applied in many engineering fields because it is accurate, effective for both infinite and semiinfinite domains, and computationally efficient.28'29 Aoki et al.ls820'30 applied the BEM to both two dimensional and threedimensional systems. Telles et al.31 improved the BEM for CP simulation by introducing the current density selfequilibrium limitation, which eliminated the need to discretize a boundary located at infinity. Brichau and Deconinck3'32 coupled the internal and external Laplace equations which were solved by coupled BEM and finite element methods (FEM). Kennelley et al.33'34 used FEM to model the influence on CP protection of discrete coating holidays that exposed bare steel. The concept of discrete holidays was extended to three dimensions using BEM by Orazem et al.35 37 Riemer and Orazem1'38 40 combined the advantages of the previous work in their development of the CP3D model. A distinguishing point of this software is that it can account for localized defects and yet is suitable for long pipelines and pipe networks. Some commercial programs are available such as PROCAT2 and OKAPPI.3 All the simulations described above solve Laplace's equation coupled with linear or nonlinear boundary conditions to obtain the potential distribution and current distribution. The analytic models can only be used for twodimensional domains with simple geometries. The numer ical models can be applied to almost any complex domain; thus, they are widely used. While many numerical models are based on the FEM, the newer generation of models use either BEM or a combination of BEM and FEM. The use of computer programs to interpret pipe survey data in terms of the un known condition of the pipe requires solution of the inverse problem, in which the properties of the pipe coating are inferred from measurements of current and po tential distributions. This approach allows interpretation of field data in a manner that would take into account the physical laws that constrain the flow of elec trical current from anode to pipe. Aoki et al. has made significant contributions to the application of inverse model to CP. Their studies included simplifying the unknown parameters in the polarization curves, changing the form of objective function, and trying different kinds of minimization methods. To minimize the ob jective function, they have employed a variety of techniques conjugate gradient, fuzzy a priori,20 and genetic algorithm methods.41'42 Their work involves mod eling protection of pipes, ships, and reinforcing steel in concrete structures.41 42 Wrobel and Miltiadou describe the application of genetic algorithms to inverse problems, including identification of coating holidays.44,45 Aoki et al. used an in verse BEM analysis to eliminate Ohmic error from measurement of polarization curves.46 The concept behind use of the inverse model for cathodic protection is related to the approach used to reconstruct the distribution of electrical resistance for med ical, chemical process, and geological applications. For medical applications, the readings from an array of sensing electrodes are used to construct an image as sociated with conductivity variations within a body.47 50 The technique is called electricalimpedance tomography.51'52 Electricalimpedance tomography has also been used to determine the interfacial area for twophase gasliquid and particu late flows in chemical processes53 56 and, through identifying the distribution of electrical resistivity, to identify composition distributions in laboratory and plant scale process equipment.57 Electricalresistivity tomography is used to interpret crossborehole resistivity measurements to obtain electrical resistivity distribu tions of geological formations.58 64Electrical resistivity is an important petroleum reservoir parameter because it is sensitive to porosity, type of pore fluid, and de gree of saturation. While use of neural networks have been suggested,65,66 most Figure 41: Schematic illustration of the pipe segment and anode used to test the inverse model. approaches are based on nonlinear regression. In this chapter, the construction of the forward model, inverse model for CP system, theoretical basis and development of the threedimensional BEM will be introduced. 4.2 Construction of the Forward Model 4.2.1 Cathodic Protection System For the threedimensional problem, a pipe was considered to be buried hori zontally under the ground surface. A cylindrical vertical anode was placed a dis tance away from the pipe. The underground region was considered to be bounded by the soil surface and to extend infinitely in the other directions in the soil. The anode was connected to the pipe with a wire, as shown in Figure 41. The pipe was placed horizontally 1.0 m (4.75 feet) below the soil surface and had a diame ter of 0.457 m (1.5 feet). The anode was placed in a vertical position and a distance away from the center of the pipe. The anode diameter was 0.152 m and the length was 1.0 m (3.28 feet). The soil resistivity was 100 Qm. A 0.5 mmthick coating was assumed to cover the side area of the pipe, with the exception that the two ends were assumed to be insulated. The potential of the pipe steel was assumed to be uniform. This assumption, valid for the short pipe segments considered in this preliminary work, will need to be extended in future work. The forward model was used to calculate the potential that one might measure with a voltmeter connected to the pipeline and a Cu/CuSO4 reference electrode (CSE). 4.2.2 Governing Equation and Boundary Conditions The current density in a dilute electrolyte solution can be described as8 i= F2V z z ujcj F zjDjVcj + Fvz zjci (41) where i is the current density vector, v is the bulk velocity, F is Faraday's constant, zi is the charge for species j, uj is the mobility for species j, cj is the concentration for species j, and Dj is the diffusion coefficient for species j. In this study, only the electrode surfaces, such as pipe walls are considered. Since the velocity v is small in the soil, the convection can be neglected. In addition, concentration gradients only have a significant effect in the diffusion layer, which is close to the struc ture and is relatively thin comparing with the characteristic length of the system. Thus, concentration gradients generally are neglected in largescale simulations. Therefore, the current density in the soil is i= F2V zYujcj (42) The conductivity of the electrolyte is defined to be S= F2V Z UjCj (43) The current density equation can be reduced to i = Vo (44) which is called Ohms Law. The conservation of charge in the bulk yields V i = 0 (45) Substituting equation(44) into equation (45) V (KVO)= 0 (46) Thus KV20 = 0 (47) Equation (47) is known as Laplace's equation for the electrochemical potential 0. The boundary conditions are shown in Figure 41. The corrosion potential of pipe steel is 0.6V. The sacrificial anode has a given potential of 1.2V. 4.2.3 Theoretical development As mentioned previously, the forward model has to be efficient in terms of com puter calculations due to its coupling with the inverse model. To achieve this goal, a linear relation, considering the coating covered pipe with various resistivity, is proposed in this section. Linear Boundary Conditions The polarization curve, which describes the relationship between the potential and the current density, indicates the corrosion condition on the pipe surface and is normally used as the boundary condition. It is not an easy task to determine the polarization curve since it strongly depends on a number of phenomena. Further more, the polarization curve can also be a function of time and history. The polar ization curve used by Riemer and Orazem1'39 had eight parameters. Amaya and Aoki et al.20'21 simplified the polarization curve from the ButlerVolmer equation and reduced the parameters to four for one cathode and one anode system. How ever, if the corrosion condition along the pipe is not uniform, more polarization curves are needed. In present study, the polarization curve was substituted by a linear relationship between the potential drop and the current density over the pipe coating for threedimensional simulation, shown in Figure 42. Coating Vouter Resistance SPipe Steel V Inner i= (1/ p5)[(V0'outer)(V, )inner)] Figure 42: Linear relationship between potential drop and current density over the pipe coating. 1 i = [(V Oouter) (V inner)] (48) PcO where i is the current density along the pipe radical direction, pc is the unit resis tance of coating material, 6 is the coating thickness, V is the potential of the pipe steel, assumed to be constant, ( is the potential in the electrolyte, V outer is the potential of the pipe referred to a reference electrode placed at the outside coating surface, and V inner is the potential of the pipe with respect to a reference elec trode placed at the inside coating surface. Under the assumption that the current densities through the coating and those at the pipe steel surface are equal, [(V Oouter) (V Oinner)] [(V inner) (V Ocorr)] ( Rc Rkinetic where Rc and Rkinetic represent the coating and polarization resistance, respec tively, and V Ocorr is the corrosion potential. As Rc > Rkinetic, (V Oinner) (V corr)] 0 (410) Therefore, 1 i = (V Oouter) (V Ocorr)] (411) PcO Equation (411) provides a linear relationship between the potential and the cur rent density. This linear boundary condition, replaces the nonlinear polarization curve and greatly simplifies the calculations of forward model. Description of the Pipe Surface According to the coatingcovered special feature of pipe, a coating resistance approximation was made. It allows the coating resistance to vary along the length of the pipe according to: P = po + pke (xxk)2/22 (412) k wherein po refers to the average good coating resistance of the pipe, pk represents the resistance reduction associated with coating defect k, xk is the center point of the coating defect, and Ok is the half width of the defect region. Equation (412) shows that if the coating position is close to a defect center xk, i.e., falling into the defect region 2ak, the coating resistance decay is evident; whereas, if the coating position is far away from the defect center, the coating resistance decay is insignif icant. The present approach has two significant advantages over assigning a resis tivity to each cylindrical element. For instance, use of equation (412) to represent two coating defects on a 150 feet pipeline is shown in Figure 43. If the pipeline is discretized with a 10 foot space interval, only seven parameters are required Po,P1,x1,a1,p2,x2 and 2 for the equation (412), as compared to 15 parameters are needed for an elementbyelement coating resistivity estimation. If the pipeline is discretized with a 5 foot spacing, still these seven parameters are required, while it is now thirty parameters for the regular method. Likewise, if the space interval is decreased to be so small that it makes the resistivity value continuous, no more than these seven parameters are enough to describe resistivity along the whole pipeline. The degree of freedom for the problem is increased dramatically, and correspondingly, the required computational time is reduced. 1200 . .. 1000 800 600 Model 60 10foot elements 400 5foot elements 200 0 50 100 150 x,ft Figure 43: Schematic illustration of the pipe surface resistivity model. 4.3 ThreeDimensional Boundary Element Method The theory for application of the boundary element method to the threedimensional halfinfinity multiconnected domain is introduced in this section. The linear cylin drical elements are adopted to discretize the pipe/anode surface in present study. 4.3.1 Infinity Domain In the twodimensional system, a closed rectangular domain was investigated. This kind of system, called an interior problem, is shown in Figure 44 and where S is the boundary of the volume V, the shaded area is the interior domain. f is the normal vector of the boundary S at that point. According to Brebbia et al. ,29 the basic equation of boundary element method corresponding to the Laplace's equation is c(xs)u(xs) = [q(y)u*(xs, y) q*(xs, y)u(y)]dS(y) (413) Equation (413) is called boundary integral equation. Wherein xs is the source point, y is the field point. u*(xs, y) and q*(xs, y) are the fundamental solutions of Green's function for potential and flux. u and q are the general solution of Laplace's function for potential and flux. c(xs) depends on the position of the n Figure 44: Interior problem. source point xs, as shown in equation (414), Xs X V, c(xs) = 0 c(x) = xs V, c(x) = 1 (414) Xs e S, c(xs) = 4 w is the solid angle subtended by V at xs, which is on the boundary S. For the threedimensional system, the studied domain is outside the pipeline and the an ode surface, and is between the soil at infinity. It is a multiconnected region of an interior problem, as shown in Figure 45. The boundary S can be considered to be the boundary of the pipeline or the anode. The boundary SR can be the boundary at infinity. When R oc, equation (413) is also valid for this case. For xs C S, c(xs)u(xs) = [q(y)u*(xs, y) q(xs, y)u(y)]dS(y) + lim [q(y)u*(xs, y) q(x, y)u(y)]dS(y) (415) Roo JSR The fundamental solutions to the Green's function for threedimensional problem at the boundary SR are 1 u*(x, y) = 47 R Ou*(xs, y) 1 q*(xs, y) = (416) On 47R2 0. SR n n Figure 45: Multi connected region of interior problem. And the general solution to the problem over the infinity boundary SR are 1 u(y) = Ku*(xs, y) + u, O( ) K 1 q(y) = 4 0( ) (417) 47R2 R2 where K and u, are two undetermined constants at the infinity boundary. K is the sum of source term (e.g., electric charge) distributed over S, K =j q(y)dS(y) 0(1) (418) Js Equation (416) and (417) are introduced into equation (415). Since limRo JSR [q(y)u*(xs, y) q*(xs, y)u(y)]dS(y) = limR fS J[q(y)u*(xs, y) q*(xs, y)(Ku*(xs, y) + uo)]dS(y) = limR,, Is {u*(xs, y)[q(y) Kq*(xs, y)] q*(xs, y)u,}dS(y) (419) where the first term inside the bracket is limRoo sR u*(xs, y)[q(y) Kq*(xs, y)]dS(y) = lim s s, 4 [O() ( )]dS(y) limn ,[O(T)) 47rR2]  0 (420) and the second term is f 1 lim q*(xs, y)udS(y) = lim u,( 47R 2) R4oo JSR Roo 47rR2 = u_ (421) Then equation (415) becomes, c(xs)u(xs) + [q(y)u*(xs y) + q*(xs, y)u(y)]dS(y) = u. (422) Equation (422) indicates that the value of K has no important effect on the bound ary S. To satisfy the pipeline and the anode CP system, the unit outward normal n on S in Figure 45 needs to change its direction. Correspondingly, equation (422) will be c(xs)u(xs) + {q(y)u*(x, y) + [q*(xs, y)u(y)]} + dS(y) = u (423) In addition, according to the physical condition, there is no current flow out of the soil surface, therefore, K = 0 is given as boundary condition and u. is unknown. 4.3.2 HalfInfinity Domain The halfinfinity domain is obtained by using a plane to split an infinite do main. Two half spaces lie on each side of the plane. Either the Dirichlet or Neu mann boundary condition is satisfied at the plane. For the buried pipelines sys tem, the underground soil domain is the halfspace being studied. Since there is no net current flowing out of the soil into the air, the Neumann boundary condi tion vanishes at the soil surface. The Green's function can satisfy that condition by using the reflection properties.29 70 If we set Q to be the region x < 0, P E Q and let P'(x, y, z) be the mirror image of a given source point P(x, y, z) with respect to the plane x = 0 (or 9OD), P' becomes the reflected source point, and r =  Q P, r' =  Q P' ( see Figure 46). Q is the field point. Consequently, halfspace fun damental solutions satisfying a homogeneous Neumann condition HN(P, Q) = 0 on O9 can be deduced by adequate superposition of freespace fundamental solu tions for sources at P and P', that is aGN(P, Q) G n(Q) = HN(, Q) = 0 (424) aQ Here, the superscript N is for the Newmann boundary conditions; n(Q) is the normal vector at point Q. Assuming two full freespace problems with VQ C Q , equation (424) can be written as GC (P, Q) 1 ( Q) n(Qr) = H(P, Q) = )n(Qr) (425) OQ 47r2 G+(P" Q) n(Qr,) = H(P', Q)= ( ) n(Qr') (426) OQ 4ir'2 From Figure 47, the normal vectors n(Qr) and n(Qr') of Fand r' at point Q are n(Qr) = 1 (427) n(Qr') = 1 (428) Meanwhile, according to the Neumann boundary condition given in equation (4 24) HN(, Q) = H(P, Q) + H(P', Q) = 0 (429) When Q e O9d, r = r', the fundamental solution in two full freespace have to be 1 G(P, Q) = (430) 47rr 1 G+(P', Q) = 1 (431) 47r' to satisfy equation (424), indicating that the two source intensities P and P' are equal and have the same sign. Therefore, the final form of the Green's function is GN(P, Q) = G(P, Q) + G (P', Q) = 1 + (432) 47 r r' ,t 1~ S\r 0\ R=Radium x Figure 46: Schematic illustration of mirror reflection technique. Figure 47: The Neuman b.c's. Fundamental solution to the halfinfinity domain satisfying the x eo eR 0 1 2 1 3 3..... n2 n1 n n+1  .    ', ,'  .   Sne i i nel+ S0 S2 S3 Snel2 Snel1 Figure 48: Pipe discretision and collocation points. 4.3.3 Boundary Discretization Rectangular and triangular elements are commonly used in the literature for BEM meshes on the surface of a threedimensional object. Aoki et al.71 used 352 constant triangular elements to discretize the side area of a 160mm long, 155mm diameter cylinder. Riemer and Orazem's model CP3D may involve thousands of elements to discretize a pipeline with discrete coating holidays. This level of de tail is appropriate for a forward model, but it is inappropriate for an inverse model where the amount of data may be insufficient to extract a high level of details con cerning the pipe condition. Besides, the model needs to be calculated not once, but thousands of times. To decrease the number of elements, a series of linear cylindrical elements were used to take advantage of the cylindrical shape of the pipe and the anode. The source points (also called collocation points) were local ized on the cylindrical side area. At the two ends of the objectpipe/anode, the two discs were set as two constant circle elements. The center of each circle was a collocation point (see Figure 48). For the cylindrical elements, there was no varia tion along the circumference. The variation was only along the axial direction72 75 such that the collocation points could be chosen at the reference line, which was on the top of the cylinder. The elements on the side area of the cylinder were from So to (Sne11), and, at the two ends of the pipe, were the two circle areas Snei and S(ne, 1). The collocation points on the side area of cylinder were from 0 to (n 1), while two collocation points at the center of circles are counted as n and (n + 1). The relation between the number of elements and the number of nodes depends on the order of the interpolation functions, which can be constant, firstorder or even higher order. The identical mesh was used for the element separation because in the follow ing inverse model, the position of the coating defects was distributed randomly along the pipeline. If the pipeline meshing depended on the position of coating defects, every time when the new coating defect positions would be assumed, the system would have to be remeshed. To save computational time, by using the identical mesh, the matrices G and H for the BEM need to be calculated only once. 4.3.4 Coordinates Definition and Transformation In this program, two kinds of objects were defined, one is pipe, and the other is anode. The user needs to input the beginning point of the pipe/anode's axle to the program. Both the global coordinate and the local coordinate were used in this study. The global coordinates was set as: z axial parallel to the soil surface and directed to north direction; x axial perpendicular to the soil surface and pointed toward the sky. The y axial and z axial build the soil surface plane. The origin is on the soil surface. The angle between the pipe/anode's axle and soil surface is defined as the elevation angle ,O. e, = 0 means that the object is parallel to the soil surface. 8O = 90 (or e, = 90) means that the object is perpendicular to the soil surface. The direction angle Od refers the angle between the pipe/anode's axle and z axial positive direction. The local coordinate is used when the collocation point and the field point are on the same object (i.e., the pipe/anode), while the global coordinate is used when the collocation point and field point are on the different objects to clarify the relative positions of these objects. Transformation between two different coordinates, such as the global and the y Figure 49: Coordinate rotation. local coordinates or two different local coordinates, usually consists of two steps. One is origin transformation, which is to move from one original point to a new origin. The other is rotation, which is to rotate previous axials around the new origin. Coordinate rotation actually includes two rotations. Figure 49 shows the procedure to rotate a coordinate around a point and one axis. The first rotation is to rotate the x axis and the z axis about the y axis with elevation angle 0e to get new x' axis and z' axis. The second rotation is to fix x', and to rotate y and z' axes with direction angle Od to get a new axes y" and z". Finally, coordinates (x, y, z) is converted to (x', y", z"). 4.3.5 Discretization of Boundary Element Method Rewriting equation (423) using P as source point, Q as field point and O0Q as boundary, we get cu(P) + u(Q)qod(O0) = q(Q)upod(O) + u (433) I P where u* and q* are the fundamental solution of the Green function, and they have the relation OU*po q (434) Discretizing equation (433), nel f u* nel ciui(P) + uj(Q) dS, = I qj(Q)Ou*odSj + u, (435) jl Jsj Ono j1 S, where nel means the number of element. Since u and q of each linear cylindrical element vary linearly along the z axial direction, then U(z) = i(l z) +2Z q(z) = q(1 z)+q2z (436) where 1 is the length of the element. Let t = ', then u and q, the functions of z, can be transferred to be the functions of t, u(t) = l[u(l t) u2t] = l[lNi(t) + u2N2(t) q(t) = l[q(1 t) + q2t] = l[qiNi(t) + q2N2(t)] (437) The basis functions can be extracted from equation (437), Ni(t) = t; N2(t) = t. (438) Upon substituting equation (437) into equation (435), the discretized boundary integral equation becomes nel &u* ciui(P) + l lj[uj,Nil(t) + uj,2N2(t)] 1jdSj j1 Sj Onj nel = lj[qj,lNl(t) + qj,2N2(t)]uidSj + um (439) j=1 j Since the fundamental solutions are 1 S= 47rij qi* O = 4 j (440) zy 9ni 47r 9O{ni According to equation (439), we define Gij = lqjkN dS (441) Jik=l Hij = lukNktdS (442) Jik=l If the field point Q is on the cylindrical element side area, dSj = RdOdt, R is the radium of cylinder, then Gij becomes j l J27 2 2 Gi = 1 lqjkNj(t)u Rd0dt k=1 1R 1 2 22 k 1 = li ) (lkNk(t))dt ii dO (443) and Hij becomes H11 l j*27 2T Hii = ljujkNk(t) qiRdOdt k=1 1 2 27 0(1) = l )(4) (I ukNk(t))dt j dO (444) 47 Jo k1 JO Onj If the field point Q is on the circle elements, dSj = rdrdO, 0 < r < R, Gij changes its form to 1 jR j27, Gij = qj(4i) rdrdO (445) 47 JO Jo rij and Hii is in form 1 R 2 (1) Hij = uij() j / rdrdO (446) 47r Jo Jo 9n The boundary integral for Gij and Hii with the different positions of the source point P and the field point Q is shown in Appendix A. In addition, adding the unknown potential at infinity u, in equation (439) to Hij, H matrix becomes H= p H, 1 (447) Ha,p Ha,a,, 1 [ H (np+na)x(np+n,+l1) where 1 column corresponds to u,. And the G matrix is G= GP Gy,a (448) Ga, p Ga,a (np+na)x(np+na) The subscript p and a indicate the pipe and the anode respectively. The subscript (p, p) and (a, a) indicate that the collocation points and the field points are on the same object. While (p, a) means that the collocation points are on the pipe and the field points are on the anode. Likewise (a, p) means that the collocation points are on the anode and the field points are on the pipe. The number of node on the pipe and the anode are represented by n, and na. 4.3.6 Row Sum Elimination The diagonal elements of the matrix H are usually difficult to calculat since the linear or higher order elements have been used and the constant ci in equation (414) involves the calculation of the solid angle subtended by the region V at xs on S. In order to overcome this difficulty, Gibbs's theorem can be applied to find the values of the diagonal.76 For a domain governed by Laplace's equation, if the potential is uniform throughout the domain, and its gradient of the potential is zero at infinity if it exists, the gradient of the potential will be equal to zero within the domain. For the threedimensional potential problem, u(y) 1 is defined in the interior region. Since q(y) = aY = 0 on the boundary, q(y) = 0 in the domain. Thus equation (413) becomes c(s)= q*(x, y)dS(y) (449) JIs It means that, YHi,i = 0 (450) ij where the diagonal terms are unknown. The diagonal terms are now specified by the negative sum of the offdiagonal terms of each row H,i = Hi,i (451) j7i where j goes from one to the number of terms in a row. The same technique used for the interior domain is applied to present work. Assuming u(y) 1, since q(y) = (y) = 0 on the boundary and in the domain, On equation (415) becomes c(x) = *(xs, y)dS(y) lim q(xs, y)dS(y) (452) Js Rao JSR Additionally, since lim q*(xs,y)dS(y)= ( 1)(4R2) = 1 (453) Roo JSR 47r2i c(xs) in equation (452) becomes c(xs) = q(xs, y)dS(y) + 1 (454) which means that, H, = 1 Hi, (455) ii Therefore, by using Gibbs's theorem, the diagonal elements of the matrix H can be obtained from the offdiagonal elements and the calculation difficulties are avoided. 4.3.7 SelfEquilibrium Cathodic protection systems do not lose or gain current from their surround ings. In other words, the current is conserved. Therefore, if the flux between the anode and the pipe is selfequilibrium, the flux at infinity will be zero. To imple ment this condition, an extra equation, equation (456) is added to the system.31'77 i qdS =0 (456) It indicates the intensity of an equivalent source distributed over S. It also can be written as i' = ip,kAp,k + Y ia,kAa,k = 0 (457) k k where Ak= dSk, k c [0, np] or [0, n] (458) Adding equation (457) to equation (447) and (448), and defining the potential vector u as Op = (459) 0. (np+na)xl the flux vector q as n Vp = n (pn(460) S(np+na)xl we get Hp,p Hp,, 1 p GP,p Gya Ha,p Ha,a 1 = Ga,p Ga,a (461) n V^a 0 0 0 0, Ap Aa V Here H matrix is singular because there is still a row of zeros in it. However, it will only be the case when Neumann type boundary conditions are specified everywhere. The Neumann problem results in an infinite number of solutions. Therefore, at least one element in the system must have a Dirichlet boundary con dition to make the H matrix nonsingular and result in a unique solution. In addi tion n VO has a minus sign because the normal vector direction is outward of boundary S. 4.4 Forward Model 4.4.1 Constant Steel Potential Assumption Since the axial variation of potential in pipe steel is trivial for a short section or low resistance pipelines, the potential of pipe steel is assumed to be a constant in this chapter. Rewriting equation (461) by making a variable substitution = V 9, we get H,,, H,,, 1 V Gpp Gp,a G V H,,p H,,, 1 V ,', = G, Ga, n (462) n (V '.,) 0 0 0 V o A, Aa where V is the potential of pipe steel. Decomposing the vector V 0 and n * V(V ), equation (462) becomes H,,,p H,,a 1 V H,,,p H,,,a 1 H,,, H,, 1 V H,,, H,,, 1 0 00 V 0 0 0 Gpp GPa = G,,p Ga, n (463) n 7 ... AP A Since the potential of steel V is assumed to be a constant, n VV = 0. Moreover, according to equation (455) Thus, the first term can readily obtain on the left hand side of equation (463) is vanished, and we Hp,p Hp,a 1 Ha,p Ha,a 1 ', 0 0 0 Equation (411) can be rewrite as Gp,p Gp,a Ga,p Ga,a Ap Aa n Psoil cor V ',' o (',= '1' corr) PcoatO Substituting equation(467) into equation (466), we can get Hp,p Hp,a 1 GY Gp,a Ha,p Ha,a 1 a Ga, p Ga,a I 0 0 0 Aoo Ap Aa Equation (468) will be applied in the two cases as below. Sacrificial anode case (466) For the sacrificial anode protection method, we set the potential of the anode as a constant and the current density of the anode is unknown except that it is zero at the two ends of the anode. The potential of the pipe is unknown. The current density of pipe is also unknown except the zero value is given at the ends of pipe. In addition, the potential and the current density of the pipe satisfy the relation of equation (467). H,p + Hp,a Ha,p + Ha,a,, (464) (465) (467) ,''. Vcorr) (468) .V,7 Rearranging the equation (468), HP'P + PPpco Gp,a 11 "1. PPpcoilS Hpa H,p+ G,psoi Ga,a 1 n V' PsoiG Ha,a S~ap + Pcoat6 P1 n a Pcoat6 +A Psoi Aa 0 J A Pil 0 P coatO PPcoat (469) it can be used to model sacrificial anode CP system. Impressed current case To simulate the impressed current protection method, we set the current den sity of the anode as a constant value, but the potential of the anode is unknown. The conditions of current density and potential of pipe are same as those of the sacrificial anode case. Likewise, the following equation can be obtained. Hp,p Gpp Psoi Hp,a 1 G'1' PPpcoS Gp,a H, + G Psoil H, 1 { G= Psoil Ga,a +AP 0 0 Ap so Aa P PcoatO PcoatO (470) 4.4.2 Simulation Results Typical simulation results for potentials on soil surface directly above a 10m pipeline with one coating defect are shown in Figure 410. The potential change at the position corresponding to the coating defects is evident, but the magnitude of the change is small, in agreement with previous calculations.35'38 4.5 Inverse Model The purpose of the inverse model is to construct an objective function and to compare different types of regression methods in order to minimize the objective function. The simulated annealing method was found appropriate in the present study since it can escape from the local minima and thereby avoid the influence of the initial guess. Potential on the Soil Surface 1.109 1.11 1.112 o_. 1 0 2 4 6 8 10 Axial Position/m 1 11 1.114 1.115 Figure 410: False color image of the onpotential on the soil surface that was gen erated by the forward model corresponding to Figure 41. 4.5.1 Objective Function An objective function, which describes the difference between the measured potential and calculated potential on the soil surface, is given as N g(x,p, ) = j(Vj j)2 (471) j1 where 4j and 4j represent measured and calculated soil surface potentials, re spectively. The calculated soil surface potential depends on the pipe coating con ditions, such as position, resistivity and width of coating defect. Thus, once the minimum of the objective function is reached, the fitted coating parameters will reflect the real physical conditions of pipeline. 4.5.2 Analysis of Regression Methods Several regression methods were tested, including the conjugate gradient, Levenberg Marquart, and quasiNewton methods. A summary of these methods is provided by Press et al.78 These methods, distinguished by the need to evaluate both the function value and the function gradient, were rejected for the pipeline inverse model. The extra calculations required for evaluation of the function gradient added significantly to the computational requirement. In addition, these methods were extremely sensitive to the value of an initial guess and had a tendency to find local rather than global minima. The downhill simplex method,78 which requires evaluation of only the objec tive function and not its gradient, was used successfully for preliminary two dimensional inverse analysis. The downhill simplex method did not work very well for threedimensional problems because it was strongly affected by the initial guess. Genetic minimization algorithms, intended to mimic the course of natural se lection, have been applied to inverse problems for cathodic protection.42'44,45 Pa rameters are normally coded as binary strings to reduce the searching population. Procedures of selection, mating and mutation are used repeatedly to create the new generation until the specified stop criterion is satisfied. This method is not very sensitive to the values of initial guesses and can escape from the local min ima. However, it has difficulty selecting between close function values.79 After trying the alternatives, the simulated annealing optimization approach was selected for the present work. This method is attractive because it is suit able for largescale problems and can search for a global minimum, which may be hidden among many local minima. Simulated annealing was better than down hill simplex because the simplex method accepts only downhill steps during its searching; whereas, simulated annealing can accept both the downhill and uphill steps. In this way, simulated annealing method can step out of the local optima and successfully locate the global minimum. 4.5.3 Simulated Annealing Method The term "simulated annealing" comes from a physical process analogy.78 When a material is heated and then is slowly cooled down, a strong crystalline structure will be obtained. This crystal is the minimum energy state of the system. In a sim ulation process, a minimum of the objective function corresponds to the ground state of the substance. The Bolzmann probability distribution, which describes the different energy state in a thermal equilibrium under given temperature, can also be used to mimic the different objective function value at certain searching region. The simulated annealing method starts with a brief view of the searching domain by making large moves and then it focuses on the most possible region. The inher ent random fluctuation permits the annealing procedure avoid the local minima. The Metropolis rule p = e (g2)/kT (472) is used to control whether the uphill step is accepted. In equation (472), p refers to the probability, gi is the objective function value and T is an important param eter in the simulated annealing method, which resembles the temperature in the thermal system. For g2 < gi, the probability p is greater than 1, which means that the downhill steps are always accepted. For g2 > gi, the probability p is compared with a uniformly distributed random number from [0, 1] to decide whether uphill steps are accepted. Since each parameter has its own limit, once the parameter is out of its limit, the following equation can be used to correct the parameter. x = XL + a(xu XL) (473) where a is a uniformly distributed random number among [0, 1], XL is the lower limit of x and xu is the upper limit of x. The parameter x can thus be guaranteed to lie within its bounds. 4.5.4 Simulation Results and Discussions Examples of the application of the proposed forward and inverse models are given in this section. In addition, several inverse strategies are introduced. The re Figure 411: Flow chart for the inverse model calculations. gression algorithm employed by the inverse model is shown in the Figure 411. The set parameters were used as inputs to the forward model to get the potential on soil surface. These results were treated as the measured data. These "measured" po tentials belong to each point of a 3 x 101grid area, as shown in Figure 412. There are three lines in the grid area: the middle line is right above the pipe centerline, and two other lines of calculated values were placed +1 m from the centerline of the pipe. The grid, therefore, comprised 303 data points. Initialize the parameters, the number of defects, each defect's center position, resistivity reduction and the width of the defects, and input them into the inverse model. Among all these pa rameters, the number of defects has the most significant impact on the regression results. In this study, we purpose a method to decide the number of defects which can be obtained from the measurements. The inverse analysis result by using sim ulated annealing method is shown in Table 4.1 for a 10meter long pipe segment x=O,y=1,Z=0~10 I x=0,y=1,Z=O0lO              S10m 0 100 Figure 412: Grid showing the location of 303 surface onpotentials calculated us ing the threedimensional forward model developed in the present work. The grid shown is for a 10m pipe segment. A scaled version of the grid was used for a 100m pipe segment. with two coating defects. The calculated data were for defects located at the 3.2 and 6.0 m positions with characteristic dimensions of 0.05 and 0.5 m, respectively. The intact coating resistivity had a value of 5.0 x 107 / Qm, which was fixed for the regression procedure. The input values for surface potentials were calculated at the limits of the precision of the forward model, i.e., no noise was added. The regression yielded the correct locations, characteristic defect dimensions, and re sistivity changes associated with the coating defects. The process of finding the minimum of the object function is shown in Figure 413. For this test problem, with the forward model and no noise added to the syn thetic data, the minimum value of the cost function was 10 15. Of the techniques used, only the simulated annealing method could find this global minimum. The best of the other techniques were able to identify local minima with values on the order of 104 to 10 8 preceding the plateau. A comparison of the set and fitted results for pipe coating resistivity, poten tial and current density, respectively, is shown in Figures 414. Corresponding to the coating resistivity decay, the potential and current density distributions have significant changes. A good agreement between the synthetic data and regression results can be observed. Table 4.1: Parameter values obtained using the threedimensional inverse model developed in the present work for a 10m pipe segment with two coating defects. Coating Position Resistivity Characteristic dimension Defect Xk / m pk / Q m 0k / m Set Values 1 3.2 4.0 x 107 0.22 2 6.0 3.5 x 107 0.71 Initial Values 1 2.5 3.0 x 107 0.1 2 5.0 3.0 x 107 0.1 Regression Result 1 3.206 3.81 x 107 0.26 2 6.004 3.51 x 107 0.70 102 0 > 106 U 0) 1o' 0 1014 1018 0 50 100 150 200 250 300 3 # Evaluations Figure 413: The regression objective function as a function of the number of evalu ations for a pipe coating with one defect region. The simulated annealing method was used for this regression. I I I I I ' I I Initial Guess \ Evaluation of Local Minima Best Minimum Identified I I I I I 5x109 cm      Set data 0 Regression result 0 2 4 6 8 10 Position /m 1 148  Set data o Regression result 1 152 1 156 1 160 I I I 4 6 Position /m 8 10 0 E < lx10l (3) C c S2X10l 3x1 3x104 Set data O Regression result S 2 4 6 8 10 Position Im (c) Figure 414: Comparison of the set and fitted results for pipe coating with two coating defects: a) coating resistivity; b) potentials; and c) current density. E a 50 S40 D 30 0) S20 8 10 10 0 . . 1 :il 1 Table 4.2: Test case parameters with five coating defects on the pipe used to demonstrate the method for determination of the number of statistically signif icant parameters (see Figure 415). The intact coating resistivity had a value of 5.0 x 107 / Qm. Coating Defect Position Resistivity Characteristic dimension Xk / m pk / Qm Ok / m Set Values 1 10 3.5 x 107 0.45 2 30 4.5 x 107 0.32 3 40 3.0 x 106 0.55 4 70 2.0 x 107 0.45 5 85 4.0 x 107 0.63 4.5.5 Inverse Strategies The strategies, including the determination of the number of significant pa rameters, the procedures used to reduce the influence of initial guess and the eval uation of the role of noise in the measured data, were explored to assess the confi dence level of the inverse model results. Determination of the Number of Statistically Significant Parameters When the collected data are scattered, it is difficult to decide how many pos sible defects or coating anomalies should be included in the regression. The ap proach taken to address this issue was to increase sequentially the number of de fects, using the regression statistics to determine when the number of statistically significant parameters was exceeded. This approach is similar to that used to as sess the number of statistically significant measurement model parameters can be obtained by regression to impedance spectroscopy data.808,s An example with five defects, shown in Table 4.2, was used to illustrate the approach taken to assess the correct number of statistically significant parameters. The relation between the regression statistic log(x2/v) and the number of coating defects assumed in the regression is presented in Figure 415. Here, v = N M represents the degree of freedom of the problem which is reduced as the number of I '' I I I 7.2 o S7.6 _o 0 \ /0 oT / 0 8.0 \ / \ / 8.4 o 0 2 4 6 8 10 # defects Figure 415: The regression statistic as a function of the number of coating de fects assumed for the model. The minimum in this value is used to identify the maximum number of coating defects that can be justified on statistical grounds. parameters is increased. N represents the number of data point, and M represents the number of parameters. x2 is the weighted objective function. The lowest point of the curve denotes the number of statistically significant coating defects obtain able by regression. Here, it is four. The fifth coating anomaly, identified in Table 4.2 as defect 3, could not be identified by the regression procedure. The number of coating anomalies identified by this procedure will depend on the amount and quality of data and on the sensitivity of the data to coating condition. Testing the Robustness of the Inverse Model To assess the influence of uncertainty in the measured data, normally dis tributed stochastic errors were added to the soil surface potential generated by the forward model. The noise added had standard deviations Onoise of 0.1, 1.0, and 2.0 mV respectively. The set values, initial guesses, and regression results are presented in Table 4.3. Three coating defects placed at 30, 40, and 70 m positions were used to generate synthetic data. For the lowest noise level, noise = 0.1 mV, the sequential procedure described in the previous section allowed only two statistically significant defects. 64 Table 4.3: Regression results from the threedimensional inverse model for a 100m pipe segment with three coating defects. The sequential regression procedure was used to identify the number of defects that were statistically significant. The intact coating resistivity had a value of 5.0 x 107 / Qm. Coating Position Resistivity Characteristic dimension Defect Xk / m pk / Q m 0k / m Set Values 1 30 4.5 x 107 0.32 2 40 3.0 x 106 0.55 3 70 2.0 x 107 0.45 Initial Guess 1 50 3.5 x 107 0.92 2 50 3.5 x 106 0.92 3 50 3.5 x 106 0.92 Regression Results noise = 0.1 mV 1 29.85 4.25 x 107 0.66 2 69.83 7.57 x 106 1.30 Noise = 1.0 mV 1 32.01 4.67 x 107 0.087 Noise = 1.0 mV 1 32.01 4.68 x 107 0.057 2 90.29 1.22 x 106 0.86 Noise = 2.0 mV 1 29.95 3.33 x 107 1.50 The initial guesses placed the defects at the midpoint of the pipe segment. The re gression results suggested that the defects exist near the 30 and 70 m locations. The missing defect is the one with the smallest deviation in coating resistivity. Thus, the regression procedure identified the correct location of the two most significant reductions in coating resistivity. For noise = 1.0 mV, the sequential procedure using the minimization of the Sx2/V criterion allowed two coating defects. The defect located at 32.01 m was consistent with the most significant defect located at 30 m, but the defect identified at 90 m did not conform to the input data. In addition, the regression failed to identify the second most significant reduction in coating resistivity at 70 m. The problem here may be an inadequate sensitivity of the X2/V statistic for identifying overfitting of data. Other criteria, such as the Akaiki information cri teria,82'83,84 provide additional penalties for adding parameters to a model. For noise = 1.0 mV, the Akaiki performance index Ap = X21 ob (474) 1 np/nob suggested that only one defect could be identified; whereas, the Akaiki informa tion criterion Aic = log (I X2 (1 + 2np/nob)) (475) suggests that two defects could be identified. Regression for a single coating defect identified a defect in the close vicinity of the most significant coating reduction. By any of the statistical measures tested, only the largest defect could be identified for noise = 2.0 mV. The location of the defect was correctly identified; however, the width of the defect was incorrectly determined. These results can be explained by examination of the synthetic surface poten tial data used for the inverse model analysis. The corresponding results of the inverse analysis for noise = 0.1 mV are shown in Figure 416 As is shown in Fig ure 4.16(a), the level of added noise did not obscure the surfacepotential features introduced by the presence of the major coating defects. The regressed and noise free target values for potential cannot be distinguished; thus, the absence of the minor coating defect did not influence the fit of the model to the synthetic data. As seen in Figure 4.16(b), the regression procedure identified the the two most significant reductions in coating resistivity at 30 and 70 m locations. In contrast, the random noise added with noise = 1.0 mV, shown in Figure 4.17(a), almost completely obscures the influence of the coating defects. Neverthe less, a major defect can be resolved by the regression procedure in the vicinity of the largest defect, as shown in Figure 4.17(b). The anomalous defect introduced by the regression procedure at 90 m has associated with it a small reduction on coat ing resistivity. There is a question as to whether this defect can be considered to 66 0.988 U 0.990 0 0.992 0.994 . 0 20 40 60 80 100 Position /m (a) 60 E C 50 40 3 30 c n S20 10 I I I lIl 0 20 40 60 80 100 Position /m (b) Figure 416: Comparison between the input values and regression results for noise =0.1 mV: a) soilsurface potential (at the centerline of the pipe); b) coating resistivity. 0.986 0 ' 0 00 0 > 0.988 Q 0 0 00 000 o 00o o000 C 0.990 0 OO 0 o 0 0Q 0.992 0.994 I , I , I , 0 20 40 60 Position /m 0 20 40 60 Position /m 80 100 80 100 Figure 417: Comparison between the input values and regression results for noise =1.0 mV: a) soilsurface potential (at the centerline of the pipe); b) coating resistivity. be statistically significant, but, nevertheless, the regression procedure would have given adequate guidance for excavation of the pipe. By visual inspection, the random noise added with noise = 2.0 mV, shown in Figure 4.18(a), would appear to obscure the influence of even the major coating defect. The regression procedure identified a single statistically significant coating defect located near the 30 m defect. The agreement between the regressed poten tial profile and the noisefree target values was surprisingly good. As shown in Figure 4.17(b), the correct location for the principal defect was correctly identified, even though the breadth of the defect was not correctly determined. The results suggest that, while it is difficult to extract coating conditions from the data when the influence of the coating defect on the surface potential is compa rable with noise, the regression procedure showed a surprising ability to identify the location of the most significant coating defects from noisy data. This result suggests that an inverse model is feasible, in particular when other types of data are included. Sensitivity to Initial Guess The values used for an initial guess have a significant influence on most of the conventional nonlinear regression methods. To reduce this effect and to test the robustness of the present inverse model, identical parameters for each defect were used. For example, as the initial guess in Table 4.3, the positions of the defects were all set to the middle of the pipeline. This kind of initial guess created significant difficulties for conventional regression methods, but posed no problems for the simulated annealing method used in the present inverse model. 69 0.984 c i p  0 i0000    i  0.984 ' 0 0 of 0 00 > 0 00 0 0 0.988 qp 0 Oo 0O00 n ^  t3 O OO, O O OOOOO 0 0 0 O0 0.996  0 20 40 60 80 100 Position /m (a) 60  E 50o S40 > 30  S20  10 10  0 20 40 60 80 100 Position /m (b) Figure 418: Comparison between the input values and regression results for Noise =2.0 mV: a) soilsurface potential (at the centerline of the pipe); b) coating resistivity. 4.6 Conclusions In this chapter, both the forward model and the inverse model for the three dimensional CP system was introduced. The simplified forward model has the following advantages. 1. The Laplace's equation with the simplified linear boundary condition was applied. This approach not only simplifies the forward calculation, but also reduces the computational effort to the inverse calculation. 2. The line shape approximation for the pipe coating resistivity makes it possi ble to limit the number of parameters needed to describe the coating condi tions along the whole pipe. It dramatically increases the degree of freedom of the problem as compared to the method that represents the coating condi tions element by element. 3. The BEM method with cylindrical element was developed. The basic theory of BEM, such as use of the mirror reflection technique for the halfinfinity domain, row sum elimination, and selfequilibrium, were introduced. In ad dition, the implementation techniques, for example, the coordinates trans formation, element discritization, were also included. 4. The simplified forward model was used to obtain the potentials of the pipe on the soil surface. The potential and the current density of the pipe and the anode could also be calculated. The inverse part of work represents an ambitious research effort with a poten tially large payoff for the oil and gas transmission industry. The central question is whether a regression approach could be used to assess pipeline coating conditions from field data in a way that is consistent with the laws of physics. The rationale, methodology, and results are presented, which demonstrate the feasibility of the inverse model for pipelines. 1. This work has demonstrated that it is possible to couple a boundary ele ment forward model with a nonlinear regression algorithm to obtain pipe sur face properties from measured soilsurface potentials. The resulting model is called an inverse model. 2. The technique identifies the location of coating anomalies as well as the breadth of the anomaly and the amount that the local resistivity has changed. 3. The performance of the inverse model is sensitive to the regression procedure used. The simulated annealing algorithm proved to be the most robust and had greatest capability to seek the global minimum for this problem. 4. An algorithm was developed that could be used to identify the maximum number of coating anomalies that can be detected. This number is sensitive to the quality of data as well as to the actual coating condition. 5. If the number of coating anomalies detected is smaller than the actual num ber of coating defects, the technique will identify the most serious anomalies. The results of this effort, limited to a single pipe in a rightofway, demonstrate the feasibility of a program to interpret survey data in terms of the state of protection of the pipe. CHAPTER 5 DEVELOPMENT OF THREEDIMENSIONAL FORWARD AND INVERSE MODELS FOR PIPELINE WITH BOTH POTENTIAL AND CURRENT DENSITY SURVEY DATA 5.1 Forward Model 5.1.1 Introduction In Riemer's work, there were two separate domains for the flow of current. The first was the soil domain up to the surfaces of the pipes and anodes. The boundary conditions for the soil domain were the kinetics of the corrosion reactions. The second domain was the internal pipe metal, anode metal and connecting wires for the return path of the cathodic protection current.38 The BEM method was used for the first domain. The FEM method was applied to the second domain. Brichau also coupled the BEM and FEM to solve the two domains.72 In this chapter, to simplify the calculations, not only the potential but also the current density along the pipe steel will be obtained by using the BEM method. 5.1.2 Pipeline with Varying Steel Potential In chapter 4, the potential V on the pipe steel was assumed to be a constant. Actually, the potential on the pipe steel varies along the pipe since the pipeline has its own internal resistance. For a short section of pipe, since the potential drop due to the pipe resistance is very small, it is reasonable to assumed it to be constant. However, for a long pipeline or a high resistance pipe, the potential drop along the pipe steel is significant and cannot be ignored. The model in chapter 4 was extended to account for the electrical potential of the pipeline steel. i0 in I q  __ y Pipeline +iqa Lc2 L2 r: Anode rim a Figure 51: Cathodic protection system with variant potential along the pipe and anode. The method of calculation is illustrated schematically in Figure 51, where a horizontally placed pipeline is connected by a wire to a vertically placed sacrifice anode. The connection points are cl on the pipe and c2 on the anode. In Figure 51, i represents the current density along the axial direction of pipeline or anode, and subscripts p and a designate pipeline and anode, respectively. The current density entering the pipe at the ends is designated by i and i;, and the current density entering the anode at the ends is designated by i, and i. The current density entering the pipe coating in the radial direction is given by q. As shown in Figure 51, the protecting current flows away from the anode to the soil, then flows to the pipeline. On the pipeline, the currents flow from the two opposite directions to the connection point cl and they are concurrent as current I. The current I flows back to the anode through the wire. The current density in the pipeline steel along axial direction is given by 1 1 dV i = ( n VV) 1 dV (51) Psteel Psteel dz where V represents the potential of the steel in the pipe. Equation (51) can be Vk1 Vk As Figure 52: Relation of potential V and current density. integrated in the axial direction along a segment with length dz, shown in Figure 52, such that V1 dV= psteelZk idz (52) The potential change across the segment is given by Vk Vk = Psteel k ( zk k1 1 Zk )dz (53) Jkl Zk z1 Z Zk 1 or (ik1 + k) Vk Vk1 = Psteel(Zk Zk1) k (54) where the potential difference between two points is found as the product of the average current density (ik1 + ik)/2 and the resistance psteel(Zk Zk). Conservation of charge requires that the current flowing into the pipe segment is equal to the current flowing out. Thus, ik = k + (55) where q is the average current density entering the coated surface of the element in the radial direction. Ac and As represent the areas of the steel crosssection and side walls, respectively. For the cathodic protected system, the protecting current flow away from the anode to the soil, and then flow to the pipeline. On the pipeline, the currents from the two sides of the connection point cl flow head Figure 53: Axial direction current flow along the pipeline. to head and are concurrent as current I flows back to the anode through the wire which is connecting the pipe and the anode. Therefore, for each segment of the pipe before connection point cl, the potential difference is, Vo V = (io + i)R = [2io (A'. (q0 + q ) ]R V1 V2 = (i 2)R = [2i(A (qo + 2ql + 2) ,A) 2 V, V, = ( = (ic'+ic) = [2io (A (qo + 2q + + 2q1 + ql)]R For each segment of the pipe after the connection point cl SA (q'+2q" P^1+...2qcl q1+ql) Vc+l+V, = (ic+l +ic)R=[2in+(Ac. +2q" +( +2qc + qc] Vn+Vn 1 = (in+i 1)R=[2in+ (A (q + )]R There are n equations if the pipeline is discritized as n segments, while there are n + 1 the unknown variables Vo, V1, Vn. One more equation is needed in order to get a unique solution. The potential V = 0 can be set at any position along the pipe. In this study, zero potential was set at the pipe and anode connection point cl, that is, Vc, = 0. The current flow in the anode through the wire will flow to the different direction along the anode. Likewise, for each segment of the anode Pipeline I Anode Figure 54: Axial direction before the connection point c2, Vo V = [2io V V [2 () A current flow along the anode. S(q+0 qi) 2 (q0 + 2q + q2) 2 S[2i (Aj+2 (o+2q1+ +2q2+ qC) ( C) For each segment of the anode after connection point C2 S(Aa (qm + 2qm1 + .. .+2qC2+1 + qC2 VC+2i+V = [2im+ A 2 Vm + Vm1 = [2im+ A (m + q2 ] Corresponding to the m segments of the anode, there will be m + 1 unknown vari ables Vo, V1, Vm. Another equation is needed. If the wire connecting the pipe and the anode has resistance Rwire, and the current flow through the wire is I, then V, = Vc2 IRire = 0 (56) In order to satisfy the selfequilibrium limitation, the current I needs to obey the relations as below. For the pipeline, as shown in Figure 53, ic = io _A\ (q + 2q1 +.+ 2qcl + qc) (57) FO ~AP 2 ic A A, (qCl +2qcl+ +.. +2q +q) (5 Fn AP 2 (58) where ifO is the current density coming from left end of pipe to the point cl, and ic is the current density coming from right end of pipe to the point cl. By combining these two equations, (io ic )AP = I (59) We have (io i")A ( 2 2q ) AP = I (510) For the anode, as shown in Figure 54, C2 = i A(As' (qO + 2q + .+ 2qC21 + qC2) iO = 2 (511) ro\ A 2 S= m A (C2 + 2qC21 + ... + 2qm + qm) (512) Tm + kA" 2 where iC0 is the current density enter at C2 and flow to the upper end of anode, and iT 2 is the current density enter at c2 and flow to the lower end of anode. Again, by combining these two equations, (iT 0 iT2)A = I (513) We have ( i im)A (q 2' (514) Summarizing the n + 1 equations for pipeline in matrix form, and defining K,, F,, Vp, qp and iend as 1 0 1 1 0 1 0 1 0 1 1 ... ... ... ... 0 1 .. 1 0 ...... 010** S0 . 0 0 0 0 0 0 2R A, 2( IA?2 ) 2 0 0 0 1 0 1 2 q0 qCl 1  S= Fpq KV =Ti ends + Fpqp p P p Vo V, VI, Vc +1 Vn1 Vn We get (515) Likewise, define Ka, Fa, V,, q, and i"d as Ka Ta = 2R 1 1 0 S. 0 0 1 1 0 0 1 1 0 0 0 1 1 0 *** 0 0 1 1 . .... ... .. 1  0 1 0 AcR AcR 2R wire 2R Rwire As Aaj 2 c Va = The m + 1 1 1 1 0 0 A Rwire Vo V1 C2 e mindss = 0 ,qa qC2 ia VC2+ qC2 Vm t qM Vn equations of the anode in matrix form are, KaVa = Tai d + Faaq Introducing the equations (515), (516) into equation (461) in the previous chap 0 0 I U 2 0 1 0 1 2 0 1 1 A Rwire (516) ter, we obtain (n V7), K, 0 0 0 0 Vp Fp Tp 0 0 0 Ka 0 0 0 Va 0 0 F, T, in 0 0 Hy, Hpa 1 < p = Gpp 0 Gpa 0 (n V )a 0 0 Hap Haa 1 Oa Gap 0 Gaa 0 1a 0 0 0 0 0 Ap Aa i" (517) The forward model, equation (517), can be used to calculate both the potential on the soil surface, and pipe steel potentials. Furthermore, pipe steel potentials allows calculation of the current flowing in the pipe, which can be compared with the field measurements. 5.1.3 Simulation Results In order to show the simulation results with variant pipe steel potentials, a test case of CP system was studied. A 500 m long pipe is assumed to be buried 1.45 m (4.75 feet) below the soil surface. Its diameter is 0.457 m (1.5 feet). A anode is placed 50 m away from the center of the pipe. The diameter of the anode is 0.152 m and its length is 1.22 m (4 feet). The soil resistivity is 100 1Qm. The boundary conditions are shown in Figure 41. A 0.5 mmthick coating is assumed to cover the side area of the pipe, with the exception at the two end, which are assumed to be insulated. A wire connected the pipeline to the anode at 50 m position and 0.5 m position, respectively. The potential of the pipe steel varies along the pipe. A coating defect is set on the middle of the length of pipeline. The simulation results of anode and pipeline are shown in Figure 55 and Fig ure 56 respectively. I I I I I I I I ' 24 22 20 1 8 1 6 1 4 Position /m 24 20 Position /m 1.0x107 75 2.0x107 3.0x10' 2.4 2.2 2.0 1.8 1.6 1.4 1.2 Position /m 012 24 22 20 1 8 1 6 Position /m 14 1 2 Figure 55: Potential and current density along the anode: a) potential Y along the anode (given value); b) potential V along the anode; c) radical direction of current density along the anode; d) axial direction of current density along the anode. 0020 E 0016 .i  0012 0 0008 I I, E 60 C r 50 S40 "3 30 S20 co o 10 0 0 100 200 300 Position / m (a) 0 E 20 40 0) " 40 O 60 500 0 100 200 300 Position / m I I I I I 0 100 200 300 400 50 Position / m I I I I I I I I 0.020  0.016 c 0.012 () 0 a_ 0.008 () (D g 0.004 0 0 100 200 300 Position / m c, U.4 I I E location Q E 0   location of bond 'c 0.4 (D 0 C S1.2 1.6 I 0 100 200 300 Position / m 400 500 400 500 Figure 56: Simulated axial distributions along the pipeline: a) input value for coating resistivity; b) radial component of current density; c) calculated value for potential 4; d) calculated value for steel potential V; and e) axial component of current density in the pipe. 5.1.4 Analysis of the Simulation Results Figure 5.5(a) shows that 4a = 1.1V, because the anode is a sacrificial anode and its potential is a given boundary condition. Figure 5.5(b) indicates the poten tials along the anode steel are variant. At the position where the potential equals zero, a wire connects the anode to the pipeline. In order to protect the pipeline, the current of the anode flows away from the anode, and flows back to the anode through the wire. Therefore, the farther the position is from the wire connected position, the lower the anode steel potential than that of the potential at the wire connection point. Figure 5.5(c) exhibits the radial direction of current density. Fig ure 5.5(d) shows that the axial direction of current density is equal to zero at the two ends of the anode. This is a given boundary condition. The current densities at the two sides of the wire connection point have different signs, which means that the current flows away from the wire connecting position to two opposite directions. The coating resistivity of the pipeline is shown in Figure 5.6(a) where the nom inal resistivity of defectfree coating is seen to be 5.0 x 107'2m. The resistivity was assumed to decrease abruptly at a position of 250 m, corresponding to the posi tion of a significant coating defect. The cathodic radial current density increased dramatically at the position of the coating defect, as shown in Figure 5.6(b). Figure 5.6(c) shows that the potential is less negative at the coating defect position than other places. It means that pipe is less protected at the coating defect. The variation of potential within the pipe steel is presented in Figure 5.6(d). Since the length of the 500 m pipe was not large, the potential drop over the pipe was very small, on the order of 10 3 10 2 /V. The calculation of steel potential was primarily useful for allowing calculation of the current flow in the pipe. The steel potential was set to a value of zero at the 50 m position where the anode and 