ECUACIONES DIFERENCIALES ADIMENSIONALIZADAS Y NÚMEROS ADIMENSIONALES

Introducción

El artículo que hoy me permito presentarles me parece interesante y pertinente porque en los magníficos libros de que dispongo lo más común es que las ecuaciones se presenten ya adimensionalizadas, y que no se expliqué qué hacer con ellas; ni de qué manera, y hasta qué extremo pueden ser útiles ni cuales son las limitaciones de los resultados de una integración numérica de éste tipo.

Para los fines consiguientes tomaré el caso de conducción unidimensional en un sólido semiinfinito, que es muy conocida y que por esa razón no voy a deducir, porque ya lo he hecho en algún artículo anterior, y en mi libro de Control Automático.

La ecuación es:

Figura 1. Ecuación de transferencai de calor unidimensional transitoria

Figura 1. Ecuación de transferencia de calor unidimensional transitoria

Si se decide adimensionalizarla, lo que la vuelve relativamente flexible, se debe decidir la manera en que se ha de definir la variable adimensional que re presente a la temperatura, la dirección de conducción, que debe ser el espesor del sólido Siendo el tiempo adimensional el Número de Fourier (αt/W^2).

En nuestro caso la pared izquierda del sólido se calienta a una temperatura T0 a tiempo cero, lo que en la práctica no es posible y adopté la siguiente expresión que es la que utilizan casi todos los autores, que es la siguiente en la que Ti es la temperatura inicial de todo el sólido antes que la perturbación que sufre en su pared izquierda y, a la vez, es la temperatura del medio que circunda al sólido.

T es el valor de la temperatura que va variando con el tiempo t y en el sentido x y T0 es el valor de la temperatura que adquiere la pared izquierda del sólido a tiempo cero. La perturbación en la pared izquierda del sólido produce su calentamiento, desde el lado más caliente (pared izquierda) hacia el lado más frío que es todo el sólido por virtud de la Segunda Ley de la Termodinámica (J. Smith y H.C. Van Ness 1959).

La Expresión para la temperatura adimensional es:

Figura 2. Definición de la temperatura adimensional

Figura 2. Definición de la temperatura adimensional

Para resolver la ecuación numéricamente, debemos elegir las células de cálculo. En este caso la célula que representa la derivada parcial respecto del tiempo es una célula de cálculo hacia delante (Ames 1969), lo que significa que –en cierta manera “viajamos hacia el futuro”.

Figura 3. Células de calculo para la integración numérica

Figura 3. Células de calculo para la integración numérica

En la fórmula anterior el subíndice i denota posición y el subíndice j denota el tiempo. El significado práctico de la utilización de las células ( de temperatura la una, y de tiempo la otra) es que teniendo tres valores de temperatura a tiempo j, es posible calcular el valor de la temperatura central a tiempo j+1 (por eso la denominación de “calculo hacia delante”).

Como el cálculo se puede hacer en Excel™, en ese caso cada fila subsiguiente (hacia abajo) corresponde a un tiempo j+1, j+2, etc.

El cálculo en Excel™ es muy práctico porque se puede usar para probar el algoritmo sin perjuicio de que aún así, si el intervalo de tiempo es muy pequeño, por ejemplo si Δt=0.001s, se requieren 1000 líneas para avanzar 1 s en la integración y la labor se vuelve muy estresante, “manual” e impráctica, debido a que, si la mano no se mueve correctamente, la integración se malogra debiéndose ubicar el sitio donde se descarrió, para volver a empezar desde allí.

La expresión anterior consta de dos células de cálculo en donde la temperatura se expresa en grados Celsius, Fahrenheit o Kelvin, la posición se expresa en metros, y el tiempo se expresa en segundos; es decir las células son dimensionales, de la manera en que han sido escritas y l -en este caso particular- queremos que sean adimensionales.

Para estos fines se procede a escribir la expresión anterior de manera que la variable temperatura se transforme en una variable adimensional, para lo que se procede a aplicar la definición que mencionamos anteriormente a cada una de las variables dimensionales de cada miembro de la ecuación anterior, dividiendo cada variable dimensional para T0-Ti, y restándole Ti (no confundir la “i” de la temperatura inicial con la “i” de posición).

Si se lleva a cabo lo que queda descrito, se obtiene la siguiente expresión que está escrita en función de las temperaturas adimensionales, que se reemplazan por el símbolo θ para obtener la expresión adimensional, lo que se muestra en la ecuación subsiguiente.

Figura 4. Ecuación de adimensionalización

Figura 4. Ecuación de adimensionalización

Si procedemos entonces a escribir la ecuación anterior mediante la notación indicada obtendremos la expresión siguiente despejando Ti,j+1i (método explícito (Ames 1969).

Figura 6. Ecuación adimensionalizada

Figura 5. Ecuación adimensionalizada

Y aquí vemos que el primer paréntesis aparece el módulo incremental de Fourier (J. Alan Adams y David F. Rogers 1973), que es adimensional y que, si se le quiere ver de alguna forma asimilable al entendimiento matemático común, puede –salvo mejores criterios, interpretarse como un un Δt adimensional.

Como decidí que la pared derecha sea convectiva, necesitamos una ecuación especial para el nodo superficial de ella.

Para obtener la ecuación de interés escribimos la ecuación de convección en términos de las variables adimensionales. Como yo dividí en cuerpo en diez espacios iguales entonces escribí las ecuaciones siguientes con sus subíndices reales de posición, y son las siguientes.

Figura 7. Condición de borde convectiva

Figura 6. Condición de borde convectiva

Obsérvese que los nodos Θ11,j y Θ9,j son temperaturas adimensionales separadas por dos espacios Δx en la malla de cálculo, y que esta es la razón para dividir su diferencia por 2Δx.

También es pertinente observar que las temperaturas que sirven para construir esta condición de borde ocurren todas al mismo tiempo j. Igualmente, es muy importante notar que el nodo Θ11,j es ficticio, por estar realmente fuera del sólido.

Para que el nodo ficticio desaparezca se procede a combinar la condición de borde anterior con la ecuación de cálculo para Θ10,j+1 que incluya el nodo ficticio.

La ecuación en cuestión es la siguiente:

Figura 8. Ecuación del nodo convectivo

Figura 7. Ecuación del nodo convectivo

No es necesario el reemplazo algébrico ya que el computador calcula primero la una expresión y automáticamente la reemplaza en la segunda.

Esta segunda ecuación es igual a las que se utilizan para el cálculo de los nodos interiores apareciendo en ella el módulo de Fourier.

Es importante resaltar que en la primera ecuación, que es una condición de borde aparece el Número de Biot, hΔx/k, que –adimensional y todo- contiene un parámetro característico del material y otro de su “circunstancia”, como diría Ortega y Gasset, que son, en su orden, el coeficiente de conductividad térmica k y el coeficiente de convección de calor h.

El Número de Biot es exactamente igual en estructura al Número de Nusselt, pero su aplicación y significado son diferentes, debiéndose esto a que el Numero de Biot (Bi) se aplica entre dos fases (en este caso entre el sólido y el aire que lo circunda), y el Número de Nusselt se aplica en una sola fase (caso de transferencia de calor en un reactor tanque-agitado, por ejemplo.

Una vez dilucidado todo lo que queda indicado, resta mostrar el programa, que inevitablemente debe contener variables subindicadas, que –como sabemos- en programación se conocen como “arreglos” y en matemática se conocen como vectores o como matrices, según sea el caso.

El Programa

A continuación se muestra el código del programa VBA, que corre muy rápidamente.

Option Base 0

Dim T(13, 13), Deltaprint, Tiempo, alfa, Deltat, Deltax As Double

Dim TT As Single

Dim jjj As Integer

Sub ecuaciones4()

‘Ingreso de variables operacionales

Ti = Worksheets(“Programa”).Cells(2, 3)

T0 = Worksheets(“Programa”).Cells(3, 3) ‘Temperatura sobre la cara izquierda del sólido a Fo cero

Deltat = Worksheets(“Programa”).Cells(4, 3)

k = Worksheets(“Programa”).Cells(5, 3)

Ro = Worksheets(“Programa”).Cells(6, 3)

Cp = Worksheets(“Programa”).Cells(7, 3)

Deltax = Worksheets(“Programa”).Cells(8, 3)

h = Worksheets(“Programa”).Cells(9, 3)

W = Worksheets(“Programa”).Cells(10, 3)

‘ffo = Worksheets(“Programa”).Cells(12, 3)

alfa = k / (Ro * Cp)

ffo = alfa * Deltat / Deltax ^ 2

W = Worksheets(“Programa”).Cells(10, 3) ‘Espesor del cuerpo semiinfinito

Deltax_bar = Deltax / W

DDeltaprint = 0.1

‘Enceramiento de temperaturas

For i = 0 To 10

For j = 0 To 10

T(i, j) = 0

T(i, j) = 0

Next j

Next i

Worksheets(“Programa”).Range(“A14:n500”).Clear

Call formateo

Fila = 16

‘Condición necesaria para convergencia del incremento de Fo, fo en CAHTA

Worksheets(“Programa”).Range(“A14:A14”).Clear

Bi = h * W / k

Worksheets(“Programa”).Cells(Fila, 1) = “Número de Biot= ” & Bi

Call Escritura_de_Titulos(Fila, Deltax, W)

principio = True

Incremento = Deltax_bar ^ 2 * ffo

verificacion = True

otra_vuelta:

Deltaprint = 0.1

Worksheets(“programa”).Range(“a14:a14”).Clear

For Fo = 0 To 1 Step Incremento

If Fo = 0 And principio = True Or Fo >= Deltaprint Then

principio = False

Call Escritura_de_Resultados(T, Fo, Deltax, Deltat, alfa, Fila, T0, Ti, Deltaprint, h, TTmax, k, W, TT, Bi, Incremento, jjj, Deltax_bar, DDeltaprint, ffo)

Else

Call Calculo_de_temperaturas(alfa, Deltat, Deltax, Deltax_bar, L, h, Ti, Tini_ref, Tcond_borde, T, k, W, Bi, Incremento, TT, Fo, T0, ffo)

End If

Next Fo

Call Escritura_de_Resultados(T, Fo, Deltax, Deltat, alfa, Fila, T0, Ti, Deltaprint, h, TTmax, k, W, TT, Bi, Incremento, jjj, Deltax_bar, DDeltaprint, ffo)

End Sub

Sub Escritura_de_Titulos(Fila, Deltax, W)

‘Worksheets(“Programa”).Cells(Fila, 1) = “Tiempo, s”

Worksheets(“Programa”).Cells(Fila, 2) = “Fo”

Call formateo

For j = 3 To 13

Worksheets(“Programa”).Cells(Fila, j) = ((j – 3)) * (Deltax) / W ‘Imprime la localización de los espacios

Next j

Fila = Fila + 1

End Sub

Sub Escritura_de_Resultados(T, Fo, Deltax, Deltat, alfa, Fila, T0, Ti, Deltaprint, h, TTmax, k, W, TT, Bi, Incremento, jjj, Deltax_bar, DDeltaprint, ffo)

If Fo = 0 Then

‘Worksheets(“Programa”).Cells(Fila, 1) = TT

Worksheets(“Programa”).Cells(Fila, 2) = Fo

For i = 0 To 10

If i = 0 Then

T(i, 1) = (T0 – Ti) / (T0 – Ti)

Worksheets(“Programa”).Cells(Fila, 3) = T(i, 1)

Else

Worksheets(“Programa”).Cells(Fila, i + 3) = T(i, 1)

End If

Next i

Fila = Fila + 1

Deltaprint = Deltaprint + DDeltaprint

Else

‘TT = Fo * W ^ 2 / alfa

‘Worksheets(“Programa”).Cells(Fila, 1) = TT

Worksheets(“Programa”).Cells(Fila, 2) = Fo

If i = 0 Then

T(0, 2) = (T0 – Ti) / (T0 – Ti)

Worksheets(“Programa”).Cells(Fila, 3) = T(0, 2)

End If

‘If i = 0 Then

‘T(0, 2) = (T0 – Ti) / (T0 – Ti)

‘End If

For i = 1 To 9

T(i, 2) = ffo * (T(i + 1, 1) + T(i – 1, 1)) + T(i, 1) * (1 – 2 * ffo)

Worksheets(“Programa”).Cells(Fila, i + 3) = T(i, 2)

Next i

T(11, 1) = T(9, 1) – 2 * Bi * Deltax_bar * T(10, 1)

T(10, 2) = ffo * (T(11, 1) + T(9, 1)) + T(10, 1) * (1 – 2 * ffo)

Worksheets(“Programa”).Cells(Fila, 13) = T(10, 2)

Fila = Fila + 1

Deltaprint = Deltaprint + DDeltaprint

End If

End Sub

Sub Calculo_de_temperaturas(alfa, Deltat, Deltax, Deltax_bar, L, h, Ti, Tini_ref, Tcond_borde, T, k, W, Bi, Incremento, TT, Fo, T0, ffo)

If i = 0 Then

T(0, 2) = (T0 – Ti) / (T0 – Ti)

End If

For i = 1 To 9

T(i, 2) = ffo * (T(i + 1, 1) + T(i – 1, 1)) + T(i, 1) * (1 – 2 * ffo)

Next i

T(11, 1) = T(9, 1) – 2 * Bi * Deltax_bar * T(10, 1)

T(10, 2) = ffo * (T(11, 1) + T(9, 1)) + T(10, 1) * (1 – 2 * ffo)

‘Se colocan las temperaturas a tiempo 2 en las temperaturas 1 con excepción de la sub 9

‘para que sea posible el próximo cálculo

For jj = 1 To 10

T(jj, 1) = T(jj, 2)

Next jj

‘TT = (Fo * W ^ 2) / alfa

End Sub

Interfase de usuario

A continuación se muestra la interface de usuario

Figura 9. Interfase de usuario

Figura 8. Interfase de usuario

Comentarios al Programa

El primer aspecto que se debe resaltar es el de Option Base 0, que permite que existan subíndices cero.

En las respectivas consultas que hice en Internet respecto de este asunto se indicaba que esa declaración no era necesaria, porque -por defecto- el VBA acepta subíndices cero (como Apple™ realiza –gracias a Dios- tantas actualizaciones a su Office™ es posible que esto sea cierto, pero preferí la senda segura, para no tener qué preocuparme de aquello.

También es interesante observar que el programa en sí es muy corto, y para resaltar esto lo he puesto en caracteres verdes.

Esta brevedad en el programa se debe al uso de subrutinas, que facilita y orden de la escritura, estructura el programa, y -como cada subrutina lleva a cabo una sola tarea el chequeo del programa se facilita muchísimo porque, porque es fácil encontrar el origen del error o anomalía.

Lo dicho se facilita aún más porque se puede utilizar la opción “Add Watch”, que con su ventana respectiva permite chequear, en tiempo real, el valor de cada variable cada vez que se ejecuta cada línea ejecutando el programa mediante la combinación (ante le llamaban “acorde” fn+F5”, en Mac.

Para que el gráfico sea completamente adimensional, las temperaturas y la posición deben variar entre cero y uno.

Resultados

Para que los resultados sean utilizables se debe transponer la matriz de resultados que se obtiene. Cuando se ejecuta el programa que en Excel™ para Mac esto es sumamente fácil ya que solamente se requiere hacer un “Copy” y un “Paste Special” en donde simplemente se escoge la opción “Transpose”

Para clarificar la explicación les incluyo una de las matrices resultado, que –como queda enunciado- debe transponerse para poder generar el gráfico.

A continuación se muestra la matriz resultado.

matriz-resultado

Figura 9. Matriz Resultado

Seguidamente se muestra la matriz transpuesta

Figura 10. Matriz transpuesta

Figura 10. Matriz transpuesta

Como les comentaba, si se quiere que la posición (x/W) sr grafique en el eje de abscisas del gráfico, siendo W el espesor a través del que se conduce el calor; y si se desea que en el eje de ordenadas aparezca la temperatura adimensional ((T-Ti)/(T0-Ti)), la matriz original debe transponerse.

En ésta matriz transpuesta se puede apreciar que la primera columna representa los tiempos adimensionales de Fourier, y que debajo de cada uno de ellos están colocadas las temperaturas adimensionales que corresponden a ese tiempo Fo={0, 0.1, o.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0}

Lo dicho significa que, al graficar cada columna se obtiene el trazo que corresponde a cada tiempo que –a su vez- se relaciona, en cada uno de sus puntos, con las temperaturas adimensionales graficadas en abscisas, entre 0 y 1..

Dispuesto el gráfico de esta manera, cualquier línea perpendicular al eje de abscisas intersecará a una línea de tiempo de Fourier, trazada o interpolada gráficamente, y la intersección permitirá obtener el valor de la temperatura adimensional respectiva.

Si se invoca el principio Cartesiano ( a cada punto de un plano corresponde un valor y vice-versa), se puede decir que a cada punto del gráfico corresponde una posición, un tiempo y una temperatura.

A continuación se muestra el resultado gráfico de la distribución de temperaturas para las condiciones de borde: isotérmica para la cara izquierda y adiabática para la cara derecha (h=0 w/m2-oC).

Figura 11. Resultados para condición de borde adiabática en cara derecha

Figura 11. Resultados para condición de borde adiabática en cara derecha

Figura 11. Resultados para condición de bordee convectiva en cara derecha

Figura 12. Resultados para condición de borde convectiva en cara derecha

Comentario General

Es mi parecer que estos gráficos son muy útiles para determinar el tiempo que discurre entre que la pared isotérmica adquiere su temperatura, de manera “instantánea” (la que sea) y la temperatura llegue a un valor dado, para SS-304 con el coeficiente de convección indicado, porque –si usted no puede cambiar los dos parámetros del número de Biot h, y k, que se han usado para la integración numérica de las ecuaciones que están en el programa, no hay mucho más que usted pueda hacer, salvo mejores y más ilustrados criterios.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Acerca de

Professor of modeling and simulation, and process design at Escuela Politécnica Nacional, in Quito, Ecuador. . In the past I was a P4, P5, and D1 at the Organization for the Prohibition of Chemical Weapons, located in the Kingdom pf the Netherlands

Publicado en Sin categoría

Responder

Introduce tus datos o haz clic en un icono para iniciar sesión:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s

Enter your email address to follow this blog and receive notifications of new posts by email.

Únete a otros 1.775 seguidores

Categorías
Artículos y comentarios sobre modelado y simulación de plantas, y equipos de la industria química, con ejemplos
PREGUNTAS O INQUIETUDES
PUEDEN ENVIAR SUS PREGUNTAS/INQUIETUDES RESPECTO DE ARTÍCULOS DEL BLOG, O TEMAS AFINES A gasteaux@hotmail.com ASUNTO: BLOG DE INGENIERÍA QUÍMICA
A %d blogueros les gusta esto: