Isothermes de van der Waals

Isothermes_vdWaals.svg

vanDerWaals.py

#!/usr/bin/env python
# Copyright: This document has been placed in the public domain.

__author__ = "Yannick Copin <[email protected]>"

import numpy as N
import scipy.optimize as O
import scipy.integrate as I
import matplotlib.pyplot as P

P.matplotlib.style.use('classic')

def isotherme(v, t=1):
    """Equation des isothermes p=P/Pc--v=V/Vc en fonction de t=T/Tc"""

    return 8*t/(3*v - 1) - 3/v**2

def dpdv(v, t=1):

    return -24*t/(3*v-1)**2 + 6/v**3

def spinodale(v):
    """Equation de la spinodale: dp/dv=0"""

    #return 2*(3*v - 1)/v**3 - 3/v**2
    return 3/v**2 - 2/v**3

def objfunc(p0, t, roots=False):

    # Solutions de p(v) = p0
    p0 = N.squeeze(p0)
    sols = N.roots([3*p0, -(p0+8*t), 9, -3])

    if not N.isreal(sols).all():
        sols = sols.real

    v1,v0,v2 = N.sort(sols)
    if roots:
        return v1,v2

    i1,err = I.quad(lambda v:isotherme(v,t=t) - p0, v1,v0)
    i2,err = I.quad(lambda v:isotherme(v,t=t) - p0, v0,v2)

    return i1 + i2

def pSat(t):
    """Calcule de la pression de vapeur saturante a t"""

    if t > 1:
        raise ValueError("La pression de vapeur saturante est definie pour t<1")

    # Intersection de la spinodale avec l'isotherme
    sols = N.roots([4*t, -9, 6, -1])
    tmp,v1,v2 = N.sort(sols)
    p1,p2 = isotherme(N.array([v1,v2]), t=t)
    # zero de objfunc, pour p1 < p0 < p2
    try:
        psat = O.brentq(objfunc, p1, p2, args=(t,))
    except ValueError:
        psat = N.nan

    return psat

def courbeSat(t=None):

    if t is None:
        t = N.linspace(0.85,1-1e-6)

    psats = N.array([ pSat(tt) for tt in t ])
    good = N.isfinite(psats)
    vsats = N.array([ objfunc(psat,tt, roots=True)
                      for psat,tt in zip(psats[good],t[good]) ])

    vs = N.concatenate((vsats[:,0],vsats[:,1][::-1]))
    ps = N.concatenate((psats[good],psats[good][::-1]))

    return vs,ps


#fig,ax = P.subplots(1,1)
fig = P.figure()
ax = fig.add_subplot(1,1,1)

ax.set_title("Isothermes du gaz de Van der Waals")
ax.set_xlabel("V/Vc")
ax.set_ylabel("P/Pc")

v = N.linspace(0.4,5,1000)
for t in N.linspace(0.8,1.2,21):
    ax.plot(v,isotherme(v,t), c='0.8', label="_")
ax.plot(v,isotherme(v,0.9), label="T/Tc=0.9")
ax.plot(v,isotherme(v,1), label="T/Tc=1")
ax.plot(v,isotherme(v,1.1), label="T/Tc=1.1")
ax.plot([1],[1],'go', label='_')
ax.annotate('Point critique', (1,1), (5,5), textcoords='offset points')

ax.plot(v,spinodale(v), 'k--', label='Spinodale')
vsats,psats = courbeSat()
ax.plot(vsats,psats, 'k-', label='Courbe de saturation')

psat = pSat(0.9)
v1,v2 = objfunc(psat,0.9,roots=True)
ax.plot([v1,v2],[psat]*2, 'bo-', label='_')
vs = N.linspace(v1,v2,100)
ax.fill_between(vs,vs*0+psat,isotherme(vs,t=0.9),
                facecolor='b', edgecolor='none', alpha=0.2)

ax.legend()
#ax.semilogx()
#xticks = [0.4,0.6,0.8,1,2,3,4,5]
#ax.set(xticks=xticks, xticklabels=map(lambda v:'%.1f'%v,xticks))
#ax.set_xlim(0.4,5)
ax.set_xlim(0.4,3)
ax.set_ylim(0,2.5)

P.show()