energies
Article
Fired Heaters Optimization by Estimating Real-Time
Combustion Products Using Numerical Methods
Ricardo Sánchez 1 , Argemiro Palencia-Díaz 1,* , Jonathan Fábregas-Villegas 1,2 and Wilmer Velilla-Díaz 3,*
1 Department of Mechanical Engineering, Universidad Tecnológica de Bolivar, Cartagena 130001, Colombia;
sanchezr@utb.edu.co (R.S.); jfabregas@utb.edu.co (J.F.-V.)
2 Industrial Engineering Program, Universidad Sergio Arboleda, Barranquilla 081007, Colombia
3 Department of Mechanical Engineering, Universidad de La Serena, La Serena 1720170, Chile
* Correspondence: argpalencia@utb.edu.co (A.P.-D.); wilmer.velilla@userena.cl (W.V.-D.)
Abstract: Fired heaters upstream of distillation towers, despite their optimal thermal efficiency, often
suffer from performance decline due to fluctuations in fuel composition and unpredictable operational
parameters. These heaters have high energy consumption, as fuel properties vary depending on the
source of the crude oil. This study aims to optimize the combustion process of a three-gas mixture,
mainly refinery gas, by incorporating more stable fuels such as natural gas and liquefied petroleum gas
(LPG) to improve energy efficiency and reduce LPG consumption. Using real-time gas chromatography-
mass spectrometry (GC-MS) data, we accurately calculate the mass fractions of individual compounds,
allowing for more precise burner flow rate determinations. Thermochemical data are used to calculate
equilibrium constants as a function of temperature, with the least squares method, while the Newton–
Raphson method solves the resulting nonlinear equations. Four key variables (X4, X6, X8, and X11),
representing H2, CO, O2, and N2, respectively, are defined, and a Jacobian matrix is constructed to ensure
convergence within a tolerance of 1 × 10−6 over a maximum of 200 iterations, implemented via Python
3.10.4 and the scipy.optimize library. The optimization resulted in a reduction in LPG consumption
by over 50%. By tailoring the fuel supply to the specific thermal needs of each processing unit, we
achieved substantial energy savings. For instance, furnaces in the hydrocracking unit, which handle
cleaner subproducts and benefit from hydrogen’s adiabatic reactions, require much less energy than
those in the primary distillation unit, where high-impurity crude oil is processed.
Citation: Sánchez, R.; Palencia-Díaz,
A.; Fábregas-Villegas, J.; Velilla-Díaz, Keywords: optimizing combustion; adiabatic flame; Newton–Raphson; fired heaters; refinery gas
W. Fired Heaters Optimization by
Estimating Real-Time Combustion
Products Using Numerical Methods.
Energies 2024, 17, 6190. https:// 1. Introduction
doi.org/10.3390/en17236190 In a wide range of refinery configurations, heaters are strategically positioned up-
Academic Editor: Qingsong Wang stream of the crude distillation tower to ensure that the incoming crude oil is heated to
a specific, required temperature. This temperature regulation is imperative for optimal
Received: 28 August 2024
operation. In instances where the crude oil does not attain this essential temperature, it
Revised: 20 September 2024
becomes necessary to either increase the combustion of fuel gas or adjust the calorific value
Accepted: 29 September 2024
Published: 9 December 2024 of the gas mixture. This adjustment is essential to achieve the desired process temperature,
underscoring the importance of thermal management and energy efficiency in refinery
operations. Fired heaters represent the primary energy consumers in both the refining
and petrochemical sectors, with approximately 55% of total energy consumption being at-
Copyright: © 2024 by the authors. tributed to their function as heat exchangers [1]. Although most were designed for thermal
Licensee MDPI, Basel, Switzerland. efficiencies of up to 80%, current operating efficiencies often fall short of expectations [2].
This article is an open access article This is thought to result from the inconsistent operation of the furnaces in line with their
distributed under the terms and design conditions.
conditions of the Creative Commons Hydrocarbons have played a crucial role in the global economy for centuries, acting
Attribution (CC BY) license (https:// as key drivers of industrial development [3–6]. Hydrocarbons in the gaseous state serve as
creativecommons.org/licenses/by/ a common fuel burned in this equipment [7,8]. However, exclusive reliance on one type of
4.0/).
Energies 2024, 17, 6190. https://doi.org/10.3390/en17236190 https://www.mdpi.com/journal/energies
Energies 2024, 17, 6190 2 of 14
gas is not feasible due to source limitations. This leads to increased costs and insufficient
combustion efficiency. Currently, a mixture of natural gas (NG) and refinery gas (RG) is
used as the primary fuel source in industrial applications [9]. RG is usually recovered from
various processes, including cracking, desulphurization, and catalytic reforming units [10].
Liquefied petroleum gas (LPG), a by-product of crude oil refining, is occasionally used
as a third gas due to its high propane and butane content, which increases the calorific
value of the mixture. This is particularly useful when the mixture of NG and RG has a
significantly low calorific value or when considering the fluctuating costs of these fuels in
the market [11]. These procedures generate a diverse matrix of fuel blends, resulting in a
wide range of compositions suitable for combustion processes [12]. Highlighting that NG
and LPG are relatively more environmentally friendly, their combustion produces lower
pollutant emissions into the atmosphere [13].
The optimization of thermal efficiency in fired heaters can be achieved through differ-
ent strategies [14–16]. Thorat and Garg, both emphasize the importance of optimal design
and application of heat tracing and management systems, which can significantly reduce
energy consumption [2,17]. Masoumi further underscores the potential for improvement in
refinery furnaces’ efficiency through mathematical modeling, particularly by considering
ambient and operational conditions [18]. These studies collectively highlight the need for a
holistic approach to addressing the challenges of energy usage in petrochemical refineries.
It has also been emphasized that net-zero global carbon dioxide CO2 emissions need to
be reached to achieve this goal by 2050 [19]. Thus, reducing greenhouse gas emissions,
especially CO2, is one of the main global challenges to achieving a more sustainable fu-
ture. Combustion furnaces within oil refineries generally account for over 65% of CO2
emissions [20]. Single-burner research furnaces have undergone testing, revealing staged
combustion air injection and flue gas recirculation as the most promising combustion mod-
ifications to decrease NO emissions from the gas-fired process heater. The modifications
have a potential for up to a 67% reduction [21].
Research also highlights the significance of fuel composition, particularly hydrogen
enrichment, in enhancing combustion efficiency. Additionally, experts are developing
simulation methods for optimizing gaseous fuel mixtures [22–26]. Saifullin’s specialized
techniques have significantly improved combustion efficiency in thermal power plants
that use variable fuel compositions [27]. Some researchers have studied real-time opti-
mization [28]. However, there is still a gap in the real-time optimization of fuel mixtures,
particularly for RG, when responding to their changing compositions. Gas chromatography-
mass spectrometry (GC-MS) data play a vital role in comprehending RG behavior. In an
industrial context, it is important to maintain a constant calorific value, which is critical for
cost management as outlined by Cote [29]. The imperative to maintain a consistent calorific
value drives many end-users to monitor and regulate their flue gas quality [30]. Chomiak
observes that the utilization of gas with higher calorific value can minimize expenses linked
to fuel consumption, thereby affecting the profitability per cubic meter of refined crude
oil [31,32].
This study aims to bridge the knowledge gap by developing a method to dynamically
optimize the combustion of a three-gas mixture. We utilize GC-MS data to analyze RG
in real time, determining the mass fraction of individual compounds. These data enable
the selection of optimal flow rate combinations for burners, ensuring efficient combustion.
Equilibrium constants are also calculated as a function of temperature by applying the
Newton–Raphson method in Python to solve the resulting nonlinear equations [33]. The
research extends beyond theoretical analysis, incorporating simulations validated with in
situ chromatography data. The novel aspect of this research lies in its real-time approach
to optimizing fuel mixtures, examining the inherent variability of RG. This approach is
expected to contribute to a significant reduction in LPG consumption, thereby offering
substantial economic benefits to oil companies. Additionally, the study contributes to
the broader goal of enhancing the overall efficiency and performance of furnaces in the
petrochemical industry.
Energies 2024, 17, 6190 3 of 14
The paper is structured as follows: Section 2 details the setup of chemical reactions
involved in combustion and the dynamic model of the gas mixture. Section 2.2 discusses
the reduction in equations from twelve variables to four, which is essential for accurately
computing combustion product mass fractions. Section 3 presents the method’s numeri-
cal performance, supported by simulations and real-time chromatography data. Finally,
Section 4 offers conclusions, highlighting the study’s contributions to efficient energy
management in fired heaters.
2. Methodology and Simulation Model
2.1. Composition and Properties of Gas Mixture
This research addresses the combustion of a three-gas mixture, with a particular focus
on how the varying composition of RG (shown as a dashed yellow line in Figure 1) affects
the process, alongside NG and a standard gas mixture (depicted by the green and red
lines, respectively).
1050 Peak at sample RG#87: 1071.28 Btu/ft
3
Sample RG#176: 952.62 Btu/ft3
1000
950
900
850
100% Refinery gas (RG)
800 100% Natural gas (NG)
750 Min at sample RG#29: 743.52 Btu/ft3 70%RG-30%NG Mixture
0 50 100 150 200 250
Chromatographic daily test
Figure 1. Daily Gas Sample Data Recorded and Analyses.
One of the key properties of gases is their ability to mix uniformly with each other,
resulting in a solution where each gas component can be analyzed independently while
maintaining the same temperature and volume within the mixture [34]. Understanding gas
mixtures becomes more straightforward with knowledge of their composition, which we
investigated through a review of statistical reports from approximately 250 chromatography
analyses of RG streams collected every 24 h. Table 1 outlines the composition of three
specific refinery gases at 29, 87 and 176 days, thus: RG29, RG87, and RG176 (See Figure 1),
which represent a range of refinery gas compositions with varying calorific values and
sulfur content. These gases were selected due to their distinct low heating value (LHV),
thus, the minimum LHV RG29, the maximum LHV RG87, and the last one with a typical
sulfur content in RG176. In addition, the composition of natural gas (NG) is predominantly
methane (98%), while LPG is primarily composed of butane (98%), making both ideal
for numerical modeling purposes due to their stable and predictable compositions. Our
code is developed to flexibly process chromatographic data for these gases across various
process parameters.
For the calculation of the thermodynamic properties of the air–fuel mixture, the
engineering equations solver (EES) software (https://fchartsoftware.com/ees/, accessed
on 19 September 2024) was used, where NASA’s ideal gas data were taken. These data
consist of specific heat, specific enthalpy and density, among others, at standard pressure as
a function of temperature. Another advantage is the ease of numerically solving thousands
of coupled non-linear algebraic and differential equations.
But first, as we observe in Table 1, none of the refinery gas samples yield the sum of the
volumetric composition equal to 1 as the other two gases, so it must be readjusted to give
exactly 1. Implementing them guarantees the precision of the solution for a great variety of
systems, in a wide range of conditions with a high degree of efficiency and reliability.
Lower heating value (Btu/ft3)
Energies 2024, 17, 6190 4 of 14
Table 1. Selecting gases for the numerical model (MF = Mole Fraction).
Fuels RG29 RG87 RG176 NG LPG
Comp. MF MF MF MF MF
H2 0.16983526 0.28839122 0.15382659 0.0 0.0
CO2 0.00355668 0.0018357 0.00180876 0.0 0.0
CH4 0.50036226 0.49431918 0.63173331 0.9831 0.0
C2H6 0.1371265 0.04462025 0.11880152 0.00258124 0.9838
C2H4 0.05565108 0.03992771 0.04821289 0.0 0.0
C3H8 0.0549727 0.0147482 0.00479593 0.002 0.0054
propylene 0.02094179 0.00305931 0.00649243 0.0 0.0098
H2S 0.0 0.0 0.00327537 0.0 0.001
C2H2 7.54 × 10−6 0.0 6.62 × 10−6 0.0 0.0
isobutane 0.00558021 0.004937739 0.00230704 0.0 0.0
propadiene 0.0 0.0 0.0 0.0 0.0
n-butane 0.00679208 0.00343208 0.00306676 0.0 0.0
O2 0.0004077 0.0153073 0.00022166 0.001 0.0
trans-2-butene 0.0032849 0.0023618 0.0013005 0.0 0.0
N2 0.0212296 0.080129 0.0109625 0.011 0.0
1-butene 0.0032146 0.00032 0.0010858 0.0 0.0
isobutene 0.0042796 0.0008848 0.0009842 0.0 0.0
2-butene 0.0022651 0.0015048 1.32 × 10−5 0.0 0.0
isopentene 0.0020144 0.000221 0.0022329 0.0002048 0.0
n-pentane 0.00162679 5.559 × 10−5 0.0021111 6.879 × 10−5 0.0
1,3-butadiene 0.0001166 0.0 4.06 × 10−5 0.0 0.0
CO 0.00653259 0.0039111 0.0067198 0.0 0.0
Hexans+ 0.0002015 3.27 × 10−5 0.0 3.03 × 10−5 0.0
∑ Xi ≡ 1 0.9999994 0.99999949 0.99999948 0.99998514 1.0
2.2. Chemical Reactions in the Combustion Process
Considering a fuel with a composition of C, H, O, N and S, mixed with air at an equiv-
alence ratio denoted as ϕ, and examining its reaction in the framework of equilibrium
thermodynamics applied to the system, we observe the formation of products at a tem-
perature set as T and a pressure set as P [35]. Under conditions where high temperatures
induce dissociation, up to 12 combustion products can be generated [36].
The entire representation of the chemical reaction can be formulated as follows,
X14[Cn1 + Hm1 + Ol1 + Nk1 + Sj1] + X15[Cn2 + Hm2 + Ol2 + Nk2 +[ Sj2] ]
X C H O N S (n+j+
m − l )
+ 16[ n3 + m3 + 4 2l3 + k3 + j3] + ϕ (O2 + 3.76N2) ⇒
X1H + X2O + X3N + X4H2 + X5OH + X6CO + X7NO
+X8O2 + X9H2O + X10CO2 + X11N2 + X12SO2 (1)
where j, k, l, m and n are the quantities of S, N, O, H and C, respectively, for the gases
1 = X14, 2 = X15, and 3 = X16. In this context, X14, X15, and X16 represent the quantities
of refinery gas, natural gas, and liquefied petroleum gas, respectively, constituting the
entirety of the fuel blend. The values ranging from X1 to X12 are the mole fraction of
the combustion products. For the sake of simplicity, we treat air as a mixture consisting
solely of O2 and 3.76N2, disregarding other components. The quantity of air utilized
in the combustion process depends directly on the fuel, and the concentration of air
is defined by [(n + j + m/4 − l/2)/ϕ](O2 + 3.76N2), where, j = X14 j1 + X15 j2 + X16 j3;
k = X14k1 + X15k2 + X16k3; l = X14l1 + X15l2 + X16l3; m = X14m1 + X15m2 + X16m3 and
n = X14n1 + X15n2 + X16n3.
Energies 2024, 17, 6190 5 of 14
The equivalence ratio between air and fuel, denoted as ϕ, is defined as follows:> 1, for fuel-rich mixtures,
ϕ = = 1, for a stoichiometric mixture,
< 1, for mixtures that are fuel-lean.
The excess air coefficient is another commonly used parameter for characterizing
the stoichiometry of the fuel-air mixture, and it is related to the equivalence ratio as
(1 − ϕ)/ϕ × 100%. By simplifying the expression, we used X13 as the sum of all the
reactants, r0 = (n + j + m/4 − l/2)/ϕ, r1 = l/2 + r0, and r2 = k/2 + 3.76r0. Therefore, we
arrive at the following equation:
X13(nC +m H +r1 O +r2 N +j S) ⇒
X1H + X2O + X3N + X4H2 + X5OH + X6CO + X7NO + X8O2 + X9H2O + X10CO2 + X11N2 + X12SO2 (2)
This represents the atom balance for the general equation.
C : nX13 = X6 + X10 (3)
H : mX13 = X1 + 2X4 + X5 + 2X9 (4)
O : 2r1X13 = X2 + X5 + X6 + X7 + 2X8 + X9 + 2X10 + 2X12 (5)
N : 2r2X13 = X3 + X7 + 2X11 (6)
S : jX13 = X12 (7)
In addition, a condition is introduced in this system which requires that the total sum
of the mole fractions of all products be equal to one mole. Therefore, this implies:
12
∑ Xi = 1 (8)
i=1
To solve the system of six equations (from Equations (3)–(8)) with 13 unknowns, an
additional set of seven equations is required (from Equations (10)–(16)). These equations
are sourced from the equilibrium reactions.
Gas Equilibrium Constants
The equilibrium constant (Kp) is determined by the ratio of the molar concentra-
tions (mol/L) of reactants and products in a chemical reaction. Its value is temperature-
dependent and must always be specified [37]. We begin by examining the general chemical
reaction to determine the equilibrium constants. The equation that defines equilibrium
constants as a function of partial pressures for a given combustion reaction is as follows:
Πi(Pi/P)µiKp = (9)Πj(Pj/P)
µj
where the partial pressures of flue gases are represented by Pi/P = Xi and the respective
stoichiometric coefficients by µi. The seven missing equations to solve the system corre-
spond to the seven chemical reactions associated with natural gas, according to Olikara as
follows [34].
1
H ⇒ XH K 1
2 2 p1
= 1 (10)
X 24
Energies 2024, 17, 6190 6 of 14
1 X
O2 ⇒O Kp2 = 22 1 (11)X 28
1 X
N 3
2 2
⇒N Kp3 = 1 (12)
X 211
1 1 X
O2 + H
5
2 2 2
⇒OH Kp5 =
X1/2X1/2
(13)
4 8
1 1 X
O + N 7
2 2 2 2
⇒NO Kp7 = 1/2 1/2 (14)X8 X11
1 ⇒ XH2 + O2 H2O K = 92 p9 X X1/2
(15)
4 8
1 X
CO + O ⇒CO K = 10 (16)
2 2 2 p10 X 1/26X8
Applying Equation (9), the equilibrium constants data were taken from JANAF Ther-
mochemical Tables [34], where log10 Kp formation for all species are tabulated as functions
of the absolute temperature (K). Theoretical investigations [38] propose a functional rela-
tionship of this nature to compute the Kp as follows:
B
log Kp = A ln TA + + C + DTA + ET2A (17)TA
where a transformed temperature TA defined as 0.005 T/9 (T is in °R) was used for fitting.
The constants A, B, C, D and E are listed in Table 2. This model was used to fit tabulated
data through a least squares fitting method. To strike a balance between precision and
temperature range, we chose the 1500 to 3000 K (2700 to 5400 °R) range for studying
combustion purposes [39]. The adiabatic flame temperature, representing the maximum
possible temperature without heat exchange with the surroundings, is crucial for optimizing
combustion. Understanding this temperature helps minimize harmful emissions, maximize
efficiency, and determine the ideal air–fuel mixture and fuel blends. Figure 2 illustrates
a case of combustion wherein the adiabatic flame would attain a peak temperature of
approximately 2000 K (3600 °R). Consequently, Kp values are computed simultaneously for
each of the reactions.
106 H
O
kp = 3.54 × 103
103 kp = 7.85 × 102
N
OH
NO
100 kp = 5.59 × 10 1 H2O
CO
kp = 2.00 × 10 2 2
3
10 3 kp = 6.74 × 10 4
kp = 1.64 × 10
10 6
10 9 kp = 9.49 × 10
10
10 12
3000 3500 4000 4500 5000 5500
Temperature (°R)
Figure 2. Equilibrium Constants Kp at 3600 °R.
Kp
Energies 2024, 17, 6190 7 of 14
Table 2. Constant coefficients of Reactions.
Equation A B C D E
(10) 0.432168 −0.112464 × 102 0.267269 × 101 −0.7457 × 10−1 0.242484 × 10−2
(11) 0.310805 −0.129540 × 102 0.321779 × 101 −0.7383 × 10−1 0.344645 × 10−2
(12) 0.389716 −0.245828 × 102 0.314505 × 101 −0.9637 × 10−1 0.585643 × 10−2
(13) −0.141784 −0.213308 × 101 0.853461 0.355015 × 10−1 −0.3102 × 10−2
(14) 0.150879 × 10−1 −0.470959 × 101 0.646096 0.272805 × 10−2 −0.1544 × 10−2
(15) −0.752364 0.124210 × 102 −0.260286 × 101 0.259556 −0.1626 × 10−1
(16) −0.4153 × 10−2 0.148627 × 102 −0.475746 × 101 0.124699 −0.9002 × 10−2
The log Kp values computed from the equations were compared with the original
data, showing deviations of less than 0.0009, these deviations are considered negligible.
In the process of data collection, we observe the following express√ions, where the√val-
ues of the√variables are in√fluenced by spec√ific constants: X1 = √Kp1 X4, X2 = Kp2√X8,
X3 = Kp3 X11, X5 = Kp5 X4X8, X7 = Kp7 X8X11, X9 = Kp9X4 X8, X10 = Kp10X6 X8.
When expressing the molar fractions using the equilibrium constants in Equations (3)–(7),
we decrease the number of constraints, resulting in a fresh system of four nonlinear equa-
tions with four unknowns. This means that every fraction is depending only on X4, X6, X8,
and X11 (H2√, CO, O2, and N2)√exclusively as foll√ows: √
Kp1√ X4 + 2X4√+ Kp5 X4X8 + 2K√p9X4 X8 − d1(X6 + Kp√10X6 X8) = 0 (18)
Kp2 X8 + Kp5 X4X8 + X6 + Kp7 X8X11 + 2X8 + Kp9X4 X8 + AA = 0 (19)
√ √
where AA = 2Kp10√X6 X8 − d√2(X6 + Kp10X6 X8) √
√ Kp√3 X11 + K√p7 X8X11 + 2X11√− d3(X6 + Kp10X6 √X8) = 0 (20)
Kp1 X4 + Kp2 X8 + Kp3 X11 + X4 + Kp5 X4X8 + X6 + Kp7 X8X11 + BB = 0 (21)
√ √ √
where BB = X8 + Kp9X4 X8 + Kp10X6 X8 + X11 + d4(X6 + Kp10X6 X8) − 1, where
d1 = m/n, d2 = 2r0/n, d3 = 2r1/n, and d4 = r2/n. This set of four interrelated non-
linear equations can be expressed as a function involving four variables: fi(X4, X6, X8, X11),
here i = 1, 2, 3, 4.
To linearize Equations (18)–(21), a Taylor series expansion is applied, yielding the
following generalized expression.
∂ f ∂ f ∂ f ∂ f
fi +
i ∆X + i ∆X i i
X 4 X 6
+ ∆X
X 8
+ ∆X11 = 0 (22)∂ 4 ∂ 6 ∂ 8 ∂X11
2.3. Numerical Modeling of the Combustion Process
Using real-time gas chromatography–mass spectrometry (GC-MS) data, we accurately
calculate the mole fractions of individual compounds, allowing for more precise burner
flow rate determinations. Thermochemical data are used to calculate equilibrium constants
as a function of temperature, with the least squares method, while the Newton–Raphson
method solves the resulting nonlinear equations. The system of four nonlinear equations
is solved using a 4 × 4 Jacobian matrix of first-order derivatives to ensure convergence
within a tolerance of 1 × 10−6 over a maximum of 200 iterations (See Equation (23)). The
optimization method is implemented in Python using the scipy.optimize library and the
newton() function.  ∂ f1 ∂ f ∂ f ∂ f 1 1 1     ∂X4 ∂X6 ∂X8 ∂X11 ∆X − f ∂ f2 ∂ f2 ∂ f
4 1
2 ∂ f2 
 ∂X4 ∂X X X  ∆X6   6 ∂ 8 ∂  − f11 2∂ f3 ∂ f ∂ f ∂ f  = (23)3 3 3 ∆X∂X ∂X ∂X ∂X 8  − f34 6 8 11 
∂ f4 ∂ f4 ∂ f4 ∂ f4 ∆X11 − f4
∂X4 ∂X6 ∂X8 ∂X11
Energies 2024, 17, 6190 8 of 14
The matrix equation in the expanded form is as follows:
 √ √

 √kp1 kp5√ X8 √ √ d k X k X k X  + √ + 2kp9 X 2 −d k X − d −
1 √p10 6 p5√ 4 √p9 4+ + + 0 
 2 X 2 X 8 1 p10 8 1  4 4
2 X8 √2 X8 √ X8 √ 
k 
 p5√
X8 √ √ d k X k k X k X k X k X 
 + kp9 X8 d d k X
6 p10 6 p2 p5 4 p7 11 p9 4 p7 8
+
2 X 5 6 √p10 8
√ + √ + √ + √√ + √ + 2 √  4 2 X8 2 X8 2 X8 2 X8 2 X8 2 X 11√

 d k X k X k k X  √ 0 −d3kp10 X8 − d −
3 √p10 √6 p7√ √11 p3 p7 8 3 + √ + √ + 2 √ √ 2 X8 2 X8 2 X 2 X  √kp1 kp5√ X8 d7k√p10 X6 √kp2 kp5√ X4 kp7√ X11 kp√
11 √11 
9 X4 √k k X + + kp9 X8 + 1 d k X d 1 p3 p7 8+ + + + + + + √ 2 X4 2 X 7 p10 8 74 2 X8 2 X8 2 X8 2 X8 2 X8 2 X11 2 X 11
 ∆X4 

 
 − (X1 + 2X4 + X5 + 2X9) + d1(X6 + X10)
 ∆X6  =  −(X2 + X5 + X6 + X7 + 2X8 + X9 + 2X10)− d2(X6 + X10)∆X8 −(X3 + X7 + 2X11) + d3(X6 + X10) 
∆X 1 − (∑1111 i=1 Xi)− d4(X6 + X10)
Initial Value Estimation
To estimate initial values, we focus on the products generated during complete com-
bustion. This involves narrowing the scope from the original 12 products to focus on only
four key ones, specifically H2O, CO2, N2, and SO2. Additionally, we consider the products
H2, CO, O2, and N2, which correspond to the fractions X4, X6, X8, and X11, respectively.
This approa(ch establishes the followin)g products within the combustion equation:
X13 nC +m H +r O +r′ N +j S ⇒ X4H2 + X6CO + X8O2 + X9H2O + CC (24)
where CC = X10CO2 + X11N2 + X12SO2.
When performing the elemental balance for each constituent of the fuel in this par-
ticular case, we have C : nX13 = X6 + X10, H : mX13 = 2X4 + 2X9, O : 2r1X13 =
X6 + 2X8 + X9 + 2X10 + 2X12, N : 2r2X13 = 2X11, S : jX13 = X12
By substituting the fractions using the constants from Section 2.2, we obtain
the following:
nC X
X 10 136 = √ (25)C10 + X8
1
2 mC5√XX 134 = (26)C5 + X8
X11 = r2X13 (27)
X12 = jX13 (28)
√ √
n(C10 +√2 X )
1 m X
0 8 2 √ 8 2X= + + 8 + 2j − 2r (29)
C10 + X8 C
1
5 + X8 X13
The quantity of X13 can be accurately estimated to ensure that the sum of the molar
fractions equals one.
12
X4 + X6 + ∑ Xi = 1 (30)
i=8
When the air–fuel equivalence ratio is less than or equal to 1, we can obtain a reliable
estimate of X13 through complete combustion.
1
X13 = m (31)
4 + r1 + 2r2
Energies 2024, 17, 6190 9 of 14
Conversely, when the air–fuel equivalence ratio exceeds 1, a dependable estimation of
X13 can be derived from incomplete combustion.
1
X13 = m (32)n + 2 + r2 + j
By replacing the estimated value of X13 into Equation (29), we obtain X8. Once
Equation (29) is solved, the remaining unknowns can be readily determined through
substitution into Equations (25)–(28). This approximation is used as the initial value for
the model.
X(1), X(1), X(1), X(1)[ 4 6 8 11 ] (33)
The previous vector approximates closely to the solution vector:
X(∗)[ 4 , X
(∗)
6 , X
(∗), X(∗)8 11 ] (34)
At each iteration, the updated vector is used to compute the partial derivatives and
evaluate the functions until the two criteria of convergences (200 iterations or tolerance of
1 × 10−6).
3. Results
3.1. Numerical Model Validation
In the research, the numerical model of the combustion process, and its validation
through operating data, considered crucial parameters such as LHV, density, adiabatic
flame temperature and CO emissions. The model postulates that emissions of CO2 and
other greenhouse gases are inextricably linked to ϕ. The algorithm begins with ϕ = 1,
which represents a stoichiometric mixture. It then optimizes between stoichiometric and
fuel-rich mixtures (ϕ ≤ 1) with the objective of reducing emissions while maximizing
energy efficiency and ensuring regulatory compliance. In the case of ideal combustion, CO2
emissions are at their highest, while CO emissions are minimized, as illustrated in Figure 3.
0.30 H2
H
0.25 O
O2
0.20 OH
H2O
0.15 CO
CO2
0.10 N
NO
0.05
0.00
0.5 1.0 1.5 2.0 2.5 3.0 3.5
Equivalence ratio, 
Figure 3. Levels fraction of harmful emissions.
Accordingly, the model’s assumptions are constrained within a range of ϕ = 1 to
ϕ = 1.5, where excess air is introduced to ensure more complete combustion, thereby
balancing CO2 production with lower levels of other harmful emissions, such as CO.
As shown in Figure 4, the density information for each chromatographic sample clearly
shows consistent alignment between the model results (blue circles) and the data reported
(red circles), thus affirming the reliability of our thermodynamic property calculations.
The error is below 0.014% among the 247 samples, providing substantial validation of the
computational model proposed in this research.
Mass fractions
Energies 2024, 17, 6190 10 of 14
0.052
0.050
0.048
0.046
0.044
0.042 Numerical model
0.040 Refinery data
0 50 100 150 200 250
Daily Chromatography analysis of RG at 15°C
Figure 4. RG Density chromatography and model results.
On the other hand, when it comes to the computed calorific values of each specific
chromatographic sample, a disparity becomes evident when compared to the reported
data in Figure 5. This examination indicates that this difference can be attributed to the
temperature conditions during their laboratory analysis, which plays a pivotal role. It is
worth noting that each sample may exhibit a temperature falling within the range of 90 °C
to 120 °C. If the exact temperatures employed for this purpose were used, this validation
would likely yield results much closer to the reported values.
1050
1000
950
900
850
800 Numerical model
750 Refinery data
0 50 100 150 200 250
Daily Chromatography Index for Composition Testing in Lab
Figure 5. RG LHV chromatography and model results.
3.2. Approaches for Achieving a 50%+ Reduction in LPG Use
The algorithm for determining optimal burner flow rates is based on the principle that,
regardless of the volume of crude oil processed at any given time (e.g., 100,000 barrels per
day), approximately 70% of the refinery’s furnace fuel demand is supplied by recovered
refinery gas. This volume, in conjunction with its physicochemical properties, which are
obtained through chromatographic analysis, serves as a fixed parameter at the outset
of the optimization process. The objective is to achieve an LHV as close as possible to
1000 BTU/ft3, in some cases without the need to incorporate LPG into the fuel mixture.
Initially, the physicochemical characterization of gases used for combustion included
determining certain properties. Table 3 presents the initial calculations of these properties
such as molecular mass, density, and LHV.
NG is the usual source for the rest unless there is a significant price difference or LHV.
In such cases, we introduce LPG as a third component to ensure optimal combustion [40].
Under these conditions, the numerical model of the mixing of the three gases must take
into account the following considerations:
1. Calculate the mass and density of the gas mixture.
2. The LHV of the mixture must be ≥ 1000 Btu/ft3.
3. Obtain the adiabatic flame temperature.
4. Reducing CO2 emissions and minimizing excess air.
LHV (Btu/ft3) Density (lb/ft3)
Energies 2024, 17, 6190 11 of 14
Once these conditions are satisfied, the resulting mixture is optimized for combustion.
Table 4 presents the blend that meets these optimization criteria.
Table 3. Physicochemical properties of gases.
Gas Mass Density Density LHV
gmol kg/m3 lb/ft3 BTU/ft3
RG29 15.48 0.6558 0.04094 743.52
RG87 20.06 0.8513 0.05314 1071.28
RG176 17.38 0.7371 0.04601 952.62
NG 16.29 0.6907 0.04311 903.98
LPG 44.23 1.9040 0.11886 2225.47
Table 4. Properties of modeled mixtures.
Mix RG-NG-LPG Mass Density LHV
%-%-% gmol lb/ft3 Btu/ft3
MRG29 70-14-16 20.19 0.0537 1003.09
MRG87 70-30-00 18.92 0.0501 1021.09
MRG176 70-25-05 18.45 0.0489 1004.10
For the mixture MRG29 including the RG gas containing the lowest calorific value, the
mixture simulation resulted in the following flow ratio: 70% RG, 14% NG and 16% LPG
(see Table 4). It was necessary to use LPG gas due to the low calorific value contained in the
refinery gas sample. In flow terms, this indicates that when the refinery gas has a calorific
value greater than 1050 BTU/ft3, it is not necessary to add LPG to the mixture. In such
cases, only 70% RG and 30% NG are kept in the mixture.
As we can see in Figure 6, when applying the Newton–Raphson method to the initial
values for the MRG29 mixture, the initial values converge after 70 iterations, although the
blue line stops acquiring negative values and the line attempts to stabilize. Consequently,
these values are replaced in the other combustion products to compute the total composition
of the molecule.
10−1
10−3 𝑋4 = 𝐻2 
𝑋6 = 𝐶𝑂 
10−5 𝑋8 = 𝑂2 
𝑋11 = 𝑁2 
10−7
10−9  
0              20              40             60             80            100           120           140
Iteration    
Figure 6. Newton Rhapson applied to MRG29 .
The lines displayed in Figure 7 reveal that gradual changes occur after each iteration.
During the initial 50 iterations, the lines extend towards the bottom of the graph in search
of negative values, which is illogical as they calculate mass or volume fractions. However,
they eventually stabilize, extending horizontally with a straight line after 60 iterations,
indicating accurate convergence of the method. The results show that the dynamic behavior
was accurate in both cases. After solving the numerical method for the three different
mixture combinations and obtaining the real values of X4, X6, X8, and X10, these values
Real Component of 𝑋𝑖   
Energies 2024, 17, 6190 12 of 14
are automatically substituted into Equations (11) to (16) to determine the exact values
for the remaining combustion products. This process completes the calculation of all
12 variables. As shown in Table 5, the sum of the molar fractions confirms that they total to
one, as expected.
Table 5. Combustion products, Numerical mixtures results MRG29 , MRG87 , MRG176 .
Comb. MRG29 MRG87 MRG176
Xi Prod. MF MF MF
X1 H 4.441 × 10−4 3.068 × 10−4 4.200 × 10−4
X2 O 1.222 × 10−6 8.294 × 10−7 1.090 × 10−6
X3 N 9.137 × 10−10 7.540 × 10−10 9.574 × 10−10
X4 H2 4.616 × 10−2 1.513 × 10−2 4.019 × 10−2
X5 OH 1.529 × 10−4 1.988 × 10−4 1.238 × 10−4
X6 CO 5.841 × 10−2 6.150 × 10−2 5.742 × 10−2
X7 NO 2.739 × 10−5 1.818 × 10−5 2.461 × 10−5
X O 1.439 × 10−6 3.496 × 10−78 2 1.034 × 10−6
X H O 1.022 × 10−1 2.133 × 10−19 2 7.422 × 10−2
X10 CO2 4.302 × 10−2 6.817 × 10−2 3.564 × 10−2
X −111 N2 7.496 × 10 6.414 × 10−1 7.913 × 10−1
X12 SO2 0.000 × 10+0 0.000 × 10+0 6.621 × 10−4
SUM 1.00 1.00 1.00
10−1
10−3
10−5
𝑋4 = 𝐻2 
𝑋6 = 𝐶𝑂 
10−7 𝑋8 = 𝑂2 
𝑋11 = 𝑁2 
10−9  
0                   20                  40                  60                  80                 100
Iteration    
Figure 7. Newton Rhapson applied to MRG87 .
4. Conclusions
This study offers a significant advance in the optimization of combustion processes in
fired heaters by integrating real-time analysis of refinery gas composition. This research
aims to establish a foundation for the integration of tools such as AI with immediate
control in equipment that regulates the flow rates of a three-gas mixture, focusing on the
variable composition of RG along with NG and LPG. This approach effectively improves
combustion efficiency and reduces fuel consumption, in particular, by reducing LPG use by
more than 50%, which offers significant economic benefits. The proposed improvement
has the potential to reduce LPG consumption by over 50% due to the current uniform
distribution of a single fuel mix (a blend of three gases) to all furnaces across the refinery,
without consideration of their specific heating requirements. By customizing the fuel mix
for each furnace based on its heating requirements, substantial savings are possible. For
instance, furnaces in the hydrocracking unit require less energy compared to those in the
primary distillation unit, which handles unrefined crude and demands higher temperatures.
Tailoring the fuel supply to match these requirements improves efficiency and can lead to
a reduction in LPG consumption of over 50%, particularly when NG prices are relatively
higher than LPG costs
Real Component of 𝑋𝑖   
Energies 2024, 17, 6190 13 of 14
The application of GC-MS data in real-time allowed for the accurate calculation of
individual compounds’ mass fractions. This innovation in monitoring and adjusting the
composition of the fuel mixture ensures a more consistent and efficient combustion process.
The successful implementation of the Newton–Raphson method in Python to solve the non-
linear equations derived from the study demonstrates the practical utility of the approach.
Furthermore, this research focus on reducing greenhouse gas emissions conforms to
worldwide sustainability targets. The numerical model helps to reduce operational costs
and contributes to improved efficiency and performance within petrochemical furnace
systems. This research represents a noteworthy advancement in achieving more sustainable
and efficient energy management within industrial processes. Additionally, it enables the
storage of mass flow data on pollutant gas emissions to report these metrics in annual
projections for reducing CO2.
Author Contributions: Conceptualization, A.P.-D., R.S. and W.V.-D.; methodology, R.S., A.P.-D. and
J.F.-V.; software, R.S. and W.V.-D.; validation, R.S., A.P.-D., W.V.-D. and J.F.-V.; formal analysis, W.V.-D.
and J.F.-V.; investigation, A.P.-D.; resources, A.P.-D.; data curation, R.S. and A.P.-D.; writing—original
draft preparation, A.P.-D., W.V.-D., R.S. and J.F.-V.; writing—review and editing, W.V.-D., A.P.-D. and
J.F.-V.; visualization, R.S. and W.V.-D.; supervision, A.P.-D.; project administration, A.P.-D.; funding
acquisition, A.P.-D. All authors have read and agreed to the published version of the manuscript.
Funding: The APC was funded by Universidad Tecnológica de Bolivar.
Data Availability Statement: Data are contained within the present article.
Acknowledgments: This work was supported by the Universidad Tecnológica de Bolívar under the
EOLITO research group.
Conflicts of Interest: The authors declare no conflicts of interest.
References
1. Garg, A. A New Approach to Optimizing Fired Heaters. 2010. Available online: https://oaktrust.library.tamu.edu/handle/1969
.1/94037 (accessed on 19 September 2024).
2. Garg, A. Optimize fired heater operations to save money. Hydrocarb. Process. 1997, 76, 97–112.
3. Litvinenko, V. The role of hydrocarbons in the global energy agenda: The focus on liquefied natural gas. Resources 2020, 9, 59.
[CrossRef]
4. Filimonova, I.V.; Cherepanova, D.; Provornaya, I.; Kozhevin, V.; Nemov, V. The dependence of sustainable economic growth on
the complex of factors in hydrocarbons-exporting countries. Energy Rep. 2020, 6, 68–73. [CrossRef]
5. Tcvetkov, P.; Cherepovitsyn, A.; Makhovikov, A. Economic assessment of heat and power generation from small-scale liquefied
natural gas in Russia. Energy Rep. 2020, 6, 391–402. [CrossRef]
6. Garba, M.D.; Usman, M.; Khan, S.; Shehzad, F.; Galadima, A.; Ehsan, M.F.; Ghanem, A.S.; Humayun, M. CO2 towards fuels:
A review of catalytic conversion of carbon dioxide to hydrocarbons. J. Environ. Chem. Eng. 2021, 9, 104756. [CrossRef]
7. Miller, S.A.; Wilkinson, J.D.; Lynch, J.T.; Hudson, H.M.; Cuellar, K.T.; Johnke, A.F.; Lewis, W.L. Hydrocarbon Gas Processing.
U.S. Patent 10,227,273, 12 March 2019.
8. Qamar, R.A.; Mushtaq, A.; Ulla, A.; Ali, Z.U. Simulation of Liquefied Petroleum Gas Recovery from Off Gases in a Fuel Oil
Refinery. J. Adv. Res. Fluid Mech. Therm. Sci. 2020, 73, 109–130. [CrossRef]
9. Eshaghi, S.; Hamrang, F. An innovative techno-economic analysis for the selection of an integrated ejector system in the flare gas
recovery of a refinery plant. Energy 2021, 228, 120594. [CrossRef]
10. Cala, O.M.; Stand, L.M.; Kafarov, V.; Rueda, J.S. Efecto de la composición del gas de refinería sobre las características del proceso
de combustión. Rev. Ing. Univ. Medellín 2013, 12, 101–111. [CrossRef]
11. Arefin, M.A.; Nabi, M.N.; Akram, M.W.; Islam, M.T.; Chowdhury, M.W. A review on liquefied natural gas as fuels for dual fuel
engines: Opportunities, challenges and responses. Energies 2020, 13, 6127. [CrossRef]
12. Amell, A.; Barraza, L.; Gómez, E. Tecnología de la Combustión de Gases y Quemadores Atmosféricos de Premezcla. Línea,
Disponible. 1996. Available online: http://es.scribd.com/doc/73707395/3-Quemadores-Atmosfericos-1 (accessed on 19 Septem-
ber 2024).
13. Murshed, M.; Alam, R.; Ansarin, A. The environmental Kuznets curve hypothesis for Bangladesh: The importance of natural gas,
liquefied petroleum gas, and hydropower consumption. Environ. Sci. Pollut. Res. 2021, 28, 17208–17227. [CrossRef]
14. Kim, S.Y.; Costa, A.L.; Bagajewicz, M.J. New robust approach for the globally optimal design of fired heaters. Chem. Eng. Res.
Des. 2023, 197, 434–448. [CrossRef]
15. Ghorashi, S.A.; Khandelwal, B. Toward the ultra-clean and highly efficient biomass-fired heaters. A review. Renew. Energy 2023,
205, 631–647. [CrossRef]
Energies 2024, 17, 6190 14 of 14
16. Qasim, F.; Lee, D.H.; Won, J.; Ha, J.K.; Park, S.J. Development of Advanced Advisory System for Anomalies (AAA) to Predict
and Detect the Abnormal Operation in Fired Heaters for Real Time Process Safety and Optimization. Energies 2021, 14, 7183.
[CrossRef]
17. Thorat, S.; McQueen, G.; Luzunaris, P.T. The role of optimal design and application of heat tracing systems to improve the energy
conservation in petrochemical facilities. IEEE Trans. Ind. Appl. 2013, 50, 163–173. [CrossRef]
18. Masoumi, M.E.; Izakmehri, Z. Improving of refinery furnaces efficiency using mathematical modeling. Int. J. Model. Optim. 2011,
1, 74. [CrossRef]
19. EU. A Clean Planet for All, a European Long-Term Strategic Vision for a Prosperous, Modern, Competitive and Climate Neutral Economy;
Depth Analysis in Support of the Commission Communication Com; European Commission: Brussels, Belgium, 2018; Volume 773.
20. Wildy, F.; Instruments, A.P. Fired heater optimization. In Technical Sales Support Manager; AMETEK Process Instruments:
Pittsburgh, PA, USA, 2000.
21. Rico, J.C.S.; Sánchez, Y.A.C. Análisis teórico de la combustión en quemadores de gas natural. Sci. Tech. 2005, 3, 139–143.
22. Bouras, F.; Khaldi, F. Optimization of Combustion Efficiency Using a Fuel Composition Based on CH4 and/or H2. Russ. J. Appl.
Chem. 2020, 93, 1954–1959. [CrossRef]
23. Horbaj, P. Simulation method for optimization of a mixture of fuel gases. Gas Waerme Int. 2000, 49, 250–252.
24. Kazi, S.R.; Sundar, K.; Srinivasan, S.; Zlotnik, A. Modeling and optimization of steady flow of natural gas and hydrogen mixtures
in pipeline networks. Int. J. Hydrogen Energy 2024, 54, 14–24. [CrossRef]
25. Simsek, S.; Uslu, S.; Simsek, H.; Uslu, G. Improving the combustion process by determining the optimum percentage of liquefied
petroleum gas (LPG) via response surface methodology (RSM) in a spark ignition (SI) engine running on gasoline-LPG blends.
Fuel Process. Technol. 2021, 221, 106947. [CrossRef]
26. Derikvand, H.; Dehaj, M.S.; Taghavifar, H. The effect of different sampling method integrated in NSGA II optimization on
performance and emission of diesel/hydrogen dual-fuel CI engine. Appl. Soft Comput. 2022, 128, 109434. [CrossRef]
27. Saifullin, E.; Larionov, V.; Busarov, A.; Busarov, V. Optimization of hydrocarbon fuels combustion variable composition in thermal
power plants. J. Phys. Conf. Ser. 2016, 669, 012037. [CrossRef]
28. Patrón, G.D.; Ricardez-Sandoval, L. An integrated real-time optimization, control, and estimation scheme for post-combustion
CO2 capture. Appl. Energy 2022, 308, 118302. [CrossRef]
29. Cote Florez, M.S. Revision de Alternativas Operativas Para el Cumplimiento de Especificaciones de Calidad de Gas en un Campo
Productor. Ph.D. Thesis, Universidad Industrial de Santander, Escuela De Ingeniería de Petróleos, Bucaramanga, Colombia, 2018.
30. Arroyo, H.; Nichiren, S.; Alonso, M.; Delfin, F. Calidad y Medición del Gas Natural. 2006. Available online: https://es.scribd.
com/document/462864821/Calidad-y-medicion-del-gas-natural (accessed on 19 September 2024).
31. Chomiak, J.; Longwell, J.; Sarofim, A. Combustion of low calorific value gases; problems and prospects. Prog. Energy Combust.
Sci. 1989, 15, 109–129. [CrossRef]
32. Pinzón Vargas, B.L.; Plazas Puentes, M.I. Evaluación Técnico-Financiera de las Tecnologías de Construcción Modular para
la Refinación de Petróleo Crudo en el Proyecto RefiBoyacá. Bachelor’s Thesis, Fundación Universidad de América, Bogotá,
Colombia, 2018.
33. Qin, C.; Li, J.; Yang, C.; Ai, B.; Zhou, Y. Comparative Study of Parameter Extraction from a Solar Cell or a Photovoltaic Module by
Combining Metaheuristic Algorithms with Different Simulation Current Calculation Methods. Energies 2024, 17, 2284. [CrossRef]
34. Olikara, C.; Borman, G.L. A Computer Program for Calculating Properties of Equilibrium Combustion Products with Some Applications to
IC Engines; Technical Report; SAE Technical Paper; SAE International: Warrendale, PA, USA, 1975.
35. Lide, D.R. CRC Handbook of Chemistry and Physics; CRC Press: Boca Raton, FL, USA, 2004; Volume 85.
36. Turns, S.R. Introdução à Combustão-: Conceitos e Aplicações; AMGH Editora: Porto Alegre, Brazil, 2013.
37. Worters, M.; Millard, D.; Hunter, G.; Helling, C.; Woitke, P. Comparison Catalogue of Gas-Equilibrium Constants, Kp. 2017.
Available online: https://research-repository.st-andrews.ac.uk/handle/10023/12242 (accessed on 19 September 2024).
38. Klotz, I.M. Chemical Thermodynamics; WA Benjamin Inc.: New York, NY, USA, 1964.
39. Yıldız, M. Chemical equilibrium based combustion model to evaluate the effects of H2 addition to biogases with different CO2
contents. Int. J. Hydrogen Energy 2024, 52, 1334–1344. [CrossRef]
40. Way, R. Methods for determination of composition and thermodynamic properties of combustion products for internal combustion
engine calculations. Proc. Inst. Mech. Eng. 1976, 190, 687–697. [CrossRef]
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual
author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to
people or property resulting from any ideas, methods, instructions or products referred to in the content.