(* Problem 4.8 *)(* VELOCITY and ACCELERATION ANALYSIS - input angle phi *)Apply[Clear,Names["Global`*"]];Off[General::spell];Off[General::spell1];(* Input data *)AB = .120 ;AC = .060 ;BD = .240 ;DE = .330 ;EF = .190 ;La = .300 ;Lb = .070 ;n = 100 ; (* rpm *)omega = n*N[Pi]/30 ; (* rad/s *)initdata = {f[t]->N[Pi]150/180, f'[t]->omega, f''[t]->0};(* Position of joint A *)xA = yA = 0;(* Position of joint C *)xC = -AC ;yC = 0 ;(* Position of joint F *)xF = La ;yF = Lb ;(* Position of joint B *)XB = AB Cos[f[t]] ;YB = AB Sin[f[t]] ;xBs=XB/.initdata;yBs=YB/.initdata;Print["xB[t] = ", XB," = ", xBs, " m" ];Print["yB[t] = ", YB," = ", yBs, " m" ];(* Linear velocity of joint B *)vBx = D[XB,t];vBy = D[YB,t];vBxs = vBx/.initdata;vBys = vBy/.initdata;Print["vBx = ", vBx ," = ", vBxs, " m/s" ];Print["vBy = ", vBy ," = ", vBys, " m/s" ];(* Linear acceleration of joint B *)aBx = D[vBx,t];aBy = D[vBy,t];aBxs = aBx/.initdata;aBys = aBy/.initdata;Print["aBx = ", aBx ," = ", aBxs, " m/s^2" ];Print["aBy = ", aBy ," = ", aBys, " m/s^2" ];ruleB={xB[t]->xBs,yB[t]->yBs, xB'[t]->vBxs,yB'[t]->vBys, xB''[t]ØaBxs,yB''[t]ØaBys}; (* Angular velocity and acceleration of the link BD (link 2)*)phi2 = ArcTan[ (yB[t]-yC)/(xB[t]-xC) ] ;omega2 = D[ phi2 , t ] ;alpha2 = D[ omega2, t ] ;Print["phi2 = phi3 = ", phi2/.ruleB , " rad = ",(phi2/.ruleB)*180/N[Pi]," deg" ];Print["omega2 = omega3 = ", omega2/.ruleB , " rad/s" ];Print["alpha2 = alpha3 = ", alpha2/.ruleB , " rad/s^2" ];rB={xBs,yBs,0};rC={xC,yC,0};rBC=rC-rB;Print["rB = ", rB, " m" ];Print["rC = ", rC, " m" ];Print["rBC = rC - rB = ", rBC, " m" ];vB={vBxs,vBys,0};aB={aBxs,aBys,0};Print["vB = ", vB , " m/s" ];Print["aB = ", aB , " m/s" ];w2={0,0,omega2/.ruleB};a2={0,0,alpha2/.ruleB};Print["w2 = ", w2, " rad/s" ];Print["a2 = ", a2, " rad/s^2" ];vC2=vB+Cross[w2,rBC];Print["vC2 = vB + w2 x rBC = ",vC2, " m/s" ];aC2=aB+Cross[a2,rBC]-w2.w2 rBC;Print["aC2 = aB + a2 x rBC - w2^2 rBC = ", aC2, " m/s^2" ];vC32=-vC2;Print["vC32 = vC3 - vC2 = ",vC32, " m/s" ];aC32cor=2 Cross[w2,vC32];Print["aC32cor = 2 w2 x vC32 = ",aC32cor, " m/s^2" ];aC32=-aC2-aC32cor;Print["aC32 = aC3 - aC2 - aC32cor = ",aC32, " m/s^2" ]; (* Position of joint D *)eqnD1 = ( xD[t]- xB[t] )^2 + ( yD[t] - yB[t] )^2 - BD^2 == 0 ;eqnD2 = ( yD[t] - yB[t] ) / ( xD[t]- xB[t] ) == ( yB[t] - yC ) / ( xB[t] - xC );solutionD = Solve [ { eqnD1/.ruleB, eqnD2/.ruleB } , { xD[t] , yD[t] } ];(* Two solutions for D *)xD1 = xD[t] /. solutionD[[1]];yD1 = yD[t] /. solutionD[[1]];xD2 = xD[t] /. solutionD[[2]];yD2 = yD[t] /. solutionD[[2]];If [ yD1 <= 0 , xDs = xD1 ; yDs = yD1 , xDs = xD2 ; yDs=yD2 ] ;Print["xD = ",xDs," m"];Print["yD = ",yDs," m"];(* Velocity of joint D *)PD={xD[t]ØxDs,yD[t]ØyDs};eqnD1V = D[eqnD1,t] ;eqnD2V = D[eqnD2,t] ;solutionDV = Solve [ { (eqnD1V/.ruleB/.PD ),( eqnD2V/.ruleB/.PD )} , { xD'[t] , yD'[t] } ];vDxs = xD'[t] /. solutionDV[[1]];vDys = yD'[t] /. solutionDV[[1]];Print["vDx = ",vDxs," m/s"];Print["vDy = ",vDys," m/s"];(* Acceleration of joint D *)PVD={xD[t]ØxDs,yD[t]ØyDs,xD'[t]ØvDxs,yD'[t]ØvDys};eqnD1A = D[eqnD1V,t] ;eqnD2A = D[eqnD2V,t] ;solutionDA = Solve [ { (eqnD1A/.ruleB/.PVD ),( eqnD2A/.ruleB/.PVD )} , { xD''[t] , yD''[t] } ];aDxs = xD''[t] /. solutionDA[[1]];aDys = yD''[t] /. solutionDA[[1]];Print["aDx = ",aDxs," m/s^2"];Print["aDy = ",aDys," m/s^2"];ruleD={xD[t]ØxDs,yD[t]ØyDs,xD'[t]ØvDxs,yD'[t]ØvDys, xD''[t]ØaDxs,yD''[t]ØaDys};rD={xDs,yDs,0};rBD=rD-rB;Print["rD = ", rD, " m" ];Print["rBD = rD - rB = ", rBD, " m" ];vD=vB+Cross[w2,rBD];Print["vD = vB + w2 x rBD = ",vD, " m/s" ];aD=aB+Cross[a2,rBD]-w2.w2 rBD;Print["aD = aB + a2 x rBD - w2^2 rBD = ", aD, " m/s^2" ];(* Position of joint E *)eqnE1=(xE[t]-xD[t])^2+(yE[t]-yD[t])^2-DE^2ã0;eqnE2=(xE[t]-xF)^2+(yE[t]-yF)^2-EF^2ã0;solutionE=Solve[{eqnE1/.ruleD,eqnE2/.ruleD},{xE[t],yE[t]}];(*Two solutions for E*)xE1=xE[t]/.solutionE[[1]];yE1=yE[t]/.solutionE[[1]];xE2=xE[t]/.solutionE[[2]];yE2=yE[t]/.solutionE[[2]];(*Select the correct position for E*)If[yE1 ≤ 0,xEs=xE1;yEs=yE1,xEs=xE2;yEs=yE2];Print["xE = ",xEs," m"];Print["yE = ",yEs," m"];(* Velocity of joint E *)PE={xE[t]ØxEs,yE[t]ØyEs};eqnE1V = D[eqnE1,t] ;eqnE2V = D[eqnE2,t] ;solutionEV = Solve [ { (eqnE1V/.ruleD/.PE ),( eqnE2V/.ruleD/.PE )} , { xE'[t] , yE'[t] } ];vExs = xE'[t] /. solutionEV[[1]];vEys = yE'[t] /. solutionEV[[1]];Print["vEx = ",vExs," m/s"];Print["vEy = ",vEys," m/s"];(* Acceleration of joint E *)PVE={xE[t]ØxEs,yE[t]ØyEs,xE'[t]ØvExs,yE'[t]ØvEys};eqnE1A = D[eqnE1V,t] ;eqnE2A = D[eqnE2V,t] ;solutionEA = Solve [ { (eqnE1A/.ruleD/.PVE ),( eqnE2A/.ruleD/.PVE )} , { xE''[t] , yE''[t] } ];aExs = xE''[t] /. solutionEA[[1]];aEys = yE''[t] /. solutionEA[[1]];Print["aEx = ",aExs," m/s^2"];Print["aEy = ",aEys," m/s^2"];ruleE={xE[t]ØxEs,yE[t]ØyEs,xE'[t]ØvExs,yE'[t]ØvEys, xE''[t]ØaExs,yE''[t]ØaEys};(* Angular velocity and acceleration of the link DE (link 4)*)phi4 = ArcTan[ (yD[t]-yE[t])/(xD[t]-xE[t]) ] ;omega4 = D[ phi4 , t ] ;alpha4 = D[ omega4, t ] ;Print["phi4 = ", phi4/.ruleD/.ruleE , " rad = ",(phi4/.ruleD/.ruleE)*180/N[Pi]," deg" ];Print["omega4 = ", omega4/.ruleD/.ruleE , " rad/s" ];Print["alpha4 = ", alpha4/.ruleD/.ruleE , " rad/s^2" ];(* Angular velocity and acceleration of the link EF (link 5)*)phi5 = ArcTan[ (yE[t]-yF)/(xE[t]-xF) ] ;omega5 = D[ phi5 , t ] ;alpha5 = D[ omega5, t ] ;Print["phi5 = ", phi5/.ruleE , " rad = ",(phi5/.ruleE)*180/N[Pi]," deg" ];Print["omega5 = ", omega5/.ruleE , " rad/s" ];Print["alpha5 = ", alpha5/.ruleE , " rad/s^2" ];Problem4_8_II.nb 1(* Problem 4.8 *)(* VELOCITY and ACCELERATION ANALYSIS - input angle phi *)Apply[Clear,Names["Global`*"]];Off[General::spell];Off[General::spell1];(* Input data *)AB = .120 ;AC = .060 ;BD = .240 ;DE = .330 ;EF = .190 ;La = .300 ;Lb = .070 ;n = 100 ; (* rpm *)omega = n*N[Pi]/30 ; (* rad/s *)initdata = {f[t]->N[Pi]150/180, f'[t]->omega, f''[t]->0};(* Position of joint A *)xA = yA = 0;(* Position of joint C *)xC = -AC ;yC = 0 ;(* Position of joint F *)xF = La ;yF = Lb ;(* Position of joint B *)XB = AB Cos[f[t]] ;YB = AB Sin[f[t]] ;xBs=XB/.initdata;yBs=YB/.initdata;Print["xB[t] = ", XB," = ", xBs, " m" ];Print["yB[t] = ", YB," = ", yBs, " m" ];(* Linear velocity of joint B *)vBx = D[XB,t];vBy = D[YB,t];vBxs = vBx/.initdata;vBys = vBy/.initdata;Print["vBx = ", vBx ," = ", vBxs, " m/s" ];Print["vBy = ", vBy ," = ", vBys, " m/s" ];(* Linear acceleration of joint B *)aBx = D[vBx,t];aBy = D[vBy,t];aBxs =