Effects of Mesh Generation on Modelling Aluminium Anode Baking Furnaces

Abstract

Anode baking is critical in carbon anode production for aluminium extraction. Operational and geometrical parameters have a direct impact on the performance of anode baking furnaces (ABF), and hence on the resulting anode quality. Gas flow patterns, velocity field, pressure drop, shear stress and turbulent dissipation rate are the main operational parameters to be optimised, considering a specific geometry that is discretised as a mesh. Therefore, this paper aims to establish the need to generate an appropriate mesh to perform accurate numerical simulations of three-dimensional turbulent flow in a single section of an ABF. Two geometries are considered for generating three meshes, using COMSOL and cfMesh, with different refinement zones. The three meshes are used for creating nine incompressible isothermal turbulent flow models, with varying operational parameters. Velocity field, convergence and turbulent viscosity ratio in the outlet of fuel inlet pipes are the quantification criteria. Quantification criteria have shown that a better physical representation is obtained by refining in the whole combustion zone. COMSOL Multiphysics’ built-in mesh generator allows quadrilateral, tetrahedron and hexahedron shapes. Adaptive cell sizes and shapes have a place within modelling, since refining a mesh in appropriate zones brings the Peclet number down when the incompressible isothermal turbulent flow is simulated.

1. Introduction

Aluminium anodes are the significant components in the extraction process of aluminium from bauxite ore at 15%, and are therefore considered essential [1]. After being baked (heat-treated) in open-top ring-type furnaces, referred to as Anode Baking Furnaces (ABF), anodes reveal particular mechanical, thermal and electrical properties that determine their suitability for aluminium production [2]. However, this process utilises vast amounts of energy and releases undesired gases, such as NOx [3]. Therefore, anode baking requires optimisation to ensure reduced NOx emission, energy consumption and soot-free combustion while improving anode quality [4].

Baking optimisation can be carried out by tuning operational and geometrical parameters. However, in situ research is expensive, disruptive, challenging and time-consuming. Instead, numerical mathematical modelling may be considered for this purpose. Computational fluid dynamics (CFD) is the preferred technique for modelling.

CFD modelling of ABF is based on a typical ABF installation. That is: fuel is injected according to a specified time-dependent law into burners located at the top of the furnace module. The rate and frequency at which fuel is fed into the burner depend on the specific burner’s characteristics, and injection frequency is higher for high-speed burners. Combustion air enters into the furnace through air inlets. Fluid is convected and diffused by the general incompressible isothermal turbulent flow. The main goal is to ensure flow uniformity inside the furnace for the correct baking of anodes, considering subsequent phenomena.

CFD modelling uses a software-assisted design of the furnace at hand. This representation is called geometry. Since computing has a discrete nature, the designed geometry has to be divided into unit controls (elements) to perform calculations of the phenomena at each specific location. The result of the division is a discrete grid or mesh, which is an input to solve Partial Differential Equations (PDEs). PDEs, representing the flow of a fluid, are well known in the industry and academia for their high complexity, due to the presence of non-linear terms.

An adequate solution of the PDEs, with minimal error, requires a mesh with small element (cell) sizes and shapes. However, how fine a mesh has to be remains an open problem [4,5,6,7]. Since each control unit (element or cell)—in all the Navier–Stokes equations—has to be solved, a fine mesh increases the computational load. This paper aims to establish the need for an appropriate mesh generation to perform accurate numerical simulations of three-dimensional turbulent flow in a single section of an ABF. Thus, two geometries are considered for generating three meshes—two meshes using cfMesh version 3.2.1 [8], and one mesh using COMSOL Multiphysics version 5.5 [9]—with different refinement zones and cell shapes and sizes. The three meshes are used for simulations using nine isothermal turbulent flow models, varying operational parameters with a nested model scheme in the last three models. Moreover, COMSOL Multiphysics version 5.5 [9] is used for solving the Navier–Stokes κ–ε turbulence model with the finite element method and a linear solver within each Newton Raphson iteration. Simulations’ accuracy is validated with the lowest error reached, along with the turbulent viscosity ratio and physical interpretation by experts.

Initially, a naive mesh, without taking advantage of symmetry, is used for illustrating the increment of computational load even with a coarse mesh—when considering the entire length of the ABF in the z-axis. Then, a second mesh, taking advantage of symmetry by halving the z-axis and a refinement in the fuel jet stream zone, is used to perform numerical simulations of the three-dimensional turbulent flow in a single section of an ABF. Finally, a third mesh, taking advantage of symmetry by halving the z-axis and a refinement in the whole combustion zone, is used to perform numerical simulations of the three-dimensional turbulent flow in a single section of an ABF. Thus, mesh 1 is a naïve mesh with no refined zones, in a testing phase, and generated using cfMesh. Mesh 2 is refined in the fuel jet stream zone and generated using cfMesh with proper (user-informed) settings. Mesh 3 is refined in the whole combustion zone and generated using COMSOL Multiphysics. Moreover, an artificial diffusion scheme is used for accelerating convergence; however, it solves a nearby problem. To overcome oscillations in the Newton-Raphson method around the optimal value, a nested scheme is used by setting the results of one computation as the initial value of another in the Newton Raphson method. As an experimental evaluation, nine models are created to evaluate the need for an appropriate mesh generation to perform accurate numerical simulations of the three-dimensional turbulent flow in a single section of an ABF.

The first four models use mesh 1. Although the numerical simulations converge, the velocity streams do not represent the three-dimensional incompressible isothermal turbulent flow in a single section of the ABF. Thus, the use of symmetry and an appropriate mesh—refining in combustion zones—are needed.

Models 5 and 6 use mesh 2. Convergence is achieved using the lowest contribution of artificial diffusion in model 6. However, artificial diffusion may solve a nearby problem. Models 5 and 6 showed that there is not enough refining in the fuel jet stream zone for representing the incompressible isothermal turbulent flow.

Models 7 and 8 use mesh 3. Model 8 has as initial values the results of model 7, in order to provide a better initialisation to the Newton Raphson method. Models 7 and 8 showed that a better physical representation is obtained by refining the whole combustion zone.

Model 9 uses mesh 2. Model 9 has as initial values the results of model 8, in order to provide a better initialisation to the Newton Raphson method, unlike model 5. It shows that there is not enough refining in the fuel jet stream zone for correctly representing the incompressible isothermal turbulent flow.

The nested scheme in model 9—initialised with the results of model 8, which was initialised with the results of model 7—did not yield the expected convergence, as it reached the lowest error at 10⁻². However, the nested scheme in model 8—it was initialised with the results of model 7—reached convergence at 10⁻³. Meshes were created with a smaller cell size only in critical areas, and a larger cell size in the rest of the geometry areas with the aim of reducing the computational load. Refining a mesh in the appropriate zones brings the Peclet number down. After observing the dissipation results, mesh 2—which is refined in the fuel jet stream zone—adequately models the split stream at the meeting of the fuel jet stream with the first tie-brick, in the incompressible isothermal turbulent flow. Regarding convergence and the turbulent viscosity ratio, mesh 3—which is refined in the whole combustion zone—adequately models the incompressible isothermal turbulent flow at the outlet of the fuel inlet pipes in the combustion zone.

Convergence, velocity field and turbulent viscosity ratio in the outlet of the fuel inlet pipes are the quantification criteria. Mesh 3, generated with a refinement in the whole combustion zone, provides a better physical representation of the incompressible isothermal turbulent flow in a single section of an ABF. Mesh 3 converged after around 50 iterations. Different cell shapes, in a mesh, may align better with a convective turbulent flow, especially in the most critical areas, such as the combustion zone. Thus, an adaptive mesh, sizes and shapes have a place when modelling a single section of an ABF to perform numerical simulations of the three-dimensional turbulent flow.

Mesh generation is a fundamental task for simulating isothermal turbulence dissipation, that should be done mainly in the whole combustion zone for achieving convergence and representing an isothermal turbulent flow according to reality—such as the flow splitting when meeting the tie-bricks and avoiding penetrating downwards, following a homogeneous direction of flow across the subsections of an ABF by distributing the flow into the whole single section.

2. Related Works

The effects of the operational and geometrical parameters on the performance of the furnace have been addressed in the literature by mathematical simulation. Specific aspects of simulations are reported in [10,11,12]. Moreover, the mathematical modelling of the anode baking process has been developed and improved significantly in the past years. In 1983, Bui et al. simulated a horizontal flue ring furnace in which they treated the furnace as a counter-flow heat exchanger [13]. Many models that have been developed at the later stage are based on these early developed models. A more detailed 3D modelling of the ABF started in the mid-90s. Kocaefe et al. presented a model in which a commercial CFD code, CFDS-FLOW3D, was used for solving the governing differential equations [14]. However, this model used simplified combustion and radiation models, failing to comment on the pollutants.

Severo and Gusberti established boundary conditions to be able to properly bake all brands of raw materials that may be expected [15]. Moreover, they developed a user-friendly software to analyse furnace energy efficiency and the minimum oxygen concentration in different sections [1]. However, the tool cannot be considered for obtaining more specific data related to soot or NOx formation with higher accuracy.

Ping et al. [16] have reported the effect of baffles and tie-brick arrangements on the flow characteristics of the anode baking process. From their report, baffle and tie-brick positioning has a significant impact on flow homogeneity. Ordronneau et al. [17] demonstrated the necessity of employing different simulation tools for meeting the challenge of increasing anode baking furnace productivity.

Other studies have explored deformations of geometry. For instance, Baiteche et al. [18] studied the effect of flue-wall deformation on anode temperature distribution. Comparing the temperature profiles on a line in the pit transverse direction for straight and deformed flue-walls, it was observed that after flue-wall deformation, the temperature profile is no longer symmetric, which indicates a non-uniform baking process. Later on, Zaidani et al. [19] have also studied the effect of flue-wall deformation on anode temperature distribution.

From their experiments, Kocaefe et al. [20] have provided an enhanced physical understanding of ABF performance. An enhanced performance may be reached using different computational tools with different levels of complexity. In this way, the κ–ε incompressible isothermal turbulent flow model has a nonlinear behaviour. It represents, in a way, a “worst-case” with demands to nonlinear iterations [21]. Valen-Sendstad et al. have shown that the Reynolds Averaged turbulent Navier–Stokes equations can be solved by a Newton Raphson iterative process after finite element discretisation with the distinct advantage of a superlinear convergence over traditionally used SIMPLE-based approaches.

A dynamic process model was developed by Oumarou et al. to investigate the effect of temperature variation in the vertical component by considering a vertical component of flue gas [22,23,24]. This allows the 2D temperature distribution boundary condition for the pit sub model. However, the model fails to optimise reducing emissions, saving energy and maintaining anode quality.

A similar work is carried out by Tajik et al., in which the effect of flue-wall design on the flow field, combustion and temperature has been modelled in Ansys Fluent [25,26,27]. The finite volume method is used in Ansys Fluent, whereas COMSOL^® Multiphysics is based on the finite element method (FEM). Comparing the two approaches, Ansys Fluent—finite volume method—and COMSOL^® Multiphysics—finite element method—may give an insight on the problem solution. In conclusion, vast modelling approaches are developed for anode baking furnace. However, the model for NOx reduction still requires significant attention.

Chaodong et al. [28] have used an FEM-based model with the aim of developing a large-scale, high-efficiency and energy-saving baking furnace. They reported results with two designs that have been optimised for the flue-wall and exhaust ramp. Gaoui et al. [29] have examined the influence of baffles. The idea was to remove the baffles and simulate using a finite element method. Some pitfalls were found on the road, like uneven heat distribution, too fast degassing or flue-wall pinching.

Nakate et al. [6] have developed a model that focuses on reducing NOx emissions. First, they modelled the incompressible isothermal turbulent flow. Then, the model has been extended by adding combustion reaction [4]. As a third step, heat transfer has been considered in the general model. However, the software used is constrained by the basic eddy dissipation model. They concluded that it is necessary to have detailed combustion models based on a probability density function. Later, Nakate et al. reported a discussion using meshes generated by COMSOL Multiphysics as input to the finite element method [4,5]. They used the turbulent viscosity ratio to demonstrate a physically unrealistic computation from the COMSOL mesh. Although Nakate et al. [4,5,6] addressed the reduction of NOx emissions, heat transfer was not considered in the models, and meshes were refined in combustion zones, in this paper.

A bottom-up study was performed by Talice [7]. Firstly, models were simulated using 2D geometry. This implies no consideration of the z-component. Modelling in two dimensions allows for a familiarisation with the model and the extraction of general flow features inside the furnace. Modelling in two dimensions allows for the description of turbulence behaviour in one planar surface. Talice used the Spalart Allmaras one-equation model. Later, Talice developed the model using three dimensions [30]. Moreover, Talice analysed the fluid flow using a more realistically represented geometry at the burner zone. Unlike [7], the two-equation Realizable κ–ε model was considered. Standard Wall Functions were used at all the solid walls. In this study, Talice concluded that using two dimensions entails a high contribution of the z-component to a well-described fluid flow [30]. Talice proposed that the conflict between the physical behaviour of the flow and the prescribed uniform value for pressure can have detrimental numerical effects. Numerical impact may occur when slowing down or simply making it impossible to reach the full convergence of the numerical method. For that reason, Talice proposed that outlet zones be redesigned to ensure flow uniformity.

Additionally, there are reports about other physical phenomena. Grégoire et al. [31] conducted a comparative study on two different model approaches for ABF combustion modelling. Later on, Grégoire and Gosselin reported a comparison of three combustion models for ABF using Ansys Fluent [32].

Table 1. Summary of literature review on ABF design modelling.

Authors	Year	Objectives	Combustion Model	Detailed Kinetics	Radiation Model
Ping. et al.	2002	Influence of the baffles on the flowing field	Non-reactive flow	Not included	Not included
Severo et al.	2005	Developing a 3D CFD model for flue-wall design modification	EDM	Not included	P1
Ordronneau et al.	2006	Application of CFD simulation for crossover design off-gas cleaning system optimisation training purposes	Not specified	Not specified	Not specified
Gregoire et al.	2011	Comparison of two modelling approaches to predict variability	Hot air jet approximation	Not specified	DO method
Kocaefe et al.	2013	Different modelling approaches on anode baking furnace	Not mentioned	Not included	Not specified
Baiteche et al.	2015	Effects of flue-wall deformation and employing different radiation models	Empirical kinetic expression	Not included	- P1 - Monte Carlo method
Ghaui et al.	2016	Implementation of baffleless flue-wall technology	Not mentioned	Not included	Not specified
Zaidani et al.	2017	Effects of flue-wall deformation	Non-reactive flow	Not included	Not specified
Chaodong et al.	2018	Optimisation and development of the furnace structures, process parameters and firing control system	Not specified	Not specified	Not specified
Nakate et al.	2018	Develop a mathematical 2D model to reduce NOx emissions considering turbulent flow, combustion model and radiation	EDM	- κ–ε - Spalart Allmaras	- P1 - DO
Talice	2018	Develop a 2D model to analyse flow behaviour	Not used	Spalart Allmaras	Not used
Nakate et al.	2019	Develop a 3D model to analyse flow behaviour	Not used	κ–ε	Not used
Talice	2019	Develop a 3D model to analyse flow behaviour	Not used	κ–ε	Not used
Nakate et al.	2021	Establish an analysis in 3D flow with a high rate of fuel injection	Energy equation	- Standard κ–ε - Realizable κ–ε	Not used

provides a summary of the recent literature. The closest publication is Nakate et al. [6], that focused on reducing the NOx and considered the heat contribution to the velocity field. Although Nakate et al. [6] used mesh 2—for simulating a non-isothermal turbulent flow as the main objective—that achieved convergence and showed that by increasing temperature the density is reduced—implying a deeper penetration of the jet—in this paper mesh 2 is used for simulating an isothermal turbulent flow model.

3. Model Description

In this section, mathematical fundamentals of the incompressible non-reactive isothermal turbulent flow are presented. More details can be found in Wilcox [33].

Standard κ–ε Turbulence Model

The following equations are solved for the six unknown parameters (pressure p, velocity at each component: $\overline{u_1}$, $\overline{u_2}$ and $\overline{u_3}$; the first transported variable as the turbulent kinetic energy κ and the second transported variable as the rate of dissipation of turbulent kinetic energy ε) during the incompressible isothermal flow computation:

$$\frac{\partial \overline{u_j}}{\partial x_j}=0,$$

(1)

$$\overline{u_j}\frac{\partial \overline{u_i}}{\partial x_j}=-\frac{\partial \overline{p}}{\partial x_i}+\frac{\partial }{\partial x_j}\left[{\nu }_{eff}\left(\frac{\partial \overline{u_i}}{\partial x_j}+\frac{\partial \overline{u_j}}{\partial x_i}\right)\right],fori,j=1,2,3$$

(2)

$$\overline{u_j}\frac{\partial k}{\partial x_j}=\frac{\partial }{\partial x_j}\left[{\nu }_{eff}^k\frac{\partial k}{\partial x_j}\right]{+\nu }_T\frac{\partial \overline{u_i}}{\partial x_j}\left(\frac{\partial \overline{u_i}}{\partial x_j}+\frac{\partial \overline{u_j}}{\partial x_i}\right)-\epsilon ,$$

(3)

$$\overline{u_j}\frac{\partial ϵ}{\partial x_j}=\frac{\partial }{\partial x_j}\left[{\nu }_{eff}^{ϵ}\frac{\partial ϵ}{\partial x_j}\right]{+C}_{1ϵ}\frac{ϵ}{k}{\nu }_T\frac{\partial \overline{u_i}}{\partial x_j}\left(\frac{\partial \overline{u_i}}{\partial x_j}+\frac{\partial \overline{u_j}}{\partial x_i}\right){-C}_{2ϵ}\frac{{ϵ}^2}{k},$$

(4)

$${\nu }_{eff}{=\nu +\nu }_T,$$

(5)

$${\nu }_{eff}^k=\nu +\frac{{\nu }_T}{{\sigma }_k},$$

(6)

$${\nu }_{eff}^{ϵ}=\nu +\frac{{\nu }_T}{{\sigma }_{ϵ}},$$

(7)

$${\nu }_T{=C}_{\mu }\frac{k^2}{ϵ+{ϵ}_s},$$

(8)

$${\sigma }_k=1,{\sigma }_{ϵ}=1.3,C_{\mu }=0.09,C_{1ϵ}=1.44\mathrm{and}{C}_{2ϵ}=1.92.$$

(9)

The laminar viscosity μ is calculated using Sutherlands’s law ${\mu =A}_s\frac{T^{3/2}}{{T+T}_s}$, where $A_s=1.67212\times {10}^{-6}$ and $T_s=170.672$ are constants.

The over-bar denotes the ordinary Reynolds averaging. The following notation is used:

• ρ is the density of the fluid (SI units: kg/m³).
• μ is the dynamic viscosity of the fluid (Pa·s or N·s/m² or kg/(m·s)).
• ν is the kinematic viscosity of the fluid (m²/s).
• ε is a small number added to avoid the division by zero.
• σk and σε are the turbulent Prandtl numbers for κ and ε.
• $\overline{p}=\overline{p}+\frac{2}{3}\rho k$ and the $\frac{1}{\rho }$ factor in front of the pressure term in the RANS equations are dropped. Then, if the true mean pressure field is sought, one has to take this into consideration.
• The default values of the model constants, ${\sigma }_k{,\sigma }_{ϵ}{,C}_{\mu }{,C}_{1ϵ}{,\mathrm{and}C}_{2ϵ},$ have been determined from experiments with air and water for fundamental turbulent shear flows, including homogeneous shear flows and decaying isotropic grid turbulence. They have been found to work fairly well for a wide range of wall-bounded and free shear flows.
• Although the default values of the model constants are the standard ones, and the most widely accepted, one can change them (if needed).

The standard κ–ε turbulence model is focused on the velocity field, as well as the turbulent flow. Hence, an isothermal flow model is used in the next section.

4. Model Configurations

Nine models are created with varying parameters as described in

Table 2. Parameters used for the different models.

Parameter	Model 1	Model 2	Model 3	Model 4	Model 5	Model 6	Model 7	Model 8	Model 9
Fluid properties
Density (k/m3)	1.2	1.2	1.2	1.2	1.2	1.2	1.2	1.2	1.2
Dynamic viscosity (Pa·s)	8.9 × 10−4	1.8 × 10−5	8.9 × 10−4	1.8 × 10−5	1.8 × 10−5	1.8 × 10−5	1.8 × 10−5	1.8 × 10−5	1.8 × 10−5
Initial values for Newton’s iteration
Ux (m/s)	70	0	70	0	0	0	0	0	0
Uy (m/s)	0	0	0	0	0	0	0	0	0
Uz (m/s)	0	0	0	0	0	0	0	0	0
Pressure (Pa)	0	0	0	0	0	0	0	0	0
Boundary conditions
Wall	No slip
Inlet	Fully developed flow
Outlet	Pressure
Geometry and mesh
Geometry	1	1	1	1	2	2	2	2	2
Mesh	1	1	1	1	2	2	3	3	2
Mesh generator tool	cfM	cfM	cfM	cfM	cfM	cfM	COMS	COMS	cfM
Artificial diffusion scheme
δ(u,p)	Off	Off	Off	Off	Off	0.005	0.25	Off	Off
δ(κ,ε)	Off	Off	Off	Off	Off	Off	0.25	Off	Off
Results used as initial value for the Newton Raphson method
Initial value	0	0	0	0	0	δ(u,p) = 0.01	0	Model 7	Model 8

. The models are implemented using COMSOL^® Multiphysics version 5.5 for solving the Navier–Stokes κ–ε turbulence model with the finite element method [9]. All solver parameters are set as default except for the linear solver. GMRES (as the Krylov subspace method), Algebraic Multigrid (as preconditioner) and Vanka (as pre- and post-smoother within Algebraic Multigrid) are selected for the linear solver.

4.1. Geometry and Mesh

Two geometries are used for representing a single section of an ABF, in Figure 1. The first geometry has a full representation of the z-axis without any symmetry consideration, in Figure 1a; and the second geometry takes advantage of symmetry using half of the z-axis, in Figure 1b, for reducing computational load. Figure 2 illustrates the fuel inlet pipes in the z-axis of both geometries. It can be observed, in Figure 2a, the fuel inlet pipe at the centre of geometry 1 at z = 0.54 m, whilst in Figure 2b the fuel inlet pipe is at the symmetry plane at z = 0.27 m. These geometries are used for generating three meshes.

Two meshes were created using cfMesh version 3.2.1 [8] by Prajakta Nakate, one for each geometry, and the third one was created using COMSOL Multiphysics version 5.5 [9] and geometry 2. A fourth coarse mesh was created using COMSOL Multiphysics version 5.5 and geometry 2, included in the Appendix A. The fourth mesh has a refinement at the fuel jet stream zone. However, the fourth mesh is coarser than meshes 2 and 3 and only has tetrahedron and triangle cells.

Mesh 1 is externally generated using cfMesh 3.2.1 [8] by Prajakta Nakate—without user input settings; and it is not intended to be used for ABF simulations—and geometry 1; it has no refinement at the combustion zone, and is shown in Figure 3a. Mesh 2 is externally generated using cfMesh 3.2.1 [8] by Prajakta Nakate with a proper (user-informed) setting and geometry 2; it has a refinement in the fuel jet stream zone and is shown in Figure 3b. Both mesh 1 and mesh 2 are represented using quadrilateral, tetrahedron, pyramid, prism, hexahedron and triangle cells. Mesh 3 was created with COMSOL Multiphysics version 5.5 [9] and geometry 2, has a refinement in the whole combustion zone, is represented using tetrahedron and triangle cells and is shown in Figure 3c.

Table 3. Description of 3D geometry and the three meshes.

Mesh	Mesh 1	Mesh 2	Mesh 3
Generator	cfMesh	cfMesh	COMSOL
Symmetry	No	Yes	Yes
Length x-axis (m)	5.5	5.5	5.5
Length y-axis (m)	5.0	5.0	5.0
Length z-axis (m)	0.54	0.27	0.27
Location of fuel inlet pipes (z-axis) (m)	0.27	0.27	0.27
Cell shape	Cartesian	Cartesian	Tetrahedral

presents a description of the meshes, along with the 3D geometries, that are used for generating nine models for modelling the incompressible isothermal turbulent flow.

Magnifications of the fuel inlet zones in meshes 2 and 3 are presented in Figure 4. Mesh 2 has a refinement along the fuel jet stream, in Figure 4a, whilst mesh 3 has a refinement at the whole combustion zone—where air inlet and fuel inlet meet—mainly at the leftmost top zone, in Figure 4b.

Skewness penalises cells with large or small angles compared to an ideal cell size and is used as a quality measure, by COMSOL Multiphysics version 5.5 [9], where values close to one are the best, in Figure 5. Mesh 2 has the best average skewness and mesh 3 the worst.

Table 4. Skewness quality measure values.

Mesh	Number of Cells	Minimum Skewness	Average Skewness
Mesh 1	2,424,973	0.00	0.79
Mesh 2	545,694	0.23	0.86
Mesh 3	4,924,080	0.08	0.66

presents a description of each mesh in terms of number of cells, as well as minimum and average skewness. Mesh 2 has a superior mesh quality histogram, two and a half time less cells than mesh 1 and nine time less cells than mesh 3.

4.2. Simulations

Four incompressible isothermal turbulent flow simulations are conducted using mesh 1. Parameters are described in Table 2 as models 1 to 4. Then, three incompressible isothermal turbulent flow simulations are conducted using mesh 2—models 5, 6 and 9. Model 6 used artificial diffusion to achieve convergence, and the diffusion parameter was tuning at δ(u,p) = 0.005. Finally, two incompressible isothermal turbulent flow simulations are conducted using mesh 3—models 7 and 8. Model 7 used artificial diffusion to achieve convergence, and the diffusion parameters were tuning at δ(u,p) = 0.25 and δ(κ,ε) = 0.25.

Meshes are used in a nested manner by setting the results of one computation as the initial value of another. The incompressible isothermal turbulent flow—using mesh 3 and artificial diffusion δ(u,p) = 0.25—was modelled as model 7, achieving convergence at 10⁻³. Then, the incompressible isothermal turbulent flow—using mesh 3 without artificial diffusion—was modelled as model 8, with initial values set from the results of model 7, achieving convergence at 10⁻³. On the other hand, the incompressible isothermal turbulent flow—using mesh 2 without artificial diffusion—were modelled as model 9, with initial values set from the results of model 8, achieving the lowest error at 10⁻².

4.3. Numerical Implementation

The numerical implementation was done using COMSOL Multiphysics version 5.5 [9], in a heterogeneous HPC cluster with Intel^® Xeon^® CPU E5-2630 v3 @2.40GHz x32 cores, and 129GB RAM. Convergence was achieved when errors reached at least 10⁻³. Additionally, simulations’ accuracy were evaluated using the graphical results of the turbulent viscosity ratio plots, defined as μ_T/μ₀, at each element, as well as the physical interpretation by experts and the recent literature [4]. This paper is focused on the physical interpretation of the resulted velocity and the isothermal turbulent behaviour when using different types of mesh cells, refinement and parameters of artificial diffusion. In particular, this paper is focused on the behaviour of these phenomena in the combustion zones—where the air and fuel meet.

4.3.1. Finite Element Method in CFD

The discretisation in the finite element method consists of the following: the solution of partial differential equations requires a discretisation of the computational domain. A visual representation is shown in Figure 6. Some examples of shape element are triangles, quadrangles, tetrahedra, prisms or hexahedra [34].

Each entity is called element or cell. Each element has to have at least another element as a neighbour. Mesh size is an important factor that determines the complexity of the solution to be calculated.

The finite element method has been approached in [35,36] for fluid problems, long time before the effectiveness of implementing FEM in computational calculations began to be evaluated.

FEM has the following advantages:

• It is a very general method,
• There is more facility to increment the element order,
• Physical fields may be reproduced more accurately,
• Physics and mathematics often require different type of functions for a phenomenon. Different phenomena can be represented at the same time with FEM,
• To reach more accuracy, increase order of polynomials and refine the mesh.

As weaknesses, FEM has the requirement to have more computations to represent phenomena at each basic unit [37].

In this work, FEM is used in COMSOL Multiphysics version 5.5 [9]. For that reason, the next subsection presents a brief introduction about FEM.

4.3.2. Theoretical Definition of FEM

Let Ω be the problem domain (i.e., the area limited by the geometry). In 3D, the curvilinear polygon $\Omega \subset {ℝ}^3$ with piecewise analytic boundary

$$\partial \Omega =\underset{i=1}{\overset{M}{{\displaystyle \cup }}}\overline{{\Gamma }_i},$$

(10)

with $\Omega \subset {ℝ}^3$, will be a bounded domain with piecewise analytic boundary $\partial \Omega$, that consists of faces ${\Gamma }_1$, ${\Gamma }_2$, ..., ${\Gamma }_M$, which are curved polygons in${ℝ}^3$, joined by edges ${\gamma }_1{,\text{…},\gamma }_{n_e}$(curves in ${ℝ}^3$) and vertices ${\mathrm{A}}_1$, ${\mathrm{A}}_2$, ..., ${\mathrm{A}}_{nw}$.

The weak formulation of the Navier–Stokes equations is:

Find $\mathbf{u}\in {\mathrm{L}}^2{(\mathrm{R}}^{+}{{[\mathrm{H}}^1(\Omega )]}^3)\cap$ ${\mathrm{C}}^0{(\mathrm{R}}^{+}{{[\mathrm{L}}^2(\Omega )]}^3)$ such that:

$$\{\begin{array}{c}{{\displaystyle \int }}_{\Omega }^{}\rho \frac{\partial \mathbf{u}}{\partial t}\mathrm{v}+{{\displaystyle \int }}_{\Omega }^{}\mu \nabla \mathbf{u}\nabla \mathrm{v}+{{\displaystyle \int }}_{\Omega }^{}\rho (\mathbf{u}·\nabla )\mathbf{u}·\mathbf{v}-{{\displaystyle \int }}_{\Omega }^{}p\nabla ·\mathbf{v}={{\displaystyle \int }}_{\Omega }^{}\mathit{f}·\mathbf{v}+{{\displaystyle \int }}_{{\Gamma }_N}^{}\mathit{h}·\mathbf{v},\forall \mathbf{v}\in \mathrm{V}\\ {{\displaystyle \int }}_{\Omega }^{}q\nabla ·\mathbf{u}=0,\forall q\in \mathrm{Q}.\end{array}$$

(11)

Function h is given Neumann boundaries. f is external force. v and q are test functions in the space V and Q, respectively. Additional details can be found at [34].

The problem has two variables to approximate, velocity u and pressure p. It is known as mixed variational formulation. A solution may be addressed using Lagrange multipliers to determine the value of each variable. However, it is more efficient to use a penalisation model of p, simplifying the discrete problem into an equations system that only depends on u. This system allows us to determine p once calculated u.

In the discretisation problem, using FEM, e will be the set of resulted elements by dividing the domain Ω into the subregions in (10). An approximation of unknowns u = (u, v, w) and p at each element will be given by:

The weak formulation of the Navier–Stokes equations is:

$$\mathrm{u}={\displaystyle \sum }_{\mathrm{i}=1}^{{\mathrm{m}}_{\mathrm{u}}}{\mathrm{N}}_{{\mathrm{u}}_i}{q}_{{\mathrm{u}}_i}{=\mathrm{N}}_{\mathrm{u}}^{\mathrm{T}}{\mathbf{q}}_{\mathrm{u}},$$

(12)

$$\mathrm{v}={\displaystyle \sum }_{\mathrm{i}=1}^{{\mathrm{m}}_{\mathrm{v}}}{\mathrm{N}}_{{\mathrm{v}}_i}{q}_{{\mathrm{v}}_i}{=\mathrm{N}}_{\mathrm{v}}^{\mathrm{T}}{\mathbf{q}}_{\mathrm{v}},$$

(13)

$$\mathrm{w}={\displaystyle \sum }_{\mathrm{i}=1}^{{\mathrm{m}}_{\mathrm{w}}}{\mathrm{N}}_{{\mathrm{w}}_i}{q}_{{\mathrm{w}}_i}{=\mathrm{N}}_{\mathrm{w}}^{\mathrm{T}}{\mathbf{q}}_{\mathrm{w}},$$

(14)

$$\mathrm{p}={\displaystyle \sum }_{\mathrm{i}=1}^{{\mathrm{m}}_{\mathrm{p}}}{\mathrm{N}}_{p_i}{q}_{{\mathrm{p}}_i}{=\mathrm{N}}_p^{\mathrm{T}}{\mathbf{q}}_p,$$

(15)

where vectors ${\mathbf{q}}_{\mathbf{u}}$,${\mathbf{q}}_{\mathbf{v}}$,${\mathbf{q}}_{\mathbf{w}}$ denote local values inside the velocities field; u = (u, v, w), and ${\mathbf{q}}_{\mathit{p}}$ denote local values inside the pressure field. ${\mathbf{N}}_{\mathbf{u}}$, ${\mathbf{N}}_{\mathbf{v}}$,${\mathbf{N}}_{\mathbf{w}}$ and ${\mathbf{N}}_{\mathit{p}}$ are shape functions of the velocity and pressure and the unknown total vector of element e is given by ${\mathbf{q}}_{\mathbf{e}}^{\mathbf{T}}=\left[{\mathbf{q}}_{\mathbf{u}}^{\mathbf{T}},{\mathbf{q}}_{\mathbf{v}}^{\mathbf{T}},{\mathbf{q}}_{\mathbf{w}}^{\mathbf{T}},{\mathbf{q}}_{\mathit{p}}^{\mathbf{T}}\right]$.

Unknowns u = (u, v, w) and p, in the problem, are no longer mathematical functions and become the values of these functions at the nodes. The complete problem solution follows the rules for discrete problems.

5. Results and Discussion

According to quantification criteria convergence, velocity field and turbulent viscosity ratio, in the outlet of the fuel inlet pipes, are the aspects to evaluate in the simulation. A model reaches convergence when the residual error value is around 10⁻³.

Table 5. Lowest error reached by the models.

Model	Model 1	Model 2	Model 3	Model 4	Model 5	Model 6	Model 7	Model 8	Model 9
Lowest error reached (u,p)	10−3	10−3	10−3	10−2	10−1	10−3	10−3	10−3	10−2
Lowest error reached (κ,ε)	10−3	10−3	10−3	10−2	10−3	10−3	10−3	10−3	10−2

shows the lowest error reached by the nine models. Models 4, 5 and 9 did not reach convergence according to the convergence criterium.

Model 1: The velocity field is presented in Figure 7a. The incompressible isothermal turbulent flow model describes the fluid flow appropriately according to the set of parameters considered in model 1. It is a simple model, and a finer mesh is required for a higher velocity and lower viscosity.

Model 2: The velocity field plot is shown in Figure 7b. In the velocity field plot, higher velocity values are modelled. Small streams with high velocity can be observed near the baffles due to the effect of lowering the viscosity. The convergence value suggests that model 2 is feasible when using low viscosities. Nevertheless, the effect of the velocity of the fuel is not considered at this stage.

Models 3 and 4: In the same way, Figure 7c,d shows the velocity field. Model 3 has reached convergence around 10⁻³, whilst model 4 has presented periodic oscillations around 10⁻². There is an incorrect representation of the velocity field in both models. In fact, in Figure 7c, the fuel jet stream penetrates the furnace downwards, avoiding the first tie-brick obstacle. This implies that the flow will not be split and distributed along the furnace. Hence, the second subsection—when the flow goes up—will not have a uniform velocity with respect to the first subsection. This physical behaviour can be explained by the use of a higher viscosity in model 3. Thus, model 4 considers the lowest viscosity and the highest velocity required, as shown in Figure 7d. Nevertheless, the effect of having no refinement at the combustion zone is observed as the fuel jet stream goes to the right side in the leftmost fuel stream. This implies a remaining non-uniform velocity distribution in the subsequent section—where the flow goes up.

When lowering the viscosity and increasing the velocity, turbulent phenomena have long eddies, which represent chaotic flows. This is not in line with what it is expected. Additionally, a refinement is required in the areas near the fuel inlet pipes and in the convergence between the two flows from the left inlet pipe and the air inlets.

Model 5: the segregated solver convergence plot is presented in Figure 8, showing the turbulent variables and the κ–ε variables. The model has not reached the desired convergence of 10⁻³, after 400 iterations. A lower viscosity and a higher velocity are required along with using a refined mesh that may increment the computer load. Additionally, the convergence plots showed repeated oscillations without signs of reaching a minimum error. This indicates that oscillations may occur [38] due to an incorrect meshing of the combustion zone. The lowest error reached was 10⁻¹ for the fluid flow variables and 10⁻³ for the turbulence variables. Mesh 2 has to be refined in the whole combustion zone for modelling the isothermal turbulent flow.

Model 6: artificial diffusion was used to achieve convergence. The isotropic diffusion model was used with the parameter δ = 0.005, as the lowest value for velocity-pressure variables. Since artificial diffusion added linearity to the model, convergence may be forced. Artificial diffusion represents a false diffusion that is not in accordance with the actual isothermal turbulent flow. The effect of artificial diffusion can be graphically observed in Figure 9a, that shows the velocity field magnitude at the symmetry plane. Convergence was reached with the lowest error of 10⁻³. Figure 9b shows the wall resolution in viscous units. Figure 9c shows the residual of velocity field calculation, and Figure 9d illustrates the turbulent viscosity ratio. Low residuals can be observed from the velocity field calculation. Nevertheless, the fuel inlet velocity is low, but not according to the incompressible isothermal turbulent flow. In particular, the turbulent viscosity ratio showed a dissipation in the right side of the combustion zone in Figure 9d. Therefore, the fuel jet stream does not penetrate downwards properly, as observed in Figure 9d. The fluid at the farthest down locations of the section of the ABF has not the desired velocities. Despite the above, the artificial diffusion scheme is an alternative to cases of high complexity and high computational load. Thus, the results of a simulation with artificial diffusion can be used as the initial value of a simulation without artificial diffusion [39]. This allows the solver to reach a solution starting from a value closer to the solution, as the Newton Raphson method assumed.

Model 7: artificial diffusion is used with the aim of achieving convergence. Figure 10 presents the velocity slide plot. The velocity plot shows that there are thin streams with high velocity. Few high-velocity streams are produced by two factors: (1) intersection of faces—as error message—was obtained when refining the whole combustion zone when using the mesh generation techniques inside COMSOL Multiphysics; (2) low convection due to the presence of false-added diffusion. Thus, the physical interpretation indicates that the solution corresponds to a nearby problem due to the use of artificial diffusion.

Model 8: The results of model 7 were used as initial values to model 8—without any artificial diffusion. Figure 11 shows the velocity field Figure 11a and the segregated convergence plot showed a convergence around 50 iterations, in Figure 11b, where thick high-velocity streams can be observed. In Figure 11a, more uniform velocity streams are produced by the calculations. Fuel jet streams are captured more appropriately, with respect to the last models—model 7. However, there is a remaining velocity jet stream that is not captured by the effect of splitting after the first tie-brick into the left part by model 8.

On the other hand, there is a stronger high-velocity fuel stream in the right part, in Figure 11a, that is better represented by model 8. Thus, it produces a better flow velocity uniformity in the second and the fourth subsections of the ABF—where the flow goes up.

Model 9: uses mesh 2 and sets the results of model 8 as the initial values for the Newton Raphson method. The lowest error value reached is 2 × 10⁻² for the velocity-pressure variables and 5 × 10⁻³ for κ, ε. Figure 12 shows the velocity slide plot at the iteration with the lowest error value reached and the convergence solver plot. There is an increment of the velocity in all locations compared to model 8. Figure 12a shows the flow stream penetrating with high velocities in ever deeper locations. After 400 iterations, model 9 has not achieved convergence, with a minimum error of 10⁻². Using an approximate solution as the initial value, without false linearity, is still useful for complex studies.

Figure 13 shows a comparison between the obtained velocity field results using mesh 2 and mesh 3. A better representation of the velocity field is yielded by mesh 3.

Appendix B shows a comparison among results from different simulations using different values of the isotropic diffusion parameter to enhance the discussion.

The turbulent viscosity ratio is used to enhance the discussion about the interpretation of results. Figure 14 shows a comparison between the turbulent viscosity ratio at the outlet of the fuel inlet pipes using meshes 2 and 3. A better physical representation of the incompressible isothermal turbulent flow is obtained by refining the whole combustion zone, mesh 3.

As a summary: models 1 to 4 use the non-symmetric mesh 1, that cannot be refined because the computational load increases dramatically. They are based on increasing fuel inlet pipe velocities and decreasing viscosity. This implied a higher computational load along with oscillations in the Newton Raphson method. Although the numerical simulations converged, the velocity streams do not represent the incompressible isothermal turbulent flow. Thus, the use of symmetry and an appropriate mesh refined in combustion zones are needed.

Models 5 and 6 use the symmetric mesh 2, which is refined along the fuel jet stream. Convergence is achieved using the lowest contribution of artificial diffusion, in model 6. Models 5 and 6 showed that there is not enough to refine along the fuel jet stream zone for representing an incompressible isothermal turbulent flow.

Models 7 and 8 use the symmetric mesh 3, which is refined in the whole combustion zone. Model 7 uses artificial diffusion for achieving convergence. Model 8 has as initial values results of model 7 in order to provide a better initialisation to the Newton Raphson method. Model 8 showed that it is possible to obtain a correct physical representation of the incompressible isothermal turbulent flow by refining the whole combustion zone. Although mesh 3 provided a better representation of the outlet of the fuel inlet pipe zones, it did not represent well the turbulence dissipation in other combustion zones.

Model 9 uses the symmetric mesh 2, which is refined along the fuel jet stream zone. Model 9 has as initial values the results of model 8, in order to provide a better initialisation to the Newton Raphson method, without reaching the convergence criterium. However, model 9 correctly represented the dissipation of turbulent flow. The large variety of cell shapes, the majority of which are hexahedral, may align the mesh in the direction of the turbulent flow. However, the fuel jet stream zone refinement is insufficient when modelling an incompressible isothermal turbulent flow.

6. Conclusions

In this paper, nine models were created in order to evaluate the need for an appropriate mesh generation to perform accurate numerical simulations of the three-dimensional incompressible isothermal turbulent flow in a single section of an ABF. Convergence and turbulent viscosity ratio, in the outlet of fuel inlet pipes, are the quantification criteria, along with a physical representation of the incompressible isothermal turbulent flows.

Models 8 and 9 allowed us to compare meshes 2 and 3, and mesh 2 was used in model 9, whilst mesh 3 was used in model 8. The comparison of mesh 2 and mesh 3 shows that mesh 3 has nine times the number of cells, but an adequate refinement zone underneath the whole combustion zone showed that mesh 3 is better than mesh 2.

However, the built-in COMSOL Multiphysics mesh generation tool is not easy to use for constructing an appropriate mesh. The lack of sufficient versatile algorithms for mesh generation inside of COMSOL Multiphysics may limit the use of different cell shapes with lower sizes in non-trivial structures when there are difficult zones to represent corners and edges, as well as the handling of face intersection errors.

Although, mesh 2 has better skewness statistics than mesh 3, the latter is optimal according to the quantification criteria to perform accurate numerical simulations of the three-dimensional incompressible isothermal turbulent flow in a single section of an ABF.

In this way, a mesh, with refinement in the whole combustion zone using a variety of cell sizes and shapes, is better for modelling the three-dimensional incompressible isothermal turbulent flow in a single section of an ABF.