Funciones de estado límite

Atendiendo los modos de falla reportados en la literatura se propone la evaluación de la probabilidad de falla bajo dos escenarios principales, independientes entre sí: el primero es el agrietamiento en el cuerpo de la cortina y el segundo es el deslizamiento en su base, aunque este último puede ser utilizado para evaluar el deslizamiento en cualquier plano débil de falla en la roca de cimentación.

Sea Gckx- la función de estado límite expresada en función de las variables inciertas x- involucradas en cierto modelo de la cortina, entonces, el margen de seguridad para el agrietamiento puede expresarse en términos del índice de reducción de rigidez secante de la estructura ( ISSR ) presentado por ^{Esteva et al. (2011)} como sigue:

donde ISSR=K0-KK0 . El término K0 es la rigidez tangente global inicial de la cortina obtenida a través de un análisis dinámico considerando una intensidad sísmica poco significativa (que no provoque comportamientos no lineales) y K es la rigidez secante calculada con el valor máximo absoluto del desplazamiento relativo de la corona de la cortina con respecto a su base, durante la respuesta no lineal de la cortina y su cortante basal asociado, expresada mediante ciclos de histéresis. En este trabajo los ciclos de histéresis se obtuvieron mediante modelos bidimensionales utilizando el método de Elemento Finito (MEF2D) considerando un comportamiento no lineal en el concreto.

La condición de falla se expresa como ISSR~1 . La confiabilidad se determina ajustando una función lognormal de probabilidad a todos los valores mínimos de la intensidad sísmica de falla ( YF ) que satisfacen esa condición. Haciendo ZF=ln⁡YF la función de confiabilidad puede expresarse como (^{Esteva et al., 2011}):

En la ecuación anterior E(∙) y σ(∙) representan el valor esperado y la desviación estándar, respectivamente. Para el deslizamiento de la cortina se consideran dos criterios para establecer la función de estado límite. El primero está basado en ^{Bernier et al. (2016)} quienes mencionaron los niveles de daños relacionados con la magnitud del deslizamiento en la base de una cortina de concreto, basado en el juicio de expertos. Estos investigadores señalaron que un deslizamiento de 150 mm implicaría un nivel de daño severo que causaría movimientos diferenciales inaceptables en relación con los monolitos adyacentes y podría, eventualmente, causar el colapso total y, por consiguiente, la liberación descontrolada del agua de la presa. Bajo estas consideraciones, la función de estado límite se expresa con la ecuación siguiente:

Donde Slb es el deslizamiento en la base de la cortina obtenido del análisis MEF2D.

El segundo criterio para establecer la función de estado límite para el deslizamiento de la cortina es mediante la consideración de un modelo de cuerpo rígido (MCR), más simple que el modelo MEF2D, pero usado ampliamente en la práctica profesional durante las etapas de diseño preliminar, basado en un modelo de falla tipo Mohr-Coulomb. Se incluye una discretización de la base de la cortina, para considerar variaciones en la cohesión y el ángulo de fricción conforme a los hallazgos de ^{Westberg y Johansson (2013)} quienes coordinaron la perforación y obtención de 8 núcleos del contacto cortina-cimentación en la presa Laxeden, Suecia. De las 8 muestras, 5 resultaron intactas (con cohesión) y 3 mostraron el contacto concreto-roca con agrietamiento (sin cohesión). En este caso, la función de estado límite para el deslizamiento Gsl2x- queda expresado, partiendo de un modelo Mohr-Coulomb, como sigue:

En la ecuación anterior k representa el número total de discretizaciones consideradas en la base de la cortina, ci es la cohesión, Ai es el área de la discretización, Ni es la fuerza normal actuando en la discretización, μi es el coeficiente de fricción calculado como μi=tan⁡ϕi , siendo ϕi el ángulo de fricción y el término ∑Fh indica la sumatoria de todas las fuerzas horizontales actuando en la cortina (presión hidrostática, presión de azolves, fuerza inercial horizontal, presión hidrodinámica y presión dinámica de azolves); el subíndice i indica la discretización correspondiente. La congruencia con las observaciones de ^{Westberg y Johansson (2013)} se dará si se considera que cuando más el 60% de los tramos resultantes de la discretización de la base de la cortina tiene cohesión y fricción y el 30% restante únicamente fricción.

Modelos de comportamiento estructural

La geometría de la cortina y las solicitaciones impuestas se muestran en la figura 5. Se consideró una altura de cortina de 40 metros y un nivel fijo en el embalse (h) a 36 metros sobre el nivel del fondo del cauce y sus respectivas masas adheridas para representar los efectos hidrodinámicos. Se supuso una eficiencia del 55% en el sistema de drenaje y un espesor de 10 metros de azolves. La sección geométrica de la cortina fue diseñada satisfaciendo los factores de seguridad ante volteo, deslizamiento y esfuerzos excesivos en su base, conforme a las recomendaciones del Manual de Obras Civiles de la Comisión Federal de Electricidad (^{MOC, 2015}). Se utilizó como coeficiente sísmico horizontal kh=a0/(1+2a0) actuando en la dirección aguas abajo, con a0 definida como la aceleración máxima en roca (PGA, por sus siglas en inglés), asociada al nivel de intensidad sísmica de diseño (probabilidad anual de excedencias 1/10,000). Adicionalmente se consideró un coeficiente sísmico vertical tomado como kv=2kh/3 .

Figura 5 Características de la cortina analizada en este estudio 

La evaluación de las funciones de estado límite fue desarrollada mediante dos modelos de comportamiento estructural: MEF2D y MCR. Con el modelo MEF2D se desarrolló un análisis dinámico paso a paso en el dominio del tiempo a través del software ABAQUS, considerando un comportamiento no lineal en el concreto mediante un modelo continuo de daño por falla a tensión y compresión que se inicia cuando los esfuerzos principales alcanzan las resistencias ante ambas condiciones. Durante el comportamiento no lineal a tensión se asumió el criterio de agrietamiento mediante la energía de fractura propuesto por Hillerborg’s en 1976, considerando una curva bilineal de ablandamiento por deformación. Para compresión, se utilizó la parábola de Hognestad para representar la curva esfuerzo-deformación en el rango no lineal. Las propiedades mecánicas y de resistencia de los materiales de la cortina para la construcción de las relaciones constitutivas mencionadas anteriormente se presentan en la tabla 1, descrita en los apartados siguientes. Las historias de aceleración utilizadas en los análisis tuvieron un tiempo total de 80 segundos.

Por otro lado, la base de la cortina fue discretizada en 10 secciones de las que 6 incluyeron cohesión y fricción (C+F) en el contacto con la cimentación y 4 incluyeron únicamente fricción (F), tomando como base las observaciones efectuadas por ^{Westberg y Johansson (2013)}. Esta consideración fue ejecutada tanto en el modelo MEF2D como en el modelo MCR (ver figura 6).

Figura 6 Ejemplo de discretización de la base de la cortina en el modelo de comportamiento estructural 

Caracterización de la amenaza sísmica

La estimación del desempeño estructural de una cortina rígida, de acuerdo con las funciones de estado límite descritas anteriormente, requiere el uso de historias de aceleración en el tiempo (acelerogramas) relacionadas con el movimiento sísmico. Los valores numéricos de estos acelerogramas pueden obtenerse de una base de datos de registros reales o generarse sintéticamente. Es necesario además que los acelerogramas estén asociados a un nivel de intensidad sísmica, pero en varias ocasiones las bases de datos no cuentan con registros suficientes para relacionarlos con el nivel de intensidad de interés. Por lo anterior, es común emplear técnicas para generar acelerogramas sintéticos mediante el escalamiento de registros reales y llevarlos hasta el nivel de intensidad sísmica requerida.

A diferencia del método de escalamiento de registros sísmicos a través de un escalar descrito por ^{Vamvatsikos y Cornell (2002)}, en este trabajo se considera la aplicación del método híbrido de simulación de historias del movimiento del terreno propuesto por ^{Ismael y Esteva (2006)} que aprovecha las ventajas las funciones generalizadas de atenuación (^{Alamilla et al., 2001}) y las funciones de Green empíricas (^{Ordaz et al., 1995}). Con lo anterior se pretende que los registros sísmicos sean escalados tanto en términos de intensidad sísmica como en contenido de frecuencias. Para lograr lo anterior, se procesaron 455 registros sísmicos de 7 estaciones ubicadas en el estado de Chiapas, México, que ocurrieron entre los años 1976 a 2016. Los registros sísmicos están contenidos en la base de datos de registros acelerográficos de la red sísmica del Instituto de Ingeniería de la UNAM.

En la figura 7 se muestra la ubicación de los epicentros, las estaciones acelerográficas y un sitio elegido para caracterizar la amenaza sísmica mediante el método híbrido mencionado anteriormente. Con los datos de magnitud, distancia y aceleraciones de los 455 registros sísmicos y la integración numérica de la ecuación 5 se obtuvo la tasa de excedencias de las intensidades sísmicas vYy en el sitio de interés (ver figura 8).

Figura 7 Epicentros y estaciones acelerográficas utilizadas para caracterizar la amenaza sísmica. Imagen de Global Mapper 

Figura 8 Tasa anual de excedencias de PGA como resultado del procesamiento de los 455 registros sísmicos en el estado de Chiapas 

De acuerdo con el procedimiento propuesto por ^{Alamilla et al. (2001)}, se estimaron la magnitud M y distancia R más probables del evento que generaría la intensidad asociada a cierta probabilidad anual de excedencias, partiendo de la función de densidad de probabilidad conjunta de M y R para un valor dado de la intensidad sísmica Y = y, misma que fue obtenida a partir del teorema de Bayes y que está representada en la ecuación 11.

En la ecuación anterior, k es una constante de normalización tal que la integral de fM,Rm,ry sobre el dominio de valores posibles de M y R es igual a la unidad, fYym,r es la función de densidad de probabilidad condicional de la intensidad sísmica Y para M=m y R=r y fM,Rm,r es la función de densidad de probabilidad conjunta de M y R. De forma inicial, se obtuvo un valor aleatorio de M condicionado a la intensidad sísmica de interés mediante la ecuación 12.

Para obtener la densidad de probabilidad f _Y (y|m,r) se consideró la función de densidad de probabilidad del error estadístico de la función de atenuación considerada en la ecuación 5. Este error estadístico se obtiene a través de la relación y1/ yc, donde y1 es la aceleración observada y yc la aceleración obtenida con la regresión. Esta relación se ajustó a una función de probabilidad lognormal Log[-0.0076, 0.917], como se muestra en el histograma de la figura 9.

Figura 9 Función de densidad de probabilidad de la relación y1/yc[/p]  

La función de densidad de probabilidad conjunta fM,Rm,r se tomó como el producto de las probabilidades marginales fM(m)fR(r) . La función de densidad de probabilidad fM(m) se representa mediante la expresión dada por ^{Alamilla et al. (2001)} presentada en la ecuación 13, considerando Mo, Mu y β tal y como fueron definidas en la ecuación 5.

La densidad de probabilidad de la distancia, fR(r) se obtuvo mediante inferencia estadística, a través del ajuste de una función de probabilidad a las distancias focales observadas, de acuerdo con el procesamiento geográfico de esas distancias al sitio de interés, para los 455 registros sísmicos mostrados en la figura 7. Para este caso, la función de probabilidad encontrada es una función lognormal con parámetros Log[5.114, 0.335]. En la figura 10 se muestra el histograma con la función de densidad de probabilidad ajustada para fR(r) , a partir de los datos observados.

Figura 10 Función de densidad de probabilidad fR(r) de las distancias focales al sitio de interés, obtenida con los datos observados procesados mediante el uso de un Sistema de Información Geográfica 

Para cada valor supuesto de la intensidad y se simuló un valor de M, de acuerdo con la ecuación 12, y se simuló un valor de la distancia focal R, considerando la función de densidad de probabilidad de la distancia condicionada a la magnitud y la intensidad sísmica de interés, expresada matemáticamente en la ecuación 14.

Considerando el procedimiento anterior se llevaron a cabo 1,000 simulaciones de la magnitud y distancia, asociados a distintos niveles de intensidad sísmica. Se escribió un código de programación en el lenguaje Wolfram-Mathematica con el fin de realizar la integración numérica del método híbrido. En la figura 11 se muestra la función de probabilidad conjunta fM,R(m,r|y) para las intensidades asociadas a probabilidades de excedencia de 1/10,000, 1/3,000, 1/1,000 y 1/145. Finalmente, para la generación de los acelerogramas sintéticos incluyendo las tres direcciones ortogonales ( ahmax , ahmin y avert ) se utilizó el método de las funciones de Green empíricas descrito en ^{Ordaz et al. (1995)}, utilizando 6 registros semilla. En la figura 12 se muestra el ejemplo de los acelerogramas sintéticos generados con el método híbrido para las tres componentes ortogonales, asociados al nivel de intensidad sísmica correspondiente a una probabilidad de excedencias del 50% en 100 años. En última instancia, en la figura 13 se presentan los espectros de respuesta elástica de 10 acelerogramas sintéticos generados para un nivel de intensidad sísmica asociada a una probabilidad de excedencias del 50% en 100 años.

Figura 11 Función de densidad de probabilidad conjunta fM,R(m,r|y) como resultado del procesamiento de 455 sismos en el estado de Chiapas 

Figura 12 Acelerogramas sintéticos de las tres componentes ortogonales: ahmax , ahmin y. avert a) Semilla, b) Función de Green en sitio y c) Acelerograma escalado a la intensidad sísmica de interés 

Figura 13 Espectros de respuesta elástica de 10 acelerogramas sintéticos generados con el método híbrido propuesto por ^{Ismael y Esteva (2006)}, asociados a una probabilidad de excedencias del 50% en 100 años y un porcentaje de amortiguamiento del 5% con respecto al crítico 

Representación probabilista de las propiedades mecánicas y de resistencia de los materiales de la cortina.

Durante el análisis estructural se consideraron 7 parámetros de forma probabilista de los que 5 corresponden a la caracterización del concreto y 2 al contacto cortina-cimentación. Todos los parámetros fueron representados por medio de funciones lognormales de probabilidad Log[λ,ζ]. La tabla 1 resume los parámetros de las funciones lognormales utilizadas. Tanto la relación de Poisson υ como el ángulo de dilatancia ψ del concreto se tomaron de forma determinista, con valores de 0.20 y 36°, respectivamente.

También se tomaron de forma determinista las propiedades de la roca de cimentación como su densidad, módulo de elasticidad y relación de Poisson. El módulo de corte (G) del contacto cortina-cimentación se obtuvo conforme a la expresión G=Ec/2(1+υ) , donde E _c es el módulo de elasticidad del concreto.

La resistencia a compresión del concreto ( f'c ) fue correlacionada con su resistencia a tensión (ft) y su módulo de elasticidad (E _c ). Para encontrar la correlación entre las resistencias a compresión y tensión se utilizó la ecuación 15. La correlación se estableció a través de 459 pruebas de resistencia a compresión simple y tensión indirecta de una presa construida en el norte de México. Se obtuvo que ϕ _t se caracteriza mediante una función de densidad de probabilidades lognormal Log[-1.03,0.29]. Para el módulo de elasticidad se utilizó una expresión similar a la de la ecuación 15, con ϕ _E ajustada a una función Log[8.63,0.16].

Para cada simulación se eligió aleatoriamente un arreglo para las discretizaciones en la base de la cortina, por lo que las zonas con cohesión y fricción (C+F) y las zonas con fricción (F) tuvieron posiciones distintas a las que se muestran en la figura 6. Lo anterior permite considerar la variación espacial de la cohesión en el contacto analizado conforme a las observaciones efectuadas por ^{Westberg y Johansson (2013)}.

El muestreo de las propiedades mecánicas y de resistencia de la cortina se efectuó con el método del Hipercubo Latino (HL). A través de este método puede reducirse el número de simulaciones utilizadas para la estimación de la probabilidad de falla pues, a diferencia del método de Monte Carlo, el método del Hipercubo Latino (HL) involucra un esquema estratificado para mejorar la cobertura del espacio de la muestra; la estratificación se cumple dividiendo el eje vertical de la gráfica de la función de distribución de probabilidad de una variable aleatoria X _j en n intervalos de igual longitud.

Probabilidad de falla por agrietamiento de la cortina. Modelo MEF2D

La estimación de la probabilidad de falla ante el agrietamiento del cuerpo de la cortina se llevó a cabo mediante un análisis dinámico incremental (IDA), utilizando 3 grupos de 10 vectores de propiedades mecánicas y de resistencia conformados aleatoriamente con el método de muestreo del Hipercubo Latino (HL), utilizando las funciones lognormales de probabilidad con los parámetros presentados en la tabla 1. En total se desarrollaron 176 análisis dinámicos de los que se obtuvieron 30 valores de la intensidad de falla; 19 correspondieron a una falla por agrietamiento en la dirección aguas abajo y 11 en la dirección aguas arriba.

Los resultados de interés del modelo MEF2D son la historia de desplazamientos δx(t) en los puntos A y B, así como la fuerza cortante en la base VX(t) representados en la figura 14. En la figura 15 se muestra un ejemplo de los desplazamientos relativos entre los puntos A y B, calculados como la sustracción definida por δxA(t)-δxB(t) , así como entre los puntos B y C, estimados como δxB(t)-δxC(t) , ambos en milímetros, para cada intervalo de tiempo t. Los desplazamientos relativos positivos son en la dirección aguas debajo de la cortina mientas que los negativos son en la dirección aguas arriba. El comportamiento estructural de la cortina puede resumirse a partir de la figura 16 y figura 17 que muestran el resultado de los 176 análisis desarrollados. En estas figuras ∆xA-xB y ∆xB-xC son el logaritmo natural de los desplazamientos relativos máximos encontrados en cada una de las series δxA(t)-δxB(t) y δxB(t)-δxC(t) , respectivamente. La variable η que se muestra en las figuras 16 y 17, representa la demanda estructural, pues está calculada como la relación entre las pseudoaceleraciones espectrales correspondientes al periodo fundamental de la cortina ( Sa ) de cada simulación y la pseudoaceleración espectral media de falla ( Saμ ) obtenida de igual manera, pero considerando los valores medios de las propiedades mecánicas y de resistencia presentados en la tabla 1.

Figura 14 Resultados de interés en el modelo MEF2D para la estimación de la probabilidad de falla ante el agrietamiento de la cortina y el deslizamiento de su base 

Figura 15 Ejemplo de historia de desplazamientos relativos entre los puntos A-B y B-C mostrados en la figura 14  

Figura 16 Resultados del desplazamiento relativo de la corona de la cortina (A) con respecto a su base (B) y la demanda sísmica (η) para los 176 análisis desarrollados 

Figura 17 Resultados del desplazamiento relativo de la base de la cortina (B) con respecto a la cimentación (C) y la demanda sísmica (η) para los 176 análisis desarrollados 

Para cada análisis se obtuvo el índice de reducción de rigidez secante expresado como IRRS=(K0-K)/K0 . El ejemplo de los ciclos de histéresis cuando la cortina exhibe ambos comportamientos, lineal y no lineal, se presenta en la figura 18. La rigidez inicial K0 está calculada como el cociente entre la fuerza cortante y el desplazamiento histerético máximo asociado, cuando la cortina exhibe un comportamiento lineal, ante intensidades sísmicas bajas. La rigidez secante K está calculada de la misma forma, pero cuando la cortina exhibe un comportamiento no lineal ante intensidades sísmicas significativas y los desplazamientos han producido grietas. Para diferenciar las rigideces de la cortina en las direcciones aguas abajo y aguas arriba se incluyó un signo + y - como subíndice, respectivamente. Nótese en el ejemplo de la figura 18 que la cortina exhibe un comportamiento no lineal más significativo en el dirección aguas abajo, por lo que las grietas están formadas en el paramento aguas arriba, mientras que en la dirección aguas arriba la rigidez secante K(-) es muy parecida a la rigidez secante K0(-) , exhibiendo un comportamiento lineal en esta dirección. Los valores del índice IRRS calculados para los 176 análisis dinámicos se muestran en la figura 19.

Figura 18 Ejemplo de dos ciclos de histéresis. Uno exhibiendo comportamiento lineal y otro exhibiendo comportamiento no lineal en el cuerpo de la cortina 

Figura 19 Valores del índice de reducción de rigidez secante ISSR para las direcciones aguas arriba (+) y aguas abajo (-) de la cortina 

Teóricamente la condición de falla se cumple cuando la rigidez secante K=0, lo que implica que IRRS=1 . Sin embargo, para valores de IRRS ligeramente menores que 1 se observaron daños considerables en el modelo MEF2D, con patrones de agrietamiento que atraviesan completamente el cuerpo de la cortina, como se muestra en la figura 20. Considerando lo anterior, se registraron las intensidades sísmicas mínimas de falla ( YF ) que produjeron un valor de IRRS mayor o igual a 0.85.

Figura 20 Ejemplo de patrones de agrietamiento y valores calculados de ISSR con base en los ciclos de histéresis obtenidos con el modelo MEF2D 

Para la estimación del índice de confiabilidad se requiere la aplicación de la ecuación 8, para las intensidades sísmicas o parámetros de demanda que satisfacen la condición de IRRS≥0.85 . En la tabla 2 se presentan los parámetros EZF y σZF obtenidos mediante un análisis de máxima verosimilitud, pero haciendo que YF=η , por lo que ZF=ln⁡[η] . Se consideró a PGA (aceleración máxima en roca) como una medida de intensidad sísmica y η como un parámetro de demanda. También, se evaluaron tres casos de falla, a saber:

Las funciones de confiablidad se presentan en figura 21 para los tres casos de falla descritos anteriormente. En esta figura β es el índice de confiabilidad propuesto por Cornell en 1969 y h representa el nivel fijo en el embalse considerado durante los análisis dinámicos. Finalmente, la función de probabilidad de falla pF(η|h) se presenta también en la figura 21, calculada como pF(η|h)=Φ(-β) , donde Φ(∙) se refiere a la función de distribución de probabilidad normal estándar.

Figura 21 Funciones de confiabilidad sísmica (izquierda) y curvas de fragilidad (derecha) ante el agrietamiento de la cortina 

Probabilidad de falla por deslizamiento en la base de la cortina. Modelo MEF2D

Se evaluó el deslizamiento en la base de la cortina ( Slb_MEF2D) mediante el valor máximo de la historia en el tiempo t de los desplazamientos relativos entre los puntos B y C mostrados en la figura 14. El deslizamiento puede entonces expresarse matemáticamente conforme a la ecuación 16.

La estimación de la probabilidad de falla considera que el deslizamiento ocurre exclusivamente en la dirección aguas abajo (ver figura 17). Esta suposición tiene un sentido físico al considerar que las solicitaciones en la cortina que producen el deslizamiento actúan mayormente en la dirección aguas abajo. De acuerdo con la función de estado límite presentada en la ecuación 9, la falla se alcanza cuando el deslizamiento Slb_MEF2D es mayor que 150 mm. Sin embargo, en ninguno de los 176 análisis dinámicos desarrollados se alcanzó esa condición. Por lo anterior, se procedió de la manera siguiente:

Figura 22 Resultados del ajuste de los valores de η en función de los desplazamientos relativos máximos entre la base de la cortina y la cimentación (deslizamientos) obtenidos con el modelo MEF2D 

La función lognormal de probabilidad obtenida en el paso 5 proporciona los parámetros EZF y σZF necesarios para obtener el índice de confiabilidad β conforme a la ecuación 8 (ver tabla 3). Tanto la función de confiabilidad como la curva de fragilidad ante el deslizamiento se muestran en la figura 23.

Figura 23 Funciones de confiabilidad sísmica (izquierda) y curvas de fragilidad (derecha) ante el deslizamiento en la base de la cortina utilizando el modelo MEF2D 

Probabilidad de falla por deslizamiento en la base de la cortina. Modelo MCR

Si se considera que la cortina de la presa es un cuerpo rígido y existe variación espacial de la cohesión y el ángulo de fricción en el contacto de la base con la cimentación mediante k discretizaciones, entonces el factor de seguridad ante el deslizamiento (FSS) puede calcularse como:

donde c _i es la cohesión, A _i es el área de la discretización, N _i es la fuerza normal actuando en la discretización, ϕ _i es el ángulo de fricción de la discretización y el término ∑Fh indica la sumatoria de todas las fuerzas horizontales actuando en la cortina: presión hidrostática, presión de azolves, fuerza inercial horizontal, presión hidrodinámica y presión dinámica de azolves.

Tomando la ecuación 17 en un contexto probabilista, tanto c _i como μ _i han sido considerados aleatorios con funciones de probabilidad lognormales definidas por los parámetros mostrados en la tabla 1. El término A _i es determinista pues depende del tamaño asignado a las discretizaciones; en un análisis bidimensional se considera que la sección de la cortina tiene un espesor unitario, por lo que Ai=Li , siendo Li la longitud de cada discretización. La fuerza normal resultante depende de tres factores principales: el peso propio de la cortina, el empuje vertical provocado por la subpresión y las fuerzas inerciales verticales. De estos tres factores tanto el peso propio W como las fuerzas inerciales verticales Qv pueden ser consideradas de forma probabilista. Por un lado, W=ADγc donde AD es el área total de la sección transversal de la cortina y yc es la densidad del concreto que se caracteriza mediante una función lognormal de probabilidad conforme a la tabla 1. Por su parte, Qv=Wkv donde kv es el coeficiente sísmico vertical; la caracterización probabilista de W ya fue señalada anteriormente. La subpresión depende del nivel del agua y de la eficiencia del sistema de drenaje de la cortina, cuyos valores se han considerado como deterministas. También, las acciones que determinan ∑Fh son deterministas, pues dependen exclusivamente del nivel del agua y del nivel de azolves, excepto las fuerzas inerciales horizontales Qh que se definen probabilísticamente de forma similar a Qv .

Para obtener las funciones de confiabilidad en este caso se utilizó el método de Montecarlo generando un muestreo aleatorio crudo (MAC) mediante la obtención de 10,000 valores aleatorios de ci , μi y γc con los parámetros mostrados en la tabla 1. La base de la cortina fue dividida en 10 discretizaciones, de las que 6 incluyeron cohesión y fricción (C+F) y 4 únicamente fricción (F), como se describe en la figura 6. Por ello, fue necesario generar 6 grupos de valores aleatorios de ci y 10 grupos de μi . Para las 10,000 simulaciones se desarrolló un análisis incremental del coeficiente sísmico horizontal kh , hasta que el factor de seguridad FSS obtuvo un valor de 0.999, implicando la falla por deslizamiento. Conviene recordar que kh=a0/(1+2a0) , actuando en la dirección aguas abajo, donde a0 está definida como la aceleración máxima en roca (PGA) y kv=2kh/3 es el coeficiente sísmico vertical; de esta forma se obtuvieron 10,000 coeficientes sísmicos horizontales de falla denotados como khF .

Se definió η como una medida de la demanda estructural, calculada mediante la relación entre cada uno de los coeficientes sísmicos horizontales de falla khF y el coeficiente sísmico asociado a la falla utilizando los valores medios de las propiedades mecánicas y de resistencia de la cortina presentados en la tabla 1. Los valores de η fueron ajustados a función lognormal de probabilidad mediante el método de máxima verosimilitud, el cual condujo a los parámetros Log[0.856,0.371], que corresponden a la función de densidad de probabilidades mostrada en la figura 24.

Figura 24 Función de densidad de probabilidad de η 

Los parámetros necesarios para aplicar la ecuación 8 y obtener las funciones de confiabilidad se presentan en la tabla 4. Se ha incluido un escenario adicional donde se ignora la cohesión en la zona de contacto de la base de la cortina con la cimentación pues algunas recomendaciones de diseño especifican una revisión adicional bajo estas condiciones (^{MOC, 2015}). Tanto la función de confiabilidad como la curva de fragilidad por el deslizamiento de la cortina en su base se presentan en la figura 25.

Figura 25 Funciones de confiabilidad (izquierda) y curvas de fragilidad sísmica (derecha) ante el deslizamiento en la base de la cortina utilizando el modelo MCR