Chapter 14 — Core saturation and inrush

Chapter 13 left one branch of the transformer suspiciously simple. The magnetizing inductance $L_m$ was treated as a constant, which it is not: iron saturates, and when it does the inductance can fall by three orders of magnitude in the space of a tenth of a per-unit of flux. This chapter is about that collapse — how to describe it, how to put it into a solver without wrecking the factorize-once economics of Chapter 3, and what it does when you close a breaker onto an unloaded transformer.

It is also the first genuinely nonlinear element in this course. The technique used here — freeze the stamped conductance and carry the nonlinearity in a current injection — is worth learning on its own account. It recurs in the machine chapters and in every EMT model whose nonlinearity is smooth.

Learning objectives

By the end of this chapter you should be able to:

• Sketch the $\lambda$–$i$ characteristic of an iron core and identify the knee and the air-core slope.
• Explain how saturation is implemented as a Norton current injection that leaves the conductance matrix untouched, and why that matters.
• Choose between the exponential, SE-point, and two-slope saturation definitions, and say which studies each is honest for.
• Predict the flux offset that follows from closing a breaker at a given point on the voltage wave.
• Read an inrush current peak off the saturation curve arithmetically, and explain why the decay is not a single exponential.

14.1 The core is nonlinear

A winding's flux linkage $\lambda$ and its current $i$ are related by the iron's magnetization curve. At low flux the iron's domains align easily and a tiny current produces a great deal of flux, so the curve is steep and $L_m=\lambda /i$ is enormous. As the domains run out of room to align, the curve bends over — the knee — and past it the iron contributes essentially nothing more. What is left is the flux the coil would produce with no iron at all, so the slope settles at the air-core inductance $L_{\text{air}}$, set purely by winding geometry.

The two slopes are separated by a startling factor. For the transformer used throughout this module, the unsaturated magnetizing reactance is 100 pu, and its air-core inductance is 0.1 pu — a thousand to one. Below the knee the branch is nearly an open circuit; above it, it is close to a short.

The magnetizing characteristic used by the energization lab: unsaturated $L_m$ below the knee at 1.2 pu, air-core slope above it. Below the knee the branch draws 3.55 A at 1.0 pu flux and 4.3 A at the knee, which on this axis is indistinguishable from zero. Above the knee every extra 0.1 pu of flux costs another 355 A. The marked point is the lab's own crest, 1.69 pu and 1745 A.

Two things a real core does that this curve does not: hysteresis (the rising and falling branches differ, which is where hysteresis loss comes from) and remanence (a de-energized core retains flux, so the next energization starts somewhere other than zero). NumaSim models neither. The core-loss conductance $G_c$ from Chapter 13 accounts for the energy lost to hysteresis and eddy currents, but the characteristic itself is single-valued and the flux state starts at zero every run. Where that limitation bites is flagged in §14.4.

14.2 Saturation without refactorizing

Here is the tempting implementation, and the reason nobody uses it. Saturation makes $L_m$ a function of flux; the magnetizing companion conductance is $G_m=\mathrm{\Delta }t/2L_m$; so recompute $G_m$ each step from the present flux and restamp. It is correct, it is obvious, and it destroys the central economy of the whole method. Chapter 3's argument was that $\mathbf{G}$ is constant between topology changes, so the expensive LU factorization happens once and every step is a cheap substitution. A conductance that changes every step means a factorization every step, on every transformer in the network.

The standard fix is to leave the conductance alone and put the entire nonlinearity into the current injection, which is on the right-hand side and therefore free. Write the saturation as a multiplicative derating of the unsaturated inductance,

$$L_{\text{sat}}(\psi )=\frac{L_m}{1+S(\psi )},$$

where $S(\psi )\ge 0$ is the saturation function and $\psi$ is the magnetizing flux in per unit. The true branch current is then $\psi /L_{\text{sat}}=(1+S)\psi /L_m$, which splits into

$$i_{\text{mag}}=\underset{\text{linear companion}}{\underbrace{\frac{\psi }{L_m}}}+\underset{\text{extra injection }{I}_{\text{sat}}}{\underbrace{S(\psi )\frac{\psi }{L_m}}}.$$

The first term is exactly the companion of Chapter 2 with the unsaturated $L_m$ — a constant conductance, stamped once. The second is a number computed from the previous step's flux and added to the injection vector. The matrix never changes. Saturation costs one exponential (or one comparison and a multiply) per leg per step and nothing else.

Verify that the split is exact rather than approximate, using the two-slope law of §14.3, whose $S$ is $S(\psi )=\frac{\psi -{\psi }_{\text{knee}}}{\psi }(\frac{L_m}{L_{\text{air}}}-1)$ above the knee:

$$I_{\text{sat}}=S\frac{\psi }{L_m}=(\psi -{\psi }_{\text{knee}})(\frac{1}{L_{\text{air}}}-\frac{1}{L_m}),$$$$\frac{\psi }{L_m}+I_{\text{sat}}=\frac{{\psi }_{\text{knee}}}{L_m}+\frac{\psi -{\psi }_{\text{knee}}}{L_{\text{air}}},$$

which is the piecewise-linear characteristic, exactly. The injection is not a linearization or a correction term; it is an algebraic rearrangement that moves the whole nonlinearity from the left-hand side to the right.

Contrast this with the switches of Chapter 5. A switch changes the network's topology, so its conductance change cannot be absorbed into an injection and a refactorization is unavoidable. Saturation is a smooth change in an element's value along a fixed topology, and smooth changes can always be shifted to the right-hand side. That is the dividing line between the two treatments.

A trap worth knowing about

The injection must be added at the point of use and must never be folded back into the companion's history state. The companion recursion feeds its own output forward as history, $I^{\text{hist}}=\alpha i+Gv$, so a term added to the stored current re-enters on every subsequent step and integrates instead of acting once. Because that spurious accumulation is per step rather than per unit time, the resulting error scales like $1/\mathrm{\Delta }t$: the simulated inrush amplitude changes when you refine the step and never converges. NumaSim keeps the linear companion current and the saturation injection in separate variables for exactly this reason, and the lab below checks convergence under step refinement as a matter of routine.

14.3 Three ways to describe the curve, and when each is honest

The component's sat_definition parameter offers three ways to specify $S(\psi )$. They are not interchangeable, and picking the wrong one is the most common way to get a plausible-looking but badly wrong answer.

Exponential. The classical form, taking coefficients directly:

$$S(\psi )=A_{\text{sat}}{e}^{B_{\text{sat}}(\psi -{\psi }_{T1})}(\psi >{\psi }_{T1}),0\text{ otherwise.}$$

${\psi }_{T1}$ is a fitting origin, not a physical knee — the standard 0.8 pu sits below rated flux and is harmless there because $S(0.8)=0$ and grows smoothly from it.

SE points. The same exponential, specified the way a machine or transformer datasheet does it. $SE(V)$ is the fractional extra magnetizing current needed to reach flux $V$ on the real open-circuit curve, relative to the straight air-gap line. Two points fix the two coefficients:

$$B_{\text{sat}}=\frac{\ln (SE(1.2)/SE(1.0))}{0.2},A_{\text{sat}}=\frac{SE(1.0)}{e^{B_{\text{sat}}(1.0-{\psi }_{T1})}}.$$

Typical values are $SE(1.0)$ around 0.05–0.15 and $SE(1.2)$ around 0.2–0.5. This is the same law as the exponential option, entered in datasheet coordinates.

Two-slope (air-core). The piecewise-linear characteristic, specified by a hard knee ${\psi }_{\text{knee}}$ and the air-core inductance $L_{\text{air}}$:

$$i(\psi )=\{\begin{array}{ll}\psi /L_m,& |\psi |\le {\psi }_{\text{knee}},\\ \mathrm{sign}(\psi )[{\displaystyle \frac{{\psi }_{\text{knee}}}{L_m}}+{\displaystyle \frac{|\psi |-{\psi }_{\text{knee}}}{L_{\text{air}}}}],& |\psi |>{\psi }_{\text{knee}}.\end{array}$$

Note that ${\psi }_{\text{knee}}$ here is a physical break, unlike the exponential's ${\psi }_{T1}$: it must sit above the transformer's normal operating flux, or the unit saturates every cycle in ordinary service. Real cores knee over between 1.1 and 1.25 pu.

So which one? The answer follows from where each is fitted. The exponential is fitted between 1.0 and 1.2 pu, which is where over-excited steady-state operation lives, and it has no air-core asymptote — extrapolated far past its fitting range it keeps reporting an inductance tens of times the leakage, because nothing in the functional form knows the iron eventually runs out.

Energization drives the flux to nearly 2 pu, far outside that range, and the answer there is governed entirely by how close the saturated inductance gets to air core. The lab below gives the numbers: on the two-slope curve the same energization peaks at 1745 A, and on an exponential curve it peaks at 347 A — a factor of five, from nothing but the choice of functional form. Neither number is a bug. Use the exponential (or SE points) for steady-state over-excitation near 1.0–1.2 pu; use two-slope for energization, inrush, and anything that pushes the core deep.

14.4 Flux is the integral of voltage

Everything about inrush follows from one relation:

$$\psi (t)=\psi (0)+{\int }_0^{t}v(\tau )d\tau .$$

In steady state, with $v=V_{pk}\cos \omega t$, the flux is $\psi =(V_{pk}/\omega )\sin \omega t$ — a sinusoid centred on zero, peaking at 1.0 pu by the definition of the flux base. Crucially, the flux is 90° behind the voltage, so the flux is at zero when the voltage is at its crest, and at its crest when the voltage is at zero.

Now close a breaker. The flux cannot jump, so it starts from whatever it was, which for a de-energized core the model treats as zero. Two cases:

• Close at the voltage crest. The steady-state flux at that instant is zero, so the initial condition already matches. The flux slides straight into its steady-state sinusoid, never exceeds 1.0 pu, and there is no inrush at all.
• Close at a voltage zero. The steady-state flux at that instant is at its own crest, $\pm 1.0$ pu, but the actual flux is zero. The difference — a full 1.0 pu — appears as a DC offset that rides on the sinusoid, so the flux swings from 0 toward 2.0 pu. This is the worst case, and it is not marginal: the core is asked to carry twice its design flux.

Top: the three flux linkages over the first three cycles of the energization lab, with the $\pm 1.2$ pu knee marked. Phase A carries an offset near $+0.83$ pu, so its crest reaches 1.69 pu and spends a few milliseconds per cycle above the knee — precisely the window in which current flows. Bottom: the per-cycle crest, trough and midpoint of ${\psi }_a$ over the full 8 s, showing the offset stalling once the crest lands on the knee.

Two per-unit numbers are worth internalizing. A 2.0 pu flux crest on a core that knees at 1.2 pu is 0.8 pu past the knee, and past the knee the inductance is $L_{\text{air}}=0.1$ pu, so the current is $0.8/0.1=8$ pu — eight times rated peak, from a healthy transformer with nothing wrong with it. And on a three-phase system you cannot avoid it by timing: the three voltages are 120° apart, so a simultaneous close puts at most one phase near its crest and the other two somewhere bad.

What NumaSim cannot reproduce

The flux state starts at zero. There is no remanent-flux parameter, so the genuine worst case — a core left with, say, 0.7 pu of residual flux in the same direction as the offset, which drives the crest toward 2.7 pu — cannot be set up. Real inrush studies for controlled switching do model remanence, and their peaks are correspondingly higher than anything this chapter produces. Read the numbers below as the no-remanence worst case.

14.5 Lab: energizing an unloaded transformer

Open the transformer energization sample in simulator →

A 230 kV source snaps onto the unloaded 100 MVA transformer at $t=0$, at the instant phase A's voltage passes through zero — the worst case of §14.4. The LV winding sits on 1 MΩ per phase, so every ampere measured is going into the core. The saturation curve is two-slope with ${\psi }_{\text{knee}}=1.2$ pu and $L_{\text{air}}=0.1$ pu. Rated HV current for this unit is 251 A RMS, i.e. 355 A peak, and every multiple below is against that peak.

The first three cycles of magnetizing current. Phase A crests at 1745 A — 4.9 times rated peak — 8.4 ms after closing, then repeats every cycle. The shape is the tell: current is essentially zero for most of the cycle and appears as a one-sided pulse a few milliseconds wide, because the branch only conducts while the flux is over the knee.

Read the peak off the curve. Run it and note the flux crest: ${\psi }_a$ reaches 1.690 pu. That is 0.490 pu past the 1.2 pu knee. One per-unit of flux across a 0.1 pu inductance is 10 pu of current, and 1 pu of current is 355 A peak, so the slope past the knee is 3550 A per pu of flux:

$$0.490\times 3550=1740\text{A},$$

on top of the 4 A the core was already drawing at the knee. The simulator reports 1745 A. You have just predicted a nonlinear transient's peak with two multiplications, because past the knee the circuit is linear again — just with a different inductance.

It is a converged number, not a step-size artefact. Refine from 50 µs to 12.5 µs and the peak moves from 1745.0 A to 1744.8 A. That is the check §14.2's warning demands: an injection that had leaked into the companion history would have moved this number by a factor of four.

Confirm what saturation is contributing. Turn the transformer's saturation parameter off and rerun. The flux still doubles — it reaches 1.998 pu, because flux is the integral of voltage and knows nothing about iron — but a linear core answers a doubled flux with a doubled current, so the peak is 7.1 A instead of 1745 A.

The same energization with saturation on and off, on a log axis. A factor of 250 in peak current, from one nonlinearity. The two traces converge at the right-hand edge because by then the saturated run has fallen back onto the linear magnetizing branch.

Now work the four variations:

1. Close at the voltage crest. Set the source phase to 90°. Phase A's inrush vanishes completely — 3.6 A, indistinguishable from ordinary magnetizing current, because its flux starts at its own crest and never needs an offset. Phases B and C take the hit instead, at 1446 A and 1466 A. You cannot null all three with one simultaneous close, which is the entire reason controlled-switching schemes stagger the poles.
2. Ramp the source. A ramp_time of 40 ms brings the voltage up over a couple of cycles and the flux tracks it without ever building an offset: 3.7 A on every phase. This is the standard way to suppress inrush in a study where it is not the subject.
3. Switch the curve. Set sat_definition to Exponential. The same energization now peaks at 347 A with the flux climbing to 1.94 pu — §14.3's factor of five, made concrete.
4. Sweep the air-core inductance. l_air_pu is the single number that sets the peak, and the relationship is close to inverse: 0.2 pu gives 1083 A, 0.1 pu gives 1745 A, 0.08 pu gives 1989 A. For a power transformer $L_{\text{air}}$ typically lands at one to three times the short-circuit reactance, so for this unit's 0.1 pu leakage the plausible range is about 0.1–0.3 pu.

14.6 Why the decay is not a single exponential

Per-cycle peak of each magnetizing current over the whole 8 s window. Phase A halves in 0.38 s, is under a tenth of its first peak by 2.0 s, and finishes at 8 A. The decay is visibly faster than exponential early and slower than exponential late, because the mechanism that damps it switches itself off.

What removes the flux offset is $I\cdot R$ — resistance in the source and the winding dissipating the DC component. That only bites while there is a large $I$, which is to say while the crest is over the knee and the branch is drawing hundreds of amperes. As the offset shrinks, the crest sinks toward the knee, the pulses narrow, and the very mechanism damping the offset weakens along with them.

The decay is therefore self-arresting. By 3 s the crest is down to 1.22 pu and the offset has all but stopped at 0.23 pu; five seconds later they are 1.20 and 0.20. The residue is not going away on this time scale, because once the crest is on the knee the only thing left to drain it is the unsaturated $L_m/R$ of the loop — roughly 140 H over 1.8 Ω, a time constant near 80 s. On an 8 s window that looks frozen. It is physics, not a numerical artefact, and it is the same effect the open-circuit test in Chapter 13 ran into.

Source impedance is the lever here, and it works the opposite way to intuition. Stiffen the source to 0.05 Ω + 1.4 mH and the first peak goes up slightly, to 1856 A, because there is less impedance between the source and the core — but the run is still at 22 A after 8 s instead of 8 A, because the resistance that was draining the offset is gone. A stiffer system means a bigger inrush that lasts longer.

14.7 Why anyone cares

Nothing in the lab is a fault. It is a healthy transformer drawing nearly five times rated current on its first cycle and staying above rated for a full second. Three practical consequences follow.

Differential protection. A transformer differential relay compares primary and secondary currents and trips when they disagree. During inrush they disagree completely — the primary carries 1745 A and the secondary carries nothing, because the difference current is magnetizing current. Every differential relay therefore needs inrush restraint, classically by detecting the second-harmonic content that the one-sided pulse shape generates in abundance. The pulse shape in §14.5's first figure is what that logic is looking at.

Controlled switching. If you can choose the closing instant per pole, you can land each phase near its own voltage crest and suppress the offset. That is variation 1 of the lab, and the reason it must be done pole by pole is variation 1's own result: nulling phase A moved the problem to B and C. Real controlled-switching schemes also compensate for remanent flux, which §14.4 notes is outside NumaSim's model today.

Voltage quality and mechanical stress. The inrush current flows through the source impedance, so it depresses and distorts the local voltage for as long as it lasts — a visible dip on a weak system. And the same current produces mechanical forces in the windings proportional to its square, so a transformer that is switched frequently accumulates mechanical fatigue that one switched rarely does not.

Beyond this chapter lie two phenomena built from the same nonlinearity: ferroresonance, where a saturable core exchanges energy with system capacitance and settles into a sustained, often chaotic overvoltage, and core-saturation instability, where a small DC bias — from a geomagnetic disturbance or a nearby HVDC ground return — pushes the core into saturation on alternate half cycles and floods the system with even harmonics. Both require exactly the model this chapter built, and both are studied with EMT for the same reason: there is no phasor description of a waveform that spends most of its cycle at zero.

14.8 Summary

• The $\lambda$–$i$ curve of an iron core is steep below the knee and falls to the air-core slope above it; for a power transformer the two inductances differ by about a factor of a thousand.
• Saturation is implemented as a Norton current injection, $I_{\text{sat}}=S(\psi )\psi /L_m$, added to the right-hand side while the stamped conductance stays at its unsaturated value. The rearrangement is algebraically exact and preserves the factorize-once structure of the solver.
• The injection must never be folded into the companion's history state, or the error accumulates per step and scales like $1/\mathrm{\Delta }t$.
• The exponential and SE-point definitions describe the same curve and are fitted for 1.0–1.2 pu; the two-slope definition carries an air-core asymptote and is the one to use for energization. On the same energization they differ by a factor of five.
• Flux is the integral of voltage, so closing at a voltage zero forces a 1.0 pu DC offset and drives the flux toward 2.0 pu; closing at the crest produces no inrush at all.
• Inrush decay is self-arresting: the offset only drains while the crest is above the knee, so it stalls at a residue that then decays on the unsaturated $L_m/R$ time constant of tens of seconds.

14.9 Problems

Problem 14.1. A transformer's core knees at 1.15 pu and its air-core inductance is 0.15 pu. An energization drives the flux crest to 1.85 pu. Estimate the peak magnetizing current in per unit and as a multiple of rated peak.

Solution 14.1

The flux is $1.85-1.15=0.70$ pu past the knee. Above the knee the inductance is 0.15 pu, so the current is $0.70/0.15=4.67$ pu, plus the small pre-knee current ($1.15/100\approx 0.012$ pu), which is negligible.

The answer is 4.67 pu, i.e. 4.67 times rated peak. Note that the peak scales inversely with $L_{\text{air}}$: halving it to 0.075 pu would give 9.3 pu.

Problem 14.2. Show that adding the saturation injection $I_{\text{sat}}=S(\psi )\psi /L_m$ to the linear companion current reproduces the exponential law's intended branch current $\psi /L_{\text{sat}}$ with $L_{\text{sat}}=L_m/(1+S)$.

Solution 14.2

The intended current is

$$i=\frac{\psi }{L_{\text{sat}}}=\frac{\psi (1+S)}{L_m}=\frac{\psi }{L_m}+S\frac{\psi }{L_m}.$$

The first term is the unsaturated companion's current and the second is $I_{\text{sat}}$ by definition, so the sum is exact for any $S(\psi )$ — the derivation never used the functional form. That is why the same code path serves the exponential, SE-point, and two-slope laws: they differ only in how $S$ is evaluated.

Problem 14.3. A datasheet gives $SE(1.0)=0.08$ and $SE(1.2)=0.26$. Compute $A_{\text{sat}}$ and $B_{\text{sat}}$ for ${\psi }_{T1}=0.8$, then evaluate $S$ at 1.2 pu and at 1.9 pu. Comment on the second result.

Solution 14.3

$B_{\text{sat}}=\ln (0.26/0.08)/0.2=\ln (3.25)/0.2=1.1787/0.2=5.894$.

$A_{\text{sat}}=0.08/e^{5.894(1.0-0.8)}=0.08/e^{1.1787}=0.08/3.250=0.0246$.

At $\psi =1.2$: $S=0.0246e^{5.894(0.4)}=0.0246e^{2.358}=0.0246\times 10.57=0.260$ — recovering $SE(1.2)$, as the fit requires.

At $\psi =1.9$: $S=0.0246e^{5.894(1.1)}=0.0246e^{6.483}=0.0246\times 654=16.1$, so $L_{\text{sat}}=L_m/17.1$. With $L_m=100$ pu that is 5.8 pu — still fifty-eight times the 0.1 pu air-core inductance the iron should have collapsed to. The exponential has no asymptote, so it understates deep saturation badly. This is the arithmetic behind §14.3's warning.

Problem 14.4. A 60 Hz transformer is energized at the instant its voltage passes through zero, from a de-energized state. Sketch the flux over the first two cycles and state the peak in per unit. How would the answer change if the core had 0.6 pu of remanent flux in the same direction?

Solution 14.4

Steady-state flux lags voltage by 90°, so at a voltage zero the steady-state flux is at $\pm 1.0$ pu. The actual flux is 0, so a $\mp 1.0$ pu DC offset is required to reconcile them. The flux is therefore $\psi (t)=1.0-\cos \omega t$ (taking the positive case), which starts at 0, peaks at 2.0 pu half a cycle later at $t=8.33$ ms, returns to 0, and repeats.

With 0.6 pu of remanence aiding, the flux starts at 0.6 and the same 1.0 pu swing carries it to 2.6 pu. Since the peak current goes as (crest − knee)/$L_{\text{air}}$, moving the crest from 2.0 to 2.6 pu with a 1.2 pu knee and 0.1 pu air core raises the current from 8 pu to 14 pu — a 75% increase from remanence alone. This is why controlled-switching schemes measure residual flux rather than just timing the close.

Problem 14.5. Why does making the source stiffer increase the inrush peak and lengthen the decay? Answer in terms of the two roles source impedance plays.

Solution 14.5

Source impedance does two separate jobs. Its reactance limits the peak: above the knee the loop is source inductance plus leakage plus $L_{\text{air}}$, so removing source inductance leaves less total inductance and admits more current. Its resistance drains the DC flux offset through $I\cdot R$ losses, so removing source resistance removes the main damping path.

A stiffer source therefore gives a larger first peak (less limiting reactance) and a slower decay (less damping resistance) — the two effects push in the same unhelpful direction. The lab measures both: 1856 A instead of 1745 A at the first peak, and 22 A instead of 8 A at $t=8$ s.

Problem 14.6. An engineer models saturation by recomputing $L_m$ each step and restamping the conductance matrix. The results are correct but the simulation is fifty times slower than expected. Explain, and give the fix.

Solution 14.6

Restamping changes $\mathbf{G}$, which invalidates the LU factorization computed at the start of the run. The solver must refactorize every step instead of doing a cheap forward/backward substitution, and factorization is the expensive operation — for a network of $n$ nodes it costs $O(n^{1.5})$ or worse against $O(n)$ for a substitution. On a network with several transformers, essentially every step becomes a factorization.

The fix is §14.2: keep the stamped conductance at the unsaturated$L_m$ and add $I_{\text{sat}}=S(\psi )\psi /L_m$ to the injection vector. The result is algebraically identical and the matrix never changes, so the factorization survives the whole run.

14.10 References

• P. Kundur, Power System Stability and Control, McGraw-Hill — Section 3.8 for the open-circuit saturation curve and the $SE$ point convention; Section 6.2 for transformer core representation.
• J. Arrillaga and N. R. Watson, Power Systems Electromagnetic Transients Simulation, IET Power and Energy Series 39 — Chapter 7 for saturation as a compensating current source and the reason it is preferred to varying the inductance.
• H. W. Dommel, Electromagnetic Transients Program (EMTP) Theory Book, Bonneville Power Administration — nonlinear inductances, piecewise-linear characteristics, and hysteresis models.

---

Previous: Chapter 13 — Transformer models and parameters · Next: Chapter 15 — The synchronous machine.