机械原理3-28MATLAB平面连杆机构运动分析,解三角函数超越方程

% while norm(f)>epsilon
J=[0 -x(6)*sin(theta2) x(7)*sin(theta3) -x(8)*sin(theta4);
0 x(6)*cos(theta2) -x(7)*cos(theta3) x(8)*cos(theta4);
cos(x(1)) 0 0 x(11)*sin(theta4); sin(x(1)) 0 0 -x(11)*cos(theta4)]; dth=inv(J)*(-1.0*f); lA=lA+dth(1); theta2=theta2+dth(2); theta3=theta3+dth(3); theta4=theta4+dth(4);
-0.0059 0.0105
0.0247 -1.2302
-0.0059 0.0109
0.0250 -1.2704
-0.0058 0.0114
0.0255 -1.3137
Columns 9 through 15
-0.0058 -0.0059 -0.0060 -0.0062 -0.0065 -0.0069 -0.0078 0.0119 0.0125 0.0131 0.0139 0.0148 0.0159 0.0175 0.0259 0.0265 0.0272 0.0281 0.0292 0.0306 0.0327 -1.3610 -1.4136 -1.4734 -1.5431 -1.6273 -1.7337 -1.8767 程序运行得到 q 矩阵，第一行到第三行分别是 a2、a3、a4 的值，第四行是杆 AF 上 滑块运动的速度，即 F 点的速度。
L2cosθ2
注意到，关于位置的四个方程组成的方程组是关于三角函数的非线性超越方 程。这里用牛顿——辛普森方法来求解。
第一步对位置方程进行求解： 用 SolidWorks 建模对机构中 AF 杆的极限位置进行分析，如图：
Hale Waihona Puke Baidu
得出 AF 杆的运动范围是-57.91°~59°
根据位置方程式编制如下 rrrposi.m 函数： function y=rrrposi(x) % %script used to implement Newton-Raphson mechod for %solving nonlinear position of RRR bar group % %Input parameters %x(1)=theta-1 %x(2)=theta-2 guess value %x(3)=theta-3 guess value %x(4)=theta-4 guess value %x(5)=l1 %x(6)=l2 %x(7)=l3 %x(8)=l4 %x(9)=l6 %x(10)=lA guess value %x(11)=lB % %Output parameters % %y(1)=lA %y(2)=theta-2 %y(3)=theta-3 %y(4)=theta-4 % theta2=x(2); theta3=x(3); theta4=x(4);
(1)位置方程
L2cosθ2-L3cosθ3-L4cos(θ4+180°)=-L1cos(θ1+180°)+L6
L2sinθ2-L3sinθ3-L4sin(θ4+180°)=-L1sin(θ1+180°)
-L’4cosθ4+ L’1cosθ1=L6
- L’4sinθ4+ L’1sinθ1=0
Xc=L1cos(θ1+180°)+L2cosθ2
L4cosθ4
0
a3
0
0
L’4sinθ4
cosθ1
a4
0
0
-L’4cosθ4
sinθ1
L’’1
=-
- w2 L2cosθ2 w3 L3cosθ3 -w4L4cosθ4
-w2 L2sinθ2 - w3 L3sinθ3
-w4L4sinθ4
0
0
w4L’4cosθ4
0
0
w4L’4sinθ4
w2
w1L1cosθ1
平面连杆机构运动学分析
3-15
已知：LAE=70mm,LAB=40mm,LEF=60mm,LDE=35mm,LCD=75mm,LBC=50mm,
原动件以等角速度 w1=10rad/s 回转。试以图解法求在θ1=50°时 C 点的速度 VC
和
加
速
度
ac.
先对机构进行位置分析： 由封闭形 ABCDEA 与 AEFA 有： L1+L2=L6+L3+L4 L’1=L6+L’4 即 L2-L3-L4=-L1+L6 -L’4+L’1=L6
x(10)*cos(x(2)) -x(11)*cos(x(3)) x(14)*cos(x(4)) 0; 0 0 x(15)*sin(x(4)) cos(x(1)); 0 0 -x(15)*cos(x(4)) sin(x(1))]; B=[-x(10)*x(6)*cos(x(2)) x(11)*x(7)*cos(x(3)) -x(12)*x(8)*cos(x(4)) 0; -x(10)*x(6)*sin(x(2)) -x(11)*x(7)*sin(x(3)) -x(12)*x(8)*sin(x(4)) 0; 0 0 x(15)*x(8)*cos(x(4)) -x(5)*sin(x(1)); 0 0 x(15)*x(8)*sin(x(4)) x(5)*cos(x(1))];
Yc= L1sin(θ1+180°) +L2 sinθ2
（2）速度方程
-L2sinθ2 L3sinθ3 L4cos(θ4+180°)
0
w2
L2cosθ2 -L3cosθ3 -L4cos(θ4+180°)
0
w3
0
0
L’4sinθ4
cosθ1
w4
0
0
-L’4cosθ4
sinθ1
L’1
L1sin(θ1+180°)
第三步进行加速度计算： 编写加速度计算函数 rrra.m：
function y=rrra(x) % %Input parameters % %x(1)=th1 %x(2)=th2 %x(3)=th3 %x(4)=th4 %x(5)=dth1 %x(6)=dth2 %x(7)=dth3 %x(8)=dth4 %x(9)=r1 %x(11)=r3 %x(12)=r4 %x(13)=r6 %x(14)=lA %x(15)=lB %x(16)=V % %y(1)=ddth2 %y(2)=ddth3 %y(3)=ddth4 %y(4)=a % A=[-x(10)*sin(x(2)) x(11)*sin(x(3)) -x(14)*sin(x(4)) 0;
x(7)*cos(x(2)) -x(8)*cos(x(3)) -x(9)*cos(x(4)+pi) 0; 0 0 x(12)*sin(x(4)) cos(x(1)); 0 0 -x(12)*cos(x(4)) sin(x(1))]; B=[x(6)*sin(x(1)+pi);-x(6)*cos(x(1)+pi);x(11)*sin(x(1));x(11)*cos(x(1))]*x(5); y=inv(A)*B;
lA=x(10) % epsilon=1.0E-6; % f=[x(6)*cos(theta2)-x(7)*cos(theta3)-x(8)*cos(pi+theta4)+ x(5)...
*cos(x(1)+pi)-x(9);
x(6)*sin(theta2)-x(7)*sin(theta3)-x(8)*sin(theta4+pi)+... x(5)*sin(x(1)+pi); -x(11)*cos(theta4)+lA*cos(x(1))-x(9); -x(11)*sin(theta4)+lA*sin(x(1))];
%x(8)=l3 %x(9)=l4 %x(10)=l6 %x(11)=lA %x(12)=lB % %Outout parameters % %y(1)=V %y(2)=dtheta-2 %y(3)=dtheta-3 %y(4)=dtheta-4 % A=[-x(7)*sin(x(2)) x(8)*sin(x(3)) x(9)*sin(pi+x(4)) 0;
根据第一步得到的数据进行数据输入，运行程序计算各速度值。程序如下：
x2=[x1' p(:,2) p(:,3) p(:,4) 10*ones(15,1) 40*ones(15,1) 50*ones(15,1)... 75*ones(15,1) 35*ones(15,1) 70*ones(15,1) p(:,1) 60*ones(15,1)]; q=zeros(4,15); for m=1:15 y2=rrrvel(x2(m,:)); q(:,m)=y2; end q
p= 93.3149 0.7163 2.5455 1.5461 91.3071 0.7045 2.5617 1.5902 89.2387 0.6929 2.5786 1.6347 87.1076 0.6815 2.5963 1.6796 84.9113 0.6703 2.6147 1.7250 82.6463 0.6592 2.6339 1.7709 80.3086 0.6482 2.6539 1.8174 77.8931 0.6372 2.6747 1.8646 75.3930 0.6263 2.6965 1.9126 72.7998 0.6154 2.7192 1.9616 70.1019 0.6043 2.7431 2.0118 67.2833 0.5930 2.7683 2.0635 64.3217 0.5812 2.7950 2.1169 61.1835 0.5687 2.8237 2.1728 57.8153 0.5551 2.8549 2.2319
w3
w1L1 sinθ1
w4
+ w1 w1L’1cosθ1+v sinθ1
L’1
w1L’1 sinθ1+v cosθ1
0 0 -w1sinθ1 -w1cosθ1
acx= w12L1cosθ1 -
w2 2
L2cosθ2
-
a2L2sinθ2
acy= w12L1 sinθ1 -
w2 2
L2
sinθ2
+a2
norm(f); end; y(1)=lA; y(2)=theta2; y(3)=theta3; y(4)=theta4;
再进行数据输入，运行程序进行运算。这里我们根据上面分析的θ1 的极限 位置取θ1 的范围为 40°~55°并均分成 15 个元素： clc
clear x1=linspace(40*pi/180,55*pi/180,15); x=zeros(length(x1),11); for n=1:15 x(n,:)=[x1(:,n) pi/6 8*pi/9 2*pi/3 40 50 75 35 70 75 60]; end p=zeros(length(x1),4); for k=1:15 y= rrrposi(x(k,:)); p(k,:)=y; end >> p
f=[x(6)*cos(theta2)-x(7)*cos(theta3)-x(8)*cos(pi+theta4)+ x(5)...
*cos(x(1)+pi)-x(9);
x(6)*sin(theta2)-x(7)*sin(theta3)-x(8)*sin(theta4+pi)+... x(5)*sin(x(1)+pi); -x(11)*cos(theta4)+lA*cos(x(1))-x(9); -x(11)*sin(theta4)+lA*sin(x(1))];
输出的 P、矩阵的第二列到第四列分别是θ2 、θ3 、θ4 的值，第一列是 AF 杆的长度 L1’。
第二步进行速度计算： 根据速度方程式编写如下 rrrvel.m 函数： function y=rrrvel(x) % %Input parameters % %x(1)=theta-1 %x(2)=theta-2 %x(3)=theta-3 %x(4)=theta-4 %x(5)=dtheta-1 %x(6)=l1 %x(7)=l2
=
-L1cos(θ1+180°)
L’1 sinθ1
-L’1cosθ1
Vcx= -L1w1sin(θ1+180°)-w2 L2sinθ2 Vcy=Lw2cos(θ1+180°)+ w2 L2cosθ2
（3）加速度方程
-L2sinθ2 L3sinθ3
- L4sinθ4
0
a2
L2cosθ2 -L3cosθ3
q=
1.0e+003 *
Columns 1 through 8
-0.0064 -0.0062 -0.0061 0.0085 0.0089 0.0092 0.0235 0.0237 0.0239 -1.0578 -1.0897 -1.1226
-0.0061 0.0096 0.0241 -1.1568
-0.0060 0.0101 0.0244 -1.1926