Simulación de estados transitorios en mecánica de fluidos

Introducción

La simulación numérica de procesos físicos se ha convertido actualmente, gracias a los avances de los ordenadores, en una herramienta de diseño de gran utilidad para el ingeniero.

Las técnicas de simulación numérica presentan ciertas ventajas respecto a los ensayos experimentales, como gran flexibilidad en las variables a medir y en los cambios de la geometría del problema, inversión económica reducida, resultados en un plazo corto de tiempo, etc.

Sin embargo, los resultados conseguidos mediante simulaciones numéricas deben ser interpretados siempre con cautela, siendo conscientes de que su exactitud dependerá en gran medida del modelo matemático que representa el proceso físico y de la técnica numérica empleada para resolver las ecuaciones. Las predicciones de simulaciones numéricas deberán ser siempre contrastadas con resultados experimentales.

Objetivos de la práctica El objetivo de la práctica es resolver numéricamente varios problemas de Mecánica de Fluidos mediante una hoja de cálculo.

La hoja de cálculo que se va a manejar en la práctica resuelve las ecuaciones diferenciales que describen los fenómenos transitorios a partir de una serie de simplificaciones que evitan el uso de técnicas matemáticas numéricas complejas, no abordables mediante una hoja de cálculo.

Los objetivos de esta práctica son los siguientes:

Estudiar la influencia de los parámetros de distintos problemas de mecánica de fluidos mediante simulaciones numéricas.

Aprender a escribir las ecuaciones de un problema de mecánica de fluidos de forma que se pueda realizar una simulación numérica.

Los casos que se van a estudiar son:

Estado estacionario en tuberías

Oscilaciones en un tubo en U

Flujo por gravedad

Cálculo de chimeneas de equilibrio

Con el objeto de no saturar el guión de prácticas, el desarrollo teórico de cada apartado se incluye en un capítulo final denominado Ecuaciones.

Estado estacionario en tuberías

Descripción del problema En este apartado se va a calcular el tiempo que tarda un flujo en alcanzar un régimen estacionario en una tubería cuando se abre una válvula repentinamente al final del conducto. Para este cálculo se tendrán en cuenta las pérdidas singulares y las de fricción.

Inmediatamente después de la apertura de la válvula, la altura del depósito está disponible para que, mediante la presión hidrostática, empiece a fluir el fluido. A medida que la velocidad aumenta también lo hacen las pérdidas singulares y por fricción, disminuyendo así la carga original disponible.

El depósito tiene en su interior una altura de líquido H y descarga a través de una tubería horizontal de longitud L cuyo coeficiente de pérdidas es K. La simulación mediante este modelo nos permite seguir la evolución del transitorio hasta la obtención de un estado pseudoestacionario del fluido circulando por la tubería de descarga (no es realmente estacionario puesto que el caudal de descarga es función de la altura de líquido en el depósito).

Se supondrá que el volumen de líquido en el tanque es lo suficientemente grande como para que se pueda suponer que la altura de líquido en su interior permanece constante durante el transitorio. Contesta a las siguientes preguntas La hoja de cálculo permite introducir valores de ciertas variables del problema (celdas amarillas) y observar cómo afectan estos cambios al fenómeno físico estudiado.

1) Indicar cuáles son las variables independientes (parámetros) y las variables dependientes (resultados) del problema.

2) Realizar un análisis paramétrico del problema. Este estudio debe incluir gráficas que muestren la influencia de una variable sobre el fenómeno físico y comentarios sobre las mismas.

Oscilaciones en un tubo en U

Descripción del problema En este apartado se simulará lo que ocurre en un tubo en U cuando el líquido en su interior es desplazado de su posición de equilibrio. Existen tres casos interesantes de oscilaciones de un líquido en un tubo en U:

Líquido sin fricción

Líquido con fricción, en régimen laminar

Líquido con fricción, en régimen turbulento

No se tratará el último de estos tres casos debido a la complejidad de las ecuaciones que gobiernan este fenómeno, imposibles de resolver mediante una hoja de cálculo simple.

Las variables del problema son: longitud total del líquido dentro del tubo (L), posición inicial de la superficie libre de la rama derecha (zo) y la velocidad inicial de dicho punto (vo). El origen de coordenadas z está en la posición de equilibrio y z toma valores positivos hacia arriba (ver figura).

El flujo transitorio en el interior de este tubo se tratará como flujo unidimensional a lo largo del recorrido del tubo. Se emplearán las ecuaciones de continuidad y de cantidad de movimiento para flujos no permanentes.

Oscilaciones sin fricción

La hoja de cálculo tiene tres campos amarillos en los que se deben introducir los datos del problema: longitud de líquido (L), altura inicial (zo) y velocidad inicial (vo) Una vez introducidos los datos, la hoja calcula las siguientes variables: altura máxima (zmax), velocidad máxima del fluido en el tubo (vmax) y el periodo de oscilación (T) (campos marcados como "Resultados de la simulación"). Además, se presenta una gráfica de evolución de la altura del líquido respecto a la posición de equilibrio y otra superpuesta, correspondiente a velocidad del fluido, ambas respecto al tiempo.

Contesta a las siguientes preguntas

1) Introduce los siguientes valores: L=5, zo=2 y vo=0. Este caso significa que hemos elevado 2 m la rama derecha y que la soltamos sin imprimirle velocidad. ¿Qué valores obtienes para zmax y para T? Observa las gráficas de evolución de z y v. ¿Es lógico que zmax sea igual a zo? ¿Por qué?

2) ¿Qué tipo de movimiento se obtiene? ¿Cuánto tiempo se tardará en alcanzar el estado estacionario?

3) ¿Qué valor toma z cuando v es cero? ¿Y qué le ocurre a v cuando pasa por la posición de equilibrio (z=0)? Razona estos resultados.

4) Realiza un análisis paramétrico de la influencia de los parámetros del problema sobre zmax, y sobre T. Este análisis debe incluir gráficas y comentarios.

5) Discute sobre si los siguientes factores pueden afectar al resultado de la simulación: anchura del tubo, densidad del líquido y radio de curvatura del tubo en U.

Oscilaciones con fricción

La hoja de cálculo tiene cuatro campos amarillos en los que se deben introducir los datos del problema: longitud de líquido (L), altura inicial (zo), velocidad inicial (vo) y viscosidad cinemática del líquido (().

Una vez introducidos los datos, la hoja calcula el parámetro m (ver hoja de cálculo), el cual permite saber si dominan las fuerzas viscosas sobre las inerciales. También se calcula la amplitud máxima de la oscilación (Z).

Finalmente, se presenta una gráfica de evolución de la altura del líquido respecto a la posición de equilibrio y otra correspondiente a velocidad del fluido, ambas respecto al tiempo.

Contesta a las siguientes preguntas

1) Analizar las oscilaciones para la siguiente situación: L=20 cm, zo=5 cm, vo=2 m/s y (=1.08e-6 m2/s (agua), utilizando tres diámetros internos: D=1 cm, 5 mm y 1 mm.

2) Calcula el diámetro crítico para el cual la viscosidad empieza a tener importancia en las oscilaciones.

3) Repetir las preguntas 1 y 2 para glicerina ((=0.77 Pa s y (=1276 kg/m3) y un diámetro D=3 cm. ¿Cambian los resultados respecto al agua?

4) En el caso de agua y D=5 mm, las gráficas muestran que la altura va disminuyendo exponencialmente respecto al tiempo. Piensa un método gráfico o numérico para obtener la constante de decaimiento k (en segundos) que aparecería en la exponencial negativa de z(t) exp(-t/k). Describe el método y calcula la constante k para el caso indicado.

Flujo por gravedad

Descripción del problema Este problema estudia el proceso transitorio que se genera al abrir la válvula de descarga de un depósito que alimenta a una línea de tuberías de cota variable.

A partir del diámetro del depósito y del diámetro de la tubería de descarga, la hoja de cálculo permite calcular la evolución temporal de la velocidad de descarga, el nivel del depósito y de la distribución de presiones en la línea de descarga.

Las cotas Z0-Z5 que se indican en la hoja de cálculo tienen su origen en una referencia inferior dibujada como un "suelo". Dichas cotas Z0-Z5 corresponden, respectivamente, con la superficie libre del depósito, el suelo del depósito y cuatro cotas más de la tubería de descarga.

Para simplificar la resolución numérica de este problema, se ha supuesto un estado cuasi-estacionario y se han despreciado las pérdidas de carga. El fluido transportado es agua.

Contesta a las siguientes preguntas Se va a trabajar con los siguientes valores: Ddepósito=50 m y Dtubería=1 m.

1) Comenta brevemente la evolución de la velocidad, el nivel del depósito y la distribución de presiones en la tubería de descarga.

2) Las gráficas muestran que la presión disminuye paulatinamente con el tiempo. Según esto, llegaría un momento que la presión disminuiría por debajo de la presión de vapor del líquido, y se produciría cavitación. Argumenta por qué esto no es posible en la situación que estamos estudiando (descarga de un depósito a la atmósfera).

3) Explica dos situaciones de un flujo en una tubería en las que sí se pueda producir cavitación.

4) Introducir z4=500 m. Se puede observar que P4=P1, que es la presión al comienzo de la tubería de descarga. ¿Es esto lógico? ¿Por qué?

5) Responde a la pregunta anterior en el caso de que se tuvieran en cuenta las pérdidas de carga.

Cálculo de chimeneas de equilibrio

Descripción del problema Las chimeneas de equilibrio (surge tanks) son conductos verticales que se comunican con la tubería de suministro de agua a las turbinas de una central hidroeléctrica con el objeto de amortiguar las variaciones de presión en los conductos. Cuando la válvula previa a las turbinas se cierra, se produce un aumento de presión en la tubería que se propaga aguas arriba. Este fenómeno se conoce con el nombre de golpe de ariete (en inglés, water hammer).

Si no se controla el golpe de ariete, el sistema puede sufrir daños serios[1]Por este motivo el aumento de presión es controlado mediante la chimenea de equilibrio. Cuando se produce este aumento de presión, parte del flujo se conduce por el camino más fácil que tiene: la chimenea. De este modo se evita que todo el aumento de la presión se produzca en una tubería cerrada, lo cual podría romper dicha tubería. Cuando se termina la onda de presión del golpe de ariete, el flujo vuelve a descender por la chimenea, lo que provoca que aumente la presión. Esto provoca de nuevo unas ondas de presión que hacen que el fluido vuelva a ascender por la chimenea, y así sucesivamente. Esto genera un movimiento oscilatorio del agua en la chimenea.

Esta simulación pretende visualizar el comportamiento de las chimeneas de compensación de presiones, analizando la evolución del agua en su interior a partir de los parámetros geométricos que definen el sistema.

Los parámetros de nuestro modelo se encuentran marcados en amarillo en la hoja de cálculo. Las variables calculadas son: altura del nivel de la chimenea (Z), velocidad de la superficie libre de la chimenea (V) y pérdidas de carga en la chimenea (hf).

Contesta a las siguientes preguntas

1) Introduce los siguientes valores: L=200 m, Dconducto=1.25 m, f=0.01, Q0=2 m3/s y Dchimenea=4 m.

Discute los resultados que se observan en la gráfica de Z, V y hf como función del tiempo. ¿Crees que las oscilaciones son amortiguadas?

2) ¿Qué relación hay entre V y hf? ¿Por qué?

3) Varía cada uno de los parámetros del problema y analiza cualitativamente su influencia en las tres variables Z, V y hf.

4) Las chimeneas de equilibrio están diseñadas para disipar parte de la energía cinética del flujo y, por tanto, amortiguar los efectos del golpe de ariete. A la vista del análisis realizado en la pregunta anterior, ¿cuál de las tres variables Z, V y hf crees que indica la energía disipada? ¿Qué parámetros de la chimenea modificarías para aumentar la eficiencia de la misma? Da algunos ejemplos numéricos obtenidos con la hoja de cálculo.

5) Calcula la velocidad de propagación de las ondas de golpe de ariete (ver datos más abajo).

6) Usando los datos de la pregunta 1, compara la energía disipada en la chimenea (ver pregunta 4) con el pico de presión de golpe de ariete dado por la fórmula de Allievi: (p=((a((u, donde a es la velocidad de propagación de las ondas en la tubería y (u es la variación de velocidad producida por el cierre de la válvula. Indica el porcentaje de energía del golpe de ariete que se disipa en la chimenea.

La velocidad de propagación de las ondas viene dada por la fórmula

donde ( es la densidad del líquido, K es el módulo de compresibilidad del líquido, E es el módulo de elasticidad del material de la tubería, D y e el diámetro y espesor de la tubería.

El valor de la constante c depende del anclaje de la tubería y del módulo de poisson (() del material de la tubería:

Tubería anclada aguas arriba y libre aguas abajo: c=5/4 – (

Tubería anclada en ambos extremos: c=1-(2

Tubería anclada en ambos extremos y con juntas de expansión: c=1-(/2

Se pueden usar los siguientes datos para agua y tubería de acero: (=1000 kg/m3, K=2.2e9 Pa, E=2e11 Pa , (=0.3, D=180 mm y e=4mm.

Suponer que la tubería está anclada en ambos extremos.

Ecuaciones

Estado estacionario en tuberías Las pérdidas de carga se pueden expresar (en unidades de altura) como:

donde K, coeficiente de pérdidas totales, se expresa como

ki siendo f el coeficiente de Darcy, L la longitud de la tubería, D su diámetro y ki las pérdidas singulares. Escribiendo la ecuación integrada de la segunda ley de Newton entre los niveles 1 y 2, y teniendo en cuenta que p1 = p2 = V1 = 0 y z1-z2 = H obtenemos:

(A)

sustituyendo el valor de hp

(B)

Si representamos la velocidad de flujo estacionario como Vest y tenemos en cuenta que para flujo estacionario dV/dt = 0, tenemos que

Sustituyendo H en la ecuación (B) y podemos despejar dt:

Integrando la ecuación anterior, y teniendo en cuenta que la constante de integración es cero, ya que V es cero cuando t es cero, obtenemos

Esta expresión también se puede poner del modo:

Esta ecuación indica que la velocidad tiende asintóticamente a Vest y que en teoría le cuesta un tiempo infinito. En la realidad, ondas elásticas y amortiguamiento provocaran que el equilibrio se alcance en un tiempo finito.

Cuando V = 0,99 Vest

ó

Oscilaciones de un tubo en U

Sin fricción En el caso sin fricción, aplicamos la ecuación de Euler-Bernoulli para flujo no permanente.

integrando y considerando p1 = p2 y V1 = V2 se obtiene:

donde L es la longitud de la columna de líquido.

Cambiando el nivel de referencia a la posición de equilibrio la ecuación anterior se transforma:

y considerando que V es únicamente función del tiempo se puede escribir

La solución general de la ecuación es

Dadas las condiciones iniciales de que cuando t = 0, z = Z y dz/dt = 0, se obtiene C1 = Z y C2 = 0, resultando la siguiente ecuación:

Esta ecuación representa el movimiento armónico simple para cualquiera de los meniscos, con periodo para una oscilación completa.

La velocidad de la columna de líquido puede hallarse derivando la coordenada z respecto del tiempo: dz/dt

Fricción laminar

Cuando un esfuerzo cortante (0 en la pared del tubo resiste el movimiento de la columna de líquido, se modifica la ecuación de Euler para el movimiento a lo largo de una línea de corriente.

Esta ecuación es válida tanto para resistencia laminar como turbulenta.

Utilizando la hipótesis de que la resistencia debida a la fricción en un flujo no permanente es la misma que la resistencia que en un flujo permanente con la misma velocidad, el esfuerzo cortante en la pared del tubo es:

Sustituyendo este valor en la ecuación anterior, integrando con respecto a s, cambiando a derivadas totales y considerando que

queda:

Las soluciones de esta ecuación son del tipo:

donde

y

Al sustituir las condiciones iniciales z = 0 y para t = 0, se obtiene:

y a partir de la derivada dz/dt, puede hallarse la velocidad v, que tiene la siguiente forma:

En la solución del problema se pueden dar dos casos diferentes:

Caso 1: Dominio de la viscosidad si

En el caso de que la viscosidad domine, ésta resulta ser tan grande que el movimiento se amortigua durante el primer período de oscilación sin que se lleguen a alcanzar valores negativos del valor de la altura del líquido con respecto a la posición de equilibrio, es decir, z.

Es posible calcular el máximo desplazamiento de la columna de líquido (z):

Caso 2: No hay dominio de la viscosidad si

En el caso de que la viscosidad no domine, el transitorio se convierte en un movimiento oscilatorio alrededor de la posición de equilibrio de amplitud decreciente, pero a lo largo de varios períodos de oscilación. El valor más grande que puede alcanzar la variable z es:

Fricción turbulenta

En la mayoría de los casos prácticos de oscilación en sistemas de tuberías, se tiene una resistencia de tipo turbulenta. En efecto, cuando se tienen tuberías grandes o túneles, el número de Reynolds suele ser grande, excepto cuando la velocidad es próxima a cero. En nuestro caso, la resistencia del fluido es proporcional al cuadrado de la velocidad promedio (se supone el factor f es constante):

Al sustituir

nos queda la siguiente ecuación:

Desarrollando esta ecuación matemáticamente y conociendo las condiciones iniciales z = z0 y para t = t0, obtenemos la siguiente expresión:

Para obtener la solución de esta ecuación se debe recurrir a métodos numéricos de resolución de ecuaciones diferenciales como el método de Runge-Kutta.

Flujo por gravedad Para este caso se puede aplicar la ecuación de Euler-Bernoulli para flujo no permanente:

Una vez definidas dos secciones, se puede integrar a lo largo de una línea de corriente:

Sin embargo en este caso como la descarga del depósito será muy prolongada en el tiempo, podemos considerar dV/dt ( 0. Por tanto, para un instante de tiempos dado, consideraremos que el sistema se encuentra en estado cuasi-estacionario.

Si tuviésemos que considerar el término de la derivada, la resolución del problema requeriría la utilización de un método numérico más complicado, como puede ser el método Runge-Kutta y su desarrollo en una hoja Excel sería más complejo. De cualquier forma, los resultados obtenidos en la simulación, corroboran la poca variación de la velocidad con el tiempo.

Con la ecuación de Euler-Bernoulli en estado estacionario y la ecuación de conservación de masa, el problema se puede resolver fácilmente tomando como referencia energética la superficie libre del depósito:

Chimeneas de equilibrio La chimenea de equilibrio se eleva por encima de la superficie de embalse para evitar sobrecargas durante la operación. En condiciones normales el nivel del agua será inferior al del embalse, en una cantidad hf igual a las pérdidas por fricción en el conducto que los une. Conforme las válvulas se cierran, el flujo se deriva hacia la chimenea de equilibrio y el nivel de agua sube por encima del embalse. A un tiempo t después de cerrar las válvulas, el nivel de la chimenea será:

Donde Z es la altura tomando como cero el nivel del agua en el embalse y hf es positiva si el flujo entra en la tubería y negativa si sale de ella.

La fuerza que se opone al movimiento del agua en el conducto de alimentación es:

Y la tasa de cambio de momento viene dada por:

Igualando las dos fuerzas, simplificando y reagrupando se tiene que:

Donde:

Si en L se incluyen como distancias equivalentes cualquier otro tipo de pérdidas en el conducto.

A la entrada de la chimenea se debe aplicar la condición de continuidad:

Donde A es la sección del conducto y Ac la sección de la chimenea y Q es el flujo hacia las turbinas.

Para resolver estas ecuaciones hay que hacer una serie de suposiciones, sobre la variación del factor de fricción y al tasa de reducción de flujo hacia las turbinas. Incluso cuando Q = 0, el factor de fricción variable todavía impide una solución analítica al problema y por eso se suele hacer una integración numérica. Sin embargo, si el factor de fricción se ignora, se puede calcular una aproximación para Zmax y el periodo de oscilación.

(Si Q = 0) Cuando t = 0, Z = 0 en el caso sin fricción se cumple que:

puesto que V es dZ/dt, integrando obtenemos:

Y el nivel máximo Zmax será:

el periodo de oscilación:

Si se incluye el factor de fricción, el nivel máximo Zmax será menor y por tanto los cálculos son conservadores. Un punto importante es que el extremo mínimo de la oscilación representado por –Zmax en el caso sin fricción pudiera ser inferior a la distancia hasta la tubería de suministro de agua y eso supondría la entrada de aire hacia la turbina. Hay que tener en cuenta este aspecto en el diseño.

Nota: Una manera de resolver las ecuaciones de manera numérica es dada un caudal inicial conocido, un factor de fricción conocido y el todos los parámetros geométricos conocidos, podemos calcular una Z y una V inicial. Así podremos resolver la ecuación:

Si fijamos un intervalos de tiempo pequeños, podemos calcular dZ y dV, que sumados a los Z y V nos permitirán ir conociendo los siguientes Z y V. Es decir, conocidos Z0, V0 y dt, podemos hallar dZ y dV de tal modo que:

Zi+1 = Zi + dZ Vi+1 = Vi + dV De este modo podemos resolver el problema incluyendo el rozamiento. Esta forma de calcularlo se basa en un método de Euler basado en diferencias hacia delante. Es necesario hacer notar que la estabilidad del método está condicionada a que los intervalos de tiempo sean menores a un intervalo de tiempo máximo, es decir:

(t < (tmáx Si el intervalo de tiempo es mayor a este intervalo máximo, la solución diverge. También es necesario hacer notar que este método numérico es muy básico, ya que es de orden 1, es decir, que el error en el resultado es del orden de (t1. Podría haberse optado por implementar un método de orden mayor, como por ejemplo el método de Runge-Kutta clásico, de orden 4, pero eso se deja como ejercicio voluntario para el alumno. De esta manera se ha operado en la hoja Microsoft( Excel para obtener valores numéricos.

1- Ver, a modo de ejemplo, el accidente de la central nuclear de Three Mile Island (Middletown, Pennsylvania, EE.UU.) en 1979, o el accidente en el Brookhaven National Laboratory (1986).

Autor:

PabloTurmero