% impact of a compound pendulum with a flat surface clear; clc; close; syms t real global L yc theta0 deltar g=9.81; L=1; R=0.01; rho=7800; m=pi*R^2*L*rho; theta = sym('theta(t)'); omega = [0 0 diff(theta,t)]; alpha = diff(omega,t); c = cos(theta); s = sin(theta); xC = (L/2)*c; yC = (L/2)*s; rC = [xC yC 0]; xA = L*c; yA = L*s; rA = [xA yA 0]; vA = diff(rA,t); G = [0 m*g 0]; IC = m*L^2/12; IO = IC + m*(L/2)^2; MO = cross(rC,G); % free fall equation eq = MO/IO; eqz = eq(3); qt = {diff(theta,t),theta}; ut = {'x(2)', 'x(1)'}; eq0 = subs(eqz, qt, ut); % first differential equation dx1 = sym('x(2)'); % second differential equation dx2 = eq0; % define right-hand side vector dx1dt = char(dx1); dx2dt = char(dx2); fid = fopen('CP.m','w+'); fprintf(fid,'function dx = CP(t,x)\n'); fprintf(fid,'dx = zeros(2,1);\n'); fprintf(fid,'dx(1) = '); fprintf(fid,dx1dt); fprintf(fid,';\n'); fprintf(fid,'dx(2) = '); fprintf(fid,dx2dt); fprintf(fid,';'); fclose(fid); option0 = odeset('RelTol',1e-3,'MaxStep',1e-3,'Events', @event0); [tt0, xs0, t0, ye0] = ode45(@CP, [0 1], [0 0], option0); theta0 = ye0(1); omega0 = ye0(2); % fprintf('t0= %g [s] \n',t0); fprintf('Initial conditions for impact\n'); fprintf('theta0= %g [rad] = %g [degrees] \n', theta0, theta0*180/pi); fprintf('omega0= %g [rad/s] \n', omega0); list0 = {omega0, theta0} ; vAn = subs(vA, qt, list0); % double(S) converts the symbolic object S to a numeric object v0 = double(vAn); % converts the symbolic object vBn to a numeric object fprintf('v0 = [ %g, %g, %g ] [m/s]\n', v0); % elastic compression E1=200e9; Sy=1.12e9; nu=.33; Ep=(2*(1-nu^2)/E1)^-1; C =1.295*exp(.736*nu); yc=(((pi*C*Sy)/(2*Ep))^2)*R; delta= L*sin(theta)-L*sin(theta0); k1=(2/(3*(1-nu^2))*E1*sqrt(R)); Pe=[0,-k1*delta^1.5,0]; eqe = (MO + cross(rA, Pe))/IO; eqze = eqe(3); % second differential equation dx2e = subs(eqze, qt, ut); % define right-hand side vector dx2dte = char(dx2e); fid = fopen('CPe.m','w+'); fprintf(fid,'function dx = CPe(t,x)\n'); fprintf(fid,'dx = zeros(2,1);\n'); fprintf(fid,'dx(1) = '); fprintf(fid,dx1dt); fprintf(fid,';\n'); fprintf(fid,'dx(2) = '); fprintf(fid,dx2dte); fprintf(fid,';'); fclose(fid); optione= odeset('RelTol',1e-6,'InitialStep',1e-10,'MaxStep',1e-8,'Events',@evente); [tte, xse, te, yee]=ode45(@CPe, [0 0.1], [theta0 omega0], optione); thetae = real(yee(1)); omegae = real(yee(2)); fprintf(' \n'); fprintf('elastic compression \n'); fprintf('te = %g [s] \n',te); fprintf('thetae = %g [rad] = %g [degrees] \n',thetae, thetae*180/pi); fprintf('omegae = %g [rad/s] \n', omegae); liste = {omegae, thetae}; deltae = double(subs(delta, qt, liste)); fprintf('deltae = %g [m] \n', deltae); pe = eval(subs(Pe(2), qt, liste)); fprintf('Pe = %g [N] \n', pe); fprintf(' \n'); thetase = real(xse(:,1)); omegase = real(xse(:,2)); % figure(1) % subplot(2,1,1),plot(tte,thetase*180/pi,'r'),... % xlabel('t (s)'),ylabel('\theta (deg)'),grid,... % subplot(2,1,2),plot(tte,omegase,'b'),... % xlabel('t (s)'),ylabel('\omega (rad/s)'),grid; % elasto-plastic compression ey = Sy/Ep; B = 0.14*exp(23*ey); A = sqrt(R*delta*(delta/(1.9*yc))^B); HS = 2.84-.92*(1-cos(pi*A / R)); Pc = (4/3)*(R/Ep)^2*(pi*C*Sy/2)^3; y = (delta/yc); Pep = Pc*(exp((-1/4)*y^(5/12))*y^(3/2)+4/C*HS *(1-exp((-1/25)*y^(5/9)))*y); fprintf(' \n'); pepe=subs(vpa(-Pep,9), theta, thetae); fprintf('Pep(te+0) = %g [N] \n', pepe); fprintf(' \n'); eqep = ( MO + cross(rA, [0 -Pep 0]) )/IO; eqzep = eqep(3); % second differential equation dx2ep = subs(eqzep, qt, ut); % define right-hand side vector dx2dtep = char(dx2ep); fid = fopen('CPep.m','w+'); fprintf(fid,'function dx = CPep(t,x)\n'); fprintf(fid,'dx = zeros(2,1);\n'); fprintf(fid,'dx(1) = '); fprintf(fid,dx1dt); fprintf(fid,';\n'); fprintf(fid,'dx(2) = '); fprintf(fid,dx2dtep); fprintf(fid,';'); fclose(fid); optionep=odeset('RelTol',1e-6,'InitialStep',1e-10,'MaxStep',1e-8,'Events',@eventep); [ttm, xsm, tm, yem]=ode45(@CPep, [te 0.1], [thetae omegae], optionep); thetam = yem(1); omegam = yem(2); fprintf(' \n'); fprintf('elasto-plastic compression \n'); fprintf('tm = %g [s] \n',tm); fprintf('thetam = %g [rad] = %g [degrees] \n',thetam, thetam*180/pi); fprintf('omegam = %g [rad/s] \n', omegam); listm = {omegam, thetam}; deltam = double(subs(delta, qt, listm)); fprintf('deltam = %g [m] \n', deltam); Pm=subs(vpa(-Pep,9), theta, thetam); fprintf('Pm = %g [N] \n', Pm); fprintf(' \n'); thetasm = xsm(:,1); omegasm = xsm(:,2); % figure(2) % subplot(2,1,1),plot(ttm,thetasm*180/pi,'r'),... % xlabel('t (s)'),ylabel('\theta (deg)'),grid,... % subplot(2,1,2),plot(ttm,omegasm,'b'),... % xlabel('t (s)'),ylabel('\omega (rad/s)'),grid % restitution deltams = deltam/yc; deltar = deltam*(1.02*(1-((deltams+5.9)/6.9)^-0.54)); Rr = 1/(deltam-deltar)^3 * (3/4 * Pm/Ep)^2; k1r = 2/(3*(1-nu^2))*E1*sqrt(Rr); Pr = [0, -k1r*(delta-deltar)^1.5, 0]; fprintf(' \n'); pr = eval(subs(Pr(2), qt, listm)); fprintf('Pm(tm+0) = %g [N] \n', pr); fprintf(' \n'); eqer = ( MO + cross(rA, Pr) )/IO; eqzr = eqer(3); % second differential equation dx2r = subs(eqzr, qt, ut); % define right-hand side vector dx2dtr = char(dx2r); fid = fopen('CPr.m','w+'); fprintf(fid,'function dx = CPr(t,x)\n'); fprintf(fid,'dx = zeros(2,1);\n'); fprintf(fid,'dx(1) = '); fprintf(fid,dx1dt); fprintf(fid,';\n'); fprintf(fid,'dx(2) = '); fprintf(fid,dx2dtr); fprintf(fid,';'); fclose(fid); optionr= odeset('RelTol',1e-6,'InitialStep',1e-10,'MaxStep',1e-8,'Events',@eventr); [ttr, xsr, tf, yer]=ode45(@CPr, [tm 0.1], [thetam 0], optionr); thetaf = real(yer(1)); omegaf = real(yer(2)); fprintf(' \n'); fprintf('restitution \n'); fprintf('tf = %g [s] \n', tf); fprintf('thetaf = %g [rad] = %g [degrees] \n',thetaf, thetaf*180/pi); fprintf('omegaf = %g [rad/s] \n', omegaf); listf = {omegaf, thetaf}; deltaf = double(subs(delta, qt, listf)); fprintf('deltaf = %g [m] \n', deltaf); fprintf(' \n'); thetasr = real(xsr(:,1)); omegasr = real(xsr(:,2)); % figure(3) % subplot(2,1,1),plot(ttr,thetasr*180/pi,'r'),... % xlabel('t (s)'),ylabel('\theta (deg)'),grid,... % subplot(2,1,2),plot(ttr,omegasr,'b'),... % xlabel('t (s)'),ylabel('\omega (rad/s)'),grid; fprintf('e = - omegaf/omega0 = %g \n', -omegaf/omega0); figure(4) plot(tte,thetase*180/pi,'r'), hold on,... plot(ttm,thetasm*180/pi,'b'), hold on,... plot(ttr,thetasr*180/pi,'g'), hold off,... xlabel('t (s)'), ylabel(' \theta (deg)'), grid; figure(5) plot(tte,omegase,'r'), hold on,... plot(ttm,omegasm,'b'), hold on,... plot(ttr,omegasr,'g'), hold off,... xlabel('t (s)'), ylabel(' \omega (rad/s)'), grid;