adds DC motor with friction

Author: ValentinKaisermayer

Project.toml

@@ -13,7 +13,7 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 IfElse = "0.1"
 ModelingToolkit = "8"
 OffsetArrays = "1"
-Symbolics = "0.1, 1, 2, 3, 4"
+Symbolics = "4.9"
 julia = "1.6"
 
 [extras]

src/Blocks/continuous.jl

@@ -111,7 +111,7 @@ Critical damping corresponds to `d=1`, which yields the fastest step response wi
 
 # Parameters:
 - `k`: Gain
-- `w`: Angular frequency
+- `w`: [`rad/s`] Angular frequency
 - `d`: Damping
 - `x_start`: Initial value of state (output)
 - `xd_start`: Initial value of derivative of state (output)
@@ -438,7 +438,7 @@ function StateSpace(; A, B, C, D = nothing, x_start = zeros(size(A, 1)), name)
     end
     @named input = RealInput(nin = nu)
     @named output = RealOutput(nout = ny)
-    @variables x[1:nx](t) = x_start
+    @variables x(t)[1:nx] = x_start
     # pars = @parameters A=A B=B C=C D=D # This is buggy
     eqs = [ # FIXME: if array equations work
         [Differential(t)(x[i]) ~ sum(A[i, k] * x[k] for k in 1:nx) +

src/Blocks/utils.jl

@@ -2,7 +2,7 @@
     if nin == 1
         @variables u(t)=u_start [input = true]
     else
-        @variables u[1:nin](t)=u_start [input = true]
+        @variables u(t)[1:nin]=u_start [input = true]
         u = collect(u)
     end
     ODESystem(Equation[], t, [u...], []; name = name)
@@ -24,7 +24,7 @@ Connector with one input signal of type Real.
     if nout == 1
         @variables u(t)=u_start [output = true]
     else
-        @variables u[1:nout](t)=u_start [output = true]
+        @variables u(t)[1:nout]=u_start [output = true]
         u = collect(u)
     end
     ODESystem(Equation[], t, [u...], []; name = name)
@@ -45,7 +45,7 @@ Connector with one output signal of type Real.
 """
     SISO(;name, u_start=0.0, y_start=0.0)
 
-Single Input Single Output continuous control block.
+Single input single output (SISO) continuous system block.
 
 # Parameters:
 - `u_start`: Initial value for the input
@@ -66,7 +66,7 @@ end
 """
     MIMO(;name, nin=1, nout=1, u_start=zeros(nin), y_start=zeros(nout))
 
-Base class for a multiple Input multiple Output continuous control block.
+Base class for a multiple input multiple output (MIMO) continuous system block.
 
 # Parameters:
 - `nin`: Input dimension
@@ -78,8 +78,8 @@ function MIMO(; name, nin = 1, nout = 1, u_start = zeros(nin), y_start = zeros(n
     @named input = RealInput(nin = nin, u_start = u_start)
     @named output = RealOutput(nout = nout, u_start = y_start)
     @variables begin
-        u[1:nin](t) = u_start
-        y[1:nout](t) = y_start
+        u(t)[1:nin] = u_start
+        y(t)[1:nout] = y_start
     end
     eqs = [
         [u[i] ~ input.u[i] for i in 1:nin]...,

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
@@ -158,5 +157,79 @@ function Short(; name)
     @named oneport = OnePort()
     @unpack v, i = oneport
     eqs = [v ~ 0]
-    extend(ODESystem(eqs, t, [], []; name = name), oneport)
+    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

@@ -6,16 +6,25 @@ module Electrical
 
 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,22 +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,7 +39,8 @@ 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()
-    @parameters J = J
+    J > 0 || throw(ArgumentError("Expected `J` to be positive"))
+    @parameters J=J
     sts = @variables begin
         phi(t) = phi_start
         w(t) = w_start
@@ -73,6 +74,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
@@ -102,9 +104,10 @@ Linear 1D rotational damper
 function Damper(; name, d)
     @named partial_comp = PartialCompliantWithRelativeStates()
     @unpack w_rel, tau = partial_comp
-    pars = @parameters d = d
-    eqs = [tau ~ d * w_rel]
-    extend(ODESystem(eqs, t, [], pars; name = name), 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)
 end
 
 """
@@ -130,7 +133,8 @@ 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
-    @parameters ratio = ratio
+    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 = [phi_a ~ flange_a.phi - phi_support
            phi_b ~ flange_b.phi - phi_support

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
 

src/Mechanical/Rotational/utils.jl

@@ -118,7 +118,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 = [support.phi ~ phi_support

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

@@ -102,7 +102,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
 
@@ -169,7 +168,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
 
@@ -365,7 +363,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

@@ -157,10 +157,10 @@ end
                                  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