Diferencia entre revisiones de «La presa de El Atazar (Grupo 1)»

r0 = 308;  % Radio base de la presa (aproximado)
b = 35;    % Curvatura del arco parabólico
H = 134;   % Altura de la presa
theta = linspace(3*pi/4, 5*pi/4, 100);  % Ángulo en radianes
z = linspace(0, H, 100);  % Altura en metros

% Creación del mallado en 2D para los ángulos y las alturas
[Z, Theta] = meshgrid(z, theta);

% Cálculo del radio para cada punto en la superficie
R = r0 + b * (1 - (Z.^2) / H^2);

% Conversión de coordenadas cilíndricas a cartesianas
X = R .* cos(Theta);
Y = R .* sin(Theta);

% Crear la figura
figure;
surf(X, Y, Z, 'FaceColor', 'b', 'EdgeColor', 'none');  % Gráfico 3D de la superficie
xlabel('Eje X (m)');
ylabel('Eje Y (m)');
zlabel('Eje Z (m)');
title('Superficie parametrizada de la presa');

% Configuración visual
axis equal;                      % Asegura la misma escala en todos los ejes
view(3);                          % Vista en 3D
grid on;                          % Muestra la cuadrícula

% Parámetros físicos
rho = 1000;  % Densidad del agua (kg/m^3)
g = 9.81;    % Aceleración de la gravedad (m/s^2)

% Cálculo de la presión en función de la altura
P = rho * g * (H - Z);  % Presión hidrostática en la superficie

% Visualización de la superficie con colores representando la presión
figure;
surf(X, Y, Z, P, 'EdgeColor', 'none');  % Superficie con mapa de colores basado en P
xlabel('Eje X (m)');
ylabel('Eje Y (m)');
zlabel('Eje Z (m)');
title('Presión sobre la superficie de la presa');

% Configuración de colores y visualización
colormap(jet);     % Paleta de colores para la presión
colorbar;          % Barra de colores para interpretar la presión
axis equal;        % Igualar las escalas de los ejes
view(3);           % Vista tridimensional
grid on;           % Mostrar la cuadrícula

% Parámetros de la presa
rho = 1000;        % Densidad del agua (kg/m^3)
g = 9.81;          % Aceleración gravitatoria (m/s^2)
H = 134;           % Altura de la presa (m)
r0 = 968/pi;       % Radio en la altura máxima (m)
b = 35;            % Factor de curvatura del arco parabólico (m)

% --- Gráfica 1: Toda la presa ---
% Coordenadas cilíndricas
theta = linspace(3*pi/4, 5*pi/4, 20);   % Ángulo θ
z = linspace(0, H, 20);                % Altura z
[Theta, Z] = meshgrid(theta, z);       % Mallado para θ y z
R = r0 + b * (1 - Z.^2 / H^2);         % Radio en función de z

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

% Derivadas parciales para las normales
dr_dz = -2 * b * Z / H^2;              % Derivada parcial de r respecto a z
n_r = 1 ./ sqrt(1 + dr_dz.^2);         % Componente radial del vector normal
n_x = n_r .* cos(Theta);               % Proyección de la normal en X
n_y = n_r .* sin(Theta);               % Proyección de la normal en Y

% Campo de presión
P = rho * g * (H - Z);                 % Presión en función de la profundidad

% Vectores de fuerza de presión
F_x = -P .* n_x;                       % Componente en X
F_y = -P .* n_y;                       % Componente en Y

% Representación gráfica de toda la presa
figure;
quiver3(X, Y, Z, F_x, F_y, zeros(size(F_x)), 1, 'b');
title('Fuerza de presión en la presa (3D)');
xlabel('X (m)');
ylabel('Y (m)');
zlabel('Z (m)');
grid on;
axis equal;

% --- Gráfica 2: Corte vertical ---
% Coordenadas cilíndricas para el corte vertical (θ = π)
theta_cut = pi;                        % Corte en θ = π
z_cut = linspace(0, H, 20);            % Altura z
R_cut = r0 + b * (1 - z_cut.^2 / H^2); % Radio en función de z

% Coordenadas cartesianas para el corte
X_cut = R_cut .* cos(theta_cut);       % Coordenadas X (constante para θ = π)
Z_cut = z_cut;                         % Coordenadas Z

% Derivadas parciales para calcular las normales
dr_dz_cut = -2 * b * Z_cut / H^2;      % Derivada parcial de r respecto a z
n_r_cut = 1 ./ sqrt(1 + dr_dz_cut.^2); % Componente radial del vector normal
n_x_cut = n_r_cut .* cos(theta_cut);   % Proyección de la normal en X
n_z_cut = -dr_dz_cut .* n_r_cut;       % Proyección de la normal en Z

% Campo de presión
P_cut = rho * g * (H - Z_cut);         % Presión en función de la profundidad

% Vectores de fuerza de presión (en el plano X-Z)
F_x_cut = -P_cut .* n_x_cut;           % Componente en X
F_z_cut = -P_cut .* n_z_cut;           % Componente en Z

% Representación gráfica del corte vertical
figure;
quiver(X_cut, Z_cut, F_x_cut, F_z_cut, 1, 'r');
title('Fuerza de presión en el corte vertical (X-Z)');
xlabel('X (m)');
ylabel('Z (m)');
grid on;
axis equal;

% Parámetros
theta = deg2rad(15); % Ángulo de salida en radianes
g = 9.81;           % Aceleración de la gravedad (m/s^2)
Hc = 25;            % Altura inicial (m)

% Velocidad inicial (previamente calculada)
v0 = sqrt(2 * g * Hc);

% Funciones de posición
t = linspace(0, 5, 1000); % Tiempo (ajustar si es necesario)
x = v0 * cos(theta) * t; 
y = Hc + v0 * sin(theta) * t - 0.5 * g * t.^2;

% Encontrar el tiempo cuando y(t) = 0 (impacto con el suelo)
t_ground = fzero(@(t) Hc + v0 * sin(theta) * t - 0.5 * g * t^2, [0, 10]);
x_ground = v0 * cos(theta) * t_ground; % Alcance horizontal

% Gráfica de la trayectoria
figure;
plot(x, y, 'b', 'LineWidth', 1.5,'DisplayName','gota');
hold on
plot(x_ground, 0, 'ro', 'MarkerSize', 8, 'DisplayName', 'Llegada');
xlabel('Distancia horizontal (m)');
ylabel('Altura (m)');
title('Trayectoria de la gota de agua');
grid on;
legend;
xlim([0, x_ground + 10]);
ylim([0, Hc + 5]);

% Parámetros iniciales
angulo = deg2rad(15);   % Ángulo de lanzamiento en radianes
gravedad = 9.81;        % Aceleración debido a la gravedad (m/s^2)
altura_inicial = 25;    % Altura inicial (m)
velocidad_inicial = sqrt(2 * gravedad * altura_inicial); % Velocidad inicial

% Tiempo y trayectorias
tiempo = linspace(0, 5, 500); % Intervalo de tiempo
pos_x = velocidad_inicial * cos(angulo) * tiempo; % Movimiento horizontal
pos_y = altura_inicial + velocidad_inicial * sin(angulo) * tiempo - 0.5 * gravedad * tiempo.^2; % Movimiento vertical

% Cálculo de las velocidades (derivadas de la posición)
vel_x = velocidad_inicial * cos(angulo) * ones(size(tiempo)); % Velocidad en x
vel_y = velocidad_inicial * sin(angulo) - gravedad * tiempo;  % Velocidad en y
velocidad_total = sqrt(vel_x.^2 + vel_y.^2);  % Magnitud de la velocidad

% Vectores unitarios tangente
tangente_x = vel_x ./ velocidad_total;
tangente_y = vel_y ./ velocidad_total;

% Cálculo de las aceleraciones (derivadas de la velocidad)
acel_x = zeros(size(tiempo)); % Aceleración en x (sin aceleración)
acel_y = -gravedad * ones(size(tiempo)); % Aceleración en y (gravedad)

aceleracion_total = sqrt(acel_x.^2 + acel_y.^2); % Magnitud de la aceleración

% Vectores unitarios normales
normal_x = (acel_x - tangente_x .* (acel_x .* tangente_x + acel_y .* tangente_y));
normal_y = (acel_y - tangente_y .* (acel_x .* tangente_x + acel_y .* tangente_y));

normal_x=normal_x./(sqrt(normal_x.^2+normal_y.^2));
normal_y=normal_y./(sqrt(normal_x.^2+normal_y.^2));
% Graficar trayectoria
figure;
plot(pos_x, pos_y, 'b', 'LineWidth', 1.5); hold on;

% Campos tangente y normal en puntos específicos a lo largo de la trayectoria
num_vectores = 20; % Número de puntos a graficar
puntos = round(linspace(1, length(tiempo), num_vectores));
quiver(pos_x(puntos), pos_y(puntos), tangente_x(puntos), tangente_y(puntos), 0.3, 'g', 'LineWidth', 1.5, 'DisplayName', 'Tangente');
quiver(pos_x(puntos), pos_y(puntos), normal_x(puntos), normal_y(puntos), 0.3, 'r', 'LineWidth', 1.5, 'DisplayName', 'Normal');

% Detalles adicionales de la gráfica
xlabel('Distancia horizontal (m)');
ylabel('Altura (m)');
title('Campos Tangente y Normal sobre la trayectoria');
legend('Trayectoria', 'Tangente', 'Normal');
grid on;
axis equal;

% Definir parámetros de la animación
t_anim = linspace(0, t_ground, 200);  % Intervalo de tiempo para la animación
x_anim = v0 * cos(theta) * t_anim;     % Posición en el eje X durante la animación
y_anim = Hc + v0 * sin(theta) * t_anim - 0.5 * g * t_anim.^2; % Posición en el eje Y

% Calcular las velocidades y aceleraciones en los puntos de la animación
vx_anim = v0 * cos(theta) * ones(size(t_anim)); % Velocidad en X
vy_anim = v0 * sin(theta) - g * t_anim;         % Velocidad en Y
speed_anim = sqrt(vx_anim.^2 + vy_anim.^2);    % Magnitud de la velocidad

% Aceleraciones constantes
ax_anim = zeros(size(t_anim));                % Aceleración en X
ay_anim = -g * ones(size(t_anim));            % Aceleración en Y

% Crear la figura para la animación
figure;
hold on;
plot(x, y, 'b', 'LineWidth', 1.5, 'DisplayName', 'Trayectoria');  % Trayectoria total
xlabel('Distancia horizontal (m)');
ylabel('Altura (m)');
title('Animación de una Gota de Agua');
grid on;
axis equal;
xlim([0, max(x)]);
ylim([0, Hc + 5]);

% Inicializar marcador para la gota
gota = plot(x_anim(1), y_anim(1), 'ro', 'MarkerSize', 10, 'DisplayName', 'Gota');

% Inicializar vectores de velocidad y aceleración
quiver_vel = quiver(x_anim(1), y_anim(1), vx_anim(1), vy_anim(1), 0.3, 'g', 'LineWidth', 1.5, 'MaxHeadSize', 1, 'DisplayName', 'Velocidad');
quiver_acc = quiver(x_anim(1), y_anim(1), ax_anim(1), ay_anim(1), 0.3, 'r', 'LineWidth', 1.5, 'MaxHeadSize', 1, 'DisplayName', 'Aceleración');

legend;

% Ejecutar la animación
for i = 1:length(t_anim)
    % Actualizar la posición de la gota
    set(gota, 'XData', x_anim(i), 'YData', y_anim(i));
    
    % Actualizar los vectores de velocidad y aceleración
    set(quiver_vel, 'XData', x_anim(i), 'YData', y_anim(i), ...
        'UData', vx_anim(i), 'VData', vy_anim(i));
    set(quiver_acc, 'XData', x_anim(i), 'YData', y_anim(i), ...
        'UData', ax_anim(i), 'VData', ay_anim(i));
    
    % Pausa para crear el efecto de animación
    pause(0.05);
end

% Definición de parámetros
area_compuerta = 1100;  % Área de la compuerta en metros cuadrados (m^2)
aceleracion_gravedad = 9.81;  % Aceleración debido a la gravedad en m/s^2
altura_salida = 25;     % Altura a la que sale el agua en metros (m)

% Cálculo de la velocidad de salida utilizando la ecuación de Torricelli
velocidad_salida = sqrt(2 * aceleracion_gravedad * altura_salida);

% Cálculo del caudal volumétrico (en m^3/s)
caudal_volumetrico = (area_compuerta * velocidad_salida);

% Mostrar los resultados
fprintf('Velocidad de salida (v0): %.2f m/s\n', velocidad_salida);
fprintf('Caudal volumétrico (Q): %.2f m^3/s\n', caudal_volumetrico);

% Datos
rho = 1000;        % Densidad del agua (kg/m^3)
g = 9.81;          % Gravedad (m/s^2)
H = 134;           % Altura de la presa (m)
r0 = 308.124;      % Radio en la altura máxima (m)

% Presa de doble curvatura (b = 35)
b = 35;                                 % Curvatura (m)
r = @(x) r0 + b * (1 - x.^2 / H^2);     % Radio
a = @(x) sqrt(r(x).^2 + (2 * b .* r(x) .* x / H^2).^2);  % Diferencial del área
P = @(x) rho * g * (H - x);             % Presión
I = @(x) P(x) .* a(x);                  % Integral de presión

A = pi/2 * integral(a, 0, H);           % Área de la superficie curva
F_total = pi/2 * integral(I, 0, H);     % Fuerza total de presión
f_porarea = F_total / A;                % Presión por unidad de área

% Presa de curvatura simple (b = 0)
b = 0;                                  % Curvatura (m)
r = @(x) r0 + b * (1 - x.^2 / H^2);     % Radio redefinido para la nueva curvatura
a = @(x) sqrt(r(x).^2 + (2 * b .* r(x) .* x / H^2).^2);  % Diferencial de área redefinido
P = @(x) rho * g * (H - x);             % Presión (sin cambios)

A0 = pi/2 * r0 * H;                     % Área de la superficie curva para curvatura simple
F_total0 = pi/2 * r0 * rho * g * H^2 / 2; % Fuerza total de presión para curvatura simple
f_porarea0 = F_total0 / A0;             % Presión por unidad de área para curvatura simple

fprintf("Fuerza total para la presa de doble curvatura: %e N\n", F_total);
fprintf("Presión por unidad de área para la presa de doble curvatura: %e N/m^2\n", f_porarea);
fprintf("Fuerza total para la presa de curvatura simple: %e N\n", F_total0);
fprintf("Presión por unidad de área para la presa de curvatura simple: %e N/m^2\n", f_porarea0);