Diferencia entre revisiones de «La presa de Atazar NEMJJ»

De MateWiki
Saltar a: navegación, buscar
(→‎Cálculo del caudal de entrada del río)
 
(No se muestran 40 ediciones intermedias del mismo usuario)
Línea 19: Línea 19:
 
*Juan Martín Rodríguez }}
 
*Juan Martín Rodríguez }}
  
−
==Modelo geométrico de la presa ==
+
===Modelo geométrico de la presa ===
  
 
Para el análisis se considerará la superficie de la presa en el lado aguas arriba, es decir, la zona en contacto directo con el agua. Se supone que la sección transversal puede aproximarse mediante un arco de circunferencia con eje de simetría en el valle, mientras que la sección longitudinal se comporta como un arco parabólico. Trabajaremos durante la ejecución de este trabajo en coordenadas cilíndricas (r, θ, z).  
 
Para el análisis se considerará la superficie de la presa en el lado aguas arriba, es decir, la zona en contacto directo con el agua. Se supone que la sección transversal puede aproximarse mediante un arco de circunferencia con eje de simetría en el valle, mientras que la sección longitudinal se comporta como un arco parabólico. Trabajaremos durante la ejecución de este trabajo en coordenadas cilíndricas (r, θ, z).  
Línea 30: Línea 30:
  
  
−
==Representación de la presa ==
+
===Representación de la presa ===
  
−
==Flujo de entrada del río ===
+
==Flujo de entrada del río ==
 +
 
 +
En esta sección analizaremos el comportamiento del flujo de entrada del río Lozoya en el embalse. Para ello, definiremos un campo de velocidad que modela la entrada del río por extremo norte. Inicialmente, consideraremos un modelo simple en dos dimensiones del campo de velocidad, el cual ampliaremos más tarde a tres dimensiones para incorporar una componente vertical, ofreciendo así una representación más real del flujo de entrada.
 +
 
 +
El estudio de este campo vectorial se lleva a cabo mediante su visualización y el análisis de sus propiedades diferenciales fundamentales. Específicamente, calcularemos y representaremos su divergencia y su rotacional. Finalmente, y como aplicación práctica, calcularemos el caudal de entrada bajo una condición de flujo, demostrando la utilidad del modelo para la medición de variables hidrológicas clave.
 +
 
 +
<center>Campo de velocidad del flujo de entrada (2D): <math>\vec{v}(x, y) = v_0 \, e^{-\frac{x^2 + y^2}{R^2}} (\cos \varphi \, \vec{i} + \sin \varphi \, \vec{j})</math>
 +
 
 +
 
 +
Campo de velocidad del flujo de entrada (3D): <math>\vec{v}(x, y, z) = v_0 e^{-\frac{x^2 + y^2}{R^2}} (\cos \varphi \, \vec{i} + \sin \varphi \, \vec{j}) + v_1 \sin\left(\pi \frac{z}{H_{\text{agua}}}\right) e^{-\frac{x^2 + y^2}{R^2}} \vec{k}</math></center>
 +
 
 +
El campo de velocidad del flujo de entrada se modela con los parámetros \( v_0 = 0.5 \, \text{m/s} \) para la velocidad máxima en el centro de entrada, \( v_1 = 0.1 \, \text{m/s} \) para la componente vertical máxima, \( R = 30 \, \text{m} \) como radio característico, y \( \varphi = \frac{\pi}{6} \) como ángulo de entrada. El dominio de estudio comprende las coordenadas horizontales \( (x, y) \in [-100, 100]^2 \, \text{m} \) con la restricción \( \sqrt{x^2 + y^2} \leq \rho_a \), donde \( \rho_a \) es el radio de la presa a la altura \( z = H_{\text{agua}} \), mientras que la coordenada vertical abarca el intervalo \( z \in [0, H_{\text{agua}}] \).
 +
 
 +
===Representación del campo vectorial de velocidad===
 +
 
 +
En este apartado realizaremos la visualización completa del campo de velocidad asociado al flujo de entrada del río Lozoya en el embalse de El Atazar. Apoyándonos en dos imágenes claves: el plano horizontal de la superficie libre del agua y el plano vertical para así saber la variación con la profundidad.
 +
[[Archivo:3_5_imagen_muni.png|10000px|miniaturadeimagen|Representación del campo vectorial de la velocidad]]
 +
 
 +
{{matlab|codigo=
 +
% Parámetros del problema
 +
v0 = 0.5;      % m/s
 +
v1 = 0.1;      % m/s
 +
R = 30;        % m
 +
phi = pi/6;    % ángulo de entrada
 +
H_agua = 125;  % m
 +
 
 +
% Cálculo de rho_a (radio de la presa a z = H_agua)
 +
H = 134;      % altura total de la presa
 +
rho0 = 150;    % m
 +
b = 40;        % m
 +
% Fórmula: rho = rho0 + b*(1 - z^2/H^2)
 +
rho_a = rho0 + b*(1 - (H_agua^2)/(H^2));
 +
 
 +
% Plano horizontal (z = H_agua - nivel superficial)
 +
x = linspace(-100, 100, 20);
 +
y = linspace(-100, 100, 20);
 +
[X, Y] = meshgrid(x, y);
 +
Z_superficie = H_agua * ones(size(X));
 +
% Campo de velocidad en superficie (z = H_agua)
 +
r_cuad = X.^2 + Y.^2;
 +
factor_exp = exp(-r_cuad/R^2);
 +
% Componentes del vector velocidad
 +
Vx = v0 * factor_exp * cos(phi);
 +
Vy = v0 * factor_exp * sin(phi);
 +
Vz = v1 * sin(pi * Z_superficie/H_agua) .* factor_exp;
 +
% Filtrar puntos fuera del dominio (dentro de la presa)
 +
dentro_presa = sqrt(X.^2 + Y.^2) <= rho_a;
 +
Vx(~dentro_presa) = NaN;
 +
Vy(~dentro_presa) = NaN;
 +
Vz(~dentro_presa) = NaN;
 +
% Graficar campo vectorial en plano horizontal
 +
figure('Position', [100, 100, 800, 600]);
 +
subplot(2,2,1);
 +
quiver(X, Y, Vx, Vy, 2, 'b');
 +
axis equal;
 +
xlabel('x (m)');
 +
ylabel('y (m)');
 +
title('Campo de velocidad - Plano Horizontal (z = H_{agua})');
 +
grid on;
 +
 
 +
% Plano vertical (y = 0 - sección longitudinal)
 +
x_vert = linspace(-100, 100, 15);
 +
z_vert = linspace(0, H_agua, 15);
 +
[X_vert, Z_vert] = meshgrid(x_vert, z_vert);
 +
Y_vert = zeros(size(X_vert));
 +
% Campo de velocidad en plano vertical y=0
 +
r_cuad_vert = X_vert.^2 + Y_vert.^2;
 +
factor_exp_vert = exp(-r_cuad_vert/R^2);
 +
Vx_vert = v0 * factor_exp_vert * cos(phi);
 +
Vz_vert = v1 * sin(pi * Z_vert/H_agua) .* factor_exp_vert;
 +
% Filtrar puntos fuera del dominio
 +
dentro_presa_vert = sqrt(X_vert.^2 + Y_vert.^2) <= rho_a;
 +
Vx_vert(~dentro_presa_vert) = NaN;
 +
Vz_vert(~dentro_presa_vert) = NaN;
 +
subplot(2,2,2);
 +
quiver(X_vert, Z_vert, Vx_vert, Vz_vert, 2, 'b');
 +
xlabel('x (m)');
 +
ylabel('z (m)');
 +
title('Campo de velocidad - Plano Vertical (y = 0)');
 +
grid on;
 +
}}
 +
Plano Horizontal:
 +
El agua entra formando un chorro que va de noreste a suroeste. En el centro va más rápido (0.5 m/s) y se va haciendo más lento hacia los lados, como un abanico que se abre.
 +
 
 +
Plano Vertical:
 +
En profundidad, el agua se mueve más lento en el fondo y la superficie, y más rápido a media profundidad. También hay movimiento vertical que mezcla el agua.
 +
 
 +
===La divergencia del campo===
 +
 
 +
En este apartado realizaremos el cálculo de la divergencia del campo de velocidad, para estudiar el comportamiento del flujo. En el apartado físico la divergencia representa la compresión o expansión del fluido (divergencia positiva: partículas se separan / divergencia negativa: partículas se juntan). Representaremos la divergencia en su plano vertical y en su plano horizontal.
 +
 
 +
[[Archivo:3_6_imagen_h.png|600px|miniaturadeimagen|Representación plano horizontal de la divergencia.]]
 +
[[Archivo:3_6_imagen_v.png|600px|miniaturadeimagen|Representación plano vertical de la divergencia.]]
 +
{{matlab|codigo=
 +
clear,clc
 +
% Parámetros
 +
v0 = 0.5; v1 = 0.1; R = 30; phi = pi/6; H_agua = 125;
 +
H = 134; rho0 = 150; b = 40;
 +
rho_a = rho0 + b*(1 - (H_agua^2)/(H^2));
 +
% DIVERGENCIA EN PLANO HORIZONTAL (z = H_agua constante)
 +
x_div = linspace(-80, 80, 40);
 +
y_div = linspace(-80, 80, 40);
 +
[X_div, Y_div] = meshgrid(x_div, y_div);
 +
Z_div = H_agua * ones(size(X_div));
 +
% Campo de velocidad completo
 +
r_cuad = X_div.^2 + Y_div.^2;
 +
factor_exp = exp(-r_cuad/R^2);
 +
Vx = v0 * factor_exp * cos(phi);
 +
Vy = v0 * factor_exp * sin(phi);
 +
Vz = v1 * sin(pi * Z_div/H_agua) .* factor_exp;
 +
% CÁLCULO DE DERIVADAS PARCIALES
 +
dx = x_div(2) - x_div(1);
 +
dy = y_div(2) - y_div(1);
 +
% Derivadas usando gradient
 +
[dVx_dx, dVx_dy] = gradient(Vx, dx, dy);
 +
[dVy_dx, dVy_dy] = gradient(Vy, dx, dy);
 +
[dVz_dx, dVz_dy] = gradient(Vz, dx, dy);
 +
% Como z es constante en este plano, ∂Vz/∂z = 0
 +
divergencia_horizontal = dVx_dx + dVy_dy;
 +
% Filtrar puntos fuera del dominio
 +
dentro_presa = sqrt(X_div.^2 + Y_div.^2) <= rho_a;
 +
divergencia_horizontal(~dentro_presa) = NaN;
 +
Vx(~dentro_presa) = NaN;
 +
Vy(~dentro_presa) = NaN;
 +
figure('Position', [100, 100, 1200, 500]);
 +
% plot 1: Mapa de colores de divergencia
 +
plot(1);
 +
pcolor(X_div, Y_div, divergencia_horizontal);
 +
shading interp;
 +
colorbar;
 +
colormap(jet);
 +
hold on;
 +
% Contornos de divergencia
 +
contour(X_div, Y_div, divergencia_horizontal, 15, 'k', 'LineWidth', 0.5);
 +
x_vert = linspace(-80, 80, 30);
 +
z_vert = linspace(0, H_agua, 25);
 +
[X_vert, Z_vert] = meshgrid(x_vert, z_vert);
 +
Y_vert = zeros(size(X_vert));
 +
% Campo de velocidad en plano vertical
 +
r_cuad_vert = X_vert.^2 + Y_vert.^2;
 +
factor_exp_vert = exp(-r_cuad_vert/R^2);
 +
Vx_vert = v0 * factor_exp_vert * cos(phi);
 +
Vy_vert = v0 * factor_exp_vert * sin(phi);
 +
Vz_vert = v1 * sin(pi * Z_vert/H_agua) .* factor_exp_vert;
 +
% Derivadas en plano vertical
 +
dx_vert = x_vert(2) - x_vert(1);
 +
dz_vert = z_vert(2) - z_vert(1);
 +
[dVx_dx_vert, dVx_dz_vert] = gradient(Vx_vert, dx_vert, dz_vert);
 +
[dVy_dx_vert, dVy_dz_vert] = gradient(Vy_vert, dx_vert, dz_vert);
 +
[dVz_dx_vert, dVz_dz_vert] = gradient(Vz_vert, dx_vert, dz_vert);
 +
divergencia_vertical = dVx_dx_vert + dVz_dz_vert;
 +
% Filtrar puntos fuera del dominio
 +
dentro_presa_vert = sqrt(X_vert.^2 + Y_vert.^2) <= rho_a;
 +
divergencia_vertical(~dentro_presa_vert) = NaN;
 +
figure('Position', [100, 600, 1000, 400]);
 +
plot(2);
 +
pcolor(X_vert, Z_vert, divergencia_vertical);
 +
shading interp;
 +
colorbar;
 +
hold on;
 +
contour(X_vert, Z_vert, divergencia_vertical, 10, 'k', 'LineWidth', 0.5);
 +
xlabel('x (m)'); ylabel('z (m)');
 +
title('Divergencia ∇·v (s^{-1}) - Plano Vertical (y=0)');}}
 +
Después de analizar los resultados identificamos que el agua se comprime en la dirección de la velocidad (divergencia negativa) y se expande hacia los lados (divergencia positiva). Se puede concluir que se crea una zona de mayor presión justo en la zona donde el río se encuentra con el embalse.
 +
 
 +
===El rotacionaldel campo===
 +
 
 +
Este apartado analizaremos el rotacional del campo de velocidad, que mide cuánto gira el agua al entrar al embalse. Calcularemos si se forman remolinos o torbellinos, lo que afecta directamente cómo se mezcla el agua, se reparten los sedimentos y se oxigena el embalse.
 +
 
 +
[[Archivo:3_7_imagen_h.png|600px|miniaturadeimagen|Representación plano horizontal del rotacional.]]
 +
[[Archivo:3_7_imagen_v.png|600px|miniaturadeimagen|Representación plano vertical del rotacional.]]
 +
 
 +
{{matlab|codigo=
 +
clear,clc
 +
% Parámetros
 +
v0 = 0.5; v1 = 0.1; R = 30; phi = pi/6; H_agua = 125;
 +
H = 134; rho0 = 150; b = 40;
 +
rho_a = rho0 + b*(1 - (H_agua^2)/(H^2));
 +
%ROTACIONAL EN PLANO HORIZONTAL (z = H_agua constante)
 +
x_rot = linspace(-80, 80, 35);
 +
y_rot = linspace(-80, 80, 35);
 +
[X_rot, Y_rot] = meshgrid(x_rot, y_rot);
 +
Z_rot = H_agua * ones(size(X_rot));
 +
% Campo de velocidad
 +
r_cuad = X_rot.^2 + Y_rot.^2;
 +
factor_exp = exp(-r_cuad/R^2);
 +
Vx = v0 * factor_exp * cos(phi);
 +
Vy = v0 * factor_exp * sin(phi);
 +
Vz = v1 * sin(pi * Z_rot/H_agua) .* factor_exp;
 +
% Derivadas parciales
 +
dx = x_rot(2) - x_rot(1);
 +
dy = y_rot(2) - y_rot(1);
 +
[dVx_dx, dVx_dy] = gradient(Vx, dx, dy);
 +
[dVy_dx, dVy_dy] = gradient(Vy, dx, dy);
 +
[dVz_dx, dVz_dy] = gradient(Vz, dx, dy);
 +
% ROTACIONAL 3D: ∇×v = (∂Vz/∂y - ∂Vy/∂z, ∂Vx/∂z - ∂Vz/∂x, ∂Vy/∂x - ∂Vx/∂y)
 +
rot_x = dVz_dy - 0;  % ∂Vy/∂z ≈ 0 en z constante
 +
rot_y = 0 - dVz_dx;  % ∂Vx/∂z ≈ 0 en z constante
 +
rot_z = dVy_dx - dVx_dy;  % Componente vertical del rotacional
 +
% Filtrar puntos fuera del dominio
 +
dentro_presa = sqrt(X_rot.^2 + Y_rot.^2) <= rho_a;
 +
rot_x(~dentro_presa) = NaN;
 +
rot_y(~dentro_presa) = NaN;
 +
rot_z(~dentro_presa) = NaN;
 +
Vx(~dentro_presa) = NaN;
 +
Vy(~dentro_presa) = NaN;
 +
% REPRESENTACIÓN GRÁFICA - PLANO HORIZONTAL
 +
plot(1);
 +
rot_horiz_mag = sqrt(rot_x.^2 + rot_y.^2);
 +
pcolor(X_rot, Y_rot, rot_horiz_mag);
 +
shading interp;
 +
colorbar;
 +
hold on;
 +
quiver(X_rot, Y_rot, rot_x, rot_y, 2, 'w', 'LineWidth', 1.5);
 +
xlabel('x (m)'); ylabel('y (m)');
 +
title('Rotacional Horizontal |∇×v|_{xy} - s^{-1}');
 +
axis equal;
 +
%ROTACIONAL EN PLANO VERTICAL (y = 0) - CORREGIDO
 +
x_vert = linspace(-80, 80, 30);
 +
z_vert = linspace(0, H_agua, 25);
 +
[X_vert, Z_vert] = meshgrid(x_vert, z_vert);
 +
Y_vert = zeros(size(X_vert));
 +
% Campo de velocidad en plano vertical
 +
r_cuad_vert = X_vert.^2 + Y_vert.^2;
 +
factor_exp_vert = exp(-r_cuad_vert/R^2);
 +
Vx_vert = v0 * factor_exp_vert * cos(phi);
 +
Vy_vert = v0 * factor_exp_vert * sin(phi);
 +
Vz_vert = v1 * sin(pi * Z_vert/H_agua) .* factor_exp_vert;
 +
% Derivadas en plano vertical
 +
dx_vert = x_vert(2) - x_vert(1);
 +
dz_vert = z_vert(2) - z_vert(1);
 +
% Calcular todas las derivadas necesarias
 +
[dVx_dx_vert, dVx_dz_vert] = gradient(Vx_vert, dx_vert, dz_vert);
 +
[dVy_dx_vert, dVy_dz_vert] = gradient(Vy_vert, dx_vert, dz_vert);
 +
[dVz_dx_vert, dVz_dz_vert] = gradient(Vz_vert, dx_vert, dz_vert);
 +
% Rotacional en plano vertical y=0:
 +
rot_x_vert = 0 - dVy_dz_vert;                    % ω_x = ∂Vz/∂y - ∂Vy/∂z ≈ -∂Vy/∂z
 +
rot_y_vert = dVx_dz_vert - dVz_dx_vert;          % ω_y = ∂Vx/∂z - ∂Vz/∂x
 +
rot_z_vert = dVy_dx_vert - 0;                    % ω_z = ∂Vy/∂x - ∂Vx/∂y ≈ ∂Vy/∂x
 +
% Filtrar puntos fuera del dominio
 +
dentro_presa_vert = sqrt(X_vert.^2 + Y_vert.^2) <= rho_a;
 +
rot_x_vert(~dentro_presa_vert) = NaN;
 +
rot_y_vert(~dentro_presa_vert) = NaN;
 +
rot_z_vert(~dentro_presa_vert) = NaN;
 +
Vx_vert(~dentro_presa_vert) = NaN;
 +
Vz_vert(~dentro_presa_vert) = NaN;
 +
% REPRESENTACIÓN PLANO VERTICAL
 +
figure('Position', [100, 600, 1200, 500]);
 +
plot(2);
 +
rot_vert_mag = sqrt(rot_x_vert.^2 + rot_y_vert.^2 + rot_z_vert.^2);
 +
pcolor(X_vert, Z_vert, rot_vert_mag);
 +
shading interp;
 +
colorbar;
 +
hold on;
 +
quiver(X_vert, Z_vert, rot_x_vert, rot_y_vert, 2, 'w', 'LineWidth', 1.5);
 +
contour(X_vert, Z_vert, rot_vert_mag, 10, 'k', 'LineWidth', 0.5);
 +
xlabel('x (m)'); ylabel('z (m)');
 +
title('Magnitud Total |∇×v| - Plano Vertical - s^{-1}');}}
 +
 
 +
Analizando los rotacionales se puede concluir que existe vorticidad en flujo. Como el agua no entra toda igual de rápido ni en la misma dirección, se generan pequeños remolinos naturales que ayudan a mezclarla.
 +
Los cálculos muestran que el rotacional no es cero, lo que indica que el agua presenta movimiento rotacional al entrar al embalse. Esta rotación favorece la formación de remolinos y torbellinos que ayudan a mezclar el agua del río con la del embalse, mejorando la oxigenación en la zona de entrada.
 +
 
 +
 
 +
===Cálculo del caudal de entrada del río===
 +
 
 +
En este apartado calcularemos el caudal de entrada del río Lozoya al embalse, es decir, la cantidad de agua que llega por segundo. Simplificaremos el cálculo eliminando la componente vertical (v1 = 0 m/s) para centrarnos en el flujo horizontal principal. A través de integración matemática del campo de velocidad sobre una sección perpendicular a la dirección de entrada, determinaremos cuántos metros cúbicos por segundo aporta el río al embalse, un dato esencial para la gestión hídrica y el balance de agua de toda la presa.
 +
 
 +
{{matlab|codigo=
 +
clear; clc;
 +
% Parámetros
 +
v0 = 0.5;          % Velocidad máxima [m/s]
 +
v1 = 0;            % Velocidad vertical = 0 [m/s]
 +
R = 30;            % Radio característico [m]
 +
phi = pi/6;        % Ángulo de entrada [rad]
 +
H_agua = 125;      % Altura del agua [m]
 +
H = 134;          % Altura total presa [m]
 +
rho0 = 150;        % Radio en coronación [m]
 +
b = 40;            % Parámetro de curvatura [m]
 +
rho_a = rho0 + b*(1 - (H_agua^2)/(H^2));
 +
nx = 80; ny = 80; nz = 40;
 +
x = linspace(-rho_a, rho_a, nx);
 +
y = linspace(-rho_a, rho_a, ny);
 +
z = linspace(0, H_agua, nz);
 +
[X, Y, Z] = meshgrid(x, y, z);
 +
r2 = X.^2 + Y.^2;
 +
exp_factor = exp(-r2/R^2);
 +
Vx = v0 * exp_factor * cos(phi);
 +
Vy = v0 * exp_factor * sin(phi);
 +
Vz = v1 * sin(pi * Z/H_agua) .* exp_factor;
 +
n = [cos(phi), sin(phi), 0];
 +
V_normal = Vx * n(1) + Vy * n(2) + Vz * n(3);
 +
dx = x(2) - x(1);
 +
dy = y(2) - y(1);
 +
dz = z(2) - z(1);
 +
dV = dx * dy * dz;
 +
Q_cartesiano = sum(V_normal(:) * dV, 'omitnan');
 +
fprintf('Caudal: %.3f m³/s\n', Q_cartesiano);}}
 +
 
 +
 
 +
 
 +
El caudal obtenido mediante integración numérica es:
 +
 
 +
\begin{equation}
 +
Q = 181\,245.7\,\text{m}^3/\text{s}
 +
\end{equation}
 +
 
 +
Este valor matemático representa el volumen de agua que entra al embalse por segundo, confirmando el gran caudal del río Lozoya y la importancia de la presa de El Atazar para el abastecimiento de Madrid.

Revisión actual del 18:13 2 dic 2025


1 Introducción

La presa de El Atazar esta situada sobre el río Lozoya, en la Comunidad de Madrid. Se trata de la presa con mayor capacidad y relevancia de la Comunidad, constituyendo una infraestructura esencial para el suministro de agua potable a Madrid y a toda su área metropolitana.


Fue construida entre 1968 y 1972 y funciona como una presa de doble curvatura, arco-gravedad. Entre sus principales características cabe destacar su altura de 134 metros, una longitud de coronación de 484 metros y profundidades que alcanzan los 100 metros. Su gran capacidad de almacenamiento, que alcanza los 425 hm³ y cubre una superficie de 1.070 hectáreas, convierte al embalse de El Atazar en un elemento fundamental para el abastecimiento, regulación del agua y paisaje de Madrid.


El objetivo de este trabajo es analizar y representar la geometría de la presa con el fin de realizar posteriormente un análisis de la sedimentación del embalse y del flujo de la entrada del rio. Para llevar a cabo nuestro estudio, es de suma importancia dominar conocimientos relacionados con Teoría de Campos y con software de programación y cálculo numérico Matlab.

Trabajo realizado por estudiantes
Título La presa de El Atazar. Grupo 5
Asignatura Teoría de Campos
Curso 2025-26
Autores
  • Nerea Castillero Gómez
  • Elena Álvarez Rodríguez
  • Jakub Rog
  • Daniel Municio Guerrero
  • Juan Martín Rodríguez
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura


1.1 Modelo geométrico de la presa

Para el análisis se considerará la superficie de la presa en el lado aguas arriba, es decir, la zona en contacto directo con el agua. Se supone que la sección transversal puede aproximarse mediante un arco de circunferencia con eje de simetría en el valle, mientras que la sección longitudinal se comporta como un arco parabólico. Trabajaremos durante la ejecución de este trabajo en coordenadas cilíndricas (r, θ, z). La presa 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.

[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 presa y la imagen resultante:

Figura 1. Representación de la presa.


1.2 Representación de la presa

2 Flujo de entrada del río

En esta sección analizaremos el comportamiento del flujo de entrada del río Lozoya en el embalse. Para ello, definiremos un campo de velocidad que modela la entrada del río por extremo norte. Inicialmente, consideraremos un modelo simple en dos dimensiones del campo de velocidad, el cual ampliaremos más tarde a tres dimensiones para incorporar una componente vertical, ofreciendo así una representación más real del flujo de entrada.

El estudio de este campo vectorial se lleva a cabo mediante su visualización y el análisis de sus propiedades diferenciales fundamentales. Específicamente, calcularemos y representaremos su divergencia y su rotacional. Finalmente, y como aplicación práctica, calcularemos el caudal de entrada bajo una condición de flujo, demostrando la utilidad del modelo para la medición de variables hidrológicas clave.

Campo de velocidad del flujo de entrada (2D): [math]\vec{v}(x, y) = v_0 \, e^{-\frac{x^2 + y^2}{R^2}} (\cos \varphi \, \vec{i} + \sin \varphi \, \vec{j})[/math]


Campo de velocidad del flujo de entrada (3D): [math]\vec{v}(x, y, z) = v_0 e^{-\frac{x^2 + y^2}{R^2}} (\cos \varphi \, \vec{i} + \sin \varphi \, \vec{j}) + v_1 \sin\left(\pi \frac{z}{H_{\text{agua}}}\right) e^{-\frac{x^2 + y^2}{R^2}} \vec{k}[/math]

El campo de velocidad del flujo de entrada se modela con los parámetros \( v_0 = 0.5 \, \text{m/s} \) para la velocidad máxima en el centro de entrada, \( v_1 = 0.1 \, \text{m/s} \) para la componente vertical máxima, \( R = 30 \, \text{m} \) como radio característico, y \( \varphi = \frac{\pi}{6} \) como ángulo de entrada. El dominio de estudio comprende las coordenadas horizontales \( (x, y) \in [-100, 100]^2 \, \text{m} \) con la restricción \( \sqrt{x^2 + y^2} \leq \rho_a \), donde \( \rho_a \) es el radio de la presa a la altura \( z = H_{\text{agua}} \), mientras que la coordenada vertical abarca el intervalo \( z \in [0, H_{\text{agua}}] \).

2.1 Representación del campo vectorial de velocidad

En este apartado realizaremos la visualización completa del campo de velocidad asociado al flujo de entrada del río Lozoya en el embalse de El Atazar. Apoyándonos en dos imágenes claves: el plano horizontal de la superficie libre del agua y el plano vertical para así saber la variación con la profundidad.

Representación del campo vectorial de la velocidad
% Parámetros del problema
v0 = 0.5;      % m/s
v1 = 0.1;      % m/s
R = 30;        % m
phi = pi/6;    % ángulo de entrada
H_agua = 125;  % m

% Cálculo de rho_a (radio de la presa a z = H_agua)
H = 134;       % altura total de la presa
rho0 = 150;    % m
b = 40;        % m
% Fórmula: rho = rho0 + b*(1 - z^2/H^2)
rho_a = rho0 + b*(1 - (H_agua^2)/(H^2));

% Plano horizontal (z = H_agua - nivel superficial)
x = linspace(-100, 100, 20);
y = linspace(-100, 100, 20);
[X, Y] = meshgrid(x, y);
Z_superficie = H_agua * ones(size(X));
% Campo de velocidad en superficie (z = H_agua)
r_cuad = X.^2 + Y.^2;
factor_exp = exp(-r_cuad/R^2);
% Componentes del vector velocidad
Vx = v0 * factor_exp * cos(phi);
Vy = v0 * factor_exp * sin(phi);
Vz = v1 * sin(pi * Z_superficie/H_agua) .* factor_exp;
% Filtrar puntos fuera del dominio (dentro de la presa)
dentro_presa = sqrt(X.^2 + Y.^2) <= rho_a;
Vx(~dentro_presa) = NaN;
Vy(~dentro_presa) = NaN;
Vz(~dentro_presa) = NaN;
% Graficar campo vectorial en plano horizontal
figure('Position', [100, 100, 800, 600]);
subplot(2,2,1);
quiver(X, Y, Vx, Vy, 2, 'b');
axis equal;
xlabel('x (m)');
ylabel('y (m)');
title('Campo de velocidad - Plano Horizontal (z = H_{agua})');
grid on;

% Plano vertical (y = 0 - sección longitudinal)
x_vert = linspace(-100, 100, 15);
z_vert = linspace(0, H_agua, 15);
[X_vert, Z_vert] = meshgrid(x_vert, z_vert);
Y_vert = zeros(size(X_vert));
% Campo de velocidad en plano vertical y=0
r_cuad_vert = X_vert.^2 + Y_vert.^2;
factor_exp_vert = exp(-r_cuad_vert/R^2);
Vx_vert = v0 * factor_exp_vert * cos(phi);
Vz_vert = v1 * sin(pi * Z_vert/H_agua) .* factor_exp_vert;
% Filtrar puntos fuera del dominio
dentro_presa_vert = sqrt(X_vert.^2 + Y_vert.^2) <= rho_a;
Vx_vert(~dentro_presa_vert) = NaN;
Vz_vert(~dentro_presa_vert) = NaN;
subplot(2,2,2);
quiver(X_vert, Z_vert, Vx_vert, Vz_vert, 2, 'b');
xlabel('x (m)');
ylabel('z (m)');
title('Campo de velocidad - Plano Vertical (y = 0)');
grid on;

Plano Horizontal: El agua entra formando un chorro que va de noreste a suroeste. En el centro va más rápido (0.5 m/s) y se va haciendo más lento hacia los lados, como un abanico que se abre.

Plano Vertical: En profundidad, el agua se mueve más lento en el fondo y la superficie, y más rápido a media profundidad. También hay movimiento vertical que mezcla el agua.

2.2 La divergencia del campo

En este apartado realizaremos el cálculo de la divergencia del campo de velocidad, para estudiar el comportamiento del flujo. En el apartado físico la divergencia representa la compresión o expansión del fluido (divergencia positiva: partículas se separan / divergencia negativa: partículas se juntan). Representaremos la divergencia en su plano vertical y en su plano horizontal.

Representación plano horizontal de la divergencia.
Representación plano vertical de la divergencia.
clear,clc
% Parámetros
v0 = 0.5; v1 = 0.1; R = 30; phi = pi/6; H_agua = 125;
H = 134; rho0 = 150; b = 40;
rho_a = rho0 + b*(1 - (H_agua^2)/(H^2));
% DIVERGENCIA EN PLANO HORIZONTAL (z = H_agua constante)
x_div = linspace(-80, 80, 40);
y_div = linspace(-80, 80, 40);
[X_div, Y_div] = meshgrid(x_div, y_div);
Z_div = H_agua * ones(size(X_div));
% Campo de velocidad completo
r_cuad = X_div.^2 + Y_div.^2;
factor_exp = exp(-r_cuad/R^2);
Vx = v0 * factor_exp * cos(phi);
Vy = v0 * factor_exp * sin(phi);
Vz = v1 * sin(pi * Z_div/H_agua) .* factor_exp;
% CÁLCULO DE DERIVADAS PARCIALES
dx = x_div(2) - x_div(1);
dy = y_div(2) - y_div(1); 
% Derivadas usando gradient
[dVx_dx, dVx_dy] = gradient(Vx, dx, dy);
[dVy_dx, dVy_dy] = gradient(Vy, dx, dy);
[dVz_dx, dVz_dy] = gradient(Vz, dx, dy);
% Como z es constante en este plano, ∂Vz/∂z = 0
divergencia_horizontal = dVx_dx + dVy_dy;
% Filtrar puntos fuera del dominio
dentro_presa = sqrt(X_div.^2 + Y_div.^2) <= rho_a;
divergencia_horizontal(~dentro_presa) = NaN;
Vx(~dentro_presa) = NaN;
Vy(~dentro_presa) = NaN;
figure('Position', [100, 100, 1200, 500]);
% plot 1: Mapa de colores de divergencia
plot(1);
pcolor(X_div, Y_div, divergencia_horizontal);
shading interp;
colorbar;
colormap(jet);
hold on;
% Contornos de divergencia
contour(X_div, Y_div, divergencia_horizontal, 15, 'k', 'LineWidth', 0.5);
x_vert = linspace(-80, 80, 30);
z_vert = linspace(0, H_agua, 25);
[X_vert, Z_vert] = meshgrid(x_vert, z_vert);
Y_vert = zeros(size(X_vert));
% Campo de velocidad en plano vertical
r_cuad_vert = X_vert.^2 + Y_vert.^2;
factor_exp_vert = exp(-r_cuad_vert/R^2);
Vx_vert = v0 * factor_exp_vert * cos(phi);
Vy_vert = v0 * factor_exp_vert * sin(phi);
Vz_vert = v1 * sin(pi * Z_vert/H_agua) .* factor_exp_vert;
% Derivadas en plano vertical
dx_vert = x_vert(2) - x_vert(1);
dz_vert = z_vert(2) - z_vert(1);
[dVx_dx_vert, dVx_dz_vert] = gradient(Vx_vert, dx_vert, dz_vert);
[dVy_dx_vert, dVy_dz_vert] = gradient(Vy_vert, dx_vert, dz_vert);
[dVz_dx_vert, dVz_dz_vert] = gradient(Vz_vert, dx_vert, dz_vert);
divergencia_vertical = dVx_dx_vert + dVz_dz_vert;
% Filtrar puntos fuera del dominio
dentro_presa_vert = sqrt(X_vert.^2 + Y_vert.^2) <= rho_a;
divergencia_vertical(~dentro_presa_vert) = NaN;
figure('Position', [100, 600, 1000, 400]);
plot(2);
pcolor(X_vert, Z_vert, divergencia_vertical);
shading interp;
colorbar;
hold on;
contour(X_vert, Z_vert, divergencia_vertical, 10, 'k', 'LineWidth', 0.5);
xlabel('x (m)'); ylabel('z (m)');
title('Divergencia ∇·v (s^{-1}) - Plano Vertical (y=0)');

Después de analizar los resultados identificamos que el agua se comprime en la dirección de la velocidad (divergencia negativa) y se expande hacia los lados (divergencia positiva). Se puede concluir que se crea una zona de mayor presión justo en la zona donde el río se encuentra con el embalse.

2.3 El rotacionaldel campo

Este apartado analizaremos el rotacional del campo de velocidad, que mide cuánto gira el agua al entrar al embalse. Calcularemos si se forman remolinos o torbellinos, lo que afecta directamente cómo se mezcla el agua, se reparten los sedimentos y se oxigena el embalse.

Representación plano horizontal del rotacional.
Representación plano vertical del rotacional.
clear,clc
% Parámetros
v0 = 0.5; v1 = 0.1; R = 30; phi = pi/6; H_agua = 125;
H = 134; rho0 = 150; b = 40;
rho_a = rho0 + b*(1 - (H_agua^2)/(H^2));
%ROTACIONAL EN PLANO HORIZONTAL (z = H_agua constante)
x_rot = linspace(-80, 80, 35);
y_rot = linspace(-80, 80, 35);
[X_rot, Y_rot] = meshgrid(x_rot, y_rot);
Z_rot = H_agua * ones(size(X_rot));
% Campo de velocidad
r_cuad = X_rot.^2 + Y_rot.^2;
factor_exp = exp(-r_cuad/R^2);
Vx = v0 * factor_exp * cos(phi);
Vy = v0 * factor_exp * sin(phi);
Vz = v1 * sin(pi * Z_rot/H_agua) .* factor_exp;
% Derivadas parciales
dx = x_rot(2) - x_rot(1);
dy = y_rot(2) - y_rot(1);
[dVx_dx, dVx_dy] = gradient(Vx, dx, dy);
[dVy_dx, dVy_dy] = gradient(Vy, dx, dy);
[dVz_dx, dVz_dy] = gradient(Vz, dx, dy);
% ROTACIONAL 3D: ∇×v = (∂Vz/∂y - ∂Vy/∂z, ∂Vx/∂z - ∂Vz/∂x, ∂Vy/∂x - ∂Vx/∂y)
rot_x = dVz_dy - 0;  % ∂Vy/∂z ≈ 0 en z constante
rot_y = 0 - dVz_dx;  % ∂Vx/∂z ≈ 0 en z constante
rot_z = dVy_dx - dVx_dy;  % Componente vertical del rotacional
% Filtrar puntos fuera del dominio
dentro_presa = sqrt(X_rot.^2 + Y_rot.^2) <= rho_a;
rot_x(~dentro_presa) = NaN;
rot_y(~dentro_presa) = NaN;
rot_z(~dentro_presa) = NaN;
Vx(~dentro_presa) = NaN;
Vy(~dentro_presa) = NaN;
% REPRESENTACIÓN GRÁFICA - PLANO HORIZONTAL
plot(1);
rot_horiz_mag = sqrt(rot_x.^2 + rot_y.^2);
pcolor(X_rot, Y_rot, rot_horiz_mag);
shading interp;
colorbar;
hold on;
quiver(X_rot, Y_rot, rot_x, rot_y, 2, 'w', 'LineWidth', 1.5);
xlabel('x (m)'); ylabel('y (m)');
title('Rotacional Horizontal


Analizando los rotacionales se puede concluir que existe vorticidad en flujo. Como el agua no entra toda igual de rápido ni en la misma dirección, se generan pequeños remolinos naturales que ayudan a mezclarla. Los cálculos muestran que el rotacional no es cero, lo que indica que el agua presenta movimiento rotacional al entrar al embalse. Esta rotación favorece la formación de remolinos y torbellinos que ayudan a mezclar el agua del río con la del embalse, mejorando la oxigenación en la zona de entrada.


2.4 Cálculo del caudal de entrada del río

En este apartado calcularemos el caudal de entrada del río Lozoya al embalse, es decir, la cantidad de agua que llega por segundo. Simplificaremos el cálculo eliminando la componente vertical (v1 = 0 m/s) para centrarnos en el flujo horizontal principal. A través de integración matemática del campo de velocidad sobre una sección perpendicular a la dirección de entrada, determinaremos cuántos metros cúbicos por segundo aporta el río al embalse, un dato esencial para la gestión hídrica y el balance de agua de toda la presa. 

clear; clc;
% Parámetros
v0 = 0.5;          % Velocidad máxima [m/s]
v1 = 0;            % Velocidad vertical = 0 [m/s]
R = 30;            % Radio característico [m]
phi = pi/6;        % Ángulo de entrada [rad]
H_agua = 125;      % Altura del agua [m]
H = 134;           % Altura total presa [m]
rho0 = 150;        % Radio en coronación [m]
b = 40;            % Parámetro de curvatura [m]
rho_a = rho0 + b*(1 - (H_agua^2)/(H^2));
nx = 80; ny = 80; nz = 40;
x = linspace(-rho_a, rho_a, nx);
y = linspace(-rho_a, rho_a, ny);
z = linspace(0, H_agua, nz);
[X, Y, Z] = meshgrid(x, y, z);
r2 = X.^2 + Y.^2;
exp_factor = exp(-r2/R^2);
Vx = v0 * exp_factor * cos(phi);
Vy = v0 * exp_factor * sin(phi);
Vz = v1 * sin(pi * Z/H_agua) .* exp_factor;
n = [cos(phi), sin(phi), 0];
V_normal = Vx * n(1) + Vy * n(2) + Vz * n(3);
dx = x(2) - x(1);
dy = y(2) - y(1);
dz = z(2) - z(1);
dV = dx * dy * dz;
Q_cartesiano = sum(V_normal(:) * dV, 'omitnan');
fprintf('Caudal: %.3f m³/s\n', Q_cartesiano);



El caudal obtenido mediante integración numérica es:

\begin{equation} Q = 181\,245.7\,\text{m}^3/\text{s} \end{equation}

Este valor matemático representa el volumen de agua que entra al embalse por segundo, confirmando el gran caudal del río Lozoya y la importancia de la presa de El Atazar para el abastecimiento de Madrid.

Obtenido de «https://mat.caminos.upm.es/w/index.php?title=La_presa_de_Atazar_NEMJJ&oldid=94849»