Application Center - Maplesoft
  Submit your work
Application Center Applications Analysis of basic equations of state (II) Preview

App Preview:

Analysis of basic equations of state (II)

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

Learn about Maple
Download Application



II. Analysis of basic equations of state

By Christian Viales M., Universidad de Costa Rica

(Made with Maple 12.00)

 

As an advanced introduction to real gases, one can calculate and later compare the accuracy of the thermodynamic data predicted by the ideal gas equation and other equations of state. In order to achieve that, one must first calculate the constants of that substance from experimental data (generally the critical point) and then calculate the volumetric coefficients. This Maple worksheet, the second of the series, pretends to facilite the calculation of the aforementioned properties for the most common two-constant fluid EOS.
NOTE: after executing the worksheet with one EOS, restart it to work with another.

 

> 

restart:
with(RealDomain):
_EnvExplicit:=true:
eqid := P = R*T/V:                                        # Ideal Gas
eqvdW:= P = R*T/(V-b)-a/V^2:                              # van der Waals
eqD  := P = R*T/(V-b)*exp(-a/(V*R*T)):                    # Dieterici
eqB  := P = R*T/(V-b)-a/(V^2*T):                          # Berthelot
eqRK := P = R*T/(V-b)-a/(sqrt(T)*V*(V+b)):                # Redlich-Kwong
eqPR := P = R*T/(V-b)-(a)/(V^2+2*b*V-b^2):                # Peng-Robinson

 

First chose the equation you want to work with.

 

> 

eq:=eqvdW;
const:=a,b:

eq := P = R*T/(V-b)-a/V^2

 (1)
> 

 

 

Now the constants can be derived from the critical point (and viceversa)...

 

> 

eq1:= subs({T=Tc,V=Vc,P=Pc},solve({0=diff(rhs(eq),V),0=diff(rhs(eq),V$2)},{const}));
eq2:= solve({eq1[1],eq1[2],subs({T=Tc,V=Vc,P=Pc},eq)},{Pc,Tc,Vc});
assign(eq2):
eq3:= Zc=simplify(Pc*Vc/(R*Tc));
eq3:= Zc=evalf(rhs(%));

eq1 := {a = 9*Vc*R*Tc/8, b = Vc/3}

eq2 := {Pc = a/(27*b^2), Tc = 8*a/(27*b*R), Vc = 3*b}

eq3 := Zc = 3/8

eq3 := Zc = .3750000000

 (2)
> 

 

 

... and the reduced form can also be computed.

NOTE: For most of the EOS the results obtained apparently differ to those calculated "by hand", nonetheless must be equivalent. For example, for van der Waals:

P[R] = (8*T[R]*V[R]^2-9*V[R]+3)/((3*V[R]-1)*(V[R]^2)) ... is equivalent to... (Pi+3/((Phi^2)))*(Phi-1/3) = 8*Tau/3

 

> 

eq4:= simplify(isolate(subs({P=P[R]*Pc,T=T[R]*Tc,V=V[R]*Vc},eq),P[R]));
unassign('Pc','Tc','Vc'):

eq4 := P[R] = (8*T[R]*V[R]^2-9*V[R]+3)/((3*V[R]-1)*V[R]^2)

 (3)
> 

 

 

Just as important, the isothermal compressibility, the isobaric expansion coefficient and the difference between the isobaric and isochoric heat capacities can be also computed (here they are computed in two equivalent expressions).

 

> 

eq5a:= kappa = (-implicitdiff(eq,V(P,T),P)/V);
eq5:= kappa = simplify(subs({isolate(eq,T)},-implicitdiff(eq,V(P,T),P)/V));
eq6a:= alpha = (implicitdiff(eq,V(P,T),T)/V);
eq6:= alpha = simplify(subs({isolate(eq,T)},implicitdiff(eq,V(P,T),T)/V));
eq7a:= C[P]-C[V]=simplify(subs({eq5a,eq6a},T*V*alpha^2/kappa));
eq7:= C[P]-C[V]=simplify(subs({eq5,eq6},T*V*alpha^2/kappa));

eq5a := kappa = V^2*(V^2-2*b*V+b^2)/(R*T*V^3-2*a*V^2+4*a*V*b-2*a*b^2)

eq5 := kappa = (V-b)*V^2/(P*V^3-a*V+2*a*b)

eq6a := alpha = R*(V-b)*V^2/(R*T*V^3-2*a*V^2+4*a*V*b-2*a*b^2)

eq6 := alpha = V^2*R/(P*V^3-a*V+2*a*b)

eq7a := C[P]-C[V] = R^2*T*V^3/(R*T*V^3-2*a*V^2+4*a*V*b-2*a*b^2)

eq7 := C[P]-C[V] = T*V^3*R^2/((P*V^3-a*V+2*a*b)*(V-b))

 (4)

 

Maple calculates the virial expansion for the compressibility factor and Boyle's temperature...

 

> 

eq8:= Z=subs(x=1/V,taylor(subs(V=1/x,rhs(eq)*V/(R*T)),x=0,5));
eq9:= Z=1+B/V+C/V^2+D/V^3;
Virial1:={B=simplify(coeff(convert(rhs(eq8),polynom),V,-1)), C=simplify(coeff(convert(rhs(eq8),polynom),V,-2)), D=simplify(coeff(convert(rhs(eq8),polynom),V,-3))};
Boyle:= T[B]=solve(0=coeff(convert(rhs(eq8),polynom),V,-1),T);

eq8 := Z = 1+(R*T*b-a)/(R*T*V)+b^2/V^2+b^3/V^3+O(1/((V^4)))

eq9 := Z = 1+B/V+C/V^2+D/V^3

Virial1 := {B = (R*T*b-a)/(R*T), C = b^2, D = b^3}

Boyle := T[B] = a/(b*R)

 (5)
> 

 

 

... and the equivalence between the pressure and volume virial forms for the compressibility factor.

 

> 

eq10:= Z=1+Beta*P+Gamma*P^2+Delta*P^3;
eq11:= Z=collect(expand(subs({eq9},subs({P=Z*R*T/V},rhs(eq10)))),1/V):
Vir2:={B=coeff(convert(rhs(eq11),polynom),V,-1), C=coeff(convert(rhs(eq11),polynom),V,-2), D=coeff(convert(rhs(eq11),polynom),V,-3)}:
Vir[1]:= Beta=solve(Vir2[1],Beta):
Vir[2]:= Gamma=solve(subs({%},Vir2[2]),Gamma):
Vir[3]:= Delta=solve(subs({%,%%},Vir2[3]),Delta):
Virial2:= {Vir[1],Vir[2],Vir[3]};
Virial3:= simplify(subs(Virial1,Virial2));

eq10 := Z = 1+Beta*P+Gamma*P^2+Delta*P^3

Virial2 := {Beta = B/(R*T), Delta = (D-3*B*C+2*B^3)/(R^3*T^3), Gamma = -(-C+B^2)/(R^2*T^2)}

Virial3 := {Beta = (R*T*b-a)/(R^2*T^2), Delta = -a*(3*b^2*R^2*T^2-6*R*T*b*a+2*a^2)/(R^6*T^6), Gamma = a*(2*R*T*b-a)/(R^4*T^4)}

 (6)
> 

 

 

The internal pressure can be calculated just as seen later in the course.

 

> 

eq12:= Diff(U,V)[T]=simplify(subs({eq,eq5a,eq6a},T*alpha/kappa-P));

eq12 := (Diff(U, V))[T] = a/V^2

 (7)
>