Estabilidad de una presa de gravedad II - Historial de revisiones

Carlos Castro en 10:38 14 dic 2017

2017-12-14T10:38:05Z

Carlos Castro en 11:54 16 mar 2017

2017-03-16T11:54:50Z

Carlos Castro: /* Programa en MATLAB */

2017-03-10T21:05:27Z

‎Programa en MATLAB

Deprecated: The each() function is deprecated. This message will be suppressed on further calls in /home/mat/public_html/w/includes/diff/DairikiDiff.php on line 434

Carlos Castro: /* Programa en MATLAB */

2017-03-10T18:02:20Z

‎Programa en MATLAB

Carlos Castro: /* Programa en MATLAB */

2017-03-10T18:00:21Z

‎Programa en MATLAB

Carlos Castro en 17:59 10 mar 2017

2017-03-10T17:59:24Z

Carlos Castro: Página creada con «En este artículo vamos a describir cómo optimizar con MATLAB el talud de aguas abajo de una presa de gravedad triangular, con talud de aguas arriba vertical, para que ten...»

2017-03-10T15:07:18Z

Página creada con «En este artículo vamos a describir cómo optimizar con MATLAB el talud de aguas abajo de una presa de gravedad triangular, con talud de aguas arriba vertical, para que ten...»

Página nueva

En este artículo vamos a describir cómo optimizar con MATLAB el talud de aguas abajo de una presa de gravedad triangular, con talud de aguas arriba vertical, para que tenga un área mínima, manteniendo la estabilidad respecto al plano presa-cimiento y con unas restricciones respecto a la tensión de compresión mínima en el pie de aguas arriba, y a la máxima en el de aguas abajo. El caso programado se corresponde con la hipótesis accidental de drenes ineficaces, despreciando la cohesión en el plano de deslizamiento. Siguiendo las notaciones y dibujos de las transparencias que envió Rafael Morán y la metodología general expuesta en el artículo [[Guía de optimización en ingeniería|Guía de optimización en ingeniería]] el código sería el siguiente:

== Programa en MATLAB ==

{{matlab|codigo=
% 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)
N = @(n) P(n)-S(n); % Normal, obtenida por el equilibrio de fuerzas verticales
% Necesitaremos también las tensiones en los pies de aguas arriba y aguas abajo: S_A y S_B
% Para ello usaremos lo siguiente:
% 1. Supondremos que la distribución de las tensiones asociadas a la
% normal crece linealmente desde el pie de aguas arriba hacia el de
% aguas abajo.
% 2. El equilibrio de momentos en el pie de aguas abajo
% Estas dos condiciones se escriben:
% Crec. lineal de ten. nor. N = 1/2*(S_A+S_B)*base(n)
% 0 = Momento del peso 0 = - P(n)*base(n)*2/3
% + momento empuje agua + E*h_NMN/3
% + momento subpresion S_2 + h_CE*base(n)*(base(n)/2)*mu_a
% + momento subpresion S_1 + 0.5*base(n)*(h_NMN-h_CE)*2/3*base(n)
% + momento de la normal N_2 + S_A*10*base(n)*(base(n)*0.5)
% + momento de la normal N_1 + 0.5*base(n)*(S_B-S_A)*(base(n)/3)
% Tenemos 2 ecuaciones con dos incógnitas. Despejando sacamos el valor de
% S_A y S_B
S_A = @(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)...
- 0.5*base(n)*(2*N(n)/base(n))*base(n)/3)/(base(n)^2*(1/2-1/3));
S_B = @(n) 2*N(n)/base(n)-S_A(n);
% Una vez definidas todas las funciones que necesitamos, calculamos la
% función coste y las restricciones

% 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.
g_rest2 = @(n) 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.
g_rest3 = @(n) S_Amin*10-S_A(n);
% la estabilidad requiere g_rest3(n) < 0

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

== Observaciones: ==

A continuación se muestra como hacer algunas de las cuestiones que surgieron en la última reunión.

* '''Funciones que dependen de varias variables.''' Es posible que tanto la función coste como las restricciones dependan de más de una variable. Por ejemplo, supongamos que nuestros parámetros de optimización son las coordenadas (x,y) de un punto y la función coste es la distancia a un punto (a,b). Entonces podemos escribir
{{matlab|codigo=
a=2; b=1;
f_coste = @(x,y) sqrt((x-a)^2+(y-b)^2);
}}
* '''Funciones que dependen de variables matriciales.''' Si los parámetros de optimización son coordenadas de varios puntos es mejor escribir esas coordenadas en una matriz que tenga por columnas las coordenadas de cada punto y usar una única variable para crear la función. Por ejemplo,
{{matlab|codigo=
% 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.
}}
Observamos que en este caso la función no devuelve un número sino un vector columna que contiene las coordenadas del centro de masas de los puntos.
* '''Definir funciones según un parámetro.''' En algún momento puede que el cálculo dependa de algún parámetro y requiera funciones diferentes según el valor del parámetro. En este caso se puede usar una condición para definir una función de una manera u otra. Ejemplo:
{{matlab|codigo=
% 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
}}
De todas formas conviene evitar este tipo situaciones ya que los métodos de optimización para problemas de esta forma son mucho menos eficientes.
* '''Restricciones adicionales.''' Además de las restricciones impuestas por la estabilidad, las geometrías que se traten pueden requerir más restricciones. Por ejemplo, si nuestras variables de optimización son varios puntos <math> (x_i, y_i), i=1,...,n, </math> y queremos que las coordenadas <math> y_i </math> de estos puntos estén ordenadas, es decir <math> y_{i+1}>y_i </math> para <math> i=1,..,n-1</math>, esto puede imponerse añadiendo una restricción más que podríamos escribir de esta forma
{{matlab|codigo=
% 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
}}

* '''Funciones complejas.''' Algunas veces las funciones que necesitamos no pueden escribirse en una sola línea fácilmente. En este caso, podemos escribirlas en un fichero aparte como en el ejemplo que mostramos para calcular el área de un polígono a partir de las coordenadas de los vértices con la fórmula que envío Rafael Morán.

{{matlab|codigo=
% 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);
}}

Es importante grabar esta función en un fichero con el mismo nombre elegido para la función. En este caso: ar_pol.m. Una vez creada la función y grabada en el mismo directorio donde está nuestro programa principal, la podemos usar en nuestro programa. Por ejemplo

{{matlab|codigo=
% 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)
}}
[[Categoría:Trabajos fin de Máster]]

@@ Línea 35: / Línea 35: @@
 % Necesitaremos también las tensiones en los pies de aguas arriba y aguas abajo: S_A y S_B
 % Para ello usaremos lo siguiente:
-%    1. Supondremos que la distribución de las tensiones asociadas a la
+%    1. Supondremos que la distribución de las tensiones verticales en la base
-%    normal varía linealmente desde el pie de aguas arriba (donde toma el valor S_A) hacia el de
+%    varía linealmente desde el pie de aguas arriba (donde toma el valor S_A) hacia el de
 %    aguas abajo (donde toma el valor S_B).
 %    2. El equilibrio de momentos en el pie de aguas abajo

@@ Línea 143: / Línea 143: @@
 }}
 De todas formas conviene evitar este tipo situaciones ya que los métodos de optimización para problemas de esta forma son mucho menos eficientes.
-* '''Restricciones adicionales.''' Además de las restricciones impuestas por la estabilidad, las geometrías que se traten pueden requerir más restricciones. Por ejemplo, si nuestras variables de optimización son varios puntos <math>  (x_i, y_i), i=1,...,n, </math>  y queremos que las coordenadas <math> y_i </math> de estos puntos estén ordenadas, es decir <math> y_{i+1}>y_i </math> para <math> i=1,..,n-1</math>, esto puede imponerse añadiendo una restricción más que podríamos escribir de esta forma
+* '''Restricciones adicionales.''' Además de las restricciones impuestas por la estabilidad, las geometrías que se traten pueden requerir más restricciones. Por ejemplo, si nuestras variables de optimización son varios puntos <math>  (x_i, y_i), i=1,...,n, </math>  y queremos que las coordenadas <math> y_i </math> de estos puntos estén ordenadas, es decir <math> y_{i+1} \gt y_i </math> para <math> i=1,..,n-1</math>, esto puede imponerse añadiendo una restricción más que podríamos escribir de esta forma
 {{matlab|codigo=
 % Llamamos Cxy la matriz que contiene por columnas las coordenadas de varios puntos y n el número de puntos, que

@@ Línea 19: / Línea 19: @@
 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_Amin    = 5;  % compresión mínima en el pie de aguas arriba (t/m2)
-S_Bmax    = 8.4;  % compresión máxima en el pie de aguas abajo (kp/cm2)
+S_Bmax    = 84;  % compresión máxima en el pie de aguas abajo (t/m2)
 % Variables de optimización: n = pendiente aguas abajo
@@ Línea 69: / Línea 69: @@
 % %% 2. Compresión en el pie de aguas abajo menor que el máximo admitido
 % por el cimiento.
-g_rest2 = @(n) S_B(n)-S_Bmax*10;
+g_rest2 = @(n) S_B(n)-S_Bmax;
 % la estabilidad requiere g_rest2(n) < 0
 % %% 3. El pie de aguas arriba debe mantenerse comprimido, con un valor superior a S_Amin.
-g_rest3 = @(n) S_Amin*10-S_A(n);
+g_rest3 = @(n) S_Amin-S_A(n);
 % la estabilidad requiere g_rest3(n) < 0

@@ Línea 63: / Línea 63: @@
 % %% 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)
+% T = N*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;
+g_rest1 = @(n) E-N(n)*tan(Phi)/F_Phi;
 % la estabilidad requiere g_rest1(n) < 0

@@ Línea 36: / Línea 36: @@
 % Para ello usaremos lo siguiente:
 %    1. Supondremos que la distribución de las tensiones asociadas a la
-%    normal crece linealmente desde el pie de aguas arriba hacia el de
+%    normal varía linealmente desde el pie de aguas arriba hacia el de
 %    aguas abajo.
 %    2. El equilibrio de momentos en el pie de aguas abajo