Second Law efficiency of a furnace#

A conventional furnace provides space heating to a building, as shown here:

Conventional furnace

The furnace combusts methane with 200% excess dry air (\(ea = 2\)). The air and methane enter the combustor at \(T_{\text{in}}\) = -10°C and \(P_{\text{in}}\) = 1 atm. The furnace heats the building (\(Q_{\text{out}}\)), which is maintained at \(T_{\text{bldg}}\) = 22°C. The combustion products leave at \(T_{\text{out}}\) = 150°C and \(P_{\text{out}}\) = 1 atm.

Assumptions: the only combustion products are carbon dioxide, water, nitrogen, and oxygen (for excess air in the reactants). The gases follow the ideal gas law.

Problems:

  • Determine the First Law efficiency of the furnace based on the higher and lower heating values of the fuel.

  • Determine the Second Law efficiency of the furnace.

This problem also demonstrates how to use Calculate enthalpies using mixture object and how to determine the Exergy of a fuel.

First, we should import the necessary modules and enter the given information.

import numpy as np
import cantera as ct
from scipy.optimize import root

from pint import UnitRegistry
ureg = UnitRegistry()
Q_ = ureg.Quantity

# for convenience:
def to_si(quant):
    '''Converts a Pint Quantity to magnitude at base SI units.
    '''
    return quant.to_base_units().magnitude

# given information
temperature_in = Q_(-10, 'degC').to('K')  # also dead state temperature
pressure_in = Q_(1, 'atm')  # also dead state pressure

temperature_building = Q_(22, 'degC').to('K')

temperature_out = Q_(150, 'degC').to('K')
pressure_out = Q_(1, 'atm')

excess_air = 2.0

First Law efficiency#

To calculate the First Law efficiency of the furnace, we need to determine its heat output, and also the lower and higher heating values of the fuel.

The heat output comes from an energy balance on the system:

\[ H_{\text{reac}} = Q_{\text{out}} + H_{\text{prod}} \;, \]

where \(H_{\text{reac}}\) and \(H_{\text{prod}}\) are the enthalpies of the reactants and products, respectively. To determine these enthalpy values, we need to determine the relative amounts of each component in the reactants and products.

Reaction stoichiometry#

The stoichiometric reaction of methane with air is

\[ \text{CH}_4 + a_s \left( 0.21 \text{O}_2 + 0.79 \text{N}2 \right) \rightarrow b_s \text{CO}_2 + c_s \text{H}_2 \text{O} + d_s \text{N}_2 \]

and we can solve for the unknown coefficients based on elemental balances of carbon, hydrogen, oxygen, and nitrogen:

\[\begin{split} 1 = b_s \\ 4 = 2 c_s \\ a_s (0.21) 2 = 2 b_s + c_s \\ a_s (0.79) 2 = 2 d_s \end{split}\]
coeffs_stoich = {}

coeffs_stoich['CO2'] = 1
coeffs_stoich['H2O'] = 4 / 2
coeffs_stoich['air'] = (
    2 * coeffs_stoich['CO2'] + coeffs_stoich['H2O']
    ) / (0.21 * 2)
coeffs_stoich['N2'] = coeffs_stoich['air'] * 0.79 * 2 / 2

Then, the actual reaction with excess air is

\[ \text{CH}_4 + a \left( 0.21 \text{O}_2 + 0.79 \text{N}2 \right) \rightarrow b \text{CO}_2 + c \text{H}_2 \text{O} + d \text{N}_2 + e \text{O}_2 \]

where \(a\) is based on the amount of excess air, with respect to the stoichiometric amount of air:

\[ a = (1 + ea) a_s \]

We can also find the actual coefficients by balancing each of the elements:

\[\begin{split} 1 = b \\ 4 = 2 c \\ a (0.21) 2 = 2 b + c + 2e \\ a (0.79) 2 = 2 d \end{split}\]
coeffs = {}

coeffs['air'] = (1 + excess_air) * coeffs_stoich['air']

coeffs['CO2'] = 1
coeffs['H2O'] = 4 / 2
coeffs['O2'] = (
    coeffs['air'] * 0.21 * 2 - 2 * coeffs['CO2'] - coeffs['H2O']
    ) / 2
coeffs['N2'] = coeffs['air'] * 0.79 * 2 / 2

print('Moles of reaction participants, per mole of fuel:')
for species, value in coeffs.items():
    print(f'{species}: {value: .2f}')

Moles of reaction participants, per mole of fuel:
air:  28.57
CO2:  1.00
H2O:  2.00
O2:  4.00
N2:  22.57

Next, to carry out the energy balance, we need to find the phase of the water in the combustion products. To do this, we can determine the dew point temperature, based on the saturation temperature of water at its partial pressure, and compare with the temperature of the products.

The partial pressure of water vapor in the products is

\[ P_v = y_v P_{\text{out}} \;, \]

where \(y_v\) is the mole fraction of water vapor (assuming no condensation):

\[ y_v = \frac{c}{b + c + d + e} \;. \]
# mole fraction of water vapor (no condensation)
mole_fraction_water = coeffs['H2O'] / (
    coeffs['CO2'] + coeffs['H2O'] + coeffs['N2'] + coeffs['O2']
    )
pressure_water = mole_fraction_water * pressure_out

water = ct.Water()

# saturated water vapor
water.PQ = to_si(pressure_water), 1.0
temperature_dewpoint = Q_(water.T, 'K')

print(f'Dew point temperature: {temperature_dewpoint: .2f}')
print(f'Products temperature: {temperature_out: .2f}')

print('Products temperature above dew point temperature: '
      f'{temperature_out > temperature_dewpoint}'
      )

Dew point temperature: 311.80 kelvin
Products temperature: 423.15 kelvin
Products temperature above dew point temperature: True

Since the temperature of the exhaust is higher than the dew point temperature, the water in the products is entirely in the vapor phase (i.e., superheated vapor).

Calculate enthalpies#

Now, we can calculate the enthalpies of the reactants and products (per mole of fuel), based on the molar specific enthalpy of the components in each:

enthalpies_reactants = {}
enthalpies_products = {}

# enthalpy of reactants
reactants = ['CH4', 'O2', 'N2']
gas = ct.Solution('gri30.yaml')

for reactant in reactants:
    gas.TPX = (
        to_si(temperature_in), to_si(pressure_in), 
        f'{reactant}:1.0'
        )
    enthalpies_reactants[reactant] = Q_(gas.enthalpy_mole, 'J/kmol')

# per kmol of fuel
enthalpy_reactants = (
    enthalpies_reactants['CH4'] + coeffs['air'] * (
        0.21 * enthalpies_reactants['O2'] + 
        0.79 * enthalpies_reactants['N2']
        )
    )

products = ['CO2', 'H2O', 'N2', 'O2']
gas = ct.Solution('gri30.yaml')

for product in products:
    gas.TPX = (
        to_si(temperature_out), to_si(pressure_out), 
        f'{product}:1.0'
        )
    enthalpies_products[product] = Q_(gas.enthalpy_mole, 'J/kmol')

# per kmol of fuel
enthalpy_products = sum([
    coeffs[p] * enthalpies_products[p] for p in products
    ])

We can now calculate the heat output of the furnace, per kmol of fuel:

heat = enthalpy_reactants - enthalpy_products
print(f'Heat output: {heat.to("MJ/kmol"): .2f}')

Heat output: 661.43 megajoule / kilomole

Calculate enthalpies using mixture object#

Rather than manually calculating the enthalpy of each component of the reactants and products separately and then combining, we could also calculate the enthalpy of the reactants and products as mixtures, using a Cantera Solution object for each.

We can specify the mixture state with the temperature, pressure, and calculated numbers of moles for each component: gas.TPX. Then, we can get the specific enthalpy of the mixture with gas.enthalpy_mole.

However, to get the extensive enthalpy of the mixture, we need to multiply this by the number of moles of the mixture. The resulting total enthalpy will be per 1 kmol of fuel, since that was the basis for our reaction stoichiometry calculations:

gas_reactants = ct.Solution('gri30.yaml')
reactant_string = (
    'CH4:1.0, '
    f"O2:{0.21 * coeffs['air']: .3f}, "
    f"N2:{0.79 * coeffs['air']: .3f}"
    )
gas.TPX = to_si(temperature_in), to_si(pressure_in), reactant_string

# per 1 kmol of fuel
moles_reactants = 1.0 + coeffs['air']
enthalpy_reactants = Q_(gas.enthalpy_mole, 'J/kmol') * moles_reactants

gas_products = ct.Solution('gri30.yaml')
product_string = ', '.join([f'{sp}:{coeffs[sp]: .3f}' for sp in products])
gas.TPX = to_si(temperature_out), to_si(pressure_out), product_string

# per 1 kmol of fuel
moles_products = sum([coeffs[p] for p in products])
enthalpy_products = Q_(gas.enthalpy_mole, 'J/kmol') * moles_products

heat = enthalpy_reactants - enthalpy_products
print(f'Heat output: {heat.to("MJ/kmol"): .2f}')

Heat output: 661.44 megajoule / kilomole

This route is a bit simpler, and we get effectively the same value.

Calculate efficiency#

Then, the First Law efficiency is based on either the lower or higher heating value of the fuel:

\[\begin{split} \eta_{\text{LHV}} = \frac{Q_{\text{out}}}{\text{LHV}} \\ \eta_{\text{HHV}} = \frac{Q_{\text{out}}}{\text{HHV}} \end{split}\]

which we can find for methane based on tabulated information: LHV = 50,032 kJ/kg and HHV = 55,516 kJ/kg. (Or, we can calculate these values, as shown later: Heating values.)

molecular_weight_methane = Q_(16.04, 'kg/kmol')

heating_value_lower = Q_(50032, 'kJ/kg') * molecular_weight_methane
heating_value_higher = Q_(55516, 'kJ/kg') * molecular_weight_methane

efficiency_lower = 100 * heat / heating_value_lower
efficiency_higher = 100 * heat / heating_value_higher

print('First Law efficiencies:')
print(f' lower heating value: {to_si(efficiency_lower): .2f} %')
print(f' higher heating value: {to_si(efficiency_higher): .2f} %')

First Law efficiencies:
 lower heating value:  82.42 %
 higher heating value:  74.28 %

The furnace appears very efficient from a first-law perspective.

Second Law efficiency#

The Second Law efficiency is based on the exergy transfer of heat out of the system and the exergy transfer in from the fuel:

\[ \eta_2 = \frac{X_{Q_{\text{out}}}}{X_{\text{fuel}}} \;, \]

where the exergy flow due to heat transfer is

\[ X_{Q_{\text{out}}} = Q_{\text{out}} \left(1 - \frac{T_0}{T_{\text{bldg}}} \right) \;. \]

The dead state temperature is the temperature of the outdoor air.

exergy_heat = heat * (1 - (temperature_in / temperature_building))
print(f'Exergy of heat transfer: {exergy_heat.to("kJ/kmol"): .2f}')

Exergy of heat transfer: 71712.94 kilojoule / kilomole

Exergy of a fuel#

The exergy of the fuel is the maximum possible work obtainable by the fuel:

\[ X_{\text{fuel}} = H_{\text{reac}} - H_{\text{prod}} - T_0 \left( S_{\text{reac}} - S_{\text{prod}} \right) \;, \]

where \(H_{\text{reac}}\) and \(H_{\text{prod}}\) are the enthalpy of the reactants and products (per mole of fuel), and \(S_{\text{reac}}\) and \(S_{\text{prod}}\) are the entropy of the reactants and products (per mole of fuel); both are evaluated at the dead state temperature (\(T_0\)) and pressure. Furthermore, the entropy of each reactant/product is evaluated at its partial pressure.

If the dead state temperature is the same as the reference temperature, \(T_0 = T_{\text{ref}}\) = 25°C, then \((H_{\text{reac}} - H_{\text{prod}})\) will be bounded by the lower and higher heating values of the fuel, with the exact value depending on the dead state definition. (The term involving entropy is much smaller.)

Thus, the fuel exergy can be found in two ways:

  1. Approximate using the heating value, or \(x_{\text{fuel}} \approx \text{HV}\).

  2. Calculate based on the dead state conditions, or

So, we can calculate the Second Law efficiency using the lower heating value (to be conservative):

efficiency_second = 100 * to_si(exergy_heat / heating_value_lower)
print(f'Second Law efficiency (LHV): {efficiency_second: .2f}%')

Second Law efficiency (LHV):  8.94%

We can also calculate the exergy of the fuel, based on the specific dead state conditions here. To do this, we need to calculate the enthalpy and entropy of the reactants and products involved, per mole of fuel. These are

\[\begin{split} H = \sum_{i=1}^C n_i \overline{h}_{i, 0} (T_0) \quad \text{and} \\ S = \sum_{i=1}^C n_i \overline{s}_{i, 0} (T_0, P_i) \;, \\ \end{split}\]

where \(C\) is the number of components in the reactants or products, \(n_i\) is the stoichiometric coefficient for component \(i\) (i.e., the number of moles per mole of fuel), \(\overline{h}_{i, 0}\) is the molar specific enthalpy for component \(i\), \(\overline{s}_{i, 0}\) is the molar specific entropy for component \(i\), and \(P_i\) is the partial pressure of component \(i\): \(P_i = y_i P\).

To calculate the exergy of the fuel, we only really need to consider the components that actively participate in the chemical reaction, meaning the fuel, oxygen, carbon dioxide, and water. All other reactants and products, including the excess air, are at the same state before and after the reaction, and so do not contribute. However, for simplicity, we can include all reactants and products.

gas = ct.Solution('gri30.yaml')

enthalpies = {}
entropies = {}

moles_reactants = 1.0 + coeffs['air']

partial_pressure = (1.0 / moles_reactants) * pressure_in
gas.TPX = to_si(temperature_in), to_si(partial_pressure), 'CH4:1.0'
enthalpies['CH4'] = Q_(gas.enthalpy_mole, 'J/kmol')
entropies['CH4'] = Q_(gas.entropy_mole, 'J/(K*kmol)')

partial_pressure = (coeffs['air'] / moles_reactants) * pressure_in
gas.TPX = to_si(temperature_in), to_si(partial_pressure), 'O2:0.21, N2:0.79'
enthalpies['air'] = Q_(gas.enthalpy_mole, 'J/kmol')
entropies['air'] = Q_(gas.entropy_mole, 'J/(K*kmol)')

# per kmol of fuel
enthalpy_reactants = (
    enthalpies['CH4'] + coeffs['air'] * enthalpies['air']
    )
entropy_reactants = (
    entropies['CH4'] + coeffs['air'] * entropies['air']
    )

moles_products = sum([coeffs[c] for c in products])
for product in products:
    partial_pressure = (coeffs[product] / moles_products) * pressure_in
    gas.TPX = (
        to_si(temperature_in), to_si(partial_pressure), 
        f'{product}:1.0'
        )
    enthalpies[product] = Q_(gas.enthalpy_mole, 'J/kmol')
    entropies[product] = Q_(gas.entropy_mole, 'J/(K*kmol)')

# per kmol of fuel
enthalpy_products = sum([
    coeffs[p] * enthalpies[p] for p in products
    ])
entropy_products = sum([
    coeffs[p] * entropies[p] for p in products
    ])

exergy_fuel = (
    enthalpy_reactants - enthalpy_products - 
    temperature_in * (entropy_reactants - entropy_products)
    )

efficiency_second = 100 * to_si(exergy_heat / exergy_fuel)
print(f'Second Law efficiency: {efficiency_second: .2f}%')

Second Law efficiency:  8.86%

We get a similar, but slightly lower, Second Law efficiency when calculating it based on the more-correct exergy of the fuel.

Either way, we see that the efficiency of this furnace is quite low from a Second Law perspective, due to the significant exergy leaving with the exhaust.