UNIVERSIDAD AUTONOMA JUAN MISEL SARACHO
MATLAB
CARACTERÍSTICAS BÁSICAS DE MATLAB
1.-ECUACIONES TERMODINAMICAS TERMODINAMICAS PARA GASES REALES 1.1.-ECUACION DE VAN DER WALLS 1.2.-ECUACION DE REDLICH-KWONG
1.3.-ECUACION DE SOAVE- REDLICH-KWONG 1.4.-ECUACION DE PENG-ROBINSON 1.5.-ECUACION DE BEATLIE BRIDGEMAN 1.6.-ECUACION VIRIAL
1.7.-ECUACION DE BENEDICT-WEBB-RUBIN BENEDICT-WEBB-RUBIN 2.-LISTADO DE PROGRAMAS
3.-EJERCICIOS RESUELTOS 4.-PROBLEMAS PROPUESTOS :
Tablas para los gases reales
Introducción: Este es el aspecto que presenta la versión R2009-b de MatLab, que será la que utilizaremos utili zaremos este este curso:
MATLAB es el nombre abreviado de “MATriz LABoratory”. Es un progra ma para realizar cálculos numéricos con vectores y matrices, y por tanto se puede trabajar también con números escalares (tanto reales como complejos), complejos), con cadenas de caracteres y con otras estructuras de información más complejas. Matlab es un lenguaje de alto rendimiento para cálculos técnicos, es al mismo tiempo un entorno y un lenguaje de programación. Uno de sus puntos fuertes es que permite construir nuestras propias herramientas reutilizables. Podemos crear fácilmente nuestras propias funciones y programas especiales (conocidos como M-archivos) en código Matlab, los podemos agrupar en Toolbox (también llamadas librerías): colección especializada de Marchivos para trabajar en clases particulares de problemas. Matlab, a parte del cálculo matricial y álgebra lineal, también puede manejar polinomios, funciones, ecuaciones diferenciales ordinarias, gráficos … Este es el aspecto que presenta la versión R2009-b de MatLab, que será la que utilizaremos utili zaremos este este curso:
CARACTERÍSTICAS BÁSICAS EL ESPACIO DE TRABAJO DE MATLAB Nada más abrir Matlab (podemos hacerlo pinchando en el icono que aparece en el escritorio o en su defecto en Inicio->Todos los programas) aparecerá una pantalla como la siguiente:
Todas las sentencias que vamos a utilizar las escribiremos en la ventana Command Window (ventana de comandos). Es la ventana de mayor tamaño. Si queremos información acerca de las variables que estamos utilizando en Matlab podemos verlas en la ventana Workspace (espacio de trabajo) o usar:
Who para obtener la lista de las variables (no de sus valores) para obtener la lista de las variables e información del tamaño, tipo y atributos Whos (tampoco da valores) Para ver esta ventana tenemos que pinchar en la pestaña que tienen este nombre. Está en la parte superior izquierda:
Si lo que queremos es conocer el valor que tiene una variable lo hacemos escribiendo el nombre de la variable y pulsando Intro. Para recordar órdenes previas usamos las flechas del teclado ↑ y ↓. También podemos verlas en la ventana Command History, ventana situada en la parte inferior izquierda:
MATEMÁTICA SENCILLA Matlab ofrece la posibilidad de realizar las siguientes operaciones básicas: Operación suma resta multiplicación división potencia
Símbolo + * / ^
Expresión en Matlab a+b a-b a*b a/b a^b
El orden de precedencia es: Orden de precedencia de operaciones 1º ^ 2º * / 3º + Matlab no tiene en cuenta los espacios. Si queremos que Matlab evalúe la línea pero que no escriba la respuesta, basta escribir punto y coma ( ;) al final de la sentencia. Si la sentencia es demasiado larga para que quepa en una sola línea podemos poner tres puntos ( …) seguido de la tecla Intro para indicar que continúa en la línea siguiente.
FORMATOS DE VISUALIZACIÓN DE NÚMEROS Matlab no cambia la representación interna de un número cuando se escogen distintos formatos, sólo se modifica la forma de visualizarlo.
Tipo format short
format long format short e format long e format short g format long g format short eng format long eng format bank format hex format rat format +
Resultado Ejemplo: >> pi Formato coma fija con 4 dígitos después de la 3.1416 coma (es el formato que viene por defecto) Formato coma fija con 14 o 15 dígitos 3.14159265358979 después de la coma Formato coma flotante con 4 dígitos después 3.1416e+000 de la coma Formato coma flotante con 14 o 15 dígitos 3.141592653589793e+000 después de la coma La mejor entre coma fija o flotante con 4 3.1416 dígitos después de la coma La mejor entre coma fija o flotante con 14 o 3.14159265358979 15 dígitos después de la coma Notación científica con 4 dígitos después de 3.1416e+000 la coma y un exponente de 3 Notación científica con 16 dígitos 3.14159265358979e+000 significantes y un exponente de 3 Formato coma fija con 2 dígitos después de la 3.14 coma Hexadecimal 400921fb54442d18 Aproximación racional 355/113 Positivo, negativo o espacio en blanco +
ACERCA DE LAS VARIABLES Matlab almacena el último resultado obtenido en la variable ans. Las variables son sensibles a las mayúsculas, deben comenzar siempre con una letra, no pueden contener espacios en blanco y pueden nombrarse hasta con 63 caracteres (en versiones anteriores no permitía tantos caracteres). Si se nombra una variable con más de 63 caracteres truncará el nombre de dicha variable. Algunas variables especiales de Matlab: Variable ans pi eps
Definición Variable usada por defecto para almacenar el último resultado Razón de una circunferencia a su diámetro Número más pequeño, tal que cuando se le suma 1, crea un número en coma flotante en el computador mayor que 1 Infinito inf Magnitud no numérica nan iyj i = j = −1 realmin El número real positivo más pequeño que es utilizable realmax El número real positivo más grande que es utilizable
Valor ??? 3.1416 2.2204e-016 Inf NaN 0 + 1.0000i 2.2251e-308 1.7977e+308
Tecleando clear podemos borrar todas las variables del espacio de trabajo, pero no borra lo de las demás ventanas, es decir, no desaparece lo que hay escrito en la ventana de comandos. Tecleando clc borramos lo que hay en la ventana de comandos pero no borra las variables de la memoria del espacio de trabajo. Algunos comandos de Matlab nos facilitan información sobre la fecha, como clock , date o calendar. OTRAS CARACTERÍSTICAS BÁSICAS Los comentarios se escriben después del símbolo de tanto por ciento ( %), de este modo todo lo que se escriba a continuación en la misma línea no será leído por Matlab. Podemos colocar varias órdenes en una línea si se separan correctamente, puede ser:
por comas (,) que hacen que se visualicen los resultados
puntos y comas (;) que suprimen la impresión en pantalla Para cerrar Matlab podemos hacerlo tecleando quit, cerrando con el aspa típico de Windows, entrando en File->Exit Matlab o con las teclas Ctrl+Q.
TRIGONOMETRÍA Función
¿Qué hace? … (x) función trigonométrica con el ángulo expresado en radianes sin (x) cos seno (radianes) (x) tan (x) coseno tangente cosecante secante cotangente csc (x) sec (x) cot (x) …d (x)
sind (x) … …h (x)
sinh (x) … a… (x)
asin (x) … a…d (x)
asind (x) … a…h (x)
función trigonométrica con el ángulo expresado en grados seno (grados) …
función trigonométrica hiperbólica con el ángulo expresado en radianes seno hiperbólico (radianes) …
inversa de la función trigonométrica con el resultado expresado en radianes arco seno (radianes) …
inversa de la función trigonométrica con el resultado expresado en grados arco seno (grados) …
inversa de la función trigonométrica hiperbólica con el resultado expresado en radianes
asinh (x)
arco seno hiperbólico (radianes)
…
…
OPERACIONES BÁSICAS CON MATRICES Símbolo + * .* / ./ \ .\ ^ .^ ' .'
Expresión Operación A + B A - Suma de matrices Resta de matrices Multiplicación de B A * B A matrices .* B A / B Multiplicación elemento a elemento de matrices A ./ B A \ B División de matrices por la derecha A .\ B A ^ n División elemento a elemento de matrices por la derecha A .^ B A ' División de matrices por la izquierda A .' División elemento a elemento de matrices por la izquierda Potenciación (n debe ser un número, no una matriz) Potenciación elemento a elemento de matrices Trasposición compleja conjugada Trasposición de matrices
TEXTO Una cadena de caracteres es texto rodeado por comillas simples ( ') y se manejan como vectores filas. Se direccionan y manipulan igual que los vectores. Son posibles las operaciones matemáticas sobre cadenas. Una vez hecha una operación matemática sobre una cadena, ésta se ve como un vector de números en ASCII. Para ver la representación ASCII de una cadena, podemos utilizar las funciones abs, double o sumamos cero. Para restaurarla y verla de nuevo como cadena de caracteres, usamos la función setstr. Si queremos cambiar a minúsculas añadiremos la diferencia entre 'a' y 'A'. Si queremos que escriba algo en pantalla podemos utilizar el comando disp.
ESTRUCTURAS CÓMO DEFINIRLAS Es una agrupación de datos de tipo diferente bajo un mismo nombre. A los datos les llamamos campos. No hace falta definir previamente el modelo de la estructura, podemos ir creando los distintos campos uno a uno o bien con el comando struct , donde los nombres de los campos se escriben entre apóstrofos ( ') seguidos del valor que se les quiere asignar. Ejemplo >> alumno.nombre = 'Pablo'; % introducimos el campo nombre en la estructura alumno >> alumno.apellido1 = 'Fernández'; % introducimos el campo apellido1 en la estructura alumno >> alumno.apellido2 = 'García'; % introducimos el campo apellido2 en la estructura alumno
>> alumno.edad = 15; % introducimos el campo edad en la estructura alumno >> alumno % escribe por pantalla la información almacenada en la estructura alumno alumno = nombre: 'Pablo' apellido1: 'Fernández' apellido2: 'García' edad: 15 >> alumno2 = struct ('nombre','Fermín','apellido1','Martínez','apellido2','Gil','edad',16) alumno2 = % otro modo de introducir los campos nombre: 'Fermín' apellido1: 'Martínez' apellido2: 'Gil' edad: 16 Pueden crearse vectores y matrices de estructuras, por ejemplo: >> alumno (1) = struct ('nombre','Pablo','apellido1','fernández','apellido2','García','edad',15); >> alumno (2) = struct ('nombre','Fermín','apellido1','Martínez','apellido2','Gil','edad',16); >> alumno % nos devuelve información sobre los campos que tiene la estructura alumno alumno = 1x2 struct array with fields: Nombre Apellido Edad >> alumno (1) % nos devuelve los datos del primer elemento del vector de la estructura ans = nombre: 'Pablo' apellido1: 'fernández' apellido2: 'García' edad: 15 >> alumno (2) % nos devuelve los datos del segundo elemento del vector de la estructura ans = nombre: 'Fermín' apellido1: 'Martínez' apellido2: 'Gil' edad: 16 Para ver un campo concreto de todos los alumnos bastaría teclear: >> alumno.nombre ans = Pablo ans = Fermín
% escribe los datos de todos los campo nombre de la estructura en orden
GRÁFICAS 2-D La orden plot genera una gráfica. Los argumentos deben ser vectores de la misma longitud. Ejemplo: >> x = [-2 -1 0 1 2 3]; y = [4 1 0 1 4 9]; >> plot (x,y)
Si queremos cambiar la apariencia de la gráfica basta pinchar en el último botón de la barra de herramientas y se abrirán unos cuadros en los laterales que nos permitirán ir haciendo los cambios deseados como darle nombre a los ejes.
La función plot nos permite otras opciones como superponer gráficas sobre los mismos ejes: >> x = [-2 -1 0 1 2 3]; y = [4 1 0 1 4 9]; z = [6 5 3 7 5 2]; >> plot (x,y,x,z)
También podemos usar distintos tipos de líneas para el dibujo de la gráfica: >> plot (x,y,'*')
Además podemos colocar etiquetas o manipular la gráfica: etiqueta sobre el eje X de la gráfica actual: eje Y de la gráfica actual: cabecera de la gráfica actual: especificado por las coordenadas: lo indicamos después con el ratón: rejilla: fija valores máximo y mínimo de los ejes: ) fija que la escala en los ejes sea igual: sea un cuadrado: equal y axis square: de gráfico: ventana de gráfico:
>> xlabel('texto') etiqueta sobre el >> ylabel('texto') título en la >> title('texto') texto en el lugar >> text(x,y, 'texto') texto, el lugar >> gtext('texto') dibujar una >> grid >> axis( [xmin xmax ymin ymax] >> axis equal fija que la gráfica >> axis square desactiva axis >> axis normal abre una ventana >> hold on borra lo que hay en la >> hold off
PROGRAMACIÓN DE MATLAB Matlab es una aplicación que permite programar fácilmente.
SENTENCIA FOR Un bloque for en cada iteración asigna a la variable la columna i-ésima de la expresión y ejecuta las órdenes. En la práctica las expresiones suelen ser del tipo escalar:escalar en cuyo caso las columnas son escalares.
for variable = expresión
…
end
SENTENCIA WHILE Un bloque while ejecuta las órdenes mientras todos los elementos de la expresión sean verdaderos.
while …
end
SENTENCIA IF Un bloque if puede escribirse de varias maneras distintas. Lo que hace es evaluar una expresión lógica y si es cierta ejecuta las órdenes que encuentre antes del end.
if <órdenes evaluadas si la expresión es verdadera> end Puede que nos interese que en caso de no ejecutar dicha orden ejecute otra distinta. Esto se lo indicaremos usando else dentro del bloque.
if <órdenes evaluadas si la expresión es verdadera> else <órdenes evaluadas si la expresión end
es
falsa
Si queremos dar una estructura mucho más completa, usaremos la más general donde sólo se evalúan las órdenes asociadas con la primera expresión verdadera de todas. En cuanto la evalúe deja de leer el resto y se dirige directamente al end.
if <órdenes evaluadas si la expresión1 es verdadera> elseif <órdenes evaluadas si la expresión2 es verdadera> elseif <órdenes evaluadas si la expresión3 es verdadera> elseif …… ……
else end <órdenes evaluadas si ninguna otra expresión es verdadera>
SENTENCIA BREAK Si queremos que en un momento dado termine la ejecución de un bucle for o un bucle while usaremos break. SENTENCIA CONTINUE La sentencia continue hace que se pase inmediatamente a la siguiente iteración del bucle for o del bucle while saltando todas las órdenes que hay entre el continue y el fin del bucle en la iteración actual. Ejemplo: Podemos mezclar en un programa varias sentencias de este estilo. Aquí podemos ver un programa que escribe por pantalla los primos del 1 al 100 usando las sentencias if, while y for. disp('Estos son los números primos menores de 100') disp(2) for i=2:100 n=2; while n <= sqrt(i) if rem(i,n)==0 n=i; else n=n+1; end end end if n~=i disp(i) end
ANÁLISIS NUMÉRICO REPRESENTACIÓN GRÁFICA Existe la función fplot que evalúa la función que se desea representar en la gráfica de salida. Como entrada, necesita conocer el nombre de la función como una cadena de caracteres y el rango de representación como un vector de dos elementos: fplot ('nombre', [ valor min, valor max] ) . Ejemplo: >> fplot ('sin', [-3*pi,3*pi] )
1.-ECUACIONES TERMODINAMICAS PARA GASES REALES 1.1.-ECUACION DE VAN DER WALLS
La ecuación de Van der Waals es una ecuación de estado de un fluido compuesto de partículas con un tamaño no despreciable y con fuerzas intermoleculares, como las fuerzas de Van der Waals. La ecuación, cuyo origen se remonta a 1873, debe su nombre a Johannes Diderik van der Waals, quien recibió el premio Nobel en 1910 por su trabajo en la ecuación de estado para gases y líquidos, la cual está basada en una modificación de la ley de los gases ideales para que se aproxime de manera más precisa al comportamiento de los gases reales al tener en cuenta su tamaño no nulo y la atracción entre sus partículas.
1.2.-ECUACION DE REDLICH-KWONG
Introducida en 1949, la ecuación de Redlich-Kwong fue una mejora considerable sobre las otras ecuaciones de la época. Aún goza de bastante interés debido a su expresión relativamente simple. Aunque es mejor que la ecuación de Van der Waals, no da buenos resultados sobre la fase líquida y por ello no puede usarse para calcular precisamente los equilibrios líquido-vapor. Sin embargo, puede usarse conjuntamente con expresiones concretas para la fase líquida en tal caso. La ecuación de Redlich-Kwong es adecuada para calcular las propiedades de la fase gaseosa cuando el cociente entre la presión y la presión crítica es menor que la mitad del cociente entre la temperatura y la temperatura crítica.
1.3.-ECUACION DE SOAVE- REDLICH-KWONG 2
En 1972 Soave reemplazó el término
a/√(T )
de la ecuación de Redlich-Kwong por una expresión
α(T,ω) función de la temperatura y del factor acéntrico. La función α fue concebida para cuadrar con los datos de las presiones de vapor de los hidrocarburos; esta ecuación describe acertadamente el comportamiento de estas sustancias.
1.4.-ECUACION DE PENG-ROBINSON
La ecuación de Peng-Robinson fue desarrollada en 1976 para cumplir los siguientes objetivos:3 1. Los parámetros habían de poder ser expresados en función de las propiedades críticas y el factor acéntrico. 2. El modelo debía ser razonablemente preciso cerca del punto crítico, particularmente para cálculos del factor de compresibilidad y la densidad líquida. 3. Las reglas de mezclado no debían emplear más que un parámetro sobre las interacciones binarias, que debía ser independiente de la presión, temperatura y composición. 4. La ecuación debía ser aplicable a todos los cálculos de todas las propiedades de los fluidos en procesos naturales de gases. Generalmente la ecuación de Peng-Robinson da unos resultados similares a la de Soave, aunque es bastante mejor para predecir las densidades de muchos compuestos en fase líquida, especialmente los apolares.
1.5.-ECUACION DE BEATLIE BRIDGEMAN
Este es un modelo de 5 constantes, cuyas ecuaciones son las siguientes
Por consiguiente, las 5 constantes son 1.6.-ECUACION VIRIAL
Aunque generalmente no es la ecuación de estado más conveniente, la ecuación del Virial es importante dado que puede ser obtenida directamente por mecánica estadística. Si se hacen las suposiciones apropiadas sobre la forma matemática de las fuerzas intermoleculares, se pueden desarrollar expresiones teóricas para cada uno de los coeficientes. En este caso
B
corresponde a
interacciones entre pares de moléculas, C a grupos de tres, y así sucesivamente
1.7.-ECUACION DE BENEDICT-WEBB-RUBIN Esta modelización realizada en 1940 especialmente para hidrocarburos livianos y las mezclas de los mismos también es denominada ecuación BWR.
2.-LISTADO DE PROGRAMAS
EL PROGRAMA CALCULA P,T,V, PARA CADA UNA DE LAS ECUACIONES 2.1.- PROGRAMA PARA LA ECUACION DE VAN DER WAALS
%Utilizando un menu para seleccionar un método para la solucion de %ecuaciones no lineales. clear, clc metodo=menu( 'Por el Metodo de Van der Waals que calcular?:' ,'Volumen','Presion','Temperatura'); syms x switch metodo case 1 %Calculo del Volumen Volar disp('Calculo del Volumen Molar de un gas' ) disp('Por la ecuacion de Van der Waals ' ) disp('Universidad Autonoma Juan Misael Saracho' ) disp('Carrera de Ing. Petroleo y Gas Natural' ) Pc=input('Introduzca la Presion critica del gas (atm):' ); R=input('Introduzca la constante R (atm*L/mol*ºK):' ); Tc=input('Introduzca la Temperatura Critica del gas (ºK):' ); T=input('Introduzca la Temperaturadel gas (ºK):' ); P=input('Introduzca la Presion del gas (atm):' ); a=0.421875*R*R*Tc*Tc/Pc; b=0.125*R*Tc/Pc; fprintf('La solucion de la constate a (atm^2*L^2*K^2/mol^2) es:%f\n' ,a); fprintf('La solucion de la constante b (L/mol) es:%f\n' ,b); disp('La ecuacion de Van der Waals' ) disp('Encontraremos el volumen molar' ) disp('Encontrando el Vo') Vo=R*T/P; fprintf('la solucion de Vo (L/mol) es:%f\n' ,Vo); disp('Por el metodo de Newton-Rapson' ) %metodo newton Rapson para calculo del Volumen Molar syms x fx=P*x^3-(P*b+R*T)*x^2+a*x-a*b; df=diff(fx,x); for i=1:100 y0=subs(fx,Vo); df0=subs(df,Vo); Vn=Vo-(y0)/(df0); fprintf('Iteracion%3d: Vn=%4.4f\n' ,i,Vn) if abs(y0/df0)<10E-9 disp('Por el método de Newton Raphson:' ) fprintf('En %3.0d iteraciones se encontró la convergencia = %4.4f\n' ,i,Vn) return else Vo=Vn; end ezplot(fx),grid on line([-20,0;20,0],[0,-20;0,20], 'color','r') end case 2 %Calculo de la Presion de un gas disp('Calculo de la Presio de un gas' ) disp('Por la ecuacion de Van der Waals ' ) disp('Universidad Autonoma Juan Misael Saracho' ) disp('Carrera de Ing. Petroleo y Gas Natural' ) Pc=input('Introduzca la Presion critica del gas (atm):' ); R=input('Introduzca la constante R (atm*L/mol*ºK):' );
desea
Tc=input('Introduzca la Temperatura Critica del gas (ºK):' ); T=input('Introduzca la Temperaturadel gas (ºK):' ); V=input('Introduzca el volumen del gas (L/mol):' ); a=0.421875*R*R*Tc*Tc/Pc; b=0.125*R*Tc/Pc; fprintf('La solucion de la constate a (atm^2*L^2*K^2/mol^2) es:%f\n' ,a); fprintf('La solucion de la constante b (L/mol) es:%f\n' ,b); P=(R*T/(V-b))-a/V^2; fprintf('La solucion de la presion es:%f\n' ,P); case 3 %Calculo de la temperatura de un gas disp('Calculo de la temperatura de un gas' ) disp('Por la ecuacion de Van der Waals ' ) disp('Universidad Autonoma Juan Misael Saracho' ) disp('Carrera de Ing. Petroleo y Gas Natural' ) Pc=input('Introduzca la Presion critica del gas (atm):' ); R=input('Introduzca la constante R (atm*L/mol*ºK):' ); Tc=input('Introduzca la Temperatura Critica del gas (ºK):' ); P=input('Introduzca la Presion del gas (atm):' ); V=input('Introduzca el volumen del gas (L/mol):' ); a=0.421875*R*R*Tc*Tc/Pc; b=0.125*R*Tc/Pc; fprintf('La solucion de la constate a (atm^2*L^2*K^2/mol^2) es:%f\n' ,a); fprintf('La solucion de la constante b (L/mol) es:%f\n' ,b); T=(P+a/(V^2))*(V-b)*1/R; fprintf('La solucion de la Temperatura (ºK) es:%f\n' ,T); end
2.2 PROGRAMA PARA LA ECUACION DE REDLICH-KWONG
%Utilizando un menu para seleccionar un método para la solucion de %ecuaciones no lineales. clear, clc metodo=menu( 'Por el Metodo de Riedlich-Kwong que calcular?:' ,'Volumen','Presion','Temperatura'); syms x switch metodo case 1 %Calculo del Volumen Molar disp('Calculo del Volumen Molar' ) disp('por el metodo de Riedlich-Kwong' ) disp('Universidad Autonoma Juan Misael Saracho' ) disp('Carrera de Ing. Petroleo y Gas Natural' ) Pc=input('Introduzca la Presion critica del gas (atm):' ); R=input('Introduzca la constante R (atm*L/mol*ºK):' ); Tc=input('Introduzca la Temperatura Critica del gas (ºK):' ); a=0.4278*(R^2)*(Tc^2.5)/(Pc); b=0.0867*R*Tc/Pc; fprintf('la solucion de la constate a (atm*L^2*K^0.5/mol^2) es:%f\n' ,a); fprintf('la solucion de la constante b (L/mol) es:%f\n' ,b); disp('Por la ecuacion de Riedlich-Kwong' ) P=input('Introduzca la Presion del gas (atm):' ); T=input('Introduzca la Temperaturadel gas (ºK):' ); disp('Encontraremos el Volumen Molar' ) disp('Encontrando el Vo') Vo=R*T/P; fprintf('la solucion de Vo (L/mol) es:%f\n' ,Vo); disp('Por el metodo de Newton-Rapson' ) %metodo newton Rapson para calculo del Volumen Molar syms x fx=P*(x^3)-R*T*(x^2)-(b^2)*P*x-R*b*T*x+(a*x)/(T^0.5)-(a*b)/(T^0.5);
desea
df=diff(fx,x); for i=1:100 y0=subs(fx,Vo); df0=subs(df,Vo); Vn=Vo-(y0)/(df0); fprintf('Iteracion%3d: Vn=%4.4f\n' ,i,Vn) if abs(y0/df0)<10E-9 disp('Por el método de Newton Raphson:' ) fprintf('En %3.0d iteraciones se encontró la convergencia = %4.4f\n' ,i,Vn) return else Vo=Vn; end ezplot(fx),grid on line([-20,0;20,0],[0,-20;0,20], 'color','r') end case 2 %Calculo de la Presion del gas disp('Calculo de la Presion del gas' ) disp('por el metodo de Riedlich-Kwong' ) disp('Universidad Autonoma Juan Misael Saracho' ) disp('Carrera de Ing. Petroleo y Gas Natural' ) Pc=input('Introduzca la Presion critica del gas (atm):' ); R=input('Introduzca la constante R (atm*L/mol*ºK):' ); Tc=input('Introduzca la Temperatura Critica del gas (ºK):' ); a=0.4278*(R^2)*(Tc^2.5)/(Pc); b=0.0867*R*Tc/Pc; fprintf('la solucion de la constate a (atm*L^2*K^0.5/mol^2) es:%f\n' ,a); fprintf('la solucion de la constante b (L/mol) es:%f\n' ,b); disp('Por la ecuacion de Riedlich-Kwong' ) V=input('Introduzca del volumen molar (L/mol):' ); T=input('Introduzca la Temperatura del gas (ºK):' ); A=R*T; U=V-b; Y=(V+b)*V*T^0.5; fprintf('la solucion A es:%f\n' ,A); fprintf('la solucion Y es:%f\n' ,Y); fprintf('la solucion U es:%f\n' ,U); P=(A/U)-(a/Y); fprintf('la solucion de la Presion del gas (atm) es:%f\n' ,P); case 3 %Calculo de la Temperatura del gas disp('Calculo de la Tempreratura del gas' ) disp('por el metodo de Riedlich-Kwong' ) disp('Universidad Autonoma Juan Misael Saracho' ) disp('Carrera de Ing. Petroleo y Gas Natural' ) Pc=input('Introduzca la Presion critica del gas (atm):' ); R=input('Introduzca la constante R (atm*L/mol*ºK):' ); Tc=input('Introduzca la Temperatura Critica del gas (ºK):' ); a=0.4278*(R^2)*(Tc^2.5)/(Pc); b=0.0867*R*Tc/Pc; fprintf('la solucion de la constate a (atm*L^2*K^0.5/mol^2) es:%f\n' ,a); fprintf('la solucion de la constante b (L/mol) es:%f\n' ,b); disp('Por la ecuacion de Riedlich-Kwong' ) V=input('Introduzca del volumen molar (L/mol):' ); P=input('Introduzca la Presion del gas (atm):' ); disp('Por la ecuacion de Riedlich-Kwong' ) disp('Encontraremos la Temperatura' ) disp('Encontrando el To') To=V*P/R; fprintf('la solucion de To (ºK) es:%f\n' ,To); disp('Por el metodo de Newton-Rapson' ) %metodo newton Rapson para calculo de la temperatura
syms x fx=P*(V^3)-R*x*(V^2)-(b^2)*P*V-R*b*x*V+(a*V)/(x^0.5)-(a*b)/(x^0.5); df=diff(fx,x); for i=1:100 y0=subs(fx,To); df0=subs(df,To); Tn=To-(y0)/(df0); fprintf('Iteracion%3d: Vn=%4.4f\n' ,i,Tn) if abs(y0/df0)<10E-9 disp('Por el método de Newton Raphson:' ) fprintf('En %3.0d iteraciones se encontró la convergencia = %4.4f\n' ,i,Tn) return else To=Tn; end ezplot(fx),grid on line([-20,0;20,0],[0,-20;0,20], 'color','r') end end
2.3 PROGRAMA PARA LA ECUACION DE SOAVE- REDLICH-KWONG
%Utilizando un menu para seleccionar un método para la solucion de %ecuaciones no lineales. clear, clc metodo=menu( 'Por el Metodo de soave-Riedlich-Kwong calcular?:' ,'Volumen','Presion','Temperatura'); syms x switch metodo case 1 %Calculo del volumen molar disp('Calculo del volumen molar' ) disp('por el metodo de soave-Riedlich-Kwong' ) disp('Universidad Autonoma Juan Misael Saracho' ) disp('Carrera de Ing. Petroleo y Gas Natural' ) Pc=input('Introduzca la Presion critica del gas (atm):' ); R=input('Introduzca la constante R (atm*L/mol*ºK):' ); Tc=input('Introduzca la Temperatura Critica del gas (ºK):' ); T=input('Introduzca la Temperaturadel gas (ºK):' ); a=0.42747*(R^2)*(Tc^2)/(Pc); b=0.086864*R*Tc/Pc; fprintf('la solucion de la constate a (atm*L^2*K/mol^2) es:%f\n' ,a); fprintf('la solucion de la constante b (L/mol) es:%f\n' ,b); w=input('introduzca el valor de w:' ); m=0.4856+1.5517*w-0.1561*w^2; fprintf('la solucion m es:%f\n' ,m); Tr=T/Tc; fprintf('la solucion Tr es:%f\n' ,Tr); alfa=(1+m*(1-Tr^0.5))^2; fprintf('la solucion de alfa es:%f\n' ,alfa); disp('Por la ecuacion de soave-Riedlich-Kwong' ) P=input('Introduzca la Presion del gas (atm):' ); disp('Encontraremos el volumen molar' ) disp('Encontrando el Vo') Vo=R*T/P; fprintf('la solucion de Vo (L/mol) es:%f\n' ,Vo); disp('Por el metodo de Newton-Rapson' ) %metodo newton Rapson para calculo del Volumen Molar syms x fx=P*(x^3)-R*T*(x^2)-(P*b^2+R*T*b-alfa*a)*x-a*b*alfa; df=diff(fx,x);
que
desea
for i=1:100 y0=subs(fx,Vo); df0=subs(df,Vo); Vn=Vo-(y0)/(df0); fprintf('Iteracion%3d: Vn=%4.4f\n' ,i,Vn) if abs(y0/df0)<10E-9 disp('Por el método de Newton Raphson:' ) fprintf('En %3.0d iteraciones se encontró la convergencia = %4.4f\n' ,i,Vn) return else Vo=Vn; end ezplot(fx),grid on line([-20,0;20,0],[0,-20;0,20], 'color','r') end case 2 %Calculo de la presion disp('Calculo de la presion' ) disp('por el metodo de soave-Riedlich-Kwong' ) disp('Universidad Autonoma Juan Misael Saracho' ) disp('Carrera de Ing. Petroleo y Gas Natural' ) Pc=input('Introduzca la Presion critica del gas (atm):' ); R=input('Introduzca la constante R (atm*L/mol*ºK):' ); Tc=input('Introduzca la Temperatura Critica del gas (ºK):' ); T=input('Introduzca la Temperaturadel gas (ºK):' ); V=input('Introduzca del volumen molar (L/mol):' ); a=0.42747*(R^2)*(Tc^2)/(Pc); b=0.086864*R*Tc/Pc; fprintf('la solucion de la constate a (atm*L^2*K/mol^2) es:%f\n' ,a); fprintf('la solucion de la constante b (L/mol) es:%f\n' ,b); w=input('introduzca el valor de w:' ); m=0.4856+1.5517*w-0.1561*w^2; fprintf('la solucion m es:%f\n' ,m); Tr=T/Tc; fprintf('la solucion Tr es:%f\n' ,Tr); alfa=(1+m*(1-Tr^0.5))^2; fprintf('la solucion de alfa es:%f\n' ,alfa); A=R*T; B=V-b; O=(V^2)+V*b; Y=a*alfa; fprintf('la solucion A es:%f\n' ,A); fprintf('la solucion B es:%f\n' ,B); fprintf('la solucion O es:%f\n' ,O); fprintf('la solucion Y es:%f\n' ,Y); disp('Por la ecuacion de soave-Riedlich-Kwong' ) P=(A/B)-(Y/O); fprintf('la solucion de la Presion del gas (atm) es:%f\n' ,P); case 3 %Calculo de la Temperatura disp('Calculo de la Temperatura' ) disp('por el metodo de soave-Riedlich-Kwong' ) disp('Universidad Autonoma Juan Misael Saracho' ) disp('Carrera de Ing. Petroleo y Gas Natural' ) Pc=input('Introduzca la Presion critica del gas (atm):' ); R=input('Introduzca la constante R (atm*L/mol*ºK):' ); Tc=input('Introduzca la Temperatura Critica del gas (ºK):' ); V=input('Introduzca del volumen molar (L/mol):' ); P=input('Introduzca la Presion del gas (atm):' ); a=0.42747*(R^2)*(Tc^2)/(Pc); b=0.086864*R*Tc/Pc; fprintf('la solucion de la constate a (atm*L^2*K/mol^2) es:%f\n' ,a); fprintf('la solucion de la constante b (L/mol) es:%f\n' ,b);
w=input('introduzca el valor de w:' ); m=0.4856+1.5517*w-0.1561*w^2; fprintf('la solucion m es:%f\n' ,m); Tr=x/Tc; alfa=(1+m*(1-Tr^0.5))^2; disp('Por la ecuacion de soave-Riedlich-Kwong' ) disp('Encontraremos la Temperatura' ) disp('Encontrando el To') To=V*P/R; fprintf('la solucion de To (ºK) es:%f\n' ,To); disp('Por el metodo de Newton-Rapson' ) %metodo newton Rapson para calculo de la temperatura syms x fx=P*(V^3)-R*x*(V^2)-(P*b^2+R*x*b-alfa*a)*V-a*b*alfa; df=diff(fx,x); for i=1:100 y0=subs(fx,To); df0=subs(df,To); Tn=To-(y0)/(df0); fprintf('Iteracion%3d: Vn=%4.4f\n' ,i,Tn) if abs(y0/df0)<10E-9 disp('Por el método de Newton Raphson:' ) fprintf('En %3.0d iteraciones se encontró la convergencia = %4.4f\n' ,i,Tn) return else To=Tn; end ezplot(fx),grid on line([-20,0;20,0],[0,-20;0,20], 'color','r') end end
2.4 PROGRAMA PARA LA ECUACION DE PENG-ROBINSON
%Utilizando un menu para seleccionar un método para la solucion de %ecuaciones no lineales. clear, clc metodo=menu( 'Por el Metodo de Peng Robinson calcular?:' ,'Volumen','Presion','Temperatura'); syms x global fv switch metodo case 1 clc,clear syms v disp ('------------------------' ) fprintf('Metodo de Peng Robinson:\n') disp ('------------------------' ) disp ('--------------------------------------------------------------' fprintf ('Ingrese las condiciones del problema\n' ) disp ('--------------------------------------------------------------' T=input('Ingrese la condicion de la Temperatura:(oK)\n T=' ); disp ('--------------------------------------------------------------' P=input('Ingrese la condicion de la Presion:(atm)\n P=' ); disp ('--------------------------------------------------------------' R=input('Ingrese la Constante Universal R:(Atm*L/oK*mol)\n R=' ); disp ('--------------------------------------------------------------' Tc=input('Ingrese la Temperatura Critica (Tabla):(oK)\n Tc=' ); disp ('--------------------------------------------------------------' Pc=input('Ingrese la Presion Critica (Tabla):(Atm)\n Pc=' ); disp ('--------------------------------------------------------------'
que
) ) ) ) ) ) )
desea
fprintf ('La Temperatura y presion reducida son:\n' ) disp ('--------------------------------------------------------------' ) disp (' Tr: Pr:') Tr=T/Tc; Pr=P/Pc; fprintf('%15.5f%15.5f%\n',Tr,Pr); fprintf ('\n-----------------------------------------------------------' ) fprintf ('\nValor de las variable "a","b" son:\n' ) a=0.45724*((R^2*Tc^2)/Pc); b=0.07780*(R*Tc/Pc); disp ('--------------------------------------------------------------' ) disp (' a [atm*L^2/mol^2]: b [Litros/mol]:' ) fprintf('%20.5f%30.5f%\n',a,b) fprintf ('\n-----------------------------------------------------------' ) w=input('\nIngrese la Constante "w" de difusividad del Gas (Tabla):\n w=' ); disp ('--------------------------------------------------------------' ) m=0.37464+1.54226*w-0.29699*w^2; fprintf('El valor de "m" es:\n' ); disp (' m:') fprintf('%15.5f%\n' ,m); fprintf ('\n-----------------------------------------------------------' ) alpha=(1+m*(1-Tr^0.5))^2; fprintf('\nEl valor de "alpha" es:\n' ); disp (' Alpha'); fprintf('%15.5f%\n' ,alpha); fprintf ('\n-----------------------------------------------------------' ) V0=R*T/P; fprintf('\nLa funcion reducida de "Peng Robinson" es: A*v^3+B*v^2+C*v+D' ); fprintf ('\n-----------------------------------------------------------' ) fv=P*v^3+(P*b-R*T)*v^2+(-3*P*b^2-2*b*R*T+a*alpha)*v+(R*T*b^2-a*b*alpha+P*b^3); df=diff(fv,v); for i=1:100 y0=subs(fv,V0); df0=subs(df,V0); Vn=V0-y0/df0; fprintf('\nIteracion%3d: Vn = %4.6f\n',i,Vn) if abs(y0/df0)<10E-8 disp('Por el método de Newton Raphson:' ) fprintf('En %3.0d iteraciones el volumen Molar encontrado es: \n\n=> Vn [litros/mol]= %4.6f\n',i,Vn) return else V0=Vn; end ezplot(fv),grid on line([-20,0;20,0],[0,-20;0,20], 'color','r') end case 2 disp ('------------------------' ) fprintf('Metodo de Peng Robinson:\n') disp ('------------------------' ) disp ('--------------------------------------------------------------' ) fprintf ('Ingrese las condiciones del problema\n' ) disp ('--------------------------------------------------------------' ) T=input('Ingrese la condicion de la Temperatura:(oK)\n T=' ); disp ('--------------------------------------------------------------' ) V=input('Ingrese el volumen molar:(L/mol)\n Vn=' ); disp ('--------------------------------------------------------------' ) R=input('Ingrese la Constante Universal R:(Atm*L/oK*mol)\n R=' ); disp ('--------------------------------------------------------------' ) Tc=input('Ingrese la Temperatura Critica (Tabla):(oK)\n Tc=' ); disp ('--------------------------------------------------------------' ) Pc=input('Ingrese la Presion Critica (Tabla):(Atm)\n Pc=' );
disp ('--------------------------------------------------------------' ) fprintf ('La Temperatura y presion reducida son:\n' ) disp ('--------------------------------------------------------------' ) disp (' Tr: ') Tr=T/Tc; fprintf('%15.5f%15.5f%\n',Tr); fprintf ('\n-----------------------------------------------------------' ) fprintf ('\nValor de las variable "a","b" son:\n' ) a=0.45724*((R^2*Tc^2)/Pc); b=0.07780*(R*Tc/Pc); disp ('--------------------------------------------------------------' ) disp (' a [atm*L^2/mol^2]: b [Litros/mol]:' ) fprintf('%20.5f%30.5f%\n',a,b) fprintf ('\n-----------------------------------------------------------' ) w=input('\nIngrese la Constante "w" de difusividad del Gas (Tabla):\n w=' ); disp ('--------------------------------------------------------------' ) m=0.37464+1.54226*w-0.29699*w^2; fprintf('El valor de "m" es:\n' ); disp (' m:') fprintf('%15.5f%\n' ,m); fprintf ('\n-----------------------------------------------------------' ) alpha=(1+m*(1-Tr^0.5))^2; fprintf('\nEl valor de "alpha" es:\n' ); disp (' Alpha'); fprintf('%15.5f%\n' ,alpha); fprintf ('\n----------------------------------------------------------' ) fprintf('\nLa ecuacion reducida es: \nP=(((R*T)*(V-b)^-1)-((a*alpha)*(V^2+2*V*bb^2)^-1))' ) fprintf ('\n-----------------------------------------------------------' ) P=(((R*T)*(V-b)^-1)-((a*alpha)*(V^2+2*V*b-b^2)^-1)); fprintf('\nEl valor de la Presion encontrado es: \n P=%5.5f\n' ,P) case 3 clc,clear syms v T fv Tr disp ('------------------------' ) fprintf('Metodo de Peng Robinson:\n') disp ('------------------------' ) disp ('--------------------------------------------------------------' ) fprintf ('Ingrese las condiciones del problema\n' ) disp ('--------------------------------------------------------------' ) v=input('Ingrese la condicion del Volumen Molara:(oK)\n Vn=' ); disp ('--------------------------------------------------------------' ) P=input('Ingrese la condicion de la Presion:(atm)\n P=' ); disp ('--------------------------------------------------------------' ) R=input('Ingrese la Constante Universal R:(Atm*L/oK*mol)\n R=' ); disp ('--------------------------------------------------------------' ) Tc=input('Ingrese la Temperatura Critica (Tabla):(oK)\n Tc=' ); disp ('--------------------------------------------------------------' ) Pc=input('Ingrese la Presion Critica (Tabla):(Atm)\n Pc=' ); disp ('--------------------------------------------------------------' ) fprintf ('La Temperatura y presion reducida son:\n' ) disp ('--------------------------------------------------------------' ) disp (' Tr: ') Pr=P/Pc; fprintf('%15.5f%15.5f%\n',Pr); fprintf ('\n-----------------------------------------------------------' ) fprintf ('\nValor de las variable "a","b" son:\n' ) a=0.45724*((R^2*Tc^2)/Pc); b=0.07780*(R*Tc/Pc); disp ('--------------------------------------------------------------' ) disp (' a [atm*L^2/mol^2]: b [Litros/mol]:' ) fprintf('%20.5f%30.5f%\n',a,b) fprintf ('\n-----------------------------------------------------------' )
w=input('\nIngrese la Constante "w" de difusividad del Gas (Tabla):\n w=' ); disp ('--------------------------------------------------------------' ) m=0.37464+1.54226*w-0.29699*w^2; fprintf('El valor de "m" es:\n' ); disp (' m:') fprintf('%15.5f%\n' ,m); fprintf ('\n-----------------------------------------------------------' ) T0=P*v/R; Tr=T/Tc; alpha=(1+m*(1-Tr^0.5))^2; fprintf('\nLa funcion reducida de "Peng Robinson" es: A*v^3+B*v^2+C*v+D' ); fprintf ('\n---------------------------------------------------------\n' ) fv=P*v^3+(P*b-R*T)*v^2+(-3*P*b^2-2*b*R*T+a*alpha)*v+(R*T*b^2-a*b*alpha+P*b^3); df=diff(fv,T); for i=1:100 y0=subs(fv,T0); df0=subs(df,T0); Tn=T0-y0/df0; fprintf('\nIteracion%3d: Tn = %4.6f\n',i,Tn) if abs(y0/df0)<10E-8 fprintf ('\n---------------------------------------------------------\n' ) disp('Por el método de Newton Raphson:' ) fprintf ('-----------------------------------------------------------\n' ) fprintf('En %3.0d iteraciones la temperatura encontrado es: \n\n=> Tn [oK]= %4.6f\n',i,Tn) return else T0=Tn; end ezplot(fv),grid on line([-20,0;20,0],[0,-20;0,20], 'color','r') end end
2.5 PROGRAMA PARA LA ECUACION DE BEATTIE-BRIDGEMAN
%Utilizando un menu para seleccionar un método para la solucion de %ecuaciones no lineales. clear, clc metodo=menu( 'Por el Metodo de BeatleBrigman calcular?:' ,'Volumen','Presion','Temperatura'); syms x switch metodo case 1 %Calculo del Volumen Molar disp('Calculo del Volumen Molar' ) disp('por el metodo de Btatle-Bridgeman ' ) disp('Universidad Autonoma Juan Misael Saracho' ) disp('Carrera de Ing. Petroleo y Gas Natural' ) Ao=input('Introduzca la constante Ao:' ); R=input('Introduzca la constante R:' ); Bo=input('Introduzca la constante Bo del gas:' ); C=input('Introduzca la constante C del gas:' ); T=input('Introduzca la Temperatura del gas :' ); a=input('Introduzca la constante a del gas :'); b=input('Introduzca la constante b del gas :' ); P=input('Introduzca la Presion del gas :' ); disp('Por la Btatle-Bridgeman ' ) disp('Encontraremos el volumen molar' ) disp('Encontrando el Vo') Vo=R*T/P;
que
desea
fprintf('La solucion de Vo (ft^3/lbm) es:%f\n' ,Vo); disp('Por el metodo de Newton-Rapson' ) %metodo newton Rapson para calculo del Volumen Molar syms x fx=R*T*(x+Bo*(1-b/x))*(1-C/(x*T^3))-Ao*(1-a/x)-P*x^2; df=diff(fx,x); for i=1:100 y0=subs(fx,Vo); df0=subs(df,Vo); Vn=Vo-(y0)/(df0); fprintf('Iteracion%3d: Vn=%4.4f\n' ,i,Vn) if abs(y0/df0)<10E-9 disp('Por el método de Newton Raphson:' ) fprintf('En %3.0d iteraciones se encontró la convergencia = %4.4f\n' ,i,Vn) return else Vo=Vn; end ezplot(fx),grid on line([-20,0;20,0],[0,-20;0,20], 'color','r') end case 2 %Calculo de la Presion de un gas disp('Calculo de la Presion de un gas' ) disp('por el metodo de Btatle-Bridgeman ' ) disp('Universidad Autonoma Juan Misael Saracho' ) disp('Carrera de Ing. Petroleo y Gas Natural' ) Ao=input('Introduzca la constante Ao :' ); R=input('Introduzca la constante R :' ); Bo=input('Introduzca la constante Bo del gas :' ); C=input('Introduzca la constante C del gas :' ); T=input('Introduzca la Temperatura del gas :' ); a=input('Introduzca la constante a del gas :' ); b=input('Introduzca la constante b del gas :' ); Vm=input('Introduzca el Volumen Molar del gas :' ); B=(R*T*Bo-Ao-R*C/T^2)/Vm^2; C=(Ao*a-R*T*Bo*b-R*Bo*C/T^2)/Vm^3; D=(R*Bo*C*b/T^2)/Vm^4; K=R*T/Vm; fprintf('La solucion de K es:%f\n' ,K); fprintf('La solucion de B es:%f\n' ,B); fprintf('La solucion de C es:%f\n' ,C); fprintf('La solucion de D es:%f\n' ,D); disp('Por la Btatle-Bridgeman ' ) P=K+B+C+D; fprintf('La solucion la Presion del gas es:%f\n' ,P); case 3 %Calculo de la Temperatura de un gas disp('Calculo de la Temperatura de un gas' ) disp('por el metodo de Btatle-Bridgeman ' ) disp('Universidad Autonoma Juan Misael Saracho' ) disp('Carrera de Ing. Petroleo y Gas Natural' ) Ao=input('Introduzca la constante Ao :' ); R=input('Introduzca la constante R :' ); Bo=input('Introduzca la constante Bo del gas :' ); C=input('Introduzca la constante C del gas :' ); Vm=input('Introduzca el Volumen Molar del gas :' ); a=input('Introduzca la constante a del gas :' ); b=input('Introduzca la constante b del gas :' ); P=input('Introduzca la Presion del gas :' ); disp('Por la Btatle-Bridgeman ' ) disp('Encontraremos la temperatura de un gas' ) disp('Encontrando el To')
To=P*Vm/R; fprintf('La solucion de To (ºR) es:%f\n' ,To); disp('Por el metodo de Newton-Rapson' ) %metodo newton Rapson para calculo de la Temperatura syms x fx=R*x*(Vm+Bo*(1-b/Vm))*(1-C/(Vm*x^3))-Ao*(1-a/Vm)-P*Vm^2; df=duiff(fx,x); for i=1:100 y0=subs(fx,To); df0=subs(df,To); Tn=To-(y0)/(df0); fprintf('Iteracion%3d: Vn=%4.4f\n' ,i,Tn) if abs(y0/df0)<10E-9 disp('Por el método de Newton Raphson:' ) fprintf('En %3.0d iteraciones se encontró la convergencia = %4.4f\n' ,i,Tn) return else To=Tn; end ezplot(fx),grid on line([-20,0;20,0],[0,-20;0,20], 'color','r') end end 2.6 PROGRAMA PARA LA ECUACION VIRIAL
%ECUASION VIRIAL %Utilizando un menu para seleccionar un método para la solucion de la ecuasion virial % Las unidads con las que se esta trabajando son: atm,lt,mol,K clear, clc metodo=menu( 'Seleccione un método de solución desde el menu:','Temperatura','volumen molar','presion'); syms x switch metodo case 1 v=input ('introduzca el valor del volumen: ' ); p=input ('introduzca el valor de la presion:' ); N=input ('introduzca el número de iteraciones:' ); R=0.082054; x0=p*v/R; A0=input('introduzca el valor de A0: ' ); B0=input('introduzca el valor de B0: ' ); b=input('introduzca el valor de b: ' ); a=input('introduzca el valor de a: ' ); c=input('introduzca el valor de c: ' ); for k= 1:N f=(((p*(v)^4)-(R*(x0)*(v)^3)-((R*(x0)*B0-((R*c)/(x0)^2)-A0)*(v)^2)-((A0*a(R*B0*c/(x0)^2)-R*(x0)*B0*b)*v)-((R)^2*B0*b*c/(x0)))); df=(-(R*(v)^3)-((R*B0+((2*R*c)/(x0)^3))*(v)^2)-(((2*R*B0*c/(x0)^3)R*B0*b)*v)+((R)^2*B0*b*c/(x0)^2)); T=x0-(f/df); x0=T fprintf(1, 'T(%d)=%f/n' ,k,T); end fprintf( 'la temperatura calculada es: %2f\n' ,T); case 2 p=input ('introduzca el valor de la presión: ' ); T=input ('introduzca el valor de la temperatura:' ); N=input ('introduzca el número de iteraciones:' ); R= 0.082054; x0=R*T/p;
A0=input('introduzca el valor de A0: ' ); B0=input('introduzca el valor de B0: ' ); b=input('introduzca el valor de b: ' ); a=input('introduzca el valor de a: ' ); c=input('introduzca el valor de c: ' ); B=(R*T*B0)-((R*c)/(T^2))-A0; C=(A0*a)-((R*B0*c)/(T^2))- (R*T*B0*b); D=((R*B0*b*c)/(T^2)); fprintf('el valor de B es: %2f\n' ,B); fprintf('el valor de C es: %2f\n' ,C); fprintf('el valor de D es: %2f\n' ,D); for k= 1:N f=(((p*(x0)^4)-(R*T*(x0)^3)-(B*(x0)^2)-(C*x0)-(D*R*T))); df=((p*4*(x0)^3)-(3*R*T*(x0)^2)-(2*B*(x0))-C); v=x0-(f/df); x0=v fprintf(1, 'v(%d)=%f/n' ,k,v); end fprintf( 'el volumen molar calculado es: %2f\n' ,v); case 3 v=input ('introduzca el valor del volumen: ' ); T=input ('introduzca el valor de la temperatura:' ); N=input ('introduzca el número de iteraciones:' ); R= 0.082054; x0=R*T/v; A0=input('introduzca el valor de A0: ' ); B0=input('introduzca el valor de B0: ' ); b=input('introduzca el valor de b: ' ); a=input('introduzca el valor de a: ' ); c=input('introduzca el valor de c: ' ); B=(R*T*B0)-((R*c)/(T^2))-A0; C=(A0*a)-((R*B0*c)/(T^2))- R*T*B0*b; D=((R*B0*b*c)/(T^2)); fprintf('el valor de B es: %2f\n' ,B); fprintf('el valor de C es: %2f\n' ,C); fprintf('el valor de D es: %2f\n' ,D); for k= 1:N f=(((x0*(v)^4)-(R*T*(v)^3)-(B*(v)^2)-(C*v)-(D*R*T))); df=((v)^4); p=x0-(f/df); x0=p fprintf(1, 'p(%d)=%f/n' ,k,p); end fprintf( 'la presion calculada es: %2f\n' ,p); end
2.7 PROGRAMA PARA LA ECUACION DE BENEDICT-WEBB-RUBIN
%Utilizando un menu para seleccionar un método para la solucion de %ecuaciones no lineales. clear, clc metodo=menu( 'Por el Metodo de Benedict-Webb-Rubin calcular?:' ,'Volumen','Presion','Temperatura'); switch metodo case 1 clc,clear %Calculo del Volumen Molar por la ecuacion de Benedict-Webb-Rubi disp('Calculo del Volumen Molar' ) disp('por el metodo de Benedict-Webb-Rubi' ) disp('Universidad Autonoma Juan Misael Saracho' ) disp('Carrera de Ing. Petroleo y Gas Natural' )
que
desea
Ao=input('Introduzca la constante Ao :' ); R=input('Introduzca la constante R :' ); Bo=input('Introduzca la constante Bo del gas :' ); Co=input('Introduzca la constante Co del gas :' ); alfa=input( 'Introduzca alfa:'); gama=input( 'Introduzca gama:'); c=input('Introduzca la constante c:' ); T=input('Introduzca la Temperatura del gas (ºK):' ); a=input('Introduzca la constante a del gas:' ); b=input('Introduzca la constante b del gas:' ); P=input('Introduzca la Presion del gas:' ); disp('Por la Benedict-Webb-Rubi') disp('Encontraremos el volumen molar' ) disp('Encontrando el Vo') Vo=R*T/P; fprintf('La solucion de Vo es:%f\n' ,Vo); syms x A=R*T; B=(Bo*R*T-Ao-(Co/T^2)); D=(b*R*T-a); F=a*alfa; J=(c/T^2); fprintf('La solucion de A es:%f\n' ,A); fprintf('La solucion de B es:%f\n' ,B); fprintf('La solucion de D es:%f\n' ,D); fprintf('La solucion de F es:%f\n' ,F); fprintf('La solucion de J es:%f\n' ,J); disp('Por el metodo de Newton-Rapson' ) %metodo newton Rapson para calculo del Volumen Molar syms x fx=-P+(A*1/x)+(B*1/x^2)+(D*1/x^3)+(F*1/x^6)+(J*(1+gama/(x^2))*exp((gama/x^2)))*(1/x^3); df=diff(fx,x); for i=1:100 y0=subs(fx,Vo); df0=subs(df,Vo); Vn=Vo-(y0)/(df0); fprintf('Iteracion%3d: Vn=%4.4f\n' ,i,Vn) if abs(y0/df0)<10E-9 disp('Por el método de Newton Raphson:' ) fprintf('En %3.0d iteraciones se encontró la convergencia Vn = %4.4f\n' ,i,Vn) return else Vo=Vn; end ezplot(fx),grid on line([-20,0;20,0],[0,-20;0,20], 'color','r') end case 2 %Calculo del Volumen Molar por la ecuacion de Benedict-Webb-Rubi disp('Calculo del Volumen Molar' ) disp('por el metodo de Benedict-Webb-Rubi' ) disp('Universidad Autonoma Juan Misael Saracho' ) disp('Carrera de Ing. Petroleo y Gas Natural' ) Ao=input('Introduzca la constante Ao :' ); R=input('Introduzca la constante R :' ); Bo=input('Introduzca la constante Bo del gas :' ); Co=input('Introduzca la constante Co del gas :' ); alfa=input( 'Introduzca alfa:'); gama=input( 'Introduzca gama:'); c=input('Introduzca la constante c:' ); T=input('Introduzca la Temperatura del gas (ºK):' ); a=input('Introduzca la constante a del gas:' );
b=input('Introduzca la constante b del gas:' ); Vm=input('Introduzca el Volumen Molar del gas:' ); disp('Por la Benedict-Webb-Rubi') A=R*T; B=(Bo*R*T-Ao-(Co/T^2)); D=(b*R*T-a); F=a*alfa; J=(c/T^2); fprintf('La solucion de A es:%f\n' ,A); fprintf('La solucion de B es:%f\n' ,B); fprintf('La solucion de D es:%f\n' ,D); fprintf('La solucion de F es:%f\n' ,F); fprintf('La solucion de J es:%f\n' ,J); disp('Por el metodo de Newton-Rapson' ) P=(A*1/Vm)+(B*1/Vm^2)+(D*1/Vm^3)+(F*1/Vm^6)+(J*(1+gama/(Vm^2))*exp((gama/Vm^2)))*(1/Vm^3); fprintf('La solucion de la Presion del gas (bar) es:%f\n' ,P); case 3 %Calculo de la Temperatura de un gas por la ecuacion de Benedict-Webb-Rubin disp('Calculo de la Tempratura de un gas' ) disp('por el metodo de Benedict-Webb-Rubi' ) disp('Universidad Autonoma Juan Misael Saracho' ) disp('Carrera de Ing. Petroleo y Gas Natural' ) Ao=input('Introduzca la constante Ao :' ); R=input('Introduzca la constante R :' ); Bo=input('Introduzca la constante Bo del gas :' ); Co=input('Introduzca la constante Co del gas :' ); alfa=input( 'Introduzca alfa:'); gama=input( 'Introduzca gama:'); c=input('Introduzca la constante c:' ); Vm=input('Introduzca el Volumen Molar del gas (m^3/Kmol):' ); a=input('Introduzca la constante a del gas:' ); b=input('Introduzca la constante b del gas:' ); P=input('Introduzca la Presion del gas:' ); disp('Por la Benedict-Webb-Rubi') disp('Encontraremos de la Temperatura del gas' ) disp('Encontrando el To') syms x To=P*Vm/R; fprintf('La solucion de Vo es:%f\n' ,To); A=R*x; B=(Bo*R*x-Ao-(Co/x^2)); D=(b*R*x-a); F=a*alfa; J=(c/x^2); fprintf('La solucion de F es:%f\n' ,F); disp('Por el metodo de Newton-Rapson' ) %metodo newton Rapson para calculo de la Temperatura del gas fx=-P+(A*1/Vm)+(B*1/Vm^2)+(D*1/Vm^3)+(F*1/Vm^6)+(J*(1+gama/(Vm^2))*exp((gama/Vm^2)))*(1/Vm^3); df=diff(fx,x); for i=1:100 y0=subs(fx,To); df0=subs(df,To); Tn=To-(y0)/(df0); fprintf('Iteracion%3d: Tn=%4.4f\n' ,i,Tn) if abs(y0/df0)<10E-9 disp('Por el método de Newton Raphson:' ) fprintf('En %3.0d iteraciones se encontró la convergencia Vn = %4.4f\n' ,i,Tn) return else To=Tn; endn
ezplot(fx),grid on line([-20,0;20,0],[0,-20;0,20], 'color','r') end end
3.-EJERCICIOS RESUELTOS: 3.1.- Se ha de transportar gas natural desde los yacimientos productores hasta los centros de
consumo, a través de un ducto, el gas es esencialmente metano puro con una presión de 204 atm y Temperatura 18.3ºC. Calcular el Volumen molar.
a) b) c) d) e) f) g) a) b) c) d) e) f)
Redlich-Kwong. Van der Waals. Benedict-Webb-Rubin Beattie-Bridgeman Peng-Robinson Soave- Redlich-Kwong
Respuestas En 5 iteraciones se encontró la convergencia = 0.0972 L/mol En 5 iteraciones se encontró la convergencia = 0.0963 L/mol En 5 iteraciones se encontró la convergencia = 0.0951 m 3/kmol En 4 iteraciones se encontró la convergencia = 0.1018 L/mol En 5 iteraciones el volumen Molar encontrado es: Vn [litros/mol]= 0.092828 En 5 iteraciones se encontró la convergencia = 0.0996 L/mol
3.2.- Se ha de transportar gas natural desde los yacimientos productores hasta los centros de consumo, a través de un ducto, el gas es esencialmente etano puro con un Volumen molar 0.096 L/mol y Temperatura 18.3ºC. Calcular la presión del gas.
a) Redlich-Kwong. b) Van der Waals. Respuestas: a) la solucion de la Presion del gas (atm) es:48.223172 b) La solucion de la presion es:173.960051 atm 3.3.- Se ha de transportar gas natural desde los yacimientos productores hasta los centros de consumo, a través de un ducto, el gas es esencialmente etano puro con un Volumen molar 0.098 L/mol y Presion 200 atm. Calcular la temperatura del gas. a) Peng-Robinson b) Soave- Redlich-Kwong
Respuestas: a) En 4 iteraciones la temperatura encontrado es: Tn [K]= 374.860321 b) En 5 iteraciones se encontró la convergencia = 358.4484 K
4.-PROBLEMAS PROPUESTOS:
1.-Determine el volumen molar del vapor de agua a 20MPa y 400°C. Utilizando la ecuación de: a) Redlich-Kwong. b) Van der Waals. c) Benedict-Webb-Rubin d) Beattie-Bridgeman e) Peng-Robinson f) Soave- Redlich-Kwong 2.- Utilizando la ecuación de: 1. Redlich-Kwong. 2. Van der Waals. 3. Benedict-Webb-Rubin 4. Beattie-Bridgeman 5. Peng-Robinson 6. Soave- Redlich-Kwong Determine el volumen, en m3, que ocupan 165 kg de metano con una presión de 200 atm y temperatura de 400 K. 3.- Utilizando la ecuación de estado de: a) Redlich-Kwong. b) Van der Waals. c) Benedict-Webb-Rubin d) Beattie-Bridgeman e) Peng-Robinson f) Soave- Redlich-Kwong Determine la presión del nitrógeno gas (N2), en MPa, a 193 K y volumen molar de 4,5 cm3/mol 4.- Utilizando las ecuaciones de: 1. 2. 3. 4. 5. 6.
Redlich-Kwong. Van der Waals. Benedict-Webb-Rubin Beattie-Bridgeman Peng-Robinson Soave- Redlich-Kwong a) El volumen molar del metano a una temperatura de 680ºC y una presión de
8 atm
b) La temperatura del metano, a un volumen molar de 900 L/mol y una presión de 101,325kpa c) La presión del metano a una temperatura 700ºF y un volumen molar de 600m3/mol
5.-Utilizando la ecuación de:
a) b) c) d) e) f)
Redlich-Kwong. Van der Waals. Benedict-Webb-Rubin Beattie-Bridgeman Peng-Robinson Soave- Redlich-Kwong 1. La presión del etano a una temperatura de 1000ºC y un volumen molar de 680m3/mol 2. El volumen molar del oxígeno a condiciones estándar de presión y te mperatura 3. La temperatura del nitrógeno a una presión de 1000 torr y un volumen molar de 1500L/mol.
6.-Utilizando la ecuacion de: a) Redlich-Kwong. b) Van der Waals. c) Benedict-Webb-Rubin d) Beattie-Bridgeman e) Peng-Robinson f) Soave- Redlich-Kwong Determine la presión, en atm, a la temperatura de 100°C si el volumen molar es 0,012 m3 /Kmol CO2. 7.- Calcular el volumen que ocupa un mol de oxígeno a 100 atm y 298 K. a) b) c) d) e) f)
Redlich-Kwong. Van der Waals. Benedict-Webb-Rubin Beattie-Bridgeman Peng-Robinson Soave- Redlich-Kwong
Bibliografía 1.- TEXTO UNIVERSITARIO INGENIERIA ASISTIDA POR COMPUTADORA
Autor: Ing. Salomón Benito Huanca
2.- http://es.wikipedia.org/wiki/Ecuaci%C3%B3n_de_estado 3.-Termodinámica Técnica - Moran Shapiro II 4.-TERMODINAMICA ED6 CENGEL 5.-Introducción a la Termodinámica en Ingeniería Química - Smith, Van Ness
TABLAS PARA LOS GASES REALES
TABLAS PARA LOS GASES REALES
Constante R = 0.0831bar*m 3/Kmol*’K
