11/9/06

萬向接頭之計算

萬向接頭之轉換角度關係可利用MATLAB函數程式unijoint()計算,其呼叫語法如下:


function [th2, omega2, alpha2]=unijoint(betax,theta3,omega3)


相關的輸入參數:

betax:兩軸間之夾角。
theta3:迴轉之角度。
omega3:迴轉之角速度,rad/s。


程式內容:

程式11.1


function [th2, omega2, alpha2]=unijoint(betax,theta3,omega3)
%
% [th2, omega2, alpha2]=unijoint(betax,theta3,omega3)
% This program calculate the displacement, velocity and acceleration of a universal joint
% Inputs:
% betax: angle between two axes, degrees
% theta3: angle of rotations, degrees
% omega3: angular velocity, rad/sec
% Example:[t,v,a]=unijoint(45,0:10:360,1)
%
d2g=pi/180;
axang=betax*d2g;
th3=theta3*d2g;
cb=cos(axang);sb=sin(axang);
st=sin(th3);ct=cos(th3);
th2=atan(tan(th3)*cb);
omega2=cb./(1-(st.*sb).^2).*omega3;
alpha2=(omega3.*omega3.*cb.*sb.*sb).*sin(2*th3)./(1-(st.*sb).^2).^2;

行星齒輪組之計算

為計算不同組合之行星齒輪組(參考機動學圖10.22),利用planetary()程式可以直接解各模組之關係,其呼叫格式如下:

function [WF,WA,WL,RV]=planetary(wf,wa,wl,gearset,mode)

其中之參數:

wf, WF:第一輸入軸之角速度,ωF(rad/s)。
wl, WL:最後齒輪軸之角速度,ωL(rad/s)。
wa, WA:旋轉臂之角速度,ωA(rad/s)。
gearset:齒輪之齒數組合(由第二個齒輪算起),列矩陣。
mode: 行星齒輪組之型式代碼(1~12, 參考圖10.22)。


執行上項函式時,輸入wf, wa, wl三個參數中,可容許一個為未知,其未知參數可用[]取代,但不能用0,因為0將認為該軸靜止不動或固定。若不知各組之齒輪數目及搭配之尺寸,則可設gearset=[],程式會自動產生一個例子,可以代入執行。然後再設法修正。程式會自檢測所輸入之組合齒數是否適當,讀者可以自行檢驗。

程式內容:



function [WF,WA,WL,RV]=planetary(wf,wa,wl,gearset,mode)
%Program for section 10.6 to solve Levai's models
% Find angular velocities of sun,arm and planetary
% Variables:
% wf,WF:angular vel. for the 1st input shaft
% wl,WL:angular vel. for the last input shaft
% wa,WA:anbular vel. for the arm shaft.
% gearset:no of gear tooth in row matrix [n1 n2,..]
% mode:types of gear set (1~12, as in Fig. 9.22)
%Note:one of wf, wl & wa can be unknown and replaced by [].
% if leave gearset in [], the program will respond
% an example that matches the specific mode.
% Example:[WF,WA,WL,RV]=planetary(60,1,[],[18 30 28 20],7)
% Designed by D. S. Fon, Bime of NTU, dated March 8,2003
% Revised: March 9, 2006
% preset signs for each set
sinx=[-1 +1 -1 +1 -1 -1 +1 +1 -1 +1 -1 +1];
example={'[15 45 105]','[18 15 12 72]','[30 15 12 18 54]',...
'[30 15 20 20 105]','[80 15 12 30 137]','[40 15 20 75]'...
,'[80 20 40 60]','[100 30 50 120]','[80 15 12 18 12 107]'...
,'[60 30 15 30 15 150]','[40 20 12 20 40 52]',...
'[60 30 12 30 162]'};
WA=0;WF=0;WL=0;RV=0;
gearmin=12;
mode(mode>12)=12;
if isempty(gearset)
disp('Example: ')
example(mode)
return
end
gearset=[mode gearset];
[code,ge]=examset(mode,gearset,gearmin);
switch code
case 2
disp('The last gear size has been modified.')
ge(2:end)
case 3
disp('Error occurs in size matching')
return
end

% set up the no of gears required in the mode.
switch mode
case 1, rv=ge(4)/ge(2)*sinx(1);
case 2, rv=ge(5)/ge(2)*sinx(2);
case 3, rv=ge(4)/ge(2)*ge(6)/ge(5)*sinx(3);
case {4,5}, rv=ge(3)/ge(2)*ge(6)/ge(4)*sinx(mode);
case {6,7,8},rv= ge(3)/ge(2)*ge(5)/ge(4)*sinx(mode);
case {9,10,11},rv= ge(3)/ge(2)*ge(5)/ge(4)*ge(7)/ge(6)*sinx(mode);
otherwise
rv= ge(4)/ge(2)*ge(6)/ge(5)*sinx(mode);
end
if isempty(wa), wa=(rv*wl-wf)/(rv-1);
elseif isempty(wf), wf=rv*(wl-wa)+wa;
else wl=1/rv*(wf-wa)+wa;
end
WA=wa;WF=wf;WL=wl;RV=rv;

function [code,ge]=examset(mode,gearset,gemin)
% Check the proper gear set imported
% set up the no of gears required for a mode.
max_item_no=[4 5 6 6 6 5 5 5 7 7 7 6];
code=1;
%Forcing the minimum required teeth in a gear.
gearset(gearset<gemin)=gemin;
itemmax=max_item_no(mode);
item_no=length(gearset);
if itemmax>item_no,
gearset=[gearset ones(1,12)*12];
gearset=gearset(1:item_no);
end
ge=abs(gearset);
switch mode
case 1, gtmp=ge(2)+2*ge(3);
case 2, gtmp=ge(2)+2*ge(3)+2*ge(4);
case 3, gtmp=ge(2)+2*ge(3)+ge(4)-ge(5);
case 4, gtmp=ge(2)+ge(3)+ge(4)+2*ge(5);
case 5, gtmp=ge(2)-ge(3)+ge(4)+2*ge(5);
case 6, gtmp=ge(2)+ge(3)+ge(4);
case 7, gtmp=ge(2)+ge(3)-ge(4);
case 8, gtmp=ge(2)-ge(3)+ge(4);
case 9, gtmp=ge(2)-ge(3)+ge(4)+ge(5)+ge(6);
case 10, gtmp=ge(2)+ge(3)+ge(4)+ge(5)+ge(6);
case 11, gtmp=ge(2)+ge(3)+ge(4)+ge(5)-ge(6);
case 12, gtmp=ge(2)+2*ge(3)+ge(4)+ge(5);
end
if ge(itemmax)~=gtmp,
ge(itemmax)=gtmp;
code=2;
end
if mode==5|mode==8|mode==9,
if ge(2)<ge(3), code=3;end
end

執行例一:


使用其內部設定之範例,設定mode=3。

>> [WF,WA,WL,RV]=planetary(60,1,[],[],3)
Example:
ans =
'[30 15 12 18 54]'
WF =
0
WA =
0
WL =
0
RV =
0


再將上值代入,執行結果:

>> [WF,WA,WL,RV]=planetary(60,1,[],[30 15 12 18 54],3)
WF =
60
WA =
1
WL =
-48.1667
RV =
-1.2000


執行例二:


如圖10.14及範例10.6所述之型式如圖10.22之第一型,其中N2=15,N3=45,N4=105。若N4內齒輪固定,旋臂(連結N3齒)每轉一圈(即ωA=1),則執行下列指令並得結果如下:

>> [WF,WA,WL,RV]=planetary([],1,0,[15 45 105],1)

WF = 8
WA = 1
WL = 0
RV = -7


執行例三:


如圖10.15及範例10.8設其旋轉臂轉速ωA為1 rad/s,而N2=A20=20,N3=B60=60,N4=D40=40,N5=E120=120,N6=C140=140。由圖10.15可知,本齒輪可由圖10.22中之第一型與第六型綜合。故開始時可利用第一型先求出輸入軸A對應旋轉臂之轉速,然後依第六型求得齒輪E之轉速。即


>> [WF,WA,WL,RV]=planetary([],1,0,[20 60 140],1)
WF = 8
WA = 1
WL = 0
RV = -7


由其結果得到ω2=8 rad/s。其次依第六型求輸出軸E如下:

>> [WF,WA,FL,rv]=planetary(WF,1,[],[20 60 40 120],6)
WF = 8
WA = 1
WL = 0.22222
RV = -9
由其結果得到ωE=0.22222 rad/s。其中WF為利用上式所得的結果。

謝福吉(Selfridge)法

由於齒數需為整數,其最大及最小齒數為Nmax、Nmin,其關係如下:

Nmin=ceil(N2N4VRlow)≤N3N5≤fix(N2N4VRhigh)=Nmax

式中,函數CEIL()表示將取接近數值之高整數;而函數FIX()則取接近數值之低整數。

利用電腦重複測試之能力可以找出合適之齒數,使其值落於Nmax與Nmin之間。
其MATLAB程式函數two_stage_trains()如下:

function [result,counter]=two_stage_trains(rv,nmin,nmax,epx)

此程式處理兩組合齒輪組。輸入參數為:

rv:速度比。
nmin:最低齒數限。
nmax:最高齒數限。
epx:容許誤差。


輸出則為RESULT及counter。前者包括齒輪2及3之齒數,該組之轉速比;齒輪4及5之齒數,該組之轉速比;總轉速比及誤差。counter為計數實際無法達成之次數,程式中每連續試一次無法達到原設定之容許誤差時,誤差值會自動加倍,直到得到較滿意的答案為止。

程式內容

程式10.2

function [RESULT,counter]=two_stage_trains(rv,nmin,nmax,epx)
%
%prog 10.2
% Find the grear teeth for the set of coaxil gear trains
% rv: velocity ratio
% nmax,nmin:max.& min. number of teeths can be used for gears
% epx: error
% RESULT:a row matrix of eight elements, including
% N2 N3 RV1 N4 N5 RV2 RV error
% Example:
% two_stage_trains(pi,15,100,3.1416e-5)
epsx=epx;
counter=0;
[' N2 N3 RV1 N4 N5 RV2 RV 誤差']
ok=1;
while ok,
ss=1;
rv_h=rv+epx; rv_l=rv-epx;
nh3=fix(nmax^2/rv_h);
nh4=fix(nmax/sqrt(rv_h));
for pin1=nmin:nh4
nhh=min(nmax,fix(nh3/pin1));
for pin2=pin1:nhh
tm=fix(pin1*pin2*rv_h);
tn=ceil(pin1*pin2*rv_l);
if tm>=tn
nm=max(nmin,fix((tn+nmax-1)/nmax));
np=sqrt(tm);
for k=tn:tm
for gear1=nm:np
if (mod(k,gear1))==0
gear2=k/gear1;
error=(rv-k/(pin1*pin2));
if error<=epx ratio1=gear1/pin1; ratio2=gear2/pin2; ratio=ratio1*ratio2; result(ss,:)=[pin1,gear1,ratio1,pin2,gear2,ratio2,ratio,abs(error)]; ss=ss+1; end end end end end end end if ss==1 epx=epx*2; counter=counter+1; else ok=0; end end RESULT=sortrows(result,8);


執行例:


設計一同軸複式齒例,其轉速比為π,最高齒數為100齒,最低為15齒,容許誤差為e-5π。執行two_stage_trains()函式,結果如下:


>> a=two_stage_trains(pi,15,100,3.1416e-5)
ans =
N2 N3 RV1 N4 N5 RV2 RV 誤差
a =
25 51 2.04 50 77 1.54 3.1416 7.3464e-006
29 88 3.0345 85 88 1.0353 3.1416 1.0503e-005
22 62 2.8182 61 68 1.1148 3.1416 1.2922e-005
33 68 2.0606 61 93 1.5246 3.1416 1.2922e-005
43 77 1.7907 57 100 1.7544 3.1416 1.7786e-005
28 85 3.0357 86 89 1.0349 3.1416 1.8642e-005
43 85 1.9767 56 89 1.5893 3.1416 1.8642e-005
23 75 3.2609 82 79 0.96341 3.1416 2.3194e-005
41 75 1.8293 46 79 1.7174 3.1416 2.3194e-005
17 54 3.1765 91 90 0.98901 3.1416 2.8336e-005
17 60 3.5294 91 81 0.89011 3.1416 2.8336e-005
16 63 3.9375 94 75 0.79787 3.1416 2.9687e-005
32 63 1.9688 47 75 1.5957 3.1416 2.9687e-005

經過排序之結果,誤差最小之組合為:

N2:N3=25:51, N4:N5=50:77。

共軸齒列之計算

題意:


設計一組共軸齒列,其減速比VR=18。

做法:

  • 1. 雖然VR不一定為整數,但若為整數則會使解題更為簡化。
  • 2. 將VR=18開平方,得4.2426,其值小於每組合應低於10之限制,故可使用兩組合串列的方式,構成一齒列。
  • 3. 若能找到將一組正好等於4.2626之轉速比時即可使用雙組同比值的方式,其軸距即可相同,並符合共軸之條件。其步驟是先以小齒輪最低齒數12開始,相乘之,利用MATLAB指令可計算如下:


>>VR=18;  %整組之轉速比設為18。
>>mr=sqrt(VR); %以開方先取得單組轉速比估計值。得mr=4.2426。
>>T1=[12:40]’; %先設定驅動輪由12齒開始至40齒可能範圍,行矩陣。
>> T2=round(T1.*mr); %以估計值計算被動齒輪之齒數。
>>VR2=(T2./T1).^2; %計算實際之總轉速比。
>>Diff=(VR2-VR)/VR*100;%求其誤差百分比。
>>[T1 T2 VR2 Diff] %列出相關值

ans =
12.0000 51.0000 18.0625 0.3472
13.0000 55.0000 17.8994 -0.5588
14.0000 59.0000 17.7602 -1.3322
15.0000 64.0000 18.2044 1.1358
16.0000 68.0000 18.0625 0.3472
17.0000 72.0000 17.9377 -0.3460
18.0000 76.0000 17.8272 -0.9602
19.0000 81.0000 18.1745 0.9695
20.0000 85.0000 18.0625 0.3472
21.0000 89.0000 17.9615 -0.2142
22.0000 93.0000 17.8698 -0.7231
23.0000 98.0000 18.1550 0.8612
24.0000 102.0000 18.0625 0.3472
25.0000 106.0000 17.9776 -0.1244
26.0000 110.0000 17.8994 -0.5588
27.0000 115.0000 18.1413 0.7849
28.0000 119.0000 18.0625 0.3472
29.0000 123.0000 17.9893 -0.0595
30.0000 127.0000 17.9211 -0.4383
31.0000 132.0000 18.1311 0.7284
32.0000 136.0000 18.0625 0.3472
33.0000 140.0000 17.9982 -0.0102
34.0000 144.0000 17.9377 -0.3460
35.0000 148.0000 17.8808 -0.6621
36.0000 153.0000 18.0625 0.3472
37.0000 157.0000 18.0051 0.0284
38.0000 161.0000 17.9508 -0.2732
39.0000 165.0000 17.8994 -0.5588
40.0000 170.0000 18.0625 0.3472


其中第一行為驅動軸由12齒至40齒之變化;第二行為轉速比維持4.2626時之被動輪齒數(整數);第三行為實際計算之轉速比;最後一行為計算誤差百分比。就所得之值觀察,T1=12之誤差為0.3472%,差強人意。最小者為T1=33,其誤差為-0.0102%,不過被動輪為140齒,顯然嫌大些。

  • 4. 採用18之質因數,並將兩組合設成不同轉速比。在每組VR值不超過10之情況下,其可能集合有[9X2]、[6X3]等兩種。在齒輪組的考慮中應選齒數較不懸殊的組合,故取[6X3]。
  • 5. 設兩組速度比之組合為[6X3]。原則上小齒輪應不小於12齒,故所選的組合必須乘一個係數。至於選取何係數,兩個組合可以不同。

共軸齒列之計算

共軸齒列之計算可利用MATLAB程式函數coaxial_trains(),其如下:

function coaxial_trains(rv,n2,n3,teeth_limit,pitch_list)


輸入參數:

rv:速度比。
n2,n3:第一齒輪組合之齒數。
teeth_limit:所選齒輪之最低齒數。
pitch_list:可被選用之齒輪周節選單。



程式內容(10.1)



function coaxial_trains(rv,n2,n3,teeth_limit,pitch_list)
%
%prog 10.1
% pitch_list is a row vector for listing pitches of gears
% rv: velocity ratio
% teeth_limit:number of teeth limit for gears
% n2,n3: no. of teeth for the first pair gears
% Example:
% pitch_list=[1:.25:2,2.5,3,3.5,4:2:18,20:4:56]
% coaxial_trains(10,20,30,15,[1:.25:2,2.5,3,3.5,4:2:18,20:4:56])
nlist=length(pitch_list);
fact=n3/n2/rv;
[' N4 N5 Pd2 Pd4']
for i=1:nlist
pd2=pitch_list(i);
for j=1:nlist
pd4=pitch_list(j);
for k=teeth_limit:200
n4=k;
n5=n4/fact;
n4f=(n2+n3)/(1+1/fact)*(pd4/pd2);
if fix(n5)==n5 & abs(n4f-n4)<0.1
[n4,n5,pd2,pd4]
end
end
end
end


執行例:


設RV=10,N2=20,N3=30 及 最低齒數=15,求其可能之組合。

>> coaxial_trains(10,20,30,15,[1:.25:2,2.5,3,3.5,4:2:18,20:4:56])
ans = N4 N5 Pd2 Pd4
ans = 87 580 1.5 20
ans = 174 1160 1.5 40
ans = 15 100 1.75 4*
ans = 87 580 3 40
ans = 15 100 3.5 8*
ans = 24 160 12 44*
ans = 15 100 14 32*
ans = 18 120 16 44*

上述中徑節比較接近的如具星號的組合。程式之輸入項中,pitch_list為業界所使用之可能徑節。

漸開線函數

漸開函數可以查表或由計算程式中得知,而計算齒厚時,可以利用這一個漸開函數計算。用MATLAB之ainv()函式可以得到相對的答案,其呼叫格式如下:

function [epsilon, epsrad]=ainv(z)

其中,z為漸開函數值,得到之對應角度Φ為epsilon, epsrad,分別以度數及弧度表示。

function [epsilon, epsrad]=ainv(z)
% Find the inv funtion of involute angle.
% Using Secant's method.
% Input: z:(<2.2)
% outputs:epsrad, angle in radians
% epsilon: angle in degrees.
e1=0.63166;z1=0.1;
e2=0.97502;z2=0.5;
ok=0;
while ~ok
e0=(z-z1)/(z2-z1)*(e2-e1)+e1;
z0=tan(e0)-e0;
if abs(z0-z)<1e-6;
break;
else
if z0>z2
z2=z0;e2=e0;
else
z1=z0;e1=e0;
end
end
if abs(z2-z1)<1e-6; break; end
end
epsrad=e0;
epsilon=epsrad*180/pi;


執行例:



>> [ep,eprad]=ainv(0.1)
ep = 36.191
eprad= 0.63166
>> [ep,eprad]=ainv(1)
ep= 64.874
eprad = 1.1323

接觸比之計算

公式9.25是接觸比之計算過程,可以用MATLAB函式contact_ratio()得到結果。此函式之輸入值為徑節、兩齒輪之齒數及壓力角,輸出項包括公式9.26至公式9.32所計算之結果。其中之參數ag為列矩陣,依序儲存兩齒輪之接近角、退遠角與作用角。其呼叫格式如下:

function [c_ratio, c_length, ad, pc, pb, d2, d3, ag] = contact_ratio(pd, n2,n3, phi)


其中輸入參數:

Pd:徑節
n2, n3:兩齒輪之齒數
phi:壓力角


輸出參數:

cr_ratio:接觸比
cr_length:接觸長度
ad:齒冠
pc, pb:周節及基周節
d2, d3:兩齒輪節圓直徑。
ag:兩齒輪之接近角、遠退角及作用角


程式內容:

 
function [c_ratio,c_length,ad,pc,pb,d2,d3,ag]=contact_ratio(pd,n2,n3, phi)
%
%Find the contact ratios
% Inputs:
% Pd: Diametrial pitch;
% n2,n4:number of both gears;
% phi: pressure angle, degrees
% Outputs:
% c_ratio, c_length: contact ratio and length
% ad:addendium
%   pc,pb: circular and basic circular pitches
%   r2, r3: radii of pitch circles
%  ag: angles of action, in matrix of
% [alpha2 beta2 theta2 alpha3 beta3 theta3]
% Example: [c_r,c_l,ad,pc,pb,d2,d3,ag]
% =contact_ratio(6,24,48,20)
% Revised: March 9, 2006
d2g=pi/180;
pangle=phi*d2g;
cosx=cos(pangle);sinx=sin(pangle);
ad=1./pd;pc=pi./pd;
pb=pc.*cosx;
r2=n2./(2*pd);r3=n3./(2*pd);d2=2*r2;d3=2*r3;
rb2=r2.*cosx;rb3=r3.*cosx;
ax=sqrt((r3+ad).^2-(r3.*cosx).^2)-r3.*sinx;
xb=sqrt((r2+ad).^2-(r2.*cosx).^2)-r2.*sinx;
c_length=ax+xb;
c_ratio=c_length./pb;
ag1=[ax./rb2 xb./rb2 c_length./rb2]/d2g;
ag2=[ax./rb3 xb./rb3 c_length./rb3]/d2g;
ag=[ag1;ag2];


執行例一:


某對齒輪之周節為4,分別有24及48齒,其壓力角為20度,求其接觸比及相關資料。

>> [c_ratio, c_length,ad,pc,pb,r2,r3,ag]=contact_ratio(4,20,30,20)
c_ratio = 1.6052
c_length = 1.1847
ad = 0.2500
pc = 0.7854
pb = 0.7380
r2 = 5
r3 = 7.5000
ag =14.8816 14.0115 28.8932
9.9211 9.3410 19.2621


所得之資訊如下:

接觸比=1.6052
接觸長度=1.1847 吋
齒冠=0.25吋
周節pc=0.7854;基周節pb=0.7380
齒輪節圓直徑r2=5吋, r3=7.5吋
齒輪2之接近角=14.8816度;遠退角=14.0115度
作用角=28.8932度
齒輪3之接近角=9.921度;遠退角=9.3410度
作用角=19.2621度


執行例二:


某對齒輪之周節為5,分別有20及30齒,其壓力角為20及25度,求其接觸比及相關資料。

>> [c_ratio, c_length,ad,pc,pb,d2,d3,ag]=contact_ratio(5,20,30,[20 25]')
c_ratio =
1.6052
1.4419
c_length =
0.9477
0.8211
ad =
0.2000
pc =
0.6283
pb =
0.5904
0.5694
d2 =
4
d3 =
6
ag =
14.8816 14.0115 28.8932
13.2629 12.6921 25.9550
9.9211 9.3410 19.2621
8.8419 8.4614 17.3033


指令中,壓力角之輸入以行型式輸入(即[20 25]' ),其結果亦以行之型式對應出現。其他變數亦可如此作對照。這種比較不限1X2行矩陣,也不限於僅使用於一個輸入變數,但需同大小之矩陣。若不需要後面之輸出項,則僅提及前面所需之項目即可。讀者可自行印證不同變數矩陣之應用,若使用矩陣輸入,以採用行向量型式為宜。

外接錐之錐角比計算

兩外錐相接,其公式如下(參閱機動學公式9.9):

tanγ3= sinΣ / [ω32 + cosΣ]

若改用電腦程式計算,其程式如下(其中參數mode=1表示外錐;mode=-1表示內錐)::

function [gama1,gama2]=bevel_angle(omega,sigma,mode)
%
% Demo9.1:
% [gama1,gama2]=bevel_angle(omega,sigma,mode)
% mode:1 for outer set; -1 for innner set
% sigma: combined angle, degrees
% omega: angular velocity ratio=w1/w2, decimal
% gama1: cone angle for gear 1, degrees
% gama2: cone angle for gear 2, degrees
% Examples:[gama1,gama2]=bevel_angle([1 2 3 4 5],45,1)
d2g=pi/180;
sa=sigma*d2g;
gama1=atan(sin(sa)./(omega+mode*cos(sa)))/d2g;
gama2=sigma-mode*gama1;

-------------------------------------------

執行例一:


作一轉速比ω2 :ω3 = 5:3比例之圓錐,設合錐角Σ為60度,其MATLAB之程式如下:
ω12 =5/3;Σ=60;

>> [gama1,gama2]=bevel_angle(5/3,[60],1)
gama1 = 21.7868
gama2 = 38.2132

執行例二:


ω12=5/3;Σ=[60 50 40 30];


>> [gama1,gama2]=bevel_angle(5/3,[60 50 40 30],1)
gama1 = 21.7868 18.3506 14.8008 11.1676
gama2 = 38.2132 31.6494 25.1992 18.8324


若改為內錐外錐相接,其公式更動如下(參閱機動學公式9.10):

tanγ3= sinΣ / [ω32 - cosΣ]


執行例三:


ω1 /ω2 =7/3;Σ=[ 30];內接錐


>> [gama1,gama2]=bevel_angle(7/3,[30],-1)
gama1 = 18.8171
gama2 = 48.8171


執行例四:


ω1 /ω2 =[7/3 5/3 4/3]7/3;Σ=[ 30];內接錐


>>[gama1,gama2]=bevel_angle([7/3 5/3 4/3],[30],-1)
gama1 = 18.8171 31.9848 46.9357
gama2 = 48.8171 61.9848 76.9357

11/8/06

凸輪從動件升程之計算

範例8.1


某凸輪開始時先在0-120∘區間滯留,其高度為零。升程落於120-180∘,並於180-210∘滯留,其總升程為8公分。設刻度區間為10∘,凸輪則以等角速度旋轉,且其升程與返程均採用等加速度運動。求其變化曲線。

由機動學課本中之演繹,可以針對升程(Demo8.1)與返程( Demo8.2)計算結果如下:

% Demo8_1
% 計算升程對應點之資料:位移,速度及加速度
theta=120:10:180;
for i=1:length(theta)
[ss(i), vv(i), aa(i)]=parabol_cam(theta(i),120,60,1,8,0);
end;
%(theta:對應角度; ss:位移, vv:速度, aa:加速度)
[theta' ss' vv' aa']

執行結果:

>> demo8_1
ans =
120.0000 0 0 29.1805
130.0000 0.4444 5.0930 29.1805
140.0000 1.7778 10.1859 29.1805
150.0000 4.0000 15.2789 -29.1805
160.0000 6.2222 10.1859 -29.1805
170.0000 7.5556 5.0930 -29.1805
180.0000 8.0000 -0.0000 -29.1805


返程之計算程式如下:

% Demo8_2
% 計算返程對應點之資料:位移,速度及加速度
theta=210:10:360;
for i=1:length(theta)
[ss(i), vv(i), aa(i)]=parabol_cam(theta(i),210,150,-1,8,0);
end;
%(theta:對應角度; ss:位移, vv:速度, aa:加速度)
[theta' ss' vv' aa']


執行結果:

>> demo8_2
ans =
210.0000 8.0000 0 -4.6689
220.0000 7.9289 -0.8149 -4.6689
230.0000 7.7156 -1.6297 -4.6689
240.0000 7.3600 -2.4446 -4.6689
250.0000 6.8622 -3.2595 -4.6689
260.0000 6.2222 -4.0744 -4.6689
270.0000 5.4400 -4.8892 -4.6689
280.0000 4.5156 -5.7041 -4.6689
290.0000 3.4844 -5.7041 4.6689
300.0000 2.5600 -4.8892 4.6689
310.0000 1.7778 -4.0744 4.6689
320.0000 1.1378 -3.2595 4.6689
330.0000 0.6400 -2.4446 4.6689
340.0000 0.2844 -1.6297 4.6689
350.0000 0.0711 -0.8149 4.6689
360.0000 0.0000 -0.0000 4.6689

拋物線運動

凸輪之從動件作上下運動,其運動型式可採拋物線運動,利用MATLAB程式可以計算抛物線之運動參數,其程式內容如下,呼叫之名稱與參數分別為:

function [y, yy, yyy]=parabol_cam(phi, phi_in, beta_range, direct, travel,rpm)


輸出入參數:

y=從動件之位移。
yy=從動件之速度。
yyy=從動件之加速度。
phi=凸輪角度
beta=運動區間度數。
phi_in=起始角度
direct=運動型式,升程為+1 ;返程為 -1 。
travel=衝程(設為1時為單位量)。
rpm=凸輪迴角速度, rpm(設為0時不考慮RPM之影響)


8-4.3.2 MATLAB計算程式



function [y, yy, yyy]=parabol_cam(phi,phi_in,beta_range,direct,travel,rpm)
% code = 2
% [y, yy, yyy]=parabol_cam(phi,phi_in,beta_range,direct,travel,rpm)
% phi=cam angle, degrees
% phi_in=starting cam angle, degrees
% beta_range=motion range, degrees
% direct=motion type; +1 for upward, -1 for downward
% travel=the follower travels(=1 for a unit of travel)
% rpm=cam rotation speed, rpm(=0 for not effective)
% Example: [y, yy, yyy]=parabol_cam(10,180,120,1,1,0)
d2r=pi/180;
th=phi*d2r; thinit=phi_in*d2r; beta=beta_range*d2r;
speed=rpm*2*pi/60;if rpm==0, speed=1;end;
theta=th-thinit;
thmed=thinit+beta/2;
thx=theta/beta;
if direct==1,
if th<thmed
y=2*thx^2;
yyy=4/beta/beta;
yy=yyy*th;
else
y=1-2*(1-thx)^2;
yy=4/beta*(1-thx);
yyy=-4/beta/beta;
end
else
if th<thmed
y=1-2*thx^2;
yyy=-4/beta/beta;
yy=yyy*th;
else
y=2*(1-thx)^2;
yy=-4/beta*(1-thx);
yyy=4/beta/beta;
end
end
y=y*travel;
yy=yy*speed*travel;
yyy=yyy*speed*speed*travel;


執行例:

1. 升程起始點120∘,β=60∘,h=0.8cm,角度為150∘時之位移、速度及加速度:

>> [y, yy, yyy]=parabol_cam(150,120,60,1,0.8,0)
y = 0.4 (cm)
yy = 1.5279 (cm/s)
yyy = -2.9181 (cm/s2)


2. 返程起始點210∘,β=150∘,h=0.8,角度為250∘時之位移、速度及加速度:

>> [y, yy, yyy]=parabol_cam(250,210,150,-1,0.8,0)
y = 0.68622(cm)
yy = -0.3259(cm/s)
yyy = -0.46689(cm/s2)

五連桿齒輪解析

機動學中之公式7.30形成一矩陣,經簡化為:AW=C之型式。再以MATLAB指令W=A\C解出桿3、4、5之角速度。程式內容參見f5link()。其中之呼叫參數說明如下:

function [omega,theta,incta]=f5link(r1,R2,R4,r5,theta2,omega2)


輸入:r1、r5為第一(地桿)及第五桿之長度。
R2、R4為成對之齒輪半徑。
   omega2為桿2之驅動角速度,rad/s。
   theta2各桿2之對應角度。
輸出:omega為桿2、3、4、5之角速度。
   theta為桿2、3、4、5之水平角度。
   Incta為桿3與桿4成一直線時,桿3、4、5之初始角度。

利用下面之程式可以直接進行求解:


function [omega,theta,incta]=f5link(r1,R2,R4,r5,tha2,w2)
% Program 7.2
% Solve the angular velocities of a 5-link system
% Input:
% r1,r5: the length of ground and the 5th links.
% R2,R4: radius of the gear pair%
% tha2: angular position of link 2, degrees
% w2: angular velocity of link 2, rad/s
% Output:
% omega:the ang. velocitis of links, rad/s
% theta:the angles of links, degrees
% incta:initial positions of links, degrees
% omega=[omega3 omega4 omega5], in rad/s
% Example [omega,theta,incta]=f5link(8,3,3.5,6,60,10)
d2r=pi/180;th2=tha2*d2r;
r3=R2+R4;
if r3+r5<=r5,
theta='Not available';
incta=0;
return
end
r34=r3+R4;
cos5=(r34*r34-r5*r5-r1*r1)/(2*r1*r5);
cos3=(r1+r5*cos5)/r34;
th50=acos(cos5);
th30=acos(cos3);th40=th30;
incta=[th30 th40 th50]/d2r;
%============================================
% To solve the angular position of each links
rr=R4/R2;
func_1=1; func_2=1;
% Force next WHILE statement to be true
th5=th50;th4=th40;
while abs(func_1)>0.01 | abs(func_2)>0.01,
% Evaluate loop equations
th3=th30+(th2+rr*(th4-th40))/(1+rr);
func_1=r3*cos(th3)+R4*cos(th4)-r5*cos(th5)-r1;
func_2=r3*sin(th3)+R4*sin(th4)-r5*sin(th5);
F=[-func_1; -func_2];
% Evaluate partial derivatives
A=[-r3*sin(th3)*rr-R4*sin(th4) r5*sin(th5);...
r3*cos(th3)*rr+R4*cos(th4) -r5*cos(th5)];
% Solve 2 equations in 2 unknowns with Cramer's rule
delta=A\F;
% Make new guess for both theta3 and theta4
th4=th4+delta(1);
th5=th5+delta(2);
end;
theta=[0 th2 th3 th4 th5]/d2r;
%============================================
% Calculate the angular velocities of links
% Set up the 3x3 maxtrix A
A=[-r3*sin(th3) -R4*sin(th4) r5*sin(th5);...
r3*cos(th3) R4*cos(th4) -r5*cos(th5);...
1+rr -rr 0];
C=[0 0 w2]';
W=A\C; %Matrix division, or solution of AW=C
omega=[w2,W'];



執行例一:


設 r1=8 cm、R2=3 cm、R4=3.5 cm、r5= 6cm,θ2= 60度,ω2=10 rad/s

>> [w,theta,incta]=f5link(8,3,3.5,6,60,10)
w =   10.0000 2.7194 -3.5211 1.8456
theta = 0 60.0000 52.1875 13.8419 95.8522
incta = 36.8699 36.8699 90.0000


故ω3=2.7194 rad/s,ω4= -3.5211 rad/s, ω5= 1.8456 rad/s

執行例二:


設r1=50 cm、R2=12 cm、R4=14 cm、r5= 50cm,θ2= 30度,ω2=10 rad/s

>> [omega ,theta,incta]=f5link(50,12,14,50,30,10)
omega = 10.0000 2.4549 -4.0123 0.4351
theta = 0 30.0000 73.5832 54.0023 133.5048
incta =   66.4218 66.4218 132.8436


故ω3= 2.4549 rad/s,ω4=-4.0123 rad/s, ω5= 0.4351 rad/s

本程式中,實際上綜合f5link_init()函式之內容,並進行作位置解,將公式7.18及7.19利用疊代法求得各桿之對應角度。其後再利公式7.30之矩陣除法求得各桿之角速度。在求解過程中,桿2之角度及角速度均必須為已知值。

齒輪組構成之五連桿分析

下面為利用MATLAB撰寫之程式函式f5link_int(),以計算一組齒輪構成之五連桿(參閱機動學範例7.1),其齒輪2與4之直徑分別為6、7cm。其他桿之相關尺寸為r1=8cm,r5=6cm。齒輪2為驅動輪。在初始位置時ABC正好成一直線,設桿2角度之初始值為零。當桿2旋轉60度角時,其他相關之θ3、θ4、θ5值之變化可由f5link_int()程式求得,其呼叫型式為:

function [th3,th4,th5]=f5link_int(r)


其輸入項為r。其中包括五桿之長度,其中第二桿為驅動齒輪之半徑r2。

function [th3,th4,th5]=f5link_int(r1,R2,R4,r5)
% Program 7.1 solving eqs 7.11-7.17
% The 5 link system consists of 2 gears with one
% holding link(r5), except the ground(r1)
% Find the intial positions of the 5-link system
% Input:
% r1,r5: the length of ground and the 5th links.
% R2,R4: radius of the gear pair
%
% Output:
% th3,th4,th5: angles of link 3,4 & 5, degrees
% Example: [th3,th4,th5]=f5link_int(8, 3, 3.5, 6)
d2r=pi/180;
r3=R2+R4;
r34=r3+R4;
cos5=(r34*r34-r5*r5-r1*r1)/(2*r1*r5);
cos3=(r1+r5*cos5)/r34;
th50=acos(cos5);
th30=acos(cos3);
th3=th30/d2r;th4=th3;th5=th50/d2r;

執行例:

>> [th3,th4,th5]=f5link_int(8, 3, 3.5, 6)
th3 = 36.8699
th4 = 36.8699
th5 = 90

11/7/06

滑塊位置分析

利用MATLAB程式slider_ solve (),可以計算任意曲柄角度之對應連結桿角度及水平距離,公式4.24與4.25,如程式4.8。其呼叫函數如下:

function [d,theta3]=slider_solve(theta2,R,L,e,mode)


其參數之定義詳如程式內容。

程式4.8
function [d,theta3]=slider_solve(theta2,R,L,e,mode)
% P4.8 [d,theta3]=slider_solve(theta2,R,L,e)
%
% Find the block coordinate and angle of connecting rod in a slider-crank system
% Input:R,L lengths of crank & connecting rod; e the offset of block
% theta2: angles of crank, cane be a matrix
% Output: d:distance of the block from the origin.
% theta3:angle of the connecting rod,
% Examples:
% [d,theta3]=slider_solve(45,5,10,3)
%
% Designed by D.S. Fon, BIME, NTU, Date:Fegruary 9,2003.
if nargin<5, mode=0;end
d2g=pi/180;
theta=theta2*d2g;
cc=(e-R.*sin(theta))./L;
if mode>=0,
theta3=asin(cc);
else
theta3=asin(-cc)+pi;
end
d=L.*cos(theta3)+R.*cos(theta);
theta3=theta3/d2g;



範例4.9


試比較一組曲桿滑塊機構R=5cm,L=10cm,e=6 cm,試求曲柄角度為30、45、60、90度時之連結桿角及滑塊之水平距離。

[解]:

>> [d,theta3]=slider_solve([30 45 60 90],5,10,6,1)

d = 13.6976 13.2271 12.3596 9.9499
theta3 = 20.4873 14.2673 9.6127 5.7392

偏置滑塊之運動

偏置滑塊之運動受制於其幾何關係,利用MATLAB程式slider_limit(),可以計算衝程及界限距離,如程式4.7。其呼叫函數如下:

function [s,theta21,theta22]=slider_limit(R,L,e)

式中,R、L之定義同前,e為偏置量。s 為衝程,theta21與theta22則分別為界限角。slider_limit()函式並未檢查上述之R、L、e間之幾何關係,但若所得之結果為虛根時,表示該關係並不存在。


程式4.7
function [s,theta21,theta22]=slider_limit(R,L,e)
% P4.7 [s,th1,th2]=slider_limit(r)
%
% Find the stroke, extreme angles of link 2 in a slider mechanism system
% Input:R,L lengths of crank & connecting rod; e the offset of block
% Output: s:stroke, theta21:angle for the right position,
% theta22:angle for the left position
% Examples:
% [s,th1,th2]=slider_limit(5,12,3)
%
% Designed by D.S. Fon, BIME, NTU, Date:October 17,2002.
d2g=pi/180;
th1=asin(e./(R+L));
th2=asin(e./abs(R-L));
s=(R+L).*cos(th1)-abs(R-L).*cos(th2);
theta21=th1/d2g;
theta22=th2/d2g;


範例4.9


試比較一組曲桿滑塊機構R=5cm,L=10cm,而偏置量e分別為2 cm、3 cm、 5 cm及7 cm之衝程及界限角度。

[解]:

>> [s,th1,th2]=slider_limit(5,10, [2 3 5 7])

s = 10.2835 10.6969 14.1421 13.2665 + 4.8990i
th1 = 7.6623 11.5370 19.4712 27.8181
th2 = 23.5782 36.8699 90.0000 90.0000 +49.6763i


當偏置量為7cm,衝程s有虛根存在,表示該點無法成立。

滑塊運動

利用MATLAB程式slider_crank(),可以計算機動學中之公式4.22與4.23。其呼叫函數如下:

function [x,V,A]=slider_crank(theta2,R,L,omega)



程式4.6
function [x,V,A]=slider_crank(theta2,R,L,omega)
%
%P4.6 function [x,V,A]=Slider_crank(theta,R,L, omega)
% Inputs: theta2: Crank angles in degrees, single value or row matrix
% R: Crank length L: length of connecting rod
% omega:constant angular velocity(=1 for default)
% Outputs: x,V,A:displacements, velocities, accelerations of block.
% Example: [x,V,A]=Slider_crank(40,5,10,10)
%
if nargin<4, omega=1;end;
d2g=pi/180;
theta=theta2*d2g;
LL=2*L;
sinx=sin(theta);
cosx=cos(theta);
x=R.*(1-cosx)+(R.*sinx).^2./LL;
V=R*omega.*(sinx+R./LL.*sin(2*theta));
A=R*omega*omega.*(cosx+R./LL.*cos(2*theta));

範例4.8


有一組曲柄滑塊機構在水平線上運轉,曲桿長為5cm,連桿為10cm,其曲桿以角速度10rad/s迴轉,若將全周以60度分隔,試求其對應位移、速度及角速度為何?

[解]:

設R=5,L=10,omega=10,呼叫函數,其結果如下:

>> [x,V,A]=Slider_crank([60:60:360],10,10)

x = 8.7500 18.7500 20.0000 18.7500 8.7500 0.0000
V =12.9904 4.3301 0 -4.3301 -12.9904 -0.0000
A =2.5000 -7.5000 -5.0000 -7.5000 2.5000 15.0000