La presa de El Atazar
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.
| Trabajo realizado por estudiantes | |
|---|---|
| Título | La presa de El Atazar. Grupo 31 |
| Asignatura | Teoría de Campos |
| Curso | 2024-25 |
| Autores |
|
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Contenido
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:
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.
% 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.
% 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
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.
% 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]);
3.2 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
% 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.
% 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.
% 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 6,04*10^10 N y su presión por unidad de superficie es de 6,95*10^5 Pa. En la de curvatura simple, b=0, la fuerza total es de 4,43*10^10 N y la presión de 6,57*10^5 Pa. Por lo tanto soportaría más presión la configuración de doble curvatura.
