Estabilidad de una presa de gravedad

En este artículo vamos a describir cómo optimizar con MATLAB el ángulo aguas abajo de una presa de gravedad triangular. Siguiendo las notaciones y dibujos de las transparencias que envió Rafa el código sería el siguiente:

% Cálculo de estabilidad de presas de gravedad
% Supondremos que el perfil es un triángulo del cual sólo podemos elegir la
% pendiente aguas abajo que llamaremos n
% 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. Equilibrio de fuerzas normales 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_A dada,
% S_B=(P-S)/(H_a*n/2)-S_A que en este caso no depende de n pero puede
% depender en general
S_B = @(n) (P(n)-S(n))/(0.5*H_a*n)-S_Amin*10;
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))