diff --git a/docs/make.jl b/docs/make.jl index 00b9cfd7d..870e1c7b6 100644 --- a/docs/make.jl +++ b/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) diff --git a/docs/pages.jl b/docs/pages.jl index 0c8967580..af2de24f3 100644 --- a/docs/pages.jl +++ b/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", ], ] diff --git a/docs/src/API/hydraulic.md b/docs/src/API/hydraulic.md new file mode 100644 index 000000000..03f22d864 --- /dev/null +++ b/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 +``` diff --git a/docs/src/index.md b/docs/src/index.md index 296c41a1b..7416e8570 100644 --- a/docs/src/index.md +++ b/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 diff --git a/src/Hydraulic/IsothermalCompressible/components.jl b/src/Hydraulic/IsothermalCompressible/components.jl index b1084267d..04705682e 100644 --- a/src/Hydraulic/IsothermalCompressible/components.jl +++ b/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,6 +491,7 @@ dm ────► │ │ area Cd = Cd Cd_reverse = Cd_reverse + minimum_area = minimum_area end vars = @variables x(t)=x_int vol(t)=x_int * area @@ -489,8 +499,8 @@ dm ────► │ │ 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 diff --git a/src/Hydraulic/IsothermalCompressible/utils.jl b/src/Hydraulic/IsothermalCompressible/utils.jl index d3b1db216..223f4646a 100644 --- a/src/Hydraulic/IsothermalCompressible/utils.jl +++ b/src/Hydraulic/IsothermalCompressible/utils.jl @@ -20,6 +20,7 @@ Connector port for hydraulic components. β μ n + let_gas ρ_gas p_gas end @@ -35,7 +36,7 @@ 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: @@ -43,17 +44,19 @@ Fluid parameter setter for isothermal compressible fluid domain. Defaults given - `Β`: [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) diff --git a/test/Hydraulic/isothermal_compressible.jl b/test/Hydraulic/isothermal_compressible.jl index f5de4c3bc..fd67f5a30 100644 --- a/test/Hydraulic/isothermal_compressible.jl +++ b/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