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 (Kays and London, 1964; Bergman et al., 2011). In addition to large–scale experimental studies (Žukauskas, 1972; Iwaki et al., 2004; Paul et al., 2007), computational technologies are in common use for theoretical considerations of key aspects of heat and mass transfer processes in tube bundles (Beale and Spalding, 1999; Wang et al., 2000; ElShaboury and Ormiston, 2005; Jayavel and Tiwari, 2009).
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, Gowda et al. (1998); Zdravistch et al. (1995); Dhaubhadel et al. (1987); Wilson and Bassiouny (2000); Khan et al. (2006); Lam et al. (2010); Benarji et al. (2008); Horvat and Mavko (2006); Wang and Jackson (2010); Tahseen et al. (2015) 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 (Landau and Lifshitz, 1987). Numerical simulation of flows around tube bundles in the presence of mass transfer is performed (Li et al., 2005) similarly to the calculation of the temperature. Various oxidation models were developed (Young, 2016) 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., Deal and Grove (1965); Xu et al. (2011)) 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 is the vertical (longitudinal) coordinate and is the horizontal (transverse) coordinate. The tube rows of a bundle can be either in–line or staggered in the direction of the fluid velocity. They form a periodic structure of cylinders of the same circular crosssection. 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 of these rectangles is occupied by the fluid, whereas the part corresponds to the tubes. The solid tube boundaries are denoted by and stands for the symmetry boundaries. The fluid with oxygen enters into the subdomains through the boundary and outflows on the boundary . Mass transfer processes are considered in the computational domain .
2.1 Hydrodynamic processes
To describe a crossflow in a tube bundle, the incompressible Navier–Stokes equations are used:
(1) 
(2) 
where and are the velocity and pressure of the fluid, respectively, while and are, respectively, the viscosity and density assumed to be constant.
Suitable boundary conditions are imposed on . These are the normal and tangential velocity components and/or components of the force applied to the boundary, which are determined by the expression . Here
is the viscous stress tensor and
stands for the outer unit normal to the boundary. The uniform velocity profile with the value is specified at the inlet:(3) 
The pressure and condition for the absence of tangential forces are prescribed at the outlet:
(4) 
The no–slip and no–permeability conditions are imposed on the rigid walls of tubes:
(5) 
On the symmetry boundaries, we put the slip conditions:
(6) 
2.2 Oxygen transport
We denote the oxygen concentration in the fluid by , 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
(7) 
where is the solute diffusion coefficient, which is assumed to be scalar and constant.
A known concentration of oxygen is specified at the inlet:
(8) 
where is constant. We prescribe zero flux of the solute at the outlet and symmetry boundaries as follows:
(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 (Young, 2016).
Following Deal and Grove (1965); Xu et al. (2011), we distinguish three stages of the oxygen transport:

oxygen transport to the external surface of the oxide film;

oxygen diffusion through the oxide film towards the oxide–metal interface;

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 , and at the oxide–metal interface as (see Fig. 2). Here 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
(10) 
where is the outer unit normal for the fluid domain.
For the mass flux in the oxide film, we have
(11) 
where 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
(12) 
where is the oxidation rate (the first–order chemical reaction).
For the quasi–stationary regime, we have
(13) 
From the second equality in (13), in view of (11), (12), we can eliminate . 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:
(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:
(15) 
where 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 (), linear kinetics is realized:
For large values of the thickness (), we have parabolic kinetics:
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
as well as for the oxygen concentration and the thickness of the oxide film, respectively:
(16) 
(17) 
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 , the inlet velocity value and the inlet concentration . Then equation (18) takes the form
(19) 
where
is the Reynolds number. The boundary condition (3) takes the form
(20) 
Taking into account the fact that the pressure is determined to within a constant, rewrite (4) in the form
(21) 
In the dimensionless variables, we get
The boundary condition (14) in the dimensionless form is written as following:
(25) 
Here is the Sherwood number that corresponds to the oxidation process based on linear kinetics:
For the Sherwood number corresponding to the approximation of parabolic kinetics, we have
Using the above reference value for the thickness of the oxide film, from (15), we arrive at the dimensionless equation
(26) 
3 Hydrodynamic processes
The numerical solution of the 2D hydrodynamic problem of a crossflow 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 is performed using the Gmsh grid generator (website gmsh.info, Geuzaine and Remacle (2009)). 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 (Gresho and Sani, 2000) is based on the variational formulation for the boundary value problem (19), (2), (5), (6), (20), (21). For the velocity , we define the function space ():
For test functions , we have
For the pressure and the corresponding test function , we set , where
Multiply equation (17) by and integrate it over the whole domain. The similar transformation of equation (2) is performed using . Taking into account the boundary conditions (6), (20), (21), we obtain the system of equations:
(27) 
(28) 
for the desired , . Here
To define the finite–element discretization, we choose finite–dimensional subspaces , and for the approximate solution and test functions. Here we use the Taylor–Hood finite element (Taylor and Hood, 1973). It consists of a continuous Lagrange element for the velocity components and a continuous 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,
Logg et al. (2012); Alnæs et al. (2015)). The convergence of the iterative Newton method for various values of the Reynolds number is presented in Table 1.Re  Iteration  Absolute residual  Relative residual 

in–line/staggered  in–line/staggered  
1  1.427e01/1.874e01  1.858e02/2.439e02  
10  2  9.574e03/3.373e02  1.246e03/4.392e03 
3  6.448e06/4.395e05  8.395e07/5.722e06  
4  1.885e11/2.581e10  2.454e12/3.360e11  
1  1.427e01/1.874e01  1.858e02/2.439e02  
50  2  3.531e02/8.926e02  4.597e03/1.162e02 
3  7.390e04/8.966e03  9.621e05/1.167e03  
4  4.896e06/2.646e04  6.374e07/3.445e05  
1  1.427e01/1.874e01  1.858e02/2.439e02  
150  2  4.975e02/9.916e02  6.477e03/1.291e02 
3  3.886e03/9.707e02  5.059e04/1.264e02  
4  3.283e04/4.778e02  4.274e05/6.220e03 
3.3 Numerical results for the stationary hydrodynamic problem
We start our numerical study with . 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 inline configuration and Fig. 8 for the staggered configuration, respectively. In these figures, the velocity components and the pressure () are given versus along the midline of the computational domain. These main plots are obtained on the fine grid. In addition, there is depicted the deviation () 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 , in fact, we can restrict ourselves to the Stokes approximation (do not take into account the convective terms in equation (19)). To show flow patterns for these 2D stationary flows, we employ the streamfunction defined by the relation
It is evaluated from the known velocity using the equation
where
and the corresponding Dirichlet boundary conditions are specified on . 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 . Define
The concentration is obtained from the equation
(29) 
where
For , from (26), we have
(30) 
For time–stepping, we employ the Crank–Nicolson scheme of second order Samarskii (2001); Ascher (2008). Let be a step–size of a uniform grid in time such that . For equation (29), we apply the following two–level scheme:
(31) 
Similarly, for (30), we have
(32) 
for the given initial conditions (see (16), (17)). We confine ourselves to the case of the homogeneous initial conditions:
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:
Time–integration process is considered up to the moment
with the time–step . The average oxygen concentration at the outlet is treated as the integral characteristic of the process:(33) 
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 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 . For a comparison, there is depicted the deviation from this solution. These two deviations obtained with and are normalized (multiplied by 100) in order to make this comparison more evident. It is easy to see that the time–grid with provides a good enough accuracy for the time–integration.
The dependence of the average concentration at the outlet defined according to (33) on the Peclet number is shown in Fig. 17. The influence of the Sherwood number is given in Fig. 18.
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 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) 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 timemoments, the growth rate of the film tends to be constant.
4.3 Effect of parabolic kinetics
In the general case, both linear and parabolic kinetics () 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.
The impact of is presented in Figs. 21–23. These calculations are performed for . As 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 , we calculate
The dependence of the oxide mass on the number is shown in Figs. 24–26. There is no essential differences between the in–line and staggered configurations.
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 CrankNicolson scheme is used for time–stepping in solving equations for the oxygen concentration and oxide film thickness.
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
 Kays and London (1964) W. M. Kays, A. L. London, Compact Heat Exchangers, McGrawHill, New York, 1964.
 Bergman et al. (2011) T. L. Bergman, A. S. Lavine, F. P. Incropera, D. P. DeWitt, Fundamentals of Heat and Mass Transfer, John Wiley & Sons, 2011.
 Žukauskas (1972) A. Žukauskas, Advances in Heat Transfer 8 (1972) 93–160.
 Iwaki et al. (2004) C. Iwaki, K. H. Cheong, H. Monji, G. Matsui, Experiments in Fluids 37 (2004) 350–363.
 Paul et al. (2007) S. S. Paul, M. F. Tachie, S. J. Ormiston, International Journal of Heat and Fluid Flow 28 (2007) 441–453.
 Beale and Spalding (1999) S. B. Beale, D. B. Spalding, Journal of Fluids and Structures 13 (1999) 723–754.
 Wang et al. (2000) Y. Q. Wang, L. A. Penner, S. J. Ormiston, Numerical Heat Transfer: Part A: Applications 38 (2000) 819–845.
 ElShaboury and Ormiston (2005) A. M. F. ElShaboury, S. J. Ormiston, Numerical Heat Transfer, Part A: Applications 48 (2005) 99–126.
 Jayavel and Tiwari (2009) S. Jayavel, S. Tiwari, International Journal of Numerical Methods for Heat & Fluid Flow 19 (2009) 931–949.
 Gowda et al. (1998) 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.
 Zdravistch et al. (1995) F. Zdravistch, C. A. Fletcher, M. Behnia, International Journal of Numerical Methods for Heat & Fluid Flow 5 (1995) 717–733.
 Dhaubhadel et al. (1987) M. N. Dhaubhadel, J. N. Reddy, D. P. Telionis, International Journal for Numerical Methods in Fluids 7 (1987) 1325–1342.
 Wilson and Bassiouny (2000) A. S. Wilson, M. K. Bassiouny, Chemical Engineering and Processing: Process Intensification 39 (2000) 1–14.
 Khan et al. (2006) W. A. Khan, J. R. Culham, M. M. Yovanovich, International Journal of Heat and Mass Transfer 49 (2006) 4831–4838.
 Lam et al. (2010) K. Lam, Y. F. Lin, L. Zou, Y. Liu, International Journal of Heat and Fluid Flow 31 (2010) 32–44.
 Benarji et al. (2008) N. Benarji, C. Balaji, S. P. Venkateshan, Heat and Mass Transfer 44 (2008) 445–461.
 Horvat and Mavko (2006) A. Horvat, B. Mavko, Numerical Heat Transfer, Part A: Applications 49 (2006) 699–715.
 Wang and Jackson (2010) Y. Q. Wang, P. L. Jackson, Journal of Thermophysics and Heat Transfer 24 (2010) 534–543.
 Tahseen et al. (2015) T. A. Tahseen, M. Ishak, M. M. Rahman, Renewable and Sustainable Energy Reviews 43 (2015) 363–380.
 Landau and Lifshitz (1987) L. D. Landau, E. Lifshitz, Fluid Mechanics, Pergamon Press, Oxford, 1987.
 Li et al. (2005) T. Li, N. G. Deen, J. A. M. Kuipers, Chemical Engineering Science 60 (2005) 1837–1847.
 Young (2016) D. J. Young, High Temperature Oxidation and Corrosion of Metals, Elsevier, Amsterdam, 2016.
 Deal and Grove (1965) B. E. Deal, A. S. Grove, Journal of Applied Physics 36 (1965) 3770–3778.
 Xu et al. (2011) Z. Xu, K. M. Rosso, S. M. Bruemmer, The Journal of Chemical Physics 135 (2011) 024108–1–024108–7.
 Geuzaine and Remacle (2009) C. Geuzaine, J.F. Remacle, International Journal for Numerical Methods in Engineering 79 (2009) 1309–1331. doi:10.1002/nme.2579.
 Gresho and Sani (2000) P. M. Gresho, R. L. Sani, Incompressible Flow and the Finite Element Method, Volume 2, Isothermal Laminar Flow, Wiley, Chichester, 2000.
 Taylor and Hood (1973) C. Taylor, P. Hood, Computers & Fluids 1 (1973) 73–100. doi:10.1016/00457930(73)900273.
 Logg et al. (2012) A. Logg, K.A. Mardal, G. N. Wells, (Eds.), Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book, SpringerVerlag, Berlin, 2012. doi:10.1007/9783642230998.
 Alnæs et al. (2015) 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.
 Samarskii (2001) A. A. Samarskii, The Theory of Difference Schemes, Marcel Dekker, Inc., New York, 2001.
 Ascher (2008) U. M. Ascher, Numerical Methods for Evolutionary Differential Equations, SIAM, Philadelphia, PA, 2008.