2 Lagrangian Mechanics

2.1        Overview

Joseph-Louis Lagrange (1736–1813)

In general, it is easier to perform engineering/technical calculations using a scalar quantity rather than a tensor/vector type quantity, mainly because a vector’s components depend on the selected coordinates system, and hence, more quantities to deal with. This was the main motivation for Joseph-Louis Lagrange (1736–1813), [7], [8] to start looking into the Newtonian mechanics close to a century after Newton developed his laws. Consequently, Lagrange developed a new formulation, so-called Lagrangian mechanics (1788).
Lagrange’s approach has advantages over that of Newton’s, specifically for analyzing complex multi-domain, multi-component systems. Lagrange’s approach releases us from having to consider a single inertia coordinates system and inter-component constraint forces. In addition, Langrangian method is faster and more efficient in terms of computation time and effort required to analyze and model engineering systems.

In Newtonian mechanics, a local condition, e.g., initial position and velocity (or momentum), is required for calculating the future states of a system. Using Newton’s law of motion, for a system or components of a system, the sum of forces (both applied, \vv{F}_a and constrained/internal, \vv{F}_c), is equal to the time rate of change of the momentum, \vv{p}.

(2.1)   \begin{equation*} \sum(\vv{F}_a + \vv{F}_c) = \frac{d\vv{p}}{dt} \end{equation*}

In order to identify the constraints, we usually isolate the components one by one from the rest of the system, while keeping the related dynamical equilibrium intact. This operation gives us the free-body diagram of each desired component, useful for analyzing the system’s motion dynamics and calculating inter-component constraint forces. However, in the Lagrangian approach, we consider a quantity that is like energy in dimension, the Lagrangian L, and use a set of partial differential equations (PDEs)—Euler-Lagrange or Lagrange’s equations— to analyze the system dynamics.

The latter is much more effective approach for analyzing the systems with many degrees of freedom and for dealing with multi-domain systems. In general, L is a function of coordinates considered and their time derivatives and, as well, could explicitly depend on time. For example, in a one-dimensional system, with designated coordinate x, the Lagrangian is written as L=L(x,\dot x,t). We can visualize L as the topography of a surface represented by L as a function of x and \dot x, as shown in Figure 2-1. This surface can vary with time, hence explicit dependence of L on time, or it could be stationary. An example of the former is the motion of a mass particle on the surface of a moving sphere. Similarly, the Lagrangian of such a system is stationary if the sphere is not moving. The visualization presented in reference [9] may help readers with understanding Lagrangian surface.

Figure 2-1 Lagrangian surface visualized in x-\dot x space

The foundation of Lagrangian mechanics rests on the principle of stationary action integral (also referred to as Hamilton’s principle) . This principle simply states that a system’s motion from a given state to another is such that a specific quantity (i.e., the system’s Lagrangian function) related to its motion is extremized (i.e., minimized or maximized); hence, the value of its integral (i.e., the action integral, \mathcal{A}) remains invariant [10].

The motion of a system from t_1 to t_2 is such that the action integral has a stationary value for the actual path of the motion

In other words, among all possible paths available for the motion of the system to go through, there exists one specific path that minimizes/maximizes (for most systems minimizes; hence, this is also referred to a principle of least action) the integral of the corresponding Lagrangian with respect to time. Mathematically, the stationary action integral can be stated as

(2.2)   \begin{equation*} \delta \mathcal{A} = \delta \left[ \int_{t_1}^{t_2} L(x,\dot x, t)dt \right] = 0 \end{equation*}

Using calculus of variations [11], [12], [13] and Equation (2.2) it can be shown (see section 2.5) that L should satisfy Lagrange’s equation, or

(2.3)   \begin{equation*} \frac{d}{dt} \left(\frac{\partial L}{\partial \dot x} \right) - \frac{\partial L}{\partial x} = 0 \end{equation*}

where L is defined as L=T-V, with T being the kinetic energy and V the potential energy functions. With reference to Figure 2-1, \dfrac{\partial L}{\partial x} \big|_{\dot x=const} is the slope at a selected point on the curve at the cross-section of surface L and a plane parallel to x-plane at desired \dot x, and \dfrac{d}{dt} \left(\dfrac{\partial L}{\partial \dot x}\right) \big|_{x=const} is the rate of change in the slope at the same selected point on the curve at the cross-section of a plane parallel to \dot x-plane drawn from and including the selected point the same point. In other words, we draw two planes parallel to the x and \dot x planes and equate their corresponding slopes at their intersectional point. Therefore, for a stationary point, these two quantities should be equal, as given by Euler’s equation (2.3). This is shown in the following sketch, see Figure 2-2.

Figure 2 2 A sketch for visualizing Euler-Lagrange’s equation

By working out a simple example, we show that the Lagrangian approach is equivalent to the Newtonian approach in terms of the system’s equation of motion.

2.2        Example: A Mass-Spring System

For this example, we show that Equation (2.3) gives the same results as that of Newton’s law of motion when applied to a simple mass-spring system, as sketched in Figure 2-3.

Figure 2-3 A frictionless mass-spring system

The kinetic energy for the mass m is T = \dfrac{1}{2} m\dot x^2 and the spring potential energy (i.e. stored elastic energy) with the spring constant k is V=\int kxdx = \dfrac{1}{2} kx^2. Therefore, using Equation (2.3), we get \dfrac{d}{dt} \left[ \dfrac{\partial}{\partial \dot x} \left(\dfrac{1}{2}m\dot x^2 - \dfrac{1}{2}kx^2 \right) \right] - \dfrac{\partial}{\partial x} \left(\dfrac{1}{2}m \dot x^2 - \dfrac{1}{2}kx^2 \right) = 0, or m\ddot x + kx = 0. Note that for this analysis we did not need to consider the free-body diagram of mass m nor the spring force as the constraining force acting on it; rather, we used the scalar quantity (T-V). However, the assumption of having a potential function V from which we can calculate the spring force is required (i.e., -\nabla V = - \dfrac{d(\dfrac{1}{2} kx^2 )}{dx} = -kx), see section 2.7.

In the following sections we expand on the Lagrangian method for discrete systems with related derivation, constraints and definitions for generalized coordinates, forces, and momenta.

2.3        Lagrange’s Equations for a Mass System in 3D Space

We consider a particle with mass m in a 3D space x_i = (x_1,x_2,x_3) \equiv (x,y,z), Cartesian system. By definition, the Lagrangian function is written as L=T-V= \dfrac{1}{2} m(\dot x^2 + \dot y^2 + \dot z^2 ) - V(x,y,z). We have assumed that the potential energy function is only a function of the space coordinates, so-called holonomic system. We now form two sets of derivatives \dfrac{\partial L}{\partial \dot x_i} = p_i and \dfrac {\partial L}{\partial x_i} = F_i of the Lagrangian function L=L( x,y,z,\dot x, \dot y, \dot z). Therefore, e.g., in 1D space, we have \frac{\partial L}{\partial \dot x} = m\dot x and \dfrac{\partial L}{\partial x} = - \dfrac{\partial V}{\partial x} = F_x. Hence F_x is a conservative force (see section 2.7). Now, using Newton’s second law, we can write the equation of motion, its x-component, as m\ddot x = F_x or \dfrac{d}{dt} \left(\dfrac{\partial L}{\partial \dot x} \right) = p_i = m \ddot x and \dfrac{\partial L}{\partial x}  = F_x. Therefore, \dfrac{d}{dt} \left(\dfrac{\partial L}{\partial \dot x} \right) = \dfrac{\partial L}{\partial x}. Similar derivation can be performed for y and z components of the equation of motion. Therefore, we get the Euler-Lagrange equations

(2.4)   \begin{equation*} \frac{d}{dt} \left(\frac{\partial L}{\partial \dot x_i} \right) - \frac{\partial L}{\partial x_i} = 0 \quad, i=1,2,3 \end{equation*}

The motion of the particle could be considered, in principle, in another coordinate system, e.g., a cylindrical or spherical system, as well. Therefore, we can define a set of coordinates q'_i = (q'_1,q'_2,q'_3) to represent arbitrary coordinate systems, including Cartesian or curvilinear, and write Equation (2.4) in terms of q'_i, as well, for generality.

2.4        Generalized Coordinates, Momenta, and Forces

As mentioned previously, one of the advantages of Lagrangian method is that we do not require consideration of the constrained forces. Therefore, we can include only those coordinates that correspond to the degrees of freedom related to a system. This consideration leads us to the concept of generalized coordinates, which is used in Lagrangian mechanics instead of inertia coordinates used in the Newtonian mechanics.

We now define the generalized coordinates. First, we expand the system discussed in section 2.3 to include N number of particles that move in 3N coordinate space, or q'_i=(q'_1,q'_2,\dots,q'_{3N}). However, in a real-world system we can have restrictions imposed on the system’s motion; hence, some of the coordinates are constrained and do not vary independently. For example, a particle moving in a plane (x-y) is constrained to move in z-direction (z=0). Or, the mass bob of a pendulum moving in (r-\theta) plane is restricted to move out of z-plane and if the pendulum rod has a fixed length, then only coordinate \theta varies during its motion. To capture these constraints, it is common and convenient to define generalized coordinates. Assume that for a 3N coordinates system we have N_c number of constraints. Therefore, the number of independent coordinates defining the motion is n=3N-N_c. By definition, for holonomic systems this is equal to the number of degrees of freedom [13]. Now we define the generalized coordinates as a subset of the original coordinates, with q_i=(q_1,q_2,\dots,q_n). Note that n<3N is the number of degrees of freedom which is equal to the number of generalized coordinates, and coordinates of the q_i system are not necessarily the same as those of the q'_i, by one-to-one comparison.

For derivation of the equations of motion of a system, using Lagrangian approach, we can calculate n number of equations for the system, one by one, related to each generalized coordinate. We can also use the generalized coordinates to define the velocity-phase space, as the combined set of generalized coordinates and their corresponding time derivatives. Therefore, the Lagrangian, as a functional, reads

(2.5)   \begin{equation*} L = L (\underbrace{q_1,q_2,\dots,q_n,\dot q_1,\dot q_2,\dots, \dot q_n}_{\text{phase space}} ,t) \end{equation*}

Note that the time dependence of Lagrangian may be explicit for some systems and implicit for others and that the phase-space coordinates do not necessarily have the same units/dimensions. For example, q_1 could be a displacement and q_2 an angle for a system like a pendulum with moving pivot point.

The fact that we can neglect the constrained coordinates in Lagrangian formulation is an advantage of this method over Newton’s because we don’t need to calculate the constrained “forces” in order to derive the equations of motion. Of course, the constrained forces can be calculated, if required, after having the solution to the system’s equations of motion.

Like the generalized coordinates, we also define associated generalized momenta and forces. As mentioned in the previous section, the definition of momentum in Lagrangian mechanics is more general than that of mass times the velocity. For example, it could be angular momentum, instead. Similarly, the definition of forces is not limited to mechanical forces; it can be applied, e.g., to voltage and temperature in electrical and thermal domains. Therefore, for each generalized coordinate we can define the corresponding generalized momentum and force. As given by Equation (2.6), we can write the generalized momenta and generalized force in terms of L, as

(2.6)   \begin{equation*} \begin{cases} p_i  = \dfrac{\partial L}{\partial \dot q_i}, \text{generalized momenta}\\ F_i = \dfrac{\partial L}{\partial q_i}, \text{generalized forces} \end{cases} \quad i = 1,\dots,n \end{equation*}

Section 2.7 discusses the topic of generalized force in terms of its types: conservative and non-conservative.

2.5        Hamilton’s Principle and Lagrange’s Equations

William Rowan Hamilton (1805–1865)

Hamilton’s principle, as given by Equation (2.2), is basically a mathematical expression of calculus of variations application for a system dynamical motion with the realization that Lagrangian functional is the function that should be extremized [12]. Therefore, Lagrange’s equations are resulted from the related calculations, naturally. This realization was first expressed by William Rowan Hamilton (1805-1865), [14], [11], [15].
Equation (2.4) can be written in terms of generalized coordinates, as

(2.7)   \begin{equation*} \frac{d}{dt} \left(\frac{\partial L}{\partial \dot q_i} \right) - \frac{\partial L}{\partial q_i} = 0 \quad, i = 1,2,\dots,n \end{equation*}

Equation (2.7) shows that Lagrange’s equation is consequence of, and necessary for, making the action integral stationary. We assume that variation \delta L results from variation in one of the arbitrarily selected coordinates, \delta q (dropping the subscript index for simplicity without losing the generality) while satisfying the fixed boundary conditions, or \delta q(t_1 ) = \delta q(t_2 ) = 0. Obviously, the same operation can be performed for all coordinates involved, i=1,\dots,n.

Figure 2-4 Sketch for variation of L for an arbitrary \delta q

Substituting Equation (2.7) into Equation (2.2), after dropping the subscript index and assuming L=L(q,\dot q) for simplicity, we get

    \begin{equation*} \delta \mathcal{A} = \delta \Big\{ \int_{t_1}^{t_2} \left[ \frac{d}{dt} \left( \frac{\partial L}{\partial \dot q} \right) - \frac{\partial L}{\partial q} \right] dt\Big\} = \Bigg\{ \int_{t_1}^{t_2} \delta \big\[\underbrace{\left[ \frac{d}{dt} \left( \frac{\partial L}{\partial \dot q} \right) - \frac{\partial L}{\partial q} \right]}_{L} \big\  dt \Bigg\} = 0 \end{equation*}

But \delta L = \dfrac{\partial L}{\partial q} \delta q + \dfrac{\partial L}{\partial \dot q} \delta \dot q and the last term can be written as \dfrac{\partial L}{\partial \dot q} \delta \dot q = \dfrac{d}{dt} \left( \dfrac{\partial L}{\partial \dot q} \delta q \right) - \dfrac{d}{dt} \left(\dfrac{\partial L}{\partial \dot q} \right) \delta q and hence, \delta L = \left[ \dfrac{\partial L}{\partial q} - \dfrac{d}{dt} \left(\dfrac{\partial L}{\partial \dot q} \right) \right] \delta q + \dfrac{d}{dt} \left( \dfrac{\partial L}{\partial \dot q} \delta q \right).

Back substituting into action integral expression, we get \delta \mathcal{A} = \int_{t_1}^{t_2} \left[ \dfrac{\partial L}{\partial q} - \dfrac{d}{dt} \left( \dfrac{\partial L}{\partial \dot q} \right) \right] \delta q \: dt + \int_{t_1}^{t_2} \dfrac{d}{dt} \left( \dfrac{\partial L}{\partial \dot q} \delta q \right) dt = 0.

But the last integral gives \int_{t_1}^{t_2} \dfrac{d}{dt} \left( \dfrac{\partial L}{\partial \dot q} \delta q \right) dt = \dfrac{\partial L}{\partial \dot q} \delta q \Big|_{t_1}^{t_2} = \dfrac{\partial L}{\partial \dot q} \bigg[ \underbrace{\delta q (t_2) - \delta q(t_1)}_{=0} \bigg] = 0.

Therefore, we have \delta \mathcal{A} = \int_{t_1}^{t_2} \left[ \dfrac{\partial L}{\partial q} - \dfrac{d}{dt} \left( \dfrac{\partial L}{\partial \dot q} \right) \right] \delta q \: dt = 0. Since \delta q is arbitrarily selected, the integrand should be equal to zero in order to have the value of the integral null, or \dfrac{\partial L}{\partial q} = \dfrac{d}{dt} \left( \dfrac{\partial L}{\partial \dot q} \right). This concludes the derivation of Lagrange’s equation using Hamilton’s principle. However, one can derive Lagrange’s equation in a more direct way using calculus of variations or virtual work principles, see [11], [13], [16].

So far, we have considered systems that do not involve energy dissipation. In practice, however, we require extra terms in Lagrange’s equation to account for friction existing in real-world systems. Therefore, we expand the discussion to include non-conservative forces, e.g., friction and dampers, and find the corresponding Lagrange equation, including related topics such as cyclic coordinates, symmetry, multi-domain, and higher-order systems, [8], [13], [17].

2.6        Cyclic Coordinates

From Equation (2.6), it can be shown that if Lagrangian function does not have explicit dependency on one of the coordinates, say q_k, among all q_i, then the conjugate momentum p_k is conserved. The proof is as follows. Writing the Lagrange’s equation for coordinate q_k, we have \dfrac{d}{dt} \left( \dfrac{\partial L}{\partial \dot q_k} \right) - \dfrac{\partial L}{\partial q_k} = 0. Since by definition, L is not a function of q_k, then \dfrac{\partial L}{\partial q_k} = 0. Therefore, \dfrac{d}{dt} \left( \dfrac{\partial L}{\partial \dot q_k} \right) = 0, and written in terms of generalized momentum p_k, we get \dfrac{d}{dt} (p_k) = 0, or p_k is invariant with respect to time, hence conserved. It is common to call the coordinate q_k, cyclic or ignorable.

2.7        Conservative and Non-Conservative Forces

The generalized forces can be conservative or non-conservative. Conservative forces are those like gravity, buoyancy, mechanical spring, electrostatic, and magnetic. Non-conservative forces are those like friction, damping, and resistance.

By definition, a conservative force is curl free, or \vv{\nabla} \times \vv{F} = 0. Writing this expression in index notation, we have \mathcal{E}_{ijk}F_{k,j} = 0, where \mathcal{E}_{ijk} is the permutation symbol [18]. For example, force under gravity is (F_x,F_y,F_z)=(0,0,az). Calculating the curl gives \mathcal{E}_{123}F_{3,2}-\mathcal{E}_{132}F_{2,3} + \mathcal{E}_{231}F_{1,3} - \mathcal{E}_{213}F_{3,1} + \mathcal{E}_{312}F_{2,1} - \mathcal{E}_{321}F_{1,2}. Each term is identically zero; hence, the force under gravity field is conservative. Now, using the vector identity \vv \nabla \times \vv \nabla(V) = 0 or \mathcal{E}_{ijk}V_{,kj} = 0; i.e., the curl of a gradient of a scalar function is identically zero, and we can write a conservative force as the gradient of a scalar, such as potential function V as \vv{F_c} = -\vv{\nabla}V. By convention, the negative sign indicates that potential energy increases when work is done against a force field and vice versa.

We now, write Equation (2.7), after dropping the index i for simplicity, for T=T(\dot q) and V=V(q). Therefore, \dfrac{d}{dt} \left(\dfrac{\partial (T-V)}{\partial \dot q} \right) - \dfrac{\partial(T-V)}{\partial q} = 0, or \dfrac{d}{dt} \left( \dfrac{\partial T}{\partial \dot q} \right) + \dfrac{\partial V}{\partial q} = 0. But \dfrac{\partial V}{\partial q} = -F_c , and we get \dfrac{d}{dt} \left( \dfrac{\partial T}{\partial \dot q} \right) = F_c . This is the equation of motion (i.e. \dfrac{d}{dt} \bigg( \dfrac{\partial (\dfrac{1}{2} m \dot q^2)}{\partial \dot q} \bigg) = \dfrac{d \dot p}{dt}). We clearly see that the conservative force is already included in the Lagrange equation given by Equation (2.7). Now, for the case that we have a non-conservative force, or that the potential function is a function of velocity \dot q and q, (i.e. V=V(q,\dot q) or V=V(\dot q)), then we can write use Equation (2.7) to write \dfrac{d}{dt} \left( \dfrac{\partial T}{\partial \dot q} \right) - \dfrac{d}{dt} \left( \dfrac{\partial V}{\partial \dot q} \right) + \dfrac{\partial V}{\partial q} = 0. Re-arranging the term in this expression, we get \dfrac{d}{dt} \left( \dfrac{\partial T}{\partial \dot q} \right) = \dfrac{d}{dt} \left( \dfrac{\partial V}{\partial \dot q} \right) - \dfrac{\partial V}{\partial q}. We define the expression on the right-hand side as the non-conservative force, as F_{nc} = \dfrac{d}{dt} \left( \dfrac{\partial V}{\partial \dot q} \right) - \dfrac{\partial V}{\partial q}. Hence, \dfrac{d}{dt} \left( \dfrac{\partial T}{\partial \dot q} \right) = \dfrac{d \dot p}{dt} = F_{nc}. Again, we have shown that the non-conservative force is already included in the Lagrange equation given by Equation (2.7), provided a modified potential function is defined, as given by F_{nc}. See reference listed at [11] for more details.

2.8        Alternative form of Lagrange’s Equation

In section 2.7, we discussed the applicability of Lagrange’s equation given by Equation (2.7) for conservative and non-conservative forces. In practice, we could benefit from a more explicit form of the Lagrange equation whose terms can be easily identified for different types of forces, including energy dissipation such as damping and resistance. In this way, we can readily calculate the related terms in the Lagrange equation for modeling and simulation of a desired system.

Jean le Rond d’Alembert (1717–1783)

There are several possible ways to derive the Lagrange equation using, e.g., principles of virtual work and d’Alembert’s principle, directly from Newton’s second law of motion and first law of thermodynamics or energy conservation (e.g., conservation of sum of kinetic and potential energies) [8], [11], [13], [15], [17].

We use the conservation of energy approach to derive the alternative form of Equation (2.7) including its expansion [17].

We consider the kinetic energy of a system with generalized coordinates q_i for (i=1,2,\dots,n) (see section 2.4) represented by T=T(q_i,\dot q_i) and its potential energy by V=V(q_i,\dot q_i). Note that, as we discussed previously, for many mechanical systems kinetic energy is a function of \dot q_i and potential energy a function of q_i, only. Therefore, the resulted Lagrange equation can be simplified, accordingly. Now, using conservation of total energy of the system, we can write

(2.8)   \begin{equation*} d(T+V)=0 \end{equation*}

But dT = \dfrac{\partial T}{\partial q_i} dq_i + \dfrac{\partial T}{\partial \dot q_i} d \dot q_i and dV = \dfrac{\partial V}{\partial q_i} dq_i + \dfrac{\partial V}{\partial \dot q_i} d \dot q_i , using their functional relationships. After substituting into Equation (2.8), we get \dfrac{\partial T}{\partial q_i} dq_i + \dfrac{\partial T}{\partial \dot q_i} d \dot q_i + \dfrac{\partial V}{\partial q_i} + \dfrac{\partial V}{\partial \dot q_i} d \dot q_i = 0. Note that the Einstein summation convention applies, or dT = \dfrac{\partial T}{\partial q_1} dq_1 + \cdots + \dfrac{\partial T}{\partial q_n} dq_n + \dfrac{\partial T}{\partial \dot q_1} d \dot q_1 + \cdots + \dfrac{\partial T}{\partial \dot q_n} d \dot q_n . Now, using the relation for the kinetic energy of the system, or

(2.9)   \begin{equation*} T=\frac{1}{2}m_{ij} \dot q_i \dot q_j \end{equation*}

where m_{ij} is defined as the generalized mass matrix, a diagonally nonzero matrix, corresponding to the generalized coordinates. Therefore, its diagonal elements could be mass or moment of inertia when the generalized coordinates are displacement and angle, respectively. For example, for a n=2 system, we have

T= \dfrac{1}{2} m_{ij} \dot q_i \dot q_j = \dfrac{1}{2} (m_{1j} \dot q_1 \dot q_j + m_{2j} \dot q_2 \dot q_j) = \dfrac{1}{2}(m_{11} \dot q_1 \dot q_1 + m_{12} \dot q_1 \dot q_2 + m_{21} \dot q_2 \dot q_1 + m_{22} \dot q_2 \dot q_2)

With having \dot q_1 = \dot x and \dot q_1 = \dot \theta, and m_{11} = m, particle mass, m_{22} = I , inertia, and m_{12} = m_{21} = 0 we get T= \dfrac{1}{2} (m \dot x^2 + I \dot \theta^2 ). Now, differentiating T with respect to \dot q_i, we get \dfrac{\partial T}{\partial \dot q_i} = m_{ij} \dot q_j and substituting into Equation (2.9), we get T = \dfrac{1}{2} \dfrac{\partial T}{\partial \dot q_i} \dot q_i. Now, we calculate total change of T using the last expression, or 2dT = d \left( \dfrac{\partial T}{\partial \dot q_i} \dot q_i \right) = d \left( \dfrac{\partial T}{\partial \dot q_i} \right) \dot q_i + \dfrac{\partial T}{\partial \dot q_i}. But we had, dT = \dfrac{\partial T}{\partial q_i} dq_i + \dfrac{\partial T}{\partial \dot q_i} d \dot q_i}. Therefore, subtracting these last two relations, gives, after simplification, dT = d \left( \dfrac{\partial T}{\partial \dot q_i} \right) \dot q_i - \dfrac{\partial T}{\partial q_i} dq_i. But we can manipulate the first term on the right-hand side as d \left( \dfrac{\partial T}{\partial \dot q_i} \right) \dot q_i = d \left( \dfrac{\partial T}{\partial \dot q_i} \right) \dfrac{dq_i}{dt} = \dfrac{d}{dt} \left( \dfrac{\partial T}{\partial \dot q_i} \right) dq_i. Substituting into the last relation for dT, we get

(2.10)   \begin{equation*} dT = \left[ \frac{d}{dt} \left( \frac{\partial T}{\partial \dot q_i} \right) - \frac{\partial T}{\partial q_i} \right] dq_i \end{equation*}

Now, substituting Equation (2.10) into (2.8), we get \left[ \dfrac{d}{dt} \left( \dfrac{\partial T}{\partial \dot q_i} \right) - \dfrac{\partial T}{\partial q_i} \right] dq_i + dV = 0. Now, if V=V(q_i), i.e. holonomic systems, then we get dV = \dfrac{\partial V}{\partial q_i}dq_i} and, after substitution, we have \left[ \dfrac{d}{dt} \left( \dfrac{\partial T}{\partial \dot q_i} \right) - \dfrac{\partial T}{\partial q_i} + \dfrac{\partial V}{\partial q_i} \right]dq_i = 0. This expression is true for any arbitrarily selected dq_i; therefore, the terms in the bracket should be identically null, or

(2.11)   \begin{equation*} \frac{d}{dt} \left( \frac{\partial T}{\partial \dot q_i} \right) - \frac{\partial T}{\partial q_i} + \frac{\partial V}{\partial q_i} = 0 \end{equation*}

Equation (2.11), is an alternative form of Lagrange’s equation and holds when forces associated with the system are conservative, included in the \dfrac{\partial V}{\partial q_i} term. Note that using Lagrangian, L=T-V and Equation (2.11) we can recover Equation (2.7). The inclusion of non-conservative generalized forces, Q_i (usually the loading associated with each coordinate) should be added to the right-hand side of Equation (2.11). Also, energy dissipation due to viscous damping or resistance is usually given as D=D(\dot q_i^2) and contributes to Lagrange equation as \dfrac{\partial D}{\partial \dot q_i}. Finally, we get the alternative form of Lagrange equation, as

(2.12)   \begin{equation*} \frac{d}{dt} \left( \frac{\partial T}{\partial \dot q_i} \right) - \frac{\partial T}{\partial q_i} + \frac{\partial D}{\partial \dot q_i} + \frac{\partial V}{\partial q_i} = Q_i \quad , \: i = 1,2,\cdots,n \end{equation*}

Recall the n is the number of generalized coordinates. In matrix form, Equation (2.12) can be written as

    \begin{equation*} \dfrac{d}{dt} \begin{Bmatrix} \dfrac{\partial T}{\partial \dot q_1}\\ \vdots \\ \dfrac{\partial T}{\partial \dot q_n} \end{Bmatrix} - \begin{Bmatrix} \dfrac{\partial T}{\partial q_1}\\ \vdots \\ \dfrac{\partial T}{\partial q_n} \end{Bmatrix} + \begin{Bmatrix} \dfrac{\partial D}{\partial \dot q_1}\\ \vdots \\ \dfrac{\partial D}{\partial \dot q_n} \end{Bmatrix} + \begin{Bmatrix} \dfrac{\partial V}{\partial q_1}\\ \vdots \\ \dfrac{\partial V}{\partial q_n} \end{Bmatrix} = \begin{Bmatrix} Q_i\\ \vdots \\ Q_n \end{Bmatrix} \end{equation*}

2.9        Multi-Domain Systems

Lagrangian method can be applied to many kinds of engineering systems, including mechanical, electrical, thermal, hydraulic, and their possible combinations as multi-domain systems. As discussed in the previous sections, the established concept of generalized coordinates, momenta, and force are key tools to model such systems.

2.10       Systems with Higher Order Equations

System equations are mostly second-order differential equations, like Newton’s second law, and Kirchohff’s law for RCL circuits. Previous sections, e.g., Equation (2.7), presented Lagrange’s equation for such systems. One may require, mostly in continuous systems, to build the Lagrangian function for higher-order systems, e.g., fourth-order bi-harmonic equation for fluid flows or plate displacements. Fortunately, the Lagrangian method can be easily extended to cover the higher-order systems by considering a Lagrangian function, as given by Equation (2.13)

(2.13)   \begin{equation*} L=L(q_i, \dot q_i, \ddot q_i, \dots, t) \end{equation*}

Using the calculus of variations and Hamilton’s principle, we can derive the corresponding Lagrange’s equation [13], [9]. This s done by:

(2.14)   \begin{equation*} \frac{\partial L}{\partial q_i} + \sum_{m=1}^{k} (-1)^m \frac{d^m}{dt^m} \left( \frac{\partial L}{\partial q_{i,m}} \right) = 0 \quad i = 1,2,\cdots,n \end{equation*}

where m is the differentiation order; e.g., for m=3, we have

    \begin{equation*} \frac{\partial L}{\partial q_i} - \frac{d}{dt} \left( \frac{\partial L}{\partial \dot q_i} \right) + \frac{d^2}{dt^2} \left( \frac{\partial L}{\partial \ddot q_i} \right) - \frac{d^3}{dt^3} \left(\frac{\partial L}{\partial \dddot q_i} \right) = 0 \end{equation*}

Worked-out examples are useful to demonstrate applications of Lagrangian method. These examples, for mechanical and electrical systems, appear below. Each example includes numerical values assigned to the parameters and presents simulation results. Selected examples include accompanying screen-recorded video files demonstrating the solution steps for related system equations using 20-sim. After learning from the related video file, the reader can modify the parameters and run the simulation for specific design cases.

2.11       Example: A Multi-Mass-Spring System

We want to find the equations governing its motion dynamics for the system sketched in Figure 2-5. For this example, we neglect the effect of gravity.

Figure 2-5 A mass-spring system with three degrees of freedom

This system has three degrees of freedom (x_1,x_2,x_3 ) associated with three masses (m_1,m_2,m_3 ). For three masses, N=3, and each can move vertically; hence, the number of constraints is N_c=2 for each mass. This gives n = 3 \times 3 - 3 \times 2 = 3. The Lagrangian method is used to find the equations of motion, or three coupled second-order differential equations. We start by writing the kinetic and potential energy expressions of the system and forming the corresponding Lagrangian. The kinetic energy of the system is T= \dfrac{1}{2} (m \dot x_1^2 + m \dot x_2^2 + m \dot x_3^2). For the potential energy, we should use the difference in displacements associated with each spring because the neutral position of the unstressed springs do not contribute to the potential energy. For example, for the spring k_5, connecting masses m_1 and m_3, we should use X =| x_1 - x_3 | as the variable, or \int k_5XdX = \dfrac{1}{2} k_5(x_1-x_3)^2. Therefore, the potential energy of the system consisting of the sum of all springs is V= \dfrac{1}{2}(k_1x_1^2 + k_2x_2^2 + k_3x_3^2) + \dfrac{1}{2}k_4(x_2 - x_3)^2 + \dfrac{1}{2}k_5(x_1 - x_3)^2 + \dfrac{1}{2}k_6(x_1 - x_2)^2. Note that for this system the kinetic energy is a function of only \dot x_i and potential energy a function of x_i. Applying Euler-Lagrange equation to each mass, or degree of freedom, we get a system of ODEs, written in matrix form,

    \begin{equation*} \begin{bmatrix} m_1 & 0 & 0\\ 0 & m_2 & 0\\ 0 & 0 & m_3 \end{bmatrix} \begin{Bmatrix} \ddot x_1\\ \ddot x_2\\ \ddot x_3 \end{Bmatrix} + \begin{bmatrix} k_1 + k_5 + k_6 & -k_6 & -k_5\\ -k_6 & k_2 + k_4 + k_6 & -k_4\\ -k_5 & -k_4 & k_3 + k_4 + k_5 \end{bmatrix} \begin{Bmatrix} x_1\\ x_2\\ x_3 \end{Bmatrix} = \begin{Bmatrix} 0\\ 0\\ 0 \end{Bmatrix} \end{equation*}

For example, the Euler-Lagrange equation associated with mass m_1 reads \dfrac{d}{dt} \left( \dfrac{\partial L}{\partial \dot x_1} \right) - \dfrac{\partial L}{\partial x_1} = 0. But we have \dfrac{\partial L}{\partial \dot x_1} = \dfrac{\partial (T-V)}{\partial \dot x_1} = \dfrac{\partial T}{\partial \dot x_1} = m \ddot x_1^2 and - \dfrac{\partial L}{\partial x_1} = - \dfrac{\partial (T-V)}{\partial x_1} = \dfrac{\partial V}{\partial x_1} = k_1x_1 + k_5x_1 + k_6x_1 - k_6x_2 - k_5x_1. Having information about initial and boundary conditions for displacements and/or velocities, we can obtain the solution of the system’s equations using 20-sim. An initial velocity of 0.2 m/s is applied to mass m_2, for example. The script code is as follows:

parameters

real m1 = 15.0 {kg};
real m2 = 30.0 {kg};
real m3 = 15.0 {kg};
real k1 = 50.0 {N/m};
real k2 = 100.0 {N/m};
real k3 = 50.0 {N/m};
real k4 = 20.0 {N/m};
real k5 = 70.0 {N/m};
real k6 = 10.0 {N/m};

variables

real x1 {m};
real x2 {m};
real x3 {m};
real x1_dot {m/s}; // velocity
real x2_dot {m/s}; // velocity
real x3_dot {m/s}; // velocity
real x1_dot_dot {m/s2}; //acceleration
real x2_dot_dot {m/s2}; //acceleration
real x3_dot_dot {m/s2}; //acceleration
real Fk1 {N}; // force spring k1
real Fk2 {N}; // force spring k2
real Fk3 {N}; // force spring k3

equations

x1_dot_dot = -(1/m1)*((k1+k5+k6)*x1-k6*x2-k5*x3);
x2_dot_dot = -(1/m2)*((k2+k4+k6)*x2-k6*x1-k4*x3);
x3_dot_dot = -(1/m3)*((k3+k4+k5)*x3-k4*x2-k5*x1);
x1_dot = int (x1_dot_dot , 0);
x2_dot = int (x2_dot_dot , 0.2); //initial velocity 0.2m/s
x3_dot = int (x3_dot_dot , 0);
x1 = int (x1_dot , 0.2); //initial displacement 0.2m
x2 = int (x2_dot , 0);
x3 = int (x3_dot , -0.1); //initial displacement -0.1m
Fk1 = k1*x1;
Fk2 = k2*x2;
Fk3 = k3*x3;

The results for displacements of the masses and velocities are shown below, see Figure 2-6.

Figure 2-6 Sample results as output from 20-sim

Here is a video showing how to build and run the model for this example in 20-sim:

 

2.12       Example: A System with Energy Dissipation and Applied External Force

We consider a system with two degrees of freedom, as shown in Figure 2-7. The damping coefficients b_1 and b_2 and spring stiffness k_1 and k_2 are used to calculate the potential and damping functions V and D, respectively.

Figure 2-7 A mass-spring-damper system with two degrees of freedom

The non-conservative Rayleigh energy dissipation function is, D= \dfrac{1}{2}b_1 \dot x_1^2 + \dfrac{1}{2} b_2( \dot x_2 - \dot x_1)^2. The derivative of this function with respect to \dot x_i, \dfrac{\partial D}{\partial \dot x_i} gives the damping forces associated with mass m_i. The kinetic energy is T= \dfrac{1}{2}(m_1 \dot x_1^2 + m_2 \dot x_2^2), and potential energy reads V = \dfrac{1}{2}k_1x_1^2 + \dfrac{1}{2}k_2(x_2 - x_1)^2

Lagrange’s equation for motion of mass m_1 reads \dfrac{d}{dt} \left( \dfrac{\partial L}{\partial \dot x_1} \right) - \dfrac{\partial L}{\partial x_1} + \dfrac{\partial D}{\partial \dot x_1} = 0 and for mass m_2 is \dfrac{d}{dt} \left( \dfrac{\partial L}{\partial \dot x_2} \right) - \dfrac{\partial L}{\partial x_2} + \dfrac{\partial D}{\partial \dot x_2} = F(t). Performing the derivatives, we get \dfrac{\partial T}{\partial \dot x_i} = m_1 \dot x_1, \dfrac{\partial T}{\partial \dot x_2} = m_2 \dot x_2, \dfrac{\partial V}{\partial x_1} = (k_1 + k_2)x_1 - k_2x_2, \dfrac{\partial V}{\partial x_2} = k_2x_2 - k_2x_1, \dfrac{\partial D}{\partial \dot x_1} = (b_1 + b_2) \dot x_1 - b_2 \dot x_2, \dfrac{\partial D}{\partial \dot x_2} = b_2 \dot x_2 - b_2 \dot x_1.

Using Lagrange’s equation, with L=T-V, we get the equations of motion of the system in matrix form as

    \begin{equation*} \begin{bmatrix} m_1 & 0\\ 0 & m_2 \end{bmatrix} \begin{Bmatrix} \ddot x_1\\ \ddot x_2 \end{Bmatrix} + \begin{bmatrix} b_1 + b_2 & -b_2\\ -b_2 & b_2 \end{bmatrix} \begin{Bmatrix} \dot x_1\\ \dot x_2 \end{Bmatrix} + \begin{bmatrix} k_1 + k_2 & -k_2\\ -k_2 & k_2 \end{bmatrix} \begin{Bmatrix} x_1\\  x_2 \end{Bmatrix} = \begin{Bmatrix} 0\\ F(t) \end{Bmatrix} \end{equation*}

We use 20-sim to solve the systems equations. A step function is used for applied force. The script code is as follows:

parameters

real m1 = 2.0 {kg};
real m2 = 1.0 {kg};
real k1 = 20.0 {N/m};
real k2 = 30.0 {N/m};
real b1 = 0.1 {N.s/m};
real b2 = 0.05 {N.s/m};
real start_time = 10 {s};
real amplitude = 0.5;

variables

real x1 {m};
real x2 {m};
real x1_dot {m/s};
real x2_dot {m/s};
real x1_dot_dot {m/s2};
real x2_dot_dot {m/s2};
real F_applied {N}; // applied force

equations

x1_dot_dot = -(b1+b2)/m1*x1_dot+b2/m1*x2_dot-(k1+k2)/m1*x1+k2/m1*x2;
x2_dot_dot = -(1/m2)*(-b2*x1_dot+b2*x2_dot-k1*x1+k2*x2+F_applied);
x1_dot = int (x1_dot_dot , 0);
x2_dot = int (x2_dot_dot , 0);
x1 = int (x1_dot , 0);
x2 = int (x2_dot , 0);
F_applied = amplitude*step (start_time);

The results for displacements of the masses and applied force are shown below, see Figure 2-8.

Figure 2-8 Sample results as output from 20-sim

Here is a video showing how to build and run the model for this example in 20-sim:

 

2.13       Example: A Two-Loop Electrical Circuit

For this example, we consider an electrical circuit with two loops/branches. For the system, we have; electric charges q_1 and q_2; resistors R_1 and R_2; inductors L_1, L_2, and L_3; and capacitors C_1 and C_2 as Figure 2-9 shows. The voltage source is u(t), connected to loop 1.

Figure 2-9 A two-loop electrical circuit with source

For comparison with a typical mechanical system, the equivalent of mass is an inductor; for spring, a capacitor; and for damper, a resistor. Therefore, using the Lagrangian method, we can write the kinetic energy of the system as T = \dfrac{1}{2}L_1 \dot q_1^2 + \dfrac{1}{2}L_2 \dot q_2^2 + \dfrac{1}{2}L_3(\dot q_1 - \dot q_2)^2. Note that electric charge is analogous to mechanical displacement and electric current to velocity, or q_i \equiv x_i and \dot q_i \equiv \dot x_i. Therefore, e.g., the term \dfrac{1}{2} L_1 \dot q_1^2 represents the stored kinetic energy in the corresponding inductor. Similarly, the potential energy is V= \dfrac{1}{2C_1}q_1^2 + \dfrac{1}{2C_2} q_2^2. Note that the capacitance is the inverse of stiffness, or C_i = \dfrac{1}{k_i}. The energy dissipation function for the system is D= \dfrac{1}{2}R_1 \dot q_1^2 + \dfrac{1}{2}R_2 \dot q_2^2. Using Langrange’s equation, \dfrac{d}{dt} \left( \dfrac{\partial (T-V)}{\partial \dot q_1} \right) - \dfrac{\partial (T-V)}{\partial q_i} + \dfrac{\partial D}{\partial q_i} = F_i, gives the electric circuit system equations as

    \begin{equation*} \left\{ \begin{array} (L_1 + L_3) \ddot q_1 - L_3 \ddot q_2 + R_1 \dot q_1 + \dfrac{q_1}{C_1} = u(t)\\ (L_2 + L_3) \ddot q_2 - L_3 \ddot q_1 + R_2 \dot q_2 + \dfrac{q_2}{C_2} = 0 \end{array} \end{equation*}

One can use rate of charge or the electric current, I as the variable by replacing I= \dot q in the system’s equations. This gives (L_1 + L_3) \dot I_1 - L_3 \dot I_2 + R_2I_1 + V_{c1} = u(t) and (L_2 + L_3) \dot I_2 - L_3 \dot I_1 + R_2I_2 + V_{c2} = 0 where V_{c1} and V_{c2} are the voltage across the capacitors, respectively.

We use 20-sim to solve the system equations. The script code is as follow

parameters

real L1 = 0.15 {H};
real L2 = 0.2 {H};
real L3 = 0.25 {H};
real C1 = 0.05 {F};
real C2 = 0.02 {F};
real R1 = 1 {ohm};
real R2 = 2 {ohm};
real omega = 3 {rad/s};
real amplitude = 1;

variables

real q1 {C};
real q2 {C};
real q1_dot {A};
real q2_dot {A};
real q1_dot_dot ;
real q2_dot_dot ;
real Voltage {V}; // applied voltage

equations // equations are manipulated

q2_dot_dot*(L1*L2+L2*L3+L1*L3)=-L3*R1*q1_dot-(L1+L3)*R2*q2_dot-L3/C1*q1-(L1+L3)/C2*q2+Voltage*L3;
q1_dot_dot*(L3) = (L2+L3)*q2_dot_dot+R2*q2_dot+(1/C2)*q2;
q1_dot = int (q1_dot_dot , 0);
q2_dot = int (q2_dot_dot , 0);
q1 = int (q1_dot , 0);
q2 = int (q2_dot , 0);
Voltage = amplitude*sin (omega*time);

Typical plots for current in each loop is shown in Figure 2-10 for a sinusoidal voltage.

Figure 2-10 Sample results as output from 20-sim

2.14       Example: A Compound Atwood’s Machine

Atwood’s machine is a collection of pulleys and masses. This example examines and models the dynamical behavior of this machine as shown in Figure 2-11.

Figure 2-11 A compound Atwood’s machine

This system has two degrees of freedom ( \underbrace{3 \times 3}_{3N} - \underbrace{3 \times 2 - 1}_{N_c} = 2) describing the motion of mass m_1 and pulley b. Therefore, two ODEs describe the system dynamical behaviour. The massless un-stretchable string length hanging over pulley a is l_a, and that of pulley b is l_b. We measure the potential energy with reference to the top of pulley a with vertical displacement designated with x and similarly from top of pulley b with y, as shown in Figure 2-11. The kinetic energy reads T = \dfrac{1}{2} (m_1 \dot x_1^2 + m_2 \dot x_2^2 + m_3 \dot x_3^2 ), where , using the geometrical constraints and string lengths, x_1 = x, x_2 = (l_a - x_1) + y, x_3 = (l_a - x_1) + (l_b - y). Therefore, \dot x _1 = \dot x, \dot x_2 = - \dot x + \dot y, \dot x_3 = -\dot x - \dot y. Substituting in kinetic energy relation, gives T= \dfrac{1}{2} [m_1 \dot x^2 + m_2 (\dot y - \dot x)^2 + m_3 (\dot y + \dot x)^2 ]. The potential energy reads V = -g(m_1x_1 + m_2x_2 + m_3x_3). After substituting for x_1,x_2, and x_3 and algebraic simplifications we get V = xg(m_2 + m_3 - m_1) + yg(m_3 - m_2 ) + C , where constant C is given by C = -g(m_2l_a + m_3l_a + m_3l_b). The Langrange equations in terms of x and y are \dfrac{d}{dt} \left( \dfrac{\partial L}{\partial \dot x} \right) - \dfrac{\partial L}{\partial x} = 0 and

\dfrac{d}{dt} \left( \dfrac{\partial L}{\partial \dot y} \right) - \dfrac{\partial L}{\partial y} = 0, having

    \begin{equation*} L = T - V = \dfrac{1}{2} [m_1 \dot x^2 + m_2 ( \dot y - \dot x)^2 + m_3 (\dot y + \dot x)^2 ] -xg(m_2 + m_3 - m_1) - yg(m_3 - m_2) \end{equation*}

We dropped C, since its differentiation is zero. Hence,

    \begin{equation*} \dfrac{d}{dt} \left( \dfrac{\partial L}{\partial \dot x} \right) = \ddot x (m_1 + m_2 + m_3 ) + \ddot y (m_3 - m_2), \end{equation*}

    \begin{equation*} \dfrac{d}{dt}\left( \dfrac{\partial L}{\partial \dot y} \right) = \ddot x (m_3 - m_2 ) + \ddot y (m_2 + m_3), \end{equation*}

    \begin{equation*} \dfrac{\partial L}{\partial x} = g(m_1 - m_2 - m_3) \end{equation*}

    \begin{equation*} \text{and} \dfrac{\partial L}{\partial y} = g(m_2 - m_3). \end{equation*}

Substituting into the corresponding Lagrange equations, we get the system’s equations of motion as

    \begin{equation*} \left\{ \begin{array} x \ddot x(m_1 + m_2 + m_3) + \ddot y(m_3 - m_32) = g(m_1 - m_2 - m_3)\\ \ddot x(m_3 - m_2) + \ddot y(m_2 + m_3) = g(m_2 - m_3) \end{array} \end{equation*}

To simplify the equations, eliminate \ddot y by multiplying the first equation by (m_2 + m_3) and the second one by (m_2 - m_3 ). After some manipulations, we get

    \begin{equation*} \left\{ \begin{array} x \ddot x(m_1m_2 + m_1m_3 + 4m_2m_3) = g(m_1 - m_2 - m_3)(m_2 - m_3)\\ \ddot y(m_2 + m_3) = (\ddot x + g)(m_2 - m_3) \end{array} \end{equation*}

We use 20-sim to solve these system equations. The script code is as follows:

parameters

real m1 = 1.0 {kg};
real m2 = 2.0 {kg};
real m3 = 4.0 {kg};
real g = 9.08 {m/s2};

variables

real x {m};
real y {m};
real x_dot {m/s};
real y_dot {m/s};
real x_dot_dot {m/s2};
real y_dot_dot {m/s2};

equations

/* x_dot_dot = (1/(m1+m2+m3))*(-y_dot_dot*(m3-m2)+g*(m1-m2-m3)); */
x_dot_dot = g*(m1-m2-m3)*(m2-m3)/(m1*m2+m1*m3+4*m2*m3);
y_dot_dot = (1/(m3+m2))*((x_dot_dot+g)*(m2-m3));
x_dot = int (x_dot_dot , 0);
y_dot = int (y_dot_dot , 0);
x = int (x_dot , 0);
y = int (y_dot , 0.1);

2.15       Example: Atwood’s Machine with Massive String and Pulley

In the analysis of Atwood’s machine, the pulley and string are usually considered massless. In this example, we include these parts, assuming the string having mass m_s, total length l, and linear mass density \rho = m_s / l_s and the pulley with mass M, radius R, and moment of inertia I [19].

Figure 2-12 Atwood’s machine

2Datum for potential energy is a horizontal plane at the level of the pulley’s centre. From the datum, the length of hanging string on the two sides of pulley is - \pi R. The potential energy is due to the masses and the string mass, or V = m_1gx - m_2g(l-x) - \dfrac{1}{2} \rho gx^2 - \dfrac{1}{2} \rho g(l - x)^2. Note that x is measured downward from the datum toward mass m_1. The kinetic energy is due to the masses, string, and the pulley’s angular kinetic energy, \dfrac{1}{2} I \omega^2 with angular velocity \omega = \dot x/R. Therefore, T = \dfrac{1}{2} (m_1 + m_2 + m_s ) \dot x^2 + \dfrac{1}{2} I ( \dot x/R)^2. The Lagrangian is written as

    \begin{equation*} L = T - V = \dfrac{1}{2} \left( m_1  + m_2  + m_s + \dfrac{I}{R^2} \right) \dot x^2 + \\ (m_1 - m_2) gx + \rho gx^2 - \rho glx + \left( \dfrac{1}{2} \rho gl^2 + m_2 gl \right) \end{equation*}

The Lagrange is equation reads \dfrac{d}{dt} \left( \dfrac{\partial L}{\partial \dot x} \right) - \dfrac{\partial L}{\partial x} = 0, or \ddot x = \dfrac{(m_1 - m_2 + \rho x  - \rho l)g}{(m_1 + m_2 + \rho l + \rho \pi R + I/R^2)}, after substituting m_s = \rho (l + \pi R). The result reduces to the familiar result of \ddot x = \dfrac{(m_1 - m_2 ) g}{(m_1  + m_2)} for massless string and pulley ( \rho = 0, I = 0).

We use 20-sim to solve these system equations. The script code is as follows:

parameters

real Ms = 2.5 {kg}; // string mass
real L = 2.0 {m}; //string length
real M = 3.0 {kg}; //mass of the pulley
real R = 30.0 {cm}; //radius of the pulley
real g = 9.08 {m/s2}; // grav. acceleration
real m1 = 4.0 {kg};
real m2 = 1.5 {kg};

variables

real x {m}; //vertical displacement
real I {kg.m2}; // pulley moment of inertia
real x_dot {m/s}; // vertical velocity
real x_dot_dot {m/s2}; // vertical acceleration

equations

I = 0.5*M*R^2;
x_dot_dot = g*(m1-m2+(Ms/L)*(x-L))/((m1+m2+(Ms/L)*(L+ pi*R)+I/R^2));
x_dot = int (x_dot_dot , 0.0);
x = int (x_dot , 0);

2.16       Example: A Complex Vibrating Mechanical System

For this example, we consider a mechanical system with three degrees of freedom, x_1, x_2, x_3, associated with three masses,m_1, m_2, m_3. The arrangement of springs and dampers is shown, with their coefficients, in Figure 2-13, with corresponding stiffness (k_1, k_2) and damping (b_1, b_2, b_3) coefficients. An applied force, f(t) acting on mass m_2 and all wall contact surfaces are considered to have negligible friction.

Figure 2-13 A complex vibrating mechanical system

The kinetic energy of the systems reads T = \dfrac{1}{2} (m_1 \dot x_1^2 + m_2 \dot x_2^2 + m_3 \dot x_3^2) and the potential energy is V = \dfrac{1}{2} k_1 x_1^2 + \dfrac{1}{2} k_2 (x_2 - x_3 )^2. Similarly, the damping function reads D = \dfrac{1}{2} b_1 \dot x_2^2 + \dfrac{1}{2} b_2 (\dot x_2 - \dot x_1 )^2 + \dfrac{1}{2} b_3 \dot x_3^2. The Lagrange’s equations are \dfrac{d}{dt} \left( \dfrac{\partial (T-V)}{\partial \dot x_i} \right) - \dfrac{\partial (T-V)}{\partial x_i} + \dfrac{\partial D}{\partial \dot ẋ_i} = F_i (t), with F_i (t) = \Bigg\{ \begin{array} x0\\ f(t)\\ 0 \end{array} \Bigg\}

because the applied force is exerted on mass m_2. Performing the differentiations, we can write the equations of the system, as

    \begin{equation*} \dfrac{d}{dt} (m_1 \dot x_1 ) + k_1 x_1 + b_2 (\dot x_1 - \dot x_2 ) = 0, \dfrac{d}{dt} (m_2 \dot x_2 ) + \\ k_2 (x_2 - x_3 ) + b_1 \dot x_2 + b_2 (\dot x_2 - \dot x_1 ) = f(t)$ \end{equation*}

and \dfrac{d}{dt} (m_3 \dot x_3 ) + k_2 (x_3  - x_2 ) + b_3 \dot x_3 = 0. In matrix form, the system’s equations are

    \begin{equation*} \begin{bmatrix} m_1 & 0 & 0\\ 0 & m_2 & 0\\ 0 & 0 & m_3 \end{bmatrix} \begin{Bmatrix} \ddot x_1\\ \ddot x_2\\ \ddot x_3 \end{Bmatrix} + \begin{bmatrix} b_2 & -b_2 & 0\\ -b_2 & b_1 + b_2 & 0\\ 0 & 0 & b_3 \end{bmatrix} \begin{Bmatrix} \dot x_1\\ \dot x_2\\ \dot x_3 \end{Bmatrix} + \begin{bmatrix} k_1 + & 0 & 0\\ 0 & k_2 & -k_2\\ 0 & -k_2 & k_2 \end{bmatrix} \begin{Bmatrix} x_1\\ x_2\\ x_3 \end{Bmatrix} = \begin{Bmatrix} 0\\ f(t)\\ 0 \end{Bmatrix} \end{equation*}

We use 20-sim to solve these system equations. The applied force is composed of three impulses applied at 5, 10, and 20 second. The script code is as follows:

parameters

real m1 = 1.0 {kg};
real m2 = 3.0 {kg};
real m3 = 2.0 {kg};
real k1 = 50.0 {N/m};
real k2 = 30.0 {N/m};
real b1 = 0.1 {N.s/m};
real b2 = 0.2 {N.s/m};
real b3 = 0.3 {N.s/m};

variables

real x1 {m};
real x2 {m};
real x3 {m};
real x1_dot {m/s};
real x2_dot {m/s};
real x3_dot {m/s};
real x1_dot_dot {m/s2};
real x2_dot_dot {m/s2};
real x3_dot_dot {m/s2};
real F_applied1 {N};
real F_applied2 {N};
real F_applied3 {N};

equations

x1_dot_dot = -b2/m1*x1_dot+b2/m1*x2_dot-k1/m1*x1;
x2_dot_dot = -(1/m2)*(-b2*x1_dot+(b1+b2)*x2_dot+k2*x2-k2*x3+F_applied1+F_applied2+F_applied3);
x3_dot_dot = -(1/m3)*(b3*x3_dot-k2*x2+k2*x3);
x1_dot = int (x1_dot_dot , 0);
x2_dot = int (x2_dot_dot , 0);
x3_dot = int (x3_dot_dot , 0);
x1 = int (x1_dot , 0);
x2 = int (x2_dot , 0);
x3 = int (x3_dot , 0);
F_applied1 = 3*impulse (5,0.1);
F_applied2 = 5*impulse (20,0.2);
F_applied3 = -10*impulse (10,0.2);

Sample results are shown in Figure 2-14.

Figure 2-14 Sample results as output from 20-sim

Here is a video showing how to build and run the model for this example in 20-sim:

 

 

2.17       Example: A Pendulum with Moving Pivot

A simple pendulum with mass m hanging from a free-moving pivot with mass M. The system has two degrees of freedom: oscillation of pivot, x = x(t) and pendulum motion about vertical designated by angle \theta = \theta (t). The pendulum string with length l is massless and unstretchable. We consider the datum at the pivot level and gravitational acceleration g pointing downwards, as in Figure 2-15.

Figure 2-15 Pendulum with oscillating pivot

Mass m coordinates read (x + l \text{sin} \theta, - l \text{cos} \theta); hence, the velocity components are (\dot x + l \dot \theta \text{cos} \theta, l \dot \theta \text{sin} \theta). We can write kinetic energy of the system as

    \begin{equation*} T = \frac{m}{2} \left[ (\dot x + l \dot \theta \text{cos} \theta)^2 + (l \dot \theta \text{sin} \theta )^2 \right] + \frac{M}{2} \dot x^2 = \frac{m}{2} (\dot x^2 + l^2 \dot \theta ^2 + 2l \dot x \dot \theta \text{cos} \theta) + \frac{M}{2} \dot x^2 \end{equation*}

Similarly, the potential energy of the systems reads V = mgy = -mgl \text{cos} \theta. Note that the pivot motion is horizontal with coordinates (x, 0). The Lagrange equation for rotational motion with respect to coordinate \theta reads \dfrac{d}{dt} \left( \dfrac{\partial (T-V)}{\partial \dot \theta} \right) - \dfrac{\partial (T-V)}{\partial \theta} = 0, or \dfrac{d}{dt} (ml^2 \dot \theta + ml \dot x \text{cos} \theta ) + ml \dot x \dot \theta \text{sin} \theta - mgl \text{sin} \theta = 0. After simplification, we get \ddot \theta - \dfrac{g}{l} \text{sin} \theta + \ddot x \text{cos} \theta = 0. Note that for fixed pivot (or x = \dot x = \ddot x = 0) we get the familiar result for a simple pendulum. The Lagrange equation for translational motion with respect to coordinate x reads \dfrac{d}{dt} \lefct( \dfrac{\partial (T-V)}{\partial \dot x} \right) - \dfrac{\partial (T-V)}{\partial x} = 0, or \dfrac{d}{dt} (m \dot x + ml \dot \theta \text{cos} \theta + M \dot x) = 0. After performing differentiation, we get (M + m) \ddot x + (ml \text{cos} \theta ) \ddot \theta - (ml \text{sin} \theta ) \dot \theta ^2 = 0. Collectively, the system’s equations of motion are

    \begin{equation*} \left\{ \begin{array} ( (ml \text{cos}^2 \theta - M - m) \ddot \theta - (ml \text{sin} \theta} \text{cos} \theta ) \dot \theta ^2 + \frac{g}{l} (M + m) \text{sin} \theta = 0 \\ \ddot x = \frac{1}{\text{cos} \theta} \left( \frac{g}{l} \text{sin} \theta - \ddot \theta \right) \end{array} \end{equation*}

We use 20-sim to solve these system equations. An initial velocity of 0.5 rad/s is applied to the pendulum. The script code is as follows:

parameters

real m = 0.5 {kg}; // pendulum/bob mass
real M = 1.0 {kg}; // pivot mass
real g = 9.08 {m/s2}; //gravity
real L = 30 {cm}; //pendulum length

variables

real x {m};
real x_dot {m/s};
real x_dot_dot {m/s2};
real theta {rad};
real theta_dot {rad/s};
real theta_dot_dot {rad/s2};

equations

x_dot_dot = (1/cos (theta))*((g/L)*sin (theta)-theta_dot_dot);
theta_dot_dot = (1/(m*L*cos (theta)^2-M-m))*(m*L*sin (theta)*cos (theta)*theta_dot^2-g/L*(m+M)*sin (theta));
x_dot = int (x_dot_dot , 0);
x = int (x_dot , 0);
theta_dot = int (theta_dot_dot , 0.5);
theta = int (theta_dot , 0);

Sample results are shown in Figure 2-16.

Figure 2-16 Sample results as output from 20-sim

Here is a video showing how to build and run the model for this example in 20-sim:

 

2.18       Example: A Pendulum Attached to a Moving Mass-Spring-Damper System

In this example we consider a system consisting of a pendulum with its pivot attached to the centre of a freely moving mass M. The mass is connected to a spring with stiffness k and a damper with damping coefficient b. The pendulum bob has a mass of m and is attached to a torsional damper with damping coefficient b_t and a torsional spring with stiffness k_t. The pendulum string is massless and has a length of l. We consider the datum at the pivot level and gravitational acceleration g pointing downwards, as in Figure 2-17. The system has two degrees of freedom; oscillation of pivot, x=x(t) and pendulum motion about vertical direction designated by angle \theta = \theta (t).

Figure 2-17 A pendulum attached to a mass-spring-damper system

The coordinates of mass m read (x + l \text{sin} \theta , - l \text{cos} \theta ), and its velocity components are (\dot x + l \dot \theta \text{cos} \theta , l \dot \theta \text{sin} \theta ). We can write kinetic energy of the system as

    \begin{equation*} T = \frac{m}{2} \left[ (\dot x + l \dot \theta \text{cos} \theta )^2 + (l \dot \theta \text{sin} \theta )^2 \right] + \frac{M}{2} \dot x^2 = \frac{m}{2} (\dot x^2 + l^2 \dot \theta ^2 + 2l \dot x \dot \theta \text{cos} \theta ) + \frac{M}{2} \dot x^2 \end{equation*}

Similarly, the potential energy of the systems reads V = \dfrac{1}{2} kx^2 + \dfrac{1}{2} k_t \theta ^2 - mgl \text{cos} \theta. The damping function of the system is D = \dfrac{1}{2} b \dot x^2 + \frac{1}{2} b_t \dot \theta ^2.

The Lagrange equation for rotational motion with respect to coordinate \theta reads \dfrac{d}{dt} \left( \dfrac{\partial (T-V)}{\partial \dot \theta} \right) - \dfrac{\partial (T-V)}{\partial \theta} + \dfrac{\partial D}{\partial \dot \theta} = 0, or \dfrac{d}{dt} (ml^2 \dot \theta + ml \dot x \text{cos} \theta ) + ml \dot x \dot \theta \text{sin⁡} \theta + k_t \theta - mgl \text{sin} ⁡\theta = 0. After simplification, we get \ddot \theta - \dfrac{g}{l} \text{sin} \theta + \dfrac{\ddot x \text{cos} \theta}{l} + \dfrac{1}{l^2} (b_t \dot \theta + k_t \theta ) = 0. The Lagrange equation for translational motion with respect to coordinate x reads \dfrac{d}{dt} \left( \dfrac{\partial (T-V)}{\partial \dot x} \right) - \dfrac{\partial (T-V)}{\partial x} + \dfrac{\partial R}{\partial \dot x} = F(t), or \dfrac{d}{dt} (m \dot x + ml \dot \theta \text{cos} \theta + M \dot x) = 0. After performing differentiation, we get (M + m) \ddot x + (ml \text{cos} \theta ) \ddot \theta - (ml \text{sin} \theta ) \dot \theta ^2 + kx + b \dot x = F(t). Collectively, the system’s equations of motion are

    \begin{equation*} \left\{ \begin{array} ( (M + m) \ddot x + ( ml \text{cos} \theta) \ddot \theta - (ml \text{sin} \theta ) \dot \theta ^2 + kx + b \dot x = F(t) \\ \ddot \theta + \dfrac{g}{l} \text{sin} \theta + \dfrac{\ddot x \text{cos} \theta}{l} + \dfrac{1}{ml^2} (b_t \dot \theta + k_t \theta ) = 0 \end{array} \end{equation*}

We use 20-sim to solve these system equations. The script code is as follows:

parameters

real m = 0.5 {kg}; // pendulum/bob mass
real M = 1.0 {kg}; // pivot mass
real g = 9.08 {m/s2}; //gravity
real L = 30 {cm}; //pendulum length
real k = 2 {N/m}; // spring stiffness
real kt = 0.5 {N.m/rad}; // torsional stiffness
real bt = 0.5 {N.m.s/rad}; // torsional damping
real b = 0.2 {N.s/m}; // damping
real amplitude = 1; // amplitude of applied force
real omega = 0.5 {rad/s}; // frequency of applied force

variables

real x {m};
real x_dot {m/s};
real x_dot_dot {m/s2};
real theta {rad};
real theta_dot {rad/s};
real theta_dot_dot {rad/s2};
real F_applied {N};
real F_spring {N}; // linear spring force
real T_spring {N.m}; // torsional spring torque
real y; // aux. variable, to help the solver

equations

x = int (x_dot , 0);
x_dot = int (x_dot_dot , 0);
theta = int (theta_dot , 0);
theta_dot = int (theta_dot_dot , 0);
y = -m*L*cos (theta)*(theta_dot_dot);
x_dot_dot = (1/(m+M))*(m*L*sin (theta)*theta_dot^2 + y -k*x-b*x_dot+F_applied);
theta_dot_dot = -g/L*sin (theta) -1/L*cos (theta)*x_dot_dot -1/(m*L^2)*(bt*theta_dot+kt*theta);
F_applied = amplitude*sin (omega*time);
F_spring = k*x;
T_spring = kt*theta;

Cart displacement, pendulum angle, and force and torque of the springs are shown in Figure 2-18.

Figure 2-18 Sample results as output from 20-sim

Here is a video showing how to build and run the model for this example in 20-sim:

 

2.19       Example: A Mass Particle Sliding on a Rotating Circular Ring

Figure 2-19 shows a particle with mass m sliding on a circular ring with radius R. The ring itself is rotating about the z-axis with a constant angular velocity \omega. We want to find the equation of motion for the mass particle.

Figure 2-19 A particle moving on a circular ring

The generalized coordinate is \theta = \theta (t), the polar angle. We can write the coordinates of the mass particle as x = R \text{sin}⁡ \theta \text{cos} \omega t, y = R \text{sin} \theta \text{sin} \omega t, and z = R \text{cos} \theta⁡. Therefore, \dot x = R \dot \theta \text{cos} \theta \text{cos} \omega t - R \omega \text{sin} \theta \text{sin} \omega t, \dot y = R \dot \theta \text{cos} \theta \text{sin} \omega t + R \omega \text{sin} \theta \text{cos} \omega t, and \dot z = -R \dot \theta \text{sin} \theta. Therefore, the kinetic energy reads T = \dfrac{m}{2} (\dot x^2 + \dot y^2 + \dot z^2) and after substitution of velocities and simplifications we get T = \dfrac{mR^2}{2} (\dot \theta^2 + \omega ^2 \text{sin}^2 \theta⁡). Similarly, the potential energy of the mass particle reads V = mgz = mgR \text{cos} \theta. Note that the kinetic energy of the particle consists of those resulted from angular velocity R \dot \theta, defined in spherical coordinates in the R-z plane due to sliding of the mass on the circular ring, and the rotational velocity (R \text{sin} \theta ) \omega, defined in a z = const plane parallel to the x-y plane at any given time during the motion.

Now we can write the Lagrange’s equations, using Equation (2.12), with the assumption that no friction and non-conservative forces exist, or D=Q=0. Hence \dfrac{d}{dt} \left( \dfrac{\partial T}{\partial \dot \theta} \right) - \dfrac{\partial T}{\partial \theta} + \dfrac{\partial V}{\partial \theta} = 0. But \dfrac{\partial T}{\partial \dot \theta} = mR^2 \dot \theta, \dfrac{\partial T}{\partial \theta} = \dfrac{1}{2} mR^2 \omega ^2 \text{sin} 2 \theta and \dfrac{\partial V}{\partial \theta} = -mgR \text{sin} \theta. After substitution and rearranging the terms, we get the equation of motion for the mass particle as

    \begin{equation*} R \ddot \theta - \frac{1}{2} R \omega ^2 \text{sin} 2 \theta = g \text{sin} \theta \end{equation*}

We use 20-sim to solve these system equations. An initial angular velocity of 0.2 rad/s is applied to the mass. The script code is as follows:

parameters

real g = 9.08 {m/s2}; //grav. acc.
real R = 40 {cm}; //ring radius
real omega = 0.8 {rad/s}; // ring angular velocity

variables

real theta {rad};
real theta_dot {rad/s};
real theta_dot_dot {rad/s2};

equations

theta_dot_dot= ((1/2)*omega^2*sin (2*theta)+g*sin (theta)/R);
theta_dot = int (theta_dot_dot , 0.2);
theta = int (theta_dot , 0);

The angular displacement, velocity and acceleration are shown in Figure 2-20.

Figure 2-20 Sample results as output from 20-sim

Here is a video showing how to build and run the model for this example in 20-sim:

 

2.20       Example: An Extensible Robotic Arm Rotating in a Plane

Figure 2-21 shows a load with mass m is carried by a robotic arm in the x-y plane. The length r of the arm and its angle \theta with respect to x-axis are functions of time t, or r=r(t) and \theta = \theta (t). The damping coefficients for radial and tangential motions are b_r and b_t, respectively.

Figure 2-21 An extensible robotic arm carrying a load

The generalized coordinates (or degrees of freedom) are q=(r,\theta), and corresponding velocities are q=(\dot r, \dot \theta), for mass m. We can write the kinetic energy as T = \dfrac{1}{2} m \dot r^2 + \dfrac{1}{2} m (r \dot \theta )^2, due to radial and tangential velocities, respectively. The potential energy, with reference to the support, is V = mgr \text{sin} \theta. The damping function is D = \dfrac{1}{2} b_r \dot r^2 + \dfrac{1}{2} b_t \dot \theta ^2. The conservative gravity force due to the load mass is accounted for through the potential function V. The force f and torque \mathcal{T} exerted by the robot-arm motor to move the mass are components of generalized force vector, or Q_i = (Q_r,Q_{\theta} ) = (f, \mathcal{T}). Now, we have \dfrac{\partial T}{\partial \dot r} = m \dot r, \dfrac{\partial T}{\partial \dot \theta} = mr^2 \dot \theta, \dfrac{\partial T}{\partial r} = mr \dot \theta ^2, \dfrac{\partial T}{\partial \theta} = 0, \dfrac{\partial D}{\partial \dot r} = b_r \dot r, \dfrac{\partial D}{\partial \dot \theta} = b_t \dot \theta, \dfrac{\partial V}{\partial r} = mg \text{sin} \theta and \dfrac{\partial V}{\partial \theta} - mgr \text{cos} \theta. Using Equation (2.12), we can write the equations of the motion for the mass m, as \dfrac{d}{dt} \left( \dfrac{\partial T}{\partial \dot r} \right) -\dfrac{\partial T}{\partial q_i} + \dfrac{\partial D}{\partial \dot q_i} + \dfrac{\partial V}{\partial q_i} = Q_i

    \begin{equation*} \left\{ \begin{array} ( m \ddot r - mr \dot \theta ^2 + b_r \dot r + mg \text{sin} \theta = f \\ mr^2 \ddot \theta + 2mr \dot r \dot \theta + b_t \dot \theta + mgr \text{cos} \theta = \mathcal{T} \end{array} \end{equation*}

We use 20-sim to solve the system equations. The script code is as follows:

parameters

real m = 0.5 {kg}; // load mass
real g = 9.08 {m/s2}; //grav. acc.
real bt = 0.5 {N.m.s/rad}; // tangential damping
real br = 0.2 {N.s/m}; // radial damping

variables

real arm {m};
real arm_dot {m/s};
real arm_dot_dot {m/s2};
real theta {rad};
real theta_dot {rad/s};
real theta_dot_dot {rad/s2};
real F {N}; //applied force
real T {N.m}; // applied torque

equations

arm_dot_dot = (arm*theta_dot^2-g*sin (theta)-br*arm_dot/m+F/m);
theta_dot_dot = (1/(m*arm^2))*(-2*m*arm*arm_dot*theta_dot-m*g*arm*cos (theta)-bt*theta_dot+T);
arm_dot = int (arm_dot_dot , 0);
arm = int (arm_dot , 0.2);
theta_dot = int (theta_dot_dot , 0);
theta = int (theta_dot , 0);
F = sin (0.2*time);
T = 0.2;

Exercise Problems for Chapter 2

Exercises

  1. Using the Equation Model tool in 20-sim, build a model for the example given in section 2.13. Using the numerical data for the parameters, run simulation and analyze the results.
  2. Using the Equation Model tool in 20-sim, build a model for the example given in section 2.14. Using the numerical data for the parameters, run simulation and analyze the results.
  3. Using the Equation Model tool in 20-sim, build a model for the example given in section 2.15. Using the numerical data for the parameters, run simulation and analyze the results.
  4. Using the Equation Model tool in 20-sim, build a model for the example given in section 2.20. Using the numerical data for the parameters, run simulation and analyze the results.
  5. Using Lagrangian method, derive the system equations for the double pendulum system shown below. Solve the resulting system of ODE’s and draw the angular displacements and velocities \theta_i and \dot \theta_i, i = 1,2) of mass m_1 and m_2 for an initial condition of m_2 at \dfrac{\pi}{12}. Also draw the phase diagram (i.e., \dot \theta vs. \theta) for each mass. Assume that the strings are massless and inextensible.
  1. For the mechanical system given, e.g., an elevator with a mass-spring-damper subsystem, verify the system equations, using Lagrangian method and solve them with 20-sim. The container could be an elevator, e.g., with a mass M and is supported by a spring k_3 and moving vertically, guided by frictionless rollers under load f(t). The subsystem is composed of a mass m, two springs k_1 and k_2, and a damper b, as shown in the figure below. The gravitational acceleration vector is directed downward, \vv{g} = (g,0,0).

    \begin{equation*} \begin{bmatrix}M & 0\\ 0 & m \end{bmatrix} \begin{Bmatrix} \ddot x_1\\ \ddot x_2 \end{Bmatrix} + b \begin{bmatrix} 1 & -1\\ -1 & 1 \end{bmatrix} \begin{Bmatrix} \dot x_1\\ \dot x_2 \end{Bmatrix} + \begin{bmatrix} k_1 + k_2 + k_3 & -k_1 -k_2 \\ -k_1 -k_2 & k_1 + k_2 \end{bmatrix} \begin{Bmatrix} x_1 \\ x_2 \end{Bmatrix} = \begin{Bmatrix} f(t) - Mg\\ -mg \end{Bmatrix} \end{equation*}

  1. Repeat the sliding mass on a rotating circular ring example given in section 2.19 assuming \omega = \omega (t). Modify the model provided for this example accordingly and run the simulation.
  2. Repeat the example given in section 2.16 after adding a mechanical spring (k_3) between mass m_3 and the wall. Modify the model provided for this example accordingly and run the simulation.
  1. Repeat the example given in section 2.17 after replacing the pendulum with a double pendulum. Modify the model provided for this example accordingly and run the simulation.
  1. Derive the system equations for the electrical circuit shown in the below sketch. Use Lagrangian method and solve the resulting system of ODEs with 20-sim.

Media Attributions

Share This Book