Revisión del 01:20 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.

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.

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 2: 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 2, 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 Curvatura

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:

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

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

Presas de gravedad

Las presas de gravedad se caracterizan por poseer un cuerpo macizo, normalmente de hormigón, cuyo propio peso es el encargado de resistir la presión del agua (resistencia al deslizamiento). Su estabilidad depende fundamentalmente de la solidez de los cimientos, por lo que requieren terrenos de roca firme y homogénea. Para evitar el vuelco, la resultante de los empujes del agua y el peso propio debe estar contenida en la base del cuerpo de presa. Se construyen con hormigón en masa prácticamente en su totalidad, armándose únicamente en puntos concretos sometidos a fuertes tracciones como las galerías. Son presas que trabajan a compresión, por lo que las tracciones se deben controlar cuidadosamente, siendo el pie de aguas arriba uno de los puntos más problemáticos.

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

Por otro lado, las presas de arco o bóveda aprovechan la curvatura de su estructura para transmitir la mayor parte del empuje hacia los laterales del valle. Esto permite reducir significativamente la cantidad de hormigón empleada, lo que convierte a estas presas en estructuras muy eficientes. El arco resulta ser antifunicular de la carga radial repartida uniformemente, es decir, un arco sometido a este tipo de cargas trabaja únicamente a axil. Dicho arco transmite esos esfuerzos de compresión a los estribos de la cerrada, por lo que estos deben tener gran resistencia.

Para que sean viables, es imprescindible que los estribos sean de roca muy resistente y que el valle sea estrecho o en forma de garganta, de manera que la curvatura trabaje adecuadamente. Las presas de arco presentan un diseño más complejo y una mayor exigencia en cuanto a calidad geológica, pero a cambio requieren menos material y pueden alcanzar grandes alturas con estructuras comparativamente delgadas. 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

También existen las presas de contrafuertes . Este tipo de presas poseen un mecanismo resistente similar al de las presas de gravedad. Se trata de presas en las que la “cabeza” está inclinada y apoyadas sobre unos elementos más largos denominados “colas” que se colocan aguas abajo sujetando la cabeza a lo largo del cuerpo de presa.

A diferencia de las de gravedad, carecen prácticamente de supresión, lo que junto a la contribución del peso de la cuña de agua sobre la cabeza, hace que tengan mucha menos masa de hormigón. En cuanto a las dimensiones, el ancho de la base es del orden de la altura de presa. Sin embargo, la complejidad de su forma, a diferencia de las presas de gravedad convencionales, hace que requieran de gran cantidad de mano de obra por lo que pese a tener menos material resultan mucho más caras. Hoy en día prácticamente no se construyen.

5.2 Presas de materiales sueltos

Las presas de materiales sueltos son presas muy versátiles que se construyen prácticamente con cualquier material, por lo que son las más abundantes en el mundo. Tienen sección trapezoidal y son mucho menos esbeltas que las presas de fábrica, siendo su principal característica la zonificación de sus materiales, es decir, cada tipo de material se coloca donde mejor ejerce su función.

5.3 Ejemplos en España

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.

La Presa de Ricobayo, en el río Esla, es un ejemplo clásico de presa de gravedad. El valle en esta zona es más ancho y dispone de cimientos graníticos adecuados para soportar el peso de un muro masivo. Su diseño es más sencillo y robusto que el de una presa de arco, aunque su construcción requiere un volumen mucho mayor de material. Representa la solución óptima cuando las condiciones geológicas del fondo del valle son favorables, pero los estribos laterales no son ideales para transmitir grandes cargas.

La Presa de Aldeadávila, en el río Duero, es una de las obras hidráulicas más importantes de España y se considera una presa de arco-gravedad. Se ubica en los Arribes del Duero, un desfiladero estrecho y rocoso donde una presa de arco sería adecuada, pero la configuración del valle y la necesidad de una gran estabilidad llevaron a emplear una tipología mixta. De este modo, parte de la presión se transmite por curvatura a los estribos y parte se resiste mediante el peso propio del muro. Este tipo de presa combina la eficiencia estructural de una presa de arco con mayor robustez y versatilidad.

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

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.

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. \text{donde } [/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]

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:

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

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:

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]

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:

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.

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: