Three pivoted beams coupled by springs

A three-degree-of-freedom oscillator: three massless beams, each pivoted at one end and carrying a point mass at the other, linked by a vertical chain of four springs. The equations of motion come from Lagrange's equations and the small-angle approximation, then from a Simulink model you can run below.

Problem statement

Derive the equations of motion of the system shown in the figure. The beams have negligible weight and the motion is restricted to the vertical plane.

Each of the three beams has length $2l$. It is pinned to a wall at one end and carries a point mass $m$ at the other. The lower and upper beams are pinned on the left, the middle beam on the right. A vertical chain of four identical springs of stiffness $k$ runs from the floor to the ceiling. It is attached to the midpoint of each beam, at distance $l$ from the pivot.

The system is conservative: there is no friction and no external input. The right-hand sides of all Lagrange equations are therefore zero, and the system performs undamped oscillations.

Figure 1. Beams 1 (bottom), 2 (middle) and 3 (top), each of length 2l with a point mass m at the free end, connected by four springs k at their midpoints

Generalized coordinates

Each beam can only rotate about its pivot, so the system has three degrees of freedom. The generalized coordinates are the rotation angles $\theta_1,\theta_2,\theta_3$ of the beams from the horizontal position. The positive directions are chosen so that a positive angle always lifts the mass:

$\theta_1$ and $\theta_3$ (beams pivoted on the left) are positive counter-clockwise,
$\theta_2$ (beam pivoted on the right) is positive clockwise.

With this convention, the rise of each mass, its speed, and the vertical displacement of each beam midpoint are:

$$\Delta h_i = 2l\sin\theta_i,\qquad v_i = \frac{d}{dt}(2l\theta_i) = 2l\dot\theta_i,\qquad y_{\text{mid},i} = l\sin\theta_i$$

Spring elongations: exact and approximate

Exact geometry gives long expressions for the spring elongations. The springs tilt as the attachment points move along arcs. For example, for the bottom and the second spring (with $e_{10}, e_{20}$ the initial lengths):

$$x_1 = \sqrt{2l^2(1-\cos\theta_1)^2 + e_{10}^2 - 2\sqrt2\,l\sqrt{1-\cos\theta_1}\,e_{10}\cos\tfrac{\theta_1}{2}},\qquad x_2 = \sqrt{\big(2l - l(\cos\theta_1+\cos\theta_2)\big)^2 + \big(e_{20} - l(\sin\theta_1-\sin\theta_2)\big)^2} - e_{20}$$

Using these directly would lead to very complicated equations that are probably impossible to solve analytically. The paper therefore assumes that springs deform only vertically, by the vertical displacement of their end points:

$$x_1 = l\sin\theta_1,\quad x_2 = l(\sin\theta_2-\sin\theta_1),\quad x_3 = l(\sin\theta_3-\sin\theta_2),\quad x_4 = l\sin\theta_3$$

This is valid only for small angles, $|\theta_i| \lesssim 5^\circ$ (≈ 0.087 rad). The bottom spring is stretched when beam 1 rises, and the top spring is compressed when beam 3 rises. Only $x_j^2$ enters the energy, so the sign does not matter.

Energies of the system

TerminologyAs in the course: attractive energy $\mathcal{A}$ = potential energy, repulsive coenergy $\mathcal{R}^*$ = kinetic coenergy, and $\mathcal{L} = \mathcal{R}^*-\mathcal{A}$.

3.1 Repulsive (kinetic) coenergy

The beams are massless, so only the three point masses contribute. Each moves on a circle of radius $2l$:

$$\mathcal{R}^* = \frac12 m\sum_{i=1}^3 v_i^2 = \frac12 m\left(4l^2\dot\theta_1^2+4l^2\dot\theta_2^2+4l^2\dot\theta_3^2\right) = 2ml^2\left(\dot\theta_1^2+\dot\theta_2^2+\dot\theta_3^2\right)$$

Typo in the paper: the paper writes the result as $ml^2(\dot\theta_1^2+\dot\theta_2^2+\dot\theta_3^2)$, dropping the factor 2. The partial derivatives that follow ($\partial\mathcal{L}/\partial\dot\theta_i = 4ml^2\dot\theta_i$) use the correct value $2ml^2$, so the final equations are not affected.

3.2 Attractive (potential) energy

$$\mathcal{A} = mg\sum_{i=1}^3\Delta h_i + \frac12 k\sum_{j=1}^4 x_j^2 = 2mgl(\sin\theta_1+\sin\theta_2+\sin\theta_3) + \frac12 kl^2\left[\sin^2\theta_1 + (\sin\theta_2-\sin\theta_1)^2 + (\sin\theta_3-\sin\theta_2)^2 + \sin^2\theta_3\right]$$

For small angles $\sin\theta_i\approx\theta_i$ and $\cos\theta_i\approx1$:

$$\mathcal{A} = 2mgl(\theta_1+\theta_2+\theta_3) + \frac12 kl^2\left[\theta_1^2 + (\theta_2-\theta_1)^2 + (\theta_3-\theta_2)^2 + \theta_3^2\right]$$

3.3 Lagrangian

$$\mathcal{L} = 2ml^2(\dot\theta_1^2+\dot\theta_2^2+\dot\theta_3^2) - 2mgl(\theta_1+\theta_2+\theta_3) - \frac12 kl^2\left[\theta_1^2 + (\theta_2-\theta_1)^2 + (\theta_3-\theta_2)^2 + \theta_3^2\right]$$

Lagrange equations

$$\frac{d}{dt}\frac{\partial\mathcal{L}}{\partial\dot\theta_i} - \frac{\partial\mathcal{L}}{\partial\theta_i} = 0,\qquad i=1,2,3$$

$$\begin{aligned}\frac{\partial\mathcal{L}}{\partial\dot\theta_i} &= 4ml^2\dot\theta_i\\ \frac{d}{dt}\frac{\partial\mathcal{L}}{\partial\dot\theta_i} &= 4ml^2\ddot\theta_i\end{aligned}$$

$$\begin{aligned}\frac{\partial\mathcal{L}}{\partial\theta_1} &= -2mgl - kl^2\theta_1 + kl^2(\theta_2-\theta_1)\\ \frac{\partial\mathcal{L}}{\partial\theta_2} &= -2mgl - kl^2(\theta_2-\theta_1) + kl^2(\theta_3-\theta_2)\\ \frac{\partial\mathcal{L}}{\partial\theta_3} &= -2mgl - kl^2(\theta_3-\theta_2) - kl^2\theta_3\end{aligned}$$

Equations of motion:

$$\begin{aligned}4ml^2\ddot\theta_1 + 2mgl + kl^2\theta_1 - kl^2(\theta_2-\theta_1) &= 0\\ 4ml^2\ddot\theta_2 + 2mgl + kl^2(\theta_2-\theta_1) - kl^2(\theta_3-\theta_2) &= 0\\ 4ml^2\ddot\theta_3 + 2mgl + kl^2\theta_3 + kl^2(\theta_3-\theta_2) &= 0\end{aligned}$$

Dividing by $4ml^2$ and introducing

$$c = \frac{g}{2l},\qquad c_1 = \frac{k}{4m}$$

gives the compact form used for Simulink:

$$\begin{aligned}\ddot\theta_1 &= -c - 2c_1\theta_1 + c_1\theta_2\\ \ddot\theta_2 &= -c - 2c_1\theta_2 + c_1\theta_1 + c_1\theta_3\\ \ddot\theta_3 &= -c - 2c_1\theta_3 + c_1\theta_2\end{aligned}$$

These are three coupled linear, non-homogeneous differential equations. The constant term $-c$ comes from gravity. In matrix form:

$$\ddot{\boldsymbol\theta} + c_1\mathbf{K}\boldsymbol\theta = -c\,\mathbf{1},\qquad \mathbf{K} = \begin{bmatrix}2&-1&0\\-1&2&-1\\0&-1&2\end{bmatrix},\quad \mathbf{1}=\begin{bmatrix}1\\1\\1\end{bmatrix}$$

Simulink model

Each angle gets a chain of two integrators ($\ddot\theta_i\to\dot\theta_i\to\theta_i$). The summing junction of each chain adds the constant $-c$, the self term $-2c_1\theta_i$, and the coupling terms $c_1\theta_j$ from the neighbouring beams. Beam 2 is coupled to both 1 and 3; beams 1 and 3 are coupled only to beam 2.

Figure 2. Structure of the Simulink model (simplified): three integrator chains, gain blocks −2c₁ (self) and +c₁ (coupling to neighbours, dashed), constant blocks −c from gravity.

The paper uses reference parameters chosen so that the motion stays inside the small-angle limit:

Parameter	Value	Derived
$m$	0.1 kg	$c = g/2l = 16.35\ \text{s}^{-2}$
$l$	0.3 m	$c_1 = k/4m = 750\ \text{s}^{-2}$
$k$	300 N/m	(the paper writes the unit as "Nm", which should be N/m)

Each parameter is then varied by ±50 % to study its influence.

Simulation results from the paper

The buttons load the same conditions into the simulator. All simulations reproduce the paper's plots.

Zero initial conditions

All angles and velocities are zero, yet the system moves. The equations are non-homogeneous: in the horizontal position the springs are unstressed and gravity is not balanced. $\theta_1=\theta_3$ oscillate between 0 and ≈ −0.065 rad, and $\theta_2$ between 0 and ≈ −0.09 rad, at about 3.3 Hz.

Non-zero initial conditions

$\theta_{01}=5^\circ$, $\theta_{02}=3^\circ$, $\theta_{03}=1^\circ$. The motion is complex-periodic (a sum of several harmonics) and undamped. It is centred around the shifted equilibrium, not around zero: the angles range from about +0.08 to −0.16 rad.

Parameter study (θ₁, zero IC)

k ↑: amplitude ↓ (∝ 1/k), frequency ↑. m ↑: amplitude ↑ (∝ m), frequency ↓. l ↑: amplitude ↓ (∝ 1/l), frequency unchanged. Values used: k = 150/300/450, m = 0.05/0.10/0.15, l = 0.15/0.30/0.45.

SymmetryWith zero initial conditions (or any with $\theta_{01}=\theta_{03}$) the responses $\theta_1$ and $\theta_3$ are identical, because the system is symmetric about the middle beam. The paper uses this to focus on $\theta_1$ only in the parameter study.

Second approach: horizontal position as equilibrium

The paper then repeats the derivation assuming that the horizontal position is the equilibrium. In that position the springs are pre-deformed by $e_{01},e_{02},e_{03},e_{04}$. The potential energy becomes

$$\mathcal{A} = 2mgl(\theta_1+\theta_2+\theta_3) + \frac12 k\left[(l\theta_1+e_{01})^2 + (l(\theta_2-\theta_1)+e_{02})^2 + (l(\theta_3-\theta_2)+e_{03})^2 + (l\theta_3+e_{04})^2\right]$$

and the equations of motion acquire constant terms from the pre-deformations:

$$\begin{aligned}4ml^2\ddot\theta_1 + 2mgl + kl^2\theta_1 - kl^2(\theta_2-\theta_1) + kl(e_{01}-e_{02}) &= 0\\ 4ml^2\ddot\theta_2 + 2mgl + kl^2(\theta_2-\theta_1) - kl^2(\theta_3-\theta_2) + kl(e_{02}-e_{03}) &= 0\\ 4ml^2\ddot\theta_3 + 2mgl + kl^2\theta_3 + kl^2(\theta_3-\theta_2) + kl(e_{03}+e_{04}) &= 0\end{aligned}$$

For $\theta=0$ to be an equilibrium, the constant terms must cancel: the spring forces must balance gravity.

$$2mgl + kl(e_{01}-e_{02}) = 0,\qquad 2mgl + kl(e_{02}-e_{03}) = 0,\qquad 2mgl + kl(e_{03}+e_{04}) = 0$$

This gives three equations for four unknowns, so one pre-deformation can be chosen freely. With $e_{04}=0$:

$$e_{01} = -\frac{6mg}{k},\qquad e_{02} = -\frac{4mg}{k},\qquad e_{03} = -\frac{2mg}{k}$$

The negative signs mean compression: the lower springs carry the weight of the masses above them. The equations become homogeneous:

$$\ddot\theta_1 + 2c_1\theta_1 - c_1\theta_2 = 0,\qquad \ddot\theta_2 + 2c_1\theta_2 - c_1\theta_1 - c_1\theta_3 = 0,\qquad \ddot\theta_3 + 2c_1\theta_3 - c_1\theta_2 = 0$$

The Simulink model is the same, but without the constant blocks $-c$. With $\theta_{01}=5^\circ,\theta_{02}=3^\circ,\theta_{03}=1^\circ$ the oscillations are now centred on zero.

Review and normal-mode analysis

The derivation is correctThe elongations, energies, partial derivatives and all equations of motion were checked, for both approaches. A numerical simulation reproduces the paper's plots: for zero initial conditions, θ₁ = θ₃ reach −0.065 rad and θ₂ reaches −0.090 rad. There are two small slips. The factor 2 is missing in the written $\mathcal{R}^*$ (section 3), but the derivatives use the correct value. The unit of $k$ is given as "Nm" instead of N/m.

8.1 The two approaches describe the same dynamics

The paper calls the two approaches "alternative solutions based on different assumptions". It is worth making precise how they relate. The first model has a constant gravity term, so its equilibrium is not at zero:

$$c_1\mathbf{K}\boldsymbol\theta_{\text{eq}} = -c\,\mathbf{1}\;\Rightarrow\;\boldsymbol\theta_{\text{eq}} = -\frac{c}{c_1}\mathbf{K}^{-1}\mathbf{1} = -\frac{c}{c_1}\begin{bmatrix}1.5\\2\\1.5\end{bmatrix} = -\frac{2mg}{kl}\begin{bmatrix}1.5\\2\\1.5\end{bmatrix}$$

For the reference parameters this is $\boldsymbol\theta_{\text{eq}} = [-0.0327,\,-0.0436,\,-0.0327]$ rad (≈ −1.9°, −2.5°, −1.9°). Substituting $\boldsymbol\theta = \boldsymbol\theta_{\text{eq}} + \boldsymbol\psi$ gives exactly the homogeneous equations of the second approach for $\boldsymbol\psi$. In a linear system, gravity only shifts the equilibrium. Frequencies, mode shapes and the response to a disturbance are identical. Which reference is physically right depends only on whether the real springs are pre-tensioned. This also explains why the paper's parameter study holds for both approaches.

8.2 Normal modes

The eigenvalues of $\mathbf{K}$ are $2-\sqrt2,\;2,\;2+\sqrt2$, so the natural frequencies are

$$\omega_1 = \sqrt{(2-\sqrt2)\,c_1},\quad \omega_2 = \sqrt{2c_1},\quad \omega_3 = \sqrt{(2+\sqrt2)\,c_1},\qquad c_1 = \frac{k}{4m}$$

Mode	Shape [θ₁, θ₂, θ₃]	Description	f (reference parameters)
1	[1, √2, 1]	all beams move together, middle one most	3.34 Hz
2	[1, 0, −1]	outer beams in opposition, middle beam at rest	6.16 Hz
3	[1, −√2, 1]	middle beam opposes the outer ones	8.05 Hz

The frequencies are not integer multiples of each other, so a general response is complex-periodic, as the paper observes. With zero initial conditions, the gravity load $\mathbf{1}$ is symmetric and excites only modes 1 and 3. Mode 1 carries about 97 % of the amplitude, which is why that response looks almost like a single sinusoid at 3.3 Hz.

8.3 Why the parameter effects come out exactly as observed

Frequency depends only on $c_1 = k/4m$: it grows with $\sqrt{k}$, falls with $\sqrt{m}$, and does not depend on $l$ or $g$. Both the mass and the spring lever arms scale with $l$, and the beams are massless, so $l$ cancels out.
Amplitude for zero initial conditions equals twice the static deflection $|\boldsymbol\theta_{\text{eq}}| \propto 2mg/(kl)$. It is exactly proportional to $m$ and inversely proportional to $k$ and $l$.

8.4 Validity of the small-angle model

The simulator also solves a model that keeps $\sin\theta_i$ and $\cos\theta_i$ (still with vertical springs). Comparing the two shows where the linear model starts to deviate. For the reference parameters, θ₂ already reaches ≈ 5.1° with zero initial conditions, right at the limit. With large initial angles (10–20°) the difference becomes clearly visible. Note also that the paper's own non-zero case (5°, 3°, 1° in the first approach) reaches about 9.4° (0.164 rad), because the oscillation is centred on the shifted equilibrium, not on zero. That is beyond the 5° limit the paper sets for its approximation. In the second approach the same initial angles stay within ±5°.

Interactive simulator

The equations are solved numerically (4th-order Runge–Kutta, 0.5 ms step). Initial angles are entered in degrees; plots show radians, as in the paper. The animation can exaggerate the angles so the small motion is visible.

Cases

Reference

Also show the model with sin θ (no small-angle approx.)

System parameters

m [kg]

l (half beam) [m]

k [N/m]

g [m/s²]

Initial conditions

θ₁(0) [°]

θ₂(0) [°]

θ₃(0) [°]

Duration [s]

Initial angular velocities are zero.

Parameter study

Vary

θ₁(t) for 50 %, 100 % and 150 % of the current value, as in the paper.

Beam angles θ₁, θ₂, θ₃

Modal coordinates q₁, q₂, q₃

Deviation from equilibrium projected on the three mode shapes. Each one is a pure sinusoid at its own frequency.

Parameter study: θ₁(t)

Energy balance

$\mathcal{R}^*+\mathcal{A}$ must stay constant: the system is conservative.

Animation

Results table (every 0.05 s)

Review questions

Why does the system move even with all initial conditions equal to zero?

In the first approach the springs are unstressed in the horizontal position, so nothing balances gravity there. The horizontal position is not an equilibrium, and the equations are non-homogeneous (constant term $-c$). The system starts falling toward $\boldsymbol\theta_{\text{eq}}$ and oscillates around it.

What assumption makes the equations linear, and when does it fail?

Two assumptions: the springs deform only vertically, and $\sin\theta\approx\theta$, $\cos\theta\approx1$. Both require small angles (roughly below 5°). For larger angles the spring tilt and the sine nonlinearity become noticeable.

Why is the natural frequency independent of the beam length $l$?

The beams are massless. The inertia ($4ml^2$) and the spring stiffness about the pivot ($kl^2$) both scale with $l^2$, so their ratio $k/4m$ does not depend on $l$.

Why are $\theta_1$ and $\theta_3$ identical for zero initial conditions?

The system and the gravity load are symmetric about the middle beam, so the antisymmetric mode 2 ([1, 0, −1]) is never excited. Without it, $\theta_1 = \theta_3$ at all times.

How would you choose initial conditions that excite only one mode?

Start from rest with the deviation from equilibrium proportional to that mode shape, e.g. [1, √2, 1]·α for mode 1. The system then oscillates at a single frequency. Try the "Mode 1/2/3" presets.

Which other generalized coordinates could be used?

Under the vertical-spring approximation, the vertical displacements of the spring attachment points ($l\theta_i$) work equally well. The dynamics are the same; only the amplitudes of these displacements no longer depend on $l$.

Parameter	Value	Derived
\(m\)	0.1 kg	\(c = g/2l = 16.35\ \text{s}^{-2}\)
\(l\)	0.3 m	\(c_1 = k/4m = 750\ \text{s}^{-2}\)
\(k\)	300 N/m	(the paper writes the unit as "Nm", which should be N/m)

Explanation of Method