Diferencia entre revisiones de «La presa de El Atazar. Grupo 16»

clear; clc;

% 1. Representar en MATLAB la superficie parametrizada de la presa en su
% cara de aguas arriba, en un unico color

% Asignacion de variables
H = 134;                              % Altura de la presa (metros)
r0 = 242;                             % Radio de la presa en la altura maxima (metros)
b = 35;                               % Curvatura del arco parabolico (metros)
theta = linspace(3*pi/4,5*pi/4,100);  % Angulo theta coordenadas cilindricas
z = linspace(0,H,100);                % Altura z coordenadas cilindricas
[Mthetha,Mz] = meshgrid(theta,z);     % Mallado
r = r0+b*(1-(Mz.^2/H^2));             % Ecuacion de la superficie de la presa

% Conversion de cilíndricas a cartesianas
X = r.*cos(Mthetha);
Y = r.*sin(Mthetha);

% Creacion del grafico
figure;
surf(X,Y,Mz,'FaceColor','#4DBEEE','EdgeColor','none');
title('Presa de El Atazar aguas arriba');
xlabel('X (m)');
ylabel('Y (m)');
zlabel('Z (m)');
axis equal;
grid on;

clear; clc;

% 2. Representar el campo escalar de presion como un mapa de colores sobre
% la superficie parametrizada de la presa (usa tambien colorbar)

% Asignacion de parametros
rho = 1000;     % Densidad del agua (Kg/(m^3))
g = 9.81;       % Gravedad (m/s^2)

% Relacion presion - altura
P = rho*g*(H-Mz);

% Creacion del gráfico
figure;
surf(X,Y,Mz,P,"EdgeColor","none");
title('Campo escalar de presion en la presa');
xlabel('X (m)');
ylabel('Y (m)');
zlabel('Z (m)');
axis equal;
grid on;
colormap("jet");
colorbar;

clear; clc;

% 3. Representación de campo de presiones 
% Derivadas y vectores normales
[Rx_theta, Rx_z] = gradient(X);
[Ry_theta, Ry_z] = gradient(Y);
[Rz_theta, Rz_z] = gradient(gridH);

% Producto cruzado para el vector normal
Nx = Ry_theta .* Rz_z - Rz_theta .* Ry_z;
Ny = Rz_theta .* Rx_z - Rx_theta .* Rz_z;
Nz = Rx_theta .* Ry_z - Ry_theta .* Rx_z;

% Magnitud del vector normal
N_magnitud = sqrt(Nx.^2 + Ny.^2 + Nz.^2);

% Vector normal unitario
Nx_unit = Nx ./ N_magnitud;
Ny_unit = Ny ./ N_magnitud;
Nz_unit = Nz ./ N_magnitud;

% Campo de fuerza de presión
Fx = -P .* Nx_unit;
Fy = -P .* Ny_unit;
Fz = -P .* Nz_unit;

% Representación del campo vectorial sobre la superficie
figure;
quiver3(X, Y, gridH, Fx, Fy, Fz, 0.5, 'Color', 'r');
hold on;
surf(X, Y, gridH, P, 'EdgeColor', 'none', 'FaceAlpha', 0.8);
colorbar;
xlabel('X (m)');
ylabel('Y (m)');
zlabel('Altura Z (m)');
title('Campo vectorial de la fuerza de presión sobre la presa');
view(-60, 20); % Cambia el ángulo de vista
axis equal;
grid on;

% Corte en planta
corte = 51; % Elegimos un corte en el plano longitudinal (aproximadamente theta = 0)
X_corte = X(:, corte);
Z_corte = gridH(:, corte);
P_corte = P(:, corte);

Nx_corte = Nx(:, corte);
Nz_corte = Nz(:, corte);
N_magnitud_corte = sqrt(Nx_corte.^2 + Nz_corte.^2);

Nx_unit_corte = Nx_corte ./ N_magnitud_corte;
Nz_unit_corte = Nz_corte ./ N_magnitud_corte;

Fx_corte = -P_corte .* Nx_unit_corte;
Fz_corte = -P_corte .* Nz_unit_corte;

% Representación del corte longitudinal
figure;
quiver(X_corte, Z_corte, Fx_corte, Fz_corte, 'r', 'LineWidth', 1.5);
hold on;
plot(X_corte, Z_corte, 'k-', 'LineWidth', 2);
xlabel('X (m)');
ylabel('Z (m)');
title('Campo de fuerza de presión en el corte longitudinal');
grid on;

clear; clc;

% Parámetros que vamos a utilizar
theta = 15 * pi / 180;  % Ángulo de salida (pasamos a radianes)
g = 9.81;               % Gravedad (m/s^2)
Hc = 25;                % Altura compuerta (m)

% Velocidad inicial, utilizamos la fórmula de Torricelli
v0 = sqrt(2 * g * Hc);

% Determinamos los factores posición
t = linspace(0, 5, 1000);                       % Tiempo
x = v0 * cos(theta) * t;                        % Posición horizontal (X)
y = Hc + v0 * sin(theta) * t - 0.5 * g * t.^2;  % Posición vertical (Y)

% Hallamos el tiempo para el cual la altura es 0, es decir, el momento en el que la gota impacta contra el suelo
t_impacto = fzero(@(t) Hc + v0 * sin(theta) * t - 0.5 * g * t^2, [0, 10]);
x_impacto = v0 * cos(theta) * t_impacto; % Distancia desde la compuerta hasta el punto de impacto con el suelo

% Representación de la gráfica de la trayectoria que describe la gota
figure;
plot(x, y, 'Color', '#4DBEEE', 'LineWidth', 1.5); % Ponemos el mismo color que hemos utilizado para la primera figura
hold on;
plot(x_impacto, 0, 'ro', 'MarkerSize', 8, 'DisplayName', 'Impacto con el suelo'); 
xlabel('Distancia horizontal de la gota (metros)');
ylabel('Altura de la gota (metros)');
title('Trayectoria de la gota de agua');
grid on;

% Insertamos una leyenda que describa lo que vemos
legend('Trayectoria gota', 'Impacto con el suelo'); 
xlim([0, x_impacto + 12]);
ylim([0, Hc + 7]);

clear; clc;

% Parámetros que vamos a utilizar
theta = 15 * pi / 180;    % Pasamos el ángulo inicial a radianes
g = 9.81;                 % Gravedad (m/s^2)
Hc = 25;                  % Altura compuerta (m)
v0 = sqrt(2 * g * Hc);    % Velocidad inicial (con Torricelli)

% Trayectoria de la curva
t = linspace(0, 5, 500);                        % Tiempo
x = v0 * cos(theta) * t;                        % Posición en el eje x
y = Hc + v0 * sin(theta) * t - 0.5 * g * t.^2;  % Posición en el eje y

% Velocidades (derivadas de la posición)
vx = v0 * cos(theta) * ones(size(t));           % Velocidad en x (dx/dt)
vy = v0 * sin(theta) - g * t;                   % Velocidad en y (dy/dt)
V = sqrt(vx.^2 + vy.^2);                        % Magnitud de la velocidad

% Vector tangente
Tx = vx ./ V;                                   % Componente tangencial en x
Ty = vy ./ V;                                   % Componente tangencial en y

% Aceleraciones (derivadas de las velocidades)
ax = zeros(size(t));                            % Aceleración en x (dvx/dt = 0)
ay = -g * ones(size(t));                        % Aceleración en y (dvy/dt = -g)
a_magnitude = sqrt(ax.^2 + ay.^2);              % Magnitud de la aceleración

% Vector normal
Nx = (ax - Tx .* (ax .* Tx + ay .* Ty)) ./ a_magnitude; % Componente normal en x
Ny = (ay - Ty .* (ax .* Tx + ay .* Ty)) ./ a_magnitude; % Componente normal en y

% Representación de la gráfica de los campos tangente y normal sobre la curva en varios instantes
figure;
plot(x, y, 'Color', '#4DBEEE', 'LineWidth', 1.5); hold on; 

% Campos tangente y normal en puntos específicos
puntos = 15; % Número de vectores a graficar
indices = round(linspace(1, length(t), puntos));  % Selección de índices

% Vectores tangentes en rojo
quiver(x(indices), y(indices), Tx(indices), Ty(indices), 0.3, 'r', 'LineWidth', 1.5, 'DisplayName', 'Tangente');

% Vectores normales en color azul oscuro
quiver(x(indices), y(indices), Nx(indices), Ny(indices), 0.3, 'Color', [0, 0, 0.5], 'LineWidth', 1.5, 'DisplayName', 'Normal');

xlabel('Distancia horizontal (m)');
ylabel('Altura (m)');
title('Campos Tangente y Normal sobre la trayectoria');
legend('Trayectoria', 'Tangente', 'Normal');
grid on;
axis equal;

% Ajustar los límites de la gráfica para asegurar que se vea la intersección con el eje x (y=0) y el eje y (x=0)
xlim([min(0, min(x) - 5), max(x) + 5]);
ylim([min(0, min(y) - 5), max(y) + 5]);

% Parámetros ajustados para la trayectoria
g = 9.81;                                            % Gravedad (m/s^2)
v0 = 30;                                             % Velocidad inicial (m/s)
theta = pi/3;                                        % Ángulo de lanzamiento (60 grados)
Hc = 25;                                             % Altura inicial (m)
t_max = 2 * v0 * sin(theta) / g + sqrt(2 * Hc / g);  % Tiempo total de vuelo
t = linspace(0, t_max, 200);                         % Intervalo de tiempo para la animación

% Posiciones en x e y durante la animación
x = v0 * cos(theta) .* t;               
y = Hc + v0 * sin(theta) .* t - 0.5 * g .* t.^2;

% Velocidades en x e y durante la animación
vx = v0 * cos(theta) * ones(size(t));    % Componente horizontal constante
vy = v0 * sin(theta) - g .* t;           % Componente vertical

% Aceleración (constante)
ax = 0;                                  % Componente horizontal de la aceleración
ay = -g;                                 % Componente vertical de la aceleración

% Ajustar límites para evitar valores negativos en y
valid_idx = y >= 0;                      % Solo considerar valores donde y >= 0
x = x(valid_idx);
y = y(valid_idx);
vx = vx(valid_idx);
vy = vy(valid_idx);

% Crear y configurar la figura
figure;
hold on;
plot(x, y, 'c', 'LineWidth', 1.5, 'DisplayName', 'Trayectoria gota'); % Trayectoria
xlabel('Distancia horizontal de la gota (metros)');
ylabel('Altura de la gota (metros)');
title('Trayectoria de la gota de agua con vectores de velocidad y aceleración');
grid on;
axis equal;
xlim([25, max(x)]);
ylim([0, Hc + 75]);

% Inicializar marcador de la gota y vectores
gota = plot(x(1), y(1), 'ro', 'MarkerSize', 10, 'DisplayName', 'Gota en movimiento');
vector_velocidad = quiver(x(1), y(1), vx(1), vy(1), 0.5, 'g', 'LineWidth', 1.5, ...
                          'MaxHeadSize', 2, 'DisplayName', 'Velocidad');
vector_aceleracion = quiver(x(1), y(1), ax, ay, 0.5, 'y', 'LineWidth', 1.5, ...
                            'MaxHeadSize', 2, 'DisplayName', 'Aceleración');

legend('Location', 'best');
hold off;

% Nombre del archivo GIF
gif_filename = 'animaciongotaG16.gif';

% Animar la trayectoria con vectores
for i = 1:length(x)
    % Actualizar la posición de la gota
    set(gota, 'XData', x(i), 'YData', y(i));
    
    % Actualizar el vector velocidad
    set(vector_velocidad, 'XData', x(i), 'YData', y(i), ...
                          'UData', vx(i), 'VData', vy(i));
    
    % Actualizar el vector aceleración
    set(vector_aceleracion, 'XData', x(i), 'YData', y(i), ...
                            'UData', ax, 'VData', ay);
    
    % Capturar el cuadro actual
    frame = getframe(gcf);
    im = frame2im(frame);
    [imind, cm] = rgb2ind(im, 256);

    % Escribir al archivo GIF
    if i == 1
        imwrite(imind, cm, gif_filename, 'gif', 'LoopCount', inf, 'DelayTime', 0.05);
    else
        imwrite(imind, cm, gif_filename, 'gif', 'WriteMode', 'append', 'DelayTime', 0.05);
    end
    
    % Pausa para crear el efecto de animación
    pause(0.05);
end

% Agregar marcador de impacto al final de la animación
hold on;
plot(x(end), y(end), 'ro', 'MarkerSize', 8, 'DisplayName', 'Impacto con el suelo');
legend('Location', 'best');
hold off;

clear; clc;

% Parámetros que vamos a utilizar
a = 1100;               % Área compuerta (m^2)
g = 9.81;               % Gravedad (m/s^2)
Hc = 25;                % Altura compuerta (m)
v0 = sqrt(2 * g * Hc);  % Velocidad inicial (con Torricelli)

% Caudal de la presa en metros cúbicos por segundo
Q = ( a * v0);

% Caudal de la presa en litros por segundo
Ql = Q ./ 1000;

% Caudales
disp(['Caudal en m³/s: ', num2str(Q), ' m³/s']);
disp(['Caudal en l/s: ', num2str(Ql), ' l/s']);

clear; clc;

rho = 1000; % Densidad del agua (kg/m^3)
g = 9.81;   % Gravedad (m/s^2)
H = 134;    % Altura de la presa (m)
r0 = 242;   % Radio de la presa en la altura máxima (metros)
b = 35;     % Factor de curvatura parabólica

% Resolución para la integración
num_theta = 100;     % Divisiones en theta
num_z = 100;         % Divisiones en z

% Rango de variables
theta = linspace(3*pi/4, 5*pi/4, num_theta);
z = linspace(0, H, num_z);

% Matriz de resultados
F_total = 0;
A_total = 0;

% Integración sobre la superficie
for i = 1:numel(theta)
    for j = 1:numel(z)
        % Coordenadas en la superficie
        r = r0 + b * (1 - (z(j)^2 / H^2));   % Radio a cada z
        h = H - z(j);                        % Profundidad
        
        % Elemento de área diferencial en coordenadas cilíndricas
        dS = r * (theta(2) - theta(1)) * (z(2) - z(1)); 
        
        % Presión a cada punto
        P = rho * g * h; 
        
        % Contribución a la fuerza total
        F_total = F_total + P * dS; 
        
        % Acumulación del área para calcular presión media (si necesaria)
        A_total = A_total + dS;
    end
end

% Fuerza total de presión
disp(['Fuerza total de presión: ', num2str(F_total), ' N']);

% Presión promedio por unidad de superficie
P_area = F_total / A_total;
disp(['Presión por unidad de superficie promedio: ', num2str(P_area), ' Pa']);