A mi esposa y compañera Ivette, a mis lindos hijos Flavia y Gianluca, por todo el cariño y apoyo entregado durante la realización de este trabajo

Agradecimientos

Deseo expresar mi gratitud y agradecimiento en primer lugar a los profesores Dr. Eugenio Oñate Ibáñez de Navarra y Dr. Juan Miquel Canet, por toda la ayuda, orientación y apropiados consejos entregados durante la dirección de esta tesis. La palabra y motivación oportuna que me han entregado, posibilitó mi comprensión de los métodos sin malla en este apasionante mundo de la investigación.

Mi total gratitud también para mi esposa e hijos, por la paciencia y comprensión demostrada, también por el constante apoyo y en especial, el cariño y todas las alegrías que me han brindado. En estos instantes especiales de mi vida, deseo expresar y agradecer a mi padre Aldo Perazzo y mi madre María Maggi, el haberme enseñado el camino correcto para llegar a tan importante meta.

Mi más sincero agradecimiento para mi amigo Dr. Carlos Sacco, por su inusual característica de escuchar y atender siempre a mis dudas, sin duda echaré de menos todos los análisis y conversaciones que mantuvimos en el bar de la escuela, en torno al tema de los métodos sin malla.

Quiero agradecer también a todos mis amigos y compañeros que me han alentado y acompañado durante este largo camino, a los doctores Alex Hanganu y Juan Carlos Cante, y finalmente a Quino, Loli, Carlos, Dudiño y Fernando, a quienes les deseo toda la suerte del mundo.

Por último quiero tener una palabra de agradecimiento para con mis colegas del departamento de Ingeniería Mecánica de la UTFSM de Chile, a todo el personal del CIMNE y del Departamento de Resistencia de Materiales de la UPC en Barcelona, por todo el apoyo y confianza que me han brindado.

1 Introducción

El método de elementos finitos (MEF) ha sido durante los últimos 25 años, y sigue siendo, la herramienta numérica más utilizada para simular diversos problemas en mecánica computacional. El buen prestigio que ostenta el método se debe en gran medida a la fiabilidad de sus resultados, cuando este se utiliza adecuadamente, y la variada gama de problemas que permite simular. Sin embargo, en la última década se ha producido un avance importante en el conocimiento y posibilidades de aplicación de un nuevo tipo de método numérico denominado 'sin malla', en inglés meshless, gridless, element-free. Este método, que presenta ciertas particularidades desde el punto de vista de su formulación, es hoy en día objeto de numerosas investigaciones y lentamente se ha comenzado a posicionar como una alternativa de solución para problemas en donde los métodos tradicionales, como el de elementos finitos y volúmens finitos, presentan ciertas desventajas.

A pesar del auge experimentado en el último tiempo por esta técnica, su utilización y primeras aplicaciones corresponden a la década de los setenta en el campo de la astrofísica [59], principalmente para modelar determinados fenómenos como fisión de estrellas en rotación, donde intervienen masas de fluído en movimiento sin la presencia de contornos. Durante los siguientes diez años el método, bautizado con el acrónimo de SPH (smoothed-particle-hydrodynamics) pasó inadvertido hasta que su versatilidad y posible utilización en otros campos fué redescubierta [64], dando lugar al desarrollo de diversas variantes que han dado origen a otros tantos nuevos métodos sin malla.

Entre éstos se destacan aquellas técnicas en las que se construye la aproximación local, utilizando un polinomio interpolante mediante mínimos cuadrados, hoy en día, se reconoce al método DEM (diffuse-element.method) [69] como precursor en utilizar este tipo de técnica en conjunto con una formulación de Galekin. Sin embargo, al igual que en otros métodos sin malla, en el DEM es necesario recurrir a una malla de fondo o auxiliar para resolver la cuadratura numérica en su forma débil. Esto produce cierta confusión respecto de la idea de método sin malla planteada en un comienzo y, aunque el uso de una malla de fondo sea para propósitos de cuadratura numérica y no para construir la aproximación, se aparta de la filosofía o motivación original de este tipo de técnica.

Una forma alternativa de eliminar el uso de esta malla de fondo es utilizar un procedimiento de Colocación Puntual, el cual también se utiliza en el método SPH, empero, el método resultante puede presentar ciertas restricciones debido a la sensibilidad en la ubicación de estos puntos de colocación. Buscar una forma de atenuar esta sensibilidad o contar con una estrategia para corregirla, permitiría aprovechar en toda su potencialidad el procedimiento de colocación en el contexto de los métodos sin malla.

1.1 Motivación

En la utilización de todo método sin malla se busca que el modelo discreto sea completamente descrito por nodos. Con esto se consigue que los datos necesarios para describir el dominio de la solución sean simplemente las coordenadas de estos nodos y cierta información del contorno. La ausencia de una malla de elementos o elementos y sus conectividades, reporta ciertas ventajas respecto de los métodos tradicionales:

  • La sencillez para discretizar el dominio de la solución en base a un conjunto de puntos distribuidos arbitrariamente, esta información por ejemplo se puede obtener a partir de cualquier sistema CAD
  • No existe una conectividad fija entre los nodos que conforman el subdominio de interpolación y que permiten obtener la aproximación en forma local
  • Menor tiempo empleado en la preparación de la información necesaria para el cálculo. Buena parte del tiempo total de una modelación compleja se invierte en realizar el modelo y la discretización por elementos o mallado, si bien los generadores de malla hoy en día son rápidos y eficientes en geometrías relativemente sencillas, lograr una malla de calidad en problemas tridimensionales sigue siendo una tarea ardua y no exenta de problemas. Por ejemplo, el proceso de generación de una malla en 3D puede tardar más tiempo que el tiempo de cálculo involucrado en la solución del sistema de ecuaciones discreto del modelo [84].
  • Facilidad para realizar un procedimiento de solución adaptable, puesto que sólo es necesario agregar más nodos en las zonas de interés sin preocupase por la cercanía con los demás. En el MEF la cercanía entre dos nodos produce un fuerte gradiente en la función de forma, lo que origina errores numéricos en la solución [71].
  • Facilidad para modelar problemas que involucren grandes deformaciones de la geometría o superfiicies en movimiento. Grandes distorsiones en los elementos de una malla hacen necesario realizar un procedimiento de remallado antes de efectuar el cálculo, esto conlleva un encarecimiento del coste computacional en aquellos problemas como conformado de metales o propagación de una grieta, que se caracterizan por un constante cambio en la geometría del dominio de análisis.
  • Pueden proporcionar soluciones cuyas derivadas sean continuas. Es fácil construir funciones del tipo , con lo cual el campo de deformaciones y tensiones en el dominio es continuo, haciendo innecesario el uso de técnicas de alisado que pueden resultar complicadas y costosas de implementar.

Para dimensionar el alcance real que pueden tener estas ventajas en un método sin malla, surge la necesidad de estudiar detenidamente su formulación matemática y evaluar las distintas alternativas que presenta su implementación para la resolución de las ecuaciones de la elasticidad lineal de sólidos. Esta ha sido la principal motivación para el desarrollo de la presente tesis, que nace dentro del las líneas de investigación en métodos sin malla que en forma conjunta llevan a cabo el Departamento de Resistencia de Materiales y Estructuras en la Ingeniería y el Centro Internacional de Métodos Numéricos en la Universidad Politécnica de Cataluña, En particular, se enfocará el estudio en el llamado 'Método de Puntos Finitos (MPF)' que ha sido la respuesta dada por estos centros de investigación a la iniciativa de formular un nuevo método sin malla.

Es pertinente mencionar en este apartado una breve aclaración respecto del nombre de Puntos Finitos, con el cual se ha denominado al método numérico objeto de investigación en esta tesis. El nombre se debe a las características particulares del método, es decir: una técnica numérica totalmente libre de malla en la cual el dominio se discretiza mediante un número finito de puntos.

1.2 Objetivos de la tesis

El objetivo principal de la tesis es obtener y desarrollar la formulación del Método de Puntos Finitos, como método totalmente libre de malla, para la resolución de las ecuaciones de la elasticidad en sólidos. Para cumplir este objetivo se pretende determinar y estudiar los aspectos relevantes de la fundamentación matemática del método, entre otros, la aproximación mediante mínimos cuadrados móviles, la función de ponderación y la técnica de colocación puntual. Para un completo análisis de las capacidades del método, se implementará un código que permita resolver diversos problemas, de índole mas bien académico, principalmente en 1D, 2D y en forma complementaria en 3D.

Como objetivos secundarios merecen ser destacados además, los siguientes:

  • estudiar y analizar las principales propiedades de la aproximación, entre otras la consistencia y convergencia
  • desarrollar una estrategia para seleccionar adecuadamente los puntos que conformarán la nube o subdominio de interpolación
  • analizar la sensibilidad del método respecto de la adecuada colocación de los puntos en el dominio, tanto para discretizaciones regulares como irregulares de nodos
  • evaluar la correcta implementación de las condiciones de contorno
  • identificar aquellos aspectos relevantes que afecten y comprometan la exactitud de la aproximación, buscando plantear alternativas de solución para estos casos

Se pretende finalmente extender el desarrollo en la formulación del MPF para abarcar, además del análisis estático, algunos problemas de la dinámica lineal de sólidos.

1.3 Contenidos de la tesis

Para cumplir con los objetivos de este trabajo, la tesis se ha estructurado en 8 capítulos, el primero de ellos describe la motivación y el marco teórico en que esta se desarrollará. En el Capítulo 2 se entregan los principales fundamentos teóricos de los métodos sin malla bajo la perspectiva de las tres caracteristicas básicas que los distinguen como son, el tipo de aproximación, la función de ponderación y la forma de resolver el sistema discreto de ecuaciones. Bajo estos mismos conceptos, en el Capítulo 3 se presentan los aspectos matemáticos necesarios para comprender el Método de Puntos Finitos como método sin malla. Es en el Capítulo 4 donde se analiza en detalle la consistencia y convergencia del MPF, a través del desarrollo de diversos ejemplos de resolución de las ecuaciones de la elasticidad en sólidos, los resultados de este capítulo son fundamentales para comprender la importancia de plantear la estrategia que se presenta posteriormente. En el Capítulo 5 se formula una nueva metodología para mejorar la solución numérica de la aproximación por el MPF, principalmente en las zonas del contorno. Con los ejemplos de la elasticidad lineal de sólidos, estáticos y dinámicos, implementados en los Capítulos 6 y 7, se pretende comprobar la validez de la modificación propuesta. Se completa el estudio y análisis de la metodología sin malla formulada, con las conclusiones y las futuras líneas de investigación presentes en el Capítulo 8.

2 Revisión de los fundamentos teóricos de los métodos sin malla

Es conveniente especificar o definir globalmente, a manera de introducción y para aclarar conceptos, las características que posee una técnica numérica para que pueda ser interpretada como sin malla o libre de malla.

Duarte [25] entrega en su trabajo una primera propuesta para definir este tipo de técnica: un método se considera sin malla si las ecuaciones básicas que gobiernan el modelo discreto del problema de contorno no dependen de la disponibilidad de una malla bien definida. Sin embargo, como podrá constatarse en este capítulo, algunos métodos son aceptados como sin malla a pesar de tener una débil dependencia respecto de una malla utilizada para calcular la cuadratura numérica de las ecuaciones integrales que gobiernan el modelo discreto.

Con posterioridad, Oñate [74] formula en su trabajo una segunda propuesta al respecto, planteando que un procedimiento sin malla debiera satisfacer las siguientes condiciones:

  1. La discretización de la función incógnita y sus derivadas deben poder ser definidas solamente en términos de la posición de los puntos localizados dentro del dominio de análisis y de los parámetros especificados en éstos.
  2. La función de ponderación y sus derivadas deben poder ser definidas solamente en términos de la posición de los puntos localizados dentro del dominio de análisis.
  3. No es necesario realizar ninguna integración de volumen o superficie, o
  4. Cualquier integración de volumen o superficie debiera ser independiente del procedimiento de interpolación escogido.

Este último planteamiento es un buen punto de partida para comprender las características y funcionamiento de los distintos métodos numéricos mesh-free existentes, ya que entrega conceptos claves que son comunes a todos ellos: discretización de la función incógnita, función de ponderación, procedimiento de interpolación, e incorpora formalmente la posible utilización de una malla de fondo o celda de integración. Estos aspectos y la similitud que presentan estas técnicas, desde el punto de vista de su formulación matemática, ha sido objeto de análisis y estudio por parte de diversos grupos de trabajo [71] [25] [61] [11].

En el presente capítulo se exponen y examinan los diversos fundamentos matemáticos de las aproximaciones numéricas sin malla, para la resolución de sistemas de ecuaciones diferenciales en derivadas parciales sometidas a condiciones iniciales y de contorno conocidas. El análisis comprende cuatro etapas estrechamente relacionadas con el Método de Puntos Finitos: proceso de interpolación o de aproximación, uso de funciones de ponderación, discretización del sistema de ecuaciones diferenciales e implementación de las condiciones de contorno prescritas.

2.1 Aproximación y funciones de interpolación

Una característica importante de todo método sin malla es, sin lugar a dudas, la forma de obtener la aproximación de la función incógnita o desconocida en el dominio de análisis. Para este propósito, y en lo que resta de este capítulo, consideraremos que la aproximación de la función en el dominio ,, o , se obtiene mediante la siguiente combinación de funciones

(2.1)

El dominio se discretiza por medio de un conjunto de nodos o puntos , con ,,, siendo el número total de puntos y el valor aproximado de la función en el punto , . La función se define localmente para cada subdominio , cumpliéndose:

(2.2)

En los métodos sin malla la función se denomina indistintamente como función de forma o función de interpolación, quedando definida en el subdominio que contiene una cantidad de puntos, tal que .

Para que la función de forma tenga un carácter local, es común en estos métodos el uso de una función de ponderación cuyo valor es distinto de cero, sobre el subdominio relativamente pequeño respecto del resto del dominio . En el argumento de la función de ponderación: indica un punto cualquiera1 del subdominio , referencia el punto donde se está efectuando la aproximación y donde la función de ponderación alcanza un máximo valor. El símbolo se reservará en esta tesis para indicar el subdominio asociado a un punto de coordenadas de la partición, sobre el cual se desea calcular o conocer los parámetros de la aproximación (figura 1). Dicho punto, en el léxico del MPF, se conoce como nodo estrella. Por otro lado, el símbolo o indicará el subdominio asociado a un punto de coordenadas genéricas o respectivamente.

Cada subdominio , dominio de influencia o soporte de la función de ponderación, se define de manera que en 1D es siempre un intervalo, mientras que en 2D y 3D adopta la forma de disco y esfera respectivamente. También es posible emplear rectángulos o paralelepípedos, incluso es posible utilizar las dos geometrías circular y rectangular en un mismo modelo [11]2. Para un caso en 2D, la figura 1 muestra el dominio de análisis , su contorno y algunos subdominios circulares asociados a su respectivo nodo estrella (sean de interior o de contorno). Como para cada nodo de la partición existirá su correspondiente subdominio, la intersección o traslape existente entre ellos será mayor a la mostrada en la figura. Esta característica, que será analizada posteriormente, posibilita el hecho de que un mismo punto usualmente pertenezca a o más subdominios de interpolación

Modelo para un metodo sin malla basado en subdominios circulares
Figura 1: Modelo para un metodo sin malla basado en subdominios circulares

A continuación se describen y examinan los principales tipos de aproximación utilizados en los métodos sin malla, se analizarán sus propiedades y la conexión existente entre las distintas funciones de interpolación utilizadas. El análisis comprende fundamentalmente aquellas técnicas que por su construcción se asemejan al MPF, lo cual facilitará posteriormente su entendimiento desde el punto de vista de su mecánica operativa.

(1) En el caso de las aproximaciones utilizadas en los metodos sin malla, este punto pertenece a la partición

(2) Las funciones de ponderación que permiten generar este tipo de subdominios se analizan en el apartado 2.3

2.1.1 Aproximación por mínimos cuadrados ponderados

Una aproximación por mínimos cuadrados busca ajustar una curva o polinomio a los valores discretos de una función en unos puntos, de forma de minimizar el error cometido en la aproximación. A pesar que este tipo de aproximación difiere de lo que se conoce tradicionalmente como una interpolación, porque el polinomio no se iguala exactamente con los datos en los puntos, en la literatura de los métodos meshless los términos aproximación e interpolación suelen tratarse como sinónimos. Esta técnica ha sido utilizada en la resolución numérica de problemas de mecánica de sólidos y fluidos [68] [47] [6] , para aproximar el campo desconocido o incógnita a través de unos valores nodales. Sin embargo, tal como se demuestra en [71], el éxito de esta aproximación presenta una importante restricción o inconveniente. La aproximación se deteriora rápidamente en la medida que el número de puntos utilizados en la interpolación local, , aumenta o crece demasiado respecto del número de términos en el polinomio base de interpolación, . Para paliar este inconveniente se recurre al uso de una función de ponderación , que permite mejorar la aproximación por ejemplo, en la vecindad del punto donde se requiere calcular la función o su derivada.

Una de las primeras aproximaciones por mínimos cuadrados, en la cual se introdujo el concepto de la función de ponderación, fue utilizada por Lancaster y Salkauskas [42] para representar o generar diversos tipos de superficies, a partir de la interpolación de una función base polinómica (interpolante) sobre un set de puntos aleatoriamente distribuidos en un dominio. Posteriormente este método, denominado por los autores como `interpolating moving least squares method' IMLS, ha sido ampliamente utilizado en el contexto de los métodos sin malla, para obtener una solución numérica aproximada a sistemas de ecuaciones en derivadas parciales [93]. Buen ejemplo de lo anterior, son las distintas familias de métodos sin malla que utilizan la técnica de interpolación tipo `moving least squares' (MLS), como: `diffuse element method' (DEM) [69], `element-free Galerkin methods' (EFGM) [9], `finite point method' (FPM) [71], `meshless local Petrov-Galerkin method' (MLPG) [3], `local boundary integral equation method' (LBIE) [105], y más recientemente `local point interpolation method' (LPIM) [49] y `least-squares collocation meshless method' (LSCMM) [103].

2.1.1.1 Aproximación tipo `moving least squares' MLS

Construcción

En un método MLS, la aproximación local de la función para cada punto , se obtiene mediante una base de funciones linealmente independientes y de un conjunto de coeficientes desconocidos , de la siguiente manera

(2.3)

siendo el operador una aplicación

(2.4)

con los vectores

(2.5)

(2.6)

La base de funciones contiene típicamente monomios que dependen de las coordenadas espaciales (por ejemplo en 2D), mientras los coeficientes dependen de la posición del punto . Como ejemplos, para el caso lineal y cuadrático respectivamente se tiene:

(2.7)

Es importante destacar la posibilidad de incluir otro tipo de funciones en la base, como por ejemplo, funciones que puedan tener singularidades o discontinuidades. Las bases así definidas, se han utilizado para la modelación de problemas en mecánica de fractura [8] [29] [91].

De acuerdo a la expresión 2.3 , para cada punto de la partición con coordenadas , siendo el conjunto de puntos pertenecientes al dominio de influencia del punto donde , existirá una diferencia o error entre el valor de la función y el valor de la aproximación, que puede ser cuantificada mediante1

(2.8)

Se debe notar de 2.8 que el valor de la aproximación, consiste en los términos de la base evaluados en y sus respectivos coeficientes en . Para obtener el vector de parámetros desconocidos , se recurre a la minimización del siguiente funcional discreto

(2.9)

cuya expresión compacta en forma matricial es

(2.10)

donde

(2.11)

(2.12)

(2.13)

La función de ponderación en 2.9, al igual que en todos los métodos `mesh-free', confiere el carácter local a la aproximación. Esto significa que se construirá únicamente con la información que aportan los valores , de los nodos que pertenezcan al subdominio para el cual tenga un valor diferente de cero. Este subdominio o soporte de la función de ponderación, suele tener forma circular (disco o esfera) centrado en y de radio .

Según [63], esta función debe cumplir las siguientes propiedades :

  1. en un subdominio ,
  2. fuera del subdominio ,
  3. Propiedad de normalidad: ,
  4. es una función monótonamente decreciente, donde ,
  5. cuando , donde es la función delta de Dirac.

Desde ahora y en lo que resta de la tesis, cada vez que se haga referencia a una función de ponderación, se entenderá que esta cumple con las cinco propiedades anteriores. En el argumento de la función de ponderación es una medida del tamaño de su soporte2. En el caso de la aproximación MLS, como consecuencia de las propiedades (i) y (ii), únicamente intervendrán en el cálculo de 2.3, aquellos puntos que cumplan . En principio, para una distribución arbitraria de puntos, distintos valores de darán lugar a diferentes funciones de ponderación , lo que se traduce en la dificultad de no contar con una única manera de definirla globalmente.

Los parámetros desconocidos , que minimizan la expresión del error, se obtienen derivando vectorialmente 2.10 respecto del vector e igualando a cero. Con esto se consigue un sistema lineal de ecuaciones

(2.14)

(2.15)

donde las matrices (matriz de momentos), y son

(2.16)

(2.17)

(2.18)

(2.19)

De la ecuación 2.15 se puede obtener el vector como

(2.20)

(2.21)

lo que permite obtener finalmente de 2.3 la expresión de la aproximación local en su forma compacta

(2.22)

y en forma desarrollada

(2.23)

La expresión de las funciones de forma se obtiene agrupando los términos que multiplican a en 2.23

(2.24)

Para extender la aproximación local 2.23 a todo el dominio, se introduce un operador global tal que [63]

(2.25)

donde

(2.26)

lo que permite obtener en el límite , siendo la aproximación global MLS en su forma compacta y desarrollada

(2.27)

(2.28)

y la expresión de las funciones de forma

(2.29)

(1) La desigualdad se justificará cuando se analice la propiedad de existencia en la aproximación MLS

(2) Generalmente el parámetro se omite como argumento de la función de ponderación

2.1.2 Interpolante de Shepard

Considérese el caso particular en que la función base de interpolación contiene como único elemento

(2.30)

en este caso las componentes de la matriz y serán respectivamente (ver 2.17, 2.19)

(2.31)

(2.32)

quedando la aproximación global definida como

(2.33)

y su correspondiente función de forma (recuérdese 2.29)

(2.34)

La aproximación global 2.33 recibe el nombre de interpolante de Shepard [90], y sus funciones de forma o funciones de Shepard cumplen las siguientes propiedades1:

En el lenguaje matemático se dice que la colección de funciones que cumplen la propiedad (ii) representan una partición de la unidad [88].

(1) Las funciones de forma que se construyen en el MEF también cumplen estas propiedades

2.1.3 Propiedades. Existencia de la aproximación.

La existencia de la aproximación MLS está condicionada a la resolución del sistema lineal de ecuaciones 2.20 , lo que se traduce en definitiva en calcular la inversa de la matriz de momentos. Recordando 2.17 , la matriz de un punto para el cual existe un subdominio y un conjunto de puntos donde , puede escribirse como

(2.35)

donde corresponde al producto interno ponderado que depende del punto , definido por [42]

(2.36)

Como muestra 2.35 la matriz es simétrica por construcción, y definida positiva por la propiedad (i) de la función de ponderación estipulada del apartado 2.1.1, si además las funciones son linealmente independientes, la matriz es invertible y la solución del sistema 2.20 existe y es única. No obstante, la independencia lineal de la base de polinomios no es condición suficiente para asegurar la no singularidad de la matriz y por consiguiente la existencia de su inversa. Para que la aproximación sea factible, deberán además respetarse unos requerimientos mínimos en la definición de los subdominios de interpolación, de cara a obtener una distribución admisible de puntos [56]. Una distribución admisible de puntos debe satisfacer los siguientes requerimientos:

  1. El número de partículas contenidas en el subdominio debe ser tal que
    (2.37)
  2. donde es un número que garantiza la existencia de las funciones de forma o regularidad de la matriz 1. Una condición necesaria para que sea una matriz regular es que

    (2.38)

    Recordando 2.13 y 2.16 se aprecia que el rango de la matriz , y en consecuencia de , puede ser como máximo igual a , por consiguiente para que sea definida positiva es necesario que . Lo interesante de la propuesta 2.38 es la definición de un posible valor para . A título de ejemplo, si se elige como base de interpolación , , la figura 2 muestra una serie de subdominios para los que no se satisface la condición .

  3. La distribución de puntos o partículas debe ser no degenerada , en concreto en 2D, significa que , un mínimo de 3 puntos deben pertenecer a y éstos no deben superponerse, es decir sus vectores de posición deben formar un elemento triangular no nulo.
  4. Distribucion de puntos inadmisible. En los subdominios marcados n=card(S(x))<3math
    Figura 2: Distribucion de puntos inadmisible. En los subdominios marcados n=card(S(x))<3math

(1) El valor de sólo tiene sentido a efectos computacionales.

2.1.4 Consistencia

Cuando se estudia la convergencia de un método sin malla, la consistencia de la aproximación utilizada, es un aspecto que debe ser analizado con detenimiento. Entenderemos por orden de consistencia de una aproximación, al grado del polinomio que debe ser representado exactamente, de esta forma la capacidad de reproducir polinomios de grado es equivalente a la consistencia de orden . Los requerimientos de consistencia dependen del orden de las ecuaciones diferenciales parciales que deben ser resueltas y del esquema de discretización empleado.

Para demostrar la consistencia de orden de una aproximación MLS si la base de interpolación está completa en el polinomio de orden , considérese que se quieren aproximar simultáneamente el siguiente conjunto de funciones de la base, agrupadas en el vector , siendo esta vez la matriz con los valores de las funciones en los puntos

(2.39)

Supóngase como antes, que es un punto arbitrario de para el cual existe un subdomino , entonces para cada punto de la partición donde , la aproximación local 2.22 tendrá esta vez la siguiente expresión (recuérdese también exp. final vec. alfa , 2.18)

(2.40)

lo que demuestra el hecho de que cualquier función que aparezca en la base puede ser reproducida exactamente. Si por ejemplo, la base contiene todos los términos constantes y monomios lineales, entonces el orden de consistencia de la aproximación MLS será lineal.

2.1.5 No interpolación

Sin duda el aspecto más destacado de las funciones de forma 2.29 que se acaban de deducir, presente además en todos los métodos sin malla que utilizan una interpolación por mínimos cuadrados, es el no cumplimiento de la propiedad de interpolación, es decir:

(2.41)

siendo el símbolo la delta de Kronecker. Como consecuencia de lo anterior, ver figura 3 , el valor de la función incógnita en el nodo es distinto del valor de la aproximación en ese punto

(2.42)
Por este motivo, los parámetros no deben ser tratados como valores nodales, sino como contribuciones nodales, en el entendido de que cada nodo aporta su contribución en la construcción de la aproximación . El subíndice que acompaña a la función de ponderación en la figura 3 indica el centro del subdominio de influencia y su argumento el punto donde esta se evalúa.
Aproximacion por minimos cuadrados ponderados tipo MLS
Figura 3: Aproximacion por minimos cuadrados ponderados tipo MLS

2.1.6 Continuidad y derivabilidad

La continuidad de la aproximación moving least squares está supeditada a la regularidad de las funciones de la base de interpolación y de la función de ponderación . Se puede comprobar sin mayor dificultad, que para un caso en que y , entonces la aproximación global MLS .

La derivada parcial de la función de forma, en la aproximación MLS, se obtiene como (recuérdese la expresión 2.29):

(2.43)

donde

(2.44)

y el subíndice que sigue a la coma, representa la derivada respecto de la coordenada espacial , es decir . Las derivadas de orden superior se pueden obtener repitiendo el mismo procedimiento, sin embargo, a pesar del grado de sistematización que se puede lograr con este proceso, el coste computacional sigue siendo bastante elevado. Al respecto, vale la pena tener presente el hecho de que para calcular las derivadas parciales de las componentes de las matrices y , es necesario calcular previamente la derivada de la función de ponderación .

Para lograr de forma más sencilla la derivada de la función de forma, Nayroles [69], primer investigador que utilizó la técnica MLS en el contexto de los métodos sin malla, propuso considerar constante el vector de parámetros desconocidos en la aproximación 2.3. Como consecuencia de este planteamiento, la derivada de las función de forma se aproxima como

(2.45)

lo que supone un considerable ahorro de cálculo en comparación con 2.43. Sin embargo, tal como se demuestra en [8], en el contexto del método sin malla (EFG), despreciar el término cuando se utiliza una interpolación por MLS puede provocar errores que afectan la exactitud de la solución. Independiente del método sin malla que se utilice, lo que si parece estar claro es que, cuando se utiliza una interpolación por MLS se debe dedicar un esfuerzo adicional para desarrollar técnicas que permitan evaluar la derivada de la función de forma de una manera razonable, es decir, sin afectar la exactitud ni la versatilidad del método. Al respecto, en trabajos como [12] [41] y [17] se pueden encontrar técnicas que permiten obtener una reducción en el coste computacional de la derivación de la aproximación MLS.

2.1.7 Aproximación con función de ponderación smooth (SPH)

El método sin malla más antiguo que recoge la literatura científica se denomina smooth particle hydrodynamics method, también conocido con el acrónimo de SPH. Desde sus comienzos [59] [34] hasta hoy en día, el método ha sido desarrollado para simular fenómenos astrofísicos como la evolución, rotación y colisiones de estrellas, a partir de un set de puntos o partículas distribuidas de manera irregular. Si bien en sus inicios el método no estaba pensado para ser utilizado en medios continuos, pues la exactitud en sus resultados estaba condicionada a determinados tipos de problemas (número reducido de partículas y ausencia de contornos), paulatinamente se ha comenzado a investigar su utilización en la simulación dinámica de materiales que presentan fragmentación o fractura [87] y en la simulación de procesos de extrusión de metales [15].

En el método SPH la aproximación , se define como [35] [63] [64] [65]:

(2.46)

donde es la aproximación de la función y es la función de ponderación smooth. Es de interés destacar, que la aproximación en el método SPH coincide a efectos prácticos con la definición del producto de convolución de dos funciones

(2.47)

en donde los límites de la integral pueden particularizarse para una zona de interés , si la función cumple los requisitos de la función de ponderación definidos en el apartado 2.1.1. En el método SPH, la función de ponderación utilizada se conoce como spline1 SPH [63].

Para aproximar numéricamente la integral 2.46 el método considera el dominio de solución , dividido en partículas elementales de masas , , ....., , siendo la contribución a la integral de una partícula , cuyo volumen es , masa y su centro de masa ubicado en

(2.48)

donde es la densidad en el nodo . De esta forma se obtiene la aproximación en término de los valores nodales de la siguiente manera

(2.49)

donde representa el conjunto de puntos o nodos de la partición pertenecientes al dominio de influencia de , tal que , y las funciones de forma de la aproximación SPH

(2.50)

Si bien es cierto, la aproximación 2.49 utilizada en el método SPH permite evaluar numéricamente la integral 2.46 sin la necesidad de una malla o una conexión fija entre los nodos, es necesario subdividir el dominio y contar con técnicas robustas que permitan asignar a cada nodo su correspondiente 2. Desarrollar estas técnicas resulta costoso y difícil de abordar, por ejemplo, en dominios 3D de geometrías irregulares.

(1) Las características de esta función de ponderación se analizan en el apartado 2.2

(2) representa una medida del dominio que rodea al nodo

2.1.8 Propiedades de la aproximación SPH

Al igual que en el caso de la aproximación MLS, el método SPH en general no cumple la propiedad de interpolación, es decir , por lo que 2.50 no puede ser entendido como un verdadero interpolante. Sin embargo, es en la propiedad de consistencia de la aproximación donde el método SPH se resiente.

Que la aproximación 2.46 sea consistente, implica por ejemplo, que las funciones constantes y lineales deberán ser representadas de manera exacta. En una dimensión, tomando y , significa que se deben verificar respectivamente las ecuaciones

(2.51)

(2.52)

La ecuación 2.51 o consistencia de orden cero, se verifica automáticamente al coincidir con la propiedad (iii) de normalidad de la función de ponderación. Para aclarar si se verifica la consistencia lineal, es decir 2.52, adviértase que 2.51 implica

(2.53)

restando 2.53 de 2.52, se obtiene

(2.54)

La ecuación 2.54 es el momento de primer orden de la función de ponderación y su cumplimiento obliga a que esta sea simétrica respecto al origen. Como se verá posteriormente, la mayoría de funciones de ponderación satisfacen esta condición y, por lo tanto, la forma continua de la aproximación SPH con función de ponderación `smooth' posee consistencia de primer orden. Sin embargo, lo anterior no garantiza la consistencia lineal, ni siquiera en 1D, de la forma discreta 2.49.

En una dimensión, utilizando la regla del trapecio para la cuadratura numérica de 2.51 y 2.54, la versión discreta de las condiciones de consistencia son

(2.55)

(2.56)

donde para un set de nodos numerados secuencialmente

(2.57)

Se puede comprobar fácilmente, para una distribución no uniforme de nodos, que la ecuación 2.56 no se satisface. Considérese, por ejemplo, la disposición de cuatro nodos que se ilustra en la figura 4(a), donde se prescriben las condiciones de consistencia en . Particularizando 2.56 para el nodo , se tiene

(2.58)

lo que equivale a

(2.59)

Para una función de ponderación como el spline SPH, la expresión 2.59 es igual a , por tanto, la condición de consistencia lineal no se satisface .

La situación se deteriora todavía más en el contorno. Para el caso de la figura 4(b), las condiciones de consistencia lineal en exigen que

(2.60)

lo que supone el cumplimiento de , sin embargo, esto es imposible de lograr con una función de ponderación que se ajuste a las condiciones (i) a (v) del apartado 2.1.1.

Consistencia lineal de la aproximacion SPH para una distribucion no uniforme de nodos
Figura 4: Consistencia lineal de la aproximacion SPH para una distribucion no uniforme de nodos

La derivada espacial de la aproximación SPH, respecto de la coordenada , se obtiene reemplazando por en 2.46 y hallando el estimador de esta nueva relación como

(2.61)

Además se puede demostrar, después de cierta manipulación algebraica [87], que

(2.62)

por lo que la derivada de la aproximación, en su forma continua, se obtiene a partir de los valores de la función y de la derivada de la función de ponderación. Como antes, la forma discreta de 2.62 se obtiene a través de una suma sobre los puntos de interpolación

(2.63)

Finalmente, en el método SPH, las expresiones 2.49 y 2.63 son colocadas en para obtener el sistema discreto de ecuaciones en derivadas parciales, cuya resolución permite obtener los coeficientes buscados .

2.1.9 Aproximación mediante el operador RK

El operador tipo reproducing kernel (RK) o núcleo generador es una clase de operador que permite, integrando sobre una función kernel , reproducir una función como

(2.64)

Un ejemplo clásico de un operador tipo núcleo generador, que sirve para comprender el estudio y desarrollo del método sin malla denominado reproducing kernel particle method (RKPM)[50] [51] [54] [55], es la transformada de Fourier. El método RKPM surge como alternativa a la aproximación SPH, para corregir los problemas de consistencia y de precisión en la solución en los contornos o cuando se utiliza un número pequeño de partículas. De esta forma, la aproximación de la función incógnita sigue un planteamiento similar al presentado en el método SPH, utilizándose esta vez una función de ponderación modificada , o reproducing kernel, que incorpora las correcciones necesarias para mejorar la aproximación.

Para comprender el funcionamiento de la técnica de aproximación utilizada en el método RKPM, a continuación se ejemplifica su uso para reproducir la función , como una suma de funciones linealmente independientes, es decir

(2.65)

donde los vectores y se definen como en 2.5 y 2.6. El vector de parámetros desconocidos , que en el apartado 2.1.1. se ha calculado mediante la técnica MLS, se obtiene multiplicando en ambos lados de 2.65 por y aplicando la integral de la función kernel , es decir1

(2.66)

donde se define la función

(2.67)

que permite obtener

(2.68)

Reemplazando 2.68 en 2.65 se tiene finalmente

(2.69)

donde la window function o función de ponderación modificada viene dada por

(2.70)

Debe destacarse, que si en la expresión 2.69 se escoge la función de corrección , entonces se recupera en su forma original la aproximación SPH. La función de ponderación utilizada en la aproximación RKPM, a diferencia del método SPH, incorpora formalmente el radio de influencia o parámetro de dilatación como argumento, siendo esta vez

(2.71)

A pesar de la diferente notación utilizada, posee las mismas propiedades y características de la función de ponderación definidas en 2.1.1.

La función depende de los distintos momentos de la función de ponderación [50] [52], siendo para el caso unidimensional con una base polinómica lineal2, ,

(2.72)

con

(2.73)

(2.74)

y

(2.75)

(2.76)

(2.77)

Incorporando las anteriores expresiones, en la forma continua de la aproximación RKPM 2.69, se obtiene finalmente para el caso unidimensional con una base de interpolación lineal

(2.78)

expresión que pone de manifiesto, el distinto papel que juegan las funciones y en el interior del dominio y en los contornos. Al respecto vale la pena destacar que, mientras en el interior y (forma original SPH), en el contorno y , lo que indica que la inclusión de la función de corrección en la aproximación RKPM juega un papel importante precisamente en los contornos, donde la consistencia la aproximación SPH se deteriora. Utilizando una cuadratura numérica como la regla del trapecio, la forma discreta de la aproximación RKPM 2.78 viene dada por

(2.79)

siendo la función de forma

(2.80)

de la misma manera, las funciones , y los distintos momentos de la función de ponderación son esta vez

(2.81)

(2.82)

(2.83)

(2.84)

(2.85)

(1) Se introduce como variable de integración

(2) La obtención de la función en un caso 2D o 3D, se consigue bajo el mismo procedimiento [51]

2.1.10 Propiedades de la aproximación RKPM

La continuidad de la aproximación RKPM, por su construcción, está supeditada a la regularidad de las funciones de la base de interpolación y de la window function, además, para el cómputo de las integrales deberá utilizarse una distribución admisible de puntos, tal como se estipuló en la aproximación MLS. La función de forma 2.80 de la aproximación RKPM no cumple con la identidad de la delta Kronecker, es decir , puesto que su valor está ponderado por una función monótona decreciente que se anula únicamente fuera del dominio de influencia del nodo . Sin embargo, la función de corrección permite introducir los cambios necesarios en la aproximación para que esta sea consistente, y así poder reproducir exactamente cualquier función que se incluya en la base de interpolación. Considérese por ejemplo, que se quiere aproximar un conjunto de funciones agrupadas en el vector , aplicando 2.69 se tiene

(2.86)

La demostración anterior supone utilizar la misma cuadratura numérica para calcular y .

La derivada de la función de forma de la aproximación RKPM, para el caso 1D y base de interpolación lineal, se obtiene diferenciando 2.80 respecto de como

(2.87)

lo que supone calcular las derivadas de las funciones y , es decir

(2.88)

y

(2.89)

además de la derivada de la función de ponderación. En las expresiones anteriores se ha utilizado la siguiente notación .

Al igual que en el caso MLS, el coste computacional asociado al cálculo de la derivada de la función de forma, incluso en el caso anterior 1D, no es para nada despreciable (notar además que es necesario calcular la derivada de los distintos momentos de la función de ponderación , y ).

2.1.11 Aproximación mediante diferencias finitas generalizada (DFG)

Una forma alternativa de obtener la aproximación de la función incógnita y sus derivadas, en el contexto de los métodos sin malla, es mediante la expresión general de su desarrollo en serie de Taylor alrededor de un punto del dominio. La técnica así desarrollada se conoce como diferencias finitas generalizadas [45] [47] y ha sido utilizada para resolver diversos tipos de problemas en mecánica aplicada [46] [96] [97]. Con posterioridad esta técnica ha sido utilizada por Liszka y Orkisz, en conjunto con el método de mínimos cuadrados ponderados, para la resolución de problemas de contorno mediante el método sin malla denominado hp-Meshless cloud method [48] [83]. Utilizando un esquema de diferencias finitas generalizadas, y con la notación empleada en esta tesis, la aproximación de la función en un entorno del punto se construye como1

(2.90)

siendo los vectores

(2.91)

(2.92)

La expresión 2.90 puede ser interpretada como la expansión en serie de Taylor, de la aproximación , alrededor del punto . El vector , que contiene los valores de la función desconocida y sus derivadas en , se calcula empleando la técnica de mínimos cuadrados ponderados, de forma semejante a lo expuesto en el apartado 2.1.1. Esta vez, el funcional discreto que minimiza el error cuadrático ponderado es (recuérdese 2.9)

(2.93)

Derivando 2.93 respecto del vector e igualando a cero, se consigue un sistema lineal de ecuaciones que permite finalmente obtener los parámetros desconocidos . Como el proceso para obtener estos parámetros es semejante al utilizado en la aproximación MLS, no se ha estimado oportuno repetirlo nuevamente para este caso (DFG), además, si en el cálculo de 2.93 intervienen todos los nodos pertenecientes al subdominio de interpolación , las aproximaciones 2.3 y 2.90 son equivalentes. Si por el contrario, no se utilizan todos los nodos que pertenecen al subdominio, las técnicas MLS y DFG entregarán distintas funciones de aproximación. Al respecto debe notarse, que en esta técnica la evaluación de los coeficientes desconocidos requieren de la conectividad de al menos 6 nodos.

(1) Se ejemplificará, sin pérdida de generalidad, la utilización de DFG para un caso 2D

2.1.12 Aproximación tipo partición de la unidad (PU)

Utilizando el concepto de partición de la unidad, también es posible construir una aproximación de la función desconocida, para que pueda ser utilizada en un método sin malla. Este planteamiento, propuesto inicialmente por [24] y [5], permite además comprender bajo un aspecto más general el funcionamiento de los distintos métodos sin malla. En una partición de la unidad (PU), el dominio es cubierto por un número finito de subdominios que se superponen, asociándose a cada uno de ellos una función diferente de cero sólo sobre (notar la semejanza entre la proposición anterior y la definición de subdominios circulares utilizadas en los métodos sin malla, figura 1). Además, si el dominio se discretiza mediante un conjunto de puntos , , siendo el número total de puntos, se dice que la familia de funciones representan una partición de la unidad respecto del conjunto de subdominios si

(2.94)

Debe destacarse que la propiedad anterior de las funciones , es idéntica a la condición de consistencia de orden cero (2.55) que deben cumplir las funciones de forma en una aproximación SPH. También las funciones de forma en una aproximación MLS son una partición de la unidad, puesto que por condición de consistencia cumplen

(2.95)

siendo el grado del monomio de la base de interpolación, en particular, si (es decir función constante) se tiene

(2.96)

De esta forma, se pueden construir particiones de la unidad a partir de las aproximaciones tipo núcleo generador y MLS analizadas anteriormente. En particular, en el método denominado hp clouds [26], se utilizan las funciones de forma MLS (2.29) para construir la partición unitaria, formando subdominios o nubes asociadas a parámetros y que permiten plantear, al igual que en el MEF, procedimientos adaptativos de la solución aumentando el número de puntos utilizados en la discretización o por el contrario aumentando el grado del polinomio de interpolación. En el método hp, la aproximación de la función viene dada por

(2.97)

donde el superíndice de las funciones de forma indica el orden del polinomio base de interpolación1, agregándose a la formulación clásica MLS, el conjunto de funciones que constituyen la base extrínseca y que contienen polinomios2 de orden superior al de las funciones (base intrínseca), o cualquier otro tipo de función que se considere adecuada para la aproximación (enhancement functions). Debe destacarse que en el caso de la aproximación utilizada en el método , además de calcular los parámetros para la base intrínseca (recuérdese el proceso de inversión de la matriz en la técnica MLS), es necesario resolver un sistema lineal de ecuaciones para obtener el conjunto de parámetros de la base extrínseca, lo que se traduce en definitiva en un aumento del coste computacional.

La idea principal de la formulación 2.97, es poder añadir elementos de forma jerárquica a la familia de funciones que representan la partición de la unidad , de manera que el nuevo set de funciones pueda reproducir polinomios de grado . Además, la base extrínseca puede ajustarse para cada nodo añadiendo términos extra, sin que las condiciones de continuidad y derivabilidad se vean afectadas [11], con ello se consigue la implementación de esquemas de refinamiento tipo . Para implementar un refinamiento tipo , los autores proponen un estimador del error en la solución a posteriori, que permite introducir nuevos nodos en aquellas zonas de interés del dominio [26].

Respecto de las propiedades de la aproximación tipo partición de la unidad, es fácil verificar que se deben respetar las mismas condiciones que garanticen la existencia de las funciones de la base intrínseca. Para el caso en que estas correspondan a las funciones de forma MLS, serán las indicadas en el apartado 2.1.1.

(1) El caso particular , corresponde a las funciones de forma del interpolanate de Shepard

(2) En [24] se utilizan polinomios de Legendre

2.2 Funciones de ponderación

Cuando se analizan las distintas características de las aproximaciones utilizadas en los métodos sin malla, es fácil identificar un aspecto común en todas ellas, esto es, la utilización de una función de ponderación . Esta función, además de conferir el carácter local a la aproximación, permite distribuir o ponderar el error cometido en la interpolación, controlando el tamaño del subdominio o nube según su radio de acción. En principio su elección es arbitraria, siempre que se respeten las condiciones (i) a (v) del apartado 2.1.1 y que la función y sus derivadas sean continuas en el grado deseado. El grado de continuidad requerido dependerá, en general, del tipo de aproximación meshless que se utilice, del orden de la ecuación diferencial del problema y del esquema de discretización empleado.

La función de ponderación, por construcción, posee una forma bastante característica, es decir presenta un valor máximo en el entorno del punto donde se requiere evaluar la aproximación, con una tendencia decreciente a medida que aumenta la distancia al punto en cuestión. Recordando que las aproximaciones analizadas se encuentran asociadas a la minimización del error a través de un funcional, la interpretación práctica de esta tendencia o forma de la función de ponderación no es otra que buscar penalizar el error cometido en la interpolación, proporcionalmente a la distancia del punto donde se está evaluando. La forma y tamaño de la función de ponderación se pueden regular a través de ciertos parámetros, propios de cada función, que entre otras cosas deben además garantizar la no singularidad del problema. Uno de estos parámetros, que además proporciona una idea del tamaño de los subdominios, es el que se conoce como radio de influencia o parámetro de dilatación . El radio de influencia juega un papel importante en la selección de los subdominios controlando el número de puntos que lo integran, de modo que , siendo la dimensión de la base de interpolación, además su valor tiene una cota máxima para conseguir una mayor eficiencia computacional.

Para conseguir una completa y adecuada definición de todos los parámetros que intervienen en una función de ponderación, existen dos alternativas:

  1. Definir a priori el número de puntos que componen el subdominio y a partir de aquí sacar los parámetros de la función (como por ejemplo el radio de influencia asociado a cada subdominio).
  2. Determinar o fijar los parámetros y a partir de ellos obtener los puntos que componen el subdominio.

Cuando se utiliza una aproximación basada en MLS, la primera alternativa resulta más segura de cara a cerciorarse que se tendrá en todos los casos, el número de puntos suficientes en el subdominio garantizando de esta forma la no singularidad del problema, acorde con la dimensión de la base de interpolación. Esta alternativa, será la utilizada en el desarrollo de los ejemplos presentados en esta tesis.

2.2.1 Propiedades y construcción

Debido al carácter local de la aproximación, las funciones de ponderación deberán ser no nulas únicamente en un subdominio cumpliéndose

(2.98)

donde es la distancia entre el punto y un punto de la partición, pertenecientes al subdominio . En lo sucesivo denotará la función de ponderación asociada a un nodo , es decir cuyo valor máximo se encuentra en . Consultando los diversos desarrollo sobre métodos sin malla, las funciones de ponderación mayoritariamente utilizadas en la práctica resultan ser del tipo [8],[93]

(2.99)

con continua, así como sus primeras derivadas. A continuación se estudian las condiciones que debe satisfacer (entero positivo) para que las primeras derivadas de la función de ponderación , con respecto a cada componente del vector , sean continuas en dicho punto.

Considérese, en primer lugar, la derivada primera de la función respecto de una componente cualquiera del vector de posición

(2.100)

Adviértase que el límite de , a medida que tiende a , no existe. Por consiguiente, la derivada primera sólo existirá para .

Para la derivada segunda, de la función de ponderación respecto de , se tiene

(2.101)

La derivada anterior existe siempre que (para valores inferiores el segundo sumando del miembro derecho de 2.101 da lugar a problemas de continuidad). Por inducción, se puede demostrar que la derivada ésima de la función de ponderación respecto de cada componente del vector existe si [9].

2.2.2 Tipos de funciones de ponderación

Habiendo revisado las principales características teóricas que han de verificar estas funciones, a continuación, se presentan las funciones de ponderación comunmente más utilizadas en las distintas aproximaciones sin malla. Considerando en primer término las funciones de ponderación cuyo soporte es circular o esférico (ver figuras 1,2), existe una función de interpolación que se aplica específicamente en la aproximación SPH, conocida como spline SPH [63], cuya expresión es

  • Spline SPH

(2.102)

donde , siendo una medida del tamaño de su soporte.

Utilizando como argumento de la función de ponderación ( se escoge igual al radio de influencia ), además del spline SPH, se destacan por su mayor utilización las siguientes funciones de ponderación

  • Función triangular

(2.103)

la que ha pesar de su sencillez tiene la desventaja que su derivada no es continua en .

  • Función cónica

(2.104)

la que constituye una generalización de la función triangular, su continuidad aumenta paralelamente con el exponente , aunque en la práctica . Según Duarte [25], esta función da lugar a funciones de forma que pueden ser integradas con mayor precisión en el caso de utilizar como esquema de discretización el método de Galerkin1.

  • Función sinusoidal

(2.105)
  • Spline de tercer orden

(2.106)
  • Spline de cuarto orden

(2.107)

tanto el spline de tercer como de cuarto orden se utilizan para problemas en los que se requiere .

  • Función exponencial o de Gauss

(2.108)

el parámetro (que se denomina factor de apuntamiento) y el exponente , determinan la forma de la función de ponderación, en definitiva los pesos relativos. En el caso del factor de apuntamiento, su valor determina por ejemplo, el que los pesos de la función de ponderación sean mayores, tanto cerca como lejos de , en la medida en que aumenta. En general, no existe un criterio o método para fijar , siendo su valor más bien arbitrario. En la literatura pueden encontrarse diversas proposiciones o recomendaciones para fijar su valor, así por ejemplo, Belytschko [8] sugiere

(2.109)

con

(2.110)

siendo el menor conjunto de puntos necesarios para determinar un polígono alrededor de . Oñate [74] propone determinar en función del radio de influencia como

(2.111)

mientras que Hegen [36] utiliza

(2.112)

Con posterioridad Atluri [4] propone la utilización de la siguiente igualdad para la selección de los parámetros y de la función de Gauss

(2.113)

siendo , , una longitud característica y el area del dominio de análisis, el número total de nodos utilizados y como antes el número de términos utilizados en la función base de interpolación, para el caso de una aproximación MLS. Reordenando la expresión 2.113, el valor para el factor resulta ser esta vez

(2.114)

La función de Gauss constituye en la práctica, junto con el spline de cuarto orden, una de las funciones de ponderación de mayor utilización en los diversos trabajos de investigación realizados sobre el desarrollo de métodos sin malla.

También es posible construir funciones de ponderación, a partir del producto tensorial de funciones como

(2.115)

donde esta vez el soporte de la función será de forma rectangular, con dimensiones , y en la dirección de los ejes coordenados , y respectivamente. La función así definida ha sido utilizada principalmente, en el caso de dominios que por su geometría favorecen una discretización mediante una distribución regular de nodos. Finalmente, a modo de comparación, en la figura 5 se muestra una representacion de las diferentes funciones de ponderación para un caso 1D. Es de interés destacar, tal como lo indica la figura, la capacidad de la función de Gauss de poder representar, mediante una adecuada selección de los parámetros y , algunas de las funciones de ponderación estudiadas anteriormente.

Funciones de ponderacion para un caso 1D
Figura 5: Funciones de ponderacion para un caso 1D

(1) Véase apartado 2.3.1

2.3 Implementación numérica

Analizados los procedimientos para obtener la función aproximada en un método sin malla, a continuación se revisan las técnicas para obtener el sistema discreto de ecuaciones diferenciales, que permitirá la resolución numérica del problema de contorno asociado. Para el planteamiento del sistema discreto de ecuaciones, se han utilizado hasta la fecha, dos tipos de formulaciones:

  • método de Galerkin
  • método de colocación puntual

El optar por una u otra estrategia, como se verá, presenta sus ventajas e inconvenietes y condiciona el desarrollo de un método sin malla. Considérese la forma general de un problema vectorial gobernado por las siguientes ecuaciones diferenciales

(2.116)

con su condición de contorno de Neumann (natural)

(2.117)

y condición de Dirichlet (esencial)

(2.118)

que debe satisfacerse en un dominio con contorno . En las expresiones anteriores, y son operadores diferenciales apropiados, será el vector de incógnitas (o campo de desplazamientos en el ámbito de la mecánica estructural) y el valor prescrito de a lo largo del contorno . Además y , representan flujos o fuerzas externas actuando sobre el dominio y a lo largo del contorno , respectivamente.

Un procedimiento general para resolver numéricamente el problema de contorno anterior es el método de los residuos ponderados [106], cuya técnica permite obtener una solución aproximada, a partir de una ecuación integral equivalente al sistema de ecuaciones diferenciales del problema. Si la solución exacta se aproxima por , por ejemplo utilizando cualquiera de las aproximaciones meshless estudiadas, es decir

(2.119)

se tiene que, en general, las ecuaciones 2.116, 2.117 y 2.118 no serán satisfechas, obteniéndose unos residuos o errores tanto en el dominio como en el contorno . La solución aproximada al problema de contorno original se consigue ponderando el error cometido en la aproximación mediante funciones de prueba o de test como sigue

(2.120)

donde , , son las denominadas funciones de test, es un subespacio finito de Sobolev y el número de incógnitas del problemas. Tomando en consideración que los integrandos , , representan el error cometido al sustituir la solución aproximada en la ecuación diferencial o en las condiciones de contorno, la expresión 2.120 puede entenderse como la integral ponderada de tales residuos.

2.3.1 Discretización mediante formulación de Galerkin

Esta formulación se basa en escoger como funciones de test, las funciones de forma utilizadas en la aproximación. Su formulación, en principio, no difiere sustancialmente de la utilizada en el método de elementos finitos, sin embargo, para poder imponer las condiciones de contorno escenciales o de Dirichlet es necesario implementar un procedimiento adicional. Para comprender las particularidades que presenta la formulación de Galerkin en una aproximación meshless, a continuación se desarrolla su implementación en un problema de contorno tipo, regido por la ecuación de Laplace en su versión escalar, es decir

(2.121)

y condiciones de contorno Neumann1

(2.122)

y condiciones de contorno Dirichlet

(2.123)

Aplicando la ecuación de residuos ponderados 2.120, en este caso particular se tiene

(2.124)

integrando por partes y aplicando el teorema de la divergencia se obtiene la forma débil de la ecuación 2.124 como [106]

(2.125)

donde . En principio, cualquiera de las aproximaciones meshless estudiadas en 2.1 puede ser utilizada en la expresión 2.125, sin embargo, el no cumplimiento de la condición de interpolación que caracteriza a estas aproximaciones, es decir

(2.126)

impide imponer las condiciones de contorno escenciales cuando y . Para hacer frente a este problema los investigadores han debido implementar diferentes soluciones, como por ejemplo, multiplicadores de Lagrange [9] [58] [66], acoplamiento con elementos finitos [10] [53] [36] [38] y métodos de penalización [104], [31]. Para tomar conciencia, del coste computacional adicional que significa tener que implementar un proceso para prescribir la condición de contorno de Dirichlet, a continuación se ejemplifica para el problema tipo, el uso de multiplicadores de Lagrange. Para ello se introduce la expresión de los multiplicadores de Lagrange, como función de los desplazamientos, es decir

(2.127)

modificándose la expresión 2.125 como sigue

(2.128)

si se adoptan funciones de test de modo que , y definida sobre se tiene

(2.129)

donde los multiplicadores de Lagrange en conjunto con la función , pasan a ser las incógnitas del problema. Utilizando una aproximación meshless para la discretizar , una aproximación para basada en funciones de forma Lagrangianas, es decir

(2.130)

donde es una coordenada que mide la longitud de arco sobre el contorno y el conjunto y seleccionando las funciones test según el método de Galerkin

(2.131)

se obtiene finalmente el siguiente sistema de ecuaciones en forma matricial

(2.132)

donde

(2.133)

(2.134)

(2.135)

(2.136)

Como primera consecuencia, el uso de multiplicadores de Lagrange para imponer las condición de contorno escencial, supone un aumento del ancho de banda de la matriz respecto de la formulación de Galerkin utilizada en el método de elementos finitos. Esto, sumado al inconveniente de que la matriz final no es definida positiva, a pesar de su simetría [8], conlleva un encarecimiento del coste computacional asociado a la resolución del sistema.

(1) representa la derivada respecto de la normal

2.3.2 Evaluación numérica de las integrales

Una vez implementada la solución para imponer la condición de contorno, surge la necesidad de evaluar numéricamente las integrales 2.133, 2.134, 2.135 y 2.136, respetando en la medida de lo posible la filosofía de los métodos sin malla. Las estrategias que se han desarrollado para solucionar este dilema, pueden clasificarse en tres tipos:

  1. Integración nodal, Integración por subdominios
  2. Celda de integración
  3. Malla auxiliar de elementos

La integración nodal es la más fácil y rápida de implementar. Al igual de lo que ocurre en las aproximaciones SPH, se aprovecha la misma partición de puntos del dominio para realizar la evaluación de las integrales, utilizándose expresiones del tipo

(2.137)

(2.138)

donde y son una medida del volúmen y superficie respectivamente, del entorno que rodea al punto . Poder seleccionar y asignar adecuadamente estos parámetros no es para nada una tarea trivial, sobre todo en geometrías irregulares y en 3D, pero quizás el mayor inconveniente de la integración nodal sea la aparición de fenómenos de inestabilidad. Este tipo de inestabilidad, conocida como tensile instabilities, fue inicialmente identificada en las aproximaciones SPH por [92]. Con posterioridad este fenómeno ha sido objeto de estudio por parte de diferentes grupos de investigación, al respecto, referencias como [DYK 95] [28] [7] [87] [23] [15] [13], reflejan el esfuerzo dedicado por los investigadores en identificar el origen de estas inestabilidades y plantear posibles soluciones.

En el procedimiento de integración por subdominios, el cómputo de las integrales se realiza a través de una cuadratura numérica sobre un subdominio definido de forma local, el cual tiene generalmente una forma geométrica sencilla como una esfera, cubo o elipsoide para facilitar la integración [3] [104] [21] [49]. Si bien el método resultante se considera libre de malla, la definición de estos subdominios de integración, especialmente en el contorno, no es una tarea sencilla y requiere de técnicas especiales.

La segunda y tercera alternativa tienen la desventaja de que el método resultante no es verdaderamente un método sin malla, puesto que en ambos se recurre a una cuadratura numérica, como las conocidas cuadraturas de Gauss o regla del trapecio, para evaluar las integrales. Previamente, se debe definir una malla auxiliar formada por una red de celdas o elementos finitos, como se muestra en la figura 6, en donde se definen los puntos de integración . Se puede también observar de la figura, que en el caso de las celdas de integración estas se definen, independiente de la posición de los puntos, como un arreglo regular, a diferencia de la malla auxiliar de elementos finitos cuyos nodos coinciden con los vértices del elemento. A pesar que la evaluación de las integrales es independiente del proceso de interpolación elegido, el hecho de particionar el dominio ya sea con un CAD o con un mallador, ensombrece notablemente las posibles ventajas de los metodos sin malla. Además, en el caso de utilizar una cuadratura de Gauss, no existe un criterio respecto del número de celdas y números de puntos