MÉTODOS NUMÉRICOS

4.1 INTRODUCCIÓN

La simulación digital es un poderoso recurso para la solución de ecuaciones que describan los sistemas de ingeniería química. Las principales dificultades son dos:

1. Solución de ecuaciones simultáneas algebraicas no lineales (usualmente realizado por algún método iterativo)
2. Integración numérica de ecuaciones diferenciales ordinarias (usando ecuaciones discretas de diferencia finita para aproximar ecuaciones diferenciales continúas).

Se debe tener en cuenta la precisión y la estabilidad numérica de estas ecuaciones de aproximación. Tanto la precisión como la estabilidad son afectadas por la ecuación de diferencia finita (o algoritmo de integración) usada. Muchos algoritmos han sido propuestos en la literatura. Algunos trabajan mejor (por ejemplo, rápidos y entonces no tan exactos para un grado de exactitud especificado) en algunos problemas que otros. Desdichadamente no hay un algoritmo que trabaje mejor para todos los problemas. Sin embargo, como discutiremos con mayor detalle posteriormente, el algoritmo de Euler de   primer orden simple explicito, es el mejor para un gran número de aplicaciones de ingeniería.

A lo largo de los años se han desarrollado un gran número de paquetes de simulación. En teoría, estos lenguajes de simulación relevan al ingeniero del conocimiento acerca de la integración numérica, facilitando al ingeniero el establecimiento y la solución de los problemas. Estos automáticamente monitorean errores y estabilidad y ajustan el intervalo de integración o tamaño de paso para estar dentro de un criterio de precisión.

4.2 USO DE MATLAB

MATLAB proviene de la contracción de los términos MATRIX LABORATORY y fue inicialmente concebido para proporcionar fácil acceso a las librerías LINPACK y EISPACK, las cuales representan hoy en día dos de las librerías más importantes en computación y cálculo matricial.

MATLAB puede considerarse como un lenguaje de programación, como FORTRAN o C, aunque sería difícil describirlo en unas cuantas palabras. He aquí algunas de sus características notables para los análisis numéricos.

• La programación es mucho más sencilla
• Hay continuidad entre valores enteros, reales y complejos
• La amplitud de intervalo y la exactitud de los números son mayores
• Cuenta con una biblioteca matemática amplia
• Abundantes herramientas gráficas, incluidas funciones de interfaz gráfica con el usuario
• Capacidad de vincularse con los lenguajes de programación tradicionales
• Transportabilidad de los programas MATLAB

Una característica extraordinaria de los números en MATLAB es que no hay distinción entre reales, complejos, enteros, de preescisión sencilla y de doble precisión. En MATLAB todos estos números están conectados continuamente, como debe ser. Esto significa que en MATLAB cualquier variable puede contener números de cualquier tipo sin una declaración especial durante la programación, con lo cual está última se hace más rápida y productiva. En FORTRAN se requiere una subrutina distinta para cada variable sencilla o doble, real o compleja, o entera, mientras que en MATLAB no hay necesidad de separarlas. La biblioteca matemática de MATLAB facilita los análisis matemáticos. Además, el usuario puede crear rutinas matemáticas adicionales con mucha mayor facilidad que en otros lenguajes de programación, gracias a la continuidad de las variables reales y complejas. Entre las numerosas funciones matemáticas; Los solucionadores de álgebra lineal desempeñan un papel crucial; de hecho, todo el sistema MATLAB se basa en estos solucionadores

Importancia de las gráficas

El análisis visual de los modelos matemáticos ayuda a comprender las matemáticas y a hacerlas más asequibles. Aunque esta ventaja es bien conocida, la presentación de resultados calculados con gráficos de computadora solía requerir un esfuerzo adicional considerable. Con MATLAB, en cambio, bastan unos cuantos comandos para producir presentaciones gráficas de material matemático. Es posible crear objetos gráficos, científicos e incluso artísticos en la pantalla mediante expresiones matemáticas.

Finalmente los lenguajes de programación como FORTRAN o C, siguen siendo importantes para computación de alto rendimiento que requiere abundante memoria o un tiempo de computo largo. La velocidad del computo con MATLAB es significativamente más baja que con FORTRAN o C porque MATLAB paga un precio elevado por sus características agradables. Por otro lado, no es necesario aprender FORTRAN ni C para entender MATLAB.

Libros de referencia para importantes para aprender MATLAB

En este texto se dan algunos ejemplos sobre la aplicación de MATLAB en simulación de procesos químicos, pero de ninguna manera pretende ser una guía del lenguaje. Se recomienda a los lectores interesados en información adicional sobre MATLAB consultar la siguiente literatura:

The Math Works, The Student Edition of MATLAB,  Versión 4, User’s guide, Prentice –Hall

MATLAB, Reference Guide, Mathworks

MATLAB, User’s Guide, Mathworks

MATLAB, Building a Graphical User Interface, Math Works

Análisis Numérico y Visualización Gráfica con MATLAB 

4.3 ÁLGEBRA LINEAL

Son dos las razones principales por las que conviene aprender álgebra lineal en esta etapa tan temprana del estudio de los métodos numéricos. En primer lugar, el álgebra lineal es fundamental para los métodos numéricos, así que cuanto más pronto la aprendamos, más fácil nos resultará el resto del estudio de los métodos numéricos. En segundo lugar, las capacidades de MATLAB se basan en las operaciones de matrices y vectores; por tanto. Podremos utilizar mucho mejor MATLAB si aprendemos álgebra lineal.

El objetivo de esta sección es comprender los fundamentos del álgebra lineal y poder resolver ecuaciones lineales, sobre todo con MATLAB. Sin embargo, esto requiere cuatro áreas de estudio. Primero debemos aprender a expresar las ecuaciones lineales y sus operaciones básicas en notación de matrices y vectores. Segundo debemos entender las ordenes de MATLAB que trabajan con ecuaciones lineales en notación de matrices. Tercero, debemos entender los problemas que son difíciles o imposibles de resolver. La cuarta área incluye temas adicionales que amplían la comprensión y que ayudan a aplicar el álgebra lineal.

En notación matemática, las matrices se encierran entra paréntesis o corchetes y siguen ciertas reglas matemáticas. En MATLAB, en cambio las matrices se escriben sin paréntesis ni corchetes. Aunque en MATLAB los términos Matriz y arreglo son sinónimo, debemos estar concientes de si un arreglo se está utilizando para una matriz en matemáticas o simplemente como variable de arreglo. En esta sección trataremos sobre las aplicaciones de matrices y vectores en el modelamiento de procesos químicos y su manipulación usando MATLAB

4.3.1 Solución de ecuaciones simultáneas

 Los problemas científicos siempre tienden a involucrar la solución de ecuaciones lineales simultáneas. Para pocas ecuaciones (dos o tres a lo sumo) esto puede hacerse a mano. Para más de tres ecuaciones es aconsejable usar una computadora.  

Ejemplo 4.1  Modelamiento con sistema de ecuaciones algebraicas

La Fig. 4.1 muestra la disposición de un problema donde las computadoras no solo son útiles, sino también prácticas debido a que hay cinco reactores interconectados o acoplados, se necesitan cinco ecuaciones de balance de masa para caracterizar el sistema.

Fig. 4-1 Cinco reactores conectados por tuberías

Si F es flujo volumétrico (volumen/tiempo) y C es la concentración (unidades de masa/volumen). Para el reactor 1, la razón de flujo de masa que entra es:

5(10) +F₃₁C₃

y la razón de flujo de masa de salida es

F₁₂C₁ + F₁₅C₁

como el sistema se encuentra en estado estacionario, los flujos de entrada y salida deben ser iguales:

5(10) + F₃₁C₃  =  F₁₂C₁ + F₁₅C₁

o sustituyendo los valores para los flujos volumétricos de la Fig. 4.1 se tiene

6C₁ – C₃ = 50

Haciendo balance de masa para los otros reactores al estado estacionario tenemos:

6C₁ – C₃ = 50

–3C₁ + 3C₂ = 0

– C₂ + 9C₃ = 160                                                                                  (4.1)

– C₂ – 8C₃ + 11C₄ – 2C₅ = 0

– 3C₁ – C₂ + 4C₅ = 0

o lo que es igual:

6C₁           –  C₃                       = 50

–3C₁ + 3C₂                                = 0

–  C₂ + 9C₃                       = 160

–  C₂ – 8C₃ + 11C₄ – 2C₅ = 0

–3C₁ –  C₂                      + 4C₅ = 0

Lo cual en notación matricial será:

Existen diferentes métodos para resolver este modelo.

1. Aplicando MATLAB, hacemos lo siguiente:

1.1    Escribimos  la matriz A y el vector b

>> A=[6 0 -1 0 0; -3 3 0 0 0;0 -1 9 0 0; 0 -1 -8 11 -2; -3 -1 0 0 4];

>> b=[5*10;0;8*20;0;0];

1.2    Damos la orden para calcular x:

>> x=A\b

4.3.2. Análisis de sistemas de reacción

Ejemplo 4.2  Una situación común en el diseño de sistemas de reacción es cuando existe un número grande de reacciones en el sistema. En este caso se deben identificar cuantas y cuales son las reacciones independientes. Estas reacciones describen al sistema y para su diseño es suficiente dirigir el estudio hacia estas reacciones. Un caso típico es el diseño de un alto horno en el cual se llevan a cabo cientos de reacciones, pero solamente unas cuantas son las que influencian al sistema.

Esta situación también se presenta en un mecanismo de reacción para el análisis de la etapa controlante del proceso.

El sistema de reacción que corresponde a la formación del HBr puede representarse por un mecanismo apropiado. En este mecanismo aparecen las sustancias Br₂, H₂, Br°, H°, HBr. Utilizando el mecanismo y representando en una matriz de coeficientes estequiométricos dicho sistema:  determine cuantas reacciones son independientes?

Solución

a)             Formulación del mecanismo

b)     Haciendo: A₁ = Br₂;  A₂ = Br°; A₃ = H₂;  A₄ = H°, A₅ = HBr. Y escribiendo las reacciones en notación generalizada convencional,

c)   Representando la matriz de coeficientes (Usando MATLAB)

>> A = [-1 2 0 0 0; 0 -1 -1 1 1; -1 1 0 -1 1; 0 1 1 -1 -1; 1 -2 0 0 0]

A =

-1 2 0 0 0

0 -1 -1 1 1

-1 1 0 -1 1

0 1 1 -1 -1

1 -2 0 0 0

d)   Aplicando el método de eliminación de Gauss

>> B = rref(A)

B =

1 0 0 2 -2

0 1 0 1 -1

0 0 1 -2 0

0 0 0 0 0

0 0 0 0 0

e)             En el mecanismo existen tres reacciones que pueden ser la etapa controlante de la reacción:  las etapas son dadas por las reacciones 1, 2 y 3 del mecanismo propuesto. Por lo tanto se debe hacer un análisis mas profundo y usando técnicas experimentales para definir la etapa controlante entre las tres reacciones dadas.

f)             Esto también se verifica aplicando el rango

» C = rank(A)

C =

3               (Tres reacciones independientes)

También se puede usar la función  rrefmovie para ver paso a paso la tranasformación

rrefmovie (A)

Original matrix

A =

-1     2    0    0    0

0    -1   -1    1    1

-1     1    0   -1    1

0     1    1   -1   -1

1    -2    0    0    0

Press any key to continue. . .

pivot = A(1,1)

A =

1   -2    0    0    0

0   -1   -1    1    1

-1    1    0   -1    1

0    1    1   -1   -1

1   -2    0    0    0

Press any key to continue. . .

eliminate in column 1

A =

1   -2    0    0    0

0   -1   -1    1    1

-1    1    0   -1    1

0    1    1   -1   -1

1   -2    0    0    0

Press any key to continue. . .

Luego de múltiples manipulaciones llagamos a la solución final

B =

1    0    0    2   -2

0    1    0    1   -1

0    0    1   -2    0

0    0    0    0    0

0    0    0    0    0

4.4 SOLUCIONES DE PRUEBA Y ERROR

Los procesos reales en ingeniería química generalmente son descritos por ecuaciones no lineales. La solución de problemas de modelos no lineales casi siempre exige al uso de un procedimiento de prueba y error para encontrar una aproximación cercana a una raíz de una ecuación. Las raíces de polinomios pueden encontrarse con la función MATLAB roots.

Ejemplo 4.3  Un proyecto en ingeniería química requiere que se calcule exactamente el volumen molar (v) del dióxido de carbono y del oxígeno para diferentes combinaciones de temperatura y presión de tal forma que los recipientes que contengan estos gases se puedan seleccionar apropiadamente. Para tal efecto se sugiere usar la ecuación de van der Waals

Solución

La ecuación de van der Waals está dada por:

donde v = V/n es el volumen molar y a y b son constantes empíricas que dependen del gas en particular

	Dióxido de carbono	Oxígeno
A	3.592	1.360
B	0.04267	0.03183
R	0.082054 Lit atm/(mol K)

Desarrollando la Ec. 4.3 tenemos un polinomio en v de la forma:

pv³ – (bp + RT)v² + av – ab = 0                                                          (4.4)

Reemplazando los valores de a y b para el dióxido de carbono se tiene:

pv³ – (0.0427p + 0.082054 T)v² + 3.592v – 0.15327064 = 0             (4.5)

Ahora podemos usar la Ec. 4.5 para cualquier juego de valores de  P y  T. En este caso usaremos P = 1 atm y T = 300 K, con lo cual tenemos la ecuación  cúbica siguiente:

v³ – 24.6589v² + 3.592v – 0.15327064 = 0

Para usar Matlab, ingresamos los coeficientes del polinomio:

>> p=[1 -24.6589 3.592 -0.15327064]

Luego la orden:

>> r = roots(p)

Con lo cual se obtiene un vector columna con las raíces del polinomio:

r =

24.5126

0.0731 + 0.0301i

0.0731 – 0.0301i

>>

Con lo cual, el volumen molar del dióxido de carbono a P = 1 atm y T = 300 K es 24.5126 lit/mol

4.5 MÉTODOS ITERATIVOS DE CONVERGENCIA

Uno de los problemas más comunes en simulación digital es la solución de ecuaciones algebraicas no lineales simultáneas. Si estas ecuaciones contienen funciones trascendentales, la solución analítica es imposible. Entonces debe usarse algún tipo de procedimiento iterativo de prueba y error. Si hay una sola incógnita se supone un valor para la solución. Este es reemplazado en la ecuación o ecuaciones para ver si satisface. Si no, se hace un nuevo supuesto y el procedimiento es repetido hasta conseguir la convergencia (el valor esperado) no el valor de la derecha.

La clave del problema es encontrar un método para hacer  que el nuevo supuesto converja rápidamente  a la respuesta correcta. Existen varias técnicas. Desdichadamente no hay un mejor método para todas las ecuaciones. Algunos métodos que convergen muy rápidamente para algunas ecuaciones divergen para otras ecuaciones; por ejemplo, la serie de nuevos valores oscilará alrededor de la solución correcta con un incremento en las desviaciones. Esta es una forma de inestabilidad numérica.

Nosotros discutiremos solamente algunos de los métodos más simples y más usados. Afortunadamente, en simulación dinámica, partimos de algún estado inicial de convergencia. A cada instante en el tiempo, las variables son cambiadas pero en rangos pequeños de los valores que tenían en un instante antes.  Así, siempre estaremos cerca de la solución correcta. Por esta razón, los métodos simples de convergencia son muy usados para simulaciones dinámicas.

El problema se entiende mejor con un ejemplo. Uno de los cálculos iterativos más comunes es un calculo del punto de burbuja en el equilibrio liquido vapor.

Ejemplo 4.4.  Se tiene la presión P y la composición x.  Queremos encontrar la temperatura del punto de burbuja y la composición del vapor como se ha discutido en la Sec. 2.2.6. Por simplicidad asumimos un sistema binario de componentes 1 y 2. El componente 1 es el más volátil, y la fracción molar de componente 1 en el liquido es x y en el vapor es y. También asumimos que el sistema es ideal. Se aplican las Leyes de Dalton y Raoult.

Las presiones parciales de dos componentes (p₁ y p₂) en las fases de vapor y liquido son:

Donde  P^s_j es la presión de vapor del componente puro j la cual es una función de la temperatura solamente y esta dada por la Ecuación de Antoine:

Las presiones parciales en las fases liquido y vapor están dadas por:

nuestro problema de convergencia es para encontrar el valor de la temperatura T que satisfaga la Ec. (4.9). El procedimiento es como sigue:

1.      Suponer una temperatura T.

2.      Calcular la presión de vapor de los componentes 1 y 2 con las Ec. (4.8).

3.      Calcular la presión total P^calc usando la Ec. (4.9).

4.      Comparar P^calc con la presión total actual dada, P. Si es suficientemente cercana a P (por ejemplo usando un criterio de convergencia relativo de 10^-6), la T supuesta es correcta. Luego la composición del vapor puede ser calculada a partir de la Ec. (4.10).

5.      Si P^calc es mucho mayor que P, el valor supuesto de la temperatura fue muy alta y debe hacerse otro supuesto de T que sea mas bajo. Si P^calc es muy bajo, debemos suponer un valor de T mas alto.

4.5.1 Método de la bisectriz

Conocido también como corte binario, de partición en dos intervalos iguales o método de Bolzano es un método de búsqueda incremental en el que el intervalo se divide siempre en dos. Si la función cambia de signo sobre un intervalo, se  evalúa el valor de la función en el punto medio. La posición de la raíz se determina situándola en el punto medio del sub intervalo dentro del cual ocurre un cambio de signo. El proceso se repite hasta obtener una mejor aproximación.

Es  simple y fácil de visualizar y programar. No es muy rápida en convergencia a la solución correcta, pero es estable. Trabaja bien en simulaciones dinámicas debido a que el tamaño del paso se puede ajustar para corresponder aproximadamente a la velocidad ala cual es cambiada la variable con respecto al tiempo durante la integración.

Fig. 4-2  Convergencia de la Bisectriz

La Fig. 4-2 muestra gráficamente el procedimiento de la bisectriz. Se hace un supuesto inicial de temperatura T₀. Se calcula P^calc a partir de la Ec. (4.11). luego P^calc es comparado con P. Un incremento fijo en la temperatura DT es adicionado o sustraído de la temperatura supuesta, dependiendo de si P^calc es mayor o menor que P.

Nos seguimos moviendo en la dirección correcta con el tamaño de este paso fijo hasta que se de un cambio en el signo del término (P – P^calc). Esto nos indica que hemos cruzado al correcto valor de T. Entonces retornamos a la mitad del trayecto, por ejemplo partimos en dos el incremento DT. Con cada iteración sucesiva partimos en dos DT, yendo siempre por sobre o debajo de la temperatura.

4.5.2 Método de Newton Raphson

Este método es probablemente el método de convergencia más popular. Es algo más complicado ya que requiere la evaluación de una derivada. También puede traer problemas de estabilidad si el supuesto inicial es pobre y si la función es altamente no lineal.

Si el valor inicial de la raíz es x_i, entonces se puede extender una tangente desde el punto [x_i, f(x_i)]. El punto donde esta tangente cruza al eje x representa una aproximación mejorada de la raíz.

Usando el problema del punto de burbuja como un ejemplo específico, definiendo la función f_(T):

Queremos encontrar el valor de T que haga f_(T) igual a cero ( queremos encontrar la raíz de f_(T)). Suponemos un valor de temperatura T_o. Entonces evaluamos a T_o, f_(To).   A continuación evaluamos la función a T_o, f ^’_{(To) = (df /dT)(To).}

Fig. 4-4   Representación gráfica de la convergencia de Newton Raphson

Luego de la geometría mostrada en la Fig. 4-4 podemos ver que

Resolviendo para T₁ se tiene

T₁ en la Ec. (4.14) es el nuevo supuesto de temperatura. Si la curva f_(T) fuese una línea recta, podríamos converger a la solución correcta solamente en una iteración.

Generalizando la Ec. (4.14), podemos obtener el algoritmo de iteración recursivo.

donde  T_n+1 = nuevo supuesto de temperatura

T_n = valor anterior de temperatura supuesta

              F_n = valor de f_(T) a T = T_n

f ^’_n = valor de la derivada de f,  df/dT a  T = T_n

La técnica requiere la evaluación de f¢, la derivada de la función f_(T)  con respecto a la temperatura. En el ejemplo del punto de burbuja, esta puede ser obtenida analíticamente.

Si la función fuese tan compleja de tal manera que no podría obtenerse una derivada explícitamente, se podría calcular una derivada aproximada numéricamente: haciendo un pequeño cambio en la temperatura D T, evaluando f a  T + D T y usando la aproximación

Si la función no es muy lineal y/o si el supuesto inicial no está cercano a la solución, el método de Newton Raphson puede divergir en lugar de converger. Las funciones que no son monotónicas son particularmente molestosas, debido a que la aproximación a cero de la derivada, hace que la Ec. (4.15) sea inestable. Por lo tanto, el método de Newton Raphson es un algoritmo muy eficiente pero que tiene problemas de convergencia. A veces estas dificultades pueden ser salvadas reduciendo el valor del cambio permitido para hacer el nuevo supuesto.

El método de Newton Raphson puede  extenderse fácilmente a problemas de iteración involucrando más de una variable. Por ejemplo, suponer que se tiene dos funciones f_{1(X1, x2)}  y f_{2(X1, x2)} que dependen de las dos variables x₁ y x₂.  Debemos encontrar los valores de x₁ y x₂ que satisfagan las dos ecuaciones

f_{1(X1, x2)}  = 0                y          f_{2(X1, x2)} = 0

Expandiendo cada una de estas funciones alrededor del punto (x_1n,  x_2n ) en una serie de Taylor y truncando antes del término de la primera derivada se tiene

Haciendo f_{1(x1, n+1, x2, n+2)}  y  f_{2(x1, n+1, x2, n+2)}  iguales a cero y resolviendo para los nuevos supuestos x_{1, n+1}  y  x_2,n+1  se tiene

todas las funciones y derivadas parciales son evaluadas a x_1n y x_2n.

Las Ecs. (4.20) y (4.21) da el algoritmo de iteración para los nuevos valores supuestos. Para cada iteración deben calcularse cuatro derivadas parciales ya sea analítica o numéricamente.

Las Ecs. (4.18) y (4.19)  pueden escribirse en forma matricial. Haciendo el lado derecho de las ecuaciones iguales a cero y denominando a los cambios en las variables supuestas como D x = x_n+1 – x_n

Todas las funciones y derivadas son evaluadas a x_1n  y x_2n . La matriz 2 x 2 de derivadas parciales en llamada matriz jacobiana 

Todos los términos de esta matriz son constantes que son calculadas en cada iteración

Los D x’s se pueden calcular resolviendo la matriz

donde J^-1 es el inverso de la matriz J

4.5.3  Método de la falsa posición

Esta técnica de convergencia es una mezcla del método de la bisectriz y de Newton Raphson, se hace un supuesto inicial T_o, y se evalúa la función f_(To) se toma un paso en la dirección correcta para evaluar una nueva temperatura T₁ y se evalúa  f_(T1). Si f_(T1)  tiene el mismo signo de f_(To), no se ha cruzado la solución y se toma otro paso (redefiniendo T₁ como el nuevo valor de T_o). Los pasos son continuados hasta que se alcanza alguna temperatura T₁  donde f_(T1) difiera en el signo de f_(To). Como se muestra en la Fig. 4.5 mediante la geometría se puede suponer un nuevo valor de temperatura T₂. de triángulos similares

Rearreglando

Generalizando, conseguimos el algoritmo recursivo:

La Ec. (4.26) es similar a la Ec. (4.15), del algoritmo de Newton Raphson. La derivada en la Ec. (4.26) se aproxima numéricamente.

Fig. 4-5 Convergencia de la Falsa Posición

4.5.4  Métodos explícitos de convergencia

Para algunos sistemas de ecuaciones es posible suponer el valor de una variable x_sup, y luego usar una de las ecuaciones para resolver explícitamente para un nuevo valor calculado de la misma variable, x_calc. Lugo se comparan el valor supuesto con el valor calculado y se hace un nuevo supuesto.

El nuevo supuesto puede ser simplemente el valor calculado (esto se denomina sustitución sucesiva). La convergencia puede ser muy baja debido a:

1)      Muy poca aproximación entre x_calc y x_sup

2)      Una oscilación de x_calc de un lado a otro alrededor de x_sup. Puede divergir.

Entonces un factor de convergencia b puede usarse acelerar o reducir la velocidad a la cual se permite el cambio de x_sup de iteración a iteración.

(x_sup)_nueva = (x_sup)_ant + b[x_calc – (x_sup)_ant]                                 (4.27)

Notar que  b = 1 corresponde a sustitución sucesiva. Este método es ilustrado en el siguiente ejemplo.

Ejemplo 4.5  Un intercambiador de calor en contracorriente es un ejemplo importante de un sistema descrito por ecuaciones que son usualmente resueltas iterativamente. La Fig. 4-6 muestra el sistema. El problema es encontrar las temperaturas de salida al estado estacionario del aceite, T_H2. y del agua de enfriamiento, T_C2, y el calor transferido Q, dadas las temperaturas de entrada, las velocidades de flujo (caudales), y el coeficiente de transferencia de calor y el área. Las ecuaciones al estado estacionario para la transferencia de calor son:

Q = U A (DT)_LM = (120)(879)(DT)_LM                                             (4.28)

Q = (70 000)(0,5)(250 – T_H2)                                                           (4.29)

         Q = (170,5)(60)(8,33)(T_C2 – 80)                                                       (4.30)

Tenemos cuatro ecuaciones y cuatro variables: Q, (DT)_LM, T_H2 y T_C2. El procedimiento iterativo es:

1. Suponer un valor para la temperatura de salida del aceite T_H2 supuesta (la cual debe ser mayor que 80 ° por razones físicas)
2. calcular Q₁ a partir de la Ec. (4.29).
3. Calcular T_C2 a partir de la Ec. (4.30), usando Q₁.
4. calcular la fuerza impulsora, diferencia media logarítmica de temperaturas (DT)_LM   a partir de la Ec. (4.31).
5. Calcular una nueva velocidad de transferencia de calor Q₂ a partir de la Ec. (4.28).
6.  Sustituir el valor de Q₂ en la Ec. (4.29) y calcular T_H2 calculada.
7. Comparar T_H2 supuesta  y T_H2 calculada
8. Si no esta dentro de los limites de error fijados volver a suponer T_H2 usando la Ec. (4.27)

En este caso tomamos un valor de b = 0.2 y hemos iniciado la iteración con  T_H2 = 90 F

Fig. 4-6  Intercambiador de calor en contracorriente

4.8  INTEGRACIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES ORDINARIAS

Como se ha discutido en la introducción a este capítulo, la solución de ecuaciones diferenciales ordinarias (EDOs) en una computadora digital involucra la integración numérica. Presentaremos algunos de los algoritmos de integración numérica más simples y populares.

Cuando se resuelven EDOs se debe tener en mente la precisión, la estabilidad numérica y la velocidad de cualquier método de integración. Mostraremos algunas comparaciones entre los diferentes tipos de algoritmos en los problemas de este capitulo y el próximo.

4.8.1 Uso de las variables de estado

Cualquier ecuación de alto orden, por decir con derivadas de orden N, puede ser reducida a N  EDOs  de primer orden con lo cual solamente necesitamos el estudio de integración numérica de EDOs de primer orden.

Ejemplo 4.6

Suponer la EDO de tercer orden:

Si definimos las nuevas variables

La Ec. (4.32 ) se escribirá:

Luego tenemos para resolver tres EDOs de primer orden

o en forma matricial (y representando las derivadas con un punto en su parte superior para así facilitar la representación)

o en forma simbólica:

que recibe el nombre de ecuación de estado

u : es el valor por el cual se cambian las entradas

4.8.2        Algoritmos explícitos de integración numérica

Los algoritmos explícitos involucran cálculos explícitos de derivadas. Dos métodos muy usados son descritos a continuación: El método de Euler y  el método de Runge y Kutta de cuarto orden. Literalmente existen cientos de métodos pero no son muy usados particularmente para problemas de ingeniería química.

A. ALGORITMO DE EULER.  El esquema posiblemente más simple de integración numérica (y el mas usado) es el de Euler, ilustrado en la Fig.4-7. Asumiendo que deseamos resolver la EDO

Donde f_{(x, t)} es, en general una función no lineal. Necesitamos conocer donde empezar. Por ejemplo, necesitamos conocer una condición inicial para x, usualmente esta es al tiempo igual a cero.

            x(o) = x₀     a    t = 0                                                                     (4.40)

Ahora si nos movemos hacia delante en el tiempo por un pequeño incremento D t, hacia t = t₁ = D t, podemos obtener un estimado del nuevo valor de x para t = D t,  x_(Dt), mediante una extrapolación lineal usando la velocidad inicial del cambio de x con el tiempo (la derivada de x a t = 0). El nuevo valor de x (x₁) es aproximadamente igual al valor anterior de x (x₀) más el producto de la derivada de x por el tamaño de paso.

Si el tamaño de paso (intervalo de integración) es pequeño, este estimado de x₁ será muy cercano al valor correcto.

Tomando otro D t para  t = t₂ = 2 D t, estimamos x_(2 D t) = x₂ a partir de

Generalizando para el paso (n + 1) en el tiempo,

La integración de Euler es extremadamente simple para programar, como se ilustrará en el Ejemplo 4-9 esta simplicidad es mantenida, aún cuando el número de EDOs se incrementa y las derivadas de las funciones sean mas complejas y no lineales.

Si tenemos dos EDOs simultáneas acopladas para integrar numéricamente

El algoritmo de integración de Euler podría ser

x₁, _n+1  = x_1n + D t f_{1 (x1n, x2n, tn)}                                                           (4.47)

x₂, _n+1  = x_2n + D t f_{2 (x1n, x2n, tn)}                                                            (4.48)

Notar que solamente se requiere la evaluación de una derivada por cada EDO a cada punto en el tiempo. Si tenemos un sistema de N ecuaciones diferenciales, deberíamos tener N ecuaciones tales como las Ecs. (4.47) y (4.48)

B.  RUNGE Y KUTTA (CUARTO ORDEN). El algoritmo de Runge y Kutta de cuarto orden es extensamente usado en ingeniería química. Para la EDO dada en la Ec. (4.39) el algoritmo de Runge y Kutta es:

Para integración numérica de dos EDOs de primer orden

Con  Runge y Kutta de cuarto orden, se evalúan cuatro k’s para cada EDO

          k₁₁ = D t f_{1(xn1 ,  xn2,   tn)}

         k₁₂ = D t f_{2(xn1 ,  xn2,   tn)}

         k₂₁ = D t f_{1(xn1 + ½ k11,  xn2 + ½ k12 ,   tn + ½ D t)}

         k₂₂ = D t f_{2(xn1 + ½ k11,  xn2 + ½ k12 ,   tn + ½ D t)}

        k₃₁ = D t f_{1(xn1 + ½ k21,  xn2 + ½ k22 ,   tn + ½ D t)}

        k₃₂ = D t f_{2(xn1 + ½ k21,  xn2 + ½ k22 ,   tn + ½ D t)}

        k₄₁ = D t f_{1(xn1 +  k31,  xn2 +  k32 ,   tn +  D t)}

        k₄₂ = D t f_{2(xn1 +  k31,  xn2 +  k32 ,   tn +  D t)}

Luego son calculados los nuevos valores de x₁ y x₂:

4.8.3 Métodos implícitos

Los métodos explícitos considerados en la sección previa involucran la evaluación de derivadas, seguida por el cálculo explicito de nuevos valores para las variables en el siguiente punto en el tiempo. Como el nombre implica, los métodos de integración implícita usan algoritmos que  dan como resultante ecuaciones que deben resolverse para los nuevos valores en el paso siguiente en el tiempo. Una simple EDO ilustra la idea.

El algoritmo implícito de Euler de primer orden es

Comparando esto cuidadosamente con el algoritmo explicito dado en la Ec. (4.44). La derivada es evaluada en el siguiente paso en el tiempo donde no conocemos la variable x_n+1. por lo tanto la incógnita x_n+1 aparece en ambos lados de la ecuación. Considerando la simple EDO

El algoritmo implícito de Euler es

             x_n+1 = x_n + (1 – x_n+1) Dt

Resolviendo para la incógnita x_n+1  se tiene

La principal ventaja del algoritmo implícito es que no se hace numéricamente inestable. Se pueden tomar tamaños de paso muy grandes sin peligro de problemas inestabilidad que si se dan en los métodos explícitos. Por lo tanto los métodos implícitos son muy usados para sistemas muy estables.

La Ec. (4.63) puede ser resuelta fácilmente para la incógnita x_n+1. Sin embargo, suponiendo que tenemos un número grande (N) de EDOs para integrar numéricamente. En general, las N derivadas dependen de todas las variables, dando N ecuaciones simultáneas (usualmente no lineales) que deben resolverse  a cada punto en el tiempo para los N valores no conocidos de las variables para el siguiente punto en el tiempo.

Considerando el conjunto de N  EDOs lineales

donde x es un vector de variables y A es una matriz cuadrada de constantes. Aplicando el algoritmo implícito de Euler a estas ecuaciones se tiene

        x_n+1 = x_n +[A x_n+1 D t]                                                                   (4.66)

Resolviendo para x_n+1 se tiene

     x_n+1 = x_n +[I – A D t] ^–1 x_n                                                                (4.67)

Donde el exponente –1 en el término significa “matriz inversa”. Esto nos lleva a deducir que el uso de algoritmos implícitos nos lleva al cálculo de matrices inversas (usualmente en cada punto en el tiempo debido a que las EDOs generalmente son no lineales).

4.9 LINEALIZACION

Como se ha mencionado anteriormente, debemos convertir las ecuaciones diferenciales no lineales rigurosas que describen el sistema en ecuaciones diferenciales lineales para poder usar las técnicas matemáticas lineales.

La primera interrogante a ser respondida es cuando una ecuación diferencial es lineal?. Básicamente esta es una que contenga variables elevadas solamente a la primera potencia en cualquiera de sus términos de la ecuación. Si aparecen en la ecuación raíces cuadradas, cuadrados, exponenciales, productos de variables, etc., esta es no lineal:

Ejemplo de lineal

donde a_o y a₁ son constantes o funciones del tiempo solamente, no de las variables dependientes o sus derivadas.

Ejemplo de no lineales

donde x₁ y x₂ son variables dependientes

Matemáticamente, una ecuación diferencial lineal es aquella para la cual se mantienen las siguientes dos propiedades.

1. Si  x_(t)  es una solución, entonces cx_(t) es también una solución, donde c es una constante.
2. Si  x₁  es una solución y x₂ es también una solución, entonces x₁ + x₂ es una solución.

La linealización se efectúa expandiendo las funciones no lineales mediante la serie de expansión de Taylor alrededor de la operación al estado estacionario y despreciando todos los términos siguientes a las primeras derivadas parciales.

Asumiendo que se tiene la función no lineal f de las variables de proceso x₁ y x₂: f_{(x1, x2)}. Por ejemplo, x₁ puede ser fracción molar, temperatura o velocidad de flujo. Nosotros denotaremos los valores al estado estacionario de estas variables usando el sub índice s:

x_1s º valor al estado estacionario de x₁

x_2s º valor al estado estacionario de x₂

expandimos la función f_{(x1, x2)}. Alrededor de sus valores al estado estacionario f_{(x1s, x2s)},

La linearización consiste en truncar la serie después de las primeras derivadas parciales

Hemos aproximado a la función real mediante una función lineal

Ejemplo 4.8

El flujo de salida de un tanque (F), depende de la raíz cuadrada del nivel de liquido (h), mediante la expresión:

La serie de expansión de Taylor alrededor del valor al estado estacionario de h, el cual es h_s en nuestra nomenclatura es

Ejemplo 4.9

La dependencia de la temperatura de la ley de Arrhenius sobre la velocidad específica de reacción k es una función altamente no lineal, la cual es linealizada como sigue:

donde  ks º k_(Ts)

Ejemplo 4.10

El producto de dos variables dependientes es una función no lineal de las dos variables:

f_{(CA, F)} = C_A F                                                                               (4.79)

Linearizando:

Notar que el proceso de linealización convierte a la función no lineal (el producto de dos variables dependientes) en una función lineal conteniendo dos términos.

Ejemplo 4.11

Considerar la ecuación diferencial ordinaria no lineal para el tanque de flujo por gravedad del Ejemplo 2.9

Linealizando el término   v²  se tiene:

v² =  + (2v_s)(v – v_s)                                                                    (4.83)

Luego la Ec. (4.82) se transforma en:

Esta EDO es ahora lineal. Los términos entre paréntesis son constantes, los cuales dependen, desde luego, del estado estacionario alrededor del cual es linealizado el sistema.

Ejemplo 4.12

La ecuación de continuidad por componente para una reacción irreversible no isotérmica de n-esimo orden ocurriendo a   volumen constante en un CSTR es:

Lineariando se tiene:

Hasta ahora hemos visto ejemplos donde la no-linealidad está en los términos derivados, es decir, en el lado derecho de la EDO. Muy a menudo el modelo de un sistema puede ser dado en una EDO la cual contenga términos no lineales dentro  del término de la derivada. Por ejemplo suponer que el modelo de un sistema no lineal es:

el procedimiento correcto para linearizar este tipo de ecuación es rearreglarlo de tal manera que las funciones no lineales aparezcan solamente en el lado derecho de la EDO y luego linearizarlo en la forma normal, por ejemplo dada la Ec. (4.87), diferenciamos el término h³ para conseguir

Luego rearreglando se tiene:

ahora estamos listos para linearizar.

Esta es una EDO lineal con coeficientes constantes