Skip to content

Add Jiles-Atherton B-H hysteresis to inductor #344

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
esaruoho wants to merge 1 commit into
base: v3-dev
Choose a base branch
from
+167 −21
Open

Add Jiles-Atherton B-H hysteresis to inductor #344

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
112 changes: 106 additions & 6 deletions src/com/lushprojects/circuitjs1/client/Inductor.java
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
/*
/*
Copyright (C) Paul Falstad and Iain Sharp

This file is part of CircuitJS1.

CircuitJS1 is free software: you can redistribute it and/or modify
Expand Down Expand Up @@ -29,6 +29,20 @@ class Inductor {
double compResistance, current;
double curSourceValue;
double saturationCurrent; // 0 = disabled (linear), >0 = saturation onset current (A)

// Jiles-Atherton hysteresis (0 = disabled). All internal state uses normalized
// dimensionless m in [-1,1] scaled by coerciveCurrent. Physics params (shape an,
// coupling alphaN) have fixed internal defaults so user only sees Ic and c.
double coerciveCurrent; // k in JA, A. 0 = hysteresis disabled.
double reversibility; // c in JA, 0..1.
static final double JA_SHAPE = 1.0/3.0; // an, gives dm/dh=1 at h=0 (virgin L = L0)
static final double JA_COUPLING = 1e-3; // alpha, small = weak interdomain
double mIrr; // irreversible magnetization, dimensionless
double mTotal; // mIrr + c*(man-mIrr)
double prevI; // I from previous startIteration
double prevM; // mTotal from previous startIteration
double lastLeff; // effective inductance last computed

Inductor(SimulationManager s) {
sim = s;
nodes = new int[2];
Expand All @@ -37,28 +51,102 @@ void setup(double ic, double cr, int f) {
inductance = ic;
current = cr;
flags = f;
lastLeff = inductance;
}
void setup(double ic, double cr, int f, double isat) {
setup(ic, cr, f);
saturationCurrent = isat;
}
void setup(double ic, double cr, int f, double isat, double icoerc, double rev) {
setup(ic, cr, f, isat);
coerciveCurrent = icoerc;
reversibility = rev;
}
boolean isTrapezoidal() { return (flags & FLAG_BACK_EULER) == 0; }
boolean hasHysteresis() { return coerciveCurrent > 0; }
void reset() { resetTo(0); }
void resetTo(double c) {
// need to set curSourceValue here in case one of inductor nodes is node 0. In that case
// calculateCurrent() may get called (from setNodeVoltage()) when analyzing circuit, before
// startIteration() gets called
curSourceValue = current = c;
mIrr = mTotal = 0;
prevI = c;
prevM = 0;
lastLeff = inductance;
}

// compute effective inductance with saturation: L(I) = L0 / (1 + (I/Isat)^2)
// smooth rolloff: at |I|=Isat, L=L0/2; at |I|=3*Isat, L=L0/10
// When hysteresis is active, effective L is derived from JA state instead
// and this function returns the last computed value.
double calcEffectiveInductance(double i) {
if (hasHysteresis()) return lastLeff;
if (saturationCurrent <= 0) return inductance;
double ratio = i / saturationCurrent;
return inductance / (1 + ratio * ratio);
}

// Advance Jiles-Atherton state by one timestep. h = I/Ic (dimensionless
// field). Magnetization m evolves according to:
// Man(he) = coth(he/an) - an/he (Langevin, normalized Ms=1)
// he = h + alpha*m
// dMirr/dh = (Man-Mirr)/(sign*1 - alpha*(Man-Mirr)) (k normalized to 1)
// m = Mirr + c*(Man-Mirr)
void advanceHysteresis(double i) {
double h = i / coerciveCurrent;
double dh = h - prevI / coerciveCurrent;
double sign = (dh >= 0) ? 1.0 : -1.0;
double he = h + JA_COUPLING * mTotal;
double man;
double heOverA = he / JA_SHAPE;
if (Math.abs(he) < 1e-6) {
man = he / (3.0 * JA_SHAPE);
} else {
// coth(x) = (e^x + e^-x) / (e^x - e^-x); for large |x| approach sign(x)
double coth;
if (Math.abs(heOverA) > 30) {
coth = (heOverA > 0) ? 1.0 : -1.0;
} else {
double ep = Math.exp(heOverA);
double em = Math.exp(-heOverA);
coth = (ep + em) / (ep - em);
}
man = coth - JA_SHAPE / he;
}
// clamp anhysteretic to [-1,1] (numerical safety near saturation)
if (man > 1) man = 1;
if (man < -1) man = -1;

double diff = man - mIrr;
double denom = sign - JA_COUPLING * diff;
// floor denominator to avoid singularity at loop tips
if (Math.abs(denom) < 1e-3) denom = sign * 1e-3;
double dMirr = diff / denom * dh;
mIrr += dMirr;
if (mIrr > 1) mIrr = 1;
if (mIrr < -1) mIrr = -1;

double newM = mIrr + reversibility * (man - mIrr);

// effective inductance from numerical slope dm/dh times L0
double dm_dh;
if (Math.abs(dh) > 1e-6) {
dm_dh = (newM - mTotal) / dh;
} else {
// no change this step: fall back to anhysteretic tangent
dm_dh = 1.0; // = 1/(3*JA_SHAPE) with an=1/3
}
// floor and cap incremental slope for convergence
if (dm_dh < 0.01) dm_dh = 0.01;
if (dm_dh > 3.0) dm_dh = 3.0;
lastLeff = inductance * dm_dh;

prevM = mTotal;
mTotal = newM;
prevI = i;
}

void stamp(int n0, int n1) {
// inductor companion model using trapezoidal or backward euler
// approximations (Norton equivalent) consists of a current
Expand All @@ -67,7 +155,7 @@ void stamp(int n0, int n1) {
// The oscillation is a real problem in circuits with switches.
nodes[0] = n0;
nodes[1] = n1;
if (saturationCurrent > 0) {
if (saturationCurrent > 0 || hasHysteresis()) {
// nonlinear: conductance changes with current, stamped in doStep()
sim.stampNonLinear(nodes[0]);
sim.stampNonLinear(nodes[1]);
Expand All @@ -82,10 +170,18 @@ void stamp(int n0, int n1) {
sim.stampRightSide(nodes[0]);
sim.stampRightSide(nodes[1]);
}
boolean nonLinear() { return saturationCurrent > 0; }
boolean nonLinear() { return saturationCurrent > 0 || hasHysteresis(); }

void startIteration(double voltdiff) {
if (saturationCurrent > 0) {
if (hasHysteresis()) {
// advance JA state once per timestep, then use resulting lastLeff
advanceHysteresis(current);
double lEff = lastLeff;
if (isTrapezoidal())
compResistance = 2*lEff/sim.timeStep;
else
compResistance = lEff/sim.timeStep;
} else if (saturationCurrent > 0) {
// recompute companion resistance from current-dependent inductance
double lEff = calcEffectiveInductance(current);
if (isTrapezoidal())
Expand All @@ -108,10 +204,14 @@ void startIteration(double voltdiff) {
return current;
}
void doStep(double voltdiff) {
if (saturationCurrent > 0) {
if (saturationCurrent > 0 || hasHysteresis()) {
// stamp companion conductance (matrix was restored to origMatrix)
sim.stampConductance(nodes[0], nodes[1], 1.0/compResistance);
}
sim.stampCurrentSource(nodes[0], nodes[1], curSourceValue);
}

// magnetization state (normalized -1..1) for info panel / debugging
double getMagnetization() { return mTotal; }
double getIrreversibleMagnetization() { return mIrr; }
}
66 changes: 51 additions & 15 deletions src/com/lushprojects/circuitjs1/client/InductorElm.java
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
/*
/*
Copyright (C) Paul Falstad and Iain Sharp

This file is part of CircuitJS1.

CircuitJS1 is free software: you can redistribute it and/or modify
Expand All @@ -27,35 +27,46 @@ class InductorElm extends CircuitElm {
double inductance;
double initialCurrent;
double saturationCurrent; // 0 = disabled (linear)
double coerciveCurrent; // 0 = hysteresis disabled
double reversibility; // JA c, 0..1
public InductorElm(int xx, int yy) {
super(xx, yy);
ind = new Inductor(sim);
inductance = 1;
ind.setup(inductance, current, flags, saturationCurrent);
reversibility = 0.1;
ind.setup(inductance, current, flags, saturationCurrent, coerciveCurrent, reversibility);
}
public InductorElm(int xa, int ya, int xb, int yb, int f,
StringTokenizer st) {
super(xa, ya, xb, yb, f);
ind = new Inductor(sim);
reversibility = 0.1;
inductance = new Double(st.nextToken()).doubleValue();
current = new Double(st.nextToken()).doubleValue();
try {
initialCurrent = new Double(st.nextToken()).doubleValue();
saturationCurrent = new Double(st.nextToken()).doubleValue();
coerciveCurrent = new Double(st.nextToken()).doubleValue();
reversibility = new Double(st.nextToken()).doubleValue();
} catch (Exception e) {}
ind.setup(inductance, current, flags, saturationCurrent);
ind.setup(inductance, current, flags, saturationCurrent, coerciveCurrent, reversibility);
}
int getDumpType() { return 'l'; }
String dump() {
return super.dump() + " " + inductance + " " + current + " " + initialCurrent + " " + saturationCurrent;
return super.dump() + " " + inductance + " " + current + " " + initialCurrent + " " + saturationCurrent
+ " " + coerciveCurrent + " " + reversibility;
}

void dumpXml(Document doc, Element elem) {
super.dumpXml(doc, elem);
XMLSerializer.dumpAttr(elem, "l", inductance);
XMLSerializer.dumpAttr(elem, "ic", initialCurrent);
if (saturationCurrent != 0)
XMLSerializer.dumpAttr(elem, "isat", saturationCurrent);
if (coerciveCurrent != 0) {
XMLSerializer.dumpAttr(elem, "ich", coerciveCurrent);
XMLSerializer.dumpAttr(elem, "hrev", reversibility);
}
}

void dumpXmlState(Document doc, Element elem) {
Expand All @@ -68,7 +79,9 @@ void undumpXml(XMLDeserializer xml) {
initialCurrent = xml.parseDoubleAttr("ic", initialCurrent);
current = xml.parseDoubleAttr("i", current);
saturationCurrent = xml.parseDoubleAttr("isat", 0);
ind.setup(inductance, current, flags, saturationCurrent);
coerciveCurrent = xml.parseDoubleAttr("ich", 0);
reversibility = xml.parseDoubleAttr("hrev", 0.1);
ind.setup(inductance, current, flags, saturationCurrent, coerciveCurrent, reversibility);
}

void setPoints() {
Expand Down Expand Up @@ -110,14 +123,26 @@ void doStep() {
ind.doStep(voltdiff);
}
void getInfo(String arr[]) {
arr[0] = (saturationCurrent > 0) ? "inductor (sat)" : "inductor";
if (ind.hasHysteresis())
arr[0] = "inductor (hyst)";
else if (saturationCurrent > 0)
arr[0] = "inductor (sat)";
else
arr[0] = "inductor";
getBasicInfo(arr);
arr[3] = "L = " + getUnitText(inductance, "H");
arr[4] = "P = " + getUnitText(getPower(), "W");
if (saturationCurrent > 0) {
int row = 5;
if (saturationCurrent > 0 && !ind.hasHysteresis()) {
double lEff = ind.calcEffectiveInductance(current);
arr[5] = "Leff = " + getUnitText(lEff, "H");
arr[6] = "Isat = " + getUnitText(saturationCurrent, "A");
arr[row++] = "Leff = " + getUnitText(lEff, "H");
arr[row++] = "Isat = " + getUnitText(saturationCurrent, "A");
}
if (ind.hasHysteresis() && row < arr.length) {
arr[row++] = "Leff = " + getUnitText(ind.calcEffectiveInductance(current), "H");
if (row < arr.length)
arr[row++] = "Ic = " + getUnitText(coerciveCurrent, "A")
+ ", M = " + showFormat.format(ind.getMagnetization());
}
}
public EditInfo getEditInfo(int n) {
Expand All @@ -133,6 +158,10 @@ public EditInfo getEditInfo(int n) {
return new EditInfo("Initial Current (on Reset) (A)", initialCurrent);
if (n == 3)
return new EditInfo("Saturation Current (A) (0=none)", saturationCurrent);
if (n == 4)
return new EditInfo("Coercive Current (A) (0=no hysteresis)", coerciveCurrent);
if (n == 5)
return new EditInfo("Hysteresis Reversibility (0-1)", reversibility);
return null;
}

Expand All @@ -149,18 +178,25 @@ public void setEditValue(int n, EditInfo ei) {
initialCurrent = ei.value;
if (n == 3 && ei.value >= 0)
saturationCurrent = ei.value;
ind.setup(inductance, current, flags, saturationCurrent);
if (n == 4 && ei.value >= 0)
coerciveCurrent = ei.value;
if (n == 5) {
if (ei.value < 0) ei.value = 0;
if (ei.value > 1) ei.value = 1;
reversibility = ei.value;
}
ind.setup(inductance, current, flags, saturationCurrent, coerciveCurrent, reversibility);
}

int getShortcut() { return 'L'; }
public double getInductance() { return inductance; }
void setInductance(double l) {
inductance = l;
ind.setup(inductance, current, flags, saturationCurrent);
ind.setup(inductance, current, flags, saturationCurrent, coerciveCurrent, reversibility);
}
void setSaturationCurrent(double isat) {
saturationCurrent = isat;
ind.setup(inductance, current, flags, saturationCurrent);
ind.setup(inductance, current, flags, saturationCurrent, coerciveCurrent, reversibility);
}
double getSaturationCurrent() { return saturationCurrent; }
}
10 changes: 10 additions & 0 deletions war/circuits/hysteresis-loop.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,10 @@
$ 1 0.000005 10.20027730826997 50 5 50 5e-11
v 192 288 192 160 0 1 50 2 0 0 0.5
w 192 160 288 160 0
l 288 160 400 160 0 0.5 0 0 0 0.25 0.1
w 400 160 400 288 0
r 400 288 320 288 0 1
w 320 288 192 288 0
368 400 160 400 288 0 Vind
o 0 32 0 4355 0.5 0.2 0 2 0 3
o 5 32 0 4099 2 0.2 1 2 5 3
Loading