Flujo de un fluido viscoso a la entrada/salida de una cavidad anular de longitud finita

Revista Internacional de Métodos Numéricos para Cálculo y Diseño en Ingeniería. Vol. 12,1, 17-29(1996) FLUJO DE UN FLUIDO VISCOSO A LA ENTRADAISALIDA DE UNA CAVIDAD ANULAR DE LONGITUD FINITA CARLOS A. VIONNET Department of Aerospace and Mechanical Engineering The University of Arizona Tucson, AZ 85721, U.S.A. RESUMEN En este trabajo se presenta un procedimiento sencillo para imponer condiciones de borde en contornos artificiales a fin de resolver las ecuaciones de Navier-Stokes en presencia de efectos centrífugos. El análisis se concentra en el flujo de un fluido viscoso atrapado en una cavidad anular de longitud finita. El algoritmo se basa en desacoplar una componente estática del campo de presiones, la cual puede ser tratada con la aproximación penalizada. Los resultados numéricos exhiben un excelente acuerdo con resultados asintóticos obtenidos mediante una expansión basada en el método de las perturbaciones. SUMMARY A simple algorithm for imposing open boundary conditions in the numerical solution of the Navier-Stokes equations in presence of centrifugal forces is presented. The analysis concentrates on the flow of a viscous fluid trapped inside an annulus of finite length. The algorithm is based on decoupling a static component of the pressure field, which can be naturally treated with the penalty function approximation. The numerical results show an excellent agreement with asymptotic results obtained with a regular perturbation expansion. Determinar las características del flujo de un fluido viscoso incompresible contenido entre dos cilindros concéntricos y rotando a diferentes velocidades es un problema que ha concitado el interés de los investigadores desde los tiempos de CouetteZ5. Problemas de esta índole pueden encontrarse en situaciones tan variadas como las que van desde el empleo de viscosímetros hasta la lubricación de cojinetes de ejes. Cuando el cilindro externo permanece en reposo, mientras el interno rota a una velocidad angular moderada, se establece un flujo ordenado o "laminar" comúnmente denominado flujo de Couette. A medida que la velocidad de rotación aumenta, este flujo deviene inestable, Recibido: Septiembre 1994 OUniversitat Politecnica de Catalunya (España) ISSN 0213-1315 18 C.A.VIONNET siendo los vórtices toroidales de Taylor el primero de una serie de estados transicionales hasta llegarse a la turbulencia ~lenal'. Para el caso de cilindros de longitud axial infinita, la solución de Couette es bien conocida1; sin embargo, el análisis no sólo pierde vigencia para velocidades angulares apreciables (el flujo es inestable ante perturbaciones de pequeña amplitud), sino que además existen dudas acerca del posible efecto de los bordes extremos sobre el comportamiento global del flujo, sean éstos abiertos o cerrados. Es claro que una solución aproximada a este problema, independientemente del régimen de flujo y del método numérico a emplear, depende en gran medida de las condiciones de borde impuestas en los extremos de la cavidad anular. La complicación es aún mayor si se tienen en cuenta oscilaciones axiales relativas entre los cilindros y propiedades que varían con la temperatura. Para el régimen estable, el caso límite de curvatura cero fue analizado por Shepherd y Di Primaz1 con ayuda de expansiones asintóticas. Barret et al.' incluyeron un factor para corregir los efectos de longitud finita en la teoría de lubricación de cilindros levemente excéntricos. En forma similar, Goenka y OHg emplearon separación de variables en un método numérico para resolver la misma clase de problemas. Para el caso más teórico de inestabilidades de Taylor-Couette, CliffeG, y más recientemente Hill15, entre otros, resolvieron el problema numéricamente para cilindros concéntricos pero sellados en sus extremos, equivalente numérico del aparato originalmente diseñado por G.I. TaylorZ2. En el presente trabajo el análisis se restringe al régimen estable de flujos tipo Couette en cavidades anulares abiertas en sus extremos, a su posterior solución numérica utilizando la formulación penalizada de la ecuación de conservación de la masa en el contexto del método de los elementos finitos (MEF) y la formulación de condiciones de borde consistente con tal propósito. FORMULACI~N DEL PROBLEMA Sea una superficie cilíndrica tal como la ilustrada en la Figura 1, en reposo con respecto a un sistema inercial, y que recubre un eje de longitud infinita que rota a velocidad angular R constante. Suponiendo que el espacio que separa ambos cuerpos está totalmente ocupado por un fluido viscoso incompresible y cuyo flujo cubre una extensión axial mucho mayor que b, la mayor dificultad aparece al truncar el dominio a fin de hacer el problema tratable desde el punto de vista computacional. Es sabido que la mejor elección es ubicar los contornos artificiales donde el flujo está totalmente desarrollado o, en su defecto, donde la acción se ha atenuado a niveles tales que es aceptable ignorar la influencia del dominio exterior. Esto no siempre es posible a menos que se cuente con información a priori sobre el movimiento del fluido14. Con referencia a la Figura 1, la geometría en la vecindad de las zonas z = O y z = b depende del problema particular considerado. Por lo tanto y dada su simplicidad, el análisis se restringe a la región O 5 z 5 b, R 5 r 5 R + d asumiendo que la superficie cilíndrica externa presenta una transición suave más allá de z < O y z > b. El objetivo es entonces establecer condiciones de borde capaces de permitir el pasaje de fluido a través de los contornos artificiales sin perturbar la solución en la región de interés. FLUJO VISCOSO EN UNA CAVIDAD ANULAR 19 Figura 1. Descripción de la cavidad anular Partiendo de las ecuaciones de Navier-Stokes escritas en coordenadas cilíndricas, en ausencia de superficies libres y asumiendo propiedades variables con la temperatura, es posible expresar las ecuaciones de gobierno como para las cuales se utilizó la siguiente normalización Además, el sistema de ecuaciones (1)-(5) está sujeto a ciertas condiciones de borde e iniciales. En (6), p* es la presión absoluta, (u(,), u(,), las componentes de la velocidad del flujo en (z, r, 8) respectivamente y T la temperatura del fluido. Las propiedades características del fluido son: densidad p, viscosidad absoluta p = p(T), calor específico a presión constante Cp y conductividad térmica k = k(T). Los valores de referencia son p8, To, /LO = p(To) y cierto incremento de temperatura lATl con respecto a To, generalmente impuesto por las condiciones de borde. Asimismo B (u, u), 2 E (x, y) , V (ax, ay) y uX = aXu (derivada parcial de u respecto a x). Los parámetros típicos son el número de Taylor Tu, el número de Reynolds Re, el número de Prandtl P,, el número de Eckert E, y la viscosidad p relativa a la viscosidad de referencia PO = pvo (donde vo es la viscosidad cinemática) definidos como R2 ~d~ RRd pCpv0 E - a2 R~ Tu = - ,Re=- , P,=- P , p*=- u: vo k ' - C~~AT~ Po (7) El sistema de ecuaciones (1)-(5) se obtiene a partir del sistema cilíndrico completo de Navier-Stokes en el límite formal1' 6 = $ -, O, cuyo rasgo dominante es su naturaleza axisimétrica donde los efectos de curvatura son ignorados excepto en el término de la aceleración centrípetaz3. Anticipando la penalización de la ecuación de continuidad, la presión es ahora desdoblada en la forma p = P + p,, donde p, es la componente estática, mientras que la componente remanente P es atribuida enteramente a la acción dinámica del fluido. La aplicación del método de los residuos ponderados conduce mediante la primera identidad de Green en un dominio w simplemente conexo y acotado por aw, a la forma débil de las ecuaciones de movimiento, a saber las que deben satisfacerse para toda función de peso l$' = (Wl, W2). Aquí y siguiendo la convención de tensores cartesianos, índices repetidos se suman tal que (u, u) (ui, u2) representa las componentes de la velocidad en las direcciones (x, y) z (xl, x2), g = i[VÜ + (vÜ)'] e1 tensor de la tasa de deformación, S,j el tensor delta de Kronecker y 6 = (nl, n2) el vector normal exterior a aw. La forma débil de la ecuación de movimiento en dirección azimutal (4) y de la ecuación de energía (5), que por brevedad no se incluyen, se obtienen en forma similar. Por último, la ecuación de conservación de la masa (1) se incorpora mediante la relación donde X es el parámetro de penali~ación~~~". De esta forma la presión es eliminada del dominio w, excepto en el contorno aw. DISCRETIZACI~N La partición del dominio w en M elementos finitos y N nodos conduce a la siguiente aproximación de las funciones incógnitas donde el supraíndice h denota la aproximación discreta de la variable dependiente f, ya sea Ü, w ó 4 , sobre la partición wh del dominio w. Las $j(Z) se refieren indistintamente FLUJO VISCOSO EN UNA CAVIDAD ANULAR a ciertas funciones de interpolación globales asociadas a cada nodo j, en cuyo caso la integración es sobre todo el dominio w, o a las conocidas funciones de interpolación locales obtenidas con el MEF (funciones sombrero en 1D o bilineales en 2D, por ejemplo) reduciéndose ahora la integración a un elemento individual de área w,. La aplicación directa del método de Galerkin (en rigor de Petrov-Galerkin) sobre la mencionada partición de elementos conduce al siguiente sistema de ecuaciones diferenciales R O ~($2) B1l(X1p*) BG(A,p*) u - & [O R] [:] + [c(%u)] + [B21(X1~*) B22(X1~*)] [E] - [oE] (11> Aquí 3 = (u1, u2,. . . , UN), u = (u1, v2,. . . , UN), = (wl, w2,. . . , wN) son los vectores incógnitas de las velocidades nodales, 4 = ($1, 42,. . . , 4~), el vector incógnita de temperaturas nodales , R la matriz de mas<^ x N, B,, la matriz de penalización más difusión viscosa N x N, para m, n = 1,2, ~(ü, g) un vector no lineal N x 1 función de la velocidad de convección U, A(a) la matriz de difusión N x N para a = p* ó a = 1/P, según corresponda, y og los vectores N x 1 de las condiciones de borde y las fuerzas de volumen y/o contacto.-LOS términos individuales se calculan a partir de donde los términos penalizados se evalúan utilizando integración reducida4. Los términos convectivos son estabilizados mediante el conocido artificio de adicionar disipación numérica en la dirección paralela al flujo3>'' (SUPG-Streamline Upwind Petrov-Galerkin) a fin de minimizar los efectos de difusión numérica en la dirección transversal1'. De esta forma, cada componente de los vectores c se obtiene por medio de la siguiente expresión U - Vgdw - donde - g designa genéricamente los vectores o 2. Cuando - g se reemplaza por o - q5 se obtiene la matriz convectiva C(i2). Arriba, IIi2II /mes el módulo de la velocidad del flujo y k. el coeficiente de disipación agregada que se calcula con ayuda de siendo Jk la llamada función de "upwind" que depende del número de celda de Peclet definido en la forma 1 1 tk(Pk) = - - - , Pk=- hk Re (no se suma sobre le), k = 1,2 (20) th(Pk) Pk 2 Los números de celda de Peclet se evalúan en los puntos de integración gaussianos y las dimensiones del elemento hl y h2 se calculan de acuerdo a lo sugerido por Yu & HeinrichZ6. Las condiciones de borde esenciales CBE o de Dirichlet se engloban en los vectores dg siguiendo el procedimiento conocido del MEFZ7. La función de disipación viscosa eñ la ecuación de energía se calcula asumiendo números de Reynolds de moderados a altos, para lo cual la contribución debida a (u, u) se considera despreciable en comparación con la disipación inducida por w. En consecuencia Finalmente, la componente crítica para la clase de problema que se quiere resolver está dada por que se obtiene a partir de la ecuación de equilibrio dps = T~wwS~~, asociada al estado básico (u,, u,, w,; 4,) (0,0,1 -y; O) del sistema (1)-(5) y al empleo de la aproximación de los productos5. Esta solución básica constituye el punto de partida para el análisis de la estabilidad del flujo ante perturbaciones de pequeña amplituds. Sin embargo, su presencia en el presente contexto es al solo fin de relajar las dificultades inherentes en la implementación de las condiciones de borde para flujos dominados por fuerzas de flotación en dominios artificialmente truncados. FLUJO VISCOSO EN UNA CAVIDAD ANULAR 23 El algoritmo de avance temporal del sistema de ecuaciones (1)-(5) se basa en el conocido método theta con matrices de masa condensadas (lump). Redisponiendo el vector ((g,~)~ u = (ul, vi, ~2~~2,. . . ,UN, vN) a fin de minimizar el semiancho de banda de la matriz de coeficientes y, asumiendo que se parte de condiciones iniciales dadas, para cada nivel de tiempo tn = nAt se evalúa y resuelve el sistema utilizando un método directo de solución basado en la eliminación de Gauss con pivote0 parcial para matrices banda asimétricas7. Aquí J = dfldu es el jacobiano del sistema de ecuaciones no lineales f (un+') = [R + ~A~B(~*)~]U"+~ + 6'Atc(un+') - {[R - (1 - 6')AtB(p*)n]un- - (1 - 6')Atc(un) + Atbt} = O (24) donde R y B son las matrices aumentadas de la ecuación (ll), At el paso temporal y theta un coeficiente de peso tal que 6' = 1 da el algoritmo totalmente implícito de Euler. El proceso descrito en (23) continúa hasta que en la v-ésima iteración se verifica 2 N "AUYII < 0.01 para 11ull = lujl IIU~+'II j=i Una vez actualizada u, se resuelve la ecuación para 3 en la forma A continuación con se evalúa b4 conforme a (21) y de manera similar a (26) - se obtiene q5n+1. Luego se actualizan p* = p*(q5) y PT = PT(q5) según las características propias delfluido considerado y se vuelve a (24) hasta que la diferencia relativa entre dos soluciones sucesivas sea menor que una tolerancia prefijada. Este esquema es válido para soluciones estacionarias. Para obtener una precisión de segundo orden en el tiempo, es decir Crank-Nicolson-Galerkin para 6' = 0.5, es necesario redistribuir el peso entre los vectores bu y b, en (24) y (26) respectivamente y ajustar los niveles temporales según una expansiónde Taylor4 en un entorno de t + At/2. A menos que p* varíe significativamente con la temperatura entre dos niveles de tiempo sucesivos, los cambios a realizar en el código numérico son mínimos. La Figura 2 muestra las dimensiones del problema y las condiciones de borde utilizadas, donde se supone que la superficie cilíndrica exterior experimenta un desplazamiento axial con velocidad unitaria. El dominio es discretizado con 800 elementos cuadrangulares (bilineales) y 858 nodos con paso uniforme Ax = .125 y Ay = .04. En los contornos artificiales dw ubicados en x = O y x = 4 se imponen las condiciones de borde siguientes: (i) se asume que p, r O en todo el dominio wU&, lo cual es equivalente a trabajar con P = p (presión total), y se ignora la integral de línea en (8). (ii) se asume p, # O en todo el dominio wUdw, p = p, y tensiones viscosas despreciables en 8,. En suma, la integral de línea (8) es nuevamente nula. Nótese que para implementar (ii) es necesario definir w, en (22). Tal como se mencionó anteriormente, esto se obtiene haciendo u, = v, = 4, = O en las ecuaciones de gobierno. Luego, el sistema (1)-(5) se reduce a wl1 = O sujeto a w,(O) = = 1, ws(1) = - "Ir=,+, = O (27) donde ' d/dy. Integrando ambas ecuaciones se obtiene Figura 2. Esquema del dominio computacional y condiciones de borde Esta solución representa el incremento necesario de presión en dirección radial para balancear los efectos de la aceleración centrífuga. El flujo principal así establecido se denomina flujo plano de Couette (invertido) para el caso límite de curvatura cero. Al emplearse la condición (ii), esta componente del flujo total es automáticamente eliminada del problema. Esto se ilustra en el ejemplo numérico donde se asume una FLUJO VISCOSO EN UNA CAVIDAD ANULAR 25 dependencia p* = exp(-€4). En el apéndice se incluye una expansión en potencias de E que permite una comparación cuantitativa con los resultados numéricos. Puede verse que el estado básico corresponde a la solución de orden cero o al caso de propiedades constantes, independientemente de la existencia de un flujo secundario u # O en la medida en que este no afecte el equilibrio hidrostático (28). En la Figura 3 se aprecian los contornos de igual presión total y función de corriente así como también el campo de velocidades obtenidos con la condición (i). A pesar de que la solución exhibe el tipo de flujo esperado en la zona media de la cavidad, hay un vórtice claramente visible a la salida y las líneas de corriente presentan una curvatura pronunciada en la vecindad de x = O. Por el contrario, la Figura 4 muestra la transición suave del flujo, tanto a la entrada como a la salida, cuando la condición (ii) es utilizada. Por último, la Figura 5 muestra el excelente acuerdo entre resultados numéricos y aproximados obtenidos con esta última condición. Los valores de los parámetros utilizados para este ejemplo son: T, = 20, P, = 0.7, E, = 10 y E = 0.2. 1, .... 1 .... c.... I ... 8. .. 0.0 1 .O 2.0 3.0 4.0 Figura 3. Flujo en la cavidad anular cuando la condición de borde (i) es utilizada. Presión total: p; líneas de corriente: $; campo de velocidades: (u, u)