# Disk Rotation Simulation inside a cylinder

### Arash Mohtat (view profile)

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');
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(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
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(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;