[19].

${\displaystyle d=1-{\frac {{\sigma }_{f}}{{\sigma }_{I}}}exp\left[A\left(1-{\frac {{\sigma }_{I}}{{\sigma }_{f}}}\right)\right]\,\,con\quad A=}$${\displaystyle {\left[{\frac {G{E}_{0}}{l{\sigma }_{f}^{2}}}-{\frac {1}{2}}\right]}^{-1}}$

(1)

donde ${\textstyle \,l}$ es la longitud de la arista del elemento finito dañado, ${\textstyle {\sigma }_{f}}$ es el límite elástico a tracción del material, ${\textstyle G}$ la energía de fractura y ${\textstyle {E}_{0}}$ el módulo de Young del material no dañado. Resulta obvio que entre dos pasos de tiempo el daño de un lado o una arista no puede disminuir, ya que ello implicaría la reparación espontánea del material, lo cual es termodinámicamente inadmisible.

3.2 Definición del daño elemental

Una vez que el daño ha sido evaluado en cada lado o arista de un elemento (es decir, un enlace virtual entre elementos discretos) el daño sobre el conjunto del elemento se calcula como el máximo daño existente en todos los planos que cortan un triángulo (2D) o un tetraedro (3D) (Fig. 1). En cada tetraedro existen cuatro planos de corte que aíslan un vértice (Fig. 1b) y otros tres que aíslan dos vértices (Fig. 1c). De esta manera, el daño sobre el plano de corte se define como el valor medio del daño de cada arista del elemento que corta dicho plano.

Si el daño elemental sobrepasa un cierto umbral entonces se elimina el triángulo o tetraedro de la malla[6-8] y se crean nuevos elementos discretos en los vértices. El tamaño de estos elementos discretos se define de manera que se conserve la masa del elemento eliminado.

Figura 1. Planos de corte. a) Elemento triangular b) Tetraedro aislando un vértice c) Tetraedro aislando dos vértices.

4.- Discretización mediante DEM

Cuando un tetraedro es totalmente eliminado (i.e. su rigidez es despreciable) se crean 4 nuevos elementos discretos en los vértices de dicho elemento finito. La creación de dichos elementos discretos queda condicionada a que no hayan sido creados con anterioridad debido a la eliminación de algún elemento finito vecino. En este trabajo se han utilizado círculos y esferas para representar los elementos discretos en 2D y 3D, respectivamente. La masa de cada elemento discreto corresponde a la masa nodal que es compatible con la obtenida mediante elementos finitos. El radio de un nuevo elemento discreto esférico se obtiene de tal manera que sea el máximo que garantice el contacto con las esferas vecinas, pero sin crear solapamientos. Existen otros algoritmos para generar elementos discretos[2,20], sin embargo, el propuesto, a pesar de su simplicidad, ha dado buenos resultados. Una vez que se crea un elemento discreto, las fuerzas en las interfaces de contacto se usan para definir la interacción de dicho elemento con los adyacentes. Estas fuerzas se deben únicamente a un mecanismo de contacto en las direcciones normal y tangencial al plano de contacto entre esferas (en 3D), considerando el radio mínimo entre las dos partículas en contacto[5].

En problemas 3D, la fuerza normal de contacto se genera en el punto de contacto entre dos esferas y viene dada por.

${\displaystyle {F}_{n}^{ij}={\frac {{A}^{ij}}{{d}_{ij}}}{E}_{0}{u}_{n}\,\,con\quad {A}^{ij}=}$${\displaystyle \pi {r}_{c}^{2}\quad \,\,(2)}$

(2)

donde ${\textstyle {A}^{ij}}$ es el área de la superficie de contacto entre las dos esferas, ${\textstyle {r}_{c}}$es el radio de la menor de las dos esferas que interactúan en la interfaz i j. ${\textstyle {u}_{n}}$ es el solapamiento en la dirección de los centros de las dos esferas y ${\textstyle {d}_{ij}}$ es la distancia entre dichos centros[5].

La fuerza tangencial ${\textstyle {F}_{s}}$ en el punto de contacto se descompone en dos direcciones ortogonales ${\textstyle {s}_{1}}$ y ${\textstyle {s}_{2}}$ contenidas en el plano normal a la dirección ${\textstyle {d}_{ij}}$. Para cada dirección ${\textstyle {s}_{i}}$ la fuerza tangencial se describe por

${\displaystyle {F}_{s}=min\left\{{\begin{matrix}{u}_{{s}_{i}}{\frac {{E}_{0}}{2(1+\upsilon )}}\\\mu {F}_{n}\left({\frac {{u}_{{s}_{i}}}{{u}_{s}}}\right)\end{matrix}}\right.\quad \quad }$

(3)

donde ${\textstyle \mu }$ es el coeficiente de rozamiento entre esferas y ${\textstyle \upsilon }$ es el coeficiente de Poisson. Más detalles de las ecuaciones (2) y (3) se pueden encontrar en Oñate et al.[5]. En general, el número total de elementos discretos generados en un análisis es solo una pequeña fracción del número total de nodos en la malla de elementos finitos[10]. Por lo tanto, los algoritmos de búsqueda para localizar las fuerzas de contacto entre elementos discretos solo representan un pequeño porcentaje del tiempo total de cálculo.

Es importante mencionar que la estrategia del uso de elementos discretos para definir el contacto entre los labios de una fisura permite resolver problemas estructurales en donde existe apertura y cierre de múltiples grietas, como se muestra en el análisis sísmico de un edificio histórico de mampostería, presentado al final de este artículo.

5.- Integración temporal en sub-pasos

Uno de los principales problemas en el análisis 3D en el tiempo es el gran número de ecuaciones que surgen de la discretización. Sumando a estos la dificultad para considerar un gran período de tiempo, se concluye que en problemas dinámicos la integración implícita en el tiempo es la mejor estrategia de solución. Sin embargo, este tipo de integraciones son complejas de utilizar en el DEM debido a que es imposible identificar y cuantificar correctamente los contactos y las fuerzas de contacto entreelementos discretos.

La implementación para el estudio de problemas transitorios seguida en esta investigación corresponde a un esquema de integración temporal de sub-etapas en el cual los elementos finitos se calculan mediante un esquema implícito de Newmark y los elementos discretos mediante un esquema explícito de diferencias centradas. Normalmente el incremento de tiempo del esquema implícito suele ser 100 veces mayor que el explícito.

La ventaja de esta estrategia es en que el número de elementos discretos suele ser mucho menor que el de elementos finitos que discretizan la estructura, por lo que la integración explícita realizada dentro de un paso de tiempo implícito es bastante rápida. Antes de comenzar un nuevo paso de tiempo ${\textstyle {\Delta }_{I}t}$ en el esquema implícito, se calcula el mismo período de tiempo con un esquema explícito sobre los elementos discretos usando un incremento de tiempo ${\textstyle {\Delta }_{E}t}$.

El contacto entre elementos discretos se cuantifica por la suma de los impulsos que se producen a lo largo del análisis explícito y se expresa como una fuerza sobre los elementos discretos, aplicada en el tiempo ${\displaystyle t+1}$ correspondiente al esquema implícito tal como se muestra en la Figura 2.

Figura2 Integración temporal en sub-etapas utilizado para el análisis conjunto de elementos finitos y discretos.

Se ha observado que no es aconsejable que los intervalos de tiempo ${\textstyle {\Delta }_{I}t}$ y ${\textstyle {\Delta }_{E}t}$ tengan una relación mayor que 1: 500 ya que pueden haber discrepancias entre la solución explícita y la implícita. Teniendo esto en cuenta, la estrategia de solución en sub-pasos implementada permite obtener excelentes resultados, como se muestra en los ejemplos del siguiente apartado.

6.- Ejemplos

En este apartado se muestran varios ejemplos a fin de mostrar el buen comportamiento de la estrategia FEM-DEM descrita. El primer ejemplo corresponde al estudio 3D de una probeta normalizada en un ensayo a tracción. El segundo ejemplo es una viga de hormigón bi-entallada donde predomina la fractura en modo mixto. El tercer ejemplo corresponde al ensayo de tracción indirecta, ampliamente usado en mecánica de rocas. El cuarto ejemplo consiste en un ensayo de cortante en hormigón propuesto por Luong[21]. El quinto ejemplo es un ensayo de una probeta a compresión simple. El sexto ejemplo es el estudio de un forjado reticular de hormigón armado, afectado por un asentamiento diferencial. El acero se considera embebido y solidario con los desplazamientos de los elementos finitos utilizados.

Finalmente, se presenta el análisis sísmico de la nave central de la iglesia del monasterio de Poblet , formada por elementos de mampostería, aplicándole un sismo de ${\textstyle 6{M}_{w}}$. En este ejemplo se puede observar como la estrategia FEM-DEM modela el efecto de la apertura y cierre de múltiples grietas.

6.1 Ensayo normalizado de tracción

El primer ejemplo corresponde al análisis de fractura de una probeta de hormigón sujeta a fuerzas de tracción. El objetivo principal es mostrar la independencia del tamaño de la malla en la generación de la grieta, de manera que la energía utilizada en la fractura sea independiente del tamaño elemental. La geometría se define de acuerdo con la norma D638 de la Sección Norteamericana de la Asociación Internacional de Ensayo de Materiales (ASTM, American Society for Testing and Materials) [22]. En la Figura 3 se muestran las tres mallas de elementos tetraédricos de 4 nodos utilizadas, así como las condiciones de contorno. La probeta se carga imponiendo un campo de velocidad constante de tracción en sus extremos, representados por la zona oscura en la figura.

Figura 3 Ensayo normalizado 3D de tracción. Malla de elementos finitos y dimensiones según la norma ASTMD638.

El estudio se ha realizado utilizando la metodología FEM-DEM en 3D antes descrita. Con el fin de localizar la fractura, sólo se permite que una banda de elementos se rompa al nivel de la tensión de fallo, usando el modelo de daño mencionado. Los resultados obtenidos se analizan dibujando los desplazamientos de los puntos PA y PB mostrados en la Figura 3.

El módulo de Young, el coeficiente de Poisson y la densidad son respectivamente ${\textstyle {E}_{0}=30x10^{9}Pa,\nu =0.2,\gamma =1.0x10^{3}N/m^{3}}$, la tensión máxima de tracción ${\textstyle {\sigma }_{f}=10x10^{3}Pa}$ y la energía de fractura ${\displaystyle G=7.5x10^{-3}J/m^{2}}$.

La probeta se deforma aplicando una velocidad constante de tracción de ${\displaystyle 0.5x10^{-7}}$ m/s en ambos extremos. La Figura 4 muestra la geometría dañada para las tres mallas del FEM consideradas.

Figure 4: Ensayo normalizado de tracción. Zona agrietada incluyendo los elementos discretos generados para las tres mallas de elementos finitos consideradas.

Con objeto de evaluar la apertura de la grieta se analiza el desplazamiento de los puntos PA y PB situados a la derecha y en el centro de la probeta, respectivamente (Figura 3). La Figura 5 muestra la relación carga-desplazamiento en estos puntos. Para las tres mallas consideradas la evolución del desplazamiento es muy similar y de acuerdo con los resultados esperados[9]. Debido a que los elementos por donde se abre la grieta tienen un tamaño diferente para cada malla, el desplazamiento del punto PB en la región elástica se hace más pequeño a medida que se reduce el tamaño del elemento.

Figure 5: Ensayo normalizado 3D de tracción. Curva carga-desplazamiento de los puntos PA y PB de la muestra.

6.2 Viga bientallada a flexión.

El siguiente ejemplo corresponde al ensayo de una viga de hormigón en masa bi-entallada. El análisis se realiza mediante las hipótesis de tensión plana y es un buen ejemplo de fractura en modo mixto. La viga se sostiene en dos puntos y se somete a flexión aplicando un desplazamiento impuesto mediante control de velocidad equivalente a ${\displaystyle 1mm/s}$ en los dos puntos representados en la Figura 6. Igualmente, en dicha figura también se muestran la geometría y las dimensiones de la probeta.

Figure 6: Viga bi-entallada a flexión. Geometría y condiciones de contorno.

La viga presenta dos puntos singulares en la punta de las dos entallas en donde las tensiones de tracción son altas y el daño comienza en esta zona. Las propiedades del material son ${\textstyle {E}_{0}=30x10^{9}Pa,\nu =0.2,{\sigma }_{f}=2MPayG=1x10^{2}J/m^{2}}$. El problema ha sido resuelto con la técnica FEM-DEM en 2D.

En la Figura 7 se muestra un detalle de las tres diferentes mallas utilizadas, formadas por ${\displaystyle 1165}$ nodos y ${\displaystyle 2202}$ elementos triangulares lineales para la malla gruesa, ${\displaystyle 1847}$ nodos y ${\displaystyle 3480}$ elementos para la malla intermedia y ${\displaystyle 5747}$ nodos y ${\displaystyle 11206}$ elementos para la malla fina. El análisis se ha realizado tanto forma cuasiestática como de forma dinámica, respetando la velocidad de aplicación de la carga. En ambos casos los resultados han sido muy similares como se observa en la Figura 9.

Figure 7: Viga bi-entallada a flexión. Mallas utilizadas.

La Figura 8 muestra la dirección de las fisuras para las tres mallas analizadas las cuales coinciden con los experimentos numéricos [14]. La Figura 9 muestra la relación entre la reacción y el desplazamiento impuesto en cualquiera de los dos puntos representados en la figura 6 (los resultados son idénticos para los dos puntos). Los gráficos son concordantes con los resultados obtenidos por Cervera et. al. [14].

Figure 8:Viga bientallada a flexión. Relación entre la fuerza y el desplazamiento en cualquiera de los puntos referenciados en la Figura 6. Los resultados obtenidos son comparados con los aquellos dados en Cervera et al. [14]

Figure 9: Viga bi-entallada a flexión. Superficies de igual desplazamiento y apertura de las grietas. a) Malla gruesa, b) malla intermedia, c) Malla fina.

6.3 Ensayo de tracción indirecta

El ensayo brasileño de tracción indirecta (BTS, de Brazilian Tensile Strength) es un procedimiento sencillo para evaluar la resistencia a la tracción de hormigón y geomateriales. La probeta de hormigón analizada es un cilindro de ${\displaystyle 0.2m}$ de diámetro (D) y ${\displaystyle 0.1m}$ de espesor (t), sujeto a una carga diametralmente opuesta (Figura 10). El valor de la resistencia a la tracción se calcula mediante la siguiente expresión[23,24].

Figure 10: Ensayo de tracción indirecta. Dimensiones de la muestra, condiciones de contorno y mallas usadas de elementos finitos tetraédricos de cuatro nodos

${\displaystyle {\sigma }_{f}={\frac {2P}{\pi tD}}\quad \,\,}$

(4)

Donde P es el valor de la carga aplicada. Las propiedades del material son ${\textstyle {E}_{0}=21x10^{9}Pa,\nu =0.2,\gamma =7.8x10^{3}N/m^{3},{\sigma }_{f}=10KPa}$ y ${\displaystyle G=1x10^{2}J/m^{2}}$. Lo que proporciona una carga máxima de fallo de ${\displaystyle P=314.16N}$.

Para realizar los análisis se han usado tres mallas de ${\displaystyle 9338,31455}$ y ${\displaystyle 61623}$ elementos tetraédricos lineales como se muestra en la Figura 10. El ensayo se realiza imponiendo una velocidad constante vertical en la parte superior de la probeta.

La Figura 11 muestra la grieta y los elementos discretos generados. Se puede observar que el patrón de fisuración es similar para las tres mallas y de acuerdo con el resultado esperado.

Figure 11: Ensayo de tracción indirecta. Grieta y elementos discretos generados. a) Malla gruesa, b) malla intermedia, c) malla fina.

La Figura 12 muestra la curva de carga-desplazamiento. Los valores obtenidos para la resistencia a la tracción máxima para las mallas gruesa, media y fina son respectivamente: ${\displaystyle 10.693KPa,10.351KPa}$ y ${\displaystyle 10.235KPa}$, correspondientes a un rango de entre ${\displaystyle 6\%}$ y ${\displaystyle 2\%}$ de error frente al valor esperado de ${\textstyle {\sigma }_{f}=10KPa.}$

Figure 12: Ensayo de tracción indirecta. Relación fuerza-desplazamiento para las tres mallas utilizadas.

Es remarcable la insensibilidad de la curva carga-desplazamiento al tipo de malla de elementos finitos utilizada. Esta “objetividad” de los resultados numéricos frente al tamaño de la malla del FEM es una de las características esenciales de la técnica FEM-DEM.

6.4 Ensayo de cortante

El ensayo de cortante se diseña para aplicar un esfuerzo de cizalladura sobre una probeta de modo que experimente una falla por deslizamiento a lo largo de un plano paralelo a las fuerzas aplicadas. Por lo general las fuerzas de cizallamiento provocan que una de las superficies de fallo de un material se mueva en una dirección y la otra superficie en dirección opuesta, de manera que el material se encuentra sometido a un estado de corte. Este ejemplo se realizó en 80 minutos usando un procesador a ${\displaystyle 2.5MHz}$.

El ensayo que aquí se considera tiene como objetivo determinar la resistencia a cortante del hormigón y ha sido propuesto por Luong[21]. La probeta tiene forma tubular y su eje coincide con el eje ${\displaystyle z=0}$. Tiene varias entallas y se somete a una carga central en una de sus caras y otra excéntrica en la cara opuesta, de manera que se generen tensiones cortantes paralelas al eje ${\displaystyle z=0}$, como se describe en la Figura 13 en la que también se muestran las condiciones de contorno impuestas. La profundidad de la entalla es de 10mm y su ancho es de 4mm. La carga se aplica imponiendo a la placa superior una velocidad constante de 1mm/s hasta llegar a la fractura.

Figure 13: Ensayo de cortante. Geometría, cargas y condiciones de contorno.

La definición geométrica y la malla de elementos finitos utilizada se muestran en la Figura 14. El plano de cizallamiento ha sido discretizado con 4 elementos finitos tetraédricos de 4 nodos con el fin de captar adecuadamente el gradiente de tensiones en esta zona.

Figure 14: Ensayo de cortante. a) Definición geométrica de la probeta b) Malla de elementos finitos.

Las propiedades del hormigón son ${\textstyle {E}_{0}=35x10^{9}Pa,\nu =0.22,\gamma =7.8x10^{3}N/m^{3},{\sigma }_{c}=30MPa,{\sigma }_{f}=3MPa}$ y ${\displaystyle G=75x10^{-3}J/m^{2}}$. En general el esfuerzo máximo de corte en el hormigón es aproximadamente 1/5 o 1/6 del esfuerzo a compresión