MTK9 Initialization Updates

Removed references to p_int and started debugging the first Hydraulics test case (Fluid Domain and Tube)

Author: mgithubk

src/Hydraulic/IsothermalCompressible/components.jl

@@ -54,7 +54,7 @@ Provides an "open" boundary condition for a hydraulic port such that mass flow `
 end
 
 """
-    TubeBase(add_inertia = true; p_int, area, length_int, head_factor = 1, perimeter = 2 * sqrt(area * pi), shape_factor = 64, name)
+    TubeBase(add_inertia = true; area, length_int, head_factor = 1, perimeter = 2 * sqrt(area * pi), shape_factor = 64, name)
 
 Variable length internal flow model of the fully developed incompressible flow friction.  Includes optional inertia term when `add_inertia = true` to model wave propagation.  Hydraulic ports have equal flow but variable pressure.  Density is averaged over the pressures, used to calculated average flow velocity and flow friction.
 
@@ -63,7 +63,6 @@ Variable length internal flow model of the fully developed incompressible flow f
 - `ddm`: [kg/s^2] Rate of change of mass flow rate in control volume.
 
 # Parameters:
-- `p_int`: [Pa] initial pressure
 - `area`: [m^2] tube cross sectional area
 - `length_int`: [m] initial tube length
 - `perimeter`: [m] perimeter of the pipe cross section (needed only for non-circular pipes)
@@ -145,12 +144,11 @@ Variable length internal flow model of the fully developed incompressible flow f
 end
 
 """
-    Tube(N, add_inertia=true; p_int, area, length, head_factor=1, perimeter = 2 * sqrt(area * pi), shape_factor = 64, name)
+    Tube(N, add_inertia=true; area, length, head_factor=1, perimeter = 2 * sqrt(area * pi), shape_factor = 64, name)
 
 Constant length internal flow model discretized by `N` (`FixedVolume`: `N`, `TubeBase`:`N-1`) which models the fully developed flow friction, compressibility (when `N>1`), and inertia effects when `add_inertia = true`.  See `TubeBase` and `FixedVolume` for more information.
 
 # Parameters:
-- `p_int`: [Pa] initial pressure
 - `area`: [m^2] tube cross sectional area
 - `length`: [m] real length of the tube
 - `perimeter`: [m] perimeter of the pipe cross section (needed only for non-circular pipes)
@@ -198,7 +196,7 @@ Constant length internal flow model discretized by `N` (`FixedVolume`: `N`, `Tub
     for i in 1:(N - 1)
         x = TubeBase(add_inertia; name = Symbol("p$i"),
             shape_factor = ParentScope(shape_factor),
-            p_int = ParentScope(p_int), area = ParentScope(area),
+            area = ParentScope(area),
             length_int = ParentScope(length) / (N - 1),
             head_factor = ParentScope(head_factor),
             perimeter = ParentScope(perimeter))
@@ -208,8 +206,7 @@ Constant length internal flow model discretized by `N` (`FixedVolume`: `N`, `Tub
     volumes = []
     for i in 1:N
         x = FixedVolume(; name = Symbol("v$i"),
-            vol = ParentScope(area) * ParentScope(length) / N,
-            p_int = ParentScope(p_int))
+            vol = ParentScope(area) * ParentScope(length) / N)
         push!(volumes, x)
     end
 
@@ -397,29 +394,27 @@ end
 end
 
 """
-    FixedVolume(; vol, p_int, name)
+    FixedVolume(; vol, name)
 
 Fixed fluid volume.
 
 # Parameters:
 - `vol`: [m^3] fixed volume
-- `p_int`: [Pa] initial pressure
 
 # Connectors:
 - `port`: hydraulic port
 """
-@component function FixedVolume(; vol, p_int, name)
+@component function FixedVolume(; vol, name)
     pars = @parameters begin
-        p_int = p_int
         vol = vol
     end
 
     systems = @named begin
-        port = HydraulicPort(; p_int)
+        port = HydraulicPort(;)
     end
 
     vars = @variables begin
-        rho(t), [guess = liquid_density(port)]
+        rho(t)#, [guess = liquid_density(port)]
         drho(t), [guess = 0]
     end
 

src/Hydraulic/IsothermalCompressible/sources.jl

@@ -55,13 +55,10 @@ end
 @deprecate Source FixedPressure
 
 """
-    Pressure(; p_int, name)
+    Pressure(; name)
 
 input pressure source
 
-# Parameters:
-- `p_int`: [Pa] initial pressure (set by `p_int` argument)
-
 # Connectors:
 - `port`: hydraulic port
 - `p`: real input 

test/Hydraulic/isothermal_compressible.jl

@@ -8,19 +8,19 @@ using ModelingToolkitStandardLibrary.Blocks: Parameter
 
 NEWTON = NLNewton(check_div = false, always_new = true, max_iter = 100, relax = 9 // 10)
 
-@testset "Fluid Domain and Tube" begin
+#@testset "Fluid Domain and Tube" begin
     function System(N; bulk_modulus, name)
         pars = @parameters begin
             bulk_modulus = bulk_modulus
         end
 
         systems = @named begin
             fluid = IC.HydraulicFluid(; bulk_modulus)
-            stp = B.Step(; height = 10e5, offset = 0, start_time = 0.005, duration = Inf,
+            stp = B.Step(; height = 2*101325, offset = 101325, start_time = 0.05, duration = Inf,
                 smooth = true)
-            src = IC.Pressure(; p_int = 0)
-            vol = IC.FixedVolume(; p_int = 0, vol = 10.0)
-            res = IC.Tube(N; p_int = 0, area = 0.01, length = 50.0)
+            src = IC.Pressure(;)
+            vol = IC.FixedVolume(; vol = 10.0)
+            res = IC.Tube(N; area = 0.01, length = 50.0)
         end
 
         eqs = [connect(stp.output, src.p)
@@ -32,35 +32,84 @@ NEWTON = NLNewton(check_div = false, always_new = true, max_iter = 100, relax =
     end
 
     @named sys1_2 = System(1; bulk_modulus = 2e9)
-    @named sys1_1 = System(1; bulk_modulus = 1e9)
+    @named sys1_1 = System(1; bulk_modulus = 1e7)
     @named sys5_1 = System(5; bulk_modulus = 1e9)
 
-    syss = structural_simplify.([sys1_2, sys1_1, sys5_1])
-    probs = [ODEProblem(sys, ModelingToolkit.missing_variable_defaults(sys), (0, 0.05))
-             for sys in syss] #
-    sols = [solve(prob, ImplicitEuler(nlsolve = NEWTON); initializealg = NoInit(),
-                dt = 1e-4, adaptive = false)
-            for prob in probs]
+#    syss = structural_simplify.([sys1_2, sys1_1, sys5_1])
 
-    s1_2 = complete(sys1_2)
+# system 1_1
+    sys1_1 = structural_simplify(sys1_1)
+    initialization_eqs = [sys1_1.vol.port.p ~ 101325, sys1_1.res.port_a.dm ~ 0]
+    initsys1_1 = ModelingToolkit.generate_initializesystem(sys1_1;initialization_eqs)
+    # here the structural_simplify of the initialization system fails, this is an indication that something needs to be fixed with the way the defaults/initial equaitons are defined
+
+    initsys1_1 = structural_simplify(initsys1_1)
+    initprob1_1 = NonlinearProblem(initsys1_1, [t=>0])
+    initsol1_1 = solve(initprob1_1)
+    initsol1_1[sys.src.port.p]
+    initsol1_1[sys.src.port.dm]
+    initsol1_1[sys.vol.port.p]
+#    probs = [ODEProblem(sys, ModelingToolkit.missing_variable_defaults(sys), (0, 0.05))
+#             for sys in syss] #
+    prob1_1 = ODEProblem(sys1_1, ModelingToolkit.missing_variable_defaults(sys), (0, 1.05))
+              #
+
+#    sols = [solve(prob, ImplicitEuler(nlsolve = NEWTON); dt = 1e-4, adaptive = false)
+#            for prob in probs]
+
+    sol1_1 = solve(prob1_1, ImplicitEuler(nlsolve = NEWTON); dt = 1e-4, adaptive = false)
+
+#    s1_2 = complete(sys1_2)
     s1_1 = complete(sys1_1)
+#    s5_1 = complete(sys5_1)
+
+# system 1_2
+    sys1_2 = structural_simplify(sys1_2)
+    initialization_eqs = [sys1_2.vol.port.p ~ 101800, sys1_2.res.port_a.dm ~ 0]
+    initsys1_2 = ModelingToolkit.generate_initializesystem(sys1_2;initialization_eqs)
+
+    initsys1_2 = structural_simplify(initsys1_2)
+    initprob1_2 = NonlinearProblem(initsys1_2, [t=>0])
+    initsol1_2 = solve(initprob1_2)
+
+    prob1_2 = ODEProblem(sys1_2, ModelingToolkit.missing_variable_defaults(sys), (0, 1.05))
+    sol1_2 = solve(prob1_2, ImplicitEuler(nlsolve = NEWTON); dt = 1e-4, adaptive = false)
+    s1_2 = complete(sys1_2)
+# system 5-1
+
+    sys5_1 = structural_simplify(sys5_1)
+    initialization_eqs = [sys5_1.vol.port.p ~ 101325,sys5_1.res.p1.port_a.dm ~ 0, 
+    sys5_1.res.p2.port_a.dm ~ 0,sys5_1.res.p3.port_a.dm ~ 0, sys5_1.res.p4.port_a.dm ~ 0,
+    sys5_1.res.p1.port_a.p ~ 101325,
+    sys5_1.res.p2.port_a.p ~ 101325,
+    sys5_1.src.port.dm~0]
+    initsys5_1 = ModelingToolkit.generate_initializesystem(sys5_1;initialization_eqs)
+
+    initsys5_1 = structural_simplify(initsys5_1)
+    initprob5_1 = NonlinearProblem(initsys5_1, [t=>0])
+    initsol5_1 = solve(initprob5_1)
+
+    prob5_1 = ODEProblem(sys5_1, ModelingToolkit.missing_variable_defaults(sys), (0, 0.05))
+    sol5_1 = solve(prob5_1, ImplicitEuler(nlsolve = NEWTON); dt = 1e-4, adaptive = false)
     s5_1 = complete(sys5_1)
 
     # higher stiffness should compress more quickly and give a higher pressure
-    @test sols[1][s1_2.vol.port.p][end] > sols[2][s1_1.vol.port.p][end]
+#    @test sols[1][s1_2.vol.port.p][end] > sols[2][s1_1.vol.port.p][end]
 
     # N=5 pipe is compressible, will pressurize more slowly
     @test sols[2][s1_1.vol.port.p][end] > sols[3][s5_1.vol.port.p][end]
 
-    # fig = Figure()
-    # ax = Axis(fig[1,1])
-    # # hlines!(ax, 10e5)
-    # lines!(ax, sols[1][s1_2.vol.port.p])
+     fig = Figure()
+     ax = Axis(fig[1,1])
+     # hlines!(ax, 10e5)
+    #  lines!(ax, sols[1][s1_2.vol.port.p])
     # lines!(ax, sols[2][s1_1.vol.port.p])
     # lines!(ax, sols[3][s5_1.vol.port.p])
-    # fig
+    lines!(ax, sol1_1[s1_1.vol.rho])
+    lines!(ax, sol1_2[s1_2.vol.rho])
+     fig
 
-end
+#end
 
 @testset "Valve" begin
     function System(; name)
@@ -353,4 +402,53 @@ end
     # fig
 end
 
+@testset "Component Flow Reversals" begin
+# Check Component Flow Reversals
+    function System(; name)
+        pars = []
+
+        systems = @named begin
+            fluid = IC.HydraulicFluid()
+            source = IC.Pressure()
+            sink = IC.FixedPressure(; p = 101325)
+            pipe = IC.Tube(1, false; area = 0.1, length =.1, head_factor = 1)
+            osc = Sine(; frequency = 0.01, amplitude = 100, offset = 101325)
+        end
+
+        eqs = [connect(fluid, pipe.port_a)
+            connect(source.port, pipe.port_a)
+            connect(pipe.port_b, sink.port)
+            connect(osc.output, source.p)]
+
+        ODESystem(eqs, t, [], []; systems)
+    end
+
+    @named sys = System()
+
+    syss = structural_simplify.([sys])
+    tspan = (0.0, 1000.0)
+    prob = ODEProblem(sys, tspan)  # u0 guess can be supplied or not
+    @time sol = solve(prob)
+    
+end
+
+@testset "Tube Discretization" begin
+    # Check Tube Discretization
+end
+
+@testset "Pressure BC" begin
+    # Ensure Pressure Boundary Condition Works
+end
+
+@testset "Massflow BC" begin
+    # Ensure Massflow Boundary Condition Works
+end
+
+@testset "Splitter Flow Test" begin
+    # Ensure FlowDivider Splits Flow Properly
+    # 1) Set flow into port A, expect reduction in port B
+
+    # 2) Set flow into port B, expect increase in port B
+end
+
 #TODO: Test Valve Inversion