Adds hydraulic library

src/Hydraulic/Hydraulic.jl

@@ -0,0 +1,16 @@
+module Hydraulic
+using ModelingToolkit, OrdinaryDiffEq
+
+@parameters t
+D = Differential(t)
+
+export HydraulicPort, FluidProperties
+include("utils.jl")
+
+export LocalRestriction, ConstantVolume, Reservoir
+include("components.jl")
+
+export PressureSource, MassFlowRateSource
+include("sources.jl")
+
+end

src/Hydraulic/components.jl

@@ -0,0 +1,109 @@
+"""
+    Reservoir(; name, p)
+
+Infinite reservoir; set BC for hydraulic system
+
+# Connectors:
+- `a` port [`HydraulicPort`](@ref)
+
+# Parameters:
+- `p`: [`Pa`] Pressure of the reservoir
+"""
+function Reservoir(; name, p)
+    @parameters p = p
+    @named a = HydraulicPort()
+    eqs = [a.p ~ p]
+    return ODESystem(eqs, t, [], [p]; systems = [a], name = name)
+end
+
+"""
+    LocalRestriction(; name, A, Cd, Re_crit)
+
+# States:
+
+# Connectors:
+- `a` left port [`HydraulicPort`](@ref)
+- `b` right port [`HydraulicPort`](@ref)
+
+# Parameters:
+- `A`: [`m^2`] Effective area of the restriction
+- `Cd`: [] Discharge coefficient
+- `Re_crit=150`: [] Critical Reynolds number
+"""
+function LocalRestriction(; name, A, Cd, Re_crit = 150)
+    @named a = HydraulicPort()
+    @named b = HydraulicPort()
+
+    pars = @parameters begin
+        A = A
+        Cd = Cd
+        p_atm# , [scope=:parent]
+        nu_atm# , [scope=:parent]
+        beta_atm# , [scope=:parent]
+        rho_atm# , [scope=:parent]
+    end
+    # lift the fluid parameters into parent scope
+    p_atm = ParentScope(p_atm)
+    nu_atm = ParentScope(nu_atm)
+    beta_atm = ParentScope(beta_atm)
+    rho_atm = ParentScope(rho_atm)
+
+    pars = [A, Cd, p_atm, nu_atm, beta_atm, rho_atm]
+
+    rho = calc_density((b.p + a.p) / 2, rho_atm, p_atm, beta_atm)
+
+    Δp = b.p - a.p
+    p_cr = pi / 4 * rho / (2 * A) * (Re_crit * nu_atm / Cd)^2
+
+    eqs = [
+        # Momentum balance
+        b.m_flow ~ Cd * A * sqrt(2 * rho) * regRoot(Δp, p_cr)
+        # Mass balance
+        a.m_flow + b.m_flow ~ 0
+    ]
+    return ODESystem(eqs, t, [], pars; name = name, systems = [a, b])
+end
+
+"""
+    ConstantVolume(; name, V, p_start)
+
+Constant volume chamber.
+
+# States:
+- `p`: [`Pa`] Pressure inside the volume
+
+# Connectors:
+- `a` port [`HydraulicPort`](@ref)
+
+# Parameters:
+- `V`: [`m^3`] Constant volume
+- `p_start=0.0`: [`Pa`] Initial pressure inside the volume
+"""
+function ConstantVolume(; name, V, p_start = 0.0)
+    @named a = HydraulicPort()
+    @variables p(t) = p_start # pressure inside the volume
+    pars = @parameters begin
+        V = V
+        p_atm# , [scope=:parent]
+        nu_atm# , [scope=:parent]
+        beta_atm# , [scope=:parent]
+        rho_atm# , [scope=:parent]
+    end
+    # lift the fluid parameters into parent scope
+    p_atm = ParentScope(p_atm)
+    nu_atm = ParentScope(nu_atm)
+    beta_atm = ParentScope(beta_atm)
+    rho_atm = ParentScope(rho_atm)
+
+    pars = [V, p_atm, nu_atm, beta_atm, rho_atm]
+
+    rho = calc_density(p, rho_atm, p_atm, beta_atm)
+
+    eqs = [
+        # Mass conservation
+        D(p) ~ 1 / (Symbolics.derivative(rho, p) * V) * a.m_flow
+        a.p ~ p # no pressure loss from inlet to volume
+    ]
+
+    return ODESystem(eqs, t, [p], pars; name = name, systems = [a])
+end

src/Hydraulic/sources.jl

@@ -0,0 +1,45 @@
+"""
+    PressureSource(; name, dp)
+
+Provides a pressure difference between port `a` and port `b`.
+
+# Connectors:
+- `a` left port [`HydraulicPort`](@ref)
+- `b` right port [`HydraulicPort`](@ref)
+
+# Parameters:
+- `dp`: [`Pa`] Pressure difference `a.p ~ b.p + dp`
+"""
+function PressureSource(; name, dp)
+    @named a = HydraulicPort()
+    @named b = HydraulicPort()
+    @parameters dp = dp
+    eqs = [
+        a.p ~ b.p + dp
+        a.m_flow + b.m_flow ~ 0
+    ]
+    return ODESystem(eqs, t, [], [dp]; systems = [a, b], name = name)
+end
+
+"""
+    MassFlowRateSource(; name, m_flow)
+
+Provides a mass flow. A positive value causes liquid to flow from `a` to `b`.
+
+# Connectors:
+- `a` left port [`HydraulicPort`](@ref)
+- `b` right port [`HydraulicPort`](@ref)
+
+# Parameters:
+- `m_flow`: [`kg/s`] Prescribed mass flow rate
+"""
+function MassFlowRateSource(; name, m_flow)
+    @named a = HydraulicPort()
+    @named b = HydraulicPort()
+    @parameters m_flow = m_flow
+    eqs = [
+        a.m_flow + b.m_flow ~ 0
+        a.m_flow ~ m_flow
+    ]
+    return ODESystem(eqs, t, [], [m_flow]; systems = [a, b], name = name)
+end

src/Hydraulic/utils.jl

@@ -0,0 +1,74 @@
+@connector function HydraulicPort(; name, p_start = 1e5, m_flow_start = 0.0)
+    @variables p(t) = p_start
+    @variables m_flow(t) = m_flow_start [connect = Flow]
+    return ODESystem(Equation[], t, [p, m_flow], []; name = name)
+end
+@doc """
+    HydraulicPort(; name, p_start=1e5, m_flow_start=0.0)
+
+Port for hydraulic systems.
+
+# States:
+- `p(t)`: [`Pa`] Pressure
+- `m_flow(t)`: [`kg/s`] Mass flow rate
+
+# Parameters:
+- `p_start`: [`Pa`] Initial pressure
+- `m_flow_start`: [`kg/s`] Initial mass flow rate
+""" HydraulicPort
+
+# smooth approximation of x * abs(x)
+regAbs(x, delta = 0.01) = x * sqrt(x^2 + delta^2)
+
+# Smoothed version of sign(x)*sqrt(abs(x)), which has a finite derivative at x=0
+regRoot(x, delta = 0.01) = x / (x^2 + delta^2)^0.25
+
+"""
+    FluidProperties(;name, p_atm=101325.0, nu_atm=1.0034e-6, beta_atm=2.1791e9, rho_atm=998.21)
+
+Defines the properties of the Liquid.
+
+# Example:
+```
+@named fluid_props = FluidProperties(;rho_atm=998.21)
+model = extend(model, fluid_props)
+```
+
+# Parameters:
+- `p_atm = 101325.0`: [Pa] Atmospheric pressure
+- `nu_atm = 1.0034e-6`: [m^2/s] Kinematic viscosity at atmospheric pressure
+- `beta_atm = 2.1791e9`: [Pa] Isothermal bulk modulus at atmospheric pressure
+- `rho_atm = 998.21`: [kg/m^3] Liquid density at atmospheric pressure
+"""
+function FluidProperties(;
+    name,
+    p_atm = 101325.0,
+    nu_atm = 1.0034e-6,
+    beta_atm = 2.1791e9,
+    rho_atm = 998.21,
+)
+    pars = @parameters p_atm = p_atm nu_atm = nu_atm beta_atm = beta_atm rho_atm = rho_atm
+    ODESystem(Equation[], t, [], pars; name = name)
+end
+
+"""
+    calc_density(p, rho_atm, p_atm, beta_atm)
+
+Computes the density of the liquid for the given conditions.
+
+Gholizadeh, Hossein, Richard Burton, and Greg Schoenau. “Fluid Bulk Modulus: Comparison of Low 
+Pressure Models.” International Journal of Fluid Power 13, no. 1 (January 2012): 7–16. 
+https://doi.org/10.1080/14399776.2012.10781042.
+
+Assumptions:
+- Zero air in the liquid
+- Constant bulk modulus
+
+# Parameters:
+- `rho_atm`: [kg/m^3] Liquid density at atmospheric pressure
+- `p_atm`: [Pa] Atmospheric pressure
+- `beta_atm`: [Pa] Isothermal bulk modulus at atmospheric pressure
+"""
+function calc_density(p, rho_atm, p_atm, beta_atm)
+    return rho_atm * exp((p - p_atm) / beta_atm)
+end

src/ModelingToolkitStandardLibrary.jl

@@ -5,5 +5,6 @@ include("Mechanical/Mechanical.jl")
 include("Electrical/Electrical.jl")
 include("Magnetic/Magnetic.jl")
 include("Thermal/Thermal.jl")
+include("Hydraulic/Hydraulic.jl")
 
 end

test/Hydraulic/hydraulic.jl

@@ -0,0 +1,38 @@
+using ModelingToolkitStandardLibrary.Hydraulic
+using ModelingToolkit, OrdinaryDiffEq, Test
+# using Plots
+
+@parameters t
+D = Differential(t)
+
+@testset "restriction and volume (shows dynamic compressibility)" begin
+    @named reservoir_left = Reservoir(; p = 2e5)
+    @named valve = LocalRestriction(; A = 1e-4, Cd = 0.64)
+    @named volume = ConstantVolume(; V = 0.001)
+    connections = [
+        connect(reservoir_left.a, valve.a)
+        connect(valve.b, volume.a)
+    ]
+    @named model = ODESystem(connections, t; systems = [reservoir_left, valve, volume])
+
+    # Water
+    @named fluid_props = FluidProperties(;
+        p_atm = 101325.0,
+        nu_atm = 1.0034e-6,
+        beta_atm = 2.1791e9,
+        rho_atm = 998.21,
+    )
+    model = extend(model, fluid_props)
+
+    sys = structural_simplify(model)
+    prob = ODEProblem(sys, [volume.p => 1e5], (0.0, 0.5e-3))
+    sol = solve(prob, Tsit5(), reltol = 1e-6)
+
+    @test sol.retcode == :Success
+    @test sol[volume.p][end] ≈ 2e5 atol = 1
+    @test sol[volume.a.m_flow][end] ≈ 0 atol = 1e-2
+
+    # p1=plot(sol; vars=[volume.p / 1e6], ylabel="p in MPa", label="Pressure of Volume")
+    # p2=plot(sol; vars=[volume.a.m_flow], ylabel="m_flow in kg/s", label="Mass Flow Rate into the volume")
+    # plot(p1, p2, layout=(2,1))
+end

test/runtests.jl

@@ -20,3 +20,6 @@ using SafeTestsets
 
 # Mechanical
 @safetestset "Mechanical" begin include("Mechanical/rotational.jl") end
+
+# Hydraulic
+@safetestset "Hydraulic" begin include("Hydraulic/hydraulic.jl") end