<a href="https://colab.research.google.com/github/OscarSantos98/Dynamics/blob/master/Mechanisms/Solution_of_the_problem_6_15_of_Michael_Rider/Solution_of_the_problem_6_15_of_Michael_Rider.ipynb" target="_parent"></a>

# Solution of the problem 6-15 of Michael Rider book, using the methods of dynamics of rigid bodies.
# by **Oscar Alberto Santos Muñoz**

This code is licensed by MIT LICENSE

Copyright (c) 2020 OscarSantos98

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

<table class="buttons" align="center">
  <td>
    <a target="_blank" href="https://drive.google.com/open?id=1TMs1iiEcBBNUgmY8XnYxq3i8kOJpYyyc">Run in Google Colab</a>
  </td>
  <td>
    <a target="_blank" href="https://github.com/OscarSantos98/Dynamics">View source on GitHub</a>
  </td>
</table>


```
import sympy as sp
```





From the problem we know that $\omega2 = 9 \; \frac{rad}{s} \; c.c.w.$ and therefore $\alpha2 = 0 \; \frac{rad}{s^{2}} $

# In the next line the values of the magnitude of the vectors are introduced.


```
a2, a3, a4, a5, a6, b3 = 70, 100, 80, 120, 100, 65
```

Print the values of the magnitudes.


```
print("Links length of the mechanism\na2 =",a2,"\na3 =",a3,"\na4 =",a4,"\na5 =",a5,"\na6 =",a6,"\nb3 =",b3)
```

    Links length of the mechanism
    a2 = 70 
    a3 = 100 
    a4 = 80 
    a5 = 120 
    a6 = 100 
    b3 = 65



```
theta2 = (95*sp.pi/180).evalf()
theta3 = (4.63498*sp.pi/180).evalf()
theta4 = (76.57626*sp.pi/180).evalf()
theta5 = (-7.2327*sp.pi/180).evalf()
theta6 = (105.74221*sp.pi/180).evalf()
beta = (35.18382*sp.pi/180).evalf()
```

Create alpha, beta and theta strings.






```
th_str = []
w_str = []
alp_str = []
for i in range(2,7):
	th_str.append(chr(952) + str(i))
	w_str.append(chr(969) + str(i))
	alp_str.append(chr(945) + str(i))
```

Print the values of the angles.


```
print("Angles\n",th_str[0]," =",theta2,"\n",th_str[1]," =",theta3,"\n",th_str[2]," =",theta4,"\n",th_str[3]," =",theta5,"\n",th_str[4]," =",theta6,"\n",chr(946)," =",beta)
```

    Angles
     θ2  = 1.65806278939461 
     θ3  = 0.0808956617640868 
     θ4  = 1.33650786586323 
     θ5  = -0.126234428808994 
     θ6  = 1.84554972283527 
     β  = 0.614073502429031


# The necessary position vectors will be computed.


```
rCO21 = sp.Matrix(1,3,[a2*sp.cos(theta2),a2*sp.sin(theta2),0])
rDC = sp.Matrix(1,3,[a3*sp.cos(theta3),a3*sp.sin(theta3),0])
rDO41 = sp.Matrix(1,3,[a4*sp.cos(theta4),a4*sp.sin(theta4),0])
rPC = sp.Matrix(1,3,[b3*sp.cos(theta3+beta),b3*sp.sin(theta3+beta),0])
rFP = sp.Matrix(1,3,[a5*sp.cos(theta5),a5*sp.sin(theta5),0])
rFO61 = sp.Matrix(1,3,[a6*sp.cos(theta6),a6*sp.sin(theta6),0])
```

Print the results of the position vectors.


```
print("Position vectors\nrCO21 =",rCO21,"\nrDC =",rDC,"\nrDO41 =",rDO41,"\nrPC =",rPC,"\nrFP =",rFP,"\nrFO61 =",rFO61)
```

    Position vectors
    rCO21 = Matrix([[-6.10090199233608, 69.7336288664222, 0]]) 
    rDC = Matrix([[99.6729729954293, 8.08074589703349, 0]]) 
    rDO41 = Matrix([[18.5720755743718, 77.8143817610846, 0]]) 
    rPC = Matrix([[49.9247741348628, 41.6235141186197, 0]]) 
    rFP = Matrix([[119.045161106793, -15.1079322561889, 0]]) 
    rFO61 = Matrix([[-27.1309590837080, 96.2492132913208, 0]])


# Velocity analysis: Definition of the angular velocities.



Declare symbolic variables.


```
w2, w3, w4, w5, w6 = sp.symbols('w2 w3 w4 w5 w6')
```

Create vectors with the known values and the symbolic variables.


```
omega2 = sp.Matrix(1,3,[0,0,9])
omega3 = sp.Matrix(1,3,[0,0,w3])
omega4 = sp.Matrix(1,3,[0,0,w4])
omega5 = sp.Matrix(1,3,[0,0,w5])
omega6 = sp.Matrix(1,3,[0,0,w6])
```

Print the vectors.


```
print("Angular velocities\n",w_str[0],"=",omega2,"\n",w_str[1],"=",omega3,"\n",w_str[2],"=",omega4,"\n",w_str[3],"=",omega5,"\n",w_str[4],"=",omega6)
```

    Angular velocities
     ω2 = Matrix([[0, 0, 9]]) 
     ω3 = Matrix([[0, 0, w3]]) 
     ω4 = Matrix([[0, 0, w4]]) 
     ω5 = Matrix([[0, 0, w5]]) 
     ω6 = Matrix([[0, 0, w6]])


Calculate velocities.


```
VC = omega2.cross(rCO21)
VD1 = VC + omega3.cross(rDC)
VD2 = omega4.cross(rDO41)
```

Print the results of the point velocities


```
print("Point velocities\nVC = ",w_str[0]," x rCO21 = ",VC,"\nVD1 = VC + ",w_str[1]," x rDC = ",VD1,"\nVD2 = ",w_str[2]," x rDO41 = ",VD2)
```

    Point velocities
    VC =  ω2  x rCO21 =  Matrix([[-627.602659797800, -54.9081179310247, 0]]) 
    VD1 = VC +  ω3  x rDC =  Matrix([[-8.08074589703349*w3 - 627.6026597978, 99.6729729954293*w3 - 54.9081179310247, 0]]) 
    VD2 =  ω4  x rDO41 =  Matrix([[-77.8143817610846*w4, 18.5720755743718*w4, 0]])


Solve for w3 and w4.


```
Eqv1 = sp.Eq(VD1[0],VD2[0])
Eqv2 = sp.Eq(VD1[1],VD2[1])
answ = sp.solve([Eqv1,Eqv2])
```

Print the answer.


```
print("Solution for w3 and w4 of matching velocity equations of point D\n",answ)
```

    Solution for w3 and w4 of matching velocity equations of point D
     {w3: 2.09422879753408, w4: 8.28285949170457}


Substitue the answer in the vector w3 and w4.


```
omega3 = omega3.subs(w3,answ[w3])
omega4 = omega4.subs(w4,answ[w4])
```


```
print("Substituting w3 and w4 in vectors",w_str[1]," and ",w_str[2],"\n",w_str[1],"=",omega3,"\n",w_str[2],"=",omega4)
```

    Substituting w3 and w4 in vectors ω3  and  ω4 
     ω3 = Matrix([[0, 0, 2.09422879753408]]) 
     ω4 = Matrix([[0, 0, 8.28285949170457]])


Calculate the rest of the velocities.


```
VP = VC + omega3.cross(rPC)
VF1 = VP + omega5.cross(rFP)
VF2 = omega6.cross(rFO61)
```

Print the result of the velocities.


```
print("Point velocities\nVP = VC + ",w_str[1]," x rPC = ",VP,"\nVF1 = VP + ",w_str[3]," x rFP = ",VF1,"\nVF2 = ",w_str[4]," x rFO61 = ",VF2)
```

    Point velocities
    VP = VC +  ω3  x rPC =  Matrix([[-714.771821719580, 49.6457817725896, 0]]) 
    VF1 = VP +  ω5  x rFP =  Matrix([[15.1079322561889*w5 - 714.77182171958, 119.045161106793*w5 + 49.6457817725896, 0]]) 
    VF2 =  ω6  x rFO61 =  Matrix([[-96.2492132913208*w6, -27.130959083708*w6, 0]])


Solve for w5 and w6.




```
Eqv3 = sp.Eq(VF1[0],VF2[0])
Eqv4 = sp.Eq(VF1[1],VF2[1])
answ2 = sp.solve([Eqv3,Eqv4])
```

Print the answer.


```
print("Solution for w5 and w6 of matching velocity equations of point F\n",answ2)
```

    Solution for w5 and w6 of matching velocity equations of point F
     {w5: -2.18777796287317, w6: 7.76966997860839}


Substitue the answer in the vector w5 and w5.


```
omega5 = omega5.subs(w5,answ2[w5])
omega6 = omega6.subs(w6,answ2[w6])
```


```
print("Substituting w5 and w6 in vectors",w_str[3]," and ",w_str[4],"\n",w_str[3],"=",omega5,"\n",w_str[4],"=",omega6)
```

    Substituting w5 and w6 in vectors ω5  and  ω6 
     ω5 = Matrix([[0, 0, -2.18777796287317]]) 
     ω6 = Matrix([[0, 0, 7.76966997860839]])


# Acceleration analysis. Definition of the angular accelerations.



Declare symbolic variables.


```
aa2, aa3, aa4, aa5, aa6 = sp.symbols('aa2 aa3 aa4 aa5 aa6')
```

Create vectors with the known values and the symbolic variables.


```
alpha2 = sp.Matrix(1,3,[0,0,0])
alpha3 = sp.Matrix(1,3,[0,0,aa3])
alpha4 = sp.Matrix(1,3,[0,0,aa4])
alpha5 = sp.Matrix(1,3,[0,0,aa5])
alpha6 = sp.Matrix(1,3,[0,0,aa6])
```

Print the vectors.


```
print("Angular accelerations\n",alp_str[0],"=",alpha2,"\n",alp_str[1],"=",alpha3,"\n",alp_str[2],"=",alpha4,"\n",alp_str[3],"=",alpha5,"\n",alp_str[4],"=",alpha6)
```

    Angular accelerations
     α2 = Matrix([[0, 0, 0]]) 
     α3 = Matrix([[0, 0, aa3]]) 
     α4 = Matrix([[0, 0, aa4]]) 
     α5 = Matrix([[0, 0, aa5]]) 
     α6 = Matrix([[0, 0, aa6]])


Calculate accelerations.


```
AC = alpha2.cross(rCO21) + omega2.cross(omega2.cross(rCO21))
AD3 = AC + alpha3.cross(rDC) + omega3.cross(omega3.cross(rDC))
AD4 = alpha4.cross(rDO41) + omega4.cross(omega4.cross(rDO41))
```

Print the results of the accelerations.


```
print("Point accelerations\nAC = ",alp_str[0]," x rCO21 + ",w_str[0]," x ",w_str[0]," x rCO21 = ",AC,"\nAD3 = AC + ",alp_str[1]," x rDC + ",w_str[1]," x ",w_str[1]," x rDC = ",AD3,"\nAD4 = ",alp_str[2]," x rDO41 + ",w_str[2]," x ",w_str[2]," x rDO41 = ",AD4)
```

    Point accelerations
    AC =  α2  x rCO21 +  ω2  x  ω2  x rCO21 =  Matrix([[494.173061379222, -5648.42393818020, 0]]) 
    AD3 = AC +  α3  x rDC +  ω3  x  ω3  x rDC =  Matrix([[-8.08074589703349*aa3 + 57.0279088954586, 99.6729729954293*aa3 - 5683.86442712301, 0]]) 
    AD4 =  α4  x rDO41 +  ω4  x  ω4  x rDO41 =  Matrix([[-77.8143817610846*aa4 - 1274.15138480262, 18.5720755743718*aa4 - 5338.51490542404, 0]])


Solve for aa3 and aa4


```
Eqa1 = sp.Eq(AD3[0],AD4[0])
Eqa2 = sp.Eq(AD3[1],AD4[1])
ansa = sp.solve([Eqa1,Eqa2])
```

Print the answer.


```
print("Solution for aa3 and aa4 of matching accelerations equations of point D\n",ansa)
```

    Solution for aa3 and aa4 of matching accelerations equations of point D
     {aa3: 0.282726957185647, aa4: -17.0777511678876}


Substitute the answer in the vectors aa3 and aa4


```
alpha3 = alpha3.subs(aa3,ansa[aa3])
alpha4 = alpha4.subs(aa4,ansa[aa4])
```


```
print("Substituting aa3 and aa4 in vectors",alp_str[1]," and ",alp_str[2],"\n",alp_str[1],"=",alpha3,"\n",alp_str[2],"=",alpha4)
```

    Substituting aa3 and aa4 in vectors α3  and  α4 
     α3 = Matrix([[0, 0, 0.282726957185647]]) 
     α4 = Matrix([[0, 0, -17.0777511678876]])


Calculate the rest of the accelerations.


```
AP3 = AC + alpha3.cross(rPC) + omega3.cross(omega3.cross(rPC))
AF5 = AP3 + alpha5.cross(rFP) + omega5.cross(omega5.cross(rFP))
AF6 = alpha6.cross(rFO61) + omega6.cross(omega6.cross(rFO61))
```

Print the results of the accelerations.


```
print("Point accelerations\nAP3 = AC + ",alp_str[1]," x rPC + ",w_str[1]," x ",w_str[1]," x rPC = ",AP3,"\nAF5 = AP3 + ",alp_str[3]," x rFP + ",w_str[3]," x ",w_str[3]," x rFP = ",AF5,"\nAF6 = ",alp_str[4]," x rFO61 + ",w_str[4]," x ",w_str[4]," x rFO61 = ",AF6)
```

    Point accelerations
    AP3 = AC +  α3  x rPC +  ω3  x  ω3  x rPC =  Matrix([[263.445184231292, -5816.86102785437, 0]]) 
    AF5 = AP3 +  α5  x rFP +  ω5  x  ω5  x rFP =  Matrix([[15.1079322561889*aa5 - 306.349291009673, 119.045161106793*aa5 - 5744.54883765817, 0]]) 
    AF6 =  α6  x rFO61 +  ω6  x  ω6  x rFO61 =  Matrix([[-96.2492132913208*aa6 + 1637.83554061634, -27.130959083708*aa6 - 5810.35052238718, 0]])


Solve for aa5 and aa6.


```
Eqa3 = sp.Eq(AF5[0],AF6[0])
Eqa4 = sp.Eq(AF5[1],AF6[1])
ansa2 = sp.solve([Eqa3,Eqa4])
```

Print the answer.


```
print("Solution for aa5 and aa6 of matching accelerations equations of point F\n",ansa2)
```

    Solution for aa5 and aa6 of matching accelerations equations of point F
     {aa5: -5.34760762788408, aa6: 21.0388849545367}


Substitute the answer in the vectors aa5 and aa6


```
alpha5 = alpha5.subs(aa5,ansa2[aa5])
alpha6 = alpha6.subs(aa6,ansa2[aa6])
```


```
print("Substituting aa5 and aa6 in vectors",alp_str[3]," and ",alp_str[4],"\n",alp_str[3],"=",alpha5,"\n",alp_str[4],"=",alpha6)
```

    Substituting aa5 and aa6 in vectors α5  and  α6 
     α5 = Matrix([[0, 0, -5.34760762788408]]) 
     α6 = Matrix([[0, 0, 21.0388849545367]])

