libro-de-metodos-numericos-en-matlab.pdf.pdf

-------------------------------------------------------------------------------------------------------------------------------------------------------1 Alberto L.Huamaní L.Huamaní Huamaní Huamaní

INDICE


INDICE


Métodos numéricos en Ingeniería Ingeniería de Alimentos

Solución numérica

CAPITULO I

SOLUCION NUMERICA ---------------------------------------------------------------------------------------------------------Objetivos: Desarrollar ejercicios de soluciones numéricas para funciones matemáticas generales y aplicaciones a la ingeniería de alimentos. ---------------------------------------------------------------------------------------------------------1.1 SOLUCION NUMERICA Uno de los problemas más antiguos y básicos del cálculo numérico es el problema de búsqueda de la solución de una ecuación, es decir encontrar los valores de la variable x que satisfacen la ecuación f(x)=0, para una función f dada. Las ecuaciones pueden ser algebraicas (la función f es un polinomio), por ejemplo: x 2+5x-4=0 o bien trascendentes puesto que están constituidas por funciones trascendentes tales como funciones exponenciales, trigonométricas, logarítmicas, etc., por ejemplo: e -x – x; x; sen x; ln x 2 – 1. 1. Solamente en casos muy simples, de ecuaciones algebraicas, existen fórmulas que permiten resolverlas en términos t érminos de sus coeficientes, para el resto de las ecuaciones se utilizan métodos aproximados que permiten mejorar la solución por simple repetición del mismo método hasta adquirir el grado de aproximación requerido. Estos métodos son apropiados para realizarlos utilizando computadoras puesto que comprenden la repetición de un proceso, es decir iteración. it eración. El objetivo principal del análisis numérico es encontrar soluciones “aproximadas” a problemas complejos utilizando sólo las operaciones operaciones más simples de la aritmética. Se requiere de una secuencia de operaciones algebraicas y lógicas que producen la aproximación al problema matemático.

Los métodos numéricos pueden ser aplicados para resolver procedimientos matemáticos en: Cálculo de derivadas, Integrales, Ecuaciones diferenciales, Operaciones con matrices, Interpolaciones, Ajuste de curvas, Polinomios Los métodos numéricos se aplican en áreas como: Ingeniería Industrial, Ingeniería Química, Ingeniería Industrias Alimentarias, Ingeniería Civil, Ingeniería Mecánica, Ingeniería eléctrica, etc... Por ello, los errores se deben:   

Identificar Cuantificar Minimizar

Cálculo Numérico es una materia de Cálculo o Matemáticas Aplicada, Aplicada, que muestra cómo a través de fórmulas e iteraciones podemos obtener resultados bastante bastante aproximados para diversos problemas que se pueden plantear. -------------------------------------------------------------------------------------------------------------------------------------------------------1 Alberto L.Huamaní L.Huamaní Huamaní Huamaní

Métodos numéricos en Ingeniería de Alimentos

Solución numérica

Se deben tener conocimientos de Cálculo Matemático, Series, Algebra Lineal, Aritmética y Trigonometría, entre otras cosas. Los errores en cálculos y medidas se pueden caracterizar con respecto a su exactitud y su precisión. Exactitud: La exactitud se refiere a qué tan cercano está el valor calculado o medido del valor verdadero. Precision: La precisión se refiere a qué tan cercanos se encuentran, unos de otros, diversos valores calculados o medidos. Errores: Los errores numéricos surgen del uso de aproximaciones para representar operaciones y cantidades matemáticas exactas. Éstas incluyen los errores de truncamiento que resultan del empleo de aproximaciones como un procedimiento matemático exacto, y los errores de redondeo que se producen cuando se usan números que tienen un límite de cifras significativas para representar números exactos. Para ambos tipos de errores, la relación entre el resultado exacto, o verdadero, y el aproximado está dada por Valor verdadero = Valor aproximado + error

(1.1)

Reordenando la ecuación (1.1) se encuentra que el error numérico es igual a la diferencia entre el valor verdadero y el valor aproximado, es decir Et = valor verdadero – valor aproximado

(1.2)

Donde Et se usa para denotar el valor exacto del error. El subíndice t indica que se trata del error “verdadero” (true). El error relativo también se puede multiplicar por 100% para expresarlo como  t 

Error verdadero

Valor verd adero

* 100%

(1.3)

Donde et denota el error relativo porcentual verdadero. Error normalizado  a 

Error aproximado

Valor aproximado

* 100%

(1.4)

Donde el subíndice a significa que el error está normalizado a un valor aproximado Por ejemplo, ciertos métodos numéricos usan un método iterativo para calcular los resultados. En tales métodos se hace una aproximación considerando la aproximación anterior. Este proceso se efectúa varias veces, o de forma iterativa, para calcular en forma sucesiva, esperando cada vez mejores aproximaciones. -------------------------------------------------------------------------------------------------------------------------------------------------------2 Alberto L.Huamaní Huamaní


Solución numérica

Error relativo porcentual En tales casos, el error a menudo se calcula como la diferencia entre la aproximación previa y la actual. Por lo tanto, el error relativo porcentual está dado por  a 

Aproximación actual  Aproximación anterior *100% Aproximación actual

(1.5)

Los signos pueden se positivos o negativos, Por lo tanto, es útil emplear el valor absoluto.  s Fijado previamente.

 a   s

(1.6)

Errores de redendeo Los números tales como pi, e ó 7 no pueden exspresarse con un número fijo de cifras significativas. Por lo tanto, no pueden ser representados exactamente por la computadora. Además, debido a que las computadoras usan una representación en base 2, no pueden representar exactamente algunos números en base 10. Esta discrepancia por la omisión de cifras significativas se llama error de redondeo. Errores de truncamiento Los errores de truncamiento son aquellos que resultan al usar una aproximación en lugar de un procedimiento matemático exacto. Por ejemplo, en el capítulo 1 aproximamos la derivada de la velocidad de caída de un paracaidista mediante una ecuación en diferencia finita dividida de la forma

 a   s

(1.7)

Metodos cerrados En los métodos cerrados, la raíz se encuentra dentro de un intervalo predeterminado por un límite inferior y otro superior. La aplicación repetida de estos métodos siempre genera aproximaciones cada vez más cercanas a la raíz. Se dice que tales métodos son convergentes porque se acercan progresivamente a la raíz a medida que se avanza en el cálculo. Metodos abiertos Los métodos abiertos se basan en fórmulas que requieren únicamente de un solo valor de inicio x o que empiecen con un par de ellos, pero que no necesariamente encierran la raíz. Éstos, algunas veces divergen o se alejan de la raíz verdadera a medida que se avanza en el cálculo. Sin embargo, cuando los métodos abiertos convergen, en general lo hacen mucho más rápido que los métodos cerrados. -------------------------------------------------------------------------------------------------------------------------------------------------------3 Alberto L.Huamaní Huamaní


Solución numérica

1.1.1 Procedimiento de solución usando programa en interfaz GUIDE Matlab 1. 2. 3. 4. 5. I.

Crear una carpeta con el nombre Raíces Guardar una imagen en jpg, dentro de la carpeta (con nombre caratula u otro) Abrir el software Matlab Hacer click en la ventana New Script Digitar el siguiente código y guardar en archivo.m (Caso para la caratula)

Programa para la CARATULA

function caratula %Autor: Ing. Alberto Luis HUAMANI HUAMANI %*************************************************************** % presentación: función que presenta la pantalla de presentación %*************************************************************** clear,clc,cla,close all %Creamos figura figdiag=figure('Units','Pixels',... 'Position',[0.06 0.06 0.9 0.9],... %Tamaño de la presentación 'Number','off',... 'Name','UNSCH/FIQM/EFP Ingenieria en Industrias Alimentarias' , ... 'Menubar','none', ... 'color',[0 0 0]); %Ubicamos ejes en figura axes('Units','Normalized',... 'Position',[0 0 1 1]); %-----Centramos la figura--------scrsz = get(0, 'ScreenSize'); pos_act=get(gcf,'Position'); xr=scrsz(3) - pos_act(3); xp=round(xr/2); yr=scrsz(4) - pos_act(4); yp=round(yr/2); set(gcf,'Position',[xp yp pos_act(3) pos_act(4)]); %--------------------------------------%Incluir imagen %Importamos imagen *.jpg,junto con su mapa de colores [x,map]=imread('caratula.jpg','jpg'); %Representamos imagen en figura, con su mapa de colores image(x),colormap(map),axis off ,hold on %Títulos sobre imagen %Título text(150,90,'E.F.P.Ingenieria en Industrias Alimentarias','Fontname','Arial','Fontsize',30,'Fontangle','Italic', ... 'Fontweight','Bold','color',[1 0 0]); text(780,1070,'Ing. Alberto HUAMANI HUAMANI','Fontname', ... 'Comic Sans MS','Fontangle','Italic','Fontweight','Bold', ... 'Fontsize',18,'color',[0 0 1]); %Botón Continuar -------------------------------------------------------------------------------------------------------------------------------------------------------4 Alberto L.Huamaní Huamaní


Solución numérica

botok=uicontrol('Style','pushbutton', ... 'Units','normalized', ... 'Position',[.84 .03 .12 .05], ... 'String','CONTINUAR',... 'Callback','clear all; close all;clc; METODOS_NUMERICOS;'); %GUI es el nombre del siguiente programa. 6. 7.

Guarda el programa en la carpeta creada como: caratula.m Ejecutando deberá salir el siguiente esquema( el fondo es la imagen guardada previamente como jpg)

II. Programa para crear el siguiente formulario

1.

Abrir la ventana de interfaz GUIDE, Haciendo click en la figurita o escribiendo guide en la entana de window

-------------------------------------------------------------------------------------------------------------------------------------------------------5 Alberto L.Huamaní Huamaní


Solución numérica

2.

Aparecerá la siguiente ventana y elegir Create New Guide y hacer click en OK

3.

Aparecerá la siguiente ventana



Solución numérica

4. Confeccionamos el siguiente formulario Hacer click en una de los comandos (Push Button) y arrastrar a la parte del formulario

5.

Hacer doble click en Push Button y escribir en la ventana de propiedades de String: SOLUCION NUMERICA

Y asi sucesivamente para todos los comandos mostrados en el siguiente formulario Que tiene todo los comandos de Push Button



6.

Solución numérica

Una vez llenado el formulario guardar el formulario dentro de la carpeta creada : Guardar como: METODOS_NUMERICOS.fig Automáticamente se creara el archivo.m

III. Programa para el siguiente formulario



Metodo Bisección

Práctica 1-1

METODO DE BISECCION 1.1 MÉTODO DE BISECCIÓN El método de bisección, conocido también como de corte binario, de partición de intervalos o de Bolzano, es un tipo de búsqueda incremental en el que el intervalo se divide siempre a la mitad. Si la función cambia de signo sobre un intervalo, se evalúa el valor de la función en el punto medio. La posición de la raíz se determina situándola en el punto medio del subintervalo, dentro del cual ocurre un cambio de signo. El proceso se repite hasta obtener una mejor aproximación (Chapra, 2007). El método de bisección se puede aplicar para resolver ecuaciones no lineales como f (x) = 0. El método de de Bisección tiene como base el Teorema de valor Intermedio, el cual a la letra dice: Teorema de valor intermedio Sea f(x) una función continua en [x 1, x 2] y sea xr un valor entre (f(x 1), f(x2)), entonces existe un valor x* entre (x1, x2), tal que f(x*) = xr Corolario Sea f(x) una función continua en [x 1, x2] y sean f(x1) y f(x2) de signos contrarios, entonces existe un valor x* entre (x 1, x2), tal que f(x*) = 0 Paso 1: Elija valores iniciales inferior, xl , y superior, x 2, que encierren la raíz, de forma tal que la función cambie de signo en el intervalo. Esto se verifica comprobando que f  x1  f  x2   0 . Paso 2: Una aproximación de la raíz xr se determina mediante: xr 

x1  x2

(1.8) 2 Paso 3: Realice las siguientes evaluaciones para determinar en qué subintervalo está la raíz: a) b)

Si f  x1  f  xr   0 , entonces la raíz se encuentra dentro del subintervalo inferior o izquierdo. Por lo tanto, haga x2 = xr y vuelva al paso 2. Si f  x1  f  xr   0 , entonces la raíz se encuentra dentro del subintervalo superior o derecho. Por lo tanto, haga x1 = x r y vuelva al paso 2.



c)

Metodo Bisección

Si f  x1  f  xr   0 , la raíz es igual a x r; termina el cálculo.

Concluyendo. El método de bisección tiene limitaciones sobre otros métodos numéricos, para obtener raíces de ecuaciones no lineales, sin embargo da resultados aproximados para la simplicidad del algoritmo. 1.

Abrir nuevamente un formulario en blanco y digitar los Push Button

Una vez digitado los comandos, guardar como: SOLUCION_NUMERICA.fig Y aparecerá nuevamente el archivo.m para este formulario. function NUMERICA_Callback(hObject, eventdata, handles) % hObject handle to NUMERICA (see GCBO) SOLUCION_NUMERICA 2. Abrir nuevamente un formulario en blanco y digitar Como en el paso 12, y confeccionamos el formulario como vemos en la figura Formulario

Procedimiento de llenado del formulario Paso 1: Realizamos el formulario copiando Static text , Edit text, table para la salida de resultados, axes1 para el grafico y push button Static text1: para texto: hacer doble click en la tabla de código y llenar: String: Funcion matematica Style: texto Edit text1: para ingreso de valores: hacer doble click en la tabla de código y llenar: String: Style:edit Tag: edit1 -------------------------------------------------------------------------------------------------------------------------------------------------------10 Alberto L.Huamaní Huamaní


Metodo Bisección

Edit text: para salida de valores: hacer doble click en la tabla de código y llenar: String: nombre a poner Style:edit Tag: edit2 o nombre de la variable Table: insertar un table y hacer anticlick y en table property poner el nombre de las variables de salida que se desea mostrar. En propiedades poner en Tag: tabla Paso 2: Completado el formulario, guardar el formulario como, guardar archivo como: METODO_BISECCION.fig 1) Inmediatamente se autogenera el archivo.m, 2) En el archivo generado buscar Pushbutton1. Digitar el código del programa en el Push button1, después de la línea siguiente para el caso y así sucesivamente. Paso 2: Digitar el Programa function pushbutton1 _Callback(hObject, eventdata, handles) f=get(handles.edit1,'string'); f=inline(f); xai=str2double(get(handles.edit2,'string')); % valor de x1 xbi=str2double(get(handles.edit3,'string')); % valor de x2 tol=str2double(get(handles.edit4,'string')); % error i=1; ea(1)=100; %%%%% Metodo Bisección %%%%%% if f(xai)*f(xbi)<0; % Comprobando que la raiz se encuentra en este intervalo xa(1)=xai; xb(1)=xbi; xr(1)=(xa(1)+xb(1))/2; %Limpiar tabla antes de mostrar resultado set(handles.tabla,'Data',{}) % Limpiar tabla, grafico en caso de que antes se haya graficado una funcion hold off cla set(handles.tabla,'Data',{}) set(handles.respuesta,'string','No hay raiz'); while abs(ea(i))>=tol if f(xa(i))*f(xr(i))<0 % Condicion de cumplimiento xa(i+1)=xa(i); xb(i+1)=xr(i); % Es la raiz(xr) si se cumple condicion end if f(xa(i))*f(xr(i))>0 % Condicion de cumplimiento xa(i+1)=xr(i); % Es la raiz(xr) si se cumple condicion xb(i+1)=xb(i); end xr(i+1)=(xa(i+1)+xb(i+1))/2; % Valor intermedio para 2° iteracion -------------------------------------------------------------------------------------------------------------------------------------------------------11 Alberto L.Huamaní Huamaní


Metodo Bisección

ea(i+1)=abs((xr(i+1)-xr(i))/(xr(i+1))*100);% error absoluto % Mostrara datos en tabla valores = {i xa(i+1) xb(i+1) xr(i+1) ea(i+1)}; temp=get(handles.tabla,'data'); valoresNuevos=[valores;temp]; set(handles.tabla,'Data',valoresNuevos) i=i+1; end % Cerramos while % Mostrando respuesta en textbox con formato coma flotante a 6 cifras decimales respuesta=sprintf('%0.6f',xr(i)); set(handles.respuesta,'string',respuesta); %Grafica de la funcion fplot(handles.axes1,f,[xai xbi]); grid on; hold on; handles.axes1=plot(xr(i),subs(f,respuesta),'r*'); else set(handles.respuesta,'string','No existe la raiz en el intervalo'); zoom on end function pushbutton2 _Callback(hObject, eventdata, handles) cla %limpiar tabla set(handles.tabla,'Data',{}) %limpiar textboxs set(handles.edit1,'string',''); set(handles.edit2,'string',''); set(handles.edit3,'string',''); set(handles.edit4,'string',''); set(handles.respuesta,'string',''); function pushbutton3 _Callback(hObject, eventdata, handles) SOLUCION_NUMERICA Paso 3: Correr el Programa, llenando los espacios en blanco con la función matemática, valor de x1, x2 y el error y hacer click en calcular Planteamiento del problema. Utilice el método gráfico para determinar el coeficiente de arrastre c necesario para que un paracaidista de masa m = 68.1 kg tenga una velocidad de 40 m/s después de una caída libre de t = 10 s. Nota: La aceleración de la gravedad es 9.8 m/s2. Solución. Este problema se resuelve determinando la raíz de la ecuación usando los parámetros t = 10, g = 9.8, v = 40 y m = 68.1: La ecuación resulta:



Metodo Bisección

 c     *t  mg  * 1  e  m   v  c 



f ( x )  f ( x) 



mg   c   * 1  exp   * t    v c   m   667,38 x

* 1  exp 0,146843* x   40

f(x)= (667.38/x)*(1-exp(-0.146843*x))-40 Resultados



Metodo Bisección

1.1.2 Usando archivo.m Paso 1: Abrir la carpeta New Script Paso 2: Digitación del código para el método bisección (caso con ingreso de función matemática) % Cálculo de ecuación matemática por el método de la bisección % Alberto HUAMANI HUAMANI % 2016 disp(' METODO DE LA BISECCION ') disp(' Metodos matemáticos en industrias alimentarias' ) disp(' ') f = input(' INGRESE LA FUNCION en x :','s'); xai = input(' INGRESE LIMITE INFERIOR DEL INTERVALO: '); xbi = input(' INGRESE LIMITE SUPERIOR DEL INTERVALO: '); tol= input(' INGRESE EL PORCENTAJE DE ERROR: '); f=inline(f); i=1; ea(1)=100; if f(xai)*f(xbi)<0; % Comprobando que la raiz se encuentra en este intervalo xa(1)=xai; xb(1)=xbi; xr(1)=(xa(1)+xb(1))/2; fprintf(' it xa xr xb error aprox\n'); fprintf('%2d\t %11.7f \t %11.7f \t %11.7f\n', i,xa(i),xr(i),xb(i)); while abs(ea(i))>=tol, if f(xa(i))*f(xr(i))<0 % Condicion de cumplimiento xa(i+1)=xa(i); xb(i+1)=xr(i); % Es la raiz(xr) si se cumple condicion end if f(xa(i))*f(xr(i))>0 % Condicion de cumplimiento xa(i+1)=xr(i); % Es la raiz(xr) si se cumple condicion xb(i+1)=xb(i); end xr(i+1)=(xa(i+1)+xb(i+1))/2; % Valor intermedio para 2° iteracion ea(i+1)=abs((xr(i+1)-xr(i))/(xr(i+1))*100);% error absolute fprintf('%2d\t %11.7f \t %11.7f \t %11.7f \t, %7.3f\n', i+1,xa(i+1),xr(i+1),xb(i+1), ea(i+1)); i=i+1; end % Cerramos while else set(handles.respuesta,'string','No existe la raiz en el intervalo'); end Paso 3: Guardar el archive como biseccion.m Paso 4: Ejecutar haciendo click en la flecha verde Paso 5: Ingresar la función, valor de x1, x2 y error Luego de ejecutar se tiene el resultado -------------------------------------------------------------------------------------------------------------------------------------------------------14 Alberto L.Huamaní Huamaní


Metodo Bisección

Resultado >> FALSA_POSICIO_2016 METODO DE LA BISECCION Metodos matemáticos en industrias alimentarias INGRESE LA FUNCION en x :(667.38/x)*(1-exp(-0.146843*x))-40 INGRESE LIMITE INFERIOR DEL INTERVALO: 12 INGRESE LIMITE SUPERIOR DEL INTERVALO: 16 INGRESE EL PORCENTAJE DE ERROR: 0.001 Resultado it 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 >>

xa 12.0000000 14.0000000 14.0000000 14.5000000 14.7500000 14.7500000 14.7500000 14.7500000 14.7656250 14.7734375 14.7773438 14.7792969 14.7792969 14.7797852 14.7800293

xr 14.0000000 15.0000000 14.5000000 14.7500000 14.8750000 14.8125000 14.7812500 14.7656250 14.7734375 14.7773438 14.7792969 14.7802734 14.7797852 14.7800293 14.7801514

xb error aprox 16.0000000 16.0000000 , 6.667 15.0000000 , 3.448 15.0000000 , 1.695 15.0000000 , 0.840 14.8750000 , 0.422 14.8125000 , 0.211 14.7812500 , 0.106 14.7812500 , 0.053 14.7812500 , 0.026 14.7812500 , 0.013 14.7812500 , 0.007 14.7802734 , 0.003 14.7802734 , 0.002 14.7802734 , 0.001



Metodo Bisección



Metodo Falsa Posición

Práctica 1-2

METODO FALSA POSICI N 1.2 MÉTODO DE LA FALSA POSICIÓN O REGULA FALSI Un inconveniente del método de bisección es que al dividir el intervalo de x l a x2 en mitades iguales, no se toman en consideración las magnitudes de f(x l) y f(x2). Por ejemplo, si f(xl) está mucho más cercana a cero que f(x 2), es lógico que la raíz se encuentre más cerca de x l que de x 2. Un método alternativo que aprovecha esta visualización gráfica consiste en unir f(x l) y f(x2) con una línea recta. La intersección de esta línea con el eje de las x representa una mejor aproximación de la raíz.También se le conoce como método de interpolacion lineal (Chapra, 2007). Usando triángulos semejantes, la intersección de la línea recta con el eje de las x se estima mediante f  x1 

 x

f  x1 



f  x2 

 x

f  x2 

(1.9)

Despejando x r x f  x2 

f  x2  x1  x2  f ( x1 )  f ( x2 )

1.2.1 Solución usando Programa GUIDE Matlab Formulario


(1.10)



Programa function varargout = pushbutton1 _Callback(h, eventdata, handles, varargin) f=get(handles.edit1,'string'); f=inline(f); xai=str2double(get(handles.edit2,'string')); % valor de x1 xbi=str2double(get(handles.edit3,'string')); % valor de x2 tol=str2double(get(handles.edit4,'string')); % error i=1; ea(1)=100; %%%%% Metodo Bisección %%%%%% if f(xai)*f(xbi)<0; % Comprobando que la raiz se encuentra en este intervalo xa(1)=xai; xb(1)=xbi; xf(1)= xb(1)-f(xb(1))*(xa(1) - xb(1))/(f(xa(1))-f(xb(1))); %Limpiar tabla antes de mostrar resultado set(handles.tabla,'Data',{}) % Limpiar tabla, grafico en caso de que antes se haya graficado una funcion hold off cla set(handles.tabla,'Data',{}) set(handles.respuesta,'string','No hay raiz'); while abs(ea(i))>=tol if f(xa(i))*f(xf(i))<0 % Condicion de cumplimiento xa(i+1)=xa(i); xb(i+1)=xf(i); % Es la raiz(xr) si se cumple condicion end if f(xa(i))*f(xf(i))>0 % Condicion de cumplimiento xa(i+1)=xf(i); % Es la raiz(xr) si se cumple condicion xb(i+1)=xb(i); end xf(i+1)= xb(i+1)-f(xb(i+1))*(xa(i+1) - xb(i+1))/(f(xa(i+1))-f(xb(i+1))); ea(i+1)=abs((xf(i+1)-xf(i))/(xf(i+1))*100);% error absoluto % Mostrara datos en tabla valores = {i xa(i+1) xb(i+1) xf(i+1) ea(i+1)}; temp=get(handles.tabla,'data'); valoresNuevos=[valores;temp]; set(handles.tabla,'Data',valoresNuevos) i=i+1; end % Cerramos while % Mostrando respuesta en textbox con formato coma flotante a 6 cifras decimales respuesta=sprintf('%0.6f',xf(i)); set(handles.respuesta,'string',respuesta); %Grafica de la funcion fplot(handles.axes1,f,[xai xbi]); grid on; hold on; handles.axes1=plot(xf(i),subs(f,respuesta),'r*'); else set(handles.respuesta,'string','No existe la raiz en el intervalo'); -------------------------------------------------------------------------------------------------------------------------------------------------------18 Alberto L.Huamaní Huamaní



zoom on end function varargout = pushbutton2 _Callback(h, eventdata, handles, varargin) cla %limpiar tabla set(handles.tabla,'Data',{}) %limpiar textboxs set(handles.edit1,'string',''); set(handles.edit2,'string',''); set(handles.edit3,'string',''); set(handles.edit4,'string',''); set(handles.respuesta,'string',''); function varargout = pushbutton3 _Callback(h, eventdata, handles, varargin) SOLUCION_NUMERICA Ejecutar: hacer click en la flecha verde Paso 5: Ingresamos la función matemática o los valores numéricos según sea el caso, ingresar los limites inferior y superior y el error; y se tiene el resultado Ingrese la función:

f ( x) 

667,38 x

* 1  exp 0,146843* x   40

f(x)= (667.38/x)*(1-exp(-0.146843*x))-40 xa1= 12 xb1=16

tol: 0.05; Hacer Click en calcular

Desventajas del método de la falsa posición Aunque el método de la falsa posición parecería ser siempre la mejor opción entre los métodos cerrados, hay casos donde funciona de manera deficiente. En efecto, como en el ejemplo siguiente, hay ciertos casos donde el método de bisección ofrece mejores resultados. -------------------------------------------------------------------------------------------------------------------------------------------------------19 Alberto L.Huamaní Huamaní





Metodo Newton Rapson

Práctica 1-3

NEWTON-RAPHSON 1.3 MÉTODO DE NEWTON-RAPHSON O DE LA TANGENTE Si el valor inicial para la raíz es xi, entonces se puede trazar una tangente desde el punto [xi, f(xi)] de la curva. Por lo común, el punto donde esta tangente cruza al eje x representa una aproximación mejorada de la raíz.

Figura 1.3. Método de Newton-Raphson El método de Newton-Raphson se deduce a partir de esta interpretación geométrica (un método alternativo basado en la serie de Taylor. De la figura 1.3, se tiene que la primera derivada en x es equivalente a la pendiente:

f '  xi  

f  xi   0 xi  xi 1

(1.11)

Reordenando la ecuación anterior xi 1  xi 

f ( xi ) f '  xi 

Deducción la fórmula de Newton-Raphson usando una serie de Taylor. La expansión de la serie de Taylor se puede expresar como: -------------------------------------------------------------------------------------------------------------------------------------------------------21 Alberto L.Huamaní Huamaní

(1.12)



2

f ( xi 1 )  f  xi   f '  xi  xi 1  xi  

f ' ' ( xi ) xi 1  xi 

2!

 ... (1.13)

Truncando la serie de Taylor después del término de la primera derivada, se obtiene una versión aproximada: f ( xi 1 )  f  xi   f '  xi  xi1  xi 

(1.14)

En la intersección con el eje x, f(x i+1) debe ser igual a cero, o 0  f  xi   f '  xi  xi1  xi  xi 1  xi 

f ( xi ) f '  xi 

(1.15)

(1.16)

Desventajas del método de Newton-Raphson Aunque en general el método de Newton-Raphson es muy eficiente, hay situaciones donde se comporta de manera deficiente. Por ejemplo en el caso especial de raíces múltiples. Sin embargo, también cuando se trata de raíces simples, se encuentran dificultades. 1.4 MÉTODO NEWTON MEJORADO Una de las condiciones para garantizar la convergencia del método de Newton es que f´(x) tiene que ser diferente de cero . Si al ejecutar el método de Newton se observa que f´(xn) se aproxima a cero, la rapidez del método disminuye y hay una posible raíz múltiple. El método de raíz múltiple también es conocido como el método de Newton mejorado, y básicamente su estructura es muy similar excepto de que se debe hallar la segunda derivada. Si en lugar de considerar los dos primeros términos de la serie de Taylor se consideran los tres primeros términos , se representa con Δxi a la diferencia entre x i+1 y xi y se iguala a cero, se tiene: f  xi    xi f '  xi  

y sustituyendo Δxi por 

 xi 2 2

f  xi  f '  xi 

f ' '  xi   0

(1.17)

(a partir de la fórmula de Newton-Raphson) queda:





1 f  xi 



2 f '  xi 

f  xi    xi  f '  xi  




f ' '  xi   0



(1.18)

Despejando ΔXi se obtiene:  xi   f '  xi  

f  xi  f  xi 

2 f '  xi 

(1.19) f ' '  xi 

De la ecuación despejando el valor de x xi 1  xi 

f '  xi  

f  xi  f  xi 

2 f '  xi 

i+1:

(1.20) f ' '  xi 

1.4.1 Procedimiento de programa en GUIDE de Matlab Formulario

Para Pop Up Menu, En string hacer doble click y escribir después de Seleccionar:

Programa function pushbutton1 _Callback(hObject, eventdata, handles) f=get(handles.edit1,'string'); x0=str2double(get(handles.edit2,'string')); -------------------------------------------------------------------------------------------------------------------------------------------------------23 Alberto L.Huamaní Huamaní



tol=str2double(get(handles.edit3,'string')); i=1; fx(i)=x0; % valor inicial de fx(i) syms x; % syms declarar la variable x ea(1)=100; f1=subs(f,x,fx(i)); % Evaluacion numerica de f en funcion de x para fx z=diff(f); % derivada de f d=subs(z,x,fx(i)); % evaluacion de z z2=diff(f,2); % calculo de la segunda derivada de f d2=subs(z,x,fx(i)); % evaluacion de z % Opciones de calculo v=get(handles.seleccionar,'value'); switch v case 2 %Limpiar tabla antes de mostrar resultado set(handles.tabla,'Data',{}) while abs(ea(i))>=tol; fx(i+1)=fx(i)-f1/d; % Expresion de Newton f1=subs(f,x,fx(i+1)); % Evalua f1 d=subs(z,x,fx(i+1)); % Evalua z ea(i+1)=abs((fx(i+1)-fx(i))/fx(i+1)*100);% Error absoluto i=i+1; end for j=1:i; %mostrara datos en tabla valores = {j-1,fx(j),ea(j)}; temp=get(handles.tabla,'data'); valoresNuevos=[valores;temp]; set(handles.tabla,'Data',valoresNuevos) end % Mostrando de raiz en textbox con formato coma flotante a 6 cifras decimales raiz=sprintf('%0.6f',fx(j)); set(handles.raiz,'string',raiz); %end %Grafica de la funcion hold off fplot(handles.axes1,f,[0 fx(j)+1]); grid on; hold on; handles.axes1=plot(fx(j),subs(f,raiz),'r*'); zoom on case 3 %Limpiar tabla antes de mostrar resultado set(handles.tabla,'Data',{}) while abs(ea(i))>=tol; fx(i+1)=fx(i)-(f1*d)/(d^2-(f1*d2)); % Expresion de Newton f1=subs(f,x,fx(i+1)); % Evaluacion numerica de f en funcion de x para fx d=subs(z,x,fx(i+1)); % Evalua z d2=subs(z,x,fx(i+1)); % evaluacion de z -------------------------------------------------------------------------------------------------------------------------------------------------------24 Alberto L.Huamaní Huamaní



ea(i+1)=abs((fx(i+1)-fx(i))/fx(i+1)*100);% Error absoluto i=i+1; end for j=1:i; % Mostrara datos en tabla valores = {j-1,fx(j),ea(j)}; temp=get(handles.tabla,'data'); valoresNuevos=[valores;temp]; set(handles.tabla,'Data',valoresNuevos) end %Mostrando de raiz en textbox con formato coma flotante a 6 cifras decimales raiz=sprintf('%0.6f',fx(j)); set(handles.raiz,'string',raiz); %Grafica de la funcion hold off fplot(handles.axes1,f,[0 fx(j)+1]); grid on; hold on; handles.axes1=plot(fx(j),subs(f,raiz),'r*'); zoom on end function pushbutton2 _Callback(hObject, eventdata, handles) % hObject handle to pushbutton2 (see GCBO) cla %limpiar tabla set(handles.tabla,'Data',{}) %limpiar textboxs set(handles.edit1,'string',''); set(handles.edit2,'string',''); set(handles.edit3,'string',''); set(handles.raiz,'string',''); set(handles.advertir,'string',''); Compilación Ingrese la función:

f ( x) 

667,38 x

* 1  exp 0,146843* x   40

f(x)= (667.38/x)*(1-exp(-0.146843*x))-40 xi= 12; tol: 0.05 a)

Newton Rapson



b)


Newton Rapson mejorado



Metodo Secante

Práctica 1-4

METODO SECANTE 1.5 MÉTODO DE LA SECANTE Surge como una variación del método de Newton-Raphson, en lugar de tomar la tangente se toma la secante. De manera que la derivada se aproxima por una diferencia finita dividida hacia atras, basada en las estimaciones sucesivas es decir, como en (figura 1.8)

Figura 1.8. Método de la Secante Esta técnica es similar a la del método de Newton-Raphson (figura 1.8) en el sentido de que una aproximación de la raíz se predice extrapolando una tangente de la función hasta el eje x. Sin embargo, el método de la secante usa una diferencia dividida en lugar de una derivada para estimar la pendiente (Chapra, 2007). f '  xi  

f  xi 1   f  xi  xi 1  xi

(1.21)

Esto puede sustituirse en la fórmula (1), quedando asi la formula de la secante: xi 1  xi 

f  xi  xi 1  xi  f  xi 1   f  xi 

(1.22)

El método requiere de dos valores iniciales pero como no se requiere que f(x) cambie de signo en el intervalo considerado, no se lo incluye dentro de los métodos que utilizan intervalos, este método no se clasifica como un método cerrado. -------------------------------------------------------------------------------------------------------------------------------------------------------27 Alberto L.Huamaní Huamaní


Metodo Secante

Observe la similitud entre los métodos de la secante y de la falsa posición. Por ejemplo, las ecuaciones son idénticas en todos los términos. Ambas usan dos valores iniciales para calcular una aproximación de la pendiente de la función que se utiliza para proyectar hacia el eje x una nueva aproximación de la raíz. Sin embargo, existe una diferencia crítica entre ambos métodos. Tal diferencia estriba en la forma en que uno de los valores iniciales se reemplaza por la nueva aproximación. 1.5.1 Procedimiento de Programa en GUIDE de Matlab Formulario

Programa function pushbutton1 _Callback(hObject, eventdata, handles) f=get(handles.edit1,'string'); x0=str2double(get(handles.edit2,'string'));% x1=str2double(get(handles.edit3,'string')); tol=str2double(get(handles.edit4,'string')); syms x; ea(1)=100; %Limpiar tabla antes de mostrar resultado set(handles.tabla,'Data',{}) %Limpiar tabla, grafico en caso de que antes se haya graficado una funcion hold off cla set(handles.tabla,'Data',{}) set(handles.respuesta,'string','No hay raiz'); i=1; -------------------------------------------------------------------------------------------------------------------------------------------------------28 Alberto L.Huamaní Huamaní


Metodo Secante

while abs(ea)>tol; x=x0; g=eval(f); x=x1; gg=eval(f); xi=x1-((gg*(x0-x1))/(g-gg)); ea=abs((xi-x1)/xi)*100; x0=x1; x1=xi; % Mostrara datos en tabla valores = {i,x xi,ea}; temp=get(handles.tabla,'data'); valoresNuevos=[valores;temp]; set(handles.tabla,'Data',valoresNuevos) i=i+1; end %Mostrando respuesta en textbox con formato coma flotante a 6 cifras decimales respuesta=sprintf('%0.6f',xi); set(handles.respuesta,'string',respuesta); %Grafica de la funcion hold off fplot(handles.axes1,f,[0 xi+1]); grid on; hold on; handles.axes1=plot(xi,subs(f,respuesta),'r*'); zoom on function pushbutton2 _Callback(hObject, eventdata, handles) cla %limpiar tabla set(handles.tabla,'Data',{}) %limpiar textboxs set(handles.edit1,'string',''); set(handles.edit2,'string',''); set(handles.edit3,'string',''); set(handles.edit4,'string',''); set(handles.respuesta,'string',''); function pushbutton3 _Callback(hObject, eventdata, handles) SOLUCION_NUMERICA Resultado



Metodo Secante



Aplicaciones a Ingeniería de Alimentos

Práctica 1-5

APLICACIONES EN INGENIERIA DE ALIMENTOS 1. Ejercicio 1: Aplicación coeficiente de friccion La siguiente relación entre el factor de fricción f y el número de Reynolds (Re) se cumple cuando hay flujo turbulento de un fluido en un tubo liso. 1 f

 0.4  1.74 Ln(Re f )

Construya una tabla de valores de f correspondientes a números de Reynolds de 10 4 hasta 106. Solución La función transformada es: 1  0.4  1.74 Ln(Re f ) f

1  0.4 f  1.74 f * Ln(Re

f )

Entonces la función a solucionar es: g ( f )  1  0.4 f  1.74 f * Ln(Re

f )

fx = -1 - (0.4 * x) + 1.74 * x * log(Re * (x ^ 0.5)) Formulario




Programa Para incorporar imagen function Problema_1_OpeningFcn(hObject, eventdata, handles, varargin) % This function has no output args, see OutputFcn. handles.output = hObject; % Update handles structure guidata(hObject, handles); axes(handles.axes1) background=imread('friccion.jpg'); axis off imshow(background) function pushbutton2 _Callback(hObject, eventdata, handles) x0=str2double(get(handles.edit1,'string'));% x1=str2double(get(handles.edit2,'string')); tol=str2double(get(handles.edit3,'string')); set(handles.tabla,'Data',{}) %Limpiar tabla, grafico en caso de que antes se haya graficado una funcion hold off %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Re = [10000; 50000;100000; 200000; 300000; 400000; 500000; 700000;900000;10^6]; n=length(Re); for k=1:n i=0;% iteracion del metodo numerico fxs=1; while abs(fxs)>tol &(i<10); fx0 = -1-(0.4*x0) + 1.74*x0*log(Re(k)*(x0^0.5)); fx1 = -1-(0.4*x1) + 1.74*x1*log(Re(k)*(x1^0.5)); xs=x1-((fx1*(x0-x1))/(fx0-fx1)); % algoritmo de metodo secante fxs = -1-(0.4*xs) + 1.74*xs*log(Re(k)*(xs^0.5)); ea=abs((xs-x1)/xs)*100; x0=x1; x1=xs; i=i+1; end % CALCULOS X(k) = xs; f=X(k); % Mostrara datos en tabla valores = {Re(k) f}; temp=get(handles.tabla,'data'); valoresNuevos=[valores;temp]; set(handles.tabla,'Data',valoresNuevos) k=k+1; %Grafica de la funcion hold on; axes(handles.axes2) -------------------------------------------------------------------------------------------------------------------------------------------------------32 Alberto L.Huamaní Huamaní



plot(Re(k),f,'b*') legend ( 'Factor de friccion en funcion del Re' ) ylabel( 'Valor de f' ) xlabel( 'Re' ) grid on; end function pushbutton2 _Callback(hObject, eventdata, handles) cla %limpiar tabla set(handles.tabla,'Data',{}) %limpiar textboxs set(handles.edit1,'string',''); set(handles.edit2,'string',''); set(handles.edit3,'string',''); set(handles.edit4,'string',''); set(handles.edit5,'string',''); function pushbutton2 _Callback(hObject, eventdata, handles) close Solucionario




2. Ejercicio 2: Aplicación en destilación Equilibrio liquido vapor y la determinación de propiedades características de este estado como son la temperatura y las composiciones. Introducción Varios procesos industriales importantes, por ejemplo, destilación, absorción y extracción, ponen en contacto a dos fases entre las que, cuando no están en equilibrio, se efectúa una transferencia de masa. La velocidad de transferencia de cada especie depende de la separación del sistema respecto al equilibrio (T, P.X, Y) del sistema. En la mayor parte de los procesos industriales las fases que coexisten son vapor y liquida, aunque también se han encontrado sistemas liquido-liquido, vapor – sólido y y liquido sólido. A continuación, haremos un planteamiento de un problema en donde se requiere conocer el comportamiento en el equilibrio para un sistema líquido vapor y los cálculos correspondientes para determinar la temperatura y composiciones de las fases de este sistema. Problema: Considere un líquido en equilibrio con su vapor. Si el líquido está formado por los componentes 1,2,3,4; con los datos dados a continuación calcule la temperatura y la composición del vapor en el equilibrio a la presión total de 75 psia . Componente

Composición del líquido (%mol)

1 2 3 4

10 54 30 6

Presión del vapor del componente puro (psia) 150 K 200 K

25 14,7 4 0,5

200 60 14,7 5

Para resolver este problema se plantean las siguientes ecuaciones: Para la presión de vapor:

   A  BT

Ln pi

0

i

i

(1)

1

Donde i =1, 2, 3, 4 y T (K). La presión total del sistema será: PT 

 Pi

(2)

Considerando que la mezcla de estos cuatro componentes, a las condiciones de presión y temperatura dadas, obedecen las leyes de Raoult y de Dalton. PT 

 p

0 i

* xi


(3)



Donde: pi0 = Presión de vapor de cada componente. PT = presión total del sistema. pi = Presión parcial de cada componente. xi = Fracción mol de cada componente en el líquido. De la ecuación de presión de vapor se tiene que 

pi  exp Ai  0



Bi 



i  1,2,3,4

T 1 

Despejando piº de 1 y reemplazándola en 3 tenemos: PT 

 



xi * exp  Ai 

Bi 



T 

Entonces despejando nos queda una ecuación la cual es función de la temperatura. La ecuación es la siguiente: f T   PT 



 

xi * exp  Ai 

Bi 

0

T 

(4)

Para obtener A i y Bi realizamos el siguiente procedimiento: Hacemos p1º, i = presión de vapor del componente i a T 1 =150 K p2º, i = presión de vapor del componente i a T 2 = 200K Entonces



0



0



Bi



Bi

Ln pi , i  Ai  Ln pi , i  Ai 

T 1 T 2

i  1,2,3,4

(5)

i  1,2,3,4

(6)

Restando estas ecuaciones se tiene  p10 , i   1 1  Ln 0   Bi     p , i   T 1 T 2   2 

De donde  p10 , i  Ln 0   p , i   2  Bi   1 1      T 1 T 2  -------------------------------------------------------------------------------------------------------------------------------------------------------35 Alberto L.Huamaní Huamaní



Reemplazando estos valores conocemos Bi y podemos obtener Ai de la ecuación (4). Valores iniciales Ahora para hallar un valor inicial de T para resolver la ecuación 4, se considera el componente dominante de la mezcla que en este caso de acuerdo a los datos dados en la tabla es el componente 2, y se usa PT en lugar de p 2º en la ecuación 1 que es la de presión de vapor. Es decir, LnPT   A2 

B2 T

De donde T 

B2 LnPT 

 A2

Con este resultado inicial y las consideraciones ya mencionadas, utilizamos el método de Newton - Raphson para hallar la temperatura del sistema (temperatura de burbuja) en el equilibrio. Método de Newton – Raphson

xi 1  xi 

f  xi  f '  xi 

Dónde: f’ (T) = - Σ xi exp ( Ai + Bi / T ) * ( - Bi / T2 ) derivada de f

Y  f T   PT 

f ' ' T   



Bi 

 x * exp A  T   0



i

 

xi * exp  Ai 

i

Bi   Bi  *   T   T 2 

Algoritmo Utilizado Para encontrar una raíz de la ecuación f (x i+1) = 0, proporcionar la función f (x i) y su derivada df(xi) y los datos: Datos: Valor inicial x 0, criterio de convergencia (ea) o error absoluto, criterio de exactitud (ea1) y número máximo de iteraciones i. Resultados: La raíz aproximada x o un mensaje de falla. PASO 1: Hacer i = 1 PASO 2: Mientras i< MAXIT, repetir los pasos 3 a 7. -------------------------------------------------------------------------------------------------------------------------------------------------------36 Alberto L.Huamaní Huamaní



PASO 3: Hacer x1 = x0 – f(x0) /df (x0) (calcula xi ) PASO 4: Si ABS (x1 – x0) < EPS, entonces IMPRIMIR x1 y TERMINAR. De otro modo CONTINUAR. PASO 5: Si ABS (f(x1)) < EPS1, entonces IMPRIMIR x1 y TERMINAR. De otro modo CONTINUAR. PASO 6: Hacer i = i + 1. PASO 7: Hacer x0 = x1 PASO 8: IMPRIMIR mensaje de falla ‘’ EL MÉTODO NO CONVERGE A UNA RAÍZ ‘’ y terminar. A)

El programa utilizado en matlab en archivo.m es el siguiente:

Función que permite calcular la temperatura de equilibrio. %METODO NEWTON RAPSON clc clear all fprintf('METODOS NUMERICOS APLICADOS A INGENIERIA\n'); fprintf('CALCULO DE TEMPERATURAS DE EQUILIBRIO \n'); fprintf('INGENIERIA EN INDUSTRIAS ALIMENTARIAS \n\n'); P1 = [25; 14.7; 4.0; 0.5]; % Presión de vapor del componente Puro (psia) a 150K P2 = [200.0; 60.0; 14.7; 5.0]; % Presión de vapor del componente Puro (psia) a 200K T1 = 150;% valor 1 de T T2 = 200; % Valor 2 de T B = log(P1./P2)/(1/T1-1/T2); ;% calculo de la constante B A = log(P1)-B/T1 ; ;% Calculo de la constante A X = [0.10; 0.54; 0.30; 0.06]; % Composición del Liquido% mol PT = 75; % composición del vapor en el equilibrio a la presión total de 75 psia i = 0;% iteracion cero f =1; tol = 0.000001; T = B(2)/(log(PT)- A(2)); fprintf (' T f(T) \n', T, f )% impresion en texto de T y f while (abs(f)>tol)&(i<10);% maxima iteracion menor de 10 (I<10) f = PT-sum(X.*exp(A + B/T));% Funcion a resolver df = sum(X.*exp(A+B/T).*(B/T^2));% derivada de la funcion T1 = T-f/df; % Algoritmo de Newton fprintf ('%10.2f %10.2e\n',T,f)% impresion de T y f en valores T = T1; % Valor de la raiz o temperatura i = i+1; end fprintf ('\n\n y(j) \n') for j = 1:i+1 y(j) = (X(j)*exp(A(j) + B(j)/T))/PT; fprintf ('%10.4f \n', y(j)); end




Solucion METODOS NUMERICOS APLICADOS A INGENIERIA CALCULO DE TEMPERATURAS DE EQUILIBRIO INGENIERIA EN INDUSTRIAS ALIMENTARIAS T f(T) 211.17 8.29e-01 211.67 -2.86e-03 211.67 -3.35e-08 y(j) 0.3761 0.5451 0.0729 0.0059 B) Usando una base de datos de P1, P2 y T en Excel 1. Primero creamos una carpeta 2. Segundo Dentro de ello guardamos en archive de excel como DESTILADO, como Libro de Excel 97- 2003 25 14.7 4 0.5

200 60 14.7 5

0.1 0.54 0.3 0.06

La primera coluna P1 segunda columna P2 y tercera columna T 3. Tercero Confeccionar el programa como archive.m clc clear all fprintf('METODOS NUMERICOS APLICADOS A INGENIERIA\n'); fprintf('CALCULO DE TEMPERATURAS DE EQUILIBRIO \n'); fprintf('INGENIERIA EN INDUSTRIAS ALIMENTARIAS \n\n'); Y=xlsread('DESTILADO','Hoja1');% importación de datos de tabla libro1 de Excel P1=Y(:,1);% datos de columna 1 base de datos P2=Y(:,2);% datos de columna 2 base de datos T1 = 150;% valor 1 de T T2 = 200; % Valor 2 de T B = log(P1./P2)/(1/T1-1/T2); % calculo de la constante B A = log(P1)-B/T1 ; % Calculo de la constante A -------------------------------------------------------------------------------------------------------------------------------------------------------38 Alberto L.Huamaní Huamaní



X=Y(:,3);% datos de columna 3 base de datos fracción molar PT = 75; % composición del vapor en el equilibrio a la presión total de 75 psia %METODO NEWTON RAPSON i = 0;% iteracion cero f =1; tol = 0.000001; T = B(2)/(log(PT)- A(2)); fprintf (' T f(T) \n', T, f )% impresion en texto de T y f while (abs(f)>tol)&(i<10);% maxima iteracion menor de 10 (I<10) f = PT-sum(X.*exp(A + B/T));% Funcion a resolver df = sum(X.*exp(A+B/T).*(B/T^2));% derivada de la funcion T1 = T-f/df; % Algoritmo de Newton fprintf ('%10.2f %10.2e\n',T,f)% impresion de T y f en valores T = T1; % Valor de la raiz o temperatura i = i+1; end fprintf ('\n\n y(j) \n') for j = 1:i+1 y(j) = (X(j)*exp(A(j) + B(j)/T))/PT; fprintf ('%10.4f \n', y(j)); end Resultados de copilacion T

f(T) 211.17 8.29e-001 211.67 -2.86e-003 211.67 -3.35e-008 y(j) 0.3761 0.5451 0.0729 0.0059




C) Programa en interfaz de GUIDE Primero hacemos de formulario

Programa function pushbutton1_Callback(hObject, eventdata, handles) % hObject handle to pushbutton1 (see GCBO) Y=xlsread('DESTILADO','Hoja1');% importación de datos de tabla libro1 de excel T1=str2double(get(handles.edit1,'string')); T2=str2double(get(handles.edit2,'string')); PT=str2double(get(handles.edit3,'string')); tol=str2double(get(handles.edit4,'string')); P1=Y(:,1); P2=Y(:,2); X=Y(:,3); B = log(P1./P2)/(1/T1-1/T2); % calculo de la constante B A = log(P1)-B/T1 ; % Calculo de la constante A %Limpiar tabla antes de mostrar resultado set(handles.tabla,'Data',{}) %METODO NEWTON RAPSON i = 0;% iteracion cero f =1; tol = 0.000001; T = B(2)/(log(PT)- A(2)); while (abs(f)>tol)&(i<10);% maxima iteracion menor de 10 (I<10) f = PT-sum(X.*exp(A + B/T));% Funcion a resolver df = sum(X.*exp(A+B/T).*(B/T^2));% derivada de la funcion T1 = T-f/df; % Algoritmo de Newton T = T1; % Valor de la raiz o temperatura i = i+1; -------------------------------------------------------------------------------------------------------------------------------------------------------40 Alberto L.Huamaní Huamaní



end for j = 1:i+1 y(j) = (X(j)*exp(A(j) + B(j)/T))/PT; % Mostrara datos en tabla valores = {T y(j)}; temp=get(handles.tabla,'data'); valoresNuevos=[valores;temp]; set(handles.tabla,'Data',valoresNuevos) end Segundo ingresamos los valores y obtenemos los resultados




3. Ejercicio 3: Aplicación en transferencia de calor: En una de las etapas de elaboración de frutos en almíbar se utiliza una solución caliente de azúcar de 40°Brix. Para realizar el calentamiento se introducen 1000 kg de dicha disolución en un tanque cilíndrico agitado de 1 m de diámetro, perfectamente aislado, provisto de un agitador tipo paleta de 30 cm de diámetro y que gira a 120 rpm. En el tanque se halla sumergido un serpentín helicoidal formado por tubos de acero inoxidable de 12 mm de diámetro interno, 1 mm de espesor de pared y 15 m de longitud total. Por el interior del serpentín circula vapor saturado de agua a 3 atm, que condensa, siendo su coeficiente de convección de calor 9300 W/(m 2 °C). Si la solución se encuentra inicialmente a 16°C; calcular: a) Coeficiente global de transmisión de calor (Ue). b) El tiempo que tarda la solución en alcanzar 60°C. (θ). c) El caudal (m/θ) y cantidad de vapor (w v) necesario para llevar a cabo este calentamiento. d) La velocidad de elevación de temperatura de la disolución cuando se halla a 50°C. Las propiedades del vapor, de Tablas de vapor saturado a 3 atm se obtienen: T = 132.9 °C H w  2721KJ / kg

 w

 2163KJ / Kg

hw  558KJ / kg

Propiedades de la solución de azúcar: Conductividad térmica: 0,814 W/mºC Calor específico = 2,85 KJ/kg ºC Viscosidad:   3,7 x107 exp(2850/ T ) Pa.s T(ºK) Densidad:   1,191  4,8 x104 T g/cc T en ºC Calentamiento de 16°C a 60°C, por tanto las propiedades de la disolución tomarán a una temperatura media t m= 38°C. Solución Datos: Producto: -------------------------------------------------------------------------------------------------------------------------------------------------------42 Alberto L.Huamaní Huamaní



Masa de producto = 100 kg Conductividad térmica: 0,814 W/mºC Calor específico = 2,85 KJ/kg ºC Viscosidad:   3,7 x107 exp(2850/ T ) Densidad:   1,191  4,8 x104 T

Pa.s T(ºK) g/cm3 T en ºC

Tanque: Diámetro = 1m N = 120 rpm Diámetro paleta = 30 cm Serpentín: Diámetro interno = 12 mm Diámetro externo = 14 mm Longitud = 15 m Vapor saturado: Presión = 3 atm a) Calculo de viscosidad y densidad Reemplazando valores se tiene Viscosidad:   3,53x103 Pa.s Densidad:   1173 kg/m3 b) Cálculo del coeficiente he para tanque agitado con calentamiento de vapor sistema serpentín es:  D p 2 N .    0,87   k   

he DT

0 , 62

1/ 3

        k    w   Cp. 

 D p 2 N .  k  he  * 0,87 *     DT  

0, 62

0.14

1/ 3

0.14

 Cp.  *   k 

    *    w 

0,3m2 2s 1 1173kg / m3 

 5,98x104

Calculo de Re, Pr 2

Re 

Pr  

D p N .



Cp. k





3

5,53 x10 Pa.s

2,85KJ / kg º C 3,53 x103 Pa.s  0,814 x10 3 KJ / s.m.º C

 12,4

Reemplazando Re, Pr y otros en la ecuación general



he 

0,814W / mº C

1m 

* 0,875,98 x10



4 0, 62


 3,53 x10 3 Pa.s   12,4   w  

0.14

1/ 3

Obteniéndose la siguiente ecuación: 0,14

he  680,6 w 

W / m 2 C

Como dato tenemos la viscosidad en función de la temperatura  2850   ( 273 tw )   

7  w  3,7 x10 exp

Para obtener he es preciso conocer tw, la temperatura en la pared del serpentín, para calcular w, viscosidad de la disolución a la temperatura de la pared. tm

he tw

Q

Tw

Espesor del tubo

hi T

La velocidad de transmisión de calor, realizando el balance de energía es:

Q  hi Ai T  T w   he Ae t w  t m  Como se supone que la pared no ofrece resistencia a la transmisión de calor: T w  tw Por tanto, de la ecuación anterior queda: t w 

hi d i T  d e T m he 2

hi Ai T  he Ae t w

t w 

hi Ai  he Ae

2

hi d i  d e he 2

2

Reemplazando valores

930012 x10 3 132,9  he 14 x10 3 38 t w  930012 x10 3   he 14 x10 3 



he  680,6 3,7 x107 exp( 2850/(273  tw)) 2

t w 



hi d i T  d e T m 680,6 * 3,7 x10 2





hi d i  d e 680,6 * 3,7 x10 2

2

7

7



0,14

0 ,14

* exp( 2850/(273  tw))

0 ,14

* exp( 2850/(273  tw))






Reemplazamos valores y tenemos 0 ,14

 2850   177 ,98  0 ,0075*5411 ,824* exp (  273 ( tw)   t w  0 ,14  2850   1 ,339  0 ,002*5411 ,824* exp (  273 ( tw)  

El cálculo de tw se realizará por iteraciones y luego se determina he  2850   he  5411,824 *  exp( (273  tw)  

0 ,14

c) Cálculo del coeficiente global : Ue 1 U e



1 he



1 hi ( d i / d e )

d) Del balance energético en el sistema dt

Término de acumulación

A  mCp

Término de entrada

E  U e Ae (T  t )

d 

T = temperatura del vapor condensante t = temperatura de la solución en el tanque  = tiempo mCp

Igualando los dos términos

dt d 

 U e Ae (T  t )

Ecuación diferencial en variables separables, que integrada con la condición límite:  =0 t=to; conduce a la expresión  T  t 0  U e Ae   T  t  mCp

Ln

Expresión que permite calcular, el tiempo de calentamiento para una determinada temperatura o viceversa: Tiempo:

 

m Cp U e Ae

 T  t 0   T t   

Ln

 U A   Temperatura: t  T  (T  t 0 ) exp  e e   m Cp 

Ae   d e L




e) Caudal de vapor y cantidad de condensado m wv  w  Cp (t  t 0 )  Masa de vapor M v  wv 

f) Velocidad de elevación de temperatura dt d 



U e Ae mCp

(T  t )

kJ  3600s   3 2 1931 x10  0,6597m  2 dt  s m º C  h   (132,9  50) d   kJ   1000kg  2,85 kg º C  

dt d 

 133,4º C / h

g) Temperatura a los 50 min     1931 x10 3 kJ  0,6597m 2  3000s      s m2 º C  t  132,9  (132,9  16) exp    kJ     1000kg  2,85 kg º C     t  102,3º C

Como se supone que se trabaja a presión atmosférica, si fuese agua no se tendría la solución acuosa, sino que podría pasar a vapor. Sin embargo, al tratarse de una solución azucarada, es posible que hierva a más de 100ºC, debido al aumento ebulloscopio que producen los sólidos solubles. a) Usando GUIDE Crear una carpeta Dentro de la carpeta debe estar la imagen en jpg Formulario




Programa Para la imagen guidata(hObject, handles); axes(handles.axes1) background=imread('tanque.jpg'); axis off imshow(background) function pushbutton1_Callback(hObject, eventdata, handles) % hObject handle to pushbutton1 (see GCBO) % PROGRAMA DE SERPENTIN m=str2double(get(handles.edit1,'string')); Ti=str2double(get(handles.edit2,'string')); Tf=str2double(get(handles.edit3,'string')); k=str2double(get(handles.edit4,'string')); Cp=str2double(get(handles.edit5,'string')); Dt=str2double(get(handles.edit6,'string')); N=str2double(get(handles.edit7,'string')); Dp=str2double(get(handles.edit8,'string')); di=str2double(get(handles.edit9,'string')); de=str2double(get(handles.edit10,'string')); L=str2double(get(handles.edit11,'string')); P=str2double(get(handles.edit12,'string')); T=str2double(get(handles.edit13,'string')); Hw=str2double(get(handles.edit14,'string')); hw=str2double(get(handles.edit15,'string')); landa=str2double(get(handles.edit16,'string')); hi=str2double(get(handles.edit17,'string')); % Calculos previos tm = (Ti + Tf) / 2; -------------------------------------------------------------------------------------------------------------------------------------------------------47 Alberto L.Huamaní Huamaní



visc = (3.7 * 10 ^ -7) * (exp(2850 / (273.15 + tm))); den = 1191 - ((4.8 * 10 ^ -4) * tm) ; Re = (Dp ^ 2 * (N / 60) * den) / visc; Pr = (Cp * visc) / (k*10^-3); b0 = (0.87 * (Re ^ 0.62) * (Pr ^ 0.333) * k) / Dt; b1=b0*(visc^0.14); a0=hi * di*T; a1 = de* tm; a2 = di* hi; Ae=3.1416*de*L; %METODOS NUMERICOS DE BISECCION x1=50; x2=150; error=0.0001; it=0; fxr=1; while abs(fxr)>error xr=(x1+x2)/2; fx1= -x1+ (a0+a1*( b1*(3.7 * (10 ^ -7) * exp(2850 / (273.15 + x1)))^-0.14))/(a2+de*( b1*(3.7 * (10 ^ -7) * exp(2850 / (273.15 + x1)))^-0.14)); fx2= -x2+ (a0+a1*( b1*(3.7 * (10 ^ -7) * exp(2850 / (273.15 + x2)))^-0.14))/(a2+de*( b1*(3.7 * (10 ^ -7) * exp(2850 / (273.15 + x2)))^-0.14)); fxr= -xr+ (a0+a1*( b1*(3.7 * (10 ^ -7) * exp(2850 / (273.15 + xr)))^-0.14))/(a2+de*( b1*(3.7 * (10 ^ -7) * exp(2850 / (273.15 + xr)))^-0.14)); if (fxr*fx2)>0; x2=xr; else x1=xr; end it=it+1; end %Calculos uw = (3.7 * 10 ^ -7) * (exp(2850 / (273.15 + xr))); he=b1*(uw^-0.14); Ue=hi*he*(di/de)/(hi*(di/de)+he); tiem = (m * Cp / ((Ue*10^-3) * Ae)) * (log((T - Ti) / (T - Tf))); wv = m * Cp * (Tf - Ti) * 3600 / (tiem * landa); mv = (tiem / 60) * wv / 60; %Grafico temperatura en funcion del tiempo syms Z % el tiempo(min) es Z que va de 0 a 25 minutos, la temperatura t del producto se estimacon la funcion t=T-(T-Ti)*exp(-Ue*Ae*60*Z/(m*Cp*1000) axes(handles.axes2) ezplot(T-(T-Ti)*exp(-Ue*Ae*60*Z/(m*Cp*1000)),[0,25]) title('CURVA DE CALENTAMIENTO'); xlabel('Tiempo (min)'),ylabel('Temperatura (ºC)') grid % Salida de resultados -------------------------------------------------------------------------------------------------------------------------------------------------------48 Alberto L.Huamaní Huamaní



set(handles.edit18,'string',visc); set(handles.edit19,'string',den); set(handles.edit20,'string',xr); set(handles.edit21,'string',Re); set(handles.edit22,'string',Pr); set(handles.edit23,'string',he); set(handles.edit24,'string',Ue); set(handles.edit25,'string',tiem); set(handles.edit26,'string',wv); set(handles.edit27,'string',mv); Solución Ingresamos las variables de ingreso y obtenemos como resultado las variables esperadas




4. Ejercicio aplicado 4: Calculo de difusividad del agua en secado de alimentos2 A partir de los datos experimentales de humedad y tiempo de secado usando la ecuación de Alvarez y Legües modificado determinar la difusividad efectiva del proceso, para los 15 términos de la serie.  (2n  1) 2  2   2 D t exp . ef   X 0  X e  n0 ( 2n  1) 2 4 L2   X  X e

x1 



8

1

8  2

 n 2 2  x  x  2 * 2 exp 2 D.t   n  r  8

1

 n 2 2  x  x  x1 * 2 exp 2 D.t  n  r  1

Donde: X = humedad de la muestra (g); Def = difusividade efectiva del água (m 2/s); t = tempo (s);

n = número de terminos de la série; L = dimensión característica (m).

Para ello los valores experimentales de humedad (%) y tiempo de secado se ingresará internamente en el programa, en la ventana se ingresará el valor inicial de la difusividad por el tiempo (Dt), tolerancia, espesor del alimento a secar (L), humedad inicial del alimento (Xo) y humedad de equilibrio (Xe). Datos experimentales de secado: Tiempo (min):0 5 10 15 20 25 30 35 40 45 50 Humedad (%):72 69,54 66,77 63,69 59,38 55,69 51,38 48 44,92 40,92 37,23 Elaboración del programa Formulario




Programa function pushbutton1 _Callback(hObject, eventdata, handles) % hObject handle to pushbutton1 (see GCBO) t=str2num(get(handles.edit1,'string'));% valores de t M=str2num(get(handles.edit2,'string'));% valores de M M0=str2double(get(handles.edit3,'string'));% Humedad inicial Meq=str2double(get(handles.edit4,'string'));% Humedad de equilibrio L=str2double(get(handles.edit5,'string')); % Espesor del alimento a secar tol=str2double(get(handles.edit6,'string')); % Error n=length(t); %Limpiar tabla antes de mostrar resultado set(handles.tabla,'Data',{}) %Limpiar tabla, grafico en caso de que antes se haya graficado hold on %Lectura de datos de tiempo t for k=1:n MR = (M(k) - Meq) / (M0 - Meq); tiempo=t(k); % METODOS NUMERICOS DE NEWTON RAPSON x1 = 0.0000001; it=0; fxp=1; while abs(fxp)>tol fx1 = 0.810566 * exp(-9.869651 * x1 / L ^ 2) + 0.090063 * exp(-88.826855 *x1 / L ^ 2) + 0.032423 * exp(-246.741264 * x1 / L ^ 2) + 0.016542 * exp(-483.612877 * x1 / L ^ 2) + 0.010007 * exp(-799.441695 * x1 / L ^ 2) + 0.006699 * exp(-1194.227718 * x1 / L ^ 2) + 0.004796 * exp(-1667.970945 * x1 / L ^ 2) + 0.003603 * exp(-2220.671376 * x1 / L ^ 2) + 0.002805 * exp(-2852.329012 * x1 / L ^ 2) + 0.002245 * exp(3562.943852 * x1 / L ^ 2)-MR; gx1 = -(8 * exp(-9.869651 * x1 / L ^ 2)) / L ^ 2 - (8 * exp(-88.826855 *x1 / L ^ 2)) / L ^ 2 - (8 * exp(-246.741264 * x1 / L ^ 2)) / L ^ 2 - (8 * exp(-483.612877 * x1 / L ^ 2)) / L ^ 2 - (8 * exp(-799.441695 * x1 / L ^ 2)) / L ^ 2 - (8 * exp(-1194.227718 * x1 / L ^ 2)) / L ^ 2 - (8 * exp(-1667.970945 * x1 / L ^ 2)) / L ^ 2 - (8 * exp(-2220.671376 * x1 / -------------------------------------------------------------------------------------------------------------------------------------------------------51 Alberto L.Huamaní Huamaní



L ^ 2)) / L ^ 2 - (8 * exp(-2852.329012 * x1 / L ^ 2)) / L ^ 2 - (8 * exp(-3562.943852 * x1 / L ^ 2)) / L ^ 2; xp = x1 -(fx1/gx1); fxp = 0.810566 * exp(-9.869651 * xp / L ^ 2) + 0.090063 * exp(-88.826855 * xp / L ^ 2) + 0.032423 * exp(-246.741264 * xp / L ^ 2) + 0.016542 * exp(-483.612877 * xp / L ^ 2) + 0.010007 * exp(-799.441695 * xp / L ^ 2) + 0.006699 * exp(-1194.227718 * xp / L ^ 2) + 0.004796 * exp(-1667.970945 * xp / L ^ 2) + 0.003603 * exp(-2220.671376 * xp / L ^ 2) + 0.002805 * exp(-2852.329012 * xp / L ^ 2) + 0.002245 * exp(3562.943852 * xp / L ^ 2)-MR; x1 = xp; it= it+1; end % CALCULOS X(k) = xp; Dt=X(k); % CALCULO DE DIFUSIVIDAD Dif = Dt/tiempo; %DIFUSIVIDAD EFECTIVA PROMEDIO SUMA=0; SUMA = (SUMA + Dif)/(n-1); % MOSTRARA DATOS EN TABLA valores ={tiempo M(k) MR Dif}; temp=get(handles.tabla,'data'); valoresNuevos=[valores;temp]; set(handles.tabla,'Data',valoresNuevos) hold on %Mostrando respuesta en edit con formato coma flotante a 6 cifras decimales respuesta=sprintf('%0.16f',SUMA); set(handles.respuesta,'string',respuesta); %Grafico de humedad en funcion del tiempo axes(handles.axes1) plot(tiempo,M(k),'r*') title('CINETICA DE SECADO DE HUMEDAD'); xlabel('Tiempo (min)'),ylabel('Humedad (g agua/100 g ms)') grid on; hold on %Grafico de razon de humedad axes(handles.axes2) plot(tiempo,MR,'b*') title('RAZON DE HUMEDAD '); xlabel('Tiempo (min)'),ylabel('MR') grid on; hold on -------------------------------------------------------------------------------------------------------------------------------------------------------52 Alberto L.Huamaní Huamaní



%Grafico de Difusividad en funcion del tiempo axes(handles.axes3) plot(tiempo,Dif,'k*') title('DIFUSIVIDAD '); xlabel('Tiempo (min)'),ylabel('Difusividad (m2/s)') grid on; k=k+1; end function pushbutton2 _Callback(hObject, eventdata, handles) close Solución



5.


Ejercicio aplicado 5: Aplicación de un modelo cinético de reacciones químicas heterogéneas gas-solido al secado por lecho fluidizado de cubos de papa.

A. Desarrollo de la ecuación integrada del modelo El paso del agua desde el centro húmedo hasta el seno del fluido de secado implica dos resistencias: la resistencia a la transferencia de masa desde el centro húmedo hasta la superficie de la partícula y la resistencia a la transferencia de masa desde la superficie hasta el seno del fluido. En todo momento, las condiciones del centro húmedo son idénticas a las condiciones iniciales y las condiciones de la coraza lo son a las condiciones en el equilibrio (Levenspiel, 1976). La velocidad de transferencia de agua desde el centro húmedo hasta la superficie de la partícula a través de la coraza seca está dada aproximadamente por: (1 / 4  r 2 )(  dM / dt )  kc(  dC / dr)  cte.

(1)

Donde: kc es el coeficiente de difusión efectivo del agua a través de la coraza, M y C son la masa y la concentración volumétrica del agua en el material y r es la posición de un punto de la partícula en coordenadas esféricas. Integrando la ecuación (1) entre los límites Cs,Rc y Ci,R se tiene: ( dM / dt )



R

Rc

dr / r 2  4 kc



Ci

Cs

dC

(2)

De donde: (  dM / dt ) 

4  RcRkc ( Cs  Ci ) ( R  Rc )

(3)

y la velocidad de transferencia de agua desde la superficie de la partícula hasta el seno del fluido está dada por: (1 / 4  R 2 )(  dM / dt )  hc (Ci  Cg )

(4)

Donde hc es el coeficiente de película. De donde: (  dM / dt )  4  R 2 hc(Ci  Cg )

(5)

Combinando las ecuaciones (3) y (5) y eliminando Ci se obtiene: ( dM / dt ) 

4  R 2 (Cs  Cg ) (1 / hc)  ( R / kc ) ( R  Rc ) / Rc

(6)

Rearreglando esta ecuación se obtiene una ecuación similar a las usadas en la cinética de reacciones químicas heterogéneas sólido-fluido: -------------------------------------------------------------------------------------------------------------------------------------------------------54 Alberto L.Huamaní Huamaní


(  dM / dt ) 

4  R 2hc(Cs  Cg ) 1  ( Rhc / kc ) ( R  Rc ) / Rc


(7)

Por otro lado la cantidad de agua a eliminar por unidad de volumen es: Ch  ( Mo  Me ) / Vo

(8)

En el tiempo t cuando el volumen del centro húmedo es Vc: M  ChVc  Ch ( 4  Rc 3 / 3)

(9)

La disminución del radio húmedo correspondiente a una disminución en la cantidad de agua a eliminar se obtiene de la derivada de la ecuación (9): dM / dRc  4 Rc 2 Ch

(10)

y dado que: dRc / dt  ( dM / dt) / ( dM / dRc)

(11)

Sustituyendo las ecuaciones (7) y (10) en la (11) se obtiene: ( dR / dt ) 

R (Cs  Cg ) 2

Ch ( Rc 2 / hc)  ( R / kc )(R ( R  Rc ))

(12)

en donde el radio R que considera el encogimiento de la partícula está relacionado con las condiciones iniciales y de equilibrio por la ecuación: R  Re 3  1  (Re/ Ro ) 3 Rc 3

1/ 3

(13)

Separando variables e integrando la ecuación (12) entre Ro,0 y Rc,t se tiene: Cs  Cg Ch

t

 dt  

Ch

2 ( Rc / hc)  ( Rc / kc) R( R  Rc)

R

Rc

0

( Cs  Cg )t

Ro



Ro  R hc 1  (Re/ Ro ) 3



2

Ro 2  Rc 2

2 kc



dRc

(14)

Ro 2  R 2

2 kc 1  (Re/ Ro )3

(15)

Reacomodando la ecuación (15): ( Cs  Cg ) 1  (Re/ Ro ) 3 t Ch ( Ro  R )

 (1 / hc)  ( 1 / kc)

R 2  Rc 2  ( Ro 2  Rc 2 )(Re/ Ro ) 3

2 ( Ro  R )

(16)

La cual es la expresión matemática del modelo propuesto y que adopta la forma de una recta. Donde: -------------------------------------------------------------------------------------------------------------------------------------------------------55 Alberto L.Huamaní Huamaní


( Cs  Cg ) 1  (Re/ Ro ) 3 t


es la variable dependiente

(17)

es la variable independiente

(18)

(1 / kc)

es la pendiente

(19)

(1 / hc)

es la ordenada al origen

(20)

Ch ( Ro  R ) R 2  Rc 2  ( Ro 2  Rc 2 )(Re/ Ro ) 3

2( Ro  R)

La ecuación (16) puede interpretarse en términos de resistencias: Resistencia total = Resistencia de la película + Resistencia de la coraza seca Donde: Resistencia total = variable dependiente. Resistencia de la película = ordenada al origen. Resistencia de la coraza seca = pendiente x variable independiente. B. Calculo de los coeficientes de transferencia del modelo Usando el método de mínimos cuadrados se calcularon los valores de la pendiente y de la ordenada al origen de la recta de mejor ajuste a las parejas de valores dados por las ecuaciones (18) y (17), evaluados a partir de los datos experimentales. De los valores de la pendiente y de la ordenada al origen se obtienen los valores de los coeficientes de transferencia kc y hc según las ecuaciones (19) y (20). C. Expresiones para estimar los coeficientes de transferencia Para estimar los coeficientes de transferencia hc y kc se obtuvieron las siguientes expresiones: hc = 2,946 (10) -4 (VR/C)1,051 e (-1955/T)

(21)

kc = 1,884 (10) -6 e (-3610/T)

(22)

Con coeficientes de correlación ajustados de 0,86 y 0,85 respectivamente. HC=(2,946*10^4)*exp(-1954,87/TEMP)*(VEL*RAD/CAM)^1,051 ; hc = 2,946 (10)-4 (VR/C) 1,051 e (-1955/T) Estructura del programa y métodos numéricos D. Programa para la simulación de la cinética de secado y del encogimiento de cubos de papa en un secador de lecho fluidizado. -------------------------------------------------------------------------------------------------------------------------------------------------------56 Alberto L.Huamaní Huamaní



Una vez estimados los valores de hc y kc se puede usar la ecuación (15) para predecir el radio del centro húmedo Rc a cualquier tiempo t. Para despejar Rc de dicha ecuación se utiliza el método numérico de aproximaciones sucesivas de Newton de primer orden (Luthe y col., 1988): Rcactual  Rcanterior 

F ( Rc ) anterior

(23)

F ' ( Rc ) anterior

Donde: 3   C s  C g    Re  3    R 1 1  2 2 2 2 e  R  RC    RO  Rc    f ( RC )   Ro  R    1     * t 2kc  hc R C  h  o      Ro  















 C s  C g   * A * t C h  

f ( x)  H  Ro  R   K R 2  x 2  B * RO  x 2   2

Y 2   R  3   1  R 2  1     R e C c        2   f ' ( x )   1   Rc          Ro    hc  R  kc   R  

Para todos los cálculos se toma como valor inicial Rc=Ro. El programa desarrollado en este trabajo, calcula el radio del centro húmedo Rc (cm) por el método de Newton, el radio de la partícula R (cm), el peso total de la muestra M (g) y la humedad en base seca X (g de agua/g s.s.) para diferentes tiempos. Requiere conocer la densidad de la muestra inicial o (g/cm3), la densidad del sólido seco ss (g/cm3), el radio equivalente inicial y en el equilibrio; Ro y Re (cm), el peso inicial, en el equilibrio y del sólido seco; Mo, Me, Mss (g), y los valores estimados de los coeficientes de transferencia en la película hc (m/s) y de difusión a través de la coraza kc (m 2/s) calculados por las ecuaciones (21) y (22) a partir de la velocidad del aire VEL (m/s), de la altura de la cama del lecho a fluidizar CAM (m), de la temperatura del aire Temp (K) y del radio de la partícula Rad (m). Rangos de aplicación El programa de simulación del proceso de secado de cubos de papa por lecho fluidizado puede ser aplicado dentro de los siguientes rangos: Aristas de los cubos de papa Radios equivalentes de los cubos Temperatura del aire de secado Velocidad del aire de secado Altura del lecho a fluidizar

0,5 a 1,2 cm. 0,4 a 0,76 cm. 50º a 100º C. 4 a 8 m/s. 7 a 13 cm.

Datos de entrada El programa requiere la siguiente información: -------------------------------------------------------------------------------------------------------------------------------------------------------57 Alberto L.Huamaní Huamaní



DO = Densidad inicial (g/cm 3) DSS = Densidad del sólido seco (g/cm 3) RO = Radio equivalente inicial (cm) RE = Radio equivalente en el equilibrio (cm) MO = Peso inicial (g) ME = Peso en el equilibrio (g) MSS = Peso del sólido seco (g) XO = Humedad inicial (g de agua/g s.s.) TEMP = Temperatura del aire de secado (ºC) VEL = Velocidad del aire de secado (m/s) CAM = Altura de la cama del lecho a fluidizar (m) Salida Para cada instante de tiempo se muestran los siguientes resultados: Numero de iteraciones para alcanzar una exactitud de 0,0001 Tiempo (min) Radio simulado de la partícula (cm) Radio simulado del centro húmedo (cm) Peso de la partícula (g) Humedad de la partícula (g agua/g s.s.) Elaboracion del programa Formulario

Programa function pushbutton1 _Callback(hObject, eventdata, handles) % hObject handle to pushbutton1 (see GCBO) % "PROGRAMA DE SIMULACION DEL PROCESO DE SECADO" % SOLO PUEDE SER APLICADO A CUBOS DE PAPA, SECADO EN LECHO FLUIDIZADO, % "======================================================" -------------------------------------------------------------------------------------------------------------------------------------------------------58 Alberto L.Huamaní Huamaní



% "CINETICAS DE SECADO" % "======================================================" t=str2num(get(handles.edit17,'string')); DO =str2double(get(handles.edit1,'string'));%Densidad inicial (g/cm3) 1.06 DSS =str2double(get(handles.edit2,'string'));% Densidad del solido seco (g/cm3) 1.245 RO =str2double(get(handles.edit3,'string')); %Radio equivalente inicial (cm) 0.4 RE =str2double(get(handles.edit4,'string')); %Radio equivalente en el equilibrio (cm) 0.216 MO =str2double(get(handles.edit5,'string')); %Peso inicial (g) 301.3 ME =str2double(get(handles.edit6,'string')); %Peso en el equilibrio (g) 61.7 MSS =str2double(get(handles.edit7,'string')); %Peso del solido seco (g) 57.254 XO =str2double(get(handles.edit8,'string')); %Humedad inicial (g de agua/g s.s.) 4.26 CAM =str2double(get(handles.edit9,'string')); %Altura de la cama del lecho (m) 0.07 TEMP=str2double(get(handles.edit10,'string')); %Temperatura del aire de secado (C) 50 VEL =str2double(get(handles.edit11,'string')); %Velocidad del aire de secado (m/s) 4 n=length(t); %n=19; % CALCULO DE LOS COEFICIENTES DE TRANSFERENCIA INTERNO Y EXTERNO RAD=RO/100; TEMP=TEMP+273; HC=.0002946*exp(-1954.87/TEMP)*(VEL*RAD/CAM)^1.051 ; KC = 1.884*(10^-6)*exp(-3610/TEMP); % CONVERSION DE UNIDADES DE LOS COEFICIENTES HC=HC*6000; KC=KC*360000; % VARIABLES AUXILIARES PARA SIMPLIFICAR EXPRESIONES CS=XO*DSS; CH=DO*(MO-ME)/MO; H=1/HC; K=1/(2*KC); J=CS/CH; B=(RE/RO)^3; A=1-B; C=RE^3; D=RE^2; %Lectura de datos de tiempo t %Limpiar tabla antes de mostrar resultado set(handles.tabla,'Data',{}) %Limpiar tabla, grafico en caso de que antes se haya graficado hold on for k=1:n tiempo=t(k); % METODOS NUMERICOS DE NEWTON RAPSON x1=0.20; it=0; fxp=1; error=0.0001; while abs(fxp)>error -------------------------------------------------------------------------------------------------------------------------------------------------------59 Alberto L.Huamaní Huamaní



fx1= H*(RO-((C+A*x1^3)^(1/3)))+K*(((C+A*x1^3)^(1/3))^2-(C/RO)-A*x1^2)(CS/CH)*A*tiempo; gx1=-A*(H*((x1^2)/((C+A*(x1)^3)^(1/3))^2)+2*K*(x1(x1^2/((C+A*(x1)^3)^(1/3))))); xp = x1 -(fx1 / gx1); fxp= H*(RO-((C+A*xp^3)^(1/3)))+K*(((C+A*xp^3)^(1/3))^2-(C/RO)-A*xp^2)(CS/CH)*A*tiempo; x1 = xp; it= it+1; end % CALCULOS XR(k)=xp; RR= XR(k); % CALCULO DEL RADIO DE LA PARTICULA" R=((C+A*(RR)^3)^(1/3)); % CALCULO DEL PESO M=ME+(MO-ME)*((RR)/RO)^3; %"CALCULO DE LA HUMEDAD BASE SECA X=(M-MSS)/MSS; % MOSTRARA DATOS EN TABLA valores ={tiempo xp R M X}; temp=get(handles.tabla,'data'); valoresNuevos=[valores;temp]; set(handles.tabla,'Data',valoresNuevos) hold on %Grafico de humedad en funcion del tiempo axes(handles.axes1) plot(tiempo,X,'r*') title('CINETICA DE SECADO'); xlabel('Tiempo (min)'),ylabel('Humedad base seca (g agua/100g ms)' ) grid on; hold on %Grafico de peso en funcion del tiempo axes(handles.axes2) plot(tiempo,M,'b*') title('CINETICA DE SECADO'); xlabel('Tiempo (min)'),ylabel('Peso (g)') grid on; hold on %Grafico de radio de particula en funcion del tiempo axes(handles.axes3) plot(tiempo,R,'k*') title('CINETICA DE SECADO'); xlabel('Tiempo (min)'),ylabel('Radio de particula (m)') grid on; k=k+1; end function pushbutton2 _Callback(hObject, eventdata, handles) close -------------------------------------------------------------------------------------------------------------------------------------------------------60 Alberto L.Huamaní Huamaní



Resultado

BIBLIOGRAFIA Alamilla, B. L. 1990. Simulación de la operación de secado de vegetales basada en un estudio de deshidratación por lecho fluidizado. Tesis de Maestría en Ciencias (Alimentos). Escuela Nacional de Ciencias Biológicas, IPN, México, D. F. Alamilla, L., Gutiérrez, G., Hernández, H. y Santiago, P. 1991. Estudio semifundamental del secado de papa en lecho fluidizado. Tec. Aliment. 25(4):24-29. Brown, G. G., Foust, A. S., Katz, D. V., Schneidewind, R., White, R. R., Wood, W. P., Brownell, L. E., Martin, J. J., Williams, G. B., banchero, J. T., and York, J. L., 1965. Fluidización de sólidos. Cap. 20, En Operaciones Básicas de la Ingeniería Química. p. 285-288. Ed. Marín, Barcelona. De Baun, R. M. 1959. Response surface design for three factors at three levels. Technometrics. 1(1):1-8. Félix, A. B., Robles, R. R. y Santiago, P. T. 1989. Estudio de Ingeniería para la deshidratación de papa por lecho fluidizado. Memorias de AMIDIQ.,México Levenspiel, O., 1976. Solid-fluid reactions. Ch. 12, In Chemical Reaction Engineering. 2da. ed., p. 357-377. Wiley International Edition, N. Y.

6. Ejercicio aplicado 5: Aplicación en mecánica de fluidos: Cuando un fluido no Newtoniano pasa por una tubería de sección circular bajo cierto régimen, se genera una pérdida de carga debido a la fricción que puede ser estimada de acuerdo a la ecuación de Darcy: -------------------------------------------------------------------------------------------------------------------------------------------------------61 Alberto L.Huamaní Huamaní


h f  f

L * v


2

2 * D * g

Dónde: hf: pérdida de carga F: Factor de fricción(s/u) L: longitud de la tubería D; diámetro de la tubería v: velocidad promedio del fluido dentro de la tubería g: aceleración de la gravedad Por otro lado, el factor de fricción”f” puede ser estimado mediante la ecuación de Colebrook:  1 2.51   2 log  , D e 3 . 7 ( / ) g f 4 f Re 4   1

n

 4n  D n v 2n Re g     n 1  3n  1  8 k La velocidad de flujo podemos estimar a través de: v

Q A



f m

 . A

Desarrollar un programa que calcule el caudal “Q”, el área transversal de la tubería “A”, la velocidad de flujo “v”, el número de Reynolds generalizado “Reg”, el factor de fricción “f”, y la pérdida de carga en la tubería “h f ”, ”, a partir de los datos de: Flujo másico: f m=2,22 kg/s Densidad del fluido: den=1165kg/m 3 Indice reológico del fluido: n=0,65 Coeficiente de consistencia del fluido: m=4,43 Rugosidad absoluta de la tubería: e=0,000005 Diámetro de la tubería: D=0,0343m D =0,0343m Longitud del fluido: L=200m

Solución Se realizaran los cálculos en el siguiente orden:

1.

Q 

f m

v 



Q A

n

 4n  D n v 2n Re g     n1  3 1 8 n k   2. -------------------------------------------------------------------------------------------------------------------------------------------------------62 Alberto L.Huamaní L.Huamaní Huamaní Huamaní



 1 2.51   2 log   4 f  3.7( D / e) Re g 4 f  se realizará por métodos métodos numéricos numéricos 1

3.

4. hf  f (a)

L * v

2

2 * D * g

Programa en matlab archivo m

Primero creamos el programa, introduciendo los datos del ejercicio, luego los cálculos previos que debe debe hacer el el programa para el cálculo cálculo de las constantes, constantes, luego luego aplicaremos la solución por el método de falsa posición. %PROBLEMA DE PERDIDA DE CARGA EN TRANSPORTE DE FLUIDO %SOLUCIÓN DE ECUACIONES NO LINEALES %"METODO DE REGLA FALSA" % Setiembre 2010 % Ingeniero Alberto HUAMANI HUAMANI clear all clc flujom = 2.22;den = 1165;n = 0.65;k = 4.43;e = 0.000005;Diam = 0.0343;lon = 200; %Calculos Q =flujom/den; Area =3.1416 * (Diam) ^ 2/ 4; vel = Q / Area; A = (4 * n / ((3 * n) + 1)) ^ n; Reg = A * ((Diam ^ n) * (vel ^ (2 - n)) * den) / (k * 8 ^ (n - 1)); b1 = 1 / (3.7 * (Diam (Diam / e)); b2 = 2.51 / Reg; Reg; % Crear la salida disp(' disp(' ')') disp(' disp(' Alberto HUAMANI HUAMANI') HUAMANI') disp(' disp(' Ingeniería en Industrias Alimentarias') Alimentarias' ) disp(' disp(' Metodos matemáticos en industrias alimentarias' ) disp(' disp(' ')') fprintf('\n' fprintf('\n'); ); x1=0.001; x2=0.1; error=0.0001; it=1; fxr=1; while abs(fxr)>error while abs(fxr)>error fx1=-1-2*(4*x1)^0.5*log(b1+b2*(4*x1)^-0.5); fx2=-1-2*(4*x2)^0.5*log(b1+b2*(4*x2)^-0.5); xr=x1-(fx1*(x2-x1)/(fx2-fx1)); -------------------------------------------------------------------------------------------------------------------------------------------------------63 Alberto L.Huamaní L.Huamaní Huamaní Huamaní



fxr=-1-2*(4*xr)^0.5*log(b1+b2*(4*xr)^-0.5); if (fxr*fx2)>0; (fxr*fx2)>0; x2=xr; else x1=xr; end it=it+1; hf = (xr * lon * (vel^2))/(2 * Diam * 9.81); end fprintf('\n fprintf('\n El coeficiente de fricción f es: %8.6f \n' ,xr); fprintf('\n fprintf('\n El número de iteraciones es: %4.0f \n' ,it); fprintf('\n fprintf('\n El caudal Q (m3/h) es: %8.4f \n' \n',Q); ,Q); fprintf('\n fprintf('\n La velocidad v(m/s) es: %8.6f \n',vel); \n' ,vel); fprintf('\n fprintf('\n El reynolds generalizado Reg es: %8.6f \n',Reg); \n' ,Reg); fprintf('\n fprintf('\n La perdida de carga (m agua) es: %8.6f \n' ,hf); Compilación: ejecutamos el programa realizado Alberto HUAMANI HUAMANI Ingeniería en Industrias Alimentarias Metodos matemáticos en industrias alimentarias El coeficiente de fricción f es: 0.010160 El número de iteraciones es: 8 El caudal Q (m3/h) es: 0.0019 La velocidad v(m/s) es: 2.062281 El reynolds generalizado Reg es: 148.820555 La pérdida de carga (m agua) es: 12.841765

7. Ejercicio 6: Secado a bajas presiones y bajas temperaturas -------------------------------------------------------------------------------------------------------------------------------------------------------64 Alberto L.Huamaní L.Huamaní Huamaní Huamaní



Un proceso alternativo entre el secado convencional y la liofilización es el secado a baja presión y baja temperatura (BPBT) que ha sido planteado por King y col. (1989), el cual compensa aspectos de costo y calidad. El método de secado BPBT consiste en evaporar el agua reduciendo simultáneamente la presión del ambiente y la temperatura del material sin llegar a congelarlo, llevándolo cerca del punto triple. El mismo autor, ha evaluado la calidad de productos secados con este método, obteniendo características cercanas a aquellas obtenidas por liofilización. Se puede establecer una analogía entre el fenómeno de secado y los fenómenos de transporte que ocurren en las partículas sólidas en las reacciones químicas heterogéneas sólido-fluido. Existen dos modelos que se aplican a las reacciones químicas heterogéneas sólido-fluido y que pueden aplicarse al secado: el modelo de conversión progresiva y el modelo del núcleo sin reaccionar (Mercado y Gutiérrez, 1995) El modelo del núcleo sin reaccionar propone que la reacción ocurre primero en la superficie exterior de la partícula sólida, después se forma un frente de reacción definido (nítido, no difuso) que se mueve a un punto más interno de la partícula (Figura 3), dejando una capa de cenizas y un núcleo en el que el reactante sólido no ha reaccionado; este núcleo irá disminuyendo de tamaño a medida que transcurre la reacción (Levenspiel, 1986) En el secado de alimentos ocurren cambios físicos en la superficie de la partícula y el frente de secado avanza hacia su interior formándose una capa de material seco la cual debe atravesar el agua para continuar el proceso de secado. Existen pocas investigaciones en las que se trata de explicar el proceso de secado a través del uso de la teoría de la cinética química y la mayor parte de los estudios reportados se basan en el análisis de la curva de secado. El modelo del núcleo sin reaccionar, tomando en cuanta el encogimiento, fue aplicado por Mercado y Gutiérrez (1995) al secado por lecho fluidizado de cubos de papa. En esta investigación se utilizará el modelo del núcleo sin reaccionar, de partícula de tamaño decreciente, es decir, considerando el encogimiento, para simular la cinética de secado de materiales en forma de placa y se aplicará al método de secado a baja presión y baja temperatura. Los datos experimentales que se utilizan en la verificación del modelo se obtuvieron de la literatura (Zazueta, 1994). Corresponden a los datos de 9 diferentes cinéticas de secado de placas de puré de papa (Solanum tuberosum) de la variedad Alpha, deshidrata por el método de secado a baja presión y baja temperatura. Los factores considerados son presión en la cámara de secado (0.67, 1.00 y 1.33 kPa) y espesor de la placa (0.26, 0.52 y 0.78 mm). El área de cada placa es de 0.01904 m2. Dado que el modelo del núcleo sin reaccionar requiere la incorporación de un coeficiente interno de transferencia de masa (D i), se propuso que este parámetro podría seguir la misma tendencia matemática que la difusividad efectiva (D eff ) de la segunda ley de Fick. Esta tendencia se determina en forma de ajuste. Obtención de la ecuación integrada de velocidad de secado -------------------------------------------------------------------------------------------------------------------------------------------------------65 Alberto L.Huamaní Huamaní



La ecuación integrada de velocidad de secado es la relación que describe el contenido de humedad en función del tiempo, X= f (t). Para obtenerla se desarrolla el modelo que representa la cinética de secado de materiales en forma de placas planas, considerando su encogimiento. El modelo original fue desarrollado por Yagi y Kunii (1955) y se encuentra descrito en Levenspiel (1986). En el presente trabajo al modelo del núcleo sin reaccionar se le denomina modelo del núcleo húmedo, cuando se aplica al proceso de secado. Para obtener el modelo, se considera que durante la deshidratación, se presentan en forma de resistencias en serie dos etapas de movimiento de agua. Una que se lleva a cabo a través de la película gaseosa que rodea a la placa y la otra a través de la capa de sólido seco que se va formando con el tiempo. La ecuación que define la velocidad de transferencia de agua desde la superficie de la placa hacia el aire, a través de la película gaseosa, es la siguiente: 1  dma 



A 

Dónde:

(1)

área de transferencia de agua (m 2) masa de agua (kg) tiempo (s) coeficiente local de transferencia de masa en la película gaseosa (m/s) concentración de agua en la superficie de la placa (kg/m 3) concentración de agua en el seno del gas (kg/m 3)

= = = = = =

A ma t k g C s C g

  K g C s  C g 

dt 

La ecuación que describe el movimiento del agua desde el núcleo húmedo hacia la superficie de la placa, a través de la capa de sólido seco, es: 1  dma 

 dC     D    A  dt  dz   i

Donde: A ma t Di

C

= = = = =

(2)

área de transferencia de agua (m 2) masa de agua (kg) tiempo (s) coeficiente interno de transferencia de masa (m 2/s) concentración de agua en la capa de sólido seco (kg/m 3)

Si el secado es controlado por condiciones externas, la ecuación (1) será la que describa el proceso. Si el sólido representa la mayor resistencia, la ecuación (2) tendrá un significado más alto. En realidad el modelo del núcleo húmedo presenta como característica ventajosa, desde el punto de vista fenomenológico, que considera una sucesión gradual de resistencia. En un principio la fase gaseosa tendrá mayor importancia y posteriormente el sólido será el factor controlante. Esto es un reflejo cercano a lo que sucede durante el proceso de secado de muchos alimentos y materiales biológicos en general. Se seleccionó el modelo exponencial para describir la difusividad efectiva en función del tiempo debido a que es altamente significativo y a que su manejo dentro del simulador es -------------------------------------------------------------------------------------------------------------------------------------------------------66 Alberto L.Huamaní Huamaní



sencillo. Por lo anterior se supone que Di seguirá un comportamiento matemático análogo al de Deff del tipo:

Di  ae bt

(3)

Obtención de la ecuación integrada de velocidad de secado La ecuación de velocidad de transferencia se convierte en: 1  dm     k gC A  dt 

Lc  L0 

k gC c C h

(4)

t

(5)

Que es el modelo que explica la disminución del espesor de la placa cuando la difusión de agua a través de la película gaseosa controla la velocidad de secado. 1  dm 

D C

  i c A  dt  L  Lc Lc  L0 

(6)

2 L0C c ebt  1 LeC h

b

(7)

Donde: Lo : el espesor inicial Le :el espesor de la placa en equilibrio Cc: la concentración de agua en el núcleo húmedo Ch : la cantidad de agua a evaporar por unidad de volumen Que es el modelo que explica la disminución del espesor de la placa cuando la difusión de agua a través de la capa de sólido seco controla la velocidad de secado. En las expresiones anteriores, aparece una resistencia combinada: capa de sólido seco y de la película gaseosa cuya importancia relativa varía a medida que se realiza el secado, controlando así el proceso global de deshidratación de la placa. Debe tomarse en cuenta que ambas resistencias actúan en serie y que ambas son lineales respecto a la concentración de agua y sólidos. Por consiguiente, se pueden combinar directamente las resistencias individuales de ambas etapas de acuerdo a: A k C  dm     dt  1  k  L  L  g

c

g

c

(8)

Di

El grupo adimensional k g(L-Lc)/Di es el número de Damkohjer (Da) para la capa de sólido seco (Levenspiel, 1986). -------------------------------------------------------------------------------------------------------------------------------------------------------67 Alberto L.Huamaní Huamaní


k g ( L  L c )

Da 

Di


(9)

En este caso Da representa el cociente de la resistencia por difusión del líquido en el interior de la capa de sólido seco y de la resistencia a la transferencia de agua de la película gaseosa. La solución representa el modelo del núcleo húmedo para placas planas que presentan encogimiento, la cual describe la disminución en tamaño del núcleo húmedo a medida que la placa se va secando, es decir, lc = f(t):

e

bCh kgCc

( Lc  L0 )

 a 

e

bt

Lek g  k gC c  Lek g 2 C c   L 0  Lc    a  h L0  bC h  bC L 0

(10)

%Cálculo de simplificación: Ch, Le, Cc, alfa y beta Di = a * exp(b * t); alfa = b * Ch / (kg * Cc); beta = Le * kg / Lo; Din = mo / (Ar * Lo); Ch = Din * (mo - me) / mo Le = me / (de * Ar); Dss = mss / (Ar * Le) xo = (mo - mss) / mss; Cc = Dss * xo: p = exp(alfa * (lcv - Lo)) q = Di + beta * (Lo - lcv + 1 / alfa) f = p * q - a - beta / alfa Masa de agua con la variación del espesor m

me  m0  me  lc l0

(11)

Humedad en base húmeda es:

X 

m  mss mss

(12)

Programa para la simulación de la cinética de secado para placas de puré de papa En la figura 1 se presenta el diagrama de flujo del programa que simula la cinética de secado para placas de puré de papa. -------------------------------------------------------------------------------------------------------------------------------------------------------68 Alberto L.Huamaní Huamaní



Parámetros de simulación (propuestos)

Dimensiones y masas de la placa

Cálculo de Ch, Le y Cc

Sí

No

Método de Newton de 2° orden

Figura 1: Diagrama de flujo del programa que simula la cinética de secado para placas de puré de papa. Los parámetros de simulación son el coeficiente local de transferencia de masa en la película gaseosa (k g) y los parámetros a y b que definen el coeficiente interno de transferencia de masa (D i) como función exponencial del tiempo. Los valores experimentales que se requieren son: A: el área, Lo : el espesor inicial mo : la masa inicial de agua en la placa, me: la masa de agua en equilibrio pe= la densidad en el equilibrio mss: la masa de sólido seco -------------------------------------------------------------------------------------------------------------------------------------------------------69 Alberto L.Huamaní Huamaní



Se requiere calcular: Ch : la cantidad de agua a evaporar por unidad de volumen Le :el espesor de la placa en equilibrio Cc: la concentración de agua en el núcleo húmedo Para calcular el espesor del núcleo húmedo (L c) a cualquier tiempo (k), se despeja de la ecuación el valor de L c utilizando el método numérico de aproximaciones sucesivas de Newton de segundo orden (James y col., 1979; Luthe y col., 1988) cuya ecuación de recurrencia es: f (lco )

lc  lco  f ' (lco )  bC h  Lc  L0 

f  Lc   e

k g C c

f (lco )

2 f ' (lco )

 bt Le k g ae  L0 

f ' ' (lco )

(13)

k C  L k C   L0  Lc  g c   a  e g c bC h  bC h L0  2

(14)

2

Le k g C c  1   f  Lc   p  Di    L0  Lc    a    bC h L0  

f  Lc   p * q  a 

 

bC h k g C c

 

: función matemática

 

p  exp *  x  L0 

f  Lc   p * q  a 

 

Le k g L0

Di  a expb * 60 * t 

 

q  Di   *  L0  x 

1 



 

: %la funcion matematica

Dónde: F’(l c) y F’’(l c) son la primera y segunda derivadas de F(lc), respectivamente.

Formulario




Programa function pushbutton1 _Callback(hObject, eventdata, handles) kg=str2double(get(handles.edit1,'string')); a=str2double(get(handles.edit2,'string')); b=str2double(get(handles.edit3,'string')); Ar=str2double(get(handles.edit4,'string')); Lo=str2double(get(handles.edit5,'string')); mo=str2double(get(handles.edit6,'string')); me=str2double(get(handles.edit7,'string')); mss=str2double(get(handles.edit8,'string')); de=str2double(get(handles.edit9,'string')); t=str2num(get(handles.edit10,'string')); % para varios valores n=length(t); %Calculus: Ch, Le, Cc, alfa y beta Din = mo/(Ar*Lo); Ch = Din*(mo-me)/mo; Le = me/(de*Ar); Dss = mss/(Ar * Le); xo = (mo - mss)/mss; Cc = Dss*xo; alfa = b*Ch/(kg * Cc); beta = Le*kg / Lo; %Cálculo del contenido de humedad en función del tiempo for k=1:n % t(k)*60; % conversión de tiempo de min a s Di = a*exp(b*t(k)*60);% Calculo de Di y conversión min a s tiempo = t(k); %METODO DE NEWTON RAPHSON DE SEGUNDO ORDEN lcv = Lo; %Aproximación inicial de la raíz e=0.000001; ea=1000; -------------------------------------------------------------------------------------------------------------------------------------------------------71 Alberto L.Huamaní Huamaní



it=0; % iteracion while ea>e p = exp(alfa*(lcv-Lo)); q = Di + beta*(Lo-lcv+1/alfa); f = p*q-a-beta/alfa; % funcion DF = p*(alfa*q-beta); % derivada de funcion DDF = p*(alfa^2*q-2*alfa*beta); % segunda derivada de funcion % lcn = lcv-f/(DF-f*DDF/(2*DF)); lcn=lcv-(f*DF)/((DF^2)-(f*DDF)); % algoritmo de Newton ea= abs(((lcn-lcv)/lcv)*100); lcv = lcn; it=it+1; end %CALCULO DE HUMEDAD EN CADA TIEMPO lc(k)=lcn; Lc=lc(k); m = me+(mo-me)*Lc/Lo;% Calculo de masa seca en cada tiempo X= (m-mss)/mss; % Humedad en cada tiempo % CALCULO DEL ENCOGIMIENTO DEL ESPESOR EN CADA TIEMPO lc(k) = lcn; L(k) = Le+(1-Le/Lo)*lc(k); Dam = kg*(L(k)-lcn)/Di; rpg(k) = 1/(1+Dam); rcss(k) = Dam/(1+Dam); % MOSTRARA DATOS EN TABLA valores ={tiempo m X L(k) }; temp=get(handles.tabla,'data'); valoresNuevos=[valores;temp]; set(handles.tabla,'Data',valoresNuevos) hold on %Grafico de humedad en funcion del tiempo axes(handles.axes1) plot(tiempo,X,'r*') title('CINETICA DE SECADO'); xlabel('Tiempo (min)'),ylabel('Humedad (g agua/ 100g ms)') grid on; hold on % Variacion de masa de agua axes(handles.axes2) plot(tiempo,m,'b*') title('CINETICA DE SECADO'); xlabel('Tiempo (min)'),ylabel('Masa de agua (g)') grid on; hold on % Grafico de variacion de espesor axes(handles.axes3) plot(tiempo,L(k),'g*') title('CINETICA DE SECADO'); xlabel('Tiempo (min)'),ylabel('espesor de nucleo (m)') grid on; -------------------------------------------------------------------------------------------------------------------------------------------------------72 Alberto L.Huamaní Huamaní



k=k+1; end function pushbutton2 _Callback(hObject, eventdata, handles) close Ejecucion del programa

Bibliografía 1. Zazueta, N.J.A. 1994. Estudio del secado a baja presión y baja temperatura y su comparación con el secado en túnel experimental. Tesis de Maestría, Universidad Autónoma de Sinaloa, México. Programa elaborado en lenguaje Visual Basic.



Ejercicios propuestos

Práctica 1-7

EJERCICIOS PROPUESTOS 1.

Problema caso de deshidratación osmótica: Determinacion de los coeficientes de difusion del agua y de los solutos.

La mayoría de los modelos existentes para el estudio de la cinética de deshidratación osmótica se basan en la segunda ley de Fick, que se basa en la ecuación de difusión, el flujo de masa es proporcional al gradiente de concentración entre el sólido y la solución. Crank (1975) presenta una serie de soluciones analíticas para la ley de la Segunda Fick, teniendo en cuenta las diferentes condiciones iniciales y de frontera. En sistema de coordenadas rectangular (x, y, z), la ecuación de difusión se expresa como:

Suponiendo que no hay forma constante y geométrica de una placa infinita con espesor 2L, donde la transferencia de humedad durante la deshidratación osmótica es predominantemente unidireccional y considerando Def constante, la ecuación 2.1 se reduce a:  X  2 X  Def t  y 2 Para la situación de humedad inicial X0 uniforme, despreciando las resistencias externas a la transferencia de materia y considerando la contracción del producto durante la deshidratación despreciable, las condiciones de contorno para la humedad X son: X  Xo

en el t = 0

0
(2.2a)



X  Xe

 X 0  y


en el y = L

t>0

(2.2b)

en el y = 0

t>0

(2.2c)

La solución de la ecuación ecuación 2.2 con las condiciones condiciones de contorno 2.2a, 2.2b e 2.2c es:  ( 2n  1) 2  2   2 exp . D t ef   X 0  X e  n0 ( 2n  1) 2 4 L2   X  X e

8



1

Para el soluto  ( 2n  1) 2  2   2 exp . D t ef   M 0  M e  n 0 ( 2n  1) 2 4 L2   M  M e

8



1

Donde:

X = humedad de la muestra (g); M = masa de sólidos (sacarosa) en la muestra (g); Def1 = difusividade efectiva del agua (m 2/s); Def2 = difusividad efetiva del soluto (m 2/s); t = tiempo (s);

n = número de terminos de la série; L = dimensión característica (m). caso placa

Realizar un programa con dos tablas t ablas de resultado y dos graficos de cinética de difusividad efectiva para el soluto y para el agua del producto. Tabla: Tabla: Valores de la Perdida de Água (PA), Ganancia de Sólidos (GS) y del Adimensional GS/PA durante la cinética de desidratación osmótica con solución de sacarosa (40, 50 y 60°Brix) a 30°C.



2.


Problema: Transporte de un fluido no Newtoniano: Se desea calcular la velocidad y caudal másico másic o con el que circula un puré de manzana que manzana que se transporta a través de una tubería de 15 cm de diámetro interno, a una temperatura de 20ºC. la distancia total que debe recorrer el puré es de 300 m, existiendo entre los puntos de salida y llegada una caída de presión de 250 kPa, estando el punto de llegada a 5 m por encima del punto punto de salida. A la temperatura de trabajo, trabajo, el puré sigue la ley de n potencia, con con un índice índice de consistencia consistencia de de 2,4 Pa.s Pa.s y un índice de comportamiento al flujo de 0,44, siendo su densidad de 1200 kg/m 3.

Solución El número de Reynold crítico correspondiente correspondiente a este flujo se obtiene de la ecuación: 6464n

Re Crítico 

1  3n2  

1    2  n 

 2  n     1 n 

Las pérdidas de energía mecánica se obtienen al aplicar la ecuación de Bernoulli entre los puntos de entrada entrada y salida es: es: Ef 

P 

   z

Se asume una velocidad a partir del cual se evalúa el valor del módulo de Reynolds generalizado n

d   4n  2n ReG  n1   v 8 k  1  3n  n

d   4n  C  n1   8 k  1  3n  n

n

Re G  C * v 2 n

Calculamos el factor de fricción f  f 

16 Re G 16 C * v

2 n

A partir de f y velocidad de circulación, se determina la velocidad 1/ 2

 2 d E f   vm   4fL  



vm  2


2 d E f 16    4*L* 2n  C * v  

2 d E * C * v 

2 n

vm  2

f

4 * L *16 

La función matemática a iterar para determinar velocidad es: f (v)  v 2 

2 d E f * C 2n v 0 4 * L *16 

Calculo de flujo de circulación 3. Problema de transferencia de calor por convección: Se desean calentar 12000 kg/h salsa de tomate desde 18°C hasta 75°C, utilizando un intercambiador de calor de tubos concéntricos. El tomate circula por el interior de un tubo de acero inoxidable AISI 304 de 2 pulgadas estandar, mientras mient ras que por el exterior condensa conden sa vapor de agua saturado a 105 °C. Si se pueden despreciar las resistencias que a la transmisión de calor ofrece la película de condensado y la pared del tubo, calcular la longitud que debe tener el intercambiador para llevar a cabo el citado calentamiento. Problema 14.1 Barboza Canovas

Datos: Propiedades del tomate triturado, en el intervalo de temperatura de operación: calor específico 3,98 KJ/Kg°C, Conductividad térmica 0,5 W/m°C, Densidad 1033 Kg/m3. La viscosidad varía con la temperatura según la expresión   1,75 x10 4 exp(4000/ T ) mPas, en la que T es la temperatura absoluta. Dimensiones del tubo de acero de 2”: diámetro interno 5,25 cm. Diámetro externo 6,03 cm. Supóngase que el coeficiente global de transmisión de calor varía linealmente con la temperatura.




De las Tablas de vapor saturado de agua para tv=105ºC, calor latente de vaporización es:  v105C  2242kJ / kg

El calor de condensación de vapor es cedido al tomate para aumentar su temperatura desde 18ºC hasta 75ºC. wc  v  w f Cp f t s  t e 



wc  2242



kJ 

  1200kg / h 3,98kJ / kg º C 75  18º C

kg 

wc  1214,2kg / h

Según el enunciado del problema se puede despreciar la resistencia que la capa de condensado y la pared ofrecen a la transmisión de calor, por lo que: U=hi Como el coeficiente global de transmisión de calor varia linealmente con la temperatura: U = a+bt, el caudal de calor que atraviesa la sección lateral del tubo metálico será: Q  AU T mtc   dL U T mtc

U T mtc 

U 2 T 1  U 1 T 2

 U 2 T 1    U 1 T 2 

Ln

U 1  h1

T 1  T v  t c  105  18  87º C

U 2  h2

T 2  T v  t s  105  75  30º C

Para el cálculo de los coeficientes individuales de transmisión de calor, se utilizará la ecuación de Sieder-Tate:    0,027k Re 0,8 Pr 0,33   h d   w 

0 ,14




Calculo de la densidad de flujo de la corriente de tomate: G

G

w1  2 d 4 4 *12000kg / h

 * 5,25 *10

m

2 2

1h 2

3600s

 1539,8

kg 2

m .s

La temperatura en la pared metálica coincidirá con la de condensación de vapor, ya que no existe resistencia de la pared metálica y de la capa de condensado: T w=tv=105ºC. Para el cálculo de Re, Pr, h i y h 2, es preciso conocer los valores de la viscosidad a las temperaturas correspondientes. Para ello se utilizará la expresión: 

 1,75 x10 4 exp( 4000 / T )

En la tabla siguiente se hallan recogidos los valores de Re, Pr y μ calculados a partir de las ecuaciones anteriores: Puede observarse que la entrada Re1=496, por lo que para el cálculo de h1 se debe utilizar una expresión para flujo laminar:   d  Nu  1,86RePr    L  

1/ 3

       w 

0,14

Luego,

  1,86k RePr d i 1/ 3    h1  d i   w 

0 ,14

L1/ 3

Para el cálculo de h1 es preciso conocer la longitud del intercambiador, lo que hace que el problema deba resolverse por iteración o tanteo. Al sustituir los valores de las variables, se obtiene: h1  892L1/ 3 W/m 2 º C

Para el cálculo de h 2 se utiliza la ecuación de Sieder.Tate. La longitud del intercambiador se calcula mediante la expresión: L 

Q

 d i U T mtc

Siendo:



Q  wc  v  756,18kW

Ejercicios propuestos 2

d i  5,25*10 m

Sustituyendo en la expresión de la longitud del intercambiador: L 

160,62 * ln(3,91 L1/ 3 ) 3,91  L1/ 3

El cálculo de L se realizará por iteraciones, para la función será: fL   L 

160,62 * ln(3,91 L1/ 3 ) 3,91  L1/ 3

0

Solución por Solución numérica de iteraciones 4. Problema de transferencia de calor: Un fluido alimentario viscoso se halla a 15ºC y se desea aumentar su temperatura hasta 40ºC para introducirlo en un pasteurizador de placas a razón de 1000 kg/h. para realizar este calentamiento se utiliza un intercambiador de tubos concéntricos de pared rasgada, con 4 paletas insertadas en un eje central que gira a 6 rpm. Por la sección anular circular 10000 kg/h de agua caliente, que se introduce a 98 ºC. si el intercambiador se halla perfectamente aislado para evitar pérdidas de calor hacia el exterior, calcular su longitud. Ejercicio 14.8 Barboza C. Datos despreciar el espesor del eje central y de las paletas

Propiedades del fluido alimentario: Conductividad térmica: 0.52 w/m ºC Viscosidad: 1.6 Pa.s Densidad: 1100 kg/m 3 Calor específico: 3.35 KJ/kg ºC Propiedades del agua: Conductividad térmica: 0.58 w/m ºC Viscosidad: 1 mPa.s Densidad: 1000 kg/m 3 Calor específico: 4.185 KJ/kg ºC Los tubos son de acero inoxidable cuya conductividad térmica es de 23 W/m ºC. El tubo interior posee un diámetro interno de 22 cm y un espesor de pared de 8 mm. El tubo exterior posee un diámetro interno de 30 cm. Solución -------------------------------------------------------------------------------------------------------------------------------------------------------80 Alberto L.Huamaní Huamaní



Sección del tubo interior: Si 

 d i

2



4

3,1416* (0,22m) 2 4

 0,0380m 2

Sección anular: S a 



 d 0  d e2 2

4

  3,1416* (0,30m)

2

4

 0,02694m 2

Diámetro equivalente sección anular: 2

 d 0  d e

De  4 * r H 

2

4  d e

d 

2

0

 d e2 

d e

0,30 

2

 0,23 62 

0,23 6

 0,1454m

Densidad de flujo másico para el agua GC 

wc S a



10000kg / h 1h kg  103,1 2 0,02694 3600s m s

El coeficiente de transmisión de calor para el agua se calcula a partir de la ecuación de Dittus-Boelter (ec. 14.14) para fluidos que se enfrían: Nu 

hDe k

 0,023Re 0,8 Pr 0,3

Módulo de Prand  kJ  3  4,185 10 Pa.s  kg º C  Cp     7, 2 Pr  k

0,58 x10 3

Módulo de Reynolds

Re 

 v De 



Gc De 

kJ

s.m.º C

103,1



kg m 2 .s

0,1454m 

3

10 Pa.s

 1,5 x10 4

Al sustituir en la ecuación Dittus-Boelter 0,58 he 

W



mº C 0,023 1,5 x10 4 0,1454m

 7,2 0 ,8

0,3

 363

W m 2 º C

El coeficiente individual de transmisión de calor para el fluido alimentario se calcula a partir de la ecuación: -------------------------------------------------------------------------------------------------------------------------------------------------------81 Alberto L.Huamaní Huamaní



1/ 2

 k  Cp P N  hi  2    

1/ 2

  kg   J   1 min   W   3   0,52  4  6  1100 3   3,35x10  m   kg º C   min 60s   mº C     hi  2         hi  98 8

W 2

m º C

Coeficiente global de transmisión de calor referido al área interna: En la que Dml es: Dml 

1 U i

d e  d i



 d  ln e   d i 

 988

0,236  0,22

1 W 2

m º C

 0,236  ln   0,22 

 0,228m

8 x10 

3

 23

1



W  0,228m 

  mº C  0,22m 

363

W  0,236m  2





m º C  0,22m 

Al operar se obtiene: U i  255

W m 2 º C

Del balance energético se obtiene 10000

kg h

4,185

kJ kg C

98  T S C  3000

kg h

3,35

kJ kg C

40  15 C

T S  92C

De la ecuación de velocidad .

Q  U i Ai Tm L  U i  diLTm L .

Q  w f Cp f t S  t e   3000

kg h

3,35

kJ

40  15C  251250

kg C

kJ h

a) Si los fluidos circulan en contracorriente



Tm Lc 

T e  t S   T S  t e  98  40   92  15   67C  T e  t S    98  40     ln  ln     92 15  T t      S e  251250

Ai 


kJ h

1h / 3600s 

kJ   3  255 *10  * 67C 2 sm C   

 4,085m 2

Siendo la longitud del intercambiador L 

4,085m 2  0,22m

 5,91m

b) Si los fluidos circulan en forma paralelo Tm Lp 

T e  t e   T S  t S  98  15  92  40   66,3C  T e  t E    98  15    ln ln    92 40   T t     S S  251250

Ai 

kJ h

1h / 3600s 

kJ   3  255 *10  * 66,3C 2 sm C   

 4,128m 2

Siendo la longitud del intercambiador L 

4,128m 2  0,22m

 5,97m

SI SE HUBIERA UTILIZADO LA SIGUIENTE EXPRESION Nu 

hd i k

 4,9 Re

0, 57

0 , 47

Pr

 di * N     v 

0,17

 d i     L 

0 ,37

Módulo de Reynolds

Re 

G f d i 



w f d i S i 

3000



kg h

1h / 3600s 0,22 m

0,038 * m  * 1,6 Pa.s  2

3

Módulo de Prand




 kJ   3,35 1,6 Pa.s  º kg C Cp    Pr    10,308 k

kJ

0,52 x10 3

s.m.º C

La velocidad lineal de circulación v

w f

 * S

d i N v

3000kg / s



kg   2 3600s 1100 3 0,038m  m  

0,22m6mi 

 1h

0,02m / s

60s

1



1h

 0,020m / s

 1,104

Al sustituir estos datos en la ecuación anterior se puede obtener el coeficiente hi hi 

0,52W / mC 0,22m

L  968 L0,37

4,93

0 ,57

10,308

0 , 47

1,104

0 ,17

 0,22     L 

0, 37

W / m 2 C

El calculo de L deberá resolverse por iteraciones 1 U i



1 hi

 2,90 4 * 10  3

Como L 

Q

1

 d i  Tm L

U i

L 

 1    2,904 * 10 3   d i Tm L  hi 

L 

1  3   2 , 904 * 10   0 , 37  d i Tm L  968 L 

Q

Q

Contracorriente: L 

Q d i Tm LC



1    2,90 4 * 10 3    0 , 37  96 8 L 

 L0,37   L   2,90 4*10 3   d i Tm LC  96 8  Q




Paralelo: L 

Q d i Tm LP



L 

Q  d i Tm LP

1    2,904 *10 3   0 , 37  968 L 

 L0 , 37    2,904*10 3   968 

2.8.2 Desarrollar el programa que permita hallar el coeficiente de fricción   1 1 2.51  2 log  4 f  3.7 D Re g 4 f  e v

Q A



  ;  

n

2  4n  D n v n Re g     n 1  3 n 1 8 k  

f m

 . A

El comportamiento reológico de una mermelada de albaricoque puede describirse mediante la ecuación de Herschel-Bulkley, presentando los siguientes datos: Flujo másico: fm = 2.2 kg/s Densidad de fluido:  = 1165 kg/m 3 Índice reológico: n = 0.65 Coeficiente de consistencia del fluido: k = 4.43 Pa.s n Rugosidad absoluta de la tubería: e = 0.00005 Diámetro de la tubería: D = 0.0343m Longitud de la tubería: L = 200m Y además que muestre los resultados de velocidad, Reynolds y coeficiente de fricción. Cuando un fluido no Newtoniano pasa por una tubería de sección circular bajo cierto régimen, se genera una pérdida de carga debido a la fricción que puede ser estimada de acuerdo a la ecuación de Darcy:



Lv


2

h f  f D2 g

Donde: hf: pérdida de carga f: Factor de fricción(s/u) L: longitud de la tubería D: diámetro de la tubería v: velocidad promedio del fluido dentro de la tubería g: aceleración de la gravedad Por otro lado, el factor de fricción “f” puede ser estimado mediante la ecuación de Colebrook:  1 2.51   2 log   4 f  3.7( D / e) Re g 4 f  , 1

n

2  4n  D n v n Re g     n 1  3n  1  8 k

Desarrollar un programa en lenguaje Matlab que calcule el caudal “Q”, el área transversal de la tubería “A”, la velocidad de flujo “v”, el número de Reynolds generalizado “Reg”, el factor de fricción “f”, y la pérdida de carga en la tubería “hf”, a partir de los datos de mermelada de albaricoque: Flujo másico: fm = 2,2 kg/s Densidad de fluido:  = 1165 kg/m3 Índice reológico: n = 0,65 Coeficiente de consistencia del fluido: m = 4.43 Pa.sn Rugosidad absoluta de la tubería: e = 0.00005 Diámetro de la tubería: D = 0.0343m Longitud de la tubería: L = 200m 2.8.3 Se desea calcular la velocidad y el caudal másico con el que circula un puré de manzana que se transporta a través de una tubería de 15 cm de diámetro interno, a una temperatura de 25°C. La distancia total que debe recorrer el puré es de 300 metros, existiendo entre los puntos de salida y llegada una caída de presión de 250 Kpa, estando el punto de llegada a 5 metros por encima del punto de salida. A ala temperatura de trabajo, el puré sigue la ley potencia, con un índice de consistencia de 2,4 Pa.s n y un índice de comportamiento al flujo de 0,44, siendo su densidad 1200 kg/m 3.




El número de Reynold crítico correspondiente a este flujo, se obtiene con la ecuación: 6464n

Re g crítico 

 2 n 

1  1 n  1  3n 2     2  n 

Las pérdidas de energía mecánica se obtienen al aplicar la ecuación de Bernoulli entre los puntos de entrada y salida

E f 

P 

 g z

Se supone una velocidad a partir de la cual se evalúa el valor del módulo de Reynold generalizado. n

n n  4n  D v 2 Re g     n 1  3n  1  8 k

El valor del coeficiente de fricción calculamos con la ecuación

f 

16 Re g

. f   2 D E  A partir de f y vm   fL 4  

1/ 2

SOLUCION



Sistema de Ecuaciones lineales




CAPITULO II

SISTEMA DE ECUACIONES LINEALES ---------------------------------------------------------------------------------------------------------Objetivos: Desarrollar ejercicios de ecuaciones lineales aplicados a la ingeniería de alimentos. ---------------------------------------------------------------------------------------------------------2.1 INTRODUCCIÓN Sistema de “m” ecuaciones lineales con “n” incógnitas

Debemos resaltar que muchos problemas prácticos y reales se reducen a un sistema de ecuaciones de este tipo, para ver esto se puede recurrir a cualquier documento de investigación de operaciones. Sistemas lineales son sistemas de ecuaciones con m ecuaciones y n incógnitas formados por ecuaciones lineales. Un sistema lineal con m ecuaciones y n incógnitas es escrito usualmente en la forma:

a11 X 1  a12 X 2  a13 X 3  .........a 1n X n  b1 a X  a X  a X  .........a X  b n 22 2 23 3 2n 2  21 1 a 31 X 1  a 32 X 2  a 33 X 3  .........a 3n X n  b3 .....  a n1 X 1  a n 2 X 2  a n3 X 3  .........a nn X n  bn

(2.1)

Donde: aij : coeficientes 1  i  m, 1  j  n x j : incógnitas j = 1, 2, ..., n bi : Termino independiente i = 1, 2, ..., m

A resolución de un sistema lineal consiste en calcular os valores de x j, j = 1, 2, ..., n, caso en los que existan, que satisfagan las m ecuaciones simultáneamente. -------------------------------------------------------------------------------------------------------------------------------------------------------89 Alberto L.Huamaní Huamaní



Usando notación matricial, el sistema linear pode ser representado por AX = B, donde  a11 a12 ... a1n  a  a22 ... a 2 n 21    . . . .  A    . . .   .  . . . .    a m1 a m 2 ... amn 

 x1   x   2 .  X    .  .     xm 

 b1  b   2 . B    . .   bm 

Este sistema de ecuaciones se simboliza como [A nxn [Xnx1 = [Bnx1, en donde A es la matriz del sistema, X es el vector incógnita y C es el vector de términos independientes, que en forma sintética se simboliza como AX=B. 2.2 CLASIFICACIÓN CUANTO AL NÚMERO DE SOLUCIONES Un sistema lineal puede ser clasificado cuanto al número de soluciones en: Sea Ax=b un sistema lineal de ecuaciones  

      

Cuando un sistema tiene al menos una solución se dice que, el sistema es consistente, en caso contrario es inconsistente. Se dice que el sistema es homogéneo cuando b = 0, es decir el vector b es el vector nulo. Estos sistemas siempre tienen solución, por lo tanto estos sistemas son consistentes. Las solución trivial de estos sistemas es también X = 0 Si b ≠ 0 se dice que el sistema es no homogéneo. Un sistema no homogéneo es consistente (tiene al menos una solución ) si rango( A) = rango(A|b) Los sistemas no homogéneos son consistentes con una solución si rango( A)=n Es consistente con más de una solución si: rango( A) = rango(A⎢b)= rrango(A)

De forma gráfica se puede clasificar los sistemas lineales en Ax=B Si r(A) ≠r(Aa) Inconsistente, no hay solución

Si r(A) = r(Aa) = k, Consistente, existe solución

Si K = n, la solución es única

Si k< n; infinitas soluciones




Cuando todos los termos independientes fueren nulos, esto es, si bi = 0, i = 0, 1, ..., m, o sistema es homogéneo. Todo sistema homogéneo é compatible, pues admitirá por lo menos la solución trivial ( x j = 0, j = 0, 1, 2, ..., n). 2.3 PROPIEDADES DE LAS MATRICES Suma y Resta de Matrices Sólo se pueden efectuar suma y resta de matrices, si tienen las mismas dimensiones. Al sumar o restar dos matrices [A] y [B], el resultado se mostrará en la matriz [C], y se calcula: cij = aij ± bij 5 9  3 2 8 11  D   B  C        1 0   6 2  7 2       Producto de Matrices Para multiplicar una matriz [A] por un escalar g, se multiplica cada elemento de [A] por g.  g * a11 X 1  g * a12 X 2  g * a13 X 3  g * a1n X n  g * b1    B  g *  A   g * a21 X 1  g * a22 X 2  g * a23 X 3  g * a2n X n  g * b2    g * a31 X 1  g * a32 X 2  g * a33 X 3  g * a3n X n  g * b3

Para multiplicar dos matrices [A] y [B], la dimensión de columnas de [A] debe ser igual a la dimensión de filas de [B].

cij 

n

a

b

ik kj

k 1

3 1  22 29 5 9    D   B* C    8 6 *   82 84  7 2 0 4  28 8 

3*5+1*7=22

8*5+6*7=82

Todo sistema de ecuación se puede escribir

 A* x  b -------------------------------------------------------------------------------------------------------------------------------------------------------91 Alberto L.Huamaní Huamaní

(2.3)



La ecuación (2.3) es para cada fila: n

 a x ij

j 1

j

 bi

(2.4)

Despejando la i-ésima incógnita de (2.4) se tiene:



  xi  b a x   i  ij j  aii  j 1  i 0   1 

n

(2.5)

2.4 MÉTODOS DIRECTOS (ALGORITMOS DIRECTOS) Un método es dicho directo cuando la solución exacta x del sistema lineal es obtenida realizando un numero finito de operaciones aritméticas. Son ejemplos conocidos la regla de Cramer, el método de eliminación de Gauss (la triangulación) y el método de Jordan. 2.4.1 Regla de Cramer Sea un sistema lineal con número de ecuaciones igual al número de incógnitas (un sistema n x n), siendo D el determinante de la matriz A, y Dx1, Dx2, Dx3, …..xn los determ inantes de las matrices obtenidas intercambiando en M, respectivamente, la columna de los coeficientes x1, x2, x3, ……xn por la columna de los térmi nos independientes tenemos que: El sistema S será compatible y será solución única si, es solamente si, D  0. En este caso la solución de S es dada por: x1 

Dx1 D

x2 

Dx2 D

x3 

Dx3 D

xn 

Dxn D

La aplicación de Regla de Cramer exige el cálculo de n+1 determinantes ( det A e det Ai, 1  i  n); para n = 20 el número total de operaciones efectuadas será 21 * 20! * 19 multiplicaciones más un número semejante de adiciones. Asimismo, un computador que efectúe cerca de 100 millones de multiplicaciones por segundo llevaría 3 x 10^5 años para efectuar las operaciones necesarias. Como eso, la regla de Cramer es inviable en función del tiempo de computación para sistemas muy grandes. 2.5 MÉTODOS ITERATIVOS Los métodos iterativos o de aproximaciones sucesivas se utilizan cuando los sistemas de ecuaciones a resolver son de grandes dimensiones o bien son sistemas dispersos (la matriz de coeficientes posee muchos ceros). Estos métodos construyen una sucesión de aproximaciones a la solución de las incógnitas hasta obtener una precisión determinada o hasta completar un número determinado de iteraciones. Son ejemplos el método de Jacobi, el de Gauss-Seidel, etc.




Practica 2-1

MÉTODO DE ELIMINACI N DE GAUSS 2.1 MÉTODO DE ELIMINACIÓN DE GAUSS El método de eliminación de Gauss es transformar el sistema lineal original en otro sistema lineal equivalente a los coeficientes de matriz triangular superior, que son de la resolución inmediata. Decimos que dos sistemas lineales son equivalentes cuando tienen la misma solución. El determinante de los sistemas lineales equivalentes es igual. Con (n - 1) pasos el sistema lineal Ax = B se transforma en un sistema triangular equivalente: UX = C, que puede ser resuelto fácilmente por sustituciones. Vamos a calcular la solución AX = B en tres pasos: 1

ª

etapa: Matriz Completa

Consiste en escribir la matriz completa o aumentada del sistema linear original. 2

ª

etapa: Triangulación

Consiste en transformar la matriz A una matriz triangular superior, mediante una secuencia de operaciones elementales en las líneas de la matriz. 3

ª

etapa: Retro-substitución

Consiste en el cálculo de los componentes x1, x2, ..., xn, solución de AX = B, a partir de la solución del último componente ( xn), y en tanto sustituimos regresivamente en las ecuaciones anteriores. Teorema: Sea AX = B un sistema linear. Aplicando sobre las ecuaciones de este sistema una secuencia de operaciones elementales escogidas entre: i) ii) iii)

Cambiar el orden de las dos ecuaciones del sistema; Multiplicar una ecuación del sistema por una constante diferente de cero; Añadir un múltiplo de una ecuación a una otra ecuación;

obtenemos un nuevo sistema

UX = C

y los sistemas AX = B e UX = C son equivalentes.




2.1.1 Resolución de Sistemas Triangulares Sea el sistema linear AX = B, donde A: matriz n x n, triangular superior, con elementos de la diagonal diferentes de cero. Escribiendo las ecuaciones de este sistema, tenemos: a11 X 1  a12 X 2  a13 X 3  .........a1n X n  b1 a 22 X 2  a 23 X 3  .........a 21n X n  b2 a 33 X 3  .........a 3n X n  b3

. . a nn X n  bn

De la última ecuación de este sistema tenemos: xn 

bn a nm

xn-1 puede en tanto obtenido de la penúltima ecuación:

xn1 

bn1  an1,n xn an1,n1

Y asimismo sucesivamente obteniéndose xn-2, ..., x2, e finalmente x1: x1 

b1  a12 x2  a13 x3  ....  a1n xn a11

2.1.2 Estrategias de Pivoteamiento El algoritmo para el método de eliminación de Gauss requiere el cálculo de los multiplicadores: mik  

aik a kk

i = k + 1, ..., n e k = 1, 2, 3, ..., n-1

cada paso del proceso k. Puesto que el coeficiente a kk llamado pivote. ¿Qué pasa si el pivote es nulo? ¿Y si el pivote está cerca de cero? Estos dos casos merecen una atención especial, ya que es imposible trabajar con un pivot e cero. Y trabajar con un pivote cerca de cero puede dar lugar a resultados totalmente inexactos. Esto es porque en ningún cálculo de la calculadora o de ordenador se llevan a cabo con precisión finita y cerca de cero pivotes son más grandes y la altura de la unidad multiplicadora, a su vez, provoca un ensanchamiento de errores de redondeo. -------------------------------------------------------------------------------------------------------------------------------------------------------94 Alberto L.Huamaní Huamaní



Para evitar estos problemas deben adoptar una estrategia de pivoteamiento, o adoptar un proceso de selección de la fila y/o columna pivote. Esta estrategia consiste en: i) ii)

En el inicio da etapa k da fase de escalonamiento, escoger para pivote y elemento de mayor módulo entre los coeficientes: aik , i = k , k + 1, ..., n; Intercambiar las líneas k y i si fuera necesario.

2.1.3 Clasificación del sistema triangular U es un sistema triangular superior escala de m ecuaciones y n incógnitas, tenemos las siguientes posibilidades: i) m = n  sistema compatible y determinado; ii) m < n  sistema compatible y indeterminado.

Si durante o escalonamiento surgir ecuaciones do tipo: 0 x1 + 0 x2 + ... + 0 xn = bm, en tanto: i) Se bm = 0, en tanto eliminaremos la ecuación y continuamos el

escalonamiento; ii) Se bm  0, en tanto concluirse que el sistema es incompatible. Ejemplo: Úsese la Eliminación Gaussiana para resolver:

0,3 x1

 0,1 x 2  7 x 2  0,2 x 2

3 x1

 0,1 x 2

 0, 2 x3

7,0033 x2

 0,2933 x3

3 x1 0,1 x1

 0,2 x3  0,3 x3 10 x3

 7,85  19,3  71,4

10,01 2 x3

 7,85  19,5617  70,08 4

Respuesta x3= 7,00003 x2 = -2,500 x1= 3,000 Desventajas:  División entre cero  Errores de redondeo  Sistemas mal condicionados El procedimiento para resolver este sistema consta de dos pasos: -------------------------------------------------------------------------------------------------------------------------------------------------------95 Alberto L.Huamaní Huamaní



1. Eliminación hacia adelante de incógnitas. 2. Sustitución hacia atrás 2.1.4 Algoritmo de Eliminación de Gauss Paso 1: Para i = 1, ..., n-1 seguir los Pasos 2 a 4. (Eliminación hacia adelante) Paso 2: Sea p el menor entero con i ≤ p ≤ n y a pi≠0. Si p no puede encontrarse, SALIDA (‘No existe solución única’) PARAR Paso 3: Si p≠i intercambiar la fila p por la fila i Paso 4: Para j= i+1, ..., n seguir los Pasos 5 a 6 Paso 5: Hacer mij 

aij aii

Paso 6: Realizar F j-mij Fi e intercambiarla por la fila F j Paso 7: Si ann=0 entonces SALIDA (‘No existe solución única’) PARAR Paso 8: Hacer X n 

C n ann

(Sustitución hacia atrás) i 1 i

b

Paso 9: Para i= n-1, ..., 1 tomar xi 

n

  aiji1 x j i 1 i 1 ij

a

Paso 10: SALIDA X (es decir (X 1, X2, ...,Xn)) PARAR 2.1.5 Procedimiento del programa en Guide Matlab Formulario




Programa function varargout = pushbutton1 _Callback(h, eventdata, handles, varargin) f=str2num(get(handles.edit1,'string')); [m,n] =size(f); salidaa =''; salida=''; i=1; k=2; for s=1:n sal=sprintf('%10.3f ',f(1,s)); salida=[salida,sal]; end while i<=m sal=sprintf('\n'); salida=[salida,sal]; divisor=f(i,i); aux=k; while k<=m j=1; numerador=f(k,i); while j<=n f(k,j)=f(k,j)-(f(i,j)*(numerador/divisor)); if aux==k sal=sprintf('%10.3f ',f(k,j)); salida=[salida,sal]; end j=j+1; end k=k+1; end i=i+1; k=i+1; end set(handles.edit2,'string',salida); % sustitucion hacia atras i=m; j=n; f(i,j)=f(i,j)/f(i,j-1); i=i-1; while i>=1 j=n-1; while j>i f(i,j)=f(i,j)*f(j,n); if f(i,j)>0 f(i,n)=f(i,n)-f(i,j); else f(i,n)=f(i,n)-f(i,j); end -------------------------------------------------------------------------------------------------------------------------------------------------------97 Alberto L.Huamaní Huamaní



j=j-1; end f(i,n)=f(i,n)/f(i,j); i=i-1;

end i=1; while i<=m sal=sprintf('X%d = %10.3f \n',i,f(i,n)); salidaa=[salidaa,sal]; i=i+1; end set(handles.edit3,'string',salidaa); function varargout = pushbutton2 _Callback(h, eventdata, handles, varargin) close Ejecución del programa Ejemplo tenemos sistema matricial o sistema de ecuaciones simples  X 1  X 2  2 X 3  3  3 X 1  X 2  X 3  1 2 X  3 X  4 X  8 2 3  1

Ingresamos los coeficientes del sistema de ecuaciones: 1 1 2 3;3 -1 1 1;2 3 -4 8 los coeficientes de la matriz El resultado es:




Practica 2-2

METODO DE GAUSS-JORDAN 2.2 MÉTODO DE GAUSS – JORDAN Consiste en aplicar operaciones elementares sobre as ecuaciones del sistema linear dado hasta que se abstenga un sistema diagonal equivalente.

LEER n , aij i=0,n,1

divisor = aii j = i , n + 1, 1

aij = aij /divisor

k=1,n,1

i~=k

pivote=a(k,i)

akj = akj - pivote* aij

ESCRIBIR ai,n+1

Fin




2.2.1 Procedimiento del programa en Guide Matlab Formulario

Programa function varargout = pushbutton1 _Callback(h, eventdata, function varargout eventdata, handles, handles, varargin) a=str2num(get(handles.edit1,'string' a=str2num(get(handles.edit1, 'string')); )); [m,n]=size(a); for i=1:m i=1:m divisor=a(i,i); for j=i:n j=i:n a(i,j)=a(i,j)/divisor; end for k=1:m k=1:m if i~=k i~=k pivote=a(k,i); for j=i:n j=i:n a(k,j)=a(k,j)-pivote*a(i,j); end end end end for i=1:m i=1:m x(i)=a(i,n); end x=x'; t=1:m; t=t'; cadena='' cadena= '';; for t=1:m t=1:m cad=sprintf('x%d=%6.2f' cad=sprintf('x%d=%6.2f',t,x(t)); ,t,x(t)); cadena=[cadena;cad]; end set(handles.edit2,'string' set(handles.edit2,'string',cadena) ,cadena) -------------------------------------------------------------------------------------------------------------------------------------------------------100 Alberto L.Huamaní L.Huamaní Huamaní Huamaní



function varargout function varargout = pushbutton2 _Callback(h, eventdata, eventdata, handles, handles, varargin) set(handles.edit1,'string' set(handles.edit1,'string',,''''); ); set(handles.edit2,'string' set(handles.edit2,'string',,''''); ); function varargout = pushbutton3 _Callback(h, eventdata, function varargout eventdata, handles, handles, varargin) close Compilación: Resolver el sistema 4 X 1  2 X 2  5 X 3  4 X 1  3 X 2  2 X 3  2 2 X 1  X 2  7 X 3  1 1 Ingresamos como: 4 2 5 4;1 3 -2 -2;2 -1 7 11







práctica 2-3

METODO DE JACOBI 2.3 METODO JACOBI Este método es muy poco utilizado debido a que el método de Gauss-Seidel converge más rápidamente a la solución y además lo hace cuando no se logra que el método de Jacobi converja. La condición suficiente para que el método de Jacobi converja es que la matriz de coeficientes sea diagonal dominante, es decir que cada elemento de la diagonal principal es mayor en valor absoluto que la suma del resto de los elementos de la misma fila en la que se encuentra el elemento en cuestión. El método de Jacobi consiste en hacer una suposición inicial de la solución e iterar con los valores previos utilizando la ecuación (2.5). La ecuación que describe la iteración de Jacobi es:

xik 1

1  aii

  n  k  b  a x  i  ij j  j 1   i 0  

(2.10)

A continuación, se presenta un algoritmo para este método iterativo. 2.3.1 Algoritmo de Jacobi Considerando la siguiente notación: n: número de ecuaciones aij: elementos de la matriz A (i indica el número de fila y j el número de columna en el que se encuentra el elemento en cuestión) Ci: elementos del vector C X0i: componentes de la primera aproximación al vector solución (esta primera aproximaciónes X0) XiK : componentes de la aproximación de orden k al vector solución, k varía de 1 a N (indica el orden de aproximación o iteración) E: cota de error o criterio de detención N: número máximo de iteraciones Paso 1: Para k = 1 -------------------------------------------------------------------------------------------------------------------------------------------------------103 Alberto L.Huamaní Huamaní



Paso 2: Mientras k ≤ N seguir con los pasos 3 a 6 Paso 3: Para i= 1, ..., n, calcular la aproximación de orden 1 mediante la fórmula:   n  1   X i  C a X  i  ij 0 j  aii  j 1  j i   Paso 4: Si |X − X0| < E, SALIDA X (es decir (X 1,X2, ...,Xn)) PARAR Paso 5: Tomar k = k+1 Paso 6: Para i= 1, ..., n tomar X0i=Xi Paso 7: SALIDA (‘Número máximo de iteraciones completado’) PARAR Algoritmo: Método de Jacobi en Matlab Entradas: matriz A, vector b, aproximaciones iniciales x0, tolerancia TOL, iteraciones máximas MAX. Salidas: valor x tal que Ax=b, residual res, iteraciones it. 2.3.2 Procedimiento del programa GUIDE MATLAB Formulario

function varargout = pushbutton1 _Callback(h, eventdata, handles, varargin) f=str2num(get(handles.edit1,'string')); error=str2num(get(handles.edit4,'string')); -------------------------------------------------------------------------------------------------------------------------------------------------------104 Alberto L.Huamaní Huamaní



err1=100; err2=100; err3=100; [m,n]=size(f); salida=''; i=1; j=1; x1=0; x2=0; x3=0; x1ant=0; x2ant=0; x3ant=0; while err1>error|err2>error|err3>error sal=sprintf('Iteracion %d\n',j); salida=[salida,sal]; %calculo de x1 x1=(f(i,n)-f(i,2)*x2ant-f(i,3)*x3ant)/f(i,i); err1=((x1-x1ant)/x1)*100; if err1<0 err1=err1*-1; end sal=sprintf('x1= %5.3f error = %5.3f\n',x1,err1); salida=[salida,sal]; i=i+1; %calculo de x2 x2=(f(i,n)-f(i,1)*x1ant-f(i,3)*x3ant)/f(i,i); err2=((x2-x2ant)/x2)*100; if err2<0 err2=err2*-1; end sal=sprintf('x2= %5.3f error = %5.3f\n',x2,err2); salida=[salida,sal]; i=i+1; %calculo de x3 x3=(f(i,n)-f(i,1)*x1ant-f(i,2)*x2ant)/f(i,i); err3=((x3-x3ant)/x3)*100; if err3<0 err3=err3*-1; end sal=sprintf('x3= %5.3f error = %5.3f\n\n\n',x3,err3); salida=[salida,sal]; i=1; j=j+1; x1ant=x1; x2ant=x2; x3ant=x3; end set(handles.edit2,'string',salida); sal=sprintf('x1= %5.3f\nx2= %5.3f\nx3= %5.3f',x1,x2,x3); -------------------------------------------------------------------------------------------------------------------------------------------------------105 Alberto L.Huamaní Huamaní



set(handles.edit3,'string',sal); function varargout = pushbutton2 _Callback(h, eventdata, handles, varargin) set(handles.edit1,'string',''); set(handles.edit2,'string',''); function varargout = pushbutton3 _Callback(h, eventdata, handles, varargin) close Resolver la matriz 2 X 1  1 X 2  0,1 X 3  3 0,1 X 1  7 X 2  0,3 X 3  20 3 X 1  2 X 2  10 0 X 3  450

2 -1 0.1 3; 0.1 7 -0.3 20; 3 -2 100 450

y el error de 0.001

Resultado




Practica 2-4

METODO DE GAUSS – SEIDEL 2.4 MÉTODO DE GAUSS – SEIDEL El método de Gauss-Seidel consiste en hacer una suposición inicial de la solución e iterar utilizando los valores previos de las incógnitas que no se han actualizado y los valores de la iteración actual de las incógnitas que ya se han iterado, todo esto utilizando la ecuación (2.5). La ecuación que describe la iteración de Gauss-Seidel es:



 k 1 k 1 k    xi  b a x a x  ij j   i  ij j aii     j 1 j i 1  i0   1 

i 1

n

(2.11)

El método de Gauss-Seidel tiene la ventaja de que converge mucho más rápido que el método de Jacobi, lo cual veremos en el Ejemplo 4.1. En el algoritmo 2.4.1 se muestra el algoritmo del método de Gauss-Seidel implementado en Matlab. 2.4.1 Algoritmo de Gauss-Seidel Considerando la siguiente notación: n: número de ecuaciones aij: elementos de la matriz A (i indica el número de fila y j el número de columna en el que se encuentra el elemento en cuestión) Ci: elementos del vector C X0i: componentes de la primera aproximación al vector solución (esta primera aproximación es X0) Xik : componentes de la aproximación de orden k al vector solución, k varía de 1 a N (indica el orden de aproximación o iteración) E: cota de error o criterio de detención N: número máximo de iteraciones Paso 1: Para k = 1 Paso 2: Mientras k ≤ N seguir con los Pasos 3 a 6 Paso 3: Para i= 1, ..., n, calcular la aproximación de orden i mediante la fórmula: -------------------------------------------------------------------------------------------------------------------------------------------------------107 Alberto L.Huamaní Huamaní



i 1 n  1   C i  aij X j  aij X 0 j  X i   aii  j 1 j i 1  Paso 4: Si X − X0 < E, SALIDA X (es decir (X1,X2, ...,Xn)) PARAR Paso 5: Tomar k = k+1 Paso 6: Para i= 1, ..., n tomar X 0i=Xi Paso 7: SALIDA (‘Número máximo de iteraciones completado’) PARAR





Las condiciones suficientes para que el método de Gauss-Seidel converja es que la matriz de coeficientes sea diagonal dominante o bien que la matriz de coeficientes sea simétrica y definida positiva. Un algoritmo del método se muestra a continuación.




2.4.2 Procedimiento del programa GUIDE Matlab Formulario

Programa function varargout = pushbutton1 _Callback(h, eventdata, handles, varargin) f=str2num(get(handles.edit1,'string')); error=str2num(get(handles.edit4,'string')); err1=100; err2=100; err3=100; [m,n]=size(f); salida=''; i=1; j=1; x1=0; x2=0; x3=0; x1ant=0; x2ant=0; x3ant=0; while err1>error|err2>error|err3>error sal=sprintf('Iteracion %d\n',j); salida=[salida,sal]; %calculo de x1 x1=(f(i,n)-f(i,2)*x2-f(i,3)*x3)/f(i,i); err1=((x1-x1ant)/x1)*100; if err1<0 err1=err1*-1; end -------------------------------------------------------------------------------------------------------------------------------------------------------109 Alberto L.Huamaní Huamaní



x1ant=x1; sal=sprintf('x1= %5.3f error = %5.3f\n',x1,err1); salida=[salida,sal]; i=i+1; %calculo de x2 x2=(f(i,n)-f(i,1)*x1-f(i,3)*x3)/f(i,i); err2=((x2-x2ant)/x2)*100; if err2<0 err2=err2*-1; end x2ant=x2; sal=sprintf('x2= %5.3f error = %5.3f\n',x2,err2); salida=[salida,sal]; i=i+1; %calculo de x3 x3=(f(i,n)-f(i,1)*x1-f(i,2)*x2)/f(i,i); err3=((x3-x3ant)/x3)*100; if err3<0 err3=err3*-1; end x3ant=x3; sal=sprintf('x3= %5.3f error = %5.3f\n\n\n',x3,err3); salida=[salida,sal]; i=1; j=j+1; end set(handles.edit2,'string',salida); sal=sprintf('x1= %5.3f\nx2= %5.3f\nx3= %5.3f',x1,x2,x3); set(handles.edit3,'string',sal); function varargout = pushbutton2 _Callback(h, eventdata, handles, varargin) set(handles.edit1,'string',''); set(handles.edit2,'string',''); function varargout = pushbutton3 _Callback(h, eventdata, handles, varargin) close (gcbf) Resolver la matriz 2 X 1  1 X 2  0,1 X 3  3 0,1 X 1  7 X 2  0,3 X 3  20 3 X 1  2 X 2  10 0 X 3  450

2 -1 0.1 3; 0.1 7 -0.3 20; 3 -2 100 450

y el error de 0.0001

Resultado -------------------------------------------------------------------------------------------------------------------------------------------------------110 Alberto L.Huamaní Huamaní



2 -1 0.1 3; 0.1 7 -0.3 20; 3 -2 100 450 y el error de 0.0001







Practica 2-5

SISTEMA DE ECUACIONES NO LINEALES Método de Newton-Raphson Este método se basa en utilizar el desarrollo de Taylor para aproximar una función derivable en las proximidades de un punto. Partimos de un sistema de la forma f 1  x1 ,......... ., xn   0 f 2  x1 ,......... ., xn   0 ………… f n  x1 ,......... ., xn   0

Del que se pretende obtener una solución. Se supone que la solución buscada, (x1, . . . , xn), se escribe de la forma x1  x10  x1 . ……..

xn  xn0  xn y, por tanto, se tiene que f 1  x10   x1 ,......... ., xn0  xn   0 . ……. f 1 x10   x1 ,......... ., xn0  xn  0

Utilizando el desarrollo de Taylor alrededor de ( x10 .......... , x n0 ) y quedándonos con el primer orden, tenemos




 0   f 1  x   0  n   x  0 f 1  x    i   i 1  xi …….  0   f 1  x   0  n f n  x     xi  0   i1  xi 0



Donde x = ( x10 ,......... ., x n0 ). Introduciendo la notación F = (f 1, . . . , f n) y

 f 1    f 1 .......   x  xn  1   DF   .......... .......   f  f n   n  .......    x1  xn  Queda    

 0    

 0

 

0

F  x   DF  x  x  x   0

Entonces 

0

 0    0  x  x  DF  x  F  x  1

Y el método de Newton consiste en calcular la sucesión k 1

x

 k    k   x  DF  x  F  x      k

1

Dado que calcular explícitamente la matriz inversa de la matriz Jacobiana no es, en general, un proceso muy eficiente desde el punto de vista numérico, a la hora de implementar el método se hace en dos pasos: 1. Se resuelve el sistema  k  k  k 1   DF  x  x   F  x     

2. Se calcula k 1

x

k

k 1

 x   x




Las siguientes funciones muestran una posible implementación en Matlab del método de Newton, para resolver el sistema

1. Programa en archivo.m Crear una carpeta y guardar los siguientes programas function y=f (x) % funcion para utilizar con newtonsi.m y(1)=x(1)^2-10*x(1)+x(2)^2+8; y(2)=x(1)*x(2)^2+x(1)-10*x(2)+8; function df= jac(x) % matriz jacobiana para usar con newtonsi.m df(1,1)=2*x(1)-10; df(1,2)=2*x(2); df(2,1)=x(2)^2+1; df(2,2)=2*x(1)*x(2)-10; function [xr,k]=newtonsi(x,tol,imax) % Metodo de Newton para sistemas de ecuaciones % Uso: [xr,k]=newtonsi(x,tol,imax) % Input: %x = vector x1,x2,...,xn inicial, %tol=tolerancia % Se ha de disponer de las funciones: % f.m funcion y=f(x) donde se define el sistema % jac.m funcion df=jac(x) donde se define la matriz % derivada del sistema. % Output: xr= raiz, k= numero de iteraciones. k=1; epi=1; x1=x; while norm(epi)>tol x=x1; fxn=f (x); axn= jac(x); epi=axn\fxn'; x1=x-epi'; k=k+1; if k>imax disp('no converge') break end end xr=x1; -------------------------------------------------------------------------------------------------------------------------------------------------------115 Alberto L.Huamaní Huamaní



2. Otro programa function [x,fx,xx] = newtons(f,x0,TolX,MaxIter,varargin) % newtons.m para resolver un sistema de ecuaciones de la forma f1(x)=0, f2(x)=0,.. %input: f = 1^st-order vector ftn equivalent to a set of equations % x0 = the initial guess of the solution % TolX = the upper limit of |x(k) - x(k - 1)| % MaxIter = the maximum # of iteration %output: x = the point which the algorithm has reached % fx = f(x(last)) % xx = the history of x h = 1e-4; TolFun = eps; EPS = 1e-6; fx = feval(f,x0,varargin{:}); Nf = length(fx); Nx = length(x0); if Nf ~= Nx, error('Incompatible dimensions of f and x0!' ); end if nargin < 4, MaxIter = 100; end if nargin < 3, TolX = EPS; end xx(1,:) = x0(:).'; %Initialize the solution as the initial row vector %fx0 = norm(fx); %(1) for k = 1: MaxIter dx = - jacob(f,xx(k,:),h,varargin{:})\fx(:);%-[dfdx]ˆ-1*fx %for l = 1: 3 %damping to avoid divergence %(2) %dx = dx/2; %(3) xx(k + 1,:) = xx(k,:) + dx.'; fx = feval(f,xx(k + 1,:),varargin{:}); fxn = norm(fx); % if fxn < fx0, break; end %(4) %end %(5) if fxn < TolFun | norm(dx) < TolX, break ; end %fx0 = fxn; %(6) end x = xx(k + 1,:); if k == MaxIter, fprintf('The best in %d iterations\n',MaxIter), end function g = jacob(f,x,h,varargin) %Jacobian of f(x) if nargin < 3, h = 1e-4; end h2 = 2*h; N = length(x); x = x(:).'; I = eye(N); for n = 1:N g(:,n) = (feval(f,x + I(n,:)*h,varargin{:})-feval(f,x - I(n,:)*h,varargin{:}))'/h2; end function y = f46(x) y(1) = x(1)*x(1) + 4*x(2)*x(2) - 5; y(2) = 2*x(1)*x(1)-2*x(1)-3*x(2) - 2.5; -------------------------------------------------------------------------------------------------------------------------------------------------------116 Alberto L.Huamaní Huamaní



Practica 2-6

APLICACIONES A LA INGENIERIA DE INDUSTRIAS ALIMENTARIAS 1. Ejercicio aplicado. Una formulación de salchicha será hecha de los siguientes ingredientes: Carne de vacuno (magra)- 14% de grasa, 67% de agua, 19% de proteína. Grasa de cerdo - 89% de grasa, 8% de agua, 3% de proteína Proteína aislada de soya - 90% de proteínas, 8% de agua. Se necesita agua para ser añadida (usualmente en forma de hielo) para conseguir la humedad deseada del contenido. La proteína aislada adicionada es el 3% del peso total de la mezcla. ¿Cuánta carne de vacuno magra, grasa de cerdo, agua y soya aislada se necesitará para obtener 100 kg de una formulación teniendo la siguiente composición? : Proteína:15% ; Humedad: 65% ; Grasa: 20% Solución

Carne Vacuno

1

Grasa=14% Proteína=19% Agua=67%

Grasa Cerdo

2

Salchicha


Proteína Soya

5

3



Agua

4

Agua=100%

X1 + X2 + X3 + X4 =100 0,19X1 + 0,03X 2 + 0,90X3 + 0X4 = 15 0,14X1 + 0,89X 2 + 0,02X3 + 0X4 = 20 0,67X1 + 0,08X 2 + 0,08X3 + X4 = 65

Total Proteína Grasa Agua

Como la proteína es 3% del total queda entonces -------------------------------------------------------------------------------------------------------------------------------------------------------117 Alberto L.Huamaní Huamaní


X1 + X2 + 3 + X4 =100 0,19X1 + 0,03X 2 + 2,7 + 0X4 = 15 0,14X1 + 0,89X 2 + 0,06 + 0X 4 = 20 0,67X1 + 0,08X 2 + 0,24+ X4 = 65 X1 + X2 + X4 =97 0,19X1 + 0,03X 2 + 0X4 = 12,3 0,14X1 + 0,89X 2 + 0X4 = 19,94 0,67X1 + 0,08X 2 + X4 = 64,76


Total Proteína Grasa Agua Total Proteína Grasa Agua

Como hay 3 variables entonces se requiere 3 ecuaciones 1 1 1 97;0.19 0.03 0 12.3;0.14 0.89 0 19.94

2. Aplicación en reactores: Se tiene un sistema de tres reactores continuos tipo tanque perfectamente agitado trabajando en serie, en donde se lleva a cabo la reacción A  Pr oductos y se opera isotérmicamente (Fig). Los volúmenes se mantienen constantes y son de 100, 50 y 50 litros, respectivamente. Un Balance de materia en cada reactor, de acuerdo con la ecuación de continuidad, conduce al siguiente sistema de ecuaciones




Calcule la concentración del A en régimen permanente en cada reactor, CA1 = , CA2 = , CA3 = . Si la reacción es de primer orden con respecto a A y la constante de velocidad de reacción es k = 0,1/min, la concentración en moles /Litro. Solucion Realizamos el balance de materia respectivo en el sistema. Entra  Sale  lo.que.reaciona  Acumulado

(1)

Primer Reactor

F .CA0  FR.CA3   ( F  FR)CA1   k 1V 1C An1  dCA1 dt

Segundo reactor

( F  FR)CA1   ( F  FR)CA2   k 2V 2 C An 2 

dCA2 dt

Tercer reactor

( F  FR)CA2   ( F  FR)CA3   k 3V 3C An3  Como no cambia con el tiempo

dCA1 dt

dCA3 dt

, entonces cero (0) el valor para cada caso.

Reemplazando los valores tenemos: Primer Reactor

10(1)  (10)(CA1 )  (10  5)CA1   (0.1)(100)CA1  0 Segundo reactor -------------------------------------------------------------------------------------------------------------------------------------------------------119 Alberto L.Huamaní Huamaní



(10  5)CA1   (10  5)CA2   (0.1)(50)CA2  0 Tercer reactor

(10  5)CA2   (15)CA3   (0.1)50CA3  0 Finalmente queda la matriz  25CA1  0CA2  5CA3  10 15CA1  20CA2  0CA3  0 0CA1  15CA2  20CA3  0

(1) (2) (3)

3. Formulación de helado: Preparar 100 kg de una mezcla de helado con las siguientes características: 11 % de grasa, 12% de proteína, 3% de gelatina y 15 % de azúcar. A partir de los siguientes ingredientes disponibles: M = mantequilla (80% grasa,0,5% SNG), LD = leche descremada (93 % SNG, 0,5% grasa, 2% lactosa), G = gelatina (95% gelatina, 2% azúcar), A = azúcar (94 % sacarosa, 3% gelatina) Tenemos el sistema matricial Mantequilla Leche Gelatina Azúcar Agua X5 Helado X1 descr X2 X3 X4 BM general 1 1 1 1 =100 1 BM Grasa 0,8 +0,005 +0 +0 +0 = 11 BM SNG 0,05 +0,93 +0 +0 +0 = 12 BM gelatina 0 +0 +0,95 +0,03 +0 =3 BM Azúcar 0 +0,02 +0,02 +0,94 +0 = 15 4. Formulación de alimento balanceado: Un granjero desea preparar una ración alimenticia para engordar ganado. Dispone de maíz, desperdicios, alfalfa y cebada, cada una con ciertas unidades de ingredientes nutritivos, de acuerdo con la tabla siguiente. Alimento Componentes

Maíz

Carbohidratos Proteínas Vitaminas Celulosa Costo $/kg

80 28 20 50 18

Desperdicio 15 72 20 10 5

Alfalfa Cebada 35 57 12 20 7

60 25 20 60 20

Requerimiento Diario Unidades/kg 230 180 80 160 --------

1. Desarrollamos los respectivos sistemas de ecuaciones (balance de materia para cada ingrediente nutritivo). 80X1 + 15X2 + 35X3 + 60X4 = 230

Carbohidratos



28X1 + 72X2 + 57X3 + 25X4 = 180 20X1 + 20X2 + 12X3 + 20X4 = 80 50X1 + 10X2 + 20X3 + 60X4 = 160 2. La matriz aumentada será 80 28  20  50

15

35 60 

72 57

25

 20 12 20   10 20 60 

 x1   x2    x3    x4


Proteína Vitaminas Celulosa 230 180 =    80    160  

5. Aplicación en extracción liquido-liquido: En una columna de 5 platos, se requiere absorber benceno contenido en una corriente de gas V, con un aceite L que circula a contracorriente del gas. Considérese que el benceno transferido no altera sustancialmente el número de moles de V y L fluyendo a contracorriente, que la relación de equilibrio está dada por la ley de Henry (y = mx) y que la columna opera a régimen permanente. Calcule la composición del benceno en cada plato. Datos: V = 100 moles/min; L = 500 moles/min y 0= 0.09 fracción molar de benceno en V, x 0= 0.0 fracción molar de benceno en L (el aceite entra por el domo sin benceno). m = 0.12.




Solución Balance de materia para el benceno en cada plato Plato

Balance de benceno

5

L ( x 0  x5 )  V ( y 4  y 5 )  0

4

L( x5  x 4 )  V ( y 3  y 4 )  0

3

L ( x 4  x 3 )  V ( y 2  y 3 )  0

2

L( x3  x 2 )  V ( y1  y 2 )  0

1

L( x 2  x1 )  V ( y 0  y1 )  9

Reemplazando los valores y reordenando se tiene el sistema de ecuaciones: 512 x5 500 x5

 12 x4  512 x4 500 x4

0

 12 x3  512 x3 500 x3

0

 12 x2  512 x2  500 x2

0

 12 x1  512 x1

0 0 0 0 9




Practica 2-6

EJERCICIOS PROPUESTOS 1. Ejercicio de formulación: En una fábrica de embutidos, se elaboran tres tipos diferentes de estos productos que responden a las formulaciones A, B y C. En el cuadro que sigue se presentan las proporciones necesarias de las materias primas principales de cada formulación. Si se dispone diariamente de 140 tn de hígado, 130 tn de carne porcina y 180 tn de tocino, indique la producción diaria de cada tipo de producto en kg/d. Tipo de embutido A B C

Hígado 3 1 3

Carne 1 3 2

Tocino 2 2 4

2. Se tienen un sistema de tres tanques cilíndricos iguales de 6 pies de diámetro, comunicados entre si por medio de tubos de 4 pulgadas de diámetro y dos pies de largo, como se muestra en la Figura. El tercer tanque tiene una salida a través de un tubo de 4 pulgadas de diámetro y 8 pies de largo. Al primer tanque llega un fluido a razón de 0,1 pies cúbicos por minuto e inicialmente su nivel tiene una altura de 20 pies, mientras que el segundo y tercer tanques están vacíos. El fluido es un aceite viscoso cuya densidad e 51,45 lb m/pie3 y cuya viscosidad es 100 centipoises. Calcule la altura del fluido en cada tanque cuando se alcance el régimen permanente. Sugerencia: use la ecuación de Poiseulle para el cálculo de la velocidad media del fluido a través de los tubos. F=0.1 pies3/min

h1

(1)

h2

L=2" D=4"

(2)

h3

L=2" D=4"

(3)

L=8" D=4"

3. Una cadena de supermercados vende carne molida del tipo popular y selecta. Un lote de molida popular contiene 3 kg de grasa y 17 kg de carne roja, un lote de molida selecta contiene 2 kg de grasa y 18 kg de carne roja. Si en un momento dado cuenta -------------------------------------------------------------------------------------------------------------------------------------------------------123 Alberto L.Huamaní Huamaní



con10 kg de grasa y 90 kg de carne roja. ¿Cuántos lotes de molida popular y selecta pueden producir utilizando toda la carne y toda la grasa sin desperdiciar nada? 4. Una granja avícola incluye en la dieta de sus aves vitaminas B, C y D, para evitar enfermedades, así como un desarrollo más rápido. En cierto mes compraron 20 cajas de vitamina B, 40 cajas de vitamina C y 50 cajas de vitamina D pagando $70000, al mes siguiente compraron 30 cajas de vitamina B, 20 de vitamina C y 50 cajas de vitamina D por un total de $51520, un mes después compraron 40 de vitamina B, 10 de vitamina C y 70 de vitamina D con un costo de 45000, sí el precio por caja no ha variado en todo ese tiempo que precio tiene cada caja de vitaminas. 5. Un granjero prepara una mezcla de avena y maíz para alimentar a su ganado. Cada kilogramo de avena contiene 0.15 kilogramos (kg) de proteína y 0.6 kg de carbohidratos, mientras que cada kilogramo de maíz contiene 0.1 kg de proteína y 0.75 kg de carbohidratos. ¿Cuántos kilogramos de cada uno pueden usarse para cumplir con los requerimientos nutricionales de 7,5 kg de proteinas y 50 kg de carbohidratos por comida? 6. Una tienda se especializa en la preparación de mezclas de café para conocedores. El dueño desea preparar bolsas de 1 kilogramo que se venderán a $120 (pesos) usando café colombiano, brasileño y de Kenia. El costo por kilogramo de estos tres tipos de café es de $140, $70 y $100, respectivamente. (a) Muestre que se debe usar al menos 21 kg de café colombiano, y no más de 72 de kg de café brasileño. (b) Determine la cantidad de cada tipo de café suponiendo que el dueño decide usar 81 de kg de café brasileño. 7. Una mezcla de jugo de naranja con 42 % de sólidos solubles es producida, mezclando un concentrado de jugo de naranja de tienda con la reciente cosecha de jugo exprimido. La relación de sólidos solubles/ ácido debe ser igual a 18. El jugo exprimido contiene 12.5% de sólidos solubles, 15.3% de sólidos totales y 0.72% de ácidos. El jugo concentrado de tienda contiene 60% de sólidos solubles, 62% de sólidos totales y 4.3% de ácidos. Calcular : La cantidad de agua que debe ser removida o aumentada para ajustar la concentración de los sólidos solubles para lograr las especificaciones indicadas. Las cantidades del jugo procesado y del concentrado de tienda necesarios para producir 100 kg de mezcla con 42% de sólidos solubles. Por ejemplo, se desea resolver el sistema X 1  X 2  2 X 3  3

3 X 1  X 2  X 3  1 2 X 1  3 X 2  4 X 3  8

Aplicando el método de eliminación gaussiana. a) Eliminación hacia adelante -------------------------------------------------------------------------------------------------------------------------------------------------------124 Alberto L.Huamaní Huamaní

libro-de-metodos-numericos-en-matlab.pdf.pdf

Recommend Documents