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. tr = temperature / tc tr # Reduced temperature # T and Tc have the units of °R. : def reduced_temperature (temperature, tc) return 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. pr = pressure / pc pr # Reduced pressure # P and Pc have the units of psia : def reduced_pressure (pressure, pc) return The reduced pressure is defined as its actual pressure divided by its critical pressure. Step 3: Calculate the Redlich-Kwong constants as: a = * (pr / tr ** ) a b = * (pr / tr) b # Redlich-Kwong constant A : def redlich_kwong_constant_a (pr, tr) 0.42748 2.5 return # Redlich-Kwong constant B : def redlich_kwong_constant_b (pr, tr) 0.08664 return 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: alpha = ( / ) * ( * (a - b - b ** ) - ) alpha beta = ( / ) * ( + ( * (a - b - b ** )) - ( * a * b)) beta # Cubic constant α : def cubic_constant_alpha (a, b) 1 3 3 2 1 return # Cubic constant β : def cubic_constant_beta (a, b) 1 27 -2 9 2 27 return Step 5: Calculate the Discriminant as: d = (beta ** / ) + (alpha ** / ) d # Discriminant : def discriminant (alpha, beta) 2 4 3 27 return Step 6: Calculate the Solution constants as: a_star = np.cbrt((-beta / ) + np.sqrt(d)) a_star b_star = np.cbrt((-beta / ) - np.sqrt(d)) b_star theta = math.acos(-sign(beta) * (math.sqrt((beta ** / ) / (-alpha ** / )))) theta # Solution Constant A* : def solution_constant_a_star (beta, d) 2 return # Solution Constant B* : def solution_constant_b_star (beta, d) 2 return # Solution Constant Theta : def solution_constant_theta (beta, alpha) 2 4 3 27 return SIGN(β) is calculated as below: beta > : sign_beta = beta < : sign_beta = : sign_beta = sign_beta # SIGN(β) # If β > 0, SIGN(β) = 1 # If β = 0, SIGN(β) = 0 # If β < 0, SIGN(β) = -1 : def sign (beta) if 0 1 elif 0 -1 else 0 return Step 7: Calculate the Trial Roots as: where the cosine term is in radians, not degrees. z1 = a_star + b_star + / z1 zi = (-( / ) * (a_star + b_star)) + (( / ) * i) zi zi1 = ( * math.sqrt(- alpha / ) * math.cos((theta / ) + (i * ((math.pi * ) / )))) + ( / ) zi1 # Trial Root Z1 : def trial_root_z1 (a_star, b_star) 1 3 return # Trial Root Zi # for i = 2, 3 : def trial_root_zi (a_star, b_star, i) 1 2 1 3 return # Trial Root Zi+1 # for i = 1, 2, 3 : def trial_root_zi1 (alpha, theta, i) 2 3 3 2 3 1 3 return Critical Constants for Selected Gases are given below: gas_critical_constants = { : [ , , ], : [ , , ], : [ , , ], : [ , , ], : [ , , ], : [ , , ], : [ , , ] } gas_critical_constants.get(gas_name) : 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 """ "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 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: values = table_gas_critical_constants(gas) tc = values[ ] pc = values[ ] 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) d < : theta = solution_constant_theta(beta, alpha) z1 = trial_root_zi1(alpha, theta, ) z2 = trial_root_zi1(alpha, theta, ) z3 = trial_root_zi1(alpha, theta, ) : a_star = solution_constant_a_star(beta, d) b_star = solution_constant_b_star(beta, d) d == : z1 = trial_root_z1(a_star, b_star) z2 = trial_root_zi(a_star, b_star, ) z3 = trial_root_zi(a_star, b_star, ) : z = trial_root_z1(a_star, b_star) z z = max(z1, z2, z3) z : 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 """ 0 1 if 0 1 2 3 else if 0 2 3 else return return Cover image from Compressibility Factor, ScienceDirect Hope this helps. Happy Coding!