UNIVERSIDAD NACIONAL DEL CALLAO

UNIVERSIDAD NACIONAL DEL CALLAO VICE- RECTORADO DE INVESTIGACION INSTITUTO DE INVESTIGACION DE LA FACULTAD DE INGENIERIA QUIMICA INFORME FINAL DEL PROYECTO DE INVESTIGACION TITULADO TEXTO " METODOS NUMERICOS UTILIZANDO POLYMATH, MATHCAD Y MATLAB APLICADOS A INGENIERIA QUIMICA” PRESENTADO POR ING. JUAN MEDINA COLLANA (Del 1 de Marzo del 2010 al 29 de Febrero 2012 Resol. N° 1312-2010 R) 1 ÍNDICE I. RESUMEN 4 II. INTRODUCCIÓN 5 III. MARCO TEÓRICO 8 IV. MATERIALES Y MÉTODOS 8 V. RESULTADOS 9 VI. DISCUSIÓN 9 VII. REFERENCIALES 10 1. ECUACIONES NO LINEALES 1.1 Método de Newton-Raphson 1.2 Método de la secante 1.3 Método de la bisección 1.4 Método de la regla falsa 1.5 Método de iteración del punto fijo 1.6 Problemas de aplicación a la Ingeniería Química 1.7 Aplicaciones en Ingeniería Química 11 11 14 16 19 22 24 32 2. SISTEMAS DE ECUACIONES NO LINEALES 2.1 Método de newton Raphson 2.1.1Aplicación 2.2 Problemas 36 36 37 40 3. INTERPOLACIÓN 3.1 Polinomios de interpolación de newton 3.1.1 Interpolación lineal 3.1.2 Interpolación cuadrática 3.1.3 Diferencias finitas divididas 3.2 Polinomios de interpolación de Lagrange 3.3 Problemas 44 45 46 47 47 48 50 4. REGRESIÓN 4.1 Regresión lineal 4.2 Regresión polinomial 4.3 Regresión lineal múltiple 4.4 Problemas 53 53 54 55 57 5. DIFERENCIACIÓN NUMÉRICA 5.1 Diferenciación mediante método Newton 5.2 Diferenciación de Lagrange: datos discretos 68 67 67 2 5.3 Problemas 69 6. INTEGRACIÓN NUMÉRICA 6.1 Método del trapezoide 6.2 Regla de Simpson 6.2.1 regla de Simpson 1/3 6.2.2Simphson 3/8 6.3 Problemas 76 7. SISTEMA DE ECUACIONES LINEALES 7.1 Métodos de Jacobi 7.2 Métodode Gauss – Seidel 7.3 Problemas 85 85 87 88 8. ECUACIONES DIFERENCIALES NUMÉRICAS 8.1 Método de Euler 8.2 Método de Euler modificado 8.3 Métodos de Runge-Kutta (rk) 8.4 Método de diferencias finitas 8.5 Sistema de ecuaciones diferenciales de primer orden 8.6 Problemas 91 92 94 96 100 101 103 9. ECUACIONES DIFERENCIALES PARCIALES (EDP) 110 9.1 Problemas 110 10. APLICACIONES DE INGENIERÍA QUÍMICA EN 76 77 78 80 81 POLYMATH- MATHCAD Y MATLAB 10.1 Problemas con polymath 10.2 Problemas con mathcad 10.3 Problemas con matlab 112 112 121 131 SILABO 147 3 I. RESUMEN La presente Investigación tuvo como propósito la elaboración de un texto universitario titulado “TEXTO: " METODOS NUMERICOS UTILIZANDO POLYMATH, MATHCAD Y MATLAB APLICADOS A INGENIERIA QUIMICA” Se trata de un texto básico que se expone brevemente los fundamentos teóricos, ilustraciones con problemas resueltos y propuestos de ingeniería química como, equilibrio químico, ecuaciones de estado, transferencia de calor, cinética química y al final se plantea problemas resueltos haciendo uso de software de polymath, mathcad y matlab.. Asi mismo en cada capitulo se prantea problemas propuestos aplicados a la ingeniería química . 4 II. INTRODUCCIÓN. El tema de investigación, está referida a la elaboración de un texto Universitario, cuyo propósito es apoyar en la labor de formación de los alumnos, en el curso de métodos numéricos aplicada a ingeniería .Quimica en la universidad Nacional del Callao facultad de Ingeniería Química. Durante mi experiencia en la docencia, en el intento de encontrar textos necesarios para el dictado del curso de métodos numéricos , se ha podido constatar que existe poca bibliografía en nuestro medio haciendo uso de mathcad, polymath y matlab en un solo texto, mediante ste trabajo se ha desarrrollado problemas aplicados a la ingeniería química haciendo uso de una técnica numérica, al mismo timpo se ha resuelto problemas con el software mathcad, polymath y mathcad, que se vienen usando con mayor intesnsidad en los últimos años a nivel de ingeniería. En este trabajo se ha resaltado el capitulo de ecuaciones no lineales ecuaciones diferenciales puesto y que gran parte de los modelos de ingenieira química se encuentran en esta area. 5 2.1 A. PLANTEAMIENTO DEL PROBLEMA DE INVESTIGACION DESCRIPCIÓN Y ANÁLISIS DEL TEMA El presente trabajo de investigación es una propuesta para la elaboración de un texto Universitario titulado: TEXTO: “METODOS NUMERICOS UTILIZANDO POLYMATH, MATHCAD Y MATLAB APLICADOS A LA INGENIERIA QUIMICA.” Dirigido a estudiantes de pre – grado de Ingeniería Química y otras especialidades afines, que presente de una manera didáctica los principios fundamentales y el uso adecuado de los software de polymath, mathcad y matlab, lo que permitirá cumplir con los propósitos de una adecuada enseñanza y formación profesional. El texto contendrá una base de teoría apropiada y práctica que van a permitir desarrollar criterios y habilidades, que resultara muy valioso para los propósitos de este texto. B. PLANTEAMIENTO DEL PROBLEMA ¿Como elaborar un TEXTO: “METODOS NUMERICOS UTILIZANDO POLYMATH, MATHCAD Y MATLAB APLICADOS A LA INGENIERIA QUIMICA.”, que oriente adecuadamente a los estudiantes de Ingeniería Química? 2.2 2.2.1 OBJETIVO Y ALCANCE DE LA INVESTIGACION OBJETIVO GENERAL Desarrollar un TEXTO: “METODOS NUMERICOS UTILIZANDO POLYMATH, MATHCAD Y MATLAB APLICADOS A LA INGENIERIA QUIMICA.” 6 2.2.2. 1. OBJETIVOS ESPECIFICOS Recopilar Información básica y actualizada, necesaria para iniciar el desarrollo del texto. 2. Analizar y procesar la información para iniciar el desarrollo del texto. 3. Desarrollar los capítulos del texto referido a los fundamentos teóricos. 4. Desarrollar los capítulos del texto referido a las aplicaciones prácticas aplicada la ingeniería química 5. Desarrollar los capítulos del texto referido a las aplicaciones prácticas aplicada la ingeniería química haciendo uso del uso polymath , mathcad y matlab 2.3. ALCANCE DE LA INVESTIGACION El presente trabajo de investigación de acuerdo a la naturaleza del problema se puede manifestar que es una investigación básica y Aplicada, dado los modelos planteados para la solución numérica proviene de fenómenos químicos, cuya solución analítica sería demasiada compleja. El aporte del trabajo de investigación estará orientado al sector académico conformado por los profesores, estudiantes y egresados de la Facultad de ingeniería química de la Universidad Nacional del Callao y otras Universidades del país. Por otro lado, este texto también podría ser utilizado por estudiantes de especialidades afines tales como ingeniería de Alimentos, ingeniería ambiental e ingeniería Industrial. 2.4. IMPORTANCIA Y JUSTIFICACION DE LA INVESTIGACION 2.4.1 IMPORTANCIA Al desarrollar el texto propuesto se facilitará el proceso de enseñanza – aprendizaje en la formación profesional del estudiante universitario a nivel de pre-grado y pueda facilitar los cálculos laboriosos de ingeniería haciendo uso de software. 2.4.2 JUSTIFICACION 7 La contribución del presente trabajo estará orientada a la preparación y el entrenamiento de los alumnos de ingeniería química en el curso de Métodos numéricos, adquiriendo fundamentos teóricos y la parte práctica que consiste en efectuar cálculos de balances de materia, energía, termodinámica, reacciones química, transferencia de calor entre otros. III. MARCO TEÓRICO En la presente investigación se ha incorporado la teoría resumida y simplificada para nueve capítulos del presente texto. Así por ejemplo, en el capítulo I se hace referencia a los diferentes métodos de solución de ecuaciones no lineales con sus respectivos ejemplos y problemas. En el capítulo III , se hace referencia de las técnicas de interpolación , resaltando el método de diferencias por newton . Asimismo en el capítulo VIII de ecuaciones diferenciales , se hace referencia de los diferentes ordenes de solución de Runge Kutta y con sus respectivos ejemplos. En el capítulo X se presenta problemas resueltos con polymath , mathcad y Matlab. IV. MATERIALES Y MÉTODOS Materiales:  Materiales De oficina  Material bibliográfico  Software Polymath  Software Mathcad  Software Matlab  Material de cómputo e impresión Métodos La elaboración del texto, propósito de la investigación le demandó al suscrito, ordenar la información disponible y complicada en función del Syllabus propuesto del curso de métodos numéricos . 8 La estructuración del texto responde a la experiencia docente en la Facultad de Ingeniería Química de la Universidad Nacional del Callao. Para la elaboración del texto, se tuvo cuidado en recurrir a la síntesis de los aspectos teóricos, selección de los problemas de ingeneiria química, elaboración de los programas en matnab, asi como revisión de tutoriales de mathcad. En cuando al planteamiento de problemas, se quiere plasmar la experiencia con el dictado del curso, que me ha permitido revisar una extensa bibliografía sobre la materia, así como volcar en el texto los resultados de mi ejercicio profesional en el campo de Ingeniería Química para satisfacer el propósito de la investigación. V. RESULTADOS El resultado de la presente investigación es la elaboración de un texto universitario titulado: TEXTO: “METODOS NUMERICOS UTILIZANDO POLYMATH, MATHCAD Y MATLAB APLICADOS A LA INGENIERIA QUIMICA.”, , el cual se adjunta al presente. El texto contiene 9 capítulos. La teoría desarrollada en el texto, responde a los aspectos básicos de métodos numericos. Los problemas resueltos en el texto, tienen el propósito de dar las pautas de la aplicación de la teoría desarrollada. Se ha logrado un texto base para el curso de métodos numéricos necesario en la formación universitaria del estudiante de Ingeniería Química. VI. DISCUSIÓN El texto universitario titulado ““METODOS NUMERICOS UTILIZANDO POLYMATH, MATHCAD Y MATLAB APLICADOS A LA INGENIERIA QUIMICA.”, aplicados a la Ingeniería Química, que es el resultado de la investigación a que se refiere el presente informe, se caracteriza por presentar la metodología de calculo de los modelos . Los problemas resueltos y planteados han sido seleccionados con la intención 9 de brindar una mayor claridad a los alumnos y puedan entender los fundamentos teóricos tratados. Los textos de METODOS NUMERICOS contienen demasiada información, muchas veces muy detallado, sin conexión directa con aplicación inmediata. Por eso, el presente texto va a tratar de desenvolver el contenido de tal manera que cada capítulo describa en forma precisa y concreta la teoría y los problemas de aplicación. Sin embargo, y por ende no sustituye el uso de la bibliografía de la especialidad, a la que deberá referirse necesariamente quienes pretendan profundizar en conocimientos de temas específicos. VII. REFERENCIALES Existen textos que se ocupan de los métodos numéricos en general entre los cuales tenemos 1. A. Constantinides and N. Mosotoufi, Numerical Methods for Chemical Engineers with MATLAB Applications Prentice Hall , Upper Saddle River, 1999. 2. Burden, R. Y Faires J. Análisis Numérico. Edit. Iberoamericana, México, 1985. 3. Carnahan, B. Luther, A. Wilkes Cálculo Numérico, Aplicaciones Editorial Rueda, Madrid, 1979 4. Carrasco Venegas Luis . Métodos Numéricos aplicados a la ingeniera 5. Nieves, A., Domínguez, F. Métodos Numéricos Aplicados a la Ingeniería Química Edit. CECSA, México 1985. 6. Nakamura, S. Métodos Numéricos aplicados con software. Edit. Prentice – Hall Hispano Americano, S.A. México, 1992. 7. S.C. Chapra y R. P. Canale,. Métodos Numéricos para Ingenieros. McGraw Hill, México,1999 10 1.ECUACIONES NO LINEALES Los métodos numéricos de resolución de ecuaciones no lineales suelen ser métodos iterativos que producen una sucesión de valores aproximados de la solución, que se espera, que converja a la raíz de la ecuación. Estos métodos van calculando las sucesivas aproximaciones en base a los anteriores, a partir de una o varias aproximaciones iníciales.. Para saber que método debemos aplicar, hay que tener en cuenta la capacidad de separar raíces cercanas, confiabilidad en el alcance de soluciones evitando errores numéricos graves y orden de convergencia. Uno de los problemas que con mayor frecuencia aparece en la ciencia y en la ingeniería es hallar las raíces de una ecuación no lineal de la forma f(x) = 0. Estudiaremos Métodos Iterativos para determinar aproximaciones a raíces reales simples de la ecuación no lineal f(x) = 0. 1.1 MÉTODO DE NEWTON-RAPHSON Este método, es uno de los más usados, a diferencia de los otros métodos, el método de Newton-Raphson no trabaja sobre un intervalo sino que basa su fórmula en un proceso iterativo. Supongamos que tenemos la aproximación xi a la raíz xr de f ( x ) , Figura Nº 1: Demostración método de newton 11 Trazamos la recta tangente a la curva en el punto  xi , f ( xi )  ; ésta cruza al eje x en un punto xi 1 que será nuestra siguiente aproximación a la raíz xr . Para calcular el punto xi 1 , calculamos primero la ecuación de la recta tangente. Sabemos que tiene pendiente m  f   xi  Y por lo tanto la ecuación de la recta tangente es: y  f  xi   f   xi  x  xi  Hacemos y  0 :  f  xi   f   xi  x  xi  Y despejamos x : x  xi  f  xi   x  xi  Que es la fórmula iterativa de Newton-Raphson para calcular la siguiente aproximación: xi 1  xi  f  xi  , si f   xi   0 f   xi  Note que el método de Newton-Raphson no trabaja con intervalos donde nos asegure que encontraremos la raíz, y de hecho no tenemos ninguna garantía de que nos aproximaremos a dicha raíz. Desde luego, existen ejemplos donde este método no converge a la raíz, en cuyo caso se dice que el método diverge. Sin embargo, en los casos donde si converge a la raíz lo hace con una rapidez impresionante, por lo cual es uno de los métodos preferidos por excelencia. También observe que en el caso de que f   xi   0 , el método no se puede aplicar. De hecho, vemos geométricamente que esto significa que la recta tangente es horizontal y por lo tanto no intersecta al eje en ningún punto, a menos que coincida con éste, en cuyo caso xi mismo es una raíz de f ( xi )! Ejemplo Usar el método de Newton-Raphson, para aproximar la raíz de f ( x)  e x  ln x , comenzando con x0  1 y hasta que a  1% . 12 Solución En este caso, tenemos que f ( x )   e  x  1 x De aquí tenemos que: Iniciamos con x0  1 y obtenemos: En este caso, el error aproximado es, Continuamos el proceso hasta reducir el error aproximado hasta donde se pidió. Resumimos los resultados en la siguiente tabla: Tabla Nº1 Aproximación de la raíz y porcentaje de error Aprox. a la raíz Error aprox. 1 1.268941421 21.19% 1.309108403 3.06% 1.309799389 0.052% De lo cual concluimos que la aproximación obtenida es: 13 1.2 MÉTODO DE LA SECANTE Uno de los objetivos de este método es eliminar el problema de la derivada de la función, ya que existen funciones que describen fenómenos físicos y químicos, cuya derivada es muy compleja El método de la secante es muy similar al de Newton con la diferencia principal que en este método de la secante no requiere de la derivada. Este método funciona por medio de dos puntos sobre el eje x, es decir un intervalo (xi-1,xi), los cuales se evalúan en la función para sacar los puntos correspondientes en el eje de la y, los puntos a obtener son f(x i-1) y f(xi), por lo que las coordenadas de los puntos que interceptan a la función son (x i1,f(xi-1)) y el (xi ,f(xi)). Se debe considerar que los puntos xi-1 y xi deben de contener a la raíz, por lo que el punto xi-1 debe estar a la izquierda y el punto xi a la derecha de la raíz. Estos dos puntos que interceptan a la función, se unen por medio de una recta, la cual al cruzar el eje de la x, genera el siguiente punto de acercamiento xi+1 , el cual quedara ubicado entre el intervalo propuesto, como se muestra en la Figura. El método de la secante parte de dos puntos (y no sólo uno como el método de Newton) y estima la tangente (es decir, la pendiente de la recta) por una aproximación de acuerdo con la expresión: Este método se basa en la fórmula de Newton-Raphson, pero evita el cálculo de la derivada usando la siguiente aproximación: f   xi   f  xi1   f  xi  xi1  xi (Recuérdese la solución numérica al problema del cuerpo en caída libre). Sustituyendo en la fórmula de Newton-Raphson, obtenemos: xi1  xi  f  xi  f  xi   xi  f  xi1   f  xi  f   xi  xi1  xi 14 Ejemplo Usar el método de la secante para aproximar la raíz de f  x   e x  x , 2 comenzando con x 0  0 , x 0  1 y hasta que a  1% . Solución Tenemos que f  x0   1 y f  x1   0, 632120558 , que sustituimos en la fórmula de la secante para calcular la aproximación x 2 : Con un error aproximado de: Como todavía no se logra el objetivo, continuamos con el proceso. Resumimos los resultados en la siguiente tabla: Tabla Nº2 Aproximación de la raíz y porcentaje de error Aprox. a la raíz Error aprox. 0 1 100% 0.612699837 63.2% 0.653442133 6.23% 0.652917265 0.08% De lo cual concluimos que la aproximación a la raíz es: x4  0,652917265 15 1.3 MÉTODO DE LA BISECCIÓN Este método es basado en el teorema de Bolzano, que establece que si una función continua cambia de signo en el intervalo (a,b), es decir, f(a)f(b) 0 se llama fórmula de diferencia progresiva, Fórmulas de diferencias divididas hacia adelante Sea los puntos x Xi Xi+1 Xi+2 XI+3 ……. f f (Xi) f (Xi+1) f (Xi+2 ) f (XI+3 ) ……. 5.1.DIFERENCIACIÓN MEDIANTE MÉTODO NEWTON Primera derivada f '  xi   f  xi 1   f  xi  h , f '  xi    f  xi  2   4 f  xi 1   3 f  xi  2h Segunda derivada f ''  xi   f ''  xi   f  xi 2   2 f  xi 1   f  xi  h2  f  xi 3   4 f  xi  2   5 f  xi 1   2 f  xi  h2 5.2. DIFERENCIACIÓN DE LAGRANGE: datos discretos Construimos el polinomio de Lagrange 68 L  x   L1  x  y1  L2  x  y2  L3  x  y3   x  x2  x  x3   x  x1  x  x3   x  x1  x  x2  y1  y2  y  x1  x2  x1  x3   x2  x1  x2  x3   x3  x2  x3  x1  3 Diferenciando la función f  x  L  x   2 x  x2  x3 y  x1  x2  x1  x3  1 2 x  x1  x3 2 x  x1  x2 y2  y  x2  x1  x2  x3   x3  x2  x3  x1  3 Asumiendo espacio constante f  x  2 x  x2  x3 2 x  x1  x3 2 x  x1  x2 y1  y2  y3 2 2 2 x x 2 x 2 f  x  2 x  x2  x3 2 x  x1  x3 2 x  x1  x2 y1  y2  y3 2 2 2 x x 2 x 2 f   x1   2 x1  x2  x3 2x  x  x 3 y1  4 y2  y3 2x  x  x y1  1 1 2 3 y2  1 1 2 2 y3  2 2 x 2 x x 2 x f   x2   2 x2  x2  x3 2 x  x1  x3 y  y1 2 x  x1  x2 y1  2 y2  2 y3  3 2 2 2 2 x 2 x x 2 x f   x3   2 x3  x2  x3 2x  x  x 2x  x  x y  4 y2  3 y3 y1  3 1 2 3 y2  3 1 2 2 y3  1 2 2 x 2 x x 2 x f  x  2 x  x2  x3 2 x  x1  x3 2 x  x1  x2 y1  y2  y3 2 2 2 x x 2 x 2 f   x   y  2 y2  y3 1 2 1 y  y2  2 y3  1 2 1 2 x x x x 2 segunda derivada 5.3 PROBLEMAS 1. El flujo de calor en la interfaz suelo-aire puede calcularse con ley dT q  z  0   k  C dz z 0 69 Z (m) T (º C) 0 13,5 1,25 12 1,75 10 Donde q = flujo de calor, k = coeficiente de difusividad térmica (3.5x10 -7),  = la densidad del suelo (1800Kg/m3), C = calor específico del suelo (840). f '  0   13.5 2  0   1.25  3.75  0  1.25 0  3.75  12 2  0   0  3.75 1.25  01.25  3.75  10 2  0   0  1.25  3.75  0 3.75  1.25  1,333 q = -3.5x10-7 . 1800 (-1,333). =70.56 2. A 400ºC y con concentraciones: CO   NO2   0,10 mol/L se obtuvieron los siguientes datos para la reacción dada: COg   NO2 g   CO2 g   NOg  Tiempo (s) 0 10 20 30 CO  mol L 0,1 0,067 0,05 0,04 Determinar el orden de la reacción dada Determinar la constante de velocidad Calcular la velocidad instantánea de la reacción en el instante t = 10 s Solución: De la ecuación: vA  - dCA  k  CA  dt n  dC  Tenemos: ln  - A   ln x  n ln  CA  ecuación lineal  dt  Diferenciación numérica con 4 puntos: x0 x1 x2 x3 Tiempo (s) 0 10 20 30 CO  mol L 0,1 0,067 0,05 0,04 y0 y1 y2 y3 h  10 70 dCA 1   11y0  18 y1  9 y2  2 y3  dt 6h Para x0:  dCA dt  1 11 0,1  18  0,067  9  0,05  2  0,04  4,4 103 6  10 dCA 1   2 y0  3y1  6 y2  y3  dt 6h Para x1:  dCA dt  1 2  0,1  3  0,067  6  0,05  0,04  2,35 103 6 10 dCA 1   - y0  6 y1  3y2  2 y3  dt 6h Para x2:  dCA dt  1 - 0,1  6  0,067  3  0,05  2  0,04  1,2 103 6  10 dCA 1   2y0  9 y1  18 y2  11y3  dt 6h Para x3:  dCA dt  1 2  0,1  9  0,067  18  0,05  11 0,04  9,5 104 6  10 Llenamos la tabla y procedemos a graficar: TABLA 11 CALCULO DIFERENCIACIÓN CO  mol L - dCA 0 0,1 10    dC A  ln dt   ln CA 4,4 x 10-3 -5,426 -2,302 0,067 2,35 x 10-3 -6,053 -2,703 20 0,05 1,2 x 10-3 -6,725 -2,995 30 0,04 9,5 x 10-4 -6,959 -3,218 Tiempo (s) dt 71 Del ajuste lineal tenemos que la ecuación es: y  1, 7 x  1, 41 Donde: Orden de reacción es: n  pendiente  1, 7  2 Constante de velocidad es: ln K  -1,41 Para t=10 s: v  k  C At 10 s   K  0,244 n v  0,244   0,067 2  1,095  103 PROBLEMAS 1. La primera ley de difusión de Fick establece que el flujo másico Flujo másico = F  D dc dx Donde flujo másico es cantidad de masa que pasa a través de una unidad de area por unidad de tiempo ( g/ cm 2/s. D= coeficiente de difusión ( cm2/s). C: concentración X: distancia en cm. Un ing. químico mide la concentración de un contaminante en los sedimentos depositados en un lago obteniéndose los siguientes datos. X (cm) 0 1 2 c. 10-6 g/cm3 0,1 0,4 0,9 72 Use la mejor técnica de derivación numérica para estimar en x =1,5 y calcule el flujo másico del contaminante cuando D= 2x10 -6cm2/s, para un lago de 3x106m2 de sedimentos ¿cuánto contaminante podría ser transportada hacia el lago en un año? 2. Para la descomposición del pentaóxido de dinitrógeno a 45ºC se obtuvieron los siguientes datos: 2N 2 O5 ( g )  O2 ( g )  4NO2 ( g ) Tiempo(s) 0  N 2 O5  200 400 600 800 mol 2,5 2,22 1,96 1,73 1,53 L determina de qué orden será la reacción 3. A 440ºC con y con concentraciones |CO| = |NO2| = 0,10 mol/L se obtuvieron los siguientes datos para la reacción dada: CO( g )  NO2 ( g )  CO2 ( g )  NO( g ) Tiempo(s) CO  mol L 0 10 20 30 40 100 1000 0,1 0,067 0,05 0,04 0,033 0,017 0,002 a) Determina el orden de la reacción dada b) Determine la constante de velocidad c) Calcula la velocidad instantánea de la reacción en el instante t= 10 s. 4. Se ha estudiado la descomposición térmica de arsenamina sobre vidrio, 2 AsH3(g) 2 As(s) + 3H2(g) La presión total del sistema varia con el tiempo, a 350C, según se indica a continuación: t/h 0 4.33 16.00 25.50 37.66 44.75 Pt/cmHg 39.2 40.30 43.65 45.35 48.05 48.80 a) Determine el orden de reacción con respecto a la arsenamina. b) Calcule la constante de velocidad. 73 5. Para poder estudiar el comportamiento físico químico de una solución de bromo en presencia de luz solar, se coloca una pequeña cantidad de solución de bromo en un matraz y se expone al sol, obteniéndose los siguientes resultados: Tiempo( min.) 10 20 30 40 50 60 Ppm Br2 a 25ºC 2,45 1,74 1,23 0,88 0,62 0,44 En base a esta información calcular la constante de velocidad en las unidades adecuadas y orden de la reacción. 6. En el estudio experimental de la velocidad de reacción para la reacción 2 NO (g) + 2 H2(g) → N2(g) + 2 H2O (g) se han obtenido los siguientes datos a 298 K: [NO] mol/L 0,02 0,02 0,03 0,06 Calcular: [H2] mol/L 0,15 0,45 0,50 0,50 Velocidad mol/l s 1,2 x 10-8 3,6 x 10-8 9,0 x 10-8 36,0 x 10-8 a) La ecuación de velocidad para esta reacción. b) El orden de reacción c) La constante de velocidad 7. Determinar el orden de reacción para la descomposición en fase gaseosa del peróxido de di-tert-butilo (CH3)3COOC(CH3)  C2H6 + 2CH3COCH3 La reacción se lleva a cabo en el laboratorio en un reactor isotérmico batch La presión total de reacción a lo largo de reacción evoluciona según la siguiente tabla 74 tiempo(min) Presión total (MmHg) 0.0 7.5 2.5 10.5 5.0 12.5 10.0 15.8 15.0 17.9 20.0 19.4 8. Determinar el orden de reacción: CH3-Cl(g) + H2O(g)  CH3-OH(g) + HCl(g) usando los datos de la tabla. v  k  [CH3 -Cl]n  [H2 O]m Experiencia [CH3-Cl] (mol/l) [H2O] (mol/l) v (mol·l–1·s–1) 1 0,25 0,25 2,83 2 0,50 0,25 5,67 3 0,25 0,5 11,35 9. La transferencia de calor está dada por la siguiente ecuación dT q   kA dy donde J   k = conductividad térmica   smK    A = área m2 T = temperatura K  y = distancia m  k  0.025 J smK A  3 m2 La temperatura en función de distancia T  1493y3  2200y 2  1076y  500 La transferencia de calor para y= 1,2m 75 6. INTEGRACIÓN NUMÉRICA En esta lección comenzamos el estudio de métodos numéricos para el cálculo numérico de integrales de la forma.  f   b  f  x  dx a Un método común para aproximar I(f) es reemplazando f(x) con un polinomio de interpolación. Este procedimiento se conoce como las reglas de Cuadratura de Newton. Examinamos los primeros dos casos de este método donde se usan polinomios de interpolación lineales y cuadráticos. 6.1 Método del trapezoide: Sea p1(x) el polinomio lineal que interpola a f(x) en x=a y x=b f(x) L(x) x0 x1 x Figura Nº9: Integración método de trapecio usando dos puntos Usando la fórmula para el área de un trapezoide o integrando p1(x) directamente se obtiene que  b  b a a f ( x)dx  h  f ( x0 )  f ( x1 ) 2 f ( x ) dx  h  f (a )  f (b) 2 Generalizando Generalizando el método de trapecio 76 f(x ) h  b a n x0 h x1 h x2 h x3 h x4 x Figura Nº10: Integración método de trapecio usando n puntos  b a x1 x2 xn x0 x1 x n 1 f(x)dx   f(x)dx   f(x)dx     f(x)dx h h h  f(x 0 )  f(x1 )   f(x1 )  f(x 2 )     f(x n 1 )  f(x n ) 2 2 2 h   f(x 0 )  2f(x1 )    2f(x i )    2 f ( xn 1 )  f ( xn )  2 b h a f(x)dx  2 f(x 0 )  2f(x1 )    2f(x i )    2 f ( xn1 )  f ( xn )  4 Evaluar la integral :  xe 2 x dx método de trapecio 0 4 I   xe 2 x dx  0 40  f (0)  f (4)  2(0  4e8 )  23847.66 2 6.2 REGLA DE SIMPSON Además de aplicar la regla trapezoidal con segmentos cada vez más finos, otra manera de obtener una estimación más exacta de una integral, es la de usar polinomios de orden superior para conectar los puntos. Por ejemplo, si hay un punto medio extra entre f(a) y f(b), entonces los tres puntos se pueden conectar con un polinomio de tercer orden. A las fórmulas resultantes de calcular la integral bajo estos polinomios se les llaman Reglas de Simpson. 77 6.2.1 REGLA DE SIMPSON 1/3 La Regla de Simpson de 1/3 proporciona una aproximación más precisa, ya que consiste en conectar grupos sucesivos de tres puntos sobre la curva mediante parábolas de segundo grado, y sumar las áreas bajo las parábolas para obtener el área aproximada bajo la curva. Por ejemplo, el área contenida en dos fajas, bajo la curva f(X) en la Figura. 2, se aproxima mediante el área sombreada bajo una parábola que pasa por los tres puntos: f(x) x0 h L(x) x1 h x2 x Figura Nº11: Integración método de Simpson usando tres puntos  b a f ( x)dx  2 c i0  i f ( x i )  c 0 f ( x 0 )  c1 f ( x1 )  c 2 f ( x 2 ) h  f ( x 0 )  4 f ( x1 )  f ( x 2 )  3 Generalizando 78 h ba n f(x) …... x0 h x1 h x2 h x3 h x4 xn-2 xn-1 xn x Figura Nº12: Integración método de Simpson usando n puntos  b a x2 x4 xn x0 x2 xn 2 f(x)dx   f(x)dx   f(x)dx     f(x)dx h h f(x 0 )  4f(x1 )  f(x 2 )  f(x 2 )  4f(x 3 )  f(x 4 ) 3 3 h     f(x n  2 )  4f(x n 1 )  f(x n )  3 h   f(x 0 )  4f(x1 )  2f(x 2 )  4f(x 3 )  2f(x 4 )   3 4f(x 2i-1 )  2 f ( x2i )  4f(x 2i 1 )    2 f ( xn  2 )  4 f ( xn 1 )  f ( xn )  Evaluar la integral  4 0 xe 2 x dx 4 I   xe 2 x dx  0  h  f (0)  4 f (2)  f (4) 3 2  0  4(2e 4 )  4e8   8240.411 3 Simpson’s 3/8 Approximate by a cubic polynomial La derivación de la Regla de los Tres Octavos de Simpson es similar a la regla de un tercio, excepto que se determina el área bajo una parábola de tercer grado que conecta 4 puntos sobre una curva dada. La forma general de la parábola de tercer grado es: 79 L(x) x0 h f(x) x1 h x2 h x3 x Figura Nº13: Integración método de Simpson usando cuatro puntos  b a 3 f ( x) dx   ci f ( xi )  c 0 f(x 0 )  c1f(x1 )  c 2 f(x 2 )  c3 f(x 3 ) i 0  h 3h  f ( x0 )  3 f ( x1 )  3 f ( x2 )  f ( x3 ) 8 b-a 3 6.2.2 SIMPHSON 3/8 4 I   xe 2 x dx  0 3h  4 8  f (0)  3 f ( )  3 f ( )  f (4)  8  3 3 3(4/3) 0  3(19.18922)  3(552.33933)  11923.832  6819.209 8 5216.926  6819.209   30.71% 5216.926  80 6.3 PROBLEMAS 1. Los datos que aparecen en la tabla siguiente corresponden a una velocidad transversal de salida en una tubería ¿cuál es el flujo volumétrico en m3/s? r (cm) 0 5 10 15 V(m/s) 50 49,5 49 48 20 25 30 35 40 45 47 50 46,5 45 43 40,5 37,5 34 25 0 Q  2  V * rdr 2. La cantidad de calor necesaria para elevar la temperatura de 1000gr de H2O desde –1000C hasta 2000C se puede evaluar de la siguiente forma:  H   m T 2 T 1 Cp  T  dT Cp  2.64 107 T 2  1.56 104 T  0.132 Determinar  empleando la técnica del trapecio y la regla de Simpson 3. La masa que entra o sale de un reactor perfectamente agitado en un periodo especifico se puede determinar mediante: t2 M   Q.cdt ti Donde t1 y t2 son el tiempo inicial y final respectivamente Q es el flujo volumétrico y c es la concentración. 81 4. T(,min ) C,mg/m3 0 10 5 22 10 35 15 47 20 55 25 58 30 52 35 40 40 37 45 32 50 34 Los datos que se muestra a continuación corresponde a la reacción de tipo. A  B en un reactor tubular X (-1/rA) 0 0.0053 0.1 0.0052 0.2 0.005 0.3 0.0045 0.4 0.0033 0.5 0.0025 0.7 0.0018 0.8 0.00125 0.85 0.001 La ecuación de diseño para este tipo de reactor es de: V  FA0   x 0 1 dx  rA V:volumen(m3)FA0: flujo de alimentación A ,X: conversión el flujo molar de alimentación es de 2 mol/ min. el 80 % es de A y la diferencia de inertes Determine el volumen del reactor para una conversión de 90% de A 82 5. Se tiene los siguientes puntos con h constante x X0 X1 X2 X3 X4 y Y0 Y1 Y2 Y3 Y4 Demostrar una fórmula de integración 6. La reacción de saponificación de acetato de etilo en medio alcalino se conduce en fase homogénea de acuerdo con la siguiente ecuación. CH3COOCH2CH3 + NaOH  CH3COONa + CH3CH2OH Donde la velocidad de reacción es de ra  5817 dCA  kCACB , k  2, 028 *109 * e T dt de reacción Constante de velocidad el tiempo de reacción para un reactor discontinuo es calculado mediante la siguiente reacción. CA t 1 dCA ( ra ) CA0  Donde : CA = CA0(1-XA) ; CB = CB0 – CA0XA Si la concentración inicial de A y B es de 0,8 mol-g/ L ¿ Cual es el tiempo para alcanzar una conversión de 95%? A 500K 7. La cantidad de calor necesaria para elevar la temperatura de 1 000 g de H2O desde 100C hasta 900C se puede evaluar de la siguiente forma: H  m  Cp T dT T2 T1 Donde: Cp  2.64 107 T 2  1.56 104 T  0.132 Determinar  empleando la técnica del trapecio y la regla de Simpson 8. Resolver la siguiente integración. 83 9. En un reactor tubular continuo se realiza la eliminación de un contaminante en fase gas mediante la reacción: A + B  R + S, a 200ºC y 7 atm. Para ello, se alimenta una corriente con A y B en relación equimolar, en las mismas condiciones de presión y temperatura, con un caudal de 500 L/min. Si se quiere alcanzar una conversión del 85%, ¿qué volumen de reactor es necesario? Datos: (rA) = 10.CA.CB [mol/L.min] 10. Una reacción A  R, de ecuación cinética r (mol/L.h) = 1,5 CA (a 50ºC), se lleva a cabo en un reactor discontinuo, de forma isoterma, siendo la temperatura de trabajo constante e igual a 50ºC. La mezcla reaccionante tiene una concentración inicial de 10 mol/L en A ¿Cuál es el tiempo de reacción necesario para alcanzar una XA de 0,90? 84 7. SISTEMA DE ECUACIONES LINEALES INTRODUCCIÓN La solución de sistemas de ecuaciones lineales es un tema de gran utilidad en diversas ramas del conocimiento como la economía, la biología, física, ingenierías, etc. La resolución de sistemas de cualquier número de ecuaciones (10, 50, 100, 500, etc.) es una realidad hoy en día, gracias a los computadores, lo cual proporciona solución directa. Un gran número de problemas prácticos de ingeniería se reduce a resolver un sistema de ecuaciones. Un sistema de ecuaciones lineales con incógnitas, tiene la forma: a11x1 + a12x2 + ….. + a1nxn = b1 a21x1 + a22x2 + ….. + a2nxn = b2  am1x1 + am2x2+ …. + amnxn = bm Con notación matricial se escribe así: a11 a  21   a m1 a12 ..... a1n  a 22 ......a 2 n    a m 2 .....a mn   x1  b1  x     2   b2           xm  bm   AX=B Donde A es la matriz de coeficientes del sistema X es el vector incógnita B es el vector de términos independientes. 7.1 METODOS DE JACOBI (Método de desplazamiento simultáneos) a11 x1 + a12 x2 + a13 x3 = b1 a21 x1 + a22 x2 + a23 x3 = b2 a31 x1 + a32 x2 + a33 x3 = b3 85 cona11, a22 y a33 distintos de cero. Se despeja x1de la primera ecuación, de la segunda y de la tercera con lo que se obtiene. =− =− que en notación matricial queda = + + − + − =− ⎡ 0 ⎢ = ⎢− ⎢ ⎢− ⎣ − − − + ⎡ ⎢ ⎢ +⎢ ⎢ ⎢ ⎣ − ⎤ ⎥ ⎥ − ⎥ 0 ⎥ ⎦ 0 ⎡ 0 ⎢ = ⎢− ⎢ ⎢− ⎣ − − 0 − ⎤ ⎥ ⎥ − ⎥ 0 ⎥ ⎦ ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ =⎢ ⎢ ⎢ ⎣ ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ Una vez que se tiene la forma en notación matricial, se propone un vector inicial ( ) que puede ser solución . ( ) = 0, o algún otro que sea aproximado al vector Para iterar existe dos variantes. Si ( ) = es el vector aproximación a la solución después de tiene para la siguiente aproximación. 86 iteraciones, entonces se ( ) 1 ( ⎡ ⎢ ⎢ 1 ( =⎢ ⎢ ⎢ 1 ( ⎣ = − − − − − )⎤ ⎥ ⎥ )⎥ ⎥ )⎥ ⎦ − 7.2 METODO DE GAUSS - SEIDEL (método de desplazamientos sucesivos) En este método los valores que se van calculando en la ( + 1) −ésima iteración se emplean para calcular los valores faltantes de esa misma iteración; es decir, con ( ) se calcula ( ) de acuerdo con: ( ⎡ ⎢ 1 ( =⎢ ⎢ 1 ⎢ ( ⎣ − − − − − 87 − ) ⎤ )⎥ ⎥ ⎥ )⎥ ⎦ 7.3 PROBLEMAS 1) Determine los valores de x1 , x2 y x3 usando el método de gaus seidel 2 x1  7 x2  11x3  6 x1  2 x2  x3   5 7 x1  5x2  2 x3  17 2) Se tiene la siguiente torre de separación realizar los balances y calcular los flujos masicos de cada corriente. H2O 0,2% de HCl 99.8% aire Masa = 2640Kg/h 43,1% de HCl 63,8% H2O 56,9% aire 36,2% de HCl Figura Nº14:Diagrama de una columna de absorción 3) Un proceso de 5 etapas en equilibrio para una extracción liquido o absorción gaseosa puede ser modelado mediante un sistema de ecuaciones lineales de acuerdo como se muestra a continuación. (1  rr ) X 1  rrX 2   F0 X i1  (1  rr) X i  rrXi1  0 para i = 2,3, ………..(n-1) X n 1  (1  rr) X n  F Para : rr = 0,9 ; F0 = 0,05 y F = 0,5 Resolver el sistema de ecuaciones lineales 88 4) Para el sgt.sistema de separación, se conoce el flujo masico(Kg/hr) y las fracciones de masa de cada sustancia de entrada de la columna de separacion1 (F0) y la salida de la columna 2(F3 ,F1Y F4)¿ Cual es el valor de F2  F1 W 1 = 0.2 W 2 = 0.6 W 3 =0.2 W 1 = 0.2 W 2 = 0.6 W 3 =0.2 F0 W 1 = 0.2 W 2 = 0.6 W 3 =0.2 F3 F2 F4 W 1 = 0.2 W 2 = 0.6 W 3 =0.2 5) Resolver la siguiente sistema de ecuación 4 x1  x 2   1  8 x1  x 2  x 3  1 3 3 x2  2 x3  4 x4   3 x 3  x 4  x 5  2 ,1 2 x 4  6 x 5  3, 4 6) La resistencia de un termistor varia con la temperatura 1 2 3  a 0  a1 ln( R )  a 2 ln R   a 3 ln R  T T es la temperatura en Kelvin, R es resistencia en ohms, y a0 , a1 , a2 , a3 son constantes R T ohm C 1101.0 25.113 911.3 30.131 636.0 40.120 451.1 50.128 Determine las constantes 89 7) En la figura se muestra la disposición de 5 tanques donde se lleva a cabo una mezcla de corrientes de entrada y salida. En base a la información determine concentración en cada en estado estacionario 2m3/min. C5 3 3m /min. 2m3/min. 5m3/min. 10mg/m3 C1 C1 C2 C4 3 3m /min. 3 1m /min. 1m3/min. C3 1m3/min. 8m3/min. 8m3/min.20mg/m3 Figura Nº15: Sistema de reactores 90 11m3/min. 8. ECUACIONES DIFERENCIALES NUMÉRICAS INTRODUCCIÓN Con mucha frecuencia aparecen problemas en ingeniería, física, química, ecología meteorología, sociología, etc., que exigen el manejo de ecuaciones diferenciales, muchas de las cuales no se pueden resolver pos los métodos convencionales, teniéndose que recurrir a métodos numéricos aproximados. Se llama ecuación diferencial aquella ecuación que contiene una variable dependiente y sus derivadas con respecto a una o más variables independientes. Se dividen en dos grandes grupos: Ordinarias, si contienen una sola variable independiente y Parciales, cuando contienen varias variables independientes. Estudiaremos únicamente las ecuaciones diferenciales ordinarias (EDO). Estas se clasifican y estudian según el orden de la mayor derivada que aparece en la respectiva ecuación diferencial. Ejemplo: dy  ky dt m d2 y  ky dt 2 es de 1er orden es de 2º orden d3 y dy  5  6y  0 es de 3er orden 3 dx dx Los problemas que encierran el uso de ecuaciones diferenciales ordinarias, constarán de: Una ecuación diferencial ordinaria: dy  y1  f(x, y) dx Las condiciones iniciales. El intervalo para la variable independiente x, para el cual debe determinarse la función y, si existe. 91 8.1. MÉTODO DE EULER Consiste en dividir el intervalo de X0 a Xn en n subintervalos de ancho h. h xn - xo obteniéndose un conjunto discreto de (n + 1) punto: x0, x1, x2 …… n xn en el intervalo [x0 , xn] generándose la sucesión de aproximaciones siguientes: y1 = y0 + hf (x0 , y0) y2 = y1 + hf (x1 , y1) y3 = y2 + hf (x2 , y2) . yn+1 = yn + hf(xn , yn) con i = 0 hasta n siendo f(xi , yi) la ecuación diferencial evaluada en xi y yi. yi+1 = yi + hf(xi, yi) Ejemplo 1: Utilice el método de Euler para integrar numéricamente la ecuación: y’=f(x,y) = -2x3 + 12x2 – 20x + 8.5 de x=1.0 hasta x=3.0 considerando 4 intervalos y el valor inicial en x=1 es y=3 x0 = 1 h f(x0 , y0) = f(1 , 3) 1.0 1.5 x0 x1 2.0 x2 2.5 x3 3.0 x4x4 3 1  0.5 4 yi+1 = yi + hf(xi , yi) i xi yi 0 1 3 F(xi , yi) yi+1=yi+hf(xi ,yi) -1.5 2.25 1 1.5 2.25 -1.25 1.625 2 2.0 1.625 0.50 1.875 3 2.5 1.875 2.25 3.000 4 3.0 3.000 2.5 4.25 92 Ejemplo 2. Dada la ecuación diferencial dy  f(x, y) = x - y , resuelva por el dx método de Euler para el intervalo [0, 1] considerando 5 intervalos y siendo y(0)=2 1 0 = 0.2 5 h Y(1)=? yi+1 = yi + hf(xi , yi) i 0 xi 0 x0 0 x1 x2 x3 0.2 0.4 0.6 yi x4 x5 0.8 00.800.80 F(xi , yi) yi+1=yi+hf(xi ,yi) 2 -2 1.6 1 0.2 1.6 -1.4 1.32 2 0.4 1.32 -0.92 1.136 3 0.6 1.136 -0.536 1.0288 4 0.8 1.0288 -0.2288 0.98304 5 1.0 0.9834 -0.01696 0.986432 El valor exacto es 1.10364 E R  1.10364  0.986432  100%  10.93% 1.10364 APLICACIÓN Un balón de acero de 1200 K es enfriado con corriente de aire a temperatura ambiente de 300K. Asumiendo la transferencia de calor en forma de radiación es mediante la ecuación diferencial . dT  2,2067  10 12 T 4  81 108 , T 0   1200 K dt Determine la temperatura para t= 480s usando el método de euler dT  2, 2067  1012 T 4  81  108  dt f  t , T    2, 2067  10 12  T 4  81  108  T1  T0  f  t0 , T0  h Ti 1  Ti  f  ti , Ti  h  1200  f  0,1200 240    1200   2.2067  10 12 1200 4  81  10 8  240 93  1200   4.5579 240 T1  106.09K  0  240 t  t1  t0  h  240 T1  T  240  106,09K 2  1  f t1 ,1  h  106.09  f  240,106.09 240    106.09  2.2067  10 12 106.09 4  81  108  240  106.09   0.017595 240  110.32K t  t2  t1  h  240  240  480 2    480  110.32K 8.2. MÉTODO DE EULER MODIFICADO Una fuente fundamental de error en el método de Euler es que la derivada al principio del intervalo se supone que se aplica a través del intervalo entero y para obtener una exactitud razonable se utiliza un intervalo muy pequeño a cambio de un error de redondeo mayor. El método de Euler modificado trata de evitar este problema utilizando un valor promedio de la derivada tomada en los dos extremos del intervalo en lugar de la derivada tomada en un solo extremo. El método de Euler modificado consta de dos pasos básicos: 1. Se parte de (x0,y0) y se utiliza el método de Euler a fin de calcular el valor de y correspondiente a x1. Este valor de y se denotará como y1 ya que solo es un valor transitorio para y1. Esta parte del proceso se conoce como paso predictor. 94 2. El segundo paso se llama corrector, pues trata de corregir la predicción. En el nuevo punto obtenido  x1 , y1  se evalúa la derivada de f  x1 , y1  usando la ecuación ordinaria que se está resolviendo. Se obtiene la media aritmética de esta derivada y la derivada en el punto inicial (x0,y0). 1  f(x 0 , y 0 )  f(x1 , y1 )   derivada promedio 2 se usa la derivada promedio para calcular un nuevo valor de yi con la ecuación de Euler yn+1 = yi + hf(xi, yi) que deberá ser más exacto que y1.  y1  y 0  h  f(x 0 , y 0 )  (f(x1 , y1 ))   2 que se tomará como valor definitivo de y1. xn  x0   h    n  Este procedimiento se repite hasta llegar a yn. El esquema iterativo para este método quedaría en general así: Usando el paso de predicción resulta: yi 1  yi  hf(xi ,yi ) Se calcula la derivada f ( x i  1 , y i  1 ) Se establece la derivada promedio (llamémosla B) B 1 f(x i , y i )  f(x i  1 , y i  1 ) 2   Se sustituye f(xi,yi) con este valor promedio en la ecuación de iteración de Euler y se obtiene y i 1  y i  h  f(x i , y i )  f(x i 1 , y i 1 )   2 o simplificado: yi+1= yi+ hB Ejemplo 1: Resuelva f(x,y) = x - y en el intervalo [0,1] Con y(0)=2 x0=0 y0=2 xn=1 yn=? y(1)=? h  n=5 y Con 5 intervalos. x n  x n 0 95  1 0 5 h = 0 .2 0 0.2 0.4 0.6 0.8 x0 x1 x2 x3 x4 1 x5 En forma iterativa, hagamos la siguiente tabla: de i=0 hasta n=5 Recuerde que: xi+1=xi+h f(xi,yi) = xi - yi yi1  yi  hf(xi ,yi ) i xi 0 0 Yi 2 f(xi,yi) y i1 B  1 f(xi , y i )  f(xi 1 , y i 1 ) 2 Xi+1 f(x i 1 , y i 1 ) B  Yi+1=yi+h b -2 1.6 0.2 -1.4 -1.7 1.66 1 0.2 1.66 -1.46 1.368 0.4 -0.968 -1.214 1.4172 2 0.4 1.4172 -1.0172 1.21376 0.6 -0.61378 -0.81548 1.254104 3 0.6 1.254104 -0.654104 1.12328 0.8 -0.32328 -0.488694 1.156365 4 0.8 1.156365 -0.356365 1.085092 1.0 -0.08509 -0.2207286 1.112220 5 1.0 1.112222 Para i=5 y5=y(1)=1.112222 El valor exacto es 1.10364 1.10364  1.112222  100% 1.10364 E R  0.78%  ER  8.3. MÉTODOS DE RUNGE-KUTTA (RK) Los métodos de Runge-Kutta tienen la exactitud del esquema de la serie de Taylor sin necesidad de cálculo de derivadas superiores y consisten en obtener un resultado que se obtendría al utilizar un número finito de términos en una serie de Taylor de la forma: yi 1  yi  f ( xi , yi )h  f ' ( xi , yi )h 2 f " ( xi , yi )h 3 f ' ' ' ( xi , yi )h 4   2! 3! 4! +.......... con una aproximación en la cual se calcula yi+1 de la fórmula: Yi 1  yi  h( 0 f ( xi , yi )  1 f ( xi  1h, yi  b1h)  ... ....   p f ( xi   p h, yi  bp h) 96 (2) donde los ,  , b se determinan de modo que si expandiera f ( xi   j h, yi  bj h) con i  j  p en series de Taylor alrededor de (xi, yi), se observaría que los coeficientes de h, h2 , h3 , etc. Coincidirían con los coeficientes correspondientes de la ecuación (1). Después de efectuar los cálculos convenientes a la ecuación (2) para establecer los valores de  ,  , b , tenemos las siguientes expresiones de Runge-Kutta (RK): 1. Método RK de segundo orden: yi+1 = yi + h [k1 +k2] 2 donde: k1=f(xi, yi) k2=f(xi+h, yi+hk1) 2. Método RK de tercer orden: yi+1 = yi + donde: h [k1 +4k2 +k3] 6 k1=f(xi, yi) k2=f(xi+ h h , yi+ k1) 2 2 k3= f(xi+h, yi -hk1+2hk2) 3. MétodoRK de cuarto orden: yi+1 = yi + h [k1 +2k2 +2k3 +k4] 6 donde: k1=f(xi, yi) k2=f(xi+ h h , yi+ k1) 2 2 k3=f(xi+ h h , yi+ k2) 2 2 k4=f(xi+h, yi+hk3) otra forma de expresar 97 k1  hf ( x n , y n ) k h k 2  hf ( x n  , y n  1 ) 2 2 k h k 3  hf ( x n  , y n  2 ) 2 2 k 4  hf ( x n  h, y n  k 3 ) 1 y n 1  y n  (k1  2k 2  2k 3  k 4 ) 6  El método más utilizado es el de cuarto orden ya que coincide con los primeros cinco términos de la serie de Taylor lo cual significa gran exactitud sin cálculo de derivadas aunque haya que evaluar la función f(x,y) cuatro veces en cada subintervalo. Ejemplo: Aplíquese el método de Runge-Kutta de cuarto orden para f(x,y)=x-y donde x=0 hasta x=1 con 5 intervalos y siendo para x=0, y=2. Ejemplo: 1. Aplíquese el método de Runge-Kutta de cuarto orden para f(x,y)= x – y desde x=0 hasta x=1 Con 5 intervalos siendo para x=0 y=2. Solución: h (x0, y0)= (0, 2) 1 0  0, 2 5 para x=1 y=? x0=0 x1=0.2 x2=0.4 Aplicamos: yi+1 = yi + x3=0.6 x4=0.8 x5=1.0 h [k1 +2k2 +2k3 +k4] 6 Se calculan los valores de k1, k2, k3, k4 en cada iteración y se halla yi+1 i xi Yi K1 K2 K3 K4 Yi+1 0 0 2 -2 -1.7 -1.73 -1.454 1.6562 1 0.2 1.6562 -1.4562 -121058 -1.235142 -1.009172 1.410973 2 0.4 1.410973 -1.010973 -0.809876 -0.829985 -0.644976 1.246451 3 0.6 1.246451 -0.746451 -0.471806 -0.49927 -0.346797 1.148004 4 0.8 1.148004 -0.348004 -0.213204 -0.226684 -0.10268 1.103656 5 1.0 1.103656 98 para x=1.0, y=1.103656 Que al comparar con el valor exacto que es 1.10364, tenemos: ER  2. 1.10364  1.103656  100%  0.001% 1.10364 Resolver F ( X i , Yi )  1 1  X i  Yi 2 aplicando el método de Runge-Kutta. 2 Solución De la condición inicial del problema se tiene que X = 0, y Y = 1; además, h = 0.1. Sustituyendo estos valores en se obtiene: Llevando estos valores a (16) y el resultante a (12) se obtiene que para X = 0.1 la solución del problema es Los valores de las kipara este punto obtenido de la solución, son: luego Continuando de la misma forma se obtiene la solución que se muestra en la siguiente tabla: 99 8.4. X Y k1 k2 k3 k4 0.0 1.0000 0.5000 0.5516 0.5544 0.6127 0.1 1.0554 0.6126 0.6782 0.6823 0.7575 0.2 1.1236 0.7575 0.8431 0.8494 0.9494 0.3 1.2085 0.9492 1.0647 1.0745 1.2121 0.4 1.3158 1.2119 1.3735 1.3896 1.5872 0.5 1.4545 1.5868 1.8234 1.8517 2.1509 MÉTODO DE DIFERENCIAS FINITAS Las diferencias divididas finitas se sustituyen por las derivadas en la ecuación original. Así una ecuación diferencial se transforma en un conjunto de ecuaciones algebraicas lineales simultaneas. Primera derivada: dy yi1  yi  dx x Segunda derivada: d 2 y yi 1  2 yi  yi 1  dx 2 x 2 Resolver la ecuación mediante diferencias finitas d2 y  2  106 y  7.5  107  (75  x) 2 dx i 1 i 1 i d 2 y yi 1  2 yi  yi 1  dx 2 (x) 2 yi 1  2 yi  yi 1  2  106 yi  7.5  107 xi (75  xi ) 2 (x) i 1 x0 i2 i 3 i4 x  25 x  50 x  75 x  25 , x  0 a x  75 con x  25 . 100 x0  0 x1  x0  x  0  25  25 x 2  x1  x  25  25  50 x3  x2  x  50  25  75 x  0 y1  0 y3  2 y 2  y1  2  10 6 y 2  7.5  10 7 x2 (75  x2 ) 2 (25) 0.0016 y1  0.003202 y 2  0.0016 y 3  7.5  10 7 ( 25)(75  25) 0.0016 y1  0.003202 y 2  0.0016 y 3  9.375  10 4 y 4  2 y3  y 2  2  10 6 y3  7.5  10 7 x3 (75  x3 ) 2 (25) 0.0016 y 2  0.003202 y 3  0.0016 y 3  7.5  10 7 (50 )(75  50 ) 0.0016 y 2  0.003202 y 3  0.0016 y 3  9.375  10 4 x  75 y4  0  0 0 0   y1  0  1  y   0.0016  0.003202  0.0016 0   2  9.375  10 4     0 0.0016  0.003202 0.0016  y3  9.375  10 4       0 0 1   y 4  0  0   y1  0  y     2    0.5852  y 3   0.5852       y 4  0 8.5. SISTEMA DE ECUACIONES DIFERENCIALES DE PRIMER ORDEN Ya sea que estemos afrontando un problema de ingeniería cuya solución implique la resolución de una ecuación diferencial de orden n, o uno que implique la resolución de un sistema de ecuaciones diferenciales ordinarias de primer orden, nos enfrentaremos a la necesidad de resolver un sistema que podremos expresar como sigue: 101 dy1  f1  x , y1 , y 2 ,............, y n  dx dy 2  f 2  x , y1 , y 2 ,............, y n  dx . . dy n  f n  x , y1 , y 2 ,............, y n  dx Por supuesto, la solución de tal sistema requiere de que se conozcan las n condiciones iniciales en el valor inicial del intervalo correspondiente a x. Todos los métodos vistos anteriormente para simples ecuaciones pueden extenderse a la resolución de sistemas como el anterior. El procedimiento para resolver un sistema de ecuaciones simplemente involucra aplicar las técnicas conocidas para cada ecuación en cada paso, antes de proceder con el siguiente. Esto quedará claramente ilustrado con el siguiente ejemplo donde hemos aplicado el método de Euler. Ejemplo 1:Resolución de un sistema de EDOs mediante el método de Euler. Resolver el siguiente conjunto de ecuaciones diferenciales ordinarias: dy1  - 0, 5  y1 dx dy 2  4 - 0, 3  y 2 - 0,1  y1 dx Resolveremos el sistema en el intervalo entre x = 0 y x = 2, con condiciones iniciales en x = 0, y1 = 4 , y2 = 6. Utilizaremos un paso h = 0.5. Se implementa el método de Euler para cada variable mediante la ya conocida expresión: yi 1  yi  f ( xi , yi )  h Primero calculamos las pendientes: dy1  f1 ( 0 , 4, 6)  - 0.5  4  - 2 dx dy 2  f 2 ( 0 , 4, 6)  4 - 0.3  6 - 0.1  4  1.8 dx Y luego los valores de la función para el primer paso: 102 y1  0.5  y1  0  f1 ( 0 , 4 , 6)  h  4   -2 0.5  3 y2  0.5  y2  0  f 2 ( 0 , 4 , 6)  h  6  1.8 0.5  6,9 Para un segundo paso volvemos a calcular las pendientes: dy1  f1 ( 0.5 , 3, 6.9)  - 0.5  3  - 1.5 dx dy 2  f 2 ( 0.5 , 3, 6.9)  4 - 0.3  6.9 - 0.1  3  1.63 dx Y luego los valores de la función para el segundo paso: y1 1.0  y1  0.5  f1 ( 0.5 , 3, 6.9)  h  3   -1.5  0.5  2.25 y2 1.0  y2  0.5  f2 ( 0.5 , 3, 6.9)  h  6.9  1.63  0.5  7.715 Y así continúa el cálculo hasta el final. Los resultados se resumen en la siguiente tabla: x y1 y2 0.0 4.000000 6.000000 0.5 3.000000 6.900000 1.0 2.250000 7.715000 1.5 1.687500 8.445250 2.0 1.265625 9.094870 8.6 PROBLEMAS 1. Dada la ecuación diferencial f(x,y)=yx2-y con x=0 a x=2 siendo y(0)=1 y n=4 intervalos resuelva por el método sencillo de Euler 2. Dada y’=f(x,y) = 2x3 –3x2 entre x=0 y x=1 siendo f(0)=1 y con n=2, evalúe utilizando a)el método sencillo de Euler. 3. Dada y’=f(x,y)= -2x3 + 12x2 –20x + 8.5 resuelva utilizando el método RK de cuarto orden desde x=1 hasta x=3 considerando 4 intervalos y siendo y(1)=3. 4. Utilice el método de Euler Modificado para resolver: a) dy/dx=2x3 –2xy con y(0)=0 y(2.5)=? 103 h=0.5 b) y’=2x2 –3y2 con y(1)=0.5 y(2)=? h=0.2 5. Resolver mediante el método euler dy  100000( y  e  x )  e  x dx y (0)  0 , use the implicit Euler method to obtain the solution from x  0 to 0,3 using the step-size of 0.1 6. Un tanque de 600 galones de capacidad contiene inicialmente 200 galones de salmuera con 25 lb. de sal. Salmuera con 2 lb. por galón entra al tanque a un flujo de 13 galones/ s. La salmuera mezclada en el tanque fluye hacia fuera a un flujo volumétrico de 8 galones/ s. ¿Cuál es la cantidad de sal y la concentración cuando el tanque se Encuentra lleno? Utilizar el método de runge kutta de segundo orden. 7. Las plantas de desalinización se usan para purificar el agua de mar , para que pueda beberse . El agua de mar contiene disuelto 8 g de sal /kg. Y es bombeada hacia un tanque de mezcla a razón de 0,5 kg./ min. Debido a una falla en el diseño del equipo, el agua se evapora del tanque a razón de 0,5kg/min. .La solución salina sale del tanque a razón de 10 kg/min. Supóngase que el balance de la disolución es agua pura. Si el tanque se llena inicialmente con 1 000 kg. de solución a la entrada. ¿ determine la concentración a la salida del tanque en función del tiempo hasta que el tanque quede vació.? 8. Las siguientes ecuaciones definen las concentraciones de tres reactores dCA  20CACC  2CB dt dCB  20CACC  2CB dt dCC  20CACC  2CB  0,2CC dt Si las condiciones iniciales son CA = 500, CB = 0 y Cc = 500 halle las Concentraciones para los tiempos que van desde 0 a 30s usando un h = 3 9. La descarga de un fluido por gravedad por un orificio por el fon deo del tanque puede ser modelado mediante el siguiente sistema d e ecuaciones diferenciales. 104 dv  0,0010 y  0,00205 v 2 dt dy 1  0 , 0624 (1  ) dt 100  v donde v es la velocidad del fluido, y nivel del liquido las condiciones iniciales son de v0 = 2,5 ; y0 = 20 cual son los valores de v, y para t = 5 con h igual 0,5 10. Par un circuito eléctrico los resistores no obedecen la ley de Ohm, y el circuito dinámico se describe por la relación siguiente l di i 3 i  R ( )   0 dt I  I i = corriente eléctrica ; I = corriente de referencia igual 1; R resistencia eléctrica 2 ; Resolver la ecuación diferencial si i(0) = 0,6 A para t= 30s con h = 3s 11. Determine usado técnicas numéricas la ecuación diferencial ordinaria de segundo orden d2x dx  3  10 x  0 2 dt dt La condición inicial para la ecuación para t = 0  x =1 Determine los valores de x entre 1,2 usando h = 0,1 12. Un tanque de 2 000 L contiene, en el inicio 400L de agua pura . Comenzando en t = 0 una solución acusa que contiene 1 g/ L de cloruro de potasio fluye hacia el tanque a razón de 8L/s y al mismo tiempo, comienza a fluir una corriente de salida a razón de 4L/s .El contenido del tanque esta mezclado perfectamente y la densidad de la corriente de alimentación y salida puede considerarse constante. Calcule la concentración de cloruro de potasio en el tanque en el momento en que este rebosa. 13. Un Reactor CSTR como se muestra en la figura, la concentración de alimentación (CA0) of 0,5 mole/m3. Determine la concentración al cabo de 0,5h 105 rA  V  k1C A (1  k2 C A ) dC A  F (C A0  C A )  (V )(rA ) dt V  2, 0m3 F  1, 0m3 / h k1  1, 0h 1 k2  1, 0m3 / mole 14. Se tiene un intercambiador de calor de tubos concéntricos en contracorriente y sin cambio de fase .Las ecuaciones que describen el intercambiador de calor en ciertas condiciones de operación son: dTB  0.03 x TS  TB  dx dTS  0.04 x TS  TB  dx CONDICIONES TB (0)  20 0 C TS (0)  1000 C TB (3)  ¿? TS (3)  ¿? Evalué la temperatura para x= 3 m 15. Resolver la siguiente ecuación mediante diferencias finitas d 2u 1 d u u   2  0, u(5)  0,08731, u(8)  0.0030769 dr 2 r dr r 106 16. El siguiente diagrama ( Figura.18) muestra tres tanques continuos con agitación conectados en serie Se ha disuelto1500g de Na 2SO4 en el primer tanque, llenando los otros tanques con solvente puro, e iniciando después un flujo de 40L/ a través del sistema. Calcular el tiempo para alcanzar una concentración de 0,01g/l en el tanque tres. Fig Nº16: Sistema de tanques de mezclado. 17. Determine el valor de y(1) resolviendo y t   0,05 y t   0,15 y t   0, y 0   0, y0  1 utilizando el método de Runge-Kutta de segundo orden, con h = 0,5 18. El circuito que se muestra en la figura E9.5 tiene una autoinductancia de L = 50H, una resistencia de R = 20 ohms y una fuente de voltaje de V = 10 volts. Si el interruptor se cierra en el instante t = 0, la corriente I(t) satisface la ecuación: L d I t   RI t   E , dt I 0   0 Determine el valor de la corriente para 0 < t  10 segundos, mediante el método de Runge-Kutta de segundo orden, con h = 0.1. L SW i E R Figura 19: Circuito electrico 107 19. El agua que sale del tanque entra a otro de 20 galones, en cual también se vierte agua pura a razón de 3 galones/minuto y se mezcla bien. La concentración de sal en el segundo tanque satisface y 2 t    3 2 y 2 t   y1 , 20 20 y 2 0   10 donde y1(t) es la concentración de sal del tanque de 50 galones del problema anterior. Utilice el método de Runge–Kutta de segundo orden para determinar cuando alcanza su máximo la concentración de sal en el tanque de 20 galones. Suponga que el segundo tanque tiene agua pura en el instante t = 0. 20. Un tanque de 50 galones de agua contiene sal con una concentración de 10 onzas/galón. Con el fin de diluir el contenido de sal, se suministra agua pura a razón de 2 galones/minuto. Si el depósito tiene una mezcla uniforme y la misma cantidad de agua que entra sale del depósito cada minuto, la concentración de sal satisface y1 t    2 y1 , 50 y1 0   10 donde y1(t) es la concentración de sal en onzas/galón y t es el tiempo en minutos. Utilice el método de Runge–Kutta de segundo orden con h = 1 minuto para determinar cuánto tiempo debe transcurrir para que la concentración de la sal sea 1/10 de su valor inicial. b) El agua que sale del tanque entra a otro de 20 galones, en cual tambien se vierte agua pura a razon de 3 galones/minuto y se mezcla bien. La concentracion de sal en el segundo tanque satisface y 2 t    3 2 y 2 t   y1 , 20 20 y 2 0   10 donde y1(t) es la concentración de sal del tanque de 50 galones del problema anterior. Utilice el método de Runge–Kutta de segundo orden para determinar cuando alcanza su máximo la concentración de sal en el tanque de 20 galones. Suponga que el segundo tanque tiene agua pura en el instante t = 0. 108 21. Resolver el sistema de EDOs mediante el método de Runge – Kutta dy1  - 0.5  y1 dx clásico de cuarto orden. dy 2  4 - 0.3  y 2 - 0.1  y1 dx Resolveremos el sistema en el intervalo entre x = 0 y x = 2, con condiciones iniciales en x = 0, y1 = 4 , y2 = 6. Utilizaremos un paso h = 0.5. 109 9. ECUACIONES DIFERENCIALES PARCIALES (EDP) INTRODUCCIÓN Una EDP es una ecuación que tiene como incógnita a una función de dos o más variables y que involucra a una o más de sus derivadas parciales. El orden de una EDP es el de la derivada con mayor orden en la ecuación. La linealidad de las ecuaciones se establece como sigue: Si los coeficientes dependen sólo de las variables independientes entonces a la ecuación se le denomina lineal. Si además dependen de la propia función o de alguna de sus derivadas parciales entonces la ecuación es no lineal. EJEMPLO  Ecuación diferencial parcial de transferencia de calor unidimensional en estado no estacionario. T T 2  2 t x T: Temperatura en K , t es el tiempo en s,  es la difusivilidad térmica en m2/s  Ecuación de transferencia de calor en estado estacionario bidimensional T 2 T 2  2 x 2 y T es la temperatura en K , x,y ejes coordenados 9.1 PROBLEMAS 1. resuelva la ecuación diferencial parcial de transferencia de calor en dos dimensiones.  2T  2T  0 X 2 Y 2 0  X 2 , 0  Y  2 , X  Y  0,5 T( 0, Y) = 100 , T( X ,0 ) = 60 , T( 2,Y) = 100 , T(X, 2) = 60 2. Resuelva la siguiente ecuación diferencial de transferencia de calor dada por 110  2T  2T  0 X 2 Y 2 T = 40 . 0  y 1 , x=0 T = 120 0 y  1 x =4 T = 40 0 X  4 y=0 T= 50 0 X  4 y=1 Utilizando y  0,25 y x  1 111 10. APLICACIONES DE INGENIERÍA QUÍMICA EN POLYMATHMATHCAD Y MATLAB 10.1 PROBLEMAS CON POLYMATH 1. Determine el volumen Molar del gas dióxido de carbono mediante la ecuación de VAN DER WAALS . a   P  V  b   RT V *V   a= 3,592 b= 0,04267 R = 0,082 PARA T= 300K , P =100 ATM f(v) = (P + a/(V2))*(V-b) –R*T = 0 POLYMATH Results NLE Solution Variable V Value 0.0793969 a 3.592 b 0.04267 T 300 R 0.082 P 100 f(x) -1.776E-13 Ini Guess 0.2505 2. En la batería de reactores de mezcla completa que se muestra en el diagrama, se lleva a cabo una reacción isotermica de Segundo orden a volumen constante si el flujo volumétrico a cada reactor es constante de 25 l/min. calcular la concentración de salida de cada reactor 112 Figura Nº17: Bateria reactores Balance de materia en cada reactor: Reactor N°1: q*CA0 –q*CA1 –K*CA12*V1 =0 Reactor N°2 q*CA1-q*CA2 –K*CA22*V2 =0 Reactor N°3 q*CA2-q*CA3 –K*CA32*V3 =0 F(CA1) = q*(CA0 –CA1) –K*CA12*V1 =0 F(CA2) = q*(CA1-q*CA2) –K*CA22*V2 = F(CA3) = q*(CA2-q*CA3) –K*CA32*V3 =0 Variable Value f(x) CA1 0.6031317 -7.105E-15 2 CA2 0.2869484 -1.332E-14 1 CA3 0.1725794 -2.944E-13 0.8 V 1200 K 0.08 CA0 2 Q 3. En s 25 un reactor químico se obtiene glucosa a partir de la hidrólisis de almidón. La ecuación diferencial que gobierna el proceso para la concentración del almidón es . V*Dca/dt = q*CA0 - q*CA + ra ra = -k*CA2V CA0: Concentración de la alimentación a la entrada = 6,3 mol/m3 113 CA: Concentración del almidon a la salida V : volumen del reactor = 500m3 Q : caudal volumetrico 100 m3/h ra : velocidad de transformación de almidon en glucosa 3 K= 0.02 m /h*mol . Calcular la composición de salida para t = 20h ODE Report (RKF45) Differential equations as entered by the user [1] d(CA)/d(t) = q*CA0 /V - q *CA/V -K*CA2 Explicit equations as entered by the user [1] q = 100 [2] V = 500 [3] CA0 = 6.3 [4] K = 0.02 Independent variable variable name : t initial value : 0 final value : 20 Variableinitial valueminimal valuemaximal valuefinal value t 0 CA 6.3 q 20 20 4.38179 6.3 4.38179 100 100 100 100 V 500 500 500 500 CA0 6.3 6.3 6.3 6.3 0.02 0.02 0.02 0.02 K 0 4. La concentración saturada de oxigeno disuelto en agua como función de temperatura y de la concentración de cloro se lista en la siguiente tabla. *Determine la ecuación que correlacione C=f(T) de grado 2 y 3  Evaluar la concentración para T = 20°C 114 Temperatura Cloro( 10 000mg/l) 5 11,6 10 10.3 15 9,1 20 8,2 25 7,4 30 6,8 REGRESION LINEAL Model: C = a0 + a1*T + a2*T2 Variable Value 95% confidence a0 13.11 0.1595539 a1 -0.3195 0.0208769 a2 0.0036429 5.839E-04 REGRESION POLINOMICAt Model: C = a0 + a1*T + a2*T2 + a3*T3 Variable Value 95% confidence a0 13.133333 0.5263102 a1 -0.3253704 0.1199362 a2 0.0040317 0.0076758 a3 -7.407E-06 1.451E-04 5. La siguiente secuencia de reacciones se llevan a cabo en estado no estacionario según la reacción: Determinar la ecuación de cada componente de la reacción así como los productos en función al tiempo. Construya la tabla tiempo vs concentración en un paso de 1 a 1 hasta 10 minutos, asuma que las concentraciones iniciales son A (0)=1 mol/L, B (0)=1 mol/L, C (0)=1 mol/L. Datos: 115 K1=0.5 min-1 K2=1.0 min-1 Solución: formando las ecuaciones diferenciales en función de la reacción. C1   k1C1 t C2  k1C1  k2C2 t C2  k 2 C2 t FIGURA 18: Datos de ingreso en Polymath FIGURA 19: Gráfica concentración vs tiempo 116 6. Resolver el siguiente sistema: y   k1 x  k2 y t x  k1 x t z  k2 y t Con valores iníciales x (0)=1, y (0)=0, z(0)=0 con t=0 a t=3 y k 1=1, k2=2 Solución: FIGURA 20: Ingreso de Datos en Polymath 117 FIGURA 21: Resultados con el Software Polymath FIGURA 22: Grafica de resultados de la Ecuación diferencial. 118 7. Resolver las siguientes ecuaciones: dC1  0,12C1  0, 02C3  1 dt dC2  0,15C1  0,15C2 dt dC3  0, 025C2  0, 225C3  4 dt dC4  0,1C3  0,1375C4  0, 025C5 dt dC5  0, 01C1  0, 01C2  0, 04C5 dt Con C1 (0)=1, C2 (0)=1, C3 (0)=1, C4 (0)=1 y C5 (0)=1 de t=0 a t=1 Solución: FIGURA 23: Ingreso de Datos en Polymath 119 FIGURA 24: Resultados Polymath FIGURA 25: Grafico Polymath 120 10.2 PROBLEMAS CON MATHCAD 1. Para flujo turbulento para todas las tuberías, el Instituto de Hidráulica de Estados Unidos y la mayoría de ingenieros consideran la ecuación de Colebrook como la más aceptable para calcular f. La ecuación es:  e 0, 251  0,86 ln    0,36d Re f f  1    Donde: Re : l número de Reynolds (adimensional) e : aspereza o rugosidad de la tubería (unidad de longitud) d :, diámetro de la tubería (unidad de longitud) Obtener el factor de fricción para un fluido con un Reynolds de 3E4 que fluye en una tubería con un diámetro de 0,1 m y una rugosidad de 2,x10-3 m. Solución con mathcad Ingresando datos: E = 0,002 D = 0,1 Re = 3104 f = 0,001 Given  E 0, 251  0,86 ln    f  0,36 D Re f 1    Calculo de la fricción F = Find (f) f = 0,055 de la figura efectuar un balance de masa y determine los flujos 121 Figura Nº26: Sistema de columna de separación 0.07x1 +0.18x2 + 0.15x3 +0.24x4 = 10.5 0.04x1 + 0.24x2 + 0.1x3 +0.65x4 =17.5 0.54x1 + 0.42x2 + 0.54x3 + 0.10x4 = 28 0.35x1 + 0.16x2 + 0.21x3 + 0.01x4 = 14  0.07 0.04 A    0.54  0.35  0.18 0.15 0.24  0.24  0.1 0.65    0.160 0.21 0.01  0.42 0.54 0.1 122  10.5     26.25  17.5     B   17.5  solucion   28   8.75   14   17.5      2. La ecuación de Antonie para el cálculo de la presión de vapor (P*) de un compuesto esta dado por: Donde: P*: presión de vapor (mmHg) T, temperatura (ºC) Las constantes para el benceno y el hexano son: Compuesto A B C benceno 6.89745 1260.350 220.237 Hexano 6.87773 1171.530 224.366 Determinar la temperatura de ebullición de la mezcla a 1 y a 5 atm. Solución log ( P*) A B C T sean las fracciones molares xa = 0,5 xb = 0,5 constantes para el benceno A1  6.89745 B1  1206.350 C1  220.237 constantes para el hexano A2  6.87773 B2  1171.530 C2  224.366 para una presión de 1atm en mmHg 123 Pt  760 asumiendo una temperatura inicial de 30ºC T  30 Given  A1 B1   A2 B2      C1  T  C2  T    Pt Xa 10  Xb 10 T  Find(T) T  73.975 3. Resolver el siguiente sistema de ecuaciones no lineales con mathcad 2 k1 k2 x y 2 ( 1  x)  ( 1  x  y ) y ( x  y) ( x  y)  ( 1  x  y) Solución Especifique las constantes de equilibrio Tome dos valores iniciales Given 2 k1 k2 x  y 2 ( 1  x)  ( 1  x  y ) y (x  y ) (x  y )(1  x  y ) vec  Find ( x  y ) vec   0.652     0.198  124 k1  0.6 k2  0.2 x  0.8 y  0.5 4.. Mediante la ecuacion de Vand Der Waals calcular el volumen molar a P = 60 atm y T = 500K.  P  a   ( V  b ) R T  2 V   ecuacion de Van der Waals Definiendo las constantes R  0.082 Para el amoniaco Tc  405 Pc  111.3 Constante a y b 27  R  Tc 2 a   64  2   Pc b  R Tc 8 Pc Especificando las constantes R  0.082 Para el amoniaco Constante a y b a  Tc  405  R 2  Tc  64  Pc 27 Evaluando las constantes a y b 2 a  Especifique la presion y temperatura b  R  Tc 8  Pc b  T  500 Asumiendo un valor inicial( comportamiento ideal) R T V  V  P Dado  P  a  (V  b )  2  V   calculo del volumen Molar V  Find    Pc  111.3 (V) V  125 R T P  60 5. Encontrar el volumen molar que ocupa el gas empleando la ecuación de Redlinh Kwong RT a P  V  b V (V  b) T  R 2TC5 / 2   RTc  a  0, 42747   b  0, 08664    Pc   PC   Atm  L    mol  gK  P en Atm V en L/mplg T en K R  0, 082  Pc: presión critica de amoniaco 111,3 atm. TC: temperatura critica de amoniaco 405K ECUACION DE ESTADO DE R-K Ingresando las contantes: Ingresando los datos: R  0.082 P  60 Pc  111.3 Tc  405 T  500 Calculo de a yb: 5   R 2  Tc 2 a  0.42747   Pc    b  0.08664   a  b  Asumiendo comportamiento ideal V  V Dado P R T (V  b)  V  Find ( V) a V ( V  b )  T V 126 R T P R  Tc Pc    6. A un sistema de separación FLASH ( Fig. 21) ingresa 100 moles/h de una mezcla VAPOR(V) 2 yi 1 Alimentacion F zi 3 xi Liquido(L) Figura Nº27:Sistema de separación de componentes Alimentación al Flash( Mol/h) Zi: Fraccion molar de los componentes de la mezcla i = 1,2,3, …NC V: Caudal de salida del vapor ( mol/h) yi : Fraccion molar de los componentes en la fase vapor L: Caudal De salida de la fase liquida ( mol/h) xi: Fraccion molar de los componentes en la fase liquida Balance de materia global: F = V + L Balance de componentes : Zi*F = xi*L + yi*V Equilibrio Liquido vapor : yi = Ki *xi V F Remplazando en las anteriores ecuaciones expresiones Fraccion vaporizada : x  zi  ( Ki  1)  1 127 Se obtiene las siguientes yi  Ki * zi  ( Ki  1)  1 Ni Zi( Ki  1)   (Ki  1)  1  0 i 1 Teniendo en cuenta los datos de la siguiente tabla Calcular la fracción vaporizada así como los flujos en cada corriente con sus respectivas composiciones COMPONENTES FRACCION MOLAR C1 0.05 16.25 C2 0.15 5.25 C3 0.25 1.99 C4 0.2 0.75 C5 0.35 0.29 Ingresando los datos: F  100 k3  1.99 k1  16.25 k4  0.75 z1  0.05 z3  0.25 K k2  5.25 k5  0.29 z2  0.15 z4  0.2 Asumiendo un valor: z5  0.35   0.5 Dado z1 ( k1  1)  ( k1  1)  1    Find( ) z2 ( k2  1)  ( k2  1)  1  z3 ( k3  1)  ( k3  1)  1  128  z4 ( k4  1)  ( k4  1)  1  z5 ( k5  1)  ( k5  1)  1 0 7. Encontrar el volumen molar que ocupa el gas empleando la ecuación de estado de Redlinh-Kwong. A una presión de 14 bar y una temperatura de 333K RT a  V  b V (V  b) T P Ingresando los datos: T  333 a b f ( V) P T b  44.897 8 a  1.5641410  R  83.14 P  14   b  2 b  R T  P    V  R T  V2  V3  P P T  a Se observa un plinomio de grado 3 Ingrese los coeficientes en una matriz de 4*1      2  b  V           P T  b  R T a    P P T   R T    P  1  a b S  polyroots ( V) 8.  71.299    229.996 S   3  1.676  10  La temperatura media logarítmica de un intercambiador de calor a contracorriente esta dado por: LMDT  T1  t2   T2  t1  T t  ln  1 2   T2  t1  129 Para este sistema la temperatura media logarítmica debe ser de 50°C y el fluido caliente se alimenta a 100°C y sale a 60 °C, mientras que el fluido frió se alimenta a 15°C ¿A qué temperatura sale del intercambiador el fluido frió? Ingrese los datos : T1  15 T3  100 T4  60 Asumiendo un valor inicial: LMTD  50 T2  50 Dado LMTD ( T3  T2)  ( T4  T1) ln T3  T2    T4  T1  T2  Find( T2) T2  9. El factor de fricción para diferentes números de Reynolds empleando la ecuación de Colebrook   1 2,51   0,86 ln    f  3, 7d Re f  Obtener el factor de friccion para Re de 10 000 , e = 0,002cm y diámetro de 2cm CALCULO DEL FACTOR DE FRICCION Ingresando datos: Re  10000 Asumiendo para f : f  0.001 Dado 0.86 ln 1 e  3.7 d f    Re f  2.51 f  Find( f ) f 130 e  0.002 d  2 10.3. PROBLEMAS CON MATLAB 1. A un reactor ingresa una mezcla de gases de 8% de dióxido de azufre y el 12% de oxigeno y la diferencia de nitrógeno, y se desarrolla la siguiente reacción. 1 SO2  O2  SO3 2 Calcular la composición en el equilibrio a presión constante de 2 atm y la constante de equilibrio KP es de 160. Algorritmo de Solucion en Matlab: clear all, close all, clc disp('H2SO4'); SO2=input('Ingrese moles de SO2: '); O2=input('Ingrese moles de O2: '); SO3=input('Ingrese moles de SO3: '); N2=input('Ingrese moles de N2: '); P=input('Ingrese Presion Total (atm): '); KP=input('Ingrese Constante de Equilibrio: '); syms x; disp('Moles en Equilibrio: '); ESO2=SO2-x; EO2=O2-0.5*x; ESO3=SO3+x; disp(ESO2); disp(EO2); disp(ESO3); ET=ESO2+EO2+ESO3+N2; f=KP-((ESO3/ET)*P)/((ESO2/ET)*((EO2/ET)^0.5)*P^1.5); 131 xai=input('Ingrese Limite Inferior: '); xbi=input('Ingrese Limite Superior: '); tol=input('Ingrese Tolerancia: '); f=inline(f); i=1; ea(1)=100; if f(xai)*f(xbi) < 0 xa(1)=xai; xb(1)=xbi; f1=f(xai); xr(1)=(xa(1)+xb(1))/2; fprintf('Iter. x1 f1 x2 Error aprox \n'); fprintf('%2d \t %11.7f \t %11.7f \t %11.7f \n',i,xa(i),f1(i),xr(i)); while abs(ea(i)) >= tol, if f(xa(i))*f(xr(i))< 0 xa(i+1)=xa(i); xb(i+1)=xr(i); end if f(xa(i))*f(xr(i))> 0 xa(i+1)=xr(i); xb(i+1)=xb(i); end xr(i+1)=(xa(i+1)+xb(i+1))/2; ea(i+1)=abs((xr(i+1)-xr(i))/(xr(i+1))*100); 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; 132 end else fprintf('No existe una raíz en este intervalo'); end Resultados: TABLA 12 Resultados en Matlab 2. Se sabe que la temperatura media logarítmica de un intercambiador de calor contracorriente está dada por:  ΔT  AQ T = LMTD = ent cal sald - T frio  - Tcalsald -T frioentr  sald  Tcalent - T frio ln  sald entr  T -T frio  cal 133    Se tiene el siguiente intercambiador: entr = 15ºC T frio INTERCAMBIADOR DE CALOR sald = T frio Tcalsald= 60ºC t Tcalent = 100ºC Se sabe que en este sistema la temperatura media logarítmica debe ser de 50ºC. El fluido caliente se alimenta al sistema a 100ºC, mientras que el fluido frío se alimenta a 15ºC,. ¿a que temperatura sale del intercambiador de fluido frío? Algoritmo de Solución en Matlab: clear all, close all, clc disp ('Temperatura en un Intercambiador de Calor') xo=input('Temperatura Inicial ='); n=input ('Numero de Iteraciones='); salida=ones(n,3); % matiz de salida de datos for i=1:n x1=xo-[(45*exp((55-xo)/50)+xo-100)]/[1+(-45/50)*(exp((55-xo)/50))]; vsal=[xo;x1]; ea=[[abs((x1-xo)/x1)]]; % error xo=x1; salida(i,1)=i; salida(i,2)=x1; salida(i,3)=ea; end disp('Iter. T(i) Error'); disp(num2str(salida)); Resultados: 134 TABLA 13 Resultados Temperatura en Matlab 3. Un gas se encuentra a una presión absoluta de 13.76 bar y una temperatura de 333 K. Encontrar el volumen molar ocupa el gas empleando la ecuación de estado de Redlich-Kwong. P= RT a - 1/2 V - b T V V + b  Para este compuesto las constantes son: P = 13.76 atm T = 333ºk a = 1.5614 x 108 (cm6 bar/(g mol)2 k1/2) b = 44.897 (cm3/ g mol) R = 83.4 (cm3 bar /g mol k) Algoritmo de Solución en Matlab: clear all, close all, clc disp('Ecuacion de Redlich-kwong 2'); a=input('Ingrese a: '); b=input('Ingrese b: '); T=input('Ingrese T(K): '); P=input('Ingrese P(bar): '); R=input('Ingrese R(cm^3.bar/K.mol): '); syms x 135 f=R*T/(x-b)-a/((sqrt(T))*x*(x+b))-P; x0=(R*T)/(P); disp('Primer punto inicial: '); disp(x0); xai=x0; xbi=input('Ingrese Limite Superior: '); tol=input('Ingrese Tolerancia: '); f=inline(f); i=1; ea(1)=100; if f(xai)*f(xbi) < 0 xa(1)=xai; xb(1)=xbi; f1=f(xai); xr(1)=(xa(1)+xb(1))/2; fprintf('Iter. V1 F1 V2 Error aprox \n'); fprintf('%2d \t %11.7f \t %11.7f \t %11.7f \n',i,xa(i),f1(i),xr(i)); while abs(ea(i)) >= tol, if f(xa(i))*f(xr(i))< 0 xa(i+1)=xa(i); xb(i+1)=xr(i); end if f(xa(i))*f(xr(i))> 0 xa(i+1)=xr(i); xb(i+1)=xb(i); end xr(i+1)=(xa(i+1)+xb(i+1))/2; 136 ea(i+1)=abs((xr(i+1)-xr(i))/(xr(i+1))*100); 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 else fprintf('No existe una raíz en este intervalo'); end Resultados: TABLA 14 Resultados volumen en Matlab 137 4. El factor de fricción (f) para el flujo turbulento en una tubería está dado por la correlación de Colebrook: ε/D 1 = - 0.86 ln f 2.54  +   3.4 Re   f  Donde Re = es el número de Reynolds (adimensional) , es la aspereza o rugosidad de la tubería (unidad de longitud) D, es el diámetro de la tubería (unidad de longitud) Obtener el factor de fricción para un fluido con un Reynolds de 3E4 que fluye en una tubería con un diámetro de 0.1 m y una rugosidad de 0.0025m. Algoritmo de Solución en Matlab: clear all, close all, clc disp('Friccion en un Flujo Turbulento'); Re=input('Ingrese el Numero de Reynolds: '); e=input('Ingrese la Rugosidad: '); D=input('Ingrese el Diametro de Tuberia: '); syms x; f=-sqrt(1/x)-0.86*log((e/D)/3.4+2.54/(Re*sqrt(x))); xai=input('Ingrese Limite Inferior del Intervalo: '); xbi=input('Ingrese Limite Superior del Intervalo: '); tol=input('Ingrese Tolerancia: '); f=inline(f); i=1; ea(1)=100; if f(xai)*f(xbi) < 0 xa(1)=xai; xb(1)=xbi; 138 f1=f(xai); xr(1)=(xa(1)+xb(1))/2; fprintf('Iter. f1 F1 f2 Error aprox \n'); fprintf('%2d \t %11.7f \t %11.7f \t %11.7f \n',i,xa(i),f1(i),xr(i)); while abs(ea(i)) >= tol, if f(xa(i))*f(xr(i))< 0 xa(i+1)=xa(i); xb(i+1)=xr(i); end if f(xa(i))*f(xr(i))> 0 xa(i+1)=xr(i); xb(i+1)=xb(i); end xr(i+1)=(xa(i+1)+xb(i+1))/2; ea(i+1)=abs((xr(i+1)-xr(i))/(xr(i+1))*100); 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 else fprintf('No existe una raíz en este intervalo'); end Resultados: TABLA 15 Resultados volumen en Matlab 139 5. Determinar volumen molar del oxigeno mediante la ecuación del VAN DER WAALS a    P  2 v  b   RT v   p = 100 Atm., T = 700 K para un gas que tiene a = 1,36 b = 0,0318 Algoritmo de Solución en Matlab: clear all, close all, clc disp ('Ecuacion de Van der Walls') T=input('Ingrese T(K): '); P=input('Ingrese P(atm): '); R=input('Ingrese R(L.atm/K.mol): '); 140 disp('Volumen Inicial: '); xo=R*T/P; disp(xo); n=input ('Numero de Iteraciones='); salida=ones(n,3); % matiz de salida de datos for i=1:n x1=xo-(100*xo^3-60.58*xo^2+1.36*xo-0.043248)/(300*xo^2121.16*xo+1.36); vsal=[xo;x1]; ea=[[abs((x1-xo)/x1)]]; % error xo=x1; salida(i,1)=i; salida(i,2)=x1; salida(i,3)=ea; end disp('Iter. V(i) Error'); disp(num2str(salida)); Resultados: TABLA 16 Resultados volumen en Matlab 141 6. La siguiente secuencia de reacciones se llevan a cabo en estado no estacionario según la reacción: k1 k2 A   B  C Determinar la ecuación de cada componente de la reacción así como los productos en función al tiempo. Construya la tabla tiempo vs conc

372 downloads 160 Views 2MB Size

Recommend Stories

Story Transcript

UNIVERSIDAD NACIONAL DEL CALLAO VICE- RECTORADO DE INVESTIGACION

INSTITUTO DE INVESTIGACION DE LA FACULTAD DE INGENIERIA QUIMICA

INFORME FINAL DEL PROYECTO DE INVESTIGACION TITULADO

TEXTO " METODOS NUMERICOS UTILIZANDO POLYMATH, MATHCAD Y MATLAB APLICADOS A INGENIERIA QUIMICA” PRESENTADO POR

ING. JUAN MEDINA COLLANA

(Del 1 de Marzo del 2010 al 29 de Febrero 2012 Resol. N° 1312-2010 R)

1

ÍNDICE I.

RESUMEN

4

II.

INTRODUCCIÓN

5

III.

MARCO TEÓRICO

8

IV.

MATERIALES Y MÉTODOS

8

V.

RESULTADOS

9

VI.

DISCUSIÓN

9

VII.

REFERENCIALES

10

1. ECUACIONES NO LINEALES 1.1 Método de Newton-Raphson 1.2 Método de la secante 1.3 Método de la bisección 1.4 Método de la regla falsa 1.5 Método de iteración del punto fijo 1.6 Problemas de aplicación a la Ingeniería Química 1.7 Aplicaciones en Ingeniería Química

11 11 14 16 19 22 24 32

2. SISTEMAS DE ECUACIONES NO LINEALES 2.1 Método de newton Raphson 2.1.1Aplicación 2.2 Problemas

36 36 37 40

3. INTERPOLACIÓN 3.1 Polinomios de interpolación de newton 3.1.1 Interpolación lineal 3.1.2 Interpolación cuadrática 3.1.3 Diferencias finitas divididas 3.2 Polinomios de interpolación de Lagrange 3.3 Problemas

44 45 46 47 47 48 50

4. REGRESIÓN 4.1 Regresión lineal 4.2 Regresión polinomial 4.3 Regresión lineal múltiple 4.4 Problemas

53 53 54 55 57

5. DIFERENCIACIÓN NUMÉRICA 5.1 Diferenciación mediante método Newton 5.2 Diferenciación de Lagrange: datos discretos

68 67 67

2

5.3

Problemas

69

6. INTEGRACIÓN NUMÉRICA 6.1 Método del trapezoide 6.2 Regla de Simpson 6.2.1 regla de Simpson 1/3 6.2.2Simphson 3/8 6.3 Problemas

76

7. SISTEMA DE ECUACIONES LINEALES 7.1 Métodos de Jacobi 7.2 Métodode Gauss – Seidel 7.3 Problemas

85 85 87 88

8. ECUACIONES DIFERENCIALES NUMÉRICAS 8.1 Método de Euler 8.2 Método de Euler modificado 8.3 Métodos de Runge-Kutta (rk) 8.4 Método de diferencias finitas 8.5 Sistema de ecuaciones diferenciales de primer orden 8.6 Problemas

91 92 94 96 100 101 103

9. ECUACIONES DIFERENCIALES PARCIALES (EDP)

110

9.1

Problemas

110

10.

APLICACIONES DE INGENIERÍA QUÍMICA EN

76

77 78 80 81

POLYMATH- MATHCAD Y MATLAB 10.1 Problemas con polymath 10.2 Problemas con mathcad 10.3 Problemas con matlab

112 112 121 131

SILABO

147

3

I.

RESUMEN

La presente Investigación tuvo como propósito la elaboración de un texto universitario titulado “TEXTO: " METODOS NUMERICOS UTILIZANDO POLYMATH, MATHCAD Y MATLAB APLICADOS A INGENIERIA QUIMICA” Se trata de un texto básico que se expone brevemente los fundamentos teóricos, ilustraciones

con problemas

resueltos

y propuestos de

ingeniería química como, equilibrio químico, ecuaciones de estado, transferencia de calor, cinética química y al final se plantea problemas resueltos haciendo uso de software de polymath, mathcad y matlab.. Asi mismo en cada capitulo se prantea problemas propuestos aplicados a la ingeniería química .

4

II.

INTRODUCCIÓN.

El tema de investigación, está referida a la elaboración de un texto Universitario, cuyo propósito es apoyar en la labor de formación de los alumnos, en el curso de métodos numéricos aplicada a ingeniería .Quimica en la universidad Nacional del Callao facultad de Ingeniería Química. Durante mi experiencia en la docencia, en el intento de encontrar textos necesarios para el dictado del curso de

métodos numéricos , se ha

podido constatar que existe poca bibliografía en nuestro medio haciendo uso de mathcad,

polymath y matlab en

un solo texto, mediante ste

trabajo se ha desarrrollado problemas aplicados a la ingeniería química haciendo uso de una técnica numérica, al mismo timpo se ha resuelto problemas con el software mathcad, polymath y mathcad, que se vienen usando con mayor intesnsidad en los últimos años a nivel de ingeniería. En este trabajo se ha resaltado el capitulo de ecuaciones no lineales ecuaciones diferenciales puesto

y

que gran parte de los modelos de

ingenieira química se encuentran en esta area.

5

2.1

A.

PLANTEAMIENTO DEL PROBLEMA DE INVESTIGACION

DESCRIPCIÓN Y ANÁLISIS DEL TEMA El presente trabajo de investigación es una propuesta para la

elaboración de un texto Universitario titulado: TEXTO:

“METODOS

NUMERICOS

UTILIZANDO

POLYMATH,

MATHCAD Y MATLAB APLICADOS A LA INGENIERIA QUIMICA.”

Dirigido a estudiantes de pre – grado de Ingeniería Química y otras especialidades afines, que presente de una manera didáctica los principios fundamentales y el uso adecuado de los software de polymath, mathcad y matlab, lo que permitirá cumplir con los propósitos de una adecuada enseñanza y formación profesional.

El texto contendrá una base de teoría apropiada y práctica que van a permitir desarrollar criterios y habilidades, que resultara muy valioso para los propósitos de este texto.

B.

PLANTEAMIENTO DEL PROBLEMA

¿Como

elaborar

un

TEXTO:

“METODOS

NUMERICOS

UTILIZANDO POLYMATH, MATHCAD Y MATLAB APLICADOS A LA INGENIERIA QUIMICA.”, que oriente adecuadamente a los estudiantes de Ingeniería Química?

2.2

2.2.1

OBJETIVO Y ALCANCE DE LA INVESTIGACION

OBJETIVO GENERAL

Desarrollar

un

TEXTO:

“METODOS

NUMERICOS

UTILIZANDO

POLYMATH, MATHCAD Y MATLAB APLICADOS A LA INGENIERIA QUIMICA.”

6

2.2.2.

1.

OBJETIVOS ESPECIFICOS

Recopilar Información básica y actualizada, necesaria para iniciar el desarrollo del texto.

2.

Analizar y procesar la información para iniciar el desarrollo del texto.

3.

Desarrollar los capítulos del texto referido a los fundamentos teóricos.

4.

Desarrollar los capítulos del texto referido a las aplicaciones prácticas aplicada la ingeniería química

5.

Desarrollar los capítulos del texto referido a las aplicaciones prácticas aplicada la ingeniería química haciendo uso del uso polymath , mathcad y matlab

2.3.

ALCANCE DE LA INVESTIGACION

El presente trabajo de investigación de acuerdo a la naturaleza del problema se puede

manifestar que es una investigación básica y

Aplicada, dado los modelos planteados para la solución numérica proviene de fenómenos químicos, cuya solución analítica sería demasiada compleja. El aporte del trabajo de investigación

estará orientado al sector

académico conformado por los profesores, estudiantes y egresados de la Facultad de ingeniería química de la Universidad Nacional del Callao y otras Universidades del país. Por otro lado, este texto también podría ser utilizado por estudiantes de especialidades afines tales como ingeniería de Alimentos, ingeniería ambiental e ingeniería Industrial.

2.4.

IMPORTANCIA Y JUSTIFICACION DE LA INVESTIGACION

2.4.1 IMPORTANCIA Al desarrollar el texto propuesto se facilitará el proceso de enseñanza – aprendizaje en la formación profesional del estudiante universitario

a nivel de pre-grado

y pueda facilitar los cálculos

laboriosos de ingeniería haciendo uso de software. 2.4.2 JUSTIFICACION 7

La contribución del presente trabajo estará orientada a la preparación y el entrenamiento de los alumnos de ingeniería química en el curso de Métodos numéricos, adquiriendo fundamentos teóricos y la parte práctica que consiste en efectuar cálculos de balances de materia, energía, termodinámica, reacciones química,

transferencia de calor

entre otros.

III.

MARCO TEÓRICO En la presente investigación se ha incorporado la teoría resumida y simplificada para nueve capítulos del presente texto. Así por ejemplo, en el capítulo I se hace referencia a los diferentes métodos de

solución de ecuaciones no lineales con

sus respectivos

ejemplos y problemas. En el capítulo

III , se hace referencia de las técnicas de

interpolación , resaltando el método de diferencias por newton . Asimismo en el capítulo VIII de ecuaciones diferenciales , se hace referencia de los diferentes ordenes de solución de Runge Kutta y con sus respectivos ejemplos. En el capítulo X se presenta problemas resueltos con polymath , mathcad y Matlab. IV.

MATERIALES Y MÉTODOS Materiales: 

Materiales De oficina



Material bibliográfico



Software Polymath



Software Mathcad



Software Matlab



Material de cómputo e impresión

Métodos La elaboración del texto, propósito de la investigación le demandó al suscrito, ordenar la información disponible y complicada en función del Syllabus propuesto del curso de métodos numéricos . 8

La estructuración del texto responde a la experiencia docente en la Facultad de Ingeniería Química de la Universidad Nacional del Callao.

Para la elaboración del texto, se tuvo cuidado en recurrir a la síntesis de los aspectos teóricos,

selección de los problemas

de

ingeneiria química, elaboración de los programas en matnab, asi como revisión de tutoriales de mathcad.

En cuando al planteamiento de problemas, se quiere plasmar la experiencia con el dictado del curso, que me ha permitido revisar una extensa bibliografía sobre la materia, así como volcar en el texto los resultados de mi ejercicio profesional en el campo de Ingeniería Química para satisfacer el propósito de la investigación.

V.

RESULTADOS El resultado de la presente investigación es la elaboración de un texto

universitario

titulado:

TEXTO:

“METODOS

NUMERICOS

UTILIZANDO POLYMATH, MATHCAD Y MATLAB APLICADOS A LA INGENIERIA QUIMICA.”, , el cual se adjunta al presente. El texto contiene 9 capítulos. La teoría desarrollada en el texto, responde a los aspectos básicos de métodos numericos. Los problemas resueltos en el texto, tienen el propósito de dar las pautas de la aplicación de la teoría desarrollada. Se ha logrado un texto base para el curso de métodos numéricos necesario

en la formación universitaria del estudiante de Ingeniería

Química. VI.

DISCUSIÓN El texto universitario titulado ““METODOS NUMERICOS UTILIZANDO POLYMATH, MATHCAD Y MATLAB APLICADOS A LA INGENIERIA QUIMICA.”, aplicados a la Ingeniería Química, que es el resultado de la investigación a que se refiere el presente informe, se caracteriza por presentar

la metodología de calculo de los modelos . Los problemas

resueltos y planteados han sido seleccionados con la intención 9

de

brindar

una mayor claridad

a

los alumnos

y puedan entender los

fundamentos teóricos tratados. Los textos de METODOS NUMERICOS

contienen demasiada

información, muchas veces muy detallado, sin conexión directa con aplicación inmediata. Por eso, el presente texto va a tratar de desenvolver el contenido de tal manera que cada capítulo describa en forma precisa y concreta la teoría y los problemas de aplicación. Sin embargo,

y por ende no

sustituye el uso de la bibliografía de la especialidad, a la que deberá referirse necesariamente quienes pretendan profundizar en conocimientos de temas específicos.

VII. REFERENCIALES

Existen textos que se ocupan de los métodos numéricos en general entre los cuales tenemos 1. A. Constantinides and N. Mosotoufi, Numerical Methods for Chemical Engineers with MATLAB Applications Prentice Hall , Upper Saddle River, 1999. 2. Burden, R. Y Faires J. Análisis Numérico. Edit. Iberoamericana, México, 1985. 3. Carnahan, B. Luther, A. Wilkes Cálculo Numérico, Aplicaciones Editorial Rueda, Madrid, 1979 4. Carrasco Venegas Luis . Métodos Numéricos aplicados a la ingeniera 5. Nieves, A., Domínguez, F. Métodos Numéricos Aplicados a la Ingeniería Química Edit. CECSA, México 1985. 6. Nakamura, S. Métodos Numéricos aplicados con

software. Edit.

Prentice – Hall Hispano Americano, S.A. México, 1992. 7. S.C. Chapra y R. P. Canale,. Métodos Numéricos para Ingenieros. McGraw Hill, México,1999

10

1.ECUACIONES NO LINEALES Los métodos numéricos de resolución de ecuaciones no lineales suelen ser métodos iterativos que producen una sucesión de valores aproximados de la solución, que se espera, que converja a la raíz de la ecuación. Estos métodos van calculando las sucesivas aproximaciones en base a los anteriores, a partir de una o varias aproximaciones iníciales.. Para saber que método debemos aplicar, hay que tener en cuenta la capacidad de separar raíces cercanas, confiabilidad en el alcance de soluciones evitando errores numéricos graves y orden de convergencia. Uno de los problemas que con mayor frecuencia aparece en la ciencia y en la ingeniería es hallar las raíces de una ecuación no lineal de la forma f(x) = 0. Estudiaremos Métodos Iterativos para determinar aproximaciones a raíces reales simples de la ecuación no lineal f(x) = 0.

1.1

MÉTODO DE NEWTON-RAPHSON

Este método, es uno de los más usados, a diferencia de los otros métodos, el método de Newton-Raphson no trabaja sobre un intervalo sino que basa su fórmula en un proceso iterativo. Supongamos que tenemos la aproximación xi a la raíz xr de f ( x ) ,

Figura Nº 1: Demostración método de newton

11

Trazamos la recta tangente a la curva en el punto  xi , f ( xi )  ; ésta cruza al eje x en un punto xi 1 que será nuestra siguiente aproximación a la raíz xr . Para calcular el punto xi 1 , calculamos primero la ecuación de la recta tangente. Sabemos que tiene pendiente m  f   xi 

Y por lo tanto la ecuación de la recta tangente es: y  f  xi   f   xi  x  xi 

Hacemos y  0 :  f  xi   f   xi  x  xi 

Y despejamos x :

x  xi 

f  xi   x  xi 

Que es la fórmula iterativa de Newton-Raphson para calcular la siguiente aproximación:

xi 1  xi 

f  xi  , si f   xi   0 f   xi 

Note que el método de Newton-Raphson no trabaja con intervalos donde nos asegure que encontraremos la raíz, y de hecho no tenemos ninguna garantía de que nos aproximaremos a dicha raíz. Desde luego, existen ejemplos donde este método no converge a la raíz, en cuyo caso se dice que el método diverge. Sin embargo, en los casos donde si converge a la raíz lo hace con una rapidez impresionante, por lo cual es uno de los métodos preferidos por excelencia. También observe que en el caso de que f   xi   0 , el método no se puede aplicar. De hecho, vemos geométricamente que esto significa que la recta tangente es horizontal y por lo tanto no intersecta al eje en ningún punto, a menos que coincida con éste, en cuyo caso xi mismo es una raíz de f ( xi )! Ejemplo Usar el método de Newton-Raphson, para aproximar la raíz de

f ( x)  e x  ln x , comenzando con x0  1 y hasta que a  1% .

12

Solución En este caso, tenemos que f ( x )   e  x 

1 x

De aquí tenemos que:

Iniciamos con x0  1 y obtenemos:

En este caso, el error aproximado es,

Continuamos el proceso hasta reducir el error aproximado hasta donde se pidió. Resumimos los resultados en la siguiente tabla: Tabla Nº1 Aproximación de la raíz y porcentaje de error Aprox. a la raíz

Error aprox.

1 1.268941421

21.19%

1.309108403

3.06%

1.309799389

0.052%

De lo cual concluimos que la aproximación obtenida es:

13

1.2

MÉTODO DE LA SECANTE

Uno de los objetivos de este método es eliminar el problema de la derivada de la función, ya que existen funciones que describen fenómenos físicos y químicos, cuya derivada es muy compleja El método de la secante es muy similar al de Newton con la diferencia principal que en este método de la secante no requiere de la derivada. Este método funciona por medio de dos puntos sobre el eje x, es decir un intervalo (xi-1,xi), los cuales se evalúan en la función para sacar los puntos correspondientes en el eje de la y, los puntos a obtener son f(x i-1) y f(xi), por lo que las coordenadas de los puntos que interceptan a la función son (x i1,f(xi-1))

y el (xi ,f(xi)).

Se debe considerar que los puntos xi-1 y xi deben de contener a la raíz, por lo que el punto xi-1 debe estar a la izquierda y el punto xi a la derecha de la raíz. Estos dos puntos que interceptan a la función, se unen por medio de una recta, la cual al cruzar el eje de la x, genera el siguiente punto de acercamiento xi+1 , el cual quedara ubicado entre el intervalo propuesto, como se muestra en la Figura. El método de la secante parte de dos puntos (y no sólo uno como el método de Newton) y estima la tangente (es decir, la pendiente de la recta) por una aproximación de acuerdo con la expresión: Este método se basa en la fórmula de Newton-Raphson, pero evita el cálculo de la derivada usando la siguiente aproximación:

f   xi  

f  xi1   f  xi  xi1  xi

(Recuérdese la solución numérica al problema del cuerpo en caída libre). Sustituyendo en la fórmula de Newton-Raphson, obtenemos:

xi1  xi 

f  xi  f  xi   xi  f  xi1   f  xi  f   xi  xi1  xi

14

Ejemplo Usar el método de la secante para aproximar la raíz de f  x   e x  x , 2

comenzando con x 0  0 , x 0  1 y hasta que a  1% . Solución Tenemos que f  x0   1 y f  x1   0, 632120558 , que sustituimos en la fórmula de la secante para calcular la aproximación x 2 :

Con un error aproximado de:

Como todavía no se logra el objetivo, continuamos con el proceso. Resumimos los resultados en la siguiente tabla: Tabla Nº2 Aproximación de la raíz y porcentaje de error

Aprox. a la raíz

Error aprox.

0 1

100%

0.612699837

63.2%

0.653442133

6.23%

0.652917265

0.08%

De lo cual concluimos que la aproximación a la raíz es:

x4  0,652917265

15

1.3

MÉTODO DE LA BISECCIÓN

Este método es basado en el teorema de Bolzano, que establece que si una función continua cambia de signo en el intervalo (a,b), es decir, f(a)f(b) 0 se llama fórmula de diferencia progresiva, Fórmulas de diferencias divididas hacia adelante Sea los puntos x

Xi

Xi+1

Xi+2

XI+3

…….

f

f (Xi)

f (Xi+1)

f (Xi+2 )

f (XI+3 )

…….

5.1.DIFERENCIACIÓN MEDIANTE MÉTODO NEWTON Primera derivada

f '  xi  

f  xi 1   f  xi  h

, f '  xi  

 f  xi  2   4 f  xi 1   3 f  xi  2h

Segunda derivada

f ''  xi   f ''  xi  

f  xi 2   2 f  xi 1   f  xi  h2  f  xi 3   4 f  xi  2   5 f  xi 1   2 f  xi  h2

5.2. DIFERENCIACIÓN DE LAGRANGE: datos discretos Construimos el polinomio de Lagrange 68

L  x   L1  x  y1  L2  x  y2  L3  x  y3 

 x  x2  x  x3   x  x1  x  x3   x  x1  x  x2  y1  y2  y  x1  x2  x1  x3   x2  x1  x2  x3   x3  x2  x3  x1  3

Diferenciando la función f  x  L  x  

2 x  x2  x3 y  x1  x2  x1  x3  1

2 x  x1  x3 2 x  x1  x2 y2  y  x2  x1  x2  x3   x3  x2  x3  x1  3

Asumiendo espacio constante f  x 

2 x  x2  x3 2 x  x1  x3 2 x  x1  x2 y1  y2  y3 2 2 2 x x 2 x 2

f  x 

2 x  x2  x3 2 x  x1  x3 2 x  x1  x2 y1  y2  y3 2 2 2 x x 2 x 2

f   x1  

2 x1  x2  x3 2x  x  x 3 y1  4 y2  y3 2x  x  x y1  1 1 2 3 y2  1 1 2 2 y3  2 2 x 2 x x 2 x

f   x2  

2 x2  x2  x3 2 x  x1  x3 y  y1 2 x  x1  x2 y1  2 y2  2 y3  3 2 2 2 2 x 2 x x 2 x

f   x3  

2 x3  x2  x3 2x  x  x 2x  x  x y  4 y2  3 y3 y1  3 1 2 3 y2  3 1 2 2 y3  1 2 2 x 2 x x 2 x

f  x 

2 x  x2  x3 2 x  x1  x3 2 x  x1  x2 y1  y2  y3 2 2 2 x x 2 x 2

f   x  

y  2 y2  y3 1 2 1 y  y2  2 y3  1 2 1 2 x x x x 2

segunda derivada

5.3 PROBLEMAS 1. El flujo de calor en la interfaz suelo-aire puede calcularse con ley dT q  z  0   k  C dz z 0

69

Z (m) T (º C)

0 13,5

1,25 12

1,75 10

Donde q = flujo de calor, k = coeficiente de difusividad térmica (3.5x10 -7),  = la densidad del suelo (1800Kg/m3), C = calor específico del suelo (840).

f '  0   13.5

2  0   1.25  3.75

 0  1.25 0  3.75

 12

2  0   0  3.75

1.25  01.25  3.75

 10

2  0   0  1.25

 3.75  0 3.75  1.25

 1,333

q = -3.5x10-7 . 1800 (-1,333). =70.56 2. A 400ºC y con concentraciones: CO   NO2   0,10 mol/L se obtuvieron los siguientes datos para la reacción dada:

COg   NO2

g  

CO2

g   NOg 

Tiempo (s)

0

10

20

30

CO  mol L

0,1

0,067

0,05

0,04

Determinar el orden de la reacción dada Determinar la constante de velocidad Calcular la velocidad instantánea de la reacción en el instante t = 10 s Solución: De la ecuación: vA  -

dCA  k  CA  dt

n

 dC  Tenemos: ln  - A   ln x  n ln  CA  ecuación lineal  dt  Diferenciación numérica con 4 puntos: x0

x1

x2

x3

Tiempo (s)

0

10

20

30

CO  mol L

0,1

0,067

0,05

0,04

y0

y1

y2

y3

h  10 70

dCA 1   11y0  18 y1  9 y2  2 y3  dt 6h

Para x0:



dCA dt



1 11 0,1  18  0,067  9  0,05  2  0,04  4,4 103 6  10

dCA 1   2 y0  3y1  6 y2  y3  dt 6h

Para x1:



dCA dt



1 2  0,1  3  0,067  6  0,05  0,04  2,35 103 6 10

dCA 1   - y0  6 y1  3y2  2 y3  dt 6h

Para x2:



dCA dt



1 - 0,1  6  0,067  3  0,05  2  0,04  1,2 103 6  10

dCA 1   2y0  9 y1  18 y2  11y3  dt 6h

Para x3:



dCA dt



1 2  0,1  9  0,067  18  0,05  11 0,04  9,5 104 6  10

Llenamos la tabla y procedemos a graficar: TABLA 11 CALCULO DIFERENCIACIÓN

CO  mol L

- dCA

0

0,1

10

 

 dC A  ln dt  

ln CA

4,4 x 10-3

-5,426

-2,302

0,067

2,35 x 10-3

-6,053

-2,703

20

0,05

1,2 x 10-3

-6,725

-2,995

30

0,04

9,5 x 10-4

-6,959

-3,218

Tiempo (s)

dt

71

Del ajuste lineal tenemos que la ecuación es: y  1, 7 x  1, 41 Donde: Orden de reacción es: n  pendiente  1, 7  2 Constante de velocidad es: ln K  -1,41 Para t=10 s: v  k  C At 10 s 



K  0,244

n

v  0,244   0,067 2  1,095  103

PROBLEMAS 1. La primera ley de difusión de Fick establece que el flujo másico Flujo másico = F  D

dc dx

Donde flujo másico es cantidad de masa que pasa a través de una unidad de area por unidad de tiempo ( g/ cm 2/s. D= coeficiente de difusión ( cm2/s). C: concentración X: distancia en cm. Un ing. químico mide la concentración de un contaminante en los sedimentos depositados en un lago obteniéndose los siguientes datos. X (cm)

0

1

2

c. 10-6 g/cm3

0,1

0,4

0,9

72

Use la mejor técnica de derivación numérica para estimar en x =1,5 y calcule el flujo másico del contaminante cuando D= 2x10 -6cm2/s, para un lago de 3x106m2 de sedimentos ¿cuánto contaminante podría ser transportada hacia el lago en un año? 2. Para la descomposición del pentaóxido de dinitrógeno a 45ºC se obtuvieron los siguientes datos:

2N 2 O5 ( g )  O2 ( g )  4NO2 ( g ) Tiempo(s) 0

 N 2 O5 

200

400

600

800

mol 2,5 2,22 1,96 1,73 1,53 L

determina de qué orden será la reacción 3. A 440ºC con y con concentraciones |CO| = |NO2| = 0,10 mol/L se obtuvieron los siguientes datos para la reacción dada:

CO( g )  NO2 ( g )  CO2 ( g )  NO( g ) Tiempo(s)

CO 

mol L

0

10

20

30

40

100

1000

0,1 0,067 0,05 0,04 0,033 0,017 0,002

a) Determina el orden de la reacción dada b) Determine la constante de velocidad c) Calcula la velocidad instantánea de la reacción en el instante t= 10 s.

4. Se ha estudiado la descomposición térmica de arsenamina sobre vidrio, 2 AsH3(g) 2 As(s) + 3H2(g) La presión total del sistema varia con el tiempo, a 350C, según se indica a continuación: t/h

0

4.33

16.00 25.50 37.66 44.75

Pt/cmHg 39.2 40.30 43.65 45.35 48.05 48.80

a) Determine el orden de reacción con respecto a la arsenamina. b) Calcule la constante de velocidad. 73

5. Para poder estudiar el comportamiento físico químico de una solución de bromo en presencia de luz solar, se coloca una pequeña cantidad de solución de bromo en un matraz y se expone al sol, obteniéndose los siguientes resultados: Tiempo( min.)

10

20

30

40

50

60

Ppm Br2 a 25ºC 2,45 1,74 1,23 0,88 0,62 0,44

En base a esta información calcular la constante de velocidad en las unidades adecuadas y orden de la reacción.

6. En el estudio experimental de la velocidad de reacción para la reacción 2 NO (g) + 2 H2(g) → N2(g) + 2 H2O (g) se han obtenido los siguientes datos a 298 K: [NO] mol/L 0,02 0,02 0,03 0,06 Calcular:

[H2] mol/L 0,15 0,45 0,50 0,50

Velocidad mol/l s 1,2 x 10-8 3,6 x 10-8 9,0 x 10-8 36,0 x 10-8

a) La ecuación de velocidad para esta reacción. b) El orden de reacción c) La constante de velocidad

7. Determinar el orden de reacción para la descomposición en fase gaseosa del peróxido de di-tert-butilo (CH3)3COOC(CH3)  C2H6 + 2CH3COCH3 La reacción se lleva a cabo en el laboratorio en un reactor isotérmico batch La presión total de reacción a lo largo de reacción evoluciona según la siguiente tabla

74

tiempo(min) Presión total (MmHg) 0.0 7.5 2.5 10.5 5.0 12.5 10.0 15.8 15.0 17.9 20.0 19.4

8. Determinar el orden de reacción: CH3-Cl(g) + H2O(g)  CH3-OH(g) + HCl(g) usando los datos de la tabla. v  k  [CH3 -Cl]n  [H2 O]m

Experiencia [CH3-Cl] (mol/l) [H2O] (mol/l) v (mol·l–1·s–1) 1

0,25

0,25

2,83

2

0,50

0,25

5,67

3

0,25

0,5

11,35

9. La transferencia de calor está dada por la siguiente ecuación dT q   kA dy donde

J   k = conductividad térmica   smK 

 

A = área m2 T = temperatura K  y = distancia m  k  0.025

J smK

A  3 m2

La temperatura en función de distancia

T  1493y3  2200y 2  1076y  500 La transferencia de calor para y= 1,2m

75

6. INTEGRACIÓN NUMÉRICA En esta lección comenzamos el estudio de métodos numéricos para el cálculo numérico de integrales de la forma.  f  

b

 f  x  dx a

Un método común para aproximar I(f) es reemplazando f(x) con un polinomio de interpolación. Este procedimiento se conoce como las reglas de Cuadratura de Newton. Examinamos los primeros dos casos de este método donde se usan polinomios de interpolación lineales y cuadráticos.

6.1

Método del trapezoide: Sea p1(x) el polinomio lineal que interpola a f(x) en x=a y x=b

f(x) L(x)

x0

x1

x

Figura Nº9: Integración método de trapecio usando dos puntos

Usando la fórmula para el área de un trapezoide o integrando p1(x) directamente se obtiene que



b



b

a

a

f ( x)dx 

h  f ( x0 )  f ( x1 ) 2

f ( x ) dx 

h  f (a )  f (b) 2

Generalizando Generalizando el método de trapecio

76

f(x )

h 

b a n x0

h

x1

h

x2

h

x3

h

x4

x

Figura Nº10: Integración método de trapecio usando n puntos



b

a

x1

x2

xn

x0

x1

x n 1

f(x)dx   f(x)dx   f(x)dx    

f(x)dx

h h h  f(x 0 )  f(x1 )   f(x1 )  f(x 2 )     f(x n 1 )  f(x n ) 2 2 2 h   f(x 0 )  2f(x1 )    2f(x i )    2 f ( xn 1 )  f ( xn )  2 b h a f(x)dx  2 f(x 0 )  2f(x1 )    2f(x i )    2 f ( xn1 )  f ( xn ) 

4

Evaluar la integral :  xe 2 x dx método de trapecio 0

4

I   xe 2 x dx  0

40  f (0)  f (4)  2(0  4e8 )  23847.66 2

6.2 REGLA DE SIMPSON Además de aplicar la regla trapezoidal con segmentos cada vez más finos, otra manera de obtener una estimación más exacta de una integral, es la de usar polinomios de orden superior para conectar los puntos. Por ejemplo, si hay un punto medio extra entre f(a) y f(b), entonces los tres puntos se pueden conectar con un polinomio de tercer orden. A las fórmulas resultantes de calcular la integral bajo estos polinomios se les llaman Reglas de Simpson.

77

6.2.1 REGLA DE SIMPSON 1/3 La Regla de Simpson de 1/3 proporciona una aproximación más precisa, ya que consiste en conectar grupos sucesivos de tres puntos sobre la curva mediante parábolas de segundo grado, y sumar las áreas bajo las parábolas para obtener el área aproximada bajo la curva. Por ejemplo, el área contenida en dos fajas, bajo la curva f(X) en la Figura. 2, se aproxima mediante el área sombreada bajo una parábola que pasa por los tres puntos:

f(x)

x0

h

L(x)

x1

h

x2

x

Figura Nº11: Integración método de Simpson usando tres puntos



b a

f ( x)dx 

2

c i0



i

f ( x i )  c 0 f ( x 0 )  c1 f ( x1 )  c 2 f ( x 2 )

h  f ( x 0 )  4 f ( x1 )  f ( x 2 )  3

Generalizando

78

h

ba n f(x)

…... x0 h x1 h x2 h x3 h

x4

xn-2 xn-1

xn

x

Figura Nº12: Integración método de Simpson usando n puntos



b

a

x2

x4

xn

x0

x2

xn 2

f(x)dx   f(x)dx   f(x)dx    

f(x)dx

h h f(x 0 )  4f(x1 )  f(x 2 )  f(x 2 )  4f(x 3 )  f(x 4 ) 3 3 h     f(x n  2 )  4f(x n 1 )  f(x n )  3 h   f(x 0 )  4f(x1 )  2f(x 2 )  4f(x 3 )  2f(x 4 )   3 4f(x 2i-1 )  2 f ( x2i )  4f(x 2i 1 )   

2 f ( xn  2 )  4 f ( xn 1 )  f ( xn )  Evaluar la integral



4

0

xe 2 x dx 4

I   xe 2 x dx  0



h  f (0)  4 f (2)  f (4) 3

2  0  4(2e 4 )  4e8   8240.411 3

Simpson’s 3/8 Approximate by a cubic polynomial La derivación de la Regla de los Tres Octavos de Simpson es similar a la regla de un tercio, excepto que se determina el área bajo una parábola de tercer grado que conecta 4 puntos sobre una curva dada. La forma general de la parábola de tercer grado es: 79

L(x)

x0

h

f(x)

x1

h

x2

h

x3

x

Figura Nº13: Integración método de Simpson usando cuatro puntos



b

a

3

f ( x) dx   ci f ( xi )  c 0 f(x 0 )  c1f(x1 )  c 2 f(x 2 )  c3 f(x 3 ) i 0



h

3h  f ( x0 )  3 f ( x1 )  3 f ( x2 )  f ( x3 ) 8

b-a 3

6.2.2 SIMPHSON 3/8 4

I   xe 2 x dx  0

3h  4 8  f (0)  3 f ( )  3 f ( )  f (4)  8  3 3

3(4/3) 0  3(19.18922)  3(552.33933)  11923.832  6819.209 8 5216.926  6819.209   30.71% 5216.926 

80

6.3 PROBLEMAS 1.

Los datos que aparecen en la tabla siguiente corresponden a una velocidad transversal de salida en una tubería ¿cuál es el flujo volumétrico en m3/s?

r (cm)

0

5

10

15

V(m/s) 50 49,5 49 48

20

25

30

35

40

45 47 50

46,5

45

43 40,5 37,5 34 25 0

Q  2  V * rdr

2.

La cantidad de calor necesaria para elevar la temperatura de 1000gr de H2O desde –1000C hasta 2000C se puede evaluar de la siguiente forma:  H



 m

T 2 T 1

Cp

 T  dT

Cp  2.64 107 T 2  1.56 104 T  0.132 Determinar  empleando la técnica del trapecio y la regla de Simpson

3.

La masa que entra

o sale de un reactor perfectamente agitado en

un periodo especifico se puede determinar mediante: t2

M   Q.cdt ti

Donde t1 y t2 son el tiempo inicial y final respectivamente Q es el flujo volumétrico y c es la concentración.

81

4.

T(,min )

C,mg/m3

0

10

5

22

10

35

15

47

20

55

25

58

30

52

35

40

40

37

45

32

50

34

Los datos que se muestra a continuación corresponde a la reacción de tipo. A  B en un reactor tubular X

(-1/rA)

0

0.0053

0.1

0.0052

0.2

0.005

0.3

0.0045

0.4

0.0033

0.5

0.0025

0.7

0.0018

0.8

0.00125

0.85

0.001

La ecuación de diseño para este tipo de reactor es de: V  FA0  

x

0

1 dx  rA

V:volumen(m3)FA0: flujo de alimentación A ,X: conversión el flujo molar de alimentación es de 2 mol/ min. el 80 % es de A y la diferencia de inertes Determine el volumen del reactor para una conversión de 90% de A 82

5.

Se tiene los siguientes puntos con h constante x

X0

X1

X2

X3

X4

y

Y0

Y1

Y2

Y3

Y4

Demostrar una fórmula de integración

6.

La reacción de saponificación de acetato de etilo en medio alcalino se conduce en fase homogénea de acuerdo con la siguiente ecuación. CH3COOCH2CH3

+ NaOH 

CH3COONa + CH3CH2OH

Donde la velocidad de reacción es de

ra 

5817 dCA  kCACB , k  2, 028 *109 * e T dt

de reacción

Constante de velocidad

el tiempo de reacción para un reactor discontinuo es

calculado mediante la siguiente reacción. CA

t

1 dCA ( ra ) CA0



Donde : CA = CA0(1-XA) ; CB = CB0 – CA0XA Si la concentración inicial de A y B es de 0,8 mol-g/ L ¿ Cual es el tiempo para alcanzar una conversión de 95%? A 500K

7.

La cantidad de calor necesaria para elevar la temperatura de 1 000 g de H2O

desde 100C hasta 900C se puede evaluar de la siguiente

forma:

H  m  Cp T dT T2

T1

Donde:

Cp  2.64 107 T 2  1.56 104 T  0.132 Determinar  empleando la técnica del trapecio y la regla de Simpson 8.

Resolver la siguiente integración.

83

9.

En un reactor tubular continuo se realiza la eliminación de un contaminante en fase gas mediante la reacción: A + B  R + S, a 200ºC y 7 atm. Para ello, se alimenta una corriente con A y B en relación equimolar, en las mismas condiciones de presión y temperatura, con un caudal de 500 L/min. Si se quiere alcanzar una conversión del 85%, ¿qué volumen de reactor es necesario? Datos: (rA) = 10.CA.CB [mol/L.min]

10.

Una reacción A  R, de ecuación cinética r (mol/L.h) = 1,5 CA (a 50ºC), se lleva a cabo en un reactor discontinuo, de forma isoterma, siendo la temperatura de trabajo constante e igual a 50ºC. La mezcla reaccionante tiene una concentración inicial de 10 mol/L en A ¿Cuál es el tiempo de reacción necesario para alcanzar una XA de 0,90?

84

7. SISTEMA DE ECUACIONES LINEALES

INTRODUCCIÓN

La solución de sistemas de ecuaciones lineales es un tema de gran utilidad en diversas ramas del conocimiento como la economía, la biología, física, ingenierías, etc. La resolución de sistemas de cualquier número de ecuaciones (10, 50, 100, 500, etc.) es una realidad hoy en día, gracias a los computadores, lo cual proporciona solución directa. Un gran número de problemas prácticos de ingeniería se reduce a resolver un sistema de ecuaciones.

Un sistema de

ecuaciones lineales con

incógnitas, tiene la forma:

a11x1 + a12x2 + ….. + a1nxn = b1 a21x1 + a22x2 + ….. + a2nxn = b2  am1x1 + am2x2+ …. + amnxn = bm Con notación matricial se escribe así: a11 a  21   a m1

a12 ..... a1n  a 22 ......a 2 n    a m 2 .....a mn 

 x1  b1  x     2   b2           xm  bm 

 AX=B Donde

A es la matriz de coeficientes del sistema X es el vector incógnita B es el vector de términos independientes.

7.1 METODOS DE JACOBI (Método de desplazamiento simultáneos)

a11 x1 + a12 x2 + a13 x3 = b1 a21 x1 + a22 x2 + a23 x3 = b2 a31 x1 + a32 x2 + a33 x3 = b3 85

cona11, a22 y a33 distintos de cero. Se despeja x1de la primera ecuación,

de la segunda y

de la tercera con lo

que se obtiene. =−

=− que en notación matricial queda

=

+

+



+



=−

⎡ 0 ⎢ = ⎢− ⎢ ⎢− ⎣



− −

+

⎡ ⎢ ⎢ +⎢ ⎢ ⎢ ⎣



⎤ ⎥ ⎥ − ⎥ 0 ⎥ ⎦

0

⎡ 0 ⎢ = ⎢− ⎢ ⎢− ⎣

− −

0



⎤ ⎥ ⎥ − ⎥ 0 ⎥ ⎦

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ =⎢ ⎢ ⎢ ⎣

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

Una vez que se tiene la forma en notación matricial, se propone un vector inicial

( )

que puede ser

solución .

( )

= 0, o algún otro que sea aproximado al vector

Para iterar existe dos variantes. Si

( )

=

es el vector aproximación a la solución

después de

tiene para la siguiente aproximación.

86

iteraciones, entonces se

(

)

1 ( ⎡ ⎢ ⎢ 1 ( =⎢ ⎢ ⎢ 1 ( ⎣

=











)⎤ ⎥ ⎥ )⎥ ⎥ )⎥ ⎦



7.2 METODO DE GAUSS - SEIDEL (método de desplazamientos sucesivos) En este método los valores que se van calculando en la (

+ 1) −ésima

iteración se emplean para calcular los valores faltantes de esa misma iteración; es decir, con

( )

se calcula

(

)

de acuerdo con:

( ⎡ ⎢ 1 ( =⎢ ⎢ 1 ⎢ ( ⎣











87



)

⎤ )⎥ ⎥ ⎥ )⎥ ⎦

7.3 PROBLEMAS 1) Determine los valores de x1 , x2 y x3 usando el método de gaus seidel

2 x1  7 x2  11x3  6 x1  2 x2  x3   5 7 x1  5x2  2 x3  17 2) Se tiene la siguiente torre de separación realizar los balances y calcular los flujos masicos de cada corriente. H2O

0,2% de HCl 99.8% aire

Masa = 2640Kg/h

43,1% de HCl

63,8% H2O

56,9% aire

36,2% de HCl

Figura Nº14:Diagrama de una columna de absorción 3) Un proceso de 5 etapas en equilibrio para una extracción liquido

o

absorción gaseosa puede ser modelado mediante un sistema

de

ecuaciones lineales de acuerdo como se muestra a continuación.

(1  rr ) X 1  rrX 2   F0

X i1  (1  rr) X i  rrXi1  0 para i = 2,3, ………..(n-1) X n 1  (1  rr) X n  F Para : rr = 0,9 ; F0 = 0,05 y F = 0,5 Resolver el sistema de ecuaciones lineales

88

4) Para el sgt.sistema de separación, se conoce el flujo masico(Kg/hr) y las fracciones de masa de cada sustancia de entrada de la columna de separacion1 (F0) y la salida de la columna 2(F3 ,F1Y F4)¿ Cual es el valor de F2 

F1

W 1 = 0.2 W 2 = 0.6 W 3 =0.2

W 1 = 0.2 W 2 = 0.6 W 3 =0.2

F0

W 1 = 0.2 W 2 = 0.6 W 3 =0.2

F3

F2

F4

W 1 = 0.2 W 2 = 0.6 W 3 =0.2

5) Resolver la siguiente sistema de ecuación

4 x1  x 2   1  8 x1  x 2  x 3  1 3 3 x2  2 x3  4 x4   3 x 3  x 4  x 5  2 ,1 2 x 4  6 x 5  3, 4 6) La resistencia de un termistor varia con la temperatura 1 2 3  a 0  a1 ln( R )  a 2 ln R   a 3 ln R  T

T es la temperatura en Kelvin, R es resistencia en ohms, y a0 , a1 , a2 , a3 son constantes R

T

ohm

C

1101.0

25.113

911.3

30.131

636.0

40.120

451.1

50.128

Determine las constantes 89

7) En la figura se muestra la disposición de 5 tanques donde se lleva a cabo una mezcla de corrientes de entrada y salida. En base a la información determine concentración en cada en estado estacionario

2m3/min. C5 3

3m /min. 2m3/min. 5m3/min. 10mg/m3

C1 C1

C2

C4

3

3m /min.

3

1m /min.

1m3/min.

C3 1m3/min.

8m3/min.

8m3/min.20mg/m3 Figura Nº15: Sistema de reactores

90

11m3/min.

8. ECUACIONES DIFERENCIALES NUMÉRICAS

INTRODUCCIÓN Con mucha frecuencia aparecen problemas en ingeniería, física, química, ecología meteorología, sociología, etc., que exigen el manejo de ecuaciones diferenciales, muchas de las cuales no se pueden resolver pos los métodos convencionales, teniéndose que recurrir a métodos numéricos aproximados. Se llama ecuación diferencial aquella ecuación que contiene una variable dependiente y sus derivadas con respecto a una o más variables independientes. Se dividen en dos grandes grupos: Ordinarias, si contienen una sola variable independiente y Parciales, cuando contienen varias variables independientes. Estudiaremos únicamente las ecuaciones diferenciales ordinarias (EDO). Estas se clasifican y estudian según el orden de la mayor derivada que aparece en la respectiva ecuación diferencial. Ejemplo:

dy  ky dt

m

d2 y  ky dt 2

es de 1er orden es de 2º orden

d3 y dy  5  6y  0 es de 3er orden 3 dx dx Los problemas que encierran el uso de ecuaciones diferenciales ordinarias, constarán de: Una ecuación diferencial ordinaria: dy  y1  f(x, y) dx

Las condiciones iniciales. El intervalo para la variable independiente x, para el cual debe determinarse la función y, si existe.

91

8.1.

MÉTODO DE EULER

Consiste en dividir el intervalo de X0 a Xn en n subintervalos de ancho h.

h

xn - xo obteniéndose un conjunto discreto de (n + 1) punto: x0, x1, x2 …… n

xn en el intervalo [x0 , xn]

generándose la sucesión de aproximaciones

siguientes: y1 = y0 + hf (x0 , y0) y2 = y1 + hf (x1 , y1) y3 = y2 + hf (x2 , y2) . yn+1 = yn + hf(xn , yn) con i = 0 hasta n siendo f(xi , yi) la ecuación diferencial evaluada en xi y yi. yi+1 = yi + hf(xi, yi)

Ejemplo 1: Utilice el método de Euler para integrar numéricamente la ecuación: y’=f(x,y) = -2x3 + 12x2 – 20x + 8.5 de x=1.0 hasta x=3.0 considerando 4 intervalos y el valor inicial en x=1 es y=3 x0 = 1

h

f(x0 , y0) = f(1 , 3) 1.0

1.5

x0

x1

2.0 x2

2.5 x3

3.0 x4x4

3 1  0.5 4

yi+1 = yi + hf(xi , yi)

i

xi

yi

0 1

3

F(xi , yi) yi+1=yi+hf(xi ,yi) -1.5

2.25

1 1.5 2.25

-1.25

1.625

2 2.0 1.625

0.50

1.875

3 2.5 1.875

2.25

3.000

4 3.0 3.000

2.5

4.25

92

Ejemplo 2. Dada la ecuación diferencial

dy  f(x, y) = x - y , resuelva por el dx

método de Euler para el intervalo [0, 1] considerando 5 intervalos y siendo y(0)=2 1 0 = 0.2 5

h

Y(1)=?

yi+1 = yi + hf(xi , yi)

i 0

xi 0

x0 0

x1

x2

x3

0.2

0.4

0.6

yi

x4

x5

0.8 00.800.80

F(xi , yi) yi+1=yi+hf(xi ,yi)

2

-2

1.6

1 0.2

1.6

-1.4

1.32

2 0.4

1.32

-0.92

1.136

3 0.6

1.136

-0.536

1.0288

4 0.8

1.0288

-0.2288

0.98304

5 1.0

0.9834

-0.01696

0.986432

El valor exacto es 1.10364 E R 

1.10364  0.986432  100%  10.93% 1.10364

APLICACIÓN Un balón

de acero

de 1200 K es enfriado con corriente de aire

a

temperatura ambiente de 300K. Asumiendo la transferencia de calor en forma de radiación es mediante la ecuación diferencial . dT  2,2067  10 12 T 4  81 108 , T 0   1200 K dt Determine la temperatura para t= 480s usando el método de euler dT  2, 2067  1012 T 4  81  108  dt f  t , T    2, 2067  10 12  T 4  81  108 

T1  T0  f  t0 , T0  h Ti 1  Ti  f  ti , Ti  h  1200  f  0,1200 240





 1200   2.2067  10 12 1200 4  81  10 8  240 93

 1200   4.5579 240

T1  106.09K  0  240

t  t1  t0  h  240 T1  T  240  106,09K

2  1  f t1 ,1  h  106.09  f  240,106.09 240





 106.09  2.2067  10 12 106.09 4  81  108  240

 106.09   0.017595 240  110.32K

t  t2  t1  h  240  240  480

2    480  110.32K

8.2.

MÉTODO DE EULER MODIFICADO

Una fuente fundamental de error en el método de Euler es que la derivada al principio del intervalo se supone que se aplica a través del intervalo entero y para obtener una exactitud razonable se utiliza un intervalo muy pequeño a cambio de un error de redondeo mayor. El método de Euler modificado trata de evitar este problema utilizando un valor promedio de la derivada tomada en los dos extremos del intervalo en lugar de la derivada tomada en un solo extremo. El método de Euler modificado consta de dos pasos básicos:

1.

Se parte de (x0,y0) y se utiliza el método de Euler a fin de calcular el valor de y correspondiente a x1. Este valor de y se denotará como y1 ya que solo es un valor transitorio para y1. Esta parte del proceso se conoce como paso predictor.

94

2.

El segundo paso se llama corrector, pues trata de corregir la predicción. En el nuevo punto obtenido  x1 , y1  se evalúa la derivada de f  x1 , y1  usando la ecuación ordinaria que se está resolviendo. Se

obtiene la media aritmética de esta derivada y la derivada en el punto inicial (x0,y0). 1  f(x 0 , y 0 )  f(x1 , y1 )   derivada promedio 2

se usa la derivada promedio para calcular un nuevo valor de yi con la ecuación de Euler

yn+1 = yi + hf(xi, yi) que deberá ser más exacto que

y1.  y1  y 0 

h  f(x 0 , y 0 )  (f(x1 , y1 ))   2

que se tomará como valor definitivo de y1.

xn  x0   h    n  Este procedimiento se repite hasta llegar a yn. El esquema iterativo para este método quedaría en general así: Usando el paso de predicción resulta:

yi 1  yi  hf(xi ,yi ) Se calcula la derivada f ( x i  1 , y i  1 ) Se establece la derivada promedio (llamémosla B) B

1 f(x i , y i )  f(x i  1 , y i  1 ) 2





Se sustituye f(xi,yi) con este valor promedio en la ecuación de iteración de Euler y se obtiene y i 1  y i 

h  f(x i , y i )  f(x i 1 , y i 1 )   2

o simplificado: yi+1= yi+ hB Ejemplo 1: Resuelva f(x,y) = x - y en el intervalo [0,1] Con y(0)=2

x0=0

y0=2

xn=1

yn=?

y(1)=?

h  n=5

y Con 5 intervalos.

x

n

 x n

0

95



1 0 5

h = 0 .2

0

0.2

0.4

0.6

0.8

x0

x1

x2

x3

x4

1 x5

En forma iterativa, hagamos la siguiente tabla: de i=0 hasta n=5 Recuerde que:

xi+1=xi+h

f(xi,yi) = xi - yi

yi1  yi  hf(xi ,yi ) i

xi

0 0

Yi

2

f(xi,yi)

y i1

B



1 f(xi , y i )  f(xi 1 , y i 1 ) 2

Xi+1

f(x i 1 , y i 1 ) B

 Yi+1=yi+h b

-2

1.6

0.2

-1.4

-1.7

1.66

1 0.2 1.66

-1.46

1.368

0.4

-0.968

-1.214

1.4172

2 0.4 1.4172

-1.0172

1.21376

0.6

-0.61378

-0.81548

1.254104

3 0.6 1.254104

-0.654104

1.12328

0.8

-0.32328

-0.488694

1.156365

4 0.8 1.156365

-0.356365

1.085092 1.0

-0.08509

-0.2207286

1.112220

5 1.0 1.112222

Para i=5

y5=y(1)=1.112222 El valor exacto es 1.10364

1.10364  1.112222  100% 1.10364 E R  0.78%

 ER 

8.3.

MÉTODOS DE RUNGE-KUTTA (RK)

Los métodos de Runge-Kutta tienen la exactitud del esquema de la serie de Taylor sin necesidad de cálculo de derivadas superiores y consisten en obtener un resultado que se obtendría al utilizar un número finito de términos en una serie de Taylor de la forma: yi 1  yi  f ( xi , yi )h 

f ' ( xi , yi )h 2 f " ( xi , yi )h 3 f ' ' ' ( xi , yi )h 4   2! 3! 4!

+.......... con una aproximación en la cual se calcula yi+1 de la fórmula:

Yi 1  yi  h( 0 f ( xi , yi )  1 f ( xi  1h, yi  b1h)  ... ....   p f ( xi   p h, yi  bp h)

96

(2)

donde los

,  ,

b

se determinan de modo que si expandiera

f ( xi   j h, yi  bj h) con i  j  p en series de Taylor alrededor de (xi, yi), se observaría

que los coeficientes de h, h2 , h3 , etc. Coincidirían con los

coeficientes correspondientes de la ecuación (1). Después de efectuar los cálculos convenientes a la ecuación (2) para establecer los valores de  ,  , b , tenemos las siguientes expresiones de Runge-Kutta (RK):

1. Método RK de segundo orden: yi+1 = yi +

h [k1 +k2] 2

donde: k1=f(xi, yi) k2=f(xi+h, yi+hk1) 2. Método RK de tercer orden: yi+1 = yi + donde:

h [k1 +4k2 +k3] 6

k1=f(xi, yi) k2=f(xi+

h h , yi+ k1) 2 2

k3= f(xi+h, yi -hk1+2hk2) 3. MétodoRK de cuarto orden: yi+1 = yi +

h [k1 +2k2 +2k3 +k4] 6

donde: k1=f(xi, yi) k2=f(xi+

h h , yi+ k1) 2 2

k3=f(xi+

h h , yi+ k2) 2 2

k4=f(xi+h, yi+hk3) otra forma de expresar

97

k1  hf ( x n , y n ) k h k 2  hf ( x n  , y n  1 ) 2 2 k h k 3  hf ( x n  , y n  2 ) 2 2 k 4  hf ( x n  h, y n  k 3 ) 1 y n 1  y n  (k1  2k 2  2k 3  k 4 ) 6



El método más utilizado es el de cuarto orden ya que coincide con los primeros cinco términos de la serie de Taylor lo cual significa gran exactitud sin cálculo de derivadas aunque haya que evaluar la función f(x,y) cuatro veces en cada subintervalo.

Ejemplo: Aplíquese el método de Runge-Kutta de cuarto orden para f(x,y)=x-y donde x=0 hasta x=1 con 5 intervalos y siendo para x=0, y=2.

Ejemplo: 1.

Aplíquese el método de Runge-Kutta de cuarto orden para f(x,y)= x – y desde x=0 hasta x=1 Con 5 intervalos siendo para x=0 y=2. Solución: h

(x0, y0)= (0, 2)

1 0  0, 2 5

para x=1 y=? x0=0

x1=0.2

x2=0.4

Aplicamos: yi+1 = yi +

x3=0.6

x4=0.8

x5=1.0

h [k1 +2k2 +2k3 +k4] 6

Se calculan los valores de k1, k2, k3, k4 en cada iteración y se halla yi+1

i xi

Yi

K1

K2

K3

K4

Yi+1

0 0

2

-2

-1.7

-1.73

-1.454

1.6562

1 0.2 1.6562

-1.4562

-121058 -1.235142 -1.009172 1.410973

2 0.4 1.410973 -1.010973 -0.809876 -0.829985 -0.644976 1.246451 3 0.6 1.246451 -0.746451 -0.471806 -0.49927 -0.346797 1.148004 4 0.8 1.148004 -0.348004 -0.213204 -0.226684 -0.10268 1.103656 5 1.0 1.103656 98

para x=1.0, y=1.103656 Que al comparar con el valor exacto que es 1.10364, tenemos: ER 

2.

1.10364  1.103656  100%  0.001% 1.10364

Resolver F ( X i , Yi ) 

1 1  X i  Yi 2 aplicando el método de Runge-Kutta. 2

Solución De la condición inicial del problema se tiene que X = 0, y Y = 1; además, h = 0.1. Sustituyendo estos valores en se obtiene:

Llevando estos valores a (16) y el resultante a (12) se obtiene que para X = 0.1 la solución del problema es

Los valores de las kipara este punto obtenido de la solución, son:

luego

Continuando de la misma forma se obtiene la solución que se muestra en la siguiente tabla: 99

8.4.

X

Y

k1

k2

k3

k4

0.0

1.0000

0.5000

0.5516

0.5544

0.6127

0.1

1.0554

0.6126

0.6782

0.6823

0.7575

0.2

1.1236

0.7575

0.8431

0.8494

0.9494

0.3

1.2085

0.9492

1.0647

1.0745

1.2121

0.4

1.3158

1.2119

1.3735

1.3896

1.5872

0.5

1.4545

1.5868

1.8234

1.8517

2.1509

MÉTODO DE DIFERENCIAS FINITAS

Las diferencias

divididas finitas se sustituyen por las derivadas en la

ecuación original. Así una ecuación diferencial se transforma en un conjunto de ecuaciones algebraicas lineales simultaneas. Primera derivada:

dy yi1  yi  dx x

Segunda derivada:

d 2 y yi 1  2 yi  yi 1  dx 2 x 2

Resolver la ecuación mediante diferencias finitas

d2 y  2  106 y  7.5  107  (75  x) 2 dx i 1

i 1

i

d 2 y yi 1  2 yi  yi 1  dx 2 (x) 2

yi 1  2 yi  yi 1  2  106 yi  7.5  107 xi (75  xi ) 2 (x) i 1 x0

i2

i 3

i4

x  25

x  50

x  75

x  25 , x  0 a x  75 con x  25 .

100

x0  0 x1  x0  x  0  25  25 x 2  x1  x  25  25  50

x3  x2  x  50  25  75 x  0 y1  0

y3  2 y 2  y1  2  10 6 y 2  7.5  10 7 x2 (75  x2 ) 2 (25) 0.0016 y1  0.003202 y 2  0.0016 y 3  7.5  10 7 ( 25)(75  25) 0.0016 y1  0.003202 y 2  0.0016 y 3  9.375  10 4

y 4  2 y3  y 2  2  10 6 y3  7.5  10 7 x3 (75  x3 ) 2 (25) 0.0016 y 2  0.003202 y 3  0.0016 y 3  7.5  10 7 (50 )(75  50 ) 0.0016 y 2  0.003202 y 3  0.0016 y 3  9.375  10 4

x  75

y4  0  0 0 0   y1  0  1  y   0.0016  0.003202  0.0016 0   2  9.375  10 4     0 0.0016  0.003202 0.0016  y3  9.375  10 4       0 0 1   y 4  0  0   y1  0  y     2    0.5852  y 3   0.5852       y 4  0

8.5.

SISTEMA DE ECUACIONES DIFERENCIALES DE PRIMER ORDEN

Ya sea que estemos afrontando un problema de ingeniería cuya solución implique la resolución de una ecuación diferencial de orden n, o uno que implique la resolución de un sistema de ecuaciones diferenciales ordinarias de primer orden, nos enfrentaremos a la necesidad de resolver un sistema que podremos expresar como sigue:

101

dy1  f1  x , y1 , y 2 ,............, y n  dx dy 2  f 2  x , y1 , y 2 ,............, y n  dx . . dy n  f n  x , y1 , y 2 ,............, y n  dx

Por supuesto, la solución de tal sistema requiere de que se conozcan las n condiciones iniciales en el valor inicial del intervalo correspondiente a x. Todos los métodos vistos anteriormente para simples ecuaciones pueden extenderse a la resolución de sistemas como el anterior. El procedimiento para resolver un sistema de ecuaciones simplemente involucra aplicar las técnicas conocidas para cada ecuación en cada paso, antes de proceder con el siguiente. Esto quedará claramente ilustrado con el siguiente ejemplo donde hemos aplicado el método de Euler.

Ejemplo 1:Resolución de un sistema de EDOs mediante el método de Euler. Resolver el siguiente conjunto de ecuaciones diferenciales ordinarias: dy1  - 0, 5  y1 dx dy 2  4 - 0, 3  y 2 - 0,1  y1 dx

Resolveremos el sistema en el intervalo entre x = 0 y x = 2, con condiciones iniciales en x = 0, y1 = 4 , y2 = 6. Utilizaremos un paso h = 0.5. Se implementa el método de Euler para cada variable mediante la ya conocida expresión:

yi 1  yi  f ( xi , yi )  h Primero calculamos las pendientes: dy1  f1 ( 0 , 4, 6)  - 0.5  4  - 2 dx dy 2  f 2 ( 0 , 4, 6)  4 - 0.3  6 - 0.1  4  1.8 dx

Y luego los valores de la función para el primer paso:

102

y1  0.5  y1  0  f1 ( 0 , 4 , 6)  h  4   -2 0.5  3 y2  0.5  y2  0  f 2 ( 0 , 4 , 6)  h  6  1.8 0.5  6,9 Para un segundo paso volvemos a calcular las pendientes: dy1  f1 ( 0.5 , 3, 6.9)  - 0.5  3  - 1.5 dx dy 2  f 2 ( 0.5 , 3, 6.9)  4 - 0.3  6.9 - 0.1  3  1.63 dx

Y luego los valores de la función para el segundo paso:

y1 1.0  y1  0.5  f1 ( 0.5 , 3, 6.9)  h  3   -1.5  0.5  2.25 y2 1.0  y2  0.5  f2 ( 0.5 , 3, 6.9)  h  6.9  1.63  0.5  7.715 Y así continúa el cálculo hasta el final. Los resultados se resumen en la siguiente tabla: x

y1

y2

0.0

4.000000

6.000000

0.5

3.000000

6.900000

1.0

2.250000

7.715000

1.5

1.687500

8.445250

2.0

1.265625

9.094870

8.6 PROBLEMAS

1. Dada la ecuación diferencial f(x,y)=yx2-y con x=0 a x=2

siendo y(0)=1 y

n=4 intervalos resuelva por el método sencillo de Euler 2. Dada y’=f(x,y) = 2x3 –3x2 entre x=0 y x=1 siendo f(0)=1 y con n=2, evalúe utilizando a)el método sencillo de Euler. 3. Dada y’=f(x,y)= -2x3 + 12x2 –20x + 8.5 resuelva utilizando el método RK de cuarto orden desde x=1 hasta x=3 considerando 4 intervalos y siendo y(1)=3. 4. Utilice el método de Euler Modificado para resolver: a) dy/dx=2x3 –2xy con y(0)=0

y(2.5)=? 103

h=0.5

b) y’=2x2 –3y2 con y(1)=0.5

y(2)=?

h=0.2

5. Resolver mediante el método euler dy  100000( y  e  x )  e  x dx

y (0)  0 , use the implicit Euler method to obtain the solution from x  0 to 0,3

using the step-size of 0.1 6. Un tanque de 600 galones de capacidad contiene inicialmente 200 galones de salmuera con 25 lb. de sal. Salmuera con 2 lb. por galón entra al tanque a un flujo de 13 galones/ s. La salmuera mezclada en el tanque fluye hacia fuera a un flujo volumétrico de 8 galones/ s. ¿Cuál es la cantidad de sal y la concentración cuando el tanque se Encuentra lleno? Utilizar el método de runge kutta de segundo orden. 7. Las plantas de desalinización se usan para purificar el agua de mar , para que pueda beberse . El agua de mar contiene disuelto 8 g de sal /kg. Y es bombeada hacia un tanque de mezcla a razón de 0,5 kg./ min. Debido a una falla en el diseño del equipo, el agua se evapora del tanque a razón de 0,5kg/min. .La solución salina sale del tanque a razón de 10 kg/min. Supóngase que el balance de la disolución es agua pura. Si el tanque se llena inicialmente con 1 000 kg. de solución a la entrada. ¿ determine la concentración a la salida del tanque en función del tiempo hasta que el tanque quede vació.? 8. Las siguientes ecuaciones definen las concentraciones de tres reactores

dCA  20CACC  2CB dt dCB  20CACC  2CB dt dCC  20CACC  2CB  0,2CC dt Si las condiciones iniciales son CA = 500, CB = 0 y Cc = 500 halle las Concentraciones para los tiempos que van desde 0 a 30s usando un h = 3

9. La descarga de un fluido por gravedad por un orificio por el fon deo del tanque puede ser modelado mediante el siguiente sistema d e ecuaciones diferenciales. 104

dv  0,0010 y  0,00205 v 2 dt dy 1  0 , 0624 (1  ) dt 100  v

donde v es la velocidad del fluido, y nivel del liquido las condiciones iniciales son de v0 = 2,5 ; y0 = 20 cual son los valores de v, y para t = 5 con h igual 0,5 10.

Par un circuito eléctrico los resistores no obedecen la ley de Ohm, y el

circuito dinámico se describe por la relación siguiente l

di i 3 i  R ( )   0 dt I  I

i = corriente eléctrica ; I = corriente de referencia igual 1; R resistencia eléctrica 2 ; Resolver la ecuación diferencial si i(0) = 0,6 A para

t=

30s con h = 3s 11.

Determine usado técnicas numéricas

la ecuación diferencial ordinaria

de segundo orden

d2x dx  3  10 x  0 2 dt dt

La condición inicial para la ecuación para t = 0

 x =1

Determine los valores de x entre 1,2 usando h = 0,1 12.

Un tanque

de 2 000 L contiene, en el inicio 400L de agua pura .

Comenzando en t = 0 una solución acusa que contiene 1 g/ L de cloruro de potasio fluye hacia el tanque a razón de 8L/s y al mismo tiempo, comienza a fluir una corriente de salida a razón de 4L/s .El contenido del tanque esta mezclado perfectamente y la densidad de la corriente de alimentación y salida puede considerarse constante.

Calcule la concentración de cloruro

de potasio en el tanque en el momento en que este rebosa.

13.

Un Reactor CSTR como se muestra en la figura, la concentración de

alimentación (CA0) of 0,5 mole/m3.

Determine la concentración al cabo de

0,5h

105

rA  V

 k1C A (1  k2 C A )

dC A  F (C A0  C A )  (V )(rA ) dt

V  2, 0m3

F  1, 0m3 / h

k1  1, 0h 1

k2  1, 0m3 / mole

14. Se tiene un intercambiador de calor de tubos concéntricos en contracorriente y sin cambio de fase .Las ecuaciones que describen el intercambiador de calor en ciertas condiciones de operación son: dTB  0.03 x TS  TB  dx dTS  0.04 x TS  TB  dx

CONDICIONES

TB (0)  20 0 C TS (0)  1000 C

TB (3)  ¿? TS (3)  ¿? Evalué la temperatura para x= 3 m

15.

Resolver la siguiente ecuación mediante diferencias finitas

d 2u 1 d u u   2  0, u(5)  0,08731, u(8)  0.0030769 dr 2 r dr r

106

16. El siguiente diagrama ( Figura.18) muestra tres tanques continuos con agitación conectados en serie

Se ha disuelto1500g de Na 2SO4 en el

primer tanque, llenando los

otros tanques con solvente puro, e

iniciando después un flujo de 40L/ a través del sistema. Calcular el tiempo para alcanzar una concentración de 0,01g/l en el tanque tres.

Fig Nº16: Sistema de tanques de mezclado.

17.

Determine el valor de y(1) resolviendo y t   0,05 y t   0,15 y t   0, y 0   0, y0  1

utilizando el método de Runge-Kutta de segundo orden, con h = 0,5

18. El circuito que se muestra en la figura E9.5 tiene una autoinductancia de L = 50H, una resistencia de R = 20 ohms y una fuente de voltaje de V = 10 volts. Si el interruptor se cierra en el instante t = 0, la corriente I(t) satisface la ecuación: L

d I t   RI t   E , dt

I 0   0

Determine el valor de la corriente para 0 < t  10 segundos, mediante el método de Runge-Kutta de segundo orden, con h = 0.1. L SW

i

E

R

Figura 19: Circuito electrico

107

19. El agua que sale del tanque entra a otro de 20 galones, en cual también se vierte agua pura a razón de 3 galones/minuto y se mezcla bien. La concentración de sal en el segundo tanque satisface y 2 t   

3 2 y 2 t   y1 , 20 20

y 2 0   10

donde y1(t) es la concentración de sal del tanque de 50 galones del problema anterior. Utilice el método de Runge–Kutta de segundo orden para determinar cuando alcanza su máximo la concentración de sal en el tanque de 20 galones. Suponga que el segundo tanque tiene agua pura en el instante t = 0.

20. Un tanque de 50 galones de agua contiene sal con una concentración de 10 onzas/galón. Con el fin de diluir el contenido de sal, se suministra agua pura a razón de 2 galones/minuto. Si el depósito tiene una mezcla uniforme y la misma cantidad de agua que entra sale del depósito cada minuto, la concentración de sal satisface y1 t   

2 y1 , 50

y1 0   10

donde y1(t) es la concentración de sal en onzas/galón y t es el tiempo en minutos. Utilice el método de Runge–Kutta de segundo orden con h = 1 minuto para determinar cuánto tiempo debe transcurrir para que la concentración de la sal sea 1/10 de su valor inicial. b) El agua que sale del tanque entra a otro de 20 galones, en cual tambien se vierte agua pura a razon de 3 galones/minuto y se mezcla bien. La concentracion de sal en el segundo tanque satisface y 2 t   

3 2 y 2 t   y1 , 20 20

y 2 0   10

donde y1(t) es la concentración de sal del tanque de 50 galones del problema anterior. Utilice el método de Runge–Kutta de segundo orden para determinar cuando alcanza su máximo la concentración de sal en el tanque de 20 galones. Suponga que el segundo tanque tiene agua pura en el instante t = 0.

108

21. Resolver el sistema de EDOs mediante el método de Runge – Kutta

dy1  - 0.5  y1 dx clásico de cuarto orden.

dy 2  4 - 0.3  y 2 - 0.1  y1 dx Resolveremos el sistema en el intervalo entre x = 0 y x = 2, con condiciones iniciales en x = 0, y1 = 4 , y2 = 6. Utilizaremos un paso h = 0.5.

109

9. ECUACIONES DIFERENCIALES PARCIALES (EDP) INTRODUCCIÓN Una EDP es una ecuación que tiene como incógnita a una función de dos o más variables y que involucra a una o más de sus derivadas parciales. El orden de una EDP es el de la derivada con mayor orden en la ecuación. La linealidad de las ecuaciones se establece como sigue: Si los coeficientes dependen sólo de las variables independientes entonces a la ecuación se le denomina lineal. Si además dependen de la propia función o de alguna de sus derivadas parciales entonces la ecuación es no lineal. EJEMPLO 

Ecuación diferencial parcial de transferencia de calor unidimensional en estado no estacionario.

T T 2  2 t x T: Temperatura en K , t es el tiempo en s,  es la difusivilidad térmica en m2/s



Ecuación de transferencia de calor en estado estacionario bidimensional T 2 T 2  2 x 2 y T es la temperatura en K , x,y ejes coordenados

9.1 PROBLEMAS 1. resuelva la ecuación diferencial parcial de transferencia de calor en dos dimensiones.

 2T  2T  0 X 2 Y 2

0  X 2 ,

0  Y  2 , X  Y  0,5

T( 0, Y) = 100 , T( X ,0 ) = 60 , T( 2,Y) = 100 , T(X, 2) = 60

2. Resuelva la siguiente ecuación diferencial de transferencia de calor dada por 110

 2T  2T  0 X 2 Y 2 T = 40 .

0  y 1 ,

x=0

T = 120

0 y  1

x =4

T = 40

0 X  4

y=0

T= 50

0 X  4

y=1

Utilizando y  0,25 y x  1

111

10. APLICACIONES DE INGENIERÍA QUÍMICA EN POLYMATHMATHCAD Y MATLAB 10.1 PROBLEMAS CON POLYMATH

1. Determine el volumen Molar del gas dióxido de carbono mediante la ecuación de VAN DER WAALS . a   P  V  b   RT V *V  

a= 3,592

b= 0,04267

R = 0,082

PARA T= 300K , P =100 ATM f(v) = (P + a/(V2))*(V-b) –R*T = 0 POLYMATH Results NLE Solution Variable V

Value 0.0793969

a

3.592

b

0.04267

T

300

R

0.082

P

100

f(x) -1.776E-13

Ini Guess 0.2505

2. En la batería de reactores de mezcla completa que se muestra en el diagrama, se lleva a cabo una reacción isotermica de Segundo orden a volumen constante si el flujo volumétrico a cada reactor es constante de 25 l/min. calcular la concentración de salida de cada reactor

112

Figura Nº17: Bateria reactores

Balance de materia en cada reactor: Reactor N°1: q*CA0 –q*CA1 –K*CA12*V1 =0 Reactor N°2 q*CA1-q*CA2 –K*CA22*V2 =0 Reactor N°3 q*CA2-q*CA3 –K*CA32*V3 =0 F(CA1) = q*(CA0 –CA1) –K*CA12*V1 =0 F(CA2) = q*(CA1-q*CA2) –K*CA22*V2 = F(CA3) = q*(CA2-q*CA3) –K*CA32*V3 =0 Variable

Value

f(x)

CA1

0.6031317

-7.105E-15

2

CA2

0.2869484

-1.332E-14

1

CA3

0.1725794

-2.944E-13

0.8

V

1200

K

0.08

CA0

2

Q

3. En

s

25

un reactor químico se obtiene

glucosa a partir de la hidrólisis de

almidón. La ecuación diferencial que gobierna el proceso para la concentración del almidón es .

V*Dca/dt = q*CA0 - q*CA + ra ra = -k*CA2V CA0: Concentración de la alimentación a la entrada = 6,3 mol/m3 113

CA: Concentración del almidon a la salida V : volumen del reactor = 500m3 Q : caudal volumetrico 100 m3/h ra : velocidad de transformación de almidon en

glucosa

3

K= 0.02 m /h*mol . Calcular la composición de salida para t = 20h ODE Report (RKF45) Differential equations as entered by the user [1] d(CA)/d(t) = q*CA0 /V - q *CA/V -K*CA2 Explicit equations as entered by the user [1] q = 100 [2] V = 500 [3] CA0 = 6.3 [4] K = 0.02 Independent variable variable name : t initial value : 0 final value : 20 Variableinitial valueminimal valuemaximal valuefinal value t

0

CA

6.3

q

20

20

4.38179

6.3

4.38179

100

100

100

100

V

500

500

500

500

CA0

6.3

6.3

6.3

6.3

0.02

0.02

0.02

0.02

K

0

4. La concentración saturada de oxigeno disuelto en agua como función de temperatura y de la concentración de cloro se lista en la siguiente tabla. *Determine la ecuación que correlacione C=f(T) de grado 2 y 3 

Evaluar la concentración para T = 20°C

114

Temperatura

Cloro( 10 000mg/l)

5

11,6

10

10.3

15

9,1

20

8,2

25

7,4

30

6,8

REGRESION LINEAL Model: C = a0 + a1*T + a2*T2 Variable Value 95% confidence a0

13.11

0.1595539

a1

-0.3195

0.0208769

a2

0.0036429

5.839E-04

REGRESION POLINOMICAt Model: C = a0 + a1*T + a2*T2 + a3*T3 Variable Value 95% confidence a0

13.133333

0.5263102

a1

-0.3253704

0.1199362

a2

0.0040317

0.0076758

a3

-7.407E-06

1.451E-04

5. La siguiente secuencia de reacciones se llevan a cabo en estado no estacionario según la reacción:

Determinar la ecuación de cada componente de la reacción así como los productos en función al tiempo. Construya la tabla tiempo vs concentración en un paso de 1 a 1 hasta 10 minutos, asuma que las concentraciones iniciales son A (0)=1 mol/L, B (0)=1 mol/L, C (0)=1 mol/L. Datos: 115

K1=0.5 min-1 K2=1.0 min-1 Solución: formando las ecuaciones diferenciales en función de la reacción. C1   k1C1 t C2  k1C1  k2C2 t C2  k 2 C2 t

FIGURA 18: Datos de ingreso en Polymath

FIGURA 19: Gráfica concentración vs tiempo 116

6. Resolver el siguiente sistema: y   k1 x  k2 y t x  k1 x t z  k2 y t

Con valores iníciales x (0)=1, y (0)=0, z(0)=0 con t=0 a t=3 y k 1=1, k2=2 Solución:

FIGURA 20: Ingreso de Datos en Polymath

117

FIGURA 21: Resultados con el Software Polymath

FIGURA 22: Grafica de resultados de la Ecuación diferencial. 118

7. Resolver las siguientes ecuaciones:

dC1  0,12C1  0, 02C3  1 dt dC2  0,15C1  0,15C2 dt dC3  0, 025C2  0, 225C3  4 dt dC4  0,1C3  0,1375C4  0, 025C5 dt dC5  0, 01C1  0, 01C2  0, 04C5 dt

Con C1 (0)=1, C2 (0)=1, C3 (0)=1, C4 (0)=1 y C5 (0)=1 de t=0 a t=1 Solución:

FIGURA 23: Ingreso de Datos en Polymath

119

FIGURA 24: Resultados Polymath

FIGURA 25: Grafico Polymath

120

10.2 PROBLEMAS CON MATHCAD

1.

Para flujo turbulento para todas las tuberías, el Instituto de Hidráulica de Estados Unidos y la mayoría de ingenieros consideran la ecuación de Colebrook como la más aceptable para calcular f. La ecuación es:

 e 0, 251  0,86 ln    0,36d Re f f 

1

  

Donde: Re

: l número de Reynolds (adimensional)

e

: aspereza o rugosidad de la tubería (unidad de longitud)

d

:, diámetro de la tubería (unidad de longitud)

Obtener el factor de fricción para un fluido con un Reynolds de 3E4 que fluye en una tubería con un diámetro de 0,1 m y una rugosidad de 2,x10-3 m. Solución con mathcad Ingresando datos: E = 0,002 D = 0,1 Re = 3104 f = 0,001 Given

 E 0, 251  0,86 ln    f  0,36 D Re f

1

  

Calculo de la fricción F = Find (f) f = 0,055 de la figura efectuar un balance de masa y determine los flujos

121

Figura Nº26: Sistema de columna de separación

0.07x1 +0.18x2 + 0.15x3 +0.24x4

= 10.5

0.04x1 + 0.24x2 + 0.1x3 +0.65x4 =17.5 0.54x1 + 0.42x2

+ 0.54x3 + 0.10x4 = 28

0.35x1 + 0.16x2 + 0.21x3 + 0.01x4 = 14  0.07 0.04 A    0.54  0.35 

0.18 0.15 0.24  0.24



0.1 0.65 

  0.160 0.21 0.01  0.42 0.54 0.1

122

 10.5     26.25  17.5     B   17.5  solucion   28   8.75   14   17.5      2.

La ecuación de Antonie para el cálculo de la presión de vapor (P*) de un compuesto esta dado por: Donde: P*: presión de vapor (mmHg) T, temperatura (ºC)

Las constantes para el benceno y el hexano son:

Compuesto

A

B

C

benceno

6.89745

1260.350

220.237

Hexano

6.87773

1171.530

224.366

Determinar la temperatura de ebullición de la mezcla a 1 y a 5 atm. Solución

log ( P*)

A

B C T

sean las fracciones molares xa = 0,5 xb = 0,5 constantes para el benceno A1  6.89745

B1  1206.350 C1  220.237

constantes para el hexano A2  6.87773

B2  1171.530 C2  224.366

para una presión de 1atm en mmHg 123

Pt  760

asumiendo una temperatura inicial de 30ºC T  30 Given

 A1 B1   A2 B2      C1  T  C2  T    Pt Xa 10  Xb 10 T  Find(T) T  73.975

3.

Resolver el siguiente sistema de ecuaciones no lineales con mathcad

2

k1

k2

x y

2

( 1  x)  ( 1  x  y ) y ( x  y) ( x  y)  ( 1  x  y)

Solución

Especifique las constantes de equilibrio Tome dos valores iniciales Given 2

k1

k2

x  y

2

( 1  x)  ( 1  x  y ) y (x  y ) (x  y )(1  x  y )

vec  Find ( x  y ) vec 

 0.652     0.198  124

k1  0.6

k2  0.2

x  0.8

y  0.5

4..

Mediante la ecuacion de Vand Der Waals calcular el volumen molar a P = 60 atm y T = 500K.

 P  a   ( V  b ) R T  2 V  

ecuacion de Van der Waals

Definiendo las constantes R  0.082 Para el amoniaco

Tc  405

Pc  111.3

Constante a y b 27  R  Tc 2

a   64 

2

 

Pc

b 

R Tc 8 Pc

Especificando las constantes R  0.082

Para el amoniaco

Constante a y b a 

Tc  405

 R 2  Tc  64  Pc 27

Evaluando las constantes a y b

2

a 

Especifique la presion y temperatura

b 

R  Tc 8  Pc b 

T  500

Asumiendo un valor inicial( comportamiento ideal) R T V  V  P Dado

 P  a  (V  b )  2  V   calculo del volumen Molar V  Find

  

Pc  111.3

(V)

V 

125

R T

P  60

5.

Encontrar el volumen molar que ocupa el gas empleando la ecuación de Redlinh Kwong RT a P  V  b V (V  b) T  R 2TC5 / 2   RTc  a  0, 42747   b  0, 08664    Pc   PC 

 Atm  L    mol  gK 

P en Atm V en L/mplg T en K R  0, 082 

Pc: presión critica de amoniaco 111,3 atm. TC: temperatura critica de amoniaco 405K ECUACION DE ESTADO DE R-K Ingresando las contantes: Ingresando los datos:

R  0.082 P  60

Pc  111.3

Tc  405

T  500

Calculo de a yb: 5   R 2  Tc 2 a  0.42747   Pc

  

b  0.08664 



a 

b 

Asumiendo comportamiento ideal

V  V

Dado P

R T (V  b)



V  Find ( V)

a V ( V  b )  T V

126

R T P

R  Tc Pc

  

6.

A un sistema de separación FLASH ( Fig. 21) ingresa 100 moles/h de una mezcla

VAPOR(V)

2

yi

1 Alimentacion F zi

3

xi

Liquido(L) Figura Nº27:Sistema de separación de componentes Alimentación al Flash( Mol/h) Zi: Fraccion molar de los componentes de la mezcla i = 1,2,3, …NC V: Caudal de salida del vapor ( mol/h) yi : Fraccion molar de los componentes en la fase vapor L: Caudal De salida de la fase liquida

( mol/h)

xi: Fraccion molar de los componentes en la fase liquida Balance de materia global: F = V + L Balance de componentes : Zi*F = xi*L + yi*V Equilibrio Liquido vapor : yi = Ki *xi V F Remplazando en las anteriores ecuaciones expresiones

Fraccion vaporizada :

x



zi  ( Ki  1)  1 127

Se obtiene las siguientes

yi 

Ki * zi  ( Ki  1)  1

Ni

Zi( Ki  1)

  (Ki  1)  1  0 i 1

Teniendo en cuenta los datos de la siguiente tabla Calcular la fracción vaporizada así como los flujos en cada corriente con sus respectivas composiciones COMPONENTES

FRACCION MOLAR

C1

0.05

16.25

C2

0.15

5.25

C3

0.25

1.99

C4

0.2

0.75

C5

0.35

0.29

Ingresando los datos:

F  100 k3  1.99

k1  16.25 k4  0.75

z1  0.05 z3  0.25

K

k2  5.25 k5  0.29

z2  0.15 z4  0.2

Asumiendo un valor:

z5  0.35

  0.5

Dado z1 ( k1  1)  ( k1  1)  1



  Find( )

z2 ( k2  1)  ( k2  1)  1



z3 ( k3  1)  ( k3  1)  1



128



z4 ( k4  1)  ( k4  1)  1



z5 ( k5  1)  ( k5  1)  1

0

7.

Encontrar el volumen molar que ocupa el gas empleando la ecuación de estado de Redlinh-Kwong.

A una presión de 14 bar y una temperatura de 333K

RT a  V  b V (V  b) T

P

Ingresando los datos: T  333 a b

f ( V)

P T

b  44.897

8

a  1.5641410 

R  83.14

P  14

  b  2

b  R T



P



  V  R T  V2  V3  P P T  a

Se observa un plinomio de grado 3 Ingrese los coeficientes en una matriz de 4*1

     2  b  V        

  P T  b  R T a    P P T   R T    P  1  a b

S  polyroots ( V)

8.

 71.299    229.996 S   3  1.676  10 

La temperatura media logarítmica de un intercambiador de calor a contracorriente esta dado por:

LMDT 

T1  t2   T2  t1  T t  ln  1 2   T2  t1 

129

Para este sistema la temperatura media logarítmica debe ser de 50°C y el fluido caliente se alimenta a 100°C y sale a 60 °C, mientras que el fluido frió se alimenta a 15°C ¿A qué temperatura sale del intercambiador el fluido frió?

Ingrese los datos : T1  15

T3  100

T4  60

Asumiendo un valor inicial:

LMTD  50

T2  50

Dado LMTD

( T3  T2)  ( T4  T1) ln

T3  T2 

  T4  T1 

T2  Find( T2) T2 

9.

El factor de fricción para diferentes números de Reynolds empleando la ecuación de Colebrook

  1 2,51   0,86 ln    f  3, 7d Re f  Obtener el factor de friccion para Re de 10 000 , e = 0,002cm y diámetro de 2cm CALCULO DEL FACTOR DE FRICCION

Ingresando datos:

Re  10000

Asumiendo para f : f  0.001 Dado 0.86 ln

1

e

 3.7 d

f



  Re f  2.51

f  Find( f ) f 130

e  0.002

d  2

10.3. PROBLEMAS CON MATLAB 1. A un reactor ingresa una mezcla de gases de 8% de dióxido de azufre y el 12% de oxigeno y la diferencia de nitrógeno, y se desarrolla la siguiente reacción. 1 SO2  O2  SO3 2

Calcular la composición en el equilibrio a presión constante de 2 atm y la constante de equilibrio KP es de 160. Algorritmo de Solucion en Matlab: clear all, close all, clc disp('H2SO4'); SO2=input('Ingrese moles de SO2: '); O2=input('Ingrese moles de O2: '); SO3=input('Ingrese moles de SO3: '); N2=input('Ingrese moles de N2: '); P=input('Ingrese Presion Total (atm): '); KP=input('Ingrese Constante de Equilibrio: '); syms x; disp('Moles en Equilibrio: '); ESO2=SO2-x; EO2=O2-0.5*x; ESO3=SO3+x; disp(ESO2); disp(EO2); disp(ESO3); ET=ESO2+EO2+ESO3+N2; f=KP-((ESO3/ET)*P)/((ESO2/ET)*((EO2/ET)^0.5)*P^1.5); 131

xai=input('Ingrese Limite Inferior: '); xbi=input('Ingrese Limite Superior: '); tol=input('Ingrese Tolerancia: '); f=inline(f); i=1; ea(1)=100; if f(xai)*f(xbi) < 0 xa(1)=xai; xb(1)=xbi; f1=f(xai); xr(1)=(xa(1)+xb(1))/2; fprintf('Iter.

x1

f1

x2

Error aprox \n');

fprintf('%2d \t %11.7f \t %11.7f \t %11.7f \n',i,xa(i),f1(i),xr(i)); while abs(ea(i)) >= tol, if f(xa(i))*f(xr(i))< 0 xa(i+1)=xa(i); xb(i+1)=xr(i); end if f(xa(i))*f(xr(i))> 0 xa(i+1)=xr(i); xb(i+1)=xb(i); end xr(i+1)=(xa(i+1)+xb(i+1))/2; ea(i+1)=abs((xr(i+1)-xr(i))/(xr(i+1))*100); 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; 132

end else fprintf('No existe una raíz en este intervalo'); end Resultados: TABLA 12 Resultados en Matlab

2. Se sabe que la temperatura media logarítmica de un intercambiador de calor contracorriente está dada por:

 ΔT  AQ

T = LMTD =

ent cal

sald - T frio  - Tcalsald -T frioentr 

sald  Tcalent - T frio ln  sald entr  T -T frio  cal

133

  

Se tiene el siguiente intercambiador: entr = 15ºC T frio

INTERCAMBIADOR DE CALOR

sald = T frio

Tcalsald= 60ºC

t

Tcalent = 100ºC

Se sabe que en este sistema la temperatura media logarítmica debe ser de 50ºC. El fluido caliente se alimenta al sistema a 100ºC, mientras que el fluido frío se alimenta a 15ºC,. ¿a que temperatura sale del intercambiador de fluido frío? Algoritmo de Solución en Matlab: clear all, close all, clc disp ('Temperatura en un Intercambiador de Calor') xo=input('Temperatura Inicial ='); n=input ('Numero de Iteraciones='); salida=ones(n,3); % matiz de salida de datos for i=1:n x1=xo-[(45*exp((55-xo)/50)+xo-100)]/[1+(-45/50)*(exp((55-xo)/50))]; vsal=[xo;x1]; ea=[[abs((x1-xo)/x1)]]; % error xo=x1; salida(i,1)=i; salida(i,2)=x1; salida(i,3)=ea; end disp('Iter.

T(i)

Error');

disp(num2str(salida)); Resultados:

134

TABLA 13 Resultados Temperatura en Matlab

3. Un gas se encuentra a una presión absoluta de 13.76 bar y una temperatura de 333 K. Encontrar el volumen molar ocupa el gas empleando la ecuación de estado de Redlich-Kwong.

P=

RT a - 1/2 V - b T V V + b 

Para este compuesto las constantes son: P = 13.76 atm T = 333ºk a = 1.5614 x 108 (cm6 bar/(g mol)2 k1/2) b = 44.897 (cm3/ g mol) R = 83.4 (cm3 bar /g mol k) Algoritmo de Solución en Matlab: clear all, close all, clc disp('Ecuacion de Redlich-kwong 2'); a=input('Ingrese a: '); b=input('Ingrese b: '); T=input('Ingrese T(K): '); P=input('Ingrese P(bar): '); R=input('Ingrese R(cm^3.bar/K.mol): '); syms x 135

f=R*T/(x-b)-a/((sqrt(T))*x*(x+b))-P; x0=(R*T)/(P); disp('Primer punto inicial: '); disp(x0); xai=x0; xbi=input('Ingrese Limite Superior: '); tol=input('Ingrese Tolerancia: '); f=inline(f); i=1; ea(1)=100; if f(xai)*f(xbi) < 0 xa(1)=xai; xb(1)=xbi; f1=f(xai); xr(1)=(xa(1)+xb(1))/2; fprintf('Iter.

V1

F1

V2

Error aprox \n');

fprintf('%2d \t %11.7f \t %11.7f \t %11.7f \n',i,xa(i),f1(i),xr(i)); while abs(ea(i)) >= tol, if f(xa(i))*f(xr(i))< 0 xa(i+1)=xa(i); xb(i+1)=xr(i); end if f(xa(i))*f(xr(i))> 0 xa(i+1)=xr(i); xb(i+1)=xb(i); end xr(i+1)=(xa(i+1)+xb(i+1))/2; 136

ea(i+1)=abs((xr(i+1)-xr(i))/(xr(i+1))*100); 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 else fprintf('No existe una raíz en este intervalo'); end Resultados: TABLA 14 Resultados volumen en Matlab

137

4. El factor de fricción (f) para el flujo turbulento en una tubería está dado por la correlación de Colebrook:

ε/D

1 = - 0.86 ln f

2.54

 +   3.4 Re

  f 

Donde Re = es el número de Reynolds (adimensional) , es la aspereza o rugosidad de la tubería (unidad de longitud) D, es el diámetro de la tubería (unidad de longitud) Obtener el factor de fricción para un fluido con un Reynolds de 3E4 que fluye en una tubería con un diámetro de 0.1 m y una rugosidad de 0.0025m. Algoritmo de Solución en Matlab: clear all, close all, clc disp('Friccion en un Flujo Turbulento'); Re=input('Ingrese el Numero de Reynolds: '); e=input('Ingrese la Rugosidad: '); D=input('Ingrese el Diametro de Tuberia: '); syms x; f=-sqrt(1/x)-0.86*log((e/D)/3.4+2.54/(Re*sqrt(x))); xai=input('Ingrese Limite Inferior del Intervalo: '); xbi=input('Ingrese Limite Superior del Intervalo: '); tol=input('Ingrese Tolerancia: '); f=inline(f); i=1; ea(1)=100; if f(xai)*f(xbi) < 0 xa(1)=xai; xb(1)=xbi; 138

f1=f(xai); xr(1)=(xa(1)+xb(1))/2; fprintf('Iter.

f1

F1

f2

Error aprox \n');

fprintf('%2d \t %11.7f \t %11.7f \t %11.7f \n',i,xa(i),f1(i),xr(i)); while abs(ea(i)) >= tol, if f(xa(i))*f(xr(i))< 0 xa(i+1)=xa(i); xb(i+1)=xr(i); end if f(xa(i))*f(xr(i))> 0 xa(i+1)=xr(i); xb(i+1)=xb(i); end xr(i+1)=(xa(i+1)+xb(i+1))/2; ea(i+1)=abs((xr(i+1)-xr(i))/(xr(i+1))*100); 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 else fprintf('No existe una raíz en este intervalo'); end Resultados: TABLA 15 Resultados volumen en Matlab

139

5. Determinar volumen molar del oxigeno mediante la ecuación del VAN DER WAALS

a    P  2 v  b   RT v   p = 100 Atm., T = 700 K para un gas que tiene a = 1,36 b = 0,0318 Algoritmo de Solución en Matlab: clear all, close all, clc disp ('Ecuacion de Van der Walls') T=input('Ingrese T(K): '); P=input('Ingrese P(atm): '); R=input('Ingrese R(L.atm/K.mol): '); 140

disp('Volumen Inicial: '); xo=R*T/P; disp(xo); n=input ('Numero de Iteraciones='); salida=ones(n,3); % matiz de salida de datos for i=1:n x1=xo-(100*xo^3-60.58*xo^2+1.36*xo-0.043248)/(300*xo^2121.16*xo+1.36); vsal=[xo;x1]; ea=[[abs((x1-xo)/x1)]]; % error xo=x1; salida(i,1)=i; salida(i,2)=x1; salida(i,3)=ea; end disp('Iter. V(i)

Error');

disp(num2str(salida));

Resultados:

TABLA 16 Resultados volumen en Matlab

141

6. La siguiente secuencia de reacciones se llevan a cabo en estado no estacionario según la reacción: k1 k2 A   B  C Determinar la ecuación de cada componente de la reacción así como los productos en función al tiempo. Construya la tabla tiempo vs concentración en un paso de 1 a 1 hasta 10 minutos, asuma que las concentraciones iníciales son A (0)=1 mol/L, B (0)=1 mol/L, C (0)=1 mol/L. Datos: K1=0,5 min-1 K2=1,0 min-1 Solución: formando las ecuaciones diferenciales en función de la reacción. C1   k1C1 t C2  k1C1  k2C2 t C2  k 2 C2 t

Metafile: function[f]=rxnserie(t,C); 142

%constantes K(1)=0.5; K(2)=1; f(1)=-K(1)*C(1); f(2)=K(1)*C(1)- K(2)*C(2); f(3)=K(2)*C(2); f=f'; file: tr=[0:1:10];C0=[1,0,1]; [t,C]=ode45('rxnserie',tr,C0); [t,C] TABLA 17 Resustado de concentracion en Matlab

143

plot(t,C); xlabel('tiempo (min)'); ylabel('concentración (mol/L)');

FIGURA 28 Resultados volumen en Matlab

144

UNIVERSIDAD NACIONAL DEL CALLAO FACUTLAD DE INGENIERIA QUIMICA ESCUELA PROFESIONAL DE ING. QUIMICA

SILABO

I.

II.

DATOS GENERALES 1.1 Asignatura

:

Métodos Numéricos

1.2 Código de la Asignatura

:

IG-301

1.3 Semestre Académico

:

2010-A

1.4 Ciclo

:

V

1.5 Créditos

:

04

1.6 Horas de teoría

:

06

1.7 Horas de práctica

:

06

1.8 Duración

:

16 semanas

1.9 Pre-requisito

:

IG 102, FM 202

1.10 Profesor de Teoría

:

Ing. Juan Medina Collana

SUMILLA 2.1 Naturaleza de la Asignatura El curso de Métodos Numéricos es un curso teórico, práctico de carácter obligatorio que sirve como base para las asignaturas de ingeniería, debido a que permite realizar cálculos complejos como resultado del modelamiento de fenómenos físicos y químicos.

2.2 Síntesis del Contenido Ecuaciones algebraicas no lineales. Sistema de Ecuaciones algebraicas no lineales. Interpolación Polinómica. Diferenciación. Integración. Ecuaciones Diferenciales Ordinarias. Sistema de ecuaciones Diferenciales Ordinarias. Sistema de Ecuaciones algebraicas lineales. Ecuaciones Diferenciales parciales. Análisis de regresión lineal, no lineal y multivariable.

145

III.

OBJETIVOS 3.1 Objetivo General Al término del curso, el alumno será capaz de resolver modelos matemáticos complejos que resulten del modelamiento de los fenómenos físicos y químicos, aplicando técnicas numéricas.

3.2

Objetivos Específicos  Realizar el estudio comparativo entre los métodos analíticos y métodos numéricos.  Adquirir destreza en el dominio de las técnicas de solución numérica.  Aplicar las técnicas numéricas para resolver problemas diversos de Ingeniería Química que resultan del modelameinto de fenómenos físicos y químicos.

IV.

PROGRAMA DE CONTENIDOS

PRIMERA SEMANA Ecuaciones algebraicas no lineales. Método de la bisección. Método de Newton Raspón de primer y segundo orden. Aplicaciones a la Ing. Química. Método de la secante. Método del punto fijo. Método de falsa posición. Métodos cuasi Newton. Métodos de la convergencia. Aplicaciones a la Ingeniería Química.

SEGUNDA SEMANA Sistema de ecuaciones no lineales. Método del punto fijo. Método de Newton Raspón estándar. Método de Newton Raspón modificado. Métodos cuasi Newton. Aplicaciones a la Ingeniería Química.

TERCERA SEMANA Interpolación polinómica. Ecuaciones en diferencias progresivas, regresiva y central de Newton. Aplicaciones a la Ingeniería Química.

146

Generación de polinomios mediante las diferencias progresiva, regresiva y central de Gauss. Interpolación con estaciones no equidistantes.

CUARTA SEMANA Análisis de Regresión Lineal y no lineal. Mínimos cuadrados. Regresión multivariable. Aplicaciones a la Ingeniería Química

QUINTA SEMANA Diferenciación numérica. Aproximación por diferencias. Diferenciación de los polinomios de Newton. Aplicaciones a la Ingeniería Química.

SEXTA SEMANA Integración numérica. Regla del Trapecio. Reglas de Simpson 1/3 y 3/8 y

SEPTIMA SEMANA Integración por cuadratura Integración múltiple. Aplicaciones a la Ingeniería Química.

OCTAVA SEMANA EXAMEN PARCIAL

NOVENA SEMANA Ecuaciones Diferenciales Ordinarias. Método de Euler. Método de EulerGauss. Método de Taylor.

DECIMA SEMANA Método de Runge-Kutta. Aplicaciones a la Ingeniería Química. Métodos Predictor-Corrector. Métodos de paso variable. Aplicaciones a la Ingeniería Química.

147

DECIMA PRIMERA SEMANA Sistema de Ecuaciones Diferenciales Ordinarias. Ecuaciones Diferenciales de Orden Superior: métodos de Ruge Kutta y matricial, explícita en la variable dependiente. Aplicaciones a la Ingeniería Química. DECIMA SEGUNDA SEMANA Sistema de ecuaciones algebraicas lineales. Métodos de Jacobi y GaussSeidal. Matrices y vectores. Sistemas especiales. Aplicaciones a la Ingeniería Química. DECIMA TERCERA SEMANA Ecuaciones Diferenciales Parciales. Ecuaciones en Diferenciales. Métodos explícito e implícito. Aplicaciones a la Ingeniería Química. Aplicaciones a la Ingeniería Química.

DECIMA CUARTA SEMANA Practica calificada

EXAMEN FINAL

EXAMEN SUSTITUTORIO

V.

PROCEDIMIENTO DIDÁCTICO Las clases teóricas se desarrollan en forma expositiva, y al mismo tiempo resolviendo ejercicios base. Con estos conocimientos se resolverán ejercicios de mayor complejidad con participación de los estudiantes. El profesor de práctica se encargará de la resolución de problemas aplicados a la Ingeniería Química con la finalidad de que los estudiantes se familiaricen con los tópicos que se estudian en Ingeniería Química.

VI.

EQUIPOS Y MATERIALES Para el desarrollo adecuado del curso, se requiere lo siguiente:  Calculadora científica.  Computadora con un software de programación de alto nivel.  Tiza blanca y de color. 148

VII. SISTEMA DE EVALUACIÓN El promedio final, se obtendrá del siguiente modo:

Peso  Examen Parcial

1.0

 Examen Final

1.0

 Promedio de Prácticas Calificadas

1.0

PF 

EP  EF  PP  11 3

Al final se tomará un Examen de Recuperación que sustituirá a la nota mas baja del Examen Parcial o Final. El examen de recuperación se tomará en base al total del avance del curso.

149

VIII. BIBLIOGRAFÍA

1. Carnahan, B. Luther, A. Wilkes

Cálculo

Numérico,

Aplicaciones

Editorial Rueda, Madrid, 1979.

2. Burden, R. Y Faires J.

Análisis Numérico. Edit. Iberoamericana, México, 1985.

3. Nieves, A., Domínguez, F.

Métodos Numéricos Aplicados a la Ingeniería Química. Edit. CECSA, México 1985.

4. Valderrama, R.

Métodos Numéricos. Edit. Trillas, UNAM, México, 1985.

5. Nakamura, S.

Métodos Numéricos aplicados con software. Edit. Prentice – Hall Hispano Americano, S.A. México, 1992.

6. Gerald, C.

Análisis Numérico Edit. Alfa – Omega, México, 1991.

7. Maron, M. y López R.

Análisis

Numérico.

Un

práctico. Edit. CECSA, México, 1998.

150

enfoque

Get in touch

Social

© Copyright 2013 - 2024 MYDOKUMENT.COM - All rights reserved.