Problem 14.08 Required Compressor Power for Isothermal and Adiabatic Compression of a Gas Mixture (CO2, O2) Using the Ideal Gas Law

A stream of 1000 mol/hour CO2 and 1000 mol/hour O2 is compressed from 290 K and 1 bar to 5 bar. Calculate the compression power for both adiabatic compression, and isothermal compression. The compression is reversible (assumed) in each case - no efficiencies are necessary.

Solution

This is a straightforward calculation. Using Thermo, working with complicated mixtures can be about as easy as pure components - if binary interaction parameters are zero. In this case, we try to load a parameter from a sample ChemSep database, but no values are available.

The values in that database are just a sample - it is entirely the user’s responsibility to provide the correct data to Thermo. If garbage is put in, garbage will come out!

The problem says to use the ideal-gas law, so we can do that too and see how the answers compare.

from scipy.constants import hour
T1 = 290
P1 = 1e5
P2 = 5e5
flow = 2000/hour # mol/s

from thermo import ChemicalConstantsPackage, PRMIX, IGMIX, FlashVL, CEOSLiquid, CEOSGas
from thermo.interaction_parameters import IPDB

constants, correlations = ChemicalConstantsPackage.from_IDs(['CO2', 'O2'])
kijs = IPDB.get_ip_asymmetric_matrix('ChemSep PR', constants.CASs, 'kij')
print(f'The PR kij matrix is {kijs}')

eos_kwargs = dict(Tcs=constants.Tcs, Pcs=constants.Pcs, omegas=constants.omegas,
                 kijs=kijs)
liquid = CEOSLiquid(PRMIX, HeatCapacityGases=correlations.HeatCapacityGases, eos_kwargs=eos_kwargs)
gas = CEOSGas(PRMIX, HeatCapacityGases=correlations.HeatCapacityGases, eos_kwargs=eos_kwargs)
flasher = FlashVL(constants, correlations, liquid=liquid, gas=gas)
zs = [.5, .5]

liquid = CEOSLiquid(IGMIX, HeatCapacityGases=correlations.HeatCapacityGases, eos_kwargs=eos_kwargs)
gas = CEOSGas(IGMIX, HeatCapacityGases=correlations.HeatCapacityGases, eos_kwargs=eos_kwargs)
flasher_ideal = FlashVL(constants, correlations, liquid=liquid, gas=gas)

The PR kij matrix is [[0.0, 0], [0, 0.0]]

Adiabatic compression

# Solve with Peng-Robinson
state_1 = flasher.flash(T=T1, P=P1, zs=zs)
state_2 = flasher.flash(S=state_1.S(), P=P2, zs=zs)
shaft_duty = (state_2.H() - state_1.H())*flow

print(f'The shaft power with Peng-Robinson is {shaft_duty:.4f} W')


state_1 = flasher_ideal.flash(T=T1, P=P1, zs=zs)
state_2 = flasher_ideal.flash(S=state_1.S(), P=P2, zs=zs)
shaft_duty = (state_2.H() - state_1.H())*flow
print(f'The shaft power with ideal-gas is {shaft_duty:.4f} W')

The shaft power with Peng-Robinson is 2632.7895 W
The shaft power with ideal-gas is 2639.9248 W

Isothermal Compression

This problem is more interesting, because there is the cooling duty as well as the compressing duty.

From theory, in an ideal gas, the cooling duty will be exactly equal to the compressing duty.

For a real-gas, it will be different as enthalpy is pressure-dependent.

In both cases, the evaluation of the following integral is required.

$$\text{duty} = \text{flow} \int_{P1}^{P2} V \partial P$$

from scipy.integrate import quad

def to_int(P, flasher):
    state = flasher.flash(T=T1, P=P, zs=zs)
    return state.V()
shaft_duty = cooling_duty = quad(to_int, P1, P2, args=(flasher_ideal,))[0]*flow

print(f'The shaft power with ideal-gas is {shaft_duty:.4f} W')
print(f'The cooling duty with ideal-gas is {cooling_duty:.4f} W')

entry = flasher.flash(T=T1, P=P1, zs=zs)
exit = flasher.flash(T=T1, P=P2, zs=zs)

shaft_duty = quad(to_int, P1, P2, args=(flasher,))[0]*flow
cooling_duty = shaft_duty - (exit.H() - entry.H())*flow

print(f'The shaft power with Peng-Robinson is {shaft_duty:.8f} W')
print(f'The cooling duty with Peng-Robinson is {cooling_duty:.8f} W')

The shaft power with ideal-gas is 2155.9263 W
The cooling duty with ideal-gas is 2155.9263 W
The shaft power with Peng-Robinson is 2139.44610002 W
The cooling duty with Peng-Robinson is 2192.57596810 W

The above shows the numerical integral calculation. That is the correct formulation.

However, it can be a little unintuitive. We can contrast this with another calculation - a series of tiny isentropic compression, then cooling steps.

cooling_duty = 0
compressing_duty = 0
increments = 3 # Number of increments
dP = (P2 - P1)/increments
old_state = entry
for i in range(increments):
    P = P1+(i+1)*dP

    # Compress another increment of pressure
    new_compressed_state = flasher.flash(S=old_state.S(), P=P, zs=zs)
    compressing_duty += (new_compressed_state.H() - old_state.H())*flow

    # Cool back to T1 at new pressure
    new_cooled_state = flasher.flash(T=T1, P=P, zs=zs)
    cooling_duty += (new_compressed_state.H() - new_cooled_state.H())*flow

    old_state = new_cooled_state

print(f'The shaft power is {compressing_duty:.8f} W')
print(f'The cooling duty is {cooling_duty:.8f} W')

The shaft power is 2322.61227046 W
The cooling duty is 2375.74213854 W