# Density calculation with cubic equation of state

I want to calculate the density of CO2 using a cubic equation of state in order to plot the P-v diagram. From theory I know that I get one real root when the substance is in single-phase region and two when it is in two-phase region. In the latter case the smaller volume is associated to liquid while the other to vapor. I'm using Matlab to solve the cubic eos, incrementing the pressure and the temperature at fixed intervals. How can I merge the two solutions to represent the actual Pv plot?

First off, in the two-phase region you have three roots not two.

Smallest volume = liquid

Largest volume = gas

Middle volume = an unstable state you can't really reach

To determine whether a given temperature and pressure gives the liquid or the gas you must compare the Gibbs free energy, $G$, between the liquid and gas roots as determined from solving the cubic equation.

The change in Gibbs free energy from liquid to gas is given by

$\Delta G = \int_{liquid}^{gas}{V dP}$

where $P$ and $V$ are pressure and volume along a path of constant temperature. Note that the need for constant temperature; you can't just hold pressure constant and change the volume because that is not constant temperature.

The proper way to do it is to render

$dP=(dP/dV) dV$

where the derivative $dP/dV$ is obtained as a function of temperature and volume by differentiating the equation of state at constant temperature. Then you have

$\Delta G = \int_{liquid}^{gas}{V(dP/dV) dV}$

with $dP/dV$ known as a function of volume at constant temperature. You carry out the integral, which requires knowing how to integrate rational functions, and check whether $\Delta G$ is positive, negative, or zero:

$\Delta G>0$ means liquid has lower free energy, so liquid is more stable. You plot the point for the liquid phase root.

$\Delta G<0$ means gas has lower free energy, so gas is more stable. You plot the point for the gas phase root.

$\Delta G=0$ means both phases have the same free energy, you hit the equilibrium between liquid and gas. You plot both points and connect them with a line segment corresponding to both phases in equilibrium. Jackpot!