How can I complete this code or reach the result according to this mechanism?
1 view (last 30 days)
Show older comments
% Define physical parameters of the mechanism
AB = 0.09; % m, length of arm 1
BC = 0.4; % m, length of arm 2
CD = 0.7; % m, length of arm 3
La = 0.75; % m
Lb = 0.35; % m
% Define simulation parameters
Nmax = 100; % Maximum number of iterations
x = [20*pi/180, 320*pi/180]; % Adjusted initial guess for th2, th3 (in radians)
xe = 0.0001 * abs(x); % Slightly relaxed error tolerance
dth = 5 * pi / 180; % Step size for thq
thQ = 0:dth:2*pi; % Range of thq (in radians)
w1 = -0.2944 * ones(size(thQ)); % Constant angular velocity w1
acc1 = zeros(size(thQ)); % Constant angular acceleration acc1
% Preallocate arrays for results
th2 = zeros(size(thQ));
th3 = zeros(size(thQ));
w2 = zeros(size(thQ));
w3 = zeros(size(thQ));
acc2 = zeros(size(thQ));
acc3 = zeros(size(thQ));
% Main loop for each th2 value
for k = 1:length(thQ)
x_temp = x; % Temporary variable for the current iteration
kerr = 1; % Flag for convergence
% Iterative solver for position analysis
for n = 1:Nmax
th2_temp = x_temp(1);
th3_temp = x_temp(2);
J = [-BC*sin(th2_temp), -CD*sin(th3_temp); BC*cos(th2_temp), CD*cos(th3_temp)];
% Check for singularity in Jacobian
if cond(J) > 1e10
disp(['Warning: Jacobian is near singular at th2 index: ', num2str(k)]);
break;
end
f = [-(AB*cos(thQ(k)) + BC*cos(th2_temp) + CD*cos(th3_temp) - La); ...
-(AB*sin(thQ(k)) + BC*sin(th2_temp) + CD*sin(th3_temp) + Lb)];
eps = J\f; % More efficient than inv(J)*f
x_temp = x_temp + eps';
if all(abs(eps) < xe)
kerr = 0;
break;
end
end
if kerr == 1
disp(['Error: Solution did not converge at th2 index: ', num2str(k), ' (th2 = ', num2str(thQ(k) * 180/pi), ' degrees)']);
else
th2(k) = x_temp(1);
th3(k) = x_temp(2);
% Velocity analysis
fv = [-AB*sin(thQ(k))*w1(k); AB*cos(thQ(k))*w1(k)];
vel = J\fv;
w2(k) = vel(1);
w3(k) = vel(2);
% Acceleration analysis
fa = [-AB*cos(thQ(k))((w1(k))^2) - BC*cos(th2(k))((w2(k))^2) - CD*cos(th3(k))*((w3(k))^2) - AB*sin(thQ(k))*acc1(k); ...
AB*cos(thQ(k))acc1(k) - BC*sin(th2(k))((w2(k))^2) - CD*sin(th3(k))((w3(k))^2) - AB*sin(thQ(k))((w1(k))^2)];
acc = J\fa;
acc2(k) = acc(1);
acc3(k) = acc(2);
end
end
% Convert angles to degrees for plotting
th2d = th2 * 180/pi;
th3d = th3 * 180/pi;
% Plotting results
figure(1);
subplot(3,3,1), plot(thQ, th2d, 'r', 'LineWidth', 2), xlabel('\theta_Q (rad)'), ylabel('\theta_2 (deg)'), grid on;
subplot(3,3,2), plot(thQ, w2, 'r', 'LineWidth', 2), xlabel('\theta_Q (rad)'), ylabel('\omega_2 (rad/s)'), grid on;
subplot(3,3,3), plot(thQ, acc2, 'r', 'LineWidth', 2), xlabel('\theta_Q (rad)'), ylabel('\alpha_2 (rad/s^2)'), grid on;
subplot(3,3,4), plot(thQ, th3d, 'r', 'LineWidth', 2), xlabel('\theta_Q (rad)'), ylabel('\theta_3 (deg)'), grid on;
subplot(3,3,5), plot(thQ, w3, 'r', 'LineWidth', 2), xlabel('\theta_Q (rad)'), ylabel('\omega_3 (rad/s)'), grid on;
subplot(3,3,6), plot(thQ, acc3, 'r', 'LineWidth', 2), xlabel('\theta_Q (rad)'), ylabel('\alpha_3 (rad/s^2)'), grid on;
2 Comments
Steven Lord
on 23 Jan 2024
It's not clear to me what help you're hoping to receive. Please ask a specific question about where you're having difficulty and we may be able to provide some guidance.
Answers (0)
See Also
Categories
Find more on Spectral Measurements in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!