Developing high-fidelity models of hydraulic systems

Why focus on hydraulics? The answer is essentially hydraulic modelling is really hard (in numerical computing terms, hydraulic models are often referred to as "stiff" ODE's, which require more rigorous solvers from standard ODE's). Solving the challenges of modeling hydraulics is applicable to the numerical modeling challenges of all other domains. Let's first start with the concept of compressibility. Often we think of a liquid as incompressible, imagine attempting to "squeeze" water, it can be done but takes some very high forces. Therefore, if the model in question won't be solving a problem with high forces, it can be assumed incompressible. However, most hydrulic industrial models will involve high forces, this is precisely the area where most hydraulic machines are used.

Compressibility

Density

Density is simply mass over volume

\[\rho = m/V\]

Given a volume and mass of liquid, if the volume were to change from $V_0$ to $V$, we know that the pressure would increase, and since the mass in this case was constant, the density will increase as well.

The change in pressure for an isothermal compressible process is typically given as

\[\Delta p = -\beta \frac{\Delta V}{V_0}\]

Calculating Density as a Function of Pressure

Substituting $\Delta p$ and $\Delta V$

\[p - p_0 = -\beta \frac{V - V_0}{V_0}\]

substituting $V = m / \rho$

\[p - p_0 = -\beta (1 - \rho/\rho_0) \]

Solving for $\rho$

\[\rho = \rho_0 (1 + (p - p_0)/\beta)\]

Taking a known $\rho_0$ when $p_0$ is 0 (at gage pressure), simplifies to

\[\rho = \rho_0 (1 + p/\beta) \]

Change in Mass

Conservation of mass gives us

\[m_{in} - m_{out} = m_s \]

The stored mass of oil is simply

\[m_s = \rho V \]

Taking the derivative gives us the rate of mass change

\[\dot{m}_{in} - \dot{m}_{out} = \frac{\delta (\rho V)}{\delta t} \]

Here is where the standard hydraulic modeling often makes a simplification.

Correct Derivation (1):

\[\frac{\delta (\rho V)}{\delta t} = \dot{\rho} V + \rho \dot{V} \]

Standard Practice^[1] (2):

\[\color{red} \frac{\delta (\rho V)}{\delta t} = \dot{\rho} V + \rho_0 \dot{V} \]

Given $\dot{\rho} = \rho_0 (\dot{p} / \beta)$, and $q = \dot{m}/\rho_0$ the above is often written as

\[\color{red} q_{in} - q_{out} = (\dot{p} / \beta) V + \dot{V} \]

Example

Problem Definition - Given:

• $M = 10,000 kg$
• $A = 0.01 m^2$
• $\rho_0 = 876 kg/m^3$
• $\beta = 1.2e9 Pa$
• $g = 9.807 m/s^2$

Find the mass flow rate ($\dot{m}$) that provides a sinusodial output of $x$:

\[x(t) = amp \cdot sin(2πtf) + x_0\]

There are 3 fundamental equations needed to solve this problem.

(1) Mass balance:

\[\dot{m} = \dot{\rho} \cdot V + \rho \cdot \dot{V}\]

where $V$ is the cylinder volume $=x \cdot A$

(2) Newton's law:

\[M \cdot \ddot{x} = p \cdot A - m \cdot g\]

(3) Density equation:

\[\rho = \rho_0 (1 + p/\beta) \]

The variables of this system are $x$, $p$, $\rho$, and $\dot{m}$. By including 1 input condition that gives 4 equations and 4 variables to be solved. We will solve the problem 3 different ways

• case 1: guess $\dot{m}$, partial mass balance
• case 2: guess $\dot{m}$, complete mass balance
• case 3: solution, solve $\dot{m}$ directly

To solve this in ModelingToolkit.jl, let's start by defining our parameters and x function

using ModelingToolkit
using DifferentialEquations
using Symbolics
using Plots
using ModelingToolkit: t_nounits as t, D_nounits as D

# parameters -------
pars = @parameters begin
    r₀ = 876 #kg/m^3
    β = 1.2e9 #Pa
    A = 0.01 #m²
    x₀ = 1.0 #m
    M = 10_000 #kg
    g = 9.807 #m/s²
    amp = 5e-2 #m
    f = 15 #Hz
end

dt = 1e-4 #s
t_end = 0.2 #s
time = 0:dt:t_end

x_fun(t,amp,f) = amp*sin(2π*t*f) + x₀

Now, to supply $\dot{m}$ we need an $\dot{x}$ function. This can be automatically generated for us with Symbolics.jl

ẋ_fun = build_function(expand_derivatives(D(x_fun(t,amp,f))), t, amp, f; expression=false)

RuntimeGeneratedFunction(#=in Symbolics=#, #=using Symbolics=#, :((t, amp, f)->begin
          #= /home/runner/.julia/packages/SymbolicUtils/qyMYa/src/code.jl:373 =#
          #= /home/runner/.julia/packages/SymbolicUtils/qyMYa/src/code.jl:374 =#
          #= /home/runner/.julia/packages/SymbolicUtils/qyMYa/src/code.jl:375 =#
          (*)((*)((*)(6.283185307179586, amp), f), NaNMath.cos((*)((*)(6.283185307179586, f), t)))
      end))

As can be seen, we get a cos function as expected taking the derivative of sin. Now let's build the variables and equations of our system. The base equations are generated in a function so we can easily compare the correct derivation of mass balance (density_type = r(t)) with the standard practice (density_type = r₀).

vars = @variables begin
    x(t) = x₀
    ẋ(t)
    ẍ(t)
    p(t) = M*g/A #Pa
    ṁ(t)
    r(t)
    ṙ(t)
end

function get_base_equations(density_type)

    eqs = [
        D(x) ~ ẋ
        D(ẋ) ~ ẍ
        D(r) ~ ṙ

        r ~ r₀*(1 + p/β)

        ṁ ~ ṙ*x*A + (density_type)*ẋ*A
        M*ẍ ~ p*A - M*g
    ]

    return eqs
end

Note: we've only specified the initial values for the known states of x and p. We will find the additional unknown initial conditions before solving. Now we have 7 variables defined and only 6 equations, missing the final driving input equation. Let's build 3 different cases:

case 1:

eqs_ṁ1 = [
    get_base_equations(r₀)...
    ṁ ~ ẋ_fun(t,amp,f)*A*r # (4) Input - mass flow guess
]

case 2:

eqs_ṁ2 = [
    get_base_equations(r)...
    ṁ ~ ẋ_fun(t,amp,f)*A*r # (4) Input - mass flow guess
]

case 3:

eqs_x = [
    get_base_equations(r)...
    x ~ x_fun(t,amp,f) # (4) Input - target x
]

Now we have 3 sets of equations, let's construct the systems and solve. If we start with case 3 with the target $x$ input, notice that the structural_simplify step outputs a system with 0 equations!

@mtkbuild odesys_x = ODESystem(eqs_x, t, vars, pars)

julia> odesys_xModel odesys_x with 0 equations
Unknowns (0):
Parameters (8):
  r₀ [defaults to 876]
  β [defaults to 1.2e9]
  A [defaults to 0.01]
  x₀ [defaults to 1.0]
⋮
Incidence matrix:0×0 SparseArrays.SparseMatrixCSC{Symbolics.Num, Int64} with 0 stored entries

What this means is ModelingToolkit.jl has found that this model can be solved entirely analytically. The full system of equations has been moved to what is called "observables", which can be obtained using the observed() function

julia> observed(odesys_x)15-element Vector{Symbolics.Equation}:
 xˍt(t) ~ 6.283185307179586amp*f*cos(6.283185307179586f*t)
 xˍtt(t) ~ -39.47841760435743amp*(f^2)*sin(6.283185307179586f*t)
 xˍttt(t) ~ -248.05021344239853amp*(f^3)*cos(6.283185307179586f*t)
 x(t) ~ x₀ + amp*sin(6.283185307179586f*t)
 ẋ(t) ~ xˍt(t)
 ẋˍt(t) ~ xˍtt(t)
 ẋˍtt(t) ~ xˍttt(t)
 ẍ(t) ~ ẋˍt(t)
 ẍˍt(t) ~ ẋˍtt(t)
 p(t) ~ (-M*g - M*ẍ(t)) / (-A)
 pˍt(t) ~ (M*ẍˍt(t)) / A
 r(t) ~ r₀*(1 + p(t) / β)
 rˍt(t) ~ (r₀*pˍt(t)) / β
 ṙ(t) ~ rˍt(t)
 ṁ(t) ~ A*r(t)*ẋ(t) + A*x(t)*ṙ(t)

This system can still be "solved" using the same steps to generate an ODESolution which allows us to easily obtain any calculated observed state.

prob_x = ODEProblem(odesys_x, [], (0, t_end))
sol_x = solve(prob_x; saveat=time)
plot(sol_x; idxs=ṁ)

Now let's solve the other system and compare the results.

@mtkbuild odesys_ṁ1 = ODESystem(eqs_ṁ1, t, vars, pars)

julia> odesys_ṁ1Model odesys_ṁ1 with 4 equations
Unknowns (4):
  x(t) [defaults to x₀]
  ẋ(t)
  r(t)
  ṙ(t)
Parameters (8):
  r₀ [defaults to 876]
  β [defaults to 1.2e9]
  A [defaults to 0.01]
  x₀ [defaults to 1.0]
⋮
Incidence matrix:4×7 SparseArrays.SparseMatrixCSC{Symbolics.Num, Int64} with 10 stored entries:
 ⋅  ×  ⋅  ×  ⋅  ⋅  ⋅
 ⋅  ⋅  ×  ⋅  ×  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  ×  ×
 ×  ×  ×  ⋅  ⋅  ⋅  ×

Notice that now, with a simple change of the system input variable, structural_simplify() outputs a system with 4 states to be solved. We can find the initial conditions needed for these states from sol_x and solve.

u0 = [sol_x[s][1] for s in unknowns(odesys_ṁ1)]
prob_ṁ1 = ODEProblem(odesys_ṁ1, u0, (0, t_end))
@time sol_ṁ1 = solve(prob_ṁ1; initializealg=NoInit());

┌ Warning: Initialization system is overdetermined. 3 equations for 0 unknowns. Initialization will default to using least squares. To suppress this warning pass warn_initialize_determined = false.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/cAhZr/src/systems/diffeqs/abstractodesystem.jl:1626
  4.755130 seconds (3.20 M allocations: 233.692 MiB, 1.42% gc time, 99.98% compilation time)

The resulting mass flow rate required to hit the target $x$ position can be seen to be completely wrong. This is the large impact that compressibility can have when high forces are involved.

plot(sol_ṁ1; idxs=ṁ, label="guess", ylabel="ṁ")
plot!(sol_x; idxs=ṁ, label="solution")

If we now solve for case 2, we can study the impact the compressibility derivation

@mtkbuild odesys_ṁ2 = ODESystem(eqs_ṁ2, t, vars, pars)
prob_ṁ2 = ODEProblem(odesys_ṁ2, u0, (0, t_end))
@time sol_ṁ2 = solve(prob_ṁ2; initializealg=NoInit());

┌ Warning: Initialization system is overdetermined. 3 equations for 0 unknowns. Initialization will default to using least squares. To suppress this warning pass warn_initialize_determined = false.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/cAhZr/src/systems/diffeqs/abstractodesystem.jl:1626
  4.599240 seconds (2.94 M allocations: 215.633 MiB, 1.35% gc time, 99.98% compilation time)

As can be seen, a significant error forms between the 2 cases. Plotting first the absolute position.

plot(sol_x; idxs=x, label="solution", ylabel="x")
plot!(sol_ṁ1; idxs=x, label="case 1: r₀")
plot!(sol_ṁ2; idxs=x, label="case 2: r")

And now plotting the difference between case 1 and 2.

plot(time, (sol_ṁ1(time)[x] .- sol_ṁ2(time)[x])/1e-3,
            label="x",
            ylabel="error (case 1 - case 2) [mm]",
            xlabel="t [s]"
        )

Also note the difference in computation.

julia> sol_ṁ1.destatsSciMLBase.DEStats
Number of function 1 evaluations:                  452
Number of function 2 evaluations:                  0
Number of W matrix evaluations:                    30
Number of linear solves:                           240
Number of Jacobians created:                       30
Number of nonlinear solver iterations:             0
Number of nonlinear solver convergence failures:   0
Number of fixed-point solver iterations:                     0
Number of fixed-point solver convergence failures:           0
Number of rootfind condition calls:                0
Number of accepted steps:                          30
Number of rejected steps:                          0

As can be seen, including the detail of full compressibility resulted in more computation: more function evaluations, Jacobians, solves, and steps.

julia> sol_ṁ2.destatsSciMLBase.DEStats
Number of function 1 evaluations:                  528
Number of function 2 evaluations:                  0
Number of W matrix evaluations:                    36
Number of linear solves:                           288
Number of Jacobians created:                       34
Number of nonlinear solver iterations:             0
Number of nonlinear solver convergence failures:   0
Number of fixed-point solver iterations:                     0
Number of fixed-point solver convergence failures:           0
Number of rootfind condition calls:                0
Number of accepted steps:                          34
Number of rejected steps:                          2

ModelingToolkitStandardLibrary.jl

Now let's re-create this example using components from the ModelingToolkitStandardLibrary.jl. It can be shown that by connecting Mass and Volume components that the same exact result is achieved. The important thing is to pay very close attention to the initial conditions.

import ModelingToolkitStandardLibrary.Mechanical.Translational as T
import ModelingToolkitStandardLibrary.Hydraulic.IsothermalCompressible as IC
import ModelingToolkitStandardLibrary.Blocks as B

using DataInterpolations
mass_flow_fun = LinearInterpolation(sol_x[ṁ], sol_x.t)

function MassVolume(; name, dx, drho, dm)

    pars = @parameters begin
        A = 0.01 #m²
        x₀ = 1.0 #m
        M = 10_000 #kg
        g = 9.807 #m/s²
        amp = 5e-2 #m
        f = 15 #Hz
        p_int=M*g/A
        dx=dx
        drho=drho
        dm=dm
    end
    vars = []
    systems = @named begin
        fluid = IC.HydraulicFluid(; density = 876, bulk_modulus = 1.2e9)
        mass = T.Mass(;v=dx,m=M,g=-g)
        vol = IC.Volume(;area=A, x=x₀, p=p_int, dx, drho, dm)
        mass_flow = IC.MassFlow(;p_int)
        mass_flow_input = B.TimeVaryingFunction(;f = mass_flow_fun)
    end

    eqs = [
        connect(mass.flange, vol.flange)
        connect(vol.port, mass_flow.port)
        connect(mass_flow.dm, mass_flow_input.output)
        connect(mass_flow.port, fluid)
    ]

    return ODESystem(eqs, t, vars, pars; systems, name)
end

dx = sol_x[ẋ][1]
drho = sol_x[ṙ][1]
dm = sol_x[ṁ][1]

@mtkbuild odesys = MassVolume(; dx, drho, dm)

prob = ODEProblem(odesys, [], (0, t_end))
sol=solve(prob)

plot(sol; idxs=odesys.vol.x, linewidth=2)
plot!(sol_x; idxs=x)

Momentum Balance

The next challenging aspect of hydraulic modeling is modeling flow through a pipe, which for compressible flow requires resolving the momentum balance equation. To derive the momentum balance we can draw a control volume (cv) in a pipe with area $A$, as shown in the figure below, and apply Newton's second law. Across this control volume from $x_1$ to $x_2$ the pressure will change from $p_1$ to $p_2$. Assuming this is written for an acausal component we put nodes at $p_1$ to $p_2$ which will have equal mass flow $\dot{m}$ entering and exiting the cv^[2].

Now taking the sum of forces acting on the cv we have the pressure forces at each end as well as the viscous drag force from the pipe wall and the body force from gravity. The sum of forces is equal to the product of mass ($\rho V$) and flow acceleration ($\dot{u}$).

\[ \rho V \dot{u} = (p_1 - p_2) A - F_{viscous} + \rho V g\]

where

\[\begin{align} F_{viscous} = A \frac{1}{2} \rho u^2 f \frac{L}{d_h} \end{align}\]

given $f$ is the fully developed flow pipe friction factor for a given shape, $L$ is the pipe length, and $d_h$ is the pipe hydraulic diameter.

The density $\rho$ is an average of $\rho_1$ and $\rho_2$. The velocity is also taken as an average of $u_1$ and $u_2$

\[u_1 = \frac{\dot{m}}{\rho_1 A}\]

\[u_2 = \frac{\dot{m}}{\rho_2 A}\]

Pipe Component

To model a pipe for compressible flow, we can combine the mass balance and momentum balance components to give both mass storage and flow resistance. Furthermore, to provide a more accurate model that allows for wave propagation we can discretize the volume connected by node of equal pressure and mass flow. The diagram below shows an example of discretizing with 3 mass balance volumes and 2 momentum balance resistive elements. Note: the Modelica Standard Library does this in a different way, by combining the mass and momentum balance in a single base class.

Dynamic Volume Component

Both Modelica and SimScape model the actuator component with simply a uniform pressure volume component. The Modelica library defines the base fluids class around the assumption of constant length (see: Object-Oriented Modeling of Thermo-Fluid Systems) and therefore adapting to a component that changes length is not possible. But in cases with long actuators with high dynamics the pressure is not at all uniform, therefore this detail cannot be ignored. Therefore, adding in the momentum balance to provide flow resistance and fluid inertia are necessary. The diagram below shows the design of a DynamicVolume component which includes both mass and momentum balance in addition to discretization by volume. The discretization is similar to the pipe, except the scheme becomes a bit more complicated with the moving wall ($x$). As the volume shrinks, the control volumes will also shrink, however not in unison, but one at a time. In this way, as the moving wall closes, the flow will come from the first volume $cv1$ and travel thru the full size remaining elements ($cv2$, $cv3$, etc.). After the first component length drops to zero, the next element will then start to shrink.

This design has a flaw unfortunately, expanding the system for N=3 gives

What happens when transitioning from one cv to the next, if the moving wall velocity is significant, then an abrupt change occurs due to the $\rho_i \dot{x}$ term. This creates an unstable condition for the solver and results in poor quality/accuracy. To resolve this problem, the mass balance equation is split into 2 parts: mass balance 1 \& 2

\[\text{mass balance 1: } \dot{m}/A = \dot{\rho} x\]

\[\text{mass balance 2: } \dot{m}/A = \rho \dot{x}\]

The below diagram explains how this component is constructed

Now the flows are simplified and are more numerically stable. The acausal connections then handle the proper summing of flows.