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.
# 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) }
# 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")
# 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"))
# 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
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.