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.

% 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
% 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 cálculo
h_NMN=NMN-Ter_arriba;
H_a=NAP-Ter_arriba;
h_CE=Ter_abajo-Ter_arriba;

% Variables de optimización
% n ángulo aguas abajo
% funciones útiles
base= @(n) H_a*n; % base del triángulo

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

% Calculamos ahora las restricciones:
% Calculamos la función peso
P = @(n) 0.5*base(n)*H_a*mu_h;

% Calculamos el empuje (no depende de n)
E = 0.5*mu_a*h_NMN^2;

% Calculamos la subpresión
S= @(n) 0.5*(h_NMN+h_CE)*base(n)*mu_a;

% Estabilidad frente al deslizamiento: E-T<0
% 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 f_rest1 < 0

% Equilibrio de fuerzas normales P-S>0
g_rest2 = @(n) S(n)-P(n);
% la estabilidad requiere f_rest2 < 0

% 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 f_rest3 < 0

% Equilibrio de momentos
% 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 f_rest4 < 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=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))