clc
clear 
%% Spherical Velocity Field for Wire Drawing
% (B. Avitzur, Metal Forming: Processes and Analysis)

%% Parameters
di=0.5;                    % Initial wire diameter
df=0.30;                   % Final wire diameter
vi=1;                      % Input wire velocity
alfa=20;                   % Die half angle

n_delta_y=10;              % Number of grid in y-direction
n_delta_x=40;              % Number of grid in x-direction
l=1;                       % Wire legth
deltat=.001;               % Time step
N=500;                     % Number of steps

%%
alfa=alfa*pi/180;
vf=di^2*vi/df^2;           % Final velocity
rf=df/sin(alfa)/2; 
ri=di/sin(alfa)/2;
deltax=l/n_delta_x;
deltay=di/n_delta_y/2;
%% Mesh generation                                                         
x=[-ri-l:deltax:-ri];     % X position
y=[-di/2:deltay:di/2];    % Y position
[x,y]=meshgrid(x,y);
z=zeros(size(x));
axes1 = axes('CameraUpVector',[0 1 0]);
axis(axes1,[-ri-l 0 (-ri-l)/2 (+ri+l)/2]);
grid(axes1,'on');
hold(axes1,'all');
f1=surf(x,y,z);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for n=1:N
    for nn=1:numel(x)
        % Check position (Zone)
        if (x(nn)^2+y(nn)^2)<rf^2        
            vu=vf;
            vv=0;
        elseif (x(nn)^2+y(nn)^2)>ri^2
            vu=vi;
            vv=0;
        else
            teta=atan(abs(y(nn)/x(nn)));
            r=(x(nn)^2+y(nn)^2)^.5;
            vu=vf*rf^2*cos(teta)^2/r^2;
            vv=-y(nn)/abs(y(nn)+eps)*vf*rf^2*cos(teta)*sin(teta)/r^2;
        end;
        x(nn)=x(nn)+vu*deltat;
        y(nn)=y(nn)+vv*deltat;
    end;

   delete(f1);
   f1=surf(x,y,z);
   daspect([1 1 1]);
   pause(0.01)
end;
delete(f1);
f1=surf(x,y,z);
