四连杆之MATLAB相关程式.docx

四連桿之 MATLAB 相關程式-26-台大生機系 馮丁樹教授4-25-2002第三章中之四連桿分析可以參考相關資料，本節則針對四連桿之動作程式加 以說明。目前所設計之程式有 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/sec2)。sigma=+1 or -1.組合模式，負值表示閉合型，正值為分支型，但有 時需視實際情況而定。driver=0 (驅動桿為第二桿); 1 (驅動桿為第三桿)輸出變數：form=組合狀態， 0 ：表示無法組合； 1：可以正確組合values=輸出矩陣，其大小為 4 X 7，各行之資料分配如下：12(deg)3(rad/s)4(rad/s2)567I桿 1 位置111VQ|VQ|VQII桿 2 位置222VP|VP|VPIII桿 3 位置333AQ|AQ|AQIV桿 4 位置444AP|AP|AP其中第一行之連桿位置向量，屬於單桿的位置向量，以格式以複數表示。第二 行為各桿之水平夾角，以度表示；第三及第四行為各桿之角度速度及角加速度， 以單位時間之弧度表示。第五至七行則為 P 點與 Q 點之速度與加速度量，第五 行為向量，第六行為絕對量，第七行為夾角，以度數表示。值得一提的是第一行之向量表示法屬於複數之型式，故若要得到其絕對值僅 需在 MATLAB 指令檔中，以 abs()這一個函數指令即可求得，而以函數 angle() 則可求得其夾角，雖然第二行與第七行之輸出亦有相對應之夾角。例一：為第二桿為驅動桿，r=3 2 4 2，theta1=0;theta2=60;td2=10; tdd2=0;sigma=-1;driver=0(crank)val,form=f4bar(3 2 4 2,0,60,10,0,-1,0)val =Columns 1 through 33001 +1.7321i60103.8682 -1.0182i-14.7465.40781.8682 +0.71389i20.91316.549Columns 4 through 601 +1.7321i201.8682 +0.71389i2-127.58173.21 -100i200-236.27364.19 -953.09i1020.3Column 76020.913-30-69.087form =1 (表示可以組合)本例中，有框線者表示其為輸入值。但第一行則已經轉換為複數型式。未來 複數型式要轉為 x-y 座表時，只要使用函數 real()及 imag()兩指令，即可進行轉 換。例二：為第三桿(coupler)為驅動桿，r=3 2 4 2，theta1=0; theta2=60; td2=10; tdd2=0; sigma=-1; driver=1(crank)val,form=f4bar(3 2 4 2,0,60,10,0,-1,1)val =Columns 1 through 33001.3321 -1.4919i-48.239-8.94872 +3.4641i60100.33205 +1.9722i80.44324.333Columns 4 through 601.3321 -1.4919i2-582.550.33205 +1.9722i20-988.55 -882.66i1325.3496.46188.64 -31.759i191.29Column 7-48.23980.443-138.24-9.5568form =1程式內容：f4bar.mfunction values,form = f4bar(r,theta1,theta2,td2,tdd2,sigma,driver)%function values,form = f4bar(r,theta1,theta2,td2,tdd2,sigma,driver)% program designed by Din-sue Fon, NTU, revised from Waldrons% This function analyzes a four-bar linkage when the crank is the% driving link. The input values are:% theta1,theta2 are angles in degrees%r(1)= length of vector 1 (frame)%r(2)= length of vector 2 (crank)%r(3)= length of vector 3 (coupler)%r(4)= length of vector 4 (rocker or slider offset)%td2= crank or coupler angular velocity (rad/sec)%tdd2= crank or coupler angular acceleration (rad/sec2)%sigma= +1 or -1. Identifies assembly mode%driver = 0 for crank as driver; 1 for coupler as driver% The results are returned in the vector values. The answers are% stored in values according to the following:%values (1:4,1)= link position%values (1:4,2) = link angles in degrees%values (1:4,3)= link angular velocities%values (1:4,4)= link angular accelerations%values (1,5) = velocity of point Q%values (2,5) = velocity of point P%values (3,5) = acceleration of point Q%values (4,5) = acceleration of point P%vakyes (4,6)=absolute values of values(:,4)%vakyes (4,7)=angles in degrees of values(:,4)%form= assembly flag. If form = 0, mechanism cannot be%assembled.%convert input datavalues=zeros(4,7);% if coupler is the driver, interchange the vetor 3 & 2% if theta2180|theta2=0)form=1;% Check for the denominator equal to zeroif abs(C-A)=epst4=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);elseif abs(C-A)=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);endendtheta(3)=t3; theta(4)=t4;%velocity calculationtd(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=AMBM;td(3)=CM(1);td(4)=CM(2);%acceleration calculationtdd(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)*c 4;.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=AMBM;tdd(3)=CM(1);tdd(4)=CM(2);%store results in array values% coordinates of P and Qif driver=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);elsec2=cx;s2=sx;endfor j=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 % position vectors values(1,5)=r(2).*exp(i*theta(2);%velocity for point Q values(2,5)=r(4).*exp(i*theta(4);%velocity for point Pvalues(3,5)=i*r(2).*(tdd(2)-td(2).*td(2).*exp(i*theta(2);%accel of Qvalues(4,5)=i*r(4).*(tdd(4)-td(4).*td(4).*exp(i*theta(4);%accel of Pfor j=1:4,values(j,6)=abs(values(j,5); %absolute values for values(:,4)values(j,7)=angle(values(j,5)/fact; %angles for values(:,4)end%find the accelerationselseform=0;if driver=1,r=r(1) r(3) r(2) r(4);for j=1:4, values(j,1)=r(j).*exp(i*theta(j);end % position vectorsendend二、drawlinks 函數有了 f4bar.m 之計算，只要連桿組可以組立(form = 1)，即表示每一個相對應 連桿的角度均可獲得，因而可以進行繪出其相關位置。程式 drawlinks 之目的即 是利用 MATLAB 繪製四連桿之相關位置。故此程式本身會呼叫 f4bar.m 函數以 計算四連桿之向量位置，然後繪圖。在其呼叫 f4bar.m 函數時，其角速度及角加 速度均設為零。其呼叫格式如下：function drawlinks(r,th1,th2,sigma,driver)其輸入各式與 f4bar.m 大體相同，茲說明如下：r(1:4) = 各桿之長度，r(1)為固定桿，其餘分別為曲桿、結合桿及被動桿。theta1 = 第一桿之水平角，或為四連桿之架構角，以角度表示。theta2 = 驅動桿之水平夾角，以角度表示。一般為曲桿角，但若為結合桿驅動，則為結合桿之水平夾角。sigma =+1 or -1.組合模式，負值表示閉合型，正值為分支型，但有時 需視實際情況而定。driver =0 (驅動桿為第二桿); 1 (驅動桿為第三桿)程式中，在呼叫 f4bar.m 過後，僅取其 values 中第一欄中之數據，即 r1、r2、 r3 及 r4 等之複數型式資料。構成基本的 O、R、P、Q 等四點。O 為原點，R 及 Q 為 r2 與 r4 的向量點，直接使用即可，而 P 點則為兩桿組合點，即為 r1+r4 或 r2+r3 的向量和，本程式取其前者。例三、第二桿為驅動桿drawlinks(3 2 4 2,0,60,-1,0)2.52P1.51Q0.50 OR-0.5-100.511.522.533.544.5圖二、四連桿之繪圖其繪出之四連桿為如圖二。黑色為第一桿，藍色為第二桿，紅色為第三桿， 綠色為第四桿。例二、第三桿為結合桿(coupler)drawlinks(3 2 4 2,0,60,-1,1)Q1.510.50OR-0.5-1P-0.500.511.522.533.5圖三、以結合桿為驅動桿(r =3 2 4 2)圖三即為所得之答案，此時四連桿為分支型(branch)，因為目前之情況無法轉 為閉合型，即使將 sigma 值變號，仍為分支型，如：clf；drawlinks(3 2 4 2,0,60,1,1)Q10.50OR-0.5-1-1.5P-0.500.511.522.533.5圖四、當 sigma=1 時並以結合桿為驅動桿(r =3 2 4 2)利用 drawlinks 亦可繪出各種角度之圖型，可以作為四連桿運動過程之觀察， 相當方便，例如：將 theta2 改以每 20 度繪一次，由繪至 360 度。此時必須使用 迴圈的繪製方式來達成：clf; for i=10:20:360,drawlinks(1 3 3 4,0,i,-1,0); endQQQQ3QPPPQQ PP2QPPQ1QPP0ORQQPP-1PPQ-2QPPQPP-3PQQQ-4-3-2-10123456圖五、多重位置之四連桿運動情形(r =1 3 3 4)其結果如圖五，當然若以第三桿(聯結桿)為驅動桿時，亦可獲得同樣的結 果，此時其他數值暫不改變，僅就參數最後一項改為 1，例如：clf;for i=10:20:360,drawlinks(1 3 3 4,0,i,-1,1);endQQQQ3QPPPPQPQPP QPQ1PQPP0ORPQQ-1PPQP-2QPQP-3PQQQ-4-3-2-10123456圖六、以第三桿為驅動桿之情形(r =1 3 3 4)圖五與圖六不同之處在於第二桿(曲桿)的角度分佈。前者係由曲桿作均勻分 佈，故可以看出 360 度的均勻分佈；而圖六則為不均均的分佈情況，但第三桿的 角度變化則應為均勻的，雖然無法立即由圖六看得出來。當 sigma 變號時，亦可看出其不同的轉動方式，如下：clf; for i=10:20:360,drawlinks(1 3 3 4,0,i,1,1); endQQQ3PQ PP2QP QP1PQPQ0ORPQP-1PQPQPPQ PPQPP-3QPQ QQQ-4-3-2-10123456圖七、當 sigma 變號(為正時) (r =1 3 3 4)drawlinks 程式內容function drawlinks(r,th1,th2,sig ma,driver)%function drawlinks(r,th1,th2,sigma,driver)%draw the positions of four-bar links%will call f4bar.m funcion%designed by Din-Sue Fon, NTU%r:row vector for four links%th1: frame angle%th2: crank angle or couple angle%sigma: assembly mode%driver: 0 for crank, 1 for coupler%clf;r b=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; if b=1plot(0 rx(1),0 ry(1),k-,LineWidth,4); hold on;if driver=0plot(0 rx(2),0 ry(2),b-,LineWidth,1.5);plot(rx(2) rx(3),ry(2) ry(3),r-,LineWidth,2); elseplot(0 rx(2),0 ry(2),r-,LineWidth,2);plot(rx(2) rx(3),ry(2) ry(3),b-,LineWidth,1.5);endplot(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); elsefprintf(Combination of links fail at degrees %6.1fn,th2);endaxis equal grid on三、drawlimits 函數四連桿之迴轉過程，能完全迴轉的情況仍然很少，有些時候無法獲得完整一 圈的迴轉。亦即依葛列夫定理四連桿之第一或第二類類型決定，前者為完整迴轉 型，後者則有迴轉角度之限制，這些限制因四連桿長度決定之。四連桿迴轉過程中，有可能其中兩桿會連成一線，或重疊成一線，前者若成 立時，即變成三角形，後者若重疊時，亦會構成另一個三角形。理論上連桿構成 三角形應不會有相對運動，故可稱為四連桿之運動極限。由這兩個極端位置，可 以知道四連桿之最終運動限制。A. 第二桿為驅動桿時在數學上，表示這兩個狀況之方法可以利用下列二種不等式進行測試：r1+r2|r3-r4|而由其不等式之方向，可以構成四種狀況，並進而求得該狀況之角度。下面 為第二桿為驅動桿時之四種情況：(1) 當 r1+r2r3+r4，|r1-r2|r3-r4| 時 min =0 ，max =2r3r4r2r1(2) 當 r1+r2r3+r4，|r1-r2|r3-r4| 時-1222min = - cosr1 +r2 -(r3+r4) /(2 r1 r2)，-1222max =cosr1 +r2 -(r3+r4) /(2 r1 r2)r3r2Amaxr4r1Amin(3) 當 r1+r2r3+r4，|r1-r2|r3-r4| 時-1222min =cosr1 +r2 -(r3-r4) /(2 r1 r2)，-1222max =cosr1 +r2 -(r3+r4) /(2 r1 r2)r3r2AmaxAminr4r1(4) 當 r1+r2r3+r4，|r1-r2|r3-r4| 時-1222min =cosr1 +r2 -(r3-r4) /(2 r1 r2) ，-1222max = 2- cosr1 +r2 -(r3-r4) /(2 r1 r2)AmaxAminr2r3r1r4B. 第三桿為驅動桿第三桿結合桿為驅動桿時，則仍然取決於四連桿屬於葛列斯荷(Grashof)一型或二 型。若屬一型連桿，則當第三桿 r3 為最短桿時，第三桿可以作 360 度迴轉。其 餘之限制條件雖不如以第二桿為驅動桿者，但其極;限狀況是當第二桿與第四桿 相平行時，變成無法繼續迴轉，除非它是處於平行四邊形。將四種情況依下列二 不等式之情況加以分類，在這些分類中，若兩式均為等號時，則應歸屬於第五類：r1+r3|r2-r4|(5) 當 r1+r3r2+r4，|r1-r3|r2-r4| 時 min =0 ，max =2r3r4r2r1(6) 當 r1+r3r2+r4，|r1-r3|r2-r4| 時-1222min = - cosr1 +r3 -(r2+r4) /(2 r1 r3) ，-1222max =cosr1 +r3 -(r2+r4) /(2 r1 r3)r2r4Aminr1r3Amax(7) 當 r1+r3r2+r4，|r1-r3|r2-r4| 時-1222min =cosr1 +r3 -(r2-r4) /(2 r1 r3) ，-1222max =cosr1 +r3 -(r2+r4) /(2 r1 r3)r2r4Aminr1r3Amax(8) 當 r1+r3r2+r4，|r1-r3|r2-r4| 時-1222min =cosr1 +r3 -(r2-r4) /(2 r1 r3) ，-1222max = 2- cosr1 +r3 -(r2-r4) /(2 r1 r3)Aminr1r3r4Amaxr2C. 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 =驅動桿(第二桿或第三桿)之最高限角度 (deg)drawlimits 函數則是呼叫 fb_angle_limits.m 函數，然後將其極限位置繪出。其輸 入項目與 drawlinks 函數相同。D.執行例(第二桿為驅動桿)(1) drawlimits(1 2 3 4,0,-1,0)Qstart =3.6e-005 Qstop = 360OsRs12=03.060.0pP1.510.50qQ-0.5-1-1.500.511.522.533.544.5圖八、第二桿驅動，屬全迴轉之情形(r = 1 2 3 4)本例中，無論第二桿或第四桿，均可完全迴轉。(2) drawlimits(3 5 2 1,0,1,0)Qstart =-33.557Qstop =33.557p 2.521.5s2=33.61q0.50OR-0.5-1Q-1.5s1=-33.6-2-2.5P-1012345圖九、第二類型，右邊有限制，左邊無限(r =3 5 2 1)圖九則因 3+52+1，|3-5|2-1|，故屬前述之第二型限制角度。 (3)drawlimits(5 4 1 3,0,1,0)Qstart =22.332Qstop =51.3183.5p 32.5q Q21.5s2=51.3P1s1=22.30.50 OR00.511.522.533.544.55圖十、第三類型，前後角度限制(r =5 4 1 3)本例中，因為各桿長度符合第三類型，即 5+41+3，|5-4|1-3|，故其角度 限制如圖十。(4)drawlimits(4 3 3 5,0,1,0)Qstart =28.955Qstop =331.04Q32P1s1=29.00ORs2=331.0-1p-2-3q-2-10123456圖十一、第四類型(r2 驅動) (r =4 3 3 5)此類型與第一類型類似，但角度限制在右邊，因為各桿之長度符合該項條 件：4+33+5，|4-3|3-5|。如圖十一。E.執行例(第三桿為驅動桿)(5) drawlimits(4 5 3 5,0,-1,1)Qstart =3.6e-005Qstop =360pPqQ4.543.532.5ss12=03.600.021.510.50OR-1012345圖十二、第三桿與第一桿平行，但可完全迴轉(r =4 5 3 5)(6) drawlimits(4 3 3 3,0,-1,1)Qstart =-117.28Qstop =117.281.5qP10.5s1=-117.30 OR-0.5s2=117.3-1Qp-1.500.511.522.533.54圖十三、 第三桿之限制角度相反，屬第六情況(r =4 3 3 3 )(7) drawlimits(4 4 6 1,0,-1,1)Qstart =26.384Qstop =55.771q0.50OR-0.5-1Q-1.5-2s1=26.4s2=55.8-2.5-3-3.5Pp-2-101234圖十四、不同角度但同一象限內之限制，屬第七狀況(r =4 4 6 1)(8) drawlimits(3 4 4 6,0,-1,1)Qstart =28.955Qstop =331.04Q 54P32s1=29.010OR-1-2s2=331.0-3-4p-5q-4-202468圖十五、上下象限之限制角度，兩者和為 360 度(r =3 4 4 6 )E.程式內容function Qstart, Qstop=fb_angle_limits(r,Q1,driver)%function Qstart, Qstop=fb_angle_limits(r,Q1,driver)% Function file that calculates the initial and final angle for four-bar% driver=0 when the crank is the driver.% driver=1 when coupler is the driver.% Variables%r=linkage row vector (cm)%Q1=frame angle(degree);%Qstart=initial crank(or coupler) angle (deg)%Qstop=final crank(or coupler) angle (deg)%Programif driver=1, r=r(1) r(3) r(2) r(4);endr1=r(1);r2=r(2);r3=r(3);r4=r(4);rmin=min(r); rmax=max(r); rtotal=sum(r);if (r1+r2)(r3+r4)& abs(r1-r2)(r3+r4)& abs(r1-r2)=abs(r3-r4) Qstart=acos(r22-(r4+r3)2+r12)/(2*r1*r2); Qstop=-Qstart;endif (r1+r2)=abs(r3-r4) Qstart=0;Qstop=2*pi;endif (r1+r2)=(r3+r4)&abs(r1-r2)abs(r3-r4) Qstart=acos(r22-(r4-r3)2+r12)/(2*r1*r2); Qstop=2*pi-Qstart;endrange = Qstop-Qstart; Qstart=Q1+(Qstart+0.0000001*range)*180/pi; Qstop=Q1+(Qstop-0.0000001*range)*180/pi; if QstopQstart;temp=Qstart; Qstart=Qstop; Qstop=temp;end%function drawlimits(r,th1,sigma,driver)%function drawlmits(r,th1,sigma,driver)%draw the positions of four-bar links%call f4