NURBS-Enhanced FEM para problemas de scattering

XX Congreso de Ecuaciones Diferenciales y Aplicaciones X Congreso de Matem´ atica Aplicada Sevilla, 24-28 septiembre 2007 (pp. 1–8) NURBS-ENHANCED FEM PARA PROBLEMAS DE SCATTERING R. Sevilla, S. Fern´ andez-M´ endez y A. Huerta Laboratori de C`alcul Num`eric (LaC`aN) Departament de Matem`atica Aplicada III Universitat Polit`ecnica de Catalunya e-mail:{ruben.sevilla,sonia.fernandez,antonio.huerta}@upc.edu web page: http://www-lacan.upc.edu . Palabras clave: NURBS, elementos finitos, CAD, modelo geom´etrico exacto, elementos isoparam´etricos de alto orden, Discontinuous Galerkin Resumen En [1] se presenta el NURBS-Enhanced Finite Element Method (NEFEM), una mejora del cl´asico m´etodo de los elementos finitos que permite trabajar con el modelo geom´etrico exacto, independientemente de la discretizaci´on espacial utilizada. En [1] se desarrollan nuevas estrategias que permiten realizar la interpolaci´on y la integraci´on num´erica para los elementos que tienen una arista definida mediante Non-Uniform Rational B-Splines (NURBS) trimadas, ver [2]. En este trabajo se plantea la aplicaci´on del NEFEM, con una formulaci´on de Galerkin discontinuo (DG), para la resoluci´on de las ecuaciones de Maxwell. La importancia de trabajar con un buen modelo geom´etrico est´a justificada en aplicaciones de scattering de ondas electromagn´eticas, donde la geometr´ıa del obst´aculo es determinante para poder obtener resultados precisos. Mediante ejemplos num´ericos se observa que el NEFEM presenta importantes ventajas respecto a los cl´asicos elementos isoparam´etricos. Para una malla fijada el NEFEM es entre 11 y 15 veces m´as preciso. Por otro lado, para una precisi´on fijada, el NEFEM tambi´en resulta m´as eficiente, ya que s´olo requiere el 38 % de grados de libertad y el 75 % del tiempo de CPU consumido por el m´etodo DG est´andar. 1 R. Sevilla, S. Fern´andez-M´endez y A. Huerta 1. Introducci´on Son muchos los autores que en los ´ultimos a˜nos destacan la importancia de disponer de un buen modelo geom´etrico para la simulaci´on num´erica de problemas de contorno. En el contexto de m´etodos DG, ver [5], la importancia del modelo geom´etrico en la resoluci´on de las ecuaciones de Euler fue claramente demostrada en [3]. Utilizando interpolaciones lineales, la p´erdida de precisi´on cerca de contornos curvos es demasiado importante y acaba afectando al comportamiento global de la soluci´on. Por otra parte, en [4] se propone el denominado an´alisis isogeom´etrico. El objetivo fundamental de este planteamiento es trabajar con el modelo geom´etrico exacto, independientemente de la discretizaci´on espacial utilizada. La estrategia adoptada consiste en utilizar las funciones base de las NURBS, para describir la geometr´ıa de todo el dominio computacional y para aproximar la soluci´on. La metodolog´ıa presentada en [1], el NEFEM, comparte el objetivo principal del an´alisis isogeom´etrico, pero es m´as natural ya que la descripci´on mediante NURBS s´olo se utiliza para el contorno del dominio, es decir, lo que usualmente se obtiene con un programa est´andar de CAD (Computed Aided Design). De esta manera, el NEFEM considera la descripci´on geom´etrica exacta pero, a diferencia del an´alisis isogeom´etrico, la soluci´on se interpola de la manera habitual en elementos finitos, mediante funciones polin´omicas en cada elemento. En la mayor parte del dominio (elementos con lados rectos) se utilizan elementos finitos cl´asicos, mientras que para los elementos con una arista definida mediante NURBS es necesario definir la interpolaci´on y dise˜nar cuadraturas num´ericas adecuadas. El uso de la interpolaci´on cl´asica del m´etodo de los elementos finitos representa una gran ventaja frente al uso de una aproximaci´on funcional mediante NURBS: el NEFEM garantiza las buenas propiedades de los elementos finitos desde un punto de vista de eficiencia computacional y de convergencia. La secci´on 2 repasa brevemente los conceptos b´asicos del NEFEM. En la secci´on 3 se presenta una aplicaci´on del NEFEM en 2D. Se combina la utilizaci´on del m´etodo DG con la descripci´on del contorno mediante NURBS para la resoluci´on de problemas de scattering de ondas electromagn´eticas en dos dimensiones. Los resultados se comparan con los obtenidos mediante elementos isoparam´etricos de alto orden. Finalmente, en la secci´on 4 se presentan las conclusiones de este trabajo. 2. Conceptos b´asicos del NEFEM En esta secci´on se recuerdan los ingredientes b´asicos del NEFEM en 2D, ver [1] para los detalles. Se considera el contorno del dominio computacional descrito mediante curvas NURBS y, por tanto, es necesario prestar una atenci´on especial a la interpolaci´on y la integraci´on en los elementos con una arista definida por una curva NURBS trimada. Por otro lado, en los elementos interiores se utilizan elementos finitos cl´asicos. En esta secci´on se considera un dominio Ω ⊂R2, donde el contorno ∂Ω (o una parte de ´el) est´a definido mediante curvas NURBS. Cada NURBS est´a definida por una parametrizaci´on C: [0,1] −→ C([0,1]) ⊆∂Ω⊂R2. Se asume, tambi´en, una triangulaci´on del dominio Ω, tal que cada elemento tiene como m´aximo una arista sobre el contorno NURBS, ver la figura 1. 2 NEFEM para problemas de scattering W Figura 1: Dominio computacional con una parte del contorno definido mediante una curva NURBS (izquierda) y una triangulaci´on v´alida para el uso del NEFEM (derecha) Para sistematizar los c´alculos se considera la transformaci´on lineal Ψ que lleva los v´ertices del elemento Ωea los v´ertices vI= (−1,−1), vII = (1,−1) y vIII = (−1,1) del elemento de referencia Ie, ver figura 2. Ψ x y ξ η − 1 1 Ie Ωe Γ e − 11 Figura 2: Transformaci´on lineal, Ψ, para un elemento con una arista descrita mediante una NURBS trimada. Para la aproximaci´on funcional se considera una interpolaci´on polin´omica de grado m en el elemento Ie= Ψ−1(Ωe), con coordenadas ξ= (ξ, η) que, mediante la transformaci´on lineal, corresponde a una interpolaci´on polin´omica de grado men coordenadas cartesianas x= (x, y). Dado un conjunto de n+ 1 nodos ξi= (ξi, ηi) en el tri´angulo Ie, se considera la base formada por los polinomios de Lagrange, Li(ξ). Para sistematizar la evaluaci´on de los polinomios Li, para cualquier distribuci´on de nodos, se adopta la implementaci´on propuesta en [7]. En el proceso de c´alculo de las integrales de la forma d´ebil es necesario definir cuadraturas adecuadas para aproximar integrales sobre el contorno definido por NURBS e integrales en el interior de elementos con una arista definida por NURBS. Para cada elemento Ωecon una arista sobre el contorno descrito por NURBS se conoce el intervalo del espacio de par´ametros [λe 1, λe 2] tal que C([λe 1, λe 2]) = Γe, donde Γees la arista de Ωedefinida por la NURBS C, ver figura 2. Para integrar una funci´on fsobre la arista Γese aproxima ZΓe f d` =Zλe 2 λe 1 f(C(λ))|JC(λ)|dλ ≈ n X i=1 f(C(λi))|JC(λi)|ωi, 3 R. Sevilla, S. Fern´andez-M´endez y A. Huerta donde JCes la derivada de la NURBS y (λi, ωi) son los puntos y pesos de integraci´on en el intervalo [λe 1, λe 2] del espacio de par´ametros. Se han comparado diferentes cuadraturas y, desde el punto de vista de eficiencia computacional, la mejor opci´on consiste en utilizar cuadraturas de Gauss-Legendre simples, ver [1] para m´as detalles. Para la integraci´on en el interior de un elemento curvo Ωese han estudiado diferentes opciones, ver [1]. La mejor alternativa (desde el punto de vista de coste computacional) consiste en definir cuadraturas de Gauss en el rect´angulo [λe 1, λe 2]×[0,1], y adaptarlas a la geometr´ıa del elemento Iemediante la aplicaci´on ϕ= (ϕ1, ϕ2) : [λe 1, λe 2]×[0,1] →Ie, definida como ϕ1(λ, ζ) := φ1(λ)(1 −ζ)−ζ, ϕ2(λ, ζ) := φ2(λ)(1 −ζ) + ζ Mediante esta transformaci´on las integrales interiores para un elemento curvo Ωese calculan como ZΩe f dx dy =|JΨ|ZIe f dξ dη ≃ |JΨ| nu X i=1 nv X j=1 f(ξij)|Jϕ(λi, ζj)|ωi$j donde JΨes el jacobiano de la transformaci´on lineal Ψ, ξij =ϕ(λi, ζj), {λi, ωi}y{ζi, $i} son los puntos y pesos de integraci´on en los intervalos [λe 1, λe 2] y [0,1] respectivamente, y Jϕes el jacobiano de la transformaci´on ϕ. 3. Ejemplos num´ericos En esta secci´on se presentan dos ejemplos num´ericos para demostrar la eficiencia del NEFEM frente a la utilizaci´on de elementos isoparam´etricos de alto orden. En todos los ejemplos se ha considerado una formulaci´on DG. Las ecuaciones de Maxwell, utilizadas para la simulaci´on de problemas de scattering de ondas electromagn´eticas, se pueden escribir como un sistema hiperb´olico de primer orden ∂U ∂t +∂Fk(U) ∂xk =0,(1) con una definici´on adecuada del vector de cantidades conservadas Uy los vectores de flujos Fk. Se utiliza la notaci´on de Einstein, es decir, ´ındices repetidos indican un sumatorio en ese ´ındice. En dos dimensiones, las ecuaciones de Maxwell se simplifican y se obtienen dos sistemas desacoplados, denominados modos TE (Transverse Electric) y TM (Transverse Magnetic). Los modos TE y TM se pueden escribir como (1) definiendo adecuadamente el vector U y los flujos Fk. Por ejemplo, para el modo TM se define U=  µH1 µH2 ²E3  ,F1=  0 −E3 −H2  yF2=  E3 0 H1  . 4 NEFEM para problemas de scattering 3.1. Scattering por un c´ırculo PEC En primer lugar se considera el scattering causado por un c´ırculo de radio unidad formado por un material conductor el´ectrico perfecto (PEC). La onda incidente viaja de izquierda a derecha a velocidad unitaria y se toma como longitud de onda λ=1. En la figura 3 se ha representado el campo scattered despu´es de 4 ciclos completos y la Radar Cross Section (RCS), para el modo TM. La RCS es la cantidad de inter´es en problemas de scattering de ondas electromagn´eticas, y proporciona informaci´on de lo visible que es un objeto ante un radar, ver [8] para una descripci´on en detalle. Se observa una excelente correspondencia entre los resultados num´ericos y anal´ıticos en la RCS utilizando el NEFEM con una interpolaci´on funcional de grado p=7. −pi −3pi/4 −pi/2 −pi/4 0 pi/4 pi/2 3pi/4 pi 4 6 8 10 12 14 16 φ RCS Analítica NEFEM Figura 3: C´ırculo PEC, soluci´on NEFEM para λ= 1 con p=7: componente E3(izquierda) y RCS (derecha) para el modo TM En la figura 4 se comparan los resultados obtenidos mediante elementos isoparam´etricos de alto orden y el NEFEM, para el modo TM. A simple vista se observa una clara ventaja en la utilizaci´on del NEFEM. Utilizando polinomios de grado p=5 el NEFEM consigue capturar con una alta precisi´on la RCS en el rango φ∈[−π/3, π/3] mientras que utilizando elementos isoparam´etricos el error cometido en φ=π/3 es muy elevado. Si se aumenta el grado de los polinomios, es decir, considerando p=6, la soluci´on obtenida con el NEFEM y la soluci´on anal´ıtica son pr´acticamente id´enticas. Sin embargo, utilizando elementos isoparam´etricos con p=6 la soluci´on mejora sustancialmente, aunque fuera del rango [−π/4, π/4] los errores son muy elevados y se pueden observar a simple vista. De manera cuantitativa, con una medida del error en la RCS en norma L2, los resultados obtenidos con el NEFEM son un orden de magnitud m´as precisos que los obtenidos con elementos isoparam´etricos, para la misma discretizaci´on espacial. Adem´as, para este ejemplo, las ventajas del NEFEM, en t´erminos de eficiencia computacional, son muy claras. Para obtener una precisi´on fijada, por ejemplo un error de 10−2, el NEFEM necesita alrededor del 50 % de grados de libertad utilizado por los elementos isoparam´etricos. Esto supone que, para obtener el mismo error, el NEFEM necesita un grado de interpolaci´on p= 6 (2688 grados de libertad) mientras que con elementos isoparam´etricos es necesario 5 R. Sevilla, S. Fern´andez-M´endez y A. Huerta utilizar un grado de interpolaci´on p= 9 (5280 grados de libertad), ver figura 5. En esta figura se observa la evoluci´on del error en la RCS, en norma L2, al aumentar el grado de la aproximaci´on funcional. Al utilizar elementos isoparam´etricos de alto orden se converge a un resultado que no es f´ısicamente correcto, ver [6] para una discusi´on en detalle de la problem´atica de los elementos isoparam´etricos de alto orden. Para evitar los errores introducidos por la transformaci´on isoparam´etrica se ha utilizado una implementaci´on de los elementos de alto orden en coordenadas cartesianas, sin tener en cuenta la transformaci´on isoparam´etrica. Los resultados obtenidos muestran que para obtener una misma precisi´on, con elementos cartesianos se necesita un grado de interpolaci´on m´as que en NEFEM, ver figura 5. Puesto que el coste de los elementos cartesianos es comparable al coste de NEFEM la ventaja de la propuesta realizada en este trabajo es clara. −pi −3pi/4 −pi/2 −pi/4 0 pi/4 pi/2 3pi/4 pi 2 4 6 8 10 12 14 16 φ RCS Analítica FEM NEFEM −pi −3pi/4 −pi/2 −pi/4 0 pi/4 pi/2 3pi/4 pi 2 4 6 8 10 12 14 16 φ RCS Analítica FEM NEFEM Figura 4: C´ırculo PEC, soluci´on del modo TM para λ= 1 con: p=5 (izquierda) y p=6 (derecha) 0 1000 2000 3000 4000 5000 6000 7000 −3 −2.5 −2 −1.5 −1 −0.5 0 0.5 ndof log10(RCS error) EF Isoparametricos EF Cartesianos NEFEM Figura 5: C´ırculo PEC, soluci´on del modo TM para λ= 1: convergencia aumentando el grado p 6 NEFEM para problemas de scattering 3.2. Scattering por un perfil NACA0012 PEC Para finalizar, se considera un segundo ejemplo num´erico con un obst´aculo de geometr´ıa m´as compleja y una onda incidente de alta frecuencia y que viaja de abajo a arriba. Concretamente, el obst´aculo corresponde a un perfil NACA0012 de ala de avi´on, y la longitud de onda considerada es λ=0.2. −pi −3pi/4 −pi/2 −pi/4 0 pi/4 pi/2 3pi/4 pi −15 −10 −5 0 5 10 15 φ RCS Referencia NEFEM Figura 6: Perfil NACA0012 PEC, soluci´on NEFEM para λ=0.2 con p=14: componente E3(izquierda) y RCS (derecha) para el modo TM Para este ejemplo no se dispone de soluci´on anal´ıtica. Como soluci´on de referencia se toma la calculada en una malla refinada con una interpolaci´on funcional de alto orden. En la figura 6 se observa una muy buena correspondencia entre los resultados num´ericos y los resultados de referencia. Se observa tambi´en una buena correspondencia con los resultados publicados por otros autores, ver por ejemplo [9]. 4. Conclusiones En [1] se presenta una mejora del cl´asico m´etodo de los elementos finitos, el NURBSEnhanced Finite Element Method (NEFEM). Se propone la utilizaci´on del modelo geom´etrico exacto que proporciona un modelo de CAD y, por tanto, es necesario modificar la interpolaci´on y proponer t´ecnicas de integraci´on num´erica para los elementos que tienen una arista definida mediante curvas NURBS. En este trabajo se propone la utilizaci´on del NEFEM para la resoluci´on num´erica de las ecuaciones de Maxwell. En concreto, se plantea la resoluci´on de problemas de scattering de ondas electromagn´eticas. Mediante los ejemplos num´ericos se demuestra la ventaja de trabajar con el NEFEM en este tipo de problemas, obteniendo una precisi´on entre 11 y 15 veces mayor que los elementos isoparam´etricos, para una discretizaci´on espacial fijada. 7 R. Sevilla, S. Fern´andez-M´endez y A. Huerta Agradecimientos Los autores agradecen el financiamiento recibido del Ministerio de Educaci´on y Ciencia (DPI2004-03000) y del Departament de Matem`atica Aplicada III de la Universitat Polit`ecnica de Catalunya. Referencias [1] R. Sevilla, A. Huerta y S. Fern´andez-M´endez , NURBS-Enhanced Finite Element Method (NEFEM) Libro de Res´umenes (CEDYA 2005). Universidad Carlos III de Madrid, 2005. [2] L. Piegl, W. Tiller, The NURBS Book, Springer-Verlag, London, 1995. [3] F. Bassi and S. Rebay, High-order accurate Discontinuous Finite Element solution of the 2D Euler equations, J. Comput. Phys., v. 138, p. 251-285, 1997. [4] T. J. R. Hughes, J. A. Cottrell and Y. Bazilevs Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Comput. Methods Appl. Mech. Eng., v. 194, p. 4135-4195, 2005. [5] B. Cockburn, Discontinuous Galerkin methods for Computational Fluid Dynamics, in: E. Stein, R. de Borst, T.J.R. Hughes (Eds.), Fluids, Encyclopedia of Computational Mechanics, vol. 3, Wiley, New York, 2004, chapter 4. [6] Barna Szab´o, Alexander D¨uster and Ernest Rank , The p-version of the Finite Element Method, in: E. Stein, R. de Borst, T.J.R. Hughes (Eds.), Fundamentals, Encyclopedia of Computational Mechanics, vol. 1, Wiley, New York, 2004, chapter 5. [7] J. S. Hesthaven and T. Warburton, Nodal high-order methods on unstructured grids - I. Time-domain solution of Maxwell’s equations, J. Comput. Phys. 181 (1), 186–221, 2002. [8] C. A. Balanis, Advanced Engineering Electromagnetics, John Wiley and Sons, New York, 1989. [9] Jie Wu and Bo-nan Jiang, A least-squares finite element mehod for electromagnetic scattering problems, NASA Technical Report (1996). 8