EL METODO DE LOS ELEMENTOS NATURALES BASADO EN FORMAS α EN ELASTICIDAD COMPRESIBLE Y CUASI INCOMPRESIBLE

E. Cueto    M. A. Martínez    M. Doblaré

División de Mecánica Estructural, Departamento de Ingeniería Mecánica. Centro Politécnico Superior. Universidad de Zaragoza. María de Luna, 3. E-50015 Zaragoza. España. E-mail: mdob1are@posta.unizar.es. Resumen

El método de los Elementos Naturales basado en formas α (MEN-α) forma parte del grupo de los llamados mιtodos sin malla. Consiste esencialmente en una variante del mιtodo de Galerkin en la que la interpolaciσn por vecinos naturales se usa para construir el conjunto de funciones de forma. El MEN-α se diferencia de otros mιtodos sin malla en su capacidad para la imposiciσn "exacta" de las condiciones de contorno esenciales o de Dirichlet y en el carαcter interpolante de sus funciones de forma. En este artículo se presenta una revisión general del método, de las condiciones que garantizan las características antes descritas y un conjunto de resultados en elasticidad compresible y cuasi-incompresible bidimensional que prueban lo anterior.

Palabras clave: Métodos sin malla, método de los elementos naturales, formas α, condiciones esenciales de contorno.

1. INTRODUCCION

Uno de los avances más importantes que se han producido en la última década dentro de la mecánica computacional ha sido la aparición y posterior desarrollo de una familia de métodos que se ha dado en llamar sin malla o de puntos. Estos métodos, basados por lo general en técnicas de Galerkin, pero también de colocación, se caracterizan por no necesitar discretizar el dominio mediante la generación de una malla en el sentido tradicional. La conectividad de los elementos se genera desde dentro del propio método, en un proceso transparente al usuario, de manera que se elimina gran parte del tiempo de preproceso del método de los elementos finitos (MEF). Las ventajas de estos métodos son patentes también en procesos con grandes desplazamientos y deformaciones, al no sufrir las limitaciones de los elementos finitos en cuanto a regularidad de la malla, ángulos interiores mínimos, etcétera.

A pesar de ello, de la gran cantidad de métodos de este tipo aparecidos (véase la revisión de Duarte [9] o la de Belytschko [2]) y del impacto que han tenido en la comunidad investigadora en pocos años, no puede decirse que supongan todavía una alternativa generalizada al MEF. Todos ellos sufren de problemas inherentes a su formulación, que han sido solucionados solo en parte. Uno de ellos, quizá el que más atención ha recibido es la imposición de condiciones de contorno [7] [12] [2] [15]. El hecho de que la mayoría de ellos posean funciones de forma aproximantes y no interpolantes hace que los coeficientes nodales en el sistema discreto de ecuaciones no sean los desplazamientos nodales.

En este artículo se presenta una revisión del último de estos métodos en ser aplicado a la elastostática: el método de los elementos naturales (MEN) [25, 3, 22] y la modificación del mismo que supone el incorporar el concepto de formas a para la restricción de la vecindad (MEN-α) [7, 6]. El MEN-α se diferencia de los demás métodos sin malla en que posee una función de forma estrictamente interpolante y en que es posible obtener una interpolación estricta en todo el contorno sin perjuicio de la exactitud del método. Se presenta a continuación una descripción del método, en su formulación estándar y basado en formas α. También se presenta su aplicación a la elastostática de materiales cuasi incompresibles, en la que muestra un comportamiento similar al elemento finito cuadrilátero bilineal. El artículo se completa con algunos ejemplos que demuestran este comportamiento

El método de los Elementos Naturales (MEN) es uno de los últimos métodos sin malla en ser aplicado a la elastostática y la elastodinámica lineales. Su planteamiento es similar a los demás métodos sin malla en el sentido de que se parte de una técnica de interpolación de datos irregularmente espaciados para construir la aproximación. Estas técnicas se escogen de forma que no exijan posiciones relativas de los nodos o no establezcan criterios acerca de ángulos mínimos en triangulaciones, como en el MEF, etc. En este caso, la técnica de interpolación se conoce como interpolación por vecinos naturales [27, 28]. Esta es usada para la construcción de una base de un subespacio finito del espacio total de búsqueda -de dimensión infinita- de la solución al problema planteado, mediante el método de Galerkin.

El método fue planteado por primera vez por Traversoni [25], que lo denominó método de Elementos Finitos de Vecindad Natural (Natural Neighbor Finite Elements). Braun y Sambridge [3, 4] lo aplicaron a diversos problemas dentro de la geofísica, campo dentro del cual la interpolación por vecinos naturales es especialmente útil, así como a la simulación de la interacción fluido-estructura. A estos autores se debe el nombre de Método de los Elementos Naturales (MEN). El análisis más profundo de su aplicación en la elastostática y en la elastodinámica lineales fue realizado por N. Sukumar [22, 20, 21, 5].

En esta sección se presenta un breve resumen de las características de este método en su aplicación a la resolución aproximada de ecuaciones en derivadas parciales (EDP).

El método de los Elementos Naturales (en adelante, el MEN) está basado, como se ha dicho, en la interpolación por vecinos naturales [18, 19]. Esta técnica hace uso de la triangulación de Delaunay y de su estructura dual, la teselación de Dirichlet o diagrama de Voronoi [26] para construir la interpolación. Estos conceptos son de excepcional importancia en la Geometría Computacional para la caracterización de una nube de puntos y la relación espacial entre éstos (O’Rourke [14]).

La triangulación de Delaunay es la única triangulación para una nube de puntos dada, que cumple el denominado criterio del circuncírculo vacío. Esto significa que dentro del círculo que pasa por los tres nodos de cada triángulo no hay ningún otro nodo de la nube (a excepción del caso degenerado en que cuatro nodos se sitúen sobre un cuadrado). En la figura 1 se puede observar la triangulación de Delaunay de un pequeño conjunto de puntos y su estructura dual, la teselación de Dirichlet.

Por simplicidad, el desarrollo que se construye a continuación se ciñe al espacio euclídeo en R² . La mayoría de los conceptos son extensibles a R³ o incluso a un espacio genérico Rⁿ , aunque frecuentemente algunas propiedades se pierden la generalización. En esos casos, se hará referencia explícita a esas particularidades.

En general, partimos de un conjunto de puntos N ={n ₁ , n₂ , ...., n_n } en R² . La teselación de Dirichlet (conocido también como diagrama de Voronoi o de Thiessen) de primer orden es la subdivisión del plano en regiones tales que todos los puntos de esa región están más cerca del nodo al que ésta está asociada que de ningún otro. Formalmente [13]

  (1)

donde T_I representa la celda de Voronoi asociada al punto I. d(·,·) representa la distancia euclídea entre dos puntos. Si definimos el bisector de dos puntos x_I , x_j como el hiperplano perpendicular a la línea éste divide al plano en dos semiplanos. Siendo h( x_I , x_j) el semiplano que contiene a x_I  y h( x_j , x_I ) el que contiene a x_j , se tiene una definición equivalente de la celda de Voronoi en la forma

       (2)

queda claro a partir de esta definición que la triangulación de Delaunay de un conjunto de puntos se construye sobre la envoltura convexa del conjunto N. Por envoltura convexa (convex hull en terminología inglesa, CH(N) o conv(N)), se entiende el más pequeño de los conjuntos convexos que contienen a N. La interpolación por vecinos naturales hace uso del diagrama de Voronoi de segundo orden. La celda de Voronoi de segundo orden T_IJ se define como el lugar geométrico de los puntos que tienen al nodo I como el más cercano y al J como el segundo en cercanía. Formalmente

   (3)

La definición de las coordenadas de vecino natural de un punto x se obtienen a partir del diagrama de Voronoi de primer orden. Si sobre éste, ya construido, se introduce x, las áreas de las celdas sufrirán una modificación por la presencia del nuevo punto. Las coordenadas de vecino natural proporcionan una cuantificación de esta variación y, por tanto, una medida de la proximidad relativa entre los nodos. Sea entonces k(x) una medida de Lebesgue (longitud, área o volumen si tratamos con R, R² o R³ , respectivamente) de la celda T _x Sea asimismo k _I (x) una medida equivalentemente definida para la celda de segundo orden T _xI . Entonces las coordenadas de vecino natural de x respecto al nodo I se definen como la relación entre el área de T_I  que ha sido transferida a T _x , respecto al área total de T _x . Esta definición se representa esquemáticamente en la figura 2 y se expresaría matemáticamente como:

     (4)

     (5)

La interpolación de la variable esencial del problema (en general, el desplazamiento) se realizaría entonces:

    (6)

u_I representa el desplazamiento del nodo I y ɸ _I (x) serían las funciones de forma asociadas a cada nodo.

Las propiedades del MEN derivan directamente de la interpolación por vecinos naturales. En esta sección se hace un repaso de ellas. En primer lugar, en contraste con la mayoría de las funciones de forma de los métodos sin malla, las coordenadas de vecino natural son funciones estrictamente interpolantes. Como tal entendemos que la superficie interpolada pasa exactamente por los valores nodales, en contraposición con el carácter aproximante de otras funciones, como las del método de mínimos cuadrados móviles (MMCM) [2]. Se cumple, por tanto, una propiedad fundamental como es

     (7)

o, dicho de otra forma, los parámetros nodales u_I son los desplazamientos nodales. La prescripción de desplazamientos nodales en el vector global de desplazamientos puede realizarse entonces por sustitución directa de estos valores. Esto supone de hecho una primera ventaja sobre otros métodos sin malla, pues evita tener que usar técnicas como interpoladores de Lagrange, etcétera. Sin embargo, debe tenerse en cuenta que esta propiedad no basta para una correcta imposición de condiciones de contorno esenciales, ya que la mera satisfacción de éstas en los nodos del contorno no asegura su cumplimiento en todo él.

Se puede observar también cómo, por construcción, (ec. 4), las funciones de forma del MEN constituyen una partición de la unidad, es decir,

    (8)

Por otro lado, es conocido que para que un método converja éste debe ser consistente y estable [2]. La consistencia del método está relacionada con el grado del polinomio que puede reproducir. Si la ecuación diferencial que rige el problema es de orden 2k, al método se le exige consistencia de orden k, es decir, que reproduzca un campo constante en la derivada de orden k de la variable esencial. En el caso de la elastostática, cuya ecuación diferencial es de grado 2, se debe poder reproducir un campo de desplazamientos lineal. Sibson [18] demostró que las coordenadas de vecino natural cumplen lo que se denomina propiedad de coordenada local, es decir:

     (9)

Esta propiedad, en conjunción con la de partición de la unidad (8), proporcionan la deseada demostración de la consistencia lineal.

En lo que se refiere a condiciones de continuidad y de derivabilidad, la función de forma del MEN es derivable, de clase C^∞ , en todos los puntos del dominio, excepto en los nodos, donde es C⁰ . Por otro lado, la aproximación que estas funciones proporcionan en determinado punto depende del número de vecinos naturales que tenga ese punto. En dominios unidimensionales, se ha demostrado [22] que la función de forma del MEN equivale a la de Elementos Finitos lineales.

En dos dimensiones, por el contrario, esto solo sucede cuando el punto en cuestión tiene tres vecinos naturales o cuatro situados en una malla regular. En el primer caso, las funciones de forma equivalen a las coordenadas baricéntricas, que como es sabido, constituyen las funciones de forma de los elementos triangulares lineales (constant strain triangles, CST), mientras que en el segundo equivalen a una aproximación bilineal. Si el número de vecinos naturales fuese cinco, o bien cuatro, pero situados irregularmente, las funciones de forma tienen una expresión racional cuártica

El MEN tiene también en lo relativo a imposición de condiciones de contorno esenciales unas características netamente diferentes del resto de métodos sin malla. La diferencia estriba en que las propiedades de aproximación del método son diferentes si tratamos con contornos cóncavos o convexos. Así, en contornos convexos, la función de forma del MEN es lineal a lo largo del contorno, mientras que en contornos no convexos no lo es. Por su interés, se reproduce a continuación la prueba que se puede encontrar en el trabajo de Sukumar [22, 20].

De acuerdo con (4), la función de forma se puede calcular como

    (10)

donde. En la misma figura se puede observar el hecho de que las celdas de Voronoi en esta parte del dominio no tienen límite superior y por tanto poseen un área no acotada. Por consiguiente, las áreas de la fórmula (10) se pueden expresar como

    (11)

donde las δ_I representan áreas finitas. Se tiene entonces que

     (12)

y resolviendo esos límites

Se puede ver entonces que debido al hecho de que las áreas asociadas a nodos en el contorno tienen un área infinita, la contribución de los nodos interiores en esta zona se hace nula. Siguiendo este mismo razonamiento, se puede observar cómo esta propiedad no es extensible a dominios no convexos. En ese caso, la contribución de nodos interiores se hace no despreciable frente a la de los nodos exteriores. En el trabajo de Sukumar [22, 20] se hace referencia a errores debidos a esta causa de un 2%.

En el trabajo de Sukumar [22, 20], la integración numérica de la forma débil del problema se realiza sobre los triángulos de la triangulación de Delaunay, tomando tres puntos por cada uno. Si bien ésta parece la forma más práctica de realizar esta integración, en ese mismo trabajo se habla de errores en el patch test del orden de 10^-3 . Dado que se ha demostrado también que la función de forma del método puede reproducir analíticamente un campo de desplazamientos constante y de variación lineal, este error debe provenir necesariamente de la integración numérica.

Este tipo de errores es común a la mayoría de los métodos sin malla, desde el momento en que se están integrando funciones sobre discretizaciones del dominio que no coinciden con el soporte de éstas. Este tipo de comportamiento, patente desde los primeros trabajos en métodos sin malla, permanece sin solución completa. Dolbow [8] trató este problema para el caso del EFGM, proponiendo una solución parcial solo válida en el caso de que se opere con soportes rectangulares, obtenidos por producto tensorial de soportes unidimensionales. Además, en la mayoría de los métodos sin malla se trabaja con funciones de forma no polinomiales, lo que complica aún más la determinación de la fórmula de integración numérica más apropiada para el conjunto de ecuaciones discretas del método de Galerkin

3. El metodo de los Elementos Naturales basado en formas α  

En el capítulo anterior se realizó un breve resumen de las características del método de los Elementos Naturales en su formulación estándar. En el presente se introducirá el concepto de forma α y se analizará su influencia en el método.

Se verá cómo la introducción de las formas α en el MEN asegura la interpolación lineal en el contorno de los dominios, con lo que se soluciona otro de los grandes problemas del método: el carácter no interpolante en contornos no convexos.

El problema podría resumirse entonces en los términos siguientes: ¿Contiene una nube de puntos información acerca de la forma del dominio sobre el que están distribuidos? Eh ojo humano es capaz de dar una respuesta afirmativa a esta pregunta pero ¿existe una forma de extraer esa información de manera automática, mediante el uso de un computador? La respuesta la dio Edelsbrunner [10, 11], estableciendo el concepto de forma α.

Un conjunto finito de puntos define un número también finito de "formas", entendiendo por tales un politopo (porción del espacio limitada  por un poliedro), que depende  del  nivel de detalle con que se quiera representar,  usualmente denominado α ∈ R. Lógicamente, ese nivel de detalle nunca podrá ser menor que la mínima distancia entre puntos. Las formas α parten de la triangulación de Delaunay^* como estructura fundamental de la nube de puntos. La triangulación de Delaunay está definida siempre sobre la envoltura convexa de la nube, es decir, sobre el más pequeño de los politopos convexos que pueden definirse englobando a todos los puntos. La forma α sería entonces una generalización del concepto de envoltura convexa a politopos no convexos. Dicho de otra forma, la envoltura convexa sería la forma α de menor nivel de detalle, esto es, α = ∞.

Sea N un conjunto finito de puntos en R³ . Sea también α ∈ R, α ∈ [0,∞) Se denomina α-esfera a la esfera abierta de radio α . Una α-esfera b estará vacía si  b ∩ N = Ø. Se denomina k-símplex, σ _T  , a la envoltura convexa de un subconjunto de N, T, tal que │T│= k + 1,con 0 ≤ k ≤ 3 (obsérvese que un 2-símplex es un triángulo y un 3-símplex es un tetraedro solo si se acepta lo que en el trabajo de Edelsbrunner se denomina la posición general de los puntos, es decir, que no hay cuatro puntos en un plano ni cinco en una esfera).

El concepto fundamental del desarrollo de la teoría de formas a viene de lo que Edelsbrunner denomina k-símplex α-expuesto. Se dice que un k-símplex  σ_T  está α-expuesto si existe una α-esfera b, con T= ∂b ∩ N, donde ∂b es la esfera o plano que limita a b. Un valor α ∈ R define entonces un conjunto F_{k ,α} de k-símphices α-expuestos, con 0 ≤ α ≤ 2. La forma α de N, S_α (N), se define entonces como el politopo cuyo contorno está formado por los triángulos de F _2,α , las aristas de F _{1, α} , y los vértices de F _{0, α} . Con esto se tendría bien definido el contorno, pero no qué parte de está dentro de S_α y cuál fuera.

Para cada triángulo α-expuesto existen dos α-esferas (no necesariamente vacνas) b₁ y b₂ tales que  b1 ≠ b2 , T Í ∂b₁ y  T Í ∂b₂ . Si las dos esferas están vacías entonces el triángulo no pertenece al contorno del interior de S_α. Si una de las esferas estuviera vacía y la otra no, entonces el triángulo limitaría el interior de S_α.

Estrictamente relacionado con lo anterior, podrían definirse los conjuntos F_k (0 ≤  k ≤  3) como los grupos de k-símphices ∂_T,│T│= k + 1, para los cuales existen esferas abiertas b tales que ∂b ∩ N = T. De esta forma, la triangulación de Delaunay, D del conjunto N se define como el complejo formado por los tetraedros en F₃, los triángulos en F₂, las aristas de F₁ y los puntos de F₀. Se cumple entonces la siguiente relación entre la triangulación de Delaunay y el contorno de S_α :

   (14)

Para cada punto p ∈ N se podría definir la celda de Voronoi V(p) como el conjunto de todos los puntos x Î R³ tales que la distancia euclídea entre  x y p es menor o igual que la distancia entre x y otro punto cualquiera de N. Esta definición, ya adelantada en el capítulo 2, pone de manifiesto que la celda de Voronoi es un poliedro convexo. Todas las celdas de Voronoi definen el diagrama de Voronoi de N, V. Cada celda del diagrama de Voronoi tridimensional sería una 3-celda de V. De la misma forma, una 2-celda sería la intersección de dos 3-celdas y una 1-celda sería la intersección de tres. Finalmente, una 0-celda sería la intersección de cuatro 3-celdas. La dualidad entre la triangulación de Delunay y el diagrama de Voronoi se establece entonces en la siguiente relación: sea T un subconjunto de N de tamaño │T│= k + 1, con 0≤ k ≤ 3 se define V_T¨=∩ _p∈N  V(p). Se tiene entonces que σ_T es un k-símplex de D si y solo si V_T es una (3 – k)-celda de V, 0 ≤  k ≤  3. Estrechamente relacionado con el concepto de forma a está el de α-complejo. De la definición de forma α se deduce que ésta es triangulada por un subconjunto de la triangulación de Delaunay, D . Un complejo simplicial C es una colección de k-símplices cerrados (0 ≤  k ≤  3) que satisface:

Considérese ahora un símplex σ_T, limitado por una esfera b_T de radio . Para 1 ≤ k ≤ 3 y 0 ≤ α ≤ ∞ se definen los conjuntos G_{k, α} como la uniσn de los k-sνmplices σ_T ∈ D para los cuales b_T está vacía y < α. G_{0, α} se define idιnticamente igual a N. Se define entonces el α-complejo de N, Cα (N) como el complejo simplicial cuyos k-símplices pertenecen a G_{k, α} o bien limitan (k + 1)-símplices de Cα.

(15)

la cual vuelve a poner de manifiesto que una forma α es triangulada por un subconjunto de la triangulación de Delaunay.

3.1.2 Formas α escaladas

En la mayoría de los trabajos de visualización científica se opera con conjuntos de puntos distribuidos de una manera más o menos regular, es decir, la separación nodal h es aproximadamente constante entre todos los puntos. Para este tipo de estructuras de puntos, las formas α han mostrado excelentes resultados. Sin embargo, ya desde los trabajos de Edelsbrunner [11] se puso de manifiesto que en el caso de conjuntos de puntos con distribución no uniforme la reproducción de la geometría no es apropiada. Se plantea entonces la posibilidad de operar con un cierto valor peso que afecte al valor del circunradio de la  α-esfera.

En la construcción de soluciones aproximadas a sistemas de ecuaciones en derivadas parciales se trata frecuentemente con conjuntos de puntos que se hacen más densos en las zonas en las que se prevé un gradiente de la variable esencial más alto, con el objeto de obtener una aproximación más exacta. El manejo de esta forma de datos plantea entonces la necesidad de reproducir esos dominios de la misma manera que una forma α lo haría con una distribución homogénea. En el trabajo de Cueto, Doblaré y Gracia [7] se plantea la posibilidad de tratar estos datos usando formas α escaladas según la densidad de puntos (Teichmann y Capps [23]). Aunque originariamente ideado para tratar casos complejos en distribuciones igualmente homogéneas, se ha demostrado cómo el uso de esta variante de las formas α ha dado excelentes resultados en la construcción de geometrías a partir de nubes de puntos.

En el trabajo de Teichmann y Capps [23] se plantean dos posibilidades complementarias: la construcción de formas  α anisótropas, en las cuales se modifica el tensor métrico localmente con vistas a que zonas próximas dejen de ser vecinas, y las formas  α escaladas, en las cuales se modifica el valor de  dependiendo de la densidad local de puntos alrededor del punto considerado. La primera forma está indicada en figuras con zonas muy próximas que no deben juntarse, intersticios, etcétera, mientras que la segunda parece mas apropiada para tratar con distribuciones no homogéneas.

Considérese de nuevo un conjunto N ={n ₁ , n₂ , ...., n_{I, ....} n_n }de puntos en R³. La densidad local de un punto n₁∈ N se define como:

   (16)

donde λ es un valor constante que indica una medida de la vecindad que rodea al punto y d(x, y) es la distancia euclídea. El valor de λ debe ser escogido por el usuario, de la misma forma que α  pero, como se verá mas adelante, representa solo una medida del orden de magnitud de la vecindad que rodea a n _I.

    (17)

entonces el valor del radio del circuncírculo se modifica, en función del valor de esa concentración nodal alrededor del punto en cuestión. El nuevo valor de α , α'  será entonces

    (18)

donde el valor β es otro nuevo valor que el usuario debe ajustar, al igual que α  y λ. µ representa la media de la densidad en todo el conjunto, como ya se ha dicho, y δ(σ_T) sería la densidad del triángulo bajo consideración. Teichmann y Capps presentan varias formas de medir esta densidad

La primera de estas formas de calcular la densidad ha mostrado excelentes resultados y parece ser además la más cóherente con aquello que se quiere obtener. Todos los ejemplos que se muestran en este trabajo se han calculado tomando la expresión (1) como referencia.

Los parámetros α, β y λ dependen del problema en consideraciσn, aunque de manera distinta. El valor de α depende del mínimo nivel de detalle con que se quiera representar la geometría y como es obvio, nunca podrá ser menor que la mínima separación nodal, h. El parámetro β depende del gradiente de la distribución de puntos y el valor de λ del tamaρo relativo de las zonas con distinto valor relativo de h o, equivalentemente, α. La figura 4 y la 5 representan dos ejemplos de formas α bi y tridimensionales, respectivamente.

En esta sección se mostrará cómo las formas α pueden modificar el esquema de interpolación por vecinos naturales. Se ha visto en la sección anterior que el hecho de que las celdas de Voronoi asociadas a nodos en contornos convexos no tengan área acotada proporciona el deseado caracter interpolante en el contorno.

Dado que esta circunstancia no se da en las celdas asociadas a nodos situados en contornos cóncavos, el valor del interpolante en esas zonas depende además del valor en nodos interiores al dominio. La interpolación por vecinos naturales pierde así el caracter interpolante en estos contornos.

La solución propuesta en Cueto, Doblaré y Gracia [7] parte de una conocida propiedad de los vecinos naturales. Si dos nodos son vecinos naturales, forman los extremos de un lado de un triángulo de Delaunay y comparten un segmento de una celda de Voronoi. De la misma forma, si se introduce un nodo x en la triangulación, sus vecinos serán aquellos que formen con él los nuevos triángulos. La forma α permite restringir esta vecindad a un determinado valor α. Se propone así una nueva definición para la celda de Voronoi:

        (19)

donde σ_T es el k-símplex formado por los nodos n_I, n_J y cualquiera del resto de los nodos de N, n_k. C_α (N) es el α-complejo asociado al dominio que se desea reproducir. Nótese que esta propuesta la vecindad ya no está asociada con la celda de Voronoi en el sentido clásico, dado que estas nuevas celdas tienen intersección no vacía fuera del dominio. De esta forma, un punto en el espacio puede estar asociado a más de una celda, contrariamente a lo que sucede en las celdas de Voronoi. Sí se mantiene la propiedad de que dos nodos que son vecinos forman un lado de un triángulo en C_α (N). Del mismo modo que en la forma estándar del método, las coordenadas de vecino natural de un punto respecto al nodo n_I se definen como la relación del área de la celda T_I transferida a T_x respecto al área de T_x. Esta forma de vecindad da lugar a lo que se denominó en la citada referencia [7] método de los Elementos Naturales basado en formas a ó MEN-α.

Para demostrar que la definición anterior da lugar a un interpolante lineal en el contorno, considérese un conjunto N de nodos distribuidos en una estructura regular (figura 6).

Supóngase un contorno cóncavo G _u  donde se han prescrito unas determinadas condiciones de contorno esenciales u = g(x). Siendo la separación entre nodos h, el valor del parámetro α apropiado sería  (distancia de la diagonal entre nodos). Un valor mayor de α podría dar lugar a una forma α en la que el nivel de detalle deseado no fuese suficiente, redondeando el ángulo que se quiere reproducir en el nodo B. Considérese también un punto   x ∈ G _u  en el que se quiere determinar el valor del interpolante. Para una imposición adecuada de las condiciones de contorno en este punto x debe conseguirse que éste tenga una celda asociada no acotada. De esta forma, la influencia de los puntos interiores sobre él se hace nula, como se vio en la sección 2.3. Para que esto ocurra, el nodo A, que en una triangulación de Delaunay sería vecino de x , y haría que su celda asociada tuviese un área finita, no debe estar asociado a x . La peor situación sería que el punto x estuviese cerca de B. En el límite, la α-esfera (círculo, en este caso) b que haría que A, B y x fuesen vecinos (es decir, que pertenecieran a C_α(N È x) tendría un radio α’ = h > α.

Esta solución es válida siempre que el ángulo formado por los segmentos   sea mayor de 90°. En el caso de la simulación de grietas, por ejemplo, este método debería ser sustituido por una triangulación conforme de Delaunay, pues resulta imposible aportar información sobre el frente de una grieta (con un radio de curvatura idealmente nulo) con nodos solamente, pues haría falta que la separación entre éstos fuese, en el límite, cero. Esto sería, de hecho, un segmento.

La función de forma resultante, asociada al nodo B puede verse en la figura 7. Puede verse la linealidad a lo largo de los segmentos  de , así como el valor unidad en B . 

En este trabajo, la función de forma se calcula tanto con arreglo al algoritmo de Watson [28] como al de Lasserre [16]. Ambos exigen almacenar los triángulos de Delaunay y su circuncentro y circunradio. Se construye una lista de los vecinos naturales para cada nodo, determinando simplemente si se cumple la ecuación

    (20)

donde v representa al circuncentro del triángulo en consideración, x el punto donde se desea conocer el valor de la función de forma y R el radio del circuncírculo. Si esta relación se diese, el punto x y los tres nodos que definen el triángulo cuyo circuncentro es v serían vecinos.

Un paso adicional de este algoritmo sería la eliminación de aquellos k-símplex que no pertenecen a Cα (N), pero como ya se ha dicho, este proceso es extremadamente sencillo, de complejidad O(n) .

Este algoritmo presenta una diferencia fundamental con respecto al original. Si bien no se ha modificado el cálculo de las funciones de forma, la definición de las celdas ya no proporciona una división unívoca del espacio Rⁿ. Como quiera que esto sólo sucede fuera del recinto de integración, no se modifica ninguna de las propiedades originales del método, sino el orden de interpolación en las distintas zonas del dominio. De hecho, el uso de un PSLG (conjunto de puntos y segmentos de recta que deben ser mantenidos en la triangulación final) para definir el contorno del dominio y la posterior eliminación de los triángulos que cayesen fuera de él también rompe esta dualidad, dado que se trabaja en todo momento con el diagrama de Voronoi a través de su estructura dual, precisamente la triangulación de Delaunay. La eliminación, por tanto, de algunos triángulos rompe inevitablemente la dualidad entre triangulación y teselación de Voronoi. Este aspecto no es puesto de manifiesto en el trabajo de Sukumar [22]. La mayoría de sus ejemplos parecen tener, en efecto, función de forma lineal en el contorno. El comportamiento del interpolante en el contorno se tratará más en profundidad en los ejemplos numéricos que completan esta sección.

A pesar de ello, el uso de un PSLG y una triangulación conforme de Delaunay no asegura el comportamiento lineal de la función de forma en el contorno, como se probará en los apartados siguientes, y presenta ciertas dificultades, como puede ser la extensión a tres dimensiones, donde la construcción de una triangulación conforme de Delaunay es un tema de investigación abierto en estos momentos.

Relacionado con estos aspectos, desde los primeros trabajos en métodos sin malla [17] se puso de manifiesto que el cálculo de la distancia entre nodos (equivalente a la noción de vecindad natural en el MEN) no debe realizarse a través de una porción del espacio Rⁿ que esté fuera del dominio del problema Ω. De esta forma, si el vector que une dos nodos atraviesa el contorno del dominio, ιste debe descomponerse en dos o más tramos que caigan enteramente en el interior del dominio, dando lugar a funciones de forma discontinuas.

Tampoco en el trabajo de Sukumar [20, 22] se ha tenido en cuenta este problema, limitándose solo al análisis de la linealidad del interpolante en el contorno. La importancia que este factor puede tener en determinados casos se pone de manifiesto en el ejemplo 5.1. Nuevamente, la definición del interpolante sobre una forma α asegura el establecimiento de la vecindad de una manera apropiada, impidiendo que los nodos estén relacionados a través de porciones del espacio no comprendidas en el dominio.

4. El MEN EN LA ELASTICIDAD INCOMPRENSIBLES

4.1 Formulación débil del problema

El trabajo que se presenta en este artículo se ciñe al problema de la elastostática lineal. Las ecuaciones que gobiernan este problema se pueden formular como sigue: sea un dominio W Ì Rⁿ, con n = 1, 2, 3, de contorno   Entonces  

         (21)

            (22)

Los términos representan la parte simétrica del operador gradiente actuando sobre el campo de desplazamientos, esto es, el tensor de pequeñas deformaciones y el tensor C es el tensor de comportamiento del material.

Las condiciones de contorno del problema se establecen en su forma mixta como

  (23)

   (24)

debiéndose cumplir que y n representa el vector normal al contorno. Sea entonces . La forma débil del problema de valores en el contorno dado por las ecuaciones (21), (23) y (24) quedaría expresada como

  (25)

La forma discreta del problema (método de Petrov-Galerkin) se tendría al escoger unas aproximaciones a u y v de la forma:

    (26)

donde v_1i son coeficientes arbitrarios, con tal de que v_i^h se anule en . Se obtiene así un sistema discreto de ecuaciones

    (28)

El método de Bubnov-Galerkin (o, simplemente, Galerkin) se obtiene al hacer coincidir las funciones . En concreto, en el Método de los Elementos Naturales, estas funciones están constituidas por las coordenadas de vecino natural, representadas en la ec. (4), del punto en cuestion.

4.2 Formulación mixta

El fenómeno del bloqueo de la solución numérica del problema elástico cuando el coeficiente de Poisson se aproxima a 0.5, esto es, cuando las características del material se aproximan a la incompresibilidad, es bien conocido. Una de las apuestas más frecuentes para la solución de este problema es el uso de principios variacionales mixtos, es decir, métodos en los cuales el desplazamiento y la presión (entre otras muchas posibilidades) son tomados como variables independientes -principio variacional de Hellinger-Reissner-.

La forma fuerte del problema de la elastostática vendría entonces dada por:

donde λ = vE / (1 + v)(1 – 2v) toma valor ∞ si v = 0.5, como puede verse. Para dar lugar a la forma débil se definen, como de costumbre, los espacios de aproximación en la forma:

Tras aplicar el teorema de la divergencia se llega a

donde μ = E/2(1 + v) es el segundo coeficiente de Lamé. Es bien sabido, asimismo, que para que el conjunto ( u , p) seasolución única del problema, la forma bilineal b( u , p) debe cumplir la condición LBB [1], también llamada condición inf-sup:

En el MEN las variables esenciales del problema se aproximan de la siguiente manera:

llegando finalmente al sistema de ecuaciones discreto

donde

siendo:

Por lo que respecta al tipo de aproximación escogida, es decir, a las funciones , se debe destacar que, como se comprobará en los ejemplos numéricos que acompañan a este artículo, tanto la aproximación estándar del MEN como la aproximación de tipo C⁰ – C⁰ –es decir, dan lugar a bloqueo de las soluciones. Debe tomarse, al igual que sucede en Elementos Finitos de tipo cuadrilátero bilineal, una aproximación de orden menor para las presiones que para los desplazamientos (de tipo C⁰– C^-1 , por ejemplo). Para lograr la aproximación de tipo C^-1 existen dos posibilidades: la primera sería tomar ψ(x) 1/n, donde n representa el número de vecinos naturales del punto considerado; la otra sería la denominada interpolación por vecino más cercano, es decir, asignar a la variable esencial (en este caso, a la presión) el valor nodal en toda la celda de Voronoi. El estudio de los resultados obtenidos con estas aproximaciones se pospone a la sección 5. 

RESULTADOS NUMERICOS

5.1 Placa con un agujero elíptico

En este apartado se probará cómo el MEN y el MEN-α no son equivalentes en dominios no convexos y cómo solamente la construcción del interpolante sobre una forma a garantiza la linealidad del mismo. Para ello, se considera una placa bidimensional con un agujero elíptico. La geometría del problema se muestra en la figura 8. La placa está sometida a una tensión unitaria aplicada en el infinito.

Figura 8. Geometría del problema de la placa con agujero elíptico

La placa se modelizó, por un lado, con un conjunto de solo 60 nodos y, por otro, con ese mismo conjunto de nodos, incluyendo la definición explícita del contorno del dominio (MEN estándar). Como propiedades del material se tomaron: módulo de elasticidad E = 1.0 y coeficiente de Poisson v = 0.25. Se simuló el problema bajo un estado de tensión plana. La solución teórica de este problema se puede encontrar, como en el caso anterior, en el libro de Timoshenko y Goodier [24]. La expresión de esta solución se suele dar en el campo complejo:

   (49)

     (50)

    (51)

donde  son las componentes del tensor de tensiones en coordenadas elípticas, u = u + iv es el campo de desplazamientos en coordenadas cartesianas y, k es, de nuevo, la constante de Kolosov. A y B son constantes, definidas como A = σ/2, B = – (1/2)a cosh  2x₀, siendo s la tensión aplicada en el infinito y x₀ el radio del agujero en coordenadas elípticas. Finalmente, representa el vector de posición en coordenadas elípticas.

Los dos modelos construidos (el basado en formas α y el construido a partir de un PSLG y una triangulación conforme) son idénticos, diferenciándose solamente en la posición de un nodo. La placa simulada con una forma α se reproduce en la figura 10, mientras que la nube de puntos original se muestra en la figura 9. La placa simulada con un PSLG y triangulación conforme se presenta a continuación (figura 11). En este último caso, se representan también los triángulos que completan la envoltura convexa del conjunto de nodos para mostrar la conformidad de la triangulación, aunque, como es lógico, se omitieron en los cálculos finales.

Figura 10. Forma a para el problema de la placa con agujero elíptico

Figura 11. Triangulación conforme de Delaunay para el problema de la placa con agujero elíptico

La razón de esta diferencia está en el deseo de poner de manifiesto que es posible (y, en general, así será) que el interpolante construido a partir del MEN estándar no sea lineal en el contorno, mientras que en el MEN-α la reproducciσn del contorno de una forma apropiada asegura esta linealidad. La importancia que este hecho puede tener numιricamente también se pone de manifiesto aquí. Así, en el caso de la placa simulada con el MEN estándar un nodo cercano al contorno se ha desplazado hacia éste, buscando que se convierta en vecino natural de otros en el mismo contorno y que se rompa la linealidad.

Como en el caso de la placa con agujero circular, se consideró solamente una cuarta parte de una porción finita de la placa, aplicándose las condiciones de contorno apropiadas, derivadas de las expresiones 49 y 50. Las normas de error se muestran en la tabla 1.

Este ejemplo muestra claramente la diferencia entre las dos variantes del método. En el MEN- a el dominio del modelo es construido exclusivamente a partir de nodos, sin definición explícita del contorno del dominio. Esta es una característica que hasta ahora solo poseían los métodos basados en colocación, razón por la cual algunos autores los denominan de puntos frente a los métodos basados en Galerkin, que necesitan en efecto de esa definición explícita del contorno y de una discretización del dominio para la integración numérica (aunque esto no pueda considerarse en sentido estricto una malla, pues los requerimientos de ésta en cuanto a regularidad, ángulos mínimos, etc. son mucho mayores). Si la posición de los nodos es escogida de forma apropiada —lo cual quiere decir solamente que la separación nodal, h, debe ser estrictamente menor que el mínimo radio de curvatura que se desee reproducir en el modelo— la linealidad del interpolante en el contorno está asegurada y la imposición de condiciones de contorno esenciales puede hacerse de forma conveniente.

TABLA 1. Resultados para el problema de una placa con agujero elíptico

Por otro lado, si se utiliza el MEN estándar mediante la definición explícita del contorno del dominio la linealidad del interpolante no está asegurada y dependerá de la posición de los puntos de Steiner (puntos no controlados por el usuario que se introducen para mantener el caracter de Delaunay de la triangulación). Controlar la posición de los puntos de Steiner no es una tarea sencilla especialmente en casos tridimensionales, donde los algoritmos de triangulación se vuelven considerablemente más complejos.

5.2 Simulación bidimensional de la flexión en material incompresible

Se presenta el ejemplo de la flexión de una viga bidimensional en voladizo sometida a carga en su extremo libre. Para ello, se considera una discretización compuesta por una nube de 85 nodos, que, junto con las dimensiones del problema, puede verse en la figura 12. El material se toma como elástico, de módulo de Young unitario y coeficiente de Poisson variable, desde 0.3 hasta 0.4999999. Los resultados, para los distintos valores del coeficiente de Poisson, se muestran en la tabla 2.

Figura 12. Nube de puntos para ala simulación bidimensional de la viga en voladizo

TABLA 2. Resultados para el problema de viga bidimensional en flexión

Los resultados representan en este caso la flecha (en tanto por ciento respecto a la solución teórica) en el extremo libre, observándose el conocido bloqueo de la solución de Elementos Finitos, a medida que el coeficiente de Poisson se acerca a 0.5. Este bloqueo también se hace presente en la simulación por Elementos Naturales con aproximación estándar en desplazamientos. Sin embargo, la aproximación mixta C⁰ – C⁰ no muestra este bloqueo, como puede observarse. Tanto los Elementos Finitos como los Elementos Naturales presentan un comportamiento similar, con ausencia de bloqueo, a pesar de no verificar la condición LBB ninguna de las dos formulaciones. Los resultados más estables vienen dados en ambos casos (MEF y MEN) por la aproximación discontinua en presiones C⁰ – C^-1.

 5.3 Flexión tridimensional de material incompresible

Al igual que en el apartado anterior, se estudia el comportamiento del MEN en la simulación del comportamiento a flexión de materiales incompresibles, pero esta vez en tres dimensiones. El estudio del comportamiento tridimensional del MEN ha sido llevado a cabo por Cueto, Calvo y Doblaré [6]. En esencia, se concluye en esa referencia que el comportamiento tridimensional del MEN es análogo al comportamiento bidimensional cuando se trata con materiales compresibles. En este trabajo se presenta el estudio correspondiente al material incompresible.

Para la simulación se consideró de nuevo un material con módulo de Young 1.0 y coeficiente de Poisson variable, desde 0.3 hasta 0.4999999. Los resultados para esta simulación, nuevamente en relación a la flecha teórica en el extremo de la viga, que se toma como el 100%, se presentan en la tabla 3.

Nuevamente, se observa un comportamiento similar al del apartado anterior y, al igual que en elastostástica compresible, el comportamiento tridimensional es una generalización del comportamiento bidimensional. En este caso, la aproximación discontinua en presiones ofrece resultados ligeramente mejores que la correspondiente de tipo C⁰ – C⁰

TABLA 3. Resultados para el problema de viga tridimensional en flexión

6. Conclusiones

En este artículo se ha presentado un estudio del comportamiento del método de los Elementos Naturales (MEN) a la simulación del comportamiento de materiales elásticos incompresibles en dos y tres dimensiones. Al igual que en la elastostática bidimensional [22] [7] y tridimensional [6], el MEN ofrece un importante potencial de uso. Se ha demostrado que el comportamiento en la simulación de materiales incompresibles es altamente similar al proporcionado por elementos finitos cuadriláteros con interpolación bilineal en desplazamientos y constante en presiones, evitándose el bloqueo de la solución en situaciones muy próximas a la incompresibilidad.

Sin embargo, la aproximación mixta discontinua en presiones para el MEN no verifica la condición LBB, con lo que se hace necesario un estudio posterior en la búsqueda de posibles soluciones a dicho comportamiento.

Con todo, el comportamiento antes descrito hace suponer que el uso del MEN en la simulación de fluidos incompresibles, por ejemplo, muestre un gran potencial de uso, unido a las ventajas derivadas de su carácter de método sin malla. Otros aspectos del MEN, como la integración numérica están aún pendientes de una revisión en profundidad.

7. REFERENCIAS

1. I. Babuška. Error bounds for finite element method. Numer. Math., 16:322-333, 1971.2. 

2. T. Belytschko, Y. Krongauz, D. Organ, M. Fleming, and P. Krysl. Meshless methods: An overview and recent developments. Computer Methods in Applied Mechanics and Engineering, 139:3-47, 19.        [ Links ]

3. J. Braun and M. Sambridge. A numerical method for solving partial differential equations on highly irregular evolving grids. Nature, 376:655-660, 1995.        [ Links ]

4. J. Braun, M. Sambridge, and H. McQueen. Geophysical parametrization and interpolation of irregular data using natural neighbours. Geophysical Journal International, 122:837-857, 1995.        [ Links ]

5. D. Bueche, N. Sukurnar, and B. Moran. Dispersive Properties of the Natural Element Method. Computational Mechanics, In Press, 1999.        [ Links ]

6. E. Cueto, B. Calvo, and M. Doblaré. Modeling piece-wise homogeneous domains using the a-shape based Natural Element Method. International Journal for Numerical Methods in Engineering, Submitted, 2001.         [ Links ]7. E. Cueto, M. Doblaré, and L. Gracia. Imposing essential boundary conditions in the Natural Element Method by means of density-scaled a-shapes. International Journal for Numerical Methods in Engineering, 49-4:519-546, 2000.        [ Links ]

8. J. Dolbow and T. Belytschko. Numerical Integration of the Galerkin Weak Forrn in Meshfree Methods. Computational Mechanics, 23:219-230, 1999.        [ Links ]

9. C. A. M. Duarte. A Review of Some Meshless Methods to Solve Partial Differential Equations. Technical Report 95-06, TICAM, University of Texas at Austin, 1995.        [ Links ]

10. H. Edelsbrunner, D. G. Kirkpatrick, and R. Seidel. On the shape of a set of points in the plane. IEEE Transactions on Information Theory, IT-29(4):551-559, 1983.        [ Links ]

11. H. Edelsbrunner and E. Mücke. Three dimensional alpha shapes. ACM Transactions on Graphics, 13:43-72, 1994.        [ Links ]

12. J. Gosz and W. K. Liu. Adrnissible approximations for essential boundary conditions in the reproducing kernel particle method. International Journal for Numerical Methods in Engineering, 38:22-22, 1996.        [ Links ]

13. P. J. Green and R. Sibson. Computing Dirichlet tesselations in the plane. The Computer Journal, 21:168-173, 1978.        [ Links ]

14. J. O’Rourke. Computational Geometry in C. Cambridge University Press, 1994.         [ Links ]15. Y. Krongauz. Application of Meshless Methods to Solid Mechanics. PhD thesis, Northwestern University, Evanston, IL, 1996.        [ Links ]

16. J. B. Lasserre. An analytical expression and an algorithm for the volurne of a convex polyhedron in R n . Journal of Optimization Theory and Applications, 39(3):363-377, 1983.        [ Links ]

17. D. Organ, M. Fleming, T. Terry, and T. Belytschko. Continuous meshless approximations for nonconvex bodies by difraction and transparency. Computational Mechanics, 18:1-11, 1996.        [ Links ]

18. R. Sibson. A Vector Identity for the Dirichlet Tesselation. Mathematical Proceedings of the Cambridge Philosophical Society, 87:151-155, 1980.        [ Links ]

19. R. Sibson. A brief description of natural neighbour interpolation. In Interpreting Multivariate Data. V. Barnett (Editor), pages 21-36. John Wiley, 1981.        [ Links ]

20. N. Sukumar. The Natural Element Method in Solid Mechanics. PhD thesis, Northwestern University, Evanston, Illinois, 1998.        [ Links ]

21. N. Sukumar and B. Moran. C 1 Natural Neighbour Interpolant for Partial Differential Equations. Numerical Methods fon Partial Differential Equations, 15(4):417-447, 1999.        [ Links ]

22. N. Sukurnar, B. Moran, and T. Belytschko. The Natural Element Method in Solid Mechanics. International Jounnal for Numerical Methods in Engineening, 43(5):839-887, 1998.        [ Links ]

23. M. Teichrnann and M. Capps. Surface reconstruction with anisotropic density-scaled alpha shapes. In Proceedings of the 1998 IEEE Visualization Conference, 1998.        [ Links ]

24. S. Timoshenko and J.Ñ. Goodier. Teoría de la Elasticidad. Editorial Urmo, 1972.L. Traversoni. Natural neighbour finite elements. In Intl. Conference on Hydraulic Engineering Software. Hydrosoft Proceedings, pages 291-297. Computational Mechanics publications, 1994.        [ Links ]

25. G. M. Voronoi. Nouvelles Applications des Pararnétres Continus á la Théorie des Formes Quadratiques. Deuxiéme Memoire: Recherches sur les parallélloédres Primitifs. J. Reine Angew. Math., 134:198-287, 1908.        [ Links ]

26. D. Watson. Computing the n-dimensional Delaunay Tessellation with Application to Voronoi Polytopes. The Computen Journal, 24(2):162-172, 1981.        [ Links ]

27. D. Watson. Nngridr. An Implementation of Natural Neighbor Interpolation. Published by the author, 1994.        [ Links ]

PIE DE PAGUINA

*En la mayoria de la literatura existen,tanto sobre trangulaciones de Delauanay como sobre formas a, se habla por simplicidad,de la triangumación de la nube de puntos,incluso si se trata de una nube puntos en  R³ . La teoría de formas alfa es aplicable tanto a dos como tres dimensiones por lo que el desarrollo que sigue debe entenderse generalizado al caso tridimensional33