In earlier chapters, we modeled systems primarily through force balances, KCL/KVL, and state equations. In this chapter, we add a complementary perspective: energy bookkeeping. The central idea is: every physically valid lumped model must satisfy a power balance of the form \( \text{input power} = \text{rate of stored energy} + \text{dissipated/transferred power} \).

If a model violates this balance, it usually indicates a sign error, inconsistent constitutive law, or incorrect interconnection. This makes energy methods not only theoretical but also extremely practical for simulation verification.

\[ P(t) = \frac{dW}{dt}, \qquad W(t_1,t_2)=\int_{t_1}^{t_2} P(t)\,dt \]

We will repeatedly use power-conjugate variables: one variable is an effort-like quantity (force, voltage, pressure, temperature), and the other is a flow-like quantity (velocity, current, volumetric flow, entropy flow). Their product gives power.

flowchart TD
  A["Boundary source"] --> B["Interconnection constraints"]
  B --> C["Storage elements"]
  B --> D["Dissipative elements"]
  C --> E["dE/dt"]
  D --> F["Loss / transfer power"]
  E --> G["Power balance check"]
  F --> G
        

Let \( e(t) \) denote an effort variable and \( f(t) \) denote a flow variable. The instantaneous power is

\[ p(t)=e(t)f(t) \]

and the work supplied over \( [t_1,t_2] \) is

\[ W = \int_{t_1}^{t_2} e(t)f(t)\,dt. \]

In a lumped model, we partition power into storage and dissipation:

\[ p_{in} = \frac{d}{dt}E_{stored} + p_{diss}. \]

This relation is exact for mechanical and electrical lumped models with standard constitutive laws. For thermal systems, care is needed because the exact conjugate pair is \( (T,\dot{S}) \) rather than \( (T,\dot{Q}) \); we will treat both the rigorous and the engineering-lumped viewpoints in Section 6.

For translational motion, the power variables are \( F \) (force) and \( v \) (velocity):

\[ p = Fv, \qquad W = \int F\,dx = \int Fv\,dt. \]

For rotational motion, the analogous variables are \( T \) (torque) and \( \omega \) (angular velocity):

\[ p = T\omega, \qquad W = \int T\,d\theta. \]

Kinetic energy theorem (translational). Starting from Newton's law \( m\dot{v}=\sum F \), multiply both sides by \( v \):

\[ m\dot{v}v = \left(\sum F\right)v \quad \Rightarrow \quad \frac{d}{dt}\left(\tfrac{1}{2}mv^2\right)=\sum p_i. \]

Thus the time derivative of kinetic energy equals the net power into the mass.

Spring and damper energies. For a linear spring \( F_s = kx \), the stored energy is

\[ E_s = \int_0^x k\xi\,d\xi = \tfrac{1}{2}kx^2. \]

For a viscous damper \( F_d = cv \), the dissipated power is

\[ p_d = F_d v = cv^2 \ge 0. \]

Mass-spring-damper power balance. Consider \( m\ddot{x} + c\dot{x} + kx = F_{in}(t) \). Multiply by \( \dot{x} \):

\[ m\ddot{x}\dot{x} + c\dot{x}^2 + kx\dot{x} = F_{in}\dot{x} \]

\[ \frac{d}{dt}\left(\tfrac{1}{2}m\dot{x}^2 + \tfrac{1}{2}kx^2\right) = F_{in}\dot{x} - c\dot{x}^2. \]

Hence the mechanical stored energy \( E_m = \tfrac{1}{2}m\dot{x}^2 + \tfrac{1}{2}kx^2 \) satisfies the standard balance law.

For an \( n \)-DOF mechanical system with generalized coordinate vector \( q \) and positive-definite mass matrix \( M(q) \), the total mechanical energy is

\[ E(q,\dot{q})=\tfrac{1}{2}\dot{q}^T M(q)\dot{q} + V(q), \]

and in many models with viscous dissipation \( D\dot{q} \) and generalized input \( u \), one obtains

\[ \dot{E} = \dot{q}^T u - \dot{q}^T D\dot{q}. \]

In electrical circuits, the power variables are \( v \) (voltage) and \( i \) (current):

\[ p = vi, \qquad w = \int vi\,dt. \]

Resistor (dissipation). For \( v=Ri \):

\[ p_R = vi = Ri^2 = \frac{v^2}{R} \ge 0. \]

Capacitor (storage). Using \( q=Cv \) and \( i=dq/dt \),

\[ p_C = vi = v\frac{dq}{dt}. \]

Since \( q=Cv \Rightarrow dq = C\,dv \), the stored energy is

\[ E_C = \int_0^q v(q')\,dq' = \int_0^v Cv'\,dv' = \tfrac{1}{2}Cv^2. \]

Inductor (storage). Using \( v=L\,di/dt \),

\[ p_L = vi = Li\frac{di}{dt} = \frac{d}{dt}\left(\tfrac{1}{2}Li^2\right), \]

so \( E_L = \tfrac{1}{2}Li^2 \).

RLC energy balance. For a series RLC driven by source \( v_s(t) \),

\[ L\dot{i} + Ri + \frac{q}{C} = v_s, \qquad \dot{q}=i. \]

Multiplying the first equation by \( i \) gives

\[ v_s i = \frac{d}{dt}\left(\tfrac{1}{2}Li^2 + \tfrac{1}{2}\frac{q^2}{C}\right) + Ri^2. \]

This is exactly the electrical counterpart of the mass-spring-damper balance.

For lumped hydraulic/fluid models, a standard power pair is \( p \) (pressure) and \( q \) (volumetric flow rate). The instantaneous hydraulic power is

\[ P_h = p q. \]

(Units check: \( \mathrm{Pa}\cdot\mathrm{m}^3/\mathrm{s} = \mathrm{N}/\mathrm{m}^2 \cdot \mathrm{m}^3/\mathrm{s} = \mathrm{N m}/\mathrm{s} = \mathrm{W} \).)

Hydraulic resistance. For a linearized resistive element \( \Delta p = R_h q \), the dissipated power is

\[ P_R = \Delta p\,q = R_h q^2 \ge 0. \]

Hydraulic inertance. For a lumped inertive element (fluid column),

\[ \Delta p = L_h \dot{q}, \]

so

\[ P_L = \Delta p\,q = L_h q\dot{q} = \frac{d}{dt}\left(\tfrac{1}{2}L_h q^2\right). \]

The term \( \tfrac{1}{2}L_h q^2 \) is the fluid kinetic-energy storage (in lumped form).

Hydraulic capacitance (compressibility / compliance). In small-signal form around an operating point, let \( \Delta V = C_h \Delta p \). Then the incremental stored energy is

\[ E_C = \int_0^{\Delta p} \Delta p'\,d(\Delta V') = \int_0^{\Delta p} \Delta p' C_h\,d\Delta p' = \tfrac{1}{2}C_h(\Delta p)^2. \]

Chamber-line example. Consider the linearized two-state model

\[ C_h\dot{p} = q_{in} - q, \qquad L_h\dot{q} = p - R_h q. \]

Multiply the first equation by \( p \) and the second by \( q \), then add:

\[ p q_{in} = \frac{d}{dt}\left(\tfrac{1}{2}C_h p^2 + \tfrac{1}{2}L_h q^2\right) + R_h q^2. \]

Again we recover the same structure: input power equals change of stored energy plus dissipation.

In thermal modeling, the quantity commonly called “heat flow” \( \dot{Q} \) already has units of power (W). For a lumped thermal capacitance with constant specific heat and no mechanical work term,

\[ C_{th}\dot{T} = \dot{Q}_{in} - \dot{Q}_{out}, \qquad C_{th} = mc_p. \]

The stored internal energy relative to a reference temperature \( T_{ref} \) is

\[ E_{th} = C_{th}(T - T_{ref}), \]

so the energy balance is exactly

\[ \frac{dE_{th}}{dt} = \dot{Q}_{in} - \dot{Q}_{out}. \]

Thermal resistance model. For two nodes at temperatures \( T_1 \) and \( T_2 \), a lumped thermal resistor uses

\[ \dot{Q} = \frac{T_1 - T_2}{R_{th}}. \]

This is the thermal analogue of Ohm’s law. However, thermodynamics is subtler than electrical/mechanical systems: the exact power-conjugate pair is temperature and entropy flow:

\[ \dot{Q}_{rev} = T\dot{S}. \]

For irreversible conduction, entropy production is nonnegative. In this lesson we use the standard lumped-energy engineering model; a deeper entropy-based treatment appears later when we discuss Hamiltonian and port-based formulations.

A useful mnemonic is that many domains share a common “storage + dissipation + source” structure. The element names differ, but the energy algebra is the same:

• Mechanical: mass/spring store energy, damper dissipates.
• Electrical: inductor/capacitor store energy, resistor dissipates.
• Fluid: inertance/compliance store energy, hydraulic resistance dissipates.
• Thermal: thermal capacitance stores energy, thermal resistance transfers heat to another node.

This motivates a generic state-space “energy audit” for simulations. If \( x \) is the state vector and \( E(x) \) is a chosen storage function, then a physically consistent model often satisfies

\[ \dot{E}(x) = u^T y - \Phi(x), \qquad \Phi(x) \ge 0 \]

where \( u^T y \) is supplied power and \( \Phi \) collects dissipation (or, in thermal networks, transfer to another lumped node/environment). Numerically, you can integrate both sides and compare them as a verification test.

flowchart TD
  A["Write state equations"] --> B["Identify source power P_in"]
  B --> C["Identify stored energy E(x)"]
  C --> D["Identify losses / transfer terms"]
  D --> E["Form dE/dt = P_in - P_loss"]
  E --> F["Simulate states"]
  F --> G["Integrate powers to get work"]
  G --> H["Compare energy balance residual"]
        

The following script simulates one representative lumped model in each domain and checks the integrated energy balance. It is intentionally explicit and easy to modify for homework experiments (change parameters, inputs, or time step).

Code: Chapter18_Lesson1.py

# Chapter18_Lesson1.py
# Work, Energy, and Power in Mechanical, Electrical, Fluid, and Thermal Systems
# The script simulates four lumped dynamic models and checks energy balances.

import numpy as np

def trapz(y, t):
    return np.trapz(y, t)

t0, tf, dt = 0.0, 12.0, 1e-4
t = np.arange(t0, tf + dt, dt)
n = len(t)

# ==========================================================
# 1) Mechanical: mass-spring-damper
#    m xdd + c xd + k x = F(t)
# ==========================================================
m, c, k = 1.5, 0.8, 12.0

def F_in(tt):
    return 2.5 * np.sin(1.2 * tt) + 1.2 * np.cos(0.4 * tt)

x = np.zeros(n)
v = np.zeros(n)
Pm_in = np.zeros(n)
Pm_diss = np.zeros(n)
Em = np.zeros(n)

for kidx in range(n - 1):
    f = F_in(t[kidx])
    a = (f - c * v[kidx] - k * x[kidx]) / m
    v[kidx + 1] = v[kidx] + dt * a
    x[kidx + 1] = x[kidx] + dt * v[kidx]

    Pm_in[kidx] = f * v[kidx]
    Pm_diss[kidx] = c * v[kidx] ** 2
    Em[kidx] = 0.5 * m * v[kidx] ** 2 + 0.5 * k * x[kidx] ** 2

# last sample bookkeeping
Pm_in[-1] = F_in(t[-1]) * v[-1]
Pm_diss[-1] = c * v[-1] ** 2
Em[-1] = 0.5 * m * v[-1] ** 2 + 0.5 * k * x[-1] ** 2

Wm_in = trapz(Pm_in, t)
Wm_diss = trapz(Pm_diss, t)
mech_balance_error = Em[-1] - Em[0] - (Wm_in - Wm_diss)

# ==========================================================
# 2) Electrical: series RLC
#    L didt + R i + q/C = v_in(t),   dq/dt = i
# ==========================================================
R, L, C = 2.0, 0.6, 0.25

def v_in(tt):
    return 5.0 * np.sin(2.0 * tt) + 2.0 * np.cos(0.5 * tt)

q = np.zeros(n)   # capacitor charge
i = np.zeros(n)
Pe_in = np.zeros(n)
Pe_diss = np.zeros(n)
Ee = np.zeros(n)

for kidx in range(n - 1):
    vin = v_in(t[kidx])
    di = (vin - R * i[kidx] - q[kidx] / C) / L
    dq = i[kidx]

    i[kidx + 1] = i[kidx] + dt * di
    q[kidx + 1] = q[kidx] + dt * dq

    Pe_in[kidx] = vin * i[kidx]
    Pe_diss[kidx] = R * i[kidx] ** 2
    Ee[kidx] = 0.5 * L * i[kidx] ** 2 + 0.5 * (q[kidx] ** 2) / C

Pe_in[-1] = v_in(t[-1]) * i[-1]
Pe_diss[-1] = R * i[-1] ** 2
Ee[-1] = 0.5 * L * i[-1] ** 2 + 0.5 * (q[-1] ** 2) / C

We_in = trapz(Pe_in, t)
We_diss = trapz(Pe_diss, t)
elec_balance_error = Ee[-1] - Ee[0] - (We_in - We_diss)

# ==========================================================
# 3) Fluid (small-signal hydraulic chamber + line)
#    C_h dp/dt = q_in(t) - q
#    L_h dq/dt = p - R_h q
#    Power pair: p * q   (pressure * volumetric flow)
# ==========================================================
C_h, L_h, R_h = 0.08, 0.15, 1.7

def q_in(tt):
    return 0.8 * np.sin(0.9 * tt) + 0.2 * np.cos(0.3 * tt)

p = np.zeros(n)   # pressure deviation
qf = np.zeros(n)  # line flow
Pf_in = np.zeros(n)
Pf_diss = np.zeros(n)
Ef = np.zeros(n)

for kidx in range(n - 1):
    qsrc = q_in(t[kidx])
    dp = (qsrc - qf[kidx]) / C_h
    dqf = (p[kidx] - R_h * qf[kidx]) / L_h

    p[kidx + 1] = p[kidx] + dt * dp
    qf[kidx + 1] = qf[kidx] + dt * dqf

    Pf_in[kidx] = p[kidx] * qsrc
    Pf_diss[kidx] = R_h * qf[kidx] ** 2
    Ef[kidx] = 0.5 * C_h * p[kidx] ** 2 + 0.5 * L_h * qf[kidx] ** 2

Pf_in[-1] = p[-1] * q_in(t[-1])
Pf_diss[-1] = R_h * qf[-1] ** 2
Ef[-1] = 0.5 * C_h * p[-1] ** 2 + 0.5 * L_h * qf[-1] ** 2

Wf_in = trapz(Pf_in, t)
Wf_diss = trapz(Pf_diss, t)
fluid_balance_error = Ef[-1] - Ef[0] - (Wf_in - Wf_diss)

# ==========================================================
# 4) Thermal: lumped thermal capacitance with thermal resistance to ambient
#    C_th dT/dt = Q_in(t) - (T - T_env)/R_th
#    Heat flow Qdot has unit W (already power)
# ==========================================================
C_th, R_th, T_env = 600.0, 0.4, 293.15

def Q_in(tt):
    return 120.0 + 40.0 * np.sin(0.15 * tt)

T = np.ones(n) * T_env
Pt_in = np.zeros(n)
Pt_out = np.zeros(n)
Et = np.zeros(n)  # stored thermal energy relative to ambient

for kidx in range(n - 1):
    qdot_in = Q_in(t[kidx])
    qdot_out = (T[kidx] - T_env) / R_th
    dT = (qdot_in - qdot_out) / C_th

    T[kidx + 1] = T[kidx] + dt * dT

    Pt_in[kidx] = qdot_in
    Pt_out[kidx] = qdot_out
    Et[kidx] = C_th * (T[kidx] - T_env)

Pt_in[-1] = Q_in(t[-1])
Pt_out[-1] = (T[-1] - T_env) / R_th
Et[-1] = C_th * (T[-1] - T_env)

Wt_in = trapz(Pt_in, t)
Wt_out = trapz(Pt_out, t)
thermal_balance_error = Et[-1] - Et[0] - (Wt_in - Wt_out)

print("=== Energy Balance Audit (explicit-Euler residuals (shrink as dt decreases)) ===")
print(f"Mechanical residual  : {mech_balance_error:.6e} J")
print(f"Electrical residual  : {elec_balance_error:.6e} J")
print(f"Fluid residual       : {fluid_balance_error:.6e} (energy units)")
print(f"Thermal residual     : {thermal_balance_error:.6e} J")
print()
print("Final stored energies:")
print(f"E_mech(tf) = {Em[-1]:.6f}")
print(f"E_elec(tf) = {Ee[-1]:.6f}")
print(f"E_fluid(tf)= {Ef[-1]:.6f}")
print(f"E_therm(tf)= {Et[-1]:.6f}")

      

This C++ version mirrors the Python logic using only the standard library. It is suitable for embedded-style numerical experiments and for comparing floating-point behavior across languages.

Code: Chapter18_Lesson1.cpp

// Chapter18_Lesson1.cpp
// Work, Energy, and Power in Mechanical, Electrical, Fluid, and Thermal Systems
// Standard C++17, no external dependencies.

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

static double trapz(const std::vector<double>& y, const std::vector<double>& t) {
    double s = 0.0;
    for (std::size_t k = 0; k + 1 < y.size(); ++k) {
        s += 0.5 * (y[k] + y[k + 1]) * (t[k + 1] - t[k]);
    }
    return s;
}

int main() {
    const double t0 = 0.0, tf = 12.0, dt = 1e-4;
    const int N = static_cast<int>((tf - t0) / dt) + 1;

    std::vector<double> t(N);
    for (int k = 0; k < N; ++k) t[k] = t0 + k * dt;

    // ==========================================================
    // 1) Mechanical: m xdd + c xd + k x = F(t)
    // ==========================================================
    const double m = 1.5, c = 0.8, ks = 12.0;
    auto Fin = [](double tt) { return 2.5 * std::sin(1.2 * tt) + 1.2 * std::cos(0.4 * tt); };

    std::vector<double> x(N, 0.0), v(N, 0.0), Pm_in(N, 0.0), Pm_diss(N, 0.0), Em(N, 0.0);

    for (int k = 0; k < N - 1; ++k) {
        const double f = Fin(t[k]);
        const double a = (f - c * v[k] - ks * x[k]) / m;
        v[k + 1] = v[k] + dt * a;
        x[k + 1] = x[k] + dt * v[k];

        Pm_in[k] = f * v[k];
        Pm_diss[k] = c * v[k] * v[k];
        Em[k] = 0.5 * m * v[k] * v[k] + 0.5 * ks * x[k] * x[k];
    }
    Pm_in[N - 1] = Fin(t[N - 1]) * v[N - 1];
    Pm_diss[N - 1] = c * v[N - 1] * v[N - 1];
    Em[N - 1] = 0.5 * m * v[N - 1] * v[N - 1] + 0.5 * ks * x[N - 1] * x[N - 1];

    const double mechResidual = Em.back() - Em.front() - (trapz(Pm_in, t) - trapz(Pm_diss, t));

    // ==========================================================
    // 2) Electrical: series RLC
    //    L didt + R i + q/C = v_in(t), dq/dt = i
    // ==========================================================
    const double R = 2.0, L = 0.6, C = 0.25;
    auto Vin = [](double tt) { return 5.0 * std::sin(2.0 * tt) + 2.0 * std::cos(0.5 * tt); };

    std::vector<double> q(N, 0.0), i(N, 0.0), Pe_in(N, 0.0), Pe_diss(N, 0.0), Ee(N, 0.0);

    for (int k = 0; k < N - 1; ++k) {
        const double vin = Vin(t[k]);
        const double di = (vin - R * i[k] - q[k] / C) / L;
        const double dq = i[k];

        i[k + 1] = i[k] + dt * di;
        q[k + 1] = q[k] + dt * dq;

        Pe_in[k] = vin * i[k];
        Pe_diss[k] = R * i[k] * i[k];
        Ee[k] = 0.5 * L * i[k] * i[k] + 0.5 * q[k] * q[k] / C;
    }
    Pe_in[N - 1] = Vin(t[N - 1]) * i[N - 1];
    Pe_diss[N - 1] = R * i[N - 1] * i[N - 1];
    Ee[N - 1] = 0.5 * L * i[N - 1] * i[N - 1] + 0.5 * q[N - 1] * q[N - 1] / C;

    const double elecResidual = Ee.back() - Ee.front() - (trapz(Pe_in, t) - trapz(Pe_diss, t));

    // ==========================================================
    // 3) Fluid: C_h dp/dt = q_in - q,   L_h dq/dt = p - R_h q
    // ==========================================================
    const double Ch = 0.08, Lh = 0.15, Rh = 1.7;
    auto qin = [](double tt) { return 0.8 * std::sin(0.9 * tt) + 0.2 * std::cos(0.3 * tt); };

    std::vector<double> p(N, 0.0), qf(N, 0.0), Pf_in(N, 0.0), Pf_diss(N, 0.0), Ef(N, 0.0);

    for (int k = 0; k < N - 1; ++k) {
        const double qsrc = qin(t[k]);
        const double dp = (qsrc - qf[k]) / Ch;
        const double dqf = (p[k] - Rh * qf[k]) / Lh;

        p[k + 1] = p[k] + dt * dp;
        qf[k + 1] = qf[k] + dt * dqf;

        Pf_in[k] = p[k] * qsrc;
        Pf_diss[k] = Rh * qf[k] * qf[k];
        Ef[k] = 0.5 * Ch * p[k] * p[k] + 0.5 * Lh * qf[k] * qf[k];
    }
    Pf_in[N - 1] = p[N - 1] * qin(t[N - 1]);
    Pf_diss[N - 1] = Rh * qf[N - 1] * qf[N - 1];
    Ef[N - 1] = 0.5 * Ch * p[N - 1] * p[N - 1] + 0.5 * Lh * qf[N - 1] * qf[N - 1];

    const double fluidResidual = Ef.back() - Ef.front() - (trapz(Pf_in, t) - trapz(Pf_diss, t));

    // ==========================================================
    // 4) Thermal: C_th dT/dt = Q_in - (T - T_env)/R_th
    // ==========================================================
    const double Cth = 600.0, Rth = 0.4, Tenv = 293.15;
    auto Qin = [](double tt) { return 120.0 + 40.0 * std::sin(0.15 * tt); };

    std::vector<double> T(N, Tenv), Pt_in(N, 0.0), Pt_out(N, 0.0), Et(N, 0.0);

    for (int k = 0; k < N - 1; ++k) {
        const double qdotIn = Qin(t[k]);
        const double qdotOut = (T[k] - Tenv) / Rth;
        const double dT = (qdotIn - qdotOut) / Cth;

        T[k + 1] = T[k] + dt * dT;

        Pt_in[k] = qdotIn;
        Pt_out[k] = qdotOut;
        Et[k] = Cth * (T[k] - Tenv);
    }
    Pt_in[N - 1] = Qin(t[N - 1]);
    Pt_out[N - 1] = (T[N - 1] - Tenv) / Rth;
    Et[N - 1] = Cth * (T[N - 1] - Tenv);

    const double thermResidual = Et.back() - Et.front() - (trapz(Pt_in, t) - trapz(Pt_out, t));

    std::cout << std::scientific << std::setprecision(6);
    std::cout << "=== Energy Balance Audit (C++) ===\n";
    std::cout << "Mechanical residual : " << mechResidual << "\n";
    std::cout << "Electrical residual : " << elecResidual << "\n";
    std::cout << "Fluid residual      : " << fluidResidual << "\n";
    std::cout << "Thermal residual    : " << thermResidual << "\n\n";

    std::cout << std::fixed << std::setprecision(6);
    std::cout << "Final stored energies:\n";
    std::cout << "E_mech(tf)  = " << Em.back() << "\n";
    std::cout << "E_elec(tf)  = " << Ee.back() << "\n";
    std::cout << "E_fluid(tf) = " << Ef.back() << "\n";
    std::cout << "E_therm(tf) = " << Et.back() << "\n";
    return 0;
}

      

The Java implementation uses arrays and explicit loops, which makes the numerical update and energy bookkeeping easy to inspect.

Code: Chapter18_Lesson1.java

// Chapter18_Lesson1.java
// Work, Energy, and Power in Mechanical, Electrical, Fluid, and Thermal Systems

public class Chapter18_Lesson1 {

    static double trapz(double[] y, double[] t) {
        double s = 0.0;
        for (int k = 0; k < y.length - 1; k++) {
            s += 0.5 * (y[k] + y[k + 1]) * (t[k + 1] - t[k]);
        }
        return s;
    }

    static double Fm(double tt) {
        return 2.5 * Math.sin(1.2 * tt) + 1.2 * Math.cos(0.4 * tt);
    }

    static double Vin(double tt) {
        return 5.0 * Math.sin(2.0 * tt) + 2.0 * Math.cos(0.5 * tt);
    }

    static double QinHyd(double tt) {
        return 0.8 * Math.sin(0.9 * tt) + 0.2 * Math.cos(0.3 * tt);
    }

    static double QdotTherm(double tt) {
        return 120.0 + 40.0 * Math.sin(0.15 * tt);
    }

    public static void main(String[] args) {
        double t0 = 0.0, tf = 12.0, dt = 1e-4;
        int N = (int) Math.round((tf - t0) / dt) + 1;

        double[] t = new double[N];
        for (int k = 0; k < N; k++) t[k] = t0 + k * dt;

        // ==========================================================
        // 1) Mechanical: m xdd + c xd + k x = F(t)
        // ==========================================================
        double m = 1.5, c = 0.8, ks = 12.0;
        double[] x = new double[N], v = new double[N];
        double[] PmIn = new double[N], PmDiss = new double[N], Em = new double[N];

        for (int k = 0; k < N - 1; k++) {
            double f = Fm(t[k]);
            double a = (f - c * v[k] - ks * x[k]) / m;

            v[k + 1] = v[k] + dt * a;
            x[k + 1] = x[k] + dt * v[k];

            PmIn[k] = f * v[k];
            PmDiss[k] = c * v[k] * v[k];
            Em[k] = 0.5 * m * v[k] * v[k] + 0.5 * ks * x[k] * x[k];
        }
        PmIn[N - 1] = Fm(t[N - 1]) * v[N - 1];
        PmDiss[N - 1] = c * v[N - 1] * v[N - 1];
        Em[N - 1] = 0.5 * m * v[N - 1] * v[N - 1] + 0.5 * ks * x[N - 1] * x[N - 1];
        double mechResidual = Em[N - 1] - Em[0] - (trapz(PmIn, t) - trapz(PmDiss, t));

        // ==========================================================
        // 2) Electrical: series RLC
        // ==========================================================
        double R = 2.0, L = 0.6, C = 0.25;
        double[] q = new double[N], i = new double[N];
        double[] PeIn = new double[N], PeDiss = new double[N], Ee = new double[N];

        for (int k = 0; k < N - 1; k++) {
            double vin = Vin(t[k]);
            double di = (vin - R * i[k] - q[k] / C) / L;
            double dq = i[k];

            i[k + 1] = i[k] + dt * di;
            q[k + 1] = q[k] + dt * dq;

            PeIn[k] = vin * i[k];
            PeDiss[k] = R * i[k] * i[k];
            Ee[k] = 0.5 * L * i[k] * i[k] + 0.5 * q[k] * q[k] / C;
        }
        PeIn[N - 1] = Vin(t[N - 1]) * i[N - 1];
        PeDiss[N - 1] = R * i[N - 1] * i[N - 1];
        Ee[N - 1] = 0.5 * L * i[N - 1] * i[N - 1] + 0.5 * q[N - 1] * q[N - 1] / C;
        double elecResidual = Ee[N - 1] - Ee[0] - (trapz(PeIn, t) - trapz(PeDiss, t));

        // ==========================================================
        // 3) Fluid: C_h dp/dt = q_in - q,   L_h dq/dt = p - R_h q
        // ==========================================================
        double Ch = 0.08, Lh = 0.15, Rh = 1.7;
        double[] p = new double[N], qf = new double[N];
        double[] PfIn = new double[N], PfDiss = new double[N], Ef = new double[N];

        for (int k = 0; k < N - 1; k++) {
            double qsrc = QinHyd(t[k]);
            double dp = (qsrc - qf[k]) / Ch;
            double dqf = (p[k] - Rh * qf[k]) / Lh;

            p[k + 1] = p[k] + dt * dp;
            qf[k + 1] = qf[k] + dt * dqf;

            PfIn[k] = p[k] * qsrc;
            PfDiss[k] = Rh * qf[k] * qf[k];
            Ef[k] = 0.5 * Ch * p[k] * p[k] + 0.5 * Lh * qf[k] * qf[k];
        }
        PfIn[N - 1] = p[N - 1] * QinHyd(t[N - 1]);
        PfDiss[N - 1] = Rh * qf[N - 1] * qf[N - 1];
        Ef[N - 1] = 0.5 * Ch * p[N - 1] * p[N - 1] + 0.5 * Lh * qf[N - 1] * qf[N - 1];
        double fluidResidual = Ef[N - 1] - Ef[0] - (trapz(PfIn, t) - trapz(PfDiss, t));

        // ==========================================================
        // 4) Thermal: C_th dT/dt = Q_in - (T - T_env)/R_th
        // ==========================================================
        double Cth = 600.0, Rth = 0.4, Tenv = 293.15;
        double[] T = new double[N], PtIn = new double[N], PtOut = new double[N], Et = new double[N];
        for (int k = 0; k < N; k++) T[k] = Tenv;

        for (int k = 0; k < N - 1; k++) {
            double qdotIn = QdotTherm(t[k]);
            double qdotOut = (T[k] - Tenv) / Rth;
            double dT = (qdotIn - qdotOut) / Cth;

            T[k + 1] = T[k] + dt * dT;

            PtIn[k] = qdotIn;
            PtOut[k] = qdotOut;
            Et[k] = Cth * (T[k] - Tenv);
        }
        PtIn[N - 1] = QdotTherm(t[N - 1]);
        PtOut[N - 1] = (T[N - 1] - Tenv) / Rth;
        Et[N - 1] = Cth * (T[N - 1] - Tenv);
        double thermResidual = Et[N - 1] - Et[0] - (trapz(PtIn, t) - trapz(PtOut, t));

        System.out.println("=== Energy Balance Audit (Java) ===");
        System.out.printf("Mechanical residual : %.6e%n", mechResidual);
        System.out.printf("Electrical residual : %.6e%n", elecResidual);
        System.out.printf("Fluid residual      : %.6e%n", fluidResidual);
        System.out.printf("Thermal residual    : %.6e%n%n", thermResidual);

        System.out.println("Final stored energies:");
        System.out.printf("E_mech(tf)  = %.6f%n", Em[N - 1]);
        System.out.printf("E_elec(tf)  = %.6f%n", Ee[N - 1]);
        System.out.printf("E_fluid(tf) = %.6f%n", Ef[N - 1]);
        System.out.printf("E_therm(tf) = %.6f%n", Et[N - 1]);
    }
}

      

The MATLAB script computes the same four-domain energy audits. The comments at the end describe the corresponding Simulink block layout (integrators, gains, sums, products, and power/work integrators), which students can assemble graphically.

Code: Chapter18_Lesson1.m

% Chapter18_Lesson1.m
% Work, Energy, and Power in Mechanical, Electrical, Fluid, and Thermal Systems
% MATLAB script (also suitable for translation to Simulink subsystems)

clear; clc;

t0 = 0; tf = 12; dt = 1e-4;
t = t0:dt:tf;
N = numel(t);

trapz_local = @(y, tt) trapz(tt, y);

%% 1) Mechanical: m*xdd + c*xd + k*x = F(t)
m = 1.5; c = 0.8; ks = 12.0;
Fin = @(tt) 2.5*sin(1.2*tt) + 1.2*cos(0.4*tt);

x = zeros(1,N); v = zeros(1,N);
Pm_in = zeros(1,N); Pm_diss = zeros(1,N); Em = zeros(1,N);

for k = 1:N-1
    f = Fin(t(k));
    a = (f - c*v(k) - ks*x(k))/m;
    v(k+1) = v(k) + dt*a;
    x(k+1) = x(k) + dt*v(k);

    Pm_in(k) = f*v(k);
    Pm_diss(k) = c*v(k)^2;
    Em(k) = 0.5*m*v(k)^2 + 0.5*ks*x(k)^2;
end
Pm_in(N) = Fin(t(N))*v(N);
Pm_diss(N) = c*v(N)^2;
Em(N) = 0.5*m*v(N)^2 + 0.5*ks*x(N)^2;
mechResidual = Em(end) - Em(1) - (trapz_local(Pm_in,t) - trapz_local(Pm_diss,t));

%% 2) Electrical: series RLC
R = 2.0; L = 0.6; C = 0.25;
Vin = @(tt) 5.0*sin(2.0*tt) + 2.0*cos(0.5*tt);

q = zeros(1,N); i = zeros(1,N);
Pe_in = zeros(1,N); Pe_diss = zeros(1,N); Ee = zeros(1,N);

for k = 1:N-1
    vin = Vin(t(k));
    di = (vin - R*i(k) - q(k)/C)/L;
    dq = i(k);

    i(k+1) = i(k) + dt*di;
    q(k+1) = q(k) + dt*dq;

    Pe_in(k) = vin*i(k);
    Pe_diss(k) = R*i(k)^2;
    Ee(k) = 0.5*L*i(k)^2 + 0.5*q(k)^2/C;
end
Pe_in(N) = Vin(t(N))*i(N);
Pe_diss(N) = R*i(N)^2;
Ee(N) = 0.5*L*i(N)^2 + 0.5*q(N)^2/C;
elecResidual = Ee(end) - Ee(1) - (trapz_local(Pe_in,t) - trapz_local(Pe_diss,t));

%% 3) Fluid (small-signal hydraulic chamber + line)
Ch = 0.08; Lh = 0.15; Rh = 1.7;
QinHyd = @(tt) 0.8*sin(0.9*tt) + 0.2*cos(0.3*tt);

p = zeros(1,N); qf = zeros(1,N);
Pf_in = zeros(1,N); Pf_diss = zeros(1,N); Ef = zeros(1,N);

for k = 1:N-1
    qsrc = QinHyd(t(k));
    dp = (qsrc - qf(k))/Ch;
    dqf = (p(k) - Rh*qf(k))/Lh;

    p(k+1) = p(k) + dt*dp;
    qf(k+1) = qf(k) + dt*dqf;

    Pf_in(k) = p(k)*qsrc;
    Pf_diss(k) = Rh*qf(k)^2;
    Ef(k) = 0.5*Ch*p(k)^2 + 0.5*Lh*qf(k)^2;
end
Pf_in(N) = p(N)*QinHyd(t(N));
Pf_diss(N) = Rh*qf(N)^2;
Ef(N) = 0.5*Ch*p(N)^2 + 0.5*Lh*qf(N)^2;
fluidResidual = Ef(end) - Ef(1) - (trapz_local(Pf_in,t) - trapz_local(Pf_diss,t));

%% 4) Thermal: Cth*dT/dt = Qdot_in - (T - Tenv)/Rth
Cth = 600; Rth = 0.4; Tenv = 293.15;
QdotTherm = @(tt) 120 + 40*sin(0.15*tt);

T = Tenv*ones(1,N);
Pt_in = zeros(1,N); Pt_out = zeros(1,N); Et = zeros(1,N);

for k = 1:N-1
    qdot_in = QdotTherm(t(k));
    qdot_out = (T(k) - Tenv)/Rth;
    dT = (qdot_in - qdot_out)/Cth;

    T(k+1) = T(k) + dt*dT;

    Pt_in(k) = qdot_in;
    Pt_out(k) = qdot_out;
    Et(k) = Cth*(T(k) - Tenv);
end
Pt_in(N) = QdotTherm(t(N));
Pt_out(N) = (T(N) - Tenv)/Rth;
Et(N) = Cth*(T(N) - Tenv);
thermResidual = Et(end) - Et(1) - (trapz_local(Pt_in,t) - trapz_local(Pt_out,t));

%% Report
fprintf('=== Energy Balance Audit (MATLAB) ===\n');
fprintf('Mechanical residual : %.6e\n', mechResidual);
fprintf('Electrical residual : %.6e\n', elecResidual);
fprintf('Fluid residual      : %.6e\n', fluidResidual);
fprintf('Thermal residual    : %.6e\n\n', thermResidual);

fprintf('Final stored energies:\n');
fprintf('E_mech(tf)  = %.6f\n', Em(end));
fprintf('E_elec(tf)  = %.6f\n', Ee(end));
fprintf('E_fluid(tf) = %.6f\n', Ef(end));
fprintf('E_therm(tf) = %.6f\n', Et(end));

%% Simulink note (manual build)
% Mechanical subsystem: Integrator(a->v), Integrator(v->x), Sum, Gain blocks
% Electrical subsystem: Integrator(di->i), Integrator(i->q), Sum, Gain blocks
% Fluid subsystem: Integrator(dp->p), Integrator(dq->qf), Sum, Gain blocks
% Thermal subsystem: Integrator(dT->T), Sum, Gain blocks
% Add Product blocks for power and Integrator blocks for work.

      

The following notebook text defines a valid \( .nb \) notebook structure with executable input cells. It uses NDSolveValue to solve the ODE models and then checks the energy balances.

Code: Chapter18_Lesson1.nb

Notebook[{
      Cell["Chapter18_Lesson1.nb", "Title"],
      Cell["Work, Energy, and Power in Mechanical, Electrical, Fluid, and Thermal Systems", "Subtitle"],
      Cell["This notebook simulates four lumped dynamic models and checks the energy balance in each domain.", "Text"],
      Cell["ClearAll[\"Global`*\"];\n\
t0 = 0; tf = 12; dt = 0.0001;\n\
tgrid = Range[t0, tf, dt];\n\
trapz[y_List, tt_List] := Total[(Most[y] + Rest[y])*(Differences[tt])/2];", "Input"],
      Cell["(* 1) Mechanical: m x'' + c x' + k x = F(t) *)\n\
m = 1.5; c = 0.8; ks = 12.0;\n\
Fin[t_] := 2.5 Sin[1.2 t] + 1.2 Cos[0.4 t];\n\
xFun = NDSolveValue[{m x''[t] + c x'[t] + ks x[t] == Fin[t], x[0] == 0, x'[0] == 0}, x, {t, t0, tf}];\n\
vFun[t_] := xFun'[t];\n\
PmIn = (Fin[#] vFun[#]) & /@ tgrid;\n\
PmDiss = (c vFun[#]^2) & /@ tgrid;\n\
Em = (0.5 m vFun[#]^2 + 0.5 ks xFun[#]^2) & /@ tgrid;\n\
mechResidual = Last[Em] - First[Em] - (trapz[PmIn, tgrid] - trapz[PmDiss, tgrid]);", "Input"],
      Cell["(* 2) Electrical: series RLC *)\n\
R = 2.0; L = 0.6; Cc = 0.25;\n\
Vin[t_] := 5.0 Sin[2.0 t] + 2.0 Cos[0.5 t];\n\
{qFun, iFun} = NDSolveValue[{\n\
   q'[t] == i[t],\n\
   L i'[t] + R i[t] + q[t]/Cc == Vin[t],\n\
   q[0] == 0, i[0] == 0}, {q, i}, {t, t0, tf}];\n\
PeIn = (Vin[#] iFun[#]) & /@ tgrid;\n\
PeDiss = (R iFun[#]^2) & /@ tgrid;\n\
Ee = (0.5 L iFun[#]^2 + 0.5 qFun[#]^2/Cc) & /@ tgrid;\n\
elecResidual = Last[Ee] - First[Ee] - (trapz[PeIn, tgrid] - trapz[PeDiss, tgrid]);", "Input"],
      Cell["(* 3) Fluid small-signal chamber + line *)\n\
Ch = 0.08; Lh = 0.15; Rh = 1.7;\n\
QinHyd[t_] := 0.8 Sin[0.9 t] + 0.2 Cos[0.3 t];\n\
{pFun, qfFun} = NDSolveValue[{\n\
   Ch p'[t] == QinHyd[t] - qf[t],\n\
   Lh qf'[t] == p[t] - Rh qf[t],\n\
   p[0] == 0, qf[0] == 0}, {p, qf}, {t, t0, tf}];\n\
PfIn = (pFun[#] QinHyd[#]) & /@ tgrid;\n\
PfDiss = (Rh qfFun[#]^2) & /@ tgrid;\n\
Ef = (0.5 Ch pFun[#]^2 + 0.5 Lh qfFun[#]^2) & /@ tgrid;\n\
fluidResidual = Last[Ef] - First[Ef] - (trapz[PfIn, tgrid] - trapz[PfDiss, tgrid]);", "Input"],
      Cell["(* 4) Thermal lumped capacitance *)\n\
Cth = 600.0; Rth = 0.4; Tenv = 293.15;\n\
QdotTherm[t_] := 120.0 + 40.0 Sin[0.15 t];\n\
TFun = NDSolveValue[{Cth T'[t] == QdotTherm[t] - (T[t] - Tenv)/Rth, T[0] == Tenv}, T, {t, t0, tf}];\n\
PtIn = (QdotTherm[#]) & /@ tgrid;\n\
PtOut = ((TFun[#] - Tenv)/Rth) & /@ tgrid;\n\
Et = (Cth (TFun[#] - Tenv)) & /@ tgrid;\n\
thermResidual = Last[Et] - First[Et] - (trapz[PtIn, tgrid] - trapz[PtOut, tgrid]);", "Input"],
      Cell["Print[\"=== Energy Balance Audit (Mathematica) ===\"];\n\
Print[\"Mechanical residual: \", NumberForm[mechResidual, {10, 6}]];\n\
Print[\"Electrical residual: \", NumberForm[elecResidual, {10, 6}]];\n\
Print[\"Fluid residual: \", NumberForm[fluidResidual, {10, 6}]];\n\
Print[\"Thermal residual: \", NumberForm[thermResidual, {10, 6}]];\n\
Print[\"Final stored energies:\"];\n\
Print[\"E_mech(tf) = \", NumberForm[Last[Em], {10, 6}]];\n\
Print[\"E_elec(tf) = \", NumberForm[Last[Ee], {10, 6}]];\n\
Print[\"E_fluid(tf)= \", NumberForm[Last[Ef], {10, 6}]];\n\
Print[\"E_therm(tf)= \", NumberForm[Last[Et], {10, 6}]];", "Input"]
    },
    WindowSize->{1000, 800},
    WindowTitle->"Chapter18_Lesson1"
    ]

      

Problem 1 (Mechanical power balance derivation): For the mass-spring-damper equation \( m\ddot{x} + c\dot{x} + kx = F_{in}(t) \), derive the energy balance equation and identify stored and dissipated terms.

Solution: Multiply the ODE by \( \dot{x} \):

\[ m\ddot{x}\dot{x} + c\dot{x}^2 + kx\dot{x} = F_{in}\dot{x}. \]

Use the identities \( m\ddot{x}\dot{x} = d(\tfrac{1}{2}m\dot{x}^2)/dt \) and \( kx\dot{x} = d(\tfrac{1}{2}kx^2)/dt \):

\[ \frac{d}{dt}\left(\tfrac{1}{2}m\dot{x}^2 + \tfrac{1}{2}kx^2\right) = F_{in}\dot{x} - c\dot{x}^2. \]

Therefore \( E_{stored} = \tfrac{1}{2}m\dot{x}^2 + \tfrac{1}{2}kx^2 \) and \( p_{diss} = c\dot{x}^2 \).

Problem 2 (Capacitor and inductor energies): Starting from \( i=C\dot{v} \) and \( v=L\dot{i} \), derive \( E_C = \tfrac{1}{2}Cv^2 \) and \( E_L = \tfrac{1}{2}Li^2 \).

Solution: For the capacitor,

\[ p_C = vi = vC\dot{v} = \frac{d}{dt}\left(\tfrac{1}{2}Cv^2\right). \]

Integrating in time yields \( E_C = \tfrac{1}{2}Cv^2 \) (assuming zero reference energy at \( v=0 \)). For the inductor,

\[ p_L = vi = Li\dot{i} = \frac{d}{dt}\left(\tfrac{1}{2}Li^2\right), \]

hence \( E_L = \tfrac{1}{2}Li^2 \).

Problem 3 (Hydraulic chamber-line balance): For \( C_h\dot{p} = q_{in} - q \) and \( L_h\dot{q} = p - R_h q \), prove the input power balance and identify the stored energy.

Solution: Multiply the first equation by \( p \) and the second by \( q \):

\[ C_h p\dot{p} = p q_{in} - pq, \qquad L_h q\dot{q} = pq - R_h q^2. \]

Adding the two equations cancels the internal exchange term \( pq \):

\[ \frac{d}{dt}\left(\tfrac{1}{2}C_h p^2 + \tfrac{1}{2}L_h q^2\right) = p q_{in} - R_h q^2. \]

Thus the storage is \( E = \tfrac{1}{2}C_h p^2 + \tfrac{1}{2}L_h q^2 \), the input power is \( p q_{in} \), and the dissipation is \( R_h q^2 \).

Problem 4 (Thermal RC transient): A thermal node has \( C_{th} = 500\,\mathrm{J/K} \), \( R_{th} = 0.5\,\mathrm{K/W} \), and ambient \( T_{env} \). It receives a constant heat input \( \dot{Q}_{in} = 100\,\mathrm{W} \). If \( T(0)=T_{env} \), find \( T(t) \) and the steady-state temperature rise.

Solution: The ODE is

\[ C_{th}\dot{T} = 100 - \frac{T - T_{env}}{R_{th}}. \]

Let \( \Delta T = T - T_{env} \). Then

\[ C_{th}\dot{\Delta T} + \frac{1}{R_{th}}\Delta T = 100. \]

The time constant is \( R_{th}C_{th} = 250\,\mathrm{s} \), and the steady-state rise is \( \Delta T_{ss} = 100R_{th} = 50\,\mathrm{K} \). Therefore

\[ T(t) = T_{env} + 50\left(1-e^{-t/250}\right). \]

Problem 5 (Cross-domain equivalence): Show that the mathematical structure of a mass-spring-damper and a series RLC circuit is identical under the mapping \( m \leftrightarrow L \), \( c \leftrightarrow R \), and \( k \leftrightarrow 1/C \).

Solution: The mechanical equation is \( m\ddot{x} + c\dot{x} + kx = F_{in} \). The series RLC in charge coordinates is \( L\ddot{q} + R\dot{q} + q/C = v_s \). Under the variable correspondence \( x \leftrightarrow q \) and source correspondence \( F_{in} \leftrightarrow v_s \), the ODEs are identical. Their energy functions are likewise identical in form:

\[ E_{mech} = \tfrac{1}{2}m\dot{x}^2 + \tfrac{1}{2}kx^2 \quad \leftrightarrow \quad E_{elec} = \tfrac{1}{2}L i^2 + \tfrac{1}{2}\frac{q^2}{C}. \]

We established a unified work-energy-power framework across four major engineering domains and derived storage/dissipation laws for canonical lumped elements. The key recurring pattern is a domain-independent power balance: source power equals the rate of stored energy plus losses/transfers. This viewpoint is foundational for Lagrangian, Hamiltonian, and bond-graph modeling in the next lessons.

Practically, students should now be able to (i) identify conjugate variables, (ii) derive element energy functions, (iii) write a power-balance equation for a lumped model, and (iv) use energy auditing as a simulation validation test.