In f0nzie/zFactor: Calculate the Compressibility Factor 'z' for Hydrocarbon Gases

knitr::opts_chunk$set(echo=T, comment=NA, error=T, warning=F, message = F, fig.align = 'center')

The correlation was put in comptational form using two sources: the book "Equations of State and PVT Analysis" by Tarek Ahmed and the paper "New explicit correlation for the compressibility factor of natural gas" by Kareem, Iwalewa and Al-Marhoun.

We found an error in the constants in one of the sources and use the other one to reconstruct the right equation. Errors are found by comparing the results against the original data points in the Standing-Katz chart (1500 points). Usually errors are evident because the equation will not converge, give infinite values or results are just way off the real z values.

DPR correlation as a function

# Dranchuk, Purvis, and Robinson (1974) developed a correlation based on the 
# Benedict-Webb-Rubin type of equation of state. Fitting the equation to 1500 
# data points from the Standing and Katz Z-factor chart optimized the eight 
# coefficients of the proposed equations. From Traker Ahmed's book.

z.DranchukPurvisRobinson <- function(pres.pr, temp.pr, tolerance = 1E-13, 
                                verbose = FALSE) {
    F <- function(rhor) {
        1 + T1 * rhor + T2 * rhor^2 + T3 * rhor^5 + 
            (T4 * rhor^2 * (1 + A8 * rhor^2) * exp(-A8 * rhor^2)) - T5 / rhor
    }

    Fprime <- function(rhor) {
        T1 + 2 * T2 * rhor + 5 * T3 * rhor^4 + 
            2 * T4 * rhor * exp(-A8 * rhor^2) * 
            ((1 + 2 * A8 * rhor^2) - A8 * rhor^2 * (1 + A8 * rhor^2)) +
             T5 / rhor^2
    }

    A1 <- 0.31506237  
    A2 <- -1.0467099  
    A3 <- -0.57832720 
    A4 <-  0.53530771
    A5 <- -0.61232032
    A6 <- -0.10488813
    A7 <- 0.68157001
    A8 <- 0.68446549

    T1 <- A1 + A2 / temp.pr + A3 / temp.pr^3
    T2 <- A4 + A5 / temp.pr
    T3 <- A5 * A6 / temp.pr
    T4 <- A7 / temp.pr^3
    T5 <- 0.27 * pres.pr / temp.pr

    rhork0 <- 0.27 * pres.pr / temp.pr    # initial guess
    rhork <- rhork0
    i <- 1
    while (TRUE) {
        if (abs(F(rhork)) < tolerance)  break    
        rhork1 <- rhork - F(rhork) / Fprime(rhork)
        delta <- abs(rhork - rhork1)
        if (delta < tolerance) break
        if (verbose) 
            cat(sprintf("%3d %11f %11f %11f \n", i, rhork, rhork1, delta))
        rhork <- rhork1
        i <- i + 1
    }
    z <- 0.27 * pres.pr / (rhork * temp.pr)
    return(z)
}

Plotting only Dranchuk-Purvis-Robinson solution

# test DAK with  using the values from paper

ppr <- c(0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5)
tpr <- c(1.3, 1.5, 1.7, 2)

tbl <- sapply(ppr, function(x) 
    sapply(tpr, function(y) z.DranchukPurvisRobinson(pres.pr = x, temp.pr = y)))

rownames(tbl) <- tpr
colnames(tbl) <- ppr
print(tbl)

library(ggplot2)
plot(x = ppr, y = tbl[1,], type = "l", main = "z @ Tpr = 1.3", ylab = "z")
plot(x = ppr, y = tbl[2,], type = "l", main = "z @ Tpr = 1.5", ylab = "z")
plot(x = ppr, y = tbl[3,], type = "l", main = "z @ Tpr = 1.7", ylab = "z")
plot(x = ppr, y = tbl[4,], type = "l", main = "z @ Tpr = 2.0", ylab = "z")

Comparative plotting against Standing-Katz chart values

# DPR vs z read from SK chart
ppr <- c(0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5)
tpr <- c(1.3, 1.5, 1.7, 2)

# values from Standing-Katz chart
z <- c(
    c(0.92, 0.76, 0.64, 0.63, 0.68, 0.76, 0.84),
    c(0.94, 0.86, 0.79, 0.77, 0.79, 0.84, 0.89),
    c(0.97, 0.92, 0.87, 0.86, 0.865, 0.895, 0.94),
    c(0.985, 0.957, 0.941, 0.938, 0.945, 0.97, 1.01)
    )

mx <- matrix(z, nrow = 4, ncol = 7, byrow = TRUE)
rownames(mx) <- tpr
colnames(mx) <- ppr
print(mx)

library(ggplot2)
plot(x = ppr, y = mx[1,], type = "l", main = "z SK @ Tpr = 1.3", ylab = "z")
lines(x = ppr, y = tbl[1,], col = "red")
legend("topright", c("F(Y)", "chart"), lty = c(1,1), col = c("black", "red"))
plot(x = ppr, y = mx[2,], type = "l", main = "z SK @ Tpr = 1.5", ylab = "z")
lines(x = ppr, y = tbl[2,], col = "red")
legend("topright", c("F(Y)", "chart"), lty = c(1,1), col = c("black", "red"))
plot(x = ppr, y = mx[3,], type = "l", main = "z SK @ Tpr = 1.7", ylab = "z")
lines(x = ppr, y = tbl[3,], col = "red")
legend("topright", c("F(Y)", "chart"), lty = c(1,1), col = c("black", "red"))
plot(x = ppr, y = mx[4,], type = "l", main = "z SK @ Tpr = 2.0", ylab = "z")
lines(x = ppr, y = tbl[4,], col = "red")
legend("topright", c("F(Y)", "chart"), lty = c(1,1), col = c("black", "red"))

Raw form of the computational form of the DPR correlation

# Dranchuk, Purvis, and Robinson (1974) developed a correlation based on the 
# Benedict-Webb-Rubin type of equation of state. Fitting the equation to 1500 
# data points from the Standing and Katz Z-factor chart optimized the eight 
# coefficients of the proposed equations. From Traker Ahmed's book.

pres.pr <- 6.5
temp.pr <- 1.3
tolerance = 1E-13
verbose <- TRUE

A1 <- 0.31506237  
A2 <- -1.0467099  
A3 <- -0.57832720 
A4 <-  0.53530771
A5 <- -0.61232032
A6 <- -0.10488813
A7 <- 0.68157001
A8 <- 0.68446549

T1 <- A1 + A2 / temp.pr + A3 / temp.pr^3
T2 <- A4 + A5 / temp.pr
T3 <- A5 * A6 / temp.pr
T4 <- A7 / temp.pr^3
T5 <- 0.27 * pres.pr / temp.pr

F <- function(rhor) {
    1 + T1 * rhor + T2 * rhor^2 + T3 * rhor^5 + 
        (T4 * rhor^2 * (1 + A8 * rhor^2) * exp(-A8 * rhor^2)) - T5 / rhor
}

Fprime <- function(rhor) {
    T1 + 2 * T2 * rhor + 5 * T3 * rhor^4 + 
        2 * T4 * rhor * exp(-A8 * rhor^2) * 
        ((1 + 2 * A8 * rhor^2) - A8 * rhor^2 * (1 + A8 * rhor^2)) +
         T5 / rhor^2
}


rhork0 <- 0.27 * pres.pr / temp.pr
rhork <- rhork0
i <- 1
while (TRUE) {
    if (abs(F(rhork)) < tolerance)  break    
    rhork1 <- rhork - F(rhork) / Fprime(rhork)
    delta <- abs(rhork - rhork1)
    if (delta < tolerance) break
    if (verbose) 
        cat(sprintf("%3d %11f %11f %11f \n", i, rhork, rhork1, delta))
    rhork <- rhork1
    i <- i + 1
}
z <- 0.27 * pres.pr / (rhork * temp.pr)
z

0.92
0.75
0.64
0.63
0.69
0.765
0.85

(0.5, 1.3) = 0.92
(1.5, 1.3) = 0.76
(2.5, 1.3) = 0.64
(3.5, 1.3) = 0.63
(4.5, 1.3) = 0.68
(5.5, 1.3) = 0.76
(6.5, 1.3) = 0.84

f0nzie/zFactor documentation built on Aug. 2, 2019, 1:42 a.m.