Diferencia entre revisiones de «La presa de El Atazar»

Revisión actual del 15:41 5 dic 2024

La presa de El Atazar, ubicada en el río Lozoya, es la mayor y más importante de la Comunidad de Madrid. Inaugurada en 1972, es una presa de gravedad y arco con 134 metros de altura y una capacidad de 425 hectómetros cúbicos, siendo clave para el abastecimiento de agua potable en Madrid y su área metropolitana.

El presente trabajo tiene por objeto la visualización y representación de la geometría de la presa para hacer un posterior análisis de la estabilidad estructural y la interacción del agua con la presa, incluyendo presión y caudal. Para su realización se empleará el software de programación y cálculo numérico Matlab.

1 Representación de la presa

En primer lugar, se representa la superficie aguas arriba de la presa. Esta tiene una altura de 134 metros, y está definida por [math]θ ∈ [\frac{3π}{4}, \frac{5π}{4}][/math] y [math]Z ∈ [0,H][/math]. Para su realización, se ha comenzado definiendo los parámetros básicos y discretizando el dominio con una altura de paso de h=100. Con el comando "meshgrid()" se construye una malla que permita parametrizar la superficie en coordenadas cilíndricas según la siguiente ecuación.

Por último, se convierte la parametrización a coordenadas cartesianas y con el comando "surf()" obtenemos la gráfica de la superficie. A continuación, se muestra el código empleado para la obtención de la presa y la imagen resultante:

Figura 1. Representación de la presa.

r0 = 200;                                      % Radio base de la presa (aproximado)
b = 35;                                        % Curvatura del arco parabólico
H = 134;                                       % Altura de la presa
theta_min = 3*pi/4;
theta_max = 5*pi/4;
theta = linspace(theta_min, theta_max, 100);   % Ángulo en radianes
z = linspace(0, H, 100);                       % Altura
[Theta, Z] = meshgrid(theta, z);               % Creación del mallado
R = r0 + b * (1 - (Z.^2) / H^2);               % Ecuación paramétrica para r
X = R .* cos(Theta);                           % Conversión a coordenadas cartesiana
Y = R .* sin(Theta);
figure;
surf(X, Y, Z, 'FaceColor', 'cyan', 'EdgeColor', 'none');
xlabel('Eje X (m)');
ylabel('Eje Y (m)');
zlabel('Eje Z (m)');
title('Superficie parametrizada de la presa (aguas arriba)');
axis equal;                                    
view(3);                                       %Visualizamos el gráfico en tres dimensiones
grid on;

2 Presiones sobre la presa

En este análisis tratamos al agua de la presa como un fluido estático, como consecuencia, el propio peso del agua genera una presión a su alrededor, que actúa en las paredes de nuestra presa.

A continuación, analizamos y visualizamos los efectos de la presión sobre la superficie de la presa mediante la representación de los campos escalares y vectoriales de presión.

2.1 Campo escalar de presión

El campo escalar de presiones viene dado por la función: [math]P(z)=ρgh(z)[/math]. donde ρ es la densidad del agua, g es la aceleración de la gravedad, y h(z) es la profundidad del agua. Representamos este campo para visualizar cómo varía la presión a lo largo de la superficie aguas arriba de la presa, lo que nos permite identificar las zonas de mayor y menor presión. Utilizamos el comando 'surf' y 'colormap' para la diferenciación de colores en P(z) y 'colorbar' para visualizar el rango de presiones.

Figura 2. Representación del campo escalar de persión sobre la presa.

% Parámetros de presión
rho = 1000;             % Densidad del agua (kg/m^3)
g = 9.81;               % Gravedad (m/s^2)
% Cálculo de la presión
P = rho * g * (H - Z);  % Presión en función de z
% Gráfica de la superficie con mapa de colores para la presión
figure;
surf(X, Y, Z, P, 'EdgeColor', 'none');
xlabel('X (m)');
ylabel('Y (m)');
zlabel('Z (m)');
title('Presión sobre la superficie de la presa');
colorbar;               % Barra de colores para indicar magnitudes de presión
colormap(jet);          % Esquema de colores
axis equal;
view(3);
grid on;

En la imagen se observan dos gamas de colores diferenciadas: una primera, correspondiente a tonos fríos (azul), que representa las zonas de menor presión, y una segunda, compuesta por tonos cálidos (amarillo), asociada a las regiones de mayor presión. Es coherente con el comportamiento hidrostático que los valores más altos de presión se localicen en la base de la presa y disminuyan progresivamente con la altura.

2.2 Campo vectorial de la fuerza de presión

En segundo lugar, se representa el campo vectorial de la fuerza de presión tanto sobre el paramento aguas arriba de la presa como en un plano que la corta verticalmente. Esto permite visualizar cómo la presión del agua genera fuerzas perpendiculares a la superficie de la presa, cuya intensidad aumenta con la profundidad (mayor presión en las zonas más profundas). Además, la representación en el plano de corte proporciona una representación detallada de la distribución de las fuerzas, facilitando un análisis preciso de su comportamiento a lo largo de la profundidad y en diferentes secciones de la presa.

Con la fuerza como [math]\vec F=P(z)*\vec n [/math] ,determinamos n como el producto cruzado de vectores tangentes a la superficie y lo multiplicamos por todos los componentes de -P(z). Usamos 'quiver3' para representar el campo vectorial en la gráfica de la superficie y en un plano vertical.

Figura 3. Representación del campo escalar de presiones aguas arriba.

Figura 3.1. Representación del campo escalar de presiones en un plano vertical.

% Derivadas numéricas para obtener vectores tangentes
[Rx_theta, Rx_z] = gradient(X);
[Ry_theta, Ry_z] = gradient(Y);
[Rz_theta, Rz_z] = gradient(Z);
% Producto cruzado de vectores tangentes para obtener 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_magnitude = sqrt(Nx.^2 + Ny.^2 + Nz.^2);
% Vector normal unitario
Nx_unit = Nx ./ N_magnitude;
Ny_unit = Ny ./ N_magnitude;
Nz_unit = Nz ./ N_magnitude;
% 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, Z, Fx, Fy, Fz, 0.5, 'Color', 'yellow'); % Escala ajustable con el factor 2
hold on;
surf(X, Y, Z, P, 'EdgeColor', 'none', 'FaceAlpha',1); % Superficie semitransparente
colorbar;
xlabel('X (m)');
ylabel('Y (m)');
zlabel('Z (m)');
title('Campo vectorial de la fuerza de presión sobre la presa');
view(3);
grid on;
% Fijamos theta = 0 (corte longitudinal)
theta_cut = 0; % Corte en el plano x-z
X_cut = X(:, 51); % Extraemos el corte longitudinal (theta cerca de 0)
Z_cut = Z(:, 51); % Correspondiente en altura
P_cut = P(:, 51); % Presión en el plano
% Vectores normales en el corte
Nx_cut = Nx(:, 51); % Componente x del vector normal
Nz_cut = Nz(:, 51); % Componente z del vector normal
N_magnitude_cut = sqrt(Nx_cut.^2 + Nz_cut.^2);
% Normalización del vector normal
Nx_unit_cut = Nx_cut ./ N_magnitude_cut;
Nz_unit_cut = Nz_cut ./ N_magnitude_cut;
% Campo de fuerza de presión en el corte
Fx_cut = -P_cut .* Nx_unit_cut;
Fz_cut = -P_cut .* Nz_unit_cut;
% Gráfica del campo vectorial en el corte
figure;
quiver(X_cut, Z_cut, Fx_cut, Fz_cut, 'b', 'LineWidth', 1.5); % Vectores de fuerza
hold on;
plot(X_cut, Z_cut, 'k-', 'LineWidth', 2); % Contorno de la presa en el plano
xlabel('X (m)');
ylabel('Z (m)');
title('Campo de fuerza de presión en el plano longitudinal');
grid on;

Nuevamente, los resultados obtenidos coinciden con lo esperado. Las presiones más altas se registran en la parte inferior de la presa, representadas en color amarillo, mientras que las más bajas se encuentran en la parte superior, en color azul. De esta manera, se observa como la presión aumenta conforme ascendemos de arriba hacia abajo.

3 Estudio de la trayectoria de una gota

Hasta ahora únicamente hemos visualizado la imagen de la curva de la presa, y la distribución de la presión en las paredes de esta. No obstante, la trayectoria de una masa ante la curvatura de las paredes, la fuerza de gravedad y la presión del agua, se vería distorsionada. Resultando la siguiente representación una imagen más aproximada a la realidad.

3.1 Representación de la trayectoria de la gota

Con la expresión de la velocidad de la gota [math]v_{0}=\sqrt{2gH_{c}}[/math] y v=15º como dato, calculamos x(t) e y(t) teniendo en cuenta la aceleración gravitatoria.

En Matlab crearemos un vector de tiempo para calcular x e y y dispondremos la gráfica con la acción 'plot'.

Figura 4. Representación de la gota.

% Parámetros
theta = deg2rad(15); % Ángulo de salida en radianes
g = 9.81;           % Aceleración de la gravedad (m/s^2)
Hc = 25;            % Altura inicial (m)

% Velocidad inicial (previamente calculada)
v0 = sqrt(2 * g * Hc);

% Funciones de posición
t = linspace(0, 5, 1000); % Tiempo (ajustar si es necesario)
x = v0 * cos(theta) * t; 
y = Hc + v0 * sin(theta) * t - 0.5 * g * t.^2;

% Encontrar el tiempo cuando y(t) = 0 (impacto con el suelo)
t_ground = fzero(@(t) Hc + v0 * sin(theta) * t - 0.5 * g * t^2, [0, 10]);
x_ground = v0 * cos(theta) * t_ground; % Alcance horizontal

% Gráfica de la trayectoria
figure;
plot(x, y, 'b', 'LineWidth', 1.5);
hold on;
plot(x_ground, 0, 'ro', 'MarkerSize', 8, 'DisplayName', 'Impacto');
xlabel('Distancia horizontal (m)');
ylabel('Altura (m)');
title('Trayectoria de la gota de agua');
grid on;
legend;
xlim([0, x_ground + 10]);
ylim([0, Hc + 5]);

Como se puede observar, la presión con la que sale la gota de la presa hace que su altura aumente ligeramente en los primeros metros, para luego caer por la gravedad.

3.2 Representación de los campos tangente y normal

Para obtener el campo tangente derivamos la posición respecto al tiempo. Tras ello calculamos el vector normal que funciona como vector perpendicular. Por último utilizamos el comando 'quiver' para representar el gráfico de los vectores anteriores en los puntos seleccionados.

Figura 5. Representación de los campos tangente y normal.

% Parámetros
theta = deg2rad(15);
g = 9.81;
Hc = 25;
v0 = sqrt(2 * g * Hc);

% Trayectoria
t = linspace(0, 5, 500); % Tiempo
x = v0 * cos(theta) * t;
y = Hc + v0 * sin(theta) * t - 0.5 * g * t.^2;

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

% Tangente unitario
Tx = vx ./ speed;
Ty = vy ./ speed;

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

% Normal unitario
Nx = (ax - Tx .* (ax .* Tx + ay .* Ty)) ./ a_magnitude;
Ny = (ay - Ty .* (ax .* Tx + ay .* Ty)) ./ a_magnitude;

% Graficar trayectoria
figure;
plot(x, y, 'b', 'LineWidth', 1.5); hold on;

% Campos tangente y normal en puntos específicos
num_points = 20; % Número de vectores a graficar
indices = round(linspace(1, length(t), num_points));
quiver(x(indices), y(indices), Tx(indices), Ty(indices), 0.3, 'r', 'LineWidth', 1.5, 'DisplayName', 'Tangente');
quiver(x(indices), y(indices), Nx(indices), Ny(indices), 0.3, 'g', 'LineWidth', 1.5, 'DisplayName', 'Normal');

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

3.3 Animación de los vectores sobre la curva de la trayectoria

En Matlab creamos un bucle con el comando 'for' para determinar la posición, la velocidad y la aceleración en cada momento y utilizamos 'plot' y 'quiver' para mostrar la gráfica con los cambios aplicados.

Figura 6 . Representación animada de una gota

% Parámetros de la animación
t_anim = linspace(0, t_ground, 200); % Tiempo para animación
x_anim = v0 * cos(theta) * t_anim;
y_anim = Hc + v0 * sin(theta) * t_anim - 0.5 * g * t_anim.^2;

% Velocidades y aceleraciones en los puntos de la animación
vx_anim = v0 * cos(theta) * ones(size(t_anim));
vy_anim = v0 * sin(theta) - g * t_anim;
speed_anim = sqrt(vx_anim.^2 + vy_anim.^2);

% Aceleración constante en y
ax_anim = zeros(size(t_anim));
ay_anim = -g * ones(size(t_anim));

% Crear la figura
figure;
hold on;
plot(x, y, 'b', 'LineWidth', 1.5, 'DisplayName', 'Trayectoria'); % Trayectoria completa
xlabel('Distancia horizontal (m)');
ylabel('Altura (m)');
title('Animación: Gota de Agua');
grid on;
axis equal;
xlim([0, max(x)]);
ylim([0, Hc + 5]);

% Añadir marcador para la gota
gota = plot(x_anim(1), y_anim(1), 'ro', 'MarkerSize', 10, 'DisplayName', 'Gota');

% Añadir vectores de velocidad y aceleración
quiver_vel = quiver(x_anim(1), y_anim(1), vx_anim(1), vy_anim(1), 0.3, 'r', 'LineWidth', 1.5, 'MaxHeadSize', 1, 'DisplayName', 'Velocidad');
quiver_acc = quiver(x_anim(1), y_anim(1), ax_anim(1), ay_anim(1), 0.3, 'g', 'LineWidth', 1.5, 'MaxHeadSize', 1, 'DisplayName', 'Aceleración');

legend;

% Animación
for i = 1:length(t_anim)
    % Actualizar marcador
    set(gota, 'XData', x_anim(i), 'YData', y_anim(i));
    
    % Actualizar vectores
    set(quiver_vel, 'XData', x_anim(i), 'YData', y_anim(i), ...
        'UData', vx_anim(i), 'VData', vy_anim(i));
    set(quiver_acc, 'XData', x_anim(i), 'YData', y_anim(i), ...
        'UData', ax_anim(i), 'VData', ay_anim(i));
    
    % Pausa para visualizar
    pause(0.05);
end

4 Caudal volumétrico

El caudal es un concepto con mucha relevancia en la realización de obras hidráulicas, por lo tanto, consideramos definirlo primero. Conocemos al caudal de una obra hidráulica como, a la cantidad de fluido que atraviesa una superficie en el espacio en una unidad de tiempo. Esto nos permite estimar la cantidad de agua a una velocidad aproximada que circula por la presa. La fórmula del caudal es [math]Q=A*v_{0}[/math], solo sustituimos los valores de las incógnitas.

% Parámetros
A = 1100;        % Área de la compuerta (m^2)
g = 9.81;        % Aceleración de la gravedad (m/s^2)
Hc = 25;         % Altura de salida del agua (m)

% Velocidad de salida (Torricelli)
v0 = sqrt(2 * g * Hc);

% Caudal volumétrico
Q = (A * v0)/1000;

% Mostrar resultados
fprintf('Velocidad de salida (v0): %.2f m/s\n', v0);
fprintf('Caudal volumétrico (Q): %.2f m^3/s\n', Q);

Utilizando este código nos dan los resultados pedidos, la velocidad de salida es 22,15 m/s y el caudal volumétrico es 24,36 m^3/s.

5 Fuerza de presión total y la presión por unidad de superficie

La presión que ejerce el agua contenida por la presa a lo largo de esta genera una fuerza que intenta desplazar y girar nuestra estructura. Teniendo estos datos en cuenta, es posible concretar más detalles de la presa, como podrían ser, la forma, el asentamiento o los materiales que serán necesarios.

Con los datos de curvatura simple y doble curvatura creamos un código en matlab para que nos de el resultado de la fuerza y la presión y poder compararlos. Para calcular la fuerza de presión utilizamos la integral [math]F=\int P(z)dA[/math]. Discretizamos en Matlab y resolvemos la integral para cada valor de b. Para el caso de doble curvatura b=35 y para el caso de curvatura simple b=0.

% Parámetros de la presa
rho = 1000;        % Densidad del agua (kg/m^3)
g = 9.81;          % Gravedad (m/s^2)
H = 134;           % Altura de la presa (m)
r0 = 80;           % Radio base (m)
b1 = 35;     % Curvatura doble
b0 = 0;      % Curvatura simple
theta_range = linspace(-pi/4, pi/4, 100); % Ángulo

% Definir función de radio (diferentes curvaturas)
r_1 = @(z) r0 + b1 * (1 - (z / H)^2);
r_0 = @(z) r0;

% Cálculo de área y fuerza (doble curvatura)
z_range = linspace(0, H, 100); % Alturas
F_1 = 0;
A_1 = 0;
for z = z_range
    r = r_1(z);
    dA = 2 * pi * r * H / length(z_range); % Aproximación diferencial
    P = rho * g * (H - z); % Presión a esa profundidad
    F_1 = F_1 + P * dA; % Sumar contribución
    A_1 = A_1 + dA; % Sumar área
end
p_s_1 = F_1 / A_1; % Presión por unidad de área

% Cálculo de área y fuerza (curvatura simple)
F_0 = 0;
A_0 = 0;
for z = z_range
    r = r_0(z);
    dA = 2 * pi * r * H / length(z_range); % Aproximación diferencial
    P = rho * g * (H - z); % Presión a esa profundidad
    F_0 = F_0 + P * dA; % Sumar contribución
    A_0 = A_0 + dA; % Sumar área
end
p_s_0 = F_0 / A_0; % Presión por unidad de área

% Resultados
fprintf('Doble curvatura:\n');
fprintf(' - Fuerza total: %.2e N\n', F_1);
fprintf(' - Presión por unidad de superficie: %.2e Pa\n\n', p_s_1);

fprintf('Curvatura simple:\n');
fprintf(' - Fuerza total: %.2e N\n', F_0);
fprintf(' - Presión por unidad de superficie: %.2e Pa\n', p_s_0);

En el caso de la doble curvatura, b=35, la fuerza total es de [math]6,04*10^{10}[/math] N y su presión por unidad de superficie es de [math]6,95*10^{5}[/math]Pa. En la de curvatura simple, b=0, la fuerza total es de [math]4,43*10^{10}[/math] N y la presión de [math]6,57*10^{5}[/math]Pa. Por lo tanto soportaría más presión la configuración de doble curvatura.

6 Tipos de presa, estabilidad y materiales

A la hora de tratar con una presa debemos saber de cual tipo estamos hablando, para ello, procederemos a clasificarlas y señalar los detalles más relevantes de cada una.

A grandes rasgos podemos clasificar estas obras en : presas de fábrica y presas de materiales sueltos.

Las presas de fábrica engloban al grupo de presas de gravedad, presas de contrafuertes y presas de arco. Todas están caracterizadas por estar construidas con hormigón en masa principalmente, y poder soportar altas cargas.

Comenzando con las de gravedad, el principal mecanismo que usa es el rozamiento de la presa con el terreno debido al peso de esta. Además, para que la presa no vuelque la resultante de los empujes del agua y el peso propio debe estar contenida en la base del cuerpo de presa. La anchura de la base suele ser 0.8 veces la altura de la presa.

Las presas de contrafuerte son muy parecidas a las de gravedad, sin embargo, en estas el "cuerpo" se apoya sobre unos elementos más largos llamados "cola". A diferencia de las de gravedad, carecen prácticamente de subpresión, que junto al peso de la cuña de agua sobre el cuerpo, hace que tengan mucha menos masa de hormigón. La anchura de la base suele ser igual que la altura de la presa.

Aquellas con forma curva se denominarán presas de arco, se caracterizan por ser antifuniculares de la fuerza ejercida por el agua, transmitiendo esta fuerza mediante un axil hasta los extremos de la cerrada, lo que obliga al asentamiento a poseer unas buenas características resistivas.

Este tipo puede dividirse a su vez en tres:

- presa arco (planta curva, sección transversal trapezoidal)

- presa de doble curvatura o en bóveda (planta curva sección curva)

- presa arco gravedad (presa de arco con estribos que refuerzan la estructura)

Estos tipos son más esbeltas que las de gravedad, ocupando las de arcos y bóveda una base 0.3 veces la altura y las de arco-gravedad 0.5 veces.

Por otra parte, nos encontramos con las presas de materiales sueltos, dirigidas a construcciones con menos capacidades resistivas pero más abundantes por su versatilidad en asentamientos. Estas son menos esbeltas que las de fábrica, mas, su característica principal es la zonificación de materiales, es decir, cada tipo de material se ubica de forma que se optimice su función. En estas obras los materiales más usados son arcillas, gravas, arenas, limos o y polvo de roca.

Las subdividiremos en tres grupos:

- Presas homogéneas: utilizan el mismo material impermeable en todo el cuerpo de presa. Se emplean en alturas pequeñas y medias siendo el ancho de su base en torno a 6 y 8 veces la altura de presa.

- Presas de núcleo: contienen un núcleo impermeable en su interior, cubierto por un material resistente tipo escollera. Se construyen de todos los tamaños con una anchura de base entre 4 y 6 veces la altura de presa.

- Presa de pantalla: aquellas protegidas por una pantalla impermeable que se apoya sobre el talud de aguas arriba del cuerpo de presa. El ancho de la base suele ser 3 veces la altura de presa.