1
$\begingroup$

I am trying to get my plot to look like this: Z vs P

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Jim Jan 8 '17 at 22:10
  • $\begingroup$ 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 $\endgroup$ – Jim Jan 8 '17 at 23:36
9
$\begingroup$

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]]

compressibility factor for reduced van der Waals equation


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]]]

compressibility factor surface

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – Jim Jan 8 '17 at 22:34
  • 1
    $\begingroup$ 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. $\endgroup$ – J. M. is away Jan 8 '17 at 22:40
  • $\begingroup$ 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? $\endgroup$ – Jim Jan 8 '17 at 23:34
  • $\begingroup$ 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$. $\endgroup$ – J. M. is away Jan 8 '17 at 23:44
  • $\begingroup$ 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? $\endgroup$ – Jim Jan 8 '17 at 23:53
0
$\begingroup$

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]]    

Z vs P

$\endgroup$
  • 1
    $\begingroup$ 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]]. $\endgroup$ – J. M. is away Jan 9 '17 at 8:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.