DISEÑO DE COLUMNAS DE ABSORCIÓN USANDO MATLAB Prof. FABIAN LEONARDO MORENO
MATLAB •
•
Laboratorio de matrices Opera matemáticamente funciones, vectores , matrices, ecuaciones, ecuaciones diferenciales, etc…
•
•
Sirve para crear programas Importa, exporta datos
MATLAB •
•
Laboratorio de matrices Opera matemáticamente funciones, vectores , matrices, ecuaciones, ecuaciones diferenciales, etc…
•
•
Sirve para crear programas Importa, exporta datos
Ventana de comandos •
Ejecuto acciones: ejemplo.
•
Suma, resta, trigonométricas, pi, date, etc…
•
Trabaja con polinomios: polinomios: p=[1,2,3,4] p=[1,2,3,4]
•
•
•
Soluciona la ecuación polinomial: roots (p) Tambien devuelve el polinomio dadas las raíces: p2= poly([-2,-3])
•
•
•
Resuelva el sistema: 2x + 3y = 4 5x – 2y = 6
•
•
•
•
•
>> a=[2, 3; 5, -2]; >> b=[4; 6]; >> x=inv(a)*b; >> x
ingrese la matriz de coeficientes ingrese el vector columna de constantes obtenga la solución con la función para invertir muestre el vector solución
•
•
Grafique la función f(x)=sen(x) ex en el intervalo 0 x
•
•
•
>> ezplot('sin(x)*exp(x)',0, 2); >> grid on;
%escriba la función entre comillas simples %muestre cuadrículas en el gráfico
•
•
Resuelva la ecuación cúbica 5x3 + 2x2 - 3x + 1 = 0;
•
•
>> a=[5, 2, -3, 1];
ecuación •
>> x=roots(a)
ingrese los coeficientes de la obtenga y muestre las raíces de la ecuación
•
obtenga y muestre las raíces de la ecuación
>> x=roots(a)
•
•
Obtenga la solución de la ecuación diferencial ordinaria : y´-x-y = 0, y(0)=1
•
•
•
•
•
>> y=dsolve('Dy-x-y=0','y(0)=1', 'x'); %defina la ecuación, condición y variable >> y %muestre la solución analítica obtenida >> ezplot(y, 0, 2); % grafique la solución para 0 x 2 >> grid on % muestre cuadrículas
•
•
Integre la función f(x) = x sen(x)
•
•
•
obtenga el resultado analítico integre entre 0 y
>> f=int('x*sin(x)') >> s=int('x*sin(x)',0,pi)
•
•
Manejo simbólico de expresiones
•
•
•
•
•
>> syms x >> y=x^3-8 >> t=factor(y) >> e=taylor(exp(x), 5);
defina x con tipo simbólico una expresión con x factorar la espresión asignada a y expandir ex con 5 términos de la serie de Taylor
graficas •
•
• •
La función plot(x), genera una gráfica de los valores de x que han sido calculados en determinada función Las graficas se pueden editar agregando titulo (tittle), rotulos (xlabel, ylabel), cuadriculas (grid), texto (text) etc… Clear all borra todo Clc limpia pantalla
Ejemplo grafica •
•
•
•
•
•
•
•
•
Graficar la función sin(x)*exp(-0.4*x) entre 0 y 10. Opción 1: ezplot('sin(x)*exp(-0.4*x)',[0,10]) Opción 2: fplot('sin(x)*exp(-0.4*x)',[0,10]) Opción 3: x=0:0.05:10; y=sin(x).*exp(-0.4*x); plot(x,y)
•
•
•
Graficar: Y = 3x^2+2x-5 Entre x = -8 y x= 8
•
•
•
fplot('sin(x)*exp(-0.4*x)',[0,10]); title('ejemplo') grid on;
PROGRAMANDO EN MATLAB Se pueden usar instrucciones en la ventana de comandos o en archivos llamados ficheros .m Bucles: For, if, while, for i=(1:2:10) (valor inicial, incremento, valor final) y=i
•
•
•
•
•
•
•
If expresión 1 hacer algo else if expresión 2 Hacer algo Else Hacer algo end
•
•
•
While expresión Hacer algo end
Ficheros .m •
Son archivos en los que se guardan instrucciones y comandos. Es como crear un programa que se puede llamar en cualquier momento.
Funciones •
•
Son ficheros .m que devuelven un valor de una función dada. Function[salida1, salida2, …] = nombre( entrada1, entrada2, …)
•
Comandos de MATLAB Return
•
La función se evalua con: feval
•
Ejemplo funciones •
•
•
•
function p=fun1(x) p=x^3—2 Despues en la ventana de comandos se llama la función y el punto donde se quiere evaluar: fun1(3*pi/2)
Ejemplo newton raphson x=-4:0.1:4 fx=x.^3+x+16 plot(x,fx) grid on
format short; x=4; fx=x.^3+x+16; while abs(fx)>0.00001 fx=x.^3+x+16; dfx=3*x^2+1; xn=x-fx/dfx; disp ([x fx dfx fx/dfx xn]); x=xn; end
Otra rutina para Newton Raphson • • • • • • • • • • • •
function [res, it]=newton(func,dfunc,x,presic) %x0 es el valor inicial, precis la presicion requerida %func la función y dfunc es su derivada it=0; x0=0; d=feval(func,x0)/feval(dfunc,x0); while abs(d)>presic x1=x0-d; it=it+1; x0=x1; d=feval(func,x0)/feval(dfunc,x0); end; res=x0
•
•
•
•
•
•
•
Se crean ficheros para las funciones, por ejemplo para la función: x^2-x-sin(x+0.15) Fichero f1 function f=f1(x); f=x^2-x-sin(x+0.15); Fichero derf1 function d=derf1(x); d=2*x-1-cos(x+0.15);
Solución de EDO en Matlab •
•
•
•
OPCIONES: Solución exacta (directa) Solución de una función definida previamente Solución por métodos numéricos programando la rutina respectiva
ejercicio •
Resolver la ecuación:
•
y‘ = (t-y)/2 Y(0) = 1
•
•
•
•
En Matlab en forma directa. En Matlab usando una rutina para el metodo de euler En excel usando una rutina con el metodo de euler y comparando el resultado para diferentes valores de n. Graficar los resultados y compararlos con la solucion real: y(t) = 3e(-t/2) – 2 + t
Solución exacta •
Resolver: y’=(t-y)/2 en [0,3] con y(0)=1
•
Se usa la función dsolve
•
•
•
•
•
Sintaxis en la ventana de comandos: u = dsolve('Du = (t-u)/2','u(0) = 1') Si queremos graficar la solución: t=0:0.1:3 plot((-2+t+3*exp(-1/2*t)))
Solución exacta •
1.
Otro ejemplo Ej: resolver la ecuación de valor inicial: U’ =1/2u ; u(0) = ¼
2. •
•
U' = 16 u(0) = 2
•
•
•
Solución: 1. u = dsolve('Du = u/2','u(0) = 1/4') 2. u = dsolve('Du = r*u','u(0) = u0')
SOLUCION DE ODE EN MATLAB •
•
Matlab ya tiene implicto el metodo de Runge Kutta para resolver las ecuaciones diferenciales ordinarias. La instrucción es ODE45, ODE23, ODE 1213, etc
•
•
•
•
•
•
•
•
•
•
•
•
Ejemplo
Resolver la ecuación: dy/dt =a + by-y^3
%Se define la función en el editor: % Comentarios, instrucciones, etc. function [f]=ecdif(t,y) % Parámetros a=6;b=1.5; % Ecuacion f=a+b*y-y^3; f = f'; Luego se resuelve el sistema en la ventana de comandos: >>tr=1;y0=1;[t,y]=ode45('ecdif',tr,y0); plot (t,y)
Solucion de sistemas de ecuaciones diferenciales •
•
•
Resolver el sistema dx1/dt =a + bX1-x1^3 dx2/dt =cX2+a
•
%Se define la función en el editor: % Comentarios, instrucciones, etc. function [f]=ecdif(t,x) % Parámetros a=6;b=1.5;c=-14; % Ecuaciones f(1)=a+b*x(1)-x(1)^3; f(2)=c*x(2)+a; f = f'; plot (t,x) Luego se resuelve el sistema en la ventana de comandos:
•
>>tr=[0,1];x0=[1,1];[t,x]=ode45('ecdif',tr,x0);
•
•
•
•
•
•
•
•
•
•
EJEMPLO 13.4 THOMPSON PRIMERO SE CREA EL FICHERO .M EN LA PESTAÑA NEW/FUNCTION
M-file
•
•
functiondy=absorber(z,y)
•
•
dy=zeros(6,1); %VALORES INCIALES
•
•
•
•
•
•
•
•
•
•
•
•
•
XAGB=.17; %CONCENTRACIÓN DE A EN GAS ENTRADA BOTTOM XALB=.026; % CONCENTRACIÓN A EN EL AGUA QUE SALE BOTTOM XAGT=.01; % CONCENTRACIÓN DE A EN EL GAS QUE SALE TOP NGB=26.6; % FLUX DE GAS QUE ENTRA BOTTOM NLB=-164.7; % FLUX LÍQUIDO QUE SALE BOTTOM DIRECCIÓN CONTRARIA Y Z POSITIVO (BAJA) DHV=55.6; % ENTALPÍA VAP AGUA DHS=1.3E4; % ENTALPÍA DISOLUCIÓN AMONÍACO EN AGUA CPG=7.0; % CP GAS CPL=18; % CP LÍQUIDO PTOT=1.0; %PTOT SISTEMA TB=35; % TEMPERATURA FONDO TORRE BOTTOM KOGAI=19; %KOG*aI
•
•
%PARÁMETROS DEL EQUILIBRIO AMONIACO
•
•
•
A1=1.6E9; A2=6276;
•
•
%PARÁMETROS ANTOINE
•
•
B1=8.10765; B2=1750.286; B3=235;
•
•
•
•
NXAGB=4.522; %FLUX GAS DE ENTRADA DE A NXALB=-4.28; %FLUX LÍQUIDO A LA SALIDA DE A NXBGB=.146; % FLUX EN EL GAS ENTRADA 10% HUMEDAD RELATIVA DE B
•
•
%FÓRMULAS:
•
•
•
•
XAG=y(1)/y(5); % HACE EXPLÍCITO XAG y(6)=XAG; % HACE EXPLÍCITO XAG NL=-y(5)+NGB+NLB; %ÚLTIMA ECUACIÓN DE LA TABLA 13-2, NO HAY QUE DEFINIR COMO Y(I) PORQUE NO ES DIFERENCIAL
•
•
•
•
•
•
•
•
•
•
MA=A1*exp(-A2/(y(3)+273)); %EQUILIBRIO A XAL=y(2)/NL; XAGS=MA*XAL; %FRACCIÓN DE A EN EL GAS EN EQUILIBRIO CON EL LÍQUIDO PARA CUALQUIER CONCENTRACIÓN EN EL LÍQUIDO NAI=KOGAI*(XAG-XAGS); % CANTIDAD DE A QUE CRUZA LA INTERFASE ALPHA=B1-B2/(B3+y(3)); % EC ANTOINE PBST=10^ALPHA; %PRESIÓN DE VAPOR DE B POR EC.ANTOINE XBGS=PBST/(PTOT*760); % FRACCIÓN DE B EN EL GAS EN EL EQUILIBRIO RAOULT XBG=y(4)/y(5); % HACE EXPLÍCITO XBG NBI=KOGAI*(XBGS-XBG); %FLUX DE B QUE CRUZA LA INTERFASE
•
•
•
•
•
•
•
•
•
% ECUACIONES DIFERENCIALES RESUELTAS SEGÚN TABLA 13-2 dy(1)=-NAI dy(2)=NAI dy(3)=(1/1.8)*(-DHV*NBI+DHS*NAI)/(y(5)*CPG+NL*CPL) dy(4)=NBI dy(5)=NBI-NAI dy(6)=(1/y(5))*(dy(1)-y(6)*dy(5)) end
EJEMPLO 13.4 THOMPSON • • • • • • • • • •
% REALIZA LA GRÁFICA DESPUÉS DE CORRER EL M-FILE clc clear all y0(1)=4.522; y0(2)=-4.28; y0(3)=35; y0(4)=.146; y0(5)=26.6; y0(6)=.17;
• •
[z,y]=ode45('absorber',[0 5],[y0]);
• • • • • •
plotyy(z,y(:,3),z,y(:,6)) ylabel('TEMPERATURE(F)') ylabel('XAG') xlabel('z (ft)') title('{Fracción molar y temperatura en un absorbedor de gas}')
PARA EL PROYECTO: •
•
•
•
Determinar el equilibrio del sistema usando aspen: Ingresar los componentes, abrir el menu tools/analysis/binary Seleccionar el diagrama Txy Generar una correlación en excel son los datos del equilibrio para tener el valor de m
