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 (se considera en este caso igual al nivel de avenida de proyecto) (msnm)
Cota_cim  = 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-Cota_cim;                       % altura de agua a embalse lleno sobre la base de cimentación (m)
H_a   = NAP-Cota_cim;                       % altura de la presa (m)
h_CE  = Ter_abajo-Cota_cim;                 % 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(n) < 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 que al menos tenemos la compresión mínima en el pie de aguas arriba, es decir S_A >= S_Amin,
% podemos calcular una cota superior para la tensión aguas abajo
% S_B<=(P-S)/(base(n)/2)-S_Amin =Cota_S_B que en este caso no depende de n pero puede
% depender en general, así que lo definimos como función
Cota_S_B = @(n) (P(n)-S(n))/(0.5*base(n))-S_Amin*10;
g_rest2 = @(n) Cota_S_B(n)-S_Bmax*10;
% la estabilidad requiere g_rest2(n) < 0

% %% 3. El pie de aguas arriba debe mantenerse comprimido, con un valor superior a S_Amin.
% El cálculo del momento resultante permite establecer esta restricción.   
% Momento= Momento del peso 
% + momento empuje agua 
% + momento subpresion S_2
% + momento subpresion S_1
% + momento de la normal N_2 (cuando suponemos S_A=S_Amin)
% + momento de la normal N_1 (cuando suponemos S_A=S_Amin y por tanto S_B=Cota_S_B)
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*10*base(n)*(base(n)*0.5) ...
    + 0.5*base(n)*(Cota_S_B(n)-S_Amin*10)*(base(n)/3);
% la estabilidad requiere g_rest3(n) < 0
% De hecho, si S_A<S_Amin, y asumimos que la distribución de tensiones crece linealmente cuando 
% nos movemos hacia el pie de aguas abajo, el momento de las fuerzas normales (para un mismo peso y subpresión) 
% sería menor que el del caso S_A=S_Amin y la restricción no se verificaría 

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

a=2; b=1; 
f_coste = @(x,y) sqrt((x-a)^2+(y-b)^2);

% Llamamos Cxy la matriz que contiene por columnas las coordenadas de varios puntos y n el número de puntos, que 
% coincidirá con el número de columnas de la matriz Cxy. La siguiente función
% devuelve el centro de masas (suponiendo que cada punto tiene masa unidad) 
c_masas = @(Cxy,n) [sum(Cxy(1,:))/n;sum(Cxy(2,:))/n];
% El argumento : se usa para seleccionar todos los posibles índices de un vector o matriz. 
% Por ejemplo, Cxy(1,:) se refiere a la primera fila de la matriz Cxy, ya que tomo la fila 1 y todas las posibles columnas.

% Vamos a definir una función dependiendo de un parámetro aleatorio 
% Elijo un número aleatorio entre 0 y 1
a = rand(1);
% Dependiendo de si a>1/2 o a<=1/2 defino la función coste
if a<=0.5
    f_coste = @(n) 2*n+4;
else
    f_coste = @(n) 2*n-4;
end

% Llamamos Cxy la matriz que contiene por columnas las coordenadas de varios puntos y n el número de puntos, que 
% coincidirá con el número de columnas de la matriz Cxy. La siguiente restricción hace que 
% las coordenadas y de los puntos estén ordenadas 
g_rest4 = @(Cxy,n) -min(Cxy(2,2:n)-Cxy(2,1:n-1));
% las coordenadas y de los puntos están ordenadas si g_rest4(Cxy,n)<0

% Función para calcular el área de un polígono
function S=ar_pol(C)
% Se introduce una matriz C con las coordenadas de los vértices del polígono 
% por columnas y ordenados en sentido antihorario. 
% Se devuelve el área S
% Aplicamos la fórmula que nos envió Rafael Morán
% S = 1/2 abs(sum yi (x_(i-1)-x_(i+1)))
% teniendo en cuenta la siguiente regla: (x_0,y_0)=(x_n,y_n) y 
% (x_(n+1),y_(n+1))=(x_1,y_1)
% donde n es el número de vértices

% Primero calculamos el número de vértices nv
nv=length(C);
% Ahora definimos los vectores con las coordenadas x e y de los vértices
x=C(1,:);    % la primera fila contiene las coordenadas x de los vértices
y=C(2,:);    % la segunda fila contiene las coordenadas y de los vértices
% Calculamos S pero separamos el primer y último término de la suma
% para poder tener en cuenta la regla
S = y(1)*(x(nv)-x(2));           % primer término
for i=2:nv-1
   S = S+y(i)*(x(i-1)-x(i+1));   % términos del 2 al nv-1
end
S = S+y(nv)*(x(nv-1)-x(1));      % último término
S=0.5*abs(S);

% Este es el programa principal
% Definimos una matriz con las coordenadas del polígono
c=[-1 0 -1 1 1;-1 0 1 1 -1];
% Calculamos su área
S=ar_pol(c)