La presa de El Atazar. Grupo 16

De MateWiki
Revisión del 14:12 8 dic 2024 de Juan Felix Aguilar (Discusión | contribuciones) (→‎.- Cálculo de la fuerza de la presión total y la presión por unidad de superficie.)

Saltar a: navegación, buscar


Trabajo realizado por estudiantes
Título La presa de El Atazar. Grupo 16
Asignatura Teoría de Campos
Curso 2024-25
Autores
  • Jaime Moral Riquelme
  • Dionisio José García Álvarez
  • Juan Felix Aguilar Romero
  • Alejandro Santisteban Sancho
  • Jose Pablo Bonilla Hurtado
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura
Presa de El Atazar (Madrid). Fuente: Colegio de Ingenieros de Caminos, Canales y Puertos de Madrid.

La presa de El Atazar, construida entre 1968 y 1972, es la más grande de la Comunidad de Madrid y una estructura clave para el suministro de agua de la capital y de toda la región. Esta presa de doble curvatura forma un arco tanto en la vista superior como en la sección vertical, alcanzando una altura de 134 metros y una longitud de coronación de 484 metros. Consideraremos la superficie de la presa en el lado aguas arriba, que está en contacto con el agua. Suponemos que la sección transversal de la presa sea un arco de circunferencia con un eje de simetría ubicado en el valle, mientras que la sección longitudinal se comporta como un arco parabólico. En coordenadas cilíndricas (r,θ,z), la superficie puede modelarse mediante las siguientes ecuaciones:

[math]r = r_{0} + b (1 - \frac{z^2}{h^2})[/math],[math]\qquad θ ∈ [\frac{3π}{4}, \frac{5π}{4}][/math], [math]\qquad Z ∈ [0,H][/math]

Siendo:

  • H: altura de la presa,
  • [math]r_{0}[/math]: radio de la presa en la altura máxima,
  • b = 35 m: factor que determina la curvatura del arco parabólico.

Consideramos el campo escalar de la presión que ejerce el agua como:

  • [math]P(z)=ρgh(z)[/math]

Además el campo vectorial de la fuerza de presión viene dado por:

  • [math]\overrightarrow{F} = −P(z) \overrightarrow{n}[/math]


1 .- Representación de la cara interior de la presa.

En primer lugar, vamos a representar la superficie aguas arriba de la presa. Esta se representa directamente con la ecuación de la superficie mencionada anteriormente. El color que hemos elegido para la representación de este primer gráfico es azul, dado que la superficie representaría la pared contra la cual el agua ejercería la presión.

Figura 1. Representación de la cara interior de la presa, vista aguas arriba.
clear; clc;

% 1. Representar en MATLAB la superficie parametrizada de la presa en su
% cara de aguas arriba, en un único color

% Asignación de variables
H = 134;                              %Altura de la presa (metros)
r0 = 242;                             %Radio de la presa en la altura máxima (metros)
b = 35;                               %Curvatura del arco parabólico (metros)
theta = linspace(3*pi/4,5*pi/4,100);  %Angulo theta coordenadas cilíndricas
z = linspace(0,H,100);                %Altura z coordenadas cilíndricas
[Mthetha,Mz] = meshgrid(theta,z);     %Mallado
r = r0+b*(1-(Mz.^2/H^2));             %Ecuación de la superficie de la presa

% Conversión de cilíndricas a cartesianas
X = r.*cos(Mthetha);
Y = r.*sin(Mthetha);

% Creación del gráfico
figure;
surf(X,Y,Mz,'FaceColor','#4DBEEE','EdgeColor','none');
title('Presa de El Atazar aguas arriba');
xlabel('X (m)');
ylabel('Y (m)');
zlabel('Z (m)');
axis equal;
grid on;


Como se puede observar en la figura 1, tenemos una superficie curva con un radio de curvatura inferior de mayor dimensión que el radio de curvatura de mayor altura de la pared de la presa. Para representar la superficie de un único color, hemos hecho uso de la función 'surf' de MATLAB.

2 .- Representación del campo escalar de presiones.

Utilizando la fórmula fundamental de la estática de fluidos, calcularemos las presiones que sufre la cara interior de la presa. Como es lógico, la presión aumenta de manera proporcional según aumenta la profundidad (color más rojo).

Figura 2. Representación de las presiones en la cara interior de la presa.
clear; clc;

% 2. Representar el campo escalar de presión como un mapa de colores sobre
% la superficie parametrizada de la presa (usa también colorbar)

% Asignación de parámetros
rho = 1000;     %Densidad del agua (Kg/(m^3))
g = 9.81;       %Gravedad terrestre

% Relación presión - altura
P = rho*g*(H-Mz);

% Creación del gráfico
figure;
surf(X,Y,Mz,P,"EdgeColor","none");
title('Campo escalar de presion en la presa');
xlabel('X (m)');
ylabel('Y (m)');
zlabel('Z (m)');
axis equal;
grid on;
colormap("jet");
colorbar;

Para la representación de esta gráfica, hemos hecho uso de la función 'surf', al igual que en el apartado anterior, y además hemos usado la función 'colorbar' para representar dónde es más fuerte la presión.

3 .- Representación del campo escalar de la fuerza de presión.

En este apartado, vamos a trabajar con el campo vectorial de la presión que ejerce el agua contra la pared interior de la presa.

Figura 3. Representación del campo vectorial de presiones en la cara interior de la presa.
Figura 4. Representación del corte en planta del campo vectorial de presiones.
clear; clc;

% 3. Representación de campo de presiones 
% Derivadas y vectores normales
[Rx_theta, Rx_z] = gradient(X);
[Ry_theta, Ry_z] = gradient(Y);
[Rz_theta, Rz_z] = gradient(gridH);

% Producto cruzado para el vector normal
Nx = Ry_theta .* Rz_z - Rz_theta .* Ry_z;
Ny = Rz_theta .* Rx_z - Rx_theta .* Rz_z;
Nz = Rx_theta .* Ry_z - Ry_theta .* Rx_z;

% Magnitud del vector normal
N_magnitud = sqrt(Nx.^2 + Ny.^2 + Nz.^2);

% Vector normal unitario
Nx_unit = Nx ./ N_magnitud;
Ny_unit = Ny ./ N_magnitud;
Nz_unit = Nz ./ N_magnitud;

% Campo de fuerza de presión
Fx = -P .* Nx_unit;
Fy = -P .* Ny_unit;
Fz = -P .* Nz_unit;

% Representación del campo vectorial sobre la superficie
figure;
quiver3(X, Y, gridH, Fx, Fy, Fz, 0.5, 'Color', 'r');
hold on;
surf(X, Y, gridH, P, 'EdgeColor', 'none', 'FaceAlpha', 0.8);
colorbar;
xlabel('X (m)');
ylabel('Y (m)');
zlabel('Altura Z (m)');
title('Campo vectorial de la fuerza de presión sobre la presa');
view(-60, 20); % Cambia el ángulo de vista
axis equal;
grid on;

% Corte en planta
corte = 51; % Elegimos un corte en el plano longitudinal (aproximadamente theta = 0)
X_corte = X(:, corte);
Z_corte = gridH(:, corte);
P_corte = P(:, corte);

Nx_corte = Nx(:, corte);
Nz_corte = Nz(:, corte);
N_magnitud_corte = sqrt(Nx_corte.^2 + Nz_corte.^2);

Nx_unit_corte = Nx_corte ./ N_magnitud_corte;
Nz_unit_corte = Nz_corte ./ N_magnitud_corte;

Fx_corte = -P_corte .* Nx_unit_corte;
Fz_corte = -P_corte .* Nz_unit_corte;

% Representación del corte longitudinal
figure;
quiver(X_corte, Z_corte, Fx_corte, Fz_corte, 'r', 'LineWidth', 1.5);
hold on;
plot(X_corte, Z_corte, 'k-', 'LineWidth', 2);
xlabel('X (m)');
ylabel('Z (m)');
title('Campo de fuerza de presión en el corte longitudinal');
grid on;

En la figura 3, podemos ver la representación en tercera dimensión del campo vectorial de presiones. Hemos conseguido representarla haciendo uso de la función 'quiver3' la cual sirve para representar campos vectoriales tridimensionales, en el gráfico, este queda representado en rojo. Además hemos representado con la función 'surf' la superficie de la pared interior de la presa.
En la figura 4, tenemos la representación en planta de lo que sería el corte longitudinal de la presa en planta. Para ella hemos hecho uso de la función 'quiver' para representar el campo vectorial de presiones en 2D y la función 'plot' para representar la curva, es decir, la sección de la pared de la presa en planta.


4 .- Representación de la curva que describe la trayectoria de una gota de agua al salir de la compuerta.

Teniendo en cuenta que la altura de la compuerta de la presa es 25 metros, que el área de la compuerta son 1100 metros cuadrados y que la fórmula de Tonirreli es:

[math]v_{0}=\sqrt{2gH_{c}}[/math]


Podemos determinar la trayectoria que describiría una gota de agua que sale de la compuerta.

Figura 5. Trayectoria que tendría una gota que sale de la compuerta de la presa.
clear; clc;

% Parámetros que vamos a utilizar
theta = 15 * pi / 180; % Ángulo de salida (pasamos a radianes)
g = 9.81;              % Gravedad en metros partido de segundo cuadrado
Hc = 25;               % Altura de la compuerta con respecto al suelo

% Velocidad inicial, utilizamos la fórmula de Torricelli
v0 = sqrt(2 * g * Hc);

% Determinamos los factores posición
t = linspace(0, 5, 1000); % Tiempo
x = v0 * cos(theta) * t;  % Posición horizontal (X)
y = Hc + v0 * sin(theta) * t - 0.5 * g * t.^2; % Posición vertical (Y)

% Hallamos el tiempo para el cual la altura es 0, es decir, el momento en el que la gota impacta contra el suelo
t_impacto = fzero(@(t) Hc + v0 * sin(theta) * t - 0.5 * g * t^2, [0, 10]);
x_impacto = v0 * cos(theta) * t_impacto; % Distancia desde la compuerta hasta el punto de impacto con el suelo

% Representación de la gráfica de la trayectoria que describe la gota
figure;
plot(x, y, 'Color', '#4DBEEE', 'LineWidth', 1.5); % Ponemos el mismo color que hemos utilizado para la primera figura
hold on;
plot(x_impacto, 0, 'ro', 'MarkerSize', 8, 'DisplayName', 'Impacto con el suelo'); 
xlabel('Distancia horizontal de la gota (metros)');
ylabel('Altura de la gota (metros)');
title('Trayectoria de la gota de agua');
grid on;

% Insertamos una leyenda que describa lo que vemos
legend('Trayectoria gota', 'Impacto con el suelo'); 
xlim([0, x_impacto + 12]);
ylim([0, Hc + 7]);

El punto (25,0) representaría el punto en el que se encuentra la compuerta, es decir, el punto de partida de la gota, por otro lado, el punto que se encuentra entre x=60 y x=70 que está rodeado de un círculo rojo, sería el punto final de la trayectoria de la gota, es decir, el punto en el que la gota impacta contra el suelo.
La curva azul representada indica por dónde se movería la gota, es decir, su trayectoria.
Para la representación de la figura 5, al ser una función bidimensional, hemos hecho uso de la función 'plot' para la representación de la gráfica.

5 .- Representación del campo tangente y el campo normal sobre la curva anterior.

En este apartado, representamos los vectores tangentes y normales de la trayectoria que describe la gota anterior que sale de la compuerta.

Figura 6. Representación del campo tangente y el campo normal sobre la curva que define la trayectoria de la curva.
clear; clc;

% Parámetros que vamos a utilizar
theta = 15 * pi / 180;  % Pasamos el ángulo inicial a radianes
g = 9.81;               % Gravedad
Hc = 25;                % Altura compuerta
v0 = sqrt(2 * g * Hc);  % Velocidad inicial (con Tonirrelli)

% Trayectoria de la curva
t = linspace(0, 5, 500); % Tiempo
x = v0 * cos(theta) * t;   % Posición en el eje x
y = Hc + v0 * sin(theta) * t - 0.5 * g * t.^2; % Posición en el eje y

% Velocidades (derivadas de la posición)
vx = v0 * cos(theta) * ones(size(t)); % Velocidad en x (dx/dt)
vy = v0 * sin(theta) - g * t;         % Velocidad en y (dy/dt)
V = sqrt(vx.^2 + vy.^2);              % Magnitud de la velocidad

% Vector tangente
Tx = vx ./ V;  % Componente tangencial en x
Ty = vy ./ V;  % Componente tangencial en y

% Aceleraciones (derivadas de las velocidades)
ax = zeros(size(t));                % Aceleración en x (dvx/dt = 0)
ay = -g * ones(size(t));            % Aceleración en y (dvy/dt = -g)
a_magnitude = sqrt(ax.^2 + ay.^2);  % Magnitud de la aceleración

% Vector normal
Nx = (ax - Tx .* (ax .* Tx + ay .* Ty)) ./ a_magnitude; % Componente normal en x
Ny = (ay - Ty .* (ax .* Tx + ay .* Ty)) ./ a_magnitude; % Componente normal en y

% Representación de la gráfica de los campos tangente y normal sobre la curva en varios instantes
figure;
plot(x, y, 'Color', '#4DBEEE', 'LineWidth', 1.5); hold on; 

% Campos tangente y normal en puntos específicos
puntos = 15; % Número de vectores a graficar
indices = round(linspace(1, length(t), puntos));  % Selección de índices

% Vectores tangentes en rojo
quiver(x(indices), y(indices), Tx(indices), Ty(indices), 0.3, 'r', 'LineWidth', 1.5, 'DisplayName', 'Tangente');

% Vectores normales en color azul oscuro
quiver(x(indices), y(indices), Nx(indices), Ny(indices), 0.3, 'Color', [0, 0, 0.5], 'LineWidth', 1.5, 'DisplayName', 'Normal');

xlabel('Distancia horizontal (m)');
ylabel('Altura (m)');
title('Campos Tangente y Normal sobre la trayectoria');
legend('Trayectoria', 'Tangente', 'Normal');
grid on;
axis equal;

% Ajustar los límites de la gráfica para asegurar que se vea la intersección con el eje x (y=0) y el eje y (x=0)
xlim([min(0, min(x) - 5), max(x) + 5]);
ylim([min(0, min(y) - 5), max(y) + 5]);


Para la elaboración de la figura 6 hemos usado la función 'quiver' para representar los campos vectoriales en 2D (Tangencial y normal) y la función 'plot' para la representación de la representación de la curva.

6 .- Representación de la animación del vector velocidad y el vector aceleración de la gota.

Sobre la gráfica anterior, se representa el cambio de la velocidad y de la aceleración que sufre una gota al seguir la curva a lo largo del recorrido. 

clear; clc;

% Parámetros ajustados para la trayectoria
g = 9.81;                         % Gravedad (m/s^2)
v0 = 30;                          % Velocidad inicial (m/s)
theta = pi/3;                     % Ángulo de lanzamiento (60 grados)
Hc = 25;                          % Altura inicial (m)
t_max = 2 * v0 * sin(theta) / g + sqrt(2 * Hc / g);  % Tiempo total de vuelo
t = linspace(0, t_max, 200);      % Intervalo de tiempo para la animación

% Posiciones en x e y durante la animación
x = v0 * cos(theta) .* t;               
y = Hc + v0 * sin(theta) .* t - 0.5 * g .* t.^2;

% Ajustar límites para evitar valores negativos en y
valid_idx = y >= 0;               % Solo considerar valores donde y >= 0
x = x(valid_idx);
y = y(valid_idx);

% Crear y configurar la figura

figure;
hold on;
plot(x, y, 'c', 'LineWidth', 1.5, 'DisplayName', 'Trayectoria gota'); % Trayectoria
xlabel('Distancia horizontal de la gota (metros)');
ylabel('Altura de la gota (metros)');
title('Trayectoria de la gota de agua (animada)');
grid on;
axis equal;
xlim([25, max(x)]);
ylim([0, Hc + 75]);

% Inicializar marcador de la gota
gota = plot(x(1), y(1), 'ro', 'MarkerSize', 10, 'DisplayName', 'Gota en movimiento');

legend('Location', 'best');
hold off;

% Animar la trayectoria
for i = 1:length(x)
    % Actualizar la posición de la gota
    set(gota, 'XData', x(i), 'YData', y(i));
    
    % Pausa para crear el efecto de animación
    pause(0.05);
end

% Agregar marcador de impacto al final de la animación
hold on;
plot(x(end), y(end), 'ro', 'MarkerSize', 8, 'DisplayName', 'Impacto con el suelo');
legend('Location', 'best');
hold off;


7 .- Cálculo del caudal volumétrico.

La presa cuenta con una compuerta a una altura de 25m con un área de 1100m2. Para calcular el caudal ideal de agua que fluye por dicha compuerta se usa el siguiente código. 

clear; clc;

% Parámetros que vamos a utilizar
a = 1100;               % Área compuerta(m^2)
g = 9.81;               % Gravedad(m/s^2)
Hc = 25;                % Altura compuerta(m)
v0 = sqrt(2 * g * Hc);  % Velocidad inicial (con Tonirrelli)

% Caudal de la presa en metros cúbicos por segundo
Q = ( a * v0)

% Caudal de la presa en litros por segundo
Q(l/s) = Q ./ 1000

Este cálculo nos da un caudal de 24,36 litros por segundo

8 .- Cálculo de la fuerza de la presión total y la presión por unidad de superficie.

clear; clc
rho = 1000; % Densidad del agua (kg/m^3)
g = 9.81;   % Gravedad (m/s^2)
H = 134;    % Altura de la presa (m)
r0 = 242;   %Radio de la presa en la altura máxima (metros)
b = 35;     % Factor de curvatura parabólica

% Resolución para la integración
num_theta = 100; % Divisiones en theta
num_z = 100;     % Divisiones en z

% Rango de variables
theta = linspace(3*pi/4, 5*pi/4, num_theta);
z = linspace(0, H, num_z);

% Matriz de resultados
force_total = 0;
area_total = 0;

% Integración sobre la superficie
for i = 1:numel(theta)
    for j = 1:numel(z)
        % Coordenadas en la superficie
        r = r0 + b * (1 - (z(j)^2 / H^2)); % Radio a cada z
        h = H - z(j); % Profundidad
        
        % Elemento de área diferencial en coordenadas cilíndricas
        dS = r * (theta(2) - theta(1)) * (z(2) - z(1)); 
        
        % Presión a cada punto
        P = rho * g * h; 
        
        % Contribución a la fuerza total
        force_total = force_total + P * dS; 
        
        % Acumulación del área para calcular presión media (si necesaria)
        area_total = area_total + dS;
    end
end

% Fuerza total de presión
disp(['Fuerza total de presión: ', num2str(force_total), ' N']);

% Presión promedio por unidad de superficie
pressure_per_unit_area = force_total / area_total;
disp(['Presión por unidad de superficie promedio: ', num2str(pressure_per_unit_area), ' Pa']);


9 .- Diferentes tipos de presas y su estabilidad.

Hay diferentes tipos de presas, en concreto hay 5: presas de gravedad , presas de arco, presas de arco-gravedad, presas bóveda o doble arco, y por último presas de contrafuertes o aligeradas.

Las presas de gravedad se caracterizan porque resisten la fuerza de el agua por su propio peso, son adecuadas para valles anchos y terrenos sólidos. Generalmente están hechas de hormigón o mampostería, son muy buenas frente al vuelco y el deslizamiento de distintas partes. La mayoría de ellas son rectas, y es importante que el suelo no erosione. Estas presas se suelen construir para resistir algunos de los terremotos más fuertes. Aunque los cimientos de una presa de gravedad se construyen para soportar el peso de la presa y de toda el agua, son bastante flexibles, ya que absorben una gran cantidad de energía y la envían a la corteza terrestre.


Como su nombre indica, en este tipo de presas la forma predominante es la de arco y es precisamente su curvatura la que resiste el empuje del agua. Es importante que se emplee un material de muy alta resistencia en los estribos de la cerrada -los laterales de la presa- ya que es allí donde se produce el mayor esfuerzo. Por eso, este tipo de presas está limitado por las condiciones topográficas, cuanto más simétrica mejor, y orográficas del terreno.

A su vez, dentro de las presas arco existen dos tipos más:

Las presas bóveda o de doble arco, en las que, por su forma curva en planta y alzado, el esfuerzo se transmite a los laterales y el fondo. Un ejemplo de este tipo de presas es la Almendra que, situada en España, es emblemática por ser la presa más alta del país.

Las presas de arco-gravedad son unas presas que tienen forma curva que transfiere la presión del agua hacia los estribos en las paredes del valle. Están hechas de hormigón y son ideales para valles estrechos y rocosos. Consumen menos material que las presas de gravedad, pero su diseño es mucho mas complejo y necesitan estribos sólidos para funcionar. La estabilidad de este tipo de presa depende de la resistencia de los estribos, de la forma del arco y de la presión del agua. Es una mezcla de la presa de arco y de la de gravedad.

Las presas de contrafuertes o aligeradas son algo similares a las de gravedad en el sentido en que, en cuanto a su mecanismo de resistencia, pero cuenta con una serie de contrafuertes para ofrecer estabilidad frente al deslizamiento y al vuelco. Se emplea así una cantidad menor de material que en las presas de gravedad, pero tienen una mayor complejidad técnica.

Figura 8. Tipos de presas. Fuente: Iberdrola.

10 .- Referencias.

https://ingcivileng.com/2020/03/24/embalse-presa-boveda-el-atazar-rio-lozoya-madrid/

https://www.iberdrola.com/conocenos/nuestra-actividad/energia-hidroelectrica/tipos-presa

https://www.canaldeisabelsegunda.es/documents/20143/787562/El+Atazar_OK.pdf/ea6f7a49-7ee0-b69b-b274-70980bb104ad?t=1567751073532

Obtenido de «https://mat.caminos.upm.es/w/index.php?title=La_presa_de_El_Atazar._Grupo_16&oldid=80511»