Métodos sin mallas

Vista previa del material en texto

UNVERSIDAD DE CARABOBO
FACULTAD EXPERIMENTAL DE CIENCIAS Y TECNOLOGÍA
DEPARTAMENTO DE MATEMÁTICAS
MÉTODOS SIN MALLAS.
Trabajo de Ascenso presentado ante la Universidad de Carabobo por
Máximo Gonzalo Mero Barcia
Como requisito para optar a la Categoŕıa de Profesor Asociado
2014
Índice general
Portada I
Resumen V
1. Introducción y Planteamiento del Problema 1
2. Métodos sin Mallas 4
2.1. Smoothed Particle Hydrodynamics - SPH . . . . . . . . . . . . . . . . . . . . 5
2.2. Mı́nimos Cuadrados Móviles - MLS . . . . . . . . . . . . . . . . . . . . . . . 8
3. Aproximación de nubes de puntos usando métodos sin mallas 12
3.1. Mı́nimos Cuadrados (LS) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3.2. Mı́nimos Cuadrados Ponderados (WLS) . . . . . . . . . . . . . . . . . . . . . 14
3.3. Mı́nimos Cuadrados móviles (MLS) . . . . . . . . . . . . . . . . . . . . . . . 17
3.3.1. Aplicaciones de Mı́nimos Cuadrados Móviles . . . . . . . . . . . . . . 19
3.4. Implementación y resultados . . . . . . . . . . . . . . . . . . . . . . . . . . 20
4. Modelos deformables usando MLS 24
i
4.1. Teoŕıa de Elasticidad Lineal . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
4.2. Métodos de Mallas Libres MLS . . . . . . . . . . . . . . . . . . . . . . . . . 25
4.2.1. Aproximación de ∇~U empleando el método de MLS . . . . . . . . . . 26
4.2.2. Funciones de Peso . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
4.2.3. Cálculo de las Fuerzas . . . . . . . . . . . . . . . . . . . . . . . . . . 28
4.3. Dinámica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
4.3.1. Dinámica de la Superficie del objeto 3D . . . . . . . . . . . . . . . . 32
4.4. Resultados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
5. Simulación de Fluidos Basados en Part́ıculas usando SPH 39
5.1. Las Ecuaciones de Navier-Stokes . . . . . . . . . . . . . . . . . . . . . . . . . 39
5.1.1. Enfoque Euleriano . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
5.2. Smoothed Particle Hydrodynamics (SPH) . . . . . . . . . . . . . . . . . . . 41
5.3. Enfoque Lagrangiano . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
5.3.1. Masa-Densidad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
5.3.2. Fuerzas Internas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
5.3.3. Viscosidad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
5.3.4. Fuerzas Externas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
5.4. Integración Numérica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
5.5. Resultados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
6. Conclusiones 64
ii
Índice de figuras
3.1. Aproximación por Mı́nimos Cuadrados . . . . . . . . . . . . . . . . . . . . . 21
3.2. Aproximación por Mı́nimos Cuadrados ponderados . . . . . . . . . . . . . . . 21
3.3. Aproximación por Mı́nimos Cuadrados Móviles . . . . . . . . . . . . . . . . 22
3.4. Aproximación por SPH . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.5. MLS con 256 datos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
4.1. Objeto en estado de reposo y en estado de deformación. . . . . . . . . . . . . 34
4.2. Función de Lennard-Jones L(α, dij). Función W1(r, h). . . . . . . . . . . . . 36
4.3. Funciones de peso W2(r, h) y W3(r, h). . . . . . . . . . . . . . . . . . . . . . 36
4.4. Funciones de peso W4(r, h) y W5(r, h). . . . . . . . . . . . . . . . . . . . . . 36
4.5. Diferencia entre hi y α. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
4.6. Animación de un octaedro con 299 part́ıculas internas. . . . . . . . . . . . . 37
4.7. Animación de un elipsoide con 568 part́ıculas internas. . . . . . . . . . . . . 37
4.8. Animación de un elipsoide formado por 715 part́ıculas internas y 990 de la
superficie y 1976 triángulos. . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
5.1. Configuración Básica de SPH para las ecuaciones de Navier-Stokes . . . . . 46
iii
5.2. Núcleo por Defecto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
5.3. Gráfica del núcleo de Presión . . . . . . . . . . . . . . . . . . . . . . . . . . 60
5.4. Gráfica del núcleo de Viscosidad . . . . . . . . . . . . . . . . . . . . . . . . . 60
5.5. Función de Peso . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
5.6. Simulación en Progreso. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
5.7. Manejo de Colisión. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
5.8. Dinámica con 65536 Part́ıculas . . . . . . . . . . . . . . . . . . . . . . . . . 63
iv
Índice de cuadros
4.1. Equipos utilizados. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
4.2. Tiempos de las animaciones basadas en el modelo de part́ıculas internas
usando la estrategia de búsqueda de vecinos de todos contra todos en los
distintos equipos. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
4.3. Parámetros usados en las animaciones basadas en el modelo de los puntos
internos y de superficie. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
5.1. Resultados para HW1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
5.2. Resultados para HW2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
5.3. Resultados para HW3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
v
RESUMEN
Los métodos sin mallas (Mesh Free -MFree) son métodos de aproximación, es decir, métodos
para reconstruir funciones y sus derivadas a partir de los valores puntuales en una serie
de part́ıculas irregularmente distribuidas, y de cierta información geométrica de bajo nivel
acerca de la distribución de part́ıculas. La reconstrucción requiere la existencia de una
interconexión de más alto nivel entre part́ıculas, esto es, la reconstrucción es independiente
de que las part́ıculas sean o no los nodos de una red.
Por extensión, se denomina habitualmente métodos sin mallas a una serie de métodos
numéricos para resolver ecuaciones diferenciales, que trabajan directamente en el espacio
f́ısico y que obtienen su marco de representación espacial a través de un método
de aproximación sin malla. En la simulación dinámica se requiere resolver ecuaciones
diferenciales complejas que gobiernan estos fenómenos. Tradicionalmente, tales ecuaciones
diferenciales en derivadas parciales se resuelven usando métodos numéricos tales como el
Método de los Elementos Finitos (FEM) y el Método de Diferencias Finitas (FDM). En estos
métodos, el dominio espacial donde las ecuaciones diferenciales parciales están definidas,
usualmente se discretiza en mallas. En FDM, la malla es usualmente llamada cuadŕıcula
(grids). En el Método de Volumen Finito (FVM), estas mallas son llamadas volúmenes o
celdas, y en FEM las mallas son llamadas elementos. Sin embargo, todas estas cuadŕıculas,
volúmenes, celdas y elementos pueden ser calificadas como mallas. La clave es que una malla
debe estar predefinida para proveer cierta relación entre los nodos, la cual es la base de la
formulación de estos métodos numéricos convencionales.
Existen muchos métodos Mfree, tales como el método Element Free Galerkin (EFG) [5], el
método sin malla local de Petrov-Galerkin (MLPG) [3], el método de interpolación puntual
(PIM) [8], el método Smoothed Particle Hydrodynamics [16].
Esta investigación se centra en el estudio de los métodos Smoothed Particles Hydrodynamics
(SPH) [15], Moving Least-Squares (MLS),( [11], [14], [8]), Least-Squares (LS) y Weighted
Least Squares (WLS) [22], su aplicación en los modelos deformables basados en la teoŕıa de
elasticidad lineal y la simulación dinámica de los fluidos.
Caṕıtulo 1
Introducción y Planteamiento del
Problema
Si trabajamos directamenteen el espacio f́ısico un método numérico, es necesario disponer
de un marco de aproximación espacial, este ha de tener asociada una cierta representación
funcional. De esta forma, mientras los métodos de diferencias finitas trabajan en un
espacio computacional transformado, de forma que las ecuaciones se resuelvan en una malla
rectangular, otros métodos numéricos, como el método de los volúmenes finitos o el método
de los elementos finitos, toman la solución numérica de ciertos subespacios de funciones
capturadas por la malla. Por ejemplo, la representación básica en un método de volúmenes
finitos es la de una función constante a trozos, es decir, es constante en cada volumen de
control.
El principal atractivo de los métodos que trabajan en el espacio f́ısico es la posibilidad de
emplear mallas no estructuradas. Se dice que una malla es estructurada cuando se dispone
de una transformación desde el espacio f́ısico, representado por la red, hacia un espacio
computacional como el empleado con el método de diferencias finitas. Si la geometŕıa es
sencilla, se podrá generar una malla estructurada con relativa facilidad. Sin embargo, la
generación de mallas estructuradas, o estructuradas por bloques, para geometŕıas muy
complejas resulta, con frecuencia, en un problema inabordable en términos de ingenieŕıa
práctica. En estos casos, resulta mucho más eficiente el empleo de mallas no estructuradas,
a las que solamente se requiere cubrir la totalidad del dominio con una serie de celdas que
no se solapen.
Estas mallas, entendidas como conjuntos de celdas conectadas en el sentido de una red,
pueden ser empleadas por los métodos que trabajan en el espacio f́ısico para desarrollar
distintos aspectos de la formulación. En los métodos de elementos finitos o elementos de
1
contornos, la conectividad de la malla sirve para construir la aproximación espacial (definir
funciones de forma), además de permitir la discretización de la forma débil del problema
en términos de contribuciones elementales (elemento a elemento). El método de volúmenes
finitos no utiliza, en principio, unas funciones de forma (salvo, quizás, si suponemos que
utiliza funciones constantes en cada celda). Sin embargo, el papel de la discretización de
celdas (volúmenes de control) está intŕınsecamente relacionado con el propio planteamiento
del esquema, a través de conceptos como forma integral del problema en cada volumen de
control y flujo a través del contorno de cada volumen de control.
Los métodos sin mallas son métodos de aproximación, es decir, métodos para reconstruir
funciones y sus derivadas a partir de los valores puntuales en una serie de part́ıculas
irregularmente distribuidas, y de cierta información geométrica de bajo nivel acerca de la
distribución de part́ıculas (por ejemplo, la distancia entre ellas). La reconstrucción requiere la
existencia de una interconexión de más alto nivel entre part́ıculas, esto es, la reconstrucción
es independiente de que las part́ıculas sean o no los nodos de una red.
Por extensión, se denomina habitualmente métodos sin mallas a una serie de métodos
numéricos para resolver ecuaciones diferenciales, que trabajan directamente en el espacio
f́ısico (capaces de utilizar información no estructurada), y que obtienen su marco
de representación espacial a través de un método de aproximación sin malla. En la
simulación dinámica se requiere resolver ecuaciones diferenciales complejas que gobiernan
estos fenómenos. Tradicionalmente, tales ecuaciones diferenciales en derivadas parciales se
resuelven usando métodos numéricos tales como el Método de los Elementos Finitos (FEM)
y el Método de Diferencias Finitas (FDM). En estos métodos, el dominio espacial donde
las ecuaciones diferenciales parciales están definidas, usualmente se discretiza en mallas. En
FDM, la malla es usualmente llamada cuadŕıcula (grids). En el Método de Volumen Finito
(FVM), estas mallas son llamadas volúmenes o celdas, y en FEM las mallas son llamadas
elementos. Sin embargo, todas estas cuadŕıculas, volúmenes, celdas y elementos pueden ser
calificadas como mallas. La clave es que una malla debe estar predefinida para proveer cierta
relación entre los nodos, la cual es la base de la formulación de estos métodos numéricos
convencionales.
Al usar una malla predefinida apropiadamente y al aplicar un principio apropiado, las
ecuaciones diferenciales parciales que gobiernan el problema pueden ser aproximadas por un
conjunto de ecuaciones algebraicas a través de la malla. El sistema de ecuaciones algebraicas
para la totalidad del dominio del problema puede ser formado al ensamblar el conjunto de
las ecuaciones algebraicas para todas las mallas.
En los métodos sin mallas MFree (mesh free) se establece un sistema de ecuaciones algebraicas
para la totalidad del dominio del problema sin el uso de la malla predefinida. Los métodos
MFree usan un conjunto de nodos diseminado o esparcidos por el interior del dominio
del problema, aśı como también un conjunto de nodos diseminados por las fronteras para
2
representar (no discretizar) tanto el dominio como sus fronteras. Estos conjuntos de puntos
diseminados no forman una malla, lo cual significa que ninguna información de la relación
entre ellos es requerida, al menos para el campo de interpolación.
Existen muchos métodos Mfree, tales como el método Element Free Galerkin (EFG) [5], el
método sin malla local de Petrov-Galerkin (MLPG) [3], el método de interpolación puntual
(PIM) [8], el método Smoothed Particle Hydrodynamics [16].
Esta investigación se enfoca en el estudio de los métodos Smoothed Particles Hydrodynamics
(SPH) [15], Moving Least-Squares (MLS),( [11], [14], [8]), Least-Squares (LS) y Weighted
Least Squares (WLS) [22] y su aplicación en los modelos deformables basados en la teoŕıa
de elasticidad lineal y la simulación dinámica de los fluidos.
Esta trabajo combina y desarrolla las ĺıneas de investigación de métodos numéricos y
computación gráfica de la Facultad de Ciencias y Tecnoloǵıa de la Universidad de Carabobo
mediante el análisis y desarrollo de los métodos de aproximación sin mallas y sus aplicaciones
en la resolución numérica de ecuaciones diferenciales.
3
Caṕıtulo 2
Métodos sin Mallas
Los denominados métodos sin malla han atráıdo muy recientemente el interés de los investi-
gadores, debido a su flexibilidad para resolver problemas prácticos de simulación numérica.
El principal objetivo de los métodos sin malla es superar las dificultades que aparecen en
los problemas de simulación numérica al tener que remallar los dominios en estudio, ya que
en estos métodos es suficiente con añadir nodos donde sea necesario. Se trata por tanto de
enfoques novedosos en el área de computación gráfica que han cambiado la estrategias a la
hora de realizar simulación o animación de objetos deformables.
Los métodos sin mallas son aquellos donde el objeto es modelado por medio de una nu-
be de part́ıculas sin conectividad previa. Existe una diversidad de métodos sin mallas y
los más destacados son: SPH (Smoothed Particle Hydrodynamics), DEM (Method Element
Difuse),EFG (Element Free Galerkin), Mı́nimos Cuadrados (Moving Least-Squares), LBIE
(Local Boundary Integral Equation), PUM (Partition of Unity Method).
El método de SPH fue desarrollado originalmente en 1977 para la simulación de gas interes-
telar en el campo de la astro f́ısica. Fue desarrollado por Bob Gingold y Joe Monaghan [9] e
independientemente por Leon Lucy [16]. El método SPH es un método poderoso que permi-
te la solución de problemas dinámicos complejos tanto de fluidos como de otros materiales.
Desde su desarrollo ha sido usado exitosamente en mútiples áreas de estudio tales como:
rompimiento de oleaje, flujo de lava, estudio de fracturas elásticas, problemas de superficie
libre y materia granular.
4
La aproximación de los mı́nimos cuadrados (MLS - Moving Least-Squares)fue desarrollada
por Lancaster y Salkauskas en 1981 [12], quienes fueron los primeros que usaron la
aproximación de MLS para el suavizado e interpolación de datos en el área de computación
gráfica.
2.1. Smoothed Particle Hydrodynamics - SPH
La formulación de SPH se divide en dos etapas. La primera es la representación integral
o también llamada aproximación por núcleos del campo de funciones. La segunda es la
aproximación de part́ıculas.
En la primera, la integración de una función arbitraria y una función de núcleo (también
llamada función de peso) dá la función de núcleo en la forma de la representación integral de
la función. La representación integral de la función es posteriormente aproximada agregando
los valores de las part́ıculas vecinas más próximas las cuales producen la aproximación de
part́ıculas de la función a un punto discreto o part́ıcula.
Aqúı, el término representación integral o aproximación por núcleos y aproximación de
part́ıculas serán usados para lo mismo.
El Concepto de representación integral de una función f(x) usada en el método SPH comienza
de la siguiente identidad
f(x) =
∫
Ω
f(y)δ(x− y)dy (2.1)
donde f es una función del vector x y δ(x− y) es la función delta de Dirac dada por
δ(x− y) =
{
1 x = y
0 x 6= y (2.2)
En la ecuación (2.1), Ω es el volumen de la integral que contiene a x. La ecuación (2.1)
implica que una función puede ser representada en una forma integral. Puesto que la función
delta de Dirac es usada, la representación integral en la ecuación 2.1 es exacta siempre y
cuando f(x) esté definida y sea continua en Ω.
5
Si la función núcleo Delta δ(x − y) es reemplazada por una función suave W (x − y, h), la
representación integral de f(x) viene dada por
f(x)
.
=
∫
Ω
f(y)W (x− y, h)dy (2.3)
Donde W es la llamada función núcleo, o función de peso. h es la longitud de suavidad que
define el área de influencia de la función de peso W . Notar que aún cuando W no es la función
de Dirac, la representación integral de la ecuación (2.3) es solamente una aproximación. Este
es el origen del término aproximación por núcleos.
La función W es elegida de modo que cumpla las condiciones siguientes:
1. La función de peso debe estar normalizada sobre su dominio∫
Ω
W (x− y, h)dy = 1 (2.4)
2. La función de peso será de soporte compacto
W (x− y, h) = 0, | x− y |> κh (2.5)
El tamaño del soporte compacto está definido por la longitud del núcleo multiplicado
por un factor κ, donde h es la longitud del núcleo y κ determina la propagación de la
función de peso dada. | x− y |≤ κh define el dominio de soporte de la part́ıcula en el
punto x.
3. W (x − y, h) ≥ 0 para todo punto y dentro del dominio de soporte para el punto x
(Positividad).
4. La función de núcleo será una función monótona decreciente a medida que se incremente
la distancia entre part́ıculas.
5. La función de núcleo debeŕıa satisfacer la condición de la función delta de Dirac a
medida que la longitud de núcleo tienda a cero (propiedad de la función Delta)
ĺım
h→0
W (x− y, h) = δ(x− y). (2.6)
6. La función de núcleo debe ser una función par (propiedad de simetŕıa)
7. La función de núcleo debe ser lo suficientemente suave (condición de suavidad )
6
La primera propiedad de normalización asegura que la integral de la función de núcleo
sobre el dominio de soporte sea la unidad. También asegura el orden de consistencia de la
representación integral de una función continua.
La segunda propiedad transforma una aproximación SPH de forma global a una de
operaciones locales.
La tercera propiedad establece que la función núcleo seŕıa no negativa en el dominio de
soporte. Esto no es matemáticamente necesario como un requerimiento de convergencia,
pero importante para asegurar un sentido f́ısico, o estable, de la representación de algunos
fenómenos f́ısicos. Algunas funciones de núcleo usadas en la práctica son negativas en partes
del dominio soporte. Sin embargo en simulaciones hidrodinámicas, valores negativos para la
función de núcleo pueden tener serias consecuencias que resultan en algunos parámetros no
f́ısicos tales como densidad negativa y enerǵıa.
La cuarta propiedad está basada sobre consideraciones f́ısicas en que las part́ıculas más
próximas tendŕıan una influencia mayor sobre la part́ıcula en cuestión. En otras palabras, con
el incremento de la distancia entre dos part́ıculas que interactúan, las fuerzas de interacción
decrecen.
La quinta propiedad asegura que cuando la longitud de núcleo tiende a cero, los valores
aproximados se acercan a los valores reales de la función.
La sexta propiedad significa que las part́ıculas a una misma distancia pero en diferentes
posiciones tienen el mismo efecto sobre una part́ıcula dada. Esta no es una condición
muy ŕıgida, y algunas veces es violada en algunos métodos sin mallas que tienen mayor
consistencia.
La séptima propiedad tiene como objetivo obtener una mejor aproximación. Para el
aproximaciones de una función y sus derivadas, la función de núcleo necesita ser lo
suficientemente continua para obtener buenos resultados.
Cualquier función que tenga las propiedades anteriores puede ser empleado como una función
de peso o función núcleo del método SPH.
En el art́ıculo original del método SPH, Lucy [16] usó la siguiente función núcleo
W (x− y, h) = W (R, h) = αd
{
(1 + 3R)(1−R)3 R ≤ 1
0 R > 1
(2.7)
7
donde αd = 5/(4h), αd = 5/(πh
2) y αd = 105/(16h
3) en una, dos y tres dimensiones
respectivamente. r es la distancia entre los puntos x y y. R es la distancia relativa entre los
puntos x y y y la longitud de núcleo h definida por R = r
h
= |x−y|
h
.
Monagan [18] señaló que para encontrar una interpretación f́ısica de una ecuación SPH,
siempre es mejor asumir que la función de núcleo sea una Gaussiana. Gingold y Monaghan [9]
en su art́ıculo original seleccionaron el siguiente núcleo Gaussiano para simular estrellas no
esféricas
W (R, h) = αde
−R2 (2.8)
donde αd = 1/(π
1/2), αd = 1/(π
3/2h2) y αd = 1/(π
3/2h3) en una, dos y tres dimensiones
respectivamente.
El núcleo Gaussiano es suficientemente suave incluso para derivadas de orden grandes y es
considerado como una buena selección ya que es muy estable y exacta, especialmente para
part́ıculas desordenadas. Sin embargo, no es compacta porque nunca se va a cero a menos que
R se aproxime al infinito. Pero, debido a que se aproxima a cero muy rápido numéricamente
es practicamente compacta
2.2. Mı́nimos Cuadrados Móviles - MLS
La aproximación de Mı́nimos Cuadrados Móviles fue inicialmente trabajada sobre el ajuste
de datos y construcción de superficies (Lancaster y Salkauskas, [12]). La invención de la
aproximación MLS fue la clave para el desarrollo de muchos métodos sin mallas, porque
los MLS pueden dar una aproximación continua para una función de campo sobre todo el
dominio del problema. Este enfoque, en la actualidad, es extensamente usado en muchos
métodos sin mallas para construir las funciones de forma. Nayroles [21] usó la aproximación
MLS la primera vez para desarrollar los llamados métodos difusos (DEM). Muchos métodos
sin mallas basados en la aproximación MLS han sido desarrollados desde entonces, tales
como los elementos libres de Galerkin [5] y el método local sin malla Petrov-Galerkin [3].
El procedimiento del MLS es el siguiente. Sea u(x) la función del campo de variables definida
en el dominio Ω. La aproximación de u(x) en el punto x es denotada por uh(x). Aśı, la
aproximación MSL escrita en términos de campo es:
8
u(x) =
m∑
j
pj(x)aj(x) = p
T (x)a(x) (2.9)
donde m es el número de términos de los monomios (base polinomial), y a(x) es un vector
de coeficientes dado por
aT (x) = {a0(x), a1(x), · · · , am(x)} (2.10)
los cuales son funciones de x.
En la ecuación (2.9), p(x) es un vector de las funciones bases que consiste usualmente de
monomios de orden menor o igual a m. En 1D, una basecompleta de polinomios de orden
m está dado por
pT (x) = {p0(x), p1(x), · · · , pm(x)} = {1, x, · · · , xm} (2.11)
El vector de coeficientes de a(x) en la ecuación (2.9) es determinado usando los valores de
la función en el conjunto de los nodos que están incluidos en el dominio de soporte de x. Un
dominio de soporte de un punto x determina el número de nodos que son usados localmente
para aproximar el valor de la función en x.
Aśı, dada un conjunto de n nodos para la función de campo u1, u2, · · · , un en n nodos
x1, x2, · · · , xn que están en el dominio de soporte, entonces la ecuación (2.9) es usada para
calcular los valores aproximados a la función de campo en estos nodos:
uh(x, xi) = p
T (xi)a(x), i = 1, 2, · · · , n (2.12)
Notar que aqúı a(x) es una función arbitraria de x. Un funcional de residual ponderado es
construido usando los valores aproximados de la función de campo y los parámetros nodales
ui = u(xi)
J =
n∑
i
Ŵ (x− xi)[uh(x, xi)− u(xi)]2
=
n∑
i
Ŵ (x− xi)[pT (xi)a(x)− u(xi)]2
(2.13)
donde Ŵ (x−xi) es una función de peso y ui es el parámetro nodal de las variables de campo
en el nodo i. Estas funciones de pesos proveen una ponderación para los residuos en los
diferentes nodos en el dominio de soporte. Se espera que los nodos que están muy alejados
9
de x tengan un peso menor que los que están más próximos a x. También esta función de peso
debe asegurar que los nodos que entren o salgan de dominio de soporte de x sea en forma
suave. La condición de la unidad que es requerida para el SPH no es necesario aqúı para la
aproximación MLS.
En la aproximación MLS, en un punto arbitrario x, a(x) es elegido para minimizar el residual
ponderado. La condición de minimización requiere
∂J
∂a
= 0 (2.14)
lo cual resulta en el siguiente sistema de ecuaciones lineales
A(x)a(x) = B(x)Us (2.15)
Donde A es la matriz de momentos ponderada, la cual está dada por
A(x) =
n∑
i
Ŵi(x)p(xi)p
T
xi
(2.16)
donde
Ŵi(x) = Ŵ (x− xi) (2.17)
En la ecuación (2.15) la matriz B tiene la forma
B(x) = [B1,B2, . . . ,Bn] (2.18)
Bi = Ŵi(x)p(xi) (2.19)
y Us es en vector que reúne la todos los parámetros nodales de las variables de campo para
todos los nodos en el dominio de soporte:
Us = {u1, u2, · · · , un}T (2.20)
10
Resolviendo la ecuación (2.15) para a(x), obtenemos
a(x) = A(x)−1B(x)Us (2.21)
y sustituyendo este valor en la ecuación (2.12) se tiene que
uh(x, xi) =
n∑
i
m∑
j
pj(x)(A(x)
−1B(x))jiui (2.22)
o
uh(x, xi) =
n∑
i
φi(x)ui (2.23)
donde las funciones φi(x) están definidas por
φi(x) =
m∑
j
pj(x)(A(x)
−1B(x))ji = p
TA−1Bi (2.24)
Notar que m es el número de términos de la base de polinomios p(x), el cuál es mucho
menor que n, donde n es el número de nodos usados en dominio de soporte para construir
las funciones de formas. El requerimiento de que n � m evita la singularidad de la matriz
de momento ponderada, por lo tanto A−1 existe. En los siguientes cápitulos se dan detalles
de los métodos Least-Squares (LS) y Weighted Least Squares(WLS) [22].
11
Caṕıtulo 3
Aproximación de nubes de puntos
usando métodos sin mallas
3.1. Mı́nimos Cuadrados (LS)
Los Mı́nimos cuadrados es una técnica de aproximación numérica de una función u(x) a
partir de su valor en ciertos puntos {xi . . .xn} pertenecientes a algún conjunto Ω ⊂ Rd donde
i ∈ {1, . . . , N}. Se quiere obtener una función definida globalmente f(x) que se aproxime
a los valores escalares dados fi en los puntos xi en el sentido de mı́nimos cuadrados con el
error funcional
JLS =
∑
i
‖f(xi)− fi‖2
Por lo tanto, se plantea el problema de minimización siguiente
mı́n
f∈
∏d
m
∑
i
‖f(xi)− fi‖2, (3.1)
donde f es tomada de
∏d
m, el espacio de polinomios de grado total m en d dimensiones
espaciales, y puede ser escrito como
f(x) = b(x)T c = b(x) · c, (3.2)
donde b(x) = [b1(x), . . . , bk(x)]
T es el vector de base polinomial y c = [c1, . . . , ck]
T es el vector
de coeficientes desconocidos que se desean miminizar en la ecuación (3.1). Algunos ejemplos
12
de bases polinomiales son (a) para m = 2 y d = 2, b(x) = [1, x, y, x2, xy, y2]T , (b) para un
ajuste lineal en R3(m = 1, d = 3), b(x) = [1, x, y, z]T , y (c) para el montaje de una constante
en dimensiones arbitrarias b(x) = [1]. En general, el número k de elementos en b(x) (y por
tanto en c) es dado por k = (d+m)!
m!d!
Solución: Se puede minimizar (3.1) al establecer las derivadas parciales del error funcional
JLS a cero, i.e, ∇JLS = 0 donde ∇ = [ ∂∂c1 , . . . ,
∂
∂ck
]T que es una condición necesaria para un
mı́nimo. Al tomar derivadas parciales con respecto a los coeficientes desconocidos c1, . . . , ck,
se obtiene un sistema de ecuaciones lineales de donde se puede calcular c [22]
∂JLS
∂c1
= 0 :
∑
i
2b1(xi)[b(xi)
T c− fi] = 0
∂JLS
∂c2
= 0 :
∑
i
2b2(xi)[b(xi)
T c− fi] = 0
...
∂JLS
∂c3
= 0 :
∑
i
2b3(xi)[b(xi)
T c− fi] = 0
(3.3)
En notatición Matriz-Vector esto se escribe como:
∑
i
2b(xi)[b(xi)
T c− fi] = 2
∑
i
[b(xi)b(xi)
T c− b(xi)fi] = 0 (3.4)
Dividiendo por la constante y reordenando se obtiene:
∑
i
b(xi)b(xi)
T c =
∑
i
b(xi)fi, (3.5)
Que se resuelve como
c =
[∑
i
b(xi)b(xi)
T
]−1∑
i
b(xi)fi, (3.6)
13
Si la matriz cuadrada ALS =
∑
i b(xi)b(xi)
T es no singular, i.e det(ALS) 6= 0, sustituyendo
la ecuación (3.6) en la ecuación (3.2) se obtiene la función de ajuste f(x). Para k pequeños
(k < 5) la inversión de la matriz en la ecuación (3.6) puede llevarse a cabo de forma expĺıcita,
de lo contrario métodos numéricos son la herramienta preferida.
Ejemplo. Dados los puntos en R2 y se desea ajustar con una cuadrática, polinomio bivariado,
i.e. d = 2,m = 2 y por lo tanto b(x) = [1, x, y, x2, xy, y2]T entonces el Mı́nimo Cuadrado
resultante se ve aśı:
∑
i

1 xi yi x
2
i xiyi y
2
i
xi x
2
i xiyi x
3
i x
2
i yi xiy
2
i
yi xiyi y
2
i x
2
i yi xiy
2
i y
3
i
x2i x
3
i x
2
i yi x
4
i x
3
i yi x
2
i y
2
i
xiyi x
2
i yi xiy
2
i x
3
i yi x
2
i y
2
i xiy
3
i
y2i xiy
2
i y
3
i x
2
i y
2
i xiy
3
i y
4
i


c1
c2
c3
c4
c5
c6
 =
∑
i

1
xi
yi
x2i
xiyi
y2i
 fi. (3.7)
3.2. Mı́nimos Cuadrados Ponderados (WLS)
Formulación del Problema. En la formulación de Mı́nimos Cuadrados Ponderados
(Weighted Least Squares) se utiliza el error funcional JWLS =
∑
i θ(‖x̄ − xi‖)‖f(xi) − fi‖2
para un punto estimado x̄ ⊂ Rd que se minimiza
mı́n
f∈
∏d
m
∑
i
θ(‖x̄− xi‖)‖f(xi)− fi‖2, (3.8)
es similar a la ecuación (3.1) sólo que ahora el error es ponderado por θ(d) donde di son las
distancias euclidianas entre x̄ y las posiciones de los puntos de la data xi.
Los coeficientes desconocidos que se deben obtener de la solución a la ecuación (3.8) son
ponderados por las distancias a x̄ y por lo tanto una función x̄.
Aśı, la local, aproximación de mı́nimos cuadrados ponderados en x̄ se escribe como
fx̄(x) = b(x)
T c(x̄) = b(x) · c(x̄), (3.9)
14
y sólo define localmente dentro de una distancia R alrededor de x̄, i.e. ‖x− x̄‖ < R.
Función de Peso. Muchas opciones para la función de ponderación θ se han propuesto en
la literatura, tal como una gaussiana, [4]
θ(d) = e−
d2
h2 (3.10)
donde h es un parámetro de espaciado que se puede utilizar para suavizar los rasgos pequeños
en los datos. Otra función de ponderación popular, [15], con soporte compacto es la función
siguiente
θ(d) =
(
1− d
h
)4(
4d
h+ 1
)
(3.11)
Esta función está bien definida en el intervalo d ∈ [h, 0] , y, además θ(0) = 1, θ(h) = 0, θ′(h) =
0, θ′′(h) = 0. Varios autores, [12], [4], [19], sugieren usar funciones de peso de la siguiente
forma:
θ(d) =
(
1
d2 + ε2
)
(3.12)
Es de hacer notar que establecer el parámetro ε a cero da como resultado una singularidad
en d = 0, lo que obliga a la función de ajuste de la MLS a interpolar los datos.
Solución. Análogo al método de Mı́nimos Cuadrados convencional, se toman derivadas
parciales del error funcional JWLS con respecto a los coeficientes desconocidos c(x̄)15
∑
i
θ(di)2b(xi)[b(xi)
T c(x̄)− fi] =
2
∑
i
[θ(di)b(xi)b(xi)
T c(x̄)− θ(di)b(xi)fi] =0,
(3.13)
donde di = ‖x̄− xi‖. Dividiendo entre la constante y reordenando se tiene
∑
i
θ(di)b(xi)b(xi)
T c(x̄) =
∑
i
θ(di)b(xi)fi, (3.14)
y despejando los coeficientes
c(x̄) = [
∑
i
θ(di)b(xi)b(xi)
T ]−1
∑
i
θ(di)b(xi)fi], (3.15)
La única diferencia entre las ecuaciones (3.6) y (3.15) son los términos de ponderación. Ob-
serve de nuevo, sin embargo, que mientras los coeficientes c en la ecuación (3.6) son globales,
los coeficientes c(x̄) son locales y tienen que ser recalculados para cada x̄. Si la matriz cua-
drada AWLS =
∑
i θ(di)b(xi)b(xi)
T (a menudo se denomina Matriz Momento) es no singular
(i.e. det(AWLS) 6= 0), sustituyendo la ecuación (3.15) en la ecuación (3.9) se obtiene la fun-
ción de ajuste fx̄(x).
Aproximación Global usando una Partición de la Unidad PU. Ajustando polinomios
a j ∈ {1, . . . , n} los puntos estimados y discretos x̄j en el parámetro de dominio Ω se
pueden ensamblar en una aproximación global a la data garantizando que cada punto en Ω
está cubierto por al menos un polinomio de aproximación, es decir, el soporte de las funciones
de peso θj centrado en los puntos x̄j que cubren Ω
Ω = ∪jsupp(θj). (3.16)
La ponderación adecuada de estas aproximaciones se puede lograr mediante la construcción
de una partición de la unidad (PU) desde θj
ϕj(x) =
θj(x)∑n
k=1 θk(x)
(3.17)
16
donde
∑
j ϕi(x) ≡ 1 en todo Ω. La aproximación global viene a ser
f(x) =
∑
j
ϕj(x)b(x)
T c(x̄j). (3.18)
Una cuestión numérica. Para evitar inestabilidades numéricas con posibles números
grandes en AWLS puede ser beneficioso llevar a cabo el procedimiento de ajuste en un sistema
de coordenadas local con respecto a x̄, es decir, para cambiar x̄ en el origen. Por lo tanto, se
reescribe la función de ajuste local en x̄ como
fx̄(x) = b(x− x̄)T c(x̄) = b(x− x̄) · c(x̄), (3.19)
los coeficientes asociados como
c(x̄) =
[∑
i
θ(di)b(xi − x̄)b(xi − x̄)T
]−1∑
i
θ(di)b(xi − x̄)fi, (3.20)
y la aproximación global como
f(x) =
∑
j
ϕj(x)b(x− x̄j)T c(x̄j). (3.21)
3.3. Mı́nimos Cuadrados móviles (MLS)
Método. El Moving Least Squares fue propuesto por Lancaster and Salkauskas [12] para
suavizar e interpolar datos. Esta técnica de aproximación es una variante de la aproximación
por mı́nimos cuadrados ponderados, con la particularidad de que ahora la función de peso
se traslada de un punto a otro, de manera que el máximo está en el punto x donde se quiere
obtener el valor aproximado, por ello el nombre de mı́nimos cuadrados móviles [2]. Esto se
lleva a cabo fácilmente si se asocian unos pesos positivos Wi (los cuales van a depender de
x ) a cada punto y se minimiza.
17
La idea es iniciar con una formulación de mı́nimos cuadrados ponderados para un punto
fijo arbitrario en Rd, a continuación mover este punto sobre el dominio entero, donde se
calcula un mı́nimo cuadrado ponderado y se evalua para cada punto individualmente. Se
puede demostrar que la función global f(x), obtenida a partir de un conjunto de funciones
locales
f(x) = fx(x), mı́n
f∈
∏d
m
∑
i
θ(‖x− xi‖)‖fx(xi)− fi‖2, (3.22)
es continuamente diferenciable si y sólo si la función de ponderación es continuamente
diferenciable.
Además θ(‖x − xi‖) es positiva, grande para los puntos x cercanos a xi, y relativamente
pequeña para los x más distantes. diferencia θ(‖x−xi‖) decrece monótonamente al aumentar
la diferencia ‖x− xi‖ hasta valer cero (soporte compacto).
Aśı que en lugar de construir la aproximación global utilizando la ecuación (3.18), se utilizan
las ecuaciones (3.9) y (3.15) (o las ecuaciones (3.19) y (3.20)) y se construye y evalúa un
polinomio local adaptado continuamente a lo largo de todo el dominio Ω, lo que resulta en
la función de ajuste MLS.
Para el caso general, cuando el sistema es de r dimensiones x = [x1, x2, ..., xr] y se quiere
aproximar una función f(x), el funcional J pasa a ser ahora de la siguiente manera:
J =
∑
j
wj(x− xj)(P T (xj)c(x)− fj)2 (3.23)
y en forma matricial
J = (Pc(x)− f)TW (x)(Pc(x)− f) (3.24)
Donde f es el vector de nx1 valores nodales y P la matriz de nxm siendo n el número de
puntos y m el número de términos del polinomio aproximador (base polinomial), y W (x) es
la matriz diagonal de la función de ponderación
W = diag[w(x− xi), ..., w(x− xn)] (3.25)
18
derivando se obtiene
∂J
∂c
= A(x)c(x)−B(x)f = 0 (3.26)
donde
A(x) =P TW (x)P
B(x) =P TW (x)
(3.27)
y por lo tanto
c(x) = A(x)−1B(x)f (3.28)
en esta ocación los coeficientes de la matriz c(x), dependen de x, debido a la inclusión de la
función de ponderación en el funcional, de forma que la función de aproximación queda
f(x) = p(x)TA(x)−1B(x)f (3.29)
donde p(x) es la base polinomial evaluada en cada x y y perteneciente al dominio discretizado,
p(x) = [1, x, y, x2, xy, y2]T
Utilizando la ecuación (3.12) como la función de ponderación con un ε muy pequeño asigna
pesos cercanos al infinito cerca de los puntos de entrada de datos, obligando a la función
de ajuste MLS a interpolar los valores de la función prescrita en estos puntos. Por lo tanto,
mediante la variación de ε se puede influir directamente en la aproximacion / interpolación
de la naturaleza de la función de ajuste MLS. O también se puede hacer uso de la función de
ponderación exponencial (3.1012) con un radio de soporte entre 0 y 1 con la cual se pueden
obtener excelentes resultados.
3.3.1. Aplicaciones de Mı́nimos Cuadrados Móviles
Diseñada originalmente para el tratamiento de datos y la generación de superficies, la ténica
de Mı́nimos Cuadrados Móviles se ha vuelto muy popular entre los grupos de investigación
que trabajan con métodos sin malla, siendo ampliamente utilizada en aplicaciones tanto
eulerianas como lagrangianas. Esta clase de métodos de aproximación es particularmente
19
competitiva en la reconstrucción de una determinada función y sus derivadas sucesivas a
partir de los valores de dicha función en una serie de puntos dispersos. Se ha pretendido
que el empleo de Mı́nimos Cuadrados Móviles proporcione algo parecido a unas funciones de
formas para esquemas de volúmenes finitos en mallas no estructuradas, representando una
interesante alternativa a los métodos existentes en la actualidad, aśı como también hacer
uso del método en aplicaciones en problemas de flujo de fluidos comprensibles. En general,
el MLS es un método bastante útil a la hora de reconstruir superficies a partir de puntos [1],
simulaciones [4], animaciones [19], interpolación y aproximación de superficies impĺıcitas [24]
3.4. Implementación y resultados
Como se observó en la sección teórica, los métodos sin mallas MLS requieren más cálculos
matriciales que toman gran cantidad de tiempo de ejecución. Sin embargo, dado los estudios
realizados a la programación de las tarjetas gráficas se han conseguido buenos resultados en
tiempo de ejecución y en obtención de la solución.
Las primeras pruebas se realizaron con pocos puntos en la data. Se pudo notar que los
métodos LS, WLS son aproximaciones ”gruesas 2la solución no se ajusta a los datos, mientras
tanto, las soluciones obtenidas con los métodos SPH y MLS dan mejores resultados.
los resultados obtenidos con los métodos sin mallas, en este caso de 9 puntos de entrada, sobre
la CPU y sobre el dispositivo gráfico se lograron en el mismo tiempo aproximadamente ya
que como se observó en la sección de implementación de los algoritmos del álgebra matricial,
en principio, para matrices no muy grandes la multiplicación de matrices es bastante rápida
en la CPU, el problema viene a ser con matrices más grandes.
Se realizaron las pruebas usando una base polinomial b(x) = [1, x, y, x2, xy, y2]T .
Pi = {(1, 1), (1,−1), (−1, 1), (−1,−1), (0, 0), (1, 0), (−1, 0), (0, 1), (0,−1)}
con el conjunto de valores de función asociados
fi = {1.0,−1.0, 0.0, 0.0, 1.0, 0.0,−1.0,−1.0, 1.0}.
En la figura (3.1) se muestra la soluciónde la aproximación a los datos usando el Método de
Mı́nimo Cuadrado(LS). Se aprecia que existen datos muy alejados de la superficie solución.
20
Figura 3.1: Aproximación por Mı́nimos Cuadrados
En la figura 3.2 se muestra la solución al problema de ajuste de datos usando el Método
de Mı́nimos Cuadrado ponderado (WLS). Aqúı, el punto x̄ fijo es el origen. Con lo cual,
la función de peso que realiza la ponderación está centrada de este punto. Por lo tanto la
solución tiene un ajuste óptimo para los datos próximos a este punto fijo, fenómeno que se
aprecia en esta imagen.
Figura 3.2: Aproximación por Mı́nimos Cuadrados ponderados
La imagen (3.3) muestra la solución de ajuste de datos usando el Método de Mı́nimo
Cuadrado Móvil (MLS). De los métodos estudiados éste es el que mejor se ajusta a los
datos, pero el tiempo de ejecución es mayor. En esta imagen se aprecia un ajuste óptimo a
los datos usando esta estrategia.
21
Figura 3.3: Aproximación por Mı́nimos Cuadrados Móviles
En la figura (3.4)se presenta la solución del ajuste de datos usando el método SPH.
Aqúı usamos la discretización de la formulación integral y observamos que la solución
obtenida con este enfoque es óptimo para el ajuste de datos. El tiempo de ejecución usando
SPH es menor que el empleado con el método MLS.
Figura 3.4: Aproximación por SPH
También se realizaron pruebas con otro tipo de datos de entrada, aqúı mostramos un
ejemplo con 256 datos. El método utilizado en este caso fue el método sin mallas Mı́nimos
Cuadrados Móviles. En este ejemplo, la nube de puntos fue generada de forma aleatoria
y el dominio de definición fue generada con 10.201. El tiempo de ejecución fue de 44,43
segundos. La programación de estos métodos fue en lenguaje C, en un equipo de escritorio
22
con un procesador Intel Core i7 920 a 2.66 Ghz. Una GPU nVidia geForce 9500GT y una
memoria RAM de 4GB DDR3.
Figura 3.5: MLS con 256 datos
23
Caṕıtulo 4
Modelos deformables usando MLS
Este caṕıtulo contiene bases teóricas relacionadas con las ecuaciones de elasticidad lineal, la
aproximación del operador de deformación utilizando el método de los mı́nimos cuadrados,
la manera en que son obtenidas las fuerzas, las distintas funciones de peso utilizadas en esta
investigación y el cómo se realiza el movimiento de las part́ıculas del cuerpo en deformación.
4.1. Teoŕıa de Elasticidad Lineal
A continuación se muestra el modelo f́ısico-matemático que surge de la teoŕıa de elasticidad
lineal y cuyas leyes son las que gobiernan el modelo deformable en estudio, [7] . Considérese
un objeto 3D cuyas coordenadas de material son: ~Xi =
[
xi, yi, zi
]T
. Para describir
la deformación de dicho objeto se utiliza un campo vectorial de desplazamiento ~Ui =[
ui, vi, wi
]T
, donde ui = u(xi, yi, zi), vi = v(xi, yi, zi) y wi = w(xi, yi, zi) son funciones
de las coordenadas de material. En la figura 4.1 se observa que las coordenadas de una
part́ıcula originalmente se encuentran en ~Xi y en su posición deformada son ~Xi + ~Ui. Se
denota por ~Xij el vector diferencia entre los vectores posición de las part́ıculas i y j.
~Xij = ~Xj − ~Xi (4.1)
En el análisis de elasticidad lineal, el Jacobiano (J) de una transformación es fundamental
para determinar ciertas variaciones del objeto en estudio. Éste se expresa para cada part́ıcula
i como:
Ji = ∇ ~Ui
T
+ I (4.2)
24
donde I es la matriz identidad 3x3. Por comodidad en la escritura, el Jacobiano se escribe
de la siguiente manera:
Ji =
 JuTiJvTi
JwTi
 (4.3)
Para medir el grado de deformación (strain) de una part́ıcula i se empleó el tensor de Green-
Saint-Venant el cual está definido por la siguiente matriz:
εi = (J
T
i Ji)− I (4.4)
De la teoŕıa de elasticidad lineal se tiene que el grado de deformación ε y la tensión (stress)
σ están linealmente relacionados, es decir:
σi = Cεi (4.5)
donde C es un tensor de rango cuatro relacionado con la ley constitutiva de los materiales
y tanto ε como σ son tensores simétricos. Para materiales isotrópicos C tiene sólo dos
coeficientes independientes llamados módulo de Young (E) y tasa de Poisson (v). De esta
teoŕıa se deduce que el tensor de tensión es:
 σi(1,1)σi(2,2)
σi(3,3)
 = E
(1 + v)(1− 2v)
 1− v v vv 1− v v
v v 1− v
 εi(1,1)εi(2,2)
εi(3,3)
 (4.6)
σi(1,2) = σi(2,1) =
E
1+v
εi(1,2)
σi(1,3) = σi(3,1) =
E
1+v
εi(1,3)
σi(2,3) = σi(3,2) =
E
1+v
εi(2,3)
(4.7)
Para calcular tensión(stress), deformación (strain) y fuerzas elásticas es necesario determinar
el campo de desplazamiento ∇~U . Para ello se hace uso del método de mallas libres MLS
siguiendo las pautas del trabajo de Müller [20]. En este caso, las derivadas en una part́ıcula
i son estimadas a partir de sus vecinos.
4.2. Métodos de Mallas Libres MLS
Los métodos de mallas libres son métodos numéricos desarrollados para aproximar las
ecuaciones diferenciales parciales sobre una nube de part́ıculas sin conectividad previa [8].
Uno de los primeros enfoques utilizados para resolver estas ecuaciones es el método SPH
(Smoothed Particle Hydrodynamics), el cual fue empleado en computación gráfica en el año
1977 por Lucy [16]. Luego Monaghan [17], en 1982, utilizó las nociones de aproximación para
núcleos.
25
Lancaster y Salkauskas [11] en 1981, fueron los primeros que usaron la aproximación de
MLS en el área de computación gráfica, para el suavizado e interpolación de datos. Si una
función ~U ~X definida sobre un dominio Ω ⊂ Rd es suficientemente suave, se puede definir una
aproximación local alrededor de un punto fijo ~̄X ∈ Ω̄ como:
~U( ~X, ~̄X) ≈ L ~̄XU( ~X) = p
T ( ~X)a( ~̄X), (4.8)
donde
~U( ~X, ~̄X) =
{
~U( ~X) ∀ ~X ∈ Ω̄, || ~X − ~̄X||2 < ρ
0 en otro caso
(4.9)
y L ~̄X es un operador. El operador || ||2 es la norma eucĺıdea. El vector p( ~X) contiene una
base completa en dimensión d con consistencia n. En vista de que la aproximación local usa
la mejor aproximación de ~U en un cierto sentido de mı́nimos cuadrados, el vector incógnita
a( ~̄X) es seleccionado para minimizar la norma 2 de los mı́nimos cuadrados. De esta forma,
el vector coeficiente a( ~̄X) es seleccionado para satisfacer la condición H ~̄X(a(
~̄X)) ≤ H ~̄X(b)
para todo b 6= a ∈ R con
H ~̄X(a) =
N∑
i=1
W ( ~X − ~Xi)[L ~̄XU( ~X)− U( ~X,
~̄X)]2 =
N∑
i=1
W ( ~X − ~Xi)
(
pT ( ~Xi)a(
~̄X)− ~Ui
)2
(4.10)
~Xi se refiere a la posición de losN nodos dentro del dominio. Cada función de pesoWi( ~X− ~Xi)
está definida en soportes Ω̃i alrededor de cada nodo ~Xi para garantizar la aproximación local.
La función de peso Wi puede ser también elegida individualmente para cada nodo.
4.2.1. Aproximación de ∇~U empleando el método de MLS
El cálculo de las derivadas parciales ∇ui,∇vi,∇wi del vector desplazamiento ~Ui de la
part́ıcula i, siguiendo el método de mallas libres MLS, se muestran a continuación.
Haciendo uso de la expansión de Taylor en cada una de las coordenadas de ~Ui, midiendo el
error de aproximación promediado con el núcleo polinomial de soporte compacto, se tiene
que las derivadas que minimizan este error satisfacen para ∇ui lo siguiente:
26
∑
j
~Xij ~X
T
ijWij∇ui =
∑
j
(uj − ui) ~XijWij (4.11)
La expresión
Ai =
∑
j
(
~Xij ~X
T
ijWij
)
, (4.12)
es conocida como la matriz de momentos.
Mediante estas ecuaciones:
∇ui = A−1i
∑
j
(
(uj − ui) ~XijWij
)
(4.13)
De forma similar se obtiene un resultado para ∇vi,∇wi. Se deduce entonces que ∇ ~Ui
T
es:
∇ ~Ui
T
=
 ∇uTi∇vTi
∇wTi
 (4.14)
Para calcular la matriz A−1 se utilizó el algoritmo de descomposición en valores singulares
[25], (SVD - Singular Values Decomposition).
4.2.2. Funciones de Peso
Se considera en este modelo que la masa de una part́ıcula es fija y se calcula como sigue:
mi = s r̄i
3 ρ (4.15)
donde s es un factor de escala igual para todas las part́ıculas, ρ es la densidad del material
y r̄i es la distancia promedio de las k part́ıculas más cercanas a la part́ıcula i. La masa mi
está distribuidaalrededor de un núcleo polinomial ó función de peso. En esta investigación
se usaron los siguientes núcleos polinomiales, obtenidos de la literatura en computación
gráfica [8], [20]:
W1(r, h) =
{
315
64πh9
(h2 − r2)3 si r < h
0 en caso contrario
(4.16)
27
W2(r, h) =
{
e−(
r
0,4h
)2 si r
h
≤ 1
0 si r
h
> 1
(4.17)
W3(r, h) =

2
3
− 4( r
h
)2 + 4( r
h
)3 si r
h
≤ 1
2
4
3
− 4( r
h
) + 4( r
h
)2 − 4
3
( r
h
)3 si 1
2
< r
h
≤ 1
0 si r
h
> 1
(4.18)
W4(r, h) =
{
1− 6( r
h
)2 + 8( r
h
)3 − 3( r
h
)4 si r
h
≤ 1
0 si r
h
> 1
(4.19)
W5(r, h) =

2
3h
(
1− 3
2
( r
h
)2 + 3
4
( r
h
)3
)
si r
h
≤ 1
2
3h
(
1
4
(2− r
h
)3
)
si 1 < r
h
≤ 2
0 si r
h
> 2
(4.20)
de estas ecuaciones se define Wij por:
Wij = W (‖ ~Xij‖2, hi), (4.21)
siendo el operador ‖ ‖2 la norma eucĺıdea. Las gráficas de estos núcleos se pueden observar
en las figuras 4.2.b, 4.3 y 4.4.
El ı́ndice j representa la part́ıcula vecina de la part́ıcula i. Las funciones asociadas a estos
núcleos relacionan una part́ıcula con sus vecinos, donde r es la distancia entre la part́ıcula i
y j, y h es el radio de soporte de estos núcleos.
La esfera centrada en la coordenada de material de la part́ıcula i tiene como radio de soporte
hi = 3r̄i, donde r̄i es la distancia promedio entre la part́ıcula i y sus vecinos. La densidad de
una part́ıcula i viene dada por ρi =
∑
j
mjWij y el volumen de una part́ıcula i se obtiene de
la forma siguiente:
voli =
mi
ρi
(4.22)
4.2.3. Cálculo de las Fuerzas
Las fuerzas internas están formadas por fuerzas elásticas, fuerzas de conservación del volumen
y fuerzas de Lennard-Jones. Las primeras son descritas de la siguiente manera:
28
Fei = (2voliJiσi)A
−1
i
∑
j
(
~XijWij
)
(4.23)
estas fuerzas conservan momento angular y lineal. Las fuerzas de conservación del volumen
se definen como:
Fvi = volikv(|Ji| − 1)ϕiA−1i
∑
j
(
~XijWij
)
(4.24)
donde kv es una constante de conservación del volumen que penaliza la desviación del
Jacobiano y
ϕi =
 (Jvi × Jwi)T(Jwi × Jui)T
(Jui × Jvi)T
 (4.25)
Las fuerzas de Lennard-Jones (Fli) [13] son fuerzas de atracción y repulsión entre part́ıculas.
Para cada part́ıcula i esta fuerza es representada por:
Fli =
∑
j
L(α, dij)~pij (4.26)
donde L(α, dij) = 4β
(
12α
12
d13ij
− 6 α6
d7ij
)
, β es el mı́nimo del potencial de enerǵıa de Lennard-
Jones, dij es la distancia entre la part́ıcula i y j (dij = ‖ ~Xi− ~Xj‖2), ~pij = 1‖ ~Xi− ~Xj‖2 (
~Xi− ~Xj)
y α es el diámetro de la esfera que representa la part́ıcula i. En la figura 4.2.a se puede ver
la gráfica de la función L(α, dij). En la figura 5.5 se puede observar que hi define el radio de
la esfera de la vecindad de la part́ıcula i y α va asociado sólo a la part́ıcula i.
Por lo tanto, las fuerzas internas de una part́ıcula i son representadas por la siguiente
expresión:
Finti = (Fei + Fvi)A
−1
i
∑
j
(
~XijWij
)
+ Fli (4.27)
En este estudio de modelos deformables, a las fuerzas internas se le pueden agregar fuerzas
externas tales como: gravedad, presión atmosférica entre otras. La suma de las fuerzas
internas y externas es la fuerza total en el modelo.
29
4.3. Dinámica
La dinámica del modelo de part́ıculas se basa en las leyes de Newton. Para la integración de
las fuerzas totales fue necesario utilizar el esquema expĺıcito de Euler. De esta manera, de la
ecuación
f t+4ti = mia
t+4t
i (4.28)
donde f t+4ti es la fuerza total, a
t+4t
i es la aceleración y mi la masa de la part́ıcula i, se
deducen tanto la velocidad
vt+4ti = v
t
i + a
t+4t
i 4t (4.29)
como la posición de una part́ıcula i
xt+4ti = x
t
i + v
t+4t
i 4t. (4.30)
En este trabajo también se utilizó un método de integración impĺıcita similar al usado en
el trabajo de Müller [20]. En este esquema, las nuevas posiciones y velocidades de todas las
part́ıculas son calculadas de la siguiente manera:
X t+4t = X t + V t+4t4t (4.31)
MV t+4t = MV t +4tF (X t+4t) (4.32)
donde los vectores X y V contienen todas las posiciones y velocidades de las part́ıculas en el
sistema. La matriz M es diagonal y contiene todas las masas. Sustituyendo la ecuación 4.31
en la ecuación 4.32 se obtiene el siguiente sistema:
MV t+4t = MV t +4tF (X t +4tV t+4t) (4.33)
de donde se deduce:
≈MV t +4tF (X t) +4t2K |Xt .V 4t+t
30
siendo F las fuerzas de todos los puntos y K |Xt= ∇XF (X t) el Jacobiano de F . K |Xt es
una matriz 3n× 3n donde n es la cantidad total de part́ıculas. Reordenando los términos se
obtiene:
(M −4t2K |Xt)V t+4t = MV t +4tF (X t) (4.34)
Resolviendo el sistema de la ecuación anterior, se pueden obtener las nuevas velocidades.
La submatriz Kkl, k, l ∈ [1..n] tiene la forma
Kkl = ∇ulfk = (
d
dul
fk,
d
dvl
fk,
d
dwl
fk) (4.35)
donde Kkl es una matriz 3× 3, siendo fk = Fek + Fvk, Fek son las fuerzas elásticas y Fvk
son las fuerzas de conservación del volumen de la part́ıcula k. La derivada de las fuerzas
elásticas para la componente ul :
d
dul
Fek = −2V olk
 dTl0
0
σk + JkC(JukdTl + dlJuTk )
 · dk (4.36)
Para la componente vl la expresión es:
d
dvl
Fek = −2V olk
 0dTl
0
σk + JkC(JvkdTl + dlJvTk )
 · dk (4.37)
Para la componente wl la expresión es:
d
dwl
Fek = −2V olk
 00
dTl
σk + JkC(JwkdTl + dlJwTk )
 · dk (4.38)
donde
31
C =
E
1− v2
 1 v 0v 1 0
0 0 1−v
2
 (4.39)
siendo E el módulo de Young y v la tasa de Poisson. Para las fuerzas de conservación del
volumen la derivada con respecto a ul es:
d
dul
Fek = −kvV olk
∣∣∣∣∣∣
 dTlJvTk
JwTk
∣∣∣∣∣∣ ·
 (Jvk × Jwk)T(Jwk × Juk)T
(Juk × Jvk)T
+ (|Jk| − 1)
 0(Jwk × dl)T
(dl × Jvk)T
 · dk
(4.40)
Para la componente vl la expresión es:
d
dvl
Fek = −kvV olk
∣∣∣∣∣∣
 JuTkdTl
JwTk
∣∣∣∣∣∣ ·
 (Jvk × Jwk)T(Jwk × Juk)T
(Juk × Jvk)T
+ (|Jk| − 1)
 (dl × Jwk)T0
(Juk × dl)T
 · dk
(4.41)
Para la componente wl la expresión es:
d
dwl
Fek = −kvV olk
∣∣∣∣∣∣
 JuTkJvTk
dTl
∣∣∣∣∣∣ ·
 (Jvk × Jwk)T(Jwk × Juk)T
(Juk × Jvk)T
+ (|Jk| − 1)
 (Jvk × dl)T(dl × Juk)T
0
 · dk
(4.42)
donde el operador | | representa el determinante de una matriz cuadrada.
Una vez calculadas las nuevas velocidades, usando la ecuación 4.31 se pueden obtener las
nuevas posiciones de todas las part́ıculas.
4.3.1. Dinámica de la Superficie del objeto 3D
El esquema utilizado en este trabajo para representar la superficie del objeto en deformación
se basa en la investigación realizada por Pauly [23]. Una vez calculadas las nuevas posiciones,
32
velocidades y desplazamientos de las part́ıculas internas (part́ıculas del modelo) se procede
a determinar el desplazamiento de los puntos de las superficie. El desplazamiento de una
part́ıcula de la superficie con respecto a los puntos internos vecinos se obtiene de la siguiente
manera:
usurf =
1∑
ω(ri)
∑
i
ω(ri)(~Ui +∇~UTi (xsurf − ~Xi)) (4.43)
donde xsurf es la posición de la part́ıcula de la superficie, ω(ri) = e
−r2i
h2 es una función de
peso gaussiano, ri = ||xsurf − ~Xi||2 y h es la distancia del radio de la esfera que contiene a
la part́ıcula de la superficie y sus part́ıculas internas vecinas.
Una vez calculado el nuevo desplazamiento de la part́ıcula de la superficie se procede a
obtener su nueva posición, la cual tiene la siguiente expresión:
xt+4tsurf = x
t
surf + usurf (4.44)
4.4. Resultados
Se crearon varios ejemplos sobre los cuales utilizamos esta metodoloǵıa para realizar
simulación dinámica de objetos deformables. Todas las animaciones fueron realizadas dentro
de cajas y cilindros. Por esta razón fue necesario aplicar técnicas de colisión de las part́ıculas
contra las cajas y cilindros que sirvieron de recipientes. La estrategia usada para la
iluminación de la superficie de estos objetos consistió en calcular las normales promediadas
por los triángulos vecinos a cada punto de la superficie. Otro aspecto visto en las simulaciones
fue que las fuerzas de Lennard-Jones permitieron que el manejo de colisiones entre las
part́ıculas fuese el adecuado.Para lograr estas animaciones fue necesario la calibración de
los parámetros f́ısicos para estos objetos.
En la figura (4.2) se muestra el potencial de Lennard-Jones junto
Los objetos de prueba 3D que utilizamos son los siguientes. Objeto 1: es un objeto 3D en
forma de un octaedro. Está formado por 299 part́ıculas internas, no tiene asignada superficie
que cubra el objeto y por lo tanto no existe ninguna triangulación de malla de la superficie.
Objeto 2: Es un objeto en forma de elipsoide formado por 568 part́ıculas y tampoco tiene
superficie que lo cubra. Objeto 5: es un objeto en forma de elipsoide formado por 715
part́ıculas, tiene superficie asignada formada por 990 nodos y una triangulación de 1976
33
triángulos.
En la tabla (4.1) se muestran las caracteŕısticas de los equipos utilizados (hardware) y los
sistemas operativos que éstos contienen. En la tabla (4.3) se describen las caracteŕısticas de
los objetos tridimensionales. En la figura 4.6 se muestra la deformación del objeto 1 usando
los equipos 1 y 2. Se empleó una integración del tiempo expĺıcita con un diferencial de tiempo
∆t = 0.005. El radio de soporte hi para definir la vecindad fue igual para todas las part́ıculas
con el valor de 0.8. Otros atributos como la densidad de material, la constante de volumen
kv, el módulo de Young y la tasa de Poisson tienen los siguientes valores 1, 10, 10 y 0.1
respectivamente.
En la figura 4.7, se aprecia la deformación del objeto 2, ejecutado en los equipos 1 y 2. En
este modelo se utilizó una tasa de Poisson 0.1, un módulo de Young de 10, kv = 10, una
presión atmosférica de 7. Se realizó con una integración expĺıcita con ∆t = 0.009.
En la tabla (4.2) se puede apreciar tiempos en FPS (cuadros por segundo) y por iteración
o ciclo de animación (TIS - Tiempo por iteración en segundos). Estos TIS solo calculan el
tiempo de ejecución del método sin realizar la visualización, el cual implica no solo pintar el
modelo sino también la iluminación. En esta tabla se puede apreciar que los TIS del objeto 1
son mayores que los del objeto 2. Esto ocurre porque a medida que aumentan la cantidad de
part́ıculas los cálculos son mayores produciendo que los FPS disminuyan y los TIS aumenten.
En la figura 4.8 se puede apreciar la simulación del objeto 5 realizada en los equipos 1 y 4. Se
utilizó una integración del tiempo expĺıcita. La función de peso usada para esta simulación
es mostrada en la ecuación 4.16.
Figura 4.1: Objeto en estado de reposo y en estado de deformación.
La figura del lado izquierdo muestra el objeto sin deformaciones. La figura del lado derecho
representa un objeto con su deformación (X + U). Donde X es el conjunto de part́ıculas que
representan las coordenadas del material y X + U son las coordenadas de las part́ıculas del
objeto deformado.
34
Equipo Procesador Velocidad Memoria Arquitectura Aceleradora Sistema
(GHZ) RAM (bits) Gráfica Operativo
1 AMD 2.20 512 MB 64 · · · · · · · · · Windows
Athlon XP
2 Intel 2.80 752 MB 32 Intel 82852 Windows
Pentium 4 82855 GM/GME XP
Graphics Controller
3 UltraSPARC 1.2 1 GB 64 · · · · · · · · · Solaris
III Cu 10
4 2 Dual Core 2.4 7.79 GB 64 NVIDIA Quadro Suse 10.2
AMD Opteron FX3450 (Linux)
/4000 SID
Cuadro 4.1: Equipos utilizados.
Objeto FPS TIS Equipo
1 7.2944 0.09926 1
1 5.9943 0.1320 2
2 1.2751 0.5123 1
2 1.029 0.6683 2
Cuadro 4.2: Tiempos de las animaciones basadas en el modelo de part́ıculas internas usando
la estrategia de búsqueda de vecinos de todos contra todos en los distintos equipos.
Parámetros Objeto 1 Objeto 2 Objeto 5
∆t 0.001 0.001 0.001
Integración del tiempo expĺıcita expĺıcita expĺıcita
Max. de vecinos para hi 10 10 10
Módulo de Young 10000 100 100
Tasa de Poisson 0.1 0.2 0.2
Constante de Volumen kv 100 100 100
Densidad de Material 10 5 5
Presión Atmosférica 10 10 100
Radio f́ısico de las part́ıculas 0.4 0.2 0.2
Radio de vecindad de puntos de la superficie 2 1.5 2.5
Función de Peso W5 W1 W1
Cuadro 4.3: Parámetros usados en las animaciones basadas en el modelo de los puntos
internos y de superficie.
35
Figura 4.2: Función de Lennard-Jones L(α, dij). Función W1(r, h).
La imagen (a) muestra la función L(α, dij). La imagen (b) muestra la función W1.
Figura 4.3: Funciones de peso W2(r, h) y W3(r, h).
La imagen (a) muestra la función W2. La imagen (b) muestra la función W3.
Figura 4.4: Funciones de peso W4(r, h) y W5(r, h).
La imagen (a) muestra la función W4. La imagen (b) muestra la función W5.
Figura 4.5: Diferencia entre hi y α.
36
Figura 4.6: Animación de un octaedro con 299 part́ıculas internas.
Figura 4.7: Animación de un elipsoide con 568 part́ıculas internas.
Las secuencias de las dos animaciones son mostradas de izquierda a derecha.
37
Figura 4.8: Animación de un elipsoide formado por 715 part́ıculas internas y 990 de la
superficie y 1976 triángulos.
En la parte superior se muestran las part́ıculas internas y externas del objeto. La animación
de la deformación es mostrada en la parte inferior.
38
Caṕıtulo 5
Simulación de Fluidos Basados en
Part́ıculas usando SPH
En computación gráfica el flujo de fluidos es importante a la hora de simular fenómenos
presentes en el d́ıa a d́ıa, por ejemplo, agua, lluvia, olas, lodo. En este caṕıtulo se presentan
las ecuaciones de Navier-Stokes para el flujo de fluidos que servirán de basamento para el
presente trabajo.
Las ecuaciones de Navier-Stokes se pueden resolver usando métodos con mallas también
llamados métodos eulerianos, o bien usando métodos sin mallas, basados en part́ıculas,
también llamados métodos lagrangianos.
Los métodos eulerianos basados en mallas son los más usados en computación grááfica y las
ecuaciones de Navier-Stokes han sido extendidas para simular muchos efectos interesantes
como tensión superficial, visco-elasticidad y elasto-plasticidad. La mayoŕıa de las soluciones
eulerianas y sus extensiones son poco útiles a la hora de realizar simulaciones interactivas en
tiempo real.
5.1. Las Ecuaciones de Navier-Stokes
Las propiedades t́ıpicas de un fluido isotérmico viscoso son: velocidad u, masa-densidad ρ y
presión p. Las propiedades f́ısicas son consideradas como campos continuos en el fluido. La
formulación clásica de la dinámica para el flujo de un fluido incompresible en el tiempo t
viene dada por las ecuaciones de Navier-Stokes.
39
ρ(
∂
∂t
+ u · ∇)u = −∇p+ µ∇ · (∇u) + f , (5.1)
∇ · u = 0, (5.2)
donde µ es la viscosidad y f es la suma de las densidades-fuerzas externas actuando sobre
el fluido, como por ejemplo, la gravedad.
La formulación de Navier-Stokes para el flujo de fluidos está basado en una estructura de
mallas, lo que significa que no solo depende del tiempo t, sino también de la posición en
la malla r(t), que también dependen del tiempo. En tres dimensiones el vector de posición
viene dado por:
r(t) = (x(t), y(t), z(t))T . (5.3)
La ecuación (5.1) describe la conservación de momento y es básicamente la segunda ley
de Newton para un fluido. A la derecha de la ecuación están representados el total de las
densidades-fuerzas mientras que del lado izquierdo está el producto de las masas-densidades
por la densidad-aceleración. Para un fluido la densidad-aceleración es la derivada del campo
de velocidad. Usando la regla de la cadena tenemos:
d
dt
u(t, r(t)) =
∂u
∂t
+
∂u
∂x
∂x
∂t
+
∂u
∂y
∂y
∂t
+
∂u
∂z
∂z
∂t
=
∂u
∂t
+
[
∂u
∂x
· ∂u
∂y
· ∂u
∂z
]T
·
[
∂x
∂t
,
∂y
∂t
,
∂z
∂t
]T
=
∂u
∂t
+ (u ·
[
∂
∂x
· ∂
∂y
· ∂
∂z
]
)u
= (
∂
∂t
+ u · ∇)u. (5.4)
La velocidad de cambio de un campo continuo adimensional (5.4), en el fluido, es llamada
la derivada sustancial de ψ y esta definida como:
dψ
dt
=
∂ψ
∂t
+ u · ∇ψ, (5.5)
Donde u · ∇ψ es el término de advección y (5.4) justifica la aparición de éste. La derivada
sustancial es la derivada de tiempo completo del campo ψ en el enfoque euleriano.
40
La ecuación (5.2) describe la conservación de masa para un fluido incompresibley fue
formulada originalmente en las ecuaciones de Euler. En general, la ecuación de conservación
de masa para fluidos compresibles es conocida como la ecuación de continuidad. Es la
derivada sustancial del campo de densidad, lo que nos lleva a:
∂ρ
∂t
+∇ · (ρu) = 0. (5.6)
Cuando la masa-densidad es constante en el fluido ∂p
∂t
= 0 y por lo tanto ρ · ∇u = ∇ · u,
entonces el fluido es incompresible.
5.1.1. Enfoque Euleriano
Las ecuaciones de Navier-Stokes (5.1) y (5.2) formulan la dinámica de un fluido. El método
euleriano considera un fluido como un conjunto de celdas alineadas en una malla regular,
donde cada celda contiene una cantidad de moléculas de dicho fluido. El reto de este método
para la simulación consiste en superar las dificultades en la visualización debido al grosor de
la malla.
La solución anaĺıtica de las ecuaciones de Navier-Stokes es todo un desaf́ıo, pero se pueden
encontrar soluciones numéricas usando varios pasos para cada componente de las ecuaciones.
Uno de los mas importantes es la descomposición de Helmholtz-Hodge [6], [10], que propone
una técnica matemática que permite representar el gradiente de presión (−∇p) y la ecuación
(5.2) en una sola ecuación. Los diferentes pasos se resuelven usando integradores expĺıcitos,
impĺıcitos y semi-impĺıcitos. La malla provee de una solución a las derivadas usando el método
de diferencias finitas.
5.2. Smoothed Particle Hydrodynamics (SPH)
La formulación de SPH surge del campo de la astrof́ısica computacional y está diseñada para
representar problemas de flujo de fluidos compresibles. El método es ampliamente utilizado
en fenómenos del campo de la astrof́ısica, por su facilidad de representar problemas complejos
de una manera sencilla. SPH es un método de interpolación para aproximar los valores y
derivadas de campos continuos usando puntos discretos de muestreo.
Los puntos de muestreo son identificados como part́ıculas y contienen valores concretos, por
41
ejemplo, masa, velocidad, posición, etc., pero estas part́ıculas además poseen cantidades
estimadas de propiedades f́ısicas dependientes del problema como, por ejemplo, masa-
densidad, temperatura, presión, etc. Las propiedades en el método SPH son macroscópicas
y son obtenidas como los promedios pesados entre las part́ıculas vecinas.
Comparado con otros métodos conocidos para la aproximación numérica de las derivadas,
entre ellos, el método de diferencias finitas, que requiere que las part́ıculas se encuentren
alineadas en una malla regular, SPH puede aproximar derivadas de campos continuos usando
diferenciación anaĺıtica en part́ıculas con posiciones arbitrarias. Cada part́ıcula está pensada
como que ocupa una fracción del espacio del problema, y para obtener promedios de
cantidades mas precisos, las part́ıculas muestreadas deben ser numerosas.
El método SPH es básicamente un método de interpolación, esta interpolación está basada
en la teoŕıa de integrales interpolantes usando funciones de peso para aproximar una función
delta. La integral interpolante de cualquier función, A(r), está definida sobre todo el espacio,
Ω, porción
AI(r) =
∫
Ω
A(r′)W (r− r′, h)dr′, (5.7)
donde r es un punto en Ω, y W es una función de peso con h como su radio. El radio, es un
factor de escala que controla la suavidad de la función de peso.
El equivalente numérico de 5.7 es obtenido mediante la aproximación de la integral
interpolante con una sumatoria interpolante,
AI(r) =
∑
Ω
AjVjW (r− rj, h), (5.8)
donde j es iterada por todas las part́ıculas, Vj es el volumen atribuido impĺıcitamente a
la part́ıcula j, rj es la posición y Aj es la cantidad A en rj. En comparación con 5.7 el
volumen de cada part́ıcula es acumulado en 5.8 lo que puede ser justificado cuando Aj se
asume suficientemente constante en Vj. La siguiente relación entre la masa, el volumen y la
masa-densidad aplica
V =
m
ρ
, (5.9)
donde m es la masa y ρ es la masa-densidad. Sustituyendo el volumen para la part́ıcula j
42
en 5.8 con 5.9 tenemos la formulación básica del método SPH, que puede ser utilizada para
aproximar cualquier cantidad en un campo continuo, y ser evaluada en cualquier lugar del
espacio contenido,
As(r) =
∑
j
Aj
mj
ρj
W (r− rj, h). (5.10)
En SPH un interpolante diferenciable de una función puede ser construido a partir de sus
valores en las part́ıculas usando un núcleo diferenciable, y por lo tanto las derivadas de una
interpolante pueden ser obtenidas mediante la diferenciación anaĺıtica del núcleo. Por lo tanto
es un método directo que permite definir el gradiente y laplaciano de un campo, mediante la
asumpción de que el núcleo es diferenciable con derivadas continuas hasta el segundo orden.
Considerando un sistema de dos part́ıculas i y j, la derivada parcial en (5.10) para x resulta
en
∂
∂x
AS(ri) =
∂
∂x
(
Aj
mj
ρj
W (ri − rj, h)
)
(5.11)
usando la regla del producto veremos
∂
∂x
(
Aj
mj
ρj
W (ri − rj, h)
)
=
∂
∂x
(
Aj
mj
ρj
)
W (ri − rj, h) + Aj
mj
ρj
∂
∂x
W (ri − rj, h)
= 0 ·W (ri − rj, h) + Aj
mj
ρj
∂
∂x
W (ri − rj, h)
= Aj
mj
ρj
∂
∂x
W (ri − rj, h), (5.12)
donde se utiliza el hecho de que Aj
mj
ρj
en este punto no depende directamente de x ni de
ningún otro eje del espacio, por lo tanto el producto se considera constante.
Sabiendo que las derivadas de una sumatoria interpolante solamente afectan al núcleo, el
gradiente de un campo resulta en
∇AS(r) =
∑
j
Aj
mj
ρj
∇W (r− rj, h). (5.13)
43
Para obtener mayor precisión en el gradiente de un campo, la interpolante puede ser obtenida
usando
∇A = 1
ρ
(∇(ρA)− A∇ρ), (5.14)
La expresión en SPH para 5.14, usando 5.13 en el termino del gradiente resulta en
∇AS(r) =
1
ρ
(∑
j
pjAj
mj
ρj
∇W (r− rj, h)− A
∑
j
ρj
mj
ρj
∇W (r− rj, h)
)
(5.15)
=
1
ρ
(∑
j
Ajmj∇W (r− rj, h)− A
∑
j
mj∇W (r− rj, h)
)
(5.16)
Una forma particular de 5.14 puede ser obtenida reescribiendo ∇A
ρ
, de acuerdo con
∇A = ρ
(
∇
(
A
ρ
)
+
A
ρ2
∇ρ
)
. (5.17)
donde el termino del SPH se convierte en
∇AS(r) = ρ
(∑
j
Aj
ρj
mj
ρj
∇W (r− rj, h) +
A
ρ2
∑
j
ρj
mj
ρj
∇W (r− rj, h)
)
= ρ
(∑
j
Aj
ρ2j
mj∇W (r− rj, h) +
∑
j
A
ρ2
mj∇W (r− rj, h)
)
= ρ
∑
j
(
Aj
ρ2j
+
A
ρ2
)
mj∇W (r− rj, h). (5.18)
El laplaciano del campo 5.10 puede ser obtenido derivando el gradiente, de donde resulta
∇2AS(r) =
∑
j
Aj
mj
ρj
∇2W (r− rj, h). (5.19)
44
Núcleos del SPH
El uso de diferentes núcleos en SPH es equivalente al uso de diferentes esquemas de
diferenciación en diferencias finitas, por lo tanto la elección de un buen núcleo para un
problema espećıfico es sumamente importante. Los núcleos en SPH deben cumplir con dos
propiedades,
∫
Ω
W (r, h)dr = 1 (5.20)
y
ĺım
h→0
W (r, h) = δ(r). (5.21)
donde δ es la función del delta de dirac
δ(r) =
{
∞ ‖r‖ = 0
0 en caso contrario
(5.22)
En la ecuación (5.20) observamos que los núcleos deben estar normalizados. Los núcleos
además deben ser positivos
W (r, h) ≥ 0, (5.23)
para asegurar que son funciones de peso. Si el núcleo es par,
W (r, h) = W (−r, h), (5.24)
se confirma la simetŕıa rotacional, lo que es importante para evitar que vaŕıe bajo rotaciones
del sistema de coordenadas. Śı (5.24) y (5.20) se cumplen, el núcleo es par y esta normalizado,
entonces el error es O(h2). También se sugiere que los núcleos estén definidos sobre soportes
compactos del radio. Se utiliza h como el radio de soporte de los núcleos, por lo que se
cumplirá W (r, h) = 0, ‖r‖ > h.
45
5.3. Enfoque Lagrangiano
Las ecuaciones de Navier-Stokes representan la formulación euleriana básica para un fluido
isotérmico incompresible. Usando part́ıculas en vez de una malla se puede simplificar la
ecuación enormemente. Se asume que la cantidad de part́ıculas es constante durante la
simulación, y manteniendo la masa constante en las part́ıculas la conservación de masa esta
garantizada, por lo tanto la ecuación ∇ · u = 0 puede ser omitida.
Figura 5.1: Configuración Básica de SPH para las ecuacionesde Navier-Stokes
La figura muestra la configuración básica de un fluido basado en part́ıculas. En la formulación
lagrangiana de un fluido, son las part́ıculas quienes lo definen completamente, por lo tanto
las part́ıculas se mueven con el fluido.
Comparativamente el enfoque euleriano mantiene su malla estática durante la simulación y la
simulación depende de la posición de un nodo en la malla además del tiempo t, con el enfoque
lagrangiano cualquier propiedad del campo que representa al fluido depende solamente
del tiempo t. Las part́ıculas contienen masa, posición y velocidad, además de cantidades,
asociadas a propiedades de material, aproximadas mediante SPH. La aceleración para una
part́ıcula de fluido lagrangiano proviene de la derivada ordinaria d
dt
de su velocidad u(t). Esto
explica el por qué el término de advección no está presente en el enfoque lagrangiano. La
formulación lagrangiana básica de Navier-Stokes para un fluido incompresible e isotérmico
viene dada por
46
ρ
du
dt
= −∇p+ µ∇2u + f . (5.25)
El lado derecho de la ecuación consiste en las fuerzas internas y externas. Los campos de
fuerzas pueden ser combinados en una suma de fuerzas, F = f internal + f external. Para la
part́ıcula i, la aceleración se define como
ai =
dui
dt
=
Fi
ρi
. (5.26)
donde ai y ui son la aceleración y la velocidad para la part́ıcula i, respectivamente, Fi es el
total de fuerza aplicada sobre la part́ıcula y ρi es la masa-densidad evaluada en la posición
de la part́ıcula i.
Es necesario seleccionar un núcleo que nos ayude a resolver las ecuaciones de SPH, para
estas razones se ha elegido el núcleo polinomial de sexto grado, que viene dado por
Wdefault(r, h) =
315
64πh9
{
(h2 − r2)3 si r < h
0 en caso contrario
(5.27)
con gradiente
∇Wdefault(r, h) =
945
32πh9
r(h2 − ‖r‖2)2, (5.28)
y laplaciano
∇2Wdefault(r, h) =
945
32πh9
(h2 − ‖r‖2)(3h2 − 7‖r‖2). (5.29)
Existen algunas caracteŕısticas interesantes de este núcleo, tales como que preserva la
curvatura de campana de Gauss y la norma de r puede ser omitida completamente, haciéndolo
computacionalmente ligero. El núcleo por defecto Wdefault(r, h) y sus derivadas son utilizadas
para el cálculo de todas las aproximaciones de las propiedades del campo, excepto por las
fuerzas internas del campo.
47
5.3.1. Masa-Densidad
Las masas de las part́ıculas que son requeridas por las ecuaciones de SPH, son constantes
definidas por el usuario, pero la masa-densidad es un campo continuo del fluido, que debe
ser calculado. Para evitar cualquier confusión, la masa-densidad es el campo que únicamente
depende de la masa de las part́ıculas y para este caso especifico se usa la ecuación (5.10)
para calcularlo. En la part́ıcula i la masa-densidad se expresa aśı
ρi = ρ(ri) (5.30)
=
∑
j
ρj
mj
ρj
Wdefault(r, h) (5.31)
=
∑
j
mjWdefault(r, h). (5.32)
5.3.2. Fuerzas Internas
Las fuerzas internas son las fuerzas definidas por el fluido mismo, éstas se representan como
presión y viscosidad, que son el primero y segundo termino de la ecuación 5.25 en su lado
derecho. Las fuerzas internas del fluido son aproximadas usando la formulación del enfoque
SPH.
5.3.2.1. Presión
La presión p en una part́ıcula puede ser determinada usando la ley del gas ideal
pV = nRT, (5.33)
donde V = 1
ρ
es el volumen por unidad de masa, n es el numero de part́ıculas de gas en
mol, R es la constante de gas universal y T es la temperatura. Para un fluido isotérmico, con
masa constante el lado derecho de la ecuación se puede mantener como constante, y por lo
tanto lo reemplazamos con una constante de rigidez k, que en teoŕıa solamente depende de la
cantidad de part́ıculas en el fluido. El término de presión puede ser reescrito de la siguiente
manera
48
pV = k
p
1
ρ
= k
p = kρ. (5.34)
Si se conoce la presión en cada part́ıcula, la fuerza local de presión sobre la part́ıcula i en
notación SPH es
fpresioni = −∇p(ri) (5.35)
= −
∑
j 6=i
pj
mj
ρj
∇Wdefault(ri − rj, h) (5.36)
Esta fuerza no es simétrica lo que llevaŕıa a la no conservación de la ley de acción-reacción,
por lo tanto reescribiremos el término −∇p(ri) usando 5.18 para obtener una fuerza simétrica
que si conserve esta ley, quedando la ecuación anterior de la siguiente manera
fpresioni = −pi
∑
j 6=i
(
pi
ρ2i
+
pj
ρ2j
)
mj∇Wdefault(ri − rj, h). (5.37)
Usar la ecuación 5.34 para el calculo de la presión resulta en un comportamiento de gas
ideal que se expande en el espacio, pero como el objetivo es representar ĺıquidos, se debe
reformular la ecuación para que el fluido posea cohesión interna y tenga una masa-densidad
constante en reposo. Para ello se puede utilizar una versión modificada de la ecuación de gas
ideal con una presión de reposo p0, que es
(p+ p0)V = k (5.38)
p+ kρ0 = kρ (5.39)
p = k(ρ− ρ0), (5.40)
donde ρ0 es la densidad en reposo del fluido. Usando esta ecuación en (5.37) se obtendrá un
comportamiento en el cual, cuando las part́ıculas se acerquen a la densidad de reposo las
fuerzas de atracción repulsión se balancearan.
49
El gradiente del campo de presión se utiliza para calcular la fuerza de presión. Para esto
es necesario un núcleo especial, debido a que con el núcleo por defecto se obtiene una
acumulación poco realista de part́ıculas en las zonas de alta presión, esta misma situación se
le presento a Desbrum et al. Para lo que propuso otro núcleo normalizado, llamado núcleo
“spiky”, que sera utilizado como núcleo para la presión
Wpresion(r, h) =
15
πh6
{
(h− r)3 si 0 ≤ r ≤ h
0 si r > h
(5.41)
con gradiente
∇Wpresion(r, h) =
15
πh6
r
‖r‖
(h− ‖r‖)2 (5.42)
ĺım
r→0−
∇Wpresion(r, h) =
45
πh6
, ĺım
r→0+
∇Wpresion(r, h) = −
45
πh6
y laplaciano
∇2Wpresion(r, h) =
90
πh6
1
‖r‖
(h− ‖r‖)(h− 2‖r‖) (5.43)
ĺım
r→0
∇2Wpresion(r, h) = −∞,
5.3.3. Viscosidad
Un fluido es una sustancia que no posee resistencia para fuerzas de quiebre y por lo tanto fluye
a través de las deformaciones. Al mismo tiempo cuando el fluido fluye, las moléculas sufren
una fricción interna que reduce su enerǵıa cinética transformándola en calor. La resistencia
al flujo es llamada viscosidad, que está controlada por el coeficiente de viscosidad µ, que
define la fuerza de viscosidad sobre fluido. La variante de SPH para el término de la fuerza
de viscosidad se define
fviscosidadi = µ∇2u(ri) (5.44)
= µ
∑
j 6=i
uj
mj
ρj
∇2W (ri − rj, h). (5.45)
50
Similar a la fuerza de presión, la fuerza de viscosidad es también asimétrica, debido a que la
velocidad varia de part́ıcula a part́ıcula. Para contrarrestar este efecto, Mullr et al. simetriza
el campo de velocidad usando
fviscosidadi = µ
∑
j 6=i
(uj − ui)
mj
ρj
∇2W (ri − rj, h). (5.46)
Esto es posible debido al hecho de que la fuerza de viscosidad solamente depende de las
diferencias de velocidades y no de las velocidades absolutas.
El laplaciano del núcleo en la ecuación anterior debe ser siempre positivo. Esto es necesario
debido a que la fuerza de viscosidad no debe incrementar las velocidades relativas, pues
esto, introduciŕıa enerǵıa al sistema, causando inestabilidad. Solo si el laplaciano es positivo
definido entonces la fuerza de viscosidad actuará como amortiguamiento para las velocidades
relativas.
El núcleo por defecto no tiene esta propiedad ni tampoco el núcleo de presión, por lo tanto
se debe emplear un núcleo que cumpla ∇2W (r, h) ≥ 0 para ‖r‖ ≤ 0, el núcleo que cumple
estas propiedades es
Wviscosidad(r, h) =
15
2πh3
{
−‖r‖
3
2h3
+ ‖r‖
2
h2
+ h
2‖r‖ − 1 0 < ‖r‖ ≤ h
0 ‖r‖ > h
(5.47)
ĺım
r→0
Wviscosidad(r, h) =∞,
con gradiente
∇Wviscosidad(r, h) =
15
2πh3
r
(
−3‖r‖
2h3
+
2
h2
− h
2‖r‖3
)
, (5.48)
ĺım
r→0−
∇Wviscosidad(r, h) =∞, ĺım
r→0+
∇Wviscosidad(r, h) = −∞,
y el laplaciano
∇2Wviscosidad(r, h) =
45
πh6
(h− ‖r‖), (5.49)
51
5.3.4. Fuerzas Externas
Las fuerzas externas son todas aquellas fuerzas que actúan sobre el modelo, estás se balancean
con las fuerzas internas del modelo y corresponden al último término