Numerical simulation of oxidation processes in a
cross–flow around tube bundles

arXiv:1708.04138v1 [cs.CE] 10 Aug 2017

Alexander G. Churbanovb , Oleg Ilieva , Valery F. Strizhovb ,
Petr N. Vabishchevichb,c,∗
a Fraunhofer

Institute for Industrial Mathematics ITWM, Fraunhofer-Platz 1, D-67663
Kaiserslautern, Germany
b Nuclear Safety Institute, Russian Academy of Sciences, 52, B. Tulskaya, Moscow, Russia
c North–Eastern Federal University, 58, Belinskogo, Yakutsk, Russia

Abstract
An oxidation process is simulated for a bundle of metal tubes in a cross–flow.
A fluid flow is governed by the incompressible Navier–Stokes equations. To
describe the transport of oxygen, the corresponding convection–diffusion equation is applied. The key point of the model is related to the description of
oxidation processes taking into account the growth of a thin oxide film in the
quasi–stationary approximation. Mathematical modeling of oxidant transport
in a tube bundle is carried out in the 2D approximation. The numerical algorithm employed in the work is based on the finite–element discretization in
space and the fully implicit discretization in time. The tube rows of a bundle
can be either in–line or staggered in the direction of the fluid flow velocity. The
growth of the oxide film on tube walls is predicted for various bundle structures
using the developed oxidation model.
Keywords: Cross–flow around tube bundles, oxidation of metals, parabolic
kinetics, Navier–Stokes equations, convection–diffusion equation,
finite–element discretization, finite–difference schemes
2010 MSC: 76D05, 76R50, 65M32

1. Introduction
Many industrial applications involve fluid flows through tube structures. An
example is the heat exchange equipment for nuclear power plants. A study of
heat and mass transfer in these systems is of great practical interest [1, 2]. In
addition to large–scale experimental studies [3, 4, 5], computational technologies
∗ Corresponding

author
Email addresses: achur@ibrae.ac.ru (Alexander G. Churbanov),
iliev@itwm.fraunhofer.de (Oleg Iliev), vfs@ibrae.ac.ru (Valery F. Strizhov),
vabishchevich@gmail.com (Petr N. Vabishchevich)

Preprint submitted to arXiv.org

March 14, 2018

are in common use for theoretical considerations of key aspects of heat and mass
transfer processes in tube bundles [6, 7, 8, 9].
In numerical simulation of cross–flows around a bundle of tubes, the emphasis is on predicting hydrodynamic and thermal phenomena. Using finite–
difference, finite–volume or finite–element methods, laminar or turbulent flows
are studied at various technological conditions (see, for example, [10, 11, 12, 13,
14, 15, 16, 17, 18, 19] and the bibliography cited in these works).
To investigate numerically the oxidation process in tube bundles, two additional phenomena should be taken into account in the mathematical model.
The first of these is associated with a distribution of oxidant in the intertubular space. The second phenomenon describes the oxidation of tubes, namely,
the formation and growth of an oxide film on tubes in the cross–flow. As a
rule, oxidant concentrations are very small and do not affect fluid dynamics.
Under these conditions, the oxidant transport can be described using the conventional convection–diffusion model [20]. Numerical simulation of flows around
tube bundles in the presence of mass transfer is performed [21] similarly to the
calculation of the temperature. Various oxidation models were developed [22]
to describe mass transfer in these flows.
In the present work, for predicting the oxidation process in a cross–flow
around tube bundles, a new model is developed to describe the growth of oxide
films. It includes both the linear kinetics and classical parabolic kinetics of
oxidation. It is based (see, e.g., [23, 24]) on the quasi–stationary approximation
for a thin oxide film. Using this model, a boundary value problem of mass
transfer is formulated. The main feature of this problem is a nonlinear boundary
condition of Robin–type for the oxidant concentration on tube surfaces. The
local thickness of the film is considered as the desired value, which, in particular,
is explicitly included in the boundary condition for the oxidant concentration.
The paper is organized as follows. A new mathematical model of oxidation
is developed in Section 2. The fluid flow in the intertubular space is governed
by the stationary incompressible Navier–Stokes equations. The flow around circular tubes arranged in in–line or staggered bundles is considered in the 2D
formulation. In Section 3, the computational algorithm for predicting hydrodynamics is given. The numerical procedure employed in the work is based on
using triangular grids and finite–element dscretization in space. The numerical results of oxidation in the cross–flow around tube bundles are presented
in Section 4. The oxidant transport is described by the unsteady convection–
diffusion equation. To model oxidation processes, the combined model of linear
and parabolic kinetics is applied taking into account the growth of the oxide
film. The results of the work are summarized in Section 5.
2. Mathematical model
We consider a cross–flow around a bundle of circular tubes in the 2D formulation. The sketch of the problem is given in Fig. 1, where x1 is the vertical
(longitudinal) coordinate and x2 is the horizontal (transverse) coordinate. The
tube rows of a bundle can be either in–line or staggered in the direction of the
2

fluid velocity. They form a periodic structure of cylinders of the same circular
cross-section. The fluid flow contains oxygen, which oxidizes metal surfaces of
the tubes. Taking into account the symmetry of the flow, it is possible to work
only with rectangular subdomains shown in Fig. 1 by solid lines for both configurations. The part Ωf of these rectangles is occupied by the fluid, whereas the
part Ωs corresponds to the tubes. The solid tube boundaries are denoted by Γs
and Γsym stands for the symmetry boundaries. The fluid with oxygen enters
into the subdomains through the boundary Γin and outflows on the boundary
Γout . Mass transfer processes are considered in the computational domain Ωf .
Γin
Γsym

Γs

x1
Ωs
Ωf
x2

Γout

Figure 1: Cross–flow around tube bundles: in–line (left) and staggered (right) arrangement
of tubes

2.1. Hydrodynamic processes
To describe a cross-flow in a tube bundle, the incompressible Navier–Stokes
equations are used:


∂u
%
+ u · ∇u + ∇p − µ∇2 u = 0,
(1)
∂t
∇ · u = 0,

x ∈ Ωf ,

t > 0,

(2)

where u(x, t) and p(x, t) are the velocity and pressure of the fluid, respectively,
while µ > 0 and ρ > 0 are, respectively, the viscosity and density assumed to
be constant.

3

Suitable boundary conditions are imposed on ∂Ωf . These are the normal
and tangential velocity components and/or components of the force applied to
the boundary, which are determined by the expression σn. Here
σ = µ(∇u + (∇u)T )
is the viscous stress tensor and n stands for the outer unit normal to the boundary. The uniform velocity profile with the value ū is specified at the inlet:
u · n = −ū,

u × n = 0,

x ∈ Γin .

(3)

The pressure p̄ and condition for the absence of tangential forces are prescribed
at the outlet:
p − σn · n = p̄, σn × n = 0, x ∈ Γout .
(4)
The no–slip and no–permeability conditions are imposed on the rigid walls of
tubes:
u · n = 0, u × n = 0, x ∈ Γs .
(5)
On the symmetry boundaries, we put the slip conditions:
u · n = 0,

σn × n = 0,

x ∈ Γsym .

(6)

2.2. Oxygen transport
We denote the oxygen concentration in the fluid by c(x, t), which is measured
in particle number per unit volume. Assuming no fluid–phase reactions, the
spatio–temporal evolution of the oxygen concentration is given by
∂c
+ ∇(uc) − D∇2 c = 0,
∂t

x ∈ Ωf ,

t > 0,

(7)

where D > 0 is the solute diffusion coefficient, which is assumed to be scalar
and constant.
A known concentration of oxygen is specified at the inlet:
c = c̄,

x ∈ Γin ,

(8)

where c̄ > 0 is constant. We prescribe zero flux of the solute at the outlet and
symmetry boundaries as follows:
D∇c · n = 0,

x ∈ Γsym ∪ Γout .

(9)

The key point in the description of the oxygen transport is related to modeling oxidation of metal surfaces. In our case of small concentrations, it is natural
to consider a thin layer of oxide. To describe the growth of the oxide layer,
various assumptions are employed [22].
Following [23, 24], we distinguish three stages of the oxygen transport:
1. oxygen transport to the external surface of the oxide film;
2. oxygen diffusion through the oxide film towards the oxide–metal interface;
4

d
c
cs

fluid

oxide

metal

Figure 2: Sketch of the oxide film growth

3. chemical reaction of oxidant with the reactive element at the oxide–metal
interface.
The second stage of the oxygen transport is associated with parabolic kinetics
(the parabolic dependence of the oxide film thickness on time), whereas the
third one is linear kinetics (the linear time–evolution).
Let us denote the oxygen concentration at the fluid–oxide boundary as c,
and at the oxide–metal interface as cs (see Fig. 2). Here d is the oxide film
thickness.
Taking into account the fact that the thickness of the oxide film is small,
we can treat the film growth as a quasi–stationary process. In this case, the
oxidant flux through the film is continuous. In addition, this flux is continuous
both at the fluid–oxide boundary and at the oxide–metal interface due to the
mass conservation law.
The mass flux at the fluid–oxide boundary from the liquid side is
qf = −D∇c · n,

(10)

where n is the outer unit normal for the fluid domain.
For the mass flux in the oxide film, we have
qo = D0

c − cs
,
d

(11)

where D0 is the diffusion coefficient of oxygen in the oxide film.
Oxygen arrived at the oxide–metal interface participates in the chemical
reaction of oxidation. The mass of oxygen involved in oxidation is evaluated as
qm = kcs ,

(12)

where k is the oxidation rate (the first–order chemical reaction).
For the quasi–stationary regime, we have
qf = qo ,

qo = qm .

5

(13)

From the second equality in (13), in view of (11), (12), we can eliminate cs .
Then, the first equality in (13), in view of (10), (11), makes possible to formulate
the boundary condition for the concentration of oxygen at the fluid boundary:
D∇c · n +

kD0
c = 0,
D0 + kd

x ∈ Γs ,

(14)

which depends on the thickness of the oxide film.
The time–variation of the thickness of the oxide film at the boundary point
is evaluated via the reacted oxidant:
kD0
∂d
%0 =
c,
∂t
D0 + kd

x ∈ Γs ,

(15)

where %0 is the density of the oxide film. In the dependence (15) for the local
thickness of the oxide film, two limiting cases can be distinguished. For a small
thickness of the oxide film (d  D0 /k), linear kinetics is realized:
∂d
%0 ≈ kc.
∂t
For large values of the thickness (d  D0 /k), we have parabolic kinetics:
D0
∂d
%0 ≈
c.
∂t
d
To consider the dynamic processes, the system of equations (1), (2), (7),
in addition to the boundary conditions (3)–(6), (8), (9), (14), is supplemented
with equation (15) describing the time–evolution of the thickness of the oxide
film. Also, we must specify the initial conditions for the velocity
u(x, 0) = u0 (x),

x ∈ Ωf ,

as well as for the oxygen concentration and the thickness of the oxide film,
respectively:
c(x, 0) = c0 (x), x ∈ Ωf .
(16)
d(x, 0) = d0 (x),

x ∈ Γs .

(17)

It seems natural to simplify the hydrodynamic problem and assume that the
fluid flow is steady–state. In view of this, instead of (1), we solve the equation
%u · ∇u + ∇p − µ∇2 u = 0.

(18)

Thus, modeling of the oxidation process is conducted via solving the initial–
boundary value problem (2)–(9), (14)–(18).

6

2.3. Dimensionless problem
Let us formulate the above problem in the dimensionless form, using for the
dimensionless velocity, pressure and concentration the same notation as for the
dimensional ones. As the reference values, we take the diameter of tubes l, the
inlet velocity value ū and the inlet concentration c̄. Then equation (18) takes
the form
1 2
∇ u = 0,
(19)
u · ∇u + ∇p −
Re
where
%lū
Re =
µ
is the Reynolds number. The boundary condition (3) takes the form
u · n = −1,

u × n = 0,

x ∈ Γin .

(20)

Taking into account the fact that the pressure is determined to within a constant,
rewrite (4) in the form
p − σn · n = 0,

σn × n = 0,

x ∈ Γout .

(21)

In the dimensionless variables, we get
σ = ∇u + (∇u)T .
The boundary conditions (5), (6) remain unchanged.
The equation (7) in the dimensionless form is written as
∂c
1 2
+ ∇(uc) −
∇ c = 0,
∂t
Pe
where
Pe =

x ∈ Ωf ,

t > 0,

(22)

lū
D

is the Peclet number.
From (8), we have
c = 1,

x ∈ Γin ,

(23)

and the boundary condition (9) takes the form
∇c · n = 0,

x ∈ Γsym ∪ Γout .

(24)

The boundary condition (14) in the dimensionless form is written as following:

−1
1
1
∇c · n +
+
d
c = 0, x ∈ Γs .
(25)
Sh1
Sh2
Here Sh1 is the Sherwood number that corresponds to the oxidation process
based on linear kinetics:
kl
Sh1 = .
D
7

For the Sherwood number corresponding to the approximation of parabolic kinetics, we have
l D0
Dc̄
Sh2 =
, d¯ =
.
D d¯
%0 ū
Using the above reference value d¯ for the thickness of the oxide film, from (15),
we arrive at the dimensionless equation
∂d
=
∂t



1
1
+
d
Sh1
Sh2

−1
c,

x ∈ Γs .

(26)

Thus, the problem under consideration is characterized by two mass transfer
parameters, namely, Re and Pe, as well as the oxidation process parameters Sh1
and Sh2 . It should be noted that the boundary conditions for the oxidant
concentration are non–linear (see (25), (26)).
3. Hydrodynamic processes
The numerical solution of the 2D hydrodynamic problem of a cross-flow
around tube bundles is obtained on the basis of finite–element discretization in
space.
3.1. Computational domain and grids
The triangulation of the computational domain Ωf is performed using the
Gmsh grid generator (website gmsh.info, [25]). Scripts for preparing geometries
of both configurations are written in the Python programming language.
In the in–line configuration, we consider 10 circular tubes with the dimensionless diameter equals 1 arranged in the bundle with the longitudinal and
transverse pitches (measured between tube centers) equal 2. The grid independence study was conducted using the sequence of refined grids presented in
Fig. 3. The staggered configuration has the same pitches, the corresponding
computational grids are shown in Fig. 4.
3.2. Computational algorithm
The hydrodynamic problem was solved separately from mass transfer. The
finite–element discretization [26] is based on the variational formulation for the
boundary value problem (19), (2), (5), (6), (20), (21). For the velocity u, we
define the function space V (u ∈ V ):
V = {u ∈ H 1 (Ωf ) : u · n = −1, u × n = 0 on Γin ,
u = 0 on Γs , u · n = 0 on Γsym }.
For test functions v ∈ V̂ , we have
V̂ = {v ∈ H 1 (Ωf ) : v = 0 on Γin , v = 0 on Γs , v · n = 0 on Γsym }.

8

9

Figure 3: Computational grids for the in–line configuration: a — basic (medium) grid (23,493 vertices and 45,328 cells), b — coarse grid (8,370
vertices and 15,764 cells), c — fine grid (71,326 vertices and 139,760 cells)

c

b

a

10

Figure 4: Computational grids for the staggered configuration: a — basic (medium) grid (23,511 vertices and 45,366 cells), b — coarse grid (8,362
vertices and 15,750 cells), c — fine grid (71,246 vertices and 139,600 cells)

c

b

a

For the pressure p and the corresponding test function q, we set p, q ∈ Q, where
Q = {q ∈ L2 (Ωf ) : q = 0 on Γout }.
Multiply equation (17) by v and integrate it over the whole domain. The similar
transformation of equation (2) is performed using q. Taking into account the
boundary conditions (6), (20), (21), we obtain the system of equations:
a(u, v) − b(v, p) = 0 ∀v ∈ V̂ ,

(27)

b(u, q) = 0 ∀q ∈ Q

(28)

for the desired u ∈ V , p ∈ Q. Here
Z
Z
1
∇u · ∇v dx,
(u · ∇u) · v dx +
a(u, v) :=
Re Ωf
Ωf
Z
b(v, q) :=
(∇ · v)q dx.
Ωf

To define the finite–element discretization, we choose finite–dimensional subspaces Vh ⊂ V , V̂h ⊂ V̂ and Qh ⊂ Q for the approximate solution and test
functions. Here we use the Taylor–Hood P2 − P1 finite element [27]. It consists of a continuous P2 Lagrange element for the velocity components and a
continuous P1 Lagrange element for the pressure field.
To solve the nonlinear variational problem, the iterative Newton method is
applied. The computational implementation is based on the FEniCS platform
for solving partial differential equations (website fenicsproject.org, [28, 29]). The
convergence of the iterative Newton method for various values of the Reynolds
number is presented in Table 1.
3.3. Numerical results for the stationary hydrodynamic problem
We start our numerical study with Re = 10. The velocity components and
pressure calculated on the basic grid for the in–line configuration are shown in
Fig. 5. Similar data for the staggered configuration are given in Fig. 6.
The grid convergence of the numerical solution demonstrates Fig. 7 for the
in-line configuration and Fig. 8 for the staggered configuration, respectively. In
these figures, the velocity components and the pressure (u1 , u2 , p) are given versus x1 along the midline of the computational domain. These main plots are obtained on the fine grid. In addition, there is depicted the deviation (δu1 , δu2 , δp)
from this finest solution. These two deviations corresponding to the coarse and
basic (medium) grids are normalized (multiplied by 100) in order to make this
visualization more evident. These figures indicate a good accuracy of the numerical results obtained on the fine grid.
The effect of the Reynolds number on the flow is shown in Fig. 9 for the in–
line configuration and in Fig. 10 for the staggered configuration, respectively.
For Re < 1, in fact, we can restrict ourselves to the Stokes approximation (do
not take into account the convective terms in equation (19)). To show flow
11

12

Figure 5: Velocity and pressure for the in–line configuration, Re = 10: a — the longitudinal velocity component u1 , b — the transverse velocity
component u2 , c — the pressure

c

b

a

13

Figure 6: Velocity and pressure for the staggered configuration, Re = 10: a — the longitudinal velocity component u1 , b — the transverse velocity
component u2 , c — the pressure

c

b

a

14

Figure 7: Velocity components and the pressure along the midline of the computational domain for the in–line configuration presented as the solution
on the fine grid and the normalized deviations from it for the basic and coarse grids: a — the longitudinal velocity component u1 , b — the transverse
velocity component u2 , c — the pressure

c

b

a

15

Figure 8: Velocity components and the pressure along the midline of the computational domain for the staggered configuration presented as the
solution on the fine grid and the normalized deviations from it for the basic and coarse grids: a — the longitudinal velocity component u1 , b — the
transverse velocity component u2 , c — the pressure

c

b

a

Table 1: Convergence of the iterative process

Re

10

50

150

Iteration
1
2
3
4
1
2
3
4
1
2
3
4

Absolute residual
in–line/staggered
1.427e-01/1.874e-01
9.574e-03/3.373e-02
6.448e-06/4.395e-05
1.885e-11/2.581e-10
1.427e-01/1.874e-01
3.531e-02/8.926e-02
7.390e-04/8.966e-03
4.896e-06/2.646e-04
1.427e-01/1.874e-01
4.975e-02/9.916e-02
3.886e-03/9.707e-02
3.283e-04/4.778e-02

Relative residual
in–line/staggered
1.858e-02/2.439e-02
1.246e-03/4.392e-03
8.395e-07/5.722e-06
2.454e-12/3.360e-11
1.858e-02/2.439e-02
4.597e-03/1.162e-02
9.621e-05/1.167e-03
6.374e-07/3.445e-05
1.858e-02/2.439e-02
6.477e-03/1.291e-02
5.059e-04/1.264e-02
4.274e-05/6.220e-03

patterns for these 2D stationary flows, we employ the streamfunction ψ defined
by the relation


∂ψ
∂ψ
u=
,−
.
∂x2
∂x1
It is evaluated from the known velocity using the equation
−∇2 ψ = ω,
where
ω=

x ∈ Ωf ,

∂u2
∂u1
−
∂x1
∂x2

and the corresponding Dirichlet boundary conditions are specified on ∂Ωf . Figures 11 and 12 present streamlines at different Reynolds numbers for the in–line
and staggered configurations, respectively.
4. Oxidation process
The oxidation process is modeled via solving the unsteady convection–diffusion
equation with the corresponding conditions on the tube surfaces.
4.1. Solution of the mass transfer problem
The unsteady problem of oxygen transport (20)–(26) with the initial conditions (16), (17) is solved numerically using the Lagrangian finite elements P1 .
Define
S = {s ∈ H 1 (Ωf ) : s = 1 on Γin },
Ŝ = {s ∈ H 1 (Ωf ) : s = 0 on Γin }.
16

17

Figure 9: Velocity components and the pressure along the midline of the computational domain for the in–line configuration at Re = 10, 50, 150: a
— the longitudinal velocity component u1 , b — the transverse velocity component u2 , c — the pressure

c

b

a

18

Figure 10: Velocity components and the pressure along the midline of the computational domain for the staggered configuration at Re = 10, 50, 150:
a — the longitudinal velocity component u1 , b — the transverse velocity component u2 , c — the pressure

c

b

a

19

Figure 11: Streamlines for the in–line configuration at various Reynolds numbers: a — Re = 10, b — Re = 50, c — Re = 150

c

b

a

20

Figure 12: Streamlines for the staggered configuration at various Reynolds numbers: a — Re = 10, b — Re = 50, c — Re = 150

c

b

a

The concentration c ∈ S is obtained from the equation




∂c
∂d
, s + e(c, s) =
,s
∀s ∈ Ŝ,
∂t
∂t
s

(29)

where
Z
e(c, s) := −
Ωf

1
cu · ∇s dx +
Pe

Z

Z
∇c · ∇s dx +

Ωf

(u · n)cs dx,
Γout

Z
(ϕ, s)s := −

ϕs dx.
Γs

For d ∈ G = L2 (Γs ), from (26), we have


∂d
,g
∂t




−
s

1
1
+
d
Sh1
Sh2

!

−1
c, g

= 0,

g ∈ G.

(30)

s

For time–stepping, we employ the Crank–Nicolson scheme of second order
[30, 31]. Let τ be a step–size of a uniform grid in time such that cn = c(tn ), tn =
nτ, n = 0, 1, .... For equation (29), we apply the following two–level scheme:

 n+1
  n+1

 n+1
c
− cn
d
− dn
c
− cn
n = 0, 1, .... (31)
,s + e
,s =
,s
τ
2
τ
s
Similarly, for (30), we have


dn+1 − dn
,g
τ




−
s

1
1 dn+1 + dn
+
Sh1
Sh2
2

−1

cn+1 + cn
,g
2

!
= 0,

(32)

s

for the given initial conditions (see (16), (17)). We confine ourselves to the case
of the homogeneous initial conditions:
c0 = 0,

x ∈ Ωf ,

d0 = 0,

x ∈ Γs .

4.2. Linear kinetics of oxidation
First, we consider the oxidation process at the initial stage that is characterized by linear kinetics. In this simplest case, predictions are performed for
the following parameter set:
Re = 50,

Pe = 10,

Sh1 = 0.001,

Sh−1
2 = 0.

Time–integration process is considered up to the moment T = 50 with the time–
step τ = 0.1. The average oxygen concentration at the outlet is treated as the
integral characteristic of the process:
R
c(x, t)dx
R
.
(33)
cout (t) = Γout
dx
Γout
21

22

Figure 13: Distribution of the concentration c for the in–line configuration at various time–moments: a — t = 5, b — t = 10, c — t = 15

c

b

a

23

Figure 14: Distribution of the concentration c for the staggered configuration at various time–moments: a — t = 5, b — t = 10, c — t = 15

c

b

a

24

Figure 15: Concentration c along the midline of the computational domain for the in–line configuration presented as the solution obtained with
τ = 0.05 and the normalized deviations from it δc obtained with τ = 0.1 and τ = 0.2 at various time–moments: a — t = 5, b — t = 10, c — t = 15

c

b

a

25

Figure 16: Concentration c along the midline of the computational domain for the staggered configuration presented as the solution obtained with
τ = 0.05 and the normalized deviations from it δc obtained with τ = 0.1 and τ = 0.2 at various time–moments: a — t = 5, b — t = 10, c — t = 15

c

b

a

The concentration distribution at different time–moments is shown in Figs. 13
and 14 for the in–line and staggered configuration, respectively. Figures 15 and
16 demonstrate the influence of the time–step on the solution for two configurations. In these figures, the concentration c is shown along the midline of the
computational domain. These main plots are obtained on the finest grid in time
corresponding to the time step τ = 0.05. For a comparison, there is depicted the
deviation δc from this solution. These two deviations obtained with τ = 0.1 and
τ = 0.2 are normalized (multiplied by 100) in order to make this comparison
more evident. It is easy to see that the time–grid with τ = 0.1 provides a good
enough accuracy for the time–integration.
The dependence of the average concentration at the outlet cout defined according to (33) on the Peclet number is shown in Fig. 17. The influence of the
Sherwood number is given in Fig. 18.

Figure 17: Influence of Pe: left — in–line configuration, right — staggered configuration

Figure 18: Effect of Sh1 : left — in–line configuration, right — staggered configuration

It is interesting to analyze the time–evolution of the local thickness of the
oxide film on tubes. The distribution of the film thickness along tube surfaces is
presented in the dependence on the local angle θ defined in Fig. 19. The value
θ = 0 corresponds to the leading edge of the tube, whereas θ = π corresponds
to the rear edge. The film thickness on the surface of the third (from the inlet)
26

θ

Figure 19: Local angle θ

Figure 20: Distribution of the oxide film thickness on the 3-rd tube surface: left — in–line
configuration, right — staggered configuration

tube is shown in Fig. 20 at various time–moments. We observe a higher growth
rate of the film at the leading edge in comparison with the rear edge. Next, for
the staggered configuration, the maximum growth of the film takes place at the
leading edge. For the in–line configuration, the maximum growth is near the
midpoint of the tube surface. At high time-moments, the growth rate of the
film tends to be constant.

27

4.3. Effect of parabolic kinetics
In the general case, both linear and parabolic kinetics (Sh−1
2 > 0) are taken
into account. Under these conditions, the problem at the new time level (see
(31), (32)) is nonlinear. To solve it, Newton’s method is applied. As a rule,
two or three iterations are necessary for the solution convergence in calculations
presented here.

Figure 21: Distribution of the oxide film thickness for Sh2 = 10−5 : left — in-line configuration,
right — staggered configuration

Figure 22: Distribution of the oxide film thickness for Sh2 = 10−6 : left — in-line configuration,
right — staggered configuration

The impact of Sh2 is presented in Figs. 21–23. These calculations are performed for Sh1 = 0.001. As Sh2 decreases, the growth rate of the oxide film
also decreases. This is due to the effect of the film thickness on the oxidation
process.
It is interesting to study the time–evolution of the oxide mass on an individual tube. For the first five tubes in the in–line and staggered configurations
with the boundaries Γs,i , we calculate
Z
mi (t) = 2
d(x, t)dx, i = 1, 2, ..., 5.
Γs,i

28

Figure 23: Distribution of the oxide film thickness for Sh2 = 10−7 : left — in-line configuration,
right — staggered configuration

The dependence of the oxide mass on the number Sh2 is shown in Figs. 24–26.
There is no essential differences between the in–line and staggered configurations.

Figure 24: Time–evolution of the oxide mass for the first five tubes at Sh2 = 10−5 : left —
in-line configuration, right — staggered configuration

29

Figure 25: Time–evolution of the oxide mass for the first five tubes at Sh2 = 10−6 : left —
in-line configuration, right — staggered configuration

Figure 26: Time–evolution of the oxide mass for the first five tubes at Sh2 = 10−7 : left —
in-line configuration, right — staggered configuration

5. Conclusions
A new mathematical model of the oxidation process is developed for a cross–
flow around tube bundles. It is based on the incompressible Navier–Stokes equations for laminar flows. The oxidant transport is described by the convection–
diffusion equation. The peculiarity of the developed mass transfer model consists in the formulation of boundary conditions for the oxidant. They include
explicitly the thickness of the oxide film and involve both linear and parabolic
kinetics of oxidation.
The computational algorithm is based on finite–element discretization in
space using Lagrangian finite elements on triangular grids. Hydrodynamics is
calculated independently of mass transfer in the stationary formulation. The
Newton method is applied to solve governing equations in the primitive variables. The Crank-Nicolson scheme is used for time–stepping in solving equations
for the oxygen concentration and oxide film thickness.

30

Two–dimensional predictions of stationary laminar flows in in–line and staggered tube bundles are conducted for different values of the Reynolds number.
The grid independence study is performed using a sequence of refined grids.
Calculations of the oxidation processes are carried out for the initial stage
of oxidation in the linear kinetics approximation. A parametric study is done
to indicate the influence of the Peclet and Sherwood numbers. For the complete model that involves both linear and parabolic kinetics, the mass transfer
phenomenon is investigated numerically with emphasis on the time–evolution
of the local thickness of the oxide film on tube surfaces.
Acknowledgements
The work of the fourth author was supported by the grant of the Russian
Federation Government (# 14.Y26.31.0013).
References
[1] W. M. Kays, A. L. London, Compact Heat Exchangers, McGraw-Hill, New
York, 1964.
[2] T. L. Bergman, A. S. Lavine, F. P. Incropera, D. P. DeWitt, Fundamentals
of Heat and Mass Transfer, John Wiley & Sons, 2011.
[3] A. Žukauskas, Advances in Heat Transfer 8 (1972) 93–160.
[4] C. Iwaki, K. H. Cheong, H. Monji, G. Matsui, Experiments in Fluids 37
(2004) 350–363.
[5] S. S. Paul, M. F. Tachie, S. J. Ormiston, International Journal of Heat and
Fluid Flow 28 (2007) 441–453.
[6] S. B. Beale, D. B. Spalding, Journal of Fluids and Structures 13 (1999)
723–754.
[7] Y. Q. Wang, L. A. Penner, S. J. Ormiston, Numerical Heat Transfer: Part
A: Applications 38 (2000) 819–845.
[8] A. M. F. El-Shaboury, S. J. Ormiston, Numerical Heat Transfer, Part A:
Applications 48 (2005) 99–126.
[9] S. Jayavel, S. Tiwari, International Journal of Numerical Methods for Heat
& Fluid Flow 19 (2009) 931–949.
[10] Y. T. K. Gowda, B. S. V. P. Patnaik, P. A. A. Narayana, K. N. Seetharamu,
International Journal of Heat and Fluid Flow 19 (1998) 49–55.
[11] F. Zdravistch, C. A. Fletcher, M. Behnia, International Journal of Numerical Methods for Heat & Fluid Flow 5 (1995) 717–733.

31

[12] M. N. Dhaubhadel, J. N. Reddy, D. P. Telionis, International Journal for
Numerical Methods in Fluids 7 (1987) 1325–1342.
[13] A. S. Wilson, M. K. Bassiouny, Chemical Engineering and Processing: Process Intensification 39 (2000) 1–14.
[14] W. A. Khan, J. R. Culham, M. M. Yovanovich, International Journal of
Heat and Mass Transfer 49 (2006) 4831–4838.
[15] K. Lam, Y. F. Lin, L. Zou, Y. Liu, International Journal of Heat and Fluid
Flow 31 (2010) 32–44.
[16] N. Benarji, C. Balaji, S. P. Venkateshan, Heat and Mass Transfer 44 (2008)
445–461.
[17] A. Horvat, B. Mavko, Numerical Heat Transfer, Part A: Applications 49
(2006) 699–715.
[18] Y. Q. Wang, P. L. Jackson, Journal of Thermophysics and Heat Transfer
24 (2010) 534–543.
[19] T. A. Tahseen, M. Ishak, M. M. Rahman, Renewable and Sustainable Energy Reviews 43 (2015) 363–380.
[20] L. D. Landau, E. Lifshitz, Fluid Mechanics, Pergamon Press, Oxford, 1987.
[21] T. Li, N. G. Deen, J. A. M. Kuipers, Chemical Engineering Science 60
(2005) 1837–1847.
[22] D. J. Young, High Temperature Oxidation and Corrosion of Metals, Elsevier, Amsterdam, 2016.
[23] B. E. Deal, A. S. Grove, Journal of Applied Physics 36 (1965) 3770–3778.
[24] Z. Xu, K. M. Rosso, S. M. Bruemmer, The Journal of Chemical Physics
135 (2011) 024108–1–024108–7.
[25] C. Geuzaine, J.-F. Remacle, International Journal for Numerical Methods
in Engineering 79 (2009) 1309–1331. doi:10.1002/nme.2579.
[26] P. M. Gresho, R. L. Sani, Incompressible Flow and the Finite Element
Method, Volume 2, Isothermal Laminar Flow, Wiley, Chichester, 2000.
[27] C. Taylor, P. Hood, Computers & Fluids 1 (1973) 73–100. doi:10.1016/
0045-7930(73)90027-3.
[28] A. Logg, K.-A. Mardal, G. N. Wells, (Eds.), Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book,
Springer-Verlag, Berlin, 2012. doi:10.1007/978-3-642-23099-8.

32

[29] M. S. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg,
C. Richardson, J. Ring, M. E. Rognes, G. N. Wells, Archive of Numerical Software 3 (2015) 9–23. doi:10.11588/ans.2015.100.20553.
[30] A. A. Samarskii, The Theory of Difference Schemes, Marcel Dekker, Inc.,
New York, 2001.
[31] U. M. Ascher, Numerical Methods for Evolutionary Differential Equations,
SIAM, Philadelphia, PA, 2008.

33

