Merge pull request #343 from langestefan/diode

Add diode component

Author: ChrisRackauckas

docs/src/API/electrical.md

@@ -31,6 +31,8 @@ Conductor
 Capacitor
 Inductor
 IdealOpAmp
+Diode
+HeatingDiode
 ```
 
 ## Analog Sensors

src/Electrical/Analog/ideal_components.jl

@@ -218,6 +218,7 @@ Temperature dependent electrical resistor
         v ~ i * R
     end
 end
+
 """
     EMF(; name, k)
 
@@ -261,3 +262,84 @@ Electromotoric force (electric/mechanic transformer)
         flange.tau ~ -k * i
     end
 end
+
+"""
+        Diode(; name, Is = 1e-6, n = 1, T = 300.15)
+
+Ideal diode based on the Shockley diode equation.
+
+# States
+
+    - See [OnePort](@ref)
+
+# Connectors
+    
+    - `p` Positive pin
+    - `n` Negative pin
+
+# Parameters
+     
+    - `Is`: [`A`] Saturation current
+    - `n`: Ideality factor
+    - `T`: [K] Ambient temperature
+"""
+@mtkmodel Diode begin
+    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)
+    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"]
+    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
+    end
+end

src/Electrical/Electrical.jl

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

test/Electrical/analog.jl

@@ -4,6 +4,7 @@ using ModelingToolkitStandardLibrary.Blocks: Step,
                                              Constant, Sine, Cosine, ExpSine, Ramp,
                                              Square, Triangular
 using ModelingToolkitStandardLibrary.Blocks: square, triangular
+using ModelingToolkitStandardLibrary.Thermal: FixedTemperature
 using OrdinaryDiffEq: ReturnCode.Success
 
 # using Plots
@@ -372,3 +373,116 @@ end
         # savefig(plt, "test_current_$(source.name)")
     end
 end
+
+@testset "Diode component test" begin
+    # Parameter values 
+    R = 1.0
+    C = 1.0
+    V = 10.0
+    n = 1.0
+    Is = 1e-3
+    f = 1.0
+
+    # Components
+    @named resistor = Resistor(R = R)
+    @named capacitor = Capacitor(C = C, v = 0.0)
+    @named source = Voltage()
+    @named diode = Diode(n = n, Is = Is)
+    @named ac = Sine(frequency = f, amplitude = V)
+    @named ground = Ground()
+
+    # Connections
+    connections = [connect(ac.output, source.V)
+                   connect(source.p, diode.p)
+                   connect(diode.n, resistor.p)
+                   connect(resistor.n, capacitor.p)
+                   connect(capacitor.n, source.n, ground.g)]
+
+    # Model
+    @named model = ODESystem(connections, t;
+        systems = [resistor, capacitor, source, diode, ac, ground])
+    sys = structural_simplify(model)
+    prob = ODEProblem(sys, Pair[], (0.0, 10.0))
+    sol = solve(prob)
+
+    # Extract solutions for testing
+    diode_voltage = sol[diode.v]
+    diode_current = sol[diode.i]
+    resistor_current = sol[resistor.i]
+    capacitor_voltage = sol[capacitor.v]
+
+    # Tests
+    @test all(diode_current .>= -Is)
+    @test capacitor_voltage[end].≈V rtol=3e-1
+
+    # For visual inspection
+    # plt = plot(sol; vars = [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
+    # 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 = HeatingDiode(n = n, Is = Is)
+    @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)
+
+    # Extract solutions for testing
+    diode_voltage = sol[heating_diode.v]
+    diode_current = sol[heating_diode.i]
+    resistor_current = sol[resistor.i]
+    capacitor_voltage = sol[capacitor.v]
+
+    # Expected thermal voltage at given temperature
+    k = 1.380649e-23  # Boltzmann constant (J/K)
+    q = 1.602176634e-19  # Elementary charge (C)
+
+    # Tests
+    @test all(diode_current .>= -Is) # Diode current should not exceed reverse saturation
+    @test capacitor_voltage[end]≈V rtol=3e-1 # Final capacitor voltage close to input voltage
+
+    # For visual inspection
+    # plt = plot(sol; vars = [heating_diode.i, resistor.i, capacitor.v],
+    #     size = (800, 600), dpi = 300,
+    #     labels = ["HeatingDiode Current" "Resistor Current" "Capacitor Voltage"],
+    #     title = "HeatingDiode Test")
+    # savefig(plt, "heating_diode_test")
+
+    # Remake model with higher amb. temperature, final capacitor voltage should be lower
+    T = 400.0
+    model = remake(prob; p = [temp.T => T])
+    sol = solve(model)
+    @test SciMLBase.successful_retcode(sol)
+    @test sol[capacitor.v][end] < capacitor_voltage[end]
+end