四连杆之MATLAB程式讲解.docx
《四连杆之MATLAB程式讲解.docx》由会员分享，可在线阅读，更多相关《四连杆之MATLAB程式讲解.docx（30页珍藏版）》请在冰豆网上搜索。
四连杆之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=驅動桿（第二桿或第三桿）之最低限角度（d
升级会员 