PDF
Stress Corrosion
Introduction
Steel pipelines are often subjected to complex stress/strain conditions in oil and gas industry. In addition to stress from internal pressure, pipelines are subjected to significant longitudinal strain due to surrounding soil movement. The effect of elastic and plastic deformation on pipeline corrosion is demonstrated in this model example.
The elasto-plastic stress simulations are performed here using a small strain plasticity model and von-Mises yielding criterion. Iron dissolution (anodic) and hydrogen evolution (cathodic) are considered as electrochemical reactions, using kinetic expressions that account for the effect of elasto-plastic deformations.
The example is based on a paper by L. Y. Xu and Y. F. Cheng (Ref. 1).
Note: This model requires the Structural Mechanics, or the MEMS Module, with the addition of either the Nonlinear Structural Materials or the Geomechanics Module.
Model Definition
The model geometry is shown in Figure 1.
Figure 1: The model geometry consists of a pipeline with corrosion defect and surrounding soil domain.
The model geometry consists of high strength alloy steel pipeline and surrounding soil domain. The pipeline length is 2 m and wall thickness is 19.1 mm. The corrosion defect on the exterior side of the pipeline is elliptical in shape with a length of 200 mm and a depth of 11.46 mm. The electrolyte conductivity of soil domain in 0.096 S/m.
Elastoplastic Stress
An elastoplastic stress simulation is performed over pipeline domain using the small strain plasticity model. The user defined isotropic hardening model is used where the hardening function, σyhard, is defined as:
where σexp is the experimental stress-strain curve, εp is the plastic deformation, σe is the von Mises stress, E is the Young’s modulus (207·109 Pa), and σys is the yield strength of high strength alloy steel(806·106 Pa).
The experimental stress-strain curve used in the model is prescribed in terms of a piecewise cubic interpolation function and is taken from Ref. 2.
Electrochemical Reactions
The iron dissolution (anodic) and hydrogen evolution (cathodic) reactions are the two electrochemical reactions that occur at the corrosion defect surface of pipelines. The rest of pipeline surfaces are assumed to be electrochemically inactive.
An anodic Tafel expression is used to model the iron dissolution reaction, with a local anodic current density defined as
where i0,a is the exchange current density (2.353·10−3 A/m2), Aa is the Tafel slope (0.118 V) and the overpotential ηa for the anodic reaction is calculated from
The equilibrium potential for the anodic reaction is calculated from
where Eeq0,a is the standard equilibrium potential for the anodic reaction (−0.859 V), ΔPm is the excess pressure to elastic deformation (2.687·108 Pa), Vm is the molar volume of steel (7.13·10−6 m3/mol), z is the charge number for steel (2), F is the Faraday’s constant, T is the absolute temperature (298.15 K), R is the ideal gas constant, νis an orientation dependent factor (0.45), α is a coefficient (1.67·1015 m-2) and N0 is the initial dislocation density (1·1012 m-2).
A cathodic Tafel expression is used to model the iron dissolution reaction, this sets the local cathodic current density to
where i0,c is the exchange current density, Ac is the Tafel slope (−0.207 V) and the overpotential ηc (SI unit: V) for the cathodic reaction is calculated from
where Eeq0,c is the standard equilibrium potential for the cathodic reaction (−0.644 V)
The exchange current density for the cathodic reaction is calculated from
where i0,c,ref is the reference exchange current density for the cathodic reaction in the absence of external stress/strain (1.457·10−2 A/m2).
Results and Discussion
Figure 2 shows the electrolyte potential distribution (V) over the soil domain and the von Mises stress distribution (MPa) over the pipe domain, as indicated by the color bars for a prescribed displacement of 4 mm in the x direction. It can be seen that the local stresses are significantly higher near the corrosion defect than that at the rest of the pipeline. A nonuniform electrolyte potential distribution near the corrosion defect is also evident in Figure 2, as indicated by a semi-circular area.
Figure 2: The electrolyte potential distribution over the soil domain and the von Mises stress distribution over the pipeline domain for a prescribed displacement of 4 mm.
Figure 3 shows the von Mises stress distribution along the length of the corrosion defect for prescribed displacements of 1.375 mm, 2.75 mm, 3.75 mm, and 4 mm, respectively. The von Mises stress increases with an increase in the tensile strain and it is found to be maximum at the center of corrosion defect. For the tensile strain of 3.75 mm and 4 mm, it is observed that the local stress, particularly at the center of corrosion defect, exceeds the yield strength of high strength alloy steel (806·106 Pa). This results in the plastic deformation at the center of corrosion defect while deformation in the remaining area of corrosion defect remains in the elastic range. For the lower tensile strains of 1.375 mm, and 2.75 mm, the entire corrosion defect is observed to be in the elastic deformation range.
Figure 3: The von Mises stress distribution along the length of the corrosion defect for prescribed displacements of 1.375 mm, 2.75 mm, 3.75 mm, and 4 mm.
Figure 4 shows the corrosion potential distribution along the length of the corrosion defect for prescribed displacements of 1.375 mm, 2.75 mm, 3.75 mm, and 4 mm. For lower tensile strains of 1.375 mm and 2.75 mm, the variation in the corrosion potential is found to be uniform along the length of the corrosion defect. However, for higher tensile strains of 3.75 mm and 4 mm, the variation in the corrosion potential is nonuniform with the more negative corrosion potential at the center of the corrosion defect than that at both the sides of the corrosion defect.
Figure 4: The corrosion potential distribution along the length of the corrosion defect for prescribed displacements of 1.375 mm, 2.75 mm, 3.75 mm, and 4 mm.
Figure 5 shows the anodic current density distribution along the length of corrosion defect for prescribed displacements of 1.375 mm, 2.75 mm, 3.75 mm, and 4 mm. For lower tensile strains of 1.375 mm and 2.75 mm, the variation in the anodic current density is found to be uniform along the length of the corrosion defect, similar to the corrosion potential behavior. However, for higher tensile strains of 3.75 mm and 4 mm, the variation in the anodic current density is significantly nonuniform, particularly at the center of the corrosion defect. It can be seen that the anodic current density increases significantly at the center of the corrosion defect whereas it decreases slightly at both the sides of the corrosion defect for higher tensile strains. The increase in the anodic current density for tensile strains of 3.75 mm and 4 mm is attributed to the plastic deformation observed at the center of the corrosion defect (see Figure 3).
Figure 5: The anodic current density distribution along the length of the corrosion defect for prescribed displacements of 1.375 mm, 2.75 mm, 3.75 mm, and 4 mm.
Figure 6 shows the cathodic current density distribution along the length of the corrosion defect for prescribed displacements of 1.375 mm, 2.75 mm, 3.75 mm, and 4 mm. It can be seen that the cathodic current density increases negatively with an increase in the tensile strain and it is found to be the most negative at the center of the corrosion defect. The nonuniformity in the cathodic current density is also found to increase with an increase in the tensile strain. Thus, the cathodic current density distribution is found to be the most nonuniform for a tensile strain of 4 mm.
Figure 6: The cathodic current density distribution along the length of the corrosion defect for prescribed displacements of 1.375 mm, 2.75 mm, 3.75 mm, and 4 mm.
Notes About the COMSOL Implementation
The model is implemented using the Solid Mechanics interface and the Secondary Current Distribution interface. Since the Solid Mechanics interface does not depend upon the Secondary Current Distribution interface, we use a sequential solver set up with a Parametric Sweep to study the impact of elastoplastic deformations on electrochemical reactions.
References
1. L. Y. Xu and Y. F. Cheng, “Development of a finite element model for simulation and prediction of mechanoelectrochemical effect of pipeline corrosion”, Corrosion Science, vol. 73, pp. 150–160, 2013.
2. L. Xu, “Assessment of corrosion defects on high-strength steel pipelines”, PhD thesis, Department of mechanical and manufacturing engineering, University of Calgary, Alberta, August 2013.
Application Library path: Corrosion_Module/Galvanic_Corrosion/stress_corrosion
Modeling Instructions
From the File menu, choose New.
New
In the New window, click  Model Wizard.
Model Wizard
1
In the Model Wizard window, click  2D.
2
In the Select Physics tree, select Structural Mechanics>Solid Mechanics (solid).
3
Click Add.
4
In the Select Physics tree, select Electrochemistry>Primary and Secondary Current Distribution>Secondary Current Distribution (cd).
5
Click Add.
6
Click  Study.
7
In the Select Study tree, select General Studies>Stationary.
8
Click  Done.
Geometry 1
Draw the geometry for pipe with corrosion defect and surrounding soil domain.
Rectangle 1 (r1)
1
In the Geometry toolbar, click  Rectangle.
2
In the Settings window for Rectangle, locate the Size and Shape section.
3
In the Width text field, type 2.
4
In the Height text field, type 19.1[mm].
Ellipse 1 (e1)
1
In the Geometry toolbar, click  Ellipse.
2
In the Settings window for Ellipse, locate the Size and Shape section.
3
In the a-semiaxis text field, type 100[mm].
4
In the b-semiaxis text field, type 11.46[mm].
5
Locate the Position section. In the x text field, type 1.
6
In the y text field, type 19.1[mm].
Difference 1 (dif1)
1
In the Geometry toolbar, click  Booleans and Partitions and choose Difference.
2
Select the object r1 only.
3
In the Settings window for Difference, locate the Difference section.
4
Find the Objects to subtract subsection. Select the  Activate Selection toggle button.
5
Select the object e1 only.
Rectangle 2 (r2)
1
In the Geometry toolbar, click  Rectangle.
2
In the Settings window for Rectangle, locate the Size and Shape section.
3
In the Width text field, type 2.
4
Click  Build All Objects.
5
Click the  Zoom Extents button in the Graphics toolbar.
Your geometry should look like Figure 1.
Global Definitions
Parameters 1
Load the model parameters.
1
In the Model Builder window, under Global Definitions click Parameters 1.
2
In the Settings window for Parameters, locate the Parameters section.
3
Click  Load from File.
4
Browse to the model’s Application Libraries folder and double-click the file stress_corrosion_parameters.txt.
Definitions
Load the model variables.
Variables 1
1
In the Model Builder window, expand the Component 1 (comp1)>Definitions node.
2
Right-click Definitions and choose Variables.
3
In the Settings window for Variables, locate the Variables section.
4
Click  Load from File.
5
Browse to the model’s Application Libraries folder and double-click the file stress_corrosion_variables.txt.
Interpolation 1 (int1)
Load the stress strain interpolation data from a text file.
1
In the Home toolbar, click  Functions and choose Local>Interpolation.
2
In the Settings window for Interpolation, locate the Definition section.
3
In the Function name text field, type stress_strain_curve.
4
Click  Load from File.
5
Browse to the model’s Application Libraries folder and double-click the file stress_corrosion_stress_strain_curve_interpolation.txt.
6
Locate the Interpolation and Extrapolation section. From the Interpolation list, choose Piecewise cubic.
7
Locate the Units section. In the Arguments text field, type 1.
8
In the Function text field, type MPa.
Solid Mechanics (solid)
Start setting up the physics. First, set the elastoplastic deformation at the Linear Elastic Material node by adding Plasticity node.
1
In the Model Builder window, under Component 1 (comp1) click Solid Mechanics (solid).
2
In the Settings window for Solid Mechanics, locate the Domain Selection section.
3
In the list, select 2.
4
Click  Remove from Selection.
5
Select Domain 1 only.
Linear Elastic Material 1
In the Model Builder window, under Component 1 (comp1)>Solid Mechanics (solid) click Linear Elastic Material 1.
Plasticity 1
1
In the Physics toolbar, click  Attributes and choose Plasticity.
2
In the Settings window for Plasticity, locate the Plasticity Model section.
3
Find the Isotropic hardening model subsection. From the list, choose Hardening function.
Materials
Now, add high-strength alloy steel material for pipe and set the values for initial yield stress, hardening function, Young’s modulus and Poisson’s ratio.
Add Material
1
In the Home toolbar, click  Add Material to open the Add Material window.
2
Go to the Add Material window.
3
In the tree, select Built-in>High-strength alloy steel.
4
Click Add to Component in the window toolbar.
5
In the Home toolbar, click  Add Material to close the Add Material window.
Materials
High-strength alloy steel (mat1)
1
In the Settings window for Material, locate the Geometric Entity Selection section.
2
In the list, select 2.
3
Click  Remove from Selection.
4
Select Domain 1 only.
In the table under Material Contents, set the value of Initial yield stress to 806e6 [Pa] and Hardening function to hardening. Also, change the value of Young’s modulus to 207e9 [Pa] and Poisson’s ratio to 0.33.
Solid Mechanics (solid)
Now, set the initial value for displacement field and then proceed to setting up boundary conditions for Solid Mechanics physics interface.
Initial Values 1
1
In the Model Builder window, under Component 1 (comp1)>Solid Mechanics (solid) click Initial Values 1.
2
In the Settings window for Initial Values, locate the Initial Values section.
3
Specify the u vector as
0.0001*X
X
0
Y
Fixed Constraint 1
1
In the Physics toolbar, click  Boundaries and choose Fixed Constraint.
2
Select Boundary 1 only.
Prescribed Displacement 1
1
In the Physics toolbar, click  Boundaries and choose Prescribed Displacement.
2
Select Boundary 7 only.
3
In the Settings window for Prescribed Displacement, locate the Prescribed Displacement section.
4
Select the Prescribed in x direction check box.
5
In the u 0 x text field, type disp.
Prescribed Displacement 2
1
In the Physics toolbar, click  Boundaries and choose Prescribed Displacement.
2
Select Boundary 2 only.
3
In the Settings window for Prescribed Displacement, locate the Prescribed Displacement section.
4
Select the Prescribed in y direction check box.
Secondary Current Distribution (cd)
Now, set up the physics for electrochemical reactions. First, set electrolyte conductivity, initial value for electrolyte potential and then set both anodic and cathodic reactions at the corrosion defect surface of pipe.
1
In the Model Builder window, under Component 1 (comp1) click Secondary Current Distribution (cd).
2
In the Settings window for Secondary Current Distribution, locate the Domain Selection section.
3
In the list, select 1.
4
Click  Remove from Selection.
5
Select Domain 2 only.
Electrolyte 1
1
In the Model Builder window, under Component 1 (comp1)>Secondary Current Distribution (cd) click Electrolyte 1.
2
In the Settings window for Electrolyte, locate the Electrolyte section.
3
From the σl list, choose User defined. In the associated text field, type sigmal.
Initial Values 1
1
In the Model Builder window, click Initial Values 1.
2
In the Settings window for Initial Values, locate the Initial Values section.
3
In the phil text field, type -(Eeq0a+Eeq0c)/2.
Electrode Surface 1
1
In the Physics toolbar, click  Boundaries and choose Electrode Surface.
2
Select Boundaries 9 and 10 only.
Electrode Reaction 1
1
In the Model Builder window, expand the Electrode Surface 1 node, then click Electrode Reaction 1.
2
In the Settings window for Electrode Reaction, locate the Equilibrium Potential section.
3
In the Eeq text field, type Eeqa.
4
Locate the Electrode Kinetics section. From the Kinetics expression type list, choose Anodic Tafel equation.
5
In the i0 text field, type i0a.
6
In the Aa text field, type ba.
Electrode Surface 1
In the Model Builder window, click Electrode Surface 1.
Electrode Reaction 2
1
In the Physics toolbar, click  Attributes and choose Electrode Reaction.
2
In the Settings window for Electrode Reaction, locate the Equilibrium Potential section.
3
In the Eeq text field, type Eeq0c.
4
Locate the Electrode Kinetics section. From the Kinetics expression type list, choose Cathodic Tafel equation.
5
In the i0 text field, type ic.
6
In the Ac text field, type bc.
Mesh 1
Set the finer mesh near corrosion defect surface of pipe.
Free Triangular 1
1
In the Mesh toolbar, click  Free Triangular.
2
In the Settings window for Free Triangular, locate the Domain Selection section.
3
From the Geometric entity level list, choose Domain.
4
From the Selection list, choose All domains.
Size 1
1
Right-click Free Triangular 1 and choose Size.
2
In the Settings window for Size, locate the Geometric Entity Selection section.
3
From the Geometric entity level list, choose Boundary.
4
Select Boundaries 2, 9, and 10 only.
5
Locate the Element Size section. Click the Custom button.
6
Locate the Element Size Parameters section. Select the Maximum element size check box.
7
In the associated text field, type 0.002.
8
Click  Build All.
Study 1
Now, set the solver settings. Since solid mechanics physics is not dependent upon electrochemical reactions, we use a sequential solver setup with a parametric sweep.
Parametric Sweep
1
In the Study toolbar, click  Parametric Sweep.
2
In the Settings window for Parametric Sweep, locate the Study Settings section.
3
Click  Add.
4
In the table, enter the following settings:
Parameter name
Parameter value list
Parameter unit
disp (Displacement)
0.001375 0.00275 0.00375 0.004
m
Step 1: Stationary
1
In the Model Builder window, click Step 1: Stationary.
2
In the Settings window for Stationary, locate the Physics and Variables Selection section.
3
In the table, clear the Solve for check box for Secondary Current Distribution (cd).
Stationary 2
1
In the Study toolbar, click  Study Steps and choose Stationary>Stationary.
2
In the Settings window for Stationary, locate the Physics and Variables Selection section.
3
In the table, clear the Solve for check box for Solid Mechanics (solid).
Solution 1 (sol1)
Lower the relative tolerance for Secondary Current Distribution study step.
1
In the Study toolbar, click  Show Default Solver.
2
In the Model Builder window, expand the Solution 1 (sol1) node, then click Stationary Solver 2.
3
In the Settings window for Stationary Solver, locate the General section.
4
In the Relative tolerance text field, type 0.00001.
Clear the Generate default plots check box. The model is now ready to be solved.
5
In the Model Builder window, click Study 1.
6
In the Settings window for Study, locate the Study Settings section.
7
Clear the Generate default plots check box.
8
In the Study toolbar, click  Compute.
Results
Corrosion potential and von Mises stress
1
In the Home toolbar, click  Add Plot Group and choose 2D Plot Group.
2
In the Settings window for 2D Plot Group, type Corrosion potential and von Mises stress in the Label text field.
3
Locate the Data section. From the Dataset list, choose Study 1/Parametric Solutions 1 (sol3).
Surface 1
1
In the Corrosion potential and von Mises stress toolbar, click  Surface.
2
In the Settings window for Surface, click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 1 (comp1)>Solid Mechanics>Stress>solid.mises - von Mises stress - N/m².
3
Locate the Expression section. From the Unit list, choose MPa.
Corrosion potential and von Mises stress
In the Model Builder window, click Corrosion potential and von Mises stress.
Surface 2
1
In the Corrosion potential and von Mises stress toolbar, click  Surface.
2
In the Settings window for Surface, locate the Expression section.
3
In the Expression text field, type -phil.
4
Click to expand the Range section. Select the Manual color range check box.
5
In the Minimum text field, type -0.733.
6
In the Maximum text field, type -0.724.
7
Locate the Coloring and Style section. Select the Reverse color table check box.
8
In the Corrosion potential and von Mises stress toolbar, click  Plot.
9
Click the  Zoom Extents button in the Graphics toolbar.
The plot should like Figure 2. One can zoom in a region close to corrosion defect using Zoom Box icon in the Graphics window.
von Mises Stress
Plot von Mises stress along the corrosion defect for different values of prescribed displacements.
1
In the Home toolbar, click  Add Plot Group and choose 1D Plot Group.
2
In the Settings window for 1D Plot Group, type von Mises Stress in the Label text field.
3
Locate the Data section. From the Dataset list, choose Study 1/Parametric Solutions 1 (sol3).
4
Locate the Plot Settings section. Select the x-axis label check box.
5
In the associated text field, type Defect length (mm).
Line Graph 1
1
In the von Mises Stress toolbar, click  Line Graph.
2
Select Boundaries 9 and 10 only.
3
In the Settings window for Line Graph, click Replace Expression in the upper-right corner of the y-Axis Data section. From the menu, choose Component 1 (comp1)>Solid Mechanics>Stress>solid.mises - von Mises stress - N/m².
4
Locate the y-Axis Data section. From the Unit list, choose MPa.
5
Locate the x-Axis Data section. From the Parameter list, choose Expression.
6
In the Expression text field, type x.
7
From the Unit list, choose mm.
8
Click to expand the Coloring and Style section. Find the Line markers subsection. From the Marker list, choose Cycle.
9
Click to expand the Legends section. Select the Show legends check box.
10
From the Legends list, choose Manual.
11
In the table, enter the following settings:
Legends
1.375 mm
2.75 mm
3.75 mm
4 mm
12
In the von Mises Stress toolbar, click  Plot.
The plot should look like Figure 3.
von Mises Stress
Now, plot corrosion potential along defect length for different prerscribed displacements.
Corrosion potential
1
In the Model Builder window, right-click von Mises Stress and choose Duplicate.
2
In the Settings window for 1D Plot Group, type Corrosion potential in the Label text field.
3
Locate the Plot Settings section. Select the y-axis label check box.
4
In the associated text field, type Corrosion potential (V).
Line Graph 1
1
In the Model Builder window, expand the Corrosion potential node, then click Line Graph 1.
2
In the Settings window for Line Graph, click Replace Expression in the upper-right corner of the y-Axis Data section. From the menu, choose Component 1 (comp1)>Secondary Current Distribution>cd.Evsref - Electrode potential vs. adjacent reference - V.
3
In the Corrosion potential toolbar, click  Plot.
The plot should look like Figure 4.
Corrosion potential
Plot the anodic current density along the corrosion defect.
Anodic current density
1
In the Model Builder window, right-click Corrosion potential and choose Duplicate.
2
In the Settings window for 1D Plot Group, type Anodic current density in the Label text field.
3
Locate the Plot Settings section. In the y-axis label text field, type Anodic current density (A/m<sup>2</sup>).
Line Graph 1
1
In the Model Builder window, expand the Anodic current density node, then click Line Graph 1.
2
In the Settings window for Line Graph, click Replace Expression in the upper-right corner of the y-Axis Data section. From the menu, choose Component 1 (comp1)>Secondary Current Distribution>Electrode kinetics>cd.iloc_er1 - Local current density - A/m².
3
In the Anodic current density toolbar, click  Plot.
The plot should look like Figure 5.
Anodic current density
Plot the cathodic current density along the corrosion defect.
Cathodic current density
1
In the Model Builder window, right-click Anodic current density and choose Duplicate.
2
In the Settings window for 1D Plot Group, type Cathodic current density in the Label text field.
3
Locate the Plot Settings section. In the y-axis label text field, type Cathodic current density (A/m<sup>2</sup>).
Line Graph 1
1
In the Model Builder window, expand the Cathodic current density node, then click Line Graph 1.
2
In the Settings window for Line Graph, click Replace Expression in the upper-right corner of the y-Axis Data section. From the menu, choose Component 1 (comp1)>Secondary Current Distribution>Electrode kinetics>cd.iloc_er2 - Local current density - A/m².
3
In the Cathodic current density toolbar, click  Plot.
The plot should look like Figure 6.