adds DC motor with friction

adds speed sensor test
adds DC motor test in runtests.jl
adds test for speed sensor
adds rel. speed sensor
adds tests
adds comment
adds parameter checks
adds ViscousFriction
adds checks for return codes
Update analog.jl
Update rotational.jl
fix test
adds missing default value
add missing link to tutorial

Author: ValentinKaisermayer

Project.toml

@@ -7,7 +7,6 @@ version = "1.4.0"
 IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
-OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [compat]
@@ -19,8 +18,11 @@ Symbolics = "0.1, 1, 2, 3, 4"
 julia = "1.6"
 
 [extras]
+DASSL = "e993076c-0cfd-5d6b-a1ac-36489fdf7917"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["SafeTestsets", "Test"]
+test = ["DASSL", "SafeTestsets", "Sundials", "OrdinaryDiffEq", "Test"]

src/Blocks/Blocks.jl

@@ -2,7 +2,7 @@
 The module `Blocks` contains common input-output components, referred to as blocks.
 """
 module Blocks
-using ModelingToolkit, Symbolics, OrdinaryDiffEq
+using ModelingToolkit, Symbolics
 using IfElse: ifelse
 
 @parameters t

src/Blocks/continuous.jl

@@ -6,6 +6,10 @@ Outputs `y = ∫k*u dt`, corresponding to the transfer function `1/s`.
 # Connectors:
 - `input`
 - `output`
+
+# Parameters:
+- `k`: Gain of integrator
+- `x_start`: Initial value of integrator
 """
 function Integrator(;name, k=1, x_start=0.0)
     @named siso = SISO()

src/Blocks/utils.jl

@@ -1,4 +1,4 @@
-@connector function RealInput(;name, nin=1, u_start=nin > 1 ? 0.0 : zeros(nin))
+@connector function RealInput(;name, nin=1, u_start=nin > 1 ? zeros(nin) : 0.0)
     if nin == 1
         @variables u(t) = u_start
     else
@@ -20,7 +20,7 @@ Connector with one input signal of type Real.
 - `u`: Value of of the connector; if nin=1 this is a scalar
 """ RealInput
 
-@connector function RealOutput(;name, nout=1, u_start=nout > 1 ? 0.0 : zeros(nout))
+@connector function RealOutput(;name, nout=1, u_start=nout > 1 ? zeros(nout) : 0.0)
     if nout == 1
         @variables u(t) = u_start
     else

src/Electrical/Analog/ideal_components.jl

@@ -159,3 +159,78 @@ function IdealOpAmp(;name)
     ]
     ODESystem(eqs, t, sts, [], systems=[p1, p2, n1, n2], name=name)
 end
+
+
+"""
+    HeatingResistor(;name, R_ref=1.0, T_ref=300.15, alpha=0)
+
+Temperature dependent electrical resistor
+
+# States
+- See [OnePort](@ref)
+- `R(t)`: [`Ω`] Temperature dependent resistance `R ~ R_ref*(1 + alpha*(heat_port.T(t) - T_ref))`
+
+# Connectors
+- `p` Positive pin
+- `n` Negative pin
+
+# Parameters: 
+- `R_ref`: [`Ω`] Reference resistance 
+- `T_ref`: [K] Reference temperature
+"""
+function HeatingResistor(;name, R_ref=1.0, T_ref=300.15, alpha=0)
+    @named oneport = OnePort()
+    @unpack v, i = oneport
+    @named heat_port = HeatPort()
+    pars = @parameters begin
+      R_ref=R_ref
+      T_ref=T_ref
+      alpha=alpha
+    end
+    @variables R
+    eqs = [
+      R ~ R_ref*(1 + alpha*(heat_port.T - T_ref))
+      heat_port.Q_flow ~ -v * i # -LossPower
+      v ~ i * R
+    ]
+    extend(ODESystem(eqs, t, [R], pars; name=name, systems=[heat_port]), oneport)
+end
+
+"""
+    EMF(;name, k)
+
+Electromotoric force (electric/mechanic transformer)
+
+# States
+- `v(t)`: [`V`] The voltage across component `p.v - n.v`
+- `i(t)`: [`A`] The current passing through positive pin
+- `phi`: [`rad`] Rotation angle (=flange.phi - support.phi)
+- `w`: [`rad/s`] Angular velocity (= der(phi))
+
+# Connectors
+- `p` Positive pin
+- `n` Negative pin
+- `flange` Shaft of EMF shaft
+- `support` Support/housing of emf shaft
+
+# Parameters: 
+- `k`: [`N⋅m/A`] Transformation coefficient 
+"""
+function EMF(;name, k)
+    @named p = Pin()
+    @named n = Pin()
+    @named flange = Flange()
+    @named support = Support()
+    @parameters k=k
+    @variables v(t)=0.0 i(t)=0.0 phi(t)=0.0 w(t)=0.0
+    eqs = [
+        v ~ p.v - n.v
+        0 ~ p.i + n.i
+        i ~ p.i
+        phi ~ flange.phi - support.phi
+        D(phi) ~ w
+        k*w ~ v
+        flange.tau ~ -k*i
+    ]
+    ODESystem(eqs, t, [v, i, phi, w], [k]; name=name, systems=[p, n, flange, support])
+end

src/Electrical/Electrical.jl

@@ -4,18 +4,28 @@ This library contains electrical components to build up analog circuits.
 """
 module Electrical
 
-using ModelingToolkit, Symbolics, IfElse, OrdinaryDiffEq
+using ModelingToolkit, Symbolics, IfElse
 using OffsetArrays
+using ..Thermal: HeatPort
+using ..Mechanical.Rotational: Flange, Support
+using ..Blocks: RealInput, RealOutput
 
 @parameters t
 D = Differential(t)
 
-using ..Blocks: RealInput, RealOutput
-
+export Pin, OnePort
 include("utils.jl")
+
+export Capacitor, Ground, Inductor, Resistor, Short, IdealOpAmp, EMF
 include("Analog/ideal_components.jl")
+
+export CurrentSensor, PotentialSensor, VoltageSensor, PowerSensor, MultiSensor
 include("Analog/sensors.jl")
+
+export SignalVoltage, ConstantVoltage, SineVoltage, StepVoltage, RampVoltage, SquareVoltage, TriangularVoltage, CosineVoltage, ExpSineVoltage
+export SignalCurrent, ConstantCurrent, SineCurrent, StepCurrent, RampCurrent, SquareCurrent, TriangularCurrent, CosineCurrent, ExpSineCurrent
 include("Analog/sources.jl")
+
 # include("Digital/components.jl")
 # include("Digital/gates.jl")
 # include("Digital/tables.jl")
@@ -26,23 +36,4 @@ include("Analog/sources.jl")
 # - machines
 # - multi-phase
 
-export #Interface
-       Pin,
-       # Analog Components
-       Capacitor, Ground, Inductor, Resistor, Conductor,
-       Short, IdealOpAmp,
-       # Analog Sensors
-       CurrentSensor, PotentialSensor, VoltageSensor,
-       PowerSensor, MultiSensor,
-       # Analog Sources
-       Voltage, Current
-
-
-       # # Digital Gates
-       # And, Or, Not, Xor, Nand, Nor, Xnor,
-       # # Digital components
-       # HalfAdder, FullAdder, MUX, DEMUX, Encoder, Decoder,
-       # # Digital Sources
-       # DigitalPin, Pulse, PulseDiff
-
 end

src/Mechanical/Rotational/Rotational.jl

@@ -3,20 +3,23 @@ Library to model 1-dimensional, rotational mechanical systems
 """
 module Rotational
 
-using ModelingToolkit, Symbolics, IfElse, OrdinaryDiffEq
+using ModelingToolkit, Symbolics, IfElse
 using OffsetArrays
 using ...Blocks: RealInput, RealOutput
 
 @parameters t
 D = Differential(t)
 
-export Flange
+export Flange, Support
 include("utils.jl")
 
-export Fixed, Inertia, Spring, Damper, IdealGear
+export Fixed, Inertia, Spring, Damper, IdealGear, RotationalFriction, ViscousFriction
 include("components.jl")
 
-export Torque
+export Torque, Speed
 include("sources.jl")
 
+export AngleSensor, SpeedSensor, TorqueSensor, RelSpeedSensor
+include("sensors.jl")
+
 end

src/Mechanical/Rotational/components.jl

@@ -3,35 +3,43 @@
 
 Flange fixed in housing at a given angle.
 
+# Connectors:
+- `flange` flange
+
 # Parameters:
 - `phi0`: [rad] Fixed offset angle of housing
 """
 function Fixed(;name, phi0=0.0)
     @named flange = Flange()
-    @parameters phi0=phi0 
+    # @parameters phi0=phi0 
     eqs = [flange.phi ~ phi0]
-    return compose(ODESystem(eqs, t, [], [phi0]; name=name), flange)
+    return compose(ODESystem(eqs, t, [], []; name=name), flange)
 end
 
 """
     Inertia(;name, J, phi_start=0.0, w_start=0.0, a_start=0.0)
 
 1D-rotational component with inertia.
 
+# States: 
+- `phi`: [rad] Absolute rotation angle of component 
+- `w`: [rad/s] Absolute angular velocity of component (= der(phi)) 
+- `a`: [rad/s²] Absolute angular acceleration of component (= der(w)) 
+
+# Connectors:
+- `flange_a` Left flange
+- `flange_b` Right flange
+
 # Parameters:
 - `J`: [kg·m²] Moment of inertia 
 - `phi_start`: [rad] Initial value of absolute rotation angle of component 
 - `w_start`: [rad/s] Initial value of absolute angular velocity of component
 - `a_start`: [rad/s²] Initial value of absolute angular acceleration of component
-
-# States: 
-- `phi`: [rad] Absolute rotation angle of component 
-- `w`: [rad/s] Absolute angular velocity of component (= der(phi)) 
-- `a`: [rad/s²] Absolute angular acceleration of component (= der(w)) 
 """
 function Inertia(;name, J, phi_start=0.0, w_start=0.0, a_start=0.0)
     @named flange_a = Flange()
     @named flange_b = Flange()
+    J > 0 || throw(ArgumentError("Expected `J` to be positive"))
     @parameters J=J
     sts = @variables begin
         phi(t)=phi_start
@@ -53,13 +61,20 @@ end
 
 Linear 1D rotational spring
 
+# States:
+- See [PartialCompliant](@ref)
+
+# Connectors:
+- See [PartialCompliant](@ref)
+
 # Parameters:
 - `c`: [N.m/rad] Spring constant
 - `phi_rel0`: Unstretched spring angle
 """
 function Spring(;name, c, phi_rel0=0.0)
     @named partial_comp = PartialCompliant()
     @unpack phi_rel, tau = partial_comp
+    c > 0 || throw(ArgumentError("Expected `c` to be positive"))
     pars = @parameters begin
         c=c 
         phi_rel0=phi_rel0  
@@ -73,12 +88,19 @@ end
 
 Linear 1D rotational damper
 
+# States:
+- See [PartialCompliantWithRelativeStates](@ref)
+
+# Connectors:
+- See [PartialCompliantWithRelativeStates](@ref)
+
 # Parameters:
 - `d`: [N.m.s/rad] Damping constant
 """
 function Damper(;name, d) 
     @named partial_comp = PartialCompliantWithRelativeStates()
     @unpack w_rel, tau = partial_comp
+    d > 0 || throw(ArgumentError("Expected `d` to be positive"))
     pars = @parameters d=d 
     eqs = [tau ~ d*w_rel]
     extend(ODESystem(eqs, t, [], pars; name=name), partial_comp)
@@ -98,6 +120,7 @@ This element characterizes any type of gear box which is fixed in the ground and
 function IdealGear(;name, ratio, use_support=false) 
     @named partial_element = PartialElementaryTwoFlangesAndSupport2(use_support=use_support)
     @unpack phi_support, flange_a, flange_b = partial_element
+    ratio > 0 || throw(ArgumentError("Expected `ratio` to be positive"))
     @parameters ratio=ratio
     sts = @variables phi_a(t)=0.0 phi_b(t)=0.0
     eqs = [ 
@@ -107,4 +130,70 @@ function IdealGear(;name, ratio, use_support=false)
         0 ~ ratio*flange_a.tau + flange_b.tau
     ]
     extend(ODESystem(eqs, t, sts, [ratio]; name=name), partial_element)
+end
+
+"""
+    RotationalFriction(;name, f, tau_c, w_brk, tau_brk) 
+
+Models rotational friction with Stribeck effect, Coulomb friction and viscous friction between the two flanges.
+The friction torque is a function of the relative angular velocity between flange_a and flange_b.
+
+Friction model: "Armstrong, B. and C.C. de Wit, Friction Modeling and Compensation, The Control Handbook, CRC Press, 1995."
+
+# States:
+- see [`PartialCompliantWithRelativeStates`](@ref)
+
+# Connectors:
+- See [`PartialCompliantWithRelativeStates`](@ref)
+
+# Parameters:
+- `f`: [`N⋅m/(rad/s)`] Viscous friction coefficient 
+- `tau_c`: [`N⋅m`] Coulomb friction torque
+- `w_brk`: [`rad/s`] Breakaway friction velocity 
+- `tau_brk`: [`N⋅m`] Breakaway friction torque
+"""
+function RotationalFriction(;name, f, tau_c, w_brk, tau_brk) 
+    @named partial_comp = PartialCompliantWithRelativeStates()
+    @unpack w_rel, tau = partial_comp
+    w_brk > 0 || throw(ArgumentError("Expected `w_brk` to be positive"))
+    tau_brk > 0 || throw(ArgumentError("Expected `tau_brk` to be positive"))
+    tau_c > 0 || throw(ArgumentError("Expected `tau_c` to be positive"))
+    tau_c <= tau_brk || throw(ArgumentError("`tau_brk` to be greater than `tau_c`"))
+    f >= 0 || throw(ArgumentError("Expected `f` to be not negative"))
+    pars = @parameters f=f tau_c=tau_c w_brk=w_brk tau_brk=tau_brk
+
+    str_scale = sqrt(2*exp(1)) * (tau_brk - tau_c)
+    w_st = w_brk * sqrt(2)
+    w_coul = w_brk / 10                     
+
+    eqs = [ 
+        tau ~ str_scale * (exp(-(w_rel/w_st)^2) * w_rel / w_st) + tau_c * tanh(w_rel / w_coul) + f * w_rel # Stribeck friction + Coulomb friction + Viscous friction
+    ]
+    extend(ODESystem(eqs, t, [], pars; name=name), partial_comp)
+end
+
+"""
+    ViscousFriction(;name, f) 
+
+Models rotational viscous friction between the two flanges.
+The friction torque is a function of the relative angular velocity between flange_a and flange_b.
+# States:
+- see [`PartialCompliantWithRelativeStates`](@ref)
+
+# Connectors:
+- See [`PartialCompliantWithRelativeStates`](@ref)
+
+# Parameters:
+- `f`: [`N⋅m/(rad/s)`] Viscous friction coefficient 
+"""
+function ViscousFriction(;name, f) 
+    @named partial_comp = PartialCompliantWithRelativeStates()
+    @unpack w_rel, tau = partial_comp
+    f >= 0 || throw(ArgumentError("Expected `f` to be not negative"))
+    pars = @parameters f=f
+
+    eqs = [ 
+        tau ~ f * w_rel # Viscous friction
+    ]
+    extend(ODESystem(eqs, t, [], pars; name=name), partial_comp)
 end