diff --git a/docs/src/API/electrical.md b/docs/src/API/electrical.md index 9cc0204f2..7f6546abd 100644 --- a/docs/src/API/electrical.md +++ b/docs/src/API/electrical.md @@ -31,6 +31,8 @@ Conductor Capacitor Inductor IdealOpAmp +Diode +HeatingDiode ``` ## Analog Sensors diff --git a/src/Electrical/Analog/ideal_components.jl b/src/Electrical/Analog/ideal_components.jl index 4ff515fbb..a90d9258a 100644 --- a/src/Electrical/Analog/ideal_components.jl +++ b/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 \ No newline at end of file diff --git a/src/Electrical/Electrical.jl b/src/Electrical/Electrical.jl index 876858197..55c1acd93 100644 --- a/src/Electrical/Electrical.jl +++ b/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 diff --git a/test/Electrical/analog.jl b/test/Electrical/analog.jl index 1ac86ef73..befd16381 100644 --- a/test/Electrical/analog.jl +++ b/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