La presa de El Atazar (Grupo 1)
| Trabajo realizado por estudiantes | |
|---|---|
| Título | La presa de El Atazar (Grupo 1) |
| 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 | |
La presa de El Atazar, es una de las infraestructuras hidráulicas más importantes de España y la mayor de la Comunidad de Madrid. Construida entre 1968 y 1972, se encuentra en la Sierra Norte de Madrid, sobre el río Lozoya. Forma parte del sistema de abastecimiento de agua de la capital y de toda la región. Su principal función es abastecer agua potable a Madrid, especialmente en épocas de sequía, producir energía hidroeléctrica y actúa como control de crecidas para regular los caudales del río Lozoya.
Esta presa, de tipo arco de gravedad, tiene una altura de 134 metros, lo que la convierte en una de las más altas del país, y una longitud de coronación de 484 metros. Su embalse tiene una capacidad de 425 hectómetros cúbicos. La presa del Atazar destaca por ser un símbolo de la ingeniería civil española, su relevancia la posiciona como un pilar fundamental para la sostenibilidad hídrica de la Comunidad de Madrid.
El objetivo principal de este trabajo es representar y visualizar la geometría de la presa como paso previo a un análisis detallado de su estabilidad estructural y de la interacción con el agua, considerando factores como la presión y el caudal. Para ello, se utilizará el software Matlab, especializado en programación y cálculos numéricos.
Contenido
- 1 Representación
- 2 Campo escalar de presión
- 3 Campo vectorial de la fuerza de presión
- 4 Representación de la trayectoria de la gota
- 5 Representación de los campos tangente y normal
- 6 Animación del vector de velocidad y aceleración sobre la curva
- 7 Caudal volumétrico
- 8 Fuerza de presión total y la presión por unidad de superficie
1 Representación
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]θ ∈ [\frac{3π}{4}, \frac{5π}{4}][/math] y [math]Z ∈ [0,H][/math].
Parámetros iniciales de la presa y el fluido
r0 = 308; % Radio base de la presa (aproximado)
b = 35; % Curvatura del arco parabólico
H = 134; % Altura de la presa
theta = linspace(3*pi/4, 5*pi/4, 100); % Ángulo en radianes
z = linspace(0, H, 100); % Altura en metros
% Creación del mallado en 2D para los ángulos y las alturas
[Z, Theta] = meshgrid(z, theta);
% Cálculo del radio para cada punto en la superficie
R = r0 + b * (1 - (Z.^2) / H^2);
% Conversión de coordenadas cilíndricas a cartesianas
X = R .* cos(Theta);
Y = R .* sin(Theta);
% Crear la figura
figure;
surf(X, Y, Z, 'FaceColor', 'b', 'EdgeColor', 'none'); % Gráfico 3D de la superficie
xlabel('Eje X (m)');
ylabel('Eje Y (m)');
zlabel('Eje Z (m)');
title('Superficie parametrizada de la presa');
% Configuración visual
axis equal; % Asegura la misma escala en todos los ejes
view(3); % Vista en 3D
grid on; % Muestra la cuadrícula2 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. En la imagen se observan dos gamas de colores, la primera corresponde a tonos fríos, que representa las zonas de menor presión y la segunda, compuesta por tonos cálidos, representa 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.
% Parámetros físicos
rho = 1000; % Densidad del agua (kg/m^3)
g = 9.81; % Aceleración de la gravedad (m/s^2)
% Cálculo de la presión en función de la altura
P = rho * g * (H - Z); % Presión hidrostática en la superficie
% Visualización de la superficie con colores representando la presión
figure;
surf(X, Y, Z, P, 'EdgeColor', 'none'); % Superficie con mapa de colores basado en P
xlabel('Eje X (m)');
ylabel('Eje Y (m)');
zlabel('Eje Z (m)');
title('Presión sobre la superficie de la presa');
% Configuración de colores y visualización
colormap(jet); % Paleta de colores para la presión
colorbar; % Barra de colores para interpretar la presión
axis equal; % Igualar las escalas de los ejes
view(3); % Vista tridimensional
grid on; % Mostrar la cuadrícula3 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;
% 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.
4 Representación de la trayectoria de la gota
A continuación, determinamos la curva plana que describe la trayectoria de una gota de agua al salir de la compuerta con los datos proporcionados, suponiendo que el agua es un fluido ideal (sin resistenciadel aire). Teniendo en cuenta la velocidad inicial del agua y la aceleración de 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,'DisplayName','gota');
hold on
plot(x_ground, 0, 'ro', 'MarkerSize', 8, 'DisplayName', 'Llegada');
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]);5 Representación de los campos tangente y normal
Representamos el campo tangente y el campo normal en varios puntos de la curva descrita por la gota de agua durante los primeros 20 segundos.
% Parámetros iniciales
angulo = deg2rad(15); % Ángulo de lanzamiento en radianes
gravedad = 9.81; % Aceleración debido a la gravedad (m/s^2)
altura_inicial = 25; % Altura inicial (m)
velocidad_inicial = sqrt(2 * gravedad * altura_inicial); % Velocidad inicial
% Tiempo y trayectorias
tiempo = linspace(0, 5, 500); % Intervalo de tiempo
pos_x = velocidad_inicial * cos(angulo) * tiempo; % Movimiento horizontal
pos_y = altura_inicial + velocidad_inicial * sin(angulo) * tiempo - 0.5 * gravedad * tiempo.^2; % Movimiento vertical
% Cálculo de las velocidades (derivadas de la posición)
vel_x = velocidad_inicial * cos(angulo) * ones(size(tiempo)); % Velocidad en x
vel_y = velocidad_inicial * sin(angulo) - gravedad * tiempo; % Velocidad en y
velocidad_total = sqrt(vel_x.^2 + vel_y.^2); % Magnitud de la velocidad
% Vectores unitarios tangente
tangente_x = vel_x ./ velocidad_total;
tangente_y = vel_y ./ velocidad_total;
% Cálculo de las aceleraciones (derivadas de la velocidad)
acel_x = zeros(size(tiempo)); % Aceleración en x (sin aceleración)
acel_y = -gravedad * ones(size(tiempo)); % Aceleración en y (gravedad)
aceleracion_total = sqrt(acel_x.^2 + acel_y.^2); % Magnitud de la aceleración
% Vectores unitarios normales
normal_x = (acel_x - tangente_x .* (acel_x .* tangente_x + acel_y .* tangente_y));
normal_y = (acel_y - tangente_y .* (acel_x .* tangente_x + acel_y .* tangente_y));
normal_x=normal_x./(sqrt(normal_x.^2+normal_y.^2));
normal_y=normal_y./(sqrt(normal_x.^2+normal_y.^2));
% Graficar trayectoria
figure;
plot(pos_x, pos_y, 'b', 'LineWidth', 1.5); hold on;
% Campos tangente y normal en puntos específicos a lo largo de la trayectoria
num_vectores = 20; % Número de puntos a graficar
puntos = round(linspace(1, length(tiempo), num_vectores));
quiver(pos_x(puntos), pos_y(puntos), tangente_x(puntos), tangente_y(puntos), 0.3, 'g', 'LineWidth', 1.5, 'DisplayName', 'Tangente');
quiver(pos_x(puntos), pos_y(puntos), normal_x(puntos), normal_y(puntos), 0.3, 'r', 'LineWidth', 1.5, 'DisplayName', 'Normal');
% Detalles adicionales 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;6 Animación del vector de velocidad y aceleración sobre la curva
% Definir parámetros de la animación
t_anim = linspace(0, t_ground, 200); % Intervalo de tiempo para la animación
x_anim = v0 * cos(theta) * t_anim; % Posición en el eje X durante la animación
y_anim = Hc + v0 * sin(theta) * t_anim - 0.5 * g * t_anim.^2; % Posición en el eje Y
% Calcular las velocidades y aceleraciones en los puntos de la animación
vx_anim = v0 * cos(theta) * ones(size(t_anim)); % Velocidad en X
vy_anim = v0 * sin(theta) - g * t_anim; % Velocidad en Y
speed_anim = sqrt(vx_anim.^2 + vy_anim.^2); % Magnitud de la velocidad
% Aceleraciones constantes
ax_anim = zeros(size(t_anim)); % Aceleración en X
ay_anim = -g * ones(size(t_anim)); % Aceleración en Y
% Crear la figura para la animación
figure;
hold on;
plot(x, y, 'b', 'LineWidth', 1.5, 'DisplayName', 'Trayectoria'); % Trayectoria total
xlabel('Distancia horizontal (m)');
ylabel('Altura (m)');
title('Animación de una Gota de Agua');
grid on;
axis equal;
xlim([0, max(x)]);
ylim([0, Hc + 5]);
% Inicializar marcador para la gota
gota = plot(x_anim(1), y_anim(1), 'ro', 'MarkerSize', 10, 'DisplayName', 'Gota');
% Inicializar vectores de velocidad y aceleración
quiver_vel = quiver(x_anim(1), y_anim(1), vx_anim(1), vy_anim(1), 0.3, 'g', 'LineWidth', 1.5, 'MaxHeadSize', 1, 'DisplayName', 'Velocidad');
quiver_acc = quiver(x_anim(1), y_anim(1), ax_anim(1), ay_anim(1), 0.3, 'r', 'LineWidth', 1.5, 'MaxHeadSize', 1, 'DisplayName', 'Aceleración');
legend;
% Ejecutar la animación
for i = 1:length(t_anim)
% Actualizar la posición de la gota
set(gota, 'XData', x_anim(i), 'YData', y_anim(i));
% Actualizar los vectores de velocidad y aceleración
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 crear el efecto de animación
pause(0.05);
end7 Caudal volumétrico
% Definición de parámetros
area_compuerta = 1100; % Área de la compuerta en metros cuadrados (m^2)
aceleracion_gravedad = 9.81; % Aceleración debido a la gravedad en m/s^2
altura_salida = 25; % Altura a la que sale el agua en metros (m)
% Cálculo de la velocidad de salida utilizando la ecuación de Torricelli
velocidad_salida = sqrt(2 * aceleracion_gravedad * altura_salida);
% Cálculo del caudal volumétrico (en litros por segundo)
caudal_volumetrico = (area_compuerta * velocidad_salida) / 1000;
% Mostrar los resultados en consola
fprintf('Velocidad de salida (v0): %.2f m/s\n', velocidad_salida);
fprintf('Caudal volumétrico (Q): %.2f m^3/s\n', caudal_volumetrico);Con este código obtenemos los resultados solicitados, donde la velocidad de salida es de 22.15 m/s y el caudal volumétrico resulta ser 24.36 m³/s.
8 Fuerza de presión total y la presión por unidad de superficie
Utilizando los datos de curvatura simple y doble curvatura, desarrollamos un código en Matlab que nos permite calcular tanto la fuerza como la presión, facilitando así la comparación de ambos resultados.
% Datos
rho = 1000; % Densidad del agua (kg/m^3)
g = 9.81; % Gravedad (m/s^2)
H = 134; % Altura de la presa (m)
r0 = 308.124; % Radio en la altura máxima (m)
% Presa de doble curvatura (b = 35)
b = 35; % Curvatura (m)
r = @(x) r0 + b * (1 - x.^2 / H^2); % Radio
a = @(x) sqrt(r(x).^2 + (2 * b .* r(x) .* x / H^2).^2); % Diferencial del área
P = @(x) rho * g * (H - x); % Presión
I = @(x) P(x) .* a(x); % Integral de presión
A = pi/2 * integral(a, 0, H); % Área de la superficie curva
F_total = pi/2 * integral(I, 0, H); % Fuerza total de presión
f_porarea = F_total / A; % Presión por unidad de área
% Presa de curvatura simple (b = 0)
b = 0; % Curvatura (m)
r = @(x) r0 + b * (1 - x.^2 / H^2); % Radio redefinido para la nueva curvatura
a = @(x) sqrt(r(x).^2 + (2 * b .* r(x) .* x / H^2).^2); % Diferencial de área redefinido
P = @(x) rho * g * (H - x); % Presión (sin cambios)
A0 = pi/2 * r0 * H; % Área de la superficie curva para curvatura simple
F_total0 = pi/2 * r0 * rho * g * H^2 / 2; % Fuerza total de presión para curvatura simple
f_porarea0 = F_total0 / A0; % Presión por unidad de área para curvatura simple
fprintf("Fuerza total para la presa de doble curvatura: %e N\n", F_total);
fprintf("Presión por unidad de área para la presa de doble curvatura: %e N/m^2\n", f_porarea);
fprintf("Fuerza total para la presa de curvatura simple: %e N\n", F_total0);
fprintf("Presión por unidad de área para la presa de curvatura simple: %e N/m^2\n", f_porarea0);Para la configuración de doble curvatura, con b = 35, la fuerza total es de 4.76 × 10¹⁰ N y la presión por unidad de superficie alcanza los 6.55 × 10⁵ Pa. En cambio, en el caso de la curvatura simple, con b = 0, la fuerza total es de 4.26 × 10¹⁰ N y la presión es de 6.57 × 10⁵ Pa. Así que, la estructura con doble curvatura soporta una mayor presión.
