Utilization of basic multi-layer perceptron artificial neural networks to resolve turbulent fine str

Abstract

Keywords: Computational Fluid Dynamics; Artificial Neural Networks; Machine Learning; Chemical Kinetics; Combustion; Eddy Dissipation Concept

Introduction

The central topic of the current research will be the application of machine learning (artificial neural networks) and combustion fields to develop a modelling approach capable of resolving chemical kinetics and turbulencechemistry interaction for gas-phase combustion without the extensive computer resources required by the current techniques used in commercial computational fluid dynamics codes. The required computer resources are reduced by using artificial neural networks (ANNs), rather than traditional initial value problem ordinary differential equation solvers or look-up tabulation methods, to predict the fine structures' temporal species mass fraction changes. The approach is developed for future implementation in combustion modelling of industrial-scale biomass, coal, oil, and gas-fired boilers using large multi-step chemical mechanisms. The current convention in gas phase combustion modelling is the utilization of the EDM (eddy dissipation model) and finite rate chemistry (FR) combustion model hybrid with two-step global reaction mechanisms, according to [1].This model works under the assumption that the combustion rate (

rate limited (net reaction rate of the chemical reactions, usually determined as a function of the Arrhenius rate), or mixing-rate limited (rate of intermolecular mixing of oxidiser with fuel and/or radical species [2]). The EDM-FR model's low computational cost has made it a popular modelling technique for industrial combustion simulations, [3]; however, the model's empirical constants need to be adjusted depending on the application and are not universally valid [4]. The other major shortcoming of this model, is that the model will yield incorrect results when using more than two reactions, due to the infinitely fast chemistry or mixing limit assumptions. The other combustion models used in industry, such as the conserved scalar PDF approach (non-premixed combustion), steady diffusion flamelet, unsteady diffusion flamelet and flamelet generated manifolds, are not always viable options due to the fact that these models (with the exception of the unsteady diffusion flamelet) only moderately account for chemical non-equilibrium. This paper was recommended for publication in revised form by Regional Editor Ahmet Selim Dalkılıç 1 John Thompson Boiler and Environmental Solutions, Cape Town, South-Africa 2 department Of Mechanical And Mechatronic Engineering, Stellenbosch University, Stellenbosch, South-Africa *E-mail address: rlaubscher@sun.ac.za Manuscript Received 17 July 2017, Accepted 19 September 2017

Journal of Thermal Engineering, Research Article, Vol. 4, No. 2, Special Issue 7, pp. 1828-1846, February, 2018 The eddy dissipation concept (EDC) combustion model is an extension of the EDM, but is able to incorporate detailed chemical mechanisms. The EDC is based on the turbulent energy cascade (larger eddies break/cascade into smaller eddies due to viscous and inertial forces), and assumes that the species are molecularly mixed in the volume that the smallest eddies occupy (fine structures). The position and volume fraction that these fine structure regions occupy is determined from the turbulent field parameters (𝑘 and 𝜀). The temporal evolution of the species in the fine structures is then calculated using ideal reactor theory (well-stirred or plug-flow reactors) [5]. The set of differential equations used to resolve the temporal evolution of the species within the fine structure reactors are said to be numerically stiff. The wide variation in chemical time scales has a severe impact on the numerical solution of the set of differential equations that govern the reaction system. The two techniques currently used to solve these stiff initial value problem equations (in ANSYS 17.0) are the direct integration (DI) [6] and insitu adaptive tabulation (ISAT) [7] methods. Both these numerical techniques have drawback namely. The DI method has large processor requirements (and low primary memory requirements) and the ISAT technique have large required primary memory requirements (and reduced processor time requirements) [8]. An industrial boiler model (± 10 × 106 million cells) using the EDM model (two-step mechanism with 4 species) solves at a speed of roughly 4.2 iterations per minute and takes about 0.83 days to converge, whereas the EDC (same two-step mechanism used in the EDM) model using ISAT solves at roughly 0.5 iterations per minute, resulting in a solving time of approximately seven days before convergence is reached. Seven days is too long a period for industrial design CFD simulations where multiple design permutations must be analysed, various load cases solved and, in some cases, multiple fuels' performance determined. The goal of the present study is thus to develop a chemistry integrator for a steady-state EDC model capable of solving at similar speed as the EDM model, making it a viable industrial design and analysis tool. As mentioned, ANNs is used to predict the temporal species changes in the turbulent fine structures due to chemical reactions in the current research. The ANN interpolates the compositional space by taking the cell's initial mass fractions (𝑌𝑘 ), temperature (𝑇) and fine structure time scale (𝜏 ∗ , which is the time over which the species react in the fine structure regions of a cell) and passing it through the network by forward propagation to produce the predicted mass fractions for the given initial states and fine structure time scale. The CPU load running the ANN chemistry integrator is reduced over the DI method by reducing the amount of required operations that are performed. The ANN integrator replaces the time-stepping procedure of the DI method by a single forward propagation procedure where the initial states and time scale of the fine structure reactor is directly mapped to the species mass fractions after reacting over the period 𝜏 ∗ . The RAM requirement is also reduced, seeing as the ANN does not have to save the large solution tables but only the network weights, which is much smaller than the memory required to store the look-up tables. The objectives of the present research were to: 1. 2. 3. 4.

Generate training data on which to train the neural networks Implement ANN chemistry integrator in ANSYS Fluent and bypass traditional solvers Generate results and compare ANN chemistry integrator results to the results generated using the traditional integration methods and the experimental measurements Show that the ANN chemistry integrator technique for the EDC combustion model could solve at the same speed as the EDM model, thus making it an industrially viable CFD design tool for combustion

The Sandia Flame D (fully turbulent diffusion methane/air flame 𝑅𝑒 = 22400) experimental results are used as validation of the model along with the solutions generated by utilising the standard integration (DI and ISAT) methods to solve the Flame D setup. A quasi-global five-step chemical mechanism, namely the Westbrook and Dryer 1 mechanism, was used to simplify the implementation of the ANN chemistry integrator [9]. Figure 1 shows a simple sketch of the experimental setup. ANNs have previously been implemented in combustion simulations by various researchers ( [10], [8] and [11]) but not in the manner the present research applies it. The work of Kempf et al. [12] and Fard et al. [13] was focused on using ANNs for storage and interpolation of the steady laminar flamelet libraries in CFD simulations of jet flames. They found good agreement between the results generated using the ANNs and the results from the traditional chemistry integration approaches and experimental data, while observing reduced computer-memory requirements.

Journal of Thermal Engineering, Research Article, Vol. 4, No. 2, Special Issue 7, pp. 1828-1846, February, 2018 Another noticeable research contribution is by Christo et al. [14] which investigated using ANNs with a modelled velocity-scalar joint probability density function (PDF) transport equation for 𝐻2 /𝐶𝑂2 turbulent jet flames.

Figure 1. Sandia flame sketch Theory and calculations Since the model is a combination of the ANN chemistry integrator and the EDC model, a brief description of both models will first be discussed along with a short discussion of the governing equations used to model the experimental flame. Next, the methodology used to generate the training and testing datasets for the neural networks' learning phase will be given. The learning and testing phase results will be stated and discussed in order to demonstrate the accuracy of the ANNs in capturing the incremental species changes for different initial states and fine structure time scales. In the next section; results of the Sandia Flame D CFD model using the ANN chemistry integrator will be presented. Governing equations There are three main differences between conservation equations for reacting flows and for non-reacting flows [15]: 1.

A reacting gas is a non-isothermal mixture of multiple species (hydrocarbons, oxidising species and moisture) where each species has it is own continuity equation that needs to be solved. The heat capacities of the reacting gas change significantly through the computational domain, thus making the simulation much more complex and unstable Species react chemically, and the combustion/reacting rate is a function of the turbulent mixing and chemical kinetics involved at that specific spatial and temporal position in the domain Since the gas is a mixture, the transport coefficients such as heat diffusivity, species diffusion, viscosity and so forth needs special attention

This subsection will define the basic conservation equations used in modelling reactive flow. Using the Favre averaged decompositions and applying them to the steady-state mass, momentum, species and energy transport equations yields the following time-averaged conservation equations: Mass,

Journal of Thermal Engineering, Research Article, Vol. 4, No. 2, Special Issue 7, pp. 1828-1846, February, 2018 Momentum,

Where 𝑢𝑖 is the velocity component in tensor notation, 𝜌̅ is the cell Reynolds average density, 𝑃̅ is the cell Reynolds average pressure tensor, 𝜇 is the viscosity, 𝑘 is the turbulent kinetic energy, 𝑌̃𝑘 is the Favre averaged mass fraction of species k, 𝑌̅𝑘 is the cell Reynolds averaged species mass fraction, 𝑆𝑐 is the Schimdt number and ̅̅̅̅ 𝜔̇ 𝑘 is the net reaction rate of species k. 𝜇𝑡 is the turbulent viscosity and is calculated using

The kinetic and dissipation (𝑘 and 𝜀 respectively) is calculated using a turbulence model. In the current work the realizable 𝑘 − 𝜀, 𝑘 − 𝜔 SST and RSM models were used. EDC gas combustion model The EDC combustion model uses the turbulent energy cascade to model the eddies' distribution for reducing energy levels all the way to a level where the eddies' length and time scales are comparable to the Kolmogorov scales. At these fine scales the assumption is made that the fuel and oxidiser eddies are molecularly mixed and react if the temperature and species mass fraction permit it [3]. Thus, the reactions are restricted to these fine scales also named the fine structure regions. The fluid volume is therefore divided into the reacting fine structure regions with a volume fraction;

where 𝐶𝛾 is the fine structure length fraction and the model constant 𝐶𝛾 = (3𝐶𝐷2 /4𝐶𝐷1 )0.25 = 2.1377. The volume fraction of the fluid, which is inert, is calculated as;

The transfer rate of the surrounding fluid into the fine structure regions can be defined as the inverse of the fine structure time scale 𝜏 ∗ , and in turn the fine structure time scale can be defined as;

Journal of Thermal Engineering, Research Article, Vol. 4, No. 2, Special Issue 7, pp. 1828-1846, February, 2018 𝐶

where in Equation 8 the model constant is defined as 𝐶𝜏 = ( 𝐷2 )0.5 = 0.4082. In Equations 8 and 6: 𝑘 and 𝜀 are 3

the turbulent kinetic energy and the turbulent dissipation rate calculated from the turbulence model and 𝜈 is the kinematic viscosity of the fluid. Knowing the fine structure regions' time scale and length/volume fraction, the net reaction rate of species k is formulated as; ˙

where 𝜌 is the fluid volume's averaged density (in Fluent 17.0 this is just the cell's density), 𝑌𝑘 is the mass-averaged mass fraction of species k in the fluid volume (where is Fluent 17.0 this is the mass fraction of the cell at the beginning of the each iteration) and 𝑌𝑘∗ is species k mass fraction in the fine structure regions. To calculate the fine structures' mass fractions for each iteration, these regions are assumed to be tiny ideal plug-flow reactors with initial mass fractions equal to the fluid volume's mass-averaged mass fractions (𝑌𝑘 ) at the beginning of each iteration, and residence time is equated to the fine structures' time scales (𝜏 ∗ ). Therefore, the plug-flow reactor's species transport and energy equations are integrated over the residence time in the fine structures. The species transport equation for species k in the plug-flow reactor is; ˙

In Equations 11 and 10: 𝜌 is the fine structure plug-flow reactor's density based on its mass fractions and temperature, 𝑐𝑃 is the specific heat of the gas mixture, 𝑀𝑘 is the molecular weight of the species k, ℎ𝑘 is the partial ˙

molar enthalpy of the species k, 𝑇 ∗ is the fine structure plug-flow reactor temperature and 𝜔𝑘,𝐹𝑅 is the net reaction rate of species k calculated using finite rate kinetics. To calculate the net reaction rate of a species using the finite rate kinetics approach we consider a chemical system of N species reacting through Z reactions; 𝑁

′ ″ ∑ 𝜈𝑘𝑗 𝑋𝑘 ↔ ∑ 𝜈𝑘𝑗 𝑋𝑘 𝑗 = 1,2,3, . . . , 𝑍 𝑘=1 𝑘=1 ′ ″ where X is the symbol for the chemical species k, 𝜈𝑘𝑗 and 𝜈𝑘𝑗 are the molar stoichiometric coefficients of species k as reactant or product respectively in reaction j. The chemical source term (net rate of production) for species k, ˙

𝜔𝑘,𝐹𝑅 is the sum of all rates for the specific species produced by all Z reactions: 𝑍 ˙ 𝜔𝑘,𝐹𝑅 = 𝑀𝑘 ∑(𝜈𝑗″ − 𝜈𝑗′ )𝑞𝑗 𝑗=1

where 𝑞𝑗 [𝑘𝑚𝑜𝑙/𝑚3 𝑠]is the rate of progress for reaction j and is defined by:

Journal of Thermal Engineering, Research Article, Vol. 4, No. 2, Special Issue 7, pp. 1828-1846, February, 2018 𝑁

In Equation 14 [𝑋𝑘 ] is the molar concentration of species k ([𝑋𝑘 ] =

reverse specific reaction rate constants [16]. The specific reaction rate constants for the Z reactions are modelled using the Arrhenius law. The Swedish chemist/physicist Svante Arrhenius postulated that only those molecules that possess energy greater than a certain amount of 𝐸𝑎 activation energy will react [17]. Using Arrhenius' Boltzman factor, the forward and reverse specific reaction rate constants for the j reaction can be calculated using:

where in Equation 15 𝐴𝑗 [1/s] is the pre-exponential factor, which is assumed to include the effect of the collision terms and the steric factor associated with orientation of the colliding molecules [16]; 𝛽𝑗 is the temperature exponent for the j reaction; 𝑅 is the universal gas constant; and 𝑇 is the temperature of the system. Therefore a rate of progress variable will exist for each reaction in the selected reaction mechanism. The Westbrook and Dryer mechanism (WD1) was selected, which is comprised of a single fuel oxidation reaction that forms 𝐶𝑂 and𝐻2 𝑂. Next, the mechanism has a reversible 𝐶𝑂 − 𝐶𝑂2 reaction. The rate constants originated from studies of 𝐶𝐻4 and 𝐶𝑂 oxidation under fuel lean conditions in a turbulent flow reactor, and they are now widely used for hydrocarbon combustion modelling as discussed by [18], [19]. The reverse reaction step for the 𝐶𝑂2 decomposition is used in order to reproduce proper heat of reaction and pressure dependence of the 𝐶𝑂/𝐶𝑂2 equilibrium. The mechanism, in addition, includes an 𝐻2 oxidation reaction.

Table 1 contains the kinetic data for the reaction mechanism in equation 16. The traditional solution strategy to solve the chemical reaction and energy equations for the fine structure plug-flow reactors uses backward differencing formulation and variable coefficient ordinary differential equation (VODE) discretisation and solver respectively - the reader is encouraged to read [20] and [21] for further details on the numerical solution method. The present work therefore sets out to replace the integration procedures of the fine structure plug-flow reactors [numerical solution of Equations 10 and 11] with the predictions of the neural networks. Each cell's fine structure length fraction and time scale is calculated in the CFD simulation and used along with each cell's mass fractions and temperatures as input to the neural network integrator model. A neural network is configured for each of the different fine structure time scale intervals found in typical combustion simulations ( Δ𝑡 = 𝜏 ∗ = 10−7 𝑠, 10−6 𝑠, 10−5 𝑠, 10−4 𝑠). Therefore, based on the given fine structure time scale value, the corresponding neural network is used to predict the species mass fractions after reacting over 𝜏 ∗ . In the next section the working of the ANN chemistry integrator will be presented with supporting equations. Table 1. WD1 reaction mechanism kinetic data No. 1 2 3 4 5

Reaction order [𝐶𝐻4 ]0.7 [𝑂2 ]0.8 [𝑂2 ]0.25 [𝐻2 𝑂]0.5 [𝐶𝑂] [𝐶𝑂2 ] [𝐻2 ][𝐻2 ]0.5 [𝐻2 𝑂]

Journal of Thermal Engineering, Research Article, Vol. 4, No. 2, Special Issue 7, pp. 1828-1846, February, 2018 ANN chemistry integrator model Neural networks is a technique utilised in computer- and cognitive science and machine learning to approximate complex functions with a high dimensional feature (input variables) and response space (output variables). A neural network is comprised of multiple neurons which act as individual processing units that weight each of the incoming signals and sum it before passing the cumulative signal to a threshold or activation function. Once the signal exits the activation function block, it is passed to all the neurons in the next layer of the network. This process is called forward propagation. In multi-layer networks, the inputs propagate through the network. The (𝑙)

neural network is parametrised by 𝑤𝑖𝑗 , which is the weight of the network, with; (17)

where d is the amount of neurons in layer 𝑙, 𝐿 is the amount of layers and 𝑤𝑖𝑗 is the weighting variables for the 𝑙 (𝑙)

where 𝑥𝑖 is the input signal from the i neuron in the preceding layer (𝑙 − 1). Next the cumulative signal is passed into the activation function for the current layer. It was shown in [22] that the best-performing activation function configuration to predict the incremental species changes due to chemical reactions for the Westbrook and Dryer 1 mechanism is to use a hyperbolic-tangent activation function for the input and hidden layers and a linear activation function on the output layer. (𝑙)

The variable 𝜃(𝑠𝑗 ) is the output signal from the 𝑗𝑡ℎ neuron in the 𝑙 𝑡ℎ layer. From here on the process above is repeated for each neuron in the 𝑙 𝑡ℎ layer before moving on to the next layer or output layer, depending on the current position in the network. Once the input signal has propagated through the entire network, the output (𝐿)

layer signal is generated for each neuron in the final layer, 𝑥𝑖 . A performance function is used to gauge how well the neural network predicts the outputs from the input dataset. The error for the observation is calculated as; (19)

Equation 19 shows that the error is a function of the weights, therefore the weights need to optimised to minimise the error, in turn maximising Equation 19. This process is called the learning procedure of the neural network and is performed by implementing the back-propagation algorithm. Weights are adjusted using the stochastic gradient ascent algorithm with the momentum modification, as seen in the relation below. (𝑙),𝑘+1

is the adjusted weighting variable for the next iteration and 𝑤𝑖,𝑗

value. 𝜂 and 𝛾 are the learning rate and momentum relaxation parameters. In Equation 20; 𝛿𝑗

(𝐿) in which layer the algorithm is, for the output layer 𝛿𝑗 is;

Journal of Thermal Engineering, Research Article, Vol. 4, No. 2, Special Issue 7, pp. 1828-1846, February, 2018 For any other layer in the neural network; 𝑑(𝑙)

Weights are continuously adjusted by performing multiple stochastic gradient ascent iterations on randomly picked input and output observation pairs, until the error value reduces below a desired value. The input and output datasets used to optimise the weighting variables is called the training dataset. Before the training dataset generation is discussed, we will first examine the requirements of the network. The inputs to the network should define the reactor initial conditions; therefore the species mass fractions and the temperature are required. The input vector to the input layer for observation i will take the following form: Input𝑖 = [𝑇, 𝑌𝐻2 , 𝑌𝑂2 , 𝑌𝐻2𝑂 , 𝑌𝐶𝐻4 , 𝑌𝐶𝑂 , 𝑌𝐶𝑂2 , 𝑌𝑁2 ]𝑇

The response vector for the i observation should be the species mass fractions after reacting over the specified reaction time step (where each reaction time step has a neural network); therefore the vector will be: ∗ Output𝑖 = [𝑌𝐻∗2 , 𝑌𝑂∗2 , 𝑌𝐻∗2𝑂 , 𝑌𝐶𝐻 , 𝑌 ∗ , 𝑌 ∗ , 𝑌 ∗ ]𝑇 4 𝐶𝑂 𝐶𝑂2 𝑁2

where 𝑌 ∗ is the fine structure regions mass fractions after reacting over time 𝜏 ∗ . The input layer of the neural network will therefore be fixed to 8 neurons and the output layer fixed to 7 neurons. The output variables do not include the temperature as for the input layer, because as seen in equation 9, the reaction rate equation of the EDC model only requires the fine scales species mass fractions for it to be evaluated. Laubscher et al. [22] performed a study to find the best-performing network configuration (amount of hidden layers and amount of neurons per hidden layer) for the current problem at hand. It was found that for the Westbrook and Dryer 1 mechanism a single hidden layer with seven neurons yielded the best results using the stochastic gradient ascent with momentum modification algorithm. Figure 2 shows a representation of the neural network used to capture the incremental species changes of the methane/air combustion mechanism.

Figure 2. Artificial neural network chemistry integration for the Westbrook and Dryer 1 chemical mechanism

Journal of Thermal Engineering, Research Article, Vol. 4, No. 2, Special Issue 7, pp. 1828-1846, February, 2018 During the back-propagation algorithm, the routine uses a training dataset (input values and corresponding known output values) to quantify the relationship between the inputs and outputs and change the weights accordingly. The EDC model approximates the reaction zones as tiny uniform single plug-flow reactors each with a residence time which calculated from the turbulence field per global numerical iteration. The training data can therefore be generated by solving only a single plug-flow reactor model with varying initial species compositions and residence time. The training data thus is generated, using a separate computer program of a single plug-flow reactor in Python 3.0.0 (using libraries from SciPy and Cantera 2.2.1 chemical kinetics and properties database), by creating random samples of initial species mass fractions and temperatures and integrating the plug-flow reactor's species and energy equations to steady state using the traditional BDF and the VODE solver. The plug-flow reactor's species mass fractions and temperature are recorded into the training dataset for each time step, where input-output sets are created by taking a time step's corresponding species mass fractions and temperature as the input vector and the next time step's corresponding mass fractions as the output vector. The neural networks implemented into the CFD model was trained using 20 000 random initial samples along with all the calculated time step mass fractions/temperatures and, as mentioned, a neural network was trained for each of the following plug-flow reactor residence times (the same as the fine structure time scale,𝜏 ∗ ): Δ𝑡 = 𝜏 ∗ = 10−7 , 10−6 , 10−5 , 10−4 . In Figure 3 the in-sample error (which is the squared difference between the actual output value of the observation and the predicted value by the neural network) for the prediction of each species in the chemical mechanism is plotted as a function of the plug-flow reactor residence time. The graph shows that as the residence time increases, the in-sample error generally also increases, indicating that the relationship between the input vector and the output vector becomes more complex. The increased complexity in the relationship between the input and output vectors is because the changes in species mass fractions are larger as the residence time of the reactor increases, and thus the relationship is increasingly non-linear, but this observation is more pronounced for certain species such as 𝐶𝑂2 , 𝐶𝑂 and 𝐻2 .

Figure 3. In-sample error as a function of reaction time step for single hidden layer (7 neurons) neural network To ensure the neural networks trained on the plug-flow reactor VODE solver's results can be used with confidence in the EDC turbulent combustion model of the Sandia flame, the estimated out-of-sample error should be minimised. Where the out-of-sample error is the error of the neural network's predictions for observations it has not trained on, the generalisation error is the difference between the in-sample and the out-of-sample errors. Therefore, if the generalisation error is minimised, the out-of-sample error is said to track or follow the in-sample error. This means the magnitude of the out-of-sample error will be very similar to the in-sample error calculated during training of the networks. Figures 4 and 5 show the species averaged squared in-sample and out-of-sample prediction errors for increasing training dataset sizes for two residence time step durations.

Journal of Thermal Engineering, Research Article, Vol. 4, No. 2, Special Issue 7, pp. 1828-1846, February, 2018

Figure 4. In-sample and out-of-sample error for 0.5E-6 network model of increasing dataset size

Figure 5. In-sample and out-of-sample error for 0.5E-3 network model of increasing dataset size Numerical simulation of Sandia Flame D experimental setup The flame selected for validation of the ANN chemistry integrator EDC model was a fully turbulent methane/air jet diffusion flame that was measured by [9]. The burner is located in a wind tunnel, where the walls of the wind tunnel are far enough from the flame so that the wall effects are negligible. The burner consists of a small round jet feeding the fuel stream, an annulus around the main jet feeding the pilot gas and air being fed outside the annulus. The geometry was modelled as a 2D-axisymmetric domain to simplify the model geometry and reduce mesh size. Figure 6 shows the computational domain. The main jet diameter d had an inner diameter of 7.2 mm, the pilot annulus inner diameter is 7.7 mm with a wall thickness of 0.25 mm, the pilot annulus outer diameter is 18.2 mm and the burner outer wall diameter is 18.9 mm with a wall thickness of 0.35 mm.

Figure 6. 2D-axisymmetric computational domain used to model piloted jet flame

Journal of Thermal Engineering, Research Article, Vol. 4, No. 2, Special Issue 7, pp. 1828-1846, February, 2018 The boundary conditions applied to the geometry are contained in Table 2. The co-flow inlet boundary condition is the surrounding air that flows parallel to the flame and is comprised of 76.7% nitrogen and 23.3% oxygen (mass based). The ambient air is at 300 K and the inlet velocity set to 0.9 m/s (measured during experiment) which is the wind tunnel exit velocity. The fuel inlet is also modelled as a velocity inlet and it has a composition of 19.664% oxygen, 64.73% nitrogen and 15.607% methane and is assumed to be at ambient temperature before it is injected into the combustion zone at a velocity of 49.6 m/s. A pilot flame is used to ignite the diffusion flame and was modelled as a velocity inlet with a composition of 5.6% oxygen, 74.2% nitrogen, 11% carbon-dioxide and 9.2% water vapour. The gas inlet temperature was measured during the experiment and the numerical model's inlet value was set to the averaged measured value of 1880 K. The pilot gas enters the computational domain through the concentric annulus surrounding the fuel inlet nozzle. The jet and pilot nozzle walls is assumed to have a zero heat flux (𝑞 = 𝑘

= 0) due to the fact that the respective thickness of the two nozzle is very small. Table 2. Boundary conditions

Assigned values Velocity magnitude = 0.9m/s Temperature = 300 K 𝑌𝑂2 = 0.233, 𝑌𝑁2 = 0.767

Velocity magnitude = 49.6m/s Temperature = 300 K 𝑌𝑂2 = 0.19664, 𝑌𝑁2 = 0.6473, 𝑌𝐶𝐻4 = 0.15607

Velocity magnitude = 11.4m/s Temperature = 1880 K 𝑌𝑂2 = 0.056, 𝑌𝑁2 = 0.742, 𝑌𝐶𝑂2 = 0.11, 𝑌𝐻2𝑂 = 0.092

Heat flux = 0 W/m2 Heat flux = 0 W/m2 Backflow temperature = 300 K Backflow 𝑌𝑂2 = 0.233

The jet and pilot nozzle walls are assumed to be perfectly insulated and the heat flux set to zero. The axis boundary condition shown in Figure 6 corresponds to the center line of the jet nozzle (and therefore the computational domain). Seeing as the flame is located inside a large controlled area where the surrounding atmosphere is stationary, the computational domain's outer boundary can be selected as a symmetry-type boundary condition. The symmetry condition has zero normal velocity and zero normal gradients for all fluid variables. The mixture species and reaction mechanism used are shown in Table 1, the density of the mixture was resolved using the ideal gas law, and the thermal conductivity was calculated using the polynomials specified in the VDI Heat Atlas [23]. The viscosity was calculated using kinetic theory of gasses. The individual species' specific heat, molecular weight, standard state enthalpy and standard state entropy were imported from the GRI 3.0 mechanism's thermodynamic properties database, which uses the NASA polynomials. The pressure-based solver was utilised and the pressure-velocity coupling was resolved using the SIMPLE algorithm. The computational domain was discretised as follows: gradients by least squared cell based, pressure using the standard scheme and density, momentum, 𝑘, 𝜀, species and energy were all discretised using the second order upwind technique. A mesh independency study was performed, where the mesh was refined while temperature and species mass fraction averages of the symmetry axis was compared after each refinement stage. The results showed that adequate mesh independency was found for a mesh size of 2352. The largest variation was seen with the averaged mass fraction value for hydrogen. The variation is acceptable seeing as the hydrogen mass fraction varies between 0% and 1% of the mixture.

Journal of Thermal Engineering, Research Article, Vol. 4, No. 2, Special Issue 7, pp. 1828-1846, February, 2018 Table 3. Property variance (%) for sequentially refined mesh Mesh size T 𝑌𝐶𝐻4 𝑌𝐶𝑂 𝑌𝐶𝑂2 𝑌𝐻2 𝑌𝐻2𝑂 𝑌𝑂2

At first the turbulence was resolved using the realizable 𝑘 − 𝜀 turbulence model with its standard model constants, as used in ANSYS Fluent 17.0. (Table 3). Investigation of the results for the CFD model using the traditional DI solver for the fine structure chemical kinetics and the realizable 𝑘 − 𝜀 turbulence model showed that the reaction rate of the 𝐶𝐻4 and consequently the 𝐶𝑂 was over-predicted, causing a temperature peak to occur closer to the burner than what the experimental results showed due to the fact that the 𝐶𝑂 and 𝐶𝐻4 is consumed faster. To determine if the results were turbulence-model dependent; the CFD model of the Sandia flame was solved using the 𝑘 − 𝜔 SST and quadratic pressure-strain Reynolds-Stress model (RSM) to resolve the turbulence field. Figures 7 and 8 show that the realizable 𝑘 − 𝜀 and 𝑘 − 𝜔 SST models predict similar results on the centre line of the flame, whereas the RSM model's 𝐶𝑂 and temperature peaks align closer to the peaks of the experimental results. The over-prediction of the species' reaction rates by the latter two models is due to the over-prediction of the turbulence around the burner outlet, which in turn increases the chemical reaction rates calculated by the EDC model. This is due to the inability of the two equation-turbulence models to accurately capture the anisotropic region (close to the nozzle) which is a result of the high mean shear rates between the fast moving jet and the slow moving co-flow and pilot gases. The quadratic pressure-strain RSM turbulence model was selected for the validation of the ANN chemistry integrator model. The RSM model predicts the anisotropic region's turbulence field more accurately by using a higher order approximation of the turbulent Reynolds stresses. Seeing as it does not suffer the same shortcomings of the two equation-turbulence models, it does not over-predict the turbulence production and in turn the mixing and combustion rates. The disparity in magnitudes (𝐶𝑂 mass fraction and the temperature) between the experimental measurements and the CFD model results (heights) is due to the chemical mechanism used; if one were to utilise a more detailed mechanism with more species and reactions, the results and the experimental measurements would align closer to each other as shown by [3], [24], [4] and [25].

Figure 7. CO mass fraction results on center line for various turbulence models

Journal of Thermal Engineering, Research Article, Vol. 4, No. 2, Special Issue 7, pp. 1828-1846, February, 2018

Figure 8. Temperature results on center line for various turbulence models Implementation of ANN integrator The variables passed to the ANN chemistry integrator are the species mass fractions of the cell, the fine structure time scales 𝜏 ∗ and the cell temperature. To output the incremental species changes in the fine scales over 𝜏 ∗ , the implemented code should calculate the mixture fraction of the cell and determine if the cell is reacting based on the limits of the cell's fuel richness or leanness (shown to be 0.1 → 0.6 where the experimental stoichiometric mixture fraction is 0.354). If the cell is reacting, the code should perform a forward propagation routine calling the correct weights and bias matrices based on the cell's fine scale chemical reaction time. Once the fine scales mass fractions are calculated, these values are utilised to determine the reaction rates for all the species transport equations and the energy release rate in the energy equation. The reaction rates and energy release rate should then be hooked into Fluent's conservation equations. The process described here is shown in Figure 9. An alternative approach to model detailed chemistry with the EDC combustion model was proposed and implemented by [26] and [27], who used the EDC as a chemistry solver applied to a developed flow, turbulence and temperature field. This entails solving the combustion with the EDM for a two-step mechanism for the given fuel. Once the simulation is converged, the flow, turbulence and energy equations are disabled. The EDC is then activated and a detailed chemical mechanism is loaded into the simulation and solved; thus the EDC and ANN integrator is used as a chemistry post-processor to predict intermediate and pollutant species. [27] found that the EDC model resulted in slow fuel conversion rates for a 50 kW non-swirling natural gas combustor, but once the EDC calculations were superimposed on the flow-field from the EDM model, superior results were achieved. [26] and [27] showed that reasonable results could be achieved with this approach. This alternative approach was also modelled in the present research and compared to the experimental data, EDC method with traditional chemistry integration and EDC with ANN chemistry integration (complete combustion and flow model). [4] showed that the original Magnussen constants 𝐴𝐸𝐷𝑀 and 𝐵 of 4 and 0.5 did not yield good prediction accuracy for a 𝐶𝑂/𝐻2 /𝑁2 flame, in that it over-predicted the reaction rate of the 𝐶𝑂 oxidation reaction. Scharler et al. proposed that the EDM constants 𝐴𝐸𝐷𝑀 and 𝐵 should be 0.6 and 0.5 respectively. These constants were used in the EDM model for reaction 2 in Table 1. The simulation was also configured only to solve the first two reactions of the WD1 mechanism. The ANN chemistry integrator of the full mechanism was therefore used as a chemistry post-processor on the converged flow, turbulence and energy fields calculated by using the traditional EDM combustion model and a two-step mechanism.

Journal of Thermal Engineering, Research Article, Vol. 4, No. 2, Special Issue 7, pp. 1828-1846, February, 2018

Figure 9. ANN chemistry integrator program implementation into Fluent 17.0

Results And Discussion

This section is comprised of two parts, the first of which is a comparison between the experimental measurements, EDC combustion model using direct integration, EDC combustion model using the ANN chemistry integrator model and the EDC-ANN chemistry post-processor applied on the converged temperature, flow and turbulence fields from the simulation using the EDM combustion model. The second part investigates the computational performance benefits of using the ANN chemistry integrator versus the standard direct integration and ISAT techniques. Figures 10a to 10f show the species mass fractions of the gas mixture components on the center line of the flame for the three models and the experimental measurements (with its variance). The center line corresponds to the axis boundary condition in Figure 6. A comparison was also performed for various radial sections along the center line of the flame, but will not be shown here. Figure 11 shows the temperature distributions of the three models and the experimental measurements along the center line of the flame. For the remainder of this discussion, the models will be referred to as follows: EDC combustion model using the direct integration - model 1, EDC combustion model using ANN chemistry integrator - model 2 and EDC

Journal of Thermal Engineering, Research Article, Vol. 4, No. 2, Special Issue 7, pp. 1828-1846, February, 2018 combustion model with ANN chemistry integrator post-processor applied to converged flow, turbulence and energy fields from an EDM combustion model - model 3. The results show that all three models have comparable accuracy compared to the experimental measurements. From the previously mentioned figures it is seen that model 1 and model 3 produces similar results for all the shown species. Figure 10a shows that models 2 and 3 over-predict the H_2 mass fractions compared to model 1 and the experimental measurements and Figure 10d shows that the ANN integrator models predicts a higher peak in Y_CO closer to the nozzle outlet. The small error in the accuracy of the ANN integrator with regards to the H_2 fractions is due to the prediction error (in-sample and out-of-sample errors) of the ANNs, which is more pronounced for H_2 due to the small quantities involved. In Figure 3 we see that the average in-sample prediction error for H_2 across the reaction time scale range is about 0.005 (square root the graph value), which is a sizable estimated error given that the maximum measured H_2 is only approximately 0.0035. Methods of reducing this error will be discussed in the conclusions and future work section of the article. The chemical energy not yet released in the hydrogen molecules causes model 2 to predict lower temperatures at certain locations along the center line compared to models 1 and 3. This is because the energy equation is solved along with the ANN integrator in model 2, causing the chemical reaction solution to alter the temperature field and produce less water vapour. This in turn causes the CO water-shift reaction to slow down due to the lower temperature and lower amount of water vapour present in the mixture. The slower CO reaction rate therefore causes a peak in CO mass fraction for the two ANN models as seen in Figure 10d. For model 3 the chemical reaction solution does not alter the temperature field seeing as the flow, turbulence and energy equations in the CFD model are frozen but the lower amount of water present causes a slower reaction rate. The difference between the Y_CO prediction between model 2 and model 3, as seen in Figure 10d, is due to the different temperature fields predicted by the EDM and EDC models, which is a result of the method the two modelling approaches use to solve the turbulence-chemistry interaction. The EDC model at very high turbulence rates slightly over predicts the chemical reaction rates which increases temperatures and reduces the Y_CO mass fraction. The over prediction was also found by [27], and could be due to the fact that the combustion is not isolated to the fine structure regions only. To gauge the computational benefits of using the ANN chemistry integrator, multiple simulations were run for increasing mesh sizes. For each mesh size, the flame model was solved using direct integration, ISAT, EDM and ANN integrator techniques. Each simulation was solved until convergence (and with a set number of iterations to ensure the comparison between the various approaches is fair) and the solving time was recorded. The results of these simulations can be found in Figure 12. The results in Figure 12 show that the ANN chemistry integrator definitely offers speed-up versus the conventional direct integration and ISAT integration techniques. This is due to the reduced CPU time required by the ANN integrator. It does not integrate each of the species equations and the energy equation for numerous time steps till the fine structure time scale is reached through time stepping, but rather replaces this procedure with two matrix multiplication operations (forward propagation). As mentioned in the introductory section of this article, one of the goals is to develop an EDC model integrator that can solve at speeds comparable to the EDM method, thus making it an industrially viable design CFD combustion model, with the added benefit of not having tunable variables and accounting for chemical nonequilibrium. The solving times of the EDC-ANN model and the "pure" EDM model for the Westbrook and Dryer 1 mechanism was compared to determine the difference in computational resource requirements for the two turbulence-chemistry interaction models. It must be noted that the EDM model, as discussed, is unable to accurately handle chemical mechanism with more than two steps, but nonetheless the complete mechanism was solved with the EDM model to keep the comparison between the EDC-ANN and EDM models fair (same amount of species transport equations and reactions). The amount of iterations per minute solved using the two models is shown in Figure 13 for increasing mesh sizes.

Journal of Thermal Engineering, Research Article, Vol. 4, No. 2, Special Issue 7, pp. 1828-1846, February, 2018

Figure 10. Species mass fractions at axial locations along center line of flame predicted using RSM turbulence model. Ø –Experimental, Ø -EDC/ANN, Ø -EDC/ANN on EDM two-step temperature and flow field and Ø – EDC with direct integration

Figure 11. Temperature at axial locations along the center line of flame predicted using the RSM turbulence model. Ø –Experimental, Ø -EDC/ANN, Ø -EDC/ANN on EDM two-step temperature and flow field and Ø – EDC with direct integration

Journal of Thermal Engineering, Research Article, Vol. 4, No. 2, Special Issue 7, pp. 1828-1846, February, 2018

Figure 12. Solving time of different combustion models and integration techniques

Figure 13. Iterations per minute for EDM and EDC/ANN models From the results we see that the EDC-ANN and EDM turbulence-chemistry interaction models have almost identical computational resource requirements.

Conclusion

In the current work, the ANN chemistry integrator for the EDC turbulence-chemistry interaction model was developed, trained, tested and implemented into a CFD simulation. The work has proven that an ANN chemistry integrator can be used to resolve the chemical kinetics in the fine structure regions of the turbulence field in a manner that is efficient with regard to both time and computer resources. The training data and the implementation in the CFD model made use of the mixture fraction ranges of the experimental setup. The training data used the compositions of the fuel, pilot and co-flow streams to calculate and use the mixture fraction to determine possible mixture compositions in the domain and train on the randomly selected mixtures and temperatures. The CFD model in turn also used the mixture fractions to determine if the cell is reacting or not. This caused the model developed in the current research to be used only on the Sandia flame D experiment. In future research the mixture fraction limits in the training phase and the simulation phase will be discarded to develop a more generic ANN chemistry integrator approach, thus creating a one-off ANN chemistry integration model for each chemical mechanism selected (GRI 3.0, WD-mult, JL 4-step, etc.). This will entail performing a statistical investigation to determine which size (completely random) training dataset is adequate to capture the actual flame behavior. Neural networks are excellent models at interpolating the compositional space; therefore, the training datasets will never have to be excessively large.

Journal of Thermal Engineering, Research Article, Vol. 4, No. 2, Special Issue 7, pp. 1828-1846, February, 2018 The in-sample and out-of-sample errors of the neural network predictions can be reduced by altering the network structure and having a network with a single species output, rather than having multiple species outputs, as is the case for the current work. Using the single species mass fraction output approach will require a neural network for all of the individual species of the chemical mechanism and then for all the fine structure regions' residence times, greatly increasing the amount of required networks. One future strategy to reduce the amount of neural networks required to train is to use the fine structure time scale as an input variable to the network, thus using the residence time of the plug-flow reactor as a feature. By implementing these changes, one would only have to train a network for each of the species in the mechanism. This will therefore increase the flexibility of the each network due to the increased number of weighting variables per species prediction, which could reduce the estimated errors. Another strategy to reduce the estimated prediction error is to use an advanced optimisation approach with the backward-propagation algorithm. One such technique is the adaptive gradient search algorithm, which reduces the learning parameter dynamically by continuous renormalization of the error gradients, making the procedure more sensitive to small changes in the optimisation space. Further research [22] where the ANN chemistry integrator was trained on the full GRI 3.0 mechanism showed reduced in-sample and out-of-sample prediction errors. This is due to the reduced complexity between incremental steps of species and increased amount of weighting variables that increase the model flexibility. Future research will also entail implementing the ANN chemistry integrator on a full-scale watertube boiler's furnace simulation.

Share and Cite