Corrosion researcher: MQL recipes

Recipes

Corrosion rate from coupon mass loss (ASTM G1/G31)
Driving force by Nernst, penetration by Faraday
Tafel slopes and Stern-Geary corrosion current
Service life of a coating by Arrhenius extrapolation
Chloride ingress and corrosion initiation in concrete
Pourbaix diagram of iron in water
EIS fit with a CPE and inhibition efficiency
Inhibitor efficiency statistics and Langmuir adsorption
Deepest pit on a pipe by Gumbel extreme values
Galvanic couple steel-copper: Evans mixed potential

executed on MentDB

Corrosion rate from coupon mass loss (ASTM G1/G31)

Three carbon steel coupons (28.5 cm2, density 7.86 g/cm3) lost 0.1523, 0.1487 and 0.1561 g in 720 h of immersion. The script takes the mean mass loss with its 95 % Student interval, converts CR = W/(A.t.rho) directly to mm/y with dimensional analysis (no 87.6 or 8.76e4 factor to remember), and also gives mpy.

Pipeline

01 · statsMean mass loss and 95 % intervalstat conf_interval_mean→02 · extractMean and marginjson select→03 · rateCR = W / (A.t.rho) in mm/yrunit calc→04 · reportRate with its uncertaintymath round_uncertainty

Script

admin

#1. Triplicate carbon steel coupons, 30 days in aerated brine, ASTM G1 / G31;
-> "[w]" "[0.1523, 0.1487, 0.1561]";
json load "rx_ch_ml" (stat conf_interval_mean [w] 0.95);
-> "[wm]" (json select "rx_ch_ml" /mean); -> "[dw]" (json select "rx_ch_ml" /margin);
#2. CR = W / A.t.rho, with A = 28.5 cm2, t = 720 h, rho = 7.86 g/cm3;
-> "[cr]" (unit calc (concat [wm] " g / (28.5 cm2 * 720 h * 7.86 g/cm3)") "mm/yr");
-> "[dcr]" (math calc "[cr]*[dw]/[wm]");
#3. Same rate in mils per year, then report;
-> "[mpy]" (unit calc (concat [cr] " mm/yr") "in/yr");
concat "CR = " (math round_uncertainty [cr] [dcr]) " mm/y at 95 % = " (math round (math calc "[mpy]*1000") 2) " mpy"

Result

mentdb

CR = 0.083 ± 0.005 mm/y at 95 % = 3.26 mpy

Go further

Compare two inhibitors with stat ttest_2 on the coupon rates.
Convert to g/(m2*d) with unit calc for the mass-loss rate.

executed on MentDB

Driving force by Nernst, penetration by Faraday

Will steel corrode in a deaerated acid at pH 4, and how fast for a measured anodic current density of 25 uA/cm2? The script computes the Nernst potentials of Fe2+/Fe at 1e-6 mol/L (Pourbaix convention) and of H+/H2, then applies Faraday's law (Fe -> Fe2+ + 2e-) to get the penetration rate in mm/y and the yearly mass loss in kg/m2.

Pipeline

01 · nernstNernst potentials E0 + RT/nF.ln(a)math calc→02 · faradayCR = i.M / (n.F.rho) in mm/yrunit calc→03 · massMass loss per m2 and per yearunit calc→04 · reportSummaryconcat

Script

admin

#1. Nernst at 25 degC: Fe2+/Fe at 1e-6 mol/L, Pourbaix convention, and H+/H2 at pH 4;
-> "[t]" 298.15;
-> "[e_fe]" (math calc "-0.447 + _R*[t]/(2*_F)*ln(1e-6)");
-> "[e_h]" (math calc "-ln(10)*_R*[t]/_F*4");
#2. Measured anodic current density 25 uA/cm2: Faraday law, Fe to Fe2+;
-> "[cr]" (unit calc "25 uA/cm2 * 55.845 g/mol / (2 * _F * 7.86 g/cm3)" "mm/yr");
-> "[ml]" (unit calc "25 uA/cm2 * 55.845 g/mol / (2 * _F) * 1 yr" "kg/m2");
concat "E(Fe2+/Fe) = " (math round [e_fe] 3) " V, E(H+/H2) = " (math round [e_h] 3) " V, driving force " (math round (math calc "[e_h]-[e_fe]") 3) " V; CR = " (math round [cr] 3) " mm/y, mass loss " (math round [ml] 2) " kg/m2 per year"

Result

mentdb

E(Fe2+/Fe) = -0.624 V, E(H+/H2) = -0.237 V, driving force 0.388 V; CR = 0.29 mm/y, mass loss 2.28 kg/m2 per year

Go further

Draw the Fe/H2O Pourbaix lines with num table over the pH.
Use const get "F" to show the CODATA Faraday constant.

executed on MentDB

Tafel slopes and Stern-Geary corrosion current

A potentiodynamic scan of +-150 mV around Ecorr and a separate LPR scan of +-10 mV were recorded on steel. The script fits the Butler-Volmer equation to get icorr and the Tafel slopes ba and bc, measures the polarization resistance Rp from the LPR slope, applies Stern-Geary icorr = B/Rp with B = ba.bc/(2.303(ba + bc)), and converts icorr to a corrosion rate (ASTM G102, EW = 27.92).

Pipeline

01 · tafelButler-Volmer fit: icorr, ba, bcstat fit_model→02 · lprRp from the linear LPR slopestat fit→03 · stern-gearyB / Rpmath calc→04 · rateCR = icorr.EW / (F.rho) in mm/yrunit calc

Script

admin

#1. Tafel scan +-150 mV around Ecorr: E = overpotential in mV, i in uA/cm2;
-> "[eta]" "[-150, -125, -100, -75, -50, -25, -10, 10, 25, 50, 75, 100, 125, 150]";
-> "[i]" "[-213.85, -135.3, -79.92, -50.76, -29.13, -14.36, -6.35, 7.22, 21.65, 66.92, 172.04, 412.86, 1016.74, 2391.27]";
#2. Butler-Volmer fit: icorr and Tafel slopes ba, bc in mV/decade;
json load "rx_ch_bv" (stat fit_model "icorr*(10^(E/ba) - 10^(-E/bc))" "E" [eta] [i] "{\"icorr\": 10, \"ba\": 60, \"bc\": 120}");
-> "[ic]" (json select "rx_ch_bv" /parameters/icorr); -> "[ba]" (json select "rx_ch_bv" /parameters/ba); -> "[bc]" (json select "rx_ch_bv" /parameters/bc);
#3. Separate LPR scan +-10 mV: Rp = dE/di in kohm.cm2, Stern-Geary icorr = B / Rp;
json load "rx_ch_lpr" (stat fit "[-10, -5, 0, 5, 10]" "[-6.16, -3.08, -0.03, 3.38, 7.06]" "linear");
-> "[rp]" (math calc (concat "1/" (json select "rx_ch_lpr" /b)));
-> "[isg]" (math calc "[ba]*[bc]/(ln(10)*([ba]+[bc]))/[rp]");
#4. ASTM G102 corrosion rate from the Tafel icorr, steel EW = 27.92 g, rho = 7.86 g/cm3;
-> "[cr]" (unit calc (concat [ic] " uA/cm2 * 27.92 g/mol / (_F * 7.86 g/cm3)") "mm/yr");
concat "ba = " (math round [ba] 1) " mV/dec, bc = " (math round [bc] 1) " mV/dec, icorr Tafel = " (math round [ic] 2) " uA/cm2, Rp = " (math round [rp] 3) " kohm.cm2, icorr Stern-Geary = " (math round [isg] 2) " uA/cm2, CR = " (math round [cr] 3) " mm/y"

Result

mentdb

ba = 66.5 mV/dec, bc = 124.4 mV/dec, icorr Tafel = 13.27 uA/cm2, Rp = 1.52 kohm.cm2, icorr Stern-Geary = 12.38 uA/cm2, CR = 0.154 mm/y

Go further

Weight the Tafel points with the current noise (sigma argument of stat fit_model).
Check the cathodic branch for diffusion control by fitting a limiting current term.

executed on MentDB

Service life of a coating by Arrhenius extrapolation

Accelerated immersion tests of an epoxy coating gave times to blistering of 8446, 3463, 1481 and 678 h at 60, 70, 80 and 90 degC. The script fits the Arrhenius lifetime law L = L80.exp(Ea/R.(1/T - 1/353.15)) and extrapolates the service life at 25 degC in years.

Pipeline

01 · fitArrhenius lifetime fitstat fit_model→02 · extractEa, its error and L80json select→03 · extrapolateLife at 298.15 Kmath calc→04 · yearsHours to yearsunit calc

Script

admin

#1. Accelerated aging of an epoxy coating: hours to blistering at 60, 70, 80, 90 degC;
-> "[tk]" "[333.15, 343.15, 353.15, 363.15]";
-> "[life]" "[8446, 3463, 1481, 678]";
#2. Arrhenius fit of the lifetime, centred on 80 degC: L = L80.exp of Ea/R.1/T - 1/353.15;
json load "rx_ch_ar" (stat fit_model "L80*exp(Ea/8.314462618*(1/T - 1/353.15))" "T" [tk] [life] "{\"L80\": 1500, \"Ea\": 80000}");
-> "[l80]" (json select "rx_ch_ar" /parameters/L80); -> "[ea]" (json select "rx_ch_ar" /parameters/Ea); -> "[dea]" (json select "rx_ch_ar" /errors/Ea);
#3. Extrapolation to the service temperature 25 degC, in years;
-> "[l25]" (math calc "[l80]*exp([ea]/_R*(1/298.15 - 1/353.15))");
concat "Ea = " (math round_uncertainty (math calc "[ea]/1000") (math calc "[dea]/1000")) " kJ/mol; predicted life at 25 degC = " (math round (unit calc (concat [l25] " h") "yr") 1) " years"

Result

mentdb

Ea = 84.87 ± 0.17 kJ/mol; predicted life at 25 degC = 35.1 years

Go further

Use a Weibull law (stat dist_inv "weibull") to give a B10 life instead of a mean life.
Propagate the Ea error to the predicted life with num uncertainty.

executed on MentDB

Chloride ingress and corrosion initiation in concrete

Chlorides diffuse into a concrete deck following Fick's second law (erfc solution). The script computes the chloride content at the rebar after 20 years, the time to corrosion initiation in closed form with the inverse error function and with a root finder, and the cover needed for a 100-year service life.

Pipeline

01 · profileChloride at the rebarerf(x)→02 · closedInitiation with erf_invmath calc→03 · checkSame time by root findingnum solve→04 · designCover for 100 yearsmath calc

Script

admin

#1. Chloride ingress in a concrete deck: Cs = 2.5 %, C0 = 0.05 %, Ccrit = 0.40 % by binder, D = 1e-12 m2/s, cover 50 mm;
-> "[Cs]" 2.5; -> "[C0]" 0.05; -> "[Cc]" 0.40; -> "[D]" 1e-12; -> "[x]" 0.050;
#2. Fick second law: C = C0 + Cs-C0 * erfc of x / 2 sqrt Dt, chloride at the rebar after 20 years;
-> "[yr]" 31557600;
-> "[C20]" (math calc "[C0]+([Cs]-[C0])*erfc([x]/(2*sqrt([D]*20*[yr])))");
#3. Time to corrosion initiation, closed form with the inverse error function;
-> "[z]" (math calc "erf_inv(1-([Cc]-[C0])/([Cs]-[C0]))");
-> "[ti]" (math calc "([x]/(2*[z]))^2/[D]/[yr]");
#4. Check with a root finder, and the cover needed for 100 years;
-> "[ti2]" (num solve (concat [C0] "+(" [Cs] "-" [C0] ")*erfc(" [x] "/(2*sqrt(" [D] "*t*" [yr] ")))-" [Cc]) "t" 1 500);
concat "Cl at rebar after 20 y = " (math round [C20] 3) " %, initiation after " (math round [ti] 1) " y (solver " (math round [ti2] 1) " y), cover for 100 y = " (math round (math calc "2*[z]*sqrt([D]*100*[yr])*1000") 1) " mm"

Result

mentdb

Cl at rebar after 20 y = 0.44 %, initiation after 18.4 y (solver 18.4 y), cover for 100 y = 116.4 mm

Go further

Draw the profiles after 10 and 50 years with plot function and a JSON array of two expressions.
Make D decrease with time (ageing factor) and solve again with num solve.

executed on MentDB

Pourbaix diagram of iron in water

Draw the simplified E-pH diagram of the Fe-H2O system at 25 degC for 1e-6 M dissolved iron (immunity, corrosion, passivity) and place the measured state of a steel pipe on it.

Pipeline

01 · slopeNernst slope ln10.RT/Fmath calc→02 · activitylog10 of the iron activitymath log_base→03 · linesEquilibrium lines and triple pointsmath round→04 · drawE-pH diagram with the pipe pointplot chart

Script

admin

#1. Nernst slope at 25 degC and activity of the dissolved iron, 1e-6 mol/L;
-> "[s]" (math calc "ln(10)*_R*298.15/_F");
-> "[lfe]" (math log_base 0.000001 10);
#2. Boundaries in V vs SHE: Fe2+/Fe, Fe2O3/Fe3+ pH, Fe2O3/Fe2+ meets Fe2+/Fe;
-> "[ea]" (math round (math calc "-0.447+[s]/2*[lfe]") 3);
-> "[ph0]" (math round (math calc "(-0.72-[lfe])/3") 2);
-> "[ph1]" (math round (math calc "(0.728-[s]*[lfe]-[ea])/(3*[s])") 2);
-> "[e14]" (math round (math calc "[ea]-[s]*(14-[ph1])") 3);
#3. Water stability lines at pH 14;
-> "[h14]" (math round (math calc "-14*[s]") 3);
-> "[o14]" (math round (math calc "1.229-14*[s]") 3);
#4. Diagram with the measured state of a steel pipe: pH 7, -0.45 V vs SHE;
plot chart (concat "{\"title\": \"Pourbaix diagram Fe-H2O, 25 degC, 1e-6 M\", \"x_label\": \"pH\", \"y_label\": \"E (V vs SHE)\", \"series\": [{\"name\": \"Fe2+ / Fe\", \"x\": [0, " [ph1] "], \"y\": [" [ea] ", " [ea] "]}, {\"name\": \"Fe3+ / Fe2+\", \"x\": [0, " [ph0] "], \"y\": [0.771, 0.771]}, {\"name\": \"Fe2O3 / Fe3+\", \"x\": [" [ph0] ", " [ph0] "], \"y\": [0.771, 1.6]}, {\"name\": \"Fe2O3 / Fe2+\", \"x\": [" [ph0] ", " [ph1] "], \"y\": [0.771, " [ea] "]}, {\"name\": \"Fe2O3 / Fe\", \"x\": [" [ph1] ", 14], \"y\": [" [ea] ", " [e14] "]}, {\"name\": \"H2O / H2\", \"x\": [0, 14], \"y\": [0, " [h14] "]}, {\"name\": \"O2 / H2O\", \"x\": [0, 14], \"y\": [1.229, " [o14] "]}, {\"name\": \"pipe\", \"type\": \"points\", \"x\": [7], \"y\": [-0.45]}]}")

Result

mentdb

The chart is shown in the Plot tab of the editor.

Go further

Add the Fe3O4 domain (Fe3O4/Fe and Fe2O3/Fe3O4 lines) for a more complete diagram.
Redraw at 80 degC: the slope changes with T and the water lines move.

executed on MentDB

EIS fit with a CPE and inhibition efficiency

Fit the Bode modulus of inhibited mild steel with the Rs + Rct//CPE circuit, compute the Brug effective capacitance, the phase at 1 Hz from the fitted model and the inhibition efficiency against the blank.

Pipeline

01 · fitWeighted fit of the Randles-CPE modulusstat fit_model→02 · weights1 % error on each modulusvector scale→03 · capacitanceBrug effective capacitancemath power→04 · modelComplex impedance at 1 Hzcomplex log→05 · reportPhase and IE %complex arg

Script

admin

#1. Bode modulus of mild steel in 1 M HCl + 200 ppm inhibitor, f in Hz, Z in ohm.cm2;
-> "[f]" "[100000, 31622.8, 10000, 3162.28, 1000, 316.228, 100, 31.6228, 10, 3.16228, 1, 0.316228, 0.1, 0.0316228, 0.01]";
-> "[z]" "[12.5, 12.72, 13.17, 15.4, 23.64, 49.4, 121.64, 295.46, 691.24, 1252.72, 1636.23, 1787.23, 1839.02, 1857.62, 1864.64]";
#2. Fit Rs + Rct parallel to a CPE, Q in uF.s^n-1, weighted by a 1 % error on Z;
-> "[a]" "Rct*Q*1e-6*(2*pi*f)^n";
-> "[mod]" "sqrt((Rs+Rct*(1+A*cos(n*pi/2))/(1+2*A*cos(n*pi/2)+A^2))^2+(Rct*A*sin(n*pi/2)/(1+2*A*cos(n*pi/2)+A^2))^2)";
json load "rx_p4c_eis" (stat fit_model (string replace [mod] "A" (concat "(" [a] ")")) "f" [f] [z] "{\"Rs\": 10, \"Rct\": 1500, \"Q\": 50, \"n\": 0.9}" (vector scale [z] 0.01));
-> "[rs]" (json select "rx_p4c_eis" /parameters/Rs); -> "[rct]" (json select "rx_p4c_eis" /parameters/Rct);
-> "[q]" (math calc (concat (json select "rx_p4c_eis" /parameters/Q) "*1e-6")); -> "[n]" (json select "rx_p4c_eis" /parameters/n);
#3. Brug effective double layer capacitance, in uF/cm2;
-> "[ceff]" (math calc (concat (math power (math calc "[q]*(1/[rs]+1/[rct])^([n]-1)") (math calc "1/[n]")) "*1e6"));
#4. Impedance at 1 Hz with the fitted model: CPE term (jw)^n = exp of n.log jw;
-> "[jwn]" (complex exp (complex mul [n] (complex log (concat (math calc "2*pi") "i"))));
-> "[z1]" (complex add [rs] (complex div [rct] (complex add 1 (complex mul (math calc "[rct]*[q]") [jwn]))));
#5. Inhibition efficiency against the blank, Rct0 = 96.5 ohm.cm2;
concat "Rs = " (math round [rs] 2) ", Rct = " (math round [rct] 0) " ohm.cm2, n = " (math round [n] 3) ", Ceff = " (math round [ceff] 1) " uF/cm2, phase at 1 Hz = " (math round (math calc (concat (complex arg [z1]) "*180/pi")) 1) " deg, IE = " (math round (math calc "100*(1-96.5/[rct])") 1) " %"

Result

mentdb

Rs = 12.38, Rct = 1849 ohm.cm2, n = 0.836, Ceff = 8.7 uF/cm2, phase at 1 Hz = -16.4 deg, IE = 94.8 %

Go further

Fit the real and the imaginary parts separately and compare Rct.
Follow Rct with the immersion time to see the adsorption kinetics of the inhibitor.

executed on MentDB

Inhibitor efficiency statistics and Langmuir adsorption

From weight-loss replicates, test whether 50 ppm of inhibitor really reduces corrosion (non-parametric test), then fit the Langmuir isotherm on the coverage and derive the free energy of adsorption.

Pipeline

01 · testBlank vs 50 ppm, Mann-Whitney Ustat mann_whitney→02 · meansMean mass lossesmath avg→03 · coveragetheta = 1 - w/w0vector add→04 · isothermLangmuir fit, Kadsstat fit_model→05 · energydG_ads = -RT ln 55.5 Kmath calc

Script

admin

#1. Mass loss in mg of mild steel coupons, 10 cm2, 24 h in 1 M HCl at 25 degC, 5 replicates;
-> "[blank]" "[118.4, 124.9, 121.7, 116.2, 126.3]";
-> "[inh50]" "[41.2, 45.8, 38.9, 43.6, 47.1]";
#2. Is 50 ppm really protective? Mann-Whitney U test, no normality assumed;
json load "rx_p4c_mw" (stat mann_whitney [blank] [inh50]);
#3. Mean mass loss per inhibitor concentration in mol/L, M = 200 g/mol;
-> "[w0]" (math avg [blank]);
-> "[c]" "[0.00025, 0.0005, 0.001, 0.002, 0.004]";
-> "[w]" (concat "[" (math avg [inh50]) ", 25.5, 14.6, 7.9, 4.25]");
#4. Surface coverage theta = 1 - w/w0 and Langmuir isotherm theta = K.C / 1+K.C;
-> "[theta]" (vector add "[1, 1, 1, 1, 1]" (vector scale [w] (math calc "-1/[w0]")));
json load "rx_p4c_lang" (stat fit_model "K*C/(1+K*C)" "C" [c] [theta] "{\"K\": 5000}");
-> "[k]" (json select "rx_p4c_lang" /parameters/K);
#5. Standard free energy of adsorption, 55.5 mol/L is the concentration of water;
-> "[dg]" (math calc "-_R*298.15*ln(55.5*[k])/1000");
concat "CR blank = " (math round (math calc "[w0]/(10*24)") 3) " mg/cm2/h; Mann-Whitney U = " (json select "rx_p4c_mw" /u) ", p = " (math round (json select "rx_p4c_mw" /p_value) 4) "; IE at 800 ppm = " (math round (math calc "100*(1-4.25/[w0])") 1) " %; Kads = " (math round [k] 0) " L/mol, R2 = " (math round (json select "rx_p4c_lang" /r2) 4) ", dG_ads = " (math round [dg] 1) " kJ/mol"

Result

mentdb

CR blank = 0.506 mg/cm2/h; Mann-Whitney U = 25, p = 0.0079; IE at 800 ppm = 96.5 %; Kads = 7316 L/mol, R2 = 0.9994, dG_ads = -32 kJ/mol

Go further

Try the Temkin or Frumkin isotherm with stat fit_model and compare the R2.
Repeat at 35 and 45 degC: the change of Kads gives the adsorption enthalpy (van t Hoff).

executed on MentDB

Deepest pit on a pipe by Gumbel extreme values

From the deepest pit measured on small coupons, extrapolate the most probable deepest pit on 25 m2 of pipe with the Gumbel law, then estimate the time to perforate the wall with a cube-root growth law.

Pipeline

01 · momentsGumbel parameters by the momentsstat stddev→02 · extrapolateReturn period T, -ln of -ln 1-1/Tmath log1p→03 · growthDepth after 10 yearsmath cbrt→04 · lifeTime to perforationmath power

Script

admin

#1. Deepest pit in mm on each of 12 coupons of 100 cm2, carbon steel, 1 year in cooling water;
-> "[d]" "[0.82, 1.05, 0.67, 0.93, 1.21, 0.78, 0.99, 1.12, 0.71, 0.88, 1.36, 0.95]";
#2. Gumbel law of the maxima by the method of moments: alpha = s.sqrt6/pi, u = mean - 0.5772.alpha;
-> "[al]" (math calc (concat (stat stddev [d]) "*sqrt(6)/pi"));
-> "[u]" (math calc (concat (stat mean [d]) " - 0.5772156649*[al]"));
#3. Return period T = pipe area / coupon area for 25 m2 of pipe, most probable deepest pit;
-> "[t]" (math calc "25/0.01");
-> "[dmax]" (math calc (concat "[u] - [al]*ln(-(" (math log1p (math calc "-1/[t]")) "))"));
#4. Pit growth law d = d1.t^1/3: depth after 10 years and time to perforate the 6.0 mm wall;
-> "[d10]" (math calc (concat "[dmax]*" (math cbrt 10)));
-> "[life]" (math power (math calc "6.0/[dmax]") 3);
concat "Gumbel u = " (math round [u] 3) " mm, alpha = " (math round [al] 3) " mm; deepest pit on the pipe after 1 year = " (math round [dmax] 2) " mm, after 10 years = " (math round [d10] 2) " mm; perforation expected after " (math round [life] 1) " years"

Result

mentdb

Gumbel u = 0.864 mm, alpha = 0.16 mm; deepest pit on the pipe after 1 year = 2.11 mm, after 10 years = 4.56 mm; perforation expected after 22.8 years

Go further

Fit the Gumbel law by maximum likelihood or on a Gumbel probability plot with stat fit.
Give the probability to exceed a depth with the CDF exp of -exp of -(x-u)/alpha raised to T.

executed on MentDB

Galvanic couple steel-copper: Evans mixed potential

Steel bolts in a copper plate in seawater: solve the mixed potential of the couple from the polarization parameters of both metals, then get the galvanic current, the acceleration factor and the corrosion rate of the steel.

Pipeline

01 · currentsButler-Volmer currents of each metalconcat→02 · coupleMixed potential: total current = 0num solve→03 · anodeSign of the steel currentmath signum→04 · rateCorrosion rate by Faradayunit calc

Script

admin

#1. Polarization data in seawater, E in V vs SCE, i in A/cm2: steel bolts 10 cm2 in a copper plate 100 cm2;
-> "[efe]" -0.650; -> "[ife]" 0.000015; -> "[afe]" 0.060; -> "[cfe]" 0.120;
-> "[ecu]" -0.250; -> "[icu]" 0.000001; -> "[acu]" 0.050; -> "[ccu]" 0.250;
#2. Butler-Volmer current of each electrode in A as a function of E;
-> "[ife_e]" (concat "10*" [ife] "*(10^((E-(" [efe] "))/" [afe] ") - 10^(-(E-(" [efe] "))/" [cfe] "))");
-> "[icu_e]" (concat "100*" [icu] "*(10^((E-(" [ecu] "))/" [acu] ") - 10^(-(E-(" [ecu] "))/" [ccu] "))");
#3. Mixed potential of the couple: total current = 0;
-> "[eg]" (num solve (concat [ife_e] " + " [icu_e]) "E" [efe] [ecu]);
#4. Net current of the steel at Eg, its sign gives the anode, and anodic current density;
-> "[net]" (math eval [ife_e] (concat "{\"E\": " [eg] "}"));
-> "[ia]" (math calc "[ife]*10^(([eg]-[efe])/[afe])");
#5. Corrosion rate of the coupled steel by Faraday, Fe EW = 27.92 g, rho = 7.86 g/cm3;
-> "[cr]" (unit calc (concat [ia] " A/cm2 * 27.92 g/mol / (_F * 7.86 g/cm3)") "mm/yr");
concat "Ecouple = " (math round [eg] 3) " V SCE, sign of the steel current = " (math signum [net]) " so steel is the anode, galvanic current " (math round (math calc "[net]*1000") 3) " mA, anodic current " (math round (math calc "[ia]*1e6") 1) " uA/cm2 = " (math round (math calc "[ia]/[ife]") 1) " x the free corrosion, CR = " (math round [cr] 2) " mm/y"

Result

mentdb

Ecouple = -0.581 V SCE, sign of the steel current = 1 so steel is the anode, galvanic current 2.103 mA, anodic current 214.3 uA/cm2 = 14.3 x the free corrosion, CR = 2.49 mm/y

Go further

Change the area ratio to show why a small anode in a large cathode is the worst case.
Draw the Evans diagram with plot function on log10 of the absolute currents.