Merge pull request #175 from SciML/bgc/negative_pressure

Bgc/negative pressure

docs/make.jl

@@ -7,6 +7,8 @@ using ModelingToolkitStandardLibrary.Magnetic
 using ModelingToolkitStandardLibrary.Magnetic.FluxTubes
 using ModelingToolkitStandardLibrary.Electrical
 using ModelingToolkitStandardLibrary.Thermal
+using ModelingToolkitStandardLibrary.Hydraulic
+using ModelingToolkitStandardLibrary.Hydraulic.IsothermalCompressible
 
 cp("./docs/Manifest.toml", "./docs/src/assets/Manifest.toml", force = true)
 cp("./docs/Project.toml", "./docs/src/assets/Project.toml", force = true)
@@ -33,7 +35,9 @@ makedocs(sitename = "ModelingToolkitStandardLibrary.jl",
              ModelingToolkitStandardLibrary.Magnetic,
              ModelingToolkitStandardLibrary.Magnetic.FluxTubes,
              ModelingToolkitStandardLibrary.Electrical,
-             ModelingToolkitStandardLibrary.Thermal],
+             ModelingToolkitStandardLibrary.Thermal,
+             ModelingToolkitStandardLibrary.Hydraulic,
+             ModelingToolkitStandardLibrary.Hydraulic.IsothermalCompressible],
          format = Documenter.HTML(assets = ["assets/favicon.ico"],
                                   canonical = "https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/"),
          pages = pages)

docs/pages.jl

@@ -14,6 +14,7 @@ pages = [
         "Magnetic Components" => "API/magnetic.md",
         "Mechanical Components" => "API/mechanical.md",
         "Thermal Components" => "API/thermal.md",
+        "Hydraulic Components" => "API/hydraulic.md",
         "Linear Analysis" => "API/linear_analysis.md",
     ],
 ]

docs/src/API/hydraulic.md

@@ -0,0 +1,44 @@
+# ModelingToolkit Standard Library: Hydrualic Components
+
+```@contents
+Pages = ["hydraulic.md"]
+Depth = 3
+```
+
+## Index
+
+```@index
+Pages = ["hydraulic.md"]
+```
+
+## IsothermalCompressible Components
+
+```@meta
+CurrentModule = ModelingToolkitStandardLibrary.Hydraulic.IsothermalCompressible
+```
+
+### IsothermalCompressible Utils
+
+```@docs
+HydraulicPort
+HydraulicFluid
+friction_factor
+```
+
+### IsothermalCompressible Components
+
+```@docs
+Cap
+TubeBase
+Tube
+FixedVolume
+DynamicVolume
+```
+
+### IsothermalCompressible Sources
+
+```@docs
+MassFlow
+Pressure
+FixedPressure
+```

docs/src/index.md

@@ -28,6 +28,7 @@ The following are the constituant libraries of the ModelingToolkit Standard Libr
   - [Electrical Components](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/API/electrical/)
   - [Magnetic Components](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/API/magnetic/)
   - [Thermal Components](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/API/thermal/)
+  - [Hydraulic Components](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/API/hydraulic/)
 
 ## Contributing
 

src/Hydraulic/IsothermalCompressible/components.jl

@@ -42,7 +42,11 @@ end
 """
     TubeBase(add_inertia = true; p_int, area, length_int, head_factor = 1, perimeter = 2 * sqrt(area * pi), shape_factor = 64, name)
 
-Variable length internal flow model of the fully developed flow friction, ignoring any compressibility.  Includes optional inertia equation when `add_inertia = true` to model wave propogation which includes change in flow and length terms.
+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.  
+
+# States:
+- `x`: [m] length of the pipe
+- `ddm`: [kg/s^2] Rate of change of mass flow rate in control volume.
 
 # Parameters:
 - `p_int`: [Pa] initial pressure
@@ -104,7 +108,7 @@ Variable length internal flow model of the fully developed flow friction, ignori
     end
 
     eqs = [0 ~ port_a.dm + port_b.dm
-           Δp ~ ifelse(x > 0, shear + inertia, 0)]
+           Δp ~ ifelse(x > 0, shear + inertia, zero(x))]
 
     if add_inertia
         push!(eqs, D(dm) ~ ddm)
@@ -116,7 +120,7 @@ end
 """
     Tube(N, add_inertia=true; p_int, area, length, head_factor=1, perimeter = 2 * sqrt(area * pi), shape_factor = 64, name)
 
-Constant length internal flow model with volume discretized by `N` which models the fully developed flow friction, compressibility, and inertia effects when `add_inertia = true`
+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
@@ -258,7 +262,10 @@ end
         port_b = HydraulicPort(; p_int = p_b_int)
     end
 
-    vars = @variables begin area(t) = area_int end
+    vars = @variables begin
+        area(t) = area_int
+        y(t) = area_int
+    end
 
     # let
     ρ = (full_density(port_a) + full_density(port_b)) / 2
@@ -275,7 +282,8 @@ end
     Cd = ifelse(Δp > 0, Cd, Cd_reverse)
 
     eqs = [0 ~ port_a.dm + port_b.dm
-           dm ~ regRoot(2 * Δp * ρ / Cd) * x]
+           dm ~ regRoot(2 * Δp * ρ / Cd) * x
+           y ~ x]
 
     ODESystem(eqs, t, vars, pars; name, systems)
 end
@@ -459,6 +467,7 @@ dm ────►               │  │ area
                                   # Valve
                                   Cd = 1e2,
                                   Cd_reverse = Cd,
+                                  minimum_area = 0,
                                   name)
     @assert(N>0,
             "the Tube component must be defined with more than 1 segment (i.e. N>1), found N=$N")
@@ -482,15 +491,16 @@ dm ────►               │  │ area
 
         Cd = Cd
         Cd_reverse = Cd_reverse
+        minimum_area = minimum_area
     end
 
     vars = @variables x(t)=x_int vol(t)=x_int * area
 
     ports = @named begin
         port = HydraulicPort(; p_int)
         flange = MechanicalPort()
-        damper = ValveBase(true; p_a_int = p_int, p_b_int = p_int, area_int = 1, Cd,
-                           Cd_reverse)
+        damper = ValveBase(; p_a_int = p_int, p_b_int = p_int, area_int = 1, Cd,
+                           Cd_reverse, minimum_area)
     end
 
     pipe_bases = []
@@ -519,7 +529,7 @@ dm ────►               │  │ area
                         Δx,
                         ifelse(x₀ - Δx * (i - 1) > 0,
                                x₀ - Δx * (i - 1),
-                               0))
+                               zero(Δx)))
 
         comp = VolumeBase(; name = Symbol("v$i"), p_int = ParentScope(p_int), x_int = 0,
                           area = ParentScope(area),
@@ -529,9 +539,11 @@ dm ────►               │  │ area
     end
 
     ratio = (x - x_min) / (x_damp - x_min)
-    damper_area = ifelse(x >= x_damp, 1,
-                         ifelse((x < x_damp) &
-                                (x > x_min), ratio, 0))
+    damper_area = ifelse(x >= x_damp,
+                         one(x),
+                         ifelse((x < x_damp) & (x > x_min),
+                                ratio,
+                                zero(x)))
 
     eqs = [vol ~ x * area
            D(x) ~ flange.v * direction

src/Hydraulic/IsothermalCompressible/utils.jl

@@ -20,6 +20,7 @@ Connector port for hydraulic components.
         β
         μ
         n
+        let_gas
         ρ_gas
         p_gas
     end
@@ -35,25 +36,27 @@ end
 """
     HydraulicFluid(; density=997, bulk_modulus=2.09e9, viscosity=0.0010016, name)
 
-Fluid parameter setter for isothermal compressible fluid domain.  Defaults given for water at 20°C and 0Pa gage (1atm absolute) reference pressure. Density is modeled using the Tait equation of state.  For pressures below the reference pressure, density is linearly interpolated to the gas state, this helps prevent pressures from going below the reference pressure.  
+Fluid parameter setter for isothermal compressible fluid domain.  Defaults given for water at 20°C and 0Pa gage (1atm absolute) reference pressure. Density is modeled using the Tait equation of state.  For pressures below the reference pressure, density is linearly interpolated to the gas state (when `let_gas` is set to 1), this helps prevent pressures from going below the reference pressure.  
 
 # Parameters:
 
 - `ρ`: [kg/m^3] fluid density at 0Pa reference gage pressure (set by `density` argument)
 - `Β`: [Pa] fluid bulk modulus describing the compressibility (set by `bulk_modulus` argument)
 - `μ`: [Pa*s] or [kg/m-s] fluid dynamic viscosity  (set by `viscosity` argument)
 - `n`: density exponent
+- `let_gas`: set to 1 to allow fluid to transition from liquid to gas (for density calculation only)  
 - `ρ_gas`: [kg/m^3] density of fluid in gas state at reference gage pressure `p_gas` (set by `gas_density` argument)
 - `p_gas`: [Pa] reference pressure (set by `gas_pressure` argument)
 """
 @connector function HydraulicFluid(; density = 997, bulk_modulus = 2.09e9,
                                    viscosity = 0.0010016, gas_density = 0.0073955,
-                                   gas_pressure = -1000, n = 1, name)
+                                   gas_pressure = -1000, n = 1, let_gas = 1, name)
     pars = @parameters begin
         ρ = density
         β = bulk_modulus
         μ = viscosity
         n = n
+        let_gas = let_gas
         ρ_gas = gas_density
         p_gas = gas_pressure
     end
@@ -96,7 +99,7 @@ function friction_factor(dm, area, d_h, density, viscosity, shape_factor)
     Re = density * u * d_h / viscosity
     f_laminar = shape_factor * regPow(Re, -1, 1e-6)
 
-    Re = maximum([Re, 1])
+    Re = max(Re, one(Re))
     f_turbulent = (shape_factor / 64) * (0.79 * log(Re) - 1.64)^(-2)
 
     f = transition(2000, 3000, f_laminar, f_turbulent, Re)
@@ -120,15 +123,26 @@ function transition(x1, x2, y1, y2, x)
 end
 
 density_ref(port) = port.ρ
+density_exp(port) = port.n
 gas_density_ref(port) = port.ρ_gas
 gas_pressure_ref(port) = port.p_gas
 bulk_modulus(port) = port.β
 viscosity(port) = port.μ
-liquid_density(port, p) = density_ref(port) * (1 + p / bulk_modulus(port))
+
+function liquid_density(port, p)
+    density_ref(port) *
+    regPow(1 + density_exp(port) * p / bulk_modulus(port), 1 / density_exp(port))
+end #Tait-Murnaghan equation of state
 liquid_density(port) = liquid_density(port, port.p)
 function gas_density(port, p)
-    density_ref(port) -
-    p * (density_ref(port) - gas_density_ref(port)) / gas_pressure_ref(port)
+    slope = (density_ref(port) - gas_density_ref(port)) / (0 - gas_pressure_ref(port))
+    b = density_ref(port)
+
+    return b + p * slope
+end
+function full_density(port, p)
+    ifelse(port.let_gas == 1,
+           ifelse(p > 0, liquid_density(port, p), gas_density(port, p)),
+           liquid_density(port, p))
 end
-full_density(port, p) = liquid_density(port, p)  #ifelse( p > 0, liquid_density(port, p), gas_density(port, p) )
 full_density(port) = full_density(port, port.p)

test/Hydraulic/isothermal_compressible.jl

@@ -8,7 +8,7 @@ using ModelingToolkitStandardLibrary.Blocks: Parameter
 @parameters t
 D = Differential(t)
 
-NEWTON = NLNewton(check_div = false, always_new = true, max_iter = 10, relax = 9 // 10)
+NEWTON = NLNewton(check_div = false, always_new = true, max_iter = 100, relax = 9 // 10)
 
 @testset "Fluid Domain and Tube" begin
     function System(N; bulk_modulus, name)
@@ -20,7 +20,7 @@ NEWTON = NLNewton(check_div = false, always_new = true, max_iter = 10, relax = 9
                          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 = 500.0)
+            res = IC.Tube(N; p_int = 0, area = 0.01, length = 50.0)
         end
 
         eqs = [connect(stp.output, src.p)
@@ -36,8 +36,8 @@ NEWTON = NLNewton(check_div = false, always_new = true, max_iter = 10, relax = 9
     @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.2))
-             for sys in syss]
+    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]
@@ -51,6 +51,15 @@ NEWTON = NLNewton(check_div = false, always_new = true, max_iter = 10, relax = 9
 
     # 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])
+    # lines!(ax, sols[2][s1_1.vol.port.p])
+    # lines!(ax, sols[3][s5_1.vol.port.p])
+    # fig
+
 end
 
 @testset "Valve" begin
@@ -131,11 +140,16 @@ end
             @named system = System(N; damping_volume)
             s = complete(system)
             sys = structural_simplify(system)
-            prob = ODEProblem(sys, [], (0, 5),
+            prob = ODEProblem(sys, ModelingToolkit.missing_variable_defaults(sys), (0, 5),
                               [s.vol1.Cd_reverse => 0.1, s.vol2.Cd_reverse => 0.1];
                               jac = true)
-            @time sol = solve(prob, ImplicitEuler(nlsolve = NEWTON); dt = 1e-4,
-                              adaptive = false, initializealg = NoInit())
+
+            @time sol = solve(prob,
+                              ImplicitEuler(nlsolve = NLNewton(check_div = false,
+                                                               always_new = true,
+                                                               max_iter = 10,
+                                                               relax = 9 // 10));
+                              dt = 0.0001, adaptive = false, initializealg = NoInit())
 
             # begin
             #     fig = Figure()
@@ -154,7 +168,7 @@ end
             #     lines!(ax, sol.t, sol[s.vol1.damper.area]; label="area 1")
             #     lines!(ax, sol.t, sol[s.vol2.damper.area]; label="area 2")
 
-            #     fig
+            #     display(fig)
             # end
 
             i1 = round(Int, 1 / 1e-4)
@@ -264,7 +278,8 @@ end
     sys = structural_simplify(system)
     defs = ModelingToolkit.defaults(sys)
     s = complete(system)
-    prob = ODEProblem(sys, [], (0, 0.1); tofloat = false, jac = true)
+    prob = ODEProblem(sys, ModelingToolkit.missing_variable_defaults(sys), (0, 0.1);
+                      tofloat = false, jac = true)
 
     # check the fluid domain
     @test Symbol(defs[s.src.port.ρ]) == Symbol(s.fluid.ρ)
@@ -291,4 +306,52 @@ end
     @test sol[s.piston.x][end]≈0.06 atol=0.01
 end
 
+@testset "Prevent Negative Pressure" begin
+    @component function System(; name)
+        pars = @parameters let_gas = 1
+
+        systems = @named begin
+            fluid = IC.HydraulicFluid(; let_gas)
+            vol = IC.DynamicVolume(5; p_int = 100e5, area = 0.001, x_int = 0.05,
+                                   x_max = 0.1, x_damp = 0.02, x_min = 0.01, direction = +1)
+            mass = T.Mass(; m = 100, g = -9.807, s_0 = 0.05)
+            cap = IC.Cap(; p_int = 100e5)
+        end
+
+        eqs = [connect(fluid, cap.port, vol.port)
+               connect(vol.flange, mass.flange)]
+
+        ODESystem(eqs, t, [], pars; name, systems)
+    end
+
+    @named system = System()
+    s = complete(system)
+    sys = structural_simplify(system)
+    prob1 = ODEProblem(sys, ModelingToolkit.missing_variable_defaults(sys), (0, 0.05))
+    prob2 = ODEProblem(sys, ModelingToolkit.missing_variable_defaults(sys), (0, 0.05),
+                       [s.let_gas => 0])
+
+    @time sol1 = solve(prob1, ImplicitEuler(nlsolve = NEWTON); adaptive = false, dt = 1e-4)
+    @time sol2 = solve(prob2, Rodas4())
+
+    # case 1: no negative pressure will only have gravity pulling mass back down
+    # case 2: with negative pressure, added force pulling mass back down
+    # - case 1 should push the mass higher
+    @test maximum(sol1[s.mass.s]) > maximum(sol2[s.mass.s])
+
+    # case 1 should prevent negative pressure less than -1000
+    @test minimum(sol1[s.vol.port.p]) > -1000
+    @test minimum(sol2[s.vol.port.p]) < -1000
+
+    # fig = Figure()
+    # ax = Axis(fig[1,1])
+    # lines!(ax, sol1.t, sol1[s.vol.port.p])
+    # lines!(ax, sol2.t, sol2[s.vol.port.p])
+
+    # ax = Axis(fig[1,2])
+    # lines!(ax, sol1.t, sol1[s.mass.s])
+    # lines!(ax, sol2.t, sol2[s.mass.s])
+    # fig
+end
+
 #TODO: Test Valve Inversion