The effect of operating conditions on the nitric oxide formation in anode baking furnace through numerical modeling

Thermal nitric-oxide (NOx) formation in industrial furnaces due to local overheating is a 1 widely known problem. Various industries made significant investments to reduce thermal NOx by 2 varying the operating conditions and designs of the furnace. Finding optimal operating conditions or 3 design parameters by experimenting in the furnace, however, is difficult. Numerical modeling can 4 provide significant information in such cases. In this paper, a three dimensional steady state finite 5 element model for the anode baking industrial furnace is discussed. The COMSOL Multiphysics 6 software is used for modeling the non-premixed turbulent combustion and the conjugate heat transfer 7 to the insulation lining. The mesh generation using the cfMesh software allows to increase the spatial 8 resolution locally at the outlet of the fuel nozzles while maintaining the overall quality of the mesh. 9 The temperature and species mass fraction obtained from the finite element model are calibrated by 10 adjusting the amount of artificial diffusion in the transport equations for the species. The simulated 11 temperature agrees well with the measured data from our industrial partner in regions distant from 12 the flames. The model underestimates the measured oxygen mass fraction. The spatial gradients 13 in oxygen mass fraction, however, are captured well by the model. The effects of variation of the 14 fuel mass flow rate and the fuel pipe diameter on the NOx generation are studied. The results show 15 that by decreasing the fuel mass flow rate and increasing the fuel pipe diameter by 45%, the peak in 16 thermal NOx ppm generated in the furnace decreases by 42%. 17

to analyse the temperature in the furnace. This can provide significant information on the generation 40 of NOx. The anode baking process practiced in Aluchemie furnace is described in the previous paper 41 [13] that describes the aerodynamics in the furnace. The study on NOx emissions continues in this geometry of the model for the two studies remains the same.

97
The meshing of the geometries with two diameters of the fuel pipes are carried out with the 98 cfMesh software. In our previous paper, the analysis of different meshing techniques are discussed 99 [13]. The less diffusive Cartesian mesh for the complex geometry of the anode baking furnace can be 100 obtained using cfMesh software. The region under the fuel outlet is refined further for better resolution 101 of the flow in that region. This region is particularly of interest with respect to NOx generation. Figure   102 2 shows the meshes of the geometries with the two diameters. It can be observed that the region of the 103 jet development is refined in the same manner for both geometries. The difference is only with respect 104 to the diameter of the fuel pipe. Table 2 provides the details of the two meshes.  The quality of the mesh is measured by the skewness parameter. Obtaining a mesh of similar mesh 106 quality for two diameter of fuel pipes is difficult with the default mesher of COMSOL Multiphysics.
107 Figure 3 shows the mesh quality histograms in terms of skewness for different diameters of fuel pipe. 108 It can be observed that in some regions of the 13 mm diameter of fuel pipe there are more elements 109 that are of low quality as compared to the 9 mm diameter. However, the overall histogram quality is 110 comparable. cfMesh software provides control over the size of elements in the region of interest and therefore, obtaining meshes with the external software is preferred.

117
The term K is defined by Equation (2).

119
The standard k − model is used for closing the Reynolds stresses from the RANS equation. The The production term P k is as presented in Equation (5).
The definition of turbulent viscosity µ T with the standard k − model is given in Equation (6).
134 Table 3. Values of the constants for standard k-model. The transport equation for each of the species is given by Equation (9) and (10).

Constant Value
The isotropic diffusion can be added to the diffusion coefficient from Equation (10). This additional 145 diffusion can be varied using a tuning parameter δ id from Equation (11). In a later section, the effect of 146 variation of this tuning parameter is studied.
The reaction term from Equation (9) is modeled using the eddy dissipation model.The

149
regularization is implemented to ensure that the reactant is consumed only when its mass fraction 150 value is higher than zero. Moreover, the mass fraction of the product is restricted to take the maximum 151 value of 1. The eddy dissipation model equations are given by Equation (12) to (16).
155 156 r f or r f or The combustion of methane translates into the generation of heat energy. This acts as a heat source .
The heat generated in the gas domain is conducted through walls to anodes. In this model, a single 178 brick layer is considered so as to account for the heat flux conducting to a combined solid phase.
The radiation in the gas domain is accounted by the radiative source term as presented in Equation

181
(20). The incident radiation, G, is computed by the P1 approximation model as defined by the Equation
The diffusion coefficient of the P1 approximation is defined based on the absorption and scattering 187 coefficient as presented in Equation (22). The scattering is neglected by assigning value zero to σ s .
The radiative flux through boundaries is defined by Equations (23) and (24). The emissivity of the    The rate constants for these reactions are provided in Table 4.  .

217
The details of the discretization as well as the solver settings for the turbulent flow are described [13].

218
A linear type of elements are considered for discretizing all transport equations of chemical species as 219 well as radiative transport equation using the finite element method.  for each of the segregated step is also mentioned in Table 5.   The velocity for the 13 mm diameter fuel pipe model are decided based on the equivalent fuel 238 mass flow rate from 9 mm diameter. Table 7 shows the values of velocities for the two fuel pipe 239 diameters that yields equal mass flow rate of fuel injected in the furnace. As can be expected, for the 240 larger diameter of fuel pipe, the injected fuel inlet velocity is lower. design.

248
The model of the 9 mm pipe diameter is developed systematically in this study. The detailed  After developing the initial flow results, the transport of chemical species is included in the model.

260
The turbulent chemistry interaction is modeled by the eddy dissipation model. This model is based 261 on the 'mixed is burnt' phenomena. This can be observed from the mass fraction distribution of CH 4 262 (reactant) and mass fraction distribution of CO 2 (product) as shown in Figure 5 and 6, respectively.

263
The layer of CH 4 stream that is in contact with the O 2 shows the appearance of the product (CO 2 ).

264
The reaction zone is limited to the mixing zone of the reactants. This produces sharp gradients of the  to be negligible.

279
The heat energy released by the combustion reaction is modeled using the enthalpy of reaction.

280
The enthalpy of reaction is further calculated based on the stoichimetric coefficient and the enthalpy of where the CO 2 is higher. It can be observed that the central part of the jet is colder than the outer part.  The maximum temperature in the furnace is as high as 2470 0 C in case of no radiative heat transfer.

290
Moreover, the lower temperature calculated by the model is unrealistic and as low as -11 0 C. Therefore, 291 heat transfer by radiation is important phenomena in the anode baking furnace that needs to be 292 incorporated in the modeling. The radiation is modeled by the P1 approximation model. As can be 293 observed from Figure 8, the temperature distribution becomes more uniform by accounting radiation.

294
The maximum temperature is lowered to 2225 0 C due to the heat transfer. The high temperature region 295 is restricted to a layer of combustion reaction zone in the jet.  The volume percentage of chemical species such as O 2 and CO 2 are carried out at several locations in 310 the furnace using Testo 350 flue gas analyser. Figure 11 shows the location at which the measurements 311 are carried out. The volume percentage of chemical species are extracted at the same locations.  Figure 11. Location of the measurement points in the furnace. Aluchemie measured volume percentage of chemical species at the given locations. The results from the test case model motivates to study the effect of different tuning of the This effect can also be visualized from the mass fraction of CO 2 in Figure 14. The appearance of CO 2 is 348 limited in the case of no additional diffusion. While the CO 2 is spread out as we increase the diffusion 349 in the transport equations.

350
The right choice of δ id is required in order to obtain comparable results with the measurements.

351
As discussed in the earlier section, the volume percentage of chemical species such as O 2 and CO 2 352 are compared with the measured values at the locations shown in Figure 11. Table 11  if we vary either the X or Y coordinate. This observation is valid while using additional diffusion.

368
As discussed earlier, the effect of additional diffusion causes more mixing. In this paper, the eddy 369 dissipation model is used for combustion. Therefore, an increase in the mixing of the two streams 370 results in the increase in combustion process. The increase in combustion implies that more reactant 371 is consumed producing more product. This can be validated from the comparisons given in Table   372 11 and 12. As we increase the value of δ id (thereby increasing diffusion), the value of O 2 decreases. 5. The diffusion improves as the tuning parameter increases from 0 to 0.5. This leads to wider reaction zones with increase in the tuning parameter. Therefore, the lower mass fraction of O 2 at the outlet is not restricted to a narrow stream for higher tuning parameters as observed from measurements.  The diffusion improves as the tuning parameter increases from 0 to 0.5. This leads to wider reaction zones with increase in the tuning parameter. Therefore, the higher mass fraction of CO 2 at the outlet is not restricted to a narrow stream for higher tuning parameters as observed from measurements. In the earlier section, the variation in the Peclet number is based on the changes in the diffusion 388 coefficient. The decrease in the Peclet number by increasing the diffusion coefficient resulted in 389 increased mixing of fuel and oxidizer stream. In this section the effect of variation on the fuel jet 390 velocity is studied. Figure 15 shows the mass fraction distribution of CH 4 at the symmetry plane for 391 varying fuel inlet velocities. It can be seen from the figure that for higher velocity, the CH 4 does not 392 readily react with oxidizer in case of 0.0037 kg/s. This can be attributed to the reduced mixing due 393 to higher momentum as well as the increased quantity of fuel. For a fixed diameter of fuel pipe, the 394 increase in velocity also accounts for the increase in the fuel quantity. For a fixed oxidizer concentration, 395 it takes longer to consume all CH 4 . Therefore, the length of CH 4 jet increases with increase in the 396 velocity.
397 Figure 15. The CH 4 mass fraction distribution at the symmetry plane for varying fuel inlet mass flow rate of (a) 0.002 kg/s (b) 0.003 kg/s and (c) 0.0037 kg/s with 9 mm fuel pipe diameter. The higher mass flow rate results in the higher momentum of the jet as well as higher amount of CH 4 mass fraction. Therefore, the region of unreacted CH 4 is higher with increased mass flow rate. Figure 16 shows the temperature distribution at the YZ plane cutting through the first burner. As way that the fuel amount remains the same. This is carried out by increasing the fuel pipe diameter.

412
The diameter of the fuel burner pipe is increased from 9 mm to 13 mm. The velocity of the fuel inlet 413 through the 13 mm fuel pipe is such that the mass flow rate of fuel is equal to that through 9 mm 414 diameter. Earlier in Table 7 Figure 18. The comparison shows that the turbulent intensity for the 9 mm diameter is 424 higher, especially in the encircled region below the fuel outlet. Due to the higher turbulence produced 425 by the jet of 9 mm diameter, there is higher mixing of the fuel and oxidizer streams. Therefore, the 426 CH 4 from the fuel jet readily reacts with O 2 from the oxidizer stream.
427 Figure 17. Velocity magnitude [m/s] comparison at the symmetry plane with the fuel pipe diameter of (a) 9 mm and (b) 13 mm for fuel mass flow rate of 0.003 kg/s. The jet is penetrated deeper for 9 mm fuel pipe diameter due to higher momentum. Figure 18. Comparison of turbulent viscosity ratio at the symmetry plane with the fuel pipe diameter of (a) 9 mm and (b) 13 mm for fuel mass flow rate of 0.003 kg/s. Higher turbulence is observed for the 9 mm fuel pipe diameter due to higher velocity magnitude. Figure 19 shows the mass fraction distribution of CH 4 at the symmetry plane with the two fuel 428 pipe diameters. It can be seen from the comparison that the CH 4 is exhausted by the reaction earlier in 429 case of 9 mm diameter. Whereas, the CH 4 is available till higher depth in case of 13 mm diameter. This 430 aligns with our expectation that due to higher turbulence CH 4 reacts readily with oxidizer in case of 9 431 mm diameter. This further results in higher temperatures in the furnace for 9 mm fuel pipe diameter. distribution at YZ plane passing through burner 1 for 9 mm and 13 mm diameter fuel pipe, respectively.

435
The maximum temperature with 9 mm fuel diameter pipe is high for both burners 1 and 2 as compared 436 to 13 mm. This is attributed to the higher rate of reaction occurring in accumulated locally in the 437 furnace. Figure 19. Comparison of mass fraction of CH 4 at the symmetry plane with the fuel pipe diameter of (a) 9 mm and (b) 13 mm for fuel mass flow rate of 0.003 kg/s. With increased diameter, turbulence decreases due to which combustion process is slower resulting into longer unreacted CH 4 jet. Figure 20. Comparison of temperature [ o C] at the YZ plane with the fuel pipe diameter of (a) 9 mm cutting through burner 1 (b) 13 mm cutting through burner 1 (c) 9 mm cutting through burner 2 and (d) 13 mm cutting through burner 2 for fuel mass flow rate of 0.003 kg/s. Due to increased turbulence with 9 mm fuel pipe diameter, combustion is faster and therefore, temperature is higher compared to 13 mm fuel pipe diameter.

439
In the above sections, the effect of varying operating conditions on the mass fraction, temperature  Before proceeding into discussion on the calculations of NOx, a summary on the impact of fuel 445 pipe diameter and fuel injection velocity is discussed. Figure 21 shows the temperature at the XY 446 symmetry plane for the studied variations. It can be observed that by increasing the fuel pipe diameter, It can be observed that for a particular symbol, for example ' ', the red plots are always higher than 471 the blue plots. This suggest that with increase in the fuel pipe diameter, the NOx across the furnace 472 decreases. Furthermore, for a particular color, for example red, the symbol 'o' is always lower than ' ' 473 which is further lower than ' '. This implies that for a given fuel pipe diameter, the NOx decreases 474 with decrease in the fuel pipe diameter. NOx calculated from the model are over predicted as compared to measured values. With increased fuel pipe diameter, NOx decreases for a certain mass flow rate of fuel. While, NOx increases with increase mass flow rate of fuel for a certain fuel pipe diameter.

476
In this paper, a three dimensional model of the heating section of the anode baking furnace is 477 examined to analyse the thermal NOx formation. The COMSOL multi-physics software is used for the 478 modeling of non-premixed turbulent combustion along with the conjugate heat transfer through lining.

479
In the first part of the paper, the calibration of the model is carried out by comparing temperature and 480 species mass fraction with the measured data from Aluchemie. The test model with tuned diffusion 481 parameter that resembles the existing design and operating conditions provide good comparison 482 with the measured data for the temperature. The O 2 mass fraction is underestimated with the model.

483
However, the trend of O 2 mass fraction through the furnace is well predicted.