add docs and tests for constant torque

Author: baggepinnen

src/Mechanical/Rotational/Rotational.jl

@@ -15,7 +15,7 @@ include("utils.jl")
 export Fixed, Inertia, Spring, Damper, IdealGear, RotationalFriction
 include("components.jl")
 
-export Torque, Speed
+export Torque, ConstantTorque, Speed
 include("sources.jl")
 
 export AngleSensor, SpeedSensor, TorqueSensor, RelSpeedSensor

src/Mechanical/Rotational/sources.jl

@@ -16,7 +16,7 @@ Input signal acting as external torque on a flange
 
 # States:
 
-  - `phi_support(t)`: [`rad`] Absolute angle of support flange"
+  - `phi_support(t)`: [`rad`] Absolute angle of support flange
 
 # Connectors:
 
@@ -35,6 +35,24 @@ function Torque(; name, use_support = false)
     return extend(ODESystem(eqs, t, [], []; name = name, systems = [tau]), partial_element)
 end
 
+"""
+    ConstantTorque(; name, tau_constant, use_support = false)
+
+Constant torque source
+
+# State variables:
+  
+- `phi_support(t)`: [`rad`] Absolute angle of support flange, only available if `use_support = true`
+- `tau`: Accelerating torque acting at flange (= -flange.tau)
+- `w`: Angular velocity of flange with respect to support (= der(phi))
+
+# Connectors:
+- `flange` [Flange](@ref)
+
+# Arguments:
+- `tau_constant`: The constant torque applied by the source
+- `use_support`: Whether or not an internal support flange is added, defaults to false.
+"""
 function ConstantTorque(; name, tau_constant, use_support = false)
     @named partial_element = PartialTorque(; use_support)
     @unpack flange, phi = partial_element

src/Mechanical/Rotational/utils.jl

@@ -18,26 +18,26 @@ Base.@doc """
     extend(ODESystem(Equation[], t, [], [], name = name), flange)
 end
 
-Base.@doc """
-    InternalSupport(;name, tau)
-
-1-dim. rotational flange of a shaft.
-
-- `tau`: External support torque (must be computed via torque balance in model where InternalSupport is used; = flange.tau)
-
-# States:
-- `phi(t)`: [`rad`] Absolute rotation angle of flange
-- `tau(t)`: [`N.m`] Cut torque in the flange
-""" Flange
-
-@connector function InternalSupport(; name, tau)
-    @named flange = Flange()
-    @variables phi(t)=0 [description = "Rotation angle of support $name"]
-    # tau(t), [connect = Flow, description = "Cut torque in support $name"],)
-    equations = [flange.tau ~ tau
-                 flange.phi ~ phi]
-    ODESystem(equations, t, [phi], [], name = name, systems = [flange]) # NOTE: tau not included since it belongs elsewhere
-end
+# Base.@doc """
+#     InternalSupport(;name, tau)
+
+# 1-dim. rotational flange of a shaft.
+
+# - `tau`: External support torque (must be computed via torque balance in model where InternalSupport is used; = flange.tau)
+
+# # States:
+# - `phi(t)`: [`rad`] Absolute rotation angle of flange
+# - `tau(t)`: [`N.m`] Cut torque in the flange
+# """ Flange
+
+# @connector function InternalSupport(; name, tau)
+#     @named flange = Flange()
+#     @variables phi(t)=0 [description = "Rotation angle of support $name"]
+#     # tau(t), [connect = Flow, description = "Cut torque in support $name"],)
+#     equations = [flange.tau ~ tau
+#                  flange.phi ~ phi]
+#     ODESystem(equations, t, [phi], [], name = name, systems = [flange]) # NOTE: tau not included since it belongs elsewhere
+# end
 
 Base.@doc """
     Support(;name)

test/Mechanical/rotational.jl

@@ -25,15 +25,15 @@ D = Differential(t)
 
     prob = ODEProblem(sys, Pair[], (0, 10.0))
     sol = solve(prob, Rodas4())
-    @test sol.retcode == Success
+    @test SciMLBase.successful_retcode(sol)
 
     prob = ODAEProblem(sys, Pair[], (0, 10.0))
     sol = solve(prob, Rodas4())
-    @test sol.retcode == Success
+    @test SciMLBase.successful_retcode(sol)
 
     prob = DAEProblem(sys, D.(states(sys)) .=> 0.0, Pair[], (0, 10.0))
     sol = solve(prob, DFBDF())
-    @test sol.retcode == Success
+    @test SciMLBase.successful_retcode(sol)
     @test all(sol[inertia1.w] .== 0)
     @test sol[inertia2.w][end]≈0 atol=1e-3 # all energy has dissipated
 
@@ -48,7 +48,7 @@ end
     @named fixed = Fixed()
     @named torque = Torque(use_support = true)
     @named inertia1 = Inertia(J = 2, phi_start = pi / 2)
-    @named spring = Spring(c = 1e4)
+    @named spring = Rotational.Spring(c = 1e4)
     @named damper = Damper(d = 10)
     @named inertia2 = Inertia(J = 4)
     @named sine = Blocks.Sine(amplitude = amplitude, frequency = frequency)
@@ -73,17 +73,37 @@ end
     prob = DAEProblem(sys, D.(states(sys)) .=> 0.0,
                       [D(D(inertia2.phi)) => 1.0; D.(states(model)) .=> 0.0], (0, 10.0))
     sol = solve(prob, DFBDF())
-    @test sol.retcode == Success
+    @test SciMLBase.successful_retcode(sol)
 
     prob = ODAEProblem(sys, Pair[], (0, 1.0))
     sol = solve(prob, Rodas4())
-    @test sol.retcode == Success
+    @test SciMLBase.successful_retcode(sol)
 
-    @test_skip begin
-        @test sol.retcode == Success
-        @test all(isapprox.(sol[inertia1.w], -sol[inertia2.w] * 2, atol = 1)) # exact opposite oscillation with smaller amplitude J2 = 2*J1
-        @test_broken all(sol[torque.flange.tau] .== -sol[sine.output.u]) # torque source is equal to negative sine
-    end
+    @test all(isapprox.(sol[inertia1.w], -sol[inertia2.w] * 2, atol = 1)) # exact opposite oscillation with smaller amplitude J2 = 2*J1
+    @test all(sol[torque.flange.tau] .== -sol[sine.output.u]) # torque source is equal to negative sine
+
+    ## Test with constant torque source
+    @named torque = ConstantTorque(use_support = true, tau_constant = 1)
+    connections = [connect(torque.support, fixed.flange)
+                   connect(torque.flange, inertia1.flange_a)
+                   connect(inertia1.flange_b, spring.flange_a, damper.flange_a)
+                   connect(spring.flange_b, damper.flange_b, inertia2.flange_a)]
+
+    @named model = ODESystem(connections, t,
+                             systems = [
+                                 fixed,
+                                 torque,
+                                 inertia1,
+                                 inertia2,
+                                 spring,
+                                 damper,
+                             ])
+    sys = structural_simplify(model)
+
+    prob = ODEProblem(sys, Pair[], (0, 10.0))
+    sol = solve(prob, Rodas4())
+    @test SciMLBase.successful_retcode(sol)
+    @test sol(sol.t[end], idxs = inertia1.w)≈sol(sol.t[end], idxs = inertia2.w) rtol=0.1 # both inertias have same angular velocity after initial transient
 end
 
 # see: https://doc.modelica.org/Modelica%204.0.0/Resources/helpWSM/Modelica/Modelica.Mechanics.Rotational.Examples.First.html
@@ -132,7 +152,7 @@ end
         sys = structural_simplify(model) #key 7 not found
         prob = ODAEProblem(sys, Pair[], (0, 1.0))
         sol = solve(prob, Rodas4())
-        @test sol.retcode == Success
+        @test SciMLBase.successful_retcode(sol)
     end
     # Plots.plot(sol; vars=[inertia2.w, inertia3.w])
 end
@@ -181,7 +201,7 @@ end
     prob = DAEProblem(sys, D.(states(sys)) .=> 0.0, Pair[], (0, 10.0))
 
     sol = solve(prob, DFBDF())
-    @test sol.retcode == Success
+    @test SciMLBase.successful_retcode(sol)
     @test sol[angle_sensor.phi.u] == sol[inertia.flange_a.phi]
 
     # p1 = Plots.plot(sol; vars=[inertia.flange_a.phi, source.phi], title="Angular Position", labels=["Inertia" "Source"], ylabel="Angle in rad")
@@ -218,7 +238,7 @@ end
 
     prob = ODEProblem(sys, Pair[], (0, 10.0))
     sol = solve(prob, Rodas4())
-    @test sol.retcode == Success
+    @test SciMLBase.successful_retcode(sol)
     @test all(sol[inertia1.w] .== 0)
     @test all(sol[inertia1.w] .== sol[speed_sensor.w.u])
     @test sol[inertia2.w][end]≈0 atol=1e-3 # all energy has dissipated
@@ -227,7 +247,7 @@ end
 
     prob = DAEProblem(sys, D.(states(sys)) .=> 0.0, Pair[], (0, 10.0))
     sol = solve(prob, DFBDF())
-    @test sol.retcode == Success
+    @test SciMLBase.successful_retcode(sol)
     @test all(sol[inertia1.w] .== 0)
     @test all(sol[inertia1.w] .== sol[speed_sensor.w.u])
     @test sol[inertia2.w][end]≈0 atol=1e-3 # all energy has dissipated