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 discreteando 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 presay la imagen resultante:
clc
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
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 el campo escalar de presiones para visualizar como varía la presión en la superficie aguas arriba de la presa.Esto nos permite identificar las zonas de mayor y menor presión.
% 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
% 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;
%APARTADO B
% 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;
3 Estudio de la trayectoria de una gota
3.1 Representación de la trayectoria 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]);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);
end4 Caudal volumétrico
% 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.
