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	Cristopher Ardon Colindres Andrea Garcia Carrasco Aaron García Martín Lara Gutiérrez Kreutzer
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

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.

[math]r = r_{0} + b * (1 - \frac{z^2}{h^2})[/math]

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:

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

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.

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;

INTERPRETACIÓN DEL CAMPO ESCALAR DE PERSIÓN

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

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;

%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

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]);

3.2 Representación de los campos tangente y normal

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

% 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

% 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

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 [math]6,04exp(10)[/math] 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.