Diferencia entre revisiones de «La presa de El Atazar NEMJJ»

De MateWiki
Saltar a: navegación, buscar
(→‎Modelo de sedimentación en el fondo del embalse)
(→‎Cálculo de la masa, volumen y porcentaje de los sedimentos)
 
(No se muestran 20 ediciones intermedias de 5 usuarios)
Línea 284: Línea 284:
 
[[Archivo:gradientedefinitivo.png|600px|miniaturadeimagen|Mapa de sedimentos en el fondo del embalse]]
 
[[Archivo:gradientedefinitivo.png|600px|miniaturadeimagen|Mapa de sedimentos en el fondo del embalse]]
 
{{matlab|codigo=
 
{{matlab|codigo=
−
## Parámetros
+
% Parámetros
 
S0 = 50;
 
S0 = 50;
 
alpha = 3;
 
alpha = 3;
 
L = 500;
 
L = 500;
  
−
## Mallado circular
+
% Mallado circular
 
rmax = 1000;
 
rmax = 1000;
 
N = 80;
 
N = 80;
Línea 299: Línea 299:
 
mask = R <= rmax;
 
mask = R <= rmax;
  
−
## Sedimentación
+
% Sedimentación
 
S = S0 * (1 + alpha * exp(-(X.^2 + Y.^2)/L^2));
 
S = S0 * (1 + alpha * exp(-(X.^2 + Y.^2)/L^2));
 
S(~mask) = NaN;
 
S(~mask) = NaN;
  
−
## Gradiente
+
% Gradiente
 
[dSdx, dSdy] = gradient(S, x, y);
 
[dSdx, dSdy] = gradient(S, x, y);
  
Línea 309: Línea 309:
 
dSdy(~mask) = NaN;
 
dSdy(~mask) = NaN;
  
−
## Gráfico
+
% Gráfico
 
figure(1); clf;
 
figure(1); clf;
 
hold on;
 
hold on;
  
−
## Heatmap
+
% Heatmap
 
pcolor(X, Y, S);
 
pcolor(X, Y, S);
 
shading interp;
 
shading interp;
Línea 319: Línea 319:
 
colorbar;
 
colorbar;
  
−
## Campo vectorial
+
% Campo vectorial
 
quiver(X, Y, dSdx, dSdy, "k");
 
quiver(X, Y, dSdx, dSdy, "k");
  
Línea 328: Línea 328:
  
 
hold off;}}
 
hold off;}}
 +
 +
<math>
 +
\nabla S(x,y)
 +
= \left(
 +
-\frac{2S_{0}\alpha x}{L^{2}}\, e^{-\frac{x^{2}+y^{2}}{L^{2}}},
 +
\;
 +
-\frac{2S_{0}\alpha y}{L^{2}}\, e^{-\frac{x^{2}+y^{2}}{L^{2}}}
 +
\right)
 +
= -\frac{2S_{0}\alpha}{L^{2}}\,e^{-\frac{x^{2}+y^{2}}{L^{2}}}\,(x,y)
 +
</math>
 +
 +
Se puede ver que el gradiente apunta hacia el centro.
 +
 +
=== Cálculo de la masa, volumen y porcentaje de los sedimentos ===
 +
 +
A la hora de afrontar este apartado, calcularemos <math>S(x,y)</math> mediante una doble integral sobre la superficie total del embalse. Además, para facilitar los cálculos, consideramos que <math>z = 0</math> es el fondo y que el área es el sector circular de la base de la presa.
 +
El radio de la base se calcula de la siguiente manera: <math> \rho_\text{base} = \rho_0 + b = 150 + 40 = 190~\text{m}. </math>, siendo 150 el radio de coronación y 40 el parámetro de curvatura.
 +
 +
Al aplicar la igualdad <math>\rho^2 = x^2 + y^2</math>, la fórmula para la concentración de sedimentos (<math>S(x,y)</math>) se convierte en <math> S(x,y) = S_0 \left( 1 + \alpha e^{-\rho^2/L^2} \right), </math>. El ángulo horizontal esta determinado en el intervalo <math>\theta \in [2\pi/3, 4\pi/3]</math>.
 +
 +
La masa total de sedimentos se obtiene integrando la concentración sobre toda la superficie: <math> M = \iint_{A} S(x,y), dA </math>.
 +
Debido a la simetría del sector circular, la integral se puede reescribir en coordenadas polares <math>( x = r\cos\theta, y = r\sin\theta, dA = r)</math>:
 +
 +
<math> M = \int_{\theta_0}^{\theta_1} \int_0^{\rho_\text{base}} S(r), r , dr , d\theta = \int_{\theta_0}^{\theta_1} \int_0^{\rho_\text{base}} S_0 \left( 1 + \alpha e^{-r^2/L^2} \right) r </math>
 +
 +
Si resolvemos la integral la integral con respecto <math>r</math>:
 +
<math> \int_0^{\rho_\text{base}} r, dr = \frac{\rho_\text{base}^2}{2}, \quad \int_0^{\rho_\text{base}} r, e^{-r^2/L^2} dr = \frac{L^2}{2} \left( 1 - e^{-\rho_\text{base}^2/L^2} \right) </math>
 +
 +
De esta manera, la masa total queda:
 +
<math> M = S_0 \Delta \theta \left[ \frac{\rho_\text{base}^2}{2} + \frac{\alpha L^2}{2} \left( 1 - e^{-\rho_\text{base}^2/L^2} \right) \right] </math>
 +
 +
Si sustituimos por los valores dados (<math>S_0 = 50~\text{kg/m²}, \alpha = 3, L = 500~\text{m}, \rho_\text{base} = 190~\text{m}, \Delta\theta = 2\pi/3</math>):
 +
<math> M \approx 7,17 \times 10^6~\text{kg}. </math>
 +
 +
Como la densidad de los sedimentos es <math>\rho_s = 1500~\text{kg/m³}</math>, y sabemos la masa total, el volumen que ocupan es:
 +
<math> V = \frac{M}{\rho_s} \approx \frac{7,17 \times 10^6}{1500} =4780~\text{m³}. </math>
 +
 +
Por último, determinamos el porcentaje con respecto a la capacidad total del embalse (<math>V_{\text{emb}} = 425~\text{hm³} = 425 \times 10^6~\text{m³}</math>).
 +
 +
 +
<math> \text{Porcentaje} = \frac{V}{V_{\text{emb}}} \cdot 100 \approx 0,0011</math>%
 +
 +
En conclusión, vemos que el porcentaje es demasiado pequeño en comparación con la capacidad embalse total como para ser preocupante, ya que la sedimentación aun no afecta al volumen posible.
 +
 +
=== Curvas de isoconcentración ===
 +
 +
[[Archivo:isoconcentracion.png|600px|miniaturadeimagen|Mapa de sedimentos en el fondo del embalse]]
 +
{{matlab|codigo=
 +
%  Parámetros
 +
S0 = 50;
 +
alpha = 3;
 +
L = 500;
 +
R = 1500;
 +
 +
%  Mallado y máscara
 +
n = 500;
 +
x = linspace(-R, R, n);
 +
y = linspace(-R, R, n);
 +
[X, Y] = meshgrid(x, y);
 +
inside = (X.^2 + Y.^2) <= R^2;
 +
 +
%  Campo
 +
S = S0 * (1 + alpha * exp(-(X.^2 + Y.^2)/L^2));
 +
S(~inside) = NaN;
 +
 +
%  Curvas de isoconcentración
 +
figure;
 +
nlevels = 12;  % número de curvas
 +
[cc, h] = contour(X, Y, S, nlevels, 'LineWidth', 1.4);
 +
clabel(cc, h, 'FontSize', 8);
 +
 +
axis equal tight;
 +
xlabel('x (m)');
 +
ylabel('y (m)');
 +
title('Curvas de isoconcentración');
 +
colormap(jet);
 +
cb = colorbar;
 +
cb.Label.String = 'S (kg/m^2)';
 +
 +
hold on;
 +
t = linspace(0, 2*pi, 400);
 +
plot(R*cos(t), R*sin(t), 'k-', 'LineWidth', 1);
 +
hold off;}}
 +
 +
Se puede observar que las curvas tienen forma aproximadamente circular, más concentradas cerca del centro.
 +
 +
 +
==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.
 +
  
 
== Análisis de tipos de presas - Grandes presas españolas==
 
== Análisis de tipos de presas - Grandes presas españolas==
Línea 371: Línea 766:
  
 
===Comparación tipos de presas===
 
===Comparación tipos de presas===
−
A continuación, llevaremos a cabo una pequeña comparación de los tipos comentados previamente
+
A continuación, expondré lo comentado anteriormente en formato de tabla, para poder apreciar sus cualidades de una manera mas visual e intuitiva.  
−
 
+
−
Estabilidad: gravedad > arco. El arco es muy estable si las paredes laterales son buenas.
+
−
 
+
−
Uso de materiales: arco < gravedad.
+
−
 
+
−
Condiciones geológicas: arco necesita garganta y roca firme - gravedad necesita cimiento sólido y capaz de soportar mucho peso - presas de tierra admiten cimientos más      débiles.
+
  
 
=== Grandes presas españolas — Tipología ===
 
=== Grandes presas españolas — Tipología ===

Revisión actual del 18:58 4 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.

right


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


2 Modelo geométrico de la presa

Para el análisis se considerará la superficie de la presa en su cara aguas arriba, zona del río antes de la presa y en contacto directo con el agua. La estructura de la presa consta de una sección transversal constituida por un arco de circunferencia en la vista horizontal, con eje de simetría en el valle, y de una sección longitudinal constituida por un arco parabólico en la sección vertical.

En coordenadas cilíndricas (𝜌,𝜃,𝑧), con los parámetros definidos [math]θ ∈ [\frac{3π}{4}, \frac{5π}{4}][/math] y [math]Z ∈ [0,H][/math] y con referencias con origen en el valle y eje 𝑧 apuntando hacia arriba. Ecuación del modelo de la superficie de la presa:

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

3 Cáculos y representaciones obligatorias

3.1 - Representación de la superficie parametrizada de la presa en su cara aguas arriba

A continuación, expondré la representación de la superficie. Para ello, usaremos comandos como "meshgrid()", el cual construye una malla que permita parametrizar la superficie en coordenadas cilíndricas la ecuación; "surf()", el cual plasmará la gráfica de la superficie. Además de dichos comandos, es importante dominar en este apartado el manejo de parametrización de coordenadas cartesianas.

Código empleado para la obtención de la presa y la imagen resultante:

Figura 2. Representación de la presa.
% Sean los datos de la presa
h = 134;      % Altura (en metros)
rho0 = 150;   % Radio (en metros)
b = 40;       % Parámetro de curvatura
% Definimos 0 a h
z = linspace(0, h, 50);
theta = linspace(2*pi/3, 4*pi/3, 50);
% Creación de la malla
[Z, TH] = meshgrid(z, theta);
% Calculamos Rho:
Rho = rho0 + b * (1 - (Z.^2) / (h^2));
% Convertimos de cilíndricas a (X, Y, Z)
X = Rho .* cos(TH);
Y = Rho .* sin(TH);
figure;
surf(X, Y, Z);
% Estética del gráfico. Sea Z la altura.
title('Representación de la Presa de El Atazar (Cara Aguas Arriba)');
xlabel('X metros');
ylabel('Y metros');
zlabel('Z metros');
axis equal;           
shading interp;      
colormap parula;      
colorbar;
grid on;



3.2 - Presiones sobre la presa - Campo escalar de presión

A continuación, analizaremos y representaremos 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. Trataremos el agua de la presa como un fluido estático, generando el peso del agua una presión a su alrededor, que actúa en las paredes de la presa.

El campo escalar de presiones viene dado por la función:

[math]P(z)=ρgh(z)[/math]
.

Sea: ρ, la densidad del agua; g, la aceleración de la gravedad; y h(z), la profundidad del agua.

A continuación, expondré la representación del campo para visualizar cómo varía la presión a lo largo de la superficie aguas arriba de la presa, permitiéndonos ello identificar las zonas de mayor y menor presión. Para ello, usaremos comandos como "surf()", el cual plasmará la gráfica de la superficie; "colomap", que diferenciara los colores en función de la presión que varía con la profundidad P(z); y "colorbar", para visualizar el rango de presiones. Además de dichos comandos, es importante dominar en este apartado el manejo de parametrización de coordenadas cartesianas.

Código empleado para la obtención de la presa y la imagen resultante:

Figura 3. Representación de las presiones de la presa.
% Sean los datos de la presa
h = 134;      % Altura (en metros)
rho0 = 150;   % Radio (en metros)
b = 40;       % Parámetro de curvatura
% Definimos 0 a h
z = linspace(0, h, 50);
theta = linspace(2*pi/3, 4*pi/3, 50);
% Creación de la malla
[Z, TH] = meshgrid(z, theta);
% Calculamos Rho:
Rho = rho0 + b * (1 - (Z.^2) / (h^2));
% Convertimos de cilíndricas a (X, Y, Z)
X = Rho .* cos(TH);
Y = Rho .* sin(TH);
% 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 metros');
ylabel('Y metros');
zlabel('Z metros');
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;



La imagen representa la variación de presión en función de la coordenada Z, la altura. Se puede apreciar dos gamas de colores: gama de tonos fríos, que representa zonas de menor presión; y gama de tonos cálidos, que representa las regiones de mayor presión. Los valores más altos de presión se localicen en la base de la presa y disminuyan progresivamente con la altura.

3.3 Fuerza de presión total y presión por unidad de superficie

A continuación, calcularemos la fuerza de presión total y la presión por unidad de superficie para dos valores distintos del parámetro de curvatura: 𝑏 = 40 m, caso de doble curvatura; y 𝑏 = 0 m, caso de curvatura simple. Por ultimo, deberemos deducir cuál de las dos configuraciones soporta mejor la presión.

Datos empleados para la obtención de la fuerza de presión total y la presión por unidad de superficie. 

% Sean los datos generales de la presa
h = 134;      % Altura (en metros)
rho0 = 150;   % Radio (en metros)
% Parámetros físicos
rho = 1000;  % Densidad del agua (kg/m^3)
g = 9.81;    % Aceleración de la gravedad (m/s^2)         
r0 = 968/pi;       % Radio en la altura máxima (m)



3.3.1 Opción 1 - Caso de doble curvatura

Código empleado para la obtención de la fuerza de presión total y la presión por unidad de superficie en el caso de doble curvatura y su representación.

Figura 4. Fuerza de presión total.
% OPCION 1 - CASO DE DOBLE CURVATURA
%Parametro de curvatura:
b = 40
% Coordenadas cilíndricas
theta = linspace(3*pi/4, 5*pi/4, 20);   % Ángulo θ
z = linspace(0, h, 20);                % Sea z la altura 
[Theta, Z] = meshgrid(theta, z);       
R = r0 + b * (1 - Z.^2 / h^2);         % Radio en función de la alturaz
% Coordenadas cartesianas
X = R .* cos(Theta);                   
Y = R .* sin(Theta);                  
% Derivadas parciales para las normales
dr_dz = -2 * b * Z / h^2;              % Derivada parcial de r respecto a z
n_r = 1 ./ sqrt(1 + dr_dz.^2);         % Componente radial del vector normal
n_x = n_r .* cos(Theta);               % Proyección de la normal en X
n_y = n_r .* sin(Theta);               % Proyección de la normal en Y
% Campo de presión
P = rho * g * (h - Z);                 % Presión en función de la profundidad
% Vectores de fuerza de presión
F_x = -P .* n_x;                       % En X
F_y = -P .* n_y;                       % En Y
% Representación gráfica de la presa
figure;
q = quiver3(X, Y, Z, F_x, F_y, zeros(size(F_x)), 1, 'Color', 'm');
q.MaxHeadSize = 2; 
title('Fuerza de presión en la presa (3D)');
xlabel('X metros');
ylabel('Y metros');
zlabel('Z metros');
grid on;
axis equal;



3.3.2 Opción 2 - Caso de curvatura simple

Código empleado para la obtención de la fuerza de presión total y la presión por unidad de superficie en el caso de curvatura simple y su representación.

Figura 4. Fuerza de presión total.
% OPCION 2 - CASO DE CURVATURA SIMPLE
%Parametro de curvatura:
b = 0
% Coordenadas cilíndricas
theta = linspace(3*pi/4, 5*pi/4, 20);   % Ángulo θ
z = linspace(0, h, 20);                % Sea z la altura 
[Theta, Z] = meshgrid(theta, z);       
R = r0 + b * (1 - Z.^2 / h^2);         % Radio en función de la alturaz
% Coordenadas cartesianas
X = R .* cos(Theta);                   
Y = R .* sin(Theta);                  
% Derivadas parciales para las normales
dr_dz = -2 * b * Z / h^2;              % Derivada parcial de r respecto a z
n_r = 1 ./ sqrt(1 + dr_dz.^2);         % Componente radial del vector normal
n_x = n_r .* cos(Theta);               % Proyección de la normal en X
n_y = n_r .* sin(Theta);               % Proyección de la normal en Y
% Campo de presión
P = rho * g * (h - Z);                 % Presión en función de la profundidad
% Vectores de fuerza de presión
F_x = -P .* n_x;                       % En X
F_y = -P .* n_y;                       % En Y
% Representación gráfica de la presa
figure;
q = quiver3(X, Y, Z, F_x, F_y, zeros(size(F_x)), 1, 'Color', 'g');
q.MaxHeadSize = 2; 
title('Fuerza de presión en la presa (3D)');
xlabel('X metros');
ylabel('Y metros');
zlabel('Z metros');
grid on;
axis equal;



La configuración que mejor soporta la presión se trata de la presa de doble curvatura (b = 40). Esto se explica debido a que la presión hidrostática se reparte hacia los lados, siendo estos estribos o montañas. La forma de arco de la presa, tanto horizontal, como vertical, hace que la curva empuje la fuerza horizontalmente. Es la opción estructuralmente más eficiente.

La presa de curvatura simple (b = 0), no existe variación de radio con la altura y la presión se transfiere más directamente al cuerpo de la presa.

4 Modelo de sedimentación en el fondo del embalse

Los ríos transportan sedimentos que se depositan en el fondo del embalse, reduciendo gradualmente su capacidad. Modelamos la concentración de sedimentos depositados en el fondo del embalse (en kg/m²) como:

[math]\displaystyle S(x,y)=S_{0}\left(1+\alpha\,e^{-\frac{x^{2}+y^{2}}{L^{2}}}\right)[/math]

donde:

  • [math]S_{0}=50\ \mathrm{kg}\,\mathrm{m}^{-2}[/math] es la sedimentación base,
  • [math]\alpha=3[/math] modela la mayor acumulación cerca de la entrada del río,
  • [math]L=500\ \mathrm{m}[/math] es una escala característica.

Para simplificar, consideramos el fondo del embalse aproximadamente plano en [math]z=0[/math] y con forma de sector circular centrado en el valle, delimitado por el arco de la presa en su base.

4.1 Mapa de sedimentos

A continuación se muestra un mapa de la distribución de sedimentos en el fondo del embalse, generado mediante una simulación en Octave:

Mapa de sedimentos en el fondo del embalse
%  Parámetros del problema 
S0 = 50;      
alpha = 3;
L = 500;      
R = 1500;     

%  Mallado cartesiano (cuadrado) y máscara circular 
n = 500;                                     
x = linspace(-R, R, n);
y = linspace(-R, R, n);
[X, Y] = meshgrid(x, y);
inside = (X.^2 + Y.^2) <= R^2;                

%  Campo S(x,y) y enmascarar el exterior 
S = S0 * (1 + alpha * exp(-(X.^2 + Y.^2) / L^2));
S(~inside) = NaN;                             

%  Gradiente analítico (dS/dx, dS/dy) y enmascarado 
dSdx = S0 * alpha * exp(-(X.^2 + Y.^2)/L^2) .* (-2 .* X / L^2);
dSdy = S0 * alpha * exp(-(X.^2 + Y.^2)/L^2) .* (-2 .* Y / L^2);
dSdx(~inside) = NaN;
dSdy(~inside) = NaN;

%  Figura con 3 subplots: 1 y 2 arriba, 3 ocupando toda la fila inferior 
figure('Name','Sedimentación - mapas y campos','Units','normalized','Position',[0.05 0.05 0.9 0.8]);

% Mapa de color 
pcolor(X, Y, S);
shading interp;
axis equal tight;
xlabel('x (m)');
ylabel('y (m)');
title('sedimentación en el fondo','FontWeight','bold','FontSize',14);   
colormap(gca, jet);
cb1 = colorbar;
cb1.Label.String = 'S (kg/m^2)';
set(gca,'XLim',[-R R],'YLim',[-R R]);


4.2 Mapa del gradiente de los sedimentos

Mapa de sedimentos en el fondo del embalse
% Parámetros
S0 = 50;
alpha = 3;
L = 500;

% Mallado circular
rmax = 1000;
N = 80;
x = linspace(-rmax, rmax, N);
y = linspace(-rmax, rmax, N);
[X,Y] = meshgrid(x,y);

R = sqrt(X.^2 + Y.^2);
mask = R <= rmax;

% Sedimentación
S = S0 * (1 + alpha * exp(-(X.^2 + Y.^2)/L^2));
S(~mask) = NaN;

% Gradiente
[dSdx, dSdy] = gradient(S, x, y);

dSdx(~mask) = NaN;
dSdy(~mask) = NaN;

% Gráfico
figure(1); clf;
hold on;

% Heatmap
pcolor(X, Y, S);
shading interp;
colormap(turbo);
colorbar;

% Campo vectorial
quiver(X, Y, dSdx, dSdy, "k");

axis equal;
xlabel("x (m)");
ylabel("y (m)");
title("Gradiente de sedimentación y campo vectorial");

hold off;


[math] \nabla S(x,y) = \left( -\frac{2S_{0}\alpha x}{L^{2}}\, e^{-\frac{x^{2}+y^{2}}{L^{2}}}, \; -\frac{2S_{0}\alpha y}{L^{2}}\, e^{-\frac{x^{2}+y^{2}}{L^{2}}} \right) = -\frac{2S_{0}\alpha}{L^{2}}\,e^{-\frac{x^{2}+y^{2}}{L^{2}}}\,(x,y) [/math]

Se puede ver que el gradiente apunta hacia el centro.

4.3 Cálculo de la masa, volumen y porcentaje de los sedimentos

A la hora de afrontar este apartado, calcularemos [math]S(x,y)[/math] mediante una doble integral sobre la superficie total del embalse. Además, para facilitar los cálculos, consideramos que [math]z = 0[/math] es el fondo y que el área es el sector circular de la base de la presa. El radio de la base se calcula de la siguiente manera: [math] \rho_\text{base} = \rho_0 + b = 150 + 40 = 190~\text{m}. [/math], siendo 150 el radio de coronación y 40 el parámetro de curvatura.

Al aplicar la igualdad [math]\rho^2 = x^2 + y^2[/math], la fórmula para la concentración de sedimentos ([math]S(x,y)[/math]) se convierte en [math] S(x,y) = S_0 \left( 1 + \alpha e^{-\rho^2/L^2} \right), [/math]. El ángulo horizontal esta determinado en el intervalo [math]\theta \in [2\pi/3, 4\pi/3][/math].

La masa total de sedimentos se obtiene integrando la concentración sobre toda la superficie: [math] M = \iint_{A} S(x,y), dA [/math]. Debido a la simetría del sector circular, la integral se puede reescribir en coordenadas polares [math]( x = r\cos\theta, y = r\sin\theta, dA = r)[/math]:

[math] M = \int_{\theta_0}^{\theta_1} \int_0^{\rho_\text{base}} S(r), r , dr , d\theta = \int_{\theta_0}^{\theta_1} \int_0^{\rho_\text{base}} S_0 \left( 1 + \alpha e^{-r^2/L^2} \right) r [/math]

Si resolvemos la integral la integral con respecto [math]r[/math]: [math] \int_0^{\rho_\text{base}} r, dr = \frac{\rho_\text{base}^2}{2}, \quad \int_0^{\rho_\text{base}} r, e^{-r^2/L^2} dr = \frac{L^2}{2} \left( 1 - e^{-\rho_\text{base}^2/L^2} \right) [/math]

De esta manera, la masa total queda: [math] M = S_0 \Delta \theta \left[ \frac{\rho_\text{base}^2}{2} + \frac{\alpha L^2}{2} \left( 1 - e^{-\rho_\text{base}^2/L^2} \right) \right] [/math]

Si sustituimos por los valores dados ([math]S_0 = 50~\text{kg/m²}, \alpha = 3, L = 500~\text{m}, \rho_\text{base} = 190~\text{m}, \Delta\theta = 2\pi/3[/math]): [math] M \approx 7,17 \times 10^6~\text{kg}. [/math]

Como la densidad de los sedimentos es [math]\rho_s = 1500~\text{kg/m³}[/math], y sabemos la masa total, el volumen que ocupan es: [math] V = \frac{M}{\rho_s} \approx \frac{7,17 \times 10^6}{1500} =4780~\text{m³}. [/math]

Por último, determinamos el porcentaje con respecto a la capacidad total del embalse ([math]V_{\text{emb}} = 425~\text{hm³} = 425 \times 10^6~\text{m³}[/math]).


[math] \text{Porcentaje} = \frac{V}{V_{\text{emb}}} \cdot 100 \approx 0,0011[/math]%

En conclusión, vemos que el porcentaje es demasiado pequeño en comparación con la capacidad embalse total como para ser preocupante, ya que la sedimentación aun no afecta al volumen posible.

4.4 Curvas de isoconcentración

Mapa de sedimentos en el fondo del embalse
%  Parámetros 
S0 = 50;
alpha = 3;
L = 500;
R = 1500;

%  Mallado y máscara 
n = 500;
x = linspace(-R, R, n);
y = linspace(-R, R, n);
[X, Y] = meshgrid(x, y);
inside = (X.^2 + Y.^2) <= R^2;

%  Campo 
S = S0 * (1 + alpha * exp(-(X.^2 + Y.^2)/L^2));
S(~inside) = NaN;

%  Curvas de isoconcentración 
figure;
nlevels = 12;  % número de curvas
[cc, h] = contour(X, Y, S, nlevels, 'LineWidth', 1.4);
clabel(cc, h, 'FontSize', 8);

axis equal tight;
xlabel('x (m)');
ylabel('y (m)');
title('Curvas de isoconcentración');
colormap(jet);
cb = colorbar;
cb.Label.String = 'S (kg/m^2)';

hold on;
t = linspace(0, 2*pi, 400);
plot(R*cos(t), R*sin(t), 'k-', 'LineWidth', 1);
hold off;


Se puede observar que las curvas tienen forma aproximadamente circular, más concentradas cerca del centro.


5 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}}] \).

5.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.

5.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.

5.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.


5.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.


6 Análisis de tipos de presas - Grandes presas españolas

6.1 Tipos de presas

6.1.1 Presa de gravedad

Consiste en un muro grueso de piedra o hormigón que soporta el agua por su propio peso.

Entre sus ventajas se halla que es muy estable, tiene un diseño simple y resiste bien empujes grandes.

Sus inconvenientes consisten en la necesidad de mucho material, su alto coste y la necesidad de cimientos muy firmes.

Su geología debería estar compuesta por roca buena y continua donde apoyar mucho peso. Se suelen hacer en lugares anchos o aquellos donde la roca no es adecuada para recibir fuera.

6.1.2 Presa de arco o bóveda

Se caracteriza por su forma curva hacia aguas abajo. Dicha curva hace que el empuje del agua se transmita a las paredes laterales, los estribos.

Sus ventajas son su alta eficiencia, el menor uso de hormigón y es apta para paisajes con cañones y gargantas estrechas.

Sus inconvenientes están marcados por la exigencia de tener paredes laterales muy resistentes y por tener cálculos y una construcción mas compleja.

Su geología debería estar compuesta por una garganta o valle con roca sólida y paredes próximas.

6.1.3 Presa arco-gravedad

Se trata de una fusión entre la presa de arco y la presa de gravedad. Contiene una forma curva pero también con suficiente peso para apoyar.

Sus ventajas e inconvenientes combina beneficios de ambos, sumando una mayor versatilidad si la garganta no es perfecta y un consumo medio de hormigón.

6.1.4 Presa de contrafuertes - buttress

Consiste en una pantalla delgada aguas arriba sostenida por contrafuertes en el lado downstream, parecido a pilares, en vez de un bloque macizo.

Entre sus ventajas se halla la necesidad de menos material que una gravedad maciza y el económico coste de material para encofrados y estructura.

Sus inconvenientes se resumen en la alta complejidad de construcción y mantenimiento y la necesidad de cimientos firmes.

6.1.5 Presa de tierra o de relleno - embalse/diques de materiales sueltos

Se trata de un montículo de tierra y roca compactada con núcleo impermeable, sea arcilla o geomembrana.

Sus ventajas son el uso de materiales locales en construcción, siendo ello más barato, y se adapta bien a cimientos menos sólidos.

Los principales inconvenientes son la necesidad de un control de filtraciones y drenaje, baja altura y es poco estética.

Geología ideal: zonas donde no hay roca superficial pero sí materiales de relleno y espacio para la base amplia.

6.2 Comparación tipos de presas

A continuación, expondré lo comentado anteriormente en formato de tabla, para poder apreciar sus cualidades de una manera mas visual e intuitiva.

6.3 Grandes presas españolas — Tipología

6.3.1 Presa de Almendra

Situada en el río Tormes, en la provincia de Salamanca. Construida entre 1963 y 1970. Es la presa más alta de España, con 202 metros de altura; longitud de coronacion de 567 metros; y tiene una gran capacidad. Forma parte del sistema hidroeléctrico del Duero Y ES una de las mayores instalaciones hidroeléctricas de España.

Se trata de una presa de arco, debido a su entorno donde tiene lugar un valle estrecho y roca granítica adecuada, lo que permite transmitir el empuje del agua a los estribos. Esta construida en hormigón y se caracterizada por su gran esbeltez y curvatura aguas abajo.

6.3.2 Presa de Ricobayo

Situada en el río Esla, en la provincia de Zamora. Construida en la década de 1930. Tiene unos 94 metros de altura y un embalse amplio que desempeña un papel fundamental en la regulación del río Esla y en la producción hidroeléctrica de la zona.

Se trata de una presa de gravedad, ya que el entorno no presenta un valle suficientemente estrecho para una presa de arco. Su estabilidad depende principalmente de su gran peso, para lo cual se utilizó una estructura maciza de hormigón asentada sobre un cimiento rocoso adecuado.

6.3.3 Presa de Aldeadávila

Situada en el río Duero, en la frontera entre España y Portugal. Construida entre 1958 y 1964. Tiene aproximadamente 139 metros de altura y alimenta una de las centrales hidroeléctricas más potentes de España, con instalaciones subterráneas de gran importancia estratégica.

Se trata de una presa de arco-gravedad, elegida por la presencia de un profundo cañón rocoso que permite combinar la curvatura para transmitir esfuerzos a los estribos con el peso propio para garantizar estabilidad. Está construida en hormigón y aprovecha al máximo la topografía para obtener un gran salto hidráulico.

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