Diferencia entre revisiones de «La presa de El Atazar»

clc
r0 = 200;                                      % Radio base de la presa (aproximado)
b = 35;                                        % Curvatura del arco parabólico
H = 134;                                       % Altura de la presa
theta_min = 3*pi/4;
theta_max = 5*pi/4;
theta = linspace(theta_min, theta_max, 100);   % Ángulo en radianes
z = linspace(0, H, 100);                       % Altura
[Theta, Z] = meshgrid(theta, z);               % Creación del mallado
R = r0 + b * (1 - (Z.^2) / H^2);               % Ecuación paramétrica para r
X = R .* cos(Theta);                           % Conversión a coordenadas cartesiana
Y = R .* sin(Theta);
figure;
surf(X, Y, Z, 'FaceColor', 'cyan', 'EdgeColor', 'none');
xlabel('Eje X (m)');
ylabel('Eje Y (m)');
zlabel('Eje Z (m)');
title('Superficie parametrizada de la presa (aguas arriba)');
axis equal;                                    
view(3);                                       %Visualizamos el gráfico en tres dimensiones
grid on;

% Parámetros de presión
rho = 1000;             % Densidad del agua (kg/m^3)
g = 9.81;               % Gravedad (m/s^2)
% Cálculo de la presión
P = rho * g * (H - Z);  % Presión en función de z
% Gráfica de la superficie con mapa de colores para la presión
figure;
surf(X, Y, Z, P, 'EdgeColor', 'none');
xlabel('X (m)');
ylabel('Y (m)');
zlabel('Z (m)');
title('Presión sobre la superficie de la presa');
colorbar;               % Barra de colores para indicar magnitudes de presión
colormap(jet);          % Esquema de colores
axis equal;
view(3);
grid on;

%APARTADO 3

% Derivadas numéricas para obtener vectores tangentes
[Rx_theta, Rx_z] = gradient(X);
[Ry_theta, Ry_z] = gradient(Y);
[Rz_theta, Rz_z] = gradient(Z);

% Producto cruzado de vectores tangentes para obtener 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_magnitude = sqrt(Nx.^2 + Ny.^2 + Nz.^2);

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

% 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, Z, Fx, Fy, Fz, 0.5, 'Color', 'yellow'); % Escala ajustable con el factor 2
hold on;
surf(X, Y, Z, P, 'EdgeColor', 'none', 'FaceAlpha',1); % Superficie semitransparente
colorbar;
xlabel('X (m)');
ylabel('Y (m)');
zlabel('Z (m)');
title('Campo vectorial de la fuerza de presión sobre la presa');
view(3);
grid on;

%APARTADO 3B


% Fijamos theta = 0 (corte longitudinal)
theta_cut = 0; % Corte en el plano x-z
X_cut = X(:, 51); % Extraemos el corte longitudinal (theta cerca de 0)
Z_cut = Z(:, 51); % Correspondiente en altura
P_cut = P(:, 51); % Presión en el plano

% Vectores normales en el corte
Nx_cut = Nx(:, 51); % Componente x del vector normal
Nz_cut = Nz(:, 51); % Componente z del vector normal
N_magnitude_cut = sqrt(Nx_cut.^2 + Nz_cut.^2);

% Normalización del vector normal
Nx_unit_cut = Nx_cut ./ N_magnitude_cut;
Nz_unit_cut = Nz_cut ./ N_magnitude_cut;

% Campo de fuerza de presión en el corte
Fx_cut = -P_cut .* Nx_unit_cut;
Fz_cut = -P_cut .* Nz_unit_cut;

% Gráfica del campo vectorial en el corte
figure;
quiver(X_cut, Z_cut, Fx_cut, Fz_cut, 'b', 'LineWidth', 1.5); % Vectores de fuerza
hold on;
plot(X_cut, Z_cut, 'k-', 'LineWidth', 2); % Contorno de la presa en el plano
xlabel('X (m)');
ylabel('Z (m)');
title('Campo de fuerza de presión en el plano longitudinal');
grid on;

%APARTADO 4

% 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);
hold on;
plot(x_ground, 0, 'ro', 'MarkerSize', 8, 'DisplayName', 'Impacto');
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]);