peng_robinson.py 文件源码

def get_Z(self, T, p):

        kappa = 0.37464 + 1.54226*self.acentric - 0.26992*self.acentric**2

        Tr = T/self.T_crit
        alpha = (1 + kappa*(1 - Tr**0.5))**2

        A = alpha * self.a * p / self.R**2 / T**2
        B = self.b * p / self.R / T

        # Solve the cubic equation for compressibility factor z
        coeffs = [1, -(1-B), A-2*B-3*B**2, -(A*B-B**2-B**3)]
        #print(coeffs)
        roots = np.roots(coeffs)

        real_roots = roots[np.isreal(roots)].real
        valid_roots = real_roots[real_roots > p*self.b/self.R/T]

        return valid_roots