第三章-四连杆之MATLAB程式讲解.doc

第三章 四連桿之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/sec2)。 .sigma = +1 or -1. 組合模式，負值表示閉合型，正值為分支型，但有時需視實際情況而定。 .driver = 0 (驅動桿為第二桿); 1 (驅動桿為第三桿) 輸出變數：.form = 組合狀態， 0 ：表示無法組合； 1：可以正確組合 .values = 輸出矩陣，其大小為 4 X 7，各行之資料分配如下： 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(3 2 4 2,0,60,10,0,-1,0)val = Columns 1 through 3 3 0 0 1 + 1.7321i 60 10 3.8682 - 1.0182i -14.746 5.4078 1.8682 + 0.71389i 20.913 16.549 Columns 4 through 6 0 1 + 1.7321i 2 0 1.8682 + 0.71389i 2 -127.58 173.21 - 100i 200 -236.27 364.19 - 953.09i 1020.3 Column 7 60 20.913 -30 -69.087 form = 1 (表示可以組合)本例中，有框線者表示其為輸入值。但第一行則已經轉換為複數型式。未來複數型式要轉為x-y座表時，只要使用函數real()及imag()兩指令，即可進行轉換。例二：為第三桿(coupler)為驅動桿val,form=f4bar(3 2 4 2,0,60,10,0,-1,1)val = Columns 1 through 3 3 0 0 1.3321 - 1.4919i -48.239 -8.9487 2 + 3.4641i 60 10 0.33205 + 1.9722i 80.443 24.333 Columns 4 through 6 0 1.3321 - 1.4919i 2 -582.55 0.33205 + 1.9722i 2 0 -988.55 - 882.66i 1325.3 496.46 188.64 - 31.759i 191.29 Column 7 -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)% 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)*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=AMBM;tdd(3)=CM(1); tdd(4)=CM(2); %store results in array values % coordinates of P and Q if 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); else c2=cx;s2=sx; end for 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 P values(3,5)=i*r(2).*(tdd(2)-td(2).*td(2).*exp(i*theta(2);%accel of Q values(4,5)=i*r(4).*(tdd(4)-td(4).*td(4).*exp(i*theta(4);%accel of P for 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 accelerationselse form=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 vectors endend二、drawlinks函數drawlinks之目的在利用MATLAB繪製四連桿之相關位置，其本身會呼叫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 (驅動桿為第三桿) 例三、第二桿為驅動桿drawlinks(3 2 4 2,0,60,-1,0)圖二、四連桿之繪圖其繪出之四連桿為如圖二。黑色為第一桿，藍色為第二桿，紅色為第三桿，綠色為第四桿。例二、第三桿為結合桿(coupler)drawlinks(3 2 4 2,0,60,-1,1)圖三、以結合桿為驅動桿(r =3 2 4 2)圖三即為所得之答案，此時四連桿為分支型(branch)，因為目前之情況無法轉為閉合型，即使將sigma值變號，仍為分支型，如：clf；drawlinks(3 2 4 2,0,60,1,1)圖四、當sigma=1時並以結合桿為驅動桿(r =3 2 4 2)利用drawlinks亦可繪出各種角度之圖型，可以作為四連桿運動過程之觀察，相當方便，例如：clf;for i=10:20:360,drawlinks(1 3 3 4,0,i,-1,0);end圖五、多重位置之四連桿運動情形(r =1 3 3 4)其結果如圖五，當然若以第三桿為驅動桿時，亦可獲得同樣的結果，例如：clf;for i=10:20:360,drawlinks(1 3 3 4,0,i,-1,1);end圖六、以第三桿為驅動桿之情形(r =1 3 3 4)當sigma 變號時，亦可看出其不同的轉動方式，如下：clf;for i=10:20:360,drawlinks(1 3 3 4,0,i,1,1);end圖七、當sigma變號(為正時) (r =1 3 3 4)drawlinks程式內容function drawlinks(r,th1,th2,sigma,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=1 plot(0 rx(1),0 ry(1),k-,LineWidth,4); hold on; if driver=0 plot(0 rx(2),0 ry(2),b-,LineWidth,1.5); plot(rx(2) rx(3),ry(2) ry(3),r-,LineWidth,2); else plot(0 rx(2),0 ry(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(Combination of links fail at degrees %6.1fn,th2);endaxis equalgrid on三、drawlimits函數四連桿之迴轉過程，能完全迴轉的情況仍然很少，有些時候無法獲得完整一圈的迴轉。亦即依葛列夫定理四連桿之第一或第二類類型決定，前者為完整迴轉型，後者則有迴轉角度之限制，這些限制因四連桿長度決定之。四連桿迴轉過程中，有可能其中兩桿會連成一線，或重疊成一線，前者若成立時，即變成三角形，後者若重疊時，亦會構成另一個三角形。理論上連桿構成三角形應不會有相對運動。故可稱為四連桿之運動極限。由這兩個極端位置，可以知道四連桿之最終運動限制。A. 第二桿為驅動桿在數學上，表示這兩個狀況之方法可以利用下列二種不等式進行測試：r1+r2|r3-r4|而由其不等式之方向，可以構成四種狀況，並進而求得該狀況之角度。下面為第二桿為驅動桿時之四種情況：(1) 當r1+r2r3+r4，|r1-r2|r3-r4| 時min =0 ，max =2(2) 當r1+r2r3+r4，|r1-r2|r3-r4| 時min = - cos-1r12+r22-(r3+r4)2/(2 r1 r2)，max = cos-1r12+r22-(r3+r4)2/(2 r1 r2)(3) 當r1+r2r3+r4，|r1-r2|r3-r4| 時min = cos-1r12+r22-(r3-r4)2/(2 r1 r2)，max = cos-1r12+r22-(r3+r4)2/(2 r1 r2)(4) 當r1+r2r3+r4，|r1-r2|r3-r4| 時min = cos-1r12+r22-(r3-r4)2/(2 r1 r2) ，max = 2- cos-1r12+r22-(r3-r4)2/(2 r1 r2)B. 第三桿為驅動桿第三桿結合桿為驅動桿時，則仍然取決於四連桿屬於葛列斯荷(Grashof)一型或二型。若屬一型連桿，則當第三桿r3為最短桿時，第三桿可以作360度迴轉。其餘之限制條件雖不如以第二桿為驅動桿者，但其極;限狀況是當第二桿與第四桿相平行時，變成無法繼續迴轉，除非它是處於平行四邊形。將四種情況依下列二不等式之情況加以分類，在這些分類中，若兩式均為等號時，則應歸屬於第五類：r1+r3|r2-r4|(5) 當r1+r3r2+r4，|r1-r3|r2-r4| 時min =0 ，max =2(6) 當r1+r3r2+r4，|r1-r3|r2-r4| 時min = - cos-1r12+r32-(r2+r4)2/(2 r1 r3) ，max = cos-1r12+r32-(r2+r4)2/(2 r1 r3)(7) 當r1+r3r2+r4，|r1-r3|r2-r4| 時min = cos-1r12+r32-(r2-r4)2/(2 r1 r3) ，max = cos-1r12+r32-(r2+r4)2/(2 r1 r3)(8) 當r1+r3r2+r4，|r1-r3|r2-r4| 時min = cos-1r12+r32-(r2-r4)2/(2 r1 r3) ，max = 2- cos-1r12+r32-(r2-r4)2/(2 r1 r3)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 = 驅動桿(第二桿或第三桿)之最高限角度 (deg)drawlimits函數則是呼叫fb_angle_limits.m函數，然後將其極限位置繪出。其輸入項目與drawlinks函數相同。D.執行例(第二桿為驅動桿)(1) drawlimits(1 2 3 4,0,-1,0)Qstart = 3.6e-005Qstop = 360圖八、第二桿驅動，屬全迴轉之情形(r = 1 2 3 4)本例中，無論第二桿或第四桿，均可完全迴轉。(2) drawlimits(3 5 2 1,0,1,0)Qstart = -33.557Qstop = 33.557圖九、第二類型，右邊有限制，左邊無限(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.318圖十、第三類型，前後角度限制(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.04圖十一、第四類型(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 = 360圖十二、第三桿與第一桿平行，但可完全迴轉(r =4 5 3 5)(6)drawlimits(4 3 3 3,0,-1,1)Qstart = -117.28Qstop = 117.28圖十三、 第三桿之限制角度相反，屬第六情況(r =4 3 3 3 )(7)drawlimits(4 4 6 1,0,-1,1)Qstart = 26.384Qstop = 55.771圖十四、不同角度但同一象限內之限制，屬第七狀況(r =4 4 6 1)(8)drawlimits(3 4 4 6,0,-1,1)Qstart = 28.955Qstop = 331.04圖十五、上下象限之限制角度，兩者和為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 f4bar funcion%r: row vector for four links%th1: frame angle%sigma: assembly mode%driver: 0 for crank, 1 for couplerQstart, Qstop=fb_angle_limits(r,th1,dr