We work with the effort/flow variables of the hydraulic domain: pressure \( p(t) \) (Pa) and volumetric flow rate \( q(t) \) (m³/s). Across a two-terminal element we define the pressure drop \( \Delta p(t) \) and the through-variable \( q(t) \) with a consistent sign convention.

The three lumped elements are defined by linear constitutive relations: hydraulic resistance \( \Delta p = R_h q \), hydraulic inertance \( \Delta p = L_h \frac{dq}{dt} \), and hydraulic capacitance \( q = C_h \frac{dp}{dt} \) (or equivalently \( p = \frac{1}{C_h}\int q\,dt \)). These mirror the mathematical roles of R, L, C in electrical circuits (Lesson 1), while the analogy choices were formalized in Lesson 2.

flowchart TD
  A["Given: geometry, fluid properties, operating regime"] --> B["Choose variables: pressure p(t), flow q(t)"]
  B --> C["Decide lumped elements: R_h, C_h, L_h"]
  C --> D["Write element laws: dp = R*q ; dp = L*dq/dt ; q = C*dp/dt"]
  D --> E["Apply continuity at junctions: sum(q)=0"]
  E --> F["Apply pressure compatibility in loops: sum(dp)=0"]
  F --> G["Obtain ODE model"]
  G --> H["Simulate and sanity-check units/limits"]
        

Scope note. We emphasize linear lumped models consistent with laminar flow and small perturbations. Strong turbulence typically introduces nonlinear pressure-loss laws (e.g., proportional to \( q|q| \)), which we only mention qualitatively here because systematic linearization is treated later (Chapter 9).

Consider a two-terminal hydraulic element with inlet pressure \( p_1(t) \) and outlet pressure \( p_2(t) \). Define the pressure drop \( \Delta p(t) = p_1(t) - p_2(t) \) and take \( q(t) \) positive from port 1 to port 2.

The instantaneous power delivered to the element is \( \mathcal{P}(t) = \Delta p(t)\,q(t) \), with units Pa·m³/s = W. This power identity is the backbone for energy and passivity proofs below.

\[ [\Delta p] = \mathrm{Pa},\quad [q] = \mathrm{m}^3/\mathrm{s},\quad [\Delta p\,q] = \mathrm{Pa}\cdot \mathrm{m}^3/\mathrm{s} = \mathrm{N}\cdot\mathrm{m}/\mathrm{s}=\mathrm{W}. \]

A constitutive law is physically admissible if it does not violate basic thermodynamics; in particular, a pure resistance must dissipate power (\( \mathcal{P}(t)\ge 0 \) for all \( t \) under the adopted sign convention), while inertance and capacitance must store recoverable energy.

For steady, fully developed, incompressible, laminar flow of a Newtonian fluid in a circular pipe (radius \( a \), length \( \ell \)), the Hagen–Poiseuille law yields a linear relation between pressure drop and flow rate: \( \Delta p = R_h q \).

3.1 Assumptions and reduced momentum balance

In cylindrical coordinates with axial velocity \( v(r) \) and constant axial pressure gradient \( dp/dx \), the axial Navier–Stokes equation reduces to \( 0 = -\frac{dp}{dx} + \mu\left(\frac{1}{r}\frac{d}{dr}\left(r\frac{dv}{dr}\right)\right) \), where \( \mu \) is dynamic viscosity.

\[ 0 = -\frac{dp}{dx} + \mu\left(\frac{1}{r}\frac{d}{dr}\left(r\frac{dv}{dr}\right)\right), \quad v(a)=0,\quad \left|\frac{dv}{dr}\right| < \infty \text{ at } r=0. \]

3.2 Proof of Poiseuille law (outline with key steps)

Integrate twice. First multiply by \( r \) and integrate:

\[ \frac{d}{dr}\left(r\frac{dv}{dr}\right) = \frac{r}{\mu}\frac{dp}{dx} \;\Rightarrow\; r\frac{dv}{dr} = \frac{1}{2\mu}\frac{dp}{dx}r^2 + c_1. \]

Finite slope at \( r=0 \) forces \( c_1=0 \), so \( \frac{dv}{dr} = \frac{1}{2\mu}\frac{dp}{dx}r \). Integrating again:

\[ v(r) = \frac{1}{4\mu}\frac{dp}{dx}r^2 + c_2. \]

Apply no-slip \( v(a)=0 \) to obtain \( c_2 = -\frac{1}{4\mu}\frac{dp}{dx}a^2 \), hence

\[ v(r) = \frac{1}{4\mu}\frac{dp}{dx}\left(r^2-a^2\right). \]

The volumetric flow rate is \( q=\int_A v\,dA = 2\pi\int_0^a v(r)\,r\,dr \). Substituting and integrating gives

\[ q = 2\pi\int_0^a \frac{1}{4\mu}\frac{dp}{dx}\left(r^2-a^2\right)r\,dr = -\frac{\pi a^4}{8\mu}\frac{dp}{dx}. \]

With \( \Delta p = p_1-p_2 = -\frac{dp}{dx}\ell \), we conclude

\[ \Delta p = \left(\frac{8\mu\ell}{\pi a^4}\right) q \;\;\Rightarrow\;\; R_h = \frac{8\mu\ell}{\pi a^4}. \]

3.3 Dissipation and passivity of \( R_h \)

For a linear resistance \( \Delta p = R_h q \) with \( R_h > 0 \), the power absorbed is

\[ \mathcal{P}(t)=\Delta p(t)\,q(t)=R_h q^2(t)\ge 0. \]

Thus the resistance is dissipative and passive: it cannot generate net energy.

A hydraulic capacitance captures the ability of a volume to store fluid (or store pressure) such that flow is proportional to the time-rate of change of pressure: \( q = C_h \frac{dp}{dt} \).

4.1 Compressible chamber: \( C_h = \frac{V}{\beta} \)

Let a control volume of nominal fluid volume \( V \) be subject to uniform pressure \( p(t) \). Define the bulk modulus \( \beta \) by

\[ \beta \equiv -V\frac{dp}{dV}. \]

For small perturbations about an operating point (lumped approximation), linearize:

\[ dp = -\frac{\beta}{V}\,dV \;\Rightarrow\; \frac{dV}{dt} = -\frac{V}{\beta}\frac{dp}{dt}. \]

If we define \( q \) as the net inflow that increases stored volume in the chamber, then \( q = \frac{dV}{dt} \) (sign depends on the chosen convention). Taking the storage sign so that increasing pressure corresponds to net inflow, we write the standard capacitance form:

\[ q = C_h \frac{dp}{dt}, \quad C_h = \frac{V}{\beta}. \]

4.2 Free-surface tank: \( C_h = \frac{A}{\rho g} \)

For a vertical tank of constant cross-sectional area \( A \) open to atmosphere at the top, with fluid height \( h(t) \), the gauge pressure at the bottom is \( p(t)=\rho g h(t) \). The stored volume is \( V(t)=A h(t) \), so \( q = dV/dt = A\,dh/dt \). Therefore

\[ q = A\frac{dh}{dt} = A\frac{1}{\rho g}\frac{dp}{dt} \;\Rightarrow\; C_h = \frac{A}{\rho g}. \]

4.3 Energy stored in \( C_h \)

For a pure capacitance \( q = C_h \frac{dp}{dt} \), the instantaneous power is \( \mathcal{P} = p q = p C_h \frac{dp}{dt} \). Integrating shows stored energy:

\[ \mathcal{E}_C(t) = \int_0^t p(\tau)\,q(\tau)\,d\tau = \int_0^t p(\tau) C_h \frac{dp(\tau)}{d\tau}\,d\tau = \frac{1}{2}C_h p^2(t) - \frac{1}{2}C_h p^2(0). \]

Hence a capacitance stores recoverable energy and is lossless in the ideal model.

Hydraulic inertance models the pressure required to accelerate a slug of fluid. For a uniform pipe segment of length \( \ell \) and cross-sectional area \( A \), with fluid density \( \rho \), we derive the relation \( \Delta p = L_h \frac{dq}{dt} \).

5.1 Derivation from Newton’s second law (quasi-1D)

The net force on the fluid slug is pressure force: \( F = \Delta p \, A \). The slug mass is \( m=\rho A\ell \), and mean velocity is \( v=q/A \), so acceleration is \( dv/dt = \frac{1}{A}\frac{dq}{dt} \). Newton’s law gives

\[ \Delta p\,A = m\frac{dv}{dt} = (\rho A\ell)\left(\frac{1}{A}\frac{dq}{dt}\right) \;\Rightarrow\; \Delta p = \left(\frac{\rho\ell}{A}\right)\frac{dq}{dt}. \]

Therefore the hydraulic inertance is

\[ L_h = \frac{\rho\ell}{A}. \]

5.2 Energy stored in \( L_h \)

For \( \Delta p = L_h \frac{dq}{dt} \), the power is \( \mathcal{P}=\Delta p\,q = L_h q \frac{dq}{dt} \). Integrating yields

\[ \mathcal{E}_L(t) = \int_0^t L_h q(\tau)\frac{dq(\tau)}{d\tau}\,d\tau = \frac{1}{2}L_h q^2(t) - \frac{1}{2}L_h q^2(0). \]

Thus inertance stores kinetic energy and is lossless in the ideal model.

For lumped hydraulic networks, the governing equations are built from:

• Continuity (junction law): at a node (common pressure), the signed sum of flows is zero.
• Pressure compatibility (loop law): around a closed loop, the signed sum of pressure drops is zero.

6.1 Junction law (continuity)

If \( q_1,\dots,q_n \) are flows incident to a node, with signs chosen positive into the node, then

\[ \sum_{k=1}^n q_k(t)=0. \]

6.2 Loop law (pressure compatibility)

For elements in a series path, pressure drops add. For a loop with elements \( e_1,\dots,e_m \),

\[ \sum_{j=1}^m \Delta p_{e_j}(t)=0. \]

6.3 Dimensional check as a modeling guardrail

It is useful to verify that each element has consistent units:

\[ [R_h] = \frac{\mathrm{Pa}}{\mathrm{m}^3/\mathrm{s}}=\frac{\mathrm{Pa}\cdot \mathrm{s}}{\mathrm{m}^3},\quad [C_h] = \frac{\mathrm{m}^3/\mathrm{s}}{\mathrm{Pa}/\mathrm{s}}=\frac{\mathrm{m}^3}{\mathrm{Pa}},\quad [L_h] = \frac{\mathrm{Pa}}{(\mathrm{m}^3/\mathrm{s}^2)}=\frac{\mathrm{Pa}\cdot \mathrm{s}^2}{\mathrm{m}^3}. \]

Consider a pressure source \( p_s(t) \) driving a series connection of a resistance and inertance into a compliant volume modeled as a capacitance to a reference pressure (ground). Let the node pressure across the capacitance be \( p(t) \), and the series flow be \( q(t) \).

flowchart TD
  PS["Source: p_s(t)"] --> R["R_h: dp = R*q"]
  R --> L["L_h: dp = L*dq/dt"]
  L --> N["Node pressure p(t)"]
  N --> C["C_h: q = C*dp/dt"]
  C --> G["Reference pressure (0 gauge)"]
        

7.1 Governing equations

Series path implies the same flow \( q(t) \) through \( R_h \) and \( L_h \), and into the capacitance:

\[ q(t) = C_h\frac{dp(t)}{dt}. \]

Pressure compatibility along the series path gives:

\[ p_s(t) - p(t) = \Delta p_R(t) + \Delta p_L(t) = R_h q(t) + L_h\frac{dq(t)}{dt}. \]

Substitute \( q = C_h \frac{dp}{dt} \) and \( \frac{dq}{dt} = C_h \frac{d^2 p}{dt^2} \):

\[ p_s(t) - p(t) = R_h C_h\frac{dp(t)}{dt} + L_h C_h\frac{d^2 p(t)}{dt^2}. \]

Rearranging yields the standard second-order ODE in pressure:

\[ L_h C_h\frac{d^2 p(t)}{dt^2} + R_h C_h\frac{dp(t)}{dt} + p(t) = p_s(t). \]

7.2 Natural frequency and damping ratio (parameter interpretation)

Divide by \( L_h C_h \) and compare to the canonical second-order form:

\[ \frac{d^2 p}{dt^2} + \frac{R_h}{L_h}\frac{dp}{dt} + \frac{1}{L_h C_h}p = \frac{1}{L_h C_h}p_s(t). \]

Define the (undamped) natural frequency \( \omega_n \) and damping ratio \( \zeta \):

\[ \omega_n = \frac{1}{\sqrt{L_h C_h}},\quad 2\zeta\omega_n = \frac{R_h}{L_h} \;\Rightarrow\; \zeta = \frac{R_h}{2}\sqrt{\frac{C_h}{L_h}}. \]

These expressions connect physical parameters to time-domain behavior: increasing \( L_h \) or \( C_h \) lowers \( \omega_n \), while increasing \( R_h \) increases damping.

We simulate the canonical ODE \( L_h C_h p'' + R_h C_h p' + p = p_s(t) \) by rewriting it as a first-order system in \( x_1=p \), \( x_2=\frac{dp}{dt} \):

\[ \dot{x}_1 = x_2,\quad \dot{x}_2 = \frac{1}{L_h C_h}\Big(p_s(t) - x_1 - R_h C_h x_2\Big). \]

8.1 Python (SciPy)

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# Parameters
R_h = 2.0e7       # Pa*s/m^3
L_h = 5.0e5       # Pa*s^2/m^3
C_h = 2.0e-10     # m^3/Pa

def p_source(t):
    # Step in source pressure (Pa)
    return 1.0e5 if t >= 0.05 else 0.0

def rhs(t, x):
    p = x[0]
    pdot = x[1]
    p_s = p_source(t)
    pddot = (p_s - p - R_h * C_h * pdot) / (L_h * C_h)
    return [pdot, pddot]

t0, tf = 0.0, 1.0
x0 = [0.0, 0.0]
t_eval = np.linspace(t0, tf, 2000)

sol = solve_ivp(rhs, (t0, tf), x0, t_eval=t_eval, rtol=1e-8, atol=1e-10)

p = sol.y[0, :]
pdot = sol.y[1, :]
q = C_h * pdot  # flow into capacitance

plt.figure()
plt.plot(sol.t, p)
plt.xlabel("t (s)")
plt.ylabel("p(t) (Pa)")
plt.title("Node pressure response")
plt.grid(True)

plt.figure()
plt.plot(sol.t, q)
plt.xlabel("t (s)")
plt.ylabel("q(t) (m^3/s)")
plt.title("Flow response (q = C_h * dp/dt)")
plt.grid(True)
plt.show()
      

Useful libraries for subsequent expansions: sympy (symbolic manipulation), scipy (ODEs), and (later, when transfer-function tools are introduced in Chapter 6) the Python control package.

8.2 C++ (from-scratch RK4 + CSV output)

#include <iostream>
#include <fstream>
#include <cmath>

static double R_h = 2.0e7;     // Pa*s/m^3
static double L_h = 5.0e5;     // Pa*s^2/m^3
static double C_h = 2.0e-10;   // m^3/Pa

double p_source(double t) {
    return (t >= 0.05) ? 1.0e5 : 0.0;
}

// x[0]=p, x[1]=pdot
void f(double t, const double x[2], double dxdt[2]) {
    double p = x[0];
    double pdot = x[1];
    double ps = p_source(t);
    dxdt[0] = pdot;
    dxdt[1] = (ps - p - R_h * C_h * pdot) / (L_h * C_h);
}

int main() {
    double t0 = 0.0, tf = 1.0, h = 1e-4;
    double x[2] = {0.0, 0.0};

    std::ofstream out("hydraulic_rlc.csv");
    out << "t,p,q\n";

    for (double t = t0; t <= tf; t += h) {
        // output
        double q = C_h * x[1];
        out << t << "," << x[0] << "," << q << "\n";

        // RK4
        double k1[2], k2[2], k3[2], k4[2], xtmp[2];

        f(t, x, k1);

        xtmp[0] = x[0] + 0.5*h*k1[0];
        xtmp[1] = x[1] + 0.5*h*k1[1];
        f(t + 0.5*h, xtmp, k2);

        xtmp[0] = x[0] + 0.5*h*k2[0];
        xtmp[1] = x[1] + 0.5*h*k2[1];
        f(t + 0.5*h, xtmp, k3);

        xtmp[0] = x[0] + h*k3[0];
        xtmp[1] = x[1] + h*k3[1];
        f(t + h, xtmp, k4);

        x[0] += (h/6.0)*(k1[0] + 2.0*k2[0] + 2.0*k3[0] + k4[0]);
        x[1] += (h/6.0)*(k1[1] + 2.0*k2[1] + 2.0*k3[1] + k4[1]);
    }

    out.close();
    std::cout << "Wrote hydraulic_rlc.csv\n";
    return 0;
}
      

For larger projects, consider Boost.Odeint (ODE integration) and Eigen (linear algebra) as foundational C++ libraries for dynamic-system simulation infrastructure.

8.3 Java (Apache Commons Math ODE integrator)

import java.io.FileWriter;
import java.io.IOException;
import org.apache.commons.math3.ode.FirstOrderDifferentialEquations;
import org.apache.commons.math3.ode.nonstiff.DormandPrince54Integrator;

public class HydraulicRLC {

    static double R_h = 2.0e7;    // Pa*s/m^3
    static double L_h = 5.0e5;    // Pa*s^2/m^3
    static double C_h = 2.0e-10;  // m^3/Pa

    static double pSource(double t) {
        return (t >= 0.05) ? 1.0e5 : 0.0;
    }

    // x[0]=p, x[1]=pdot
    static class Model implements FirstOrderDifferentialEquations {
        public int getDimension() { return 2; }

        public void computeDerivatives(double t, double[] x, double[] xDot) {
            double p = x[0];
            double pdot = x[1];
            double ps = pSource(t);
            xDot[0] = pdot;
            xDot[1] = (ps - p - R_h * C_h * pdot) / (L_h * C_h);
        }
    }

    public static void main(String[] args) throws IOException {
        double t0 = 0.0, tf = 1.0;
        double[] x = new double[] {0.0, 0.0};

        DormandPrince54Integrator integrator =
            new DormandPrince54Integrator(1e-6, 1e-2, 1e-10, 1e-8);

        FileWriter out = new FileWriter("hydraulic_rlc_java.csv");
        out.write("t,p,q\n");

        // simple fixed sampling loop: integrate to each sample time
        int N = 2000;
        for (int i = 0; i <= N; i++) {
            double t = t0 + (tf - t0) * i / (double) N;
            integrator.integrate(new Model(), (i==0 ? t0 : (t0 + (tf - t0) * (i-1) / (double) N)), x, t, x);

            double p = x[0];
            double q = C_h * x[1];
            out.write(t + "," + p + "," + q + "\n");
        }

        out.close();
        System.out.println("Wrote hydraulic_rlc_java.csv");
    }
}
      

Apache Commons Math provides robust ODE solvers; for plotting, export CSV and visualize with your preferred tool.

8.4 MATLAB (ode45)

% Parameters
R_h = 2.0e7;      % Pa*s/m^3
L_h = 5.0e5;      % Pa*s^2/m^3
C_h = 2.0e-10;    % m^3/Pa

p_source = @(t) (t >= 0.05) * 1.0e5;

% x(1)=p, x(2)=pdot
rhs = @(t,x) [ x(2);
               (p_source(t) - x(1) - R_h*C_h*x(2)) / (L_h*C_h) ];

tspan = [0 1];
x0 = [0; 0];

opts = odeset('RelTol',1e-8,'AbsTol',1e-10);
[t,x] = ode45(rhs, tspan, x0, opts);

p = x(:,1);
q = C_h * x(:,2);

figure; plot(t,p); grid on; xlabel('t (s)'); ylabel('p(t) (Pa)');
title('Node pressure response');

figure; plot(t,q); grid on; xlabel('t (s)'); ylabel('q(t) (m^3/s)');
title('Flow response (q = C_h * dp/dt)');
      

8.5 Simulink (block-level implementation)

Implement the first-order system \( \dot{x}_1=x_2 \), \( \dot{x}_2=\frac{1}{L_h C_h}(p_s - x_1 - R_h C_h x_2) \) using two Integrator blocks in cascade:

• Integrator 1 output: \( x_2 \) (input is \( \dot{x}_2 \))
• Integrator 2 output: \( x_1=p \) (input is \( \dot{x}_1=x_2 \))
• Compute \( \dot{x}_2 \) using Sum and Gain blocks implementing \( (p_s - x_1 - R_h C_h x_2)/(L_h C_h) \)
• Compute flow as \( q = C_h x_2 \) using a Gain block

The following MATLAB script programmatically creates a minimal Simulink model with these blocks:

model = 'hydraulic_RLC';
new_system(model); open_system(model);

% Add blocks
add_block('simulink/Sources/Step', [model '/Step']);
add_block('simulink/Math Operations/Sum', [model '/Sum']);
add_block('simulink/Math Operations/Gain', [model '/Gain_RCh']);
add_block('simulink/Math Operations/Gain', [model '/Gain_invLCh']);
add_block('simulink/Continuous/Integrator', [model '/Int_pdot']); % outputs x2
add_block('simulink/Continuous/Integrator', [model '/Int_p']);    % outputs x1
add_block('simulink/Math Operations/Gain', [model '/Gain_C']);    % q = C*x2
add_block('simulink/Sinks/Scope', [model '/Scope_p']);
add_block('simulink/Sinks/Scope', [model '/Scope_q']);

% Parameters (set in workspace)
R_h = 2.0e7; L_h = 5.0e5; C_h = 2.0e-10;

% Configure gains
set_param([model '/Gain_RCh'], 'Gain', 'R_h*C_h');
set_param([model '/Gain_invLCh'], 'Gain', '1/(L_h*C_h)');
set_param([model '/Gain_C'], 'Gain', 'C_h');

% Sum block signs: ps - p - RCh*pdot
set_param([model '/Sum'], 'Inputs', '+++'); % we will wire negative via gain and sum structure

% Wire: Step -> Sum(+)
add_line(model, 'Step/1', 'Sum/1');

% Wire: p (x1) feedback with negative sign using an extra Gain(-1) via Gain_invLCh placement
% For simplicity, directly place a Gain(-1) block if you want explicit negatives.

% Minimal wiring approach (recommended to do manually for clarity):
% Step feeds Sum; p and RCh*pdot feed Sum with negative signs; Sum feeds Gain_invLCh -> Int_pdot -> Int_p.

save_system(model);
      

In practice, students typically build this diagram manually to internalize the signal flow, then use scripts for parameter sweeps and reproducibility.

8.6 Wolfram Mathematica (symbolic + numerical)

(* Parameters *)
Rh = 2.0*^7;   (* Pa*s/m^3 *)
Lh = 5.0*^5;   (* Pa*s^2/m^3 *)
Ch = 2.0*^-10; (* m^3/Pa  *)

ps[t_] := Piecewise[{ {1.0*^5, t >= 0.05} }, 0.0];

(* ODE: Lh Ch p'' + Rh Ch p' + p = ps(t) *)
ode = Lh*Ch*p''[t] + Rh*Ch*p'[t] + p[t] == ps[t];
ics = {p[0] == 0, p'[0] == 0};

sol = NDSolve[{ode, ics}, p, {t, 0, 1}, AccuracyGoal -> 12, PrecisionGoal -> 10][[1]];

pPlot = Plot[Evaluate[p[t] /. sol], {t, 0, 1}, PlotGridLines -> Automatic,
  AxesLabel -> {"t (s)", "p(t) (Pa)"}, PlotLabel -> "Node pressure response"];

q[t_] := Ch*(p'[t] /. sol);
qPlot = Plot[Evaluate[q[t]], {t, 0, 1}, PlotGridLines -> Automatic,
  AxesLabel -> {"t (s)", "q(t) (m^3/s)"}, PlotLabel -> "Flow response"];

{pPlot, qPlot}
      

Problem 1 (Derive \( R_h \) for a circular pipe): Starting from the reduced laminar-flow equation \( 0 = -\frac{dp}{dx} + \mu\left(\frac{1}{r}\frac{d}{dr}\left(r\frac{dv}{dr}\right)\right) \), derive \( \Delta p = R_h q \) and obtain \( R_h \) in terms of \( \mu,\ell,a \).

Solution: From Section 3, integrating twice and applying boundary conditions gives

\[ v(r)=\frac{1}{4\mu}\frac{dp}{dx}(r^2-a^2),\quad q=2\pi\int_0^a v(r)r\,dr=-\frac{\pi a^4}{8\mu}\frac{dp}{dx}. \]

With \( \Delta p = -\frac{dp}{dx}\ell \), we obtain

\[ \Delta p = \left(\frac{8\mu\ell}{\pi a^4}\right) q, \quad R_h = \frac{8\mu\ell}{\pi a^4}. \]

Problem 2 (Capacitance from compressibility): A sealed chamber of volume \( V=5\times 10^{-4}\,\mathrm{m}^3 \) contains a fluid with bulk modulus \( \beta = 1.6\times 10^9\,\mathrm{Pa} \). Compute \( C_h \). If the pressure changes at \( dp/dt = 2\times 10^6\,\mathrm{Pa/s} \), what net inflow \( q \) is required?

Solution:

\[ C_h=\frac{V}{\beta}=\frac{5\times 10^{-4}}{1.6\times 10^9} = 3.125\times 10^{-13}\,\mathrm{m}^3/\mathrm{Pa}. \]

Using \( q=C_h \frac{dp}{dt} \),

\[ q = (3.125\times 10^{-13})(2\times 10^6) = 6.25\times 10^{-7}\,\mathrm{m}^3/\mathrm{s}. \]

Problem 3 (Inertance and kinetic energy): A straight pipe segment has length \( \ell=2\,\mathrm{m} \) and area \( A=5\times 10^{-4}\,\mathrm{m}^2 \). The fluid density is \( \rho=1000\,\mathrm{kg/m}^3 \). (a) Compute \( L_h \). (b) If the flow is \( q=2\times 10^{-4}\,\mathrm{m}^3/\mathrm{s} \), compute stored kinetic energy \( \mathcal{E}_L=\tfrac{1)2}L_h q^2 \).

Solution:

\[ L_h=\frac{\rho\ell}{A}=\frac{(1000)(2)}{5\times 10^{-4}}=4\times 10^6\;\mathrm{Pa\cdot s^2/m^3}. \]

\[ \mathcal{E}_L=\frac{1}{2}L_h q^2 =\frac{1}{2}(4\times 10^6)(2\times 10^{-4})^2 =\frac{1}{2}(4\times 10^6)(4\times 10^{-8}) = 0.08\;\mathrm{J}. \]

Problem 4 (Derive the pressure ODE for the \( R_h \)–\( L_h \)–\( C_h \) network): Using \( q=C_h\,dp/dt \) and \( p_s-p=R_h q + L_h\,dq/dt \), derive the second-order ODE in \( p(t) \). Then identify \( \omega_n \) and \( \zeta \) in terms of \( R_h, L_h, C_h \).

Solution:

Substitute \( q=C_h p' \) and \( q' = C_h p'' \) into \( p_s-p=R_h q + L_h q' \):

\[ p_s(t)-p(t)=R_h C_h p'(t)+L_h C_h p''(t) \;\Rightarrow\; L_h C_h p'' + R_h C_h p' + p = p_s(t). \]

Divide by \( L_h C_h \):

\[ p'' + \frac{R_h}{L_h}p' + \frac{1}{L_h C_h}p = \frac{1}{L_h C_h}p_s(t). \]

Comparing with \( p'' + 2\zeta\omega_n p' + \omega_n^2 p = \omega_n^2 p_s(t) \) yields

\[ \omega_n = \frac{1}{\sqrt{L_h C_h}},\quad \zeta = \frac{R_h}{2}\sqrt{\frac{C_h}{L_h}}. \]

Problem 5 (Consistency checks via limiting cases): Consider the ODE \( L_h C_h p'' + R_h C_h p' + p = p_s(t) \). (a) Show that if \( L_h \rightarrow 0 \), the model reduces to a first-order equation. (b) Show that if \( C_h \rightarrow 0 \) with bounded signals, the node pressure tracks the source pressure (in the lumped limit).

Solution:

(a) Let \( L_h \rightarrow 0 \). The term \( L_h C_h p'' \) vanishes, giving

\[ R_h C_h p'(t) + p(t) = p_s(t), \]

which is first-order with time constant \( R_h C_h \). This corresponds to removing inertia so pressure responds without momentum lag.

(b) Let \( C_h \rightarrow 0 \). Then \( q=C_h p' \rightarrow 0 \) for bounded \( p' \), meaning negligible flow into storage. The series path must then have negligible flow, implying negligible pressure drops across \( R_h \) and \( L_h \) in the lumped linear model; hence \( p(t) \approx p_s(t) \). Formally, with \( C_h \rightarrow 0 \), the terms \( L_h C_h p'' \) and \( R_h C_h p' \) vanish and the ODE reduces to \( p=p_s \).

We defined hydraulic effort/flow variables and constructed lumped models using three primitive elements: resistance \( \Delta p=R_h q \) (derived via Hagen–Poiseuille under laminar assumptions), capacitance \( q=C_h dp/dt \) (from compressibility \( C_h=V/\beta \) or free-surface storage \( C_h=A/(\rho g) \)), and inertance \( \Delta p=L_h dq/dt \) (from momentum balance \( L_h=\rho\ell/A \)). Using continuity and pressure-compatibility laws, we assembled networks into ODEs and analyzed a canonical \( R_h \)–\( L_h \)–\( C_h \) example whose behavior is governed by \( \omega_n=1/\sqrt{L_h C_h} \) and \( \zeta=(R_h/2)\sqrt{C_h/L_h} \). We then implemented the model in Python, C++, Java, MATLAB/Simulink, and Mathematica.