Los elementos finitos de alto orden (hp-FEM) como método de cálculo en problemas de estabilidad fluidodinámica

XX Congreso de Ecuaciones Diferenciales y Aplicaciones X Congreso de Matem´ atica Aplicada Sevilla, 24-28 septiembre 2007 (pp. 1–8) Los elementos finitos de alto orden (hp-FEM) como m´etodo de c´alculo en problemas de estabilidad fluidodin´amica. Maite Pe˜ na Alcar´ az1, Leo M. Gonz´ alez2, Vassilis Theofilis3 1ETS de Ingenier´ıa ICAI, Universidad Pontificia Comillas de Madrid, Alberto Aguilera 23, 28015 Madrid. E-mail: [email protected]. 2Dpto. de Ense˜nanzas B´asicas de la Ingenier´ıa Naval, ETSIN-UPM, Arco de la Victoria s/n, 28040 Madrid. E-mail: [email protected]. 3ETS Ingenieros Aerona´uticos, Univ. Polit´ecnica de Madrid, Plaza Cardenal Cisneros 3, 28040 Madrid. E-mail: [email protected]. Palabras clave: Bi-Global instability, high order methods, hp-FEM Resumen En los problemas de estabilidad BiGlobal en un contexto fluidodin´amico se han encontrado un cierto tipo de condicionantes que todo buen m´etodo num´erico debe ser capaz de satisfacer. Por un lado y como condici´on m´as indispensable, se encuentra la libertad geom´etrica, es decir, el m´etodo debe ser capaz de discretizar cualquier dominio espacial de forma que no quede sujeto a formas simples o f´acilmente transformables a ´estas. En segundo lugar se desea que las condiciones de contorno necesarias en los cl´asicos problemas de Mec´anica de Fluidos sean f´acilmente implementables. Por ´ultimo que el orden del m´etodo elegido en la resoluci´on del problema sea f´acilmente variable en funci´on de la precisi´on exigida. En este contexto se conoce que la precisi´on necesaria crece a trav´es de un par´ametro que no es otro que el n´umero de Reynolds del problema, y por tanto es muy ventajoso que no se sea el tama˜no de la malla el ´unico grado de libertad para mejorar la precisi´on del c´alculo, ver [2] y [3]. Con todos estos condicionantes se ha construido un c´odigo basado el elementos finitos de alto orden (hp-FEM) que ha sido aplicado al c´alculo de los valores propios y los modos propios de algunos problemas cl´asicos de la Mec´anica de Fluidos como son el movimiento de fluidos en un conductos. En los casos aqu´ı presentados, hemos aplicado el m´etodo de Arnoldi para el c´alculo del espectro y los modos propios en los problemas de autovalores generalizados que aparecen. Como resultados del trabajo se muestran el c´alculo del flujo base de una tuber´ıa rectangular y el posterior an´alisis de estabilidad BiGlobal de dicho flujo base. Los resultados han sido coherentes y en buena consonancia con los obtenidos mediante otros m´etodos en la literatura cl´asica, ver [1]. 1 Maite Pe˜na Alcar´az, Leo M. Gonz´alez, Vassilis Theofilis 1. Introducci´on La descripci´on cl´asica de la teor´ıa de estabilidad BiGlobal se basa en una perturbaci´on lineal de las ecuaciones de Navier-Stokes alrededor de un flujo base (U, V , W), el cual es una soluci´on particular de las mismas. El c´alculo de esta soluci´on particular (U, V , W) se ha obtenido resolviendo simplificaciones bidimensionales o unidimensionales de las ecuaciones de Navier-Stokes mediante una discretizaci´on de alto orden hp-FEM. Una vez se ha calculado el flujo base se ha linealizado el problema entorno a dicha soluci´on mediante el ansatz: ui(x, y, z, t) = Ui(x, y) + ²bui(x, y)eiωteiβz (1) p(x, y, z, t) = P(x, y) + ²bp(x, y)eiωteiβz (2) Considerando ²¿1, tenemos que (bui,bp) representan las amplitudes de la perturbaci´on, el n´umero complejo ωcontiene como parte imaginaria la tasa de crecimiento y como parte real la frequencia de dicha perturbaci´on, por ´ultimo β=2π LzyLzson el n´umero de onda y la longitud de onda de la perturbaci´on, respectivamente. Al introducir este ansatz en las ecuaciones de Navier-Stokes y considerar el flujo base como una soluci´on particular se llega a un problema de autovalores generalizado del tipo: A    ˆu1 ˆu2 ˆu3 ˆp    =−ωB     ˆu1 ˆu2 ˆu3 ˆp    .(3) Donde el operador complejo y no sim´etrico Ase define como: A=     a11 ∂¯u1 ∂y 0∂ ∂x ∂¯u2 ∂x a22 0∂ ∂y ∂¯u3 ∂x ∂¯u3 ∂y a33 iβ ∂ ∂x ∂ ∂y iβ 0     .(4) aii =(¯uj ∂ ∂xj +∂¯ui ∂xi −1 Re Ã∂2 ∂x2 j −β2!+iβ ¯u3), j = 1,2.(5) y no se considerar´a el convenio de indices repetidos de Einstein para la ecuaci´on (5). Por otro lado el operador real sim´etrico Bse define como: B=    M0 0 0 0M0 0 0 0 M0 0 0 0 0    ,(6) donde Mrepresenta la cl´asica matriz de masas de elementos finitos. 2 Normas para la edici´on (CEDYA 2007) En la ecuaci´on (3) el valor propio viene dado por ωy el autovector (bu, bv, bw, bp) contiene las amplitudes de las perturbaciones. La resoluci´on directa del problema de autovalores generalizado mediante un m´etodo tipo QZ implicar´ıa un coste computacional excesivo para matrices AyBde gran dimensi´on, es por ello que la metodolog´ıa escogida para la resoluci´on de este tipo de problemas ha sido un m´etodo iterativo basado en la creaci´on de un subespacio de Krylov conocido como m´etodo de Arnoldi. Para ello en cada iteraci´on se ha resuelto un problema lineal que implicaba a la matriz Acomo matriz de coeficientes y que por tanto obligaba a una descomposici´on LU de la misma al principio del c´alculo. 2. Discretizaci´on espacial del problema. Elementos finitos de alto orden: hp-FEM La resoluci´on num´erica del problema de autovalores generalizado descrito anteriormente, requiere una discretizaci´on espacial que se har´a mediante el uso de una base de funciones, Φq(x). Algunas versiones de este tipo de c´alculos y discretizaciones cuando las funciones base son lagrangianas, nodales y de bajo orden se puede ver en [1]. De cualquier modo la experiencia con funciones de bajo orden muestra que para resolver estructuras a n´umeros de Reynolds moderados se requiere de mallas muy finas con el consecuente consumo de memoria y tiempo de computaci´on. En el presente trabajo se propone una discretizaci´on modal de alto orden basada en polinomios de Legendre Lp(x). El hecho de que sea modal implica que el valor de los coeficientes de la expansi´on no tiene interpretaci´on f´ısica, por otro lado sabemos que las expansiones modales tienen estructura jer´arquica, lo que implica que un conjunto de funciones hasta orden P−1 est´a incluida en una expansi´on de orden P. Una de las propiedades fundamentales de este tipo de expansiones es la ortogonalidad, que junto con el hecho de poseer estructura jer´arquica da lugar a matrices bien condicionadas [2]. Para problemas como los aqu´ı presentados que implican un segundo orden de diferenciaci´on basta con garantizar que la soluci´on se halla en H1. T´ıpicamente, el m´etodo de los elementos finitos resuelve este problema mediante la imposici´on de una condici´on de continuidad C0entre elementos vecinos, ello implica que las expansiones globales son continuas en todo el dominio mientras que la continuidad en la derivada se consigue con la convergencia, ver [2]. En un intervalo est´andar Ω = {ξ| − 1≤ξ≤1}una expansi´on modal de orden pse denota ψp(x) y se define como ψp(ξ) =                      1−ξ 2p= 0 ·1−ξ 2¸·1 + ξ 2¸Lp−1(ξ) 0 < p < P 1 + ξ 2p=P (7) En caso de tener elementos que no coinciden con el elemento est´andar [−1,1], ser´a posi3 Maite Pe˜na Alcar´az, Leo M. Gonz´alez, Vassilis Theofilis ble transformarlos a dicho intervalo mediante una transformaci´on lineal. Por lo tanto ´unicamente definiremos nuestras funciones base para el intervalo est´andar. Del mismo modo tambi´en es necesario se˜nalar que dado que los elementos finitos requieren una formulaci´on tipo Galerkin es necesario elegir un m´etodo de integraci´on num´erica que en ´este caso ser´a el de Gauss. Para una situaci´on multidimensional se puede generalizar esta metodolog´ıa, y m´as concretamente en un caso bidimensional la discretizaci´on se puede realizar tanto mediante elementos de tipo triangular o de tipo cuadril´atero. De este modo las funciones base se pueden expandir a m´as dimensiones mediante un producto tensorial de funciones base unidimensionales ψpq(ξ1, ξ2) = ψp(ξ1)ψq(ξ2). En el caso de un cuadril´atero est´andar, las coordenadas cartesianas (ξ1, ξ2) quedan definidas mediante contornos constantes Q2={(ξ1, ξ2)| − 1≤ξ1, ξ2≤1}. Por el contrario, para un elemento triangular est´andar, donde hay una dependencia entre las coordenadas cartesianas (ξ1, ξ2) a la hora de definir el dominio, es decir T2= {(ξ1, ξ2)| − 1≤ξ1, ξ2, ξ1+ξ2≤0}. Un conjunto de funciones base definidas en elemntos triangulares ha sido desarrollada por Karniadakis & Sherwin [2], presentando la ventaja de poseer un sistema de coordenadas en el que no hay dependencia entre las coordenadas para la definici´on del dominio. La gran ventaja de ello es la posibilidad de realizar expansiones multidimensionales mediante productos tensoriales de las elaboradas en una dimension. La forma de trabajar con elementos triangulares que aqu´ı se propone, ver [2], se basa en una transformaci´on de coordenadas, mediante el cambio: η1= 21 + ξ1 1−ξ2 −1, η2=ξ2.(8) Por lo tanto, un paso necesario para realizar un c´alculo mediante elementos finitos hpFEM require de diversas transformaciones de los elementos geom´etricos a los elementos est´andar [−1,1] en el caso 1-D, Q2para el caso 2-D y elementos tipo cuadril´ateros y T2 para el caso 2-D y elementos tipo triangular. 2.1. Estudio de la Forma de las Matrices de Masas y Rigidez A consecuencia de la ortogonalidad de los polinomios de Legendre Lp−1(ξ) las matrices de masas y rigidez son, respectivamente, triy pentadiagonales. Este tipo de discretizaci´on mediante elementos finitos tipo hp−(FEM) presenta en este sentido una gran ventaja frente a discretizaciones tipo espectral Chebyshev, dado que en este ´ultimo caso la matriz de rigidez es densa. Se consideran las funciones definidas por la ecuaci´on (7), y se definen las matrices sim´etricas M={Mij}i,j=0,...,P yR={Rij}i,j=0,...,P donde Mij =Z1 −1 ψiψjdξ Rij =Z1 −1 ∂ψi ∂ξ ∂ψj ∂ξ dξ Se consideran conocidas las siguientes propiedades de los polinomios de Legendre: 4 Normas para la edici´on (CEDYA 2007) Z1 −1 LpLqdξ =2 2p+ 1δp,q,(9) (k+ 1)Lk+1(ξ) = (2k+ 1)ξLk(ξ)−kLk−1(ξ),(10) (2k+ 1)Lk(ξ) = L0 k+1(ξ)−L0 k−1(ξ),(11) gr(Lk) = k(12) siendo gr(p) el grado del polinomio p(ξ) y se denotar´a ∂f ∂ξ con f0of0(ξ). Por otro lado se sabe que si Lp(ξ) es el polinomio de Legendre de orden p, entonces ser´a un polinomio par si pes par e impar si pes impar. Del mismo modo si Lp(ξ) es un polinomio par, entonces L0 p(ξ) es impar o viceversa. Teorema. Mij = 0 si i+jes impar, con i, j 6= 0, P Demostraci´on. Si i+jes impar, iser´a impar o jser´a impar. Se puede suponer por simetr´ıa en el razonamiento que ies impar y jpar. Se tiene que Mij =M2i0−1,2j0=Z1 −1µ1−ξ2 4¶2 L2i0L2j0−1dξ = 0 por ser producto de dos funciones pares y una impar, y por tanto ser una funci´on impar integrada entre −1 y 1, por lo que la integral se anula. Teorema. Rij = 0 si i+jes impar, con i, j 6= 0, P Demostraci´on.Si i, j 6= 0, P ,ψk(ξ) = 1−ξ2 2Lk−1(ξ) con k=i, j por lo que es inmediato que si kes impar, entonces ψkes impar y viceversa. Por tanto ψ0 kser´a par si y s´olo si klo es y ser´a impar en caso contrario. Por tanto un razonamiento an´alogo al anterior permite concluir que si i+jes impar, tambi´en lo ser´a el producto ψiψjy por tanto la integral se anular´a. Puesto que ψ0 0=1 2yψ0 P=−1 2se obtiene de forma inmediata que Ri0=RiP = 0 si i es impar. Para el c´alculo de las matrices de masas y rigidez a partir del grado de los polinomios, tenemos la siguiente proposici´on: Proposici´on. Para cualquier polinomio de Legendre Lp(ξ) se tiene que Z1 −1 Lp(ξ)p(ξ)dξ = 0 si p(ξ) es un polinomio de grado g < p. Demostraci´on. Si p(ξ) es un polinomio de grado menor que pentonces p(ξ) se puede expresar como combinaci´on lineal de polinomios de Legendre de grado estrictamente menor que p. Por tanto p(ξ) = PαiLi(ξ) y ello implica que Z1 −1 Lpp(ξ)dξ =Z1 −1XαiLiLpdξ =XαiZ1 −1 LiLpdξ =Xαiδi,p 2 2p+ 1 = 0 5 Maite Pe˜na Alcar´az, Leo M. Gonz´alez, Vassilis Theofilis Teorema. R1 −1ψiψjdξ = 0 si |i−j|>4, i, j 6= 0, P Demostraci´on. Se puede suponer sin p´erdida de la generalidad que i≥j. Entonces Z1 −1 ψiψjdξ =Z1 −1µ1−ξ2 2¶2 Lj−1Li−1dξ =Z1 −1 p(ξ)Li−1dξ donde p(ξ) es un polinomio de grado j−1 + 4 = j−3. Por tanto si i−j > 4 la prueba se acaba utilizando el resultado de la proposici´on anterior. En el caso particular en que iojsean igual a 0 o Pse obtiene el siguiente resultado cuya demostraci´on de este teorema es an´aloga a la del teorema anterior: Z1 −1 ψ0ψjdξ =Z1 −1 ψPψjdξ = 0 si j > 4, j 6= 0, P Ahora se ha terminado el an´alisis de todos los t´erminos de la matriz de masas exceptuando la diagonal principal, dos diagonales secundarias a distancia 2 y 4 y algunos t´erminos dispersos como i= 0, P. Para terminar un simple c´alculo permite asegurar que en general ning´un otro t´ermino de la matriz ha de ser nulo. Para la matriz de rigidez el razonamiento es muy parecido, basta considerar el siguiente: Teorema. R1 −1 ∂ψi ∂ξ ∂ψj ∂ξ dξ =R1 −1 ∂2ψi ∂ξ2ψjdξ ∀i, j 6= 0, P Demostraci´on. Basta integrar por partes sabiendo que la funci´on ψi(−1) = ψi(1) = 0. Una vez que se tiene este resultado, la forma de la matriz de rigidez se obtiene de forma inmediata. Teorema. R1 −1 ∂ψi ∂ξ ∂ψj ∂ξ dξ = 0 si i, j 6= 0, P y se cumple que i > j + 2 o viceversa. Demostraci´on. Usando que la integral del corolario es igual R1 −1 ∂2ψj ∂ξ2ψidξ y sabiendo que ∂2ψj ∂ξ2ψi= Li−1p(ξ) donde el grado de p(ξ) es j+ 1, con lo que si i−1> j + 1 la integral ser´a nula como ya se ha probado. Por ´ultimo, falta analizar la forma de los t´erminos de la matriz de rigidez cuando alguno de los ´ındices es igual a 0, P . En estos casos se cumple el siguiente resultado: Teorema. R1 −1 ∂ψi ∂ξ ∂ψj ∂ξ dξ = 0 ∀j6= 0, P ,i= 0, P Demostraci´on. Basta realizar un simple c´alculo para comprobar que |∂ψi ∂ξ |=1 2, por lo que la integral, queda como: R1 −1 ∂ψj ∂ξ dξ =k·(ψj(1) −ψj(−1)) = 0. Por tanto, los ´unicos elementos no nulos de la matriz de rigidez Ri,j son aquellos donde |(i−j)| ≤ 2, con i, j ambos pares o impares si i, j 6= 0, P yR0,0, R0,P , RP,0, RP,P . Como resumen, en las dos siguientes matrices se ha dibujado la forma de los elementos no nulos de MyR: 6 Normas para la edici´on (CEDYA 2007) M=          • • • • • • • • • • • • • • • • • • • • • • • • • • • • ••••• ••••• • • • • • • • • ••• ••• • • • • • •          R=      • • • • • • • • • ••• ••• • • • • • • • • •       3. Aspectos num´ericos Para asegurar la estabilidad de una discretizaci´on mediante el m´etodo de los elementos finitos de las ecuaciones de Navier-Stokes o problemas derivados de ellas, se sabe que la condici´on de compatibilidad inf −sup debe ser satisfecha entre aquellos espacios donde se discretiza la velocidad y la presi´on. En el presente trabajo el n´umero de modos empleados para la discretizaci´on de la presi´on p−1 ha sido rebajado en una unidad respecto al de velocidad p. La dimension mdel subespacio de Krylov para todos los calculos que se presentan a continuaci´on ha sido de 100. Por otro lado, el n´umero de elementos de matriz complejos que se reservan para el c´alculo del operador Ason 4p4N2 e. De este modo las estimaciones de memoria RAM en Mb necesaria para llevar a cabo estos c´alculos ha sido 64p4N2 e 10242= 64p4µNe 1024¶2 =µ8p2Ne 1024 ¶2 (13) A consecuencia de todo esto el problema m´as grande resoluble en un PC con 4 Gb de memoria RAM viene determinado por la condici´on Nep2≤8000. 4. Aplicaci´on. En este caso se considerar´a un flujo base unidimensional (0,0, W) creado en una tuber´ıa de secci´on triangular o cuadrada que se extiende a lo largo del eje Zcon gradiente de presi´on favorable y constante. El caso triangular el dominio estar´a formado por los puntos interiores al tri´angulo definido por los puntos Γb={(1,−1),(−1,1),(−1,−1)}. Para obtener el flujo base se resuelve num´ericamente la ecuaci´on de Poisson, ∇2¯ W(x, y) = −2,¯ W(x, y)Γb= 0,(14) El an´alisis de dicho flujo base contiene, como se puede ver en las ecuaciones (4)-(5), una pareja de par´ametros libres: Re yβ. En nuestro caso se han realizado an´alisis para un amplio rango de valores de β∈(0,25) y Re ≤1000. Por otra parte, en el caso de una tuber´ıa rectangular tambi´en se podr´ıa hablar de un tercer par´ametro geom´etrico como es el aspecto geom´etrico A, que no es otra cosa que el cociente entre los lados del conducto. La convergencia en el c´alculo de los autovalores se consigue tanto mediante un aumento del n´umero de elementos h, como del orden polin´omico p, ver tabla 1. En la figura 1 se puede puede ver la forma de las componentes de los modos asociados al autovalor menos estable en el caso de tener contorno triangular. 7 Maite Pe˜na Alcar´az, Leo M. Gonz´alez, Vassilis Theofilis h p ωrωih p ωrωi 1 8 0.5991425642 -0.1490564284 2 4 0.5982952805 -0.1374149414 1 12 0.5941967173 -0.1404610035 2 8 0.5941884517 -0.1405027982 1 16 0.5941885536 -0.1405012852 4 4 0.5941712455 -0.1404787289 1 20 0.5941885590 -0.1405012304 4 8 0.5941885554 -0.1405012287 Tabla 1: Autovalor menos estable en el caso de una tuber´ıa cuadrada A= 1 a Re = 100, β = 1, para diferentes valores de hyp. El valor de referencia para estos par´ametros [4, 5] es ω= 0.594177-0.140507 i. Figura 1: Componentes de la velocidad y presi´on del primer modo propio obtenido en un conducto triangular. Parte superior izquierda bu, superior derecha bv, inferior derecha bw, inferior izquierda bp Agradecimientos Este trabajo ha sido llevado a cabo gracias al proyecto concedido por el Ministerio de Educaci´on y Ciencia TRA2005-08983/TAIR cuyo investigador principal es D.Rafael G´omez Blanco. Referencias [1] L.Gonz´alez, V.Theofilis, Rafael G´omez Blanco. Finite-element numerical methods for viscous incompressible BiGlobal linear instability analysis on unstructured meshes. AIAA Journal 2007, vol.45 no.4 (840-854) [2] G.E.Karniadakis, S. Sherwin, Spectral/hp Element Methods for Computational Fluid Dynamics ,Oxford Science Publications, 2005. [3] Ch. Schwab. p and hp-Finite Element Methods. Theory and Applications in Solid and Fluid Mechanics. Oxford Science Publications, 2004. [4] T. Tatsumi, T. Yoshimura. Stability of the laminar flow in a rectangular duct. J. Fluid. Mech. 1990,vol.212 (437-449). [5] V. Theofilis, P. W. Duck, J. Owen. Viscous linear stability analysis of rectangular duct and cavity flows. J. Fluid. Mech. 2004,vol.505 (249-286). 8