How to Create a Compressibility Factor Calculator in Python

In thermodynamics, the compressibility factor (Z), also known as the compression factor or the gas deviation factor, is a correction factor which describes the deviation of a real gas from ideal gas behavior behavior.

The calculations below can be used to find the inlet or the outlet compressibility factor by using the values for temperature and pressure. T and P are the temperature and pressure in units of °R and psia, respectively. The value calculated for Z should not be greater than zero. The compressibility factor has no units of measure.

Below are the steps and the equations to calculate the compressibility factor for a given gas based on the Redlich-Kwong equation of state:

Step 1: Calculate the Reduced temperature as:

where T and Tc have the units of °R.

# Reduced temperature
# T and Tc have the units of °R.
def reduced_temperature(temperature, tc):
    tr = temperature / tc
    return tr

The reduced temperature of a fluid is its actual temperature, divided by its critical temperature. A critical point (or critical state) is the end point of a phase equilibrium curve. At the critical point, defined by a critical temperature Tc and a critical pressure Pc, phase boundaries vanish.

Step 2: Calculate the Reduced pressure as:

where P and Pc have the units of psia.

# Reduced pressure
# P and Pc have the units of psia
def reduced_pressure(pressure, pc):
    pr = pressure / pc
    return pr

The reduced pressure is defined as its actual pressure divided by its critical pressure.

Step 3: Calculate the Redlich-Kwong constants as:

# Redlich-Kwong constant A
def redlich_kwong_constant_a(pr, tr):
    a = 0.42748 * (pr / tr ** 2.5)
    return a


# Redlich-Kwong constant B
def redlich_kwong_constant_b(pr, tr):
    b = 0.08664 * (pr / tr)
    return b

In thermodynamics, the Redlich–Kwong equation of state is an empirical, algebraic equation that relates temperature, pressure, and volume of gases. The Redlich–Kwong equation is adequate for calculation of gas phase properties when the ratio of the pressure to the critical pressure (reduced pressure) is less than about one-half of the ratio of the temperature to the critical temperature (reduced temperature).

Step 4: Calculate the Cubic constants as:

# Cubic constant α
def cubic_constant_alpha(a, b):
    alpha = (1 / 3) * (3 * (a - b - b ** 2) - 1)
    return alpha


# Cubic constant β
def cubic_constant_beta(a, b):
    beta = (1 / 27) * (-2 + (9 * (a - b - b ** 2)) - (27 * a * b))
    return beta

Step 5: Calculate the Discriminant as:

# Discriminant
def discriminant(alpha, beta):
    d = (beta ** 2 / 4) + (alpha ** 3 / 27)
    return d

Step 6: Calculate the Solution constants as:

# Solution Constant A*
def solution_constant_a_star(beta, d):
    a_star = np.cbrt((-beta / 2) + np.sqrt(d))
    return a_star


# Solution Constant B*
def solution_constant_b_star(beta, d):
    b_star = np.cbrt((-beta / 2) - np.sqrt(d))
    return b_star


# Solution Constant Theta
def solution_constant_theta(beta, alpha):
    theta = math.acos(-sign(beta) * (math.sqrt((beta ** 2 / 4) / (-alpha ** 3 / 27))))
    return theta

SIGN(β) is calculated as below:

# SIGN(β)
# If β > 0, SIGN(β) = 1
# If β = 0, SIGN(β) = 0
# If β < 0, SIGN(β) = -1
def sign(beta):
    if beta > 0:
        sign_beta = 1
    elif beta < 0:
        sign_beta = -1
    else:
        sign_beta = 0
    return sign_beta

Step 7: Calculate the Trial Roots as:

where the cosine term is in radians, not degrees.

# Trial Root Z1
def trial_root_z1(a_star, b_star):
    z1 = a_star + b_star + 1 / 3
    return z1


# Trial Root Zi
# for i = 2, 3
def trial_root_zi(a_star, b_star, i):
    zi = (-(1 / 2) * (a_star + b_star)) + ((1 / 3) * i)
    return zi


# Trial Root Zi+1
# for i = 1, 2, 3
def trial_root_zi1(alpha, theta, i):
    zi1 = (2 * math.sqrt(- alpha / 3) * math.cos((theta / 3) + (i * ((math.pi * 2) / 3)))) + (1 / 3)
    return zi1

Critical Constants for Selected Gases are given below:

def table_gas_critical_constants(gas_name):
    """
    Retrieves gas critical constants table data
    :param gas_name: Gas name
    :type gas_name: str
    :return: Gas critical constants data
    :rtype: list
    """
    gas_critical_constants = {
        "Air": [238.56, 549.11, 28.965],
        "Carbon Dioxide": [547.60, 1070.600, 44.011],
        "Hydrogen": [59.82, 190.82, 2.02],
        "Methane": [343.90, 673.100, 16.04],
        "Nitrogen": [227.25, 492.420, 28.13],
        "Propane": [666.00, 618.700, 44.09],
        "Typical Natural Gas": [360.00, 777.373, 17.185]
    }
    return gas_critical_constants.get(gas_name)

If Discriminant (D) is lesser than zero, Solution constant (θ) and Trial roots (Zi+1 for i = 1, 2, 3) will be calculated. The max of the Trial roots will be calculated as compressibility factor (Z).

If Discriminant (D) is greater than or equal to zero, Solution constants (A* and B*) will be calculated. On that if D is equal to zero, then Trial roots (Z1 and Zi for i = 2, 3) will be calculated. The max of the Trial roots will be calculated as compressibility factor (Z). But if D is greater than zero, then Trial root (Z1) will be calculated and it will give the compressibility factor (Z).

Below is the program using all the above equations and steps for the compressibility factor calculation for a given gas based on the Redlich-Kwong equation of state:

def compressibility_factor_calc(gas, temperature, pressure):
    """
    Calculates compressibility factor
    :param gas: Gas name
    :type gas: str
    :param temperature: Temperature
    :type temperature: float
    :param pressure: Pressure
    :type pressure: float
    :return: Compressibility factor(z)
    :rtype: float
    """
    values = table_gas_critical_constants(gas)
    tc = values[0]
    pc = values[1]
    tr = reduced_temperature(temperature, tc)
    pr = reduced_pressure(pressure, pc)
    a = redlich_kwong_constant_a(pr, tr)
    b = redlich_kwong_constant_b(pr, tr)
    alpha = cubic_constant_alpha(a, b)
    beta = cubic_constant_beta(a, b)
    d = discriminant(alpha, beta)
    if d < 0:
        theta = solution_constant_theta(beta, alpha)
        z1 = trial_root_zi1(alpha, theta, 1)
        z2 = trial_root_zi1(alpha, theta, 2)
        z3 = trial_root_zi1(alpha, theta, 3)
    else:
        a_star = solution_constant_a_star(beta, d)
        b_star = solution_constant_b_star(beta, d)
        if d == 0:
            z1 = trial_root_z1(a_star, b_star)
            z2 = trial_root_zi(a_star, b_star, 2)
            z3 = trial_root_zi(a_star, b_star, 3)
        else:
            z = trial_root_z1(a_star, b_star)
            return z
    z = max(z1, z2, z3)
    return z

Cover image from Compressibility Factor, ScienceDirect

Hope this helps. Happy Coding!