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.
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.
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 [1]. En la literatura pueden encontrarse aplicaciones diversas en el campo de la ingeniería geotécnica.
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 [2].
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 [3].
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.
2Antecedentes del método del punto materialSegún Beuth et al. [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 [5].
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.
Por último, se calculan las deformaciones en los PM y se actualizan las tensiones mediante una relación constitutiva. En la 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 [6].
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 [4,7]; a aspectos relacionados con la consolidación [4,8]; a asentamientos en taludes [9,10]; a cimentaciones superficiales [1,11,12]; a problemas de penetración [12,13]; a empujes activos y pasivos [1,7,10,11,14,15]; a discontinuidades, fisuras y grietas [16,17]; a asentamientos en terraplenes [18]; a la estabilidad de los taludes y las laderas [19–21]; a cimentaciones sobre suelos blandos [8], y a problemas de penetración [14,22]. Sin embargo, cabe destacar que no abundan las validaciones rigurosas sobre los resultados obtenidos.
3El método del punto materialEs 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. [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.
3.1FormulaciónLa 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:
donde σ˜˜ es el tensor de tensiones de segundo orden; x˜ representa el vector de posición de un punto, I˜˜ expresa el tensor de identidad de segundo orden; ρ simboliza el campo de densidad escalar; b˜ son las fuerzas de cuerpo y externas que actúan en un punto, y a˜ 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.
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 w˜x˜ similar al MEF, que toma la forma siguiente, después de aplicar el teorema de la divergencia:
donde t˜ representa las fuerzas de tracción en la superficie en una región A.
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 χpx˜ 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:
La ecuación (3) establece que la tasa de cambio del momento lineal en los vértices de la malla q˙˜n 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.
Las fuerzas externas en los vértices de la malla f˜next reciben las contribuciones de las tracciones en la superficie t˜ y de las fuerzas de cuerpo b˜ 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 Sn e integrándolas sobre el área de la superficie, lo cual se representa por el primer término de la ecuación (3).
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:
Por otro lado, la transferencia de tensiones se logra multiplicándolas por el gradiente de la función de forma G˜n=dSn(x˜)/dx˜, que solo deberá calcularse dentro del dominio del PM dado por la función χpx˜. De esta manera, la extrapolación se obtiene del G¯˜np, calculado de la manera siguiente:
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 [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. [23].
4Simulaciones numéricasLas simulaciones presentadas fueron ejecutadas en condiciones bidimensionales, en estado plano de tensiones y utilizando el código abierto NairnMPM 8.1.0 [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.
4.1Ensayo de compresión no confinadaUno 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 [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.
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.
El módulo de elasticidad secante calculado suponiendo un 3% de deformación fue E=405kPa; la tensión de plastificación, σy=52kPa, y la densidad, γ=18kN/m3. Y debido a que la muestra se encontraba en condiciones saturadas, se asumió ν=0,5. En la 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).
En la 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 151s en total en alcanzar la máxima deformación registrada.
En la 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[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.
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 [1] y Konagai y Johansson [27] obtuvieron unos resultados similares.
4.2Cimentación superficial sobre un estrato de suelo finitoUno 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.
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 [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 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.
En la 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.
En la 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.
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.
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.
En la 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.
4.3Deformación de un talud artificial
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 su=0,5σy. 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.
Se eligió arbitrariamente una inclinación de 45° y una altura variable de entre 1 y 8m. De acuerdo con el ábaco de estabilidad de Taylor revisado por Steward et al. [29] en suelos puramente cohesivos, el factor de seguridad del talud a una altura de 8m estaría muy cerca de la unidad (colapso). En la 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 [9,10].
En la 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 5m, 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.
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 [24].
La 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 (5m de altura). En la 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.
Finalmente, la figura 9 muestra la deformación de desviación (¿1−¿3) 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.
5Conclusiones
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.
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.
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.
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.
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.