Clausius–Clapeyron equation: shape of phase diagrams makes no sense

Question

I am trying to model the melting point of a substance at varying pressures (ranging from very small to very very large). All I am trying to do is make an equation that relates melting temperature to pressure, so $T(P)$ is some function. To do this, I am trying to use the Clausius–Clapeyron equation (CC), which states that

$$\frac{\mathrm dP}{\mathrm dT} = \frac{L}{TΔV}.$$

In other words, the slope of the equilibrium line on the phase diagram should decrease as temperature is increased.

However, this is not the case; the curve of the equilibrium line is exponential and the slope $\mathrm dP/\mathrm dT$ increases as $T$ increases. Integrating CC we arrive at a logarithmic function, which again is not what empirical measurements reflect.

As I see it, the empirical results and the equation that is supposed to describe them are mutually exclusive. There is no way to arrive at an exponential curve from a slope that varies with $1/x.$ The CC equation and phase diagrams cannot be both be true at once and it is driving me mad.

Why is this the case? Is the CC equation valid at all because it seems to be totally false? What function do I use to model melting points at different temperatures?

The results that are so stupefying are these:

The shape of the curve is exponential. But, the supposed derivative is $1/T$, in which case the slope of each curve (red and blue here) should flatten as $T$ increases, but it steepens. Also, integrating that supposed derivative gives us $\ln (T)$ which is definitely not the shape of the phase diagram. This discrepancy is true for both the liquid/solid and liquid/gas curves. I hope this clarifies the question!

Answer 1

This is a summary of the equations to use to calculate phase transitions.

The Clapeyron equation $\displaystyle p_2-p_1=\frac{\Delta H}{\Delta V}\ln\left( \frac{T_2}{T_1} \right)$ is used for a solid-liquid transition. The changes in enthalpy and volume relate therefore to changes occurring in fusion.

The Clausius-Clapeyron equation describes solid-vapour and liquid-vapour changes because the final volume is far greater than the initial one, and is $\displaystyle \frac{dp}{dT}=p\frac{\Delta H}{RT^2}$ where $\Delta H$ the enthalpy change at the liquid–vapour or sublimation transition. Integrating this last equation from pressure $p_1 \to p_2$ and temperature $T_1 \to T_2$ gives $\displaystyle \ln\left(\frac{p_2}{p_1} \right) = -\frac{\Delta_{vap}H}{R}\left( \frac{1}{T_2}-\frac{1}{T_1} \right) $.

The change in the volume during fusion is $\displaystyle \Delta_{fus}V = m\left(\frac{1}{d_l}-\frac{1}{d_s} \right)$ where $m$ is the molar mass and $d_l$ and $d_s$ the densities of the liquid and solid. The pressure variation for the solid to liquid (melting or fusion) change is

$$\displaystyle p_2=p_1+\frac{\Delta_{fus}H}{\Delta_{fus}V}\ln\left(\frac{T_2}{T_1} \right)$$

and for evaporation and sublimation

$$\displaystyle p_2=p_1\exp\left( -\frac{\Delta_{vap}H}{R}\left(\frac{1}{T_2}-\frac{1}{T_1} \right) \right)$$

with the appropriate $\Delta H$. This is $\Delta_{vap}H$ for evaporation and $\Delta_{fus}H + \Delta_{vap}H$ for sublimation. Sublimation is treated as two steps merged into one; melting and instantaneous evaporation.

The $p$ vs $T$ plot for benzene is shown in the figure. Notice how this is different to how these are generally shown. This is sometimes due to the fact that log pressure is sometimes plotted but this is not always shown on the figure. Notice also how the solid liquid line is effectively vertical.

The data used is

R = 8.314 #( J/mol/K)
dens_sol = 981.0 #( kg/m^3)
dens_liq = 879.0 #( kg/m^3)
mol_mass = 78.0/1000.0 #( kg/mol)
DH_vap = 30.8*1000 #( J/mol)
DH_fus = 10.6*1000 #( J/mol)
p3 = 36.0/760*101325 #( triple point pressure Pa)
T3 = 5.5 + 273.16 #( triple point temperature K)
DV_fus = mol_mass*(1/dens_liq-1/dens_sol) # delta volume fusion

Python/numpy functions for the pressure are:
p_liq_vap= lambda T: p3*np.exp( (DH_vap/R)*(1/T3-1/T))
p_sol_vap= lambda T: p3*np.exp( ((DH_fus + DH_vap)/R ) * (1/T3-1/T) )
p_sol_liq= lambda T: p3 + DH_fus/DV_fus*(np.log(T) - np.log(T3))