working Pipe model

Author: Brad Carman

src/Hydraulic/IsothermalCompressible/Fluids/sae_30_oil_20C.jl

@@ -0,0 +1,7 @@
+@named sae_30_oil_20C = FluidType()
+const sae_30_oil_20C_type = typeof(get_fluid(sae_30_oil_20C))
+
+# https://wiki.anton-paar.com/us-en/engine-oil/ at 20C
+density(::sae_30_oil_20C_type) = 881.5
+bulk_modulus(::sae_30_oil_20C_type) = 1.5e9
+viscosity(::sae_30_oil_20C_type) = 0.23939

src/Hydraulic/IsothermalCompressible/Fluids/water_10C.jl

@@ -0,0 +1,6 @@
+@named water_20C = FluidType()
+const water_20C_type = typeof(get_fluid(water_20C)) 
+
+density(::water_20C_type) = 997             # [kg/m^3]
+bulk_modulus(::water_20C_type) = 2.09e9     # [Pa]
+viscosity(::water_20C_type) = 0.0010016     # [Pa*s] or [kg/m-s]

src/Hydraulic/IsothermalCompressible/Fluids/water_20C.jl

@@ -1,5 +1,5 @@
 """
-Library to model iso-thermal compressible fluid flow 
+Library to model iso-thermal compressible liquid fluid flow 
 """
 module IsothermalCompressible
 
@@ -12,7 +12,7 @@ D = Differential(t)
 export HydraulicPort
 include("utils.jl")
 
-export Source, FixedVolume, LaminarResistance
+export Source, InputSource, FixedVolume, Pipe
 include("components.jl")
 
 

src/Hydraulic/IsothermalCompressible/Fluids/water_30C.jl

@@ -31,15 +31,15 @@ end
 
 
 
-function InputSource(;name, p)
+function InputSource(;name, p_int)
     pars = @parameters begin
-        p = p
+        p_int = p_int
     end
 
     vars = []
 
     systems = @named begin
-        port = HydraulicPort(; p_int = p)
+        port = HydraulicPort(; p_int)
         input = RealInput()
     end
 
@@ -98,7 +98,7 @@ end
 """
     LaminarResistance(fluid, shape=:circle; name, p_int, area, length, perimeter=2*sqrt(area*pi))
 
-fixed fluid volume
+laminar pipe resistance
 
 # Parameters:
 
@@ -110,7 +110,7 @@ fixed fluid volume
 
   - `port`: hydraulic port
 """
-function LaminarResistance(;fluid, shape=Shapes.circle, name, p_int, area, length, perimeter=2*sqrt(area*pi))
+function PipeBase(;name, fluid, shape=Shapes.circle, p_int, area, length, perimeter=2*sqrt(area*pi))
     pars = @parameters begin
         p_int = p_int
         area = area
@@ -120,9 +120,7 @@ function LaminarResistance(;fluid, shape=Shapes.circle, name, p_int, area, lengt
         shape=shape
     end
 
-    vars = @variables begin
-        v(t) = 0
-    end
+    vars = []
 
     systems = @named begin
         port_a = HydraulicPort(; p_int)
@@ -133,19 +131,69 @@ function LaminarResistance(;fluid, shape=Shapes.circle, name, p_int, area, lengt
     # let ----------------------
     Δp = port_a.p - port_b.p
     dm = port_a.dm
+    p = (port_a.p + port_b.p)/2
 
     d_h = 4*area/perimeter
-    rho = density(fluid, port_a.p)
-    Re = rho*v*d_h/viscosity(fluid)
-    f = friction_factor(shape)/Re
-    
+
+    ρ = density(fluid, p)
+    μ = viscosity(fluid)
+    Φ = shape_factor(shape)
+    f = friction_factor(dm, area, d_h, ρ, μ, Φ)
+    u = dm/(ρ*area)
 
     eqs = [        
-        Δp ~ 1/2 * rho * v^2 * f * length
-        dm ~ rho*v*area
-        0 ~ port_a.dm - port_b.dm
+        Δp ~ 1/2 * ρ * u^2 * f * (length/d_h)
+        0 ~ port_a.dm + port_b.dm
     ]
 
     ODESystem(eqs, t, vars, pars; name, systems)
 end
 
+function Pipe(N; name, fluid, shape=Shapes.circle, p_int, area, length, perimeter=2*sqrt(area*pi))
+
+    @assert(N>1, "the pipe component must be defined with more than 1 segment (i.e. N>1), found N=$N")
+
+    pars = @parameters begin
+        p_int = p_int
+        area = area
+        length = length
+        perimeter = perimeter
+        fluid=fluid
+        shape=shape
+    end
+    
+    vars = []
+    
+    ports = @named begin
+        port_a = HydraulicPort(; p_int)
+        port_b = HydraulicPort(; p_int)
+    end
+
+
+    pipe_bases = []
+    for i=1:N-1
+        x = PipeBase(;name=Symbol("p$i"), fluid=ParentScope(fluid), shape=ParentScope(shape), p_int=ParentScope(p_int), area=ParentScope(area), length=ParentScope(length)/(N-1), perimeter=ParentScope(perimeter))
+        push!(pipe_bases, x)
+    end
+
+    volumes = []
+    for i=1:N
+        x = FixedVolume(; name=Symbol("v$i"), fluid=ParentScope(fluid), vol=ParentScope(area)*ParentScope(length)/N, p_int=ParentScope(p_int))
+        push!(volumes, x)
+    end
+    
+
+    eqs = [
+        connect(volumes[1].port, pipe_bases[1].port_a, port_a)
+        connect(volumes[end].port, pipe_bases[end].port_b, port_b)
+    ]
+
+    for i=2:N-1
+        eq = connect(volumes[i].port, pipe_bases[i-1].port_b, pipe_bases[i].port_a)
+        push!(eqs, eq)
+    end
+
+
+    ODESystem(eqs, t, vars, pars; name, systems=[ports; pipe_bases; volumes])
+end
+

src/Hydraulic/IsothermalCompressible/IsothermalCompressible.jl

@@ -12,21 +12,21 @@ Base.@doc """
 port for hydraulic components
 
 # States:
-- `p`: [Pa] pressure
+- `p`: [Pa] gage pressure
 - `dm`: [kg/s] mass flow
 """ HydraulicPort
 
 
 struct FluidType
     name::Symbol
-    value::Int
+    value::Float64
 end
-FluidType(value; name) = FluidType(name, value)
+FluidType(; name) = FluidType(name, float(hash(name)))
 
 
 Base.float(x::FluidType) = Base.float(x.value)
-Base.eltype(x::FluidType) = Int
-Base.promote_type(x::Type{FluidType}, y::Type{Int64}) = Int64
+Base.eltype(x::FluidType) = Float64
+# Base.promote_type(x::Type{FluidType}, y::Type{Int64}) = Int64
 Base.promote_type(x::Type{FluidType}, y::Type{Float64}) = Float64
 Base.convert(T::Type{Float64}, x::FluidType) = float(x)
 
@@ -35,20 +35,10 @@ function Base.show(io::IO, ::MIME"text/plain", x::FluidType)
     print(io, "$(x.name)::FluidType")
 end
 
-get_fluid(x::Float64) = Val(round(Int, x))
-get_fluid(x::Int) = Val(x)
+get_fluid(x::Float64) = Val(x)
 get_fluid(x::FluidType) = Val(x.value)
 
 
-module Fluids
-import ..IsothermalCompressible: FluidType
-using ModelingToolkit
-@named water = FluidType(0)
-@named sae_30_oil = FluidType(1)
-end
-
-
-
 density(fluid) = density(get_fluid(fluid))
 @register_symbolic density(fluid)
 
@@ -59,18 +49,12 @@ viscosity(fluid) = viscosity(get_fluid(fluid))
 @register_symbolic viscosity(fluid)
 
 
-density(::Val{Fluids.water.value}) = 997 # [kg/m^3]
-bulk_modulus(::Val{Fluids.water.value}) = 2.09e9 # [Pa]
-viscosity(::Val{Fluids.water.value}) = 0.0010016 # [Pa*s] or [kg/m-s]
-
-# https://wiki.anton-paar.com/us-en/engine-oil/ at 20C
-density(::Val{Fluids.sae_30_oil.value}) = 881.5
-bulk_modulus(::Val{Fluids.sae_30_oil.value}) = 1.5e9
-viscosity(::Val{Fluids.sae_30_oil.value}) = 0.23939
-
-
 density(fluid, p) = density(fluid)*(1 + p/bulk_modulus(fluid))
 
+# Fluids
+include("Fluids/water_20C.jl")
+include("Fluids/sae_30_oil_20C.jl")
+
 
 
 
@@ -101,11 +85,57 @@ module Shapes
 import ..IsothermalCompressible: ShapeType
 using ModelingToolkit
 @named circle = ShapeType(0)
-@named triangle = ShapeType(1)
+@named square = ShapeType(1)
+@named triangle = ShapeType(2)
+@named rectangle_r2 = ShapeType(3)
+@named rectangle_r3 = ShapeType(4)
+@named rectangle_r4 = ShapeType(5)
+@named rectangle_r8 = ShapeType(6)
+@named rectangle_rinf = ShapeType(7)
 end
 
+shape_factor(shape) = shape_factor(get_shape(shape))
+@register_symbolic shape_factor(shape)
+
+shape_factor(::Val{Shapes.circle.value}) = 64
+shape_factor(::Val{Shapes.square.value}) = 57
+shape_factor(::Val{Shapes.triangle.value}) = 53
+shape_factor(::Val{Shapes.rectangle_r2.value}) = 62
+shape_factor(::Val{Shapes.rectangle_r3.value}) = 69
+shape_factor(::Val{Shapes.rectangle_r4.value}) = 73
+shape_factor(::Val{Shapes.rectangle_r8.value}) = 82
+shape_factor(::Val{Shapes.rectangle_rinf.value}) = 96
+
+function friction_factor(dm, area, d_h, rho, mu, shape_factor)
+
+    u = abs(dm)/(rho*area)
+    Re = maximum([rho*u*d_h/mu, 1]) 
+
+    f_laminar = shape_factor/Re
+
+
+    f_turbulent = (0.79*log(Re)-1.64)^(-2)
+
+    f = transition(2000, 3000, f_laminar, f_turbulent, Re)
+
+    return f
+end
+@register_symbolic friction_factor(dm, area, d_h, rho, mu, shape_factor)
+
+
+
+function transition(x1,x2,y1,y2,x)
+
+    if x < x1
+        return y1
+    elseif x > x2
+        return y2
+    else
+        u = (x-x1)/(x2-x1)
+        blend = 3*u^2 - 2*u^3;
+        return (1-blend)*y1 + blend*y2;
+    end
+    
+end
 
-friction_factor(shape) = friction_factor(get_shape(shape))
-@register_symbolic friction_factor(shape)
 
-friction_factor(::Val{Shapes.circle.value}) = 64