Numerical Modeling and Multi-Objective Control of Air Pollutant Emissions from Natural Gas Internal Combustion Engines

At this stage, due to the increasingly serious environmental problems and energy consumption, distributed energy systems have become a major hotspot of the current energy development Chuan. As a clean and efficient low-carbon energy source, natural gas can effectively improve the environment, reduce carbon emissions, and optimize the energy structure [1-2]. Recycling through combustion of gas internal combustion units to generate electricity, and waste heat using lithium bromide refrigeration units and other equipment is one of the ways to achieve the gradient utilization of energy, and is also the main way of supplying energy to data centers, hospitals and other users at present [3-5]. And in the actual combustion process, the emission of pollutants is unavoidable. The internal combustion engine is one of the main sources of urban air pollution. During the operation of gas internal combustion units, the main air pollutants produced are NO_x (nitrogen oxides), accompanied by a small amount of SO2 and CH4 and other gases [6-7]. Nowadays, NO_x and SO2 have become important air pollutant control indicators for internal combustion engines, while CO and CH4 are gas components that characterize the combustion state of internal combustion engines [8-9]. These pollutants are very harmful to human health and environment, and their emission control can be reduced by numerical simulation techniques.

In the combustion process of an internal combustion engine, factors such as reaction temperature, concentration of oxidizers and combustibles, reaction time and reactant distribution affect the generation and emission of combustion products [10-12]. Numerical simulation allows for the simulation of combustion and pollutant emission scenarios in situations that are not convenient for humans to physically observe or cannot be naturally generated [13]. It is carried out on a computer, which can provide the change trend of various variables in the combustion process, obtain the generation and emission laws of combustion products, adjust the choice of fuel, and change the reaction temperature, rate, gas transport, particle explosion and other parameters of combustion, so as to maintain the state of combustion reaction, reduce pollutant emissions and improve combustion efficiency, so as to achieve the control of pollutant emissions [14-15]. However, there are specific standards for the emission of air pollutants from internal combustion engines, which refer to the standards that stipulate the emission limits and monitoring requirements for various pollutants [16-17]. Its formulation and implementation are of great significance for reducing pollutant emissions and improving the quality of the air environment. By setting strict emission standards, enterprises can be guided to take effective pollutant management measures to reduce the negative impact on air quality.

In this paper, the control of pollutant emissions is mainly achieved by constructing numerical simulation equations and building and solving a multi-objective control model for air pollutants in natural gas internal combustion engines. The basic equations for the three-dimensional model are constructed for the structure of a natural gas internal combustion engine. According to the principle of fluid dynamics calculation, the hydrodynamic multidimensional simulation model of the operation of natural gas internal combustion engine is established in order to apply fluent software for numerical simulation. Combining the non-dominated genetic algorithm with the multi-objective optimization problem, the TOPSIS method is selected as the decision-making method to solve the multi-objective optimization of the constructed air pollutant emission model of the natural gas internal combustion engine, and the optimal solution close to the Pareto optimal frontier is obtained. Finally, the numerical simulation results of the temperature field, NO_x and combustion characteristics under the baseline operating conditions are shown, and the multi-objective pollutant control arithmetic simulation analysis is carried out.

2

Three-dimensional modeling and numerical simulation of emissions from a natural gas internal combustion engine

In this paper, Fluent software is applied to establish a three-dimensional numerical model of a single natural gas internal combustion engine and to perform three-dimensional numerical simulation analysis of the four-stroke process.

2.1

Basic equations for 3D modeling

For numerical simulation analysis, the following basic equations are first constructed in this paper.

2.1.1

The continuity equation

Equation (1) is the basic form of the mass equation and applies to both compressible and incompressible fluids: (1) $\frac{\partial ρ}{\partial t} + \frac{\partial}{\partial x_{j}} (ρ u_{j}) = S_{m}$

where S_m is for a quality source term or a quality additive term.

2.1.2

Momentum conservation equations

The momentum conservation equation in the i-direction, for example, is as follows: (2) $\frac{\partial (ρ u_{i})}{\partial t} + \frac{\partial}{\partial x_{j}} (ρ u_{i} u_{j}) = - \frac{\partial p}{\partial x_{i}} + \frac{\partial τ_{i j}}{\partial x_{j}} + g_{i} - f_{i}$

where p is the static pressure, ρ is the fluid mixture density, u_i is the i-directional velocity. g_i and f_i are gravity and other external forces, respectively, and t_ij is the viscous stress tensor: (3) $τ_{i j} = 2 μ S_{i j} - \frac{2}{3} μ S_{k k} δ_{i j}$

where S_ij is the fluid strain rate tensor: (4) $S_{i j} = \frac{1}{2} (\frac{\partial u_{i}}{\partial x_{j}} + \frac{\partial u_{j}}{\partial x_{i}})$

S_kk is the fluid dispersion, which characterizes the volumetric expansion and compressibility of the fluid, and μ is the dynamic viscosity coefficient of the fluid. δ_ij is the second-order unit tensor, δ_ij = 1 when i = j, and δ_ij = 0 when i ≠ j.

2.1.3

Energy conservation equation

The solution of the energy equation is also essential when the fluid flow is considered for heat transfer or compression: (5) $\frac{\partial (ρ h_{0} r)}{\partial t} + \frac{\partial}{\partial x_{j}} l (ρ u_{j} h_{0} r) = \frac{\partial}{\partial x_{j}} g l (Γ_{h} \frac{\partial h_{0}}{\partial x_{j}} g r) + S_{k}$

Where h₀ is the stagnation enthalpy, F_h the enthalpy transport coefficient or exchange coefficient, and S_k the internal heat source of the fluid and the portion of the fluid’s mechanical energy that is converted to heat due to viscous action.

2.1.4

Component transport equations

In internal combustion engine simulation, the fuel injection, diffusion, and distribution of the work mass need to be taken into account, and Fluent solves for the mass fraction of localized chemicals i, Y_i by calculating the convection-diffusion equations, as shown in the model’s component transport equation (6): (6) $\frac{\partial}{\partial t} (ρ Y_{i}) + \nabla (ρ \vec{v} Y_{i}) = \nabla \cdot {\vec{J}}_{i} + R_{i} + S_{i}$

where R_i is the rate of mass production of a substance in a chemical reaction and S_i is the rate of production of a discrete phase or a custom source term.

2.2

Main models in internal combustion engine calculations

The solution of the multidimensional model, i.e., it is mainly based on the conservation equations of mass, momentum, energy and chemical components describing the phenomena in the cylinder. However, due to the complexity of the simulation environment, many equations cannot be directly solved, which requires further development of relevant solution models. In addition, due to the fact that the operation of internal combustion engine involves a series of action processes, a number of sub-models are included in the 3D simulation, such as the gas flow model, the natural gas spray mixing model, the chemical reaction model, and the heat transfer model.

2.2.1

Turbulence k − ε Modeling

In the process of solving the basic equations, considering the strong turbulent eddies in the cylinder as well as the role of the near-wall surfaces, the article adopts the RNG model [18] in the form of equations: (7) $\frac{\partial}{\partial t} (ρ k) + \frac{\partial}{\partial x_{i}} (ρ k u_{i}) = \frac{\partial}{\partial x_{j}} g (\frac{μ_{a f f}}{σ_{k}} \frac{\partial k}{\partial x_{j}} g) + G_{K} + G_{b} - ρ ε - Y_{M} + S_{k}$ (8) $\frac{\partial}{\partial t} (ρ ε) + \frac{\partial}{\partial x_{i}} (ρ ε u_{i}) = \frac{\partial}{\partial x_{j}} g (\frac{μ_{q f f}}{σ_{s}} \frac{\partial ε}{\partial x_{j}} g) + C_{1 s} \frac{ε}{k} (G_{k} + C_{3 s} G_{b}) - R_{s} + S_{s}$

where G_k is the turbulent kinetic energy k generation term due to the mean velocity gradient, G_b is the turbulent kinetic energy k generation term due to the buoyancy force, c_3ε is the coefficient of influence of the buoyancy force on the turbulence ε flow, Y_M represents the turbulent dissipation rate generated by expansion pulsations in compressive turbulence, σ_k, σ_ε is the turbulence Prandtl number representing the k, ε respectively, and S_kS_ϵ is the customized source term respectively.

2.2.2

Jet model

A typical distribution for the particle size distribution of ejected droplets is the Rosin-Rammler distribution. This distribution assumes that the diameter d of a droplet is distributed exponentially with respect to the mass fraction Y_d of the droplet at that diameter [19]: (9) $Y_{d} = e^{- {(d / \bar{d})}^{n}}$

where $\bar{d}$ is the average particle diameter and n is the diffusion coefficient, which is determined according to the Fluent settings.

2.2.3

Discrete phase model

When the fuel is injected from the nozzle and enters the cylinder, it forms droplets, and it is necessary to calculate its trajectory and its interaction with the continuous phase during the mixing process with the continuous gas in the cylinder.

The DPM model is often referred to as a discrete phase model.Fluent traces a second discrete phase in the simulated flow field in Lagrangian coordinates after solving the transport equations for the continuous phase.

The choice of phase coupling calculation in Fluent takes into account the influence of the discrete phase on the continuous phase. The presence of the discrete phase affects the flow field of the continuous phase, which in turn affects the distribution of the discrete phase, and the discrete and continuous phases are computed alternately until convergence of the computational results is reached. 1)

Equation of motion for a single particle

The particle trajectory’s are based on the force equilibrium calculations (x direction) obtained in the continuous phase: (10) $\frac{d u_{p}}{d t} = F_{D} (u - u_{p}) + \frac{g_{x} (ρ_{p} - ρ)}{ρ_{p}} + F_{x}$

where F_D(u − u_p) is the trailing force per unit mass of the particle: (11) $F_{D} = \frac{18 μ}{ρ_{p} d_{p}^{2}} \frac{C_{D} Re}{24}$

Where u is the fluid phase velocity, u_p is the particle velocity, ρ is the fluid density, μ is the hydrodynamic viscosity, ρ_p is the particle density, d_p is the particle diameter, and c_D is the trailing coefficient. 2)

Randomized orbital model of particles

Based on the strong turbulent flow in the cylinder, the turbulent diffusion of particles is not negligible. In the random orbit model, along the particle orbit, Fluent in the integration calculation process, the fluid velocity in the particle orbit equation is the instantaneous velocity u + u′(t).

The governing equations for the randomized orbits: (12) $\frac{d u_{k}}{d t} = (\bar{u} + u^{'} - u_{k}) / τ_{r k}$ (13) $\frac{d ν_{k}}{d t} = (\bar{ν} + ν^{'} - ν_{k}) / τ_{r k} + w_{k}^{2} / r_{k} + g$ (14) $\frac{d w_{k}}{d t} = (\bar{w} + w^{'} - w_{k}) / τ_{r k} - w_{k} v_{k} / r_{k}$

where u_k, ν_k, w_k is the instantaneous axial, radial and tangential velocities of the particles, $\bar{u}, \bar{ν}, \bar{w}$ and u′, ν′, w′ are the time-averaged and pulsating velocities of the gas phase in the three directions, respectively.

2.2.4

Combustion rate modeling - EDC modeling

Fluent can simulate problems involving gas-phase reactions, volumetric reactions with chemical reaction kinetics, and reactions with particulate matter.Fluent provides five models for simulating gas-phase reactions, including a general finite-rate model, a non-premixed combustion model, a premixed combustion model, a partially premixed combustion model, and a probability density function (PDF) model. The generalized finite-rate model has wide applicability and can simulate laminar or turbulent flow, premixed or non-premixed combustion systems. The generalized finite rate focuses on the mass fraction of each component by solving the mass transport equation. The transport equation reaction source term R_i is expressed through the Arrhenius chemical reaction rate.

In this paper, the vortex cluster dissipation conceptual model (EDC) will be used. The vortex cluster dissipation conceptual model is actually an extension of the vortex cluster dissipation model to consider the chemical reaction mechanism. It assumes that the reaction occurs on a fine scale with a length scaling calculation: (15) $ξ^{*} = C_{ξ} {(\frac{ν ε}{k^{2}})}^{1 / 4}$

where * represents the fine-scale quantity, C_ζ is the volume constant, and C_ξ = 2.1377, v are the kinematic viscosities.

Time scale calculations of component reactions on a fine grid: (16) $τ^{*} = C_{τ} {(\frac{ν}{ε})}^{1 / 2}$

C₁ is the time scale constant, C₁ = 0.04082, and its reaction rate source term R_i is expressed as: (17) $R_{i} = \frac{ρ {(ξ^{*})}^{2}}{τ^{*} [1 - {(ξ^{*})}^{3}]} (Y_{i}^{*} - Y_{i})$

where $Y_{i}^{*}$ is the mass fraction of component i after reaction time t^*.

Using the EDC model, a detailed kinetic model of the chemical reaction can be imported to simulate the reaction with turbulent flow. The detailed chemical reaction mechanism of natural gas will be mentioned later.

2.2.5

Emission modeling - NO_x

For the combustion of internal combustion engines, for the quality of the environment of the increasingly concerned, exhaust emissions pollution has become the focus of attention, the content of its emissions, has become an important measure of the performance of the internal combustion engine indicators, in particular, the amount of NO_x generation, NO_x, including NO and NO₂, but mainly NO-based.

In internal combustion engines, NO generation mainly comes from the following three situations: (1) thermal NO, which mainly comes from the formation of nitrogen in the air oxidized at high temperatures. (2) Rapid NO, which is generated by the rapid reaction between nitrogen in the air and the hydrocarbon components under the oxygen-rich fuel, and its generation ratio is not high. (3) Fuel NO. In comparison, thermal NO is the main source of NOx generation, which is mainly generated by nitrogen from air in a high temperature environment. Its generation mechanism is the expanded Zeldovich mechanism: (18) $\begin{matrix} N_{2} + O \Leftrightarrow_{k 2}^{k 1} N O + N \\ N + O_{2} \Leftrightarrow_{k 4}^{k 3} N O + O \\ N + O H \Leftrightarrow_{k 6}^{k 5} N O + H \end{matrix}$

3

Multi-objective control methods for air pollutant emissions from natural gas internal combustion engines

3.1

Genetic Algorithms

Genetic algorithm is a global search and optimization method derived from the principle of genetic evolution of organisms in nature, by encoding the genetic material located on the chromosomes, evaluating each chromosome based on pre-set fitness values, and selecting chromosomes by simulating the environmental laws, which ultimately leads to the optimal structure [20].

The genetic algorithm has the following advantages: 1)

Compared with traditional optimization methods, genetic algorithms perform operations such as selection, crossover and mutation on a population rather than on a particular individual, so that the computational process is potentially parallel and multiple individuals can be compared at the same time.

2)

Genetic algorithm is a global optimization algorithm, which can avoid the phenomenon of “dead loop” in the process of a single iterative loop.

3)

Genetic algorithm in the optimization process in accordance with the “survival of the fittest” principle of independent search, the process is simple. It does not need to set and describe all the characteristics of the problem, but only to determine the individual evaluation parameters after selecting the coding method.

4)

Stronger robustness. The objective function does not require whether the search domain is continuous or not, and compared to the traditional genetic algorithm, it can find the optimal value with greater probability.

3.2

Multi-objective optimization problem

A multi-objective optimization problem (MOP) is to find a solution that satisfies all the optimization objectives simultaneously. Since MOP involves multiple optimization objectives and the objectives are independent and contradictory, it is difficult to determine the degree of superiority of the unique solution during the optimization process, so the solution is usually an uncertain set of points, which is the Pareto-optimal solution set. The task of multi-objective optimization is to find out the distribution of the solution set, and then find the appropriate solution according to the conditions and weights.

The multi-objective optimization problem of this paper is described as follows: (19) $\min {f_{1} (x), f_{2} (x), ...., f_{m} (x)}$ (20) $s . t . {\begin{array}{l} l b < x < u b \\ g_{i} (x) \geq 0, i \in I \\ y_{j} (x) \geq 0, j \in J \end{array}$

In the formula,

f(x) - objective function to be optimized, consisting of m sub-objective functions, and

lb - lower bound of constraints for variable x, the

ub - the upper bound of constraints for variable x, the

g_i(x) - the inequality constraint function, the

y_j(x) - the equality constraint function, the

I-Set of inequality constraint indicators, the

J-Set of Equation Constraint Indicators.

A Pareto optimal solution set is defined as a Pareto optimal or non-inferior solution of X is said to be x if, for any solution x, there exists no x ∈ X such that F(x′) = (f₁(x′), f₂(x′), …, f_m(x′)) prevails over F(x) = f_i(x), f₂(x), …, f_m(x)) in the feasible region X. For a multi-objective feasible solution set, there are different relationships between individuals, and for the sake of convenience of description, let P be a set in which there exist n attributes for each element, and the relationship between the elements of P is. 1)

If ∀x, y ∈ P, f_i(x) ≤ f_i(y).i = {1.2…n}, ∃I ∈ {1, 2…n}, ∃I ∈ {1, 2…n}, then x is said to dominate y. In this case, x is non-dominant and y is dominated.

2)

If the domination relation between ∀x.y ∈ P, x and y does not hold, then x and y are said to be either non-dominating or play-irrelevant.

According to the above relationship can be seen, in the feasibility of the region does not have any one solution (hollow point) back to dominate the solid point, the pipe center point that is the inferior solution, the solid point of the set composed of the optimal solution set of the required graph Pareto. In general, multi-objective optimization that seeks the income value or minimum value of the objective function, in order to facilitate the analysis, the objective optimization function can be unified transformation. In this paper, the maximization multi-objective optimization problem is transformed into a minimum value problem: (21) $\max f_{i} (x) = - \min [- f_{i} (x)]$

3.3

Multi-objective genetic algorithm (NSGA-II)

Combining genetic algorithm and multi-objective optimization for solving the optimal solution set of the optimized function, the multi-objective genetic algorithm has three main objective requirements: 1)

To find as many optimal solutions in the solution set as possible.

2)

The solutions in the solution set should be as close as possible to the Pareto optimal frontier.

3)

The solutions found should be uniformly distributed on the global Pareto optimal frontier as much as possible.

There are many kinds of multi-objective genetic algorithms, mainly including random weight method, SPEA, NSGA and NSGA-II. Among them, the non-dominated genetic algorithm (NSGA-II) is widely used because of its faster convergence speed and greatly shortened running time, and the multi-objective optimization function Gamultiobj in Matlab is based on this algorithm, which can display Pareto-optimal solutions in the output. The optimal solutions are displayed in the output result graph, and at the same time, due to the introduction of the crowding operator makes these optimal solutions can be uniformly distributed in the Pareto frontier graph, which meets the solution requirements of multi-objective optimization [21].

Simple genetic algorithm is a genetic process of selection, crossover and mutation of individuals according to a certain ratio under the condition of setting the initial population and the number of termination generations, and updating iteration of the population. The NSGA-II algorithm, on the other hand, has different computational steps. Firstly, according to the dominance relationship between individuals in the population, sorting according to the non-inferiority rank, and then according to the congestion operator on the congestion distance, which can make the optimal solution can maximize the uniform distribution on the Pareto front, select the population and then according to the genetic algorithm’s computational process for the population iteration, and ultimately get the set of solutions to satisfy the conditions.

3.4

Decision-making methods

Decision making i.e. how to select the solution that meets the needs of the real situation from the Pareto optimal solution. In the multi-objective decision-making problem, TOPSIS method is a more effective solution, also known as the superior and inferior solution distance method.TOPSIS method selects the target object by measuring the distance between the target object and the contingent optimal solution and the worst solution, and the closer the distance between the target object and the optimal solution asks, the better, and vice versa, it is non-optimal [22]. The optimal solution of each objective function is the positive ideal solution, and the worst solution in the Pareto optimal set is called the negative ideal solution. In order to achieve the purpose of selecting the most suitable solution in the optimal solution set, the weights need to be selected, and at the same time, the ideal solution and the relative proximity need to be calculated as follows: 1)

Set the original decision matrix. Let the decision problem have n evaluation object and m evaluation indicators, then the jth evaluation indicator of the ird evaluation object is x_ij, and the decision matrix is X: (22) $X_{i j} = [\begin{matrix} x_{11} & x_{12} & ... & x_{1 m} \\ x_{21} & x_{22} & ... & x_{2 m} \\ ... & ... & ... & ... \\ x_{n 1} & x_{n 2} & ... & x_{n m} \end{matrix}]$

2)

Determine the weight vector y = {y_i}, i = 1, 2, …m according to the actual demand.

3)

Multiply the weight vector pointwise with the decision matrix to obtain the weighted decision matrix Z = {z_ij}: (23) $Z_{i j} = y_{i} \cdot X_{i j}$

4)

Determine the positive and negative ideal solutions of the optimization objective: (24) $\begin{matrix} p^{+} = \min (z_{j}) \\ p^{-} = \max (z_{j}) \end{matrix}$

Where j = 1, 2, …m. 5)

Calculate the distance of each scheme to the positive and negative ideal solutions d⁺, d^*: (25) $d^{+} = \sqrt{\sum {(z_{j} - x^{+})}^{2}}$ (26) $d^{-} = \sqrt{\sum {(z_{j} - x^{-})}^{2}}$

6)

Determine the similarity of the programs to the optimal program C_i: (27) $C_{i} = \frac{d^{-}}{(d^{-} + d^{+})}$

4

Results of numerical simulation and multi-objective optimization of air pollutant emissions

4.1

Analysis of simulation results for the baseline operating conditions of natural gas internal combustion engines

In this paper, the ISO condition is used as the benchmark condition, and numerical simulation is mainly carried out for the distribution of flow field, temperature distribution and the distribution characteristics of each component in the combustion chamber.

The simulation results of the combustion chamber outlet in the base case are compared with the actual data of the power plant as shown in Table 1. The simulation errors of the numerical simulation methods in this paper are controlled within 5%, and the simulation results are good enough to meet the requirements of subsequent analysis and optimization.

Table 1.

Comparison of the analog and actual data of the combustion chamber exit

Item	Analog value	Actual value	Error(%)
Combustion chamber outlet temperature(K)	1789.27	1796.57	0.41%
Mean NO_x of combustion chamber exit(ppm)	24.08	24.52	1.79%
Maximum amplitude of pressure pulsation(kPa)	12.79	12.43	2.91%
Average flow rate of exit section(m/s)	115.58	120.61	4.17%
Mean CO of combustion chamber exit(ppm)	5.88	5.95	1.18%

4.1.1

Temperature field analysis

Figure 1 shows the trend of temperature on the center axis. Fuel is ejected from the nozzle, and there is an air passage in the center, so the flame shows a gourd-like distribution. In the middle of the combustion chamber, the flame tightens on the outside, and the bundle of air from the inner air channel narrows until it is absent, constituting a violent combustion in the main combustion section, and the blowing air from the outside is used to prevent excessive temperatures from being generated. Nozzle outlet and between the nozzle part of the reflux, the reflux of high temperature natural gas coil suction so that the high temperature region in the flame tube is mainly concentrated in the middle of the position of the main combustion zone. Therefore, analyzing the steady state temperature field shows that the specific pattern of temperature distribution on the central axis in the time-averaged state is smoother and more time-averaged than the transient results.

By analyzing the temperature at z=800mm, it is possible to get an idea of the transient temperature in the most intense part of the combustion. The maximum temperature is around 1940K and is reached around v=-0.1 and v=0.1, the rest of the section shows a decreasing trend because it is more affected by the cooling air.

4.1.2

NO_x distribution analysis

Figure 2 shows the radial distribution pattern of NO_x in each cross-section, and it can be seen that the amount of NO_x generation decreases gradually along the central axis toward the exit of the combustion chamber. At the tangent line of x=0 at the same axial distance, the NO_x generation shows two hump-like distributions, both peaks are due to the combustion center region at the nozzle exit, and the intermediate trough values are due to the lower temperature of the center nozzle, and subsequently the lower NO_x generation, where the peaks at different axial distances are 2.4 × 10⁻⁵, 1.9 × 10⁻⁵, 1.4 × 10⁻⁵, respectively.

4.1.3

Analysis of combustion characteristics

(a)~(c) in Fig. 3 demonstrate the index parameters of axial distance distribution under different inlet temperature conditions. The overall distribution pattern along the axial direction remains unchanged, mainly due to the change in the area of the low, medium, and high flow velocity regions. This summary plot of the flow velocity pattern shows that as the inlet temperature increases, the flow velocities in all parts of the overall basin increase, and that the flow velocity variations between basins are smaller because of the higher doping for the conditions with high inlet temperatures. At an inlet temperature of 600K to 900K, the flow rate fluctuations gradually flatten out, but always remain at a high level.

As the inlet temperature increases, the extension area of the flame high temperature region is increasing, and the combustion temperature gradually increases in the transition section, the flame cluster formed by the horn-shaped distribution of the blowing air has a wider distribution range than that of the condition with lower inlet temperature. Flame center temperature is higher, the temperature in the combustion chamber is also generally open high, which is also due to the increase in the temperature of the blowing air inlet makes the chamber of the high temperature gas distribution is wider than before.

The graphs of the temperature patterns along the axial distribution for each condition show that the increase in inlet temperature brings about more intense flame fluctuations, which is the main reason for the increase in the amplitude of the pressure pulsations in the combustion chamber. The regularity in the time-averaged results is more distinctive, as the NO volume concentration fraction in the center longitudinal section gradually increases with the increase of inlet temperature and increases from the width of the columnar distribution. In terms of the pattern of the NO volume concentration fraction along the axial distribution, an increase in the inlet temperature brings about fluctuations in the NO production, and ultimately a larger gap in the NO volume fraction at the outlet. At an inlet temperature of 600 K, the volume concentration fraction at the outlet is 1.0 × 10⁻⁵, but when the inlet temperature increases to 900 K, the NO volume concentration fraction at the outlet is 2.2 × 10⁻⁵, which is also closely related to the distribution of the temperature field.

(d) In Fig. 3 shows the analyzed pressure pulsation spectrum at the P2 monitoring point for different inlet temperatures. It can be seen that the amplitude of pressure oscillations in the combustion chamber decreases with the increase of inlet temperature, and the main frequency of pressure oscillations increases with the increase of inlet temperature, the main frequency increases from 983 Hz to 1925 Hz, and the corresponding amplitude decreases from 5131 Pa to 3015 Pa. When the temperature increases, the premixed fuel as well as the blowing air are more mixed, so it makes the combustion more full and more stable. The combustion is more complete and stable. However, it should be noted that in the process of temperature increase, there will be a phenomenon of main frequency increase, so it is also necessary to reasonably control the inlet temperature, so as not to be too high to cause an increase in the speed of sound in the combustion chamber, resulting in the improvement of thermoacoustic instability.

4.2

Analysis of multi-objective pollutant control results

4.2.1

Experimental setup

The multi-objective control method of pollutant emission in this paper is used for experiments. Under normal circumstances, it is impossible to make each sub-objective of multi-objective optimization reach the optimal value at the same time, and can only coordinate and compromise between them to make each sub-objective optimized as much as possible. Therefore, in practical applications, it is generally impossible to obtain the exact Pareto optimal solution set. Therefore, the research focus of multi-objective optimization is to obtain an approximate solution for the Pareto optimal solution set. And the genetic algorithm can coordinate the relationship between the sub-objectives to find out the optimal solution set that makes each sub-objective reach a relatively large (or relatively small) value as far as possible, i.e., the Pareto frontier can be obtained. The structure of the algorithm used in this paper is shown in Fig. 4. The running parameters of the genetic algorithm are shown in Table 2.

Table 2.

Algorithm operation parameters

Item	value
Population size	150
Cross probability	0.8
Mutation probability	0.2
The number of individuals per generation	0.7

System efficiency, total annual cost, and CO2 emission rate have been selected as objective functions. According to the research object of this paper, 11 design parameters were selected for system optimization, and the feasible range of design parameters is given in Table 3.

Table 3.

The variation ranges of the design parameters

Design parameter	Variation range
High pressure steam pressure value	9.5MPa~12.5MPa
Medium pressure vapor pressure value	2.0MPa~2.6MPa
Low pressure steam pressure value	0.3 MPa ~0.9 MPa
High pressure steam temperature	<575.65°C
Reheat steam temperature	<575.65°C
Close temperature difference (high, medium, low pressure)	6~18°C
Node temperature difference (high, medium, low pressure)	10~22°C

4.2.2

Analysis of results

Considering the three objectives, namely system efficiency, total annual cost, and CO2 emission rate, the genetic algorithm was applied to the multi-objective optimization of the system, and the resulting Pareto frontier is shown in Figure 5. With the increase in the total annual cost, the efficiency of the system gradually increases, and the carbon dioxide emission rate gradually decreases. Economic performance is getting worse while environmental performance is getting better, which means there is a contradiction between system efficiency, total annual cost, and CO2 emission rate.

Therefore, the Pareto front obtained by the above three-objective optimization can be projected onto three planes to obtain the results of dual-objective optimization. Fig. 6~8 show the Pareto frontier of dual-objective optimization of total cost and CO2 emissions, total cost and system efficiency, and CO2 emissions and system efficiency, respectively.

Figure 6 shows that the CO2 emission rate decreases with increasing total annual costs and is approximately linear. When the CO2 emission rate gradually increases to 352.69kg/(MW·h), the total annual cost decreases rapidly, and when the CO2 emission rate continues to increase, the reduction rate of the total annual cost slows down.

As can be seen from Figure 7, the total annual cost increases with the increase of system efficiency, showing an approximate linear relationship. When the system efficiency is less than 58%, the increase in the total annual cost is small, and when the system efficiency continues to increase, the increase in the total annual cost is very fast.

Figure 8 shows that the CO2 emission rate decreases as the system increases, with an approximate linear relationship. Since the optimization goal of this paper is to maximize the total annual cost and carbon dioxide emission rate of the system, it can be found that the system efficiency and carbon dioxide emission rate are consistent in the optimization process, and the total annual cost conflicts with the other two objective functions. Therefore, due to the conflict of different objective functions, there is no single scheme to make all objective functions reach the optimal value, and the optimization scheme needs to be obtained from the Pareto frontier according to the requirements.

It can be seen that the total annual cost is inconsistent with the system efficiency and carbon dioxide emission rate, and in the case of dual-objective optimization, the carbon dioxide emission rate also tends to be optimal when the system efficiency is optimized. Therefore, according to the requirements of this paper, the optimization scheme only needs to be obtained from the Pareto frontier of the system and the total annual cost. Point C corresponds to the smaller system and the lowest total annual cost, while point A corresponds to the largest system and the highest total annual cost. Suppose there is an equilibrium point, and the objective function corresponding to this point will be optimized, i.e., the system is the highest, and the total annual cost is the lowest. In fact, the equilibrium point does not exist, so in a general multi-objective problem, point B, which is closest to the equilibrium point on the Pareto front, can be considered as the final optimization solution.

As shown in Figure 7, the Pareto boundary is relatively flat between the BC segments of the Pareto front. However, in the AB section of the Pareto front after point B, the total annual cost increases significantly in the pursuit of higher system and lower CO2 emission rates. As a result, a small change in the system of the AB section of the Pareto front obtained in this paper will lead to a sharp change in the total annual cost, and point B, which is closest to the equilibrium point, is not suitable as the final optimization scheme. This is because when the efficiency of the system increases slightly, the total annual cost increases rapidly, which cannot meet the economic performance requirements of the unit. Therefore, the point where the objective function value is better than the original design parameters can be selected from the BC segment of the Pareto front as the final optimization scheme according to preference.

The design parameters of the optimized system selected in this paper are shown in Table 4. The energy, economic, and environmental performance indicators of the optimization points are shown in Table 5. Compared with the original system, the values of these three objective functions of the optimized system are improved. Among them, the system efficiency is 56.82%, the carbon dioxide emission rate is 353.31kg/(Mwah), and the total annual cost is 6.45×10⁸¥. The optimized system has obvious advantages in terms of environmental protection and energy saving, and can provide guidance for the application of pollutant emissions.

Table 4.

Optimized system and original system design parameters

Design parameter	Unit	Before optimization	After optimization
High pressure steam pressure value	MPa	9.921	12.568
Medium pressure vapor pressure value	MPa	2.13	2.05
Low pressure steam pressure value	MPa	0.38	0.27
High pressure steam temperature	°C	567.6	573.69
Reheat steam temperature	°C	576.6	576.98
Close temperature difference (high, medium, low pressure)	°C	15/15/15	17/8/8
Node temperature difference (high, medium, low pressure)	°C	8/8/8	5/9.5/5

Table 5.

The energy, economic and environmental performance of the optimization point

Target function	Unit	Before optimization	After optimization
Cost	¥	6.52×10⁸	6.45×10⁸
H_GTCC	%	56.31	56.82
ε	kg/(MW·h)	357.61	353.31

Table 6 shows the design parameters of the optimized natural gas-internal combustion engine combined cycle system at different natural gas turbine load rates. Table 7 shows the energy, economic, and environmental performance metrics of the optimized point at different natural gas turbine load rates. The objective function values of the optimized system are improved at different natural gas turbine load rates. Particularly on the economic side, the total annual cost increases when the load factor of the natural gas turbine decreases. As the natural gas turbine load factor decreases, the difference in total annual cost between the optimized system and the original system increases, which indicates that the optimized system design parameters can greatly improve the economic performance of the system under off-design conditions.

Table 6.

Design parameters of the optimized points

Design parameter	Unit	Internal combustion turbine load rate
Design parameter	Unit	100%	75%	50%	30%
High pressure steam pressure value	MPa	12.381	12.5	12.5	12.5
Medium pressure vapor pressure value	MPa	2.11	2.126	2.11	2.11
Low pressure steam pressure value	MPa	0.284	0.265	0.225	0.21
High pressure steam temperature	°C	576.36	571	571	558.56
Reheat steam temperature	°C	576.96	576.8	576.8	5573.879

Table 7.

Comparison of energy, economic, environment performance after optimization

Target function	Unit	Before optimization
Target function	Unit	100%	75%	50%	30%
Cost	¥(×10⁸)	6.52	5.62	4.45	3.61
H_GTCC	%	56.31	53.38	50.37	45.98
ε	kg/MW·h	357.61	376.79	400.65	442.3
Target function	Unit	After optimization
Target function	Unit	100%	75%	50%	30%
Cost	¥(×10⁸)	6.45	5.52	4.36	3.24
H_GTCC	%	56.82	54.36	51.87	47.66
ε	kg/MW·h	353.31	368.39	387.54	421.74

Since the multi-objective optimization results in this paper add the environmental cost to the total annual cost, which can reflect the economic loss caused by pollutant emissions to the environment. In turn, the total annual cost can reflect the actual economic situation of the unit. Therefore, in order to further obtain the economic losses caused by pollutant emissions to the environment and to reflect the advantages of the optimized system in terms of environmental protection and energy saving, it is necessary to calculate the proportion of the environmental part of the system’s cost in the total annual cost. The increased cost due to considering the environment has been internalized and the environmental component of the system environmental thermoeconomics cost reflects the cash cost consumed by the system to generate 1kWh due to environmental impacts. For the purposes of this paper, the ratio of the environmental component of the system to the total annual cost is defined as the environmental pollution factor, or β.

β reflects the annual cash cost consumed by the natural gas-Internal Combined Cycle system due to environmental impacts as a percentage of the total cost. Figure 9 shows the variation of β. The figure shows that environmental costs cannot be ignored. When the natural gas turbine load factor decreases, β increases. And the difference of β of the system before and after optimization increases as the natural gas turbine load factor decreases, which indicates that the optimized system design parameters can greatly reduce the economic loss of the system to the environment due to pollutant emissions, which has good environmental and economic benefits. Moreover, the optimized system reduces the emission of pollutants to meet the requirements of environmental protection.

5

Conclusion

The main aim of this study is to construct a hydrodynamic model of air pollutant emissions from a natural gas internal combustion engine and solve it using a multi-objective optimization algorithm using an undominated genetic algorithm.

Under the benchmark working conditions of a natural gas internal combustion engine, the numerical simulation error of each index is controlled within 5%. Among them, the transient temperature of the combustion engine is about 1940K, and the NOx generation presents two hump-like distributions, and it is found that the inlet temperature is closely related to the natural gas flow rate, high-temperature gas distribution, flame fluctuation, pressure pulsation, etc. The results of multi-objective pollutant control show that there is a contradiction between system efficiency, total annual cost, and carbon dioxide emission rate. The system efficiency of the scheme is 56.82%, the carbon dioxide emission rate is 353.31kg/(Mwah), and the total annual cost is 6.45×10⁸¥, which achieves better economic and environmental benefits compared with the pollutant emission situation before optimization.

Therefore, this paper realizes the simulation and optimization of the air pollutant emission process through the numerical simulation of natural gas hydrodynamic model and the multi-objective control method of pollutant emission, which can help the natural gas internal combustion engine to obtain better environmental and economic benefits.

Language:: English

Publication timeframe:: 1 times per year
Journal Subjects:: Life Sciences, Life Sciences, other, Mathematics, Applied Mathematics, General Mathematics, Physics, Physics, other

Journal RSS Feed

Numerical Modeling and Multi-Objective Control of Air Pollutant Emissions from Natural Gas Internal Combustion Engines

Xiaoqing Wu

Jun Chen

Jianing Kang

Yunwei Li

Published Online: Mar 24, 2025

Received: Oct 09, 2024

Accepted: Feb 06, 2025

DOI: https://doi.org/10.2478/amns-2025-0784

KeywordsNumerical simulation, Multi-target control, Natural gas internal combustion engines, Pareto optimal, TOPSIS decision-making

© 2025 Xiaoqing Wu et al., published by Sciendo

This work is licensed under the Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License.

Keywords
Numerical simulation, Multi-target control, Natural gas internal combustion engines, Pareto optimal, TOPSIS decision-making