Numerical, theoretical and experimental validation of the material point method to solve geotechnical engineering problems



<span id='centerContent'></span>

<span id='centerInner'></span>
<span id='tit0005'></span>
=Validación numérica, teórica y experimental del método del punto material para resolver problemas geotécnicos=

Numerical, theoretical and experimental validation of the material point method to solve geotechnical engineering problems

M.A. Llano-Serna<sup>, </sup>, M.M. Farias

<ul><li><span id='aff0005'></span>
Facultad de Ingeniería Civil y Ambiental, Universidad de Brasilia, Brasilia, Brasil</li>
</ul>

Under a Creative Commons [http://creativecommons.org/licenses/by-nc-nd/4.0/ license]

<span id='ppvPlaceHolder'></span>

<span id='refersToAndreferredToBy'></span>
<span id='referredToBy'></span>

<span id='authorabs00101'></span>
==Resumen==

<span id='spar0005'></span>
Este trabajo investiga la posibilidad de utilizar el método del punto material (MPM) para resolver problemas geotécnicos cuasiestáticos de pequeñas deformaciones y problemas dinámicos que presentan grandes distorsiones. Métodos tradicionales como el método de los elementos finitos (MEF) tienen muchas dificultades para resolver los problemas de las grandes deformaciones. Por este motivo, herramientas como el MPM han adquirido importancia en los últimos años. Como cualquier herramienta nueva, es conveniente verificar su eficacia en la simulación de problemas propios de la ingeniería geotécnica. Primero se describe de forma breve la formulación matemática del MPM, y a continuación se presentan algunas simulaciones realizadas con el MPM de una cimentación superficial, un ensayo de compresión no confinada y deformaciones en un talud. A tal efecto, se ha utilizado un software de código abierto y los resultados obtenidos se han comparado con modelos de MEF, soluciones analíticas y ensayos de laboratorio reales. Los resultados muestran una concordancia cualitativa y cuantitativa, y un mejor desempeño en la simulación de las tensiones que en la de las deformaciones. El conjunto de los modelos evaluados demuestra que el MPM es una herramienta válida y útil para simular problemas propios de geotecnia de pequeñas y grandes deformaciones, aunque el MEF tradicional sigue siendo computacionalmente más eficiente para resolver problemas cuasiestáticos.

<span id='authorabs00052'></span>
==Abstract==

<span id='spar0010'></span>
This paper investigates the possibility of using the material point method (MPM) to solve small strains quasi-static problems and dynamic problems related to large distortions. Traditional methods such as the finite element method (FEM) face difficulties when large strains are involved. Therefore, tools such as the MPM have become more important in recent years. As a new tool, the MPM needs to prove its functionality for geotechnical engineering problems. In first place the MPM mathematical formulation is shortly described. Next, numerical simulations of a shallow foundation, an unconfined compression test and a slope problem are performed in an open source MPM code. The results are compared with FEM simulations, analytical solutions and real laboratory tests. The study shows qualitative and quantitative agreement when compared; a better performance of MPM for solving stresses better than strains is detected. The set of simulations validates the MPM to solve geotechnical engineering problems when dealing with small and large strains. However, the traditional FEM showed a better performance for quasi-static cases.

<span id='kwd_1'></span>
==Palabras clave==

Método del punto material; Grandes deformaciones; Compresión no confinada de arcillas; Cimentación superficial sobre suelos blandos; Deformaciones en taludes de suelos cohesivos

<span id='kwd_2'></span>
==Keywords==

Material point method; Large strains; Unconfined compression of clay; Shallow foundations in soft soil; Slope deformations in cohesive soils

<span id='embeddedPDF'></span>

<span id='SD_BA1P'></span>

<span id='sec0005'></span>
==1. Introducción==

<span id='par0005'></span>
Debido al comportamiento complejo de los suelos, los ingenieros de hoy en día se ven en la necesidad de utilizar metodologías prácticas para diseñar obras cada vez más complejas al servicio de la sociedad. En las últimas décadas se han desarrollado muchas herramientas numéricas basadas en el método de los elementos finitos (MEF), debido a su versatilidad [[#bib0150|[1]]]. En la literatura pueden encontrarse aplicaciones diversas en el campo de la ingeniería geotécnica.

<span id='par0010'></span>
A pesar del uso extendido del MEF, presenta algunas limitaciones evidentes al intentar resolver problemas de grandes deformaciones. Para extender el uso de este método a la simulación de grandes deformaciones, manteniendo su enfoque tradicional, es necesario utilizar tensores de grandes deformaciones y actualizar la configuración del objeto de análisis. Si la malla utilizada para discretizar el medio experimenta una gran distorsión, tiene que ser actualizada, lo cual genera imprecisiones numéricas a consecuencia de la pérdida de información cuando tiene que mapearse una variable de estado al nuevo esquema de discretización [[#bib0155|[2]]].

<span id='par0015'></span>
Por ello, las investigaciones más recientes se han centrado en mejorar la simulación de los problemas que implican una gran distorsión. Este tipo de problemas son muy relevantes cuando intervienen fenómenos de impacto, penetración, fragmentación e interacción entre diferentes materiales [[#bib0160|[3]]].

<span id='par0020'></span>
Entre las nuevas técnicas para resolver grandes deformaciones, destaca el método del punto material (MPM), por las posibilidades que ofrece. El MPM es una extensión del método ''particle in cell'' (PIC), muy utilizado en mecánica computacional de fluidos. Este método puede describirse como la unión de 2 esquemas de discretización: una formulación lagrangiana para representar el cuerpo continuo como un conjunto de puntos materiales y una formulación euleriana para resolver las ecuaciones de campo mediante una malla de cálculo.

<span id='sec0010'></span>
==2. Antecedentes del método del punto material==

<span id='par0025'></span>
Según Beuth et al. [[#bib0165|[4]]], el MPM surgió en la década de 1960 en el Laboratorio Nacional de Los Álamos durante el proceso de búsqueda de una solución a problemas complejos de mecánica de fluidos, mediante una representación de puntos o partículas que se movían a través de una malla fija [[#bib0170|[5]]].

<span id='par0030'></span>
La trayectoria de este conjunto de partículas o puntos materiales (PM) que componen el volumen de análisis normalmente queda registrada. Los PM almacenan toda la información relevante: velocidades, deformaciones, tensiones, etc. La información de los PM es extrapolada o transferida periódicamente a la malla de fondo (similar a una malla de MEF). Las ecuaciones de movimiento (conservación de momento linear) se resuelven de manera explícita, para un intervalo de tiempo finito, utilizando el método de las diferencias finitas (MDF). Una vez resuelto el sistema, las velocidades en la malla de fondo se interpolan de nuevo a los PM para así poder determinar su posición final.

<span id='par0035'></span>
Por último, se calculan las deformaciones en los PM y se actualizan las tensiones mediante una relación constitutiva. En la [[#fig0005|figura 1]] se ilustra el esquema de discretización. Una vez transferida la información a los PM, se descartan los valores de la malla de fondo y el proceso se repite para un nuevo intervalo de tiempo; de esta forma, la solución evoluciona en el tiempo hasta lograr una estabilización total o llegar al punto de análisis deseado. Durante este proceso, los PM pueden circular a través de diferentes células de la malla de fondo. Se observó la presencia de ruido numérico en el MPM cuando los PM cruzaban las fronteras de las células que los contenían, situación que fue resuelta utilizando funciones de interpolación suavizadas [[#bib0175|[6]]].

<span id='figure_fig0005'></span>

<span id='fig0005'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_799892250-1-s2.0-S0213131515000292-gr1.jpg|center|349px|Discretización utilizada en el MPM.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 1.

<span id='spar0015'></span>
Discretización utilizada en el MPM.

</span>
|}

<span id='par0040'></span>
En la literatura existen algunos estudios que tratan de la validación y la aplicación del MPM en la geotecnia, referidos a modelos simples, como vigas en voladizo [[#bib0165|[4]]] and [[#bib0180|[7]]]; a aspectos relacionados con la consolidación [[#bib0165|[4]]] and [[#bib0185|[8]]]; a asentamientos en taludes [[#bib0190|[9]]] and [[#bib0195|[10]]]; a cimentaciones superficiales [[#bib0150|[1]]], [[#bib0200|[11]]] and [[#bib0205|[12]]]; a problemas de penetración [[#bib0205|[12]]] and [[#bib0210|[13]]]; a empujes activos y pasivos [[#bib0150|[1]]], [[#bib0180|[7]]], [[#bib0195|[10]]], [[#bib0200|[11]]], [[#bib0215|[14]]] and [[#bib0220|[15]]]; a discontinuidades, fisuras y grietas [[#bib0225|[16]]] and [[#bib0230|[17]]]; a asentamientos en terraplenes [[#bib0235|[18]]]; a la estabilidad de los taludes y las laderas [[#bib0240|[19]]], [[#bib0245|[20]]] and [[#bib0250|[21]]]; a cimentaciones sobre suelos blandos [[#bib0185|[8]]], y a problemas de penetración [[#bib0215|[14]]] and [[#bib0255|[22]]]. Sin embargo, cabe destacar que no abundan las validaciones rigurosas sobre los resultados obtenidos.

<span id='sec0015'></span>
==3. El método del punto material==

<span id='par0045'></span>
Es posible encontrar en la literatura numerosas formulaciones que describen el MPM. En este trabajo, nos remitimos a la formulación presentada por Buzzi et al. [[#bib0260|[23]]]. La notación tensorial adoptada es la siguiente: el tensor es representado por una letra y su orden es marcado por el número de virgulillas (∼) bajo la letra; un punto en negrita (·) simboliza una contracción simple; los dos puntos (:) significan una contracción doble, y el símbolo (⊗) es el producto diádico o tensorial.

<span id='sec0020'></span>
===3.1. Formulación===

<span id='par0050'></span>
La ecuación de campo que describe el movimiento dinámico de puntos en un medio continuo ha de satisfacer las condiciones de equilibrio, las cuales pueden ser representadas mediante notación tensorial de la manera siguiente:

<span id='eq0005'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{d\underset{\sim}{\underset{\sim}{\sigma }}}{d\underset{\sim}{x}}:\underset{\sim}{\underset{\sim}{I}}+\rho \underset{\sim}{b}=\rho \underset{\sim}{a}</math>
|}
| style="text-align: right;" | (1)
|}

donde <math display="inline">\underset{\sim}{\underset{\sim}{\sigma }}</math> es el tensor de tensiones de segundo orden; <math display="inline">\underset{\sim}{x}</math> representa el vector de posición de un punto, <math display="inline">\underset{\sim}{\underset{\sim}{I}}</math> expresa el tensor de identidad de segundo orden; ''ρ  '' simboliza el campo de densidad escalar; <math display="inline">\underset{\sim}{b}</math> son las fuerzas de cuerpo y externas que actúan en un punto, y <math display="inline">\underset{\sim}{a}</math> es el vector de aceleración en un punto. La ecuación (1) es simplemente la representación de la segunda ley de Newton para los puntos de un medio continuo.

<span id='par0055'></span>
La forma discreta débil de la ecuación (1) sobre las células de la malla con volumen (''V  '') se obtiene utilizando el método de los residuos ponderados, mediante funciones de peso <math display="inline">\underset{\sim}{w}\left(\underset{\sim}{x}\right)</math> similar al MEF, que toma la forma siguiente, después de aplicar el teorema de la divergencia:

<span id='eq0010'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\int_A\underset{\sim}{w}\cdot \underset{\sim}{t}dA+\int_V\rho \underset{\sim}{w}\cdot \underset{\sim}{b}dV-\int_V\frac{d\underset{\sim}{w}}{d\underset{\sim}{x}}:\underset{\sim}{\underset{\sim}{\sigma }}dV=\int_V\rho \underset{\sim}{w}\cdot \underset{\sim}{a}dV</math>
|}
| style="text-align: right;" | (2)
|}

donde <math display="inline">\underset{\sim}{t}</math> representa las fuerzas de tracción en la superficie en una región ''A''.

<span id='par0060'></span>
Los PM ocupan solo una fracción del volumen de la célula, normalmente entre 1/4 y 1/9 del volumen total de la célula (''V  ''), y los PM pueden moverse libremente entre las células, por lo cual su dominio debe registrarse. Esto se logra mediante una función característica en la partícula <math display="inline">\chi _p\left(\underset{\sim}{x}\right)</math> que toma el valor de 1 cuando se encuentra dentro del dominio del punto material y de 0 cuando se halla fuera de sus límites. Utilizando este artificio matemático, y tras realizar algunas manipulaciones matemáticas, la ecuación (2) puede transformarse en:

<span id='eq0015'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\underset{\underset{\sim}{f}_n^{ext}}{\underbrace{\int_AS_n\underset{\sim}{t}dA+\sum_pm_p\underset{\sim}{b}\bar{S}_{np}}}-\underset{\underset{\sim}{f}_n^{int}}{\underbrace{\sum_pV_p\underset{\sim}{\underset{\sim}{\sigma }}_p\cdot \underset{\sim}{\bar{G}}_{np}}}=\underset{\underset{\sim}{\dot{q}}_n}{\underbrace{\sum_p\dot{\underset{\sim}{q}}_p\bar{S}_{np}}}</math>
|}
| style="text-align: right;" | (3)
|}

<span id='par0065'></span>
La ecuación (3) establece que la tasa de cambio del momento lineal en los vértices de la malla <math display="inline">\underset{\sim}{\dot{q}}_n</math> es igual a la resultante de la sumatoria de las fuerzas externas e internas, es decir, se cumple la ley de conservación del momento lineal.

<span id='par0070'></span>
Las fuerzas externas en los vértices de la malla <math display="inline">\underset{\sim}{f}_n^{ext}</math> reciben las contribuciones de las tracciones en la superficie <math display="inline">\left(\underset{\sim}{t}\right)</math> y de las fuerzas de cuerpo <math display="inline">\left(\underset{\sim}{b}\right)</math> aplicadas en el medio continuo. De forma similar al MEF, las tracciones en la superficie son transferidas directamente a los vértices, multiplicándolas por las funciones de forma ''S''<sub>''n''</sub> e integrándolas sobre el área de la superficie, lo cual se representa por el primer término de la ecuación (3).

<span id='par0075'></span>
Las fuerzas de cuerpo también han de transponerse de los PM internos de las células a los vértices de la malla, lo cual se logra mediante las funciones de forma y posterior integración sobre el dominio:

<span id='eq0020'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\bar{S}_{np}(\underset{\sim}{x})=\frac{1}{V_p}\int_{V\ast }S_n(\underset{\sim}{x})\chi _p(\underset{\sim}{x})dV</math>
|}
| style="text-align: right;" | (4)
|}

<span id='par0080'></span>
Por otro lado, la transferencia de tensiones se logra multiplicándolas por el gradiente de la función de forma <math display="inline">\underset{\sim}{G}_n=dS_n(\underset{\sim}{x})/d\underset{\sim}{x}</math>, que solo deberá calcularse dentro del dominio del PM dado por la función <math display="inline">\chi _p\left(\underset{\sim}{x}\right)</math>. De esta manera, la extrapolación se obtiene del <math display="inline">\underset{\sim}{\bar{G}}_{np}</math>, calculado de la manera siguiente:

<span id='eq0025'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\underset{\sim}{\bar{G}}_{np}(\underset{\sim}{x})=\frac{1}{V_p}\int_{V\ast }\frac{dS_n(\underset{\sim}{x})}{d\underset{\sim}{x}}\chi _p(\underset{\sim}{x})dV</math>
|}
| style="text-align: right;" | (5)
|}

<span id='par0085'></span>
Las integrales de las ecuaciones (4) y (5) pueden calcularse analíticamente, puesto que su forma es sencilla y, por tanto, no es necesario utilizar un esquema de integración numérica. En el trabajo de Bardenhagen y Kober [[#bib0175|[6]]] pueden encontrarse expresiones explícitas para estas matrices, y el ciclo completo de la solución es descrito en detalle por Buzzi et al. [[#bib0260|[23]]].

<span id='sec0025'></span>
==4. Simulaciones numéricas==

<span id='par0090'></span>
Las simulaciones presentadas fueron ejecutadas en condiciones bidimensionales, en estado plano de tensiones y utilizando el código abierto NairnMPM 8.1.0 [[#bib0265|[24]]]. Este código ha sido verificado para algunos problemas de madera. Debido al uso extendido del MEF, es natural comparar los resultados obtenidos en este trabajo con este método, como también es válido el uso de ensayos de laboratorio y soluciones analíticas. Para las simulaciones con el MEF, se utilizó el código comercial Plaxis.

<span id='sec0030'></span>
===4.1. Ensayo de compresión no confinada===

<span id='par0095'></span>
Uno de los problemas que permite evaluar, de forma sencilla y efectiva, la eficacia del MPM para grandes deformaciones es el ensayo de compresión no confinada o compresión simple [[#bib0270|[25]]]. El material utilizado para ejecutar el ensayo fue caolín manufacturado. Con el objetivo de alcanzar una gran distorsión, la muestra de suelo fue llevada hasta una deformación del 35 %, midiendo la tensión principal y la deformación en el eje vertical.

<span id='par0100'></span>
Una vez ejecutado el ensayo, se seleccionó un modelo de ruptura simple del tipo Von Mises, dado que el ensayo de compresión no confinada solo supone la medición de las tensiones totales y la muestra de suelo fue sometida a condiciones saturadas no drenadas (carga rápida). Además, se supuso un comportamiento elástico en la primera etapa de la carga, asumiendo que el suelo tiene 2 tipos de comportamiento: primero perfectamente elástico y, una vez lograda la ruptura, perfectamente plástico.

<span id='par0105'></span>
El módulo de elasticidad secante calculado suponiendo un 3% de deformación fue ''E'' = 405 kPa; la tensión de plastificación, ''σ''<sub>''y''</sub> = 52 kPa, y la densidad, ''γ'' = 18 kN/m<sup>3</sup>. Y debido a que la muestra se encontraba en condiciones saturadas, se asumió ''ν'' = 0,5. En la [[#fig0010|figura 2]] se presentan las etapas de deformación simuladas hasta alcanzar una deformación de un 35% con el MPM (lado izquierdo de cada etapa). En el caso del MEF (lado derecho), el problema solo se consiguió simular hasta una deformación máxima del 23% —las diferencias se acentúan en la última etapa del esquema. Es importante señalar que la simulación con el MPM puede reproducir mayores distorsiones y que la simulación que aquí se presenta se detuvo expresamente una vez superada la frontera establecida por el MEF (Plaxis).

<span id='figure_fig0010'></span>

<span id='fig0010'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_799892250-1-s2.0-S0213131515000292-gr2.jpg|center|357px|Etapas de deformación en la simulación del ensayo de compresión no confinada ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 2.

<span id='spar0020'></span>
Etapas de deformación en la simulación del ensayo de compresión no confinada utilizando MPM.

</span>
|}

<span id='par0110'></span>
En la [[#tbl0005|tabla 1]] se enumeran a título comparativo las características geométricas y numéricas de las simulaciones realizadas con ambos métodos numéricos. También se indica el tiempo computacional que invirtió cada uno de los modelos: el MPM tardó un 30% más tiempo en lograr el mismo nivel de deformación que el MEF y 151 s en total en alcanzar la máxima deformación registrada.

<span id='table_tbl0005'></span>

<span id='tbl0005'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabla 1.

Características generales de la simulación del ensayo de compresión no confinada

|-

! 
! MPM
! MEF
|-

| Elementos
| 238 cuadrados
| 226 triangulares (15 nodos)
|-

| Discretización
| 324 puntos materiales
| 2.711 puntos de Gauss
|-

| Tiempo computacional (s)
| 113
| 88
|}

<span id='par0115'></span>
En la [[#fig0015|figura 3]] se indican los resultados de las curvas obtenidas en el laboratorio, la idealización teórica y el modelo numérico en el MPM. Es posible ver que, en el intervalo de deformaciones elásticas hasta alcanzar una deformación del 10%, los resultados de la simulación tienen una tendencia muy similar a la teórica. Sin embargo, a partir de este punto se observa que la solución numérica experimenta oscilaciones hasta alcanzar la tensión de plastificación, fenómeno conocido como ''kinematic locking''[[#bib0275|[26]]]. No obstante, los resultados de la simulación alcanzan la tensión de ruptura con una diferencia del 3% con respecto al ensayo real.

<span id='figure_fig0015'></span>

<span id='fig0015'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_799892250-1-s2.0-S0213131515000292-gr3.jpg|center|345px|Comparación de las curvas tensión vertical-deformación del ensayo de compresión ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 3.

<span id='spar0025'></span>
Comparación de las curvas tensión vertical-deformación del ensayo de compresión no confinada.

</span>
|}

<span id='par0120'></span>
Al comparar la simulación con los datos experimentales puede observarse que el ensayo en el laboratorio presenta una transición más suave ya que, de hecho, muestra un comportamiento elastoplástico desde el inicio de la deformación. A pesar de las diferencias en la transición elastoplástica, las simulaciones numéricas presentaron un excelente desempeño en el tramo inicial (intervalo elástico) y el tramo final (intervalo plástico). Coetzee [[#bib0150|[1]]] y Konagai y Johansson [[#bib0280|[27]]] obtuvieron unos resultados similares.

<span id='sec0035'></span>
===4.2. Cimentación superficial sobre un estrato de suelo finito===

<span id='par0125'></span>
Uno de los problemas clásicos de la ingeniería geotécnica son las cimentaciones superficiales, las cuales se pueden representar como una carga distribuida uniforme sobre un estrato de suelo finito. Los resultados de las simulaciones utilizando el MPM se compararon con soluciones analíticas y modelos numéricos utilizando el MEF con el programa Plaxis.

<span id='par0130'></span>
La solución analítica a las tensiones considerando una interfaz rugosa entre el estrato de suelo y la base rígida es proporcionada por Poulos y Davis [[#bib0285|[28]]]. Para las simulaciones numéricas se utilizaron modelos constitutivos elásticos lineales, con las propiedades medidas en el ensayo de compresión no confinada. La geometría del modelo cuasiestático se representa en la [[#fig0020|figura 4]], donde se resalta la sección de análisis ''(T-T’)''. Debido a que las estrategias de discretización inherentes al método son diferentes, se intentó realizar una división similar del dominio.

<span id='figure_fig0020'></span>

<span id='fig0020'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_799892250-1-s2.0-S0213131515000292-gr4.jpg|center|360px|Modelo geométrico de la discretización realizada en MPM y MEF.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 4.

<span id='spar0030'></span>
Modelo geométrico de la discretización realizada en MPM y MEF.

</span>
|}

<span id='par0135'></span>
En la [[#tbl0010|tabla 2]] se establecen las características geométricas y numéricas de los modelos numéricos de cimentación superficial; la más importante de ellas, y que es preciso resaltar, es el tiempo de cómputo. En este sentido, el MEF es claramente mucho más eficiente, pues reduce 100 veces el tiempo de cálculo.

<span id='table_tbl0010'></span>

<span id='tbl0010'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabla 2.

Características generales de la simulación de una cimentación superficial

|-

! 
! MPM
! MEF
|-

| Elementos/tipo
| 1.768 cuadrados
| 1.139 triangulares (15 nodos)
|-

| Discretización
| 6.000 puntos materiales
| 9.275 puntos de Gauss
|-

| Tiempo computacional (s)
| 193
| 2
|}

<span id='par0140'></span>
En la [[#fig0025|figura 5]] se indican las tensiones a lo largo del eje de análisis. En el primer milímetro de profundidad (20% del ancho de la carga) ambos métodos numéricos presentan tensiones muy divergentes. Estos problemas en la frontera de análisis se deben a las condiciones de contorno aplicadas en el MPM y a las extrapolaciones calculadas en los puntos de Gauss para la superficie, en el caso del MEF. Además, este es un punto singular para la solución teórica, pues la sección de análisis se encuentra en el límite entre la carga y el tramo libre.

<span id='figure_fig0025'></span>

<span id='fig0025'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_799892250-1-s2.0-S0213131515000292-gr5.jpg|center|342px|Tensiones en la masa de suelo a lo largo del eje T-T’: a)tensión vertical; ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 5.

<span id='spar0035'></span>
Tensiones en la masa de suelo a lo largo del eje ''T-T’'': a) tensión vertical; b) tensión horizontal; c) tensión cortante.

</span>
|}

<span id='par0145'></span>
En el caso de la tensión vertical, el primer tramo de la curva calculada con el MEF se ajusta mejor al de la solución teórica, mientras que el MPM la subestima. Después del primer tercio, en cambio, las 3 aproximaciones arrojan prácticamente el mismo valor. Las tensiones horizontales calculadas numéricamente presentan una mayor concordancia entre ellas y conservan la misma tendencia cualitativa con la solución analítica, pero sobreestimando su valor al compararla a la solución analítica.

<span id='par0150'></span>
También se observa una cierta divergencia superficial en el caso de las tensiones cizallantes, debido a las condiciones de contorno de las soluciones. Sin embargo, un 20% por debajo de la profundidad de análisis, los resultados numéricos y teóricos coinciden.

<span id='par0155'></span>
En la [[#fig0030|figura 6]] se han grafiado los desplazamientos. Los desplazamientos verticales máximos calculados con el MPM son un 25% mayores que los calculados con el MEF. Además, el error estimado con el MPM es del 26%, mientras que el del MEF fue del 36%, comparado con la solución analítica. En el caso de los desplazamientos horizontales, cabe notar que ambas tendencias numéricas son similares y que el valor máximo de la solución del MPM es mucho más preciso que el obtenido aplicando el MEF.

<span id='figure_fig0030'></span>

<span id='fig0030'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_799892250-1-s2.0-S0213131515000292-gr6.jpg|center|331px|Desplazamientos en la masa de suelo a lo largo del eje T-T’: a)desplazamiento ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 6.

<span id='spar0040'></span>
Desplazamientos en la masa de suelo a lo largo del eje ''T-T’'': a) desplazamiento vertical; b) desplazamiento horizontal.

</span>
|}

<span id='sec0040'></span>
===4.3. Deformación de un talud artificial===

<span id='par0160'></span>
Por último, se comparan los resultados de los cálculos obtenidos utilizando el MEF y el MPM de una serie de taludes artificiales. Ambos modelos utilizaron las mismas condiciones de contorno (desplazamientos verticales y rotaciones restringidas en la base y restricción de los desplazamientos horizontales y rotaciones en el segmento vertical opuesto al talud). El material de referencia para todos los taludes correspondía a la muestra de caolín descrita anteriormente y se asumía una resistencia no drenada de ''s''<sub>''u''</sub> = 0,5 σ<sub>y</sub>. Para las simulaciones se utilizó un modelo elástico lineal perfectamente plástico, con un criterio de falla Von Mises en el modelo del MPM y Mohr-Coulomb (ϕ ≈ 0) en el modelo del MEF.

<span id='par0165'></span>
Se eligió arbitrariamente una inclinación de 45° y una altura variable de entre 1 y 8 m. De acuerdo con el ábaco de estabilidad de Taylor revisado por Steward et al. [[#bib0145|[29]]] en suelos puramente cohesivos, el factor de seguridad del talud a una altura de 8 m estaría muy cerca de la unidad (colapso). En la [[#fig0035|figura 7]] se muestra la variación del desplazamiento horizontal del borde superior del talud como producto de las deformaciones ocasionadas por el propio peso de la masa de suelo ante la variación de altura del mismo. Pueden encontrarse resultados similares en la literatura [[#bib0190|[9]]] and [[#bib0195|[10]]].

<span id='figure_fig0035'></span>

<span id='fig0035'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_799892250-1-s2.0-S0213131515000292-gr7.jpg|center|349px|Relación entre el desplazamiento horizontal en la corona del talud y la altura ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 7.

<span id='spar0045'></span>
Relación entre el desplazamiento horizontal en la corona del talud y la altura para diferentes taludes.

</span>
|}

<span id='par0170'></span>
En la [[#fig0035|figura 7]] se comparan las curvas de desplazamiento calculadas con el MPM con los resultados del programa Plaxis. Puede observarse que, aunque los esquemas de discretización utilizados son diferentes, el resultado entre ambos métodos numéricos es prácticamente el mismo. Sin embargo, destaca que, a partir de 5 m, el exceso de deformaciones impide obtener una respuesta mediante el MEF, mientras que el MPM calcula sin problemas las deformaciones hasta la altura de la falla del talud. Además, cabe señalar que, puesto que el MPM es un método explícito en el tiempo, es necesario utilizar amortiguamiento numérico para obtener una solución comparable con el MEF; de lo contrario se obtendrán valores de desplazamiento superiores en un 30%, debido al impacto generado por la carga gravitacional.

<span id='par0175'></span>
El amortiguamiento numérico es necesario ya que se pretende simular un problema cuasiestático por medio de una formulación dinámica. Para tal fin, la ecuación (3) debe amortiguarse a través del llamado ''amortiguamiento local viscoso'', buscando un equilibrio para no sobreamortiguar el modelo. Su estimación tuvo como base el primer período de oscilación de un estrato uniforme de suelo  [[#bib0265|[24]]].

<span id='par0180'></span>
La [[#fig0040|figura 8]] muestra el esquema de discretización típico adoptado por cada método, así como el estado final de deformación a escala real para la máxima altura comparable del talud (5 m de altura). En la [[#tbl0015|tabla 3]] se comparan los detalles numéricos entre ambos modelos y se observa nuevamente que, a pesar de que el MPM está preparado para resolver el problema de las grandes deformaciones, el MEF es más veloz.

<span id='figure_fig0040'></span>

<span id='fig0040'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_799892250-1-s2.0-S0213131515000292-gr8.jpg|center|357px|Discretización numérica y estado de deformación final de un talud sintético de ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 8.

<span id='spar0050'></span>
Discretización numérica y estado de deformación final de un talud sintético de 5 m de alto y 45° de inclinación por medio de métodos numéricos: a) MEF; b) MPM.

</span>
|}

<span id='table_tbl0015'></span>

<span id='tbl0015'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabla 3.

Características generales de la simulación de un talud artificial

|-

! 
! MPM
! MEF
|-

| Elementos/tipo
| 3.416 cuadrados
| 1.010 triangulares (15 nodos)
|-

| Discretización
| 8.775 puntos materiales
| 12.120 puntos de Gauss
|-

| Tiempo computacional (s)
| 134
| 30
|}

<span id='par0185'></span>
Finalmente, la [[#fig0045|figura 9]] muestra la deformación de desviación (''ɛ''<sub>1</sub> − ''ɛ''<sub>3</sub>) calculada en la masa de suelo. En este caso destaca la similitud entre la magnitud de los valores calculados y la forma general en que las deformaciones están distribuidas, sin olvidar que el fenómeno de ''kinematic locking'' puede observarse nuevamente en los resultados del MPM. En ambos gráficos también puede verse que las deformaciones de desviación se concentran en la cara del talud, lo cual da una idea del inicio de la superficie de falla.

<span id='figure_fig0045'></span>

<span id='fig0045'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_799892250-1-s2.0-S0213131515000292-gr9.jpg|center|372px|Deformación desviadora en la masa de suelo de un talud sintético utilizando ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 9.

<span id='spar0055'></span>
Deformación desviadora en la masa de suelo de un talud sintético utilizando métodos numéricos: a) MEF; b) MPM.

</span>
|}

<span id='sec0045'></span>
==5. Conclusiones==

<span id='par0190'></span>
Tanto el MPM como el software de código abierto NairnMPM se han utilizado con éxito en la simulación de problemas clásicos de la ingeniería geotécnica. Y se han obtenido con éxito soluciones cuasiestáticas como una cimentación superficial y un talud, además de resolver un problema explícito en el tiempo con solución dinámica, como la compresión no confinada.

<span id='par0195'></span>
Los resultados obtenidos han sido positivos, desde el punto de vista cualitativo y cuantitativo, y son comparables en precisión a las soluciones obtenidas con un método numérico conocido y ampliamente utilizado como es el MEF. Sin embargo, cabe señalar que, pese a que el MPM es utilizado cada vez con más frecuencia, los resultados disponibles en la literatura no siempre están acompañados de evaluaciones rigurosas de los mismos.

<span id='par0200'></span>
La limitación principal del MPM sigue siendo el tiempo de cómputo: en el caso de la solución cuasiestática, el MEF es mucho más veloz. No obstante, el MPM sigue siendo una herramienta atractiva para simular problemas de grandes deformaciones, como ha quedado demostrado y validado en el caso del ensayo de compresión no confinada y de las deformaciones en un talud.

<span id='par0205'></span>
Debido al mayor tiempo de cálculo del MPM con respecto al MEF, los dominios numéricos voluminosos deben analizarse con precisión, buscando un nivel óptimo de refinamiento de la malla entre el nivel de detalle que se requiere y el tiempo de procesamiento. Mallas muy refinadas y problemas tridimensionales pueden hacer inviable la simulación numérica utilizando los computadores convencionales.

<span id='ack001'></span>
==Agradecimientos==

<span id='par0210'></span>
Los autores desear expresar su agradecimiento a la Coordinación de Perfeccionamiento del Personal de Nivel Superior (CAPES) y al Consejo Nacional de Investigación (CNPq) del Gobierno brasileño por el apoyo económico recibido para la ejecución de este trabajo.

<span id='bibl0005'></span>
==Bibliografía==

<ol style='list-style-type: none;margin-left: 0px;'><li><span id='bib0150'></span>
[[#bib0150|[1]]] C.J. Coetzee; The Modelling of Granular Flow Using the Particle-in-Cell Method [tesis doctoral]; Universidad de Stellenbosch (2004)</li>
<li><span id='bib0155'></span>
[[#bib0155|[2]]] I.K.J. Al-Kafaji; Formulation of a Dynamic Material Point Method (MPM) for Geomechanical Problems [tesis doctoral]; Universidad de Stuttgart (2013)</li>
<li><span id='bib0160'></span>
[[#bib0160|[3]]] Z. Chen, R. Brannon; An Evaluation of the Material Point Method; Sandia National Laboratories, Albuquerque and Livermore (2002)</li>
<li><span id='bib0165'></span>
[[#bib0165|[4]]] L. Beuth, T. Benz, P.A. Vermeer, Z. Wieckowsi, C.J. Coetzee; Large deformation analysis using a quasi-static Material Point Method; Comput. Meth. Mech. (2007), pp. 1–6</li>
<li><span id='bib0170'></span>
[[#bib0170|[5]]] F.H. Harlow; The particle-in-cell computing method for fluid dynamics; Meth. Comput. Phys., 3 (1964), pp. 319–343</li>
<li><span id='bib0175'></span>
[[#bib0175|[6]]] S.G. Bardenhagen, E.M. Kober; The generalized interpolation material point method; Tech. Sci. Press, 5 (6) (2004), pp. 477–495</li>
<li><span id='bib0180'></span>
[[#bib0180|[7]]] L. Beuth, Z. Wieckowski, P.A. Vermeer; Solution of quasi-static large-strain problems by the material point method; Int. J. Numer. Anal. Meth. Geomech., 35 (2011), pp. 1451–1465</li>
<li><span id='bib0185'></span>
[[#bib0185|[8]]] F. Zabala; Modelacion de problemas geotécnicos hidromecánicos utilizando el método del punto material [tesis doctoral]; Universitat Politècnica de Catalunya (2010)</li>
<li><span id='bib0190'></span>
[[#bib0190|[9]]] L. Beuth, T. Benz, P.A. Vermeer; Large deformation analysis using a quasi-static material point method; J. Theor. Appl. Mech., 38 (1-2) (2008), pp. 45–60</li>
<li><span id='bib0195'></span>
[[#bib0195|[10]]] P.A. Vermeer, L. Beuth, T. Benz; A Quasi-Static Method for Large Deformation Problems in Geomechanics; 12th International Conference on Computer Methods and Advances in Geomechanics, Goa (2008), pp. 55–63</li>
<li><span id='bib0200'></span>
[[#bib0200|[11]]] D.Z. Zhang, X. Ma, P. Giguere; Computation of multi-material interactions using material point method; Int. Conf. Particle Based. Meth. Barcelona (2009), pp. 1–8</li>
<li><span id='bib0205'></span>
[[#bib0205|[12]]] J. Ma; Multiscale Simulation Using the Generalized Interpolation Material Point Method, Discrete Dislocations and Molecular Dynamics [tesis doctoral]; Oklahoma State University (2002)</li>
<li><span id='bib0210'></span>
[[#bib0210|[13]]] R.K. Raghav; Three Dimensional Generalized Interpolation Material Point (GIMP) Simulations in Samrai Environment [disertación de Maestría]; Oklahoma State University (2005)</li>
<li><span id='bib0215'></span>
[[#bib0215|[14]]] Z. Wieckowski; The material point method in soil mechanics problems; 21st International Congress of Theoretical and Applied Mechanics, Warsaw, Polonia (2004), pp. 3–4</li>
<li><span id='bib0220'></span>
[[#bib0220|[15]]] Z. Wieckowski; Application of dynamic relaxation in material point method; Comput. Meth. Mech (CMM-2011), Warsaw, Polonia (2011), pp. 1–2</li>
<li><span id='bib0225'></span>
[[#bib0225|[16]]] Y.J. Guo, J.A. Nairn; Three-dimensional dynamic fracture analysis using the material point method; Tech Sci. Press, 1 (1) (2006), pp. 11–25</li>
<li><span id='bib0230'></span>
[[#bib0230|[17]]] N.P. Daphalapurkar, H. Lu, D. Coker, R. Komanduri; Simulation of dynamic crack growth using the generalized interpolation material point (GIMP) method; Int. J. Fract., 143 (2007), pp. 79–102</li>
<li><span id='bib0235'></span>
[[#bib0235|[18]]] S. Zhou, J. Stormont, Z. Chen; Simulation of geomembrane response to settlement in landfills by using the material point method; Int. J. Numer. Methods Geomech., 23 (1999), pp. 1977–1994</li>
<li><span id='bib0240'></span>
[[#bib0240|[19]]] M. Numada, K. Konagai, H. Ito, J. Johansson; Material point method for run-out analysis of earthquake-induced long-traveling soil flows; JSCE J. Earthq. Eng., 27 (2003), pp. 1–4</li>
<li><span id='bib0245'></span>
[[#bib0245|[20]]] W. Shin, G.R. Miller, P. Arduino, P. Mackenzie-Helnwein; Dynamic meshing for material point method computations; Eng. Technol., 48 (9) (2010), pp. 84–92</li>
<li><span id='bib0250'></span>
[[#bib0250|[21]]] S. Andersen, L. Andersen; Modelling of landslides with the material-point method; Comput. Geosci., 14 (1) (2010), pp. 137–147</li>
<li><span id='bib0255'></span>
[[#bib0255|[22]]] M.A. Llano-Serna; Aplicação do Método do Ponto Material (MPM) a problemas geotécnicos [disertación de Maestria]; Universidad de Brasilia (2012)</li>
<li><span id='bib0260'></span>
[[#bib0260|[23]]] O. Buzzi, D.M. Pedroso, A. Giacomini; Caveats on the implementation of the generalized material point method; Tech. Sci. Press., 31 (2) (2008), pp. 85–106</li>
<li><span id='bib0265'></span>
[[#bib0265|[24]]] J. A. Nairn, Software para MPM y FEA de código abierto — NairnMPM y NairnFEA [consultado 2014]. Disponible en: http://osupdocs.forestry.oregonstate.edu/index.php/Main_Page</li>
<li><span id='bib0270'></span>
[[#bib0270|[25]]] ASTM D 2216, Standard test method for laboratory determination of water (moisture) content of soil and rock by mass; Am. Soc. Test. Mater. (2010), pp. 1–5</li>
<li><span id='bib0275'></span>
[[#bib0275|[26]]] C.M. Mast, P. Mackenzie-Helnwein, P. Arduino, G.R. Miller, W. Shin; Mitigating kinematic locking in the material point method; J. Comput. Phys., 231 (16) (2012), pp. 5351–5373</li>
<li><span id='bib0280'></span>
[[#bib0280|[27]]] K. Konagai, J. Johansson; Two dimensional lagrangian particle finite difference method for modeling large soil deformations; Struct. Eng. Eng., 18 (2) (2001), pp. 105–110</li>
<li><span id='bib0285'></span>
[[#bib0285|[28]]] H.G. Poulos, E.H. Davis; Elastic Solutions for Soil and Rock Mechanics; John Wiley, New York (1974)</li>
<li><span id='bib0145'></span>
[[#bib0145|[29]]] T. Steward, N. Sivakugan, S.K. Shukla, B.M. Das; Taylor's Slope Stability Charts revisited; Int. J. Geomech. (2011), pp. 348–352</li>
</ol>

<span id='footer-area-dummy'></span>