From 7156a5946631c35cb94947f84a6992d1338ad0ba Mon Sep 17 00:00:00 2001 From: Fredrik Ahlgren Date: Thu, 22 Oct 2026 08:51:38 +0200 Subject: [PATCH 1/4] feat(fuse): the fuse model as an importable package The IEC 60269 curves and the per-phase damage rule lived in package main, so FTW's fast closed-loop simulation (go/internal/sitesim) kept its own port of the gG curve. They move to package fuse at the module root, with no dependencies: Lookup, All, Profile.TripSeconds and Profile.Step on a Phase. The site loop now calls it; its JSON and behaviour are unchanged. Co-Authored-By: Claude Opus 5.5 --- CLAUDE.md | 3 + cmd/simulator/fuse.go | 239 ++++--------------------------------- cmd/simulator/fuse_test.go | 24 ++++ fuse/fuse.go | 196 ++++++++++++++++++++++++++++++ fuse/fuse_test.go | 116 ++++++++++++++++++ 5 files changed, 362 insertions(+), 216 deletions(-) create mode 100644 cmd/simulator/fuse_test.go create mode 100644 fuse/fuse.go create mode 100644 fuse/fuse_test.go diff --git a/CLAUDE.md b/CLAUDE.md index 98cc62d..aab2f81 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -86,6 +86,9 @@ The app supports two modes controlled by the `--desktop` flag: - `cmd/simulator/templates/index.html` - Embedded dashboard UI - `cmd/simulator/wails.json` - Wails configuration +**Public package:** +- `fuse/` - IEC 60269 fuse curves and per-phase damage; no dependencies, so an EMS's own simulation (FTW's `go/internal/sitesim`) imports the same model + **Internal packages:** - `internal/modbus/server.go` - Modbus TCP server with register storage and logging - `internal/modbus/registers.go` - Register definitions, encoding/decoding (U16, I16, U32, I32) diff --git a/cmd/simulator/fuse.go b/cmd/simulator/fuse.go index aa9f36c..feba551 100644 --- a/cmd/simulator/fuse.go +++ b/cmd/simulator/fuse.go @@ -3,38 +3,22 @@ package main import ( "math" "time" + + "github.com/srcfl/device-simulator/fuse" ) +// The fuse model lives in package fuse, which an EMS's own simulation can +// import; these names keep the site code and its JSON as they were. const ( - FuseProfileManual = "manual" - FuseProfileIecGg = "iec_gg" - FuseProfileIecGm = "iec_gm" - FuseProfileIecAm = "iec_am" - FuseProfileIecGr = "iec_gr" - FuseProfileIecGs = "iec_gs" + FuseProfileManual = fuse.Manual + FuseProfileIecGg = fuse.IecGg + FuseProfileIecGm = fuse.IecGm + FuseProfileIecAm = fuse.IecAm + FuseProfileIecGr = fuse.IecGr + FuseProfileIecGs = fuse.IecGs ) -type FuseCurvePoint struct { - Multiple float64 - Seconds float64 -} - -type FuseProfileDefinition struct { - ID string - Label string - Curve []FuseCurvePoint - CooldownSeconds float64 - MinTripSeconds float64 - TripStartRatio float64 -} - -type FusePhaseOverload struct { - Ratio float64 `json:"ratio"` - TripSeconds float64 `json:"trip_s"` - RemainingSeconds float64 `json:"remaining_s"` - Damage float64 `json:"damage"` - Active bool `json:"active"` -} +type FusePhaseOverload = fuse.Status type FuseOverload struct { L1 FusePhaseOverload `json:"l1"` @@ -42,141 +26,21 @@ type FuseOverload struct { L3 FusePhaseOverload `json:"l3"` } -var fuseProfiles = map[string]FuseProfileDefinition{ - FuseProfileManual: { - ID: FuseProfileManual, - Label: "Manual", - }, - FuseProfileIecGg: { - ID: FuseProfileIecGg, - Label: "IEC 60269 gG/gL (general purpose)", - Curve: []FuseCurvePoint{ - {Multiple: 1.25, Seconds: 3600}, - {Multiple: 1.6, Seconds: 1800}, - {Multiple: 2.0, Seconds: 600}, - {Multiple: 3.0, Seconds: 120}, - {Multiple: 5.0, Seconds: 20}, - {Multiple: 10.0, Seconds: 3}, - {Multiple: 20.0, Seconds: 0.5}, - }, - CooldownSeconds: 300, - MinTripSeconds: 0.05, - TripStartRatio: 1, - }, - FuseProfileIecGm: { - ID: FuseProfileIecGm, - Label: "IEC 60269 gM (motor, dual-rated) - UK common", - Curve: []FuseCurvePoint{ - {Multiple: 1.25, Seconds: 7200}, - {Multiple: 1.6, Seconds: 3600}, - {Multiple: 2.0, Seconds: 1200}, - {Multiple: 3.0, Seconds: 240}, - {Multiple: 5.0, Seconds: 40}, - {Multiple: 10.0, Seconds: 6}, - {Multiple: 20.0, Seconds: 1}, - }, - CooldownSeconds: 420, - MinTripSeconds: 0.05, - TripStartRatio: 1, - }, - FuseProfileIecAm: { - ID: FuseProfileIecAm, - Label: "IEC 60269 aM (motor, short-circuit)", - // Partial-range curve: trip only at higher multiples. - Curve: []FuseCurvePoint{ - {Multiple: 4.0, Seconds: 60}, - {Multiple: 6.3, Seconds: 60}, - {Multiple: 8.0, Seconds: 0.5}, - {Multiple: 10.0, Seconds: 0.2}, - {Multiple: 19.0, Seconds: 0.18}, - }, - CooldownSeconds: 600, - MinTripSeconds: 0.05, - TripStartRatio: 4, - }, - FuseProfileIecGr: { - ID: FuseProfileIecGr, - Label: "IEC 60269 gR (semiconductor, fast)", - Curve: []FuseCurvePoint{ - {Multiple: 1.25, Seconds: 1080}, - {Multiple: 1.6, Seconds: 540}, - {Multiple: 2.0, Seconds: 180}, - {Multiple: 3.0, Seconds: 36}, - {Multiple: 5.0, Seconds: 6}, - {Multiple: 10.0, Seconds: 0.9}, - {Multiple: 20.0, Seconds: 0.15}, - }, - CooldownSeconds: 180, - MinTripSeconds: 0.02, - TripStartRatio: 1, - }, - FuseProfileIecGs: { - ID: FuseProfileIecGs, - Label: "IEC 60269 gS (semiconductor, full-range)", - Curve: []FuseCurvePoint{ - {Multiple: 1.25, Seconds: 1800}, - {Multiple: 1.6, Seconds: 900}, - {Multiple: 2.0, Seconds: 300}, - {Multiple: 3.0, Seconds: 60}, - {Multiple: 5.0, Seconds: 10}, - {Multiple: 10.0, Seconds: 1.5}, - {Multiple: 20.0, Seconds: 0.25}, - }, - CooldownSeconds: 240, - MinTripSeconds: 0.03, - TripStartRatio: 1, - }, -} +var fuseProfiles = func() map[string]fuse.Profile { + m := map[string]fuse.Profile{} + for _, p := range fuse.All() { + m[p.ID] = p + } + return m +}() -func resolveFuseProfile(id string) FuseProfileDefinition { +func resolveFuseProfile(id string) fuse.Profile { if profile, ok := fuseProfiles[id]; ok { return profile } return fuseProfiles[FuseProfileManual] } -func fuseTripTimeSeconds(profile FuseProfileDefinition, multiple float64) float64 { - if multiple <= 1 { - return math.Inf(1) - } - if len(profile.Curve) == 0 { - return math.Inf(1) - } - - points := profile.Curve - - if multiple <= points[0].Multiple { - return points[0].Seconds - } - - for i := 1; i < len(points); i++ { - if multiple <= points[i].Multiple { - return logLogInterpolate(multiple, points[i-1], points[i]) - } - } - - // Extrapolate using the last segment slope for very high currents. - last := points[len(points)-1] - prev := points[len(points)-2] - return logLogInterpolate(multiple, prev, last) -} - -func logLogInterpolate(x float64, a, b FuseCurvePoint) float64 { - if x <= 0 || a.Multiple <= 0 || b.Multiple <= 0 || a.Seconds <= 0 || b.Seconds <= 0 { - return b.Seconds - } - lx := math.Log10(x) - lx1 := math.Log10(a.Multiple) - lx2 := math.Log10(b.Multiple) - if lx2 == lx1 { - return b.Seconds - } - ly1 := math.Log10(a.Seconds) - ly2 := math.Log10(b.Seconds) - t := (lx - lx1) / (lx2 - lx1) - return math.Pow(10, ly1+(ly2-ly1)*t) -} - func updateFuseState(site *Site, phasePowers PhasePowers, simDelta time.Duration) PhaseFuses { simSeconds := simDelta.Seconds() if simSeconds < 0 { @@ -211,67 +75,10 @@ func updateFuseState(site *Site, phasePowers PhasePowers, simDelta time.Duration overload := FuseOverload{} updatePhase := func(power float64, blown bool, damage *float64) (FusePhaseOverload, bool) { - if blown { - if *damage < 1 { - *damage = 1 - } - return FusePhaseOverload{Damage: *damage}, false - } - - current := math.Abs(power) / 230.0 - ratio := 0.0 - if fuseSize > 0 { - ratio = current / float64(fuseSize) - } - - startRatio := profile.TripStartRatio - if startRatio <= 0 { - startRatio = 1 - } - - if ratio <= startRatio { - if profile.CooldownSeconds > 0 && simSeconds > 0 { - *damage -= simSeconds / profile.CooldownSeconds - if *damage < 0 { - *damage = 0 - } - } - return FusePhaseOverload{Ratio: ratio, Damage: *damage, Active: false}, false - } - - tripSeconds := fuseTripTimeSeconds(profile, ratio) - if tripSeconds <= 0 || math.IsInf(tripSeconds, 1) { - return FusePhaseOverload{Ratio: ratio, Damage: *damage, Active: false}, false - } - if profile.MinTripSeconds > 0 && tripSeconds < profile.MinTripSeconds { - tripSeconds = profile.MinTripSeconds - } - - if simSeconds > 0 { - *damage += simSeconds / tripSeconds - if *damage > 1 { - *damage = 1 - } - } - - if *damage >= 1 { - return FusePhaseOverload{ - Ratio: ratio, - TripSeconds: tripSeconds, - RemainingSeconds: 0, - Damage: *damage, - Active: false, - }, true - } - - remaining := (1 - *damage) * tripSeconds - return FusePhaseOverload{ - Ratio: ratio, - TripSeconds: tripSeconds, - RemainingSeconds: remaining, - Damage: *damage, - Active: true, - }, false + ph := fuse.Phase{Damage: *damage, Blown: blown} + status := profile.Step(&ph, math.Abs(power)/230.0, float64(fuseSize), simSeconds) + *damage = ph.Damage + return status, ph.Blown && !blown } status, blew := updatePhase(phasePowers.L1, fuseBlown.L1, &damage.L1) diff --git a/cmd/simulator/fuse_test.go b/cmd/simulator/fuse_test.go new file mode 100644 index 0000000..ecd8cd2 --- /dev/null +++ b/cmd/simulator/fuse_test.go @@ -0,0 +1,24 @@ +package main + +import ( + "testing" + "time" +) + +// The site's fuse is package fuse's iec_gg: 2 × a 25 A rating on L1 blows +// it after 600 s, and the meter's other phases stay whole. +func TestSiteFuseBlowsOnTheGGCurve(t *testing.T) { + site := &Site{FuseSize: 25, FuseProfile: FuseProfileIecGg} + at := 0 + for s := 1; s <= 700 && at == 0; s++ { + if updateFuseState(site, PhasePowers{L1: 50 * 230, L2: 1000, L3: -1000}, time.Second).L1 { + at = s + } + } + if at != 600 || site.FuseBlown.L2 || site.FuseBlown.L3 { + t.Fatalf("L1 blew after %d s (want 600), L2/L3 %v/%v", at, site.FuseBlown.L2, site.FuseBlown.L3) + } + if site.FuseOverload.L1.Damage != 1 || site.FuseOverload.L2.Active { + t.Fatalf("overload %+v", site.FuseOverload) + } +} diff --git a/fuse/fuse.go b/fuse/fuse.go new file mode 100644 index 0000000..939b193 --- /dev/null +++ b/fuse/fuse.go @@ -0,0 +1,196 @@ +// Package fuse models a fuse on each phase of a site's main supply: the time +// it takes to blow at a multiple of its rating (IEC 60269 time-current +// curves) and the damage a phase builds up and sheds as the current moves. +// +// Every second above the rating spends 1/t of the fuse's life, with t the +// curve's time to blow at that multiple; below it the damage recovers over +// the profile's cooldown. A phase blows when its damage reaches 1. The curve +// blows slowly from just above the rating, so it is stricter than a real gG +// fuse below 1.25 × its rating. +// +// It has no dependencies, so an EMS's own fast simulation can import the +// same curves the device simulator uses. +package fuse + +import ( + "math" + "sort" +) + +// Profile ids. +const ( + Manual = "manual" // never blows by itself + IecGg = "iec_gg" + IecGm = "iec_gm" + IecAm = "iec_am" + IecGr = "iec_gr" + IecGs = "iec_gs" +) + +// Point is one point of a time-current curve: at Multiple × the rating the +// fuse blows after Seconds. +type Point struct { + Multiple float64 + Seconds float64 +} + +// Profile is a fuse type. +type Profile struct { + ID string + Label string + Curve []Point + // CooldownS is how long damage 1 takes to recover at or below the start ratio. + CooldownS float64 + // MinTripS is the shortest time to blow. + MinTripS float64 + // TripStartRatio is the multiple of the rating above which damage builds. + TripStartRatio float64 +} + +var profiles = map[string]Profile{ + Manual: {ID: Manual, Label: "Manual"}, + IecGg: { + ID: IecGg, + Label: "IEC 60269 gG/gL (general purpose)", + Curve: []Point{ + {1.25, 3600}, {1.6, 1800}, {2.0, 600}, {3.0, 120}, {5.0, 20}, {10.0, 3}, {20.0, 0.5}, + }, + CooldownS: 300, MinTripS: 0.05, TripStartRatio: 1, + }, + IecGm: { + ID: IecGm, + Label: "IEC 60269 gM (motor, dual-rated) - UK common", + Curve: []Point{ + {1.25, 7200}, {1.6, 3600}, {2.0, 1200}, {3.0, 240}, {5.0, 40}, {10.0, 6}, {20.0, 1}, + }, + CooldownS: 420, MinTripS: 0.05, TripStartRatio: 1, + }, + IecAm: { + ID: IecAm, + Label: "IEC 60269 aM (motor, short-circuit)", + // Partial-range curve: trips only at higher multiples. + Curve: []Point{ + {4.0, 60}, {6.3, 60}, {8.0, 0.5}, {10.0, 0.2}, {19.0, 0.18}, + }, + CooldownS: 600, MinTripS: 0.05, TripStartRatio: 4, + }, + IecGr: { + ID: IecGr, + Label: "IEC 60269 gR (semiconductor, fast)", + Curve: []Point{ + {1.25, 1080}, {1.6, 540}, {2.0, 180}, {3.0, 36}, {5.0, 6}, {10.0, 0.9}, {20.0, 0.15}, + }, + CooldownS: 180, MinTripS: 0.02, TripStartRatio: 1, + }, + IecGs: { + ID: IecGs, + Label: "IEC 60269 gS (semiconductor, full-range)", + Curve: []Point{ + {1.25, 1800}, {1.6, 900}, {2.0, 300}, {3.0, 60}, {5.0, 10}, {10.0, 1.5}, {20.0, 0.25}, + }, + CooldownS: 240, MinTripS: 0.03, TripStartRatio: 1, + }, +} + +// Lookup returns the profile with this id. +func Lookup(id string) (Profile, bool) { + p, ok := profiles[id] + return p, ok +} + +// All returns every profile, ordered by id. +func All() []Profile { + out := make([]Profile, 0, len(profiles)) + for _, p := range profiles { + out = append(out, p) + } + sort.Slice(out, func(i, j int) bool { return out[i].ID < out[j].ID }) + return out +} + +// TripSeconds is the time to blow at multiple × the rating, interpolated on +// log-log axes, with the last segment extended above its end. It is +Inf at +// or below 1 ×, and for a profile without a curve. +func (p Profile) TripSeconds(multiple float64) float64 { + if multiple <= 1 || len(p.Curve) == 0 { + return math.Inf(1) + } + pts := p.Curve + if multiple <= pts[0].Multiple { + return pts[0].Seconds + } + for i := 1; i < len(pts); i++ { + if multiple <= pts[i].Multiple { + return logLog(multiple, pts[i-1], pts[i]) + } + } + return logLog(multiple, pts[len(pts)-2], pts[len(pts)-1]) +} + +func logLog(x float64, a, b Point) float64 { + if x <= 0 || a.Multiple <= 0 || b.Multiple <= 0 || a.Seconds <= 0 || b.Seconds <= 0 { + return b.Seconds + } + lx, lx1, lx2 := math.Log10(x), math.Log10(a.Multiple), math.Log10(b.Multiple) + if lx2 == lx1 { + return b.Seconds + } + ly1, ly2 := math.Log10(a.Seconds), math.Log10(b.Seconds) + return math.Pow(10, ly1+(ly2-ly1)*(lx-lx1)/(lx2-lx1)) +} + +// Phase is the state of one phase's fuse. +type Phase struct { + Damage float64 // 0 new, 1 blown + Blown bool +} + +// Status describes a phase after a step. +type Status struct { + Ratio float64 `json:"ratio"` // current over the rating + TripSeconds float64 `json:"trip_s"` // time to blow at this ratio from new + RemainingSeconds float64 `json:"remaining_s"` + Damage float64 `json:"damage"` + Active bool `json:"active"` // building damage +} + +// Step advances one phase by dt seconds at amps through a fuse rated +// ratingA. The sign of amps does not matter: export heats a fuse too. A +// profile without a curve leaves the phase as it is. +func (p Profile) Step(ph *Phase, amps, ratingA, dt float64) Status { + if len(p.Curve) == 0 { + return Status{Damage: ph.Damage} + } + if dt < 0 { + dt = 0 + } + if ph.Blown { + ph.Damage = math.Max(ph.Damage, 1) + return Status{Damage: ph.Damage} + } + ratio := 0.0 + if ratingA > 0 { + ratio = math.Abs(amps) / ratingA + } + start := p.TripStartRatio + if start <= 0 { + start = 1 + } + if ratio <= start { + if p.CooldownS > 0 { + ph.Damage = math.Max(0, ph.Damage-dt/p.CooldownS) + } + return Status{Ratio: ratio, Damage: ph.Damage} + } + trip := p.TripSeconds(ratio) + if trip <= 0 || math.IsInf(trip, 1) { + return Status{Ratio: ratio, Damage: ph.Damage} + } + trip = math.Max(trip, p.MinTripS) + ph.Damage = math.Min(1, ph.Damage+dt/trip) + if ph.Damage >= 1 { + ph.Blown = true + return Status{Ratio: ratio, TripSeconds: trip, Damage: ph.Damage} + } + return Status{Ratio: ratio, TripSeconds: trip, RemainingSeconds: (1 - ph.Damage) * trip, Damage: ph.Damage, Active: true} +} diff --git a/fuse/fuse_test.go b/fuse/fuse_test.go new file mode 100644 index 0000000..99b666e --- /dev/null +++ b/fuse/fuse_test.go @@ -0,0 +1,116 @@ +package fuse + +import ( + "math" + "testing" +) + +func gg(t *testing.T) Profile { + t.Helper() + p, ok := Lookup(IecGg) + if !ok { + t.Fatal("no iec_gg profile") + } + return p +} + +func TestTripSecondsFollowsTheCurve(t *testing.T) { + p := gg(t) + for _, pt := range p.Curve { + if got := p.TripSeconds(pt.Multiple); math.Abs(got-pt.Seconds) > 1e-9*pt.Seconds { + t.Errorf("%.2f×: %v s, want %v s", pt.Multiple, got, pt.Seconds) + } + } + if got := p.TripSeconds(1); !math.IsInf(got, 1) { + t.Errorf("1×: %v s, want never", got) + } + if got := p.TripSeconds(1.1); got != 3600 { + t.Errorf("1.1×: %v s, want the first point's 3600 s", got) + } + if got := p.TripSeconds(1.4); got <= 1800 || got >= 3600 { + t.Errorf("1.4×: %.0f s, want between the 1.25× and 1.6× points", got) + } + if got := p.TripSeconds(40); got >= 0.5 { + t.Errorf("40×: %v s, want under the 20× point", got) + } +} + +// A phase held at a multiple of its rating blows after the curve's time, +// in one-second steps; the other phases are untouched. Export heats a fuse +// as much as import. +func TestStepBlowsOnTheCurve(t *testing.T) { + p := gg(t) + for _, tc := range []struct{ multiple, wantS float64 }{{2, 600}, {5, 20}, {10, 3}} { + var phases [3]Phase + blewAt := 0 + for s := 1; s <= 4000 && blewAt == 0; s++ { + p.Step(&phases[0], -tc.multiple*25, 25, 1) + p.Step(&phases[1], 0, 25, 1) + if phases[0].Blown { + blewAt = s + } + } + if float64(blewAt) != tc.wantS || phases[1].Damage != 0 || phases[1].Blown { + t.Errorf("%.0f×: blew after %d s, want %v s", tc.multiple, blewAt, tc.wantS) + } + } + var ph Phase + for s := 0; s < 4000; s++ { + p.Step(&ph, 25, 25, 1) + } + if ph.Blown || ph.Damage != 0 { + t.Errorf("at the rating: %+v, want no damage", ph) + } +} + +func TestDamageRecoversBelowTheRating(t *testing.T) { + p := gg(t) + ph := Phase{Damage: 0.5} + st := p.Step(&ph, 10, 25, 75) + if math.Abs(ph.Damage-0.25) > 1e-12 || st.Active { + t.Fatalf("after 75 s at 0.4×: damage %v, want 0.25 (cooldown %v s)", ph.Damage, p.CooldownS) + } + st = p.Step(&ph, 50, 25, 1) + if !st.Active || math.Abs(st.TripSeconds-600) > 1e-9 || math.Abs(st.RemainingSeconds-(1-ph.Damage)*st.TripSeconds) > 1e-9 { + t.Fatalf("at 2×: %+v", st) + } +} + +func TestBlownStaysBlownAndManualNeverBlows(t *testing.T) { + p := gg(t) + ph := Phase{Blown: true, Damage: 0.9} + if st := p.Step(&ph, 0, 25, 1000); !ph.Blown || ph.Damage != 1 || st.Damage != 1 { + t.Fatalf("blown phase after a cool second: %+v", ph) + } + m, _ := Lookup(Manual) + ph = Phase{} + m.Step(&ph, 1000, 25, 1000) + if ph.Blown || ph.Damage != 0 { + t.Fatalf("manual: %+v, want untouched", ph) + } +} + +func TestAMBuildsOnlyFromItsStartRatio(t *testing.T) { + p, _ := Lookup(IecAm) + ph := Phase{} + p.Step(&ph, 3.9*25, 25, 100) + if ph.Damage != 0 { + t.Fatalf("aM at 3.9×: damage %v, want none", ph.Damage) + } + p.Step(&ph, 5*25, 25, 1) + if ph.Damage == 0 { + t.Fatal("aM at 5×: no damage") + } +} + +func TestAllIsSortedAndComplete(t *testing.T) { + all := All() + if len(all) != 6 { + t.Fatalf("%d profiles", len(all)) + } + for i := 1; i < len(all); i++ { + if all[i-1].ID >= all[i].ID { + t.Fatalf("not sorted: %s before %s", all[i-1].ID, all[i].ID) + } + } +} From a6baf29be5116c53b536b5fc3d9792827df55e37 Mon Sep 17 00:00:00 2001 From: Fredrik Ahlgren Date: Thu, 22 Oct 2026 19:52:12 +0200 Subject: [PATCH 2/4] feat(fuse): a thermal model calibrated to IEC 60269 gates MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The fuse was a damage counter: 1/t of its life per second above the rating, recovered linearly over a fixed cooldown. A fuse that had carried its rating for an hour blew no sooner than a cold one, and 1.13 × blew within the hour where IEC 60269-1 lets a gG carry 1.25 × for it. Now each phase is three thermal bodies in a ladder (element, fuse body and sand, holder and cabinet) heated by I² and melting at tin's 232 °C, with the ambient as a parameter. Each gG band (16, 20, 25, 32, 35, 40, 50, 63 A) is fitted to IEC 60269-1's gates and pre-arcing I²t at 0.01 s and keeps the conventional currents; aM to its gates. gM, gR and gS are gG with scaled times and say they are not calibrated. Steps are exact, so the result does not depend on the step size. The site keeps each phase's heat between ticks and clears it with the fuse. Co-Authored-By: Claude Opus 5.5 --- cmd/simulator/app.go | 9 +- cmd/simulator/fuse.go | 15 +- cmd/simulator/fuse_test.go | 41 ++- cmd/simulator/main.go | 15 +- cmd/simulator/templates/index.html | 4 + fuse/fuse.go | 536 ++++++++++++++++++++++------- fuse/fuse_test.go | 284 +++++++++++---- 7 files changed, 698 insertions(+), 206 deletions(-) diff --git a/cmd/simulator/app.go b/cmd/simulator/app.go index aabbeb8..62383c7 100644 --- a/cmd/simulator/app.go +++ b/cmd/simulator/app.go @@ -7,6 +7,7 @@ import ( "math" "time" + "github.com/srcfl/device-simulator/fuse" "github.com/srcfl/device-simulator/internal/modbus" "github.com/srcfl/device-simulator/internal/modbus/devices" "github.com/srcfl/device-simulator/internal/settings" @@ -909,6 +910,7 @@ func (a *App) UpdateSiteConfig(id int, config map[string]interface{}) error { if fuseSize, ok := config["fuse_size"].(float64); ok { site.FuseSize = int(fuseSize) site.FuseDamage = PhaseLoads{} + site.FuseHeat = [3]fuse.Phase{} site.FuseOverload = FuseOverload{} } if fuseProfile, ok := config["fuse_profile"].(string); ok { @@ -918,6 +920,7 @@ func (a *App) UpdateSiteConfig(id int, config map[string]interface{}) error { if site.FuseProfile != fuseProfile { site.FuseProfile = fuseProfile site.FuseDamage = PhaseLoads{} + site.FuseHeat = [3]fuse.Phase{} site.FuseOverload = FuseOverload{} } } @@ -935,19 +938,19 @@ func (a *App) UpdateSiteConfig(id int, config map[string]interface{}) error { if fb, ok := config["fuse_blown"].(map[string]interface{}); ok { if l1, ok := parseBool(fb["l1"]); ok { if site.FuseBlown.L1 && !l1 { - site.FuseDamage.L1 = 0 + site.FuseDamage.L1, site.FuseHeat[0] = 0, fuse.Phase{} } site.FuseBlown.L1 = l1 } if l2, ok := parseBool(fb["l2"]); ok { if site.FuseBlown.L2 && !l2 { - site.FuseDamage.L2 = 0 + site.FuseDamage.L2, site.FuseHeat[1] = 0, fuse.Phase{} } site.FuseBlown.L2 = l2 } if l3, ok := parseBool(fb["l3"]); ok { if site.FuseBlown.L3 && !l3 { - site.FuseDamage.L3 = 0 + site.FuseDamage.L3, site.FuseHeat[2] = 0, fuse.Phase{} } site.FuseBlown.L3 = l3 } diff --git a/cmd/simulator/fuse.go b/cmd/simulator/fuse.go index feba551..9bda46e 100644 --- a/cmd/simulator/fuse.go +++ b/cmd/simulator/fuse.go @@ -74,28 +74,31 @@ func updateFuseState(site *Site, phasePowers PhasePowers, simDelta time.Duration fuseBlown := site.FuseBlown overload := FuseOverload{} - updatePhase := func(power float64, blown bool, damage *float64) (FusePhaseOverload, bool) { - ph := fuse.Phase{Damage: *damage, Blown: blown} - status := profile.Step(&ph, math.Abs(power)/230.0, float64(fuseSize), simSeconds) + // Each phase keeps its thermal state between ticks, so a fuse that has + // been warm blows sooner. + updatePhase := func(i int, power float64, blown bool, damage *float64) (FusePhaseOverload, bool) { + ph := &site.FuseHeat[i] + ph.Blown, ph.AmbientC = blown, site.FuseAmbientC + status := profile.Step(ph, math.Abs(power)/230.0, float64(fuseSize), simSeconds) *damage = ph.Damage return status, ph.Blown && !blown } - status, blew := updatePhase(phasePowers.L1, fuseBlown.L1, &damage.L1) + status, blew := updatePhase(0, phasePowers.L1, fuseBlown.L1, &damage.L1) overload.L1 = status if blew { fuseBlown.L1 = true changed = true } - status, blew = updatePhase(phasePowers.L2, fuseBlown.L2, &damage.L2) + status, blew = updatePhase(1, phasePowers.L2, fuseBlown.L2, &damage.L2) overload.L2 = status if blew { fuseBlown.L2 = true changed = true } - status, blew = updatePhase(phasePowers.L3, fuseBlown.L3, &damage.L3) + status, blew = updatePhase(2, phasePowers.L3, fuseBlown.L3, &damage.L3) overload.L3 = status if blew { fuseBlown.L3 = true diff --git a/cmd/simulator/fuse_test.go b/cmd/simulator/fuse_test.go index ecd8cd2..9cb3784 100644 --- a/cmd/simulator/fuse_test.go +++ b/cmd/simulator/fuse_test.go @@ -1,24 +1,59 @@ package main import ( + "math" "testing" "time" + + "github.com/srcfl/device-simulator/fuse" ) // The site's fuse is package fuse's iec_gg: 2 × a 25 A rating on L1 blows -// it after 600 s, and the meter's other phases stay whole. +// it after the band's cold time in the cabinet's 25 °C, and the meter's +// other phases stay whole. func TestSiteFuseBlowsOnTheGGCurve(t *testing.T) { site := &Site{FuseSize: 25, FuseProfile: FuseProfileIecGg} + gg, _ := fuse.Lookup(fuse.IecGg) + want := gg.MeltSeconds(2, 25, fuse.DefaultAmbientC) at := 0 for s := 1; s <= 700 && at == 0; s++ { if updateFuseState(site, PhasePowers{L1: 50 * 230, L2: 1000, L3: -1000}, time.Second).L1 { at = s } } - if at != 600 || site.FuseBlown.L2 || site.FuseBlown.L3 { - t.Fatalf("L1 blew after %d s (want 600), L2/L3 %v/%v", at, site.FuseBlown.L2, site.FuseBlown.L3) + if math.Abs(float64(at)-want) > 1 || site.FuseBlown.L2 || site.FuseBlown.L3 { + t.Fatalf("L1 blew after %d s (want %.1f), L2/L3 %v/%v", at, want, site.FuseBlown.L2, site.FuseBlown.L3) } if site.FuseOverload.L1.Damage != 1 || site.FuseOverload.L2.Active { t.Fatalf("overload %+v", site.FuseOverload) } } + +// The site keeps each phase's heat between ticks: a phase that has carried +// its rating for an hour blows at 2 × sooner than a cold one, and clearing +// a blown fuse starts it cold. +func TestSiteFuseRemembersItsHeat(t *testing.T) { + blowAt := func(site *Site) int { + for s := 1; s <= 700; s++ { + if updateFuseState(site, PhasePowers{L1: 50 * 230}, time.Second).L1 { + return s + } + } + return 0 + } + cold := blowAt(&Site{FuseSize: 25, FuseProfile: FuseProfileIecGg}) + warm := &Site{FuseSize: 25, FuseProfile: FuseProfileIecGg} + for s := 0; s < 3600; s++ { + updateFuseState(warm, PhasePowers{L1: 25 * 230}, time.Second) + } + if d := warm.FuseDamage.L1; d < 0.4 || d > 0.6 { + t.Fatalf("after 1 h at the rating damage %.2f, want about half way to melting", d) + } + if hot := blowAt(warm); hot == 0 || hot >= cold*2/3 { + t.Fatalf("2× after 1 h at the rating blew after %d s, cold after %d s", hot, cold) + } + warm.FuseBlown.L1, warm.FuseDamage.L1, warm.FuseHeat[0] = false, 0, fuse.Phase{} + if again := blowAt(warm); again != cold { + t.Fatalf("a cleared fuse blew after %d s, want the cold %d s", again, cold) + } +} diff --git a/cmd/simulator/main.go b/cmd/simulator/main.go index 257552f..1fb4932 100644 --- a/cmd/simulator/main.go +++ b/cmd/simulator/main.go @@ -20,6 +20,7 @@ import ( "syscall" "time" + "github.com/srcfl/device-simulator/fuse" "github.com/srcfl/device-simulator/internal/device" simmdns "github.com/srcfl/device-simulator/internal/mdns" "github.com/srcfl/device-simulator/internal/modbus" @@ -143,6 +144,8 @@ type Site struct { FuseBlown PhaseFuses `json:"fuse_blown"` // Per-phase blown fuse state FuseDamage PhaseLoads `json:"-"` FuseOverload FuseOverload `json:"-"` + FuseHeat [3]fuse.Phase `json:"-"` // each phase's thermal state (L1, L2, L3) + FuseAmbientC float64 `json:"-"` // the fuse cabinet's air; 0 is 25 °C BoundMeter bool `json:"bound_meter"` // True if meter was auto-created with site P1Buffer *P1IntervalBuffer // P1 port interval buffer @@ -2642,14 +2645,15 @@ func handleSiteConfig(w http.ResponseWriter, r *http.Request, site *Site) { } if config.FuseSize != nil { // Validate fuse size (common values: 10, 16, 20, 25, 32 Amps) - validFuses := map[int]bool{10: true, 16: true, 20: true, 25: true, 32: true} + validFuses := map[int]bool{10: true, 16: true, 20: true, 25: true, 32: true, 35: true, 40: true, 50: true, 63: true} if !validFuses[*config.FuseSize] { site.mu.Unlock() - http.Error(w, "Fuse size must be 10, 16, 20, 25, or 32 Amps", http.StatusBadRequest) + http.Error(w, "Fuse size must be 10, 16, 20, 25, 32, 35, 40, 50 or 63 Amps", http.StatusBadRequest) return } site.FuseSize = *config.FuseSize site.FuseDamage = PhaseLoads{} + site.FuseHeat = [3]fuse.Phase{} site.FuseOverload = FuseOverload{} } if config.FuseProfile != nil { @@ -2661,6 +2665,7 @@ func handleSiteConfig(w http.ResponseWriter, r *http.Request, site *Site) { if site.FuseProfile != *config.FuseProfile { site.FuseProfile = *config.FuseProfile site.FuseDamage = PhaseLoads{} + site.FuseHeat = [3]fuse.Phase{} site.FuseOverload = FuseOverload{} } } @@ -2684,13 +2689,13 @@ func handleSiteConfig(w http.ResponseWriter, r *http.Request, site *Site) { } if config.FuseBlown != nil { if site.FuseBlown.L1 && !config.FuseBlown.L1 { - site.FuseDamage.L1 = 0 + site.FuseDamage.L1, site.FuseHeat[0] = 0, fuse.Phase{} } if site.FuseBlown.L2 && !config.FuseBlown.L2 { - site.FuseDamage.L2 = 0 + site.FuseDamage.L2, site.FuseHeat[1] = 0, fuse.Phase{} } if site.FuseBlown.L3 && !config.FuseBlown.L3 { - site.FuseDamage.L3 = 0 + site.FuseDamage.L3, site.FuseHeat[2] = 0, fuse.Phase{} } site.FuseBlown = *config.FuseBlown } diff --git a/cmd/simulator/templates/index.html b/cmd/simulator/templates/index.html index d2571e1..ae568e9 100644 --- a/cmd/simulator/templates/index.html +++ b/cmd/simulator/templates/index.html @@ -322,6 +322,10 @@

Site Settings

+ + + +
diff --git a/fuse/fuse.go b/fuse/fuse.go index 939b193..b9cdc3f 100644 --- a/fuse/fuse.go +++ b/fuse/fuse.go @@ -1,15 +1,29 @@ -// Package fuse models a fuse on each phase of a site's main supply: the time -// it takes to blow at a multiple of its rating (IEC 60269 time-current -// curves) and the damage a phase builds up and sheds as the current moves. +// Package fuse models the fuse on each phase of a site's main supply as a +// thermal system: the current heats the fuse element (I², so export heats +// it as much as import), the heat flows through the fuse body and sand to +// the holder and the cabinet, and the element melts when it reaches its +// melting temperature. A fuse that has carried its rating for an hour is +// hot and blows at an overload sooner than a cold one; short overloads add +// up; a fuse cools along the sum of its bodies' exponentials. // -// Every second above the rating spends 1/t of the fuse's life, with t the -// curve's time to blow at that multiple; below it the damage recovers over -// the profile's cooldown. A phase blows when its damage reaches 1. The curve -// blows slowly from just above the rating, so it is stricter than a real gG -// fuse below 1.25 × its rating. +// Each phase has three thermal bodies in a ladder: // -// It has no dependencies, so an EMS's own fast simulation can import the -// same curves the device simulator uses. +// heat I² → element → fuse body and sand → holder and cabinet → ambient +// +// In per-unit terms the heat into the element is u² (u = I/In) and the +// three thermal resistances add up to 1, so at a steady current the element +// sits u² above ambient. It melts at UMin², UMin being the steady current +// (as a multiple of In) that just melts it at the reference ambient. The +// element and the body set the short times (adiabatic I²t below 0.1 s, the +// body's seconds); the holder sets the long ones. +// +// The gG bands are calibrated per rating to IEC 60269-1's conventional +// currents and gates and the pre-arcing I²t at 0.01 s (see gg below); aM to +// its gates. The other profiles are gG with scaled times and are not +// calibrated. The model is independent of any EMS's own fuse protection: it +// is the judge. +// +// It has no dependencies, so an EMS's own fast simulation can import it. package fuse import ( @@ -27,71 +41,187 @@ const ( IecGs = "iec_gs" ) -// Point is one point of a time-current curve: at Multiple × the rating the -// fuse blows after Seconds. +const ( + // MeltC is the element's melting temperature: tin, 231.93 °C, the + // solder spot (M-effect) at which a gG element melts under overload. + MeltC = 232.0 + // RefAmbientC is the ambient the calibration holds at: IEC 60269 gives + // time-current characteristics at a reference ambient of 20 °C, + // approximately valid at 30 °C (Enel GSCL004 rev. 0, 2016, §6.1). + RefAmbientC = 20.0 + // DefaultAmbientC is a meter cabinet's air when a phase gives none. + DefaultAmbientC = 25.0 +) + +// Point is one point of a time-current curve from cold: at Multiple × the +// rating the fuse melts after Seconds. type Point struct { Multiple float64 Seconds float64 } -// Profile is a fuse type. +// Gate bounds the melting time from cold at a multiple of the rating: with +// Min the fuse must not melt sooner than Seconds, otherwise it must melt +// within Seconds. +type Gate struct { + Multiple float64 + Seconds float64 + Min bool +} + +// Thermal is a band's three bodies. Re, Rb and Rh are the shares of the +// steady temperature rise across element-to-body, body-to-holder and +// holder-to-ambient (they add up to 1); each Tau is that body's heat +// capacity times its resistance onwards. +type Thermal struct { + UMin float64 // steady melting multiple at the reference ambient + Re, TauE float64 // element + Rb, TauB float64 // fuse body and sand + Rh, TauH float64 // holder and cabinet + TimeFactor float64 // scales every time constant; 0 or 1 as calibrated +} + +// Band is a fuse type at one rating: its thermal bodies and the points it +// is checked against. +type Band struct { + RatingA float64 + Thermal Thermal + // Gates are the sourced bounds; I2tMin and I2tMax bound the pre-arcing + // I²t at 0.01 s (A²s). Targets are the centre points the band is fitted + // to. + Gates []Gate + I2tMin, I2tMax float64 + Targets []Point + + sys system +} + +// Profile is a fuse type over its rating bands. type Profile struct { - ID string - Label string + ID string + Label string + Source string // where its numbers come from + Bands []Band + // Curve is the reference band's cold curve, sampled, for display. Curve []Point - // CooldownS is how long damage 1 takes to recover at or below the start ratio. - CooldownS float64 - // MinTripS is the shortest time to blow. - MinTripS float64 - // TripStartRatio is the multiple of the rating above which damage builds. - TripStartRatio float64 -} - -var profiles = map[string]Profile{ - Manual: {ID: Manual, Label: "Manual"}, - IecGg: { - ID: IecGg, - Label: "IEC 60269 gG/gL (general purpose)", - Curve: []Point{ - {1.25, 3600}, {1.6, 1800}, {2.0, 600}, {3.0, 120}, {5.0, 20}, {10.0, 3}, {20.0, 0.5}, - }, - CooldownS: 300, MinTripS: 0.05, TripStartRatio: 1, - }, - IecGm: { - ID: IecGm, - Label: "IEC 60269 gM (motor, dual-rated) - UK common", - Curve: []Point{ - {1.25, 7200}, {1.6, 3600}, {2.0, 1200}, {3.0, 240}, {5.0, 40}, {10.0, 6}, {20.0, 1}, - }, - CooldownS: 420, MinTripS: 0.05, TripStartRatio: 1, - }, - IecAm: { - ID: IecAm, - Label: "IEC 60269 aM (motor, short-circuit)", - // Partial-range curve: trips only at higher multiples. - Curve: []Point{ - {4.0, 60}, {6.3, 60}, {8.0, 0.5}, {10.0, 0.2}, {19.0, 0.18}, - }, - CooldownS: 600, MinTripS: 0.05, TripStartRatio: 4, - }, - IecGr: { - ID: IecGr, - Label: "IEC 60269 gR (semiconductor, fast)", - Curve: []Point{ - {1.25, 1080}, {1.6, 540}, {2.0, 180}, {3.0, 36}, {5.0, 6}, {10.0, 0.9}, {20.0, 0.15}, - }, - CooldownS: 180, MinTripS: 0.02, TripStartRatio: 1, - }, - IecGs: { - ID: IecGs, - Label: "IEC 60269 gS (semiconductor, full-range)", - Curve: []Point{ - {1.25, 1800}, {1.6, 900}, {2.0, 300}, {3.0, 60}, {5.0, 10}, {10.0, 1.5}, {20.0, 0.25}, - }, - CooldownS: 240, MinTripS: 0.03, TripStartRatio: 1, - }, + // RefRatingA is the band TripSeconds and Curve use. + RefRatingA float64 +} + +// The holder and cabinet: a modelling choice, not from a standard. IEC 60269 +// fixes a gG fuse only at 1.25 × In (must not melt in 1 h) and 1.6 × In (must +// melt within 1 h) in the long-time region; these two numbers put a cold +// 1.6 × melt at about 13 minutes. +const ( + holderShare = 0.4 // of the steady rise at a steady current + holderTauS = 1200.0 // s +) + +// gg are the gG bands. Gates, I²t windows and conventional currents from +// IEC 60269-1 as tabled in Schneider Electric, "Complementary technical +// information" CA908036E v3.2 (2018), gate and pre-arcing I²t tables; the +// 25, 40, 50 and 63 A rows also in Enel GSCL004 rev. 0 (2016) Tables 6, 7 +// and 14, and 25 A in Mersen EduPack TM-104 (2012) slide 6. 35 A is not in +// those tables: its gates and I²t are interpolated on log axes between 32 A +// and 40 A. Conventional currents: Inf = 1.25 In and If = 1.6 In with a +// conventional time of 1 h for 16 ≤ In ≤ 63 A. +// +// Re, TauE and TauB are fitted (least squares on log time, from cold at the +// reference ambient) to three centre points: the 0.01 s point at the +// geometric mean of the I²t window, the 0.1 s point at the geometric mean +// of the two 0.1 s gate currents, and √50 s at the geometric mean of +// Imin(10 s) and Imax(5 s). UMin = √(1.25 × 1.6) lies between Inf and If. +var gg = []struct { + ratingA float64 + imin10, imax5, imin01, imax1 float64 // A + i2tMin, i2tMax float64 // A²s at 0.01 s + re, tauE, tauB float64 // fitted +}{ + {16, 33, 65, 85, 150, 300, 1000, 0.042489, 0.040263, 15.5955}, + {20, 42, 85, 110, 200, 500, 1800, 0.038486, 0.040450, 16.8370}, + {25, 52, 110, 150, 260, 1000, 3000, 0.034171, 0.042170, 17.0309}, + {32, 75, 150, 200, 350, 1800, 5000, 0.030930, 0.040115, 21.9215}, + {35, 82.47, 164.94, 222.22, 387.17, 2209.86, 6331.18, 0.030262, 0.041014, 22.1359}, // interpolated + {40, 95, 190, 260, 450, 3000, 9000, 0.029308, 0.042411, 22.4587}, + {50, 125, 250, 350, 610, 5000, 16000, 0.024424, 0.038492, 25.1411}, + {63, 160, 320, 450, 820, 9000, 27000, 0.022440, 0.038871, 25.8209}, +} + +func ggBands(timeFactor float64, withGates bool) []Band { + uMin := math.Sqrt(1.25 * 1.6) + var out []Band + for _, r := range gg { + th := Thermal{UMin: uMin, Re: r.re, TauE: r.tauE, Rb: 1 - r.re - holderShare, TauB: r.tauB, + Rh: holderShare, TauH: holderTauS, TimeFactor: timeFactor} + b := Band{RatingA: r.ratingA, Thermal: th} + if withGates { + in := r.ratingA + b.Gates = []Gate{ + {1.25, 3600, true}, {1.6, 3600, false}, + {r.imin10 / in, 10, true}, {r.imax5 / in, 5, false}, + {r.imin01 / in, 0.1, true}, {r.imax1 / in, 0.1, false}, + } + b.I2tMin, b.I2tMax = r.i2tMin, r.i2tMax + b.Targets = []Point{ + {math.Sqrt(math.Sqrt(r.i2tMin*r.i2tMax)/0.01) / in, 0.01}, + {math.Sqrt(r.imin01*r.imax1) / in, 0.1}, + {math.Sqrt(r.imin10*r.imax5) / in, math.Sqrt(50)}, + } + } + out = append(out, b) + } + return out } +// amBand is aM at any rating. Gates from Mersen EduPack TM-104 (2012) slide +// 7, "as per IEC 60269", for any current rating: not melting within 60 s at +// 4 In, 0.5 s at 8 In and 0.2 s at 10 In; melting within 60 s at 6.3 In, +// 0.5 s at 12.5 In and 0.1 s at 19 In. Fitted like gG to √(4 × 6.3) In at +// 60 s, 10 In at √(0.2 × 0.5) s and 19 In at 0.05 s, with UMin = 4 (aM +// protects against short circuits, not overloads) and the holder at a 0.2 +// share, which the fit needs. +func amBand() Band { + return Band{RatingA: 0, + Thermal: Thermal{UMin: 4, Re: 0.2005, TauE: 0.2001, Rb: 0.5995, TauB: 45.9302, Rh: 0.2, TauH: holderTauS}, + Gates: []Gate{{4, 60, true}, {6.3, 60, false}, {8, 0.5, true}, {10, 0.2, true}, + {12.5, 0.5, false}, {19, 0.1, false}}, + Targets: []Point{{math.Sqrt(4 * 6.3), 60}, {10, math.Sqrt(0.2 * 0.5)}, {19, 0.05}}, + } +} + +var profiles = func() map[string]Profile { + m := map[string]Profile{ + Manual: {ID: Manual, Label: "Manual", Source: "never blows by itself"}, + IecGg: {ID: IecGg, Label: "IEC 60269 gG/gL (general purpose)", RefRatingA: 25, Bands: ggBands(1, true), + Source: "IEC 60269-1 gates, I²t and conventional currents per rating (see gg)"}, + IecAm: {ID: IecAm, Label: "IEC 60269 aM (motor, short-circuit)", Bands: []Band{amBand()}, + Source: "IEC 60269 aM gates (see amBand)"}, + // Not calibrated: gG with every time scaled, as the simulator's + // earlier curves had them against gG. + IecGm: {ID: IecGm, Label: "IEC 60269 gM (motor, dual-rated) - UK common", RefRatingA: 25, Bands: ggBands(2, false), + Source: "not calibrated: gG with times × 2 (a gM is a gG of its characteristic current, which the site does not give)"}, + IecGr: {ID: IecGr, Label: "IEC 60269 gR (semiconductor, fast)", RefRatingA: 25, Bands: ggBands(0.3, false), + Source: "not calibrated: gG with times × 0.3"}, + IecGs: {ID: IecGs, Label: "IEC 60269 gS (semiconductor, full-range)", RefRatingA: 25, Bands: ggBands(0.5, false), + Source: "not calibrated: gG with times × 0.5"}, + } + for id, p := range m { + for i := range p.Bands { + p.Bands[i].sys = newSystem(p.Bands[i].Thermal) + } + if len(p.Bands) > 0 { + b := p.band(p.RefRatingA) + for _, x := range []float64{1.5, 1.6, 2, 3, 5, 10, 20} { + if t := b.MeltSeconds(x, RefAmbientC); !math.IsInf(t, 1) { + p.Curve = append(p.Curve, Point{x, t}) + } + } + } + m[id] = p + } + return m +}() + // Lookup returns the profile with this id. func Lookup(id string) (Profile, bool) { p, ok := profiles[id] @@ -108,89 +238,261 @@ func All() []Profile { return out } -// TripSeconds is the time to blow at multiple × the rating, interpolated on -// log-log axes, with the last segment extended above its end. It is +Inf at -// or below 1 ×, and for a profile without a curve. -func (p Profile) TripSeconds(multiple float64) float64 { - if multiple <= 1 || len(p.Curve) == 0 { - return math.Inf(1) - } - pts := p.Curve - if multiple <= pts[0].Multiple { - return pts[0].Seconds +// band is the band nearest ratingA on log axes; a rating outside the bands +// takes the nearest end. +func (p Profile) band(ratingA float64) *Band { + if len(p.Bands) == 0 { + return nil } - for i := 1; i < len(pts); i++ { - if multiple <= pts[i].Multiple { - return logLog(multiple, pts[i-1], pts[i]) + best, dist := 0, math.Inf(1) + for i, b := range p.Bands { + if b.RatingA <= 0 || ratingA <= 0 { + return &p.Bands[i] + } + if d := math.Abs(math.Log(ratingA / b.RatingA)); d < dist { + best, dist = i, d } } - return logLog(multiple, pts[len(pts)-2], pts[len(pts)-1]) + return &p.Bands[best] } -func logLog(x float64, a, b Point) float64 { - if x <= 0 || a.Multiple <= 0 || b.Multiple <= 0 || a.Seconds <= 0 || b.Seconds <= 0 { - return b.Seconds +// Band returns the band a fuse rated ratingA uses, and false for a profile +// without bands. +func (p Profile) Band(ratingA float64) (Band, bool) { + b := p.band(ratingA) + if b == nil { + return Band{}, false } - lx, lx1, lx2 := math.Log10(x), math.Log10(a.Multiple), math.Log10(b.Multiple) - if lx2 == lx1 { - return b.Seconds + return *b, true +} + +// TripSeconds is the time to melt from cold at multiple × the rating, at the +// reference band and ambient; +Inf if it never melts there. +func (p Profile) TripSeconds(multiple float64) float64 { + return p.MeltSeconds(multiple, p.RefRatingA, RefAmbientC) +} + +// MeltSeconds is the time to melt from cold at multiple × a rating of +// ratingA in ambientC; +Inf if it never melts. +func (p Profile) MeltSeconds(multiple, ratingA, ambientC float64) float64 { + b := p.band(ratingA) + if b == nil { + return math.Inf(1) } - ly1, ly2 := math.Log10(a.Seconds), math.Log10(b.Seconds) - return math.Pow(10, ly1+(ly2-ly1)*(lx-lx1)/(lx2-lx1)) + return b.MeltSeconds(multiple, ambientC) +} + +// MeltSeconds is the band's time to melt from cold. +func (b *Band) MeltSeconds(multiple, ambientC float64) float64 { + return b.sys.timeToMelt([3]float64{}, multiple, b.sys.meltRise(ambientC)) } // Phase is the state of one phase's fuse. type Phase struct { - Damage float64 // 0 new, 1 blown + // Theta is each body's rise over ambient, in per-unit (element, body, + // holder); zero is cold. + Theta [3]float64 + // AmbientC is the cabinet's air; 0 means DefaultAmbientC. + AmbientC float64 + // Damage is the element's way from ambient to melting, 0 to 1; Step + // sets it. + Damage float64 Blown bool } // Status describes a phase after a step. type Status struct { - Ratio float64 `json:"ratio"` // current over the rating - TripSeconds float64 `json:"trip_s"` // time to blow at this ratio from new - RemainingSeconds float64 `json:"remaining_s"` - Damage float64 `json:"damage"` - Active bool `json:"active"` // building damage + Ratio float64 `json:"ratio"` // current over the rating + TripSeconds float64 `json:"trip_s"` // time to melt at this ratio from cold; 0 if never + RemainingSeconds float64 `json:"remaining_s"` // time to melt at this ratio from now; 0 if never + Damage float64 `json:"damage"` // the element's way from ambient to melting + Active bool `json:"active"` // it will melt if the current stays } // Step advances one phase by dt seconds at amps through a fuse rated -// ratingA. The sign of amps does not matter: export heats a fuse too. A -// profile without a curve leaves the phase as it is. +// ratingA. The sign of amps does not matter. The step is exact for a +// current held over dt, so the result does not depend on the step size. A +// profile without bands leaves the phase as it is. func (p Profile) Step(ph *Phase, amps, ratingA, dt float64) Status { - if len(p.Curve) == 0 { + b := p.band(ratingA) + if b == nil { return Status{Damage: ph.Damage} } if dt < 0 { dt = 0 } if ph.Blown { - ph.Damage = math.Max(ph.Damage, 1) - return Status{Damage: ph.Damage} + ph.Damage = 1 + return Status{Damage: 1} } - ratio := 0.0 + u := 0.0 if ratingA > 0 { - ratio = math.Abs(amps) / ratingA + u = math.Abs(amps) / ratingA + } + ambient := ph.AmbientC + if ambient == 0 { + ambient = DefaultAmbientC + } + melt := b.sys.meltRise(ambient) + st := Status{Ratio: u} + if cold := b.sys.timeToMelt([3]float64{}, u, melt); !math.IsInf(cold, 1) { + st.TripSeconds = cold + } + left := b.sys.timeToMelt(ph.Theta, u, melt) + if left <= dt { + ph.Theta = b.sys.advance(ph.Theta, u, left) + ph.Blown, ph.Damage = true, 1 + st.Damage = 1 + return st + } + ph.Theta = b.sys.advance(ph.Theta, u, dt) + ph.Damage = math.Max(0, math.Min(1, ph.Theta[0]/melt)) + st.Damage = ph.Damage + if !math.IsInf(left, 1) { + st.RemainingSeconds, st.Active = left-dt, true + } + return st +} + +// system is a band's ladder solved once: C dθ/dt = -G θ + (u², 0, 0). +// With y = θ/s (s = 1/√C), dy/dt = -M y + ..., M = s G s symmetric, so +// θ(t) = θ∞ + P e^{-Λt} Q (θ0 - θ∞), P = diag(s) V, Q = Vᵀ diag(1/s). +type system struct { + uMin2 float64 + inf [3]float64 // θ∞ at u = 1 + lambda [3]float64 + p, q [3][3]float64 +} + +func newSystem(t Thermal) system { + f := t.TimeFactor + if f <= 0 { + f = 1 + } + c := [3]float64{f * t.TauE / t.Re, f * t.TauB / t.Rb, f * t.TauH / t.Rh} + g1, g2, g3 := 1/t.Re, 1/t.Rb, 1/t.Rh + g := [3][3]float64{{g1, -g1, 0}, {-g1, g1 + g2, -g2}, {0, -g2, g2 + g3}} + var s [3]float64 + for i := range s { + s[i] = 1 / math.Sqrt(c[i]) + } + var m [3][3]float64 + for i := 0; i < 3; i++ { + for j := 0; j < 3; j++ { + m[i][j] = s[i] * g[i][j] * s[j] + } + } + lam, v := eigen3(m) + sys := system{uMin2: t.UMin * t.UMin, lambda: lam, + inf: [3]float64{t.Re + t.Rb + t.Rh, t.Rb + t.Rh, t.Rh}} + for i := 0; i < 3; i++ { + for j := 0; j < 3; j++ { + sys.p[i][j] = s[i] * v[i][j] + sys.q[i][j] = v[j][i] / s[j] + } + } + return sys +} + +// meltRise is the element's melting rise over ambient, in per-unit: UMin² +// at the reference ambient, less in a warmer cabinet. +func (s *system) meltRise(ambientC float64) float64 { + ambientC = math.Min(ambientC, MeltC-1) + return s.uMin2 * (MeltC - ambientC) / (MeltC - RefAmbientC) +} + +// modes returns θ∞ at u and the modal amplitudes of θ0 - θ∞. +func (s *system) modes(theta [3]float64, u float64) (inf, amp [3]float64) { + u2 := u * u + var d [3]float64 + for i := range inf { + inf[i] = s.inf[i] * u2 + d[i] = theta[i] - inf[i] + } + for k := 0; k < 3; k++ { + for j := 0; j < 3; j++ { + amp[k] += s.q[k][j] * d[j] + } } - start := p.TripStartRatio - if start <= 0 { - start = 1 + return inf, amp +} + +func (s *system) advance(theta [3]float64, u, dt float64) [3]float64 { + inf, amp := s.modes(theta, u) + var out [3]float64 + for i := 0; i < 3; i++ { + out[i] = inf[i] + for k := 0; k < 3; k++ { + out[i] += s.p[i][k] * math.Exp(-s.lambda[k]*dt) * amp[k] + } } - if ratio <= start { - if p.CooldownS > 0 { - ph.Damage = math.Max(0, ph.Damage-dt/p.CooldownS) + return out +} + +// timeToMelt is when the element first reaches melt at u from theta; +Inf +// if it never does (within ten million seconds). +func (s *system) timeToMelt(theta [3]float64, u, melt float64) float64 { + inf, amp := s.modes(theta, u) + f := func(t float64) float64 { + x := inf[0] - melt + for k := 0; k < 3; k++ { + x += s.p[0][k] * math.Exp(-s.lambda[k]*t) * amp[k] } - return Status{Ratio: ratio, Damage: ph.Damage} + return x + } + if f(0) >= 0 { + return 0 } - trip := p.TripSeconds(ratio) - if trip <= 0 || math.IsInf(trip, 1) { - return Status{Ratio: ratio, Damage: ph.Damage} + prev := 0.0 + for t := 1e-6; t <= 1e7; t *= 1.2 { + if f(t) >= 0 { + lo, hi := prev, t + for i := 0; i < 60 && hi-lo > 1e-9*hi; i++ { + mid := (lo + hi) / 2 + if f(mid) >= 0 { + hi = mid + } else { + lo = mid + } + } + return hi + } + prev = t } - trip = math.Max(trip, p.MinTripS) - ph.Damage = math.Min(1, ph.Damage+dt/trip) - if ph.Damage >= 1 { - ph.Blown = true - return Status{Ratio: ratio, TripSeconds: trip, Damage: ph.Damage} + return math.Inf(1) +} + +// eigen3 diagonalises a symmetric 3×3 matrix by Jacobi rotations: the +// eigenvalues and the eigenvectors as columns. +func eigen3(a [3][3]float64) ([3]float64, [3][3]float64) { + v := [3][3]float64{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}} + for sweep := 0; sweep < 50; sweep++ { + off := a[0][1]*a[0][1] + a[0][2]*a[0][2] + a[1][2]*a[1][2] + if off < 1e-30 { + break + } + for _, pq := range [][2]int{{0, 1}, {0, 2}, {1, 2}} { + p, q := pq[0], pq[1] + if math.Abs(a[p][q]) < 1e-300 { + continue + } + theta := (a[q][q] - a[p][p]) / (2 * a[p][q]) + t := math.Copysign(1, theta) / (math.Abs(theta) + math.Sqrt(theta*theta+1)) + c := 1 / math.Sqrt(t*t+1) + sn := t * c + for k := 0; k < 3; k++ { + akp, akq := a[k][p], a[k][q] + a[k][p], a[k][q] = c*akp-sn*akq, sn*akp+c*akq + } + for k := 0; k < 3; k++ { + apk, aqk := a[p][k], a[q][k] + a[p][k], a[q][k] = c*apk-sn*aqk, sn*apk+c*aqk + } + for k := 0; k < 3; k++ { + vkp, vkq := v[k][p], v[k][q] + v[k][p], v[k][q] = c*vkp-sn*vkq, sn*vkp+c*vkq + } + } } - return Status{Ratio: ratio, TripSeconds: trip, RemainingSeconds: (1 - ph.Damage) * trip, Damage: ph.Damage, Active: true} + return [3]float64{a[0][0], a[1][1], a[2][2]}, v } diff --git a/fuse/fuse_test.go b/fuse/fuse_test.go index 99b666e..d7545ac 100644 --- a/fuse/fuse_test.go +++ b/fuse/fuse_test.go @@ -5,112 +5,252 @@ import ( "testing" ) -func gg(t *testing.T) Profile { +func profile(t *testing.T, id string) Profile { t.Helper() - p, ok := Lookup(IecGg) + p, ok := Lookup(id) if !ok { - t.Fatal("no iec_gg profile") + t.Fatalf("no %s profile", id) } return p } -func TestTripSecondsFollowsTheCurve(t *testing.T) { - p := gg(t) - for _, pt := range p.Curve { - if got := p.TripSeconds(pt.Multiple); math.Abs(got-pt.Seconds) > 1e-9*pt.Seconds { - t.Errorf("%.2f×: %v s, want %v s", pt.Multiple, got, pt.Seconds) +// hold runs a phase at multiple × rating for seconds in steps of dt and +// says how long it took to blow (+Inf if it held). +func hold(p Profile, ph *Phase, multiple, rating, seconds, dt float64) float64 { + for t := 0.0; t < seconds; t += dt { + p.Step(ph, multiple*rating, rating, dt) + if ph.Blown { + return t + dt } } - if got := p.TripSeconds(1); !math.IsInf(got, 1) { - t.Errorf("1×: %v s, want never", got) - } - if got := p.TripSeconds(1.1); got != 3600 { - t.Errorf("1.1×: %v s, want the first point's 3600 s", got) + return math.Inf(1) +} + +// Every gG band, from cold at the reference ambient, keeps IEC 60269-1's +// gates and its pre-arcing I²t window at 0.01 s, and sits within 0.05 +// decades of the centre points it was fitted to. +func TestGGBandsMeetTheGates(t *testing.T) { + p := profile(t, IecGg) + for _, b := range p.Bands { + for _, g := range b.Gates { + got := b.MeltSeconds(g.Multiple, RefAmbientC) + if g.Min && got < g.Seconds || !g.Min && got > g.Seconds { + kind := "within" + if g.Min { + kind = "not sooner than" + } + t.Errorf("%g A at %.2f×: melts after %.3g s, want %s %g s", b.RatingA, g.Multiple, got, kind, g.Seconds) + } + } + u := multipleAt(&b, 0.01) + if i2t := math.Pow(u*b.RatingA, 2) * 0.01; i2t < b.I2tMin || i2t > b.I2tMax { + t.Errorf("%g A: pre-arcing I²t at 0.01 s %.0f A²s, want %g-%g", b.RatingA, i2t, b.I2tMin, b.I2tMax) + } + for _, pt := range b.Targets { + got := b.MeltSeconds(pt.Multiple, RefAmbientC) + if d := math.Abs(math.Log10(got / pt.Seconds)); d > 0.05 { + t.Errorf("%g A at %.2f×: %.3g s, %.2f decades from the fitted %.3g s", b.RatingA, pt.Multiple, got, d, pt.Seconds) + } + } } - if got := p.TripSeconds(1.4); got <= 1800 || got >= 3600 { - t.Errorf("1.4×: %.0f s, want between the 1.25× and 1.6× points", got) +} + +// multipleAt is the multiple at which the band melts after seconds. +func multipleAt(b *Band, seconds float64) float64 { + lo, hi := 1.0, 1000.0 + for i := 0; i < 100; i++ { + mid := math.Sqrt(lo * hi) + if b.MeltSeconds(mid, RefAmbientC) > seconds { + lo = mid + } else { + hi = mid + } } - if got := p.TripSeconds(40); got >= 0.5 { - t.Errorf("40×: %v s, want under the 20× point", got) + return hi +} + +// Through Step, as a site runs it: 1.25 × In holds for the conventional +// hour, 1.6 × In blows within it. +func TestConventionalCurrents(t *testing.T) { + p := profile(t, IecGg) + for _, b := range p.Bands { + cold := func() *Phase { return &Phase{AmbientC: RefAmbientC} } + if got := hold(p, cold(), 1.25, b.RatingA, 3600, 1); !math.IsInf(got, 1) { + t.Errorf("%g A: 1.25× blew after %.0f s, want it to hold 1 h", b.RatingA, got) + } + if got := hold(p, cold(), 1.6, b.RatingA, 3600, 1); got > 3600 { + t.Errorf("%g A: 1.6× held 1 h, want it to blow", b.RatingA) + } } } -// A phase held at a multiple of its rating blows after the curve's time, -// in one-second steps; the other phases are untouched. Export heats a fuse -// as much as import. -func TestStepBlowsOnTheCurve(t *testing.T) { - p := gg(t) - for _, tc := range []struct{ multiple, wantS float64 }{{2, 600}, {5, 20}, {10, 3}} { - var phases [3]Phase - blewAt := 0 - for s := 1; s <= 4000 && blewAt == 0; s++ { - p.Step(&phases[0], -tc.multiple*25, 25, 1) - p.Step(&phases[1], 0, 25, 1) - if phases[0].Blown { - blewAt = s +// A fuse that has carried its rating for an hour blows at 2 × in under +// half its cold time; at 1.5 × in under two thirds. +func TestPreloadBlowsSooner(t *testing.T) { + p := profile(t, IecGg) + for _, rating := range []float64{16, 25, 63} { + for _, c := range []struct{ multiple, factor float64 }{{2, 0.5}, {1.5, 2.0 / 3}} { + cold := hold(p, &Phase{AmbientC: RefAmbientC}, c.multiple, rating, 7200, 0.5) + warm := &Phase{AmbientC: RefAmbientC} + hold(p, warm, 1, rating, 3600, 1) + hot := hold(p, warm, c.multiple, rating, 7200, 0.5) + t.Logf("%g A, %g×: cold %.1f s, after 1 h at In %.1f s", rating, c.multiple, cold, hot) + if hot >= c.factor*cold { + t.Errorf("%g A, %g× after 1 h at In: %.1f s, want under %.0f %% of the cold %.1f s", rating, c.multiple, hot, 100*c.factor, cold) } } - if float64(blewAt) != tc.wantS || phases[1].Damage != 0 || phases[1].Blown { - t.Errorf("%.0f×: blew after %d s, want %v s", tc.multiple, blewAt, tc.wantS) + } +} + +// 30 s at 1.5 × every 2 min (0.5 × between) heats a fuse further each +// cycle: the element ends each burst hotter than the last, and after ten +// cycles a held 1.5 × blows sooner than from cold. +func TestShortOverloadsAccumulate(t *testing.T) { + p := profile(t, IecGg) + ph := &Phase{AmbientC: RefAmbientC} + var peaks []float64 + for cycle := 0; cycle < 10; cycle++ { + hold(p, ph, 1.5, 25, 30, 1) + peaks = append(peaks, ph.Damage) + hold(p, ph, 0.5, 25, 90, 1) + } + for i := 1; i < len(peaks); i++ { + if peaks[i] <= peaks[i-1] { + t.Fatalf("burst %d left damage %.3f, not above the last %.3f", i+1, peaks[i], peaks[i-1]) } } - var ph Phase - for s := 0; s < 4000; s++ { - p.Step(&ph, 25, 25, 1) + t.Logf("damage after each burst: %.3f", peaks) + after := hold(p, ph, 1.5, 25, 7200, 1) + cold := p.MeltSeconds(1.5, 25, RefAmbientC) + t.Logf("1.5× after ten bursts blows in %.0f s, from cold %.0f s", after, cold) + if after >= 0.9*cold { + t.Errorf("1.5× after ten bursts: %.0f s, want under 90 %% of the cold %.0f s", after, cold) + } +} + +// With no current a hot fuse cools along its slowest body's exponential: +// after an hour, each further 10 minutes take off the same share. +func TestCoolingIsExponential(t *testing.T) { + p := profile(t, IecGg) + ph := &Phase{AmbientC: RefAmbientC} + hold(p, ph, 1.3, 25, 7200, 1) + hold(p, ph, 0, 25, 3600, 1) + var ratios []float64 + for i := 0; i < 4; i++ { + before := ph.Theta[0] + hold(p, ph, 0, 25, 600, 1) + ratios = append(ratios, ph.Theta[0]/before) } - if ph.Blown || ph.Damage != 0 { - t.Errorf("at the rating: %+v, want no damage", ph) + for _, r := range ratios[1:] { + if math.Abs(r/ratios[0]-1) > 0.01 { + t.Fatalf("10-minute cooling ratios %.4f, want one constant share", ratios) + } + } + if want := math.Exp(-600 / holderTauS); math.Abs(ratios[0]-want) > 0.05 { + t.Errorf("cooling ratio %.3f per 10 min, want about e^(-600/τ) = %.3f", ratios[0], want) } } -func TestDamageRecoversBelowTheRating(t *testing.T) { - p := gg(t) - ph := Phase{Damage: 0.5} - st := p.Step(&ph, 10, 25, 75) - if math.Abs(ph.Damage-0.25) > 1e-12 || st.Active { - t.Fatalf("after 75 s at 0.4×: damage %v, want 0.25 (cooldown %v s)", ph.Damage, p.CooldownS) +func TestExportHeatsAsImport(t *testing.T) { + p := profile(t, IecGg) + var in, out Phase + for s := 0; s < 120; s++ { + a := p.Step(&in, 40, 25, 1) + b := p.Step(&out, -40, 25, 1) + if a != b || in != out { + t.Fatalf("second %d: import %+v, export %+v", s, a, b) + } } - st = p.Step(&ph, 50, 25, 1) - if !st.Active || math.Abs(st.TripSeconds-600) > 1e-9 || math.Abs(st.RemainingSeconds-(1-ph.Damage)*st.TripSeconds) > 1e-9 { - t.Fatalf("at 2×: %+v", st) +} + +// The step is exact: 0.5 s and 2 s steps blow at the same moment, to +// within the step. +func TestStepSizeDoesNotMatter(t *testing.T) { + p := profile(t, IecGg) + for _, m := range []float64{1.7, 2.5, 6} { + fine := hold(p, &Phase{AmbientC: RefAmbientC}, m, 25, 7200, 0.5) + coarse := hold(p, &Phase{AmbientC: RefAmbientC}, m, 25, 7200, 2) + cold := p.MeltSeconds(m, 25, RefAmbientC) + if math.Abs(fine-cold) > 0.5 || math.Abs(coarse-cold) > 2 { + t.Errorf("%g×: blew after %.1f s (0.5 s steps) and %.1f s (2 s steps), want %.2f s", m, fine, coarse, cold) + } } } -func TestBlownStaysBlownAndManualNeverBlows(t *testing.T) { - p := gg(t) - ph := Phase{Blown: true, Damage: 0.9} - if st := p.Step(&ph, 0, 25, 1000); !ph.Blown || ph.Damage != 1 || st.Damage != 1 { - t.Fatalf("blown phase after a cool second: %+v", ph) +// Status: trip_s is the cold time, remaining_s the time from now, damage +// the element's way from ambient to melting. +func TestStatusDescribesThePhase(t *testing.T) { + p := profile(t, IecGg) + ph := &Phase{AmbientC: RefAmbientC} + first := p.Step(ph, 50, 25, 1) + second := p.Step(ph, 50, 25, 1) + cold := p.MeltSeconds(2, 25, RefAmbientC) + if !first.Active || math.Abs(first.TripSeconds-cold) > 1e-6*cold || math.Abs(first.RemainingSeconds-(cold-1)) > 1e-3 { + t.Fatalf("first second at 2×: %+v, want active, trip %.2f s, remaining %.2f s", first, cold, cold-1) + } + if math.Abs(second.RemainingSeconds-(first.RemainingSeconds-1)) > 1e-3 || second.Damage <= first.Damage { + t.Fatalf("second second: %+v after %+v", second, first) + } + idle := p.Step(&Phase{AmbientC: RefAmbientC}, 25, 25, 1) + if idle.Active || idle.TripSeconds != 0 || idle.RemainingSeconds != 0 || idle.Damage <= 0 { + t.Fatalf("at the rating: %+v, want warming but never melting", idle) } - m, _ := Lookup(Manual) - ph = Phase{} - m.Step(&ph, 1000, 25, 1000) - if ph.Blown || ph.Damage != 0 { - t.Fatalf("manual: %+v, want untouched", ph) + if hold(p, ph, 2, 25, 7200, 1); !ph.Blown || ph.Damage != 1 { + t.Fatalf("held at 2×: %+v, want blown", ph) + } + if st := p.Step(ph, 0, 25, 1); st.Damage != 1 || !ph.Blown { + t.Fatalf("a blown fuse: %+v", st) } } -func TestAMBuildsOnlyFromItsStartRatio(t *testing.T) { - p, _ := Lookup(IecAm) - ph := Phase{} - p.Step(&ph, 3.9*25, 25, 100) - if ph.Damage != 0 { - t.Fatalf("aM at 3.9×: damage %v, want none", ph.Damage) +// A warmer cabinet leaves less margin to the melting point. +func TestWarmCabinetBlowsSooner(t *testing.T) { + p := profile(t, IecGg) + at20 := p.MeltSeconds(1.5, 25, 20) + at40 := p.MeltSeconds(1.5, 25, 40) + if at40 >= at20 { + t.Errorf("1.5× at 40 °C %.0f s, at 20 °C %.0f s: want sooner when warmer", at40, at20) + } + // IEC 60269's curves hold approximately at 30 °C as at 20 °C. + if at30 := p.MeltSeconds(5, 25, 30); math.Abs(math.Log10(at30/p.MeltSeconds(5, 25, 20))) > 0.1 { + t.Errorf("5× at 30 °C %.2f s, more than 0.1 decades from 20 °C", at30) } - p.Step(&ph, 5*25, 25, 1) - if ph.Damage == 0 { - t.Fatal("aM at 5×: no damage") +} + +// aM keeps its gates at any rating. +func TestAMGates(t *testing.T) { + p := profile(t, IecAm) + for _, rating := range []float64{16, 63} { + b, _ := p.Band(rating) + for _, g := range b.Gates { + got := b.MeltSeconds(g.Multiple, RefAmbientC) + if g.Min && got < g.Seconds || !g.Min && got > g.Seconds { + t.Errorf("aM at %g×: %.3g s, gate %g s (min %v)", g.Multiple, got, g.Seconds, g.Min) + } + } } } -func TestAllIsSortedAndComplete(t *testing.T) { - all := All() - if len(all) != 6 { - t.Fatalf("%d profiles", len(all)) +func TestManualNeverBlows(t *testing.T) { + p := profile(t, Manual) + ph := &Phase{} + if hold(p, ph, 20, 25, 600, 1); ph.Blown { + t.Fatal("a manual fuse blew") } - for i := 1; i < len(all); i++ { - if all[i-1].ID >= all[i].ID { - t.Fatalf("not sorted: %s before %s", all[i-1].ID, all[i].ID) +} + +// Every profile but manual melts above its curve's start and holds at its +// rating. +func TestEveryProfileHasACurve(t *testing.T) { + for _, p := range All() { + if p.ID == Manual { + continue + } + if len(p.Curve) == 0 || len(p.Bands) == 0 || p.Source == "" { + t.Errorf("%s: %d curve points, %d bands, source %q", p.ID, len(p.Curve), len(p.Bands), p.Source) + } + if !math.IsInf(p.MeltSeconds(1, 25, RefAmbientC), 1) { + t.Errorf("%s melts at its rating", p.ID) } } } From 5710bba787d01b0a7eafa9a2f6f10079cb6caa33 Mon Sep 17 00:00:00 2001 From: Fredrik Ahlgren Date: Fri, 23 Oct 2026 06:02:12 +0200 Subject: [PATCH 3/4] fix(fuse): follow ABB's gG curves within IEC's tube; explicit ambient; refuse bad input MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The review of a6baf29 (srcfl/device-simulator#9): - The 1.4-2 × region, where home overloads land, rested on three free choices (UMin, the holder's share and time constant) that no source pinned. All six thermal parameters are now fitted per band to ABB's pre-arcing curves for OFAF gG size 000 (1SCC317001C0201 p. 14, digitized from the PDF's vector paths), keeping every IEC 60269-1 gate, I²t window and conventional current: rms 0.08-0.095 decades, at most 0.17. A temperature-dependent element resistance made the fit worse and is left out. - Ambient 0 meant unset, so a 0 °C cabinet became 25 °C. Phase.AmbientC is now a pointer; nil is the default. - A NaN current, step or ambient, an infinite step, an ambient at the melting point or a rating without a band left NaN in the state or blew a fuse that never melts. Step now refuses them, leaves the state and says why in Status.Invalid; an infinite current blows the fuse. - 10 A has other conventional currents (1.5/1.9 In) than the 16-63 A band it borrowed. gG now models 16-63 A only; the simulator refuses other sizes, and a saved 10 A site reports the refusal. - Preheat puts a phase at the steady state of a held current, for scenarios that start mid-evening. - dt = 0 is documented as a time-to-melt probe; the step-size test checks the exactness to 1e-9; the Tau wording and the modes are corrected. Co-Authored-By: Claude Opus 5.5 --- CLAUDE.md | 2 +- cmd/simulator/app.go | 6 + cmd/simulator/fuse.go | 22 +++ cmd/simulator/fuse_test.go | 22 +++ cmd/simulator/main.go | 8 +- cmd/simulator/templates/index.html | 1 - fuse/fuse.go | 295 +++++++++++++++++++---------- fuse/fuse_test.go | 178 +++++++++++++---- internal/mcp/server.go | 2 +- 9 files changed, 390 insertions(+), 146 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index aab2f81..0226a6b 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -87,7 +87,7 @@ The app supports two modes controlled by the `--desktop` flag: - `cmd/simulator/wails.json` - Wails configuration **Public package:** -- `fuse/` - IEC 60269 fuse curves and per-phase damage; no dependencies, so an EMS's own simulation (FTW's `go/internal/sitesim`) imports the same model +- `fuse/` - each phase's fuse as a thermal model (element, body, holder; I² heating), within IEC 60269-1's gates and following ABB's gG curves for 16-63 A; no dependencies, so an EMS's own simulation (FTW's `go/internal/sitesim`) imports the same model **Internal packages:** - `internal/modbus/server.go` - Modbus TCP server with register storage and logging diff --git a/cmd/simulator/app.go b/cmd/simulator/app.go index 62383c7..139585f 100644 --- a/cmd/simulator/app.go +++ b/cmd/simulator/app.go @@ -908,6 +908,12 @@ func (a *App) UpdateSiteConfig(id int, config map[string]interface{}) error { site.BaseLoadW = baseLoad } if fuseSize, ok := config["fuse_size"].(float64); ok { + if fuseSize != math.Trunc(fuseSize) { + return fmt.Errorf("fuse size %v A is not whole amperes", fuseSize) + } + if err := checkFuseSize(int(fuseSize)); err != nil { + return err + } site.FuseSize = int(fuseSize) site.FuseDamage = PhaseLoads{} site.FuseHeat = [3]fuse.Phase{} diff --git a/cmd/simulator/fuse.go b/cmd/simulator/fuse.go index 9bda46e..89d6cb0 100644 --- a/cmd/simulator/fuse.go +++ b/cmd/simulator/fuse.go @@ -1,6 +1,7 @@ package main import ( + "fmt" "math" "time" @@ -34,6 +35,27 @@ var fuseProfiles = func() map[string]fuse.Profile { return m }() +// fuseSizes are the ratings the gG model is calibrated for (IEC 60269-1's +// 16-63 A band, plus 35 A); smaller and larger fuses have other +// conventional currents, which the model does not have. +var fuseSizes = func() []int { + gg, _ := fuse.Lookup(fuse.IecGg) + var out []int + for _, r := range gg.Ratings() { + out = append(out, int(r)) + } + return out +}() + +func checkFuseSize(a int) error { + for _, s := range fuseSizes { + if a == s { + return nil + } + } + return fmt.Errorf("fuse size must be one of %v A", fuseSizes) +} + func resolveFuseProfile(id string) fuse.Profile { if profile, ok := fuseProfiles[id]; ok { return profile diff --git a/cmd/simulator/fuse_test.go b/cmd/simulator/fuse_test.go index 9cb3784..c0105af 100644 --- a/cmd/simulator/fuse_test.go +++ b/cmd/simulator/fuse_test.go @@ -57,3 +57,25 @@ func TestSiteFuseRemembersItsHeat(t *testing.T) { t.Fatalf("a cleared fuse blew after %d s, want the cold %d s", again, cold) } } + +// A site saved with a size the model does not cover says so on every phase +// and never blows; new sizes outside the bands are refused. +func TestUnsupportedFuseSize(t *testing.T) { + site := &Site{FuseSize: 10, FuseProfile: FuseProfileIecGg} + for s := 0; s < 60; s++ { + updateFuseState(site, PhasePowers{L1: 40 * 230}, time.Second) + } + if site.FuseBlown.L1 || site.FuseOverload.L1.Invalid == "" { + t.Fatalf("10 A at 40 A: blown %v, status %+v; want refused", site.FuseBlown.L1, site.FuseOverload.L1) + } + for _, a := range []int{16, 20, 25, 32, 35, 40, 50, 63} { + if err := checkFuseSize(a); err != nil { + t.Errorf("%d A refused: %v", a, err) + } + } + for _, a := range []int{0, 10, 30, 80} { + if checkFuseSize(a) == nil { + t.Errorf("%d A accepted", a) + } + } +} diff --git a/cmd/simulator/main.go b/cmd/simulator/main.go index 1fb4932..721b182 100644 --- a/cmd/simulator/main.go +++ b/cmd/simulator/main.go @@ -145,7 +145,7 @@ type Site struct { FuseDamage PhaseLoads `json:"-"` FuseOverload FuseOverload `json:"-"` FuseHeat [3]fuse.Phase `json:"-"` // each phase's thermal state (L1, L2, L3) - FuseAmbientC float64 `json:"-"` // the fuse cabinet's air; 0 is 25 °C + FuseAmbientC *float64 `json:"-"` // the fuse cabinet's air; nil is 25 °C BoundMeter bool `json:"bound_meter"` // True if meter was auto-created with site P1Buffer *P1IntervalBuffer // P1 port interval buffer @@ -2644,11 +2644,9 @@ func handleSiteConfig(w http.ResponseWriter, r *http.Request, site *Site) { site.BaseLoadW = *config.BaseLoadW } if config.FuseSize != nil { - // Validate fuse size (common values: 10, 16, 20, 25, 32 Amps) - validFuses := map[int]bool{10: true, 16: true, 20: true, 25: true, 32: true, 35: true, 40: true, 50: true, 63: true} - if !validFuses[*config.FuseSize] { + if err := checkFuseSize(*config.FuseSize); err != nil { site.mu.Unlock() - http.Error(w, "Fuse size must be 10, 16, 20, 25, 32, 35, 40, 50 or 63 Amps", http.StatusBadRequest) + http.Error(w, err.Error(), http.StatusBadRequest) return } site.FuseSize = *config.FuseSize diff --git a/cmd/simulator/templates/index.html b/cmd/simulator/templates/index.html index ae568e9..c3b20e6 100644 --- a/cmd/simulator/templates/index.html +++ b/cmd/simulator/templates/index.html @@ -317,7 +317,6 @@

Site Settings

- - - - - - - - +
diff --git a/fuse/abb_test.go b/fuse/abb_test.go new file mode 100644 index 0000000..9097e23 --- /dev/null +++ b/fuse/abb_test.go @@ -0,0 +1,51 @@ +package fuse + +import ( + "math" + "testing" +) + +// Every gG band against every vertex of the ABB curve it was fitted to, +// not only the eleven fitted points: a three-body ladder cannot follow the +// curve's whole shape (see gg), so the bound is loose, 0.5 decades, and +// only catches a regression. +func TestGGBandsAgainstTheWholeABBCurve(t *testing.T) { + p := profile(t, IecGg) + for _, b := range p.Bands { + vs := abbVertices[b.RatingA] + if len(vs) != 48 { + t.Fatalf("%g A: %d ABB vertices", b.RatingA, len(vs)) + } + var sq, worst, worstAt float64 + for _, v := range vs { + got := b.MeltSeconds(v[0], RefAmbientC) + d := math.Log10(got / v[1]) + if math.IsInf(got, 1) { + d = math.Inf(1) + } + sq += d * d + if math.Abs(d) > math.Abs(worst) { + worst, worstAt = d, v[0] + } + } + rms := math.Sqrt(sq / float64(len(vs))) + t.Logf("%g A: rms %.3f decades, worst %+.2f at %.2f ×", b.RatingA, rms, worst, worstAt) + if math.Abs(worst) > 0.5 { + t.Errorf("%g A: %+.2f decades from ABB at %.2f ×, want within 0.5", b.RatingA, worst, worstAt) + } + } +} + +// abbVertices are the 48 vertices of each ABB OFAF 000 gG pre-arcing curve +// (1SCC317001C0201 page 14), as multiples of In and seconds, read from the +// PDF's vector paths. +var abbVertices = map[float64][][2]float64{ + 16: {{1.4324, 7990.5}, {1.4564, 5987.4}, {1.4718, 5000}, {1.4898, 3994.1}, {1.5147, 2993.7}, {1.5517, 1997}, {1.5777, 1496.8}, {1.6139, 998.51}, {1.6334, 801.12}, {1.6608, 600.3}, {1.681, 499.25}, {1.7039, 400.56}, {1.7351, 300.15}, {1.7857, 200.28}, {1.8239, 150.07}, {1.8855, 100.14}, {1.9228, 80.014}, {1.9758, 59.956}, {2.0118, 50.068}, {2.0612, 39.996}, {2.1308, 29.978}, {2.243, 20.003}, {2.3329, 14.989}, {2.478, 9.9988}, {2.5657, 7.9895}, {2.6886, 6.0112}, {2.7752, 4.9993}, {2.8909, 3.9947}, {3.0522, 3.0056}, {3.3065, 1.9973}, {3.5125, 1.5028}, {3.8396, 0.99866}, {4.0422, 0.80124}, {4.3263, 0.60039}, {4.5201, 0.49932}, {4.7722, 0.40051}, {5.1313, 0.30019}, {5.6945, 0.20031}, {6.1324, 0.15009}, {6.8258, 0.10012}, {7.251, 0.080004}, {7.8318, 0.059965}, {8.2196, 0.050076}, {8.7308, 0.040002}, {9.4583, 0.029982}, {10.593, 0.020001}, {11.511, 0.014987}, {13.008, 0.01}}, + 20: {{1.4136, 7990.5}, {1.435, 5987.4}, {1.4481, 5000}, {1.4657, 3994.1}, {1.488, 2993.7}, {1.5221, 1997}, {1.5452, 1496.8}, {1.5807, 998.51}, {1.6022, 801.12}, {1.6316, 600.3}, {1.6513, 499.25}, {1.6766, 400.56}, {1.7124, 300.15}, {1.7675, 200.28}, {1.8107, 150.07}, {1.8833, 100.14}, {1.9264, 80.014}, {1.9914, 59.956}, {2.0341, 50.068}, {2.0932, 39.996}, {2.1771, 29.978}, {2.3092, 20.003}, {2.4198, 14.989}, {2.5899, 9.9988}, {2.6978, 7.9895}, {2.8442, 6.0112}, {2.9491, 4.9993}, {3.0812, 3.9947}, {3.2681, 3.0056}, {3.5564, 1.9973}, {3.7893, 1.5028}, {4.1486, 0.99866}, {4.367, 0.80124}, {4.6672, 0.60039}, {4.8763, 0.49932}, {5.1409, 0.40051}, {5.5192, 0.30019}, {6.1062, 0.20031}, {6.5656, 0.15009}, {7.2975, 0.10012}, {7.7521, 0.080004}, {8.3722, 0.059965}, {8.8003, 0.050076}, {9.362, 0.040002}, {10.173, 0.029982}, {11.463, 0.020001}, {12.511, 0.014987}, {14.225, 0.01}}, + 25: {{1.397, 7990.5}, {1.4204, 5987.4}, {1.4355, 5000}, {1.4551, 3994.1}, {1.4795, 2993.7}, {1.5134, 1997}, {1.5412, 1496.8}, {1.5788, 998.51}, {1.6028, 801.12}, {1.6321, 600.3}, {1.6544, 499.25}, {1.6797, 400.56}, {1.7156, 300.15}, {1.7734, 200.28}, {1.8195, 150.07}, {1.8895, 100.14}, {1.9358, 80.014}, {1.9951, 59.956}, {2.0377, 50.068}, {2.0971, 39.996}, {2.1745, 29.978}, {2.3028, 20.003}, {2.406, 14.989}, {2.5674, 9.9988}, {2.6703, 7.9895}, {2.8109, 6.0112}, {2.9059, 4.9993}, {3.0314, 3.9947}, {3.2104, 3.0056}, {3.499, 1.9973}, {3.7223, 1.5028}, {4.0879, 0.99866}, {4.3164, 0.80124}, {4.6335, 0.60039}, {4.8485, 0.49932}, {5.1426, 0.40051}, {5.554, 0.30019}, {6.2112, 0.20031}, {6.7488, 0.15009}, {7.6035, 0.10012}, {8.1387, 0.080004}, {8.9105, 0.059965}, {9.4511, 0.050076}, {10.147, 0.040002}, {11.159, 0.029982}, {12.803, 0.020001}, {14.124, 0.014987}, {16.277, 0.01}}, + 32: {{1.3939, 7990.5}, {1.4216, 5987.4}, {1.4367, 5000}, {1.4586, 3994.1}, {1.4851, 2993.7}, {1.5261, 1997}, {1.554, 1496.8}, {1.5969, 998.51}, {1.6211, 801.12}, {1.6558, 600.3}, {1.6784, 499.25}, {1.7065, 400.56}, {1.7456, 300.15}, {1.8074, 200.28}, {1.857, 150.07}, {1.9315, 100.14}, {1.9786, 80.014}, {2.0455, 59.956}, {2.0892, 50.068}, {2.1501, 39.996}, {2.2328, 29.978}, {2.3683, 20.003}, {2.4744, 14.989}, {2.6483, 9.9988}, {2.7544, 7.9895}, {2.9039, 6.0112}, {3.0064, 4.9993}, {3.141, 3.9947}, {3.3316, 3.0056}, {3.6363, 1.9973}, {3.8748, 1.5028}, {4.2614, 0.99866}, {4.4996, 0.80124}, {4.8376, 0.60039}, {5.0694, 0.49932}, {5.3769, 0.40051}, {5.8166, 0.30019}, {6.5236, 0.20031}, {7.1101, 0.15009}, {8.0599, 0.10012}, {8.6521, 0.080004}, {9.5155, 0.059965}, {10.124, 0.050076}, {10.934, 0.040002}, {12.116, 0.029982}, {14.092, 0.020001}, {15.734, 0.014987}, {18.493, 0.01}}, + 35: {{1.3931, 7990.5}, {1.4208, 5987.4}, {1.4381, 5000}, {1.46, 3994.1}, {1.4889, 2993.7}, {1.5299, 1997}, {1.5602, 1496.8}, {1.6057, 998.51}, {1.6301, 801.12}, {1.6649, 600.3}, {1.6876, 499.25}, {1.7186, 400.56}, {1.758, 300.15}, {1.82, 200.28}, {1.8675, 150.07}, {1.9451, 100.14}, {1.9897, 80.014}, {2.057, 59.956}, {2.1007, 50.068}, {2.162, 39.996}, {2.2452, 29.978}, {2.3814, 20.003}, {2.488, 14.989}, {2.6629, 9.9988}, {2.7696, 7.9895}, {2.9199, 6.0112}, {3.023, 4.9993}, {3.1584, 3.9947}, {3.35, 3.0056}, {3.6568, 1.9973}, {3.8961, 1.5028}, {4.2849, 0.99866}, {4.5244, 0.80124}, {4.8643, 0.60039}, {5.0973, 0.49932}, {5.4066, 0.40051}, {5.8487, 0.30019}, {6.5697, 0.20031}, {7.1493, 0.15009}, {8.1044, 0.10012}, {8.6999, 0.080004}, {9.568, 0.059965}, {10.18, 0.050076}, {10.995, 0.040002}, {12.183, 0.029982}, {14.17, 0.020001}, {15.82, 0.014987}, {18.595, 0.01}}, + 40: {{1.3987, 7990.5}, {1.4263, 5987.4}, {1.4437, 5000}, {1.4634, 3994.1}, {1.4925, 2993.7}, {1.5335, 1997}, {1.564, 1496.8}, {1.6071, 998.51}, {1.6315, 801.12}, {1.6664, 600.3}, {1.6891, 499.25}, {1.7201, 400.56}, {1.7595, 300.15}, {1.8244, 200.28}, {1.8747, 150.07}, {1.9527, 100.14}, {2.0005, 80.014}, {2.0681, 59.956}, {2.1154, 50.068}, {2.177, 39.996}, {2.264, 29.978}, {2.4051, 20.003}, {2.5203, 14.989}, {2.7016, 9.9988}, {2.8142, 7.9895}, {2.9715, 6.0112}, {3.0764, 4.9993}, {3.2191, 3.9947}, {3.4193, 3.0056}, {3.7382, 1.9973}, {3.9948, 1.5028}, {4.3934, 0.99866}, {4.6461, 0.80124}, {4.9951, 0.60039}, {5.2425, 0.49932}, {5.5605, 0.40051}, {6.0238, 0.30019}, {6.7665, 0.20031}, {7.3642, 0.15009}, {8.3471, 0.10012}, {8.9752, 0.080004}, {9.8708, 0.059965}, {10.485, 0.050076}, {11.342, 0.040002}, {12.549, 0.029982}, {14.573, 0.020001}, {16.246, 0.014987}, {19.039, 0.01}}, + 50: {{1.4162, 7990.5}, {1.4464, 5987.4}, {1.464, 5000}, {1.4884, 3994.1}, {1.518, 2993.7}, {1.5645, 1997}, {1.598, 1496.8}, {1.6469, 998.51}, {1.6745, 801.12}, {1.7155, 600.3}, {1.7416, 499.25}, {1.7761, 400.56}, {1.8224, 300.15}, {1.8953, 200.28}, {1.9533, 150.07}, {2.0439, 100.14}, {2.0971, 80.014}, {2.1744, 59.956}, {2.2277, 50.068}, {2.2995, 39.996}, {2.3988, 29.978}, {2.5597, 20.003}, {2.6865, 14.989}, {2.8928, 9.9988}, {3.0177, 7.9895}, {3.1958, 6.0112}, {3.3137, 4.9993}, {3.4728, 3.9947}, {3.7001, 3.0056}, {4.051, 1.9973}, {4.3291, 1.5028}, {4.7758, 0.99866}, {5.0427, 0.80124}, {5.4215, 0.60039}, {5.6812, 0.49932}, {6.0259, 0.40051}, {6.5086, 0.30019}, {7.2892, 0.20031}, {7.908, 0.15009}, {8.9377, 0.10012}, {9.5659, 0.080004}, {10.489, 0.059965}, {11.142, 0.050076}, {12.016, 0.040002}, {13.274, 0.029982}, {15.369, 0.020001}, {17.109, 0.014987}, {20.017, 0.01}}, + 63: {{1.317, 7990.5}, {1.3512, 5987.4}, {1.3718, 5000}, {1.3991, 3994.1}, {1.4354, 2993.7}, {1.4861, 1997}, {1.5247, 1496.8}, {1.5858, 998.51}, {1.6197, 801.12}, {1.6695, 600.3}, {1.7025, 499.25}, {1.7467, 400.56}, {1.8058, 300.15}, {1.8979, 200.28}, {1.971, 150.07}, {2.0811, 100.14}, {2.1482, 80.014}, {2.2409, 59.956}, {2.3027, 50.068}, {2.3805, 39.996}, {2.491, 29.978}, {2.662, 20.003}, {2.7938, 14.989}, {2.9993, 9.9988}, {3.124, 7.9895}, {3.2936, 6.0112}, {3.4049, 4.9993}, {3.5519, 3.9947}, {3.7562, 3.0056}, {4.075, 1.9973}, {4.3223, 1.5028}, {4.7108, 0.99866}, {4.9517, 0.80124}, {5.2759, 0.60039}, {5.5037, 0.49932}, {5.8024, 0.40051}, {6.2293, 0.30019}, {6.9131, 0.20031}, {7.4784, 0.15009}, {8.3874, 0.10012}, {8.964, 0.080004}, {9.799, 0.059965}, {10.378, 0.050076}, {11.158, 0.040002}, {12.291, 0.029982}, {14.165, 0.020001}, {15.719, 0.014987}, {18.282, 0.01}}, +} diff --git a/fuse/fuse.go b/fuse/fuse.go index ca65beb..4fd40ab 100644 --- a/fuse/fuse.go +++ b/fuse/fuse.go @@ -74,7 +74,7 @@ type Gate struct { // holder-to-ambient (they add up to 1). Each Tau is that body's heat // capacity times its own outgoing resistance; they are parameters, not the // system's modes, which come out of the coupled ladder (at 25 A 0.12 s, -// 9.4 s and 1,315 s). +// 9.4 s and 1,308 s). type Thermal struct { UMin float64 // steady melting multiple at the reference ambient Re, TauE float64 // element @@ -129,7 +129,18 @@ type Profile struct { // parameters are fitted, least squares on log time from cold at 20 °C, to // ABB's melting times at 1.5, 1.6, 2, 2.5, 3, 4, 5, 8 and 10 × In, the // current that melts in 0.01 s and the curve's start at 7,990 s; the fit -// keeps every gate. Residuals: rms 0.08-0.095 decades, at most 0.17. +// keeps every gate. At those eleven points the residuals are rms +// 0.08-0.095 decades, at most 0.17. Against all 48 vertices of each ABB +// curve they are rms 0.13-0.20, at most 0.22-0.46: the model is too slow +// at 1.6-1.95 × (16 A at 1.82 ×: ABB 150 s, model 447 s) and too fast at +// 2.0-2.5 × (down to -0.27). Where the element's quasi-steady rise over +// the body reaches the melting rise, u = UMin/√(Re + Rb) = 2.02-2.09 × In, +// the body stage takes over from the holder and the curve drops steeply: +// its log-log slope between 1.9 and 2.3 × reaches -18 to -52, ABB's -7 to +// -11. Fitting all 48 vertices reaches only rms 0.12-0.15, so it is the +// three-body ladder's shape. A simulation judging overloads at 1.9-2.1 × +// In sits on that step: a small change of current there moves the melting +// time a lot. // Letting the element's resistance rise with temperature (silver, 0.0038 // per K) made the fit worse (rms 0.11-0.14), so it is left out. var gg = []struct { @@ -141,22 +152,22 @@ var gg = []struct { abb [9]float64 abb001A, abbStartA float64 // A melting in 0.01 s; A at the curve's 7,990 s start }{ - {16, 33, 65, 85, 150, 300, 1000, 1.428717, 0.073252, 0.074250, 0.500994, 7.03215, 1616.13, - [9]float64{3546.8, 1165.2, 53.1, 9.445, 3.29, 0.8381, 0.3328, 0.0554, 0.0246}, 208.1, 22.9}, - {20, 42, 85, 110, 200, 500, 1800, 1.414377, 0.070387, 0.091390, 0.513906, 9.47446, 1241.79, - [9]float64{2592.9, 819.5, 57.8, 12.3, 4.58, 1.177, 0.4498, 0.0711, 0.0318}, 284.5, 28.3}, - {25, 52, 110, 150, 260, 1000, 3000, 1.395152, 0.073492, 0.128142, 0.528300, 8.82109, 1302.24, - [9]float64{2341.3, 821.6, 58.7, 11.8, 4.222, 1.098, 0.445, 0.0847, 0.0419}, 406.9, 34.9}, - {32, 75, 150, 200, 350, 1800, 5000, 1.391843, 0.075659, 0.167797, 0.532303, 10.60028, 1501.37, - [9]float64{2581.2, 970.5, 72.9, 14.1, 5.056, 1.311, 0.5272, 0.1026, 0.0519}, 591.8, 44.6}, - {35, 82.47, 164.94, 222.22, 387.17, 2209.86, 6331.18, 1.391721, 0.075520, 0.169985, 0.532484, 10.92192, 1591.95, // gates interpolated - [9]float64{2680.3, 1049.8, 76.5, 14.6, 5.206, 1.342, 0.5387, 0.1044, 0.0527}, 650.8, 48.8}, - {40, 95, 190, 260, 450, 3000, 9000, 1.394819, 0.073959, 0.177486, 0.530528, 11.82141, 1606.21, - [9]float64{2778.0, 1066.7, 80.2, 15.8, 5.714, 1.494, 0.5982, 0.1149, 0.0577}, 761.5, 55.9}, - {50, 125, 250, 350, 610, 5000, 16000, 1.409631, 0.071707, 0.197219, 0.525392, 16.77180, 1951.30, - [9]float64{3566.2, 1471.4, 121.6, 23.2, 8.242, 2.115, 0.8293, 0.1445, 0.0696}, 1000.8, 70.8}, - {63, 160, 320, 450, 820, 9000, 27000, 1.314812, 0.063486, 0.165089, 0.603945, 19.24608, 1551.42, - [9]float64{1798.4, 910.0, 134.6, 29.3, 9.987, 2.193, 0.7666, 0.1183, 0.0563}, 1151.8, 83.0}, + {16, 33, 65, 85, 150, 300, 1000, 1.429959, 0.073378, 0.074249, 0.500174, 7.03102, 1604.18, + [9]float64{3546.8, 1165.2, 53.1, 9.445, 3.29, 0.8381, 0.3328, 0.0554, 0.0246}, 208.1, 22.918}, + {20, 42, 85, 110, 200, 500, 1800, 1.412945, 0.070246, 0.091392, 0.514825, 9.47660, 1251.63, + [9]float64{2592.9, 819.5, 57.8, 12.3, 4.58, 1.177, 0.4498, 0.0711, 0.0318}, 284.5, 28.272}, + {25, 52, 110, 150, 260, 1000, 3000, 1.396178, 0.073599, 0.128139, 0.527645, 8.81978, 1295.49, + [9]float64{2341.3, 821.6, 58.7, 11.8, 4.222, 1.098, 0.445, 0.0847, 0.0419}, 406.9, 34.925}, + {32, 75, 150, 200, 350, 1800, 5000, 1.392001, 0.075676, 0.167797, 0.532203, 10.60002, 1500.19, + [9]float64{2581.2, 970.5, 72.9, 14.1, 5.056, 1.311, 0.5272, 0.1026, 0.0519}, 591.8, 44.605}, + {35, 82.47, 164.94, 222.22, 387.17, 2209.86, 6331.18, 1.390453, 0.075385, 0.169990, 0.533283, 10.92408, 1601.99, // gates interpolated + [9]float64{2680.3, 1049.8, 76.5, 14.6, 5.206, 1.342, 0.5387, 0.1044, 0.0527}, 650.8, 48.759}, + {40, 95, 190, 260, 450, 3000, 9000, 1.396104, 0.074093, 0.177481, 0.529719, 11.81891, 1595.82, + [9]float64{2778.0, 1066.7, 80.2, 15.8, 5.714, 1.494, 0.5982, 0.1149, 0.0577}, 761.5, 55.948}, + {50, 125, 250, 350, 610, 5000, 16000, 1.409864, 0.071730, 0.197218, 0.525250, 16.77090, 1948.86, + [9]float64{3566.2, 1471.4, 121.6, 23.2, 8.242, 2.115, 0.8293, 0.1445, 0.0696}, 1000.8, 70.810}, + {63, 160, 320, 450, 820, 9000, 27000, 1.314326, 0.063439, 0.165090, 0.604208, 19.24866, 1554.48, + [9]float64{1798.4, 910.0, 134.6, 29.3, 9.987, 2.193, 0.7666, 0.1183, 0.0563}, 1151.8, 82.971}, } var abbMultiples = [9]float64{1.5, 1.6, 2, 2.5, 3, 4, 5, 8, 10} @@ -396,9 +407,11 @@ func (p Profile) Step(ph *Phase, amps, ratingA, dt float64) Status { } u := math.Abs(amps) / ratingA st := Status{Ratio: u} - if math.IsInf(u, 1) { + if math.IsInf(u*u, 1) { + // A current whose square overflows (or an infinite one) melts the + // fuse at once; the ratio stays finite for JSON. ph.Blown, ph.Damage = true, 1 - st.Damage = 1 + st.Ratio, st.Damage = math.MaxFloat64, 1 return st } melt := b.sys.meltRise(ambient) @@ -423,13 +436,17 @@ func (p Profile) Step(ph *Phase, amps, ratingA, dt float64) Status { // Preheat puts the phase in the steady state of a current held for ever in // its ambient, as a fuse that has carried the evening's load: a scenario -// that starts mid-evening starts warm, not cold. It refuses a current that -// would melt the fuse before reaching that state, and anything Step -// refuses, leaving the phase as it was. +// that starts mid-evening starts warm, not cold. It refuses a blown fuse +// (clear Blown and Theta to replace one), a current that would melt the +// fuse before reaching that state, and anything Step refuses, leaving the +// phase as it was. func (p Profile) Preheat(ph *Phase, amps, ratingA float64) error { if len(p.Bands) == 0 { return fmt.Errorf("%s has no fuse model", p.ID) } + if ph.Blown { + return fmt.Errorf("the fuse is blown") + } if st := p.Step(&Phase{AmbientC: ph.AmbientC}, amps, ratingA, 0); st.Invalid != "" { return fmt.Errorf("%s", st.Invalid) } @@ -439,7 +456,7 @@ func (p Profile) Preheat(ph *Phase, amps, ratingA float64) error { if inf[0] >= b.sys.meltRise(ph.ambient()) { return fmt.Errorf("%.3g × %g A melts the fuse before it settles", u, ratingA) } - ph.Theta, ph.Blown = inf, false + ph.Theta = inf ph.Damage = inf[0] / b.sys.meltRise(ph.ambient()) return nil } diff --git a/fuse/fuse_test.go b/fuse/fuse_test.go index 9ec2195..8ed3439 100644 --- a/fuse/fuse_test.go +++ b/fuse/fuse_test.go @@ -1,6 +1,7 @@ package fuse import ( + "encoding/json" "math" "testing" ) @@ -129,6 +130,13 @@ func TestPreheat(t *testing.T) { t.Errorf("preheat at %v A: err %v, state %v; want refused, cold", amps, err, ph.Theta) } } + // A blown fuse stays blown: preheating does not replace it. + blown := ref() + hold(p, blown, 3, 25, 600, 1) + before := *blown + if err := p.Preheat(blown, 25, 25); err == nil || !blown.Blown || blown.Theta != before.Theta { + t.Errorf("preheat on a blown fuse: err %v, blown %v; want refused, still blown", err, blown.Blown) + } } // 30 s at 1.5 × every 2 min (0.5 × between) heats a fuse further each @@ -300,9 +308,18 @@ func TestInvalidInputLeavesTheState(t *testing.T) { t.Errorf("%s: %+v, state %v → %v; want refused, unchanged", c.name, st, before.Theta, ph.Theta) } } - ph := *warm - if st := p.Step(&ph, math.Inf(-1), 25, 1); !ph.Blown || st.Damage != 1 || st.Invalid != "" { - t.Errorf("an infinite current: %+v, want blown", st) + for _, amps := range []float64{math.Inf(-1), 1e160 * 25} { + ph := *warm + st := p.Step(&ph, amps, 25, 1) + _, err := json.Marshal(st) + if !ph.Blown || st.Damage != 1 || st.Invalid != "" || math.IsInf(st.Ratio, 0) || err != nil { + t.Errorf("%g A: %+v (json: %v), want blown with a finite ratio", amps, st, err) + } + for _, x := range ph.Theta { + if math.IsNaN(x) || math.IsInf(x, 0) { + t.Errorf("%g A left the state %v", amps, ph.Theta) + } + } } if p.Supports(10) || p.Supports(80) || !p.Supports(16) || !p.Supports(30) || !p.Supports(63) { t.Error("gG supports 16-63 A only")