Subs Function not working properly

Hello,

In the code below, I determine and solve the equations of motion of an upright slider crank mechanism by setting up force and moment balance equations, solving for the unkown forces, and then substituting them back into one of the moment equations to arrive at the equations of motion. I can sucessfully plot the angular displacement, velocity, and acceleration of the crank and coupler linkage, but I want to also be able to plot the forces acting on the system and slider kinematics as well. It seemed to be pretty straightforward; I thought I would have just had to convert the solution struct equations for the forces to function handles, but I just cannot seem to work it out, especially because of the fact that the equations contain the symbolic functions theta(t), phi(t), diff(theta,t), diff(phi,t), diff(theta,t,t), and diff(phi,t,t). These are all plotted and known after the solver finishes the numerical analysis (as different variables, of course), but I can't seem to figure out a way to substitute them into any of the force solution equations. When I try substituting new variables for diff(theta,t), diff(phi,t), and so on, the subs function doesn't work properly:

If someone could help point me in the right direction, it would be greatly appreciated, thanks.

% Author: Jason Louison

clearvars

syms theta(t) R_C_m R_C m_C I_C 

syms phi(t) R_L_m R_L m_L I_L

syms m_S D

syms T_A T_B F_A_x F_A_y F_B_x F_B_y F_C_x F_N F_C_y F_P_x F_P_y

syms g u 

%Crankshaft Linkage

%x and y coordinates of COG (Referenced from point A)

R_C_m_x = -R_C_m*sin(theta);
R_C_m_y = R_C_m*cos(theta);

% x and y coordinates of end of link (Referenced from point A)

R_C_x = -R_C*sin(theta);
R_C_y = R_C*cos(theta);

% We can differentiate the COG coordinate equations to obtain the instantaneous 
% acceleration of the COG: 

% x and y components of COG acceleration

a_1_x = diff(R_C_m_x,2);
a_1_y = diff(R_C_m_y,2);

%Connecting Rod

R_L_m_x = -R_L_m*sin(phi);
R_L_m_y = R_L_m*cos(phi);

R_L_x = -R_L*sin(phi);
R_L_y = R_L*cos(phi);

%Connecting Rod COG Position Acceleration: 

P_COG_L_x = R_C_x + R_L_m_x;
P_COG_L_y = R_C_y + R_L_m_y; 

a_2_x = diff(P_COG_L_x,2);
a_2_y = diff(P_COG_L_y,2);

%Piston COG Position and Acceleration 

RSx = R_C_x+R_L_x;
RSy = R_C_y+R_L_y;

a_3_x = diff(RSx,2); 
a_3_y = diff(RSy,2); 

%Crankshaft Force/Moment Balance

F1x = F_A_x + F_B_x == m_C*a_1_x;
F1y = F_A_y + F_B_y + m_C*g == m_C*a_1_y;
M1 = - (R_C_m_x*F_A_y - R_C_m_y*F_A_x) + (R_C_x - R_C_m_x)*F_B_y - (R_C_y - R_C_m_y)*F_B_x + T_A == I_C*diff(theta,t,t);

%Connecting Rod Force/Moment Balance

F2x = -(F_B_x + F_C_x) == m_L*a_2_x; 
F2y = -(F_B_y + F_C_y) + m_L*g == m_L*a_2_y; 
M2 = -(-(R_L_m_x*F_B_y - R_L_m_y*F_B_x) + (R_L_x-R_L_m_x)*F_C_y - (R_L_y-R_L_m_y)*F_C_x) == I_L*diff(phi,t,t);

%Piston Force Balance 

F3x = F_C_x + F_N == 0;

F3y = F_C_y + m_S*g - u*F_N == m_S*a_3_y;

%Solving for x and y Forces at Nodes A, B, and C

eqns_set_1 = [F1x; F1y; F2x; F2y; F3x; F3y; M1; M2]; %Force and Moment 'Coefficient' Vector.

vars_set_1 = [F_A_x; F_A_y; F_B_x; F_B_y; F_C_x; F_C_y; F_N; T_A]; %Unknown Vector

nodeForces = solve(eqns_set_1,vars_set_1); %Force Solution Struct Matrix 

%derive the equation of motion for the connecting rod and crankshaft. 

M_C = subs(M1,[F_A_x, F_A_y, F_B_x, F_B_y],[nodeForces.F_A_x, nodeForces.F_A_y, nodeForces.F_B_x, nodeForces.F_B_y]);

syms S(t) % <-- Piston/Slider Displacement

Dx = -R_C*sin(theta) -R_L*sin(phi) - D == 0;
Dy = R_C*cos(theta) + R_L*cos(phi) - S == 0;

Vx = diff(Dx,1);
Vy = diff(Dy,1);

Ax = diff(Dx,2);
Ay = diff(Dy,2);

%Derive Connecting Rod Angular Displacement, Velocity, and Acceleration

phif = isolate(Dx,phi);
dphif = isolate(Vx,diff(phi,t));
ddphif = isolate(Ax,diff(phi,t,t));

%Derive Slider Displacement, Velocity, and Acceleration

Sf = isolate(Dy,S);
dSf = isolate(Vy,diff(S,t));
ddSf = isolate(Ay,diff(S,t,t));

syms Y1 Y2 Y3 Y4 Y5 Y6

Disp = subs(rhs(Sf),[theta phi],[Y1 Y2])

Disp = 

Veloce = subs(rhs(Sf),[theta phi diff(theta,1) diff(phi,1)],[Y1 Y2 Y3 Y4])

Veloce = 

%==================== Input Parameters =====================%

%Crankshaft: COG Vector, Length Vector, Mass, Inertia

R_C_m = 37.5e-3; R_C = 50e-3; m_C = 3; I_C = 0.75;

%Crankshaft: COG Vector, Length Vector, Mass, Inertia

R_L_m = 62.5e-3; R_L = 150e-3; m_L = 0.650; I_L = 0.5;

%Piston and Cylinder: Piston Mass, Bore Centerline Offset

m_S = 0.7; D = 0;

%Forces: External Torque Applied to Crankshaft 

T_A = 0; 

%Constants: Gravitational Acceleration; Cylinder Wall Friction Coefficient 

g = -9.81; u = 0; 

A = subs([M_C ddphif]);

V = odeToVectorField(A);

%Convert the differential equations to a function

M = matlabFunction(V,'vars',{'t','Y'});

%initial Conditions 

theta_0 = 0; %Initial Crank Angle (rad)

dtheta_0 = 100; %Initial Crank Angular Velocity (rad/s)

phi_0 = -asin((D+R_C*sin(theta_0))/R_L);

dphi_0 = -(R_C*cos(theta_0)*dtheta_0)/(R_L*cos(phi_0));

%Integration Function

[t,Y] =  rk4th(M,0,0.12,[theta_0 dtheta_0 phi_0 dphi_0],1000000);

[R,C] = size(Y);

%Plot State Variables

figure('Name','Slider Crank Time Response')
for i = 1:C
subplot(3,2,i)
plot(t,Y(:,i),'Color','b')
if i == 1 
    title('Crank Anglular Displacement vs. Time')
    xlabel('t (sec)')
    ylabel('\theta (rad)')
end
if i == 2 
    hold on 
    title('Coupler Anglular Displacement vs. Time')
    xlabel('t (sec)')
    ylabel('\phi (rad)')
end 
if i == 3 
    hold on 
    title('Crank Anglular Velocity vs. Time')
    xlabel('t (sec)')
    ylabel('\theta'' (rad/s)')
end
if i == 4 
    hold on 
    title('Coupler Anglular Velocity vs. Time')
    xlabel('t (sec)')
    ylabel('\phi'' (rad/s)')
end 
if i == 5 
    hold on 
    title('Crank Anglular Acceleration vs. Time')
    xlabel('t (sec)')
    ylabel('\theta'''' (rad/s^2)')
end 
if i == 6 
    hold on 
    title('Coupler Anglular Acceleration vs. Time')
    xlabel('t (sec)')
    ylabel('\phi'''' (rad/s^2)')
end
end

%Runge Kutta 4th Order Integration Function. 

function [t,Y] =  rk4th(M,t0,tf,y0,N)
h = (tf-t0)/N;
t = t0:h:tf;
y = zeros(length(t),4);
dy = zeros(length(t),4);
y(1,[1 2 3 4]) = y0;
for i = 1:(length(t)-1)
    k_1 = M(t(i),y(i,:)).';
    k_2 = M(t(i),y(i,:)+0.5*h.*k_1).';
    k_3 = M(t(i),y(i,:)+0.5*h.*k_2).';
    k_4 = M(t(i),y(i,:)+h.*k_3).';
    y(i+1,:) = y(i,:) + (h/6).*(k_1+2.*k_2+2.*k_3+k_4);
    dy(i,:) = M(t(i),y(i,:)).';
    dy(i+1,:) = M(t(i+1),y(i+1,:)).';
end
y1 = y(:,1);
y2 = y(:,2);
y3 = y(:,3);
y4 = y(:,4);
y5 = dy(:,2);
y6 = dy(:,4);
Y = [y1 y3 y2 y4 y5 y6];
end

Answer 1

[Walter Roberson]

syms theta(t) R_C_m R_C m_C I_C 
syms phi(t) R_L_m R_L m_L I_L

theta and phi are symfun, not just plain sym

Veloce = subs(rhs(Sf),[theta phi diff(theta,1) diff(phi,1)],[Y1 Y2 Y3 Y4])

diff() of a symfun is a symfun .

So you are [] four symfun. The result is a single symfun, not a vector of four symfun.

You need to use

Veloce = subs(rhs(Sf), {theta phi diff(theta,1) diff(phi,1)}, {Y1 Y2 Y3 Y4})

[Jason Louison]

Hello Walter,

I also had made an error in the subs equation for Veloce. I was trying to sub four variables into an equation of two variables. However, it seems to be there is a bigger, deeper problem here, as the substitution works completely fine for the displacement equation. Why not the diff() symfunctions as well? I believe I have successfully done this operation before without running into any issues. To make matters worse, when I corrected my error in the Veloce subs equation, I get this result:

This is, of course, a nonsensical result. I know how to plot these functions, its a matter of using 'matlabFunction' to convert the symbolic functions to function handles, and then evaluating at the required values. But I cannot do so without substituting the symbolic function variables (theta(t), diff(theta,t), etc.) with a regular sym variable (Y,X,Z1, etc). I have noticed, however, that Sf, dSf, and ddSf are all "sym" variables, when I expected them to be "symfun" variables. This is the case with the force solution equations as well. I don't know if that has anything to do with it though, because it seems to work the same way.

[Jason Louison]

Walter,

I managed to figure it out somehow. In order for it to work, you need to enter the "diff()" variables arranged from highest order to lowest order, like so:

F(i) = subs(Fsol_v(i,1),[diff(phi,2) diff(theta,2) diff(phi,1) diff(theta,1) phi theta],{Y6 Y5 Y4 Y3 Y2 Y1});

I can successfully plot anything I need to now. Thank you for pointing me in the right direction. Now, if only there were an easy way to combine my 3x2 subplot figures into a single 3x4 subplot figure. One thing after the other, haha. Results from a run are shown below.

Thank you for pointing me in the right direction, I am truly grateful.

Code:

% Author: Jason Louison

clearvars

syms theta(t) R_C_m R_C m_C I_C 

syms phi(t) R_L_m R_L m_L I_L

syms m_S D

syms T_A T_B F_A_x F_A_y F_B_x F_B_y F_C_x F_N F_C_y F_P_x F_P_y

syms g u 

%Crankshaft Linkage

%x and y coordinates of COG (Referenced from point A)

R_C_m_x = -R_C_m*sin(theta);
R_C_m_y = R_C_m*cos(theta);

% x and y coordinates of end of link (Referenced from point A)

R_C_x = -R_C*sin(theta);
R_C_y = R_C*cos(theta);

% We can differentiate the COG coordinate equations to obtain the instantaneous 
% acceleration of the COG: 

% x and y components of COG acceleration

a_1_x = diff(R_C_m_x,2);
a_1_y = diff(R_C_m_y,2);

%Connecting Rod

R_L_m_x = -R_L_m*sin(phi);
R_L_m_y = R_L_m*cos(phi);

R_L_x = -R_L*sin(phi);
R_L_y = R_L*cos(phi);

%Connecting Rod COG Position Acceleration: 

P_COG_L_x = R_C_x + R_L_m_x;
P_COG_L_y = R_C_y + R_L_m_y; 

a_2_x = diff(P_COG_L_x,2);
a_2_y = diff(P_COG_L_y,2);

%Piston COG Position and Acceleration 

RSx = R_C_x+R_L_x;
RSy = R_C_y+R_L_y;

a_3_x = diff(RSx,2); 
a_3_y = diff(RSy,2); 

%Crankshaft Force/Moment Balance

F1x = F_A_x + F_B_x == m_C*a_1_x;
F1y = F_A_y + F_B_y + m_C*g == m_C*a_1_y;
M1 = - (R_C_m_x*F_A_y - R_C_m_y*F_A_x) + (R_C_x - R_C_m_x)*F_B_y - (R_C_y - R_C_m_y)*F_B_x + T_A == I_C*diff(theta,t,t);

%Connecting Rod Force/Moment Balance

F2x = -(F_B_x + F_C_x) == m_L*a_2_x; 
F2y = -(F_B_y + F_C_y) + m_L*g == m_L*a_2_y; 
M2 = -(-(R_L_m_x*F_B_y - R_L_m_y*F_B_x) + (R_L_x-R_L_m_x)*F_C_y - (R_L_y-R_L_m_y)*F_C_x) == I_L*diff(phi,t,t);

%Piston Force Balance 

F3x = F_C_x + F_N == 0;

F3y = F_C_y + m_S*g - u*F_N == m_S*a_3_y;

%Solving for Forces at Nodes A, B, and C

K = [F1x; F1y; F2x; F2y; F3x; F3y; M1; M2];

U = [F_A_x; F_A_y; F_B_x; F_B_y; F_C_x; F_C_y; F_N; T_A];

[G,F] = equationsToMatrix(K,U);

Fsol = G\F;

%Now that we have equations for all the dynamic forces acting on the
%mechanism, we can substitute them into the moment balance equations to
%derive the equation of motion for the connecting rod and crankshaft. 

M_C = subs(M1,[F_A_x, F_A_y, F_B_x, F_B_y],[Fsol(1,1), Fsol(2,1), Fsol(3,1), Fsol(4,1)]);

syms S(t)

Dx = -R_C*sin(theta) -R_L*sin(phi) - D == 0;
Dy = R_C*cos(theta) + R_L*cos(phi) - S == 0;

Vx = diff(Dx,1);
Vy = diff(Dy,1);

Ax = diff(Dx,2);
Ay = diff(Dy,2);


phif = isolate(Dx,phi);
dphif = isolate(Vx,diff(phi,t));
ddphif = isolate(Ax,diff(phi,t,t));

Sf = isolate(Dy,S);
dSf = isolate(Vy,diff(S,t));
ddSf = isolate(Ay,diff(S,t,t));

%==================== Input Parameters: =====================%

%Crankshaft

R_C_m = 0; R_C = 50e-3; m_C = 3; I_C = 0.75;

%Connecting Rod

R_L_m = 75e-3; R_L = 150e-3; m_L = 0.650; I_L = 0.5;

%Piston and Cylinder 

m_S = 0.7; D = 0; %D = Bore Centerline Offset

%Forces 

T_A = 0; 

%Constants 

g = -0; u = 0; %Gravitational Acceleration; Cylinder Wall Friction Coefficient

A = subs([M_C ddphif]);

Fsol_v = subs(Fsol);

V = odeToVectorField(A);

M = matlabFunction(V,'vars',{'t','Y'});

syms Y1 Y2 Y3 Y4 Y5 Y6

idx = size(Fsol_v);

for i = 1:idx
F(i) = subs(Fsol_v(i,1),[diff(phi,2) diff(theta,2) diff(phi,1) diff(theta,1) phi theta],{Y6 Y5 Y4 Y3 Y2 Y1});
Ff{i} = matlabFunction(F(i,1),'vars',{t,[Y1 Y2 Y3 Y4 Y5 Y6]});
end

%initial Conditions 

theta_0 = 0; %Initial Crank Angle (rad)

dtheta_0 = 105; %Initial Crank Angular Velocity (rad/s)

phi_0 = -asin((D+R_C*sin(theta_0))/R_L);

dphi_0 = -(R_C*cos(theta_0)*dtheta_0)/(R_L*cos(phi_0));

[t,Y] =  rk4th(M,0,0.06,[theta_0 dtheta_0 phi_0 dphi_0],1000000);

[R,C] = size(Y);

for i = 1:idx
fresult{i} = Ff{i}(t,Y(:,[1 2 3 4 5 6]));
end

figure('Name','Slider Crank Time Response')
for i = 1:C
subplot(3,2,i)
plot(t,Y(:,i),'Color','b')
if i == 1 
    title('Crank Anglular Displacement vs. Time')
    xlabel('Time (sec)')
    ylabel('Anglular Displacement \theta (rad)')
end
if i == 2 
    hold on 
    title('Coupler Anglular Displacement vs. Time')
    xlabel('Time (sec)')
    ylabel('Anglular Displacement \phi (rad)')
end 
if i == 3 
    hold on 
    title('Crank Anglular Velocity vs. Time')
    xlabel('Time (sec)')
    ylabel('Anglular Velocity \theta'' (rad/s)')
end
if i == 4 
    hold on 
    title('Coupler Anglular Velocity vs. Time')
    xlabel('Time (sec)')
    ylabel('Anglular Velocity \phi'' (rad/s)')
end 
if i == 5 
    hold on 
    title('Crank Anglular Acceleration vs. Time')
    xlabel('Time (sec)')
    ylabel('Anglular Acceleration \theta'''' (rad/s^2)')
end 
if i == 6 
    hold on 
    title('Coupler Anglular Acceleration vs. Time')
    xlabel('Time (sec)')
    ylabel('Anglular Acceleration \phi'''' (rad/s^2)')
end
end

figure('Name','Forces Acting on Mechanism')
for i = 1:C
subplot(3,2,i)
plot(t,fresult{1,i},'Color','b')
if i == 1 
    title('x-Direction Force at Crank Journal vs. Time')
    xlabel('Time (sec)')
    ylabel('Force (N)')
end
if i == 2 
    hold on 
    title('y-Direction Force at Crank Journal vs. Time')
    xlabel('Time (sec)')
    ylabel('Force (N)')
end 
if i == 3 
    hold on 
    title('x-Direction Force at Rod Journal vs. Time')
    xlabel('Time (sec)')
    ylabel('Force (N)')
end
if i == 4 
    hold on 
    title('y-Direction Force at Rod Journal vs. Time')
    xlabel('Time (sec)')
    ylabel('Force (N)')
end 
if i == 5 
    hold on 
    title('x-Direction Force at Piston vs. Time')
    xlabel('Time (sec)')
    ylabel('Force (N)')
end 
if i == 6 
    hold on 
    title('y-Direction Force at Piston vs. Time')
    xlabel('Time (sec)')
    ylabel('Force (N)')
end
end
%Runge Kutta 4th Order Integration Function. 

function [t,Y] =  rk4th(M,t0,tf,y0,N)
h = (tf-t0)/N;
t = t0:h:tf;
y = zeros(length(t),4);
dy = zeros(length(t),4);
y(1,[1 2 3 4]) = y0;
for i = 1:(length(t)-1)
    k_1 = M(t(i),y(i,:)).';
    k_2 = M(t(i),y(i,:)+0.5*h.*k_1).';
    k_3 = M(t(i),y(i,:)+0.5*h.*k_2).';
    k_4 = M(t(i),y(i,:)+h.*k_3).';
    y(i+1,:) = y(i,:) + (h/6).*(k_1+2.*k_2+2.*k_3+k_4);
    dy(i,:) = M(t(i),y(i,:)).';
    dy(i+1,:) = M(t(i+1),y(i+1,:)).';
end
y1 = y(:,1);
y2 = y(:,2);
y3 = y(:,3);
y4 = y(:,4);
y5 = dy(:,2);
y6 = dy(:,4);
Y = [y1 y3 y2 y4 y5 y6];
end