Simplify temperature dependency for analog components

Author: langestefan

src/Electrical/Analog/ideal_components.jl

@@ -18,30 +18,57 @@ node.
 end
 
 """
-    Resistor(; name, R)
+    Resistor(; name, R_ref = 1.0, T_ref = 300.15, alpha = 0, T_dep = false)
 
-Creates an ideal Resistor following Ohm's Law.
+Generic resistor with optional temperature dependency.
 
 # States:
 
-See [OnePort](@ref)
+  - See [OnePort](@ref)
+  - `R(t)`: [`Ω`] Resistance (temperature dependent if `T_dep = true`)
 
 # Connectors:
 
   - `p` Positive pin
   - `n` Negative pin
+  - `heat_port` [HeatPort](@ref) (only if `T_dep = true`) Heat port to model the temperature dependency
 
 # Parameters:
 
-  - `R`: [`Ohm`] Resistance
+  - `R`: [`Ω`] Reference resistance
+  - `T_ref`: [K] Reference temperature
+  - `alpha`: [K⁻¹] Temperature coefficient of resistance
+  - `T_dep`: [bool] Temperature dependency
 """
 @mtkmodel Resistor begin
     @extend v, i = oneport = OnePort()
+
+    @structural_parameters begin
+        T_dep = false
+    end
+
     @parameters begin
-        R, [description = "Resistance"]
+        R = 1.0, [description = "Reference resistance", unit = "Ω"]
+        T_ref = 300.15, [description = "Reference temperature", unit = "K"]
+        alpha = 0.0, [description = "Temperature coefficient of resistance", unit = "K⁻¹"]
     end
-    @equations begin
-        v ~ i * R
+
+    if T_dep
+        @components begin
+            heat_port = HeatPort()
+        end
+        @variables begin
+            R_T(t), [description = "Temperature-dependent resistance", unit = "Ω"]
+        end
+        @equations begin
+            R_T ~ R * (1 + alpha * (heat_port.T - T_ref))  # Temperature-dependent resistance
+            heat_port.Q_flow ~ -v * i  # -LossPower
+            v ~ i * R_T  # Ohm's Law
+        end
+    else
+        @equations begin
+            v ~ i * R  # Ohm's Law for constant resistance
+        end
     end
 end
 
@@ -264,9 +291,9 @@ Electromotoric force (electric/mechanic transformer)
 end
 
 """
-        Diode(; name, Is = 1e-6, n = 1, T = 300.15)
+    Diode(; name, Is = 1e-6, n = 1, T_ref = 300.15, T_dep = false)
 
-Ideal diode based on the Shockley diode equation.
+Generic diode with optional temperature dependency.
 
 # States
 
@@ -276,70 +303,49 @@ Ideal diode based on the Shockley diode equation.
 
     - `p` Positive pin
     - `n` Negative pin
+    - `port` [HeatPort](@ref) (only if `T_dep = true`) Heat port to model the temperature dependency
 
-# Parameters
-
+# Parameters:
+    
     - `Is`: [`A`] Saturation current
     - `n`: Ideality factor
-    - `T`: [K] Ambient temperature
+    - `T_ref`: [K] Reference temperature
+    - `T_dep`: [bool] Temperature dependency
 """
 @mtkmodel Diode begin
     @constants begin
         k = 1.380649e-23 # Boltzmann constant (J/K)
         q = 1.602176634e-19 # Elementary charge (C)
     end
     @extend v, i = oneport = OnePort(; v = 0.0)
-    @parameters begin
-        Is = 1e-6, [description = "Saturation current (A)"]
-        n = 1, [description = "Ideality factor"]
-        T = 300.15, [description = "Ambient temperature"]
-    end
-    @equations begin
-        i ~ Is * (exp(v * q / (n * k * T)) - 1)
-    end
-end
-
-"""
-    HeatingDiode(; name, Is = 1e-6, n = 1)
 
-Temperature dependent diode based on the Shockley diode equation.
-
-# States
-
-    - See [OnePort](@ref)
-
-# Connectors
-
-    - `p` Positive pin
-    - `n` Negative pin
-    - `port` [HeatPort](@ref) Heat port to model the temperature dependency
-
-# Parameters:
-
-    - `Is`: [`A`] Saturation current
-    - `n`: Ideality factor
-"""
-@mtkmodel HeatingDiode begin
-    begin
-        k = 1.380649e-23 # Boltzmann constant (J/K)
-        q = 1.602176634e-19 # Elementary charge (C)
+    @structural_parameters begin
+        T_dep = false
     end
 
-    @extend v, i = oneport = OnePort(; v = 0.0)
-    @components begin
-        port = HeatPort()
-    end
     @parameters begin
         Is = 1e-6, [description = "Saturation current (A)"]
         n = 1, [description = "Ideality factor"]
+        T_ref = 300.15, [description = "Reference temperature (K)"]
+        Vt_const = k * T_ref / q, [description = "Constant thermal voltage"]
     end
-    @variables begin
-        Vt(t), [description = "Thermal voltage"]
-    end
-    @equations begin
-        Vt ~ k * port.T / q  # Thermal voltage equation
-        i ~ Is * (exp(v / (n * Vt)) - 1)  # Shockley diode equation
-        port.Q_flow ~ -v * i  # -LossPower
+
+    if T_dep
+        @components begin
+            port = HeatPort()
+        end
+        @variables begin
+            Vt(t), [description = "Thermal voltage"]
+        end
+        @equations begin
+            Vt ~ k * port.T / q  # Thermal voltage equation
+            i ~ Is * (exp(v / (n * Vt)) - 1)  # Shockley diode equation with temperature dependence
+            port.Q_flow ~ -v * i  # -LossPower
+        end
+    else
+        @equations begin
+            i ~ Is * (exp(v / (n * Vt_const)) - 1)  # Shockley diode equation
+        end
     end
 end
 
@@ -366,20 +372,20 @@ R = R_const + pos * R_ref * (1 + alpha * (port.T - T_ref))
     - `R(t)`: Resistance
 
 # Connectors
-
-        - `p` Positive pin
-        - `n` Negative pin
-        - `position` RealInput to set the position of the wiper
-        - `port` [HeatPort](@ref) Heat port to model the temperature dependency
+        
+    - `p` Positive pin
+    - `n` Negative pin
+    - `position` RealInput to set the position of the wiper
+    - `port` [HeatPort](@ref) Heat port to model the temperature dependency
 
 # Parameters
-
-        - `R_ref`: [`Ω`] Resistance at temperature T_ref when fully closed (pos=1.0)
-        - `T_ref`: [K] Reference temperature
-        - `R_const`: [`Ω`] Constant resistance between p and n
-        - `T_dep`: Temperature dependency
-        - `alpha`: [K⁻¹] Temperature coefficient of resistance
-        - `enforce_bounds`: Enforce bounds for the position of the wiper (0-1)
+    
+    - `R_ref`: [`Ω`] Resistance at temperature T_ref when fully closed (pos=1.0)
+    - `T_ref`: [K] Reference temperature
+    - `R_const`: [`Ω`] Constant resistance between p and n
+    - `T_dep`: [bool] Temperature dependency
+    - `alpha`: [K⁻¹] Temperature coefficient of resistance
+    - `enforce_bounds`: Enforce bounds for the position of the wiper (0-1)
 """
 @mtkmodel VariableResistor begin
     @extend v, i = oneport = OnePort()

src/Electrical/Electrical.jl

@@ -15,7 +15,7 @@ include("utils.jl")
 
 export Capacitor,
        Ground, Inductor, Resistor, Conductor, Short, IdealOpAmp, EMF,
-       HeatingResistor, Diode, HeatingDiode, VariableResistor
+       Diode, VariableResistor
 include("Analog/ideal_components.jl")
 
 export CurrentSensor, PotentialSensor, VoltageSensor, PowerSensor, MultiSensor

test/Electrical/analog.jl

@@ -416,45 +416,45 @@ end
     @test capacitor_voltage[end].≈8.26 rtol=3e-1
 
     # For visual inspection
-    # plt = plot(sol; vars = [diode.i, resistor.i, capacitor.v],
+    # plt = plot(sol; idxs = [diode.i, resistor.i, capacitor.v],
     #     size = (800, 600), dpi = 300,
     #     labels = ["Diode Current" "Resistor Current" "Capacitor Voltage"],
     #     title = "Diode Test")
     # savefig(plt, "diode_test")
 end
 
-@testset "HeatingDiode component test" begin
-    @mtkmodel HeatingDiodeTest begin
-        @parameters begin
-            R = 1.0
-            C = 1.0
-            V = 10.0
-            T = 300.0 # Ambient temperature in Kelvin
-            n = 2.0
-            Is = 1e-6
-            f = 1.0
-        end
-        @components begin
-            resistor = Resistor(R = R)
-            capacitor = Capacitor(C = C, v = 0.0)
-            source = Voltage()
-            heating_diode = HeatingDiode(n = n, Is = Is)
-            ac = Sine(frequency = f, amplitude = V)
-            ground = Ground()
-            temp = FixedTemperature(T = T)
-        end
-        @equations begin
-            connect(ac.output, source.V)
-            connect(source.p, heating_diode.p)
-            connect(heating_diode.n, resistor.p)
-            connect(resistor.n, capacitor.p)
-            connect(capacitor.n, ground.g)
-            connect(source.n, ground.g)
-            connect(temp.port, heating_diode.port)
-        end
-    end
-
-    @mtkbuild sys = HeatingDiodeTest()
+@testset "Diode with temperature dependency component test" begin
+    # Parameter values
+    R = 1.0
+    C = 1.0
+    V = 10.0
+    T = 300.0 # Ambient temperature in Kelvin
+    n = 2.0
+    Is = 1e-6
+    f = 1.0
+
+    # Components
+    @named resistor = Resistor(R = R)
+    @named capacitor = Capacitor(C = C, v = 0.0)
+    @named source = Voltage()
+    @named heating_diode = Diode(n = n, Is = Is, T_dep = true)
+    @named ac = Sine(frequency = f, amplitude = V)
+    @named ground = Ground()
+    @named temp = FixedTemperature(T = T)
+
+    # Connections
+    connections = [connect(ac.output, source.V),
+        connect(source.p, heating_diode.p),
+        connect(heating_diode.n, resistor.p),
+        connect(resistor.n, capacitor.p),
+        connect(capacitor.n, ground.g),
+        connect(source.n, ground.g),
+        connect(temp.port, heating_diode.port)]
+
+    # Model
+    @named model = ODESystem(connections, t;
+        systems = [resistor, capacitor, source, heating_diode, ac, ground, temp])
+    sys = structural_simplify(model)
     prob = ODEProblem(sys, Pair[], (0.0, 10.0))
     sol = solve(prob, Rodas4())
 
@@ -473,10 +473,10 @@ end
     @test capacitor_voltage[end]≈7.75 rtol=3e-1 # Final capacitor voltage close to input voltage
 
     # For visual inspection
-    # plt = plot(sol; vars = [heating_diode.i, resistor.i, capacitor.v],
+    # plt = plot(sol; idxs = [heating_diode.i, resistor.i, capacitor.v],
     #     size = (800, 600), dpi = 300,
-    #     labels = ["HeatingDiode Current" "Resistor Current" "Capacitor Voltage"],
-    #     title = "HeatingDiode Test")
+    #     labels = ["Diode Current" "Resistor Current" "Capacitor Voltage"],
+    #     title = "Diode Test")
     # savefig(plt, "heating_diode_test")
 
     # Remake model with higher amb. temperature, final capacitor voltage should be lower

test/multi_domain.jl

@@ -181,7 +181,8 @@ end
             ground = Ground()
             source = Voltage()
             voltage_sine = Blocks.Sine(amplitude = 220, frequency = 1)
-            heating_resistor = HeatingResistor(R_ref = 100, alpha = 1e-3, T_ref = 293.15)
+            heating_resistor = Resistor(R = 100, alpha = 1e-3,
+                T_ref = 293.15, T_dep = true)
             thermal_conductor = ThermalConductor(G = 50)
             env = FixedTemperature(T = 273.15 + 20)
         end