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

Revisión del 21:26 4 dic 2025

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 alzado— la convierte en un ejemplo destacado 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.

El objetivo principal de este trabajo es analizar, mediante herramientas de cálculo, geometría diferencial y visualización en MATLAB, distintos fenómenos físicos asociados al funcionamiento de una presa de estas características. Para ello se parte de un modelo geométrico paramétrico del paramento aguas arriba y del campo de presión hidrostática generado por el agua embalsada.

2 Representación geométrica de la presa

2.1 -Modelo paramétrico en coordenadas cilíndricas

Para realizar el siguiente análisis se considerará la superficie de la presa en su cara aguas arriba, es decir, la zona en contacto directo con el agua del embalse. La presa presenta una doble curvatura, es decir, una curvatura en planta (arco circular) y otra en alzado (arco parábolico). Este tipo de geometría optimiza la transmisión de esfuerzos hacia los estribos del valle, reduciendo notablemente la presión neta soportada. Para describir la superficie aguas arriba de la presa se emplean coordenadas cilíndricas \((\rho, \theta, z)\), tomando como origen el fondo del valle y considerando el eje vertical \(z\). Los parámetros recorren los intervalos:

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:

donde:

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

2.2 -Superficie parametrizada en MATLAB

Para representar la superficie en MATLAB se transforma la parametrización cilíndrica a coordenadas cartesianas mediante:

Este cambio de coordenadas permite generar una malla tridimensional que precisa del paramento mojado, y el cuál servirá posteriormente para superponer el campo de presiones y analizar las fuerzas hidrostáticas.

A continuación se muestra el código utilizado para generar la superficie paramétrica:

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 una presión hidrostática que aumenta linealmente con la profundidad. Para una columna de agua de altura [math]H_{\text{agua}}[/math], y la presión absoluta a una altura \(z\) medida desde el fondo se expresa como:

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]

3.2 - Mapa de colores sobre la superficie

Una vez calculado P(z), se puede asignar un color a cada punto de la superficie para visualizar la distribución de presión. El código MATLAB que parte de la superficie construida en la sección anterior:

Figura 2. Campo escalar de presión

% Parámetros geométricos de la presa
H = 134; 
rho0= 150;  
b = 40; 
theta = linspace(2*013, 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.A2)/(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, la presencia de un gradiente vertical pronunciado. En las zonas situadas a mayor profundidad, la presión alcanzada es significativamente superior, lo que se refleja en la aparición de colores cálidos, fundamentalmente amarillos y rojos. Por otro lado, las regiones próximas a la coronación de la presa muestran valores mínimos de presión, representados mediante tonalidades frías, principalmente azules. En cuanto a la variación horizontal, podemos decir que es prácticamente nula debido a que la presión hidrostática únicamente depende de la coordenada vertical \(z\); es decir, está determinada exclusivamente por la profundidad de la columna de agua.

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

En este aparado 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 cual sería la que mejor resiste la presión ejercida por el agua.

Para calcular cómo actúa la presión hidrostática sobre la presa, es necesario considerar primero la distribución de la presión. Sabemos que el agua transmite presión de manera isotrópica y que ésta depende únicamente de la profundidad. Por lo tanto, en cualquier punto de la superficie de la presa, la presión se puede expresar como \; [math]p(z) = ρ g (H − z)[/math], donde \(ρ\) es la densidad del agua, \(g\) la gravedad, \(H\) 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, y que no depende de la coordenada horizontal.

Para integrar esta presión sobre la superficie de la presa, es necesario parametrizar geométricamente toda su cara aguas arriba.

Para poder analizar la acción del agua sobre la presa, es fundamental considerar que la fuerza total como la distribución de tensiones sobre la superficie. Se utilizan las coordenadas cilíndricas \((\rho, \theta, z)\) para poder modelizar la geometría de la presa, siendo:

• \(z\): altura desde el fondo \((z=0)\) hasta la coronación \((z=H)\)
• [math]θ[/math]: ángulo horizontal sobre el arco de la presa
• [math]ρ[/math]: distancia al eje del valle

La curvatura de la presa se describe con la ecuación de apartado 2:

Sabiendo que [math]b = 40 \text{ m}[/math] corresponde a la doble curvatura, y [math]b = 0 \text{ m}[/math] a la curvatura simple (cilíndrica). El ángulo total de la presa es [math]Δθ=2π/3.[/math]

La superficie del paramento aguas arriba se puede parametrizar en coordenadas cilíndricas como: [math]\mathbf{S}(\theta, z) = (\rho(z)\cos\theta,, \rho(z)\sin\theta,, z)[/math] con derivadas parciales respecto a [math]z[/math] y [math]\theta[/math]: [math]\frac{\partial \mathbf{S}}{\partial z} = (\rho'(z)\cos\theta,, \rho'(z)\sin\theta,, 1),[/math] [math]\frac{\partial \mathbf{S}}{\partial \theta} = (-\rho(z)\sin\theta,, \rho(z)\cos\theta,, 0)[/math] donde la derivada radial es: [math]\rho'(z) = -\frac{2bz}{H^2}[/math]

El elemento diferencial de área se obtiene a partir del módulo del producto vectorial de las derivadas parciales: [math]\mathrm{d}A = \left|\frac{\partial \mathbf{S}}{\partial z} \times \frac{\partial \mathbf{S}}{\partial \theta}\right| , \mathrm{d}\theta , \mathrm{d}z[/math] Evaluando el producto vectorial se obtiene: [math]\mathrm{d}A = \rho(z)\sqrt{1 + (\rho'(z))^2}, \mathrm{d}\theta , \mathrm{d}z[/math] Este elemento diferencial es el que se utiliza posteriormente para calcular la fuerza total y la tensión media sobre la presa.

El vector de fuerza de presión sobre la superficie se obtiene como: [math]\overrightarrow{F} = −P(z) \overrightarrow{n}[/math], donde[math] \overrightarrow{n}[/math] es el vector normal a la superficie (apuntando hacia el agua).

La presa de doble curvatura distribuye la presión de manera más uniforme a lo largo de toda la superficie, reduciendo picos de carga local y aumentando la estabilidad estructural. Por ello, la configuración con [math]b = 40 \text{ m}[/math] es más eficiente y soporta mejor la presión que la curvatura simple.

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

4.1 -Sedimentación en el embalse

La sedimentación es un proceso natural que afecta gradualmente la capacidad de almacenamiento del embalse, ya que los sedimentos transportados por el río Lozoya tienden a acumularse en el fondo del embalse de El Atazar. Para analizar este fenómeno, consideramos un modelo simplificado de la concentración de sedimentos depositados sobre el fondo, expresada en kg/m²:

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.

El fondo se aproxima a un plano horizontal [math]z = 0[/math] y con forma de sector circular centrado en el valle, delimitado por el arco de la presa en su base.

4.1.1 Representación del mapa de sedimentos en MATLAB

El campo escalar [math]S(x, y)[/math] se puede visualizar como un mapa de colores sobre el fondo del embalse. El código creado en MATLAB sería:

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;
axis equal;
xlabel('x (m)');
ylabel('y (m)');
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 define como: [math]\nabla S(x, y) = \left( \frac{\partial S}{\partial x}, \frac{\partial S}{\partial y} \right)[/math] donde las componentes son las derivadas parciales de [math]S[/math] respecto a las coordenadas horizontales [math]x[/math] e [math]y[/math].

Dado el modelo de sedimentación: [math]S(x, y) = S_0 \left( 1 + \alpha e^{- \frac{x^2 + y^2}{L^2}} \right)[/math] [math]S_0[/math] es la sedimentación base (kg/m²). [math]\alpha[/math] modela la mayor acumulación de sedimentos cerca de la entrada del río. [math]L[/math] es una escala característica (m). Las derivadas parciales se calculan paso a paso: Derivada parcial respecto a x: [math]\frac{\partial S}{\partial x} = S_0 \cdot \alpha \cdot \frac{\partial}{\partial x} \left( e^{- \frac{x^2 + y^2}{L^2}} \right)[/math] Usando la regla de la cadena: [math]\frac{\partial}{\partial x} \left( e^{- \frac{x^2 + y^2}{L^2}} \right) = e^{- \frac{x^2 + y^2}{L^2}} \cdot \frac{\partial}{\partial x} \left( -\frac{x^2 + y^2}{L^2} \right) = -\frac{2x}{L^2} e^{- \frac{x^2 + y^2}{L^2}}[/math] Por lo tanto: [math]\frac{\partial S}{\partial x} = -\frac{2 S_0 \alpha x}{L^2} e^{- \frac{x^2 + y^2}{L^2}}[/math] Derivada parcial respecto a y: De manera análoga: [math]\frac{\partial S}{\partial y} = -\frac{2 S_0 \alpha y}{L^2} e^{- \frac{x^2 + y^2}{L^2}}[/math]

La magnitud del gradiente, que indica el cambio máximo de concentración, se obtiene como: [math]|\nabla S(x, y)| = \sqrt{\left(\frac{\partial S}{\partial x}\right)^2 + \left(\frac{\partial S}{\partial y}\right)^2} = \sqrt{\left(-\frac{2 S_0 \alpha x}{L^2} e^{- \frac{x^2 + y^2}{L^2}} \right)^2 + \left(-\frac{2 S_0 \alpha y}{L^2} e^{- \frac{x^2 + y^2}{L^2}} \right)^2}[/math] Factorizando términos comunes: [math]|\nabla S(x, y)| = \frac{2 S_0 \alpha}{L^2} \sqrt{x^2 + y^2} , e^{- \frac{x^2 + y^2}{L^2}}[/math]

La magnitud del gradiente nos da la rapidez con que varía la concentración: [math]|\nabla S(x,y)| = \sqrt{\left(\frac{\partial S}{\partial x}\right)^2 + \left(\frac{\partial S}{\partial y}\right)^2}[/math] Para entender cómo cambia la concentración de sedimentos en el fondo del embalse, se utiliza el gradiente del campo escalar S(x,y). El gradiente es un vector que indica: Dirección en la que la concentración aumenta más rápidamente. Magnitud del cambio máximo en esa dirección. Matemáticamente, el gradiente se define como:

son las derivadas parciales respecto a las coordenadas horizontales x e y. El gradiente indica la dirección de máximo incremento de sedimentos y puede representarse con vectores:

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')
step = 8;
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')
plot(0,0,'r.','MarkerSize',18)
rmax = L/sqrt(2);
theta = linspace(0,2*pi,300);
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 puede observar en la imagen dada por MATLAB, los vectores apuntan hacia las zonas con mayor concentración de sedimentos, generalmente cerca de la entrada del río.

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 simplificar, 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 abarca [math]\theta \in [2\pi/3, 4\pi/3][/math]. El campo de sedimentos se modela como: [math] 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] 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 = r\cos\theta, y = r\sin\theta, dA = r,dr,d\theta)[/math]:

[math] M = \int_{\theta_0}^{\theta_1} \int_0^{\rho_\text{base}} S(r)\, r \, dr \, d\theta = \int_{\theta_0}^{\theta_1} \int_0^{\rho_\text{base}} S_0 \left( 1 + \alpha e^{-r^2/L^2} \right) r \, dr \, d\theta [/math] Resolviendo la integral respecto a [math]r[/math]: [math] \int_0^{\rho_\text{base}} r\, dr = \frac{\rho_\text{base}^2}{2}, \quad \int_0^{\rho_\text{base}} r\, e^{-r^2/L^2} dr = \frac{L^2}{2} \left( 1 - e^{-\rho_\text{base}^2/L^2} \right) [/math] Por lo tanto, la masa total es: [math] M = S_0 \Delta \theta \left[ \frac{\rho_\text{base}^2}{2} + \frac{\alpha L^2}{2} \left( 1 - e^{-\rho_\text{base}^2/L^2} \right) \right] [/math] Sustituyendo los valores proporcionados por el enunciado ([math]S_0 = 50~\text{kg/m²}, \alpha = 3, L = 500~\text{m}, \rho_\text{base} = 190~\text{m}, \Delta\theta = 2\pi/3[/math]): [math] M \approx 7,17 \times 10^6~\text{kg}. [/math]

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

Por último, para conocer el \textbf{porcentaje respecto a la capacidad total del embalse} es [math]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] \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

Las curvas de isoconcentración representan las líneas en el fondo del embalse donde la concentración de sedimentos es constante, es decir, donde [math]S(x,y) = \text{constante}[/math]. Estas curvas permiten visualizar cómo se distribuyen los sedimentos y facilitan la identificación de zonas de mayor acumulación. Dado que el campo de sedimentos tiene simetría radial respecto al centro del embalse, se puede expresar 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 obtener las curvas de isoconcentración para un valor dado [math]S_c[/math], se resuelve la ecuación: [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] Esto proporciona el radio de la curva de isoconcentración correspondiente a cada concentración específica [math]S_c[/math]. Luego, en coordenadas cartesianas, las 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]

Para visualizar cómo se distribuyen los sedimentos en el fondo del embalse, se representan las \textbf{curvas de isoconcentración}, que son las líneas donde la concentración $S(x,y)$ permanece constante. Estas curvas permiten identificar zonas de mayor o menor acumulación de sedimentos.

• • 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 y facilita la planificación de posibles medidas de limpieza o gestión del embalse.

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, 5, 15, 'LineWidth",1.2);
colorbar;
axis equal;
title('Curvas de isoconcentración de 5(x,y)');
xlabel('x (m)');
ylabel('y (m)');

4.2 -Flujo de entrada del río Lozoya

El estudio del flujo de entrada del río Lozoya permite analizar cómo el agua se distribuye y cómo transporta sedimentos o calor dentro del embalse. Para simplificar, se modela el flujo como un campo de velocidad tridimensional dependiente de la posición [math](x, y, z)[/math].

4.2.1 Modelo del campo de velocidad

El río Lozoya entra al embalse por el extremo norte. Modelamos el campo de velocidad del flujo de entrada como:

Campo 2D simplificado: [math] \vec{v}(x,y) = v_0 \, e^{\frac{-x^2 + y^2}{R^2}} (\cos\phi \, \hat{i} + \sin\phi \, \hat{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.

Se considera el dominio: [math](x,y) \in [-100,100]^2 \text{ m}, \quad \sqrt{x^2 + y^2} \le \rho_a[/math], donde [math]\rho_a[/math] es el radio de la presa a la altura [math]z = H_{\text{agua}}[/math].

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

donde:

• [math]v_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) se mantiene igual: [math](x,y) \in [-100,100]^2 \text{ m}, \quad \sqrt{x^2 + y^2} \le \rho_a[/math].

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

Para visualizar el flujo en MATLAB, se construyen planos horizontales (a nivel superficial) y planos verticales (secciones perpendiculares al eje del embalse). El campo vectorial se representa mediante vectores [math]\vec{v}(x,y,z)[/math]: La longitud del vector indica la magnitud de la velocidad. La dirección indica hacia dónde se mueve el agua en cada punto del plano. Esta representación permite identificar zonas de mayor velocidad, así como regiones de convergencia o dispersión.

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

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

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;

4.2.5 Cálculo del caudal de entrada integrando la velocidad

5 Clasificación de las presas - LAS PRESAS ESPAÑOLAS

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

• 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 - En España:

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

• https://www.iberdrola.com/conocenos/nuestra-actividad/energia-hidroelectrica/tipos-presa

• https://masqueingenieria.com/blog/tipos-de-presas-y-su-clasificacion/

• https://www.miteco.gob.es/es/agua/temas/seguridad-de-presas-y-embalses/principales-documentos-gestion-seguridad-presas-embalses/clasificacion.html

• https://www.iagua.es/data/infraestructuras/presas/almendra