1. Introducción

El Método de los Elementos Finitos (FEM) permite obtener una solución numérica aproximada sobre un dominio (medio continuo), sobre el que están definidas ciertas ecuaciones diferenciales en forma débil o integral, que caracterizan el comportamiento físico del problema modelado. Es un método muy usado debido a su generalidad y a la facilidad de introducir dominios de cálculo complejos. Además, es fácilmente adaptable a problemas de diversos tipos, como los problemas de transmisión de calor.

Aunque el FEM es un método bastante versátil y ha sido ampliamente desarrollado desde mediados de la década de los 50 del siglo pasado ^{[1,2,3,4,5,6],} sus características no le permiten representar adecuadamente el comportamiento real de ciertos tipos de problemas. Por ejemplo, aquellos que involucran discontinuidades o singularidades, como el problema de transmisión de calor en un dominio no homogéneo, en donde los subdominios o fases, están separadas por una interfase, y poseen condiciones (o conductividades) distintas. En estos casos, es necesario utilizar una malla de elementos finitos que esté ajustada a las aristas de los elementos y a la interfase, a fin de poder capturar dicha singularidad; sin embargo, esto implica un mayor costo computacional (ya que, dependiendo de la disposición de la interfase, varía la malla a usar), sobretodo en problemas de tipo transitorio. En base a esto, surgen nuevas ideas para resolver este problema, evitando el remallado y ahorrando costo computacional. Es entonces, cuando a finales de la década de los 90 y principios del año 2000 ^{[7,8,9,10,11,12,13],} surge un nuevo método llamado el Método de los Elementos Finitos Extendido (XFEM). Este método resulta útil para aproximar de manera eficiente singularidades o discontinuidades asociadas a fracturas, fisuras o cambio de fases, mediante la adición de funciones especiales, que enriquecen el espacio solución de la aproximación de FEM, todo esto independientemente de la malla de elementos finitos que se esté usando. De esta forma, el XFEM constituye un importante avance en la resolución de problemas con interfase, y es hoy en día un área donde se realiza un amplio trabajo, tanto en definir estimadores de error y procesos adaptativos ^[14], o en la utilización de elementos no tradicionales ^[15]. No obstante, cuando se aplica XFEM a los problemas de difusión en un sistema de dos fases con condiciones de materiales muy distintas, éste produce una representación inexacta de los flujos en la vecindad de la interfase, debido a que el enriquecimiento XFEM mejora la calidad global de la solución, pero no satisface algunos rasgos locales de los flujos. Con el fin de remediar este inconveniente, surge un nuevo método, llamado XFEM+, el cual agrega una restricción adicional a la formulación XFEM permitiendo que se reproduzcan con mayor precisión los rasgos de los flujos locales en la zona de transición de las fases. Este nuevo método fue publicado en ^[16], y su efectividad queda comprobada mediante varios ejemplos numéricos.

El aporte de este trabajo, radica en la implementación de dichos métodos para elementos cuadriláteros. Los resultados numéricos obtenidos para este tipo de elementos, y distintas interfases, dejan ver la notable mejora del XFEM al aumentar el orden de interpolación, en comparación con los resultados obtenidos al implementar XFEM+ con elementos triangulares.

2. Planteamiento del problema

Se considera un dominio Ω ⊂ ℝ², abierto y acotado con frontera suave dividida en dos partes por la interfase a, que será una curva suave que completará el dominio Ω1 ∪ Ω2 ∪ α (ver Figura 1). En este dominio se desea resolver el siguiente problema:

Figura 1
Dominio Ω = Ω₁ ∪ Ω₂ ∪α fases o subdominios Ω_i, interfase α, y n vector normal a α.

Donde: u incógnita, f є L²(Ω) el término fuente, u₀ condición de frontera tipo Dirichlet establecida sobre ∂Ω, es el salto del flujo normal y v el parámetro que define los materiales

Se define el flujo mediante la expresión siguiente:

Asimismo, se define para cualquier función suave u en Ω el salto de u sobre α como:

donde es la restricción de u sobre Ω_i, con i = 1,2. Ahora, para cualquier punto arbitrario P є α, se introducen los flujos a ambos lados de la interfase:

Al denotar como la componente normal del gradiente, en donde n =(n₁, n₂) es un vector unitario normal a P є α, se tiene que el flujo normal está dado por:

Se debe resaltar que, la solución de la ecuación posee una característica de discontinuidad del gradiente sobre la zona de la interfase, es por esto que dicha ecuación diferencial no está definida sobre la interfase por no tener sentido en esos puntos donde el parámetro v es discontinuo. Sin embargo, sobre la interfase se ha de verificar la continuidad de la componente normal del flujo, lo cual se obtiene con la última condición impuesta al problema (1). Además, a mayor diferencia entre los valores de v₁ y v₂, mayor será el salto del gradiente; por ello la discontinuidad del gradiente (en la dirección normal a la interfase) es un rasgo característico que posee la solución exacta del problema que se está estudiando, y por tanto, cualquier aproximación debería contener esta propiedad.

3. El método de los elementos finitos extendido: XFEM

3.1. Descripción de fases usando las funciones: Level Set y Ridge

Para facilitar el desarrollo teórico de dicho método, se hace necesario introducir las siguientes funciones¹:

3.1.1. Level Set

Denotada por ϕ, esta función permite establecer donde está ubicado cada punto del dominio, respecto a la interfase α. Para efectos prácticos, ϕ se toma como la función distancia a la interfase, con signo:

A menudo esta función ϕ se trunca a una cierta distancia, de tal manera que a partir de ese momento ϕ se toma constante (ver Figura 2).

Figura 2
La función Level Set es, en general, la función distancia signada a la interfase y a partir de cierta distancia se toma constante.

Recuerde que el dominio Ω se divide en dos subdominios Ω₁ y Ω₂. A los elementos Ωⁱ de la malla de elementos finitos, que contienen la interfase en su interior, se les llamará elementos enriquecidos, y a los demás elementos estándar (ver Figura 3).

Figura 3
Bosquejo en ID, de los elementos estándar, elementos enriquecidos y Level Set.

3.1.2. Función de Ridge

La función de Ridge a utilizar fue introducida en ^[11], denotada de ahora en adelante R, y definida como:

donde:

• Nnodo es el número total de nodos que posee la malla a usar.

• ϕ_i es la función Level Set evaluada en el nodo x_¡.

• N_¡ son las funciones de forma de FEM.

Es de acotar que esta función se anula por completo sobre los elementos estándar de la malla; y en los elementos enriquecidos, se anula solo en los nodos, además presenta una discontinuidad del gradiente, exactamente sobre la interfase (Ver Figura 4), lo cual es exactamente el comportamiento de la solución analítica del problema. Por ello, la necesidad de introducir esta función para enriquecer el espacio de las funciones de forma, y así obtener una solución mucho más cercana a la exacta.

Figura 4
Ridge y Level Set en 2D (para elementos cuadriláteros).

donde, a_j son los valores nodales en los nodos asociados a elementos enriquecidos. Por otra parte, los valores nodales ϕ_i , ¡ = 1, ..,N_nodo permiten interpolar ϕ en cada elemento, de la siguiente forma:

De esta manera, la descripción discreta de la interfase es una curva continua, suave en el interior de los elementos y con discontinuidades en la pendiente en los lados de los elementos enriquecidos (ver Figura 4).

Esto permite también calcular el vector normal unitario a la interfase, en un punto P є α interior a un elemento Ωⁱ:

El XFEM enriquece la solución introduciendo la posibilidad de que aparezca un salto en la componente normal del gradiente para intentar satisfacer sobre α. Así, el gradiente de u_X resulta ser

Note que

Es decir, por tanto el salto del gradiente de R es paralelo a la normal n ya que esta es paralela a Gracias a esta propiedad, se espera que el XFEM sea capaz de obtener una mejora significativa en la aproximación de la solución exacta.

3.2. Aproximación del flujo en la zona de la interfase

El FEM no puede satisfacer la condición de continuidad del flujo normal en los puntos de la interfase, ya que el gradiente es continuo en el interior de los elementos. Sin embargo, en la implementación numérica de XFEM, se observa un fenómeno: a pesar de la capacidad que tiene el método para capturar el salto del gradiente, cuando los parámetros v₁ y v₂ son altamente distintos, el método es poco preciso y arroja una solución no tan buena de la aproximación del flujo en la zona cercana a la interfase.

3.3. Continuidad del flujo normal sobre la interfase

Para imponer de manera explícita la continuidad del flujo normal al sistema XFEM, se propone añadir en algunos de los puntos P є α la ecuación. Para esto, en cada punto fronterizo se añade una condición al sistema general. De ahora en adelante, se va a considerar la posición de P (en el segmento de interfase dentro del elemento enriquecido Ω^k) de forma parametrizada, dada por un parámetro escalar λ є [0,1], de la siguiente manera:

donde P₁y P₂ son los puntos de intersección entre la interfase α y los lados del elemento enriquecido Ω^k (Figura 5).

Figura 5
Un elemento cuadrilátero enriquecido Ω^k con un punto P en la interfase, donde se va a imponer de manera fuerte la continuidad del flujo normal.

Para lograr dicha imposición, se presenta a continuación la expresión algebraica en el contexto de la aproximación u_X de XFEM.

Denotado por

Se obtiene que:

Por otro lado, la expresión del salto de flujos a ambos lados de un punto está dada por

donde y la diferencia de los flujos a ambos lados de la interfase q₁- q₂ se puede escribir como:

Observe que la condición (2) se puede imponer en cada punto P є α (en un elemento enriquecido Ω^k ) y representa una restricción algebraica sobre las incógnitas u_iy Recuerde también que la solución XFEM, u_X, no verifica esta restricción puesto que sólo se impone de manera débil.

3.4. Restricción algebraica

Como se dijo anteriormente, la restricción se impone sobre puntos Estratégicamente, se tomaran estos puntos P como los puntos de Gauss para aproximar la integral de linea sobre Además en la práctica, dependiendo del tipo de elemento (triangular o cuadrilátero) será más apropiado tomar 2 o 3 puntos P: Recordemos que, para elementos triangulares, las funciones de forma a usar son lineales:

Por tanto, es adecuado tomar solo 2 puntos ya que, como es sabido, una recta viene determinada por únicamente dos puntos.

Para los elementos cuadriláteros, las funciones de forma a usar son bilineales:

Por tal motivo, resulta apropiado tomar solo 3 puntos pues, en este caso, una parábola viene determinada por únicamente tres puntos.

Por otro lado, es importante destacar que la interfase puede atravesar de diversas maneras un elemento enriquecido (bien sea triangular o cuadrilátero), y según la forma en que lo haga, varían las coordenadas de los puntos de intersección P₁ y P₂, y por tanto, varían los vectores de funciones de forma asociados a estos puntos (y denotados N^[k] (P₁) y N^[k] (P₂) respectivamente). Para encontrar la forma explícita de estos vectores, basta con utilizar la condición de P₁ y P₂ de pertenencia a la interfase α:

4.1. Resultados numéricos

4.1. Para una interfase recta

Sea Ω = ((0,1) x (0,1)), Γ_b el lado inferior de Ω, y Γ_t, el lado superior de Ω. El problema fuerte a resolver es:

donde v₁ = 1 y v₂ 1000, con una interfase α inclinada 10°.

Es importante acotar que, la efectividad y precisión del método numérico se llega a visualizar mejor al momento de graficar los flujos (ya que permite ver si se cumple la continuidad del flujo normal sobre α), es por ello que se muestra solamente los flujos obtenidos de la solución numérica.

En este ejemplo, se considera un caso bastante extremo puesto que los valores de v₁ y v₂ son altamente distintos. Observando la Figura 6 se puede apreciar que (tanto en la malla triangular como en la malla de cuadriláteros) FEM falla por completo, perdiendo totalmente la precisión de los flujos sobre a (Figuras 6a y 6b), y haciendo la comparación con los flujos obtenidos de XFEM (Figuras 6c y 6d), se aprecia que dicho método empieza a fallar, perdiendo la precisión en cuanto a la continuidad del flujo normal sobre la interfase a y sus alrededores.

Figura 6
Flujos numéricos obtenidos con los valores v₁= 1 y v₂ = 1000.

Sin embargo, es de acotar que, al utilizar elementos cuadriláteros, se está aumentando el orden de la aproximación (de funciones lineales a funciones cuadráticas), por ello, el XFEM con cuadriláteros tiene menores errores en la aproximación de flujos globales a comparación con el XFEM con elementos triangulares. Finalmente, esta incapacidad de XFEM de reproducir con exactitud los flujos locales, se debe a que, la condición sobre α del problema original se impuso de manera débil en la formulación variacional, y recordemos que los métodos numéricos que se están usando, resuelven el problema débil. Por tanto, se hace necesario imponer de manera fuerte sobre el método, dicha condición de continuidad del flujo normal sobre α, y esta necesidad es la que da origen a una modificación de XFEM, conocida como XFEM+.

Comparando los resultados obtenidos con XFEM+ (con malla triangular, Figura 6e), y XFEM (con malla de cuadriláteros, Figura 6d), podemos observar que XFEM con cuadriláteros alcanza una precisión muy similar a la de XFEM+ con triángulos, en cuanto a la aproximación de los flujos globales, esto se debe a que, al utilizar cuadriláteros se está aumentando el orden de aproximación en todo el dominio Ω al utilizar funciones cuadráticas. Sin embargo, aunque se está aumentando el orden, XFEM con cuadriláteros no va a evitar la imprecisión de los flujos sobre la interfase, sólo que para que XFEM falle, los valores de v₁ y v₂ deben ser de mucho mayor contraste (ver Figura 7).

Figura 7
Malla de cuadriláteros. Flujos numéricos obtenidos con FEM y XFEM, con valores v₁= 1 y v₂ = 5000.

Es de acotar que, para este ejemplo que se está tratando y para una disposición de interfase recta inclinada 10°, el método de XFEM+ (con malla de cuadriláteros) no resulta adecuado de implementar puesto que se estaría sobreestimando la interfase recta, ya que el Level Set sería aproximado con funciones cuadráticas, lo cual no tiene sentido. De aquí que, para un problema con una interfase curva, si sea adecuada la implementación del XFEM+.

4.2. Para una interfase curva

Sea Ω = (0,1) x (0,1), y Γ el lado derecho y superior de Ω. El problema fuerte a resolver es:

donde v₁= 1 y v₂ = 5000, y la interfase α está dada por un semicírculo de radio y centro (0,0).

En la Figura 8 se muestran las aproximaciones de los flujos para los tres métodos al igual que los flujos exactos.

Figura 8
Flujos numéricos obtenidos con los valores v₁= 1 y v₂ = 5000.

Se puede observar que el FEM falla completamente al aproximar los flujos cerca de la interfase, y XFEM mejora notablemente los flujos en la zona de transición. En este caso, como la interfase es una curva, se puede estimar de manera correcta usando el Level Set, con las funciones cuadráticas, observando que XFEM+ no desmejora los flujos globales y locales arrojados por XFEM.

5. Conclusiones

Se propuso un nuevo método de los elementos finitos para problemas de difusión con dos fases, haciendo énfasis en problemas de tipo elíptico, con condiciones de materiales distintas. Entonces,

• Dada la imposibilidad de FEM, de satisfacer la propiedad de discontinuidad del gradiente sobre la interfase, se modifica el método, enriqueciendo el espacio de soluciones, a fin de poder recrear esta característica de la solución exacta. Dicho método es llamado XFEM.

• Luego, para problemas con condiciones (conductividades) altamente distintas, XFEM es incapaz de recrear con precisión dicha continuidad en los alrededores de la interfase. Por ello, se recurre a imponer (de manera fuerte) esta condición sobre algunos puntos de la interfase. Esto es lo que se conoce como XFEM+.

• Para lograr esta imposición, se construyen unas ecuaciones algebraicas que serán añadidas (posteriormente) al sistema solución de XFEM; en ellas, se hace que el salto del flujo normal se anule en un conjunto finito de puntos sobre la interfase. De esta forma, se ha conseguido mejorar la calidad de los flujos numéricos en el entorno de la interfase.

En la etapa de implementación numérica de estos métodos (FEM, XFEM y XFEM+), se puede concluir que:

• Tanto en la utilización de mallas con elementos triangulares, o cuadriláteros: El FEM resulta un método bastante inadecuado para aproximar problemas de dos fases, esto debido a la incapacidad de la aproximación de tener gradiente discontinuo sobre la interfase.

• Tanto en la utilización de mallas con elementos triangulares, o cuadriláteros: Al enriquecer el espacio de soluciones para incorporar a la aproximación la propiedad de discontinuidad del gradiente sobre la interfase, el método XFEM resulta satisfactorio para problemas con conductividades distintas. Acotando que, al utilizar mallas con cuadriláteros, el método arroja mejores aproximaciones globales de los flujos, esto se debe a que se está aumentando el orden de la aproximación en todo el dominio, de funciones lineales (con los elementos triangulares) a funciones bilineales (con los elementos cuadriláteros).

• Cuando se tienen dos conductividades altamente distintas y una interfase recta inclinada: tanto en la utilización de mallas con elementos triangulares o cuadriláteros, se observa que XFEM empieza a perder precisión en cuanto a la aproximación de los flujos locales a la interfase, por ello la necesidad de mejorar dicha aproximación aplicando el método de XFEM+. En el caso de las mallas con elementos cuadriláteros, es de acotar que, no sería adecuado implementar el XFEM+, puesto que se estaría sobreestimando la interfase, ya que se estaría aproximando el Level Set (recuerde que el Level Set describe de forma discreta la interfase) con funciones cuadráticas (parábolas), lo cual no tiene sentido. De aquí que, para una interfase curva, XFEM+ resulte adecuado de implementar.

• En la utilización de mallas con elementos cuadriláteros, cuando se tienen dos conductividades altamente distintas y una interfase curva: al implementar XFEM y XFEM+ (para el ejemplo planteado) y comparar los flujos numéricos obtenidos con cada método, se pudo observar que XFEM es bastante preciso, arrojando una buena aproximación de los flujos globales y locales; en tanto que XFEM+ no desmejoró los flujos obtenidos por XFEM.

6. Recomendaciones

Es importante destacar que el enriquecimiento XFEM tiene un costo computacional, que viene dado por la extensión del espacio de soluciones. De aquí que, dependiendo de las necesidades del usuario que esté resolviendo el problema, resulte más adecuado utilizar mallas con elementos triangulares o cuadriláteros, sabiendo que con las mallas triangulares se puede tener un menor tiempo computacional, pero si se desea una mejor precisión de los flujos locales, va a ser inevitable tener que implementar el XFEM, e inclusive el XFEM+ (dependiendo de las características del problema). En tanto que la utilización de las mallas con cuadriláteros, acarrea un mayor tiempo computacional (mayor número de evaluaciones en las funciones de forma) pero se obtiene una mejor precisión de los flujos globales.

Referencias

[1] A. S. of Civil Engineers, Conference Papers2ndConference on Electronic Computation, Pittsburgh, Pa., September 8-9, 1960. Sponsored by the Committee on Electronic Computation, Structural Division, ASCE and the Pittsburgh Section, ASCE. American Society of Civil Engineers, 1960.

[2] R. Clough and Y. Rashid, "Finite element analysis of axisymmetric solids," Journal of the Engineering Mechanics Division, 1965,Vol.91, Issue 1, Pg.71-86, vol. 91, no. 1, pp. 71-86, 1965.

[3] H. Martin, Introduction to Matrix Methods of Structural Analysis. McGrawHill Book Company, New York, 1966.

[4] S. Pawsey and R. Clough, "Improved numerical integration of thick shell finite elements," International Journal for Numerical Methods in Engineering, vol. 3, no. 4, pp. 575-586, 1971.

[5] M. Turner, R. Clough ,H. Martin, and L. Topp, "Stiffness and deflection analysis of complex structures," Journal of Aerospace Science and Technologies, vol. 23, pp. 805-824, 1956.

[6] O. Zienkiewicz and R. Taylor, El Método de los Elementos Finitos: Formulación básica y problemas lineales, 4fh ed. McGraw-Hill Book Company, 1994, vol. 1.

[7] T. Belystchko and T. Black, "Elastic crack growth in finite elements with minimal remeshing," International Journal for Numerical Methods in Engineering, no. 45 (5), pp. 601-620, 1999.

[8] J. Chessa and T. Belystchko, "An extended finite element method for two-phase fluids," Transactions of the ASME, vol. 70, pp. 10-17, 2003.

[9] C. Daux, N. Moes, J. Dolbow, N. Sukumar, and T. Belystchko , "Arbitrary branched and intersecting cracks with the extended finite element method," International Journal for Numerical Methods in Engineering, no. 48, pp. 1741-1760, 2000.

[10] T. P. Fries and T. Belytschko, "The extended/generalized finite element method: An overview of the method and its applications," International Journal for Numerical Methods in Engineering, vol. 84, no. 3, pp. 253-304, 2010.

[11] N. Möes,J. Dolbow , and T. Belystchko , "A finite element method for crack growth without remeshing," International Journal for Numerical Methods in Engineering, no. 46, pp. 131-150, 1999.

[12] Q. Xiao and B. Karihaloo, "Improving the accuracy of xfem crack tip fields using higher order quadrature and statically admissible stress recovery," International Journal for Numerical Methods in Engineering, no. 66 (9), pp. 1378-1410, 2006.

[13] S. Zlotnik, P. Diez, M. Fernández, and J. Verges, "Numerical modelling of tectonic plates subduction using X-FEM," Computer Methods in Applied Mechanics and Engineering, no. 196, pp.4283-4293, 2007.

[14] J. Rodenas, O. González-Estrada, F. Fuenmayor, and F. Chinesta, "Enhanced error estimator based on a nearly equilibrated moving least squares recovery technique for fem and xfem," Computational Mechanics, no. 52 (2), pp. 321-344, 2013.

[15] O. González-Estrada , S. Natarajan, and C. Graciano, "Reconstrucción de tensiones para el método de elementos finitos con mallas poligonales," Revista UIS Ingenierías, no. 16 (1), pp. 23-34, 2017.

[16] F. Cordero and P. Diez , "XFEM+: una modificación de XFEM para mejorar la precisión de los flujos locales en problemas de difusión con conductividades muy distintas." Revista Internacional Métodos Numéricos para Cálculo y Diseño en Ingeniería, vol. 2, no. 26, pp. 121-133, 2010.

Notas

Notas