added Actuator

Author: Brad Carman

src/Hydraulic/IsothermalCompressible/IsothermalCompressible.jl

@@ -5,6 +5,7 @@ module IsothermalCompressible
 
 using ModelingToolkit, Symbolics
 using ...Blocks: RealInput, RealOutput
+using ...Mechanical.Translational: MechanicalPort
 
 @parameters t
 D = Differential(t)

src/Hydraulic/IsothermalCompressible/components.jl

@@ -197,3 +197,43 @@ function Pipe(N; name, fluid, shape=Shapes.circle, p_int, area, length, perimete
     ODESystem(eqs, t, vars, pars; name, systems=[ports; pipe_bases; volumes])
 end
 
+
+function Actuator(direction; name, fluid, p_int, x_int, area, dead_volume)
+
+    pars = @parameters begin
+        fluid = fluid
+        p_int = p_int
+        x_int = x_int
+        area = area
+        dead_volume = dead_volume
+    end
+
+    vars = @variables begin
+        x(t) = x_int
+        dx(t) = 0
+        rho(t) = density(fluid, p_int)
+        drho(t) = 0
+    end
+
+    systems = @named begin
+        port = HydraulicPort(; p_int)
+        flange = MechanicalPort()
+    end
+
+    # let -------------
+    vol = dead_volume + area*x
+
+    eqs = [
+        D(x) ~ dx
+        D(rho) ~ drho
+
+        dx ~ flange.v*direction
+        rho ~ density(fluid, port.p)
+        
+        port.dm ~ drho*vol + rho*area*dx
+        flange.f ~ -port.p*area*direction
+    ]
+
+    ODESystem(eqs, t, vars, pars; name, systems)
+end
+

src/Hydraulic/IsothermalCompressible/utils.jl

@@ -106,6 +106,8 @@ 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)

test/Hydraulic/isothermal_compressible.jl

@@ -1,6 +1,7 @@
 using ModelingToolkit, OrdinaryDiffEq, Test
 import ModelingToolkitStandardLibrary.Hydraulic.IsothermalCompressible as IC
 import ModelingToolkitStandardLibrary.Blocks as B
+import ModelingToolkitStandardLibrary.Mechanical.Translational as T
 
 
 
@@ -67,23 +68,89 @@ prob25 = ODEProblem(cc25, [], (0,t_end))
 
 dt = 1e-4
 NEWTON = NLNewton(check_div = false, always_new = true, max_iter=100, relax=4//10)
-sol0 = solve(prob0, ImplicitEuler(nlsolve=NEWTON); adaptive=false, dt, initializealg=NoInit())  
-sol1 = solve(prob1, ImplicitEuler(nlsolve=NEWTON); adaptive=false, dt, initializealg=NoInit())  
-sol2 = solve(prob2, ImplicitEuler(nlsolve=NEWTON); adaptive=false, dt, initializealg=NoInit())  
-sol5 = solve(prob5, ImplicitEuler(nlsolve=NEWTON); adaptive=false, dt, initializealg=NoInit())  
-sol25 = solve(prob25, ImplicitEuler(nlsolve=NEWTON); adaptive=false, dt, initializealg=NoInit())  
+# sol0 = solve(prob0, ImplicitEuler(nlsolve=NEWTON); adaptive=false, dt, initializealg=NoInit())  
+# sol1 = solve(prob1, ImplicitEuler(nlsolve=NEWTON); adaptive=false, dt, initializealg=NoInit())  
+# sol2 = solve(prob2, ImplicitEuler(nlsolve=NEWTON); adaptive=false, dt, initializealg=NoInit())  
+# sol5 = solve(prob5, ImplicitEuler(nlsolve=NEWTON); adaptive=false, dt, initializealg=NoInit())  
+# sol25 = solve(prob25, ImplicitEuler(nlsolve=NEWTON); adaptive=false, dt, initializealg=NoInit())  
+
+
+
 
+# using GLMakie
+# fig = Figure()
+# ax = Axis(fig[1,1])
+# lines!(ax, sol0.t, sol0[sys1.res.port_a.p] .- sol0[sys1.res.port_b.p]; label="pipe base, oil", linestyle=:dash)
+# lines!(ax, sol1.t, sol1[sys1.res.port_a.p] .- sol1[sys1.res.port_b.p]; label="pipe base")
+# lines!(ax, sol2.t, sol2[sys2.res.port_a.p] .- sol2[sys2.res.port_b.p]; label="pipe N=2")
+# lines!(ax, sol5.t, sol5[sys5.res.port_a.p] .- sol5[sys5.res.port_b.p]; label="pipe N=5")
+# lines!(ax, sol25.t, sol25[sys25.res.port_a.p] .- sol25[sys25.res.port_b.p]; label="pipe N=25")
+# axislegend(ax)
+# fig
 
 
 
+
+function actuator_system(; name, fluid_a, fluid_b = fluid_a)
+
+    pars = @parameters begin
+        fluid_a = fluid_a
+        fluid_b = fluid_b
+    end
+
+    systems = @named begin
+        stp = B.Step(;height = 10e5, offset = 0, start_time = 0.005, duration = Inf, smooth = 0)
+
+        src_a = IC.InputSource(;p_int=0)
+        res_a = IC.Pipe(5; p_int=0, area=0.01, length=500.0, fluid=fluid_a)
+        act_a = IC.Actuator(+1; fluid=fluid_a, p_int=0, x_int=0, area=0.1, dead_volume=100)
+
+        src_b = IC.InputSource(;p_int=0)
+        res_b = IC.Pipe(5; p_int=0, area=0.01, length=500.0, fluid=fluid_b)
+        act_b = IC.Actuator(-1; fluid=fluid_b, p_int=0, x_int=0, area=0.1, dead_volume=100)
+
+        m = T.Mass(; m=1000)
+        
+    end
+
+
+
+
+    eqs = [
+        connect(stp.output, src_a.input, src_b.input)
+
+        connect(src_a.port, res_a.port_a)
+        connect(res_a.port_b, act_a.port)
+
+        connect(src_b.port, res_b.port_a)
+        connect(res_b.port_b, act_b.port)
+
+        connect(act_a.flange, act_b.flange, m.flange)
+    ]
+
+    ODESystem(eqs, t, [], pars; name, systems)
+end
+
+@named act_sys = actuator_system(; fluid_a = IC.water_20C)
+cc = structural_simplify(act_sys; allow_parameter=false)
+
+t_end = 0.5
+prob1 = ODEProblem(cc, [], (0,t_end), [complete(act_sys).fluid_b => user_oil])
+prob2 = ODEProblem(cc, [], (0,t_end))
+
+sol1 = solve(prob1, ImplicitEuler(nlsolve=NEWTON); adaptive=false, dt, initializealg=NoInit())  
+sol2 = solve(prob2, ImplicitEuler(nlsolve=NEWTON); adaptive=false, dt, initializealg=NoInit())  
+
+
 using GLMakie
 fig = Figure()
 ax = Axis(fig[1,1])
-lines!(ax, sol0.t, sol0[sys1.res.port_a.p] .- sol0[sys1.res.port_b.p]; label="pipe base, oil", linestyle=:dash)
-lines!(ax, sol1.t, sol1[sys1.res.port_a.p] .- sol1[sys1.res.port_b.p]; label="pipe base")
-lines!(ax, sol2.t, sol2[sys2.res.port_a.p] .- sol2[sys2.res.port_b.p]; label="pipe N=2")
-lines!(ax, sol5.t, sol5[sys5.res.port_a.p] .- sol5[sys5.res.port_b.p]; label="pipe N=5")
-lines!(ax, sol25.t, sol25[sys25.res.port_a.p] .- sol25[sys25.res.port_b.p]; label="pipe N=25")
-axislegend(ax)
-fig
-
+lines!(ax, sol1.t, sol1[act_sys.act_a.x]; color=:red)
+lines!(ax, sol1.t, sol1[act_sys.act_b.x]; color=:red, linestyle=:dash)
+lines!(ax, sol2.t, sol2[act_sys.act_a.x]; color=:blue)
+lines!(ax, sol2.t, sol2[act_sys.act_b.x]; color=:blue, linestyle=:dash)
+
+ax = Axis(fig[2,1])
+lines!(ax, sol1.t, (sol1[act_sys.act_a.port.p] .- sol1[act_sys.act_b.port.p])/1e5; color=:red)
+lines!(ax, sol2.t, (sol2[act_sys.act_a.port.p] .- sol2[act_sys.act_b.port.p])/1e5; color=:blue)
+fig