Numerical modelling of landslides in reservoirs via the Particle Finite Element Method (PFEM)

Resumen

Se presentan resultados obtenidos mediante la aplicación del Método de Partículas y Elementos Finitos (PFEM) en la simulación de deslizamientos de ladera en embalses. Es un fenómeno complejo, por la interacción entre el material deslizado, la masa de agua del embalse, y el conjunto formado por el vaso y la presa. PFEM es un esquema numérico original con el que se ha afrontado con éxito la resolución de problemas de interacción fluido-estructura. Combina un enfoque Lagrangiano con la resolución de las ecuaciones de elementos finitos mediante la generación de una malla, que se actualiza en cada paso de tiempo. Se presentan resultados de casos de validación en los que se han reproducido ensayos en laboratorio existentes en la bibliografía técnica. Se muestran también otros cálculos más complejos, sobre la cartografía a escala real de un embalse, donde se aprecia el fenómeno de generación de la ola, su propagación por el embalse y la afección a la presa. Por último, se ha modelado el deslizamiento ocurrido en 1958 en la bahía de Lituya (Alaska), en el que la caída de 90 millones de toneladas de roca produjo una ola que alcanzó una sobreelevación máxima de 524 m en la ladera opuesta. Los resultados permiten afirmar que PFEM puede ser una herramienta útil en el análisis de riesgos frente a este tipo de fenómenos, ofreciendo una buena aproximación de las afecciones potenciales.

Abstract

The paper presents the results of the application of the Particle Finite Element Method (PFEM) to the analysis of landslides in reservoirs. This is a complex phenomenon, because of the interaction between the landslide, the still water in the reservoir and the dam. PFEM combines a Lagrangian approach with the solution of the governing equations of the problem using the FEM. A mesh connecting the initial set of particles (nodes) is re-generated in every time step. Some validation cases are presented, in which PFEM results are compared with experimental data. More complex calculations have been made over the actual geometry of reservoirs taken from the cartographic information of the sites. In these cases the wave generation, its propagation and dam overtopping are reproduced. Finally, Lituya bay rock slide in which 90 × 10⁶ tons of rocks fell on the bay, generating a huge wave that caused a maximum run-up of 524 m on the opposite shore, has been simulated in 3D. The results show that PFEM is a useful tool in risk assessment related with landslides in reservoirs as it gives a good approximation to the potential affections, thus allowing the appropriate design of protection measures.

Palabras clave

Deslizamiento ; Embalse ; Olas ; Presas ; Seguridad ; PFEM ; Riesgo

Keywords

Landslide ; Reservoir ; Wave ; Dam ; Safety ; PFEM ; Risk

1. Introducción

Suele utilizarse el término deslizamiento de ladera para referirse a un buen número de fenómenos de diferente naturaleza. Pueden caracterizarse por la velocidad con que se producen (desde centímetros por año hasta metros por segundo), por la naturaleza de los materiales involucrados (rocas, material granular, coladas de barro, etc.), y por las causas que lo provocan (meteorización, erosión interna, variación del nivel freático, etc.). Varnes [1] define los movimientos de ladera como «un movimiento hacia abajo y hacia afuera de los materiales que forman una ladera bajo la influencia de la gravedad». Ayala [2] completa la anterior definición añadiendo que «se acompañan a veces por otras fuerzas naturales como las sísmicas, volcánicas, o la presión de gases, y representa la materia sólida en porcentaje del peso más del 70%».

Los deslizamientos de ladera son fenómenos que a lo largo de la historia han provocado la rotura de varias presas y la pérdida de un gran número de vidas humanas. Aunque la probabilidad de ocurrencia es muy baja, su importancia en lo referente a pérdidas humanas y materiales es importante, constituyendo el tercer riesgo natural por víctimas, tras terremotos e inundaciones [2] .

Los deslizamientos en embalses tienen un interés especial, porque las variaciones de nivel de agua (especialmente los desembalses rápidos) pueden favorecer su ocurrencia. La presencia del embalse durante un periodo de tiempo suficiente hace que el material de las márgenes se sature, con lo que aumenta la presión intersticial y por tanto se reduce la tensión efectiva. Este efecto desestabilizador se ve parcialmente compensado por el incremento de la tensión total producido por la presión hidrostática ejercida por el agua del embalse [3] . Un desembalse rápido elimina este efecto estabilizador en un lapso de tiempo que en general no permite que se disipe la presión intersticial (esto depende de la permeabilidad del material y de la velocidad de descenso del nivel de embalse, pero es muy frecuente), con lo que aumenta la probabilidad de ocurrencia de un deslizamiento.

En este campo, existe un importante esfuerzo investigador en todo el mundo, principalmente concentrado en el desarrollo de tecnología para la predicción de los deslizamientos. Con mayor frecuencia se afronta el problema analizando las características de los materiales que forman la ladera, la presencia de agua, el grado de fracturación, etc. En ocasiones, se han utilizado herramientas matemáticas tan singulares como las técnicas de análisis fractal [4] para analizar la topografía de las zonas susceptibles, así como la distribución geográfica de los deslizamientos ocurridos previamente.

El presente trabajo se centra en los deslizamientos rápidos, que son los más peligrosos, por la ausencia de tiempo de aviso y evacuación de la población potencialmente afectada. Concretamente, se ha estudiado el caso de los deslizamientos producidos en laderas de embalses, analizando la formación de la ola, su propagación, y la afección a la presa y a las márgenes. No se analizan ni las causas que los provocan, ni se trata de predecirlos ni de estimar su probabilidad de ocurrencia. Se estudian las posibles consecuencias de un deslizamiento con unas características básicas (geometría, material) definidas.

La herramienta utilizada para tal fin es un código de cálculo desarrollado en CIMNE durante los últimos años basado en la técnica de partículas y elementos finitos (PFEM), que constituye un método innovador para el análisis de problemas de fluidos con grandes movimientos de la superficie libre, así como para estudiar la interacción entre fluidos y estructuras. Puede encontrarse información detallada sobre los fundamentos del método en [5] , [6]  and [7] . El código de cálculo se ha validado para su aplicación a diferentes problemas de ingeniería hidráulica [8] . Con esa base, se ha aplicado al análisis de varios casos de deslizamientos de ladera documentados en la bibliografía técnica, tanto de ensayos de laboratorio como de casos reales.

2. Estado del arte

2.1. Clasificación de los deslizamientos

Existen diversas clasificaciones de los deslizamientos de ladera, atendiendo a diferentes aspectos, como la causa que lo provoca, la velocidad, o la naturaleza del material deslizado. Introducciones más detalladas al fenómeno pueden encontrarse en Ayala [2] , y en Haddad [9] . Desde el punto de vista de las consecuencias, los más importantes son sin duda los deslizamientos rápidos. La tabla 1 muestra la clasificación elaborada por el UNESCO Working Party del World Landslide Inventory de 1995, basada en la clasificación de Varnes de 1978.

Tabla 1. Severidad de los deslizamientos en función de la velocidad. Fuente: Working Party of the World Landslide Inventory, tomado de Ayala [2] .

Velocidad	Designación	Consecuencias
5 m/s	Extremadamente rápida	Catástrofe de violencia mayor; edificaciones expuestas totalmente destruidas, población muerta por el impacto o el material desplazado o por disgregación de la masa deslizada
3 m/min	Muy rápida	Algunas vidas perdidas porque la velocidad es demasiado grande para permitir el escape; destrucción mayor
1,8 m/h	Rápida	Escape y evacuación posible; estructura, posesiones y equipo destruidos
13 m/mes	Moderada	Estructuras invulnerables pueden ser mantenidas si están localizadas a corta distancia del pie del movimiento; daño extenso en las situadas sobre el movimiento
1,6 m/año	Lenta	Carreteras y estructuras vulnerables pueden ser mantenidas con trabajos fuertes e intensos si el movimiento no dura demasiado y si los movimientos diferenciales en los márgenes del movimientos se distribuyen en una zona amplia
0,16 m/año	Muy lenta	Algunas estructuras permanentes no sufren daño, o si lo sufren, pueden ser reparadas
> 0,16 m/año	Extremadamente lenta	Sin daños a estructuras construidas con cuidado

Otras clasificaciones se basan en la naturaleza del material movilizado, o bien en las causas desencadenantes. El trabajo presentado no se centra en estos aspectos de los movimientos de ladera, por lo que no se presentan estas clasificaciones. Sin embargo, existe una clara relación entre lo anterior y la cinemática del deslizamiento, característica clave en la interacción con el agua del embalse, y por tanto en la formación de la ola y en sus consecuencias. Por esta razón se describen a continuación las características principales de los diferentes tipos de deslizamientos rápidos.

2.1.1. Flujos de derrubios y coladas de barro

Estos tipos de deslizamientos se caracterizan por la presencia importante de agua. Se diferencian principalmente por el diámetro de las partículas sólidas de la mezcla, que en las coladas de barro presentan al menos un 50% de materiales del tamaño de arenas, limos y arcillas, mientras que en los flujos de derrubios hay una proporción mayor de material grueso [1] . La principal característica de ambos, desde el punto de vista de su modelación numérica, es que la velocidad relativa entre la fase sólida y la fase líquida es pequeña, con lo que es admisible tratar la mezcla como un continuo [10] . Además, la presencia de agua hace que su comportamiento pueda asemejarse al de un fluido denso.

2.1.2. Avalanchas granulares

Se producen a partir de la pérdida de estabilidad de grandes bloques de roca, que, a medida que descienden y aumenta su velocidad, se van disgregando hasta convertirse en una masa de material granular (derrubios) con partículas de tamaño variable. De este modo el proceso pasa de tener características de sólido rígido a poder asemejarse a un fluido denso [10] . Otra posible causa es la movilización de depósitos de rocas debido a un terremoto. La gran diferencia con las coladas de barro o los flujos de derrubios es la ausencia de agua.

2.1.3. Flujos deslizantes

Este término se refiere a los movimientos de materiales procedentes de escombreras, es decir, de rellenos antrópicos. La naturaleza de estos materiales puede variar, pero en general tienen la característica común de estar poco compactados. Diferentes causas pueden provocar la licuefacción del material, con la consecuente pérdida de estabilidad que genera el deslizamiento. Por las características mencionadas, este tipo de movimientos se han simulado mediante técnicas de mecánica de fluidos [9] .

2.2. Deslizamientos en embalses

Hasta este punto se ha prestado atención únicamente al deslizamiento de ladera. El trabajo se centra sin embargo en los deslizamientos ocurridos en laderas de embalses, o más propiamente, en aquellos casos en que la masa deslizada interacciona con una masa de agua inicialmente en reposo (todos los desarrollos pueden aplicarse a deslizamientos ocurridos en lagos, o junto al mar, como se verá en alguno de los ejemplos). Si el fenómeno de generación y propagación del deslizamiento es complejo, más aún lo es la interacción con una masa de agua en reposo, y la formación y propagación de la consiguiente ola.

Este es quizá el punto más difícil de afrontar desde el punto de vista de la modelación. Hay numerosos trabajos relacionados con el tema, muchos de ellos basados en ensayos de laboratorio. El mecanismo de generación de la ola utilizado varía de un experimento a otro. En algunos casos se han utilizado sólidos rígidos, como una cuña que desliza, una placa que gira [11] , un bloque que cae verticalmente, un sólido en forma de medio elipsoide [12] , o bloques paralelepípedos [13] . En otros, la masa deslizada se ha simulado como un material granular, para lo cual se han desarrollado equipos de laboratorio capaces de mantener la forma de la masa, y controlar su velocidad [14] . Existen también aproximaciones analíticas, como las de Noda [15] o Di Risio y Sammarco [16] . Por último, el problema se ha afrontado también previamente mediante diferentes técnicas de cálculo numérico, como Monaghan y Kos [17] con técnicas SPH, o Pastor et al. [18] , con elementos finitos y SPH.

Fritz et al. [19] realizaron un estudio paramétrico mediante investigación experimental, consistente en una batería de casos de deslizamientos granulares, en los que estudiaron la relación entre las variables propias del deslizamiento y el tipo de ola generada. A la vista de los resultados de esta campaña propusieron caracterizar los deslizamientos en base a tres variables: a) espesor, b) volumen y c) número de Froude. En su opinión este último parámetro es el que gobierna el fenómeno, y se define como el cociente entre la velocidad del centro de masas del deslizamiento en el momento del impacto y la raíz cuadrada del producto de la aceleración de la gravedad por la profundidad inicial del agua.

La fase de propagación de la ola en el embalse puede ser estudiada con diferentes herramientas, siendo las más comunes los modelos numéricos integrados en profundidad.

En este campo también se han realizado campañas de ensayos en modelo físico, destacando la de Huber y Hager [20] , que extrajeron una expresión para la amplitud máxima de la ola en función de la distancia al punto de deslizamiento y del ángulo respecto de la dirección del movimiento de la masa.

El otro aspecto de interés del fenómeno es la caracterización de la afección de la ola tanto a la presa como a las márgenes del embalse. Este campo, como los anteriores, se ha estudiado mediante modelación física, habiéndose obtenido expresiones que estiman el valor de la altura de ascenso de la ola sobre el paramento y el volumen de vertido en función de las características de la ola. En Vischer y Hager [21] pueden encontrarse estas expresiones, además de otras relativas a otros aspectos del fenómeno.

Garrote y Laguna [22] realizaron una revisión del estado del arte en el campo de los deslizamientos de ladera en embalses.

2.3. Modelación numérica de deslizamientos de ladera

En los últimos años se ha invertido un considerable esfuerzo en todo el mundo en el desarrollo de métodos robustos y eficientes para el análisis de problemas complejos de dinámica de fluidos con grandes movimientos de la superficie libre del fluido y la posible existencia de cuerpos sumergidos total o parcialmente. Los ejemplos de este tipo son comunes en problemas de hidrodinámica de barcos, estructuras portuarias off-shore, flujos en lámina libre en canales de descarga de aliviaderos de presas, tanques de agua y de gas licuado, reactores de mezclas, procesos de llenado de moldes, etc.

Las dificultades típicas del análisis numérico de los problemas de dinámica de fluidos con superficie libre muy inestable e irregular son fundamentalmente el tratamiento de los términos convectivos y de la condición de incompresibilidad en las ecuaciones de fluido, el modelado y seguimiento de la superficie libre, la transferencia de la información entre los dominios del fluido y el contorno vía las interfases de contacto, el modelado de la rotura de olas, la actualización eficiente de las mallas de elementos finitos para la estructura y el fluido, etc.

La mayoría de estos problemas desaparecen si se utiliza una descripción Lagrangiana para formular las ecuaciones de gobierno de los dominios del fluido y del contorno. En la formulación Lagrangiana se sigue el movimiento de cada una de las partículas de líquido o del sólido de forma individual y, consecuentemente, los nodos en una malla de elementos finitos pueden considerarse como «partículas» en movimiento. Por consiguiente, el movimiento de la malla que discretiza el dominio total (incluyendo los dominios del fluido y de la estructura) se sigue durante la solución en el tiempo.

Una formulación Lagrangiana muy utilizada en la comunidad científica es la denominada Smooth Particle Hydrodynamics (SPH) [40] . Este método se ha desarrollado mucho en los últimos años, por su relativa sencillez de programación. Las dificultades principales de este método son la aplicación de las condiciones de contorno de una manera eficiente, y el cumplimiento de la condición de incompresibilidad.

En los últimos años, CIMNE ha desarrollado un tipo particular de formulación Lagrangiana para resolver problemas en los que interviene la interacción entre fluidos y sólidos. El método se denomina Método de Partículas y Elementos Finitos (PFEM). El PFEM trata los nodos en la malla, tanto en los dominios del fluido como de la estructura, como partículas que pueden moverse libremente e incluso separarse del dominio principal del fluido representando, por ejemplo, el efecto de gotas o chorreones de agua. Una malla de elementos finitos conecta los nodos que definen el dominio discretizado donde se resuelven las ecuaciones de gobierno de la mecánica de fluidos (para el líquido) y de la mecánica de sólidos (para la estructura) en la forma estándar del MEF [6] , [7] , [8] , [23] , [24] , [25] , [26]  and [27] .

Una ventaja de la formulación Lagrangiana es que los términos convectivos desaparecen de las ecuaciones del fluido. La dificultad, sin embargo, se transfiere al problema de mover adecuadamente (y eficientemente) los nodos de la malla. En general, suele ser necesario remallar a lo largo de la solución en cada paso de tiempo. La técnica desarrollada por CIMNE utiliza un procedimiento de regeneración de la malla que mezcla elementos de diferentes formas mediante un método extendido de Delaunay [28] . Estos elementos finitos poliédricos necesitan funciones de forma especiales. En los desarrollos realizados por CIMNE se utilizan las funciones de forma del denominado método de elementos finitos sin malla (MFEM) [23] .

En la formulación Lagrangiana existe todavía la necesidad de tratar adecuadamente la condición de la incompresibilidad en el fluido. El uso de interpolaciones de elementos finitos estándar puede conducir al bloqueo de la solución por deformación volumétrica, a menos que se tomen algunas precauciones. En la literatura se encuentran diferentes procedimientos de elementos finitos para aliviar el problema del bloqueo en fluidos incompresibles [29] . Un objetivo general de todos estos métodos es poder utilizar elementos de bajo orden con el mismo tipo de interpolaciones para las variables de velocidad y presión. En el método se utiliza una técnica de estabilización basada en el método de cálculo finito (FIC). En las referencias [30] , [31]  and [32] se pueden encontrar diferentes aplicaciones de la técnica FIC para problemas de fluidos incompresibles utilizando elementos triangulares y tetraédricos lineales.

El PFEM tiene muchas ventajas para seguir el movimiento de las partículas del fluido en flujos en donde existen grandes desplazamientos de la superficie libre como en el caso de olas que rompen sobre una estructura. Este es el caso también de los deslizamientos de ladera en embalses, razón por la cual se ha desarrollado el presente trabajo. Recordamos que la información en el método PFEM es típicamente nodal, es decir, la malla de elementos finitos se utiliza fundamentalmente para obtener los valores de las variables de estado (por ejemplo las velocidades, presiones, etc.) en los nodos.

3. Ideas básicas del Método de Partículas y Elementos Finitos (PFEM)

Consideremos un dominio que contiene subdominios de fluido y de sólido. Las partículas de fluido en movimiento interaccionan con los contornos del sólido induciendo por tanto la deformación del sólido que a su vez afecta al movimiento del fluido y, por consiguiente, el problema está totalmente acoplado.

En la técnica PFEM, tanto los dominios del fluido como del sólido se modelan utilizando una formulación Lagrangiana actualizada (fig. 1 ). Se utiliza el método de los elementos finitos para resolver las ecuaciones de gobierno en ambos dominios. Por tanto debe generarse una malla que discretiza estos dominios para resolver las ecuaciones de gobierno para el fluido y el sólido con el método de elementos finitos tradicional. Se destaca de nuevo que los nodos que discretizan los dominios del fluido y del sólido pueden interpretarse como partículas cuyo movimiento se sigue durante la solución transitoria.

Figura 1. Descripción Lagrangiana actualizada para un contorno que contiene un dominio de fluido y otro de sólido.

La calidad de la solución numérica depende lógicamente de la discretización utilizada, como sucede en el MEF estándar. Con el PFEM es posible definir diferentes tamaños de malla en diversas zonas del dominio de análisis, para mejorar la solución en zonas donde ocurran grandes movimientos del fluido o de la estructura, o bien donde se quiera obtener mayor precisión.

Como se ha explicado, la formulación Lagrangiana permite seguir el movimiento de cada partícula individual del fluido (un nodo). Esto es útil para modelar la separación de las partículas del líquido del dominio principal del fluido y para seguir su movimiento como partículas individuales con una velocidad inicial y sometidas a la fuerza de la gravedad.

En resumen, una solución típica con el PFEM involucra las etapas siguientes (fig. 2 ):

• Discretizar los dominios del fluido y de la estructura con una malla de elementos finitos. En nuestro trabajo utilizamos un esquema de generación de malla basado en una teselación extendida de Delaunay [23] , [24]  and [25] .
• Identificar los contornos externos para los dominios del fluido y del sólido. Esta es una etapa esencial ya que algunos contornos (como la superficie libre en los fluidos) pueden haberse distorsionado severamente durante el proceso de solución, incluyendo la separación y reentrada de nodos. En nuestro trabajo utilizamos el método Alpha Shape [33] para la definición de los nodos del contorno.
• Resolver las ecuaciones acopladas del movimiento Lagrangiano para los dominios del fluido y de la estructura. Calcular las variables de estado relevante en ambos dominios en cada paso de tiempo: las velocidades, la presión y las tensiones viscosas en el fluido, y los desplazamientos, las deformaciones y las tensiones en la estructura.
• Mover los nodos de la malla a una nueva posición en función del paso de tiempo. Esta etapa es típicamente una consecuencia del proceso de solución de la etapa 3.
• Generar una nueva malla en caso necesario.
• Volver a la etapa 2 y repetir el proceso de solución para el paso de tiempo siguiente.

Figura 2. Secuencia de etapas para actualizar una «nube» de nodos que representa un continuo que contiene un fluido y un sólido del tiempo n (t  = tn ) al n  + 2 (t  = tn  + 2Δt ).

Para completar esta descripción, en la figura 3 se muestra un ejemplo típico de una aplicación bi-dimensional del PFEM [6] . Las figuras corresponden al análisis del problema de la rotura de una columna de agua en el interior de un recipiente. La figura 3 a muestra la malla inicial que discretiza el dominio del fluido y las paredes del recipiente. Se aprecia en las figuras 3 b y 3 c la malla en los dominios del fluido y del recipiente en dos instantes de tiempo diferentes.

Figura 3. Rotura de una columna de agua. (a) Discretización del dominio del fluido y de las paredes del recipiente. (b) y (c) Mallas en el fluido y en las paredes del recipiente en dos tiempos distintos.

En la página web www.cimne.com/pfem pueden encontrarse publicaciones sobre PFEM y videos de diferentes aplicaciones.

4. Formulación del PFEM para un medio continuo Lagrangiano

Consideremos un medio continuo compuesto de materiales con distintas propiedades. Estos materiales pueden ser fluidos incompresibles (agua), medios granulares o rocosos y estructuras con materiales tradicionales.

Las ecuaciones a resolver sobre el dominio «compuesto» son las de un continuo Lagrangiano (excluyendo efectos térmicos), es decir:

4.1. Conservación de cantidad de movimiento

$\rho {\frac {\partial v_{i}}{\partial t}}={\frac {\partial {\sigma }_{ij}}{\partial x_{j}}}+$$b_{i}\quad en\quad \Omega$

( 1)

4.2. Relación entre presión y velocidad

${\frac {1}{K}}{\frac {\partial p}{\partial t}}-{\frac {\partial v_{i}}{\partial x_{i}}}=$$0\quad en\quad \Omega$

( 2)

En esas ecuaciones ${\textstyle v_{i}}$ es la velocidad a lo largo del eje cartesiano global i , p es la presión, ρ y K son la densidad y el módulo de compresibilidad del material, respectivamente, b_i las fuerzas de masa y σ_ij las tensiones (de Cauchy) que están relacionadas con las velocidades por la ecuación constitutiva[37] :

4.2.1. Medio incompresible

${}^{t+1}{\sigma }_{ij}=2\mu {\dot {\varepsilon }}_{ij}-$${}^{t+1}ps_{ij}$

( 3a)

4.2.2. Medio compresible o cuasi-incompresible

${}^{t+1}s_{ij}={}^{t}{\overset {\mbox{ˆ}}{\sigma }}_{ij}+$$2\mu {\dot {\varepsilon }}_{ij}+\lambda {\dot {\varepsilon }}_{ii}{\delta }_{ij}$

( 3b)

siendo ${\textstyle {\overset {\mbox{ˆ}}{\sigma }}_{ij}}$ son las componentes del tensor

$[{\overset {\mbox{ˆ}}{\sigma }}]={\frac {1}{J}}{\boldsymbol {\mbox{F}}}^{T}{\boldsymbol {\mbox{S}}}{\boldsymbol {\mbox{F}}}$

( 3c)

donde S es el segundo tensor de tensiones de Piola-Kirchhoff, F el tensor gradiente de deformación y J  = det F .

En las ecuaciones anteriores ${\textstyle {\dot {\varepsilon }}_{ij}}$ es la tasa de deformación ${\textstyle \left({\dot {\varepsilon }}_{ij}=\right.}$$\left.{\frac {1}{2}}\left({\frac {\partial v_{i}}{\partial x_{j}}}+\right.\right.$$\left.\left.{\frac {\partial v_{j}}{\partial x_{i}}}\right)\right)$ y ^t (·) indica valores en el tiempo t . Los parámetros μ y λ toman los siguientes valores para un fluido o un sólido.

• μ   : viscosidad; ${\textstyle \lambda =\Delta tK-{\frac {2\mu }{3}}}$
• ${\textstyle \mu ={\frac {\Delta tG}{J}}}$ ; ${\textstyle \lambda ={\frac {2G\nu \Delta t}{J(1-2\nu )}}}$ , siendo G el módulo de rigidez, ν el coeficiente de Poison y Δt el incremento de tiempo.

Las ecuaciones (1) –((3a) , (3b)  and (3c) ) se completan con las condiciones de contorno estándar en velocidades y fuerzas de superficie en el contorno [37] .

4.3. Discretización de las ecuaciones

La forma variacional de las ecs. (1) –((3a) , (3b)  and (3c) ) se obtiene por el método de residuos ponderados (o por la expresión de trabajos virtuales equivalente). Un problema esencial en la solución de estas ecuaciones es la satisfacción de la ec. (2) para material incompresible, como el agua, o algunos suelos. En estos casos K  =∞ y la ec. (2) se transforma en ${\textstyle {\frac {\partial v_{i}}{\partial x_{i}}}=}$$0$ .

En la literatura existen diversos métodos numéricos para tratar el problema de la incompresibilidad. En nuestro caso se utiliza una formulación estabilizada basada en la técnica de cálculo finito (FIC) [30] , [31]  and [32] . La esencia de este método es la solución de una ecuación presión-velocidad modificada que se escribe como

${\frac {1}{K}}{\frac {\partial p}{\partial t}}-{\frac {\partial v_{i}}{\partial x_{i}}}-$$\sum _{i=1}^{3}\tau {\frac {\partial q}{\partial x_{i}}}\left[{\frac {\partial p}{\partial x_{i}}}+\right.$$\left.{\pi }_{i}\right]=0$

( 4)

donde q son funciones de peso y τ es un parámetro de estabilización dado por [30]

$\tau ={\left({\frac {2\rho \vert {\boldsymbol {\mbox{v}}}\vert }{h}}+{\frac {8\mu }{3h^{2}}}\right)}^{-1}$

( 5)

En la expresión anterior, h es una longitud característica de cada elemento finito y |v | es el módulo del vector de velocidad. En la ec. (4)π_i son variables auxiliares que se interpretan como el gradiente de presión proyectado sobre la malla de elementos finitos. Estas variables garantizan que el término de estabilización sea pequeño. De hecho ese término puede interpretarse como una suma ponderada de las ecuaciones de conservación de la cantidad de movimiento [30] . El conjunto de ecuaciones de gobierno se completa añadiendo una ecuación adicional que establece que los términos de estabilización se anulan en el equilibrio [7] , [30]  and [32]

${\int }_{V}\tau w_{i}\left({\frac {\partial p}{\partial x_{i}}}+\right.$$\left.{\pi }_{i}\right)dV=0\quad i=1,n_{d}\quad \left(no\quad suma\quad en\right){\mbox{ }}i$

( 6)

donde ${\textstyle w_{i}}$ son funciones de peso arbitrarias.

El proceso de discretización implica interpolar las incógnitas del problema (las velocidades, ${\textstyle v_{i}}$ , la presión p y los gradientes de presión proyectados π_i ). En nuestro trabajo se utiliza una interpolación lineal de todas las variables sobre triángulos de 3 nodos (en 2D) y tetraedros de 4 nodos (en 3D). El conjunto de ecuaciones discretizadas utilizando el método de Galerkin tiene la forma siguiente

4.4. Conservación de cantidad de movimiento

${\boldsymbol {\mbox{M}}}{\dot {\overline {\boldsymbol {\mbox{v}}}}}+$${\boldsymbol {\mbox{K}}}{\overline {\boldsymbol {\mbox{v}}}}-$${\boldsymbol {\mbox{G}}}{\overline {\boldsymbol {\mbox{p}}}}=$${\boldsymbol {\mbox{f}}}$

( 7)

4.5. Presión-velocidad

${\overline {\boldsymbol {\mbox{M}}}}{\dot {\overline {\boldsymbol {\mbox{p}}}}}-$${\boldsymbol {\mbox{G}}}{\overline {\boldsymbol {\mbox{v}}}}-$${\boldsymbol {\mbox{L}}}{\overline {\boldsymbol {\mbox{p}}}}-$${\boldsymbol {\mbox{Q}}}{\overline {\boldsymbol {\pi }}}=$${\boldsymbol {0}}$

( 8)

4.6. Proyección del gradiente de presión

${\overset {\mbox{ˆ}}{\boldsymbol {\mbox{M}}}}{\overline {\boldsymbol {\pi }}}+$${\boldsymbol {\mbox{Q}}}^{T}{\overline {\boldsymbol {\mbox{p}}}}=$${\boldsymbol {0}}$

( 9)

En las ecs. (7) –(9)${\textstyle ({\overline {\cdot }})}$ denota variables nodales, ej. ${\textstyle ({\dot {\overline {\cdot }}})={\frac {\partial }{\partial t}}({\overline {\cdot }})}$ . Las matrices y vectores de estas ecuaciones pueden verse en [7]  and [37] .

En el cuadro 1 se muestra el algoritmo de integración temporal utilizado para resolver las ecs. (7) –(9) . Para más detalle consultar [7] .

Cuadro 1.

Algoritmo básico PFEM para un dominio fluido

1. BUCLE SOBRE PASOS DEL TIEMPO, ${\textstyle t=1}$ , NTIME

Valores conocidos

${\textstyle {}^{t}{\overline {\boldsymbol {\mbox{x}}}},{}^{t}{\overline {\boldsymbol {\mbox{v}}}},{}^{t}{\overline {\boldsymbol {\mbox{p}}}},{}^{t}{\overline {\boldsymbol {\pi }}},{}^{t}{\overline {T}},{}^{t}\mu ,{}^{t}{\boldsymbol {\mbox{f}}},{}^{t}{\boldsymbol {\mbox{q}}},{}^{t}C,{}^{t}V,{}^{t}M}$

2. BUCLE SOBRE NÚMERO DE ITERACIONES, ${\textstyle i=1}$ , NITER

${\textstyle \bullet }$ Cálculo de velocidades (ec. [[[#eq0045|7]] ])

${\textstyle \left[{\frac {1}{\Delta t}}M+K\right]{}^{t+1}{\overline {v}}^{i+1}=}$${}^{t+1}f+G{}^{t+1}{\overline {p}}^{i}+{\frac {1}{\Delta t}}M{}^{t}{\overline {v}}$

${\textstyle \bullet }$ Cálculo de presiones (ec. [(8 )]

${\textstyle \left[{\boldsymbol {\mbox{L}}}+\right.}$$\left.{\frac {1}{\Delta t}}{\overline {\boldsymbol {\mbox{M}}}}\right]{}^{t+1}{\overline {p}}^{i+1}=$$-G{}^{T}{}^{t+1}{\overline {v}}^{i+1}-Q{}^{t+1}{\overline {\pi }}^{i}+$${\frac {1}{\Delta t}}{\overline {\boldsymbol {\mbox{M}}}}{}^{t}{\overline {\boldsymbol {\mbox{p}}}}$

${\textstyle \bullet }$ Cálculo de gradientes de presión proyectados (ec. [[[#eq0055|9]] ])

${\textstyle {}^{n+1}{\overline {\pi }}^{i+1}=}$$-{\overset {\mbox{ˆ}}{M}}_{D}^{-1}\left[Q^{T}\right]{}^{t+1}{\overline {p}}^{i+1}{\begin{array}{ccc}&,&{\overset {\mbox{ˆ}}{M}}_{D}=diag\left[{\overset {\mbox{ˆ}}{M}}_{D}\right]\end{array}}$

${\textstyle \bullet }$ Actualizar posición de nodos:

${\textstyle {}^{t+\Delta t}{\overline {x}}^{i+1}=}$${}^{t}x^{i}+{}^{t+\Delta t}v^{i+1}\Delta t$

Nueva «nube» de nodos ${\textstyle {}^{t+\Delta t}C^{i+1}}$

Comprueba convergencia ${\textstyle \rightarrow }$ NO ${\textstyle \rightarrow }$ Nueva iteración ${\textstyle i\rightarrow i+1}$

${\textstyle \downarrow }$ SÍ

Nuevo paso de tiempo ${\textstyle t\rightarrow t+1}$

${\textstyle \bullet }$ Identifica nuevos contornos: ${\textstyle {}^{t+1}V}$

${\textstyle \bullet }$ Genera malla:${\textstyle {}^{t+1}M}$

Regresa a 1

5. Particularidades del PFEM

5.1. Generación de la malla

Uno de los puntos clave para el éxito del PFEM es la rápida regeneración de una malla en cada paso de tiempo a partir de la posición de los nodos en el dominio espacial. En nuestro trabajo hemos generado la malla utilizando la denominada Teselación Extendida de Delaunay (EDT) [24] . El método EDT permite generar mallas no estándar combinando elementos de formas poliédricas arbitrarias (triángulos, quadriláteros y otros polígonos en 2D y tetraedros, hexaedros y poliedros arbitrarios en 3D) en un tiempo de cálculo n , donde n es el número total de nodos en la malla ( figura 4 ). Las funciones de forma de clase C^∘ de los elementos pueden obtenerse sencillamente utilizando la denominada formulación del FEM sin malla [28] .

Figura 4. Generación de mallas no estándar combinando diferentes polígonos (en 2D) y poliedros (en 3D) utilizando la técnica de Delaunay extendida [24] .

Una vez que se ha generado la malla en cada paso de tiempo se aborda la solución numérica utilizando el algoritmo descrito en el apartado previo.

5.2. Identificación de los continuos

Una de las principales etapas en el PFEM es la definición correcta de los contornos. A veces los nodos de contorno se definen explícitamente como tales de manera diferente de los contornos en el interior de dominio. En otros casos, solo se dispone de la información del conjunto total de nodos y el algoritmo debe reconocer los nodos de contorno.

La utilización de la partición de Delaunay extendida facilita reconocer los nodos de contorno.

Considerando que los nodos siguen una distribución variable h (x ), siendo h (x ) la distancia mínima entre dos nodos, se ha utilizado el siguiente criterio para identificar nodos de contorno. Para todos los nodos se busca el círculo (la esfera en 3D) de radio αh que los contiene. Todo par de nodos contenido en un círculo que no contenga otro nodo en su interior se consideran nodos de contorno. En la práctica, α es un parámetro cercano a, pero mayor que, uno. Este criterio es muy similar al método Alpha Shape [33] .

Una vez que se ha decidido qué nodos estén sobre los contornos, debe definirse la superficie de contorno. Dicha superficie se define por todas las superficies poliédricas (o polígonos en 2D) que tienen todos sus nodos sobre el contorno y que pertenecen a un único poliedro.

La figura 5 muestra ejemplos de la técnica de reconocimiento de contornos utilizando el método Alpha Shape.

Figura 5. Ejemplos de reconocimiento del contorno con el método de Alpha Shape. Los círculos de radio mayor que αh (x ) y que no contienen nodos en su interior definen los nodos de contorno.

La definición correcta del contorno es muy importante para definir la normal externa a la superficie de contorno. Además, en formas débiles (Galerkin) como las aquí utilizadas es muy importante una evaluación correcta del volumen del dominio. En el criterio empleado, el error en la definición del contorno es proporcional a h lo que es un error aceptable. La única de forma de obtener una superficie de contorno más precisa es reduciendo la distancia entre los nodos.

El método descrito permite también identificar las partículas de fluido aisladas fuera del dominio principal del fluido. Estas partículas se tratan como parte del contorno exterior donde la presión se fija al valor atmosférico.

La figura 6 muestra un ejemplo esquemático del proceso para identificar las partículas individuales (o un grupo de partículas) comenzando de una colección dada de nodos. Una aplicación práctica del método para identificar las partículas de la superficie libre se muestra en la figura 7 . El ejemplo corresponde al análisis del movimiento de un fluido dentro de un recipiente de forma elipsoide oscilante. Adviértase que el método captura las gotas de agua individuales que se separan de la superficie libre durante la oscilación del recipiente.

Figura 6. Identificación de partículas individuales (o un grupo de partículas) partiendo de una colección de nodos.

Figura 7. Movimiento de un líquido sobre un recipiente oscilante. (a) Distribución original de partículas (nodos) antes de la oscilación. (b) Posición de las partículas del líquido en dos tiempos diferentes. Las partículas de contorno que representan la superficie libre, las gotas de líquido y las paredes del recipiente se representan con un color más pálido. Las flechas denotan los vectores de velocidad para cada partícula.

5.3. Tratamiento del contacto entre la interfase del fluido y los sólidos

La condición de velocidades o presiones prescritas en los contornos del sólido puede aplicarse en el PFEM en forma fuerte a los nodos del contorno. Estos nodos pueden pertenecer a contornos exteriores fijos, o a contornos móviles relacionados o vinculados a los sólidos que intereactúan con el fluido. En algunos problemas es útil definir una capa de nodos adyacente al contorno exterior en el fluido donse se impone la condición de velocidad prescrita. Estos nodos típicamente permanecen fijos durante el proceso de solución. El contacto entre las partículas de agua y los contornos de sólido se tiene en cuenta a través de la condición de incompresibilidad, que previene de forma natural que los nodos de agua penetren en los contornos de sólido . Esta forma sencilla de tratar el contacto entre el agua y las paredes del sólido es otra característica atractiva de la formulación PFEM [7] , [30]  and [40] .

5.4. Modelización de un sólido deformable como un fluido viscoso

Una forma simple y efectiva para analizar un movimiento de sólidos rígidos (o cuasi-rígidos) en interior de fluidos con el PFEM, es modelar el sólido como un fluido con una viscosidad mucho mayor que la del fluido que le rodea. Puede aplicarse entonces el esquema del cuadro 1 y resolver al mismo tiempo para el movimiento simultáneo de ambos dominios (el fluido real y el fluido ficticio que modela el cuerpo cuasi-rígido).

Este es ciertamente un procedimiento muy interesante para modelar un deslizamiento considerándolo como un cuasi-fluido. Esta técnica, ya utilizada por otros autores [18]  and [36] , ha sido la escogida en este trabajo para modelar el material del deslizamiento.

Este procedimiento puede extenderse para incluir la deformación elástica del material deslizante considerado como un fluido visco-elástico.

6. Aplicación de PFEM al análisis de deslizamientos de ladera en embalses

6.1. Descripción de los estudios realizados

Por las características explicadas del PFEM, parecía razonable utilizarlo para estudiar el fenómeno de generación de olas por deslizamientos, su propagación y la afección producida sobre las márgenes del embalse y la presa. El método ya ha sido aplicado con éxito para analizar el efecto del oleaje en obras portuarias [38]  and [39] . Se trata de un fenómeno en el que hay un movimiento de una gran masa, que impacta contra un fluido inicialmente en reposo, con lo que se genera una ola que puede tener grandes dimensiones.

En primer lugar, se pretendía comprobar la validez del método para analizar la generación de olas por deslizamiento, tanto producidos por bloques rígidos como por masas granulares. Para ello se seleccionaron sendos ensayos de laboratorio de los existentes en la bibliografía, y se resolvieron con el PFEM.

El objetivo final es aplicar el método a casos reales, por lo que el siguiente paso ha sido simular un deslizamiento sobre la geometría a escala 1:1 de un embalse. Se trata de un caso puramente académico, del que por tanto no se pueden predecir los resultados ni realizar validaciones. El objetivo es únicamente comprobar que no hay limitaciones técnicas a priori.

Por último, se ha estudiado el caso del deslizamiento de Lituya Bay en tres dimensiones. Se trata de un ejercicio que hasta ahora no se ha realizado, ni con modelos físicos ni numéricos. Fritz et al. [34] , en uno de los artículos más recientes relativos al deslizamiento de Lituya Bay, describen algunos casos estudiados en un modelo cuasi-tridimensional, pero concluyen que está pendiente realizar un análisis completamente tridimensional (3D) del evento, con la batimetría completa de la bahía (el modelo presentado se realizaba sobre una piscina rectangular), que contemple la generación de la ola, su propagación y su ascenso sobre las laderas.

Conviene recordar las características especiales del fenómeno de los deslizamientos de ladera en embalses. Se trata de un problema en el que existe una importante incertidumbre sobre muchos de los aspectos que intervienen. En primer lugar, la aplicación práctica trata de evaluar las consecuencias de un deslizamiento posible, que por lo tanto no ha ocurrido. Es necesario definir la geometría, las características físicas (densidad, grado de fracturación, presencia de agua, etc.), y cinemáticas del fenómeno, para lo que la información geológico-geotécnica es fundamental. Si no es fácil predecir la ocurrencia de un deslizamiento, mucho más difícil es estimar cuáles van a ser sus características.

Por lo tanto, al aplicar un método numérico al fenómeno del deslizamiento de laderas en embalses, hay que tener presente que las condiciones del caso real van a ser muy diferentes de las que se puedan reproducir, tanto numéricamente como en laboratorio. Los resultados que se obtengan serán estimaciones de las posibles consecuencias de un deslizamiento potencial, a pesar de lo cual pueden resultar muy útiles para tomar las medidas de protección y aviso a la población, en su caso.

6.2. Validación bidimensional con deslizamiento granular

Para la primera validación del método se ha seleccionado la modelación del deslizamiento ocurrido en la bahía de Lituya (Alaska) en su versión bidimensional (2D). Como se explica más adelante, se trata del evento más estudiado en todo el mundo, con modelos tanto numéricos como físicos. Se ha seleccionado el modelo físico de Fritz et al. [35] , y el modelo numérico de Quecedo et al. [36] para comparar los resultados obtenidos con PFEM.

En el ensayo de Fritz [35] se analizó el deslizamiento mediante un canal prismático a escala 1:675. Teniendo en cuenta la geometría de la bahía, se consideró que la propagación lateral de la onda era poco importante, con lo que se podían esperar buenos resultados de un modelo bidimensional. La modelación numérica de Quecedo et al. [36] se desarrollo así mismo con un modelo en dos dimensiones. La figura 8 muestra la geometría esquemática en que se basó el modelo físico, que es la que se ha reproducido para el análisis con PFEM.

Figura 8. Geometría del análisis en dos dimensiones (fuente: Fritz et al. [35] ).

En el modelo físico, el deslizamiento fue modelado mediante un material granular artificial de 4 mm de diámetro medio, con una densidad seca de 2.640 kg/m³ , de acuerdo con la densidad estimada del material real.

El modelo numérico de Quecedo et al. [36] considera los tres fluidos involucrados (aire, agua y deslizamiento) por separado, como si fueran inmiscibles. Se repitió el cálculo modificando los valores de la densidad y del coeficiente de fricción, concluyendo que los resultados no varían de manera sustancial.

El cálculo con PFEM es más simplificado, porque, tal y como se menciona en el apartado 5,4, la masa deslizada se trata como un fluido viscoso (incompresible) , y no se tiene en cuenta la presencia del aire. El único parámetro con influencia significativa en los resultados es la densidad del material deslizado (se ha comprobado que la influencia de la viscosidad es despreciable). Se ha probado con valores de a) 1.000 kg/m³ (igual a la del agua), b) 1.600 kg/m³ (densidad del material de ensayo, teniendo en cuenta los huecos), y c) 2.640 kg/m³ (densidad del material de ensayo, supuesto sin huecos). Los resultados más acordes con los del modelo físico corresponden al segundo caso, como cabía esperar. En la tabla 2 se resumen los resultados de los tres modelos comparados. Las variables que se incluyen son las que se midieron en el modelo físico.

Tabla 2. Comparación de los resultados de la modelación bidimensional del deslizamiento de Lituya Bay

	Medido (Fritz)	Quecedo et al.	PFEM
Duración del deslizamiento (s)	7	9	9, 5
Velocidad de impacto (m/s)	110	85	83
Longitud del deslizamiento en el momento del impacto (m)	748	1.092	1.052
Altura de ola máxima (m)	> 200	226	234
Tiempo de ola máxima (s)	11	21	26
Coordenada x de ola máxima (m)	600	600	814
Altura de ola en x = 885 m (m)	152	266	232
Tiempo de ola máxima en x = 885 (s)	16	26,8	27,1

La figura 9 muestra distintas imágenes de los modelos comparados. Puede apreciarse que el fenómeno se reproduce cualitativamente con el modelo numérico: la formación de la ola, la creación y posterior colapso de una cavidad de aire en el trasdós de la ola, y el ascenso por la margen opuesta.

Figura 9. Resultados de los modelos físico[35] y numérico (PFEM) en instantes de tiempo correspondientes.

Las alturas de ola obtenidas con los modelos numéricos son mayores en todos los casos que las obtenidas en el ensayo experimental. Este hecho podría ser consecuencia de la no consideración de la mezcla agua-deslizamiento-aire, lo que hace que no se tenga en cuenta la disipación de energía que se produce.

Otra explicación de la discrepancia podría ser la cinemática del deslizamiento. En el modelo físico se utilizó un mecanismo para acelerar la masa y controlar su forma, mientras que en nuestro cálculo únicamente se deja caer la masa deslizada por acción de la gravedad. Ello se traduce en una diferencia en la forma y velocidad del deslizamiento en el momento del impacto.

6.3. Validación bidimensional con deslizamiento de bloques

El segundo experimento seleccionado para validación se ha tomado de Sælevik et al. [13] . Consiste en un canal prismático de 25 m de longitud, 0,51 m de anchura y 1 m de profundidad. En el instante inicial se encuentra lleno de agua hasta una altura de 0,6 m. La diferencia principal con el caso anterior estriba en la masa deslizada, que en este caso está formada por bloques rígidos que se aceleran sobre unos raíles. En la figura 10 se muestra un esquema del experimento.

Figura 10. Esquema del experimento de Sælevik et al. [13] , con la situación de los sensores.

Como en el caso anterior, se muestran los resultados mediante instantáneas comparativas (fig. 11 ).

Figura 11. Comparación entre el resultado del modelo numérico (izquierda) y el físico (derecha) en la fase de formación de la ola.

En este caso, además, se dispone del registro de la evolución de la superficie libre en tres puntos, donde se instalaron sensores (fig. 10 ). En la figura 12 se muestran los gráficos comparativos de los resultados del experimento y los obtenidos mediante PFEM.

Figura 12. Comparación entre los resultados del modelo físico (tomado de Sælevik et al. [13] , línea continua), y el cálculo con PFEM (línea discontinua) De izquierda a derecha, sensores WG1, WG2 y WG3 (ver situación en la figura 10 ).

Puede apreciarse que el modelo numérico reproduce los resultados del experimento, aunque la altura de ola máxima se subestima hasta en un 15% en el primer sensor. Hay que aclarar que, como en el caso anterior, en este experimento la velocidad de los bloques se controla con un mecanismo que los acelera, mientras que en el modelo numérico se dejan caer por acción de la gravedad. Se ha comprobado que la velocidad en el momento del impacto, que es el parámetro más determinante en la altura de ola máxima, coincide en ambos modelos. A pesar de ello, el movimiento de los bloques no es exactamente el mismo. Esta es la razón más probable de la discrepancia, especialmente la que se observa en el segundo máximo (ver gráfico correspondiente al sensor WG1 de la figura 10 ).

6.4. Caso de prueba a escala real

El siguiente ejemplo es un caso sobre la geometría real de un embalse. En esta ocasión no se dispone de datos sobre el deslizamiento ni sus consecuencias. Se ha tomado la cartografía real de un embalse, incluida la presa, y se ha definido una zona inestable de un tamaño del orden de magnitud de la presa, de modo que se produjera una ola que generase sobrevertido.

El objetivo era comprobar que en un caso de estas dimensiones era posible aplicar el método, salvando los problemas de coste computacional que normalmente presentan los códigos de cálculo tridimensionales. Por otra parte, se quería estudiar el mejor método para trabajar con la información topográfica en los formatos más comúnmente utilizados (shp, dxf, etc.), para convertirla en la malla de elementos finitos que sirve de base para el cálculo con PFEM. Este proceso se ha realizado de una manera rápida y eficaz con el progama GiD (www.gidhome.com ).

En la figura 13 se muestran algunas imágenes del resultado obtenido.

Figura 13. Deslizamiento hipotético a escala real. Se observa la generación de la ola, su propagación, y el sobrevertido sobre la presa, que en este caso es de materiales sueltos (derecha).

Uno de los objetivos actuales es poder realizar cálculos acoplados en los que se considere simultáneamente el deslizamiento, con toda la complejidad que se ha descrito, y la afección a la presa. PFEM ya permite la posibilidad de estudiar la erosión de un material por la acción del agua [7]  and [39] . Este módulo está siendo desarrollado para relacionar los parámetros con los que se define la resistencia a la erosión en el modelo numérico con aquellos que se utilizan habitualmente para la caracterización de cada tipo de material.

6.5. Simulación tridimensional del deslizamiento de la bahía de Lituya

De entre todos los eventos de generación de olas por deslizamiento, uno de los más conocidos en la comunidad científica y que ha sido objeto de más profundo análisis es el deslizamiento de Lituya Bay (Alaska). El 9 de julio de 1958, se produjo un deslizamiento de unos 90 millones de toneladas de rocas, generado por un gran terremoto, que se precipitaron sobre la bahía y formaron una ola que alcanzó una altura máxima sobre la ladera opuesta de 524 m. En la figura 14 se muestra un esquema de la geometría del deslizamiento. Muchos investigadores han estudiado este caso, además de por su magnitud, por el hecho de que puede estimarse con cierta precisión el alcance de la ola en los diferentes puntos de la bahía a partir de la marca de destrucción que dejó en las laderas, al arrasar la capa vegetal que alcanzó. Además, por la geometría del deslizamiento y de la bahía en la zona de impacto, en general se ha considerado aceptable estudiar el problema en dos dimensiones. Por último, la realización de ensayos bidimensionales en modelo reducido bien documentados ha permitido a los desarrolladores de códigos numéricos utilizar los resultados para validación.

Figura 14. Esquema del deslizamiento de Lituya Bay. La destrucción del bosque de la ladera izquierda corresponde con la cota alcanzada por la ola generada (Fritz et al. [35] ).

En la revisión bibliográfica que se ha realizado no se ha encontrado un trabajo en el que se describa una modelación tridimensional del deslizamiento de Lituya Bay. Este es el objetivo planteado en esta investigación.

Como base de partida se tomó el trabajo de Fritz et al. [35] , en el que se recopila con detalle gran cantidad de información sobre el citado evento.

La cartografía se ha obtenido de la página web del U. S. Geological Survey (http://eros.usgs.gov/ ). Tras algunas transformaciones utilizando herramientas GIS, se obtienen las curvas de nivel de la figura 15 , a partir de las cuales se crea la malla de elementos finitos mediante el programa GiD. En este caso, resulta inútil buscar una precisión mayor, porque hay otras incertidumbres sobre la situación en el momento del deslizamiento, como la propia geometría de la masa deslizada, la forma del frente del glaciar en el momento del suceso, o la batimetría. Por ejemplo, Fritz et al. [35] desprecian el efecto de los depósitos existentes junto al frente del glaciar.

Figura 15. Vista en planta del deslizamiento de Lituya Bay. Izquierda: Cartografía utilizada en el modelo; curvas de nivel cada 100 m. Derecha: imagen de Fritz[35] en la que se marca el deslizamiento.

La geometría inicial del deslizamiento se ha generado a partir de la información de las figuras 8 , 14 y 15 . Hay que resaltar que el único parámetro del modelo sobre el que se ha trabajado aparte de la geometría del deslizamiento ya comentada, es el tamaño de malla. En los resultados presentados se ha utilizado una malla de 30 m de dimensión característica. Todas las simulaciones se han realizado con una densidad de la masa deslizada de 1.600 kg/m³ , que es la que mejores resultados dio en el caso bidimensional.

En las figuras 16 y 17 se muestran algunas instantáneas en las que se observan los resultados del cálculo. Se han comparado con los datos existentes relativos a la máxima cota alcanzada por la ola, así como con el plano de la afección producida (fig. 18 ).

Figura 16. Formación de la ola, propagación, y alcance máximo en la margen norte. Curvas de nivel cada 100 m.

Figura 17. Propagación de la ola en la bahía y alcance máximo en la margen sur. Curvas de nivel cada 100 m.

Figura 18. Comparación de los datos de Fritz[35] (izquierda) con los resultados del cálculo (derecha; curvas de nivel cada 100 m). Se muestra una superposición de imágenes de los resultados en distintos pasos de tiempo.

El alcance máximo obtenido sobre la ladera opuesta con PFEM es de 551 m (fig. 16 ), algo superior al observado (524 m), lo que supone un error relativo del 5%. El punto en el que se alcanza este máximo difiere en unos 300 m del indicado por Fritz. En dicha zona, justo en la arista de la colina situada frente al deslizamiento, el alcance obtenido con el modelo numérico es de 492 m (un 6% menor). En cuanto a la ladera sur de la bahía, el máximo alcance fue de 208 m, mientras que con el PFEM se ha obtenido un valor de 195 m (figura 17 ). Esto supone un error relativo del 6%. Conviene señalar que las discrepancias entre lo observado y los resultados del modelo son del mismo orden de magnitud que la malla de elementos finitos empleada.

7. Conclusiones y líneas de investigación futuras

Se han presentado los resultados obtenidos en la aplicación del PFEM al análisis de deslizamientos de ladera en embalses. La masa deslizada se ha considerado, dependiendo del caso, como un sólido rígido o como un fluido denso. Hay que tener en cuenta que se trata de un fenómeno que, además de ser muy complejo, presenta una gran incertidumbre en muchas de las variables que lo gobiernan. Se trata de una característica común a todos los modelos que tratan de reproducir el terreno, puesto que en general la información de que se dispone es, cuando menos, discreta, y siempre debe extrapolarse a un continuo. En este caso, la geometría y características mecánicas y cinemáticas del deslizamiento son dos aspectos clave que deben estimarse.

Los modelos numéricos, como el PFEM, además deben basarse en simplificaciones del fenómeno, por lo que sus resultados sólo pueden considerarse aproximaciones a la realidad. En los casos presentados no se ha tenido en cuenta la presencia del aire ni su mezcla con el agua. Por otra parte, al analizar casos a escala real, se desprecian las irregularidades del terreno, así como el fenómeno de erosión y el arrastre de nuevo material originado por el propio deslizamiento, lo que en general debería conducir a errores en el alcance de la afección.

Otras fuentes de error provienen del desconocimiento del comportamiento real del material desestabilizado. En general, se desconoce hasta qué punto un deslizamiento en bloque puede romperse antes del impacto. Es decir, si su movimiento puede pasar de un deslizamiento puro a una avalancha completamente desarrollada antes de impactar con la masa de agua.

No obstante todo lo anterior, los resultados del análisis 3D del deslizamiento de Lituya Bay permiten concluir que PFEM puede ser útil en la estimación de las afecciones de deslizamientos potenciales.

En la actualidad, los autores están trabajando en la optimización del método en relación con el tiempo de cálculo, lo que permitirá en el futuro afrontar problemas de gran tamaño con mallas más finas[38]  and [39] . Así mismo, se está perfeccionando el módulo de erosión, con el objetivo de simular fenómenos como el aumento de masa que se incorpora al deslizamiento durante la caída, o el efecto del sobrevertido producido por la ola sobre la estabilidad de la presa, especialmente si esta es de materiales sueltos. Todos estos fenómenos pueden reproducirse con el PFEM y serán objeto de futuras investigaciones.

Agradecimientos

Los autores expresan su agradecimiento por la ayuda prestada para el desarrollo del presente trabajo a Daniel Rodríguez Martínez, Miguel Ángel Celigueta y Pedro Arnau.

Este trabajo se ha realizado en el marco de los proyectos SEDUREC del Programa Consolider, XLIDE del Programa INNPACTO, ambos del MICINN y SAFECON del European Research Council de la Comisión Europea.

References

1. [1] D.J. Varnes; Slope movements types and processes; Landslides: analysis and control, Transportation Research Board, Special report (1978) 11–33
2. [2] F.J. Ayala, Análisis de riesgos por movimientos de ladera. Riesgos Naturales, Ayala, F. J. y Olcina, J. (coordinadores). ISBN 84-344-8034-4, (2002) 379–410.
3. [3] E.E. Alonso, N.M. Pinyol, Landslides in reservoirs and dam operation, en: R. Romeo García et al. (Eds.), Dam Maintenance and Rehabilitation II, CRC Press, Londres, 2010, pp. 3–27.
4. [4] T. Kubota; Study on the fractal dimensions and geological condition of landslides; Chacón, Iriagaray, Fernández (Eds.), Landslides, Balkema, Rotterdam (1996), pp. 385–392
5. [5] S.R. Idelsohn, E. Oñate, F. Del Pin, N. Calvo; Fluid-structure interaction using the particle finite element method; Computer Methods in Applied Mechanics and Engineering, 195 (2006), pp. 2100–2123
6. [6] E. Oñate, S.R. Idelsohn, F. Del Pin, R. Aubry; The particle finite element method: an overview; International Journal on Computational Methods, 1 (2004), pp. 267–307
7. [7] E. Oñate, S.R. Idelsohn, M.A. Celigueta, R. Rossi; Advances in the particle finite element method for the analysis of fluid-multibody interaction and bed erosion in free surface flows; Comp. Methods in Appl. Mech. and Eng., 197 (2008), pp. 1777–1800
8. [8] A. Larese, R. Rossi, E. Oñate, S.R. Idelsohn; Validation of the particle finite element method (PFEM) for simulation of free surface flows; Int. J. for Computer-aided Engineering and Software, 25 (4) (2008), pp. 385–425
9. [9] B. Haddad, Modelización numérica mediante elementos finitos y SPH de los geomateriales fluidificados: aplicación a los deslizamientos rápidos de ladera. Tesis doctoral. ISBN: 978-84-669-3128-1 (2007).
10. [10] Pierson; Erosion and deposition by debris flows at Mt. Thomas, North Canterbury, New Zealand; Earth Surf. Processes, 5 (1986), pp. 227–247
11. [11] J. Sander, K. Hutter; Evolution of weakly non-linear shallow water waves generated by a moving boundary; Acta Mechanica, 91 (1992), pp. 119–155
12. [12] M. Di Risio, G. Bellotti, A. Panizzo, P. De Girolamo; Three-dimensional experiments on landslide generated impulse waves at a sloping coast; Coastal Engineering, 56 (2009), pp. 659–671
13. [13] G. Sælevik, A. Jensen, G. Pedersen; Experimental investigation of impact generated tsunami; related to a potential rock slide, Western Norway; Coastal Engineering., 56 (2009), pp. 897–906
14. [14] H.M. Fritz, P. Moser; Pneumatic landslide generator; Int. J. Fluid Power, 4 (1) (2003), pp. 49–57
15. [15] E. Noda; Water waves generated by landslides; J. Waterways, Harbors and Coastal Engineering, Div, 96 (WW4) (1970), pp. 835–855
16. [16] M. Di Risio, P. Sammarco; Analytical modeling of landslide generated waves; Journal of Waterway, Port, Coastal and Ocean Engineering-ASCE, 134 (1) (2008), pp. 53–60
17. [17] J.J. Monaghan, A. Kos; Scott Russells wave generator; Physics of Fluids, 12 (3) (2000), pp. 622–630
18. [18] M. Pastor, I. Herreros, J.A. Fernández Merodo, P. Mira, B. Haddad, M. Quecedo, E. González, C. Álvarez-Cedrón, V. Drempetic; Modelling of fast catastrophic landslides and impulse waves induced by them in fjords, lakes and reservoirs; Engineering Geology, 109 (2009), pp. 124–134
19. [19] H.M. Fritz, W.H. Hager, H.E. Minor; Near field characteristics of landslide generated impulse waves; Journal of Waterway, Port, Coastal and Ocean Engineering. ASCE, 130 (6) (2004), pp. 287–302
20. [20] A. Huber, W.H. Hager; Forecasting impulse waves in reservoirs. 18 ICOLD Congress; Florencia, C31 (1997), pp. 993–1005
21. [21] D.L. Vischer, W.H. Hager; Dam hydraulics; Ed. Wiley & Sons (1998) ISBN 978 0471972891
22. [22] L. Garrote, F. Laguna, Generación de olas por deslizamiento. Jornadas técnicas sobre estabilidad de laderas en embalses. Zaragoza (2007).
23. [23] S.R. Idelsohn, E. Oñate, F. Del Pin; A Lagrangian meshless finite element method applied to fluid-structure interaction problems; Comput. Struct., 81 (2003), pp. 655–671
24. [24] S.R. Idelsohn, N. Calvo, E. Oñate; Polyhedrization of an arbitrary point set; Comput. Meth. Appl. Mech. Engrg., 192 (22–24) (2003), pp. 2649–2668
25. [25] S.R. Idelsohn, E. Oñate, F. Del Pin; The particle finite element method: a powerful tool to solve incompressible flows with free-surfaces and breaking waves; Int. J. Numer. Methods Engrg., 61 (2004), pp. 964–989
26. [26] R. Aubry, S.R. Idelsohn, E. Oñate; Particle finite element method in fluid mechanics including thermal convection-diffusion; Computer and Structures, 83 (2005), pp. 1459–1475
27. [27] E. Oñate, S.R. Idelsohn, F. Del Pin; Lagrangian formulation for incompressible fluids using finite calculus and the finite element method; Y. Kuznetsov, P. Neittanmaki, O. Pironneau (Eds.), Numerical Methods for Scientific Computing Variational Problems and Applications, CIMNE, Barcelona (2003)
28. [28] S.R. Idelsohn, E. Oñate, E.N. Calvo, F. Del Pin; The meshless finite element method; Int. J. Num. Meth. Engng., 58 (6) (2003), pp. 893–912
29. [29] R. Codina, O.C. Zienkiewicz; CBS versus GLS stabilization of the incompressible Navier-Stokes equations and the role of the time step as stabilization parameter; Communications in Numerical Methods in Engineering, 18 (2) (2002), pp. 99–112
30. [30] E. Oñate, J. García, S.R. Idelsohn, F. Del Pin; Finite calculus formulations for finite element analysis of incompressible flows. Eulerian, ALE and Lagrangian approaches; Computer Methods in Applied Mechanics and Engineering, 195 (23–24) (2006), pp. 3001–3037
31. [31] E. Oñate; Possibilities of finite calculus in computational mechanics; Int. J. Num. Meth. Engng., 60 (1) (2004), pp. 255–281
32. [32] E. Oñate, J. García; A finite element method for fluid-structure interaction with surface waves using a finite calculus formulation; Comput. Meth. Appl. Mech. Engrg., 191 (2001), pp. 635–660
33. [33] H. Edelsbrunner, E.P. Mucke; Three dimensional alpha shapes; ACM Trans. Graphics, 13 (1999), pp. 43–72
34. [34] H.M. Fritz, F. Mohammed, J. Yoo; Lituya Bay landslide impact generated mega-tsunami 50th anniversary; Pure appl. Geophys., 166 (2009), pp. 153–175
35. [35] H.M. Fritz, W.H. Hager, H.E. Minor; Lituya Bay Case: Rockslide impact and wave run-up; Science of Tsunami Hazards, 19 (1) (2001)
36. [36] M. Quecedo, M. Pastor, M.I. Herreros; Numerical modelling of impulse wave generated by fast landslides; International Journal for Numerical Methods in Engineering, 59 (2004), pp. 1633–1656
37. [37] S.R. Idelsohn, J. Marti, A. Limache, E. Oñate; A unified Lagrangian formulations for both: the elastic solid and the incomopressible fluid; Computer Methods in Applied Mechanics and Engineering, 197 (19–20) (2008), pp. 1762–1776
38. [38] E. Oñate, S.R. Idelsohn, M.A. Celigueta, R. Rossi, J. Marti, J.M. Carbonell, P. Ryzhakov, B. Suárez; Advances in the Particle Finite Element method (PFEM) for solving coupled problems in engineering; E. Oñate, R. Owen (Eds.), Particle-Based Methods, Vol. 25, Springer (2011), pp. 1–49
39. [39] E. Oñate, S.R. Idelsohn, M.A. Celigueta, R. Rossi, S. Latorre, Possibilities of the particle finite element method in fluid-soil-structure interaction problems. Submitted to Computational Mechanics (2011).
40. [40] M.B. Liu, G.-R. Liu; Smoothed particle hydrodynamics (SPH): an overview and recent developments; Archives of Computational Methods in Engineering, 17 (1) (2010), pp. 25–76