Incorrect calculation of real density for compressible medium

[dad616610]

Describe the bug
As of now the real density is calculated this way $\rho = \frac{P \cdot T_n}{P_n \cdot T \cdot Z}$, however the correct formula is this one: $\rho = \frac{P \cdot T_n \cdot \boldsymbol Z_n}{P_n \cdot T \cdot Z}$, where $Z_n$ is a compressibility factor taken at normal conditions ($P_n = 101325 \ Pa$ and $T = 273.15 \ K$)

Reasoning
Compressibility factor is $Z = \frac{P \cdot V}{n \cdot R \cdot T} \Rightarrow 1 = \frac{P \cdot V}{n \cdot R \cdot T \cdot Z}$, which should hold true for any $P$, $V$, $T$ and $Z$ of the same medium.

Let's consider 2 states: state 1 is at given conditions and state 2 is at normal conditions. Then:

$$\frac{P_1 \cdot V_1}{n \cdot R \cdot T_1 \cdot Z_1} = \frac{P_2 \cdot V_2}{n \cdot R \cdot T_2 \cdot Z_2} \xRightarrow[\text{cancelling out n, R and m}]{V=\frac{m}{\rho}} \frac{P_1}{T_1 \cdot Z_1 \cdot \rho_1} = \frac{P_2}{T_2 \cdot Z_2 \cdot \rho_2} \xRightarrow[]{} \rho_1 = \rho_2 \frac{P_1 \cdot T_2 \cdot Z_2}{P_2 \cdot T_1 \cdot Z_1} = \rho_n \frac{P \cdot T_n \cdot Z_n}{P_n \cdot T \cdot Z}$$

So basically, for state 2 we're using a real gas, instead of ideal gas, which $Z = 1$

Data-driven proof
Calculate real density for Methane

Using NIST webbook

1. Let $P_n = 101325\ Pa$, $T_n = 273.15\ K$
2. Obtain $V_n = 0.022361\ m^3/mol$ and $\rho_n = 44.722\ mol/m^3$ at given conditions via: webbook link
3. Calculate $Z_n = 0.997631$ using $Z = \frac{P \cdot V}{n \cdot R \cdot T}$
4. Let $P = 60 \cdot P_n \ Pa$, $T = 320\ K$
5. Obtain $V = 0.00040562\ m^3/mol$ and $\rho = 2465.3 mol/m^3$ at given conditions via: webbook link
6. Calculate $Z = 0.926831$
7. Calculate $\rho_{\text{wrong}} = \rho_n \frac{P \cdot T_n}{P_n \cdot T \cdot Z} = 2471.286570\ mol/m^3$
8. Calculate $\rho_{\text{right}} = \rho_n \frac{P \cdot T_n \cdot Z_n}{P_n \cdot T \cdot Z} = 2465.432281\ mol/m^3$

Copy-pastable python script

R = 8.31451
calc_z = lambda _p, _t, _v: _p * _v / (R * _t)

pn = 101325  # Pa
tn = 273.15  # K
vn = 0.022361 # m^3/mol, obtained through link[1]
rhon = 44.722 # mol/m^3, obtained through link[1]
zn = calc_z(pn, tn, vn)  # 0.997631

p = 60 * pn  # Pa
t = 320  # K
v = 0.00040562  # m^3/mol, obtained through link[2]
rho = 2465.3 # mol/m^3, obtained through link[2]
z = calc_z(p, t, v)  # 0.926831

wrong_rho = rhon * (p * tn) / (pn * t * z)  # 2471.286570 mol/m^3
right_rho = wrong_rho * zn  # 2465.432281 mol/m^3

rel_err = lambda orig, approx: abs(approx - orig) / orig
print(f"{rel_err(rho, wrong_rho):.5%}")  # 0.24283 %
print(f"{rel_err(rho, right_rho):.5%}") # 0.00537 %

Using CoolProp library

import CoolProp.CoolProp

# create equation of state for Methane
eos = CoolProp.AbstractState("HEOS", "Methane")

pn = 101325  # Pa
tn = 273.15  # K
eos.update(CoolProp.PT_INPUTS, pn, tn) # update internal state to be at given conditions
rhon = eos.rhomolar() # 44.721543 mol/m^3, obtained through link[1]
zn = eos.compressibility_factor()  # 0.997613

p = 60 * pn  # Pa
t = 320  # K
eos.update(CoolProp.PT_INPUTS, p, t) # update internal state to be at given conditions
rho = eos.rhomolar()  # 2465.335806 mol/m^3, obtained through link[2]
z = eos.compressibility_factor()  # 0.926841

wrong_rho = rhon * (p * tn) / (pn * t * z)  # 2471.235416 mol/m^3
right_rho = wrong_rho * zn  # 2465.335806 mol/m^3

rel_err = lambda orig, approx: abs(approx - orig) / orig
print(f"{rel_err(rho, wrong_rho):.5%}")  # 0.23930 %
print(f"{rel_err(rho, right_rho):.5%}") # 0.00000 %