Diferencia entre revisiones de «El Vortice de Rankine (Grupo 11)»

%% Pregunta 1 - Vórtice de Rankine (versión simplificada)
clear; clc; close all;

% Parámetros EF4
R = 250; vR = 90;
Gamma = 2*pi*R*vR;
fprintf('Circulacion Gamma = %.3e m^2/s\n', Gamma);

% ===== FUNCIÓN VELOCIDAD TANGENCIAL =====
vtheta_fun = @(rho) (rho<=R).* (Gamma/(2*pi*R^2).*rho) + ...
                    (rho> R).* (Gamma./(2*pi*rho));

%% FIGURA 1: Perfil v_theta(rho)
rho = linspace(0,1000,1000);
vtheta = vtheta_fun(rho);

figure; hold on; grid on;
plot(rho(rho<=R), vtheta(rho<=R), 'r', 'LineWidth',2);
plot(rho(rho> R), vtheta(rho> R), 'b', 'LineWidth',2);
yl = ylim; plot([R R],[yl(1) yl(2)],'k--','LineWidth',1.3);

xlabel('\rho (m)'); ylabel('v_\theta(m/s)');
title('Perfil de velocidad tangencial del Tornado categoria EF4');
leg1 = legend('Núcleo EF4','Exterior EF4','Radio del ojo R','Location','northeast');
set(leg1,'FontSize',7);

%% FIGURA 2: Campo vectorial (vista superior) en x,y ∈ [-800,800]
L = 800; N = 25;
[x,y] = meshgrid(linspace(-L,L,N));
rho_xy = sqrt(x.^2+y.^2);  rho_safe = max(rho_xy,1e-6);

v_xy = vtheta_fun(rho_xy);
Ux = -v_xy .* (y ./ rho_safe);
Uy =  v_xy .* (x ./ rho_safe);

mask     = rho_xy <= L;      % dominio circular dentro del cuadrado [-L,L]^2
mask_in  = mask & (rho_xy <= R);
mask_out = mask & (rho_xy >  R);

figure; hold on; grid on; axis equal;
axis([-L L -L L]);   % <- aquí fuerzas el dominio x,y ∈ [-800,800]

quiver(x(mask_in),  y(mask_in),  Ux(mask_in),  Uy(mask_in), 0.7, 'r');
quiver(x(mask_out), y(mask_out), Ux(mask_out), Uy(mask_out), 0.7, 'b');

theta = linspace(0,2*pi,200);
plot(R*cos(theta), R*sin(theta),'k--','LineWidth',1.3);

xlabel('x (m)'); ylabel('y (m)');
title('Campo de velocidad tangencial del Tornado categoria EF4 (plano horizontal)');
leg2 = legend('Núcleo EF4','Exterior EF4','Radio del ojo R','Location','northeast');
set(leg2,'FontSize',7);

%% Rotacional y magnitud de la vorticidad del vórtice de Rankine
% Parámetros del tornado EF4
R    = 250;        % radio del núcleo (m)
vR   = 90;         % velocidad tangencial en rho = R (m/s)
Gamma = 2*pi*R*vR; % circulación (m^2/s)

% Dominio horizontal
L = 800;           % extensión máxima (m)
N = 200;           % puntos por dirección
[x,y] = meshgrid(linspace(-L, L, N));
rho   = sqrt(x.^2 + y.^2);

% Rotacional analítico (solo componente vertical):
%   (nabla x v)_z = Gamma/(pi*R^2) si rho <= R, y 0 si rho > R
wz = (rho <= R) .* (Gamma/(pi*R^2));   % campo escalar del rotacional

%% FIGURA Rotacional 

figure;
imagesc(linspace(-L,L,N), linspace(-L,L,N), wz);
set(gca,'YDir','normal');
axis equal tight;
colorbar;
xlabel('x (m)');
ylabel('y (m)');
title('Rotacional del campo de velocidad');

hold on;
theta = linspace(0,2*pi,300);
plot(R*cos(theta), R*sin(theta), 'k--', 'LineWidth', 1.5);

%% Campo de presión del vórtice (versión mínima con unidad en colorbar)
clear; clc; close all;

% Parámetros
R=250; vR=90; z0=2800;
P0=920e2; Pinf=1013e2;
rho=1.225; g=9.81;

Gamma=2*pi*R*vR;
v=@(r) (r<=R).*(Gamma/(2*pi*R^2).*r)+(r>R).*(Gamma./(2*pi*r));

% Mallado y presión
[x,z]=meshgrid(linspace(0,1000,200), linspace(0,z0,200));
p = (x<=R).*(P0+0.5*rho*v(x).^2 - rho*g.*z) + ...
    (x> R).*(Pinf-0.5*rho*v(x).^2 - rho*g.*z);

% Gráfica
figure;
[~,h]=contourf(x(1,:), z(:,1), p/100, 20);
set(h,'LineColor','k'); set(gca,'YDir','normal');

xlabel('\rho (m)');
ylabel('z (m)');
title('Campo de presión sobre el sección vertical');

c = colorbar;
c.Label.String = 'Presión (mbar)';   % ← UNIDAD AÑADIDA AQUÍ

axis tight;

<syntaxhighlight lang="matlab">
%% Campo de presión p(rho,z) en un vórtice de Rankine
clear; close all; clc;

%% Parámetros físicos
rho_f = 1.2;      % densidad del fluido [kg/m^3]
g     = 9.81;     % gravedad [m/s^2]

%% Parámetros del vórtice
Gamma = 400;      % circulación total [m^2/s]
a     = 100;      % radio del núcleo [m]
Omega = Gamma/(2*pi*a^2);

%% Dominio
rho_max = 1000;   % [m]
z0      = 1000;   % [m]

Nr = 500;         % resolución radial
Nz = 300;         % resolución vertical

rho = linspace(0,rho_max,Nr);
z   = linspace(0,z0,Nz);
[R,Z] = meshgrid(rho,z);

%% Perfil de velocidad v_theta(rho)
vtheta = zeros(size(rho));
inside  = rho <= a;
outside = rho >  a;

vtheta(inside)  = Omega .* rho(inside);
vtheta(outside) = Gamma ./ (2*pi*rho(outside));

%% Cálculo de presión radial p(rho)
p = zeros(size(rho));
p0 = 101325;  % presión de referencia (atm)

% región interna
p(inside) = p0 - rho_f*(Omega^2/2).*rho(inside).^2;

% valor en el borde para continuidad
p_a = p(find(inside,1,'last'));

% región externa
for i = find(outside)
    r = rho(i);
    p(i) = p_a - rho_f*(Gamma^2/(8*pi^2))*(1/a^2 - 1/r^2);
end

%% Extender p(rho) a p(rho,z): añadir hidrostática
P = p + rho_f*g*(z0 - Z);   % matriz completa NxM

%% --- Gráfica ---
figure('Position',[200 100 900 700]);
imagesc