%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% %%
%% Doswell Vortex %%
%% %%
%% Analytical Solution %%
%% %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The motion is a pure rotation with
% vanishing radial velocity and
% tangential velocity:
%
% V = (sech r)^2 * tanh r
%
% The initial temperature depends only
% on the `North-south' direction:
%
% T0 = - tanh y
%
% The study is purely kinematic. The
% velocity is held constant in time.
%
% The evolving temperature field is
% obtained by integration the simple
% advection equation
%
% dT/dt = 0
%
% The analytical solution is easily found.
% Each fluid parcel is carried round a
% circle with constant angular velocity
% omega depending only on the radius.
% So the value of T(x,y) at time t is equal
% to its value at a point (x0,y0) displaced
% through an angle omega*t.
%
% [See Doswell, C. A., 1984: A kinematic
% analysis of frontogenesis associated
% with a nondivergent vortex. J. Atmos.
% Sci., 41, 1242-1248.]
%
% [See also EIgil Kaas, Tellus, 2007 or 8]
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clf
colormap('jet')
rmax = 4;
timemax = 30;
Ntime = 1*timemax;
dtime = timemax/Ntime;
Nx = 101; Ny = 101;
Nx = 201; Ny = 201;
xmax = 0.7*rmax; ymax = xmax;
dx = 2*xmax/(Nx-1); dy = dx;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Get fields on Cartesian grid.
for nx=1:Nx
X(nx) = -xmax+(nx-1)*dx;
end
for ny=1:Ny
Y(ny) = -ymax+(ny-1)*dy;
end
[XX,YY] = meshgrid(X,Y);
XX = XX'; YY = YY';
for ny=1:Ny
for nx=1:Nx
x = X(nx);
y = Y(ny);
r = sqrt(x^2+y^2);
theta = atan2(y,x);
Vrot = (sech(r))^2 * tanh(r);
TT(nx,ny) = -tanh(y);
UU(nx,ny) = -Vrot*sin(theta);
VV(nx,ny) = +Vrot*cos(theta);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Plot the initial temperature
contourf(XX,YY,TT)
axis('square'); colorbar
drawnow
pause
%% Plot the vortex flow
quiver(XX,YY,UU,VV)
axis('square')
hold on
SS = sqrt(UU.^2+VV.^2);
contour(XX,YY,SS);
hold off
drawnow
pause
surf(XX,YY,SS)
view(0,90); shading('interp')
axis('square'); colorbar
drawnow
pause
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Compute the solution of the
%% advection equation for each time.
for ntime=1:Ntime
time = ntime*dtime;
% Compute the departure point for each (x,y)
for ny=1:Ny
for nx=1:Nx
x = X(nx);
y = Y(ny);
r = sqrt(x^2+y^2);
theta = atan2(y,x);
omega = ((sech(r))^2*tanh(r))/(r+eps);
omt = omega*time;
r0 = r;
theta0 = theta - omt;
x0 = r0*cos(theta0);
y0 = r0*sin(theta0);
TT(nx,ny) = -tanh(y0*cos(omt)-x0*sin(omt));
end
end
% Plot the temperature at each time
surf(XX,YY,TT)
view(0,90)
% view(-90,60)
axis('square')
axis('off')
shading('interp')
colorbar
heading = ['Time = ' num2str(time) ];
title(heading)
drawnow
% pause
day = ntime*dtime;
if(day==1|day==3|day==5|day==10|day==20) pause; end
end % End of time loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%