Application Center - Maplesoft

App Preview:

Van der Waals equation of state (I)

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application




I. The isotherms of van der Waals

By Christian Viales M., Universidad de Costa Rica

(made with Maple 12.00)

 

Because of the great educational value of van der Waals' equation of state in teaching thermodynamics of real substances, this Maple worksheet is the first of a series that intends to facilite the visualization and calculation of properties based on it. Here the reader will find four independent procedures that will plot the vdW isotherms around the critical point, and whose input consists in: A) the vdW constants in SI units (a = [m^6*Pa/((mol^2))], b = [m^3/mol]), B) the temperature window around the critical point that will be plotted (in kelvins), and C) the number of isotherms wanted.

 

The first procedure, isothermsW, plots the equation in the traditional way.

The second, isothermsMW, plots the van der Waals equation using Maxwell's area equalization method to determine the vapor pressure lines. It requires a fifth input variable: D) the tolerance (in pascals) the program will assume as acceptable for equalization.

The third, compresW, will plot the compressibility factor based on the virial expansion at low pressure and plotted around Boyle's temperature (instead of the critical temperature).

And finally, compresMW will plot the compressibility factor based on Maxwell's method around the critical temperature (it also requires the tolerance).

 

restart:

NOTE: For hydrogen, helium and other substances with very low critical temperatures do not use Delta*T < T[c]/2.

isothermsW:=proc(a::positive,b::positive,DT::positive,n::posint)
             local TC,PC,P,T,V,i,plo,pl,R:
             description "This procedure plots the n+1 isotherms of van der Waals in a range of DT around the critical temperature and of the gas described by the constants (a,b) in SI units":
             uses plots:
               pl:=printlevel: printlevel:=0: R:=8.31447:
               TC:=8*a/(27*R*b): PC:=a/(27*b^2): T[0]:=TC-DT/2:
               if T[0]<0 then error "Invalid DT: %1" end if;
               if n<1 then error "Invalid n: %1" end if;
               P:=(V,T)->R*T/(V-b)-a/V^2:
               for i from 0 to n do
                 T[i]:=T[0]+i*DT/n:
                 plo[i]:=plot(P(V,T[i]),V=1.1*b..10*b,color=COLOR(RGB,i/n,1-i/n,1),legend=cat(trunc(T[i])," K"),view=[b..10*b, 0..2*PC]):
               end do:
               printlevel:=pl:
               display({seq(plo[i],i=0..n)},labels=["V(m3/mol)","P(Pa)"]);
           end proc:

isothermsMW:=proc(a::positive,b::positive,DT::positive,n::posint,tol::positive)
               local TC,PC,p0,p1,q,A,B,V1,V2,V3,P,m1,m2,f,v,T,V,i,c,plo,pl,R:
               description "This procedure plots the n+1 isotherms of van der Waals, after Maxwell's construction, in a range of DT around the critical temperature and of the gas described by the constants (a,b) in SI units. tol defines the difference in pressure the procedure will tolerate as criteria of area equalization":
               uses plots:
                 pl:=printlevel: printlevel:=0: R:=8.31447:
                 TC:=8*a/(27*R*b): PC:=a/(27*b^2): T[0]:=TC-DT/2:
                 if T[0]<0 then error "Invalid DT: %1" end if;
                 if n<1 then error "Invalid n: %1" end if;
                 P:=(V,T)->R*T/(V-b)-a/V^2:
                 for i from 0 to n do
                   T[i]:=T[0]+i*DT/n: p0:=P(3.375*b,T[i]): p1:=p0: q:=p0/10*tol:
                   A,B:=0,0:
                   if T[i]<TC then
                     for c from 0 to 1000 do
                       V1,V2,V3:=solve(p0=R*T[i]/(V-b)-a/V^2,V):
                       m1:=A-B:
                       A:=int(p0-P(V,T[i]),V=V1..V2):
                       B:=int(P(V,T[i])-p0,V=V2..V3):
                       m2:=A-B:
                       if A>B then p0:=p0-q: end if:
                       if A<B then p0:=p0+q: end if:
                       if abs(p0-p1)<tol then break end if:
                       q:=abs(p0-p1)/2:
                       p1:=p0:
                     end do: f[i]:=V->piecewise(V1<=V and V<=V3,p0,P(V,T[i])):
                     else f[i]:=V->P(V,T[i]):
                   end if:
                   plo[i]:=plot(f[i](v),v=1.1*b..10*b,color=COLOR(RGB,i/n,1-i/n,0),legend=cat(trunc(T[i])," K"),view=[b..10*b, 0..2*PC]):
                 end do:
                 printlevel:=pl:
                 display({seq(plo[i],i=0..n)},labels=["V(m3/mol)","P(Pa)"]);
             end proc:

compresW:=proc(a::positive,b::positive,DT::positive,n::posint)
            local PC,TB,P,T,V,Z,eq,x,i,B,c,d,plo,pl,R:
            description "This procedure plots the n+1 isotherms of the compressibility factor calculated from van der Waals' equation in a range of DT around the Boyle's temperature and of the gas described by the constants (a,b) in SI units":
            uses plots:
              pl:=printlevel: printlevel:=0: R:=8.31447:
              TB:=a/(R*b): PC:=a/(27*b^2): T[0]:=TB-DT/2:
              if T[0]<0 then error "Invalid DT: %1" end if;
              if n<1 then error "Invalid n: %1" end if;
              for i from 0 to n do
                T[i]:=T[0]+i*DT/n:
                eq:=convert(eval(subs(x=1/V,taylor(subs(V=1/x,(R*T[i]/(V-b)-a/V^2)*V/(R*T[i])),x=0,4))),polynom);
                B:=coeff(eq,V,-1): c:=coeff(eq,V,-2): d:=coeff(eq,V,-3):
                Z:=(P,T)->1+B/(R*T)*P+(c-B^2)/(R*T)^2*P^2+(d-3*B*c-2*B^3)/(R*T)^3*P^3:
                plo[i]:=plot(Z(P,T[i]),P=0..PC,color=COLOR(RGB,i/n,0,1-i/n),legend=cat(trunc(T[i])," K")):
              end do:
              printlevel:=pl:
              display({seq(plo[i],i=0..n)},labels=["P(Pa)","Z"]);
          end proc:

compresMW:=proc(a::positive,b::positive,DT::positive,n::posint,tol::positive)
             local PC,TC,P,T,V,Z,eq,x,i,B,c,d,plo,pl,R:
             description "This procedure plots the n+1 isotherms of the compressibility factor calculated from van der Waals' equation in a range of DT around the Boyle's temperature and of the gas described by the constants (a,b) in SI units. In case the datum are within the vdW waves, it invokes the Maxwell procedure to generate a single answer":
             uses plots,RealDomain:
               pl:=printlevel: printlevel:=0: R:=8.31447:
               TC:=8*a/(27*R*b): PC:=a/(27*b^2): T[0]:=TC-DT/2:
               if T[0]<0 then error "Invalid DT: %1" end if;
               if n<1 then error "Invalid n: %1" end if;
               V:=proc(P,T)
                    local f,y,v1,v2,v3,p0,p1,p2,q,A,B,m1,m2,c:
                    y:=solve(P=R*T/(v-b)-a/v^2,v):
                    if nops({y})=3 then
                      A,B,p0,p1,q:=0,0,P,P,P/10*tol:
                      for c from 0 to 1000 do
                        v1,v2,v3:=solve(p0=R*T/(v-b)-a/v^2,v):
                        m1:=A-B:
                        A:=int(p0-(R*T/(v-b)-a/v^2),v=v1..v2):
                        B:=int((R*T/(v-b)-a/v^2)-p0,v=v2..v3):
                        m2:=A-B:
                        if A>B then p0:=p0-q: end if:
                        if A<B then p0:=p0+q: end if:
                        if abs(p0-p1)<tol then break end if:
                        q:=abs(p0-p1)/2:
                        p1:=p0:
                      end do:
                      f:=p2->piecewise(p2<p0,max(v1,v2,v3),p0<p2,min(v1,v2,v3),v1+v2+v3-min(v1,v2,v3)-max(v1,v2,v3)):
                      else f:=p2->y:
                    end if:
                    f(P);
               end proc:
               Z:=(P,T)->V(P,T)/(V(P,T)-b)-a/(R*T*V(P,T)):
               for i from 0 to n do
                 T[i]:=T[0]+i*DT/n:
                 plo[i]:=plot(Z(P,T[i]),P=0..5*PC,color=COLOR(RGB,0.5*(1-i/n),i/n,1-i/n),legend=cat(trunc(T[i])," K")):
               end do:
               printlevel:=pl:
               display({seq(plo[i],i=0..n)},labels=["P(Pa)","Z"],view=[0..5*PC,0..1]);
           end proc:

 

And as an example, the isotherms for carbon dioxide.

a,b:=0.3640,4.27e-5:

isothermsW(a,b,140,10);

 

isothermsMW(a,b,140,10,0.5e-3);

 

compresW(a,b,140,10);

 

compresMW(a,b,140,10,0.5e-3);

 

NOTE: Either the procedure written or Maple fails to recognize the discontinous nature of the liquid-gas transition, so the plot is not as "pretty" as one would like.