Bond graphs provide a unified modeling language for mechanical, electrical, fluid, and thermal systems. Instead of starting from domain-specific laws only, we begin with the universal concept of power transfer. A bond carries two variables: \( e \) (effort) and \( f \) (flow).

\[ P = e f \]

The sign of \( P \) depends on the bond orientation (reference power direction). Positive power means energy flows in the direction of the bond arrow.

flowchart TD
  A["Physical subsystem"] --> B["Choose effort/flow variables"]
  B --> C["Assign bonds and directions"]
  C --> D["Pick elements: Se, Sf, R, C, I, TF, GY"]
  D --> E["Connect with 0/1 junctions"]
  E --> F["Write junction constraints"]
  F --> G["Write constitutive laws"]
  G --> H["Derive state equations"]
  H --> I["Check power and energy balance"]
        

The power-product structure \( P=e f \) makes bond graphs domain independent. Common pairings include:

• Translational mechanics: \( e=F \), \( f=v \)
• Rotational mechanics: \( e=T \), \( f=\omega \)
• Electrical systems: \( e=V \), \( f=i \)
• Hydraulic systems: \( e=p \), \( f=Q \)
• Thermal systems (strict power form): \( e=T \), \( f=\dot{S} \)

Storage elements introduce integrated variables. For a \( C \)-element the state is \( q \), and for an \( I \)-element the state is \( p \):

\[ q(t)=q(t_0)+\int_{t_0}^{t} f(\xi)\,d\xi, \qquad p(t)=p(t_0)+\int_{t_0}^{t} e(\xi)\,d\xi \]

These variables will become the natural state variables in bond-graph-based state-space models.

3.1 Sources \(Se\) and \(Sf\)

An effort source imposes effort and a flow source imposes flow:

\[ Se:\; e(t)=e_s(t), \qquad Sf:\; f(t)=f_s(t) \]

Their instantaneous powers are \( P_{Se}=e_s f \) and \( P_{Sf}=e f_s \). Depending on sign, sources can inject or absorb energy.

3.2 Resistive Element \(R\)

A resistive element models dissipation (damper, resistor, hydraulic resistance, etc.). In the linear case:

\[ e = R f \qquad \text{or} \qquad f = \frac{1}{R}e \]

Passivity proof: for \( R > 0 \),

\[ P_R = e f = R f^2 = \frac{e^2}{R} \ge 0 \]

Therefore a positive linear \( R \) element cannot generate net power.

3.3 Capacitive Element \(C\)

A \( C \)-element stores potential-type energy. Its constitutive law is naturally expressed through an energy function \( H_C(q) \):

\[ \dot{q}=f, \qquad e=\frac{\partial H_C(q)}{\partial q} \]

For a linear \( C \)-element with compliance parameter \( C \):

\[ H_C(q)=\frac{q^2}{2C}, \qquad e=\frac{q}{C} \]

Energy-rate proof:

\[ \frac{dH_C}{dt} = \frac{\partial H_C}{\partial q}\dot{q} = e f \]

3.4 Inertial Element \(I\)

An \( I \)-element stores kinetic-type energy:

\[ \dot{p}=e, \qquad f=\frac{\partial H_I(p)}{\partial p} \]

For a linear \( I \)-element with inertance parameter \( I \):

\[ H_I(p)=\frac{p^2}{2I}, \qquad f=\frac{p}{I} \]

Energy-rate proof:

\[ \frac{dH_I}{dt} = \frac{\partial H_I}{\partial p}\dot{p} = f e \]

Junctions encode ideal interconnection constraints. Let \( \sigma_i \in \{-1,+1\} \) denote the sign of bond \( i \) with respect to a junction orientation convention.

4.1 0-Junction (Common Effort)

\[ e_1=e_2=\cdots=e_m=e_J, \qquad \sum_{i=1}^{m}\sigma_i f_i = 0 \]

Power-conservation proof:

\[ \sum_{i=1}^{m}\sigma_i P_i = \sum_{i=1}^{m}\sigma_i e_i f_i = e_J \sum_{i=1}^{m}\sigma_i f_i = 0 \]

4.2 1-Junction (Common Flow)

\[ f_1=f_2=\cdots=f_m=f_J, \qquad \sum_{i=1}^{m}\sigma_i e_i = 0 \]

Power-conservation proof:

\[ \sum_{i=1}^{m}\sigma_i P_i = \sum_{i=1}^{m}\sigma_i e_i f_i = f_J \sum_{i=1}^{m}\sigma_i e_i = 0 \]

This is the bond-graph counterpart of ideal Kirchhoff-type interconnection constraints.

Bond graphs use two ideal power-conserving two-port elements for domain conversion and power scaling.

5.1 Transformer \(TF\)

For modulus \( n \), a standard constitutive form is

\[ e_2 = n e_1, \qquad f_1 = n f_2 \]

Power invariance:

\[ e_1 f_1 = e_1 (n f_2) = (n e_1)f_2 = e_2 f_2 \]

This represents ideal gears, levers, and other kinematic transformations.

5.2 Gyrator \(GY\)

For modulus \( r \), one common form is

\[ e_2 = r f_1, \qquad e_1 = r f_2 \]

Power invariance:

\[ e_1 f_1 = (r f_2)f_1 = (r f_1)f_2 = e_2 f_2 \]

Gyrators are important in electromechanical energy conversion models.

Consider a force-driven mass–spring–damper system with source force \( F_s(t) \), mass \( m \), damping \( b \), and stiffness \( k \). In bond-graph notation, the elements are \( Se \), \( I=m \), \( R=b \), and \( C=1/k \).

flowchart LR
  Se["Se: F_s(t)"] --> J1["1-junction (common v)"]
  J1 --> I["I: m"]
  J1 --> R["R: b"]
  J1 --> C["C: 1/k"]
        

At the 1-junction, all flows are equal: \( f=v \). Choose states \( q \) (for \( C \)) and \( p \) (for \( I \)).

\[ v=\frac{p}{m}, \qquad e_R=bv=b\frac{p}{m}, \qquad e_C=\frac{q}{C}=kq \]

Effort balance at the 1-junction gives

\[ F_s(t) - e_I - e_R - e_C = 0 \]

Since \( e_I=\dot{p} \), we obtain

\[ \dot{p} = F_s(t) - b\frac{p}{m} - \frac{q}{C} \]

The \( C \)-element equation gives

\[ \dot{q}=v=\frac{p}{m} \]

Therefore the bond-graph state model is

\[ \begin{bmatrix}\dot{q}\\ \dot{p}\end{bmatrix} = \begin{bmatrix} 0 & \frac{1}{m}\\ -\frac{1}{C} & -\frac{b}{m} \end{bmatrix} \begin{bmatrix}q\\p\end{bmatrix} + \begin{bmatrix}0\\1\end{bmatrix}F_s(t) \]

Using \( C=1/k \) and \( p=m\dot{q} \), this reduces to the familiar ODE

\[ m\ddot{q} + b\dot{q} + kq = F_s(t) \]

Stored energy is

\[ H(q,p)=\frac{p^2}{2m}+\frac{q^2}{2C} \]

and its rate satisfies the balance law

\[ \dot{H}=F_s(t)\,v - b v^2 \]

which explicitly separates source power and dissipated power.

The following Python program implements the bond-graph state equations in variables \( (q,p) \), performs RK4 integration, and checks the numerical power balance.

Code: Chapter18_Lesson4.py

# Chapter18_Lesson4.py
# Bond graph simulation: Se - 1 - (I,R,C) for a mass-spring-damper

import numpy as np
import matplotlib.pyplot as plt

def Se(t):
    return 2.0*np.sin(2*np.pi*0.8*t) + (1.0 if t >= 1.0 else 0.0)

def f_bg(t, x, m, b, C):
    q, p = x                     # states: q (C-state), p (I-state)
    v = p / m                    # common flow at 1-junction
    eR = b * v
    eC = q / C                   # = k q because C = 1/k
    qdot = v
    pdot = Se(t) - eR - eC       # 1-junction effort balance
    return np.array([qdot, pdot], dtype=float)

def rk4(fun, t, x, h, *params):
    k1 = fun(t, x, *params)
    k2 = fun(t + h/2, x + h*k1/2, *params)
    k3 = fun(t + h/2, x + h*k2/2, *params)
    k4 = fun(t + h, x + h*k3, *params)
    return x + (h/6)*(k1 + 2*k2 + 2*k3 + k4)

def main():
    m = 1.5
    k = 12.0
    b = 1.2
    C = 1.0 / k

    h = 1e-3
    t = np.arange(0.0, 10.0 + h, h)
    x = np.zeros((len(t), 2))
    x[0] = [0.10, 0.0]  # q(0), p(0)

    for i in range(len(t)-1):
        x[i+1] = rk4(f_bg, t[i], x[i], h, m, b, C)

    q = x[:, 0]
    p = x[:, 1]
    v = p / m
    eS = np.array([Se(ti) for ti in t])

    H = 0.5*p**2/m + 0.5*q**2/C
    Pdiss = b*v**2
    Ediss = np.cumsum(Pdiss) * h
    dH = np.gradient(H, h)
    residual = eS*v - Pdiss - dH

    print('Max |power residual| =', float(np.max(np.abs(residual))))

    np.savetxt(
        'Chapter18_Lesson4_python_output.csv',
        np.c_[t, q, v, H, Ediss, residual],
        delimiter=',',
        header='t,q,v,H,Ediss,residual',
        comments=''
    )

    plt.figure()
    plt.plot(t, q, label='q(t)')
    plt.plot(t, v, label='v(t)')
    plt.grid(True); plt.legend(); plt.title('States (bond-graph form)')

    plt.figure()
    plt.plot(t, H, label='Stored energy H')
    plt.plot(t, Ediss, label='Dissipated energy')
    plt.grid(True); plt.legend(); plt.title('Energy accounting')

    plt.figure()
    plt.plot(t, residual)
    plt.grid(True); plt.title('Power-balance residual')
    plt.show()

if __name__ == '__main__':
    main()

These versions implement exactly the same bond-graph dynamics and export CSV outputs for comparison.

Code: Chapter18_Lesson4.cpp

// Chapter18_Lesson4.cpp
// Bond graph simulation: Se - 1 - (I,R,C) for a mass-spring-damper
// Compile: g++ -std=c++17 Chapter18_Lesson4.cpp -o Chapter18_Lesson4

#include <array>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <vector>

struct Params { double m, b, C; };

double Se(double t) {
    return 2.0 * std::sin(2.0 * M_PI * 0.8 * t) + (t >= 1.0 ? 1.0 : 0.0);
}

std::array<double,2> rhs(double t, const std::array<double,2>& x, const Params& p) {
    double q = x[0];
    double mom = x[1];
    double v = mom / p.m;
    double eR = p.b * v;
    double eC = q / p.C;
    return {v, Se(t) - eR - eC};
}

std::array<double,2> rk4(double t, const std::array<double,2>& x, double h, const Params& p) {
    auto k1 = rhs(t, x, p);
    auto x2 = std::array<double,2>{x[0] + 0.5*h*k1[0], x[1] + 0.5*h*k1[1]};
    auto k2 = rhs(t + 0.5*h, x2, p);
    auto x3 = std::array<double,2>{x[0] + 0.5*h*k2[0], x[1] + 0.5*h*k2[1]};
    auto k3 = rhs(t + 0.5*h, x3, p);
    auto x4 = std::array<double,2>{x[0] + h*k3[0], x[1] + h*k3[1]};
    auto k4 = rhs(t + h, x4, p);

    return {
        x[0] + (h/6.0)*(k1[0] + 2*k2[0] + 2*k3[0] + k4[0]),
        x[1] + (h/6.0)*(k1[1] + 2*k2[1] + 2*k3[1] + k4[1])
    };
}

int main() {
    Params p{1.5, 1.2, 1.0/12.0};
    double h = 1e-3, t0 = 0.0, tf = 10.0;
    int N = static_cast<int>((tf - t0)/h) + 1;

    std::vector<double> t(N), q(N), mom(N), v(N), H(N), Ediss(N), res(N);
    q[0] = 0.10; mom[0] = 0.0; Ediss[0] = 0.0;
    for (int i = 0; i < N; ++i) t[i] = t0 + i*h;

    for (int i = 0; i < N - 1; ++i) {
        auto xn = rk4(t[i], {q[i], mom[i]}, h, p);
        q[i+1] = xn[0];
        mom[i+1] = xn[1];
    }

    for (int i = 0; i < N; ++i) {
        v[i] = mom[i]/p.m;
        H[i] = 0.5*mom[i]*mom[i]/p.m + 0.5*q[i]*q[i]/p.C;
    }

    for (int i = 1; i < N; ++i) {
        Ediss[i] = Ediss[i-1] + h * (p.b * v[i-1] * v[i-1]);
    }

    double maxRes = 0.0;
    for (int i = 1; i < N - 1; ++i) {
        double dH = (H[i+1] - H[i-1])/(2*h);
        res[i] = Se(t[i])*v[i] - p.b*v[i]*v[i] - dH;
        if (std::abs(res[i]) > maxRes) maxRes = std::abs(res[i]);
    }
    std::cout << "Max |power residual| = " << maxRes << "\n";

    std::ofstream out("Chapter18_Lesson4_cpp_output.csv");
    out << "t,q,v,H,Ediss,residual\n";
    out << std::setprecision(10);
    for (int i = 0; i < N; ++i) {
        out << t[i] << "," << q[i] << "," << v[i] << "," << H[i] << "," << Ediss[i] << "," << res[i] << "\n";
    }
    out.close();
    return 0;
}

Code: Chapter18_Lesson4.java

// Chapter18_Lesson4.java
// Bond graph simulation: Se - 1 - (I,R,C) for a mass-spring-damper
// Compile: javac Chapter18_Lesson4.java
// Run:     java Chapter18_Lesson4

import java.io.FileWriter;
import java.io.PrintWriter;

public class Chapter18_Lesson4 {
    static class Params {
        double m, b, C;
        Params(double m, double b, double C) { this.m = m; this.b = b; this.C = C; }
    }

    static double Se(double t) {
        return 2.0 * Math.sin(2.0 * Math.PI * 0.8 * t) + (t >= 1.0 ? 1.0 : 0.0);
    }

    static double[] rhs(double t, double[] x, Params p) {
        double q = x[0], mom = x[1];
        double v = mom / p.m;
        double eR = p.b * v;
        double eC = q / p.C;
        return new double[]{v, Se(t) - eR - eC};
    }

    static double[] rk4(double t, double[] x, double h, Params p) {
        double[] k1 = rhs(t, x, p);
        double[] x2 = new double[]{x[0] + 0.5*h*k1[0], x[1] + 0.5*h*k1[1]};
        double[] k2 = rhs(t + 0.5*h, x2, p);
        double[] x3 = new double[]{x[0] + 0.5*h*k2[0], x[1] + 0.5*h*k2[1]};
        double[] k3 = rhs(t + 0.5*h, x3, p);
        double[] x4 = new double[]{x[0] + h*k3[0], x[1] + h*k3[1]};
        double[] k4 = rhs(t + h, x4, p);

        return new double[]{
            x[0] + (h/6.0)*(k1[0] + 2*k2[0] + 2*k3[0] + k4[0]),
            x[1] + (h/6.0)*(k1[1] + 2*k2[1] + 2*k3[1] + k4[1])
        };
    }

    public static void main(String[] args) throws Exception {
        Params p = new Params(1.5, 1.2, 1.0/12.0);
        double h = 1e-3, t0 = 0.0, tf = 10.0;
        int N = (int)((tf - t0)/h) + 1;

        double[] t = new double[N], q = new double[N], mom = new double[N], v = new double[N];
        double[] H = new double[N], Ediss = new double[N], res = new double[N];

        q[0] = 0.10; mom[0] = 0.0; Ediss[0] = 0.0;
        for (int i = 0; i < N; i++) t[i] = t0 + i*h;

        for (int i = 0; i < N - 1; i++) {
            double[] xn = rk4(t[i], new double[]{q[i], mom[i]}, h, p);
            q[i+1] = xn[0];
            mom[i+1] = xn[1];
        }

        for (int i = 0; i < N; i++) {
            v[i] = mom[i] / p.m;
            H[i] = 0.5*mom[i]*mom[i]/p.m + 0.5*q[i]*q[i]/p.C;
        }

        for (int i = 1; i < N; i++) {
            Ediss[i] = Ediss[i-1] + h * (p.b * v[i-1] * v[i-1]);
        }

        double maxRes = 0.0;
        for (int i = 1; i < N - 1; i++) {
            double dH = (H[i+1] - H[i-1]) / (2*h);
            res[i] = Se(t[i]) * v[i] - p.b*v[i]*v[i] - dH;
            maxRes = Math.max(maxRes, Math.abs(res[i]));
        }
        System.out.println("Max |power residual| = " + maxRes);

        PrintWriter out = new PrintWriter(new FileWriter("Chapter18_Lesson4_java_output.csv"));
        out.println("t,q,v,H,Ediss,residual");
        for (int i = 0; i < N; i++) {
            out.printf(java.util.Locale.US, "%.10f,%.10f,%.10f,%.10f,%.10f,%.10f%n",
                    t[i], q[i], v[i], H[i], Ediss[i], res[i]);
        }
        out.close();
    }
}

The MATLAB version uses ode45 and the same bond-graph state variables \( q \) and \( p \). In Simulink, the equivalent model is obtained with two integrators and algebraic gain/sum blocks. The Mathematica version uses NDSolveValue.

Code: Chapter18_Lesson4.m

% Chapter18_Lesson4.m
% Bond graph simulation: Se - 1 - (I,R,C) for a mass-spring-damper

clear; clc; close all;

m = 1.5; b = 1.2; k = 12.0; C = 1/k;
x0 = [0.10; 0.0]; % [q0; p0]

[t, X] = ode45(@(t,x) rhs_bg(t, x, m, b, C), [0 10], x0);

q = X(:,1);
p = X(:,2);
v = p/m;
eS = arrayfun(@Se, t);

H = 0.5*(p.^2)/m + 0.5*(q.^2)/C;
Pdiss = b*(v.^2);
Ediss = cumtrapz(t, Pdiss);
dH = gradient(H, t);
residual = eS.*v - Pdiss - dH;

disp(['Max |power residual| = ', num2str(max(abs(residual)))]);

T = table(t, q, v, H, Ediss, residual);
writetable(T, 'Chapter18_Lesson4_matlab_output.csv');

figure; plot(t,q,t,v); grid on; legend('q','v'); title('Bond-graph states');
figure; plot(t,H,t,Ediss); grid on; legend('H','E_d_i_s_s'); title('Energy accounting');
figure; plot(t,residual); grid on; title('Power-balance residual');

function dx = rhs_bg(t, x, m, b, C)
    q = x(1); p = x(2);
    v = p/m;
    eR = b*v;
    eC = q/C;
    dx = [v; Se(t) - eR - eC];
end

function y = Se(t)
    y = 2*sin(2*pi*0.8*t) + double(t >= 1.0);
end

Code: Chapter18_Lesson4.nb

(* Chapter18_Lesson4.nb *)
(* Wolfram Language code saved with .nb extension *)

ClearAll["Global`*"];
m = 1.5; b = 1.2; k = 12.0; C = 1/k;
Se[t_] := 2 Sin[2 Pi 0.8 t] + If[t >= 1.0, 1.0, 0.0];

eqns = {
  q'[t] == p[t]/m,
  p'[t] == Se[t] - b (p[t]/m) - q[t]/C,
  q[0] == 0.10,
  p[0] == 0.0
};

sol = NDSolveValue[eqns, {q, p}, {t, 0, 10}];
qf[t_] := sol[[1]][t];
pf[t_] := sol[[2]][t];
vf[t_] := pf[t]/m;

H[t_] := 0.5 pf[t]^2/m + 0.5 qf[t]^2/C;
Pdiss[t_] := b vf[t]^2;
Ediss[t_] := NIntegrate[Pdiss[s], {s, 0, t}];

times = Range[0, 10, 0.01];
data = Table[
  Module[{dH = (H[tt + 10^-4] - H[tt - 10^-4])/(2*10^-4)},
    {tt, qf[tt], vf[tt], H[tt], Ediss[tt], Se[tt] vf[tt] - Pdiss[tt] - dH}
  ],
  {tt, times}
];

Export["Chapter18_Lesson4_mathematica_output.csv",
  Prepend[data, {"t","q","v","H","Ediss","residual"}]
];

Plot[{qf[t], vf[t]}, {t, 0, 10}, GridLines -> Automatic, PlotLegends -> {"q","v"}]
Plot[{H[t], Ediss[t]}, {t, 0, 10}, GridLines -> Automatic, PlotLegends -> {"H","Ediss"}]
Plot[Se[t] vf[t] - Pdiss[t] - D[H[s], s] /. s -> t, {t, 0.01, 9.99}, GridLines -> Automatic]

Problem 1 (0-Junction Conservation): Show that an ideal 0-junction is power-conserving.

Solution: At a 0-junction all efforts are equal: \( e_i=e_J \), and flows satisfy \( \sum_i \sigma_i f_i = 0 \). Therefore

\[ \sum_i \sigma_i P_i = \sum_i \sigma_i e_i f_i = e_J \sum_i \sigma_i f_i = 0 \]

Hence no net power is created or destroyed at an ideal 0-junction.

Problem 2 (1-Junction Conservation): Show that an ideal 1-junction is power-conserving.

Solution: At a 1-junction all flows are equal: \( f_i=f_J \), and efforts satisfy \( \sum_i \sigma_i e_i = 0 \). Then

\[ \sum_i \sigma_i P_i = \sum_i \sigma_i e_i f_i = f_J \sum_i \sigma_i e_i = 0 \]

Therefore the ideal 1-junction is also lossless.

Problem 3 (Linear \(C\)-Element Energy Function): Given \( e=q/C \), derive \( H_C(q) \).

Solution: Using \( e=\partial H_C/\partial q \),

\[ \frac{\partial H_C}{\partial q} = \frac{q}{C} \]

Integrating with respect to \( q \):

\[ H_C(q)=\int \frac{q}{C}\,dq = \frac{q^2}{2C} + \text{constant} \]

With zero reference energy at \( q=0 \), we obtain \( H_C(q)=q^2/(2C) \).

Problem 4 (Mass–Spring–Damper from Bond Graph): Starting from the 1-junction bond graph in Section 6, derive the state equations and show equivalence to the classical ODE.

Solution: Common flow is \( v=p/m \), so \( \dot{q}=v=p/m \). Resistive and capacitive efforts are \( e_R=b p/m \) and \( e_C=q/C \). Junction effort balance gives

\[ \dot{p} = F_s(t) - b\frac{p}{m} - \frac{q}{C} \]

Since \( p=m\dot{q} \) and \( C=1/k \),

\[ m\ddot{q} + b\dot{q} + kq = F_s(t) \]

which matches the standard second-order mechanical model.

Problem 5 (Transformer Losslessness): For \( e_2=n e_1 \) and \( f_1=n f_2 \), prove that an ideal transformer is lossless.

Solution:

\[ P_1 = e_1 f_1 = e_1(n f_2)=n e_1 f_2, \qquad P_2 = e_2 f_2 = (n e_1)f_2 = n e_1 f_2 \]

Hence \( P_1=P_2 \), so the transformer preserves power.

Bond graphs unify system dynamics modeling through effort/flow variables and power bonds. In this lesson, we defined the core one-port and two-port elements, proved junction conservation laws, and derived a complete state-space model of a mass–spring–damper system in bond-graph coordinates \( (q,p) \). This provides a strong foundation for later multi-domain mechatronic and robotics models, where energy consistency is crucial.