Report on the Heavy Spring With Double Pendulum On T junction Demonstration

Nasser M. Abbasi, updated Sept 10, 2012

This report is a description of the model used and mathematical derivation of equations of motion of the following dynamic system. It is a spring pendulum which carries as its “bob” another spring pendulum

The top spring represents the pendulum. It is a heavy spring with mass m and attached to the end of it is massless spring which holds a bob at the end of it with mass M. Both springs are assumed infinitely rigid against rotation and can only move along their length. If it helps visualize this, we can imagine the lower spring to be wrapped around an infinitely rigid thin massless rod with the bob being a bead that slides along the rod attached to this spring.
The bob has three degrees of freedom: the angle of rotation of the pendulum, the distance of the bob along the lower rod and the current length of the top pendulum. The top spring is assumed to have mass, while the lower spring is massless. 
To derive the equations of motion (3 equations of motions are needed since there are 3 degrees of freedom) Lagrangian method is used. First the kinetic energy and then the potential energy are derived.

Before deriving the Lagrangian, we need to discuss the top spring more. This spring will have relaxed length of (Subscript[L, 0]) which is the length when the spring is laying on a horizontal surface, hence its weight has no effect. When this spring is held upright, it will now extend a little more due to its own weight. We use Rayleigh factor of (1/3) to find this additional extension. When we attach the bob, there will be an additional extension, called the static extension. It is at this state that the simulation starts. Any additional pulling of the spring at this time is what will cause the vibration.

Displacement of the top spring is measured from the static deflection position as shown in this diagram

To model the top spring, since it has mass m and a length, its kinetic energy comes from rotational motion and motion due to translation of its center of mass in the top-down direction as it vibrates. For rotation it is modeled as rigid rod of moment of inertia ((m (L^2))/3)where L is the total current length of the top spring measured from the origin where the pendulum is attached. For kinetic energy due to translational vibration motion up and down, it is modeled as point mass located at the bottom of the spring of mass m/3

For the lower spring, it is assumed to have zero relaxed length and the extension of the lower spring is measured as x starting from where the lower spring is attached to the top spring. In the above diagram, all terms are assumed to be in the positive sense of motion.
Given the above diagram, we can now find the kinetic and potential energy of the system and derive the equations of motion.

Derivations of the equations of motion

KE

bobHasMass=False;
If[bobHasMass,
L0=len0+((m0 g)/(3 k1))+((M0 g)/k1),
L0=len0+((m0 g)/(3 k1))
];

This is the position vector of the bob mass as measured in the inertial frame of reference. The origin of this frame of reference is where the top pendulum is attached. The position x direction is towards the right, and positive y direction is downwards.

r[t]={(L0+y[t])Sin[θ[t]]+x[t]Cos[θ[t]],(L0+y[t])Cos[θ[t]]-x[t] Sin[θ[t]]};
bobVelocity=D[r[t],t]

This is speed of the bob after squaring it. It comes from the above expression.

bobVelocitySquared=((x'[t]+(L0+y[t]) θ'[t])^2+(x[t] θ'[t]-y'[t])^2);

If[bobHasMass,
ke=(1/2)(M0) bobVelocitySquared+(1/2)m0 (((L0+y[t])^2)/3) θ'[t]^2+(1/2)(m0/3)y'[t]^2
,
ke=(1/2)m0 (L0+y[t])^2/3 θ'[t]^2+(1/2)(m0/3)y'[t]^2
]

PE

Assume PE is zero at the bottom of the pendulum when it hangs in the vertical position.

If[bobHasMass,
pe= m0 g ((L0+y[t])/2)(1-Cos[θ[t]])+ M0 g((L0+y[t]) (1-Cos[θ[t]])+x[t] Sin[θ[t]])+(1/2) k x[t]^2+(1/2) k1 (y[t])^2
,
pe= m0 g ((L0+y[t])/2)(1-Cos[θ[t]])+(1/2) k1 (y[t])^2

]

Lagrangian

lagrangian=ke-pe

Equations of motion

x equation of motion

If[bobHasMass,
eq1=D[D[lagrangian,x'[t]],t]-D[lagrangian,x[t]]==0
]

y equation of motion

eq2=D[D[lagrangian,θ'[t]],t]-D[lagrangian,θ[t]]==0

θ equation of motion

eq3=D[D[lagrangian,y'[t]],t]-D[lagrangian,y[t]]==0

Use NDSolve to solve the equations of motion using typical values for parameters

If[bobHasMass,
params={k->500,k1->1000,M0->10,g->9.8,len0->0.5,m0->5},
params={k1->1000,g->9.8,len0->0.5,m0->5},
]
If[bobHasMass,
ic={x[0]==0.35,x'[0]==0,θ[0]==Pi/8,θ'[0]==0,y[0]==0,y'[0]==0}
,
ic={θ[0]==0,θ'[0]==0.1,y[0]==0.1,y'[0]==0}
];

If[bobHasMass,
sol=First@NDSolve[Flatten[{eq1,eq2,eq3,ic}]//.params,{x,y,θ,x',y',θ'},{t,0,9},Method->{"BDF"}]
,
sol=First@NDSolve[Flatten[{eq2,eq3,ic}]//.params,{y,θ,y',θ'},{t,0,9},Method->{"BDF"}]
]

Plot the solution

data=Table[(pe+ke)//.params/.sol,{t,0,9,.01}];
ListPlot[data,Joined->True,PlotRange->All]

Plot[180/Pi θ[t]/.sol,{t,0,9},PlotRange->All,AxesOrigin->{0,0}]

Plot[y[t]/.sol,{t,0,9},PlotRange->All,AxesOrigin->{0,0}]

Plot[x[t]/.sol,{t,0,9},PlotRange->All,AxesOrigin->{0,0}]