In [1]:
# Initial import statements
%matplotlib notebook
import matplotlib.pyplot as plt
import numpy as np
import sympy as sy
from matplotlib.pyplot import *
from numpy import *
from numpy.linalg import *
from sympy import *
Note: I assumed there is a mistake in the assignment and instead I copied the equations from the book. $µ_4$ is not defined but it was used in the equations. The results with the given expressions made no sense or I was just not able to calculate them accordingly.
Four blocks of different masses $m_i$ are connected by ropes of negligble mass. Three of the blocks lie on a inclined plane, the coefficients of friction between the blocks and the plane being $µ_i$. The equations of motion for the blocks can be shown to be:
where $T_i$ denotes the tensile forces in the ropes and a is the acceleration of the system.
(a) Determine $a$ and $T_i$, when $θ = 45^{\circ}$ and $g = 9.81 m/s^2$, $m = [10.0, 4.0, 5.0, 6.0]$ kg, and $µ = [0.25, 0.30, 0.20]$.
(b) What the angle should be in order that the system is in balance? Try a couple of different values for angle and find out what are the values for $a$ and $T_i$. Based on these values make a graph (x-axis = angle = θ, y-axis = acceleration = a) and based on the graph estimate the angle giving the acceleration $a = 0.0 m/s^2$.
Before I actually start to create the code, it might be a good idea to convert the equations to matrix form:
$$ \begin{bmatrix} 1.0 & 0.0 & 0.0 \\ -1.0 & 1.0 & 0.0 \\ 0.0 & -1.0 & 1.0 \end{bmatrix} \begin{bmatrix} T_1 \\ T_2 \\ T_3 \end{bmatrix} = \begin{bmatrix} m_1g(sin θ - µ_1 cos θ) - m_1a\\ m_2g(sin θ - µ_2 cos θ) - m_2a\\ m_3g(sin θ - µ_3 cos θ) - m_3a \end{bmatrix} $$...but from that point forward I have no idea how to continue with matrices, so I decided to use the sympy.solve() function to find out the equations for t1, t2 and t3 to be used when calculating a. Once I have a numerical value for a, I can use it to find out the values for t1, t2 and t3 using basic matrix calculations.
In [10]:
# The value for gravity doesn't change as long as we're in on earth, so it can be defined "globally"
gravity = 9.81
# Define characteristics of our application
mass = np.array(([10.0, 4.0, 5.0, 6.0]))
friction = np.array(([0.25, 0.30, 0.20]))
# Calculate value for acceleration
def calculate_acc(theta):
a = Symbol('a')
t1 = Symbol('t1')
temp = sy.solve([t1 + (mass[0] * a) - (mass[0] * gravity * (sin(theta) - friction[0] * (cos(theta))))], [t1, a])
t1 = temp[t1]
t2 = Symbol('t2')
temp = sy.solve([-t1 + t2 + (mass[1] * a) - (mass[1] * gravity * (sin(theta) - friction[1] * (cos(theta))))], [t2, a])
t2 = temp[t2]
t3 = Symbol('t3')
temp = sy.solve([-t2 + t3 + (mass[2] * a) - (mass[2] * gravity * (sin(theta) - friction[2] * (cos(theta))))], [t3, a])
t3 = temp[t3]
acc = np.array(zeros(1)) # Create array as the result sympy solve is in list format
acc = sy.solve(-(t3) + mass[3]*a + mass[3]*gravity, a)
return acc[0]
a = calculate_acc(np.deg2rad(45.0))
def calculate_t(theta):
A = np.array(([1.0, 0.0, 0.0],
[-1.0, 1.0, 0.0],
[0.0, -1.0, 1.0]))
b = np.array(([mass[0] * gravity * (sin(theta) - friction[0] * (cos(theta))) - (mass[0]*a)],
[mass[1] * gravity * (sin(theta) - friction[1] * (cos(theta))) - (mass[0]*a)],
[mass[2] * gravity * (sin(theta) - friction[2] * (cos(theta))) - (mass[0]*a)]))
t = dot(inv(A), b)
return t
t = calculate_t(np.deg2rad(45.0))
print("Value of a is:", a)
print("Values for t are:\nt1 = {}\nt2 = {}\nt3 = {}".format(t[0], t[1], t[2]))
Now that we have values for a, t1, t2 and t3 with the angle of 45$^{\circ}$, we can start to iterate to find out the angle in which is the application is in balance. The application is in balance when acceleration equals 0.0 (e.g. it doesn't move one way or another).
In [9]:
%matplotlib notebook
iterations = 100
step = -1
theta = np.deg2rad(45.0)
# Create empty arrays to store values into
x = []
y = []
for i in range(iterations):
aOld = a.copy
# Increment theta by step
theta = degrees(theta)
theta += step
theta = deg2rad(theta)
# Perform calculations with new theta
aOld = calculate_acc(theta)
t = calculate_t(theta)
x.append(degrees(theta))
y.append(aOld)
# Break if we go under 0.0
if aOld <= 0.0:
break
a = aOld
print("Value of a at break point is: {}".format(a))
print("Values for t are:\nt1 = {}\nt2 = {}\nt3 = {}".format(t[0], t[1], t[2]))
# Create plot
figure('Problem 1')
plt.plot(x, y)
plt.ylabel("Acceleration (m/s^2)")
plt.xlabel("Angle (degrees)")
plt.grid()
plt.show()
From the graph we can see to get acceleration of $0.0 m/s^2$ the angle must be roughly 31.7457 degrees. We can confirm this by using the functions previously created.
In [11]:
angle = deg2rad(31.7457)
a = calculate_acc(angle)
calculate_t(angle)
print("Value of a is:", a)
print("Values for t are:\nt1 = {}\nt2 = {}\nt3 = {}".format(t[0], t[1], t[2]))
So as we can see, the value of a is close to zero at the given angle.