第三章 四连杆之MATLAB程式讲解.docx
《第三章 四连杆之MATLAB程式讲解.docx》由会员分享,可在线阅读,更多相关《第三章 四连杆之MATLAB程式讲解.docx(32页珍藏版)》请在冰豆网上搜索。
第三章四连杆之MATLAB程式讲解
第三章四連桿之MATLAB程式
第三章中之四連桿分析可以參考相關資料。
本節則針對四連桿之動作程式加以說明。
目前所設計之程式有f4bar.m、drawlinks.m、fb_angle_limits.m、drawlimits.m等四個程式,茲分別說明如下:
圖一、四連桿之關係位置及各桿名稱
一、f4bar函數:
f4bar函數之呼叫格式如下:
function[values,form]=f4bar(r,theta1,theta2,td2,tdd2,sigma,driver)
輸入變數:
.r(1:
4)=各桿之長度,r
(1)為固定桿,其餘分別為曲桿、結合桿及被動桿。
.theta1=第一桿之水平角,或為四連桿之架構角,以角度表示。
.theta2=驅動桿之水平夾角,以角度表示。
一般為曲桿角,但若為結合桿驅動,則為結合桿之水平夾角。
.td2=驅動桿(第二桿或第三桿)之角速度(rad/sec)。
.tdd2=驅動桿(第二桿或第三桿)之角加速度(rad/sec^2)。
.sigma=+1or-1.組合模式,負值表示閉合型,正值為分支型,但有時需視實際情況而定。
.driver=0(驅動桿為第二桿);1(驅動桿為第三桿)
輸出變數:
.form=組合狀態,0:
表示無法組合;1:
可以正確組合
.values=輸出矩陣,其大小為4X7,各行之資料分配如下:
1
2(deg)
3(rad/s)
4(rad/s2)
5
6
7
I
桿1位置
θ1
ω1
α1
VQ
|VQ|
∠VQ
II
桿2位置
θ2
ω2
α2
VP
|VP|
∠VP
III
桿3位置
θ3
ω3
α3
AQ
|AQ|
∠AQ
IV
桿4位置
θ4
ω4
α4
AP
|AP|
∠AP
其中第一行之連桿位置向量,屬於單桿的位置向量。
第二行為各桿之水平夾角,第三及第四行為各桿之角度速度及角加速度。
第五至七行則為P點與Q點之速度與加速度量,第五行為向量,第六行為絕對量,第七行為夾角。
值得一提的是第一行、三行、四行及五行之向量表示法屬於複數之型式。
故若要得到其絕對值僅需在MATLAB指令檔中,以abs()這一個函數指令即可求得,而以函數angle()則可求得其夾角,雖然第二行與第七行之輸出亦有相對應之夾角。
例一:
為第二桿為驅動桿
[val,form]=f4bar([3242],0,60,10,0,-1,0)
val=Columns1through3
300
1+1.7321i6010
3.8682-1.0182i-14.7465.4078
1.8682+0.71389i20.91316.549
Columns4through6
01+1.7321i2
01.8682+0.71389i2
-127.58173.21-100i200
-236.27364.19-953.09i1020.3
Column7
60
20.913
-30
-69.087
form=1(表示可以組合)
本例中,有框線者表示其為輸入值。
但第一行則已經轉換為複數型式。
未來複數型式要轉為x-y座表時,只要使用函數real()及imag()兩指令,即可進行轉換。
例二:
為第三桿(coupler)為驅動桿
[val,form]=f4bar([3242],0,60,10,0,-1,1)
val=
Columns1through3
300
1.3321-1.4919i-48.239-8.9487
2+3.4641i6010
0.33205+1.9722i80.44324.333
Columns4through6
01.3321-1.4919i2
-582.550.33205+1.9722i2
0-988.55-882.66i1325.3
496.46188.64-31.759i191.29
Column7
-48.239
80.443
-138.24
-9.5568
form=1
程式內容:
function[values,form]=f4bar(r,theta1,theta2,td2,tdd2,sigma,driver)
%
%function[values,form]=f4bar(r,theta1,theta2,td2,tdd2,sigma,driver)
%programdesignedbyDin-sueFon,NTU,revisedfromWaldron's
%Thisfunctionanalyzesafour-barlinkagewhenthecrankisthe
%drivinglink.Theinputvaluesare:
%theta1,theta2areanglesindegrees
%r
(1)=lengthofvector1(frame)
%r
(2)=lengthofvector2(crank)
%r(3)=lengthofvector3(coupler)
%r(4)=lengthofvector4(rockerorslideroffset)
%td2=crankorcouplerangularvelocity(rad/sec)
%tdd2=crankorcouplerangularacceleration(rad/sec^2)
%sigma=+1or-1.Identifiesassemblymode
%driver=0forcrankasdriver;1forcouplerasdriver
%Theresultsarereturnedinthevector"values".Theanswersare
%storedinvaluesaccordingtothefollowing:
%values(1:
4,1)=linkposition
%values(1:
4,2)=linkanglesindegrees
%values(1:
4,3)=linkangularvelocities
%values(1:
4,4)=linkangularaccelerations
%values(1,5)=velocityofpointQ
%values(2,5)=velocityofpointP
%values(3,5)=accelerationofpointQ
%values(4,5)=accelerationofpointP
%vakyes(4,6)=absolutevaluesofvalues(:
4)
%vakyes(4,7)=anglesindegreesofvalues(:
4)
%form=assemblyflag.Ifform=0,mechanismcannotbe
%assembled.
%convertinputdata
values=zeros(4,7);
%ifcoupleristhedriver,interchangethevetor3&2
%iftheta2>180|theta2<0,sigma=-sigma;end
ifdriver==1,
r=[r
(1)r(3)r
(2)r(4)];
end
rr=r.*r;
fact=pi/180;
theta=zeros(4,1);
td=zeros(4,1);
tdd=zeros(4,1);
theta(1:
2)=[theta1theta2]*fact;
t1=theta
(1);
tx=theta
(2);
s1=sin(t1);
c1=cos(t1);
sx=sin(tx);
cx=cos(tx);
%positioncalculations
A=2*r
(1)*r(4)*c1-2*r
(2)*r(4)*cx;
C=rr
(1)+rr
(2)+rr(4)-rr(3)-2*r
(1)*r
(2)*(c1*cx+s1*sx);
B=2*r
(1)*r(4)*s1-2*r
(2)*r(4)*sx;
arg=B*B-C*C+A*A;
if(arg>=0)
form=1;
%Checkforthedenominatorequaltozero
ifabs(C-A)>=eps
t4=2*atan((-B+sigma*sqrt(arg))/(C-A));
s4=sin(t4);
c4=cos(t4);
t3=atan2((r
(1)*s1+r(4)*s4-r
(2)*sx),(r
(1)*c1+r(4)*c4-r
(2)*cx));
s3=sin(t3);
c3=cos(t3);
elseifabs(C-A)%Ifthedenominatoriszero,computetheta(3)first
A=-2*r
(1)*r(3)*c1+2*r
(2)*r(3)*cx;
B=-2*r
(1)*r(3)*s1+2*r
(2)*r(3)*sx;
C=rr
(1)+rr
(2)+rr(3)-rr(4)-2*r
(1)*r
(2)*(c1*cx+s1*sx);
arg=B*B-C*C+A*A;
if(arg>=0)
t3=2*atan((-B-sigma*sqrt(arg))/(C-A));
s3=sin(t3);
c3=cos(t3);
t4=atan2((-r
(1)*s1+r(3)*s3+r
(2)*sx),(-r
(1)*c1+r(3)*c3+r
(2)*cx));
s4=sin(t4);
c4=cos(t4);
end
end
theta(3)=t3;
theta(4)=t4;
%velocitycalculation
td
(2)=td2;
AM=[-r(3)*s3,r(4)*s4;-r(3)*c3,r(4)*c4];
BM=[r
(2)*td
(2)*sx;r
(2)*td
(2)*cx];
CM=AM\BM;
td(3)=CM
(1);
td(4)=CM
(2);
%accelerationcalculation
tdd
(2)=tdd2;
BM=[r
(2)*tdd
(2)*sx+r
(2)*td
(2)*td
(2)*cx+r(3)*td(3)*td(3)*c3-r(4)*td(4)*td(4)*c4;...
r
(2)*tdd
(2)*cx-r
(2)*td
(2)*td
(2)*sx-r(3)*td(3)*td(3)*s3+r(4)*td(4)*td(4)*s4];
CM=AM\BM;
tdd(3)=CM
(1);
tdd(4)=CM
(2);
%storeresultsinarrayvalues
%coordinatesofPandQ
ifdriver==1,
r=[r
(1)r(3)r
(2)r(4)];
c2=c3;c3=cx;s2=s3;s3=sx;
td=[td
(1)td(3)td
(2)td(4)];
tdd=[tdd
(1)tdd(3)tdd
(2)tdd(4)];
theta=[theta
(1)theta(3)theta
(2)theta(4)];
else
c2=cx;s2=sx;
end
forj=1:
4,
values(j,1)=r(j).*exp(i*theta(j));
values(j,2)=theta(j)/fact;
values(j,3)=td(j);
values(j,4)=tdd(j);
end%positionvectors
values(1,5)=r
(2).*exp(i*theta
(2));%velocityforpointQ
values(2,5)=r(4).*exp(i*theta(4));%velocityforpointP
values(3,5)=i*r
(2).*(tdd
(2)-td
(2).*td
(2)).*exp(i*theta
(2));%accelofQ
values(4,5)=i*r(4).*(tdd(4)-td(4).*td(4)).*exp(i*theta(4));%accelofP
forj=1:
4,
values(j,6)=abs(values(j,5));%absolutevaluesforvalues(:
4)
values(j,7)=angle(values(j,5))/fact;%anglesforvalues(:
4)
end
%findtheaccelerations
else
form=0;
ifdriver==1,
r=[r
(1)r(3)r
(2)r(4)];
forj=1:
4,values(j,1)=r(j).*exp(i*theta(j));end%positionvectors
end
end
二、drawlinks函數
drawlinks之目的在利用MATLAB繪製四連桿之相關位置,其本身會呼叫f4bar.m函數以計算四連桿之向量位置,然後繪圖。
其呼叫格式如下:
functiondrawlinks(r,th1,th2,sigma,driver)
其輸入各式與f4bar.m大體相同,茲說明如下:
.r(1:
4)=各桿之長度,r
(1)為固定桿,其餘分別為曲桿、結合桿及被動桿。
.theta1=第一桿之水平角,或為四連桿之架構角,以角度表示。
.theta2=驅動桿之水平夾角,以角度表示。
一般為曲桿角,但若為結合桿驅動,則為結合桿之水平夾角。
.sigma=+1or-1.組合模式,負值表示閉合型,正值為分支型,但有時需視實際情況而定。
.driver=0(驅動桿為第二桿);1(驅動桿為第三桿)
例三、第二桿為驅動桿
drawlinks([3242],0,60,-1,0)
圖二、四連桿之繪圖
其繪出之四連桿為如圖二。
黑色為第一桿,藍色為第二桿,紅色為第三桿,綠色為第四桿。
例二、第三桿為結合桿(coupler)
drawlinks([3242],0,60,-1,1)
圖三、以結合桿為驅動桿(r=[3242])
圖三即為所得之答案,此時四連桿為分支型(branch),因為目前之情況無法轉為閉合型,即使將sigma值變號,仍為分支型,如:
clf;drawlinks([3242],0,60,1,1)
圖四、當sigma=1時並以結合桿為驅動桿(r=[3242])
利用drawlinks亦可繪出各種角度之圖型,可以作為四連桿運動過程之觀察,相當方便,例如:
clf;fori=10:
20:
360,drawlinks([1334],0,i,-1,0);end
圖五、多重位置之四連桿運動情形(r=[1334])
其結果如圖五,當然若以第三桿為驅動桿時,亦可獲得同樣的結果,例如:
clf;fori=10:
20:
360,drawlinks([1334],0,i,-1,1);end
圖六、以第三桿為驅動桿之情形(r=[1334])
當sigma變號時,亦可看出其不同的轉動方式,如下:
clf;fori=10:
20:
360,drawlinks([1334],0,i,1,1);end
圖七、當sigma變號(為正時)(r=[1334])
drawlinks程式內容
functiondrawlinks(r,th1,th2,sigma,driver)
%
%functiondrawlinks(r,th1,th2,sigma,driver)
%drawthepositionsoffour-barlinks
%willcallf4bar.mfuncion
%designedbyDin-SueFon,NTU
%
%r:
rowvectorforfourlinks
%th1:
frameangle
%th2:
crankangleorcoupleangle
%sigma:
assemblymode
%driver:
0forcrank,1forcoupler
%clf;
[rb]=f4bar(r,th1,th2,0,0,sigma,driver);
r(3,1)=r(1,1)+r(4,1);
rx=real(r(:
1));rx(4)=0;
ry=imag(r(:
1));ry(4)=0;
ifb==1
plot([0rx
(1)],[0ry
(1)],'k-','LineWidth',4);
holdon;
ifdriver==0
plot([0rx
(2)],[0ry
(2)],'b-','LineWidth',1.5);
plot([rx
(2)rx(3)],[ry
(2)ry(3)],'r-','LineWidth',2);
else
plot([0rx
(2)],[0ry
(2)],'r-','LineWidth',2);
plot([rx
(2)rx(3)],[ry
(2)ry(3)],'b-','LineWidth',1.5);
end
plot([rx
(1)rx(3)],[ry
(1)ry(3)],'-g');
plot(rx,ry,'bo');
text(0,0,'O');text(rx
(1),ry
(1),'R');
text(rx
(2),ry
(2),'P');text(rx(3),ry(3),'Q');
else
fprintf('Combinationoflinksfailatdegrees%6.1f\n',th2);
end
axisequal
gridon
三、drawlimits函數
四連桿之迴轉過程,能完全迴轉的情況仍然很少,有些時候無法獲得完整一圈的迴轉。
亦即依葛列夫定理四連桿之第一或第二類類型決定,前者為完整迴轉型,後者則有迴轉角度之限制,這些限制因四連桿長度決定之。
四連桿迴轉過程中,有可能其中兩桿會連成一線,或重疊成一線,前者若成立時,即變成三角形,後者若重疊時,亦會構成另一個三角形。
理論上連桿構成三角形應不會有相對運動。
故可稱為四連桿之運動極限。
由這兩個極端位置,可以知道四連桿之最終運動限制。
A.第二桿為驅動桿
在數學上,表示這兩個狀況之方法可以利用下列二種不等式進行測試:
r1+r2|r1-r2|>|r3-r4|
而由其不等式之方向,可以構成四種狀況,並進而求得該狀況之角度。
下面為第二桿為驅動桿時之四種情況:
(1)當r1+r2≦r3+r4,|r1-r2|≧|r3-r4|時
θmin=0,θmax=2π
(2)當r1+r2≧r3+r4,|r1-r2|≧|r3-r4|時
θmin=-cos-1{[r12+r22-(r3+r4)2]/(2r1r2)},
θmax=cos-1{[r12+r22-(r3+r4)2]/(2r1r2)}
(3)當r1+r2≧r3+r4,|r1-r2|≦|r3-r4|時
θmin=cos-1{[r12+r22-(r3-r4)2]/(2r1r2)},
θmax=cos-1{[r12+r22-(r3+r4)2]/(2r1r2)}
(4)當r1+r2≦r3+r4,|r1-r2|≦|r3-r4|時
θmin=cos-1{[r12+r22-(r3-r4)2]/(2r1r2)},
θmax=2π-cos-1{[r12+r22-(r3-r4)2]/(2r1r2)}
B.第三桿為驅動桿
第三桿結合桿為驅動桿時,則仍然取決於四連桿屬於葛列斯荷(Grashof)一型或二型。
若屬一型連桿,則當第三桿r3為最短桿時,第三桿可以作360度迴轉。
其餘之限制條件雖不如以第二桿為驅動桿者,但其極;限狀況是當第二桿與第四桿相平行時,變成無法繼續迴轉,除非它是處於平行四邊形。
將四種情況依下列二不等式之情況加以分類,在這些分類中,若兩式均為等號時,則應歸屬於第五類:
r1+r3|r1-r3|>|r2-r4|
(5)當r1+r3≦r2+r4,|r1-r3|≧|r2-r4|時
θmin=0,θmax=2π
(6)當r1+r3≧r2+r4,|r1-r3|≧|r2-r4|時
θmin=-cos-1{[r12+r32-(r2+r4)2]/(2r1r3)},
θmax=cos-1{[r12+r32-(r2+r4)2]/(2r1r3)}
(7)當r1+r3≧r2+r4,|r1-r3|≦|r2-r4|時
θmin=cos-1{[r12+r32-(r2-r4)2]/(2r1r3)},
θmax=cos-1{[r12+r32-(r2+r4)2]/(2r1r3)}
(8)當r1+r3≦r2+r4,|r1-r3|≦|r2-r4|時
θmin=cos-1{[r12+r32-(r2-r4)2]/(2r1r3)},
θmax=2π-cos-1{[r12+r32-(r2-r4)2]/(2r1r3)}
C.fb_angle_limits函數
觀察上面討論之四個極限角度,可以寫一組程式進行計算。
由於以第三桿驅動與第二桿驅動,在計算上僅是將其中之r2與r3之位置對調即可。
為尋找上述極限角度θmin、θmax,可用函數fb_angle_limits進行尋找,其格式如下:
function[Qstart,Qstop]=fb_angle_limits(r,Q1,driver)
其中輸入項目有:
r=四連桿之長度向量,其定義與前函數相同。
Q1=第一桿之夾角,角度表示(deg)。
driver=驅動模式(=0第二桿驅動;=1第三桿驅動)。
而輸出項為兩個角度:
Qstart=驅動桿(第二桿或第三桿)之最低限角度(deg)
Qstop=