Figure	Section
Figure 1‑1 Modern design process for a system or component	1
Figure 2‑1 Lagrangian surface visualized in space	2.1
Figure 2‑2 A sketch for visualizing Euler-Lagrange’s equation	2.1
Figure 2‑3 A frictionless mass-spring system	2.2
Figure 2‑4 Sketch for variation of for an arbitrary	2.5
Figure 2‑5 A mass-spring system with three degrees of freedom	2.11
Figure 2‑6 Sample results as output from 20-sim	2.11
Figure 2‑7 A mass-spring-damper system with two degrees of freedom	2.12
Figure 2‑8 Sample results as output from 20-sim	2.12
Figure 2‑9 A two-loop electrical circuit with source	2.13
Figure 2‑10 Sample results as output from 20-sim	2.13
Figure 2‑11 A compound Atwood’s machine.	2.14
Figure 2‑12 Atwood’s machine	2.15
Figure 2‑13 A complex vibrating mechanical system	2.16
Figure 2‑14 Sample results as output from 20-sim	2.16
Figure 2‑15 Pendulum with oscillating pivot	2.17
Figure 2‑16 Sample results as output from 20-sim	2.17
Figure 2‑17 A pendulum attached to a mass-spring-damper system	2.18
Figure 2‑18 Sample results as output from 20-sim	2.18
Figure 2‑19 A particle moving on a circular ring	2.19
Figure 2‑20 Sample results as output from 20-sim	2.19
Figure 2‑21 An extensible robotic arm carrying a load	2.20
Figure 3‑1 Sketch of a mechanical system with its components category types	3.2
Figure 3‑2 Bond graph power direction and associated effort and flow definitions, B receives power from A	3.2
Figure 3‑3 Causality assignment definition and directions of effort and flow between elements A and B	3.3
Figure 3‑4 Causality assignments for an I-element, with preferred integral causality indicated by dashed circle (left) and derivative causality (right)	3.4.1
Figure 3‑5 Block diagram (left) and equivalent bond graph for an I-element with assigned integral causality and state variable	3.4.1
Figure 3‑6 Block diagram (left) and equivalent bond graph for an I-element with assigned derivative causality	3.4.1
Figure 3‑7 Causality assignments for a C-element, with preferred one indicated by dashed circle, integral causality (right) and derivative causality (left)	3.4.2
Figure 3‑8 Block diagram (left) and equivalent bond graph for a C-element with assigned integral causality and state variable	3.4.2
Figure 3‑9 Block diagram (left) and equivalent bond graph for a C-element with assigned derivative causality	3.4.2
Figure 3‑10 Causality assignments for an R-element	3.4.3
Figure 3‑11 Block diagrams (left) and equivalent bond graph for an R-element with assigned causality	3.4.3
Figure 3‑12 Bond graph symbols for effort source (left) and flow source (right) with their assigned causalities	3.4.4
Figure 3‑13 Bond graph symbol for 1-junction element with four connecting bonds, corresponding causalities, and strong bond identified with thick half-arrow	3.4.5
Figure 3‑14 Bond graph symbol for 0-junction element with four connecting bonds, corresponding causalities, and strong bond identified with thick half-arrow	3.4.5
Figure 3‑15 Block diagrams (left) and equivalent bond graphs for a TF-element with related assigned causalities—inputs are shown with thick arrows	3.4.6
Figure 3‑16 Block diagrams (left) and equivalent bond graphs for a GY-element with related assigned causalities—inputs are shown with thick arrows	3.4.6
Figure 3‑17 Consistency of integral causality assignment and state variable for an I-element with parameter m	3.5.1
Figure 3‑18 Consistency of integral causality assignment and state variable for a C-element with parameter c	3.5.1
Figure 4‑1 A mass-spring-damper mechanical system	4.3
Figure 4‑2 Bond graph model for a one-DOF mass-spring-damper system	4.3
Figure 4‑3 Simplified bond graph model for a one-DOF mass-spring-damper system	4.3
Figure 4‑4 A two-DOF mass-spring-damper mechanical system	4.4
Figure 4-5 Bond graph model for a two-DOF mass-spring-damper mechanical system	4.4
Figure 4-6 A three-DOF mass-spring-damper mechanical system	4.5
Figure 4-7 Bond graph model for a three-DOF mass-spring-damper mechanical system	4.5
Figure 4-8 Kinetics of a one-DOF mechanical system with showing the stream of efforts with its bond graph model	4.6
Figure 4-9 Kinematics of a one-DOF mechanical system with showing the stream of flows with its bond graph model	4.6
Figure 4-10 Traditional approach for system simulation and design	4.7.1
Figure 4-11 Modern approach for system simulation and design	4.7.2
Figure 5‑1 The 20-sim Editor interface	5.1
Figure 5‑2 The 20-sim Simulator interface	5.1
Figure 5‑3 Process steps for design of a system using the modeling, simulation and analysis	5.1
Figure 5‑4 Sketch for a car seat mechanical system	5.2
Figure 5‑5 Bode graph model for a car seat	5.2
Figure 5‑6 Sketch for a cart system carrying a load	5.3
Figure 5‑7 Bond graph model for the cart carrying a load	5.3
Figure 6‑1 Decomposition of 2D rigid-body motion into translation and rotation	6.3
Figure 6‑2 Gear and shaft mechanical system sketch	6.4
Figure 6‑3 Bond graph model for a gear shaft system, built in 20-sim	6.4
Figure 6‑4 Bond graph model for a gear shaft system, derivative causalities removed	6.4
Figure 6‑5 A double pinion-rack mechanical system	6.5
Figure 6‑6 Bond graph model for the double pinion-rack system	6.5
Figure 6‑7 A mass-spring-damper system on an inclined plane	6.6
Figure 6‑8 Bond graph model for the mass-spring-damper system on an inclined plane	6.6
Figure 6‑9 Half-car mechanical system sketch	6.7
Figure 6‑10 Bond graph model for a half-car mechanical system	6.7
Figure 6‑11 A mass-spring-damper mechanical system attached to a lever	6.8
Figure 6‑12 A mass-spring-damper mechanical system attached to a beam	6.9
Figure 6‑13 Two moving mass-spring system attached to a lever	6.10
Figure 6‑14 A two-pulley mechanical system	6.11
Figure 7‑1 Sign convention for electrical current through passive elements, passive sign convention	7.2
Figure 7‑2 Electrical power sign for several elements according to passive sign convention	7.2
Figure 7‑3 Sketch for a RCL electrical circuit in series	7.4
Figure 7‑4 Bond graph model for a RCL electrical circuit in series	7.4
Figure 7‑5 Simplified bond graph model for a RCL electrical circuit in series	7.4
Figure 7‑6 Sketch for a RCL electrical circuit in parallel	7.5
Figure 7‑7 Bond graph model for a RCL electrical circuit in parallel	7.5
Figure 7‑8 A two-loop RCL electrical circuit	7.6
Figure 7‑9 Bond graph model for the two-loop RCL electrical circuit	7.6
Figure 7‑10 A three-loop electrical circuit	7.7
Figure 7‑11 Bond graph model for the three-loop RCL electrical circuit	7.7
Figure 7‑12 A Wheatstone bridge electrical circuit	7.8
Figure 7‑13 Bond graph model for the Wheatstone bridge circuit	7.8
Figure 7‑14 A multi-loop electrical circuit	7.9
Figure 7‑15 Bond graph model for the multi-loop electrical circuit	7.9
Figure 7‑16 A multi-loop electrical circuit with transformer	7.10
Figure 7‑17 Bond graph model for the multi-loop electric circuit with transformer	7.10
Figure 8‑1 Sketches for pressure drop in a pipe and a storage tank	8.1
Figure 8‑2 Sketch for a control volume of flowing fluid in a pipe	8.4
Figure 8‑3 Sketch of velocity profile for a Hagen-Poiseuille flow in a pipe	8.5
Figure 8‑4 Sketch for a two-tank hydraulic system	8.8
Figure 8‑5 Bond graph model for a two-tank hydraulic system	8.8
Figure 8‑6 A hydraulic system with pump	8.10
Figure 8‑7 Bond graph model for the pump-reservoir-valve hydraulic system	8.10
Figure 8‑8 A hydraulic lift system	8.11
Figure 8‑9 Bond graph model for hydraulic lift system	8.11
Figure 9‑1 A car brake hydro-mechanical system	9.2
Figure 9‑2 Bond graph model for the car brake hydro-mechanical system	9.2
Figure 9‑3 An Electro-mechanical system with load	9.3
Figure 9‑4 Bond graph model for an electro-mechanical system	9.3
Figure 9‑5 A drive shaft mechanical system carrying a torsional load	9.4
Figure 9‑6 Bond graph model for the drive shaft mechanical system carrying a torsional load	9.4
Figure 9‑7 An inverted double pendulum system	9.5
Figure 9‑8 Bond graph model for the inverted double pendulum	9.5
Figure 10‑1 Linear system sketch for processing inputs and outputs	10.3
Figure 10‑2 Linear System Editor interface in 20-sim	10.4
Figure 10‑3 Transfer Function Editor interface in 20-sim	10.4
Figure 10‑4 Bode plots for a PI controller	10.4
Figure 10‑5 Model Linearization interface in 20-sim	10.4.1
Figure 10‑6 Typical Bode plots for a system	10.4.1
Figure 10‑7 Bode plots for the transfer function	10.5
Figure 10‑8 Mechanical system sketch for Example given in section 10-6	10.6
Figure 10‑9 Bode plots for mechanical system given in section 10-6	10.6
Figure 11‑1 Implementing gravity force for a bond graph model, in 20-sim	11.3
Figure 11‑2 A mechanical system sketch	11.4.1
Figure 11‑3 Bond graph model for the mechanical system with labelled power bonds	11.4.1
Figure 11‑4 Bond graph model for the electrical system with labelled power bonds	11.4.2
Figure 11‑5 A b ond graph model with derivative causality, colour-coded	11.5.1
Figure 11‑6 A b ond graph model with algebraic loop causality	11.5.2
Figure 11‑7 The b ond graph model with removed algebraic loop-selecting R3	11.5.2
Figure 11‑8 The b ond graph model with removed algebraic loop-selecting R2	11.5.2
Figure 11‑9 The b ond graph model with removed algebraic loop-selecting R1	11.5.2

Table	Section
Table 3‑1 Analogous quantity definitions in BG modelling method related to different systems	3.2
Table 3‑2 Typical mechanical components and their corresponding BG elements	3.4
Table 5‑1 Data for car seat example and BG model	5.2
Table 5‑2 Data for cart with load example and BG model	5.3
Table 6‑1 Typical mechanical rotational components and their corresponding BG elements	6.2
Table 6‑2 Data for the gear-shaft mechanical system	6.4
Table 7‑1 Typical electrical components and their corresponding BG elements	7.1
Table 8‑1 Typical hydraulic components and their corresponding BG elements	8.5
Table 8‑2 Data for system components in example 8-9.	8.9
Table 10‑1 Response of a PI controller to a sinusoidal signal, $/sin /omega t$	10.3
Table 11‑1 Inputs to the mechanical system by all elements- Q1 question.	11.4.1
Table 11‑2 Outputs from the mechanical system to all storage elements- Q2 question.	11.4.1
Table 11‑3 Inputs to the electrical system by all elements- Q1 question	11.4.2
Table 11‑4 Outputs from the electrical system to all storage elements- Q2 question	11.4.2
Table 11‑5 Definitions for effort and flow for BG and pseudo BG models, thermal systems	11.6

Introduction

The foundations of engineering practice are mathematical models, the principles of physics, and empirical results obtained from experiments for defining design criteria. An engineer must know the laws of physics very well and use the relevant mathematical models and their solutions, either exact or numerical, to design parts, systems, and complex machines which function with certain reliability for an assumed lifetime. To help with this task, an engineer may use modelling tools to simulate the behavior of systems and their components. Modelling and the application of software tools are becoming increasingly common in modern engineering practice. As shown in Figure 1‑1, modeling and simulation results can help optimize and refine a design before the physical prototype is built. This minimizes the time required for the design process. In addition, application of modelling can minimize the final cost of a prototype or a product.

Figure 1-1 Modern design process for a system or component

Isaac Newton (1643–1727)

Modelling has a long history starting from ancient times when scientists used “equations” to relate variables or parameters to one another (e.g., Archimedes, Thales, Khawrazmi). Later, scientists and mathematicians developed “equations” which could represent the way that natural phenomena work and materials behave. These “equations” are sometimes referred to as laws of physics and constitutive equations since they are validated through time and the obtained results match with what we experience or measure in the real world with some approximations, of course. For example, Newton’s second law is given as a model which predicts the behaviour of material bodies under given forces applied to them, i.e., the relationship between forces applied to a body mass and the change of its momentum with respect to time.

Similarly, Ohm’s law is a model which relates the voltage across a resistor to the electrical current using the resistor’s material property. These models, and many other similar ones (e.g., Hooke’s, Fick’s, Fourier’s) related to different engineering disciplines, form the foundation of engineering. It is through their application that we trust the behavior and responses of our designs in the real world. Assume that we are flying in an airplane which is designed based on laws and governing equations or models applied to fluid mechanics and solid mechanics, among others. If we don’t trust and accept these laws and models, then it would not be logical to ride in an airplane!

Real-world phenomena are complex and usually involve many types of physics. For application in engineering, we usually simplify these phenomena and consider the dominant physics involved. For example, the length of a simple spring linearly changes under a given load according to Hooke’s law. But it becomes a more complex problem if the spring’s material behaves non-linearly, or if for example, electrical charges flow through it. Traditionally, the simplification of a problem is/was due to lack of tools for finding a solution which could represent more accurately that problem’s real world behaviour. It is at this point that modelling methods, e.g., Lagrangian and BG, and advanced modelling software tools, e.g., 20-sim, are valuable resources for finding solutions to complex engineering systems and optimizing our designs to have more economical, reliable, and durable products as end results. Although this book focuses on using bond graphs as a modelling method, we also emphasize the importance of learning and, hence, understanding the foundation and mathematics behind an energy-based approach for system analysis. For this purpose, we summarize Lagrangian mechanics in chapter 2 and provide some references for further reading.

The main body of the text is devoted to the BG method. This graphical (i.e., it can be sketched similar to engineering drawings) method translates the physical laws relevant to a desired system at hand into graphical interactions of interconnected assigned elements. The method uses laws of thermodynamics and the principle of cause and effect (in an acausal way) with the inclusion of constitutive relations relevant to system components.

Henry Paynter (1923–2002). Courtesy MIT Museum.

In 1959, Henry M. Paynter at the MIT Department of Mechanical Engineering developed the bond graph method [1]. This method has fluctuated in application and popularity in the industry, with a recent rise due to its strength in modelling multi-energy-domain systems and the widespread availability of economically viable computer power [2].

In this book, we make use of facilities available in 20-sim, as a software tool for building, among others, BG models. 20-sim also offers solvers for finding solutions for the resulting system equations for simulation and design of systems. We use these solvers, with the modern script language SIDOPS++ included, to solve system equations as ordinary differential equations (ODEs). The system equations could be extracted from BG models or using Lagrangian method. The script language SIDOPS++ is suitable for complex system modelling and solving the relevant equations [3].

The reader may come across or already be familiar with other available methods/tools for modelling engineering systems, including block diagram, a signal processing graphical method; icon-diagram, a component-iconic graphical method; and advanced script languages/tools, e.g., Dymola, Smile, and recently Modelica [4], [5], [6].

Media Attributions

Isaac Newton © J. MacArdell after E. Seeman is licensed under a CC BY (Attribution) license
Henry Paynter © MIT Museum is licensed under a All Rights Reserved license

Adapted and modified, with permission from Mercury Learning and Information LLC.

Acausal method, like bond graph, allows the user to select input and ouput ports, in contrast to causal method, for which the ports are fixed in terms of input and output signals, e.g., block diagram method. Acausal methods can be interpreted as two-way streets vs. causal methods as one-way streets.

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)$ .

After expanding the expression in the bracket on the R.H.S., we get:

$T=\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_1 \dot x_1^2 + m_2 \dot x_2^2 + m_3 \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_1 \dot x_1$ 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_3$ . 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:

One or more interactive elements has been excluded from this version of the text. You can view them online here: https://pressbooks.bccampus.ca/engineeringsystems/?p=58#oembed-7

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:

One or more interactive elements has been excluded from this version of the text. You can view them online here: https://pressbooks.bccampus.ca/engineeringsystems/?p=58#oembed-8

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_i} \right) - \dfrac{\partial (T-V)}{\partial q_i} + \dfrac{\partial D}{\partial \dot 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.81 {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.81 {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 ẋ_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:

One or more interactive elements has been excluded from this version of the text. You can view them online here: https://pressbooks.bccampus.ca/engineeringsystems/?p=58#oembed-9

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.81 {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:

One or more interactive elements has been excluded from this version of the text. You can view them online here: https://pressbooks.bccampus.ca/engineeringsystems/?p=58#oembed-10

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 system 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 + b_t \dot \theta = 0$ . After simplification, we get $\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$ . 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 D}{\partial \dot x} = F(t)$ , or $\dfrac{d}{dt} (m \dot x + ml \dot \theta \text{cos} \theta + M \dot x) +kx + b \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.81 {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:

One or more interactive elements has been excluded from this version of the text. You can view them online here: https://pressbooks.bccampus.ca/engineeringsystems/?p=58#oembed-11

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.81 {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.0);
theta = int (theta_dot , 0.2);

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:

One or more interactive elements has been excluded from this version of the text. You can view them online here: https://pressbooks.bccampus.ca/engineeringsystems/?p=58#oembed-12

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.81 {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

Derive Lagrange equation for the system given in example 2.13. Using the Equation Model tool in 20-sim, build a model for this example. Using the numerical data for the parameters, run simulation and analyze the results.
Derive Lagrange equation for the system given in example 2.14. Using the Equation Model tool in 20-sim, build a model for this example. Using the numerical data for the parameters, run simulation and analyze the results.
Derive Lagrange equation for the system given in example 2.15. Using the Equation Model tool in 20-sim, build a model for this example. Using the numerical data for the parameters, run simulation and analyze the results.
Derive Lagrange equation for the system given in example 2.18. Using the Equation Model tool in 20-sim, build a model for this example. Using the numerical data for the parameters, run simulation and analyze the results.
Derive Lagrange equation for the system given in example 2.20. Using the Equation Model tool in 20-sim, build a model for this example. Using the numerical data for the parameters, run simulation and analyze the results.
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.

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*}$

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.
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.

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.

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

Joseph Louis Lagrange © Zephirin Belliard is licensed under a CC BY (Attribution) license
William Rowan Hamilton adapted by Quikbik is licensed under a Public Domain license
Jean le Rond d’Alembert is licensed under a Public Domain license
fig 2-11_edit

Bond Graph Modelling Method

3.1 Overview

All engineering systems share the physical phenomenon of the transfer and distribution of energy among their corresponding components while converting one form of energy to another. The balance of energy “flowing” through a system should be maintained. The total amount of energy remains constant—energy is conserved—according to the first law of thermodynamics. In 1959, Henry Paynter used the first law and common system features to create a general graphical method for analyzing and modelling multi-domain engineering systems. His objectives were mainly to have a unified graphical method for modelling single- and multi-domain systems as well as a common procedural algorithm to develop such models and obtain their relevant systems’ equations. Hence, the bond graph (BG) method was created [1], [20], [21], [22].

BG represents a system through graphical modelling. The BG method assigns ports (the communication point) for each component of a system and connects each port to the adjacent component through bonds (the communication path and direction) for a two-way energy/power exchange. At any instant of time, each component either receives (sends) a quantity called effort and simultaneously sends (receives) another quantity called flow. The product of the quantities of effort and flow has the dimension of power—or time rate of energy change. In a mechanical system, force is the effort and velocity is the flow; in an electrical system, voltage is the effort and current is the flow. The collection of bonds—with the inclusion of the related system components’ constitutive laws, constraints, and boundary conditions—forms the system BG model. Building a BG model requires nine basic elements, defined as follows. See section 3.4 for full description.

The nine basic BG elements, along with the principle of causality, can be employed for building a BG model representing a given system’s dynamical behaviour (see further sections for detailed explanation). The resulting BG model, then, would clearly show the kinematic (i.e., continuous stream of flow) and kinetic (i.e., continuous stream of effort) of the system and can be used to extract the equations governing the dynamical behaviour of the whole system. In addition, the insights provided by a BG model are valuable for understanding the physics/dynamics of the system and provide a powerful tool for simulation, design, and optimization of the system. The procedure for building a BG model is similar for analogous engineering systems. For example, when the parameters of the pertinent components are used, the BG model for a mechanical mass-spring-damper system is identical to those of an electrical resistor-capacitor-inductor (RCL) system.

In this chapter, we discuss, among other topics, the definitions for basic BG elements, the causality principle and assignments, and the concept of state variables.

3.2 Categorizing System Components—Generalized Effort and Flow

The components of a system can be categorized according to energy transfer through the system into three types. These are kinetic energy storages, potential energy storages, and energy dissipaters. In addition, we have energy source/sink components acting with the surroundings at the boundary of the system. There may also exist components that simply transfer energy without storing or dissipating it. Finally, a system may include components, such as a distributor, that perform as junctions. Figure 3‑1 shows a sketch of a mechanical system with examples of component categories, as mentioned above. All these types of components can be modelled using nine basic BG elements, as discussed in further sections.

Figure 3-1 Sketch of a mechanical system with components’ category types

The dynamical behaviour of a system changes with time. Therefore, time rate of energy or power is the quantity of interest in BG models. The relation between energy $\mathcal{E}$ and power $\mathcal{P}$ can mathematically be written as

(3.1) $\begin{equation*} \mathcal{P} (t) = \frac{d \mathcal{E} (t)}{dt} \quad or \quad \mathcal{E} (t) = \int \mathcal{P} (t)dt \end{equation*}$

We can identify the components of a given system as “lumped” entities that exchange energy with one another. Using the first law of thermodynamics, we can write the change in energy as the sum of work $\mathcal{W}$ and heat $\mathcal{Q}$ exchanges, or $d \mathcal{E} = \mathcal{P} dt = \delta \mathcal{W} + \delta \mathcal{Q}$ . Summing up the energy changes of lumped components in a system gives the total energy change of the system. For example, without losing generality, we consider a mechanical system component receiving power and exhibiting a displacement $x(t)$ and velocity $v(t) = \dfrac{dx}{dt} = \dot x(t)$ . Using Equation (3.1), the amount of energy in terms of work input can be written as $\mathcal{W} = \int \mathcal{P} dt$ . But the work is also equal to the force times the displacement; hence, $\mathcal{W} = \int F(t)dx$ . Substituting for $dx = v(t)dt$ , we get $\mathcal{W} = \int Fvdt$ . Considering $t_1$ and $t_2$ as the time limits associated with the duration of energy transfer, we can, after integrating, write the work as

(3.2) $\begin{equation*} \begin{dcases} \Delta \mathcal{W} |_\textit{potential} = \int_{t_1}^{t_2} F \left( \underbrace{vdt}_{dx} \right) = Fx|_{t_2}-Fx|_{t_1} \\ \Delta \mathcal{W} |_\textit{kinetic} = \int_{t_1}^{t_2} v \left( \underbrace{Fdt}_{dp} \right) = \int_{t_1}^{t_2} v(mdv) = \frac{m}{2} (v^2 |_{t_2} - v^2 |_{t_1}) \\ \Delta \mathcal{W} |_\textit{dissipation} = \int_{t_1}^{t_2} (Fv)dt = \mathcal{P}|_{t_2} - \mathcal{P}|_{t_1} \end{dcases} \end{equation*}$

where $m$ is mass. In a BG model, each system component is designated by a suitable basic element and associated port(s). Depending on the type of element used, the number of ports could be one, two, or more. The power direction is designated by a half-arrow ( $\longrightharpoondown$ ) which shows the direction of power to or from the port for each element. Traditionally, half-arrows are used in BG models to keep the full-arrow shape for one-way signal data, as in block diagram graphs.

As mentioned above and by Equation (3.2), for mechanical systems, the power is composed of two quantities: force and velocity. In BG method, we generalize this concept and show the power with the product of $e$ and $f$ , the effort and flow, respectively. Hence, the product of effort and flow has the dimension of power, or $\textit{power} = e \cdot f$ . For example, for a rotational motion, $e$ is the torque and $f$ is the angular velocity (see Table 3‑1). In other words, in a BG model, the kinetics of a system is modelled by transfer of the efforts of its components according to the equilibrium, and the kinematics by transfer of components’ flows according to compatibility requirement. We will discuss this feature of BG method, using some examples, in section 4.6.

Figure 3‑2 shows the definition of power direction for element A sending power to B, and the associated effort, shown above the half-arrow, and flow, shown, by definition, below the half-arrow.

Figure 3-2 BG power direction and associated effort and flow definitions: B receives power from A

Table 3‑1 Analogous quantity definitions in BG modelling method related to different systems
Systems	Effort ( $e$ )	Flow ( $f$ )	Displacement $q = \int fdt$	Momentum $p= \int edt$
mechanical-translational	force [N]	velocity [m/s]	distance [m]	[kg.m/s]
rotational mechanical	torque [N.m]	angular velocity [rad/s]	angle [rad]	angular momentum [kg.m2/s]
hydraulic	pressure [Pa]	volume flow rate [m³/s]	volume [m³]	hydraulic momentum [Pa.s]
thermal/thermodynamics	temperature [K]	entropy change rate [J/ K.s]	entropy [J/K]	—
thermo-fluid	enthalpy (specific) [J]	mass flow rate [kg/s]	mass flow [kg]	flow momentum
electrical	voltage [V]	current [A]	charge [C]	flux linkage [V.s]
magnetics	magnetic force [A]	magnetic flux rate [Wb/s]	magnetic flux [Wb]	—
chemical	chemical potential [J/mol]	mole flow rate [mol/s]	mole flow [mol]	—

3.3 Causality Principle and Assignment

To establish the principle of cause and effect relationship in BG method, we use the definition of causality assignment. The cause signal brings all the history data to the system/element, and through the dynamical behaviour of the system, the present signal effect is decided and provided as output.

As mentioned, in BG method, the half-arrow indicates the direction of power between related elements in a BG model. However, the half-arrow does not provide information about the direction of power constituents, i.e., effort or of the flow. In principle, we can arbitrarily define these directions. For example, in Figure 3‑2, we can assign $e$ direction pointing from component A to B (hence, $f$ should be directing from B to A) or vice versa. In other words, the causality assignment is a symmetrical one. By definition, a small transverse/vertical line, a causality stroke, is drawn close to one of the ports at the power bond to show the direction of effort toward it, hence the direction of flow away from it, as shown in Figure 3‑3. This operation is critical for building BG models and, in terms of providing a definite solution, has consequences in the resulting equations of the system. After the causality is assigned, then the signal received by the element is the cause, and the returning signal—or the element response—is the effect.

The preferred causality assignment is called integral causality, and the alternative option is the derivative/differentiate causality. We will discuss the details further in section 3.5.

Figure 3-3 Causality assignment definition and directions of effort and flow between elements A and B

3.4 Nine Basic Elements of Bond Graph Method

As mentioned in the previous section, building a BG model of a physical system involves consideration of the energy conservation, transfer, and conversion through the system. In a BG model, we focus on the rate of energy or power as the quantity to deal with.

For energy storage, we define two elements, represented by $I$ (inertial element) for kinetic energy and $C$ (capacity element) for potential energy storages. For energy dissipation, we define one element, represented by $R$ (friction or resistor element). We represent the energy source/sink acting at the boundary of the system by two elements, one for effort $S_e$ and one for flow $S_f$ . To manage the distribution of energy through the system, we define two elements as junctions, represented by junction 1 and junction 0. For energy transfer/conversion, we define two elements, represented by transformer $TF$ and gyrator $GY$ . Therefore, in total, we have nine elements available and sufficient for building a BG for any given physical system, with the inclusion of their modulated versions ( $MS_e$ , $MS_f$ , etc.) for when a signal is input to the corresponding element from an external source. Examples of physical/engineering systems are mechanical, electrical, thermal, hydraulic systems, or some hybrid systems composed of subsystems assembled of different energy media.

Portrait of a Mathematician, thought to be of Robert Hooke (1635–1703)

Each one of the BG elements mentioned above should behave according to the relevant physical laws represented by their constitutive relations—a mathematical model. For example, a linear mechanical spring is modelled by element $C$ , whose governing equation should comply with Hooke’s law. However, a given spring can go under deformation either by receiving an effort (i.e., force) or a flow (i.e., displacement rate/velocity). Depending on the system and computational preferences, we can assign causality strokes to the element $C$ to specify that the desired spring receives effort or flow. This rule, the causality assignment, must be applied to all bonds in a BG model. Examples of typical translational mechanical elements are shown in Table 3‑2.

Table 3‑2 Typical mechanical components and their corresponding BG elements
$R$ -element (damper)	$C$ -element (spring)	$I$ -element (mass)	$TF$ -element (lever)	$GY$ -element (gyroscope)

In the next sections, we will define the constitutive equations, preferred causality, and physical representation examples for all nine BG elements.

3.4.1 Inertia Element I: Kinetic Energy Storage

In BG modelling, the $I$ -element is a passive element; it should receive power to return a signal. This requirement means that the half-arrow power bond should be drawn toward this element. An $I$ -element has only one port for communicating to the rest of the system. Examples are mass bodies in mechanical systems and inductors in electrical systems.

As shown in Figure 3‑4, the input quantity for the $I$ -element can be either effort ( $e$ ) or flow ( $f$ ); consequently, the response is flow or effort, respectively. Note that the causality stroke (the vertical/transverse line) specifies the direction of effort defined to be toward the stroke; hence, the direction of flow is to be away from it. We use red colour for specifying non-integral causality strokes.

Figure 3-4 Causality assignments for an $I$ -element, with preferred integral causality indicated by dashed circle (left) and derivative causality (right)

Now the question is, how do we choose between these two possible options when building a BG model? What are the implications when choosing one option versus the other? The short answer is that both options are legitimate, but there is a preference for having the $I$ -element receiving the effort and sending the flow out—integral causality—hence, the causality stroke is placed at the half-arrow head at the port close to the element. The effort is the cause, and the flow is the effect relevant to $I$ -element when it is integrally causalled.

In a bond graph model for an $I$ -element, the preferred causality assignment is effort-in, so-called integral causality.

Mathematically, the statement given in the box can be analyzed as follows. In a mechanical system, for example, we consider a point mass $m$ and apply Newton’s second law to the motion of that point mass. Therefore, we can write ${F}={m}\dfrac{dv}{dt}={m}\dot{v}$ ( $F$ is net applied force, and $v$ is the velocity of the mass), or in BG generalized notation, ${e}={I}\dfrac{df}{dt}={I}\dot{f}$ . Recall that the $e$ symbol represents effort (force) and $f$ represents flow (velocity) in a mechanical system (see Table 3‑1). We also use the symbol $I$ , representing mass $m$ , or inductance for electrical systems. Now, for the effort-in option that we have, since the input should be $e$ , $e ={I}\dot{f}\Rightarrow\underbrace{\dot f}_{output}=\underbrace{\dfrac{e}{I}}_{input}$ or after integration, $\int\dot{f}{dt}=\dfrac{1}{I}\underbrace{\smallint e\:dt}_{momentum}\Rightarrow{f}=\dfrac{p}{I}.$ Note that the integral of force with respect to time is the momentum $p$ , ( $\smallint{F}{dt}=\smallint{edt}$ ). The equation ${f}=\dfrac{p}{I}$ is the key point here. Let’s see what it means. The flow (velocity) is equal to momentum divided by the mass. This is well-known! In BG method, however, it has an important meaning: for the $I$ -element, the input effort quantity, after being integrated, is divided by the $I$ -element parameter $m$ , and the output quantity is flow or velocity. This can be shown in a block/signal diagram along with equivalent BG model diagram (see Figure 3‑5). Since the integration of effort is involved, we call the related causality assignment an integral causality which is preferred for $I$ -elements. From the physical point of view, the integration of effort collects all the input data and hence represents a more comprehensive description of the system in terms of modelling. In addition, the resulting system’s equations (see section 3.5) are first-order ODEs when integral causality is assigned.

Figure 3-5 Block diagram (left) and equivalent bond graph for $I$ -element with assigned integral causality and state variable $p$

The constitutive equation for the $I$ -element in a BG model is given as

(3.3) $\begin{equation*} f={\frac{p}{I}}={\frac{1}{I}}\underbrace{\int e\:dt}_{state\:var.} \end{equation*}$

The momentum $p$ , which is the result of input/effort integration, is a state variable (see section 3.5).

Now, we consider the second possible option with flow-in signal (see Figure 3‑4). We have $e=I\dot{f}=I\dfrac{df}{dt}$ . This equation matches with the input and output data, since the time derivative of input flow, given on the right-hand side of the relation, multiplied by the $I$ -element parameter is the element output or effort, given on the left-hand side. This is the derivative causality assignment since the derivative of input data is involved. This case can be shown in a block diagram along with equivalent BG model diagram (see Figure 3‑6).

Figure 3-6 Block diagram (left) and equivalent bond graph for an $I$ -element with assigned derivative causality

3.4.2 Capacity Element C: Potential Energy Storage Element

In BG modelling method, the $C$ -element is a passive element because it should receive power to react to. This requirement means that the half-arrow power bond should be drawn toward this element. A $C$ -element has only one port for communicating to the rest of the system. Examples are springs in mechanical and capacitors in electrical systems. As shown in Figure 3‑7, the input quantity can be either effort ( $e$ ) or flow ( $f$ ); consequently, the response is flow or effort, respectively. Note that the causality stroke (the vertical line) specifies the direction of effort defined to be toward the stroke; hence, the direction of flow is to be away from it.

Figure 3-7 Causality assignments for a $C$ -element, with preferred one indicated by dashed circle, integral causality (right) and derivative causality (left)

Now the question is, how do we choose between these two possible options when building a model? What are the implications when choosing one option versus the other? The short answer is that both options are legitimate, but there is a preference for having the $C$ -element sending the effort and receiving the flow—integral causality—hence, the causality stroke is placed at the opposite end of the half-arrow head away from the element’s port.

In a bond graph model for a $C$ -element the preferred causality assignment is effort-out, so-called integral causality.

Mathematically, the statement given in the box can be analyzed as follows. In a mechanical system, e.g., we consider a linear mechanical spring with stiffness $k$ and apply Hooke’s law to its motion. Therefore, we can write $F=kx$ ( $F$ is net applied force, and $x$ is the displacement) or, in generalized BG notation, $e=kq=\dfrac{q}{c}$ , where $c=1/k$ the spring compliance. Recall that the $e$ symbol represents effort (force) and $f$ represents flow (velocity) in, e.g., a mechanical system, (see Table 3‑1). We use the symbol $c$ , representing spring compliance or capacitance in electrical systems as well.

Now, for the effort-out option having the flow $f$ as the input, we can write $e=\dfrac{q}{c}=\dfrac{1}{c}{\int f\:dt}$ . That is, for the $C$ -element, the input flow quantity, after integration, is divided by the $C$ -element’s compliance and gives the output quantity as effort $e$ . For $C$ -element, the displacement $q$ , which is the result of input/flow integration, is the state variable.

(3.4) $\begin{equation*} e=\frac{q}{c}=\frac{1}{c}\underbrace{\int{f}\:dt}_{state\:var.} \end{equation*}$

This can be shown in a block/signal diagram along with equivalent BG model diagram (see Figure 3‑8).

Figure 3-8 Block diagram (left) and equivalent bond graph for a $C$ -element with assigned integral causality and state variable $q$

Now, we consider the second possible option with effort-in signal (see Figure 3‑7). We can write $q=ce$ , with effort being the input and displacement as the output data; hence, the time derivative of output displacement ( $q$ ) is required to get the flow/velocity. This is derivative causality assignment, since the derivative/differential operation is needed to get the output signal involved. This case can be shown in a block diagram along with equivalent BG model diagram (see Figure 3‑9).

Figure 3-9 Block diagram (left) and equivalent bond graph for a $C$ -element with assigned derivative causality

3.4.3 Friction Element R: Energy Dissipation Element

In BG modelling method, the $R$ -element is a passive element since it should receive power to return a signal. This requirement means that the half-arrow power bond should be drawn toward this element. An $R$ -element has only one port for communicating to the rest of the system. Examples are dampers in mechanical and resistors in electrical systems.

As Figure 3‑10 shows, the input quantity for the $R$ -element can be either effort ( $e$ ) or flow ( $f$ ); consequently, the response is flow or effort, respectively. Note that the causality stroke (the vertical line) specifies the direction of effort defined to be toward the stroke; hence, the direction of flow is to be away from it.

Figure 3-10 Causality assignments for an $R$ -element

There is no preference for having the $R$ -element receiving the effort or the flow. Therefore, the causality stroke can be placed at either end of the half-arrow power connection, according to the causality requirement for the adjacent elements.

In a bond graph model for an $R$ -element, there is no preferred causality assignment- i.e., it is neutrally causalled.

Mathematically, the statement given in the box can be analyzed as follows. In a mechanical system, for example, we consider a damper with viscous damping coefficient $R$ . The constitutive equation gives the force applied on the damper proportional to the rate of displacement. Hence, we can write ${F}=Rv$ ( $F$ is net applied force, and $v$ is the velocity). Writing in BG generalized notation, $e=Rf$ . Now, for the effort-in option we have, since the input should be $e$ , $\underbrace{f}_{output}=\underbrace{\dfrac{e}{R}}_{input}$

Now, we consider the option with flow-in data (see Figure 3‑10). We have $e=Rf$ . Since the constitutive equation for a linear viscous damper is algebraic, we do not need to integrate or differentiate the input signal to obtain the output signal for an $R$ -element. Therefore, there is no preference, and $R$ -element is neutrally causalled. Figure 3‑11 shows block diagrams along with equivalent BG model diagram with causality assignments for an $R$ -element.

Figure 3-11 Block diagrams (left) and equivalent bond graph for an $R$ -element with assigned causality

3.4.4 Source Elements S_e and S_f : System Boundary Input Elements

In BG modelling method, the boundary source elements are of two types. The sources for effort (such as force, voltage) and flow (such as velocity, current) are represented by $S_e$ and $S_f$ respectively. These elements are active, and the half-arrow power bond should be drawn from these sources to the connecting elements in the system. Source elements have only one port each, for communicating to the rest of the system. As shown in Figure 3‑12, the causality assignments are uniquely assigned for these elements.

Figure 3‑12 BG symbols for effort source (left) and flow source (right) with their assigned causalities

3.4.5 1- and 0-junctions: Distribution Constraint Elements

In BG modelling method, system-required constraints for distribution of energy are applied using two elements. These are multi-port elements with symbols “1” and “0” that can receive or send power to the elements connecting to them. This requirement means that the half-arrow power bond can be drawn toward or from these elements.

A 1-junction is a flow equalizer or an effort summator element. For example, in a mechanical system, a common node with connecting system components exhibits the same value of velocity, or the elements of an electrical circuit experience the same value of current. The causality assignment for a 1-junction element must comply with its definition of distributing the flow received through one of the connecting bonds to the rest of bonds. Therefore, only one bond can send flow to a 1-junction—the strong bond—and the remaining connecting bonds should send the same flow to connecting elements; hence, the causalities are assigned accordingly, as shown in Figure 3‑13.

Figure 3‑13 BG symbol for 1-junction element with four connecting bonds, corresponding causalities, and strong bond identified with thick half-arrow

After labelling the bonds with arbitrary numbers, we can write the conservation energy law, in terms of its rate, as $e_1 f_1 + e_2 f_2 - e_3 f_3 - e_4 f_4 = 0$ . But the 1-junction distributes the flow received from the strong bond (i.e., the bond labelled “1”) equally to bonds 2, 3, and 4. Hence, $f_1 = f_2 = f_3 = f_4$ . From these relations, after substitution, we get $e_1 + e_2 - e_3 - e_4 = 0$ . Similarly, for $N$ number of bonds connecting to a 1-junction, we have the constraint relations for the 1-junction as

(3.5) $\begin{equation*} \begin{dcases} f_1 = f_2 = \cdots = f_N \\ \sum_{i=1}^{N} e_i= 0 \end{dcases} \end{equation*}$

In Equation (3.5), the summation for efforts received by 1-junction is algebraic, or the input power is considered to be positive, and the output power has a negative sign.

A 0-junction is an effort equalizer or a flow summator element. For example, in a mechanical system, a common node with connecting system components experiences relative velocity values or the nodes in an electrical circuit with common voltage. The causality assignment for a 0-junction element must comply with its definition of distributing the effort received through one of the connecting bonds to the rest of bonds. Therefore, only one bond can send effort to a 0-junction—the strong bond—and the remaining connecting bonds should send the same effort to connecting elements, hence, the causalities are assigned accordingly, as shown in Figure 3‑14.

Figure 3‑14 BG symbol for a 0-junction element with four connecting bonds, corresponding causalities, and strong bond identified with thick half-arrow

After labelling the bonds with arbitrary numbers, we can write the conservation energy law, in terms of power or energy rate, as $e_1f_1+e_2f_2-e_3f_3-e_4f_4=0$ . But the 0-junction distributes the effort received from the strong bond (i.e., the bond labelled “1”) equally to bonds 2, 3, and 4. Hence, $e_1=e_2=e_3=e_4$ . From these relations, after substitution, we get $f_1+f_2-f_3-f_4=0$ . Similarly, for $N$ number of bonds connecting to a 0-junction, we have the constraint relations for the 0-junction as

(3.6) $\begin{equation*} \begin{dcases} e_1=e_2=\cdots=e_N\\ \sum_{i=1}^{N} f_i=0 \end{dcases} \end{equation*}$

In Equation (3.6), the summation for flows received by the 0-junction is algebraic, or the input power is considered to be positive and the output power has a negative sign.

3.4.6 Transformer TF and Gyrator GY: Energy Conversion Elements

In physical engineering systems, energy may be converted by some components while its conservation is maintained. Examples are levers and gearbox in mechanical systems or electrical transformers and motors in electrical systems. In BG modelling method, there exist two elements for modelling convertors: transformer $TF$ and gyrator $GY$ . These elements are two-port elements and can receive power through one of their ports as input and deliver a converted power from the other port as output, in terms of the power variables effort and flow. The causality assignments determine the directions of flows and efforts as being inputs or outputs. In this section, we present the details of $TF$ -element followed by those of $GY$ -element.

A transformer element, $TF$ represents the converter that receives the same type of physical quantity as the type it delivers, after conversion. For example, a force applied at one end of a lever is converted to a magnified/reduced force at the other end, or the velocity of the lever’s end point is converted to another velocity value related to another point proportional to their distances from the lever’s pivot.

As shown in Figure 3‑15, a $TF$ -element can have one effort and one flow as inputs, and consequently, delivers converted corresponding effort and flow as outputs. The conversion parameter $\lambda$ should be defined, based on the physical system data. For example, for the case that flow $f_1$ is the input and flow $f_2$ the output, we can write $f_2 = \lambda f_1$ to define $\lambda$ . But from energy conservation we have $e_1 f_1 = e_2 f_2$ , or the output effort $e_1 = \underbrace{(f_2/f_1)}_{\lambda}e_2 = \lambda e_2$ . Similarly, for the case that effort $e_1$ is the input and effort $e_2$ the output, we can write $e_2 = \lambda^{-1} e_1$ , using $\lambda$ . But from energy conservation, we have, $e_1 f_1 = e_2 f_2$ or the output flow $f_1 = \underbrace{(e_2/e_1)}_{{1}/\lambda}f_2 = \lambda^{-1} f_2$ .

These relations constitute the $TF$ -element equations and are shown in Figure 3‑15, for each case where the inputs to the $TF$ -element are identified with thick arrows.

Figure 3-15 Block diagrams (left) and equivalent bond graphs for a $TF$ -element with related assigned causalities—inputs are shown with thick arrows

Note that the $TF$ -element should have only one of the two required causality strokes near it for either cases, as shown in Figure 3‑15.

In a bond graph model for the $TF$ -element, there should be only one causality stroke close to it and another one away from it. A $TF$ -element converts flows to flows and efforts to efforts.

A gyrator element, $GY$ , represents the converter that receives a type of physical quantity and delivers a different type after conversion. Examples are a DC motor which converts voltage (effort) to angular velocity (flow) of the rotor or the attached shaft. The reverse operation is that of an electric generator.

As Figure 3‑16 shows, a $GY$ -element can have one effort and one flow as inputs and, consequently, delivers corresponding flow and effort as outputs. The conversion parameter $\kappa$ should be defined, based on the physical system data. For example, for the case with flow $f_1$ as the input and effort $e_2$ being the output, we can write $e_2 =\kappa f_1$ to define $\kappa$ . But from energy conservation we have $e_1 f_1 = e_2 f_2$ , or the output effort $e_1 = \underbrace{(e_2/f_1)}_{\kappa}f_2 = \kappa f_2$ . Similarly, for the case with effort $e_1$ as the input and flow $f_2$ being the output, we can write $f_2 = \kappa^{-1}e_1$ , using $\kappa$ . But from energy conservation we have $e_1 f_1 = e_2 f_2$ , or the output flow $f_1 = \underbrace{(f_2/e_1)}_{1/ \kappa} e_2 = \kappa^{-1}e_2$ . These relations constitute the $GY$ -element equations and are shown for each case where the inputs for the $GY$ -element are identified with thick arrows in Figure 3‑16.

Figure 3-16 Block diagrams (left) and equivalent bond graphs for a GY-element with related assigned causalities—inputs are shown with thick arrows

Note that the $GY$ -element should have both required causality strokes near it or away from it, as shown in Figure 3‑16.

In a bond graph model for the $GY$ -element, there should be two causality strokes close to it or both away from it. A $GY$ -element converts flows to efforts and efforts to flows.

Now that we have all nine elements of BG method defined, in the following sections we discuss the state variables, their definitions, and relation with integral causality. State variables are key quantities in analyzing engineering system dynamics and behaviour and are a critical part of BG method. A sound understanding of the state variables will help in developing a high level of competency in BG method and its applications to engineering systems.

3.5 System State Variables

The main objective of BG models is to derive system equations that describe the behaviour of the system and to follow up by solving these equations for simulation and design purposes.

The system equations may be ODEs of second order or higher. However, when writing these governing system equations in terms of state variables—those variables that uniquely and sufficiently describe the system dynamics—we end up having first-order ODEs, a huge advantage when using numerical/analytical solution methods. In addition, when we extract system equations from the corresponding BG model (see chapter 11), additional algebraic equations are involved; hence, we have a system of differential-algebraic equations (DAEs) that could benefit from having the related ODEs written as first-order equations.

In this section, we define the state variables that relate themselves to the storage elements in BG method i.e., $I$ -element and $C$ -element. Other BG elements correspond to the algebraic equations of the system DAEs and do not possess state variables of their own.

We now consider the kinetic energy storage element or inertia $I$ -element. The energy stored can be written as the integral of power (i.e., effort multiplied by flow) with respect to time, $\int e \cdot f\: dt$ or as $\int e \cdot f\: dt=\int(edt)f$ . But $edt=dp$ , the generalized momentum differential/change. Hence, $\int (edt) f = \int f\: dp$ , or the energy stored in an inertia element $I$ is the integral of flow (e.g., velocity) with respect to momentum as the independent variable. Therefore, a functional form of the type $f = \mathcal{F} (p)$ is required to perform the integral operation. In other words, the area under the curve of the flow $f$ in the $f-p$ coordinate system is equal to the energy stored. Recall that, e.g., in mechanical systems, this function (i.e., $f = \mathcal{F} (p)$ ), is derived from Newton’s second law, or $f = p/I$ (the parameter $I$ is mass or inductance, for example). Therefore, we have $\textit{Energy} =\int f\: dp=\dfrac{1}{I}\int p\: dp=\dfrac{p^2}{2I}$ , or

(3.7) $\begin{equation*} \textit{Energy stored as an I element} = \frac{p^2}{2I} \end{equation*}$

Equation (3.7) clearly shows that the energy stored by an $I$ -element is uniquely defined by its generalized momentum. Therefore, the momentum of an $I$ -element is identified as a state variable of the system.

The generalized momentum associated with Inertia element in the bond graph model is a system state variable, so-called $p$ on $I$ .

Similarly, we consider the potential energy storage element, or $C$ -element. The energy stored can be written as the integral of power with respect to time, $\int e \cdot f\: dt$ or as $\int e \cdot f\: dt=\int(fdt)e$ . But $fdt=dq$ , the generalized displacement differential/change. Hence, $\int(fdt)e=\int e\: dq$ or the energy stored in a $C$ -element is the integral of effort (e.g., force) with respect to displacement as the independent variable. Therefore, a functional form of the type $e= \mathcal{H} (q)$ is required to perform the integral operation. In other words, the area under the curve of $e$ as a function of $q$ in the $e-q$ coordinate system is equal to the energy stored. Recall that, e.g., in mechanical systems, this function (i.e., $e=\mathcal{H} (q)$ ) is derived from Hooke’s law, or $e=q/c$ (the parameter $c$ is spring compliance or capacitor capacitance, for example). Therefore, we have $\textit{Energy} = \int e\:dq=\dfrac{1}{c}\int q\:dq=\dfrac{q^2}{2c}$ , or

(3.8) $\begin{equation*} \textit{Energy stored in a C-element} = \frac{q^2}{2c} \end{equation*}$

Equation (3.8) clearly shows that the energy related to a $C$ -element is uniquely defined by its generalized displacement. Therefore, the displacement of a $C$ -element is identified as a state variable of the system.

The generalized displacement associated with a $C$ -element in the bond graph model is a system state variable, so-called $q$ on $C$ .

These two state variables ( $p$ and $q$ ) are key variables when extracting system equations from the corresponding bond graph (see chapter 11). The total number of independent system equations is equal to the total number of state variables, or $p$ on $I$ and $q$ on $C$ .

The reader should also note that the assumed governing equations for these two elements (i.e., Newton’s second law for $I$ -elements and Hooke’s law for $C$ -elements) determine the functional forms of $f= \mathcal{F}(p)$ for an $I$ -element and $e=\mathcal{H} (q)$ for a $C$ -element, respectively. Other constitutive equations: e.g., non-linear relations could be used if desirable, but the uniqueness of energy stored on the $p$ and $q$ remains for each of these two elements.

3.5.1 Integral Causality and State Variables: I– and C-elements

The main objective of assigning a causality stroke to an element is to make the element definite in terms of its inputs and outputs (i.e., either effort or flow). Since we have two choices (either effort or flow being the input or the output), the preferred causality is the one that, when assigned, allows the input to the element such that the element-related laws of physics are satisfied and the state variable is concluded as well. For example, if an element receives effort, then it should respond with flow, and the related state variable should be the outcome of the application of the laws of physics to this element. These objectives are met when we use the integral causality strokes for $I$ -element and $C$ -element. In other words, when the integral of the cause signal is equal to the state variable of the corresponding storage element, then that element is integrally causalled.

In the previous sections (see sections 3.4.1 and 3.4.2), we discussed the preferred causalities for $I$ – and $C$ – elements as being the integral causality types. Having defined the state variables for $I$ – and $C$ – elements (see section 3.5), we can expand the discussion on why the integral causality is the preferred one for these elements.

The state variable corresponding to elements $I$ and $C$ specifies the assignment of causality strokes for the integral causality is defined such that the integral of input quantity (either effort or flow) for $I$ – or $C$ -elements result in the corresponding state variable.

Recall that generalized momentum $p$ is the state variable for an $I$ -element. Now, we consider the choice of having the flow $f$ or effort $e$ as the input for $I$ -element according to the causality stroke assignment (see Figure 3‑4). When the effort is selected as the input, we can integrate it (hence, the designation of integral causality for this choice), and get the momentum, i.e., the state variable, as well as the flow for the element response. This is consistent with the $I$ -element governing equation (i.e., Newton’s second law). Therefore, having the causality stroke at the port of $I$ -element, or the preferred causality assignment (see Figure 3‑4), satisfies all the mathematical requirements and provides the flow as the response and the momentum as the state variable. The whole process is shown in Figure 3‑17. The choice of having flow as the input for $I$ -element—the derivative causality—does not fulfill all the objectives mentioned above; hence, it is not preferred. Note that when derivative causality is assigned, Newton’s second law still is satisfied, but the state variable is not explicitly involved.

Figure 3-17 Consistency of integral causality assignment and state variable for an $I$ -element with parameter $m$

Similarly, for a $C$ -element, we can have a similar argument. Recall that generalized displacement $q$ is the state variable for a $C$ -element. Now, we consider the choice of having the flow $f$ or effort $e$ as the input for $C$ -element according to the causality stroke assignment (see Figure 3‑7). When the flow is selected as the input, we can integrate it (hence the designation of integral causality for this choice) and get the displacement, i.e., the state variable, as well as the effort as the element’s response. This is consistent with the $C$ -element governing equation, i.e., Hooke’s law. Therefore, having the causality stroke away from the port of $C$ -element, or the preferred causality assignment (see Figure 3‑7) satisfies all the mathematical requirements and provides the effort as the response and the displacement as the state variable. The whole process is shown in Figure 3‑18. The choice of having effort as the input for $C$ -element—the derivative causality—does not fulfill all the objectives mentioned above; hence, it is not preferred. Note that when derivative causality is assigned Hooke’s law still is satisfied but the state variable is not explicitly involved.

Figure 3-18 Consistency of integral causality assignment and state variable for a $C$ -element with parameter $c$

Exercise Problems for Chapter 3

Exercises

Using Figure 3‑1, identify each component in terms of its type related to energy storage, dissipation, converter, and source.
Using Figure 3‑3, explain if the power bond direction and causality stroke assignment are independent from each other or dependent.
List nine basic bond graph elements and sketch them with their preferred causalities, where applicable.
For each bond graph sketch, perform the operations given below:
1. Write the energy rate balance equation at each junction
2. Identify strong power bond.
3. Assign all remaining causality strokes, using red colour to distinguish them
The massless lever rotates about the pivot point with angular velocity , as shown in the below sketch. Draw the bond graph model using TF- element along with equation model for each case:
1. Velocity magnitude at point A is given, $V_A$ . Calculate the transformer parameter $\lambda_V$ .
2. Force magnitude at point A is given, $F_A$ . Calculate the transformer parameter $\lambda_F$
3. Discuss the relation between $\lambda_V$ and $\lambda_F$ .
Describe system state variables and explain their significance related to a system’s equations. Identify BG elements associated with these variables.
Discuss the principle of cause and effect in relation to causality assignment in BG method. For the following elements, assign the causalities and identify the cause and effect for each one. Also identify the integral vs. the derivative causality.

$\longrightharpoondown I \longrightharpoondown C \longrightharpoondown R$

Media Attributions

Portrait of a Mathemetician © Mary Beale is licensed under a Public Domain license
Figure 3-16

force per unit displacement of the spring

displacement per unit force applied to the spring

Building Bond Graph Models: General Procedure and Application

4.1 Overview

To demonstrate applications of BG method, we discuss the procedure for building BG models for physical systems, using the material presented in chapter 3. We use examples related to mechanical systems to establish the guidelines and steps required to build a BG model. In further chapters, we present more worked-out examples for several engineering systems and disciplines, including electrical, hydraulic systems.

4.2 Steps for Building Bond Graph Models: General Guidelines

As mentioned, BG method can be used to build models for single- and multi-domain physical systems. The building blocks are the nine basic BG elements, including their modulated versions, and causality assignment rules (see section 3.4). A model for any specific system also requires definitions of relevant sign conventions for general displacement and forces. This chapter will discuss the latter and will present some worked-out examples. The following are the steps for building a BG model, in general, including for mechanical, electrical, and hydraulic systems:

Identify the physical system components in terms of their type (energy storage, source, dissipater, etc.).
Identify the DOF (degrees of freedom) of the system. This step is optional but recommended.
Identify and list the required BG elements.
Identify distinct physical points/nodes of the physical systems:
- velocity or force (mechanical systems): translational
- angular velocity or torque (mechanical systems): rotational
- voltage or current (electrical systems): electrical circuits
- pressure or flow rate (hydraulic systems): fluid network
Assign proper BG multi-port junction elements to items from step 4:
- “1” for velocity, angular velocity, and electrical and flow currents
- “0” for force, voltage, and pressure
Connect associated elements, using BG elements and power bonds, to the items from step 5.
Assign proper BG multi-port junction elements in between those items from step 5:
- “0” for relative velocity and angular velocity
- “1” for voltage drop and pressure drop
- $TF$ and $GY$ for energy conversion
Connect associated elements to items from step 7, using BG elements and power bonds.
Define sign convention and connect all remaining power bonds.
Apply all causality assignments (integral causalities must be given priority).
Draw and build the BG model in 20-sim (when available).
Perform simulation and design, using the obtained BG model (when required).

In further sections, we will demonstrate implementation of the procedure/algorithm mentioned above with some worked-out examples, including power bond direction, causality assignment, and sign convention.

4.2.1 Guidelines for Power Bond Direction

Connecting elements in a BG model with power bonds requires compliance with the direction of energy flow in the physical system. Therefore, the directions of half-arrows are critical. The following guidelines may be helpful:

Draw power bonds from BG source elements ( $S_e$ and $S_f$ ) toward the system, connecting to the adjacent elements.
Draw power bonds toward BG passive elements (i.e., $I$ , $C$ , and $R$ )
Draw power bonds to and from BG junction elements (“1” and “0”) according to a previously defined sign convention (see section 3.4.5).
Draw remaining power bonds to have all BG elements connected.
Some simplifications of the BG model may be justified, but not required.

After drawing all power bonds for the model, assign the causality strokes. The next section provides a list of guidelines for causality assignments.

4.2.2 Guidelines for Assigning Causality Strokes

The assignment of causality strokes is a required step in building any BG model. The following steps help with achieving this requirement.

Assign causality to BG source elements.
Assign causality assignments with preferred integral causality strokes to $I$ – and $C$ – elements.
As far as possible, extend the causality assignments to other power bonds, using the causality requirements for connecting elements (e.g., 1, 0, $TF$ and $GY$ )
Assign causality assignments to $R$ -elements that accept neutral causality stroke assignment.
As far as possible, using the causality requirements for connecting elements, extend the causality assignments to all remaining power bonds in the model.

If execution of step 5 from the above list cannot be completed, then the BG model contains some specific mathematical properties—algebraic loop or differential/derivative causality (see chapter 11). The application 20-sim automatically assigns the causality strokes with prioritizing integral causalities, and if present, identifies the derivative or algebraic loops causalities in the model with red-colour strokes. In further sections, we will explore these features, with some examples.

4.3 Example: BG Model for a One-DOF Mass-Spring-Damper Mechanical System

A mechanical system consists of mass $m$ [kg], spring $k$ [N/m], and damper $b$ [N.s/m]. The applied force on mass is $F(t)$ . Build a BG model for this system as shown in Figure 4‑1, neglecting friction of the rollers.

Figure 4-1 A mass-spring-damper mechanical system

Solution:

DOF = 1 (1D translational motion of one mass) and the required BG elements are: $I$ (representing the mass), $C$ (representing the spring), $R$ (representing the damper), $S_e$ (representing force $F$ ) and $S_f$ (representing the wall velocity). Also, we are required to have junctions “1” and “0.”

Distinct velocity nodes are the mass and the wall (although the wall usually is stationary). Hence, we need two “1” junctions to represent common velocity for all elements attached to the mass and the wall.

We draw them as

As well, for each junction, it is useful to assign a name related to its representation.

We draw all elements connecting to the junctions which have the same distinct velocity values and connect them with power bonds. Therefore, for the wall-velocity junction, we use flow source $S_f$ and for mass-velocity junction and inertial element $I$ , representing the mass $m$ and a $S_e$ representing the applied force $F(t).$ Note that $I$ -element should receive power (passive element), and sources send power to the system (active elements).

The spring and damper experience the same value of relative velocity, | $V_m-V_W$ |, which is represented by 0-junctions. Recall that 0-junction element is a flow summator. According to the power bonds connecting the spring (or damper) to the 0-junction, we can have ( $V_m>V_W$ ) or ( $V_m<V_W$ ), considering the $x$ -coordinate as given in Figure 4‑1. Therefore, to specify the associated power bond directions we should define a sign convention. The common practice is to consider the spring (or damper) from the BG model and define either tension force as being positive (+T) or the compression force being positive (+C). For this example, we use the spring displacement/velocity to demonstrate the sign convention. A similar argument applies for the damper’s displacement.

To represent the relative velocity, we add two 0-junctions and $C$ – and $R$ – elements to the model and use, e.g., (+T) sign convention, as shown below:

Note that $C$ – and $R$ – elements are passive and should receive power from the system. After labelling the bonds connecting to the 0-junction associated with the $C$ – element, we can write the power balance as $e_2f_2-e_3f_3-e_1f_1=0$ . But $e_1=e_2=e_3$ . Hence, $f_2-f_3=f_1$ , or $V_m-V_W=V_C$ where $V_C$ is the spring displacement rate or velocity equal to the relative velocity. Now, to have the displacement of the spring in the + $x$ direction, we should have $V_C \ \textgreater \ 0$ or $V_m>V_W$ . This implies that the displacement/velocity of the mass should be larger than that of the wall, for the spring is experiencing a positive tension force. Therefore, the spring is under tension and the assigned sign convention (+T) is satisfied, considering the + $x$ direction.

Now, if we change the power direction of the bonds connecting to the 0-junction, as shown in the sketch below, we have $V_W-V_m=V_C>0$ , or $V_W>V_m$ ; hence, we have the spring under compression. Therefore, the (+C) sign convention is satisfied. Note that only the power direction of the four bonds associated with the two 0-juctions can change their directions since the rest are associated with source or passive elements and are unique in their directions, as shown below.

Both sign conventions are legitimate, but only one should be selected and used consistently for building a BG model. We continue, using the BG model with (+T) sign convention.

Causality assignments are now applied, according to the rules discussed in chapter 3. Following the guidelines given in section 4.2.2, we start applying the causality to the source elements, followed by those for $I$ – and $C$ – elements. Recall that integral causalities are preferred for elements $I$ (i.e., $I$ receives effort) and $C$ (i.e., C sends effort). The causality strokes are shown with transvers lines, as shown below.

Extend the causality assignments to the remaining bonds, using the rules for 1-junction (can receive only one flow signal through its strong bond) and 0-junction (can receive only one effort signal through its strong bond), as shown below in the model sketch (see Figure 4‑2).

Figure 4‑2 Bond graph model for a one-DOF mass-spring-damper system

To a reasonable extent, we can simplify the BG model such that it clearly resembles the physical system. The two 0-junctions represent the same value of relative velocity, $V_m-V_W$ . Therefore, we can combine them into a single 0-junction and share the relative velocity value through a 1-junction element with the $C$ – and $R$ – elements. This simplification becomes very useful for building large BG models for more complex systems. Figure 4‑3 shows the resulting BG model. Note that the causality strokes should be adjusted after simplifications are made.

Figure 4‑3 Simplified bond graph model for a one-DOF mass-spring-damper system

4.4 Example: BG Model for a Two-DOF Mass-Spring-Damper Mechanical System

Build the BG model for the mechanical system as shown in Figure 4-4. Consider the (+C) to be the sign convention for internal forces.

Figure 4-4 A two-DOF mass-spring-damper mechanical system

Solution:

This system has two DOF and four distinct velocity points, corresponding to mass $m_1$ and $m_2$ and the two walls. Therefore, we lay out four 1-junctions to represent them in the model. The remaining required BG elements are $I$ , $C$ , $R$ , $S_e$ , $S_f$ , and 1- and 0- junctions.

We follow the same guidelines demonstrated in the previous example (see section 4.3) and build the BG model as shown in Figure 4-5.

Figure 4-5 BG model for a two-DOF mass-spring-damper mechanical system

The reader is encouraged to build this BG model and to compare the results with those provided in Figure 4-5. The 20-sim BG model and a screen recording are available as companion resources describing the process to build the equation model and typical results for this example.

One or more interactive elements has been excluded from this version of the text. You can view them online here: https://pressbooks.bccampus.ca/engineeringsystems/?p=62#oembed-3

4.5 Example: BG Model for a Three-DOF Mass-Spring-Damper Mechanical System

Build the BG model for the mechanical system as shown in Figure 4-6. Consider the (+C) to be the sign convention for internal forces.

Figure 4-6 A three-DOF mass-spring-damper mechanical system

Solution:

This system has three DOF and five distinct velocity points corresponding to mass $m_1$ , $m_2$ , and $m_3$ and the two walls. Therefore, we lay out five 1-junctions to represent them in the model. The remaining required BG elements are $I$ , $C$ , $R$ , $S_e$ , $S_f$ , and 1- and 0-junctions.

We follow the same guidelines demonstrated in the previous example (see section 4.3) and build the BG model, as shown in Figure 4-7.

Figure 4-7 BG model for a three-DOF mass-spring-damper mechanical system

The reader is encouraged to build this BG model and compare the results with those provided in Figure 4-7. The 20-sim BG model and a screen recording are available as companion resources describing the process to build the equation model and typical results for this example.

One or more interactive elements has been excluded from this version of the text. You can view them online here: https://pressbooks.bccampus.ca/engineeringsystems/?p=62#oembed-4

4.6 Example: Kinetics and Kinematics of a Mechanical System Using BG Model

As mentioned previously (see section 3.2), one of the advantages of BG modelling method is that a BG model allows us to gain insights by visual inspection of the BG model. This can be achieved by drawing the streams of efforts (kinetics) and flows (kinematics) for a BG model.

In this example, we use the results from the example given in section 4.3 to demonstrate this property by explicitly drawing the effort and flow associated with each power bond in the model. First, we look at the kinetics of the system by drawing the efforts, as shown in Figure 4-8. As shown, the efforts/forces associated with the spring $e_C$ and dumper $e_R$ are collected as force $e$ and transferred to the mass $m$ in addition to the applied force $F$ shown as $e_F$ . Clearly, the wall receives the collected force $e=e_C + e_R$ .

Figure 4-8 Kinetics of a one-DOF mechanical system showing the stream of efforts with its BG model

Similarly, by drawing the flows, as shown in Figure 4-9 the kinematics of the system can be visualized. As shown, the flows/velocities associated with the mass $f_m$ and wall $f_w$ are collected as velocity $f$ and transferred to the spring and damper. Clearly, these elements receive the relative velocity $f=f_m+ f_w$ due to the motion of mass and the wall (if stationary, $f_w=0$ ).

Figure 4-9 Kinematics of a one-DOF mechanical system showing the stream of flows with its BG model

4.7 Modelling and Simulation Approaches in Engineering: Modern vs. Traditional

Considering BG—our focus in this textbook—as the modelling method, once we have the corresponding BG model, we can proceed to simulation, and hence, design of a system. One can take two approaches to perform this task: traditional or modern. As mentioned, the main objective of modelling and simulation is to help with more effective design of the systems in terms of their cost, function, material consumption, etc. Therefore, any modelling method, including bond graph, should result in a mathematical model consisting of the systems’ equations. The solution of the equations can be used for system simulation and analysis to support effective system design. In the following sections, we briefly describe possible approaches and the reason for choosing the modern approach over the more traditional one.

4.7.1 Traditional Approach

Once the BG model is available for a system, we can derive/extract the system equations from the BG model—usually a laborious task—and then using numerical methods, we can obtain their solutions. This task is usually achieved with the help of computer programs (usually developed from scratch) based on a selected numerical method. This approach, along with both its system equation extraction and especially the computer coding, is limited in practical application, being specific from one problem to another one, and is laboriously time consuming. Figure 4-10 shows the major steps of the traditional approach.

Figure 4-10 Traditional approach for system simulation and design

In practice, it is inefficient to develop computer codes for each specific design: the amount of person power, computer power, and other resources become overwhelming for the fast-paced engineering design needs of today’s industries. Therefore, a huge effort has been made to develop commercially available software tools to help meet these objectives and to make the whole process of system design more effective, economically viable, and efficient.

4.7.2 Modern Approach

Alternatively—or rather, preferably—the modern approach in engineering and system design employs related software tools. These tools provide opportunities to perform systems simulation immediately after obtaining the BG model. The software tool that we introduce and use in this textbook, 20-sim, helps with extracting the system equations from the BG model seamlessly and provides facilities for system simulation and parametric analysis. This modern approach is more effective in engineering practice and provides more and quicker insights into engineering systems design. In addition, the modern approach helps to respond more effectively to the fast-paced engineering demands in industry and is recommended for engineers in practice. Figure 4-11 shows the major steps of the modern approach. Note that verification and validation should be considered in the modelling step as well.

Figure 4‑11 Modern approach for system simulation and design

In the next section, we introduce the 20-sim software package with a focus on bond graph modelling, simulation, and time and frequency analysis for engineering systems and design. In further sections, we use 20-sim to build BG models and their simulations and to study their dynamical behaviour.

Exercise Problems for Chapter 4

Exercises

Build the BG model, including causality assignment, for the example given in section 4.4 considering (+T) as the sign convention for internal forces.
1. Draw a kinetic map of the system, using the stream of efforts.
2. Draw a kinematic map of the system, using the stream of flows.
Build the BG model, including causality assignment, for the example given in section 4.5 considering (+T) as the sign convention for internal forces.
1. Draw a kinetic map of the system, using the stream of efforts.
2. Draw a kinematic map of the system, using the stream of flows.
Discuss the benefits of modern vs. traditional approaches for simulation and design of systems, from a practical point of view. Include speed of calculations, and economical aspects of the two methods in your discussion.

Media Attributions

fig4-3

Recall that 1- junction is a flow equalizer (or effort summator) and 0-junction is an effort equalizer (or flow summator).

Introduction to 20-sim Software Tool

5.1 Overview

In this chapter, we introduce the engineering software tool 20-sim, focusing on its BG modelling and simulation facilities. Originally developed in the Netherlands at the University of Twente and released in 1995, 20-sim was the first commercially released software package for bond graph modelling. 20-sim was further developed and currently is available through Controllab Products [3].

An advanced modelling and simulation software tool, 20-sim uses Microsoft Windows™ operating system. 20-sim enables users to model the behaviour of systems, such as mechanical, electrical, hydraulic, thermal, or a combination of these systems (i.e., multi-domain systems). Below is a list of 20-sim’s modelling and simulation tools. The first four items in the list are the main tools:

bond graph
block diagrams
iconic diagrams
equation models
3D animation
3D mechanics
code generation
controller design
frequency and time domain
multi-domain systems
virtual reality: Unity Toolbox

20-sim allows users create models graphically, similar to drawing an engineering sketch. With these models you can simulate and analyze the behaviour of single-domain and multi-domain dynamic systems and create control systems. One can even generate C-code and run this code on hardware for rapid prototyping and hardware-in-the-loop (HIL) simulation.

20-sim has two main window interfaces: Editor and Simulator, shown in Figure 5‑1 and Figure 5‑2, respectively. To open the Simulator window from the keyboard, select Ctrl + R. To open the Simulator window from the Editor window, choose Tools > Simulator.

Figure 5-1 20-sim Editor interface

Figure 5-2 20-sim Simulator interface

As described in the manual Getting started with 20-sim 4.8, the Editor consists of four parts:

Model tab/Library tab: This is the part at the middle left. The Model tab shows the model hierarchy, i.e., the hierarchical composition (all the elements) of the model that is created in the Editor. The Library tab shows the 20-sim library.
Graphical Editor/Equation Editor: This is the big white space, with or without grids, at the middle right. The Editor is used to create graphical models and enter equations.
Output tab/Process tab/Find tab: This is the part at the bottom right. The Output tab shows the files that are opened and stored. The Process tab shows the compiler messages. The Find tab shows the search results.
Interface tab/Icon tab: This is the part at the bottom left. The Interface tab shows the interface of a selected model. Double-click to open the Interface Editor.
The Icon tab shows the icon of a selected model. Double-click to open the Icon Editor.

Here is a video guiding the reader through the features of 20-sim:

One or more interactive elements has been excluded from this version of the text. You can view them online here: https://pressbooks.bccampus.ca/engineeringsystems/?p=64#oembed-5

More information about using 20-sim is available at 20-sim webhelp and at 20-sim Help Manuals.

Figure 5‑3 shows the main steps for system design. In further sections, we will follow these steps with several worked-out examples to demonstrate the application of this process.

Figure 5-3 Process steps for design of a system using modelling, simulation, and analysis

5.2 Example: BG Model for a Car Seat Mechanical System

To demonstrate 20-sim application, using the data provided in Table 5‑1, we build a BG model of a car seat. The car seat system schematic is shown in Figure 5‑4. Consider the (+C) to be the sign convention for internal forces.

Figure 5-4 Sketch for a mechanical system for a car seat

Table 5‑1 Data for car-seat example and BG model
Component/element	Value	Unit
Mass, $m_1$	27	kg
Mass, $m_2$	80	kg
Spring, $k_1$	1800	N/m
Spring, $k_2$	19 x 10⁴	N/m
Damper, $b_1$	400	N.s/m
Damper, $b_2$	900	N.s/m
Car floor speed, $V_floor$	signal	step-function

Solution:

Here are two videos showing how to build and run the model for this example in 20-sim:

One or more interactive elements has been excluded from this version of the text. You can view them online here: https://pressbooks.bccampus.ca/engineeringsystems/?p=64#oembed-6

One or more interactive elements has been excluded from this version of the text. You can view them online here: https://pressbooks.bccampus.ca/engineeringsystems/?p=64#oembed-7

Launch 20-sim and follow the solution steps provided in the videos. Figure 5‑5 shows the resulting BG model.

Figure 5-5 Bond graph model for a car seat

5.3 Example: BG Model for a Cart Carrying a Load

To demonstrate 20-sim, we build a BG model of a cart carrying a load. Consider the (+C) to be the sign convention for internal forces.

Table 5‑2 shows the data for the system components. Figure 5‑6 shows the cart system schematic.

Table 5‑2 Data for example of cart with load and BG model
Component/element	Value	Unit
Mass, $m_1$	40	kg
Mass, $m_2$	20	kg
Spring, $k_1$	2000	N/m
Spring, $k_2$	1500	N/m
Damper, $b_1$	400	N.s/m
Damper, $b_2$	900	N.s/m
Force, $F(t)$	signal	sinusoidal

Figure 5‑6 Sketch for a cart system carrying a load

Solution:

The process for building a BG model and related simulation for this example is given through a screen recording that guides the readers through the whole process. Launch 20-sim, and to build and run the 20-sim model for this example, follow the solution steps provided in the following video:

One or more interactive elements has been excluded from this version of the text. You can view them online here: https://pressbooks.bccampus.ca/engineeringsystems/?p=64#oembed-8

Figure 5‑7 shows the resulting BG model.

Figure 5‑7 Bond graph model for a cart carrying a load

Exercise Problems for Chapter 5

Exercises

Familiarize yourself with the 20-sim tools and features, using the screen recording in section 5.1.
Build the BG model using 20-sim for the example given in section 5.2 considering:.
1. (+T) as the sign convention for internal forces
2. (+C) as the sign convention for internal forces

Run the model for simulation and create and report typical graphs.

Build the BG model using 20-sim for the example given in section 5.3 considering:
1. (+T) as the sign convention for internal forces
2. (+C) as the sign convention for internal forces

Run the model for simulation and create and report typical graphs.

Bond Graph Models for Complex Mechanical Systems

6.1 Overview

In the previous chapters, we established concepts such as the basic elements of bond graph method and the algorithm for building BG models. We now continue with more worked-out examples for selected complex mechanical systems. These systems may have many components, involve many degrees of freedom, and exhibit translational and rotational motions in one-dimensional (1D) or two-dimensional (2D) space. So far, we have used translational mechanical systems and demonstrated how to build their related BG models (see chapters 4 and 5). In this chapter, we expand the discussion to rotational mechanical systems with rotational and/or 2D/plane rigid-body motions, including their related BG model examples. First, we establish the theories and related equations and then use those for building the BG models.

6.2 Mechanical Systems—Rotational

A mechanical system may consist of rotational components, e.g., shafts, discs, gears, pulleys, and levers. The generalized BG elements and relations apply to the motion of rotational components in a similar way that the translational motion was treated; i.e., they are analogous (see Table 3‑1). In other words, rotation angle $\theta$ is equivalent to the generalized displacement $q$ , angular velocity $\omega=\dot\theta$ to the flow $f$ , and torque $T$ to the effort $e$ . The polar moment of inertia $J$ is represented by $I$ -element, the shaft by $C$ -element, and bearing by $R$ -element. The generalized momentum is the integral of $T$ with respect to time. Therefore, we can write

$\begin{equation*} \begin{dcases} \textit{e}\equiv\textit{T},\textit{torque}\\ \textit{f}\equiv\omega=\dot\theta,\textit{angular velocity}\\ \textit{e}\cdot\textit{f}\equiv T\cdot\omega,\textit{power}\\ \textit{q}\equiv\theta=\int{\omega}\:{dt},\textit{twist angle}\\ \textit{p}\equiv{p_T}=\int\textit{T}\:{dt},\textit{angular momentum}\\ \end{dcases} \end{equation*}$

where $p_T$ is the rotational momentum or so-called angular momentum. Using the constitutive relations, for an $I$ -element we have $e=I\cdot\dot{f}$ or $T=J\dot\omega$ , for a $C$ -element we have $e = q / c$ or, $T=\theta/c$ where $c$ represents the torsional compliance or inverse of torsional stiffness $k_T$ , $c = 1 / k_T$ . Similarly, for an $R$ -element we have $e=Rf$ or $T={b_T}\omega$ where $b_T$ is the friction of the torsional bearing. The energy associated with storage elements can be written using Equations (3.7) and (3.8), or for elements $I$ and $C$ , as $\dfrac{1}{2J}{p}{_T^2}$ and $\dfrac{1}{2k_T} {T^2}$ , respectively. One advantage of bond graph method is its analogous applicability to different domains using the common constitutive relations, as described above for rotational motion.

Table 6‑1 shows typical rotational components.

Table 6-1 Typical rotational mechanical components and their corresponding BG elements
$R$ -element (bearings)	$C$ -element (shaft)	$I$ -element (disc)	$TF$ -element (gears)	$GY$ -element (gyrator)

6.3 Mechanical Systems—Two-Dimensional Rigid Plane Motion

The components of mechanical systems that we considered so far are assumed as point masses. In other words, they are point elements but can have motions either in translation and/or rotation. However, two-dimensional components such as rigid plates, car chassis, and thin rods can have relative 2D motion and cannot be treated as point elements.

In general, a 3D solid component/body has six degrees of freedom; i.e., its centre of mass can move in three translational directions and through three associated rotational angles. In many mechanical systems, however, we can assume components as two-dimensional planes with negligible deformations, or as 2D rigid bodies having three degrees of freedom: two in-plane translations and one rotation about the perpendicular axis to the plane of motion.

Building a BG model requires transforming the velocities and angular velocities associated with the rigid plane and making them available to the contact points with other components of the system. For example, a car’s chassis moving forward on a wavy road may experience rotations like pitch and roll (i.e., rotations about the axes parallel to the ground) in addition to translational motion. Considering the chassis as a 2D rigid body, we need to know how the linear and angular velocities are transmitted to the suspensions connecting to it.

We present here an analysis of 2D rigid-body motion, with focus on applications to BG modelling. For further readings on this topic, consult with available references [13], [18], [20], [23].

The general motion of a 2D rigid body can be decomposed into translation of the whole body and a rotation about a fixed point of the body. This is the result of the principle of superposition and can be shown using the geometry of the motion. As Figure 6‑1 shows, we assume a rigid body going through a planar motion with reference to a fixed coordinate system $x-y$ . We identify a line/vector on the body connecting two arbitrarily selected points A and B, with point B taken as a reference, usually the centre of mass. We then capture a picture of the body at a later time, $\Delta{t}$ during its motion, and find out the line BA in its new orientation and position, as shown in the sketch on the right in Figure 6‑1. Since the body does not deform, the length of the line BA (or the magnitude of vector $\overrightarrow{BA}, BA=|\overrightarrow{BA}|=|\overrightarrow{r}|=r$ ) remains constant. Using this property, we can draw a circle with its centre at the new position of point B and radius of BA. Then, we draw the radial line BA’ parallel to the line BA at its initial position. To orient BA’ according to the new position of BA, we then rotate BA’ about point B through angle $\theta=\omega\triangle{t}$ , where $\omega=|\overrightarrow{\omega}|$ is the magnitude of the angular velocity vector of the rigid body, perpendicular to the plane of motion. Consequently, we can claim that original point A is translated (not rotated) by the velocity of point B, $\overrightarrow{V}_B$ from its initial position to a new position A’ and subsequently rotated about point B by angle $\theta$ to orient in its final new position of point A. The initial position of point A is arbitrarily selected; therefore, the argument equally applies to all points of the rigid body.

Figure 6-1 Decomposition of 2D rigid-body motion into translation and rotation

The general planar motion of a 2D rigid body can be decomposed into the translation of an arbitrary point on the body followed by a rotation about the point.

Mathematically, we can write $\overrightarrow{V}_A=\overrightarrow{V}_B+\overrightarrow{V}_A/_B$ . The relative velocity $\overrightarrow{V}_A/_B$ is the tangential velocity due to rotation and can be written as $\overrightarrow{\omega}\times\overrightarrow{BA}$ . Therefore, $\overrightarrow{V}_A/_B = \overrightarrow{\omega}\times\overrightarrow{r}=(0,0,\omega)\times (r_x,r_y,0)=(-{\omega}{r_y},{\omega}{r_x},0)$ . Therefore, we can write the velocity components of point A resulted from rigid-body motion as

(6.1) $\begin{equation*} \left\{\begin{array}{rcl} {V}_A_x=V_B_x-\omega{r_y}\\ V_A_y=V_B_y+\omega{r_x}\\ \end{array}\right \end{equation*}$

But $r_x=r\cos\theta$ and $r_y=r\sin\theta$ where, $\theta$ is the angle between vector $\overrightarrow{r}$ and positive direction of $x$ -axis. After substituting, we get

(6.2) $\begin{equation*} \left\{\begin{array}{rcl} V_A_x=V_B_x-\omega{r}\sin}\theta}\\ V_A_y=V_B_y+\omega{r}\cos\theta\\ \end{array}\right \textit\textbf{{,large\:rotation}} \end{equation*}$

We can use Equations (6.2) for large rotations. However, for small rotations (i.e., $\theta\ll1$ ), we can linearize these relations by substituting for $\sin\theta\approx{0}$ and $\cos\theta\approx{1}$ , or

(6.3) $\begin{equation*} \begin{cases} V_A_x=V_B_x\\ V_A_y=V_B_y+\omega{r}\\ \end{cases} \quad , small \: rotation \end{equation*}$

For example, for a car chassis, the front and rear ends are moving with the same speed as that of the centre of mass, and their velocities in the vertical direction are the algebraic sum of the vertical speed and that of due to the pitch rotation (i.e., rotation about the axis parallel to the ground and perpendicular to the direction of the motion). See the example in section 6.7.

In the following sections, we present examples of several mechanical systems along with their BG models.

6.4 Example: Gear-Shaft Mechanical System—Rotational

The sketch in Figure 6‑2 shows a system composed of four gears, three shafts, six bearings, and two discs.

Build the BG model for this system, including bearings. Use 20-sim.
Identify derivative causalities and the related elements. Discuss the reasoning and how to remove the derivative causalities.
Remove all bearings from the model and perform some analysis using the data provided.

For gears, use velocity ratio, $VR$ equal to the inverse ratio of the number of teeth, $N$ or gears’ diameters, $D$ given by $VR=\dfrac{\omega1}{\omega2}=\dfrac{D2}{D1}=\dfrac{N2}{N1}$ . Angular velocity of gears is represented by the symbol $\omega$ . System data is given in Table 6‑2.

Figure 6-2 Gear and shaft mechanical system

Table 6‑2 Data for the gear-shaft mechanical system
Shaft torsional stiffness (MN.m/rad)	Gear/disc rotational inertia (kg.m2)	Gear number of teeth
K1=2500	J1=40	N/A
K1=2500	J2=15	120
K2=1000	J3=80	400
K2=1000	J4=20	200
K3=700	J5=25	150
K3=700	J6=35	N/A

Solution:

This system has six distinct angular velocities related to gears and discs. The torsional shafts are potential energy storages, and the torsional inertia are kinetic energy storages. The BG elements required are $I$ , $C$ , $R$ , $S_e$ , $S_f$ , and 1- and 0-junctions.

The following video shows how to build and run the model for this example in 20-sim.

One or more interactive elements has been excluded from this version of the text. You can view them online here: https://pressbooks.bccampus.ca/engineeringsystems/?p=66#oembed-12

The resulted BG models are shown in Figure 6‑3 and Figure 6‑4.

Figure 6-3 BG model for a gear shaft system, built in 20-sim

Figure 6-4 BG model for a gear shaft system, derivative causalities removed

6.5 Example: Double Rack-and-Pinion Mechanical System—Rotational

Figure 6‑5 shows a double rack-and-pinion mechanical system. Build a BG model for this system using 20-sim. A torque is applied on the disc connected to the two shafts.

Figure 6-5 A double rack-and-pinion mechanical system

The following video shows how to build and run the model for this example in 20-sim.

One or more interactive elements has been excluded from this version of the text. You can view them online here: https://pressbooks.bccampus.ca/engineeringsystems/?p=66#oembed-13

Figure 6‑6 shows the BG model for this system.

Figure 6-6 BG model for the double rack-and-pinion system

6.6 Example: Mass-Spring-Damper System on an Inclined Plane—Translational

Figure 6‑7 show a mass-spring-damper system on an inclined plane. Build a BG model for this system using 20-sim. Build a BG model for this system using 20-sim.

Figure 6-7 A mass-spring-damper system on an inclined plane

The following video shows how to build and run the model for this example in 20-sim.

One or more interactive elements has been excluded from this version of the text. You can view them online here: https://pressbooks.bccampus.ca/engineeringsystems/?p=66#oembed-14

Figure 6‑8 shows the BG model for this system.

Figure 6-8 BG model for the mass-spring-damper system on an inclined plane

6.7 Example: 2D Rigid-Body Motion—Half-Car Model

In this example, we demonstrate how to build a BG model for a half-car model as shown in Figure 6‑9. The chassis of the car is modelled as a rigid body with two degrees of freedom. The vertical displacement of the centre of mass is the heave, and its angular velocity is the pitch rate. In the BG model, transformer elements are used to transfer the front and rear velocities to the corresponding connecting points between the suspensions and the chassis. The suspension are modelled as spring-dampers and the tires as mass-spring subsystems.

Figure 6-9 Half-car mechanical system sketch

Below are two videos (parts 1 and 2) showing how to build and run the model for this example in 20-sim, including the implementation of the BG transformer element and the setup of the equation model.

Watch the videos and practice building the model on your own, with modified parameters and input signals.

One or more interactive elements has been excluded from this version of the text. You can view them online here: https://pressbooks.bccampus.ca/engineeringsystems/?p=66#oembed-15

One or more interactive elements has been excluded from this version of the text. You can view them online here: https://pressbooks.bccampus.ca/engineeringsystems/?p=66#oembed-16

Figure 6‑10 shows the corresponding BG model. The compression force is considered to be positive.

Figure 6-10 BG model for a half-car mechanical system

6.8 Example: Mass-Spring-Damper System Connected to a Massless Lever

In this example, we demonstrate how to build a BG model for the mechanical system shown in Figure 6‑11. The lever is represented with a $TF$ -element.

Figure 6-11 A mass-spring-damper mechanical system attached to a lever

Below are two videos (parts 1 and 2) showing how to build and run the model for this example in 20-sim, including the implementation of the BG transformer element and the setup of the equation model.

One or more interactive elements has been excluded from this version of the text. You can view them online here: https://pressbooks.bccampus.ca/engineeringsystems/?p=66#oembed-17

One or more interactive elements has been excluded from this version of the text. You can view them online here: https://pressbooks.bccampus.ca/engineeringsystems/?p=66#oembed-18

6.9 Example: Mass-Spring-Damper System Connected to a Lever

For this example we discuss and demonstrate how to build a BG model for the mechanical system as shown in Figure 6‑12. The lever is represented with a $TF$ -element.

Figure 6-12 A mass-spring-damper mechanical system attached to a beam

Below are two videos (parts 1 and 2) showing how to build and run the model for this example in 20-sim, including the implementation of the BG transformer element and the setup of the equation model.

One or more interactive elements has been excluded from this version of the text. You can view them online here: https://pressbooks.bccampus.ca/engineeringsystems/?p=66#oembed-19

One or more interactive elements has been excluded from this version of the text. You can view them online here: https://pressbooks.bccampus.ca/engineeringsystems/?p=66#oembed-20

6.10 Example: Inclined Lever and Mass-Spring-Damper System

In this example, we demonstrate how to build a BG model for a mechanical system consisting of two moving masses attached to a rod, as shown in Figure 6‑13. The rod can rotate as a lever and is represented with a $TF$ -element.

Figure 6-13 Two moving mass-spring system attached to a lever

The following video shows how to build and run the model for this example in 20-sim.

One or more interactive elements has been excluded from this version of the text. You can view them online here: https://pressbooks.bccampus.ca/engineeringsystems/?p=66#oembed-21

6.11 Example: A Pulley-Mass-Spring System

In this example, we demonstrate how to build a BG model for a mechanical system consisting of two pulleys and three masses, as shown in Figure 6‑14.

Figure 6-14 A two-pulley mechanical system

The video below shows how to build and run the model for this example in 20-sim.

One or more interactive elements has been excluded from this version of the text. You can view them online here: https://pressbooks.bccampus.ca/engineeringsystems/?p=66#oembed-22

Exercise Problems For Chapter 6

Exercises

Repeat the BG model using 20-sim for the half-car system given in section 6.7, considering following cases, inclusively and suggested data: , (kg), , (m), (kg.m²/rad), , , (kN/m), , (kN.s/m)
1. ground velocity signal as a pulse signal with start time at 1.5 sec., stop time at 3 sec., and amplitude of 10 cm/s
2. ground velocity signal as a step signal with start time at 1.5 sec. and amplitude of 10 cm/s
3. ground velocity signals such that the car hits a trapezoidal-shaped bump at 1.5s of time, reaches the top of the ramp at 3.5s, stays on top of the plateaued bump for 2s, and comes back to ground level with a similar slope. The amplitude of the velocity signals is 10 cm/s.
4. ground velocity signal as a pulse wave signal with interval of 1 sec., pulse length of 0.1 sec., and amplitude of 5 cm.
5. For all inputs, graph the displacements of the tires and the heave and pitch of the car chassis.
Using 20-sim, build a BG model for the mechanical system shown in the sketch below. Use the data to graph the displacement of the mass $M_s$ , $M_4$ , and spring $k_6$ . Consider the floor velocity input as a pulse signal with start time of 3 sec., stop time of 4.5 sec., and amplitude of 10 cm. Compression forces are considered to be positive (+C). Gravity direction and positive displacements are shown in the sketch.

Masses (kg)

$M_s=30$

$M_1=35$

$M_2=6$

$M_3=15$

$M_4=4$

—

Springs (N/m)

$K_1=18000$

$K_2=50000$

$K_3=190000$

$K_4=135000$

$K_5=20000$

$K_6=10000$

Dampers (N.s/m)

$B_1=400$

$B_2=2500$

$B_3=900$

$B_4=250$

$B_5=350$

$B_6=200$

Bond Graph Models for Electrical Systems

7.1 Overview

An electrical system may consist of components such as resistors, inductors, capacitors, transformers, and batteries/source. The generalized BG elements and relations apply to the analysis of dynamics of electrical systems and components in a similar way that the mechanical components were treated; i.e., they are analogous (see Table 3-1). In other words, electric charge $Q$ is equivalent to the generalized displacement $q$ , electrical current to $i=\dfrac{dQ}{dt}=\.{Q}$ to the flow $f$ , and voltage $V$ to the effort $e$ . The inductor (with inductance $L$ ) is analogous to point mass and is represented by $I$ -element; the capacitor (with capacitance $c$ ) is analogous to a mechanical spring and is represented by $C$ -element; and resistor (with resistance $R$ ) is analogous to mechanical damper and is represented by $R$ -element. The generalized momentum $p_E$ or flux linkage is the integral of $V$ with respect to time. Therefore, we can write

$\begin{equation*} \begin{dcases} e\equiv V, voltage\\ f \equiv i=\.{Q}, current\\ e \cdot{f} \equiv{V} \cdot{i}, power\\ q \equiv{Q}=\int\textit{i dt}, charge\\ p \equiv{P_E}=\int{V}{dt}, \textit{electrical momentum} \end{dcases} \end{equation*}$

Using the constitutive relations, for an $I$ -element we have $e=I\.{f}$ or $V=L\dfrac{di}{dt}$ , and for a C-element we have $e=q/c$ or $V=Q/c$ . Similarly, for an R-element we have $e=Rf$ or $V=Ri$ . The energy associated with storage elements can be written using Equations (3.7) and (3.8), or for elements I and C, as $\dfrac{1}{2L} p^2_E=\dfrac{L}{2} i^2$ and $\dfrac{1}{2c} Q^2=\dfrac{c}{2} V^2$ , respectively. An advantage of the bond graph method is its analogous applicability to different domains using the common constitutive relations, as described above for electrical systems.
The analogy between mechanical and electrical systems can be summarized as follows:

For a system with series-connected components, we have equal effort for mechanical and equal flow for electrical systems. For example, when a spring and a damper are connected in series, they experience the same force, and when a capacitor and a resistor are connected in series, they experience the same current.
For a system with parallel-connected components, we have equal flow for mechanical and equal effort for electrical systems. For example, when a spring and a damper are connected in parallel, they experience the same velocity (or rate of displacement), and when a capacitor and a resistor are connected in parallel, they experience the same voltage.

In other words, the relations of efforts and flows are swapped according to the type of the physical system between mechanical and electrical systems.

$\begin{multiline} mechanical|_{parallel} \: \equiv \: electrical|_{series}\\ \: \\ mechanical|_{series}\:\equiv \: electrical|_{parallel} \end{multiline}$

Table 7-1 shows typical components for resistor/ $R$ , capacitor/ $C$ , Inductor/ $I$ , Transformer/ $TF$ , Electric motor/ $GY$ .

Table 7‑1 Typical electrical components and their corresponding BG elements
$R$ -Element (resistor)	$C$ -Element (capacitor)	$I$ -Element (inductor)	$TF$ -Element (transformer)	$GY$ -Element (electric motor)

7.2 Example: Sign Convention for Electrical Systems

Like mechanical systems, for which we defined +C and +T sign convention for internal forces, we require to define a sign convention for electrical systems. It is customary to use the passive sign convention (PSC) for defining the direction of electrical current ( $i=\dfrac{dQ}{dt}$ ) passing through the elements of an electrical circuit. The background for the PSC is to have power being positive when absorbed by passive elements, e.g., $R$ -, $C$ -, and $I$ -elements in BG method. Therefore, for a typical passive element, by definition, the electrical current is considered as being positive when input into the element from its higher-voltage node (i.e., positive voltage/ $+V$ ) and output from the relatively lower-voltage node (i.e., negative voltage/ $-V$ ). Otherwise, the current is negative. See Figure 7‑1.

Figure 7-1 Sign convention for electrical current through passive elements, passive sign convention

Using the PSC, we have power defined positive for positive current and negative for negative current, or $power=Vi>0$ when $V>0$ and $i>0$ ; hence, power is absorbed by the element. Otherwise, power is generated, when $\textit{power}=Vi<0$ when $V>0$ and $i<0$ . Figure 7‑2 shows the electrical power sign convention for passive elements ( $R$ , $C$ , $L$ ) and active elements (voltage and current sources).

Figure 7‑2 Electrical power sign for several elements according to passive sign convention

In the next section, we use the PSC for defining the current and voltage signs and discuss the step-by-step procedure for building BG models for electrical systems.

7.3 Guidelines for Drawing BG for Electrical Systems

Gustav Robert Kirchhoff (1824–1887)

As mentioned in chapter 4, the general guidelines for drawing BG model can be applied to electrical systems along with causality assignment rules. For electrical systems, we follow these guidelines, along with Kirchhoff’s circuit laws [24] and the PSC for building their BG models, as described in the following steps:

Assign voltage polarity ( $+V, -V$ ) for each element in the electrical circuit.
Assign current direction based on PSC for each element (see Figure 7‑1 and Figure 7‑2).
Assign 0-junction for each distinct voltage node in the circuit, according to Kirchhoff’s voltage law (KVL)—the algebraic sum of all voltage drops around a closed circuit is equal to zero.
Assign 1-junction for each element in the circuit, according to Kirchhoff’s current law (KCL)—the algebraic sum of all electrical currents entering and leaving a node is equal to zero). This is for taking care of relative voltage or drops related to each element located between two 0-junctions, since 1-junction is an effort summator.
Select a node in the circuit as a reference, i.e., the grounding, with zero voltage.
Assign $C$ -element for capacitors, $R$ -element for resistors, $I$ -element for inductors, $S_e$ for voltage, and $S_f$ for current sources.
Assign $TF$ -element for electrical transformers and $GY$ -element for electric motors.
Connect the elements with power bonds, assign causalities, and simplify by neglecting the bonds and the 0-junction which are connected to the ground source.

The above steps are based on KVL, and the process starts with assigning 0-junctions for each distinct voltage node. It is also possible to start with KCL and assign 1-junctions for the current in each closed-circuit loop and use 0-junctions in between for distribution of the current to corresponding circuit loops. The latter will result in a more simplified BG model and is recommended for complex circuits that involve several electric loops. In practice, we sometimes use a combination of these two approaches for building BG model for electrical systems.

In the following sections, we demonstrate the application of the procedure discussed above, with some worked-out examples.

7.4 Example: An RCL Circuit—in Series

Figure 7‑3 shows an RCL circuit consisting of a resistor, a capacitor, and an inductor connected in series. To build the BG model, we apply the PSC and use the procedure listed in section 7.3. The four nodes identified by solid circles have distinct voltages. Therefore, four 0-junctions are assigned at the four corners of the circuit. For voltage drop across each element, we assign 1-juction and connect it to the corresponding element with a power bond. Note that the current direction in the circuit is consistent with the PSC convention. The resulting BG model is shown in Figure 7‑4 after being simplified with deleted ground-connecting bonds shown in the dashed circle. Alternatively, we can simplify the BG model and use a 1-juction for the current in the circuit loop according to KCL. In other words, the electrical current flowing through all elements should be identical. The resulting simplified BG model is shown in Figure 7‑5.

It is useful to discuss the analogy between the RCL circuit and mechanical mass-spring-damper systems (see Figure 4‑1) and their identical BG model. Assuming a ground connection for the circuit is analogous to a wall with zero velocity for the mass-spring-damper system, the current through the inductor is analogous to the velocity of the mass. The same current flows through the resistor and the capacitors, analogous to the velocity of the spring and damper components. Therefore, the simplified BG model (see Figure 7‑5) is identical for both electrical and mechanical systems. In other words, the BG model is identical to the one for a mass-spring-damper connected in parallel.

Figure 7-3 Sketch for a RCL electrical circuit in series

Figure 7‑4 BG model for a RCL electrical circuit in series

Figure 7‑5 Simplified BG model for an RCL electrical circuit in series

7.5 Example: An RCL Circuit—in Parallel

Figure 7‑6 shows an RCL circuit consisting of two inductors, a resistor, and a capacitor connected in parallel. We use the KCL approach to build the BG model for this example. Because the voltages across all components are identical, we can, using power bonds, apply a 0-junction (i.e., voltage equalizer) and connect the $R$ , $C$ , and $I$ components to it. This can be obtained by simplifying the BG model shown in Figure 7‑7.

Figure 7‑6 Sketch for a RCL electrical circuit in parallel

The following video shows how to build and run the model for this example in 20-sim.

One or more interactive elements has been excluded from this version of the text. You can view them online here: https://pressbooks.bccampus.ca/engineeringsystems/?p=68#oembed-7

Figure 7‑7 BG model for a RCL electrical circuit in parallel

7.6 Example: An Electrical Circuit—Two Loops

Figure 7‑8 shows an RCL two-loop circuit consisting of resistors, inductors, and a capacitor connected in parallel. We use the KCL approach to build the BG model for this example.

Figure 7‑8 A two-loop RCL electrical circuit

The following video shows how to build and run the model for this example in 20-sim.

One or more interactive elements has been excluded from this version of the text. You can view them online here: https://pressbooks.bccampus.ca/engineeringsystems/?p=68#oembed-8

The simplified BG model with a supplied voltage signal as a square wave is shown in Figure 7‑9.

Figure 7‑9 BG model for the two-loop RCL electrical circuit

7.7 An Electrical Circuit—Three Loops

Figure 7‑10 shows an RCL three-loop circuit consisting of resistors, inductors, and a capacitor connected in parallel. We use the KCL approach to build the BG model for this example.

Figure 7‑10 A three-loop electrical circuit

The following video shows how to build and run the model for this example in 20-sim.

One or more interactive elements has been excluded from this version of the text. You can view them online here: https://pressbooks.bccampus.ca/engineeringsystems/?p=68#oembed-9

The simplified BG model with a supplied voltage signal as a block wave is shown in Figure 7‑11.

Figure 7‑11 BG model for the three-loop RCL electrical circuit

7.8 An Electrical Circuit—Wheatstone Bridge

Georg Simon Ohm (1789–1854)

Figure 7‑12 shows a Wheatstone circuit consisting of resistors. This circuit is usually used to measure an unknown resistor, e.g., placed in the system as $R_{4}$ , by adjusting the variable $R_{2}$ such that the current through $R_{L}$ is null, i.e., the balanced point. Using Kirchhoff’s and Ohm’s laws [25], we can calculate the currents going through the branch $abd$ as $\dfrac{V}{R_{1}+R_{2}}$ and branch $acd$ as $\dfrac{V}{R_{3}+R_{4}}.$ Therefore, the voltages at nodes $b$ and $c$ , with reference to the ground, are $V_{b}=\dfrac{V}{R_{1}+R_{2}}R_{2}$ and $V_{c}=\dfrac{V}{R_{3}+R_{4}}R_{4}$ , respectively. For having null voltage across $R_{L}$ , we let $V_{b}=V_{c}$ or after some manipulations, we get $R_{4}=\dfrac{R_{2}}{R_{1}}R_{3}$ . As shown, the balanced point is independent of the voltage supplied. We use the KCL approach to build the BG model for this example.

Figure 7‑12 A Wheatstone bridge electrical circuit

The following video shows how to build and run the model for this example in 20-sim.

One or more interactive elements has been excluded from this version of the text. You can view them online here: https://pressbooks.bccampus.ca/engineeringsystems/?p=68#oembed-10

The simplified BG model with a supplied voltage is shown in Figure 7‑13.

Figure 7‑13 BG model for the Wheatstone bridge circuit

7.9 An Electrical Circuit—Multi-loop

Figure 7‑14 shows an RCL multi-loop circuit consisting of resistors, inductors, and capacitors connected in series and parallel. We use the KCL/KVL approach to build the BG model for this example.

Figure 7‑14 A multi-loop electrical circuit

The following video shows how to build and run the model for this example in 20-sim.

One or more interactive elements has been excluded from this version of the text. You can view them online here: https://pressbooks.bccampus.ca/engineeringsystems/?p=68#oembed-11

The simplified BG model with a supplied voltage is shown in Figure 7‑15.

Figure 7‑15 BG model for the multi-loop electrical circuit

7.10 An Electrical Circuit—Multi-loop with Transformer

Figure 7‑16 shows an RCL multi-loop circuit consisting of resistors, inductors, capacitors, and a transformer connected in series and parallel. We use the KCL/KVL approach to build the BG model for this example.

Figure 7‑16 A multi-loop electrical circuit with transformer

The following video shows how to build and run the model for this example in 20-sim.

One or more interactive elements has been excluded from this version of the text. You can view them online here: https://pressbooks.bccampus.ca/engineeringsystems/?p=68#oembed-12

The simplified BG model with a supplied voltage is shown in Figure 7‑17.

Figure 7‑17 BG model for the multi-loop electric circuit with transformer

Exercise Problems for Chapter 7

Exercises

Build the BG model for the electrical system as shown in the sketch. Run the model and report the following quantities:
1. charge accumulated on capacitors
2. current across resistors
3. voltage drop across resistor $R_2$
4. momentum (flux linkage) for the inductor.

Use following data: $V(t)=100\sin(20t)$ , $R_1=0.1k\Omega$ , $R_2=0.5k\Omega$ , $R_3=0.8k\Omega$ , $C_1=200\mu{F}$ , $C_2=400\mu{F}$ , $L_1=5mH$ , and transformer parameter 2:1. Perform Parameter Sweep on a range of $TF$ parameter values, 0.5-3 and graph the results for electric charge on capacitor $C_2$ .

Build a BG model for the electrical circuit shown. Use $R=3k\Omega$ , $C=1\mu{F}$ , $L=100mH$ for simulation. Report voltages across each element for a direct source voltage of $5V$ . Also, run the model for a range of capacitance $10nF$ , $100nF$ , $1000nF$ using Parameter Sweep and report the across the inductance for these values. Draw the sketch.

For the electrical system shown in the sketch, build the BG model.

A modified Wheatstone bridge circuit is shown in the sketch. Build a BG model and show that the voltage across the bridge resistor (R5) is null when the bridge is balanced.

An electrical circuit is shown in the below sketch below. The circuit consists of two capacitors, two inductors, and one resistor. Build the corresponding BG model.

Modify the example given in section 7.7 by making the components branching from node $a$ to be laid out in parallel. Build the BG model for the modified circuit.
For the example given in above section 7-9, use the corresponding BG model and the following data to simulate the system: $V(t)=60\sin(15t)$ , $R_1=0.1k\Omega$ , $R_2=0.5k\Omega$ , $C_1=200\mu{F}$ , $C_2=400\mu{F}$ , $C_3=300\mu{F}$ , $L_1=5mH$ , $L_2=10mH$ .
Build the BG model for the electrical circuit shown below. After building the model in 20-sim, simplify it and interpret the simplified model. Perform a parametric sweep analysis for the capacitor and inductor.

Media Attributions

Gustav Robert Kirchhoff is licensed under a Public Domain license
Exercise-7-1

Bond Graph Models for Hydraulic Systems

8.1 Overview

The generalized BG elements and relations apply to the modelling of dynamics of hydraulic systems in a similar way that the mechanical or electrical systems were treated; i.e.; they are analogous (see Table 3‑1). In this chapter, we define the effort and flow for hydraulic systems and derive the relations for hydraulic capacitance, inertance, and resistance corresponding to BG elements $C$ , $I$ , and $R$ , respectively. Note that the complexity of fluid behaviour in static or dynamic flow conditions require us to pay more attention to identify these quantities and relations as compared to those for mechanical and electrical systems.

For modelling hydraulic systems, we are usually interested in having a relationship between pressure and fluid volume in static conditions and between pressure drop and fluid volume flow rate in dynamic conditions. For example, we might be interested to know the pressure drop for a given flow rate in a pipe, or we might want to know the pressure at a given depth in a storage tank, as sketched in Figure 8‑1.

Figure 8-1 Sketches for pressure drop in a pipe and in a storage tank

8.2 Definitions of Effort, Flow, and Momentum for Hydraulic Systems

Recall that power is the quantity of interest in BG method (see section 3.2). Considering a control volume $\mathbb{V}$ of an incompressible fluid flowing under pressure $\mathbb{p}$ , we can write the power $\mathcal{P}$ as the product of the force $F$ exerting on the fluid, resulted from applied pressure, and the velocity $v$ of the fluid flowing through the volume, or $\mathcal{P}=Fv.$ But the velocity of the fluid can be written as $v={\dot{\mathbb{V}}}/{\mathbb{A}}$ , using the continuity relation, where $\dot{\mathbb{V}}$ is volume flow rate of the fluid and $\mathbb{A}$ is the cross-sectional area of the control volume. Therefore, after substitution, we get $\mathcal{P} = F \dfrac{\dot{\mathbb{V}}}{\mathbb{A}} = \dfrac{F}{A} \dot{\mathbb{V}} = \mathbb{p}\dot{\mathbb{V}}$ , or equivalently rate of energy $\mathbb{pV}$ . Comparing the relation $\mathcal{P} = \mathbb{p}\dot{\mathbb{V}}$ with the BG generalized relation for power, i.e., $\mathcal{P} = ef$ , we can write $e \equiv \mathbb{p}$ and $f \equiv \dot{\mathbb{V}}$ . In other words, for hydraulic systems, pressure is equivalent to BG effort, and fluid volume flow rate is the BG flow. Similarly, we can write the generalized BG displacement $q$ as the volume of the fluid, or $q = \int f dt = \int \dot{\mathbb{V}} dt = \mathbb{V}$ .

In BG method, the generalized momentum is the integral of effort. Therefore, we can write the fluid momentum $p_{\mathbb{P}}$ as the integral of pressure, or $p_{\mathbb{p}} = \int \mathbb{p} dt$ . Summarizing these relations, we have

$\begin{equation*} \begin{dcases} e \equiv \mathbb{p}, pressure \\ f \equiv \dot{\mathbb{V}} , volume \: flow \: rate \\ e \cdot f \equiv \mathbb{p} \dot{\mathbb{V}}, power \\ q \equiv \mathbb{V} = \int \dot{\mathbb{V}}dt, volume \\ p \equiv p_{\mathbb{p}} = \int \mathbb{p} dt, pressure \: momentum \\ \end{dcases} \end{equation*}$

For hydraulic systems, fluid pressure is equivalent to generalized BG effort, and fluid volume flow rate is the generalized BG flow.

8.3 Fluid Compliance: C-element

Fluid compliance or hydraulic capacitance describes potential energy storage with a fluid, e.g., the height of fluid in a tank. It is equivalent to mechanical spring compliance or electrical capacitor capacitance. For a $C$ -element in BG method, we have $e=q/c$ . Using equivalent quantities for hydraulic systems, we can write $\mathbb{p} = \mathbb{V}/c$ , or fluid compliance is volume change per unit of pressure acting on the fluid volume. For an incompressible fluid with density $\rho$ , the hydrostatic pressure at depth $h$ is $\mathbb{p} = \rho gh$ and the volume of the fluid is $\mathbb{V} = \mathbb{A}h$ . After substitution, we get $\mathbb{p} = \mathbb{A}h/c$ , or after rearranging and simplifying, the hydraulic capacitance for incompressible fluid is

(8.1) $\begin{equation*} c = \dfrac{\mathbb{A}}{{\rho}g} , incompressible \end{equation*}$

where $g$ is the gravitational acceleration. The dimension of fluid compliance can be worked out as $[\dfrac{m^4s^2}{kg}]$ = [ $\dfrac{m^5}{N}$ ] = $\dfrac{volume}{pressure}$ .

Note that the pressure could be replaced by total dynamic pressure for fluid in motion.

If the fluid is compressible, we use the bulk modulus of elasticity $\beta$ for calculating the change in volume. By definition, $\beta$ is pressure needed to change fluid volume per unit of volume, or $\beta = \vartriangle \mathbb{p}/(\vartriangle\mathbb{V}/\mathbb{V})$ . Therefore, $c = \dfrac{\vartriangle\mathbb{v}}{\vartriangle\mathbb{p}} = \dfrac{\mathbb{v}}{\beta}$ .

(8.2) $\begin{equation*} c = \dfrac{\mathbb{V}}{\beta} , compressible \end{equation*}$

For more complex flow and non-uniform, flexible tubes, consult with chapter 4 of Dean, Karnopp, Margolis, and Rosenberg [20].

Having the hydraulic capacitance, we can write the relation between the flow rate and the pressure as $\vartriangle\dot{\mathbb{V}} = \dot{\mathbb{V}}_2 - \dot{\mathbb{V}}_1 = c\dfrac{d\mathbb{p}}{dt}$ , useful to calculate the flow rate for given pressure. Similarly, we can write $\mathbb{p} = \dfrac{1}{c}\int\vartriangle\overset{.}{\mathbb{V}}dt = \dfrac{1}{c}\int( \dot{\mathbb{V}}_2 - \dot{\mathbb{V}}_1)dt$ , useful for calculating pressure for given flow rates.

(8.3) $\begin{equation*} \begin{dcases} \vartriangle\dot{\mathbb{V}} = \dot{\mathbb{V}}_2 - \dot{\mathbb{V}}_1 = c\dfrac{d\mathbb{p}}{dt}\\ \mathbb{p} = \dfrac{1}{c}\int( \dot{\mathbb{V}}_2 - \dot{\mathbb{V}}_1)dt\\ \end{dcases} \end{equation*}$

Note the similarity between relations given by Equation (8.3) and those given for mechanical spring (when pressure is replaced by force and fluid volume flow rate by velocity) and electrical capacitance (when pressure is replaced by voltage and fluid volume flow rate by current).

8.4 Fluid Inertia: I-element

Fluid inertia, or hydraulic inertance, describes kinetic energy storage with a fluid or the inertia, e.g., of a fluid flowing in a pipe. It is equivalent to inertia related to mass in mechanical or inductance in electrical systems. For an $I$ -element in BG method, $p = If$ describes the relation between generalized momentum and flow. Using equivalent quantities for hydraulic systems, we can write $p_\mathbb{p} = I\dot{\mathbb{V}}$ , or fluid pressure momentum is the product of fluid inertia by its volume flow rate. To derive the relation for $I$ , we require to have the relationship between the momentum and volume flow rate of the fluid flow. For derivation, we consider a control volume with length $\ell$ and cross-sectional area $\mathbb{A}$ of the fluid with density $\rho$ , as shown in Figure 8-2.

Figure 8-2 Sketch for a control volume of flowing fluid in a pipe

Assuming a pressure difference $\vartriangle\mathbb{p}$ between two ends of the control volume acting on the fluid, we can write Newton’s second law for the fluid motion as $\underbrace{\vartriangle\mathbb{pA}}_{force} = \underbrace{\rho\mathbb{V}}_{mass}\dfrac{d}{dt}$ $\underbrace{\left(\dfrac{\dot{\mathbb{V}}}{\mathbb{A}}\right)}_{velocity}$ . Rearranging the terms and integrating the pressure, we get the pressure momentum $p_{\mathbb{p}} = \int {\mathbb{p}dt = \dfrac{\rho{\mathbb{V}}}{\mathbb{A}}\dfrac{\dot{\mathbb{V}}}{{\mathbb{A}}}$ . But volume is $\mathbb{V} = \mathbb{A}\ell$ . After substituting, we get the relationship between pressure momentum and the volume flow rate as

(8.4) $\begin{equation*} p_\mathbb{p} = \dfrac{\rho\ell}{\mathbb{A}}\dot{\mathbb{V}} \end{equation*}$

Comparing Equation (8.4) with the generalized momentum equation for $I$ -element, $p = If$ we can write the fluid inertia as

(8.5) $\begin{equation*} I= \dfrac{\rho\ell}{\mathbb{A}}, incompressible \end{equation*}$

The dimension of fluid inertia can be worked out as $[\,\dfrac{kg}{m^4}]\, \equiv \dfrac{pressure\: momentum}{volume\: flow \:rate}$ .

From Equation (8.5), we can conclude that a fluid has larger inertia when flowing in small diameter tubes, compared to in larger tubes because $I$ is inversely proportional to $\mathbb{A}$ . This effect is counterintuitive and is sometimes misinterpreted with the wrong assumption that large-size tubes should exhibit larger inertia effects. Note that we consider only force due to pressure, and not friction due to viscosity, assuming an ideal fluid.

For more complex flow and non-uniform and/or flexible tubes, consult with chapters 4 and 12 of Dean, Karnopp, Margolis, and Rosenberg [20]. For example, if the cross-section of the pipe and the density of fluid change along its $x$ -axis, then we get

(8.6) $\begin{equation*} I = \int\dfrac{pdx}{{\mathbb{A}}(x)} \end{equation*}$

Having the hydraulic inertance, we can write the relation between the flow rate and the pressure as $\vartriangle\dot{\mathbb{V}} = \dot{\mathbb{V}}_2 - \dot{\mathbb{V}}_1 = \dfrac{1}{I}\int\mathbb{p}dt = \dfrac{\mathbb{A}}{\rho\ell}\int\mathbb{p}dt$ , useful to calculate the flow rate for given pressure. Similarly, we can write $\vartriangle\mathbb{p} = I\dfrac{d\dot{\mathbb{V}}}{dt} = \dfrac{\rho\ell}{\mathbb{A}}\dfrac{d\dot{\mathbb{V}}}{dt}$ , useful for calculating pressure for given flow rates.

(8.7) $\begin{equation*} \begin{dcases} \vartriangle\dot{\mathbb{V}} = \dfrac{1}{I}\int\mathbb{p}dt\\ \vartriangle\mathbb{p} = I\dfrac{d\dot{\mathbb{V}}}{dt}\\ \end{dcases} \end{equation*}$

Note the similarity between relations given by Equation (8.7) and those given for mechanical systems (when pressure is replaced by force, fluid flow rate by velocity, and inertance by mass) and electrical systems (when pressure is replaced by voltage, fluid flow rate by current, and inertance by inductance).

8.5 Fluid Resistance: R-Element

Gotthilf Hagen, (1797–1884)

Fluid, or hydraulic, resistance describes energy dissipation with a fluid, e.g., friction of a fluid flowing in a pipe. Fluid resistance is equivalent to dampers in mechanical or resistors in electrical systems. For an $R$ -element in BG method, we have $e=Rf.$ Using equivalent quantities for hydraulic systems, we can write

(8.8) $\begin{equation*} \mathbb{p} = R\dot{\mathbb{V}} \end{equation*}$

or fluid resistance is equal to pressure change per unit volume flow rate. This relationship depends on the state of the flow (e.g., laminar, turbulent) and the fluid properties (e.g., ideal, viscous,) [26], [27], [28]. Note the similarity between relations given by Equation (8.8) and those given for mechanical systems (when pressure is replaced by force and fluid flow rate by velocity) and electrical systems (when pressure is replaced by voltage and fluid flow rate by current).

Jean-Léonard-Marie Poiseuille, (1797–1869)

To demonstrate the derivation of the relation for $R$ , we consider a laminar flow of a viscous incompressible fluid in a pipe (so-called Hagen-Poiseuille flow) and write Newton’s second law for a cylindrical differential control volume of the fluid with length $dz$ along the pipe axis and a cross-section with radius $r$ , as shown in Figure 8‑3. This flow is axisymmetric, and the velocity profile changes along the radius related to a cylindrical coordinate system $(r, \theta, z).$

Figure 8-3 Sketch of velocity profile for a Hagen-Poiseuille flow in a pipe

For a steady flow (i.e., non-transient), we have $\sum$ $F_z = 0$ . But the forces applied on the fluid are due to pressures $\mathbb{p}_1$ at point $z_1$ and $\mathbb{p}_2 = (\mathbb{p}_1 - d\mathbb{p}) < \mathbb{p}_1$ at point $z_2$ and the viscous-induced shear stress $\tau$ . Hence, $(\mathbb{p}_1 - \mathbb{p}_1 + d\mathbb{p}) \pi r^2 - \tau(2\pi{r})dz = 0$ where $r$ is the radial dimension in the $x-y$ plane parallel to the pipe cross-section. We need the relation for fluid friction effect due to viscosity. According to Newton’s law for a viscous fluid, we have $\tau = \mu\dfrac{du(r)}{dr}$ , assuming the shear stress due to the fluid’s viscosity be proportional to the velocity gradient along the radius with the proportionality constant being the dynamic viscosity $\mu$ . Note that velocity profile at any cross-section of the pipe is only a function of radius, or velocity vector is $\overrightarrow{u} = (0, 0, u(r))\equiv (0, 0, u)$ . After substitution, we have $\pi r^2 d\mathbb{p} = \mu\dfrac{du}{dr} (2\pi r)dz$ . After simplifying and rearranging the terms, we have $\dfrac{r}{2\mu}d\mathbb{p} = \dfrac{du}{dr}dz$ . Integrating the latter relation, noting that velocity is not a function of $z$ , gives $\dfrac{r}{2\mu}\displaystyle\int_{\mathbb{p}_1}^{\mathbb{p}_2}d\mathbb{p} = \dfrac{du}{dr}\displaystyle\int_{z_1}^{z_2}dz$ , or $\dfrac{r}{2\mu}({\mathbb{p}_2}-{\mathbb{p}_1}) = \dfrac{du}{dr}(z_2 - z_1)$ . Now, we rearrange the terms and let $z_2 - z_1 = \ell$ , the length of the control volume, and use the pressure difference as a positive constant quantity $\mathbb{p} = (\mathbb{p}_1 - \mathbb{p}_2) > 0$ in the direction of the fluid flow. Hence, $du = -\dfrac {rdr}{2\mu\ell}\mathbb{p}$ . Integrating both sides (the left-hand side with respect to $u$ and the right-hand side with respect to $r$ ) gives $u = -\dfrac{r^2}{4\mu\ell}\mathbb{p} + constant$ . The constant of integration can be obtained using the information at the boundary of the pipe assuming the no-slip condition, or $u|_{r=D/2} = 0$ , where D is the pipe diameter. Therefore, $constant = \dfrac{\mathbb{p}D^2}{16\mu\ell}$ . Hence, after back substitution, we get the velocity profile $u = - \dfrac{\mathbb{p}r^2}{4\mu\ell} + \dfrac{\mathbb{p}D^2}{16\mu\ell}$ . This relation is the famous parabolic velocity profile for the flow in a pipe and can be written in its functional form as

(8.9) $\begin{equation*} u(r) = \dfrac{\mathbb{p}}{4\mu\ell}\left({\dfrac{D^2}{4} - r^2}\right) \end{equation*}$

Using Equation (8.9), we can calculate the velocity for any given value of $r$ , e.g., $u|_{r=0} = \dfrac{\mathbb{p}D^2}{16\mu\ell}$ at the centre-line of the pipe or $u|_{r=D/2} = 0$ at the interior wall of the pipe.

Now, to find the volume flow rate, we integrate the velocity $u(r)$ over the whole cross-section of the pipe using a differential area element $d\mathbb{A} = rdrd\theta$ . Or, the volume of the fluid passing through the whole cross-section of the pipe per unit of time is given by

$\dot{\mathbb{V}}=\iint ud{\mathbb{A}}=\int_o^{2\pi} d\theta \int_0^{D/2}urdr = 2\pi \dfrac{\mathbb{p}}{4\mu\ell}\int_o^{D/2}\left({\dfrac{D^2}{4} - r^2}\right)rdr = \dfrac{\pi{D^4}}{128\mu\ell}\mathbb{p}$ , or

(8.10) $\begin{equation*} \dot{\mathbb{V}} = \dfrac{\pi{D^4}}{128\mu\ell}\mathbb{p} \end{equation*}$

Comparing Equation (8.10) with Equation (8.8), we can write the fluid resistance as

(8.11) $\begin{equation*} R = \dfrac{128\mu\ell}{\pi{D^4}} \end{equation*}$

The fluid resistance can be interpreted as the amount of pressure drop per unit of volume flow rate of the fluid in the pipe. The dimension of fluid resistance can be worked out as $\Big[\dfrac{kg}{s.m^4}\Big] = \Big[\dfrac{N.s}{m^5}\Big]\equiv \dfrac{pressure}{volume\: flow\: rate}$ .

Other BG elements for hydraulic systems are sources of flow (e.g., centrifugal pumps) and efforts (e.g., reservoirs, tanks, displacement pumps). Pumps provide flow of a fluid at a certain flow rate according to their types and specifications. Reservoirs or pressure chambers provide certain pressure to the system as an effort source. The transformers elements are those like piston-cylinder (plunger), and gyrators are those elements like reaction turbines or hydraulic motors. Sketches below show some related elements.

Typical hydraulic components are shown in Table 8‑1.

Table 8-1 Typical hydraulic components and their corresponding BG elements
$R$ -element (valve)	$C$ -element (storage)	$I$ -element (fluid mass)	$TF$ -element (plunger)	$GY$ -element (pump)

8.6 Sign Convention for BG Modelling of Hydraulic Systems

The sign convention for hydraulic systems can be defined by specifying the relative high/low pressure points in the system and, hence, the positive fluid flow direction along the pressure drop. The pressure reference is commonly taken to be the atmospheric pressure (i.e., one atm for absolute and zero for gauge pressures). For BG modelling, it is recommended to have all pressures in gauge and define a zero-pressure point for reference atmospheric pressure. If the results are required in absolute pressure units, then one unit of atmospheric pressure can be added to the obtained values from the BG model.

8.7 Guidelines for Drawing BG for Hydraulic Systems

As mentioned in chapter 4, the general guidelines for drawing BG models can be applied to hydraulic systems, along with causality assignment rules. For hydraulic systems, we follow the guidelines given for electrical systems (see section 7.3) as described in the following steps:

1) Assign sign convention for fluid flow directions.

2) Assign 0-junction for each distinct pressure point in the system.

3) Assign 1-junction for each element in the system. This is for taking care of relative pressure drops related to each element located between two adjacent 0-junctions, since 1-junction is effort summator.

4) Select a node in the system as a reference, i.e., the atmospheric pressure point, and assign a 0-junction element to it. If gauge pressures are used, then this 0-junction and all connected power bonds can be eliminated to simplify the model.

5) Assign $C$ -element for storage/capacitors, $R$ -element for friction, $I$ -element for fluid mass, and $S_e$ for pressure and $S_f$ for flow sources.

6) Assign $TF$ -element for hydraulic transformers and $GY$ -element for hydraulic gyrators.

7) Connect the elements with power bonds and assign causalities. Simplify by neglecting the bonds and the 0-junction which are connected to the 0-junction representing the atmospheric pressure.

Similarly, a 1-junction-based approach can be used for distinct flow rates and hence simplifying the BG model, as we demonstrated in the previous chapter with electrical systems.

8.8 Example: Hydraulic Reservoir-Valve System

Figure 8‑4 shows a hydraulic system consisting of two tanks, pipes, and valves. Build a BG model for this system.

Figure 8-4 Sketch for a two-tank hydraulic system

Solution:

We follow the systematic procedure for building the BG model, listed in section 8.7. For this system, we can easily identify two pressure points located at the bottom of tanks. We assign two 0-junctions for each. For flow input, we assign a flow source element, and for the output, an effort source element to define the atmospheric pressure at that location. For the tanks, we only consider capacitance, assuming slow fluid motion and neglect inertia and friction (i.e., no inertance nor resistance). For the pipe sections, we consider inertance and resistance. As well, we assign 1-junctions for flows in the pipes that represent the pressure changes for these components. Figure 8-5 shows the resulting BG model.

Figure 8-5 Bond graph model for a two-tank hydraulic system

8.9 Example: Hydraulic Reservoir-Valve System Simulation

In this example, we use the BG model developed in section 8.8, along with data assigned to parameters for simulation. Considering water as the fluid ( $\rho = 10^3kg/m^3 , \mu = 10^-3Pa.s$ ) and the data given in Table 8‑2, we can calculate the related $C$ , $R$ , and $I$ of the elements in the system. The diameter of the pipes is 15 cm, and $g = 9.81\dfrac{m}{s^2}$ .

The following video shows how to build and run the model for this example in 20-sim.

One or more interactive elements has been excluded from this version of the text. You can view them online here: https://pressbooks.bccampus.ca/engineeringsystems/?p=70#oembed-4

Table 8‑2 Data for system components in example 8-9
Component	$X$ -section area [ $m^2$ ]	Length [ $m$ ]	$C$ [ $m^5/N$ ] Eq. (8.1)	$I$ [ $kg/m^4$ ] Eq. (8.5)	$R$ [ $N.s/m^5$ ] Eq. (8.11)
Storage tanks	2	–	$2.04 X 10^4$	–	–
Pipe1	0.01767	4	–	226372.4	322
Pipe2	0.01767	2	–	113186.2	161

8.10 Example: Hydraulic Pump-Reservoir-Valve System

In this example, we use the BG model developed in section 8.8, adding a pump to the system as shown in Figure 8‑6. In this example we discuss in more detail the BG model of a pump. For further details related to BG modelling of pumps, consult with references cited as [21] and [29].

Figure 8-6 A hydraulic system with a pump

The following video shows how to build and run the model for this example in 20-sim. The resulting BG model is shown in Figure 8‑7.

One or more interactive elements has been excluded from this version of the text. You can view them online here: https://pressbooks.bccampus.ca/engineeringsystems/?p=70#oembed-5

Figure 8-7 Bond graph model for the pump-reservoir-valve hydraulic system

8.11 Example: A Hydraulic Lift System

Blaise Pascal (1623-1662)

In this example, we consider a hydraulic lift, as sketched in Figure 8‑8. We build a BG model for this hydraulic system. The continuity relation applies to the fluid motion and Pascal’s law defines the pressure distribution of the fluid in the cylinders. Two transformer elements are used in the BG model to convert linear velocities of the pistons to/from volume flow rate and convert forces to pressures (see Figure 8‑9). The transformers’ parameters are explained in the video clip.

Figure 8‑8 A hydraulic lift system

The following video shows how to build and run the model for this example in 20-sim.

One or more interactive elements has been excluded from this version of the text. You can view them online here: https://pressbooks.bccampus.ca/engineeringsystems/?p=70#oembed-6

The BG model is shown in Figure 8‑9, along with the detail of the transformers’ inputs and outputs.

Figure 8-9 Bond graph model for hydraulic lift system

Exercise Problems for Chapter 8

Exercises

Build the bond graph for a two-way safety valve.
Repeat the example in section 8-9 and perform a parametric sweep for some parameters in the simulation, for example pipe diameters and lengths.
Expand the BG model given in section 8.10 with running simulation with some data for the system parameters, similar to those given in section 8.9. Also, expand the model of the pump using some pump-chart (H-Q).
Use some data and run simulation for the example given in section 8.11, the hydraulic lift.

Media Attributions

Jean-Léonard-Marie Poiseuille is licensed under a Public Domain license
fig-8-3_edits
Fig-8-4

Bond Graph Models for Multi-Domain Systems

9.1 Overview

As mentioned in previous sections (see chapter 3), the universality of BG elements for application to analogous quantities is the main strength of this method for modelling and simulation of multi-domain systems. In practice, many systems are composed of multi-energy domains, e.g., an electric motor turning a mechanical shaft, a car engine generating and transferring power to the torsion shaft, a wind turbine, a robot. The BG modelling method can serve as a powerful tool for modelling these real-world systems.

In this chapter, we present several examples of multi-domain systems and build their BG models. We emphasize that the control sub-system is a major part of any engineering system. In this textbook, however, we focus on BG method and how to build BG models for systems. The full treatment of the topic of control is left for a possible future volume. For information about control theory and modelling, consult with references such as those cited as [20], [21], [30], and [31].

9.2 Example: Car Brake System

For this example we consider a car brake system as shown in Figure 9‑1. The driver applies a force on the brake pedal, which is transferred to the brake discs through the hydro-mechanical system. The process of force transfer is modelled with using several transformer elements ( $TF$ ).

Figure 9-1 A car brake hydro-mechanical system

The following video shows how to build and run the model for this example in 20-sim.

One or more interactive elements has been excluded from this version of the text. You can view them online here: https://pressbooks.bccampus.ca/engineeringsystems/?p=72#oembed-5

The BG model for this system is shown in Figure 9‑2. For further reading, consult with the reference cited as [32].

Figure 9‑2 BG model for the car brake hydro-mechanical system

9.3 Example: Electro-mechanical Hoist System

For this example, we consider an electro-mechanical hoist system as shown in Figure 9‑3. The electric motor is connected to a shaft-drum mechanical system. The load is represented by a mass connected by an elastic extensible string to the hoist drum. A gyrator ( $GY$ ) and a transformer ( $TF$ ) elements are used in this model. The $GY$ -element models the electric motor by transforming the motor voltage (current) to the angular velocity (torque) of the drum/shaft and the $TF$ -element transforms the angular velocity of the drum to the linear velocity of the mass.

Figure 9‑3 An electro-mechanical system with load

The following video shows how to build and run the model for this example in 20-sim.

One or more interactive elements has been excluded from this version of the text. You can view them online here: https://pressbooks.bccampus.ca/engineeringsystems/?p=72#oembed-6

Figure 9‑4 shows the BG model for this system. For further reading, consult with the references cited as [21] and [33].

Figure 9‑4 BG model for an electro-mechanical system

9.4 Example: Drive Shaft-Load Mechanical System

For this example, we consider a mechanical drive shaft system carrying a torsional load as shown in Figure 9‑5. The applied torque is transferred to the gear-shaft system. The load can be applied through an electric motor (not shown). This model uses several transformer elements ( $TF$ ). The $TF$ -elements exchange the angular velocity of the gears, using compatibility requirement.

Figure 9‑5 A drive shaft mechanical system carrying a torsional load

The following video shows how to build and run the model for this example in 20-sim.

One or more interactive elements has been excluded from this version of the text. You can view them online here: https://pressbooks.bccampus.ca/engineeringsystems/?p=72#oembed-7

The BG model for this system is shown in Figure 9‑6. For further reading, consult with the references cited as [20] and [34].

Figure 9‑6 BG model for drive shaft mechanical system carrying a torsional load

9.5 Example: Inverted Double Pendulum

For this example, we consider an inverted double pendulum system pivoted to a moving mass, $M$ as shown in Figure 9‑7. The rods are pinned together, and each has one rotational DOF, represented by $\theta_1$ and $\theta_2$ . The centre of mass is located at the mid-points $\dfrac {l_1}{2}$ and $\dfrac {l_2}{2}$ of each rod, represented by $x_1$ and $x_2$ ; associated mass by $m_1$ and $m_2$ ; and rotational moment of inertia by $J_1$ and $J_2$ . The motion is measured with reference to a fixed coordinate system, $x-y$ , initially at origin. Applied force $F(t)$ is exerted on mass $M$ , moving horizontally with its displacement designated by $x$ .

The coordinates of the centre of mass for the rods are given as $x_1=x+\dfrac {l_1}{2} \sin \theta_1$ , $x_2=x+l_1 \sin \theta_1 + \dfrac {l_2}{2} \sin \theta_2$ , and $y_1=\dfrac {l_1}{2} \cos\theta_1$ , $y_2=l_1\cos\theta_1+ \dfrac {l_2}{2} \cos\theta_2$ . Therefore, the corresponding velocities are $\dot {x}_1=\dot {x} + \dfrac {l_1}{2} \dot {\theta}_1 \cos \theta_1$ , $\dot {x}_2=\dot {x}+l_1\dot {\theta}_1 \cos\theta_1+\dfrac {l_2}{2}\dot {\theta}_2\cos\theta_2$ , $\dot {y}_1=-\dfrac {l_1}{2}\dot {\theta}_1 \sin\theta_1$ , and $\dot {y}_2=-l_1\dot {\theta}_1\sin\theta_1-\dfrac {l_2}{2}\dot {\theta}_2\sin\theta_2$ . Assuming small angles, or $\sin\theta_i\cong0$ and $\cos\theta_i\cong1$ , we get $\dot {y}_1= \dot {y}_2\cong0$ and

$\begin{equation*} \begin{dcases} \dot {x}_1=\dot {x}+\dfrac {l_1}{2}\dot {\theta}_1 \\ \dot {x}_2=\dot {x}+l_1\dot {\theta}_1+\dfrac {l_2}{2}\dot {\theta}_2 \end{dcases} \end{equation*}$

Figure 9-7 An inverted double pendulum system

The following video shows how to build and run the model for this example in 20-sim.

One or more interactive elements has been excluded from this version of the text. You can view them online here: https://pressbooks.bccampus.ca/engineeringsystems/?p=72#oembed-8

Figure 9‑8 shows the BG model for this system. For further reading, consult with the references cited as [20] and [32].

Figure 9‑8 BG model for the inverted double pendulum

Exercise Problems for Chapter 9

Exercises

Use the BG model given in section 9.2 to complete the following:
1. Build the BG model.
2. Reset the causalities and manually apply them to identify the algebraic loop and related power bonds. List all the options that might exist.
3. Using some typical numerical values for the car-brake system parameters, build the model and run some simulation scenarios.
Use the BG model given in section 9.3 to complete the following:
1. Build the BG model.
2. Reset the causalities and manually apply them to related power bonds.
3. Draw the arrows for showing the streams of flow and effort in the whole system.
4. Using some typical numerical values (as shown below) for the hoist system parameters, build the model and run some simulation scenarios. Make graphs of mass velocity and study the effects of string elasticity and drum moment of inertia using the Parameter Sweep tool in 20-sim.

Frequency Analysis: Bode Plots and Transfer Function

10.1 Overview

Studying the behaviour of systems with respect to time is the primary approach for systems modelling and analysis. However, when a system behaves in a repetitive mode under some applied load and/or boundary conditions—a quasi-static mode—we are interested in the changes in terms of the inherently-involved frequencies of the system rather than the details of instantaneous variations with respect to time. Therefore, transforming the domain of analysis from time to frequency provides us with useful and important information about the behaviour of systems. For example, identifying characteristics of a system—such as its natural frequency, behaviour at large and small frequencies, and magnitude of certain quantities at specific frequencies—provides useful insights in terms of system analysis, design, and control.

In this section, we present a brief background of frequency analysis and methods with focus on Bode plot method and transfer function, with worked-out examples. However, this textbook does not present a full discussion of control theory and related methods for system analysis.

For further reading, consult with references cited in this chapter. 20-sim has tools for performing frequency analysis using BG models and for drawing Bode plots for systems. Through some examples, we will demonstrate how to use these tools.

10.2 Background

Pierre-Simon Laplace (1749 – 1827)

Performing analysis in frequency domain requires a transformation from time domain to frequency domain. Having a mathematical model describing the behaviour of a system, we can use a transformer to convert the governing equations from time domain to frequency domain $\omega$ . Using applied engineering mathematics, we usually employ Laplace and/or Fourier transforms for such an operation. The Laplace transform (defined in complex $s$ -domain) is a more general case of Fourier’s transform (defined in $\omega$ -domain), as given below for transforming a function of time $F(t)$ , [12], [35].

(10.1) $\begin{equation*} \begin{dcases} \tilde{F}(s)=\mathcal{L}(F(t))={\int_{0}^{t}}e^{-st}F(t)dt\:\textit{,Laplace}\\ \check{F}(\omega)=\mathcal{F}(F(t))= {\int_{-\infty}^{+\infty}}e^{-j\omega{t}}F(t){dt}\:,\textit{Fourier}\\ \end{dcases} \end{equation*}$

where $j=\sqrt{-1}$ and for $s=j\omega$ these two transforms are comparable. Note that, in principle, $s=\sigma+{j\omega}$ but the real part, $\sigma$ of complex variable $s$ , is not included here since we are interested in equilibrium at a steady state in frequency analysis. Application of these transforms greatly simplifies the solutions of system equations, both for ODEs and PDEs. The original functions in time domain can be calculated back using the inverse transforms of Laplace and Fourier, given as

Jean Baptiste Joseph Fourier (1768–1830)

(10.2) $\begin{equation*} \begin{dcases} \textit{F}(t)=\mathcal{L}^{-1}({\tilde{F}(s))=\dfrac{1}{2\pi{j}}{\int_{-j\omega}^{j\omega}}e^{st}\tilde{F}(s)ds\textit{, inverse Laplace}\\ \textit{F}(t)=\mathcal{F}^{-1}(\check{F}(\omega))=\dfrac{1}{2\pi} {\int_{-\infty}^{+\infty}}e^{j\omega{t}}\check{F}(\omega){d}\omega\:,\textit{inverse Fourier}\\ \end{dcases} \end{equation*}$

In control theory for systems, several methods are used for studying system behaviour and design, including Bode and Nichols plots requiring frequency-domain response representation, root-locus method requiring complex-domain pole-zero representation, and polynomial-domain design requiring transfer-function representation [31], [36]. A transfer function, by definition, is the ratio of the output signal (i.e., magnitude, power) of a system over selected input signal values. Among these methods, we focus on Bode plots for application in system design using BG models. Bode plots are used for linear systems or linearized non-linear systems. For more detail, see Reference Manual 20-sim 4.6. To access the manual, from the 20-sim Editor window, go to Help, and then select Manual (PDF).

10.3 Motivational Example: A Linear System

A linear system, or LTI (linear time invariant), is a system for which the linear combination of output of a set of inputs is equal to the sum of outputs resulting from each input, e.g., as shown in Figure 10‑1.

Figure 10-1 Linear system sketch for processing inputs and outputs

Considering a proportional integral (PI) controller (for which the output is equal to the input multiplied by a constant and added to the integral of the input), we provide a sinusoidal signal input with frequency of 0.5 rad/s (or $\dfrac{0.5}{2\pi}Hz$ ) given as 0.5 $t$ . Assuming the constant multiplier to be 2, we get the output as $F(t)=2\sin0.5t+\int\sin0.5t\:{dt}=2(\sin0.5t-\cos0.5t$ ). From the output, we observe that the frequency remains the same as that of the input. Now, by multiplying $\dfrac{\sqrt{2}}{\sqrt{2}}=1$ to the output, we rewrite it as a single sinusoidal function, or:

$F(t)=2\sqrt{2}(\dfrac{1}{\sqrt{2}}\sin0.5t-\dfrac{1}{\sqrt{2}}\cos0.5t)=2\sqrt{2}(\cos\dfrac{\pi}{4}\sin0.5t-\sin\dfrac{\pi}{4}\cos0.5t)=2\sqrt{2}\sin(0.5t-\dfrac{\pi}{4}).$

Therefore, the input amplitude is magnified by a factor of $2\sqrt{2}$ , frequency remains the same as mentioned, and a phase change of $-\dfrac{\pi}{4}=-45^{\circ}$ is introduced to the output signal by the PI controller. But one can ask the question: what would be the controller response to a similar signal with a different frequency? For example, if we repeat the same calculation for an input signal given as $\sin0.8t$ , we get the response $\dfrac{\sqrt{89}}{4}\sin(0.8t-32.01^{\circ})$ . To find the response to a spectrum of input frequencies, in principle we can repeat similar calculations and analyze the system behaviour. Table 10‑1 shows some typical response result for input signals of the form $\sin\omega{t}$ into a PI controller.

Table 10‑1 Response of a PI controller to a sinusoidal signal, $\sin\omega{t}$
Frequency (rad/s), $\omega$	Response Amplitude	Response Phase (deg.)
0.1	10.198	-78.690
0.2	5.385	-68.199
0.3	3.887	-59.036
0.4	3.202	-51.340
0.5	2.828	-45.000
0.6	2.603	-39.806
0.7	2.458	-35.538
0.8	2.358	-32.005
0.9	2.288	-29.055
1	2.236	-26.565
2	2.062	-14.036
3	2.028	-9.462
4	2.016	-7.125
5	2.010	-5.711
6	2.007	-4.764
7	2.005	-4.086
8	2.004	-3.576
9	2.003	-3.180
10	2.002	-2.862

One can make graphs of the response amplitude and phase changes versus input frequency to study the behaviour of the PI controller used in this example.

10.4 Bode Plots and Cutoff Frequency

Hendrik Wade Bode (1905–1982)

The results of the calculations mentioned in section 10.3 motivate us to look for a more general and practical method of frequency analysis of systems. The calculations presented are laborious, although new computer tools, e.g., Excel or even computer coding can speed up the process. However, Hendrik Wade Bode, working at Bell Labs in the 1930s, suggested a more practical and now commonly used graphical method—Bode plots, [21], [31], [37], [38]. This method has a wide range of applications in system dynamics, control, and design. An important outcome of having Bode plots for a system is a quick visual insight into the system’s dynamical behaviour for a wide range of frequencies.

In this section, we present the basic idea and some formulas related to Bode plots and what they intend to represent when applied to a system.

Following the example presented in section 10.3, we assume, without losing generality, to have an input sinusoidal signal to the PI controller with frequency $\omega$ given as $A_0\sin\omega{t}$ . The output is, then, $A\sin(\omega{t}+\varphi)$ with amplitude $A$ and phase $\varphi$ . The relations for the output amplitude and phase angle depend on the PI controller specifications. For our example having a proportionality factor of two and an integration, using similar manipulations as those given in section 10.3, we find $A=A_0\sqrt{4+1/\omega^{2}$ and $\varphi=\tan^{-1}(-\dfrac{1}{2\omega})$ . Now, after transforming the input and output signals to the Laplace $s$ -domain, using Laplace transform (see Equation (10.1)), we get $\mathcal{L}\underbrace{(A_0\sin\omega{t})}_{input}=A_0\dfrac{\omega}{s^{2}+\omega^{2}}$ and $\mathcal{L}\underbrace{(2A_0\sin\omega{t}+{\smallint}A_0\sin\omega{t})}_{output}=A_0(\dfrac{2\omega}{s^{2}+\omega^{2}}+\dfrac{\omega/s}{s^{2}+\omega^{2}})$ . The transfer function $G(s)$ (also referred to as gain function) is defined as the ratio of the output amplitude over input amplitude. For our examples, we get, after some simplifications, $G(s)=2+\dfrac{1}{s}$ . Having the transfer function, we can substitute $s=j\omega$ , to transform from $s$ -domain to $\omega$ -domain, guided by Equation (10.1). Therefore, $G(\omega)=G(s)|_{s=j\omega}=2+\dfrac{1}{j\omega}$ , or the transfer function in frequency domain. This is a function of a complex variable and can be written, in general, as its real and imaginary parts or $G(\omega)=\mathfrak{R}+j\mathfrak{J}$ . Therefore, after some manipulations, we get $G(\omega)=2+\dfrac{1}{j\omega}=2-\dfrac{j}{\omega}$ . This gives the real part $\mathfrak{R}=2$ and imaginary part $\mathfrak{J}=-1/\omega$ . Therefore, the magnitude of the transfer function is $|{G}|=\sqrt{4+1/\omega^{2}}$ . The relations for transfer function including its magnitude and phase are summarized in Equation (10.3).

(10.3) $\begin{equation*} \begin{dcases} \textit{G}(\omega)=\mathfrak{R}+j\mathfrak{J}\\ |\textit{G}(\omega)|=\sqrt{\mathfrak{R}^{2}+\mathfrak{J}^{2}}\\ \varphi=\tan^{-1}(\frac{\mathfrak{J}}{\mathfrak{R}})\\ \end{dcases} \end{equation*}$

A set of plots consisting of magnitude $|G(\omega)|$ and $\varphi$ versus logarithm (at base 10) of $\omega$ is called Bode plots. However, magnitude is traditionally measured in decibel (dB), phase angle in degrees, and frequency as logarithm of frequency in rad/s (or in Hz).

Recall that dB (one tenth of a bel) is a unit for measuring the power of a signal with reference to a threshold. For example, the threshold for human hearing is $10^{-12}W/m^{2}$ , given as power intensity; dB is measured in logarithm of the power ratios at base ten, or $dB=10\log_{10}\dfrac{power\:out}{power\:in}$ . But since power of a wave signal is proportional to its amplitude squared, then we get $dB=20\log_{10}\dfrac{amplitude\:out}{amplitude\:in}$ , or

(10.4) $\begin{equation*} \textit{d}B=20\log_{10}|G(\omega)|\\ \end{dcases} \end{equation*}$

Therefore, Bode plots are composed of two charts: signal gain in dB and phase in degrees versus logarithm of $\omega$ , usually given in a single graph chart. The graphs in Figure 10‑4 show the Bode plots for the PI controller, generated using 20-sim. The tools available for drawing Bode plots in 20-sim can also be used when a transfer function is available or calculated and also after a BG model is built for systems. For this example, we calculated the transfer function and used it to draw the corresponding Bode plots. For drawing the corresponding Bode plots, follow these steps:

From the Editor window, go to Tools, select Frequency Domain Toolbox, and then select Linear System Editor. The 20-sim Linear System Editor window opens, as shown in Figure 10‑2.

Figure 10‑2 Linear System Editor interface in 20-sim

Select Transfer Function and click on the Edit A window opens in which you can enter the coefficients of the transfer function, as shown in Figure 10‑3.
Enter the desired transfer function as a polynomial fraction with numerator and denominator polynomials with their corresponding coefficients in descending power of $s$ . For this example, having $G(s)=2+\dfrac{1}{s}=\dfrac{2{s}+1}{s},$ the coefficients for the numerator polynomial are (2, 1) and for the denominator polynomial are (1, 0). Note the zero term, i.e., the coefficient for the term $s^{0}$ , or the constant term. A space can be used instead of a comma, to separate the coefficients.
Click on Apply and then OK This takes you back to the Linear System Editor window with the transfer function listed. Double check the resulted transfer function to make sure it is entered correctly into 20-sim.

Figure 10‑3 Transfer Function Editor interface in 20-sim
From the 20-sim System Editor window (see Figure 10‑2) under Plots, select Bode. The Bode Plot window opens with the corresponding Bode plots, as shown in Figure 10‑4. From these plots, we can conclude that the PI controller is a high-pass filter system because it passes through the high frequency signals but attenuates low frequency signals. The gain for low frequency signals decreases linearly with a slope of 20 dB per decade from 40 dB. The phase change is from close to $-90^{\circ}$ at low frequencies to null at high frequencies. Therefore, the PI controller is in phase with the input signals at high frequencies and out of phase, by about $-90^{\circ}$ , at low frequencies.
The asymptotes to the high frequency and low frequency gains intersect at a point defined as the cutoff frequency, $\omega_{b}$ (also referred to as corner or break frequency). This frequency is defined when the output power reaches to 50% of the input signal power (so-called half-power point), or $|dB|=\Big|10\log_{10}\dfrac{power\:out}{power\:in}\Big|=|10\log_{10}0.5|=\Big|20\log_{10}\underbrace{0.5^{0.5}}_{=1/\sqrt{2}}\Big|\cong3$ .
The 3dB-point is the standard method of finding cut off frequency from the Bode plot gain chart. For this example, the asymptote to the high frequency gain (i.e., the horizontal line as the frequency $\rightarrow\infty$ ) is at about 6.02 dB. Hence, the cutoff frequency corresponds to the point at 6.02+3=9.02 dB, or $|G(\omega)|=2.825$ using Equation (10.4). This gives the cutoff frequency of $\omega_{b}=0.5013$ rad/s, using $|G(\omega)|=2.825=\sqrt{4+1/\omega^{2}$ . At the cutoff frequency the phase reads $-45$ deg, using $\varphi=\tan^{-1}(-\dfrac{1}{2\times0.5013})$ , after conversion.

Figure 10-4 Bode plots for a PI controller

By following the steps presented in the following section, we can also draw Bode plots using 20-sim when a BG model of a system is available.

10.4.1 Guideline for Drawing Bode Plots for BG Models

After using 20-sim to build a BG model, to draw related Bode plots, click on Tools and select the options in the Frequency Domain Toolbox to draw related Bode plots. We can draw Bode plots using transfer function, either manually or using computer graphing tools. The following steps can be used for drawing Bode plots, by either (A) using 20-sim software tool or (B) manually:

1. Drawing Bode plots using 20-sim (See the 20-sim Reference Manual.)
  1. Build BG model. Include data for related variables.
  2. From the Editor window, go to Tools, select Frequency Domain Toolbox, and then select Model Linearization. The 20-sim Model Linearization window opens, as shown in Figure 10‑5.

Figure 10‑5 Model Linearization interface in 20-sim

1. Select the input and output variables for which you want to draw the Bode plots. The resulting transfer function is related to this selection. Click on the Variable Chooser icon to get a list of model variables to choose from. Leave the rest of options as selected by default. Note that unless output is used as feedback, usually Open Loop is selected. Select OK.
2. The 20-sim System Editor interface window opens (see Figure 10‑2). The transfer function, based on input/output selection appears. From the selections under Plots, choose the Bode The corresponding Bode plots appear in a new window as shown in Figure 10‑6.

Figure 10‑6 Typical Bode plots for a system

1. 1. Go to Properties and select Plots to edit plots for title, axes scales, labels, legends, etc.

1. Manual drawing of Bode plots using transfer function

1. Derive transfer function and transform it to $s$ -domain, $G(s)$ , using Laplace transform.
2. Plug in $s=j\omega$ into transfer function, to get $G(j\omega)$ .
3. Calculate the real and imaginary parts of the $G(j\omega)$ .
4. Calculate magnitude $|G(j\omega)|$ and power, using Equation (10.4).
5. Calculate the phase angle in degrees, using Equation (10.3).

Alternatively, Bode plots can be drawn using 20-sim after having the desired transfer function $G(s)$ from step 1 above, by following the guideline given in section 10.4.

10.5 Example: Bode Plots Using Transfer Function

A system’s transfer function $G(s)=\dfrac{1}{s^{2}+6s+8}=\dfrac{1}{(s+2)(s+4)}$ is given. Find expressions for magnitude and phase angle and draw the corresponding Bode plots. Consider the frequency range 0.01–1000 rad/s. Discuss the system dynamical behaviour at low and high frequencies, including cutoff frequency.

Solution:

Using the guidelines, we substitute $s=j\omega$ in $G(s)$ . Therefore, $G(j\omega)=\dfrac{1}{(j\omega+2)(j\omega+4)}=\dfrac{1}{8-\omega^{2}+6j\omega}$ After multiplying and dividing by the conjugate of the denominator, we get $G(j\omega)=\dfrac{1}{8-\omega^{2}+6j\omega}\times\dfrac{8-\omega^{2}-6j\omega}{8-\omega^{2}-6j\omega}=\dfrac{8-\omega^{2}}{\omega^{4}+20\omega^{2}+64}-j\dfrac{6\omega}{\omega^{4}+20\omega^{2}+64}$ . From this expression, we get the real and imaginary parts as $\mathfrak{R}=\dfrac{8-\omega^{2}}{\omega^{4}+20\omega^{2}+64}$ and $\mathfrak{J}=-\dfrac{6\omega}{\omega^{4}+20\omega^{2}+64}$ . Using the real and imaginary parts and Equation (10.3), we can calculate the magnitude and phase as $|G(j\omega)|=\dfrac{1}{\sqrt{\omega^{4}+20\omega^{2}+64}}$ and $\varphi=\tan^{-1}(\frac{6\omega}{\omega^{2}-8})$ . Note that at $\omega=2\sqrt{2},\varphi=\pi/2$ and at $\omega=0,\varphi=\pi$ . We can draw the Bode plots, manually or using 20-sim, for example. The cutoff frequency can be calculated as follows. The gain magnitude at $\omega=0$ reads $|G(j\omega)|=\frac{1}{\sqrt{64}}=0.125$ . Therefore, the corresponding power is $20\log_{10}0.125=-18.0618\: dB$ . The cutoff frequency corresponds to $-18.0618+(-3)=-21.0618\:dB$ . Therefore, the corresponding gain is $G=0.0885$ , using $-21.0618=20\log{G}$ . After substituting for $G$ and using $G=\dfrac{1}{\sqrt{\omega^{4}+20\omega^{2}+64}}$ , we get cutoff frequency $\omega_{b}=1.672\:rad/s$ . The phase at the cutoff frequency can be calculated using $\varphi=\tan^{-1}(\dfrac{6\times1.672}{1.672^{2}-8})=-1.092\:rad$ or $\varphi=-62.57^{\circ}$ .

Following the steps given in section 10.4, we can use 20-sim and the transfer function $G(s)=\dfrac{1}{s^{2}+6s+8}$ to draw the Bode plots. The coefficients for the polynomials are (1) for numerator, and (1, 6, 8) for denominator. The resulting Bode plots are shown in Figure 10‑7.

The following video shows how to build and run the model for this example in 20-sim.

One or more interactive elements has been excluded from this version of the text. You can view them online here: https://pressbooks.bccampus.ca/engineeringsystems/?p=74#oembed-3

Figure 10‑7 Bode plots for the transfer function

10.6 Example: Bode Plots Using a BG Model

Build the BG model and its Bode plots considering transfer function based on force $F$ as input and displacement of the spring $k_{2}$ as output. Use the mechanical system as shown in Figure 10‑8 with

$m_{1}=20\:kg$ , $m_{2}=10\:kg$ , $k_{1}=35\:kN/m$ , $k_{2}=38kN/m$ , $k_{3}=15kN/m$ , $b_{1}=0.4\:N.s/m$ , $b_{2}=0.5\:N.s/m,$ $F=200\:N$ , and $b_{3}=0.1\:N.s/m$ .

The damper $b_{3}$ connects mass $m_{1}$ and the wall at the right side. Repeat the simulation for the following cases:

1. 1. Use given damper coefficient values $b_{3}=0.1,0.2,0.3,0.4,0.5,0.6 N.s/m$ to study its effect on the system with Parameter Sweep tool in 20-sim (available at the 20-sim Simulator window: select Tools > Time Domain Toolbox >Parameter Sweep). During the sweep, monitor the displacement of mass $m_{1}$ .
  2. Use a pulse-type signal as the applied force with amplitude 200 $N$ , start time 2 sec, and stop time 3.5 sec.

Figure 10-8 Mechanical system sketch for example given in section 10-6

Solution:

The following video shows how to build and run the model for this example in 20-sim.

One or more interactive elements has been excluded from this version of the text. You can view them online here: https://pressbooks.bccampus.ca/engineeringsystems/?p=74#oembed-4

Figure 10‑9 shows the resulting Bode plots. Note that in this video, the force is applied to mass $m_{2}$ , in the first try and then moved to mass $m_{1}$ according to the sketch shown in Figure 10‑8.

Figure 10-9 Bode plots for mechanical system given in section 10-6

Media Attributions

Hendrik Wade Bode is licensed under a Public Domain license

Miscellaneous Topics

11.1 Overview

This chapter discusses several topics pertinent to modelling and simulation of systems dynamics with focus on BG method. The reader can consult with these sections as independent topics or as supplementary to previously discussed topics.

11.2 Energy and Power Conjugate Variables

The product of conjugate variables related to a physical quantity should, by definition, give the corresponding quantity. In chapter 3, sections 3.2 and 3.5, we discussed the power variables (effort $e$ and flow $f$ ) and their relations to the state variables (momentum $p$ and displacement $q$ ). The product ( $e.f$ ) is power—so-called power variables—and integrating them, with respect to time, gives $p$ and $q$ , respectively. Hence the $p$ and $q$ are called energy variables. However, further attention to the dimension of the product of the variables, i.e., ( $p.q$ ), indicates that the dimension of this quantity is not equivalent to that of energy. For example, in the SI system of units, we get $[p.q] = [N.s.m] = [J.s]$ . Therefore, we can call $e$ and $f$ conjugate power variables because the product ( $e.f$ ) is power. But this definition does not apply to the product ( $p.q$ ). Instead, the product ( $\dot p . q$ ) or product ( $p. \dot q$ ) can be defined as true conjugate energy variables.

One can ask the question, what is the relationship among these variables? To investigate further and provide a possible answer, we consider the system energy function $\mathcal{E} = \mathcal{E} (p,q)$ . This functional form is legitimate, since we know that the energy stored, in storage elements $I$ and $C$ , is uniquely defined by state variables (see section 3.5). The total differential change of energy is then $d \mathcal{E} = \dfrac{\partial \mathcal{E}}{\partial p} dp + \dfrac{\partial \mathcal{E}}{\partial q} dq$ . In principle, $\mathcal{E} = \mathcal{E}(p,q,t)$ . But the partial derivative $\dfrac{\partial \mathcal{E}}{\partial t} |_{p,q}$ can be ignored at equilibrium state of the system, or when $p$ and $q$ are kept constant and energy is not changing with time. Now, the total differential change of power $d \mathcal{P}$ can be written as the time derivative of energy, or $d \mathcal{P} = \dfrac{d \mathcal{E}}{dt} = \dfrac{\partial \mathcal{E}}{\partial p} \dfrac{dp}{dt} + \dfrac{\partial \mathcal{E}}{\partial q} \dfrac{dq}{dt}$ . But $\dfrac{dp}{dt} = e$ and $\dfrac{dq}{dt} = f$ , the effort and flow in BG method, respectively. Therefore, we can write

(11.1) $\begin{equation*} d \mathcal{P} = e \dfrac{\partial \mathcal{E}}{\partial p} + f \dfrac{\partial \mathcal{E}}{\partial q} \end{equation*}$

Equation (11.1) clearly defines the relationship among power variables (effort $e$ and flow $f$ ) and state variables (momentum $p$ and displacement $q$ ) with respect to the power ( $\mathcal{P}$ ) and energy ( $\mathcal{E}$ ) of the system. Note that for a system, we can write this equation for each component and sum it up to get the relation for the system.

From Equation (11.1), we can conclude that by dimensional analysis, $\dfrac{\partial \mathcal{E}}{\partial p}$ has the dimension of flow and $\dfrac{\partial \mathcal{E}}{\partial q}$ the dimension of effort, according to BG terminology. We define these quantities as system flow and effort, or (note that the partial derivative definition is explicitly shown here)

(11.2) $\begin{equation*} \begin{dcases} \dfrac{\partial \mathcal{E}}{\partial p} |_{q=constant} , \: system \: flow \\ \dfrac{\partial \mathcal{E}}{\partial q} |_{p=constant} , \: system \: effort \end{dcases} \end{equation*}$

For example, considering a mass, represented by an $I$ -element, its kinetic energy can be written as $\dfrac{p^2}{2m}$ , but $\dfrac{\partial}{\partial p} \left( \dfrac{p^2}{2m} \right) = \dfrac{p}{m} = \dfrac{mv}{m} = v$ , the velocity of the mass or the flow associated with the $I$ -element. Similarly, considering a mechanical spring with stiffness $k$ , represented by a $C$ -element, its potential energy can be written as $k \frac{q^2}{2}$ , but $\dfrac{\partial}{\partial q} \left( k \dfrac{q^2}{2} \right) = kq = kx = F$ , the force acting on the spring or the effort associated with the $C$ -element. Equation (11.2) applies analogously to mechanical, electrical, hydraulic, or thermal systems. Note that the energy function $\mathcal{E}$ relates to the Hamiltonian of the system, [13].

11.3 Including Gravity in BG Models

Forces due to gravitational acceleration $g$ can be implemented in a BG model using effort source element $S_e$ . But it feeds power to the system when the assigned positive direction of displacement/velocity and the gravitational acceleration vector are in the same direction. Otherwise, the system feeds power into the gravity field, or loses power. In manually drawn BG, the direction of the power bond is drawn from $S_e$ toward the system when the gravity and displacement directions are matched; otherwise, the power bond should be toward the $S_e$ element. In 20-sim, however, it is not possible to draw a power bond pointing toward an $S_e$ element because it is defined as a source element. Hence, a negative value of $g$ is assigned to the associated parameter in the corresponding equation model, i.e., −9.81 $m/s^2$ (see Figure 11-1). Note that $g$ is interpreted as gravitational force per unit mass, here.

Figure 11-1 Implementing gravity force for a BG model in 20-sim.

11.4 Extracting System Equations from BG Models

Although graphics of BG models can give us some insights into the dynamics of systems before actually solving the systems equations, having the mathematical model consisting of a system of equations is the ultimate goal resulting from a BG model, or for that matter, any model. In the previous sections, we discussed that for a BG model with integral causalities, we can derive a system of ODEs in terms of state variables, i.e., $p$ and $q$ (see section 3.5). For such a system, the number of equations, $n$ is equal to the sum of number of $p$ associated with $I$ – elements and number of $q$ associated with the $C$ – elements. In other words, the system equations are coupled first-order ODEs explicitly given in terms of state variables and can be solved simultaneously using numerical schemes. The solution of these coupled equations is, obviously, easier to obtain for linear systems rather than it is for non-linear systems.

In this section we discuss the procedure and guidelines of how to extract system equations from a BG model. For now, we assume that all storage elements (i.e., $I$ – and $C$ – elements) in the system are assigned with integral causality strokes. In the next section, we discuss cases when derivative causality and/or algebraic-loop situations may exist in a BG model, along with the definitions and resulting consequences for such systems equations.

Referring to the energy management of a system, we have already discussed that the energy input to a system is partially stored in storages and the rest is dissipated through dissipaters. Therefore, the storage elements are key elements for studying the dynamics of any system. This principle also applies to the BG models, and we use it to extract system equations.

The system equations can be given in a general form as $\dot x_i = \mathcal{F}_i (x_i, u_j, t)$ where $\dot x_i$ , the vector representing time derivative of state variables, $u_j$ , the vector representing the inputs, with indices $i = 1, \cdots , n,$ and $j = 1, \cdots, r$ . One can expand this functional to get

(11.3) $\begin{equation*} \begin{dcases} \dot x_1 = \mathcal{F}_1 (x_1, x_2, \dots, x_n, u_1, u_2, \dots, u_r, t) \\ \dot x_2 = \mathcal{F}_2 (x_1, x_2, \dots, x_n, u_1, u_2, \dots, u_r, t) \\ \quad \quad \quad \quad \quad \quad \quad \quad \vdots \\ \dot x_n = \mathcal{F}_n (x_1, x_2, \dots, x_n, u_1, u_2, \dots, u_r, t) \\ \end{dcases} \end{equation*}$

After building the BG model for a system, we ask the following two key questions, [21]:

Q1: What does each component/element send to the system?

Q2: What does the system send back to the storage components/elements?

Each question guides us to list the corresponding relations involving the state variables and performing some manipulations to write the system equations in terms of the state variables, using the general guideline of $p$ on $I$ and $q$ on $C$ . The following steps can be used to extract the system equations:

Build the BG model and assign causalities, with preferred integral causalities.
Simplify the BG model, as far as possible.
Assign labels as numbers to each power bond in the BG model. The order of numbering is arbitrary.
Answer Q1, in terms of all efforts or flows input into the system.
Answer Q2, in terms of all $p$ ’s on $I$ ’s and all $q$ ’s on $C$ ’s, or momenta (or their derivatives) associated with the inertia elements and displacements (or their derivatives) associated with the capacity elements.
Apply constraints implemented by 1- and 0- junction elements and perform algebraic manipulations to write the system equations only in terms of state variables, as independent variables with inputs and components’ properties as parameters.

Now, we present some worked-out examples to demonstrate the application of the above mentioned procedure.

11.4.1 Example: System Equations for a Mechanical System

For this example, we consider a mechanical system as shown in Figure 11-2. This systems has three DOF associated with three masses. The number of state variables is seven, associated with the storage elements, three masses, and four springs. The positive displacements are considered as shown, and tension forces are positive (+T).

Figure 11-2 A mechanical system sketch

Solution:

We first build the BG model, either manually or using the tools available through 20-sim, with all power bonds labelled, as shown in Figure 11-3. As mentioned, the labels are arbitrary; they merely help bookkeeping of the variables involved without affecting final solution of the system equations. Therefore, the state variables are: $p_2, p_6, p_{13}$ associated with the momenta of the masses and $q_4, q_{10}, q_{15}, q_{19}$ displacements of the mechanical springs. The corresponding power bonds to these elements are colour coded (see Figure 11-3). The system equations are ODEs consisting of these state variables. To extract the system equations, we ask two key questions (see section 11.4) and list their corresponding answers, using the power bond labels and causality assignments in terms of BG symbols $e$ and $f$ , as follows:

Figure 11-3 BG model for the mechanical system with labelled power bonds

Q1: What does each component/element send to the system?

The inputs to the system are from elements at the boundary of the system associated with bonds 1, 2, 4, 6, 10, 11, 13, 15, 17, and 19, as listed in Table 11-1.

Table 11-1 Inputs to the mechanical system by all elements—Q1 question
Power bond label #	System input/relation
1	$e_1 = F$
2	$f_2 = \dot q_2 = p_2/m_3$
4	$e_4 = k_4q_4$
6	$f_6 = \dpt q_6 = p_6/m_2$
10	$e_{10} = k_3q_{10}$
11	$e_{11} = b_1f_{11}$
13	$f_{13} = \dot q_{13} = p_{13}/m_1$
15	$e_{15} = k_2q_{15}$
17	$f_{17} = V_{wall}$
19	$e_{19} = k_1q_{19}$

Next, we list the relations to answer the second question:

Q2: What does the system send back to the storage components/elements?

Here, we are only interested in storage elements, as shown by colour-coded bonds in Figure 11-3. For example, considering power bond number 2, the system sends an effort signal (or rate of momentum) to the $I$ -element representing the mass $m_3$ . Therefore, we can write $e_2 = \dot p_2 = m_3 \dot f_2$ . Similarly, considering power bond number 15, the system sends a flow signal (or rate of displacement) to the $C$ -element representing the spring $k_2$ . Therefore, we can write $f_{15} = \dot q_{15}$ . Table 11-2 lists all relations corresponding to storage elements.

Table 11-2 Outputs from the mechanical system to all storage elements—Q2 question
Power bond label #(element)	Storage element/relation	Relations, using constraints
2 ( $I:m_3$ )	$e_2 = \dot p_2$	$e_2 = \pmb{\dot p_2} = e_1 - e_3 = e_1 - e_4 = \pmb{F - k_4q_4}$
4 ( $C:k_4$ )	$f_4 = \dot q_4$	$f_4 = \pmb{\dot q_4} = f_3 - f_5 = f_2 - f_6 = \pmb{\dfrac{p_2}{m_3} - \dfrac{p_6}{m_2}}$
6 ( $I:m_2$ )	$e_6 = \dot p_6$	$e_6 = \pmb{\dot p_6} = e_5 - e_7 - e_8 = \pmb{k_4q_4 - k_3 q_{10}} - b_1 f_{11} - \pmb{k_1 q_{19}}$
10 ( $C:k_3$ )	$f_{10} = \dot q_{10}$	$f_{10} = \pmb{\dot q_{10}} = f_9 = f_7 - f_{12} = f_6 - f_{13} = \pmb{\dfrac{p_6}{m_2} - \dfrac{p_{13}}{m_1}}$
13 ( $I:m_1$ )	$e_{13} = \dot p_{13}$	$e_{13} = \pmb{\dot p_{13}} = e_{12} - e_{14} = e_9 - e_{15} = (e_{10} + e_{11}) - e_{15} = \pmb{k_3q_{10}} + \underbrace{b_1f_{11}}_{=b_1f_{10}} - \pmb{k_2q_{15} = k_3q_{10}} + b_1 \pmb{ \left( \dfrac{p_6}{m_2} - \dfrac{p_{13}}{m_1} \right) - k_2q_{15}}$
15 ( $C:k_2$ )	$f_{15} = \dot q_{15}$	$f_{15} = \pmb{\dot q_{15}} = f_{14} - f_{16} = f_{13} - f_{17} = \pmb{\dfrac{p_{13}}{m_1} - V_{wall}}$
19 ( $C:k_1$ )	$f_{19} = \dot q_{19}$	$f_{19} = \pmb{\dot q_{19}} = f_8 - f_{18} = f_6 - f_{17} = \pmb{\dfrac{p_6}{m_2} - V_{wall}}$

In the relations listed in the third column of Table 11-2, we used the constraints resulted from 1- and 0-junctions and elements’ constitutive equations, and included the relations from Table 11-1. The final relations are given in bold. Finally, after some manipulations, we can write the system equations in terms of seven state variables in matrix form, as given in Equation (11.4).

(11.4) $\begin{equation*} \begin{matrix} \begin{Bmatrix} \dot p_2\\ \dot q_4 \\ \dot p_6 \\ \dot q_{10} \\ \dot p_{13} \\ \dot q_{15} \\ \dot q_{19} \end{Bmatrix} = \begin{bmatrix} 0 & -k_4 & 0 & 0 & 0 & 0 & 0 \\ 1/m_3 & 0 & -1/m_2 & 0 & 0 & 0 & 0 \\ 0 & k_4 & -b_1/m_2 & -k_3 & b_1/m_1 & 0 & -k_1 \\ 0 & 0 & 1/m_2 & 0 & -1/m_1 & 0 & 0 \\ 0 & 0 & b_1/m_2 & k_3 & -b_1/m_1 & -k_2 & 0 \\ 0 & 0 & 0 & 0 & 1/m_1 & 0 & 0 \\ 0 & 0 & 1/m_2 & 0 & 0 & 0 & 0 \end{bmatrix} \\ \\ \begin{Bmatrix} p_2\\ q_4 \\ p_6 \\ q_{10} \\ p_{13} \\ q_{15} \\ q_{19} \end{Bmatrix} + \begin{Bmatrix} F\\ 0 \\ 0 \\ 0 \\ 0 \\ -V_{wall} \\ -V_{wall} \end{Bmatrix}} \end{matrix} \end{equation*}$

11.4.2 Example: System Equations for an Electrical System

For this example, we consider the electrical system and BG model given in section 7.6 and shown in Figure 7-8. The system has two loops. The number of state variables is three, associated with the storage elements, two inductors and one capacitor.

Solution:

Figure 11-4 shows the BG model with labelled power bonds for this system. The power bonds associated with storage elements, hence state variables, are colour coded.

Figure 11-4 BG model for the electrical system with labelled power bonds

As mentioned, the labels are arbitrary and they merely help bookkeeping of the variables involved without affecting final solution of the system equations. The state variables are: $p_2,p_5,$ associated with the momenta of the inductors, and $q_4$ , electrical charge of the capacitor. The system equations are ODEs consisting of these variables. To extract the system equations, we ask two key questions (see section 11.4) and list their corresponding answers, using the power bond labels and causality assignments in terms of BG symbols e and f, as follows:

Q1: What does each component/element send to the system?

The inputs to the system are from elements at the boundary of the system associated with bonds
1, 2, 3, 4, 5, and 6, as listed in Table 11-3.

Table 11-3 Inputs to the electrical system by all elements—Q1 question
Power bond label	System input/relation
1	$e_1 = V$
2	$f_2 = \dot q_2 = p_2/L_2$
3	$e_3 = R_2f_3$
4	$e_4 = q_4/c$
5	$f_5 = \dot q_5 = p_5/L_1$
6	$e_6 = R_1f_6$

Next, we list the relations to answer the second question:

Q2: What does the system send back to the storage components/elements?

Here, we are only interested in storage elements, as shown by colour-coded bonds in Figure 11-4. For example, considering power bond number 2, the system sends an effort signal (or rate of momentum/flux linkage) to the $I$ -element representing the inductor $L_2$ . Therefore, we can write $e_2 = \dot p_2 = L_2 \dot f_2$ . Similarly, considering power bond number 4, the system sends a flow signal (i.e., current or rate of charge) to the $C$ -element representing the capacitor $c$ . Therefore, we can write $f_4 = \dot q_4$ . Table 11-4 lists all relations corresponding to storage elements.

Table 11-4 Outputs from the electrical system to all storage elements—Q2 question
Power bond label (element)	Storage element/relation	Relations, using constraints
2 ( $I:L_2$ )	$e_2 = \dot p_2$	$e_2 = \dot p_2 = e_9 - e_3 = e_7 - e_3 = e_1 - e_6 - e_3 = V - R_1 \underbrace{f_6}_{f_8 + f_9} - R_2 \underbrace{f_3}_{=f_2} = V - R_1 \left( \dfrac{p_5}{L_1} + \dfrac{p_2}{L_2} \right) - R_2 \dfrac{p_2}{L_2}$
4 ( $C:c$ )	$f_4 = \dot q_4$	$f_4 = \dot q_4 = f_5 = \dfrac{p_5}{L_1}$
5 ( $I:L_1$ )	$e_5 = \dot p_5$	$e_5 = \dot p_5 = \underbrace{e_8}_{=e_7} - e_4 = e_1 - e_6 - \dfrac{q_4}{c} = V- R_1 \underbrace{f_6}_{f_8+f_9} - \dfrac{q_4}{c} = V - R_1 \left( \dfrac{p_5}{L_1} + \dfrac{p_2}{L_2} \right) - \dfrac{q_4}{c}$

Note that $f_6 = f_7 = f_8 + f_9 = f_5 + f_2 - \dfrac{p_5}{L_1} + \dfrac{p_2}{L_2}$ . In the relations listed in the third column of Table 11-4, we use the constraints resulted from 1- and 0-junctions and elements’ constitutive equations, and included the relations from Table 11-3. The final relations are given in bold. Finally, we can write the system equations, after some manipulations, in terms of three state variables in matrix form, as given in Equation (11.5).

(11.5) $\begin{equation*} \begin{Bmatrix} \dot p_2\\ \dot q_4 \\ \dot p_5 \end{Bmatrix} = \begin{bmatrix} - \left( \dfrac{R_1}{L_2} + \dfrac{R_2}{L_1} \right) & 0 & - \dfrac{R_2}{L_1}\\ 0 & 0 & \dfrac{1}{L_1} \\ - \dfrac{R_1}{L_2} & - \dfrac{1}{c} & - \dfrac{R_1}{L_1} \end{bmatrix} \begin{Bmatrix} p_2 \\ q_4 \\ p_5 \end{Bmatrix} + \begin{Bmatrix} V \\ 0 \\ V \end{Bmatrix} \end{equation*}$

11.5 Derivative Causality and Algebraic Loop: Implicit System Equations

When assigning causality strokes to a BG model, three scenarios may occur [20], [21]. These scenarios are:

all causality strokes are possible to be assigned as preferred integral causalities;
at least one storage element ( $I$ – or $C$ – element) cannot be assigned with integral
causality; instead derivative causality is forced upon the element; or
there is more than one option for having a BG model with assigned integral causalities for the given system, a so-called algebraic loop.

All of these scenarios are legitimate in terms of BG modelling rules, and with red colour-coded causality strokes, 20-sim provides warnings, not errors, for scenarios 2 and 3 as listed above. As mentioned, scenario 1 is desirable and preferred. For this scenario, the system equations can be derived explicitly in terms of state variables as a system of coupled first-order ODEs. The solution of the system of equations provides the answers for the generalized momenta and displacements associated with storage elements. After having the solutions for system equations, we can calculate other desired quantities related to the system.

For the second scenario, the system equations are not independent, i.e., implicit. In other words, the equations related to the storage elements with integral causalities are independent; pertinent state variables can be uniquely calculated by solving these equations, simultaneously. But the equation related to the storage element with derivative causality is not independent; it relates algebraically to the independent equations. Therefore, the related state variable for the element with derivative causality can be calculated using the solution of those of the independent variables. The mathematical dependence of the equation for the derivatively causalled element indicates that the dynamics of the element are related to and defined by the other storage elements in the system. In addition, for the derivative causality case, we may have to force more than one element having non-integral causalities.

For the third scenario, after assigning integral causalities for some of the elements in a definite manner, we face one or more options and have to make choices for one or more elements and arbitrarily assign them with causalities. The resulting system equations are again independent coupled ODEs and can be simultaneously solved to find their solutions. However, during the manipulation, we encounter an implicit equation (or equations) for a particular state variable (or its derivative) which is a function of inputs and itself as well—an algebraic loop. For linear systems/elements, the algebraic-loop situation does not pose a problem; however, for non-linear systems/components, this makes it more difficult to manipulate the equations to find their solutions.

Note that an algebraic loop may occur at several levels and, hence, make it harder to manipulate the equations. This is when, after selecting an assigning causality stroke for an element arbitrarily, further selection(s) should be made to proceed.

For both scenarios, derivative causalities and algebraic loop, the fact that they appear in a BG model is useful information about the system and/or assumptions made for the model, even before solving the resulting equations. In many cases, the modeller can improve the system with additional components, usually of the types of $R$ – element.

As mentioned, 20-sim sends warning messages, but no error messages when derivative causalities and algebraic loop exist. The resulted system equations can be solved by solvers available, and the system simulation and design can proceed as usual.

In further sections, we present some examples to demonstrate the derivative causality and algebraic-loop cases.

11.5.1 Example: BG Model with Derivative Causality

We consider the BG model shown in Figure 11-5. The colour-coded bonds associated with two inertia $I$ -elements are forced, having non-integral causalities or derivative causality. This can be vividly explored by starting assigning causality from $S_f$ -element sending a flow signal to the adjacent 1-junction. Therefore, all other connecting bonds to this 1-junction should have their causality strokes at the ports connecting to this 1-junction. This requirement forces the associated $I$ -element to have a derivative causality, or sending effort signal, instead of receiving it for being integrally causalled, to the 1-junction. Similarly, the $I$ -element $I_1$ is forced to have derivative causality. As a result, the momenta of these two $I$ -elements (i.e. $I$ and $I_1$ ) will depend on the remaining storage elements, i.e., $I$ -element $I_2$ and $C$ -element $c$ , and can be calculated having the corresponding solutions for $I_2$ – and $C$ – elements. For further reading, consult with references cited as [5], [20], [21], and [30].

As mentioned, in 20-sim a warning message appears on the screen when derivative causality exists in the BG model. But the software solver tools take care of this and solve the system equations. This capability is useful for practical applications in design of the systems.

Figure 11-5 A BG model with derivative causality, colour coded

11.5.2 Example: BG Model with Algebraic Loop

We consider the BG model shown in Figure 11-6. Starting with the process of assigning causalities, we realize that the colour-coded bonds can’t be integrally assigned and present options. In other words, to proceed further with assigning causalities to the whole system’s BG model, we require to make a selection arbitrarily. Usually, the selection can be more easily made with $R$ -elements, since it can be assigned with neutral causality. For example, if we select the element $R:R3$ and assign causality stroke to it such that it sends a flow signal to the adjacent 1-junction, then the remaining bonds’ causality strokes can be assigned, as shown colour coded in Figure 11-7. Alternatively, element $R:R2$ or $R:R1$ could be selected as an option. The corresponding BG models and related assigned causalities are shown in Figure 11-8 and Figure 11-9, respectively.
For further reading, consult with references cited as [5], [20], [21], and [30].

Figure 11-6 A BG model with algebraic-loop causality

Figure 11-7 The BG model with removed algebraic loop—selecting R3

Figure 11-8 The BG model with removed algebraic loop—selecting R2

Figure 11-9 The BG model with removed algebraic loop—selecting R1

As mentioned, in 20-sim, a warning message appears on the screen when algebraic loop causality exists in the BG model. But the software solver tools take care of this and solve the system equations. This capability is useful for practical applications in design of the systems.

11.6 Thermal Systems: Pseudo Bond Graph

Thermal systems are unique in terms of BG modelling. They are closely analogous to electrical systems. In practice, to analyze thermal systems, we use this analogy to build models as thermal networks, similar to electrical networks, [39], [40].

Considering heat conduction through a solid slab, e.g., the temperature difference $\Delta T$ (analogous to voltage) and heat flow rate $\dot Q_{th}$ (analogous to electrical current) satisfy $\Delta T = R_{th} \dot Q_{th}$ . The quantity $R_{th} = \dfrac{\Delta L}{k_TA}$ (analogous to electrical resistance) is thermal resistance with dimension [K/W], in $SI$ , where, $\Delta L$ is the length across which the $\Delta T$ applies; $A$ is the cross-sectional area for heat flow, and $k_T$ is the thermal conductivity of the material. Intuitively, this analogy suggests to consider temperature as the effort and the heat flow rate as the flow when building BG models for thermal systems. With this assumption at hand, the product ( $T \dot Q_{th}$ ) would be considered as power variables, or the dimension of this quantity should be as of that for power. However, close examination of the dimension for the product ( $T \dot Q_{th}$ ) gives its dimension as $[(T\dot Q_{th})] = K \cdot W$ , or $Kelvin \times Watt$ . Therefore, according to BG modelling rules, we cannot accept temperature and heat flow rate as conjugate power variables. This discrepancy leads us to call such a BG model a pseudo bond graph (pseudo BG), for which temperature is the effort and heat flow rate is the flow. This is acceptable as far as the thermal system is the only system involved and the modeller is aware of the fact that the “power” variables are not defined fully and correctly in the related thermal pseudo BG model. However, for multi-domain systems, the exchange of power between a thermal sub-system part and other domains of the system becomes problematic.

The solution to this challenge comes from the second law of thermodynamics, which guides us to consider the flow in the BG model as the entropy rate, instead of heat flow rate. Recall that entropy is a state function, and for a reversible system, we have the relationship $\Delta Q_{th} = T \Delta S$ for exchange of heat $\Delta Q_{th}$ and entropy $\Delta S$ using temperature $T$ , [41]. Therefore, the heat flow rate is $\Delta \dot Q_{th} = T \Delta \dot S$ . The dimension of $\Delta \dot S$ , the entropy rate, is then [ $W/K$ ]. Therefore the dimension of the power variables, or ( $T \Delta \dot S$ ) is [ $W$ ], as it should be; hence, we can accept the temperature as the effort and entropy rate as the flow (or entropy as displacement) in BG model for a thermal system. Table 11-5 shows the variables involved in Pseudo BG and BG models for thermal systems.

Table 11-5 Definitions for effort and flow for BG and pseudo BG models, thermal systems
Model	Effort	Flow	Power	System
Pseudo BG	temperature	heat flow rate	ill-defined	single thermal domain
BG	temperature	entropy rate	defined	multi-domain

$d \mathcal{E} (p,q,t) = \dfrac{\partial \mathcal{E}}{\partial p} dp + \dfrac{\partial \mathcal{E}}{\partial q} dq + \dfrac{\partial \mathcal{E}}{\partial t} dt$

	References
[1]	H. M. Paynter, Analysis and Design of Engineering Systems, Boston, MA: MIT Press, 1961.
[2]	Ronald Zamora-Musa and Jesus Coba, “Trends in modeling and simulation in the automotive industry concerning the bond graph framework,” International Journal of Engineering & Technology, vol. 7, pp. 160-163, 2018.
[3]	Controllab, Introduction to 20-sim, Accessed: 2020. [Online]. Available: https://www.20sim.com/overview/introduction/
[4]	M. Tiller, Introduction to Physical Modeling with Modelica, Kluwer Academic Publishers, 2001.
[5]	Charles M. Close and Dean K. Frederick, Modeling and Analysis of Dynamic Systems, 3rd ed., Wiley, 2015.
[6]	Modelica Association, Modelica Tools, Accessed: 2020. [Online]. Available: https://www.modelica.org/tools.
[7]	J. Lagrange, Méchanique Analytique; English, translated and edited by Auguste Boissonnade and Victor N. Vagliente., Boston: Kluwer Academic Publishers, 1997.
[8]	P. Hamill, A Student’s Guide to Lagrangians and Hamiltonians, Cambridge University Press, 2013.
[9]	E. Khutoryansky, Euler-Lagrange equation explained intuitively – Lagrangian Mechanics (Nov. 18 2018).[Online Video]. Available: https://www.youtube.com/watch?v=EceVJJGAFFI
[10]	J.E. Marsden and T.S. Ratiu, Introduction to Mechanics and Symmetry – A Basic Exposition of Classical Mechanical Systems, 2nd edition, Springer, 2002.
[11]	Goldstein, H., Poole Jr, C.P., and Charles, P, Classical Mechanics, 3rd ed., Pearson, 2013.
[12]	E. Kreyszig, Advanced Engineering Mathematics, 10th ed., Wiley, 2011.
[13]	Tom W.B Kibble and Frank H. Berkshire, Classical Mechanics, 5th ed., Imperial College Press, 2014.
[14]	W. Hamilton, “On a General Method in Dynamics,” Philosophical Transactions of the Royal Society, pp. Part I, pp. 95-144, Part II, pp. 247-308, 1834,1835.
[15]	K. Cassel, Variational Methods with Applications in Science and Engineering, Cambridge University Press, 2013.
[16]	W. Greiner, Classical Mechanics, Springer-Verlag, 2010.
[17]	R. d. Callafon, , MAE courses. [Online]. UCSD Available: http://maecourses.ucsd.edu/callafon/labcourse/lecturenotes/Lagrange_Handout.pdf. [Accessed Oct., 2019].
[18]	M. Tabatabaian, Tensor Analysis for Engineers, 2nd ed. Mercury and Learning Information, 2021.
[19]	N. A. Lemos, “Atwood’s machine with a massive string,” European Journal of Physics, vol. 38, no. 6, Sep., 2017. doi: arXiv:1710.01263V1.
[20]	Dean C. Karnopp, Donald L. Margolis, and Ronald C. Rosenberg, System Dynamics: Modeling and Simulation of Mechatronic Systems, 4th ed. Wiley Publishing, 2006.
[21]	S. Das, Mechatronic Modeling and Simulation Using Bond Graphs, 1st ed. Boca Raton: CRC Press, 2009.
[22]	B. T. Forbes, Engineering Systems Dynamics: A Unified Graph-Centered Approach, 2nd ed. CRC Press., 2007.
[23]	James L. Meriam, L. G. Kraige, J. N. Bolton, Engineering mechanics: Dynamics, 8th Edition, Wiley, 2015.
[24]	Charles Alexander and Matthew Sadiku, Fundamentals of Electric Circuits, 7th ed. McGraw-Hill Education, 2020.
ED.[25]	NASA, “Wheatstone Bridge Circuit,” N. Hall, Ed. 2015. [Online]. Available: https://www.grc.nasa.gov/WWW/K-12/airplane/tunwheat.html. [Accessed 2020].
[26]	I. Shames, Mechanics of Fluids, 4th ed. McGraw-Hill, 2002.
[27]	F. M. White, Fluid Mechanics, McGraw-Hill, 2015.
[28]	M. Tabatabaian, CFD Module: Turbulent Flow Modeling, Mercury Learning and Information, 2015.
[29]	Yilei Li, Zhencai Zhu, and Guoan Chen, “Bond Graph Modeling and Validation of an Energy Regenerative System for Emulsion Pump Tests,” The Scientific World Journal, vol. 2014, pp. 1-12, 2014.
[30]	J. A. Kypuros, System Dynamics and Control with Bond Graph Modeling, CRC Press, 2013.
[31]	N. S. Nise, Control Systems Engineering, 8th ed., John Wiley & Sons, 2019.
[32]	P. Gawthrop and L. Smith, Metamodelling: For Bond Graphs and Dynamic Systems, Prentice Hall, 1996.
[33]	W. Borutzky, “Bond Graph Modelling and Simulation of Mechatronic Systems- An Introduction into the Methology,” in 20th Europeann Conference on Modelling and Simulation, Bonn, 2006.
[34]	Fan, Zhun, et. al., “Exploring Multiple Design Topologies Using Genetic Programming And Bond Graphs,” in GECCO 2002: Proceedings of the Genetic and Evolutionary Computation, New York, 2002.
[35]	G. Strang, Introduction to Applied Mathematics, Wellesley-Cambridge Press, 1986.
[36]	K. Ogata, Modern Control Engineering, 5th ed. Pearson, 2009.
[37]	H. W. Bode, Network Analysis and Feedback Amplifier Design, D. Van Nostrand, 1945.
[38]	B. York, “Frequency Response and Bode Plots,” 2009. [Online]. Available: https://web.njit.edu/~levkov/classes_files/ECE232/Handouts/Frequency%20Response.pdf. [Accessed 2020].
[39]	John H. Lienhard V and John H. Lienhard IV, A Heat Transfer Textbook, 5th ed. Cambridge, Massachusetts: Phlogiston Press, 2019.
[40]	Y. A. Cengel, Heat Transfer: A Practical Approach, McGraw-Hill, 2002.
[41]	M. Tabatabaian and ER.R.K. Rajput, Advanced Thermodynamics: Fundamentals, Mathematics, Applications, Mercury Learning and Information, 2017.
[42]	G. Romero, et. al., “Kinematic Analysis of Mechanism By Using Bond-Graph Language,” in 20th European Conference on Modelling and Simulation, 2006.

Version	Date	Change
1.0 1.1	October 2021 May 2022	PressBook Created Typos, Figures, some Exercises

Engineering Systems Dynamics, Modelling, Simulation, and Design

Engineering Systems Dynamics, Modelling, Simulation, and Design

Lagrangian and Bond Graph Methods

Epigraph

Contents

Acknowledgements

About the Author

Preface

List of Figures

List of Tables

Introduction

Media Attributions

Lagrangian Mechanics

2.1 Overview

2.2 Example: A Mass-Spring System

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

2.4 Generalized Coordinates, Momenta, and Forces

2.5 Hamilton’s Principle and Lagrange’s Equations

2.6 Cyclic Coordinates

2.7 Conservative and Non-Conservative Forces

2.8 Alternative form of Lagrange’s Equation

2.9 Multi-Domain Systems

2.10 Systems with Higher Order Equations

2.11 Example: A Multi-Mass-Spring System

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

2.13 Example: A Two-Loop Electrical Circuit

2.14 Example: A Compound Atwood’s Machine

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

2.16 Example: A Complex Vibrating Mechanical System

2.17 Example: A Pendulum with Moving Pivot

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

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

2.20 Example: An Extensible Robotic Arm Rotating in a Plane

Exercise Problems for Chapter 2

Media Attributions

Bond Graph Modelling Method

3.1 Overview

3.2 Categorizing System Components—Generalized Effort and Flow

3.3 Causality Principle and Assignment

3.4 Nine Basic Elements of Bond Graph Method

3.4.1 Inertia Element I: Kinetic Energy Storage

3.4.2 Capacity Element C: Potential Energy Storage Element

3.4.3 Friction Element R: Energy Dissipation Element

3.4.4 Source Elements Se and Sf : System Boundary Input Elements

3.4.5 1- and 0-junctions: Distribution Constraint Elements

3.4.6 Transformer TF and Gyrator GY: Energy Conversion Elements

3.5 System State Variables

3.5.1 Integral Causality and State Variables: I– and C-elements

Exercise Problems for Chapter 3

Media Attributions

Building Bond Graph Models: General Procedure and Application

4.1 Overview

4.2 Steps for Building Bond Graph Models: General Guidelines

4.2.1 Guidelines for Power Bond Direction

4.2.2 Guidelines for Assigning Causality Strokes

4.3 Example: BG Model for a One-DOF Mass-Spring-Damper Mechanical System

4.4 Example: BG Model for a Two-DOF Mass-Spring-Damper Mechanical System

4.5 Example: BG Model for a Three-DOF Mass-Spring-Damper Mechanical System

4.6 Example: Kinetics and Kinematics of a Mechanical System Using BG Model

4.7 Modelling and Simulation Approaches in Engineering: Modern vs. Traditional

4.7.1 Traditional Approach

4.7.2 Modern Approach

Exercise Problems for Chapter 4

Media Attributions

Introduction to 20-sim Software Tool

5.1 Overview

5.2 Example: BG Model for a Car Seat Mechanical System

5.3 Example: BG Model for a Cart Carrying a Load

Exercise Problems for Chapter 5

Bond Graph Models for Complex Mechanical Systems

6.1 Overview

6.2 Mechanical Systems—Rotational

6.3 Mechanical Systems—Two-Dimensional Rigid Plane Motion

6.4 Example: Gear-Shaft Mechanical System—Rotational

6.5 Example: Double Rack-and-Pinion Mechanical System—Rotational

6.6 Example: Mass-Spring-Damper System on an Inclined Plane—Translational

6.7 Example: 2D Rigid-Body Motion—Half-Car Model

6.8 Example: Mass-Spring-Damper System Connected to a Massless Lever

6.9 Example: Mass-Spring-Damper System Connected to a Lever

6.10 Example: Inclined Lever and Mass-Spring-Damper System

3.4.4 Source Elements S_e and S_f : System Boundary Input Elements