JULIO / DICIEMBRE
2019 -
VOLUMEN 29 (2)
10.15517/ri.v29i2.36140
Ingeniería 29 (2): 44-58, julio-diciembre, 2019. ISSN: 2215-2652. San José, Costa Rica.
Esta obra está bajo una Licencia de Creative Commons. Atribución - No Comercial - Compartir Igual
Método EBE de los elementos nitos para sistemas de gran tamaño
y su aplicación en la Física Digital de Rocas
EBE nite element methdology for large scale systems applied to
digital rock Physics
Francisco José Benavides Murillo
Universidad Nacional de Costa Rica, Costa Rica.
Correos: fbenavides@id.uff.br
Recibido: 2 de febrero 2019 Aceptado: 5 de julio 2019
_________________________________________________________
Resumen
El método de los elementos nitos (FEM) es una técnica numérica ampliamente utilizada en física e
ingeniería, que aproxima soluciones de ecuaciones diferenciales en derivadas parciales en dominios arbitrarios.
En general, el paso nal de esta técnica consiste en la resolución de un sistema lineal de la forma Ax=b, en
el que la matriz A es dispersa y cuyo ancho de banda (número de las no nulas en la matriz) depende del
soporte de las funciones de forma denidas en los elementos. Cuando este número de elementos es muy
grande, inclusive las representaciones en estructuras de datos ecientes para matrices dispersas son capaces
de consumir en su totalidad la memoria de un computador poderoso. Por tanto, en este artículo se describe
el método EBE (Element by Element), un truco computacional diseñado para implementar el método FEM
sin utilizar la representación explícita de esta matriz. Se describe también un contexto en el que es necesaria
esta implementación, en el área de la petrofísica digital de rocas, para estimar las propiedades elásticas de
una muestra de roca a partir de imágenes tomográcas.
Palabras clave:
Análisis Numérico, Elemento Finito, Elemento por Elemento, Petrofísica, Propiedades elásticas
Abstract
The nite element method (FEM) is a numerical technique that estimates solutions of partial differential
equations on arbitrary domains. It has been widely used to solve problems in Physics and Engineering. In
general, the nal step of this technique consists of a linear system of equations Ax=b in which the matrix
A is sparse and its bandwidth depends on the nite element shape functions support. When the number of
elements is large, even the efcient data structure sparse matrix representations can consume the entire
computer memory. In this article, we describe a technique to solve these large-scale problems without
explicitly representing this matrix. This computational trick is known as EBE (Element By Element). We
also describe an application in which such kind of implementation is necessary, in the eld of digital rock
physics, to estimate the elastic coefcient of rock samples using micro-tomographic images.
Keywords:
Numerical Analysis, Finite Element, Element By Element, Petrophysics, Elastic properties
Ingeniería 29 (2): 44-58, julio-diciembre, 2019. ISSN: 2215-2652. San José, Costa Rica. DOI 10.15517/ri.v29i2.36140
45
1. INTRODUCCIÓN
El método de los elementos nitos (FEM) es una técnica proyectiva, en la que una función
ϕ
∈
X de un espacio de Hilbert X, es aproximada de manera óptima en un sub-espacio U ⊂ X mediante
su única proyección
ϕ
U
(Kress, 1999). Por una parte, la función
ϕ
es la solución de una ecuación
diferencial en derivadas parciales sobre un dominio arbitrario, dominio que puede modelar un objeto
físico real. Por otra parte U es, en general, un sub-espacio de dimensión nita cuya base está formada
por funciones de soporte compacto llamadas funciones de forma.
Tales funciones suelen ser denidas con base en subdominios o elementos geométricos (triángulos
o hexaedros, por ejemplo). La aproximación minimiza la norma ||
ϕ
-
ϕ
U
|| asociada a la forma variacional
del problema, la cual se calcula en términos de productos internos en el espacio de funciones U
(Hughes, 2000). La forma variacional representa la energía asociada al problema físico modelado;
por esta razón, el FEM es ampliamente utilizado en aplicaciones de ingeniería, en áreas tales como
Mecánica de Materiales y Teoría de Elasticidad (Slaughter, 2002).
También, el FEM es utilizado en la simulación de procesos físicos en la Física de Rocas; ciencia
que relaciona observaciones geofísicas obtenidas remotamente con las propiedades físicas subyacentes
de las rocas, tales como composición, porosidad, permeabilidad, coecientes elásticos, etc (Mavko
et al., 2003).
Las propiedades mecánicas de las rocas son útiles para estimar sus procesos de deformación y
son aplicadas en la Ingeniería de Petróleo para optimizar los métodos de perforación y extracción
(Tiab and Donaldson, 2012). Estas propiedades pueden ser obtenidas a partir de ensayos de
laboratorio, pero recientemente esto se ha hecho a partir de imágenes tomográcas de muestras
de las rocas del reservatorio. Obtener tales imágenes para simular, en un computador, los procesos
físicos del espacio poroso y mineral de la estructura rocosa es el paradigma de la Física Digital
de las Rocas.
Cabe aclarar que este artículo se basa en el enfoque descrito en Garboczi y Day (1995) para
estimar el módulo de Young o módulo de Elasticidad de una muestra de roca, digitalizada mediante
una secuencia de imágenes tomográcas. Asimismo, utiliza la versión más simple de la técnica de
los elementos nitos, utilizando apenas interpolaciones lineales entre los nodos y métodos clásicos
de gradiente para resolver el sistema lineal asociado.
El interés principal de la investigación no está dirigido hacia la justicación del uso de los
elementos nitos para resolver el problema, ni en proponer algún método novedoso para el FEM. Se
considera que ese tipo de justicaciones ya han sido abordadas, por ejemplo en (Arns et al., 2002;
Garboczi and Day, 1995; Makarynska et al., 2008; y Roberts y Garboczi, 2000). En dichos artículos,
se utiliza FEM para estimar el módulo de Young, utilizando una cantidad muy grande de grados de
libertad (deniendo un elemento por voxel), con sistemas de ecuaciones lineales que pueden llegar
fácilmente a más de cien millones de incógnitas.
Asimismo, las investigaciones mencionadas no describen con detalle los aspectos computacionales
que hay que tener en cuenta a la hora de resolver sistemas de ese tamaño, lo que puede llevar a la
BENAVIDES MURILLO: Método EBE de los elementos nitos...
46
conclusión errónea de que son necesarios super computadores de gran tamaño para su tratamiento.
En realidad, sistemas de este orden pueden ser resueltos en Workstations comunes con la metodología
que se describe en este artículo, la cual, de hecho, ha sido utilizada con éxito en (Leiderman et al.,
2017). Cabe resaltar que el rango de aplicaciones se puede extender a otros tipos de problemas
mecánicos (Betancur et al., 2018).
El truco computacional descrito en este trabajo no es novedoso, pero representa una marcada
diferencia entre la utilización óptima de los recursos computacionales y la falsa necesidad de aumentar
la capacidad del hardware con el n de abordar este tipo de problemas.
Vale la pena remarcar que, debido a la existencia actual de una gran cantidad de implementaciones
computacionales del FEM, tanto comerciales como libres, su utilización para estimar el módulo de
Young a partir de imágenes parece una escogencia obvia, natural y preferible. Sin embargo, como se
verá a lo largo de este artículo, la mayor parte de estos programas hacen implementaciones genéricas
de los elementos nitos que no optimizan las características especiales del caso aquí presentado. De
hecho, su utilización para simular un ensayo de laboratorio no es posible para tamaños relativamente
grandes, que son normales en el contexto de la Física Digital de Rocas, y una implementación genérica
del FEM agota fácilmente los recursos de la máquina.
De este modo, la implementación descrita en el presente artículo surge a raíz de la necesidad de
simular numérica y computacionalmente un experimento de laboratorio (especícamente un ensayo
de compresión uniaxial, como se describe en las siguientes secciones) con base en una secuencia de
imágenes de una muestra de roca.
Tales imágenes representan, con precisión, la estructura del medio poroso, y permiten
realizar múltiples ensayos virtuales no destructivos en la muestra. La utilización de imágenes
para simular experimentos de laboratorio tiene algunas ventajas de rapidez, eciencia y costo,
pero también tiene el interés cientíco de conducir a una mejor comprensión del fenómeno
físico modelado (Andrä et al., 2013). Entre los muchos ejemplos que se pueden citar, se puede
considerar el trabajo de (Benavides et al., 2017), donde se propone una metodología computacional
para mejorar la estimación de un parámetro de la roca, utilizando experimentos virtuales en
imágenes tomográcas.
El resto del artículo se distribuye de la siguiente manera: en el marco teórico se hace una
descripción sucinta del ensayo de compresión uniaxial; luego, se describe el método FEM desde
el punto de vista de la implementación para este tipo de ensayo. Se aborda la resolución, por
medio de FEM, de la ecuación elastostática, describiendo los operadores y la forma de calcularlos.
Posteriormente, se describe el método EBE de manera detallada, tratando de definir, paso por
paso, un algoritmo implementable. Se tratan de llenar los vacíos de otras descripciones menos
detalladas, por ejemplo, (Leiderman et al., 2017; Ulrich et al., 1998. Luego se describen otros
detalles computacionales, como la forma de paralelizar y validar; finalmente, se presentan las
conclusiones.
Ingeniería 29 (2): 44-58, julio-diciembre, 2019. ISSN: 2215-2652. San José, Costa Rica. DOI 10.15517/ri.v29i2.36140
47
2. MARCO TEÓRICO
2.1 Ecuación elastostática y el ensayo de compresión uniaxial
La prueba de compresión uniaxial es un ensayo de laboratorio en el que se estima el coeciente
de Young de una muestra de roca (Hudson and Harrison, 2000). Este mide la compresibilidad
homogénea equivalente a una muestra heterogénea de roca, compuesta por diferentes tipos de
materiales con diferentes propiedades mecánicas y por un espacio de poros vacíos con geometría
arbitraria. La prueba consiste en someter la muestra a una deformación ja u_y a lo largo del eje
Y, manteniendo libertad de movimiento en los otros ejes (Figura 1, izquierda).
Para lograr tal deformación, se ha necesitado aplicar una presión en el área superior, dando
como resultado una fuerza total F. Tal fuerza se relaciona con la deformación mediante la
ecuación: , donde E es el coeciente de Young efectivo, h es la altura de la muestra y A
es el área de la parte superior.
Dicha ecuación dene la ley de Hooke en una dimensión. Para calcular esta fuerza se utiliza la
solución del campo de desplazamiento de la ecuación elastostática. Con este campo, es posible calcular
las presiones en la supercie superior de la muestra y, mediante integración numérica, calcular la
fuerza en toda la supercie. Las condiciones de contorno son descritas mediante el desplazamiento
u
y
y la jación de la supercie inferior de la muestra. No hay otras presiones externas. Un modelo
de muestra de roca (Figura 1, derecha) incluye gran cantidad de detalles del medio poroso, por lo
que la malla debe ser na. Habitualmente se utiliza un voxel por elemento, lo que produce sistemas
de gran tamaño. Esto se aborda con más detalle en las siguientes secciones.
Figura 1. Izquierda: esquema de un ensayo de compresión uniaxial; derecha: Textura 3D de una muestra
de roca, obtenido como una malla 3D. Nota: La resolución es de 1 micrómetro por voxel y la imagen es
de 200x200x200 voxels
BENAVIDES MURILLO: Método EBE de los elementos nitos...
48
2.2 Resolución de la ecuación elastostática
En esta sección se hace una breve presentación del método de los elementos nitos (FEM) en el
contexto de la ecuación elastostática. El objetivo de esta sección es describir algunas ideas básicas,
sucientes para explicar la técnica EBE. Se sigue de cerca la presentación de Hughes (2000), así
como la de Garboczi and Day (1995). Lo aquí expuesto puede utilizarse para una implementación
clásica del FEM.
Quizás la forma más directa de visualizar el FEM es como una técnica proyectiva, usada para
resolver ecuaciones diferenciales y que culmina, en general, con un sistema lineal Kc=f.
Se puede suponer es un campo vectorial con producto interno denido a(
̇
,
̇
).
En el método de los elementos nitos, cada una de las entradas de u, denotada por ui,i=1,2,3 es
aproximada mediante una combinación lineal de funciones de forma
y una función g
i
asociada a las condiciones de contorno del dominio. Si ũ es tal aproximación,
entonces cada una de sus entradas ũ
i
está dada por:
(1)
Esto produce 3N incógnitas c
Ai
. Si e
1
, e
2
, e
3
son los vectores de la base Euclídea en , entonces,
al aplicar la propiedad ortográca a esta proyección, se obtiene:
a(N
A
e
j
,u-ũ)=0 j=1,2,3 (2)
La forma de elegir el operador a(
̇
,
̇
) está asociada con la ecuación diferencial por resolver. En
este caso, se trata de la ecuación elastostática, que en notación inicial tiene la forma:
σ
ij, j
=0 (3)
Aquí, el doble índice denota sumatoria y la coma, diferenciación con respecto a la variable
indizada. Los valores de σ
ij
conforman el tensor de tensiones, simétrico, de tamaño 3×3. Este tipo
de notación es empleada por Hughes (2000), y se caracteriza por ser compacta y breve.
Ahora bien, hallar la solución consiste en obtener el campo de dislocamiento u=(u
1
, u
2
, u
3
) en
todo el dominio. Este campo está relacionado con el tensor de tensiones por medio de la ley de
Hooke (Slaughter, 2002):
σ
ij
=c
ijkl
∈
kl
Donde las constantes c
ijkl
corresponden al tensor de rigidez, dependiente del tipo de material y
∈
kl
es el tensor simétrico de desplazamientos, denido en términos de los desplazamientos u
1
como:
Ingeniería 29 (2): 44-58, julio-diciembre, 2019. ISSN: 2215-2652. San José, Costa Rica. DOI 10.15517/ri.v29i2.36140
49
En materiales isotrópicos el tensor c
ijkl
es simplicado en una matriz 6×6 sacando provecho
de la simetría de los tensores de desplazamiento y el de tensiones. La forma de la Ley de Hooke
para este caso es:
(4)
Donde μ y λ son las constantes de Lamé del material modelado. Para que el problema esté bien
planteado, se deben denir las condiciones de contorno que se clasican en dos tipos: condiciones
de Neuman y condiciones de Drichlet (Hughes, 2000). En mecánica, estas condiciones pueden ser
interpretadas de la siguiente forma:
1. Condiciones de contorno en las que el desplazamiento es conocido T
G
2. Condiciones de contorno en las que es conocida la tensión externa aplicada T
H
En este caso, no existen tensiones externas, pero sí un desplazamiento conocido, tal y como se
corresponde en la prueba de compresión uniaxial. De esta manera, las condiciones de contorno son:
u
i
=q
i
para los desplazamientos conocidos y σ
ij
n
j
= 0 en los puntos de frontera donde no se aplican
presiones externas. El vector n = (n
1
, n
2
, n
3
) es un vector normal supercial. En la expansión de la
ecuación (1), la función g
i
estará asociada a las condiciones de contorno T
G
.
Para la implementación computacional del operador a(
̇
,
̇
), se dene la matriz de derivadas para la
función de forma indizada por A como D
A
:
La matriz de integrales para los índices A y B, K
AB
= se calcula en todo el dominio
de integración. La matriz K
e
está denida por la ecuación (4). El producto interno a(N
A
e
i
,N
B
e
j
),
BENAVIDES MURILLO: Método EBE de los elementos nitos...
50
que depende de los índices A y B de las funciones de forma y de los índices i y j contenidos
en el conjunto {1, 2, 3}, se dene como la entrada (i,j) de la matriz K
AB
. Esto puede escribirse
como . Con este montaje, lo que ocurre es que para cualesquier campo
vectorial W=W
1
e
1
+W
2
e
2
+W
3
e
3
y V=V
1
e
1
+V_2 e
2
+V
3
e
3
se tendrá .
Volviendo a la propiedad ortográca de la ecuación (2), se tiene que:
a(N
A
e
j
, u-ũ ) = a(N
A
e
j
, u) - a(N
A
e
j
, ũ )=0
Donde a(N
A
e
j
, u ) = ʃ (N
A
e
j
)
(p, q)
c
pqkl
u
(k, l)
= ʃ (N
A
e
j
)
(p, q)
σ
p q
Por el Teorema de Green, se obtiene la siguiente separación (Hughes, 2000):
La integral de contorno, por su parte, es evaluada en T=T
H
UT
G
. Suponiendo que las funciones
de forma N
A
son iguales a cero en T
G
y que no existen tensiones externas, el resultado de esta suma
es igual a 0, debido a la ecuación (3). La ecuación (2) se simplica de la siguiente manera:
Utilizando ahora la expansión de la ecuación (1) en este resultado, se obtiene:
(5)
En el FEM, se asume que g es combinación lineal de funciones de soporte compacto (si-
milares a las de forma) que no se anulan en el contorno T
G
. La resolución de este sistema lineal
K
c
=f, donde K es la matriz de productos internos a(N
A
e
j
, N
B
e
i
) y f tiene por entradas los productos
-a(N
A
e
j
,g), es el paso nal del método. K es de tamaño 3N×3N. Las funciones de forma son cero
en la mayoría de los elementos del dominio, de modo que para cada A, existe solo un número
limitado de índices B para los cuales el producto a(N
A
e
j
, N
B
e
i
) no es cero. El número máximo
de índices para el que esto ocurre es el ancho de banda de la matriz de rigidez K. El sistema de
ecuaciones es resuelto mediante la técnica iterativa del gradiente conjugado (Trefethen and Bau,
1997).
Ingeniería 29 (2): 44-58, julio-diciembre, 2019. ISSN: 2215-2652. San José, Costa Rica. DOI 10.15517/ri.v29i2.36140
51
3. DESCRIPCIÓN DEL ALGORITMO EBE
3.1 Implementación del método EBE
Para evitar los costos computacionales de representar la matriz de productos internos
K
ij
=
a(N
A
e
j
, N
B
e
i
)
, se utiliza el método EBE. Esta idea ya fue utilizada en la estimación de los
coecientes elásticos de huesos trabeculares (van Rietbergen et al., 1995). En su investigación, los
autores mencionan que la matriz de productos internos K no precisa ser montada, y que se puede
aprovechar la regularidad del modelo (todos los elementos tienen el mismo tamaño) para ahorrar
memoria y procesamiento. Asumen que las propiedades de la fase sólida de todos los elementos
son idénticos, lo que, como se verá a continuación, no es necesario. Los autores tampoco aplican
ninguna precondición a esta matriz. En el algoritmo que se describe a continuación también se
aplica, de manera opcional, un precondicionador de Jacobi (Trefethen and Bau, 1997).
La estrategia EBE está ligada con el proceso iterativo que obtiene la solución del sistema lineal
K
C
=f. El método de los gradientes conjugados (así como cualquier otra técnica iterativa aplica-
da a este caso), se basa en sucesivos productos de matriz por vector s=Kr que van mejorando la
aproximación de la solución en cada iteración. Aunque en esta técnica numérica se efectúan otros
cálculos vectoriales, este producto matricial es el proceso más pesado computacionalmente. En
lugar de calcular esta matriz de una sola vez, la estrategia EBE toma como entrada el vector r y lo
va multiplicando por los productos internos a(N
A
e
j
, N
A
e
i
) elemento por elemento, actualizando al
mismo tiempo, el vector resultante s.
En otras palabras, la matriz K es ensamblada, por partes, en cada iteración, actuando inme-
diatamente sobre el vector resultante, sin almacenar nunca en memoria todos los valores de la
estructura matricial.
La tarea de calcular los productos internos de las funciones de forma en cada iteración del
método del gradiente conjugado es más optimizable cuando la malla es regular o parcialmente
regular. Es decir, cuando está conformada por un conjunto limitado de clases de elementos. El
almacenamiento en memoria de los productos internos asociados a cada elemento requiere única-
mente una matriz 24×24.
Dicho producto depende de las constantes de Lamé y la geometría del elemento. En física
de rocas, las muestras en general se consideran bifásicas (solo tienen fase sólida y poro), pues el
mayor interés se centra en el estudio del espacio de poro. Los elementos son cubos del mismo
tamaño, uno por cada voxel. Por esa razón, la matriz de productos internos sólo debe ser calculada
una vez. Finalmente, para mejorar la velocidad de convergencia, se utiliza el precondicionador de
Jacobi, de modo que los elementos de la diagonal de K son unitarios. Resumiendo, el método EBE
requiere del cálculo previo de la matriz de productos internos P
k
= a(N
A
e
j
, N
A
e
i
)
ij
para cada tipo de
elemento. Cada una de estas matrices tiene tamaño 24×24. El producto vectorial Kr donde K es
la matriz de rigidez del sistema, se calcula iterando en cada elemento E
w
con base en el siguiente
algoritmo, que aprovecha también la simetría de la matriz de rigidez K:
BENAVIDES MURILLO: Método EBE de los elementos nitos...
52
Entrada: vector r, matrices de productos internos P
w
(para cada elemento E
w
)
Salida: vector s=Kr. Vector precondicionador d.
Inicializar s= 0.
Iterar sobre todos los elementos E
w
Para A=1,2,…,8, índices de los nodos de E
w
Para B=A,…,8, índices de los nodos de E
w
Para i
1
=1,2,3 e i
2
= 1,2,3
Sea p el índice de la la de K correspondiente al nodo A y al grado de libertad
i
1
.
Sea q el índice de la la de K correspondiente al nodo B y al grado de libertad
i
2
.
Sea v=P
w
[3(A-1)+i
1
, 3(B-1)+ i
2
]
s[p]=s[p]+r[q]
̇
v
Si A≠B entonces s[q]=s[q]+r[p]
̇
v
Si p=q se actualiza el precondicionador d[p]=d[p]+v (esto únicamente se
realiza durante la primera iteración, guardando en memoria el precondicio-
nador d).
Multiplicar cada entrada s[i] por 1/d[i] (precondicionador de Jacobi).
3.2 Paralelización
Con el objeto de acelerar las operaciones computacionales, el trabajo que se realiza en cada
elemento en el algoritmo descrito en la sección anterior debe ser ejecutado en paralelo. No es
posible distribuir los elementos entre cada procesador de manera arbitraria, ya que dos elementos
distintos, siendo procesados en procesadores diferentes, pueden intentar actualizar la misma entra-
da del vector resultante s.
Para evitar lo anterior, es necesario separar los elementos en grupos diferentes, de modo que
no exista intersección de los nodos de los elementos de cada grupo. Los grupos se procesan se-
cuencialmente pero cada grupo se procesa en paralelo, pues no existen colisiones. En una malla
regular tridimensional en la que cada voxel es un elemento, son necesarios únicamente 8 grupos
Ingeniería 29 (2): 44-58, julio-diciembre, 2019. ISSN: 2215-2652. San José, Costa Rica. DOI 10.15517/ri.v29i2.36140
53
diferentes. Esta división se lleva a cabo numerando cada voxel de la imagen de manera secuencial.
Para el elemento asociado al voxel j, el grupo al que pertenece está dado simplemente por j (mod
8). Este tipo de paralelización puede ser llevada a cabo tanto en la CPU como en un procesador
gráco (GPU) con suciente memoria.
En la Figura 2 se muestra el esquema de división de grupos y numeración para la paraleliza-
ción. En el caso ilustrado se actualizan los productos internos correspondientes a los elementos
numerados con los valores de j tales que j (mod 8)=1, coloreados en rojo. Los productos internos
de las funciones de forma dentro de esos elementos no actualizan nunca la misma la del vector
resultante, porque los nodos nunca coinciden en cada elemento. Por esta razón, todos esos elemen-
tos en rojo pueden ser procesados en paralelo, sin necesidad de utilizar comandos de exclusión o
MUTEXes (Reif and Rajasekaran, 2008).
Figura 2. Distribución paralela de los elementos durante una iteración
Se muestran tres grupos de 8 voxels donde se realizan los productos internos de las funciones
de forma que no se nulican en los voxels coloreados en rojo. Es de notar que estos voxels colo-
reados no tienen nodos en común, por tanto no actualizan nunca la misma la del vector resultante.
Los otros 7 grupos son trabajados en paralelo en pasos subsecuentes.
Es importante recalcar que la metodología propuesta aquí simplemente posibilita la resolución
de la ecuación elastostática en casos donde el número de elementos es tan grande que no puede ser
administrado por un software FEM comercial. Pero el resultado nal entre montar la matriz K y no
hacerlo debería ser exactamente el mismo.
Por esa razón, una forma de validar esta metodología consiste en comparar el campo de dislo-
camiento obtenido, con el de un software comercial (en una malla pequeña, Figura 3). Una prueba
de este tipo sería suciente para validar la implementación, pero también es posible utilizar una
malla grande (por ejemplo, 400×400×400) sin poros, comparando el módulo de Young obtenido
con el que fue utilizando en cada elemento. El resultado debería ser el mismo.
BENAVIDES MURILLO: Método EBE de los elementos nitos...
54
Figura 3. Izquierda: malla pequeña, porosa y regular; derecha: ensayo de compresión
uniaxial aplicado a la malla, usando FEM
4. VALIDACIÓN
La eciencia y velocidad del truco computacional expuesto en el presente trabajo depende
de detalles de implementación tales como el hardware utilizado, el lenguaje de programación en
el que se desarrolla el método, la forma en que se organiza la disposición de los elementos y los
nodos en la memoria, las particularidades del medio poroso representado.
Se ha pretendido que el algoritmo descrito en la sección 3-1 no dependa de estos factores.
Asimismo, se considera relevante observar que la implementación, realizada en C++ y ejecutada
en una Workstation Dell XPS (con 16 GB de RAM y 8 procesadores en hyperhtread), procesa
imágenes de tamaño 400×400×400 con más de 190 millones de grados de libertad en poco más de
25 minutos.
Dicho experimento sería imposible de resolver utilizando el método FEM representando explí-
citamente la matriz de rigidez sin hacer uso de una memoria RAM de más de 120 GB. Por el con-
trario, con el método EBE es necesario apenas un 5% de esta memoria. No obstante, se considera
que es posible aún mejorar múltiples detalles de implementación y obtener mejores rendimientos.
Las consideraciones antes expuestas, de tipo cualitativo, no deberían ser una guía para la
validación o la conrmación de una implementación especíca, sino apenas una orientación para
quien necesite crear un software de elemento nito con características similares a las expuestas en
el presente trabajo.
Además, es importante notar que lo que se ha descrito aquí es apenas es una técnica de
implementación, cuyos resultados deberían ser totalmente equivalentes a los obtenidos con
un método FEM tradicional, con software comercial o libre, bajo las mismas condiciones. Por
ejemplo, el número de iteraciones del gradiente conjugado usando EBE, es el mismo que se
obtiene con la representación explícita de la matriz, pues la novedad de este algoritmo es que se
está alterando es la forma de representar la matriz de forma y la técnica para calcular el producto
matriz-vector. En general, el método EBE es más lento que el explícito, pero solo pueden realizarse
Ingeniería 29 (2): 44-58, julio-diciembre, 2019. ISSN: 2215-2652. San José, Costa Rica. DOI 10.15517/ri.v29i2.36140
55
tales comparaciones en imágenes de tamaño máximo de 100×100×100 en el computador descrito.
La diferencia en velocidad es aproximadamente de 1 a 5.
En el contexto de la estimación del módulo de Young para medios porosos, es posible validar
la implementación utilizando el trabajo de (Hashin, 1983). Suponiendo una muestra bifásica con
constantes de compresibilidad y elasticidad transversal K
S
, μ
S
para la matriz y K
P
, μ
P
para las in-
clusiones porosas, entonces los módulos de compresibilidad y elasticidad transversal K, μ de la
muestra están dados por:
K=K
S
+pP
SP
(K
P
- K
S
)
μ=μ
S
+ pQ
SP
(μ
P
- μ
S
) (6)
Donde
y p es la porosidad o razón volumétrica de las inclusiones. El módulo de Young equivalente
se obtiene utilizando estas constantes (Slaughter, 2002):
(7)
Estas ecuaciones son válidas bajo la suposición de que las inclusiones porosas dentro de la
muestra son esféricas y están diluidas dentro de la matriz. Esto signica que la distribución de los
poros es uniforme (isotropía) y que la tensión interna en un poro no está inuenciada por los otros
poros. Este es el comportamiento de una matriz porosa innita (Roberts and Garboczi, 2002).
Estas fuertes condiciones pueden ser reproducidas en el espacio discreto de las imágenes,
usando esferas discretas de Bressenham (Bresenham, 1965), distribuidas simétricamente en la ma-
triz (como se ilustra en la Figura 4, izquierda). Los resultados del ensayo de compresión en esta
estructura se ajustan muy bien a las ecuaciones teóricas (6) y (7) (Figura 4, derecha).
Para esta validación se ha utilizado en unidades normalizadas
para las fases sólida y porosa. El error relativo se incrementa con la porosidad, hasta que la hipó-
tesis de poros diluidos de Hashin deja de cumplirse (porosidad mayor a 0.6). Los errores relativos
de los otros casos, para porosidades menores a 0.6, aunque muy pequeños (menores al 1%), no son
nulos. Esto se debe a los errores de discretización del modelado de las inclusiones esféricas en un
espacio basado en voxels.
Dicho esquema de validación es útil cuando es imposible realizar un ensayo de compresión
basado en una implementación comercial tradicional, debido al tamaño y número de las imágenes.
Los resultados mostrados se han ejecutado en 200 imágenes de tamaño 200×200, lo que produce
8 millones de elementos.
BENAVIDES MURILLO: Método EBE de los elementos nitos...
56
Figura 4. Izquierda: distribución uniforme de un arreglo de 27 poros esféricos dentro de una matriz sólida;
derecha: esferas.
Nota: Módulo de Young normalizado obtenido a partir de un ensayo de compresión aplicado
a la muestra, alterando el radio de las esferas para obtenerlo como una función de la porosidad.
Línea continua, módulo de Young teórico, obtenido como función de la porosidad utilizando las
ecuaciones (6) y (7).
5. CONCLUSIONES
Como se ha dicho anteriormente, este tipo de implementación de elementos nitos es útil en
contextos similares a lo que se ha descrito aquí: Física Digital de Rocas para la estimación de pará-
metros de modelos muy heterogéneos, utilizando FEM con muchos elementos iguales y regulares.
Además, se adapta fácilmente para condiciones de contorno periódicas, en donde se supone
que la muestra está sumergida en un medio innito y regular. Su aplicación en computadores de
mayor poder permite resolver imágenes de orden 1000
3
o inclusive2000
3
(en los que el número de
grados de libertad está en el orden de billones).
Sin duda, el método EBE es la única opción viable para implementar la técnica clásica de los
elementos nitos cuando el número de elementos y de grados de libertad producen una matriz de
rigidez que sobrepasa las capacidades físicas de un computador.
La gran cantidad de implementaciones comerciales existentes para resolver problemas numé-
ricos mediante el FEM puede causar la impresión de que no es necesario programar otro tipo de
implementaciones y apenas utilizar las que ya existen, pues de todos modos ya están optimizadas
y depuradas.
En este artículo se ha descrito un contexto en el que una implementación genérica comercial
no es aplicable. Bajo estas circunstancias, un usuario de estas aplicaciones puede verse en una
situación en la que los requerimientos de hardware son imposibles o tienen un costo muy alto, por
lo que el esfuerzo de “reinventar la rueda” y crear una implementación propia tiene una ventaja
importante desde el punto de vista del costo-benecio.
Ingeniería 29 (2): 44-58, julio-diciembre, 2019. ISSN: 2215-2652. San José, Costa Rica. DOI 10.15517/ri.v29i2.36140
57
REFERENCIA BIBLIOGRÁFICAS
Andrä, H., et al. (2013). Digital rock physics benchmarks-Part I: Imaging and segmentation. Comput. Geosci.
50, 25–32. https://doi.org/10.1016/j.cageo.2012.09.005
Arns, C.H., Knackstedt, M.A., Pinczewski, W.V. y Garboczi, E.J. (2002). Computation of linear elastic pro-
perties from microtomographic images: Methodology and agreement between theory and experiment.
Geophysics 67, 1396–1405. https://doi.org/10.1190/1.1512785
Benavides, F., Leiderman, R., Souza, A., Carneiro, G. y Bagueira, R. (2017). Estimating the surface relaxivity
as a function of pore size from NMR T2 distributions and micro-tomographic images. Comput. Geosci.
106. https://doi.org/10.1016/j.cageo.2017.06.016
Betancur, A., Anor, C., Pereira, A. y Leiderman, R. (2018). Determination of the Effective Elastic Modu-
lus for Nodular Cast Iron Using the Boundary Element Method. Metals (Basel). 8, 641. https://doi.
org/10.3390/met8080641
Bresenham, J.E. (1965). Algorithm for computer control of a digital plotter. IBM Syst. J. 4, 25–30. https://
doi.org/10.1147/sj.41.0025
Garboczi, E.J. y Day, A.R. (1995). An algorithm for computing the effective linear elastic properties of hete-
rogeneous materials: Three-dimensional results for composites with equal phase poisson ratios. J. Mech.
Phys. Solids 43, 1349–1362. https://doi.org/10.1016/0022-5096(95)00050-S
Hashin, Z. (1983). Analysis of Composite Materials---A Survey. J. Appl. Mech. Asme 50. https://doi.
org/10.1115/1.3167081
Hudson, J.A. y Harrison, J.P. (2000). Engineering Rock Mechanics: An Introduction to the Principles. Else-
vier Science, 86–96.
Hughes, T.J.R. (2000). The Finite Element Method: Linear Static and Dynamic Finite Element Analysis.
Dover Civil and Mechanical Engineering, 57–107.
Kress, R. (1999). Linear Integral Equations, in: Applied Mathematical Sciences. Springer New York, 218–245.
Leiderman, R., Pereira, A.M.B., Benavides, F.M.J., Silveira, C.S., Almeida, R.M.R. y Bagueira, R.A. (2017).
PERSONAL COMPUTER-BASED DIGITAL PETROPHYSICS. Rev. Bras. Geofísica 35, 95–107.
https://doi.org/10.22564/rbgf.v35i2.891
Makarynska, D., Gurevich, B., Ciz, R., Arns, C.H. y Knackstedt, M.A. (2008). Finite element modelling
of the effective elastic properties of partially saturated rocks. Comput. Geosci. 34, 647–657. https://doi.
org/http://dx.doi.org/10.1016/j.cageo.2007.06.009
Mavko, G., Mukerji, T. y Dvorkin, J. (2003). The Rock Physics Handbook: Tools for Seismic Analysis of
Porous Media, in: Stanford-Cambridge Program. Cambridge University Press, 21–81.
Reif, J. y Rajasekaran, S. (2008). Handbook of Parallel Computing: Models, Algorithms, and Applications,
First Edit. ed. Chapman & Hall/CRC Computer and information science series.
BENAVIDES MURILLO: Método EBE de los elementos nitos...
58
Roberts, A.P. y Garboczi, E. (2000). Elastic Properties of Model Porous Ceramics. J. Am. Ceram. Soc.
Roberts, A.P. y Garboczi, E.J. (2002). Computation of the linear elastic properties of random porous materials
with a wide variety of microstructure. Proc. R. Soc. A Math. Phys. Eng. Sci. 458, 1033–1054. https://
doi.org/10.1098/rspa.2001.0900
Slaughter, W.S. (2002). The Linearized Theory of Elasticity. Birkhäuser Boston, 387–426.
Tiab, D. y Donaldson, E.C. (2012). Petrophysics: Theory and Practice of Measuring Reservoir Rock and
Fluid Transport Properties. EngineeringPro Collection. Elsevier Science, 554–660.
Trefethen, L.N. y Bau, D. (1997). Numerical Linear Algebra. Other Titles in Applied Mathematics. Society
for Industrial and Applied Mathematics (SIAM, 3600 Market Street, Floor 6, Philadelphia, PA 19104),
293–303.
Ulrich, D., Van Rietbergen, B. y Weinans, H., R P. (1998). Finite element analysis of trabecular bone structure:
A comparison of image-based meshing techniques. J. Biomech. 31, 1187–1192. https://doi.org/10.1016/
S0021-9290(98)00118-3
Van Rietbergen, B., Weinans, H., Huiskes, R. y Odgaard, A. (1995). A new method to determine trabecular
bone elastic properties and loading using micromechanical nite-element models. J. Biomech. 28, 69–81.
https://doi.org/http://dx.doi.org/10.1016/0021-9290(95)80008-5
