Free Body Diagram for particles

Renato Naville Watanabe

In [2]:
import numpy as np
import matplotlib.pyplot as plt

%matplotlib notebook

Steps to draw a free-body diagram

1 - Draw separately each object considered in the problem. How you separate depends on what questions you want to answer.

2 - Identify the forces acting on each object. If you are analyzing more than one object, remember the third Newton law (action and reaction), and identify where the reaction of a force is being applied.

3 - Draw all the identified forces, representing them as vectors. The vectors should be represented with the origin in the object. In the case of particles, the origin should be in the center of the particle.

4 - If necessary, you should represent the reference frame in the free-body diagram.

5 - After this, you can solve the problem using the Second Newton Law (see, e.g, Newton Laws) to find the motion of the particle.

Basic element and forces

Spring

Is an element used to represent a force proportional to some length or displacement. It produces a force in the same direction of the vector linking the spring extremities and opposite to its length or displacement from an equilibrium length. Frequently it has a linear relation, but it could be nonlinear as well. The force exerted by the spring in one of the extremities is:

\begin{equation} \vec{\bf{F}} = - k(||\vec{\bf{r}}||-l_0)\frac{\vec{\bf{r}}}{||\vec{\bf{r}}||} = -k\vec{\bf{r}} +kl_0\frac{\vec{\bf{r}}}{||\vec{\bf{r}}||} = -k\left(1-\frac{l_0}{||\vec{\bf{r}}||}\right)\vec{\bf{r}} \end{equation}

where $\vec{\bf{r}}$ is the vector linking the extremity applying the force to the other extremity and $l_0$ is the equilibrium length of the spring.

Since the spring element is a massless element, the force in both extremities have the same absolute value and opposite directions.

Damping

Is an element used to represent a force proportional to the velocity of displacement. It produces a force in the same direction of the vector linking the element extremities and opposite to its velocity.

Frequently it has a linear relation, but it could be nonlinear as well. The force exerted by the damping element in one of its extremities is:

\begin{equation} \vec{\bf{F}} = - b||\vec{\bf{v}}||\frac{\vec{\bf{v}}}{||\vec{\bf{v}}||} = -b\vec{\bf{v}} = -b\frac{d\vec{\bf{r}}}{dt} \end{equation}

where $\vec{\bf{r}}$ is the vector linking the extremity applying the force to the other extremity.

For the same reason of the spring, the force in both extremities have the same absolute value and opposite directions.

Gravity

The gravity force acts on two masses, each one atracting each other:

\begin{equation} \vec{{\bf{F}}} = - G\frac{m_1m_2}{||\vec{\bf{r}}||^2}\frac{\vec{\bf{r}}}{||\vec{\bf{r}}||} \end{equation}

where $G = 6.67.10^{−11} Nm^2/kg^2$ and $\vec{\bf{r}}$ is a vector with length equal to the distance between the masses and directing towards the other mass. Note the forces acting on each mass have the same absolute value.

Since the mass of the Earth is $m1=5.9736×10^24kg$ and its radius is 6.371×10^6m (see this notebook, for details), the gravity force on Earth is:

\begin{equation} \vec{{\bf{F}}} = m\vec{\bf{g}} \end{equation}

with the absolute value of $\vec{\bf{g}}$ approximately equal to 9.81 $m/s^2$, and point towards the center of Earth.

Below, will be shown some examples of how to draw the free-body diagram and obtain the motion equations to solve the problems.

1) No force acting on the particle

The most trivial situation is a particle with no force acting on it.

The free-body diagram is below, with no force vectors acting on the particle.

In this situation, the resultant force is:

\begin{equation} \vec{\bf{F}} = \vec{\bf{0}} \end{equation}

And the second Newton law for this particle is:

\begin{equation} m\frac{d^2\vec{\bf{r}}}{dt^2} = \vec{\bf{0}} \rightarrow \frac{d^2\vec{\bf{r}}}{dt^2} = \vec{\bf{0}} \end{equation}

The motion of of the particle can be found by integrating twice both times, getting the following:

\begin{equation} \vec{\bf{r}} = \vec{\bf{v_0}}t + \vec{\bf{r0}} \end{equation}

The particle continues to change its position with the same velocity it was at the beginning of the analysis. This could be predicted by the first Newton law.

2) Gravity force acting on the particle

Now, we will consider a ball with the gravity force acting on it. The free-body diagram is depicted below.

The only force acting on the ball is the gravitational force:

\begin{equation} \hat{\bf{F_g}} = - mg \hat{\bf{j}} \end{equation}

So, we write the Second Newton Law:

\begin{equation} \hat{\bf{F_g}} = m \frac{d^2\vec{\bf{r}}}{dt^2} \rightarrow - mg \hat{\bf{j}} = m \frac{d^2\vec{\bf{r}}}{dt^2} \rightarrow - g \hat{\bf{j}} =\frac{d^2\vec{\bf{r}}}{dt^2} \end{equation}

Now, we can separate the equation in two components (x and y):

\begin{equation} 0 = \frac{d^2x}{dt^2} \end{equation}

and

\begin{equation}

  • g = \frac{d^2y}{dt^2} \end{equation}

These equations were solved in detail in this Notebook about Newton laws.

3) Ground reaction force

Now, we will analyze the situation of a particle in contact with the ground. To simplify the analysis, only the vertical movement will be considered.

The forces acting on the particle are the ground reaction force (often called as normal force) and the gravity force. The free-body diagram of the particle is below:

So, the resultant force in the particle is:

\begin{equation} \vec{\bf{F}} = \overrightarrow{\bf{GRF}} + m\vec{\bf{g}} = \overrightarrow{\bf{GRF}} - mg\hat{\bf{j}} \end{equation}

Considering only the y direction:

\begin{equation} F = GRF - mg \end{equation}

Applying the second Newton law to the particle:

\begin{equation} m \frac{d^2y}{dt^2} = GRF - mg \end{equation}

Note that since we have no information about how the force GRF varies along time, we cannot solve this equation. To find the position of the particle along time, one would have to measure the ground reaction force. See the notebook on Vertical jump for an application of this model.

4) Linear spring in horizontal movement

The example below is mass attached to a spring. The other extremity of the spring is fixed.

The only force force acting on the mass is from the spring. Below is the free-body diagram from the mass.

Since the movement is horizontal, we can neglect the gravity force.

\begin{equation} \vec{\bf{F}} = -k\left(1-\frac{l_0}{||\vec{\bf{r}}||}\right)\vec{\bf{r}} \end{equation}

Applying the second Newton law to the mass:

\begin{equation} m\frac{d^2\vec{\bf{r}}}{dt^2} = -k\left(1-\frac{l_0}{||\vec{\bf{r}}||}\right)\vec{\bf{r}} \rightarrow \frac{d^2\vec{\bf{r}}}{dt^2} = -\frac{k}{m}\left(1-\frac{l_0}{||\vec{\bf{r}}||}\right)\vec{\bf{r}} \end{equation}

Since the movement is unidimensional, we can deal with it scalarly:

\begin{equation} \frac{d^2x}{dt^2} = -\frac{k}{m}\left(1-\frac{l_0}{x}\right)x = -\frac{k}{m}(x-l_0) \end{equation}

To solve this equation numerically, we must break the equations into two first-order differential equation:

\begin{equation} \frac{dv_x}{dt} = -\frac{k}{m}(x-l_0) \end{equation}

\begin{equation} \frac{dx}{dt} = v_x \end{equation}

In the numerical solution below, we will use $k = 40 N/m$, $m = 2 kg$, $l_0 = 0.5 m$ and the mass starts from the position $x = 0.8m$ and at rest.

In [3]:
k = 40
m = 2
l0 = 0.5
x0 = 0.8
v0 = 0

x = x0
v = v0

dt = 0.001
t = np.arange(0, 3, dt)

r = np.array([x])

for i in t[1:]:
    dxdt = v
    dvxdt = -k/m*(x-l0)
    x = x + dt*dxdt
    v = v + dt*dvxdt
    r = np.vstack((r,np.array([x])))
    

plt.figure()
plt.plot(t,r)
plt.show()

5) Linear spring in bidimensional movement in horizontal plane

This example has two masses attached to the spring. To solve the motion of both masses, we must draw a free-body diagram for each one of the masses.

The only force acting on each mass is the force due to the spring. Since the movement is happening in the horizontal plane, the gravity force can be neglected.

So, the forces acting on mass 1 is:

\begin{equation} \vec{\bf{F_1}} = k\left(||\vec{\bf{r_2}}-\vec{\bf{r_1}}||-l_0\right)\frac{(\vec{\bf{r_2}}-\vec{\bf{r_1}})}{||\vec{\bf{r_2}}-\vec{\bf{r_1}}||} \end{equation}

and the forces acting on mass 2 is:

\begin{equation} \vec{\bf{F_2}} =k\left(||\vec{\bf{r_2}}-\vec{\bf{r_1}}||-l_0\right)\frac{(\vec{\bf{r_1}}-\vec{\bf{r_2}})}{||\vec{\bf{r_2}}-\vec{\bf{r_1}}||} \end{equation}

Applying the second Newton law for the masses:

\begin{equation} m_1\frac{d^2\vec{\bf{r_1}}}{dt^2} = k\left(||\vec{\bf{r_2}}-\vec{\bf{r_1}}||-l_0\right)\frac{(\vec{\bf{r_2}}-\vec{\bf{r_1}})}{||\vec{\bf{r_2}}-\vec{\bf{r_1}}||} \rightarrow \frac{d^2\vec{\bf{r_1}}}{dt^2} = -\frac{k}{m_1}\left(1-\frac{l_0}{||\vec{\bf{r_2}}-\vec{\bf{r_1}}||}\right)\vec{\bf{r_1}}+\frac{k}{m_1}\left(1-\frac{l_0}{||\vec{\bf{r_2}}-\vec{\bf{r_1}}||}\right)\vec{\bf{r_2}} \rightarrow \frac{d^2x_1\hat{\bf{i}}}{dt^2}+\frac{d^2y_1\hat{\bf{j}}}{dt^2} = -\frac{k}{m_1}\left(1-\frac{l_0}{||\vec{\bf{r_2}}-\vec{\bf{r_1}}||}\right)(x_1\hat{\bf{i}}+y_1\hat{\bf{j}})+\frac{k}{m_1}\left(1-\frac{l_0}{||\vec{\bf{r_2}}-\vec{\bf{r_1}}||}\right)(x_2\hat{\bf{i}}+y_2\hat{\bf{j}}) \end{equation}

\begin{equation} m_2\frac{d^2\vec{\bf{r_2}}}{dt^2} = k\left(||\vec{\bf{r_2}}-\vec{\bf{r_1}}||-l_0\right)\frac{(\vec{\bf{r_1}}-\vec{\bf{r_2}})}{||\vec{\bf{r_2}}-\vec{\bf{r_1}}||} \rightarrow \frac{d^2\vec{\bf{r_2}}}{dt^2} = -\frac{k}{m_2}\left(1-\frac{l_0}{||\vec{\bf{r_2}}-\vec{\bf{r_1}}||}\right)\vec{\bf{r_2}}+\frac{k}{m_2}\left(1-\frac{l_0}{||\vec{\bf{r_2}}-\vec{\bf{r_1}}||}\right)\vec{\bf{r_1}} \rightarrow \frac{d^2x_2\hat{\bf{i}}}{dt^2}+\frac{d^2y_2\hat{\bf{j}}}{dt^2} = -\frac{k}{m_2}\left(1-\frac{l_0}{||\vec{\bf{r_2}}-\vec{\bf{r_1}}||}\right)(x_2\hat{\bf{i}}+y_2\hat{\bf{j}})+\frac{k}{m_2}\left(1-\frac{l_0}{||\vec{\bf{r_2}}-\vec{\bf{r_1}}||}\right)(x_1\hat{\bf{i}}+y_1\hat{\bf{j}}) \end{equation}

Now, we can separate the equations for each of the coordinates:

\begin{equation} \frac{d^2x_1}{dt^2} = -\frac{k}{m_1}\left(1-\frac{l_0}{||\vec{\bf{r_2}}-\vec{\bf{r_1}}||}\right)x_1+\frac{k}{m_1}\left(1-\frac{l_0}{||\vec{\bf{r_2}}-\vec{\bf{r_1}}||}\right)x_2=-\frac{k}{m_1}\left(1-\frac{l_0}{\sqrt{(x_2-x_1)^2+(y_2-y_1)^2}}\right)(x_1-x_2) \end{equation}

\begin{equation} \frac{d^2y_1}{dt^2} = -\frac{k}{m_1}\left(1-\frac{l_0}{||\vec{\bf{r_2}}-\vec{\bf{r_1}}||}\right)y_1+\frac{k}{m_1}\left(1-\frac{l_0}{||\vec{\bf{r_2}}-\vec{\bf{r_1}}||}\right)y_2=-\frac{k}{m_1}\left(1-\frac{l_0}{\sqrt{(x_2-x_1)^2+(y_2-y_1)^2}}\right)(y_1-y_2) \end{equation}

\begin{equation} \frac{d^2x_2}{dt^2} = -\frac{k}{m_2}\left(1-\frac{l_0}{||\vec{\bf{r_2}}-\vec{\bf{r_1}}||}\right)x_2+\frac{k}{m_2}\left(1-\frac{l_0}{||\vec{\bf{r_2}}-\vec{\bf{r_1}}||}\right)x_1=-\frac{k}{m_2}\left(1-\frac{l_0}{\sqrt{(x_2-x_1)^2+(y_2-y_1)^2}}\right)(x_2-x_1) \end{equation}

\begin{equation} \frac{d^2y_2}{dt^2} = -\frac{k}{m_2}\left(1-\frac{l_0}{||\vec{\bf{r_2}}-\vec{\bf{r_1}}||}\right)y_2+\frac{k}{m_2}\left(1-\frac{l_0}{||\vec{\bf{r_2}}-\vec{\bf{r_1}}||}\right)y_1=-\frac{k}{m_2}\left(1-\frac{l_0}{\sqrt{(x_2-x_1)^2+(y_2-y_1)^2}}\right)(y_2-y_1) \end{equation}

To solve these equations numerically, you must break these equations into first-order equations:

\begin{equation} \frac{dv_{x_1}}{dt} = -\frac{k}{m_1}\left(1-\frac{l_0}{\sqrt{(x_2-x_1)^2+(y_2-y_1)^2}}\right)(x_1-x_2) \end{equation}

\begin{equation} \frac{dv_{y_1}}{dt} = -\frac{k}{m_1}\left(1-\frac{l_0}{\sqrt{(x_2-x_1)^2+(y_2-y_1)^2}}\right)(y_1-y_2) \end{equation}

\begin{equation} \frac{dv_{x_2}}{dt} = -\frac{k}{m_2}\left(1-\frac{l_0}{\sqrt{(x_2-x_1)^2+(y_2-y_1)^2}}\right)(x_2-x_1) \end{equation}

\begin{equation} \frac{dv_{y_2}}{dt} = -\frac{k}{m_2}\left(1-\frac{l_0}{\sqrt{(x_2-x_1)^2+(y_2-y_1)^2}}\right)(y_2-y_1) \end{equation}

\begin{equation} \frac{dx_1}{dt} = v_{x_1} \end{equation}

\begin{equation} \frac{dy_1}{dt} = v_{y_1} \end{equation}

\begin{equation} \frac{dx_2}{dt} = v_{x_2} \end{equation}

\begin{equation} \frac{dy_2}{dt} = v_{y_2} \end{equation}

Note that if you did not wanted to know the details about the motion of each mass, but only the motion of the center of mass of the masses-spring system, you could have modeled the whole system as a single particle.

To solve the equations numerically, we will use the $m_1=1 kg$, $m_2 = 2 kg$, $l_0 = 0.5 m$, $k = 90 N/m$ and $x_{1_0} = 0 m$, $x_{2_0} = 0 m$, $y_{1_0} = 1 m$, $y_{2_0} = -1 m$, $v_{x1_0} = -2 m/s$, $v_{x2_0} = 0 m$, $v_{y1_0} = 0 m$, $v_{y2_0} = 0 m$.

In [4]:
x01 = 0
y01= 0.5
x02 = 0
y02 = -0.5
vx01 = 0.1
vy01 = 0
vx02 = -0.1
vy02 = 0

x1= x01
y1 = y01
x2= x02
y2 = y02
vx1= vx01
vy1 = vy01
vx2= vx02
vy2 = vy02
r1 = np.array([x1,y1])
r2 = np.array([x2,y2])

k = 30
m1 = 1
m2 = 1
l0 = 0.5

dt = 0.0001
t = np.arange(0,5,dt)

for i in t[1:]:
    
    dvx1dt = -k/m1*(x1-x2)*(1-l0/np.sqrt((x2-x1)**2+(y2-y1)**2))
    dvx2dt = -k/m2*(x2-x1)*(1-l0/np.sqrt((x2-x1)**2+(y2-y1)**2))
    dvy1dt = -k/m1*(y1-y2)*(1-l0/np.sqrt((x2-x1)**2+(y2-y1)**2))
    dvy2dt = -k/m2*(y2-y1)*(1-l0/np.sqrt((x2-x1)**2+(y2-y1)**2))
    
    dx1dt = vx1
    dx2dt = vx2
    dy1dt = vy1
    dy2dt = vy2
    
    x1 = x1 + dt*dx1dt
    x2 = x2 + dt*dx2dt
    y1 = y1 + dt*dy1dt
    y2 = y2 + dt*dy2dt
    
    vx1 = vx1 + dt*dvx1dt
    vx2 = vx2 + dt*dvx2dt
    vy1 = vy1 + dt*dvy1dt
    vy2 = vy2 + dt*dvy2dt
    
    r1 = np.vstack((r1,np.array([x1,y1])))
    r2 = np.vstack((r2,np.array([x2,y2])))

springLength = np.sqrt((r1[:,0]-r2[:,0])**2+(r1[:,1]-r2[:,1])**2)    
plt.figure()
plt.plot(t, springLength)
plt.show()

6)Particle with gravity and linear air resistance

Below is the free-body diagram of a particle with the gravity force and a linear drag force due to the air resistance.

the forces being applied in the ball are:

\begin{equation} \vec{\bf{F}} = -mg \hat{\bf{j}} - b\vec{\bf{v}} = -mg \hat{\bf{j}} - b\frac{d\vec{\bf{r}}}{dt} = -mg \hat{\bf{j}} - b\left(\frac{dx}{dt}\hat{\bf{i}}+\frac{dy}{dt}\hat{\bf{j}}\right) = - b\frac{dx}{dt}\hat{\bf{i}} - \left(mg + b\frac{dy}{dt}\right)\hat{\bf{j}} \end{equation}

Writing down the Second Newton Law:

\begin{equation} \vec{\bf{F}} = m \frac{d^2\vec{\bf{r}}}{dt^2} \rightarrow - b\frac{dx}{dt}\hat{\bf{i}} - \left(mg + b\frac{dy}{dt}\right)\hat{\bf{j}} = m\left(\frac{d^2x}{dt^2}\hat{\bf{i}}+\frac{d^2y}{dt^2}\hat{\bf{j}}\right) \end{equation}

Now, we can separate into one equation for each coordinate:

\begin{equation}

  • b\frac{dx}{dt} = m\frac{d^2x}{dt^2} -\rightarrow \frac{d^2x}{dt^2} = -\frac{b}{m} \frac{dx}{dt} \end{equation}

\begin{equation} -mg - b\frac{dy}{dt} = m\frac{d^2y}{dt^2} \rightarrow \frac{d^2y}{dt^2} = -\frac{b}{m}\frac{dy}{dt} - g \end{equation}

These equations were solved in this notebook.

7) Particle with gravity and nonlinear air resistance

Below, is the free-body diagram of a particle with the gravity force and a drag force due to the air resistance proportional to the square of the particle velocity.

The forces being applied in the ball are:

\begin{equation} \vec{\bf{F}} = -mg \hat{\bf{j}} - bv^2\hat{\bf{e_t}} = -mg \hat{\bf{j}} - b (v_x^2+v_y^2) \frac{v_x\hat{\bf{i}}+v_y\hat{\bf{j}}}{\sqrt{v_x^2+v_y^2}} = -mg \hat{\bf{j}} - b \sqrt{v_x^2+v_y^2} \,(v_x\hat{\bf{i}}+v_y\hat{\bf{j}}) = -mg \hat{\bf{j}} - b \sqrt{\left(\frac{dx}{dt} \right)^2+\left(\frac{dy}{dt} \right)^2} \,\left(\frac{dx}{dt} \hat{\bf{i}}+\frac{dy}{dt}\hat{\bf{j}}\right) \end{equation}

Writing down the Second Newton Law:

\begin{equation} \vec{\bf{F}} = m \frac{d^2\vec{\bf{r}}}{dt^2} \rightarrow -mg \hat{\bf{j}} - b \sqrt{\left(\frac{dx}{dt} \right)^2+\left(\frac{dy}{dt} \right)^2} \,\left(\frac{dx}{dt} \hat{\bf{i}}+\frac{dy}{dt}\hat{\bf{j}}\right) = m\left(\frac{d^2x}{dt^2}\hat{\bf{i}}+\frac{d^2y}{dt^2}\hat{\bf{j}}\right) \end{equation}

Now, we can separate into one equation for each coordinate:

\begin{equation}

  • b \sqrt{\left(\frac{dx}{dt} \right)^2+\left(\frac{dy}{dt} \right)^2} \,\frac{dx}{dt} = m\frac{d^2x}{dt^2} \rightarrow \frac{d^2x}{dt^2} = - \frac{b}{m} \sqrt{\left(\frac{dx}{dt} \right)^2+\left(\frac{dy}{dt} \right)^2} \,\frac{dx}{dt} \end{equation}

\begin{equation} -mg - b \sqrt{\left(\frac{dx}{dt} \right)^2+\left(\frac{dy}{dt} \right)^2} \,\frac{dy}{dt} = m\frac{d^2y}{dt^2} \rightarrow \frac{d^2y}{dt^2} = - \frac{b}{m} \sqrt{\left(\frac{dx}{dt} \right)^2+\left(\frac{dy}{dt} \right)^2} \,\frac{dy}{dt} -g \end{equation}

These equations were solved numerically in this notebook.

8) Linear spring and damping on bidimensional horizontal movement

This situation is very similar to the example of horizontal movement with one spring and two masses, with a damper added in parallel to the spring.

Now, the forces acting on each mass are the force due to the spring and the force due to the damper.

So, the forces acting on mass 1 is:

\begin{equation} \vec{\bf{F_1}} = b\frac{d(\vec{\bf{r_2}}-\vec{\bf{r_1}})}{dt} + k\left(||\vec{\bf{r_2}}-\vec{\bf{r_1}}||-l_0\right)\frac{(\vec{\bf{r_2}}-\vec{\bf{r_1}})}{||\vec{\bf{r_2}}-\vec{\bf{r_1}}||} = b\frac{d(\vec{\bf{r_2}}-\vec{\bf{r_1}})}{dt} + k\left(1-\frac{l_0}{||\vec{\bf{r_2}}-\vec{\bf{r_1}}||}\right)(\vec{\bf{r_2}}-\vec{\bf{r_1}}) \end{equation}

and the forces acting on mass 2 is:

\begin{equation} \vec{\bf{F_2}} = b\frac{d(\vec{\bf{r_1}}-\vec{\bf{r_2}})}{dt} + k\left(||\vec{\bf{r_2}}-\vec{\bf{r_1}}||-l_0\right)\frac{(\vec{\bf{r_1}}-\vec{\bf{r_2}})}{||\vec{\bf{r_1}}-\vec{\bf{r_2}}||}= b\frac{d(\vec{\bf{r_1}}-\vec{\bf{r_2}})}{dt} + k\left(1-\frac{l_0}{||\vec{\bf{r_2}}-\vec{\bf{r_1}}||}\right)(\vec{\bf{r_1}}-\vec{\bf{r_2}}) \end{equation}

Applying the second Newton law for the masses:

\begin{equation} m_1\frac{d^2\vec{\bf{r_1}}}{dt^2} = b\frac{d(\vec{\bf{r_2}}-\vec{\bf{r_1}})}{dt}+k\left(1-\frac{l_0}{||\vec{\bf{r_2}}-\vec{\bf{r_1}}||}\right)(\vec{\bf{r_2}}-\vec{\bf{r_1}}) \rightarrow \frac{d^2\vec{\bf{r_1}}}{dt^2} = -\frac{b}{m_1}\frac{d\vec{\bf{r_1}}}{dt} -\frac{k}{m_1}\left(1-\frac{l_0}{||\vec{\bf{r_2}}-\vec{\bf{r_1}}||}\right)\vec{\bf{r_1}} + \frac{b}{m_1}\frac{d\vec{\bf{r_2}}}{dt}+\frac{k}{m_1}\left(1-\frac{l_0}{||\vec{\bf{r_2}}-\vec{\bf{r_1}}||}\right)\vec{\bf{r_2}} \rightarrow \frac{d^2x_1\hat{\bf{i}}}{dt^2}+\frac{d^2y_1\hat{\bf{j}}}{dt^2} = -\frac{b}{m_1}\left(\frac{dx_1\hat{\bf{i}}}{dt}+\frac{dy_1\hat{\bf{j}}}{dt}\right)-\frac{k}{m_1}\left(1-\frac{l_0}{||\vec{\bf{r_2}}-\vec{\bf{r_1}}||}\right)(x_1\hat{\bf{i}}+y_1\hat{\bf{j}})+\frac{b}{m_1}\left(\frac{dx_2\hat{\bf{i}}}{dt}+\frac{dy_2\hat{\bf{j}}}{dt}\right)+\frac{k}{m_1}\left(1-\frac{l_0}{||\vec{\bf{r_2}}-\vec{\bf{r_1}}||}\right)(x_2\hat{\bf{i}}+y_2\hat{\bf{j}}) = -\frac{b}{m_1}\left(\frac{dx_1\hat{\bf{i}}}{dt}+\frac{dy_1\hat{\bf{j}}}{dt}-\frac{dx_2\hat{\bf{i}}}{dt}-\frac{dy_2\hat{\bf{j}}}{dt}\right)-\frac{k}{m_1}\left(1-\frac{l_0}{\sqrt{(x_2-x_1)^2+(y_2-y_1)^2}}\right)(x_1\hat{\bf{i}}+y_1\hat{\bf{j}}-x_2\hat{\bf{i}}-y_2\hat{\bf{j}}) \end{equation}

\begin{equation} m_2\frac{d^2\vec{\bf{r_2}}}{dt^2} = b\frac{d(\vec{\bf{r_1}}-\vec{\bf{r_2}})}{dt}+k\left(1-\frac{l_0}{||\vec{\bf{r_2}}-\vec{\bf{r_1}}||}\right)(\vec{\bf{r_1}}-\vec{\bf{r_2}}) \rightarrow \frac{d^2\vec{\bf{r_2}}}{dt^2} = -\frac{b}{m_2}\frac{d\vec{\bf{r_2}}}{dt} -\frac{k}{m_2}\left(1-\frac{l_0}{||\vec{\bf{r_2}}-\vec{\bf{r_1}}||}\right)\vec{\bf{r_2}} + \frac{b}{m_2}\frac{d\vec{\bf{r_1}}}{dt}+\frac{k}{m_2}\left(1-\frac{l_0}{||\vec{\bf{r_2}}-\vec{\bf{r_1}}||}\right)\vec{\bf{r_1}} \rightarrow \frac{d^2x_2\hat{\bf{i}}}{dt^2}+\frac{d^2y_2\hat{\bf{j}}}{dt^2} = -\frac{b}{m_2}\left(\frac{dx_2\hat{\bf{i}}}{dt}+\frac{dy_2\hat{\bf{j}}}{dt}\right)-\frac{k}{m_2}\left(1-\frac{l_0}{||\vec{\bf{r_2}}-\vec{\bf{r_1}}||}\right)(x_2\hat{\bf{i}}+y_2\hat{\bf{j}})+\frac{b}{m_2}\left(\frac{dx_1\hat{\bf{i}}}{dt}+\frac{dy_1\hat{\bf{j}}}{dt}\right)+\frac{k}{m_2}\left(1-\frac{l_0}{||\vec{\bf{r_2}}-\vec{\bf{r_1}}||}\right)(x_1\hat{\bf{i}}+y_1\hat{\bf{j}})=-\frac{b}{m_2}\left(\frac{dx_2\hat{\bf{i}}}{dt}+\frac{dy_2\hat{\bf{j}}}{dt}-\frac{dx_1\hat{\bf{i}}}{dt}-\frac{dy_1\hat{\bf{j}}}{dt}\right)-\frac{k}{m_2}\left(1-\frac{l_0}{\sqrt{(x_2-x_1)^2+(y_2-y_1)^2}}\right)(x_2\hat{\bf{i}}+y_2\hat{\bf{j}}-x_1\hat{\bf{i}}-y_1\hat{\bf{j}}) \end{equation}

Now, we can separate the equations for each of the coordinates:

\begin{equation} \frac{d^2x_1}{dt^2} = -\frac{b}{m_1}\left(\frac{dx_1}{dt}-\frac{dx_2}{dt}\right)-\frac{k}{m_1}\left(1-\frac{l_0}{\sqrt{(x_2-x_1)^2+(y_2-y_1)^2}}\right)(x_1-x_2) \end{equation}

\begin{equation} \frac{d^2y_1}{dt^2} = -\frac{b}{m_1}\left(\frac{dy_1}{dt}-\frac{dy_2}{dt}\right)-\frac{k}{m_1}\left(1-\frac{l_0}{\sqrt{(x_2-x_1)^2+(y_2-y_1)^2}}\right)(y_1-y_2) \end{equation}

\begin{equation} \frac{d^2x_2}{dt^2} = -\frac{b}{m_2}\left(\frac{dx_2}{dt}-\frac{dx_1}{dt}\right)-\frac{k}{m_2}\left(1-\frac{l_0}{\sqrt{(x_2-x_1)^2+(y_2-y_1)^2}}\right)(x_2-x_1) \end{equation}

\begin{equation} \frac{d^2y_2}{dt^2} = -\frac{b}{m_2}\left(\frac{dy_2}{dt}-\frac{dy_1}{dt}\right)-\frac{k}{m_2}\left(1-\frac{l_0}{\sqrt{(x_2-x_1)^2+(y_2-y_1)^2}}\right)(y_2-y_1) \end{equation}

If you want to solve these equations numerically, you must break these equations into first-order equations:

\begin{equation} \frac{dv_{x_1}}{dt} = -\frac{b}{m_1}\left(v_{x_1}-v_{x_2}\right)-\frac{k}{m_1}\left(1-\frac{l_0}{\sqrt{(x_2-x_1)^2+(y_2-y_1)^2}}\right)(x_1-x_2) \end{equation}

\begin{equation} \frac{dv_{y_1}}{dt} = -\frac{b}{m_1}\left(v_{y_1}-v_{y_2}\right)-\frac{k}{m_1}\left(1-\frac{l_0}{\sqrt{(x_2-x_1)^2+(y_2-y_1)^2}}\right)(y_1-y_2) \end{equation}

\begin{equation} \frac{dv_{x_2}}{dt} = -\frac{b}{m_2}\left(v_{x_2}-v_{x_1}\right)-\frac{k}{m_2}\left(1-\frac{l_0}{\sqrt{(x_2-x_1)^2+(y_2-y_1)^2}}\right)(x_2-x_1) \end{equation}

\begin{equation} \frac{dv_{y_2}}{dt} = -\frac{b}{m_2}\left(v_{y_2}-v_{y_1}\right)-\frac{k}{m_2}\left(1-\frac{l_0}{\sqrt{(x_2-x_1)^2+(y_2-y_1)^2}}\right)(y_2-y_1) \end{equation}

\begin{equation} \frac{dx_1}{dt} = v_{x_1} \end{equation}

\begin{equation} \frac{dy_1}{dt} = v_{y_1} \end{equation}

\begin{equation} \frac{dx_2}{dt} = v_{x_2} \end{equation}

\begin{equation} \frac{dy_2}{dt} = v_{y_2} \end{equation}

To solve the equations numerically, we will use the $m_1=1 kg$, $m_2 = 2 kg$, $l_0 = 0.5 m$, $k = 10 N/m$, $b = 0.6 Ns/m$ and $x_{1_0} = 0 m$, $x_{2_0} = 0 m$, $y_{1_0} = 1 m$, $y_{2_0} = -1 m$, $v_{x1_0} = -2 m/s$, $v_{x2_0} = 1 m$, $v_{y1_0} = 0 m$, $v_{y2_0} = 0 m$.

In [5]:
x01 = 0
y01= 1
x02 = 0
y02 = -1

vx01 = -2
vy01 = 0
vx02 = 1
vy02 = 0

x1= x01
y1 = y01
x2= x02
y2 = y02
vx1= vx01
vy1 = vy01
vx2= vx02
vy2 = vy02
r1 = np.array([x1,y1])
r2 = np.array([x2,y2])

k = 10
m1 = 1
m2 = 2
b = 0.6
l0 = 0.5

dt = 0.001
t = np.arange(0,5,dt)

for i in t[1:]:
    dvx1dt = -b/m1*(vx1-vx2) -k/m1*(1-l0/np.sqrt((x2-x1)**2+(y2-y1)**2))*(x1-x2)
    dvx2dt = -b/m2*(vx2-vx1) -k/m2*(1-l0/np.sqrt((x2-x1)**2+(y2-y1)**2))*(x2-x1)
    dvy1dt = -b/m1*(vy1-vy2) -k/m1*(1-l0/np.sqrt((x2-x1)**2+(y2-y1)**2))*(y1-y2)
    dvy2dt = -b/m2*(vy2-vy1) -k/m2*(1-l0/np.sqrt((x2-x1)**2+(y2-y1)**2))*(y2-y1)
    dx1dt = vx1
    dx2dt = vx2
    dy1dt = vy1
    dy2dt = vy2
    x1 = x1 + dt*dx1dt
    x2 = x2 + dt*dx2dt
    y1 = y1 + dt*dy1dt
    y2 = y2 + dt*dy2dt
    vx1 = vx1 + dt*dvx1dt
    vx2 = vx2 + dt*dvx2dt
    vy1 = vy1 + dt*dvy1dt
    vy2 = vy2 + dt*dvy2dt
    r1 = np.vstack((r1,np.array([x1,y1])))
    r2 = np.vstack((r2,np.array([x2,y2])))

springDampLength = np.sqrt((r1[:,0]-r2[:,0])**2+(r1[:,1]-r2[:,1])**2)    
plt.figure()
plt.plot(t, springDampLength)
plt.show()

plt.figure()
plt.plot(r1[:,0], r1[:,1],'r.')
plt.plot(r2[:,0], r2[:,1],'b.')
plt.plot((m1*r1[:,0]+m2*r2[:,0])/(m1+m2), (m1*r1[:,1]+m2*r2[:,1])/(m1+m2),'g.')
plt.xlim(-2,2)
plt.ylim(-2,2)
plt.show()

9) Simple muscle model

The diagram below shows a simple muscle model. The spring in the left represents the tendinous tissues and the spring in the right represents the elastic properties of the muscle fibers. The damping is present to model the viscous properties of the muscle fibers, the element CE is the contractile element (force production) and the mass $m$ is the muscle mass.

The length $L_{MT}$ is the length of the muscle plus the tendon. In our model $L_{MT}$ is constant, but it could be a function of the joint angle.

The length of the tendon will be denoted by $l_T(t)$ and the muscle length, by $l_{m}(t)$. Both lengths are related by each other by the following expression:

\begin{equation} l_t(t) + l_m(t) = L_{MT} \end{equation}

The free-body diagram of the muscle mass is depicted below.

The resultant force being applied in the muscle mass is:

$\vec{\bf{F}} = -k_T(||\vec{\bf{r_m}}||-l_{t_0})\frac{\vec{\bf{r_m}}}{||\vec{\bf{r_m}}||} + b\frac{d(L_{MT}\hat{\bf{i}} - \vec{\bf{r_{m}}})}{dt} + k_m (||L_{MT}\hat{\bf{i}} - \vec{\bf{r_{m}}}||-l_{m_0})\frac{L_{MT}\hat{\bf{i}} - \vec{\bf{r_{m}}}}{||L_{MT}\hat{\bf{i}} - \vec{\bf{r_{m}}}||} +\vec{\bf{F}}{\bf{_{CE}}}(t)$

Since the model is unidimensional, we can assume that the force $\vec{\bf{F}}\bf{_{CE}}(t)$ is in the x direction, so the analysis will be done only in this direction.

$F = -k_T(l_t-l_{t_0}) + b\frac{d(L_{MT} - l_t)}{dt} + k_m (l_m-l_{m_0}) + F_{CE}(t) = -k_T(l_t-l_{t_0}) -b\frac{dl_t}{dt} + k_m (L_{MT}-l_t-l_{m_0}) + F_{CE}(t) = -b\frac{dl_t}{dt}-(k_T+k_m)l_t+F_{CE}(t)+k_Tl_{t_0}+k_m(L_{MT}-l_{m_0})$

Applying the second Newton law:

$m\frac{d^2l_t}{dt^2} = -b\frac{dl_t}{dt}-(k_T+k_m)l_t+F_{CE}(t)+k_Tl_{t_0}+k_m(L_{MT}-l_{m_0})$

To solve this equation, we must break the equation into two first-order differential equations:

\begin{equation} \frac{dvt}{dt} = - \frac{b}{m}v_t - \frac{k_T+k_m}{m}l_t +\frac{F_{CE}(t)}{m} + \frac{k_T}{m}l_{t_0}+\frac{k_m}{m}(L_{MT}-l_{m_0}) \end{equation}

\begin{equation} \frac{d l_t}{dt} = v_t \end{equation}

Now, we can solve these equations by using some numerical method. To obtain the solution, we will use the damping factor of the muscle as $b = 10\,Ns/m$, the muscle mass is $m = 2 kg$, the stiffness of the tendon as $k_t=1000\,N/m$ and the elastic element of the muscle as $k_m=1500\,N/m$. The tendon-length is $L_{MT} = 0.35\,m$, and the tendon equilibrium length is $l_{t0} = 0.28\,m$ and the muscle fiber equilibrium length is $l_{m0} = 0.07\,m$. Both the tendon and the muscle fiber are at their equilibrium lengths and at rest.

Also, we will model the force of the contractile element as a Heaviside step of $90\,N$ (90 N beginning at $t=0$), but normally it is modeled as a function of $l_m$ and $v_m$ having a neural activation signal as input.

In [6]:
m = 2
b = 10
km = 1500
kt = 1000
lt0 = 0.28
lm0 = 0.07

Lmt = 0.35
vt0 = 0

dt = 0.0001
t = np.arange(0, 10, dt)

Fce = 90

lt = lt0
vt = vt0
ltp = np.array([lt0])
lmp = np.array([lm0])
Ft = np.array([0])
for i in range(1,len(t)):
    dvtdt = -b/m*vt-(kt+km)/m*lt  + Fce/m + kt/m*lt0 +km/m*(Lmt-lm0)
    dltdt = vt
    vt = vt + dt*dvtdt
    lt = lt + dt*dltdt
    Ft = np.vstack((Ft,np.array(kt*(lt-lt0))))
    ltp = np.vstack((ltp,np.array(lt)))
    lmp = np.vstack((lmp,np.array(Lmt - lt)))

plt.figure()
plt.plot(t,Ft)
plt.xlabel('t')
plt.ylabel('Tendon force (N)')
plt.show()
plt.figure()
plt.plot(t,ltp)
plt.plot(t,lmp)
plt.xlabel('t')
plt.ylabel('Length (m)')
plt.show()

Problems

Solve the problems 3.3.9, 3.3.20, 10.1.6, 12.1.6(a,b,c,d,f), 12.1.7, 12.1.10 (a,b) from Ruina and Pratap book

References