Estabilidad de una presa de gravedad

En este artículo vamos a describir cómo optimizar con MATLAB la pendiente aguas abajo de una presa de gravedad triangular, con pendiente aguas arriba vertical, para que tenga área mínima manteniendo la estabilidad. Siguiendo las notaciones y dibujos de las transparencias que envió Rafa el código sería el siguiente:

% Este programa calcula la pendiente óptima aguas abajo de una presa de gravedad 
% manteniendo la estabilidad. Supondremos que el perfil es un triángulo del cual 
% sólo podemos elegir la pendiente aguas abajo que llamaremos n. La pendiente 
% aguas arriba supondremos que es vertical. Buscaremos el perfil con área mínima.

% Parámetros
NMN       =170; % Nivel máximo normal (m)
NAP       =174; % Nivel de avenida m
Ter_arriba=108; % unidad m
Ter_abajo =118; % unidad m
mu_h      =2.4; % densidad del hormigón en t/m3
mu_a      =1.0; % densidad del agua en t/m3
Phi       =pi/4; % ángulo de rozamiento
F_Phi     =1.2; % coeficiente seguridad respecto al rozamiento
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
h_NMN = NMN-Ter_arriba;                     % altura normal máxima
H_a   = NAP-Ter_arriba;                     % altura
h_CE  = Ter_abajo-Ter_arriba;               % altura del terreno aguas abajo
base  = @(n) H_a*n;                         % base del triángulo
P     = @(n) 0.5*base(n)*H_a*mu_h;          % Peso total
E     = 0.5*h_NMN^2*mu_a;                   % Empuje (no depende de n)
S     = @(n) 0.5*(h_NMN+h_CE)*base(n)*mu_a; % subpresión

% 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. Estabilidad de fuerzas verticales S-P<0
g_rest2 = @(n) S(n)-P(n);
% la estabilidad requiere g_rest2 < 0

% %% 3. Tensión en el pie de aguas abajo menor que el máximo permitido
% N=P-S=1/2*H_a*n*(S_A+S_B)
% Suponiendo la compresión mínima S_Amin dada,
% S_B=(P-S)/(H_a*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*H_a*n);
g_rest3=@(n) S_B(n)-S_Bmax*10;
% la estabilidad requiere g_rest3 < 0

% %% 4. Equilibrio de momentos: el momento en sentido antihorario debe ser suficiente para que 
% la tensión en el pie de aguas arriba sea superior a S_Amin. De forma equivalente planteamos el 
% equilibrio haciendo que cuando la tensión en el pie de aguas arriba es S_Amin, el momento
% resultante vaya en sentido antihorario
% 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_rest4=@(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_rest4 < 0

% Pasamos ahora a optimizar la pendiente n de manera que se cumplan las restricciones
%% 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  =20;   % 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)&(g_rest4(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
sprintf('El óptimo se alcanza en n = %f',n(indice))