add SpringDamper component (#162)

* add SpringDamper component

* format

Author: baggepinnen

docs/src/API/mechanical.md

@@ -35,6 +35,7 @@ Fixed
 Inertia
 Spring
 Damper
+SpringDamper
 IdealGear
 RotationalFriction
 ```

src/Mechanical/Rotational/Rotational.jl

@@ -12,7 +12,7 @@ D = Differential(t)
 export Flange, Support
 include("utils.jl")
 
-export Fixed, Inertia, Spring, Damper, IdealGear, RotationalFriction
+export Fixed, Inertia, Spring, Damper, SpringDamper, IdealGear, RotationalFriction
 include("components.jl")
 
 export Torque, ConstantTorque, Speed

src/Mechanical/Rotational/components.jl

@@ -119,6 +119,42 @@ Linear 1D rotational damper
     extend(ODESystem(eqs, t, [], pars; name = name), partial_comp)
 end
 
+"""
+    SpringDamper(;name, d)
+
+Linear 1D rotational spring and damper
+
+# States:
+
+  - `phi_rel(t)`: [`rad`] Relative rotation angle (= flange_b.phi - flange_a.phi)
+  - `w_rel(t)`: [`rad/s`] Relative angular velocity (= D(phi_rel))
+  - `a_rel(t)`: [`rad/s²`] Relative angular acceleration (= D(w_rel))
+  - `tau(t)`: [`N.m`] Torque between flanges (= flange_b.tau)
+
+# Connectors:
+
+  - `flange_a` [Flange](@ref)
+  - `flange_b` [Flange](@ref)
+
+# Parameters:
+
+  - `d`: [`N.m.s/rad`] Damping constant
+  - `c`: [`N.m/rad`] Spring constant
+"""
+@component function SpringDamper(; name, c, d, phi_rel0 = 0.0)
+    @named partial_comp = PartialCompliantWithRelativeStates()
+    @unpack phi_rel, w_rel, tau = partial_comp
+    @variables tau_c(t) [description = "Spring torque"]
+    @variables tau_d(t) [description = "Damper torque"]
+    @parameters d=d [description = "Damping constant"]
+    @parameters c=c [description = "Spring constant"]
+    @parameters phi_rel0=phi_rel0 [description = "Unstretched spring angle"]
+    eqs = [tau_c ~ c * (phi_rel - phi_rel0)
+           tau_d ~ d * w_rel
+           tau ~ tau_c + tau_d]
+    extend(ODESystem(eqs, t; name = name), partial_comp)
+end
+
 """
     IdealGear(;name, ratio, use_support=false)
 

test/Mechanical/rotational.jl

@@ -24,8 +24,8 @@ D = Differential(t)
     sys = structural_simplify(model)
 
     prob = ODEProblem(sys, Pair[], (0, 10.0))
-    sol = solve(prob, Rodas4())
-    @test SciMLBase.successful_retcode(sol)
+    sol1 = solve(prob, Rodas4())
+    @test SciMLBase.successful_retcode(sol1)
 
     prob = ODAEProblem(sys, Pair[], (0, 10.0))
     sol = solve(prob, Rodas4())
@@ -37,6 +37,21 @@ D = Differential(t)
     @test all(sol[inertia1.w] .== 0)
     @test sol[inertia2.w][end]≈0 atol=1e-3 # all energy has dissipated
 
+    @named springdamper = SpringDamper(; c = 1e4, d = 10)
+    connections = [connect(fixed.flange, inertia1.flange_b)
+                   connect(inertia1.flange_b, springdamper.flange_a)
+                   connect(springdamper.flange_b, inertia2.flange_a)]
+
+    @named model = ODESystem(connections, t,
+                             systems = [fixed, inertia1, inertia2, springdamper])
+    sys = structural_simplify(model)
+
+    prob = ODEProblem(sys, Pair[], (0, 10.0))
+    sol2 = solve(prob, Rodas4())
+    @test SciMLBase.successful_retcode(sol)
+
+    @test sol2(0:1:10, idxs = inertia2.w).u≈sol1(0:1:10, idxs = inertia2.w).u atol=1e-3
+
     # Plots.plot(sol; vars=[inertia1.w, inertia2.w])
 end