Laplace Transform and Bessel Functions in Solving an Electrochemical Corrosion Problem

Electrochemical corrosion protection of metallic structures, particularly, in oil storage tanks, is a problem of great industrial importance. This paper proposes a mathematical model based on partial differential equations, which can be simplified to an ordinary differential equation with variable coefficients, under certain conditions. The latter is solved using Laplace transform, resulting in a solution expressed in terms of modified Bessel functions.

Introduction

A problem of great importance to industry—specifically the oil industry—is the protection of its metallic structures against electrochemical corrosion. In recent literature, various approaches to solve this problem can be found. These include the use of corrosion inhibitors [1-3], including those derived from natural extracts such as Caesalpinia spinosa and Ilex guayusa, which have shown promising results in mitigating steel corrosion [3]. Other strategies involve constructing structures from specially engineered alloys with enhanced resistance to corrosion, such as nitrogen-alloyed VCoNi systems developed for tribocorrosion conditions [4]. Additionally, the use of smart materials has gained attention for their self-healing and responsive capabilities, offering innovative possibilities in corrosion protection [5].

One of the most widely used techniques to achieve this objective is the use of sacrificial anodes, which must be placed in precise locations to maximize protection and minimize installation costs.

This paper examines the corrosion protection of ANCAP’s oil storage tank roofs (the state-owned oil company of Uruguay), which are located at the Eastern Terminal, near Punta José Ignacio, a cape on the Uruguayan Atlantic coast [6].

Although oil does not conduct ions –therefore it is not a corrosive medium– it should be noted that crude oil is typically accompanied by water, which accumulates at the bottom of the aforementioned storage tanks. The design of anodic protection involves choosing the appropriate material, calculating their mass, the position, and number of anodes to be used. This paper summarizes the key steps of the mathematical modeling of the problem. This model corresponds to an elliptic PDE in cylindrical coordinates, accompanied by nonlinear boundary conditions [7,8]. After several simplifications, it results in a second-order ODE with variable coefficients.

The simplified equation is solved by Laplace Transform, resulting in a general solution, expressed as a linear combination of Bessel functions plus a constant. This solution will serve as the fundamental input for an analytical-numerical method, in which the radii at which the sacrificial anodes will be located are calculated using an implicit iterative method, to be developed in further research.

The mathematical model

In this section, we will succinctly describe how the mathematical model used to determine the distribution of the anodes within the tank was obtained.

Firstly, a mass balance is performed in a control volume [9] within the electrolyte. In the case of species i this balance gives:

$\frac{\partial c_i}{\partial t}=-\nabla \cdot J_i+R_i\text{ (1)}$

Similarly, the overall balance for all species involved is:

$\frac{\partial {{\displaystyle \sum }}_i{c}_i}{\partial t}=-\nabla \cdot \left({{\displaystyle \sum }}_i{J}_i\right)+{{\displaystyle \sum }}_i{R}_i\text{ (2)}$

Equation (2) is transformed into a charge balance by multiplying by the charges of the species z_i and the Faraday constant F. This charge balance is presented in Equation 3.

$F \frac{\partial {{\displaystyle \sum }}_i{z}_{i }{c}_i}{\partial t}=-\nabla \cdot \left(F{{\displaystyle \sum }}_i{z}_{i }{J}_i\right)+F{{\displaystyle \sum }}_i{z}_{i }{R}_i\text{ (3)}$

Since in this case we study an electrochemical system with reaction at the interface, it can be assumed that there is no generation in the control volume (then R_i = 0). Also, $F \frac{\partial {{\displaystyle \sum }}_i{z}_{i }{c}_i}{\partial t}=0$ due to the steady-state condition, and so, we obtain:

$\nabla \cdot \left(F{{\displaystyle \sum }}_i{z}_{i }{J}_i\right)=\nabla \cdot \overrightarrow{j}=0\text{ (4)}$

where $\overrightarrow{j}$ is the current density vector.

Secondly, $\overrightarrow{j}$ is developed considering migration, diffusion, and advection [10], which initially results in a three-summand equation [8], like in the following equation:

$\nabla \cdot \overrightarrow{j}=-\nabla \cdot \left({{\displaystyle \sum }}_i \left|z_i\right| u_{i }{c}_{i }\nabla E\right)-\nabla \cdot \left(F {{\displaystyle \sum }}_i D_{i }{z}_{i }\nabla c_i\right)+\nabla \cdot \left(Fv{{\displaystyle \sum }}_i z_{i }{c}_{i }\right)=0\text{ (5)}$

Where v is the local velocity of the fluid, D_i is the diffusivity of species i, and u_i is the ionic mobility.

Equation (5) is simplified under the assumptions that the solution is an isotropic medium (so the ionic conductivity and diffusivity of the species have zero gradients) and applying the electroneutrality condition (thus the convective term vanishes). Because of this, Equation (5) further simplifies to:

$\nabla \cdot \overrightarrow{j}=-\text{ χ }{\nabla }^2\text{E }-F{{\displaystyle \sum }}_i{D}_{i }{z}_{i }{\nabla }^2{c}_i+\nabla \cdot \left(Fv{{\displaystyle \sum }}_i z_{i }{c}_{i }\right)=0\text{ (6)}$

Where x is the electrolyte conductivity?

Due to the condition of electroneutrality, the convective term is zero, and then, Equation (6) can be simplified one more time to give:

$\nabla \cdot \overrightarrow{j}=-\text{ χ }{\nabla }^2\text{E }-F{{\displaystyle \sum }}_i{D}_{i }{z}_{i }{\nabla }^2{c}_i=0\text{ (7)}$

Where x is the electrolyte conductivity?

Due to the condition of electroneutrality, the convective term is zero, and then, Equation (6) can be simplified one more time to give:

$\nabla \cdot \overrightarrow{j}=-\text{ χ }{\nabla }^2\text{E }=0\text{ (8)}$

Equation (8) gives rise to the following elliptic PDE:

${\nabla }^{2}E=0\text{ (9)}$

Since our system is a cylindrical tank with D = 60 m and h = 0.015 m (average height of the conductive electrolyte present, the heigh tank is 20 m), it is more convenient to use cylindrical coordinates [11], and we obtain:

$\frac{{\partial }^{2}E}{\partial r^2}+\frac{1}{r} \frac{\partial E}{\partial r}+\frac{1}{r^2} \frac{{\partial }^{2}E}{\partial {\theta }^2}+\frac{{\partial }^{2}E}{\partial z^2}=0\text{ (10)}$

Due to symmetry reasons, the derivatives with respect to the polar angle θ can be canceled, obtaining:

$\frac{{\partial }^{2}E}{\partial r^2}+\frac{1}{r} \frac{\partial E}{\partial r}+\frac{{\partial }^{2}E}{\partial z^2}=0\text{ (11)}$

Furthermore, for symmetry reasons, there is no flow in the center of the tank; therefore,

${\frac{\partial E}{\partial r}|}_{r=0}=0\text{ (12)}$

The other boundary condition is given by the experimental polarization curves (electrochemical potential vs. current density), obtained in the laboratory and piecewise linearized.

Finally, a new boundary condition is obtained at the oil/water interface, where there is no flow, and therefore,

$\frac{\partial E}{\partial r}=0$ Evaluated at z = L (13)

In this system, the Frumkin condition for unidirectional flow is satisfied [12] for unidirectional flow is satisfied, which simplifies the differential Equation (11) to give

$\frac{{\partial }^{2}E}{\partial r^2}+\frac{1}{r} \frac{\partial E}{\partial r}=0\text{ (14)}$

If Equation (14) is written in terms of the electrochemical potential ( $\overline{\mu }$ ), the equation results:

$\frac{{\partial }^2\overline{\mu }}{\partial r^2}+\frac{1}{r} \frac{\partial \overline{\mu }}{\partial r}=\pm i_S \frac{2\rho }{L}\text{ (15)}$

being i_s the normal current to the metallic surface.

The boundary condition can be piecewise linearized, using appropriate software, or manually. The only condition is that the curve can be locally considered as a straight line segment with an approximation error below the predetermined error margins.

The result of this process is given by

${\frac{\partial \overline{\mu }}{\partial z}|}_{z=0}=i_S=a+b\overline{\mu }\text{ (16)}$

In Equation (16), a and b are the parameters corresponding to the piecewise linearization of the boundary condition.

After a change of variable x = r / R, Equation (15) becomes:

$\frac{{\partial }^2\overline{\mu }}{\partial x^2}+\frac{1}{x} \frac{\partial \overline{\mu }}{\partial x}=\pm \left(\frac{a}{b}+\overline{\mu }\right) \frac{b{R}^{2}2\rho }{L}\text{ (17)}$

Equation (17) is of the form

$y^{\text{'}\text{'}}+ \frac{1}{x} y^{\text{'}}-Ky=AK\text{ (18)}$

Where $\frac{b{R}^{2}2\rho }{L}=\text{K}$ and $\frac{a}{b}=A$ depend on the piecewise linearization of the experimental polarization curves, the geometric characteristics, and the conductivity of the medium.

The general solution

The homogeneous equation associated with Equation (18) is:

$y^{″}+ \frac{1}{x} y^{\prime }-Ky=0\text{ (19)}$

If we put K = a² in equation (19), and multiplying by x we obtain:

$x{y}^{″}+ y^{\prime }-a^{2}xy=0\text{ (20)}$

Let be Y(s) the Laplace transform of y(x), then:

$ℒ\left\{y\left(x\right)\right\}=Y\left(s\right),ℒ\left\{y\text{'}\left(x\right)\right\}=sY\left(s\right)-\alpha ,ℒ\left\{y^{″}\left(x\right)\right\}=s^{2}Y\left(s\right)-\alpha s-\beta \text{ (21)}$

Notably, since x = r / R, then x = ∈ [0,1], therefore, the functions y(x), y′ (x), y′′ (x) possess compact support, thereby ensuring convergence of their Laplace transforms .

Therefore, the transformed equation becomes:

$ℒ\left\{x y\left(x\right)\right\}=-\frac{d}{ds} ℒ\left\{y\left(x\right)\right\}= -Y^{\prime }\left(s\right)\text{ (22)}$

And,

$ℒ\left\{x y^{″}\left(x\right)\right\}=-\frac{d}{ds} ℒ\left\{{y^{″}}^{\left(x\right)}\right\}= -s^2{Y}^{\prime }\left(s\right)-2s Y\left(s\right)+\alpha \text{ (23)}$

Then by applying Laplace transform, the ODE (20) becomes:

$\left(a^2-s^2\right)Y^{\prime }\left(s\right)-s Y\left(s\right)=0\text{ (24)}$

Separating variables and integrating we obtain:

$Y\left(s\right)=C.\frac{1}{\sqrt{s^2-a^2}}\text{ (25)}$

the inverse Laplace transform of which is C.I_o (a x), where I_o(a x) is the modified Bessel function of the first kind, with order zero [13].

From Equation (25) it is easy to obtain the general solution of the homogeneous ODE as:

$c_1.I_0\left(a x\right)+c_2.K_0\left(a x\right)\text{ (26)}$

where I_o(x) and K_o x) are modified Bessel functions of order zero.

It can be shown that y(x) = −A is a particular solution of the non-homogeneous ODE (18) and therefore the general solution of (18) is of the form:

$y\left(x\right)=c_1.I_0\left(a x\right)+c_2.K_0\left(a x\right)-A\text{ (27)}$

Some considerations on the use of the general solution

To determine the optimal electrode placement radii for placing the sacrificial electrodes, it is necessary to divide the domain of variable x = r / R into several intervals

$\left[x_0 , x_1\right],\left[x_1 , x_2\right] , \dots , \left[x_{n-1} , x_n\right]\text{ (28)}$

where the coefficients of the ODE (18), K and A take different values, depending on the piecewise linearization of the experimental polarization curves.

As a consequence, the constants c₁, c₂ of the general solution, take also different values in each interval.

Considering that there are more unknowns than equations, it is necessary to add other equalities, which can arise from deriving (27), to obtain:

$y^{\prime }\left(x\right)=c_1.a.I_1\left(a x\right)-c_2.a.K_1\left(a x\right)\text{ (29)}$

since

$\frac{d}{dx}{I}_0\left( x\right)=I_1\left( x\right)$ and [13] $\frac{d}{dx}{K}_0\left( x\right)=-K_1\left( x\right)$ (30)

The derivative at the internal nodes (x_i, i = 1, 2, … , n − 1) can be equated to the bilateral numerical derivative to obtain more equations, which would allow obtaining an analytical-numerical method to find the positions of the electrodes.

However, if this iterative process is performed forward, it becomes unstable. Therefore, it is better to use these equations to obtain a more complicated, but more stable, implicit, backward method.

Conclusion

In this work, a model was developed for calculating the positions of the sacrificial electrodes for tank protection in the petroleum industry. The proposed equation admits an analytical solution using the Laplace transform and modified Bessel functions. This improves previous results [8], presented in the form of power series, which are not easy to implement in the solution of experimental problems.

The proposed solution is both simpler and more elegant; however, the number of unknowns, which exceeds the number of equations, requires the addition of new numerical-analytical equalities.

Since forward methods proved unstable, in order to use those equations, the iterative method must be a backward analytical-numerical method. This implicit backward method must be carried out by dividing the interval (x = ∈ [0,1]) into subintervals, where the obtained analytical solution and its derivative will be applied. The analytical derivative will be matched to the bilateral numerical derivative in all the internal nodes of the partition.

The proposed method could be applied to other cylindrical geometries where the unidirectional flow condition is satisfied, making it readily adaptable to similar cylindrical corrosion problems.

The authors wish to thank the anonymous reviewers for their helpful suggestions for improving the chapter and also thank Marjorie Chaves for her contribution in reviewing this manuscript.

1. Hidalgo J, Hidalgo L, Serrano Aguiar CD, García Madroñero DB, Galambos I, Vilasó-Cadre JE, et al. Study of Caesalpinia spinosa extracts as green corrosion inhibitor for mild steel. Langmuir. 2025;41(14):9406–21. Available from: https://doi.org/10.1021/acs.langmuir.5c00224
2. Hidalgo J, Cotolan N, Lupan A, Brânzanic AMV, Turdean GL. A poly(methyl methacrylate)–ibuprofen composite film as anticorrosive coating of Ti–6Al–4V surface. J Solid State Electrochem. 2024;28:479–94. Available from: https://doi.org/10.1007/s10008-023-05681-w
3. Hidalgo J, Hidalgo L, Serrano C, Punina D, Rea E, Ilbay M, et al. A study of the inhibition capacity of a novel Ilex guayusa green extract for preventing corrosion in mild steel exposed to different conditions. Analytica. 2025;6(1). Available from: https://doi.org/10.3390/analytica6010001
4. Jiao Z, Dong Y, Li Q, Zhou Q. Enhancing tribocorrosion resistance of VCoNi alloys in artificial seawater via nitrogen alloying. Corros Sci. 2025;243:112600. Available from: http://dx.doi.org/10.1016/j.corsci.2024.112600
5. Siddiqui NA, Muzamil M, Jamil T, Hussain G. Smart materials in manufacturing. Smart Mater Manuf. 2025;3:100080. Available from: https://www.keaipublishing.com/en/journals/smart-materials-in-manufacturing/
6. Campo L, Martinez-Luaces V, Ohanian M. Sistema de protección catódica interna en tanques de crudo. Montevideo (UY): ANCAP-UdelaR; 2017. Technical report.
7. Ohanian M, Martinez-Luaces V. Corrosion potential profile simulation in a tube under cathodic protection. Int J Corros. 2014;2014:Article ID 102363. Available from: http://dx.doi.org/10.1155/2014/102363
8. Martinez-Luaces V, Campo L, Ohanian M. Calculus of the distribution of sacrificial anodes for roof oil tanks cathodic corrosion protection. Adv Math Model Appl. 2019;4(3):224–31. Available from: https://jomardpublishing.com/UploadFiles/Files/journals/AMMAV1N1/V4N3/Martinez-Luaces.pdf
9. Ibl N. Current distribution. In: Yeager E, Bockris JO, Conway BE, Sarangapani S, editors. Comprehensive Treatise of Electrochemistry. Boston (USA): Springer. 1983;6:239–315. Available from: https://link.springer.com/chapter/10.1007/978-1-4615-6690-8_4
10. Newman J. Engineering design of electrochemical systems. Ind Eng Chem. 1968;60(4):12–27. Available from: https://pubs.acs.org/doi/10.1021/ie50700a005
11. Zill DG, Cullen MR. Differential equations with boundary-value problems. Belmont (USA): Cengage Learning; 2009. p. 480–1. Available from: http://www.uop.edu.pk/ocontents/DIFFERENTIAL%20EQUATIONS%20with%20Boundary-Value%20Problemsa%20Zill%20Cullen(1).pdf
12. Frumkin AN. Distribution of the corrosion process along the tube length. J Phys Chem. 1949;23:1477–82.
13. Spiegel M, Lipschutz S, Liu J. Schaum's outline of mathematical handbook of formulas and tables. 3rd ed. London (UK): McGraw-Hill Education; 2008. Available from: http://mcb111.org/w01/Mathematical_Handbook_of_Formulas_and_Ta.pdf