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