CARDANO VIETTA
Este programa fue diseñado para aplicar el método Cardano-Vietta con el fin
de hallar las raíces de una función cubica de volumen, la cual fue
determinada a partir de la ecuación de gas real de Redlich-Wkong.
Aquí se recorre las temperaturas manteniendo una presión fija, una vez
termina un recorrido de temperaturas, cambia el valor de la presión en
aumento; las unidades fueron tomadas consistentes para R con unidades de
bar*m3/mol*K.
Este programa es la versión final obtenida, a diferencia del anterior
intento, este no contiene restricciones para la creación de una matriz (en la
que se tabulan las tres raíces a correspondientes presiones y temperaturas),
sin embargo, estas restricciones no fueron necesarias puesto que todas las
raíces obtenidas fueron positivas. El cambio efectuado para esta última
versión fue definir a “a” y a “b” como constantes numéricas precalculadas y
definir a “R’ con las unidades anteriormente mencionadas, en este sentido, no
fue necesario mantener la definición de P crítica y T critica dentro del
programa.
Al final se presenta una matriz de 5 columnas mostrando, respectivamente, las
tres raíces halladas (X1, X2 y X3), la presión correspondiente a las raíces,
y las temperaturas correspondientes; todos los valores se presentan en
notación científica de (10-2).
Para su graficación, se tomaron en el eje “x” los valores de las dos raíces
más convenientes y en el eje “y” los valores de las presiones hasta 204.5
bar.
clc;
i=0;
for P=20:4.1:225
i=(i+1);
for T=298:4.02:700
R=0.00008314472;
a= 0.000005526914954;
b=0.00002081832548;
A1=(-((R*T)/P));
A2=(-(((P*((b)^2))+(R*T*b)-a)/P));
A3=(-((a*b)/P));
Q=(((A1^2)-(3*A2))/9);
M=(((2*(A1^3))-(9*A1*A2)+(27*A3))/54);
c= ((Q^3)-(M^2));
k=(atan(sqrt(((Q^3)/(M^2))-1)));
if c>=0
if M>0
X1=((-2)*sqrt(Q)*(cos(k/3))-(A1/3));
X2=((-2)*sqrt(Q)*(cos((k+(2*pi))/3))-(A1/3));
X3=((-2)*sqrt(Q)*(cos((k+(4*pi))/3))-(A1/3));
dT=(T*1);
end
if M<0
X1=(2*sqrt(Q)*(cos(k/3))-(A1/3));
X2=(2*sqrt(Q)*(cos((k+(2*pi))/3))-(A1/3));
X3=(2*sqrt(Q)*(cos((k+(4*pi))/3))-(A1/3));
dT=(T*1);
end
else
if M>0
q=1;
end
if M<0
q=(-1);
end
X= (((-q)*(((sqrt((M^2)-(Q^3))+(abs(M)))^(1/3))+(Q/(((sqrt((M^2)-
(Q^3)))+(abs(M)))^(1/3)))))-(A1/3));
end
end
m(i,1)=X1;
m(i,2)=X2;
m(i,3)=X3;
m(i,4)=P;
m(i,5)=dT;
end
format long
m
2.5
1.5
0.5
0
0 0.000005 0.00001 0.000015 0.00002 0.000025
2.5
1.5
0.5
0
0 0.000005 0.00001 0.000015 0.00002 0.000025
Aquí se recorre las temperaturas manteniendo una presión fija, una vez
termina un recorrido de temperaturas, cambia el valor de la presión en
aumento; las unidades fueron tomadas del SI.
Este programa solo acepta producir la matriz final (en la que se tabulan los
pares de raíces y sus presiones) si se obtienen por lo menos dos raíces
positivas, en tal caso, se definiría el valor de la isoterma para las raíces
que cumplen la condición en esta presión. Sin embargo, no se arroja una
matriz debido a que no se encuentran dos raíces positivas.
clc;
Tc=658.416405;
Pc=22500;
i=0;
for P=2000:410:22500
i=(i+1);
for T=298:4.02:700
R=8.3144598;
a=((0.42748*(((R^2)*(Tc^2.5))/(Pc*(T^2.5)))));
b=0.02108;
A1=(-((R*T)/P));
A2=(-(((P*((b)^2))+(R*T*b)-a)/P));
A3=(-((a*b)/P));
Q=(((A1^2)-(3*A2))/9);
M=(((2*(A1^3))-(9*A1*A2)+(27*A3))/54);
c= ((Q^3)-(M^2));
k=(atan(sqrt(((Q^3)/(M^2))-1)));
if c>=0
if M>0
X1=((-2)*sqrt(Q)*(cos(k/3))-(A1/3));
X2=((-2)*sqrt(Q)*(cos((k+(2*pi))/3))-(A1/3));
X3=((-2)*sqrt(Q)*(cos((k+(4*pi))/3))-(A1/3));
if X2>0 && X3>0
m(i,1)=X1;
m(i,2)=X2;
m(i,3)=X3;
m(i,4)=P;
dT=(T*1);
m(i,5)=dT;
end
if X1>=0 && X3>=0
m(i,1)=X1;
m(i,2)=X2;
m(i,3)=X3;
m(i,4)=P;
dT=(T*1);
m(i,5)=dT;
end
if X1>=0 && X2>=0
m(i,1)=X1;
m(i,2)=X2;
m(i,3)=X3;
m(i,4)=P;
dT=(T*1);
m(i,5)=dT;
end
end
if M<0
X1=(2*sqrt(Q)*(cos(k/3))-(A1/3));
X2=(2*sqrt(Q)*(cos((k+(2*pi))/3))-(A1/3));
X3=(2*sqrt(Q)*(cos((k+(4*pi))/3))-(A1/3));
if X2>0 && X3>0
m(i,1)=X1;
m(i,2)=X2;
m(i,3)=X3;
m(i,4)=P;
dT=(T*1);
m(i,5)=dT;
end
if X1>=0 && X3>=0
m(i,1)=X1;
m(i,2)=X2;
m(i,3)=X3;
m(i,4)=P;
dT=(T*1);
m(i,5)=dT;
end
if X1>=0 && X2>=0
m(i,1)=X1;
m(i,2)=X2;
m(i,3)=X3;
m(i,4)=P;
dT=(T*1);
m(i,5)=dT;
end
end
end
if c<0
if M>0
q=1;
end
if M<0
q=(-1);
end
X= (((-q)*(((sqrt((M^2)-(Q^3))+(abs(M)))^(1/3))+(Q/(((sqrt((M^2)-
(Q^3)))+(abs(M)))^(1/3)))))-(A1/3));
end
end
end
format long
m
Este programa fue el último intento funcional para la aplicación del método Cardano-Vietta; ya se han
hecho aquí las correcciones correspondientes al uso de la función atan. Este programa hace un recorrido
de presiones con recorrido interno de temperaturas, pero toma siempre la última temperatura del
recorrido(Tc) para hallar las raíces de la función cuadrática, dos de las cuales son negativas. Esto ocurre
al no tener criterios que discriminen raíces tomadas para la matriz.
clc;
Tc=658.416405;
Pc=22500;
i=0;
for P=2000:410:22500
i=(i+1);
for T=298:4.02:700
R=8.3144598;
a=((0.42748*(((R^2)*(Tc^2.5))/(Pc*(T^2.5)))));
b=0.02108;
A1=(-((R*T)/P));
A2=(-(((P*((b)^2))+(R*T*b)-a)/P));
A3=(-((a*b)/P));
Q=(((A1^2)-(3*A2))/9);
M=(((2*(A1^3))-(9*A1*A2)+(27*A3))/54);
c= ((Q^3)-(M^2));
k=(atan(sqrt(((Q^3)/(M^2))-1)));
if c>=0
if M>0
X1=((-2)*sqrt(Q)*(cos(k/3))-(A1/3));
X2=((-2)*sqrt(Q)*(cos((k+(2*pi))/3))-(A1/3));
X3=((-2)*sqrt(Q)*(cos((k+(4*pi))/3))-(A1/3));
end
if M<0
X1=(2*sqrt(Q)*(cos(k/3))-(A1/3));
X2=(2*sqrt(Q)*(cos((k+(2*pi))/3))-(A1/3));
X3=(2*sqrt(Q)*(cos((k+(4*pi))/3))-(A1/3));
end
end
if c<0
if M>0
q=1;
end
if M<0
q=(-1);
end
X= (((-q)*(((sqrt((M^2)-(Q^3))+(abs(M)))^(1/3))+(Q/(((sqrt((M^2)-
(Q^3)))+(abs(M)))^(1/3)))))-(A1/3));
end
end
m(i,1)=X1;
m(i,2)=X2;
m(i,3)=X3;
m(i,4)=P;
end
format long
m
A continuacion, la ultima columna de la matriz reportada en command window,
registra que la temperatura guarada para los calculos de las raices es la Tc,
ya que no se tienen restricciones para la aceptacion de raices.
----------------------------------------------------------------------------------------------------------------------------------
Inicialmente y por error, se usó la función acot para el cálculo de del ángulo teta; a continuación, se
muestran las gráficas que se obtuvieron y la posterior corrección del cálculo de teta en el programa.
El siguiente es el primer error en el proceso del diseño del programa, del cual se tiene registro. En un
principio, se definió “a” como una constante del agua, pero, al tener en cuenta el recorrido de las
temperaturas se determinó la dependencia de “a” hacia T y se reemplazó su valor fijo por la fórmula que
lo define.
