Revisión del 20:04 6 dic 2025

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.

Figura 1: La presa de El Atazar vista desde arriba

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:

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:

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:

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

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.

Figura 2: 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

3 Presiones sobre la presa

La presa responde a un campo de presiones modelizado por la función:

donde:

• P(z) = presión a la altura z
• P₀ = presión atmosférica
• Q_w = densidad del agua, 1000kg/m³
• g = gravedad, 9,81g/m²
• H_w = altura a la que se encuentra llena la presa, 125m
• z = altura

Figura 3: Representación de presiones de la presa realizada con el programa MATLAB

% Parámetros geométricos 
H = 134;      % Altura de la presa
Hw = 125;     % Altura agua
rho0 = 150;   % Radio en la coronacion
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);

% Datos densidad y gravedad
Po = 1000;
g = 9.81;

% Relación presión-altura
P = Po * g * (Hw - Z);

% Gráfica de la superficie
figure;
surf(X, Y, Z, P, "EdgeColor", "none");
title('Campo escalar de presión en la presa');
xlabel('x');
ylabel('y');
zlabel('z');
axis equal;
grid on;
colormap("jet");
colorbar;

Observando la Figura 3, vemos que la presión es lineal, aumentando conforme vamos perdiendo altura. Esto sigue la lógica de la fórmula fundamental de la estática de fluidos.

4 Tensiones

En este apartado se determinará la fuerza total debida a la presión hidrostática que actúa sobre la cara aguas arriba de la presa, así como la presión media (o tensión superficial equivalente) para los dos modelos geométricos analizados:

• Presa doble curvatura: [math]b=40m[/math]
• Presa de curvatura simple (cilíndrica): [math]b=0m[/math]

El objetivo es comparar ambos modelos y analizar cuál presenta una mejor resistencia frente a la presión ejercida por el agua.

Para determinar cómo actúa la presión hidrostática sobre la presa, se parte de su distribución en función de la profundidad. El agua transmite presión en todas direcciones y ésta depende únicamente de la profundidad. Así, en cualquier punto de la superficie de la presa, la presión viene dada por:

Quitamos [math]p_{0}[/math] ya que, para no ir arrastrando la presión atmosférica, tomamos manométricas.

Donde [math]\rho[/math] es la densidad del agua, [math]g[/math] es la aceleración de la gravedad, [math]H[/math] es la altura máxima de la lámina de agua y [math]z[/math] es la distancia medida desde la base de la presa.

Para calcular la fuerza, utilizaremos la siguiente fórmula: [math]F=\iint_{}^{}P(z)dA[/math]

Como mencionamos en la introducción, la forma del paramento aguas arriba se define como:

Donde:

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

Y además sabemos que nuestro vector de posición es: [math]\vec{r}(\theta ,z)=(\rho cos\theta, \rho sin\theta, z)[/math].

El área:

Ya con esto, podemos resolver:

Cada elemento diferencial de fuerza actúa perpendicularmente a la superficie:

y la fuerza total se obtiene integrando sobre toda la superficie de la presa: [math] \quad F = \int_{\theta_1}^{\theta_2} \int_0^H p(z) \, \rho(z) \sqrt{1 + (\rho'(z))^2} \, dz \, d\theta [/math]

Evaluando numéricamente con los parámetros del enunciado:

• Para [math]b = 40~\text{m}[/math] (doble curvatura): [math]F = 3,027 \times 10^{10}~\text{N}[/math]
• Para [math]b = 0~\text{m}[/math] (curvatura simple): [math]F = 2,407 \times 10^{10}~\text{N}[/math]

A continuación, para conocer la tensión promedio que soporta la presa, calculamos: [math] \sigma = \frac{F}{A} [/math]
Donde [math]A[/math] es el área de la superficie: [math] A = \int_{\theta_1}^{\theta_2} \int_0^H \rho(z) \sqrt{1 + (\rho'(z))^2} \, dz \, d\theta [/math]
Resultados numéricos son:

• [math]b = 40~\text{m}[/math]: [math]A = 4,89\times10^4~\text{m²}\;\Longrightarrow\;\sigma = 0,613~\text{MPa}[/math]
• [math]b = 0~\text{m}[/math]: [math]A = 3,92\times10^4~\text{m²}\;\Longrightarrow\; \sigma = 0,619~\text{MPa}[/math]

En resumen, aunque la presa con doble curvatura recibe una fuerza global más alta, la tensión promedio que soporta es menor. Esto significa que los esfuerzos se reparten de manera más equilibrada sobre toda la superficie, haciendo que la estructura responda mejor frente al empuje del agua. Por ello, la configuración de doble curvatura es la que ofrece un comportamiento más seguro, ya que evita acumulaciones de tensión y disminuye la probabilidad de fallos en el material.

5 Tipos de presas

Figura 4: Tipos de presas

Las presas pueden clasificarse de diversas maneras según la forma en que resisten el empuje del agua y los materiales empleados en su construcción. Cada tipología presenta ventajas e inconvenientes que condicionan su idoneidad para un emplazamiento concreto. En líneas generales, podemos clasificarlas en dos grupos: presas de fábrica (de hormigón) y presas de materiales sueltos, siendo estas últimas las más comunes por su versatilidad.

5.1 Presas de fábrica

Se tratan de presas relativamente esbeltas, construidas con hormigón. Este tipo de presas podemos subclasificarlas a su vez en tres grandes grupos:

5.1.1 Presas de gravedad

Figura 5: Presas de gravedad

Las presas de gravedad se caracterizan por tener un cuerpo macizo —generalmente de hormigón— cuyo propio peso es el que resiste la presión del agua, evitando el deslizamiento. Su estabilidad depende en gran medida de la solidez y homogeneidad del terreno de cimentación, por lo que suelen construirse sobre macizos rocosos firmes. Para impedir el vuelco, la resultante de los empujes (peso propio + presión del agua) debe quedar siempre dentro de la base de la presa. Estas estructuras se ejecutan casi por completo con hormigón en masa, utilizando armaduras solo en zonas puntuales que puedan estar sometidas a tracción, como las galerías interiores. Como trabajan principalmente a compresión, es esencial controlar las tensiones de tracción, especialmente en el pie de aguas arriba, uno de los puntos más críticos desde el punto de vista estructural. Su principal ventaja es la robustez y sencillez estructural, ya que no dependen de la forma del valle ni de transmitir cargas a los estribos. Sin embargo, necesitan un gran volumen de materiales, lo que eleva los costes y las hace menos eficientes en términos económicos y ambientales. Son adecuadas cuando el valle es ancho y existe buena calidad del terreno en el fondo.

5.1.2 Presas de arco o bóveda

Figura 6: Presas de arco o bóveda

Por otro lado, las presas de arco o bóveda aprovechan la curvatura de su estructura para transmitir la mayor parte del empuje del agua hacia los estribos del valle. Gracias a este comportamiento, es posible reducir significativamente la cantidad de hormigón, lo que las convierte en soluciones especialmente eficientes. El arco actúa como antifunicular de la carga radial uniformemente repartida; es decir, trabaja esencialmente a esfuerzos axiles de compresión, sin tracciones relevantes. Estos esfuerzos se transmiten directamente a los estribos, que deben ser rocas muy resistentes capaces de soportar las elevadas cargas concentradas. Para que este tipo de presa sea viable, es fundamental que el valle sea estrecho o en forma de garganta, y que los laterales presenten una alta calidad geológica. Aunque su diseño estructural es más complejo y exige un análisis geotécnico más riguroso, las presas de arco permiten alcanzar grandes alturas con secciones delgadas y un uso mucho más eficiente de materiales. La sección de este tipo de presas puede ser:

• Presa arco : sección trapezoidal de planta curva.
• De doble curvatura : sección curva de planta curva. También se conocen como presas en bóveda.
• Arco-gravedad : combinan la curvatura de una presa de arco con parte del peso propio de una presa de gravedad. Esto les permite adaptarse a valles que no son tan estrechos como los necesarios para una presa de arco pura, pero que tampoco justifican un muro macizo típico de gravedad. Su diseño permite ahorrar materiales respecto a una presa de gravedad, aunque exige buen apoyo lateral y estudios geológicos rigurosos. Son un compromiso que aprovecha parcialmente las ventajas de ambos sistemas.

5.1.3 Presas de contrafuertes

Figura 7: Presas de contrafuertes

También existen las presas de contrafuertes. Este tipo de presa presenta un mecanismo resistente similar al de las presas de gravedad, aunque con un diseño más ligero. La parte principal o cabeza de la presa se encuentra inclinada y apoyada sobre una serie de contrafuertes —elementos longitudinales o “colas”— situados aguas abajo, que transmiten los esfuerzos al terreno a lo largo de toda la estructura.

A diferencia de las presas de gravedad, estas estructuras prácticamente no requieren supresión, y junto con la contribución del peso de la cuña de agua sobre la cabeza, permiten reducir notablemente la masa de hormigón empleada. En cuanto a proporciones, el ancho de la base suele ser del orden de la altura de la presa, lo que las hace esbeltas en comparación con las de gravedad maciza.

Sin embargo, su forma más compleja exige mayor cantidad de mano de obra y precisión constructiva, lo que incrementa los costes pese al menor volumen de material. Por este motivo, y dado que existen alternativas más eficientes, hoy en día apenas se construyen presas de contrafuertes..

5.2 Presas de materiales sueltos

Figura 8: Presas de materiales sueltos

Las presas de materiales sueltos son estructuras muy versátiles que pueden construirse prácticamente con cualquier tipo de material disponible en la zona —gravas, arenas, arcillas o materiales seleccionados—, lo que explica que sean las presas más abundantes del mundo. Presentan una sección trapezoidal y son mucho menos esbeltas que las presas de fábrica (hormigón), debido a que necesitan taludes amplios para garantizar su estabilidad.Su característica más distintiva es la zonificación, es decir, la colocación deliberada de cada tipo de material en la zona donde mejor cumple su función: un núcleo impermeable que minimiza filtraciones, espaldones de materiales más permeables que aportan estabilidad y filtros o drenes que controlan el flujo de agua. Esta organización interna permite optimizar recursos locales y adaptar la presa a condiciones geotécnicas muy diversas..

5.3 Ejemplos en España

Figura 9: Presa de Almendra

En España existen varios ejemplos destacados que ilustran la aplicación de estas tipologías:

La Presa de Almendra, situada sobre el río Tormes, es una de las más altas de España y se clasifica como presa de bóveda o doble curvatura. Su ubicación, en un valle profundo con paredes rocosas muy estables, es ideal para este tipo de presa, ya que permite transmitir el empuje hacia los estribos y reducir considerablemente el volumen de hormigón. Esta tipología resulta especialmente adecuada para grandes alturas y para ubicar centrales hidroeléctricas de alta potencia aprovechando saltos muy profundos.

Figura 10: Presa de Ricobayo

La Presa de Ricobayo, ubicada en el río Esla, es un ejemplo clásico de presa de gravedad. En este tramo el valle es más ancho y cuenta con cimientos graníticos muy estables, adecuados para soportar el peso de un muro masivo. Su diseño resulta más sencillo y robusto que el de una presa de arco, aunque su construcción implica emplear un volumen de hormigón considerablemente mayor. Este tipo de solución es óptima cuando las condiciones geológicas del fondo del valle son favorables, pero los estribos laterales no presentan la resistencia necesaria como para transmitir grandes cargas, lo que descarta la tipología de arco.

Figura 11: Presa de Aldeadávila

La Presa de Aldeadávila, situada en el río Duero, es una de las obras hidráulicas más relevantes de España y se clasifica como presa de arco-gravedad. Su emplazamiento en los Arribes del Duero, un desfiladero estrecho y profundamente encajado en roca, hace que la tipología de arco resulte adecuada; sin embargo, la geometría del valle y las altas exigencias de estabilidad llevaron a optar por una solución mixta. En este diseño, una parte del empuje del agua se transmite por curvatura a los estribos, mientras que otra parte se resiste mediante el peso propio del muro, combinando así la eficiencia estructural del arco con la robustez de una presa de gravedad. Gracias a ello, Aldeadávila logra un excelente comportamiento frente a solicitaciones muy elevadas y permite alojar uno de los complejos hidroeléctricos más potentes del país.

Figura 12: Presa de El Atazar

Otro ejemplo es mismamente la presa sobre la que va este trabajo, la Presa de El Atazar. Es una de las presas más singulares de España y también corresponde a una tipología de arco-gravedad. Su diseño curvado y su considerable altura aprovechan la estrechez del valle y la calidad de la roca de los estribos, al mismo tiempo que la masa central de la presa contribuye a resistir el empuje del agua. Esta combinación permitió optimizar el gasto de material y obtener una estructura muy estable. Su forma y funcionamiento demuestran cómo en algunos valles de anchura intermedia resulta más adecuado emplear un diseño mixto que uno puramente de arco o puramente de gravedad.

En conjunto, la elección del tipo de presa depende principalmente de la forma del valle, la resistencia de los materiales geológicos, la altura deseada y los recursos económicos disponibles.

En resumen

• Las presas de gravedad aportan gran estabilidad pero requieren mucho material.
• Las de arco y bóveda son más eficientes, pero dependen de un terreno lateral muy resistente.
• Las de arco-gravedad son una solución intermedia que proporciona flexibilidad.
• Las de contrafuertes permiten ahorrar material a costa de mayor complejidad constructiva.
• Las de materiales sueltos son muy versátiles, depende del material que utilicen.
  

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:

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

Figura 13: 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:

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:

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

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.

Figura 14: 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 aplicamos la fórmula: [math]\quad M = \iint_{A} S(x,y)\, dA [/math].

Ahora para facilitar los cálculos podemos usar coordenadas polares debido a que el dominio de integración es un sector circular:
Polares: \begin{cases} x = r\cos\theta & ,\text{donde } r^2 = x^2 + y^2 \\ y = r\sin\theta \\ \end{cases}
Resolvemos la integral:

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

Separando términos:

[math] M = S_0 \theta \int_0^R r\,dr + S_0 \alpha \theta \int_0^R r e^{-r^2/L^2}\,dr \qquad \text{donde } \theta =\theta _{1}-\theta _{0} [/math]

Las integrales son

[math] \int_0^R r\,dr = \frac{R^2}{2}, \qquad \int_0^R r e^{-r^2/L^2}\,dr = \frac{L^2}{2}\left(1-e^{-R^2/L^2}\right). [/math]

Por tanto,

[math] M = \frac{S_0 \theta}{2} \left[ R^2 + \alpha L^2 \left(1-e^{-R^2/L^2}\right) \right]. [/math]

Finalmente, el volumen de sedimentos es

[math] V = \frac{M}{\rho_s}. [/math]

Usamos la expresión obtenida: [math] M=\frac{S_0\theta}{2}\left[R^2+\alpha L^2\left(1-e^{-R^2/L^2}\right)\right] [/math]

Sustituimos los valores: \(S_0=50\), \(\alpha=3\), \(L=500\), \(R=190\), \(\theta=2\pi/3\).

Entonces: [math] M=\frac{50\cdot (2\pi/3)}{2}\left[190^2+3\cdot 500^2\left(1-e^{-190^2/500^2}\right)\right] [/math]

Calculamos cada término: [math] 190^2=36100,\qquad 500^2=250000,\qquad \frac{190^2}{500^2}=0.1444 [/math] [math] e^{-0.1444}\approx 0.8652,\qquad 1-0.8652=0.1348 [/math]

Entonces: [math] M=\frac{50\cdot (2\pi/3)}{2}\left[36100+3\cdot 250000\cdot 0.1348\right] [/math] [math] 3\cdot 250000\cdot 0.1348=101100 [/math]

Así: [math] M=\frac{50\cdot (2\pi/3)}{2}\left[36100+101100\right] =\frac{50\cdot (2\pi/3)}{2}\cdot 137200 [/math]

[math] \frac{50\cdot (2\pi/3)}{2}=\frac{50\pi}{3}\approx 52.3599 [/math]

Finalmente: [math] M\approx 52.3599\cdot 137200\approx 7.17\times 10^6\;\text{kg} [/math]

El volumen de sedimentos es: [math] V=\frac{M}{\rho_s}=\frac{7.17\times 10^6}{1500}\approx 4.78\times 10^3\;\text{m}^3 [/math] y representa una fracción de la capacidad: [math] \frac{4.78\times 10^3}{4.25\times 10^8}\times 100\approx 0.00112\%. [/math]

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.

Figura 15: 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:

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:

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

Figura 16: Plano horizontal para Z = 125 hecho con el programa MATLAB

Figura 17: 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:

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:

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]

Figura 18: Divergencia en 2D hecha en MATLAB

Figura 19: 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:

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:

Y sustituyendo f(x,y) y g(z):

El cálculo del rotacional muestra que el flujo de entrada al embalse presenta vorticidad, lo que significa que el agua no entra de forma laminar, sino con tendencia a girar y formar remolinos. Las componentes horizontal y vertical del flujo generan rotación tanto en el plano horizontal como en planos verticales, produciendo torbellinos y zonas de recirculación a los lados del chorro. Esto implica que el flujo entrante induce una mezcla importante del agua, altera la distribución de velocidades y puede influir en el transporte de sedimentos y en la dinámica interna del embalse. En conjunto, los resultados indican que la entrada del río produce un patrón de circulación complejo y tridimensional dentro del embalse.

Figura 20: Representación horizontal de la divergencia hecha en MATLAB

Figura 21:Representación vertical con x=0 de la divergencia hecha en MATLAB

clear; clc; close all;

% Parámetros físicos y geométricos
v0 = 0.5; 
v1 = 0.1; 
R  = 30; 
phi = pi/6; 
Hagua = 125;

H  = 134; 
rho0 = 150;
b = 40;

% Radio de la presa a la altura del agua
rho_lim = rho0 + b*(1 - (Hagua^2)/(H^2));

%   PLANO HORIZONTAL (z = Hagua)

nx = 25;  
ny = 25;

x = linspace(-50, 50, nx);
y = linspace(-50, 50, ny);

[X, Y] = meshgrid(x, y);
Z = Hagua * ones(size(X));

% Campo en z = constante
dist2 = X.^2 + Y.^2;  
exp_term = exp(-dist2 / R^2);

Vx = v0 * cos(phi) * exp_term;
Vy = v0 * sin(phi) * exp_term;
Vz = v1 * exp_term .* sin(pi * Z / Hagua);

% Derivadas espaciales
dx = x(2) - x(1);
dy = y(2) - y(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 en el plano horizontal
omega_x = dVz_dy;
omega_y = -dVz_dx;
omega_z = dVy_dx - dVx_dy;

omega_mag = sqrt(omega_x.^2 + omega_y.^2);

figure;
pcolor(X, Y, omega_mag); 
shading interp;
colorbar;
hold on;
quiver(X, Y, omega_x, omega_y, 1.8,'k'); 
title('Rotacional en el plano horizontal (z = H_{agua})');
xlabel('x (m)'); ylabel('y (m)');
axis equal tight;
hold off;

%   PLANO VERTICAL (x = 0)

nz = 25;
ny = 25;

x0 = 0;        

z = linspace(0, Hagua, nz);
y = linspace(-50, 50, ny);

[Yv, Zv] = meshgrid(y, z);
Xv = x0 * ones(size(Yv));

% Campo de velocidades en el plano vertical
dist2v = Xv.^2 + Yv.^2;
exp_term_v = exp(-dist2v / R^2);

Vxv = v0 * cos(phi) * exp_term_v;
Vyv = v0 * sin(phi) * exp_term_v;
Vzv = v1 * exp_term_v .* sin(pi * Zv / Hagua);

% Derivadas espaciales en el plano vertical (∂/∂y y ∂/∂z)
dyv = y(2) - y(1);
dzv = z(2) - z(1);

[dVx_dy_v, dVx_dz_v] = gradient(Vxv, dyv, dzv);
[dVy_dy_v, dVy_dz_v] = gradient(Vyv, dyv, dzv);
[dVz_dy_v, dVz_dz_v] = gradient(Vzv, dyv, dzv);

% Rotacional completo:

omega_x_v = dVz_dy_v - dVy_dz_v;
omega_y_v = dVx_dz_v;          
omega_z_v = - dVx_dy_v;         

omega_mag_v = sqrt(omega_x_v.^2 + omega_z_v.^2);

figure;
pcolor(Yv, Zv, omega_mag_v);
shading interp;
colorbar;
hold on;

% Campo de rotacional proyectado (y,z)
quiver(Yv, Zv, omega_z_v, omega_x_v, 1.5, 'k');

title('Rotacional en el plano vertical (x = 0)');
xlabel('y (m)'); 
ylabel('z (m)');
axis equal tight;
hold off;

7.4 Caudal de entrada

El caudal nos dice cuánta agua o líquido está fluyendo y qué tan rápido lo hace. Conocer el caudal del río es fundamental para una presa porque permite: asegurar la estructura, dimensionar correctamente la presa y sus compuertas, gestionar el agua para riego o energía, y prever riesgos como inundaciones o sequías.

Sabemos que nuestro campo de velocidad de entrada es el siguiente:

Pues tenemos como condición que [math]v_{1}=0[/math], es decir, que no tenga componente vertical.
El caudal se define con la siguiente fórmula: [math]Q=v S[/math], siendo [math]v[/math] la velocidad y [math]S[/math] la sección transversal.
Sin embargo, al no tener una sección fácil de calcular (rectángulo, circunferencia, triángulo...), usaremos la siguiente: [math]Q=\iint_{seccion}^{}v\cdot ndS[/math].

Sabemos que el caudal entra en dirección [math]\left( \cos\varphi\, \hat{\imath} + \sin\varphi\, \hat{\jmath} \right)\;[/math], entonces, un vector normal de la sección transversal. ([math]\vec n = \left( \cos\varphi\, \hat{\imath} + \sin\varphi\, \hat{\jmath} \right)\;[/math]).

Calculemos el producto [math]v\cdot n[/math]: [math]v_{0}\, e^{-\frac{x^{2}+y^{2}}{R^{2}}}\left( \cos\varphi\, \hat{\imath} + \sin\varphi\, \hat{\jmath} \right)\;\cdot\left( \cos\varphi\, \hat{\imath} + \sin\varphi\, \hat{\jmath} \right)\;=v_{0}\, e^{-\frac{x^{2}+y^{2}}{R^{2}}}\left ( cos\varphi cos\varphi +sin\varphi sin\varphi \right )=v_{0}\, e^{-\frac{x^{2}+y^{2}}{R^{2}}}\left ( cos^2\varphi+sin^2\varphi\right )=v_{0}e^{-\frac{x^{2}+y^{2}}{R^{2}}}[/math]

Como tenemos información sobre el radio, vamos a pasarlo a coordenadas cilíndricas para simplificar la integral:

Ahora resolvemos la integral, para ello haremos un cambio de variable:

8 Póster

Imagen póster

9 Bibliografía y webgrafía

Todo el código y gráficas realizadas con MATLAB.
Hecho con la asistencia de ChatGPT.
Apuntes de la asignatura de Teoría de Campos.
Apuntes de la asignatura de Física de Sólidos y Fluidos.
Tipos de presas:
https://masqueingenieria.com/blog/tipos-de-presas-y-su-clasificacion/
https://webaero.net/
https://www.iagua.es/data/infraestructuras/presas/almendra
https://aldeadavila.es/presa-de-iberdrola/
https://es.linkedin.com/pulse/la-presa-de-el-atazar-fundaci%C3%B3n-canal
https://www.cronicanorte.es/cifras-curiosidades-del-embalse-del-atazar-canal-isabel-ii/98414 Póster hecho con canva