How to plot Compressibility factor Z vs Pressure P using ParametricPlot?

I am trying to get my plot to look like this:

but I cannot seem to get it right when I try to plot each curve as a parametric equation. Here is what I have been trying:

 ParametricPlot[{{v/(v - 0.0043067) - 0.23026/(8.3144598*(180)*v), (
8.3144598*(180))/(v - 0.0043067) - 0.23026/v^2}, {v/(
v - 0.0043067) - 0.23026/(8.3144598*(189)*v), (8.3144598*(189))/(
v - 0.0043067) - 0.23026/v^2}, {v/(v - 0.0043067) - 0.23026/(
8.3144598*(190)*v), (8.3144598*(190))/(v - 0.0043067) - 0.23026/
v^2}, {v/(v - 0.0043067) - 0.23026/(8.3144598*(200)*v), (
8.3144598*(200))/(v - 0.0043067) - 0.23026/v^2}, {v/(
v - 0.0043067) - 0.23026/(8.3144598*(250)*v), (8.3144598*(250))/(
v - 0.0043067) - 0.23026/v^2}}, {v, 1, 2}]


In all honesty, I don't know what the range on v should be, since it is the molar volume. On each equation, the temperature in Kelvin is already in.

• For compressibility I used the van der waals equation and solved for Z in terms of molar volume (Z=v/(v-b)-a/RTv) and for pressure I also used the van der waals equation and solved for pressure in terms of molar volume (P=RT/(v-b)-a/v^2). I tried to use both equations in the parametric plot but I might be going about it all wrong. – Jim Jan 8 '17 at 22:10
• Also this plot is for methane at 180 K, 189 K, 190 K, 200 K, and 250 K. The van der Waals parameters for methane are a=2.3026 and b=0.043067 – Jim Jan 8 '17 at 23:36

The plot you were trying to reproduce is a plot of the reduced van der Waals equation; that is, it is using the reduced pressure $P_r=\frac{P}{P_c}$ and reduced temperature $T_r=\frac{T}{T_c}$ as the independent variables. Thus you need to use the reduced form

$$Z=\frac38\left(\frac{8V_r}{3V_r-1}-\frac3{V_r T_r}\right)$$

where $V_r$ is obtained from solving the equation

$$P_r=\frac{8T_r}{3V_r-1}-\frac3{V_r^2}$$

Note that the constants $a$ and $b$ in the conventional presentation are absent, having been absorbed into the critical constants.

Thus,

With[{tlist = {1., 1.1, 1.2, 1.3, 1.5, 2.}},
Plot[Table[With[{v = Root[-3 + 9 #1 - (p + 8 t) #1^2 + 3 p #1^3 &, 1]},
3/8 ((8 v)/(3 v - 1) - 3/(t v))], {t, tlist}] // Evaluate,
{p, 0, 7}, Axes -> None, Frame -> True,
FrameLabel -> {"Reduced Pressure (\!$$\*SubscriptBox[\(p$$, $$r$$]\))",
"Compressibility Factor (Z)"}, GridLines -> Automatic,
PlotLegends -> tlist, PlotRange -> All]]


Here is a 3D version as a bonus:

Plot3D[With[{v = Root[-3 + 9 #1 - (p + 8 t) #1^2 + 3 p #1^3 &, 1]},
3/8 ((8 v)/(3 v - 1) - 3/(t v))], {p, 0, 7}, {t, 0.9, 2.1},
AxesLabel -> {"(\!$$\*SubscriptBox[\(p$$, $$r$$]\))",
"(\!$$\*SubscriptBox[\(t$$, $$r$$]\))", "Z"},
BoundaryStyle -> None, Mesh -> {{1., 1.1, 1.2, 1.3, 1.5, 2.}},
MeshFunctions -> (#2 &), PlotStyle -> Opacity[1/2, ColorData[97, 1]]]


• Awesome feedback! Is there any way to get a plot of Z vs P without using the reduced van der waals equation? Or in other words, plot Z against P for methane (a=2.3026 dm^6 bar mol^-2 and b=0.043067 dm^3 mol^-1) using only the van der waals equation by solving for Z and P? – Jim Jan 8 '17 at 22:34
• Next time, you could be more descriptive with your questions; I shouldn't have had to guess that you were using vdW. Anyway: an explicit expression for molar volume would go something like With[{r = QuantityMagnitude[UnitConvert[Quantity[1, "MolarGasConstant"]]]/100, a = 2.3026, b = 0.043067}, Root[-a b + a # - (b p + r t) #^2 + p #^3 &, 1]], which you can now plug into vdW. – J. M. is away Jan 8 '17 at 22:40
• Ok, sorry for the ambiguity, I am still new to all of this. So using With[{r = QuantityMagnitude[UnitConvert[Quantity[1, "MolarGasConstant"]]]/ 100, a = 2.3026, b = 0.043067}, Root[-a b + a # - (b p + r t) #^2 + p #^3 &, 1]] how would one end up with the plot of Z vs P? – Jim Jan 8 '17 at 23:34
• Did you notice that the function supplied is now only dependent on p and t? Replace the molar volume in your compressibility factor expression with this, and you now have $Z$ as a function of $P$ and $T$. – J. M. is away Jan 8 '17 at 23:44
• Makes sense. The only problem I still have is translating that into mathematica language. Could I just plot that function on its own as a list of functions with different values of T and have it depend on P? – Jim Jan 8 '17 at 23:53

Well it's messy but I figured it out. I used solve to solve the van der Waals equation for z and then plotted it.

With[{tlist = {180, 189, 190, 200, 250}},
Plot[Table[
With[{r =
QuantityMagnitude[
UnitConvert[Quantity[1, "MolarGasConstant"]]]/100,
a = 2.3026, b = 0.043067}, (-b p r^2 t^2 - r^3 t^3)/(
3 r^3 t^3) - (2^(
1/3) (3 a p r^4 t^4 - (-b p r^2 t^2 -
r^3 t^3)^2))/(3 r^3 t^3 (18 a b p^2 r^6 t^6 +
2 b^3 p^3 r^6 t^6 - 9 a p r^7 t^7 + 6 b^2 p^2 r^7 t^7 +
6 b p r^8 t^8 +
2 r^9 t^9 + \[Sqrt]((18 a b p^2 r^6 t^6 +
2 b^3 p^3 r^6 t^6 - 9 a p r^7 t^7 +
6 b^2 p^2 r^7 t^7 + 6 b p r^8 t^8 + 2 r^9 t^9)^2 +
4 (3 a p r^4 t^4 - (-b p r^2 t^2 - r^3 t^3)^2)^3))^(
1/3)) + 1/(
3 2^(1/3)
r^3 t^3) (18 a b p^2 r^6 t^6 + 2 b^3 p^3 r^6 t^6 -
9 a p r^7 t^7 + 6 b^2 p^2 r^7 t^7 + 6 b p r^8 t^8 +
2 r^9 t^9 + \[Sqrt]((18 a b p^2 r^6 t^6 +
2 b^3 p^3 r^6 t^6 - 9 a p r^7 t^7 + 6 b^2 p^2 r^7 t^7 +
6 b p r^8 t^8 + 2 r^9 t^9)^2 +
4 (3 a p r^4 t^4 - (-b p r^2 t^2 - r^3 t^3)^2)^3))^(
1/3)], {t, tlist}] // Evaluate, {p, 0, 700}, Axes -> None,
Frame -> True, PlotLegends -> tlist, PlotRange -> All]]


• Try this: With[{r = QuantityMagnitude[UnitConvert[Quantity[1, "MolarGasConstant"]]]/100, a = 2.3026, b = 0.043067, tlist = {180, 189, 190, 200, 250}}, Plot[Table[With[{v = Root[-a b + a # - (b p + r t) #^2 + p #^3 &, 1]}, v/(v - b) - a/(r t v)], {t, tlist}] // Evaluate, {p, 0, 700}, Axes -> None, Frame -> True, PlotLegends -> tlist, PlotRange -> All]]. – J. M. is away Jan 9 '17 at 8:52