Problem 14.03 Reversible and Isothermal Compression of Liquid Water

A flow of 2000 kg/h liquid water at 25 °C and 1 bar is pumped to a pressure of 100 bar. The pump is “cooled”, so the process is reversible and isothermal. What is the duty of the pump shaft, and the energy that must be removed from the water being compressed?


We can use the high-accuracy IAPWS-95 implementation of the properties of water to easily and extremely accurately calculate these values.

from scipy.constants import bar, hour
import numpy as np
from thermo import FlashPureVLS, IAPWS95Liquid, IAPWS95Gas, iapws_constants, iapws_correlations
from scipy.integrate import quad
from chemicals import property_molar_to_mass
import numpy as np

T1 = T2 = 25 + 273.15
P1 = 1*bar
P2 = 100*bar

liquid = IAPWS95Liquid(T=T1, P=P1, zs=[1])
gas = IAPWS95Gas(T=T1, P=P1, zs=[1])
flasher = FlashPureVLS(iapws_constants, iapws_correlations, gas, [liquid], [])

mass_flow = 2000/hour
mole_flow = property_molar_to_mass(mass_flow, MW=iapws_constants.MWs[0])

entry = flasher.flash(T=T1, P=P1)
leaving = flasher.flash(T=T2, P=P2)

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

cooling_duty = shaft_duty - (leaving.H() - entry.H())*mole_flow

print(f'The shaft power is {shaft_duty:.8f} W')
print(f'The cooling duty is {cooling_duty:.4f} W')
The shaft power is 5504.05633851 W
The cooling duty is 431.1770 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 = 30 # 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)
    compressing_duty += (new_compressed_state.H() - old_state.H())*mole_flow

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

    old_state = new_cooled_state

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