No BSD License  

Highlights from
Disk Rotation Simulation inside a cylinder

Disk Rotation Simulation inside a cylinder

by

 

21 Feb 2006 (Updated )

The program animates and simulates the rotation od a disk inside a half cylinder

final_out=semicircle_rotation
function final_out=semicircle_rotation
% DISK ROTATION on semi-circular path
%        This program Simulates & Animates the motion of a disk on the
% inner surface of a lying cylinder.The disk CM oscillates inside a
% half-circle on the vertical plane.
%    The program generates a real-time animation(simultaneous calc/plot
% in each loop-turn).The solving-algorithm used is a step-by-step analysis
% in general,however in case of Steady Pure Rotation a direct DE-solving
% approach is adapted which makes use of Extended Symbolic Math Toolbox.
%    The Steady-state of P.R. occures from the beginning or after transient
% energy-wasting(damping) oscillations of the disk CM.And when occures,
% is recognized at the beginning by a simple sufficiency-criterion for mus
% or in the middle by a comparison-test of the successive max-altitudes of
% the CM.
%    The program is terminated when 'Done' button is pressed so the Results
% Diagrams are graphed & the final_out output struct is returned.
%
% Arash Mohtat   March-2005

%clear;
clear maplemex;
close all;

% Input dialog-box (+variable assignments)
Uinput=inputdlg({'Enter the radius of the cylindrical path:',...
                 'Enter the radius of the disk:',...
                 'Enter the initial postion-angle of disk CM:',...
                 'Enter the disk mass:',...
                 'Enter the coefficient of static friction:',...
                 'Enter the coefficient of dynamic friction:',...
                 'Step accuracy coefficient (near .01 recommended)',...
                 'Precision of subroutine for f (near .01 recommended)'},...
                 'input',1,{'5','1','0','1','0.5','0.3','0.01','0.025'});
g=10;
h=eval(Uinput{1});
R=eval(Uinput{2});
th0=eval(Uinput{3});
m=eval(Uinput{4});
mus=eval(Uinput{5});
mud=eval(Uinput{6});
step_coeff=eval(Uinput{7});
epsilon_f=eval(Uinput{8});
b=h-R;

% Graphic Settings (figure,axes,path-background,done-button)
fig=figure('name','Disk Rotation on a semi-circular path','numbertitle','off','color','white','renderer','OpenGL');
[xpath,ypath]=arc(h,0,0,pi,2*pi);
patch([-1.1*h,-1.1*h,xpath,1.1*h,1.1*h],[-1.2*h,0,ypath,0,-1.2*h],'k');
axis([-1.1*h 1.1*h -1.2*h 1.05*R]);axis off;
set(gca,'dataaspectratio',[1 1 1]);
hold on;
uicontrol('string','Done',...
    'callback','evalin(''caller'',''endloop1=true;endloop2=true;'')');

% Disk-related objects & chronograph
[xc,yc]=pol2cart(pi+th0,b);
[xp,yp]=pol2cart(pi+th0,h);
[Xdisk,Ydisk]=circle(R,xc,yc);
disk_h=patch(Xdisk,Ydisk,'y','edgecolor','r');
radius_h=plot([xc,xp],[yc,yp],'r');
text(0,R,'t=','HorizontalAlignment','right');
chrono=text(0,R,'0','HorizontalAlignment','left');

% Loop-initializing assignments
delta_t=step_coeff/sqrt(2*g/b);
t=-delta_t;
th=th0;delta_th=0;
phi=0;delta_phi=0;
v=0;w=0;
aT=0;alpha=0;
n=0;
endloop1=false;
endloop2=false;
fs=[0,0];fd=[0,0];
ymax=-R*sin(th0);
results(1,:)=[0,th0,0,0,0,...
    sign(pi/2-th0)*min([mus*m*g*sin(th0),abs(m*g*cos(th0)/3)]),NaN];

% Steady Pure Rotation from beginning?
Pure_Rotation=false;
if mus*sin(th)>=cos(th)/3
    Pure_Rotation=true;
    uiwait(msgbox('PURE ROTATION'));
    endloop1=true;
end

% Real-time animation/calculation loop1 (Step-by-step Algorithm)
while ishandle(fig)
    drawnow;
    n=n+1;
    t=t+delta_t;
    th=th+delta_th;
    phi=phi+delta_phi;
    %-----Disk plotting in position at t-----
    [xc,yc]=pol2cart(pi+th,b);
    xp=xc-R*cos(th-phi);
    yp=yc-R*sin(th-phi);
    [Xdisk,Ydisk]=circle(R,xc,yc);
    set(disk_h,'Xdata',Xdisk,'Ydata',Ydisk);
    set(radius_h,'Xdata',[xc,xp],'Ydata',[yc,yp]);
    set(chrono,'string',num2str(t));
    %----------------------------------------
    v=v+aT*delta_t;
    w=w+alpha*delta_t;
    N=m*v^2/b+m*g*sin(th);
    %-----Decision making subroutine for f-----
    if abs(v-R*w)<epsilon_f  
                       %epsilon_f can have a drastic effect on the results.
                       % try defaultinput with: 
                       %  mud=.03 & epsilon_f= 0.01/.001/.0001 ! 
        fs=sign(pi/2-th)*min([mus*N,abs(m*g*cos(th)/3)]);
        fd=NaN;
        f=fs;
    else
        fd=sign(v-R*w)*mud*N;
        fs=NaN;
        f=fd;
    end
    %-------------------------------------------
    aT=g*cos(th)-f/m;
    alpha=2*f/(m*R);
    delta_th=(0.5*aT*delta_t^2+v*delta_t)/b;
    delta_phi=w*delta_t+delta_th;
    results(n+1,:)=[t,th,phi,w,v,fs,fd];
    %-----Max altitude saving/comparing subroutine-----
    if results(end,5)*results(end-1,5)<=0 & n>3
        ymax=[ymax,-b*sin(th)];
        if abs(ymax(end)-ymax(end-1))<b/100
            uiwait(msgbox('Steady phase of PURE ROTATION'));
            Pure_Rotation=true;
            endloop1=true;
        end   
    end
    %---------------------------------------------------
    if endloop1
        break;
    end 
end

% Steady Pure Rotation solution by Maple:dsolve,numeric
if Pure_Rotation
    delta_t=3*delta_t;        % Because each loop-turn in loop2 takes
                              % almost 3times more time than loop1
                              % (note that in this algorithm increasing                             
    th0=th;                   %  delta_t does not affect accuracy!)
    phi0=phi;
    eval(['maple t0:=',num2str(t)]);
    eval(['maple h:=',Uinput{1}]);
    eval(['maple R:=',Uinput{2}]);
    maple b:=h-R;
    maple g:=10;
    eval(['maple th0:=',num2str(th0)]);
    maple de:=diff(theta(t),t,t)=2*g*cos(theta(t))/(3*b);
    maple sol:=dsolve({de,theta(t0)=th0,D(theta)(t0)=0},theta(t),numeric,output=listprocedure);
    maple th:=eval(theta(t),sol);
    maple thdot:=eval(diff(theta(t),t),sol);
    % Real-time animation/calculation loop2 (Direct DE-solving Algorithm) 
    while ishandle(fig)
        if endloop2
            break;
        end 
        drawnow;
        n=n+1;
        t=t+delta_t;
        th=mfun('th',t);
        phi=h/R*(th-th0)+phi0;
        [xc,yc]=pol2cart(pi+th,b);
        xp=xc-R*cos(th-phi);
        yp=yc-R*sin(th-phi);
        [Xdisk,Ydisk]=circle(R,xc,yc);
        set(disk_h,'Xdata',Xdisk,'Ydata',Ydisk);
        set(radius_h,'Xdata',[xc,xp],'Ydata',[yc,yp]);
        set(chrono,'string',num2str(t));
        v=b*mfun('thdot',t);
        w=(h/R-1)*v/b;
        fs=m*g*cos(th)/3;
        results(n+1,:)=[t,th,phi,w,v,fs,NaN];
    end   
end     


% Results Diagrams
figure('windowstyle','docked','numbertitle','off','name','theta vs time');plot(results(:,1),results(:,2));grid on;       
figure('windowstyle','docked','numbertitle','off','name','phi vs time');plot(results(:,1),results(:,3));grid on;
figure('windowstyle','docked','numbertitle','off','name','omega vs time');plot(results(:,1),results(:,4));grid on;
figure('windowstyle','docked','numbertitle','off','name','CM speed vs time');plot(results(:,1),results(:,5)); grid on;
figure('windowstyle','docked','numbertitle','off','name','touch-point speed vs time');plot(results(:,1),results(:,5)-R*results(:,4));grid on; 
figure('windowstyle','docked','numbertitle','off','name','friction vs time');hold on;
h1=plot(results(:,1),results(:,6),'b.');grid on;
h2=plot(results(:,1),results(:,7),'r.');grid on;
legend([h1,h2],'f_s','f_d'); 


final_out=struct('t',results(:,1),'theta',results(:,2),...
         'phi',results(:,3),'omega',results(:,4),'v',results(:,5),...
         'fs',results(:,6),'fd',results(:,7));
     
disp('theta defines the position of CM  &')
disp('phi is the angle that the marked radius makes with the radial direction')
disp('alpha is not brought in the output.However,you can easily obtain it by:')
disp('alpha=R*f/I=2*f/(m*R)')

%--------------------------------------------------------------------------
function [x,y]=arc(R,x0,y0,th1,th2)
th=linspace(th1,th2);
[x,y]=pol2cart(th,R*ones(1,100));
x=x+x0;
y=y+y0;
%--------------------------------------------------------------------------
function [x,y]=circle(R,x0,y0)
th=linspace(0,2*pi);
x=R*cos(th)+x0;
y=R*sin(th)+y0;

Contact us