5/16- 第九次作業

作業九 廖婉婷 B94611040


*本週(5/3)有來上課*

由課本上定義,滑塊之滑動方向若與水平方向有一距離e時
稱為偏置機構,而其中e稱之偏置量
而此偏置機構可視為一個四連桿結構
其中偏置量將影響滑桿機構之活動範圍
<本題著重在極限位置之分析>

首先依題目給定R=50，L=55，e=10
由課本之slider_limit函式
計算滑塊水平滑行最大距離(即衝程)和界限角度
以下為程式內容:

function [ s theta21,theta22 ] = slider_limit( R,L,e )
%R= 曲桿長
%L= 連桿長
%e= 偏置量
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;

執行後得到以下結果:

 slider_limit( 50,55,10 )
abs(s)
abs(theta21)
abs(theta22)

ans =

  1.0452e+002 -8.6603e+000i


ans =

  104.8809


ans =

    5.4650


ans =

  117.4463

直接對各項取絕對值
我們得到衝程為104.88
界限角度為5.47和117.45

討論:
在這邊我先發現了個問題
由於篇置量固定
最大衝程應該出現在曲桿和連桿為一直線時
由畢氏定理計算出值應為104.52?





接著利用老師講義上所提供繪製極限位置之函式(drawsldlimits)
此程式除呼叫sld_angle_limits獲得迴轉之範圍外
呼叫drawsldlink函數繪出各桿之位置。
以下為函式內容:

function [Qstart, Qstop]=sld_angle_limits(r,theta1,linkdrive)
%輸入函數:
%r = 四連桿之長度向量
%theta1 = 第一桿之夾角，角度表示(deg)
%linkdrive = 驅動模式(=0 第二桿驅動； =1 第三桿驅動 =2 滑塊驅動)
%輸出(兩個角度)
%Qstart = 驅動桿(第二桿或第三桿)之最低限角度 (deg)
%Qstop = 驅動桿(第二桿或第三桿)之最高限角度 (deg)
r1=r(1);r2=r(2);r3=r(3);r4=r(4);
g2d=180/pi;
switch linkdrive
case 0 %crank
  if r3+r4=0 & r3>r4
     Qstart=asin((r4-r3)/r2);Qstop=asin((r4+r3)/r2);
  elseif r3+r4>=r2 & r3>=r4 & r3>=0
     Qstart=asin((r4-r3)/r2);Qstop=pi-asin((r4-r3)/r2);
  elseif r3-r4<=r2 & r4<0 & r3>=-r4
     Qstart=asin((r4-r3)/r2);Qstop=asin((r4+r3)/r2);
  elseif r3-r4>=r2 & r3>=0 & -r4>=r3
     Qstart=-pi-asin((r4+r3)/r2);Qstop=asin((r4+r3)/r2);
  else
     Qstart=0;Qstop=2*pi;
  end
case 1 %coupler
   if r2-r4<=r3 & r4>=0 & r2>=r4
       Qstart=asin((r4-r2)/r3);Qstop=pi-asin((r4-r2)/r3);
   elseif r2+r4=0
       Qstart=asin((r4-r2)/r3);Qstop=asin((r4+r2)/r3);
   elseif r2+r4<=r3 & r4<=0 & r2+r4>=0
       Qstart=-pi-asin((r4+r2)/r3);
       Qstop=asin((r4+r2)/r3);
   elseif r2-r4 < r3 & r4 < =0
       Qstart=asin((r4-r2)/r3);Qstop=asin((r4+r2)/r3);
   else %r2>=(r3+abs(r4))
       Qstart=0;Qstop=2*pi;
   end
case 2 %slider displacement
  Qstart=0;Qstop=0;
  arg2=(r2+r3)^2-r4^2;
  if abs(r2-r3)>=r4
     arg1=(r2-r3)^2-r4^2;
     if arg1>0,Qstart=sqrt(arg1);end;
     Qstop=sqrt(arg2);
  else
     if arg2<0, return; end
     Qstart=sqrt(arg2);Qstop=-sqrt(arg2);
  end
  theta1=0;g2d=1;
end %case
if Qstop < Qstart,TT=Qstart;Qstart=Qstop;Qstop=TT;end
adjust=(Qstop-Qstart)*1e-8;
Qstart=theta1+(Qstart+adjust)*g2d;
Qstop=theta1+(Qstop-adjust)*g2d;

function [values]=drawlinks(r,th1,th2,mode,linkdrive)
%r = 四桿之向量位置
%th1 = 第一桿水平角度
%th2 = 驅動桿水平角度
%mode = +1 or -1 組合模數,負值表閉合型;正值表分支型
%linkdriver = 0 (驅動桿為第二桿); 1 (驅動桿為第三桿) 
if nargin<5,linkdrive=0;end
if nargin<4,mode=1;end
[values b]=f4bar(r,th1,th2,0,0,mode,linkdrive);
rr=values(:,1);
rr(3,1)=rr(1,1)+rr(4,1);
rx=real(rr(:,1));rx(4)=0;
ry=imag(rr(:,1));ry(4)=0;
if b==1
plot([0 rx(1)],[0 ry(1)],'k-','LineWidth',4);
hold on;
if linkdrive==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-','LineWidth',1.5);
plot(rx,ry,'bo');
text(0,0,' O');text(rx(1),ry(1),' R');
text(rx(2),ry(2),' Q');text(rx(3),ry(3),' P');
else
fprintf('Combination of links fail at degrees %6.1f\n',th2);
end
axis equal
grid on

function drawsldlimits(r,th1,sigma,driver)
%r = 四桿之向量位置
%th1 = 第一桿水平角度
%sigma = +1 or -1 組合模數,負值表閉合型;正值表分支型
%driver = 0 (驅動桿為第二桿); 1 (驅動桿為第三桿) ;2(滑塊驅動)
g2d=180/pi;
[Qstart, Qstop]=sld_angle_limits(r,th1,driver)
if Qstart==0 & Qstop==0, return; end
switch driver
case 0 %Link 2 as driving link
  [rr]=drawsldlinks(r,th1,Qstart,sigma,driver);
  text(real(rr(2,1)/2),imag(rr(2,1)/2),['s1=' num2str(Qstart,'%6.1f')]);
  [rr]=drawsldlinks(r,th1,Qstop,sigma,driver);
  text(real(rr(2,1)/2),imag(rr(2,1)/2),['s2=' num2str(Qstop,'%6.1f')]);
case 1 %Link 3 as driving link
  [rr]=drawsldlinks(r,th1,Qstart,sigma,driver);
  text(real(rr(2,1)+rr(3,1)/2),imag(rr(2,1)+rr(3,1)/2),...
            ['s1=' num2str(Qstart,'%6.1f')]);
  [rr]=drawsldlinks(r,th1,Qstop,sigma,driver);
  text(real(rr(2,1)+rr(3,1)/2),imag(rr(2,1)+rr(3,1)/2),...
            ['s1=' num2str(Qstop,'%6.1f')]);
case 2 %Slide as driving link
  if abs(r(2)-r(3)) < r(4);
     angx=asin((r(2)+r(3))/r(4))*g2d;
     tmax=th1+angx;
     tmin=th1+180-angx;
  else
     tmin=th1+asin((r(2)+r(3))/r(4))*g2d;
     tmax=th1+asin((r(2)-r(3))/r(4))*g2d;
  end
  r(1)=Qstart;offr=0.1*Qstart;
  [rr]=drawsldlinks(r,th1,tmin,sigma,driver);
  text(real(rr(1,1)/2),imag(rr(1,1)/2)+offr,['r1=' num2str(Qstart,'%6.1f')]);
  r(1)=Qstop;
  [rr]=drawsldlinks(r,th1,tmax,sigma,driver);
  text(real(rr(1,1)/2),imag(rr(1,1)/2)-offr,['r1=' num2str(Qstop,'%6.1f')]);
end
axis equal
grid on

雖然題目只給定R,L,e之值 也就是相當於四連桿之r2,r3,r4
由於前面slider_limit已求得衝程為104.88 也就是相當於四連桿之r1
因此在drawsldlimits函式中
即可直接代入四桿之值

(1)曲桿驅動(以第二桿為驅動桿)
因為曲桿驅動
driver為0

>> axis([-110 110 -40 40])
drawsldlimits([104.88 50 55 10],0,-1,0)

Qstart =

  3.6000e-006


Qstop =

  360.0000

討論:
在這邊我發現了一個問題
照程式結果跑出來的角度可作360度旋轉
但是根據課本上所述:
"若曲柄須能完成360度的迴轉，則其相關連桿之長度須依其所能夠成三角形關係而定
而若有偏置的情形時，則連結桿與曲柄之長度差應大於偏置距離e"
與本題條件所給R-L<0條件不符,顯然有錯誤
在請教過歐陽太閒同學後
發現是sld_angle_limits程式碼中少了r3>=r4>0這種情況(也就是講義所列第二情形)
因此將原程式碼中:

   elseif r3+r4>=r2 & r4>=r3 & r3>=0
Qstart=asin((r4-r3)/r2);Qstop=pi-asin((r4-r3)/r2);

修正為:

   elseif r3+r4>=r2 & r3>=r4 & r3>=0
Qstart=asin((r4-r3)/r2);Qstop=pi-asin((r4-r3)/r2);

重新跑過一次後:

axis([-110 110 -40 40])
drawsldlimits([104.88 50 55 10],0,-1,0)
title('第二桿驅動極限位置 -b94611040')

Qstart =

  -64.1581


Qstop =

  244.1581

得到了正確的結果





(2)連桿驅動(以第三桿為驅動桿)

%drive為1
axis([-110 110 -40 40])
drawsldlimits([104.88 50 55 10],0,-1,1)
title('第三桿驅動極限位置 -b94611040')

Qstart =

  -46.6582


Qstop =

  226.6582

(3)滑塊驅動

%driver為2
>>  axis([-110 110 -40 40])
drawsldlimits([104.88 50 55 10],0,-1,2)
title('滑塊驅動極限位置 -b94611040')

Qstart =

 -104.5227


Qstop =

  104.5227

需特別注意
由於滑塊驅動屬於直線運動
故在程式中之Qstart及Qstop值為第一桿位移之最大與最小量。

討論
做到這邊 一直在思考r1的長度問題
原本是使用slider_limit計算出之最大衝程直接當作r1
但是後來考慮到此為滑塊驅動情況
當為曲稈驅動時會有角度的限制
因此試的代入其他值至r1 發現結果一樣
後來觀察sld_angle_limits程式碼
發現限制角度運算都是使用r2,r3和r4
r1會隨運動角度不同而改變 因此設任何值似乎並不影響


接著分別就第二及第三桿驅動之運動狀態作成動畫:
(1)曲桿驅動(以第二桿為驅動桿)

for i=[-64.1581:244.1581]
title('第二桿驅動極限位置 -b94611040')
pause(0.005);
clf;
drawsldlinks([0 50 55 10],0,i,1,0)
axis([-60 120 -80 80]);
axis equal;
end;

(2)連桿驅動(以第三桿為驅動桿)

for i=[-46.6582:226.6582]
title('第三桿驅動極限位置 -b94611040')
pause(0.005);
clf;
drawsldlinks([0 50 55 10],0,i,1,1)
axis([-100 100 -100 100]);
axis equal;
end;

如果要分析各桿之位置,速度,加速度
我們使用講義中的function sldlinks作分析
以下為程式碼內容:
由前面 已知r2=50, r3=55, r4=10,
r1值隨滑塊運動作變動,在此先設為1
另設驅動桿角速度為5, 角加速度為0

function [data,form] = sldlink(r,theta1,theta2,td2,tdd2,sigma,driver)
% r= [r1 r2 r3 r4],四桿之長
% td2= 驅動桿之角速度(rad/sec)
% tdd2= 驅動桿之角加速度(rad/sec^2)
% sigma= +1 or -1.  組合模數, 負值表示閉合型，正值為分支型
% driver= 0(驅動桿為第二桿); 1 (驅動桿為第三桿) ;2(滑塊驅動)
f=@(num,ndg) round(num*10^ndg)/10^ndg;
data=zeros(4,8);
if driver==1, r(2:3)=[r(3) r(2)];end
rr=r.*r;d2g=pi/180;
theta=[theta1 theta2 0 theta1+90]*d2g;
t1=theta(1);t4=theta(4);
[td,tdd,rd,rdd]=deal(zeros(4,1));
s1=sin(t1);c1=cos(t1);
s4=sin(t4);c4=cos(t4);
switch driver
case 2 % for the sliding block as driver
  rd(1)=td2;
  rdd(1)=tdd2;
  A=-2*r(1)*r(2)*c1-2*r(2)*r(4)*c4;
  B=-2*r(1)*r(2)*s1-2*r(2)*r(4)*s4;
  C=rr(1)+rr(2)+rr(4)-rr(3)+2*r(1)*r(4)*(c1*c4+s1*s4);
  arg=B*B-C*C+A*A;
  arg=f(arg,5);
  if (arg>=0)
    form=1; %assembly flag
    t2=2*atan((-B+sigma*sqrt(arg))/(C-A));
    s2=sin(t2);   c2=cos(t2);
    t3=atan2((r(1)*s1+r(4)*s4-r(2)*s2),(r(1)*c1+r(4)*c4-r(2)*c2));
    theta(2)=t2;  theta(3)=t3;
    s3=sin(t3);   c3=cos(t3);
    %velocity calculation
    AM=[-r(2)*s2, -r(3)*s3; r(2)*c2, r(3)*c3];
    BM=[rd(1)*c1;rd(1)*s1];
    CM=AM\BM;
    td(2)=CM(1); td(3)=CM(2);
    %acceleration calculation
    BM=[r(2)*td(2)*td(2)*c2+r(3)*td(3)*td(3)*c3+rdd(1)*c1;...
       r(2)*td(2)*td(2)*s2+r(3)*td(3)*td(3)*s3+rdd(1)*s1];
    CM=AM\BM;
    tdd(2)=CM(1); tdd(3)=CM(2);
  else
    form=0; return;
  end
case {0,1}  % for link 2 or 3 as a crank
  td(2)=td2; tdd(2)=tdd2; tx=theta(2);
  sx=sin(tx); cx=cos(tx);
  % position calculations
  A=2*r(4)*(c1*c4+s1*s4)-2*r(2)*(c1*cx+s1*sx);
  B=rr(2)+rr(4)-rr(3)-2*r(2)*r(4)*(cx*c4+sx*s4);
  arg=A*A-4*B;
  arg=f(arg,5);
  if (arg>=0)
    form=1; %assembly flag
    r(1)=(-A+sigma*sqrt(arg))/2;
    t3=atan2((r(1)*s1+r(4)*s4-r(2)*sx),(r(1)*c1+r(4)*c4-r(2)*cx));
    theta(3)=t3;
    s3=sin(t3);  c3=cos(t3);
    %velocity calculation
    AM=[c1, r(3)*s3; s1, -r(3)*c3];
    BM=[-r(2)*td(2)*sx;r(2)*td(2)*cx];
    CM=AM\BM;
    rd(1)=CM(1);  td(3)=CM(2);
    %acceleration calculation
    BM=[-r(2)*tdd(2)*sx-r(2)*td(2)*td(2)*cx-r(3)*td(3)*td(3)*c3;...
       r(2)*tdd(2)*cx-r(2)*td(2)*td(2)*sx-r(3)*td(3)*td(3)*s3];
    CM=AM\BM;
    rdd(1)=CM(1);  tdd(3)=CM(2);
  else
    form=0;
    if driver==1,
      r(2:3)=[r(3) r(2)];
      for j=1:4, data(j,1)=r(j).*exp(i*theta(j));end % position vectors
    end
    return;
  end
    if driver==1,
      r(2:3)=[r(3) r(2)];
      c2=c3;c3=cx;s2=s3;s3=sx;
      td(2:3)=[td(3) td(2)];
      tdd(2:3)=[tdd(3) tdd(2)];
      theta(2:3)=[theta(3) theta(2)];
    elseif driver==0
      c2=cx;s2=sx;
    end
end % end of switch
for j=1:4,
  data(j,1)=r(j).*exp(i*theta(j));
  data(j,2:4)=[theta(j)/d2g,td(j),tdd(j)];
end
data(:,5)=[rd(1);rdd(1);data(2,1);data(2,1)+data(3,1)];
data(1,6)=i*r(2)*td(2)*exp(i*theta(2));
data(2,6)=data(1,6)+i*r(3)*td(3)*exp(i*theta(3));
data(3,6)=r(2)*(i*tdd(2)-td(2)*td(2))*exp(i*theta(2));
data(4,6)=data(3,6)+r(3)*(i*tdd(3)-td(3)*td(3))*exp(i*theta(3));
%[Vq;Vp;Aq;Ap]
data(1:4,7)=abs(data(1:4,6));data(1:4,8)=angle(data(1:4,6))/d2g;

由預設條件代入函式執行得:

[values,form] = sldlink([0 50 55 10],0,100,5,0,-1,0)
     abs(values)

ans =
%  column1   column2   column3   column4
    0.0472         0         0         0    
    0.0500    0.1000    0.0050         0    
    0.0550    0.1345    0.0011    0.0307    
    0.0100    0.0900         0         0    

   column5   column6   column7    column8
    0.2904    0.2500    0.2500    0.1700
    0.9368    0.2904    0.2904    0.1800
    0.0500    1.2500    1.2500    0.0800
    0.0483    0.9368    0.9368    0.1800

其中第一行為各連桿位置向量。
第二行為各桿之水平夾角，第三及第四行為各桿之角度速度及角加速度。
第五行為第一桿之軸向速度與加速度、P及Q點之位置向量。
第六至八行則為P點與Q點之速度與加速度量，第六行為向量，第七行為絕對量，第八行為夾角。

此為以第二桿為軀動桿 且為預設角度之情況
因此若要以第三桿,連桿驅動, 則可比照此作法執行
另外,在嘗試將r1代入不同值時
我們發現在此處r1的值也不影響執行結果


此外,若要作連續性的速度/加速度分析
則可利用之前講義上的draw_4paths函式
由於在此不需作桿上特定點分析,因此r6,th6之值設為0
function draw_4paths(r6,th6,nlink,r,th1,theta,td2,tdd2,sigma,npts,driver,mode)
%Inputs:
% r= 四連桿之長度，以向量表示
% th1,theta: 第一桿及驅動桿之水平角度
% td2,tdd2: 驅動桿之角速度及角加速度
% sigma: 組合模數
% driver: 0為第二桿驅動, 1為第三桿驅動
% npts: 設定分割的點數或位置
% r6,rh6,nlink: 桿上特定點之位置，包括桿長，與桿之夾角及附於何桿
% mode=0 畫簡單位置圖
% =1 畫所有圖表
% =2 畫簡單位置圖+速度圖
% =3 畫所有圖表，但用簡單位置圖
代入預設值執行程式後:

>> draw_4paths(0,0,3,[0 50 55 10],0,30,2,5,1,100,0,1)

不過由於尚無法完全掌握此函式之內容
因此在瞭解清楚之前 還不敢亂加猜測其中意義