Diferencia entre revisiones de «La Presa de El Atazar (GRUPO 14)»

De MateWiki
Saltar a: navegación, buscar
(→‎Representación de curvas de isoconcentración y análisis geométrico)
(→‎Representación de curvas de isoconcentración y análisis geométrico)
Línea 479: Línea 479:
 
%Curvas de nivel
 
%Curvas de nivel
 
figure
 
figure
−
contour (X, Y, 5, 15, 'LineWidth",1.2);
+
contour (X, Y, S, 10, 'LineWidth', 1.2);
 
colorbar;
 
colorbar;
−
axis equal;
+
axis equal
−
title('Curvas de isoconcentración de 5(x,y)');
+
title('Curvas de isoconcentración de S(x,y)');
 
xlabel('x (m)');
 
xlabel('x (m)');
 
ylabel('y (m)');
 
ylabel('y (m)');

Revisión del 21:40 6 dic 2025

Trabajo realizado por estudiantes
Título La Presa de El Atazar. Grupo 14
Asignatura Teoría de Campos
Curso 2025-26
Autores René Orbea Malla
David Holguin González
Selena Ajenjo Jiménez
Cristina Ojog Timofti
Julia Almudena Carrasco Gonzales
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura


1 Introducción

La Presa de El Atazar, situada en el río Lozoya y construida entre 1968 y 1972, es la mayor infraestructura hidráulica de la Comunidad de Madrid y una de las más relevantes de España. Su diseño de doble curvatura —arco en planta y paramento parabólico en el alzado— la convierte en un ejemplo muy interesante para el estudio matemático de la estabilidad estructural de las presas de curvatura y de la interacción entre el agua embalsada y la superficie resistente de la presa.


La presa de El Atazar


El objetivo de este trabajo es analizar varios fenómenos físicos asociados al funcionamiento de una presa de este tipo mediante herramientas de cálculo, geometría diferencial y visualización en MATLAB. Para ello se parte de un modelo paramétrico que describe la geometría del paramento aguas arriba y del campo de presión hidrostática generado por el embalse. A partir de este modelo se desarrollan distintas representaciones y cálculos que permiten estudiar la distribución de presiones, las fuerzas resultantes sobre la superficie y otros procesos relevantes del embalse, como la sedimentación o el flujo de entrada del río.


2 Representación geométrica de la presa

2.1 -Modelo paramétrico en coordenadas cilíndricas

Para estudiar la acción del agua sobre la presa se considera el paramento aguas arriba, es decir, la superficie en contacto directo con el embalse. La presa de El Atazar presenta una doble curvatura, combinando un arco circular en planta y un perfil parabólico en alzado.

Este tipo de geometría permite transmitir gran parte de los esfuerzos hacia los estribos del valle, reduciendo la presión neta soportada por el hormigón y aumentando la estabilidad global de la estructura. La superficie se describe mediante coordenadas cilíndricas \((\rho, \theta, z)\), tomando como origen un punto situado en el fondo del valle y considerando el eje \(z\) como vertical ascendente. Los parámetros recorren los intervalos:

[math]\theta \in \left[\frac{2\pi}{3}, \frac{4\pi}{3}\right],\quad z \in [0, H][/math]


donde [math]H = 134 \text{ m}[/math] es la altura total de la presa. En este sistema, la geometría del paramento se modela mediante la expresión:

[math]\rho = \rho_0 + b\left(1 - \frac{z^2}{H^2}\right)[/math]

con:

  • [math]\rho_0 = 150 \text{ m}[/math]: radio en la coronación (altura máxima)
  • [math]b = 40 \text{ m}[/math]: parámetro que controla la curvatura parabólica en el alzado


2.2 -Superficie parametrizada en MATLAB

Para representar la superficie en MATLAB se transforma la parametrización expresada en coordenadas cilíndricas a coordenadas cartesianas mediante el cambio:


[math] \begin{cases} x(\theta,z)\ = \rho(\theta,z)\,\cos\theta,\\[2mm] y(\theta,z)\ = \rho(\theta,z)\,\sin\theta,\\[1mm] z = z. \end{cases} [/math]


Este cambio de coordenadas permite generar una malla tridimensional del paramento aguas arriba, que será esencial para superponer posteriormente el campo de presión hidrostática y analizar las fuerzas ejercidas por el agua. A continuación se muestra el código utilizado para generar la superficie paramétrica en MATLAB:

Figura 1. Representación de la presa parametrizada
% Parámetros geométricos de la presa
H = 134;         
rho0 = 150;       
b = 40;           
theta_min = 2*pi/3;
theta_max = 4*pi/3;
Ntheta = 200;
Nz = 200;

% Discretización de la malla
theta = linspace(theta_min, theta_max, Ntheta);
z = linspace(0, H, Nz);
[Theta, Z] = meshgrid(theta, z);

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

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

% Representación de la superficie
figure;
surf(X, Y, Z, 'EdgeColor', 'none', 'FaceColor', [0.2 0.4 1] );
xlabel('Ancho (m)')
ylabel('Largo (m)')
zlabel('Altura (m)')
title('Superficie aguas arriba de la presa de El Atazar');
axis equal;
view(45, 30);


3 Campo de presión hidrostática

3.1 - Modelo matemático de la presión

El agua del embalse ejerce sobre el paramento aguas arriba una presión hidrostática que aumenta linealmente con la profundidad. Asumiendo que la superficie libre del agua se encuentra a una altura [math]H_{\text{agua}}[/math]respecto al fondo de la presa, la presión absoluta a una altura \(z\) medida desde la base se expresa como:


[math] P(z) = P_0 + \rho_{\text{agua}}* g \left(H_{\text{agua}} - z\right), \quad z \in [0, H_{\text{agua}}][/math]

donde:

  • [math] P_0 : \text{presión atmosférica}[/math]
  • [math] \rho_{\text{agua}} = 1000~\text{kg/m}^3[/math]
  • [math]g = 9.81~\text{m/s}^2[/math]
  • [math]H_{\text{agua}} = 125~\text{m}[/math], altura del nivel del agua en el embalse


La fórmula refleja que la presión es máxima en la base de la presa \( (z=0) \) y disminuye linealmente hasta alcanzar \(P_0\) en la superficie del embalse.


3.2 - Representación del mapa de presión sobre la superficie

Una vez obtenido el campo de presiones \(P(z)\), es posible asignar un color a cada punto del paramento aguas arriba para visualizar cómo varía la presión con la profundidad. Para ello, se utiliza la misma malla \(( \theta, z)\) construida en la sección anterior y se evalúa la expresión de \(P(z)\) sobre toda la superficie. El resultado se representa mediante un mapa de colores que permite identificar de forma inmediata las zonas sometidas a mayor presión. El siguiente código de MATLAB genera dicha representación:

Figura 2. Campo escalar de presión
% Parámetros geométricos de la presa
H = 134; 
rho0= 150;  
b = 40; 
theta = linspace(2*pi/3, 4*pi/3, 300); 
z = linspace(0, H, 300);  
[TH, Z] = meshgrid(theta, z); 

% Radio parabólico en función de z 
R = rho0 + b*(1 - (Z.^2)/(H^2)); 

% Conversión a coordenadas cartesianas 
X = R .* cos(TH); 
Y = R .* sin(TH); 
Zs = Z;

% Campo de presión hidrostática 
P0 = 101325; 
rho_agua = 1000;  
g = 9.81;
Hagua = 125;

P = P0 + rho_agua * g .* (Hagua - Zs); 

% Representación de la presión sobre la presa
figure 
surf(X, Y, Zs, P)
shading interp 
colormap jet 
colorbar 
title('Campo escalar de presión P(z) sobre la presa del Atazar') 
xlabel('x (m)') 
ylabel('y (m)') 
zlabel('z (m)') 
axis equal 
view(45, 20)

Podemos observar como en la representación gráfica del campo de presiones, se identifica claramente cómo la presión aumenta en profundidad, alcanzando sus valores máximos en la base de la presa y disminuyendo progresivamente hasta la superficie libre del embalse.


3.3 - Cálculo de la fuerza total y presión por unidad de superficie

En este apartado se calculará la fuerza total de la presión hidrostática que actúa sobre la cara aguas arriba de la presa, así como la presión por unidad de superficie (tensión media) para los dos modelos geométricos considerados:

  • Doble curvatura: [math]b = 40 \text{ m}[/math]
  • Curvatura simple: [math]b = 0 \text{ m}[/math] (cilindro)

Para, a continuación, poder comparar y analizar cuál de las dos geometrías resiste mejor la presión ejercida por el agua.


La presión hidrostática depende únicamente de la profundidad, ya que el agua transmite presión de manera isotrópica. Por lo tanto, en cualquier punto de la superficie de la presa, la presión se puede expresar como:

[math] p(z) = \rho g (H_{\text{agua}} - z) [/math]


donde \(\rho\) es la densidad del agua, \(g\) la gravedad, \(H_{\text{agua}} = 125~\text{m}\) la altura máxima del agua y \(z\) la profundidad desde la base. Este planteamiento nos permite entender que la presión será máxima en la base y mínima en la coronación.


Para integrar esta presión sobre la superficie de la presa, es necesario parametrizar geométricamente toda su cara aguas arriba, utilizando las coordenadas cilíndricas ([math]\rho, \theta, z[/math]):

[math]\rho(z) = \rho_0 + b \left(1 - \frac{z^2}{H^2}\right), \quad \rho_0 = 150~\text{m}[/math]

La derivada necesaria para calcular el área diferencial sería : [math] \quad \rho'(z) = -\frac{2 b z}{H^2}[/math]


El ángulo horizontal que abarca la presa es: [math] \quad \theta \in \left[ \frac{2\pi}{3}, \frac{4\pi}{3}\right][/math]


La superficie parametrizada, la cual vimos en el apartado 2.2, se define como:

[math] S(z,\theta) =(\rho(z)\cos\theta, \rho(z)\sin\theta, z)[/math]


El elemento del área diferencial de la superficie se obtiene a partir de las derivadas parciales:

[math] dA = \left\| \frac{\partial \mathbf{s}}{\partial z} \times \frac{\partial \mathbf{s}}{\partial \theta} \right\| \, dz\, d\theta \;\Longrightarrow\; dA = \rho(z)\,\sqrt{1 + \left(\rho'(z)\right)^2}\, dz\, d\theta [/math]



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

[math] d\overrightarrow{F} = p(z) \, dA \, \overrightarrow{n} [/math]

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, \quad \theta_1 = \frac{2\pi}{3}, \, \theta_2 = \frac{4\pi}{3} [/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]


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 = 48898,98~\text{m²}\;\Longrightarrow\;\sigma = 0,613~\text{MPa}[/math]
  • [math]b = 0~\text{m}[/math]: [math]A = 39269,19~\text{m²}\;\Longrightarrow\; \sigma = 0,619~\text{MPa}[/math]


  • Conclusión:

Aunque la fuerza total es mayor en la presa de doble curvatura, la tensión media es menor, lo que indica que la distribución de esfuerzos es más uniforme. Por lo tanto, la estructura de doble curvatura soporta mejor el empuje hidrostático, reduciendo riesgos de concentración de esfuerzos y posibles fallos estructurales.


4 Modelización y análisis de procesos en el embalse

4.1 -Sedimentación en el embalse

La sedimentación en un embalse es un proceso natural que consiste en el depósito gradual de partículas sólidas transportadas por el río. Estas partículas tienden a acumularse en el fondo, reduciendo progresivamente la capacidad de almacenamiento del embalse y afectando potencialmente su funcionamiento y eficiencia. Para analizar este fenómeno, consideramos un modelo simplificado en el que el fondo se trata como un plano horizontal en \(z=0\) y el dominio está limitado por un sector circular definido por la base de la presa.

En este modelo, la concentración de sedimentos depositados sobre el fondo se expresa en kg/m² y se define mediante la función::

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

donde:

  • [math]S_0 = 50 \text{ kg/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 \text{ m}[/math] es una escala característica.
  • [math](x,y)[/math] son coordenadas horizontales sobre el fondo del embalse.


4.1.1 Representación del mapa de sedimentos en MATLAB

El campo escalar de sedimentos [math]S(x, y)[/math] se puede visualizar mediante un mapa de colores sobre el fondo del embalse, lo que permite identificar las zonas de mayor y menor acumulación de sedimentos.

A continuación se muestra un ejemplo de código MATLAB para generar esta representación:

Figura 3. Representación de la presa parametrizada
%opcional 1 apartado 1.5
%Parámetros
S0 = 50;
alpha = 3;
L = 500;

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

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

%Figura: mapa de colores
figure('Units', 'normalized', 'Position', [0.1 0.1 0.65 0.65])
p = pcolor(X, Y, S);
set (p, 'EdgeColor', 'none');
colorbar;
xlabel('x (m)');
ylabel('y (m)');
axis equal
title('Campo escalar S(x,y) (sedimentacion en el fondo)') ;


4.1.2 Cálculo y representación del gradiente de sedimentación

Para entender cómo cambia la concentración de sedimentos en el fondo del embalse, se utiliza el gradiente del campo escalar [math]S(x,y)[/math]. El gradiente es un vector que indica:

  • La dirección en la que la concentración aumenta más rápidamente.
  • La magnitud del cambio máximo en esa dirección.

Matemáticamente, el gradiente se obtiene derivando la función [math]S(x,y)[/math] respecto a las coordenadas horizontales [math]x[/math] e [math]y[/math]:


[math] \left. \begin{array}{l} * \displaystyle \frac{\partial S}{\partial x} = -\frac{2 S_0 \alpha x}{L^2} e^{- \frac{x^2 + y^2}{L^2}} \\[0.9cm] * \displaystyle \frac{\partial S}{\partial y} = -\frac{2 S_0 \alpha y}{L^2} e^{- \frac{x^2 + y^2}{L^2}} \end{array} \right\} \nabla S = \left( \frac{\partial S}{\partial x}, \frac{\partial S}{\partial y} \right) = -\frac{2 S_0 \alpha}{L^2} e^{- \frac{x^2 + y^2}{L^2}} \left( x \overrightarrow{i} + y \overrightarrow{j} \right) [/math]


De estas derivadas se construye el vector gradiente [math]\nabla S = \left( \frac{\partial S}{\partial x}, \frac{\partial S}{\partial y} \right)[/math], cuya magnitud indica qué tan rápido cambia la concentración en un punto, y cuya dirección señala hacia dónde se incrementa más rápidamente la sedimentación. Por lo tanto se puede interpretar que :

  • Cerca del centro del embalse, donde [math]\rho = \sqrt{x^2 + y^2} [/math] es pequeña, la magnitud del gradiente es grande. Esto indica que los sedimentos tienden a acumularse rápidamente hacia el eje de entrada del río.
  • Hacia los bordes del embalse, [math]ρ[/math] aumenta, el gradiente disminuye, lo que significa que el cambio de concentración es más lento y la acumulación de sedimentos es menor.


En MATLAB, la magnitud del gradiente puede visualizarse mediante un mapa de colores, donde las zonas de color más intenso representan cambios rápidos en la concentración. La dirección del gradiente se representa con flechas usando un quiver plot, indicando hacia dónde se desplazan los sedimentos hacia las zonas de mayor acumulación.

Figura 4. Representación de la presa parametrizada
%Opcional 1 apartado 1.6

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

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

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

% Gradiente
coef = -2 * S0 * alpha / (L^2);
Sx = coef * X .* exp(-R2 / L^2);
Sy = coef * Y .* exp(-R2 / L^2);

% Magnitud del gradiente
mag = sqrt(Sx.^2 + Sy.^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('Magnitud del gradiente y campo vectorial de sedimentacion')

% Paso de muestreo para quiver
step = 8;

% Escala de flechas quiver para visualización
quiver(X(1:step:end,1:step:end), ...
            Y(1:step:end,1:step:end), ...
	    Sx(1:step:end,1:step:end), ...
	    Sy(1:step:end,1:step:end), ...
	    'k')
% Origen
plot(0,0,'r.','MarkerSize',18)

% Radio r = L/sqrt(2)
rmax = L/sqrt(2);
theta = linspace(0,2*pi,300);

% Leyenda
plot(rmax*cos(theta), rmax*sin(theta),'w--','LineWidth',1.5)
legend({'Magnitud del gradiente','Campo vectorial','Origen','Radio r=L/sqrt(2)'},
	      'Location','northeastoutside')
hold off


Como se observa en la representación generada por MATLAB, los vectores del gradiente señalan la dirección hacia la cual la concentración de sedimentos aumenta más rápidamente, apuntando mayoritariamente hacia el eje del ingreso del río. La magnitud de los vectores indica la rapidez con que cambia la concentración: cerca del eje, los vectores son más largos, reflejando un aumento rápido de sedimentos, mientras que hacia los bordes del embalse, los vectores son más cortos, indicando cambios más suaves.


4.1.3 Cálculo de la masa, volumen y porcentaje de sedimentos

Para cuantificar la cantidad de sedimentos depositados en el fondo del embalse, se integra el campo de concentración [math]S(x,y)[/math] sobre la superficie del embalse. Para simplificarlo, consideramos que el fondo del embalse es un plano horizontal [math]z = 0[/math] y que su extensión está limitada por el sector circular definido por la presa en la base. El radio de la base se obtiene sumando el radio de coronación con el parámetro de curvatura:

[math] \rho_\text{base} = \rho_0 + b = 150 + 40 = 190~\text{m} [/math]


El ángulo horizontal que abarca la presa es: [math]\theta \in \left[ \frac{2\pi}{3}, \frac{4\pi}{3}\right][/math]


El campo de sedimentos se modela como:[math] \quad S(x,y) = S_0 \left( 1 + \alpha e^{-(x^2+y^2)/L^2} \right) = S_0 \left( 1 + \alpha e^{-\rho^2/L^2} \right), [/math] donde [math]\rho^2 = x^2 + y^2[/math]


La masa total de sedimentos se calcula integrando la concentración sobre toda el área del fondo: [math]\quad M = \iint_{A} S(x,y)\, dA [/math]


Aprovechando la simetría radial del sector circular, se transforma la integral a coordenadas polares [math]( x = \rho\cos\theta, y = \rho\sin\theta) \;\Longrightarrow\; dA = \rho \, d \rho \, d\theta[/math]:


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


Resolviendo la integral respecto a [math] \rho [/math]:

[math] \quad \int_0^{\rho_\text{base}} \rho\, dr + \int_0^{\rho_\text{base}} \rho\, e^{- \rho^2/L^2} dr =\frac{\rho_\text{base}^2}{2}+ \frac{L^2}{2} \left( 1 - e^{-\rho_\text{base}^2/L^2} \right) [/math]


Por lo tanto, la masa total es:[math] \quad 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]


Sustituyendo los valores proporcionados por el enunciado:

[math] \begin{array}{lcl} \alpha &=& 3, \\ L &=& 500~\text{m}, \\ \rho_\text{base} &=& 190~\text{m}, \\ \Delta\theta &=& \dfrac{2\pi}{3} \end{array} \quad\Longrightarrow\quad M \approx 7{,}17 \times 10^6~\text{kg} [/math]


Si la densidad de los sedimentos es [math]\rho_s = 1500~\text{kg/m³}[/math], el volumen total ocupado será: [math]\quad V = \frac{M}{\rho_s} \approx \frac{7,17 \times 10^6}{1500} \approx 4780~\text{m³} [/math]


Por último, para conocer el porcentaje respecto a la capacidad total del embalse: [math]\quad V_{\text{emb}} = 425~\text{hm³} = 425 \times 10^6~\text{m³}[/math]


Entonces, el porcentaje de capacidad ocupado por sedimentos se calcula como: [math] \quad \text{Porcentaje} = \frac{V}{V_{\text{emb}}} \cdot 100 \approx 0,0011\% [/math]


En conclusión podemos decir que la masa de sedimentos depositada hasta el momento es muy pequeña comparada con la capacidad total del embalse. Por lo tanto, la sedimentación aún no afecta significativamente al volumen disponible, aunque es importante monitorizarla para el mantenimiento futuro del embalse.


4.1.4 Representación de curvas de isoconcentración y análisis geométrico

Para visualizar cómo se distribuyen los sedimentos en el fondo del embalse, se emplean las curvas de isoconcentración, que son líneas sobre el plano [math]z=0[/math] donde la concentración [math]S(x,y)[/math] permanece constante. Estas curvas permiten identificar regiones de alta y baja acumulación de sedimentos y analizar su forma geométrica.

Dado que el campo de sedimentos tiene simetría radial respecto al centro del embalse, es conveniente expresarlo en coordenadas polares:


[math] S(\rho) = S_0 \left( 1 + \alpha e^{-\rho^2/L^2} \right)[/math] , donde [math]\rho = \sqrt{x^2 + y^2}[/math]


Para un valor específico de concentración [math]S_c[/math], el radio de la curva de isoconcentración correspondiente se obtiene resolviendo:


[math] S_c = S_0 \left( 1 + \alpha e^{-\rho^2/L^2} \right) \quad \Rightarrow \quad \rho = L \sqrt{-\ln\left(\frac{S_c/S_0 - 1}{\alpha}\right)} [/math]


En coordenadas cartesianas, estas curvas se representan como circunferencias concéntricas centradas en el origen:

[math] x = \rho \cos \theta, \quad y = \rho \sin \theta, \quad \theta \in [2\pi/3, 4\pi/3]. [/math]

Interpretación geométrica:

  • Las curvas más pequeñas (con menor radio) corresponden a zonas de mayor concentración de sedimentos, generalmente cerca de la entrada del río.
  • Las curvas más grandes (mayor radio) indican que la concentración disminuye hacia los bordes del embalse.

Este análisis permite identificar visualmente las áreas donde los sedimentos tienden a acumularse, facilitando la planificación de limpiezas, dragados o modificaciones operativas del embalse para mejorar su funcionamiento y prevenir acumulaciones localizadas.

Figura 5. Representación de la presa parametrizada
%opcional 1 apartado 1.8

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

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

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

%Curvas de nivel
figure
contour (X, Y, S, 10, 'LineWidth', 1.2);
colorbar;
axis equal
title('Curvas de isoconcentración de S(x,y)');
xlabel('x (m)');
ylabel('y (m)');


4.2 -Flujo de entrada del río Lozoya

El río Lozoya entra al embalse desde el extremo norte, generando un flujo de agua que se distribuye en el sector de entrada. Para entender y analizar el comportamiento de este flujo, se define un campo vectorial de velocidad, el cual describe tanto la dirección como la magnitud del movimiento del agua en cada punto del dominio.Esto permite estudiar cómo el agua se distribuye y como transporta sedimentos o calor dentro del embalse.


4.2.1 Modelo del campo de velocidad

El campo de velocidad se expresa como:

Campo 2D (plano horizontal):[math] \quad \overrightarrow{v}(x,y) = v_0 \, e^{-\frac{x^2 + y^2}{R^2}} (\cos\phi \, \overrightarrow{i} + \sin\phi \, \overrightarrow{j})[/math]

donde:

  • [math]v_0 = 0.5 \text{ m/s}[/math] es la velocidad máxima en el centro de entrada
  • [math]R = 30 \text{ m}[/math] es el radio característico
  • [math]\phi = \pi/6[/math] es el ángulo de entrada.

El dominio considerado es un cuadrado de \( [-100\,\text{m},\, 100\,\text{m}]\) en \(x\) y \(y\), limitado por el radio de la presa a la altura del agua.


Extensión 3D (con componente vertical):[math] \quad \overrightarrow{v}(x,y,z) = v_0 \, e^{-\frac{x^2 + y^2}{R^2}} (\cos\phi \, \overrightarrow{i} + \sin\phi \, \overrightarrow{j}) + v_1 \, \sin\left(\frac{\pi z}{H_{\text{agua}}}\right) \, e^{\frac{-x^2 + y^2}{R^2}} \,\overrightarrow{k}[/math]

donde:

  • [math]v_1 = 0.1 \text{ m/s}[/math] es la velocidad máxima de la componente vertical,
  • [math]z \in [0, H_{\text{agua}}][/math] es la altura sobre la base de la presa,

El dominio en \(x\) y \(y\) se mantiene igual al del Campo 2D


4.2.2 Representación vectorial del flujo en plano horizontal y vertical

Una vez definido el campo de velocidad del flujo de entrada del río Lozoya, el siguiente paso es visualizarlo para interpretar su comportamiento dentro del embalse. La herramienta habitual para ello es la representación mediante campos vectoriales, que muestran en cada punto la dirección y magnitud del flujo.

Figura 6. Representación vectorial del flujo
%Parámetros
v0 = 0.5;
v1 = 0.1;
R = 30;
phi = pi/6; 
H = 125;

%Dominio (modifica resolución si quieres más detalle)
xvec = -100:4:100;
yvec = -100:4:100;
[x,y] = meshgrid(xvec, yvec); 

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

%Función E
E = exp(-(x.^2+y.^2)/R^2);

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

% 3.5
%representacion plano horizontal (z = H)
figure('Name', 'Campo en superficie (z = H)', 'NumberTitle', 'off');
quiver(x, y, vx, vy, 'AutoScale', 'on'); 
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, 22] meshgrid(xz, zz);
E_vert = exp(-(XZ.^2)/R^2); 

VX_vert = v0 * cos(phi) .*E_vert;
VZ_vert = v1 *  sin(pi*ZZ/H) .*E_vert;

figure( 'Name', 'Campo en plano vertical y=0', 'NumberTitle', 'off');
quiver(XZ, ZZ, VX vert, VZ vert, Autoscale', 'an');
xlabel('x (m)'); ylabel('z (m)');
title('Plano vertical (y = 0): componentes v_x (h) y v_z (v)');
axis tight;
grid on;


4.2.3 Cálculo y representación de la divergencia

La divergencia de un campo vectorial sirve para analizar cómo se comporta el flujo dentro de un sistema hidrodinámico. En este trabajo, estudiar la divergencia del campo de velocidad del río Lozoya permite identificar zonas donde el agua tiende a acumularse, dispersarse o mantenerse en equilibrio local.

La divergencia se define como la suma de las derivadas parciales de cada componente del vector velocidad respecto a su coordenada correspondiente: [math] \quad \nabla \cdot \overrightarrow{v}= \frac{\partial v_x}{\partial x}+ \frac{\partial v_y}{\partial y}+ \frac{\partial v_z}{\partial z}[/math]

El campo de velocidad que tomaremos como referencia será : [math] \quad \overrightarrow{v}(x,y,z) = v_0 \, e^{-\frac{x^2 + y^2}{R^2}} (\cos\phi \, \overrightarrow{i} + \sin\phi \, \overrightarrow{j}) + v_1 \, \sin\left(\frac{\pi z}{H_{\text{agua}}}\right) \, e^{\frac{-x^2 + y^2}{R^2}} \, \overrightarrow{k} [/math]


Para ayudarnos en los cálculos definimos: [math]\quad E = e^{-\frac{x^2 + y^2}{R^2}} [/math]

Entonces:

[math] \begin{align*} v_x &= v_0 \cos\phi \, E &\Longrightarrow\quad& \frac{\partial v_x}{\partial x} = v_0 \cos\phi \, E \left(-\frac{2x}{R^2}\right) \\ v_y &= v_0 \sin\phi \, E &\Longrightarrow\quad& \frac{\partial v_y}{\partial y} = v_0 \sin\phi \, E \left(-\frac{2y}{R^2}\right) \\ v_z &= v_1 \sin\left(\frac{\pi z}{H_{\text{agua}}}\right) E &\Longrightarrow\quad& \frac{\partial v_z}{\partial z} = v_1 \frac{\pi}{H_{\text{agua}}} \cos\left(\frac{\pi z}{H_{\text{agua}}}\right) E \end{align*} [/math]


Sustituyendo:

[math] \nabla\cdot\vec{v}= v_0\cos\phi\,E\left(-\frac{2x}{R^2}\right) + v_0\sin\phi\,E\left(-\frac{2y}{R^2}\right) + v_1\frac{\pi}{H_{\text{agua}}} \cos\left(\frac{\pi z}{H_{\text{agua}}}\right) E [/math]

Factorizando:

[math] \nabla \cdot \vec{v} = E\left[ -\frac{2v_0}{R^2}\left(x\cos\phi + y\sin\phi\right) + v_1\frac{\pi}{H_{\text{agua}}} \cos\left(\frac{\pi z}{H_{\text{agua}}}\right) \right] = e^{-\frac{x^2 + y^2}{R^2}} \left[ -\frac{2v_0}{R^2}\left(x\cos\phi + y\sin\phi\right) + v_1\frac{\pi}{H_{\text{agua}}} \cos\left(\frac{\pi z}{H_{\text{agua}}}\right) \right] [/math]


El valor máximo de la divergencia se alcanza en las condiciones en las que tanto el factor dentro del corchete como el término exponencial son máximos. Dado que [math] e^{-(x^2+y^2)/R^2}[/math], el exponencial alcanza su valor máximo en [math] x = y = 0 [/math], lo cual conlleva que se encuentra en el centro.

Desde el punto de vista físico, la divergencia proporciona información sobre la presencia de fuentes y sumideros locales en el flujo: una divergencia positiva indica que el fluido se aleja de esa zona (fuente), mientras que una divergencia negativa indica convergencia del flujo hacia esa región (sumidero). En este caso particular, la divergencia máxima en el centro de la superficie sugiere que el fluido se eleva verticalmente desde el centro hacia arriba, con una magnitud que depende de \(v_1\) y de la altura total del agua.


Figura 7. Representación de la divergencia
% 3.6: Divergencia en el plano superficial z = H

 coeff x = 2*v0*cos(phi)/R^2;
 coeff y = 2*v0*sin(phi)/R^2;
 coeff_zterm=v1 * (pi/ H);

 divergence = E .* (coeff_x  .* x + coeff_y .* y + coeff_zterm .* cos(pi*z_surface / H) );

 figure('Name', 'Divergencia en superficie', 'NumberTitle', 'off');
 surf(x, y, divergence);
 shading interp; hold on;
 contour3(x,y, divergence, 12, 'k'); hold off;
 xlabel('x (m)'); ylabel('y (m)'); zlabel('div (1/s)');
 title('Divergencia (plano superficial z = H)');
 colorbar;
 view(45,30);
 axis tight;

% Representación en mapa de colores 2D (vista desde arriba)

figure('Name', 'Mapa de colores divergencia (superficie)', 'NumberTitle', 'off');
imagesc(xvec, yvec, divergence);
axis xy; axis equal;
xlabel('x (m)'); ylabel('y (m)');
title('Mapa de colores: divergencia en superficie (z = H)');
colorbar;


4.2.4 Cálculo y representación del rotacional

El rotacional del campo de velocidad matemáticamente describe la tendencia del flujo a girar o generar vorticidad. En el contexto del río Lozoya entrando por el norte, el rotacional permite medir cuánto gira el agua al entrar al embalse. Se define como:

[math]\nabla \times \overrightarrow{v} = \begin{vmatrix} \overrightarrow{i} & \overrightarrow{j} & \overrightarrow{k} \\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \\ v_x & v_y & v_z \end{vmatrix} [/math]

Las componentes son:

[math] v_x= v_0 \cos \phi \, e^{-\frac{x^2+y^2}{R^2}}\qquad \ v_y= v_0 \sin \phi \, e^{-\frac{x^2+y^2}{R^2}} \qquad \ v_z= v_1 \sin\Big(\frac{\pi z}{H_{\text{agua}}}\Big)\, e^{-\frac{x^2+y^2}{R^2}} [/math]


Para ayudarnos en los cálculos definimos: [math]\quad E = e^{-\frac{x^2 + y^2}{R^2}} [/math]


Componente \( x \): ​

  • [math] (\nabla \times \overrightarrow{v})_x = \frac{\partial v_z}{\partial y} - \frac{\partial v_y}{\partial z} \quad \Longrightarrow \quad \frac{\partial v_z}{\partial y} = v_1\sin\Big(\frac{\pi z}{H_{\text{agua}}}\Big) \;E\left(-\frac{2y}{R^2}\right) ,\quad \frac{\partial v_x}{\partial z} = 0 \quad \Longrightarrow \quad (\nabla \times \overrightarrow{v})_y = v_1\sin\Big(\frac{\pi z}{H_{\text{agua}}}\Big)\;E\left(-\frac{2y}{R^2}\right) [/math]


Componente \( y \): ​

  • [math] (\nabla \times \overrightarrow{v})_y = \frac{\partial v_x}{\partial z} - \frac{\partial v_z}{\partial x} \quad \Longrightarrow \quad \frac{\partial v_x}{\partial z} = 0, \quad \frac{\partial v_z}{\partial x} =v_1\sin\Big(\frac{\pi z}{H_{\text{agua}}}\Big)\,E \left(-\frac{2x}{R^2}\right) \quad \Longrightarrow \quad (\nabla \times \overrightarrow{v})_y = v_1\sin (\frac{ \pi z}{H_{\text{agua}} })\, E\left(\frac{2x}{R^2}\right) [/math]


Componente \( z \): ​

  • [math] (\nabla \times \overrightarrow{v})_z = \frac{\partial v_y}{\partial x} - \frac{\partial v_x}{\partial y} \quad \Longrightarrow \quad (\nabla \times \overrightarrow{v})_z = \frac{2y}{R^2} v_0 \cos\phi \, E - \frac{2x}{R^2} v_0 \sin\phi \, E [/math]


Derivadas de \( v_y \) y \( v_z \):

  • [math]\frac{\partial v_y}{\partial x} = - \frac{2x}{R^2} v_0 \sin\phi \, E[/math]
  • [math]\frac{\partial v_x}{\partial y} = - \frac{2y}{R^2} v_0 \cos\phi \, E[/math]


El resultado rotacional sería: [math] \nabla \times \overrightarrow{v} = \begin{pmatrix} - \frac{2}{R^2} \, y \, v_1 \, \sin\Big(\frac{\pi z}{H_{\text{agua}}}\Big) \, E \\[1em] \frac{2}{R^2} \, x \, v_1 \, \sin\Big(\frac{\pi z}{H_{\text{agua}}}\Big) \, E\\[1em] \frac{2}{R^2} \, y \, v_0 \, \cos\phi \, E - \frac{2}{R^2} \, x \, v_0 \, \sin\phi \, E \end{pmatrix} [/math]


Figura 8. Representación del Rotacional
% 3.8: Rotacional en plano superficial z = H

%coef numericos

c1 = 2 * v1 / R^2;
c2 = 2 * v0 / R^2;

curl_x = -c1 .* y .* sin(piz_surface / H) .* Ε;
curl_y = -c1 .* x .* sin(piz_surface / H) .* Ε;
curl_z = -c2 .* E .* (x.*sin(phi) - y.*cos(phi));

% Parte horizontal del rotacional

figure('Name', 'Rotacional (horizontal) en superficie', 'NumberTitle', 'off');
quiver(x, y, curl_x, curl_y, 'AutoScale', 'on');
axis equal;
xlabel('x (m)'); ylabel('y (m)');
title('Rotacional (componentes horizontales) en superficie');
grid on;

% Componente vertical del rotacional en mapa de colores

figure('Name', 'Componente z del rotacional', 'NumberTitle', 'off');
surf (x, y, curl_z);
shading interp; colorbar;
xlabel('x (m)'); ylabel('y (m)'); zlabel('(\omega_z)');
title('Componente vertical del rotacional \omega_z en superficie (z=H)');
view(45,30);
axis tight;


Las componentes x e y del rotacional dependen de la parte vertical oscilatoria del flujo (la componente v z). La componente z del rotacional depende únicamente del flujo horizontal y describe la aparición de un giro en el plano. El factor E =atenúa la vorticidad rápidamente al alejarnos del punto de entrada del río, lo que concuerda con el comportamiento físico esperado.


4.2.5 Cálculo del caudal de entrada mediante integración del flujo

El caudal es la cantidad de volumen de agua que atraviesa una superficie por unidad de tiempo. Matemáticamente, el caudal \(Q\) a través de una superficie \(S\) perpendicular al flujo se define como:

[math] Q = \iint_S \overrightarrow{v}\cdot \overrightarrow{n}=\, dA [/math]

Donde:

  • [math]\overrightarrow{v}[/math] es el campo de velocidad
  • [math]\overrightarrow{n}[/math] es el vector normal a la superficie
  • [math]dA[/math]es el elemento diferencial de área.


Para este apartado, se indica que se debe considerar \(v_1=0\), es decir, se anula la componente vertical del flujo. De este modo, el campo de velocidad queda reducido a un flujo completamente horizontal: [math] \quad \overrightarrow{v}(x,y) = v_0 \, e^{-\frac{x^2 + y^2}{R^2}} (\cos\phi \, \overrightarrow{i} + \sin\phi \, \overrightarrow{j})[/math]


Como el río entra por el norte, modelizamos la entrada mediante una superficie perpendicular a la dirección del flujo, cuya dirección viene dada por: [math] \overrightarrow{u} = (\cos\phi, \sin\phi) [/math] [math]\quad \phi = \frac{\pi}{6}[/math]


Si la superficie es perpendicular al flujo, su vector normal [math]\overrightarrow{n}[/math]será precisamente el vector [math]\overrightarrow{u}[/math]. Para el cálculo teórico, se toma una superficie plana lo suficientemente grande para que contenga la región donde el flujo es significativo, es decir, donde el término exponencial no es despreciable: [math]x,y ∈[−100,100][/math]


Como la superficie es perpendicular al flujo: [math]\overrightarrow{v}\cdot\overrightarrow{u}= v_0 e^{-\frac{x^2 + y^2}{R^2}} [/math]


Así el caudal se expresa como: [math] Q = \iint_S v_0 e^{-\frac{x^2 + y^2}{R^2}}\, dA. [/math]

Aunque el dominio es un cuadrado, el campo es radial, por lo que resulta más natural utilizar coordenadas polares:

[math] x = r\cos\theta,\qquad y = r\sin\theta,\qquad dA = r\,dr\,d\theta [/math]


Como el flujo solo existe de manera significativa en una región circular de radio máximo unos pocos múltiplos de R=30 m, se integrará sobre un disco grande donde el flujo prácticamente ha decaído:

[math] 0 \le r \le 100,\qquad 0 \le \theta \le 2\pi [/math]


Sustituyendo obtenemos la integral del caudal:

[math] Q = \int_0^{2\pi} \int_0^{100} v_0\, e^{-r^2/R^2}\, r\, dr\, d\theta \quad \Longrightarrow \quad Q = v_0 \int_0^{2\pi} d\theta \int_0^{100} r\, e^{-r^2/R^2}\, dr [/math]


Para resolverla aplicaremos un cambio de variable : [math] \int_0^{100} r\, e^\frac{-r^2}{R^2}\, dr \quad \Longrightarrow \quad u = -\frac{r^2}{R^2} ,\quad du = -\frac{2r}{R^2} dr \quad \Longrightarrow \quad\int_0^{100} r\, e^\frac{-r^2}{R^2}dr = \frac{R^2}{2} \left[ 1 - e^{-100^2/R^2} \right][/math]


Se puede observar que el término exponencial es extremadamente pequeño por lo que se permite despreciarla:[math] \quad\ e^\frac{-100^2}{30^2} \approx 1.5\times 10^{-5}. [/math]


Nos queda: [math] \int_0^{100} r\, e^{-r^2/R^2}\, dr \approx \frac{R^2}{2} [/math]


Volvemos a sustituir : [math]Q \approx v_0 \int_0^{2\pi} d\theta \, \left(\frac{R^2}{2}\right) \quad \Longrightarrow \quad Q \approx v_0 \left(\frac{R^2}{2}\right) (2\pi) \quad \Longrightarrow \quad Q = v_0 \pi R^2 [/math]


Tomando los valores del enunciado:[math]\quad v_0 = 0.5~\text{m/s} \qquad R = 30~\text{m}. [/math]


Y finalmente obtenemos el valor del caudal : [math] Q = 0.5 \cdot \pi \cdot 30^2 \quad \Longrightarrow \quad Q \approx 1413,72\, \text{m}^3/\text{s} [/math]


El valor obtenido es bastante alto, y esto se debe principalmente a que el modelo que estamos utilizando es una idealización del flujo real. En concreto, estamos suponiendo un perfil gaussiano bastante amplio y con una velocidad relativamente grande en toda la zona de entrada, lo que aumenta mucho el caudal calculado. Es importante remarcar que este resultado no representa el caudal real del río Lozoya. Lo que estamos obteniendo es el caudal asociado al modelo matemático simplificado que plantea la práctica. Por tanto, su utilidad no está en describir el comportamiento real del río, sino en comprobar cómo funciona el modelo y si es coherente con las hipótesis que hemos hecho.


5 Clasificación de las presas y ejemplos en España

5.1 - Según su material de construcción

  • Presas de hormigón: Sus buenas prestaciones en lo que a durabilidad, resistencia, impermeabilidad y facilidad de construcción se refiere, junto a un precio relativamente bajo las convierten en las más utilizadas. Dentro de las presas de hormigón, este puede ser de tipo convencional o de consistencia seca, que permite ser compactado con rodillos.
  • Presas de mampostería: Normalmente son presas más pequeñas, que emplean principalmente una estructura de piedra, arena y cemento. Su construcción suele ser perpendicular a las cárcavas (socavones en el terreno) para así reducir la velocidad del escurrimiento del agua al formar escalones que reducen la erosión.
  • Presas de materiales sueltos: Debido a su gran versatilidad y bajo coste son las más empleadas en países en desarrollo. Cuentan con un relleno de tierras sin cementar, normalmente piedras y gravas, que aportan la resistencia necesaria ante el empuje de las aguas. A esto se añaden pantallas impermeables de otros materiales.


Presa de hormigón
Presas de mampostería
Presas de materiales sueltos


5.2 - Según su forma

Tipos de presas

  • Presas de gravedad:


El mecanismo resistente de este tipo de presas es el rozamiento del cuerpo de presa con el terreno sobre el que se apoya debido a su gran peso. Para evitar posibles deslizamientos son construidas 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. Al depender de la orografía, se requiere de un suelo que sea muy estable. Su diseño es triangular, ya que es más eficiente en cuanto a resistencia y desde el punto de vista económico. La anchura de la base suele ser alrededor de un 80% de la altura, y pueden ser de eje recto o curvo. Algunos ejemplos en España son las presas de Torrejón, Villalcampo, Saucelle o Castro.


  • Presas arco:


La curvatura de estas presas resiste el empuje del agua. Es importante que se emplee un material de muy alta resistencia en los laterales de la presa, ya que es allí donde se produce el mayor esfuerzo. Por eso, este tipo de presas está limitado por las condiciones topográficas (cuanto más simétrica mejor) y orográficas del terreno. A su vez, dentro de las presas arco existen dos tipos más: Las presas bóveda o de doble arco. Su forma curva en planta y alzado trasmite el esfuerzo a los laterales y al fondo. Presas arco-gravedad, que combinan elementos de la presa arco y la de gravedad. Es curva pero su peso ayuda a los estribos a resistir los esfuerzos. Grandes presas en España como las de Cedillo o la de Aldeadávila son de arco-gravedad.


  • Presas de contrafuertes:


Este tipo de presas presenta similitudes con las de gravedad en su mecanismo de resistencia, pero cuenta con una serie de contrafuertes para ofrecer estabilidad frente al deslizamiento y al vuelco. Se emplea así una cantidad menor de material que en las presas de gravedad, pero tienen una mayor complejidad técnica. La presa de José María de Oriol, en el río Tajo, también conocida como “La Catedral”, es de tipo contrafuertes.


5.3 - Las presas españolas mas relevantes :

Presa de Almendra: es una presa de bóveda o también conocida como doble arco y es la más alta de España llegando a alcanzar una altura de 202 metros.
Presa de Aldeadávila: es una presa de hormigón del tipo arco-gravedad, ubicado en el Parque Natural de Arribes del Duero.
Presa de Ricobayo: es una presa de gravedad construida con hormigón. Su función es retener el agua del río Esla para la generación de energía hidroeléctrica, a través de las centrales que están a pie de presa y de forma subterránea.


Presa de Almendra
Presa de Aldeadávila
Presa de Ricobayo

6 Bibliografía / Referencias

Obtenido de «https://mat.caminos.upm.es/w/index.php?title=La_Presa_de_El_Atazar_(GRUPO_14)&oldid=102718»