En los problemas de estabilidad BiGlobal en un contexto fluidodinámico se han encontrado un cierto tipo de condicionantes que todo buen método numérico debe ser capaz de satisfacer. Por un lado y como condición más indispensable, se encuentra la libertad geométrica, es decir, el método debe ser capaz de discretizar cualquier dominio espacial de forma que no quede sujeto a formas simples o fácilmente transformables a estas. En segundo lugar se desea que las condiciones de contorno necesarias en los clásicos problemas de Mecánica de Fluidos sean fácilmente implementables. Por
 último que el orden del método elegido en la resolución del problema sea fácilmente variable en función de la precisión exigida. En este contexto se conoce que la precisión necesaria crece a través de un parámetro que no es otro que el número de Reynolds del problema, y por tanto es muy ventajoso que no se sea el tamaño de la malla el único grado de libertad para mejorar la precisión del cálculo, ver [2] G.E.Karniadakis, S. Sherwin, Spectral/hp Element Methods for Computational Fluid Dynamics, Oxford
 Science Publications,2005 y [3] Ch. Schwab. p and hp-Finite Element Methods. Theory and Applications in Solid and Fluid Mechanics.
 Oxford Science Publications, 2004 Con todos estos
 condicionantes se ha construido un código basado el elementos finitos de alto orden (hp-FEM) que ha sido aplicado al cálculo de los valores propios y los modos propios de algunos problemas clásicos de la Mecánica de Fluidos como son el movimiento de fluidos en un conductos.
En los casos aquí presentados, hemos aplicado el método de Arnoldi para el cálculo del espectro y los modos propios en los problemas de autovalores generalizados que aparecen. Como resultados del trabajo se muestran el cálculo del flujo base de una tubería rectangular y el posterior análisis de estabilidad BiGlobal de dicho flujo base.
 Los resultados han sido coherentes y en buena consonancia con los obtenidos mediante otros métodos en la literatura clásica, ver [1] L.González, V.Theofilis, Rafael Gómez Blanco. Finite-element numerical methods for viscous incompressible BiGlobal linear instability analysis on unstructured meshes. AIAA Journal 2007, vol.45 no.4 (840-854). 
Eliminar seleccionado

González Gutiérrez, Leo Miguel

Peña Alcaraz, Maite

Theofilis, Vassilis

idUS. Depósito de Investigación Universidad de Sevilla

XX Congreso de Ecuaciones Diferenciales y AplicacionesX Congreso de Matema´tica AplicadaSevilla, 24-28 septiembre 2007(pp. 1–8)Los elementos finitos de alto orden (hp-FEM) como me´todode ca´lculo en problemas de estabilidad fluidodina´mica.Maite Pen˜a Alcara´z1, Leo M. Gonza´lez2, VassilisTheofilis31 ETS de Ingenier´ıa ICAI, Universidad Pontificia Comillas de Madrid, Alberto Aguilera 23, 28015Madrid. E-mail: maitepalcaraz@hotmail.com.2 Dpto. de Ensen˜anzas Ba´sicas de la Ingenier´ıa Naval, ETSIN-UPM, Arco de la Victoria s/n, 28040Madrid. E-mail: leo.gonzalez@upm.es.3 ETS Ingenieros Aeronau´ticos, Univ. Polite´cnica de Madrid, Plaza Cardenal Cisneros 3, 28040 Madrid.E-mail: vassilis@torroja.dmt.upm.es.Palabras clave: Bi-Global instability, high order methods, hp-FEMResumenEn los problemas de estabilidad BiGlobal en un contexto fluidodina´mico se hanencontrado un cierto tipo de condicionantes que todo buen me´todo nume´rico debe sercapaz de satisfacer. Por un lado y como condicio´n ma´s indispensable, se encuentra lalibertad geome´trica, es decir, el me´todo debe ser capaz de discretizar cualquier dominioespacial de forma que no quede sujeto a formas simples o fa´cilmente transformablesa e´stas. En segundo lugar se desea que las condiciones de contorno necesarias enlos cla´sicos problemas de Meca´nica de Fluidos sean fa´cilmente implementables. Poru´ltimo que el orden del me´todo elegido en la resolucio´n del problema sea fa´cilmentevariable en funcio´n de la precisio´n exigida. En este contexto se conoce que la precisio´nnecesaria crece a trave´s de un para´metro que no es otro que el nu´mero de Reynolds delproblema, y por tanto es muy ventajoso que no se sea el taman˜o de la malla el u´nicogrado de libertad para mejorar la precisio´n del ca´lculo, ver [2] y [3]. Con todos estoscondicionantes se ha construido un co´digo basado el elementos finitos de alto orden(hp-FEM) que ha sido aplicado al ca´lculo de los valores propios y los modos propiosde algunos problemas cla´sicos de la Meca´nica de Fluidos como son el movimiento defluidos en un conductos.En los casos aqu´ı presentados, hemos aplicado el me´todo de Arnoldi para el ca´lculodel espectro y los modos propios en los problemas de autovalores generalizados queaparecen. Como resultados del trabajo se muestran el ca´lculo del flujo base de unatuber´ıa rectangular y el posterior ana´lisis de estabilidad BiGlobal de dicho flujo base.Los resultados han sido coherentes y en buena consonancia con los obtenidos medianteotros me´todos en la literatura cla´sica, ver [1].1Maite Pen˜a Alcara´z, Leo M. Gonza´lez, Vassilis Theofilis1. Introduccio´nLa descripcio´n cla´sica de la teor´ıa de estabilidad BiGlobal se basa en una perturbacio´nlineal de las ecuaciones de Navier-Stokes alrededor de un flujo base (U, V ,W ), el cual esuna solucio´n particular de las mismas. El ca´lculo de esta solucio´n particular (U, V ,W ) seha obtenido resolviendo simplificaciones bidimensionales o unidimensionales de las ecua-ciones de Navier-Stokes mediante una discretizacio´n de alto orden hp-FEM. Una vez seha calculado el flujo base se ha linealizado el problema entorno a dicha solucio´n medianteel ansatz:ui(x, y, z, t) = U i(x, y) + ²ûi(x, y)eiωteiβz (1)p(x, y, z, t) = P (x, y) + ²p̂(x, y)eiωteiβz (2)Considerando ²¿ 1, tenemos que (ûi, p̂) representan las amplitudes de la perturbacio´n,el nu´mero complejo ω contiene como parte imaginaria la tasa de crecimiento y como partereal la frequencia de dicha perturbacio´n, por u´ltimo β = 2piLz y Lz son el nu´mero de onday la longitud de onda de la perturbacio´n, respectivamente. Al introducir este ansatz enlas ecuaciones de Navier-Stokes y considerar el flujo base como una solucio´n particular sellega a un problema de autovalores generalizado del tipo:Auˆ1uˆ2uˆ3pˆ = −ωBuˆ1uˆ2uˆ3pˆ . (3)Donde el operador complejo y no sime´trico A se define como:A =a11∂u¯1∂y 0∂∂x∂u¯2∂x a22 0∂∂y∂u¯3∂x∂u¯3∂y a33 iβ∂∂x∂∂y iβ 0 . (4)aii ={u¯j∂∂xj+∂u¯i∂xi− 1Re(∂2∂x2j− β2)+ iβu¯3}, j = 1, 2. (5)y no se considerara´ el convenio de indices repetidos de Einstein para la ecuacio´n (5).Por otro lado el operador real sime´trico B se define como:B =M 0 0 00 M 0 00 0 M 00 0 0 0 , (6)donde M representa la cla´sica matriz de masas de elementos finitos.2Normas para la edicio´n (CEDYA 2007)En la ecuacio´n (3) el valor propio viene dado por ω y el autovector (û, v̂, ŵ, p̂) contienelas amplitudes de las perturbaciones.La resolucio´n directa del problema de autovalores generalizado mediante un me´todotipo QZ implicar´ıa un coste computacional excesivo para matrices A y B de gran dimen-sio´n, es por ello que la metodolog´ıa escogida para la resolucio´n de este tipo de problemasha sido un me´todo iterativo basado en la creacio´n de un subespacio de Krylov conocidocomo me´todo de Arnoldi. Para ello en cada iteracio´n se ha resuelto un problema linealque implicaba a la matriz A como matriz de coeficientes y que por tanto obligaba a unadescomposicio´n LU de la misma al principio del ca´lculo.2. Discretizacio´n espacial del problema. Elementos finitosde alto orden: hp-FEMLa resolucio´n nume´rica del problema de autovalores generalizado descrito anterior-mente, requiere una discretizacio´n espacial que se hara´ mediante el uso de una base defunciones, Φq(x). Algunas versiones de este tipo de ca´lculos y discretizaciones cuando lasfunciones base son lagrangianas, nodales y de bajo orden se puede ver en [1]. De cualquiermodo la experiencia con funciones de bajo orden muestra que para resolver estructurasa nu´meros de Reynolds moderados se requiere de mallas muy finas con el consecuenteconsumo de memoria y tiempo de computacio´n.En el presente trabajo se propone una discretizacio´n modal de alto orden basada enpolinomios de Legendre Lp(x). El hecho de que sea modal implica que el valor de loscoeficientes de la expansio´n no tiene interpretacio´n f´ısica, por otro lado sabemos quelas expansiones modales tienen estructura jera´rquica, lo que implica que un conjunto defunciones hasta orden P − 1 esta´ incluida en una expansio´n de orden P . Una de laspropiedades fundamentales de este tipo de expansiones es la ortogonalidad, que junto conel hecho de poseer estructura jera´rquica da lugar a matrices bien condicionadas [2].Para problemas como los aqu´ı presentados que implican un segundo orden de difer-enciacio´n basta con garantizar que la solucio´n se halla en H1. T´ıpicamente, el me´todode los elementos finitos resuelve este problema mediante la imposicio´n de una condicio´nde continuidad C0 entre elementos vecinos, ello implica que las expansiones globales soncontinuas en todo el dominio mientras que la continuidad en la derivada se consigue conla convergencia, ver [2].En un intervalo esta´ndar Ω = {ξ| − 1 ≤ ξ ≤ 1} una expansio´n modal de orden p sedenota ψp (x) y se define comoψp (ξ) =1− ξ2p = 0[1− ξ2] [1 + ξ2]Lp−1(ξ) 0 < p < P1 + ξ2p = P(7)En caso de tener elementos que no coinciden con el elemento esta´ndar [−1, 1], sera´ posi-3Maite Pen˜a Alcara´z, Leo M. Gonza´lez, Vassilis Theofilisble transformarlos a dicho intervalo mediante una transformacio´n lineal. Por lo tantou´nicamente definiremos nuestras funciones base para el intervalo esta´ndar. Del mismomodo tambie´n es necesario sen˜alar que dado que los elementos finitos requieren una for-mulacio´n tipo Galerkin es necesario elegir un me´todo de integracio´n nume´rica que ene´ste caso sera´ el de Gauss. Para una situacio´n multidimensional se puede generalizar estametodolog´ıa, y ma´s concretamente en un caso bidimensional la discretizacio´n se puederealizar tanto mediante elementos de tipo triangular o de tipo cuadrila´tero. De este modolas funciones base se pueden expandir a ma´s dimensiones mediante un producto tensorialde funciones base unidimensionales ψpq(ξ1, ξ2) = ψp(ξ1)ψq(ξ2).En el caso de un cuadrila´tero esta´ndar, las coordenadas cartesianas (ξ1, ξ2) quedandefinidas mediante contornos constantes Q2 = {(ξ1, ξ2)| − 1 ≤ ξ1, ξ2 ≤ 1}.Por el contrario, para un elemento triangular esta´ndar, donde hay una dependenciaentre las coordenadas cartesianas (ξ1, ξ2) a la hora de definir el dominio, es decir T 2 ={(ξ1, ξ2)| − 1 ≤ ξ1, ξ2, ξ1 + ξ2 ≤ 0}. Un conjunto de funciones base definidas en elemntostriangulares ha sido desarrollada por Karniadakis & Sherwin [2], presentando la ventaja deposeer un sistema de coordenadas en el que no hay dependencia entre las coordenadas parala definicio´n del dominio. La gran ventaja de ello es la posibilidad de realizar expansionesmultidimensionales mediante productos tensoriales de las elaboradas en una dimension.La forma de trabajar con elementos triangulares que aqu´ı se propone, ver [2], se basaen una transformacio´n de coordenadas, mediante el cambio:η1 = 21 + ξ11− ξ2 − 1, η2 = ξ2. (8)Por lo tanto, un paso necesario para realizar un ca´lculo mediante elementos finitos hp-FEM require de diversas transformaciones de los elementos geome´tricos a los elementosesta´ndar [−1, 1] en el caso 1-D, Q2 para el caso 2-D y elementos tipo cuadrila´teros y T 2para el caso 2-D y elementos tipo triangular.2.1. Estudio de la Forma de las Matrices de Masas y RigidezA consecuencia de la ortogonalidad de los polinomios de Legendre Lp−1 (ξ) las matricesde masas y rigidez son, respectivamente, tri- y penta- diagonales. Este tipo de discretizacio´nmediante elementos finitos tipo hp−(FEM) presenta en este sentido una gran ventajafrente a discretizaciones tipo espectral Chebyshev, dado que en este u´ltimo caso la matrizde rigidez es densa.Se consideran las funciones definidas por la ecuacio´n (7), y se definen las matricessime´tricas M = {Mij}i,j=0,...,P y R = {Rij}i,j=0,...,P dondeMij =∫ 1−1ψiψjdξ Rij =∫ 1−1∂ψi∂ξ∂ψj∂ξdξSe consideran conocidas las siguientes propiedades de los polinomios de Legendre:4Normas para la edicio´n (CEDYA 2007)∫ 1−1LpLqdξ =22p+ 1δp,q, (9)(k + 1)Lk+1(ξ) = (2k + 1)ξLk(ξ)− kLk−1(ξ), (10)(2k + 1)Lk(ξ) = L′k+1(ξ)− L′k−1(ξ), (11)gr(Lk) = k (12)siendo gr(p) el grado del polinomio p(ξ) y se denotara´ ∂f∂ξ con f′ o f ′(ξ).Por otro lado se sabe que si Lp(ξ) es el polinomio de Legendre de orden p, entoncessera´ un polinomio par si p es par e impar si p es impar. Del mismo modo si Lp(ξ) es unpolinomio par, entonces L′p(ξ) es impar o viceversa.Teorema. Mij = 0 si i+ j es impar, con i, j 6= 0, PDemostracio´n. Si i+ j es impar, i sera´ impar o j sera´ impar. Se puede suponer porsimetr´ıa en el razonamiento que i es impar y j par. Se tiene queMij =M2i′−1,2j′ =∫ 1−1(1− ξ24)2L2i′L2j′−1dξ = 0por ser producto de dos funciones pares y una impar, y por tanto ser una funcio´n imparintegrada entre −1 y 1, por lo que la integral se anula.Teorema. Rij = 0 si i+ j es impar, con i, j 6= 0, PDemostracio´n.Si i, j 6= 0, P , ψk(ξ) = 1−ξ22 Lk−1(ξ) con k = i, j por lo que es inmediatoque si k es impar, entonces ψk es impar y viceversa. Por tanto ψ′k sera´ par si y so´lo si k loes y sera´ impar en caso contrario. Por tanto un razonamiento ana´logo al anterior permiteconcluir que si i+ j es impar, tambie´n lo sera´ el producto ψiψj y por tanto la integral seanulara´.Puesto que ψ′0 =12 y ψ′P = −12 se obtiene de forma inmediata que Ri0 = RiP = 0 si ies impar.Para el ca´lculo de las matrices de masas y rigidez a partir del grado de los polinomios,tenemos la siguiente proposicio´n:Proposicio´n. Para cualquier polinomio de Legendre Lp(ξ) se tiene que∫ 1−1Lp(ξ)p(ξ)dξ = 0si p(ξ) es un polinomio de grado g < p.Demostracio´n.Si p(ξ) es un polinomio de grado menor que p entonces p(ξ) se puede expresar comocombinacio´n lineal de polinomios de Legendre de grado estrictamente menor que p. Portanto p(ξ) =∑αiLi(ξ) y ello implica que∫ 1−1Lpp(ξ)dξ =∫ 1−1∑αiLiLpdξ =∑αi∫ 1−1LiLpdξ =∑αiδi,p22p+ 1= 05Maite Pen˜a Alcara´z, Leo M. Gonza´lez, Vassilis TheofilisTeorema.∫ 1−1 ψiψjdξ = 0 si |i− j| > 4, i, j 6= 0, PDemostracio´n.Se puede suponer sin pe´rdida de la generalidad que i ≥ j. Entonces∫ 1−1ψiψjdξ =∫ 1−1(1− ξ22)2Lj−1Li−1dξ =∫ 1−1p(ξ)Li−1dξdonde p(ξ) es un polinomio de grado j − 1+ 4 = j − 3. Por tanto si i− j > 4 la prueba seacaba utilizando el resultado de la proposicio´n anterior.En el caso particular en que i o j sean igual a 0 o P se obtiene el siguiente resultadocuya demostracio´n de este teorema es ana´loga a la del teorema anterior:∫ 1−1ψ0ψjdξ =∫ 1−1ψPψjdξ = 0 si j > 4, j 6= 0, PAhora se ha terminado el ana´lisis de todos los te´rminos de la matriz de masas ex-ceptuando la diagonal principal, dos diagonales secundarias a distancia 2 y 4 y algunoste´rminos dispersos como i = 0, P . Para terminar un simple ca´lculo permite asegurar queen general ningu´n otro te´rmino de la matriz ha de ser nulo.Para la matriz de rigidez el razonamiento es muy parecido, basta considerar el siguiente:Teorema.∫ 1−1∂ψi∂ξ∂ψj∂ξ dξ =∫ 1−1∂2ψi∂ξ2ψjdξ ∀i, j 6= 0, PDemostracio´n.Basta integrar por partes sabiendo que la funcio´n ψi(−1) = ψi(1) = 0.Una vez que se tiene este resultado, la forma de la matriz de rigidez se obtiene deforma inmediata.Teorema.∫ 1−1∂ψi∂ξ∂ψj∂ξ dξ = 0 si i, j 6= 0, P y se cumple que i > j + 2 o viceversa.Demostracio´n.Usando que la integral del corolario es igual∫ 1−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 sera´ nulacomo ya se ha probado.Por u´ltimo, falta analizar la forma de los te´rminos de la matriz de rigidez cuandoalguno de los ı´ndices es igual a 0, P . En estos casos se cumple el siguiente resultado:Teorema.∫ 1−1∂ψi∂ξ∂ψj∂ξ dξ = 0 ∀j 6= 0, P , i = 0, PDemostracio´n.Basta realizar un simple ca´lculo para comprobar que |∂ψi∂ξ | = 12 , por lo que la integral,queda como:∫ 1−1∂ψj∂ξ dξ = k · (ψj(1)− ψj(−1)) = 0.Por tanto, los u´nicos 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 y R0,0, R0,P , RP,0, RP,P .Como resumen, en las dos siguientes matrices se ha dibujado la forma de los elementosno nulos de M y R:6Normas para la edicio´n (CEDYA 2007)M =• • • • • •• • • • •• • • • •• • • • • •• • • • • •• • • • •• • • • •• • • •• • • •• • •• • •• • • • • •R =• •• •• •• • •• • •• • •• • •• •• •• •3. Aspectos nume´ricosPara asegurar la estabilidad de una discretizacio´n mediante el me´todo de los elementosfinitos de las ecuaciones de Navier-Stokes o problemas derivados de ellas, se sabe que lacondicio´n de compatibilidad inf −sup debe ser satisfecha entre aquellos espacios donde sediscretiza la velocidad y la presio´n. En el presente trabajo el nu´mero de modos empleadospara la discretizacio´n de la presio´n p − 1 ha sido rebajado en una unidad respecto al develocidad p.La dimension m del subespacio de Krylov para todos los calculos que se presentan acontinuacio´n ha sido de 100. Por otro lado, el nu´mero de elementos de matriz complejosque se reservan para el ca´lculo del operador A son 4p4N2e . De este modo las estimacionesde memoria RAM en Mb necesaria para llevar a cabo estos ca´lculos ha sido64p4N2e10242= 64p4(Ne1024)2=(8p2Ne1024)2(13)A consecuencia de todo esto el problema ma´s grande resoluble en un PC con 4 Gb dememoria RAM viene determinado por la condicio´n Nep2 ≤ 8000.4. Aplicacio´n.En este caso se considerara´ un flujo base unidimensional (0, 0,W ) creado en una tu-ber´ıa de seccio´n triangular o cuadrada que se extiende a lo largo del eje Z con gradiente depresio´n favorable y constante. El caso triangular el dominio estara´ formado por los pun-tos interiores al tria´ngulo definido por los puntos Γb = {(1,−1), (−1, 1), (−1,−1)}. Paraobtener el flujo base se resuelve nume´ricamente la ecuacio´n de Poisson,∇2W¯ (x, y) = −2, W¯ (x, y)Γb = 0, (14)El ana´lisis de dicho flujo base contiene, como se puede ver en las ecuaciones (4)-(5),una pareja de para´metros libres: Re y β. En nuestro caso se han realizado ana´lisis paraun amplio rango de valores de β ∈ (0, 25) y Re ≤ 1000. Por otra parte, en el caso de unatuber´ıa rectangular tambie´n se podr´ıa hablar de un tercer para´metro geome´trico como esel aspecto geome´trico A, que no es otra cosa que el cociente entre los lados del conducto.La convergencia en el ca´lculo de los autovalores se consigue tanto mediante un aumentodel nu´mero de elementos h, como del orden polino´mico p, ver tabla 1.En la figura 1 se puede puede ver la forma de las componentes de los modos asociadosal autovalor menos estable en el caso de tener contorno triangular.7Maite Pen˜a Alcara´z, Leo M. Gonza´lez, Vassilis Theofilish p ωr ωi h p ωr ωi1 8 0.5991425642 -0.1490564284 2 4 0.5982952805 -0.13741494141 12 0.5941967173 -0.1404610035 2 8 0.5941884517 -0.14050279821 16 0.5941885536 -0.1405012852 4 4 0.5941712455 -0.14047872891 20 0.5941885590 -0.1405012304 4 8 0.5941885554 -0.1405012287Tabla 1: Autovalor menos estable en el caso de una tuber´ıa cuadrada A = 1 a Re =100, β = 1, para diferentes valores de h y p. El valor de referencia para estos para´metros[4, 5] es ω = 0.594177-0.140507 i.Figura 1: Componentes de la velocidad y presio´n del primer modo propio obtenido en unconducto triangular. Parte superior izquierda û, superior derecha v̂, inferior derecha ŵ,inferior izquierda p̂AgradecimientosEste trabajo ha sido llevado a cabo gracias al proyecto concedido por el Ministeriode Educacio´n y Ciencia TRA2005-08983/TAIR cuyo investigador principal es D.RafaelGo´mez Blanco.Referencias[1] L.Gonza´lez, V.Theofilis, Rafael Go´mez Blanco. Finite-element numerical methods for viscous incom-pressible 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 ,Ox-ford Science Publications, 2005.[3] Ch. Schwab. p and hp-Finite Element Methods. Theory and Applications in Solid and Fluid Mechan-ics. 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 cavityflows. J. Fluid. Mech. 2004,vol.505 (249-286).8

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

https://idus.us.es/xmlui/bitstream/handle/11441/35798/Los%20elementos%20finitos%20de%20alto%20orden%20%28hp-FEM%29%20como%20m%c3%a9todo%20de%20c%c3%a1lculo.pdf?sequence=1&isAllowed=y

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

Abstract

Similar works

Full text

Available Versions

idUS. Depósito de Investigación Universidad de Sevilla