Detalles técnicos

A la hora de aplicar estos métodos a una simulación de partículas hay que afrontar una serie de problemas técnicos.

El principal es que, a diferencia de las distribuciones de nodos consideradas en las Refs. [31,33], aplicadas a problemas de elasticidad, las configuraciones de las partículas en una simulación pueden ser bastante desordenadas. Por ello es preciso determinar adecuadamente el valor de $g$ . Como se ve en la figura 7.1, en un sistema con dos zonas de distinto valor de $h$ un cálculo en el que no se ha tenido en cuenta este hecho produce unos resultados poco aceptables: la partícula señalada en rojo a la izquierda tiene muchos vecinos, pero la de la derecha tiene pocos.

Para ser más precisos, es $h$ lo relevante; en una implementación numérica los exponentes $f_a$ de la Ec. 6.10 se evaluarían así:

$$f_a= (\lambda h_a) \frac{x-x_a}{h_a} - (\beta h_a)^2 \left[ \frac{(x-x_a)^2}{h_a^2}- \frac{g}{h_a^2}\right] ,$$ (7.5)

de modo que todas las distancia se dividen por el valor de $h_a$ local, y el valor de $g^*=g/h_a^2$ , que es adimensional, está fijo para todo el sistema. Es decir, el valor de $g$ en cada punto es $g_a=g^* h_a^2$ : más pequeño si las partículas están juntas y más grande si están separadas. Esta elección solucionaría la situación de la figura 7.1.

En una dimensión $g^*$ puede tomarse alrededor de $0.25$ , como se sugiere en [33]; en 2D, puede ser sencillamente $0.25 I$ , proporcional al tensor identidad. Una posibilidad, mencionada en [33] y examinada por nosotros (con poco éxito, de momento), es considerar tensores más complejos, que se adapten a la anisotropía local.

Figura 7.1: Funciones base SME para un sistema con dos valores de $h$ muy distintos. De [33].

Así pues, una etapa previa en cada paso temporal de la simulación es la determinación de un valor de $h$ local. Una opción ingenua sería tomar el valor del vecino más próximo. Esta idea hay que descartarla, ya que en general los otros vecinos pueden ser mucho más lejanos, y $h$ debe ser lo bastante grande para englobar un número ``suficiente'' de vecinos.

De hecho, en una dimensión un criterio que da mejores resultados es tomar la distancia del primer vecino (a la derecha o izquierda) que esté más alejado. Esto se generalizaría en 2D de manera elegante: tomar $h$ como la distancia al vecino de Delaunay más lejano; de esta manera, los más cercanos también quedan englobados. El parentesco entre las funciones base y los elementos finitos sobre la triangulación de Delaunay (aunque haya quedado más difuminado que para el método LME) hace que esta elección sea natural. Hemos comprobado que este criterio da buenos resultados; por otro lado, esto exige el cálculo de la triangulación de Delaunay, que justamente se pretendía evitar.

Así pues, en 2D hemos utilizado un método más sencillo, que produce resultados análogos: tomar $h$ como la distancia de la sexta partícula más cercana. Que sea la sexta justamente está motivado porque, en promedio, una partícula en una configuración aleatoria tiene seis vecinos; en todo caso, con 5 o 7 vecinos no cambian los resultados apenas. Esta distancia se determina, en la implementación actual, de una manera poco refinada: se evalúan, para una partícula dada, las distancias a todas las demás, se ordenan y se toma el sexto valor. Al menos, la lista de distancias puede aprovecharse para almacenar una lista de vecinos, descartando los que estén más allá de, aproximadamente $3h$ .

En cuanto al