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

# Find the z factor using Hall-Yarborough correlation
# This procedure follows Tarek Ahmed method. Pag 156 of his book
pres.pr <- 0.5
temp.pr <- 1.3    # pseudo-reduce temperature
t <- 1 / temp.pr  # reciprocal

yk0 <- 0.0125 * pres.pr * t * exp(-1.2 * (1 - t)^2)

    X1 <- -0.06125 * pres.pr * t * exp(-1.2 * (1 - t)^2)
    X2 <- t * (14.76 - 9.76 * t + 4.58 * t^2)
    X3 <- t * (90.7 - 242.2 * t + 42.4 * t^2)
    X4 <- 2.18 + 2.82 * t

    yk <- yk0

while (TRUE) {       
    y <- yk

    fy <- X1 + (y + y^2 + y^3 - y^4) / (1 - y)^3  - X2 * y^2 + X3 * y^X4  

    fydot <- (1 + 4 * y + 4 * y^2 - 4 * y^3 + y^4 ) / (1 - y)^4 - 
        2 * X2 * y + X3 * X4 * y^(X4-1)

    if (is.na(fy) || is.na(fydot)) {
        # yk <- 0.0125 * pres.pr * t * exp(-1.2 * (1 - t)^2)
        yk <- yk + 0.001
    } else {

        yk1 <- yk - fy / fydot

        error <- abs(yk - yk1)

        if (error < 1E-13) break

        cat(yk, yk1, fydot, error, "\n")

        yk <- yk1
    }
}

cat("\n")   
 y <- yk
 z <- -X1 / y
 print(z)    

Values from Standing and Katz chart

(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

tol <- 0.0000001
tol <- 1E-13

f <- function(y) {
    - A * pres.pr + (y + y^2 + y^3 - y^4) / (1 - y)^3  - B * y^2 + C * y^D    
}

fdot <- function(y) {
    (1 + 4 * y + 4 * y^2 - 4 * y^3 + y^4 ) / (1 - y)^4 - 2 * B * y + C * D * y^(D-1)
}

pres.pr <- 0.5
temp.pr <- 1.3

t <- 1 / temp.pr
yk <- 0.0125 * pres.pr * t * exp(-1.2 * (1 - t)^2)
delta <- 1
i <- 1
while (delta >= tol) {

    A <- 0.06125 * t * exp(-1.2 * (1 - t)^2)
    B <- t * (14.76 - 9.76 * t + 4.58 * t^2)
    C <- t * (90.7 - 242.2 * t + 42.4 * t^2)
    D <- 2.18 + 2.82 * t

    fyk <- f(yk)
    if (abs(fyk) < tol) break

    yk1 <- yk - f(yk) / fdot(yk)
    delta <- abs(yk - yk1)
    cat(sprintf("%3d %10f %10f %10f \n", i, delta, yk, fyk))
    yk <- yk1
    i <- i + 1
}    
cat("\n")   
 y <- yk
 z <- A * pres.pr / y
 print(z)

Values from Standing and Katz chart

(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.