La presa de El Atazar (Grupo 44)

Trabajo realizado por estudiantes
Título	La presa de El Atazar (Grupo 44)
Asignatura	Teoría de Campos
Curso	2025-26
Autores	Marina Lacho Mora José Luis Leines Almeida Rocío Martín Renzini Jaime García Suárez Pablo Fernández Arce
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

La presa de El Atazar, situada en el valle del río Lozoya al norte de la Comunidad de Madrid, es la infraestructura hidráulica más relevante del sistema de abastecimiento regional. Construida entre 1968 y 1972, es una presa de arco-gravedad que combina la acción del peso propio con el trabajo de la doble curvatura, aprovechando al máximo la geometría del valle. Con sus 134 m de altura y aproximadamente 484 m de coronación, genera un embalse de más de 425 hm³, fundamental para garantizar el suministro de agua a Madrid. Además de su capacidad, destaca por su diseño estructural y por el papel clave que desempeña dentro del conjunto de presas gestionadas en la cuenca del Lozoya.

El objetivo de este trabajo es estudiar la geometría y el comportamiento hidrostático de la presa mediante herramientas de modelado y visualización en MATLAB. Para ello se representará la superficie del paramento de aguas arriba, se analizará el campo de presiones y se compararán dos configuraciones geométricas en términos de la fuerza que soportan. Además, se realizará una breve revisión de los principales tipos de presas y de algunos ejemplos destacados en España.

1 Introducción

1.1 Modelo geométrico de la presa

Se analiza la superficie de la presa en su paramento de aguas arriba, es decir, la zona en contacto directo con el embalse. Dado que la presa presenta doble curvatura, su trazado adopta una forma de arco circular en planta y un arco de tipo parabólico en la sección vertical. Utilizando un sistema de coordenadas cilíndricas (𝜌,𝜃,𝑧), cuyo origen se sitúa en el centro del valle y con el eje 𝑧 dirigido hacia arriba, la geometría se describe mediante la siguiente superficie:

[math]\rho = \rho_{0} + b (1 - \frac{z^2}{h^2})[/math],[math]\qquad θ ∈ [\frac{2π}{3}, \frac{4π}{3}][/math], [math]\qquad Z ∈ [0,H][/math]

Donde:

• 𝐻 = 134 m: altura de la presa.
• [math]\rho_{0}[/math]= 150 m: radio en la coronación (altura máxima).
• 𝑏 = 40 m: parámetro de curvatura parabólica.

1.2 Campo de presión hidrostática

Cuando el embalse está lleno hasta la altura [math]H_{agua}[/math] medida desde la base de la presa, la presión que ejerce el agua sobre el paramento de aguas arriba viene dada por:

[math]P(z)=P_{0}+g\varrho_{agua}(H_{agua}-z)[/math] [math]\qquad z\in \left[0,H_{agua}\right][/math]

Donde P₀ es la presión representa la presión atmosférica, [math]\varrho_{agua}[/math] es la densidad del agua y g la aceleración de la gravedad. En este estudio se toma [math]H_{agua}=125m[/math].

El campo de fuerzas asociado a dicha presión sobre la superficie puede expresarse como:

[math]\overrightarrow{F} = −P(z) \overrightarrow{n}[/math]

Siendo [math]\overrightarrow{n}[/math] el vector normal a la superficie, orientado hacia el interior del embalse.

2 Superficie parametrizada

A la derecha nos encontramos con la superficie parametrizada de la presa en su cara aguas arriba.
La presa se modela mediante una parametrización en coordenadas cilíndricas (𝜌,𝜃,𝑧), donde:
1. Curvatura Horizontal (Arco Circular): El ángulo 𝜃 define la forma de arco circular del plano horizontal, que permite que la presa resista la presión del agua por compresión
2. Curvatura Vertical (Arco Parabólico): El radio 𝜌 se hace dependiente de la altura 𝑧 mediante una función parabólica

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

Esto provoca que el radio sea máximo en la base (𝑧=0) y mínimo en la coronación (𝑧=H).Esta variación optimiza la distribución de tensiones, ya que las presas son más gruesas y necesitan un radio mayor en la base donde la presión del agua es máxima.

Superficie parametrizada de la presa realizada con el programa MATLAB

% Parámetros
H = 134;           % Altura de la presa
rho0 = 150;        % Radio en la coronación
b = 40;            % Parámetro de curvatura parabólica

% Rango de parámetros
theta = linspace(2*pi/3, 4*pi/3, 200);    % Ángulo (cara aguas arriba)
z = linspace(0, H, 200);                  % Altura

% Crear malla
[Theta, Z] = meshgrid(theta, z);

% Expresión del radio en función de z
Rho = rho0 + b*(1 - (Z.^2)/(H^2));

% Convertir a coordenadas cartesianas
X = Rho .* cos(Theta);
Y = Rho .* sin(Theta);

% Gráfica de la superficie
figure
surf(X, Y, Z)
shading interp
colormap('turbo')
xlabel('x')
ylabel('y')
zlabel('z')
title('Cara aguas arriba de la presa (parámetro cilíndrico)')
axis equal
grid on

6 Sedimentación en el embalse

Los sedimentos transportados por el agua del río se depositan en el fondo del embalse, reduciendo el volumen, y por lo tanto, reduciendo su capacidad. Modelamos la concentración de sedimentos depositados en el fondo del embalse (en kg/m²) como:

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

Donde [math]S_{0}=50kg/m^2[/math] es la sedimentación base, [math]\alpha = 3[/math] modela la mayor acumulación cerca de la entrada del río, y [math]L=500m[/math] es una escala característica. Para simplificar, considera el fondo del embalse como 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.

6.1 Representación

Sedimentación en el fondo hecho con el programa MATLAB

clear, clc, close all
%Parámetros
S0 = 50;
alpha = 3;
L = 500;

%Dominio del fondo del embalse
x = linspace(-700,700,300);
y = linspace(-700,700,300);
[X,Y] = meshgrid(x,y);

%Campo escalar S(x,y)
R = X.^2 + Y.^2;
S = S0*(1+ alpha*exp(-R/L^2));

%Figura: mapa de colores
figure('Units', 'normalized')
p = pcolor(X, Y, S);
set (p, 'EdgeColor', 'none');
colorbar;
axis equal;
title('Campo escalar S(x,y)');

6.2 Gradiente de sedimentación

Gracias al gradiente, sabremos hacia qué dirección aumenta más rápido la sedimentación y qué tan rápido cambia.

Dado:

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

La derivada en función de x: [math] \frac{\partial S}{\partial x}=S_{0}\alpha e^{(-\frac{x^2+y^2}{L^2})}\left ( -\frac{2x}{L^2} \right )[/math]

La derivada en función de y: [math] \frac{\partial S}{\partial y}=S_{0}\alpha e^{(-\frac{x^2+y^2}{L^2})}\left ( -\frac{2y}{L^2} \right )[/math]

Por lo tanto, el resultado de mi gradiente es:

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

Veamos de una forma sencilla ahora hacia qué lado apunta. Para ello, simplificaremos el gradiente de la siguiente manera:

[math]∇S(x,y)=\left ( -\frac{2x}{L^2}(algo positivo),-\frac{2y}{L^2}(algo positivo) \right )[/math]

Vamos a centrarnos en lo que no es positivo.

Supongamos que [math]x\gt0[/math], entonces [math]-\frac{2x}{L^2}\lt0[/math]. Por lo tanto, la flecha apunta hacia la izquierda, porque los valores negativos de x están a la izquierda.

Ahora supongamos [math]x\lt0[/math], entonces [math]-\frac{2x}{L^2}\gt0[/math]. Por lo tanto, la flecha apunta hacia la derecha, porque los valores positivos de x están a la derecha.

Como vemos, sin importar que x tomemos, apuntará al origen.

Ocurre lo mismo si suponemos [math]y\gt0, y\lt0[/math].

Por lo tanto, el gradiente apunta al centro de la presa, [math](0,0)[/math]

¿Cómo interpretamos estos resultados? En las zonas cercanas al centro del embalse, el gradiente es grande. Esto refleja que la sedimentación crece de forma pronunciada hacia el punto por donde entra el río. En las regiones más alejadas, el gradiente disminuye, lo que implica que la concentración cambia más lentamente y la acumulación es menos notable.

La magnitud del gradiente puede visualizarse mediante un mapa de colores: los tonos más fuertes corresponden a variaciones intensas de la concentración. Por su parte, la dirección del gradiente se representa con flechas negras, que muestran hacia dónde se desplaza el flujo de sedimentos hacia áreas de mayor concentración.

Gradiente de sedimentación hecho con el programa MATLAB

clc, clear, close all
% Parámetros
S0 = 50;
alpha = 3;
L = 500;

%Dominio
x = linspace(-400,400,200);
y = linspace(-400,400,200);
[X,Y] = meshgrid(x,y);

% Radio y campo S
R = X.^2 + Y.^2;
S = S0*(1 + alpha*exp(-R / L^2));

% Gradiente
coef = -2 * S0 * alpha / (L^2);
dx = coef * X .* exp(-R / L^2);
dy = coef * Y .* exp(-R / L^2);

% Magnitud del gradiente
mag = sqrt(dx.^2 + dy.^2);

% Figura
figure('Units','normalized','Position',[0.1 0.1 0.6 0.6])
p = pcolor(X, Y, mag);
set(p,'EdgeColor','none')
colorbar

hold on
axis equal
xlabel('x (m)');
ylabel('y (m)');
title('Representación gradiente de sedimentación')
step = 6;
quiver(X(1:step:end,1:step:end), Y(1:step:end,1:step:end), dx(1:step:end,1:step:end), dy(1:step:end,1:step:end), 'k')
plot(0,0,'y.','MarkerSize',10)
legend({'Magnitud del gradiente','Campo vectorial','Origen'})     
hold off

6.3 Masa sedimentos

Para calcular la masa total de sedimentos depositados utilizamos la fórmula:

6.4 Curvas de isoconcentración

Una curva de isoconcentración es una línea que une todos los puntos donde la concentración es la misma. [math](S(x,y)=cte)[/math]

En mi función, [math]S_{0}[/math], [math]\alpha[/math] y [math]1[/math] son constantes, por lo que el único término que varía es [math]e^{-\frac{x^2+y^2}{L^2}}[/math]. Como [math]L^2[/math] también es constante, mi función depende solo de [math]x^2+y^2=r^2[/math]. Eso significa que S depende solo de la distancia al origen, así que mis curvas de nivel serán círculos.

Gradiente de sedimentación hecho con el programa MATLAB

clc, clear, close all

% Parámetros
S0 = 50;
alpha = 3;
L = 500;

% Dominio 
x = linspace(-800,800,300);
y = linspace(-800,800,300);
[X,Y] = meshgrid(x,y);

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

% Representación de curvas de nivel
figure;
contour(X, Y, S, 10, 'LineWidth', 1.5);
axis equal;
xlabel('x (m)');
ylabel('y (m)');
title('Curvas de isoconcentración de sedimentación S(x,y)');
grid on;

7 Flujo de entrada del río

El río Lozoya entra al embalse por el extremo norte. Los modelos a continuación describen cómo se mueve el agua al entrar en el embalse, mostrando que el flujo es más intenso cerca de la entrada, y que disminuye rápidamente con la distancia. Además, incluye una pequeña componente vertical que representa la mezcla interna del agua.

Todo ello permite analizar el transporte de sedimentos, la mezcla de temperaturas y la dinámica interna del embalse.

Modelado en 2D:

[math]\vec{v}(x,y)=v_{0}e^{-\frac{x^{2+y^2}}{R^2}}(cos\varphi \vec{i}+sen\varphi \vec{j})[/math]

Donde \(v_{0} = 0.5\,\text{m/s}\) es la velocidad máxima en el centro de entrada, \(R = 30\,\text{m}\) es el radio característico, y \(\varphi = \pi/6\) es el ángulo de entrada. Se considera el dominio \((x,y) \in [-100\,\text{m},\, 100\,\text{m}]^{2}\), \(\sqrt{x^{2}+y^{2}}\leq \rho _{a}\), donde \(\rho_{a}\) es el radio de la presa a la altura \(z = H_{\text{agua}}\). El factor exponencial indica que el agua entra concentrada y luego se abre perdiendo fuerza.

Modelado en 3D:

[math]\vec{v}(x,y,z)= v_{0}\, e^{-\frac{x^{2}+y^{2}}{R^{2}}}\left( \cos\varphi\, \hat{\imath} + \sin\varphi\, \hat{\jmath} \right)\;+\;v_{1}\,\sin\left( \frac{\pi z}{H_{\text{agua}}} \right)\, e^{-\frac{x^{2}+y^{2}}{R^{2}}} \hat{k}.[/math]

Donde \(v_{1} = 0.1\,\text{m/s}\) es la velocidad máxima de la componente vertical, y \(z \in [0, H_{\text{agua}}]\).

7.1 Campo vectorial de la velocidad

Plano horizontal para Z = 125 hecho con el programa MATLAB

Plano vertical para y = 0 hecho con el programa MATLAB

clear, clc, close all
% Parámetros
v0 = 0.5;        % Velocidad horizontal máxima (m/s)
v1 = 0.1;        % Velocidad vertical máxima (m/s)
R = 30;          % Radio característico (m)
phi = pi/6;      % Ángulo de entrada (rad)
H = 125;         % Altura del agua (m)

% Dominio horizontal (plano superficial z = H)
xvec = -100:4:100;
yvec = -100:4:100;
[x, y] = meshgrid(xvec, yvec);
z_surface = H;   % Plano superficial

% Factor exponencial
E = exp(-(x.^2 + y.^2)/R^2);

% Componentes del campo en superficie
vx = v0 * cos(phi) .* E;
vy = v0 * sin(phi) .* E;
vz_surface = v1 * sin(pi*z_surface/H) .* E;

% Representación plano horizontal (z = H)
figure('Name', 'Campo en superficie (z = H)');
quiver(x, y, vx, vy,'k');
axis equal;
xlim([min(xvec) max(xvec)]);
ylim([min(yvec) max(yvec)]);
xlabel('x (m)');
ylabel('y (m)');
title(sprintf('Campo de velocidad en superficie (z = %d m)', H));
grid on;

% Representación en plano vertical (y = 0)
xz = linspace(-100, 100, 61);
zz = linspace(0, H, 41);
[XZ, ZZ] = meshgrid(xz, zz);

% Factor exponencial solo depende de X en este plano
E_vert = exp(-(XZ.^2)/R^2);

% Componentes del campo en el plano vertical
VX_vert = v0 * cos(phi) .* E_vert;
VZ_vert = v1 * sin(pi*ZZ/H) .* E_vert;

% Figura vertical
figure('Name', 'Campo en plano vertical y=0', 'NumberTitle', 'off');
quiver(XZ, ZZ, VX_vert, VZ_vert,'k');
xlabel('x (m)');
ylabel('z (m)');
title('Plano vertical (y = 0)');
axis tight;
grid on;

7.2 Divergencia del campo

La divergencia de un campo vectorial [math]\vec{v}=(v_{x},v_{y},v_{z})[/math] se define: [math] ∇\cdot \vec{v}=\frac{\partial v_{x}}{\partial x}+\frac{\partial v_{y}}{\partial y}+\frac{\partial v_{z}}{\partial z}[/math]

La divergencia mide si en un punto el flujo se está expandido o concentrando:

• [math] ∇\cdot \vec{v}\gt0[/math] indica que el flujo se abre o se dispersa.
• [math] ∇\cdot \vec{v}\lt0[/math] indica que el flujo se acumula.

En el caso del embalse, calcular la divergencia del campo de velocidad del agua que entra por el río nos permite:

• Identificar zonas de acumulación de agua y sedimentos, donde el flujo converge ([math] ∇\cdot \vec{v}\lt0[/math])
• Detectar regiones donde el flujo se dispersa, con menor sedimentación ([math] ∇\cdot \vec{v}\gt0[/math])
• Comprender la dinámica interna del embalse, cómo se mezcla el agua y cómo se distribuyen los sedimentos.

Mi campo es:

[math]\vec{v}(x,y,z)= v_{0}\, e^{-\frac{x^{2}+y^{2}}{R^{2}}}\left( \cos\varphi\, \hat{\imath} + \sin\varphi\, \hat{\jmath} \right)\;+\;v_{1}\,\sin\left( \frac{\pi z}{H_{\text{agua}}} \right)\, e^{-\frac{x^{2}+y^{2}}{R^{2}}} \hat{k}.[/math]

Donde [math]v_{x}=v_{0}cos(\varphi)e^{-\frac{x^{2}+y^{2}}{R^2}},\qquad v_{y}=v_{0}sin(\varphi)e^{-\frac{x^{2}+y^{2}}{R^2}},\qquad v_{z}=v_{1}sin\left ( \frac{\pi z}{H_{agua}} \right )e^{-\frac{x^2+y^2}{R^2}}[/math]

Empecemos entonces haciendo las derivadas parciales:

• Para [math]v_{x}[/math]: [math]\frac{\partial v_{x}}{\partial x}=v_{0}cos(\varphi)\cdot \frac{\partial}{\partial x}e^{-\frac{x^{2}+y^{2}}{R^2}}=v_{0}cos(\varphi)e^{-\frac{x^{2}+y^{2}}{R^2}}\cdot \left ( -\frac{2x}{R^2} \right )= -\frac{2xv_{0}cos\varphi }{R^{2}}e^{-\frac{x^{2}+y^{2}}{R^2}}[/math]

• Para [math]v_{y}[/math]: [math]\frac{\partial v_{y}}{\partial y}=v_{0}sin(\varphi)\cdot \frac{\partial}{\partial y}e^{-\frac{x^{2}+y^{2}}{R^2}}=v_{0}sin(\varphi)e^{-\frac{x^{2}+y^{2}}{R^2}}\cdot \left ( -\frac{2y}{R^2} \right )= -\frac{2yv_{0}sin\varphi }{R^{2}}e^{-\frac{x^{2}+y^{2}}{R^2}}[/math]

• Para [math]v_{z}[/math]: [math]\frac{\partial v_{z}}{\partial z}=v_{1}e^{-\frac{x^2+y^2}{R^2}}\frac{\partial }{\partial z}\left [ sin\left ( \frac{\pi z}{H_{agua}} \right )\right ]=v_{1}e^{-\frac{x^2+y^2}{R^2}}\cdot \frac{\pi }{H_{agua}}cos\left ( \frac{\pi z}{H_{agua}} \right )[/math]

Y ahora sumamos:

[math]∇\cdot \vec{v}=\frac{\partial v_{x}}{\partial x}+\frac{\partial v_{y}}{\partial y}+\frac{\partial v_{z}}{\partial z}=e^{-\frac{x^2+y^2}{R^2}}\left [-\frac{2xv_{0}cos\varphi }{R^{2}}-\frac{2yv_{0}sin\varphi }{R^{2}}+\frac{\pi v_{1}}{H_{agua}}cos(\frac{\pi z}{H_{agua}}) \right ] [/math]

Por lo tanto, la divergencia máxima se da en el centro de la entrada del río y en la superficie [math](0, 0, 0)[/math]

¿Dónde es máximo? Los términos [math]-\frac{2xv_{0}cos\varphi }{R^{2}}[/math] y [math]-\frac{2yv_{0}sin\varphi }{R^{2}}[/math] son [math]0[/math] cuando [math]x=0, y=0[/math]. El término vertical [math]\frac{\pi v_{1}}{H_{agua}}cos\left (\frac{\pi z}{H_{agua}}\right )[/math] es máximo cuando [math]cos\left (\frac{\pi z}{H_{agua}}\right )=1\Rightarrow z=0[/math]. El término exponencial [math]e^{-\frac{x^2+y^2}{R^2}}[/math] es máximo en [math]x=0, y=0[/math]

Por lo tanto, la divergencia máxima se da en el centro de la entrada del río y en la superficie [math](0, 0, 0)[/math]

Divergencia en 2D hecha en MATLAB

Divergencia en 3D hecha en MATLAB

clear, clc, close all
% Parámetros
v0 = 0.5;      % m/s, velocidad horizontal máxima
v1 = 0.1;      % m/s, velocidad vertical máxima
R  = 30;       % m, radio característico
phi = pi/6;    % ángulo de entrada
H = 125;       % m, altura del agua

% Dominio horizontal (plano z = H)
xv = linspace(-100,100,100);
yv = linspace(-100,100,100);
[x, y] = meshgrid(xv, yv);

z_surface = H;  % Plano superficial (z = H)

% Factor exponencial
E = exp(-(x.^2 + y.^2)/R^2);

% Coeficientes para la divergencia
coefx = -2 * v0 * cos(phi) / R^2;
coefy = -2 * v0 * sin(phi) / R^2;
coefz = v1 * (pi / H);

% Divergencia
diver = E .* (coefx .* x + coefy .* y + coefz .* cos(pi * z_surface / H));

% Representación 3D 
figure('Name', 'Divergencia en superficie');
surf(x, y, diver);
shading interp
hold on;
contour3(x, y, diver, 20, 'b'); 
hold off;
title('Divergencia (plano superficial z = H)');
colorbar;
view(30,25);
axis tight;

% Representación en mapa de colores 2D 
figure('Name', 'Mapa de colores divergencia (superficie)');
imagesc(xv, yv, diver);
axis xy; axis equal;
xlabel('x (m)'); ylabel('y (m)');
title('Mapa de colores: divergencia en superficie (z = H)');
colorbar;

7.3 Rotacional del campo

En este estudio, el rotacional del campo de velocidad permite identificar si el flujo del río Lozoya entra al embalse con vorticidad, es decir, con tendencia a girar o formar remolinos. Calcular el rotacional nos indica dónde el agua presenta zonas de giro tanto en el plano horizontal como en el vertical. Esto es importante porque la vorticidad influye en la estabilidad del flujo, en la mezcla del agua, en el transporte de sedimentos y en el modo en que el chorro de entrada interactúa con el embalse.

En este caso, como demostraremos a continuación, el rotacional no es cero, por lo que el flujo presenta remolinos y rotación, especialmente cerca de los bordes del chorro de entrada, lo que confirma que la entrada de agua no es uniforme ni completamente laminar.

Definiremos: [math]f(x,y)=e^{-\frac{x^2+y^2}{R^2}}, \qquad g(z)=sin\left ( \frac{\pi z}{H_{agua}} \right )[/math]

De esta forma: [math]v(x,y,z)=v_{0}f(x,y)(cos(\varphi) \vec{i}+sin(\varphi )\vec{j})+v_{1}f(x,y)g(z)\vec{k}[/math]

En coordenadas cartesianas:

[math]∇\times v=\begin{vmatrix} \vec{i}& \vec{j} & \vec{k} \\ \frac{\partial }{\partial x}& \frac{\partial }{\partial y} & \frac{\partial }{\partial z}\\ v_{x}& v_{y} & v_{z} \\ \end{vmatrix}=\left ( \frac{\partial v_{z}}{\partial y}-\frac{\partial v_{y}}{\partial z} \right )\vec{i}+\left ( \frac{\partial v_{x}}{\partial z}-\frac{\partial v_{z}}{\partial x} \right )\vec{j}+\left ( \frac{\partial v_{y}}{\partial x}-\frac{\partial v_{x}}{\partial y} \right )\vec{k}[/math]

Siendo: [math]v_{x}=v_{0}f(x,y)cos(\varphi) \qquad v_{y}=v_{0}f(x,y)sin(\varphi) v_{z}=v_{1}f(x,y)g(z) [/math]

Vemos que:

• [math]v_{x}, v_{y}[/math] dependen solo de [math]x, y[/math].
• [math]v_{z}[/math] depende de [math]x, y, z[/math].

Calculamos:

[math]\left ( -\frac{2y}{R^2}v_{1}fg(z)-0\right )\vec{i}+\left ( 0--\frac{2x}{R^2}v_{1}fg(z)\right )\vec{j}+\left ( \frac{\partial v_{y}}{\partial x}-\frac{\partial v_{x}}{\partial y} \right )\vec{k}[/math]

[math][/math]