Basic models and approximation for the engineering description of the kinetics of the oxide layer of steel in a flow of heavy liquid metal coolant under various oxygen conditions
expand article infoAlexandr V. Avdeenkov, Oleg I. Achakovsky, Vladimir V. Ketlerov, Vladimir Ya. Kumaev, Alexander I. Orlov§
‡ JSC SSC RF-IPPE, Obninsk, Russia
§ Proryv JSC, Moscow, Russia
Open Access


The article presents the results of corrosion processes, kinetics and changes in the oxide layer modeling using MASKA-LM software complex. The complex is intended for a numerical simulation of three-dimensional non-stationary processes of mass transfer and interaction of impurity components in a heavy liquid metal coolant (HLMC: lead, lead-bismuth). The software complex is based on the numerical solution of coupled three-dimensional equations of hydrodynamics, heat transfer, formation and convective-diffusive transport of chemically interacting components of impurities.

Examples of calculations of mass transfer processes and interaction of impurity components in HLMC, formation of protective oxide films on the surfaces of steels are given to justify the coolant technology.


coolant, corrosion, films, HLMC, heavy liquid metal coolant, impurities, lead, lead-bismuth, mass transfer


With the advent of liquid metal coolants necessary for nuclear power plants with fast neutron reactors, the key elements of the safety justification were attributed to the problems associated with maintaining the quality of the coolant. Maintaining the necessary physicochemical regime of the coolant (or coolant technology) has turned from an important but additional operation into one of the main components in substantiating the design and operating parameters of the reactor. There are two main tasks of the heavy liquid metal coolant (lead or lead-bismuth) technology:

  • ensuring the cleanliness of the coolant and the surfaces of the circulation circuit to maintain the design thermohydraulic characteristics with long service life (several decades when the reactor plant operates at a power of up to 100%);
  • prevention of corrosion and erosion of structural materials with long service life (several decades, when the reactor plant is operating at a power of up to 100%).

It should be noted that the main processes of the lead-bismuth and lead coolant technology have a common physicochemical nature. In the lead-bismuth eutectic, lead exhibits the dominant chemical activity. Considering also that the construction materials for both coolants are identical, the methods and means of the coolant technology are also similar in the principle of operation.

Therefore understanding the process of generation of the oxide layer of steel in a flow of heavy liquid metal coolant is highly important for justification of circuit integrity and reactor safety.

Kinetics of corrosion processes of steel components in HLMC is determined by the associated processes of hydrodynamics and the interaction of dissolved components, primarily such as oxygen and iron. The justification and appropriate modeling of these processes is an important component for substantiating heavy liquid metal coolant technology (Robertson et al. 1988; Fazio et al. 2001; Aiello et al. 2004; Ivanov et al. 2005, 2017; Zhang et al. 2005, 2007; Lim et al. 2010; Martinelli et al. 2008a; Weisenburger 2008; Schroer and Konys 2009; Steiner 2009; Hwang and Lim 2010; Askhadullin et al. 2011; Schroer et al. 2012; Zhang 2013; Nuclear Energy Agency 2015; Tsisar et al. 2017; Salaev et al. 2018).

The MASKA-LM software complex implements basic models and approximations with the following features:

  • modeling of the distribution of dissolved impurities along the circulation path of lead coolant in various operating modes of a reactor unit, taking into account the condition of the surfaces of structural steels;
  • simulation of a transient distribution of impurities when changing operating modes;
  • predicting the deviation of the quality indicators of lead coolant from normal operating conditions for various initiating events;
  • simulation of the distribution of the gas phase in lead coolant;
  • simulation of the transfer of solid-phase impurities along the circulation circuit.

The use of commercial CFD codes in this area has not yet been widely applied because of the need to build a complex apparatus of user functions for coordinated calculations of thermal hydraulics and processes of interaction and formation of impurities. For example, the works (Marino 2015; Marino et al. 2018) on modeling physicochemical processes in lead coolant show the possibility of using this approach. However, even for the simplest cases (individual sections, areas, etc.), the dimension of computational grids reaches excessively large values ​​at which the use of only a CFD code will require significant machine resources. At the same time, the general accuracy of the calculations is mainly determined by the accuracy of the equations of physicochemical processes. "Massive" calculations of three-dimensional non-stationary distributions of impurity components in the primary circuit in lead coolant at "life cycle" mode using a CFD code with fine grids are problematic due to unreasonable calculation time, and the use of coarse grids eliminates the "hydrodynamic" advantages of CFD.

Mathematical model of the processes of mass transfer and physico-chemical interaction of impurities in the volume of the coolant and on the surfaces of structural steels

A system of three-dimensional equations of an incompressible multicomponent medium in a Cartesian coordinate system was used to describe the processes of formation, transformation, and transfer of impurities in the primary circuit of a reactor with HLMC (Alekseev et al. 2010a). The system includes hydrodynamic equations in the Boussinesq approximation, equations of conservation and transfer of thermal energy, equations of physicochemical kinetics of impurities (Kumaev et al. 2002), equations of conservation and transfer of components of a multicomponent medium. The system of equations is solved by finite-difference methods.

The medium is assumed to be incompressible with temperature-dependent thermophysical properties. The presence of dissolved and dispersed particulate impurities does not affect the hydrodynamics, heat transfer, and thermophysical properties of lead. Concentration distributions of dissolved impurities and particles are affected by molecular, Brownian, and turbulent diffusion. On the surfaces of structural materials, oxide films are formed that impede the diffusion exchange of oxygen and iron between the coolant and structural materials.

The system of equations describing the non-stationary thermohydrodynamics of the components of a multicomponent medium is presented in the form of equations of motion, conservation of mass and thermal energy of an inhomogeneous continuous medium (Kumaev et al. 2002, 2005; Alekseev et al. 2010a, b, c). The components of the medium are fuel, structural and other materials that are in a solid state and do not move in space, liquid media that can move, as well as components of various impurities that undergo transformations in physicochemical reactions. The viscosity and thermal conductivity coefficients of a liquid are presented as the sum of molecular and turbulent components.

The equations of motion for the characteristics of a multicomponent medium averaged over components that are not impurities were obtained in (Kumaev et al. 2002, 2005) and in the Cartesian coordinate system in the Boussinesq approximation and have the following form:

Uiτ+UiUkxk=-1ρpxi+1ρxkμUixk+giρ~(t)ρ-ΛikUk , (1)


where: Ui – projections (components) of the medium velocity vector in the Cartesian coordinate system; p – pressure; ρ – density of medium averaged over temperature and components; μ – dynamic effective viscosity coefficient; ῀ρ (t) – averaged only over components and temperature-dependent density of the medium; gi – projections of the gravitational acceleration vector on the coordinate axes; Λik – tensor of drag coefficients.

The drag coefficients of the medium are also used to simulate the solid and liquid states of the medium. Zero values of the resistance coefficients correspond to the liquid state of the components; large values (tending to infinity) correspond to the solid state of the components.

Environmental properties are obtained by averaging the initial properties of the components and are determined by the following dependencies:

ρ=nρ0nεn ; where ρn0 – density of n component;

μ=nμ0nεn ; where μn0 – n component dynamic viscosity coefficient;

Uk=nuknεn – projections of the velocity vector of the averaged motion of the medium;

εn – volume fraction of component with number n;

ukn=Uk+Δukn – projections of the velocity vector of the medium component with number n.

The drag coefficients of the medium Λik are used to simulate the solid and liquid states of the medium. For components in the liquid state, Λik is set equal to zero, for components in the solid state it is set to a larger value, which ensures that the velocity of the medium is equal to zero.

The equations of transport of non-impurity components of the medium, taking into account stratification and taking into account the sources of components that can occur during chemical reactions, have the form:

εnτ+xkεnUk=-xkεfnΔukn+p=1NJnp , (2)

For components that are impurities, diffusion processes are additionally taken into account in the transport equations:

Cnτ+xkCnUk=-xkCfnΔukn+xkDnxkCn+p=1NGnp , (3)

where Δukn – components of the velocity vectors of the relative motion of the components of the medium (stratification velocity);

ε fn – volume fraction of component n, that is in a liquid state;

Jnp – volume source of component n due to loss of component р;

N fn – mass concentration of the impurity component n, that is in the liquid;

Gnp – mass source of the impurity component n due to a decrease in the impurity component р;

Dn – component n diffusion coefficient.

The left-hand sides of equations (2.3) describe the transfer of medium components due to averaged motion. The right-hand sides of equations (2) describe the process of separation of a mixture of components in a liquid state, sources of components due to chemical reactions or other processes, as well as diffusion processes of component transfer. Solid components do not delaminate. The initial and boundary conditions that describe the flows of impurities onto solid surfaces, as well as the values of the concentrations of impurities at the input, are added to the transport equation.

Modeling of heat transfer and heat exchange processes in a multicomponent medium taking into account stratification is carried out by solving the energy equation in the following form:

(cρt)τ+cρUktxk+nxkcnρnΔuknεfnt=xkλtxk+qv , (4)


с – the specific heat of the medium averaged over the components;

λ – thermal conductivity coefficient averaged over the components of the medium;

Uk – medium volume velocity;

Δukn – projections of the component stratification velocity vector n;

ε fn – volume fraction of a component n in a liquid state.

The first two terms on the left side of equation (4) describe the convective transfer of thermal energy. The third term is the transfer of thermal energy due to the separation of components.

Description of the processes of physicochemical interaction of impurities in the coolant volume and on the surfaces of structural steels

Mass transfer processes in nonisothermal circuits with HLMC are closely related to the interaction of the coolant with structural steels. Structural steels are the main sources of metallic impurities entering the coolant (iron, chromium, nickel, etc.). In turn, the coolant is a supplier of dissolved oxygen for structural steels, which is necessary for the formation of protective oxide films on steel surfaces that prevent the development of corrosion processes.

Dissolved oxygen in the coolant is both in free form, and in the form of various compounds with a coolant and metal impurities. Moreover, in quantitative terms, free oxygen in the solution is negligible relative to the bound forms of dissolved oxygen.

The main metal impurity, which should be taken into account when calculating the simulation of mass transfer processes in the circulation circuit with HLMC, is an iron impurity. In lead melt, iron can be both free and in oxygen-bound form. The relationship between these forms of existence of dissolved iron in the coolant depends on the oxygen content and can vary widely.

In the model of interaction of dissolved impurities in the coolant, only those forms of existence of oxide compounds that play a decisive role should be taken into account: PbO, Fe3O4.

The most important in modeling physicochemical interaction of impurities are the processes of formation and growth of protective oxide films on the surfaces of structural steels. The composition of oxide film, as a rule, varies from chrome spinel near the boundary with the coolant to chromium spinel at the border with steel. In this work, only the process of magnetite formation at the "coolant – oxide film" interface was taken into account, and the formation of chromium spinel is not modeled separately, but is completely determined by the model of magnetite film formation. This approach can be considered as simplified, but it gives quite adequate description of the evolution of the oxide film.

Oxide film growth kinetics model

Oxidative kinetics leading to the parabolic law of growth of an oxide film is usually used when the effective diffusion of reagents is independent of time and size of the film. The oxidation mechanism of ferritic-martensitic steels is determined by the diffusion of Fe ions in the direction of the coolant and variations in the oxygen potential. In this approach, oxidation of the surface of steels is described by the parabolic velocity equation with the rate constant predicted by Wagner’s theory if the resulting oxide film has an "ideal" crystal lattice, and ion diffusion is dominant.

General Provisions

For simplicity, we consider the growth of a single-layer film. In the developed (Kumaev et al. 2002, 2005; Alekseev et al. 2010c) model of the growth of a single-layer film, at the first stage, the iron flux "remaining" in the wall layer is determined as the difference between the fluxes through the oxide film and the flux that has gone into the coolant. The "oncoming" oxygen flow from the coolant is determined according to the stoichiometry of the resulting magnetite.

In this approximation, it is unambiguously assumed that the growth of the film (magnetite) occurs due to the reaction on the surface of the oxide layer:

3Fe + 4PbO ⇔ Fe3O4 + 4Pb, (5)

and at the same time, all four components of the reaction are in chemical equilibrium. Since the rate of chemical reactions is much higher than the characteristic rates of all other physical processes, the chemical equilibrium can be considered steady in the volume of the coolant at each moment of time and the distribution of components substantially depends on thermohydraulic processes.

Let us consider a model assuming the presence of a diffusion "reaction" layer at the interface between the oxide film and the coolant. The main components interacting with oxygen are iron and lead. We consider the oxide film as a single layer and consisting of magnetite. At the steel-coolant boundary, we distinguish four subregions: steel (conventionally Fe), oxide, an intermediate diffusion layer, and HLMC (Fig. 1).

Figure 1. 

Layers near the interface.

Intermediate layer is a virtual layer where the above-described surface reaction is carried out (5). For flows of components and sources of impurities in the media presented in Fig. 1, we use the notation: J βO – O (PbO) oxygen flow in the coolant; J βFe – Fe oxygen flow in the coolant; J δO – oxygen flow in oxide film; J δFe – iron flow in oxide film; JOγ~vOγt – oxygen flow (source) in the intermediate layer; JFeγ~vFeγt – iron flow (source) in the intermediate layer; v γO,Fe – the amount of oxygen and iron in the intermediate layer, for impurity concentrations , where n = O, PbO, Fe, Fe3O4 and activity an=CnCnS , C Sn – saturation concentration.

The chemical equilibrium of the oxygen-lead system with the formation of the [PbO] complex is described in the standard way:

kβox=aPbOaPbaOβ=exp-ΔGPbORT (6)

Lead activity is equal to one; therefore, this ratio shows the ratio of oxygen concentrations and the dynamic complex of PbO. For the liquid phase, the existence of a separate oxygen phase [O] does not have deep physical meaning, and its concentration is much lower than the concentration of [PbO]. The concentration of [O] can be considered as the concentration of "instantly free" oxygen.

Using the data (Nuclear Energy Agency 2015) for the Gibbs energy parameters for the formation of lead oxide and relation (6), it is easy to estimate that aPbO / aO ~108, so for oxygen activity it is enough to understand the activity of the complex [PbO].

For HLMC medium with impurities (Pb, PbO, Fe, Fe3O4) we introduce the notation β. In the presence of iron in the coolant, three reactions are formally possible:

Pb + O = PbO

3Fe + 4O = Fe3O4 (7)

3Fe + 4PbO = Fe3O4 + 4Pb

Chemical equilibrium constant for the second reaction:

kβO=1aFeβ3a0β4=exp-ΔGFe3O4RT (8)

The chemical equilibrium constant for the third reaction, provided that the activity of magnetite is equal to unity:

kβm=aFeβ3a0β4=exp-4ΔGPbO-ΔGFe3O4RT (9)

where anβ=CnβCnS – activities corresponding to the coolant component.

The equilibrium constants (6, 8, 9) in the dynamic and chemical equilibrium of all components are obviously related:

kβm=kβox4/kβO (10)

Obviously, in conditions of dynamic equilibrium, it is sufficient to consider any two of the three equations. From the point of view of the thermodynamics of the lead-iron-oxygen system, it is more convenient to consider the first two equations. And the area of existence of iron oxides of various compositions in terms of oxygen potentials overlaps with the region of existence of oxygen solutions in the lead melt.

Let’s consider the surface layer at the media interface of oxide film-coolant (medium γ). As we discussed above, if we formally approach this region as well as the coolant region, then there will be no differences between them and, moreover, there will not even be a need to introduce this region for consideration. Since a possible film growth is assumed in this medium, that is, a kind of phase transition:

Fe3O4 (solid phase) + Pb (liquid phase) →

dissolved (Fe, PbO) in Pb (liquid phase), (11)

then in this environment we assume the possibility of only the third reaction from (11).

Equilibrium constant in this environment:

kγ=aFeγ3aOγ4=exp-4ΔGPbO-ΔGFe3O4RT (12)

Once again, we note that in general case, knowledge of the equilibrium constant is not enough for this medium, since we will assume that the existence of equilibrium is not a prerequisite. In this case, it is necessary to know the parameters of forward and reverse reactions (7). In addition, the third reaction (7), which occurs on the surface of magnetite, is in reality not "one-step" as it is written, but contains steps for adsorption and desorption of the elements involved. Therefore, the reaction constants do not have to be identical with the constants in a homogeneous medium. The oxygen-coolant system is assumed in chemical equilibrium, similarly to formula (6). It should also be noted that the heterogeneous reaction rate is directly proportional to the area of contacting reagents. Therefore, "rough" surface, obviously, will provide a more intensive reaction process. Therefore, oxygen consumption by the steel surface will depend not only on its type, but also on the method of its manufacture. The latter is almost impossible to take into account when modeling, which can significantly reduce the physical justification for the direct application of the Wagner model (Robertson et al. 1988).

Let’s consider the layer of the film itself, which at this stage is assumed to be composed only of magnetite Fe3O4 (δ medium). In general, film growth also occurs on the metal-magnetite interface, but since the oxygen diffusion rate is much lower than the rate of the chemical reaction of incorporation of oxygen into the magnetite lattice, it is sufficient to know the concentration at the metal-magnetite interface to calculate the oxygen flux.

The constant of chemical equilibrium in this medium, taking into account the fact that the activity of magnetite and iron are equal to unity:

kOδ=1aOδ4=exp-ΔGFe3O4RTδ , (13)

where a αO – oxygen activity in the medium δ, Τ δ – film temperature.

The condition of chemical equilibrium in the intermediate layer during the formation of magnetite:

νOγνFeγ=νOγνFeγ=43=JOγJFeγ (14)

The formulas for the sources of oxygen and iron in the surface layer γ, J γO and J γFe, are proportional to the growth rate (or degradation) of this "reaction" layer, which in turn is proportional to the amount of oxygen and iron substance in it. Since the system in question is assumed in chemical equilibrium, the velocities of direct, R = K (C γFe)3 (C γO=PbO)4 and reverse, R = K (CPb)4CO=Fe3O4 reactions are the same and from the chemical equilibrium condition (12) it is possible to determine only the ratio of concentrations and "fluxes" of iron and oxygen. That is, the condition of chemical equilibrium and the assumption of the instantaneous establishment of this chemical equilibrium in comparison with the slow diffusion process does not make it possible to determine the slow growth of the film and additional assumptions about the mechanism of film growth are needed, namely, that the molar fluxes of oxygen and iron are treated as 4/3, in accordance with the kinetics of reaction (7), with instant subsequent formation of magnetite on the surface.

The source of oxygen in the system is the coolant (oxygen content is regulated through external sources). Therefore, "diffusion" of oxygen from the coolant to the steel-coolant surface section and the formation of a conventional oxide film (magnetite in our model) are initially assumed. That is, oxygen is "consumed" on the growth of the inner and outer parts of the film, respectively:

J0βx=δ+γ=J0γ+J0δδ (15)

Similarly for iron that diffuses from steel:

JFeδδ=JFeγ+JFeβx=δ+γ (16)

Note that the model considered in (Alekseev et al. 2010b, c) is a special case of the proposed model. According to (Alekseev et al. 2010b, c), iron fluxes at the interface can be calculated as follows:

JFeδδ=DFeδCFeδ-CFeγδ (17)

JFeβx=δ+γ=χFeCFeγ-CFeβ , (18)

where χFe – mass transfer coefficient of iron at the border.

Similarly, oxygen flows in media can be represented by δ (Alekseev et al. 2010b, c) and β respectively:

JOδδ=2DOδCOδ-COγδ (19)

JOβx=δ+γ=χOCOβ-COγ , (20)

where χO – oxygen mass transfer coefficient.

Using ratios (16) and (17) we get:

θJFeδ-JFeβ=JOδ-JOβ , (21)

where θ=434MO3MFe – for molar (mass) flows, respectively, or

4MO3MFeDFeδCFeδ-CFeγδ-χFeCFeγ-CFeβ=χOCOβ-COγ-2DOδCOδ-COγδ , (22)

where JO,Feβ,δ ‒ flows of oxygen and iron in the coolant stream and the film, respectively, C δFe – iron concentration in the layer δ, C βFe – iron concentration in the coolant flow, C γFe – iron concentration in the layer γ, C βO – oxygen concentration in the coolant flow, C δO – oxygen concentration in the layer δ.

Using relations (17–20) and the relationship between the activities of oxygen and iron in the wall layer γ, we can obtain a nonlinear equation for the concentration (or activity) of oxygen in the wall layer in which the concentrations of iron and oxygen, C βFe and C βO, are calculated from equations (3–5). Concentrations C δFe and C δO are calculated based on the properties of magnetite, namely, the activity of iron in magnetite is assumed equal to one. The activity of oxygen in magnetite is found from the condition of chemical equilibrium (13). Thus, the equation obtained from relation (20) will contain one unknown parameter, namely, the film thickness δ. Note that the obtained equation for the concentration C γO (or the activity corresponding to it) parametrically depends only on the concentrations in other regions (β and δ), which in turn depend on time. The thickness of the film cannot be expressed only through these parameters (concentration or fluxes), and the film growth will directly depend on time.

Film growth model within the Wagner approach

As was determined above, the diffusion processes of iron and oxygen in the oxide film are slow, and therefore the determining processes of its growth. Therefore, formulas (17) and (19) are a rather crude approximation.

According to Wagner’s theory, the flow of iron through the oxide layer as a function of coordinates can be represented (Zhang 2005):

JFeδ=-CFeδDFeRTdμFedx=12αDFeCFeδdlnPO2dx=23DFeCFeδdlnPO2dx , (23)

where μFe – chemical potential of iron, PO2 – partial pressure of oxygen, x – coordinate, α=⁴/₃ for Fe3O4.

The fulfillment of relation (23) gives a constant flow of iron [10]:

JFeδδ=23CFeδaO2(0)aO2(δ)DFedlnaO2 (24)

where aO2 – oxygen activity in oxide at a point with a coordinate x, determined as aO2=PO2PO20 , where P 0O2 = 1 atm., CFeδ=3MFeρoxMO [kg/m3] ‒ the average density of iron in magnetite, where MOx = 3MFe + MO, MFe and MO – molar masses of iron and oxygen, respectively, ρox ‒ magnetite density, aO20=3MFeρoxMoxCO2S0 , C SO2 (0) – oxygen saturation concentration in magnetite.

Using relation (24) and the dependence of the diffusion coefficient of iron cations, which takes into account the transport of iron ions over vacancies and along the grain boundaries of the metal in the framework of the theory of point defects (Topfer et al. 1995; Zhang 2005; Martinelli et al. 2008a, 2008b; Mikityuk 2010):

DFe=DV0exp-EV/RTPO22/3+DI0exp-EI/RTPO2-2/3 , (25)

D 0V and D 0I – constants characterizing the transport of vacancies and between grains, respectively), we can obtain the value of the iron flux from the layer δ as a function of oxygen activity at the layer boundary (x = δ):

JFeδ=CFeδδAaO2(δ)2/3-aO2(0)2/3-BaO2(δ)-2/3-aO2(0)-2/3 , (26)

A=DV0exp-EV/RTPO20-2/3 ,

B=DI0exp-EI/RTPO202/3 ,

In the accepted notation aO2 (δ) = a γO2, aO2 (0) = a δO2.

Assuming that in the coolant flow C βFe = 0 and neglecting the last term in equation (22) we obtain that:

4MO3MFeJFeδ-χFeCFeγχ0C0β-C0γ (27)

Using the relations (26, 27) and kγ = (a γFe)3(a γ0)4, we can obtain the equation for C γO (which does not have a simple analytical solution, but can be solved numerically (the only solution)).

The rate of change of the film is determined by the flow of oxygen J βO:

dδdt=Joβρox1+34MFeMo (28)

From the general formulas (22) and (28), one can derive particular approximations of Martinelli (Maruyama et al. 2004) and Zhang (Zhang and Li 2005).

Martinelli Approximation

Neglecting both the oxygen flow in the film and the iron flow in the coolant flow, we can obtain the approximation of film growth according to Martinelli (Martinelli et al. 2008a). In this case

4MO3MFeJFeδ=JOβ (29)

Then the film growth equation (28) takes the form:

CFeδdδdt=JFeδ (30)

Taking into account expression (26), we obtain an expression for the dependence of film growth on time:

δ2 = t . Kp, (31)


Zhang’s approach

Neglecting the oxygen flux in the film, we can obtain the equation of film growth in the Zhang approximation (Zhang and Li 2005):

JOβ=4MO3MFeJFeδ-JFeβ4MO3MFeJFeδ-χFeCFeγ-CFeβ (32)

Then the film growth equation takes the following form:

dδdt=Kp2δ-χFeCFeγ-CFeβCFeδ (33)

Denote the value χFe'=χFeCFeγ-CFeβCFeδ , which reflects the dissolution rate (corrosion) of the film. Then the equation for changing the film takes an even simpler form:

dδdt=Kp2δ-χFe' (34)

This equation corresponds to the Tedmon model (Martinelli et al. 2008b). If the parameters of the equation are considered as constants, then the equation has an analytical solution:

t=-δ-δ0χFe'-Kp2χFe'2ln1-2χFe'Kpδ-ln1-2χFe'Kpδ0 (35)

where δ0 – initial film thickness.

Expression (35) shows that the film thickness at t → ∞ reaches its limit δc=Kp2χFe' .

At δ0 << δc and 2χFe'2Kp1 expression (35) can be approximately written as

δ-δ0=Kpt1/2-23χFe't=Kpt1/21-232χFe'2Kpt (36)

In this case, we can assume that the parabolic film growth is quite well performed and the film corrosion can be neglected at times tKp2χFe' .

Approximation of a semi-empirical account of diffusion through an oxide film

The derivation of the film growth rate constant of the parabolic law (31) is based on the model (Topfer et al. 1995), where, when modeling the diffusion of iron cations, the assumption was made that the concentration of vacancies and interstitial ion concentrations are independent of porosity and grain size of magnetite, that is, ion diffusion within the framework of the vacancy model for ideal "magnetite" (Hallstrom et al. 2011). Such models are quite fundamental in nature and allow one to calculate the diffusion coefficient of iron ions, which is necessary for calculating the film growth rate constant in the framework of the Wagner model. The main result of the model (Topfer et al. 1995) and similar ones (Diekmann 1998) is that in the region of low partial oxygen pressures Kp~PO22/3~dOγ4/3 . Experimental studies of film growth on various alloys show slightly different dependences of the film growth constant on the partial pressure of oxygen than was obtained for pure magnetite using models of the type (Topfer 2005; Hallstrom et al. 2011).

Since in reality the film for the overwhelming majority of materials under consideration is not pure magnetite, but is at least two-layer, the models developed in (Topfer 2005; Hallstrom et al. 2011) it’s not necessarily to describe the dependence of the film growth rate constant (conditional name) on the partial oxygen pressure like Kp~PO22/3 . Since in the works (Surman 1973; Smith 1982; Saito. et al. 1985; Sato et al. 2002) the dependencies PO20.145, PO20.135, PO20.141 and PO20.279;0.313 have been experimentally found for different types of steel correspondingly. These empirical estimates imply the fulfillment of the parabolic law of film growth (31), i.e., without taking into account possible mass transfer processes leading to its dissolution (corrosion). The latter was taken into account in (33) and for large times provides a nonparabolic change in the growth of the oxide film. Such a different dependence of the film velocity constant is determined by both the composition of the metal so that by the properties of the grains, which leads to the failure of the theoretical dependence Kp~PO22/3 , characteristic to the "ideal" magnetite (Diekmann 1998). In (Sato et al. 2002), it was assumed that the film growth rate constant is generally determined by the diffusion of chromium ions in chromium oxide, the theoretical dependence of which Kp~PO23/16, which, in principle, is somewhat closer to the experimental values.

Note that in (Surman 1973; Smith 1982; Saito et al. 1988; Diekmann 1998; Sato et al. 2002) the oxidation of steel surface was carried out at low partial oxygen pressures typical for partial pressures in liquid lead, but in a gaseous medium. Therefore, the processes of film dissolution and erosion, in principle, could not play a significant role. The results of experiments in a stationary coolant (Steiner 2009, Askhadullin et al. 2011) led to an empirical dependence PO20.125 for steels EP823-(Sh) and E302-(Sh). For T91 steel, the available experimental data shows a dependence of the type PO20.11 .

The parabolic dependence of the film growth (31), in principle, is quite well confirmed experimentally (Aiello et al. 2004), if corrosion-erosion processes have low intensity (for example, in a gaseous medium), and the Wagner approach is currently the only reasonably substantiated approximation. Deviations from the parabolic law can be explained within the framework of the approach (Zhang and Li 2005; Steiner et al. 2008), in which the film growth is described by equation (34).

In the general case, the film growth constant can be written as:

Kp=Aexp-QRTPO2γn-PO2δn (37а)

P γO2, P δO2 – partial oxygen pressure at the interface between the coolant-oxide and steel-oxide media, respectively, while unknown parameters are determined empirically.

The partial pressure at the coolant-oxide interface, using relation (10), is defined as:

PO2γ=COCOS2exp2ΔGPbORT (37b)

The value of the partial pressure at the oxide-steel interface is much smaller than the value of the partial pressure of oxygen at the surface and can be neglected (Robertson and Manning 1988). In the case of lead – bismuth coolant, it can be estimated that since the binding energy of bismuth oxide is much lower than the binding energy of lead oxide, the formation of the former does not occur until a significant amount of lead oxide is formed. Therefore, in all calculations, the formation of bismuth oxide is neglected for its smallness.

When describing the film growth constant using formula (37), the ratio of fluxes at the boundary (22) takes the following form:

4MOCFeδ3MFeKp2δ-χFeCFeγ-CFeβCFeδχ0C0β-C0γ (38)

Given the discussion of relation (31), this expression means that δcδ-1~χ0C0β-C0γ, where δc=Kp2χFeCFeγ-CFeβCFeδ .

That is, upon reaching the maximum value of the film δc, the oxygen concentrations in the film and in the volume are equal. That is, the global isoconcentration regime is established.

Oxide film growth criterion

Based on the foregoing, in the previous section, we rewrite equation (38) in the following form:

4MO3MFeδcδ-1χ0C0β-C0γχFeCFeγ-CFeβ (39)

The expression for the maximum value of the film, provided that there is no iron in the coolant flow and the saturation concentration of iron is reached at its surface, takes the form:

δc=Kp2χFeCFesaFeγρ3MFeMoxρox0.2KpχFeCFesaFeγ (40)

Here, the densities of lead and magnetite are taken for temperature 550 °С, C SFe – fractional concentration of iron saturation. Thus, the film growth condition from the initial value δ0, δc=0.2KpχFeCFeSaFeγ>δ0 , which is accompanied by simultaneous natural conditions (39) that the concentration of oxygen in the coolant is greater than at the wall and vice versa for iron. Note that the parameter δcδc (T,P γO2,υ,dh) will characterize either the growth or dissolution of the film with respect to the initial thickness to an equilibrium value δc.

Parabolic law for the growth of an oxide film

From relations (22) and (38) it follows that the condition of the parabolic law of film growth corresponds to the condition

JFeδJFeβ=δcδ1 (41)

Based on the smallness parameter for equation (39), parabolic behavior is expected for times ttk,tk=Kp2χFe'2 . The relation also holds δc2=Kp2tk . Thus, if we conditionally consider the fulfillment of the parabolic law with time less than ttk / 20, which is due to the smallness parameter of equation (36), then this is expected for film thickness no more than δδc/10 . The latter meets the condition JFeδJFeβ10 .

Let us estimate the order of times at which parabolic growth of the oxide film is expected. Based on the available experimental information (Nuclear Energy Agency 2015), we determine the expected order of the film thickness in the region of several tens of microns, for example, 50 μm, and suppose that the film grew parabolic to this value. Then the maximum film thickness is about 150 microns. Further, the expected order of the mass transfer member can be estimated:

χFe'=χFeCFeγ-CFeβCFeδχFecFesafeγρMox3MFeρox10-10afeγ , [m/s] (42)

Let us take the activity of iron at the surface of the film of the order 10-1–10-2, and mass transfer constant χFe ~ 10-4[m/s] at characteristic parameters of temperatures and coolant velocity. Therefore, tk=δcχFe' ~ 300–3000 hours. The latter means that parabolic behavior is expected at a time of less than 15–150 hours. Thus, if the experimentally observed film thicknesses are relatively close to the saturation film, then a parabolic approximation to this value will take place over a relatively short period of time, up to several hundred hours.

In order for the parabolic film growth to last for several thousand hours, the saturation parameter of the film growth δc must be in the order of 10-2 m.

This corresponds to the fact that the film growth constant should be of the order of 10-11 m2/s at Т = 600 °С and an order of magnitude or two less at temperatures up to 400 °С. Such a film growth constant is very large and gives a huge value of the film thickness, which does not correspond to the known experimental data.

When the activity of iron at the surface of the film is of the order 10-3 parabolic film growth is provided for times from several thousand hours to several tens of thousands of hours depending on temperature.

Double layer film

An oxide film, as a rule, consists of two different layers and a transition region between the inner oxide layer and steel (Robertson and Manning 1988; Martinelli et al. 2008a; Alekseev et al. 2010a; Nuclear Energy Agency 2015). The inner layer is chromium-rich spinel, in which the concentration of chromium is approximately equal to its concentration in the metal, while the outer layer is a porous iron-rich oxide. A thin transition layer is observed in which the iron concentration gradually changes from the concentration in the oxide to the concentration in the alloy, the chromium concentration is practically unchanged, and the oxygen concentration gradually changes from the concentration in the oxide to zero.

Numerous experimental data show that the layers of magnetite and spinel are comparable in their thicknesses, which is explained in principle in the framework of the "accessible space" model (Robertson and Manning 1988; Martinelli et al. 2008a). The model assumes that the spinel layer Fe-Cr grows into an accessible space created by vacancies formed in steel due to diffusion of iron, while the mobility of iron ions is several orders of magnitude higher than the mobility of chromium ions. Vacancies can accumulate upto the formation of nanocavities and nanochannels, which facilitates the access of oxygen and provides the formation of spinel. But the rate of vacancy formation is completely determined by the diffusion of iron cations, and it is this process that determines the rate of both magnetite growth and the spinel formation rate. That is why the ratio of the thicknesses of magnetite and spinel Fe13Cr0.6O4 is approximately 1.16 (Robertson and Manning 1988; Martinelli et al. 2008a). In turn, the degree of diffusion of iron cations depends on the partial pressure at the interface between the media and determines the rate constant of the (parabolic) film growth Кр (40).

Based on the foregoing and following the ideas developed in (Robertson and Manning 1988; Martinelli et al. 2008a), it does not make much sense to complicate the engineering approach for calculating the changes not only in the magnetite film, but also in its spinel part, since the latter is completely determined by the diffusion (departure) of iron cations from the metal into magnetite. Let’s consider, in general terms, an equation similar to equation (28) for film growth, taking into account spinel formation and magnetite growth:

1Sνoxt1SνSpt+1Sνmagt=ρmagMmagδmagt+ρSpMSpδSpt=Joβ4Mo ,

where υox, υmag, υSp – the amount of the substance of the oxide, which is the sum of the substance of magnetite and spinel. Based on the above accepted model, spinel growth is proportional to magnetite growth, δSp = αδmag and the film growth equation can be rewritten as:

ρoxMoxδtρmagMmag1+αρSpMmagρmagMSp1+αδt=Joβ4Mo ,


1+αρSpMmagρmagMSp1+αδt=Joβ4MoMmagρmag , (43)

where ρmag,Sp, Mmag,Sp – densities and molar masses of magnetite and spinel Fe3-xCrxO4, where x reflects the stoichiometry of the compound, respectively. Note that the densities of spinel and magnetite are fairly close to each other (Zhang 2013) and since ρSpMmagρmagMSp1 with accuracy of less than 3%, then 1+αρSpMmagρmagMSp1+α~1 and equation (43) takes the form of (31) with high accuracy. Thus, within the framework of this model, a separate calculation of the growth of magnetite and spinel films does not make much sense, but only significantly complicates the calculations.

The interface between the outer and inner layers with sufficient accuracy coincides with the original metal surface (Robertson and Manning 1988; Schroer et al. 2012; Tsisar et al. 2017). Based on the mass balance and on the assumption that iron does not enter the coolant, but remains in the oxide outer and inner layers, and that all oxidized chromium is in the inner layer, we can write the ratio of the thicknesses of the outer layer and the inner (Robertson and Manning 1988):

δoutδin=x1+CFe/CCr3-1 , (44)

where CFe and CCr determine the concentration of iron and chromium in the alloy (mol /m3).

It should be noted that the steel studied in sufficient detail for the structure of the oxide film does not give an unambiguous answer to the question of the presence or absence of an external layer of magnetite. But, in general, it can be argued that under nominal and similar modes, and in the case of dynamic experiments (coolant velocity 1–3 m/s)) (Zhang and Li 2007; Weisenburger et al. 2008; Schroer et al. 2009; Tsisar et al. 2017), the magnetite layer is most likely absent due to erosion or dissolution of iron oxide films in the coolant. Although in some cases the same steels (for example, T91) show the presence of a layer of magnetite in dynamic experiments (Alekseev et al. 2010a; Schroer et al. 2012).

In static experiments, as a rule, a sufficiently bright two-layer structure of the film is observed and, as a rule, it is significantly thicker compared to dynamic experiments, in equal conditions.

Model of volume sources in a coolant flow

For an individual unit volume of coolant, the activities of oxygen, iron and magnetite in a state of chemical equilibrium are connected by the relation (the coolant index β in this section is omitted):

k=aFe3O4τ+ΔτaOτ+Δτ4aFeτ+Δτ3 . (45)

Chemical equilibrium occurs due to the loss of oxygen and iron, the accumulation of magnetite:

aFe3O4τ+Δτ=aFe3O4τ+ΔaFe3O4;aFeτ+Δτ=aFeτ-ΔaFe;aOτ+Δτ=aOτ-ΔaO. (46)

The accumulation of magnetite is associated with a decrease in oxygen and iron ratio:

ΔaFe=αΔaFe3O4CsFe3O4CsFe;ΔaO=βΔaFe3O4CsFe3O4CsO;α+β=1. (47)

Considering the stoichiometric coefficients of the reaction of magnetite formation and the ratio of the atomic weights of oxygen and iron, we can obtain the coefficients α = 0.726 and β = 0.274.

Combining relations (45–47), we can obtain an equation for determining the profit of magnetite in the process of establishing equilibrium:

aFeτ-αΔaFe3O4CsFe3O4csFe=aFe3o4τ+ΔaFe3O41/3kaOτ-βΔaFe304CsFe3O4CsO41/3. (48)

Equation (48) is solved at each time step by the tangent method. According to the found values of the magnetite profit, the values of the loss of oxygen and iron are found at each step of solving the transport equations. Thus, the volume sources of impurities in the coolant are calculated by the relations:

GFe3O4=ΔaFe3O4CsFe3O4/Δτ; (49)



For the near-wall cells of the computational grid, surface flows of oxygen and iron are added to volume sources.

Parabolic film growth constant

As a rule, the parabolic film growth constant is determined experimentally in a gaseous medium or under the condition that the processes of corrosion (dissolution) of the film and erosion can be neglected. Obviously, these two effects will be quite small in a static experiment. Under the assumption that the yield of iron from steel to the coolant is small compared to its consumption by an oxide film, the Martinelli approximation (29) should work well, which ensures parabolic film growth when solving equation (28). In this case, the oxygen flux is inversely proportional to the thickness of the oxide film. Thus, to determine the parabolic film growth constant, it is necessary to measure the oxidation rate of the surface at various temperatures and oxygen concentrations.

In the empirical determination of the parameters of the parabolic constant, the assumption of a parabolic growth of the oxide film in the static mode is fundamental. This assumption works fine provided that the diffusion yield of steel components in the coolant is small. In passivation mode in the initial period of formation of oxide films (Ivanov et al. 2005), the diffusion yield of metal components into the coolant is very intense, and its fraction can reach ~ 50% and the diffusion flux does not have to be proportional to the oxygen flux. According to the data of (Ivanov et al. 2005), for an oxidation duration of more than 400 hours, the proportion of steel components entering the coolant of the total amount of oxidized is ~ 1% and less, with an oxidation duration of 200 to 400 hours, this proportion is ~ 10%. Therefore, during long campaigns, at least more than several hundred hours, oxygen consumption occurs only on the surface and parabolic film growth in static experiments should be manifested to a greater extent than, for example, in dynamic ones. Thus, the empirical determination of the parabolic constant in the best way can be made by the thickness of the film after a sufficiently long oxidation (more than ~ 400 hours) and when using the static mode. Unfortunately, there are extremely few such experimental data.

Such experiments were carried out at JSC "SSC RF-IPPE" (Askhadullin et al. 2011), where a technique was developed that made it possible to estimate the flows of oxygen consumed for the oxidation of the steel surface under variations in temperature and oxygen conditions of the coolant. There is a matrix of various experimental values of the flow of oxygen consumed for steel oxidation, and the temperature and partial pressure of oxygen in the coolant (Askhadullin et al. 2011; Ivanov et al. 2017).

Based on the available experimental data (Askhadullin et al. 2011; Ivanov et al. 2017), the following dependence was proposed for EP-832 steel:

Kp=exp39.5pO20.125exp-30000T,[μm2/h], (50)

Similar experimental studies, by definition, were also carried out for EP-302 steel. The following dependence is proposed for a parabolic constant:

Kp=exp31pO20.125exp-20000T,[μm2/h] (51)

For comparison, we present the dependences of parabolic constants for the basic oxygen regime, obtained in (Hwang and Lim 2010), based on an analysis of the available experimental data (Fazio et al. 2001; Barbier and Rusanov 2001; Aiello et al. 2004; Balbaud-Celerier 2004; Zhang et al. 2005; Lim et al. 2010; Martinelli et al. 2008a; Hwang and Lim 2010; Nuclear Energy Agency 2015).

For steels of ferritic-martensitic class (9–12Cr) НТ9, Т91 the following dependence is obtained:

Kp=1.156×10-3exp-135693RT,cm2/s (52)

For austenitic steel AISI 316L:

Kp=1.156×10-3exp-135693RT,cm2/s (53)

In (Zhang 2013), the following dependences were also obtained for HT9 and T91 steels that take into account the partial pressure at the interface between media, namely, steel Т91:

Kp=2.096×10-5co0.27exp-157802RT,m2/s, (54)

And for steel HT9:

Kp=9.735×10-8co0.27exp-126802RT,m2/s, (55)

where с0(ppm, millionth share) ‒ mass fraction of oxygen.

An analysis of the empirical dependences of parabolic constants (Fig. 2) shows that among ferritic-martensitic steels, oxygen consumption will be the smallest for EP-823 steel, which, under similar other conditions, will provide the smallest film growth. Austenitic steels (for example, the AISI 316L shown here) generally have lower oxygen consumption. Note that in (Barbier and Rusanov 2001; Balbaud-Celerier et al. 2004) the geometric data of the circulation loops are not presented in sufficient detail, on the basis of which the temperature dependences for the parabolic constant are obtained, which introduces additional calculated errors.

Figure 2. 

Empirical dependences of parabolic constants for various types of steels (Zhang and Li 2007; Hwang and Lim 2010).

Corrosion processes. Mass transfer coefficient

It was also noted in (Zhang and Li 2007; Steiner 2009; Alekseev et al. 2010c) that the parabolic dependence does not necessarily describe the available experimental data. Further, when taking into account corrosion processes, we will adhere to the general classification of corrosion oxidation regimes adopted in (Zhang and Li 2007), namely, a comparison of Kp,c ~ Kp / δ (Kp ‒ parabolic constant) and Qc,o ‒ corrosion rates during corrosion oxidation:

Kp,с ~ 0, metal corrosion without oxide formation;

Kp,с < Qc,o, oxidation affects corrosion;

Kp,с > Qc,o, corrosion oxidation process;

– Qc,o ~ 0, oxide corrosion, slight erosion;

Kp,с = Qc,o, no new oxide formations.

The general concept of corrosion oxidation due to mass transfer effects is the assumption that the rate of corrosion oxidation is proportional to the corrosion of the metal itself (Zhang and Li 2007):

Qc,0 ~ Qc = χ (cscb) ~ χo (ν,D,d)Vn exp(q (ν,D,d)/RT), (56)

where Qc ‒ corrosion rate of the actual metal surface, V ‒ coolant speed, χo and q ‒ constants depending on the properties of the coolant and flow geometry, n – determined by type of flow (Alekseev et al. 2010c).

At low coolant speeds, the corrosion rate of a metal surface is fully or partially controlled by the mass transfer effect (Chen et al. 1992; Balbaud-Célérier and Barbier 2001). The corrosion rate increases with increasing coolant velocity. At higher coolant velocities, the reactions of diffusion transfer of metal ions become equal to the diffusion rates of metal particles in the coolant or even become lower than them. As a result, the overall corrosion rate begins to be controlled by the reaction rates in the metal (activated controlled), and the corrosion rate ceases to depend on the speed of the coolant. At even higher speeds, mechanical erosion begins.

Thus, using the approach of (Zhang and Li 2007) for a pure metal, the corrosion rate is defined as QcχFeC sFe, while Qc,oχFeCFeγMox3MFeρox , that is

Qc,oQc=Moxρmetall3MFeρoxaFeγaFeγaFeγd,T,C0 (57)

Thus, the activity of iron at the surface of the oxide film gives an understanding of how much the corrosion process on the film is slower compared to corrosion of the metal surface. As a rule, this value is about 100–1000. Based on the chemical equilibrium on the wall, the greater the oxygen activity in the coolant, the lower is the activity of iron and the less corrosion of the surface. Of course, it must be kept in mind that a large (close to saturation) oxygen activity can lead to excessive and undesirably large film growth, that is, the yield of iron can be minimized by an increase in oxygen in the coolant, but thereby its output will increase, in fact, in the oxide film, which also is not always desirable.

There are several expressions of mass transfer coefficients (Figs 2, 8) obtained on the basis of experimental data (Harriott and Hamilton 1965; Berger and Hau 1977; Silverman 1984):

kb₋h = 0.0165υ⁻0.53D0.67V0.86d⁻0.14

ks = 0.0177υ⁻0.579D0.704V0.875d⁻0.125 (58)

kh₋h = 0.0096υ⁻0.567D0.654V0.913d⁻0.087

where υ ‒ viscosity, D ‒ diffusion coefficient, V ‒ speed and d ‒ hydraulic diameter.

Using formulas (58), numerical values are obtained for the mass transfer coefficients, but these differences do not exceed 20% at a speed of 2 m/s and less at lower speeds. In further calculations, we use the expression kb-h. Fig. 3 shows typical values of the rates of dissolution (corrosion) of the film taking into account this mass transfer coefficient (58) depending on temperature.

Figure 3. 

Typical dissolution (corrosion) rates of a film (see equation (34)).

When using experimental data to verify the calculation models, it is not always possible to extract data on the geometry of the working sections with sufficient accuracy, but from general considerations, the limits of the hydraulic diameters of these sections can be estimated. As a rule, hydraulic diameters are within the range of 10–100 mm. Corrosion rate parameter change χFe (34) when changing the hydraulic diameter from 10 to 100 mm lies within 20%, as can be seen in Fig. 4.

Figure 4. 

The value of the parameter of dissolution (corrosion) of the film (see equation (34)) depending on the hydraulic diameter.

Justification of the calculation method

The main goal of a computer program is to calculate the thickness of an oxide film and the associated concentrations of oxygen and iron. Thermohydraulic calculations are a necessary basis, but the accuracy of the calculated thermohydraulic characteristics (in our case, this is the field of temperatures and velocities) does not have to exceed the accuracy of the calculation of physicochemical processes (in our case, this is the thickness of the oxide film and the concentrations accompanying it). The main errors in the calculation are introduced by the uncertainty of empirical parameters.

The film growth in the static case is parabolic in nature. It is static experiments that serve as the basis for determining the parameters of the parabolic constant (37a), which is determined by the oxygen consumption of the film. Based on the fact that in the static case there is a parabolic dependence of the film thickness growth δ2 ~ Kpt, then the relation between the relative errors of the film thickness and the uncertainty of the parabolic constant is written as εKp = 2εδ.

The experiments on oxygen consumption (Askhadullin et al. 2011; Ivanov et al. 2017; Salaev et al. 2018) make it possible to relatively accurately determine only the parameters of the dependence on temperature and partial pressure (37b). Therefore, the absolute values of the parabolic constant can be determined only by reference experiments. In our case, these are experiments (Barbier and Rusanov 2001; Balbaud-Celerier 2004), since the experiments were performed for different types of steels on the same installation. Based on the data in Figs 810 (Balbaud-Celerier and Terlain 2004), we accept εKp ≈ 0.6. Since the accepted empirical dependence (37a) is determined by three parameters, A, Q, and n, the relative uncertainty of the parabolic constant is determined by the relative errors for these parameters. The dependence on the uncertainty of the degree of partial pressure can be neglected due to its smallness, then the uncertainty of the parameters A and Q, εA, εQ contribute to the uncertainty of the parabolic constant as

εKp=εA2+Q/Tεq2 (60)

Since, based on the available data, it is impossible to unambiguously indicate which of the uncertainties makes the largest contribution, we accept this contribution to be the same, i.e.

εA=|Q/T|εQ=0.6/20.43 (61)

That is, uncertainty εA may be high enough while uncertainty εQ should be small since |Q/T| >> 1.

Thus, during variational calculations of the film thickness and impurity concentrations, the boundaries of the varied parameters are estimated А and Q.

The errors of partial pressure and average flow rate, temperatures are taken from the available experimental data.

In dynamic experiments, the flow rate and hydraulic diameter were additionally varied. The speed varied based on the data of the error of the flow meter, this value, as a rule, lies within ±10%. The hydraulic diameter ranged from 10–100 mm (see previous section).

Based on the selected intervals, the average deviation of the calculated value from the experimental ones was calculated to get the relative error

ε=1Ni=1Nabsxicalc-xiexpxiexp (62)

where xicalc, xiexp calculated and experimental values of the simulated value, respectively.

Subsequent numerical calculations showed that the errors in the calculation of dissolved oxygen concentrations are not more than 15% when comparing the results obtained using the MASKA-LM and STAR-CCM + codes, and the oxide film thickness is not more than 30% for the static case (stationary coolant) and not more than 55% for the dynamic case (coolant with average speeds of 0.1–3 m/s) when comparing the calculated and experimental data.

When considering experimental errors for oxygen concentration, it should be noted that at the moment in the literature there is no described well-established measurement procedure with a proven error in a wide range of oxygen concentrations and it, in principle, can reach significant values. When considering experimental errors for oxide film thicknesses, we note that experimental results are very often presented with no errors at all, and those data presented with errors should be considered as a spread of experimental data without evaluating the experimental errors themselves (see, for example, Fig. 11).

Simulation of film growth in the framework of the simplified 0D-model

The joint solution of equations (34) and (38) within the framework of the point model, that is, without taking into account the exact geometry (Tadmon’s approach), allows a reasonable accuracy to estimate the film growth rate depending on temperature and coolant velocity. For this, we used the film growth parameters selected for EP-823 steel, since for this grade of steel there is a fairly large amount of experimental data.

The growth of an oxide film as a function of time is considered in detail in a subsequent chapter. The most indicative characteristics parameters of film growth are its maximum value (see equations (34–36)) and the characteristic film growth time according to the parabolic law ttk,tk=Kp2χFe'2 . Based on (40), the saturation yield corresponds to the equality of the fluxes of iron through the film and in the coolant, J δFe = J βFe, which corresponds to a zero flow of oxygen. This, of course, is performed provided that a constant concentration is maintained at the selected oxygen regime.

Figs 57 show the typical model dependences for the film saturation thickness on the oxygen velocity and concentration, as well as the characteristic time of transition to parabolic growth. The simulation parameters are taken for lead coolant and steel EP-823.

The growth of the oxide film upon saturation, depending on the coolant speed, is shown in Fig. 5. Obviously, at low speeds, a substantial, almost exponential, increase in the thickness of the saturation film is observed. That is, the saturation thickness substantially depends on the velocity of the coolant.

Also, calculations showed that at temperatures of 450–650 °C, parabolic film growth is expected for times from less than an hour to several hundred hours (Fig. 6).

The calculations shown in Fig. 7 demonstrate that under the basic oxygen regime (oxygen concentration in the region of 10-8 kg/kg) and at standard coolant velocities (1–2 m/s), a very small film growth is expected (up to 1 μm). If you deviate, in the direction of increasing the oxygen content, from the basic mode, you can expect film growth to several tens of micrometers. At a low oxygen content (dash dotted line in Fig. 7), the formation of an oxide film is practically unrealistic, which can lead to corrosion.

Figs 810 show the calculations for steels EP-823, T91 and 316L for which both static and dynamic experiments were performed (Barbier and Rusanov 2001; Balbaud-Celerier and Terlain 2004), which allows a consistent analysis of these data, that is, the same parameterization of the parabolic constant is used for both types of tests. The coolant speed in the dynamic test is 1.9 m/s. Below is an analysis for T = 600 °C. Taking into account the numerical errors from the variation of the parameters, we obtain fairly good agreement with the experimental data. An interesting feature of the above calculations and experimental data is that although the steels EP-823 and T91 are quite close to each other in composition, the growth and thickness of the film for them is significantly different, which determines the adequacy of our engineering approach.

The authors of (Zhang and Li 2007; Weisenburger et al. 2011) analyzed the obtained experimental information for T91 steel. The experiments were carried out on the same installation of JSC "SSC RF-IPPE" as in (Balbaud-Celerier and Terlain 2004; Zhang et al. 2005), but with an average coolant speed of 1 m/s and under the basic oxygen regime. As discussed earlier, it is not possible to reconstruct the exact geometric parameters. Therefore, a certain spread in the values of the hydraulic diameter in the range of 0.005–0.1 m, which gives no more than 20% error in the value of the film thickness, is accepted. As an example, a calculation with an average coolant velocity of 0.3 m/s is presented to demonstrate sensitivity to this parameter. This speed was chosen intentionally to place most of the experimental data between the two calculated curves (Fig. 11). The calculated estimate demonstrates that having a spread of experimental data can be, inter alia, related to a "spread" in speed from 0.3 m/s to 1 m/s.

Figure 5. 

Dependence of the saturation thickness of the oxide film on the coolant velocity at a basic oxygen concentration of 10-8 kg/kg.

Figure 6. 

Dependence of the transition time to parabolic film growth on the coolant velocity at a basic oxygen concentration 10-8 kg/kg.

Figure 7. 

Dependence of the saturation thickness of the oxide film on the coolant velocity for various oxygen concentrations at T=550°C.

Figure 8. 

Calculated results and experimental data for oxide film growth based on static and dynamic experiments for steel EP-823 (Balbaud-Celerier et al. 2004). The calculated error is 55% for a dynamic experiment and 30% for a static.

Figure 9. 

Calculated results and experimental data based on static and dynamic experiments for T91 steel (Balbaud-Celerier et al. 2004). The calculated error is 55% for a dynamic experiment and 30% for a static.

Figure 10. 

Calculated results and experimental data based on static and dynamic experiments for 316L steel (Balbaud-Celerier et al. 2004). The calculated error is 55% for a dynamic experiment and 30% for a static.

Figure 11. 

Calculated results and experimental data for steel T91 (Weisenburger et al. 2011). The estimated error is 55%.

Numerical analysis of oxide film growth on surfaces of simple geometries using MASKA-LM and STAR-CCM+ codes

Distribution of oxygen at a given drain on the pipe wall

This section provides calculations of model oxygen distributions at various temperatures. Test calculations of the elementary model of HLMC flow in a pipe were carried out using the commercial StarCCM+ code and MASKA-LM code of own design.

The STAR CCM+ code interface allows you to use user-defined functions to model certain physical and chemical processes. The system of equations (33), (38), functions dependent on them, and semi-empirical dependencies of type (2, 20–25) were incorporated into the code for joint solution with the STAR CCM+ thermohydraulic models.

For numerical analysis, a simple model was chosen for the flow of lead coolant of a given temperature in the pipe: the length of the computational domain is 1 m; diameter 50 mm; design grid 1 million cells; Reynolds number is 105; input mass concentration of oxygen 2∙10-8 kg/kg. To simulate oxygen corrosion, the following value was selected for the drain (consumption) of oxygen to the wall 4∙10-8 kg/m2∙s. The selected amount of oxygen consumption is very significant compared to that usually found in real calculations and was chosen only to demonstrate the distribution of oxygen.

When using the MASKA-LM code, a rectangular computational grid with constant spatial steps in the Cartesian coordinate system, consisting of 275,000 cells (50 × 50 × 110) (along the directions of the coordinate axes) was entered into the computational domain. That is, the number of cells and the quality of the grid is noticeably inferior to the StarCCM+ code.

The results of calculations using the StarCCM+ and MASKA-LM codes are shown in Figs 12, 13 (for coolant temperature of 400 °C and 500 °C, respectively, the differences are minimal). Calculations are given for various turbulence models (StarCCM+) and for various turbulent diffusion coefficients used in the MASKA-LM code. Changing the coefficient of turbulent diffusion by an order of magnitude gives the maximum error in calculating the oxygen concentration of less than 15%. It is worth noting that at the moment, the methodology for determining the measurement error of the oxygen activity sensor has not been developed, and according to qualitative estimates it can be at least 50%. In this case, the choice of a turbulence model or a simplified description of turbulent motion will not affect the description of available experimental results at all.

Figure 12. 

Distribution of oxygen dissolved in HLMC in the cross section of the pipe at Т = 400 °С.

Figure 13. 

Distribution of oxygen dissolved in HLMC in the cross section of the pipe at Т = 500 °С.

Thus, a turbulent diffusion coefficient of the order of (1–3)*10-5 gives a completely adequate description of the oxygen distribution with an error of no more than 15% at the surface. With smaller oxygen sinks, the scatter of values is even smaller. The scatter in the values of the coefficient of turbulent diffusion gives a very small error in the thickness of the oxide film.

Numerical analysis of oxide film growth

This section presents model calculations of the growth of an oxide film in an EP-823 steel pipe and a beam of rods at various flow rates of HLMC and temperatures.

Figs 1416 show the calculations for different temperatures, oxygen concentrations and average speeds in the pipe using the MASKA-LM and STAR CCM+ codes.

Fig. 14 shows film growth versus time. All models, in accordance with Zhang’s approach, demonstrate the approximation of the film thickness to the maximum possible value. The calculation shown by the black solid line is made taking into account only the parabolic constant (there is no dependence on the coolant speed through taking into account mass transfer). As can be seen, for calculations at a speed of 1.8 m/s, the film growth rate close to the parabolic behavior is observed only from about 1 to 50 hours, while at a speed of 0.1 m/s it is observed up to several thousand hours. Although it is more correct to say that film growth is closer to the paralinear law in these areas. The Tedmon’s approach fairly well describes only areas close to saturation in film growth.

A similar behavior of the increase in the thickness of the oxide film is also observed at a temperature of T = 500 °C. The main difference from the previous case is that the approximation to the asymptotic value of the film thickness occurs much later (Fig. 15). All models, including the simple 0D model, show very similar results. Iron concentration distributions also show fairly similar behavior, and the main difference, again, is observed near the surface.

Fig. 16 shows the calculations for an increased oxygen concentration; in general, the concentration of the basic oxygen regime is more than an order of magnitude higher. In this case, for a rather long time, a pronounced parabolic growth of the film is observed. Comparison of calculations using the codes in question was carried out only for temperatures Т = 500 °С and Т = 650 °С. The maximum difference in the thickness of the oxide film is not more than 20%.

For the following geometry, a bunch of rods was chosen with a diameter of 10 mm, a length of 0.97 m, in a triangular package, with a step of 14 mm; the simulated area is cut off along the planes of symmetry, and the grid contains 3.1 × 106 cells when using the STAR CCM+ code; mass oxygen content of 10-8 kg/kg; temperature 500 °С, 650 °С; coolant velocity 1.7 m/s (Re~105), turbulence model K-ε when using STAR CCM+ code.

Fig. 17 shows the calculation results for the thickness of the oxide film at various temperatures. As in the previous case, the difference in the obtained calculations does not exceed 20%.

Note that in the above examples only isothermal calculations were considered.

Figure 14. 

Results of calculations of film growth carried out with MASKA-LM, STAR-CCM+ and the solution of the Tedmon’s approach at a temperature of T = 650oC. Mass oxygen concentration - 10-8 kg/kg.

Figure 15. 

Results of calculations of film growth carried out in the framework of MASKA-LM, STAR-CCM+ and the analytical solution of the Tedmon’s approach at a temperature T = 500oC. Mass oxygen content – 10-8 kg/kg.

Figure 16. 

Results of film growth calculations performed within the framework of MASKA-LM, STAR-CCM+ for various temperatures and mass oxygen content 4*10-7 kg/kg.

Figure 17. 

The results of calculations of film growth carried out in the framework of MASKA-LM, STAR-CCM+ at temperatures T=650oC and Т=500оС. Mass oxygen content is – 10-8 kg/kg. The curves are trimmed after the thickness of the films reaches the saturation value.


An engineering model is presented for a self-consistent calculation of the growth of an oxide film in circulation loops with a heavy liquid metal coolant and concentrations of impurities (oxygen, iron, magnetite). The modeling of thermohydraulic and physicochemical processes is based on solving the associated three-dimensional equations of hydrodynamics, heat transfer, convective-diffusive transport, and the formation of chemically interacting impurity components in the coolant volume and on the surface of steels.

The change in the thickness of the oxide film on the metal surface and the corresponding exit of steel components into the coolant substantially depends on the steel grade. For a more adequate justification of the evolution of the oxide film, a semi-empirical model is proposed for using the empirical parameterization of the parabolic constant not only in the equation for changing the thickness of the oxide film, but also in the mass balance equation associated with it. The parabolic constant, which is determined by the degree of oxygen consumption by a steel, obviously significantly depends on the type of this steel and on the method of its preparation. Therefore, the direct application of the Wagner approach is unlikely to adequately describe such differences in steels, while the experimental parametric dependences on oxygen consumption by steel depending on temperature and oxygen partial pressure are obviously unique for each steel.

For an adequate experimental determination of the parametric dependences for a parabolic constant, at least two conditions are necessary: the experiment should be carried out for a static case or the average coolant velocity should be low in order to exclude the influence of corrosion-erosion processes; steel should have a sufficient initial oxide film to exclude the influence of the initial large flow of iron, and possibly corrosion of steel, since with a large yield of iron the film growth does not have to obey the parabolic law even in the static case.

The developed approach regarding accounting for the main physical and chemical processes is used not only in our MASKA-LM code, but is also implemented as user-defined functions in the STAR CCM+ code. The cross-verification calculations of the two test models showed a fairly good agreement between the results in terms of describing the increase in the thickness of the oxide film.

In the future, the developed approach will be applied to more complex systems in terms of geometry, expanded for non-isothermal loops, and applied to simultaneously account for the formation of deposits in loops.


  • Aiello A, Azzati M, Benamati G, Gessi A, Long B, Scaddozzo G (2004) Corrosion behaviour of stainless steels in flowing LBE at low and high oxygen concentration. Journal of Nuclear Materials 335(2): 169–173.
  • Alekseev VV, Kozlov FA, Kumaev VYa et al. (2010a) Modeling of mass transfer and corrosion of steels in nuclear power plants with lead coolant. Obninsk, preprint "FEI-3179".
  • Alekseev VV, Orlova EA, Kozlov FA et al. (2010b) Modeling of mass transfer and corrosion of steels in nuclear power plants with lead coolant. Preprint SSC RF-IPPE No. 3128. Obninsk, 2008.
  • Alekseev VV, Orlova EA, Kozlov FA, Torbenkova IYu, Kondratyev AS (2010c) Calculation and theoretical analysis of the process of steel oxidation in a lead coolant. Issues of Atomic Science and Technology, Series: Nuclear Constants 1–2: 56–56.
  • Askhadullin RSh, Ivanov KD, Shelemetyev VM, Sadovnichiy RP (2011) Evaluation of the intensity of oxidation processes of structural steels of the first circuit of a nuclear power plant with heavy coolants. Izvestiya vuzov – Yadernaya energetika 4: 121–146. [in Russian]
  • Balbaud-Celerier F, Terlain A (2004) Influence of the Pb–Bi hydrodynamics on the corrosion of T91 martensitic steel and pure iron. Journal of Nuclear Materials 335(2): 204–209.
  • Chen TY, Moccari AA, Macdonald DD (1992) Development of controlled hydrodynamic techniques for corrosion testing. Corrosion 48(3): 239–255.
  • Diekmann R (1998) Point defects and transport in non-stoichiometric oxides: Solved and unsolved problems. Journal of physics and chemistry of solids. Journal of Physics and Chemistry of Solids 59(4): 507–525.
  • FBU "STC NRS" (2012) Design ratios and methods for calculating the hydrodynamic and thermal characteristics of elements and equipment of nuclear power plants with a liquid metal coolant. Safety Guide for the Use of Atomic Energy RB-075-12.
  • Hwang IS, Lim J (2010) Structural Developments for Lead-Bismuth Cooled Fast Reactors, PEACER and PASCAR. In: The 25th KAIF/KNS Annual Conference, World Trade Center, COEX, Seoul, South Korea, April 14–16, 2010, 10 pp.
  • Ivanov KD, Lavrova OV, Salaev SV (2005) The use of the developed technique for estimation of metal components diffusion yield from steels to study corrosion resistance of these steels in heavy coolants. Book of Abstracts of the Conference "Thermo-hydraulic safety aspects of NPP’s with fast reactors". SSC RF-IPPE Publ., Obninsk, 117–117. [in Russian]
  • Ivanov KD, Niyazov SA, Lavrova OV, Salaev SV, Askhadullin RSh (2017) Development of a method for determining the oxidation rate of structural steels in heavy liquid metal coolants. Izvestiya vuzov – Yadernaya energetika 4: 127–137.
  • Kumaev VYa, Lebezov AA, Pyshin IV, Alekseev VV (2002) MASKA-LM – code for calculating the mass transfer of impurities in liquid metal circuits. In: Heat and mass transfer and properties of liquid metals: materials conference – Obninsk, 2002.Vol. 1, 295–298.
  • Kumaev VYa, Lebezov A, Alexeev V (2005) Development and application of MASKA-LM code for calculation of thermal hydraulics and mass transfer of lead cooled fast reactors. In: The 11th International Topical Meeting on Nuclear Reactor Thermal-Hydraulics (NURETH-11) (Avignon, France, October 2-6, 2005). Paper 191: 1–6.
  • Laptev AG (2007) Boundary layer models and calculation of heat and mass transfer processes. Kazan University Press, Kazan, 500 pp. [in Russian]
  • Lim J, Nam H-O, Hwang I-S (2010) Evaluation of the corrosion resistance of FeCrAl alloy in stagnant LBE with active oxygen control. In: Actinide and fission product partitioning and transmutation – 10th Information Exchange Meeting, Mito, Japan (6–10 October 2008), 321–328.
  • Marino A, Lim J, Keijers S, Deconinck J, Aerts A (2018) Numerical modeling of oxygen mass transfer in a wire wrapped fuel assembly underflowing lead bismuth eutectic. Journal of Nuclear Materials 506: 53–62.
  • Martinelli L, Balbaud-Célérier F, Terlain A, Bosonnet S, Picard G, Santarini G (2008a) Oxidation mechanism of Fe–9Cr–1Mo steel by liquid Pb–Bi eutectic alloy at 470 C (Part II). Corrosion Science 50(9): 2537–2548.
  • Martinelli L, Balbaud-Célérier F, Picard G, Santarini G (2008b) Oxidation mechanism of a Fe-9Cr-1Mo steel by liquid Pb–Bi eutectic alloy (Part III). Corrosion Science 50(9): 2549–2559.
  • Nuclear Energy Agency (2015) Handbook on Lead-bismuth Eutectic Alloy and Lead Properties, Materials Compatibility, Thermal-hydraulics and Technologies, 2015 Edition.
  • Reynolds AD (1979) Turbulent flows in engineering applications. Energiya, Moscow, 408 pp.
  • Salaev SV, Askhadullin RSh, Legkih AYu, Ivanov KD, Lavrova OV, Niyazov S-AS (2018) Experimental test of exploratory procedure of study of the oxidation structural steels kinetics and main results. VANT: Nuclear-reactor constants 4: 104–111. [in Russian]
  • Sato I, Takaki M, Arima T, Furuya H, Idemitsu K, Inagaki Y, Momoda M, Namekawa T (2002) Oxidation behavior of modified SUS316 (PNC316) stainless steel under low oxygen partial pressure. Journal of nuclear materials 304(1): 21–28.
  • Schroer C, Konys J (2009) Quantification of the Long-Term Performance of Steels T91 and 316L in Oxygen-Containing Flowing Lead-Bismuth Eutectic at 550 °C . In: Proceedings of the 17th International Conference on Nuclear Engineering. Volume 1: Plant Operations, Maintenance, Engineering, Modifications and Life Cycle; Component Reliability and Materials Issues; Next Generation Systems [Brussels, Belgium. July 12–16, 2009], 671–679.
  • Schroer C, Wedemeyer O, Skrypnik A, Novotny J, Konys J (2012) Corrosion kinetics of Steel T91 in flowing oxygen-containing lead–bismuth eutectic at 450 °C. Journal of Nuclear Materials 431(1–3): 105–112.
  • Tsisar V, Schroer C, Wedemeyer O, Skrypnik A, Konys J (2017) Characterization of corrosion phenomena and kinetics on T91 ferritic/martensitic steel exposed at 450 and 550 °C to flowing Pb-Bi eutectic with 10−7 mass% dissolved oxygen. Journal of Nuclear Materials 494: 422–438.
  • Weisenburger A, Heinzel A, Mueller G, Muscher H, Rousanov A (2008) T91 cladding tubes with and without modified FeCrAlY coatings exposed in LBE at different flow, stress and temperature conditions. Journal of Nuclear Materials 376(3): 274–281.
  • Weisenburger A, Schroer C, Jianu A, Heinzel A, Konys J, Steiner H, Müller G, Fazio C, Gessi A, Babayan S, Kobzova A, Martinelli L, Ginestar K, Balbaud-Célerier F, Martín-Muñoz FJ, Soler Crespo L (2011) Long-term corrosion on T91 and AISI1 316L steel in flowing lead alloy and corrosion protection barrier development: Experiments and models. Journal of Nuclear Materials. 415(3): 260–269.
  • Zhang J (2013) Long-Term Behaviors of Oxide Layer in Liquid Lead–Bismuth Eutectic (LBE), Part I: Model Development and Validation. Materials Science. Oxidation of Metals 80(6): 669–685.


Other Necessary Closing Relations

As the closing relations for the system of equations of motion and equations for the transport of impurities, we used dependences to estimate the values of the coefficients of turbulent thermal conductivity, turbulent viscosity, and turbulent diffusion (Reynolds 1979).

To calculate the coefficient of turbulent viscosity, we use the relations of the three-layer Karman model (Laptev 2007; FBU "STC NRS" 2012):

vTv=0 for y+<5vTv=y+5-1 for 5y+30vTv=y+2,5-1 for y+>30

where νТ – turbulent viscosity,

y+=u*yν – dimensionless coordinate in the boundary layer, ν – kinematic viscosity,

u*=τρ – dynamic speed, τ – wall shear stress, ρ – density. For pipes with a circular cross-section and the average shear stress along the perimeter of the bar in the assembly of smooth bars, the shear stress is calculated by the formula τ=ξρu-28 (FBU "STC NRS" 2012), i.e. u*=u¯ξ8 , where ξ – hydraulic resistance coefficient.

For pipes with circular cross-section at Re = 4 × 103 ₋ 108:

ξ = 1/(1.821∙lgRe₋1.64)2

For assembling smooth rods in triangular grid with Re = 6 × 103 ₋ 108:

ξ = ξ0 (1+(h/d₋1)0.32), ξ0 = 0.21/Re0.25, where h – grid period, d – rod diameter.

In (FBU "STC NRS" 2012) other correlations are also presented for calculating the hydraulic resistance of a rod assembly, which can be parametrically used in the code.

According to the Prandtl model, the average turbulent viscosity at y ≥ δ (that is, in the flow core) is calculated by the formula (Laptev 2007):

νT = 0.4u*δ , where the thickness of the boundary layer is calculated in the formulation of the three-layer Karman model (derivation of the relation in [51]):


Thus, the turbulent viscosity in the flow is defined as:


The formulas used to determine the turbulent viscosity are not the only possible ones for use in calculations. In (Laptev 2007) alternative formulas are given for calculations both in the framework of the three-layer and in the framework of the two-layer model. In (Laptev 2007), a sufficient similarity of the results of using different models is demonstrated.

Estimates of the coefficient of turbulent diffusion of an admixture are based on the analogy of transport processes, the turbulent Schmidt number ScT = νT / DT ≈ 1.