La ecuación de van der Waals

La ecuación de estado de un gas ideal es PV=nRT

P es la presión del gas, V es el volumen, T es la temperatura, n es el número de moles y R la constante de los gases
Describe aproximadamente la conducta de los gases reales a muy bajas presiones.

La ecuación de van der Waals tiene en cuenta el volumen finito de las moléculas y las fuerzas atractivas que una molécula ejerce sobre otra a distancias muy cercanas entre ellas.

( P+ n 2 a V 2 )( Vnb )=nRT

Las constantes a y b son característicos de cada gas y se obtienen a partir de los datos de la presión, Pc, volumen Vc y la temperatura Tc crítica. El punto crítico es un punto de inflexión de la isoterma Tc en el diagrama P-V de modo que se cumple que

( P V ) T =0 nRT ( Vnb ) 2 + 2 n 2 a V 3 =0 ( 2 P V 2 ) T =0 2nRT ( Vnb ) 3 6 n 2 a V 4 =0

De estas dos ecuaciones obtenemos, el volumen Vc y la temperatura Tc crítica.

V c =3nb T c = 8a 27Rb

Sustituyendo Vc y Tc en la ecuación de van der Waals obtenemos la presión crítica, Pc,

P c = a 27 b 2

Escribimos la ecuación de van der Waals universal más simple en términos de nuevas variables: p=P/Pc, v=V/Vc y t=T/Tc

( p+ 3 v 2 )( 3v1 )=8t

En el punto crítico pc=1, tc=1 y vc=1.

Como vamos a dibujar las isotermas en un diagrama P-V, despejamos la presión p.

p= 8t 3v1 3 v 2

Creamos un script para representar las isotermas de temperaturas t=0.8,0.9,1.0,1.1 y 1.2. Señalamos con una marca el punto crítico tal como se aprecia en la figura

f=@(t,x) 8*t./(3*x-1)-3./(x.^2);
hold on
axis([0.5 4 0 2])
for t=0.8:0.1:1.2;
    v=linspace(0.5,4,100);
    p=f(t,v);
    plot(v,p)
end
plot(1,1,'ro')
xlabel('V_R')
ylabel('P_R')
grid on
hold off

Examinamos las isotermas por debajo de la temperatura crítica t<1.

En la figura vemos que para una determinada presión p, la recta p=cte corta a la curva en tres puntos A, B y C. Los volúmenes v1, v2 y v3 que corresponden a estos tres puntos se obtienen resolviendo la ecuación cúbica

v 3 8t+p 3p v 2 + 3 p v 1 p =0

Para obtener las raíces, utilizamos la función roots de MATLAB o la función raices_3 que nos calcula las tres raíces exactas de una ecuación de tercer grado.

function x = raices_3(p)
    Q=(p(2)*p(2)-3*p(3))/9;
    R=(2*p(2)^3-9*p(2)*p(3)+27*p(4))/54;
    x=zeros(3,1);
    if (R*R)<(Q^3)
        tetha=acos(R/sqrt(Q^3));
        x(1)=-2*sqrt(Q)*cos(tetha/3)-p(2)/3;
        x(2)=-2*sqrt(Q)*cos((tetha+2*pi)/3)-p(2)/3;
        x(3)=-2*sqrt(Q)*cos((tetha-2*pi)/3)-p(2)/3;
    else
        A=-sign(R)*nthroot(abs(R)+sqrt(R*R-Q^3),3);
        if A==0
            B=0;
        else
            B=Q/A;
        end
        x(1)=(A+B)-p(2)/3;
        x(2)=-(A+B)/2-p(2)/3+(sqrt(3)*(A-B)/2)*sqrt(-1); 
        x(3)=-(A+B)/2-p(2)/3-(sqrt(3)*(A-B)/2)*sqrt(-1);
    end
end

Ordenamos los valores de las raíces de menor a mayor mediante la función sort de MATLAB.

La isoterma de temperatura t, tiene un mínimo y un máximo local. Determinamos sus coordenadas (vm1, pm1) y (vm2, pm2) calculando la derivada primera e igualando a cero.

p= 8t 3v1 3 v 2 ( p v ) t =0 24t ( 3v1 ) 2 6 v 3 =0 v 3 9 4t v 2 + 3 2t v 1 4t =0

Tenemos de nuevo, que calcular las raíces de una ecuación cúbica.

En la zona de transición entre la fase gaseosa y la fase líquida la presión no oscila, se mantiene constante y el potencial químico se mantiene constante. La regla de Maxwell elimina el comportamiento oscilante de la isoterma de la ecuación de van der Waals y la sustituye por un segmento horizontal de presión pr tal que el área comprendida entre el segmento horizontal AB y la isoterma es igual al área comprendida entre la isoterma y el segmento horizontal BC, tal como se ve en la figura más abajo.

En la figura, vemos la isoterma de temperatura t=0.9,

El el área A1comprendida entre el segmento horizontal AB y la isoterma se calcula integrando

p ( v 2 v 1 ) v 1 v 2 p · d v = p ( v 2 v 1 ) v 1 v 2 ( 8 t 3 v 1 3 v 2 ) · d v = p ( v 2 v 1 ) 8 3 t ln 3 v 2 1 3 v 1 1 3 ( 1 v 2 1 v 1 )

El el área A2 comprendida entre el segmento horizontal BC y la isoterma se calcula integrando

v 2 v 3 p · d v p ( v 3 v 2 ) = 8 3 t ln 3 v 3 1 3 v 2 1 + 3 ( 1 v 3 1 v 2 ) p ( v 3 v 2 )

Creamos una función denominada igualArea que calcule la diferencia entre las dos áreas A1-A2.

function res=igualArea(p,t)
    pol=[1 -(p+8*t)/(3*p) 3/p -1/p];
    v = sort(raices_3(pol));
    area_1 = p*(v(2)-v(1))-8/3*t*log((3*v(2)-1)
/(3*v(1)-1))-3*(v(1)-v(2))/(v(1)*v(2));
    area_2 =8/3*t*log((3*v(3)-1)/(3*v(2)-1))+3*(v(2)-v(3))
/(v(2)*v(3))-p*(v(3)-v(2));
    res=area_1-area_2;
end

Creamos un script para realizar las siguientes tareas:.

  1. Establezca el valor de la temperatura t reducida, menor que la unidad
  2. Represente la isoterma de color rojo en el diagrama P-V.
  3. Calcule el mínimo pm1 y máximo pm2 local resolviendo la ecuación de tercer grado en v.
  4. La presión pr buscada estará en el intervalo ( pm1, pm2 ), o bien en el intervalo ( 0, pm2 ) si pm1 fuese negativa (la presión es siempre positiva). El manejador (handle) de esta función igualArea se lo pasaremos a la función fzero de MATLAB para que calcule el valor pr de la presión que hace la que diferencia de áreas sea nula.
  5. Una vez obtenido el valor de la presión pr, se calcula las raíces de la ecuación cúbica para determinar los volúmenes v1, v2 y v3. El primero corresponde al volumen de la fase líquida y el tercero al del gas, para la presión pr y la temperatura t
  6. Rellene de color amarillo las dos regiones, para mostrar visualmente que las áreas son iguales.
t=input('temperatura reducida (0.8-1): ');
%gráfica de la isoterma
f=@(x) 8*t./(3*x-1)-3./(x.^2);
Vr=linspace(0.5,4,100);
Pr=f(Vr);
hold on
axis([0 4 0 1.4])
plot(Vr,Pr,'r')
xlabel('V_R')
ylabel('P_R')

%máximo y mínimo
pol=[1 -9/(4*t) 3/(2*t) -1/(4*t)];
z=sort(raices_3(pol));
plot(z(1),f(z()),z(2),f(z(2)),'k*',z(3),f(z(3)),'k*')
if f(z(2))>0
    a=f(z(2))+0.01;
else
    a=0.01;
end
b=f(z(3))-0.01;
%calcula p para que las áreas sean iguales
f1=@(p)igualArea(p,t);
p=fzero(f1,[a b]);
%igualArea(p,t) %para comprobar que las areas son iguales
plot([0 4],[p p],'k--') 
pol=[1 -(p+8*t)/(3*p) 3/p -1/p];
v = sort(raices_3(pol));
plot(v(1),p,'bo',v(2),p,'bo',v(3),p,'bo')

%rellenar la areas
xx = [v(1) Vr(Vr >= v(1) & Vr <= v(2))  v(2)];
yy = [f(v(1)) Pr(Vr >= v(1) & Vr <= v(2)) f(v(2))];
fill(xx,yy,'y')

xx = [v(2) Vr(Vr >= v(2) & Vr <= v(3)) v(3)];
yy = [f(v(2)) Pr(Vr >= v(2) & Vr <= v(3)) f(v(3))];
fill(xx,yy,'y')
plot(Vr,Pr,'r')
%imprimir datos
text(1,1.35,sprintf('presión %1.3f',p))
text(1,1.25,sprintf('volumen gas %1.3f',v(3)))
text(1,1.15,sprintf('volumen líquido %1.3f',v(1)))

hold off

Ahora vamos a dibujar varias isotermas por encima y por debajo del punto crítico t=1. Sustituiremos la región ondulada entre vliquido y vgas a la presión pr por un segmento y uniremos los puntos que marcan la región de coexistencia de las dos fases, tal como se muestra en la figura.

Creamos un script para que realice las siguientes tareas.

  1. Dibuje las isotermas t≥1 (para t=1,1.05,1.1,1.15 y 1.2)
  2. Dibuje las isotermas para t<1 (para t=0.8, 0.85, 0.9, 0.95).
  3. Se añade a la matriz bidimensional las coordendas del punto crítico v=1 y p=1. Ordenamos la matriz de datos (v,p) utilizando la función MATLAB sortrows, y unimos los puntos (v,p) mediante líneas de color azul tal como se muestra en la figura.
Vr=linspace(0.1,4,100);
hold on
i=1;
PV=zeros(9,2);
axis([0.4 4 0.2 1.4])
for t=0.8:0.05:1.2;
    f=@(x) 8*t./(3*x-1)-3./(x.^2);
    if t<1
        %calcula máximo y mínimo local de p
        pol=[1 -9/(4*t) 3/(2*t) -1/(4*t)];
        z=sort(raices_3(pol));
        if f(z(2))>0
            a=f(z(2))+0.01;
        else
            a=0.01;
        end
        b=f(z(3))-0.01;
        %calcula p para que las áreas sean iguales
        f1=@(p)igualArea(p,t);
        p=fzero(f1,[a b]);
        pol=[1 -(p+8*t)/(3*p) 3/p -1/p];
        v = sort(raices_3(pol));
        %guarda Vlíquido, Vgas y la presión constante p 
        PV(i,1)=v(1);PV(i,2)=p;
       PV(i+1,1)=v(3);PV(i+1,2)=p;
       i=i+2;
 %gráfica
        Pr=f(Vr);
        Vf = [Vr(Vr <= v(1)) v(1) v(3) Vr(Vr>= v(3))];
        Pf = [Pr(Vr <= v(1)) f(v(1)) f(v(3)) Pr(Vr>= v(3))];
        plot(Vf,Pf,'r')
    else %en el punto crítico y por encima
        Pr=f(Vr);
        plot(Vr,Pr,'r')
    end 
end
%une los puntos (p Vlíquido), (p Vgas) cambio de fase
PV(i,1)=1; PV(i,2)=1; %punto crítico
PVs=sortrows(PV,1);
plot(PVs(:,1),PVs(:,2),'b')
plot(1,1,'bo') %punto crítico
xlabel('V_R')
ylabel('P_R')
hold off