萬花尺matlab仿真（圓內旋輪線，異形齒輪）

• 0 前言
• 1 圓形齒輪，單點
• 2 圓形齒輪，圖形孔
• 3 異形齒輪，單點
• 4 異形齒輪，異形孔

0 前言

萬花尺是一種常見的玩具，通常由兩個齒輪組成。大齒輪固定，小齒輪在大齒輪內側運動。小齒輪上有許多空隙點，將筆尖插入空隙點轉動小齒輪，即可畫出很漂亮的圖形。

（網上隨便找的圖，如有問題請告知↑）

看著圖形很好看，所以準備借此水一篇。接下來嘗試用matlab進行仿真，模擬萬花尺生成的曲線

1 圓形齒輪，單點

首先簡單的繪制一個草圖，建立數學模型：

內部小齒輪的公轉角度設為$\theta$，小齒輪的自轉角度設為$?$。假設小齒輪的自轉角速度為$\omega$，則根據理論力學知識，綠色的點為瞬心，小齒輪質心處的速度為$\omega r$。
其中，R為大齒輪半徑，r為小齒輪半徑。
小齒輪質心處的公轉角速度可以計算為：${\omega }_c=\frac{\omega r}{R?r}$

從地面坐標系上看，小圓質心可以用公轉角度$\theta$和公轉角速度${\omega }_c$表示，${\theta }_0$為初始時刻小圓的角度。
$${\omega }_c=\frac{\omega r}{R?r}$$
$$\theta ={\theta }_0+{\omega }_{c}t$$
因此，以大圓的圓心作為坐標原點，建立xy直角坐標系：
$$x_c=R\cos (\theta )$$
$$y_c=R\sin (\theta )$$
小圓初始的自轉角度為${?}_0$，則時間t后，自轉的角度$?$為：
$$?={?}_0?\omega t$$
因此，圓盤上的一點，可以用小圓上的質心位置，加上一個旋轉后的相對位置得到：
$$x_b=x_c+b\cos (\theta )$$
$$y_b=y_c+b\sin (\theta )$$

b為圓上一點到圓盤質心的距離。

利用得到的$x_b$，$y_b$，我們可以得到一個參數方程，用來繪制這個點的軌跡。我們把這個軌跡稱為圓內旋輪線。

旋輪線的周期與大圓和小圓之間的比例相關：
如果大圓比小圓等于2：1，則代表一個公轉周期內，小圓自轉2圈，此時軌跡圖案為一條直線；
如果大圓比小圓等于3：1，則代表一個公轉周期內，小圓自轉3圈；
如果大圓比小圓等于5：2，則代表2個公轉周期內，小圓自轉5圈。

當偏心率等于(R-r)/r時，每個自轉周期內，旋輪線都會通過原點。

接下來用matlab編寫程序演示：

clear clc close all%設定參數 %圖1，周期演示 lamda=0.7;%偏心比例,等于b/r period=[20,41];%第一個數代表公轉周期，第二個數代表自轉周期，互質。等于r/R。 [xbl,ybl]=hypocycloid(lamda,period); figure(1);plot(xbl,ybl) xlim([-1,1]);ylim([-1,1]) axis equal %圖2,全部通過中心演示 lamda=(13-11)/11;%偏心比例,等于b/r period=[11,13];%第一個數代表公轉周期，第二個數代表自轉周期，互質。等效于[r,R]。 [xbl,ybl]=hypocycloid(lamda,period); figure(2);plot(xbl,ybl) %xlim([-1,1]);ylim([-1,1]) axis equal %圖3,圖1的彩色版 lamda=0.7; period=[20,41]; [xbl,ybl]=hypocycloid(lamda,period); figure(3); colormap(spring) patch([xbl;NaN],[ybl;NaN],-[sqrt(xbl.^2+ybl.^2);NaN],'EdgeColor','interp','Marker','none',...'MarkerFaceColor','flat','LineWidth',1); xlim([-1,1]);ylim([-1,1]) axis equal box on %圖4，偏心率大于1時 lamda=1.1;%偏心比例,等于b/r period=[20,41];%第一個數代表公轉周期，第二個數代表自轉周期，互質。等于r/R。 [xbl,ybl]=hypocycloid(lamda,period); figure(4);plot(xbl,ybl) xlim([-1,1]);ylim([-1,1]) axis equalfunction [xbl,ybl]=hypocycloid(lamda,period) %繪制內旋 %幾何參數 R=1;%大圓半徑 r=R*period(1)/period(2);%小圓半徑 b=lamda*r;%點的偏心位置 w=1;%角速度（隨意） theta0=0;%小圓初始位置 phi0=0;%小圓初始旋轉角度tl=0:0.01:period(2)*2*pi; N=length(tl); xbl=zeros(N,1); ybl=zeros(N,1);for k=1:Nt=tl(k);%小圓圓心角度wc=w*r/(R-r);theta=theta0+wc*t;%小圓圓心位置xc=(R-r)*cos(theta);yc=(R-r)*sin(theta);%小圓上一點的角度phi=phi0-w*t;%小圓上一點的位置xb=xc+b*cos(phi);yb=yc+b*sin(phi);%記錄xbl(k)=xb;ybl(k)=yb; end end

下圖為輸出的結果，可以看到效果是很不錯的。當然如果仔細設計，再利用不同的顏色組合繪制，應該還可以達到更好的效果。

2 圓形齒輪，圖形孔

有時候，小齒輪內部并不是單純的幾個點，而是存在很多其它圖形。這時候如果再去仿真筆尖下的線條，就不能直接套用旋輪線公式去分析。

因此，提出一個假設，認為無論什么樣的圖形，在用筆繪制并使得小齒輪旋轉的時候，筆尖位置一定在最遠離大圓中心的位置附近。比如下圖，紅色的點是認為最遠離中心的點，也是筆尖所在的點。

因此，我修改程序，得到了小齒輪上，任意圖形下的萬花尺筆尖軌跡仿真。當然，也不是任意圖形，需要保證滿足“最遠端”的點在任意時刻具有唯一性，因為不滿足的話，在實際使用中會因為受力原因，出現小齒輪反轉的情況，這樣實際繪圖的話也會出現困擾。

clear clc close all%萬花尺上的圖形繪圖 %設定參數 period=[3,11];%第一個數代表公轉周期，第二個數代表自轉周期，互質。等于r/R。%幾何參數 R=1;%大圓半徑 r=R*period(1)/period(2);%小圓半徑 w=1;%角速度（隨意） theta0=0;%小圓初始位置 phi0=0;%小圓初始旋轉角度%設定圖形，對于小圓圓心的坐標 xf=r/4+r/3*cos(0:0.01:2*pi); yf=0+r/8*sin(0:0.01:2*pi);tl=0:0.01:period(1)*2*pi*(R-r)/r; N=length(tl); xbl=zeros(N,1); ybl=zeros(N,1);figure() for k=1:Nt=tl(k);%小圓圓心角度wc=w*r/(R-r);theta=theta0+wc*t;%小圓圓心位置xc=(R-r)*cos(theta);yc=(R-r)*sin(theta);%小圓上一點的角度phi=phi0-w*t;%小圓上圖形的坐標xyf2=[cos(phi),-sin(phi);sin(phi),cos(phi)]*[xf;yf]+[xc;yc];xf2=xyf2(1,:);yf2=xyf2(2,:);%求圖形上最遠離大圓圓心的那一個點[~,I]=max(xf2.^2+yf2.^2);xf2_i=xf2(I);yf2_i=yf2(I);%記錄xfl(k)=xf2_i;yfl(k)=yf2_i;%小圓標記線xyf_temp=[cos(phi),-sin(phi);sin(phi),cos(phi)]*[0,r;0,0]+[xc;yc];%繪圖if mod(k,5)==0clfxlim([-1,1])ylim([-1,1])axis equalhold onrectangle('Position',[-1,-1,2,2],'Curvature',[1 1])%繪制大齒輪rectangle('Position',[xc-r,yc-r,2*r,2*r],'Curvature',[1 1])%繪制小齒輪plot(xf2,yf2);%繪制圖形plot(xfl,yfl)%繪制點的軌跡plot(xyf_temp(1,:),xyf_temp(2,:))%標記線，證明內齒輪在轉動hold offpause(0.01)end end

輸出圖像如下：

3 異形齒輪，單點

有時候，萬花尺的某些小齒輪甚至不是圓的，比如下面這些：

其中，凸的光滑圖形可以在圓上滾動，當局部曲率半徑小于大圓曲率半徑時，每一個點都可以與大圓接觸。但是對于凹的圖形，凹陷的部分一定不可能與大圓接觸。

因此，對于不能接觸的部分，在運動過程中等效于不接觸面積補全后，補全圖形的運動。比如上圖十字圖形，接觸補全后如右圖所示。比如上圖的三角形，下面直線部分曲率半徑無窮大，所以運動式也等效于在下面補全圓弧的圖形。

因此，我們擁有把任意形狀的齒輪等效變形為，可以讓每個點都與大圓接觸并一一對應的新圖形。之后，利用接觸弧長l作為參數，對異形齒輪進行建模。

為了驗證思路的正確性，以橢圓在直線上滾動為例。橢圓滾動的距離利用長度l進行計算，橢圓滾動的角度利用接觸點切線來計算。

結果看上去沒有問題。

之后，對圓形齒輪內部的運動進行建模，思路同樣，但是接觸角度利用大齒輪的切線和小齒輪的切線共同確定。

clear clc close all%異形齒輪圓周 %接觸面x^2+y^2=R^2 R=1;%外側固定齒輪半徑%首先，把齒輪轉換為填充后的等效形式%齒輪上采點（沿著運動方向） % xr=0.12970964*cos(0:0.001:2*pi); % yr=0.25941928*sin(0:0.001:2*pi);%周期1/5%xr=0.25941928*cos(0:0.001:2*pi); %yr=0.51883856*sin(0:0.001:2*pi);%周期2/5xr=0.185299485714286*cos(0:0.001:2*pi); yr=0.370598971428571*sin(0:0.001:2*pi);%周期2/7%定義繪圖點 xyp=[0.05;0.1]; xyp_l=[];%計算每一點對應的長度 lr=sqrt((xr(2:end)-xr(1:end-1)).^2+(yr(2:end)-yr(1:end-1)).^2); lr=cumsum(lr); lr=[lr,lr(end)+sqrt((xr(end)-xr(1)).^2 + (xr(end)-xr(1)).^2)]; %2*pi/lr(end)-5%初始接觸點(對于齒輪來說） x0=xr(1); y0=yr(1); %初始接觸點旋轉角度對x軸 phi0=angle((xr(2)-xr(end))+(yr(2)-yr(end))*1i);figure(1) for l=0:0.005:4*pil_G=l;%地面坐標系下，運動的長度%移動周長l后的接觸點位置%l=0.01;if l>lr(end)l=l-fix(l/lr(end))*lr(end);endk=find(l<lr,1);if k==1phik=angle((xr(1)-xr(end))+(yr(1)-yr(end))*1i);phi=-(phik-phi0);else%接觸點的旋轉角度phik=angle((xr(k)-xr(k-1))+(yr(k)-yr(k-1))*1i);phi=-(phik-phi0);end%接觸點壁面角度theta=l_G/R;%接觸點位置（地面坐標，利用固定接觸面計算）xk_G=R*cos(theta);yk_G=R*sin(theta);%齒輪繞接觸點旋轉phi2=phi+theta;xyrk_G=[cos(phi2),-sin(phi2);sin(phi2),cos(phi2)]*([xr;yr]-[xr(k);yr(k)])+[xk_G;yk_G];xr_G=xyrk_G(1,:);yr_G=xyrk_G(2,:);%繪圖點做同樣的旋轉xyp_G=[cos(phi2),-sin(phi2);sin(phi2),cos(phi2)]*(xyp-[xr(k);yr(k)])+[xk_G;yk_G];%儲存繪圖點xyp_l=[xyp_l,xyp_G];if mod(k,5)==0clfhold onplot(xr_G,yr_G)rectangle('Position',[-1,-1,2,2],'Curvature',[1 1])%大圓plot(xyp_l(1,:),xyp_l(2,:));%異形齒輪set(gca,'Units','centimeters')set(gca,'position',[1,1,9,9])ylim([-1,1]);xlim([-1,1])hold offbox on%axis equalpause(0.01)end endfunction L=Dis1(xy1,xy2) L=sqrt((xy1(1)-xy2(1))^2+(xy1(2)-xy2(2))^2); end

最終的繪圖效果如下：

4 異形齒輪，異形孔

到了最復雜的組合，齒輪是異形齒輪的，孔也異形的，怎么辦？

這就需要結合第2章和第3章的代碼，來解決問題：

matlab代碼如下：

clear clc close all%異型齒輪圓周 %接觸面x^2+y^2=R^2 R=1;%外側固定齒輪半徑%首先，把齒輪轉換為填充后的等效外形%齒輪上采點（沿著運動方向） xr=0.12970964*cos(0:0.001:2*pi); yr=0.25941928*sin(0:0.001:2*pi);%周期1/5%設定圖形，對于異型齒輪的坐標 xf=0+0.04*cos(0:0.001:2*pi); yf=0.1+0.08*sin(0:0.001:2*pi); xyf=[xf;yf];%定義繪圖點 xyp_l=[];%計算每一點對應的長度 lr=sqrt((xr(2:end)-xr(1:end-1)).^2+(yr(2:end)-yr(1:end-1)).^2); lr=cumsum(lr); lr=[lr,lr(end)+sqrt((xr(end)-xr(1)).^2 + (xr(end)-xr(1)).^2)]; %2*pi/lr(end)-5%初始接觸點(對于齒輪來說） x0=xr(1); y0=yr(1); %初始接觸點旋轉角度對x軸 phi0=angle((xr(2)-xr(end))+(yr(2)-yr(end))*1i);figure(1) for l=0:0.005:4*pil_G=l;%地面坐標系下，運動的長度%移動周長l后的接觸點位置%l=0.01;if l>lr(end)l=l-fix(l/lr(end))*lr(end);endk=find(l<lr,1);if k==1phik=angle((xr(1)-xr(end))+(yr(1)-yr(end))*1i);phi=-(phik-phi0);else%接觸點的旋轉角度phik=angle((xr(k)-xr(k-1))+(yr(k)-yr(k-1))*1i);phi=-(phik-phi0);end%接觸點壁面角度theta=l_G/R;%接觸點位置（地面坐標，利用固定接觸面計算）xk_G=R*cos(theta);yk_G=R*sin(theta);%齒輪繞接觸點旋轉phi2=phi+theta;xyrk_G=[cos(phi2),-sin(phi2);sin(phi2),cos(phi2)]*([xr;yr]-[xr(k);yr(k)])+[xk_G;yk_G];xr_G=xyrk_G(1,:);yr_G=xyrk_G(2,:);%繪圖圖形做同樣的旋轉，圖形的坐標xyf_G=[cos(phi2),-sin(phi2);sin(phi2),cos(phi2)]*(xyf-[xr(k);yr(k)])+[xk_G;yk_G];xf_G=xyf_G(1,:);yf_G=xyf_G(2,:);%求圖形上最遠離大圓圓心的那一個點[~,I]=max(xf_G.^2+yf_G.^2);xf_G_i=xf_G(I);yf_G_i=yf_G(I);%儲存繪圖點xy_G_i=[xf_G(I);yf_G(I)];xyp_l=[xyp_l,xy_G_i];if true%mod(k,5)==0clfhold onplot(xr_G,yr_G)rectangle('Position',[-1,-1,2,2],'Curvature',[1 1])plot(xyp_l(1,:),xyp_l(2,:));plot(xf_G,yf_G);set(gca,'Units','centimeters')set(gca,'position',[1,1,9,9])ylim([-1,1]);xlim([-1,1])hold offbox on%axis equalpause(0.01)end endfunction L=Dis1(xy1,xy2) L=sqrt((xy1(1)-xy2(1))^2+(xy1(2)-xy2(2))^2); end

總結

以上是生活随笔為你收集整理的万花尺matlab仿真（圆内旋轮线，异形齿轮）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。