knitr::opts_chunk$set(echo=T, comment=NA, error=T, warning=F, message = F, fig.align = 'center', results="hold")
Kenneth Hall and Lyman Yarborough used the hard-sphere equation as the basis for the equation of state. They tested the correlation with 12 reservoir gas reservoir systems up to Ppr
as high as 20.5. The Standing-Katz chart only extends to Ppr=15
. At that moment the Standing-Katz chart had 30 years of existance. See [@Hall1973].
z
at selected Ppr
and Tpr=2.0
Use the the corelation to calculate z
. From the Standing-Katz chart we obtain a digitized point at the given Tpr
and Ppr
.
# get a z value using HY library(zFactor) ppr <- 1.5 tpr <- 2.0 z.calc <- z.HallYarborough(pres.pr = ppr, temp.pr = tpr) # get a z value from the SK chart at the same Ppr and Tpr z.chart <- getStandingKatzMatrix(tpr_vector = tpr, pprRange = "lp")[1, as.character(ppr)] # calculate the APE (Average Percentage Error) ape <- abs((z.calc - z.chart) / z.chart) * 100 df <- as.data.frame(list(Ppr = ppr, z.calc =z.calc, z.chart = z.chart, ape=ape)) rownames(df) <- tpr df
z
at selected Ppr
and Tpr = 1.1
From the Standing-Katz chart we obtain a digitized point:
library(zFactor) ppr <- 1.5 tpr <- 1.1 z.calc <- z.HallYarborough(pres.pr = ppr, temp.pr = tpr) # From the Standing-Katz chart we obtain a digitized point: z.chart <- getStandingKatzMatrix(tpr_vector = tpr, pprRange = "lp")[1, as.character(ppr)] # calculate the APE ape <- abs((z.calc - z.chart) / z.chart) * 100 df <- as.data.frame(list(Ppr = ppr, z.calc =z.calc, z.chart = z.chart, ape=ape)) rownames(df) <- tpr df
We see here a noticeable difference between the values of
z
from theHY
correlation and the value read from the Standing-Katz chart. This is expected at these low pressures and temperatures. This area is a challenge for all of the correlations as well.
z
for combinations of Ppr
and Tpr
In this example we provide vectors instead of a single point.
With the same ppr
and tpr
vectors that we use for the correlation, we do the same for the Standing-Katz
chart. We want to compare both and find the absolute percentage error
.
library(zFactor) ppr <- c(0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5) tpr <- c(1.05, 1.1, 1.7, 2) # calculate using the correlation z.calc <- z.HallYarborough(ppr, tpr) # With the same ppr and tpr vector, we do the same for the Standing-Katz chart z.chart <- getStandingKatzMatrix(ppr_vector = ppr, tpr_vector = tpr) ape <- abs((z.calc - z.chart) / z.chart) * 100 # calculate the APE cat("z.correlation \n"); print(z.calc) cat("\n z.chart \n"); print(z.chart) cat("\n APE \n"); print(ape)
You can see errors of 22.48% and 12.09% in the isotherm
Tpr=1.05
atPpr=1.5 and 2.5, respectively. Other errors, greater than one, can also be found at the isotherm
Tpr=1.1`. Then, the rest of the curves are fine.
isotherms
We apply the function summary
over the transpose of the matrix. Applying the transpose will alow us to see the statistics at each of the isotherms.
sum_t_ape <- summary(t(ape)) sum_t_ape
We see that the errors in
z
are considerable with ar sum_t_ape[1,1]
% andr sum_t_ape[6,1]
% forTpr=1.05
, and ar sum_t_ape[1,2]
% andr sum_t_ape[6,2]
% forTpr=1.10
.
The Hall-Yarborough correlation shows a very high error at values of Tpr
lower or equal than 1.1
, being Tpr=1.05
the worst curve to calculate z values from. Keep that in mind. We will explore later a comparative tile chart where we confirm these early calculations.
Tpr
Let's see the numbers now for higher values of Tpr
at various values of Ppr
.
library(zFactor) # enter vectors for Tpr and Ppr tpr2 <- c(1.2, 1.3, 1.5, 2.0, 3.0) ppr2 <- c(0.5, 1.5, 2.5, 3.5, 4.5, 5.5) # get z values from the SK chart z.chart <- getStandingKatzMatrix(ppr_vector = ppr2, tpr_vector = tpr2, pprRange = "lp") # We do the same with the HY correlation: # calculate z values at lower values of Tpr z.calc <- z.HallYarborough(pres.pr = ppr2, temp.pr = tpr2) ape <- abs((z.calc - z.chart) / z.chart) * 100 # calculate the APE cat("z.correlation \n"); print(z.calc) cat("\n z.chart \n"); print(z.chart) cat("\n APE \n"); print(ape)
At
Tpr
above or equal to 1.2 theHY
correlation behaves very well.
isotherms
Now, let's apply the function summary
over the transpose of the matrix to observe the error of the correlation at each isotherm.
sum_t_ape <- summary(t(ape)) sum_t_ape
In all, the error of
z
at these isotherms is less than 0.82%. Pretty good.
SK
chart values vs HY
correlationNow, we will be plotting the difference between the z
values in the Standing-Katz and the z
values calculated by Hall-Yarborough. For the moment let's not pay attention to the numerical error but to the error bars
in the plot (the orange bars).
library(zFactor) library(tibble) library(ggplot2) tpr2 <- c(1.05, 1.1, 1.2, 1.3) ppr2 <- c(0.5, 1.0, 1.5, 2, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5) sk_corr_2 <- createTidyFromMatrix(ppr2, tpr2, correlation = "HY") as_tibble(sk_corr_2) p <- ggplot(sk_corr_2, aes(x=Ppr, y=z.calc, group=Tpr, color=Tpr)) + geom_line() + geom_point() + geom_errorbar(aes(ymin=z.calc-abs(dif), ymax=z.calc+abs(dif)), width=.4, position=position_dodge(0.05)) print(p)
As we were expecting, the errors can be found on Tpr=1.05 and Tpr=1.1
Tpr
curvesIn this last example, we compare the values of z
at all the isotherms.
We use the function getCurvesDigitized
in the zfactor
package to obtain all the isotherms or Tpr
curves in the Standing-Katz chart that have been digitized. The next function createTidyFromMatrix
calculates z
using the correlation and prepares a tidy dataset ready to plot.
library(ggplot2) library(tibble) # get all `lp` Tpr curves tpr_all <- getStandingKatzTpr(pprRange = "lp") ppr <- c(0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5) sk_corr_all <- createTidyFromMatrix(ppr, tpr_all, correlation = "HY") as_tibble(sk_corr_all) p <- ggplot(sk_corr_all, aes(x=Ppr, y=z.calc, group=Tpr, color=Tpr)) + geom_line() + geom_point() + geom_errorbar(aes(ymin=z.calc-dif, ymax=z.calc+dif), width=.4, position=position_dodge(0.05)) print(p)
As we saw before, we confirm that the greatest errors are localized in two of the
Tpr
curves: at 1.05 and 1.1.
What we will see here is the distribution of the Mean Average Percentage Error or MAPE accross all the Tpr
isotherms and Ppr
grids. Remember from the README
: red and yellow are bad.
# MSE: Mean Squared Error # RMSE: Root Mean Sqyared Error # RSS: residual sum of square # RMSLE: Root Mean Squared Logarithmic Error. Penalizes understimation. # MAPE: Mean Absolute Percentage Error = AARE # MPE: Mean Percentage error = ARE # MAE: Mean Absolute Error library(dplyr) grouped <- group_by(sk_corr_all, Tpr, Ppr) smry_tpr_ppr <- summarise(grouped, RMSE= sqrt(mean((z.chart-z.calc)^2)), MPE = sum((z.calc - z.chart) / z.chart) * 100 / n(), MAPE = sum(abs((z.calc - z.chart) / z.chart)) * 100 / n(), MSE = sum((z.calc - z.chart)^2) / n(), RSS = sum((z.calc - z.chart)^2), MAE = sum(abs(z.calc - z.chart)) / n(), RMLSE = sqrt(1/n()*sum((log(z.calc +1)-log(z.chart +1))^2)) ) ggplot(smry_tpr_ppr, aes(Ppr, Tpr)) + geom_tile(data=smry_tpr_ppr, aes(fill=MAPE), color="white") + scale_fill_gradient2(low="blue", high="red", mid="yellow", na.value = "pink", midpoint=12.5, limit=c(0, 25), name="MAPE") + theme(axis.text.x = element_text(angle=45, vjust=1, size=11, hjust=1)) + coord_equal() + ggtitle("Hall-Yarborough", subtitle = "HY")
Tpr
and Ppr
values that show more errorThis is just plotting the couple of isotherms where we see the largest errors.
library(dplyr) sk_corr_all %>% filter(Tpr %in% c("1.05", "1.1")) %>% ggplot(aes(x = z.chart, y=z.calc, group = Tpr, color = Tpr)) + geom_point(size = 3) + geom_line(aes(x = z.chart, y = z.chart), color = "black") + facet_grid(. ~ Tpr) + geom_errorbar(aes(ymin=z.calc-abs(dif), ymax=z.calc+abs(dif)), position=position_dodge(0.05))
Finally, the dataframe with the calculated errors between the z
from the correlation and the z
read from the chart:
as_tibble(smry_tpr_ppr)
$$RMSE = \sqrt {\frac {1}{n} \sum_{i=1}^{n} (y_i - \hat y_i)^2}$$
getStandingKatzMatrix
is equivalent to using the sapply
function with the internal function .z.HallYarborough
(the dot means it's internal), which we call adding the prefix zFactor:::
. That is, the package name and three colons.
# test HY with 1st-derivative using the values from the 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) hy <- sapply(ppr, function(x) sapply(tpr, function(y) zFactor:::.z.HallYarborough(pres.pr = x, temp.pr = y))) rownames(hy) <- tpr colnames(hy) <- ppr print(hy)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.