% Push Rod Simulation
% ME 244L Dynamic Systems and Controls
% R.G. Longoria: revised 1/26/98
% Solution for Problem 7-6 from
% "Physical System Dynamics" by Rosenberg and Karnopp (1983)
%
% This routine uses a fixed-step RK4 algorithm (rk4fixed)
% (Alternatively, ode23 or ode45 can also be used; see below.)
%
% Note: this routine was updated for Matlab 5.0.1, July 1998 by rgl.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Define global parameters
%
global k2 m34 k6 m5 b7 omega
%
% stiffness of pushrod, k2 N/m
k2 = 100000;
% rocker modulus
m34 = 1;
% valve spring stiffness, k6 N/m
k6 = 1000;
% valve mass, m5 kg
m5 = 0.1;
% valve guide/lifter damping, b7 N*s/m
b7 = 100;
%
% frequency of cam rotation
% 240 rad/sec is about 2292 rpm
% 110 rad/sec is about 1050 rpm
% 10 rad/sec is about 95 rpm
omega = 240;
%
%
% 
% Note: in matlab 5.0 and higher, move global declaration
% to the beginning - before you define the variables.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define the initial conditions for the simulation.
%
% Pre-loads in both pushrod and valve spring of 1 N
% Initial velocity of valve is zero, so p5(0) = 0
%
preload = 1;
%
% Initial condition vector:
%
x0=[preload/k2;preload/k6;0];
%
% Note: the cam input velocity is defined in the function file
% pushrod.m as: u=2.4*sin(omega*t), where omega = 240 rad/sec
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Time integration parameters
%
% time interval for fixed-step simulation: dt = 0.001 sec
%
dt=0.001;
%
% start time
t0 = 0.0;
% go to final time of tf
tf = 0.08;
%
N = floor((tf-t0)/dt);
%
% Call rk4fixed.m
%
[t,x]=rk4fixed('pushrod',[t0 tf],x0,N);
% 
% If you want to use one of the solvers from the Matlab ODE suite:
%
%[t,x]=ode45('pushrod',[t0 tf],x0);
%
% the states are returned in the vector x = [x(1);x(2);x(3)]
% x(1) = axial deflection of pushrod
% x(2) = deflection of valve spring
% x(3) = momentum of valve mass
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% We are interested in the output variables:
%
% y1 = x6 = x(2), valve spring deflection
% y2 = F2 = k2*x2 = k2*x(1), pushrod force
%
% The cam input displacement is the integral of the velocity
% which is
xcam = 2.4*(1 - cos(omega*t))/omega;
%
% We'll need xcam to compare to the simulated response
%
% Plotting:
%
subplot(2,1,1), plot(t,1000*(x(:,2)-x(1,2)),t,1000*xcam,'-.')
title('Plot Contrasting x6(t) to Cam Motion')
legend('x6','cam')
ylabel('Displacement, mm')
xlabel('Time, sec')
subplot(2,1,2), plot(t,x(:,1))
title('Plot of contact force F(t) = k2*x2(t)')
ylabel('Force, N')
xlabel('Time, sec')
