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 cilindricas 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;

%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 terrestre

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

%Creacion del grafico
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;

% 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 longitudinal
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;