Estabilidad de una presa de gravedad

En este artículo vamos a describir cómo optimizar con MATLAB el talud de aguas abajo de una presa de gravedad triangular, con talud de aguas arriba vertical, para que tenga un área mínima, manteniendo la estabilidad respecto al plano presa-cimiento y con unas restricciones respecto a la tensión de compresión mínima en el pie de aguas arriba, y a la máxima en el de aguas abajo. El caso programado se corresponde con la hipótesis accidental de drenes ineficaces, despreciando la cohesión en el plano de deslizamiento. Siguiendo las notaciones y dibujos de las transparencias que envió Rafael Morán y la metodología general expuesta en el artículo Guía de optimización en ingeniería el código sería el siguiente:

% Este programa calcula el talud óptimo de aguas abajo de una presa de gravedad 
% manteniendo la estabilidad y las condiciones de tensiones en cimentación 
% en la hipótesis de tensión plana. Supondremos que el perfil es un triángulo del cual 
% sólo podemos elegir el talud de aguas abajo que llamaremos n. El talud de 
% aguas arriba supondremos que es vertical. Buscaremos el perfil con área mínima.

% Parámetros específicos del problema
NMN       =170;  % Nivel máximo normal de embalse (msnm)
NAP       =174;  % Cota del vértice resistente (se considera en este caso igual al nivel de avenida de proyecto) (msnm)
Ter_arriba=108;  % Cota de la base de cimentación, que en este caso es horizontal (msnm)
Ter_abajo =118;  % Cota del terreno natural aguas abajo (se considera que el nivel del río aguas abajo está a esa cota) (msnm)
mu_h      =2.4;  % peso específico del hormigón en t/m3
mu_a      =1.0;  % peso específico del agua en t/m3
Phi       =pi/4; % ángulo de rozamiento en el contacto presa-cimiento (radianes)
F_Phi     =1.2;  % coeficiente de seguridad al rozamiento exigido por la normativa en la hipótesis de drenes ineficaces
S_Amin    =0.5;  % compresión mínima en el pie de aguas arriba (kp/cm2)
S_Bmax    =8.4;  % compresión máxima en el pie de aguas abajo (kp/cm2)     

% Variables de optimización: n = pendiente aguas abajo

% variables y funciones útiles para el cálculo (las acciones se evalúan en toneladas-fuerza por metro lineal de sección)
h_NMN = NMN-Ter_arriba;                     % altura de agua a embalse lleno sobre la base de cimentación (m)
H_a   = NAP-Ter_arriba;                     % altura de la presa (m)
h_CE  = Ter_abajo-Ter_arriba;               % subpresión en el pie de aguas abajo (mca)
base  = @(n) H_a*n;                         % ancho de la base del triángulo (m)
P     = @(n) 0.5*base(n)*H_a*mu_h;          % peso de la sección de hormigón (t/m)
E     = 0.5*h_NMN^2*mu_a;                   % empuje hidrostático con embalse lleno (no depende de n)
S     = @(n) 0.5*(h_NMN+h_CE)*base(n)*mu_a; % empuje debido a la subpresión en cimentación (t/m)

% Función coste: área del triángulo
f_coste = @(n) 0.5*base(n)*H_a;

% Calculamos ahora las restricciones:
% %% 1. Estabilidad frente al deslizamiento: E-T<0
% donde T es la fuerza necesaria para vencer el rozamiento:
% T = (P-S)*tan(Phi)/F_Phi (despreciando la cohesión y añadiendo el coef. seg. al rozamiento)
g_rest1 = @(n) E-(P(n)-S(n))*tan(Phi)/F_Phi; 
% la estabilidad requiere g_rest1 < 0

% %% 2. Compresión en el pie de aguas abajo menor que el máximo admitido por el cimiento
% N=P-S=1/2*(S_A+S_B)*base(n)
% Suponiendo la compresión mínima S_Amin dada,
% S_B=(P-S)/(base(n)/2)-S_Amin que en este caso no depende de n pero puede
% depender en general, así que lo definimos como función
S_B = @(n) (P(n)-S(n))/(0.5*base(n))-S_Amin*10;
g_rest2=@(n) S_B(n)-S_Bmax*10;
% la estabilidad requiere g_rest2 < 0

% %% 3. Equilibrio de momentos: el momento resultante debe ser tal que 
% el pie de aguas arriba se mantenga comprimido, con un valor superior a S_Amin. 
% Momento= Momento del peso 
% + momento empuje agua 
% + momento subpresion S_2
% + momento subpresion S_1
% + momento de la normal N_2
% + momento de la normal N_1
g_rest3=@(n) -P(n)*base(n)*2/3 ...
    + E*h_NMN/3 ...
    + h_CE*base(n)*(base(n)/2)*mu_a ...
    + 0.5*base(n)*(h_NMN-h_CE)*2/3*base(n)...
    + S_Amin*base(n)*(base(n)*0.5) ...
    + 0.5*base(n)*(S_B(n)-S_Amin)*(base(n)/3);
% la estabilidad requiere g_rest3 < 0

% %% Planteamiento matemático del problema de optimización:
% Buscamos n en el intervalo [0.1,80] que minimice la función coste:
% f_coste(n)
% con las restricciones:
% g_rest1(n)<0, g_rest2(n)<0, g_rest3(n)<0.

%% Optimización por fuerza bruta
% Intervalo de parámetros:
h  =0.01; % distancia entre dos valores del parámetro
L  =0.1;  % extremo izquierdo
M  =80;   % extremo derecho 
n  =L:h:M;% vector con los parámetros
% Inicialización
Optimo=10000; % Inicializamos con un valor muy alto 
indice=0;     % inicializamos el indice
% Probamos todos los valores n(i)
for i=1:length(n)
    % creamos una variable lógica que sea cierta sólo si se verifican todas las restricciones
    restriccion=(g_rest1(n(i))<0)&(g_rest2(n(i))<0)&(g_rest3(n(i))<0);
    if restriccion              % solo en este caso n(i) es un parámetro admisible
        coste = f_coste(n(i));  % Calculo el coste
        if coste<Optimo         % en este caso mejoramos!
            Optimo = coste;     % actualizo el coste con el nuevo parámetro
            indice=i;           % actualizo el índice en el que está el óptimo
        end
    end        
end
% Escribimos el resultado de la optimización por pantalla
if indice == 0
    sprintf('Ningún valor de n cumple las restricciones')
else
    sprintf('El óptimo se alcanza en n = %f',n(indice))
end