Diferencia entre revisiones de «Estabilidad de una presa de gravedad»

% 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 (en este caso igual al Nivel de avenida de proyecto (msnm)
Ter_arriba=108;  % Nivel terreno aguas arriba (m)
Ter_abajo =118;  % Nivel terreno aguas abajo (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 con embalse lleno (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*(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));
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))