diff --git a/src/Hydraulic/IsothermalCompressible/components.jl b/src/Hydraulic/IsothermalCompressible/components.jl index be891ee4c..56228c8a7 100644 --- a/src/Hydraulic/IsothermalCompressible/components.jl +++ b/src/Hydraulic/IsothermalCompressible/components.jl @@ -53,7 +53,7 @@ Provides an "open" boundary condition for a hydraulic port such that mass flow ` end """ - TubeBase(add_inertia = true; area, length_int, head_factor = 1, perimeter = 2 * sqrt(area * pi), shape_factor = 64, name) + TubeBase(add_inertia = true, variable_length = 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. @@ -89,7 +89,7 @@ Variable length internal flow model of the fully developed incompressible flow f @variables begin x(t), [guess = length_int] - ddm(t), [guess = 0] + ddm(t) = 0 end vars = [] @@ -119,7 +119,7 @@ Variable length internal flow model of the fully developed incompressible flow f f = friction_factor(dm, area, d_h, μ, shape_factor) u = dm / (ρ * area) - shear = (1 / 2) * ρ * u * abs(u) * f * head_factor * (c / d_h) + shear = (1 / 2) * ρ * regPow(u, 2) * f * head_factor * (c / d_h) inertia = if add_inertia (c / area) * ddm else @@ -160,21 +160,10 @@ Constant length internal flow model discretized by `N` (`FixedVolume`: `N`, `Tub """ @component function Tube(N, add_inertia = true; area, length, head_factor = 1, perimeter = 2 * sqrt(area * pi), - shape_factor = 64, name) + shape_factor = 64, p_int, name) @assert(N>0, "the Tube component must be defined with at least 1 segment (i.e. N>0), found N=$N") - if N == 1 - return TubeBase(add_inertia, - false; - shape_factor, - area, - length_int = length, - head_factor, - perimeter, - name) - end - #TODO: How to set an assert effective_length >= length ?? pars = @parameters begin area = area @@ -182,6 +171,7 @@ Constant length internal flow model discretized by `N` (`FixedVolume`: `N`, `Tub head_factor = head_factor perimeter = perimeter shape_factor = shape_factor + p_int = p_int end vars = [] @@ -192,33 +182,30 @@ Constant length internal flow model discretized by `N` (`FixedVolume`: `N`, `Tub end pipe_bases = [] - for i in 1:(N - 1) - x = TubeBase(add_inertia; name = Symbol("p$i"), + for i in 1:N + x = TubeBase(add_inertia, false; name = Symbol("p$i"), shape_factor = ParentScope(shape_factor), area = ParentScope(area), - length_int = ParentScope(length) / (N - 1), + length_int = N > 1 ? ParentScope(length) / (N - 1) : ParentScope(length), head_factor = ParentScope(head_factor), perimeter = ParentScope(perimeter)) push!(pipe_bases, x) end + + eqs = [ + connect(pipe_bases[1].port_a, port_a) + connect(pipe_bases[end].port_b, port_b) + ] + volumes = [] - for i in 1:N + for i in 1:(N-1) x = FixedVolume(; name = Symbol("v$i"), - vol = ParentScope(area) * ParentScope(length) / N) + vol = ParentScope(area) * ParentScope(length) / (N-1), + 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 in 2:(N - 1) push!(eqs, - connect(volumes[i].port, pipe_bases[i - 1].port_b, pipe_bases[i].port_a)) - end - - for i in 1:(N - 1) - push!(eqs, pipe_bases[i].x ~ length / (N - 1)) + connect(x.port, pipe_bases[i].port_b, pipe_bases[i+1].port_a)) end return ODESystem(eqs, t, vars, pars; name, systems = [ports; pipe_bases; volumes]) @@ -281,8 +268,8 @@ end end vars = @variables begin - area(t) - y(t) + area(t), [guess=0] + y(t), [guess=0] end # let @@ -306,7 +293,7 @@ end eqs = [0 ~ port_a.dm + port_b.dm domain_connect(port_a, port_b) - dm ~ regRoot(2 * Δp * ρ / c) * x + dm ~ regRoot(2 * Δp * ρ / c) * x # I think this should be reformulated as: regRoot(2 DP rho) c x y ~ x] ODESystem(eqs, t, vars, pars; name, systems) @@ -357,11 +344,12 @@ Valve with `area` input and discharge coefficient `Cd` defined by https://en.wik ODESystem(eqs, t, vars, pars; name, systems) end -@component function VolumeBase(; area, dead_volume = 0, Χ1 = 1, Χ2 = 1, +@component function VolumeBase(; area, dead_volume = 0, p_int, x_int, name) pars = @parameters begin area = area dead_volume = dead_volume + p_int = p_int end systems = @named begin @@ -369,11 +357,11 @@ end end vars = @variables begin - x(t) + x(t) = x_int dx(t), [guess = 0] rho(t), [guess = liquid_density(port)] - drho(t), [guess = 0] - vol(t) + m(t), [guess = 0] + vol(t) end # let @@ -382,11 +370,14 @@ end eqs = [vol ~ dead_volume + area * x D(x) ~ dx - D(rho) ~ drho - rho ~ full_density(port, p) - dm ~ drho * vol * Χ1 + rho * area * dx * Χ2] + D(m) ~ dm + #rho ~ full_density(port, p) + p ~ full_pressure(port, rho) # see https://github.com/SciML/OrdinaryDiffEq.jl/issues/2561 + m ~ rho * vol] - ODESystem(eqs, t, vars, pars; name, systems) + initialization_eqs = [p ~ p_int] + + ODESystem(eqs, t, vars, pars; name, systems, initialization_eqs) end """ @@ -400,9 +391,10 @@ Fixed fluid volume. # Connectors: - `port`: hydraulic port """ -@component function FixedVolume(; vol, name) +@component function FixedVolume(; vol, name, p_int) pars = @parameters begin vol = vol + p_int = p_int end systems = @named begin @@ -411,16 +403,19 @@ Fixed fluid volume. vars = @variables begin rho(t), [guess = liquid_density(port)] - drho(t), [guess = 0] + m(t), [guess=vol*liquid_density(port)] + p(t)=p_int end # let dm = port.dm - p = port.p + - eqs = [D(rho) ~ drho - rho ~ full_density(port, p) - dm ~ drho * vol] + eqs = [D(m) ~ dm + # rho ~ full_density(port, p) + p ~ full_pressure(port, rho) # see https://github.com/SciML/OrdinaryDiffEq.jl/issues/2561 + p ~ port.p + m ~ rho * vol] ODESystem(eqs, t, vars, pars; name, systems) end @@ -467,24 +462,32 @@ See also [`FixedVolume`](@ref), [`DynamicVolume`](@ref) #parameters area, - direction = +1, name) + direction = +1, + x_int, + name) + pars = @parameters begin area = area + x_int = x_int end vars = @variables begin - x(t) - dx(t) - p(t) - f(t) - rho(t) - drho(t) - dm(t) + x(t)=x_int + dx(t), [guess=0] + p(t), [guess=0] + f(t), [guess=0] + rho(t), [guess=0] + m(t), [guess=0] + dm(t), [guess=0] end systems = @named begin port = HydraulicPort() flange = MechanicalPort() + damper = ValveBase(reversible; + Cd, + Cd_reverse, + minimum_area) end eqs = [ @@ -496,18 +499,20 @@ See also [`FixedVolume`](@ref), [`DynamicVolume`](@ref) # differentials D(x) ~ dx - D(rho) ~ drho + D(m) ~ dm # physics - rho ~ liquid_density(port, p) + # rho ~ full_density(port, p) + p ~ full_pressure(port, rho) # see https://github.com/SciML/OrdinaryDiffEq.jl/issues/2561 + f ~ p * area - dm ~ drho * x * area + rho * dx * area] + m ~ rho * x * area] - ODESystem(eqs, t, vars, pars; name, systems, defaults = [rho => liquid_density(port)]) + ODESystem(eqs, t, vars, pars; name, systems) end """ - DynamicVolume(N, add_inertia=true; p_int, area, x_int = 0, x_max, x_min = 0, x_damp = x_min, direction = +1, perimeter = 2 * sqrt(area * pi), shape_factor = 64, head_factor = 1, Cd = 1e2, Cd_reverse = Cd, name) + DynamicVolume(; p_int, area, x_int = 0, x_max, x_min = 0, x_damp = x_min, direction = +1, perimeter = 2 * sqrt(area * pi), shape_factor = 64, head_factor = 1, Cd = 1e2, Cd_reverse = Cd, name) Volume with moving wall with `flange` connector for converting hydraulic energy to 1D mechanical. The `direction` argument aligns the mechanical port with the hydraulic port, useful when connecting two dynamic volumes together in oppsing directions to create an actuator. @@ -524,7 +529,6 @@ dm ────► │ │ area ``` # Features: -- volume discretization with flow resistance and inertia: use `N` to control number of volume and resistance elements. Set `N=0` to turn off volume discretization. See `TubeBase` for more information about flow resistance. - minimum volume flow shutoff with damping and directional resistance. Use `reversible=false` when problem defines volume position `x` and solves for `dm` to prevent numerical instability. # Parameters: @@ -550,7 +554,7 @@ dm ────► │ │ area - `port`: hydraulic port - `flange`: mechanical translational port """ -@component function DynamicVolume(N, add_inertia = true, reversible = false; +@component function DynamicVolume(reversible = false; area, x_int = 0, x_max, @@ -562,14 +566,18 @@ dm ────► │ │ area perimeter = 2 * sqrt(area * pi), shape_factor = 64, head_factor = 1, + p_int, # Valve Cd = 1e2, Cd_reverse = Cd, minimum_area = 0, + + # Damping + d = 0, + name) - @assert(N>=0, - "the Tube component must be defined with 0 or more segments (i.e. N>=0), found N=$N") + @assert (direction == +1)||(direction == -1) "direction argument must be +/-1, found $direction" #TODO: How to set an assert effective_length >= length ?? @@ -584,42 +592,34 @@ dm ────► │ │ area perimeter = perimeter shape_factor = shape_factor head_factor = head_factor + p_int = p_int Cd = Cd Cd_reverse = Cd_reverse minimum_area = minimum_area + + d = d end - vars = @variables x(t)=x_int vol(t)=x_int * area + vars = @variables begin + x(t)=x_int + vol(t), [guess=x_int * area] + end - ports = @named begin + systems = @named begin port = HydraulicPort(;) flange = MechanicalPort(;) damper = ValveBase(reversible; Cd, Cd_reverse, minimum_area) + moving_volume = VolumeBase(; + area, + dead_volume = area * x_int, + p_int, + x_int=0) end - pipe_bases = [] - for i in 1:N - comp = TubeBase(add_inertia; name = Symbol("p$i"), - shape_factor = ParentScope(shape_factor), - area = ParentScope(area), - length_int = 0, #set in equations - head_factor = ParentScope(head_factor), - perimeter = ParentScope(perimeter)) - push!(pipe_bases, comp) - end - - #TODO: How to handle x_int? - #TODO: Handle direction - @named moving_volume = VolumeBase(; - area, - dead_volume = N == 0 ? area * x_int : 0, - Χ1 = N == 0 ? 1 : 0, - Χ2 = 1) # changed x_int to x_min - ratio = (x - x_min) / (x_damp - x_min) damper_area = if reversible @@ -628,55 +628,20 @@ dm ────► │ │ area ifelse(x >= x_damp, one(x), ifelse((x < x_damp) & (x > x_min), ratio, zero(x))) end + dx = moving_volume.dx + p = moving_volume.port.p + eqs = [vol ~ x * area D(x) ~ flange.v * direction damper.area ~ damper_area - connect(port, damper.port_b)] - - volumes = [] - if N > 0 - Δx = ParentScope(x_max) / N - x₀ = ParentScope(x_int) - - for i in 1:N - length = ifelse(x₀ > Δx * i, - Δx, - ifelse(x₀ - Δx * (i - 1) > 0, - x₀ - Δx * (i - 1), - zero(Δx))) - - comp = VolumeBase(; name = Symbol("v$i"), - area = ParentScope(area), - dead_volume = ParentScope(area) * length, Χ1 = 1, Χ2 = 0) - - push!(volumes, comp) - end - - push!(eqs, connect(moving_volume.port, volumes[end].port, pipe_bases[end].port_a)) - push!(eqs, connect(pipe_bases[1].port_b, damper.port_a)) - for i in 1:(N - 1) - push!(eqs, - connect(volumes[i].port, pipe_bases[i + 1].port_b, pipe_bases[i].port_a)) - end - - for i in 1:N - push!(eqs, - volumes[i].dx ~ ifelse( - (vol >= (i - 1) * (x_max / N) * area) & - (vol < (i) * (x_max / N) * area), - flange.v * direction, 0)) - push!(eqs, pipe_bases[i].x ~ volumes[i].vol / volumes[i].area) - end - else - push!(eqs, connect(moving_volume.port, damper.port_a)) - end + connect(port, damper.port_b) + connect(moving_volume.port, damper.port_a) + dx ~ flange.v * direction + p * area - dx*d ~ -flange.f * direction + ] - push!(eqs, moving_volume.dx ~ flange.v * direction) - push!(eqs, -moving_volume.port.p * area * direction ~ flange.f) - ODESystem(eqs, t, vars, pars; name, - systems = [ports; pipe_bases; volumes; moving_volume], - defaults = [flange.v => 0]) + return ODESystem(eqs, t, vars, pars; name, systems) end """ @@ -698,10 +663,11 @@ Spool valve with `x` valve opening input as mechanical flange port and `d` diame See [`Valve`](@ref) for more information. """ -@component function SpoolValve(reversible = false; Cd, d, name) +@component function SpoolValve(reversible = false; Cd, d, x_int, name) pars = @parameters begin d = d Cd = Cd + x_int = x_int end systems = @named begin @@ -712,8 +678,8 @@ See [`Valve`](@ref) for more information. end vars = @variables begin - x(t) - dx(t) + x(t)=x_int + dx(t), [guess=0] end eqs = [D(x) ~ dx @@ -751,7 +717,7 @@ end See [`SpoolValve`](@ref) for more information. """ -@component function SpoolValve2Way(reversible = false; m, g, Cd, d, name) +@component function SpoolValve2Way(reversible = false; m, g, Cd, d, x_int, dx_int, name) pars = @parameters begin m = m g = g @@ -759,13 +725,16 @@ See [`SpoolValve`](@ref) for more information. d = d Cd = Cd + + x_int = x_int + dx_int = dx_int end vars = [] systems = @named begin - vSA = SpoolValve(reversible; Cd, d) - vBR = SpoolValve(reversible; Cd, d) + vSA = SpoolValve(reversible; Cd, d, x_int) + vBR = SpoolValve(reversible; Cd, d, x_int) port_s = HydraulicPort(;) port_a = HydraulicPort(;) @@ -782,8 +751,13 @@ See [`SpoolValve`](@ref) for more information. connect(vBR.port_a, port_b) connect(vBR.port_b, port_r) connect(vSA.flange, vBR.flange, mass.flange, flange)] + + initialization_eqs = [ + mass.s ~ x_int + # mass.v ~ dx_int + ] - ODESystem(eqs, t, vars, pars; name, systems) + ODESystem(eqs, t, vars, pars; name, systems, initialization_eqs) end """ @@ -853,7 +827,7 @@ Actuator made of two DynamicVolumes connected in opposite direction with body ma - `port_b`: hydraulic port - `flange`: mechanical translational port """ -@component function Actuator(N, add_inertia = true, reversible = false; +@component function Actuator(reversible = false; area_a, area_b, perimeter_a = 2 * sqrt(area_a * pi), @@ -867,12 +841,16 @@ Actuator made of two DynamicVolumes connected in opposite direction with body ma m, g, x_int = 0, + dx_int = 0, minimum_volume_a = 0, minimum_volume_b = 0, damping_volume_a = minimum_volume_a, damping_volume_b = minimum_volume_b, Cd = 1e4, Cd_reverse = Cd, + d=0, + p_a_int, + p_b_int, name) pars = @parameters begin area_a = area_a @@ -884,6 +862,7 @@ Actuator made of two DynamicVolumes connected in opposite direction with body ma head_factor_a = head_factor_a head_factor_b = head_factor_b x_int = x_int + dx_int = dx_int length_a_int = length_a_int length_b_int = length_b_int minimum_volume_a = minimum_volume_a @@ -894,18 +873,21 @@ Actuator made of two DynamicVolumes connected in opposite direction with body ma Cd_reverse = Cd_reverse m = m g = g + d = d + p_a_int = p_a_int + p_b_int = p_b_int end vars = @variables begin x(t) = x_int - dx(t) + dx(t) = dx_int end total_length = length_a_int + length_b_int #TODO: include effective_length systems = @named begin - vol_a = DynamicVolume(N, add_inertia, reversible; direction = +1, + vol_a = DynamicVolume(reversible; direction = +1, area = area_a, x_int = length_a_int, x_max = total_length, @@ -915,9 +897,11 @@ Actuator made of two DynamicVolumes connected in opposite direction with body ma shape_factor = shape_factor_a, head_factor = head_factor_a, Cd, - Cd_reverse) + Cd_reverse, + d, + p_int = p_a_int) - vol_b = DynamicVolume(N, add_inertia, reversible; direction = -1, + vol_b = DynamicVolume(reversible; direction = -1, area = area_b, x_int = length_b_int, x_max = total_length, @@ -927,7 +911,9 @@ Actuator made of two DynamicVolumes connected in opposite direction with body ma shape_factor = shape_factor_b, head_factor = head_factor_b, Cd, - Cd_reverse) + Cd_reverse, + d, + p_int = p_b_int) mass = Mass(; m, g) port_a = HydraulicPort() port_b = HydraulicPort() @@ -940,5 +926,57 @@ Actuator made of two DynamicVolumes connected in opposite direction with body ma D(x) ~ dx dx ~ vol_a.flange.v] - ODESystem(eqs, t, vars, pars; name, systems) + initialization_eqs = [ + mass.s ~ x_int + ] + + ODESystem(eqs, t, vars, pars; name, systems, initialization_eqs) +end + +""" + Orifice() + +A valve in fixed position, with parameters for area and the discharge coefficient (fitting the form Effective Area = area x Cd) + +``` + ┌ + │ + ▲ +dm ────► effective area + ▼ + │ + └ +``` + +# Features: +- + +# Parameters: +## volume +- `area`: [m^2] physical area +- `cd`: [unitless] discharge coefficient + +# Connectors: +- `port_a`: hydraulic port +- `port_b`: hydraulic port +""" + +@mtkmodel Orifice begin + @parameters begin + orifice_area=0.00094 + Cd=0.6 # TODO Cd here is defined differently from Valve(). + # Here it follows the form Effective Orifice Area = Cd x Physical Orifice Area + # The Valve component should be updated too. + end + @components begin + area = Constant(k=orifice_area) + valve = Valve(Cd= 1 / (Cd * Cd)) + port_a = HydraulicPort() + port_b = HydraulicPort() + end + @equations begin + connect(valve.area, area.output) + connect(valve.port_a, port_a) + connect(valve.port_b, port_b) + end end diff --git a/src/Hydraulic/IsothermalCompressible/utils.jl b/src/Hydraulic/IsothermalCompressible/utils.jl index 425686d1f..2ba047552 100644 --- a/src/Hydraulic/IsothermalCompressible/utils.jl +++ b/src/Hydraulic/IsothermalCompressible/utils.jl @@ -1,6 +1,5 @@ -import ChainRulesCore - -regPow(x, a, delta = 0.01) = x * (x * x + delta * delta)^((a - 1) / 2); +# regPow(x, a, delta = 0.01) = x * (x * x + delta * delta)^((a - 1) / 2); +regPow(x, a, delta = 0.01) = ifelse(abs(x/delta) >= 1, sign(x)*abs(x/delta)^a * delta^a, (delta^a*x)/delta) regRoot(x, delta = 0.01) = regPow(x, 0.5, delta) """ @@ -74,6 +73,13 @@ Fluid parameter setter for isothermal compressible fluid domain. Defaults given ODESystem(eqs, t, vars, pars; name) end + +function transition(x1, x2, y1, y2, x) + u = (x - x1) / (x2 - x1) + blend = u^2 * (3 - 2 * u) + return (1 - blend) * y1 + blend * y2 +end + f_laminar(shape_factor, Re) = shape_factor * regPow(Re, -1, 0.1) #regPow used to avoid dividing by 0, min value is 0.1 f_turbulent(shape_factor, Re) = (shape_factor / 64) / (0.79 * log(Re) - 1.64)^2 @@ -118,15 +124,6 @@ end @register_symbolic friction_factor(dm, area, d_h, viscosity, shape_factor) Symbolics.derivative(::typeof(friction_factor), args, ::Val{1}) = 0 Symbolics.derivative(::typeof(friction_factor), args, ::Val{4}) = 0 -function ChainRulesCore.frule(_, ::typeof(friction_factor), args...) - (friction_factor(args...), ChainRulesCore.ZeroTangent()) -end - -function transition(x1, x2, y1, y2, x) - u = (x - x1) / (x2 - x1) - blend = u^2 * (3 - 2 * u) - return (1 - blend) * y1 + blend * y2 -end density_ref(port) = port.ρ density_exp(port) = port.n @@ -141,19 +138,41 @@ function liquid_density(port, p) end #Tait-Murnaghan equation of state liquid_density(port) = liquid_density(port, port.p) -function liquid_density_differential(port, dp) - density_ref(port) * dp / bulk_modulus(port) +# p = beta*(rho/rho_0 - 1) +# (p/beta + 1)*rho_0 = rho + +function liquid_pressure(port, rho) + + (rho/density_ref(port) - 1)*bulk_modulus(port) + end + function gas_density(port, p) slope = (density_ref(port) - gas_density_ref(port)) / (0 - gas_pressure_ref(port)) b = density_ref(port) return b + p * slope end + +function gas_pressure(port, rho) + + slope = (0 - gas_pressure_ref(port)) / (density_ref(port) - gas_density_ref(port)) + b = 0 + + return b + rho * slope +end + function full_density(port, p) ifelse(port.let_gas == 1, - ifelse(p > 0, liquid_density(port, p), gas_density(port, p)), + ifelse(p >= 0, liquid_density(port, p), gas_density(port, p)), liquid_density(port, p)) end full_density(port) = full_density(port, port.p) + +function full_pressure(port, rho) + ifelse(port.let_gas == 1, + ifelse( rho >= density_ref(port), liquid_pressure(port, rho), gas_pressure(port, rho)), + liquid_pressure(port, rho) + ) +end \ No newline at end of file diff --git a/test/Hydraulic/isothermal_compressible.jl b/test/Hydraulic/isothermal_compressible.jl index 6d79c833d..e4bda0652 100644 --- a/test/Hydraulic/isothermal_compressible.jl +++ b/test/Hydraulic/isothermal_compressible.jl @@ -9,72 +9,75 @@ using ModelingToolkitStandardLibrary.Blocks: Parameter NEWTON = NLNewton( check_div = false, always_new = true, max_iter = 100, relax = 9 // 10, κ = 1e-6) -#TODO: add initialization test for N=5 @testset "Fluid Domain and Tube" begin function System(N; bulk_modulus, name) pars = @parameters begin bulk_modulus = bulk_modulus + p_int = 0 end systems = @named begin fluid = IC.HydraulicFluid(; bulk_modulus) - stp = B.Step(; height = 10e5, offset = 0, start_time = 0.005, duration = Inf, - smooth = true) - src = IC.Pressure(;) - vol = IC.FixedVolume(; vol = 10.0) - res = IC.Tube(N; area = 0.01, length = 50.0) + stp = B.Step(; height = 10e5, offset = 0, start_time = 0.005, duration = Inf, smooth = true) + src = IC.Pressure() + vol = IC.FixedVolume(; vol = 10.0, p_int) + res = IC.Tube(N; area = 0.01, length = 50.0, p_int) end - eqs = [connect(stp.output, src.p) + eqs = [ + connect(stp.output, src.p) connect(fluid, src.port) + connect(src.port, res.port_a) - connect(res.port_b, vol.port)] + connect(res.port_b, vol.port) + + ] ODESystem(eqs, t, [], pars; name, systems) end - @named sys1_2 = System(1; bulk_modulus = 1e9) - @named sys1_1 = System(1; bulk_modulus = 1e7) - @named sys5_1 = System(5; bulk_modulus = 1e9) + @mtkbuild s1_1 = System(1; bulk_modulus = 1e9) + @mtkbuild s1_2 = System(1; bulk_modulus = 2e9) + @mtkbuild s5_1 = System(5; bulk_modulus = 1e9) - syss = structural_simplify.([sys1_2, sys1_1]) #removed sys5_1 for now - probs = [ODEProblem(sys, [], (0, 0.05); - initialization_eqs = [sys.vol.port.p ~ 0, sys.res.port_a.dm ~ 0]) - for sys in syss] # - sols = [solve(prob, Rodas5P()) - for prob in probs] + + + p1_1 = ODEProblem(s1_1, [], (0, 0.05)) + p1_2 = ODEProblem(s1_2, [], (0, 0.05)) + p5_1 = ODEProblem(s5_1, [], (0, 0.05)) - s1_2 = complete(sys1_2) - s1_1 = complete(sys1_1) - s5_1 = complete(sys5_1) + sol1_1 = solve(p1_1, Rodas5P()) + sol1_2 = solve(p1_2, Rodas5P()) + sol5_1 = solve(p5_1, Rodas5P()) + # fig = Figure() + # tm = 0:0.001:0.05 |> collect + # ax = Axis(fig[1,1]) + # lines!(ax, tm, sol1_1.(tm; idxs=s1_2.vol.port.p)); fig + # lines!(ax, tm, sol1_2.(tm; idxs=s1_1.vol.port.p)); fig + # lines!(ax, tm, sol5_1.(tm; idxs=s5_1.vol.port.p)); fig + # fig + # 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 sol1_2[s1_2.vol.port.p][end] > sol1_1[s1_1.vol.port.p][end] - #TODO: bring back after implementing N=5 # 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] + @test sol1_1[s1_1.vol.port.p][end] > sol5_1[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 + function System(; name) pars = [] systems = @named begin fluid = IC.HydraulicFluid() sink = IC.FixedPressure(; p = 10e5) - vol = IC.FixedVolume(; vol = 0.1) + vol = IC.FixedVolume(; vol = 5, p_int=1e5) valve = IC.Valve(; Cd = 1e5, minimum_area = 0) - ramp = B.Ramp(; height = 1, duration = 0.001, offset = 0, start_time = 0.001, - smooth = true) + ramp = B.Ramp(; height = 0.1, duration = 0.1, offset = 0, start_time = 0.1, smooth = true) end eqs = [connect(fluid, sink.port) @@ -86,23 +89,23 @@ end end @named valve_system = System() - sys = structural_simplify(valve_system) - initialization_eqs = [sys.vol.port.p ~ 101325] - initsys = ModelingToolkit.generate_initializesystem(sys; initialization_eqs) - # initsys = structural_simplify(initsys) - # initprob = NonlinearProblem(initsys, [t=>0]) - # initsol = solve(initprob) - prob = ODEProblem(sys, [], (0, 0.01)) - sol = solve(prob, Rodas5P()) + sys = structural_simplify(valve_system) + prob = ODEProblem(sys, [], (0, 1)) + sol = solve(prob, Rodas5P(); abstol=1e-6, reltol=1e-9) s = complete(valve_system) # the volume should discharge to 10bar @test sol[s.vol.port.p][end]≈10e5 atol=1e5 + + # fig = Figure() + # tm = 0:0.01:1 |> collect + # ax = Axis(fig[1,1]) + # lines!(ax, tm, sol.(tm; idxs=sys.vol.port.p)); + # fig end -#TODO: implement initialization system, currently seems to have an issue with DynamicVolume and issue: https://github.com/SciML/ModelingToolkit.jl/issues/2952 @testset "DynamicVolume and minimum_volume feature" begin # Need help here - function System(N; name, area = 0.01, length = 0.1, damping_volume) + function System(; name, area = 0.01, length = 0.1, damping_volume=length*area*0.1) pars = [] # DynamicVolume values @@ -112,24 +115,30 @@ end src1 = IC.Pressure(;) src2 = IC.Pressure(;) - vol1 = IC.DynamicVolume(N; direction = +1, + vol1 = IC.DynamicVolume(; direction = +1, area, x_int = length, x_max = length * 2, x_min = length * 0.1, - x_damp = damping_volume / area + length * 0.1) + x_damp = damping_volume / area + length * 0.1, + d = 1e3, + p_int = 10e5) + # vol1 = IC.Volume(;area, direction = +1, x_int=length) - vol2 = IC.DynamicVolume(N; direction = -1, + vol2 = IC.DynamicVolume(; direction = -1, area, x_int = length, x_max = length * 2, x_min = length * 0.1, - x_damp = damping_volume / area + length * 0.1) + x_damp = damping_volume / area + length * 0.1, + d = 1e3, + p_int = 10e5) + # vol2 = IC.Volume(;area, direction = -1, x_int=length) mass = T.Mass(; m = 10) - sin1 = B.Sine(; frequency = 0.5, amplitude = +0.5e5, offset = 10e5) - sin2 = B.Sine(; frequency = 0.5, amplitude = -0.5e5, offset = 10e5) + sin1 = B.Sine(; frequency = 0.5, amplitude = +1e5, offset = 10e5) + sin2 = B.Sine(; frequency = 0.5, amplitude = -1e5, offset = 10e5) end eqs = [connect(fluid, src1.port) @@ -140,88 +149,58 @@ end connect(src1.p, sin1.output) connect(src2.p, sin2.output)] - ODESystem(eqs, t, [], pars; name, systems) - end + initialization_eqs = [ + mass.s ~ 0.0 + mass.v ~ 0.0 + ] - for N in [1, 2] - for damping_volume in [0.01 * 0.1 * 0.25] - @named system = System(N; damping_volume) - s = complete(system) - sys = structural_simplify(system) - - u0 = ModelingToolkit.missing_variable_defaults(sys) |> Dict{Num, Num} - - u0[sys.vol1.x] = 0.1 - u0[sys.vol1.v1.rho] = IC.liquid_density(sys.fluid, 10e5) - u0[sys.vol1.v1.port.p] = 10e5 - u0[sys.vol1.damper.port_a.p] = 10e5 - - u0[sys.vol2.x] = 0.1 - u0[sys.vol2.v1.rho] = IC.liquid_density(sys.fluid, 10e5) - u0[sys.vol2.v1.port.p] = 10e5 - u0[sys.vol2.damper.port_a.p] = 10e5 - - if N == 2 - u0[sys.vol1.v2.rho] = IC.liquid_density(sys.fluid, 10e5) - u0[sys.vol1.v2.port.p] = 10e5 - u0[sys.vol2.v2.rho] = IC.liquid_density(sys.fluid, 10e5) - u0[sys.vol2.v2.port.p] = 10e5 - end - - prob = ODEProblem(sys, u0, (0, 5), - [s.vol1.Cd_reverse => 0.1, s.vol2.Cd_reverse => 0.1]; - jac = true) - - @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() - - # ax = Axis(fig[1,1], ylabel="position [m]", xlabel="time [s]") - # lines!(ax, sol.t, sol[s.vol1.x]; label="vol1") - # lines!(ax, sol.t, sol[s.vol2.x]; label="vol2") - # Legend(fig[1,2], ax) - - # ax = Axis(fig[2,1], ylabel="pressure [bar]", xlabel="time [s]") - # lines!(ax, sol.t, sol[s.vol1.damper.port_a.p]/1e5; label="vol1") - # lines!(ax, sol.t, sol[s.vol2.damper.port_a.p]/1e5; label="vol2") - # ylims!(ax, 10-0.5, 10+0.5) - - # ax = Axis(fig[3,1], ylabel="area", xlabel="time [s]") - # lines!(ax, sol.t, sol[s.vol1.damper.area]; label="area 1") - # lines!(ax, sol.t, sol[s.vol2.damper.area]; label="area 2") - - # display(fig) - # end - - i1 = round(Int, 1 / 1e-4) - i2 = round(Int, 2 / 1e-4) - i3 = round(Int, 3 / 1e-4) - i4 = round(Int, 4 / 1e-4) - - # volume/mass should stop moving at opposite ends - @test round(sol[s.vol1.x][1]; digits = 2) == 0.1 - @test round(sol[s.vol2.x][1]; digits = 2) == 0.1 - @test round(sol[s.vol1.x][i1]; digits = 2) == +0.19 - @test round(sol[s.vol2.x][i1]; digits = 2) == +0.01 - @test round(sol[s.vol1.x][i2]; digits = 2) == +0.01 - @test round(sol[s.vol2.x][i2]; digits = 2) == +0.19 - @test round(sol[s.vol1.x][i3]; digits = 2) == +0.19 - @test round(sol[s.vol2.x][i3]; digits = 2) == +0.01 - @test round(sol[s.vol1.x][i4]; digits = 2) == +0.01 - @test round(sol[s.vol2.x][i4]; digits = 2) == +0.19 - end + ODESystem(eqs, t, [], pars; name, systems, initialization_eqs) end + + @mtkbuild sys = System() + prob = ODEProblem(sys, [], (0, 5)) + sol = solve(prob, Rodas5P(); abstol=1e-6, reltol=1e-9) + + # begin + # fig = Figure() + + # ax = Axis(fig[1,1], ylabel="position [m]", xlabel="time [s]") + # lines!(ax, sol.t, sol[sys.vol1.x]; label="vol1") + # lines!(ax, sol.t, sol[sys.vol2.x]; label="vol2") + # Legend(fig[1,2], ax) + + # ax = Axis(fig[2,1], ylabel="pressure [bar]", xlabel="time [s]") + # lines!(ax, sol.t, sol[sys.vol1.damper.port_a.p]/1e5; label="vol1") + # lines!(ax, sol.t, sol[sys.vol2.damper.port_a.p]/1e5; label="vol2") + # ylims!(ax, 10-2, 10+2) + + # ax = Axis(fig[3,1], ylabel="area", xlabel="time [s]") + # lines!(ax, sol.t, sol[sys.vol1.damper.area]; label="area 1") + # lines!(ax, sol.t, sol[sys.vol2.damper.area]; label="area 2") + + # display(fig) + # end + + # volume/mass should stop moving at opposite ends + @test sol(0; idxs=sys.vol1.x) == 0.1 + @test sol(0; idxs=sys.vol2.x) == 0.1 + + @test round(sol(1; idxs=sys.vol1.x); digits=2) == 0.19 + @test round(sol(1; idxs=sys.vol2.x); digits=2) == 0.01 + + @test round(sol(2; idxs=sys.vol1.x); digits=2) == 0.01 + @test round(sol(2; idxs=sys.vol2.x); digits=2) == 0.19 + + @test round(sol(3; idxs=sys.vol1.x); digits=2) == 0.19 + @test round(sol(3; idxs=sys.vol2.x); digits=2) == 0.01 + + @test round(sol(4; idxs=sys.vol1.x); digits=2) == 0.01 + @test round(sol(4; idxs=sys.vol2.x); digits=2) == 0.19 + end -#TODO: implement initialization system, currently seems to have an issue with DynamicVolume and issue: https://github.com/SciML/ModelingToolkit.jl/issues/2952 @testset "Actuator System" begin - function System(use_input, f; name) + function System(use_input; name) pars = @parameters begin p_s = 200e5 p_r = 5e5 @@ -250,8 +229,8 @@ end systems = @named begin src = IC.FixedPressure(; p = p_s) - valve = IC.SpoolValve2Way(; g, m = m_f, d, Cd) - piston = IC.Actuator(0; + valve = IC.SpoolValve2Way(; g, m = m_f, d, Cd, x_int = 0, dx_int = 0) + piston = IC.Actuator(; length_a_int = l_1, length_b_int = l_2, area_a = A_1, @@ -261,14 +240,16 @@ end minimum_volume_a = A_1 * 1e-3, minimum_volume_b = A_2 * 1e-3, damping_volume_a = A_1 * 5e-3, - damping_volume_b = A_2 * 5e-3) - body = T.Mass(; m = 1500) - pipe = IC.Tube(1; area = A_2, length = 2.0) + damping_volume_b = A_2 * 5e-3, + p_a_int = p_1, + p_b_int = p_2) + # body = T.Mass(; m = 1500) + # pipe = IC.Tube(1; area = A_2, length = 2.0, p_int = p_2) snk = IC.FixedPressure(; p = p_r) pos = T.Position() - m1 = IC.FlowDivider(; n = 3) - m2 = IC.FlowDivider(; n = 3) + # m1 = IC.FlowDivider(; n = 3) + # m2 = IC.FlowDivider(; n = 3) fluid = IC.HydraulicFluid() end @@ -276,7 +257,8 @@ end if use_input @named input = B.SampledData(Float64) else - @named input = B.TimeVaryingFunction(f) + #@named input = B.TimeVaryingFunction(f) + @named input = B.Constant(k=0) end push!(systems, input) @@ -284,24 +266,38 @@ end eqs = [connect(input.output, pos.s) connect(valve.flange, pos.flange) connect(valve.port_a, piston.port_a) - connect(piston.flange, body.flange) - connect(piston.port_b, m1.port_a) - connect(m1.port_b, pipe.port_b) - connect(pipe.port_a, m2.port_b) - connect(m2.port_a, valve.port_b) + # connect(piston.flange, body.flange) + + connect(piston.port_b, valve.port_b) + + # connect(piston.port_b, pipe.port_b) + # # connect(piston.port_b, m1.port_a) + # # connect(m1.port_b, pipe.port_b) + + # connect(pipe.port_a, valve.port_b) + # # connect(pipe.port_a, m2.port_b) + # # connect(m2.port_a, valve.port_b) + connect(src.port, valve.port_s) connect(snk.port, valve.port_r) connect(fluid, src.port, snk.port) - D(body.v) ~ ddx] - - ODESystem(eqs, t, vars, pars; name, systems) + D(piston.mass.v) ~ ddx + ] + + initialization_eqs = [ + # body.s ~ 0 + ] + + ODESystem(eqs, t, vars, pars; name, systems, initialization_eqs) end - @named system = System(true, nothing) + @mtkbuild initsys = System(false) + + initprob = ODEProblem(initsys, [], (0,0)) + initsol = solve(initprob, Rodas5P()) + + @mtkbuild sys = System(true) - sys = structural_simplify(system) - defs = ModelingToolkit.defaults(sys) - s = complete(system) dt = 1e-4 time = 0:dt:0.1 @@ -309,47 +305,41 @@ end x = @. (time - 0.015)^2 - 10 * (time - 0.02)^3 x[1:150] = zeros(150) - defs[s.input.buffer] = Parameter(0.5 * x, dt) - - u0 = ModelingToolkit.missing_variable_defaults(sys) |> Dict{Num, Num} - - u0[sys.piston.vol_a.moving_volume.rho] = IC.liquid_density(sys.fluid, 45e5) - u0[sys.piston.vol_b.moving_volume.rho] = IC.liquid_density(sys.fluid, 45e5) - u0[sys.valve.vBR.valve.port_a.p] = 45e5 - u0[sys.piston.vol_a.moving_volume.port.p] = 45e5 - u0[sys.piston.vol_b.moving_volume.port.p] = 45e5 - u0[sys.piston.vol_a.damper.port_b.p] = 45e5 - u0[sys.piston.vol_b.damper.port_b.p] = 45e5 + defs = ModelingToolkit.defaults(sys) + defs[sys.input.buffer] = Parameter(0.5 * x, dt) - prob = ODEProblem(sys, u0, (0, 0.1); - tofloat = false) + # NOTE: bypassing initialization system: https://github.com/SciML/ModelingToolkit.jl/issues/3312 + prob = ODEProblem(sys, initsol[1], (0, 0.1); build_initializeprob = false) + + #TODO: Implement proper initialization system after issue is resolved + #TODO: How to bring the body back and not have an overdetermined system? # check the fluid domain - @test Symbol(defs[s.src.port.ρ]) == Symbol(s.fluid.ρ) - @test Symbol(defs[s.valve.port_s.ρ]) == Symbol(s.fluid.ρ) - @test Symbol(defs[s.valve.port_a.ρ]) == Symbol(s.fluid.ρ) - @test Symbol(defs[s.valve.port_b.ρ]) == Symbol(s.fluid.ρ) - @test Symbol(defs[s.valve.port_r.ρ]) == Symbol(s.fluid.ρ) - @test Symbol(defs[s.snk.port.ρ]) == Symbol(s.fluid.ρ) - - prob = remake(prob; tspan = (0, time[end])) - @time sol = solve(prob, ImplicitEuler(nlsolve = NEWTON); adaptive = false, dt, - initializealg = NoInit()) - + @test Symbol(defs[sys.src.port.ρ]) == Symbol(sys.fluid.ρ) + @test Symbol(defs[sys.valve.port_s.ρ]) == Symbol(sys.fluid.ρ) + @test Symbol(defs[sys.valve.port_a.ρ]) == Symbol(sys.fluid.ρ) + @test Symbol(defs[sys.valve.port_b.ρ]) == Symbol(sys.fluid.ρ) + @test Symbol(defs[sys.valve.port_r.ρ]) == Symbol(sys.fluid.ρ) + @test Symbol(defs[sys.snk.port.ρ]) == Symbol(sys.fluid.ρ) + + @time sol = solve(prob, Rodas5P(); initializealg = NoInit()); + @test sol[sys.ddx][1] == 0.0 @test maximum(sol[sys.ddx]) > 200 - @test sol[s.piston.x][end] > 0.6 + @test sol[sys.piston.x][end] > 0.6 + + end -#TODO: implement initialization system, currently seems to have an issue with DynamicVolume and issue: https://github.com/SciML/ModelingToolkit.jl/issues/2952 + @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(1; area = 0.001, x_int = 0.05, - x_max = 0.1, x_damp = 0.02, x_min = 0.01, direction = +1) + vol = IC.DynamicVolume(; area = 0.001, x_int = 0.05, + x_max = 0.1, x_damp = 0.02, x_min = 0.01, direction = +1, p_int=100e5) mass = T.Mass(; m = 100, g = -9.807) # s = 0.05 cap = IC.Cap() end @@ -357,46 +347,42 @@ 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) + initialization_eqs = [ + mass.s ~ 0.05 + mass.v ~ 0 + ] + - u0 = ModelingToolkit.missing_variable_defaults(sys) |> Dict{Num, Num} - - u0[sys.mass.s] = 0.05 - u0[sys.vol.v1.port.p] = 100e5 - u0[sys.vol.v1.rho] = IC.liquid_density(sys.fluid, 100e5) - u0[sys.vol.damper.port_b.p] = 100e5 + return ODESystem(eqs, t, [], pars; name, systems, initialization_eqs) + end - prob1 = ODEProblem(sys, u0, (0, 0.05)) - prob2 = ODEProblem(sys, u0, (0, 0.05), - [s.let_gas => 0]) + @mtkbuild sys = System() + + prob1 = ODEProblem(sys, [], (0, 0.05)) + # prob1 = remake(prob1; u0 = BigFloat.(prob1.u0)) + prob2 = ODEProblem(sys, [], (0, 0.05), [sys.let_gas => 0]) - @time sol1 = solve(prob1, ImplicitEuler(nlsolve = NEWTON); - adaptive = false, dt = 1e-4, initializealg = NoInit()) - @time sol2 = solve(prob2, ImplicitEuler(nlsolve = NEWTON); - adaptive = false, dt = 1e-4, initializealg = NoInit()) + # @time sol1 = solve(prob1, Rodas5P(); abstol=1e-9, reltol=1e-9) #BUG: Using BigFloat gives... ERROR: MethodError: no method matching getindex(::Missing, ::Int64) + @time sol1 = solve(prob1, Rodas5P(); adaptive=false, dt=1e-6) #TODO: fix BigFloat to implment abstol=1e-9, reltol=1e-9 + @time sol2 = solve(prob2, Rodas5P()) # 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 sol1[s.mass.s][end] > sol2[s.mass.s][end] + @test sol1[sys.mass.s][end] > sol2[sys.mass.s][end] # 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 + @test minimum(sol1[sys.vol.port.p]) > -5000 + @test minimum(sol2[sys.vol.port.p]) < -5000 # 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]) + # lines!(ax, sol1.t, sol1[sys.vol.port.p]); fig + # lines!(ax, sol2.t, sol2[sys.vol.port.p]); fig # ax = Axis(fig[1,2]) - # lines!(ax, sol1.t, sol1[s.mass.s]) - # lines!(ax, sol2.t, sol2[s.mass.s]) + # lines!(ax, sol1.t, sol1[sys.mass.s]) + # lines!(ax, sol2.t, sol2[sys.mass.s]) # fig end