adds DC motor with friction

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]
@@ -20,7 +19,9 @@ julia = "1.6"
 
 [extras]
 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 = ["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/Electrical/Analog/ideal_components.jl

@@ -26,7 +26,7 @@ See [OnePort](@ref)
 - `n` Negative pin
 
 # Parameters:
-- `R`: [`Ω`] Resistance
+- `R`: [`Ohm`] Resistance
 """
 function Resistor(;name, R)
     @named oneport = OnePort()
@@ -70,8 +70,7 @@ end
 Creates an ideal capacitor.
 
 # States:
-- `v(t)`: [`V`]
-  The voltage across the capacitor, given by `D(v) ~ p.i / C`
+- `v(t)`: [`V`] The voltage across the capacitor, given by `D(v) ~ p.i / C`
 
 # Connectors:
 - `p` Positive pin
@@ -161,4 +160,78 @@ function Short(; name)
     @unpack v, i = oneport
     eqs = [v ~ 0]
     extend(ODESystem(eqs, t, [], []; name=name), oneport)
+end
+
+"""
+    HeatingResistor(;name, R_ref=1.0, T_ref=300.15, alpha=0)
+
+Temperature dependent electrical resistor
+
+# States
+- See [OnePort](@ref)
+- `R(t)`: [`Ohm`] 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` [Pin](@ref) Positive pin
+- `n` [Pin](@ref) Negative pin
+- `flange` [Flange](@ref) Shaft of EMF shaft
+- `support` [Support](@ref) 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,27 @@ 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 Voltage, Current
 include("Analog/sources.jl")
+
 # include("Digital/components.jl")
 # include("Digital/gates.jl")
 # include("Digital/tables.jl")
@@ -26,23 +35,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/components.jl

@@ -39,6 +39,7 @@ end
 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
@@ -75,6 +76,7 @@ Linear 1D rotational spring
 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  
@@ -104,6 +106,7 @@ Linear 1D rotational damper
 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)
@@ -132,6 +135,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 = [ 
@@ -143,7 +147,6 @@ function IdealGear(;name, ratio, use_support=false)
     extend(ODESystem(eqs, t, sts, [ratio]; name=name), partial_element)
 end
 
-
 """
     RotationalFriction(;name, f, tau_c, w_brk, tau_brk) 
 

src/Mechanical/Rotational/sources.jl

@@ -1,5 +1,5 @@
 """
-    Torque(;name) 
+    Torque(; name, use_support=false) 
 
 Input signal acting as external torque on a flange
 
@@ -22,7 +22,6 @@ function Torque(;name, use_support=false)
 end
 
 
-
 """
     Speed(; name, use_support=false, exact=false, f_crit=50) 
 

src/Mechanical/Rotational/utils.jl

@@ -121,7 +121,7 @@ Partial model for a component with one rotational 1-dim. shaft flange and a supp
 function PartialElementaryOneFlangeAndSupport2(;name, use_support=false)
     @named flange = Flange()
     sys = [flange]
-    @variables phi_support(t)
+    @variables phi_support(t)=0.0
     if use_support
         @named support = Support() 
         eqs = [

src/ModelingToolkitStandardLibrary.jl

@@ -2,8 +2,8 @@ module ModelingToolkitStandardLibrary
 
 include("Blocks/Blocks.jl")
 include("Mechanical/Mechanical.jl")
+include("Thermal/Thermal.jl")
 include("Electrical/Electrical.jl")
 include("Magnetic/Magnetic.jl")
-include("Thermal/Thermal.jl")
 
 end

test/Blocks/sources.jl

@@ -97,7 +97,6 @@ end
     sys = structural_simplify(iosys)
     prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 10.0))
     sol = solve(prob, Rodas4())
-
     @test sol.retcode == :Success
     @test sol[src.output.u] ≈ cosine.(sol.t, frequency, amplitude, phase, offset, start_time) atol=1e-3
 
@@ -161,7 +160,6 @@ end
     sys = structural_simplify(iosys)
     prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 10.0))
     sol = solve(prob, Rodas4())
-
     @test sol.retcode == :Success
     @test sol[src.output.u] ≈ ramp.(sol.t, offset, height, duration, start_time) atol=1e-3
 
@@ -350,7 +348,6 @@ end
     sys = structural_simplify(iosys)
     prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 10.0))
     sol = solve(prob, Rodas4())
-
     @test sol.retcode == :Success
     @test sol[src.output.u] ≈ exp_sine.(sol.t, amplitude, frequency, damping, phase, start_time) atol=1e-3
 

test/Electrical/analog.jl

@@ -149,10 +149,10 @@ end
         @named model = ODESystem(connections, t; systems=[resistor, capacitor, source, ground, voltage])
         sys = structural_simplify(model)
         prob = ODAEProblem(sys, [capacitor.v => 10.0], (0.0, 10.0))
-        sol = solve(prob, Rodas5())
-        @test sol.retcode == :Success
         sol = solve(prob, Tsit5())
         @test sol.retcode == :Success
+        sol = solve(prob, Rodas4())
+        @test sol.retcode == :Success
 
         # Plots.plot(sol; vars=[voltage.v, capacitor.v])
     end
@@ -272,7 +272,6 @@ _damped_sine_wave(x, f, A, st, ϕ, d) = exp((st-x)*d)*A*sin(2*π*f*(x-st) + ϕ)
 
         prob = ODAEProblem(vsys, u0, (0, 10.0))
         sol = solve(prob, dt=0.1, Tsit5())
-
         @test sol.retcode == :Success
         @test sol[voltage.V.u] ≈ waveforms(i, sol.t) atol=1e-1
         @test sol[voltage.p.v] ≈ sol[voltage.V.u]
@@ -326,12 +325,11 @@ end
 
         prob = ODAEProblem(isys, u0, (0, 10.0))
         sol = solve(prob, dt=0.1, Tsit5())
-
         @test sol.retcode == :Success
         @test sol[current.I.u] ≈ waveforms(i, sol.t) atol=1e-1
         @test sol[current.I.u] ≈ sol[current.p.i] atol=1e-1
         # For visual inspection
         # plt = plot(sol)
         # savefig(plt, "test_current_$(source.name)")
     end
-end
+end